source: sans/XOP_Dev/MonteCarlo/MonteCarlo.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: 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.