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
|
/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright 1997-2003, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: Guide_curved
*
* %I
* Written by: Ross Stewart
* Date: November 18 2003
* Origin: <a href="http://www.ill.fr">ILL (France)</a>.
* Modified by: E. Farhi, uniformize parameter names (Jul 2008)
*
* Non-focusing curved neutron guide.
*
* %D
* Models a rectangular curved guide tube with entrance centered on the Z axis.
* The entrance lies in the X-Y plane. Draws a true depiction
* of the guide, and trajectories. Guide is not focusing.
*
* Example: Guide_curved(w1=0.1, h1=0.1, l=2.0, R0=0.99, Qc=0.021,
* alpha=6.07, m=2, W=0.003, curvature=2700)
*
* %BUGS
* This component does not work with gravitation on. Use component Guide_gravity then.
* Systematic error on transmitted flux is found to be about 10%.
*
* %P
* INPUT PARAMETERS:
*
* w1: [m] Width at the guide entry
* h1: [m] Height at the guide entry
* l: [m] length of guide
* 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
* curvature: [m] Radius of curvature of the guide
*
* %L
* <a href="../optics/Bender.comp.html">Bender</a>
*
* %E
*******************************************************************************/
DEFINE COMPONENT Guide_curved
SETTING PARAMETERS (w1, h1, l, R0=0.995, Qc=0.0218, alpha=4.38, m=2, W=0.003, curvature=2700)
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
SHARE
%{
%include "ref-lib"
%}
INITIALIZE
%{
if (mcgravitation) fprintf(stderr,"WARNING: Guide_curved: %s: "
"This component produces wrong results with gravitation !\n"
"Use Guide_gravity.\n",
NAME_CURRENT_COMP);
%}
TRACE
%{
double t11, t12, t21, t22, theta, alphaAng, endtime, phi;
double time, time1, time2, q, R;
int ii, i_bounce;
double whalf = 0.5*w1, hhalf = 0.5*h1; /* half width and height of guide */
double z_off = curvature*sin(l/curvature); /* z-component of total guide length */
double R1 = curvature - whalf; /* radius of curvature of inside mirror */
double R2 = curvature + whalf; /* radius of curvature of outside mirror */
double vel = sqrt(vx*vx + vy*vy + vz*vz); /* neutron velocity */
double vel_xz = sqrt(vx*vx + vz*vz); /* in plane velocity */
double K = V2K*vel; /* neutron wavevector */
double lambda = 2.0*PI/K; /* neutron wavelength */
/* Propagate neutron to guide entrance. */
PROP_Z0;
if(x <= -whalf || x >= whalf || y <= -hhalf || y >= hhalf)
ABSORB;
SCATTER;
for(;;)
{
double par[]={R0, Qc, alpha, m, W};
/* Find itersection points of neutron with inside and outside guide walls */
ii = cylinder_intersect(&t11, &t12 ,x - curvature, y, z, vx, vy, vz, R1, h1);
ii = cylinder_intersect(&t21, &t22 ,x - curvature, y, z, vx, vy, vz, R2, h1);
/* Choose appropriate reflection time */
time1 = (t11 < 1e-7) ? t12 : t11;
time2 = (t21 < 1e-7) ? t22 : t21;
time = (time1 < 1e-7 || time2 < time1) ? time2 : time1;
/* Has neutron left the guide? */
endtime = (z_off - z)/vz;
if (time > endtime || time <= 1e-7) break;
PROP_DT(time);
/* Find reflection surface */
R = (time == time1) ? R1 : R2;
i_bounce = (fabs(y - hhalf) < 1e-7 || fabs(y + hhalf) < 1e-7) ? 2 : 1;
switch(i_bounce) {
case 1: /* Inside or Outside wall */
phi = atan(vx/vz); /* angle of neutron trajectory */
alphaAng = asin(z/R); /* angle of guide wall */
theta = fabs(phi - alphaAng); /* angle of reflection */
vz = vel_xz*cos(2.0*alphaAng - phi);
vx = vel_xz*sin(2.0*alphaAng - phi);
break;
case 2: /* Top or Bottom wall */
theta = fabs(atan(vy/vz));
vy = -vy;
break;
}
/* Now compute reflectivity. */
if (m == 0 || !R0) ABSORB;
q = 4.0*PI*sin(theta)/lambda;
StdReflecFunc(q, par, &R);
if (R >= 0) p *= R; else ABSORB;
SCATTER;
}
%}
MCDISPLAY
%{
double x1, x2, z1, z2;
double xplot1[100], xplot2[100], zplot1[100], zplot2[100];
int n = 100;
int j = 1;
double R1 = (curvature - 0.5*w1); /* radius of inside arc */
double R2 = (curvature + 0.5*w1); /* radius of outside arc */
for(j=0;j<n;j++) {
z1 = ((double)j)*(R1*l/curvature)/(double)(n - 1);
z2 = ((double)j)*(R2*l/curvature)/(double)(n - 1);
x1 = curvature - sqrt(R1*R1 - z1*z1);
x2 = curvature - sqrt(R2*R2 - z2*z2);
xplot1[j] = x1;
xplot2[j] = x2;
zplot1[j] = z1;
zplot2[j] = z2;
}
line(xplot1[0], 0.5*h1,zplot1[0],xplot2[0], 0.5*h1,zplot2[0]);
line(xplot1[0], 0.5*h1,zplot1[0],xplot1[0],-0.5*h1,zplot1[0]);
line(xplot2[0],-0.5*h1,zplot2[0],xplot2[0], 0.5*h1,zplot2[0]);
line(xplot1[0],-0.5*h1,zplot1[0],xplot2[0],-0.5*h1,zplot2[0]);
for(j=0;j<n-1;j++) {
line(xplot1[j], 0.5*h1, zplot1[j], xplot1[j+1], 0.5*h1, zplot1[j+1]);
line(xplot2[j], 0.5*h1, zplot2[j], xplot2[j+1], 0.5*h1, zplot2[j+1]);
line(xplot1[j], -0.5*h1, zplot1[j], xplot1[j+1], -0.5*h1, zplot1[j+1]);
line(xplot2[j], -0.5*h1, zplot2[j], xplot2[j+1], -0.5*h1, zplot2[j+1]);
}
line(xplot1[n-1], 0.5*h1,zplot1[n-1],xplot2[n-1], 0.5*h1,zplot2[n-1]);
line(xplot1[n-1], 0.5*h1,zplot1[n-1],xplot1[n-1],-0.5*h1,zplot1[n-1]);
line(xplot2[n-1],-0.5*h1,zplot2[n-1],xplot2[n-1], 0.5*h1,zplot2[n-1]);
line(xplot1[n-1],-0.5*h1,zplot1[n-1],xplot2[n-1],-0.5*h1,zplot2[n-1]);
%}
END
|
