Changeset 711 for sans/XOP_Dev


Ignore:
Timestamp:
Jul 1, 2010 12:18:44 PM (12 years ago)
Author:
srkline
Message:

Added wavelength distribution as part of the simulation

Fixed bug where pixel locations of (-1,0) row or column were incorrectly rounded towards zero, and ended up on the detector in row 0 or column 0, effectively doubling the count there.

Separated the 4 functions into 4 separate files, for ease in checking that they are all the same. Now the only difference in each file should be the RNG used. And they had better be different - I'm almost positive now that cross-calling of the same RNG by different threads is responsible for the crashing.

Location:
sans/XOP_Dev/MonteCarlo
Files:
3 added
2 edited

Legend:

Unmodified
Added
Removed
  • sans/XOP_Dev/MonteCarlo/MonteCarlo.c

    r677 r711  
    88 */ 
    99 
     10// 
     11// uses ran3() 
     12// 
    1013 
    1114#include "XOPStandardHeaders.h"                 // Include ANSI headers, Mac headers, IgorXOP.h, XOP.h and XOPSupport.h 
     
    1518//static int gCallSpinProcess = 1;              // Set to 1 to all user abort (cmd dot) and background processing. 
    1619 
    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 // 
    4220 
    4321 
     
    7048        long NDoubleCoherent,NMultipleScatter,isOn,xCtr_long,yCtr_long; 
    7149        long NMultipleCoherent,NCoherentEvents; 
     50        double deltaLam,v1,v2,currWavelength,rsq,fac;           //for simulating wavelength distribution 
    7251 
    7352        // for accessing the 2D wave data, direct method (see the WaveAccess example XOP) 
     
    156135        sig_incoh = inputWave[9]; 
    157136        sig_sas = inputWave[10]; 
     137        deltaLam = inputWave[11]; 
    158138        xCtr_long = (long)(xCtr+0.5);   //round() not defined in VC8 
    159139        yCtr_long = (long)(yCtr+0.5);   // and is probably wrong to use anyways (returns double!) 
     
    257237                } while(rr>r1); 
    258238 
     239                //pick the wavelength out of the wavelength spread, approximate as a gaussian 
     240                // from NR - pg 288. Needs random # from [0,1]. del is deltaLam/lam (as FWHM) and the 
     241                // 2.35 converts to a gaussian std dev. 
     242                do { 
     243                        v1=2.0*ran3(&seed)-1.0; 
     244                        v2=2.0*ran3(&seed)-1.0; 
     245                        rsq=v1*v1+v2*v2; 
     246                } while (rsq >= 1.0 || rsq == 0.0); 
     247                fac=sqrt(-2.0*log(rsq)/rsq); 
     248                                 
     249                //              gset=v1*fac             //technically, I'm throwing away one of the two values 
     250                currWavelength = (v2*fac)*deltaLam*wavelength/2.35 + wavelength; 
     251                 
    259252                do {   //Scattering Loop, will exit when "done" == 1 
    260253                                // keep scattering multiple times until the neutron exits the sample 
     
    297290                                                 
    298291                                                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 
     292                                                theta = q0/2.0/pi*currWavelength;               //SAS approximation. 1% error at theta=30 degrees 
    300293                                                 
    301294                                                find_theta = 1;         //always accept 
     
    339332                                if( index != 0) {               //neutron was scattered, figure out where it went 
    340333                                        theta_z = acos(vz);             // Angle (= 2theta) WITH respect to z axis. 
    341                                         testQ = 2.0*pi*sin(theta_z)/wavelength; 
     334                                        testQ = 2.0*pi*sin(theta_z)/currWavelength; 
    342335                                         
    343336                                        // pick a random phi angle, and see if it lands on the detector 
     
    347340                                         
    348341                                        // is it on the detector?        
    349                                         FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 
     342                                        FindPixel(testQ,testPhi,currWavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 
    350343                                                                                                         
    351344                                        if(xPixel != -1 && yPixel != -1) { 
     
    471464} 
    472465////////        end of main function for calculating multiple scattering 
    473  
    474 int 
    475 FindPixel(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 
    507 int 
    508 NewDirection(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  
    543 double 
    544 path_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  
    564 float 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  
    614 float 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  
    660 float 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  
    708 float 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 
    752 double 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 */ 
    788 static long 
    789 RegisterFunction() 
    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 */ 
    819 static void 
    820 XOPEntry(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 */ 
    839 HOST_IMPORT void 
    840 main(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 } 
  • sans/XOP_Dev/MonteCarlo/MonteCarlo2.c

    r673 r711  
    11/* 
    2  *  MonteCarlo.c 
    3  *  SANSAnalysis 
     2 *  MonteCarlo2.c 
     3 *  SANSMonteCarlo 
    44 * 
    5  *  Created by Steve Kline on 10/16/08. 
    6  *  Copyright 2008 __MyCompanyName__. All rights reserved. 
     5 *  Created by Steve Kline on 7/1/10. 
     6 *  Copyright 2010 NCNR. All rights reserved. 
    77 * 
    88 */ 
     
    5151        long NDoubleCoherent,NMultipleScatter,isOn,xCtr_long,yCtr_long; 
    5252        long NMultipleCoherent,NCoherentEvents; 
     53        double deltaLam,v1,v2,currWavelength,rsq,fac;           //for simulating wavelength distribution 
    5354 
    5455 
     
    138139        sig_incoh = inputWave[9]; 
    139140        sig_sas = inputWave[10]; 
     141        deltaLam = inputWave[11]; 
    140142        xCtr_long = (long)(xCtr+0.5); 
    141143        yCtr_long = (long)(yCtr+0.5); 
     
    239241                } while(rr>r1); 
    240242 
     243                //pick the wavelength out of the wavelength spread, approximate as a gaussian 
     244                // from NR - pg 288. Needs random # from [0,1]. del is deltaLam/lam (as FWHM) and the 
     245                // 2.35 converts to a gaussian std dev. 
     246                do { 
     247                        v1=2.0*ran1(&seed)-1.0; 
     248                        v2=2.0*ran1(&seed)-1.0; 
     249                        rsq=v1*v1+v2*v2; 
     250                } while (rsq >= 1.0 || rsq == 0.0); 
     251                fac=sqrt(-2.0*log(rsq)/rsq); 
     252                 
     253                //              gset=v1*fac             //technically, I'm throwing away one of the two values 
     254                currWavelength = (v2*fac)*deltaLam*wavelength/2.35 + wavelength; 
     255                 
    241256                do {   //Scattering Loop, will exit when "done" == 1 
    242257                                // keep scattering multiple times until the neutron exits the sample 
     
    279294                                                 
    280295                                                q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran1(&seed))*delta; 
    281                                                 theta = q0/2/pi*wavelength;             //SAS approximation 
     296                                                theta = q0/2/pi*currWavelength;         //SAS approximation 
    282297                                                 
    283298                                                find_theta = 1;         //always accept 
     
    321336                                if( index != 0) {               //neutron was scattered, figure out where it went 
    322337                                        theta_z = acos(vz);             // Angle (= 2theta) WITH respect to z axis. 
    323                                         testQ = 2*pi*sin(theta_z)/wavelength; 
     338                                        testQ = 2*pi*sin(theta_z)/currWavelength; 
    324339                                         
    325340                                        // pick a random phi angle, and see if it lands on the detector 
     
    329344                                         
    330345                                        // is it on the detector?        
    331                                         FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 
     346                                        FindPixel(testQ,testPhi,currWavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 
    332347                                                                                                         
    333348                                        if(xPixel != -1 && yPixel != -1) { 
     
    451466        return(0); 
    452467} 
    453 ////////        end of X2 
    454  
    455  
    456  
    457 ////////////////                X3 using ran3a 
    458 int 
    459 Monte_SANSX3(MC_ParamsPtr p) { 
    460         double *inputWave;                              /* pointer to double precision wave data */ 
    461         double *ran_dev;                                /* pointer to double precision wave data */ 
    462         double *nt;                             /* pointer to double precision wave data */ 
    463         double *j1;                             /* pointer to double precision wave data */ 
    464         double *j2;                             /* pointer to double precision wave data */ 
    465         double *nn;                             /* pointer to double precision wave data */ 
    466 //      double *MC_linear_data;                         /* pointer to double precision wave data */ 
    467         double *results;                                /* pointer to double precision wave data */ 
    468         double retVal;                          //return value 
    469  
    470         long imon; 
    471         double r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas; 
    472         long ind,index,n_index; 
    473         double qmax,theta_max,q0,zpow; 
    474         long n1,n2,n3; 
    475         double dth,zz,xx,yy,phi; 
    476         double theta,ran,ll,rr; 
    477         long done,find_theta,err;               //used as logicals 
    478         long xPixel,yPixel; 
    479         double vx,vy,vz,theta_z; 
    480         double sig_abs,ratio,sig_total; 
    481         double testQ,testPhi,left,delta,dummy,pi; 
    482         double sigabs_0,num_bins; 
    483         long NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent; 
    484         long NDoubleCoherent,NMultipleScatter,isOn,xCtr_long,yCtr_long; 
    485         long NMultipleCoherent,NCoherentEvents; 
    486  
    487  
    488         // for accessing the 2D wave data, direct method (see the WaveAccess example XOP) 
    489         waveHndl wavH; 
    490 //      int waveType,hState; 
    491         long numDimensions; 
    492         long dimensionSizes[MAX_DIMENSIONS+1]; 
    493 //      char* dataStartPtr; 
    494 //      long dataOffset; 
    495 //      long numRows, numColumns; 
    496         long numRows_ran_dev; 
    497 //      double *dp0, *dp; 
    498         double value[2];                                // Pointers used for double data. 
    499         long seed; 
    500         long indices[MAX_DIMENSIONS]; 
    501          
    502 //      char buf[256]; 
    503                  
    504         /* check that wave handles are all valid */ 
    505         if (p->inputWaveH == NIL) { 
    506                 SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */ 
    507                 return(NON_EXISTENT_WAVE); 
    508         } 
    509         if (p->ran_devH == NIL) { 
    510                 SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */ 
    511                 return(NON_EXISTENT_WAVE); 
    512         }        
    513         if (p->ntH == NIL) { 
    514                 SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */ 
    515                 return(NON_EXISTENT_WAVE); 
    516         } 
    517         if (p->j1H == NIL) { 
    518                 SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */ 
    519                 return(NON_EXISTENT_WAVE); 
    520         } 
    521         if (p->j2H == NIL) { 
    522                 SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */ 
    523                 return(NON_EXISTENT_WAVE); 
    524         } 
    525         if (p->nnH == NIL) { 
    526                 SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */ 
    527                 return(NON_EXISTENT_WAVE); 
    528         } 
    529         if (p->MC_linear_dataH == NIL) { 
    530                 SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */ 
    531                 return(NON_EXISTENT_WAVE); 
    532         } 
    533         if (p->resultsH == NIL) { 
    534                 SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */ 
    535                 return(NON_EXISTENT_WAVE); 
    536         } 
    537          
    538         p->retVal = 0; 
    539          
    540 // trusting that all inputs are DOUBLE PRECISION WAVES!!! 
    541         inputWave = WaveData(p->inputWaveH); 
    542         ran_dev = WaveData(p->ran_devH); 
    543         nt = WaveData(p->ntH); 
    544         j1 = WaveData(p->j1H); 
    545         j2 = WaveData(p->j2H); 
    546         nn = WaveData(p->nnH); 
    547 //      MC_linear_data = WaveData(p->MC_linear_dataH); 
    548         results = WaveData(p->resultsH); 
    549          
    550         seed = (long)results[0]; 
    551          
    552 //      sprintf(buf, "input seed = %ld\r", seed); 
    553 //      XOPNotice(buf); 
    554          
    555         if(seed >= 0) { 
    556                 seed = -1234509876; 
    557         } 
    558  
    559         dummy = ran3a(&seed);           //initialize the random sequence by passing in a negative value 
    560         seed = 12348765;                //non-negative after that does nothing 
    561  
    562         imon = (int)inputWave[0]; 
    563         r1 = inputWave[1]; 
    564         r2 = inputWave[2]; 
    565         xCtr = inputWave[3]; 
    566         yCtr = inputWave[4]; 
    567         sdd = inputWave[5]; 
    568         pixSize = inputWave[6]; 
    569         thick = inputWave[7]; 
    570         wavelength = inputWave[8]; 
    571         sig_incoh = inputWave[9]; 
    572         sig_sas = inputWave[10]; 
    573         xCtr_long = (long)(xCtr+0.5); 
    574         yCtr_long = (long)(yCtr+0.5); 
    575          
    576         dummy = MDGetWaveScaling(p->ran_devH, 0, &delta, &left);                //0 is the rows 
    577         if (retVal = MDGetWaveDimensions(p->ran_devH, &numDimensions, dimensionSizes)) 
    578                 return retVal; 
    579         numRows_ran_dev = dimensionSizes[0]; 
    580          
    581         pi = 4.0*atan(1.0);      
    582          
    583         // access the 2D wave data for writing using the direct method 
    584         wavH = p->MC_linear_dataH; 
    585         if (wavH == NIL) 
    586                 return NOWAV; 
    587  
    588 //      waveType = WaveType(wavH); 
    589 //      if (waveType & NT_CMPLX) 
    590 //              return NO_COMPLEX_WAVE; 
    591 //      if (waveType==TEXT_WAVE_TYPE) 
    592 //              return NUMERIC_ACCESS_ON_TEXT_WAVE; 
    593 //      if (retVal = MDGetWaveDimensions(wavH, &numDimensions, dimensionSizes)) 
    594 //              return retVal; 
    595 //      numRows = dimensionSizes[0]; 
    596 //      numColumns = dimensionSizes[1]; 
    597          
    598 //      if (retVal = MDAccessNumericWaveData(wavH, kMDWaveAccessMode0, &dataOffset)) 
    599 //              return retVal; 
    600                  
    601 //      hState = MoveLockHandle(wavH);          // So wave data can't move. Remember to call HSetState when done. 
    602 //      dataStartPtr = (char*)(*wavH) + dataOffset; 
    603 //      dp0 = (double*)dataStartPtr;                    // Pointer to the start of the 2D wave data. 
    604          
    605 //scattering power and maximum qvalue to bin 
    606 //      zpow = .1               //scattering power, calculated below 
    607         qmax = 4.0*pi/wavelength;       //maximum Q to bin 1D data. (A-1) (not really used) 
    608         sigabs_0 = 0.0;         // ignore absorption cross section/wavelength [1/(cm A)] 
    609         n_index = 50;   // maximum number of scattering events per neutron 
    610         num_bins = 200;         //number of 1-D bins (not really used) 
    611          
    612 //c       total SAS cross-section 
    613 // 
    614         zpow = sig_sas*thick;                   //since I now calculate the sig_sas from the model 
    615         sig_abs = sigabs_0 * wavelength; 
    616         sig_total = sig_abs + sig_sas + sig_incoh; 
    617 //      Print "The TOTAL XSECTION. (CM-1) is ",sig_total 
    618 //      Print "The TOTAL SAS XSECTION. (CM-1) is ",sig_sas 
    619 //      results[0] = sig_total; 
    620 //      results[1] = sig_sas; 
    621 //      RATIO = SIG_ABS / SIG_TOTAL 
    622         ratio = sig_incoh / sig_total; 
    623          
    624         theta_max = wavelength*qmax/(2*pi); 
    625 //C     SET Theta-STEP SIZE. 
    626         dth = theta_max/num_bins; 
    627 //      Print "theta bin size = dth = ",dth 
    628  
    629 //C     INITIALIZE COUNTERS. 
    630         n1 = 0; 
    631         n2 = 0; 
    632         n3 = 0; 
    633         NSingleIncoherent = 0; 
    634         NSingleCoherent = 0; 
    635         NDoubleCoherent = 0; 
    636         NMultipleScatter = 0; 
    637         NScatterEvents = 0; 
    638         NMultipleCoherent = 0; 
    639         NCoherentEvents = 0; 
    640          
    641         isOn = 0; 
    642          
    643 //C     MONITOR LOOP - looping over the number of incedent neutrons 
    644 //note that zz, is the z-position in the sample - NOT the scattering power 
    645 // NOW, start the loop, throwing neutrons at the sample. 
    646         do { 
    647                 ////SpinProcess() IS A CALLBACK, and not good for Threading! 
    648 //              if ((n1 % 1000 == 0) && gCallSpinProcess && SpinProcess()) {            // Spins cursor and allows background processing. 
    649 //                              retVal = -1;                                                            // User aborted. 
    650 //                              break; 
    651 //              } 
    652          
    653                 vx = 0.0;                       // Initialize direction vector. 
    654                 vy = 0.0; 
    655                 vz = 1.0; 
    656                  
    657                 theta = 0.0;            //      Initialize scattering angle. 
    658                 phi = 0.0;                      //      Intialize azimuthal angle. 
    659                 n1 += 1;                        //      Increment total number neutrons counter. 
    660                 done = 0;                       //      True when neutron is absorbed or when  scattered out of the sample. 
    661                 index = 0;                      //      Set counter for number of scattering events. 
    662                 zz = 0.0;                       //      Set entering dimension of sample. 
    663                 incoherentEvent = 0; 
    664                 coherentEvent = 0; 
    665                  
    666                 do      {                               //      Makes sure position is within circle. 
    667                         ran = ran3a(&seed);             //[0,1] 
    668                         xx = 2.0*r1*(ran-0.5);          //X beam position of neutron entering sample. 
    669                         ran = ran3a(&seed);             //[0,1] 
    670                         yy = 2.0*r1*(ran-0.5);          //Y beam position ... 
    671                         rr = sqrt(xx*xx+yy*yy);         //Radial position of neutron in incident beam. 
    672                 } while(rr>r1); 
    673  
    674                 do {   //Scattering Loop, will exit when "done" == 1 
    675                                 // keep scattering multiple times until the neutron exits the sample 
    676                         ran = ran3a(&seed);             //[0,1]  RANDOM NUMBER FOR DETERMINING PATH LENGTH 
    677                         ll = path_len(ran,sig_total); 
    678                         //Determine new scattering direction vector. 
    679                         err = NewDirection(&vx,&vy,&vz,theta,phi);              //vx,vy,vz updated, theta, phi unchanged by function 
    680                                                                          
    681                         //X,Y,Z-POSITION OF SCATTERING EVENT. 
    682                         xx += ll*vx; 
    683                         yy += ll*vy; 
    684                         zz += ll*vz; 
    685                         rr = sqrt(xx*xx+yy*yy);         //radial position of scattering event. 
    686  
    687                         //sprintf(buf, "xx,yy,zz,vx,vy,vz,ll = %g %g %g %g %g %g %g\r",xx,yy,zz,vx,vy,vz,ll); 
    688                         //XOPNotice(buf); 
    689                                                  
    690                         //Check whether interaction occurred within sample volume. 
    691                         if (((zz > 0.0) && (zz < thick)) && (rr < r2)) { 
    692                                 //NEUTRON INTERACTED. 
    693                                 //sprintf(buf,"neutron interacted\r"); 
    694                                 //XOPNotice(buf); 
    695                                  
    696                                 index += 1;                     //Increment counter of scattering events. 
    697                                 if (index == 1) { 
    698                                         n2 += 1;                //Increment # of scat. neutrons 
    699                                 } 
    700                                 ran = ran3a(&seed);             //[0,1] 
    701                                 //Split neutron interactions into scattering and absorption events 
    702                                 if (ran > ratio ) {             //C             NEUTRON SCATTERED coherently 
    703                                         //sprintf(buf,"neutron scatters coherently\r"); 
    704                                         //XOPNotice(buf); 
    705                                         coherentEvent += 1; 
    706                                         find_theta = 0;                 //false 
    707                                         do { 
    708                                                 // pick a q-value from the deviate function 
    709                                                 // pnt2x truncates the point to an integer before returning the x 
    710                                                 // so get it from the wave scaling instead 
    711 //                                              q0 =left + binarysearchinterp(ran_dev,ran3a(seed))*delta; 
    712                                                  
    713                                                 q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran3a(&seed))*delta; 
    714                                                 theta = q0/2/pi*wavelength;             //SAS approximation 
    715                                                  
    716                                                 find_theta = 1;         //always accept 
    717  
    718                                                 //sprintf(buf, "after locate_interp call q0 = %g, theta = %g,left = %g,delta = %g\r",q0,theta,left,delta); 
    719                                                 //XOPNotice(buf); 
    720  
    721                                         } while(!find_theta); 
    722                                          
    723                                         ran = ran3a(&seed);             //[0,1] 
    724                                         phi = 2.0*pi*ran;                       //Chooses azimuthal scattering angle. 
    725                                 } else { 
    726                                         //NEUTRON scattered incoherently 
    727                                         //sprintf(buf,"neutron scatters incoherent\r"); 
    728                                         //XOPNotice(buf); 
    729                                         incoherentEvent += 1; 
    730                                   // phi and theta are random over the entire sphere of scattering 
    731                                         // !can't just choose random theta and phi, won't be random over sphere solid angle 
    732  
    733                                         ran = ran3a(&seed);             //[0,1] 
    734                                         theta = acos(2.0*ran-1); 
    735                          
    736                                         ran = ran3a(&seed);             //[0,1] 
    737                                         phi = 2.0*pi*ran;                       //Chooses azimuthal scattering angle. 
    738                                 }               //(ran > ratio) 
    739                         } else { 
    740                                 //NEUTRON ESCAPES FROM SAMPLE -- bin it somewhere                                                                
    741                                 done = 1;               //done = true, will exit from loop 
    742                                 //Increment #scattering events array 
    743                                 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 
    744                                 indices[0] =index;                      //this sets access to nn[index] 
    745                                 if (index <= n_index) { 
    746                                         if (retVal = MDGetNumericWavePointValue(p->nnH, indices, value)) 
    747                                                 return retVal; 
    748                                         value[0] += 1; // add one to the value 
    749                                         if (retVal = MDSetNumericWavePointValue(p->nnH, indices, value)) 
    750                                                 return retVal; 
    751                                 //      nn[index] += 1; 
    752                                 } 
    753                                                                                                  
    754                                 if( index != 0) {               //neutron was scattered, figure out where it went 
    755                                         theta_z = acos(vz);             // Angle (= 2theta) WITH respect to z axis. 
    756                                         testQ = 2*pi*sin(theta_z)/wavelength; 
    757                                          
    758                                         // pick a random phi angle, and see if it lands on the detector 
    759                                         // since the scattering is isotropic, I can safely pick a new, random value 
    760                                         // this would not be true if simulating anisotropic scattering. 
    761                                         testPhi = ran3a(&seed)*2*pi; 
    762                                          
    763                                         // is it on the detector?        
    764                                         FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 
    765                                                                                                          
    766                                         if(xPixel != -1 && yPixel != -1) { 
    767                                                 isOn += 1; 
    768                                                 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 
    769                                                 indices[0] = xPixel; 
    770                                                 indices[1] = yPixel; 
    771                                                 if (retVal = MDGetNumericWavePointValue(wavH, indices, value)) 
    772                                                         return retVal; 
    773                                                 value[0] += 1; // Real part 
    774                                                 if (retVal = MDSetNumericWavePointValue(wavH, indices, value)) 
    775                                                         return retVal; 
    776                                                 //if(index==1)  // only the single scattering events 
    777                                                         //dp = dp0 + xPixel + yPixel*numColumns;                //offset the pointer to the exact memory location 
    778                                                         //*dp += 1;             //increment the value there 
    779                                                 //endif 
    780                                         } 
    781                                          
    782  
    783         /*              is this causing me a problem since I'm not locking these? Probably not, since it crashes even if I comment these out... */ 
    784                                         if(theta_z < theta_max) { 
    785                                                 //Choose index for scattering angle array. 
    786                                                 //IND = NINT(THETA_z/DTH + 0.4999999) 
    787                                                 ind = (long)(theta_z/dth + 0.4999999);          //round is eqivalent to nint() 
    788                                                 nt[ind] += 1;                   //Increment bin for angle. 
    789                                                 //Increment angle array for single scattering events. 
    790                                                 if (index == 1) { 
    791                                                         j1[ind] += 1; 
    792                                                 } 
    793                                                 //Increment angle array for double scattering events. 
    794                                                 if (index == 2) { 
    795                                                         j2[ind] += 1; 
    796                                                 } 
    797                                         } 
    798         /**/ 
    799                                          
    800                                         // increment all of the counters now since done==1 here and I'm sure to exit and get another neutron 
    801                                         NScatterEvents += index;                //total number of scattering events 
    802                                         if(index == 1 && incoherentEvent == 1) { 
    803                                                 NSingleIncoherent += 1; 
    804                                         } 
    805                                         if(index == 1 && coherentEvent == 1) { 
    806                                                 NSingleCoherent += 1; 
    807                                         } 
    808                                         if(index == 2 && coherentEvent == 1 && incoherentEvent == 0) { 
    809                                                 NDoubleCoherent += 1; 
    810                                         } 
    811                                         if(index > 1) { 
    812                                                 NMultipleScatter += 1; 
    813                                         } 
    814                                         if(coherentEvent >= 1 && incoherentEvent == 0) { 
    815                                                 NCoherentEvents += 1; 
    816                                         } 
    817                                         if(coherentEvent > 1 && incoherentEvent == 0) { 
    818                                                 NMultipleCoherent += 1; 
    819                                         } 
    820  
    821                                 } else {        // index was zero, neutron must be transmitted, so just increment the proper counters and data 
    822                                                 isOn += 1; 
    823                                                 nt[0] += 1; 
    824                                                 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 
    825                                                 //indices[0] = xCtr_long;               //don't put everything in one pixel 
    826                                                 //indices[1] = yCtr_long; 
    827                                                 indices[0] = (long)(xCtr+xx/pixSize+0.5); 
    828                                                 indices[1] = (long)(yCtr+yy/pixSize+0.5); 
    829                                                 // check for valid indices - got an XOP error, probably from here 
    830                                                 if(indices[0] > 127) indices[0] = 127; 
    831                                                 if(indices[0] < 0) indices[0] = 0; 
    832                                                 if(indices[1] > 127) indices[1] = 127; 
    833                                                 if(indices[1] < 0) indices[1] = 0; 
    834                                                  
    835                                                 if (retVal = MDGetNumericWavePointValue(wavH, indices, value)) 
    836                                                         return retVal; 
    837                                                 value[0] += 1; // Real part 
    838                                                 if (retVal = MDSetNumericWavePointValue(wavH, indices, value)) 
    839                                                         return retVal; 
    840                                         }        
    841                         } 
    842                  } while (!done); 
    843         } while(n1 < imon); 
    844          
    845 // assign the results to the wave 
    846  
    847         MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 
    848         value[0] = (double)n1; 
    849         indices[0] = 0; 
    850         if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 
    851                 return retVal; 
    852         value[0] = (double)n2; 
    853         indices[0] = 1; 
    854         if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 
    855                 return retVal; 
    856         value[0] = (double)isOn; 
    857         indices[0] = 2; 
    858         if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 
    859                 return retVal; 
    860         value[0] = (double)NScatterEvents; 
    861         indices[0] = 3; 
    862         if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 
    863                 return retVal; 
    864         value[0] = (double)NSingleCoherent; 
    865         indices[0] = 4; 
    866         if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 
    867                 return retVal; 
    868         value[0] = (double)NMultipleCoherent; 
    869         indices[0] = 5; 
    870         if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 
    871                 return retVal; 
    872         value[0] = (double)NMultipleScatter; 
    873         indices[0] = 6; 
    874         if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 
    875                 return retVal;   
    876         value[0] = (double)NCoherentEvents; 
    877         indices[0] = 7; 
    878         if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 
    879                 return retVal; 
    880  
    881 //      HSetState((Handle)wavH, hState);                //release the handle of the 2D data wave 
    882 //      WaveHandleModified(wavH);                       // Inform Igor that we have changed the wave. (CALLBACK! needed, but not allowed in Threading) 
    883          
    884         return(0); 
    885 } 
    886 ////////        end of X3 
    887  
    888  
    889 ///////////////// X4 using ran1a() 
    890 int 
    891 Monte_SANSX4(MC_ParamsPtr p) { 
    892         double *inputWave;                              /* pointer to double precision wave data */ 
    893         double *ran_dev;                                /* pointer to double precision wave data */ 
    894         double *nt;                             /* pointer to double precision wave data */ 
    895         double *j1;                             /* pointer to double precision wave data */ 
    896         double *j2;                             /* pointer to double precision wave data */ 
    897         double *nn;                             /* pointer to double precision wave data */ 
    898 //      double *MC_linear_data;                         /* pointer to double precision wave data */ 
    899         double *results;                                /* pointer to double precision wave data */ 
    900         double retVal;                          //return value 
    901  
    902         long imon; 
    903         double r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas; 
    904         long ind,index,n_index; 
    905         double qmax,theta_max,q0,zpow; 
    906         long n1,n2,n3; 
    907         double dth,zz,xx,yy,phi; 
    908         double theta,ran,ll,rr; 
    909         long done,find_theta,err;               //used as logicals 
    910         long xPixel,yPixel; 
    911         double vx,vy,vz,theta_z; 
    912         double sig_abs,ratio,sig_total; 
    913         double testQ,testPhi,left,delta,dummy,pi; 
    914         double sigabs_0,num_bins; 
    915         long NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent; 
    916         long NDoubleCoherent,NMultipleScatter,isOn,xCtr_long,yCtr_long; 
    917         long NMultipleCoherent,NCoherentEvents; 
    918  
    919  
    920         // for accessing the 2D wave data, direct method (see the WaveAccess example XOP) 
    921         waveHndl wavH; 
    922 //      int waveType,hState; 
    923         long numDimensions; 
    924         long dimensionSizes[MAX_DIMENSIONS+1]; 
    925 //      char* dataStartPtr; 
    926 //      long dataOffset; 
    927 //      long numRows, numColumns; 
    928         long numRows_ran_dev; 
    929 //      double *dp0, *dp; 
    930         double value[2];                                // Pointers used for double data. 
    931         long seed; 
    932         long indices[MAX_DIMENSIONS]; 
    933          
    934 //      char buf[256]; 
    935                  
    936         /* check that wave handles are all valid */ 
    937         if (p->inputWaveH == NIL) { 
    938                 SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */ 
    939                 return(NON_EXISTENT_WAVE); 
    940         } 
    941         if (p->ran_devH == NIL) { 
    942                 SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */ 
    943                 return(NON_EXISTENT_WAVE); 
    944         }        
    945         if (p->ntH == NIL) { 
    946                 SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */ 
    947                 return(NON_EXISTENT_WAVE); 
    948         } 
    949         if (p->j1H == NIL) { 
    950                 SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */ 
    951                 return(NON_EXISTENT_WAVE); 
    952         } 
    953         if (p->j2H == NIL) { 
    954                 SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */ 
    955                 return(NON_EXISTENT_WAVE); 
    956         } 
    957         if (p->nnH == NIL) { 
    958                 SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */ 
    959                 return(NON_EXISTENT_WAVE); 
    960         } 
    961         if (p->MC_linear_dataH == NIL) { 
    962                 SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */ 
    963                 return(NON_EXISTENT_WAVE); 
    964         } 
    965         if (p->resultsH == NIL) { 
    966                 SetNaN64(&p->retVal);                                   /* return NaN if wave is not valid */ 
    967                 return(NON_EXISTENT_WAVE); 
    968         } 
    969          
    970         p->retVal = 0; 
    971          
    972 // trusting that all inputs are DOUBLE PRECISION WAVES!!! 
    973         inputWave = WaveData(p->inputWaveH); 
    974         ran_dev = WaveData(p->ran_devH); 
    975         nt = WaveData(p->ntH); 
    976         j1 = WaveData(p->j1H); 
    977         j2 = WaveData(p->j2H); 
    978         nn = WaveData(p->nnH); 
    979 //      MC_linear_data = WaveData(p->MC_linear_dataH); 
    980         results = WaveData(p->resultsH); 
    981          
    982         seed = (long)results[0]; 
    983          
    984 //      sprintf(buf, "input seed = %ld\r", seed); 
    985 //      XOPNotice(buf); 
    986          
    987         if(seed >= 0) { 
    988                 seed = -1234509876; 
    989         } 
    990  
    991         dummy = ran1a(&seed);           //initialize the random sequence by passing in a negative value 
    992         seed = 12348765;                //non-negative after that does nothing 
    993  
    994         imon = (int)inputWave[0]; 
    995         r1 = inputWave[1]; 
    996         r2 = inputWave[2]; 
    997         xCtr = inputWave[3]; 
    998         yCtr = inputWave[4]; 
    999         sdd = inputWave[5]; 
    1000         pixSize = inputWave[6]; 
    1001         thick = inputWave[7]; 
    1002         wavelength = inputWave[8]; 
    1003         sig_incoh = inputWave[9]; 
    1004         sig_sas = inputWave[10]; 
    1005         xCtr_long = (long)(xCtr+0.5); 
    1006         yCtr_long = (long)(yCtr+0.5); 
    1007          
    1008         dummy = MDGetWaveScaling(p->ran_devH, 0, &delta, &left);                //0 is the rows 
    1009         if (retVal = MDGetWaveDimensions(p->ran_devH, &numDimensions, dimensionSizes)) 
    1010                 return retVal; 
    1011         numRows_ran_dev = dimensionSizes[0]; 
    1012          
    1013         pi = 4.0*atan(1.0);      
    1014          
    1015         // access the 2D wave data for writing using the direct method 
    1016         wavH = p->MC_linear_dataH; 
    1017         if (wavH == NIL) 
    1018                 return NOWAV; 
    1019  
    1020 //      waveType = WaveType(wavH); 
    1021 //      if (waveType & NT_CMPLX) 
    1022 //              return NO_COMPLEX_WAVE; 
    1023 //      if (waveType==TEXT_WAVE_TYPE) 
    1024 //              return NUMERIC_ACCESS_ON_TEXT_WAVE; 
    1025 //      if (retVal = MDGetWaveDimensions(wavH, &numDimensions, dimensionSizes)) 
    1026 //              return retVal; 
    1027 //      numRows = dimensionSizes[0]; 
    1028 //      numColumns = dimensionSizes[1]; 
    1029          
    1030 //      if (retVal = MDAccessNumericWaveData(wavH, kMDWaveAccessMode0, &dataOffset)) 
    1031 //              return retVal; 
    1032                  
    1033 //      hState = MoveLockHandle(wavH);          // So wave data can't move. Remember to call HSetState when done. 
    1034 //      dataStartPtr = (char*)(*wavH) + dataOffset; 
    1035 //      dp0 = (double*)dataStartPtr;                    // Pointer to the start of the 2D wave data. 
    1036          
    1037 //scattering power and maximum qvalue to bin 
    1038 //      zpow = .1               //scattering power, calculated below 
    1039         qmax = 4.0*pi/wavelength;       //maximum Q to bin 1D data. (A-1) (not really used) 
    1040         sigabs_0 = 0.0;         // ignore absorption cross section/wavelength [1/(cm A)] 
    1041         n_index = 50;   // maximum number of scattering events per neutron 
    1042         num_bins = 200;         //number of 1-D bins (not really used) 
    1043          
    1044 //c       total SAS cross-section 
    1045 // 
    1046         zpow = sig_sas*thick;                   //since I now calculate the sig_sas from the model 
    1047         sig_abs = sigabs_0 * wavelength; 
    1048         sig_total = sig_abs + sig_sas + sig_incoh; 
    1049 //      Print "The TOTAL XSECTION. (CM-1) is ",sig_total 
    1050 //      Print "The TOTAL SAS XSECTION. (CM-1) is ",sig_sas 
    1051 //      results[0] = sig_total; 
    1052 //      results[1] = sig_sas; 
    1053 //      RATIO = SIG_ABS / SIG_TOTAL 
    1054         ratio = sig_incoh / sig_total; 
    1055          
    1056         theta_max = wavelength*qmax/(2*pi); 
    1057 //C     SET Theta-STEP SIZE. 
    1058         dth = theta_max/num_bins; 
    1059 //      Print "theta bin size = dth = ",dth 
    1060  
    1061 //C     INITIALIZE COUNTERS. 
    1062         n1 = 0; 
    1063         n2 = 0; 
    1064         n3 = 0; 
    1065         NSingleIncoherent = 0; 
    1066         NSingleCoherent = 0; 
    1067         NDoubleCoherent = 0; 
    1068         NMultipleScatter = 0; 
    1069         NScatterEvents = 0; 
    1070         NMultipleCoherent = 0; 
    1071         NCoherentEvents = 0; 
    1072          
    1073         isOn = 0; 
    1074          
    1075 //C     MONITOR LOOP - looping over the number of incedent neutrons 
    1076 //note that zz, is the z-position in the sample - NOT the scattering power 
    1077 // NOW, start the loop, throwing neutrons at the sample. 
    1078         do { 
    1079                 ////SpinProcess() IS A CALLBACK, and not good for Threading! 
    1080 //              if ((n1 % 1000 == 0) && gCallSpinProcess && SpinProcess()) {            // Spins cursor and allows background processing. 
    1081 //                              retVal = -1;                                                            // User aborted. 
    1082 //                              break; 
    1083 //              } 
    1084          
    1085                 vx = 0.0;                       // Initialize direction vector. 
    1086                 vy = 0.0; 
    1087                 vz = 1.0; 
    1088                  
    1089                 theta = 0.0;            //      Initialize scattering angle. 
    1090                 phi = 0.0;                      //      Intialize azimuthal angle. 
    1091                 n1 += 1;                        //      Increment total number neutrons counter. 
    1092                 done = 0;                       //      True when neutron is absorbed or when  scattered out of the sample. 
    1093                 index = 0;                      //      Set counter for number of scattering events. 
    1094                 zz = 0.0;                       //      Set entering dimension of sample. 
    1095                 incoherentEvent = 0; 
    1096                 coherentEvent = 0; 
    1097                  
    1098                 do      {                               //      Makes sure position is within circle. 
    1099                         ran = ran1a(&seed);             //[0,1] 
    1100                         xx = 2.0*r1*(ran-0.5);          //X beam position of neutron entering sample. 
    1101                         ran = ran1a(&seed);             //[0,1] 
    1102                         yy = 2.0*r1*(ran-0.5);          //Y beam position ... 
    1103                         rr = sqrt(xx*xx+yy*yy);         //Radial position of neutron in incident beam. 
    1104                 } while(rr>r1); 
    1105  
    1106                 do {   //Scattering Loop, will exit when "done" == 1 
    1107                                 // keep scattering multiple times until the neutron exits the sample 
    1108                         ran = ran1a(&seed);             //[0,1]  RANDOM NUMBER FOR DETERMINING PATH LENGTH 
    1109                         ll = path_len(ran,sig_total); 
    1110                         //Determine new scattering direction vector. 
    1111                         err = NewDirection(&vx,&vy,&vz,theta,phi);              //vx,vy,vz updated, theta, phi unchanged by function 
    1112                                                                          
    1113                         //X,Y,Z-POSITION OF SCATTERING EVENT. 
    1114                         xx += ll*vx; 
    1115                         yy += ll*vy; 
    1116                         zz += ll*vz; 
    1117                         rr = sqrt(xx*xx+yy*yy);         //radial position of scattering event. 
    1118  
    1119                         //sprintf(buf, "xx,yy,zz,vx,vy,vz,ll = %g %g %g %g %g %g %g\r",xx,yy,zz,vx,vy,vz,ll); 
    1120                         //XOPNotice(buf); 
    1121                                                  
    1122                         //Check whether interaction occurred within sample volume. 
    1123                         if (((zz > 0.0) && (zz < thick)) && (rr < r2)) { 
    1124                                 //NEUTRON INTERACTED. 
    1125                                 //sprintf(buf,"neutron interacted\r"); 
    1126                                 //XOPNotice(buf); 
    1127                                  
    1128                                 index += 1;                     //Increment counter of scattering events. 
    1129                                 if (index == 1) { 
    1130                                         n2 += 1;                //Increment # of scat. neutrons 
    1131                                 } 
    1132                                 ran = ran1a(&seed);             //[0,1] 
    1133                                 //Split neutron interactions into scattering and absorption events 
    1134                                 if (ran > ratio ) {             //C             NEUTRON SCATTERED coherently 
    1135                                         //sprintf(buf,"neutron scatters coherently\r"); 
    1136                                         //XOPNotice(buf); 
    1137                                         coherentEvent += 1; 
    1138                                         find_theta = 0;                 //false 
    1139                                         do { 
    1140                                                 // pick a q-value from the deviate function 
    1141                                                 // pnt2x truncates the point to an integer before returning the x 
    1142                                                 // so get it from the wave scaling instead 
    1143 //                                              q0 =left + binarysearchinterp(ran_dev,ran1a(seed))*delta; 
    1144                                                  
    1145                                                 q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran1a(&seed))*delta; 
    1146                                                 theta = q0/2/pi*wavelength;             //SAS approximation 
    1147                                                  
    1148                                                 find_theta = 1;         //always accept 
    1149  
    1150                                                 //sprintf(buf, "after locate_interp call q0 = %g, theta = %g,left = %g,delta = %g\r",q0,theta,left,delta); 
    1151                                                 //XOPNotice(buf); 
    1152  
    1153                                         } while(!find_theta); 
    1154                                          
    1155                                         ran = ran1a(&seed);             //[0,1] 
    1156                                         phi = 2.0*pi*ran;                       //Chooses azimuthal scattering angle. 
    1157                                 } else { 
    1158                                         //NEUTRON scattered incoherently 
    1159                                         //sprintf(buf,"neutron scatters incoherent\r"); 
    1160                                         //XOPNotice(buf); 
    1161                                         incoherentEvent += 1; 
    1162                                   // phi and theta are random over the entire sphere of scattering 
    1163                                         // !can't just choose random theta and phi, won't be random over sphere solid angle 
    1164  
    1165                                         ran = ran1a(&seed);             //[0,1] 
    1166                                         theta = acos(2.0*ran-1); 
    1167                          
    1168                                         ran = ran1a(&seed);             //[0,1] 
    1169                                         phi = 2.0*pi*ran;                       //Chooses azimuthal scattering angle. 
    1170                                 }               //(ran > ratio) 
    1171                         } else { 
    1172                                 //NEUTRON ESCAPES FROM SAMPLE -- bin it somewhere                                                                
    1173                                 done = 1;               //done = true, will exit from loop 
    1174                                 //Increment #scattering events array 
    1175                                 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 
    1176                                 indices[0] =index;                      //this sets access to nn[index] 
    1177                                 if (index <= n_index) { 
    1178                                         if (retVal = MDGetNumericWavePointValue(p->nnH, indices, value)) 
    1179                                                 return retVal; 
    1180                                         value[0] += 1; // add one to the value 
    1181                                         if (retVal = MDSetNumericWavePointValue(p->nnH, indices, value)) 
    1182                                                 return retVal; 
    1183                                 //      nn[index] += 1; 
    1184                                 } 
    1185                                                                                                  
    1186                                 if( index != 0) {               //neutron was scattered, figure out where it went 
    1187                                         theta_z = acos(vz);             // Angle (= 2theta) WITH respect to z axis. 
    1188                                         testQ = 2*pi*sin(theta_z)/wavelength; 
    1189                                          
    1190                                         // pick a random phi angle, and see if it lands on the detector 
    1191                                         // since the scattering is isotropic, I can safely pick a new, random value 
    1192                                         // this would not be true if simulating anisotropic scattering. 
    1193                                         testPhi = ran1a(&seed)*2*pi; 
    1194                                          
    1195                                         // is it on the detector?        
    1196                                         FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 
    1197                                                                                                          
    1198                                         if(xPixel != -1 && yPixel != -1) { 
    1199                                                 isOn += 1; 
    1200                                                 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 
    1201                                                 indices[0] = xPixel; 
    1202                                                 indices[1] = yPixel; 
    1203                                                 if (retVal = MDGetNumericWavePointValue(wavH, indices, value)) 
    1204                                                         return retVal; 
    1205                                                 value[0] += 1; // Real part 
    1206                                                 if (retVal = MDSetNumericWavePointValue(wavH, indices, value)) 
    1207                                                         return retVal; 
    1208                                                 //if(index==1)  // only the single scattering events 
    1209                                                         //dp = dp0 + xPixel + yPixel*numColumns;                //offset the pointer to the exact memory location 
    1210                                                         //*dp += 1;             //increment the value there 
    1211                                                 //endif 
    1212                                         } 
    1213                                          
    1214  
    1215         /*              is this causing me a problem since I'm not locking these? Probably not, since it crashes even if I comment these out... */ 
    1216                                         if(theta_z < theta_max) { 
    1217                                                 //Choose index for scattering angle array. 
    1218                                                 //IND = NINT(THETA_z/DTH + 0.4999999) 
    1219                                                 ind = (long)(theta_z/dth + 0.4999999);          //round is eqivalent to nint() 
    1220                                                 nt[ind] += 1;                   //Increment bin for angle. 
    1221                                                 //Increment angle array for single scattering events. 
    1222                                                 if (index == 1) { 
    1223                                                         j1[ind] += 1; 
    1224                                                 } 
    1225                                                 //Increment angle array for double scattering events. 
    1226                                                 if (index == 2) { 
    1227                                                         j2[ind] += 1; 
    1228                                                 } 
    1229                                         } 
    1230         /**/ 
    1231                                          
    1232                                         // increment all of the counters now since done==1 here and I'm sure to exit and get another neutron 
    1233                                         NScatterEvents += index;                //total number of scattering events 
    1234                                         if(index == 1 && incoherentEvent == 1) { 
    1235                                                 NSingleIncoherent += 1; 
    1236                                         } 
    1237                                         if(index == 1 && coherentEvent == 1) { 
    1238                                                 NSingleCoherent += 1; 
    1239                                         } 
    1240                                         if(index == 2 && coherentEvent == 1 && incoherentEvent == 0) { 
    1241                                                 NDoubleCoherent += 1; 
    1242                                         } 
    1243                                         if(index > 1) { 
    1244                                                 NMultipleScatter += 1; 
    1245                                         } 
    1246                                         if(coherentEvent >= 1 && incoherentEvent == 0) { 
    1247                                                 NCoherentEvents += 1; 
    1248                                         } 
    1249                                         if(coherentEvent > 1 && incoherentEvent == 0) { 
    1250                                                 NMultipleCoherent += 1; 
    1251                                         } 
    1252  
    1253                                 } else {        // index was zero, neutron must be transmitted, so just increment the proper counters and data 
    1254                                                 isOn += 1; 
    1255                                                 nt[0] += 1; 
    1256                                                 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 
    1257                                                 //indices[0] = xCtr_long;               //don't put everything in one pixel 
    1258                                                 //indices[1] = yCtr_long; 
    1259                                                 indices[0] = (long)(xCtr+xx/pixSize+0.5); 
    1260                                                 indices[1] = (long)(yCtr+yy/pixSize+0.5); 
    1261                                                 // check for valid indices - got an XOP error, probably from here 
    1262                                                 if(indices[0] > 127) indices[0] = 127; 
    1263                                                 if(indices[0] < 0) indices[0] = 0; 
    1264                                                 if(indices[1] > 127) indices[1] = 127; 
    1265                                                 if(indices[1] < 0) indices[1] = 0; 
    1266                                                  
    1267                                                 if (retVal = MDGetNumericWavePointValue(wavH, indices, value)) 
    1268                                                         return retVal; 
    1269                                                 value[0] += 1; // Real part 
    1270                                                 if (retVal = MDSetNumericWavePointValue(wavH, indices, value)) 
    1271                                                         return retVal; 
    1272                                         }        
    1273                         } 
    1274                  } while (!done); 
    1275         } while(n1 < imon); 
    1276          
    1277 // assign the results to the wave 
    1278  
    1279         MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 
    1280         value[0] = (double)n1; 
    1281         indices[0] = 0; 
    1282         if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 
    1283                 return retVal; 
    1284         value[0] = (double)n2; 
    1285         indices[0] = 1; 
    1286         if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 
    1287                 return retVal; 
    1288         value[0] = (double)isOn; 
    1289         indices[0] = 2; 
    1290         if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 
    1291                 return retVal; 
    1292         value[0] = (double)NScatterEvents; 
    1293         indices[0] = 3; 
    1294         if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 
    1295                 return retVal; 
    1296         value[0] = (double)NSingleCoherent; 
    1297         indices[0] = 4; 
    1298         if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 
    1299                 return retVal; 
    1300         value[0] = (double)NMultipleCoherent; 
    1301         indices[0] = 5; 
    1302         if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 
    1303                 return retVal; 
    1304         value[0] = (double)NMultipleScatter; 
    1305         indices[0] = 6; 
    1306         if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 
    1307                 return retVal;   
    1308         value[0] = (double)NCoherentEvents; 
    1309         indices[0] = 7; 
    1310         if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 
    1311                 return retVal; 
    1312  
    1313 //      HSetState((Handle)wavH, hState);                //release the handle of the 2D data wave 
    1314 //      WaveHandleModified(wavH);                       // Inform Igor that we have changed the wave. (CALLBACK! needed, but not allowed in Threading) 
    1315          
    1316         return(0); 
    1317 } 
    1318 ////////        end of X4 
    1319  
    1320  
    1321  
     468 
Note: See TracChangeset for help on using the changeset viewer.