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