1/**************************************************************************
2*
3* McStas, neutron ray-tracing package
4*         Copyright 1997-2006, All rights reserved
5*         Risoe National Laboratory, Roskilde, Denmark
6*         Institut Laue Langevin, Grenoble, France
7*
8* Component: Pol_triafield
9*
10* %I
11* Written by: Morten Sales, based on Pol_constBfield by Peter Christiansen
12* Date: 2013
13* Origin: Helmholtz-Zentrum Berlin
14*
15* Constant magnetic field in a isosceles triangular coil
16*
17* %D
18*
19* Rectangular box with constant B field along y-axis (up) in a isosceles triangle.
20* There is a guide (or precession) field as well. It is along y in the entire rectangular box.
21* A neutron hitting outside the box opening or the box sides is absorbed.
22*
23*
24*     __________________
25*    |        /\        |
26*    | Bguide/  \Bguide |      x
27*    |      /    \      |      ^
28*    |     /      \     |      |
29*    |    /   B    \    |      |-----> z
30*    |   /   and    \   |
31*    |  /   Bguide   \  |
32*    | /              \ |
33*    |/________________\|
34*
35* The angle of the inclination of the triangular field boundary is given by the arctangent to xwidth/(0.5*zdepth)
36*
37* This component does NOT take gravity into account.
38*
39* Example: Pol_triafield(xwidth=0.1, yheight=0.1, zdepth=0.2, B=1e-3, Bguide=0.0)
40*
41* %P
42* INPUT PARAMETERS:
43*
44* xwidth: [m]   Width of opening
45* yheight: [m]  Height of opening
46* zdepth: [m]   zdepth of field
47* B: [T]        Magnetic field along y-direction inside triangle
48* Bguide: [T]   Magnetic field along y-direction inside entire box
49*
50* OUTPUT PARAMETERS:
51*
52* %E
53*******************************************************************************/
54
55DEFINE COMPONENT Pol_triafield
56SETTING PARAMETERS (xwidth, yheight, zdepth, B=0, Bguide=0)
57OUTPUT PARAMETERS(omegaL, omegaLguide)
58/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
59
60SHARE
61%{
62double IntersectWall(double pos, double vel, double wallpos) {
63    /* Function to calculate where the neutron hit the wall */
64
65    if(vel==0)
66      return -1;
67
68    if(vel>0)
69      return (wallpos-pos)/vel;
70    else
71      return (-wallpos-pos)/vel;
72  }
73%}
74
75DECLARE
76%{
77  /* Larmor frequency */
78  double omegaL;
79  double omegaLguide;
80%}
81
82INITIALIZE
83%{
84  omegaL = 0;
85  omegaLguide = 0;
86
87
88  double velocity = 0, time = 0;
89
90  if ((xwidth<=0) || (yheight<=0) || (zdepth<=0)) {
91    fprintf(stderr, "Pol_filter: %s: Null or negative volume!\n"
92	    "ERROR      (xwidth, yheight, zdepth). Exiting\n",
93	    NAME_CURRENT_COMP);
94    exit(1);
95  }
96
97  omegaL	  = -1.832472e8 * (B - Bguide); // B and Bguide is in Tesla
98  omegaLguide = -1.832472e8 * Bguide;       // Bguide is in Tesla
99  %}
100
101TRACE
102%{
103  double deltaT, deltaTx, deltaTy, sx_in1, sz_in1, sx_in2, sz_in2, iz1, iz2, denom1, denom2, deltaTtria;
104
105  PROP_Z0;
106  if (!inside_rectangle(x, y, xwidth, yheight))
107    ABSORB;
108
109  // Time spent in Bguide-field
110  deltaT = zdepth/vz;
111
112  // This calculates the intersections on the xz-plane between the neutron trajectory and the triangular field boundaries
113  // The neutron trajectory is given by the points    (        x, 0,       0) and (     x+vx, 0,       vz)
114  // The first field boundary is given by the points  (-xwidth/2, 0,       0) and ( xwidth/2, 0, zdepth/2)
115  // The second field boundary is given by the points ( xwidth/2, 0,zdepth/2) and (-xwidth/2, 0,   zdepth)
116  // iz1 and iz2 are the z-values for the intersection
117    denom1 = (-vz)*((-xwidth/2)-xwidth/2)-(x-(x+vx))*(-zdepth/2);
118    iz1    = ((-x*vz)*(-zdepth/2)-(-vz)*(-(-xwidth/2)*zdepth/2))/denom1;
119
120    denom2 = (-vz)*(xwidth/2-(-xwidth/2))-(x-(x+vx))*(zdepth/2-zdepth);
121    iz2    = ((-x*vz)*(zdepth/2-zdepth)-(-vz)*(zdepth/2*(-xwidth/2)-xwidth/2*zdepth))/denom2;
122  // Time spent in triangular B-field
123    deltaTtria	= (iz2-iz1)/vz;
124
125  // check that track goes throgh without hitting the walls
126  if (!inside_rectangle(x+vx*deltaT, y+vy*deltaT, xwidth, yheight)) {
127
128    // Propagate to the wall and absorb
129    deltaTx = IntersectWall(x, vx, xwidth/2);
130    deltaTy = IntersectWall(y, vy, yheight/2);
131
132    if (deltaTx>=0 && deltaTx<deltaTy)
133      deltaT = deltaTx;
134    else
135      deltaT = deltaTy;
136
137    PROP_DT(deltaT);
138
139    ABSORB;
140  }
141
142  PROP_DT(deltaT);
143
144  // These are the incoming spin directions
145  sx_in1 = sx;
146  sz_in1 = sz;
147
148  // This calculates the spin rotation caused by the guide/precession field
149  sz_in2 = cos(omegaLguide*deltaT)*sz_in1 - sin(omegaLguide*deltaT)*sx_in1;
150  sx_in2 = sin(omegaLguide*deltaT)*sz_in1 + cos(omegaLguide*deltaT)*sx_in1;
151
152  // This calculated the spin rotation caused by the triangular field
153  sz = cos(omegaL*deltaTtria)*sz_in2 - sin(omegaL*deltaTtria)*sx_in2;
154  sx = sin(omegaL*deltaTtria)*sz_in2 + cos(omegaL*deltaTtria)*sx_in2;
155  %}
156
157MCDISPLAY
158%{
159
160
161  box(0, 0, 0, xwidth, yheight, zdepth);
162  %}
163
164END
165