/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2008, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: ViewModISIS
*
* %I
* Written by: G. Skoro, based on ViewModerator4 from S. Ansell
* Date: February 2016
* Origin: ISIS
* 
* ISIS Moderators
*
* %D
* Produces a neutron distribution at the ISIS TS1 or TS2 beamline shutter front face (or corresponding moderator) position.  
* The Face argument determines which TS1 or TS2 beamline is to be sampled by using corresponding Etable file.
* Neutrons are created having a range of energies determined by the E0 and E1 arguments.
* Trajectories are produced such that they pass through the shutter front face (RECOMENDED) or moderator face (defined by 
* xw and yh) and a focusing rectangle (defined by focus_xw, focus_yh and dist). 
*
* Example 1:   ViewModISISver1(Face="TS1_S04_Merlin.mcstas", E0 = E_min, E1 = E_max,
*                               dist = 1.7, focus_xw = 0.094, focus_yh = 0.094, modPosition=0,
*                               xw=0.12,yh=0.12)
*
* 	MERLIN simulation. 
*	modPosition=0 so initial surface is at (near to) the moderator face.
*	In this example,  focus_xw and focus_yh are chosen to be identical to the shutter opening dimension.
*	dist = 1.7 is the 'real' distance to the shutter front face => This is TimeOffset value (=170 [cm])
*	from TS1_S04_Merlin.mcstas file.
*
*
* Example 2:   ViewModISISver1(Face="Let_timeTest_155.mcstas", E0 = E_min, E1 = E_max,
*              modPosition=1, xw=0.196, yh = 0.12, focus_xw = 0.04, focus_yh = 0.094, dist = 0.5)
*
* 	LET simulation. 
*	modPosition=1 so initial surface is at front face of shutter insert.
*
*	(IMPORTANT) If modPosition=1, the xw and yh values are obsolete. The dimensions of
*	initial surface are automatically calculated using "RDUM" values in Let_timeTest_155.mcstas file.
*
*	In this example,  focus_xw and focus_yh are arbitrary chosen to be identical to the shutter opening dimension.
* 
*
*
* N.B. Absolute normalization: The result of the Mc-Stas simulation will show neutron intensity for beam current of 1 uA.
*
*
* %Parameters
* INPUT PARAMETERS:
*
* Face:        [string]  Etable filename
* E0:          [meV]     Lower edge of energy distribution
* E1:          [meV]     Upper edge of energy distribution
* modPosition: [int]     Defines the initial surface for neutron distribution production. Possible values = 0 (at moderator) or 1 (at shutter insert)
* xw:          [m]       Moderator width
* yh:          [m]       Moderator height
* focus_xw:    [m]       Width of focusing rectangle
* focus_yh:    [m]       Height of focusing rectangle
* dist:        [m]       Distance from source surface to the focusing rectangle
* verbose:     [int]     Flag to output debugging information
* beamcurrent: [uA]      ISIS beam current
* 
* %E
*******************************************************************************/
DEFINE COMPONENT ViewModISIS
SETTING PARAMETERS (string Face="TS1_S04_Merlin.mcstas",E0 , E1, modPosition=0,
  xw=0.12, yh=0.12, focus_xw=0.094, focus_yh=0.094,
  dist=1.7, int verbose=0, beamcurrent=1)

