/*******************************************************************************
*
* McStas, the neutron ray-tracing package
*         Maintained by Kristian Nielsen and Kim Lefmann,
*         Copyright 1997-2000 Risoe National Laboratory, Roskilde, Denmark
*
*Component: SNS_source
*
* %I
* Written by: G. Granroth
* Date: July 2004
* Origin: SNS Project Oak Ridge National Laboratory
* Modified: G. E. Granroth  Nov 2014 Move to Mcstas V2
* Modified: J.Y.Y. Lin April 2021 to Fix race condition
* Modified: P. WIllendrup 2021 Move to Move to Mcstas Version 3 and enable GPus
*
* A source that produces a time and energy distribution from the SNS moderator files
*
* %D
* Produces a time-of-flight spectrum from SNS moderator files
* moderator files can be obtained from the <A href="http://www.sns.gov/users/instrument_systems/components/moderator/">SNS website </A>.
* <strong>IMPORTANT: The output units of this component are N/pulse</strong>
* <strong>IMPORTANT: The component needs a FULL PATH to the source input file</strong>
* Notes:
* (1) the raw moderator files are per Sr.  The focusing parameters provide the solid
* angle accepted by the guide to remove the Sr dependence from the output. Therefore
* the best practice is to set focus_xw and focus_yh to the width and height of the next beam
* component, respectively. The dist parameter should then be set as the distance
* from the moderator to the first component.
* (2) This component works purely by interpolation. Therefore be sure that Emin and
* Emax are within the limits of the moderator file
*
*
* %P
* Input parameters:
* filename: []   Filename of source data
* xwidth: [m]    width of moderator
* yheight: [m]   height of moderator
* dist: [m]      Distance from source to the focusing rectangle
* xw: [m]        Width of focusing rectangle
* focus_yh: [m]  Height of focusing rectangle
* focus_xw: [m]  Width of focusing rectangle
* Emin: [meV]    minimum energy of neutron to generate
* Emax: [meV]    maximum energy of neutron to generate
* %E
*******************************************************************************/

DEFINE COMPONENT SNS_source

SETTING PARAMETERS (string filename, xwidth=0.1, yheight=0.12, dist, focus_xw, focus_yh, Emin, Emax)


