Changeset 812 for sans/XOP_Dev


Ignore:
Timestamp:
Jun 21, 2011 1:09:46 PM (11 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.

Location:
sans/XOP_Dev/MonteCarlo
Files:
6 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                        } 
  • sans/XOP_Dev/MonteCarlo/MonteCarlo.h

    r758 r812  
    5656 
    5757long MC_round(double x); 
    58 int FindPixel(double testQ, double testPhi, double lam, double sdd, double pixSize, double xCtr, double yCtr, long *xPixel, long *yPixel); 
     58int FindPixel(double testQ, double testPhi, double lam, double yg_d, double sdd, double pixSize, double xCtr, double yCtr, long *xPixel, long *yPixel); 
    5959int NewDirection(double *vx, double *vy, double *vz, double theta, double phi); 
    6060double path_len(double aval, double sig_tot); 
  • 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                        } 
  • sans/XOP_Dev/MonteCarlo/MonteCarlo3.c

    r711 r812  
    5151        long NMultipleCoherent,NCoherentEvents; 
    5252        double deltaLam,v1,v2,currWavelength,rsq,fac;           //for simulating wavelength distribution 
    53          
     53        double ssd, sourAp, souXX, souYY, magn;         //source-to-sample, and source Ap radius for initlal trajectory 
     54        double vz_1,g,yg_d;                             //gravity terms 
     55         
     56        vz_1 = 3.956e5;         //velocity [cm/s] of 1 A neutron 
     57        g = 981.0;                              //gravity acceleration [cm/s^2]  
    5458         
    5559        // for accessing the 2D wave data, direct method (see the WaveAccess example XOP) 
     
    139143        sig_sas = inputWave[10]; 
    140144        deltaLam = inputWave[11]; 
     145        ssd = inputWave[12];                    // in cm, like SDD 
     146        sourAp = inputWave[13];         // radius, in cm, like r1 and r2         
     147         
    141148        xCtr_long = (long)(xCtr+0.5); 
    142149        yCtr_long = (long)(yCtr+0.5); 
     
    190197        ratio = sig_incoh / sig_total; 
    191198         
    192         theta_max = wavelength*qmax/(2*pi); 
     199        theta_max = wavelength*qmax/(2.0*pi); 
    193200        //C     SET Theta-STEP SIZE. 
    194201        dth = theta_max/num_bins; 
     
    219226                //              } 
    220227                 
    221                 vx = 0.0;                       // Initialize direction vector. 
    222                 vy = 0.0; 
    223                 vz = 1.0; 
     228//              vx = 0.0;                       // Initialize direction vector. 
     229//              vy = 0.0; 
     230//              vz = 1.0; 
    224231                 
    225232                theta = 0.0;            //      Initialize scattering angle. 
     
    231238                incoherentEvent = 0; 
    232239                coherentEvent = 0; 
    233                  
     240 
     241                // pick point in source aperture area            
     242                do      {                               //      Makes sure position is within circle. 
     243                        ran = ran3a(&seed);             //[0,1] 
     244                        souXX = 2.0*sourAp*(ran-0.5);           //X beam position of neutron entering sample. 
     245                        ran = ran3a(&seed);             //[0,1] 
     246                        souYY = 2.0*sourAp*(ran-0.5);           //Y beam position ... 
     247                        rr = sqrt(souXX*souXX+souYY*souYY);             //Radial position of neutron in incident beam. 
     248                } while(rr>sourAp); 
     249 
     250                // pick point in sample aperture 
    234251                do      {                               //      Makes sure position is within circle. 
    235252                        ran = ran3a(&seed);             //[0,1] 
     
    248265                        rsq=v1*v1+v2*v2; 
    249266                } while (rsq >= 1.0 || rsq == 0.0); 
    250                 fac=sqrt(-2.0*log(rsq)/rsq); 
     267                fac=sqrt(-2.0*log10(rsq)/rsq);          //be sure to use log10() here, to duplicate the Igor code 
    251268                 
    252269                //              gset=v1*fac             //technically, I'm throwing away one of the two values 
    253270                currWavelength = (v2*fac)*deltaLam*wavelength/2.35 + wavelength; 
     271                 
     272                magn = sqrt((souXX - xx)*(souXX - xx) + (souYY - yy)*(souYY - yy) + ssd*ssd); 
     273                vx = (souXX - xx)/magn;         // Initialize direction vector. 
     274                vy = (souYY - yy)/magn; 
     275                vz = (ssd - 0.)/magn; 
    254276                 
    255277                do {   //Scattering Loop, will exit when "done" == 1 
     
    293315                                                 
    294316                                                q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran3a(&seed))*delta; 
    295                                                 theta = q0/2/pi*currWavelength;         //SAS approximation 
     317                                                theta = q0/2.0/pi*currWavelength;               //SAS approximation 
    296318                                                 
    297319                                                find_theta = 1;         //always accept 
     
    333355                                } 
    334356                                 
     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                                 
    335360                                if( index != 0) {               //neutron was scattered, figure out where it went 
    336361                                        theta_z = acos(vz);             // Angle (= 2theta) WITH respect to z axis. 
    337                                         testQ = 2*pi*sin(theta_z)/currWavelength; 
     362                                        testQ = 2.0*pi*sin(theta_z)/currWavelength; 
    338363                                         
    339364                                        // pick a random phi angle, and see if it lands on the detector 
    340365                                        // since the scattering is isotropic, I can safely pick a new, random value 
    341366                                        // this would not be true if simulating anisotropic scattering. 
    342                                         testPhi = ran3a(&seed)*2*pi; 
     367                                        testPhi = ran3a(&seed)*2.0*pi; 
    343368                                         
    344369                                        // is it on the detector?        
    345                                         FindPixel(testQ,testPhi,currWavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 
     370                                        FindPixel(testQ,testPhi,currWavelength,yg_d,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 
    346371                                         
    347372                                        if(xPixel != -1 && yPixel != -1) { 
     
    403428                                        isOn += 1; 
    404429                                        nt[0] += 1; 
    405                                         MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 
    406                                         //indices[0] = xCtr_long;               //don't put everything in one pixel 
    407                                         //indices[1] = yCtr_long; 
    408                                         indices[0] = (long)(xCtr+xx/pixSize+0.5); 
    409                                         indices[1] = (long)(yCtr+yy/pixSize+0.5); 
    410                                         // check for valid indices - got an XOP error, probably from here 
    411                                         if(indices[0] > 127) indices[0] = 127; 
    412                                         if(indices[0] < 0) indices[0] = 0; 
    413                                         if(indices[1] > 127) indices[1] = 127; 
    414                                         if(indices[1] < 0) indices[1] = 0; 
    415                                          
    416                                         if (retVal = MDGetNumericWavePointValue(wavH, indices, value)) 
    417                                                 return retVal; 
    418                                         value[0] += 1; // Real part 
    419                                         if (retVal = MDSetNumericWavePointValue(wavH, indices, value)) 
    420                                                 return retVal; 
     430                                         
     431                                        //figure out where it landed 
     432                                        theta_z = acos(vz);             // Angle (= 2theta) WITH respect to z axis. 
     433                                        testQ = 2.0*pi*sin(theta_z)/currWavelength; 
     434                                         
     435                                        // pick a random phi angle, and see if it lands on the detector 
     436                                        // since the scattering is isotropic, I can safely pick a new, random value 
     437                                        // this would not be true if simulating anisotropic scattering. 
     438                                        testPhi = ran3a(&seed)*2.0*pi; 
     439                                         
     440                                        // is it on the detector?        
     441                                        FindPixel(testQ,testPhi,currWavelength,yg_d,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 
     442                                         
     443                                        if(xPixel != -1 && yPixel != -1) { 
     444                                                isOn += 1; 
     445                                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 
     446                                                indices[0] = xPixel; 
     447                                                indices[1] = yPixel; 
     448                                                if (retVal = MDGetNumericWavePointValue(wavH, indices, value)) 
     449                                                        return retVal; 
     450                                                value[0] += 1; // Real part 
     451                                                if (retVal = MDSetNumericWavePointValue(wavH, indices, value)) 
     452                                                        return retVal; 
     453                                                //if(index==1)  // only the single scattering events 
     454                                                //dp = dp0 + xPixel + yPixel*numColumns;                //offset the pointer to the exact memory location 
     455                                                //*dp += 1;             //increment the value there 
     456                                                //endif 
     457                                        } 
    421458                                }        
    422459                        } 
  • sans/XOP_Dev/MonteCarlo/MonteCarlo4.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]  
    5559         
    5660        // for accessing the 2D wave data, direct method (see the WaveAccess example XOP) 
     
    140144        sig_sas = inputWave[10]; 
    141145        deltaLam = inputWave[11]; 
     146        ssd = inputWave[12];                    // in cm, like SDD 
     147        sourAp = inputWave[13];         // radius, in cm, like r1 and r2 
     148         
    142149        xCtr_long = (long)(xCtr+0.5); 
    143150        yCtr_long = (long)(yCtr+0.5); 
     
    191198        ratio = sig_incoh / sig_total; 
    192199         
    193         theta_max = wavelength*qmax/(2*pi); 
     200        theta_max = wavelength*qmax/(2.0*pi); 
    194201        //C     SET Theta-STEP SIZE. 
    195202        dth = theta_max/num_bins; 
     
    220227                //              } 
    221228                 
    222                 vx = 0.0;                       // Initialize direction vector. 
    223                 vy = 0.0; 
    224                 vz = 1.0; 
     229//              vx = 0.0;                       // Initialize direction vector. 
     230//              vy = 0.0; 
     231//              vz = 1.0; 
    225232                 
    226233                theta = 0.0;            //      Initialize scattering angle. 
     
    233240                coherentEvent = 0; 
    234241                 
     242                // pick point in source aperture area            
     243                do      {                               //      Makes sure position is within circle. 
     244                        ran = ran1a(&seed);             //[0,1] 
     245                        souXX = 2.0*sourAp*(ran-0.5);           //X beam position of neutron entering sample. 
     246                        ran = ran1a(&seed);             //[0,1] 
     247                        souYY = 2.0*sourAp*(ran-0.5);           //Y beam position ... 
     248                        rr = sqrt(souXX*souXX+souYY*souYY);             //Radial position of neutron in incident beam. 
     249                } while(rr>sourAp); 
     250                 
     251                 
     252                // pick point in sample aperture                 
    235253                do      {                               //      Makes sure position is within circle. 
    236254                        ran = ran1a(&seed);             //[0,1] 
     
    249267                        rsq=v1*v1+v2*v2; 
    250268                } while (rsq >= 1.0 || rsq == 0.0); 
    251                 fac=sqrt(-2.0*log(rsq)/rsq); 
     269                fac=sqrt(-2.0*log10(rsq)/rsq);          //be sure to use log10() here, to duplicate the Igor code 
    252270                 
    253271                //              gset=v1*fac             //technically, I'm throwing away one of the two values 
    254272                currWavelength = (v2*fac)*deltaLam*wavelength/2.35 + wavelength; 
     273                 
     274                magn = sqrt((souXX - xx)*(souXX - xx) + (souYY - yy)*(souYY - yy) + ssd*ssd); 
     275                vx = (souXX - xx)/magn;         // Initialize direction vector. 
     276                vy = (souYY - yy)/magn; 
     277                vz = (ssd - 0.)/magn; 
    255278                 
    256279                do {   //Scattering Loop, will exit when "done" == 1 
     
    294317                                                 
    295318                                                q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran1a(&seed))*delta; 
    296                                                 theta = q0/2/pi*currWavelength;         //SAS approximation 
     319                                                theta = q0/2.0/pi*currWavelength;               //SAS approximation 
    297320                                                 
    298321                                                find_theta = 1;         //always accept 
     
    333356                                        //      nn[index] += 1; 
    334357                                } 
     358 
     359                                // calculate fall due to gravity (in cm) (note that it is negative) 
     360                                yg_d = -0.5*g*sdd*(ssd+sdd)*(currWavelength/vz_1)*(currWavelength/vz_1); 
    335361                                 
    336362                                if( index != 0) {               //neutron was scattered, figure out where it went 
    337363                                        theta_z = acos(vz);             // Angle (= 2theta) WITH respect to z axis. 
    338                                         testQ = 2*pi*sin(theta_z)/currWavelength; 
     364                                        testQ = 2.0*pi*sin(theta_z)/currWavelength; 
    339365                                         
    340366                                        // pick a random phi angle, and see if it lands on the detector 
    341367                                        // since the scattering is isotropic, I can safely pick a new, random value 
    342368                                        // this would not be true if simulating anisotropic scattering. 
    343                                         testPhi = ran1a(&seed)*2*pi; 
     369                                        testPhi = ran1a(&seed)*2.0*pi; 
    344370                                         
    345371                                        // is it on the detector?        
    346                                         FindPixel(testQ,testPhi,currWavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 
     372                                        FindPixel(testQ,testPhi,currWavelength,yg_d,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 
    347373                                         
    348374                                        if(xPixel != -1 && yPixel != -1) { 
     
    404430                                        isOn += 1; 
    405431                                        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; 
     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 = ran1a(&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                                        } 
    422459                                }        
    423460                        } 
  • sans/XOP_Dev/MonteCarlo/MonteCarlo_Main.c

    r789 r812  
    3434 
    3535int 
    36 FindPixel(double testQ, double testPhi, double lam, double sdd, 
     36FindPixel(double testQ, double testPhi, double lam, double yg_d, double sdd, 
    3737                  double pixSize, double xCtr, double yCtr, long *xPixel, long *yPixel) { 
    3838         
     
    4242        qx = testQ*cos(testPhi); 
    4343        qy = testQ*sin(testPhi); 
     44         
     45        //correct qy for gravity 
     46        qy += 4.0*pi/lam*(yg_d/sdd/2.0); 
     47 
    4448         
    4549        //convert qx,qy to pixel locations relative to # of pixels x, y from center 
Note: See TracChangeset for help on using the changeset viewer.