source: sans/XOP_Dev/MonteCarlo/MonteCarlo2.c @ 591

Last change on this file since 591 was 591, checked in by srkline, 13 years ago

Cosmetic changes to monteCarlo code. Added some ".0" here and there to make sure that all of the numerical values are really interpreted as floating point. Also changed the variable "result" to be "retVal" to avoid confusion with the wave of returned values "results"

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