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