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