[neutron-mc] Bender: how it works
Emmanuel Farhi
farhi at ill.fr
Tue Sep 17 17:12:02 CEST 2002
Hello Markus and neutron-mc guys !
If we look at the usual formulas for curved guide computation we can
write (for the example given)
m=2 (coating)
H=0.06 [m] (width)
Gamma_c = m*0.1 deg/Angs (critical angle)
r=1640.0 [m] (radius)
Then the critical wavelength is simply (from Mildner and Coppley papers)
Lambda_c = 180*sqrt(2*H/r)/2/pi/Gamma_c = 1.22 Angs for m=2 (2.45 for
m=1)
This means that the transmitted intensity versus the wavelength should
vary as Lambda^2 below Lambda_c, tends to be linear after, and then
saturates far away.
I've been running the 'test.instr' given by Markus with
mcstas-1.6-ill-alpha, and new developpment versions. They all give the
same result, and it seems reasonable (see mcstas.ps plot from
test2lam.dat with 1e6 neutrons)
I also attach the files for components Source_flat_lambda and Bender
that I used. There may be differences with your library. Please check
(specially the source).
Cheers, Emmanuel.
--
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 20 76 48 \__U_/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman2.mcstas.org/pipermail/mcstas-users/attachments/20020917/8c0adfac/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: mcstas.ps
Type: application/postscript
Size: 18639 bytes
Desc: not available
URL: <http://mailman2.mcstas.org/pipermail/mcstas-users/attachments/20020917/8c0adfac/attachment.ps>
-------------- next part --------------
/*******************************************************************************
*
* McStas, the neutron ray-tracing package: Bender component
* Copyright 2000-2001 Risoe National Laboratory, Roskilde, Denmark
* Component: Bender
*
* %Identification
* Written by: Philipp Bernhardt
* Date: Februar 7 1999
* Origin: McStas 1.5/Uni. Erlangen (Germany)
* Version: 1.0
*
* Models a curved neutron guide.
*
* %Description
* Models a curved neutron guide.The entrance lies in the X-Y plane, centered
* on the Z axis. The neutrons will also leave the bender in the X-Y plane at
* the z-value l=r*Win, so you don't have to calculate the real exit coordinates
* and you do not need a new arm. The bender is bent to the negative X axis;
* it behaves like a parallel guide in the Y axis.
* You may either enter the deviation angle 'Win' or the length 'l'.
* The bender is shown straight, even if the implemented model is curved.
* Thus, you don't have to shift the following component by the curved deviation.
*
* Example: Bender(w=0.05,h=0.12,r=250,d=0.001,Win=0.04,k=1,
* R0a=0.99,Qca=0.021,alphaa=6.07,ma=2,Wa=0.003,
* R0i=0.99,Qci=0.021,alphai=6.07,mi=2,Wi=0.003,
* R0s=0.99,Qcs=0.021,alphas=6.07,ms=2,Ws=0.003)
*
* See also <a href="http://neutron.risoe.dk/mcstas/support/misc/Bender.html">Additional note</a> from <a href="mailto:philipp.bernhardt at krist.uni-erlangen.de">Philipp Bernhardt</a>.
*
*
* %Parameters
* INPUT PARAMETERS:
*
* w: (m) Width at the bender entry and exit
* h: (m) Height at the bender entry and exit
* r: (m) Radius of the bender
* d: (m) Thickness of the partition, which separates the channels
* Win: (rad) Angle of the deflection of the whole bender
* k: (1) Number of channels inside the bender
* R0a: (1) Low-angle reflectivity at the bender <b>concave</b> side
* Qca: (AA-1) Critical scattering vector
* alphaa: (AA) Slope of reflectivity
* ma: (1) m-value of material
* Wa: (AA-1) Width of supermirror cut-off
* R0i: (1) Low-angle reflectivity at the bender <b>convex</b> side
* Qci: (AA-1) Critical scattering vector
* alphai: (AA) Slope of reflectivity
* mi: (1) m-value of material
* Wi: (AA-1) Width of supermirror cut-off
* R0s: (1) Low-angle reflectivity at bender <b>top and bottom</b> side
* Qcs: (AA-1) Critical scattering vector
* alphas: (AA) Slope of reflectivity
* ms: (1) m-value of material
* Ws: (AA-1) Width of supermirror cut-off
*
* Optional parameters:
* l: (m) length of bender l=r*Win
*
* %End
*******************************************************************************/
DEFINE COMPONENT Bender
DEFINITION PARAMETERS
(w,h,r,d=0.001,Win=0.04,k=1,R0a=0.99,Qca=0.021,alphaa=6.07,ma=2,Wa=0.003,R0i=0.99,Qci=0.021,alphai=6.07,mi=2,Wi=0.003,R0s=0.99,Qcs=0.021,alphas=6.07,ms=2,Ws=0.003,l=0)
SETTING PARAMETERS ()
OUTPUT PARAMETERS (bk, mWin)
STATE PARAMETERS (x, y, z, vx, vy, vz, t, s1, s2, p)
SHARE
%{
#ifndef BENDER_DECLARE
#define BENDER_DECLARE
double bd_sgn(double x) {
if (x>=0)
return 1.0;
else
return -1.0; }
#endif
%}
DECLARE
%{
double bk, mWin;
%}
INITIALIZE
%{
if (r <0)
{ fprintf(stderr,"Bender: error: %s: to bend in the other direction\n", NAME_CURRENT_COMP);
fprintf(stderr," rotate comp on z-axis by 180 deg.\n"); exit(-1); }
if (k*d > w)
{ fprintf(stderr,"Bender: error: %s has (k*d > w).\n", NAME_CURRENT_COMP);
exit(-1); }
if (w*h*r*Win*k == 0)
{ fprintf(stderr,"Bender: error: %s has one of w,h,r,Win,k null.\n", NAME_CURRENT_COMP);
exit(-1); }
/* width of one channel + thickness d of partition */
mWin = Win;
if (l!= 0 && r != 0) mWin = (double)l/(double)r;
bk=(w+d)/k;
%}
TRACE
%{
int i,num,numa,numi;
double dru,ab,dab,R,Q,arg,arga,argi,Ta,vpl;
double einmWin,ausmWin,zykmWin,aeumWin,innmWin,ref,innref,aeuref;
double einzei,auszei,zykzei;
/* does the neutron hit the bender at the entrance? */
PROP_Z0;
if ((fabs(x)<w/2) && (fabs(y)<h/2))
{
/*** reflections in the XZ-plane ***/
/* distance between neutron and concave side of the channel at the entrance */
dru=floor((w/2-x)/bk)*bk;
ab=w/2.0-x-dru;
/* radius of the channel */
R=r-dru;
/* does the neutron hit the partition at the entrance? */
if (ab<bk-d)
{
/* velocity in the XZ-plane */
vpl=sqrt(vx*vx+vz*vz);
/* divergence of the neutron at the entrance */
einmWin=atan(vx/vz);
/* maximal distance between neutron and concave side of the channel */
dab=R-cos(einmWin)*(R-ab);
/* reflection angle at the concave side */
aeumWin=acos((R-dab)/R);
/* reflection coefficient at the concave side */
arga=0.0;
Q=2.0*V2K*vpl*sin(aeumWin);
if (Q<=Qca)
aeuref=R0a;
else {
arga=(Q-ma*Qca)/Wa;
aeuref=0.5*R0a*(1.0-tanh(arga))*(1.0-alphaa*(Q-Qca)); }
/* does the neutron hit the convex side of the channel? */
innmWin=0.0;
innref=1.0;
argi=0.0;
if (dab>bk-d)
{
/* reflection coefficient at the convex side */
innmWin=acos((R-dab)/(R-bk+d));
Q=2.0*V2K*vpl*sin(innmWin);
if (Q<=Qci)
innref=R0i;
else {
argi=(Q-mi*Qci)/Wi;
innref=0.5*R0i*(1.0-tanh(argi))*(1.0-alphai*(Q-Qci)); }
}
/* divergence of the neutron at the exit */
zykmWin=2.0*(aeumWin-innmWin);
ausmWin=fmod(mWin+einmWin+aeumWin-innmWin
*(1.0+bd_sgn(einmWin)),zykmWin)-zykmWin/2.0;
ausmWin+=innmWin*bd_sgn(ausmWin);
/* number of reflections at the concave side */
numa=(mWin+einmWin+aeumWin-innmWin*(1.0+bd_sgn(einmWin)))/zykmWin;
/* number of reflections at the convex side */
numi=numa;
if (ausmWin*einmWin<0)
{
if (ausmWin-einmWin>0)
numi++;
else
numi--;
}
/* is the reflection coefficient too small? */
if (((numa>0) && (arga>10.0)) || ((numi>0) && (argi>10.0)))
ABSORB;
/* calculation of the neutron probability weight p */
for (i=1;i<=numa;i++)
p*=aeuref;
for (i=1;i<=numi;i++)
p*=innref;
/* time to cross the bender */
Ta=(2*numa*(tan(aeumWin)-tan(innmWin))
+tan(ausmWin)-tan(einmWin)
-tan(innmWin)*(bd_sgn(ausmWin)-bd_sgn(einmWin)))
*(R-dab)/vpl;
t+=Ta;
/* distance between neutron and concave side of channel at the exit */
ab=R-(R-dab)/cos(ausmWin);
/* calculation of the exit coordinates in the XZ-plane */
x=w/2.0-ab-dru;
z=r*mWin;
vx=sin(ausmWin)*vpl;
vz=cos(ausmWin)*vpl;
/*** reflections at top and bottom side (Y axis) ***/
if (vy!=0.0)
{
/* reflection coefficent at the top and bottom side */
Q=2.0*V2K*fabs(vy);
if (Q<=Qcs)
ref=R0s;
else {
arg=(Q-ms*Qcs)/Ws;
ref=0.5*R0s*(1.0-tanh(arg))*(1.0-alphas*(Q-Qcs)); }
/* number of reflections at top and bottom */
einzei=h/2.0/fabs(vy)+y/vy;
zykzei=h/fabs(vy);
num=(Ta+einzei)/zykzei;
/* time between the last reflection and the exit */
auszei=fmod(Ta+einzei,zykzei);
/* is the reflection coefficient too small? */
if ((num>0) && (arg>10.0))
ABSORB;
/* calculation of the probability weight p */
for (i=1;i<=num;i++) {
p*=ref;
vy*=-1.0; }
/* calculation of the exit coordinate */
y=auszei*vy-vy*h/fabs(vy)/2.0;
} /* if (vy!=0.0) */
SCATTER;
} /* if (dab>bk-d) */
else
ABSORB; /* hit separating walls */
}
else /* if ((fabs(x)<w/2) && (fabs(y)<h/2)) */
ABSORB; /* miss entry window */
%}
MCDISPLAY
%{
int i;
double w1c, w2c, h1, h2, L, w1, w2;
w1c = (w + d)/(double)k;
w2c = w1c; h1 = h; h2 = h;
L = r*mWin; w1 = w; w2 = w;
magnify("xy");
for(i = 0; i < k; i++)
{
multiline(5,
i*w1c - w1/2.0, -h1/2.0, 0.0,
i*w2c - w2/2.0, -h2/2.0, (double)L,
i*w2c - w2/2.0, h2/2.0, (double)L,
i*w1c - w1/2.0, h1/2.0, 0.0,
i*w1c - w1/2.0, -h1/2.0, 0.0);
multiline(5,
(i+1)*w1c - d - w1/2.0, -h1/2.0, 0.0,
(i+1)*w2c - d - w2/2.0, -h2/2.0, (double)L,
(i+1)*w2c - d - w2/2.0, h2/2.0, (double)L,
(i+1)*w1c - d - w1/2.0, h1/2.0, 0.0,
(i+1)*w1c - d - w1/2.0, -h1/2.0, 0.0);
}
line(-w1/2.0, -h1/2.0, 0.0, w1/2.0, -h1/2.0, 0.0);
line(-w2/2.0, -h2/2.0, (double)L, w2/2.0, -h2/2.0, (double)L);
%}
END
-------------- next part --------------
/*******************************************************************************
*
* McStas, the neutron ray-tracing package: Source_flat_lambda.comp
* Copyright 1997-2001 Risoe National Laboratory, Roskilde, Denmark
*
* Component: Source_flat_lambda
*
* %I
* Written by: Kristian Nielsen and Kim Lefmann
* Date: 1999
* Version: $Revision: 1.2 $
* Origin: McStas release 1.5.1
* Modified by: KN, 1999 from Source_flat.comp
* Modified by: KL, October 8, 2001
* Modified by: Emmanuel Farhi, October 30, 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 --------------
# Instrument-source: test.instr
# Date: Tue Sep 17 16:49:46 2002
# Ncount: 1e+06/1e+06
# Trace: no
# type: array_1d(100)
# component: lmonben1
# title: Wavelength monitor
# filename: 'test2lam.dat'
# variables: L I I_err N
# xvar: L
# yvar: (I,I_err)
# xlabel: 'Wavelength [AA]'
# ylabel: 'Intensity'
# xlimits: 1 13
# statistics: Min_y=3.39545e-09, Max_y=3.92813e-07, Mean_y= 1.5587e-07
# statistics: Sum_y=1.5587e-05, X0=9.59903, dX=2.54913
1.06 3.39545e-09 8.48919e-10 49
1.18 4.37996e-09 9.6054e-10 54
1.3 4.84361e-09 9.86483e-10 73
1.42 6.1298e-09 1.12775e-09 84
1.54 6.59692e-09 1.16944e-09 110
1.66 6.10644e-09 1.10794e-09 107
1.78 8.44359e-09 1.29007e-09 135
1.9 9.57924e-09 1.38572e-09 170
2.02 8.92885e-09 1.32436e-09 195
2.14 1.35994e-08 1.5967e-09 252
2.26 1.27096e-08 1.58203e-09 254
2.38 1.35136e-08 1.60806e-09 285
2.5 1.55953e-08 1.73185e-09 315
2.62 1.9543e-08 1.93376e-09 340
2.74 2.23384e-08 2.07486e-09 375
2.86 2.0597e-08 2.00105e-09 394
2.98 2.27889e-08 2.06995e-09 427
3.1 2.59097e-08 2.21029e-09 485
3.22 2.83344e-08 2.32127e-09 527
3.34 3.53963e-08 2.60966e-09 638
3.46 3.54004e-08 2.59811e-09 654
3.58 3.61231e-08 2.62728e-09 647
3.7 4.447e-08 2.90773e-09 764
3.82 3.98292e-08 2.76219e-09 757
3.94 3.73348e-08 2.66667e-09 772
4.06 4.73015e-08 3.01006e-09 909
4.18 4.83861e-08 3.03497e-09 908
4.3 5.72317e-08 3.31391e-09 963
4.42 5.77996e-08 3.32194e-09 1070
4.54 6.1714e-08 3.42511e-09 1116
4.66 6.25438e-08 3.46906e-09 1147
4.78 7.32873e-08 3.72548e-09 1267
4.9 7.11333e-08 3.66417e-09 1328
5.02 6.82532e-08 3.5747e-09 1344
5.14 7.63744e-08 3.8211e-09 1450
5.26 8.45252e-08 3.98365e-09 1520
5.38 8.37023e-08 3.96293e-09 1559
5.5 8.61337e-08 4.00925e-09 1648
5.62 8.41357e-08 3.96257e-09 1651
5.74 9.18248e-08 4.16195e-09 1772
5.86 9.82229e-08 4.29541e-09 1852
5.98 1.00559e-07 4.37385e-09 1931
6.1 1.06835e-07 4.45238e-09 1962
6.22 1.12234e-07 4.58235e-09 2085
6.34 1.21032e-07 4.77881e-09 2135
6.46 1.13003e-07 4.60067e-09 2124
6.58 1.25013e-07 4.83136e-09 2336
6.7 1.23905e-07 4.80728e-09 2379
6.82 1.24167e-07 4.81742e-09 2345
6.94 1.30358e-07 4.93352e-09 2465
7.06 1.29637e-07 4.88723e-09 2564
7.18 1.41628e-07 5.1429e-09 2646
7.3 1.43255e-07 5.18424e-09 2688
7.42 1.40465e-07 5.11929e-09 2643
7.54 1.5506e-07 5.37111e-09 2861
7.66 1.5638e-07 5.38142e-09 2900
7.78 1.67634e-07 5.5465e-09 3052
7.9 1.62123e-07 5.46052e-09 3031
8.02 1.65955e-07 5.54e-09 3099
8.14 1.81474e-07 5.82729e-09 3212
8.26 1.85497e-07 5.86463e-09 3249
8.38 1.87456e-07 5.86351e-09 3319
8.5 1.88236e-07 5.89105e-09 3355
8.62 1.92801e-07 5.95464e-09 3359
8.74 2.01544e-07 6.10109e-09 3412
8.86 2.04073e-07 6.14177e-09 3495
8.98 2.15846e-07 6.29408e-09 3542
9.1 2.20208e-07 6.38566e-09 3507
9.22 2.23723e-07 6.42285e-09 3673
9.34 2.34697e-07 6.58221e-09 3726
9.46 2.26025e-07 6.46345e-09 3580
9.58 2.34388e-07 6.55377e-09 3720
9.7 2.46912e-07 6.74562e-09 3858
9.82 2.5763e-07 6.88126e-09 3918
9.94 2.68701e-07 7.01331e-09 3920
10.06 2.57266e-07 6.87365e-09 3864
10.18 2.70653e-07 7.07253e-09 3832
10.3 2.7881e-07 7.14479e-09 3928
10.42 2.815e-07 7.19768e-09 4023
10.54 2.65846e-07 6.96374e-09 3947
10.66 2.87487e-07 7.23674e-09 3967
10.78 2.99664e-07 7.38873e-09 4029
10.9 2.93714e-07 7.31519e-09 3980
11.02 3.01298e-07 7.41165e-09 3984
11.14 3.0429e-07 7.48426e-09 4061
11.26 3.07364e-07 7.47734e-09 4040
11.38 3.13543e-07 7.6085e-09 3944
11.5 3.18987e-07 7.62315e-09 4031
11.62 3.2914e-07 7.7488e-09 4130
11.74 3.26451e-07 7.7207e-09 4047
11.86 3.33047e-07 7.80729e-09 4087
11.98 3.52774e-07 8.02349e-09 4111
12.1 3.27641e-07 7.69304e-09 3963
12.22 3.49116e-07 8.00001e-09 4034
12.34 3.6878e-07 8.21574e-09 4112
12.46 3.75126e-07 8.30722e-09 4120
12.58 3.75791e-07 8.29981e-09 4148
12.7 3.77995e-07 8.30043e-09 4111
12.82 3.74977e-07 8.28818e-09 4074
12.94 3.92813e-07 8.48698e-09 4115
More information about the mcstas-users
mailing list