Repatch (from linux) of OSX IRAF

1 files changed, 2530 insertions, 0 deletions

diff --git a/noao/obsutil/src/sptime/t_sptime.x b/noao/obsutil/src/sptime/t_sptime.x
new file mode 100644
index 00000000..6d052acf
--- /dev/null
+++ b/noao/obsutil/src/sptime/t_sptime.x
@@ -0,0 +1,2530 @@
+include	<mach.h>
+include	<error.h>
+include <math.h>
+include <gset.h>
+include	<ctype.h>
+include	"sptime.h"
+
+
+# T_SPTIME -- Spectroscopic exposure time calculator.
+
+procedure t_sptime ()
+
+bool	interactive
+int	i, j, nexp, niter, npix, nw, fd, outlist
+real	nobj[4], nsky[4], time, minexp, maxexp, snr, sngoal
+real	wave, x, x1, dx, thruput, sat, dnmax
+pointer sp, str, err, st, tab, waves, counts, snrs, gp
+
+bool	streq(), tabexists(), fp_equalr()
+int	stgeti(), strdic()
+int	clpopnu(), fntopnb(), fntgfnb(), nowhite(), open(), tabgeti()
+real	stgetr(), tabgetr(), tabinterp1(), gr_mag(), gr_getr()
+real	st_x2w()
+pointer tabopen(), gopen(), un_open()
+errchk	tabopen, tabgeti, tabgetr, tabinterp1
+errchk	st_gtable, st_gtable1, stgeti, stgetr, st_snr, st_disperser
+errchk	open, gopen
+
+begin
+	call smark (sp)
+	call salloc (st, ST_LEN, TY_STRUCT)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (err, SZ_FNAME, TY_CHAR)
+
+	# Load tables.
+	ST_SEARCH(st) = clpopnu ("search")
+	ST_TAB(st) = tabopen ()
+	tab = ST_TAB(st)
+	call st_gtable (st, "spectrograph", "")
+	call st_gtable (st, "spectrum", "spectrograph")
+	call st_gtable (st, "sky", "spectrograph")
+	call st_gtable (st, "extinction", "spectrograph")
+	call st_gtable (st, "telescope", "spectrograph")
+	call st_gtable (st, "emissivity", "spectrograph")
+	call st_gtable (st, "adc", "spectrograph")
+	call st_gtable (st, "filter", "spectrograph")
+	call st_gtable (st, "filter2", "spectrograph")
+	call st_gtable (st, "aperture", "spectrograph")
+	call st_gtable (st, "fiber", "spectrograph")
+	call st_gtable (st, "aperture", "fiber")
+	call st_gtable (st, "collimator", "spectrograph")
+	call st_gtable (st, "disperser", "spectrograph")
+	call st_gtable (st, "xdisperser", "spectrograph")
+	call st_gtable (st, "corrector", "spectrograph")
+	call st_gtable (st, "camera", "spectrograph")
+	call st_gtable (st, "vignetting", "spectrograph")
+	call st_gtable (st, "vignetting", "camera")
+	call st_gtable (st, "detector", "spectrograph")
+	call st_gtable (st, "sensfunc", "spectrograph")
+	
+	call st_gtable1 (st, "abjohnson", "abjohnson")
+
+	# Set dispersion units.
+	call stgstr (st, "units", "spectrograph", "Angstroms",
+	    Memc[str], SZ_LINE)
+	call strcpy (Memc[str], ST_DUNITS(st), ST_SZSTRING) 
+	ST_DUNANG(st) = un_open ("angstroms")
+	ST_DUN(st) = un_open (ST_DUNITS(st))
+
+	# Set spectrum.
+	ST_REFW(st) = stgetr (st, "refwave", "spectrum", INDEFR)
+	if (!IS_INDEFR(ST_REFW(st)))
+	    call un_ctranr (ST_DUN(st), ST_DUNANG(st), ST_REFW(st),
+		ST_REFW(st), 1)
+	ST_REFF(st) = stgetr (st, "refflux", "spectrograph", 10.)
+	call stgstr (st, "funits", "spectrograph", "AB", Memc[str], SZ_LINE)
+	ST_FUNITS(st) = strdic (Memc[str], Memc[str], SZ_LINE, FUNITS)
+	switch (ST_SPEC(st)) {
+	case SPEC_BB:
+	    ST_PARAM(st) = stgetr (st, "temperature", "spectrograph", 6000.)
+	case SPEC_FL, SPEC_FN:
+	    ST_PARAM(st) = stgetr (st, "index", "spectrograph", 0.)
+	}
+	ST_RV(st) = stgetr (st, "R", "spectrum", 3.1)
+	ST_AV(st) = ST_RV(st) * stgetr (st, "E", "spectrum", 0.)
+
+	# Set observing conditions.
+	ST_SEEING(st) = stgetr (st, "seeing", "spectrograph", 1.)
+	ST_AIRMASS(st) = stgetr (st, "airmass", "spectrograph", 1.)
+	ST_PHASE(st) = stgetr (st, "phase", "spectrograph", 0.)
+
+	# Set thermal background.
+	ST_THERMAL(st) = stgetr (st, "thermal", "telescope", 0.)
+
+	# Set instrument.
+	ST_CW(st) = stgetr (st, "wave", "spectrograph", INDEFR)
+	if (!IS_INDEFR(ST_CW(st)))
+	    call un_ctranr (ST_DUN(st), ST_DUNANG(st), ST_CW(st), ST_CW(st), 1)
+	ST_ORDER(st,1) = stgeti (st, "order", "spectrograph", INDEFI)
+	ST_ORDER(st,2) = stgeti (st, "xorder", "spectrograph", INDEFI)
+
+	# Aperture.
+	if (!tabexists (tab, "aperture")) {
+	    if (tabexists (tab, "fiber"))
+		call st_gtable1 (st, "aperture", "circle")
+	    else if (!IS_INDEFR(stgetr(st, "diameter", "spectrograph", INDEFR)))
+		call st_gtable1 (st, "aperture", "circle")
+	    else
+		call st_gtable1 (st, "aperture", "slit")
+	}
+	call stgstr (st, "aptype", "aperture", "", Memc[str], SZ_LINE)
+	if (Memc[str] == EOS) {
+	    iferr (call tabgstr (tab, "aperture", "", "type",
+		Memc[str], SZ_LINE)) {
+		if (tabgeti (tab, "aperture", "", "table.ndim") == 2)
+		    call strcpy ("circular", Memc[str], SZ_LINE)
+		else
+		    call strcpy ("rectangular", Memc[str], SZ_LINE)
+	    }
+	}
+	ST_APTYPE(st) = strdic (Memc[str], Memc[str], SZ_LINE, APTYPES)
+	switch (ST_APTYPE(st)) {
+	case CIRCULAR:
+	    if (tabexists (tab, "fiber")) {
+		ST_APSIZE(st,1) = stgetr (st, "diameter", "fiber", INDEFR)
+		if (!IS_INDEFR(ST_APSIZE(st,1)) && ST_APSIZE(st,1) > 0.) {
+		    ST_TELSCALE(st) = stgetr (st, "scale", "telescope", 10.)
+		    ST_APSIZE(st,1) = ST_APSIZE(st,1) * ST_TELSCALE(st)
+		}
+	    }
+	    if (IS_INDEFR(ST_APSIZE(st,1)))
+		ST_APSIZE(st,1) = stgetr (st, "diameter", "aperture", -2.)
+	    ST_APSIZE(st,2) = ST_APSIZE(st,1)
+	case RECTANGULAR:
+	    ST_APSIZE(st,1) = stgetr (st, "width", "aperture", -2.)
+	    ST_APSIZE(st,2) = stgetr (st, "length", "aperture", -100.)
+	default:
+	    call sprintf (Memc[err], SZ_FNAME,
+		"Unknown aperture type (%s)")
+		call pargstr (Memc[str])
+	    call error (1, Memc[err])
+	}
+
+	ST_INOUTA(st,1) = stgetr (st, "inoutangle", "spectrograph", INDEFR)
+	ST_INOUTA(st,2) = stgetr (st, "xinoutangle", "spectrograph", INDEFR)
+	ST_BIN(st,1) = stgeti (st, "xbin", "detector", 1)
+	ST_BIN(st,2) = stgeti (st, "ybin", "detector", 1)
+	ST_GAIN(st) = stgetr (st, "gain", "detector", 1.)
+	ST_RDNOISE(st) = stgetr (st, "rdnoise", "detector", 0.)
+	ST_DARK(st) = stgetr (st, "dark", "detector", 0.)
+
+	# Set filter flag.
+	ST_FILTBLK(st) = NO
+	call stgstr (st, "block", "filter", "no", Memc[str], SZ_LINE) 
+	if (streq (Memc[str], "yes"))
+	    ST_FILTBLK (st) = YES
+
+	# Set sky subtraction parameters.
+	if (tabexists (tab, "sky") ||
+	    (tabexists (tab, "emissivity") && ST_THERMAL(st) > 0.)) {
+	    switch (ST_APTYPE(st)) {
+	    case CIRCULAR:
+		call stgstr (st, "skysub", "spectrograph", "multiap",
+		    Memc[str], SZ_LINE)
+		ST_NSKYAPS(st) = stgeti (st, "nskyaps", "spectrograph", 10)
+	    case RECTANGULAR:
+		call stgstr (st, "skysub", "spectrograph", "longslit",
+		    Memc[str], SZ_LINE)
+	    }
+	} else
+	    call stgstr (st, "skysub", "spectrograph", "none", Memc[str],
+		SZ_LINE)
+	    
+	ST_SKYSUB(st) = strdic (Memc[str], Memc[str], SZ_LINE, SKYSUB)
+
+	# Set calculation parameters.
+	ST_MINEXP(st) = stgetr (st, "minexp", "spectrograph", MINEXP)
+	ST_MAXEXP(st) = stgetr (st, "maxexp", "spectrograph", 3600.)
+	if (ST_MINEXP(st) <= 0.)
+	    ST_MINEXP(st) = MINEXP
+	if (ST_MAXEXP(st) <= ST_MINEXP(st))
+	    ST_MAXEXP(st) = ST_MINEXP(st)
+	time = stgetr (st, "time", "spectrograph", INDEFR)
+	sngoal = stgetr (st, "sn", "spectrograph", 25.)
+	if (IS_INDEF(time) && IS_INDEF(sngoal))
+	    call error (1,
+		"Either an exposure time or a desired S/N must be specified")
+	ST_SUBPIXELS(st) = stgeti (st, "subpixels", "spectrograph", 1)
+
+	# Set output parameters.
+	gp = NULL; fd = NULL
+	call stgstr (st, "output", "spectrograph", "", Memc[str], SZ_LINE)
+	outlist = fntopnb (Memc[str], NO)
+	if (fntgfnb (outlist, Memc[str], SZ_LINE) != EOF) {
+	    if (streq (Memc[str], "ALL")) {
+		call strcpy (OUTTYPES, Memc[str], SZ_LINE)
+		j = str+1
+		for (i=j; Memc[i] != EOS; i=i+1) { 
+		    if (IS_WHITE(Memc[i]) || Memc[i] == '\n')
+			next
+		    if (Memc[i] == Memc[str])
+			Memc[j] = ','
+		    else
+			Memc[j] = Memc[i]
+		    j = j + 1
+		}
+		Memc[j] = EOS
+		call fntclsb (outlist)
+		outlist = fntopnb (Memc[str+1], NO)
+	    } else
+		call fntrewb (outlist)
+	    nw = stgeti (st, "nw", "spectrograph", 101)
+	    call stgstr (st, "graphics", "spectrograph", "stdgraph",
+		Memc[str], SZ_LINE)
+	    if (nowhite (Memc[str], Memc[str], SZ_LINE) > 0) {
+		gp = gopen (Memc[str], NEW_FILE+AW_DEFER, STDGRAPH)
+		call stgstr (st, "interactive", "spectrograph", "yes",
+		    Memc[str], SZ_LINE)
+		interactive = streq (Memc[str], "yes")
+	    }
+	    call stgstr (st, "list", "spectrograph", "", Memc[str], SZ_LINE)
+	    if (nowhite (Memc[str], Memc[str], SZ_LINE) > 0)
+		fd = open (Memc[str], APPEND, TEXT_FILE)
+	}
+
+	# Focal lengths.
+	ST_COLFL(st) = stgetr (st, "colfl", "collimator", INDEFR)
+	if (IS_INDEFR(ST_COLFL(st))) {
+	    iferr (ST_COLFL(st) = tabgetr (tab, "collimator", "", "fl"))
+		ST_COLFL(st) = 1.
+	}
+	ST_CAMFL(st) = stgetr (st, "camfl", "camera", INDEFR)
+	if (IS_INDEFR(ST_CAMFL(st))) {
+	    iferr (ST_CAMFL(st) = tabgetr (tab, "camera", "", "fl"))
+		ST_CAMFL(st) = 1.
+	}
+
+	call st_disperser (st, "disperser", 1)
+	if (ST_DISPTYPE(st,1) == 0)
+	    call error (1, "No disperser specified")
+	call st_disperser (st, "xdisperser", 2)
+
+	ST_AREA(st) = 10000 * stgetr (st, "area", "telescope", 1.)
+
+	# Scales.
+	ST_TELSCALE(st) = stgetr (st, "scale", "telescope", 10.)
+	ST_SCALE(st,1) = ST_TELSCALE(st)
+	ST_SCALE(st,2) = ST_TELSCALE(st)
+	x = gr_mag (ST_GR(st,1), ST_CW(st), ST_ORDER(st,1))
+	if (!IS_INDEF(x))
+	    ST_SCALE(st,1) = ST_SCALE(st,1) / x
+	x = gr_mag (ST_GR(st,2), ST_CW(st), ST_ORDER(st,2))
+	if (!IS_INDEF(x))
+	    ST_SCALE(st,2) = ST_SCALE(st,2) / x
+	x = ST_COLFL(st) / ST_CAMFL(st)
+	ST_SCALE(st,1) = ST_SCALE(st,1) * x
+	ST_SCALE(st,2) = ST_SCALE(st,2) * x
+	ST_SCALE(st,1) = ST_SCALE(st,1) *
+	    stgetr (st, "apmagdisp", "spectrograph", 1.)
+	ST_SCALE(st,2) = ST_SCALE(st,2) *
+	    stgetr (st, "apmagxdisp", "spectrograph", 1.)
+	ST_PIXSIZE(st) = stgetr (st, "pixsize", "detector", 0.02)
+	ST_SCALE(st,1) = ST_SCALE(st,1) * ST_PIXSIZE(st) * ST_BIN(st,1)
+	ST_SCALE(st,2) = ST_SCALE(st,2) * ST_PIXSIZE(st) * ST_BIN(st,2)
+
+	# Convert aperture sizes to arcsec.
+	if (ST_APSIZE(st,1) < 0.)
+	    ST_APSIZE(st,1) = -ST_APSIZE(st,1) * ST_SCALE(st,1)
+	else
+	    ST_APSIZE(st,1) = ST_APSIZE(st,1)
+	if (ST_APSIZE(st,2) < 0.)
+	    ST_APSIZE(st,2) = -ST_APSIZE(st,2) * ST_SCALE(st,2)
+	else
+	    ST_APSIZE(st,2) = ST_APSIZE(st,2)
+
+	# Set dispersion per pixel and per resolution element.
+	ST_RES(st,1) = stgetr (st, "resolution", "camera", INDEFR)
+	if (IS_INDEFR(ST_RES(st,1)))
+	    ST_RES(st,1) = 2
+	else
+	    ST_RES(st,1) = ST_RES(st,1) / ST_PIXSIZE(st)
+	ST_RES(st,2) = ST_RES(st,1)
+	ST_DISP(st,1) = abs (gr_getr (ST_GR(st,1), "dispersion"))
+	ST_DISP(st,1) = ST_DISP(st,1) * ST_PIXSIZE(st) * ST_BIN(st,1)
+	x = 1 + min (ST_SEEING(st), ST_APSIZE(st,1)) / ST_SCALE(st,1)
+	ST_DISP(st,2) = ST_DISP(st,1) * max (2., ST_RES(st,1), x)
+
+	# Set number of pixels in object.
+	switch (ST_APTYPE(st)) {
+	case CIRCULAR:
+	    x = ST_APSIZE(st,2) / ST_SCALE(st,2) + ST_RES(st,2)
+	    npix = max (1, int (x + 0.999))
+	    ST_NOBJPIX(st) = npix
+	    ST_APLIM(st) = YES
+	case RECTANGULAR:
+	    x = ST_APSIZE(st,2) / ST_SCALE(st,2) + ST_RES(st,2)
+	    npix = max (1, int (x + 0.999))
+	    x = min (ST_APSIZE(st,2), 3*ST_SEEING(st)) / ST_SCALE(st,2) +
+		ST_RES(st,2)
+	    ST_NOBJPIX(st) = min (npix, int (x + 0.999))
+	    if (ST_NOBJPIX(st) > npix)
+		ST_APLIM(st) = NO
+	    else
+		ST_APLIM(st) = YES
+	}
+
+	# Set number of pixels in sky.
+	switch (ST_SKYSUB(st)) {
+	case SKY_NONE:
+	    ST_NSKYPIX(st) = 0
+	case SKY_LONGSLIT:
+	    ST_NSKYPIX(st) = max (0, npix - ST_NOBJPIX(st))
+	case SKY_MULTIAP:
+	    ST_NSKYPIX(st) = npix * ST_NSKYAPS(st)
+	case SKY_SHUFFLE:
+	    ST_NSKYPIX(st) = npix
+	}
+
+	# Compute exposure time and S/N.
+	if (!tabexists (tab, "spectrum")) {
+	    if (IS_INDEF(ST_REFW(st)))
+		ST_REFW(st) = ST_CW(st)
+	    switch (ST_FUNITS(st)) {
+	    case AB:
+		ST_REFFL(st) = 10. ** (-0.4 *
+		    (ST_REFF(st) + 5*log10(ST_REFW(st)) + 2.40))
+	    case FLAMBDA:
+		if (ST_REFF(st) < 0.)
+		    call error (1,
+			"Monochromatic flux (F-lambda) must be greater than 0")
+		ST_REFFL(st) = ST_REFF(st)
+	    case FNU:
+		if (ST_REFF(st) < 0.)
+		    call error (1,
+			"Monochromatic flux (F-nu) must be greater than 0")
+		ST_REFFL(st) = ST_REFF(st) * C / (ST_REFW(st) * ST_REFW(st))
+	    case UMAG,BMAG,VMAG,RMAG,IMAG,JMAG,HMAG,KSMAG,KMAG,LMAG,LPMAG,MMAG:
+		switch (ST_FUNITS(st)) {
+		case UMAG:
+		    ST_REFW(st) = 3650.
+		case BMAG:
+		    ST_REFW(st) = 4400.
+		case VMAG:
+		    ST_REFW(st) = 5500.
+		case RMAG:
+		    ST_REFW(st) = 7000.
+		case IMAG:
+		    ST_REFW(st) = 9000.
+		case JMAG:
+		    ST_REFW(st) = 12150.
+		case HMAG:
+		    ST_REFW(st) = 16540.
+		case KSMAG:
+		    ST_REFW(st) = 21570.
+		case KMAG:
+		    ST_REFW(st) = 21790.
+		case LMAG:
+		    ST_REFW(st) = 35470.
+		case LPMAG:
+		    ST_REFW(st) = 37610.
+		case MMAG:
+		    ST_REFW(st) = 47690.
+		}
+
+		ST_REFFL(st) = ST_REFF(st) +
+		    tabinterp1 (tab, "abjohnson", ST_REFW(st))
+		ST_REFFL(st) = 10. ** (-0.4 *
+		    (ST_REFFL(st) + 5*log10(ST_REFW(st)) + 2.40))
+	    }
+	}
+
+	# Check saturation.
+	sat = stgetr (st, "saturation", "detector", MAX_REAL)
+	dnmax = stgetr (st, "dnmax", "detector", MAX_REAL)
+
+	wave = ST_CW(st)
+	if (!IS_INDEF(time)) {
+	    if (time <= 0.) {
+		call eprintf ("Total Exposure Time must be greater than 0.\n")
+		return
+	    }
+
+	    if (ST_MAXEXP(st) > 0.) {
+		nexp = max (1, int (time / ST_MAXEXP(st) + 0.99))
+		time = time / nexp
+	    } else
+		nexp = 1
+
+	} else {
+	    if (sngoal <= 0.) {
+		call printf ("Desired S/N must be greater than 0.\n")
+		return
+	    }
+
+	    nexp = 1
+	    minexp = ST_MINEXP(st)
+	    maxexp = ST_MAXEXP(st)
+	    time = maxexp
+	    snr = 0.
+
+	    # Iterate to try and achieve the requested SNR.
+	    do niter = 1, MAXITER {
+
+		if (snr > 0.) {
+		    x = time
+		    i = nexp
+
+		    # After the first pass we use the last calculated SNR to
+		    # estimate a new time per exposure and number of exposures.
+		    time = time*sngoal*sngoal/(nexp*snr*snr)
+
+		    # Divide into multiple exposures if the time per exposure
+		    # exceeds a maximum.  Note the maximum may be reset by
+		    # saturation criteria.
+		    if (time > maxexp) {
+			nexp = nexp * max (1, int (time / maxexp + 0.99))
+			time = x*sngoal*sngoal/(nexp*snr*snr)
+		    } 
+
+		    # Apply a minimum time per exposure if possible.
+		    if (time < minexp && nexp > 1) {
+			nexp = max (1, nexp * time / minexp)
+			time = x*sngoal*sngoal/(nexp*snr*snr)
+		    }
+
+		    # New time per exposure to try.
+		    time = max (minexp, min (maxexp, time))
+		    if (fp_equalr (time, x) && nexp == i)
+			break
+		}
+
+		# Compute SNR.
+		call st_snr (st, NULL, wave, nexp, time, nobj, nsky,
+		    snr, thruput)
+
+		# Reset maximum time per exposure to avoid saturation.
+		if (nobj[1]+nsky[1] > sat && time > minexp && snr < sngoal) {
+		    time = time * nexp
+		    snr = snr * snr * nexp
+		    nexp = nexp * (1 + (nobj[1] + nsky[1]) / sat)
+		    time = time / nexp
+		    snr = sqrt (snr / nexp)
+		    maxexp = max (minexp, time)
+		    next
+		}
+		if ((nobj[1]+nsky[1])/ST_GAIN(st) > dnmax && time > minexp &&
+			snr < sngoal) {
+		    time = time * nexp
+		    snr = snr * snr * nexp
+		    nexp = nexp * (1 + (nobj[1] + nsky[1]) /
+			(ST_GAIN(st) * dnmax))
+		    time = time / nexp
+		    snr = sqrt (snr / nexp)
+		    maxexp = max (minexp, time)
+		    next
+		}
+
+		if (abs ((sngoal-sqrt(real(nexp))*snr)/sngoal) < 0.001)
+		    break
+	    }
+	}
+	ST_NEXP(st) = nexp
+	ST_TIME(st) = time
+
+	# Output.
+	npix = stgetr (st, "ndisp", "detector", 2048.)
+	nw = max (1, min (nw, npix))
+	x1 = npix * ST_PIXSIZE(st)
+
+	call salloc (waves, nw, TY_REAL)
+	call salloc (counts, nw, TY_REAL)
+	call salloc (snrs, nw, TY_REAL)
+
+	if (nw > 1) {
+	    dx = x1 / (nw - 1)
+	    x1 = -x1 / 2
+	    do i = 0, nw-1 {
+		x = x1 + dx * i
+		Memr[waves+i] = st_x2w (st, 1, x)
+	    }
+	} else
+	    Memr[waves] = wave
+
+	# Output result summary.
+	call st_results (st, STDOUT)
+	call st_check (st, STDOUT, Memr[waves], nw)
+	call st_snr (st, STDOUT, wave, nexp, time, nobj, nsky, snr, thruput)
+
+	while (fntgfnb (outlist, Memc[str], SZ_LINE) != EOF)
+	    call st_output (st, gp, fd, interactive, Memc[str], Memr[waves], nw)
+	    
+	# Finish up.
+	call un_close (ST_DUN(st))
+	call un_close (ST_DUNANG(st))
+	call clpcls (ST_SEARCH(st))
+	call tabclose (ST_TAB(st))
+	call gr_close (ST_GR(st,1))
+	call gr_close (ST_GR(st,2))
+	if (gp != NULL)
+	    call gclose (gp)
+	if (fd != NULL)
+	    call close (fd)
+	call fntclsb (outlist)
+	call sfree (sp)
+end
+
+
+# ST_RESULTS -- Print result summary.
+
+procedure st_results (st, fd)
+
+pointer st		#I SPECTIME structure
+int	fd		#I Output file descriptor
+
+char	eff[SZ_FNAME]
+real	x, y
+int	npix, order, tabgeti(), stgeti()
+pointer	tab
+real	gr_getr()
+bool	tabexists()
+
+begin
+	tab = ST_TAB(st)
+
+	call fprintf (fd, "\n")
+	if (tabexists (tab, "spectrum"))
+	    call st_description (st, fd, "Object spectrum: ", "spectitle",
+		"spectrum")
+	else {
+	    call fprintf (fd, "Object spectrum: ")
+	    switch (ST_SPEC(st)) {
+	    case SPEC_BB:
+		call fprintf (fd, "Blackbody spectrum of temperature %g K\n")
+		    call pargr (ST_PARAM(st))
+	    case SPEC_FL:
+		call fprintf (fd, "F(lambda) power law of index %g\n")
+		    call pargr (ST_PARAM(st))
+	    case SPEC_FN:
+		call fprintf (fd, "F(nu) power law of index %g\n")
+		    call pargr (ST_PARAM(st))
+	    }
+	    if (ST_AV(st) > 0.) {
+		call fprintf (fd,
+		    "Reddening: E(B-V) of %g with A(V)/E(B-V) of %g\n")
+		    call pargr (ST_AV(st) / ST_RV(st))
+		    call pargr (ST_RV(st))
+	    }
+	    call un_ctranr (ST_DUNANG(st), ST_DUN(st), ST_REFW(st), x, 1)
+	    call fprintf (fd, "Reference wavelength: %.4g %s\n")
+		call pargr (x)
+		call pargstr (ST_DUNITS(st))
+	    call fprintf (fd, "Reference flux: ")
+	    switch (ST_FUNITS(st)) {
+	    case AB:
+		call fprintf (fd, "AB = %.3g (%.3g ergs/s/cm^2/A)\n")
+		    call pargr (ST_REFF(st))
+		    call pargr (ST_REFFL(st))
+	    case FLAMBDA:
+		call fprintf (fd, "%.3g ergs/s/cm^2/A\n")
+		    call pargr (ST_REFF(st))
+	    case FNU:
+		call fprintf (fd, "%.3g ergs/s/cm^2/Hz\n")
+		    call pargr (ST_REFF(st))
+	    case UMAG, BMAG, VMAG, RMAG, IMAG, JMAG:
+		switch (ST_FUNITS(st)) {
+		case UMAG:
+		    call fprintf (fd, "U = %.3g ")
+			call pargr (ST_REFF(st))
+		case BMAG:
+		    call fprintf (fd, "B = %.3g ")
+			call pargr (ST_REFF(st))
+		case VMAG:
+		    call fprintf (fd, "V = %.3g ")
+			call pargr (ST_REFF(st))
+		case RMAG:
+		    call fprintf (fd, "R = %.3g ")
+			call pargr (ST_REFF(st))
+		case IMAG:
+		    call fprintf (fd, "I = %.3g ")
+			call pargr (ST_REFF(st))
+		case JMAG:
+		    call fprintf (fd, "J = %.3g ")
+			call pargr (ST_REFF(st))
+		}
+		call fprintf (fd, "(%.3g ergs/s/cm^2/A)\n")
+		    call pargr (ST_REFFL(st))
+	    }
+	}
+	if (tabexists (tab, "sky")) {
+	    call st_description (st, fd, "Sky spectrum: ", "skytitle", "sky")
+	    if (tabgeti (tab, "sky", "", "table.ndim") == 2) {
+		call fprintf (fd, "\tMoon phase: %d\n")
+		    call pargr (ST_PHASE(st))
+	    }
+	}
+	if (ST_AIRMASS(st) > 0) {
+	    call st_description (st, fd, "Extinction: ", "exttitle",
+		"extinction")
+	    call fprintf (fd, "\tAirmass: %.3g\n")
+		call pargr (ST_AIRMASS(st))
+	}
+	call fprintf (fd, "Seeing: %.2g\" (FWHM)\n")
+	    call pargr (ST_SEEING(st))
+	if (tabexists (tab, "emissivity") && ST_THERMAL(st) > 0.) {
+	    call fprintf (fd, "Thermal Background:\n")
+	    call st_description (st, fd, "\tEmissivity: ",
+		"emistitle", "emissivity")
+	    call fprintf (fd, "\tTemperature: %.1g K\n")
+		call pargr (ST_THERMAL(st))
+	}
+
+	call fprintf (fd, "\n")
+	call st_description (st, fd, "Telescope: ", "teltitle", "telescope")
+	call fprintf (fd, "\tArea: %.1f m^2, Scale: %.4g arcsec/mm\n")
+	    call pargr (ST_AREA(st) / 10000.)
+	    call pargr (ST_TELSCALE(st))
+	if (tabexists (tab, "adc"))
+	    call st_description (st, fd, "ADC: ", "adctitle", "adc")
+	if (tabexists (tab, "spectrograph"))
+	    call st_description (st, fd, "Spectrograph: ", "title",
+		"spectrograph")
+	call st_description (st, fd, "Collimator: ", "coltitle", "collimator")
+	call fprintf (fd, "\tFocal length = %.4g m\n")
+	    call pargr (ST_COLFL(st))
+	if (tabexists (tab, "aperture")) {
+	    call st_description (st, fd, "Apertures: ", "aptitle", "aperture")
+	    switch (ST_APTYPE(st)) {
+	    case CIRCULAR:
+		call fprintf (fd, "\tSize: %.2f\", %.3g mm, %.1f pixels\n")
+		    call pargr (ST_APSIZE(st,1))
+		    call pargr (ST_APSIZE(st,1) / ST_TELSCALE(st))
+		    call pargr (ST_APSIZE(st,1) / ST_SCALE(st,1))
+	    case RECTANGULAR:
+		call fprintf (fd,
+		"\tSize: %.2f\" x %.2f\", %.3g x %.3g mm, %.1f x %.1f pixels\n")
+		    call pargr (ST_APSIZE(st,1))
+		    call pargr (ST_APSIZE(st,2))
+		    call pargr (ST_APSIZE(st,1) / ST_TELSCALE(st))
+		    call pargr (ST_APSIZE(st,2) / ST_TELSCALE(st))
+		    call pargr (ST_APSIZE(st,1) / ST_SCALE(st,1))
+		    call pargr (ST_APSIZE(st,2) / ST_SCALE(st,2))
+	    }
+	}
+	if (tabexists (tab, "fiber"))
+	    call st_description (st, fd, "Fibers: ", "fibtitle", "fiber")
+	if (tabexists (tab, "filter"))
+	    call st_description (st, fd, "Filter: ", "ftitle", "filter")
+	if (tabexists (tab, "filter2"))
+	    call st_description (st, fd, "Filter: ", "f2title", "filter2")
+	if (ST_DISPTYPE(st,1) != 0) {
+	    call st_description (st, fd, "Disperser: ", "disptitle",
+		"disperser")
+	    order = nint (gr_getr (ST_GR(st,1), "order"))
+	    call fprintf (fd, "\tCentral order = %d\n")
+		call pargi (order)
+	    x = gr_getr (ST_GR(st,1), "wavelength")
+	    call un_ctranr (ST_DUNANG(st), ST_DUN(st), x, x, 1)
+	    call fprintf (fd, "\tCentral wavelength = %.6g %s\n")
+		call pargr (x)
+		call pargstr (ST_DUNITS(st))
+	    x = abs(gr_getr(ST_GR(st,1), "dispersion"))
+	    call un_ctranr (ST_DUNANG(st), ST_DUN(st), x, x, 1)
+	    call un_ctranr (ST_DUNANG(st), ST_DUN(st), ST_DISP(st,1), y, 1)
+	    call fprintf (fd,
+	    "\tCentral dispersion = %.3g %s/mm, %.3g %s/pixel\n")
+		call pargr (x)
+		call pargstr (ST_DUNITS(st))
+		call pargr (y)
+		call pargstr (ST_DUNITS(st))
+	    if (ST_DISPTYPE(st,1) == GRATING &&
+		int(gr_getr(ST_GR(st,1),"full"))==YES) {
+		call fprintf (fd, "\tRuling = %d lines/mm\n")
+		    call pargr (gr_getr (ST_GR(st,1), "g"))
+		call fprintf (fd, "\tBlaze = %.1f deg\n")
+		    call pargr (gr_getr (ST_GR(st,1), "blaze"))
+		x = gr_getr (ST_GR(st,1), "tilt")
+		if (abs(x) > 0.1) {
+		    call fprintf (fd, "\tGrating tilt = %.1f degrees\n")
+			call pargr (gr_getr (ST_GR(st,1), "tilt"))
+		    x = gr_getr (ST_GR(st,1), "mag")
+		    if (abs(x) < 0.99) {
+			call fprintf (fd, "\tGrating magnification = %.2f\n")
+			    call pargr (x)
+		    }
+		}
+		call sprintf (eff, SZ_FNAME, "eff%d")
+		    call pargi (order)
+		if (!tabexists (tab, eff)) {
+		    if (order == 1 && tabexists (tab, "disperser")) {
+			if (tabgeti (tab, "disperser", "", "table.ndim") != 0)
+			    call strcpy ("disperser", eff, SZ_FNAME)
+		    }
+		}
+		if (tabexists (tab, eff))
+		    call fprintf (fd, "\tUsing tabulated efficiencies\n")
+		else
+		    call fprintf (fd, "\tUsing predicted efficiencies\n")
+	    }
+	}
+	if (ST_DISPTYPE(st,2) != 0) {
+	    call st_description (st, fd, "Crossdisperser: ", "xdisptitle",
+		"xdisperser")
+	    order = nint (gr_getr (ST_GR(st,2), "order"))
+	    call fprintf (fd, "\tCentral order = %d\n")
+		call pargi (order)
+	    x = gr_getr (ST_GR(st,2), "wavelength")
+	    call un_ctranr (ST_DUNANG(st), ST_DUN(st), x, x, 1)
+	    call fprintf (fd, "\tCentral wavelength = %.6g %s\n")
+		call pargr (x)
+		call pargstr (ST_DUNITS(st))
+	    x = abs(gr_getr(ST_GR(st,1), "dispersion"))
+	    call un_ctranr (ST_DUNANG(st), ST_DUN(st), x, x, 1)
+	    call fprintf (fd, "\tCentral dispersion = %.3g %s/mm\n")
+		call pargr (x)
+		call pargstr (ST_DUNITS(st))
+	    if (ST_DISPTYPE(st,2) == GRATING &&
+		int(gr_getr(ST_GR(st,2),"full"))==YES) {
+		call fprintf (fd, "\tRuling = %d lines/mm\n")
+		    call pargr (gr_getr (ST_GR(st,2), "g"))
+		call fprintf (fd, "\tBlaze = %d deg\n")
+		    call pargr (gr_getr (ST_GR(st,2), "blaze"))
+		x = gr_getr (ST_GR(st,2), "tilt")
+		if (abs(x) > 0.1) {
+		    call fprintf (fd, "\tGrating tilt = %.1f degrees\n")
+			call pargr (gr_getr (ST_GR(st,2), "tilt"))
+		    x = gr_getr (ST_GR(st,2), "mag")
+		    if (abs(x) < 0.99) {
+			call fprintf (fd, "\tGrating magnification = %.2f\n")
+			    call pargr (x)
+		    }
+		}
+		call sprintf (eff, 10, "xeff%d")
+		    call pargi (order)
+		if (!tabexists (tab, eff)) {
+		    if (order == 1 && tabexists (tab, "xdisperser")) {
+			if (tabgeti (tab, "xdisperser", "", "table.ndim") != 0)
+			    call strcpy ("xdisperser", eff, SZ_FNAME)
+		    }
+		}
+		if (tabexists (tab, eff))
+		    call fprintf (fd, "\tUsing tabulated efficiencies\n")
+		else
+		    call fprintf (fd, "\tUsing predicted efficiencies\n")
+	    }
+	}
+	if (tabexists (tab, "corrector"))
+	    call st_description (st, fd, "Corrector: ", "cortitle", "corrector")
+	call st_description (st, fd, "Camera: ", "camtitle", "camera")
+	call fprintf (fd, "\tFocal length = %.4g m, Resolution = %.1g pixels\n")
+	    call pargr (ST_CAMFL(st))
+	    call pargr (ST_RES(st,1))
+	call st_description (st, fd, "Detector: ", "dettitle", "detector")
+	call fprintf (fd, "\tPixel size: %d microns, %.2f\"\n")
+	    call pargr (1000 * ST_PIXSIZE(st))
+	    call pargr (ST_SCALE(st,1))
+	if (ST_BIN(st,1) != 1 || ST_BIN(st,2) != 1) {
+	    call fprintf (fd, "\tBinning: %dx%d (XxY)\n")
+		call pargi (ST_BIN(st,1))
+		call pargi (ST_BIN(st,2))
+	}
+	npix = stgeti (st, "ndisp", "detector", 2048) / ST_BIN(st,1)
+	call fprintf (fd, "\tNumber of pixels = %d\n")
+	    call pargi (npix)
+	call fprintf (fd,
+	    "\tRead noise = %.1f e-, Gain = %.1f e-/DN, Dark = %.3g e-/s\n") 
+	    call pargr (ST_RDNOISE(st))
+	    call pargr (ST_GAIN(st))
+	    call pargr (ST_DARK(st))
+
+	call fprintf (fd, "\n")
+	call fprintf (fd,
+	    "Source pixels: %d pixels")
+	    call pargi (ST_NOBJPIX(st))
+	if (ST_APLIM(st) == YES)
+	    call fprintf (fd, " (aperture & camera resolution limited)\n")
+	else
+	    call fprintf (fd, " (3xFWHM seeing disk & camera resolution)\n")
+	switch (ST_SKYSUB(st)) {
+	case SKY_NONE:
+	    call fprintf (fd,
+		"Background pixels: no background subtraction done\n")
+	case SKY_LONGSLIT:
+	    call fprintf (fd, "Background pixels: %3d pixels\n")
+		call pargi (ST_NSKYPIX(st))
+	case SKY_MULTIAP:
+	    call fprintf (fd,
+		"Background pixels: %d apertures each with %d pixels\n")
+		call pargi (ST_NSKYAPS(st))
+		call pargi (ST_NSKYPIX(st) / ST_NSKYAPS(st))
+	case SKY_SHUFFLE:
+	    call fprintf (fd,
+		"Background pixels: shuffle with %d pixels (same as source)\n")
+		call pargi (ST_NSKYPIX(st))
+	}
+	call fprintf (fd, "\n")
+end
+
+
+# ST_DESCRIPTION -- Print description of table.
+
+procedure st_description (st, fd, label, title, table)
+
+pointer st		#I SPECTIME structure
+int	fd		#I Output file descriptor
+char	label[ARB]	#I Label
+char	title[ARB]	#I Title parameter name
+char	table[ARB]	#I Table name
+
+pointer sp, str
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Print label and title.
+	call st_gtitle (st, title, table, Memc[str], SZ_LINE)
+	call fprintf (fd, "%s%s\n")
+	    call pargstr (label)
+	    call pargstr (Memc[str])
+
+	# Print description if available.
+	ifnoerr (call tabgstr (ST_TAB(st), table, "", "description", Memc[str],
+	    SZ_LINE)) {
+	    call fprintf (fd, "%s\n")
+		call pargstr (Memc[str])
+	}
+
+	call sfree (sp)
+end
+
+
+# ST_GTITLE -- Get title.
+
+procedure st_gtitle (st, param, table, title, maxchar)
+
+pointer st		#I SPECTIME structure
+char	param[ARB]	#I Title parameter name
+char	table[ARB]	#I Table name
+char	title[ARB]	#O Title
+int	maxchar		#I Maximum string length
+
+char	dummy[1]
+int	nowhite()
+
+begin
+	call clgstr (param, title, maxchar)
+	if (nowhite (title, dummy, 1) > 0)
+	    return
+
+	iferr (call tabgstr (ST_TAB(st), table, "spectrograph", param, title,
+	    maxchar)) {
+	    iferr (call tabgstr (ST_TAB(st), table, "", "title",
+		title, maxchar)) {
+		iferr (call tabgstr (ST_TAB(st), table, "",
+		    "table.filename", title, maxchar))
+		    title[1] = EOS
+	    }
+	}
+end
+
+
+# ST_SNR -- Source, sky, SNR, and thruput including all orders.
+
+procedure st_snr (st, fd, wave, nexp, time, nobj, nsky, snr, thruput)
+
+pointer st		#I SPECTIME structure
+int	fd		#I Output descriptor
+real	wave		#I Wavelengths
+int	nexp		#I Number of exposures
+real	time		#I Exposure time
+real	nobj[4]		#O Number of object photons (1=total)
+real	nsky[4]		#O Number of sky photons (1=total)
+real	snr		#O S/N per pixel for total
+real	thruput		#O System thruput at primary wavelength
+
+int	i, j, order
+real	w, f, st_spectrum()
+errchk	st_spectrum, st_snr1
+
+begin
+	# Initialize.
+	do i = 1, 4 {
+	    nobj[i] = 0
+	    nsky[i] = 0
+	}
+
+	# If not x-dispersed and disperser is a grating do overlapping orders.
+	if (ST_DISPTYPE(st,2) == 0 && ST_FILTBLK(st) == NO &&
+	    (ST_DISPTYPE(st,1) == GRATING || ST_DISPTYPE(st,1) == GRISM)) {
+	    order = ST_ORDER(st,1)
+	    do j = 2, 4 {
+		i = order + j - 3
+		if (i < 1)
+		    next
+		w = wave * order / i
+		f = st_spectrum (st, w)
+		ST_ORDER(st,1) = i
+		call st_snr1 (st, NULL, w, f, nexp, time, nobj[j], nsky[j],
+		    snr, thruput)
+		if (i != order) {
+		    nobj[1] = nobj[1] + nobj[j]
+		    nsky[1] = nsky[1] + nsky[j]
+		}
+	    }
+	    ST_ORDER(st,1) = order
+
+	    w = wave
+	    f = st_spectrum (st, w)
+	    call st_snr1 (st, fd, w, f, nexp, time, nobj, nsky, snr, thruput)
+	} else {
+	    w = wave
+	    f = st_spectrum (st, w)
+	    call st_snr1 (st, fd, w, f, nexp, time, nobj, nsky, snr, thruput)
+	    nobj[3] = nobj[1]
+	    nsky[3] = nsky[1]
+	}
+end
+
+
+# ST_SNR1 -- Compute and print photons and S/N for given exposure time.
+
+procedure st_snr1 (st, fd, wave, flux, nexp, time, nobj, nbkg, snr, thruput)
+
+pointer st		#I SPECTIME structure
+int	fd		#I Output descriptor
+real	wave		#I Wavelength
+real	flux		#I Flux
+int	nexp		#I Number of exposures
+real	time		#I Exposure time
+real	nobj		#U Number of object photons
+real	nbkg		#U Number of bkg photons
+real	snr		#O S/N per pixel
+real	thruput		#O System thruput
+
+int	i, ncols
+real	tobs, n, ndark
+real	w, w1, nmag, sky, thermal, bkg, dobj, dbkg, ddark, dreadout, tnoise
+real	ext, tel, adc, ap, fib, inst, fltr1, fltr2, col, eff1, eff2, disp
+real	cor, cam, vig, dqe, cum, sqnexp
+
+real	tabinterp1(), tabinterp2(), st_dispeff(), st_w2dw(), st_w2x()
+errchk	tabinterp1, tabinterp2, st_dispeff
+
+begin
+	# Check for reasonable wavlength and source flux.
+	if (wave < 0. || flux < 0.) {
+	    nobj = 0.
+	    nbkg = 0.
+	    snr = 0.
+	    thruput = 0.
+	    return
+	}
+
+	# Set observation time.
+	switch (ST_SKYSUB(st)) {
+	case SKY_SHUFFLE:
+	    tobs = time / 2
+	default:
+	    tobs = time
+	}
+
+	# Compute pixel counts over subsampled pixels.
+	disp = st_w2dw (st, 1, wave) * ST_PIXSIZE(st) * ST_BIN(st,1) /
+	    ST_SUBPIXELS(st)
+	w1 = wave - disp * (ST_SUBPIXELS(st) + 1) / 2.
+	do i = 1, ST_SUBPIXELS(st) {
+	    w = w1 + disp * i
+
+	    # Atmospheric transmission.
+	    iferr {
+		ext = tabinterp1 (ST_TAB(st), "extinction", w)
+		ext = 10 ** (-0.4 * ext * ST_AIRMASS(st))
+	    } then
+		ext = INDEF
+
+	    # Telescope transmission.
+	    iferr (tel = tabinterp1 (ST_TAB(st), "telescope", w))
+		tel = 1
+
+	    # ADC transmission.
+	    iferr (adc = tabinterp1 (ST_TAB(st), "adc", w))
+		adc = INDEF
+
+	    # Aperture thruput.
+	    iferr {
+		switch (ST_APTYPE(st)) {
+		case CIRCULAR:
+		    ap = tabinterp1 (ST_TAB(st), "aperture",
+			ST_APSIZE(st,1) / ST_SEEING(st))
+		case RECTANGULAR:
+		    ap = tabinterp2 (ST_TAB(st), "aperture",
+			ST_APSIZE(st,1) / ST_SEEING(st),
+			ST_APSIZE(st,2) / ST_SEEING(st))
+		}
+	    } then
+		ap = 1
+
+	    # Fiber transmission.
+	    iferr (fib = tabinterp1 (ST_TAB(st), "fiber", w))
+		fib = INDEF
+
+	    # Spectrograph transmission.
+	    iferr (inst = tabinterp1 (ST_TAB(st), "spectrograph", w))
+		inst = INDEF
+
+	    # Filter transmission.
+	    iferr (fltr1 = tabinterp1 (ST_TAB(st), "filter", w))
+		fltr1 = INDEF
+	    iferr (fltr2 = tabinterp1 (ST_TAB(st), "filter2", w))
+		fltr2 = INDEF
+
+	    # Collimator transmission.
+	    iferr (col = tabinterp1 (ST_TAB(st), "collimator", w))
+		col = 1
+
+	    # Disperser efficiency.
+	    eff1 = st_dispeff (st, "disperser", w, ST_ORDER(st,1))
+	    eff2 = st_dispeff (st, "xdisperser", w, ST_ORDER(st,2))
+
+	    # Corrector transmission.
+	    iferr (cor = tabinterp1 (ST_TAB(st), "corrector", w))
+		cor = INDEF
+
+	    # Camera transmission.
+	    iferr (cam = tabinterp1 (ST_TAB(st), "camera", w))
+		cam = 1
+
+	    # Vignetting.
+	    iferr (vig = tabinterp1 (ST_TAB(st), "vignetting",
+		st_w2x (st, 1, w)))
+		vig = INDEF
+
+	    # Detector DQE.
+	    iferr (dqe = tabinterp1 (ST_TAB(st), "detector", w))
+		dqe = 1
+
+	    # Cumulative transmission.
+	    thruput = 1
+	    if (!IS_INDEF(tel)) {
+		tel = max (0., tel)
+		thruput = thruput * tel
+	    }
+	    if (!IS_INDEF(adc)) {
+		adc = max (0., adc)
+		thruput = thruput * adc
+	    }
+	    if (!IS_INDEF(ap)) {
+		ap = max (0., ap)
+		thruput = thruput * ap
+	    }
+	    if (!IS_INDEF(fib)) {
+		fib = max (0., fib)
+		thruput = thruput * fib
+	    }
+	    if (!IS_INDEF(inst)) {
+		inst = max (0., inst)
+		thruput = thruput * inst
+	    }
+	    if (!IS_INDEF(fltr1)) {
+		fltr1 = max (0., fltr1)
+		thruput = thruput * fltr1
+	    }
+	    if (!IS_INDEF(fltr2)) {
+		fltr2 = max (0., fltr2)
+		thruput = thruput * fltr2
+	    }
+	    if (!IS_INDEF(eff1)) {
+		eff1 = max (0., eff1)
+		thruput = thruput * eff1
+	    }
+	    if (!IS_INDEF(eff2)) {
+		eff2 = max (0., eff2)
+		thruput = thruput * eff2
+	    }
+	    if (!IS_INDEF(cor)) {
+		cor = max (0., cor)
+		thruput = thruput * cor
+	    }
+	    if (!IS_INDEF(col)) {
+		col = max (0., col)
+		thruput = thruput * col
+	    }
+	    if (!IS_INDEF(cam)) {
+		cam = max (0., cam)
+		thruput = thruput * cam
+	    }
+	    if (!IS_INDEF(vig)) {
+		vig = max (0., vig)
+		thruput = thruput * vig
+	    }
+	    if (!IS_INDEF(dqe)) {
+		dqe = max (0., dqe)
+		thruput = thruput * dqe
+	    }
+
+	    # Source photons detected.
+	    nmag = flux / (C * H / w)
+	    n = nmag * ST_AREA(st) * thruput * tobs * disp
+	    if (!IS_INDEF(ext))
+		n = n * ext
+	    n = max (0., n)
+	    if (n < 100000.)
+		n = int (n)
+	    nobj = nobj + n
+
+	    # Sky photon flux.
+	    iferr (sky = tabinterp2 (ST_TAB(st), "sky", w, ST_PHASE(st)))
+		iferr (sky = tabinterp1 (ST_TAB(st), "sky", w))
+		    sky = 0
+	    sky = sky / (C * H / w)
+
+	    # Thermal photon flux.
+	    thermal = 0.
+	    if (!IS_INDEF(ST_THERMAL(st))) {
+		iferr {
+		    thermal = tabinterp1 (ST_TAB(st), "emissivity", w)
+		    # 1.41e24 = 2 * c * (arcsec/rad)^2 * (A/cm) * (A/cm)^-4
+		    thermal = thermal * 1.41e24 / (w ** 4) /
+			(exp (min (30., 1.43877e8 / (w * ST_THERMAL(st)))) - 1)
+		} then
+		    thermal = 0.
+	    }
+
+	    # Total background.
+	    bkg = sky + thermal
+
+	    switch (ST_APTYPE(st)) {
+	    case CIRCULAR:
+		bkg = bkg * PI * (ST_APSIZE(st,1) / 2) ** 2
+	    case RECTANGULAR:
+		bkg = bkg * ST_APSIZE(st,1) * ST_SCALE(st,1) * ST_NOBJPIX(st)
+	    }
+	    n = bkg * ST_AREA(st) * thruput / ap * tobs * disp
+	    n = max (0., n)
+	    if (n < 100000.)
+		n = int (n)
+	    nbkg = nbkg + n
+	}
+
+	# Dark counts.
+	ndark = ST_NOBJPIX(st) * ST_DARK(st) * time
+	ndark = max (0., ndark)
+	if (ndark <100000.)
+	    ndark = int (ndark)
+
+	# Noise.
+	dobj = sqrt (nobj)
+	dbkg = sqrt (nbkg)
+	ddark = sqrt (ndark)
+	dreadout = sqrt (ST_NOBJPIX(st) * ST_RDNOISE(st)**2)
+	tnoise = nobj + nbkg + ndark + dreadout**2
+
+	# Background subtraction statistics.
+	switch (ST_SKYSUB(st)) {
+	case SKY_NONE:
+	    nobj = nobj + nbkg
+	    nbkg = 0.
+	default:
+	    if (ST_NSKYPIX(st) > 0)
+		tnoise = tnoise + ((nbkg + ndark) / ST_NOBJPIX(st) +
+		    ST_RDNOISE(st)**2) / ST_NSKYPIX(st)
+	}
+
+	# Final S/N.
+	snr = nobj
+	tnoise = sqrt (tnoise)
+	if (tnoise > 0.)
+	    snr = snr / tnoise
+
+	if (fd == NULL)
+	    return
+
+	# Print results.
+	sqnexp = sqrt (real(nexp))
+	ncols = nint (ST_DISP(st,2) / ST_DISP(st,1))
+
+	call un_ctranr (ST_DUNANG(st), ST_DUN(st), wave, w, 1)
+	if (nexp > 1) {
+	    call fprintf (fd,
+	    "---- Results for %d exposures of %.2fs at %.4g %s ----\n")
+		call pargi (nexp)
+		call pargr (time)
+		call pargr (w)
+		call pargstr (ST_DUNITS(st))
+	} else {
+	    call fprintf (fd,
+	    "---- Results for %.2fs exposure at %.4g %s ----\n")
+		call pargr (time)
+		call pargr (w)
+		call pargstr (ST_DUNITS(st))
+	}
+
+	call fprintf (fd, "\nSource flux: %.3g photons/s/cm^2/A (AB=%.3g)\n")
+	    call pargr (nmag)
+	    call pargr (-2.5*log10(flux) - 5*log10(wave) - 2.40)
+	call fprintf (fd, "Background flux: %.3g photons/s/cm^2/A (over source pixels)\n")
+	    call pargr (bkg)
+
+	call fprintf (fd, "\nTransmision/Efficiencies:%40t%10s %10s\n")
+	    call pargstr ("individual")
+	    call pargstr ("cumulative")
+	cum = 1
+	if (!IS_INDEF(ext) && ext < 1) {
+	    cum = cum * ext
+	    call fprintf (fd,
+		"    Atmosphere (%.2g mag/airmass)%40t%9.1f%% %9.1f%%\n")
+		call pargr (-2.5 * log10 (ext) / ST_AIRMASS(st))
+		call pargr (100 * ext)
+		call pargr (100 * cum)
+	}
+	if (!IS_INDEF(tel) && tel < 1) {
+	    cum = cum * tel
+	    call fprintf (fd, "    Telescope%40t%9.1f%% %9.1f%%\n")
+		call pargr (100 * tel)
+		call pargr (100 * cum)
+	}
+	if (!IS_INDEF(adc) && adc < 1) {
+	    cum = cum * adc
+	    call fprintf (fd, "    ADC%40t%9.1f%% %9.1f%%\n")
+		call pargr (100 * adc)
+		call pargr (100 * cum)
+	}
+	if (!IS_INDEF(ap) && ap < 1) {
+	    cum = cum * ap
+	    call fprintf (fd,
+		"    Aperture (seeing=%.2g\")%40t%9.1f%% %9.1f%%\n")
+		call pargr (ST_SEEING(st))
+		call pargr (100 * ap)
+		call pargr (100 * cum)
+	}
+	if (!IS_INDEF(fib) && fib < 1) {
+	    cum = cum * fib
+	    call fprintf (fd, "    Fiber%40t%9.1f%% %9.1f%%\n")
+		call pargr (100 * fib)
+		call pargr (100 * cum)
+	}
+	if (!IS_INDEF(fltr1) && fltr1 < 1) {
+	    cum = cum * fltr1
+	    call fprintf (fd, "    Filter%40t%9.1f%% %9.1f%%\n")
+		call pargr (100 * fltr1)
+		call pargr (100 * cum)
+	}
+	if (!IS_INDEF(fltr2) && fltr2 < 1) {
+	    cum = cum * fltr2
+	    call fprintf (fd, "    Filter%40t%9.1f%% %9.1f%%\n")
+		call pargr (100 * fltr2)
+		call pargr (100 * cum)
+	}
+	if (!IS_INDEF(inst) && inst < 1) {
+	    cum = cum * inst
+	    call fprintf (fd, "    Spectrograph%40t%9.1f%% %9.1f%%\n")
+		call pargr (100 * inst)
+		call pargr (100 * cum)
+	}
+	if (!IS_INDEF(col) && col < 1) {
+	    cum = cum * col
+	    call fprintf (fd, "    Collimator%40t%9.1f%% %9.1f%%\n")
+		call pargr (100 * col)
+		call pargr (100 * cum)
+	}
+	if (!IS_INDEF(eff1) && eff1 < 1) {
+	    cum = cum * eff1
+	    call fprintf (fd, "    Disperser%40t%9.1f%% %9.1f%%\n")
+		call pargr (100 * eff1)
+		call pargr (100 * cum)
+	}
+	if (!IS_INDEF(eff2) && eff2 < 1) {
+	    cum = cum * eff2
+	    call fprintf (fd, "    Cross disperser%40t%9.1f%% %9.1f%%\n")
+		call pargr (100 * eff2)
+		call pargr (100 * cum)
+	}
+	if (!IS_INDEF(cor) && cor < 1) {
+	    cum = cum * cor
+	    call fprintf (fd, "    Corrector%40t%9.1f%% %9.1f%%\n")
+		call pargr (100 * cor)
+		call pargr (100 * cum)
+	}
+	if (!IS_INDEF(cam) && cam < 1) {
+	    cum = cum * cam
+	    call fprintf (fd, "    Camera%40t%9.1f%% %9.1f%%\n")
+		call pargr (100 * cam)
+		call pargr (100 * cum)
+	}
+	if (!IS_INDEF(vig) && vig < 0.999) {
+	    cum = cum * vig
+	    call fprintf (fd, "    Vignetting%40t%9.1f%% %9.1f%%\n")
+		call pargr (100 * vig)
+		call pargr (100 * vig)
+	}
+	if (!IS_INDEF(dqe) && dqe < 1) {
+	    cum = cum * dqe
+	    call fprintf (fd, "    Detector DQE%40t%9.1f%% %9.1f%%\n")
+		call pargr (100 * dqe)
+		call pargr (100 * cum)
+	}
+
+	call fprintf (fd, "\nStatistics per exposure per pixel:%40t%10s %10s\n")
+	    call pargstr ("e-")
+	    call pargstr ("sigma")
+	call fprintf (fd, "    Source%40t%10d %10.1f\n")
+	    call pargr (nobj)
+	    call pargr (dobj)
+	if (nbkg > 0.) {
+	    call fprintf (fd, "    Background%40t%10d %10.1f\n")
+		call pargr (nbkg)
+		call pargr (dbkg)
+	}
+	call fprintf (fd, "    Dark%40t%10d %10.1f\n")
+	    call pargr (ndark)
+	    call pargr (ddark)
+	call fprintf (fd, "    Readout%40t%10w %10.1f\n")
+	    call pargr (dreadout)
+	call fprintf (fd, "    Net%40t%10d %10.1f\n")
+	    call pargr (nobj + nbkg + ndark)
+	    call pargr (tnoise)
+
+	call fprintf (fd, "\nSignal-to-Noise Statistics:\n")
+	call fprintf (fd, "  %3d%7tper exposure per pixel\n")
+	    call pargr (snr)
+	    call pargr (ST_DISP(st,1))
+	call un_ctranr (ST_DUNANG(st), ST_DUN(st), ST_DISP(st,2), w, 1)
+	call fprintf (fd,
+	    "  %3d%7tper exposure per %.3g %s (%d pixel) resolution element\n")
+	    call pargr (snr * sqrt (real (ncols)))
+	    call pargr (w)
+	    call pargstr (ST_DUNITS(st))
+	    call pargi (ncols)
+	if (nexp > 1.) {
+	    call fprintf (fd,
+		"  %3d%10tper %d exposures per pixel\n")
+		call pargr (snr * sqnexp)
+		call pargi (nexp)
+		call pargr (ST_DISP(st,1))
+	    call un_ctranr (ST_DUNANG(st), ST_DUN(st), ST_DISP(st,2), w, 1)
+	    call fprintf (fd,
+	"  %3d%10tper %d exposures per %.3g %s (%d pixel) resolution element\n")
+		call pargr (snr * sqrt (real (ncols)) * sqnexp)
+		call pargi (nexp)
+		call pargr (w)
+		call pargstr (ST_DUNITS(st))
+		call pargi (ncols)
+	}
+
+	if (nexp > 1) {
+	    call fprintf (fd,
+		"\nExposure time: %d exposures of %.2fs")
+		call pargi (nexp)
+		call pargr (time)
+	    if (nexp * time > 60.) {
+		call fprintf (fd, " (%.1h)")
+		    call pargr (nexp * time / 3600.)
+	    } else {
+		call fprintf (fd, " (%.1f)")
+		call pargr (nexp * time)
+	    }
+	} else {
+	    call fprintf (fd, "\nExposure time: %.2fs")
+		call pargr (time)
+	    if (time > 60.) {
+		call fprintf (fd, " (%.1h)")
+		    call pargr (time / 3600.)
+	    }
+	}
+	if (ST_SKYSUB(st) == SKY_SHUFFLE)
+	    call fprintf (fd, " (shuffled)\n")
+	else
+	    call fprintf (fd, "\n")
+	call fprintf (fd, "\n")
+end
+
+
+# ST_CHECK -- Check for possible errors such as lack of proper order blocking.
+
+procedure st_check (st, fd, waves, npix)
+
+pointer st		#I SPECTIME structure
+int	fd		#I Output file descriptor
+real	waves[ARB]	#I Wavelengths
+int	npix		#I Number of pixels
+
+int	i, n, step
+real	nobj[4], nsky[4]
+real	w1, w2, w, dw, snr, thruput, sat, dnmax, maxcount, maxfrac
+pointer	tab
+real	stgetr(), tabgetr(), gr_getr()
+
+begin
+	tab = ST_TAB(st)
+
+	# Check for extrapolations.  This has to be done before the order
+	# overlap because that may reference wavelength which are extrapolated.
+
+	ifnoerr (w1 = tabgetr (tab, "spectrum", "", "table.xmin")) {
+	    w2 = tabgetr (tab, "spectrum", "", "table.xmax")
+	    if (w1 > waves[1] || w2 < waves[npix])
+		call fprintf (fd,
+		    "WARNING: Spectrum table extrapolated in wavelength\n")
+	}
+	ifnoerr (w1 = tabgetr (tab, "sky", "", "table.xmin")) {
+	    w2 = tabgetr (tab, "sky", "", "table.xmax")
+	    if (w1 > waves[1] || w2 < waves[npix])
+		call fprintf (fd,
+		    "WARNING: Sky table extrapolated in wavelength\n")
+	    ifnoerr (w1 = tabgetr (tab, "sky", "", "table.ymin")) {
+		w2 = tabgetr (tab, "sky", "", "table.ymax")
+		if (w1 > ST_PHASE(st) || w2 < ST_PHASE(st))
+		    call fprintf (fd,
+			"WARNING: Sky table extrapolated in lunar phase\n")
+	    }
+	}
+	ifnoerr (w1 = tabgetr (tab, "sensfunc", "", "table.xmin")) {
+	    w2 = tabgetr (tab, "sensfunc", "", "table.xmax")
+	    if (w1 > waves[1] || w2 < waves[npix])
+		call fprintf (fd,
+		    "WARNING: Sensitivity table extrapolated in wavelength\n")
+	}
+	ifnoerr (w1 = tabgetr (tab, "extinction", "", "table.xmin")) {
+	    w2 = tabgetr (tab, "extinction", "", "table.xmax")
+	    if (w1 > waves[1] || w2 < waves[npix])
+		call fprintf (fd,
+		    "WARNING: Extinction table extrapolated in wavelength\n")
+	}
+	ifnoerr (w1 = tabgetr (tab, "telescope", "", "table.xmin")) {
+	    w2 = tabgetr (tab, "telescope", "", "table.xmax")
+	    if (w1 > waves[1] || w2 < waves[npix])
+		call fprintf (fd,
+		    "WARNING: Telescope table extrapolated in wavelength\n")
+	}
+	ifnoerr (w1 = tabgetr (tab, "adc", "", "table.xmin")) {
+	    w2 = tabgetr (tab, "adc", "", "table.xmax")
+	    if (w1 > waves[1] || w2 < waves[npix])
+		call fprintf (fd,
+		    "WARNING: ADC table extrapolated in wavelength\n")
+	}
+	ifnoerr (w1 = tabgetr (tab, "filter", "", "table.xmin")) {
+	    w2 = tabgetr (tab, "filter", "", "table.xmax")
+	    if (w1 > waves[1] || w2 < waves[npix])
+		call fprintf (fd,
+		    "WARNING: Filter table extrapolated in wavelength\n")
+	}
+	ifnoerr (w1 = tabgetr (tab, "filter2", "", "table.xmin")) {
+	    w2 = tabgetr (tab, "filter2", "", "table.xmax")
+	    if (w1 > waves[1] || w2 < waves[npix])
+	    call fprintf (fd,
+		"WARNING: Second filter table extrapolated in wavelength\n")
+	}
+	ifnoerr (w1 = tabgetr (tab, "fiber", "", "table.xmin")) {
+	    w2 = tabgetr (tab, "fiber", "", "table.xmax")
+	    if (w1 > waves[1] || w2 < waves[npix])
+		call fprintf (fd,
+		    "WARNING: Fiber table extrapolated in wavelength\n")
+	}
+	ifnoerr (w1 = tabgetr (tab, "collimator", "", "table.xmin")) {
+	    w2 = tabgetr (tab, "collimator", "", "table.xmax")
+	    if (w1 > waves[1] || w2 < waves[npix])
+		call fprintf (fd,
+		    "WARNING: Collimator table extrapolated in wavelength\n")
+	}
+	ifnoerr (w1 = tabgetr (tab,"disperser","","table.xmin")/ST_ORDER(st,1)) {
+	    w2 = tabgetr (tab, "disperser", "", "table.xmax") / ST_ORDER(st,1)
+	    if (w1 > waves[1] || w2 < waves[npix])
+		call fprintf (fd,
+		    "WARNING: Disperser table extrapolated in wavelength\n")
+	    ifnoerr (w1 = tabgetr (tab, "disperser", "", "table.ymin")) {
+		w2 = tabgetr (tab, "disperser", "", "table.ymax")
+		if (w1 > ST_ORDER(st,1) || w2 < ST_ORDER(st,1))
+		    call fprintf (fd,
+			"WARNING: Disperser table extrapolated in order\n")
+	    }
+	}
+	ifnoerr (w1 = tabgetr (tab,"xdisperser", "","table.xmin")/ST_ORDER(st,2)) {
+	    w2 = tabgetr (tab, "xdisperser", "", "table.xmax") / ST_ORDER(st,2)
+	    if (w1 > waves[1] || w2 < waves[npix])
+		call fprintf (fd,
+		"WARNING: Cross disperser table extrapolated in wavelength\n")
+	    ifnoerr (w1 = tabgetr (tab, "xdisperser", "", "table.ymin")) {
+		w2 = tabgetr (tab, "xdisperser", "", "table.ymax")
+		if (w1 > ST_ORDER(st,2) || w2 < ST_ORDER(st,2))
+		    call fprintf (fd,
+		    "WARNING: Cross disperser table extrapolated in order\n")
+	    }
+	}
+	ifnoerr (w1 = tabgetr (tab, "corrector", "", "table.xmin")) {
+	    w2 = tabgetr (tab, "corrector", "", "table.xmax")
+	    if (w1 > waves[1] || w2 < waves[npix])
+		call fprintf (fd,
+		    "WARNING: Corrector table extrapolated in wavelength\n")
+	}
+	ifnoerr (w1 = tabgetr (tab, "camera", "", "table.xmin")) {
+	    w2 = tabgetr (tab, "camera", "", "table.xmax")
+	    if (w1 > waves[1] || w2 < waves[npix])
+		call fprintf (fd,
+		    "WARNING: Camera table extrapolated in wavelength\n")
+	}
+	ifnoerr (w1 = tabgetr (tab, "detector", "", "table.xmin")) {
+	    w2 = tabgetr (tab, "detector", "", "table.xmax")
+	    if (w1 > waves[1] || w2 < waves[npix])
+		call fprintf (fd,
+		    "WARNING: Detector table extrapolated in wavelength\n")
+	}
+
+	# Check for saturation and DN limits.
+	sat = stgetr (st, "saturation", "detector", MAX_REAL)
+	dnmax = stgetr (st, "dnmax", "detector", MAX_REAL)
+
+	# Check for insufficient sky pixels.
+	switch (ST_SKYSUB(st)) {
+	case SKY_LONGSLIT:
+	    if (ST_NSKYPIX(st) <= 0)
+		call fprintf (fd,
+		    "WARNING: Slit is too short for sky subtraction\n")
+	}
+
+	# Check for order overlaps and saturation.
+	if (npix > 3) {
+	    step = max (1, npix / 21)
+	    maxfrac = 0.
+	    maxcount = 0
+	    do i = 3*step, npix-3*step, step {
+		w = waves[i]
+		call st_snr (st, NULL, w, ST_NEXP(st), ST_TIME(st), nobj, nsky,
+		    snr, thruput)
+		maxcount = max (nobj[1]+nsky[3], maxcount)
+		if (nobj[1] > 0.) {
+		    maxfrac = max (nobj[2]/nobj[1], maxfrac)
+		    maxfrac = max (nobj[4]/nobj[1], maxfrac)
+		}
+	    }
+	} else {
+	    w = waves[1]
+	    call st_snr (st, NULL, w, ST_NEXP(st), ST_TIME(st), nobj, nsky,
+		snr, thruput)
+	    maxcount = max (nobj[1]+nsky[3], maxcount)
+	    if (nobj[1] > 0.) {
+		maxfrac = max (nobj[2]/nobj[1], maxfrac)
+		maxfrac = max (nobj[4]/nobj[1], maxfrac)
+	    }
+	}
+
+	if (maxcount > sat)
+	    call fprintf (fd, "WARNING: Exposure may saturate.\n")
+	if (maxcount / ST_GAIN(st) > dnmax)
+	    call fprintf (fd, "WARNING: Exposure may overflow DN maximum.\n")
+	if (maxfrac > 0.1)
+	    call fprintf (fd,
+		"WARNING: More than 10%% contribution from other orders.\n")
+
+	if (ST_DISPTYPE(st,2) != 0 && !IS_INDEFI(ST_ORDER(st,1))) {
+	    w = ST_CW(st)
+	    i = ST_ORDER(st,1)
+	    dw = abs (gr_getr (ST_GR(st,2), "dispersion"))
+	    dw = dw * ST_PIXSIZE(st) * ST_BIN(st,2)
+	    n = w * (1 - real (i) / real (i + 1)) / dw
+	    if (n < ST_APSIZE(st,2) / ST_SCALE(st,2)) 
+		call fprintf (fd, "WARNING: Orders overlap\n")
+	}
+
+	call fprintf (fd, "\n")
+end
+
+
+# ST_GTABLE -- Get table name from CL and load if a name is given.
+# Otherwise get table name from another table, if specified, and load if
+# a name is found.
+
+procedure st_gtable (st, name, table)
+
+pointer st		#I SPECTIME structure
+char	name[ARB]	#I Table to load (CL parameter name and table name)
+char	table[ARB]	#I Table to use if not defined by CL parameter
+
+pointer sp, dir, path, fname
+int	i, nowhite(), strdic()
+bool	streq()
+errchk	st_gtable1
+
+define	done_	10
+
+begin
+	call smark (sp)
+	call salloc (dir, SZ_FNAME, TY_CHAR)
+	call salloc (path, SZ_FNAME, TY_CHAR)
+	call salloc (fname, SZ_FNAME, TY_CHAR)
+
+	# Determine table name from CL parameter or table.
+	call stgstr (st, name, table, "", Memc[fname], SZ_FNAME)
+	i = nowhite (Memc[fname], Memc[fname], SZ_FNAME)
+
+	# Special cases.
+	if (streq (Memc[fname], "none"))
+	    goto done_
+	if (streq (name, "spectrum")) {
+	    ST_SPEC(st) = strdic (Memc[fname], Memc[path], SZ_FNAME, SPECTYPES)
+	    if (ST_SPEC(st) != SPEC_TAB && streq (Memc[fname], Memc[path]))
+		goto done_
+	    ST_SPEC(st) = SPEC_TAB
+	}
+
+	# If a filename is given find it and load it.
+	if (Memc[fname] != EOS)
+	    call st_gtable1 (st, name, Memc[fname])
+
+done_	call sfree (sp)
+end
+
+
+# ST_GTABLE1 -- Load table with search path.
+# Otherwise get table name from another table, if specified, and load if
+# a name is found.
+
+procedure st_gtable1 (st, name, fname)
+
+pointer st		#I SPECTIME structure
+char	name[ARB]	#I Table name
+char	fname[ARB]	#I File
+
+pointer sp, dir, path
+real	value
+int	i, ctor(), strlen(), clgfil(), access()
+errchk	tabload
+
+define	done_	10
+
+begin
+	call smark (sp)
+	call salloc (dir, SZ_FNAME, TY_CHAR)
+	call salloc (path, SZ_FNAME, TY_CHAR)
+
+	# Check for constant value.
+	i = 1
+	if (ctor (fname, i, value) == strlen (fname)) {
+	    call tabload (ST_TAB(st), name, fname)
+	    goto done_
+	}
+
+	# Absolute path or current directory.
+	if (access (fname, READ_ONLY, 0) == YES) {
+	    call tabload (ST_TAB(st), name, fname)
+	    goto done_
+	}
+
+	# Search directories.
+	call clprew (ST_SEARCH(st))
+	while (clgfil (ST_SEARCH(st), Memc[dir], SZ_FNAME) != EOF) {
+	    call sprintf (Memc[path], SZ_FNAME, "%s/%s")
+		call pargstr (Memc[dir])
+		call pargstr (fname)
+	    if (access (Memc[path], READ_ONLY, 0) == NO)
+		next
+	    call tabload (ST_TAB(st), name, Memc[path])
+	    goto done_
+	}
+
+	# Search sptimelib$.
+	call sprintf (Memc[path], SZ_FNAME, "%s/%s")
+	    call pargstr ("sptimelib$")
+	    call pargstr (fname)
+	if (access (Memc[path], READ_ONLY, 0) == YES) {
+	    call tabload (ST_TAB(st), name, Memc[path])
+	    goto done_
+	}
+
+	call sprintf (Memc[path], SZ_FNAME, "Table `%s' not found")
+	    call pargstr (fname)
+	call error (1, Memc[path])
+
+done_	call sfree (sp)
+end
+
+
+# ST_SPECTRUM -- Return spectrum flux.
+
+real procedure st_spectrum (st, wave)
+
+pointer st		#I SPECTIME pointer
+real	wave		#I Wavelength
+real	flux		#O Flux
+
+real	tabinterp1(), ccm()
+errchk	tabinterp1, ccm
+
+begin
+	switch (ST_SPEC(st)) {
+	case SPEC_TAB:
+	    flux = tabinterp1 (ST_TAB(st), "spectrum", wave)
+	case SPEC_BB:
+	    flux = (exp (min (30., 1.4338e8/(ST_REFW(st)*ST_PARAM(st)))) - 1) /
+		(exp (min (30., 1.4338e8/(wave*ST_PARAM(st)))) - 1)
+	    flux = ST_REFFL(st) * (ST_REFW(st) / wave) ** 5 * flux
+	case SPEC_FL:
+	    flux = ST_REFFL(st) * (wave/ST_REFW(st)) ** ST_PARAM(st)
+	case SPEC_FN:
+	    flux = ST_REFFL(st) * (wave/ST_REFW(st)) ** (-2.-ST_PARAM(st))
+	}
+	flux = flux * 10. ** (0.4 * ST_AV(st) *
+	    (ccm (ST_REFW(st), ST_RV(st)) - ccm (wave, ST_RV(st))))
+
+	return (flux)
+end
+
+
+# ST_OUTPUT -- Output graphs and/or lists.
+
+procedure st_output (st, gp, fd, interactive, output, w, npts)
+
+pointer st			#I SPECTIME structure
+pointer gp			#U GIO pointer
+int	fd			#I FIO pointer
+bool	interactive		#I Interactive?
+char	output[ARB]		#I Output type
+real	w[npts]			#I Wavelengths
+int	npts			#I Number of points
+
+int	i, j, ndata, type
+real	nobj[4], nsky[4]
+real	wx1, wx2, wy1, wy2, wmin, wmax, a, b, buf, f, snr, thruput, airmass
+pointer sp, spec, name
+pointer	title, xlabel, ylabel, id[4], xdata, ydata, str, tab
+
+bool	tabexists()
+int	strdic(), clgcur()
+real	st_dispeff(), tabinterp1(), tabinterp2(), st_spectrum(), tabgetr()
+errchk	st_snr, st_dispeff, tabinterp1, tabinterp2, st_spectrum
+
+begin
+	if (gp == NULL && fd == NULL)
+	    return
+
+	call smark (sp)
+	call salloc (spec, SZ_LINE, TY_CHAR)
+	call salloc (name, SZ_LINE, TY_CHAR)
+	call salloc (title, SZ_LINE, TY_CHAR)
+	call salloc (xlabel, SZ_LINE, TY_CHAR)
+	call salloc (ylabel, SZ_LINE, TY_CHAR)
+	call salloc (id[1], SZ_LINE, TY_CHAR)
+	call salloc (id[2], SZ_LINE, TY_CHAR)
+	call salloc (id[3], SZ_LINE, TY_CHAR)
+	call salloc (id[4], SZ_LINE, TY_CHAR)
+	call salloc (xdata, npts, TY_REAL)
+	call salloc (ydata, 4*npts, TY_REAL)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	tab = ST_TAB(st)
+
+	# Set title, data, and data limits.
+	call st_gtitle (st, "title", "spectrograph", Memc[title], SZ_LINE)
+	call st_gtitle (st, "spectitle", "spectrum", Memc[spec], SZ_LINE)
+	if (Memc[spec] == EOS) {
+	    switch (ST_SPEC(st)) {
+	    case SPEC_BB:
+		call sprintf (Memc[spec], SZ_LINE,
+		    "Blackbody spectrum of temperature %g K")
+		    call pargr (ST_PARAM(st))
+	    case SPEC_FL:
+		call sprintf (Memc[spec], SZ_LINE,
+		    "F(lambda) power law spectrum of index %g")
+		    call pargr (ST_PARAM(st))
+	    case SPEC_FN:
+		call sprintf (Memc[spec], SZ_LINE,
+		    "F(nu) power law spectrum of index %g")
+		    call pargr (ST_PARAM(st))
+	    }
+	}
+
+	call sprintf (Memc[xlabel], SZ_LINE, "Dispersion (%s)")
+	    call pargstr (ST_DUNITS(st))
+
+	ndata = npts
+	#call amovr (w, Memr[xdata], ndata)
+	call un_ctranr (ST_DUNANG(st), ST_DUN(st), w, Memr[xdata], ndata)
+	Memr[ydata] = INDEF
+	Memr[ydata+npts] = INDEF
+	Memr[ydata+2*npts] = INDEF
+	Memr[ydata+3*npts] = INDEF
+	#wx1 = INDEF; wx2 = INDEF; wy1 = -4.; wy2 = 104.
+	wx1 = INDEF; wx2 = INDEF; wy1 = INDEF; wy2 = INDEF
+
+	type = strdic (output, Memc[name], SZ_LINE, OUTTYPES)
+	switch (type) {
+	case OUT_RATE:
+	    if (Memc[title] != EOS)
+		call strcat ("\n", Memc[title], SZ_LINE)
+	    call strcat (Memc[spec], Memc[title], SZ_LINE)
+	    call strcpy ("Photons/s/A", Memc[ylabel], SZ_LINE)
+	    call strcpy ("Object (ph/s/A)", Memc[id[1]], SZ_LINE) 
+	    call strcpy ("Background (ph/s/A)", Memc[id[2]], SZ_LINE) 
+	    call strcpy ("Total (ph/s/A)", Memc[id[3]], SZ_LINE) 
+
+	    do i = 1, npts {
+		call st_snr (st, NULL, w[i], ST_NEXP(st), ST_TIME(st),
+		    nobj, nsky, snr, thruput)
+		Memr[ydata+i-1] = nobj[1] / ST_TIME(st) / ST_DISP(st,1)
+		Memr[ydata+npts+i-1] = nsky[1] / ST_TIME(st) / ST_DISP(st,1)
+		Memr[ydata+2*npts+i-1] = (nobj[1]+nsky[1]) /
+		    ST_TIME(st) / ST_DISP(st,1)
+	    }
+	    call alimr (Memr[ydata+npts], npts, a, b)
+	    if (b <= 0.) {
+		Memr[ydata+npts] = INDEF
+		Memr[ydata+2*npts] = INDEF
+		if (Memc[title] != EOS)
+		    call strcat ("\n", Memc[title], SZ_LINE)
+		call sprintf (Memc[str], SZ_LINE, "Object photon rates\n")
+		call strcat (Memc[str], Memc[title], SZ_LINE)
+	    } else {
+		if (Memc[title] != EOS)
+		    call strcat ("\n", Memc[title], SZ_LINE)
+		call sprintf (Memc[str], SZ_LINE,
+		    "Object, background, and total photon rates\n")
+		call strcat (Memc[str], Memc[title], SZ_LINE)
+	    }
+
+	    wy1 = INDEF; wy2 = INDEF
+	case OUT_OBJ:
+	    if (Memc[title] != EOS)
+		call strcat ("\n", Memc[title], SZ_LINE)
+	    call strcat (Memc[spec], Memc[title], SZ_LINE)
+	    call sprintf (Memc[str], SZ_LINE, "\nExposure time = %d seconds")
+		call pargr (ST_NEXP(st) * max (ST_TIME(st),1.))
+	    call strcat (Memc[str], Memc[title], SZ_LINE)
+	    if (ST_NEXP(st) > 1) {
+		call sprintf (Memc[str], SZ_LINE, "(%d exposures)")
+		    call pargi (ST_NEXP(st))
+		call strcat (Memc[str], Memc[title], SZ_LINE)
+	    }
+	    call strcpy ("Object DN", Memc[ylabel], SZ_LINE)
+	    call strcpy ("Total Object (counts)", Memc[id[1]], SZ_LINE) 
+	    call sprintf (Memc[id[2]], SZ_LINE, "Order %d (counts)")
+		call pargi (ST_ORDER(st,1)-1)
+	    call sprintf (Memc[id[3]], SZ_LINE, "Order %d (counts)")
+		call pargi (ST_ORDER(st,1))
+	    call sprintf (Memc[id[4]], SZ_LINE, "Order %d (counts)")
+		call pargi (ST_ORDER(st,1)+1)
+
+	    do i = 1, npts {
+		call st_snr (st, NULL, w[i], 1, ST_NEXP(st) * ST_TIME(st),
+		    nobj, nsky, snr, thruput)
+		Memr[ydata+i-1] = nobj[1] / ST_GAIN(st)
+		Memr[ydata+npts+i-1] = nobj[2] / ST_GAIN(st)
+		Memr[ydata+2*npts+i-1] = nobj[3] / ST_GAIN(st)
+		Memr[ydata+3*npts+i-1] = nobj[4] / ST_GAIN(st)
+	    }
+	    call alimr (Memr[ydata+npts], npts, a, b)
+	    if (b <= 0.)
+		Memr[ydata+npts] = INDEF
+	    call alimr (Memr[ydata+3*npts], npts, a, b)
+	    if (b <= 0.)
+		Memr[ydata+3*npts] = INDEF
+	    if (IS_INDEF(Memr[ydata+npts]) && IS_INDEF(Memr[ydata+3*npts])) {
+		Memr[ydata+2*npts] = INDEF
+		if (Memc[title] != EOS)
+		    call strcat ("\n", Memc[title], SZ_LINE)
+		call sprintf (Memc[str], SZ_LINE,
+		    "Object DN at gain of %.2g\n")
+		    call pargr (ST_GAIN(st))
+		call strcat (Memc[str], Memc[title], SZ_LINE)
+	    } else {
+		if (Memc[title] != EOS)
+		    call strcat ("\n", Memc[title], SZ_LINE)
+		call sprintf (Memc[str], SZ_LINE,
+		    "Object DN at gain of %.2g with order contributions\n")
+		    call pargr (ST_GAIN(st))
+		call strcat (Memc[str], Memc[title], SZ_LINE)
+	    }
+
+	    wy1 = INDEF; wy2 = INDEF
+	case OUT_SNR:
+	    if (Memc[title] != EOS)
+		call strcat ("\n", Memc[title], SZ_LINE)
+	    call strcat (Memc[spec], Memc[title], SZ_LINE)
+	    call sprintf (Memc[str], SZ_LINE, "\nExposure time = %d seconds")
+		call pargr (ST_NEXP(st) * max(ST_TIME(st),1.))
+	    call strcat (Memc[str], Memc[title], SZ_LINE)
+	    if (ST_NEXP(st) > 1) {
+		call sprintf (Memc[str], SZ_LINE, "(%d exposures)")
+		    call pargi (ST_NEXP(st))
+		call strcat (Memc[str], Memc[title], SZ_LINE)
+	    }
+	    call strcpy ("S/N", Memc[ylabel], SZ_LINE)
+	    call strcpy ("S/N", Memc[id[1]], SZ_LINE) 
+
+	    a = sqrt (real(ST_NEXP(st)))
+	    do i = 1, npts {
+		call st_snr (st, NULL, w[i], ST_NEXP(st), ST_TIME(st),
+		    nobj, nsky, snr, thruput)
+		Memr[ydata+i-1] = a * snr
+	    }
+	    wy1 = INDEF; wy2 = INDEF
+	case OUT_ATM:
+	    call st_gtitle (st, "exttitle", "extinction", Memc[str], SZ_LINE)
+	    if (Memc[str] == EOS)
+		call strcpy ("Extinction", Memc[str], SZ_LINE)
+	    if (Memc[title] != EOS)
+		call strcat ("\n", Memc[title], SZ_LINE)
+	    call strcat (Memc[str], Memc[title], SZ_LINE)
+	    call sprintf (Memc[str], SZ_LINE, "\nAirmass of %g")
+		call pargr (ST_AIRMASS(st))
+	    call strcat (Memc[str], Memc[title], SZ_LINE)
+	    call strcpy ("Transmission (%)", Memc[ylabel], SZ_LINE)
+	    call strcpy ("ATM Transmission (%)", Memc[id[1]], SZ_LINE) 
+
+	    iferr  {
+		do i = 1, npts {
+		    f = tabinterp1 (tab, "extinction", w[i])
+		    Memr[ydata+i-1] = 100 * 10 ** (-0.4 * f * ST_AIRMASS(st))
+		}
+	    } then
+		call amovkr (100., Memr[ydata], npts)
+	case OUT_TEL:
+	    call st_gtitle (st, "teltitle", "telescope", Memc[str], SZ_LINE)
+	    if (Memc[str] == EOS)
+		call strcpy ("Telescope", Memc[str], SZ_LINE)
+	    if (Memc[title] != EOS)
+		call strcat ("\n", Memc[title], SZ_LINE)
+	    call strcat (Memc[str], Memc[title], SZ_LINE)
+	    call strcpy ("Transmission (%)", Memc[ylabel], SZ_LINE)
+	    call strcpy ("Telescope Transmission (%)", Memc[id[1]], SZ_LINE) 
+
+	    iferr  {
+		do i = 1, npts
+		    Memr[ydata+i-1] = 100 *
+			tabinterp1 (tab, Memc[name], w[i])
+	    } then
+		call amovkr (100., Memr[ydata], npts)
+	case OUT_AP:
+	    call st_gtitle (st, "aptitle", "aperture", Memc[str], SZ_LINE)
+	    if (Memc[str] == EOS)
+		call strcpy ("Aperture", Memc[str], SZ_LINE)
+	    if (Memc[title] != EOS)
+		call strcat ("\n", Memc[title], SZ_LINE)
+	    call strcat (Memc[str], Memc[title], SZ_LINE)
+	    call sprintf (Memc[str], SZ_LINE, "\nSeeing of %.2f\"")
+		call pargr (ST_SEEING(st))
+	    call strcat (Memc[str], Memc[title], SZ_LINE)
+	    call strcpy ("Transmission (%)", Memc[ylabel], SZ_LINE)
+	    call strcpy ("Aperture (%)", Memc[id[1]], SZ_LINE) 
+
+	    iferr  {
+		switch (ST_APTYPE(st)) {
+		case CIRCULAR:
+		    f = 100 * tabinterp1 (tab, Memc[name],
+			ST_APSIZE(st,1) / ST_SEEING(st))
+		case RECTANGULAR:
+		    f = 100 * tabinterp2 (tab, Memc[name],
+			ST_APSIZE(st,1) / ST_SEEING(st),
+			ST_APSIZE(st,2) / ST_SEEING(st))
+		}
+	    } then
+		f = 100
+	    call amovkr (f, Memr[ydata], npts)
+	case OUT_FIB:
+	    if (tabexists (tab, Memc[name])) {
+		call st_gtitle (st, "fibtitle", "fiber", Memc[str], SZ_LINE)
+		if (Memc[str] == EOS)
+		    call strcpy ("Fiber", Memc[str], SZ_FNAME)
+		if (Memc[title] != EOS)
+		    call strcat ("\n", Memc[title], SZ_LINE)
+		call strcat (Memc[str], Memc[title], SZ_LINE)
+		call strcpy ("Transmission (%)", Memc[ylabel], SZ_LINE)
+
+		iferr  {
+		    do i = 1, npts
+			Memr[ydata+i-1] = 100 *
+			    tabinterp1 (tab, Memc[name], w[i])
+		} then
+		    call amovkr (100., Memr[ydata], npts)
+	    }
+	case OUT_FILT, OUT_FILT2:
+	    if (tabexists (tab, Memc[name])) {
+		if (type == OUT_FILT)
+		    call st_gtitle (st, "ftitle", "filter", Memc[str], SZ_LINE)
+		else
+		    call st_gtitle (st, "f2title", "filter2", Memc[str],SZ_LINE)
+		if (Memc[str] == EOS)
+		    call strcpy ("Filter", Memc[str], SZ_LINE)
+		if (Memc[title] != EOS)
+		    call strcat ("\n", Memc[title], SZ_LINE)
+		call strcat (Memc[str], Memc[title], SZ_LINE)
+		call strcpy ("Transmission (%)", Memc[ylabel], SZ_LINE)
+		call strcpy ("Fiber (%)", Memc[id[1]], SZ_LINE) 
+
+		iferr  {
+		    do i = 1, npts
+			Memr[ydata+i-1] = 100 *
+			    tabinterp1 (tab, Memc[name], w[i])
+		} then
+		    call amovkr (100., Memr[ydata], npts)
+	    }
+	case OUT_COL:
+	    call st_gtitle (st, "coltitle", "collimator", Memc[str], SZ_LINE)
+	    if (Memc[str] == EOS)
+		call strcpy ("Collimator", Memc[str], SZ_LINE)
+	    if (Memc[title] != EOS)
+		call strcat ("\n", Memc[title], SZ_LINE)
+	    call strcat (Memc[str], Memc[title], SZ_LINE)
+	    call strcpy ("Transmission (%)", Memc[ylabel], SZ_LINE)
+	    call strcpy ("Collimator (%)", Memc[id[1]], SZ_LINE) 
+
+	    iferr  {
+		do i = 1, npts
+		    Memr[ydata+i-1] = 100 *
+			tabinterp1 (tab, Memc[name], w[i])
+	    } then
+		call amovkr (100., Memr[ydata], npts)
+	case OUT_DISP:
+	    if (ST_DISPTYPE(st,1) != 0) {
+		call st_gtitle (st, "disptitle", "disperser", Memc[str],SZ_LINE)
+		if (Memc[str] == EOS)
+		    call strcpy ("Disperser", Memc[str], SZ_LINE)
+		if (Memc[title] != EOS)
+		    call strcat ("\n", Memc[title], SZ_LINE)
+		call strcat (Memc[str], Memc[title], SZ_LINE)
+		call strcpy ("Efficiency (%)", Memc[ylabel], SZ_LINE)
+		call strcpy ("Disperser (%)", Memc[id[1]], SZ_LINE) 
+
+		do i = 1, npts
+		    Memr[ydata+i-1] = 100 * st_dispeff (st, Memc[name], w[i],
+			ST_ORDER(st,1))
+	    }
+	case OUT_XDISP:
+	    if (ST_DISPTYPE(st,2) != 0) {
+		call st_gtitle (st, "xdisptitle","xdisperser",Memc[str],SZ_LINE)
+		if (Memc[str] == EOS)
+		    call strcpy ("Crossdisperser", Memc[str], SZ_LINE)
+		if (Memc[title] != EOS)
+		    call strcat ("\n", Memc[title], SZ_LINE)
+		call strcat (Memc[str], Memc[title], SZ_LINE)
+		call strcpy ("Efficiency (%)", Memc[ylabel], SZ_LINE)
+		call strcpy ("Cross-disperser (%)", Memc[id[1]], SZ_LINE) 
+
+		do i = 1, npts
+		    Memr[ydata+i-1] = 100 * st_dispeff (st, Memc[name], w[i],
+			ST_ORDER(st,2))
+	    }
+	case OUT_COR:
+	    if (tabexists (tab, Memc[name])) {
+		call st_gtitle (st, "cortitle", "corrector", Memc[str], SZ_LINE)
+		if (Memc[str] == EOS)
+		    call strcpy ("Corrector", Memc[str], SZ_LINE)
+		if (Memc[title] != EOS)
+		    call strcat ("\n", Memc[title], SZ_LINE)
+		call strcat (Memc[str], Memc[title], SZ_LINE)
+		call strcpy ("Transmission (%)", Memc[ylabel], SZ_LINE)
+		call strcpy ("Corrector (%)", Memc[id[1]], SZ_LINE) 
+
+		iferr  {
+		    do i = 1, npts
+			Memr[ydata+i-1] = 100 *
+			    tabinterp1 (tab, Memc[name], w[i])
+		} then
+		    call amovkr (100., Memr[ydata], npts)
+	    }
+	case OUT_CAM:
+	    call st_gtitle (st, "camtitle", "camera", Memc[str], SZ_LINE)
+	    if (Memc[str] == EOS)
+		call strcpy ("Camera", Memc[str], SZ_LINE)
+	    if (Memc[title] != EOS)
+		call strcat ("\n", Memc[title], SZ_LINE)
+	    call strcat (Memc[str], Memc[title], SZ_LINE)
+	    call strcpy ("Transmission (%)", Memc[ylabel], SZ_LINE)
+	    call strcpy ("Camera (%)", Memc[id[1]], SZ_LINE) 
+
+	    iferr  {
+		do i = 1, npts
+		    Memr[ydata+i-1] = 100 *
+			tabinterp1 (tab, Memc[name], w[i])
+	    } then
+		call amovkr (100., Memr[ydata], npts)
+	case OUT_DET:
+	    call st_gtitle (st, "dettitle", "detector", Memc[str], SZ_LINE)
+	    if (Memc[str] == EOS)
+		call strcpy ("Detector", Memc[str], SZ_LINE)
+	    if (Memc[title] != EOS)
+		call strcat ("\n", Memc[title], SZ_LINE)
+	    call strcat (Memc[str], Memc[title], SZ_LINE)
+	    call strcpy ("DQE (%)", Memc[ylabel], SZ_LINE)
+	    call strcpy ("DQE (%)", Memc[id[1]], SZ_LINE) 
+
+	    iferr  {
+		do i = 1, npts
+		    Memr[ydata+i-1] = 100 *
+			tabinterp1 (tab, Memc[name], w[i])
+	    } then
+		call amovkr (100., Memr[ydata], npts)
+	case OUT_SPEC:
+	    call strcpy ("Spectrograph Other", Memc[str], SZ_LINE)
+	    if (Memc[title] != EOS)
+		call strcat ("\n", Memc[title], SZ_LINE)
+	    call strcat (Memc[str], Memc[title], SZ_LINE)
+	    call strcpy ("Transmission (%)", Memc[ylabel], SZ_LINE)
+	    call strcpy ("Spectrograph Other (%s)", Memc[id[1]], SZ_LINE) 
+
+	    iferr  {
+		do i = 1, npts
+		    Memr[ydata+i-1] = 100 *
+			tabinterp1 (tab, Memc[name], w[i])
+	    } then
+		Memr[ydata] = INDEF
+	case OUT_ADC:
+	    if (tabexists (tab, Memc[name])) {
+		call st_gtitle (st, "adctitle", "adc", Memc[str], SZ_LINE)
+		if (Memc[str] == EOS)
+		    call strcpy ("ADC", Memc[str], SZ_LINE)
+		if (Memc[title] != EOS)
+		    call strcat ("\n", Memc[title], SZ_LINE)
+		call strcat (Memc[str], Memc[title], SZ_LINE)
+		call strcpy ("Transmission (%)", Memc[ylabel], SZ_LINE)
+		call strcpy ("ADC (%)", Memc[id[1]], SZ_LINE) 
+
+		iferr  {
+		    do i = 1, npts
+			Memr[ydata+i-1] = 100 *
+			    tabinterp1 (tab, Memc[name], w[i])
+		} then
+		    call amovkr (100., Memr[ydata], npts)
+	    }
+	case OUT_EMIS:
+	    if (tabexists (tab, Memc[name]) && ST_THERMAL(st) > 0.) {
+		if (Memc[title] != EOS)
+		    call strcat ("\n", Memc[title], SZ_LINE)
+		call st_gtitle (st, "emistitle", "emissivity",Memc[str],SZ_LINE)
+		if (Memc[str] == EOS)
+		    call strcpy ("Emissivity", Memc[str], SZ_LINE)
+		call strcat (Memc[str], Memc[title], SZ_LINE)
+		call sprintf (Memc[str], SZ_LINE, "\nTemperature  = %.1f K")
+		    call pargr (ST_THERMAL(st))
+		call strcat (Memc[str], Memc[title], SZ_LINE)
+		call strcpy ("Emissivity (%)", Memc[ylabel], SZ_LINE)
+		call strcpy ("Emissivity (%)", Memc[id[1]], SZ_LINE) 
+
+		iferr  {
+		    do i = 1, npts
+			Memr[ydata+i-1] = 100 *
+			    tabinterp1 (tab, Memc[name], w[i])
+		} then
+		    call amovkr (100., Memr[ydata], npts)
+	    }
+	case OUT_THRU:
+	    if (Memc[title] != EOS)
+		call strcat ("\n", Memc[title], SZ_LINE)
+	    call strcat ("System Thruput (Telescope to Detected Photons)",
+		Memc[title], SZ_LINE)
+	    call strcpy ("Thruput (%)", Memc[ylabel], SZ_LINE)
+	    call strcpy ("Thruput (%)", Memc[id[1]], SZ_LINE) 
+
+	    do i = 1, npts {
+		call st_snr (st, NULL, w[i], ST_NEXP(st), ST_TIME(st),
+		    nobj, nsky, snr, thruput)
+		Memr[ydata+i-1] = 100 * thruput
+	    }
+	case OUT_COUNTS:
+	    if (Memc[title] != EOS)
+		call strcat ("\n", Memc[title], SZ_LINE)
+	    call strcat (Memc[spec], Memc[title], SZ_LINE)
+	    call sprintf (Memc[str], SZ_LINE, "\nExposure time = %d seconds")
+		call pargr (ST_NEXP(st) * max(ST_TIME(st),1.))
+	    call strcat (Memc[str], Memc[title], SZ_LINE)
+	    if (ST_NEXP(st) > 1) {
+		call sprintf (Memc[str], SZ_LINE, "(%d exposures)")
+		    call pargi (ST_NEXP(st))
+		call strcat (Memc[str], Memc[title], SZ_LINE)
+	    }
+	    call strcpy ("DN", Memc[ylabel], SZ_LINE)
+	    call strcpy ("Object (counts)", Memc[id[1]], SZ_LINE) 
+	    call strcpy ("Background (counts)", Memc[id[2]], SZ_LINE) 
+	    call strcpy ("Total (counts)", Memc[id[3]], SZ_LINE) 
+
+	    do i = 1, npts {
+		call st_snr (st, NULL, w[i], 1, ST_NEXP(st) * ST_TIME(st),
+		    nobj, nsky, snr, thruput)
+		Memr[ydata+i-1] = nobj[1] / ST_GAIN(st)
+		Memr[ydata+npts+i-1] = nsky[1] / ST_GAIN(st)
+		Memr[ydata+2*npts+i-1] = (nobj[1]+nsky[1]) / ST_GAIN(st)
+	    }
+	    call alimr (Memr[ydata+npts], npts, a, b)
+	    if (b <= 0.) {
+		Memr[ydata+npts] = INDEF
+		Memr[ydata+2*npts] = INDEF
+		if (Memc[title] != EOS)
+		    call strcat ("\n", Memc[title], SZ_LINE)
+		call sprintf (Memc[str], SZ_LINE,
+		    "Object DN at gain of %.2g\n")
+		    call pargr (ST_GAIN(st))
+		call strcat (Memc[str], Memc[title], SZ_LINE)
+	    } else {
+		if (Memc[title] != EOS)
+		    call strcat ("\n", Memc[title], SZ_LINE)
+		call sprintf (Memc[str], SZ_LINE,
+		    "Object, background, and total DN at gain of %.2g\n")
+		    call pargr (ST_GAIN(st))
+		call strcat (Memc[str], Memc[title], SZ_LINE)
+	    }
+
+	    wy1 = INDEF; wy2 = INDEF
+	case OUT_SENS:
+	    airmass = ST_AIRMASS(st)
+	    ST_AIRMASS(st) = 0.
+
+	    if (Memc[title] != EOS)
+		call strcat ("\n", Memc[title], SZ_LINE)
+	    call strcat ("Sensitivity Function", Memc[title], SZ_LINE)
+	    call strcpy ("Magnitudes", Memc[ylabel], SZ_LINE)
+	    call strcpy ("Sensitivity (mag)", Memc[id[1]], SZ_LINE) 
+
+	    wy1 = MAX_REAL
+	    wy2 = -MAX_REAL
+	    do i = 1, npts {
+		call st_snr (st, NULL, w[i], ST_NEXP(st), ST_TIME(st),
+		    nobj, nsky, snr, thruput)
+		a = nobj[1] / ST_GAIN(st) / ST_TIME(st) / ST_DISP(st,1)
+		b = st_spectrum (st, w[i])
+		if (a > 0. && b > 0.) {
+		    a = 2.5 * (log10 (a) - log10 (b))
+		    wy1 = min (wy1, a)
+		    wy2 = max (wy2, a)
+		} else
+		    a = 0
+		Memr[ydata+i-1] = a
+	    }
+	    if (wy1 < wy2) {
+		buf = 0.1 * (wy2 - wy1)
+		wy1 = wy1 - buf
+		wy2 = wy2 + buf
+	    } else {
+		wy1 = INDEF
+		wy2 = INDEF
+	    }
+
+	    ST_AIRMASS(st) = airmass
+	case OUT_CORRECT:
+	    if (tabexists (tab, "sensfunc")) {
+		airmass = ST_AIRMASS(st)
+		ST_AIRMASS(st) = 0.
+
+		iferr (call tabgstr (tab, "sensfunc", "",
+		    "title", Memc[str], SZ_LINE)) {
+		    call tabgstr (tab, "sensfunc", "", "table.filename",
+			Memc[str], SZ_LINE)
+		}
+		if (Memc[title] != EOS)
+		    call strcat ("\n", Memc[title], SZ_LINE)
+		call strcat ("Correction from: ", Memc[title], SZ_LINE)
+		call strcat (Memc[str], Memc[title], SZ_LINE)
+		call strcpy ("Correction factor", Memc[ylabel], SZ_LINE)
+		call strcpy ("Correction factor", Memc[id[1]], SZ_LINE) 
+
+		wmin = tabgetr (tab, "sensfunc", "", "table.xmin")
+		wmax = tabgetr (tab, "sensfunc", "", "table.xmax")
+
+		wy1 = MAX_REAL
+		wy2 = -MAX_REAL
+		ndata = 0
+		do i = 1, npts {
+		    if (w[i] < wmin || w[i] > wmax)
+			next
+		    call st_snr (st, NULL, w[i], ST_NEXP(st), ST_TIME(st),
+			nobj, nsky, snr, thruput)
+		    a = nobj[1] / ST_GAIN(st) / ST_TIME(st) / ST_DISP(st,1)
+		    b = st_spectrum (st, w[i])
+		    if (a <= 0. || b <= 0.)
+			next
+		    a = 2.5 * (log10 (a) - log10 (b))
+		    b = tabinterp1 (tab, "sensfunc", w[i])
+		    a = 10. ** (0.4 * (b - a))
+		    wy1 = min (wy1, a)
+		    wy2 = max (wy2, a)
+		    Memr[xdata+ndata] = w[i]
+		    Memr[ydata+ndata] = a
+		    ndata = ndata + 1
+		}
+		if (wy1 < wy2) {
+		    buf = 0.1 * max (0.1, wy2 - wy1)
+		    wy1 = wy1 - buf
+		    wy2 = wy2 + buf
+		} else {
+		    wy1 = INDEF
+		    wy2 = INDEF
+		}
+
+		ST_AIRMASS(st) = airmass
+	    }
+	}
+
+	# Draw graph.
+	if (gp != NULL && !IS_INDEF(Memr[ydata])) {
+	    call greactivate (gp, 0)
+	    if (IS_INDEF(wx1) || IS_INDEF(wx2)) {
+		call alimr (Memr[xdata], ndata, wx1, wx2)
+		buf = 0.1 * (wx2 - wx1)
+		wx1 = wx1 - buf
+		wx2 = wx2 + buf
+	    }
+	    if (IS_INDEF(wy1) || IS_INDEF(wy2)) {
+		call alimr (Memr[ydata], ndata, wy1, wy2)
+		do i = 1, 3 {
+		    if (!IS_INDEF(Memr[ydata+i*npts])) { 
+			call alimr (Memr[ydata+i*npts], ndata, a, b)
+			wy1 = min (wy1, a)
+			wy2 = max (wy2, b)
+		    }
+		}
+		buf = 0.1 * (wy2 - wy1)
+		wy1 = wy1 - buf
+		wy2 = wy2 + buf
+	    }
+
+	    call gclear (gp)
+	    call gsview (gp, 0.15, 0.95, 0.15, 0.85)
+	    call gswind (gp, wx1, wx2, wy1, wy2)
+	    call glabax (gp, Memc[title], Memc[xlabel], Memc[ylabel])
+	    do i = 1, 4 {
+		if (!IS_INDEF(Memr[ydata+(i-1)*npts])) { 
+		    call gseti (gp, G_PLTYPE, i)
+		    call gseti (gp, G_PLCOLOR, i)
+		    call gpline (gp, Memr[xdata], Memr[ydata+(i-1)*npts], ndata)
+		}
+	    }
+
+	    if (interactive) {
+		call printf (
+		    "Type 'q' to quit graphs or any other key to continue")
+		i = clgcur ("gcur", a, b, i, j, Memc[str], SZ_LINE)
+		if (j == 'q')
+		    call gclose (gp)
+	    }
+	    if (gp != NULL)
+		call gdeactivate (gp, 0)
+	}
+
+	# Output list.
+	if (fd != NULL && !IS_INDEF(Memr[ydata])) {
+	    j = 0
+	    do i = 0, ARB {
+		if (Memc[title+i] == EOS || Memc[title+i] == '\n') {
+		    if (j != 0) {
+			Memc[str+j] = EOS
+			call fprintf (fd, "# %s\n")
+			    call pargstr (Memc[str])
+		    }
+		    if (Memc[title+i] == EOS)
+			break
+		    j = 0
+		} else {
+		    Memc[str+j] = Memc[title+i]
+		    j = j + 1
+		}
+	    }
+	    call fprintf (fd, "# Column 1: %s\n")
+		call pargstr (Memc[xlabel])
+	    do j = 1, 4 {
+		if (IS_INDEF(Memr[ydata+(j-1)*npts]))
+		    next
+		call fprintf (fd, "# Column %d: %s\n")
+		    call pargi (j+1)
+		    call pargstr (Memc[id[j]])
+	    }
+	    do i = 1, ndata {
+		call fprintf (fd, "%12.8g")
+		    call pargr (Memr[xdata+i-1])
+		do j = 1, 4 {
+		    if (IS_INDEF(Memr[ydata+(j-1)*npts]))
+			next
+		    call fprintf (fd, "  %12.8g")
+			call pargr (Memr[ydata+(j-1)*npts+i-1])
+		}
+		call fprintf (fd, "\n")
+	    }
+	    call fprintf (fd, "\n")
+	}
+
+	call sfree (sp)
+end
+
+
+# CCM -- Compute CCM Extinction Law
+
+real procedure ccm (wavelength, rv) 
+
+real	wavelength		# Wavelength in Angstroms
+real	rv			# A(V) / E(B-V)
+
+real	x, y, a, b
+
+begin
+	# Convert to inverse microns
+	x = 10000. / wavelength
+	x = max (0.3, min (10., x))
+
+	# Compute a(x) and b(x)
+	if (x < 0.3) {
+	    call error (1, "Wavelength out of range of extinction function")
+
+	} else if (x < 1.1) {
+	    y = x ** 1.61
+	    a = 0.574 * y
+	    b = -0.527 * y
+
+	} else if (x < 3.3) {
+	    y = x - 1.82
+	    a = 1 + y * (0.17699 + y * (-0.50447 + y * (-0.02427 +
+		y * (0.72085 + y * (0.01979 + y * (-0.77530 + y * 0.32999))))))
+	    b = y * (1.41338 + y * (2.28305 + y * (1.07233 + y * (-5.38434 +
+		y * (-0.62251 + y * (5.30260 + y * (-2.09002)))))))
+
+	} else if (x < 5.9) {
+	    y = (x - 4.67) ** 2
+	    a = 1.752 - 0.316 * x - 0.104 / (y + 0.341)
+	    b = -3.090 + 1.825 * x + 1.206 / (y + 0.263)
+
+	} else if (x < 8.0) {
+	    y = (x - 4.67) ** 2
+	    a = 1.752 - 0.316 * x - 0.104 / (y + 0.341)
+	    b = -3.090 + 1.825 * x + 1.206 / (y + 0.263)
+
+	    y = x - 5.9
+	    a = a - 0.04473 * y**2 - 0.009779 * y**3
+	    b = b + 0.2130 * y**2 + 0.1207 * y**3
+
+	} else if (x <= 10.0) {
+	    y = x - 8
+	    a = -1.072 - 0.628 * y + 0.137 * y**2 - 0.070 * y**3
+	    b = 13.670 + 4.257 * y - 0.420 * y**2 + 0.374 * y**3
+
+	} else {
+	    call error (1, "Wavelength out of range of extinction function")
+	    
+	}
+
+	# Compute A(lambda)/A(V)
+	y = a + b / rv
+	return (y)
+end
+
+
+# STGETR -- Get real parameter.
+# Returns INDEFR if parameter is not defined.
+
+real procedure stgetr (st, param, table, default)
+
+pointer st		#I SPECTIME structure
+char	param[ARB]	#I Parameter name
+char	table[ARB]	#I Name of table
+real	default		#I Default value
+real	val		#O Returned value
+
+real	clgetr(), tabgetr()
+errchk	tabgetr
+
+begin
+	# Get parameter from CL.
+	val = clgetr (param)
+
+	# If the value is INDEF get the value from the table.
+	if (IS_INDEFR(val)) {
+	    iferr (val = tabgetr (ST_TAB(st), table, "spectrograph", param))
+		val = INDEFR
+	}
+
+	# If the value is INDEF set default.
+	if (IS_INDEFR(val))
+	    val = default
+
+	return (val)
+end
+
+
+# STGETI -- Get integer parameter.
+# Returns INDEFI if parameter is not defined.
+
+int procedure stgeti (st, param, table, default)
+
+pointer st		#I SPECTIME structure
+char	param[ARB]	#I Parameter name
+char	table[ARB]	#I Name of table
+int	default		#I Default value
+int	val		#O Returned value
+
+int	clgeti(), tabgeti()
+errchk	tabgeti
+
+begin
+	# Get parameter from CL.
+	val = clgeti (param)
+
+	# If the value is INDEF get the value from the table.
+	if (IS_INDEFI(val)) {
+	    iferr (val = tabgeti (ST_TAB(st), table, "spectrograph", param))
+		val = INDEFI
+	}
+
+	# If the value is INDEF set the default.
+	if (IS_INDEFI(val))
+	    val = default
+
+	return (val)
+end
+
+
+# STGSTR -- Get string parameter.
+# Returns null string if parameter is not defined.
+
+procedure stgstr (st, param, table, default, val, maxchar)
+
+pointer st		#I SPECTIME structure
+char	param[ARB]	#I Parameter name
+char	table[ARB]	#I Name of table
+char	default[ARB]	#I Default value
+char	val[ARB]	#O Returned value
+int	maxchar		#I Maximum string length
+
+char	tmp[10]
+int	nowhite()
+errchk	tabgste
+
+begin
+	# Get parameter from CL.
+	call clgstr (param, val, maxchar)
+
+	# If the value is a null string get the value from a table.
+	if (nowhite (val, tmp, 10) == 0) {
+	    iferr (call tabgstr (ST_TAB(st), table, "spectrograph", param,
+		val, maxchar))
+		val[1] = EOS
+	}
+
+	# If the value is null set the default value.
+	if (val[1] == EOS)
+	    call strcpy (default, val, maxchar)
+end