/******************************************************************************* * * 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 * Version: $Revision: 1.1 $ * 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 SNS website . * 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 * 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 * %E *******************************************************************************/ DEFINE COMPONENT SNS_source DEFINITION PARAMETERS (S_filename) SETTING PARAMETERS (width, height, dist, xw, yh, Emin, Emax) OUTPUT PARAMETERS (hdiv,vdiv,p_in) 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((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("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]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;idxerr*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(idx10){ 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(idx10){ 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;idx10.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