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

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

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