Rebin large set of x-y data on log-scale

//This routine will rebin approximately linearly distributed data into approximately logarithmically spaced data. 
//It will replace the input Wx and Wy with new Wx and Wy with the required NumberOfPoints
//Code will replace all other input data also, so make sure you pass in only copies of the data. 
//If MinStep > 0 it will try to set the values so the minimum step of the log scale is MinStep
//note: there are some limits on step size and number of points in and out. Basically, the step size needs to be 
//approximately equivalent the linear step. It cannot be much larger. And number of new points must be significantly smaller than input num points. 
//optional Wsdev and Wxsdev are standard deviation for each Wy (and Wx) value, they will be propagated through - sum(sdev^2)/numpnts for each bin. 
//optional Wxwidth will generate width of each new bin in x. NOTE: the edge is half linear distance between the two points, no log  
//skewing is done for edges. Therefore the width is really half of the distance between p-1 and p+1 points.  
//optional W1-W5 will be averaged for each bin, so this is way to propagate other data one may need to. 
//**********************************************************************************************************
Function RebinLogData(Wx,Wy,NumberOfPoints,MinStep,[Wsdev,Wxsdev,Wxwidth,W1, W2, W3, W4, W5])
        Wave Wx, Wy
        Variable NumberOfPoints, MinStep
        Wave Wsdev, Wxsdev
        Wave Wxwidth
        Wave W1, W2, W3, W4, W5
        variable CalcSdev, CalcxSdev, CalcWidth, CalcW1, CalcW2, CalcW3, CalcW4, CalcW5
        CalcSdev = ParamIsDefault(Wsdev) ?  0 : 1
        CalcxSdev = ParamIsDefault(Wxsdev) ?  0 : 1
        CalcWidth = ParamIsDefault(Wxwidth) ?  0 : 1
        CalcW1 = ParamIsDefault(W1) ?  0 : 1
        CalcW2 = ParamIsDefault(W2) ?  0 : 1
        CalcW3 = ParamIsDefault(W3) ?  0 : 1
        CalcW4 = ParamIsDefault(W4) ?  0 : 1
        CalcW5 = ParamIsDefault(W5) ?  0 : 1
        
        variable OldNumPnts=numpnts(Wx)
        if(3*NumberOfPoints>OldNumPnts)
            print "User requested rebinning of data, but old number of points is less than 3*requested number of points, no rebinning done"
            return 0
        endif
        variable StartX, EndX, iii, isGrowing, CorrectStart, logStartX, logEndX
        if(Wx[0]<=0)                //log scale cannot start at 0, so let's pick something close to what user wanted...  
            Wx[0] = Wx[1]/2
        endif
        CorrectStart = Wx[0]
        if(MinStep>0)
            StartX = FindCorrectLogScaleStart(Wx[0],Wx[numpnts(Wx)-1],NumberOfPoints,MinStep)
        else
            StartX = CorrectStart
        endif
        Endx = StartX +abs(Wx[numpnts(Wx)-1] - Wx[0])
        isGrowing = (Wx[0] < Wx[numpnts(Wx)-1]) ? 1 : 0
        make/O/D/FREE/N=(NumberOfPoints) tempNewLogDist, tempNewLogDistBinWidth
        logstartX=log(startX)
        logendX=log(endX)
        tempNewLogDist = logstartX + p*(logendX-logstartX)/numpnts(tempNewLogDist)
        tempNewLogDist = 10^(tempNewLogDist)
        startX = tempNewLogDist[0]
        tempNewLogDist += CorrectStart - StartX
    
        tempNewLogDistBinWidth[1,numpnts(tempNewLogDist)-2] = tempNewLogDist[p+1] - tempNewLogDist[p-1]
        tempNewLogDistBinWidth[0] = tempNewLogDistBinWidth[1]
        tempNewLogDistBinWidth[numpnts(tempNewLogDist)-1] = tempNewLogDistBinWidth[numpnts(tempNewLogDist)-2]
        make/O/D/FREE/N=(NumberOfPoints) Rebinned_WvX, Rebinned_WvY, Rebinned_Wv1, Rebinned_Wv2,Rebinned_Wv3, Rebinned_Wv4, Rebinned_Wv5, Rebinned_Wsdev, Rebinned_Wxsdev
        Rebinned_WvX=0
        Rebinned_WvY=0
        Rebinned_Wv1=0  
        Rebinned_Wv2=0  
        Rebinned_Wv3=0  
        Rebinned_Wv4=0  
        Rebinned_Wv5=0  
        Rebinned_Wsdev=0    
                Rebinned_Wxsdev=0
 
        variable i, j
        variable cntPoints, BinHighEdge
        //variable i will be from 0 to number of new points, moving through destination waves
        j=0     //this variable goes through data to be reduced, therefore it goes from 0 to numpnts(Wx)
        For(i=0;i<NumberOfPoints;i+=1)
            cntPoints=0
            BinHighEdge = tempNewLogDist[i]+tempNewLogDistBinWidth[i]/2
            if(isGrowing)
                Do
                    Rebinned_WvX[i]     += Wx[j]
                    Rebinned_WvY[i] += Wy[j]
                    if(CalcW1)
                        Rebinned_Wv1[i] += W1[j]
                    endif
                    if(CalcW2)
                        Rebinned_Wv2[i] += W2[j]
                    endif
                    if(CalcW3)
                        Rebinned_Wv3[i] += W3[j]
                    endif
                    if(CalcW4)
                        Rebinned_Wv4[i] += W4[j]
                    endif
                    if(CalcW5)
                        Rebinned_Wv5[i] += W5[j]
                    endif
                    if(CalcSdev)
                        Rebinned_Wsdev[i] += Wsdev[j]^2
                    endif
                    if(CalcxSdev)
                        Rebinned_Wxsdev[i] += Wxsdev[j]^2
                    endif
                    cntPoints+=1
                    j+=1
                While(Wx[j-1]<BinHighEdge && j<OldNumPnts)
            else
                Do
                    Rebinned_WvX[i]     += Wx[j]
                    Rebinned_WvY[i] += Wy[j]
                    if(CalcW1)
                        Rebinned_Wv1[i] += W1[j]
                    endif
                    if(CalcW2)
                        Rebinned_Wv2[i] += W2[j]
                    endif
                    if(CalcW3)
                        Rebinned_Wv3[i] += W3[j]
                    endif
                    if(CalcW4)
                        Rebinned_Wv4[i] += W4[j]
                    endif
                    if(CalcW5)
                        Rebinned_Wv5[i] += W5[j]
                    endif
                    if(CalcSdev)
                        Rebinned_Wsdev[i] += Wsdev[j]^2
                    endif
                    if(CalcxSdev)
                        Rebinned_Wxsdev[i] += Wxsdev[j]^2
                    endif
                    cntPoints+=1
                    j+=1
                While((Wx[j-1]>BinHighEdge) && (j<OldNumPnts))
            endif
            Rebinned_WvX[i]/=cntPoints   
            Rebinned_WvY[i]/=cntPoints
            if(CalcW1)
                Rebinned_Wv1[i]/=cntPoints
            endif
            if(CalcW2)
                Rebinned_Wv2[i]/=cntPoints
            endif
            if(CalcW3)
                Rebinned_Wv3[i]/=cntPoints
            endif
            if(CalcW4)
                Rebinned_Wv4[i]/=cntPoints
            endif
            if(CalcW5)
                Rebinned_Wv5[i]/=cntPoints
            endif
            if(CalcSdev)
                Rebinned_Wsdev[i]=sqrt(Rebinned_Wsdev[i])/cntPoints
            endif
            if(CalcxSdev)
                Rebinned_Wxsdev[i]=sqrt(Rebinned_Wxsdev[i])/cntPoints
            endif
        endfor
 
    Redimension/N=(numpnts(Rebinned_WvX))/D Wx, Wy
    Wx=Rebinned_WvX
    Wy=Rebinned_WvY
 
    if(CalcW1)
        Redimension/N=(numpnts(Rebinned_WvX))/D W1
        W1=Rebinned_Wv1
    endif
    if(CalcW2)
        Redimension/N=(numpnts(Rebinned_WvX))/D W2
        W2=Rebinned_Wv2
    endif
    if(CalcW3)
        Redimension/N=(numpnts(Rebinned_WvX))/D W3
        W3=Rebinned_Wv3
    endif
    if(CalcW4)
        Redimension/N=(numpnts(Rebinned_WvX))/D W4
        W4=Rebinned_Wv4
    endif
    if(CalcW5)
        Redimension/N=(numpnts(Rebinned_WvX))/D W5
        W5=Rebinned_Wv5
    endif
 
    if(CalcSdev)
        Redimension/N=(numpnts(Rebinned_WvX))/D Wsdev
        Wsdev = Rebinned_Wsdev
    endif
 
    if(CalcxSdev)
        Redimension/N=(numpnts(Rebinned_WvX))/D Wsdev
        Wxsdev = Rebinned_Wxsdev
    endif
    
    if(CalcWidth)
        Redimension/N=(numpnts(Rebinned_WvX))/D Wxwidth
        Wxwidth = tempNewLogDistBinWidth
    endif
end     
//**********************************************************************************************************
Function FindCorrectLogScaleStart(StartValue,EndValue,NumPoints,MinStep)
    variable StartValue,EndValue,NumPoints,MinStep
    Make/Free/N=3 w
    w={EndValue-StartValue, NumPoints,MinStep}
    Optimize /H=100/L=1e-5/I=100/T=(MinStep/10)/Q myFindStartValueFunc, w
        //the optimize parameters: 10^H is last reasonable number, 10^L is really close to 1, other are open for discussion.  
    //Test this works if you want to make sure...
        //  variable startX=log(V_minloc)
        //  variable endX=log(V_minloc+range)
        //  variable LastMinStep = 10^(startX + (endX-startX)/NumPoints) - 10^(startX)
        //  print LastMinStep
    return V_minloc
end
//**********************************************************************************************************
static Function myFindStartValueFunc(w,x1)     //static to keep this form polluting names others may use. 
    Wave w      //this is {totalRange, NumSteps,MinStep}
    Variable x1 //this is startValue where we need to start with log stepping...
    variable LastMinStep = 10^(log(X1) + (log(X1+w[0])-log(X1))/w[1]) - 10^(log(X1))
    return abs(LastMinStep-w[2])
End
 
//*****************************************************************************************************************