Ignore:
Timestamp:
Jan 28, 2010 6:09:34 PM (13 years ago)
Author:
srkline
Message:

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

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

File:
1 edited

Legend:

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

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