# 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)                //log scale cannot start at 0, so let's pick something close to what user wanted...
Wx = Wx/2
endif
CorrectStart = Wx
if(MinStep>0)
StartX = FindCorrectLogScaleStart(Wx,Wx[numpnts(Wx)-1],NumberOfPoints,MinStep)
else
StartX = CorrectStart
endif
Endx = StartX +abs(Wx[numpnts(Wx)-1] - Wx)
isGrowing = (Wx < 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
tempNewLogDist += CorrectStart - StartX

tempNewLogDistBinWidth[1,numpnts(tempNewLogDist)-2] = tempNewLogDist[p+1] - tempNewLogDist[p-1]
tempNewLogDistBinWidth = tempNewLogDistBinWidth
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)-log(X1))/w) - 10^(log(X1))
return abs(LastMinStep-w)
End

//***************************************************************************************************************** Forum Support Gallery