SHARE
%{

#include <ctype.h>

typedef struct
{
  int nEnergy;        ///< Number of energy bins
  int nTime;          ///< number of time bins

  double XAxis;
  double ZAxis;

  double rdumMid;     ///< tally time Window mid point
  double timeOffset;  ///< Time separation
  double* TimeBin;    ///< Time bins
  double* EnergyBin;  ///< Energy bins

  double** Flux;       ///< Flux per bin (integrated)
  double* EInt;        ///< Integrated Energy point
  double Total;        ///< Integrated Total

} Source_ViewMod;


double**
matrix(const int m,const int n)
 /*!
   Determine a double matrix
 */
{
  int i;
  double* pv;
  double** pd;

  if (m<1) return 0;
  if (n<1) return 0;
  pv = (double*) malloc(m*n*sizeof(double));
  pd = (double**) malloc(m*sizeof(double*));
  if (!pd)
    {
      fprintf(stderr,"No room for matrix!\n");
      exit(1);
    }
  for (i=0;i<m;i++)
    pd[i]=pv + (i*n);
  return pd;
} /* matrix */

#pragma acc routine seq
double polInterp(double* X,double* Y,int Psize,double Aim)
  /*!
    returns the interpolated polynomial between Epnts
    and the integration
    \param X :: X coordinates
    \param Y :: Y coordinates
    \param Psize :: number of valid point in array to use
    \param Aim :: Aim point to intepolate result (X coordinate)
    \returns Energy point
  */
{
  double out,errOut;         /* out put variables */
  double *C = malloc(Psize*sizeof(double));
  double *D = malloc(Psize*sizeof(double));
  double testDiff,diff;

  double w,den,ho,hp;           /* intermediate variables */
  int i,m,ns;


  ns=0;
  diff=fabs(Aim-X[0]);
  C[0]=Y[0];
  D[0]=Y[0];
  for(i=1;i<Psize;i++)
    {
      testDiff=fabs(Aim-X[i]);
      if (diff>testDiff)
	{
	  ns=i;
	  diff=testDiff;
	}
      C[i]=Y[i];
      D[i]=Y[i];
    }

  out=Y[ns];
  ns--;              /* Now can be -1 !!!! */

  for(m=1;m<Psize;m++)
    {
      for(i=0;i<Psize-m;i++)
	{
	  ho=X[i]-Aim;
	  hp=X[i+m]-Aim;
	  w=C[i+1]-D[i];
	  /*	  den=ho-hp;  -- test !=0.0 */
	  den=w/(ho-hp);
	  D[i]=hp*den;
	  C[i]=ho*den;
	}

      errOut= (2*(ns+1)<(Psize-m)) ? C[ns+1] : D[ns--];
      out+=errOut;
    }
  free(C);
  free(D);
  return out;
} /* polInterp */

#pragma acc routine seq
int binSearch(int Npts,double* AR,double V)
  /*!
    Object is to find the point in
    array AR, closest to the value V
    Checked for ordered array returns lower of backeting objects
  */
{
  int klo,khi,k;
  if (Npts<=0)
    return 0;
  if (V>AR[Npts-1])
    return Npts;

  if(AR[0]>0.0)AR[0]=0.0;

  if (V<AR[0])
    {
      // if(AR[0]>0.0)AR[0]=0.0;
      return 0;
    }
  klo=0;
  khi= Npts-1;
  while (khi-klo >1)
    {
      k=(khi+klo) >> 1;    // quick division by 2
      if (AR[k]>V)
	khi=k;
      else
	klo=k;
    }
  return khi;
} /* binSearch */

int
cmdnumberD(char *mc,double* num)
 /*!
   \returns 1 on success 0 on failure
 */
{
  int i,j;
  char* ss;
  char **endptr;
  double nmb;
  int len;

  len=strlen(mc);
  j=0;

  for(i=0;i<len && mc[i] &&
	(mc[i]=='\t' || mc[i]==' '  || mc[i]==',');i++);
  if(i==len || !mc[i]) return 0;
  ss=malloc(sizeof(char)*(len+1));

  for(;i<len && mc[i]!='\n' && mc[i]
	&& mc[i]!='\t' && mc[i]!=' ' && mc[i]!=',';i++)
    {
      ss[j]=mc[i];
      j++;
    }
  if (!j)
    {
      free(ss);
      return 0;         //This should be impossible
    }
  ss[j]=0;
  endptr=malloc(sizeof(char*));
  nmb = strtod(ss,endptr);
  if (*endptr != ss+j)
    {
      free(endptr);
      free(ss);
      return 0;
    }
  *num = (double) nmb;
  for(j=0;j<i && mc[j];j++)
    mc[j]=' ';
  free(endptr);
  free(ss);
  return 1;
} /* cmdnumberD */

int
notComment(char* Line)
 /*!
   \returns 0 on a comment, 1 on a non-comment
 */
{
  int len,i;

  len=strlen(Line);
  for(i=0;i<len && isspace(Line[i]);i++);

  if (!Line[i] || Line[i]=='c' || Line[i]=='C' ||
      Line[i]=='!' || Line[i]=='#')
    return 0;
  return 1;
} /* notComment */

void
orderEnergy(double* A,double *B)
{
  double tmp;
  if (*A>*B)
    {
      tmp=*A;
      *A=*B;
      *B=tmp;
    }
  return;
} /* orderEnergy */

int
timeStart(char* Line)
 /*!
   Search for a word time at the start of
   the line.
   \param Line :: Line to search
   \returns 1 on success 0 on failure
 */
{
  int len,i;

  len=strlen(Line);
  for(i=0;i<len && isspace(Line[i]);i++);
  if (len-i<4) return 0;
  return (strncmp(Line+i,"time",4)) ? 0 : 1;
} /* timeStart */

void
cutToNumber(char* Line)
{
  int len,i;

  len=strlen(Line);
  for(i=0;i<len && !isdigit(Line[i]) && Line[i]!='-';i++)
    Line[i]=' ';
  return;
} /* cutToNumber */

int
findTimeOffset(char* Line,double* TO)
  /*!
    Determine the time offset from the file
    \return 1 on success
   */
{
  int len,i;
  double D;

  len=strlen(Line);
  for(i=0;i<len && Line[i]!='T';i++);
  if (len-i<11) return 0;

  if (strncmp(Line+i,"TimeOffset",11) &&
      cmdnumberD(Line+i+11,&D))
    {
      *TO=D/100.0;
      return 1;
    }
  return 0;
} /* findTimeOffset */

int
findRDUM(char* Line)
 /*!
   Search for a word rdum at the start of
   the line.
   \param Line :: Line to search
   \returns 1 on success 0 on failure
 */
{
  int len,i;

  len=strlen(Line);
  for(i=0;i<len && Line[i]!='r';i++);
  if (len-i<4) return 0;
  return (strncmp(Line+i,"rdum",4)) ? 0 : 1;
} /* findRDUM */

double 
convertEnergy(double E)
  /*!
    Convert the energy from eV [not change] or 
    from negative -ve which is angstrom
  */
{
  return (E>0.0) ? E : 81.793936/(E*E);
} /* convertEnergy */

double 
EtoLambda(double E)
  /*!
    Convert the energy from eV [not change] o
    to lambda [A]
  */
{
  return sqrt(81.793936/E);
} /* EtoLambda */
  

int
timeEnd(char* Line)
 /*!
   Search for a word time at the start of
   the line.
   \param Line :: Line to search
   \returns 1 on success 0 on failure
 */
{
  int len,i;

  len=strlen(Line);
  for(i=0;i<len && isspace(Line[i]);i++);
  if (len-i<5) return 0;
  return (strncmp(Line+i,"total",5)) ? 0 : 1;
} /* timeEnd */

int
energyBin(char* Line,double Einit,double Eend,double* Ea,double* Eb)
     /*!
       Search for a word "energy bin:" at the start of
       the line. Then separte off the energy bin values
       \param Line :: Line to search
       \param Ea :: first energy bin [meV]
       \param Eb :: second energy bin [meV]
       \returns 1 on success 0 on failure
     */
{
  int len,i;
  double A,B;

  len=strlen(Line);
  for(i=0;i<len && isspace(Line[i]);i++);
  if (len-i<11) return 0;


  if (strncmp(Line+i,"energy bin:",11))
    return 0;

  i+=11;
  if (!cmdnumberD(Line+i,&A))
    return 0;
  // remove 'to'
  for(;i<len-1 && Line[i]!='o';i++);
  i++;
  if (!cmdnumberD(Line+i,&B))
    return 0;
  A*=1e9;
  B*=1e9;
  *Ea=A;
  *Eb=B;
  if (*Eb>Einit && *Ea<Eend)
    return 1;
  return 0;
} /* energyBin */

double
calcFraction(double EI,double EE,double Ea,double Eb)
 /*!
   Calculate the fraction of the bin between Ea -> Eb
   that is encompassed by EI->EE
 */
{
  double frac;
  double dRange;

  if (EI>Eb)
    return 0.0;
  if (EE<Ea)
    return 0.0;

  dRange=Eb-Ea;
  frac=(EI>Ea) ? (Eb-EI)/dRange : 1.0;


  frac-=(EE<Eb) ? (Eb-EE)/dRange : 0.0;

  //  if(frac != 1.0)
  //  fprintf(stderr,"frac %g, Ea %g,Eb %g, EI %g, EE %g\n",frac,Ea,Eb,EI,EE);

  return frac;
} /* calcFraction */

double
calcRDum(double* RPts, Source_ViewMod TS)
  /*!
    Caluclate the mean distance from the window centre point
   */
{
  double sum,zMid;
  int i;

  sum=0.0;
  for(i=0;i<4;i++)
    {
      fprintf(stderr,"RDUM %g %g %g\n",RPts[i*3],RPts[i*3+1],RPts[i*3+2]);
      sum+=sqrt(RPts[i*3]*RPts[i*3]+
		RPts[i*3+1]*RPts[i*3+1]);
    }
  sum/=400.0;   // Convert to metres from cm
  return sum;
} /* calcRDum */

void
processHeader(FILE* TFile, Source_ViewMod *TS, double modPosition, double xw, double yh)
  /*!
    Determine the window and the time offset
    \param TFile :: file [rewound on exit]
  */
{
  char ss[255];          /* BIG space for line */
  int rdumCnt;
  double Ea,Eb;
  double RPts[12];
  int timeOffsetFlag;
  int i,j;
  double D;
  
  rdumCnt=0;
  timeOffsetFlag=0;

  while(fgets(ss,255,TFile))
    {
      if (!timeOffsetFlag) {
	double toff;
	timeOffsetFlag=findTimeOffset(ss,&toff);
	TS->timeOffset=toff;
      }
      if (!rdumCnt)
	rdumCnt=findRDUM(ss);
      if (rdumCnt && rdumCnt<5)
	{
	  cutToNumber(ss);
	  for(i=0;i<3 && cmdnumberD(ss,&D);i++)
	    RPts[(rdumCnt-1)*3+i]=D;
	  rdumCnt++;
	}
      // EXIT CONDITION:
      if (rdumCnt*timeOffsetFlag==5)
	{
	  for(j=0;j<3;j++)
	    {
	      TS->ZAxis+=(RPts[3+j]-RPts[j])*(RPts[3+j]-RPts[j]);
	      TS->XAxis+=(RPts[6+j]-RPts[3+j])*(RPts[6+j]-RPts[3+j]);
	    }	


	  TS->ZAxis=sqrt(TS->ZAxis)/100.0;   // convert to metres from cm
	  TS->XAxis=sqrt(TS->XAxis)/100.0;
              if (!modPosition)
                 {
	          TS->ZAxis=yh;
	          TS->XAxis=xw;
                  }
	  fprintf(stderr,"Time off sec == %g \n",TS->timeOffset);
	  fprintf(stderr,"mod size z == %g \n",TS->ZAxis);
	  // TS.rdumMid=calcRDum(RPts, TS);
	  TS->rdumMid=TS->timeOffset; // Goran
	  return;
	}
    }
  return;
} /* processHeader */



int
readHtable(FILE* TFile,const double Einit,const double Eend, 
  Source_ViewMod *TS, double modPosition, double xw, double yh, int verbose)
/*!
  Process a general h.o file to create an integrated
  table of results from Einit -> Eend
  \param Einit :: inital Energy
  \param Eend  :: final energy
*/
{
  char ss[255];          /* BIG space for line */
  double Ea,Eb;
  double T,D;
  double Efrac;          // Fraction of an Energy Bin
  int Ftime;             // time Flag
  int eIndex;             // energy Index
  int tIndex;             // time Index
  double Tsum;           // Running integration
  double Efraction;      // Amount to use for an energy/time bin

  // extern Source TS;

  int DebugCnt;
  int i;
  /*!
    Status Flag::
    Ftime=1 :: [time ] Reading Time : Data : Err [Exit on Total]


    /*
    Double Read File to determine how many bins and
    memery size
  */
  if (!TFile) return(0);
  Ea=0.0;
  Eb=0.0;
  fprintf(stderr,"Energy == %g %g\n",Einit,Eend);
  eIndex= -1;
  DebugCnt=0;
  Ftime=0;
  tIndex=0;
  TS->nTime=0;
  TS->nEnergy=0;
  processHeader(TFile, TS, modPosition, xw, yh);

  // Read file and get time bins:
  while(fgets(ss,255,TFile) && Eend>Ea)
    {
      if (notComment(ss))
	{
	  DebugCnt++;
          if (!Ftime)
	    {
	      if (energyBin(ss,Einit,Eend,&Ea,&Eb))
		{
		  if (eIndex==0)
		    TS->nTime=tIndex;
		  eIndex++;
		}
	      else if (timeStart(ss))
		{
		  Ftime=1;
		  tIndex=0;
		}
	    }
	  else  // In the time section
	    {
	      if (timeEnd(ss))     // Found "total"
		Ftime=0;
	      else
		{
		  // Need to read the line in the case of first run
		  if (TS->nTime==0)
		    {
		      if (cmdnumberD(ss,&T) &&
			  cmdnumberD(ss,&D))
			tIndex++;
		    }
		}
	    }
	}
    }
  // Plus 2 since we have a 0 counter and we have missed the last line.
  TS->nEnergy=eIndex+2;
  if (!TS->nTime && tIndex)
    TS->nTime=tIndex;
  // printf("tIndex %d %d %d %d \n",tIndex,eIndex,TS->nEnergy,TS->nTime);

  /* SECOND TIME THROUGH:: */
  rewind(TFile);

  TS->Flux=matrix(TS->nEnergy,TS->nTime);
  TS->EInt=(double*) malloc(TS->nEnergy*sizeof(double));
  TS->TimeBin=(double*) malloc(TS->nTime*sizeof(double));
  TS->EnergyBin=(double*) malloc(TS->nEnergy*sizeof(double));

  Tsum=0.0;
  Ea=0.0;
  Eb=0.0;
  eIndex=-1;
  DebugCnt=0;
  Ftime=0;
  tIndex=0;
  TS->EInt[0]=0.0;
  // Read file and get time bins
  while(fgets(ss,255,TFile) && Eend>Ea)
    {
      if (notComment(ss))
	{
	  DebugCnt++;
          if (!Ftime)
	    {
	      if (energyBin(ss,Einit,Eend,&Ea,&Eb))
		{
		  eIndex++;
		  TS->EnergyBin[eIndex]=(Einit>Ea) ? Einit : Ea;
		  Efraction=calcFraction(Einit,Eend,Ea,Eb);
		  Ftime++;
		}
	    }
	  else if (Ftime==1)
	    {
	      if (timeStart(ss))
		{
		  Ftime=2;
		  tIndex=0;
		}
	    }

	  else           // In the time section
	    {
	      if (timeEnd(ss))     // Found "total"
		{
		  Ftime=0;
		  TS->EInt[eIndex+1]=Tsum;
		}
	      else
		{
		  // Need to read the line in the case of first run
		  if (cmdnumberD(ss,&T) &&
		      cmdnumberD(ss,&D))
		    {
		      TS->TimeBin[tIndex]=T/1e8;     // convert Time into second (from shakes)
		      Tsum+=D*Efraction;
		      TS->Flux[eIndex][tIndex]=Tsum;
		      tIndex++;
		    }
		}
	    }
	}
    }

  TS->EnergyBin[eIndex+1]=Eend;
  TS->Total=Tsum;

  //  printf("tIndex %d %d %d \n",tIndex,eIndex,TS.nTime);
  //printf("Tsum %g \n",Tsum);
  //fprintf(stderr,"ebin1 ebinN %g %g\n",TS.EnergyBin[0],TS.EnergyBin[TS.nEnergy-1]);

  return 1;
} // readHtable


#pragma acc routine seq
void getPoint(double* TV,double* EV,double* lim1, double* lim2, Source_ViewMod TS, _class_particle *_particle)
 /*!
   Calculate the Time and Energy
   by sampling the file.
   Uses TS table to find the point
   \param TV ::
   \param EV ::
   \param lim1 ::
   \param lim2 ::
 */
{
  int i;

  double R0,R1,R,Rend;
  int Epnt;       ///< Points to the next higher index of the neutron integral
  int Tpnt;
  int iStart,iEnd;
  double TRange,Tspread;
  double Espread,Estart;
  double *EX;

  // So that lowPoly+highPoly==maxPoly
  const int maxPoly=6;
  const int highPoly=maxPoly/2;
  const int lowPoly=maxPoly-highPoly;

  // static int testVar=0;

  R0=rand01();
  /* if (testVar==0)
    {
    R0=1.0e-8;
    testVar=1;
    }
  */
  Rend=R=TS.Total*R0;
  // This gives Eint[Epnt-1] > R > Eint[Epnt]
  Epnt=binSearch(TS.nEnergy-1,TS.EInt,R);

  //      if (Epnt < 0)
  //   Epnt=1;
  Tpnt=binSearch(TS.nTime-1,TS.Flux[Epnt-1],R);
  //  fprintf(stderr,"TBoundaryX == %12.6e %12.6e \n",TS.TimeBin[Tpnt-1],TS.TimeBin[Tpnt]);
  //  fprintf(stderr,"TFlux == %12.6e %12.6e %12.6e \n\n",TS.Flux[Epnt-1][Tpnt-1],R,TS.Flux[Epnt-1][Tpnt]);
  //  if (Epnt == -1)
  //{
  //    Epnt=0;
  // fprintf(stderr,"\n Rvals == %g %d %d %g\n",R,Epnt,Tpnt,TS.TimeBin[0]);
  //  fprintf(stderr,"EInt == %d %12.6e %12.6e %12.6e %12.6e \n",Epnt,TS.EInt[Epnt-1],R,TS.EInt[Epnt],TS.EInt[Epnt+1]);
  // printf("EBoundary == %12.6e %12.6e \n",TS.EnergyBin[Epnt],TS.EnergyBin[Epnt+1]);

  //  fprintf(stderr,"TFlux == %12.6e %12.6e %12.6e \n\n",TS.Flux[Epnt+1][Tpnt],R,TS.Flux[Epnt+1][Tpnt+1]);
  // }

  if(R < TS.Flux[Epnt-1][Tpnt-1] || R >TS.Flux[Epnt-1][Tpnt] )
    {
#ifndef OPENACC
      fprintf(stderr, "outside bin limits Tpnt/Epnt problem  %12.6e %12.6e %12.6e \n",TS.Flux[Epnt-1][Tpnt-1],R,TS.Flux[Epnt-1][Tpnt]);
#endif
    }

  if(Epnt == 0)
    {
      Estart=0.0;
      Espread=TS.EInt[0];
      *EV=TS.EnergyBin[1];
    }
  else
    {
      Estart=TS.EInt[Epnt-1];
      Espread=TS.EInt[Epnt]-TS.EInt[Epnt-1];
      *EV=TS.EnergyBin[Epnt+1];
    }

  if (Tpnt==0 || Epnt==0)
    {
#ifndef OPENACC
      fprintf(stderr,"BIG ERROR WITH Tpnt: %d and Epnt: %d\n",Tpnt,Epnt);
      exit(1);
#endif
    }
  if (Tpnt==TS.nTime)
    {
#ifndef OPENACC
      fprintf(stderr,"BIG ERROR WITH Tpnt and Epnt\n");
      exit(1);
#endif

      *TV=0.0;
      Tspread=TS.Flux[Epnt-1][0]-TS.EInt[Epnt-1];
      TRange=TS.TimeBin[0];
      R-=TS.EInt[Epnt-1];
    }
  else
    {
      *TV=TS.TimeBin[Tpnt-1];
      TRange=TS.TimeBin[Tpnt]-TS.TimeBin[Tpnt-1];
      Tspread=TS.Flux[Epnt-1][Tpnt]-TS.Flux[Epnt-1][Tpnt-1];
      R-=TS.Flux[Epnt-1][Tpnt-1];
    }
  //  printf("R == %12.6e\n",R);
  R/=Tspread;
  *TV+=TRange*R;


  R1=TS.EInt[Epnt-1]+Espread*rand01();
  iStart=Epnt>lowPoly ? Epnt-lowPoly : 0;                  // max(Epnt-halfPoly,0)
  iEnd=TS.nEnergy>Epnt+highPoly ? Epnt+highPoly : TS.nEnergy-1;  // min(nEnergy-1,Epnt+highPoly

  *EV=polInterp(TS.EInt+iStart,TS.EnergyBin+iStart,1+iEnd-iStart,R1);

  //  fprintf(stderr,"Energy == %d %d %12.6e %12.6e \n",iStart,iEnd,R1,*EV);
  //  fprintf(stderr,"bins == %12.6e %12.6e %12.6e %12.6e \n",TS.EnergyBin[iStart],TS.EnergyBin[iEnd],
  //	  TS.EInt[Epnt],TS.EInt[Epnt-1]);

    if(*TV < TS.TimeBin[Tpnt-1] || *TV > TS.TimeBin[Tpnt])
    {
#ifndef OPENACC
      fprintf(stderr,"%d Tpnt %d Tval %g Epnt %d \n",TS.nTime,Tpnt,*TV,Epnt);
      fprintf(stderr,"TBoundary == %12.6e,%g , %12.6e \n\n",TS.TimeBin[Tpnt-1],*TV,TS.TimeBin[Tpnt]);
#endif
    }


  if(*EV < *lim1 || *EV > *lim2)
    {
#ifndef OPENACC
      fprintf(stderr,"outside boundaries\n Epnt= %d, Tpnt= %d binlo %g|%g| binhi %g \n",Epnt,Tpnt,TS.EnergyBin[Epnt-1],*EV,TS.EnergyBin[Epnt]);

      fprintf(stderr,"TS == %g %g :: %d %d \n",TS.EInt[Epnt-1],TS.EInt[Epnt],iStart,iEnd);
      fprintf(stderr,"Points (%g) == ",R1);

      for(i=0;i<iEnd-iStart;i++)
	fprintf(stderr," %g %g",*(TS.EInt+i+iStart),*(TS.EnergyBin+iStart+i));
      fprintf(stderr,"\n");
#endif
    }
  return;
} /* getPoint */


int
cmdnumberI(char *mc,int* num,const int len)
  /*!
    \param mc == character string to use
    \param num :: Place to put output
    \param len == length of the character string to process
    returns 1 on success and 0 on failure
  */
{
  int i,j;
  char* ss;
  char **endptr;
  double nmb;

  if (len<1)
    return 0;
  j=0;
  
  for(i=0;i<len && mc[i] &&
	(mc[i]=='\t' || mc[i]==' '  || mc[i]==',');i++);
  if(i==len || !mc[i]) return 0;
  ss=malloc(sizeof(char)*(len+1));
  /*  char *ss=new char[len+1]; */
  for(;i<len && mc[i]!='\n' && mc[i]
	&& mc[i]!='\t' && mc[i]!=' ' && mc[i]!=',';i++)
    {
      ss[j]=mc[i];
      j++;
    }
  if (!j)
    {
      free(ss);
      return 0;         //This should be impossible
    }
  ss[j]=0;
  endptr=malloc(sizeof(char*));
  nmb = strtod(ss,endptr);
  if (*endptr != ss+j)
    {
      free(endptr);
      free(ss);
      return 0;
    }
  *num = (double) nmb;
  for(j=0;j<i && mc[j];j++)
    mc[j]=' ';
  free(endptr);
  free(ss);
  return 1;
} /* cmdnumberI */


FILE* 
openFile(char* FileName)
/*
  Tries to open the file:
  (i) In current working directory 
  (ii) In MC_Path directory + ISIS_tables extension
*/
{
  FILE* efile=0;
  int fLen;
  char extFileName[1024];
  char testFileName[2048];
  char mct[2048];
  
  fLen=strlen(FileName);

  if (fLen>512) 
    {
      fprintf(stderr,"Filename excessivley long [EXIT]:\n %s\n",FileName);
      exit(1);
    }
  
  
  strcpy(extFileName,FileName);
  strcpy(extFileName+fLen,".mcstas");
  
  /* Now check for the requested moderator file. In terms of precedence, 
     1) Use MCTABLES location if available
     2) Is a local file available in PWD?
     3) Is there an ISIS_tables in PWD?
     4) Is the file available from the MCSTAS/contrib/ISIS_tables folder?
     - otherwise fail */

  fprintf(stderr,"Searching for %s in multiple locations... -- \n",FileName);

  /* 1) Is MCTABLES set and file located there? */
  if (getenv("MCTABLES"))
    {
      strcpy(mct, getenv("MCTABLES"));
      sprintf(testFileName, "%s%c%s", mct, MC_PATHSEP_C, FileName);
      efile=fopen(testFileName,"r");
      if (efile) {
	fprintf(stderr," - Found in MCTABLES folder %s!\n",mct);
	return efile;
      }
    }

  /* 2) Is the file located in working dir? */
  efile=fopen(FileName,"r");
  if (efile) {
    fprintf(stderr," - Found in current directory!\n");
    return efile;
  }
  
  efile=fopen(extFileName,"r");
  if (efile) return efile; 
  
  /* 3) Is the file in a local 'tablefiles' folder? */
  sprintf(testFileName,"%s%c%s","ISIS_tables",MC_PATHSEP_C,FileName);
  efile=fopen(testFileName,"r");
  if (efile) {
    fprintf(stderr," - Found in local ISIS_tables directory!\n");
    return efile;
  }

  /* 4) Is the file available within the MCSTAS install dir? */
  sprintf(testFileName,"%s%c%s%c%s%c%s",
    MCSTAS,MC_PATHSEP_C,"data",MC_PATHSEP_C,"ISIS_tables",MC_PATHSEP_C,FileName);
  efile=fopen(testFileName,"r");
  if (efile) {
    fprintf(stderr," - Found in MCSTAS system dir: \n %s%c%s%c%s\n",
      MCSTAS,MC_PATHSEP_C,"contrib",MC_PATHSEP_C,"ISIS_tables");
    return efile;
  }

  /* If we reach here, the file was not found - raise fatal error */
  fprintf(stderr,"%s - Not found! \nPlease check your McStas installation and/or MCTABLES environment variable!\n",FileName);
  exit(1);
} /* openFile */

double strArea(Source_ViewMod TS, double focus_xw, double focus_yh, double dist)
{
  /*
    Returns the mean Str view of the viewport
    This integrates over each point on the window focus_xw to focus_yh
    View port is symmetric so use only 1/4 of the view
    for the calcuation.
    Control Values TS.XAxis TS.ZAxis focus_xw focus_yh
  */
  
  double A;
  double Vx,Vy;        // view temp points
  double Mx,My;        // moderator x,y
  double D2;           // Distance ^2
  double projArea;     // viewport projection to moderator
  int i,j,aa,bb;       // loop variables
  int NStep;

  NStep=50;
  D2=dist*dist;
  A=0.0;
  fprintf(stderr,"TS axis == %g %g\n",TS.XAxis,TS.ZAxis);
  fprintf(stderr,"AW axis == %g %g\n",focus_xw,focus_yh);
  for(i=0;i<NStep;i++)              // Mod X
    {
      Mx=(i+0.5)*TS.XAxis/(NStep*2.0);
      for(j=0;j<NStep;j++)         // Mod Y
	{
	  My=(j+0.5)*TS.ZAxis/(NStep*2.0);
	  // Position on moderator == (Mx,My)
	  for(aa= -NStep;aa<=NStep;aa++)  //view port
	    for(bb= -NStep;bb<=NStep;bb++)
	      {
		Vx=aa*focus_xw/(2.0*NStep+1.0);
		Vy=bb*focus_yh/(2.0*NStep+1.0);
		A+=1.0/((Mx-Vx)*(Mx-Vx)+(My-Vy)*(My-Vy)+D2);
	      }
	}
    }
  //change to Mx*My 
  A*= (TS.XAxis*TS.ZAxis)/(NStep*NStep*(2.0*NStep+1.0)*(2.0*NStep+1.0));
  // Correct for the area of the viewport. (tables are per cm^2)
  A*=focus_xw*focus_yh*10000;

  //   testing only - Goran              
    // if (!modPosition)
        // {
            // projArea=focus_xw*focus_yh/tally/tally*(tally+dist)*(tally+dist);
	    // A*=TS.XAxis*TS.ZAxis/(TS.XAxis*TS.ZAxis-projArea);
         // }  

  fprintf(stderr,"Viewport == %g %g Moderator size == (%g * %g) m^2 \n",focus_xw,focus_yh,TS.XAxis,TS.ZAxis);
  fprintf(stderr,"Dist == %g (metres) \n",dist);
  fprintf(stderr,"Viewport Solid angle == %g str\n",A/(focus_xw*focus_yh*10000.0));
  fprintf(stderr,"Solid angle used == %g str\n",A);
  return A;
} /* strArea */

%}

