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
|
/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright (C) 1997-2010, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: Guide
*
* %I
* Written by: Emmanuel Farhi, adapted by Peter Link and Gaetano Mangiapia
* Original Date: August 4th 2010
* Revised for McStas 3.1 on: May 31st 2022
* Origin: ILL/MLZ
*
* Reflecting surface (guide and mirror) with any shape, defined from an OFF file
* and a patch to allow m-, alpha- and W-values added per polygon face.
* Derived from Guide_anyshape / Guide_anyshape_r .
*
* %D
*
* This is a reflecting object component, derived from Guide_anyshape.
* Its shape, as well as m-, alpha- and W-values (IN THIS ORDER) for each face, are defined from
* an OFF file, given with its file name. The object size is set as given from the file, where
* dimensions should be in meters. The bounding box may be re-scaled by specifying
* xwidth,yheight,zdepth parameters. The object may optionally be centered when
* using 'center=1'.
*
* If in input only the values of R0 and Qc are specified, leaving m, alpha and W all null, then
* Guide_anyshape_g will use the values of these last three parameters that are specified in the OFF
* file for each face: floating point based values may be provided for them.
* In case at least one among m, alpha and W is not null, then Guide_anyshape_g will use these
* values as input (common for all the faces), ignoring those specified in the OFF file
*
* The complex OFF/PLY geometry option handles any closed non-convex polyhedra.
* It supports the OFF and NOFF file format but not COFF (colored faces).
* Standard .OFF files may be generated from XYZ data using:
* qhull < coordinates.xyz Qx Qv Tv o > geomview.off
* or
* powercrust coordinates.xyz
* and viewed with geomview or java -jar jroff.jar (see below).
* The default size of the object depends of the OFF file data, but its
* bounding box may be resized using xwidth,yheight and zdepth.
* PLY geometry files are also supported.
*
* %P
* INPUT PARAMETERS:
*
* geometry: [str] Name of the OFF/PLY file describing the guide geometry
* xwidth: [m] Redimension the object bounding box on X axis is non-zero
* yheight: [m] Redimension the object bounding box on Y axis is non-zero
* zdepth: [m] Redimension the object bounding box on Z axis is non-zero
* center: [1] When set to 1, the object will be centered w.r.t. the local coordinate frame
* R0: [1] Low-angle reflectivity
* Qc: [AA-1] Critical scattering vector
* alpha: [AA] Slope of reflectivity
* m: [1] m-value of material. Zero means completely absorbing.
* W: [AA-1] Width of supermirror cut-off
* transmit: [1] When true, non reflected neutrons are transmitted through the surfaces, instead of being absorbed. No material absorption is taken into account though
*
* CALCULATED PARAMETERS:
* SCATTERED: [] number of reflected events
*
* %D
* Example values: m=4 Qc=0.0219 W=1/300 alpha=6.49 R0=1
*
* %L
* <a href="http://www.geomview.org">Geomview and Object File Format (OFF)</a>
* %L
* Java version of Geomview (display only) <a href="http://www.holmes3d.net/graphics/roffview/">jroff.jar</a>
* %L
* <a href="http://qhull.org">qhull</a>
* %L
* <a href="http://www.cs.ucdavis.edu/~amenta/powercrust.html">powercrust</a>
*
* %E
*******************************************************************************/
DEFINE COMPONENT Guide_anyshape_r
SETTING PARAMETERS (xwidth=0, yheight=0, zdepth=0, center=0, transmit=0,
R0=0.99, Qc=0.0219, alpha=0, m=0, W=0, string geometry=0)
DEPENDENCY " -DUSE_OFF "
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
SHARE
%{
%include "read_table-lib"
%include "r-interoff-lib"
%include "ref-lib"
%}
DECLARE
%{
r_off_struct offdata;
int use_file_coatings; // Flag to specify whether to use (1) or not (0) the mirror parameters given in the .OFF file
%}
INITIALIZE
%{
/* initialize OFF object from the file(s) */
if (!r_off_init( geometry, xwidth, yheight, zdepth, !center, &offdata )) exit(-1);
if (m==0 && alpha==0 && W==0) {
use_file_coatings = 1;
printf("Guide_anyshape_r: %s: All of m, alpha, W assigned 0: \n - We are using H. Jacobsen reflectivity model if m-values are given in %s\n", NAME_CURRENT_COMP, geometry);
} else {
printf("Guide_anyshape_r: %s: alpha / m / W values provided in input: \n The corresponding values specified in %s will be ignored!\n", NAME_CURRENT_COMP, geometry);
use_file_coatings = 0;
}
%}
TRACE
%{
int intersect=0;
int counter=0;
/* main loop for multiple reflections */
#ifdef OPENACC
r_off_struct thread_offdata = offdata;
#else
#define thread_offdata offdata
#endif
do {
double t0=0, t3=0, dt=0;
unsigned long faceindex0=0, faceindex3=0, fi=0;
Coords n0, n3, n={0,0,0};
/* determine intersections with object; PL: get index of the intersected face */
double mc_gx=0.0, mc_gy=0.0, mc_gz=0.0;
if (mcgravitation) {
Coords locgrav;
locgrav = rot_apply(_comp->_rotation_absolute, coords_set(0,-GRAVITY,0));
coords_get(locgrav, &mc_gx, &mc_gy, &mc_gz);
}
intersect = r_off_intersect(&t0, &t3, &n0, &n3, &faceindex0, &faceindex3,
x,y,z, vx, vy, vz, mc_gx, mc_gy, mc_gz, thread_offdata );
/* get the smallest positive */
if (t0 > 0) { dt = t0; n=n0; fi=faceindex0;}
if (intersect > 1 && dt <= 0 && t3 > dt) { dt = t3; n=n3; fi=faceindex3;}
/* exit loop when no intersection forward */
if (dt <= 0 || !intersect) break;
double nx,ny,nz;
coords_get(n, &nx, &ny, &nz);
/* test if the angle is large in case the object has an internal coating */
double n2 = nx*nx+ny*ny+nz*nz;
double n_dot_v = scalar_prod(vx,vy,vz,nx,ny,nz);
double q = 2*fabs(n_dot_v)*V2K/sqrt(n2);
/* propagate neutron to reflection point */
PROP_DT(dt);
/* handle surface intersection */
double R=0;
double m_value = m;
double alpha_value = alpha;
double W_value = W;
if (use_file_coatings == 1) {
m_value = offdata.face_m_Array[fi];
alpha_value = offdata.face_alpha_Array[fi];
W_value = offdata.face_W_Array[fi];
}
double par[] = {R0, Qc, alpha_value, m_value, W_value};
StdReflecFunc(q, par, &R);
if (R > 1) {
fprintf(stderr,"Guide_anyshape_r: %s: Warning: Reflectivity R=%g > 1 lowered to R=1.\n", NAME_CURRENT_COMP, R);
R=1;
}
/* now handle either probability when transmit or reflect */
if (R > 0) {
/* when allowing transmission, we should check if indeed we reflect */
if (!transmit || (transmit && rand01() < R)) {
/* reflect velocity: -q -> -2*n*n.v/|n|^2 */
if (!transmit) p *= R;
n_dot_v *= 2/n2;
vx -= nx*n_dot_v;
vy -= ny*n_dot_v;
vz -= nz*n_dot_v;
SCATTER;
} else {
if (transmit) {
p *= (1-R); /* transmitted beam has non reflected weight */
} else ABSORB;
}
} else {
/* R=0: no reflection: absorb or transmit through when allowed */
if (!transmit) ABSORB;
}
/* leave surface by a small amount so that next intersection is not the same one */
PROP_DT(1e-9);
} while (intersect && counter++<CHAR_BUF_LENGTH);
/* end of main loop */
%}
MCDISPLAY
%{
/* display the object */
r_off_display(offdata);
%}
END
|
