procedure fusecombine (inlist, fitsout)

file	inlist	{"",prompt="Name of file with list of input FUSE files"}
file	fitsout	{"",prompt="Name of output summed flux file"}

begin

int	n, naxis2
real	texp, ttot
file	fitsfile, outfile, flxtmp, errtmp, flx, err, flxtab, outtmp
string	line

# Get input file names
    list = inlist
    outfile = fitsout

# 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
    flxtmp = mktemp ("tmp$flxtmp.")
    errtmp = mktemp ("tmp$errtmp.")
    outtmp = mktemp ("tmp$outtmp.")

# Loop through list of images
    n = 0
    while (fscan (list, line) != EOF) {
	fitsfile = line

	print(n," ",fitsfile)
	if ( n == 0 ) {

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

	    # Copy first input image to the output table file
	    # to preserve the header info
	    tcopy(fitsfile, outtmp)

	    # Get exposure time from file header
	    imgets (fitsfile//"[0]", "exptime")
	    ttot = real(imgets.value)

	    # Extract flux and error arrays into IRAF images
	    tabim(fitsfile, flxtmp, "FLUX", 1, naxis2)
	    tabim(fitsfile, errtmp, "ERROR", 1, naxis2)

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

	} else {	# For n > 0, go through the summation process

	# Get exposure time from file header
	    imgets (fitsfile//"[0]", "exptime")
	    texp = real(imgets.value)
	    ttot = ttot + texp

	# Create temporary files
	    flx = mktemp ("tmp$flx.")
	    err = mktemp ("tmp$err.")

	# Extract flux and error columns into IRAF image files and sum
	    tabim(fitsfile, flx, "FLUX", 1, naxis2)
	    imarith(flxtmp, "+", flx, flxtmp)

	    tabim(fitsfile, err, "ERROR", 1, naxis2)
	    imarith(err, "*", err, err)
	    imarith(errtmp, "+", err, errtmp)

	# Delete temporary files flx1 and err1 for use in next pass
	    imdelete(flx)
	    imdelete(err)

	}
	n = n + 1
    }

# Convert summed variance into output error array
    imfunc(errtmp, errtmp, "sqrt")

# Normalize by the number of combined spectra
    imarith(flxtmp, "/", n, flxtmp)
    imarith(errtmp, "/", n, errtmp)

# Create temporary file for summed tables
    flxtab = mktemp ("tmp$flxtab.")

# Put accumulated flux and error images into an STSDAS table file.
    imtab (flxtmp, flxtab, "FLUX", pname="", wcs="logical", formats="",\
	   tbltype="row")
    print("Added flux column to flxtab.")
    imtab (errtmp, flxtab, "ERROR", pname="", wcs="logical", formats="",\
	   tbltype="row")
    print("Added error column to flxtab.")

# Merge summed table with the original.  (Merging overwrites original
# columns with the summed fluxes and errors. Wavelengths stay the same.)
    tmerge (outtmp//","//flxtab, outfile, "merge", allcols=yes, tbltype="row",\
	    extracol=0)

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

# Update total exposure time in the output file
    hedit (outfile//"[0]", "EXPTIME", ttot,\
           add=no, delete=no, verify=no, show-, update=yes)

# Delete temporary images
    imdelete(flxtmp)
    imdelete(errtmp)
    delete(flxtab)
    list = ""

end