Basic Rebinning

#pragma rtGlobals=1		// Use modern global access method.

Function Pla_intRebin(x_init, y_init,s_init, x_rebin)
	Wave x_init, y_init, s_init,x_rebin

//Rebins histogrammed data into boundaries set by x_rebin.
//makes the waves W_rebin and W_RebinSD.

//it uses linear interpolation to work out the proportions of which original cells should be placed
//in the new bin boundaries.

// Precision will normally be lost.  
// It does not make sense to rebin to smaller bin boundaries.
	
	if(checkSorted(x_rebin) || checkSorted(x_init))
		print "The x_rebin and x_init must be monotonically increasing (Pla_intRebin)"
		return 1
	endif
	
	if(wavedims(X_init)!=1 || wavedims(y_init)!=1 || wavedims(s_init)!=1 || wavedims(X_rebin)!=1)
		print "All supplied waves must be 1D (Pla_intrebin)"	
		return 1
	endif
	
	if(numpnts(X_init)-1!= numpnts(y_init) || numpnts(y_init)!=numpnts(s_init))
		print "y_init and s_init must have one less point than x_init (Pla_intRebin)"
		return 1
	endif

	
	make/o/d/n =(numpnts(x_rebin)-1) W_rebin=0,W_RebinSD=0
	
	variable ii=0, kk = 0
	variable lowlim,upperlim
	
	for(ii=0; ii< numpnts(x_rebin)-1 ; ii+=1)

		//this gives the approximate position of where the new bin would start in the old bin		
		lowlim = binarysearchinterp(x_init,x_rebin[ii])
		upperlim = binarysearchinterp(x_init,x_rebin[ii+1])	
		
		//if your rebin x boundaries are set outisde those of the initial data then you won't get any counts.
		if(numtype(lowlim) && numtype(upperlim))
			W_rebin[ii] = 0
			W_RebinSD[ii] = 1
			continue
		endif
		
		//lower limit for your rebinned data may be outside the original histogram boundary
		//set it to the lowest point in this case
		if(numtype(lowlim) && numtype(upperlim) == 0)
			lowlim = 0
		endif
		
		//lower limit for rebinned boundary is in the original boundaries
		//but upperlimit has escaped, so set to the highest from the original data.
		if(numtype(lowlim)==0 && numtype(upperlim) )
			upperlim = numpnts(x_init)-1
		endif
		
		//now need to add the counts together
		
		//both upperlimit and lower limit rebin boundaries aren't the same unbinned cell
		//need to take a proportion of a lower and upper cell 
		if(trunc(lowlim) != trunc(upperlim))
			W_rebin[ii]  =  y_init[trunc(lowlim)]*(ceil(lowlim) - lowlim)
			W_rebin[ii] += y_init[trunc(upperlim)]*(upperlim - trunc(upperlim))
			
			W_RebinSD[ii]  = (s_init[trunc(lowlim)]*(ceil(lowlim) - lowlim))^2
			W_RebinSD[ii] += (s_init[trunc(upperlim)]*(upperlim - trunc(upperlim)))^2
			W_RebinSD[ii] = sqrt(W_RebinSD[ii])
		else
		//the upper and lower limits are in the same binned cell.  Need to work out
		//what proportion of the original cell is occupied by the difference between the limits
			W_rebin[ii] =	y_init[trunc(lowlim)]*(upperlim-lowlim)
			W_RebinSD[ii] =	s_init[trunc(lowlim)]*(upperlim-lowlim)
		endif
		
		//if the upper and lower limits span several of the original data, then you need to add counts 
		//from each individual cell.
		if((ceil(lowlim)<trunc(upperlim)) &&(trunc(upperlim) - ceil(lowlim) >=1))
			for(kk = 0 ; kk < trunc(upperlim) - ceil(lowlim) ;kk += 1)
				W_rebin[ii] += y_init[ceil(lowlim) + kk]
				W_RebinSD[ii] += (s_init[ceil(lowlim) + kk])^2
			endfor
			W_RebinSD[ii] = sqrt(W_RebinSD[ii])
		endif
	endfor
	
	return 0
End

Function checkSorted(aWave)
	Wave aWave
	variable ii
	for(ii=dimsize(awave,0)-1 ; ii>= 0 ; ii-=1)
		if(awave[ii] < awave[ii-1])
			return 1
		endif
	endfor
	return 0
End