/*******************************************************************************
*
* McXtrace, X-ray tracing package
*         Copyright, All rights reserved
*         DTU Physics, Kgs. Lyngby, Denmark
*         Synchrotron SOLEIL, Saint-Aubin, France
*
* Component: Source_spectra
*
* %Identification
* Written by: Erik Knudsen 
* Date: November 11, 2019
* Origin: Risoe
* Release: McXtrace 1.5
*
* Specialized X-ray source for reading in SPECTRA 10 source definitions
*
* %Description
*
* This is a source component for connecting SPECTRA 10-output files with McXtrace.
* json-style SPECTRA 11 output files are not yet supported.
* 
* SPECTRA is an application software to calculate optical properties of synchrotron 
* radiation (SR) emitted from bending magnets, wigglers (conventional and elliptical) 
* and undulators (conventional, helical, elliptical and figure-8). Calculations 
* of radiation from an arbitrary magnetic field distribution are also available. 
* Parameters on the electron beam and the source can be edited completely on 
* graphical user interfaces (GUIs) and it is possible to show the calculation 
* result graphically. The energy spectrum and radiation power after transmitting 
* various filters and convolution of detector's resolution are also available. 
* See <a href="http://spectrax.org/spectra/">SPECTRA</a>.
*
* If the source is symmetric in x and/or y it is possible to speed up the spectra
* calculations by only including one half-plane or quadrant. The other side/quadrants will then
* be mirrored by McXtrace.
*
* %BUGS
* Absolute intensity of 4D (x,y,x',y') is nor correctly normalized.
*
* %Parameters
* E0:             [keV] Mean energy of X-rays.
* dE:             [keV] Energy spread of X-rays.
* Emin:           [keV] Energy of low end of the Spectra-calculated data.
* Emax:           [keV] Energy of high end of the Spectra-calculated data.
* nE:             [int] Number of steps in the spectra-calculations.
* randomphase:    [0/1] If !=0 the photon phase is chosen randomly.
* nx:             [int] Number of grid points along x in datafiles. If zero this is computed from the files.
* ny:             [int] Number of grid points along y in datafiles. If zero this is computed from the files.
* npx:            [int] Number of grid points along x' in datafiles. If zero this is computed from the files.
* npy:            [int] Number of grid points along y' in datafiles. If zero this is computed from the files.
* phase:          [rad] Value of the photon phase (only used if randomphase==0).
* verbose:        [0/1] If non-zero output more warning messages.
* initial_serial: [int] First serial number of the series of spectra files.
* symmetricx:     [0/1] If nonzero the source is mirrored in the x-axis. This to allow smaller spectra-calculations.
* symmetricy:     [0/1] If nonzero the source is mirrored in the y-axis. This to allow smaller spectra-calculations.
* spectra_stem_x: [str] Filename stem of x-projection of source distribution. -n.xxx will be added where n is a serial number and xxx spectra_suffix.
* spectra_stem_y: [str] Filename stem of y-projection of source distribution. -n.xxx will be added where n is a serial number and xxx spectra_suffix.
* spectra_stem:   [str] Filename stem of x,x',y,y'-source distribution distribution. -n.xxx will be added where n is a serial number and xxx spectra_suffix.
* spectra_suffix: [str] Suffix of spectra output files.
* flag4d:         [0/1] Use either (0) x,y-projections or (1) full 4D x,y,x',y' datafiles.
* noinit:         [0/1] Do no initialize the component. Can be usefiul in conjunction with a deactivating WHEN-clause.
*
* CALCULATED PARAMETERS:
*
* %Link
* Tanaka, J. Synchrotron Rad. (2001). 8, 1221-1228. https://doi.org/10.1107/S090904950101425X
* http://spectrax.org/spectra/
*
* %End
*******************************************************************************/

DEFINE COMPONENT Source_spectra

SETTING PARAMETERS (
    string spectra_stem_x="", string spectra_stem_y="", string spectra_stem="", string spectra_suffix="dsc",
    E0=0, dE=0, Emin,Emax, int nE, int randomphase=1, phase=0,
    int nx=0, int ny=0, int npx=0, int npy=0, int initial_serial=1, int symmetricx=0, int symmetricy=0, int verbose=0,
    int flag4d=0, int noinit=0)

/* X-ray parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */ 

