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

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

Added poor man's threading to the MonteCarlo? calculation.

My guess is that the ran() function from NR is not thread safe (it is non-reentrant). So I simply duplicated Monte_SANSX to Monte_SANSX2, where each incarnation uses a different random number generator, either ran1() or ran3(). This means that currently only two processors are supported. Not a big deal. At least it works.

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