Changeset 1207 for sans/Dev/branches/elliptical_averaging/NCNR_User_Procedures/Reduction/SANS/RectAnnulAvg.ipf
- Timestamp:
- Sep 6, 2019 2:29:24 PM (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sans/Dev/branches/elliptical_averaging/NCNR_User_Procedures/Reduction/SANS/RectAnnulAvg.ipf
r961 r1207 660 660 Return 0 661 661 End 662 663 //performs an elliptical average using isointensity contours, centered on a 664 //specified q-value (Intensity vs. angle) 665 //the parameters in the global keyword-string must have already been set somewhere 666 //either directly or from the protocol 667 // 668 //the input (data in the "type" folder) must be on linear scale - the calling routine is 669 //responsible for this 670 //averaged data is written to the data folder and plotted. data is not written 671 //to disk from this routine. 672 // 673 Function EllipticalAverageTo1D(type) 674 String type 675 676 SVAR keyListStr = root:myGlobals:Protocols:gAvgInfoStr 677 678 //type is the data type to do the averaging on, and will be set as the current folder 679 //get the current displayed data (so the correct folder is used) 680 String destPath = "root:Packages:NIST:"+type 681 682 Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,dtsize,dtdist 683 Variable rcentr,large_num,small_num,dtdis2,nq,xoffst,xbm,ybm,ii 684 Variable rc,delr,rlo,rhi,dphi,nphi,dr 685 Variable lambda,trans 686 Wave reals = $(destPath + ":RealsRead") 687 688 // center of detector, for non-linear corrections 689 NVAR pixelsX = root:myGlobals:gNPixelsX 690 NVAR pixelsY = root:myGlobals:gNPixelsY 691 692 xcenter = pixelsX/2 + 0.5 // == 64.5 for 128x128 Ordela 693 ycenter = pixelsY/2 + 0.5 // == 64.5 for 128x128 Ordela 694 695 // beam center, in pixels 696 x0 = reals[16] 697 y0 = reals[17] 698 //detector calibration constants 699 sx = reals[10] //mm/pixel (x) 700 sx3 = reals[11] //nonlinear coeff 701 sy = reals[13] //mm/pixel (y) 702 sy3 = reals[14] //nonlinear coeff 703 704 dtsize = 10*reals[20] //det size in mm 705 dtdist = 1000*reals[18] // det distance in mm 706 lambda = reals[26] 707 708 Variable qc = NumberByKey("QCENTER",keyListStr,"=",";") 709 Variable nw = NumberByKey("QDELTA",keyListStr,"=",";") 710 711 dr = 1 //minimum annulus width, keep this fixed at one 712 NVAR numPhiSteps = root:Packages:NIST:gNPhiSteps 713 nphi = numPhiSteps //number of anular sectors is set by users 714 715 rc = 2*dtdist*asin(qc*lambda/4/Pi) //in mm 716 delr = nw*sx/2 717 rlo = rc-delr 718 rhi = rc + delr 719 dphi = 360/nphi 720 721 /// data wave is data in the current folder which was set at the top of the function 722 Wave data=$(destPath + ":data") 723 //Check for the existence of the mask, if not, make one (local to this folder) that is null 724 725 if(WaveExists($"root:Packages:NIST:MSK:data") == 0) 726 Print "There is no mask file loaded (WaveExists)- the data is not masked" 727 Make/O/N=(pixelsX,pixelsY) $(destPath + ":mask") 728 WAVE mask = $(destPath + ":mask") 729 mask = 0 730 else 731 Wave mask=$"root:Packages:NIST:MSK:data" 732 Endif 733 734 rcentr = 150 //pixels within rcentr of beam center are broken into 9 parts 735 // values for error if unable to estimate value 736 //large_num = 1e10 737 large_num = 1 //1e10 value (typically sig of last data point) plots poorly, arb set to 1 738 small_num = 1e-10 739 740 // output wave are expected to exist (?) initialized to zero, what length? 741 // 300 points on VAX --- 742 Variable wavePts=500 743 Make/O/N=(wavePts) $(destPath + ":phival"),$(destPath + ":aveint") 744 Make/O/N=(wavePts) $(destPath + ":ncells"),$(destPath + ":sig"),$(destPath + ":sigave") 745 WAVE phival = $(destPath + ":phival") 746 WAVE aveint = $(destPath + ":aveint") 747 WAVE ncells = $(destPath + ":ncells") 748 WAVE sig = $(destPath + ":sig") 749 WAVE sigave = $(destPath + ":sigave") 750 751 phival = 0 752 aveint = 0 753 ncells = 0 754 sig = 0 755 sigave = 0 756 757 dtdis2 = dtdist^2 758 nq = 1 759 xoffst=0 760 //distance of beam center from detector center 761 xbm = FX(x0,sx3,xcenter,sx) 762 ybm = FY(y0,sy3,ycenter,sy) 763 764 //BEGIN AVERAGE ********** 765 Variable xi,xd,x,y,yd,yj,nd,fd,nd2,iphi,ntotal,var 766 Variable jj,data_pixel,xx,yy,ll,kk,rij,phiij,avesq,aveisq 767 768 // IGOR arrays are indexed from [0][0], FORTAN from (1,1) (and the detector too) 769 // loop index corresponds to FORTRAN (old code) 770 // and the IGOR array indices must be adjusted (-1) to the correct address 771 ntotal = 0 772 ii=1 773 do 774 xi = ii 775 xd = FX(xi,sx3,xcenter,sx) 776 x = xoffst + xd -xbm //x and y are in mm 777 778 jj = 1 779 do 780 data_pixel = data[ii-1][jj-1] //assign to local variable 781 yj = jj 782 yd = FY(yj,sy3,ycenter,sy) 783 y = yd - ybm 784 if(!(mask[ii-1][jj-1])) //masked pixels = 1, skip if masked (this way works...) 785 nd = 1 786 fd = 1 787 if( (abs(x) > rcentr) || (abs(y) > rcentr)) //break pixel into 9 equal parts 788 nd = 3 789 fd = 2 790 Endif 791 nd2 = nd^2 792 ll = 1 //"el-el" loop index 793 do 794 xx = x + (ll - fd)*sx/3 795 kk = 1 796 do 797 yy = y + (kk - fd)*sy/3 798 //test to see if center of pixel (i,j) lies in annulus 799 rij = sqrt(x*x + y*y)/dr + 1.001 800 //check whether pixel lies within width band 801 if((rij > rlo) && (rij < rhi)) 802 //in the annulus, do something 803 if (yy >= 0) 804 //phiij is in degrees 805 phiij = atan2(yy,xx)*180/Pi //0 to 180 deg 806 else 807 phiij = 360 + atan2(yy,xx)*180/Pi //180 to 360 deg 808 Endif 809 if (phiij > (360-0.5*dphi)) 810 phiij -= 360 811 Endif 812 iphi = trunc(phiij/dphi + 1.501) 813 aveint[iphi-1] += 9*data_pixel/nd2 814 sig[iphi-1] += 9*data_pixel*data_pixel/nd2 815 ncells[iphi-1] += 9/nd2 816 ntotal += 9/nd2 817 Endif //check if in annulus 818 kk+=1 819 while(kk<=nd) 820 ll += 1 821 while(ll<=nd) 822 Endif //masked pixel check 823 jj += 1 824 while (jj<=pixelsY) 825 ii += 1 826 while(ii<=pixelsX) //end of the averaging 827 828 //compute phi-values and errors 829 830 ntotal /=9 831 832 kk = 1 833 do 834 phival[kk-1] = dphi*(kk-1) 835 if(ncells[kk-1] != 0) 836 aveint[kk-1] = aveint[kk-1]/ncells[kk-1] 837 avesq = aveint[kk-1]*aveint[kk-1] 838 aveisq = sig[kk-1]/ncells[kk-1] 839 var = aveisq - avesq 840 if (var <=0 ) 841 sig[kk-1] = 0 842 sigave[kk-1] = 0 843 ncells[kk-1] /=9 844 else 845 if(ncells[kk-1] > 9) 846 sigave[kk-1] = sqrt(9*var/(ncells[kk-1]-9)) 847 sig[kk-1] = sqrt( abs(aveint[kk-1])/(ncells[kk-1]/9) ) 848 ncells[kk-1] /=9 849 else 850 sig[kk-1] = 0 851 sigave[kk-1] = 0 852 ncells[kk-1] /=9 853 Endif 854 Endif 855 Endif 856 kk+=1 857 while(kk<=nphi) 858 859 // data waves were defined as 200 points (=wavePts), but now have less than that (nphi) points 860 // use DeletePoints to remove junk from end of waves 861 Variable startElement,numElements 862 startElement = nphi 863 numElements = wavePts - startElement 864 DeletePoints startElement,numElements, phival,aveint,ncells,sig,sigave 865 866 //////////////end of VAX Phibin.for 867 868 //angle dependent transmission correction is not done in phiave 869 Ann_1D_Graph(aveint,phival,sigave) 870 871 //get rid of the default mask, if one was created (it is in the current folder) 872 //don't just kill "mask" since it might be pointing to the one in the MSK folder 873 Killwaves/z $(destPath+":mask") 874 875 //return to root folder (redundant) 876 SetDataFolder root: 877 878 Return 0 879 End
Note: See TracChangeset
for help on using the changeset viewer.
