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.

File:
1 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 } 
Note: See TracChangeset for help on using the changeset viewer.