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