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

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

Added wavelength distribution as part of the simulation

Fixed bug where pixel locations of (-1,0) row or column were incorrectly rounded towards zero, and ended up on the detector in row 0 or column 0, effectively doubling the count there.

Separated the 4 functions into 4 separate files, for ease in checking that they are all the same. Now the only difference in each file should be the RNG used. And they had better be different - I'm almost positive now that cross-calling of the same RNG by different threads is responsible for the crashing.

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