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