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
|
/*******************************************************************************
*
* 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
* Date: August 4th 2010
* Origin: ILL
*
* Reflecting surface (guide and mirror) with any shape, defined from an OFF file.
*
* %D
*
* This is a reflecting object component. Its shape is 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'.
*
* The reflectivity is specified either from the usual parametric description
* R0,Qc,alpha,W,m, or from a reflectivity file 'reflect' with a 2 column
* format [q(Angs-1) R(0-1)].
*
* 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).
* Such 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
* reflect: [q(Angs-1) R(0-1)] (str) Reflectivity file name. Format
* 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
SETTING PARAMETERS (xwidth=0, yheight=0, zdepth=0, center=0, transmit=0,
R0=0.99, Qc=0.0219, alpha=3, m=2, W=0.003, string reflect=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 "interoff-lib"
%include "ref-lib"
%}
DECLARE
%{
off_struct offdata;
t_Table pTable;
%}
INITIALIZE
%{
/* initialize OFF object from the file(s) */
if (!off_init( geometry, xwidth, yheight, zdepth, !center, &offdata )) exit(-1);
/* optionally initialize reflectivity curves */
if (reflect && strlen(reflect) && strcmp(reflect,"NULL") && strcmp(reflect,"0")) {
if (Table_Read(&pTable, reflect, 1) <= 0) /* read 1st block data from file into pTable */
exit(fprintf(stderr,"Guide_anyshape: %s: can not read file %s. Aborting.\n", NAME_CURRENT_COMP, reflect));
else
Table_Rebin(&pTable); /* rebin as evenly, increasing array */
}
%}
TRACE
%{
int intersect=0;
int counter=0;
/* main loop for multiple reflections */
#ifdef OPENACC
off_struct thread_offdata = offdata;
#else
#define thread_offdata offdata
#endif
do {
double t0=0, t3=0, dt=0;
Coords n0, n3, n={0,0,0};
/* determine intersections with object */
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 = off_intersect(&t0, &t3, &n0, &n3, 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; }
if (intersect > 1 && dt <= 0 && t3 > dt) { dt = t3; n=n3; }
/* exit loop when no intersection forward */
if (dt <= 0 || !intersect) break;
double nx,ny,nz;
coords_get(n, &nx, &ny, &nz);
/* propagate neutron to reflection point, must be before n_dot_v with gravity */
PROP_DT(dt);
/* 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);
/* handle surface intersection */
double R=0;
/* Reflectivity (see component Guide). */
if (reflect && strlen(reflect) && strcmp(reflect,"NULL") && strcmp(reflect,"0"))
TableReflecFunc(q, &pTable, &R);
else {
double par[] = {R0, Qc, alpha, m, W};
StdReflecFunc(q, par, &R);
}
if (R > 1) {
fprintf(stderr,"Guide_anyshape: %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 */
off_display(offdata);
%}
END
|
