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
|
/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright 1997-2002, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: V_selector
*
* %I
* Written by: Kim Lefmann
* Date: Nov 25, 1998
* Origin: Risoe
*
* Velocity selector.
*
* %D
* Velocity selector consisting of rotating Soller-like blades
* defining a helically twisted passage.
* Geometry defined by two identical, centered apertures at 12 o'clock
* position, Origo is at the centre of the selector (input is at -zdepth/2).
* Transmission is analytical assuming a continuous source.
*
* Example: V_selector(xwidth=0.03, yheight=0.05, zdepth=0.30, radius=0.12, alpha=48.298, length=0.25, d=0.0004, nu=20000, nslit=72)
* These are values for the D11@ILL Dornier 'Dolores' Velocity Selector (NVS 023)
*
* %VALIDATION
* Jun 2005: extensive external test, no problems found
* Validated by: K. Lieutenant
*
* %P
* INPUT PARAMETERS:
*
* xwidth: [m] Width of entry aperture
* yheight: [m] Height of entry aperture
* zdepth: [m] Distance between apertures, for housing containing the rotor
* radius: [m] Height from aperture centre to rotation axis
* alpha: [deg] Twist angle along the cylinder
* length: [m] Length of cylinder/rotor (less than zdepth)
* d: [m] Thickness of blades
* nu: [Hz] Cylinder rotation speed, counter-clockwise, which is ideally 3956*alpha*DEG2RAD/2/PI/lambda/length
* nslit: [1] Number of Soller blades
*
* %E
*******************************************************************************/
DEFINE COMPONENT V_selector
SETTING PARAMETERS (xwidth=0.03, yheight=0.05, zdepth=0.30, radius=0.12, alpha=48.298, length=0.25, d=0.0004, nu=300, nslit=72)
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
DECLARE
%{
double omega;
double alpha_rad;
%}
INITIALIZE
%{
omega=nu*2*PI;
alpha_rad = alpha*DEG2RAD;
if (zdepth < length) zdepth=length;
%}
TRACE
%{
double dt0,dt1;
double r_i,r_f,r_mean;
double theta_i,theta_f;
double A;
double d_s_alpha;
if (vz == 0)
ABSORB;
dt1= (-zdepth/2.0 - z)/vz;
PROP_DT(dt1); /* Propagate to the entry aperture */
if (x<(-xwidth/2.0) || x>(xwidth/2.0) || y<(-yheight/2.0) || y>(yheight/2.0))
ABSORB;
dt0 = (zdepth-length)/(2.0*vz); /* Propagate to the cylinder start */
PROP_DT(dt0);
r_i = sqrt(x*x+(y+radius)*(y+radius));
theta_i = atan2(x,y+radius);
dt1 = length/vz; /* Propagate along the cylinder length */
PROP_DT(dt1);
r_f = sqrt(x*x+(y+radius)*(y+radius));
theta_f = atan2(x,y+radius);
dt0 = (zdepth-length)/(2.0*vz); /* Propagate to the exit aperture */
PROP_DT(dt0);
if (x<(-xwidth/2.0) || x>(xwidth/2.0) || y<(-yheight/2.0) || y>(yheight/2.0))
ABSORB;
/* Calculate analytical transmission assuming continuous source */
r_mean = (r_i + r_f)/2.0; /* Approximation using mean radius */
d_s_alpha = theta_f-theta_i;
A = nslit/(2*PI)*( d/r_mean + fabs(alpha_rad+d_s_alpha-omega*length/vz) );
if (A >= 1)
ABSORB;
p*= (1-A);
SCATTER;
%}
MCDISPLAY
%{
double r = radius + yheight;
double x0 = -xwidth/2.0;
double x1 = xwidth/2.0;
double y0 = -yheight/2.0;
double y1 = yheight/2.0;
double z0 = -zdepth/2.0;
double z1 = -length/2.0;
double z2 = length/2.0;
double z3 = zdepth/2.0;
double a;
double xw, yh;
xw = xwidth/2.0;
yh = yheight/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, -radius, z1, r);
circle("xy", 0, -radius, z2, r);
line(0, -radius, z1, 0, -radius, z2);
for(a = 0; a < 2*PI; a += PI/8)
{
multiline(4,
0.0, -radius, z1,
r*cos(a), r*sin(a) - radius, z1,
r*cos(a + DEG2RAD*alpha), r*sin(a + DEG2RAD*alpha) - radius, z2,
0.0, -radius, z2);
}
%}
END
|
