# RESUM -- Rebinning using a partial pixel summation technique to 
#          preserve total flux.

procedure resum (pixin, pixout, invert, npin, npout)

real	pixin[ARB]	# Input spectrum
real	pixout[ARB]	# Output spectrum (returned)
double	invert[ARB]	# Output to input pixel map
int 	npin		# Number of pixels in input spectrum
int 	npout		# Number of pixels in output spectrum

int 	i
double	x1, x2, xlo, xhi, dx, xmin, xmax

real	pixel_parts()

begin
	# Initialize
	xmin = 0.5d0
	xmax = double (npin) + 0.5d0
	x1 = invert [1]
	x2 = invert [2]
	dx = x2 - x1
	if (dx < 0) {
	    dx = x1 - x2
	    xlo = x2 - (dx / 2.0)
	    xhi = x1 + dx/2
	    }
	else {
	    xlo = x1 - (dx / 2.0)
	    xhi = x2 + dx/2
	    }
	xlo = max (xlo, xmin)
	xhi = min (xhi, xmax)
	pixout[1] = pixel_parts (pixin, xlo, xhi)

	do i = 2, npout {
	    x1 = x2
	    x2 = invert [i + 1]
	    dx = x2 - x1
	    if (dx < 0) {
		dx = x1 - x2
		xhi = xlo
		xlo = x2 - (dx / 2.0)
		}
	    else {
		xlo = xhi
		xhi = x2 + dx/2
		}
	    xlo = max (xlo, xmin)
	    xhi = min (xhi, xmax)

	    # Integrate over the partial pixels from xlo->xhi
	    pixout[i] = pixel_parts (pixin, xlo, xhi)
	    }
end