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