source: sans/XOP_Dev/MonteCarlo/MonteCarlo2.c @ 623

Last change on this file since 623 was 623, checked in by srkline, 13 years ago

Updated the MonteCarlo? code to allow 4 processors, but simply copying the function 4 times, and defining 4 different random number generators. Still can't figure out what the problem is with threading a single version, but not worth the effort. Copy/paste is way faster.

Also added some simple (non-optimized) calculations for using Debye's sphere method. These are largely undocumented at this point - so see the code. These are XOP versions of the old ipf code I've used in the past, and stripped of the now-obsolete AltiVec? code (I now lose the 4x speedup from the vectorization...)

File size: 44.7 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#include "XOPStandardHeaders.h"                 // Include ANSI headers, Mac headers, IgorXOP.h, XOP.h and XOPSupport.h
12#include "MonteCarlo.h"
13
14static 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,ttot;
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
54
55        // for accessing the 2D wave data, direct method (see the WaveAccess example XOP)
56        waveHndl wavH;
57        int waveType,hState;
58        long numDimensions;
59        long dimensionSizes[MAX_DIMENSIONS+1];
60        char* dataStartPtr;
61        long dataOffset;
62        long numRows, numColumns,numRows_ran_dev;
63        double *dp0, *dp, value[2];                             // Pointers used for double data.
64        long seed;
65        long indices[MAX_DIMENSIONS];
66       
67        char buf[256];
68               
69        /* check that wave handles are all valid */
70        if (p->inputWaveH == NIL) {
71                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
72                return(NON_EXISTENT_WAVE);
73        }
74        if (p->ran_devH == NIL) {
75                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
76                return(NON_EXISTENT_WAVE);
77        }       
78        if (p->ntH == NIL) {
79                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
80                return(NON_EXISTENT_WAVE);
81        }
82        if (p->j1H == NIL) {
83                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
84                return(NON_EXISTENT_WAVE);
85        }
86        if (p->j2H == NIL) {
87                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
88                return(NON_EXISTENT_WAVE);
89        }
90        if (p->nnH == NIL) {
91                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
92                return(NON_EXISTENT_WAVE);
93        }
94        if (p->MC_linear_dataH == NIL) {
95                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
96                return(NON_EXISTENT_WAVE);
97        }
98        if (p->resultsH == NIL) {
99                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
100                return(NON_EXISTENT_WAVE);
101        }
102       
103        p->retVal = 0;
104       
105// trusting that all inputs are DOUBLE PRECISION WAVES!!!
106        inputWave = WaveData(p->inputWaveH);
107        ran_dev = WaveData(p->ran_devH);
108        nt = WaveData(p->ntH);
109        j1 = WaveData(p->j1H);
110        j2 = WaveData(p->j2H);
111        nn = WaveData(p->nnH);
112//      MC_linear_data = WaveData(p->MC_linear_dataH);
113        results = WaveData(p->resultsH);
114       
115        seed = (long)results[0];
116       
117//      sprintf(buf, "input seed = %ld\r", seed);
118//      XOPNotice(buf);
119       
120        if(seed >= 0) {
121                seed = -1234509876;
122        }
123
124        dummy = ran1(&seed);            //initialize the random sequence by passing in a negative value
125        seed = 12348765;                //non-negative after that does nothing
126
127        imon = (int)inputWave[0];
128        r1 = inputWave[1];
129        r2 = inputWave[2];
130        xCtr = inputWave[3];
131        yCtr = inputWave[4];
132        sdd = inputWave[5];
133        pixSize = inputWave[6];
134        thick = inputWave[7];
135        wavelength = inputWave[8];
136        sig_incoh = inputWave[9];
137        sig_sas = inputWave[10];
138        xCtr_long = round(xCtr);
139        yCtr_long = round(yCtr);
140       
141        dummy = MDGetWaveScaling(p->ran_devH, 0, &delta, &left);                //0 is the rows
142        if (retVal = MDGetWaveDimensions(p->ran_devH, &numDimensions, dimensionSizes))
143                return retVal;
144        numRows_ran_dev = dimensionSizes[0];
145       
146        pi = 4.0*atan(1.0);     
147       
148        // access the 2D wave data for writing using the direct method
149        wavH = p->MC_linear_dataH;
150        if (wavH == NIL)
151                return NOWAV;
152
153//      waveType = WaveType(wavH);
154//      if (waveType & NT_CMPLX)
155//              return NO_COMPLEX_WAVE;
156//      if (waveType==TEXT_WAVE_TYPE)
157//              return NUMERIC_ACCESS_ON_TEXT_WAVE;
158//      if (retVal = MDGetWaveDimensions(wavH, &numDimensions, dimensionSizes))
159//              return retVal;
160//      numRows = dimensionSizes[0];
161//      numColumns = dimensionSizes[1];
162       
163//      if (retVal = MDAccessNumericWaveData(wavH, kMDWaveAccessMode0, &dataOffset))
164//              return retVal;
165               
166//      hState = MoveLockHandle(wavH);          // So wave data can't move. Remember to call HSetState when done.
167//      dataStartPtr = (char*)(*wavH) + dataOffset;
168//      dp0 = (double*)dataStartPtr;                    // Pointer to the start of the 2D wave data.
169       
170//scattering power and maximum qvalue to bin
171//      zpow = .1               //scattering power, calculated below
172        qmax = 4.0*pi/wavelength;       //maximum Q to bin 1D data. (A-1) (not really used)
173        sigabs_0 = 0.0;         // ignore absorption cross section/wavelength [1/(cm A)]
174        n_index = 50;   // maximum number of scattering events per neutron
175        num_bins = 200;         //number of 1-D bins (not really used)
176       
177//c       total SAS cross-section
178//
179        zpow = sig_sas*thick;                   //since I now calculate the sig_sas from the model
180        sig_abs = sigabs_0 * wavelength;
181        sig_total = sig_abs + sig_sas + sig_incoh;
182//      Print "The TOTAL XSECTION. (CM-1) is ",sig_total
183//      Print "The TOTAL SAS XSECTION. (CM-1) is ",sig_sas
184//      results[0] = sig_total;
185//      results[1] = sig_sas;
186//      RATIO = SIG_ABS / SIG_TOTAL
187        ratio = sig_incoh / sig_total;
188       
189        theta_max = wavelength*qmax/(2*pi);
190//C     SET Theta-STEP SIZE.
191        dth = theta_max/num_bins;
192//      Print "theta bin size = dth = ",dth
193
194//C     INITIALIZE COUNTERS.
195        n1 = 0;
196        n2 = 0;
197        n3 = 0;
198        NSingleIncoherent = 0;
199        NSingleCoherent = 0;
200        NDoubleCoherent = 0;
201        NMultipleScatter = 0;
202        NScatterEvents = 0;
203        NMultipleCoherent = 0;
204        NCoherentEvents = 0;
205       
206        isOn = 0;
207       
208//C     MONITOR LOOP - looping over the number of incedent neutrons
209//note that zz, is the z-position in the sample - NOT the scattering power
210// NOW, start the loop, throwing neutrons at the sample.
211        do {
212                ////SpinProcess() IS A CALLBACK, and not good for Threading!
213//              if ((n1 % 1000 == 0) && gCallSpinProcess && SpinProcess()) {            // Spins cursor and allows background processing.
214//                              retVal = -1;                                                            // User aborted.
215//                              break;
216//              }
217       
218                vx = 0.0;                       // Initialize direction vector.
219                vy = 0.0;
220                vz = 1.0;
221               
222                theta = 0.0;            //      Initialize scattering angle.
223                phi = 0.0;                      //      Intialize azimuthal angle.
224                n1 += 1;                        //      Increment total number neutrons counter.
225                done = 0;                       //      True when neutron is absorbed or when  scattered out of the sample.
226                index = 0;                      //      Set counter for number of scattering events.
227                zz = 0.0;                       //      Set entering dimension of sample.
228                incoherentEvent = 0;
229                coherentEvent = 0;
230               
231                do      {                               //      Makes sure position is within circle.
232                        ran = ran1(&seed);              //[0,1]
233                        xx = 2.0*r1*(ran-0.5);          //X beam position of neutron entering sample.
234                        ran = ran1(&seed);              //[0,1]
235                        yy = 2.0*r1*(ran-0.5);          //Y beam position ...
236                        rr = sqrt(xx*xx+yy*yy);         //Radial position of neutron in incident beam.
237                } while(rr>r1);
238
239                do {   //Scattering Loop, will exit when "done" == 1
240                                // keep scattering multiple times until the neutron exits the sample
241                        ran = ran1(&seed);              //[0,1]  RANDOM NUMBER FOR DETERMINING PATH LENGTH
242                        ll = path_len(ran,sig_total);
243                        //Determine new scattering direction vector.
244                        err = NewDirection(&vx,&vy,&vz,theta,phi);              //vx,vy,vz updated, theta, phi unchanged by function
245                                                                       
246                        //X,Y,Z-POSITION OF SCATTERING EVENT.
247                        xx += ll*vx;
248                        yy += ll*vy;
249                        zz += ll*vz;
250                        rr = sqrt(xx*xx+yy*yy);         //radial position of scattering event.
251
252                        //sprintf(buf, "xx,yy,zz,vx,vy,vz,ll = %g %g %g %g %g %g %g\r",xx,yy,zz,vx,vy,vz,ll);
253                        //XOPNotice(buf);
254                                               
255                        //Check whether interaction occurred within sample volume.
256                        if (((zz > 0.0) && (zz < thick)) && (rr < r2)) {
257                                //NEUTRON INTERACTED.
258                                //sprintf(buf,"neutron interacted\r");
259                                //XOPNotice(buf);
260                               
261                                index += 1;                     //Increment counter of scattering events.
262                                if (index == 1) {
263                                        n2 += 1;                //Increment # of scat. neutrons
264                                }
265                                ran = ran1(&seed);              //[0,1]
266                                //Split neutron interactions into scattering and absorption events
267                                if (ran > ratio ) {             //C             NEUTRON SCATTERED coherently
268                                        //sprintf(buf,"neutron scatters coherently\r");
269                                        //XOPNotice(buf);
270                                        coherentEvent += 1;
271                                        find_theta = 0;                 //false
272                                        do {
273                                                // pick a q-value from the deviate function
274                                                // pnt2x truncates the point to an integer before returning the x
275                                                // so get it from the wave scaling instead
276//                                              q0 =left + binarysearchinterp(ran_dev,ran1(seed))*delta;
277                                               
278                                                q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran1(&seed))*delta;
279                                                theta = q0/2/pi*wavelength;             //SAS approximation
280                                               
281                                                find_theta = 1;         //always accept
282
283                                                //sprintf(buf, "after locate_interp call q0 = %g, theta = %g,left = %g,delta = %g\r",q0,theta,left,delta);
284                                                //XOPNotice(buf);
285
286                                        } while(!find_theta);
287                                       
288                                        ran = ran1(&seed);              //[0,1]
289                                        phi = 2.0*pi*ran;                       //Chooses azimuthal scattering angle.
290                                } else {
291                                        //NEUTRON scattered incoherently
292                                        //sprintf(buf,"neutron scatters incoherent\r");
293                                        //XOPNotice(buf);
294                                        incoherentEvent += 1;
295                                  // phi and theta are random over the entire sphere of scattering
296                                        // !can't just choose random theta and phi, won't be random over sphere solid angle
297
298                                        ran = ran1(&seed);              //[0,1]
299                                        theta = acos(2.0*ran-1);
300                       
301                                        ran = ran1(&seed);              //[0,1]
302                                        phi = 2.0*pi*ran;                       //Chooses azimuthal scattering angle.
303                                }               //(ran > ratio)
304                        } else {
305                                //NEUTRON ESCAPES FROM SAMPLE -- bin it somewhere                                                               
306                                done = 1;               //done = true, will exit from loop
307                                //Increment #scattering events array
308                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
309                                indices[0] =index;                      //this sets access to nn[index]
310                                if (index <= n_index) {
311                                        if (retVal = MDGetNumericWavePointValue(p->nnH, indices, value))
312                                                return retVal;
313                                        value[0] += 1; // add one to the value
314                                        if (retVal = MDSetNumericWavePointValue(p->nnH, indices, value))
315                                                return retVal;
316                                //      nn[index] += 1;
317                                }
318                                                                                               
319                                if( index != 0) {               //neutron was scattered, figure out where it went
320                                        theta_z = acos(vz);             // Angle (= 2theta) WITH respect to z axis.
321                                        testQ = 2*pi*sin(theta_z)/wavelength;
322                                       
323                                        // pick a random phi angle, and see if it lands on the detector
324                                        // since the scattering is isotropic, I can safely pick a new, random value
325                                        // this would not be true if simulating anisotropic scattering.
326                                        testPhi = ran1(&seed)*2*pi;
327                                       
328                                        // is it on the detector?       
329                                        FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel);
330                                                                                                       
331                                        if(xPixel != -1 && yPixel != -1) {
332                                                isOn += 1;
333                                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
334                                                indices[0] = xPixel;
335                                                indices[1] = yPixel;
336                                                if (retVal = MDGetNumericWavePointValue(wavH, indices, value))
337                                                        return retVal;
338                                                value[0] += 1; // Real part
339                                                if (retVal = MDSetNumericWavePointValue(wavH, indices, value))
340                                                        return retVal;
341                                                //if(index==1)  // only the single scattering events
342                                                        //dp = dp0 + xPixel + yPixel*numColumns;                //offset the pointer to the exact memory location
343                                                        //*dp += 1;             //increment the value there
344                                                //endif
345                                        }
346                                       
347
348        /*              is this causing me a problem since I'm not locking these? Probably not, since it crashes even if I comment these out... */
349                                        if(theta_z < theta_max) {
350                                                //Choose index for scattering angle array.
351                                                //IND = NINT(THETA_z/DTH + 0.4999999)
352                                                ind = round(theta_z/dth + 0.4999999);           //round is eqivalent to nint()
353                                                nt[ind] += 1;                   //Increment bin for angle.
354                                                //Increment angle array for single scattering events.
355                                                if (index == 1) {
356                                                        j1[ind] += 1;
357                                                }
358                                                //Increment angle array for double scattering events.
359                                                if (index == 2) {
360                                                        j2[ind] += 1;
361                                                }
362                                        }
363        /**/
364                                       
365                                        // increment all of the counters now since done==1 here and I'm sure to exit and get another neutron
366                                        NScatterEvents += index;                //total number of scattering events
367                                        if(index == 1 && incoherentEvent == 1) {
368                                                NSingleIncoherent += 1;
369                                        }
370                                        if(index == 1 && coherentEvent == 1) {
371                                                NSingleCoherent += 1;
372                                        }
373                                        if(index == 2 && coherentEvent == 1 && incoherentEvent == 0) {
374                                                NDoubleCoherent += 1;
375                                        }
376                                        if(index > 1) {
377                                                NMultipleScatter += 1;
378                                        }
379                                        if(coherentEvent >= 1 && incoherentEvent == 0) {
380                                                NCoherentEvents += 1;
381                                        }
382                                        if(coherentEvent > 1 && incoherentEvent == 0) {
383                                                NMultipleCoherent += 1;
384                                        }
385
386                                } else {        // index was zero, neutron must be transmitted, so just increment the proper counters and data
387                                                isOn += 1;
388                                                nt[0] += 1;
389                                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
390                                                //indices[0] = xCtr_long;               //don't put everything in one pixel
391                                                //indices[1] = yCtr_long;
392                                                indices[0] = (long)round(xCtr+xx/pixSize);
393                                                indices[1] = (long)round(yCtr+yy/pixSize);
394                                                // check for valid indices - got an XOP error, probably from here
395                                                if(indices[0] > 127) indices[0] = 127;
396                                                if(indices[0] < 0) indices[0] = 0;
397                                                if(indices[1] > 127) indices[1] = 127;
398                                                if(indices[1] < 0) indices[1] = 0;
399                                               
400                                                if (retVal = MDGetNumericWavePointValue(wavH, indices, value))
401                                                        return retVal;
402                                                value[0] += 1; // Real part
403                                                if (retVal = MDSetNumericWavePointValue(wavH, indices, value))
404                                                        return retVal;
405                                        }       
406                        }
407                 } while (!done);
408        } while(n1 < imon);
409       
410// assign the results to the wave
411
412        MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
413        value[0] = (double)n1;
414        indices[0] = 0;
415        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
416                return retVal;
417        value[0] = (double)n2;
418        indices[0] = 1;
419        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
420                return retVal;
421        value[0] = (double)isOn;
422        indices[0] = 2;
423        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
424                return retVal;
425        value[0] = (double)NScatterEvents;
426        indices[0] = 3;
427        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
428                return retVal;
429        value[0] = (double)NSingleCoherent;
430        indices[0] = 4;
431        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
432                return retVal;
433        value[0] = (double)NMultipleCoherent;
434        indices[0] = 5;
435        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
436                return retVal;
437        value[0] = (double)NMultipleScatter;
438        indices[0] = 6;
439        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
440                return retVal; 
441        value[0] = (double)NCoherentEvents;
442        indices[0] = 7;
443        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
444                return retVal;
445
446//      HSetState((Handle)wavH, hState);                //release the handle of the 2D data wave
447//      WaveHandleModified(wavH);                       // Inform Igor that we have changed the wave. (CALLBACK! needed, but not allowed in Threading)
448       
449        return(0);
450}
451////////        end of X2
452
453
454
455////////////////                X3 using ran3a
456int
457Monte_SANSX3(MC_ParamsPtr p) {
458        double *inputWave;                              /* pointer to double precision wave data */
459        double *ran_dev;                                /* pointer to double precision wave data */
460        double *nt;                             /* pointer to double precision wave data */
461        double *j1;                             /* pointer to double precision wave data */
462        double *j2;                             /* pointer to double precision wave data */
463        double *nn;                             /* pointer to double precision wave data */
464//      double *MC_linear_data;                         /* pointer to double precision wave data */
465        double *results;                                /* pointer to double precision wave data */
466        double retVal;                          //return value
467
468        long imon;
469        double r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas;
470        long ind,index,n_index;
471        double qmax,theta_max,q0,zpow;
472        long n1,n2,n3;
473        double dth,zz,xx,yy,phi;
474        double theta,ran,ll,rr,ttot;
475        long done,find_theta,err;               //used as logicals
476        long xPixel,yPixel;
477        double vx,vy,vz,theta_z;
478        double sig_abs,ratio,sig_total;
479        double testQ,testPhi,left,delta,dummy,pi;
480        double sigabs_0,num_bins;
481        long NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent;
482        long NDoubleCoherent,NMultipleScatter,isOn,xCtr_long,yCtr_long;
483        long NMultipleCoherent,NCoherentEvents;
484
485
486        // for accessing the 2D wave data, direct method (see the WaveAccess example XOP)
487        waveHndl wavH;
488        int waveType,hState;
489        long numDimensions;
490        long dimensionSizes[MAX_DIMENSIONS+1];
491        char* dataStartPtr;
492        long dataOffset;
493        long numRows, numColumns,numRows_ran_dev;
494        double *dp0, *dp, value[2];                             // Pointers used for double data.
495        long seed;
496        long indices[MAX_DIMENSIONS];
497       
498        char buf[256];
499               
500        /* check that wave handles are all valid */
501        if (p->inputWaveH == NIL) {
502                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
503                return(NON_EXISTENT_WAVE);
504        }
505        if (p->ran_devH == NIL) {
506                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
507                return(NON_EXISTENT_WAVE);
508        }       
509        if (p->ntH == NIL) {
510                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
511                return(NON_EXISTENT_WAVE);
512        }
513        if (p->j1H == NIL) {
514                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
515                return(NON_EXISTENT_WAVE);
516        }
517        if (p->j2H == NIL) {
518                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
519                return(NON_EXISTENT_WAVE);
520        }
521        if (p->nnH == NIL) {
522                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
523                return(NON_EXISTENT_WAVE);
524        }
525        if (p->MC_linear_dataH == NIL) {
526                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
527                return(NON_EXISTENT_WAVE);
528        }
529        if (p->resultsH == NIL) {
530                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
531                return(NON_EXISTENT_WAVE);
532        }
533       
534        p->retVal = 0;
535       
536// trusting that all inputs are DOUBLE PRECISION WAVES!!!
537        inputWave = WaveData(p->inputWaveH);
538        ran_dev = WaveData(p->ran_devH);
539        nt = WaveData(p->ntH);
540        j1 = WaveData(p->j1H);
541        j2 = WaveData(p->j2H);
542        nn = WaveData(p->nnH);
543//      MC_linear_data = WaveData(p->MC_linear_dataH);
544        results = WaveData(p->resultsH);
545       
546        seed = (long)results[0];
547       
548//      sprintf(buf, "input seed = %ld\r", seed);
549//      XOPNotice(buf);
550       
551        if(seed >= 0) {
552                seed = -1234509876;
553        }
554
555        dummy = ran3a(&seed);           //initialize the random sequence by passing in a negative value
556        seed = 12348765;                //non-negative after that does nothing
557
558        imon = (int)inputWave[0];
559        r1 = inputWave[1];
560        r2 = inputWave[2];
561        xCtr = inputWave[3];
562        yCtr = inputWave[4];
563        sdd = inputWave[5];
564        pixSize = inputWave[6];
565        thick = inputWave[7];
566        wavelength = inputWave[8];
567        sig_incoh = inputWave[9];
568        sig_sas = inputWave[10];
569        xCtr_long = round(xCtr);
570        yCtr_long = round(yCtr);
571       
572        dummy = MDGetWaveScaling(p->ran_devH, 0, &delta, &left);                //0 is the rows
573        if (retVal = MDGetWaveDimensions(p->ran_devH, &numDimensions, dimensionSizes))
574                return retVal;
575        numRows_ran_dev = dimensionSizes[0];
576       
577        pi = 4.0*atan(1.0);     
578       
579        // access the 2D wave data for writing using the direct method
580        wavH = p->MC_linear_dataH;
581        if (wavH == NIL)
582                return NOWAV;
583
584//      waveType = WaveType(wavH);
585//      if (waveType & NT_CMPLX)
586//              return NO_COMPLEX_WAVE;
587//      if (waveType==TEXT_WAVE_TYPE)
588//              return NUMERIC_ACCESS_ON_TEXT_WAVE;
589//      if (retVal = MDGetWaveDimensions(wavH, &numDimensions, dimensionSizes))
590//              return retVal;
591//      numRows = dimensionSizes[0];
592//      numColumns = dimensionSizes[1];
593       
594//      if (retVal = MDAccessNumericWaveData(wavH, kMDWaveAccessMode0, &dataOffset))
595//              return retVal;
596               
597//      hState = MoveLockHandle(wavH);          // So wave data can't move. Remember to call HSetState when done.
598//      dataStartPtr = (char*)(*wavH) + dataOffset;
599//      dp0 = (double*)dataStartPtr;                    // Pointer to the start of the 2D wave data.
600       
601//scattering power and maximum qvalue to bin
602//      zpow = .1               //scattering power, calculated below
603        qmax = 4.0*pi/wavelength;       //maximum Q to bin 1D data. (A-1) (not really used)
604        sigabs_0 = 0.0;         // ignore absorption cross section/wavelength [1/(cm A)]
605        n_index = 50;   // maximum number of scattering events per neutron
606        num_bins = 200;         //number of 1-D bins (not really used)
607       
608//c       total SAS cross-section
609//
610        zpow = sig_sas*thick;                   //since I now calculate the sig_sas from the model
611        sig_abs = sigabs_0 * wavelength;
612        sig_total = sig_abs + sig_sas + sig_incoh;
613//      Print "The TOTAL XSECTION. (CM-1) is ",sig_total
614//      Print "The TOTAL SAS XSECTION. (CM-1) is ",sig_sas
615//      results[0] = sig_total;
616//      results[1] = sig_sas;
617//      RATIO = SIG_ABS / SIG_TOTAL
618        ratio = sig_incoh / sig_total;
619       
620        theta_max = wavelength*qmax/(2*pi);
621//C     SET Theta-STEP SIZE.
622        dth = theta_max/num_bins;
623//      Print "theta bin size = dth = ",dth
624
625//C     INITIALIZE COUNTERS.
626        n1 = 0;
627        n2 = 0;
628        n3 = 0;
629        NSingleIncoherent = 0;
630        NSingleCoherent = 0;
631        NDoubleCoherent = 0;
632        NMultipleScatter = 0;
633        NScatterEvents = 0;
634        NMultipleCoherent = 0;
635        NCoherentEvents = 0;
636       
637        isOn = 0;
638       
639//C     MONITOR LOOP - looping over the number of incedent neutrons
640//note that zz, is the z-position in the sample - NOT the scattering power
641// NOW, start the loop, throwing neutrons at the sample.
642        do {
643                ////SpinProcess() IS A CALLBACK, and not good for Threading!
644//              if ((n1 % 1000 == 0) && gCallSpinProcess && SpinProcess()) {            // Spins cursor and allows background processing.
645//                              retVal = -1;                                                            // User aborted.
646//                              break;
647//              }
648       
649                vx = 0.0;                       // Initialize direction vector.
650                vy = 0.0;
651                vz = 1.0;
652               
653                theta = 0.0;            //      Initialize scattering angle.
654                phi = 0.0;                      //      Intialize azimuthal angle.
655                n1 += 1;                        //      Increment total number neutrons counter.
656                done = 0;                       //      True when neutron is absorbed or when  scattered out of the sample.
657                index = 0;                      //      Set counter for number of scattering events.
658                zz = 0.0;                       //      Set entering dimension of sample.
659                incoherentEvent = 0;
660                coherentEvent = 0;
661               
662                do      {                               //      Makes sure position is within circle.
663                        ran = ran3a(&seed);             //[0,1]
664                        xx = 2.0*r1*(ran-0.5);          //X beam position of neutron entering sample.
665                        ran = ran3a(&seed);             //[0,1]
666                        yy = 2.0*r1*(ran-0.5);          //Y beam position ...
667                        rr = sqrt(xx*xx+yy*yy);         //Radial position of neutron in incident beam.
668                } while(rr>r1);
669
670                do {   //Scattering Loop, will exit when "done" == 1
671                                // keep scattering multiple times until the neutron exits the sample
672                        ran = ran3a(&seed);             //[0,1]  RANDOM NUMBER FOR DETERMINING PATH LENGTH
673                        ll = path_len(ran,sig_total);
674                        //Determine new scattering direction vector.
675                        err = NewDirection(&vx,&vy,&vz,theta,phi);              //vx,vy,vz updated, theta, phi unchanged by function
676                                                                       
677                        //X,Y,Z-POSITION OF SCATTERING EVENT.
678                        xx += ll*vx;
679                        yy += ll*vy;
680                        zz += ll*vz;
681                        rr = sqrt(xx*xx+yy*yy);         //radial position of scattering event.
682
683                        //sprintf(buf, "xx,yy,zz,vx,vy,vz,ll = %g %g %g %g %g %g %g\r",xx,yy,zz,vx,vy,vz,ll);
684                        //XOPNotice(buf);
685                                               
686                        //Check whether interaction occurred within sample volume.
687                        if (((zz > 0.0) && (zz < thick)) && (rr < r2)) {
688                                //NEUTRON INTERACTED.
689                                //sprintf(buf,"neutron interacted\r");
690                                //XOPNotice(buf);
691                               
692                                index += 1;                     //Increment counter of scattering events.
693                                if (index == 1) {
694                                        n2 += 1;                //Increment # of scat. neutrons
695                                }
696                                ran = ran3a(&seed);             //[0,1]
697                                //Split neutron interactions into scattering and absorption events
698                                if (ran > ratio ) {             //C             NEUTRON SCATTERED coherently
699                                        //sprintf(buf,"neutron scatters coherently\r");
700                                        //XOPNotice(buf);
701                                        coherentEvent += 1;
702                                        find_theta = 0;                 //false
703                                        do {
704                                                // pick a q-value from the deviate function
705                                                // pnt2x truncates the point to an integer before returning the x
706                                                // so get it from the wave scaling instead
707//                                              q0 =left + binarysearchinterp(ran_dev,ran3a(seed))*delta;
708                                               
709                                                q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran3a(&seed))*delta;
710                                                theta = q0/2/pi*wavelength;             //SAS approximation
711                                               
712                                                find_theta = 1;         //always accept
713
714                                                //sprintf(buf, "after locate_interp call q0 = %g, theta = %g,left = %g,delta = %g\r",q0,theta,left,delta);
715                                                //XOPNotice(buf);
716
717                                        } while(!find_theta);
718                                       
719                                        ran = ran3a(&seed);             //[0,1]
720                                        phi = 2.0*pi*ran;                       //Chooses azimuthal scattering angle.
721                                } else {
722                                        //NEUTRON scattered incoherently
723                                        //sprintf(buf,"neutron scatters incoherent\r");
724                                        //XOPNotice(buf);
725                                        incoherentEvent += 1;
726                                  // phi and theta are random over the entire sphere of scattering
727                                        // !can't just choose random theta and phi, won't be random over sphere solid angle
728
729                                        ran = ran3a(&seed);             //[0,1]
730                                        theta = acos(2.0*ran-1);
731                       
732                                        ran = ran3a(&seed);             //[0,1]
733                                        phi = 2.0*pi*ran;                       //Chooses azimuthal scattering angle.
734                                }               //(ran > ratio)
735                        } else {
736                                //NEUTRON ESCAPES FROM SAMPLE -- bin it somewhere                                                               
737                                done = 1;               //done = true, will exit from loop
738                                //Increment #scattering events array
739                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
740                                indices[0] =index;                      //this sets access to nn[index]
741                                if (index <= n_index) {
742                                        if (retVal = MDGetNumericWavePointValue(p->nnH, indices, value))
743                                                return retVal;
744                                        value[0] += 1; // add one to the value
745                                        if (retVal = MDSetNumericWavePointValue(p->nnH, indices, value))
746                                                return retVal;
747                                //      nn[index] += 1;
748                                }
749                                                                                               
750                                if( index != 0) {               //neutron was scattered, figure out where it went
751                                        theta_z = acos(vz);             // Angle (= 2theta) WITH respect to z axis.
752                                        testQ = 2*pi*sin(theta_z)/wavelength;
753                                       
754                                        // pick a random phi angle, and see if it lands on the detector
755                                        // since the scattering is isotropic, I can safely pick a new, random value
756                                        // this would not be true if simulating anisotropic scattering.
757                                        testPhi = ran3a(&seed)*2*pi;
758                                       
759                                        // is it on the detector?       
760                                        FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel);
761                                                                                                       
762                                        if(xPixel != -1 && yPixel != -1) {
763                                                isOn += 1;
764                                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
765                                                indices[0] = xPixel;
766                                                indices[1] = yPixel;
767                                                if (retVal = MDGetNumericWavePointValue(wavH, indices, value))
768                                                        return retVal;
769                                                value[0] += 1; // Real part
770                                                if (retVal = MDSetNumericWavePointValue(wavH, indices, value))
771                                                        return retVal;
772                                                //if(index==1)  // only the single scattering events
773                                                        //dp = dp0 + xPixel + yPixel*numColumns;                //offset the pointer to the exact memory location
774                                                        //*dp += 1;             //increment the value there
775                                                //endif
776                                        }
777                                       
778
779        /*              is this causing me a problem since I'm not locking these? Probably not, since it crashes even if I comment these out... */
780                                        if(theta_z < theta_max) {
781                                                //Choose index for scattering angle array.
782                                                //IND = NINT(THETA_z/DTH + 0.4999999)
783                                                ind = round(theta_z/dth + 0.4999999);           //round is eqivalent to nint()
784                                                nt[ind] += 1;                   //Increment bin for angle.
785                                                //Increment angle array for single scattering events.
786                                                if (index == 1) {
787                                                        j1[ind] += 1;
788                                                }
789                                                //Increment angle array for double scattering events.
790                                                if (index == 2) {
791                                                        j2[ind] += 1;
792                                                }
793                                        }
794        /**/
795                                       
796                                        // increment all of the counters now since done==1 here and I'm sure to exit and get another neutron
797                                        NScatterEvents += index;                //total number of scattering events
798                                        if(index == 1 && incoherentEvent == 1) {
799                                                NSingleIncoherent += 1;
800                                        }
801                                        if(index == 1 && coherentEvent == 1) {
802                                                NSingleCoherent += 1;
803                                        }
804                                        if(index == 2 && coherentEvent == 1 && incoherentEvent == 0) {
805                                                NDoubleCoherent += 1;
806                                        }
807                                        if(index > 1) {
808                                                NMultipleScatter += 1;
809                                        }
810                                        if(coherentEvent >= 1 && incoherentEvent == 0) {
811                                                NCoherentEvents += 1;
812                                        }
813                                        if(coherentEvent > 1 && incoherentEvent == 0) {
814                                                NMultipleCoherent += 1;
815                                        }
816
817                                } else {        // index was zero, neutron must be transmitted, so just increment the proper counters and data
818                                                isOn += 1;
819                                                nt[0] += 1;
820                                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
821                                                //indices[0] = xCtr_long;               //don't put everything in one pixel
822                                                //indices[1] = yCtr_long;
823                                                indices[0] = (long)round(xCtr+xx/pixSize);
824                                                indices[1] = (long)round(yCtr+yy/pixSize);
825                                                // check for valid indices - got an XOP error, probably from here
826                                                if(indices[0] > 127) indices[0] = 127;
827                                                if(indices[0] < 0) indices[0] = 0;
828                                                if(indices[1] > 127) indices[1] = 127;
829                                                if(indices[1] < 0) indices[1] = 0;
830                                               
831                                                if (retVal = MDGetNumericWavePointValue(wavH, indices, value))
832                                                        return retVal;
833                                                value[0] += 1; // Real part
834                                                if (retVal = MDSetNumericWavePointValue(wavH, indices, value))
835                                                        return retVal;
836                                        }       
837                        }
838                 } while (!done);
839        } while(n1 < imon);
840       
841// assign the results to the wave
842
843        MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
844        value[0] = (double)n1;
845        indices[0] = 0;
846        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
847                return retVal;
848        value[0] = (double)n2;
849        indices[0] = 1;
850        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
851                return retVal;
852        value[0] = (double)isOn;
853        indices[0] = 2;
854        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
855                return retVal;
856        value[0] = (double)NScatterEvents;
857        indices[0] = 3;
858        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
859                return retVal;
860        value[0] = (double)NSingleCoherent;
861        indices[0] = 4;
862        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
863                return retVal;
864        value[0] = (double)NMultipleCoherent;
865        indices[0] = 5;
866        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
867                return retVal;
868        value[0] = (double)NMultipleScatter;
869        indices[0] = 6;
870        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
871                return retVal; 
872        value[0] = (double)NCoherentEvents;
873        indices[0] = 7;
874        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
875                return retVal;
876
877//      HSetState((Handle)wavH, hState);                //release the handle of the 2D data wave
878//      WaveHandleModified(wavH);                       // Inform Igor that we have changed the wave. (CALLBACK! needed, but not allowed in Threading)
879       
880        return(0);
881}
882////////        end of X3
883
884
885///////////////// X4 using ran1a()
886int
887Monte_SANSX4(MC_ParamsPtr p) {
888        double *inputWave;                              /* pointer to double precision wave data */
889        double *ran_dev;                                /* pointer to double precision wave data */
890        double *nt;                             /* pointer to double precision wave data */
891        double *j1;                             /* pointer to double precision wave data */
892        double *j2;                             /* pointer to double precision wave data */
893        double *nn;                             /* pointer to double precision wave data */
894//      double *MC_linear_data;                         /* pointer to double precision wave data */
895        double *results;                                /* pointer to double precision wave data */
896        double retVal;                          //return value
897
898        long imon;
899        double r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas;
900        long ind,index,n_index;
901        double qmax,theta_max,q0,zpow;
902        long n1,n2,n3;
903        double dth,zz,xx,yy,phi;
904        double theta,ran,ll,rr,ttot;
905        long done,find_theta,err;               //used as logicals
906        long xPixel,yPixel;
907        double vx,vy,vz,theta_z;
908        double sig_abs,ratio,sig_total;
909        double testQ,testPhi,left,delta,dummy,pi;
910        double sigabs_0,num_bins;
911        long NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent;
912        long NDoubleCoherent,NMultipleScatter,isOn,xCtr_long,yCtr_long;
913        long NMultipleCoherent,NCoherentEvents;
914
915
916        // for accessing the 2D wave data, direct method (see the WaveAccess example XOP)
917        waveHndl wavH;
918        int waveType,hState;
919        long numDimensions;
920        long dimensionSizes[MAX_DIMENSIONS+1];
921        char* dataStartPtr;
922        long dataOffset;
923        long numRows, numColumns,numRows_ran_dev;
924        double *dp0, *dp, value[2];                             // Pointers used for double data.
925        long seed;
926        long indices[MAX_DIMENSIONS];
927       
928        char buf[256];
929               
930        /* check that wave handles are all valid */
931        if (p->inputWaveH == NIL) {
932                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
933                return(NON_EXISTENT_WAVE);
934        }
935        if (p->ran_devH == NIL) {
936                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
937                return(NON_EXISTENT_WAVE);
938        }       
939        if (p->ntH == NIL) {
940                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
941                return(NON_EXISTENT_WAVE);
942        }
943        if (p->j1H == NIL) {
944                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
945                return(NON_EXISTENT_WAVE);
946        }
947        if (p->j2H == NIL) {
948                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
949                return(NON_EXISTENT_WAVE);
950        }
951        if (p->nnH == NIL) {
952                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
953                return(NON_EXISTENT_WAVE);
954        }
955        if (p->MC_linear_dataH == NIL) {
956                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
957                return(NON_EXISTENT_WAVE);
958        }
959        if (p->resultsH == NIL) {
960                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
961                return(NON_EXISTENT_WAVE);
962        }
963       
964        p->retVal = 0;
965       
966// trusting that all inputs are DOUBLE PRECISION WAVES!!!
967        inputWave = WaveData(p->inputWaveH);
968        ran_dev = WaveData(p->ran_devH);
969        nt = WaveData(p->ntH);
970        j1 = WaveData(p->j1H);
971        j2 = WaveData(p->j2H);
972        nn = WaveData(p->nnH);
973//      MC_linear_data = WaveData(p->MC_linear_dataH);
974        results = WaveData(p->resultsH);
975       
976        seed = (long)results[0];
977       
978//      sprintf(buf, "input seed = %ld\r", seed);
979//      XOPNotice(buf);
980       
981        if(seed >= 0) {
982                seed = -1234509876;
983        }
984
985        dummy = ran1a(&seed);           //initialize the random sequence by passing in a negative value
986        seed = 12348765;                //non-negative after that does nothing
987
988        imon = (int)inputWave[0];
989        r1 = inputWave[1];
990        r2 = inputWave[2];
991        xCtr = inputWave[3];
992        yCtr = inputWave[4];
993        sdd = inputWave[5];
994        pixSize = inputWave[6];
995        thick = inputWave[7];
996        wavelength = inputWave[8];
997        sig_incoh = inputWave[9];
998        sig_sas = inputWave[10];
999        xCtr_long = round(xCtr);
1000        yCtr_long = round(yCtr);
1001       
1002        dummy = MDGetWaveScaling(p->ran_devH, 0, &delta, &left);                //0 is the rows
1003        if (retVal = MDGetWaveDimensions(p->ran_devH, &numDimensions, dimensionSizes))
1004                return retVal;
1005        numRows_ran_dev = dimensionSizes[0];
1006       
1007        pi = 4.0*atan(1.0);     
1008       
1009        // access the 2D wave data for writing using the direct method
1010        wavH = p->MC_linear_dataH;
1011        if (wavH == NIL)
1012                return NOWAV;
1013
1014//      waveType = WaveType(wavH);
1015//      if (waveType & NT_CMPLX)
1016//              return NO_COMPLEX_WAVE;
1017//      if (waveType==TEXT_WAVE_TYPE)
1018//              return NUMERIC_ACCESS_ON_TEXT_WAVE;
1019//      if (retVal = MDGetWaveDimensions(wavH, &numDimensions, dimensionSizes))
1020//              return retVal;
1021//      numRows = dimensionSizes[0];
1022//      numColumns = dimensionSizes[1];
1023       
1024//      if (retVal = MDAccessNumericWaveData(wavH, kMDWaveAccessMode0, &dataOffset))
1025//              return retVal;
1026               
1027//      hState = MoveLockHandle(wavH);          // So wave data can't move. Remember to call HSetState when done.
1028//      dataStartPtr = (char*)(*wavH) + dataOffset;
1029//      dp0 = (double*)dataStartPtr;                    // Pointer to the start of the 2D wave data.
1030       
1031//scattering power and maximum qvalue to bin
1032//      zpow = .1               //scattering power, calculated below
1033        qmax = 4.0*pi/wavelength;       //maximum Q to bin 1D data. (A-1) (not really used)
1034        sigabs_0 = 0.0;         // ignore absorption cross section/wavelength [1/(cm A)]
1035        n_index = 50;   // maximum number of scattering events per neutron
1036        num_bins = 200;         //number of 1-D bins (not really used)
1037       
1038//c       total SAS cross-section
1039//
1040        zpow = sig_sas*thick;                   //since I now calculate the sig_sas from the model
1041        sig_abs = sigabs_0 * wavelength;
1042        sig_total = sig_abs + sig_sas + sig_incoh;
1043//      Print "The TOTAL XSECTION. (CM-1) is ",sig_total
1044//      Print "The TOTAL SAS XSECTION. (CM-1) is ",sig_sas
1045//      results[0] = sig_total;
1046//      results[1] = sig_sas;
1047//      RATIO = SIG_ABS / SIG_TOTAL
1048        ratio = sig_incoh / sig_total;
1049       
1050        theta_max = wavelength*qmax/(2*pi);
1051//C     SET Theta-STEP SIZE.
1052        dth = theta_max/num_bins;
1053//      Print "theta bin size = dth = ",dth
1054
1055//C     INITIALIZE COUNTERS.
1056        n1 = 0;
1057        n2 = 0;
1058        n3 = 0;
1059        NSingleIncoherent = 0;
1060        NSingleCoherent = 0;
1061        NDoubleCoherent = 0;
1062        NMultipleScatter = 0;
1063        NScatterEvents = 0;
1064        NMultipleCoherent = 0;
1065        NCoherentEvents = 0;
1066       
1067        isOn = 0;
1068       
1069//C     MONITOR LOOP - looping over the number of incedent neutrons
1070//note that zz, is the z-position in the sample - NOT the scattering power
1071// NOW, start the loop, throwing neutrons at the sample.
1072        do {
1073                ////SpinProcess() IS A CALLBACK, and not good for Threading!
1074//              if ((n1 % 1000 == 0) && gCallSpinProcess && SpinProcess()) {            // Spins cursor and allows background processing.
1075//                              retVal = -1;                                                            // User aborted.
1076//                              break;
1077//              }
1078       
1079                vx = 0.0;                       // Initialize direction vector.
1080                vy = 0.0;
1081                vz = 1.0;
1082               
1083                theta = 0.0;            //      Initialize scattering angle.
1084                phi = 0.0;                      //      Intialize azimuthal angle.
1085                n1 += 1;                        //      Increment total number neutrons counter.
1086                done = 0;                       //      True when neutron is absorbed or when  scattered out of the sample.
1087                index = 0;                      //      Set counter for number of scattering events.
1088                zz = 0.0;                       //      Set entering dimension of sample.
1089                incoherentEvent = 0;
1090                coherentEvent = 0;
1091               
1092                do      {                               //      Makes sure position is within circle.
1093                        ran = ran1a(&seed);             //[0,1]
1094                        xx = 2.0*r1*(ran-0.5);          //X beam position of neutron entering sample.
1095                        ran = ran1a(&seed);             //[0,1]
1096                        yy = 2.0*r1*(ran-0.5);          //Y beam position ...
1097                        rr = sqrt(xx*xx+yy*yy);         //Radial position of neutron in incident beam.
1098                } while(rr>r1);
1099
1100                do {   //Scattering Loop, will exit when "done" == 1
1101                                // keep scattering multiple times until the neutron exits the sample
1102                        ran = ran1a(&seed);             //[0,1]  RANDOM NUMBER FOR DETERMINING PATH LENGTH
1103                        ll = path_len(ran,sig_total);
1104                        //Determine new scattering direction vector.
1105                        err = NewDirection(&vx,&vy,&vz,theta,phi);              //vx,vy,vz updated, theta, phi unchanged by function
1106                                                                       
1107                        //X,Y,Z-POSITION OF SCATTERING EVENT.
1108                        xx += ll*vx;
1109                        yy += ll*vy;
1110                        zz += ll*vz;
1111                        rr = sqrt(xx*xx+yy*yy);         //radial position of scattering event.
1112
1113                        //sprintf(buf, "xx,yy,zz,vx,vy,vz,ll = %g %g %g %g %g %g %g\r",xx,yy,zz,vx,vy,vz,ll);
1114                        //XOPNotice(buf);
1115                                               
1116                        //Check whether interaction occurred within sample volume.
1117                        if (((zz > 0.0) && (zz < thick)) && (rr < r2)) {
1118                                //NEUTRON INTERACTED.
1119                                //sprintf(buf,"neutron interacted\r");
1120                                //XOPNotice(buf);
1121                               
1122                                index += 1;                     //Increment counter of scattering events.
1123                                if (index == 1) {
1124                                        n2 += 1;                //Increment # of scat. neutrons
1125                                }
1126                                ran = ran1a(&seed);             //[0,1]
1127                                //Split neutron interactions into scattering and absorption events
1128                                if (ran > ratio ) {             //C             NEUTRON SCATTERED coherently
1129                                        //sprintf(buf,"neutron scatters coherently\r");
1130                                        //XOPNotice(buf);
1131                                        coherentEvent += 1;
1132                                        find_theta = 0;                 //false
1133                                        do {
1134                                                // pick a q-value from the deviate function
1135                                                // pnt2x truncates the point to an integer before returning the x
1136                                                // so get it from the wave scaling instead
1137//                                              q0 =left + binarysearchinterp(ran_dev,ran1a(seed))*delta;
1138                                               
1139                                                q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran1a(&seed))*delta;
1140                                                theta = q0/2/pi*wavelength;             //SAS approximation
1141                                               
1142                                                find_theta = 1;         //always accept
1143
1144                                                //sprintf(buf, "after locate_interp call q0 = %g, theta = %g,left = %g,delta = %g\r",q0,theta,left,delta);
1145                                                //XOPNotice(buf);
1146
1147                                        } while(!find_theta);
1148                                       
1149                                        ran = ran1a(&seed);             //[0,1]
1150                                        phi = 2.0*pi*ran;                       //Chooses azimuthal scattering angle.
1151                                } else {
1152                                        //NEUTRON scattered incoherently
1153                                        //sprintf(buf,"neutron scatters incoherent\r");
1154                                        //XOPNotice(buf);
1155                                        incoherentEvent += 1;
1156                                  // phi and theta are random over the entire sphere of scattering
1157                                        // !can't just choose random theta and phi, won't be random over sphere solid angle
1158
1159                                        ran = ran1a(&seed);             //[0,1]
1160                                        theta = acos(2.0*ran-1);
1161                       
1162                                        ran = ran1a(&seed);             //[0,1]
1163                                        phi = 2.0*pi*ran;                       //Chooses azimuthal scattering angle.
1164                                }               //(ran > ratio)
1165                        } else {
1166                                //NEUTRON ESCAPES FROM SAMPLE -- bin it somewhere                                                               
1167                                done = 1;               //done = true, will exit from loop
1168                                //Increment #scattering events array
1169                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
1170                                indices[0] =index;                      //this sets access to nn[index]
1171                                if (index <= n_index) {
1172                                        if (retVal = MDGetNumericWavePointValue(p->nnH, indices, value))
1173                                                return retVal;
1174                                        value[0] += 1; // add one to the value
1175                                        if (retVal = MDSetNumericWavePointValue(p->nnH, indices, value))
1176                                                return retVal;
1177                                //      nn[index] += 1;
1178                                }
1179                                                                                               
1180                                if( index != 0) {               //neutron was scattered, figure out where it went
1181                                        theta_z = acos(vz);             // Angle (= 2theta) WITH respect to z axis.
1182                                        testQ = 2*pi*sin(theta_z)/wavelength;
1183                                       
1184                                        // pick a random phi angle, and see if it lands on the detector
1185                                        // since the scattering is isotropic, I can safely pick a new, random value
1186                                        // this would not be true if simulating anisotropic scattering.
1187                                        testPhi = ran1a(&seed)*2*pi;
1188                                       
1189                                        // is it on the detector?       
1190                                        FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel);
1191                                                                                                       
1192                                        if(xPixel != -1 && yPixel != -1) {
1193                                                isOn += 1;
1194                                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
1195                                                indices[0] = xPixel;
1196                                                indices[1] = yPixel;
1197                                                if (retVal = MDGetNumericWavePointValue(wavH, indices, value))
1198                                                        return retVal;
1199                                                value[0] += 1; // Real part
1200                                                if (retVal = MDSetNumericWavePointValue(wavH, indices, value))
1201                                                        return retVal;
1202                                                //if(index==1)  // only the single scattering events
1203                                                        //dp = dp0 + xPixel + yPixel*numColumns;                //offset the pointer to the exact memory location
1204                                                        //*dp += 1;             //increment the value there
1205                                                //endif
1206                                        }
1207                                       
1208
1209        /*              is this causing me a problem since I'm not locking these? Probably not, since it crashes even if I comment these out... */
1210                                        if(theta_z < theta_max) {
1211                                                //Choose index for scattering angle array.
1212                                                //IND = NINT(THETA_z/DTH + 0.4999999)
1213                                                ind = round(theta_z/dth + 0.4999999);           //round is eqivalent to nint()
1214                                                nt[ind] += 1;                   //Increment bin for angle.
1215                                                //Increment angle array for single scattering events.
1216                                                if (index == 1) {
1217                                                        j1[ind] += 1;
1218                                                }
1219                                                //Increment angle array for double scattering events.
1220                                                if (index == 2) {
1221                                                        j2[ind] += 1;
1222                                                }
1223                                        }
1224        /**/
1225                                       
1226                                        // increment all of the counters now since done==1 here and I'm sure to exit and get another neutron
1227                                        NScatterEvents += index;                //total number of scattering events
1228                                        if(index == 1 && incoherentEvent == 1) {
1229                                                NSingleIncoherent += 1;
1230                                        }
1231                                        if(index == 1 && coherentEvent == 1) {
1232                                                NSingleCoherent += 1;
1233                                        }
1234                                        if(index == 2 && coherentEvent == 1 && incoherentEvent == 0) {
1235                                                NDoubleCoherent += 1;
1236                                        }
1237                                        if(index > 1) {
1238                                                NMultipleScatter += 1;
1239                                        }
1240                                        if(coherentEvent >= 1 && incoherentEvent == 0) {
1241                                                NCoherentEvents += 1;
1242                                        }
1243                                        if(coherentEvent > 1 && incoherentEvent == 0) {
1244                                                NMultipleCoherent += 1;
1245                                        }
1246
1247                                } else {        // index was zero, neutron must be transmitted, so just increment the proper counters and data
1248                                                isOn += 1;
1249                                                nt[0] += 1;
1250                                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
1251                                                //indices[0] = xCtr_long;               //don't put everything in one pixel
1252                                                //indices[1] = yCtr_long;
1253                                                indices[0] = (long)round(xCtr+xx/pixSize);
1254                                                indices[1] = (long)round(yCtr+yy/pixSize);
1255                                                // check for valid indices - got an XOP error, probably from here
1256                                                if(indices[0] > 127) indices[0] = 127;
1257                                                if(indices[0] < 0) indices[0] = 0;
1258                                                if(indices[1] > 127) indices[1] = 127;
1259                                                if(indices[1] < 0) indices[1] = 0;
1260                                               
1261                                                if (retVal = MDGetNumericWavePointValue(wavH, indices, value))
1262                                                        return retVal;
1263                                                value[0] += 1; // Real part
1264                                                if (retVal = MDSetNumericWavePointValue(wavH, indices, value))
1265                                                        return retVal;
1266                                        }       
1267                        }
1268                 } while (!done);
1269        } while(n1 < imon);
1270       
1271// assign the results to the wave
1272
1273        MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
1274        value[0] = (double)n1;
1275        indices[0] = 0;
1276        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
1277                return retVal;
1278        value[0] = (double)n2;
1279        indices[0] = 1;
1280        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
1281                return retVal;
1282        value[0] = (double)isOn;
1283        indices[0] = 2;
1284        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
1285                return retVal;
1286        value[0] = (double)NScatterEvents;
1287        indices[0] = 3;
1288        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
1289                return retVal;
1290        value[0] = (double)NSingleCoherent;
1291        indices[0] = 4;
1292        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
1293                return retVal;
1294        value[0] = (double)NMultipleCoherent;
1295        indices[0] = 5;
1296        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
1297                return retVal;
1298        value[0] = (double)NMultipleScatter;
1299        indices[0] = 6;
1300        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
1301                return retVal; 
1302        value[0] = (double)NCoherentEvents;
1303        indices[0] = 7;
1304        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
1305                return retVal;
1306
1307//      HSetState((Handle)wavH, hState);                //release the handle of the 2D data wave
1308//      WaveHandleModified(wavH);                       // Inform Igor that we have changed the wave. (CALLBACK! needed, but not allowed in Threading)
1309       
1310        return(0);
1311}
1312////////        end of X4
1313
1314
1315
Note: See TracBrowser for help on using the repository browser.