/******************************************************************************* * * McStas, neutron ray-tracing package * Copyright (C) 1997-2006, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: ISIS_moderator * * * %I * Written by: S. Ansell and D. Champion * Date: AUGUST 2005 * Version: $Revision: 1.4 $ * Origin: ISIS * Release: McStas 1.9 * * ISIS Moderators * * %D * Produces a TS1 or TS2 ISIS moderator distribution. The Face argument determines which moderator * is to be sampled. Each face uses a different Etable file (the location of which is determined by * the MCTABLES environment variable). Neutrons are created having a range of energies * determined by the E0 and E1 arguments. Trajectories are produced such that they pass * through the moderator face (defined by modXsixe and modYsize) and a focusing rectangle * (defined by xh,yh and dist). Please download the documentation for precise instructions for use. * * Example: ISIS_moderator(Face ="water", E0 = 49.0,E1 = 51.0, dist = 1.0, xw = 0.01, * yh = 0.01, modXsize = 0.074, modYsize = 0.074, CAngle = 0.0,SAC= 1) * * %P * INPUT PARAMETERS: * * Face: (word) Name of the face - TS2:groove,hydrogen,narrow,broad - * - TS1:Water,ch4,h2,merlin or instrument name eg maps, crisp etc. * - TS2: W1-9 E1-9 * E0: (meV) Lower edge of energy distribution * E1: (meV) Upper edge of energy distribution * dist: (m) Distance from source to the focusing rectangle * xw: (m) Width of focusing rectangle * yh: (m) Height of focusing rectangle * modXsize: (m) Moderator vertical size * modYsize: (m) Moderator Horizontal size * CAngle: (radians) Angle from the centre line * SAC: (boolean) Apply solid angle correction or not (1/0) * * %L * Further information should be * downloaded from the ISIS MC website. * %E * * * * * *******************************************************************************/ DEFINE COMPONENT ISIS_moderator DEFINITION PARAMETERS (Face) SETTING PARAMETERS (E0, E1,dist,xw,yh,modXsize,modYsize,CAngle,SAC) STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p) POLARISATION PARAMETERS (sx,sy,sz) DECLARE %{ #include typedef struct { int nEnergy; ///< Number of energy bins int nTime; ///< number of time bins double* TimeBin; ///< Time bins double* EnergyBin; ///< Energy bins double** Flux; ///< Flux per bin (integrated) double* EInt; ///< Integrated Energy point double Total; ///< Integrated Total } Source; /* New functions */ int cmdnumberD(char *,double*); int cmdnumberI(char *,int*,const int); double polInterp(double*,double*,int,double); FILE* openFile(char*); FILE* openFileTest(char*); int readHtable(FILE*,const double,const double); int timeStart(char*); int timeEnd(char*); int energyBin(char*,double,double,double*,double*); int notComment(char*); double strArea(); /* 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 TS; /* runtime variables*/ double rtE0,rtE1; /* runtime Energy minima and maxima so we can use angstroms as negative input */ double rtmodX,rtmodY; /* runtime moderator sizes, so that a negative argument may give a default size */ 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;itestDiff) { ns=i; diff=testDiff; } C[i]=Y[i]; D[i]=Y[i]; } out=Y[ns]; ns--; /* Now can be -1 !!!! */ for(m=1;mAR[Npts-1]) return Npts; if(AR[0]>0.0)AR[0]=0.0; if (V0.0)AR[0]=0.0; fprintf(stderr,"here"); 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; } 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;iEinit && *Ea Eb that is encompassed by EI->EE */ { double frac; double dRange; if (EI>Eb) return 0.0; if (EEEa) ? (Eb-EI)/dRange : 1.0; frac-=(EE Eend \param Einit :: inital Energy \parma 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 */ 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; // 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; } void getPoint(double* TV,double* EV,double* lim1, double* lim2) /*! 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; extern Source TS; 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] ) { 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]); } 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) { fprintf(stderr,"BIG ERROR WITH Tpnt: %d and Epnt: %d\n",Tpnt,Epnt); exit(1); } if (Tpnt==TS.nTime) { fprintf(stderr,"BIG ERROR WITH Tpnt and Epnt\n"); exit(1); *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; // printf("R == %12.6e\n",R); *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]) { 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]); } if(*EV < *lim1 || *EV > *lim2) { 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;i0) | (rtE0>0 && E1<0)) { fprintf(stderr,"Cannot have differing signs for E0 and E1, choose Angstroms or meV!\n"); exit(1); } if (rtE0<0 && E1<0) { fprintf (stderr,"converting Angstroms to meV\n"); rtE0=81.793936/(rtE0*rtE0); rtE1=81.793936/(rtE1*rtE1); } if (rtE0>rtE1) { tmp=rtE1; rtE1=rtE0; rtE0=tmp; fprintf (stderr,"%g A -> %g A => %g meV -> %g meV\n",-E0,-E1,rtE0,rtE1); } /**********************************************************************/ Tnpts=0; Ncount=0; fprintf(stderr,"Face == %s \n",Face); for(i=0;Face[i] && Face[i]!=' ';i++) lowerFace[i]=tolower(Face[i]); lowerFace[i]=0; for(i=0;i \n"); for(i=0;i0.0) ? strArea() : 2*3.141592654; else angleArea=1.0; %} TRACE %{ double v,r,E; double xf,yf,dx,dy,w_focus; /* mxp ->max var in param space */ double Ival,Tval,Eval; double Ddist; /* Temp versions of dist */ Ncount++; p=p_in; p=1.0; /* forcing */ z=0; x = 0.5*rtmodX*randpm1(); /* Get point +/-0.5 * */ y = 0.5*rtmodY*randpm1(); xf = 0.5*xw*randpm1(); /* Choose focusing position uniformly */ yf = 0.5*yh*randpm1(); dx = xf-x; dy = yf-y; if (dist>0.0) { r = sqrt(dx*dx+dy*dy+dist*dist); /* Actual distance to point */ Ddist=dist; w_focus = (SAC) ? angleArea : scaleSize*(dist*dist)/(r*r); } else /* Assume that we have a window 1metre infront of the moderator */ /* with size area of detector and solid angle 1.0 */ { r=1.0; w_focus=scaleSize; Ddist=1.0; } getPoint(&Tval,&Eval,&rtE0,&rtE1); //fprintf(stderr,"outside %g mev\n", TS.Total ); if(Eval>rtE1 || Eval %g %g %g %g \n",Ncount,Eval,Tval,TS.Total,Ival); t=Tval; p=w_focus*Ival/Nsim; %} MCDISPLAY %{ double cirp=0.0,cirq=0.3,pi=3.141592654; int pp=0; /* circle drawing parameter*/ magnify("xy"); multiline(5,-0.5*rtmodX,-0.5*rtmodY,0.0, 0.5*rtmodX,-0.5*rtmodY,0.0, 0.5*rtmodX,0.5*rtmodY,0.0, -0.5*rtmodX,0.5*rtmodY,0.0, -0.5*rtmodX,-0.5*rtmodY,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