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
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
|
/**************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright 1997-2006, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: Pol_Bfield
*
* %I
* Written by: Erik B Knudsen, Peter Christiansen and Peter Willendrup
* Date: July 2011
* Origin: RISOE
*
* Magnetic field component.
*
* %D
*
* Region with a definable magnetic field.
*
* The component is nestable if the concentric flag is set (the default). This means that it may have a
* construction like the following:
* // START MAGNETIC FIELD
* COMPONENT msf =
* Pol_Bfield(xwidth=0.08, yheight=0.08, zdepth=0.2, Bx=0, By=-0.678332e-4, Bz=0)
* AT (0, 0, 0) RELATIVE armMSF
*
* // OTHER COMPONENTS INSIDE THE MAGNETIC FIELD CAN BE PLACED HERE
*
* // STOP MAGNETIC FIELD
* COMPONENT msf_stop = Pol_simpleBfield_stop(magnet_comp_stop=msf)
* AT (0, 0, 0) RELATIVE armMSF
*
* Note that if you have objects within the magnetic field that extend outside it you may get
* wrong results, because propagation within the field will then possibly extend outside, e.g.
* when using a tabled field. The evaluated field will then use the nearest defined field point
* _also_ outside the defintion area. If these outside-points have a non-zero field precession will
* continue - even after the neutron has left the field region.
*
* In between the two component instances the propagation routine
* PROP_DT also handles spin propagation.
* The current algorithm used for spin propagation is:
* SimpleNumMagnetPrecession
* in pol-lib.
*
* Example: Pol_Bfield(xwidth=0.1, yheight=0.1, zdepth=0.2, Bx=0, By=1, Bz=0)
* Pol_Bfield(xwidth=0.1, yheight=0.1, zdepth=0.2, field_type=1)
*
* Functions supplied by the system are (the number is the ID of the function to be given as the field_type parameter:
* 1. Constant magnetic field: Constant field (Bx,By,Bz) within the region
* 2. Rotating magnetic_field: Field is initially (0,By,0) but after a length of zdepth
* has rotated to (By,0,0)
* 3. Majorana magnetic_field: Field is initially (Bx,By,0) with By<<Bx, then linearly transforms to
* (-Bx,By,0) at z=zdepth.
* 4. MSF field:
* 5. RF field: A radio frequency field is modeled by an implcit use of a roating frame.
* 6. Gradient field: Similar to Majorana by without an x-component to the field. I.e. it (0,By,0) at z=0, and
* becomes (0,-By,0) AT z=zdepth.
*
* If the concentric parameter is set to 1 the magnetic field is turned on by an
* instance of this component, and turned off by an instance of Pol_simpleBfield_stop.
* Anything in between is considered inside the field.
* If concentric is zero the field is considered a closed component - neutrons are propagated
* through the field region, and no other components are permitted inside the field.
*
* %P
* INPUT PARAMETERS:
*
* xwidth: [m] Width of opening.
* yheight: [m] Height of opening.
* zdepth: [m] Length of field.
* radius: [m] Radius of field if it is cylindrical or spherical.
* Bx: [T] Parameter used for x composant of field.
* By: [T] Parameter used for y composant of field.
* Bz: [T] Parameter used for z composant of field.
* field_type: [] ID of function describing the magnetic field. 1:constant field, 2: rotating fiel, 3: majorana type field, 4: MSF, 5: RF-field, 6: gradient field.
* concentric: [ ] Allow components and/or other fields inside field. Requires a subsequent Pol_simpleBfield_stop component.
*
* %E
****************************************************************************/
DEFINE COMPONENT Pol_Bfield
SETTING PARAMETERS (int field_type=0, xwidth=0, yheight=0,zdepth=0, radius=0, Bx=0, By=0, Bz=0, concentric=1)
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
SHARE
%{
%include "pol-lib"
double fmax(double, double);
double fmin(double, double);
%}
DECLARE
%{
/* Larmor frequency and scalar product threshold*/
double Bprms[8];
int shape;
%}
INITIALIZE
%{
enum shapes {NONE=0, BOX, WINDOW, CYLINDER, SPHERE, ANY};
enum field_functions ff=field_type;
/* these default field functions are part of pol-lib */
if (ff==constant){
double *t=Bprms;
t[0]=Bx;t[1]=By;t[2]=Bz;
} else if (ff==rotating){
double *t=Bprms;
t[0]=By;t[1]=zdepth;
} else if (ff==majorana){
double *t=Bprms;
t[0]=Bx; t[1]=By; t[2]=zdepth;
}
if(xwidth && yheight && zdepth){
shape=BOX;
}else if (xwidth && yheight && !zdepth){
shape=WINDOW;
}else if(radius && yheight){
shape=CYLINDER;
}else if (radius) {
shape=SPHERE;
}else{
shape=NONE;
}
if(shape == WINDOW && !concentric){
printf("Warning (%s): Magnetic field has window shape and concentric==0 => Field volume == 0.\n",NAME_CURRENT_COMP);
}
if(shape == NONE) {
printf("Warning (%s): Magnetic field has no geometry. Full Z=0-plane is used as boundary.\n",NAME_CURRENT_COMP);
}
%}
TRACE
%{
double t0,t1;
int hit;
enum shapes {NONE=0, BOX, WINDOW, CYLINDER, SPHERE, ANY};
/*enter through whatever object we are*/
switch (shape){
case BOX:
hit=box_intersect(&t0,&t1,x,y,z,vx,vy,vz,xwidth,yheight,zdepth);
/*terminate neutrons which miss the component*/
if(!hit) ABSORB;
/*If we do hit - propagate to the start of the field unless the neutron is already there.*/
if(t0>FLT_EPSILON) {
PROP_DT(t0);
t1-=t0;
}else if (t0<-FLT_EPSILON && concentric){
PROP_DT(t1);
}
break;
case CYLINDER:
hit=cylinder_intersect(&t0,&t1,x,y,z,vx,vy,vz,radius,yheight);
/*terminate neutrons which miss the component*/
if(!hit)ABSORB;
/*If we do hit - propagate to the start of the field unless the neutron is already there.*/
if(t0>FLT_EPSILON) {
PROP_DT(t0);
t1-=t0;
}else if (t0<-FLT_EPSILON && concentric){
PROP_DT(t1);
}
break;
case WINDOW:
PROP_Z0;
if (2*x>xwidth || 2*x<-xwidth || 2*y>yheight || 2*y<-yheight){
/*terminate neutrons which miss the component*/
ABSORB;
}
break;
default:
PROP_Z0;
}
mcmagnet_push(_particle, field_type,&(ROT_A_CURRENT_COMP),&(POS_A_CURRENT_COMP),0,Bprms);
#ifdef MCDEBUG
mcmagnet_print_stack();
#endif
/*does the field "close/stop" itself*/
if (!concentric){
switch (shape){
case BOX:
PROP_DT(t1);
break;
case CYLINDER:
PROP_DT(t1);
break;
case WINDOW:
PROP_Z0;
/*terminate neutrons which miss the exit window*/
if (2*x>xwidth || 2*x<-xwidth || 2*y>yheight || 2*y<-yheight){
ABSORB;
}
break;
default:
PROP_Z0;
}
mcmagnet_pop(_particle);
}
%}
MCDISPLAY
%{
enum shapes {NONE=0, BOX, WINDOW, CYLINDER, SPHERE, ANY};
switch (shape){
case BOX:
box(0,0,0,xwidth,yheight,zdepth,0, 0, 1, 0);
break;
case CYLINDER:
circle("xz",0, yheight/2.0,0,radius);
circle("xz",0,-yheight/2.0,0,radius);
line(-radius,-yheight/2.0,0,-radius,yheight/2.0,0);
line( radius,-yheight/2.0,0, radius,yheight/2.0,0);
line(0,-yheight/2.0,-radius,0,yheight/2.0,-radius);
line(0,-yheight/2.0, radius,0,yheight/2.0, radius);
break;
case SPHERE:
circle("xz",0,0,0,radius);
circle("xy",0,0,0,radius);
circle("yz",0,0,0,radius);
break;
case WINDOW:
rectangle("xy",0,0,0,xwidth,yheight);
break;
}
%}
END
|
