Changeset 1207


Ignore:
Timestamp:
Sep 6, 2019 2:29:24 PM (3 years ago)
Author:
krzywon
Message:

Add elliptical averaging as an option to the average parameters window. The method called currently performs an annular average.

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

Legend:

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

    r1165 r1207  
    3838        Variable/G root:myGlobals:Drawing:gDrawQCtr = 0 
    3939        Variable/G root:myGlobals:Drawing:gDrawQDelta = 1 
     40        Variable/G root:myGlobals:Drawing:gDrawRAxes = 1 
    4041         
    4142        //return to root 
     
    131132                        AnnularAverageTo1D(type) 
    132133                        break 
     134                case "Elliptical": 
     135                        EllipticalAverageTo1D(type) 
     136                        break 
    133137                case "Circular": 
    134138                case "Sector": 
    135139                        //circular or sector 
    136140                        CircularAverageTo1D(type)               //graph is drawn here 
     141                        break 
     142                case "2D_NXcanSAS": 
     143                        WriteNxCanSAS2D(type,"",1) 
    137144                        break 
    138145                case "2D ASCII": 
     
    156163                        case "2D ASCII": 
    157164                        case "QxQy ASCII": 
     165                        case "2D_NXcanSAS": 
    158166                                break 
    159167                        default: 
     
    225233                case "2D ASCII":                // execute if case matches expression 
    226234                case "QxQy ASCII": 
     235                case "2D_NXcanSAS": 
    227236                        String/G root:myGlobals:Drawing:gDrawInfoStr = ReplaceStringByKey("AVTYPE", tempStr, choice, "=", ";") 
    228237                        Button P_DoAvg,title="Save ASCII" 
     
    232241                case "Sector_PlusMinus": 
    233242                case "Rectangular": 
     243                case "Elliptical": 
    234244                case "Annular": 
    235245                        String/G root:myGlobals:Drawing:gDrawInfoStr = ReplaceStringByKey("AVTYPE", tempStr, choice, "=", ";") 
     
    253263        strswitch(choice)       // string switch 
    254264                case "2D ASCII": 
    255                 case "QxQy ASCII":                                       
     265                case "QxQy ASCII":       
     266                case "2D_NXcanSAS":      
    256267                case "Circular":                //disable everything for these three choices 
    257268                        SetVariable Phi_p,disable=yes 
     
    260271                        SetVariable DPhi_p,disable=yes 
    261272                        SetVariable width_p,disable=yes 
     273                        SetVariable RAxes_p,disable=yes 
    262274                        popupmenu sides,disable=yes 
    263275                        break 
     
    269281                        SetVariable DPhi_p,disable=no 
    270282                        SetVariable width_p,disable=yes 
     283                        SetVariable RAxes_p,disable=yes 
    271284                        popupmenu sides,disable=no 
    272285                        break 
     
    277290                        SetVariable DPhi_p,disable=yes 
    278291                        SetVariable width_p,disable=no 
     292                        SetVariable RAxes_p,disable=yes 
    279293                        popupmenu sides,disable=no 
    280294                        break 
     
    285299                        SetVariable DPhi_p,disable=yes 
    286300                        SetVariable width_p,disable=yes 
     301                        SetVariable RAxes_p,disable=yes 
    287302                        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 
    288312                        break 
    289313                default:                                                        // optional default expression executed 
     
    307331        PopupMenu av_choice,pos={61,7},size={144,20},proc=AvTypePopMenuProc,title="AverageType" 
    308332        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;\"" 
    310334        Button ave_help,pos={260,7},size={25,20},proc=ShowAvePanelHelp,title="?" 
    311335        Button ave_help,help={"Show the help file for averaging options"} 
     
    326350        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."} 
    327351        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 
    328355        SetVariable width_p,pos={15,155},size={115,15},proc=WidthSetVarProc,title="Width (pixels)" 
    329356        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  
    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// 
     673Function 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 
     879End 
Note: See TracChangeset for help on using the changeset viewer.