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
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
|
/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright (C) 1997-2011, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: Guide
*
* %I
* Written by: Kristian Nielsen
* Modified by: Nikolaos Tsapatsaris to accept differently coated mirror planes
* Date: September 2 1998, Modified 2013
* Version: $Revision: 1.32 $
* Origin: Risoe
* Release: McStas 1.12c
*
* Neutron guide.
*
* %D
* Models a rectangular guide tube centered on the Z axis. The entrance lies
* in the X-Y plane.
* For details on the geometry calculation see the description in the McStas
* reference manual.
* The reflectivity profile may either use an analytical mode (see Component
* Manual) or a 2-columns reflectivity free text file with format
* [q(Angs-1) R(0-1)].
*
* Example: Guide(w1=0.1, h1=0.1, w2=0.1, h2=0.1, l=2.0,
* R0=0.99, Qc=0.021, alpha=6.07, m=2, W=0.003
*
* %VALIDATION
* May 2005: extensive internal test, no bugs found
* Validated by: K. Lieutenant
*
* %BUGS
* This component does not work with gravitation on. Use component Guide_gravity then.
*
* %P
* INPUT PARAMETERS:
*
* w1: [m] Width at the guide entry
* h1: [m] Height at the guide entry
* w2: [m] Width at the guide exit
* h2: [m] Height at the guide exit
* l: [m] Length of guide
* R0_left: [1] Low-angle reflectivity, left mirrror.
* Qc_left: [AA-1] Critical scattering vector, left mirrror.
* alpha_left: [AA] Slope of reflectivity, left mirrror.
* m_left: [1] m-value of material. Zero means completely absorbing, left mirrror.
* W_left: [AA-1] Width of supermirror cut-off, left mirrror.
* R0_right: [1] Low-angle reflectivity, right mirror.
* Qc_right: [AA-1] Critical scattering vector, right mirror.
* alpha_right: [AA] Slope of reflectivity, right mirror.
* m_right: [1] m-value of material. Zero means completely absorbing, right mirror.
* W_right: [AA-1] Width of supermirror cut-off, right mirror.
* R0_top: [1] Low-angle reflectivity, top mirror.
* Qc_top: [AA-1] Critical scattering vector, top mirror.
* alpha_top: [AA] Slope of reflectivity, top mirror.
* m_top: [1] m-value of material. Zero means completely absorbing, top mirror.
* W_top: [AA-1] Width of supermirror cut-off, top mirror.
* R0_bottom: [1] Low-angle reflectivity, bottom mirror.
* Qc_bottom: [AA-1] Critical scattering vector, bottom mirror.
* alpha_bottom: [AA] Slope of reflectivity, bottom mirror.
* m_bottom: [1] m-value of material. Zero means completely absorbing, bottom mirror.
* W_bottom: [AA-1] Width of supermirror cut-off, bottom mirror.
* reflect: [str] Reflectivity file name. Format [q(Angs-1) R(0-1)]
*
* %D
* Example values: m=4 Qc=0.0219 W=1/300 alpha=6.49 R0=1
*
* %E
*
* Mirror reflectivity files
* For the simulation of the guide surfaces Mirrotron reflectivity files [6] were used as a template for creating reflectivity files in VITESS using the following mathematical description:
* R = 1/2 R0(1 - tanh((Q - m Qc,Ni)/W) (1 - α(Q - Qc))
* using
* α = (Rm - R0) / (m Qc,Ni - Qc).
*
*
* m = Qmax/Qmax(Ni) Reflectivity at Q=m*Qc(Ni) Width of cut-off Alpha
* 1 0.99 0.003 0.000
* 1.5 0.97 0.003 1.904
* 2 0.95 0.003 1.904
* 2.5 0.92 0.003 2.222
* 3.6 0.78 0.003 3.846
* 4 0.75 0.003 3.333
*******************************************************************************/
DEFINE COMPONENT Guide_m
SETTING PARAMETERS (string reflect=0, w1, h1, w2, h2, l, R0_left=0.99, R0_right=0.99, R0_top=0.99, R0_bottom=0.99, Qc_left=0.0219, Qc_right=0.0219, Qc_top=0.0219, Qc_bottom=0.0219, alpha_left=6.07, alpha_right=6.07, alpha_top=6.07, alpha_bottom=6.07, m_left=2, m_right=2, m_top=2, m_bottom=2, W_left=0.003, W_right=0.003 , W_top=0.003 , W_bottom=0.003)
SHARE
%{
%include "read_table-lib"
%}
DECLARE
%{
t_Table pTable;
double m;
double alpha;
double Qc;
double R0;
double W;
%}
INITIALIZE
%{
if (mcgravitation) fprintf(stderr,"WARNING: Guide: %s: "
"This component produces wrong results with gravitation !\n"
"Use Guide_gravity.\n",
NAME_CURRENT_COMP);
if (reflect && strlen(reflect)) {
if (Table_Read(&pTable, reflect, 1) <= 0) /* read 1st block data from file into pTable */
exit(fprintf(stderr,"Guide: %s: can not read file %s\n", NAME_CURRENT_COMP, reflect));
} else {
if (W_left < 0 || W_right < 0 || W_top < 0 || W_bottom < 0 || R0_left < 0 ||R0_right < 0 || R0_top < 0 || R0_bottom < 0 || Qc_left < 0 || Qc_right < 0 || Qc_top < 0 || Qc_bottom < 0 || m_left < 0 || m_right < 0 || m_top < 0 || m_bottom < 0)
{ fprintf(stderr,"Guide: %s: W R0 Qc must be >0.\n", NAME_CURRENT_COMP);
exit(-1); }
if (m_left < 1 && m_left != 0) fprintf(stderr,"WARNING: Guide: %s: m_left < 1 behaves as if m=1.\n",
NAME_CURRENT_COMP);
if (m_right < 1 && m_right != 0) fprintf(stderr,"WARNING: Guide: %s: m_right < 1 behaves as if m=1.\n",
NAME_CURRENT_COMP);
if (m_top < 1 && m_top != 0) fprintf(stderr,"WARNING: Guide: %s: m_top < 1 behaves as if m=1.\n",
NAME_CURRENT_COMP);
if (m_bottom < 1 && m_bottom != 0) fprintf(stderr,"WARNING: Guide: %s: m_bottom < 1 behaves as if m=1.\n",
NAME_CURRENT_COMP);
}
%}
TRACE
%{
double t1,t2; /* Intersection times. */
double av,ah,bv,bh,cv1,cv2,ch1,ch2,d; /* Intermediate values */
double weight; /* Internal probability weight */
double vdotn_v1,vdotn_v2,vdotn_h1,vdotn_h2; /* Dot products. */
int i; /* Which mirror hit? */
double q; /* Q [1/AA] of reflection */
double nlen2; /* Vector lengths squared */
/* ToDo: These could be precalculated. */
double ww = .5*(w2 - w1), hh = .5*(h2 - h1);
double whalf = .5*w1, hhalf = .5*h1;
/* Propagate neutron to guide entrance. */
PROP_Z0;
/* Scatter here to ensure that fully transmitted neutrons will not be
absorbed in a GROUP construction, e.g. all neutrons - even the
later absorbed ones are scattered at the guide entry. */
SCATTER;
if(x <= -whalf || x >= whalf || y <= -hhalf || y >= hhalf)
ABSORB;
for(;;)
{
/* Compute the dot products of v and n for the four mirrors. */
av = l*vx; bv = ww*vz;
ah = l*vy; bh = hh*vz;
vdotn_v1 = bv + av; /* Left vertical */
vdotn_v2 = bv - av; /* Right vertical */
vdotn_h1 = bh + ah; /* Lower horizontal */
vdotn_h2 = bh - ah; /* Upper horizontal */
/* Compute the dot products of (O - r) and n as c1+c2 and c1-c2 */
cv1 = -whalf*l - z*ww; cv2 = x*l;
ch1 = -hhalf*l - z*hh; ch2 = y*l;
/* Compute intersection times. */
t1 = (l - z)/vz;
i = 0;
if(vdotn_v1 < 0 && (t2 = (cv1 - cv2)/vdotn_v1) < t1)
{
t1 = t2;
i = 1;
}
if(vdotn_v2 < 0 && (t2 = (cv1 + cv2)/vdotn_v2) < t1)
{
t1 = t2;
i = 2;
}
if(vdotn_h1 < 0 && (t2 = (ch1 - ch2)/vdotn_h1) < t1)
{
t1 = t2;
i = 3;
}
if(vdotn_h2 < 0 && (t2 = (ch1 + ch2)/vdotn_h2) < t1)
{
t1 = t2;
i = 4;
}
if(i == 0)
break; /* Neutron left guide. */
PROP_DT(t1);
switch(i)
{
case 1: /* Left vertical mirror */
nlen2 = l*l + ww*ww;
q = V2Q*(-2)*vdotn_v1/sqrt(nlen2);
d = 2*vdotn_v1/nlen2;
vx = vx - d*l;
vz = vz - d*ww;
m = m_left;
Qc = Qc_left;
W = W_left;
alpha= alpha_left;
R0= R0_left;
break;
case 2: /* Right vertical mirror */
nlen2 = l*l + ww*ww;
q = V2Q*(-2)*vdotn_v2/sqrt(nlen2);
d = 2*vdotn_v2/nlen2;
vx = vx + d*l;
vz = vz - d*ww;
m = m_right;
Qc = Qc_right;
W = W_right;
alpha= alpha_right;
R0= R0_right;
break;
case 3: /* Lower horizontal mirror */
nlen2 = l*l + hh*hh;
q = V2Q*(-2)*vdotn_h1/sqrt(nlen2);
d = 2*vdotn_h1/nlen2;
vy = vy - d*l;
vz = vz - d*hh;
m = m_bottom;
Qc = Qc_bottom;
W = W_bottom;
alpha= alpha_bottom;
R0= R0_bottom;
break;
case 4: /* Upper horizontal mirror */
nlen2 = l*l + hh*hh;
q = V2Q*(-2)*vdotn_h2/sqrt(nlen2);
d = 2*vdotn_h2/nlen2;
vy = vy + d*l;
vz = vz - d*hh;
m = m_top;
Qc = Qc_top;
W = W_top;
alpha= alpha_top;
R0= R0_top;
break;
}
/* Now compute reflectivity. */
weight = 1.0; /* Initial internal weight factor */
if(m == 0)
ABSORB;
if (reflect && strlen(reflect))
weight = Table_Value(pTable, q, 1);
else if(q > Qc)
{
double arg = (q-m*Qc)/W;
if(arg < 10)
weight = .5*(1-tanh(arg))*(1-alpha*(q-Qc));
else
ABSORB; /* Cutoff ~ 1E-10 */
weight *= R0;
} else { /* q <= Qc */
weight *= R0;
}
p *= weight;
SCATTER;
}
%}
MCDISPLAY
%{
multiline(5,
-w1/2.0, -h1/2.0, 0.0,
w1/2.0, -h1/2.0, 0.0,
w1/2.0, h1/2.0, 0.0,
-w1/2.0, h1/2.0, 0.0,
-w1/2.0, -h1/2.0, 0.0);
multiline(5,
-w2/2.0, -h2/2.0, (double)l,
w2/2.0, -h2/2.0, (double)l,
w2/2.0, h2/2.0, (double)l,
-w2/2.0, h2/2.0, (double)l,
-w2/2.0, -h2/2.0, (double)l);
line(-w1/2.0, -h1/2.0, 0, -w2/2.0, -h2/2.0, (double)l);
line( w1/2.0, -h1/2.0, 0, w2/2.0, -h2/2.0, (double)l);
line( w1/2.0, h1/2.0, 0, w2/2.0, h2/2.0, (double)l);
line(-w1/2.0, h1/2.0, 0, -w2/2.0, h2/2.0, (double)l);
%}
END
|
