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

    r789 r812  
    4949        long NMultipleCoherent,NCoherentEvents; 
    5050        double deltaLam,v1,v2,currWavelength,rsq,fac;           //for simulating wavelength distribution 
    51  
     51        double ssd, sourAp, souXX, souYY, magn;         //source-to-sample, and source Ap radius for initlal trajectory 
     52        double vz_1,g,yg_d;                             //gravity terms 
     53         
     54        vz_1 = 3.956e5;         //velocity [cm/s] of 1 A neutron 
     55        g = 981.0;                              //gravity acceleration [cm/s^2] 
     56         
    5257        // for accessing the 2D wave data, direct method (see the WaveAccess example XOP) 
    5358        waveHndl wavH; 
     
    136141        sig_sas = inputWave[10]; 
    137142        deltaLam = inputWave[11]; 
     143        ssd = inputWave[12];                    // in cm, like SDD 
     144        sourAp = inputWave[13];         // radius, in cm, like r1 and r2         
     145         
    138146        xCtr_long = (long)(xCtr+0.5);   //round() not defined in VC8 
    139147        yCtr_long = (long)(yCtr+0.5);   // and is probably wrong to use anyways (returns double!) 
     
    216224//              } 
    217225         
    218                 vx = 0.0;                       // Initialize direction vector. 
    219                 vy = 0.0; 
    220                 vz = 1.0; 
     226//              vx = 0.0;                       // Initialize direction vector. 
     227//              vy = 0.0; 
     228//              vz = 1.0; 
    221229                 
    222230                theta = 0.0;            //      Initialize scattering angle. 
     
    229237                coherentEvent = 0; 
    230238                 
     239                // pick point in source aperture area            
     240                do      {                               //      Makes sure position is within circle. 
     241                        ran = ran3(&seed);              //[0,1] 
     242                        souXX = 2.0*sourAp*(ran-0.5);           //X beam position of neutron entering sample. 
     243                        ran = ran3(&seed);              //[0,1] 
     244                        souYY = 2.0*sourAp*(ran-0.5);           //Y beam position ... 
     245                        rr = sqrt(souXX*souXX+souYY*souYY);             //Radial position of neutron in incident beam. 
     246                } while(rr>sourAp); 
     247                                 
     248                                 
     249                // pick point in sample aperture 
    231250                do      {                               //      Makes sure position is within circle. 
    232251                        ran = ran3(&seed);              //[0,1] 
     
    245264                        rsq=v1*v1+v2*v2; 
    246265                } while (rsq >= 1.0 || rsq == 0.0); 
    247                 fac=sqrt(-2.0*log(rsq)/rsq); 
     266                fac=sqrt(-2.0*log10(rsq)/rsq);          //be sure to use log10() here, to duplicate the Igor code 
    248267                                 
    249268                //              gset=v1*fac             //technically, I'm throwing away one of the two values 
    250269                currWavelength = (v2*fac)*deltaLam*wavelength/2.35 + wavelength; 
     270                 
     271                magn = sqrt((souXX - xx)*(souXX - xx) + (souYY - yy)*(souYY - yy) + ssd*ssd); 
     272                vx = (souXX - xx)/magn;         // Initialize direction vector. 
     273                vy = (souYY - yy)/magn; 
     274                vz = (ssd - 0.)/magn; 
     275                 
    251276                 
    252277                do {   //Scattering Loop, will exit when "done" == 1 
     
    329354                                //      nn[index] += 1; 
    330355                                } 
    331                                                                                                  
     356                                 
     357                                // calculate fall due to gravity (in cm) (note that it is negative) 
     358                                yg_d = -0.5*g*sdd*(ssd+sdd)*(currWavelength/vz_1)*(currWavelength/vz_1); 
     359                                 
     360                                 
    332361                                if( index != 0) {               //neutron was scattered, figure out where it went 
    333362                                        theta_z = acos(vz);             // Angle (= 2theta) WITH respect to z axis. 
     
    340369                                         
    341370                                        // is it on the detector?        
    342                                         FindPixel(testQ,testPhi,currWavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 
     371                                        FindPixel(testQ,testPhi,currWavelength,yg_d,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 
    343372                                                                                                         
    344373                                        if(xPixel != -1 && yPixel != -1) { 
     
    400429                                                isOn += 1; 
    401430                                                nt[0] += 1; 
    402                                                 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 
    403                                                 //indices[0] = xCtr_long;               //don't put everything in one pixel 
    404                                                 //indices[1] = yCtr_long; 
    405                                                 indices[0] = (long)(xCtr+xx/pixSize+0.5); 
    406                                                 indices[1] = (long)(yCtr+yy/pixSize+0.5); 
    407                                                 // check for valid indices - got an XOP error, probably from here 
    408                                                 if(indices[0] > 127) indices[0] = 127; 
    409                                                 if(indices[0] < 0) indices[0] = 0; 
    410                                                 if(indices[1] > 127) indices[1] = 127; 
    411                                                 if(indices[1] < 0) indices[1] = 0; 
    412                                                  
    413                                                 if (retVal = MDGetNumericWavePointValue(wavH, indices, value)) 
    414                                                         return retVal; 
    415                                                 value[0] += 1; // Real part 
    416                                                 if (retVal = MDSetNumericWavePointValue(wavH, indices, value)) 
    417                                                         return retVal; 
     431                                         
     432                                                //figure out where it landed 
     433                                                theta_z = acos(vz);             // Angle (= 2theta) WITH respect to z axis. 
     434                                                testQ = 2.0*pi*sin(theta_z)/currWavelength; 
     435                                                 
     436                                                // pick a random phi angle, and see if it lands on the detector 
     437                                                // since the scattering is isotropic, I can safely pick a new, random value 
     438                                                // this would not be true if simulating anisotropic scattering. 
     439                                                testPhi = ran3(&seed)*2.0*pi; 
     440                                                 
     441                                                // is it on the detector?        
     442                                                FindPixel(testQ,testPhi,currWavelength,yg_d,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 
     443                                                 
     444                                                if(xPixel != -1 && yPixel != -1) { 
     445                                                        isOn += 1; 
     446                                                        MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 
     447                                                        indices[0] = xPixel; 
     448                                                        indices[1] = yPixel; 
     449                                                        if (retVal = MDGetNumericWavePointValue(wavH, indices, value)) 
     450                                                                return retVal; 
     451                                                        value[0] += 1; // Real part 
     452                                                        if (retVal = MDSetNumericWavePointValue(wavH, indices, value)) 
     453                                                                return retVal; 
     454                                                        //if(index==1)  // only the single scattering events 
     455                                                        //dp = dp0 + xPixel + yPixel*numColumns;                //offset the pointer to the exact memory location 
     456                                                        //*dp += 1;             //increment the value there 
     457                                                        //endif 
     458                                                } 
     459                                                 
    418460                                        }        
    419461                        } 
Note: See TracChangeset for help on using the changeset viewer.