source: sans/XOP_Dev/MonteCarlo/MonteCarlo3.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: 17.0 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        //changed for TK6
61        int numDimensions;
62        CountInt 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 = (double*)WaveData(p->inputWaveH);
115        ran_dev = (double*)WaveData(p->ran_devH);
116        nt = (double*)WaveData(p->ntH);
117        j1 = (double*)WaveData(p->j1H);
118        j2 = (double*)WaveData(p->j2H);
119        nn = (double*)WaveData(p->nnH);
120        //      MC_linear_data = (double*)WaveData(p->MC_linear_dataH);
121        results = (double*)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 = ran3a(&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       
166        //scattering power and maximum qvalue to bin
167        //      zpow = .1               //scattering power, calculated below
168        qmax = 4.0*pi/wavelength;       //maximum Q to bin 1D data. (A-1) (not really used)
169        sigabs_0 = 0.0;         // ignore absorption cross section/wavelength [1/(cm A)]
170        n_index = 50;   // maximum number of scattering events per neutron
171        num_bins = 200;         //number of 1-D bins (not really used)
172       
173        //c       total SAS cross-section
174        //
175        zpow = sig_sas*thick;                   //since I now calculate the sig_sas from the model
176        sig_abs = sigabs_0 * wavelength;
177        sig_total = sig_abs + sig_sas + sig_incoh;
178        //      Print "The TOTAL XSECTION. (CM-1) is ",sig_total
179        //      Print "The TOTAL SAS XSECTION. (CM-1) is ",sig_sas
180        //      results[0] = sig_total;
181        //      results[1] = sig_sas;
182        //      RATIO = SIG_ABS / SIG_TOTAL
183        ratio = sig_incoh / sig_total;
184       
185        theta_max = wavelength*qmax/(2.0*pi);
186        //C     SET Theta-STEP SIZE.
187        dth = theta_max/num_bins;
188        //      Print "theta bin size = dth = ",dth
189       
190        //C     INITIALIZE COUNTERS.
191        n1 = 0;
192        n2 = 0;
193        n3 = 0;
194        NSingleIncoherent = 0;
195        NSingleCoherent = 0;
196        NDoubleCoherent = 0;
197        NMultipleScatter = 0;
198        NScatterEvents = 0;
199        NMultipleCoherent = 0;
200        NCoherentEvents = 0;
201       
202        isOn = 0;
203       
204        //C     MONITOR LOOP - looping over the number of incedent neutrons
205        //note that zz, is the z-position in the sample - NOT the scattering power
206        // NOW, start the loop, throwing neutrons at the sample.
207        do {
208                ////SpinProcess() IS A CALLBACK, and not good for Threading!
209                //              if ((n1 % 1000 == 0) && gCallSpinProcess && SpinProcess()) {            // Spins cursor and allows background processing.
210                //                              retVal = -1;                                                            // User aborted.
211                //                              break;
212                //              }
213               
214//              vx = 0.0;                       // Initialize direction vector.
215//              vy = 0.0;
216//              vz = 1.0;
217               
218                theta = 0.0;            //      Initialize scattering angle.
219                phi = 0.0;                      //      Intialize azimuthal angle.
220                n1 += 1;                        //      Increment total number neutrons counter.
221                done = 0;                       //      True when neutron is absorbed or when  scattered out of the sample.
222                index = 0;                      //      Set counter for number of scattering events.
223                zz = 0.0;                       //      Set entering dimension of sample.
224                incoherentEvent = 0;
225                coherentEvent = 0;
226
227                // pick point in source aperture area           
228                do      {                               //      Makes sure position is within circle.
229                        ran = ran3a(&seed);             //[0,1]
230                        souXX = 2.0*sourAp*(ran-0.5);           //X beam position of neutron entering sample.
231                        ran = ran3a(&seed);             //[0,1]
232                        souYY = 2.0*sourAp*(ran-0.5);           //Y beam position ...
233                        rr = sqrt(souXX*souXX+souYY*souYY);             //Radial position of neutron in incident beam.
234                } while(rr>sourAp);
235
236                // pick point in sample aperture
237                do      {                               //      Makes sure position is within circle.
238                        ran = ran3a(&seed);             //[0,1]
239                        xx = 2.0*r1*(ran-0.5);          //X beam position of neutron entering sample.
240                        ran = ran3a(&seed);             //[0,1]
241                        yy = 2.0*r1*(ran-0.5);          //Y beam position ...
242                        rr = sqrt(xx*xx+yy*yy);         //Radial position of neutron in incident beam.
243                } while(rr>r1);
244               
245                //pick the wavelength out of the wavelength spread, approximate as a gaussian
246                // from NR - pg 288. Needs random # from [0,1]. del is deltaLam/lam (as FWHM) and the
247                // 2.35 converts to a gaussian std dev.
248                do {
249                        v1=2.0*ran3a(&seed)-1.0;
250                        v2=2.0*ran3a(&seed)-1.0;
251                        rsq=v1*v1+v2*v2;
252                } while (rsq >= 1.0 || rsq == 0.0);
253                fac=sqrt(-2.0*log10(rsq)/rsq);          //be sure to use log10() here, to duplicate the Igor code
254               
255                //              gset=v1*fac             //technically, I'm throwing away one of the two values
256                currWavelength = (v2*fac)*deltaLam*wavelength/2.35 + wavelength;
257               
258                magn = sqrt((souXX - xx)*(souXX - xx) + (souYY - yy)*(souYY - yy) + ssd*ssd);
259                vx = (souXX - xx)/magn;         // Initialize direction vector.
260                vy = (souYY - yy)/magn;
261                vz = (ssd - 0.)/magn;
262               
263                do {   //Scattering Loop, will exit when "done" == 1
264                        // keep scattering multiple times until the neutron exits the sample
265                        ran = ran3a(&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 = ran3a(&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,ran3a(seed))*delta;
301                                               
302                                                q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran3a(&seed))*delta;
303                                                theta = q0/2.0/pi*currWavelength;               //SAS approximation
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 = ran3a(&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 = ran3a(&seed);             //[0,1]
323                                        theta = acos(2.0*ran-1);
324                                       
325                                        ran = ran3a(&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                                if( index != 0) {               //neutron was scattered, figure out where it went
347                                        theta_z = acos(vz);             // Angle (= 2theta) WITH respect to z axis.
348                                        testQ = 2.0*pi*sin(theta_z)/currWavelength;
349                                       
350                                        // pick a random phi angle, and see if it lands on the detector
351                                        // since the scattering is isotropic, I can safely pick a new, random value
352                                        // this would not be true if simulating anisotropic scattering.
353                                        testPhi = ran3a(&seed)*2.0*pi;
354                                       
355                                        // is it on the detector?       
356                                        FindPixel(testQ,testPhi,currWavelength,yg_d,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel);
357                                       
358                                        if(xPixel != -1 && yPixel != -1) {
359                                                isOn += 1;
360                                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
361                                                indices[0] = xPixel;
362                                                indices[1] = yPixel;
363                                                if (retVal = MDGetNumericWavePointValue(wavH, indices, value))
364                                                        return retVal;
365                                                value[0] += 1; // Real part
366                                                if (retVal = MDSetNumericWavePointValue(wavH, indices, value))
367                                                        return retVal;
368                                                //if(index==1)  // only the single scattering events
369                                                //dp = dp0 + xPixel + yPixel*numColumns;                //offset the pointer to the exact memory location
370                                                //*dp += 1;             //increment the value there
371                                                //endif
372                                        }
373                                       
374                                       
375                                        /*              is this causing me a problem since I'm not locking these? Probably not, since it crashes even if I comment these out... */
376                                        if(theta_z < theta_max) {
377                                                //Choose index for scattering angle array.
378                                                //IND = NINT(THETA_z/DTH + 0.4999999)
379                                                ind = (long)(theta_z/dth + 0.4999999);          //round is eqivalent to nint()
380                                                nt[ind] += 1;                   //Increment bin for angle.
381                                                //Increment angle array for single scattering events.
382                                                if (index == 1) {
383                                                        j1[ind] += 1;
384                                                }
385                                                //Increment angle array for double scattering events.
386                                                if (index == 2) {
387                                                        j2[ind] += 1;
388                                                }
389                                        }
390                                        /**/
391                                       
392                                        // increment all of the counters now since done==1 here and I'm sure to exit and get another neutron
393                                        NScatterEvents += index;                //total number of scattering events
394                                        if(index == 1 && incoherentEvent == 1) {
395                                                NSingleIncoherent += 1;
396                                        }
397                                        if(index == 1 && coherentEvent == 1) {
398                                                NSingleCoherent += 1;
399                                        }
400                                        if(index == 2 && coherentEvent == 1 && incoherentEvent == 0) {
401                                                NDoubleCoherent += 1;
402                                        }
403                                        if(index > 1) {
404                                                NMultipleScatter += 1;
405                                        }
406                                        if(coherentEvent >= 1 && incoherentEvent == 0) {
407                                                NCoherentEvents += 1;
408                                        }
409                                        if(coherentEvent > 1 && incoherentEvent == 0) {
410                                                NMultipleCoherent += 1;
411                                        }
412                                       
413                                } else {        // index was zero, neutron must be transmitted, so just increment the proper counters and data
414                                        isOn += 1;
415                                        nt[0] += 1;
416                                       
417                                        //figure out where it landed
418                                        theta_z = acos(vz);             // Angle (= 2theta) WITH respect to z axis.
419                                        testQ = 2.0*pi*sin(theta_z)/currWavelength;
420                                       
421                                        // pick a random phi angle, and see if it lands on the detector
422                                        // since the scattering is isotropic, I can safely pick a new, random value
423                                        // this would not be true if simulating anisotropic scattering.
424                                        testPhi = ran3a(&seed)*2.0*pi;
425                                       
426                                        // is it on the detector?       
427                                        FindPixel(testQ,testPhi,currWavelength,yg_d,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel);
428                                       
429                                        if(xPixel != -1 && yPixel != -1) {
430                                                isOn += 1;
431                                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
432                                                indices[0] = xPixel;
433                                                indices[1] = yPixel;
434                                                if (retVal = MDGetNumericWavePointValue(wavH, indices, value))
435                                                        return retVal;
436                                                value[0] += 1; // Real part
437                                                if (retVal = MDSetNumericWavePointValue(wavH, indices, value))
438                                                        return retVal;
439                                                //if(index==1)  // only the single scattering events
440                                                //dp = dp0 + xPixel + yPixel*numColumns;                //offset the pointer to the exact memory location
441                                                //*dp += 1;             //increment the value there
442                                                //endif
443                                        }
444                                }       
445                        }
446                } while (!done);
447        } while(n1 < imon);
448       
449        // assign the results to the wave
450       
451        MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
452        value[0] = (double)n1;
453        indices[0] = 0;
454        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
455                return retVal;
456        value[0] = (double)n2;
457        indices[0] = 1;
458        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
459                return retVal;
460        value[0] = (double)isOn;
461        indices[0] = 2;
462        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
463                return retVal;
464        value[0] = (double)NScatterEvents;
465        indices[0] = 3;
466        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
467                return retVal;
468        value[0] = (double)NSingleCoherent;
469        indices[0] = 4;
470        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
471                return retVal;
472        value[0] = (double)NMultipleCoherent;
473        indices[0] = 5;
474        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
475                return retVal;
476        value[0] = (double)NMultipleScatter;
477        indices[0] = 6;
478        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
479                return retVal; 
480        value[0] = (double)NCoherentEvents;
481        indices[0] = 7;
482        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
483                return retVal;
484       
485        //      WaveHandleModified(wavH);                       // Inform Igor that we have changed the wave. (CALLBACK! needed, but not allowed in Threading)
486       
487        return(0);
488}
Note: See TracBrowser for help on using the repository browser.