Ticket #7: Collimator_radial.comp

File Collimator_radial.comp, 10.3 KB (added by farhi, 2 years ago)

Fixed Collimator_radial

Line 
1/*******************************************************************************
2*
3* McStas, neutron ray-tracing package
4*         Copyright 1997-2002, All rights reserved
5*         Risoe National Laboratory, Roskilde, Denmark
6*         Institut Laue Langevin, Grenoble, France
7*
8* Component: Collimator_radial
9*
10* %I
11* Written by: Emmanuel Farhi <farhi@ill.fr>
12* Date: July 2005
13* Version: $Revision: 1.11 $
14* Origin: ILL
15* Release: McStas 1.9
16*
17* A radial Soller collimator.
18*
19* %D
20* Radial Soller collimator with rectangular opening and specified length.
21* The collimator is made of many rectangular channels stacked radially.
22* Each channel is a set of transmitting layers (nslit), separated by an absorbing
23* material (infinitely thin), the whole stuff is inside an absorbing housing.
24* When using zero as the number of channels (nchan), the collimator is continuous,
25* whithout shadowing effect.
26* The component should be positioned at the radius center.
27* The component can be made oscillating (usual on diffractometers and TOF
28* machines) with the 'roc' parameter.
29* The neutron beam outside the collimator angular area is transmitted unaffected.
30*
31* Example: Channelled radial collimator with shadow parts
32*        Collimator_radial(xwidth=0.015, yheight=.3, length=.35,
33*                   divergence=40,transmission=1,
34*                   theta_min=5, theta_max=165, nchan=128, radius=0.9)
35*   A continuous radial collimator
36*        Collimator_radial(yheight=.3, length=.35,
37*                   divergence=40,transmission=1,
38*                   theta_min=5, theta_max=165, radius=0.9)
39*
40* %P
41* INPUT PARAMETERS:
42*
43* xwidth:     (m)          Soller  window width, filled with nslit slits.
44*                            Use 0 value for continuous collimator.
45* yheight:    (m)          Collimator height.
46* length:     (m)          Length/Distance between inner and outer slits.
47* divergence: (min of arc) Divergence angle. May also be specified with the
48*                            nslit parameter
49* theta_min:  (deg)        Minimum Theta angle for the radial setting.
50* theta_max:  (deg)        Maximum Theta angle for the radial setting.
51* nchan:      (1)          Number of Soller channels in the theta range.
52*                            Use 0 value for continuous collimator.
53* radius:     (m)          Radius of the collimator (to entry window).
54*
55* Optional parameters
56* transmission:(1)         Transmission of Soller (0<=t<=1).
57* nslit:       (1)         Number of blades composing each Soller. Overrides
58*                            the divergence parameter.
59* roc:         (deg)       Amplitude of oscillation of collimator. 0=fixed.
60* verbose:     (0/1)       Gives additional information.
61* approx:      (0/1)       Use Soller triangular transmission approximation.
62*
63* %E
64*******************************************************************************/
65DEFINE COMPONENT Collimator_radial
66DEFINITION PARAMETERS ()
67SETTING PARAMETERS (xwidth=0, yheight=.3, length=.35,
68                    divergence=0,transmission=1,
69                    theta_min=5, theta_max=165, nchan=0, radius=1.3, nslit=0,
70                    roc=0, verbose=0, approx=0)
71OUTPUT PARAMETERS ()
72STATE PARAMETERS (x,y,z,vx,vy,vz,t,sx,sy,sz,p)
73DECLARE
74%{
75  double width_of_slit=0, width_of_Soller=0, slit_theta=0;
76%}
77INITIALIZE
78%{
79  if (radius <= 0)
80    exit(printf("Collimator_radial: %s: incorrect radius=%g\n",
81      NAME_CURRENT_COMP, radius));
82  if (length <= 0)
83    exit(printf("Collimator_radial: %s: invalid collimator length=%g\n",
84      NAME_CURRENT_COMP, length));
85  if (transmission <= 0 || transmission >1)
86    exit(printf("Collimator_radial: %s: invalid transmission=%g\n",
87      NAME_CURRENT_COMP, transmission));
88 
89  theta_max *= DEG2RAD;
90  theta_min *= DEG2RAD;
91  roc       *= DEG2RAD;
92  divergence*= MIN2RAD;
93 
94  if (xwidth && !nchan)
95    nchan  = ceil(radius*fabs(theta_max-theta_min)/xwidth);
96  else if (!xwidth && nchan)
97    xwidth = radius*fabs(theta_max-theta_min)/nchan;
98 
99  /* determine total width [m] of Soller channels, containing nslit in xwidth */
100  if (nchan) width_of_Soller = radius*fabs(theta_max-theta_min)/nchan;
101  else       width_of_Soller = 0;
102 
103  if (!nchan || !xwidth || xwidth > width_of_Soller)
104    nchan=xwidth=width_of_Soller=0; /* continuous collimator */
105 
106  /* determine width [m] of slits */
107  if (divergence) {
108    width_of_slit = length*tan(divergence);
109    if (xwidth)           /* Soller */
110      nslit = ceil(xwidth/width_of_slit);
111    else if (!nchan)      /* continuous collimator */
112      nslit = ceil(radius*fabs(theta_max-theta_min)/width_of_slit);
113  } else {
114    if (!nchan && nslit)  /* continuous collimator */
115      width_of_slit = radius*fabs(theta_max-theta_min)/nslit;
116    else if (nchan && nslit)  /* Soller */
117      width_of_slit = xwidth/nslit;
118    divergence = atan2(width_of_slit,length);
119  }
120 
121  if (nslit <= 0)
122    exit(printf("Collimator_radial: %s: number of channels must be positive nslit=%g.\n"
123                "                   Specify alternatively divergence, xwidth and nchan.\n",
124                NAME_CURRENT_COMP, nslit));
125 
126  if (verbose) {
127    printf("Collimator_radial: %s: divergence %g [min] %s"
128           ". Total opening [%g:%g] [deg]\n",
129           NAME_CURRENT_COMP, divergence*RAD2MIN,
130           (roc ? "oscillating" : ""),
131           theta_min*RAD2DEG, theta_max*RAD2DEG);
132    if (approx)
133      printf("    Using triangular approximation model");
134    else if (!nchan)
135      printf("    Using continuous model");
136    else
137      printf("    Using %i Soller channels of width %g [cm]",
138        (int)nchan, width_of_Soller*100);
139
140    printf(" with %i slits of width %g [mm] pitch %g [deg].\n",
141      (int)nslit, width_of_slit*1000, atan2(width_of_slit, radius)*RAD2DEG);
142  }
143
144%}
145
146TRACE
147%{
148  double intersect=0;
149  double t0, t1, t2, t3;
150  /* determine intersection with inner and outer cylinders */
151  intersect=cylinder_intersect(&t0,&t3,x,y,z,vx,vy,vz,radius,yheight);
152  if (!intersect) ABSORB;
153  else if (t3 > t0) t0 = t3;
154
155  intersect=cylinder_intersect(&t1,&t2,x,y,z,vx,vy,vz,radius+length,yheight);
156  if (!intersect) ABSORB;
157  else if (t2 > t1) t1 = t2;
158 
159  /* propagate/determine neutron position at inner cylinder */
160  if (t0 > 0 && t1 > t0) {
161    double input_chan=0;
162    double input_theta=0, output_theta=0;
163    double roc_theta=0;
164    double input_slit=0,  output_slit=0;
165   
166    PROP_DT(t0);
167 
168    /* apply ROC oscillation with a linear random distribution */
169    if (roc) roc_theta = roc*randpm1()/2; else roc_theta=0;
170   
171    /* angle on inner cylinder */
172    input_theta = atan2(x, z) + roc_theta; 
173   
174    /* check if we are within min/max collimator input bounds */
175    if (input_theta >= theta_min && input_theta <= theta_max) {
176      SCATTER;
177      /* input Soller channel index */
178      if (width_of_Soller) {
179        input_chan = radius*(input_theta-theta_min)/width_of_Soller;
180        input_chan = input_chan-floor(input_chan); /* position within Soller [0:1] */
181       
182        /* check if we hit an absorbing housing (between Sollers): ABSORB
183         *   total Soller aperture is width_of_Soller,
184         *   containg a slit pack  of aperture xwidth
185         */
186        if (input_chan < (1-xwidth/width_of_Soller)/2
187         || input_chan > (1+xwidth/width_of_Soller)/2) ABSORB;
188      }
189
190      /* determine input slit index */
191      input_slit = floor(input_theta*radius/width_of_slit);
192     
193    } else /* neutron missed collimator input range */
194      input_theta=4*PI;
195 
196    /* propagate to outer cylinder */
197    PROP_DT(t1-t0);
198   
199    /* angle on outer cylinder */
200    output_theta = atan2(x, z) + roc_theta;
201   
202    /* check if we are within min/max collimator output bounds */
203    if (output_theta >= theta_min && output_theta <= theta_max) {
204      /* check if we come from sides:   ABSORB */
205      if (input_theta > 2*PI) ABSORB; /* input_theta=4*PI when missed input */
206     
207      if (approx) {
208        double phi=atan2(x, z)-atan2(vx, vz); /* difference between positional slit angle and velocity */
209        if (fabs(phi) > divergence)
210          ABSORB; /* get outside transmission */
211        else
212          p *= (1.0 - phi/divergence);
213      } else {
214        /* check if we have changed slit: ABSORB */
215        /* slits are considered radial so that their output size is:
216           width_of_slit*(radius+length)/radius and it turns out this the same exp as for input
217         */
218        output_slit = floor(output_theta*radius/width_of_slit);
219        if (input_slit != output_slit) ABSORB;
220      }
221      SCATTER;
222      p *= transmission;
223     
224    } /* else neutron missed collimator output range */
225 
226  } /* else did not encounter collimator cylinders */
227
228%}
229
230MCDISPLAY
231%{
232  double Soller_theta;
233  double height=yheight/2;
234  double theta1, theta2;
235  double x_in_l,  z_in_l,  x_in_r,  z_in_r;
236  double x_out_l, z_out_l, x_out_r, z_out_r;
237  int i;
238 
239  /* display collimator radial geometry:
240     in order to avoid too many lines, we shown main housing and channels
241     but no slit */
242     
243  if (!nchan || nchan > 20)     nchan=20;
244  if (nchan > 64)               nchan=64;
245  Soller_theta=fabs(theta_max-theta_min)/nchan; /* angular width of Soller */
246 
247  /* draw all channels, which also show housing */
248  magnify("xy");
249  for (i = 0; i < nchan; i++) {
250   
251    theta1 = i*Soller_theta+theta_min;
252    theta2 = theta1+Soller_theta;
253
254    z_in_l = radius*cos(theta1);
255    x_in_l = radius*sin(theta1);
256    z_in_r = radius*cos(theta2);
257    x_in_r = radius*sin(theta2);
258
259    z_out_l = (radius+length)*cos(theta1);
260    x_out_l = (radius+length)*sin(theta1);
261    z_out_r = (radius+length)*cos(theta2);
262    x_out_r = (radius+length)*sin(theta2);
263    /* left side */
264    multiline(6,
265      x_in_l, -height, z_in_l,
266      x_in_l,  height, z_in_l,
267      x_out_l, height, z_out_l,
268      x_out_l,-height, z_out_l,
269      x_in_l, -height, z_in_l,
270      x_in_r, -height, z_in_r);
271   /* left -> right lines */
272   line(x_in_l,   height, z_in_l,  x_in_r,  height, z_in_r);
273   line(x_out_l,  height, z_out_l, x_out_r, height, z_out_r);
274   line(x_out_l, -height, z_out_l, x_out_r,-height, z_out_r);
275  }
276  /* remaining bits */
277  theta1 = nchan*Soller_theta+theta_min;
278  z_in_l = radius*cos(theta1);
279  x_in_l = radius*sin(theta1);
280  z_out_l = (radius+length)*cos(theta1);
281  x_out_l = (radius+length)*sin(theta1);
282  multiline(5,
283      x_in_l, -height, z_in_l,
284      x_in_l,  height, z_in_l,
285      x_out_l, height, z_out_l,
286      x_out_l,-height, z_out_l,
287      x_in_l, -height, z_in_l);
288%}
289
290END