DECLARE %{
  /* global variables */

  double p_in;            /* Polorization term (from McSTAS) */
  int Tnpts;              /* Number of points in parameteriation */

  double scaleSize;        /* correction for the actual area of the moderator viewed */

  double angleArea;       /* Area seen by the window  */

  double Nsim;		/* Total number of neutrons to be simulated */

  int Ncount;          /* Number of neutron simulate so far*/
  Source_ViewMod TS;

  /* runtime variables*/
  double fullAngle;
  double rtE0;       /* runtime Energy minima and maxima so we can use angstroms as negative input */
  double rtE1;
%}

INITIALIZE
%{
  /* READ IN THE ENERGY FILE */
  FILE* Tfile;
  

  Nsim=mcget_ncount();  // Number of points requested.

  Tfile=openFile(Face);	              // Get open file
  rtE0=convertEnergy(E0);
  rtE1=convertEnergy(E1);
  orderEnergy(&rtE0,&rtE1);

  readHtable(Tfile,rtE0,rtE1, &TS, modPosition, xw, yh, verbose);
  fclose(Tfile);
  // Below pragma was needed with PGI 19.x, compilation fails with NVC 20.7 
  //#pragma acc declare create( TS )
  /**********************************************************************/

  Tnpts=0;
  Ncount=0;

  fprintf(stderr,"Face == %s \n",Face);
  fprintf(stderr,"Number of Energy Points == %d\n",TS.nEnergy);
  if (dist<0)
    {
      dist=TS.rdumMid;
      fprintf(stderr,"Setting distance to moderator surface == %g\n",
	      dist);
    }
  /* Do solid angle correction */
  angleArea= strArea(TS, focus_xw, focus_yh, dist);
  // Below pragma was needed with PGI 19.x, compilation fails with NVC 20.7 
  //#pragma acc update host( TS )
  fprintf(stderr,"Totals:: %g %d %d \n",TS.Total,TS.nEnergy,TS.nTime);

%}

