# File Util/rebin.x
# February 12, 2001
# Adapted by Doug Mink to read MWCS data
# Adapted by Stephen Levine , University of Wisconsin (madraf::levine)
# from AL_REBIN in T_REBIN (in onedspec)

# REBIN rebins a spectrum to a specified range in wavelength
# REBINL rebins a spectrum to a specified range in log wavelength

include <mach.h>
include <ctype.h>
include <error.h>
include <imhdr.h>
include <fset.h>
include <math/iminterp.h>
include <math/curfit.h>
include <smw.h>

#rebinning defs	
define RB_NEAREST	1
define RB_LINEAR	2
define RB_SPLINE3	3
define RB_POLY3		4
define RB_POLY5		5
define RB_SUMS		6
define RB_FUNCTIONS	"|nearest|linear|spline3|poly3|poly5|sums|"

procedure rebinl (imin, sh, indlog, imout, outpix, outw0, outdw)

real	imin[ARB]	# Input spectrum
pointer	sh		# Spectrum header structure
double	indlog		# Log wavelength shift for input spectrum
real	imout[ARB]	# Output spectrum
int	outpix		# Number of pixels in output spectrum
double	outw0		# Output starting log wavelength
double	outdw		# Output log wavelength increment

pointer	invert, sp
char	interp_mode[SZ_LINE]
int	inpix		# Number of pixels in input spectrum
double	lw		# Log wavelength
double	wl		# Wavelength in Angstroms
double	px		# Pixel number (1=center of first)
int	mode, user_mode
int	i
int	clgwrd()
double	wcs_l2p()
#int	fd, open()
#double	inw0		# Input starting wavelength
#double	inw1		# Input final wavelength

begin

# Get rebinning method
	user_mode = clgwrd ("interp_mode", interp_mode, SZ_LINE, RB_FUNCTIONS)
	call cfit_mode (user_mode,mode)

# Make room for inverted solution
	call smark (sp)
	call salloc (invert, outpix+1, TY_DOUBLE)

#	fd = open ("rebin.log", APPEND, TEXT_FILE)
#	call fprintf (fd,"\n%s wavelengths %11.5f - %11.5f by %11.5f %d\n")
#	    call pargstr (SPECTRUM(sh))
#	    call pargd (W0(sh))
#	    call pargd (W1(sh))
#	    call pargd (WP(sh))
#	    call pargi (DC(sh))
#	inw0 = log10 (W0(sh)) + indlog
#	inw1 = log10 (W1(sh)) + indlog
#	call fprintf (fd,"%s log wavelengths %11.9f - %11.9f shifted by %11.9f\n")
#	    call pargstr (TITLE(sh))
#	    call pargd (inw0)
#	    call pargd (inw1)
#	    call pargd (indlog)

# Compute pixel position as a function of log lambda.
	call wcs_set (sh)
	do i = 1, outpix {
	    lw = outw0 + (outdw * double (i - 1)) - indlog
	    wl = 10.d0 ** lw
	    px = wcs_l2p (lw)
	    Memd[invert+i-1] = px
#	    call fprintf (fd," %d wl=%11.5f pixel=%11.5f\n")
#		call pargi (i)
#		call pargd (wl)
#		call pargd (px)
	    }
#	call close (fd)

# Interpolate input spectrum to output log lambda spectrum
	inpix = SN(sh)
	if (mode == RB_SUMS)
	    call resum (imin, imout, Memd[invert], inpix, outpix)
	else
	    call reinterp (imin, imout, Memd[invert], inpix, outpix, mode)

	call sfree (sp)

end


procedure rebin (imin, sh, indlog, imout, outpix, outw0, outdw)

real	imin[ARB]	# Input spectrum
pointer	sh		# Spectrum header structure
double	indlog		# Log wavelength shift for input spectrum
real	imout[ARB]	# Output spectrum
int	outpix		# Number of pixels in output spectrum
double	outw0		# Output starting wavelength
double	outdw		# Output wavelength increment

pointer	invert, sp
char	interp_mode[SZ_LINE]
int	inpix		# Number of pixels in input spectrum
double	lw		# Log wavelength
double	wl,wl0		# Wavelength in Angstroms
double	px		# Pixel number (1=center of first)
int	mode, user_mode
int	i
int	clgwrd()
double	wcs_w2p()

begin

