1/*******************************************************************************
2*
3* McStas, neutron ray-tracing package
4*         Copyright (C) 1997-2006, All rights reserved
5*         Risoe National Laboratory, Roskilde, Denmark
6*         Institut Laue Langevin, Grenoble, France
7*
8* Component: Exact_radial_coll
9*
10* %I
11* Written by: Roland Schedler <roland.schedler at hmi.de>
12* Modified by: using Collimator_radial Component
13* Date: October 2006
14* Origin: HMI
15* Modified by: E. Farhi, uniformize parameter names (Jul 2008)
16*
17* An exact radial Soller collimator.
18*
19* %D
20* Radial Soller collimator with rectangular opening, specified length and
21* specified foil thickness.
22* The collimator is made of many trapezium shaped nslit stacked radially.
23* The nslit are separated by absorbing foils, the whole stuff is inside
24* an absorbing housing.
25* The component should be positioned at the radius center. The model is exact.
26* The neutron beam outside the collimator area is transmitted unaffected.
27*
28* Example: Exact_radial_coll(theta_min=-5, theta_max=5, nslit=100,
29*          radius=1.0, length=.3, h_in=.2, h_out=.3, d=0.0001)
30*
31*
32* %P
33* INPUT PARAMETERS:
34*
35* theta_min: [deg]  Minimum Theta angle for the radial setting
36* theta_max: [deg]  Maximum Theta angle for the radial setting
37* nslit: [1]        Number of channels in the theta range
38* radius: [m]       Inner radius (focus point to foil start point).
39* length: [m]       Length of the foils / collimator
40* h_in: [m]         Input  window height
41* h_out: [m]        Output window height
42* d: [m]            Thickness of the absorbing foils
43* verbose: [0/1]    Gives additional information
44*
45* %E
46*******************************************************************************/
47
48
49DEFINE COMPONENT Exact_radial_coll
50DEFINITION PARAMETERS ()
51SETTING PARAMETERS (theta_min=-5, theta_max=5, nslit=100,
52radius=1.0, length=.5, h_in=.3, h_out=.4,
53d=0.0001, verbose=0)
54OUTPUT PARAMETERS (alpha_in, alpha_out, out_radius, beta_in, beta_out, iw, ow,divergence,theta)
55/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
56DECLARE
57%{
58double alpha_in, alpha_out, beta_in, beta_out, theta;
59double out_radius, iw, ow, divergence;
60%}
61INITIALIZE
62%{
63/* check for input parameters */
64if (radius <= 0) exit(printf("Exact_radial_coll: %s: radius must be positive\n", NAME_CURRENT_COMP));
65  if (h_in <= 0) exit(printf("Exact_radial_coll: %s: h_in must be positive\n", NAME_CURRENT_COMP));
66  if (h_out <= 0) exit(printf("Exact_radial_coll: %s: h_out must be positive\n", NAME_CURRENT_COMP));
67  if (d <= 0) exit(printf("Exact_radial_coll: %s: d must be positive\n", NAME_CURRENT_COMP));
68  if (nslit <= 0)  exit(printf("Exact_radial_coll: %s: number of channels must be positive\n", NAME_CURRENT_COMP));
69  if ((nslit - floor (nslit)) > 0) exit(printf("Exact_radial_coll: %s: number of channels must be an integer\n", NAME_CURRENT_COMP));
70  if (length <= 0)    exit(printf("Exact_radial_coll: %s: collimator length must be positive\n", NAME_CURRENT_COMP));
71  if (theta_max <= theta_min) exit(printf("Exact_radial_coll: %s: theta_max must be greater than theta_min\n", NAME_CURRENT_COMP));
72
73  theta_max *= DEG2RAD;
74  theta_min *= DEG2RAD;
75  theta = theta_max - theta_min;
76  out_radius = radius + length;
77  beta_in =  2*asin(d / (2 * radius));
78  beta_out=  2*asin(d / (2 * out_radius));
79  if (theta < nslit*beta_in) exit(printf("Exact_radial_coll: %s: the %6.0f foils of %g [meter]\n"
80                                             "do not fit within the angular range theta = %4.2f [deg]\n",
81                                             NAME_CURRENT_COMP, nslit, d, theta*RAD2DEG));
82  alpha_in = (theta - nslit*beta_in)/nslit;
83  alpha_out = (theta - nslit*beta_out)/nslit;
84  iw = 2*radius*sin((alpha_in/2));
85  ow = 2*out_radius*sin((alpha_out/2));
86  divergence=(iw+ow)/(sqrt(4*length*length-(ow-iw)*(ow-iw)));
87
88  if (verbose) {
89    printf("Exact_radial_coll: %s: foil thickness is %.2g [millimeter]\n", NAME_CURRENT_COMP, d*1000);
90    printf("                  opening each  input slit [%.3g:%.0f] [millimeter]\n", iw*1000, h_in*1000);
91    printf("                  opening each output slit [%.3g:%.0f] [millimeter]\n", ow*1000, h_out*1000);
92    printf("                  divergence per channel is %g [min] \n", divergence*RAD2MIN);
93  }
94%}
95
96TRACE
97%{
98  double phi, t0, t1, t2, t3;
99  int    intersect;
100  long   input_chan,  output_chan;
101  double input_theta, output_theta;
102  double input_center,output_center;
103  double window_theta;
104  char   ok=0;
105
106  /* first compute intersection time with input cylinder */
107  intersect=cylinder_intersect(&t0,&t3,x,y,z,vx,vy,vz,radius,h_in);
108  if (!intersect) ABSORB;
109  else if (t3 > t0) t0 = t3;
110
111  intersect=cylinder_intersect(&t1,&t2,x,y,z,vx,vy,vz,out_radius,h_out);
112  if (!intersect) ABSORB;
113  else if (t2 > t1) t1 = t2;
114
115  /* get index of input slit */
116  if (t0 > 0 && t1 > t0) {
117      PROP_DT(t0);
118      input_theta = atan2(x, z);
119     /* channel number (start at 0) */
120      window_theta = (theta_max - theta_min)/nslit;
121      input_chan  = floor((input_theta - theta_min)/window_theta);
122      if (input_chan >= 0 && input_chan < nslit && fabs(y) < h_in/2) ok=1;
123    if (ok) {
124        input_center= theta_min + input_chan*window_theta + (window_theta)/2;
125        /* are we outside the soller or in the foil? */
126        phi = input_theta - input_center;
127        if (fabs(phi) > alpha_in/2) ABSORB; /* inside the foil*/
128        SCATTER;
129
130      /* propagate to output radius */
131      PROP_DT(t1-t0);
132      SCATTER;
133        output_theta = atan2(x, z);
134        /* channel number (start at 0) */
135        output_chan  = floor((output_theta - theta_min)/window_theta);
136        /* did we change channel ? */
137        if (output_chan != input_chan) ABSORB; /* changed slit */
138        output_center= theta_min + output_chan*window_theta
139                    + (window_theta)/2;
140        /* are we outside the soller */
141        phi = output_theta -output_center;
142        if (fabs(phi) > alpha_out/2 || fabs(y) > h_out/2) ABSORB; /* outside output slit */
143
144    } /* else we pass aside the entrance window of radial collimator */
145    else {
146      /* propagate to output radius */
147      PROP_DT(t1-t0);
148      SCATTER;
149        output_theta = atan2(x, z);
150        /* channel number (start at 0) */
151        output_chan  = floor((output_theta - theta_min)/window_theta);
152       /* are we come from outside into the soller or in the foil?*/
153        if (output_chan >= 0 || output_chan < nslit) ABSORB;
154    } /* else we pass aside the exit window of radial collimator */
155  }   /* else did not encounter collimator */
156
157%}
158
159MCDISPLAY
160%{
161  int i;
162  double theta1, theta2, theta3, theta4;
163  double x_in_l,  z_in_l,  x_in_r,  z_in_r;
164  double x_out_l, z_out_l, x_out_r, z_out_r;
165  double window_theta, y1, y2;
166
167
168  window_theta = alpha_in + beta_in;
169  y1 = h_in/2;
170  y2 = h_out/2;
171
172  theta1 = theta_min;
173  theta3 = theta1+beta_in/2;
174  theta4 = theta1+beta_out/2;
175
176  z_in_l = radius*cos(theta1);
177  x_in_l = radius*sin(theta1);
178  z_in_r = radius*cos(theta3);
179  x_in_r = radius*sin(theta3);
180
181  z_out_l = out_radius*cos(theta1);
182  x_out_l = out_radius*sin(theta1);
183  z_out_r = out_radius*cos(theta4);
184  x_out_r = out_radius*sin(theta4);
185
186    multiline(5,
187      x_in_l, -y1, z_in_l,
188      x_in_l,  y1, z_in_l,
189      x_out_l, y2, z_out_l,
190      x_out_l,-y2, z_out_l,
191      x_in_l, -y1, z_in_l);
192
193   line(x_in_l,   y1, z_in_l,  x_in_r,  y1, z_in_r);
194   line(x_in_l,  -y1, z_in_l,  x_in_r, -y1, z_in_r);
195   line(x_out_l,  y2, z_out_l, x_out_r, y2, z_out_r);
196   line(x_out_l, -y2, z_out_l, x_out_r,-y2, z_out_r);
197
198   multiline(5,
199      x_in_r, -y1, z_in_r,
200      x_in_r,  y1, z_in_r,
201      x_out_r, y2, z_out_r,
202      x_out_r,-y2, z_out_r,
203      x_in_r, -y1, z_in_r);
204
205  for (i = 1; i < nslit; i++) {
206    theta1 = i*window_theta+theta_min-beta_in/2;
207    theta2 = i*window_theta+theta_min+beta_in/2;
208    theta3 = i*window_theta+theta_min-beta_out/2;
209    theta4 = i*window_theta+theta_min+beta_out/2;
210
211    z_in_l = radius*cos(theta1);
212    x_in_l = radius*sin(theta1);
213    z_in_r = radius*cos(theta2);
214    x_in_r = radius*sin(theta2);
215
216    z_out_l = out_radius*cos(theta3);
217    x_out_l = out_radius*sin(theta3);
218    z_out_r = out_radius*cos(theta4);
219    x_out_r = out_radius*sin(theta4);
220    /* left side */
221    multiline(5,
222      x_in_l, -y1, z_in_l,
223      x_in_l,  y1, z_in_l,
224      x_out_l, y2, z_out_l,
225      x_out_l,-y2, z_out_l,
226      x_in_l, -y1, z_in_l);
227   /* left -> right lines */
228   line(x_in_l,   y1, z_in_l,  x_in_r,  y1, z_in_r);
229   line(x_in_l,  -y1, z_in_l,  x_in_r, -y1, z_in_r);
230   line(x_out_l,  y2, z_out_l, x_out_r, y2, z_out_r);
231   line(x_out_l, -y2, z_out_l, x_out_r,-y2, z_out_r);
232   /* right side */
233   multiline(5,
234      x_in_r, -y1, z_in_r,
235      x_in_r,  y1, z_in_r,
236      x_out_r, y2, z_out_r,
237      x_out_r,-y2, z_out_r,
238      x_in_r, -y1, z_in_r);
239  }
240
241  /* remaining bits */
242
243  theta1 = theta_max;
244  theta3 = theta1-beta_in/2;
245  theta4 = theta1-beta_out/2;
246
247  z_in_l = radius*cos(theta1);
248  x_in_l = radius*sin(theta1);
249  z_in_r = radius*cos(theta3);
250  x_in_r = radius*sin(theta3);
251
252  z_out_l = out_radius*cos(theta1);
253  x_out_l = out_radius*sin(theta1);
254  z_out_r = out_radius*cos(theta4);
255  x_out_r = out_radius*sin(theta4);
256
257    multiline(5,
258      x_in_l, -y1, z_in_l,
259      x_in_l,  y1, z_in_l,
260      x_out_l, y2, z_out_l,
261      x_out_l,-y2, z_out_l,
262      x_in_l, -y1, z_in_l);
263
264   line(x_in_l,   y1, z_in_l,  x_in_r,  y1, z_in_r);
265   line(x_in_l,  -y1, z_in_l,  x_in_r, -y1, z_in_r);
266   line(x_out_l,  y2, z_out_l, x_out_r, y2, z_out_r);
267   line(x_out_l, -y2, z_out_l, x_out_r,-y2, z_out_r);
268
269   multiline(5,
270      x_in_r, -y1, z_in_r,
271      x_in_r,  y1, z_in_r,
272      x_out_r, y2, z_out_r,
273      x_out_r,-y2, z_out_r,
274      x_in_r, -y1, z_in_r);
275
276%}
277
278END
279