Changeset 1212


Ignore:
Timestamp:
Sep 18, 2019 2:31:55 PM (3 years ago)
Author:
krzywon
Message:

Combine elliptical averaging with circular averaging and wrote the binning routine. Need to verify calculations and fix issues with elliptical drawing routine.

Location:
sans/Dev/branches/elliptical_averaging/NCNR_User_Procedures/Reduction/SANS
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • sans/Dev/branches/elliptical_averaging/NCNR_User_Procedures/Reduction/SANS/AvgGraphics.ipf

    r1208 r1212  
    3232        //ok, create the globals, fill the keyword string with all possible values (default) 
    3333        String/G root:myGlobals:Drawing:gDrawInfoStr = "AVTYPE=Circular;PHI=0;DPHI=0;WIDTH=0;SIDE=both;" 
    34                 root:myGlobals:Drawing:gDrawInfoStr += "QCENTER=0;QDELTA=0;" 
     34                root:myGlobals:Drawing:gDrawInfoStr += "QCENTER=0;QDELTA=0;RATIOAXES=1;" 
    3535        Variable/G root:myGlobals:Drawing:gDrawPhi =0 
    3636        Variable/G root:myGlobals:Drawing:gDrawWidth = 1 
     
    133133                        break 
    134134                case "Elliptical": 
    135                         EllipticalAverageTo1D(type) 
    136                         break 
    137135                case "Circular": 
    138136                case "Sector": 
    139                         //circular or sector 
     137                        //circular, elliptical, or sector 
    140138                        CircularAverageTo1D(type)               //graph is drawn here 
    141139                        break 
     
    309307                        SetVariable width_p,disable=yes 
    310308                        SetVariable RAxes_p,disable=no 
    311                         popupmenu sides,disable=no 
     309                        popupmenu sides,disable=yes 
    312310                        break 
    313311                default:                                                        // optional default expression executed 
     
    350348        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."} 
    351349        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" 
     350        SetVariable RAxes_p,pos={166,176},size={110,15},proc=AxisRatioSetVarProc,title="Axis Ratio" 
    353351        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."} 
    354352        SetVariable RAxes_p,limits={0,1,0.001},value= root:myGlobals:Drawing:gDrawRAxes 
     
    491489                variable RAxes = NumberByKey("RATIOAXES",drawStr,"=",";") 
    492490                 
     491                rr=0 
     492                gg=50000 
     493                bb=0 
     494                 
    493495                Make/O/N=(128,128) ellipsewave 
    494496                variable b = 64 
    495497                variable a = RAxes * b 
    496                 DrawAnEllipsoid("ellipsewave",x0,y0,a,b,phi) 
     498                DrawAnEllipsoid("ellipsewave",x0,y0,a,b,phi,rr,gg,bb) 
    497499          
    498500                Return 0                //exit the Draw routine 
     
    826828// OUTPUTS: 
    827829//   make all changes directly in xywave 
    828 Function DrawAnEllipsoid(wname,xo,yo,a,b,phi) 
     830Function DrawAnEllipsoid(wname,xo,yo,a,b,phi,rr,gg,bb) 
    829831    string wname 
    830     variable xo,yo,a,b,phi 
     832    variable xo,yo,a,b,phi,rr,gg,bb 
    831833    int ii 
    832834 
    833835    wave xywave = $wname 
    834836     
    835     make/FREE/n=361 urx,ury 
    836     make/n=361 xw,yw = 0 
     837    make/O/FREE/N=361 urx,ury 
     838    make/O/N=361 xw,yw = 0 
    837839    
    838840    variable cosalpha = cos(phi*Pi/180), sinalpha = sin(phi*Pi/180) 
     
    848850    xywave[xw,yw] = 1 
    849851     
    850     rr=0 
    851     gg=50000 
    852     bb=0 
    853      
    854852         //SetDrawEnv/W=SANS_Data xcoord= bottom,ycoord= left,linefgc= (rr,gg,bb),linethick= (thick),fillpat=0 
    855853         //AppendImage/W=SANS_Data ellipsewave 
     
    944942        ControlInfo/W=Average_Panel RAxes_p 
    945943        Variable val = V_Value 
    946         SVAR tempStr = root:myGlobals:Drawing:gDrawRAxes 
     944        SVAR tempStr = root:myGlobals:Drawing:gDrawInfoStr 
    947945        String newStr = ReplaceNumberByKey("RATIOAXES", tempStr, val, "=", ";") 
    948946         
    949         String/G root:myGlobals:Drawing:gDrawRAxes = newStr 
     947        String/G root:myGlobals:Drawing:gDrawInfoStr = newStr 
    950948         
    951949        //redraw the angles 
  • sans/Dev/branches/elliptical_averaging/NCNR_User_Procedures/Reduction/SANS/CircSectAve.ipf

    r902 r1212  
    3232         
    3333        SVAR keyListStr = root:myGlobals:Protocols:gAvgInfoStr          //this is the list that has it all 
    34         Variable isCircular = 0 
     34        Variable isCircular = 0, isElliptical = 0 
    3535         
    3636        if( cmpstr("Circular",StringByKey("AVTYPE",keyListStr,"=",";")) ==0) 
    3737                isCircular = 1          //set a switch for later 
     38        Endif 
     39        if( cmpstr("Elliptical",StringByKey("AVTYPE",keyListStr,"=",";")) ==0) 
     40                isElliptical = 1                //set a switch for later 
    3841        Endif 
    3942         
     
    4447        // 
    4548        Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,dtsize,dtdist,dr,ddr 
    46         Variable lambda,trans 
     49        Variable lambda,trans,raxes 
    4750        WAVE reals = $(destPath + ":RealsRead") 
    4851        WAVE/T textread = $(destPath + ":TextRead") 
     
    8891                phi_y = sin(phi_rad) 
    8992        Endif 
     93        if(isElliptical) 
     94                raxes = NumberByKey("RATIOAXES",keyListStr,"=",";") 
     95        EndIf 
    9096         
    9197        /// data wave is data in the current folder which was set at the top of the function 
     
    145151        //BEGIN AVERAGE ********** 
    146152        Variable xi,dxi,dx,jj,data_pixel,yj,dyj,dy,mask_val=0.1 
    147         Variable dr2,nd,fd,nd2,ll,kk,dxx,dyy,ir,dphi_p 
     153        Variable dr2,nd,fd,nd2,ll,kk,dxx,dyy,ir,dphi_p,rho 
    148154         
    149155        // IGOR arrays are indexed from [0][0], FORTAN from (1,1) (and the detector too) 
     
    178184                                        do 
    179185                                                dyy = dy + (kk - fd)*sy/3 
    180                                                 if(isCircular) 
    181                                                         //circular average, use all pixels 
    182                                                         //(increment)  
    183                                                         nq = IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2) 
     186                                                if(isCircular || isElliptical) 
     187                                                        //use all pixels 
     188                                                        if (isElliptical) 
     189                                                                rho = atan(dyy/dxx) - phi_rad 
     190                                                                nq = IncrementEllipticalPixel(data_pixel,ddr,dxx,dyy,rho,raxes, aveint,dsq,ncells,nq,nd2) 
     191                                                        else 
     192                                                                //circular average 
     193                                                                nq = IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2) 
     194                                                        EndIf 
    184195                                                else 
    185196                                                        //a sector average - determine azimuthal angle 
     
    396407End 
    397408 
     409//returns nq, new number of q-values 
     410//arrays aveint,dsq,ncells are also changed by this function 
     411// 
     412Function IncrementEllipticalPixel(dataPixel,ddr,dxx,dyy,rho,raxes,aveint,dsq,ncells,nq,nd2) 
     413        Variable dataPixel,ddr,dxx,dyy,rho,raxes 
     414        Wave aveint,dsq,ncells 
     415        Variable nq,nd2 
     416         
     417        Variable irCircular,ir 
     418         
     419        irCircular = (sqrt(dxx*dxx+dyy*dyy)/ddr) 
     420        ir = irCircular*sqrt(cos(rho)*cos(rho) + raxes*raxes*sin(rho)*sin(rho))+1 
     421        if (ir>nq) 
     422                nq = ir         //resets maximum number of q-values 
     423        endif 
     424        aveint[ir-1] += dataPixel/nd2           //ir-1 must be used, since ir is physical 
     425        dsq[ir-1] += dataPixel*dataPixel/nd2 
     426        ncells[ir-1] += 1/nd2 
     427         
     428        Return nq 
     429End 
     430 
    398431//function determines azimuthal angle dphi that a vector connecting 
    399432//center of detector to pixel makes with respect to vector 
  • sans/Dev/branches/elliptical_averaging/NCNR_User_Procedures/Reduction/SANS/RectAnnulAvg.ipf

    r1208 r1212  
    660660        Return 0 
    661661End 
    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 phi = NumberByKey("PHI",keyListStr,"=",";") 
    709         Variable ratio = NumberByKey("RATIOAXES",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.