Changeset 812 for sans/XOP_Dev/MonteCarlo
- Timestamp:
- Jun 21, 2011 1:09:46 PM (12 years ago)
- Location:
- sans/XOP_Dev/MonteCarlo
- Files:
-
- 6 edited
-
MonteCarlo.c (modified) (8 diffs)
-
MonteCarlo.h (modified) (1 diff)
-
MonteCarlo2.c (modified) (9 diffs)
-
MonteCarlo3.c (modified) (9 diffs)
-
MonteCarlo4.c (modified) (9 diffs)
-
MonteCarlo_Main.c (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
sans/XOP_Dev/MonteCarlo/MonteCarlo.c
r789 r812 49 49 long NMultipleCoherent,NCoherentEvents; 50 50 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 52 57 // for accessing the 2D wave data, direct method (see the WaveAccess example XOP) 53 58 waveHndl wavH; … … 136 141 sig_sas = inputWave[10]; 137 142 deltaLam = inputWave[11]; 143 ssd = inputWave[12]; // in cm, like SDD 144 sourAp = inputWave[13]; // radius, in cm, like r1 and r2 145 138 146 xCtr_long = (long)(xCtr+0.5); //round() not defined in VC8 139 147 yCtr_long = (long)(yCtr+0.5); // and is probably wrong to use anyways (returns double!) … … 216 224 // } 217 225 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; 221 229 222 230 theta = 0.0; // Initialize scattering angle. … … 229 237 coherentEvent = 0; 230 238 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 231 250 do { // Makes sure position is within circle. 232 251 ran = ran3(&seed); //[0,1] … … 245 264 rsq=v1*v1+v2*v2; 246 265 } 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 248 267 249 268 // gset=v1*fac //technically, I'm throwing away one of the two values 250 269 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 251 276 252 277 do { //Scattering Loop, will exit when "done" == 1 … … 329 354 // nn[index] += 1; 330 355 } 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 332 361 if( index != 0) { //neutron was scattered, figure out where it went 333 362 theta_z = acos(vz); // Angle (= 2theta) WITH respect to z axis. … … 340 369 341 370 // 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); 343 372 344 373 if(xPixel != -1 && yPixel != -1) { … … 400 429 isOn += 1; 401 430 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 418 460 } 419 461 } -
sans/XOP_Dev/MonteCarlo/MonteCarlo.h
r758 r812 56 56 57 57 long 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);58 int FindPixel(double testQ, double testPhi, double lam, double yg_d, double sdd, double pixSize, double xCtr, double yCtr, long *xPixel, long *yPixel); 59 59 int NewDirection(double *vx, double *vy, double *vz, double theta, double phi); 60 60 double path_len(double aval, double sig_tot); -
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 } -
sans/XOP_Dev/MonteCarlo/MonteCarlo3.c
r711 r812 51 51 long NMultipleCoherent,NCoherentEvents; 52 52 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] 54 58 55 59 // for accessing the 2D wave data, direct method (see the WaveAccess example XOP) … … 139 143 sig_sas = inputWave[10]; 140 144 deltaLam = inputWave[11]; 145 ssd = inputWave[12]; // in cm, like SDD 146 sourAp = inputWave[13]; // radius, in cm, like r1 and r2 147 141 148 xCtr_long = (long)(xCtr+0.5); 142 149 yCtr_long = (long)(yCtr+0.5); … … 190 197 ratio = sig_incoh / sig_total; 191 198 192 theta_max = wavelength*qmax/(2 *pi);199 theta_max = wavelength*qmax/(2.0*pi); 193 200 //C SET Theta-STEP SIZE. 194 201 dth = theta_max/num_bins; … … 219 226 // } 220 227 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; 224 231 225 232 theta = 0.0; // Initialize scattering angle. … … 231 238 incoherentEvent = 0; 232 239 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 234 251 do { // Makes sure position is within circle. 235 252 ran = ran3a(&seed); //[0,1] … … 248 265 rsq=v1*v1+v2*v2; 249 266 } 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 251 268 252 269 // gset=v1*fac //technically, I'm throwing away one of the two values 253 270 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; 254 276 255 277 do { //Scattering Loop, will exit when "done" == 1 … … 293 315 294 316 q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran3a(&seed))*delta; 295 theta = q0/2 /pi*currWavelength; //SAS approximation317 theta = q0/2.0/pi*currWavelength; //SAS approximation 296 318 297 319 find_theta = 1; //always accept … … 333 355 } 334 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 335 360 if( index != 0) { //neutron was scattered, figure out where it went 336 361 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; 338 363 339 364 // pick a random phi angle, and see if it lands on the detector 340 365 // since the scattering is isotropic, I can safely pick a new, random value 341 366 // this would not be true if simulating anisotropic scattering. 342 testPhi = ran3a(&seed)*2 *pi;367 testPhi = ran3a(&seed)*2.0*pi; 343 368 344 369 // 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); 346 371 347 372 if(xPixel != -1 && yPixel != -1) { … … 403 428 isOn += 1; 404 429 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 } 421 458 } 422 459 } -
sans/XOP_Dev/MonteCarlo/MonteCarlo4.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] 55 59 56 60 // for accessing the 2D wave data, direct method (see the WaveAccess example XOP) … … 140 144 sig_sas = inputWave[10]; 141 145 deltaLam = inputWave[11]; 146 ssd = inputWave[12]; // in cm, like SDD 147 sourAp = inputWave[13]; // radius, in cm, like r1 and r2 148 142 149 xCtr_long = (long)(xCtr+0.5); 143 150 yCtr_long = (long)(yCtr+0.5); … … 191 198 ratio = sig_incoh / sig_total; 192 199 193 theta_max = wavelength*qmax/(2 *pi);200 theta_max = wavelength*qmax/(2.0*pi); 194 201 //C SET Theta-STEP SIZE. 195 202 dth = theta_max/num_bins; … … 220 227 // } 221 228 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; 225 232 226 233 theta = 0.0; // Initialize scattering angle. … … 233 240 coherentEvent = 0; 234 241 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 235 253 do { // Makes sure position is within circle. 236 254 ran = ran1a(&seed); //[0,1] … … 249 267 rsq=v1*v1+v2*v2; 250 268 } 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 252 270 253 271 // gset=v1*fac //technically, I'm throwing away one of the two values 254 272 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; 255 278 256 279 do { //Scattering Loop, will exit when "done" == 1 … … 294 317 295 318 q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran1a(&seed))*delta; 296 theta = q0/2 /pi*currWavelength; //SAS approximation319 theta = q0/2.0/pi*currWavelength; //SAS approximation 297 320 298 321 find_theta = 1; //always accept … … 333 356 // nn[index] += 1; 334 357 } 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); 335 361 336 362 if( index != 0) { //neutron was scattered, figure out where it went 337 363 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; 339 365 340 366 // pick a random phi angle, and see if it lands on the detector 341 367 // since the scattering is isotropic, I can safely pick a new, random value 342 368 // this would not be true if simulating anisotropic scattering. 343 testPhi = ran1a(&seed)*2 *pi;369 testPhi = ran1a(&seed)*2.0*pi; 344 370 345 371 // 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); 347 373 348 374 if(xPixel != -1 && yPixel != -1) { … … 404 430 isOn += 1; 405 431 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 } 422 459 } 423 460 } -
sans/XOP_Dev/MonteCarlo/MonteCarlo_Main.c
r789 r812 34 34 35 35 int 36 FindPixel(double testQ, double testPhi, double lam, double sdd,36 FindPixel(double testQ, double testPhi, double lam, double yg_d, double sdd, 37 37 double pixSize, double xCtr, double yCtr, long *xPixel, long *yPixel) { 38 38 … … 42 42 qx = testQ*cos(testPhi); 43 43 qy = testQ*sin(testPhi); 44 45 //correct qy for gravity 46 qy += 4.0*pi/lam*(yg_d/sdd/2.0); 47 44 48 45 49 //convert qx,qy to pixel locations relative to # of pixels x, y from center
Note: See TracChangeset
for help on using the changeset viewer.