# Get rebinning method
	user_mode = clgwrd ("interp_mode", interp_mode, SZ_LINE, RB_FUNCTIONS)
	call cfit_mode (user_mode,mode)

# Make room for inverted solution
	call smark (sp)
	call salloc (invert, outpix+1, TY_DOUBLE)

# Compute wavelength at each pixel position
	call wcs_set (sh)
	wl0 = 0.d0
	do i = 1, outpix {
	    if (indlog != 0.d0) {
		wl0 = outw0 + (outdw * double (i - 1))
		lw = dlog10 (wl0) - indlog
		wl = 10.d0 ** lw
		}
	    else
		wl = outw0 + (outdw * double (i - 1))
	    px = wcs_w2p (wl)
	    Memd[invert+i-1] = px
#	    call printf ("%d: %.4fA -> %.4fA = %.4f: %.4f\n")
#		call pargi (i)
#		call pargd (wl0)
#		call pargd (wl)
#		call pargd (lw)
#		call pargd (px)
	    }

# Interpolate input spectrum to output spectrum
	inpix = SN(sh)
	if (mode == RB_SUMS)
	    call resum (imin, imout, Memd[invert], inpix, outpix)
	else
	    call reinterp (imin, imout, Memd[invert], inpix, outpix, mode)

	call sfree (sp)

end


# CFIT_MODE -- Transform users input mode to CURFIT package mode.

procedure cfit_mode (user_mode, interp_mode)

int	user_mode, interp_mode

begin
	switch (user_mode) {
	    case RB_NEAREST:
		interp_mode = II_NEAREST

	    case RB_LINEAR:
		interp_mode = II_LINEAR

	    case RB_SPLINE3:
		interp_mode = II_SPLINE3

	    case RB_POLY3:
		interp_mode = II_POLY3

	    case RB_POLY5:
		interp_mode = II_POLY5

	    case RB_SUMS:
		interp_mode = RB_SUMS
	    }
end

# 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

# PIXEL_PARTS -- Integrate over partial pixels to obtain total flux
# 	         over specified region.

real procedure pixel_parts (y, x1, x2)

real	y[ARB]		# Spectrum vector
double	x1		# Starting fractional pixel
double	x2		# Ending fractional pixel

int 	i, i1, i2
double	cx1, cx2, frac1, frac2, sum

begin
	# Remember that pixel centers occur at intgral values
	# so a pixel extends from i-0.5 to i+0.5

	if (x1 >= x2 || x2 <= x1)
	    return (0.d0)
	cx1 = x1 - 0.5
	cx2 = x2 - 0.5

	i1 = int (cx1) + 1
	i2 = int (cx2) + 1

	if (i1 == i2) {
	    frac1 = x2 - x1
	    frac2 = 0.0
	    }
	else {
	    frac1 = int (cx1) + 1.0 - cx1
	    frac2 = cx2 - int (cx2)
	    }

	sum = frac1 * y[i1] + frac2 * y[i2]

	# Include inclusive whole pixels
	do i = i1+1,i2-1
	    sum = sum + y[i]

	return (sum)
end

# REINTERP -- Rebin the vector by interpolation.

procedure reinterp (pixin, pixout, invert, npin, npout, mode)

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 	mode		# Interpoloation mode

int	j
real	xpos

real	arieval()

begin
	do j = 1, npout {
	    xpos = invert[j]
	    if (xpos < 1.0 || xpos > real (npin))
		pixout[j] = 0.0
	    else
		pixout[j] = arieval (xpos, pixin, npin, mode)
	    }
end

# Jun  4 1993	Transform wavelength to pixel in subroutine
# Jun 15 1993	Clarify subroutines; add nearest value option
# Jun 16 1993	Change computations to double
# Jun 30 1993	Add spectrum header structure argument
# Jun 30 1993	Perform log-wavelength shift (for templates)
# Jul  1 1993	Keep subroutine arguments same type for Decstations
# Oct  6 1995	Change SHDR_* calls to SPHD_*

# Aug  7 1996	Use WCS_L2P instead of SHDR_PL

# May 16 1997	Add REBIN to rebin to linear wavelength; change file name
# Jul 17 1997	Fix REBIN linear wavelength bug: log10, not log
# Oct  6 1997	Use smw.h from IRAF, not local version

# Feb 12 2001	Fix SUMS interpolation option so it works with -dispersion