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 * Library: share/pol-lib.c
9 *
10 * %Identification
11 * Written by: Erik Knudsen, Astrid Rømer & Peter Christiansen
12 * Date: Oct 08
13 * Origin: RISOE
14 * Release: McStas 1.12
15 * Version: $Revision: 4466 $
16 *
17 * This file is to be imported by polarisation components.
18 * It handles some shared functions.
19 * Embedded within instrument in runtime mode.
20 * Variable names have prefix 'mc_pol_' for 'McStas Polarisation'
21 * to avoid conflicts
22 *
23 * Usage: within SHARE
24 * %include "pol-lib"
25 *
26 ****************************************************************************/
27
28 #ifndef POL_LIB_H
29 #include "pol-lib.h"
30 #endif
31
32 #include<sys/stat.h>
33
34 %include "read_table-lib"
35 %include "interpolation-lib"
36
37 enum {MCMAGNET_STACKSIZE=12} mcmagnet_constants;
38
39 /*definition of the magnetic stack*/
40 static mcmagnet_field_info *stack[MCMAGNET_STACKSIZE];
41 /*definition of the precession function*/
42 #ifdef MC_POL_COMPAT
43 extern mcmagnet_prec_func *mcMagnetPrecession;
44 extern Coords mcMagnetPos;
45 extern Rotation mcMagnetRot;
46 extern double* mcMagnetData;
47 /* mcMagneticField(x, y, z, t, Bx, By, Bz) */
48 extern int (*mcMagneticField) (double, double, double, double,
49 double*, double*, double*, void *);
50 #else
51 #ifndef POL_LIB_C
52 static mcmagnet_prec_func *mcMagnetPrecession=SimpleNumMagnetPrecession;
53 static Coords mcMagnetPos;
54 static Rotation mcMagnetRot;
55 static double* mcMagnetData = NULL;
56 static int (*mcMagneticField) (double, double, double, double,
57 double*, double*, double*, void *);
58 /*Threshold below which two magnetic fields are considered to be
59 * in the same direction.*/
60 static double mc_pol_angular_accuracy = 1.0*DEG2RAD; /*rad.*/
61 /*The maximal timestep taken by neutrons in a const field*/
62 static double mc_pol_initial_timestep = 1e-5;
63 #define POL_LIB_C 1
64 #endif
65 #endif
66
mcmagnet_init()67 int mcmagnet_init(){
68 mcMagnetPrecession=SimpleNumMagnetPrecession;
69 return 1;
70 }
71
mc_pol_set_angular_accuracy(double domega)72 void mc_pol_set_angular_accuracy(double domega){
73 mc_pol_angular_accuracy = domega;
74 }
75
mc_pol_set_timestep(double dt)76 void mc_pol_set_timestep(double dt){
77 mc_pol_initial_timestep=dt;
78 }
79
80 #ifdef PROP_MAGNET
81 #undef PROP_MAGNET
82 #define PROP_MAGNET(dt) \
83 do { \
84 /* change coordinates from local system to magnet system */ \
85 Rotation rotLM; \
86 Coords posLM = POS_A_CURRENT_COMP; \
87 rot_transpose(ROT_A_CURRENT_COMP, rotLM); \
88 mcMagnetPrecession(mcnlx, mcnly, mcnlz, mcnlt, mcnlvx, mcnlvy, mcnlvz, \
89 &mcnlsx, &mcnlsy, &mcnlsz, dt, posLM, rotLM); \
90 } while(0)
91 #endif
92
93 /*traverse the stack and return the magnetic field*/
mcmagnet_get_field(double x,double y,double z,double t,double * bx,double * by,double * bz,void * dummy)94 int mcmagnet_get_field(double x, double y, double z, double t, double *bx,double *by, double *bz, void *dummy){
95 mcmagnet_field_info *p=stack[0];
96 Coords in,loc,b,bsum={0,0,0},zero={0,0,0};
97 Rotation r;
98
99 /*PROP_MAGNET takes care of transforming local "PROP" coordinates to lab system*/
100 in.x=x;in.y=y;in.z=z;
101
102 int i=0,stat=1;
103 p=stack[i];
104 *bx=0;*by=0;*bz=0;
105 if (!p) return 0;
106 //mcmagnet_print_stack();
107 //printf("getfield_(lab):_(xyz,t)=( %g %g %g %g )\n",x,y,z,t);
108 while(p){
109 /*transform to the coordinate system of the particular magnetic function*/
110 loc=coords_sub(rot_apply(*(p->rot),in),*(p->pos));
111 stat=(p->func) (loc.x,loc.y,loc.z,t,&(b.x),&(b.y),&(b.z),p->data);
112 /*check if the field function should be garbage collected*/
113 //printf("getfield_(loc):_(xyz,t)=( %g %g %g %g )\n",loc.x,loc.y,loc.z,t);
114 if (stat){
115 /*transform to the lab system and add up. (resusing loc variable - to now contain the field in lab coords)*/
116 rot_transpose(*(p->rot),r);
117 loc=rot_apply(r,b);
118 bsum.x+=loc.x;bsum.y+=loc.y;bsum.z+=loc.z;
119 //printf("Bs=(%g %g %g), B=(%g %g %g)\n",bsum.x,bsum.y,bsum.z,loc.x,loc.y,loc.z);
120 }
121 if (p->stop) break;
122 p=stack[++i];
123 }
124 /*we now have the magnetic field in lab coords in loc., transfer it back to caller*/
125 *bx=bsum.x;
126 *by=bsum.y;
127 *bz=bsum.z;
128 return 1;
129 }
130
131 /*void mcmagnet_init(void){
132 mcmagnet_field_info *p;
133 for (p=&(stack[0]);p<&(stack[MCMAGNET_STACKSIZE]);p++){
134 *p = malloc (sizeof(mcmagnet_field_info));
135 }
136 }
137 */
mcmagnet_push(mcmagnet_field_func * func,Rotation * magnet_rot,Coords * magnet_pos,int stopbit,void * prms)138 void *mcmagnet_push(mcmagnet_field_func *func, Rotation *magnet_rot, Coords *magnet_pos, int stopbit, void * prms){
139 mcmagnet_field_info *p;
140 int i;
141 /*move the stack one step down start from -2 since we have 0-indexing (i.e. last item is stacksize-1) */
142 for (i=MCMAGNET_STACKSIZE-2;i>=0;i--){
143 stack[i+1]=stack[i];
144 }
145 stack[0]=(mcmagnet_field_info *)malloc(sizeof(mcmagnet_field_info));
146 mcmagnet_pack(stack[0],func,magnet_rot,magnet_pos,stopbit,prms);
147 mcmagnet_set_active(stack[0]);
148 if(stack[0] && stack[0]->func){
149 MAGNET_ON;
150 }
151 return (void *) stack[0];
152 }
153
mcmagnet_pop(void)154 void *mcmagnet_pop(void) {
155 mcmagnet_field_info **p,*t;
156 /*move the stack one step up*/
157 int i;
158 t=stack[0];
159 for (i=0;i<MCMAGNET_STACKSIZE-2;i++){
160 stack[i]=stack[i+1];
161 }
162 stack[MCMAGNET_STACKSIZE-1]=NULL;
163 mcmagnet_set_active(stack[0]);
164 if(stack[0] && stack[0]->func){
165 MAGNET_ON;
166 }else{
167 MAGNET_OFF;
168 }
169 return (void*) t;
170 }
171
mcmagnet_free_stack(void)172 void mcmagnet_free_stack(void){
173 mcmagnet_field_info **p;
174 for (p=&(stack[0]);p<&(stack[MCMAGNET_STACKSIZE]);p++){
175 free(*p);
176 }
177 }
178
mcmagnet_init_par_backend(int dummy,...)179 void *mcmagnet_init_par_backend(int dummy, ...){
180 void * data;
181 unsigned char *p=NULL;
182 int q,dp=0;
183 va_list arg_list;
184
185 va_start(arg_list,dummy);
186 p=(unsigned char *)arg_list;
187 q=va_arg(arg_list,int);
188 while (q!=MCMAGNET_STOP_ARG){
189 q=va_arg(arg_list,int);
190 }
191 dp=(unsigned char *)arg_list-p;
192 data=(void *) malloc(sizeof(int)*dp);
193 memcpy(data,p,sizeof(int)*dp);
194 return data;
195 }
196
mcmagnet_print_active()197 void mcmagnet_print_active(){
198 Rotation *p;
199 printf("address of magnetic field function:%p\n",mcMagneticField);
200 p=&mcMagnetRot;
201 printf("rotation matrix of magnetic field:[%g %g %g; %g %g %g; %g %g %g]\n",(*p)[0][0],(*p)[0][1],(*p)[0][2],(*p)[1][0],(*p)[1][1],(*p)[1][2],(*p)[2][0],(*p)[2][1],(*p)[2][2]);
202 printf("origin position of magnet (x,y,z) :[%g %g %g]\n",mcMagnetPos.x,mcMagnetPos.y,mcMagnetPos.z);
203 printf("address of magnetic field parameters: %p\n",mcMagnetData);
204 }
205
mcmagnet_print_field(mcmagnet_field_info * magnet)206 void mcmagnet_print_field(mcmagnet_field_info *magnet){
207 Rotation *p;
208 if (magnet!=NULL){
209 printf("address of magnetic field function:%p\n",magnet->func);
210 p=magnet->rot;
211 printf("rotation matrix of magnetic field:[%g %g %g; %g %g %g; %g %g %g]\n",(*p)[0][0],(*p)[0][1],(*p)[0][2],(*p)[1][0],(*p)[1][1],(*p)[1][2],(*p)[2][0],(*p)[2][1],(*p)[2][2]);
212 printf("origin position of magnet (x,y,z) :[%g %g %g]\n",magnet->pos->x,magnet->pos->y,magnet->pos->z);
213 printf("address of magnetic field parameters: %p\n",magnet->data);
214 } else {
215 printf("magnet is NULL\n");
216 }
217 }
218
mcmagnet_print_stack()219 void mcmagnet_print_stack(){
220 mcmagnet_field_info *p=stack[0];
221 int i=0;
222 p=stack[i];
223 printf("magnetic stack info:\n");
224 if (!p) return;
225 while(p) {
226 printf("magnet %d:\n",i);
227 mcmagnet_print_field(p);
228 if (p->stop) break;
229 p=stack[++i];
230 }
231 }
232
233
234 /*Example magnetic field functions*/
const_magnetic_field(double x,double y,double z,double t,double * bx,double * by,double * bz,void * data)235 int const_magnetic_field(double x, double y, double z, double t,
236 double *bx, double *by, double *bz, void *data) {
237 int stat=1;
238 if (!data) return 0;
239 *bx=((double *)data)[0];
240 *by=((double *)data)[1];
241 *bz=((double *)data)[2];
242 return stat;
243 }
244
rot_magnetic_field(double x,double y,double z,double t,double * bx,double * by,double * bz,void * data)245 int rot_magnetic_field(double x, double y, double z, double t,
246 double *bx, double *by, double *bz, void *data) {
247 /* Field of magnitude By that rotates to x in magnetLength m*/
248
249 if (!data) return 0;
250 double Bmagnitude=((double *)data)[0];// = mcMagnetData[1];
251 double magnetLength=((double *)data)[1];// = mcMagnetData[5];
252 *bx = Bmagnitude * sin(PI/2*(z+magnetLength/2.0)/magnetLength);
253 *by = Bmagnitude * cos(PI/2*(z+magnetLength/2.0)/magnetLength);
254 *bz = 0;
255 //printf("mag field at (x,y,z)=( %g %g %g ) t=%g is B=( %g %g %g )\n",x,y,z,t,*bx,*by,*bz);
256 return 1;
257 }
258
majorana_magnetic_field(double x,double y,double z,double t,double * bx,double * by,double * bz,void * data)259 int majorana_magnetic_field(double x, double y, double z, double t,
260 double *bx, double *by, double *bz, void *data) {
261 /* Large linearly decreasing (from +Bx to -Bx in magnetLength) component along x axis,
262 * small constant component along y axis
263 */
264 if (!data) return 0;
265 double Blarge = ((double *)data)[0];
266 double Bsmall = ((double *)data)[1];
267 double magnetLength = ((double *)data)[2];
268 *bx = Blarge -2*Blarge*z/magnetLength;
269 *by = Bsmall;
270 *bz = 0;
271 return 1;
272 }
273
table_magnetic_field(double x,double y,double z,double t,double * bx,double * by,double * bz,void * data)274 int table_magnetic_field(double x, double y, double z, double t,
275 double *bx, double *by, double *bz,
276 void *data)
277 {
278 if (!data) return 0;
279 struct interpolator_struct *interpolator = (struct interpolator_struct*)data;
280 return(interpolator_interpolate3_3(interpolator, x,y,z, bx,by,bz) != NULL);
281 }
282
283
284 /****************************************************************************
285 * void GetMonoPolFNFM(double Rup, double Rdown, double *FN, double *FM)
286 *
287 * ACTION: Calculate FN and FM from reflectivities Rup and Rdown
288 *
289 * For a monochromator (nuclear and magnetic scattering), the
290 * polarisation is done by defining the reflectivity for spin up (Rup)
291 * and spin down (Rdown) (which can be negative, see now!) and based on
292 * this the nuclear and magnetic structure factors are calculated:
293 * FM = sign(Rup)*sqrt(|Rup|) + sign(Rdown)*sqrt(|Rdown|)
294 * FN = sign(Rup)*sqrt(|Rup|) - sign(Rdown)*sqrt(|Rdown|)
295 *****************************************************************************/
GetMonoPolFNFM(double mc_pol_Rup,double mc_pol_Rdown,double * mc_pol_FN,double * mc_pol_FM)296 void GetMonoPolFNFM(double mc_pol_Rup, double mc_pol_Rdown,
297 double *mc_pol_FN, double *mc_pol_FM) {
298 if (mc_pol_Rup>0)
299 mc_pol_Rup = sqrt(fabs(mc_pol_Rup));
300 else
301 mc_pol_Rup = -sqrt(fabs(mc_pol_Rup));
302
303 if (mc_pol_Rdown>0)
304 mc_pol_Rdown = sqrt(fabs(mc_pol_Rdown));
305 else
306 mc_pol_Rdown = -sqrt(fabs(mc_pol_Rdown));
307
308 *mc_pol_FN = 0.5*(mc_pol_Rup + mc_pol_Rdown);
309 *mc_pol_FM = 0.5*(mc_pol_Rup - mc_pol_Rdown);
310 return;
311 }
312
313 /****************************************************************************
314 * void GetMonoPolRefProb(double FN, double FM, double sy, double *prob)
315 *
316 * ACTION: Calculate reflection probability from sy, FN and FM
317 *
318 * For a monochromator with up direction along y the reflection
319 * probability is given as:
320 * prob = FN*FN + 2*FN*FM*sy_in + FM*FM
321 * (= |Rup| + |Rdown| (for sy_in=0))
322 * where FN and FM are calculated from Rup and Rdown by GetMonoPolFNFM
323 *****************************************************************************/
GetMonoPolRefProb(double mc_pol_FN,double mc_pol_FM,double mc_pol_sy,double * mc_pol_prob)324 void GetMonoPolRefProb(double mc_pol_FN, double mc_pol_FM,
325 double mc_pol_sy, double *mc_pol_prob) {
326 *mc_pol_prob = mc_pol_FN*mc_pol_FN + mc_pol_FM*mc_pol_FM
327 + 2*mc_pol_FN*mc_pol_FM*mc_pol_sy;
328 return;
329 }
330
331 /****************************************************************************
332 * void SetMonoPolRefOut(double FN, double FM, double refProb,
333 * double* sx, double* sy, double* sz) {
334 *
335 * ACTION: Set the outgoing polarisation vector of the reflected neutrons
336 * given FN, FM and the reflection probability.
337 *
338 * For a monochromator with up direction along y the outgoing polarisation
339 * is given as:
340 * sx = (FN*FN - FM*FM)*sx_in/R0;
341 * sy = ((FN*FN - FM*FM)*sy_in + 2*FN*FM + FM*FM*sy_in)/R0;
342 * sz = (FN*FN - FM*FM)*sz_in/R0;
343 * where sx_in, sy_in, and sz_in is the incoming polarisation, and
344 * FN and FM are calculated from Rup and Rdown by GetMonoPolFNFM
345 *****************************************************************************/
SetMonoPolRefOut(double mc_pol_FN,double mc_pol_FM,double mc_pol_refProb,double * mc_pol_sx,double * mc_pol_sy,double * mc_pol_sz)346 void SetMonoPolRefOut(double mc_pol_FN, double mc_pol_FM,
347 double mc_pol_refProb, double* mc_pol_sx,
348 double* mc_pol_sy, double* mc_pol_sz) {
349 *mc_pol_sx = (mc_pol_FN*mc_pol_FN - mc_pol_FM*mc_pol_FM)*(*mc_pol_sx)
350 /mc_pol_refProb;
351 *mc_pol_sy = ((mc_pol_FN*mc_pol_FN - mc_pol_FM*mc_pol_FM)*(*mc_pol_sy)
352 + 2*mc_pol_FN*mc_pol_FM + 2*mc_pol_FM*mc_pol_FM*(*mc_pol_sy))
353 /mc_pol_refProb;
354 *mc_pol_sz = (mc_pol_FN*mc_pol_FN - mc_pol_FM*mc_pol_FM)*(*mc_pol_sz)
355 /mc_pol_refProb;
356 return;
357 }
358
359 /****************************************************************************
360 * void SetMonoPolTransOut(double FN, double FM, double refProb,
361 * double* sx, double* sy, double* sz) {
362 *
363 * ACTION: Set the outgoing polarisation vector of the transmitted neutrons
364 * given FN, FM and the REFLECTION probability.
365 *
366 * We use that the polarization is conserved so:
367 * s_in = refProb*s_ref+(1-refProb)*s_trans, and then
368 * s_trans = (s_in-refProb*s_ref)/(1-refProb)
369 * where refProb is calculated using the routine GetMonoPolRefProb
370 * and s_ref is calculated by SetMonoPolRefOut
371 *****************************************************************************/
SetMonoPolTransOut(double mc_pol_FN,double mc_pol_FM,double mc_pol_refProb,double * mc_pol_sx,double * mc_pol_sy,double * mc_pol_sz)372 void SetMonoPolTransOut(double mc_pol_FN, double mc_pol_FM,
373 double mc_pol_refProb, double* mc_pol_sx,
374 double* mc_pol_sy, double* mc_pol_sz) {
375 double mc_pol_sx_ref = *mc_pol_sx, mc_pol_sy_ref = *mc_pol_sy;
376 double mc_pol_sz_ref = *mc_pol_sz;
377
378 // By passing 1 as probability we get mc_pol_refProb*s_out_ref
379 SetMonoPolRefOut(mc_pol_FN, mc_pol_FM, 1,
380 &mc_pol_sx_ref, &mc_pol_sy_ref, &mc_pol_sz_ref);
381 *mc_pol_sx = (*mc_pol_sx - mc_pol_sx_ref)/(1 - mc_pol_refProb);
382 *mc_pol_sy = (*mc_pol_sy - mc_pol_sy_ref)/(1 - mc_pol_refProb);
383 *mc_pol_sz = (*mc_pol_sz - mc_pol_sz_ref)/(1 - mc_pol_refProb);
384 return;
385 }
386
387 /****************************************************************************
388 * void SimpleNumMagnetPrecession(double x, double y, double z, double t,
389 * double vx, double vy, double vz,
390 * double* sx, double* sy, double* sz, double dt)
391 *
392 *****************************************************************************/
SimpleNumMagnetPrecession(double mc_pol_x,double mc_pol_y,double mc_pol_z,double mc_pol_time,double mc_pol_vx,double mc_pol_vy,double mc_pol_vz,double * mc_pol_sx,double * mc_pol_sy,double * mc_pol_sz,double mc_pol_deltaT,Coords mc_pol_posLM,Rotation mc_pol_rotLM)393 void SimpleNumMagnetPrecession(double mc_pol_x, double mc_pol_y,
394 double mc_pol_z, double mc_pol_time,
395 double mc_pol_vx, double mc_pol_vy,
396 double mc_pol_vz,
397 double* mc_pol_sx, double* mc_pol_sy,
398 double* mc_pol_sz, double mc_pol_deltaT,
399 Coords mc_pol_posLM, Rotation mc_pol_rotLM) {
400
401 double mc_pol_Bx, mc_pol_By, mc_pol_Bz, mc_pol_phiz;
402 double mc_pol_BxStart, mc_pol_ByStart, mc_pol_BzStart, mc_pol_Bstart;
403 double mc_pol_BxTemp, mc_pol_ByTemp, mc_pol_BzTemp, mc_pol_Btemp;
404 double mc_pol_Bstep, mc_pol_timeStep, mc_pol_sp;
405 const double mc_pol_spThreshold = cos(mc_pol_angular_accuracy);
406 const double mc_pol_startTimeStep = mc_pol_initial_timestep; // s
407 double dummy1, dummy2;
408 Rotation mc_pol_rotBack;
409
410 mcMagneticField=mcmagnet_get_field;
411
412 //printf("pos_at_caller(xyz)( %g %g %g )\n", mc_pol_x,mc_pol_y,mc_pol_z);
413 // change coordinates from current local system to lab system
414 mccoordschange(mc_pol_posLM, mc_pol_rotLM,
415 &mc_pol_x, &mc_pol_y, &mc_pol_z,
416 &mc_pol_vx, &mc_pol_vy, &mc_pol_vz, mc_pol_sx, mc_pol_sy, mc_pol_sz);
417 //printf("pos_at_labaftertranformation(xyz)( %g %g %g )\n", mc_pol_x,mc_pol_y,mc_pol_z);
418
419 // get initial B-field value
420 mcMagneticField(mc_pol_x, mc_pol_y, mc_pol_z, mc_pol_time,
421 &mc_pol_BxTemp, &mc_pol_ByTemp, &mc_pol_BzTemp,NULL);
422 do {
423 mc_pol_Bx = 0; mc_pol_By = 0; mc_pol_Bz = 0; mc_pol_phiz = 0;
424 mc_pol_BxStart = mc_pol_BxTemp; mc_pol_ByStart = mc_pol_ByTemp;
425 mc_pol_BzStart = mc_pol_BzTemp;
426 mc_pol_Bstart =
427 sqrt(mc_pol_BxStart*mc_pol_BxStart + mc_pol_ByStart*mc_pol_ByStart
428 + mc_pol_BzStart*mc_pol_BzStart);
429 mc_pol_timeStep = mc_pol_startTimeStep;
430
431 if(mc_pol_deltaT<mc_pol_timeStep)
432 mc_pol_timeStep = mc_pol_deltaT;
433
434 do {
435
436 mcMagneticField(mc_pol_x+mc_pol_vx*mc_pol_timeStep,
437 mc_pol_y+mc_pol_vy*mc_pol_timeStep,
438 mc_pol_z+mc_pol_vz*mc_pol_timeStep,
439 mc_pol_time+mc_pol_timeStep,
440 &mc_pol_BxTemp, &mc_pol_ByTemp, &mc_pol_BzTemp,NULL);
441 // not so elegant, but this is how we make sure that the steps decrease
442 // when the WHILE condition is not met
443 mc_pol_timeStep *= 0.5;
444
445 mc_pol_Btemp =
446 sqrt(mc_pol_BxTemp*mc_pol_BxTemp + mc_pol_ByTemp*mc_pol_ByTemp
447 + mc_pol_BzTemp*mc_pol_BzTemp);
448
449 mc_pol_sp =
450 scalar_prod(mc_pol_BxStart, mc_pol_ByStart, mc_pol_BzStart,
451 mc_pol_BxTemp, mc_pol_ByTemp, mc_pol_BzTemp);
452 mc_pol_sp /= mc_pol_Bstart*mc_pol_Btemp;
453
454 } while (mc_pol_sp<mc_pol_spThreshold && mc_pol_timeStep>FLT_EPSILON);
455
456 mc_pol_timeStep*=2;
457
458 // update coordinate values
459 mc_pol_x += mc_pol_vx*mc_pol_timeStep;
460 mc_pol_y += mc_pol_vy*mc_pol_timeStep;
461 mc_pol_z += mc_pol_vz*mc_pol_timeStep;
462 mc_pol_time += mc_pol_timeStep;
463 mc_pol_deltaT -= mc_pol_timeStep;
464
465 mc_pol_Bx = 0.5 * (mc_pol_BxStart + mc_pol_BxTemp);
466 mc_pol_By = 0.5 * (mc_pol_ByStart + mc_pol_ByTemp);
467 mc_pol_Bz = 0.5 * (mc_pol_BzStart + mc_pol_BzTemp);
468 mc_pol_phiz = fmod(sqrt(mc_pol_Bx*mc_pol_Bx+
469 mc_pol_By*mc_pol_By+
470 mc_pol_Bz*mc_pol_Bz)
471 *mc_pol_timeStep*mc_pol_omegaL, 2*PI);
472
473 // Do the neutron spin precession
474
475 if(!(mc_pol_Bx==0 && mc_pol_By==0 && mc_pol_Bz==0)) {
476
477 double mc_pol_sx_in = *mc_pol_sx;
478 double mc_pol_sy_in = *mc_pol_sy;
479 double mc_pol_sz_in = *mc_pol_sz;
480
481 rotate(*mc_pol_sx, *mc_pol_sy, *mc_pol_sz,
482 mc_pol_sx_in, mc_pol_sy_in, mc_pol_sz_in,
483 mc_pol_phiz, mc_pol_Bx, mc_pol_By, mc_pol_Bz);
484 }
485 } while (mc_pol_deltaT>0);
486
487 // change back spin coordinates from lab system to local system
488 rot_transpose(mc_pol_rotLM, mc_pol_rotBack);
489 mccoordschange_polarisation(mc_pol_rotBack, mc_pol_sx, mc_pol_sy, mc_pol_sz);
490
491 }
492
493 /****************************************************************************
494 * double GetConstantField(double length, double lambda, double angle)
495 *
496 * Return the magnetic field in Tesla required to flip a neutron with
497 * wavelength lambda(1/velocity), angle degrees, over the specified
498 * length(=time*velocity).
499 *
500 *****************************************************************************/
GetConstantField(double mc_pol_length,double mc_pol_lambda,double mc_pol_angle)501 double GetConstantField(double mc_pol_length, double mc_pol_lambda,
502 double mc_pol_angle)
503 {
504 const double mc_pol_velocity = K2V*2*PI/mc_pol_lambda;
505 const double mc_pol_time = mc_pol_length/mc_pol_velocity;
506
507 // B*omegaL*time = angle
508 return mc_pol_angle*DEG2RAD/mc_pol_omegaL/mc_pol_time; // T
509 }
510
511 /* end of regular pol-lib.c */
512