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