Changeset 623 for sans/XOP_Dev/MonteCarlo/MonteCarlo2.c
- Timestamp:
- Jan 28, 2010 6:09:34 PM (13 years ago)
- File:
-
- 1 edited
-
sans/XOP_Dev/MonteCarlo/MonteCarlo2.c (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
sans/XOP_Dev/MonteCarlo/MonteCarlo2.c
r591 r623 14 14 static int gCallSpinProcess = 1; // Set to 1 to all user abort (cmd dot) and background processing. 15 15 16 ////////// 17 // PROGRAM Monte_SANS 18 // PROGRAM simulates multiple SANS. 19 // revised 2/12/99 JGB 20 // added calculation of random deviate, and 2D 10/2008 SRK 21 22 // N1 = NUMBER OF INCIDENT NEUTRONS. 23 // N2 = NUMBER INTERACTED IN THE SAMPLE. 24 // N3 = NUMBER ABSORBED. 25 // THETA = SCATTERING ANGLE. 26 27 // works fine in the single-threaded case. 16 // these versions are DIRECT COPIES of the main version in MonteCarlo.c 17 // make changes there and copy them here. All that changes here is that the random 18 // number calls are different. 28 19 // 29 /// currently crashes if threaded. apparently something here is either doing an unknown callback, or is accessing30 // a bad place in memory.31 //32 // the operations SpinProcess() and WaveHandleModified() are callbacks and MUST be commented out before threading.33 // - some locks are non-existent34 // - supposedly safe wave access routines are used35 //36 // random number generators are not thread-safe, and can give less than random results, but is this enough to crash?37 // -- a possible workaround is to define multiple versions (poor man's threading)38 //39 //40 //41 42 20 // version X uses ran3 43 21 // version X2 uses ran1 22 // version X3 uses ran3a 23 // version X4 usus ran1a 44 24 45 25 int … … 469 449 return(0); 470 450 } 471 //////// end of main function for calculating multiple scattering 472 451 //////// end of X2 452 453 454 455 //////////////// X3 using ran3a 456 int 457 Monte_SANSX3(MC_ParamsPtr p) { 458 double *inputWave; /* pointer to double precision wave data */ 459 double *ran_dev; /* pointer to double precision wave data */ 460 double *nt; /* pointer to double precision wave data */ 461 double *j1; /* pointer to double precision wave data */ 462 double *j2; /* pointer to double precision wave data */ 463 double *nn; /* pointer to double precision wave data */ 464 // double *MC_linear_data; /* pointer to double precision wave data */ 465 double *results; /* pointer to double precision wave data */ 466 double retVal; //return value 467 468 long imon; 469 double r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas; 470 long ind,index,n_index; 471 double qmax,theta_max,q0,zpow; 472 long n1,n2,n3; 473 double dth,zz,xx,yy,phi; 474 double theta,ran,ll,rr,ttot; 475 long done,find_theta,err; //used as logicals 476 long xPixel,yPixel; 477 double vx,vy,vz,theta_z; 478 double sig_abs,ratio,sig_total; 479 double testQ,testPhi,left,delta,dummy,pi; 480 double sigabs_0,num_bins; 481 long NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent; 482 long NDoubleCoherent,NMultipleScatter,isOn,xCtr_long,yCtr_long; 483 long NMultipleCoherent,NCoherentEvents; 484 485 486 // for accessing the 2D wave data, direct method (see the WaveAccess example XOP) 487 waveHndl wavH; 488 int waveType,hState; 489 long numDimensions; 490 long dimensionSizes[MAX_DIMENSIONS+1]; 491 char* dataStartPtr; 492 long dataOffset; 493 long numRows, numColumns,numRows_ran_dev; 494 double *dp0, *dp, value[2]; // Pointers used for double data. 495 long seed; 496 long indices[MAX_DIMENSIONS]; 497 498 char buf[256]; 499 500 /* check that wave handles are all valid */ 501 if (p->inputWaveH == NIL) { 502 SetNaN64(&p->retVal); /* return NaN if wave is not valid */ 503 return(NON_EXISTENT_WAVE); 504 } 505 if (p->ran_devH == NIL) { 506 SetNaN64(&p->retVal); /* return NaN if wave is not valid */ 507 return(NON_EXISTENT_WAVE); 508 } 509 if (p->ntH == NIL) { 510 SetNaN64(&p->retVal); /* return NaN if wave is not valid */ 511 return(NON_EXISTENT_WAVE); 512 } 513 if (p->j1H == NIL) { 514 SetNaN64(&p->retVal); /* return NaN if wave is not valid */ 515 return(NON_EXISTENT_WAVE); 516 } 517 if (p->j2H == NIL) { 518 SetNaN64(&p->retVal); /* return NaN if wave is not valid */ 519 return(NON_EXISTENT_WAVE); 520 } 521 if (p->nnH == NIL) { 522 SetNaN64(&p->retVal); /* return NaN if wave is not valid */ 523 return(NON_EXISTENT_WAVE); 524 } 525 if (p->MC_linear_dataH == NIL) { 526 SetNaN64(&p->retVal); /* return NaN if wave is not valid */ 527 return(NON_EXISTENT_WAVE); 528 } 529 if (p->resultsH == NIL) { 530 SetNaN64(&p->retVal); /* return NaN if wave is not valid */ 531 return(NON_EXISTENT_WAVE); 532 } 533 534 p->retVal = 0; 535 536 // trusting that all inputs are DOUBLE PRECISION WAVES!!! 537 inputWave = WaveData(p->inputWaveH); 538 ran_dev = WaveData(p->ran_devH); 539 nt = WaveData(p->ntH); 540 j1 = WaveData(p->j1H); 541 j2 = WaveData(p->j2H); 542 nn = WaveData(p->nnH); 543 // MC_linear_data = WaveData(p->MC_linear_dataH); 544 results = WaveData(p->resultsH); 545 546 seed = (long)results[0]; 547 548 // sprintf(buf, "input seed = %ld\r", seed); 549 // XOPNotice(buf); 550 551 if(seed >= 0) { 552 seed = -1234509876; 553 } 554 555 dummy = ran3a(&seed); //initialize the random sequence by passing in a negative value 556 seed = 12348765; //non-negative after that does nothing 557 558 imon = (int)inputWave[0]; 559 r1 = inputWave[1]; 560 r2 = inputWave[2]; 561 xCtr = inputWave[3]; 562 yCtr = inputWave[4]; 563 sdd = inputWave[5]; 564 pixSize = inputWave[6]; 565 thick = inputWave[7]; 566 wavelength = inputWave[8]; 567 sig_incoh = inputWave[9]; 568 sig_sas = inputWave[10]; 569 xCtr_long = round(xCtr); 570 yCtr_long = round(yCtr); 571 572 dummy = MDGetWaveScaling(p->ran_devH, 0, &delta, &left); //0 is the rows 573 if (retVal = MDGetWaveDimensions(p->ran_devH, &numDimensions, dimensionSizes)) 574 return retVal; 575 numRows_ran_dev = dimensionSizes[0]; 576 577 pi = 4.0*atan(1.0); 578 579 // access the 2D wave data for writing using the direct method 580 wavH = p->MC_linear_dataH; 581 if (wavH == NIL) 582 return NOWAV; 583 584 // waveType = WaveType(wavH); 585 // if (waveType & NT_CMPLX) 586 // return NO_COMPLEX_WAVE; 587 // if (waveType==TEXT_WAVE_TYPE) 588 // return NUMERIC_ACCESS_ON_TEXT_WAVE; 589 // if (retVal = MDGetWaveDimensions(wavH, &numDimensions, dimensionSizes)) 590 // return retVal; 591 // numRows = dimensionSizes[0]; 592 // numColumns = dimensionSizes[1]; 593 594 // if (retVal = MDAccessNumericWaveData(wavH, kMDWaveAccessMode0, &dataOffset)) 595 // return retVal; 596 597 // hState = MoveLockHandle(wavH); // So wave data can't move. Remember to call HSetState when done. 598 // dataStartPtr = (char*)(*wavH) + dataOffset; 599 // dp0 = (double*)dataStartPtr; // Pointer to the start of the 2D wave data. 600 601 //scattering power and maximum qvalue to bin 602 // zpow = .1 //scattering power, calculated below 603 qmax = 4.0*pi/wavelength; //maximum Q to bin 1D data. (A-1) (not really used) 604 sigabs_0 = 0.0; // ignore absorption cross section/wavelength [1/(cm A)] 605 n_index = 50; // maximum number of scattering events per neutron 606 num_bins = 200; //number of 1-D bins (not really used) 607 608 //c total SAS cross-section 609 // 610 zpow = sig_sas*thick; //since I now calculate the sig_sas from the model 611 sig_abs = sigabs_0 * wavelength; 612 sig_total = sig_abs + sig_sas + sig_incoh; 613 // Print "The TOTAL XSECTION. (CM-1) is ",sig_total 614 // Print "The TOTAL SAS XSECTION. (CM-1) is ",sig_sas 615 // results[0] = sig_total; 616 // results[1] = sig_sas; 617 // RATIO = SIG_ABS / SIG_TOTAL 618 ratio = sig_incoh / sig_total; 619 620 theta_max = wavelength*qmax/(2*pi); 621 //C SET Theta-STEP SIZE. 622 dth = theta_max/num_bins; 623 // Print "theta bin size = dth = ",dth 624 625 //C INITIALIZE COUNTERS. 626 n1 = 0; 627 n2 = 0; 628 n3 = 0; 629 NSingleIncoherent = 0; 630 NSingleCoherent = 0; 631 NDoubleCoherent = 0; 632 NMultipleScatter = 0; 633 NScatterEvents = 0; 634 NMultipleCoherent = 0; 635 NCoherentEvents = 0; 636 637 isOn = 0; 638 639 //C MONITOR LOOP - looping over the number of incedent neutrons 640 //note that zz, is the z-position in the sample - NOT the scattering power 641 // NOW, start the loop, throwing neutrons at the sample. 642 do { 643 ////SpinProcess() IS A CALLBACK, and not good for Threading! 644 // if ((n1 % 1000 == 0) && gCallSpinProcess && SpinProcess()) { // Spins cursor and allows background processing. 645 // retVal = -1; // User aborted. 646 // break; 647 // } 648 649 vx = 0.0; // Initialize direction vector. 650 vy = 0.0; 651 vz = 1.0; 652 653 theta = 0.0; // Initialize scattering angle. 654 phi = 0.0; // Intialize azimuthal angle. 655 n1 += 1; // Increment total number neutrons counter. 656 done = 0; // True when neutron is absorbed or when scattered out of the sample. 657 index = 0; // Set counter for number of scattering events. 658 zz = 0.0; // Set entering dimension of sample. 659 incoherentEvent = 0; 660 coherentEvent = 0; 661 662 do { // Makes sure position is within circle. 663 ran = ran3a(&seed); //[0,1] 664 xx = 2.0*r1*(ran-0.5); //X beam position of neutron entering sample. 665 ran = ran3a(&seed); //[0,1] 666 yy = 2.0*r1*(ran-0.5); //Y beam position ... 667 rr = sqrt(xx*xx+yy*yy); //Radial position of neutron in incident beam. 668 } while(rr>r1); 669 670 do { //Scattering Loop, will exit when "done" == 1 671 // keep scattering multiple times until the neutron exits the sample 672 ran = ran3a(&seed); //[0,1] RANDOM NUMBER FOR DETERMINING PATH LENGTH 673 ll = path_len(ran,sig_total); 674 //Determine new scattering direction vector. 675 err = NewDirection(&vx,&vy,&vz,theta,phi); //vx,vy,vz updated, theta, phi unchanged by function 676 677 //X,Y,Z-POSITION OF SCATTERING EVENT. 678 xx += ll*vx; 679 yy += ll*vy; 680 zz += ll*vz; 681 rr = sqrt(xx*xx+yy*yy); //radial position of scattering event. 682 683 //sprintf(buf, "xx,yy,zz,vx,vy,vz,ll = %g %g %g %g %g %g %g\r",xx,yy,zz,vx,vy,vz,ll); 684 //XOPNotice(buf); 685 686 //Check whether interaction occurred within sample volume. 687 if (((zz > 0.0) && (zz < thick)) && (rr < r2)) { 688 //NEUTRON INTERACTED. 689 //sprintf(buf,"neutron interacted\r"); 690 //XOPNotice(buf); 691 692 index += 1; //Increment counter of scattering events. 693 if (index == 1) { 694 n2 += 1; //Increment # of scat. neutrons 695 } 696 ran = ran3a(&seed); //[0,1] 697 //Split neutron interactions into scattering and absorption events 698 if (ran > ratio ) { //C NEUTRON SCATTERED coherently 699 //sprintf(buf,"neutron scatters coherently\r"); 700 //XOPNotice(buf); 701 coherentEvent += 1; 702 find_theta = 0; //false 703 do { 704 // pick a q-value from the deviate function 705 // pnt2x truncates the point to an integer before returning the x 706 // so get it from the wave scaling instead 707 // q0 =left + binarysearchinterp(ran_dev,ran3a(seed))*delta; 708 709 q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran3a(&seed))*delta; 710 theta = q0/2/pi*wavelength; //SAS approximation 711 712 find_theta = 1; //always accept 713 714 //sprintf(buf, "after locate_interp call q0 = %g, theta = %g,left = %g,delta = %g\r",q0,theta,left,delta); 715 //XOPNotice(buf); 716 717 } while(!find_theta); 718 719 ran = ran3a(&seed); //[0,1] 720 phi = 2.0*pi*ran; //Chooses azimuthal scattering angle. 721 } else { 722 //NEUTRON scattered incoherently 723 //sprintf(buf,"neutron scatters incoherent\r"); 724 //XOPNotice(buf); 725 incoherentEvent += 1; 726 // phi and theta are random over the entire sphere of scattering 727 // !can't just choose random theta and phi, won't be random over sphere solid angle 728 729 ran = ran3a(&seed); //[0,1] 730 theta = acos(2.0*ran-1); 731 732 ran = ran3a(&seed); //[0,1] 733 phi = 2.0*pi*ran; //Chooses azimuthal scattering angle. 734 } //(ran > ratio) 735 } else { 736 //NEUTRON ESCAPES FROM SAMPLE -- bin it somewhere 737 done = 1; //done = true, will exit from loop 738 //Increment #scattering events array 739 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 740 indices[0] =index; //this sets access to nn[index] 741 if (index <= n_index) { 742 if (retVal = MDGetNumericWavePointValue(p->nnH, indices, value)) 743 return retVal; 744 value[0] += 1; // add one to the value 745 if (retVal = MDSetNumericWavePointValue(p->nnH, indices, value)) 746 return retVal; 747 // nn[index] += 1; 748 } 749 750 if( index != 0) { //neutron was scattered, figure out where it went 751 theta_z = acos(vz); // Angle (= 2theta) WITH respect to z axis. 752 testQ = 2*pi*sin(theta_z)/wavelength; 753 754 // pick a random phi angle, and see if it lands on the detector 755 // since the scattering is isotropic, I can safely pick a new, random value 756 // this would not be true if simulating anisotropic scattering. 757 testPhi = ran3a(&seed)*2*pi; 758 759 // is it on the detector? 760 FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 761 762 if(xPixel != -1 && yPixel != -1) { 763 isOn += 1; 764 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 765 indices[0] = xPixel; 766 indices[1] = yPixel; 767 if (retVal = MDGetNumericWavePointValue(wavH, indices, value)) 768 return retVal; 769 value[0] += 1; // Real part 770 if (retVal = MDSetNumericWavePointValue(wavH, indices, value)) 771 return retVal; 772 //if(index==1) // only the single scattering events 773 //dp = dp0 + xPixel + yPixel*numColumns; //offset the pointer to the exact memory location 774 //*dp += 1; //increment the value there 775 //endif 776 } 777 778 779 /* is this causing me a problem since I'm not locking these? Probably not, since it crashes even if I comment these out... */ 780 if(theta_z < theta_max) { 781 //Choose index for scattering angle array. 782 //IND = NINT(THETA_z/DTH + 0.4999999) 783 ind = round(theta_z/dth + 0.4999999); //round is eqivalent to nint() 784 nt[ind] += 1; //Increment bin for angle. 785 //Increment angle array for single scattering events. 786 if (index == 1) { 787 j1[ind] += 1; 788 } 789 //Increment angle array for double scattering events. 790 if (index == 2) { 791 j2[ind] += 1; 792 } 793 } 794 /**/ 795 796 // increment all of the counters now since done==1 here and I'm sure to exit and get another neutron 797 NScatterEvents += index; //total number of scattering events 798 if(index == 1 && incoherentEvent == 1) { 799 NSingleIncoherent += 1; 800 } 801 if(index == 1 && coherentEvent == 1) { 802 NSingleCoherent += 1; 803 } 804 if(index == 2 && coherentEvent == 1 && incoherentEvent == 0) { 805 NDoubleCoherent += 1; 806 } 807 if(index > 1) { 808 NMultipleScatter += 1; 809 } 810 if(coherentEvent >= 1 && incoherentEvent == 0) { 811 NCoherentEvents += 1; 812 } 813 if(coherentEvent > 1 && incoherentEvent == 0) { 814 NMultipleCoherent += 1; 815 } 816 817 } else { // index was zero, neutron must be transmitted, so just increment the proper counters and data 818 isOn += 1; 819 nt[0] += 1; 820 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 821 //indices[0] = xCtr_long; //don't put everything in one pixel 822 //indices[1] = yCtr_long; 823 indices[0] = (long)round(xCtr+xx/pixSize); 824 indices[1] = (long)round(yCtr+yy/pixSize); 825 // check for valid indices - got an XOP error, probably from here 826 if(indices[0] > 127) indices[0] = 127; 827 if(indices[0] < 0) indices[0] = 0; 828 if(indices[1] > 127) indices[1] = 127; 829 if(indices[1] < 0) indices[1] = 0; 830 831 if (retVal = MDGetNumericWavePointValue(wavH, indices, value)) 832 return retVal; 833 value[0] += 1; // Real part 834 if (retVal = MDSetNumericWavePointValue(wavH, indices, value)) 835 return retVal; 836 } 837 } 838 } while (!done); 839 } while(n1 < imon); 840 841 // assign the results to the wave 842 843 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 844 value[0] = (double)n1; 845 indices[0] = 0; 846 if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 847 return retVal; 848 value[0] = (double)n2; 849 indices[0] = 1; 850 if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 851 return retVal; 852 value[0] = (double)isOn; 853 indices[0] = 2; 854 if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 855 return retVal; 856 value[0] = (double)NScatterEvents; 857 indices[0] = 3; 858 if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 859 return retVal; 860 value[0] = (double)NSingleCoherent; 861 indices[0] = 4; 862 if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 863 return retVal; 864 value[0] = (double)NMultipleCoherent; 865 indices[0] = 5; 866 if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 867 return retVal; 868 value[0] = (double)NMultipleScatter; 869 indices[0] = 6; 870 if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 871 return retVal; 872 value[0] = (double)NCoherentEvents; 873 indices[0] = 7; 874 if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 875 return retVal; 876 877 // HSetState((Handle)wavH, hState); //release the handle of the 2D data wave 878 // WaveHandleModified(wavH); // Inform Igor that we have changed the wave. (CALLBACK! needed, but not allowed in Threading) 879 880 return(0); 881 } 882 //////// end of X3 883 884 885 ///////////////// X4 using ran1a() 886 int 887 Monte_SANSX4(MC_ParamsPtr p) { 888 double *inputWave; /* pointer to double precision wave data */ 889 double *ran_dev; /* pointer to double precision wave data */ 890 double *nt; /* pointer to double precision wave data */ 891 double *j1; /* pointer to double precision wave data */ 892 double *j2; /* pointer to double precision wave data */ 893 double *nn; /* pointer to double precision wave data */ 894 // double *MC_linear_data; /* pointer to double precision wave data */ 895 double *results; /* pointer to double precision wave data */ 896 double retVal; //return value 897 898 long imon; 899 double r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas; 900 long ind,index,n_index; 901 double qmax,theta_max,q0,zpow; 902 long n1,n2,n3; 903 double dth,zz,xx,yy,phi; 904 double theta,ran,ll,rr,ttot; 905 long done,find_theta,err; //used as logicals 906 long xPixel,yPixel; 907 double vx,vy,vz,theta_z; 908 double sig_abs,ratio,sig_total; 909 double testQ,testPhi,left,delta,dummy,pi; 910 double sigabs_0,num_bins; 911 long NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent; 912 long NDoubleCoherent,NMultipleScatter,isOn,xCtr_long,yCtr_long; 913 long NMultipleCoherent,NCoherentEvents; 914 915 916 // for accessing the 2D wave data, direct method (see the WaveAccess example XOP) 917 waveHndl wavH; 918 int waveType,hState; 919 long numDimensions; 920 long dimensionSizes[MAX_DIMENSIONS+1]; 921 char* dataStartPtr; 922 long dataOffset; 923 long numRows, numColumns,numRows_ran_dev; 924 double *dp0, *dp, value[2]; // Pointers used for double data. 925 long seed; 926 long indices[MAX_DIMENSIONS]; 927 928 char buf[256]; 929 930 /* check that wave handles are all valid */ 931 if (p->inputWaveH == NIL) { 932 SetNaN64(&p->retVal); /* return NaN if wave is not valid */ 933 return(NON_EXISTENT_WAVE); 934 } 935 if (p->ran_devH == NIL) { 936 SetNaN64(&p->retVal); /* return NaN if wave is not valid */ 937 return(NON_EXISTENT_WAVE); 938 } 939 if (p->ntH == NIL) { 940 SetNaN64(&p->retVal); /* return NaN if wave is not valid */ 941 return(NON_EXISTENT_WAVE); 942 } 943 if (p->j1H == NIL) { 944 SetNaN64(&p->retVal); /* return NaN if wave is not valid */ 945 return(NON_EXISTENT_WAVE); 946 } 947 if (p->j2H == NIL) { 948 SetNaN64(&p->retVal); /* return NaN if wave is not valid */ 949 return(NON_EXISTENT_WAVE); 950 } 951 if (p->nnH == NIL) { 952 SetNaN64(&p->retVal); /* return NaN if wave is not valid */ 953 return(NON_EXISTENT_WAVE); 954 } 955 if (p->MC_linear_dataH == NIL) { 956 SetNaN64(&p->retVal); /* return NaN if wave is not valid */ 957 return(NON_EXISTENT_WAVE); 958 } 959 if (p->resultsH == NIL) { 960 SetNaN64(&p->retVal); /* return NaN if wave is not valid */ 961 return(NON_EXISTENT_WAVE); 962 } 963 964 p->retVal = 0; 965 966 // trusting that all inputs are DOUBLE PRECISION WAVES!!! 967 inputWave = WaveData(p->inputWaveH); 968 ran_dev = WaveData(p->ran_devH); 969 nt = WaveData(p->ntH); 970 j1 = WaveData(p->j1H); 971 j2 = WaveData(p->j2H); 972 nn = WaveData(p->nnH); 973 // MC_linear_data = WaveData(p->MC_linear_dataH); 974 results = WaveData(p->resultsH); 975 976 seed = (long)results[0]; 977 978 // sprintf(buf, "input seed = %ld\r", seed); 979 // XOPNotice(buf); 980 981 if(seed >= 0) { 982 seed = -1234509876; 983 } 984 985 dummy = ran1a(&seed); //initialize the random sequence by passing in a negative value 986 seed = 12348765; //non-negative after that does nothing 987 988 imon = (int)inputWave[0]; 989 r1 = inputWave[1]; 990 r2 = inputWave[2]; 991 xCtr = inputWave[3]; 992 yCtr = inputWave[4]; 993 sdd = inputWave[5]; 994 pixSize = inputWave[6]; 995 thick = inputWave[7]; 996 wavelength = inputWave[8]; 997 sig_incoh = inputWave[9]; 998 sig_sas = inputWave[10]; 999 xCtr_long = round(xCtr); 1000 yCtr_long = round(yCtr); 1001 1002 dummy = MDGetWaveScaling(p->ran_devH, 0, &delta, &left); //0 is the rows 1003 if (retVal = MDGetWaveDimensions(p->ran_devH, &numDimensions, dimensionSizes)) 1004 return retVal; 1005 numRows_ran_dev = dimensionSizes[0]; 1006 1007 pi = 4.0*atan(1.0); 1008 1009 // access the 2D wave data for writing using the direct method 1010 wavH = p->MC_linear_dataH; 1011 if (wavH == NIL) 1012 return NOWAV; 1013 1014 // waveType = WaveType(wavH); 1015 // if (waveType & NT_CMPLX) 1016 // return NO_COMPLEX_WAVE; 1017 // if (waveType==TEXT_WAVE_TYPE) 1018 // return NUMERIC_ACCESS_ON_TEXT_WAVE; 1019 // if (retVal = MDGetWaveDimensions(wavH, &numDimensions, dimensionSizes)) 1020 // return retVal; 1021 // numRows = dimensionSizes[0]; 1022 // numColumns = dimensionSizes[1]; 1023 1024 // if (retVal = MDAccessNumericWaveData(wavH, kMDWaveAccessMode0, &dataOffset)) 1025 // return retVal; 1026 1027 // hState = MoveLockHandle(wavH); // So wave data can't move. Remember to call HSetState when done. 1028 // dataStartPtr = (char*)(*wavH) + dataOffset; 1029 // dp0 = (double*)dataStartPtr; // Pointer to the start of the 2D wave data. 1030 1031 //scattering power and maximum qvalue to bin 1032 // zpow = .1 //scattering power, calculated below 1033 qmax = 4.0*pi/wavelength; //maximum Q to bin 1D data. (A-1) (not really used) 1034 sigabs_0 = 0.0; // ignore absorption cross section/wavelength [1/(cm A)] 1035 n_index = 50; // maximum number of scattering events per neutron 1036 num_bins = 200; //number of 1-D bins (not really used) 1037 1038 //c total SAS cross-section 1039 // 1040 zpow = sig_sas*thick; //since I now calculate the sig_sas from the model 1041 sig_abs = sigabs_0 * wavelength; 1042 sig_total = sig_abs + sig_sas + sig_incoh; 1043 // Print "The TOTAL XSECTION. (CM-1) is ",sig_total 1044 // Print "The TOTAL SAS XSECTION. (CM-1) is ",sig_sas 1045 // results[0] = sig_total; 1046 // results[1] = sig_sas; 1047 // RATIO = SIG_ABS / SIG_TOTAL 1048 ratio = sig_incoh / sig_total; 1049 1050 theta_max = wavelength*qmax/(2*pi); 1051 //C SET Theta-STEP SIZE. 1052 dth = theta_max/num_bins; 1053 // Print "theta bin size = dth = ",dth 1054 1055 //C INITIALIZE COUNTERS. 1056 n1 = 0; 1057 n2 = 0; 1058 n3 = 0; 1059 NSingleIncoherent = 0; 1060 NSingleCoherent = 0; 1061 NDoubleCoherent = 0; 1062 NMultipleScatter = 0; 1063 NScatterEvents = 0; 1064 NMultipleCoherent = 0; 1065 NCoherentEvents = 0; 1066 1067 isOn = 0; 1068 1069 //C MONITOR LOOP - looping over the number of incedent neutrons 1070 //note that zz, is the z-position in the sample - NOT the scattering power 1071 // NOW, start the loop, throwing neutrons at the sample. 1072 do { 1073 ////SpinProcess() IS A CALLBACK, and not good for Threading! 1074 // if ((n1 % 1000 == 0) && gCallSpinProcess && SpinProcess()) { // Spins cursor and allows background processing. 1075 // retVal = -1; // User aborted. 1076 // break; 1077 // } 1078 1079 vx = 0.0; // Initialize direction vector. 1080 vy = 0.0; 1081 vz = 1.0; 1082 1083 theta = 0.0; // Initialize scattering angle. 1084 phi = 0.0; // Intialize azimuthal angle. 1085 n1 += 1; // Increment total number neutrons counter. 1086 done = 0; // True when neutron is absorbed or when scattered out of the sample. 1087 index = 0; // Set counter for number of scattering events. 1088 zz = 0.0; // Set entering dimension of sample. 1089 incoherentEvent = 0; 1090 coherentEvent = 0; 1091 1092 do { // Makes sure position is within circle. 1093 ran = ran1a(&seed); //[0,1] 1094 xx = 2.0*r1*(ran-0.5); //X beam position of neutron entering sample. 1095 ran = ran1a(&seed); //[0,1] 1096 yy = 2.0*r1*(ran-0.5); //Y beam position ... 1097 rr = sqrt(xx*xx+yy*yy); //Radial position of neutron in incident beam. 1098 } while(rr>r1); 1099 1100 do { //Scattering Loop, will exit when "done" == 1 1101 // keep scattering multiple times until the neutron exits the sample 1102 ran = ran1a(&seed); //[0,1] RANDOM NUMBER FOR DETERMINING PATH LENGTH 1103 ll = path_len(ran,sig_total); 1104 //Determine new scattering direction vector. 1105 err = NewDirection(&vx,&vy,&vz,theta,phi); //vx,vy,vz updated, theta, phi unchanged by function 1106 1107 //X,Y,Z-POSITION OF SCATTERING EVENT. 1108 xx += ll*vx; 1109 yy += ll*vy; 1110 zz += ll*vz; 1111 rr = sqrt(xx*xx+yy*yy); //radial position of scattering event. 1112 1113 //sprintf(buf, "xx,yy,zz,vx,vy,vz,ll = %g %g %g %g %g %g %g\r",xx,yy,zz,vx,vy,vz,ll); 1114 //XOPNotice(buf); 1115 1116 //Check whether interaction occurred within sample volume. 1117 if (((zz > 0.0) && (zz < thick)) && (rr < r2)) { 1118 //NEUTRON INTERACTED. 1119 //sprintf(buf,"neutron interacted\r"); 1120 //XOPNotice(buf); 1121 1122 index += 1; //Increment counter of scattering events. 1123 if (index == 1) { 1124 n2 += 1; //Increment # of scat. neutrons 1125 } 1126 ran = ran1a(&seed); //[0,1] 1127 //Split neutron interactions into scattering and absorption events 1128 if (ran > ratio ) { //C NEUTRON SCATTERED coherently 1129 //sprintf(buf,"neutron scatters coherently\r"); 1130 //XOPNotice(buf); 1131 coherentEvent += 1; 1132 find_theta = 0; //false 1133 do { 1134 // pick a q-value from the deviate function 1135 // pnt2x truncates the point to an integer before returning the x 1136 // so get it from the wave scaling instead 1137 // q0 =left + binarysearchinterp(ran_dev,ran1a(seed))*delta; 1138 1139 q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran1a(&seed))*delta; 1140 theta = q0/2/pi*wavelength; //SAS approximation 1141 1142 find_theta = 1; //always accept 1143 1144 //sprintf(buf, "after locate_interp call q0 = %g, theta = %g,left = %g,delta = %g\r",q0,theta,left,delta); 1145 //XOPNotice(buf); 1146 1147 } while(!find_theta); 1148 1149 ran = ran1a(&seed); //[0,1] 1150 phi = 2.0*pi*ran; //Chooses azimuthal scattering angle. 1151 } else { 1152 //NEUTRON scattered incoherently 1153 //sprintf(buf,"neutron scatters incoherent\r"); 1154 //XOPNotice(buf); 1155 incoherentEvent += 1; 1156 // phi and theta are random over the entire sphere of scattering 1157 // !can't just choose random theta and phi, won't be random over sphere solid angle 1158 1159 ran = ran1a(&seed); //[0,1] 1160 theta = acos(2.0*ran-1); 1161 1162 ran = ran1a(&seed); //[0,1] 1163 phi = 2.0*pi*ran; //Chooses azimuthal scattering angle. 1164 } //(ran > ratio) 1165 } else { 1166 //NEUTRON ESCAPES FROM SAMPLE -- bin it somewhere 1167 done = 1; //done = true, will exit from loop 1168 //Increment #scattering events array 1169 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 1170 indices[0] =index; //this sets access to nn[index] 1171 if (index <= n_index) { 1172 if (retVal = MDGetNumericWavePointValue(p->nnH, indices, value)) 1173 return retVal; 1174 value[0] += 1; // add one to the value 1175 if (retVal = MDSetNumericWavePointValue(p->nnH, indices, value)) 1176 return retVal; 1177 // nn[index] += 1; 1178 } 1179 1180 if( index != 0) { //neutron was scattered, figure out where it went 1181 theta_z = acos(vz); // Angle (= 2theta) WITH respect to z axis. 1182 testQ = 2*pi*sin(theta_z)/wavelength; 1183 1184 // pick a random phi angle, and see if it lands on the detector 1185 // since the scattering is isotropic, I can safely pick a new, random value 1186 // this would not be true if simulating anisotropic scattering. 1187 testPhi = ran1a(&seed)*2*pi; 1188 1189 // is it on the detector? 1190 FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel); 1191 1192 if(xPixel != -1 && yPixel != -1) { 1193 isOn += 1; 1194 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 1195 indices[0] = xPixel; 1196 indices[1] = yPixel; 1197 if (retVal = MDGetNumericWavePointValue(wavH, indices, value)) 1198 return retVal; 1199 value[0] += 1; // Real part 1200 if (retVal = MDSetNumericWavePointValue(wavH, indices, value)) 1201 return retVal; 1202 //if(index==1) // only the single scattering events 1203 //dp = dp0 + xPixel + yPixel*numColumns; //offset the pointer to the exact memory location 1204 //*dp += 1; //increment the value there 1205 //endif 1206 } 1207 1208 1209 /* is this causing me a problem since I'm not locking these? Probably not, since it crashes even if I comment these out... */ 1210 if(theta_z < theta_max) { 1211 //Choose index for scattering angle array. 1212 //IND = NINT(THETA_z/DTH + 0.4999999) 1213 ind = round(theta_z/dth + 0.4999999); //round is eqivalent to nint() 1214 nt[ind] += 1; //Increment bin for angle. 1215 //Increment angle array for single scattering events. 1216 if (index == 1) { 1217 j1[ind] += 1; 1218 } 1219 //Increment angle array for double scattering events. 1220 if (index == 2) { 1221 j2[ind] += 1; 1222 } 1223 } 1224 /**/ 1225 1226 // increment all of the counters now since done==1 here and I'm sure to exit and get another neutron 1227 NScatterEvents += index; //total number of scattering events 1228 if(index == 1 && incoherentEvent == 1) { 1229 NSingleIncoherent += 1; 1230 } 1231 if(index == 1 && coherentEvent == 1) { 1232 NSingleCoherent += 1; 1233 } 1234 if(index == 2 && coherentEvent == 1 && incoherentEvent == 0) { 1235 NDoubleCoherent += 1; 1236 } 1237 if(index > 1) { 1238 NMultipleScatter += 1; 1239 } 1240 if(coherentEvent >= 1 && incoherentEvent == 0) { 1241 NCoherentEvents += 1; 1242 } 1243 if(coherentEvent > 1 && incoherentEvent == 0) { 1244 NMultipleCoherent += 1; 1245 } 1246 1247 } else { // index was zero, neutron must be transmitted, so just increment the proper counters and data 1248 isOn += 1; 1249 nt[0] += 1; 1250 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 1251 //indices[0] = xCtr_long; //don't put everything in one pixel 1252 //indices[1] = yCtr_long; 1253 indices[0] = (long)round(xCtr+xx/pixSize); 1254 indices[1] = (long)round(yCtr+yy/pixSize); 1255 // check for valid indices - got an XOP error, probably from here 1256 if(indices[0] > 127) indices[0] = 127; 1257 if(indices[0] < 0) indices[0] = 0; 1258 if(indices[1] > 127) indices[1] = 127; 1259 if(indices[1] < 0) indices[1] = 0; 1260 1261 if (retVal = MDGetNumericWavePointValue(wavH, indices, value)) 1262 return retVal; 1263 value[0] += 1; // Real part 1264 if (retVal = MDSetNumericWavePointValue(wavH, indices, value)) 1265 return retVal; 1266 } 1267 } 1268 } while (!done); 1269 } while(n1 < imon); 1270 1271 // assign the results to the wave 1272 1273 MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 1274 value[0] = (double)n1; 1275 indices[0] = 0; 1276 if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 1277 return retVal; 1278 value[0] = (double)n2; 1279 indices[0] = 1; 1280 if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 1281 return retVal; 1282 value[0] = (double)isOn; 1283 indices[0] = 2; 1284 if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 1285 return retVal; 1286 value[0] = (double)NScatterEvents; 1287 indices[0] = 3; 1288 if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 1289 return retVal; 1290 value[0] = (double)NSingleCoherent; 1291 indices[0] = 4; 1292 if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 1293 return retVal; 1294 value[0] = (double)NMultipleCoherent; 1295 indices[0] = 5; 1296 if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 1297 return retVal; 1298 value[0] = (double)NMultipleScatter; 1299 indices[0] = 6; 1300 if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 1301 return retVal; 1302 value[0] = (double)NCoherentEvents; 1303 indices[0] = 7; 1304 if (retVal = MDSetNumericWavePointValue(p->resultsH, indices, value)) 1305 return retVal; 1306 1307 // HSetState((Handle)wavH, hState); //release the handle of the 2D data wave 1308 // WaveHandleModified(wavH); // Inform Igor that we have changed the wave. (CALLBACK! needed, but not allowed in Threading) 1309 1310 return(0); 1311 } 1312 //////// end of X4 1313 1314 1315
Note: See TracChangeset
for help on using the changeset viewer.
