1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
|
/***************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright 1997-2002, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: He3_cell
*
* %Identification
* Written by: Trefor Roberts & Erik B Knudsen
* Date: March 1999
* Origin: ILL/DTU Physics
* version: $Revision$
*
* Polarised 3He cell
*
* %Description
* Simple cylindrical polarised 3He neutron spin filter cell.
* The glass container for the cell is not included in the model.
*
* This component has been validated against:
* Batz, M, Baessler, S, Heil, W, et al., J Res Natl Inst Stand Technol. 2005;110(3):293–298.
*
* Example: He3_cell(radius=0.1,length=.2,pressure=3,p3he=0.7,bx=0,by=1e-3,bz=0)
*
* %Parameters
* Input parameters:
*
* radius: [m] radius of the cylinder
* length: [m] length of the cylinder
* pressure: [bar] pressure of the gas in the cell
* p3he: polarisation of the 3He gas [-1 to +1]
* bx: [tesla] x component of field at the cell
* by: [tesla] y component of field at the cell
* bz: [tesla] z component of field at the cell
*
* %End
***************************************************************************/
DEFINE COMPONENT He3_cell
SETTING PARAMETERS (radius,length,pressure=3,p3he=0.7,bx=0,by=1e-3,bz=0)
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
DECLARE
%{
%}
INITIALIZE
%{
%}
TRACE
%{
double t0,t1; /* time that neutron enters and leaves gas (s) */
double v,lambda; /* neutron velocity and wavelength (ms-1, Anstroms) */
double l_full; /* path length of neutron through gas (m) */
double dt0; /* time neutron spends in the gas (s) */
double opacity; /* opacity of the gas for this neutron (dimless) */
double omega; /* angle through which polarisation precesses
over the path through the cell (radians) */
double px,py,pz; /* auxilliary vector polarisation (-1 to +1) */
double bnx,bny,bnz; /* normal vector parallel to the magnetic field*/
const double gyro = 1.832e8; /*absolute value of the gyromagnetic ratio of the neutron (1/(s.T))*/
const double opac_cnv = 7.33; /*opacity conversion factor to express opacity in (bar.m.angstrom)*/
/* calculate the intersection times with the volume of gas, if the neutron
goes through the cell, continue with calculation otherwise all done.
Note that y and z are swapped - this is because the cylindrical axis of
a 3He cell lies along the beam. */
if(cylinder_intersect(&t0,&t1,x,z,y,vx,vz,vy,radius,length))
{
/* Calculate the neutron velocity and wavelength */
v=sqrt(vx*vx+vy*vy+vz*vz);
lambda=2*PI/(V2K*v);
/* Calculate the path length of the neutron through the gas */
dt0=t1-t0;
l_full=v*dt0;
/* Calculate the opacity of the cell for the path length travelled */
opacity=pressure*l_full*lambda*opac_cnv;
/* propagate the polarisation accross the cell (assuming a constant
magnetic field). The actual interaction point is not taken into account
as only the parallel and perpendicular components are important. */
omega=dt0*gyro*sqrt(bx*bx+by*by+bz*bz);
rotate(px,py,pz,sx,sy,sz,omega,bx,by,bz);
sx=px;
sy=py;
sz=pz;
/* adjust the neutron weight according to spin state relative to the
3He nuclei - antiparallel spins are as good as absorbed, whereas parallel
spins are transmitted (depending upon degree of polarisation of gas */
bnx=bx;bny=by;bnz=bz;
NORM(bnx,bny,bnz);
double pm1=scalar_prod(sx,sy,sz,bnx,bny,bnz);
double phiu,phid,T, Tup,Tdown;
phiu=((pm1)+1)/2.0;
phid=(1-(pm1))/2.0;
Tup=exp(-opacity*(1.0-p3he));
Tdown=exp(-opacity*(1.0+p3he));
p*=phiu*Tup + phid*Tdown;
/*set the outgoing polarization vector without touching the part perpendicular
to the polarisation of the gas .
Perpendicular comp: P_=S - P// = S - (S.b_n)b_n
Parallel comp: P// = Pn bn, where Pn is the polarisation efficiency.*/
px=sx - pm1*bnx;
py=sy - pm1*bny;
pz=sz - pm1*bnz;
double Pn=(Tup-Tdown)/(Tup+Tdown);
sx=px + Pn*bnx;
sy=py + Pn*bny;
sz=pz + Pn*bnz;
SCATTER;
}
%}
MCDISPLAY
%{
circle("xy",0.0,0.0,-length/2.0,radius);
circle("xy",0.0,0.0,length/2.0,radius);
line(0.0,radius,-length/2.0,0.0,radius,length/2.0);
line(radius,0.0,-length/2.0,radius,0.0,length/2.0);
line(0.0,-radius,-length/2.0,0.0,-radius,length/2.0);
line(-radius,0.0,-length/2.0,-radius,0.0,length/2.0);
%}
END
|
