source: sans/XOP_Dev/MonteCarlo/MonteCarlo.c @ 591

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

Cosmetic changes to monteCarlo code. Added some ".0" here and there to make sure that all of the numerical values are really interpreted as floating point. Also changed the variable "result" to be "retVal" to avoid confusion with the wave of returned values "results"

File size: 21.2 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//////////
17//    PROGRAM Monte_SANS
18//    PROGRAM simulates multiple SANS.
19//       revised 2/12/99  JGB
20//            added calculation of random deviate, and 2D 10/2008 SRK
21
22//    N1 = NUMBER OF INCIDENT NEUTRONS.
23//    N2 = NUMBER INTERACTED IN THE SAMPLE.
24//    N3 = NUMBER ABSORBED.
25//    THETA = SCATTERING ANGLE.
26
27// works fine in the single-threaded case.
28//
29/// currently crashes if threaded. apparently something here is either doing an unknown callback, or is accessing
30// a bad place in memory.
31//
32// the operations SpinProcess() and WaveHandleModified() are callbacks and MUST be commented out before threading.
33// - some locks are non-existent
34// - supposedly safe wave access routines are used
35//
36// random number generators are not thread-safe, and can give less than random results, but is this enough to crash?
37// -- a possible workaround is to define multiple versions (poor man's threading)
38//
39//
40//
41
42
43int
44Monte_SANSX(MC_ParamsPtr p) {
45        double *inputWave;                              /* pointer to double precision wave data */
46        double *ran_dev;                                /* pointer to double precision wave data */
47        double *nt;                             /* pointer to double precision wave data */
48        double *j1;                             /* pointer to double precision wave data */
49        double *j2;                             /* pointer to double precision wave data */
50        double *nn;                             /* pointer to double precision wave data */
51//      double *MC_linear_data;                         /* pointer to double precision wave data */
52        double *results;                                /* pointer to double precision wave data */
53        double retVal;                          //return value
54
55        long imon;
56        double r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas;
57        long ind,index,n_index;
58        double qmax,theta_max,q0,zpow;
59        long n1,n2,n3;
60        double dth,zz,xx,yy,phi;
61        double theta,ran,ll,rr,ttot;
62        long done,find_theta,err;               //used as logicals
63        long xPixel,yPixel;
64        double vx,vy,vz,theta_z;
65        double sig_abs,ratio,sig_total;
66        double testQ,testPhi,left,delta,dummy,pi;
67        double sigabs_0,num_bins;
68        long NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent;
69        long NDoubleCoherent,NMultipleScatter,isOn,xCtr_long,yCtr_long;
70        long NMultipleCoherent,NCoherentEvents;
71
72        // for accessing the 2D wave data, direct method (see the WaveAccess example XOP)
73        waveHndl wavH;
74        int waveType,hState;
75        long numDimensions;
76        long dimensionSizes[MAX_DIMENSIONS+1];
77        char* dataStartPtr;
78        long dataOffset;
79        long numRows, numColumns,numRows_ran_dev;
80        double *dp0, *dp, value[2];                             // Pointers used for double data.
81        long seed;
82        long indices[MAX_DIMENSIONS];
83       
84        char buf[256];
85               
86        /* check that wave handles are all valid */
87        if (p->inputWaveH == NIL) {
88                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
89                return(NON_EXISTENT_WAVE);
90        }
91        if (p->ran_devH == NIL) {
92                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
93                return(NON_EXISTENT_WAVE);
94        }       
95        if (p->ntH == NIL) {
96                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
97                return(NON_EXISTENT_WAVE);
98        }
99        if (p->j1H == NIL) {
100                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
101                return(NON_EXISTENT_WAVE);
102        }
103        if (p->j2H == NIL) {
104                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
105                return(NON_EXISTENT_WAVE);
106        }
107        if (p->nnH == NIL) {
108                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
109                return(NON_EXISTENT_WAVE);
110        }
111        if (p->MC_linear_dataH == NIL) {
112                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
113                return(NON_EXISTENT_WAVE);
114        }
115        if (p->resultsH == NIL) {
116                SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */
117                return(NON_EXISTENT_WAVE);
118        }
119       
120        p->retVal = 0.0;
121       
122// trusting that all inputs are DOUBLE PRECISION WAVES!!!
123        inputWave = WaveData(p->inputWaveH);
124        ran_dev = WaveData(p->ran_devH);
125        nt = WaveData(p->ntH);
126        j1 = WaveData(p->j1H);
127        j2 = WaveData(p->j2H);
128        nn = WaveData(p->nnH);
129//      MC_linear_data = WaveData(p->MC_linear_dataH);
130        results = WaveData(p->resultsH);
131       
132        seed = (long)results[0];
133       
134//      sprintf(buf, "input seed = %ld\r", seed);
135//      XOPNotice(buf);
136       
137        if(seed >= 0) {
138                seed = -1234509876;
139        }
140
141        dummy = ran3(&seed);            //initialize the random sequence by passing in a negative value
142        seed = 12348765;                //non-negative after that does nothing
143
144        imon = (int)inputWave[0];
145        r1 = inputWave[1];
146        r2 = inputWave[2];
147        xCtr = inputWave[3];
148        yCtr = inputWave[4];
149        sdd = inputWave[5];
150        pixSize = inputWave[6];
151        thick = inputWave[7];
152        wavelength = inputWave[8];
153        sig_incoh = inputWave[9];
154        sig_sas = inputWave[10];
155        xCtr_long = round(xCtr);
156        yCtr_long = round(yCtr);
157       
158        dummy = MDGetWaveScaling(p->ran_devH, 0, &delta, &left);                //0 is the rows
159        if (retVal = MDGetWaveDimensions(p->ran_devH, &numDimensions, dimensionSizes))
160                return retVal;
161        numRows_ran_dev = dimensionSizes[0];
162       
163        pi = 4.0*atan(1.0);     
164       
165        // access the 2D wave data for writing using the direct method
166        wavH = p->MC_linear_dataH;
167        if (wavH == NIL)
168                return NOWAV;
169
170//      waveType = WaveType(wavH);
171//      if (waveType & NT_CMPLX)
172//              return NO_COMPLEX_WAVE;
173//      if (waveType==TEXT_WAVE_TYPE)
174//              return NUMERIC_ACCESS_ON_TEXT_WAVE;
175//      if (retVal = MDGetWaveDimensions(wavH, &numDimensions, dimensionSizes))
176//              return retVal;
177//      numRows = dimensionSizes[0];
178//      numColumns = dimensionSizes[1];
179       
180//      if (retVal = MDAccessNumericWaveData(wavH, kMDWaveAccessMode0, &dataOffset))
181//              return retVal;
182               
183//      hState = MoveLockHandle(wavH);          // So wave data can't move. Remember to call HSetState when done.
184//      dataStartPtr = (char*)(*wavH) + dataOffset;
185//      dp0 = (double*)dataStartPtr;                    // Pointer to the start of the 2D wave data.
186       
187//scattering power and maximum qvalue to bin
188//      zpow = .1               //scattering power, calculated below
189        qmax = 4.0*pi/wavelength;       //maximum Q to bin 1D data. (A-1) (not really used)
190        sigabs_0 = 0.0;         // ignore absorption cross section/wavelength [1/(cm A)]
191        n_index = 50;   // maximum number of scattering events per neutron
192        num_bins = 200;         //number of 1-D bins (not really used)
193       
194//c       total SAS cross-section
195//
196        zpow = sig_sas*thick;                   //since I now calculate the sig_sas from the model
197        sig_abs = sigabs_0 * wavelength;
198        sig_total = sig_abs + sig_sas + sig_incoh;
199//      Print "The TOTAL XSECTION. (CM-1) is ",sig_total
200//      Print "The TOTAL SAS XSECTION. (CM-1) is ",sig_sas
201//      results[0] = sig_total;
202//      results[1] = sig_sas;
203//      RATIO = SIG_ABS / SIG_TOTAL
204        ratio = sig_incoh / sig_total;
205       
206        theta_max = wavelength*qmax/(2.0*pi);
207//C     SET Theta-STEP SIZE.
208        dth = theta_max/num_bins;
209//      Print "theta bin size = dth = ",dth
210
211//C     INITIALIZE COUNTERS.
212        n1 = 0;
213        n2 = 0;
214        n3 = 0;
215        NSingleIncoherent = 0;
216        NSingleCoherent = 0;
217        NDoubleCoherent = 0;
218        NMultipleScatter = 0;
219        NScatterEvents = 0;
220        NMultipleCoherent = 0;
221        NCoherentEvents = 0;
222
223        isOn = 0;
224       
225//C     MONITOR LOOP - looping over the number of incedent neutrons
226//note that zz, is the z-position in the sample - NOT the scattering power
227// NOW, start the loop, throwing neutrons at the sample.
228        do {
229                ////SpinProcess() IS A CALLBACK, and not good for Threading!
230//              if ((n1 % 1000 == 0) && gCallSpinProcess && SpinProcess()) {            // Spins cursor and allows background processing.
231//                              retVal = -1;                                                            // User aborted.
232//                              break;
233//              }
234       
235                vx = 0.0;                       // Initialize direction vector.
236                vy = 0.0;
237                vz = 1.0;
238               
239                theta = 0.0;            //      Initialize scattering angle.
240                phi = 0.0;                      //      Intialize azimuthal angle.
241                n1 += 1;                        //      Increment total number neutrons counter.
242                done = 0;                       //      True when neutron is absorbed or when  scattered out of the sample.
243                index = 0;                      //      Set counter for number of scattering events.
244                zz = 0.0;                       //      Set entering dimension of sample.
245                incoherentEvent = 0;
246                coherentEvent = 0;
247               
248                do      {                               //      Makes sure position is within circle.
249                        ran = ran3(&seed);              //[0,1]
250                        xx = 2.0*r1*(ran-0.5);          //X beam position of neutron entering sample.
251                        ran = ran3(&seed);              //[0,1]
252                        yy = 2.0*r1*(ran-0.5);          //Y beam position ...
253                        rr = sqrt(xx*xx+yy*yy);         //Radial position of neutron in incident beam.
254                } while(rr>r1);
255
256                do {   //Scattering Loop, will exit when "done" == 1
257                                // keep scattering multiple times until the neutron exits the sample
258                        ran = ran3(&seed);              //[0,1]  RANDOM NUMBER FOR DETERMINING PATH LENGTH
259                        ll = path_len(ran,sig_total);
260                        //Determine new scattering direction vector.
261                        err = NewDirection(&vx,&vy,&vz,theta,phi);              //vx,vy,vz updated, theta, phi unchanged by function
262                                                                       
263                        //X,Y,Z-POSITION OF SCATTERING EVENT.
264                        xx += ll*vx;
265                        yy += ll*vy;
266                        zz += ll*vz;
267                        rr = sqrt(xx*xx+yy*yy);         //radial position of scattering event.
268
269                        //sprintf(buf, "xx,yy,zz,vx,vy,vz,ll = %g %g %g %g %g %g %g\r",xx,yy,zz,vx,vy,vz,ll);
270                        //XOPNotice(buf);
271                                               
272                        //Check whether interaction occurred within sample volume.
273                        if (((zz > 0.0) && (zz < thick)) && (rr < r2)) {
274                                //NEUTRON INTERACTED.
275                                //sprintf(buf,"neutron interacted\r");
276                                //XOPNotice(buf);
277                               
278                                index += 1;                     //Increment counter of scattering events.
279                                if (index == 1) {
280                                        n2 += 1;                //Increment # of scat. neutrons
281                                }
282                                ran = ran3(&seed);              //[0,1]
283                                //Split neutron interactions into scattering and absorption events
284                                if (ran > ratio ) {             //C             NEUTRON SCATTERED coherently
285                                        //sprintf(buf,"neutron scatters coherently\r");
286                                        //XOPNotice(buf);
287                                        coherentEvent += 1;
288                                        find_theta = 0;                 //false
289                                        do {
290                                                // pick a q-value from the deviate function
291                                                // pnt2x truncates the point to an integer before returning the x
292                                                // so get it from the wave scaling instead
293//                                              q0 =left + binarysearchinterp(ran_dev,ran1(seed))*delta;
294                                               
295                                                q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran3(&seed))*delta;
296                                                theta = q0/2.0/pi*wavelength;           //SAS approximation. 1% error at theta=30 degrees
297                                               
298                                                find_theta = 1;         //always accept
299
300                                                //sprintf(buf, "after locate_interp call q0 = %g, theta = %g,left = %g,delta = %g\r",q0,theta,left,delta);
301                                                //XOPNotice(buf);
302
303                                        } while(!find_theta);
304                                       
305                                        ran = ran3(&seed);              //[0,1]
306                                        phi = 2.0*pi*ran;                       //Chooses azimuthal scattering angle.
307                                } else {
308                                        //NEUTRON scattered incoherently
309                                        //sprintf(buf,"neutron scatters incoherent\r");
310                                        //XOPNotice(buf);
311                                        incoherentEvent += 1;
312                                  // phi and theta are random over the entire sphere of scattering
313                                        // !can't just choose random theta and phi, won't be random over sphere solid angle
314
315                                        ran = ran3(&seed);              //[0,1]
316                                        theta = acos(2.0*ran-1);
317                       
318                                        ran = ran3(&seed);              //[0,1]
319                                        phi = 2.0*pi*ran;                       //Chooses azimuthal scattering angle.
320                                }               //(ran > ratio)
321                        } else {
322                                //NEUTRON ESCAPES FROM SAMPLE -- bin it somewhere                                                               
323                                done = 1;               //done = true, will exit from loop
324                                //Increment #scattering events array
325                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
326                                indices[0] =index;                      //this sets access to nn[index]
327                                if (index <= n_index) {
328                                        if (retVal = MDGetNumericWavePointValue(p->nnH, indices, value))
329                                                return retVal;
330                                        value[0] += 1; // add one to the value
331                                        if (retVal = MDSetNumericWavePointValue(p->nnH, indices, value))
332                                                return retVal;
333                                //      nn[index] += 1;
334                                }
335                                                                                               
336                                if( index != 0) {               //neutron was scattered, figure out where it went
337                                        theta_z = acos(vz);             // Angle (= 2theta) WITH respect to z axis.
338                                        testQ = 2.0*pi*sin(theta_z)/wavelength;
339                                       
340                                        // pick a random phi angle, and see if it lands on the detector
341                                        // since the scattering is isotropic, I can safely pick a new, random value
342                                        // this would not be true if simulating anisotropic scattering.
343                                        testPhi = ran3(&seed)*2.0*pi;
344                                       
345                                        // is it on the detector?       
346                                        FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel);
347                                                                                                       
348                                        if(xPixel != -1 && yPixel != -1) {
349                                                isOn += 1;
350                                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
351                                                indices[0] = xPixel;
352                                                indices[1] = yPixel;
353                                                if (retVal = MDGetNumericWavePointValue(wavH, indices, value))
354                                                        return retVal;
355                                                value[0] += 1; // Real part
356                                                if (retVal = MDSetNumericWavePointValue(wavH, indices, value))
357                                                        return retVal;
358                                                //if(index==1)  // only the single scattering events
359                                                        //dp = dp0 + xPixel + yPixel*numColumns;                //offset the pointer to the exact memory location
360                                                        //*dp += 1;             //increment the value there
361                                                //endif
362                                        }
363                                       
364
365        /*              is this causing me a problem since I'm not locking these? Probably not, since it crashes even if I comment these out... */
366                                        if(theta_z < theta_max) {
367                                                //Choose index for scattering angle array.
368                                                //IND = NINT(THETA_z/DTH + 0.4999999)
369                                                ind = round(theta_z/dth + 0.4999999);           //round is eqivalent to nint()
370                                                nt[ind] += 1;                   //Increment bin for angle.
371                                                //Increment angle array for single scattering events.
372                                                if (index == 1) {
373                                                        j1[ind] += 1;
374                                                }
375                                                //Increment angle array for double scattering events.
376                                                if (index == 2) {
377                                                        j2[ind] += 1;
378                                                }
379                                        }
380        /**/
381                                       
382                                        // increment all of the counters now since done==1 here and I'm sure to exit and get another neutron
383                                        NScatterEvents += index;                //total number of scattering events
384                                        if(index == 1 && incoherentEvent == 1) {
385                                                NSingleIncoherent += 1;
386                                        }
387                                        if(index == 1 && coherentEvent == 1) {
388                                                NSingleCoherent += 1;
389                                        }
390                                        if(index == 2 && coherentEvent == 1 && incoherentEvent == 0) {
391                                                NDoubleCoherent += 1;
392                                        }
393                                        if(index > 1) {
394                                                NMultipleScatter += 1;
395                                        }
396                                        if(coherentEvent >= 1 && incoherentEvent == 0) {
397                                                NCoherentEvents += 1;
398                                        }
399                                        if(coherentEvent > 1 && incoherentEvent == 0) {
400                                                NMultipleCoherent += 1;
401                                        }
402
403                                } else {        // index was zero, neutron must be transmitted, so just increment the proper counters and data
404                                                isOn += 1;
405                                                nt[0] += 1;
406                                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
407                                                //indices[0] = xCtr_long;               //don't put everything in one pixel
408                                                //indices[1] = yCtr_long;
409                                                indices[0] = (long)round(xCtr+xx/pixSize);
410                                                indices[1] = (long)round(yCtr+yy/pixSize);
411                                                // check for valid indices - got an XOP error, probably from here
412                                                if(indices[0] > 127) indices[0] = 127;
413                                                if(indices[0] < 0) indices[0] = 0;
414                                                if(indices[1] > 127) indices[1] = 127;
415                                                if(indices[1] < 0) indices[1] = 0;
416                                               
417                                                if (retVal = MDGetNumericWavePointValue(wavH, indices, value))
418                                                        return retVal;
419                                                value[0] += 1; // Real part
420                                                if (retVal = MDSetNumericWavePointValue(wavH, indices, value))
421                                                        return retVal;
422                                        }       
423                        }
424                 } while (!done);
425        } while(n1 < imon);
426       
427// assign the results to the wave
428
429        MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
430        value[0] = (double)n1;
431        indices[0] = 0;
432        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
433                return retVal;
434        value[0] = (double)n2;
435        indices[0] = 1;
436        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
437                return retVal;
438        value[0] = (double)isOn;
439        indices[0] = 2;
440        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
441                return retVal;
442        value[0] = (double)NScatterEvents;
443        indices[0] = 3;
444        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
445                return retVal;
446        value[0] = (double)NSingleCoherent;
447        indices[0] = 4;
448        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
449                return retVal;
450        value[0] = (double)NMultipleCoherent;
451        indices[0] = 5;
452        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
453                return retVal;
454        value[0] = (double)NMultipleScatter;
455        indices[0] = 6;
456        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
457                return retVal; 
458        value[0] = (double)NCoherentEvents;
459        indices[0] = 7;
460        if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value))
461                return retVal; 
462
463
464//      HSetState((Handle)wavH, hState);                //release the handle of the 2D data wave
465//      WaveHandleModified(wavH);                       // Inform Igor that we have changed the wave. (CALLBACK! needed, but not allowed in Threading)
466       
467        return(0);
468}
469////////        end of main function for calculating multiple scattering
470
471int
472FindPixel(double testQ, double testPhi, double lam, double sdd,
473                double pixSize, double xCtr, double yCtr, long *xPixel, long *yPixel) {
474
475        double theta,dy,dx,qx,qy,pi;
476        pi = 4.0*atan(1.0);     
477        //decompose to qx,qy
478        qx = testQ*cos(testPhi);
479        qy = testQ*sin(testPhi);
480
481        //convert qx,qy to pixel locations relative to # of pixels x, y from center
482        theta = 2.0*asin(qy*lam/4.0/pi);
483        dy = sdd*tan(theta);
484        *yPixel = round(yCtr + dy/pixSize);
485       
486        theta = 2.0*asin(qx*lam/4.0/pi);
487        dx = sdd*tan(theta);
488        *xPixel = round(xCtr + dx/pixSize);
489
490        //if on detector, return xPix and yPix values, otherwise -1
491        if(*yPixel > 127 || *yPixel < 0) {
492                *yPixel = -1;
493        }
494        if(*xPixel > 127 || *xPixel < 0) {
495                *xPixel = -1;
496        }
497       
498        return(0);
499}
500
501
502//calculates new direction (xyz) from an old direction
503//theta and phi don't change
504int
505NewDirection(double *vx, double *vy, double *vz, double theta, double phi) {
506       
507        int err=0;
508        double vx0,vy0,vz0;
509        double nx,ny,mag_xy,tx,ty,tz;
510       
511        //store old direction vector
512        vx0 = *vx;
513        vy0 = *vy;
514        vz0 = *vz;
515       
516        mag_xy = sqrt(vx0*vx0 + vy0*vy0);
517        if(mag_xy < 1e-12) {
518                //old vector lies along beam direction
519                nx = 0;
520                ny = 1;
521                tx = 1;
522                ty = 0;
523                tz = 0;
524        } else {
525                nx = -vy0 / mag_xy;
526                ny = vx0 / mag_xy;
527                tx = -vz0*vx0 / mag_xy;
528                ty = -vz0*vy0 / mag_xy;
529                tz = mag_xy ;
530        }
531       
532        //new scattered direction vector
533        *vx = cos(phi)*sin(theta)*tx + sin(phi)*sin(theta)*nx + cos(theta)*vx0;
534        *vy = cos(phi)*sin(theta)*ty + sin(phi)*sin(theta)*ny + cos(theta)*vy0;
535        *vz = cos(phi)*sin(theta)*tz + cos(theta)*vz0;
536       
537        return(err);
538}
539
540double
541path_len(double aval, double sig_tot) {
542       
543        double retval;
544       
545        retval = -1.0*log(1.0-aval)/sig_tot;
546       
547        return(retval);
548}
549
550
551#define IA 16807
552#define IM 2147483647
553#define AM (1.0/IM)
554#define IQ 127773
555#define IR 2836
556#define NTAB 32
557#define NDIV (1+(IM-1)/NTAB)
558#define EPS 1.2e-7
559#define RNMX (1.0-EPS)
560
561float ran1(long *idum)
562{
563        int j;
564        long k;
565        static long iy=0;
566        static long iv[NTAB];
567        float temp;
568
569        if (*idum <= 0 || !iy) {
570                if (-(*idum) < 1) *idum=1;
571                else *idum = -(*idum);
572                for (j=NTAB+7;j>=0;j--) {
573                        k=(*idum)/IQ;
574                        *idum=IA*(*idum-k*IQ)-IR*k;
575                        if (*idum < 0) *idum += IM;
576                        if (j < NTAB) iv[j] = *idum;
577                }
578                iy=iv[0];
579        }
580        k=(*idum)/IQ;
581        *idum=IA*(*idum-k*IQ)-IR*k;
582        if (*idum < 0) *idum += IM;
583        j=iy/NDIV;
584        iy=iv[j];
585        iv[j] = *idum;
586        if ((temp=AM*iy) > RNMX) return RNMX;
587        else return temp;
588}
589#undef IA
590#undef IM
591#undef AM
592#undef IQ
593#undef IR
594#undef NTAB
595#undef NDIV
596#undef EPS
597#undef RNMX
598
599////////////////////////
600#define MBIG 1000000000
601#define MSEED 161803398
602#define MZ 0
603#define FAC (1.0/MBIG)
604
605float ran3(long *idum)
606{
607        static int inext,inextp;
608        static long ma[56];
609        static int iff=0;
610        long mj,mk;
611        int i,ii,k;
612
613        if (*idum < 0 || iff == 0) {
614                iff=1;
615                mj=MSEED-(*idum < 0 ? -*idum : *idum);
616                mj %= MBIG;
617                ma[55]=mj;
618                mk=1;
619                for (i=1;i<=54;i++) {
620                        ii=(21*i) % 55;
621                        ma[ii]=mk;
622                        mk=mj-mk;
623                        if (mk < MZ) mk += MBIG;
624                        mj=ma[ii];
625                }
626                for (k=1;k<=4;k++)
627                        for (i=1;i<=55;i++) {
628                                ma[i] -= ma[1+(i+30) % 55];
629                                if (ma[i] < MZ) ma[i] += MBIG;
630                        }
631                inext=0;
632                inextp=31;
633                *idum=1;
634        }
635        if (++inext == 56) inext=1;
636        if (++inextp == 56) inextp=1;
637        mj=ma[inext]-ma[inextp];
638        if (mj < MZ) mj += MBIG;
639        ma[inext]=mj;
640        return mj*FAC;
641}
642#undef MBIG
643#undef MSEED
644#undef MZ
645#undef FAC
646
647
648// returns the interpolated point value in xx[0,n-1] that has the value x
649double locate_interp(double xx[], long n, double x)
650{
651        unsigned long ju,jm,jl,j;
652        int ascnd;
653        double pt;
654
655//      char buf[256];
656       
657        jl=0;
658        ju=n-1;
659        ascnd=(xx[n-1] > xx[0]);
660        while (ju-jl > 1) {
661                jm=(ju+jl) >> 1;
662                if (x > xx[jm] == ascnd)
663                        jl=jm;
664                else
665                        ju=jm;
666        }
667        j=jl;           // the point I want is between xx[j] and xx[j+1]
668        pt = (x- xx[j])/(xx[j+1] - xx[j]);              //fractional distance, using linear interpolation
669       
670//      sprintf(buf, "x = %g, j= %ld, pt = %g\r",x,j,pt);
671//      XOPNotice(buf);
672       
673        return(pt+(double)j);
674}
675
676
677
678
679/////////////////////////////
680/*      RegisterFunction()
681       
682        Igor calls this at startup time to find the address of the
683        XFUNCs added by this XOP. See XOP manual regarding "Direct XFUNCs".
684*/
685static long
686RegisterFunction()
687{
688        int funcIndex;
689
690        funcIndex = GetXOPItem(0);              // Which function is Igor asking about?
691        switch (funcIndex) {
692                case 0:                                         //
693                        return((long)Monte_SANSX);
694                        break;
695                case 1:                                         //
696                        return((long)Monte_SANSX2);
697                        break;
698
699        }
700        return(NIL);
701}
702
703/*      XOPEntry()
704
705        This is the entry point from the host application to the XOP for all messages after the
706        INIT message.
707*/
708static void
709XOPEntry(void)
710{       
711        long result = 0;
712
713        switch (GetXOPMessage()) {
714                case FUNCADDRS:
715                        result = RegisterFunction();
716                        break;
717        }
718        SetXOPResult(result);
719}
720
721/*      main(ioRecHandle)
722
723        This is the initial entry point at which the host application calls XOP.
724        The message sent by the host must be INIT.
725        main() does any necessary initialization and then sets the XOPEntry field of the
726        ioRecHandle to the address to be called for future messages.
727*/
728HOST_IMPORT void
729main(IORecHandle ioRecHandle)
730{       
731        XOPInit(ioRecHandle);                                                   // Do standard XOP initialization.
732        SetXOPEntry(XOPEntry);                                                  // Set entry point for future calls.
733       
734        if (igorVersion < 600)                          // Requires Igor Pro 6.00 or later.
735                SetXOPResult(OLD_IGOR);                 // OLD_IGOR is defined in WaveAccess.h and there are corresponding error strings in WaveAccess.r and WaveAccessWinCustom.rc.
736        else
737                SetXOPResult(0L);
738}
Note: See TracBrowser for help on using the repository browser.