a TOF source
H.Kadowaki
kadowaki at phys.metro-u.ac.jp
Thu Sep 14 19:01:44 CEST 2000
Hello McStas users,
I have written a TOF source component producing pulse shapes of
the Ikeda-Carpenter function.
It emulates the energy spectrum and pulse shape of a KEK-JAERI 1MW
pulse-source to be constructed, decoupled supercritical hydrogen
moderator (12x12x5 cm^3) at 20 K with an H2O premoderator in a Pb
reflector (decoupling energy = 1 eV).
I am sending 3 files.
(1) Modr_dcHPb_flux.comp
TOF source component
(2) modr_dcHPb.c
functional forms of energy spectrum and pulse shape are
described.
(3) T_Modr_dcHPb_flux.instr
example of testing Modr_dcHPb_flux.comp
Its results are not so bad.
Regards
H. Kadowaki
Department of Physics, Tokyo Metropolitan University,
Hachioji-shi, Tokyo 192-0397, Japan
E-mail: kadowaki at phys.metro-u.ac.jp
############### remove this line ##################
/******************begin Modr_dcHPb_flux.comp ******/
/*
* McStas, the neutron ray-tracing package
* Maintained by Kristian Nielsen and Kim Lefmann,
* Copyright 1997-2000 Risoe National Laboratory, Roskilde, Denmark
*
* %I
* Written by: H KADOWAKI kadowaki at phys.metro-u.ac.jp
* Date: 2000/9/8
* Version: 0.2
* Origin: modified from Moderator.comp of McStas1.3 release
*
* A pulsed source of a decoupled H_2 moderator for time-of-flight.
* This emulates the energy spectrum and pulse shape of
* a KEK-JAERI 1MW pulse-source to be constructed,
* decoupled supercritical hydrogen moderator (12x12x5 cm^3) at 20 K
* with an H2O premoderator in a Pb reflector (decoupling energy = 1 eV).
*
* %D
* Produces a time-of-flight spectrum, with a energy distribution
* using fitted functions written in the Ikead & Carpenter paper.
* The functions espectr(E) and pulse_shape(E, t) are defined
* in modr_dcHPb.c
*
* Since the flux of the source is given in an absolute scale and
* the weight p is adjusted properly, the computed intensity I is
* counting rate in unit time [counts/(sec MW)].
*
* To use this component, compile with modr_dcHPb.c as follows.
* mcstas test.instr
* gcc -O -o test.ex test.c modr_dcHPb.c -lm
*
* %P
* Input parameters:
*
* xws: (m) Width of source
* yhs: (m) Height of source
* angle (degree) Angle of rotation by which the front face is rotated
* around the y-axis (angle=0 : xy-plane)
* Emin: (meV) Lower edge of energy distribution
* Emax (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
*
* %E
***************************************************************************/
DEFINE COMPONENT Modr_dcHPb_flux
DEFINITION PARAMETERS (xws, yhs, angle, Emin, Emax, dist, xw, yh)
SETTING PARAMETERS ()
OUTPUT PARAMETERS (hdiv,vdiv,p_in,angle_rad,xxx_low_limit,xx,del_E,n_E
,tauEtab)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
DECLARE
%{
#define TABLE_LENGTH_E 1024
double hdiv,vdiv;
double p_in,angle_rad,xxx_low_limit;
double xx[TABLE_LENGTH_E],del_E;
int n_E;
double tauEtab[TABLE_LENGTH_E];
%}
INITIALIZE
%{
int i;
double ee,denom,ang;
extern double espectr( double ee);
extern double pulse_shape(double ee, double ttt);
hdiv = atan(xw/(2.0*dist));
vdiv = atan(yh/(2.0*dist));
ang=angle; angle_rad = ang/45.0*atan(1.0);
/* set xx[i] */
n_E=TABLE_LENGTH_E-1; /* number of division */
del_E=(Emax-Emin)/n_E; /* step of E */
denom=0.0;
for(i=0; i<n_E; i++)
{
ee=Emin+del_E*(i+0.5);
denom=denom + espectr(ee);
xx[i+1]=denom;
}
xx[0]=0.0; xx[n_E]=1.0;
for(i=1; i<n_E; i++)
{ xx[i]=xx[i]/denom; }
/* xx[i] = [\int_{Emin}^{E } espectr(e) de ]
/[\int_{Emin}^{Emax} espectr(e) de ]
E = Emin + i*del_E */
denom=denom*del_E; /*[\int_{Emin}^{Emax} espectr(E) dE ]*/
/* set p_in
p_in = p_0
= [\int_{Emin}^{Emax} espectr(E) dE ]
* [(dA) (d\Omega) (dt=1 sec)]/N_{sim}
dA = (xws*yhs*1.0e4) [cm^2]
d\Omega = (4.0*hdiv*vdiv) [str]
N_{sim} = mcget_ncount() = number of neutron generated
\int_{Emin}^{Emax} espectr(E) dE = (denom * 1.0e-3)
[neutrons/(cm^2 sr sec MW)]
see manual of Source_flux_lambda
*/
p_in = (xws*yhs*1.0e4)*(4.0*hdiv*vdiv); /* cm^2 str Small angle approx. */
p_in = p_in * (denom * 1.0e-3) / mcget_ncount();
xxx_low_limit= 5.0/ mcget_ncount();
/* set tauEtab[i]= decay time of neutrons with E = Emin+del_E*i
= tau [microsec] defined by
pulse_shape(E, tau) = max[ pulse_shape(E, tt) ]/2.71 */
for(i=0; i<=n_E; i++)
{
double tt,peak,tpeak,yy,tau;
ee=Emin+del_E*i; /* Emin<=ee<=Emax */
peak=0.0; tpeak=0.0;
for(tt=0.25; tt<10000.0; tt=tt+0.25)
{
yy=pulse_shape(ee, tt);
if(peak < yy )
{peak=yy; tpeak=tt;}
else
{
if( peak > (yy*2.71) )
{tau=tt; break;}
}
}
/* printf(" %g %g %g \n",ee,tpeak,tau); */
tauEtab[i]=tau;
}
%}
TRACE
%{
double theta0,phi0,theta,phi,v,E;
double distp,hdivp,vdivp;
double xxx,tauE;
int il,ih,im;
extern double pulse_shape(double ee, double ttt);
p=p_in;
x = xws*(-0.5+rand01()) ; /* Choose point on source */
y = yhs*(-0.5+rand01()) ; /* with uniform distribution. */
z = - x * sin(angle_rad); /* rotation around the y-axis */
x = x * cos(angle_rad);
distp = dist-z;
theta0 = -atan(x/distp); /* Angles to aim at target centre */
phi0 = -atan(y/distp);
hdivp = atan(xw/(2.0*distp));
vdivp = atan(yh/(2.0*distp));
theta = theta0 + hdivp*randpm1(); /* Small angle approx. */
phi = phi0 + vdivp*randpm1();
/* set E by energy distribution of espectr(E) */
xxx=rand01(); /* xxx -> E */
/* xxx = [\int_{Emin}^{E } espectr(e) de ]
/[\int_{Emin}^{Emax} espectr(e) de ] */
/* xx[il] <= xxx < xx[ih=il+1] find out il,ih */
il=0; ih=n_E;
for(im=il+(ih-il)/2; (ih-il)>1; im=il+(ih-il)/2)
{
if(xxx<xx[im]) { ih=im; }
else { il=im; }
}
/* set E by linear interpolation using xx[i]
xx[i] = [\int_{Emin}^{E } espectr(e) de ]
/[\int_{Emin}^{Emax} espectr(e) de ]
E = Emin + i*del_E, i=il+(xxx-xx[il])/(xx[ih]-xx[il]) */
E=il+(xxx-xx[il])/(xx[ih]-xx[il]);
E=Emin+del_E*E;
v = SE2V*sqrt(E);
vz = v*cos(phi)*cos(theta); /* Small angle approx. */
vy = v*sin(phi);
vx = v*cos(phi)*sin(theta);
/* set tauE = decay time of neutron with energy E */
il=(E-Emin)/del_E;
if(il >= n_E) {il=n_E-1;}
if(il < 0) {il=0;}
xxx=(E-Emin-del_E*il)/del_E;
xxx=(tauEtab[il+1]-tauEtab[il])*xxx;
tauE=tauEtab[il]+xxx; /* linear interpolation */
/* set t using distribution exp(-t/tauE)/tauE */
do{
xxx=rand01(); /* 1-rand01() */
/* if(xxx <= xxx_low_limit )
{printf(" xxx <= %g, xxx= %g \n",xxx_low_limit,xxx);}
discard small xxx , i.e., very large t, p */
}
while( xxx <= xxx_low_limit );
t=-tauE*log(xxx); /* xxx -> t */
/* set p by weight multiplication
weight = \pi_j = pulse_shape(E, t) / [exp(-t/tauE)/tauE]
= pulse_shape(E, t)*tauE /(xxx=1-rand01()) */
xxx=tauE/xxx;
xxx= pulse_shape(E, t)*xxx; /* weight \pi_j */
p=p*xxx;
t=t*1.0e-6; /* microsec -> sec */
%}
/* MCDISPLAY
%{
double xmin,xmax,ymin,ymax,zmin,zmax;
xmax=xws*0.5*cos(angle_rad);xmin=-xmax;
ymax=yhs*0.5;ymin=-ymax;
zmax=-xws*0.5*sin(angle_rad);zmin=-zmax
magnify("xy");
multiline(5, (double)xmin, (double)ymin, (double)zmin,
(double)xmax, (double)ymin, (double)zmax,
(double)xmax, (double)ymax, (double)zmax,
(double)xmin, (double)ymax, (double)zmin,
(double)xmin, (double)ymin, (double)zmin);
%} */
END
/******************end Modr_dcHPb_flux.comp ******/
############### remove this line ##################
/** begin modr_dcHPb.c **/
/* written by H KADOWAKI kadowaki at phys.metro-u.ac.jp
version 0.2 2000/9/8
*/
/*
Energy spectrum and pulse shape of
Decoupled supercritical hydrogen moderator (12 x 12 x 5 cm^3) at 20 K
with an H2O premoderator in a Pb reflector (decoupling energy = 1 eV)
Data were provided by simulations by TESHIGAWARA(2000/8/3).
To emulate these data in MCSTAS, they are fitted
using functions of the Ikeda & Carpenter (IC) paper, and
resulting functional forms are described as.
espectr(E) and pulse_shape(E, t)
Ref: S. Ikeda and J.M. Carpenter: Nucl. Instr. Meth. A239 (1985) 536.
espectr(E) = neutron flux as a function of energy
eq. (18) of IC [neutrons/(cm^2 sr eV sec MW)]
E [meV]
pulse_shape(E, t) = pulse shape of neutrons with energy=E
eq. (10') (11) (14) of IC
E [meV] t [microsec]
\begin{eqnarray*}
espectr(E) &=& \bar{i}(E) \\
&=&
\bar{i}_{Th} \frac{E}{E_{T}^{2}} \exp(-E/E_{T})
+ \bar{i}_{epi} \frac{\Delta(E)}{\Delta(E_0)}
\frac{1}{E} \left( \frac{E}{E_0} \right)^{\alpha} \\
\Delta(E) &=& [1 + \exp (a \lambda -b ) ]^{-1}
\end{eqnarray*}
\begin{eqnarray*}
time\_struct(E, t) &=& \psi(v, t) \\
&=& \frac{1}{2} (1-R) \alpha^3 t^2 \exp(- \alpha t) \\
&+& R \frac{\alpha^3 \beta}{(\alpha-\beta)^3}
\{
\exp(- \beta t) -
\exp(- \alpha t)
[1+(\alpha-\beta) t + \frac{1}{2} (\alpha-\beta)^2 t^2]
\} \\
\alpha &=& v \Sigma \\
\Sigma &=& (S_1^2 + S_2^2 \lambda^2 )^{1/2} \\
R &=& \exp (-E/E_0)
\end{eqnarray*}
To use these in MCSTAS, they should be declared
as external functions and compile, for example, as follows.
$mcstas test.instr
$gcc -O -o test.ex test.c modr_dcHPb.c -lm
*/
#include <stdio.h>
#include <math.h>
/* Energy spectrum */
double espectr(double ee)
{
double e0=1000.0,wlen0=0.2860057;
/* e0 = E_0 = 1000 [meV],
wlen0 = \lambda_0 [Angstrom] of E_0 */
double pa1=0.9695810E+15; /* \bar{i}_{Th} */
double pa2= 4.100054 ; /* E_T [meV] */
double pa3=0.2948123E+15; /* \bar{i}_{epi} */
double pa4= 10.29586 ; /* a [Angstrom^{-1}] */
double pa5= 26.28562 ; /* b [1] */
double pa6=0.2965801E-01; /* \alpha [1] */
/* pa1-pa6 are based on fit in range 1< E < 1000 meV */
double xx,xxx,yyy,wlen;
/* ee = energy [meV], wlen = wave length [Angstrom]
xx,xxx,yyy = work variables */
xxx=pa1/(pa2*pa2);
xx=-ee/pa2;
xxx=xxx*ee*exp(xx);
wlen=9.04429/sqrt(ee);
yyy=pa3*(1.0+exp(pa4*wlen0-pa5));
yyy=yyy/(1.0+exp(pa4*wlen -pa5));
yyy=yyy/ee;
yyy=yyy*pow((ee/e0),pa6);
return (xxx+yyy);
}
/* pulse shape Ikeda & Cartpenter function */
double pulse_shape(double ee, double ttt)
/* ee = energy [meV], ttt = time [microsec]*/
{
double pa1= .9669852 ; /* S_1 [cm^{-1}] */
double pa2= .1525358 ; /* S_2 [cm^{-1} Angstrom^{-1}] */
double pa3=0.1200498E-01; /* \beta [microsec^{-1}]*/
double pa4= 20.94545 ; /* E_0 of R=exp(-E/E_0) [meV]*/
/* based on fit in range 1< E < 400 meV */
double pi=3.141592654;
double velc,wlen,sigma,alpha,beta,amb,RR;
double xxx,zz,yyy,yyy1,uu,vv;
/* velc = velocity, wlen= wave length, sigma=\Sigma,
alpha=\alpha, beta=\beta, amb=\alpha-\beta, RR=R,
xxx,zz,yyy,yyy1,uu,vv= work variables */
if(ttt < 0.0)
{ xxx=0.0; }
else
{
xxx=sqrt(ee/2.072);
wlen=(2.0*pi)/xxx; /* wlen=\lambda [Angstrom]*/
velc=xxx*0.06296; /* velc=v=velocity [cm/micro sec] */
sigma=sqrt(pa1*pa1+(pa2*pa2*wlen*wlen)); /* sigma=\Sigma [cm^{-1}]*/
alpha=sigma*velc;
RR=-ee/pa4;
RR=exp(RR); /* R=exp(-E/E_0) */
xxx=0.5*alpha*(1.0-RR);
zz=alpha*ttt;
xxx=xxx*(zz*zz)*exp(-zz); /* 1st term of eq.(10') */
beta=pa3;
amb=alpha-beta;
yyy=alpha/amb;
yyy=yyy*yyy*yyy;
yyy=yyy*RR*beta;
uu=beta*ttt;
vv=amb*ttt;
yyy1=1.0+vv+0.5*vv*vv;
yyy1=exp(-uu)-exp(-zz)*yyy1;
yyy=yyy*yyy1; /* 2nd term of eq.(10') */
xxx=xxx+yyy;
}
return (xxx);
}
/** end modr_dcHPb.c **/
############### remove this line ##################
DEFINE INSTRUMENT T_Modr_dcHPb_flux (Elow,Ehigh)
/* T_Modr_dcHPb_flux.instr
test of Modr_dcHPb_flux.comp
mcstas T_Modr_dcHPb_flux.instr
gcc -O -o T_Modr_dcHPb_flux.ex T_Modr_dcHPb_flux.c modr_dcHPb.c -lm
*/
TRACE
COMPONENT arm1=Arm()
AT (0,0,0) ABSOLUTE
COMPONENT mod1=Modr_dcHPb_flux(
xws=0.01, yhs=0.01, angle=0.0,
Emin=Elow, Emax=Ehigh,
dist=1.0, xw=0.01, yh=0.01)
AT (0,0,0) RELATIVE arm1
COMPONENT mon1= TOF_monitor(
xmin=-0.005,xmax=0.005,ymin=-0.005,ymax=0.005,
nchan=1000,dt=0.5,
filename="tstr.T.Mod.dH.dat3")
AT (0,0,0.0001) RELATIVE arm1
COMPONENT mon2= E_monitor(
xmin=-0.005,xmax=0.005,ymin=-0.005,ymax=0.005,
Emin=Elow,Emax=Ehigh,
nchan=2000,
filename="espc.T.Mod.dH.dat3")
AT (0,0,1.001) RELATIVE arm1
END
More information about the mcstas-users
mailing list