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