source: sans/XOP_Dev/MonteCarlo/MonteCarlo2.c @ 825

Last change on this file since 825 was 825, checked in by srkline, 11 years ago

Added new DebyeSphere? functions to the Win.rc file

Changed the MonteCarlo?.c files to do all of the declarations before any assignments. This was generating a heap of errors in VC10, but no problems in Xcode.

XOPs all seem to compile fine on the newly set up sansvm2.

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