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