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

//*****************************************************************************************************************