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

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

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

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

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