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