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
|
/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright 1997-2002, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: disp_sample
*
* %I
* Written by: Garrett Granroth
* Date: 14.11.14
* Origin: Oak Ridge National Laboratory
*
* Spot sample.
*
* %D
* A sample that is a series of delta functions in omega and Q. The Q is
* determined by a two theta value and an equal number of spots around the
* beam center. This component is a variation of the V sample written
* by Kim Lefmann and Kristian Nielsen. The following text comes from their
* original component.
* A Double-cylinder shaped incoherent scatterer (a V-sample)
* No multiple scattering. Absorbtion included.
* The shape of the sample may be a box with dimensions xwidth, yheight, zthick.
* The area to scatter to is a disk of radius 'focus_r' situated at the target.
* This target area may also be rectangular if specified focus_xw and focus_yh
* or focus_aw and focus_ah, respectively in meters and degrees.
* The target itself is either situated according to given coordinates (x,y,z), or
* setting the relative target_index of the component to focus at (next is +1).
* This target position will be set to its AT position. When targeting to centered
* components, such as spheres or cylinders, define an Arm component where to
* focus at.
*
* Example: spot_sample(radius_o=0.01, h=0.05, pack = 1,
* xwidth=0, yheight=0, zthick=0, Eideal=100.0,w=50.0,two_theta=25.0,n_spots=4)
*
* %P
* INPUT PARAMETERS:
*
* radius_o: [m] Outer radius of sample in (x,z) plane
* h: [m] Height of sample y direction
* pack: [1] Packing factor
* Eideal: [meV] The presumed incident energy
* w: [meV] The energy transfer of the delta function
* two_theta: [degrees] the scattering angle of the spot
* n_spots: [1] The number of spots to generate symmetrically around the beam
*
* Optional parameters
* xwidth: [m] horiz. dimension of sample, as a width
* yheight: [m] vert.. dimension of sample, as a height
* zthick: [m] thickness of sample
*
* Variables calculated in the component
*
* V_my_s: [] Attenuation factor due to scattering (m^-1)
* V_my_a: [] Attenuation factor due to absorbtion (m^-1)
*
* %L
*
*
*
* %E
*******************************************************************************/
DEFINE COMPONENT Spot_sample
SETTING PARAMETERS (radius_o=0.01, h=0.05, pack = 1,
xwidth=0, yheight=0, zthick=0, Eideal=100.0,w=50.0,two_theta=25.0,n_spots=4)
//STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p) removed to go to mcstas version 2
SHARE
%{
struct StructVarsVspot
{
double sigma_a; /* Absorption cross section per atom (barns) */
double sigma_i; /* Incoherent scattering cross section per atom (barns) */
double rho; /* Density of atoms (AA-3) */
double my_s;
double my_a_v;
char isrect; /* true when sample is a box */
};
%}
DECLARE %{
struct StructVarsVspot VarsV;
%}
INITIALIZE
%{
if (!radius_o || !h) {
if (!xwidth || !yheight || !zthick) exit(fprintf(stderr,"V_sample: %s: sample has no volume (zero dimensions)\n", NAME_CURRENT_COMP));
else VarsV.isrect=1; }
else VarsV.isrect=0;
VarsV.sigma_a=5.08; /* in barns */
VarsV.sigma_i=4.935;
VarsV.rho = (2*pack/(3.024*3.024*3.024));
VarsV.my_s=(VarsV.rho * 100 * VarsV.sigma_i);
VarsV.my_a_v=(VarsV.rho * 100 * VarsV.sigma_a * 2200);
/* now compute target coords if a component index is supplied */
%}
TRACE
%{
double t0, t3; /* Entry/exit time for outer cylinder */
double t1, t2; /* Entry/exit time for inner cylinder */
double v; /* Neutron velocity */
double dt0, dt1, dt2, dt; /* Flight times through sample */
double l_full; /* Flight path length for non-scattered neutron */
double l_i, l_o=0; /* Flight path lenght in/out for scattered neutron */
double my_a; /* Velocity-dependent attenuation factor */
double solid_angle=0; /* Solid angle of target as seen from scattering point */
double aim_x, aim_y, aim_z; /* Position of target relative to scattering point */
double kix,kiy,kiz,qx,qy,qz;
double kf,kfx,kfy,kfz,kiideal,kfideal,Efideal;
double Ef,Ei,pol,rbool;
int spot;
int intersect=0;
if (VarsV.isrect)
intersect = box_intersect(&t0, &t3, x, y, z, vx, vy, vz, xwidth, yheight, zthick);
else
intersect = cylinder_intersect(&t0, &t3, x, y, z, vx, vy, vz, radius_o, h);
if(intersect)
{
if(t0 < 0) ABSORB; /* we already passed the sample */
/* Neutron enters at t=t0. */
dt0 = t3-t0; /* Time in sample, */
v = sqrt(vx*vx + vy*vy + vz*vz);
kix=vx*V2K;kiy=vy*V2K;kiz=vz*V2K;
Ei=v*v*VS2E;
l_full = v * (dt0); /* Length of full path through sample */
dt = rand01()*(dt0); /* Time of scattering (relative to t0) */
l_i = v*dt; /* Penetration in sample */
PROP_DT(dt+t0); /* Point of scattering */
Efideal=Eideal-w;
kfideal=sqrt(Efideal/2.0723);
kiideal=SE2V*sqrt(Eideal)*V2K;
spot=floor(n_spots*rand01())+1;
pol=(spot-1)*2.0*PI/n_spots;
qz=kiideal-kfideal*cos(two_theta*DEG2RAD);
qx=-kfideal*cos(pol)*sin(two_theta*DEG2RAD);
qy=-kfideal*sin(pol)*sin(two_theta*DEG2RAD);
kfx=kix-qx;
kfy=kiy-qy;
kfz=kiz-qz;
if(!VarsV.isrect) {
if(!cylinder_intersect(&t0, &t3, x, y, z, vx, vy, vz, radius_o, h))
{
/* ??? did not hit cylinder */
printf("FATAL ERROR: Did not hit cylinder from inside.\n");
exit(1);
}
dt = t3;
}
vx = kfx*K2V;
vy = kfy*K2V;
vz = kfz*K2V;
/*printf("vx:%g vy:%g vz:%g \n", vx,vy,vz);*/
my_a = VarsV.my_a_v/v;
p*=1;
SCATTER;
}
%}
MCDISPLAY
%{
if (!VarsV.isrect) {
circle("xz", 0, h/2.0, 0, radius_o);
circle("xz", 0, -h/2.0, 0, radius_o);
line(-radius_o, -h/2.0, 0, -radius_o, +h/2.0, 0);
line(+radius_o, -h/2.0, 0, +radius_o, +h/2.0, 0);
line(0, -h/2.0, -radius_o, 0, +h/2.0, -radius_o);
line(0, -h/2.0, +radius_o, 0, +h/2.0, +radius_o);
}
else
{
double xmin = -0.5*xwidth;
double xmax = 0.5*xwidth;
double ymin = -0.5*yheight;
double ymax = 0.5*yheight;
double zmin = -0.5*zthick;
double zmax = 0.5*zthick;
multiline(5, xmin, ymin, zmin,
xmax, ymin, zmin,
xmax, ymax, zmin,
xmin, ymax, zmin,
xmin, ymin, zmin);
multiline(5, xmin, ymin, zmax,
xmax, ymin, zmax,
xmax, ymax, zmax,
xmin, ymax, zmax,
xmin, ymin, zmax);
line(xmin, ymin, zmin, xmin, ymin, zmax);
line(xmax, ymin, zmin, xmax, ymin, zmax);
line(xmin, ymax, zmin, xmin, ymax, zmax);
line(xmax, ymax, zmin, xmax, ymax, zmax);
}
%}
END
|
