1/****************************************************************************
2*
3* McStas, neutron ray-tracing package
4*         Copyright 1997-2003, All rights reserved
5*         Risoe National Laboratory, Roskilde, Denmark
6*         Institut Laue Langevin, Grenoble, France
7*
8* Component: Pol_guide_vmirror
9*
10* %I
11* Written by: Peter Christiansen and Erik B Knudsen
12* Date: July 2006
13* Origin: RISOE
14*
15* Polarising guide with a supermirror along its diagonal.
16*
17* %D
18* Models a rectangular guide with entrance centered on the Z axis and
19* with one supermirros sitting on the diagonal inside.
20* The entrance lies in the X-Y plane.  Draws a true depiction
21* of the guide with mirror and trajectories.
22* The polarisation is handled similar to in Monochromator_pol.
23* The reflec functions are handled similar to Pol_mirror.
24* The up direction is hardcoded to be along the y-axis (0, 1, 0)
25*
26* Note that this component can also be used as a frame overlap-mirror
27* if the up and down reflectivities are set equal. In this case the wall
28* refletivity (rPar) should probably be set to absorb.
29*
30* The parameters can either be
31* double pointer initializations (e.g. {R0, Qc, alpha, m, W})
32* or table names (e.g."supermirror_m2.rfl" AND useTables=1).
33* NB! This might cause warnings by the compiler that can be ignored.
34*
35* GRAVITY: YES
36*
37* %BUGS
38* No absorption by mirror.
39*
40* %P
41* INPUT PARAMETERS:
42*
43* xwidth: [m]      Width at the guide entry
44* yheight: [m]     Height at the guide entry
45* length: [m]      length of guide
46* rFunc: [1]       Guide Reflection function
47* rPar: [1]        Guide Parameters for rFunc
48* rUpFunc: [1]     Mirror Reflection function for spin up
49* rUpPar: [1]      Mirror Parameters for rUpFunc
50* rDownFunc: [1]   Mirror Reflection function for spin down
51* rDownPar: [1]    Mirror Parameters for rDownFunc
52* useTables: [1]   Parameters are 0: Values, 1: Table names
53* debug: [1]       if debug > 0 print out some internal runtime parameters
54*
55* OUTPUT PARAMETERS:
56*
57* localG: [m/s/s]  Gravity vector in guide reference system
58* normalTop: [1]   One of several normal vectors used for defining the geometry
59* pointTop: [1]    One of several points used for defining the geometry
60* rParPtr:      One of several pointers to reflection parameters used with the ref. functions. []
61* SCATTERED: []    is 1 for reflected, and 2 for transmitted neutrons
62*
63*
64* %L
65*
66* %E
67*******************************************************************************/
68
69DEFINE COMPONENT Pol_guide_mirror
70DEFINITION PARAMETERS (xwidth, yheight, length,
71rFunc=StdReflecFunc, rUpFunc=StdReflecFunc,
72rDownFunc=StdReflecFunc,
73rPar    ={1.0, 0.0219, 4.07, 3.2, 0.003},
74rUpPar  ={1.0, 0.0219, 4.07, 3.2, 0.003},
75rDownPar={0.1, 0.0219, 4.07, 3.2, 0.003},
76useTables=0)
77SETTING PARAMETERS (int debug=0)
78OUTPUT PARAMETERS (localG,
79normalTop, normalBot, normalLeft, normalRight, normalInOut,
80pointTop, pointBot, pointLeft, pointRight,
81pointIn, pointOut,
82rParPtr, rUpParPtr, rDownParPtr)
83/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
84
85SHARE
86%{
87%include "pol-lib"
88%include "ref-lib"
89%}
90
91DECLARE
92%{
93Coords localG;
94Coords normalTop, normalBot, normalLeft, normalRight, normalInOut;
95Coords pointTop, pointBot, pointLeft, pointRight, pointIn, pointOut;
96
97#if (useTables)
98  t_Table *rParPtr     = 0;
99  t_Table *rUpParPtr   = 0;
100  t_Table *rDownParPtr = 0;
101#else
102  double rParPtr[]     = (double [])rPar;
103  double rUpParPtr[]   = (double [])rUpPar;
104  double rDownParPtr[] = (double [])rDownPar;
105#endif
106  %}
107
108INITIALIZE
109%{
110#if (useTables)
111  rParPtr     = (t_Table*) malloc(sizeof(t_Table));
112  rUpParPtr   = (t_Table*) malloc(sizeof(t_Table));
113  rDownParPtr = (t_Table*) malloc(sizeof(t_Table));
114  if (Table_Read(rParPtr, rPar, 1) <= 0) {
115    fprintf(stderr,"Pol_guide_vmirror: %s: can not read file %s\n",
116	    NAME_CURRENT_COMP, rPar);
117    exit(1);
118  }
119  if (Table_Read(rUpParPtr, rUpPar, 1) <= 0) {
120    fprintf(stderr,"Pol_guide_vmirror: %s: can not read file %s\n",
121	    NAME_CURRENT_COMP, rUpPar);
122    exit(1);
123  }
124  if (Table_Read(rDownParPtr, rDownPar, 1) <= 0) {
125    fprintf(stderr,"Pol_guide_vmirror: %s: can not read file %s\n",
126	    NAME_CURRENT_COMP, rDownPar);
127    exit(1);
128  }
129#endif
130
131  if ((xwidth<=0) || (yheight<= 0) || (length<=0)) {
132    fprintf(stderr, "Pol_guide_vmirror: %s: NULL or negative length scale!\n"
133	                  "ERROR      (xwidth,yheight,length). Exiting\n",
134	    NAME_CURRENT_COMP);
135    exit(1);
136  }
137
138  if (mcgravitation) {
139
140    localG = rot_apply(ROT_A_CURRENT_COMP, coords_set(0,-GRAVITY,0));
141    fprintf(stdout,"Pol_guide_vmirror: %s: Gravity is on!\n",
142	    NAME_CURRENT_COMP);
143  } else
144    localG = coords_set(0, 0, 0);
145
146  // To be able to handle the situation properly where a component of
147  // the gravity is along the z-axis we also define entrance (in) and
148  // exit (out) planes
149  // The entrance and exit plane are defined by the normal vector
150  // (0, 0, 1)
151  // and the two points pointIn=(0, 0, 0) and pointOut=(0, 0, length)
152
153  normalInOut = coords_set(0, 0, 1);
154  pointIn   = coords_set(0, 0, 0);
155  pointOut  = coords_set(0, 0, length);
156
157  // Top plane (+y dir) can be spanned by (1, 0, 0) & (0, 0, 1)
158  // and the point (0, yheight/2, 0)
159  // A normal vector is (0, 1, 0)
160  normalTop  = coords_set(0, 1, 0);
161  pointTop = coords_set(0, yheight/2, 0);
162
163  // Bottom plane (-y dir) can be spanned by (1, 0, 0) & (0, 0, 1)
164  // and the point (0, -yheight/2, 0)
165  // A normal vector is (0, 1, 0)
166  normalBot  = coords_set(0, 1, 0);
167  pointBot = coords_set(0, -yheight/2, 0);
168
169  // Left plane (+x dir) can be spanned by (0, 1, 0) & (0, 0, 1)
170  // and the point (xwidth/2, 0, 0)
171  // A normal vector is (1, 0, 0)
172  normalLeft  = coords_set(1, 0, 0);
173  pointLeft = coords_set(xwidth/2, 0, 0);
174
175  // Right plane (-x dir) can be spanned by (0, 1, 0) & (0, 0, 1)
176  // and the point (-xwidth/2, 0, 0)
177  // A normal vector is (1, 0, 0)
178  normalRight  = coords_set(1, 0, 0);
179  pointRight = coords_set(-xwidth/2, 0, 0);
180
181  %}
182
183TRACE
184%{
185  /* time threshold */
186  const double tThreshold = 1e-10/sqrt(vx*vx + vy*vy + vz*vz);
187  const double xwhalf = xwidth/2;
188  const double norm = 1.0/sqrt(xwidth*xwidth + length*length);
189  double R;
190
191  Coords normalMirror, pointMirror;
192  Coords* normalPointer = 0;
193
194  // Pol variables
195  double FN, FM, Rup, Rdown, refWeight;
196
197  /* Propagate neutron to guide entrance. */
198  PROP_Z0;
199
200  if (!inside_rectangle(x, y, xwidth, yheight))
201    ABSORB;
202
203  normalMirror  = coords_set(norm*length, 0, -norm*xwidth);
204  pointMirror = coords_set(-xwhalf, 0, 0);
205
206  for(;;) {
207    double tLeft, tRight, tTop, tBot, tIn, tOut, tMirror;
208    double tUp, tSide, time, endtime;
209    double Q; //, dummy1, dummy2, dummy3;
210    Coords vVec, xVec;
211    int mirrorReflect;
212
213    mirrorReflect = 0;
214    xVec = coords_set(x, y, z);
215    vVec = coords_set(vx, vy, vz);
216
217    solve_2nd_order(&tTop, NULL, 0.5*coords_sp(normalTop,localG),
218		    coords_sp(normalTop, vVec),
219		    coords_sp(normalTop, coords_sub(xVec, pointTop)));
220
221    solve_2nd_order(&tBot, NULL, 0.5*coords_sp(normalBot,localG),
222		    coords_sp(normalBot, vVec),
223		    coords_sp(normalBot, coords_sub(xVec, pointBot)));
224
225    solve_2nd_order(&tRight, NULL, 0.5*coords_sp(normalRight,localG),
226		    coords_sp(normalRight, vVec),
227		    coords_sp(normalRight, coords_sub(xVec, pointRight)));
228
229    solve_2nd_order(&tLeft, NULL, 0.5*coords_sp(normalLeft,localG),
230		    coords_sp(normalLeft, vVec),
231		    coords_sp(normalLeft, coords_sub(xVec, pointLeft)));
232
233    solve_2nd_order(&tIn, NULL, 0.5*coords_sp(normalInOut,localG),
234		    coords_sp(normalInOut, vVec),
235		    coords_sp(normalInOut, coords_sub(xVec, pointIn)));
236
237    solve_2nd_order(&tOut, NULL, 0.5*coords_sp(normalInOut,localG),
238		    coords_sp(normalInOut, vVec),
239		    coords_sp(normalInOut, coords_sub(xVec, pointOut)));
240
241    solve_2nd_order(&tMirror, NULL, 0.5*coords_sp(normalMirror,localG),
242		    coords_sp(normalMirror, vVec),
243		    coords_sp(normalMirror, coords_sub(xVec, pointMirror)));
244
245    /* Choose appropriate reflection time */
246    if (tTop>tThreshold && (tTop<tBot || tBot<=tThreshold))
247      tUp=tTop;
248    else
249      tUp=tBot;
250
251    if (tLeft>tThreshold && (tLeft<tRight || tRight<=tThreshold))
252      tSide=tLeft;
253    else
254      tSide=tRight;
255
256    if (tUp>tThreshold && (tUp<tSide || tSide<=tThreshold))
257      time=tUp;
258    else
259      time=tSide;
260
261    if (tMirror>tThreshold && tMirror<time) {
262
263      time=tMirror;
264      mirrorReflect = 1; // flag to show which reflection function to use
265    }
266
267    if (time<=tThreshold)
268      fprintf(stdout, "Pol_guide_vmirror: %s: tTop: %f, tBot:%f, tRight: %f, tLeft: %f\n"
269	      "tUp: %f, tSide: %f, time: %f\n", NAME_CURRENT_COMP,
270	      tTop, tBot, tRight, tLeft, tUp, tSide, time);
271
272    /* Has neutron left the guide? */
273    if (tOut>tThreshold && (tOut<tIn || tIn<=tThreshold))
274      endtime=tOut;
275    else
276      endtime=tIn;
277
278    if (time > endtime)
279      break;
280
281    if(time <= tThreshold) {
282
283      printf("Time below threshold!\n");
284      fprintf(stdout, "Pol_guide_vmirror: %s: tTop: %f, tBot:%f, tRight: %f, tLeft: %f\n"
285	      "tUp: %f, tSide: %f, time: %f\n", NAME_CURRENT_COMP,
286	      tTop, tBot, tRight, tLeft, tUp, tSide, time);
287      break;
288    }
289
290    if(debug>0 && time==tLeft) {
291
292      fprintf(stdout, "\nPol_guide_vmirror: %s: Left side hit: x, v, normal, point, gravity\n", NAME_CURRENT_COMP);
293      coords_print(xVec);
294      coords_print(vVec);
295      coords_print(normalLeft);
296      coords_print(pointLeft);
297      coords_print(localG);
298
299      fprintf(stdout, "\nA: %f, B: %f, C: %f, tLeft: %f\n",
300	      0.5*coords_sp(normalLeft,localG),coords_sp(normalLeft, vVec),
301	      coords_sp(normalLeft, coords_sub(xVec, pointLeft)), tLeft);
302    }
303
304    if(debug>0)
305      fprintf(stdout, "Pol_guide_vmirror: %s: tTop: %f, tBot:%f, tRight: %f, tLeft: %f\n"
306	      "tUp: %f, tSide: %f, time: %f\n", NAME_CURRENT_COMP,
307	      tTop, tBot, tRight, tLeft, tUp, tSide, time);
308
309    if(debug>0)
310      fprintf(stdout, "Pol_guide_vmirror: %s: Start v: (%f, %f, %f)\n", NAME_CURRENT_COMP, vx, vy, vz);
311
312    PROP_DT(time);
313    if (mcgravitation)
314      vVec = coords_set(vx, vy, vz);
315    SCATTER;
316
317    if(time==tTop)
318      normalPointer = &normalTop;
319    else if(time==tBot)
320      normalPointer = &normalBot;
321    else if(time==tRight)
322      normalPointer = &normalRight;
323    else if(time==tLeft)
324      normalPointer = &normalLeft;
325    else if(time==tMirror)
326      normalPointer = &normalMirror;
327    else
328      fprintf(stderr, "Pol_guide_vmirror: %s: This should never happen!!!!\n", NAME_CURRENT_COMP);
329
330    Q = 2*coords_sp(vVec, *normalPointer)*V2K;
331
332    if(!mirrorReflect) {
333      // we have hit one of the sides. Always reflect.
334      vVec = coords_add(vVec, coords_scale(*normalPointer, -Q*K2V));
335      rFunc(fabs(Q), rParPtr, &refWeight);
336      p *= refWeight;
337
338    } else {
339      // we have hit one of the mirrors
340      rUpFunc(fabs(Q), rUpParPtr, &Rup);
341      rDownFunc(fabs(Q), rDownParPtr, &Rdown);
342      if (Rup <  0)   ABSORB;
343      if (Rup >  1)   Rup =1 ;
344      if (Rdown <  0) ABSORB;
345      if (Rdown >  1) Rdown =1 ;
346
347      GetMonoPolFNFM(Rup, Rdown, &FN, &FM);
348      GetMonoPolRefProb(FN, FM, sy, &refWeight);
349      // check that refWeight is meaningfull
350      if (refWeight <  0) ABSORB;
351      if (refWeight >  1) refWeight =1 ;
352
353      if (rand01()<refWeight) {
354	      //reflect: SCATTERED==1 for reflection
355
356	      vVec = coords_add(vVec, coords_scale(*normalPointer, -Q*K2V));
357	      SetMonoPolRefOut(FN, FM, refWeight, &sx, &sy, &sz);
358      } else {
359	      // transmit: SCATTERED==2 for transmission
360	      SCATTER;
361	      SetMonoPolTransOut(FN, FM, refWeight, &sx, &sy, &sz);
362      }
363
364      if (sx*sx+sy*sy+sz*sz>1.000001) { // check that polarisation is meaningfull
365        fprintf(stderr, "Pol_guide_vmirror: %s: polarisation |s|=%g > 1 s=[%g,%g,%g]\n",
366          NAME_CURRENT_COMP, sx*sx+sy*sy+sz*sz, sx, sy, sz);
367      }
368    }
369
370    if(p==0) {
371      ABSORB;
372      break;
373    }
374
375    // set new velocity vector
376    coords_get(vVec, &vx, &vy, &vz);
377
378    if(debug>0)
379      fprintf(stdout, "Pol_guide_vmirror: %s: End v: (%f, %f, %f)\n", NAME_CURRENT_COMP, vx, vy, vz);
380  }
381
382  %}
383
384MCDISPLAY
385%{
386  int i, j;
387
388
389
390  // draw box
391  box(0, 0, length/2.0, xwidth, yheight, length);
392
393  for(j = -1; j<=1; j+=2)
394      line(-xwidth/2.0, j*yheight/2, 0, xwidth/2.0, j*yheight/2, length);
395
396  %}
397
398END
399