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
|
/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright 1997-2002, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: Source_div
*
* %I
* Written by: KL
* Date: November 20, 1998
* Modified by: KL, 8 October 2001
* Origin: Risoe
*
* Neutron source with Gaussian or uniform divergence
*
* %D
* The routine is a rectangular neutron source, which has a gaussian or uniform
* divergent output in the forward direction.
* The neutron energy is 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).
* In the case of uniform distribution (gauss=0), angles are uniformly distributed
* between -focus_aw and +focus_aw as well as -focus_ah and +focus_ah.
* For Gaussian distribution (gauss=1), 'focus_aw' and 'focus_ah' define the
* FWHM of a Gaussian distribution. Energy/wavelength distribution is also
* Gaussian.
*
* Example: Source_div(xwidth=0.1, yheight=0.1, focus_aw=2, focus_ah=2, E0=14, dE=2, gauss=0)
*
* %VALIDATION
* Feb 2005: tested by Kim Lefmann (o.k.)
* Apr 2005: energy distribution used in external tests of Fermi choppers (o.k.)
* Jun 2005: wavelength distribution used in external tests of velocity selectors (o.k.)
* Validated by: K. Lieutenant
*
* %BUGS
* distribution is uniform in (hor. and vert.) angle (relative to moderator normal),
* therefore not suited for large angles
*
* %P
* xwidth: [m] Width of source
* yheight: [m] Height of source
* focus_aw: [deg] FWHM (Gaussian) or maximal (uniform) horz. width divergence
* focus_ah: [deg] FWHM (Gaussian) or maximal (uniform) vert. height divergence
* E0: [meV] Mean energy of neutrons.
* dE: [meV] Energy half spread of neutrons.
* lambda0: [Ang] Mean wavelength of neutrons (only relevant for E0=0)
* dlambda: [Ang] Wavelength half spread of neutrons.
* gauss: [0|1] Criterion: 0: uniform, 1: Gaussian distributions
* flux: [1/(s cm 2 st energy_unit)] flux per energy unit, Angs or meV
*
* CALCULATED PARAMETERS:
* sigmah: [rad] parameter 'sigma' of the Gaussian distribution for horizontal divergence
* sigmav: [rad] parameter 'sigma' of the Gaussian distribution for vertical divergence
* p_init: [1] normalisation factor 1/'neutron_count'
*
* %E
*******************************************************************************/
DEFINE COMPONENT Source_div
SETTING PARAMETERS (xwidth, yheight, focus_aw, focus_ah, E0=0.0, dE=0.0, lambda0=0.0, dlambda=0.0, gauss=0, flux=1)
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
DECLARE
%{
double sigmah;
double sigmav;
double p_init;
double dist;
double focus_xw;
double focus_yh;
%}
INITIALIZE
%{
sigmah = DEG2RAD*focus_aw/(sqrt(8.0*log(2.0)));
sigmav = DEG2RAD*focus_ah/(sqrt(8.0*log(2.0)));
if (xwidth < 0 || yheight < 0 || focus_aw < 0 || focus_ah < 0) {
printf("Source_div: %s: Error in input parameter values!\n"
"ERROR Exiting\n",
NAME_CURRENT_COMP);
exit(0);
}
if ((!lambda0 && !E0 && !dE && !dlambda)) {
printf("Source_div: %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_div: %s: Unmeaningful definition of wavelength or energy range!\n ERROR - Exiting\n",
NAME_CURRENT_COMP);
exit(0);
}
/* compute distance to next component */
Coords ToTarget;
double tx,ty,tz;
ToTarget = coords_sub(POS_A_COMP_INDEX(INDEX_CURRENT_COMP+1),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);
/* compute target area */
if (dist) {
focus_xw=dist*tan(focus_aw*DEG2RAD);
focus_yh=dist*tan(focus_ah*DEG2RAD);
}
p_init = flux*1e4*xwidth*yheight/mcget_ncount();
if (!focus_aw || !focus_ah)
exit(printf("Source_div: %s: Zero divergence defined. \n"
"ERROR Use non zero values for focus_aw and focus_ah.\n",
NAME_CURRENT_COMP));
p_init *= 2*fabs(DEG2RAD*focus_aw*sin(DEG2RAD*focus_ah/2)); /* solid angle */
if (dlambda)
p_init *= 2*dlambda;
else if (dE)
p_init *= 2*dE;
%}
TRACE
%{
double E,lambda,v;
double tan_h;
double tan_v;
double thetah;
double thetav;
p=p_init;
z=0;
t=0;
x=randpm1()*xwidth/2.0;
y=randpm1()*yheight/2.0;
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);
}
if (gauss==1) {
thetah = randnorm()*sigmah;
thetav = randnorm()*sigmav;
} else {
thetah = randpm1()*focus_aw*DEG2RAD/2;
thetav = randpm1()*focus_ah*DEG2RAD/2;
}
tan_h = tan(thetah);
tan_v = tan(thetav);
/* Perform the correct treatment - no small angle approx. here! */
vz = v / sqrt(1 + tan_v*tan_v + tan_h*tan_h);
vy = tan_v * vz;
vx = tan_h * vz;
%}
MCDISPLAY
%{
multiline(5, -xwidth/2.0, -yheight/2.0, 0.0,
xwidth/2.0, -yheight/2.0, 0.0,
xwidth/2.0, yheight/2.0, 0.0,
-xwidth/2.0, yheight/2.0, 0.0,
-xwidth/2.0, -yheight/2.0, 0.0);
if (dist) {
dashed_line(0,0,0, -focus_xw/2,-focus_yh/2,dist, 4);
dashed_line(0,0,0, focus_xw/2,-focus_yh/2,dist, 4);
dashed_line(0,0,0, focus_xw/2, focus_yh/2,dist, 4);
dashed_line(0,0,0, -focus_xw/2, focus_yh/2,dist, 4);
}
%}
END
|
