- Timestamp:
- Sep 6, 2019 2:29:24 PM (3 years ago)
- Location:
- sans/Dev/branches/elliptical_averaging/NCNR_User_Procedures/Reduction/SANS
- Files:
-
- 2 edited
-
AvgGraphics.ipf (modified) (12 diffs)
-
RectAnnulAvg.ipf (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
sans/Dev/branches/elliptical_averaging/NCNR_User_Procedures/Reduction/SANS/AvgGraphics.ipf
r1165 r1207 38 38 Variable/G root:myGlobals:Drawing:gDrawQCtr = 0 39 39 Variable/G root:myGlobals:Drawing:gDrawQDelta = 1 40 Variable/G root:myGlobals:Drawing:gDrawRAxes = 1 40 41 41 42 //return to root … … 131 132 AnnularAverageTo1D(type) 132 133 break 134 case "Elliptical": 135 EllipticalAverageTo1D(type) 136 break 133 137 case "Circular": 134 138 case "Sector": 135 139 //circular or sector 136 140 CircularAverageTo1D(type) //graph is drawn here 141 break 142 case "2D_NXcanSAS": 143 WriteNxCanSAS2D(type,"",1) 137 144 break 138 145 case "2D ASCII": … … 156 163 case "2D ASCII": 157 164 case "QxQy ASCII": 165 case "2D_NXcanSAS": 158 166 break 159 167 default: … … 225 233 case "2D ASCII": // execute if case matches expression 226 234 case "QxQy ASCII": 235 case "2D_NXcanSAS": 227 236 String/G root:myGlobals:Drawing:gDrawInfoStr = ReplaceStringByKey("AVTYPE", tempStr, choice, "=", ";") 228 237 Button P_DoAvg,title="Save ASCII" … … 232 241 case "Sector_PlusMinus": 233 242 case "Rectangular": 243 case "Elliptical": 234 244 case "Annular": 235 245 String/G root:myGlobals:Drawing:gDrawInfoStr = ReplaceStringByKey("AVTYPE", tempStr, choice, "=", ";") … … 253 263 strswitch(choice) // string switch 254 264 case "2D ASCII": 255 case "QxQy ASCII": 265 case "QxQy ASCII": 266 case "2D_NXcanSAS": 256 267 case "Circular": //disable everything for these three choices 257 268 SetVariable Phi_p,disable=yes … … 260 271 SetVariable DPhi_p,disable=yes 261 272 SetVariable width_p,disable=yes 273 SetVariable RAxes_p,disable=yes 262 274 popupmenu sides,disable=yes 263 275 break … … 269 281 SetVariable DPhi_p,disable=no 270 282 SetVariable width_p,disable=yes 283 SetVariable RAxes_p,disable=yes 271 284 popupmenu sides,disable=no 272 285 break … … 277 290 SetVariable DPhi_p,disable=yes 278 291 SetVariable width_p,disable=no 292 SetVariable RAxes_p,disable=yes 279 293 popupmenu sides,disable=no 280 294 break … … 285 299 SetVariable DPhi_p,disable=yes 286 300 SetVariable width_p,disable=yes 301 SetVariable RAxes_p,disable=yes 287 302 popupmenu sides,disable=yes 303 break 304 case "Elliptical": 305 SetVariable Phi_p,disable=no 306 SetVariable Qctr_p,disable=yes 307 SetVariable Qdelta_p,disable=yes 308 SetVariable DPhi_p,disable=yes 309 SetVariable width_p,disable=yes 310 SetVariable RAxes_p,disable=no 311 popupmenu sides,disable=no 288 312 break 289 313 default: // optional default expression executed … … 307 331 PopupMenu av_choice,pos={61,7},size={144,20},proc=AvTypePopMenuProc,title="AverageType" 308 332 PopupMenu av_choice,help={"Select the type of average to perform, then make the required selections below and click \"DoAverage\" to plot the results"} 309 PopupMenu av_choice,mode=1,popvalue="Circular",value= #"\"Circular;Sector;Annular;Rectangular; 2D ASCII;QxQy ASCII;Sector_PlusMinus;\""333 PopupMenu av_choice,mode=1,popvalue="Circular",value= #"\"Circular;Sector;Annular;Rectangular;Elliptical;2D_NXcanSAS;2D ASCII;QxQy ASCII;Sector_PlusMinus;\"" 310 334 Button ave_help,pos={260,7},size={25,20},proc=ShowAvePanelHelp,title="?" 311 335 Button ave_help,help={"Show the help file for averaging options"} … … 326 350 SetVariable DPhi_p,help={"Enter the +/- range (0,45) of azimuthal angles to be included in the average. The bounding angles will be drawin in blue."} 327 351 SetVariable DPhi_p,limits={0,90,1},value= root:myGlobals:Drawing:gDrawDPhi 352 SetVariable RAxes_p,pos={166,176},size={110,15},proc=DeltaPhiSetVarProc,title="Axis Ratio" 353 SetVariable RAxes_p,help={"Enter the ratio of the minor to major axes for an isointensity contour. By definition, this value should be less than 1."} 354 SetVariable RAxes_p,limits={0,1,0.001},value= root:myGlobals:Drawing:gDrawRAxes 328 355 SetVariable width_p,pos={15,155},size={115,15},proc=WidthSetVarProc,title="Width (pixels)" 329 356 SetVariable width_p,help={"Enter the total width of the rectangular section in pixels. The bounding lines will be drawn in blue."} -
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.
