1 /******************************************************************************* 2 * 3 * McStas, neutron ray-tracing package 4 * Copyright (C) 1997-2009, All rights reserved 5 * Risoe National Laboratory, Roskilde, Denmark 6 * Institut Laue Langevin, Grenoble, France 7 * 8 * Runtime: share/mcstas-r.h 9 * 10 * %Identification 11 * Written by: KN 12 * Date: Aug 29, 1997 13 * Release: McStas X.Y 14 * Version: $Revision$ 15 * 16 * Runtime system header for McStas. 17 * 18 * In order to use this library as an external library, the following variables 19 * and macros must be declared (see details in the code) 20 * 21 * struct mcinputtable_struct mcinputtable[]; 22 * int mcnumipar; 23 * char mcinstrument_name[], mcinstrument_source[]; 24 * int mctraceenabled, mcdefaultmain; 25 * extern MCNUM mccomp_storein[]; 26 * extern MCNUM mcAbsorbProp[]; 27 * extern MCNUM mcScattered; 28 * #define MCCODE_STRING "the McStas version" 29 * 30 * Usage: Automatically embbeded in the c code. 31 * 32 * $Id$ 33 * 34 *******************************************************************************/ 35 36 #ifndef MCSTAS_R_H 37 #define MCSTAS_R_H "$Revision$" 38 39 /* Following part is only embedded when not redundent with mcstas.h ========= */ 40 41 #ifndef MCCODE_H 42 43 #define AA2MS 629.622368 /* Convert k[1/AA] to v[m/s] */ 44 #define MS2AA 1.58825361e-3 /* Convert v[m/s] to k[1/AA] */ 45 #define K2V AA2MS 46 #define V2K MS2AA 47 #define Q2V AA2MS 48 #define V2Q MS2AA 49 #define SE2V 437.393377 /* Convert sqrt(E)[meV] to v[m/s] */ 50 #define VS2E 5.22703725e-6 /* Convert (v[m/s])**2 to E[meV] */ 51 52 #define SCATTER do {mcDEBUG_SCATTER(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz, \ 53 mcnlt,mcnlsx,mcnlsy,mcnlsz, mcnlp); mcScattered++;} while(0) 54 #define ABSORB do {mcDEBUG_STATE(mcnlx, mcnly, mcnlz, mcnlvx, mcnlvy, mcnlvz, \ 55 mcnlt,mcnlsx,mcnlsy,mcnlsz, mcnlp); mcDEBUG_ABSORB(); MAGNET_OFF; goto mcabsorb;} while(0) 56 57 #define STORE_NEUTRON(index, x, y, z, vx, vy, vz, t, sx, sy, sz, p) \ 58 mcstore_neutron(mccomp_storein,index, x, y, z, vx, vy, vz, t, sx, sy, sz, p); 59 #define RESTORE_NEUTRON(index, x, y, z, vx, vy, vz, t, sx, sy, sz, p) \ 60 mcrestore_neutron(mccomp_storein,index, &x, &y, &z, &vx, &vy, &vz, &t, &sx, &sy, &sz, &p); 61 62 #define MAGNET_ON \ 63 do { \ 64 mcMagnet = 1; \ 65 } while(0) 66 67 #define MAGNET_OFF \ 68 do { \ 69 mcMagnet = 0; \ 70 } while(0) 71 72 #define ALLOW_BACKPROP \ 73 do { \ 74 mcallowbackprop = 1; \ 75 } while(0) 76 77 #define DISALLOW_BACKPROP \ 78 do { \ 79 mcallowbackprop = 0; \ 80 } while(0) 81 82 #define PROP_MAGNET(dt) \ 83 do { \ 84 }while (0) 85 /* change coordinates from local system to magnet system */ 86 /* Rotation rotLM, rotTemp; \ 87 Coords posLM = coords_sub(POS_A_CURRENT_COMP, mcMagnetPos); \ 88 rot_transpose(ROT_A_CURRENT_COMP, rotTemp); \ 89 rot_mul(rotTemp, mcMagnetRot, rotLM); \ 90 mcMagnetPrecession(mcnlx, mcnly, mcnlz, mcnlt, mcnlvx, mcnlvy, mcnlvz, \ 91 &mcnlsx, &mcnlsy, &mcnlsz, dt, posLM, rotLM); \ 92 } while(0) 93 */ 94 95 #define mcPROP_DT(dt) \ 96 do { \ 97 if (mcMagnet && dt > 0) PROP_MAGNET(dt);\ 98 mcnlx += mcnlvx*(dt); \ 99 mcnly += mcnlvy*(dt); \ 100 mcnlz += mcnlvz*(dt); \ 101 mcnlt += (dt); \ 102 if (isnan(p) || isinf(p)) { mcAbsorbProp[INDEX_CURRENT_COMP]++; ABSORB; }\ 103 } while(0) 104 105 /* ADD: E. Farhi, Aug 6th, 2001 PROP_GRAV_DT propagation with acceleration */ 106 #define PROP_GRAV_DT(dt, Ax, Ay, Az) \ 107 do { \ 108 if(dt < 0 && mcallowbackprop == 0) { mcAbsorbProp[INDEX_CURRENT_COMP]++; ABSORB; }\ 109 if (mcMagnet) printf("Spin precession gravity\n"); \ 110 mcnlx += mcnlvx*(dt) + (Ax)*(dt)*(dt)/2; \ 111 mcnly += mcnlvy*(dt) + (Ay)*(dt)*(dt)/2; \ 112 mcnlz += mcnlvz*(dt) + (Az)*(dt)*(dt)/2; \ 113 mcnlvx += (Ax)*(dt); \ 114 mcnlvy += (Ay)*(dt); \ 115 mcnlvz += (Az)*(dt); \ 116 mcnlt += (dt); \ 117 DISALLOW_BACKPROP;\ 118 } while(0) 119 120 121 #define PROP_DT(dt) \ 122 do { \ 123 if(dt < 0) { RESTORE=1; goto mcabsorbComp; }; \ 124 if (mcgravitation) { Coords mcLocG; double mc_gx, mc_gy, mc_gz; \ 125 mcLocG = rot_apply(ROT_A_CURRENT_COMP, coords_set(0,-GRAVITY,0)); \ 126 coords_get(mcLocG, &mc_gx, &mc_gy, &mc_gz); \ 127 PROP_GRAV_DT(dt, mc_gx, mc_gy, mc_gz); } \ 128 else mcPROP_DT(dt); \ 129 DISALLOW_BACKPROP;\ 130 } while(0) 131 132 133 #define PROP_Z0 \ 134 do { \ 135 if (mcgravitation) { Coords mcLocG; int mc_ret; \ 136 double mc_dt, mc_gx, mc_gy, mc_gz; \ 137 mcLocG = rot_apply(ROT_A_CURRENT_COMP, coords_set(0,-GRAVITY,0)); \ 138 coords_get(mcLocG, &mc_gx, &mc_gy, &mc_gz); \ 139 mc_ret = solve_2nd_order(&mc_dt, NULL, -mc_gz/2, -mcnlvz, -mcnlz); \ 140 if (mc_ret && mc_dt>=0) {PROP_GRAV_DT(mc_dt, mc_gx, mc_gy, mc_gz); mcnlz=0;}\ 141 else { if (mcallowbackprop ==0) {mcAbsorbProp[INDEX_CURRENT_COMP]++; ABSORB; }}; }\ 142 else mcPROP_Z0; \ 143 DISALLOW_BACKPROP;\ 144 } while(0) 145 146 #define mcPROP_Z0 \ 147 do { \ 148 double mc_dt; \ 149 if(mcnlvz == 0) { mcAbsorbProp[INDEX_CURRENT_COMP]++; ABSORB; }; \ 150 mc_dt = -mcnlz/mcnlvz; \ 151 if(mc_dt < 0 && mcallowbackprop == 0) { mcAbsorbProp[INDEX_CURRENT_COMP]++; ABSORB; }; \ 152 mcPROP_DT(mc_dt); \ 153 mcnlz = 0; \ 154 DISALLOW_BACKPROP;\ 155 } while(0) 156 157 #define PROP_X0 \ 158 do { \ 159 if (mcgravitation) { Coords mcLocG; int mc_ret; \ 160 double mc_dt, mc_gx, mc_gy, mc_gz; \ 161 mcLocG = rot_apply(ROT_A_CURRENT_COMP, coords_set(0,-GRAVITY,0)); \ 162 coords_get(mcLocG, &mc_gx, &mc_gy, &mc_gz); \ 163 mc_ret = solve_2nd_order(&mc_dt, NULL, -mc_gx/2, -mcnlvx, -mcnlx); \ 164 if (mc_ret && mc_dt>=0) PROP_GRAV_DT(mc_dt, mc_gx, mc_gy, mc_gz); \ 165 else { if (mcallowbackprop ==0) {mcAbsorbProp[INDEX_CURRENT_COMP]++; ABSORB; }}; }\ 166 else mcPROP_X0; \ 167 DISALLOW_BACKPROP;\ 168 } while(0) 169 170 #define mcPROP_X0 \ 171 do { \ 172 double mc_dt; \ 173 if(mcnlvx == 0) { mcAbsorbProp[INDEX_CURRENT_COMP]++; ABSORB; }; \ 174 mc_dt = -mcnlx/mcnlvx; \ 175 if(mc_dt < 0 && mcallowbackprop == 0) { mcAbsorbProp[INDEX_CURRENT_COMP]++; ABSORB; }; \ 176 mcPROP_DT(mc_dt); \ 177 mcnlx = 0; \ 178 DISALLOW_BACKPROP;\ 179 } while(0) 180 181 #define PROP_Y0 \ 182 do { \ 183 if (mcgravitation) { Coords mcLocG; int mc_ret; \ 184 double mc_dt, mc_gx, mc_gy, mc_gz; \ 185 mcLocG = rot_apply(ROT_A_CURRENT_COMP, coords_set(0,-GRAVITY,0)); \ 186 coords_get(mcLocG, &mc_gx, &mc_gy, &mc_gz); \ 187 mc_ret = solve_2nd_order(&mc_dt, NULL, -mc_gy/2, -mcnlvy, -mcnly); \ 188 if (mc_ret && mc_dt>=0) PROP_GRAV_DT(mc_dt, mc_gx, mc_gy, mc_gz); \ 189 else { if (mcallowbackprop ==0) {mcAbsorbProp[INDEX_CURRENT_COMP]++; ABSORB; }}; }\ 190 else mcPROP_Y0; \ 191 DISALLOW_BACKPROP;\ 192 } while(0) 193 194 195 #define mcPROP_Y0 \ 196 do { \ 197 double mc_dt; \ 198 if(mcnlvy == 0) { mcAbsorbProp[INDEX_CURRENT_COMP]++; ABSORB; }; \ 199 mc_dt = -mcnly/mcnlvy; \ 200 if(mc_dt < 0 && mcallowbackprop == 0) { mcAbsorbProp[INDEX_CURRENT_COMP]++; ABSORB; }; \ 201 mcPROP_DT(mc_dt); \ 202 mcnly = 0; \ 203 DISALLOW_BACKPROP; \ 204 } while(0) 205 206 /*moved from mccode-r.h*/ 207 void mcsetstate(double x, double y, double z, double vx, double vy, double vz, 208 double t, double sx, double sy, double sz, double p); 209 210 #ifdef DEBUG 211 212 #define mcDEBUG_STATE(x,y,z,vx,vy,vz,t,sx,sy,sz,p) if(!mcdotrace); else \ 213 printf("STATE: %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g\n", \ 214 x,y,z,vx,vy,vz,t,sx,sy,sz,p); 215 #define mcDEBUG_SCATTER(x,y,z,vx,vy,vz,t,sx,sy,sz,p) if(!mcdotrace); else \ 216 printf("SCATTER: %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g\n", \ 217 x,y,z,vx,vy,vz,t,sx,sy,sz,p); 218 219 #else 220 221 #define mcDEBUG_STATE(x,y,z,vx,vy,vz,t,sx,sy,sz,p) 222 #define mcDEBUG_SCATTER(x,y,z,vx,vy,vz,t,sx,sy,sz,p) 223 224 #endif 225 226 #endif /* !MCCODE_H */ 227 228 #endif /* MCSTAS_R_H */ 229