/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2002, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: SNS_source
*
* %I
* Written by: G. Granroth
* Date: July 2004
* Version: $Revision: 1.4 $
* Origin: SNS Project Oak Ridge National Laboratory
*
* 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>.
* The output units of this component are N/pulse
* 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 xw and 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:
* S_filename:   Filename of source data (str)
* width:  (m)   width of moderator
* height: (m)   height of moderator
* dist:   (m)   Distance from source to the focusing rectangle
* xw:     (m)   Width of focusing rectangle
* yh:     (m)   Height of focusing rectangle
* Emin: (meV)   minimum energy of neutron to generate
* Emax: (meV)   maximum energy of neutron to generate
*
* %L
* <A href="http://www.sns.gov/users/instrument_systems/components/moderator/">Source files </A>
* SNS_test example instrument
* %E
*******************************************************************************/

DEFINE COMPONENT SNS_source
DEFINITION PARAMETERS (S_filename)
SETTING PARAMETERS (width, height, dist, xw, yh, Emin, Emax)
OUTPUT PARAMETERS (hdiv,vdiv,p_in,
  inxvec,inyvec,Pvec,xylength,tcol,Ecol,txval,tyval,tPvec,Ptmat,
  EPmax, EPmin,INorm,INorm,
  ntvals,idxstart,idxstop,tidxstart,tidxstop,nEvals,
  xonly,Pfunc,txonly,tPfunc)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
DECLARE
%{
  double xonly(double);
  double Pfunc(double,double);
  double txonly(double);
  double tPfunc(double,double);
  double hdiv,vdiv;
  double p_in;
  double *inxvec,*inyvec,*Pvec;
  int xylength;
  double *tcol, *Ecol;
  double *txval, *tyval;
  double *tPvec;
  double **Ptmat;
  double EPmax, EPmin,INorm,INorm2;
  int ntvals,idxstart,idxstop,tidxstart,tidxstop,nEvals;
  #define Maxlength 200
  #define MAXCOLS 500
  /* ----------------------------------------------------------------
      routine to load E, I and t I data from SNS source files
  -----------------------------------------------------------------*/

  void sns_source_load(char filename[], 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 *lntoken, *cp;
    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",filename);
    fp=fopen(filename,"r");
    if (fp==NULL){
      printf("Error opening file");
      }
    /* 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++;
        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 */
    fgets(line,Maxlength,fp);
    while(strrchr(line,'#')==NULL){
      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);
  }
/*-------------------------------------------------------------
        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 */
double linint(double xdes,double x1, double x2, double y1, double y2)
{
  double m;
  m=(y2-y1)/(x2-x1);
  return (m*(xdes-x1)+y1);
}

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("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]);
     }

  }
/*------------------------------------------------------------------------
            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 */
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;
}

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("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]);
     }

  }
/*-------------------------------------------------------------------
  integration routines
---------------------------------------------------------------------*/
double integtrap(double (*func)(double),double prev,double low,double high, int step)
{
 double s,npts,stpsze,sum,x;
 int pw2, idx;
 if (step==1){
   return(s=0.5*(high-low)*((*func)(high)+(*func)(low)));
 }
 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);
   }
   s=0.5*(s+(high-low)*sum/npts);
   return s;
 }
}
double integ1(double (*func)(double),double low, double high, double err)
{
  double out,outprev;
  int idx;
  out=integtrap(*func,0.0,low,high,1);
  outprev=out;
  out=integtrap(*func,out,low,high,2);
  /*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);
   /* printf("out %g outprev %g \n",out,outprev);*/
  }
  return out;
}
/*---------------------------------------------------------------------------
   Routine for finding zeros.
  Modified version of rtbis from "Numerical Recipes in C: pg 354
 -----------------------------------------------------------------------------*/

double zero_find(double (*func)(double, double),double yval,double xmin,double xmax, double tol)
{
  double xl,xh,f,fmid,xmid,dx,rtb;
  xl=xmin;
  xh=pow(10,(log10(xmin)+yval*(log10(xmax)-log10(xmin))));
  f=(*func)(xl,yval);
  fmid=(*func)(xh,yval);
  while (fmid*f>=0.0){
    xh=xh+(xh-xl)*2.0;
    fmid=(*func)(xh,yval);
  }
  dx=xh-xl;
  rtb=xl;
  while(fabs((*func)(rtb,yval))>tol){
    dx=dx*0.5;
    xmid=rtb+dx;
    fmid=(*func)(xmid,yval);
    if (fmid<0){
      rtb=xmid;
    }
  }
 return rtb;
}
/*----------------------------------------------------------------------------
Routine for calculating Probability distribution
----------------------------------------------------------------------------*/
void Pcalc(double (*func)(double),double llim, double hlim, double *xvec, double *Prob, int veclen, int *idxstart, int *idxstop)
{
  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);
   idx1++;
  while(xvec[idx1]<=hlim){
     junk=integ1((*func),xvec[idx1-1],xvec[idx1],0.001);
     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;
    }
  }
}
/*----------------------------------------------------------------------------
Routine for calculating t Probability distribution
----------------------------------------------------------------------------*/
void tPcalc(double (*func)(double),double llim, double hlim, double *xvec, double *Prob, int veclen, int *idxstart, int *idxstop)
{
  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);
  while(xvec[idx1]<=hlim){
     junk=integ1((*func),xvec[idx1-1],xvec[idx1],0.001);
     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])*/;
    }
  }
}
/*-----------------------------------------------------------------
Functions for random energy generation
------------------------------------------------------------------*/
double xonly(double x)
{
  return linfuncint(x,xylength,inxvec,inyvec);
}
double Pfunc(double x, double y)
{
  return quadfuncint(x,xylength,inxvec,Pvec)-y;
}
/*----------------------------------------------------------------
Functions for random time generation
------------------------------------------------------------------*/
double txonly(double t)
{
  return linfuncint(t,ntvals,txval,tyval);
}
double tPfunc(double t,double y)
{
  return quadfuncint(t,ntvals,txval,tyval)-y;
}

%}
INITIALIZE
%{
  FILE *fp;
  double llim, hlim,ltlim,htlim,junk;
  double tycol[200];
  double **Imat;
  int idx1,idx2;
  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));
  tyval=malloc(500*sizeof(double));
  txval=malloc(500*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 ",S_filename);
   sns_source_load(S_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);
  Pcalc(xonly,llim,hlim,inxvec,Pvec,xylength,&idxstart,&idxstop);
  /*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);
      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 = xw*yh/(dist*dist);
  for(idx1=0;idx1<200;idx1++){
    free(Imat[idx1]);
  }
  free(Imat);
printf("Finished calculating probability distribution\n");

%}
TRACE
%{
  double theta,phi,v,tauv,E,Eval,tval,randp;
  double hdivmin,hdivmax,vdivmin,vdivmax;
  int idx1,idx3,Eidxl,Eidxh;
  p=p_in;
  z=0;

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

 hdivmax=atan((xw/2.0-x)/dist);
 hdivmin=atan(-(xw/2.0+x)/dist);
 vdivmax=atan((yh/2.0-y)/dist);
 vdivmin=atan(-(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 */
    Eval=zero_find(Pfunc,randp,inxvec[idxstart],inxvec[idxstop],1e-5);
   /* 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){
      tval=zero_find(tPfunc,randp,txval[tidxstart],txval[tidxstop-1],1e-5);}
    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*=INorm2/mcget_ncount();
%}
FINALLY
%{
  int idxf;
  free(txval);free(tyval);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=-width/2.0;y1=-height/2.0;x2=width/2.0;y2=height/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