Ignore:
Timestamp:
Jun 21, 2011 1:09:46 PM (12 years ago)
Author:
srkline
Message:

Added source aperture size and SSD to simulation to get realistic initial trajectory, rather than only along the z-axis.

Added gravity fall to properly account for fall of neutrons and its effect on the primary transmitted beam and on the scattering pattern.

File:
1 edited

Legend:

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

    r711 r812  
    5252        long NMultipleCoherent,NCoherentEvents; 
    5353        double deltaLam,v1,v2,currWavelength,rsq,fac;           //for simulating wavelength distribution 
    54  
     54        double ssd, sourAp, souXX, souYY, magn;         //source-to-sample, and source Ap radius for initlal trajectory 
     55        double vz_1,g,yg_d;                             //gravity terms 
     56         
     57        vz_1 = 3.956e5;         //velocity [cm/s] of 1 A neutron 
     58        g = 981.0;                              //gravity acceleration [cm/s^2] 
     59         
    5560 
    5661        // for accessing the 2D wave data, direct method (see the WaveAccess example XOP) 
     
    140145        sig_sas = inputWave[10]; 
    141146        deltaLam = inputWave[11]; 
     147        ssd = inputWave[12];                    // in cm, like SDD 
     148        sourAp = inputWave[13];         // radius, in cm, like r1 and r2         
     149         
    142150        xCtr_long = (long)(xCtr+0.5); 
    143151        yCtr_long = (long)(yCtr+0.5); 
     
    191199        ratio = sig_incoh / sig_total; 
    192200         
    193         theta_max = wavelength*qmax/(2*pi); 
     201        theta_max = wavelength*qmax/(2.0*pi); 
    194202//C     SET Theta-STEP SIZE. 
    195203        dth = theta_max/num_bins; 
     
    220228//              } 
    221229         
    222                 vx = 0.0;                       // Initialize direction vector. 
    223                 vy = 0.0; 
    224                 vz = 1.0; 
     230//              vx = 0.0;                       // Initialize direction vector. 
     231//              vy = 0.0; 
     232//              vz = 1.0; 
    225233                 
    226234                theta = 0.0;            //      Initialize scattering angle. 
     
    233241                coherentEvent = 0; 
    234242                 
     243                 
     244                // pick point in source aperture area            
     245                do      {                               //      Makes sure position is within circle. 
     246                        ran = ran1(&seed);              //[0,1] 
     247                        souXX = 2.0*sourAp*(ran-0.5);           //X beam position of neutron entering sample. 
     248                        ran = ran1(&seed);              //[0,1] 
     249                        souYY = 2.0*sourAp*(ran-0.5);           //Y beam position ... 
     250                        rr = sqrt(souXX*souXX+souYY*souYY);             //Radial position of neutron in incident beam. 
     251                } while(rr>sourAp); 
     252                 
     253                // pick point in sample aperture 
    235254                do      {                               //      Makes sure position is within circle. 
    236255                        ran = ran1(&seed);              //[0,1] 
     
    249268                        rsq=v1*v1+v2*v2; 
    250269                } while (rsq >= 1.0 || rsq == 0.0); 
    251                 fac=sqrt(-2.0*log(rsq)/rsq); 
     270                fac=sqrt(-2.0*log10(rsq)/rsq);          //be sure to use log10() here, to duplicate the Igor code 
    252271                 
    253272                //              gset=v1*fac             //technically, I'm throwing away one of the two values 
    254273                currWavelength = (v2*fac)*deltaLam*wavelength/2.35 + wavelength; 
     274 
     275                magn = sqrt((souXX - xx)*(souXX - xx) + (souYY - yy)*(souYY - yy) + ssd*ssd); 
     276                vx = (souXX - xx)/magn;         // Initialize direction vector. 
     277                vy = (souYY - yy)/magn; 
     278                vz = (ssd - 0.)/magn; 
    255279                 
    256280                do {   //Scattering Loop, will exit when "done" == 1 
     
    294318                                                 
    295319                                                q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran1(&seed))*delta; 
    296                                                 theta = q0/2/pi*currWavelength;         //SAS approximation 
     320                                                theta = q0/2.0/pi*currWavelength;               //SAS approximation 
    297321                                                 
    298322                                                find_theta = 1;         //always accept 
     
    333357                                //      nn[index] += 1; 
    334358                                } 
    335                                                                                                  
     359                                 
     360                                // calculate fall due to gravity (in cm) (note that it is negative) 
     361                                yg_d = -0.5*g*sdd*(ssd+sdd)*(currWavelength/vz_1)*(currWavelength/vz_1); 
     362                                 
    336363                                if( index != 0) {               //neutron was scattered, figure out where it went 
    337364                                        theta_z = acos(vz);             // Angle (= 2theta) WITH respect to z axis. 
    338                                         testQ = 2*pi*sin(theta_z)/currWavelength; 
     365                                        testQ = 2.0*pi*sin(theta_z)/currWavelength; 
    339366                                         
    340367                                        // pick a random phi angle, and see if it lands on the detector 
    341368                                        // since the scattering is isotropic, I can safely pick a new, random value 
    342369                                        // this would not be true if simulating anisotropic scattering. 
    343                                         testPhi = ran1(&seed)*2*pi; 
     370                                        testPhi = ran1(&seed)*2.0*pi; 
    344371                                         
    345372                                        // is it on the detector?        
    346                                         FindPixel(testQ,testPhi,currWavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 
     373                                        FindPixel(testQ,testPhi,currWavelength,yg_d,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 
    347374                                                                                                         
    348375                                        if(xPixel != -1 && yPixel != -1) { 
     
    404431                                                isOn += 1; 
    405432                                                nt[0] += 1; 
    406                                                 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 
    407                                                 //indices[0] = xCtr_long;               //don't put everything in one pixel 
    408                                                 //indices[1] = yCtr_long; 
    409                                                 indices[0] = (long)(xCtr+xx/pixSize+0.5); 
    410                                                 indices[1] = (long)(yCtr+yy/pixSize+0.5); 
    411                                                 // check for valid indices - got an XOP error, probably from here 
    412                                                 if(indices[0] > 127) indices[0] = 127; 
    413                                                 if(indices[0] < 0) indices[0] = 0; 
    414                                                 if(indices[1] > 127) indices[1] = 127; 
    415                                                 if(indices[1] < 0) indices[1] = 0; 
    416                                                  
    417                                                 if (retVal = MDGetNumericWavePointValue(wavH, indices, value)) 
    418                                                         return retVal; 
    419                                                 value[0] += 1; // Real part 
    420                                                 if (retVal = MDSetNumericWavePointValue(wavH, indices, value)) 
    421                                                         return retVal; 
     433                                                 
     434                                                //figure out where it landed 
     435                                                theta_z = acos(vz);             // Angle (= 2theta) WITH respect to z axis. 
     436                                                testQ = 2.0*pi*sin(theta_z)/currWavelength; 
     437                                                 
     438                                                // pick a random phi angle, and see if it lands on the detector 
     439                                                // since the scattering is isotropic, I can safely pick a new, random value 
     440                                                // this would not be true if simulating anisotropic scattering. 
     441                                                testPhi = ran1(&seed)*2.0*pi; 
     442                                                 
     443                                                // is it on the detector?        
     444                                                FindPixel(testQ,testPhi,currWavelength,yg_d,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 
     445                                                 
     446                                                if(xPixel != -1 && yPixel != -1) { 
     447                                                        isOn += 1; 
     448                                                        MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 
     449                                                        indices[0] = xPixel; 
     450                                                        indices[1] = yPixel; 
     451                                                        if (retVal = MDGetNumericWavePointValue(wavH, indices, value)) 
     452                                                                return retVal; 
     453                                                        value[0] += 1; // Real part 
     454                                                        if (retVal = MDSetNumericWavePointValue(wavH, indices, value)) 
     455                                                                return retVal; 
     456                                                        //if(index==1)  // only the single scattering events 
     457                                                        //dp = dp0 + xPixel + yPixel*numColumns;                //offset the pointer to the exact memory location 
     458                                                        //*dp += 1;             //increment the value there 
     459                                                        //endif 
     460                                                }                                        
     461                                         
    422462                                        }        
    423463                        } 
Note: See TracChangeset for help on using the changeset viewer.