source: sans/XOP_Dev/MonteCarlo/MonteCarlo_Main.c

Last change on this file was 997, checked in by srkline, 7 years ago

Updates for Xcode 7, Toolkit 7, and 64-bit version

File size: 9.3 KB
Line 
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
23int MC_round(double x)
24{
25        if(x == 0)      {
26                return (int)(0);
27        }
28        if(x > 0)       {
29                return (int)(x + 0.5f);
30        } else {                // x < 0
31                return (int)(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
38int
39FindPixel(double testQ, double testPhi, double lam, double yg_d, double sdd,
40                  double pixSize, double xCtr, double yCtr, int *xPixel, int *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
97int
98NewDirection(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
133double
134path_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
154float ran1(SInt32 *idum)
155{
156        int j;
157        SInt32 k;
158        static SInt32 iy=0;
159        static SInt32 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
204float ran1a(SInt32 *idum)
205{
206        int j;
207        SInt32 k;
208        static SInt32 iy=0;
209        static SInt32 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
250float ran3(SInt32 *idum)
251{
252        static int inext,inextp;
253        static SInt32 ma[56];
254        static int iff=0;
255        SInt32 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
298float ran3a(SInt32 *idum)
299{
300        static int inext,inextp;
301        static SInt32 ma[56];
302        static int iff=0;
303        SInt32 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
342double locate_interp(double xx[], IndexInt n, double x)
343{
344        UInt32 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 */
378static XOPIORecResult
379RegisterFunction()
380{
381        int funcIndex;
382       
383        funcIndex = (int)GetXOPItem(0);         // Which function is Igor asking about?
384        switch (funcIndex) {
385                case 0:                                         //
386                        return((XOPIORecResult)Monte_SANSX);
387                        break;
388                case 1:                                         //
389                        return((XOPIORecResult)Monte_SANSX2);
390                        break;
391                case 2:                                         //
392                        return((XOPIORecResult)DebyeSpheresX);
393                        break;
394                case 3:                                         //
395                        return((XOPIORecResult)Monte_SANSX3);
396                        break;
397                case 4:                                         //
398                        return((XOPIORecResult)Monte_SANSX4);
399                        break;
400                case 5:                                         //
401                        return((XOPIORecResult)maxDistanceX);
402                        break;
403                case 6:                                         //
404                        return((XOPIORecResult)binDistanceX);
405                        break;
406                case 7:                                         //
407                        return((XOPIORecResult)SobolX);
408                        break;
409                case 8:                                         //
410                        return((XOPIORecResult)binSLDDistanceX);
411                        break;
412                case 9:                                         //
413                        return((XOPIORecResult)MetropolisX);
414                        break;
415        }
416        return(0);
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 updated for TK7 May2016
425 */
426//extern "C" void       //only use this if file is .cpp
427static void
428XOPEntry(void)
429{       
430    XOPIORecResult result = 0;
431       
432        switch (GetXOPMessage()) {
433                case FUNCADDRS:
434                        result = RegisterFunction();
435                        break;
436        }
437        SetXOPResult(result);
438}
439
440/*      main(ioRecHandle)
441 
442 This is the initial entry point at which the host application calls XOP.
443 The message sent by the host must be INIT.
444 main() does any necessary initialization and then sets the XOPEntry field of the
445 ioRecHandle to the address to be called for future messages.
446 
447 changed for TK6 30JAN2012 SRK
448 */
449HOST_IMPORT int
450XOPMain(IORecHandle ioRecHandle)
451{       
452        XOPInit(ioRecHandle);                                                   // Do standard XOP initialization.
453        SetXOPEntry(XOPEntry);                                                  // Set entry point for future calls.
454       
455        if (igorVersion < 620)  {                                       // Requires Igor Pro 6.00 or later.
456                SetXOPResult(OLD_IGOR);
457                return EXIT_FAILURE;
458        }
459
460        SetXOPResult(0);
461        return EXIT_SUCCESS;
462       
463}
Note: See TracBrowser for help on using the repository browser.