Changeset 812 for sans/XOP_Dev/MonteCarlo/MonteCarlo2.c
- Timestamp:
- Jun 21, 2011 1:09:46 PM (12 years ago)
- File:
-
- 1 edited
-
sans/XOP_Dev/MonteCarlo/MonteCarlo2.c (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
sans/XOP_Dev/MonteCarlo/MonteCarlo2.c
r711 r812 52 52 long NMultipleCoherent,NCoherentEvents; 53 53 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 55 60 56 61 // for accessing the 2D wave data, direct method (see the WaveAccess example XOP) … … 140 145 sig_sas = inputWave[10]; 141 146 deltaLam = inputWave[11]; 147 ssd = inputWave[12]; // in cm, like SDD 148 sourAp = inputWave[13]; // radius, in cm, like r1 and r2 149 142 150 xCtr_long = (long)(xCtr+0.5); 143 151 yCtr_long = (long)(yCtr+0.5); … … 191 199 ratio = sig_incoh / sig_total; 192 200 193 theta_max = wavelength*qmax/(2 *pi);201 theta_max = wavelength*qmax/(2.0*pi); 194 202 //C SET Theta-STEP SIZE. 195 203 dth = theta_max/num_bins; … … 220 228 // } 221 229 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; 225 233 226 234 theta = 0.0; // Initialize scattering angle. … … 233 241 coherentEvent = 0; 234 242 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 235 254 do { // Makes sure position is within circle. 236 255 ran = ran1(&seed); //[0,1] … … 249 268 rsq=v1*v1+v2*v2; 250 269 } 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 252 271 253 272 // gset=v1*fac //technically, I'm throwing away one of the two values 254 273 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; 255 279 256 280 do { //Scattering Loop, will exit when "done" == 1 … … 294 318 295 319 q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran1(&seed))*delta; 296 theta = q0/2 /pi*currWavelength; //SAS approximation320 theta = q0/2.0/pi*currWavelength; //SAS approximation 297 321 298 322 find_theta = 1; //always accept … … 333 357 // nn[index] += 1; 334 358 } 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 336 363 if( index != 0) { //neutron was scattered, figure out where it went 337 364 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; 339 366 340 367 // pick a random phi angle, and see if it lands on the detector 341 368 // since the scattering is isotropic, I can safely pick a new, random value 342 369 // this would not be true if simulating anisotropic scattering. 343 testPhi = ran1(&seed)*2 *pi;370 testPhi = ran1(&seed)*2.0*pi; 344 371 345 372 // 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); 347 374 348 375 if(xPixel != -1 && yPixel != -1) { … … 404 431 isOn += 1; 405 432 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 422 462 } 423 463 }
Note: See TracChangeset
for help on using the changeset viewer.
