Changeset 711 for sans/XOP_Dev/MonteCarlo
 Timestamp:
 Jul 1, 2010 12:18:44 PM (13 years ago)
 Location:
 sans/XOP_Dev/MonteCarlo
 Files:

 3 added
 2 edited
Legend:
 Unmodified
 Added
 Removed

sans/XOP_Dev/MonteCarlo/MonteCarlo.c
r677 r711 8 8 */ 9 9 10 // 11 // uses ran3() 12 // 10 13 11 14 #include "XOPStandardHeaders.h" // Include ANSI headers, Mac headers, IgorXOP.h, XOP.h and XOPSupport.h … … 15 18 //static int gCallSpinProcess = 1; // Set to 1 to all user abort (cmd dot) and background processing. 16 19 17 //////////18 // PROGRAM Monte_SANS19 // PROGRAM simulates multiple SANS.20 // revised 2/12/99 JGB21 // added calculation of random deviate, and 2D 10/2008 SRK22 23 // N1 = NUMBER OF INCIDENT NEUTRONS.24 // N2 = NUMBER INTERACTED IN THE SAMPLE.25 // N3 = NUMBER ABSORBED.26 // THETA = SCATTERING ANGLE.27 28 // works fine in the singlethreaded case.29 //30 /// currently crashes if threaded. apparently something here is either doing an unknown callback, or is accessing31 // a bad place in memory.32 //33 // the operations SpinProcess() and WaveHandleModified() are callbacks and MUST be commented out before threading.34 //  some locks are nonexistent35 //  supposedly safe wave access routines are used36 //37 // random number generators are not threadsafe, and can give less than random results, but is this enough to crash?38 //  a possible workaround is to define multiple versions (poor man's threading)39 //40 //41 //42 20 43 21 … … 70 48 long NDoubleCoherent,NMultipleScatter,isOn,xCtr_long,yCtr_long; 71 49 long NMultipleCoherent,NCoherentEvents; 50 double deltaLam,v1,v2,currWavelength,rsq,fac; //for simulating wavelength distribution 72 51 73 52 // for accessing the 2D wave data, direct method (see the WaveAccess example XOP) … … 156 135 sig_incoh = inputWave[9]; 157 136 sig_sas = inputWave[10]; 137 deltaLam = inputWave[11]; 158 138 xCtr_long = (long)(xCtr+0.5); //round() not defined in VC8 159 139 yCtr_long = (long)(yCtr+0.5); // and is probably wrong to use anyways (returns double!) … … 257 237 } while(rr>r1); 258 238 239 //pick the wavelength out of the wavelength spread, approximate as a gaussian 240 // from NR  pg 288. Needs random # from [0,1]. del is deltaLam/lam (as FWHM) and the 241 // 2.35 converts to a gaussian std dev. 242 do { 243 v1=2.0*ran3(&seed)1.0; 244 v2=2.0*ran3(&seed)1.0; 245 rsq=v1*v1+v2*v2; 246 } while (rsq >= 1.0  rsq == 0.0); 247 fac=sqrt(2.0*log(rsq)/rsq); 248 249 // gset=v1*fac //technically, I'm throwing away one of the two values 250 currWavelength = (v2*fac)*deltaLam*wavelength/2.35 + wavelength; 251 259 252 do { //Scattering Loop, will exit when "done" == 1 260 253 // keep scattering multiple times until the neutron exits the sample … … 297 290 298 291 q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran3(&seed))*delta; 299 theta = q0/2.0/pi* wavelength; //SAS approximation. 1% error at theta=30 degrees292 theta = q0/2.0/pi*currWavelength; //SAS approximation. 1% error at theta=30 degrees 300 293 301 294 find_theta = 1; //always accept … … 339 332 if( index != 0) { //neutron was scattered, figure out where it went 340 333 theta_z = acos(vz); // Angle (= 2theta) WITH respect to z axis. 341 testQ = 2.0*pi*sin(theta_z)/ wavelength;334 testQ = 2.0*pi*sin(theta_z)/currWavelength; 342 335 343 336 // pick a random phi angle, and see if it lands on the detector … … 347 340 348 341 // is it on the detector? 349 FindPixel(testQ,testPhi, wavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel);342 FindPixel(testQ,testPhi,currWavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 350 343 351 344 if(xPixel != 1 && yPixel != 1) { … … 471 464 } 472 465 //////// end of main function for calculating multiple scattering 473 474 int475 FindPixel(double testQ, double testPhi, double lam, double sdd,476 double pixSize, double xCtr, double yCtr, long *xPixel, long *yPixel) {477 478 double theta,dy,dx,qx,qy,pi;479 pi = 4.0*atan(1.0);480 //decompose to qx,qy481 qx = testQ*cos(testPhi);482 qy = testQ*sin(testPhi);483 484 //convert qx,qy to pixel locations relative to # of pixels x, y from center485 theta = 2.0*asin(qy*lam/4.0/pi);486 dy = sdd*tan(theta);487 *yPixel = (long)(yCtr + dy/pixSize+0.5);488 489 theta = 2.0*asin(qx*lam/4.0/pi);490 dx = sdd*tan(theta);491 *xPixel = (long)(xCtr + dx/pixSize+0.5);492 493 //if on detector, return xPix and yPix values, otherwise 1494 if(*yPixel > 127  *yPixel < 0) {495 *yPixel = 1;496 }497 if(*xPixel > 127  *xPixel < 0) {498 *xPixel = 1;499 }500 501 return(0);502 }503 504 505 //calculates new direction (xyz) from an old direction506 //theta and phi don't change507 int508 NewDirection(double *vx, double *vy, double *vz, double theta, double phi) {509 510 int err=0;511 double vx0,vy0,vz0;512 double nx,ny,mag_xy,tx,ty,tz;513 514 //store old direction vector515 vx0 = *vx;516 vy0 = *vy;517 vz0 = *vz;518 519 mag_xy = sqrt(vx0*vx0 + vy0*vy0);520 if(mag_xy < 1e12) {521 //old vector lies along beam direction522 nx = 0;523 ny = 1;524 tx = 1;525 ty = 0;526 tz = 0;527 } else {528 nx = vy0 / mag_xy;529 ny = vx0 / mag_xy;530 tx = vz0*vx0 / mag_xy;531 ty = vz0*vy0 / mag_xy;532 tz = mag_xy ;533 }534 535 //new scattered direction vector536 *vx = cos(phi)*sin(theta)*tx + sin(phi)*sin(theta)*nx + cos(theta)*vx0;537 *vy = cos(phi)*sin(theta)*ty + sin(phi)*sin(theta)*ny + cos(theta)*vy0;538 *vz = cos(phi)*sin(theta)*tz + cos(theta)*vz0;539 540 return(err);541 }542 543 double544 path_len(double aval, double sig_tot) {545 546 double retval;547 548 retval = 1.0*log(1.0aval)/sig_tot;549 550 return(retval);551 }552 553 554 #define IA 16807555 #define IM 2147483647556 #define AM (1.0/IM)557 #define IQ 127773558 #define IR 2836559 #define NTAB 32560 #define NDIV (1+(IM1)/NTAB)561 #define EPS 1.2e7562 #define RNMX (1.0EPS)563 564 float ran1(long *idum)565 {566 int j;567 long k;568 static long iy=0;569 static long iv[NTAB];570 float temp;571 572 if (*idum <= 0  !iy) {573 if ((*idum) < 1) *idum=1;574 else *idum = (*idum);575 for (j=NTAB+7;j>=0;j) {576 k=(*idum)/IQ;577 *idum=IA*(*idumk*IQ)IR*k;578 if (*idum < 0) *idum += IM;579 if (j < NTAB) iv[j] = *idum;580 }581 iy=iv[0];582 }583 k=(*idum)/IQ;584 *idum=IA*(*idumk*IQ)IR*k;585 if (*idum < 0) *idum += IM;586 j=iy/NDIV;587 iy=iv[j];588 iv[j] = *idum;589 if ((temp=AM*iy) > RNMX) return RNMX;590 else return temp;591 }592 #undef IA593 #undef IM594 #undef AM595 #undef IQ596 #undef IR597 #undef NTAB598 #undef NDIV599 #undef EPS600 #undef RNMX601 602 603 //////////a complete copy of ran1(), simply renamed ran1a()604 #define IA 16807605 #define IM 2147483647606 #define AM (1.0/IM)607 #define IQ 127773608 #define IR 2836609 #define NTAB 32610 #define NDIV (1+(IM1)/NTAB)611 #define EPS 1.2e7612 #define RNMX (1.0EPS)613 614 float ran1a(long *idum)615 {616 int j;617 long k;618 static long iy=0;619 static long iv[NTAB];620 float temp;621 622 if (*idum <= 0  !iy) {623 if ((*idum) < 1) *idum=1;624 else *idum = (*idum);625 for (j=NTAB+7;j>=0;j) {626 k=(*idum)/IQ;627 *idum=IA*(*idumk*IQ)IR*k;628 if (*idum < 0) *idum += IM;629 if (j < NTAB) iv[j] = *idum;630 }631 iy=iv[0];632 }633 k=(*idum)/IQ;634 *idum=IA*(*idumk*IQ)IR*k;635 if (*idum < 0) *idum += IM;636 j=iy/NDIV;637 iy=iv[j];638 iv[j] = *idum;639 if ((temp=AM*iy) > RNMX) return RNMX;640 else return temp;641 }642 #undef IA643 #undef IM644 #undef AM645 #undef IQ646 #undef IR647 #undef NTAB648 #undef NDIV649 #undef EPS650 #undef RNMX651 ///////////////652 653 654 ////////////////////////655 #define MBIG 1000000000656 #define MSEED 161803398657 #define MZ 0658 #define FAC (1.0/MBIG)659 660 float ran3(long *idum)661 {662 static int inext,inextp;663 static long ma[56];664 static int iff=0;665 long mj,mk;666 int i,ii,k;667 668 if (*idum < 0  iff == 0) {669 iff=1;670 mj=MSEED(*idum < 0 ? *idum : *idum);671 mj %= MBIG;672 ma[55]=mj;673 mk=1;674 for (i=1;i<=54;i++) {675 ii=(21*i) % 55;676 ma[ii]=mk;677 mk=mjmk;678 if (mk < MZ) mk += MBIG;679 mj=ma[ii];680 }681 for (k=1;k<=4;k++)682 for (i=1;i<=55;i++) {683 ma[i] = ma[1+(i+30) % 55];684 if (ma[i] < MZ) ma[i] += MBIG;685 }686 inext=0;687 inextp=31;688 *idum=1;689 }690 if (++inext == 56) inext=1;691 if (++inextp == 56) inextp=1;692 mj=ma[inext]ma[inextp];693 if (mj < MZ) mj += MBIG;694 ma[inext]=mj;695 return mj*FAC;696 }697 #undef MBIG698 #undef MSEED699 #undef MZ700 #undef FAC701 702 //////////////////////// a complete copy of ran3() renamed ran3a()703 #define MBIG 1000000000704 #define MSEED 161803398705 #define MZ 0706 #define FAC (1.0/MBIG)707 708 float ran3a(long *idum)709 {710 static int inext,inextp;711 static long ma[56];712 static int iff=0;713 long mj,mk;714 int i,ii,k;715 716 if (*idum < 0  iff == 0) {717 iff=1;718 mj=MSEED(*idum < 0 ? *idum : *idum);719 mj %= MBIG;720 ma[55]=mj;721 mk=1;722 for (i=1;i<=54;i++) {723 ii=(21*i) % 55;724 ma[ii]=mk;725 mk=mjmk;726 if (mk < MZ) mk += MBIG;727 mj=ma[ii];728 }729 for (k=1;k<=4;k++)730 for (i=1;i<=55;i++) {731 ma[i] = ma[1+(i+30) % 55];732 if (ma[i] < MZ) ma[i] += MBIG;733 }734 inext=0;735 inextp=31;736 *idum=1;737 }738 if (++inext == 56) inext=1;739 if (++inextp == 56) inextp=1;740 mj=ma[inext]ma[inextp];741 if (mj < MZ) mj += MBIG;742 ma[inext]=mj;743 return mj*FAC;744 }745 #undef MBIG746 #undef MSEED747 #undef MZ748 #undef FAC749 750 751 // returns the interpolated point value in xx[0,n1] that has the value x752 double locate_interp(double xx[], long n, double x)753 {754 unsigned long ju,jm,jl,j;755 int ascnd;756 double pt;757 758 // char buf[256];759 760 jl=0;761 ju=n1;762 ascnd=(xx[n1] > xx[0]);763 while (jujl > 1) {764 jm=(ju+jl) >> 1;765 if (x > xx[jm] == ascnd)766 jl=jm;767 else768 ju=jm;769 }770 j=jl; // the point I want is between xx[j] and xx[j+1]771 pt = (x xx[j])/(xx[j+1]  xx[j]); //fractional distance, using linear interpolation772 773 // sprintf(buf, "x = %g, j= %ld, pt = %g\r",x,j,pt);774 // XOPNotice(buf);775 776 return(pt+(double)j);777 }778 779 780 781 782 /////////////////////////////783 /* RegisterFunction()784 785 Igor calls this at startup time to find the address of the786 XFUNCs added by this XOP. See XOP manual regarding "Direct XFUNCs".787 */788 static long789 RegisterFunction()790 {791 int funcIndex;792 793 funcIndex = GetXOPItem(0); // Which function is Igor asking about?794 switch (funcIndex) {795 case 0: //796 return((long)Monte_SANSX);797 break;798 case 1: //799 return((long)Monte_SANSX2);800 break;801 case 2: //802 return((long)DebyeSpheresX);803 break;804 case 3: //805 return((long)Monte_SANSX3);806 break;807 case 4: //808 return((long)Monte_SANSX4);809 break;810 }811 return(NIL);812 }813 814 /* XOPEntry()815 816 This is the entry point from the host application to the XOP for all messages after the817 INIT message.818 */819 static void820 XOPEntry(void)821 {822 long result = 0;823 824 switch (GetXOPMessage()) {825 case FUNCADDRS:826 result = RegisterFunction();827 break;828 }829 SetXOPResult(result);830 }831 832 /* main(ioRecHandle)833 834 This is the initial entry point at which the host application calls XOP.835 The message sent by the host must be INIT.836 main() does any necessary initialization and then sets the XOPEntry field of the837 ioRecHandle to the address to be called for future messages.838 */839 HOST_IMPORT void840 main(IORecHandle ioRecHandle)841 {842 XOPInit(ioRecHandle); // Do standard XOP initialization.843 SetXOPEntry(XOPEntry); // Set entry point for future calls.844 845 if (igorVersion < 600) // Requires Igor Pro 6.00 or later.846 SetXOPResult(OLD_IGOR); // OLD_IGOR is defined in WaveAccess.h and there are corresponding error strings in WaveAccess.r and WaveAccessWinCustom.rc.847 else848 SetXOPResult(0L);849 } 
sans/XOP_Dev/MonteCarlo/MonteCarlo2.c
r673 r711 1 1 /* 2 * MonteCarlo .c3 * SANS Analysis2 * MonteCarlo2.c 3 * SANSMonteCarlo 4 4 * 5 * Created by Steve Kline on 10/16/08.6 * Copyright 20 08 __MyCompanyName__. All rights reserved.5 * Created by Steve Kline on 7/1/10. 6 * Copyright 2010 NCNR. All rights reserved. 7 7 * 8 8 */ … … 51 51 long NDoubleCoherent,NMultipleScatter,isOn,xCtr_long,yCtr_long; 52 52 long NMultipleCoherent,NCoherentEvents; 53 double deltaLam,v1,v2,currWavelength,rsq,fac; //for simulating wavelength distribution 53 54 54 55 … … 138 139 sig_incoh = inputWave[9]; 139 140 sig_sas = inputWave[10]; 141 deltaLam = inputWave[11]; 140 142 xCtr_long = (long)(xCtr+0.5); 141 143 yCtr_long = (long)(yCtr+0.5); … … 239 241 } while(rr>r1); 240 242 243 //pick the wavelength out of the wavelength spread, approximate as a gaussian 244 // from NR  pg 288. Needs random # from [0,1]. del is deltaLam/lam (as FWHM) and the 245 // 2.35 converts to a gaussian std dev. 246 do { 247 v1=2.0*ran1(&seed)1.0; 248 v2=2.0*ran1(&seed)1.0; 249 rsq=v1*v1+v2*v2; 250 } while (rsq >= 1.0  rsq == 0.0); 251 fac=sqrt(2.0*log(rsq)/rsq); 252 253 // gset=v1*fac //technically, I'm throwing away one of the two values 254 currWavelength = (v2*fac)*deltaLam*wavelength/2.35 + wavelength; 255 241 256 do { //Scattering Loop, will exit when "done" == 1 242 257 // keep scattering multiple times until the neutron exits the sample … … 279 294 280 295 q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran1(&seed))*delta; 281 theta = q0/2/pi* wavelength; //SAS approximation296 theta = q0/2/pi*currWavelength; //SAS approximation 282 297 283 298 find_theta = 1; //always accept … … 321 336 if( index != 0) { //neutron was scattered, figure out where it went 322 337 theta_z = acos(vz); // Angle (= 2theta) WITH respect to z axis. 323 testQ = 2*pi*sin(theta_z)/ wavelength;338 testQ = 2*pi*sin(theta_z)/currWavelength; 324 339 325 340 // pick a random phi angle, and see if it lands on the detector … … 329 344 330 345 // is it on the detector? 331 FindPixel(testQ,testPhi, wavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel);346 FindPixel(testQ,testPhi,currWavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 332 347 333 348 if(xPixel != 1 && yPixel != 1) { … … 451 466 return(0); 452 467 } 453 //////// end of X2 454 455 456 457 //////////////// X3 using ran3a 458 int 459 Monte_SANSX3(MC_ParamsPtr p) { 460 double *inputWave; /* pointer to double precision wave data */ 461 double *ran_dev; /* pointer to double precision wave data */ 462 double *nt; /* pointer to double precision wave data */ 463 double *j1; /* pointer to double precision wave data */ 464 double *j2; /* pointer to double precision wave data */ 465 double *nn; /* pointer to double precision wave data */ 466 // double *MC_linear_data; /* pointer to double precision wave data */ 467 double *results; /* pointer to double precision wave data */ 468 double retVal; //return value 469 470 long imon; 471 double r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas; 472 long ind,index,n_index; 473 double qmax,theta_max,q0,zpow; 474 long n1,n2,n3; 475 double dth,zz,xx,yy,phi; 476 double theta,ran,ll,rr; 477 long done,find_theta,err; //used as logicals 478 long xPixel,yPixel; 479 double vx,vy,vz,theta_z; 480 double sig_abs,ratio,sig_total; 481 double testQ,testPhi,left,delta,dummy,pi; 482 double sigabs_0,num_bins; 483 long NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent; 484 long NDoubleCoherent,NMultipleScatter,isOn,xCtr_long,yCtr_long; 485 long NMultipleCoherent,NCoherentEvents; 486 487 488 // for accessing the 2D wave data, direct method (see the WaveAccess example XOP) 489 waveHndl wavH; 490 // int waveType,hState; 491 long numDimensions; 492 long dimensionSizes[MAX_DIMENSIONS+1]; 493 // char* dataStartPtr; 494 // long dataOffset; 495 // long numRows, numColumns; 496 long numRows_ran_dev; 497 // double *dp0, *dp; 498 double value[2]; // Pointers used for double data. 499 long seed; 500 long indices[MAX_DIMENSIONS]; 501 502 // char buf[256]; 503 504 /* check that wave handles are all valid */ 505 if (p>inputWaveH == NIL) { 506 SetNaN64(&p>retVal); /* return NaN if wave is not valid */ 507 return(NON_EXISTENT_WAVE); 508 } 509 if (p>ran_devH == NIL) { 510 SetNaN64(&p>retVal); /* return NaN if wave is not valid */ 511 return(NON_EXISTENT_WAVE); 512 } 513 if (p>ntH == NIL) { 514 SetNaN64(&p>retVal); /* return NaN if wave is not valid */ 515 return(NON_EXISTENT_WAVE); 516 } 517 if (p>j1H == NIL) { 518 SetNaN64(&p>retVal); /* return NaN if wave is not valid */ 519 return(NON_EXISTENT_WAVE); 520 } 521 if (p>j2H == NIL) { 522 SetNaN64(&p>retVal); /* return NaN if wave is not valid */ 523 return(NON_EXISTENT_WAVE); 524 } 525 if (p>nnH == NIL) { 526 SetNaN64(&p>retVal); /* return NaN if wave is not valid */ 527 return(NON_EXISTENT_WAVE); 528 } 529 if (p>MC_linear_dataH == NIL) { 530 SetNaN64(&p>retVal); /* return NaN if wave is not valid */ 531 return(NON_EXISTENT_WAVE); 532 } 533 if (p>resultsH == NIL) { 534 SetNaN64(&p>retVal); /* return NaN if wave is not valid */ 535 return(NON_EXISTENT_WAVE); 536 } 537 538 p>retVal = 0; 539 540 // trusting that all inputs are DOUBLE PRECISION WAVES!!! 541 inputWave = WaveData(p>inputWaveH); 542 ran_dev = WaveData(p>ran_devH); 543 nt = WaveData(p>ntH); 544 j1 = WaveData(p>j1H); 545 j2 = WaveData(p>j2H); 546 nn = WaveData(p>nnH); 547 // MC_linear_data = WaveData(p>MC_linear_dataH); 548 results = WaveData(p>resultsH); 549 550 seed = (long)results[0]; 551 552 // sprintf(buf, "input seed = %ld\r", seed); 553 // XOPNotice(buf); 554 555 if(seed >= 0) { 556 seed = 1234509876; 557 } 558 559 dummy = ran3a(&seed); //initialize the random sequence by passing in a negative value 560 seed = 12348765; //nonnegative after that does nothing 561 562 imon = (int)inputWave[0]; 563 r1 = inputWave[1]; 564 r2 = inputWave[2]; 565 xCtr = inputWave[3]; 566 yCtr = inputWave[4]; 567 sdd = inputWave[5]; 568 pixSize = inputWave[6]; 569 thick = inputWave[7]; 570 wavelength = inputWave[8]; 571 sig_incoh = inputWave[9]; 572 sig_sas = inputWave[10]; 573 xCtr_long = (long)(xCtr+0.5); 574 yCtr_long = (long)(yCtr+0.5); 575 576 dummy = MDGetWaveScaling(p>ran_devH, 0, &delta, &left); //0 is the rows 577 if (retVal = MDGetWaveDimensions(p>ran_devH, &numDimensions, dimensionSizes)) 578 return retVal; 579 numRows_ran_dev = dimensionSizes[0]; 580 581 pi = 4.0*atan(1.0); 582 583 // access the 2D wave data for writing using the direct method 584 wavH = p>MC_linear_dataH; 585 if (wavH == NIL) 586 return NOWAV; 587 588 // waveType = WaveType(wavH); 589 // if (waveType & NT_CMPLX) 590 // return NO_COMPLEX_WAVE; 591 // if (waveType==TEXT_WAVE_TYPE) 592 // return NUMERIC_ACCESS_ON_TEXT_WAVE; 593 // if (retVal = MDGetWaveDimensions(wavH, &numDimensions, dimensionSizes)) 594 // return retVal; 595 // numRows = dimensionSizes[0]; 596 // numColumns = dimensionSizes[1]; 597 598 // if (retVal = MDAccessNumericWaveData(wavH, kMDWaveAccessMode0, &dataOffset)) 599 // return retVal; 600 601 // hState = MoveLockHandle(wavH); // So wave data can't move. Remember to call HSetState when done. 602 // dataStartPtr = (char*)(*wavH) + dataOffset; 603 // dp0 = (double*)dataStartPtr; // Pointer to the start of the 2D wave data. 604 605 //scattering power and maximum qvalue to bin 606 // zpow = .1 //scattering power, calculated below 607 qmax = 4.0*pi/wavelength; //maximum Q to bin 1D data. (A1) (not really used) 608 sigabs_0 = 0.0; // ignore absorption cross section/wavelength [1/(cm A)] 609 n_index = 50; // maximum number of scattering events per neutron 610 num_bins = 200; //number of 1D bins (not really used) 611 612 //c total SAS crosssection 613 // 614 zpow = sig_sas*thick; //since I now calculate the sig_sas from the model 615 sig_abs = sigabs_0 * wavelength; 616 sig_total = sig_abs + sig_sas + sig_incoh; 617 // Print "The TOTAL XSECTION. (CM1) is ",sig_total 618 // Print "The TOTAL SAS XSECTION. (CM1) is ",sig_sas 619 // results[0] = sig_total; 620 // results[1] = sig_sas; 621 // RATIO = SIG_ABS / SIG_TOTAL 622 ratio = sig_incoh / sig_total; 623 624 theta_max = wavelength*qmax/(2*pi); 625 //C SET ThetaSTEP SIZE. 626 dth = theta_max/num_bins; 627 // Print "theta bin size = dth = ",dth 628 629 //C INITIALIZE COUNTERS. 630 n1 = 0; 631 n2 = 0; 632 n3 = 0; 633 NSingleIncoherent = 0; 634 NSingleCoherent = 0; 635 NDoubleCoherent = 0; 636 NMultipleScatter = 0; 637 NScatterEvents = 0; 638 NMultipleCoherent = 0; 639 NCoherentEvents = 0; 640 641 isOn = 0; 642 643 //C MONITOR LOOP  looping over the number of incedent neutrons 644 //note that zz, is the zposition in the sample  NOT the scattering power 645 // NOW, start the loop, throwing neutrons at the sample. 646 do { 647 ////SpinProcess() IS A CALLBACK, and not good for Threading! 648 // if ((n1 % 1000 == 0) && gCallSpinProcess && SpinProcess()) { // Spins cursor and allows background processing. 649 // retVal = 1; // User aborted. 650 // break; 651 // } 652 653 vx = 0.0; // Initialize direction vector. 654 vy = 0.0; 655 vz = 1.0; 656 657 theta = 0.0; // Initialize scattering angle. 658 phi = 0.0; // Intialize azimuthal angle. 659 n1 += 1; // Increment total number neutrons counter. 660 done = 0; // True when neutron is absorbed or when scattered out of the sample. 661 index = 0; // Set counter for number of scattering events. 662 zz = 0.0; // Set entering dimension of sample. 663 incoherentEvent = 0; 664 coherentEvent = 0; 665 666 do { // Makes sure position is within circle. 667 ran = ran3a(&seed); //[0,1] 668 xx = 2.0*r1*(ran0.5); //X beam position of neutron entering sample. 669 ran = ran3a(&seed); //[0,1] 670 yy = 2.0*r1*(ran0.5); //Y beam position ... 671 rr = sqrt(xx*xx+yy*yy); //Radial position of neutron in incident beam. 672 } while(rr>r1); 673 674 do { //Scattering Loop, will exit when "done" == 1 675 // keep scattering multiple times until the neutron exits the sample 676 ran = ran3a(&seed); //[0,1] RANDOM NUMBER FOR DETERMINING PATH LENGTH 677 ll = path_len(ran,sig_total); 678 //Determine new scattering direction vector. 679 err = NewDirection(&vx,&vy,&vz,theta,phi); //vx,vy,vz updated, theta, phi unchanged by function 680 681 //X,Y,ZPOSITION OF SCATTERING EVENT. 682 xx += ll*vx; 683 yy += ll*vy; 684 zz += ll*vz; 685 rr = sqrt(xx*xx+yy*yy); //radial position of scattering event. 686 687 //sprintf(buf, "xx,yy,zz,vx,vy,vz,ll = %g %g %g %g %g %g %g\r",xx,yy,zz,vx,vy,vz,ll); 688 //XOPNotice(buf); 689 690 //Check whether interaction occurred within sample volume. 691 if (((zz > 0.0) && (zz < thick)) && (rr < r2)) { 692 //NEUTRON INTERACTED. 693 //sprintf(buf,"neutron interacted\r"); 694 //XOPNotice(buf); 695 696 index += 1; //Increment counter of scattering events. 697 if (index == 1) { 698 n2 += 1; //Increment # of scat. neutrons 699 } 700 ran = ran3a(&seed); //[0,1] 701 //Split neutron interactions into scattering and absorption events 702 if (ran > ratio ) { //C NEUTRON SCATTERED coherently 703 //sprintf(buf,"neutron scatters coherently\r"); 704 //XOPNotice(buf); 705 coherentEvent += 1; 706 find_theta = 0; //false 707 do { 708 // pick a qvalue from the deviate function 709 // pnt2x truncates the point to an integer before returning the x 710 // so get it from the wave scaling instead 711 // q0 =left + binarysearchinterp(ran_dev,ran3a(seed))*delta; 712 713 q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran3a(&seed))*delta; 714 theta = q0/2/pi*wavelength; //SAS approximation 715 716 find_theta = 1; //always accept 717 718 //sprintf(buf, "after locate_interp call q0 = %g, theta = %g,left = %g,delta = %g\r",q0,theta,left,delta); 719 //XOPNotice(buf); 720 721 } while(!find_theta); 722 723 ran = ran3a(&seed); //[0,1] 724 phi = 2.0*pi*ran; //Chooses azimuthal scattering angle. 725 } else { 726 //NEUTRON scattered incoherently 727 //sprintf(buf,"neutron scatters incoherent\r"); 728 //XOPNotice(buf); 729 incoherentEvent += 1; 730 // phi and theta are random over the entire sphere of scattering 731 // !can't just choose random theta and phi, won't be random over sphere solid angle 732 733 ran = ran3a(&seed); //[0,1] 734 theta = acos(2.0*ran1); 735 736 ran = ran3a(&seed); //[0,1] 737 phi = 2.0*pi*ran; //Chooses azimuthal scattering angle. 738 } //(ran > ratio) 739 } else { 740 //NEUTRON ESCAPES FROM SAMPLE  bin it somewhere 741 done = 1; //done = true, will exit from loop 742 //Increment #scattering events array 743 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 744 indices[0] =index; //this sets access to nn[index] 745 if (index <= n_index) { 746 if (retVal = MDGetNumericWavePointValue(p>nnH, indices, value)) 747 return retVal; 748 value[0] += 1; // add one to the value 749 if (retVal = MDSetNumericWavePointValue(p>nnH, indices, value)) 750 return retVal; 751 // nn[index] += 1; 752 } 753 754 if( index != 0) { //neutron was scattered, figure out where it went 755 theta_z = acos(vz); // Angle (= 2theta) WITH respect to z axis. 756 testQ = 2*pi*sin(theta_z)/wavelength; 757 758 // pick a random phi angle, and see if it lands on the detector 759 // since the scattering is isotropic, I can safely pick a new, random value 760 // this would not be true if simulating anisotropic scattering. 761 testPhi = ran3a(&seed)*2*pi; 762 763 // is it on the detector? 764 FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 765 766 if(xPixel != 1 && yPixel != 1) { 767 isOn += 1; 768 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 769 indices[0] = xPixel; 770 indices[1] = yPixel; 771 if (retVal = MDGetNumericWavePointValue(wavH, indices, value)) 772 return retVal; 773 value[0] += 1; // Real part 774 if (retVal = MDSetNumericWavePointValue(wavH, indices, value)) 775 return retVal; 776 //if(index==1) // only the single scattering events 777 //dp = dp0 + xPixel + yPixel*numColumns; //offset the pointer to the exact memory location 778 //*dp += 1; //increment the value there 779 //endif 780 } 781 782 783 /* is this causing me a problem since I'm not locking these? Probably not, since it crashes even if I comment these out... */ 784 if(theta_z < theta_max) { 785 //Choose index for scattering angle array. 786 //IND = NINT(THETA_z/DTH + 0.4999999) 787 ind = (long)(theta_z/dth + 0.4999999); //round is eqivalent to nint() 788 nt[ind] += 1; //Increment bin for angle. 789 //Increment angle array for single scattering events. 790 if (index == 1) { 791 j1[ind] += 1; 792 } 793 //Increment angle array for double scattering events. 794 if (index == 2) { 795 j2[ind] += 1; 796 } 797 } 798 /**/ 799 800 // increment all of the counters now since done==1 here and I'm sure to exit and get another neutron 801 NScatterEvents += index; //total number of scattering events 802 if(index == 1 && incoherentEvent == 1) { 803 NSingleIncoherent += 1; 804 } 805 if(index == 1 && coherentEvent == 1) { 806 NSingleCoherent += 1; 807 } 808 if(index == 2 && coherentEvent == 1 && incoherentEvent == 0) { 809 NDoubleCoherent += 1; 810 } 811 if(index > 1) { 812 NMultipleScatter += 1; 813 } 814 if(coherentEvent >= 1 && incoherentEvent == 0) { 815 NCoherentEvents += 1; 816 } 817 if(coherentEvent > 1 && incoherentEvent == 0) { 818 NMultipleCoherent += 1; 819 } 820 821 } else { // index was zero, neutron must be transmitted, so just increment the proper counters and data 822 isOn += 1; 823 nt[0] += 1; 824 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 825 //indices[0] = xCtr_long; //don't put everything in one pixel 826 //indices[1] = yCtr_long; 827 indices[0] = (long)(xCtr+xx/pixSize+0.5); 828 indices[1] = (long)(yCtr+yy/pixSize+0.5); 829 // check for valid indices  got an XOP error, probably from here 830 if(indices[0] > 127) indices[0] = 127; 831 if(indices[0] < 0) indices[0] = 0; 832 if(indices[1] > 127) indices[1] = 127; 833 if(indices[1] < 0) indices[1] = 0; 834 835 if (retVal = MDGetNumericWavePointValue(wavH, indices, value)) 836 return retVal; 837 value[0] += 1; // Real part 838 if (retVal = MDSetNumericWavePointValue(wavH, indices, value)) 839 return retVal; 840 } 841 } 842 } while (!done); 843 } while(n1 < imon); 844 845 // assign the results to the wave 846 847 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 848 value[0] = (double)n1; 849 indices[0] = 0; 850 if (retVal = MDSetNumericWavePointValue(p>resultsH, indices, value)) 851 return retVal; 852 value[0] = (double)n2; 853 indices[0] = 1; 854 if (retVal = MDSetNumericWavePointValue(p>resultsH, indices, value)) 855 return retVal; 856 value[0] = (double)isOn; 857 indices[0] = 2; 858 if (retVal = MDSetNumericWavePointValue(p>resultsH, indices, value)) 859 return retVal; 860 value[0] = (double)NScatterEvents; 861 indices[0] = 3; 862 if (retVal = MDSetNumericWavePointValue(p>resultsH, indices, value)) 863 return retVal; 864 value[0] = (double)NSingleCoherent; 865 indices[0] = 4; 866 if (retVal = MDSetNumericWavePointValue(p>resultsH, indices, value)) 867 return retVal; 868 value[0] = (double)NMultipleCoherent; 869 indices[0] = 5; 870 if (retVal = MDSetNumericWavePointValue(p>resultsH, indices, value)) 871 return retVal; 872 value[0] = (double)NMultipleScatter; 873 indices[0] = 6; 874 if (retVal = MDSetNumericWavePointValue(p>resultsH, indices, value)) 875 return retVal; 876 value[0] = (double)NCoherentEvents; 877 indices[0] = 7; 878 if (retVal = MDSetNumericWavePointValue(p>resultsH, indices, value)) 879 return retVal; 880 881 // HSetState((Handle)wavH, hState); //release the handle of the 2D data wave 882 // WaveHandleModified(wavH); // Inform Igor that we have changed the wave. (CALLBACK! needed, but not allowed in Threading) 883 884 return(0); 885 } 886 //////// end of X3 887 888 889 ///////////////// X4 using ran1a() 890 int 891 Monte_SANSX4(MC_ParamsPtr p) { 892 double *inputWave; /* pointer to double precision wave data */ 893 double *ran_dev; /* pointer to double precision wave data */ 894 double *nt; /* pointer to double precision wave data */ 895 double *j1; /* pointer to double precision wave data */ 896 double *j2; /* pointer to double precision wave data */ 897 double *nn; /* pointer to double precision wave data */ 898 // double *MC_linear_data; /* pointer to double precision wave data */ 899 double *results; /* pointer to double precision wave data */ 900 double retVal; //return value 901 902 long imon; 903 double r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas; 904 long ind,index,n_index; 905 double qmax,theta_max,q0,zpow; 906 long n1,n2,n3; 907 double dth,zz,xx,yy,phi; 908 double theta,ran,ll,rr; 909 long done,find_theta,err; //used as logicals 910 long xPixel,yPixel; 911 double vx,vy,vz,theta_z; 912 double sig_abs,ratio,sig_total; 913 double testQ,testPhi,left,delta,dummy,pi; 914 double sigabs_0,num_bins; 915 long NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent; 916 long NDoubleCoherent,NMultipleScatter,isOn,xCtr_long,yCtr_long; 917 long NMultipleCoherent,NCoherentEvents; 918 919 920 // for accessing the 2D wave data, direct method (see the WaveAccess example XOP) 921 waveHndl wavH; 922 // int waveType,hState; 923 long numDimensions; 924 long dimensionSizes[MAX_DIMENSIONS+1]; 925 // char* dataStartPtr; 926 // long dataOffset; 927 // long numRows, numColumns; 928 long numRows_ran_dev; 929 // double *dp0, *dp; 930 double value[2]; // Pointers used for double data. 931 long seed; 932 long indices[MAX_DIMENSIONS]; 933 934 // char buf[256]; 935 936 /* check that wave handles are all valid */ 937 if (p>inputWaveH == NIL) { 938 SetNaN64(&p>retVal); /* return NaN if wave is not valid */ 939 return(NON_EXISTENT_WAVE); 940 } 941 if (p>ran_devH == NIL) { 942 SetNaN64(&p>retVal); /* return NaN if wave is not valid */ 943 return(NON_EXISTENT_WAVE); 944 } 945 if (p>ntH == NIL) { 946 SetNaN64(&p>retVal); /* return NaN if wave is not valid */ 947 return(NON_EXISTENT_WAVE); 948 } 949 if (p>j1H == NIL) { 950 SetNaN64(&p>retVal); /* return NaN if wave is not valid */ 951 return(NON_EXISTENT_WAVE); 952 } 953 if (p>j2H == NIL) { 954 SetNaN64(&p>retVal); /* return NaN if wave is not valid */ 955 return(NON_EXISTENT_WAVE); 956 } 957 if (p>nnH == NIL) { 958 SetNaN64(&p>retVal); /* return NaN if wave is not valid */ 959 return(NON_EXISTENT_WAVE); 960 } 961 if (p>MC_linear_dataH == NIL) { 962 SetNaN64(&p>retVal); /* return NaN if wave is not valid */ 963 return(NON_EXISTENT_WAVE); 964 } 965 if (p>resultsH == NIL) { 966 SetNaN64(&p>retVal); /* return NaN if wave is not valid */ 967 return(NON_EXISTENT_WAVE); 968 } 969 970 p>retVal = 0; 971 972 // trusting that all inputs are DOUBLE PRECISION WAVES!!! 973 inputWave = WaveData(p>inputWaveH); 974 ran_dev = WaveData(p>ran_devH); 975 nt = WaveData(p>ntH); 976 j1 = WaveData(p>j1H); 977 j2 = WaveData(p>j2H); 978 nn = WaveData(p>nnH); 979 // MC_linear_data = WaveData(p>MC_linear_dataH); 980 results = WaveData(p>resultsH); 981 982 seed = (long)results[0]; 983 984 // sprintf(buf, "input seed = %ld\r", seed); 985 // XOPNotice(buf); 986 987 if(seed >= 0) { 988 seed = 1234509876; 989 } 990 991 dummy = ran1a(&seed); //initialize the random sequence by passing in a negative value 992 seed = 12348765; //nonnegative after that does nothing 993 994 imon = (int)inputWave[0]; 995 r1 = inputWave[1]; 996 r2 = inputWave[2]; 997 xCtr = inputWave[3]; 998 yCtr = inputWave[4]; 999 sdd = inputWave[5]; 1000 pixSize = inputWave[6]; 1001 thick = inputWave[7]; 1002 wavelength = inputWave[8]; 1003 sig_incoh = inputWave[9]; 1004 sig_sas = inputWave[10]; 1005 xCtr_long = (long)(xCtr+0.5); 1006 yCtr_long = (long)(yCtr+0.5); 1007 1008 dummy = MDGetWaveScaling(p>ran_devH, 0, &delta, &left); //0 is the rows 1009 if (retVal = MDGetWaveDimensions(p>ran_devH, &numDimensions, dimensionSizes)) 1010 return retVal; 1011 numRows_ran_dev = dimensionSizes[0]; 1012 1013 pi = 4.0*atan(1.0); 1014 1015 // access the 2D wave data for writing using the direct method 1016 wavH = p>MC_linear_dataH; 1017 if (wavH == NIL) 1018 return NOWAV; 1019 1020 // waveType = WaveType(wavH); 1021 // if (waveType & NT_CMPLX) 1022 // return NO_COMPLEX_WAVE; 1023 // if (waveType==TEXT_WAVE_TYPE) 1024 // return NUMERIC_ACCESS_ON_TEXT_WAVE; 1025 // if (retVal = MDGetWaveDimensions(wavH, &numDimensions, dimensionSizes)) 1026 // return retVal; 1027 // numRows = dimensionSizes[0]; 1028 // numColumns = dimensionSizes[1]; 1029 1030 // if (retVal = MDAccessNumericWaveData(wavH, kMDWaveAccessMode0, &dataOffset)) 1031 // return retVal; 1032 1033 // hState = MoveLockHandle(wavH); // So wave data can't move. Remember to call HSetState when done. 1034 // dataStartPtr = (char*)(*wavH) + dataOffset; 1035 // dp0 = (double*)dataStartPtr; // Pointer to the start of the 2D wave data. 1036 1037 //scattering power and maximum qvalue to bin 1038 // zpow = .1 //scattering power, calculated below 1039 qmax = 4.0*pi/wavelength; //maximum Q to bin 1D data. (A1) (not really used) 1040 sigabs_0 = 0.0; // ignore absorption cross section/wavelength [1/(cm A)] 1041 n_index = 50; // maximum number of scattering events per neutron 1042 num_bins = 200; //number of 1D bins (not really used) 1043 1044 //c total SAS crosssection 1045 // 1046 zpow = sig_sas*thick; //since I now calculate the sig_sas from the model 1047 sig_abs = sigabs_0 * wavelength; 1048 sig_total = sig_abs + sig_sas + sig_incoh; 1049 // Print "The TOTAL XSECTION. (CM1) is ",sig_total 1050 // Print "The TOTAL SAS XSECTION. (CM1) is ",sig_sas 1051 // results[0] = sig_total; 1052 // results[1] = sig_sas; 1053 // RATIO = SIG_ABS / SIG_TOTAL 1054 ratio = sig_incoh / sig_total; 1055 1056 theta_max = wavelength*qmax/(2*pi); 1057 //C SET ThetaSTEP SIZE. 1058 dth = theta_max/num_bins; 1059 // Print "theta bin size = dth = ",dth 1060 1061 //C INITIALIZE COUNTERS. 1062 n1 = 0; 1063 n2 = 0; 1064 n3 = 0; 1065 NSingleIncoherent = 0; 1066 NSingleCoherent = 0; 1067 NDoubleCoherent = 0; 1068 NMultipleScatter = 0; 1069 NScatterEvents = 0; 1070 NMultipleCoherent = 0; 1071 NCoherentEvents = 0; 1072 1073 isOn = 0; 1074 1075 //C MONITOR LOOP  looping over the number of incedent neutrons 1076 //note that zz, is the zposition in the sample  NOT the scattering power 1077 // NOW, start the loop, throwing neutrons at the sample. 1078 do { 1079 ////SpinProcess() IS A CALLBACK, and not good for Threading! 1080 // if ((n1 % 1000 == 0) && gCallSpinProcess && SpinProcess()) { // Spins cursor and allows background processing. 1081 // retVal = 1; // User aborted. 1082 // break; 1083 // } 1084 1085 vx = 0.0; // Initialize direction vector. 1086 vy = 0.0; 1087 vz = 1.0; 1088 1089 theta = 0.0; // Initialize scattering angle. 1090 phi = 0.0; // Intialize azimuthal angle. 1091 n1 += 1; // Increment total number neutrons counter. 1092 done = 0; // True when neutron is absorbed or when scattered out of the sample. 1093 index = 0; // Set counter for number of scattering events. 1094 zz = 0.0; // Set entering dimension of sample. 1095 incoherentEvent = 0; 1096 coherentEvent = 0; 1097 1098 do { // Makes sure position is within circle. 1099 ran = ran1a(&seed); //[0,1] 1100 xx = 2.0*r1*(ran0.5); //X beam position of neutron entering sample. 1101 ran = ran1a(&seed); //[0,1] 1102 yy = 2.0*r1*(ran0.5); //Y beam position ... 1103 rr = sqrt(xx*xx+yy*yy); //Radial position of neutron in incident beam. 1104 } while(rr>r1); 1105 1106 do { //Scattering Loop, will exit when "done" == 1 1107 // keep scattering multiple times until the neutron exits the sample 1108 ran = ran1a(&seed); //[0,1] RANDOM NUMBER FOR DETERMINING PATH LENGTH 1109 ll = path_len(ran,sig_total); 1110 //Determine new scattering direction vector. 1111 err = NewDirection(&vx,&vy,&vz,theta,phi); //vx,vy,vz updated, theta, phi unchanged by function 1112 1113 //X,Y,ZPOSITION OF SCATTERING EVENT. 1114 xx += ll*vx; 1115 yy += ll*vy; 1116 zz += ll*vz; 1117 rr = sqrt(xx*xx+yy*yy); //radial position of scattering event. 1118 1119 //sprintf(buf, "xx,yy,zz,vx,vy,vz,ll = %g %g %g %g %g %g %g\r",xx,yy,zz,vx,vy,vz,ll); 1120 //XOPNotice(buf); 1121 1122 //Check whether interaction occurred within sample volume. 1123 if (((zz > 0.0) && (zz < thick)) && (rr < r2)) { 1124 //NEUTRON INTERACTED. 1125 //sprintf(buf,"neutron interacted\r"); 1126 //XOPNotice(buf); 1127 1128 index += 1; //Increment counter of scattering events. 1129 if (index == 1) { 1130 n2 += 1; //Increment # of scat. neutrons 1131 } 1132 ran = ran1a(&seed); //[0,1] 1133 //Split neutron interactions into scattering and absorption events 1134 if (ran > ratio ) { //C NEUTRON SCATTERED coherently 1135 //sprintf(buf,"neutron scatters coherently\r"); 1136 //XOPNotice(buf); 1137 coherentEvent += 1; 1138 find_theta = 0; //false 1139 do { 1140 // pick a qvalue from the deviate function 1141 // pnt2x truncates the point to an integer before returning the x 1142 // so get it from the wave scaling instead 1143 // q0 =left + binarysearchinterp(ran_dev,ran1a(seed))*delta; 1144 1145 q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran1a(&seed))*delta; 1146 theta = q0/2/pi*wavelength; //SAS approximation 1147 1148 find_theta = 1; //always accept 1149 1150 //sprintf(buf, "after locate_interp call q0 = %g, theta = %g,left = %g,delta = %g\r",q0,theta,left,delta); 1151 //XOPNotice(buf); 1152 1153 } while(!find_theta); 1154 1155 ran = ran1a(&seed); //[0,1] 1156 phi = 2.0*pi*ran; //Chooses azimuthal scattering angle. 1157 } else { 1158 //NEUTRON scattered incoherently 1159 //sprintf(buf,"neutron scatters incoherent\r"); 1160 //XOPNotice(buf); 1161 incoherentEvent += 1; 1162 // phi and theta are random over the entire sphere of scattering 1163 // !can't just choose random theta and phi, won't be random over sphere solid angle 1164 1165 ran = ran1a(&seed); //[0,1] 1166 theta = acos(2.0*ran1); 1167 1168 ran = ran1a(&seed); //[0,1] 1169 phi = 2.0*pi*ran; //Chooses azimuthal scattering angle. 1170 } //(ran > ratio) 1171 } else { 1172 //NEUTRON ESCAPES FROM SAMPLE  bin it somewhere 1173 done = 1; //done = true, will exit from loop 1174 //Increment #scattering events array 1175 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 1176 indices[0] =index; //this sets access to nn[index] 1177 if (index <= n_index) { 1178 if (retVal = MDGetNumericWavePointValue(p>nnH, indices, value)) 1179 return retVal; 1180 value[0] += 1; // add one to the value 1181 if (retVal = MDSetNumericWavePointValue(p>nnH, indices, value)) 1182 return retVal; 1183 // nn[index] += 1; 1184 } 1185 1186 if( index != 0) { //neutron was scattered, figure out where it went 1187 theta_z = acos(vz); // Angle (= 2theta) WITH respect to z axis. 1188 testQ = 2*pi*sin(theta_z)/wavelength; 1189 1190 // pick a random phi angle, and see if it lands on the detector 1191 // since the scattering is isotropic, I can safely pick a new, random value 1192 // this would not be true if simulating anisotropic scattering. 1193 testPhi = ran1a(&seed)*2*pi; 1194 1195 // is it on the detector? 1196 FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 1197 1198 if(xPixel != 1 && yPixel != 1) { 1199 isOn += 1; 1200 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 1201 indices[0] = xPixel; 1202 indices[1] = yPixel; 1203 if (retVal = MDGetNumericWavePointValue(wavH, indices, value)) 1204 return retVal; 1205 value[0] += 1; // Real part 1206 if (retVal = MDSetNumericWavePointValue(wavH, indices, value)) 1207 return retVal; 1208 //if(index==1) // only the single scattering events 1209 //dp = dp0 + xPixel + yPixel*numColumns; //offset the pointer to the exact memory location 1210 //*dp += 1; //increment the value there 1211 //endif 1212 } 1213 1214 1215 /* is this causing me a problem since I'm not locking these? Probably not, since it crashes even if I comment these out... */ 1216 if(theta_z < theta_max) { 1217 //Choose index for scattering angle array. 1218 //IND = NINT(THETA_z/DTH + 0.4999999) 1219 ind = (long)(theta_z/dth + 0.4999999); //round is eqivalent to nint() 1220 nt[ind] += 1; //Increment bin for angle. 1221 //Increment angle array for single scattering events. 1222 if (index == 1) { 1223 j1[ind] += 1; 1224 } 1225 //Increment angle array for double scattering events. 1226 if (index == 2) { 1227 j2[ind] += 1; 1228 } 1229 } 1230 /**/ 1231 1232 // increment all of the counters now since done==1 here and I'm sure to exit and get another neutron 1233 NScatterEvents += index; //total number of scattering events 1234 if(index == 1 && incoherentEvent == 1) { 1235 NSingleIncoherent += 1; 1236 } 1237 if(index == 1 && coherentEvent == 1) { 1238 NSingleCoherent += 1; 1239 } 1240 if(index == 2 && coherentEvent == 1 && incoherentEvent == 0) { 1241 NDoubleCoherent += 1; 1242 } 1243 if(index > 1) { 1244 NMultipleScatter += 1; 1245 } 1246 if(coherentEvent >= 1 && incoherentEvent == 0) { 1247 NCoherentEvents += 1; 1248 } 1249 if(coherentEvent > 1 && incoherentEvent == 0) { 1250 NMultipleCoherent += 1; 1251 } 1252 1253 } else { // index was zero, neutron must be transmitted, so just increment the proper counters and data 1254 isOn += 1; 1255 nt[0] += 1; 1256 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 1257 //indices[0] = xCtr_long; //don't put everything in one pixel 1258 //indices[1] = yCtr_long; 1259 indices[0] = (long)(xCtr+xx/pixSize+0.5); 1260 indices[1] = (long)(yCtr+yy/pixSize+0.5); 1261 // check for valid indices  got an XOP error, probably from here 1262 if(indices[0] > 127) indices[0] = 127; 1263 if(indices[0] < 0) indices[0] = 0; 1264 if(indices[1] > 127) indices[1] = 127; 1265 if(indices[1] < 0) indices[1] = 0; 1266 1267 if (retVal = MDGetNumericWavePointValue(wavH, indices, value)) 1268 return retVal; 1269 value[0] += 1; // Real part 1270 if (retVal = MDSetNumericWavePointValue(wavH, indices, value)) 1271 return retVal; 1272 } 1273 } 1274 } while (!done); 1275 } while(n1 < imon); 1276 1277 // assign the results to the wave 1278 1279 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 1280 value[0] = (double)n1; 1281 indices[0] = 0; 1282 if (retVal = MDSetNumericWavePointValue(p>resultsH, indices, value)) 1283 return retVal; 1284 value[0] = (double)n2; 1285 indices[0] = 1; 1286 if (retVal = MDSetNumericWavePointValue(p>resultsH, indices, value)) 1287 return retVal; 1288 value[0] = (double)isOn; 1289 indices[0] = 2; 1290 if (retVal = MDSetNumericWavePointValue(p>resultsH, indices, value)) 1291 return retVal; 1292 value[0] = (double)NScatterEvents; 1293 indices[0] = 3; 1294 if (retVal = MDSetNumericWavePointValue(p>resultsH, indices, value)) 1295 return retVal; 1296 value[0] = (double)NSingleCoherent; 1297 indices[0] = 4; 1298 if (retVal = MDSetNumericWavePointValue(p>resultsH, indices, value)) 1299 return retVal; 1300 value[0] = (double)NMultipleCoherent; 1301 indices[0] = 5; 1302 if (retVal = MDSetNumericWavePointValue(p>resultsH, indices, value)) 1303 return retVal; 1304 value[0] = (double)NMultipleScatter; 1305 indices[0] = 6; 1306 if (retVal = MDSetNumericWavePointValue(p>resultsH, indices, value)) 1307 return retVal; 1308 value[0] = (double)NCoherentEvents; 1309 indices[0] = 7; 1310 if (retVal = MDSetNumericWavePointValue(p>resultsH, indices, value)) 1311 return retVal; 1312 1313 // HSetState((Handle)wavH, hState); //release the handle of the 2D data wave 1314 // WaveHandleModified(wavH); // Inform Igor that we have changed the wave. (CALLBACK! needed, but not allowed in Threading) 1315 1316 return(0); 1317 } 1318 //////// end of X4 1319 1320 1321 468
Note: See TracChangeset
for help on using the changeset viewer.