SHARE
%{
  %include "read_table-lib";

  int source_spectra_find_offset(char * fn){
      /*find the first line that starts with [-0-9], i.e. can be considered a number*/
      char line[3][512];
      int linecheck[3],done=0;
      long pos[3]={0,-1,-2};
      double buf[6];
      FILE *fs;

      if( (fs=Open_File(fn,"rb",NULL))==NULL){
        /*Open_File from read_table-lib searches the McXtrace library. Will report error on failure, so just exit.*/
        exit(-1);
      }

      /*read lines and save position 3 lines back, when three consecutive line have 3
        columns we have an offset*/
      line[0][0]='\0';line[1][0]='\0';line[2][0]='\0';
      line[0][511]='\0';line[1][511]='\0';line[2][511]='\0';

      do {
          pos[2]=pos[1];pos[1]=pos[0];pos[0]=ftell(fs);
          strncpy(line[2],line[1],511);
          strncpy(line[1],line[0],511);
          fgets(line[0],512,fs);

          /*check for file overrun*/
          if (feof(fs)){
              fprintf(stderr,"ERROR (Source_spectra): Could not strip header from file %s\n",fn);
              exit(-1);
          }

          int i;
          for (i=0;i<3;i++){
              linecheck[i]=sscanf(line[i],"%lf %lf %lf %lf %lf %lf", buf, buf+1, buf+2, buf+3, buf+4, buf+5);
          }
          if(linecheck[0]==linecheck[1] && linecheck[2]==linecheck[1] && (linecheck[0]==3 || linecheck[0]==5) ){
              done=1;
          }
      } while(!done);
      return pos[2];
  }

  double interpolate_4d_spectra(t_Table data, int *idx,int *N, double *alpha){
    /*interpolate in the 4d spectra data structure at the points idx;idx+1 according
      the normalized coordinates alpha*/
    double first[16],second[8],third[4],fourth[2];
    int i,j,k,l,m;
    double result;
    /*pick the "corner" values to interpolate in*/
    for (m=0;m<16;m++){
      i=idx[0]+ m%2;
      j=idx[1]+ (m/2)%2;
      k=idx[2]+ (m/4)%4;
      l=idx[3]+ (m/8)%2;
      first[m]=Table_Index(data,i+ j*N[0] + k*N[0]*N[1] + l*N[0]*N[1]*N[2],4);
    }
    /*reduce and interpolate the values by successive pairwise weighting, i.e. alpha*p{i,j,k,l} + (1-alpha*p_{i+1,j,k,l} etc.*/
    for (m=0;m<8;m++){
      second[m]=(1-alpha[0])*first[m*2] + alpha[0]*first[m*2+1];
    }
    for (m=0;m<4;m++){
      third[m]=(1-alpha[1])*second[m*2] + alpha[1]*second[m*2+1];
    }
    for (m=0;m<2;m++){
      fourth[m]=(1-alpha[2])*third[m*2] + alpha[2]*third[m*2+1];
    }
    result=(1-alpha[3])*fourth[0] + alpha[3]*fourth[1];
    return result;

  }

%}

DECLARE
%{
  double K;
  double dK;
  double pmul;
  double pint;
  t_Table *xproj;
  t_Table *yproj;
  t_Table *map;
  double *Ix;
  double *Iy;
  double *Imap;
  double xmin;
  double xmax;
  double ymin;
  double ymax;
  double xpmin;
  double xpmax;
  double ypmin;
  double ypmax;
  double xstep;
  double ystep;
  double xpstep;
  double ypstep;
  int brilliance_column;
%}

