source: sans/XOP_Dev/MonteCarlo/MonteCarlo_Main.c @ 790

Last change on this file since 790 was 789, checked in by srkline, 12 years ago

Added XOP-ized functions for doing the distance-binned Debye Spheres calculation. Two versions are added, one that includes the SLDs and one that does not.

Also added a pseudo-random sequence generator, SobolX, that will generate pseudo-random 2D or 3D distributions. This fills space more uniformly than a random generation.

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