/******************************************************************************* * * 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 SNS website . * IMPORTANT: The output units of this component are N/pulse * IMPORTANT: The component needs a FULL PATH to the source input file * 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 DEFINITION PARAMETERS () SETTING PARAMETERS (string filename, xwidth=0.1, yheight=0.12, dist, focus_xw, focus_yh, Emin, Emax) OUTPUT PARAMETERS () 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((idx20)){ idx2++; } if(idx23)){ 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(idx1xylen){ 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]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;idxerr*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++=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++=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++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(idx10){ 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;idx10.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