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/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.