source: sans/XOP_Dev/MonteCarlo/MonteCarlo.c @ 834

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

Changes to the XOP code to upgrade to ToolKit? v6. Changes are the ones outlined in the Appendix A of the TK6 manual. SOme of the XOP support routines and the #pragma for 2-byte structures have changed. Per Howard Rodstein, there is no need to change the c files to cpp. c should work and compile just fine.

These changes work correctly on my mac. Next is to make sure that they work correctly on the two build machines.

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