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
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
|
/**************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright 1997-2006, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: Pol_triafield
*
* %I
* Written by: Morten Sales, based on Pol_constBfield by Peter Christiansen
* Date: 2013
* Origin: Helmholtz-Zentrum Berlin
*
* Constant magnetic field in a isosceles triangular coil
*
* %D
*
* Rectangular box with constant B field along y-axis (up) in a isosceles triangle.
* There is a guide (or precession) field as well. It is along y in the entire rectangular box.
* A neutron hitting outside the box opening or the box sides is absorbed.
*
*
* __________________
* | /\ |
* | Bguide/ \Bguide | x
* | / \ | ^
* | / \ | |
* | / B \ | |-----> z
* | / and \ |
* | / Bguide \ |
* | / \ |
* |/________________\|
*
* The angle of the inclination of the triangular field boundary is given by the arctangent to xwidth/(0.5*zdepth)
*
* This component does NOT take gravity into account.
*
* Example: Pol_triafield(xwidth=0.1, yheight=0.1, zdepth=0.2, B=1e-3, Bguide=0.0)
*
* %P
* INPUT PARAMETERS:
*
* xwidth: [m] Width of opening
* yheight: [m] Height of opening
* zdepth: [m] zdepth of field
* B: [T] Magnetic field along y-direction inside triangle
* Bguide: [T] Magnetic field along y-direction inside entire box
*
* CALCULATED PARAMETERS:
*
* %E
*******************************************************************************/
DEFINE COMPONENT Pol_triafield
SETTING PARAMETERS (xwidth, yheight, zdepth, B=0, Bguide=0)
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
SHARE
%{
double IntersectWall(double pos, double vel, double wallpos) {
/* Function to calculate where the neutron hit the wall */
if(vel==0)
return -1;
if(vel>0)
return (wallpos-pos)/vel;
else
return (-wallpos-pos)/vel;
}
%}
DECLARE
%{
/* Larmor frequency */
double omegaL;
double omegaLguide;
%}
INITIALIZE
%{
omegaL = 0;
omegaLguide = 0;
double velocity = 0, time = 0;
if ((xwidth<=0) || (yheight<=0) || (zdepth<=0)) {
fprintf(stderr, "Pol_filter: %s: Null or negative volume!\n"
"ERROR (xwidth, yheight, zdepth). Exiting\n",
NAME_CURRENT_COMP);
exit(1);
}
omegaL = -1.832472e8 * (B - Bguide); // B and Bguide is in Tesla
omegaLguide = -1.832472e8 * Bguide; // Bguide is in Tesla
%}
TRACE
%{
double deltaT, deltaTx, deltaTy, sx_in1, sz_in1, sx_in2, sz_in2, iz1, iz2, denom1, denom2, deltaTtria;
PROP_Z0;
if (!inside_rectangle(x, y, xwidth, yheight))
ABSORB;
// Time spent in Bguide-field
deltaT = zdepth/vz;
// This calculates the intersections on the xz-plane between the neutron trajectory and the triangular field boundaries
// The neutron trajectory is given by the points ( x, 0, 0) and ( x+vx, 0, vz)
// The first field boundary is given by the points (-xwidth/2, 0, 0) and ( xwidth/2, 0, zdepth/2)
// The second field boundary is given by the points ( xwidth/2, 0,zdepth/2) and (-xwidth/2, 0, zdepth)
// iz1 and iz2 are the z-values for the intersection
denom1 = (-vz)*((-xwidth/2)-xwidth/2)-(x-(x+vx))*(-zdepth/2);
iz1 = ((-x*vz)*(-zdepth/2)-(-vz)*(-(-xwidth/2)*zdepth/2))/denom1;
denom2 = (-vz)*(xwidth/2-(-xwidth/2))-(x-(x+vx))*(zdepth/2-zdepth);
iz2 = ((-x*vz)*(zdepth/2-zdepth)-(-vz)*(zdepth/2*(-xwidth/2)-xwidth/2*zdepth))/denom2;
// Time spent in triangular B-field
deltaTtria = (iz2-iz1)/vz;
// check that track goes throgh without hitting the walls
if (!inside_rectangle(x+vx*deltaT, y+vy*deltaT, xwidth, yheight)) {
// Propagate to the wall and absorb
deltaTx = IntersectWall(x, vx, xwidth/2);
deltaTy = IntersectWall(y, vy, yheight/2);
if (deltaTx>=0 && deltaTx<deltaTy)
deltaT = deltaTx;
else
deltaT = deltaTy;
PROP_DT(deltaT);
ABSORB;
}
PROP_DT(deltaT);
// These are the incoming spin directions
sx_in1 = sx;
sz_in1 = sz;
// This calculates the spin rotation caused by the guide/precession field
sz_in2 = cos(omegaLguide*deltaT)*sz_in1 - sin(omegaLguide*deltaT)*sx_in1;
sx_in2 = sin(omegaLguide*deltaT)*sz_in1 + cos(omegaLguide*deltaT)*sx_in1;
// This calculated the spin rotation caused by the triangular field
sz = cos(omegaL*deltaTtria)*sz_in2 - sin(omegaL*deltaTtria)*sx_in2;
sx = sin(omegaL*deltaTtria)*sz_in2 + cos(omegaL*deltaTtria)*sx_in2;
%}
MCDISPLAY
%{
box(0, 0, 0, xwidth, yheight, zdepth,0, 0, 1, 0);
%}
END
|
