Changeset 812 for sans/XOP_Dev/MonteCarlo/MonteCarlo.c

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:

• sans/XOP_Dev/MonteCarlo/MonteCarlo.c (modified) (8 diffs)

sans/XOP_Dev/MonteCarlo/MonteCarlo.c

@@ -49 +49 @@
         long NMultipleCoherent,NCoherentEvents;
         double deltaLam,v1,v2,currWavelength,rsq,fac;           //for simulating wavelength distribution
-
+        double ssd, sourAp, souXX, souYY, magn;         //source-to-sample, and source Ap radius for initlal trajectory
+        double vz_1,g,yg_d;                             //gravity terms
+        
+        vz_1 = 3.956e5;         //velocity [cm/s] of 1 A neutron
+        g = 981.0;                              //gravity acceleration [cm/s^2]
+        
         // for accessing the 2D wave data, direct method (see the WaveAccess example XOP)
         waveHndl wavH;
@@ -136 +141 @@
         sig_sas = inputWave[10];
         deltaLam = inputWave[11];
+        ssd = inputWave[12];                    // in cm, like SDD
+        sourAp = inputWave[13];         // radius, in cm, like r1 and r2        
+        
         xCtr_long = (long)(xCtr+0.5);   //round() not defined in VC8
         yCtr_long = (long)(yCtr+0.5);   // and is probably wrong to use anyways (returns double!)
@@ -216 +224 @@
 //              }
         
-                vx = 0.0;                       // Initialize direction vector.
-                vy = 0.0;
-                vz = 1.0;
+//              vx = 0.0;                       // Initialize direction vector.
+//              vy = 0.0;
+//              vz = 1.0;
                 
                 theta = 0.0;            //      Initialize scattering angle.
@@ -229 +237 @@
                 coherentEvent = 0;
                 
+                // pick point in source aperture area           
+                do      {                               //      Makes sure position is within circle.
+                        ran = ran3(&seed);              //[0,1]
+                        souXX = 2.0*sourAp*(ran-0.5);           //X beam position of neutron entering sample.
+                        ran = ran3(&seed);              //[0,1]
+                        souYY = 2.0*sourAp*(ran-0.5);           //Y beam position ...
+                        rr = sqrt(souXX*souXX+souYY*souYY);             //Radial position of neutron in incident beam.
+                } while(rr>sourAp);
+                                
+                                
+                // pick point in sample aperture
                 do      {                               //      Makes sure position is within circle.
                         ran = ran3(&seed);              //[0,1]
@@ -245 +264 @@
                         rsq=v1*v1+v2*v2;
                 } while (rsq >= 1.0 || rsq == 0.0);
-                fac=sqrt(-2.0*log(rsq)/rsq);
+                fac=sqrt(-2.0*log10(rsq)/rsq);          //be sure to use log10() here, to duplicate the Igor code
                                 
                 //              gset=v1*fac             //technically, I'm throwing away one of the two values
                 currWavelength = (v2*fac)*deltaLam*wavelength/2.35 + wavelength;
+                
+                magn = sqrt((souXX - xx)*(souXX - xx) + (souYY - yy)*(souYY - yy) + ssd*ssd);
+                vx = (souXX - xx)/magn;         // Initialize direction vector.
+                vy = (souYY - yy)/magn;
+                vz = (ssd - 0.)/magn;
+                
                 
                 do {   //Scattering Loop, will exit when "done" == 1
@@ -329 +354 @@
                                 //      nn[index] += 1;
                                 }
-                                                                                                
+                                
+                                // calculate fall due to gravity (in cm) (note that it is negative)
+                                yg_d = -0.5*g*sdd*(ssd+sdd)*(currWavelength/vz_1)*(currWavelength/vz_1);
+                                
+                                
                                 if( index != 0) {               //neutron was scattered, figure out where it went
                                         theta_z = acos(vz);             // Angle (= 2theta) WITH respect to z axis.
@@ -340 +369 @@
                                         
                                         // is it on the detector?       
-                                        FindPixel(testQ,testPhi,currWavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel);
+                                        FindPixel(testQ,testPhi,currWavelength,yg_d,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel);
                                                                                                         
                                         if(xPixel != -1 && yPixel != -1) {
@@ -400 +429 @@
                                                 isOn += 1;
                                                 nt[0] += 1;
-                                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
-                                                //indices[0] = xCtr_long;               //don't put everything in one pixel
-                                                //indices[1] = yCtr_long;
-                                                indices[0] = (long)(xCtr+xx/pixSize+0.5);
-                                                indices[1] = (long)(yCtr+yy/pixSize+0.5);
-                                                // check for valid indices - got an XOP error, probably from here
-                                                if(indices[0] > 127) indices[0] = 127;
-                                                if(indices[0] < 0) indices[0] = 0;
-                                                if(indices[1] > 127) indices[1] = 127;
-                                                if(indices[1] < 0) indices[1] = 0;
-                                                
-                                                if (retVal = MDGetNumericWavePointValue(wavH, indices, value))
-                                                        return retVal;
-                                                value[0] += 1; // Real part
-                                                if (retVal = MDSetNumericWavePointValue(wavH, indices, value))
-                                                        return retVal;
+                                        
+                                                //figure out where it landed
+                                                theta_z = acos(vz);             // Angle (= 2theta) WITH respect to z axis.
+                                                testQ = 2.0*pi*sin(theta_z)/currWavelength;
+                                                
+                                                // pick a random phi angle, and see if it lands on the detector
+                                                // since the scattering is isotropic, I can safely pick a new, random value
+                                                // this would not be true if simulating anisotropic scattering.
+                                                testPhi = ran3(&seed)*2.0*pi;
+                                                
+                                                // is it on the detector?       
+                                                FindPixel(testQ,testPhi,currWavelength,yg_d,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel);
+                                                
+                                                if(xPixel != -1 && yPixel != -1) {
+                                                        isOn += 1;
+                                                        MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
+                                                        indices[0] = xPixel;
+                                                        indices[1] = yPixel;
+                                                        if (retVal = MDGetNumericWavePointValue(wavH, indices, value))
+                                                                return retVal;
+                                                        value[0] += 1; // Real part
+                                                        if (retVal = MDSetNumericWavePointValue(wavH, indices, value))
+                                                                return retVal;
+                                                        //if(index==1)  // only the single scattering events
+                                                        //dp = dp0 + xPixel + yPixel*numColumns;                //offset the pointer to the exact memory location
+                                                        //*dp += 1;             //increment the value there
+                                                        //endif
+                                                }
+                                                
                                         }       
                         }