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