INITIALIZE
%{
  /*if noinit is set - return early from the init function*/
  if (noinit){
    return(_comp);
  }

  int num,status;
  long offset, orig_offset;
  char fnx[256]="";
  char fny[256]="";

  if(flag4d==0){
    /*flag4d not set - the datafiles are x,x',p and y,y',p projections*/
    brilliance_column=2;
    xproj=calloc(nE,sizeof(t_Table));
    yproj=calloc(nE,sizeof(t_Table));
    Ix=calloc(nE,sizeof(double));
    Iy=calloc(nE,sizeof(double));
    if(xproj==NULL || yproj==NULL || Ix==NULL || Iy==NULL){
      fprintf(stderr,"ERROR (%s): Memory allocation error\n",NAME_CURRENT_COMP);
      exit(-1);
    }
    /*find the offset of the datafiles. Assume to be identical for all of them.*/
    snprintf(fnx,255,"%s-%d.%s",spectra_stem_x,initial_serial,spectra_suffix);
    orig_offset=source_spectra_find_offset(fnx);

    for (num=0;num<nE;num++){
      snprintf(fnx,255,"%s-%d.%s",spectra_stem_x,num+initial_serial,spectra_suffix);
      offset=orig_offset;/*Have to do this every time since Table_Read_Offset overwrites offset*/
      if ( (status=Table_Read_Offset(&(xproj[num]),fnx,0,&offset,0))==-1){
        fprintf(stderr,"ERROR (%s): Could not parse file \"%s\"\n",NAME_CURRENT_COMP,fnx);
        exit(-1);
      }
      snprintf(fny,255,"%s-%d.%s",spectra_stem_y,num+initial_serial,spectra_suffix);
      offset=orig_offset;/*Have to do this every time since Table_Read_Offset overwrites offset*/
      if ( (status=Table_Read_Offset(&(yproj[num]),fny,0,&offset,0))==-1){
        fprintf(stderr,"ERROR (%s): Could not parse file \"%s\"\n",NAME_CURRENT_COMP,fny);
        exit(-1);
      }
      Ix[num]=Iy[num]=0;
      /*sum the brilliances to get something to normalize to*/
      int r;
      for (r=0;r<xproj[num].rows;r++){
        Ix[num]+=Table_Index(xproj[num],r,brilliance_column);//xproj.data[r*xproj.columns+ brilliance_column];
      }

      for (r=0;r<yproj[num].rows;r++){
        Iy[num]+=Table_Index(yproj[num],r,brilliance_column);//yproj.data[r*yproj.columns+ brilliance_column];
      }
      if (verbose && Ix[num]!=Iy[num]){
        fprintf(stderr,"WARNING (%s): Integrated intensities do not match up for x and y projections at num %d\n",NAME_CURRENT_COMP,num);
      }
      if (verbose) printf("INFO (%s): Integrated intensity for projections I [%d] = (%g,%g)\n",NAME_CURRENT_COMP,num,Ix[num],Iy[num]);

      if (num==0){
        /*check the data structure for the first two input files*/
        /*if not given deduce the number of sample-points in datafiles*/
        if (nx==0){
          int r;
          double p1,p2;
          for (r=0;r<xproj[0].rows;r++){
            if ( nx==0 && (p1=Table_Index(xproj[0],r,0))>(p2=Table_Index(xproj[0],r+1,0)) ){
              /*this means we have found where the first coordinate starts over*/
              nx=r+1;
              break;
            }
          }
          if (nx==0){
            nx=1;
          }
          npx/=nx;
        }
        if (npx==0){
          npx=xproj[0].rows/nx;
        }

        if( nx*npx != xproj[0].rows){
          fprintf(stderr,"Error (%s): number of read rows (%d) in %s does not match nx*npx = ( %d * %d). Please check the input files. Aborting.\n",NAME_CURRENT_COMP,xproj[0].rows,fnx,nx,npx);
          exit(-1);
        }

        if(ny==0){
          int r;
          for (r=0;r<yproj[0].rows;r++){
            if ( ny==0 && Table_Index(yproj[0],r,0)>Table_Index(yproj[0],r+1,0) ){
              /*this means we have found where the first coordinate starts over*/
              ny=r+1;
              break;
            }
          }
          if (ny==0){
            ny=1;
          }
        }
        if (npy==0){
          npy=yproj[0].rows/ny;
        }
        if( ny*npy != yproj[0].rows){
          fprintf(stderr,"ERROR (%s): number of read rows (%d) in %s does not match ny*npy = ( %d * %d). Please check the input files. Aborting.\n",NAME_CURRENT_COMP,yproj[0].rows,fny,ny,npy);
          exit(-1);
        }
        if(verbose) printf("INFO (%s): (nx,nxp) = ( %d %d ), (ny,npy) = ( %d %d )\n",NAME_CURRENT_COMP,nx,npx,ny,npy);
      }/*if num==0*/
    }
    /*find limits in x,x',y, and y', assuming they're the same across all source files.*/
    /*these would be relevant for a search*/

    t_Table *xptr=&(xproj[0]);
    t_Table *yptr=&(yproj[0]);
    xmin=Table_Index(*xptr,0,0);
    xpmin=Table_Index(*xptr,0,1);
    ymin=Table_Index(*yptr,0,0);
    ypmin=Table_Index(*yptr,0,1);
    xmax=Table_Index(*xptr,nx-1,0);
    xpmax=Table_Index(*xptr,nx*npx-1,1);
    ymax=Table_Index(*yptr,ny-1,0);
    ypmax=Table_Index(*yptr,ny*npy-1,1);
    xstep=Table_Index(*xptr,1,0)-Table_Index(*xptr,0,0);
    xpstep=Table_Index(*xptr,nx,1)-Table_Index(*xptr,0,1);
    ystep=Table_Index(*yptr,1,0)-Table_Index(*yptr,0,0);
    ypstep=Table_Index(*yptr,ny,1)-Table_Index(*yptr,0,1);

    if(verbose && xmin==0 && symmetricx==0){
    fprintf(stderr,"WARNING (%s): Minimum x-value in datafile is 0 but symmetricx is not set.\n",NAME_CURRENT_COMP);
    }
    if(verbose && xmin==0 && symmetricy==0){
      fprintf(stderr,"WARNING (%s): Minimum y-value in datafile is 0 but symmetricy is not set.\n",NAME_CURRENT_COMP);
    }
  } else {
    /*The datafiles are 4D: i.e. x,x',y,y' and p columns;*/
    brilliance_column=4;
    map=calloc(nE,sizeof(t_Table));
    Imap=calloc(nE,sizeof(double));
    sprintf(fnx,"%s-%d.%s",spectra_stem,initial_serial,spectra_suffix);
    orig_offset=source_spectra_find_offset(fnx);
    for (num=0;num<nE;num++){
      sprintf(fnx,"%s-%d.%s",spectra_stem,num+initial_serial,spectra_suffix);
      offset=orig_offset;/*Table_Read_Offset overwrites offset*/
      if ( (status=Table_Read_Offset(&(map[num]),fnx,0,&offset,0))==-1){
        fprintf(stderr,"Source_spectra(%s) Error: Could not parse file \"%s\"\n",NAME_CURRENT_COMP,fnx);
        exit(-1);
      }
      Imap[num]=0;
      /*sum the brilliances to get something to normalize to*/
      int r;
      for (r=0;r<map[num].rows;r++){
        Imap[num]+=Table_Index(map[num],r,brilliance_column);
      }
    }
    /*if the limits are not given, try to deduce the number of steps from the data files*/
    int r;
    double p1,p2;
    for (r=0;r<map[0].rows;r++){
      if ( nx==0 && (p1=Table_Index(map[0],r,0))>(p2=Table_Index(map[0],r+1,0)) ){
        /*this means we have found where the first coordinate starts over*/
        nx=r+1;
      }
      if ( ny==0 && (p1=Table_Index(map[0],r,1))>(p2=Table_Index(map[0],r+1,1)) ){
        /*this means we have found where the 2nd coordinate starts over*/
        ny=(r+1)/nx;
      }
      if ( npx==0 && (p1=Table_Index(map[0],r,2))>(p2=Table_Index(map[0],r+1,2)) ){
        /*this means we have found where the 3rd coordinate starts over*/
        npx=(r+1)/(nx*ny);
      }
      /*if first 3 are set stop searching - set the last later*/
      if(nx && ny && npx ){
        break;
      }
    }
    if (nx==0){
      nx=1;
    }
    if (ny==0){
      ny=1;
    }
    if (npx==0){
      npx=1;
    }
    /*last one is set frm th evalues of the other ones*/
    if (npy==0){
      npy=map[0].rows/(nx*ny*npx);
    }
    if (verbose) printf("INFO (%s): [Nx,Ny,Nx,Ny]= [ %d %d %d %d ] for files: %s-X.%s\n",NAME_CURRENT_COMP,nx,ny,npx,npy,spectra_stem,spectra_suffix);

    /*set limit values. Use that the last row contains the maximum value for all 4 dimensions.*/
    xmin=Table_Index(map[0],0,0);
    xmax=Table_Index(map[0],map[0].rows,0);
    xstep=(xmax-xmin)/nx;
    ymin=Table_Index(map[0],0,1);
    ymax=Table_Index(map[0],map[0].rows,1);
    ystep=(ymax-ymin)/ny;
    xpmin=Table_Index(map[0],0,2);
    xpmax=Table_Index(map[0],map[0].rows,2);
    xpstep=(xpmax-xpmin)/npx;
    ypmin=Table_Index(map[0],0,3);
    ypmax=Table_Index(map[0],map[0].rows,3);
    ypstep=(ypmax-ypmin)/npy;

    if (verbose) printf("INFO (%s): (xmin,xmax)=( %g %g ), (ymin,ymax)=( %g %g ), (xpmin,xpmax)=( %g %g ), (ypmin,ypmax)=( %g %g )\n",NAME_CURRENT_COMP,xmin,xmax,ymin,ymax,xpmin,xpmax,ypmin,ypmax);
  }

  if(E0-dE-Emin <-FLT_MAX || E0+dE-Emax>FLT_MAX){
    fprintf(stderr,"WARNING(%s): Sampled energy interval (%g+-%g keV) reaches outside what\'s defined by datafiles (%g+-%g keV)\n",NAME_CURRENT_COMP,E0,dE,(Emin+Emax)*0.5,(Emax-Emin)*0.5);
  }

  /*downweight accoring to number of rays*/
  pmul=1.0/((double) mcget_ncount());

  /*downweight for not using the full energy window. Don't do this for deterministic energy (dE=0).*/
  if(dE && ( (E0-dE>Emin) || E0+dE<Emax) ){
    pmul*=2*dE/(Emax-Emin);
  }

%}

