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
|
/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright 1997-2002, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: Source_simple
*
* %I
* Written by: Kim Lefmann
* Date: October 30, 1997
* Modified by: KL, October 4, 2001
* Modified by: Emmanuel Farhi, October 30, 2001. Serious bug corrected.
* Origin: Risoe
*
* A circular neutron source with flat energy spectrum and arbitrary flux
*
* %D
* The routine is a circular neutron source, which aims at a square target
* centered at the beam (in order to improve MC-acceptance rate). The angular
* divergence is then given by the dimensions of the target.
* The neutron energy is uniformly distributed between lambda0-dlambda and
* lambda0+dlambda or between E0-dE and E0+dE.
* The flux unit is specified in n/cm2/s/st/energy unit (meV or Angs).
*
* This component replaces Source_flat, Source_flat_lambda,
* Source_flux and Source_flux_lambda.
*
* Example: Source_simple(radius=0.1, dist=2, focus_xw=.1, focus_yh=.1, E0=14, dE=2)
*
* %P
* radius: [m] Radius of circle in (x,y,0) plane where neutrons are generated.
* yheight: [m] Height of rectangle in (x,y,0) plane where neutrons are generated.
* xwidth: [m] Width of rectangle in (x,y,0) plane where neutrons are generated.
* target_index: [1] relative index of component to focus at, e.g. next is +1 this is used to compute 'dist' automatically.
* dist: [m] Distance to target along z axis.
* focus_xw: [m] Width of target
* focus_yh: [m] Height of target
* E0: [meV] Mean energy of neutrons.
* dE: [meV] Energy half spread of neutrons (flat or gaussian sigma).
* lambda0: [AA] Mean wavelength of neutrons.
* dlambda: [AA] Wavelength half spread of neutrons.
* flux: [1/(s*cm**2*st*energy unit)] flux per energy unit, Angs or meV if flux=0, the source emits 1 in 4*PI whole space.
* gauss: [1] Gaussian (1) or Flat (0) energy/wavelength distribution
*
* %E
*******************************************************************************/
DEFINE COMPONENT Source_simple
SETTING PARAMETERS (radius=0.1, yheight=0, xwidth=0,
dist=0, focus_xw=.045, focus_yh=.12,
E0=0, dE=0, lambda0=0, dlambda=0,
flux=1, gauss=0, int target_index=+1)
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
DECLARE
%{
double pmul;
double srcArea;
int square;
double tx;
double ty;
double tz;
%}
INITIALIZE
%{
square = 0;
/* Determine source area */
if (radius && !yheight && !xwidth ) {
square = 0;
srcArea = PI*radius*radius;
} else if(yheight && xwidth) {
square = 1;
srcArea = xwidth * yheight;
}
if (flux) {
pmul=flux*1e4*srcArea/mcget_ncount();
if (dlambda)
pmul *= 2*dlambda;
else if (dE)
pmul *= 2*dE;
} else {
gauss = 0;
pmul=1.0/(mcget_ncount()*4*PI);
}
if (target_index && !dist)
{
Coords ToTarget;
ToTarget = coords_sub(POS_A_COMP_INDEX(INDEX_CURRENT_COMP+target_index),POS_A_CURRENT_COMP);
ToTarget = rot_apply(ROT_A_CURRENT_COMP, ToTarget);
coords_get(ToTarget, &tx, &ty, &tz);
dist=sqrt(tx*tx+ty*ty+tz*tz);
} else if (dist) {
tx = 0;
ty = 0;
tz = dist;
}
if (srcArea <= 0) {
printf("Source_simple: %s: Source area is <= 0 !\n ERROR - Exiting\n",
NAME_CURRENT_COMP);
exit(0);
}
if (dist <= 0 || focus_xw <= 0 || focus_yh <= 0) {
printf("Source_simple: %s: Target area unmeaningful! (negative dist / focus_xw / focus_yh)\n ERROR - Exiting\n",
NAME_CURRENT_COMP);
exit(0);
}
if ((!lambda0 && !E0 && !dE && !dlambda)) {
printf("Source_simple: %s: You must specify either a wavelength or energy range!\n ERROR - Exiting\n",
NAME_CURRENT_COMP);
exit(0);
}
if ((!lambda0 && !dlambda && (E0 <= 0 || dE < 0 || E0-dE <= 0))
|| (!E0 && !dE && (lambda0 <= 0 || dlambda < 0 || lambda0-dlambda <= 0))) {
printf("Source_simple: %s: Unmeaningful definition of wavelength or energy range!\n ERROR - Exiting\n",
NAME_CURRENT_COMP);
exit(0);
}
%}
TRACE
%{
double chi,E,lambda,v,r, xf, yf, rf, dx, dy, pdir;
t=0;
z=0;
if (square == 1) {
x = xwidth * (rand01() - 0.5);
y = yheight * (rand01() - 0.5);
} else {
chi=2*PI*rand01(); /* Choose point on source */
r=sqrt(rand01())*radius; /* with uniform distribution. */
x=r*cos(chi);
y=r*sin(chi);
}
randvec_target_rect_real(&xf, &yf, &rf, &pdir,
tx, ty, tz, focus_xw, focus_yh, ROT_A_CURRENT_COMP, x, y, z, 2);
dx = xf-x;
dy = yf-y;
rf = sqrt(dx*dx+dy*dy+rf*rf);
p = pdir*pmul;
if(lambda0==0) {
if (!gauss) {
E=E0+dE*randpm1(); /* Choose from uniform distribution */
} else {
E=E0+randnorm()*dE;
}
v=sqrt(E)*SE2V;
} else {
if (!gauss) {
lambda=lambda0+dlambda*randpm1();
} else {
lambda=lambda0+randnorm()*dlambda;
}
v = K2V*(2*PI/lambda);
}
vz=v*dist/rf;
vy=v*dy/rf;
vx=v*dx/rf;
%}
MCDISPLAY
%{
if (square == 1) {
rectangle("xy",0,0,0,xwidth,yheight);
} else {
circle("xy",0,0,0,radius);
}
if (dist) {
dashed_line(0,0,0, -focus_xw/2+tx,-focus_yh/2+ty,tz, 4);
dashed_line(0,0,0, focus_xw/2+tx,-focus_yh/2+ty,tz, 4);
dashed_line(0,0,0, focus_xw/2+tx, focus_yh/2+ty,tz, 4);
dashed_line(0,0,0, -focus_xw/2+tx, focus_yh/2+ty,tz, 4);
}
%}
END
|
