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
|
/****************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright 1997-2003, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: Pol_mirror
*
* %I
* Written by: Peter Christiansen
* Date: July 2006
* Origin: RISOE
*
* Polarising mirror.
*
* %D
* This component models a rectangular infinitely thin mirror.
* For an unrotated component, the mirror surface lies in the Y-Z
* plane (ie. parallel to the beam).
* It relies on similar physics as the Monochromator_pol.
* The reflectivity function (see e.g. share/ref-lib for examples) and parameters
* are passed to this component to give a bigger freedom.
* The up direction is hardcoded to be along the y-axis (0, 1, 0)
* IMPORTANT: At present the component only works correctly for polarization along the up/down
* direction and for completely unpolarized beams, i.e. sx=sy=sz=0 for all rays.
*
* For now we assume:
* P(Transmit|Q) = 1 - P(Reflect|Q)
* i.e. NO ABSORPTION!
*
*
* The component can both reflect and transmit neutrons with a respective proportion
* depending on the p_reflect parameter:
* p_reflect=-1 Reflect and transmit (proportions given from reflectivity) [default]
* p_reflect=1 Only handle reflected events
* p_reflect=0 Only handle transmitted events (reduce weight)
* p_reflect=0-1 Both transmit and reflect with fixed statistics proportions
*
* The parameters can either be
* double pointer initializations (e.g. {R0, Qc, alpha, m, W})
* or table names (e.g."supermirror_m2.rfl" AND useTables=1).
* NB! This might cause warnings by the compiler that can be ignored.
*
* Examples:
* Reflection function parametrization
* Pol_mirror(zwidth = 0.40, yheight = 0.40, rUpFunc=StdReflecFunc, rUpPar={1.0, 0.0219, 6.07, 2.0, 0.003})
*
* Table function
* Pol_mirror(zwidth = 0.40, yheight = 0.40, rUpFunc=TableReflecFunc, rUpPar="supermirror_m2.rfl", rDownFunc=TableReflecFunc, rDownPar="supermirror_m3.rfl", useTables=1)
*
* See also the example instrument Test_Pol_Mirror (under tests).
*
* GRAVITY: YES
*
* %BUGS
* NO ABSORPTION
*
* %P
* INPUT PARAMETERS:
*
* zwidth: [m] Width of the mirror
* yheight: [m] Height of the mirror
* rUpFunc: [1] Reflection function for spin up (q, *par, *r)
* rUpPar: [1] Parameters for rUpFunc
* rDownFunc: [1] Reflection function for spin down (q, *par, *r)
* rDownPar: [1] Parameters for rDownFunc
* rUpData: [str] Mirror Reflectivity data file for spin up
* rDownData: [str] Mirror Reflectivity data file for spin down
* useTables: [1] Parameters are 0: Values, 1: Table names
* p_reflect: [1] Proportion of reflected events. Use 0 to only get the transmitted beam, and 1 to get only the reflected beam. Use -1 to use the mirror reflectivity. This value is purely computational and is not related to the actual reflectivity
*
* CALCULATED PARAMETERS:
* SCATTERED: [] is 1 for reflected, and 2 for transmitted neutrons
*
* %L
*
* %E
*******************************************************************************/
DEFINE COMPONENT Pol_mirror
SETTING PARAMETERS (
vector rUpPar={0.99,0.0219,6.07,2.0,0.003},
vector rDownPar={0.99,0.0219,6.07,2.0,0.003},
string rUpData="", string rDownData="",
p_reflect=-1,
zwidth,
yheight)
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
SHARE
%{
%include "pol-lib"
%include "ref-lib"
%}
DECLARE
%{
t_Table rUpTable;
int rUpTableFlag;
t_Table rDownTable;
int rDownTableFlag;
%}
INITIALIZE
%{
if (strlen(rUpData) && strcmp(rUpData,"NULL")){
if (Table_Read(&rUpTable, rUpData, 1) <= 0) {
fprintf(stderr,"Pol_mirror: %s: can not read file %s\n",
NAME_CURRENT_COMP, rUpPar);
exit(1);
}
rUpTableFlag=1;
}else{
rUpTableFlag=0;
}
if (strlen(rUpData) && strcmp(rUpData,"NULL")){
if (Table_Read(&rDownTable, rDownData, 1) <= 0) {
fprintf(stderr,"Pol_mirror: %s: can not read file %s\n",
NAME_CURRENT_COMP, rDownPar);
exit(1);
}
rDownTableFlag=1;
}else{
rDownTableFlag=0;
}
if ((zwidth<=0) || (yheight <= 0)) {
fprintf(stderr, "Pol_mirror: %s: NULL or negative length scale!\n"
"ERROR (zwidth=%g,yheight=%g). Exiting\n",
NAME_CURRENT_COMP, zwidth, yheight);
exit(1);
}
%}
TRACE
%{
double Q, Rup, Rdown, FN, FM, refWeight;
int reflect = -1;
int isPolarising = 0;
// propagate to mirror plane
PROP_X0;
if (inside_rectangle(z, y, zwidth, yheight)) {/* Intersect the crystal? */
// calculate scattering vector magnitude
Q = fabs(2*vx*V2K);
// calculate reflection probability
// downgraded from defpars version
//rUpFunc(Q, rUpParPtr, &Rup);
if(rUpTableFlag){
Rup=Table_Value(rUpTable,Q,1);
}else{
StdReflecFunc(Q, rUpPar, &Rup);
}
// downgraded from defpars version
//rDownFunc(Q, rDownParPtr, &Rdown);
if(rDownTableFlag){
Rdown=Table_Value(rDownTable,Q,1);
}else{
StdReflecFunc(Q, rDownPar, &Rdown);
}
if(Rup != Rdown) {
isPolarising = 1;
GetMonoPolFNFM(Rup, Rdown, &FN, &FM);
GetMonoPolRefProb(FN, FM, sy, &refWeight);
/* Output of PW discussions with Hal Lee 2024/03/08
We have now done our QM "measurement", thus
forcing the spin to assume up/down: */
sx=0; sz=0;
} else
refWeight = Rup;
// check that refWeight is meaningfull
if (refWeight < 0) ABSORB;
if (refWeight > 1) refWeight =1 ;
// find out if neutrons is reflected or transmitted
if (p_reflect < 0 || p_reflect > 1) { // reflect OR transmit using mirror reflectivity
if (rand01()<refWeight) //reflect
reflect = 1;
else
reflect = 0;
} else { // reflect OR transmit using a specified weighting
if (rand01()<p_reflect) {
reflect = 1; // reflect
p *= refWeight/p_reflect;
} else {
reflect = 0; // transmit
p *= (1-refWeight)/(1-p_reflect);
}
}
// set outgoing velocity and polarisation
if (reflect==1) { // reflect: SCATTERED==1 for reflection
vx = -vx;
if(isPolarising)
SetMonoPolRefOut(FN, FM, refWeight, &sx, &sy, &sz);
} else { // transmit: SCATTERED==2 for transmission
if(isPolarising)
SetMonoPolTransOut(FN, FM, refWeight, &sx, &sy, &sz);
SCATTER;
}
if(isPolarising)
if(sx*sx+sy*sy+sz*sz>1.000001)
fprintf(stderr,"Pol_mirror: %s: Warning: polarisation |s| = %g > 1\n",
NAME_CURRENT_COMP, sx*sx+sy*sy+sz*sz); // check that polarisation is meaningfull
SCATTER;
} /* End intersect the mirror */
else
{
/* neutron will miss the mirror, so don't propagate i.e restore it */
RESTORE_NEUTRON(INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p);
}
%}
MCDISPLAY
%{
rectangle("yz", 0, 0, 0, zwidth, yheight);
%}
END
|