TRACE
%{
  double v,r,E;
  double xf,yf,dx,dy;    /* mxp ->max var in param space */
  double Ival,Tval,Eval;

  Ncount++;

  x = TS.XAxis*(0.5-rand01());            /* Get point on shutter enterance */
  y = TS.ZAxis*(0.5-rand01());            /* Get point on shutter enterance */
    
  
  xf = focus_xw*(0.5-rand01());          /* Choose focusing position uniformly */
  yf = focus_yh*(0.5-rand01());          /* Choose focusing position uniformly */

   getPoint(&Tval,&Eval,&rtE0,&rtE1, TS, _particle);
  
  Ival=TS.Total*6.2415093e+12; // Number of proton in 1uAmp

  dx = xf-x;
  dy = yf-y;
  r = sqrt(dx*dx+dy*dy+dist*dist);      // Actual distance to point 
  v = SE2V*sqrt(Eval);                  // Calculate the velocity 
  vz = v*fabs(dist)/r;
  vy = v*dy/r;
  vx = v*dx/r;


   t=Tval-(TS.rdumMid-TS.timeOffset)/vz;
   
     if (modPosition)
     {
     t+=TS.rdumMid/vz;
     }

  p=beamcurrent*angleArea*Ival/Nsim;


#ifndef OPENACC
  if (verbose && !(Ncount % 100000))
    {
      fprintf(stderr,"FPos[%d]=> %g %g %g  \n",Ncount,x,y,z);
      fprintf(stderr,"FDir[%d]=> %g %g %g  \n",Ncount,vx,vy,vz);
      fprintf(stderr,"PlaneAxis %g %g \n",TS.XAxis,fullAngle);
      fprintf(stderr,"RMID %g \n",TS.rdumMid);
      fprintf(stderr,"TimeMid[%d]=> %g\n",Ncount,TS.rdumMid);
      fprintf(stderr,"TimeOffset[%d]=> %g\n",Ncount,TS.timeOffset);
      fprintf(stderr,"TimeTval[%d]=> %g\n",Ncount,Tval);
      fprintf(stderr,"TimeZero[%d]=> %g\n",Ncount,t);
    }
#endif  

%}

