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

Last change on this file since 437 was 437, checked in by srkline, 14 years ago

Change to calculation to avoid unnecessary calculations for transmitted neutrons. Marginal speedup. Done to be consistent with non-XOP version.

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