SHARE %{

  #define Maxlength 200
  #define MAXCOLS 500

  /* ----------------------------------------------------------------
      routine to load E, I and t I data from SNS source files
  -----------------------------------------------------------------*/
  void sns_source_load(char* fname, double *xvec, double *yvec, int xcol, int ycol,
    int *veclenptr, double *tcol, double *Ecol, double **Imat,int *ntvals, int *nEvals)
  {
    FILE *fp;
    int idx1,idx2,idx3; /* counter for number of x, y values */
    int jk,idx3max;
    int numtvals;
    int totalvals;
    float indat[6];
    double *Icoltmp, *tcoltmp, *Ecoltmp;
    char *line;
    char *ret;
    Icoltmp=malloc(100000*sizeof(double));
    tcoltmp=malloc(100000*sizeof(double));
    Ecoltmp=malloc(100000*sizeof(double));
    line=malloc(200*sizeof(char));
    /* open file */
    printf("%s\n",fname);
    fp=fopen(fname,"r");
    if (fp==NULL){
      exit(printf("Error opening file: %s. Check existence/permission. Aborting.\n", fname));
      }
    /* skip header lines any line that begin with # */
    while((fgets(line,Maxlength,fp)!=NULL)&&(strchr(line,'#')!=NULL)){
    }
    idx1=0;
    /* read all lines that fit the format for the time integrated data*/
    while(sscanf(line," %f %f %f %f %f %f",&indat[0], &indat[1], &indat[2], &indat[3],&indat[4],&indat[5])>0){
        xvec[idx1]=indat[xcol];
        yvec[idx1]=indat[ycol];
        //printf("idx1 %i xvec %g yvec %g\n",idx1,xvec[idx1],yvec[idx1]);
        idx1++;
        ret=fgets(line,Maxlength,fp);
    }
    idx1--;   // correct index since it counts one line past useful data
   // printf("idx1 %i\n",idx1);
    idx2=floor(idx1/2);
    while((idx2<idx1)&&(yvec[idx2]>0)){
      idx2++;
    }
    if(idx2<idx1){
      *veclenptr=idx2;
    }
    else{
      *veclenptr=idx1-2;
    }
  /* find t data header */
    ret=fgets(line,Maxlength,fp);
    while(strrchr(line,'#')==NULL){
      ret=fgets(line,Maxlength,fp);
    }
  /*find end of t data header */
    while((fgets(line,Maxlength,fp)!=NULL)&&(strchr(line,'#')!=NULL)){
    }
  /* read in t data */
 /*printf("t data read start\n");*/
    idx2=0;
   // while(fgets(line,Maxlength,fp)!=NULL){
   do {
        jk=sscanf(line," %f %f %f %f",&indat[0], &indat[1], &indat[2], &indat[3]);
          if ((line!=NULL)&&(jk>3)){
            tcoltmp[idx2]=indat[0];
            Ecoltmp[idx2]=indat[1];
            Icoltmp[idx2]=indat[2];
           // printf("%d %d %g %g %g %g\n",idx2,jk,tcoltmp[idx2],Ecoltmp[idx2],Icoltmp[idx2],indat[3]);
            idx2++;
          }
    }
    while(fgets(line,Maxlength,fp)!=NULL);
    fclose(fp);
    totalvals=idx2+1;
    printf("total vals: %d\n",totalvals);
    /* reformat data into an Ecol, a tcol, and an I matrix*/
    idx1=0;idx2=0;idx3=0;
    Ecol[idx2]=Ecoltmp[idx1];
    tcol[idx3]=tcoltmp[idx1];
    Imat[idx3][idx2]=Icoltmp[idx1];
    idx1++;idx3++;
    while(idx1<totalvals){
      jk=idx1-1;
      if(Ecoltmp[idx1]==Ecoltmp[jk]){
        tcol[idx3]=tcoltmp[idx1];
         Imat[idx3][idx2]=Icoltmp[idx1];
         idx1++;idx3++;
      }
      else{
        idx3max=idx3;
        idx2++;idx3=0;
        Ecol[idx2]=Ecoltmp[idx1];
        tcol[idx3]=tcoltmp[idx1];
        Imat[idx3][idx2]=Icoltmp[idx1];
        idx1++;
      }
    }
   *ntvals=idx3max+1;*nEvals=idx2;
   printf("ntvals: %i idx: %i\n",*ntvals, idx3);
   free(Icoltmp);free(tcoltmp);free(Ecoltmp);free(line);
  } /* sns_source_load */
/*-------------------------------------------------------------
        End load file routines
--------------------------------------------------------------*/

/*----------------------------------------------------------------------
                routine to do a 1D linear interpolation
------------------------------------------------------------------------*/
/* given a point (x1,y1) on the low side of xdes and one (x2,y2) on the
high side of xdes return the interpolated y values */
#pragma acc routine
double linint(double xdes,double x1, double x2, double y1, double y2)
{
  double m;
  m=(y2-y1)/(x2-x1);
  return (m*(xdes-x1)+y1);
} /*linint */

#pragma acc routine
double linfuncint(double xdes, double xylen, double *vecx, double *vecy)
  {
  int idx;
  idx=0;
  while((vecx[idx]<xdes)&&idx<xylen){
     idx++;
    }
  if (idx>xylen){
      printf("SNS_source: linfuncint: error exceeded vector length");
    }
  if (vecx[idx]==xdes){
      return vecy[idx];
      }
   else
     {
     return linint(xdes,vecx[idx-1],vecx[idx],vecy[idx-1],vecy[idx]);
     }

  } /* linfuncint */
/*------------------------------------------------------------------------
            routine to perform a 1 d quadratic interpolation
--------------------------------------------------------------------------*/
/* given 2 points on the low side of xdes and one on the high side, return
a quadratically interpolated result */
#pragma acc routine
double quadint(double xdes,double x1, double x2,double x3, double y1, double
y2, double y3)
{
  double t1, t2, t3;
  t1=((xdes-x2)*(xdes-x3)*y1)/((x1-x2)*(x1-x3));
  t2=((xdes-x1)*(xdes-x3)*y2)/((x2-x1)*(x2-x3));
  t3=((xdes-x1)*(xdes-x2)*y3)/((x3-x1)*(x3-x2));
  return t1+t2+t3;
} /* quadint */

#pragma acc routine
double quadfuncint(double xdes, double xylen, double *vecx, double *vecy)
  {
  int idx;
  idx=1;
  while((vecx[idx]<xdes)&&idx<xylen){
     idx++;
    }
  if (idx>xylen){
      printf("SNS_source: quadfuncint: error exceeded vector length");
    }

   if (vecx[idx]==xdes){
      return vecy[idx];
      }
   else
     {
     return quadint(xdes,vecx[idx-2],vecx[idx-1],vecx[idx],vecy[idx-2],vecy[idx-1],vecy[idx]);
     }

  } /* quadfuncint */

/*-----------------------------------------------------------------
Functions for random energy generation
------------------------------------------------------------------*/
#pragma acc routine
double xonly(double x,double xylength,double *inxvec,double *inyvec)
{
  return linfuncint(x,xylength,inxvec,inyvec);
}
#pragma acc routine
double Pfunc(double x, double y,double xylength,double *inxvec,double *Pvec)
{
  return quadfuncint(x,xylength,inxvec,Pvec)-y;
}
  
/*----------------------------------------------------------------
Functions for random time generation
------------------------------------------------------------------*/
#pragma acc routine
double txonly(double t,double ntvals,double *txval,double *tyval)
{
  return linfuncint(t,ntvals,txval,tyval);
}
#pragma acc routine
double tPfunc(double t,double y,double ntvals,double *txval,double *tyval)
{
  return quadfuncint(t,ntvals,txval,tyval)-y;
}

 
/*-------------------------------------------------------------------
  integration routines
---------------------------------------------------------------------*/
double integtrap(double (*func)(double,double,double*,double*),double prev,double low,double high, int step,
  double len, double *xvec, double *yvec)
{
 double s,npts,stpsze,sum,x;
 int pw2, idx;
 if (step==1){
   return(s=0.5*(high-low)*((*func)(high, len, xvec,yvec)+(*func)(low, len, xvec,yvec)));
 }
 else{
   s=prev;
   for(pw2=1,idx=1;idx<step-1;idx++){
     pw2<<=1;
   }
   npts=pw2;
   stpsze=(high-low)/npts;
   x=low+0.5*stpsze;
   for(sum=0.0,idx=1;idx<=pw2;idx++,x+=stpsze){
     sum+=(*func)(x, len, xvec,yvec);
   }
   s=0.5*(s+(high-low)*sum/npts);
   return s;
 }
} /* integtrap */

double integ1(double (*func)(double,double,double*,double*),double low, double high, double err,
  double len, double *xvec, double *yvec)
{
  double out,outprev;
  int idx;
  out=integtrap(*func,0.0,low,high,1, len, xvec, yvec);
  outprev=out;
  out=integtrap(*func,out,low,high,2, len, xvec, yvec);
  /*printf("out %g outprev %g \n",out,outprev);*/
  idx=2;
  while(fabs(out-outprev)>err*fabs(out)){
    idx++;
    outprev=out;
    out=integtrap(*func,out,low,high,idx, len, xvec, yvec);
   /* printf("out %g outprev %g \n",out,outprev);*/
  }
  return out;
}
/*---------------------------------------------------------------------------
   Routine for finding zeros.
  Algorithm equivalent to rtbis from "Numerical Recipes in C: pg 354
 -----------------------------------------------------------------------------*/


double zero_find(double (*func)(double, double, double, double*, double*),
  double yval,double xmin,double xmax, double tol, double arg1, double *arg2, double *arg3)
{
  double xl,xh,f,fmid,xmid,dx,rtb;
  int max_iter;
  int iter_no;
  xl=xmin;
  xh=pow(10,(log10(xmin)+yval*(log10(xmax)-log10(xmin))));
  f=(*func)(xl,yval, arg1, arg2, arg3);
  fmid=(*func)(xh,yval, arg1, arg2, arg3);
  while (fmid*f>=0.0){
    xh=xh+(xh-xl)*2.0;
    fmid=(*func)(xh,yval, arg1, arg2, arg3);
  }
  dx=xh-xl;
  rtb=xl;
  max_iter = 10000;
  iter_no=0;
  while(fabs((*func)(rtb,yval, arg1, arg2, arg3))>tol && iter_no++<max_iter){
    dx=dx*0.5;
    xmid=rtb+dx;
    fmid=(*func)(xmid,yval, arg1, arg2, arg3);
    if (fmid<0){
      rtb=xmid;
    }
  }
 return rtb;
}

#pragma acc routine
double zero_find_gpu1(double yval,double xmin,double xmax, double tol, double arg1, double *arg2, double *arg3)
{
  double xl,xh,f,fmid,xmid,dx,rtb;
  int max_iter;
  int iter_no;
  xl=xmin;
  xh=pow(10,(log10(xmin)+yval*(log10(xmax)-log10(xmin))));
  f=Pfunc(xl,yval, arg1, arg2, arg3);
  fmid=Pfunc(xh,yval, arg1, arg2, arg3);
  while (fmid*f>=0.0){
    xh=xh+(xh-xl)*2.0;
    fmid=Pfunc(xh,yval, arg1, arg2, arg3);
  }
  dx=xh-xl;
  rtb=xl;
  max_iter = 10000;
  iter_no=0;
  while(fabs(Pfunc(rtb,yval, arg1, arg2, arg3))>tol && iter_no++<max_iter){
    dx=dx*0.5;
    xmid=rtb+dx;
    fmid=Pfunc(xmid,yval, arg1, arg2, arg3);
    if (fmid<0){
      rtb=xmid;
    }
  }
 return rtb;
}
#pragma acc routine
double zero_find_gpu2(double yval,double xmin,double xmax, double tol, double arg1, double *arg2, double *arg3)
{
  double xl,xh,f,fmid,xmid,dx,rtb;
  int max_iter;
  int iter_no;
  xl=xmin;
  xh=pow(10,(log10(xmin)+yval*(log10(xmax)-log10(xmin))));
  f=tPfunc(xl,yval, arg1, arg2, arg3);
  fmid=tPfunc(xh,yval, arg1, arg2, arg3);
  while (fmid*f>=0.0){
    xh=xh+(xh-xl)*2.0;
    fmid=tPfunc(xh,yval, arg1, arg2, arg3);
  }
  dx=xh-xl;
  rtb=xl;
  max_iter = 10000;
  iter_no=0;
  while(fabs(tPfunc(rtb,yval, arg1, arg2, arg3))>tol && iter_no++<max_iter){
    dx=dx*0.5;
    xmid=rtb+dx;
    fmid=tPfunc(xmid,yval, arg1, arg2, arg3);
    if (fmid<0){
      rtb=xmid;
    }
  }
 return rtb;
}
/*----------------------------------------------------------------------------
Routine for calculating Probability distribution
----------------------------------------------------------------------------*/

// Pcalc(xonly,llim,hlim,inxvec,Pvec,xylength,&idxstart,&idxstop, inyvec);
// xonly(double x,double xylength,double *inxvec,double *inyvec)
void Pcalc(double (*func)(double,double,double*,double*),
  double llim, double hlim, double *xvec, double *Prob, int veclen,
  int *idxstart, int *idxstop, double *inyvec)
{
  int idx1,idx2;
  double junk,Norm;
  idx1=0;
  while(xvec[idx1]<=llim){
     Prob[idx1]=0;
     idx1++;

  }
  if (idx1<1){
   printf("Error: lower energy limit is out of bounds\n");
   exit(0);
  }
  *idxstart=idx1;
   Prob[idx1]=integ1((*func),llim,xvec[idx1],0.001, veclen, xvec, inyvec);
   idx1++;
  while(xvec[idx1]<=hlim){
     junk=integ1((*func),xvec[idx1-1],xvec[idx1],0.001, veclen, xvec, inyvec);
     Prob[idx1]=(Prob[idx1-1]+junk);
     idx1++;
   }
   *idxstop=idx1;
   while(idx1<veclen){
     Prob[idx1]=1;
     idx1++;
   }

  /*Normalize all Probability values*/
    Norm=Prob[*idxstop-1];
  if (Norm>0){
    for(idx2=*idxstart;idx2<*idxstop;idx2++){
      Prob[idx2]=Prob[idx2]/Norm;
    }
  }
} /* Pcalc */
/*----------------------------------------------------------------------------
Routine for calculating t Probability distribution
----------------------------------------------------------------------------*/
// tPcalc(txonly,ltlim,htlim,tcol,tPvec,ntvals,&tidxstart,&tidxstop, txval,tyval);
// txonly(double t,double ntvals,double *txval,double *tyval)
void tPcalc(double (*func)(double,double,double*,double*),double llim, double hlim,
  double *xvec, double *Prob, int veclen, int *idxstart, int *idxstop,
  double *txval,double *tyval)
{
  int idx1,idx2;
  double junk,Norm;
  idx1=0;
  while(xvec[idx1]<=llim){
     Prob[idx1]=0;
     idx1++;

  }
  *idxstart=idx1;
   Prob[idx1]=integ1((*func),llim,xvec[idx1],0.001, veclen, txval, tyval);
  while(xvec[idx1]<=hlim){
     junk=integ1((*func),xvec[idx1-1],xvec[idx1],0.001, veclen, txval, tyval);
     Prob[idx1]=(Prob[idx1-1]+junk);
     idx1++;
   }
   *idxstop=idx1;
   while(idx1<veclen){
     Prob[idx1]=1;
     idx1++;
   }
   /* calculate normalization*/
   Norm=Prob[*idxstop-1];
   /*printf("Norm %f\n",Norm); */
  /*Normalize all Probability values*/
  if (Norm>0){
    for(idx2=*idxstart;idx2<*idxstop;idx2++){
      Prob[idx2]=Prob[idx2]/Norm;
      /*printf("%d %g \n",idx2,Prob[idx2])*/;
    }
  }
} /* tPcalc */

%}

