source: sans/XOP_Dev/MonteCarlo/MonteCarlo3.c

Last change on this file was 997, checked in by srkline, 7 years ago

Updates for Xcode 7, Toolkit 7, and 64-bit version

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