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