DECLARE
%{
  double p_in;
  double *inxvec;
  #pragma acc shape(inxvec[0:500])
  double *inyvec;
  #pragma acc shape(inyvec[0:500])
  double *Pvec;
  #pragma acc shape(Pvec[0:500])
  int xylength;
  double *tcol;
  #pragma acc shape(tcol[0:200])
  double *Ecol;
  #pragma acc shape(Ecol[0:200])
  double *tPvec;
  #pragma acc shape(tPvec[0:500])
  double **Ptmat;
  #pragma acc shape(Ptmat[0:200][0:200])
  double EPmax;
  double EPmin;
  double INorm;
  double INorm2;
  int ntvals;
  int idxstart;
  int idxstop;
  int tidxstart;
  int tidxstop;
  int nEvals;
  double pmul;
%}

INITIALIZE
%{
  FILE *fp;
  double llim, hlim,ltlim,htlim,junk;
  double **Imat;
  int idx1,idx2;

  double tyval[500];
  double txval[500];

  Pvec=malloc(500*sizeof(double));
  inxvec=malloc(500*sizeof(double));
  inyvec=malloc(500*sizeof(double));
  tcol=malloc(200*sizeof(double));
  Ecol=malloc(200*sizeof(double));
  tPvec=malloc(500*sizeof(double));
  Ptmat=malloc(200*sizeof(double *));
  for(idx1=0;idx1<200;idx1++){
    Ptmat[idx1]=malloc(200*sizeof(double));
  }
  Imat=malloc(200*sizeof(double*));
  for(idx1=0;idx1<200;idx1++){
    Imat[idx1]=malloc(500*sizeof(double));
  }
  ltlim=0.1;
  htlim=1.8e3;
 /* read file */
  printf("%s%s\n","Loading moderator file ",filename);
   sns_source_load(filename,inxvec,inyvec,0,2,&xylength,tcol,Ecol,Imat,&ntvals,&nEvals);
 /* calculate probabilty distribution function points for use in interpolation routine */

 llim=inxvec[1];hlim=inxvec[xylength];
 printf("Start calculating probability distribution\n");
 /* calculate total number of neutrons in specified energy window */
 INorm2=integ1(xonly,Emin/1000.0,Emax/1000.0,0.001, xylength,inxvec,inyvec);
  Pcalc(xonly,llim,hlim,inxvec,Pvec,xylength,&idxstart,&idxstop,
    inyvec);
  /*calculate probability distribution as a function of t for each energy value */
  tyval[0]=Imat[0][0];
  //printf("outntvals %i\n",ntvals);
  //printf("%g \n",tyval[0]);
  for(idx1=0;idx1<nEvals;idx1++){
      for(idx2=0;idx2<ntvals;idx2++){
        //printf("idx2 %d idx1 %d tyval: %g\n",idx2,idx1,Imat[idx2][idx1]);
        tyval[idx2]=Imat[idx2][idx1];
        txval[idx2]=tcol[idx2];
      }
      tPcalc(txonly,ltlim,htlim,tcol,tPvec,ntvals,&tidxstart,&tidxstop,
        txval,tyval);
      for(idx2=0;idx2<ntvals;idx2++){
        Ptmat[idx2][idx1]=tPvec[idx2];
      }
   }
  /* Generate min and max Prob values corresponding to specified energy window*/
    /* use quadratic interpolation  */
   EPmax=quadfuncint(Emax/1000.0,xylength,inxvec,Pvec);
   EPmin=quadfuncint(Emin/1000.0,xylength,inxvec,Pvec);
  /* Calculate solid Angle */
  p_in = focus_xw*focus_yh/(dist*dist);
  for(idx1=0;idx1<200;idx1++){
    free(Imat[idx1]);
  }
  free(Imat);

  printf("Finished calculating probability distribution\n");
  pmul=1;
  pmul/=mcget_ncount();
%}
TRACE
%{
  double theta,phi,v,tauv,E,Eval,tval,randp;  
  double hdiv;
  double vdiv;
  double hdivmin,hdivmax,vdivmin,vdivmax;
  int idx1,idx3,Eidxl,Eidxh;
  double tyval[500];
  double txval[500];
  p=p_in;
  z=0;

  x = (rand01()-0.5)*xwidth; /* choose points uniformly distributed on the source */
  y = (rand01()-0.5)*yheight;

 hdivmax=atan((focus_xw/2.0-x)/dist);
 hdivmin=atan(-(focus_xw/2.0+x)/dist);
 vdivmax=atan((focus_yh/2.0-y)/dist);
 vdivmin=atan(-(focus_yh/2.0+y)/dist);

  theta = hdivmin + (hdivmax-hdivmin)*rand01(); /* Small angle approx. */
  phi = vdivmin + (vdivmax-vdivmin)*rand01();
  hdiv=theta;
  vdiv=phi;
  /* generate random numbers*/

   /*First generate E random value */
    randp=EPmin+rand01()*(EPmax-EPmin);
    /* find E value corresponding to random probability */
    #ifndef OPENACC
    Eval=zero_find(Pfunc,randp,inxvec[idxstart],inxvec[idxstop],1e-7,
      xylength,inxvec,Pvec);
    #else
    Eval=zero_find_gpu1(randp,inxvec[idxstart],inxvec[idxstop],1e-7,
      xylength,inxvec,Pvec);
    #endif
   /* from a known E value generate an emission time value */
   /* find the index of the E values that bracket the random E value */
   idx1=0;
   while((idx1<nEvals)&&(Ecol[idx1]<Eval)){
     idx1++;
   }
   Eidxh=idx1;
   Eidxl=idx1-1;
   /*calculate the interpolated t, P curve to determine t value from uniform random value*/
   for (idx3=0;idx3<ntvals;idx3++){
     tyval[idx3]=linint(Eval,Ecol[Eidxl],Ecol[Eidxh],Ptmat[idx3][Eidxl],Ptmat[idx3][Eidxh]);
     txval[idx3]=tcol[idx3];
   }
   randp=tyval[tidxstart]+rand01()*(tyval[tidxstop-1]-tyval[tidxstart]);
   if (randp>0.0){
     #ifndef OPENACC
     tval=zero_find(tPfunc,randp,txval[tidxstart],txval[tidxstop-1],1e-5,
        ntvals,txval,tyval); }
     #else
     tval=zero_find_gpu2(randp,txval[tidxstart],txval[tidxstop-1],1e-5,
        ntvals,txval,tyval); }
     #endif
    else{
       tval=0;}
  E = Eval*1000.0;  /* Convert Energy from Ev to meV */
  t = tval*1e-6;      /* Convert time from mus to S */
  v = SE2V*sqrt(E);
  /* Calculate components of velocity vector such that the neutron is within the focusing rectangle */
  vz = v*cos(phi)*cos(theta);   /* Small angle approx. */
  vy = v*sin(phi);
  vx = v*cos(phi)*sin(theta);

  p*=(xwidth*yheight/(0.1*0.12))*INorm2*pmul;

%}
FINALLY
%{
  int idxf;
  free(tPvec);
  free(inxvec);free(inyvec);free(Pvec);free(tcol);free(Ecol);
  for(idxf=0;idxf<200;idxf++){
    free(Ptmat[idxf]);
  }
  free(Ptmat);

%}

MCDISPLAY
%{
  double x1,y1,x2,y2;
  x1=-xwidth/2.0;y1=-yheight/2.0;x2=xwidth/2.0;y2=yheight/2.0;
  multiline(4,(double)x1,(double)y1,0.0,(double)x1,(double)y2,0.0,(double)x2,(double)y2,0.0,(double)x2,(double)y1,0.0,(double)x1,(double)y1,0.0);
%}

END