MCDISPLAY
%{
  double cirp=0.0,cirq=0.3,pi=3.141592654;
  int pp=0; /* circle drawing parameter*/



  magnify("xy");
  multiline(5,-0.5*TS.XAxis,-0.5*TS.ZAxis,0.0,
	    0.5*TS.XAxis,-0.5*TS.ZAxis,0.0,
	    0.5*TS.XAxis,0.5*TS.ZAxis,0.0,
	    -0.5*TS.XAxis,0.5*TS.ZAxis,0.0,
	    -0.5*TS.XAxis,-0.5*TS.ZAxis,0.0);
  /* circle("xy",0.0,0.0,0.0,cos(cirp)); */

  /*line(0.5*sin(cirp),0.0,0.5*cos(cirp),0.5*sin(cirq),0.0,0.5*cos(cirq));*/

  /*line(-0.5,0.0,0.0,0.0,0.0,0.5);*/

  for (pp=0;pp<=20;pp=pp+2)
    {
      cirp= (pp*(pi/21.0))-(0.5*pi);
      cirq= ((pp+1)*(pi/21.0))-(0.5*pi);
      line(0.5*sin(cirp),0.0,0.5*cos(cirp),0.5*sin(cirq),0.0,0.5*cos(cirq));
    }

  %}

END