source: sans/XOP_Dev/MonteCarlo/MonteCarlo.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: 17.0 KB
Line 
1/*
2 *  MonteCarlo.c
3 *  SANSAnalysis
4 *
5 *  Created by Steve Kline on 10/16/08.
6 *  Copyright 2008 __MyCompanyName__. All rights reserved.
7 *
8 */
9
10//
11// uses ran3()
12//
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//static int gCallSpinProcess = 1;              // Set to 1 to all user abort (cmd dot) and background processing.
19
20
21
22int
23Monte_SANSX(MC_ParamsPtr p) {
24        double *inputWave;                              /* pointer to double precision wave data */
25        double *ran_dev;                                /* pointer to double precision wave data */
26        double *nt;                             /* pointer to double precision wave data */
27        double *j1;                             /* pointer to double precision wave data */
28        double *j2;                             /* pointer to double precision wave data */
29        double *nn;                             /* pointer to double precision wave data */
30//      double *MC_linear_data;                         /* pointer to double precision wave data */
31        double *results;                                /* pointer to double precision wave data */
32        double retVal;                          //return value
33
34        SInt64 imon;
35        double r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas;
36        IndexInt ind,index,n_index;
37        double qmax,theta_max,q0,zpow;
38        CountInt n1,n2,n3;
39        double dth,zz,xx,yy,phi;
40        double theta,ran,ll,rr;
41        int done,find_theta,err;                //used as logicals
42        int xPixel,yPixel;
43        double vx,vy,vz,theta_z;
44        double sig_abs,ratio,sig_total;
45        double testQ,testPhi,left,delta,dummy,pi;
46        double sigabs_0,num_bins;
47        CountInt NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent;
48        CountInt NDoubleCoherent,NMultipleScatter,isOn,xCtr_long,yCtr_long;
49        CountInt NMultipleCoherent,NCoherentEvents;
50        double deltaLam,v1,v2,currWavelength,rsq,fac;           //for simulating wavelength distribution
51        double ssd, sourAp, souXX, souYY, magn;         //source-to-sample, and source Ap radius for initlal trajectory
52        double vz_1,g,yg_d;                             //gravity terms
53       
54
55       
56        // for accessing the 2D wave data, direct method (see the WaveAccess example XOP)
57        waveHndl wavH;
58//      int waveType,hState;
59        //changed for TK6
60        int numDimensions;
61        CountInt dimensionSizes[MAX_DIMENSIONS+1];
62//      char* dataStartPtr;
63//      long dataOffset;
64//      long numRows, numColumns;
65        CountInt numRows_ran_dev;
66//      double *dp0, *dp;
67        double value[2];                                // Pointers used for double data.
68        SInt32 seed;
69        IndexInt indices[MAX_DIMENSIONS];
70       
71//      char buf[256];
72               
73        vz_1 = 3.956e5;         //velocity [cm/s] of 1 A neutron
74        g = 981.0;                              //gravity acceleration [cm/s^2]
75
76        /* check that wave handles are all valid */
77        if (p->inputWaveH == NIL) {
78                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
79                return(NON_EXISTENT_WAVE);
80        }
81        if (p->ran_devH == NIL) {
82                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
83                return(NON_EXISTENT_WAVE);
84        }       
85        if (p->ntH == NIL) {
86                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
87                return(NON_EXISTENT_WAVE);
88        }
89        if (p->j1H == NIL) {
90                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
91                return(NON_EXISTENT_WAVE);
92        }
93        if (p->j2H == NIL) {
94                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
95                return(NON_EXISTENT_WAVE);
96        }
97        if (p->nnH == NIL) {
98                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
99                return(NON_EXISTENT_WAVE);
100        }
101        if (p->MC_linear_dataH == NIL) {
102                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
103                return(NON_EXISTENT_WAVE);
104        }
105        if (p->resultsH == NIL) {
106                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
107                return(NON_EXISTENT_WAVE);
108        }
109       
110        p->retVal = 0.0;
111       
112// trusting that all inputs are double PRECISION WAVES!!!
113        inputWave = (double*)WaveData(p->inputWaveH);
114        ran_dev = (double*)WaveData(p->ran_devH);
115        nt = (double*)WaveData(p->ntH);
116        j1 = (double*)WaveData(p->j1H);
117        j2 = (double*)WaveData(p->j2H);
118        nn = (double*)WaveData(p->nnH);
119//      MC_linear_data = (double*)WaveData(p->MC_linear_dataH);
120        results = (double*)WaveData(p->resultsH);
121       
122        seed = (SInt32)results[0];
123       
124        //      sprintf(buf, "input seed = %ld\r", seed);
125        //      XOPNotice(buf);
126       
127        if(seed >= 0) {
128                seed = -1234509876;
129        }
130
131        dummy = ran3(&seed);            //initialize the random sequence by passing in a negative value
132        seed = 12348765;                //non-negative after that does nothing
133
134        imon = (SInt64)inputWave[0];
135        r1 = inputWave[1];
136        r2 = inputWave[2];
137        xCtr = inputWave[3];
138        yCtr = inputWave[4];
139        sdd = inputWave[5];
140        pixSize = inputWave[6];
141        thick = inputWave[7];
142        wavelength = inputWave[8];
143        sig_incoh = inputWave[9];
144        sig_sas = inputWave[10];
145        deltaLam = inputWave[11];
146        ssd = inputWave[12];                    // in cm, like SDD
147        sourAp = inputWave[13];         // radius, in cm, like r1 and r2       
148       
149        xCtr_long = (CountInt)(xCtr+0.5);       //round() not defined in VC8
150        yCtr_long = (CountInt)(yCtr+0.5);       // and is probably wrong to use anyways (returns double!)
151       
152        dummy = MDGetWaveScaling(p->ran_devH, 0, &delta, &left);                //0 is the rows
153        if (retVal = MDGetWaveDimensions(p->ran_devH, &numDimensions, dimensionSizes))
154                return retVal;
155        numRows_ran_dev = dimensionSizes[0];
156       
157        pi = 4.0*atan(1.0);     
158       
159        // access the 2D wave data for writing using the direct method
160        wavH = p->MC_linear_dataH;
161        if (wavH == NIL)
162                return NOWAV;
163       
164//scattering power and maximum qvalue to bin
165//      zpow = .1               //scattering power, calculated below
166        qmax = 4.0*pi/wavelength;       //maximum Q to bin 1D data. (A-1) (not really used)
167        sigabs_0 = 0.0;         // ignore absorption cross section/wavelength [1/(cm A)]
168        n_index = 50;   // maximum number of scattering events per neutron
169        num_bins = 200;         //number of 1-D bins (not really used)
170       
171//c       total SAS cross-section
172//
173        zpow = sig_sas*thick;                   //since I now calculate the sig_sas from the model
174        sig_abs = sigabs_0 * wavelength;
175        sig_total = sig_abs + sig_sas + sig_incoh;
176//      Print "The TOTAL XSECTION. (CM-1) is ",sig_total
177//      Print "The TOTAL SAS XSECTION. (CM-1) is ",sig_sas
178//      results[0] = sig_total;
179//      results[1] = sig_sas;
180//      RATIO = SIG_ABS / SIG_TOTAL
181        ratio = sig_incoh / sig_total;
182       
183        theta_max = wavelength*qmax/(2.0*pi);
184//C     SET Theta-STEP SIZE.
185        dth = theta_max/num_bins;
186//      Print "theta bin size = dth = ",dth
187
188//C     INITIALIZE COUNTERS.
189        n1 = 0;
190        n2 = 0;
191        n3 = 0;
192        NSingleIncoherent = 0;
193        NSingleCoherent = 0;
194        NDoubleCoherent = 0;
195        NMultipleScatter = 0;
196        NScatterEvents = 0;
197        NMultipleCoherent = 0;
198        NCoherentEvents = 0;
199
200        isOn = 0;
201       
202//C     MONITOR LOOP - looping over the number of incedent neutrons
203//note that zz, is the z-position in the sample - NOT the scattering power
204// NOW, start the loop, throwing neutrons at the sample.
205        do {
206                ////SpinProcess() IS A CALLBACK, and not good for Threading!
207//              if ((n1 % 1000 == 0) && gCallSpinProcess && SpinProcess()) {            // Spins cursor and allows background processing.
208//                              retVal = -1;                                                            // User aborted.
209//                              break;
210//              }
211       
212//              vx = 0.0;                       // Initialize direction vector.
213//              vy = 0.0;
214//              vz = 1.0;
215               
216                theta = 0.0;            //      Initialize scattering angle.
217                phi = 0.0;                      //      Intialize azimuthal angle.
218                n1 += 1;                        //      Increment total number neutrons counter.
219                done = 0;                       //      True when neutron is absorbed or when  scattered out of the sample.
220                index = 0;                      //      Set counter for number of scattering events.
221                zz = 0.0;                       //      Set entering dimension of sample.
222                incoherentEvent = 0;
223                coherentEvent = 0;
224               
225                // pick point in source aperture area           
226                do      {                               //      Makes sure position is within circle.
227                        ran = ran3(&seed);              //[0,1]
228                        souXX = 2.0*sourAp*(ran-0.5);           //X beam position of neutron entering sample.
229                        ran = ran3(&seed);              //[0,1]
230                        souYY = 2.0*sourAp*(ran-0.5);           //Y beam position ...
231                        rr = sqrt(souXX*souXX+souYY*souYY);             //Radial position of neutron in incident beam.
232                } while(rr>sourAp);
233                               
234                               
235                // pick point in sample aperture
236                do      {                               //      Makes sure position is within circle.
237                        ran = ran3(&seed);              //[0,1]
238                        xx = 2.0*r1*(ran-0.5);          //X beam position of neutron entering sample.
239                        ran = ran3(&seed);              //[0,1]
240                        yy = 2.0*r1*(ran-0.5);          //Y beam position ...
241                        rr = sqrt(xx*xx+yy*yy);         //Radial position of neutron in incident beam.
242                } while(rr>r1);
243
244                //pick the wavelength out of the wavelength spread, approximate as a gaussian
245                // from NR - pg 288. Needs random # from [0,1]. del is deltaLam/lam (as FWHM) and the
246                // 2.35 converts to a gaussian std dev.
247                do {
248                        v1=2.0*ran3(&seed)-1.0;
249                        v2=2.0*ran3(&seed)-1.0;
250                        rsq=v1*v1+v2*v2;
251                } while (rsq >= 1.0 || rsq == 0.0);
252                fac=sqrt(-2.0*log10(rsq)/rsq);          //be sure to use log10() here, to duplicate the Igor code
253                               
254                //              gset=v1*fac             //technically, I'm throwing away one of the two values
255                currWavelength = (v2*fac)*deltaLam*wavelength/2.35 + wavelength;
256               
257                magn = sqrt((souXX - xx)*(souXX - xx) + (souYY - yy)*(souYY - yy) + ssd*ssd);
258                vx = (souXX - xx)/magn;         // Initialize direction vector.
259                vy = (souYY - yy)/magn;
260                vz = (ssd - 0.)/magn;
261               
262               
263                do {   //Scattering Loop, will exit when "done" == 1
264                                // keep scattering multiple times until the neutron exits the sample
265                        ran = ran3(&seed);              //[0,1]  RANDOM NUMBER FOR DETERMINING PATH LENGTH
266                        ll = path_len(ran,sig_total);
267                        //Determine new scattering direction vector.
268                        err = NewDirection(&vx,&vy,&vz,theta,phi);              //vx,vy,vz updated, theta, phi unchanged by function
269                                                                       
270                        //X,Y,Z-POSITION OF SCATTERING EVENT.
271                        xx += ll*vx;
272                        yy += ll*vy;
273                        zz += ll*vz;
274                        rr = sqrt(xx*xx+yy*yy);         //radial position of scattering event.
275
276                        //sprintf(buf, "xx,yy,zz,vx,vy,vz,ll = %g %g %g %g %g %g %g\r",xx,yy,zz,vx,vy,vz,ll);
277                        //XOPNotice(buf);
278                                               
279                        //Check whether interaction occurred within sample volume.
280                        if (((zz > 0.0) && (zz < thick)) && (rr < r2)) {
281                                //NEUTRON INTERACTED.
282                                //sprintf(buf,"neutron interacted\r");
283                                //XOPNotice(buf);
284                               
285                                index += 1;                     //Increment counter of scattering events.
286                                if (index == 1) {
287                                        n2 += 1;                //Increment # of scat. neutrons
288                                }
289                                ran = ran3(&seed);              //[0,1]
290                                //Split neutron interactions into scattering and absorption events
291                                if (ran > ratio ) {             //C             NEUTRON SCATTERED coherently
292                                        //sprintf(buf,"neutron scatters coherently\r");
293                                        //XOPNotice(buf);
294                                        coherentEvent += 1;
295                                        find_theta = 0;                 //false
296                                        do {
297                                                // pick a q-value from the deviate function
298                                                // pnt2x truncates the point to an integer before returning the x
299                                                // so get it from the wave scaling instead
300//                                              q0 =left + binarysearchinterp(ran_dev,ran1(seed))*delta;
301                                               
302                                                q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran3(&seed))*delta;
303                                                theta = q0/2.0/pi*currWavelength;               //SAS approximation. 1% error at theta=30 degrees
304                                               
305                                                find_theta = 1;         //always accept
306
307                                                //sprintf(buf, "after locate_interp call q0 = %g, theta = %g,left = %g,delta = %g\r",q0,theta,left,delta);
308                                                //XOPNotice(buf);
309
310                                        } while(!find_theta);
311                                       
312                                        ran = ran3(&seed);              //[0,1]
313                                        phi = 2.0*pi*ran;                       //Chooses azimuthal scattering angle.
314                                } else {
315                                        //NEUTRON scattered incoherently
316                                        //sprintf(buf,"neutron scatters incoherent\r");
317                                        //XOPNotice(buf);
318                                        incoherentEvent += 1;
319                                  // phi and theta are random over the entire sphere of scattering
320                                        // !can't just choose random theta and phi, won't be random over sphere solid angle
321
322                                        ran = ran3(&seed);              //[0,1]
323                                        theta = acos(2.0*ran-1);
324                       
325                                        ran = ran3(&seed);              //[0,1]
326                                        phi = 2.0*pi*ran;                       //Chooses azimuthal scattering angle.
327                                }               //(ran > ratio)
328                        } else {
329                                //NEUTRON ESCAPES FROM SAMPLE -- bin it somewhere                                                               
330                                done = 1;               //done = true, will exit from loop
331                                //Increment #scattering events array
332                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
333                                indices[0] =index;                      //this sets access to nn[index]
334                                if (index <= n_index) {
335                                        if (retVal = MDGetNumericWavePointValue(p->nnH, indices, value))
336                                                return retVal;
337                                        value[0] += 1; // add one to the value
338                                        if (retVal = MDSetNumericWavePointValue(p->nnH, indices, value))
339                                                return retVal;
340                                //      nn[index] += 1;
341                                }
342                               
343                                // calculate fall due to gravity (in cm) (note that it is negative)
344                                yg_d = -0.5*g*sdd*(ssd+sdd)*(currWavelength/vz_1)*(currWavelength/vz_1);
345                               
346                               
347                                if( index != 0) {               //neutron was scattered, figure out where it went
348                                        theta_z = acos(vz);             // Angle (= 2theta) WITH respect to z axis.
349                                        testQ = 2.0*pi*sin(theta_z)/currWavelength;
350                                       
351                                        // pick a random phi angle, and see if it lands on the detector
352                                        // since the scattering is isotropic, I can safely pick a new, random value
353                                        // this would not be true if simulating anisotropic scattering.
354                                        testPhi = ran3(&seed)*2.0*pi;
355                                       
356                                        // is it on the detector?       
357                                        FindPixel(testQ,testPhi,currWavelength,yg_d,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel);
358                                                                                                       
359                                        if(xPixel != -1 && yPixel != -1) {
360                                                isOn += 1;
361                                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
362                                                indices[0] = xPixel;
363                                                indices[1] = yPixel;
364                                                if (retVal = MDGetNumericWavePointValue(wavH, indices, value))
365                                                        return retVal;
366                                                value[0] += 1; // Real part
367                                                if (retVal = MDSetNumericWavePointValue(wavH, indices, value))
368                                                        return retVal;
369                                                //if(index==1)  // only the single scattering events
370                                                        //dp = dp0 + xPixel + yPixel*numColumns;                //offset the pointer to the exact memory location
371                                                        //*dp += 1;             //increment the value there
372                                                //endif
373                                        }
374                                       
375
376        /*              is this causing me a problem since I'm not locking these? Probably not, since it crashes even if I comment these out... */
377                                        if(theta_z < theta_max) {
378                                                //Choose index for scattering angle array.
379                                                //IND = NINT(THETA_z/DTH + 0.4999999)
380                                                ind = (IndexInt)(theta_z/dth + 0.4999999);              //round is eqivalent to nint()
381                                                nt[ind] += 1;                   //Increment bin for angle.
382                                                //Increment angle array for single scattering events.
383                                                if (index == 1) {
384                                                        j1[ind] += 1;
385                                                }
386                                                //Increment angle array for double scattering events.
387                                                if (index == 2) {
388                                                        j2[ind] += 1;
389                                                }
390                                        }
391        /**/
392                                       
393                                        // increment all of the counters now since done==1 here and I'm sure to exit and get another neutron
394                                        NScatterEvents += index;                //total number of scattering events
395                                        if(index == 1 && incoherentEvent == 1) {
396                                                NSingleIncoherent += 1;
397                                        }
398                                        if(index == 1 && coherentEvent == 1) {
399                                                NSingleCoherent += 1;
400                                        }
401                                        if(index == 2 && coherentEvent == 1 && incoherentEvent == 0) {
402                                                NDoubleCoherent += 1;
403                                        }
404                                        if(index > 1) {
405                                                NMultipleScatter += 1;
406                                        }
407                                        if(coherentEvent >= 1 && incoherentEvent == 0) {
408                                                NCoherentEvents += 1;
409                                        }
410                                        if(coherentEvent > 1 && incoherentEvent == 0) {
411                                                NMultipleCoherent += 1;
412                                        }
413
414                                } else {        // index was zero, neutron must be transmitted, so just increment the proper counters and data
415                                                isOn += 1;
416                                                nt[0] += 1;
417                                       
418                                                //figure out where it landed
419                                                theta_z = acos(vz);             // Angle (= 2theta) WITH respect to z axis.
420                                                testQ = 2.0*pi*sin(theta_z)/currWavelength;
421                                               
422                                                // pick a random phi angle, and see if it lands on the detector
423                                                // since the scattering is isotropic, I can safely pick a new, random value
424                                                // this would not be true if simulating anisotropic scattering.
425                                                testPhi = ran3(&seed)*2.0*pi;
426                                               
427                                                // is it on the detector?       
428                                                FindPixel(testQ,testPhi,currWavelength,yg_d,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel);
429                                               
430                                                if(xPixel != -1 && yPixel != -1) {
431                                                        isOn += 1;
432                                                        MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
433                                                        indices[0] = xPixel;
434                                                        indices[1] = yPixel;
435                                                        if (retVal = MDGetNumericWavePointValue(wavH, indices, value))
436                                                                return retVal;
437                                                        value[0] += 1; // Real part
438                                                        if (retVal = MDSetNumericWavePointValue(wavH, indices, value))
439                                                                return retVal;
440                                                        //if(index==1)  // only the single scattering events
441                                                        //dp = dp0 + xPixel + yPixel*numColumns;                //offset the pointer to the exact memory location
442                                                        //*dp += 1;             //increment the value there
443                                                        //endif
444                                                }
445                                               
446                                        }       
447                        }
448                 } while (!done);
449        } while(n1 < imon);
450       
451// assign the results to the wave
452
453        MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
454        value[0] = (double)n1;
455        indices[0] = 0;
456        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
457                return retVal;
458        value[0] = (double)n2;
459        indices[0] = 1;
460        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
461                return retVal;
462        value[0] = (double)isOn;
463        indices[0] = 2;
464        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
465                return retVal;
466        value[0] = (double)NScatterEvents;
467        indices[0] = 3;
468        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
469                return retVal;
470        value[0] = (double)NSingleCoherent;
471        indices[0] = 4;
472        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
473                return retVal;
474        value[0] = (double)NMultipleCoherent;
475        indices[0] = 5;
476        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
477                return retVal;
478        value[0] = (double)NMultipleScatter;
479        indices[0] = 6;
480        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
481                return retVal; 
482        value[0] = (double)NCoherentEvents;
483        indices[0] = 7;
484        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
485                return retVal; 
486
487
488//      WaveHandleModified(wavH);                       // Inform Igor that we have changed the wave. (CALLBACK! needed, but not allowed in Threading)
489       
490        return(0);
491}
492////////        end of main function for calculating multiple scattering
Note: See TracBrowser for help on using the repository browser.