McStas 1.5: bad sources
Emmanuel Farhi
farhi at ill.fr
Mon Oct 29 15:17:38 CET 2001
Hello McStas users,
I've found an error in all the 4 'standard' mcstas sources. None of them
work. Bug was really obvious.
They generate neutrons with speed about 1 m/s. I like UCN sources, but
usualy I prefer them warmer.
I here attach 'good' sources. You may also use the 'Source_gen' that can
make the jobs of all these (and works perfectly).
The Source_Maxwell_3 is ok. Source_Maxwell does not focus correctly on
the target.
Invitation:
These modifications are part of the McStas version 1.6-ill
<http://www.ill.fr/tas/mcstas/> (that one from which Per-Olof says "The
parallel development of McStas at ILL is very unfortunate". Wow ! you're
not very kind Per-Olof, why ?)
--
What's up Doc ?
--------------------------------------------
Emmanuel FARHI, http://www.ill.fr/tas/people/Farhi.html \|/ ____ \|/
CS-Group ILL4/156, Institut Laue-Langevin (ILL) Grenoble ~@-/ oO \-@~
6 rue J. Horowitz, BP 156, 38042 Grenoble Cedex 9,France /_( \__/ )_\
Work :Tel (33/0) 4 76 20 71 35. Fax (33/0) 4 76 48 39 06 \__U_/
-------------- next part --------------
/*******************************************************************************
*
* McStas, the neutron ray-tracing package: Source_flat_lambda.comp
* Copyright 1997-2001 Risoe National Laboratory, Roskilde, Denmark
*
* %I
* Written by: Kristian Nielsen and Kim Lefmann
* Date: 1999
* Version: $Revision: 1.5 $
* Origin: McStas 1.5
* Modified by: Kim Lefmann, October 8, 2001
*
* Neutron source with flat wavelength spectrum and arbitrary flux.
*
* %D
* The routine is a circular neutron source, which aims at a square target
* centered at the beam (in order to improve MC-acceptance rate). The angular
* divergence is then given by the dimensions of the target. The neutron
* wavelength is uniformly distributed between lambda_0 - d_lambda and
* lambda_0 + d_lambda.
*
* Example: Source_flat_lambda(radius=0.1, dist=2, xw=.1, yh=.1,
* lambda_0=2.36, d_lambda=.3)
*
* %P
* radius: (m) Radius of circle in (x,y,0) plane where neutrons
* are generated.
* dist: (m) Distance to target along z axis.
* xw: (m) Width(x) of target
* yh: (m) Height(y) of target
* lambda_0: (AA) Mean wavelength of neutrons.
* d_lambda: (AA) Wavelength spread of neutrons.
*
* %E
*******************************************************************************/
DEFINE COMPONENT Source_flat_lambda
DEFINITION PARAMETERS ()
SETTING PARAMETERS (radius, dist, xw, yh, lambda_0, d_lambda)
OUTPUT PARAMETERS (pmul)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
DECLARE
%{
double pmul;
%}
INITIALIZE
%{
pmul=1.0/(mcget_ncount()*4*PI);
%}
TRACE
%{
double chi,lambda,v,r,xf,yf,rf,dx,dy;
t=0;
z=0;
chi=2*PI*rand01(); /* Choose point on source */
r=sqrt(rand01())*radius; /* with uniform distribution. */
x=r*cos(chi);
y=r*sin(chi);
xf = 0.5*xw*randpm1(); /* Choose focusing position uniformly */
yf = 0.5*yh*randpm1();
dx = xf-x;
dy = yf-y;
rf = sqrt(dx*dx+dy*dy+dist*dist);
p = pmul*xw*yh*dist/(rf*rf*rf); /* target area * cos(phi)/rf^2 */
lambda = lambda_0+d_lambda*randpm1();
v = K2V*(2*PI/lambda);
vz=v*dist/rf;
vy=v*dy/rf;
vx=v*dx/rf;
%}
MCDISPLAY
%{
magnify("xy");
circle("xy",0,0,0,radius);
%}
END
-------------- next part --------------
/*******************************************************************************
*
* McStas, the neutron ray-tracing package: Source_flat.comp
* Copyright 1997-2001 Risoe National Laboratory, Roskilde, Denmark
*
* %I
* Written by: Kim Lefmann
* Date: October 30, 1997
* Modified by: KL, October 4, 2001
* Version: $Revision: 1.9 $
* Origin: McStas 1.5
*
* A circular neutron source with flat energy spectrum and arbitrary flux
*
* %D
* The routine is a circular neutron source, which aims at a square target
* centered at the beam (in order to improve MC-acceptance rate). The angular
* divergence is then given by the dimensions of the target.
* The neutron energy is uniformly distributed between E0-dE and E0+dE.
*
* Example: Source_flat(radius=0.1, dist=2, xw=.1, yh=.1, E0=14, dE=2)
*
* %P
* radius: (m) Radius of circle in (x,y,0) plane where neutrons
* are generated.
* dist: (m) Distance to target along z axis.
* xw: (m) Width(x) of target
* yh: (m) Height(y) of target
* E0: (meV) Mean energy of neutrons.
* dE: (meV) Energy spread of neutrons.
*
* %E
*******************************************************************************/
DEFINE COMPONENT Source_flat
DEFINITION PARAMETERS ()
SETTING PARAMETERS (radius, dist, xw, yh, E0, dE)
OUTPUT PARAMETERS ()
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
DECLARE
%{
double pmul;
%}
INITIALIZE
%{
pmul=1.0/(mcget_ncount()*4*PI);
%}
TRACE
%{
double chi,E,v,r, xf, yf, rf, dx, dy;
t=0;
z=0;
chi=2*PI*rand01(); /* Choose point on source */
r=sqrt(rand01())*radius; /* with uniform distribution. */
x=r*cos(chi);
y=r*sin(chi);
xf = 0.5*xw*randpm1(); /* Choose focusing position uniformly */
yf = 0.5*yh*randpm1();
dx = xf-x;
dy = yf-y;
rf = sqrt(dx*dx+dy*dy+dist*dist);
p = pmul*xw*yh*dist/(rf*rf*rf); /* target area * cos(phi)/rf^2 */
E=E0+dE*randpm1(); /* Choose from uniform distribution */
v=sqrt(E)*SE2V;
vz=v*dist/rf;
vy=v*dy/rf;
vx=v*dx/rf;
%}
MCDISPLAY
%{
magnify("xy");
circle("xy",0,0,0,radius);
%}
END
-------------- next part --------------
/**********************************************************************
*
* McStas, the neutron ray-tracing package: Source_flux.comp
* Copyright 1997-2001 Risoe National Laboratory, Roskilde, Denmark
*
* Component: Source_flux
*
* %Identification
* Written by: Kristian Nielsen
* Date: 1998
* Version: $Revision: 1.5 $
* Origin: McStas 1.5
*
* An old variant of the official <b>Source_flux_lambda</b> component.
*
* %Description
* Models a reactor source with a flat energy distribution and a
* given neutron flux.
* This is useful for simulations where the absolute value of neutron flux,
* detector counts, etc. is needed for comparison with real instruments and
* experiments.
*
* Example: Source_flat(radius=0.1, dist=2, xw=.1, yh=.1, E0=14, dE=2, flux=1e13)
*
* %Parameters
* INPUT PARAMETERS
*
* radius: (m) Radius of circle in (x,y,0) plane where neutrons
* are generated.
* dist: (m) Distance to target along z axis.
* xw: (m) Width(x) of target
* yh: (m) Height(y) of target
* E0: (meV) Mean energy of neutrons.
* dE: (meV) Energy spread of neutrons.
* flux: (n/s/cm^2/Angs) Neutron flux
*
* %End
***********************************************************************/
DEFINE COMPONENT Source_flux
DEFINITION PARAMETERS (radius=0.1, dist=2.33, xw=0.05, yh=0.1, E0=14.68, dE=2, flux=1e13)
SETTING PARAMETERS ()
OUTPUT PARAMETERS (p_in)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
DECLARE
%{
double hdiv,vdiv;
double p_in;
%}
INITIALIZE
%{
double factor, lambda_min, lambda_max, delta_lambda, source_area;
lambda_min = sqrt(81.81/(E0+dE)); /* AAngstroem */
lambda_max = sqrt(81.81/(E0-dE));
delta_lambda = lambda_max - lambda_min;
source_area = radius*radius*PI*1e4; /* cm^2 */
p_in = flux/mcget_ncount()*delta_lambda*source_area;
%}
TRACE
%{
double chi,v,r, E;
double xf, yf, dx, dy, rf;
z=0;
chi=2*PI*rand01(); /* Choose point on source */
r=sqrt(rand01())*radius; /* with uniform distribution. */
x=r*cos(chi);
y=r*sin(chi);
xf = 0.5*xw*randpm1(); /* Choose focusing position uniformly */
yf = 0.5*yh*randpm1();
dx = xf-x;
dy = yf-y;
rf = sqrt(dx*dx+dy*dy+dist*dist);
p = p_in*xw*yh*dist/(rf*rf*rf); /* target area * cos(phi)/rf^2 */
E=E0+dE*randpm1(); /* Assume linear distribution */
v=sqrt(E)*SE2V;
vz=v*dist/rf;
vy=v*dy/rf;
vx=v*dx/rf;
%}
MCDISPLAY
%{
magnify("xy");
circle("xy",0,0,0,radius);
%}
END
-------------- next part --------------
/*******************************************************************************
*
* McStas, the neutron ray-tracing package: Source_flux_lambda.comp
* Copyright 1997-2001 Risoe National Laboratory, Roskilde, Denmark
*
* %I
* Written by: Kristian Nielsen
* Date: 1998
* Version: $Revision: 1.6 $
* Origin: McStas 1.5
*
* Neutron source with flat wavelength spectrum and user-specified flux.
*
* %D
* The routine is a circular neutron source, which aims at a square target
* centered at the beam (in order to improve MC-acceptance rate). The angular
* divergence is then given by the dimensions of the target. The neutron
* wavelength is uniformly distributed between lambda_0 - d_lambda and
* lambda_0 + d_lambda. The source flux is specified in neutrons per steradian
* per square cm per AAngstroem.
*
* Example: Source_flux_lambda(radius=0.1, dist=2, xw=.1, yh=.1,
* lambda_0=2.36, d_lambda=.3, flux=1e13)
*
* %P
* radius: (m) Radius of circle in (x,y,0) plane where neutrons
* are generated.
* dist: (m) Distance to target along z axis.
* xw: (m) Width(x) of target
* yh: (m) Height(y) of target
* lambda_0: (AA) Mean wavelength of neutrons.
* d_lambda: (AA) Wavelength spread of neutrons.
* flux: (1/(cm**2*st*AA)) Source flux
*
* %E
*******************************************************************************/
DEFINE COMPONENT Source_flux_lambda
DEFINITION PARAMETERS ()
SETTING PARAMETERS (radius, dist, xw, yh, lambda_0, d_lambda, flux)
OUTPUT PARAMETERS (p_in)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
DECLARE
%{
double p_in;
%}
INITIALIZE
%{
double factor, delta_lambda, source_area;
delta_lambda = 2*d_lambda;
source_area = radius*radius*PI*1e4; /* cm^2 */
p_in = flux/mcget_ncount()*delta_lambda*source_area;
%}
TRACE
%{
double chi,lambda,v,r;
double xf, yf, dx, dy, rf;
z=0;
chi=2*PI*rand01(); /* Choose point on source */
r=sqrt(rand01())*radius; /* with uniform distribution. */
x=r*cos(chi);
y=r*sin(chi);
xf = 0.5*xw*randpm1(); /* Choose focusing position uniformly */
yf = 0.5*yh*randpm1();
dx = xf-x;
dy = yf-y;
rf = sqrt(dx*dx+dy*dy+dist*dist);
p = p_in*xw*yh*dist/(rf*rf*rf); /* target area * cos(phi)/rf^2 */
lambda = lambda_0+d_lambda*randpm1();
v = K2V*(2*PI/lambda);
vz=v*dist/rf;
vy=v*dy/rf;
vx=v*dx/rf;
%}
MCDISPLAY
%{
magnify("xy");
circle("xy",0,0,0,radius);
%}
END
-------------- next part --------------
/*******************************************************************************
*
* McStas, the neutron ray-tracing package: Source_gen.comp
* Copyright 1997-2001 Risoe National Laboratory, Roskilde, Denmark
*
* Component: Source_gen
*
* %I
* Written by: Kristian Nielsen, Kim Lefmann and Emmanuel Farhi
* Date: October 30, 1997
* Version: $Revision: 1.7 $
* Origin: McStas 1.5
* Modified by: EF, Aug 27, 2001 ; can use Energy/wavelength and I1
* Modified by: EF, Sep 18, 2001 ; corrected illumination bug. Add options
*
* Circular/squared neutron source with flat or Maxwellian energy spectrum
* (possibly spatially gaussian)
*
* %D
* This routine is a neutron source (rectangular or circular), which aims at
* a square target centered at the beam (in order to improve MC-acceptance
* rate). The angular divergence is then given by the dimensions of the
* target.
* The neutron energy/wavelength is distributed between Emin=E0-dE and
* Emax=E0+dE or Lmin=Lambda0-dLambda and Lmax=Lambda0+dLambda. The I1 may
* be either arbitrary (I1=0), or specified in neutrons per steradian per
* square cm per Angstrom. A Maxwellian spectra may be selected if you give
* the source temperatures (up to 3). The source shape is defined by its radius,
* or can alternatively be squared if you specify non-zero h and w parameters.
* The beam is spatially uniform, but becomes gaussian if one of the source
* dimensions (radius, h or w) is negative, or you set the gaussian flag.
* Divergence profiles are triangular. The source may have a thickness.
* For the Maxwellian spectra, the generated intensity is dPhi/dLambda (n/s/AA)
* For flat spectra, the generated intensity is Phi (n/s).
*
* Usage example:
* Source_gen(radius=0.1,Lambda0=2.36,dLambda=0.16, T1 = 20, T2 = 38, I1=1e13)
* Source_gen(h=0.1,w=0.1,Emin=1,Emax=3, I1=1e13, verbose=1, gaussian=1)
*
* Some sources parameters:
* PSI cold source T1=150.42 I1=3.67e11 T2=38.74 I2=3.64e11
* T3=14.84 I3=0.95e11
* ILL VCS cold source T1=216.8 I1=1.24e+13 T2=33.9 I2=1.02e+13
* (H1) T3=16.7 I3=3.0423e+12
* ILL HCS cold source T1=213.4 I1=1.47e+13 T2=83.1 I2=3.26e+12
* (H5) T3=26.5 I3=1.21e13
* ILL Thermal tube T1=683.7 I1=1.7278e+13 T2=257.7 I2=7.3823e+13
* (H12) T3=6.71 I3=6.7125e+09
*
*
* %P
* radius: (m) Radius of circle in (x,y,0) plane where neutrons
* are generated. You may also use 'h' and 'w' for a square source
* dist: (m) Distance to target along z axis.
* xw: (m) Width(x) of target. If dist is not given, xw = horz. div (deg)
* yh: (m) Height(y) of target. If dist is not given, yh = vert. div (deg)
* E0: (meV) Mean energy of neutrons. You may also use Lambda0.
* dE: (meV) Energy spread of neutrons (half width).
*
* Optional parameters:
* Lambda0: (AA) Mean wavelength of neutrons.
* dLambda: (AA) Wavelength spread of neutrons (half width).
* Emin: (meV) Minimum energy of neutrons
* Emax: (meV) Maximum energy of neutrons
* Lmin: (meV) Minimum wavelength of neutrons
* Lmax: (meV) Maximum wavelength of neutrons
* h (m) Source y-height (then does not use radius parameter)
* w (m) Source x-width (then does not use radius parameter)
* length (m) Source z-length (not anymore flat)
* T1 (K) Temperature of the Maxwellian source (0=none)
* I1: (1/(cm**2*st*AA)) Source flux per solid angle, area and Angstrom
* T2 (K) Second Maxwellian source Temperature (0=none)
* I2: (1/(cm**2*st*AA)) Second Maxwellian Source flux
* T3 (K) Third Maxwellian source Temperature (0=none)
* I3: (1/(cm**2*st*AA)) Third Maxwellian Source flux
* gaussian (0/1) Use gaussian beam (normal distributions)
* verbose (0/1) display info about the source. -1 unactivate source.
*
* %E
*******************************************************************************/
DEFINE COMPONENT Source_gen
DEFINITION PARAMETERS ()
SETTING PARAMETERS (radius=0.1, dist=57.296, xw=2, yh=2, E0=14.68, dE=0, Lambda0=2.36, dLambda=0.1627, I1=0, h=0, w=0, gaussian=0, verbose=0, T1=0,Lmin=0,Lmax=0,Emin=0,Emax=0,T2=0,I2=0,T3=0,I3=0,length=0)
OUTPUT PARAMETERS (p_in,lambda0, lambda02, L2P,lambda_0,d_lambda,lambda0b, lambda02b, L2Pb,lambda0c, lambda02c, L2Pc)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
DECLARE
%{
double p_in;
double lambda_0=0,d_lambda=0;
double lambda0, lambda02, L2P; /* first Maxwellian source */
double lambda0b, lambda02b, L2Pb; /* second Maxwellian source */
double lambda0c, lambda02c, L2Pc; /* second Maxwellian source */
%}
INITIALIZE
%{
double factor, delta_lambda, source_area, k;
if ((radius < 0) || (h < 0) || (w < 0))
gaussian = 1;
else
gaussian = 0;
if ((Emin != 0) && (Emax != 0))
{ E0 = (Emin+Emax)/2; dE=fabs(Emin-Emax)/2; }
if ((E0 != 0) && (dE != 0))
{
Lmin = sqrt(81.81/(E0+dE)); /* AAngstroem */
Lmax = sqrt(81.81/(E0-dE));
}
if ((Lmin != 0) && (Lmax != 0))
{ Lambda0 = (Lmin+Lmax)/2; dLambda=fabs(Lmin-Lmax)/2; }
if ((Lambda0 != 0) || (dLambda != 0))
{ lambda_0 = Lambda0; d_lambda=dLambda; }
radius = fabs(radius); w=fabs(w); h=fabs(h); I1=fabs(I1);
lambda_0=fabs(lambda_0); d_lambda=fabs(d_lambda);
xw = fabs(xw); yh=fabs(yh); dist=fabs(dist);
if ((lambda_0-d_lambda <= 0) || (lambda_0+d_lambda <= 0))
{
fprintf(stderr,"Source_gen: Warning: Wavelength will reach negative values\n");
}
if (dist == 0)
{
fprintf(stderr,"Source_gen: warning: focusing distance is null. Set to 1 rad\n");
dist = 57.296;
}
delta_lambda = 2*d_lambda;
Lmin = lambda_0 - d_lambda; /* AAngstroem */
Lmax = lambda_0 + d_lambda;
if (I1 != 0)
{
if ((h == 0) || (w == 0))
source_area = radius*radius*PI*1e4; /* circular cm^2 */
else
source_area = h*w*1e4; /* square cm^2 */
factor = I1*source_area*delta_lambda/mcget_ncount();
p_in = 4*factor; /* p_in = (4*hdiv*vdiv)*factor; */
}
else
p_in = 1/PI; /* Small angle approx. */
k = 1.38066e-23; /* k_B */
if (T1 > 0)
{
lambda0 = 1.0e10*sqrt(HBAR*HBAR*4.0*PI*PI/2.0/MNEUTRON/k/T1);
lambda02 = lambda0*lambda0;
L2P = 2*lambda02*lambda02;
}
else
{ lambda0 = lambda_0; }
if (T2 > 0)
{
lambda0b = 1.0e10*sqrt(HBAR*HBAR*4.0*PI*PI/2.0/MNEUTRON/k/T2);
lambda02b = lambda0b*lambda0b;
L2Pb = 2*lambda02b*lambda02b;
}
else
{ lambda0b = lambda_0; }
if (T3 > 0)
{
lambda0c = 1.0e10*sqrt(HBAR*HBAR*4.0*PI*PI/2.0/MNEUTRON/k/T3);
lambda02c = lambda0c*lambda0c;
L2Pc = 2*lambda02c*lambda02c;
if (I3 == 0) I3 = I1;
}
else
{ lambda0b = lambda_0; }
if (verbose == 1)
{
printf("Source_gen: component %s ", NAME_CURRENT_COMP);
if ((h == 0) || (w == 0))
printf("(square)");
else
printf("(disk)");
printf("\n spectra ");
printf("%.2f to %.2f AA (%.2f to %.2f meV)", Lmin, Lmax, 81.81/Lmax/Lmax, 81.81/Lmin/Lmin);
if (gaussian)
printf(", gaussian beam");
printf("\n");
if (T1 != 0)
printf(" T1=%.1f K (%.2f AA)", T1, lambda0);
if (T2 != 0)
printf(", T2=%.1f K (%.2f AA)", T2, lambda0b);
if (T3 != 0)
printf(", T3=%.1f K (%.2f AA)", T3, lambda0c);
if (T1) printf("\n");
}
else
if (verbose == -1)
printf("Source_gen: component %s unactivated", NAME_CURRENT_COMP);
%}
TRACE
%{
double theta0,phi0,theta1,phi1,chi,theta,phi,v,r, lambda;
double Maxwell, lambda2, lambda5;
if (verbose >= 0)
{
z=0;
if ((h == 0) || (w == 0))
{
chi=2*PI*rand01(); /* Choose point on source */
r=sqrt(rand01())*radius; /* with uniform distribution. */
x=r*cos(chi);
y=r*sin(chi);
}
else
{
x = w*randpm1()/2; /* select point on source (uniform) */
y = h*randpm1()/2;
}
if (length != 0)
z = length*randpm1()/2;
theta0= -atan((x-xw/2.0)/dist); /* Angles to aim at target */
phi0 = -atan((y-yh/2.0)/dist);
theta1= -atan((x+xw/2.0)/dist);
phi1 = -atan((y+yh/2.0)/dist);
/* shot towards target : flat distribution */
if (gaussian)
{
theta= theta0+(theta1- theta0)*(randnorm()*FWHM2RMS+0.5);
phi = phi0 +(phi1 - phi0) *(randnorm()*FWHM2RMS+0.5);
}
else
{
theta= theta0+(theta1- theta0)*rand01();
phi = phi0 +(phi1 - phi0) *rand01();
}
/* Assume linear distribution */
lambda = lambda_0+d_lambda*randpm1();
if (lambda <= 0) ABSORB;
lambda2 = lambda*lambda;
lambda5 = lambda2*lambda2*lambda;
v = K2V*(2*PI/lambda);
p *= p_in*fabs((theta1- theta0)*(phi1 - phi0));
if (T1 > 0)
{
Maxwell= L2P/lambda5*exp(-lambda02/lambda2); /* 1/AA */
if ((T2 > 0) && (I1 != 0))
{
if (I2 == 0) I2 = I1;
Maxwell += (I2/I1)*L2Pb/lambda5*exp(-lambda02b/lambda2);
}
if ((T3 > 0) && (I1 != 0))
{
if (I3 == 0) I3 = I1;
Maxwell += (I3/I1)*L2Pc/lambda5*exp(-lambda02c/lambda2);
}
p *= Maxwell;
}
vz=v*cos(phi)*cos(theta);
vy=v*sin(phi);
vx=v*cos(phi)*sin(theta);
SCATTER;
}
%}
MCDISPLAY
%{
double xmin;
double xmax;
double ymin;
double ymax;
if ((h == 0) || (w == 0))
{
magnify("xy");
circle("xy",0,0,0,radius);
if (gaussian)
circle("xy",0,0,0,radius/2);
}
else
{
xmin = -w/2; xmax = w/2;
ymin = -h/2; ymax = h/2;
magnify("xy");
multiline(5, (double)xmin, (double)ymin, 0.0,
(double)xmax, (double)ymin, 0.0,
(double)xmax, (double)ymax, 0.0,
(double)xmin, (double)ymax, 0.0,
(double)xmin, (double)ymin, 0.0);
if (gaussian)
circle("xy",0,0,0,sqrt(w*w+h*h)/4);
}
%}
END
More information about the mcstas-users
mailing list