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

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

Updated the MonteCarlo? code to allow 4 processors, but simply copying the function 4 times, and defining 4 different random number generators. Still can't figure out what the problem is with threading a single version, but not worth the effort. Copy/paste is way faster.

Also added some simple (non-optimized) calculations for using Debye's sphere method. These are largely undocumented at this point - so see the code. These are XOP versions of the old ipf code I've used in the past, and stripped of the now-obsolete AltiVec? code (I now lose the 4x speedup from the vectorization...)

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