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

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

I think this resolved my conflicts...

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