| 1 | /* |
|---|
| 2 | * MonteCarlo_Main.c |
|---|
| 3 | * SANSMonteCarlo |
|---|
| 4 | * |
|---|
| 5 | * Created by Steve Kline on 7/1/10. |
|---|
| 6 | * Copyright 2010 NIST. All rights reserved. |
|---|
| 7 | * |
|---|
| 8 | */ |
|---|
| 9 | |
|---|
| 10 | // Contains all of the XOP entry points |
|---|
| 11 | // contains all of the RNG definitions |
|---|
| 12 | // contains all of the utility functions required by Monte Carlo instances |
|---|
| 13 | |
|---|
| 14 | #include "XOPStandardHeaders.h" // Include ANSI headers, Mac headers, IgorXOP.h, XOP.h and XOPSupport.h |
|---|
| 15 | #include "MonteCarlo.h" |
|---|
| 16 | #include "DebyeSpheres.h" |
|---|
| 17 | #include "Metropolis.h" |
|---|
| 18 | |
|---|
| 19 | |
|---|
| 20 | // this function is here simply because MS visual studio does not contain lround() or round() in math.h |
|---|
| 21 | // rounds away from zero |
|---|
| 22 | // -- only used in FindPixel |
|---|
| 23 | long MC_round(double x) |
|---|
| 24 | { |
|---|
| 25 | if(x == 0) { |
|---|
| 26 | return (long)(0); |
|---|
| 27 | } |
|---|
| 28 | if(x > 0) { |
|---|
| 29 | return (long)(x + 0.5f); |
|---|
| 30 | } else { // x < 0 |
|---|
| 31 | return (long)(x - 0.5f); |
|---|
| 32 | } |
|---|
| 33 | } |
|---|
| 34 | |
|---|
| 35 | |
|---|
| 36 | // corrected Feb 2013 to return values in [0,127] rather than detector coordinates (which are passed in) |
|---|
| 37 | // -- also corrected the pixel calculation |
|---|
| 38 | int |
|---|
| 39 | FindPixel(double testQ, double testPhi, double lam, double yg_d, double sdd, |
|---|
| 40 | double pixSize, double xCtr, double yCtr, long *xPixel, long *yPixel) { |
|---|
| 41 | |
|---|
| 42 | // double theta,dy,dx,qx,qy,pi; |
|---|
| 43 | double pi,two_theta,dist,dx,dy; |
|---|
| 44 | |
|---|
| 45 | pi = 4.0*atan(1.0); |
|---|
| 46 | |
|---|
| 47 | /* OLD way |
|---|
| 48 | //decompose to qx,qy |
|---|
| 49 | qx = testQ*cos(testPhi); |
|---|
| 50 | qy = testQ*sin(testPhi); |
|---|
| 51 | |
|---|
| 52 | //correct qy for gravity |
|---|
| 53 | qy += 4.0*pi/lam*(yg_d/sdd/2.0); |
|---|
| 54 | |
|---|
| 55 | |
|---|
| 56 | //convert qx,qy to pixel locations relative to # of pixels x, y from center |
|---|
| 57 | theta = 2.0*asin(qy*lam/4.0/pi); |
|---|
| 58 | dy = sdd*tan(theta); |
|---|
| 59 | // *yPixel = lround(yCtr + dy/pixSize); //corrected 7/2010 to round away from zero, to avoid 2x counts in row 0 and column 0 |
|---|
| 60 | *yPixel = MC_round(yCtr + dy/pixSize); //corrected 7/2010 to round away from zero, to avoid 2x counts in row 0 and column 0 |
|---|
| 61 | |
|---|
| 62 | theta = 2.0*asin(qx*lam/4.0/pi); |
|---|
| 63 | dx = sdd*tan(theta); |
|---|
| 64 | // *xPixel = lround(xCtr + dx/pixSize); |
|---|
| 65 | *xPixel = MC_round(xCtr + dx/pixSize); |
|---|
| 66 | |
|---|
| 67 | */ |
|---|
| 68 | |
|---|
| 69 | // corrected way Feb 2013 |
|---|
| 70 | two_theta = 2.0*asin(testQ*lam/4.0/pi); |
|---|
| 71 | dist = sdd*tan(two_theta); //hypot in det plane |
|---|
| 72 | |
|---|
| 73 | dx = dist*cos(testPhi); |
|---|
| 74 | dy = dist*sin(testPhi); |
|---|
| 75 | *xPixel = MC_round(dx/pixSize + xCtr); |
|---|
| 76 | *yPixel = MC_round(dy/pixSize + yCtr + yg_d/pixSize); //shift down due to gravity |
|---|
| 77 | |
|---|
| 78 | // to array coordinates rather than detector coordinates |
|---|
| 79 | *xPixel -= 1; |
|---|
| 80 | *yPixel -= 1; |
|---|
| 81 | |
|---|
| 82 | |
|---|
| 83 | //if on detector, return xPix and yPix values, otherwise -1 |
|---|
| 84 | if(*yPixel > 127 || *yPixel < 0) { |
|---|
| 85 | *yPixel = -1; |
|---|
| 86 | } |
|---|
| 87 | if(*xPixel > 127 || *xPixel < 0) { |
|---|
| 88 | *xPixel = -1; |
|---|
| 89 | } |
|---|
| 90 | |
|---|
| 91 | return(0); |
|---|
| 92 | } |
|---|
| 93 | |
|---|
| 94 | |
|---|
| 95 | //calculates new direction (xyz) from an old direction |
|---|
| 96 | //theta and phi don't change |
|---|
| 97 | int |
|---|
| 98 | NewDirection(double *vx, double *vy, double *vz, double theta, double phi) { |
|---|
| 99 | |
|---|
| 100 | int err=0; |
|---|
| 101 | double vx0,vy0,vz0; |
|---|
| 102 | double nx,ny,mag_xy,tx,ty,tz; |
|---|
| 103 | |
|---|
| 104 | //store old direction vector |
|---|
| 105 | vx0 = *vx; |
|---|
| 106 | vy0 = *vy; |
|---|
| 107 | vz0 = *vz; |
|---|
| 108 | |
|---|
| 109 | mag_xy = sqrt(vx0*vx0 + vy0*vy0); |
|---|
| 110 | if(mag_xy < 1e-12) { |
|---|
| 111 | //old vector lies along beam direction |
|---|
| 112 | nx = 0; |
|---|
| 113 | ny = 1; |
|---|
| 114 | tx = 1; |
|---|
| 115 | ty = 0; |
|---|
| 116 | tz = 0; |
|---|
| 117 | } else { |
|---|
| 118 | nx = -vy0 / mag_xy; |
|---|
| 119 | ny = vx0 / mag_xy; |
|---|
| 120 | tx = -vz0*vx0 / mag_xy; |
|---|
| 121 | ty = -vz0*vy0 / mag_xy; |
|---|
| 122 | tz = mag_xy ; |
|---|
| 123 | } |
|---|
| 124 | |
|---|
| 125 | //new scattered direction vector |
|---|
| 126 | *vx = cos(phi)*sin(theta)*tx + sin(phi)*sin(theta)*nx + cos(theta)*vx0; |
|---|
| 127 | *vy = cos(phi)*sin(theta)*ty + sin(phi)*sin(theta)*ny + cos(theta)*vy0; |
|---|
| 128 | *vz = cos(phi)*sin(theta)*tz + cos(theta)*vz0; |
|---|
| 129 | |
|---|
| 130 | return(err); |
|---|
| 131 | } |
|---|
| 132 | |
|---|
| 133 | double |
|---|
| 134 | path_len(double aval, double sig_tot) { |
|---|
| 135 | |
|---|
| 136 | double retval; |
|---|
| 137 | |
|---|
| 138 | retval = -1.0*log(1.0-aval)/sig_tot; |
|---|
| 139 | |
|---|
| 140 | return(retval); |
|---|
| 141 | } |
|---|
| 142 | |
|---|
| 143 | |
|---|
| 144 | #define IA 16807 |
|---|
| 145 | #define IM 2147483647 |
|---|
| 146 | #define AM (1.0/IM) |
|---|
| 147 | #define IQ 127773 |
|---|
| 148 | #define IR 2836 |
|---|
| 149 | #define NTAB 32 |
|---|
| 150 | #define NDIV (1+(IM-1)/NTAB) |
|---|
| 151 | #define EPS 1.2e-7 |
|---|
| 152 | #define RNMX (1.0-EPS) |
|---|
| 153 | |
|---|
| 154 | float ran1(long *idum) |
|---|
| 155 | { |
|---|
| 156 | int j; |
|---|
| 157 | long k; |
|---|
| 158 | static long iy=0; |
|---|
| 159 | static long iv[NTAB]; |
|---|
| 160 | float temp; |
|---|
| 161 | |
|---|
| 162 | if (*idum <= 0 || !iy) { |
|---|
| 163 | if (-(*idum) < 1) *idum=1; |
|---|
| 164 | else *idum = -(*idum); |
|---|
| 165 | for (j=NTAB+7;j>=0;j--) { |
|---|
| 166 | k=(*idum)/IQ; |
|---|
| 167 | *idum=IA*(*idum-k*IQ)-IR*k; |
|---|
| 168 | if (*idum < 0) *idum += IM; |
|---|
| 169 | if (j < NTAB) iv[j] = *idum; |
|---|
| 170 | } |
|---|
| 171 | iy=iv[0]; |
|---|
| 172 | } |
|---|
| 173 | k=(*idum)/IQ; |
|---|
| 174 | *idum=IA*(*idum-k*IQ)-IR*k; |
|---|
| 175 | if (*idum < 0) *idum += IM; |
|---|
| 176 | j=iy/NDIV; |
|---|
| 177 | iy=iv[j]; |
|---|
| 178 | iv[j] = *idum; |
|---|
| 179 | if ((temp=AM*iy) > RNMX) return RNMX; |
|---|
| 180 | else return temp; |
|---|
| 181 | } |
|---|
| 182 | #undef IA |
|---|
| 183 | #undef IM |
|---|
| 184 | #undef AM |
|---|
| 185 | #undef IQ |
|---|
| 186 | #undef IR |
|---|
| 187 | #undef NTAB |
|---|
| 188 | #undef NDIV |
|---|
| 189 | #undef EPS |
|---|
| 190 | #undef RNMX |
|---|
| 191 | |
|---|
| 192 | |
|---|
| 193 | //////////a complete copy of ran1(), simply renamed ran1a() |
|---|
| 194 | #define IA2 16807 |
|---|
| 195 | #define IM2 2147483647 |
|---|
| 196 | #define AM2 (1.0/IM2) |
|---|
| 197 | #define IQ2 127773 |
|---|
| 198 | #define IR2 2836 |
|---|
| 199 | #define NTAB2 32 |
|---|
| 200 | #define NDIV2 (1+(IM2-1)/NTAB2) |
|---|
| 201 | #define EPS2 1.2e-7 |
|---|
| 202 | #define RNMX2 (1.0-EPS2) |
|---|
| 203 | |
|---|
| 204 | float ran1a(long *idum) |
|---|
| 205 | { |
|---|
| 206 | int j; |
|---|
| 207 | long k; |
|---|
| 208 | static long iy=0; |
|---|
| 209 | static long iv[NTAB2]; |
|---|
| 210 | float temp; |
|---|
| 211 | |
|---|
| 212 | if (*idum <= 0 || !iy) { |
|---|
| 213 | if (-(*idum) < 1) *idum=1; |
|---|
| 214 | else *idum = -(*idum); |
|---|
| 215 | for (j=NTAB2+7;j>=0;j--) { |
|---|
| 216 | k=(*idum)/IQ2; |
|---|
| 217 | *idum=IA2*(*idum-k*IQ2)-IR2*k; |
|---|
| 218 | if (*idum < 0) *idum += IM2; |
|---|
| 219 | if (j < NTAB2) iv[j] = *idum; |
|---|
| 220 | } |
|---|
| 221 | iy=iv[0]; |
|---|
| 222 | } |
|---|
| 223 | k=(*idum)/IQ2; |
|---|
| 224 | *idum=IA2*(*idum-k*IQ2)-IR2*k; |
|---|
| 225 | if (*idum < 0) *idum += IM2; |
|---|
| 226 | j=iy/NDIV2; |
|---|
| 227 | iy=iv[j]; |
|---|
| 228 | iv[j] = *idum; |
|---|
| 229 | if ((temp=AM2*iy) > RNMX2) return RNMX2; |
|---|
| 230 | else return temp; |
|---|
| 231 | } |
|---|
| 232 | #undef IA2 |
|---|
| 233 | #undef IM2 |
|---|
| 234 | #undef AM2 |
|---|
| 235 | #undef IQ2 |
|---|
| 236 | #undef IR2 |
|---|
| 237 | #undef NTAB2 |
|---|
| 238 | #undef NDIV2 |
|---|
| 239 | #undef EPS2 |
|---|
| 240 | #undef RNMX2 |
|---|
| 241 | /////////////// |
|---|
| 242 | |
|---|
| 243 | |
|---|
| 244 | //////////////////////// |
|---|
| 245 | #define MBIG 1000000000 |
|---|
| 246 | #define MSEED 161803398 |
|---|
| 247 | #define MZ 0 |
|---|
| 248 | #define FAC (1.0/MBIG) |
|---|
| 249 | |
|---|
| 250 | float ran3(long *idum) |
|---|
| 251 | { |
|---|
| 252 | static int inext,inextp; |
|---|
| 253 | static long ma[56]; |
|---|
| 254 | static int iff=0; |
|---|
| 255 | long mj,mk; |
|---|
| 256 | int i,ii,k; |
|---|
| 257 | |
|---|
| 258 | if (*idum < 0 || iff == 0) { |
|---|
| 259 | iff=1; |
|---|
| 260 | mj=MSEED-(*idum < 0 ? -*idum : *idum); |
|---|
| 261 | mj %= MBIG; |
|---|
| 262 | ma[55]=mj; |
|---|
| 263 | mk=1; |
|---|
| 264 | for (i=1;i<=54;i++) { |
|---|
| 265 | ii=(21*i) % 55; |
|---|
| 266 | ma[ii]=mk; |
|---|
| 267 | mk=mj-mk; |
|---|
| 268 | if (mk < MZ) mk += MBIG; |
|---|
| 269 | mj=ma[ii]; |
|---|
| 270 | } |
|---|
| 271 | for (k=1;k<=4;k++) |
|---|
| 272 | for (i=1;i<=55;i++) { |
|---|
| 273 | ma[i] -= ma[1+(i+30) % 55]; |
|---|
| 274 | if (ma[i] < MZ) ma[i] += MBIG; |
|---|
| 275 | } |
|---|
| 276 | inext=0; |
|---|
| 277 | inextp=31; |
|---|
| 278 | *idum=1; |
|---|
| 279 | } |
|---|
| 280 | if (++inext == 56) inext=1; |
|---|
| 281 | if (++inextp == 56) inextp=1; |
|---|
| 282 | mj=ma[inext]-ma[inextp]; |
|---|
| 283 | if (mj < MZ) mj += MBIG; |
|---|
| 284 | ma[inext]=mj; |
|---|
| 285 | return mj*FAC; |
|---|
| 286 | } |
|---|
| 287 | #undef MBIG |
|---|
| 288 | #undef MSEED |
|---|
| 289 | #undef MZ |
|---|
| 290 | #undef FAC |
|---|
| 291 | |
|---|
| 292 | //////////////////////// a complete copy of ran3() renamed ran3a() |
|---|
| 293 | #define MBIG2 1000000000 |
|---|
| 294 | #define MSEED2 161803398 |
|---|
| 295 | #define MZ2 0 |
|---|
| 296 | #define FAC2 (1.0/MBIG2) |
|---|
| 297 | |
|---|
| 298 | float ran3a(long *idum) |
|---|
| 299 | { |
|---|
| 300 | static int inext,inextp; |
|---|
| 301 | static long ma[56]; |
|---|
| 302 | static int iff=0; |
|---|
| 303 | long mj,mk; |
|---|
| 304 | int i,ii,k; |
|---|
| 305 | |
|---|
| 306 | if (*idum < 0 || iff == 0) { |
|---|
| 307 | iff=1; |
|---|
| 308 | mj=MSEED2-(*idum < 0 ? -*idum : *idum); |
|---|
| 309 | mj %= MBIG2; |
|---|
| 310 | ma[55]=mj; |
|---|
| 311 | mk=1; |
|---|
| 312 | for (i=1;i<=54;i++) { |
|---|
| 313 | ii=(21*i) % 55; |
|---|
| 314 | ma[ii]=mk; |
|---|
| 315 | mk=mj-mk; |
|---|
| 316 | if (mk < MZ2) mk += MBIG2; |
|---|
| 317 | mj=ma[ii]; |
|---|
| 318 | } |
|---|
| 319 | for (k=1;k<=4;k++) |
|---|
| 320 | for (i=1;i<=55;i++) { |
|---|
| 321 | ma[i] -= ma[1+(i+30) % 55]; |
|---|
| 322 | if (ma[i] < MZ2) ma[i] += MBIG2; |
|---|
| 323 | } |
|---|
| 324 | inext=0; |
|---|
| 325 | inextp=31; |
|---|
| 326 | *idum=1; |
|---|
| 327 | } |
|---|
| 328 | if (++inext == 56) inext=1; |
|---|
| 329 | if (++inextp == 56) inextp=1; |
|---|
| 330 | mj=ma[inext]-ma[inextp]; |
|---|
| 331 | if (mj < MZ2) mj += MBIG2; |
|---|
| 332 | ma[inext]=mj; |
|---|
| 333 | return mj*FAC2; |
|---|
| 334 | } |
|---|
| 335 | #undef MBIG2 |
|---|
| 336 | #undef MSEED2 |
|---|
| 337 | #undef MZ2 |
|---|
| 338 | #undef FAC2 |
|---|
| 339 | |
|---|
| 340 | |
|---|
| 341 | // returns the interpolated point value in xx[0,n-1] that has the value x |
|---|
| 342 | double locate_interp(double xx[], long n, double x) |
|---|
| 343 | { |
|---|
| 344 | unsigned long ju,jm,jl,j; |
|---|
| 345 | int ascnd; |
|---|
| 346 | double pt; |
|---|
| 347 | |
|---|
| 348 | // char buf[256]; |
|---|
| 349 | |
|---|
| 350 | jl=0; |
|---|
| 351 | ju=n-1; |
|---|
| 352 | ascnd=(xx[n-1] > xx[0]); |
|---|
| 353 | while (ju-jl > 1) { |
|---|
| 354 | jm=(ju+jl) >> 1; |
|---|
| 355 | if (x > xx[jm] == ascnd) |
|---|
| 356 | jl=jm; |
|---|
| 357 | else |
|---|
| 358 | ju=jm; |
|---|
| 359 | } |
|---|
| 360 | j=jl; // the point I want is between xx[j] and xx[j+1] |
|---|
| 361 | pt = (x- xx[j])/(xx[j+1] - xx[j]); //fractional distance, using linear interpolation |
|---|
| 362 | |
|---|
| 363 | // sprintf(buf, "x = %g, j= %ld, pt = %g\r",x,j,pt); |
|---|
| 364 | // XOPNotice(buf); |
|---|
| 365 | |
|---|
| 366 | return(pt+(double)j); |
|---|
| 367 | } |
|---|
| 368 | |
|---|
| 369 | |
|---|
| 370 | |
|---|
| 371 | |
|---|
| 372 | ///////////////////////////// |
|---|
| 373 | /* RegisterFunction() |
|---|
| 374 | |
|---|
| 375 | Igor calls this at startup time to find the address of the |
|---|
| 376 | XFUNCs added by this XOP. See XOP manual regarding "Direct XFUNCs". |
|---|
| 377 | */ |
|---|
| 378 | static long |
|---|
| 379 | RegisterFunction() |
|---|
| 380 | { |
|---|
| 381 | int funcIndex; |
|---|
| 382 | |
|---|
| 383 | funcIndex = GetXOPItem(0); // Which function is Igor asking about? |
|---|
| 384 | switch (funcIndex) { |
|---|
| 385 | case 0: // |
|---|
| 386 | return((long)Monte_SANSX); |
|---|
| 387 | break; |
|---|
| 388 | case 1: // |
|---|
| 389 | return((long)Monte_SANSX2); |
|---|
| 390 | break; |
|---|
| 391 | case 2: // |
|---|
| 392 | return((long)DebyeSpheresX); |
|---|
| 393 | break; |
|---|
| 394 | case 3: // |
|---|
| 395 | return((long)Monte_SANSX3); |
|---|
| 396 | break; |
|---|
| 397 | case 4: // |
|---|
| 398 | return((long)Monte_SANSX4); |
|---|
| 399 | break; |
|---|
| 400 | case 5: // |
|---|
| 401 | return((long)maxDistanceX); |
|---|
| 402 | break; |
|---|
| 403 | case 6: // |
|---|
| 404 | return((long)binDistanceX); |
|---|
| 405 | break; |
|---|
| 406 | case 7: // |
|---|
| 407 | return((long)SobolX); |
|---|
| 408 | break; |
|---|
| 409 | case 8: // |
|---|
| 410 | return((long)binSLDDistanceX); |
|---|
| 411 | break; |
|---|
| 412 | case 9: // |
|---|
| 413 | return((long)MetropolisX); |
|---|
| 414 | break; |
|---|
| 415 | } |
|---|
| 416 | return(NIL); |
|---|
| 417 | } |
|---|
| 418 | |
|---|
| 419 | /* XOPEntry() |
|---|
| 420 | |
|---|
| 421 | This is the entry point from the host application to the XOP for all messages after the |
|---|
| 422 | INIT message. |
|---|
| 423 | */ |
|---|
| 424 | static void |
|---|
| 425 | XOPEntry(void) |
|---|
| 426 | { |
|---|
| 427 | long result = 0; |
|---|
| 428 | |
|---|
| 429 | switch (GetXOPMessage()) { |
|---|
| 430 | case FUNCADDRS: |
|---|
| 431 | result = RegisterFunction(); |
|---|
| 432 | break; |
|---|
| 433 | } |
|---|
| 434 | SetXOPResult(result); |
|---|
| 435 | } |
|---|
| 436 | |
|---|
| 437 | /* main(ioRecHandle) |
|---|
| 438 | |
|---|
| 439 | This is the initial entry point at which the host application calls XOP. |
|---|
| 440 | The message sent by the host must be INIT. |
|---|
| 441 | main() does any necessary initialization and then sets the XOPEntry field of the |
|---|
| 442 | ioRecHandle to the address to be called for future messages. |
|---|
| 443 | |
|---|
| 444 | changed for TK6 30JAN2012 SRK |
|---|
| 445 | */ |
|---|
| 446 | HOST_IMPORT int |
|---|
| 447 | main(IORecHandle ioRecHandle) |
|---|
| 448 | { |
|---|
| 449 | XOPInit(ioRecHandle); // Do standard XOP initialization. |
|---|
| 450 | SetXOPEntry(XOPEntry); // Set entry point for future calls. |
|---|
| 451 | |
|---|
| 452 | if (igorVersion < 600) { // Requires Igor Pro 6.00 or later. |
|---|
| 453 | SetXOPResult(IGOR_OBSOLETE); |
|---|
| 454 | return EXIT_FAILURE; |
|---|
| 455 | } |
|---|
| 456 | |
|---|
| 457 | SetXOPResult(0); |
|---|
| 458 | return EXIT_SUCCESS; |
|---|
| 459 | |
|---|
| 460 | } |
|---|
