procedure fusebin (infile, fitsout, binsize)

file	infile	{"",prompt="Name of file with list of input FUSE files"}
file	fitsout	{"",prompt="Name of output summed flux file"}
int	binsize {"",prompt="Number of pixels to bin"}

begin

int	nbin, naxis2
file	fitsfile, outfile, wave1, flx1, err1, wave2, flx2, err2

# Get input parameters
    fitsfile = infile
    outfile = fitsout
    nbin = binsize

# Test for existence of output file
    if ( access(outfile) ) {
	print("Operation would overwrite existing file ", outfile)
	bye
    }

# Load necessary packages
    if (! defpac ("ttools")) {
	print("Please load the 'ttools' package before running this task.")
	bye
    }

# Create temporary files
    wave1 = mktemp ("tmp$wave1.")//".fits"
    flx1 = mktemp ("tmp$flx1.")//".fits"
    err1 = mktemp ("tmp$err1.")//".fits"
    wave2 = mktemp ("tmp$wave2.")//".fits"
    flx2 = mktemp ("tmp$flx2.")//".fits"
    err2 = mktemp ("tmp$err2.")//".fits"

# Get dimension of input table
    keypar (fitsfile//"[1]", "NAXIS2")
    naxis2 = int(keypar.value)

# Extract wavelength, flux and error arrays into IRAF images
    tabim(fitsfile, wave1, "WAVE", 1, naxis2)
    tabim(fitsfile, flx1, "FLUX", 1, naxis2)
    tabim(fitsfile, err1, "ERROR", 1, naxis2)

# Square error array for quadrature sum
    imarith(err1, "*", err1, err1)

# Bin the data
    blkavg(wave1, wave2, nbin, 1, 1, 1, 1, 1, 1, option="average")
    blkavg(flx1, flx2, nbin, 1, 1, 1, 1, 1, 1, option="average")
    blkavg(err1, err2, nbin, 1, 1, 1, 1, 1, 1, option="average")

# Convert summed variance into output error array
    imarith(err2, "/", nbin, err2)
    imfunc(err2, err2, "sqrt")

# Convert image files back into FUSE-format table
    imtab (wave2, outfile, "WAVE", pname="", formats="", tbltype="row")
    imtab (flx2, outfile, "FLUX", pname="", formats="", tbltype="row")
    imtab (err2, outfile, "ERROR", pname="", formats="", tbltype="row")

# Copy over header of original file
    imcopy(fitsfile//"[0]", outfile//"[0,over+]")

# Delete temporary images
    imdelete(wave1)
    imdelete(flx1)
    imdelete(err1)
    imdelete(wave2)
    imdelete(flx2)
    imdelete(err2)

end