TRACE
%{
  double kk,theta_x,theta_y,l,e,k,xp,yp;
  int num,ix,ipx,iy,ipy;
  double alpha,beta,Iinterpx,Iinterpy;
  t_Table *xptr,*yptr;
  
  p=pmul;
  theta_x=(xpmin + rand01()*(xpmax-xpmin))*1e-3;
  theta_y=(ypmin + rand01()*(ypmax-ypmin))*1e-3;
  
  x=(xmin+rand01()*(xmax-xmin))*1e-3;
  y=(ymin+rand01()*(ymax-ymin))*1e-3;

  
  /*So now interpolate to get at Brilliance values*/
  /*Need to normalize to something*/
  
  /*pick an energy randomly*/
  e=rand01()*2*dE+(E0-dE);
  if(e<Emin || e>Emax){
    ABSORB;
  }
  k=E2K*e;

  kx=tan(theta_x);
  ky=tan(theta_y);
  kz=1;
  NORM(kx,ky,kz);

  kx*=k;
  ky*=k;
  kz*=k;
  /*compute xp and yp*/
  xp=kx/kz*1e3;/*spectra output is in millirad*/
  yp=ky/kz*1e3;
  double xx=x*1e3;/*spectra output is in mm*/
  double yy=y*1e3;

  ix  = (int)floor((xx - xmin)*(nx-1)/(xmax - xmin));
  ipx = (int)floor((xp- xpmin)*(npx-1)/(xpmax-xpmin));
  iy  = (int)floor((yy - ymin)*(ny-1)/(ymax - ymin));
  ipy = (int)floor((yp- ypmin)*(npy-1)/(ypmax-ypmin));

  int ie;
  double pe[2];
  double ealpha,estep;
  estep=(Emax-Emin)/(nE-1);
  num=(int)floor( (e-Emin)/estep);
  if (num<0) num=0;
  if (num>nE-1) num=nE-1;
  ealpha = (e- (Emin+estep*num))/estep;

  if(!flag4d){
    xptr=&(xproj[num]);
    yptr=&(yproj[num]);
    for (ie=0;ie<2;ie++){
      alpha=( (xx - Table_Index(*xptr,ix,0)) /xstep);/*regular grid so no need to do ix + ipx*nx*/
      beta=( (xp - Table_Index(*xptr,ipx*nx,1)) /xpstep) ;

      double t0,t1;
      t0=(1-alpha)*fabs(Table_Index(*xptr,ix+ipx*nx,2)) + alpha*fabs(Table_Index(*xptr,(ix+1)+ipx*nx,2));
      t1=(1-alpha)*fabs(Table_Index(*xptr,ix+(ipx+1)*nx,2)) + alpha*fabs(Table_Index(*xptr,(ix+1)+(ipx+1)*nx,2));
      Iinterpx = (1-beta)*t0+beta*t1 *xstep*1e-3*xpstep*1e-3;

      alpha=( (yy - Table_Index(*yptr,iy,0)) /ystep);/*regular grid so no need to do iy + ipy*ny*/
      beta=( (yp - Table_Index(*yptr,ipy*ny,1)) /ypstep) ;

      t0=(1-alpha)*fabs(Table_Index(*yptr,iy+ipy*ny,2)) + alpha*fabs(Table_Index(*yptr,(iy+1)+ipy*ny,2));
      t1=(1-alpha)*fabs(Table_Index(*yptr,iy+(ipy+1)*ny,2)) + alpha*fabs(Table_Index(*yptr,(iy+1)+(ipy+1)*ny,2));
      Iinterpy = (1-beta)*t0+beta*t1 * ystep*1e-3*ypstep*1e-3;

      pe[ie]=Iinterpx/Ix[num+ie] * Iinterpy/Iy[num+ie] * (Ix[num+ie]+Iy[num+ie])*0.5;
    }
    /* Set the photon ray weight to the mean of the two multiplied by initial weight
     * due to energy interval and ray count*/

    //p=pmul * Iinterpx/Ix[num] * Iinterpy/Iy[num] * (Ix[num]+Iy[num])*0.5;
    /*now we interpolate in energy*/
    p=pmul * ( (1-ealpha)*pe[0] + ealpha*pe[1]);
  }else{
    /*each pt in x,x',y,y' space is surrounded by 16 points. Interpolate between them*/
    double alx,alxp,aly,alyp;
    alx=( (xx - Table_Index(map[ie],ix,0)) /xstep);
    aly=( (yy - Table_Index(map[ie],iy*nx,1)) /ystep);
    alxp=( (xp - Table_Index(map[ie],ipx*nx*ny,2)) /xpstep);
    alyp=( (yp - Table_Index(map[ie],ipy*nx*ny*npx,3)) /ypstep);

    for (ie=0;ie<2;ie++){
      int idx[]={ix,iy,ipx,ipy};
      double alpha[]={alx,aly,alxp,alyp};
      int NN[]={nx,ny,npx,npy};
      pe[ie]=interpolate_4d_spectra(map[num+ie],idx,NN,alpha)*xstep*ystep*1e-6 *xpstep*ypstep*1e-6;
    }
    /*lastly interpolate over energy*/
    p=pmul * ( (1-ealpha)*pe[0] + ealpha*pe[1]);
  }
  /*if symmetric source possibly reflect x*/
  if (symmetricx){
    if (rand01()<0.5){
      x=-x;
    }
    if (rand01()<0.5){
      kx=-kx;
    }
  }
  if (symmetricy){
    if (rand01()<0.5){
      y=-y;
    }
    if (rand01()<0.5){
      ky=-ky;
    }
  }

  /*spectra output is in brilliance (unit photons/s/mm^2/mrad^2/0.1%BW), so scale to get at raw flux in photons/s */

  /*set polarization and phase to something known*/
  Ex=0;Ey=0;Ez=0;
  if (!randomphase){
    phi=0;
  }else{
    phi=rand01()*M_2_PI;
  }

  /*set polarization vector*/
  Ex=0;Ey=0;Ez=0;

%}

FINALLY
%{
  if(!noinit){
    if(!flag4d){
      free(Ix);
      free(Iy);
      Table_Free(xproj);
      Table_Free(yproj);
      free(yproj);
      free(xproj);
    }else{
      free(Imap);
      Table_Free(map);
      free(map);
    }
  }
%}

MCDISPLAY
%{
  double dist=1;
  multiline(5, xmin, ymin, 0.0,
      xmax, ymin, 0.0,
      xmax, ymax, 0.0,
      xmin, ymax, 0.0,
      ymin, ymin, 0.0);

  dashed_line(0,0,0, tan(xpmax)*dist,0,dist,4);
  dashed_line(0,0,0, tan(xpmin)*dist,0,dist,4);
    
  dashed_line(0,0,0, 0,tan(ypmax)*dist,dist,4);
  dashed_line(0,0,0, 0,tan(ypmin)*dist,dist,4);
%}

END