path: root/noao/obsutil/src/sptime/t_sptime.x

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