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.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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.