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

Last change on this file since 812 was 812, checked in by srkline, 11 years ago

Added source aperture size and SSD to simulation to get realistic initial trajectory, rather than only along the z-axis.

Added gravity fall to properly account for fall of neutrons and its effect on the primary transmitted beam and on the scattering pattern.

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