/*******************************************************************************
*
* McStas, the neutron ray-tracing package
* Maintained by Kristian Nielsen and Kim Lefmann,
* Copyright 1997-2000 Risoe National Laboratory, Roskilde, Denmark
*
* - Adjusted time intervals to be considered for sampling time August 18, 2010
* - reworked and fixed sample_t=0 option September 24, 2010
* - added analog sampling of energy for sample_E=0 October 15, 2010
* the former sample_E=0 moved to sample_E=1 (uniform in energy)
* the former sample_E=1 moved to sample_E=2 (uniform in wavelength)
* the former sample_E=2 moved to sample_E=3 (uniform in log(energy))
* - corrected for weight adjustment with sample_t=1 October 15, 2010
* in cases of using only a time slice of the full pulse
* - corrected for setting ttmin=tinmin if tinmin>ttmin October 15, 2010
* - added normalization parameters n_pulses and p_power October 18, 2010
* - normalizes with October 18, 2010
* (xwidth x yheight/0.1/0.12) *n_pulses *p_power/2.0
*
*
* Component: SNS_source_analytic
*
* %I
* Written by: F. X. Gallmeier
* Date: October 18, 2010
* Origin: SNS Oak Ridge National Laboratory
*
* A source that produces a time and energy distribution from
* parameterized SNS moderator files
*
* %D
* Produces a time-of-flight spectrum from SNS moderator files
* moderator files can be obtained from the author
* 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 parameters focus_xw focus_yh and dist
* provide the solid angle with which the beam line is served.
* The best practice is to set focus_xw and focus_yh to the width and height of the
* closest beam component, and dist as the distance from the moderator
* location to this component.
* (2) Be sure that Emin and Emax are within the limits of 1e-5 to 100 eV
* (3) the proton pulse length T determines short- and long-pulse mode
*
* Beamport definitions:
* BL11: a1Gw2-11-f5_fit_fit.dat
* BL14: a1Gw2-14-f5_fit_fit.dat
* BL17: a1Gw2-17-f5_fit_fit.dat
* BL02: a1Gw2-2-f5_fit_fit.dat
* BL05: a1Gw2-5-f5_fit_fit.dat
* BL08: a1Gw2-8-f5_fit_fit.dat
*
* %P
* Input parameters:
* filename: [] Filename of source data
* xwidth: [m] width of moderator (default 0.1 m)
* yheight: [m] height of moderator (default 0.12 m)
* dist: [m] Distance from source to the focusing rectangle (no default)
* focus_xw: [m] Width of focusing rectangle (no default)
* focus_yh: [m] Height of focusing rectangle (no default)
* Emin: [meV] minimum energy of neutrons (default 0.01 meV)
* Emax: [meV] maximum energy of neutrons (default 1e5 meV)
* tinmin: [us] minimum time of neutrons (default 0 us)
* tinmax: [us] maximum time of neutrons (default 2000 us)
* sample_E: [] sample energy from: (default 0)
* 0="probability distribution",
* 1="energy range",
* 2="lambda range",
* 3="log(energy) range"
* sample_t: [] sample t from: (default 0)
* 0="probability distribution",
* 1="time range"
* proton_T: [us] proton pulse length (0 us)
* p_power: [MW] proton power (default 2 MW)
* n_pulses: [] number of pulses (default 1)
*
* %E
*******************************************************************************/
DEFINE COMPONENT SNS_source_analytic
DEFINITION PARAMETERS ()
SETTING PARAMETERS (
string filename,
xwidth=0.10,
yheight=0.12,
dist,
focus_xw,
focus_yh,
Emin=0.01,
Emax=1.0e5,
tinmin=0.0,
tinmax=2000.0,
sample_E=0,
sample_t=0,
proton_T=0.0,
p_power=2.0,
n_pulses=1.0)
OUTPUT PARAMETERS ()
SHARE
%{
#ifdef OPENACC
/* there's no abort() on the GPU */
#define _ABORT()
#else
#define _ABORT() abort()
#endif
/*############################################################################################
#
# slowing-down spectrum and two Maxwellians joined by a modified Wescott function
#
# I(E) = I*1e12 * exp(-c/sqrt(E))
# * ( R1*E/(kT1)**2*exp(-E/kT1) + R2*E/(kT2)**2*exp(-E/kT2)
# + R3*E/(kT3)**2*exp(-(E/kT3)**b) + D(E)*rho(E)/E**(1-a) )
# with
# D(E) = 1/(1+(Ecut/E)**s)
# rho(E) = 1 + delta*exp(-x)(1 + x +0.5*x**2)
# x(E) = g*(E-2B); for E>2B
# = 0; for E<=2B
#
# constants:
# k = 1.3805e-23 J/K = 8.617e-5 eV/K
# B = 7.36e-3 eV
#
# parameters:
# I1
# c
# R1
# T1
# R2
# T2
# R3
# T3
# a
# b
# Ecut
# s
# delta
# g
#*/
#pragma acc routine seq
double spectral_function( double para[14], double E )
{
double c, R1, T1, R2, T2, R3, T3, a, b, Ecut, s, delta, g, Io;
double D, x, B, k, rho, arg1, arg2, arg3, arg4, arg5, arg6;
/* constants */
k = 8.617e-5;
B = 7.36e-3;
/* initialization of parameters */
c = para[0]; ;
R1 = para[1];
T1 = para[2];
R2 = para[3];
T2 = para[4];
R3 = para[5];
T3 = para[6];
a = para[7];
b = para[8];
Ecut = para[9];
s = para[10];
delta= para[11];
g = para[12];
Io = para[13];
/* evaluation of spectral function */
D = 1.0/(1+pow(Ecut/E,s));
x = 0.0;
if(E>2.0*B) {x = g*(E-2.0*B);}
rho = 1.0 + delta*exp(-x)*(1 + x +0.5*x*x);
arg1 = Io*1.0e12 * exp(-c/sqrt(E));
arg2 = R1*E/pow(k*T1,2) *exp(-E/(k*T1));
arg3 = R2*E/pow(k*T2,2) *exp(-E/(k*T2));
arg4 = R3*E/pow(k*T3,2) *exp(-pow(E/(k*T3),b));
arg5 = D*rho/pow(E,1-a);
arg6 =(arg1 * ( arg2 + arg3 + arg4 + arg5 ));
return(arg6);
}
/*############################################################################################
#
# prepares a vector of 1000 equiprobable energies in the range of Emin to Emax
# from the spectral distribution function
# returns the 1001 bin boundaries in array csfE
# returns the specctral integral in the range Emin to Emax in CItot
#
*/
double prepare_cumulative_spectral_function ( double *csfE, double Emin, double Emax,
double para[14] )
{
double E, E0, E1, CI, dCI, CItot, I0, I1, arg;
double u, umax, umin, du;
int i;
umax = log(Emax*1e-3);
umin = log(Emin*1e-3);
du = (umax - umin)*1e-5;
CItot=0.0;
E0=0.0;
I0 = spectral_function( para, Emin );
for(u=umin+du; u<=umax; u+=du) {
E1=exp(u);
I1 = spectral_function( para, E1 );
CItot += (E1-E0)*(I1+I0)*0.5;
I0 = I1;
E0 = E1;
}
printf("\n CItot = %12.4e\n", CItot);
CI = 0.0;
dCI = 1.0e-3;
E0 = Emin*1e-3;
u = umin;
I0 = spectral_function( para, E0 );
i=0;
csfE[0]= Emin*1e-3;
while (i<1000) {
E1 = exp(u+du);
I1 = spectral_function( para, E1 );
arg = (I1+I0)*0.5/CItot;
if(CI+arg*(E1-E0)>dCI) {
E = E0 + (dCI-CI)/arg;
i++;
csfE[i] = E;
// printf("E-%d = %12.4e\n", i, csfE[i]);
CI=0.0;
I0 = spectral_function( para, E );
E0 = E;
u = log(E);
}
else {
u = u+du;
CI += arg*(E1-E0);
I0 = I1;
E0 = E1;
}
}
return(CItot);
}
/*#############################################################################################
#
# Short-pulse Ikeda Carpenter energy-time brightness function
#
# f(E,t) = a/2 *( (1-R)(a*t)**2 *exp(-a*t)
# +2*R*a**2*b/(a-b)**3 *[ exp(-b*t) - exp(-a*t)*(1+(a-b)*t+0.5*(a-b)**2*t**2) ] }
# with
#
# t=t-to
#
# parameters:
#
# a
# b
# R
# to
#*/
#pragma acc routine seq
double f_sp (double para[4], double t)
{
double a, b, R, to;
double tt, eat, ebt, arg1, arg2, arg3, arg4;
a = para[0];
b = para[1];
R = para[2];
to= para[3];
tt = t - to*10.0;
if(tt<0.0) {tt=0.0;}
eat = exp(-a*tt);
ebt = exp(-b*tt);
arg1 = (1-R) *pow(a*tt,2) *eat;
arg2 = 2*R*a*a *b/pow(a-b,3);
arg3 = 1 + (a-b)*tt + 0.5*pow(a-b,2) *tt*tt;
arg4 = 0.5*a *( arg1 + arg2 *( ebt - eat*arg3 ) );
return(arg4);
}
/*#############################################################################################
#
# cummulative f_sp short pulse distribution function
# derived anlytically
#
# F(t) = (1-R)*(1-(gamma*exp(-a*(t-to))) for t-to >0
# *R *(1-delta *(exp(-b*(t-to)) -exp(-a*(t-to)) *(b/a)
# *(1 + (a-b)/a *(a*(t-to) +1) +((a-b)/a)**2 *gamma))
# = 0 for t-to <=0
#
# with
# gamma = 0.5*((a*(t-to))**2 +a*(t-to) +1
# delta = a**3/(a-b)**3
# tm = max(to,(t-T))
#
# f_sp(t) = F(tm) -F(t)
#
# parameters:
# a = proportial to scattering cross section
# b = decay constant of longest living flux eigenfunction
# R = fraction of slowing down term
# to = delay time
# T = proton pulse length
#*/
#pragma acc routine seq
double cummulative_f_sp (double para[4], double t)
{
double a, b, R, to;
double tt, eat, ebt, arg1, arg2, arg3, arg4, arg5;
double delta, gt, ambda, ambda2;
a = para[0];
b = para[1];
R = para[2];
to= para[3];
ambda = (a-b)/a;
ambda2= ambda*ambda;
delta = pow(1/ambda,3);
tt=t-to*10.0;
if(tt<0.0) {tt=0.0;}
gt = 0.5*pow(a*tt,2) + a*tt + 1.0;
eat = exp(-a*tt);
ebt = exp(-b*tt);
arg1 = 1.0-gt*eat;
arg2 = delta*(1 -ebt);
arg3 = delta*b/a*eat*(1.0 + ambda*(a*tt+1.0) + ambda2*gt);
arg4 = delta*b/a*(1+ambda+ambda2);
arg5 = (1.0-R)*arg1 + R*(arg2+arg3-arg4);
return(arg5);
}
/*#############################################################################################
#
# samples t from f_sp within interval [tmin,tmax]
#
# F(t) = (1-R)*(1-(gamma*exp(-a*(t-to))) for t-to >0
# *R *(1-delta *(exp(-b*(t-to)) -exp(-a*(t-to)) *(b/a)
# *(1 + (a-b)/a *(a*(t-to) +1) +((a-b)/a)**2 *gamma))
# = 0 for t-to <=0
#
# with
# gamma = 0.5*((a*(t-to))**2 +a*(t-to) +1
# delta = a**3/(a-b)**3
# tm = max(to,(t-T))
#
# f_sp(t) = F(tm) -F(t)
#
# parameters:
# a = proportial to scattering cross section
# b = decay constant of longest living flux eigenfunction
# R = fraction of slowing down term
# to = delay time
# T = proton pulse length
#*/
#pragma acc routine seq
double sample_t_from_f_sp (double para[4], double tmin, double tmax, double randd, double *p_t)
{
double arg;
double t, tm, tp, fm, fp, renorm, Imin, Imax;
int n;
double diff, eps;
eps = 1.0e-4;
Imax = cummulative_f_sp(para, tmax);
/*
printf(" tmax=%e",tmax);
printf(" Imax=%e",Imax);
printf("\n");
*/
Imin = cummulative_f_sp(para, tmin);
/*
printf(" tmin=%e",tmin);
printf(" Imin=%e",Imin);
printf("\n");
*/
renorm = 1.0 /(Imax-Imin);
tm = tmin;
fm = 0.0;
tp = tmax;
fp = 1.0;
n = 0;
diff = 1.0;
while(fabs(diff)>eps) {
n++;
if(n<10) {
t = 0.5*(tp+tm);
}
else {
t = tm +(tp-tm)/(fp-fm)*(randd-fm);
}
arg = cummulative_f_sp(para, t);
arg = (arg -Imin)*renorm;
diff = arg - randd;
if(n>50){
printf(" sample_t_from_f_sp exeeds 50 iterations!\n");
printf(" tmax=%01.4e tmin=%10.4e " , tmax, tmin);
printf(" Imax=%01.4e Imin=%10.4e arg4=%10.4e" , Imax, Imin, arg);
printf(" randd=%10.4e diff=%10.4e t=%10.4e\n", randd, diff, t);
_ABORT();
}
if(argrandd) {
tp = t;
fp = arg;
}
}
/* printf("tmin=%f tmax=%f Imin=%e Imax=%e t=%e n=%d\n", ttmin, ttmax, Imin, Imax, t, n); */
*p_t = Imax-Imin;
return(t);
}
/*#############################################################################################
#
# Long-pulse Ikeda Carpenter function
#
# proton pulse is a atep function with pulse length T (Heavyside function)
#
# folded with the Ikeada Carpenter Function f_sp
#
# F(t) = (-1/T)*(1-R)*(gamma*exp(-a*(t-to)) for t-to >0
# (-1/T)*delta *(exp(-b*(t-to)) -exp(-a*(t-to)) *(b/a)
# *(1 + (a-b)/a *(a*(t-to) +1) +((a-b)/a)**2 *gamma)
# = 0 for t-to <=0
#
# with
# gamma = 0.5*((a*(t-to))**2 +a*(t-to) +1
# delta = R*a**3/(a-b)**3
# tm = max(to,(t-T))
#
# f_lp(t) = F(tm) - F(t)
#
# parameters:
# a = proportial to scattering cross section
# b = decay constant of longest living flux eigenfunction
# R = fraction of slowing down term
# to = delay time
# T = proton pulse length
#*/
#pragma acc routine seq
double f_lp (double para[4], double t, double T)
{
double a, b, R, to;
double tt, eat, ebt, arg1, arg2, Imax, Imin, Itot, delta, ambda, ambda2, gt;
a = para[0];
b = para[1];
R = para[2];
to= para[3];
ambda = (a-b)/a;
ambda2= ambda*ambda;
delta = pow(ambda,-3);
/* upper time boundary */
tt = t - to*10.0;
if(tt<0.0) {tt=0.0;}
eat = exp(-a*tt);
ebt = exp(-b*tt);
gt = 0.5*pow(a*tt,2) + a*tt + 1.0;
arg1 = gt*eat;
arg2 = delta*(ebt - eat*b/a*(1.0 + ambda*(a*tt+1.0) + ambda2*gt));
Imax = (1.0-R)*arg1 + R*arg2;
/*
printf(" tt =%e",tt);
printf(" eatm=%e",eat);
printf(" ebtm=%e",ebt);
printf(" gtm=%e",gt);
printf(" arg1=%e",arg1);
printf(" arg2=%e",arg2);
printf(" Imax=%e",Imax);
printf("\n");
*/
/* lower time boundary */
tt = tt - T;
if(tt<0.0) {tt=0.0;}
eat = exp(-a*tt);
ebt = exp(-b*tt);
gt = 0.5*pow(a*tt,2) + a*tt + 1.0;
arg1 = gt*eat;
arg2 = delta*(ebt - eat*b/a*(1.0 + ambda*(a*tt+1.0) + ambda2*gt));
Imin = (1.0-R)*arg1 + R*arg2;
Itot = (Imin - Imax)/T;
/*
printf(" tto =%e",tt);
printf(" eat =%e",eat);
printf(" ebt =%e",ebt);
printf(" gt =%e",gt);
printf(" arg1=%e",arg1);
printf(" arg2=%e",arg2);
printf(" Imin=%e",Imin);
printf(" Itot=%e",Itot);
printf("\n");
*/
return(Itot);
}
/*#############################################################################################
#
# Ikeda carpenter function integrated twice over time
# base function used for calculating cumulative long-pulse time distribution function
#
# FF(t) = (-1/T)*(1-R)*(-1/a)*gamma*exp(-a*(t-to)) for t-to >0
# +(-1/T)*R*delta*(-1/b)*(exp(-b*(t-to))
# +(-1/T)*R*delta*(b/a**2)
# *(1+2*ambda + ambda*a*t + ambda**2*gamma)
# *exp(-a*(t-to))
# = 0 for t-to <=0
#
# with
# gamma = 0.5*(a*(t-to))**2 +2*a*(t-to) +3
# ambda =(a-b)/a
# delta = ambda**(-3)
#
# remember:
# F(t) = integral(f_sp(t))
# f_sp(t) = short-pulse emission time distribution
# f_lp(t) = integral_(T-t^t(f_sp(t))
# = (F(t)-F((t-T))/T
# hence:
# FF(t) = integral(F(t))
#
# parameters:
# a = proportial to scattering cross section
# b = decay constant of longest living flux eigenfunction
# R = fraction of slowing down term
# to = delay time
# T = proton pulse length
#*/
#pragma acc routine seq
double integral_integral_IC (double para[4], double t, double T)
{
double a, b, R, to;
double eat, ebt, arg1, arg2, arg3, arg4;
double tt, delta, gamma, ambda, ambda2;
a = para[0];
b = para[1];
R = para[2];
to= para[3];
tt = t-to*10.0;
if(tt<0.0) {tt=0.0;}
ambda = (a-b)/a;
ambda2 = ambda*ambda;
delta = pow(ambda,-3);
/* integrate 0 to T */
eat = exp(-a*tt);
ebt = exp(-b*tt);
gamma = 0.5*pow(a*tt,2) + 2.0*a*tt + 3.0;
arg1 = gamma/a*eat;
arg2 = delta/b*ebt;
arg3 = delta*b/a/a*(1.0+2.0*ambda+ambda*a*tt+ambda2*gamma)*eat;
arg4 = ((1.0-R)*arg1 + R*arg2 -R*arg3)/T;
/*
printf(" t=%e", t);
printf(" eatm=%e",eat);
printf(" ebtm=%e",ebt);
printf(" gtm=%e",gamma);
printf(" arg1=%e",arg1);
printf(" arg2=%e",arg2);
printf(" arg3=%e",arg3);
printf(" arg4=%e",arg4);
printf("\n");
*/
return(arg4);
}
/*#############################################################################################
#
# cummulative distribution function f_lp
# integrated analytically:from f_lp
#
# FF(t) = (-1/T)*(1-R)*(-1/a)*gamma*exp(-a*(t-to)) for t-to >0
# +(-1/T)*R*delta*(-1/b)*(exp(-b*(t-to))
# +(-1/T)*R*delta*(b/a**2)
# *(1+2*ambda + ambda*a*t + ambda**2*gamma)
# *exp(-a*(t-to))
# = 0 for t-to <=0
#
# with
# gamma = 0.5*(a*(t-to))**2 +2*a*(t-to) +3
# ambda =(a-b)/a
# delta = ambda**(-3)
# tm = max(to,(t-T))
#
# remember:
# F(t) = integral(f_sp(t))
# f_sp(t) = short-pulse emission time distribution
# f_lp(t) = integral_(T-t^t(f_sp(t))
# = (F(t)-F((t-T))/T
# hence:
# FF(t) = integral(F(t))
# CI(t) = integral_0^t(f_lp(t))
# = FF(t) - FF(tm) + f_lp(0)*(t-to)
# with
# tm=max(0,t-T)
#
# parameters:
# a = proportial to scattering cross section
# b = decay constant of longest living flux eigenfunction
# R = fraction of slowing down term
# to = delay time
# T = proton pulse length
#*/
#pragma acc routine seq
double cummulative_f_lp (double para[4], double t, double T)
{
double a, b, R, to;
double arg;
double tt, delta, ambda, ambda2, phi;
a = para[0];
b = para[1];
R = para[2];
to= para[3];
tt = t-to*10.0;
if(tt<0.0) {tt=0.0;}
ambda = (a-b)/a;
ambda2 = ambda*ambda;
delta = pow(ambda,-3);
phi = (1.0-R)+ R*delta*(1.0 - b/a*(1.0+ambda+ambda2));
/* integrate 0 to T */
tt=t-to*10.;
if(tt<0.0) {tt=0.0;}
if(tt-1e-5 ) {CImin = 0.0;}
if(CImax>0.0 && CImin>=0.0) {
renorm = 1.0 /(CImax-CImin);
}
else {
printf("troubble in sample_t_from_f_lp\n");
printf(" E=%e tmax=%e tmin=%e CImax=%e CImin=%e\n", E, tmax, tmin, CImax, CImin);
_ABORT();
}
tm = tmin;
fm = 0.0;
tp = tmax;
fp = 1.0;
n = 0;
diff = 1.0;
while(fabs(diff)>eps) {
n++;
if(n<10) {
t = 0.5*(tp+tm);
}
else {
t = tm +(tp-tm)/(fp-fm)*(randd-fm);
}
arg = cummulative_f_lp(para,t,T);
arg = (arg -CImin)*renorm;
diff = arg - randd;
if(n<50){
}
else {
if(fabs(diff)randd) {
tp = t;
fp = arg;
}
}
/* printf("tmin=%f tmax=%f Imin=%e Imax=%e t=%e n=%d\n", ttmin, ttmax, CImin, CImax, t, n); */
*p_t = CImax-CImin;
return(t);
}
/*#############################################################################################
#
# Pade type fitting function
#
# f(x) = log(a*x**b*(1+c*x+d*x**2+(x/f)**g)/(1+h*x+i*x**2+(x/j)**k))
#
# with
# a,b,c,d,f,g,h,i,j,k
#
#
#*/
#pragma acc routine seq
double pade_function (double para[10], double E)
{
double a, b, c, d, f, g, h, i, j, k;
a = para[0];
b = para[1];
c = para[2];
d = para[3];
f = para[4];
g = para[5];
h = para[6];
i = para[7];
j = para[8];
k = para[9];
return(a*pow(E,b)*(1+c*E+d*E*E+pow(E/f,g))/(1+h*E+i*E*E+pow(E/j,k)));
}
%}
DECLARE
%{
double hdiv;
double vdiv;
double p_in;
double CL;
double CU;
double para_sp[14];
double para_a[10];
double para_b[10];
double para_R[10];
double para_to[10];
double csfE[1001];
double CItot;
%}
INITIALIZE
%{
char line[1000];
int linelength=1000;
int k, kk;
FILE *fp, *fopen();
fp = fopen(filename,"r");
if (fp==NULL){
printf("Error opening file\n");
}
else {
printf("%s\n","File opened...");
}
/* spectral parameters */
while( (fgets(line,linelength,fp) != NULL) && (strchr(line,'#') != NULL) ){
printf("%s",line);
}
kk=sscanf(line," %le %le %le %le %le %le %le %le %le %le %le %le %le %le", ¶_sp[0], ¶_sp[1], ¶_sp[2],
¶_sp[3], ¶_sp[4], ¶_sp[5], ¶_sp[6], ¶_sp[7], ¶_sp[8], ¶_sp[9], ¶_sp[10],
¶_sp[11], ¶_sp[12], ¶_sp[13]);
/* a parameter of emission time distribuation */
while( (fgets(line,linelength,fp) != NULL) && (strchr(line,'#') != NULL) ){
printf("%s",line);
}
kk=sscanf(line," %le %le %le %le %le %le %le %le %le %le", ¶_a[0], ¶_a[1], ¶_a[2],
¶_a[3], ¶_a[4], ¶_a[5], ¶_a[6], ¶_a[7], ¶_a[8], ¶_a[9]);
if(kk!=10) {
printf("para_a: insufficient number of data entries read\n");
_ABORT();
}
/* b parameter of emission time distribuation */
while((fgets(line,linelength,fp)!=NULL)&&(strchr(line,'#')!=NULL)){
printf("%s",line);
}
kk=sscanf(line," %le %le %le %le %le %le %le %le %le %le", ¶_b[0], ¶_b[1], ¶_b[2],
¶_b[3], ¶_b[4], ¶_b[5], ¶_b[6], ¶_b[7], ¶_b[8], ¶_b[9]);
if(kk!=10) {
printf("para_b: insufficient number of data entries read\n");
_ABORT();
}
/* R parameter of emission time distribuation */
while((fgets(line,linelength,fp)!=NULL)&&(strchr(line,'#')!=NULL)){
printf("%s",line);
}
kk=sscanf(line," %le %le %le %le %le %le %le %le %le %le", ¶_R[0], ¶_R[1], ¶_R[2],
¶_R[3], ¶_R[4], ¶_R[5], ¶_R[6], ¶_R[7], ¶_R[8], ¶_R[9]);
if(kk!=10) {
printf("para_R: kk=%d insufficient number of data entries read\n", kk);
_ABORT();
}
/* to parameter of emission time distribuation */
while((fgets(line,linelength,fp)!=NULL)&&(strchr(line,'#')!=NULL)){
printf("%s",line);
}
kk=sscanf(line," %le %le %le %le %le %le %le %le %le %le", ¶_to[0], ¶_to[1], ¶_to[2],
¶_to[3], ¶_to[4], ¶_to[5], ¶_to[6], ¶_to[7], ¶_to[8], ¶_to[9]);
if(kk!=10) {
printf("para_R: insufficient number of data entries read\n");
_ABORT();
}
printf("\n");
for(k=0; k<14; k++){
printf(" para_sp(%d)=%e\n",k,para_sp[k]);
}
printf("\n");
for(k=0; k<10; k++){
printf(" para_a(%d)=%e\n",k,para_a[k]);
}
printf("\n");
for(k=0; k<10; k++){
printf(" para_b(%d)=%e\n",k,para_b[k]);
}
printf("\n");
for(k=0; k<8; k++){
printf(" para_R(%d)=%e\n",k,para_R[k]);
}
printf("\n");
for(k=0; k<10; k++){
printf(" para_to(%d)=%e\n",k,para_to[k]);
}
if(sample_E==0) {
CItot = prepare_cumulative_spectral_function ( csfE, Emin, Emax, para_sp );
if(CItot<=0.0) {
printf(" choice of Emin and Emax gives zero neutron intensity\n");
_ABORT();
}
// for(i=0;i<=1000;i++){
// printf("E-%d = %12.4e\n", i, csfE[i]);
// }
}
/* some checks */
if(dist<=0.0) {
printf(" dist must be greater zero\n");
_ABORT();
}
else if(focus_xw<=0.0) {
printf(" focus_xw must be greater zero\n");
_ABORT();
}
else if(focus_yh<=0.0) {
printf(" focus_yh must be greater zero\n");
_ABORT();
}
else if(p_power<=0) {
printf(" p_power must be greater zero\n");
_ABORT();
}
else if(n_pulses<=0) {
printf(" n_pulses must be greater zero\n");
_ABORT();
}
/* Normalization */
// Calculate solid Angle
p_in = focus_xw*focus_yh/(dist*dist);
// Normalize to viewed area
p_in *= xwidth*yheight/0.1/0.12;
// Normalize to proton power
p_in *= p_power/2.0;
// Normalize to number of pulses
p_in *= n_pulses;
/* constants for conversion into wavelength and lethargy units */
CL = (1.0/sqrt(Emin*1e-3)-1.0/sqrt(Emax*1e-3))/(Emax-Emin)*1e3;
CU = (log(Emax*1e-3)-log(Emin*1e-3))/(Emax-Emin)*1e3;
#ifdef USE_MPI
p_in /= mpi_node_count;
#endif
p_in /= mcget_ncount();
%}
TRACE
%{
double theta,phi,v,E,Eval,tval;
double hdivmin,hdivmax,vdivmin,vdivmax;
double p_E, p_t, para_t[4], p_sa, randd, ttmin, ttmax;
int i;
p=p_in;
z=0;
x = (rand01()-0.5)*xwidth; /* choose points uniformly distributed on the source */
y = (rand01()-0.5)*yheight;
hdivmax=atan((focus_xw/2.0-x)/dist);
hdivmin=atan(-(focus_xw/2.0+x)/dist);
vdivmax=atan((focus_yh/2.0-y)/dist);
vdivmin=atan(-(focus_yh/2.0+y)/dist);
theta = hdivmin + (hdivmax-hdivmin)*rand01(); /* Small angle approx. */
phi = vdivmin + (vdivmax-vdivmin)*rand01();
hdiv=theta;
vdiv=phi;
/* find E value corresponding to randomly in E range
sample_E=0 draws uniformly from probability distribution
sample_E=1 draws uniformly in energy range
sample_E=2 draws uniformly from wavelength range
sample_E=3 draws uniformly from lethargy range
*/
if(sample_E==0) {
randd = rand01();
i = randd*1e3;
Eval = csfE[i] + (csfE[i+1]-csfE[i])*(randd*1e3-i);
p_sa = 1.0;
p_E = CItot;
}
else if(sample_E==1){
Eval=(Emin+(Emax-Emin)*rand01())*1e-3;
p_sa = 1.0;
p_E = spectral_function(para_sp, Eval) * (Emax-Emin)*1e-3;
}
else if(sample_E==2) {
Eval = pow(1.0/(1.0/sqrt(Emax) +(1.0/sqrt(Emin)-1.0/sqrt(Emax))*rand01()),2)*1e-3;
p_sa = CL*2.0*pow(Eval,1.5);
p_E = spectral_function(para_sp, Eval) * (Emax-Emin)*1e-3;
}
else if(sample_E==3) {
Eval = exp(log(Emin) +(log(Emax)-log(Emin))*rand01())*1e-3;
p_sa = CU*Eval;
p_E = spectral_function(para_sp, Eval) * (Emax-Emin)*1e-3;
}
else {
printf("sample_E allows only values 0, 1, 2 or 3!\n");
_ABORT();
}
/* determine tval from emisstime distribution */
para_t[0] = pade_function(para_a, Eval);
para_t[1] = pade_function(para_b, Eval);
para_t[2] = pade_function(para_R, Eval);
para_t[3] = pade_function(para_to, Eval);
// printf("para1=%f para_b=%f para_R=%f para_to=%f\n", para_t[0], para_t[1], para_t[2], para_t[3]);
/* find t value corresponding to random probability
sample_t=0: draw from probability interval [0,1] and calculate corresponding tval
sample_t=1: draw tval from interval [0,20/beta] for short-pulse
from interval [0,20/beta+proton_T] for long-pulse
*/
if(sample_t==0) {
randd = rand01();
if(proton_T<=0.0) {
ttmax = 20.0/para_t[1]+para_t[3]*10.;
if(ttmax>tinmax) {ttmax=tinmax;}
ttmin = para_t[3]*10.;
if(tinmin>ttmin) {ttmin = tinmin;}
tval = sample_t_from_f_sp (para_t, ttmin, ttmax, randd, &p_t);
}
else {
ttmax = proton_T+20.0/para_t[1]+para_t[3]*10.;
if(ttmax>tinmax) {ttmax=tinmax;}
ttmin = para_t[3]*10.;
if(tinmin>ttmin) {ttmin = tinmin;}
tval = sample_t_from_f_lp (para_t, Eval, ttmin, ttmax, proton_T, randd, &p_t);
}
}
else if(sample_t==1) {
if(proton_T<=0.0) {
ttmax = 20.0/para_t[1]+para_t[3]*10.;
if(tinmaxttmin) {ttmin = tinmin;}
tval=ttmin+(ttmax-tinmin)*rand01();
p_t = f_sp(para_t, tval);
}
else {
ttmax = proton_T+20.0/para_t[1]+para_t[3]*10.;
if(tinmaxttmin) {ttmin = tinmin;}
tval=ttmin+(ttmax-ttmin)*rand01();
p_t = f_lp(para_t, tval, proton_T);
}
p_t *= (ttmax-ttmin);
}
else {
printf("sample_t allows only values 0 or 1!\n");
}
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*=p_E*p_sa*p_t;
//printf("Eval=%e tval=%e p_E=%e p_t=%e p_sa=%e p=%e\n", Eval, tval, p_E, p_t, p_sa, p);
%}
FINALLY
%{
%}
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