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
|
/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright 1997-2002, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: Selector
*
* %Identification
* Written by: Peter Link, <a href="mailto:Andreas.Ostermann@frm2.tum.de">Andreas Ostermann</a>
* Date: MARCH 1999
* Origin: Uni. Gottingen (Germany)
*
* velocity selector (helical lamella type) such as <b>V_selector</b> component
*
* %Description
* Velocity selector consisting of rotating Soller-like blades
* defining a helically twisted passage.
* Geometry is defined by two identical apertures at 12 o'clock position,
* The origin is at the ENTRANCE of the selector.
*
* Example: Selector(xmin=-0.015, xmax=0.015, ymin=-0.025, ymax=0.025, length=0.25,
* nslit=72,d=0.0004, radius=0.12, alpha=48.298, nu=500)
* These are values for the D11@ILL Dornier 'Dolores' Velocity Selector (NVS 023)
*
* %VALIDATION
* Jun 2005: extensive external test, one minor problem
* Jan 2006: problem solved (for McStas-1.9.1)
* Validated by: K. Lieutenant
*
* %BUGS
* for transmission calculation, each neutron is supposed to be in the guide centre
*
* %Parameters
* INPUT PARAMETERS:
*
* xmin: [m] Lower x bound of entry aperture
* xmax: [m] Upper x bound of entry aperture
* ymin: [m] Lower y bound of entry aperture
* ymax: [m] Upper y bound of entry aperture
* xwidth: [m] Width of entry. Overrides xmin,xmax.
* yheight: [m] Height of entry. Overrides ymin,ymax.
* length: [m] rotor length
* nslit: [1] number of absorbing subdivinding spokes/lamella
* d: [m] width of spokes at beam-center
* radius: [m] radius of beam-center
* alpha: [deg] angle of torsion
* nu: [Hz] frequency of rotation, which is ideally 3956*alpha*DEG2RAD/2/PI/lambda/length
*
* %Links
* See also Additional notes <a href="http://mcstas.risoe.dk/pipermail/neutron-mc/1999q1/000134.html">March 1999</a> and <a href="http://mcstas.risoe.dk/pipermail/neutron-mc/1999q2/000136.html">Jan 2000</a>.
*
* %End
*******************************************************************************/
DEFINE COMPONENT Selector
SETTING PARAMETERS (xmin=-0.015, xmax=0.015, ymin=-0.025, ymax=0.025, length=0.25,
xwidth=0, yheight=0, nslit=72, d=0.0004, radius=0.12, alpha=48.298, nu=500)
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
INITIALIZE
%{
if (xwidth > 0) { xmax=xwidth/2; xmin=-xmax; }
if (yheight > 0) { ymax=yheight/2; ymin=-ymax; }
%}
TRACE
%{
double E;
double dt;
double open_angle, closed_angle, n_angle_in, n_angle_out;
double sel_phase, act_radius; // distance between neutron and selector axle
PROP_Z0;
E=VS2E*(vx*vx+vy*vy+vz*vz);
if (x<xmin || x>xmax || y<ymin || y>ymax)
ABSORB; /* because outside frame */
dt = length/vz;
/* get phase angle of selector rotor as MonteCarlo choice
only the free space between two neighboring spokes is taken
p is adjusted to transmission for parallel beam
*/
n_angle_in = atan2( x,y+radius)*RAD2DEG;
act_radius = sqrt(x*x + pow(radius+y,2));
closed_angle = d/act_radius*RAD2DEG;
open_angle = 360/nslit-closed_angle;
sel_phase = open_angle*rand01();
p *= (open_angle/(closed_angle+open_angle));
PROP_DT(dt);
if (x<xmin || x>xmax || y<ymin || y>ymax)
ABSORB; /* because outside frame */
/* now let's look whether the neutron is still
between the same two spokes or absorbed meanwhile */
n_angle_out = atan2(x,y+radius)*RAD2DEG; /* neutron beam might be divergent */
sel_phase = sel_phase + nu*dt*360 - alpha; /* rotor turned, but spokes are torsaded */
if (n_angle_out<(n_angle_in-sel_phase) || n_angle_out>(n_angle_in-sel_phase+open_angle) )
ABSORB; /* because must have passed absorber */
else
SCATTER;
%}
MCDISPLAY
%{
double phi, r0, Width, height, l0, l1;
double r;
double x0;
double x1;
double y0;
double y1;
double z0;
double z1;
double z2;
double z3;
double a;
double xw, yh;
phi = alpha;
Width = (xmax-xmin)/2;
height = (ymax-ymin)/2;
x0 = xmin; x1 = xmax;
y0 = ymin; y1 = ymax;
l0 = length; l1 = l0;
r0 = radius;
r = r0 + height;
x0 = -Width/2.0;
x1 = Width/2.0;
y0 = -height/2.0;
y1 = height/2.0;
z0 = 0;
z1 = 0;
z2 = l1;
z3 = l0;
xw = Width/2.0;
yh = height/2.0;
/* Draw apertures */
for(a = z0;;)
{
multiline(3, x0-xw, (double)y1, a,
(double)x0, (double)y1, a,
(double)x0, y1+yh, a);
multiline(3, x1+xw, (double)y1, a,
(double)x1, (double)y1, a,
(double)x1, y1+yh, a);
multiline(3, x0-xw, (double)y0, a,
(double)x0, (double)y0, a,
(double)x0, y0-yh, a);
multiline(3, x1+xw, (double)y0, a,
(double)x1, (double)y0, a,
(double)x1, y0-yh, a);
if(a == z3)
break;
else
a = z3;
}
/* Draw cylinder. */
circle("xy", 0, -r0, z1, r);
circle("xy", 0, -r0, z2, r);
line(0, -r0, z1, 0, -r0, z2);
for(a = 0; a < 2*PI; a += PI/8)
{
multiline(4,
0.0, -r0, z1,
r*cos(a), r*sin(a) - r0, z1,
r*cos(a + DEG2RAD*phi), r*sin(a + DEG2RAD*phi) - r0, z2,
0.0, -r0, z2);
}
/*
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); */
%}
END
|
