aboutsummaryrefslogtreecommitdiff
path: root/noao/obsutil/src/sptime/t_sptime.x
diff options
context:
space:
mode:
Diffstat (limited to 'noao/obsutil/src/sptime/t_sptime.x')
-rw-r--r--noao/obsutil/src/sptime/t_sptime.x2530
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
generated by cgit (git 2.34.1) at 2025-09-20 10:44:08 -0400