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
|
/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright 1997-2002, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: Sqq_w_monitor
*
* %I
* Written by: Peter Willendrup
* Date: June-July, 2018
* Origin: DTU
*
* Monitor outputting a series of energy-planes in a subset of reciprocal space, spanned by
* scattering vectors qa(x,z) and qb(x,z) in the component x-z plane.
*
* %D
*
* Cylindrical monitor on the x-z plane outputting a series of energy-planes in a subset of reciprocal
* space, spanned by scattering vectors qa(x,z) and qb(x,z).
*
* The radius and yheight parameters are not used for propagation, but only to define the outgoing divergence
* limit considered when estimating k_f
*
* The assumption is that the "current" neutron represents the final state, whereas the incoming state
* is found by restoring the neutron state "index" components earlier.
*
* Example: Sqq_w_monitor(filename="output",radius=1, yheight=0.05, Emin=0,Emax=5,nE=11,nqa=100,nqb=100,qamin=1,qamax=10,qbmin=1qbmax=10, vix="vix", viy="viy", viz="viz")
* AT (0,0,0) RELATIVE sample
*
* %P
* INPUT PARAMETERS:
*
* radius: [m] Cylinder radius
* yheight: [m] Cylinder height
* qax: [1] x-component of 1st q-vector
* qaz: [1] z-component of 1st q-vector
* qbx: [1] x-component of 2nd q-vector
* qbz: [1] z-component of 2nd q-vector
* qamin: [AA^-1] Defines interval (qamin,qamax) where monitor measures in nqa bins
* qamax: [AA^-1] Defines interval (qamin,qamax) where monitor measures in nqa bins
* qbmin: [AA^-1] Defines interval (qbmin,qbmax) where monitor measures in nqb bins
* qbmax: [AA^-1] Defines interval (qbmin,qbmax) where monitor measures in nqb bins
* nqa: [int] Number of bins along qa direction
* nqb: [int] Number of bins along qb direction
* Emin: [meV] Defines the energy-transfer [Emin,Emax] window to monitor in nE bins
* Emax: [meV] Defines the energy-transfer [Emax,Emax] window to monitor in nE bins
* nE: [int] Number of energy slices
* vix: [string] Points to instrument-level USERVAR for reading an earlier x-velocity
* viy: [string] Points to instrument-level USERVAR for reading an earlier y-velocity
* viz: [string] Points to instrument-level USERVAR for reading an earlier z-velocity
* filename: [string] Base filename to use, nE+1 files will be output
* nowritefile: [1] If set, monitor will skip writing to disk
* nosum: [1] If set, monitor will skip writing the energy-summed array to disk
*
* CALCULATED PARAMETERS:
*
* M_N: [] 3D array of neutron counts
* M_p: [] 3D array of neutron weight counts
* N_p2: [] 3D array of second moments
* M_Ns: [] 2D array of neutron counts
* M_ps: [] 2D array of neutron weight counts
* N_p2s: [] 2D array of second moments
*
* %E
*******************************************************************************/
DEFINE COMPONENT Sqq_w_monitor
SETTING PARAMETERS (radius=1, yheight=0.05, qax=1,qaz=0,qbx=0,qbz=1,qamin=0,qamax=2,qbmin=0,qbmax=2,Emin=0,Emax=5,int nqa=90, int nqb=90, int nE=10, string filename=0, int nowritefile=0, int nosum=0, string vix="", string viy="", string viz="")
DECLARE
%{
DArray3d M_N;
DArray3d M_p;
DArray3d M_p2;
DArray2d M_Ns;
DArray2d M_ps;
DArray2d M_p2s;
double dE;
%}
INITIALIZE
%{
/* Make checks for limits on qa, qb, w grid */
M_N=create_darr3d(nE,nqa,nqb);
M_p=create_darr3d(nE,nqa,nqb);
M_p2=create_darr3d(nE,nqa,nqb);
M_Ns=create_darr2d(nqa,nqb);
M_ps=create_darr2d(nqa,nqb);
M_p2s=create_darr2d(nqa,nqb);
dE=(Emax-Emin)/(1.0*nE-1);
// Use instance name for monitor output if no input was given
if (!strcmp(filename,"\0")) sprintf(filename,"%s",NAME_CURRENT_COMP);
%}
TRACE
%{
int i,j,k;
//double rx,ry,rz,rvx,rvy,rvz,rt,rsx,rsy,rsz,rp;
double rvx,rvy,rvz;
double Ei,Ef,E,Ki,Kf;
double qx,qy,qz;
double qqa, qqb;
double kix,kiy,kiz;
double kfx,kfy,kfz;
double t0,t1;
double nx,ny,nz;
nx=vx;
ny=vy;
nz=vz;
NORM(nx,ny,nz);
double v_div=scalar_prod(nx, ny, nz, 0, 1, 0);
/* Initial check to see if this neutron should be counted or not... */
if (fabs(v_div)<=yheight/radius && cylinder_intersect(&t0, &t1, x, y, z, vx, vy, vz, radius, yheight)) {
if(t0<0 && t1>0) {
/* This one hits our cylindrical band, arriving from inside */
int fail;
rvx = particle_getvar(_particle,vix,&fail); if(fail) rvx=0;
rvy = particle_getvar(_particle,viy,&fail); if(fail) rvy=0;
rvz = particle_getvar(_particle,viz,&fail); if(fail) rvz=0;
/* Calculate energy transfer */
Ei = VS2E*(rvx*rvx + rvy*rvy + rvz*rvz);
Ef = VS2E*( vx*vx + vy*vy + vz*vz);
E=Ei-Ef;
/* calculate k vectors and momentum transfer*/
kix=rvx;
kiy=rvy;
kiz=rvz;
kfx=vx;
kfy=vy;
kfz=vz;
NORM(kix, kiy, kiz);
NORM(kfx, kfy, kfz);
/* K-vector lengths */
Ki=V2K*sqrt((rvx*rvx)+(rvy*rvy)+(rvz*rvz));
Kf=V2K*sqrt((vx*vx)+(vy*vy)+(vz*vz));
kix=Ki*kix; kiy=Ki*kiy; kiz=Ki*kiz;
kfx=Kf*kfx; kfy=Kf*kfy; kfz=Kf*kfz;
qx=kix-kfx;
qy=kiy-kfy;
qz=kiz-kfz;
/* Calculate projections along qa and qb */
qqa = qx*qax + qz*qaz;
qqb = qx*qbx + qz*qbz;
/* Check if we are within the selected q/e range */
if (qqa <= qamax && qqa >=qamin && qqb <=qbmax && qqb>=qbmin && E<= Emax && E>=Emin) {
i = floor((qqa - qamin)*nqa/(qamax - qamin));
j = floor((qqb - qbmin)*nqb/(qbmax - qbmin));
k = floor((E - Emin)*nE/(Emax - Emin));
double p2 = p*p;
#pragma acc atomic
M_N[k][i][j] = M_N[k][i][j] + 1;
#pragma acc atomic
M_p[k][i][j] = M_p[k][i][j] + p;
#pragma acc atomic
M_p2[k][i][j] = M_p2[k][i][j] + p2;
#pragma acc atomic
M_Ns[i][j] = M_Ns[i][j] + 1 ;
#pragma acc atomic
M_ps[i][j] = M_ps[i][j] + p;
#pragma acc atomic
M_p2s[i][j] = M_p2s[i][j] + p2;
SCATTER;
}
}
}
/* Stuff that don't hit cylinder properly is disregarded */
%}
SAVE
%{
if (!nowritefile) {
int kk,ll;
char ff[256];
char tt[256];
if (!nosum) {
sprintf(ff, "%s_Sum",filename);
DETECTOR_OUT_2D(
"qa vs qb monitor E Sum",
"qa [AA^-1]",
"qb [AA^-1]",
qamin, qamax, qbmin, qbmax,
nqa, nqb,
&M_Ns[0][0],&M_ps[0][0],&M_p2s[0][0],
ff);
}
for (kk=0; kk<nE; kk++) {
sprintf(ff, "%s_%i",filename,kk);
sprintf(tt, "qa vs qb monitor E slice %i ~ %g meV",kk,dE*kk+Emin);
DETECTOR_OUT_2D(
tt,
"qa [AA^-1]",
"qb [AA^-1]",
qamin, qamax, qbmin, qbmax,
nqa, nqb,
&M_N[kk][0][0],&M_p[kk][0][0],&M_p2[kk][0][0],
ff);
}
}
%}
FINALLY %{
destroy_darr3d(M_N);
destroy_darr3d(M_p);
destroy_darr3d(M_p2);
destroy_darr2d(M_Ns);
destroy_darr2d(M_ps);
destroy_darr2d(M_p2s);
%}
MCDISPLAY
%{
circle("xz", 0, -(yheight/2.0), 0, radius );
circle("xz", 0, (yheight/2.0), 0, radius );
%}
END
|
