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
|
/*****************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright 1997-2002, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: Monochromator_pol
*
* %I
*
* Written by: Peter Christiansen
* Date: 2006
* Origin: RISOE
*
* Flat polarizaing monochromator crystal.
*
* %D
* Based on Monochromator_flat.
* Flat, infinitely thin mosaic crystal, useful as a monochromator or analyzer.
* For an unrotated monochromator component, the crystal surface lies in the Y-Z
* plane (ie. parallel to the beam).
* The mosaic and d-spread distributions are both Gaussian.
* Neutrons are just reflected (billard ball like). No correction is done for
* mosaicity of reflecting crystal.
* The crystal is assumed to be a ferromagnet with spin pointing up
* eta-tilde = (0, 1, 0) (along y-axis), so that the magnetic field is
* pointing opposite (0, -|B|, 0).
*
* The polarisation is done by defining the reflectivity for spin up
* (Rup) and spin down (Rdown) (which can be negative, see now!) and
* based on this the nuclear and magnetic structure factors are
* calculated:
* FM = sign(Rup)*sqrt(|Rup|) + sign(Rdown)*sqrt(|Rdown|)
* FN = sign(Rup)*sqrt(|Rup|) - sign(Rdown)*sqrt(|Rdown|)
* and the physics is calculated as
* Pol in = (sx_in, sy_in, sz_in)
* Reflectivity R0 = FN*FN + 2*FN*FM*sy_in + FM*FM
* (= |Rup| + |Rdown| (for sy_in=0))
* Pol out:
* sx = (FN*FN - FM*FM)*sx_in/R0;
* sy = ((FN*FN - FM*FM)*sy_in + 2*FN*FM + FM*FM*sy_in)/R0;
* sz = (FN*FN - FM*FM)*sz_in/R0;
*
* These equations are taken from:
* Lovesey: "Theory of neutron scattering from condensed matter, Volume
* 2", Eq. 10.96 and Eq. 10.110
*
* This component works with gravity (uses PROP_X0).
*
* Example: Monochromator_pol(zwidth=0.2, yheight=0.2, mosaic=30.0, dspread=0.0025, Rup=1.0, Rdown=0.0, Q=1.8734)
*
* Monochromator lattice parameter
* PG 002 DM=3.355 AA (Highly Oriented Pyrolythic Graphite)
* PG 004 DM=1.677 AA
* Heusler 111 DM=3.362 AA (Cu2MnAl)
* CoFe DM=1.771 AA (Co0.92Fe0.08)
* Ge 111 DM=3.266 AA
* Ge 311 DM=1.714 AA
* Ge 511 DM=1.089 AA
* Ge 533 DM=0.863 AA
* Si 111 DM=3.135 AA
* Cu 111 DM=2.087 AA
* Cu 002 DM=1.807 AA
* Cu 220 DM=1.278 AA
* Cu 111 DM=2.095 AA
*
* %P
* INPUT PARAMETERS:
*
* zwidth: [m] Width of crystal
* yheight: [m] Height of crystal
* mosaic: Mosaicity (FWHM) [arc minutes]
* dspread: [1] Relative d-spread (FWHM)
* Q: [AA-1] Magnitude of scattering vector
* Rup: [1] Reflectivity of neutrons with polarization up
* Rdown: [1] Reflectivity of neutrons with polarization down
*
* optional parameters
* DM: [AA] monochromator d-spacing instead of Q = 2*pi/DM
* pThreshold: [] if probability>pThreshold then accept and weight else random
* debug: [] if debug > 0, print out some info about the calculations
*
* %E
*****************************************************************************/
DEFINE COMPONENT Monochromator_pol
SETTING PARAMETERS (zwidth=0.1, yheight=0.1, mosaic=30.0, dspread=0, Q=1.8734, DM=0, pThreshold=0, Rup = 1, Rdown =1, int debug=0)
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
SHARE
%{
%include "pol-lib"
%}
DECLARE
%{
double mos_rms; /* root-mean-square of mosaic in radians */
double d_rms; /* root-mean-square of d-spread in AA */
double mono_Q;
double FN; /* Unit cell nuclear structure factor */
double FM; /* Unit cell magnetic structure factor */
%}
INITIALIZE
%{
mos_rms = MIN2RAD*mosaic/sqrt(8*log(2));
mono_Q = Q;
if (DM != 0)
mono_Q = 2*PI/DM;
DM = 2*PI/mono_Q;
d_rms = dspread*DM/sqrt(8*log(2));
// calculate the unit cell nuclear and magnetic structure factors
if(debug > 0)
printf("Rup: %f, Rdown: %f\n", Rup, Rdown);
GetMonoPolFNFM(Rup, Rdown, &FN, &FM);
if(debug > 0)
printf("FN: %f, FM: %f\n", FN, FM);
%}
TRACE
%{
double y1, z1, t1, dt, vel;
double sinTheta, lambdaBragg, lambda, dlambda2, sigmaLambda2, p_reflect;
double R0; /* reflection probability based on FN and FM */
double sx_in, sy_in, sz_in;
int i;
/* Propagate to crystal */
PROP_X0;
if (inside_rectangle(z, y, zwidth, yheight)) {/* Intersect the crystal? */
// calculate sin(Bragg angle)
vel = sqrt(vx*vx + vy*vy + vz*vz);
sinTheta = abs(vx)/vel;
// calculate lambdaBragg
lambdaBragg = 2.0*DM*sinTheta;
// calculate lambda of neutron
lambda = 2*PI/(V2K*vel);
// calculate deltalambda squared and sigmaLambda squared
dlambda2 = (lambda-lambdaBragg)*(lambda-lambdaBragg);
// The sigmaLambda is propagated by differentiating the Bragg
// condition: lambda = 2*d*sinTheta
sigmaLambda2 = 2.0*2.0 * sinTheta*sinTheta * d_rms*d_rms+
2.0*2.0 * DM*DM * (1.0-sinTheta*sinTheta) * mos_rms*mos_rms;
// calculate peak reflection probability
GetMonoPolRefProb(FN, FM, sy, &R0);
// calculate reflection probability
p_reflect = R0*exp(-dlambda2/(2.0*sigmaLambda2));
if(debug > 0) {
printf("\n lambda: %f, Lambda_Bragg: %f\n", lambda, lambdaBragg);
printf("sigmaLambda: %f, R0: %f, p_reflect: %f\n",
sqrt(sigmaLambda2), R0, p_reflect);
printf("S_in: (%f, %f, %f)\n", sx, sy, sz);
}
if((pThreshold>0 && p_reflect>pThreshold) || rand01()<p_reflect) {
/* Reflect */
// scale weight if neutron was accepted because of threshold
if(pThreshold>0 && p_reflect>pThreshold)
p*=p_reflect;
vx = -vx;
// Outgoing polarisation
SetMonoPolRefOut(FN, FM, R0, &sx, &sy, &sz);
if(debug > 0)
printf("S_out: (%f, %f, %f)\n", sx, sy, sz);
if(sx*sx+sy*sy+sz*sz>1)
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 MC choice to reflect or transmit neutron */
} /* End intersect the crystal */
%}
MCDISPLAY
%{
rectangle("yz", 0, 0, 0, zwidth, yheight);
%}
END
|
