aboutsummaryrefslogtreecommitdiff
path: root/noao/obsutil/src/sptime
diff options
context:
space:
mode:
Diffstat (limited to 'noao/obsutil/src/sptime')
-rw-r--r--noao/obsutil/src/sptime/Revisions81
-rw-r--r--noao/obsutil/src/sptime/abzero.cl10
-rw-r--r--noao/obsutil/src/sptime/blazeang.cl24
-rw-r--r--noao/obsutil/src/sptime/blazefunc.cl76
-rw-r--r--noao/obsutil/src/sptime/grating.x1107
-rw-r--r--noao/obsutil/src/sptime/lib/abjohnson17
-rw-r--r--noao/obsutil/src/sptime/lib/circle21
-rw-r--r--noao/obsutil/src/sptime/lib/slit103
-rw-r--r--noao/obsutil/src/sptime/mkcircle.cl16
-rw-r--r--noao/obsutil/src/sptime/mkpkg20
-rw-r--r--noao/obsutil/src/sptime/mkslit.cl37
-rw-r--r--noao/obsutil/src/sptime/rates.cl74
-rw-r--r--noao/obsutil/src/sptime/specpars.par85
-rw-r--r--noao/obsutil/src/sptime/sptime.h132
-rw-r--r--noao/obsutil/src/sptime/sptime.par53
-rw-r--r--noao/obsutil/src/sptime/stdisperser.x455
-rw-r--r--noao/obsutil/src/sptime/t_cgiparse.x110
-rw-r--r--noao/obsutil/src/sptime/t_sptime.x2530
-rw-r--r--noao/obsutil/src/sptime/tabinterp.x698
-rw-r--r--noao/obsutil/src/sptime/x_spectime.x2
20 files changed, 5651 insertions, 0 deletions
diff --git a/noao/obsutil/src/sptime/Revisions b/noao/obsutil/src/sptime/Revisions
new file mode 100644
index 00000000..67cc787b
--- /dev/null
+++ b/noao/obsutil/src/sptime/Revisions
@@ -0,0 +1,81 @@
+.help revisions Jun01 spectime
+.nf
+
+t_sptime.x
+ 1. The use of the fiber diameter to set the aperture was not working.
+ 2. The listing of the fiber information was not working because it was
+ looking for a table "fibers" instead of "fiber".
+ 3. Fixed typo in label "indivdual" -> "individual".
+ (8/13/04, Valdes)
+
+============================
+SPECTIME V2.2: May 15, 2003
+============================
+
+t_sptime.x
+sptime.par
+sptime.h
+ 1. Improved algorithm for handling saturation.
+ 2. Added a minimum exposure parameter.
+ (5/15/03, Valdes)
+
+t_sptime.x
+ The thermal background calculation was wrong. (3/17/03, Valdes)
+
+============================
+SPECTIME V2.1: March 3, 2003
+============================
+
+sptime.h
+t_sptime.x
+sptime.par
+../doc/sptime.hlp
+ Added a "shuffle" option for the sky subtraction. This assumes half
+ the cycle is object and half is sky and the same number of pixels are
+ used for both. The times in the calculator are the total time
+ (object+sky). (2/10/03)
+
+specpars.par
+ The choices for aperture type listed "cicle" instead of the correct
+ "circular". (2/10/03, Valdes)
+
+t_cgiparse.x +
+mkpkg
+x_spectime.x
+../../obsutil.cl
+ Task to parse the CGI QUERY_STRING and set task parameters.
+ (4/19/02, Valdes)
+
+============================
+SPECTIME V2.0: August 15, 2001
+============================
+
+Greatly revised version.
+
+t_sptime.x
+ Fixed a type mismatch in a max() function. (6/13/02, Valdes)
+
+============================
+SPECTIME V1.2: June 13, 2001
+============================
+
+mkpkg
+ Added missing dependencies for t_sptime.x and tabinterp.x (12/13/01, MJF)
+
+t_sptime.x
+sptime.par
+sptime.h
+mkpkg
+ 1. Added "units" parameter to allow different dispersion units.
+ This requires linking with libsmw.a from the noao package.
+ 2. Added IR bands.
+ 3. Changed so that transmisions of 1 are not reported.
+ 4. Added a parameter to the filter to override order overlap message.
+ The idea is that a single filter function can be used without requiring
+ the user to chose the filter.
+ (6/13/01, Valdes)
+
+================================
+SPECTIME V1.1: February 11, 2000
+================================
+.endhelp
diff --git a/noao/obsutil/src/sptime/abzero.cl b/noao/obsutil/src/sptime/abzero.cl
new file mode 100644
index 00000000..c8b3e76b
--- /dev/null
+++ b/noao/obsutil/src/sptime/abzero.cl
@@ -0,0 +1,10 @@
+procedure abzero (w, logf)
+
+real w {prompt="Wavelength (microns)"}
+real logf {prompt="Log f_lambda (W/cm^2/micron)"}
+
+begin
+ x = 10000 * w
+ y = -2.5 * (logf + 3) - 5 * log10 (x) - 2.4
+ printf ("%d %.3f\n", x, y)
+end
diff --git a/noao/obsutil/src/sptime/blazeang.cl b/noao/obsutil/src/sptime/blazeang.cl
new file mode 100644
index 00000000..2c55b3a2
--- /dev/null
+++ b/noao/obsutil/src/sptime/blazeang.cl
@@ -0,0 +1,24 @@
+procedure blazeang (g, w)
+
+real g = 316 {prompt="l/mm"}
+real w = 7500 {prompt="Blaze wavelength (A)"}
+real phi = 46. {prompt="Camera-collimator angle (deg)"}
+real m = 1 {prompt="Order"}
+real n = 1. {prompt="Index of refraction"}
+real prism = 22 {prompt="Prism angle (deg)"}
+
+begin
+ real dtor, val
+
+ dtor = 3.14159 / 180.
+
+ if (n <= 1.) {
+ val = g * w * m / cos (phi/2*dtor) / 2e7
+ val = atan2 (val, sqrt (1 - val**2)) / dtor
+ } else {
+# val = g * w * m / 1e7 / (n - 1.)
+# val = atan2 (val, sqrt (1 - val**2)) / dtor
+ val = g * w * m / 1e7 / sin (dtor * prism) + 1
+ }
+ printf ("%.4g\n", val)
+end
diff --git a/noao/obsutil/src/sptime/blazefunc.cl b/noao/obsutil/src/sptime/blazefunc.cl
new file mode 100644
index 00000000..31416fc7
--- /dev/null
+++ b/noao/obsutil/src/sptime/blazefunc.cl
@@ -0,0 +1,76 @@
+procedure blazefunc (grating, order)
+
+file grating {prompt="Grating"}
+int order = INDEF {prompt="Order"}
+real camcolangle = 45. {prompt="Camera-grating-collimator angle (deg)"}
+string search = "spectimedb$KPNO/Gratings" {prompt="Directory search list\n"}
+
+string title = "" {prompt="Title"}
+real w1 = 3000. {prompt="Lower wavelength to plot"}
+real w2 = 12000. {prompt="Upper wavelength to plot\n"}
+
+real x1 = 3000. {prompt="Left graph wavelength"}
+real x2 = 12000. {prompt="Right graph wavelength"}
+real y1 = -5. {prompt="Bottom graph efficiency"}
+real y2 = 105. {prompt="Top graph efficiency"}
+string ltype = "1" {prompt="Line type"}
+string color = "1" {prompt="Color"}
+bool append = no {prompt="Append?"}
+
+struct *fd
+
+begin
+ file tmp1, tmp2
+ real x, y
+
+ tmp1 = mktemp ("tmp$iraf")
+ tmp2 = mktemp ("tmp$iraf")
+
+ # Spectrograph.
+ print ("# area = 1", >> tmp1)
+ print ("# scale = 1", >> tmp1)
+ print ("# fl = 1", >> tmp1)
+ print ("# ndisp = 2000", >> tmp1)
+ print ("# pixsize = 1", >> tmp1)
+
+ sptime (time=1., maxexp=3600., sn=25., spectrum="blackbody",
+ sky="", sensfunc="none", airmass=1., seeing=1., phase=0.,
+ temperature=6000., index=0., refwave=INDEF, refflux=10.,
+ funits="AB", abjohnson="none", wave=INDEF, order=order,
+ xorder=INDEF, width=1., length=1., diameter=1.,
+ inoutangle=camcolangle, xinoutangle=INDEF, xbin=1, ybin=1,
+ search=search, spectrograph=tmp1, filter="none",
+ filter2="none", disperser=grating, xdisperser="none",
+ fiber="none", telescope=tmp1, adc="none", collimator=tmp1,
+ corrector="none", camera=tmp1, detector=tmp1, aperture=tmp1,
+ extinction="", gain=1., rdnoise=0., dark=0., skysub="none",
+ nskyaps=10, output="disperser", list=tmp2, graphics="",
+ interactive=no, nw=1000, > "dev$null")
+
+ delete (tmp1, verify-)
+
+ fd = tmp2
+ while (fscan (fd, x, y) != EOF) {
+ if (nscan() != 2)
+ next
+ if (x < w1 || x > w2)
+ next
+ print (x, y, >> tmp1)
+ }
+ fd = ""
+
+ if (title == "")
+ title = grating
+
+ graph (tmp1, wx1=x1, wx2=x2, wy1=y1, wy2=y2, wcs="logical", axis=1,
+ transpose=no, pointmode=no, marker="box", szmarker=0.005,
+ ltypes=ltype, colors=color, logx=no, logy=no, box=yes,
+ ticklabels=yes, xlabel="Wavelength (A)", ylabel="Efficiency (%)",
+ xformat="wcsformat", yformat="", title=title,
+ lintran=no, p1=0., p2=0., q1=0., q2=1., vx1=0., vx2=0., vy1=0.,
+ vy2=0., majrx=5, minrx=5, majry=7, minry=5, overplot=no,
+ append=append, device="stdgraph", round=no, fill=yes)
+
+ delete (tmp1, verify-)
+ delete (tmp2, verify-)
+end
diff --git a/noao/obsutil/src/sptime/grating.x b/noao/obsutil/src/sptime/grating.x
new file mode 100644
index 00000000..373b335d
--- /dev/null
+++ b/noao/obsutil/src/sptime/grating.x
@@ -0,0 +1,1107 @@
+include <error.h>
+include <math.h>
+
+define GR_LEN 24
+define GR_W Memr[P2R($1)] # Wavelength
+define GR_O Memi[$1+1] # Order
+define GR_P Memr[P2R($1+2)] # Blaze peak scale factor
+define GR_WB Memr[P2R($1+3)] # First order wavelength at blaze (A)
+define GR_DB Memr[P2R($1+4)] # First order dispersion at blaze (A/mm)
+define GR_OREF Memi[$1+5] # Reference order
+define GR_F Memr[P2R($1+6)] # Focal length (mm)
+define GR_G Memr[P2R($1+7)] # Ruling (groves/A)
+define GR_BLAZE Memr[P2R($1+8)] # Blaze angle (rad)
+define GR_N Memr[P2R($1+9)] # Index of refraction
+define GR_PHI Memr[P2R($1+10)] # Alpha - beta (rad)
+define GR_ALPHA Memr[P2R($1+11)] # Incident angle (rad)
+define GR_BETA Memr[P2R($1+12)] # Diffraction angle (rad)
+define GR_TYPE Memr[P2R($1+13)] # 1=Reflection, -1=Transmission
+define GR_FULL Memi[$1+14] # Full solution?
+
+define GR_PIS Memr[P2R($1+15)] # PI/G*S
+define GR_CA Memr[P2R($1+16)] # cos (ALPHA)
+define GR_SA Memr[P2R($1+17)] # sin (ALPHA)
+define GR_CB Memr[P2R($1+18)] # cos (BETA)
+define GR_TB Memr[P2R($1+19)] # tan (BETA)
+define GR_SE Memr[P2R($1+20)] # sin (ALPHA - BLAZE)
+define GR_CE Memr[P2R($1+21)] # cos (ALPHA - BLAZE)
+define GR_CBLZ Memr[P2R($1+22)] # cos (BLAZE)
+define GR_T2BLZ Memr[P2R($1+23)] # tan (2 * BLAZE)
+
+
+# Definitions of INDEF parameter flags.
+define F 1B # Focal length
+define G 2B # Groves
+define T 4B # Blaze angle
+define A 10B # Incident angle
+define B 20B # Diffracted angle
+define P 40B # Incident - diffracted
+define W 100B # Wavelength
+define D 200B # Dispersions
+define N 400B # Index of refraction
+
+# Combinations
+define FG 3B
+define FT 5B
+define FA 11B
+define FW 101B
+define GT 6B
+define GA 12B
+define GW 102B
+define GD 202B
+define TA 14B
+define TAW 114B
+define TAD 214B
+define TB 24B
+define TP 44B
+define TW 104B
+define TD 204B
+define AB 30B
+define AP 50B
+define AW 110B
+define AD 210B
+define BP 140B
+define WD 300B
+define ABP 70B
+define GTA 16B
+
+
+# GRATING_OPEN -- Open grating structure.
+# Check and derive grating parameters.
+#
+# This procedure hasn't yet been fixed up for grisms (index of refraction
+# is not accounted for) so all input parameters should be defined with
+# alpha=beta=blaze=prism angle.
+
+pointer procedure gr_open (w, o, p, wb, db, oref, f, gmm, blaze, n, phi,
+ alpha, beta, mode, full)
+
+real w #I Wavelength (A)
+int o #I Order
+real p #I Blaze peak scale factor
+real f #I Focal length (mm)
+real wb #I Blaze wavelength (A)
+real db #I Blaze dispersion (A/mm)
+int oref #I Reference order
+real gmm #I Groves (groves/mm)
+real blaze #I Blaze angle (deg)
+real n #I Index of refraction for grism
+real phi #I Incident - diffracted (deg)
+real alpha #I Incident angle (deg)
+real beta #I Diffracted angle (deg)
+int mode #I 1 = incident > diffracted
+int full #I Do full solution?
+pointer gr #O Grating pointer
+
+int flags
+real x
+define err_ 10
+
+begin
+ call malloc (gr, GR_LEN, TY_STRUCT)
+
+ GR_W(gr) = w
+ GR_O(gr) = o
+ GR_P(gr) = p
+
+ GR_F(gr) = f
+ GR_G(gr) = gmm
+ GR_BLAZE(gr) = blaze
+ GR_N(gr) = n
+ GR_PHI(gr) = phi
+ GR_ALPHA(gr) = alpha
+ GR_BETA(gr) = beta
+
+ GR_OREF(gr) = oref
+ GR_WB(gr) = wb
+ GR_DB(gr) = db
+
+ # The grating is reflection unless the index of refraction is not
+ # 1, the blaze dispersion is negative, beta is greater than
+ # 180 degrees, or the incident and diffraction angles are the same.
+
+ if (GR_N(gr) == 1.)
+ GR_TYPE(gr) = 1
+ else if (!IS_INDEF(GR_N(gr)))
+ GR_TYPE(gr) = -1
+ else {
+ if (!IS_INDEF(GR_BETA(gr))
+ && (GR_BETA(gr) > 180. || GR_BETA(gr) < -180))
+ GR_TYPE(gr) = -1
+ else if (!IS_INDEF(GR_DB(gr)) && GR_DB(gr) < 0.)
+ GR_TYPE(gr) = -1
+ else if (!IS_INDEF(GR_ALPHA(gr)) && !IS_INDEF(GR_BETA(gr)) &&
+ GR_ALPHA(gr) == GR_BETA(gr))
+ GR_TYPE(gr) = -1
+ else if (GR_PHI(gr) == 0.)
+ GR_TYPE(gr) = -1
+ else
+ GR_TYPE(gr) = 1
+ }
+
+ # Set INDEF values to reasonable defaults. Convert degrees to radians.
+ if (IS_INDEF(GR_P(gr)))
+ GR_P(gr) = 1
+ if (!IS_INDEF(GR_WB(gr))) {
+ if (GR_WB(gr) <= 0.)
+ GR_WB(gr) = INDEF
+ else
+ GR_WB(gr) = GR_WB(gr) * GR_OREF(gr)
+ }
+ if (!IS_INDEF(GR_DB(gr)))
+ GR_DB(gr) = GR_TYPE(gr) * GR_OREF(gr) * abs (GR_DB(gr))
+ if (!IS_INDEF(GR_F(gr))) {
+ if (GR_F(gr) <= 0.)
+ GR_F(gr) = INDEF
+ }
+ if (!IS_INDEF(GR_G(gr))) {
+ if (GR_G(gr) <= 0.)
+ GR_G(gr) = INDEF
+ else
+ GR_G(gr) = GR_G(gr) / 1e7
+ }
+ if (!IS_INDEF(GR_PHI(gr)))
+ GR_PHI(gr) = DEGTORAD (GR_PHI(gr))
+ if (!IS_INDEF(GR_ALPHA(gr)))
+ GR_ALPHA(gr) = DEGTORAD (GR_ALPHA(gr))
+ if (!IS_INDEF(GR_BETA(gr))) {
+ GR_BETA(gr) = DEGTORAD (GR_BETA(gr))
+ if (GR_BETA(gr) > HALFPI)
+ GR_BETA(gr) = GR_BETA(gr) - PI
+ else if (GR_BETA(gr) < -HALFPI)
+ GR_BETA(gr) = GR_BETA(gr) + PI
+ }
+ if (!IS_INDEF(GR_BLAZE(gr)))
+ GR_BLAZE(gr) = DEGTORAD (GR_BLAZE(gr))
+
+ # Compute missing angles, if possible, based on the other angles.
+ # This assumes the given values are for the blaze peak.
+
+ flags = 0
+ if (IS_INDEF(GR_BLAZE(gr)))
+ flags = flags + T
+ if (IS_INDEF(GR_ALPHA(gr)))
+ flags = flags + A
+ if (IS_INDEF(GR_BETA(gr)))
+ flags = flags + B
+ if (IS_INDEF(GR_PHI(gr)))
+ flags = flags + P
+
+ switch (flags) {
+ case T, P, TP:
+ GR_PHI(gr) = GR_ALPHA(gr) - GR_BETA(gr)
+ GR_BLAZE(gr) = (GR_ALPHA(gr) + GR_BETA(gr)) / 2.
+ case A, TA:
+ GR_ALPHA(gr) = GR_BETA(gr) + GR_PHI(gr)
+ GR_BLAZE(gr) = (GR_ALPHA(gr) + GR_BETA(gr)) / 2.
+ case AB:
+ if (mode == 1) {
+ GR_ALPHA(gr) = GR_BLAZE(gr) + GR_PHI(gr)/2.
+ GR_BETA(gr) = GR_BLAZE(gr) - GR_PHI(gr)/2.
+ } else {
+ GR_ALPHA(gr) = GR_BLAZE(gr) - GR_PHI(gr)/2.
+ GR_BETA(gr) = GR_BLAZE(gr) + GR_PHI(gr)/2.
+ }
+ case AP:
+ GR_ALPHA(gr) = 2 * GR_BLAZE(gr) - GR_BETA(gr)
+ GR_PHI(gr) = GR_ALPHA(gr) - GR_BETA(gr)
+ case B, TB:
+ GR_BETA(gr) = GR_ALPHA(gr) - GR_PHI(gr)
+ GR_BLAZE(gr) = (GR_ALPHA(gr) + GR_BETA(gr)) / 2.
+ case BP:
+ GR_BETA(gr) = 2 * GR_BLAZE(gr) - GR_ALPHA(gr)
+ GR_PHI(gr) = GR_ALPHA(gr) - GR_BETA(gr)
+ case ABP:
+ GR_ALPHA(gr) = GR_BLAZE(gr)
+ GR_BETA(gr) = GR_BLAZE(gr)
+ GR_PHI(gr) = 0.
+ }
+
+ # Compute index of refraction if possible.
+ if (IS_INDEF(GR_N(gr))) {
+ if (GR_TYPE(gr) == 1)
+ GR_N(gr) = 1
+ else if (!IS_INDEF(GR_WB(gr)) && !IS_INDEF(GR_G(gr)) &&
+ !IS_INDEF(GR_BLAZE(gr)))
+ GR_N(gr) = (GR_G(gr) * GR_WB(gr)) / sin (GR_BLAZE(gr)) + 1
+ }
+
+ # Compute other parameters if possible.
+ flags = 0
+ if (IS_INDEF(GR_F(gr)))
+ flags = flags + F
+ if (IS_INDEF(GR_G(gr)))
+ flags = flags + G
+ if (IS_INDEF(GR_BLAZE(gr)))
+ flags = flags + T
+ if (IS_INDEF(GR_ALPHA(gr)))
+ flags = flags + A
+ if (IS_INDEF(GR_WB(gr)))
+ flags = flags + W
+ if (IS_INDEF(GR_DB(gr)))
+ flags = flags + D
+ if (IS_INDEF(GR_N(gr)))
+ flags = flags + N
+
+ switch (flags) {
+ case 0, F, G, T, A, N, W, D:
+ switch (flags) {
+ case F:
+ GR_F(gr) = GR_TYPE(gr) * cos (GR_BETA(gr)) /
+ (GR_G(gr) * GR_DB(gr))
+ case G:
+ GR_G(gr) = (GR_N(gr) * sin (GR_ALPHA(gr)) +
+ GR_TYPE(gr) * sin (GR_BETA(gr))) / GR_WB(gr)
+ if (GR_G(gr) == 0.)
+ GR_G(gr) = INDEF
+ case T:
+ if (GR_ALPHA(gr) > PI) {
+ x = GR_G(gr) * GR_WB(gr) / (2 * cos (GR_ALPHA(gr)))
+ if (abs (x) > 1.)
+ goto err_
+ GR_BLAZE(gr) = asin (x)
+ GR_ALPHA(gr) = GR_ALPHA(gr) - TWOPI + GR_BLAZE(gr)
+ } else {
+ x = GR_G(gr) * GR_WB(gr) - GR_N(gr) * sin (GR_ALPHA(gr))
+ if (abs (x) > 1.)
+ goto err_
+ GR_BLAZE(gr) = (GR_ALPHA(gr) + asin (x)) / 2
+ }
+ GR_BETA(gr) = 2 * GR_BLAZE(gr) - GR_ALPHA(gr)
+ GR_PHI(gr) = GR_ALPHA(gr) - GR_BETA(gr)
+ case A:
+ x = GR_TYPE(gr) * GR_G(gr) * GR_WB(gr) / (2 * sin(GR_BLAZE(gr)))
+ if (abs (x) > 1.)
+ goto err_
+ if (mode == 1) {
+ GR_ALPHA(gr) = GR_BLAZE(gr) + acos (x)
+ GR_BETA(gr) = 2 * GR_BLAZE(gr) - GR_ALPHA(gr)
+ } else {
+ GR_BETA(gr) = GR_BLAZE(gr) + acos (x)
+ GR_ALPHA(gr) = 2 * GR_BLAZE(gr) - GR_BETA(gr)
+ }
+ GR_PHI(gr) = GR_ALPHA(gr) - GR_BETA(gr)
+ case N:
+ GR_N(gr) = (GR_G(gr) * GR_WB(gr)) / sin (GR_BLAZE(gr)) + 1
+ }
+ GR_WB(gr) = (GR_N(gr) * sin (GR_ALPHA(gr)) +
+ GR_TYPE(gr) * sin (GR_BETA(gr))) / GR_G(gr)
+ GR_DB(gr) = GR_TYPE(gr) * cos (GR_BETA(gr)) / (GR_F(gr)*GR_G(gr))
+ case FG:
+ x = (GR_N(gr) * sin (GR_ALPHA(gr)) +
+ GR_TYPE(gr) * sin (GR_BETA(gr))) / GR_WB(gr)
+ if (x == 0.)
+ goto err_
+ GR_G(gr) = x
+ GR_F(gr) = GR_TYPE(gr) * cos (GR_BETA(gr)) / (GR_G(gr) * GR_DB(gr))
+ case FT:
+ if (GR_ALPHA(gr) > PI) {
+ x = GR_TYPE(gr) * GR_G(gr) * GR_WB(gr) /
+ (2 * cos (GR_ALPHA(gr)))
+ if (abs (x) > 1.)
+ goto err_
+ GR_BLAZE(gr) = asin (x)
+ GR_ALPHA(gr) = GR_ALPHA(gr) - TWOPI + GR_BLAZE(gr)
+ } else {
+ x = GR_TYPE(gr) * GR_G(gr) * GR_WB(gr) -
+ GR_N(gr) * sin (GR_ALPHA(gr))
+ if (abs (x) > 1.)
+ goto err_
+ GR_BLAZE(gr) = (GR_ALPHA(gr) + asin (x)) / 2
+ }
+ GR_BETA(gr) = 2 * GR_BLAZE(gr) - GR_ALPHA(gr)
+ GR_PHI(gr) = GR_ALPHA(gr) - GR_BETA(gr)
+ GR_F(gr) = GR_TYPE(gr) * cos (GR_BETA(gr)) / (GR_G(gr) * GR_DB(gr))
+ case FA:
+ x = GR_TYPE(gr) * GR_G(gr) * GR_WB(gr) / (2 * sin (GR_BLAZE(gr)))
+ if (abs (x) > 1.)
+ goto err_
+ if (mode == 1) {
+ GR_ALPHA(gr) = GR_BLAZE(gr) + acos (x)
+ GR_BETA(gr) = 2 * GR_BLAZE(gr) - GR_ALPHA(gr)
+ } else {
+ GR_BETA(gr) = GR_BLAZE(gr) + acos (x)
+ GR_ALPHA(gr) = 2 * GR_BLAZE(gr) - GR_BETA(gr)
+ }
+ GR_PHI(gr) = GR_ALPHA(gr) - GR_BETA(gr)
+ GR_F(gr) = GR_TYPE(gr) * cos (GR_BETA(gr)) / (GR_G(gr) * GR_DB(gr))
+ case FW:
+ GR_WB(gr) = (GR_N(gr) * sin (GR_ALPHA(gr)) +
+ GR_TYPE(gr) * sin (GR_BETA(gr))) / GR_G(gr)
+ GR_F(gr) = GR_TYPE(gr) * cos (GR_BETA(gr)) / (GR_G(gr) * GR_DB(gr))
+ case GT:
+ x = GR_TYPE(gr) * GR_F(gr) * GR_DB(gr) / GR_WB(gr)
+ if (GR_ALPHA(gr) > PI) {
+ GR_BLAZE(gr) = atan (1 / (2 * x - tan (GR_ALPHA(gr))))
+ GR_ALPHA(gr) = GR_ALPHA(gr) - TWOPI + GR_BLAZE(gr)
+ } else {
+ x = (tan (GR_ALPHA(gr)) - x) / (1 + 2 * x * tan (GR_ALPHA(gr)))
+ GR_BLAZE(gr) = atan (x + sqrt (1 + x * x))
+ }
+ GR_BETA(gr) = 2 * GR_BLAZE(gr) - GR_ALPHA(gr)
+ GR_PHI(gr) = GR_ALPHA(gr) - GR_BETA(gr)
+ GR_G(gr) = (GR_N(gr) * sin (GR_ALPHA(gr)) +
+ GR_TYPE(gr) * sin (GR_BETA(gr))) / GR_WB(gr)
+ case GA:
+ GR_ALPHA(gr) = GR_BLAZE(gr) +
+ atan (2 * GR_TYPE(gr) * GR_F(gr) * GR_DB(gr) /
+ GR_WB(gr) - 1 / tan (GR_BLAZE(gr)))
+ GR_BETA(gr) = 2 * GR_BLAZE(gr) - GR_ALPHA(gr)
+ GR_PHI(gr) = GR_ALPHA(gr) - GR_BETA(gr)
+ GR_G(gr) = (GR_N(gr) * sin (GR_ALPHA(gr)) +
+ GR_TYPE(gr) * sin (GR_BETA(gr))) / GR_WB(gr)
+ case GW:
+ GR_G(gr) = GR_TYPE(gr) * cos (GR_BETA(gr)) / (GR_F(gr) * GR_DB(gr))
+ if (GR_G(gr) == 0.)
+ GR_G(gr) = INDEF
+ else
+ GR_WB(gr) = (GR_N(gr) * sin (GR_ALPHA(gr)) +
+ GR_TYPE(gr) * sin (GR_BETA(gr))) / GR_G(gr)
+ case GD:
+ x = (GR_N(gr) * sin (GR_ALPHA(gr)) +
+ GR_TYPE(gr) * sin (GR_BETA(gr))) / GR_WB(gr)
+ if (x == 0.)
+ goto err_
+ GR_G(gr) = x
+ GR_DB(gr) = GR_TYPE(gr) * cos (GR_BETA(gr)) / (GR_F(gr) * GR_G(gr))
+ case TAD:
+ if (IS_INDEF(GR_PHI(gr)))
+ GR_PHI(gr) = 0.
+ x = GR_G(gr) * GR_WB(gr) / (2 * cos (GR_PHI(gr)/2))
+ if (mode == 1) {
+ GR_ALPHA(gr) = GR_PHI(gr)/2 + asin (x)
+ GR_BETA(gr) = GR_ALPHA(gr) - GR_PHI(gr)
+ } else {
+ GR_BETA(gr) = GR_PHI(gr)/2 + asin (x)
+ GR_ALPHA(gr) = GR_BETA(gr) - GR_PHI(gr)
+ }
+ GR_BLAZE(gr) = (GR_ALPHA(gr) + GR_BETA(gr)) / 2
+ GR_DB(gr) = GR_TYPE(gr) * cos (GR_BETA(gr)) /
+ (GR_F(gr) * GR_G(gr))
+ case TAW:
+ if (!IS_INDEF(GR_PHI(gr))) {
+ x = GR_TYPE(gr) * GR_F(gr) * GR_G(gr) * GR_DB(gr)
+ if (abs (x) > 1.) {
+ if (abs (x) > 1.1)
+ goto err_
+ else {
+ x = 1.
+ GR_DB(gr) = x / (GR_F(gr) * GR_G(gr))
+ }
+ }
+ if (mode == 1) {
+ GR_BETA(gr) = acos (x)
+ GR_ALPHA(gr) = GR_BETA(gr) + GR_PHI(gr)
+ } else {
+ GR_ALPHA(gr) = acos (x)
+ GR_BETA(gr) = GR_ALPHA(gr) + GR_PHI(gr)
+ }
+ GR_BLAZE(gr) = (GR_ALPHA(gr) + GR_BETA(gr)) / 2
+ GR_WB(gr) = (GR_N(gr) * sin (GR_ALPHA(gr)) +
+ GR_TYPE(gr) * sin (GR_BETA(gr))) / GR_G(gr)
+ }
+ case TA:
+ if (!IS_INDEF(GR_PHI(gr))) {
+ x = GR_G(gr) * GR_WB(gr) / (2 * cos (GR_PHI(gr)/2))
+ if (mode == 1) {
+ GR_ALPHA(gr) = GR_PHI(gr)/2 + asin (x)
+ GR_BETA(gr) = GR_ALPHA(gr) - GR_PHI(gr)
+ } else {
+ GR_BETA(gr) = GR_PHI(gr)/2 + asin (x)
+ GR_ALPHA(gr) = GR_BETA(gr) - GR_PHI(gr)
+ }
+ GR_BLAZE(gr) = (GR_ALPHA(gr) + GR_BETA(gr)) / 2
+ GR_DB(gr) = GR_TYPE(gr) * cos (GR_BETA(gr)) /
+ (GR_F(gr) * GR_G(gr))
+ } else {
+ x = GR_TYPE(gr) * GR_F(gr) * GR_G(gr) * GR_DB(gr)
+ if (abs (x) > 1.) {
+ if (abs (x) > 1.1)
+ goto err_
+ else {
+ x = 1.
+ GR_DB(gr) = x / (GR_F(gr) * GR_G(gr))
+ }
+ }
+ GR_BETA(gr) = acos (x)
+ x = GR_G(gr) * GR_WB(gr) - GR_TYPE(gr) * sin (GR_BETA(gr))
+ if (abs (x) > 1.)
+ goto err_
+ GR_ALPHA(gr) = asin (x)
+ GR_BLAZE(gr) = (acos (GR_TYPE(gr) * GR_F(gr) * GR_G(gr) *
+ GR_DB(gr)) + GR_ALPHA(gr)) / 2
+ GR_BETA(gr) = 2 * GR_BLAZE(gr) - GR_ALPHA(gr)
+ GR_PHI(gr) = GR_ALPHA(gr) - GR_BETA(gr)
+ }
+ case TW:
+ GR_BLAZE(gr) = (GR_ALPHA(gr) +
+ acos (GR_TYPE(gr) * GR_F(gr) * GR_G(gr) * GR_DB(gr))) / 2
+ GR_BETA(gr) = 2 * GR_BLAZE(gr) - GR_ALPHA(gr)
+ GR_PHI(gr) = GR_ALPHA(gr) - GR_BETA(gr)
+ GR_WB(gr) = (GR_N(gr) * sin (GR_ALPHA(gr)) +
+ GR_TYPE(gr) * sin (GR_BETA(gr))) / GR_G(gr)
+ case TD:
+ if (GR_ALPHA(gr) > PI) {
+ x = GR_G(gr) * GR_WB(gr) / (2 * cos (GR_ALPHA(gr)))
+ if (abs (x) > 1.)
+ goto err_
+ GR_BLAZE(gr) = asin (x)
+ GR_ALPHA(gr) = GR_ALPHA(gr) - TWOPI + GR_BLAZE(gr)
+ } else {
+ x = GR_G(gr) * GR_WB(gr) - GR_N(gr) * sin (GR_ALPHA(gr))
+ if (abs (x) > 1.)
+ goto err_
+ GR_BLAZE(gr) = (GR_ALPHA(gr) + asin (x)) / 2
+ }
+ GR_BETA(gr) = 2 * GR_BLAZE(gr) - GR_ALPHA(gr)
+ GR_PHI(gr) = GR_ALPHA(gr) - GR_BETA(gr)
+ GR_DB(gr) = GR_TYPE(gr) * cos (GR_BETA(gr)) / (GR_F(gr) * GR_G(gr))
+ case AW:
+ x = GR_TYPE(gr) * GR_F(gr) * GR_G(gr) * GR_DB(gr)
+ if (abs (x) > 1.)
+ goto err_
+ GR_ALPHA(gr) = 2 * GR_BLAZE(gr) - acos (x)
+ GR_BETA(gr) = 2 * GR_BLAZE(gr) - GR_ALPHA(gr)
+ GR_PHI(gr) = GR_ALPHA(gr) - GR_BETA(gr)
+ GR_WB(gr) = (GR_N(gr) * sin (GR_ALPHA(gr)) + sin (GR_BETA(gr))) /
+ GR_G(gr)
+ case AD:
+ x = GR_G(gr) * GR_WB(gr) / (2 * sin (GR_BLAZE(gr)))
+ if (abs (x) > 1.)
+ goto err_
+ if (mode == 1) {
+ GR_ALPHA(gr) = GR_BLAZE(gr) + acos (x)
+ GR_BETA(gr) = 2 * GR_BLAZE(gr) - GR_ALPHA(gr)
+ } else {
+ GR_BETA(gr) = GR_BLAZE(gr) + acos (x)
+ GR_ALPHA(gr) = 2 * GR_BLAZE(gr) - GR_BETA(gr)
+ }
+ GR_PHI(gr) = GR_ALPHA(gr) - GR_BETA(gr)
+ GR_DB(gr) = GR_TYPE(gr) * cos (GR_BETA(gr)) / (GR_F(gr) * GR_G(gr))
+ case WD:
+ GR_WB(gr) = (GR_N(gr) * sin (GR_ALPHA(gr)) +
+ GR_TYPE(gr) * sin (GR_BETA(gr))) / GR_G(gr)
+ GR_DB(gr) = GR_TYPE(gr) * cos (GR_BETA(gr)) / (GR_F(gr) * GR_G(gr))
+ case GTA:
+ # Assume beta=alpha=blaze.
+ x = (GR_TYPE(gr) * GR_WB(gr)) /
+ ((GR_N(gr) + GR_TYPE(gr)) * GR_F(gr) * GR_DB(gr))
+ GR_BETA(gr) = atan (x)
+ GR_ALPHA(gr) = GR_BETA(gr)
+ GR_BLAZE(gr) = GR_BETA(gr)
+ GR_PHI(gr) = 0
+ GR_G(gr) = (GR_N(gr) * sin (GR_ALPHA(gr)) +
+ GR_TYPE(gr) * sin (GR_BETA(gr))) / GR_WB(gr)
+ }
+
+ # The result should give the blaze wavelength and dispersion.
+ # If this cannot be computed then it is an error.
+
+ if (IS_INDEF(GR_WB(gr)) || IS_INDEF(GR_DB(gr))) {
+ call gr_close (gr)
+ call mfree (gr, TY_STRUCT)
+ call error (1,
+ "Insufficient information to to resolve grating parameters")
+ }
+
+ # If all the other parameters cannot be computed then use linear.
+ if (full==NO || IS_INDEF(GR_F(gr)) || IS_INDEF(GR_G(gr)) ||
+ IS_INDEF(GR_BETA(gr)) || IS_INDEF(GR_PHI(gr)) ||
+ IS_INDEF(GR_N(gr))) {
+ GR_FULL(gr) = NO
+ if (IS_INDEF(GR_F(gr)))
+ GR_F(gr) = 1
+ GR_W(gr) = GR_WB(gr)
+ GR_O(gr) = GR_OREF(gr)
+ GR_CE(gr) = PI * GR_DB(gr)
+ GR_CA(gr) = 1.
+ GR_CB(gr) = 1.
+
+ # A full grating solution is possible.
+ } else {
+ GR_FULL(gr) = YES
+
+ # Set the order and wavelength to the blaze values if not given.
+ if (IS_INDEF(GR_W(gr))) {
+ if (!IS_INDEFI(GR_O(gr)))
+ GR_W(gr) = GR_WB(gr) / GR_O(gr)
+ else
+ GR_W(gr) = GR_WB(gr) / GR_OREF(gr)
+ }
+ if (IS_INDEFI(GR_O(gr)))
+ GR_O(gr) = max (1, nint (GR_WB(gr) / GR_W(gr)))
+
+ # Convert to incident angle at desired wavelength.
+ if (GR_PHI(gr) != 0.) {
+ x = GR_G(gr) * GR_O(gr) * GR_W(gr) /
+ (2 * cos (GR_PHI(gr)/2))
+ if (abs (x) > 1.)
+ goto err_
+ if (mode == 1) {
+ GR_ALPHA(gr) = asin (x) + GR_PHI(gr) / 2
+ GR_BETA(gr) = asin (x) - GR_PHI(gr) / 2
+ } else {
+ GR_ALPHA(gr) = asin (x) - GR_PHI(gr) / 2
+ GR_BETA(gr) = asin (x) + GR_PHI(gr) / 2
+ }
+ } else {
+ x = (GR_G(gr) * GR_O(gr) * GR_W(gr) -
+ GR_TYPE(gr) * sin (GR_BETA(gr))) / GR_N(gr)
+ if (abs (x) > 1.)
+ goto err_
+ GR_ALPHA(gr) = asin (x)
+ }
+
+ # The following parameters are for efficiency. Beta terms are
+ # for the blaze peak diffraction.
+
+ x = 2 * GR_BLAZE(gr) - GR_ALPHA(gr)
+ GR_CA(gr) = cos (GR_ALPHA(gr))
+ GR_SA(gr) = sin (GR_ALPHA(gr))
+ GR_CB(gr) = GR_TYPE(gr) * cos (x)
+ GR_TB(gr) = tan (x)
+ GR_SE(gr) = sin (GR_ALPHA(gr) - GR_BLAZE(gr))
+ GR_CE(gr) = cos (GR_ALPHA(gr) - GR_BLAZE(gr))
+ GR_CBLZ(gr) = cos (GR_BLAZE(gr))
+ GR_T2BLZ(gr) = tan (2 * GR_BLAZE(gr))
+ if (GR_ALPHA(gr) > x)
+ GR_PIS(gr) = PI / GR_G(gr) * GR_CA(gr)
+ else
+ GR_PIS(gr) = PI / GR_G(gr) * GR_CBLZ(gr)
+ }
+
+ return (gr)
+
+err_ call error (2, "Impossible combination of grating parameters")
+end
+
+
+# GR_CLOSE -- Free grating structure.
+
+procedure gr_close (gr)
+
+pointer gr #I Grating pointer
+
+begin
+ call mfree (gr, TY_STRUCT)
+end
+
+# GR_GETR -- Get grating parameter.
+
+real procedure gr_getr (gr, param)
+
+pointer gr #I Grating pointer
+char param[ARB] #I Parameter name
+
+bool streq()
+
+begin
+ if (gr == NULL)
+ return (INDEF)
+
+ switch (param[1]) {
+ case 'a':
+ if (streq (param, "alpha")) {
+ if (IS_INDEFR(GR_ALPHA(gr)))
+ return (GR_ALPHA(gr))
+ else
+ return (RADTODEG(GR_ALPHA(gr)))
+ }
+ case 'b':
+ if (streq (param, "blaze")) {
+ if (IS_INDEFR(GR_BLAZE(gr)))
+ return (GR_BLAZE(gr))
+ else
+ return (RADTODEG(GR_BLAZE(gr)))
+ } else if (streq (param, "beta")) {
+ if (IS_INDEFR(GR_BETA(gr)))
+ return (GR_BETA(gr))
+ else
+ return (RADTODEG(GR_BETA(gr)))
+ }
+ case 'd':
+ if (streq (param, "dblaze"))
+ return (GR_DB(gr) / GR_OREF(gr))
+ if (streq (param, "dispersion")) {
+ if (GR_FULL(gr) == NO)
+ return (GR_DB(gr) / GR_O(gr))
+ else
+ return (GR_CB(gr) / (GR_G(gr) * GR_O(gr) *GR_F(gr)))
+ }
+ case 'f':
+ if (streq (param, "full"))
+ return (real (GR_FULL(gr)))
+ else if (streq (param, "f"))
+ return (GR_F(gr))
+ case 'g':
+ if (streq (param, "g")) {
+ if (IS_INDEFR(GR_G(gr)))
+ return (GR_G(gr))
+ else
+ return (GR_G(gr) * 1E7)
+ }
+ case 'm':
+ if (streq (param, "mag"))
+ return (GR_CA(gr) / GR_CB(gr))
+ case 'n':
+ if (streq (param, "n"))
+ return (GR_N(gr))
+ case 'o':
+ if (streq (param, "order"))
+ return (real (GR_O(gr)))
+ if (streq (param, "oref"))
+ return (real (GR_OREF(gr)))
+ case 'p':
+ if (streq (param, "phi")) {
+ if (IS_INDEFR(GR_PHI(gr)))
+ return (GR_PHI(gr))
+ else
+ return (RADTODEG(GR_PHI(gr)))
+ }
+ case 't':
+ if (streq (param, "tilt")) {
+ if (GR_FULL(gr) == NO)
+ return (INDEFR)
+ else
+ return (RADTODEG((GR_ALPHA(gr)+GR_TYPE(gr)*GR_BETA(gr))/2))
+ }
+ case 'w':
+ if (streq (param, "wavelength"))
+ return (GR_W(gr))
+ if (streq (param, "wblaze"))
+ return (GR_WB(gr) / GR_OREF(gr))
+ }
+ call error (1, "gr_getr: unknown parameter")
+end
+
+
+# GR_LIST -- List grating parameters.
+
+procedure gr_list (gr, order, col)
+
+pointer gr #I Grating pointer
+int order #I Order
+int col #I Column to indent
+
+begin
+ if (gr == NULL)
+ return
+
+ if (GR_FULL(gr) == NO) {
+ call printf ("%*tReference order = %d\n")
+ call pargi (col)
+ call pargi (order)
+ call printf (
+ "%*tBlaze wavelength of reference order = %.6g Angstroms\n")
+ call pargi (col)
+ call pargr (GR_WB(gr) / order)
+ call printf (
+ "%*tBlaze dispersion of reference order = %.4g Angstroms/mm\n")
+ call pargi (col)
+ call pargr (GR_DB(gr) / order)
+ } else {
+ call printf ("%*tFocal length = %d mm\n")
+ call pargi (col)
+ call pargr (GR_F(gr))
+ call printf ("%*tGrating = %.1f grooves/mm\n")
+ call pargi (col)
+ call pargr (GR_G(gr) * 1e7)
+ call printf ("%*tBlaze angle = %.1f degrees\n")
+ call pargi (col)
+ call pargr (RADTODEG(GR_BLAZE(gr)))
+ call printf ("%*tIncident to diffracted angle = %.1f degrees\n")
+ call pargi (col)
+ call pargr (RADTODEG(GR_PHI(gr)))
+ call printf ("%*tReference order = %d\n")
+ call pargi (col)
+ call pargi (order)
+ call printf (
+ "%*tBlaze wavelength = %.6g Angstroms\n")
+ call pargi (col)
+ call pargr (GR_WB(gr) / order)
+ call printf (
+ "%*tBlaze dispersion = %.4g Angstroms/mm\n")
+ call pargi (col)
+ call pargr (GR_DB(gr) / order)
+
+ call printf ("\n%*tCentral wavelength = %.6g Angstroms\n")
+ call pargi (col)
+ call pargr (GR_W(gr))
+ call printf ("%*tCentral dispersion = %.4g Angstroms/mm\n")
+ call pargi (col)
+ call pargr (GR_TYPE(gr) * cos (GR_BETA(gr)) /
+ (GR_G(gr) * order *GR_F(gr)))
+ call printf ("%*tOrder = %d\n")
+ call pargi (col)
+ call pargi (GR_O(gr))
+ call printf ("%*tTilt = %.1f degrees\n")
+ call pargi (col)
+ call pargr (RADTODEG(GR_ALPHA(gr) - GR_PHI(gr) / 2))
+ call printf ("%*tGrating magnification = %.2f\n")
+ call pargi (col)
+ call pargr (GR_CA(gr)/GR_CB(gr))
+ call printf ("%*tIncidence angle = %.1f degrees\n")
+ call pargi (col)
+ call pargr (RADTODEG(GR_ALPHA(gr)))
+ call printf ("%*tDiffracted angle = %.1f degrees\n")
+ call pargi (col)
+ call pargr (RADTODEG(GR_BETA(gr)))
+ }
+end
+
+
+# GR_W2PSI -- Given wavelength return tan(psi).
+
+real procedure gr_w2psi (gr, w, m)
+
+pointer gr #I Grating pointer
+real w #I Wavelength
+int m #I Order
+real x #O Pixel position
+
+real gr_sinbeta()
+
+begin
+ if (GR_FULL(gr) == NO)
+ return ((m*w - GR_WB(gr)) / GR_DB(gr) * GR_F(gr))
+
+ x = gr_sinbeta (gr, w, m)
+ if (!IS_INDEF(x)) {
+ x = x / sqrt (1 - x * x)
+ x = (x - GR_TB(gr)) / (1 + x * GR_TB(gr))
+ }
+ return (x)
+
+end
+
+
+# GR_W2PSIR -- Given wavelength return tan(psi) of reflected component.
+
+real procedure gr_w2psir (gr, w, m)
+
+pointer gr #I Grating pointer
+real w #I Wavelength
+int m #I Order
+real x #O Pixel position
+
+real gr_sinbeta()
+
+begin
+ x = gr_sinbeta (gr, w, m)
+ if (!IS_INDEF(x)) {
+ #x = x / sqrt (1 - x * x)
+ #x = (x - GR_TB(gr)) / (1 + x * GR_TB(gr))
+ #x = (x + GR_T2BLZ(gr)) / (1 - x * GR_T2BLZ(gr))
+ x = PI - asin(x)
+ x = 2 * GR_BLAZE(gr) - x - GR_BETA(gr)
+ x = tan (x)
+ }
+ return (x)
+end
+
+
+# GR_DELTA -- Given pixel position and wavelength return blaze function
+# phase angle.
+
+real procedure gr_delta (gr, w, m)
+
+pointer gr #I Grating pointer
+real w #I Wavelength
+int m #I Order
+real d #O Delta
+
+real tanpsi, cospsi, sinpsi, gr_w2psi()
+
+begin
+ tanpsi = gr_w2psi (gr, w, m)
+ if (IS_INDEF(tanpsi))
+ return (INDEF)
+
+ if (GR_FULL(gr) == NO)
+ d = PI * GR_DB(gr) / w * tanpsi
+ else {
+ cospsi = 1 / sqrt (1 + tanpsi * tanpsi)
+ sinpsi = tanpsi * cospsi
+ d = GR_PIS(gr) / w * (GR_CE(gr) * sinpsi + GR_SE(gr) * (1 - cospsi))
+ }
+ if (abs(d) < 0.01) {
+ if (d < 0)
+ d = -0.01
+ else
+ d = 0.01
+ }
+ return (d)
+end
+
+
+# GR_DELTAR -- Blaze function phase angle for reflected component.
+
+real procedure gr_deltar (gr, w, m)
+
+pointer gr #I Grating pointer
+real w #I Wavelength
+int m #I Order
+real d #O Delta
+
+real tanpsi, cospsi, sinpsi, gr_w2psir()
+
+begin
+ tanpsi = gr_w2psir (gr, w, m)
+ if (IS_INDEF(tanpsi))
+ return (INDEF)
+
+ if (GR_FULL(gr) == NO)
+ d = PI * GR_DB(gr) / w * tanpsi
+ else {
+ cospsi = 1 / sqrt (1 + tanpsi * tanpsi)
+ sinpsi = tanpsi * cospsi
+ d = GR_PIS(gr) / w * (GR_CE(gr) * sinpsi + GR_SE(gr) * (1 - cospsi))
+ }
+ if (abs(d) < 0.01) {
+ if (d < 0)
+ d = -0.01
+ else
+ d = 0.01
+ }
+ return (d)
+end
+
+
+# GR_BLAZE -- Blaze pattern at given wavelength.
+
+real procedure gr_blaze (gr, w, m)
+
+pointer gr #I Grating pointer
+real w #I Wavelength
+int m #I Order
+real val #O Blaze pattern
+
+real d, gr_delta()
+
+begin
+ d = gr_delta (gr, w, m)
+ if (IS_INDEF(d))
+ val = 0.
+ else
+ val = (sin (d) / d) ** 2
+ return (val)
+end
+
+
+# GR_PEAK -- Blaze peak corrected for light defracted into other orders.
+
+real procedure gr_peak (gr, w, m)
+
+pointer gr #I Grating pointer
+real w #I Wavelength (A)
+int m #I Order
+real val #O Blaze peak
+
+int i, j
+real frac, p
+real gr_delta(), gr_deltar()
+
+begin
+ # Find the absolute response of the gratings at the reference
+ # blaze peak.
+
+ p = gr_delta (gr, w, m)
+ if (IS_INDEF(p))
+ return (INDEF)
+ val = (sin (p) / p) ** 2
+ frac = 0.
+ do i = m - 1, 1, -1 {
+ p = gr_delta (gr, w, i)
+ if (IS_INDEF(p))
+ break
+ frac = frac + (sin (p) / p) ** 2
+ if (abs (p) > 1000.)
+ break
+ }
+ do i = m + 1, ARB {
+ p = gr_delta (gr, w, i)
+ if (IS_INDEF(p))
+ break
+ frac = frac + (sin (p) / p) ** 2
+ if (abs (p) > 1000.)
+ break
+ }
+
+ if (GR_FULL(gr) == YES && GR_TYPE(gr) > 0) {
+ j = (GR_N(gr) * GR_SA(gr) + GR_TYPE(gr)) / GR_G(gr) / w
+ do i = j+1, ARB, 1 {
+ p = gr_deltar (gr, w, i)
+ if (IS_INDEF(p))
+ break
+ frac = frac + (sin (p) / p) ** 2
+ if (abs (p) > 1000.)
+ break
+ }
+ }
+
+ val = val / (val + frac)
+
+ # Shadowing
+ if (GR_ALPHA(gr) < GR_BETA(gr) && GR_TYPE(gr) > 0)
+ val = val * abs (GR_CA(gr) / GR_CB(gr))
+
+ return (val)
+end
+
+
+# GR_EFF -- Efficiency at specified wavelength and order.
+
+real procedure gr_eff (gr, w, m)
+
+pointer gr #I Grating pointer
+real w #I Wavelength
+int m #I Order
+real eff #O Efficiency
+
+real gr_blaze(), gr_peak()
+
+begin
+ if (gr == NULL)
+ return (INDEF)
+
+ if (GR_FULL(gr) == NO)
+ return (GR_P(gr))
+
+ eff = gr_blaze (gr, w, m)
+ if (eff > 0.)
+ eff = eff * GR_P(gr) * gr_peak (gr, GR_WB(gr)/m, m)
+
+ return (eff)
+end
+
+
+# GR_W2DW -- Grating dispersion = cos (beta(w,m)) / (g * m * f)
+# This is corrected to a detector plane.
+
+real procedure gr_w2dw (gr, w, m)
+
+pointer gr #I Grating pointer
+real w #I Wavelength
+int m #I Order
+real disp #O Dispersion (A/mm)
+
+real gr_sinbeta(), gr_w2x()
+
+begin
+ if (GR_FULL(gr) == NO)
+ return (GR_DB(gr) / m)
+
+ disp = gr_sinbeta (gr, w, m)
+ if (!IS_INDEF(disp)) {
+ disp = sqrt (1 - disp * disp) / (GR_G(gr) * m * GR_F(gr))
+ disp = disp / (1 + (gr_w2x (gr, w, m) / GR_F(gr))**2)
+ }
+ return (disp)
+end
+
+
+# GR_X2W -- Wavelength at given position relative to center.
+
+real procedure gr_x2w (gr, x, m)
+
+pointer gr #I Grating pointer
+real x #I Pixel position (mm from center)
+int m #I Order
+real w #O Wavelength (Angstroms)
+
+begin
+ if (GR_FULL(gr) == NO) {
+ w = GR_WB(gr) + GR_DB(gr) / m * x
+ return (w)
+ }
+
+ w = x / GR_F(gr)
+ w = atan (w) + GR_BETA(gr)
+ w = (GR_N(gr) * GR_SA(gr) + GR_TYPE(gr) * sin (w)) /
+ (GR_G(gr) * m)
+ return (w)
+end
+
+
+# GR_W2X -- Position at given wavelength.
+
+real procedure gr_w2x (gr, w, m)
+
+pointer gr #I Grating pointer
+real w #I Wavelength (Angstroms)
+int m #I Order
+real x #I Pixel position (mm from center)
+
+begin
+ if (GR_FULL(gr) == NO) {
+ x = (w - GR_WB(gr)) * m / GR_DB(gr)
+ return (x)
+ }
+
+ x = (w * m * GR_G(gr) - GR_N(gr) * GR_SA(gr)) / GR_TYPE(gr)
+ x = asin (x) - GR_BETA(gr)
+ x = x * GR_F(gr)
+ return (x)
+end
+
+
+# GR_MAG -- Grating magnification = cos (alpha) / cos (beta(w,m))
+
+real procedure gr_mag (gr, w, m)
+
+pointer gr #I Grating pointer
+real w #I Wavelength
+int m #I Order
+real mag #O mag
+
+real gr_sinbeta()
+
+begin
+ mag = gr_sinbeta (gr, w, m)
+ if (!IS_INDEF(mag))
+ mag = GR_CA(gr) / sqrt (1 - mag * mag)
+ return (mag)
+end
+
+
+# GR_TILT -- Grating tilt = (alpha + beta) / 2
+
+real procedure gr_tilt (gr, w, m)
+
+pointer gr #I Grating pointer
+real w #I Wavelength
+int m #I Order
+real tilt #O tilt
+
+real gr_sinbeta()
+
+begin
+ tilt = gr_sinbeta (gr, w, m)
+ if (!IS_INDEF(tilt))
+ tilt = (GR_ALPHA(gr) + GR_TYPE(gr)*asin(tilt)) / 2
+ return (tilt)
+end
+
+
+# GR_SINBETA -- sin(beta(w,m)) = g * m * w - n * sin(alpha)
+# n is index of refraction which is different from 1 for a grism.
+
+real procedure gr_sinbeta (gr, w, m)
+
+pointer gr #I Grating pointer
+real w #I Wavelength
+int m #I Order
+real sb #O sin(beta)
+
+begin
+ if (gr == NULL)
+ return (INDEF)
+ if (GR_FULL(gr) == NO)
+ return (INDEF)
+
+ sb = (GR_G(gr) * m * w - GR_N(gr) * GR_SA(gr)) / GR_TYPE(gr)
+ if (abs(sb) >= 1.)
+ sb = INDEF
+
+ return (sb)
+end
diff --git a/noao/obsutil/src/sptime/lib/abjohnson b/noao/obsutil/src/sptime/lib/abjohnson
new file mode 100644
index 00000000..b2bd6d2a
--- /dev/null
+++ b/noao/obsutil/src/sptime/lib/abjohnson
@@ -0,0 +1,17 @@
+# AB zeropoints.
+# Values are from Astrophysical Quantities converted from
+# log f (in W/cm^2/micron) to AB = -2.5 * log (F_lambda) - 5 log (lambda) - 2.4
+# where lambda is in Angstroms and F_lambda is in ergs/s/cm^2/A. The
+# values are for U, B, V, R, I, J, H, Ks, K, L, L', M.
+
+ 3650 0.714
+ 4400 -0.167
+ 5500 -0.052
+ 7000 0.275
+ 9000 0.529
+12150 0.878
+16540 1.356
+21570 1.857
+35470 2.803
+37610 2.921
+47690 3.397
diff --git a/noao/obsutil/src/sptime/lib/circle b/noao/obsutil/src/sptime/lib/circle
new file mode 100644
index 00000000..ff3ca32e
--- /dev/null
+++ b/noao/obsutil/src/sptime/lib/circle
@@ -0,0 +1,21 @@
+# title = "Circular aperture"
+# type = "circular"
+
+0.00 0.000
+0.10 0.007
+0.13 0.011
+0.16 0.017
+0.20 0.027
+0.25 0.043
+0.32 0.067
+0.40 0.104
+0.50 0.160
+0.63 0.241
+0.79 0.354
+1.00 0.500
+1.26 0.667
+1.58 0.825
+2.00 0.937
+2.51 0.987
+3.16 0.999
+3.98 1.000
diff --git a/noao/obsutil/src/sptime/lib/slit b/noao/obsutil/src/sptime/lib/slit
new file mode 100644
index 00000000..2c39234a
--- /dev/null
+++ b/noao/obsutil/src/sptime/lib/slit
@@ -0,0 +1,103 @@
+# title = "Slit with Gaussian profile"
+# type = "rectangular"
+
+0.00 0.00 0.000
+0.13 0.00 0.000
+0.21 0.00 0.000
+0.32 0.00 0.000
+0.52 0.00 0.000
+0.82 0.00 0.000
+1.30 0.00 0.000
+2.05 0.00 0.000
+3.26 0.00 0.000
+5.17 0.00 0.000
+0.00 0.13 0.000
+0.13 0.13 0.009
+0.21 0.13 0.014
+0.32 0.13 0.022
+0.52 0.13 0.034
+0.82 0.13 0.051
+1.30 0.13 0.071
+2.05 0.13 0.088
+3.26 0.13 0.093
+5.17 0.13 0.094
+0.00 0.21 0.000
+0.13 0.21 0.014
+0.21 0.21 0.022
+0.32 0.21 0.034
+0.52 0.21 0.053
+0.82 0.21 0.080
+1.30 0.21 0.113
+2.05 0.21 0.139
+3.26 0.21 0.148
+5.17 0.21 0.148
+0.00 0.32 0.000
+0.13 0.32 0.022
+0.21 0.32 0.034
+0.32 0.32 0.054
+0.52 0.32 0.084
+0.82 0.32 0.126
+1.30 0.32 0.177
+2.05 0.32 0.218
+3.26 0.32 0.232
+5.17 0.32 0.233
+0.00 0.52 0.000
+0.13 0.52 0.034
+0.21 0.52 0.053
+0.32 0.52 0.084
+0.52 0.52 0.130
+0.82 0.52 0.196
+1.30 0.52 0.275
+2.05 0.52 0.338
+3.26 0.52 0.360
+5.17 0.52 0.361
+0.00 0.82 0.000
+0.13 0.82 0.051
+0.21 0.82 0.080
+0.32 0.82 0.126
+0.52 0.82 0.196
+0.82 0.82 0.294
+1.30 0.82 0.413
+2.05 0.82 0.509
+3.26 0.82 0.541
+5.17 0.82 0.542
+0.00 1.30 0.000
+0.13 1.30 0.071
+0.21 1.30 0.113
+0.32 1.30 0.177
+0.52 1.30 0.275
+0.82 1.30 0.413
+1.30 1.30 0.579
+2.05 1.30 0.714
+3.26 1.30 0.759
+5.17 1.30 0.761
+0.00 2.05 0.000
+0.13 2.05 0.088
+0.21 2.05 0.139
+0.32 2.05 0.218
+0.52 2.05 0.338
+0.82 2.05 0.509
+1.30 2.05 0.714
+2.05 2.05 0.880
+3.26 2.05 0.935
+5.17 2.05 0.938
+0.00 3.26 0.000
+0.13 3.26 0.093
+0.21 3.26 0.148
+0.32 3.26 0.232
+0.52 3.26 0.360
+0.82 3.26 0.541
+1.30 3.26 0.759
+2.05 3.26 0.935
+3.26 3.26 0.994
+5.17 3.26 0.997
+0.00 5.17 0.000
+0.13 5.17 0.094
+0.21 5.17 0.148
+0.32 5.17 0.233
+0.52 5.17 0.361
+0.82 5.17 0.542
+1.30 5.17 0.761
+2.05 5.17 0.938
+3.26 5.17 0.997
+5.17 5.17 1.000
diff --git a/noao/obsutil/src/sptime/mkcircle.cl b/noao/obsutil/src/sptime/mkcircle.cl
new file mode 100644
index 00000000..2fc61210
--- /dev/null
+++ b/noao/obsutil/src/sptime/mkcircle.cl
@@ -0,0 +1,16 @@
+# MKCIRCLE -- Fraction of Gaussian profile going through a circular aperture
+# of specified diameter in units of the profile FWHM.
+
+real d, logd
+
+printf ("## CIRCLE -- Fraction of Gaussian profile going through a")
+printf (" cicular\n## aperture as a function of the diameter in units")
+printf (" of FWHM.\n\n")
+
+printf ("%4.2f %5.3f\n", 0, 0)
+for (logd=-1; logd<=0.6; logd=logd+0.1) {
+ d = 10.**logd
+ z = d * 0.8325546111577
+ z = 1 - exp (-(z * z))
+ printf ("%4.2f %5.3f\n", d, z)
+}
diff --git a/noao/obsutil/src/sptime/mkpkg b/noao/obsutil/src/sptime/mkpkg
new file mode 100644
index 00000000..21f5845d
--- /dev/null
+++ b/noao/obsutil/src/sptime/mkpkg
@@ -0,0 +1,20 @@
+# Make SPECTIME.
+
+$checkout libpkg.a ../
+$update libpkg.a
+$checkin libpkg.a ../
+$exit
+
+standalone:
+ $update libpkg.a
+ $omake x_spectime.x
+ $link x_spectime.o libpkg.a -lsmw -liminterp -o xx_spectime.e
+ ;
+
+libpkg.a:
+ grating.x <error.h> <math.h>
+ stdisperser.x sptime.h
+ t_cgiparse.x
+ t_sptime.x sptime.h <error.h> <gset.h> <math.h> <ctype.h> <mach.h>
+ tabinterp.x <error.h> <math/iminterp.h> <mach.h>
+ ;
diff --git a/noao/obsutil/src/sptime/mkslit.cl b/noao/obsutil/src/sptime/mkslit.cl
new file mode 100644
index 00000000..240f851e
--- /dev/null
+++ b/noao/obsutil/src/sptime/mkslit.cl
@@ -0,0 +1,37 @@
+# MKSLIT -- Fraction of Gaussian profile going through a slit aperture
+# of specified width and height in units of the profile FWHM.
+
+real logx, logy, t, erfcc, xval, yval
+
+printf ("## SLIT -- Fraction of Gaussian profile going through a")
+printf (" slit\n## aperture as a function of the width and height")
+printf (" in units of FWHM.\n\n")
+
+printf ("%4.2f %4.2f %5.3f\n", 0, 0, 0)
+for (logx=-1; logx<=0.6; logx=logx+0.2) {
+ x = 10.**logx
+ printf ("%4.2f %4.2f %5.3f\n", x, 0, 0)
+}
+for (logy=-1; logy<=0.6; logy=logy+0.2) {
+ y = 10.**logy
+ z = y * 0.8325546111577
+ t = 1. / (1. + 0.5 * z)
+ erfcc = t * exp (-z * z - 1.26551223 + t * (1.00002368 +
+ t * (0.37409196 + t * (0.09678418 + t * (-0.18628806 +
+ t * (0.27886807 + t * (-1.13520398 + t * (1.48851587 +
+ t * (-0.82215223 + t * 0.17087277)))))))))
+ yval = (1 - erfcc)
+ printf ("%4.2f %4.2f %5.3f\n", 0, y, 0)
+ for (logx=-1; logx<=0.6; logx=logx+0.2) {
+ x = 10.**logx
+ z = x * 0.8325546111577
+ t = 1. / (1. + 0.5 * z)
+ erfcc = t * exp (-z * z - 1.26551223 + t * (1.00002368 +
+ t * (0.37409196 + t * (0.09678418 + t * (-0.18628806 +
+ t * (0.27886807 + t * (-1.13520398 + t * (1.48851587 +
+ t * (-0.82215223 + t * 0.17087277)))))))))
+ xval = (1 - erfcc)
+ z = xval * yval
+ printf ("%4.2f %4.2f %5.3f\n", x, y, z)
+ }
+}
diff --git a/noao/obsutil/src/sptime/rates.cl b/noao/obsutil/src/sptime/rates.cl
new file mode 100644
index 00000000..f664b097
--- /dev/null
+++ b/noao/obsutil/src/sptime/rates.cl
@@ -0,0 +1,74 @@
+procedure rates (grating, filter, order)
+
+file grating {prompt="Grating"}
+file filter {prompt="Filter"}
+int order = INDEF {prompt="Order"}
+real wave = INDEF {prompt="Central wavelength"}
+file spectrograph = "rcspec" {prompt="Spectrograph"}
+real width = 10. {prompt="Slit width"}
+
+string title = "" {prompt="Title"}
+real w1 = 3000. {prompt="Lower wavelength to plot"}
+real w2 = 12000. {prompt="Upper wavelength to plot\n"}
+
+real x1 = 3000. {prompt="Left graph wavelength"}
+real x2 = 12000. {prompt="Right graph wavelength"}
+real y1 = -5. {prompt="Bottom graph efficiency"}
+real y2 = 105. {prompt="Top graph efficiency"}
+string ltype = "1" {prompt="Line type"}
+string color = "1" {prompt="Color"}
+bool append = no {prompt="Append?"}
+
+struct *fd
+
+begin
+ string search
+ file tmp1, tmp2
+ real x, y
+
+ tmp1 = mktemp ("tmp$iraf")
+ tmp2 = mktemp ("tmp$iraf")
+
+ search = "spectimedb$,sptimeKPNO$"
+ search = search // ",sptimeKPNO$Spectrographs"
+ search = search // ",sptimeKPNO$Gratings"
+ search = search // ",sptimeKPNO$Grisms"
+ search = search // ",sptimeKPNO$Filters/RCCRYO"
+ search = search // ",sptimeKPNO$Filters/GCAM"
+
+ sptime (time=1., maxexp=3600., sn=25., spectrum="fnu_power",
+ sky="none", sensfunc="none", airmass=1., seeing=1.5, phase=0.,
+ temperature=6000., index=0., refwave=INDEF, refflux=10.,
+ funits="AB", abjohnson="none", wave=wave, order=order,
+ xorder=INDEF, width=width, length=INDEF, diameter=INDEF,
+ inoutangle=INDEF, xinoutangle=INDEF, xbin=1, ybin=1,
+ search=search, spectrograph=spectrograph, filter=filter,
+ filter2="none", disperser=grating, xdisperser="none",
+ fiber="none", telescope="", adc="", collimator="",
+ corrector="", camera="", detector="", aperture="",
+ extinction="", gain=INDEF, rdnoise=INDEF, dark=INDEF, skysub="none",
+ nskyaps=10, output="rate", list=tmp2, graphics="",
+ interactive=no, nw=1000, > "dev$null")
+
+ fd = tmp2
+ while (fscan (fd, x, y) != EOF) {
+ if (nscan() != 2)
+ next
+ if (x < w1 || x > w2)
+ next
+ print (x, y, >> tmp1)
+ }
+ fd = ""
+
+ graph (tmp1, wx1=x1, wx2=x2, wy1=y1, wy2=y2, wcs="logical", axis=1,
+ transpose=no, pointmode=no, marker="box", szmarker=0.005,
+ ltypes=ltype, colors=color, logx=no, logy=no, box=yes,
+ ticklabels=yes, xlabel="Wavelength (A)", ylabel="Rate",
+ xformat="wcsformat", yformat="", title=title,
+ lintran=no, p1=0., p2=0., q1=0., q2=1., vx1=0., vx2=0., vy1=0.,
+ vy2=0., majrx=7, minrx=3, majry=7, minry=3, overplot=no,
+ append=append, device="stdgraph", round=no, fill=yes)
+
+ delete (tmp1, verify-)
+ delete (tmp2, verify-)
+end
diff --git a/noao/obsutil/src/sptime/specpars.par b/noao/obsutil/src/sptime/specpars.par
new file mode 100644
index 00000000..79c6c7f0
--- /dev/null
+++ b/noao/obsutil/src/sptime/specpars.par
@@ -0,0 +1,85 @@
+spectrograph,s,h,"",,,Spectrograph transmission or table
+title,s,h,"",,,Spectrograph title
+apmagdisp,r,h,INDEF,,,Aperture magnification along dispersion
+apmagxdisp,r,h,INDEF,,,Aperture magnification across dispersion
+inoutangle,r,h,INDEF,,,Incident to diffracted angle (deg)
+xinoutangle,r,h,INDEF,,,"Incident to diffracted cross disperser angle (deg)
+
+# Telescope"
+telescope,s,h,"",,,Telescope transmission or table
+teltitle,s,h,"",,,Telescope title
+area,r,h,INDEF,,,Telescope effective area (m^2)
+scale,r,h,INDEF,,,Telescope scale at entrance aperture (arcsec/mm)
+emissivity,s,h,"",,,Emissivity value or table
+emistitle,s,h,"",,,"Emissivity title
+
+# Correctors"
+corrector,s,h,"",,,Corrector transmission or table
+cortitle,s,h,"",,,Corrector title
+adc,s,h,"",,,ADC transmission or table
+adctitle,s,h,"",,,"ADC title
+
+# Disperser"
+disperser,s,h,"",,,Disperser efficiency or table
+disptitle,s,h,"",,,Disperser title
+disptype,s,h,"",,,Disperser type (grating|grism|generic)
+gmm,r,h,INDEF,,,Grating (grooves/mm)
+blaze,r,h,INDEF,,,Blaze or prism angle (deg)
+oref,r,h,INDEF,,,Reference order
+wavelength,r,h,INDEF,,,Central wavelength in ref order (A)
+dispersion,r,h,INDEF,,,Dispersion at central wavelength in ref order (A/mm)
+indexref,r,h,INDEF,,,Index of refraction at first order
+eff,r,h,INDEF,,,"Peak efficiency
+
+# Crossdisperser"
+xdisperser,s,h,"",,,Crossdisperser efficiency or table
+xdisptitle,s,h,"",,,Crossdisperser title
+xdisptype,s,h,"",,,Disperser type (grating|grism|generic)
+xgmm,r,h,INDEF,,,Grating (grooves/mm)
+xblaze,r,h,INDEF,,,Blaze or prism angle (deg)
+xoref,r,h,INDEF,,,Reference order
+xwavelength,r,h,INDEF,,,Central wavelength in ref order (A)
+xindexref,r,h,INDEF,,,Index of refraction at first order
+xdispersion,r,h,INDEF,,,Dispersion at central wavelength in ref order (A/mm)
+xeff,r,h,INDEF,,,"Peak efficiency
+
+# Aperture"
+aperture,s,h,"",,,Aperture transmission or table
+aptitle,s,h,"",,,Aperture title
+aptype,s,h,"",,,"Aperture type (rectanglular|circular)
+
+# Fiber"
+fiber,s,h,"",,,Fiber transmission or table
+fibtitle,s,h,"",,,"Fiber title
+
+# Filters"
+filter,s,h,"",,,Filter transmission or table
+ftitle,s,h,"",,,Filter title
+filter2,s,h,"",,,Filter transmission or table
+f2title,s,h,"",,,Filter title
+block,s,h,"",,,"Assume other orders are blocked (yes|no)
+
+# Collimator"
+collimator,s,h,"",,,Collimator transmission or table
+coltitle,s,h,"",,,Collimator title
+colfl,r,h,INDEF,,,"Collimator focal length (m)
+
+# Camera"
+camera,s,h,"",,,Camera transmission or table
+camtitle,s,h,"",,,Camera title
+camfl,r,h,INDEF,,,Camera focal length (m)
+resolution,r,h,INDEF,,,Camera resolution (mm)
+vignetting,s,h,"",,,"Vignetting table
+
+# Detector"
+detector,s,h,"",,,Detector DQE or table
+dettitle,s,h,"",,,Detector title
+ndisp,r,h,INDEF,,,Detector pixels along dispersion
+pixsize,r,h,INDEF,,,Detector pixel size (mm)
+gain,r,h,INDEF,,,Detector gain (photons/ADU)
+rdnoise,r,h,INDEF,,,Detector read noise (photons)
+dark,r,h,INDEF,,,Detector dark count rate (photons/s)
+saturation,r,h,INDEF,,,Detector saturation (photons)
+dnmax,r,h,INDEF,,,Detector maximum counts (ADU)
+xbin,i,h,1,1,,Detector binning (dispersion)
+ybin,i,h,1,1,,Detector binning (spatial)
diff --git a/noao/obsutil/src/sptime/sptime.h b/noao/obsutil/src/sptime/sptime.h
new file mode 100644
index 00000000..c17b377a
--- /dev/null
+++ b/noao/obsutil/src/sptime/sptime.h
@@ -0,0 +1,132 @@
+# Definitions for SPTIME.
+
+# Spectral distribution types.
+define SPECTYPES "|blackbody|flambda_power|fnu_power|"
+define SPEC_TAB 0
+define SPEC_BB 1
+define SPEC_FL 2
+define SPEC_FN 3
+
+# Flux units.
+define FUNITS "|AB|F_lambda|F_nu|U|B|V|R|I|J|H|Ks|K|L|L'|M|"
+define AB 1
+define FLAMBDA 2
+define FNU 3
+define UMAG 4
+define BMAG 5
+define VMAG 6
+define RMAG 7
+define IMAG 8
+define JMAG 9
+define HMAG 10
+define KSMAG 11
+define KMAG 12
+define LMAG 13
+define LPMAG 14
+define MMAG 15
+
+# Sky subtraction options.
+define SKYSUB "|none|longslit|multiap|shuffle|"
+define SKY_NONE 1
+define SKY_LONGSLIT 2
+define SKY_MULTIAP 3
+define SKY_SHUFFLE 4
+
+# Aperture types.
+define APTYPES "|circular|rectangular|"
+define CIRCULAR 1
+define RECTANGULAR 2
+
+# Output types.
+define OUTTYPES "|counts|snr|object|rate|atmosphere|telescope|adc|\
+ |aperture|fiber|filter|filter2|collimator|disperser|\
+ |xdisperser|corrector|camera|detector|spectrograph|\
+ |emissivity|thruput|sensfunc|correction|"
+define OUT_COUNTS 1
+define OUT_SNR 2
+define OUT_OBJ 3
+define OUT_RATE 4
+define OUT_ATM 5
+define OUT_TEL 6
+define OUT_ADC 7
+define OUT_AP 9
+define OUT_FIB 10
+define OUT_FILT 11
+define OUT_FILT2 12
+define OUT_COL 13
+define OUT_DISP 14
+define OUT_XDISP 16
+define OUT_COR 17
+define OUT_CAM 18
+define OUT_DET 19
+define OUT_SPEC 20
+define OUT_EMIS 22
+define OUT_THRU 23
+define OUT_SENS 24
+define OUT_CORRECT 25
+
+# Grating types.
+define DISPTYPES "|grating|grism|generic"
+define GRATING 1
+define GRISM 2
+define GENERIC 3
+
+# Macros.
+define MINEXP 0.01 # Minimum exposure time.
+define MAXNEXP 10000 # Maximum number of exposures.
+define MAXITER 50
+define H 6.626E-27
+define C 2.99792456E18
+
+# Data structure.
+define ST_SZSTRING 99 # Length of strings
+define ST_LEN 154 # Structure length
+
+define ST_TAB Memi[$1] # Tables pointer
+define ST_SEARCH Memi[$1+1] # Search list
+define ST_SPEC Memi[$1+2] # Spectrum type
+define ST_PARAM Memr[P2R($1+3)] # Spectrum parameter
+define ST_THERMAL Memr[P2R($1+4)] # Thermal background temperature
+define ST_AV Memr[P2R($1+5)] # A(V)
+define ST_RV Memr[P2R($1+6)] # A(V)/E(B-V)
+define ST_TIME Memr[P2R($1+7)] # Exposure time
+define ST_NEXP Memi[$1+8] # Number of exposures
+define ST_MINEXP Memr[P2R($1+9)] # Minimum time per integration
+define ST_MAXEXP Memr[P2R($1+10)] # Maximum time per integration
+define ST_AIRMASS Memr[P2R($1+11)] # Airmass
+define ST_SEEING Memr[P2R($1+12)] # Seeing (FWHM in arcsec)
+define ST_PHASE Memr[P2R($1+13)] # Moon phase
+define ST_REFW Memr[P2R($1+14)] # Reference wavelength
+define ST_REFF Memr[P2R($1+15)] # Reference flux
+define ST_REFFL Memr[P2R($1+16)] # Reference flambda
+define ST_CW Memr[P2R($1+17)] # Central wavelength
+define ST_ORDER Memi[$1+17+$2] # Grating orders (2)
+define ST_APSIZE Memr[P2R($1+19+$2)] # Aperture sizes (2)
+define ST_BIN Memi[$1+21+$2] # Binning (2)
+define ST_GAIN Memr[P2R($1+24)] # Detector gain
+define ST_RDNOISE Memr[P2R($1+25)] # Detector read noise
+define ST_DARK Memr[P2R($1+26)] # Detector dark counts
+define ST_INOUTA Memr[P2R($1+26+$2)] # Incident-diffracted angle(deg)
+define ST_SKYSUB Memi[$1+29] # Sky subtraction type
+define ST_NSKYAPS Memi[$1+30] # Number of sky apertures
+
+define ST_DISPTYPE Memi[$1+30+$2] # Disperser type
+define ST_GR Memi[$1+32+$2] # Grating pointer (2)
+define ST_DISP Memr[P2R($1+34+$2)] # Dispersion (2)
+define ST_SCALE Memr[P2R($1+36+$2)] # Scale at detector (2)
+define ST_RES Memr[P2R($1+38+$2)] # Resolution at detector (2)
+define ST_AREA Memr[P2R($1+41)] # Effective collecting area
+define ST_TELSCALE Memr[P2R($1+42)] # Telescope scale at aperture
+define ST_NOBJPIX Memi[$1+43] # Number of object pixels
+define ST_NSKYPIX Memi[$1+44] # Number of sky pixels
+define ST_APLIM Memi[$1+45] # Aperture limited?
+define ST_APTYPE Memi[$1+46] # Aperture type
+define ST_PIXSIZE Memr[P2R($1+47)] # Pixel weighted disp. size
+define ST_FILTBLK Memi[$1+48] # Block filter flag
+define ST_DUNANG Memi[$1+49] # Angstrom units pointer
+define ST_DUN Memi[$1+50] # User units pointer
+define ST_SUBPIXELS Memi[$1+51] # Number of subpixels
+define ST_COLFL Memr[P2R($1+52)] # Collimator focal length
+define ST_CAMFL Memr[P2R($1+53)] # Collimator focal length
+define ST_DUNITS Memc[P2C($1+54)] # Dispersion units
+define ST_FUNITS Memc[P2C($1+104)] # Flux units
diff --git a/noao/obsutil/src/sptime/sptime.par b/noao/obsutil/src/sptime/sptime.par
new file mode 100644
index 00000000..0342e4f2
--- /dev/null
+++ b/noao/obsutil/src/sptime/sptime.par
@@ -0,0 +1,53 @@
+time,r,h,INDEF,,,Total exposure time (sec)
+sn,r,h,25.,,,"S/N per pixel if time is INDEF
+
+# Source and background parameters"
+spectrum,s,h,"blackbody",,,Source spectrum
+spectitle,s,h,"",,,Spectrum title
+E,r,h,0.,,,E(B-V)
+R,r,h,3.1,,,A(V)/E(B-V)
+sky,s,h,"",,,Sky spectrum
+skytitle,s,h,"",,,Sky title
+extinction,s,h,"",,,Extinction value or table
+exttitle,s,h,"",,,"Extinction title
+
+# Spectral distribution parameters"
+refwave,r,h,INDEF,,,Reference wavelength
+refflux,r,h,10,,,Source flux at reference wavelength
+funits,s,h,"AB","AB|F_lambda|F_nu|U|B|V|R|I|J|H|Ks|K|L|L'|M|",,Reference flux units
+temperature,r,h,6000.,,,Blackbody temperature (K)
+index,r,h,0.,,,"Power law index
+
+# Observation parameters"
+seeing,r,h,1.,,,Seeing (arcsec)
+airmass,r,h,1.,0.,,Airmass
+phase,r,h,0.,0.,14.,Moon phase (0-14)
+thermal,r,h,0.,,,Thermal background temperature (K)
+wave,r,h,INDEF,,,Central wavelength
+order,i,h,INDEF,,,Order for grating
+xorder,i,h,INDEF,,,Order for crossdisperser grating
+width,r,h,INDEF,,,Slit width (-=pixels +=arcsec)
+length,r,h,INDEF,,,Slit length (-=pixels +=arcsec)
+diameter,r,h,INDEF,,,Circular/fiber aperture diameter (-=pixels +=arcsec)
+xbin,i,h,1,1,,Detector binning (dispersion)
+ybin,i,h,1,1,,"Detector binning (spatial)
+
+# Calculation parameters"
+search,s,h,"spectimedb$",,,List of directories to search for tables
+minexp,r,h,0.01,,,Minimum time per exposure (sec)
+maxexp,r,h,3600.,,,Maximum time per exposure (sec)
+units,s,h,"Angstroms",,,Wavelength units
+skysub,s,h,"",,,Type of sky subtaction (none|longslit|multiap|shuffle)
+nskyaps,i,h,10,0,,Number of multiap sky apertures
+subpixels,i,h,INDEF,,,Number of subpixels in calculation
+sensfunc,s,h,"",,,"Sensitivity function value or table
+
+# Output parameters"
+output,s,h,"",,,List of output quantities
+nw,i,h,101,,,Number of output dispersion points
+list,s,h,"",,,Output list file
+graphics,s,h,"stdgraph",,,Graphics device
+interactive,b,h,yes,,,"Interactive pause after graphs?
+
+# Spectrograph parameters"
+specpars,pset,h,"",,,Spectrograph parameters
diff --git a/noao/obsutil/src/sptime/stdisperser.x b/noao/obsutil/src/sptime/stdisperser.x
new file mode 100644
index 00000000..090cf010
--- /dev/null
+++ b/noao/obsutil/src/sptime/stdisperser.x
@@ -0,0 +1,455 @@
+include "sptime.h"
+
+# These routines interface between any type of disperser and the
+# disperser type specific routines (such as those for gratings).
+
+
+# ST_DISPERSER -- Initialize disperser data.
+
+procedure st_disperser (st, name, index)
+
+pointer st #I SPECTIME pointer
+char name[ARB] #I Table name
+int index #I Grating index
+
+int order, oref, stgeti(), tabgeti(), strdic()
+real f, phi, g, blaze, wb, db, ref, stgetr(), tabgetr(), gr_getr()
+pointer sp, str, fname, gr, gr_open()
+bool streq()
+errchk gr_open, tab_getr, gr_getr, st_gtable1
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+ call salloc (fname, SZ_LINE, TY_CHAR)
+
+ ST_GR(st,index) = NULL
+ ST_DISPTYPE(st,index) = 0
+
+ # Get disperser type.
+ if (streq (name, "disperser")) {
+ call stgstr (st, "disptype", "disperser", "", Memc[str], SZ_LINE)
+ if (Memc[str] == EOS) {
+ iferr (call tabgstr (ST_TAB(st), "disperser", "spectrograph",
+ "type", Memc[str], SZ_LINE))
+ call strcpy ("generic", Memc[str], SZ_LINE)
+ }
+ } else if (streq (name, "xdisperser")) {
+ call stgstr (st, "xdisptype", "xdisperser", "", Memc[str], SZ_LINE)
+ if (Memc[str] == EOS) {
+ iferr (call tabgstr (ST_TAB(st), "xdisperser", "spectrograph",
+ "type", Memc[str], SZ_LINE))
+ Memc[str] = EOS
+ }
+ } else
+ Memc[str] = EOS
+ ST_DISPTYPE(st,index) = strdic (Memc[str],Memc[str],SZ_LINE,DISPTYPES)
+
+ switch (ST_DISPTYPE(st,index)) {
+ case GRATING:
+ f = ST_CAMFL(st) * 1000
+ switch (index) {
+ case 1:
+ g = stgetr (st, "gmm", name, INDEFR)
+ blaze = stgetr (st, "blaze", name, INDEFR)
+ oref = stgeti (st, "oref", name, 1)
+ wb = stgetr (st, "wavelength", name, INDEFR)
+ db = stgetr (st, "dispersion", name, INDEFR)
+ ref = stgetr (st, "eff", name, INDEFR)
+ case 2:
+ g = stgetr (st, "xgmm", name, INDEFR)
+ blaze = stgetr (st, "xblaze", name, INDEFR)
+ oref = stgeti (st, "xoref", name, 1)
+ wb = stgetr (st, "xwavelength", name, INDEFR)
+ db = stgetr (st, "xdispersion", name, INDEFR)
+ ref = stgetr (st, "xeff", name, INDEFR)
+
+ # Check old names.
+ if (IS_INDEFR(g)) {
+ iferr (g = tabgetr (ST_TAB(st), name, "spectrograph",
+ "gmm"))
+ g = INDEFR
+ }
+ }
+
+ # Old names.
+ if (IS_INDEFR(g)) {
+ iferr (g = tabgetr (ST_TAB(st), name, "spectrograph",
+ "gmm"))
+ g = INDEFR
+ }
+ if (IS_INDEFR(blaze)) {
+ iferr (blaze = tabgetr (ST_TAB(st), name, "spectrograph",
+ "blaze"))
+ blaze = INDEFR
+ }
+ if (IS_INDEFI(oref)) {
+ iferr (oref = tabgeti (ST_TAB(st), name, "spectrograph",
+ "oref"))
+ oref = 1
+ }
+ if (IS_INDEFR(wb)) {
+ iferr (wb = tabgetr (ST_TAB(st), name, "spectrograph",
+ "wavelength"))
+ wb = INDEFR
+ }
+ if (IS_INDEFR(db)) {
+ iferr (db = tabgetr (ST_TAB(st), name, "spectrograph",
+ "dispersion"))
+ db = INDEFR
+ }
+ if (IS_INDEFR(ref)) {
+ iferr (ref = tabgetr (ST_TAB(st), name, "spectrograph",
+ "reflectance"))
+ ref = 1.
+ }
+
+ phi = ST_INOUTA(st,index)
+# if (!IS_INDEFR(db))
+# db = db / f
+
+ iferr (gr = gr_open (ST_CW(st), ST_ORDER(st,index), ref, wb, db,
+ oref, f, g, blaze, 1., phi, INDEF, INDEF, 1, YES)) {
+ g = 300.
+ blaze = 6
+ gr = gr_open (ST_CW(st), ST_ORDER(st,index), ref, wb, db,
+ oref, f, g, blaze, 1., phi, INDEF, INDEF, 1, YES)
+ }
+ ST_GR(st,index) = gr
+
+ if (IS_INDEF(ST_CW(st)))
+ ST_CW(st,index) = gr_getr (gr, "wavelength")
+ if (IS_INDEFI(ST_ORDER(st,index)))
+ ST_ORDER(st,index) = nint (gr_getr (gr, "order"))
+ ST_DISP(st,index) = gr_getr (gr, "dblaze") * oref * f
+
+ # Look for explicit blaze functions.
+ do order = ST_ORDER(st,index)-1, ST_ORDER(st,index)+1 {
+ call sprintf (Memc[str], SZ_LINE, "eff%d")
+ call pargi (order)
+ ifnoerr (call tabgstr (ST_TAB(st), name, "spectrograph",
+ Memc[str], Memc[fname], SZ_LINE)) {
+ if (streq (name, "disperser"))
+ call st_gtable1 (st, Memc[str], Memc[fname])
+ else if (streq (name, "xdisperser")) {
+ call sprintf (Memc[str], SZ_LINE, "xeff%d")
+ call pargi (order)
+ call st_gtable1 (st, Memc[str], Memc[fname])
+ }
+ }
+ }
+ case GRISM:
+ f = ST_CAMFL(st) * 1000
+ switch (index) {
+ case 1:
+ g = stgetr (st, "gmm", name, INDEFR)
+ blaze = stgetr (st, "blaze", name, INDEFR)
+ ref = stgetr (st, "eff", name, 1.)
+ db = stgetr (st, "indexref", name, INDEFR)
+ if (!IS_INDEFI(ST_ORDER(st,index)) && ST_ORDER(st,index)!=1) {
+ call sprintf (Memc[str], SZ_LINE, "index%d")
+ call pargi (ST_ORDER(st,index))
+ iferr (wb = tabgetr (ST_TAB(st), name, "spectrograph",
+ Memc[str]))
+ wb = db
+ db = wb
+ }
+ case 2:
+ g = stgetr (st, "xgmm", name, INDEFR)
+ blaze = stgetr (st, "xblaze", name, INDEFR)
+ ref = stgetr (st, "xeff", name, 1.)
+ db = stgetr (st, "xindexref", name, INDEFR)
+ if (!IS_INDEFI(ST_ORDER(st,index)) && ST_ORDER(st,index)!=1) {
+ call sprintf (Memc[str], SZ_LINE, "index%d")
+ call pargi (ST_ORDER(st,index))
+ iferr (wb = tabgetr (ST_TAB(st), name, "spectrograph",
+ Memc[str]))
+ wb = db
+ db = wb
+ }
+ }
+
+ # Old names.
+ if (IS_INDEFR(g)) {
+ iferr (g = tabgetr (ST_TAB(st), name, "spectrograph",
+ "gmm"))
+ g = INDEFR
+ }
+ if (IS_INDEFR(blaze)) {
+ iferr (blaze = tabgetr (ST_TAB(st), name, "spectrograph",
+ "prism"))
+ blaze = INDEFR
+ }
+ if (IS_INDEFR(ref)) {
+ iferr (ref = tabgetr (ST_TAB(st), name, "spectrograph",
+ "transmission"))
+ ref = 1
+ }
+ if (IS_INDEFR(db)) {
+ if (!IS_INDEFI(ST_ORDER(st,index))) {
+ call sprintf (Memc[str], SZ_LINE, "index%d")
+ call pargi (ST_ORDER(st,index))
+ iferr (db = tabgetr (ST_TAB(st), name, "spectrograph",
+ Memc[str]))
+ db = tabgetr (ST_TAB(st), name, "spectrograph",
+ "index1")
+ } else
+ db = tabgetr (ST_TAB(st), name, "spectrograph", "index1")
+ }
+ oref = 1
+
+ iferr (gr = gr_open (ST_CW(st), ST_ORDER(st,index), ref, INDEF,
+ INDEF, oref, f, g, blaze, db, 0., blaze, blaze, 1, YES)) {
+ g = 300.
+ blaze = 6.
+ gr = gr_open (ST_CW(st), ST_ORDER(st,index), ref, INDEF,
+ INDEF, oref, f, g, blaze, db, 0., blaze, blaze, 1, YES)
+ }
+
+ ST_GR(st,index) = gr
+
+ if (IS_INDEF(ST_CW(st)))
+ ST_CW(st,index) = gr_getr (gr, "wavelength")
+ if (IS_INDEFI(ST_ORDER(st,index)))
+ ST_ORDER(st,index) = nint (gr_getr (gr, "order"))
+ ST_DISP(st,index) = gr_getr (gr, "dblaze") * oref * f
+
+ # Look for explicit blaze functions.
+ do order = ST_ORDER(st,index)-1, ST_ORDER(st,index)+1 {
+ call sprintf (Memc[str], SZ_LINE, "eff%d")
+ call pargi (order)
+ ifnoerr (call tabgstr (ST_TAB(st), name, "spectrograph",
+ Memc[str], Memc[fname], SZ_LINE)) {
+ if (streq (name, "disperser"))
+ call st_gtable1 (st, Memc[str], Memc[fname])
+ else if (streq (name, "xdisperser")) {
+ call sprintf (Memc[str], SZ_LINE, "xeff%d")
+ call pargi (order)
+ call st_gtable1 (st, Memc[str], Memc[fname])
+ }
+ }
+ }
+ case GENERIC:
+ f = ST_CAMFL(st) * 1000
+ g = INDEFR
+ blaze = INDEFR
+ oref = 1
+
+ switch (index) {
+ case 1:
+ g = INDEFR
+ blaze = INDEFR
+ oref = 1
+ wb = stgetr (st, "wavelength", name, INDEFR)
+ db = stgetr (st, "dispersion", name, INDEFR)
+ ref = stgetr (st, "eff", name, INDEFR)
+ case 2:
+ g = INDEFR
+ blaze = INDEFR
+ oref = 1
+ wb = stgetr (st, "xwavelength", name, INDEFR)
+ db = stgetr (st, "xdispersion", name, INDEFR)
+ ref = stgetr (st, "xeff", name, INDEFR)
+ }
+
+ if (IS_INDEFR(wb)) {
+ iferr (wb = tabgetr (ST_TAB(st), name, "spectrograph",
+ "wavelength"))
+ wb = INDEFR
+ }
+ if (IS_INDEFR(db)) {
+ iferr (db = tabgetr (ST_TAB(st), name, "spectrograph",
+ "dispersion"))
+ db = INDEFR
+ }
+ if (IS_INDEFR(ref)) {
+ iferr (ref = tabgetr (ST_TAB(st), name, "spectrograph",
+ "reflectance"))
+ ref = 1.
+ }
+
+ phi = ST_INOUTA(st,index)
+
+ gr = gr_open (ST_CW(st), ST_ORDER(st,index), ref, wb, db,
+ oref, f, g, blaze, 1., phi, INDEF, INDEF, 1, NO)
+ ST_GR(st,index) = gr
+
+ if (IS_INDEF(ST_CW(st)))
+ ST_CW(st,index) = gr_getr (gr, "wavelength")
+ if (IS_INDEFI(ST_ORDER(st,index)))
+ ST_ORDER(st,index) = nint (gr_getr (gr, "order"))
+ ST_DISP(st,index) = gr_getr (gr, "dblaze") * oref * f
+
+ # Look for explicit blaze functions.
+ do order = ST_ORDER(st,index)-1, ST_ORDER(st,index)+1 {
+ call sprintf (Memc[str], SZ_LINE, "eff%d")
+ call pargi (order)
+ ifnoerr (call tabgstr (ST_TAB(st), name, "spectrograph",
+ Memc[str], Memc[fname], SZ_LINE)) {
+ if (streq (name, "disperser"))
+ call st_gtable1 (st, Memc[str], Memc[fname])
+ else if (streq (name, "xdisperser")) {
+ call sprintf (Memc[str], SZ_LINE, "xeff%d")
+ call pargi (order)
+ call st_gtable1 (st, Memc[str], Memc[fname])
+ }
+ }
+ }
+ }
+
+ call sfree (sp)
+end
+
+
+# ST_DISPEFF -- Return disperser efficiency.
+
+real procedure st_dispeff (st, name, wave, order)
+
+pointer st #I SPECTIME pointer
+char name[ARB] #I Table name
+real wave #I Wavelength
+int order #I Order
+real eff #O Efficiency
+
+real w
+pointer sp, tab, str
+
+int tabgeti()
+real tabinterp1(), tabgetr(), gr_eff()
+bool tabexists(), streq()
+errchk tabinterp1, gr_eff
+
+begin
+ tab = ST_TAB(st)
+ eff = INDEF
+
+ if (streq (name, "disperser") && ST_DISPTYPE(st,1) == 0)
+ return (eff)
+ if (streq (name, "xdisperser") && ST_DISPTYPE(st,2) == 0)
+ return (eff)
+
+ call smark (sp)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+ if (streq (name, "disperser")) {
+ call sprintf (Memc[str], SZ_FNAME, "eff%d")
+ call pargi (order)
+ if (!tabexists (tab, Memc[str])) {
+ if (order == 1 && tabexists (tab, name)) {
+ if (tabgeti (tab, name, "", "table.ndim") != 0)
+ call strcpy (name, Memc[str], SZ_FNAME)
+ }
+ }
+ if (tabexists (tab, Memc[str])) {
+ eff = tabinterp1 (tab, Memc[str], wave)
+ w = max (tabgetr (tab, Memc[str], "", "table.xmin"), wave)
+ w = min (tabgetr (tab, Memc[str], "", "table.xmax"), w)
+ if (w != wave) {
+ eff = tabinterp1 (tab, Memc[str], w) /
+ gr_eff (ST_GR(st,1), w, order)
+ eff = eff * gr_eff (ST_GR(st,1), wave, order)
+ }
+ } else
+ eff = gr_eff (ST_GR(st,1), wave, order)
+ } else if (streq (name, "xdisperser")) {
+ call sprintf (Memc[str], SZ_FNAME, "xeff%d")
+ call pargi (order)
+ if (!tabexists (tab, Memc[str])) {
+ if (order == 1 && tabexists (tab, name)) {
+ if (tabgeti (tab, name, "", "table.ndim") != 0)
+ call strcpy (name, Memc[str], SZ_FNAME)
+ }
+ }
+ if (tabexists (tab, Memc[str])) {
+ eff = tabinterp1 (tab, Memc[str], wave)
+ w = max (tabgetr (tab, Memc[str], "", "table.xmin"), wave)
+ w = min (tabgetr (tab, Memc[str], "", "table.xmax"), w)
+ if (w != wave) {
+ eff = tabinterp1 (tab, Memc[str], w) /
+ gr_eff (ST_GR(st,2), w, order)
+ eff = eff * gr_eff (ST_GR(st,2), wave, order)
+ }
+ } else
+ eff = gr_eff (ST_GR(st,2), wave, order)
+ }
+
+ if (IS_INDEF(eff))
+ eff = 0.
+
+ call sfree (sp)
+ return (eff)
+end
+
+
+# ST_X2W -- Return wavelength at given position on detector.
+
+real procedure st_x2w (st, index, x)
+
+pointer st #I SPECTIME pointer
+int index #I Grating index
+real x #I Detector position (mm from center)
+real w #O Wavelength (Angstroms)
+
+real gr_x2w()
+
+begin
+ switch (ST_DISPTYPE(st,index)) {
+ case GRATING:
+ w = gr_x2w (ST_GR(st,index), x, ST_ORDER(st,index))
+ case GRISM:
+ w = gr_x2w (ST_GR(st,index), x, ST_ORDER(st,index))
+ case GENERIC:
+ w = gr_x2w (ST_GR(st,index), x, ST_ORDER(st,index))
+ }
+
+ return (w)
+end
+
+
+# ST_W2X -- Return wavelength at given position on detector.
+
+real procedure st_w2x (st, index, w)
+
+pointer st #I SPECTIME pointer
+int index #I Grating index
+real w #I Wavelength (Angstroms)
+real x #O Detector position (mm from center)
+
+real gr_w2x()
+
+begin
+ switch (ST_DISPTYPE(st,index)) {
+ case GRATING:
+ x = gr_w2x (ST_GR(st,index), w, ST_ORDER(st,index))
+ case GRISM:
+ x = gr_w2x (ST_GR(st,index), w, ST_ORDER(st,index))
+ case GENERIC:
+ x = gr_w2x (ST_GR(st,index), w, ST_ORDER(st,index))
+ }
+
+ return (x)
+end
+
+
+# ST_W2DW -- Return dispersion on detector at given wavelength.
+
+real procedure st_w2dw (st, index, w)
+
+pointer st #I SPECTIME pointer
+int index #I Grating index
+real w #I Wavelength (Angstroms)
+real d #I Dispersion (Angstroms/mm)
+
+real gr_w2dw()
+
+begin
+ switch (ST_DISPTYPE(st,index)) {
+ case GRATING:
+ d = gr_w2dw (ST_GR(st,index), w, ST_ORDER(st,index))
+ case GRISM:
+ d = gr_w2dw (ST_GR(st,index), w, ST_ORDER(st,index))
+ case GENERIC:
+ d = gr_w2dw (ST_GR(st,index), w, ST_ORDER(st,index))
+ }
+
+ return (d)
+end
diff --git a/noao/obsutil/src/sptime/t_cgiparse.x b/noao/obsutil/src/sptime/t_cgiparse.x
new file mode 100644
index 00000000..8985a8f0
--- /dev/null
+++ b/noao/obsutil/src/sptime/t_cgiparse.x
@@ -0,0 +1,110 @@
+define SZ_QUERY 2048
+
+
+# T_CGIPARSE -- Parse the CGI QUERY_STRING environment variable.
+# The string is expected to be a set of task.param=value pairs.
+# The task parameter values are set.
+
+procedure t_cgiparse ()
+
+char c, hex[2]
+int i
+long hexval
+pointer sp, str, par, val
+pointer ip, op
+
+int envfind(), gctol()
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_QUERY, TY_CHAR)
+ call salloc (par, SZ_LINE, TY_CHAR)
+ call salloc (val, SZ_LINE, TY_CHAR)
+
+ # Get the query string. If there isn't one then do nothing.
+ if (envfind ("QUERY_STRING", Memc[str], SZ_QUERY) <= 0)
+ return
+
+ # Parse the query string into parameters and values.
+ # For each pair set the task parameter.
+
+ op = par
+ for (ip=str;; ip=ip+1) {
+ c = Memc[ip]
+ switch (c) {
+ case '&', EOS: # End of parameter=value
+ Memc[op] = EOS
+ call cgi_setpar (Memc[par], Memc[val])
+ if (c == EOS)
+ break
+ Memc[par] = EOS
+ Memc[val] = EOS
+ op = par
+ next
+ case '=': # Separator between parameters and value
+ Memc[op] = EOS
+ op = val
+ next
+ case '+': # Space character
+ c = ' '
+ case '%': # Special characters in hex
+ call strcpy (Memc[ip+1], hex, 2)
+ i = 1
+ if (gctol (hex, i, hexval, 16) > 0) {
+ c = hexval
+ ip = ip + 2
+ }
+ }
+
+ Memc[op] = c
+ op = op + 1
+ }
+
+ call sfree (sp)
+end
+
+
+# CGI_SETPAR -- Set parameter value.
+# There is no error check for undefined tasks or parameters.
+# Values of the wrong type are ignored.
+
+procedure cgi_setpar (param, val)
+
+char param[ARB] #I Task parameter
+char val[ARB] #I Value string
+
+bool bval, streq()
+int ival, ip, ctoi(), ctod()
+double dval
+pointer sp, str, type
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+ call salloc (type, 10, TY_CHAR)
+
+ # Determine if parameter type.
+ call sprintf (Memc[str], SZ_LINE, "%s.p_type")
+ call pargstr (param)
+ call clgstr (Memc[str], Memc[type], 10)
+
+ # Set parameter in the approriate type.
+ ip = 1
+ switch (Memc[type]) {
+ case 'i':
+ if (ctoi (val, ip, ival) > 0)
+ call clputi (param, ival)
+ case 'r':
+ if (ctod (val, ip, dval) > 0)
+ call clputd (param, dval)
+ case 'b':
+ if (streq (val, "no") || streq (val, "yes")) {
+ bval = streq (val, "yes")
+ call clputb (param, bval)
+ }
+ default:
+ call clpstr (param, val)
+ }
+
+ call sfree (sp)
+end
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
diff --git a/noao/obsutil/src/sptime/tabinterp.x b/noao/obsutil/src/sptime/tabinterp.x
new file mode 100644
index 00000000..518f7b20
--- /dev/null
+++ b/noao/obsutil/src/sptime/tabinterp.x
@@ -0,0 +1,698 @@
+include <error.h>
+include <mach.h>
+include <math/iminterp.h>
+
+# Table structure.
+define TAB_DIM 2 # Maximum dimension of table
+define TAB_SZFNAME 99 # Size of file name
+define TAB_SZPARAMS (10*SZ_LINE) # Size of parameter string
+define TAB_SIZE (55+2*TAB_DIM) # Size of table structure
+
+define TAB_FNAME Memc[P2C($1)] # File name
+define TAB_VALUE Memr[P2R($1+50)] # Constant table value
+define TAB_PARAMS Memi[$1+51] # Pointer to parameters
+define TAB_NDIM Memi[$1+52] # Dimension of table
+define TAB_INTERP Memi[$1+53] # Interpolation pointer
+define TAB_NEXTRAP Memi[$1+54] # Number of extrapolations
+define TAB_LEN Memi[$1+54+$2] # Length of axes in table
+define TAB_COORD Memi[$1+54+$2+TAB_DIM] # Pointer to axes coordinates
+
+
+procedure tabtest ()
+
+char table[SZ_FNAME], param[SZ_FNAME], str[SZ_FNAME]
+int clglpr()
+real x, y, z, tabinterp1(), tabinterp2(), tabgetr()
+pointer tab, tabopen()
+errchk tabopen, tabinterp1, tabinterp2, tabgetr, tabgstr
+
+begin
+ tab = tabopen ()
+
+ call clgstr ("table", table, SZ_FNAME)
+
+ call clgstr ("param", param, SZ_FNAME)
+ z = tabgetr (tab, table, "", param)
+ call tabgstr (tab, table, "", param, str, SZ_FNAME)
+ call printf ("%s = %g = %s\n")
+ call pargstr (param)
+ call pargr (z)
+ call pargstr (str)
+
+ while (clglpr ("x", x) != EOF) {
+ iferr {
+ z = tabinterp1 (tab, table, x)
+ call printf ("%g %g\n")
+ call pargr (x)
+ call pargr (z)
+ } then {
+ if (clglpr ("y", y) == EOF)
+ break
+ z = tabinterp2 (tab, table, x, y)
+ call printf ("%g %g %g\n")
+ call pargr (x)
+ call pargr (y)
+ call pargr (z)
+ }
+ }
+
+ call tabclose (tab)
+end
+
+
+# TABOPEN -- Open table interpolation package.
+
+pointer procedure tabopen ()
+
+pointer stp, stopen()
+
+begin
+ stp = stopen ("tables", 10, 1000, 1000)
+ return (stp)
+end
+
+
+# TABCLOSE -- Close table interpolation package.
+
+procedure tabclose (stp)
+
+pointer stp #I Symbol table pointer
+
+pointer sym, sthead(), stnext()
+
+begin
+ for (sym = sthead(stp); sym != NULL; sym = stnext(stp, sym))
+ call tabfree (sym)
+ call stclose (stp)
+end
+
+
+procedure tabfree (sym)
+
+pointer sym #I Table structure
+
+int i
+
+begin
+ call mfree (TAB_PARAMS(sym), TY_CHAR)
+ if (TAB_INTERP(sym) != NULL) {
+ if (TAB_LEN(sym,1) > 1 && TAB_LEN(sym,2) > 1)
+ call msifree (TAB_INTERP(sym))
+ else
+ call asifree (TAB_INTERP(sym))
+ do i = 1, TAB_NDIM(sym)
+ call mfree (TAB_COORD(sym,i), TY_REAL)
+ }
+end
+
+
+
+# TABLOAD -- Load a table in the symbol table.
+
+procedure tabload (stp, name, fname)
+
+pointer stp # Symbol table pointer
+char name[ARB] # Name of table
+char fname[ARB] # File name of table
+
+int i, fd, ndim, npts, len[2], open(), errget(), nowhite(), ctor(), strlen()
+real value
+pointer sp, str, sym, params, data, coords[2], stfind(), stenter()
+bool streq()
+errchk open, tabread
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ # If no value is specified then don't enter into symbol table.
+ if (nowhite (fname, Memc[str], SZ_LINE) == 0) {
+ call sfree (sp)
+ return
+ }
+
+ # The special string "none" is equivalent to no table.
+ if (streq (Memc[str], "none")) {
+ call sfree (sp)
+ return
+ }
+
+ # Check if table has already been loaded.
+ sym = stfind (stp, name)
+ if (sym != NULL) {
+ if (streq (fname, TAB_FNAME(sym)))
+ return
+ }
+
+ # Check if constant is specified.
+ i = 1
+ if (ctor (Memc[str], i, value) == strlen (Memc[str])) {
+ sym = stenter (stp, name, TAB_SIZE)
+ call strcpy (fname, TAB_FNAME(sym), TAB_SZFNAME)
+ TAB_VALUE(sym) = value
+ TAB_NDIM(sym) = -1
+ TAB_PARAMS(sym) = NULL
+ TAB_INTERP(sym) = NULL
+
+ call sfree (sp)
+ return
+ }
+
+ # Read the table.
+ fd = open (Memc[str], READ_ONLY, TEXT_FILE)
+ iferr (call tabread (fd, params, data, npts, len, coords, ndim)) {
+ ndim = errget (Memc[str], SZ_LINE)
+ call strcat (" (", Memc[str], SZ_LINE)
+ call strcat (name, Memc[str], SZ_LINE)
+ call strcat (")", Memc[str], SZ_LINE)
+ call error (1, Memc[str])
+ }
+ call close (fd)
+
+ if (sym == NULL)
+ sym = stenter (stp, name, TAB_SIZE)
+ else
+ call tabfree (sym)
+
+ if (data != NULL) {
+ if (len[1] > 1 && len[2] > 1) {
+ call msiinit (TAB_INTERP(sym), II_BILINEAR)
+ call msifit (TAB_INTERP(sym), Memr[data], len[1], len[2],
+ len[1])
+ } else if (len[2] == 1) {
+ if (len[1] == 1)
+ call asiinit (TAB_INTERP(sym), II_NEAREST)
+ else
+ call asiinit (TAB_INTERP(sym), II_LINEAR)
+ call asifit (TAB_INTERP(sym), Memr[data], len[1])
+ } else {
+ call asiinit (TAB_INTERP(sym), II_LINEAR)
+ call asifit (TAB_INTERP(sym), Memr[data], len[2])
+ }
+ } else
+ TAB_INTERP(sym) = NULL
+
+ call strcpy (fname, TAB_FNAME(sym), TAB_SZFNAME)
+ TAB_PARAMS(sym) = params
+ TAB_NDIM(sym) = ndim
+ TAB_NEXTRAP(sym) = 0
+ call amovi (len, TAB_LEN(sym,1), 2)
+ call amovi (coords, TAB_COORD(sym,1), 2)
+ call mfree (data, TY_REAL)
+ call sfree (sp)
+end
+
+
+# TABREAD -- Read data from a table.
+
+procedure tabread (fd, params, data, npts, len, coords, ndim)
+
+int fd #I File descriptor
+pointer params #O Pointer to parameters
+pointer data #O Pointer to data
+int npts #O Number of data points
+int len[ARB] #O Number of points along each dimension
+pointer coords[ARB] #O Pointers to coordinates along each dimension
+int ndim #O Dimension of table
+
+int i, j, nread, fdparams, fscan(), nscan(), stropen(), ctor()
+real coord[4], scale
+pointer sp, cmt, key, eq, val
+bool streq()
+errchk stropen()
+
+begin
+ iferr {
+ call smark (sp)
+ call salloc (cmt, 2, TY_CHAR)
+ call salloc (key, SZ_FNAME, TY_CHAR)
+ call salloc (eq, 1, TY_CHAR)
+ call salloc (val, SZ_LINE, TY_CHAR)
+
+ npts = 0
+ ndim = 0
+ params = NULL
+ data = NULL
+ scale = INDEF
+ call aclri (len, 2)
+ call aclri (coords, 2)
+ while (fscan (fd) != EOF) {
+ call gargr (coord[1])
+ call gargr (coord[2])
+ call gargr (coord[3])
+ call gargr (coord[4])
+ nread = nscan()
+
+ if (nread == 0) { # Check for parameters.
+ call reset_scan()
+ call gargwrd (Memc[cmt], 2)
+ call gargwrd (Memc[key], SZ_LINE)
+ call gargwrd (Memc[eq], 1)
+ call gargwrd (Memc[val], SZ_LINE)
+ if (nscan() != 4 || Memc[cmt] != '#' || Memc[eq] != '=' ||
+ Memc[cmt+1] != EOS)
+ next
+ if (streq (Memc[key], "tablescale")) {
+ i = 1
+ if (ctor (Memc[val], i, scale) == 0)
+ call error (1, "Syntax error in table scale factor")
+ }
+ if (params == NULL) {
+ call malloc (params, TAB_SZPARAMS, TY_CHAR)
+ fdparams = stropen (Memc[params], TAB_SZPARAMS,
+ WRITE_ONLY)
+ }
+ call fprintf (fdparams, "%s %s\n")
+ call pargstr (Memc[key])
+ call pargstr (Memc[val])
+ next
+ } else if (nread == 1)
+ next
+
+ if (ndim == 0) {
+ ndim = nread - 1
+ if (ndim > 2)
+ call error (1, "Table dimension is too high")
+ do i = ndim+1, 2
+ len[i] = 1
+ }
+ if (nread-1 != ndim)
+ call error (2, "Table has variable number of columns")
+
+ do i = 1, ndim {
+ if (len[i] == 0) {
+ call malloc (coords[i], 100, TY_REAL)
+ Memr[coords[i]] = coord[i]
+ len[i] = len[i] + 1
+ } else {
+ do j = 0, len[i] - 1
+ if (coord[i] == Memr[coords[i]+j])
+ break
+ if (j >= len[i]) {
+ if (mod (len[i], 100) == 0)
+ call realloc (coords[i], len[i]+100, TY_REAL)
+ Memr[coords[i]+len[i]] = coord[i]
+ len[i] = len[i] + 1
+ }
+ }
+ }
+
+ if (npts == 0)
+ call malloc (data, 100, TY_REAL)
+ else if (mod (npts, 100) == 0)
+ call realloc (data, npts+100, TY_REAL)
+
+ Memr[data+npts] = coord[nread]
+ npts = npts + 1
+ }
+
+ if (npts > 0) {
+ j = 1
+ do i = 1, ndim
+ j = j * len[i]
+ if (j != npts)
+ call error (4, "Table is not regular")
+ }
+
+ if (!IS_INDEF(scale))
+ call amulkr (Memr[data], scale, Memr[data], npts)
+
+ call close (fdparams)
+ call sfree (sp)
+ } then {
+ if (params != NULL) {
+ call close (fdparams)
+ call mfree (params, TY_CHAR)
+ }
+ do i = 1, 2
+ call mfree (coords[i], TY_REAL)
+ call mfree (data, TY_REAL)
+ call sfree (sp)
+ call erract (EA_ERROR)
+ }
+end
+
+
+# TABEXISTS -- Determine if table exists.
+
+bool procedure tabexists (stp, name)
+
+pointer stp #I Symbol table pointer
+char name[ARB] #I Name of table
+
+pointer stfind()
+
+begin
+ return (stfind (stp, name) != NULL)
+end
+
+
+# TABGETR -- Get real parameter from table.
+
+real procedure tabgetr (stp, name, alt, param)
+
+pointer stp #I Symbol table pointer
+char name[ARB] #I Name of table
+char alt[ARB] #I Name of alternate table
+char param[ARB] #I Parameter name
+real val #O Returned value
+
+real rval
+int i, fd, strncmp(), strdic(), stropen(), fscan(), nscan()
+bool streq()
+pointer sp, key, sym[2], stfind()
+errchk stfind, stropen
+
+begin
+ # Return if no table.
+ if (name[1] == EOS && alt[1] == EOS)
+ return (INDEFR)
+
+ call smark (sp)
+ call salloc (key, SZ_FNAME, TY_CHAR)
+
+ # Find tables.
+ sym[1] = stfind (stp, name)
+ if (alt[1] != EOS)
+ sym[2] = stfind (stp, alt)
+ else
+ sym[2] = NULL
+
+ # Return an error if there is no table.
+ if (sym[1] == NULL && sym[2] == NULL) {
+ call sprintf (Memc[key], SZ_FNAME, "Table `%s' not found")
+ call pargstr (name)
+ call error (1, Memc[key])
+ }
+
+ # Get the parameter value.
+ val = INDEFR
+ do i = 1, 2 {
+ if (sym[i] == NULL)
+ next
+ if (strncmp (param, "table.", 6) == 0) {
+ switch (strdic (param[7], Memc[key], SZ_FNAME,
+ "|ndim|xmin|xmax|ymin|ymax|nextrap|")) {
+ case 1:
+ val = TAB_NDIM(sym[i])
+ case 2:
+ if (TAB_NDIM(sym[i]) >= 1)
+ val = Memr[TAB_COORD(sym[i],1)]
+ else if (TAB_NDIM(sym[i]) == -1)
+ val = -MAX_REAL
+ case 3:
+ if (TAB_NDIM(sym[i]) >= 1)
+ val = Memr[TAB_COORD(sym[i],1)+TAB_LEN(sym[i],1)-1]
+ else if (TAB_NDIM(sym[i]) == -1)
+ val = MAX_REAL
+ case 4:
+ if (TAB_NDIM(sym[i]) >= 2)
+ val = Memr[TAB_COORD(sym[i],2)]
+ else if (TAB_NDIM(sym[i]) == -1)
+ val = -MAX_REAL
+ case 5:
+ if (TAB_NDIM(sym[i]) >= 2)
+ val = Memr[TAB_COORD(sym[i],2)+TAB_LEN(sym[i],1)-1]
+ else if (TAB_NDIM(sym[i]) == -1)
+ val = MAX_REAL
+ case 6:
+ val = TAB_NEXTRAP(sym[i])
+ }
+ break
+
+ } else if (TAB_PARAMS(sym[i]) != NULL) {
+ fd = stropen (Memc[TAB_PARAMS(sym[i])], TAB_SZPARAMS, READ_ONLY)
+ while (fscan (fd) != EOF) {
+ call gargwrd (Memc[key], SZ_FNAME)
+ call gargr (rval)
+ if (nscan() != 2)
+ next
+ if (streq (Memc[key], param)) {
+ val = rval
+ break
+ }
+ }
+ call close (fd)
+ }
+
+ if (!IS_INDEF(val))
+ break
+ }
+
+ # Return error if no value was found.
+ if (IS_INDEF(val)) {
+ call sprintf (Memc[key], SZ_FNAME,
+ "Table parameter `%s' not found (%s)")
+ call pargstr (param)
+ call pargstr (name)
+ call error (1, Memc[key])
+ }
+
+ call sfree (sp)
+ return (val)
+
+end
+
+
+# TABGETI -- Get integer paraemter from table.
+
+int procedure tabgeti (stp, name, alt, param)
+
+pointer stp #I Symbol table pointer
+char name[ARB] #I Name of table
+char alt[ARB] #I Name of alternate table
+char param[ARB] #I Parameter name
+
+real val, tabgetr()
+errchk tabgetr
+
+begin
+ val = tabgetr (stp, name, alt, param)
+ if (IS_INDEFR(val))
+ return (INDEFI)
+
+ return (nint(val))
+end
+
+
+# TABGSTR -- Get string parameter from table.
+
+procedure tabgstr (stp, name, alt, param, val, maxchar)
+
+pointer stp #I Symbol table pointer
+char name[ARB] #I Name of table
+char alt[ARB] #I Name of alternate table
+char param[ARB] #I Parameter name
+char val[ARB] #O Returned value
+int maxchar #I Maximum string length
+
+int i, fd, strncmp(), strdic(), stropen(), fscan(), nscan()
+bool streq()
+pointer sp, key, str, sym[2], stfind()
+errchk stfind, stropen()
+
+begin
+ # Return if no table.
+ if (name[1] == EOS && alt[1] == EOS) {
+ val[1] = EOS
+ return
+ }
+
+ call smark (sp)
+ call salloc (key, SZ_FNAME, TY_CHAR)
+ call salloc (str, maxchar, TY_CHAR)
+
+ # Find tables.
+ sym[1] = stfind (stp, name)
+ if (alt[1] != EOS)
+ sym[2] = stfind (stp, alt)
+ else
+ sym[2] = NULL
+
+ # Return an error if there is no table.
+ if (sym[1] == NULL && sym[2] == NULL) {
+ call sprintf (Memc[key], SZ_FNAME, "Table `%s' not found")
+ call pargstr (name)
+ call error (1, Memc[key])
+ }
+
+ # Get the parameter value.
+ val[1] = EOS
+ do i = 1, 2 {
+ if (sym[i] == NULL)
+ next
+ if (strncmp (param, "table.", 6) == 0) {
+ switch (strdic (param[7], Memc[key], SZ_FNAME, "|filename|")) {
+ case 1:
+ call strcpy (TAB_FNAME(sym[i]), val, maxchar)
+ }
+ break
+
+ } else if (TAB_PARAMS(sym[i]) != NULL) {
+ fd = stropen (Memc[TAB_PARAMS(sym[i])], TAB_SZPARAMS, READ_ONLY)
+ while (fscan (fd) != EOF) {
+ call gargwrd (Memc[key], SZ_FNAME)
+ call gargstr (Memc[str], SZ_LINE)
+ if (nscan() != 2)
+ next
+ if (streq (Memc[key], param)) {
+ if (val[1] == EOS)
+ call strcpy (Memc[str+1], val, maxchar)
+ else {
+ call strcat ("\n", val, maxchar)
+ call strcat (Memc[str+1], val, maxchar)
+ }
+ }
+ }
+ call close (fd)
+
+ if (val[1] != EOS)
+ break
+ }
+ }
+
+ # Return error if no value was found.
+ if (val[1] == EOS) {
+ call sprintf (Memc[key], SZ_FNAME,
+ "Table parameter `%s' not found (%s)")
+ call pargstr (param)
+ call pargstr (name)
+ call error (1, Memc[key])
+ }
+
+ call sfree (sp)
+end
+
+
+# TABINTERP1 -- Interpolate a named 1D table.
+
+real procedure tabinterp1 (stp, name, x)
+
+pointer stp # Symbol table pointer
+char name[ARB] # Name of table
+real x # Interpolation coordinate
+
+char err[SZ_FNAME]
+int i, nx
+real xi, y, asieval()
+pointer sym, xp, stfind()
+errchk stfind
+
+begin
+ # Find the table.
+ sym = stfind (stp, name)
+ if (sym == NULL) {
+ call sprintf (err, SZ_FNAME, "Table `%s' not found")
+ call pargstr (name)
+ call error (1, err)
+ }
+
+ # If a constant table return the value.
+ if (TAB_NDIM(sym) == -1)
+ return (TAB_VALUE(sym))
+
+ # Check if the table is of the proper dimensionality.
+ if (TAB_NDIM(sym) != 1) {
+ call sprintf (err, SZ_FNAME, "Table is not one dimensional (%s)")
+ call pargstr (name)
+ call error (1, err)
+ }
+
+ nx = TAB_LEN(sym,1)
+ xp = TAB_COORD(sym,1)
+ if (x < Memr[xp] || x > Memr[xp+nx-1])
+ TAB_NEXTRAP(sym) = TAB_NEXTRAP(sym) + 1
+
+ if (nx == 1)
+ xi = 1
+ else {
+ do i = 1, nx-2
+ if (x < Memr[xp+i])
+ break
+
+ xi = i + (x - Memr[xp+i-1]) / (Memr[xp+i] - Memr[xp+i-1])
+ }
+ xi = max (1., min (real(nx), xi))
+ y = asieval (TAB_INTERP(sym), xi)
+
+
+ return (y)
+end
+
+
+# TABINTERP2 -- Interpolate a named 2D table.
+
+real procedure tabinterp2 (stp, name, x, y)
+
+pointer stp # Symbol table pointer
+char name[ARB] # Name of table
+real x, y # Interpolation coordinate
+
+char err[SZ_FNAME]
+int i, nx, ny
+real xi, yi, z, asieval(), msieval()
+pointer sym, xp, yp, stfind()
+errchk stfind
+
+begin
+ # Find the table.
+ sym = stfind (stp, name)
+ if (sym == NULL) {
+ call sprintf (err, SZ_FNAME, "Table `%s' not found")
+ call pargstr (name)
+ call error (1, err)
+ }
+
+ # If a constant table return the value.
+ if (TAB_NDIM(sym) == -1)
+ return (TAB_VALUE(sym))
+
+ # Check if the table is of the proper dimensionality.
+ if (TAB_NDIM(sym) != 2) {
+ call sprintf (err, SZ_FNAME, "Table is not two dimensional (%s)")
+ call pargstr (name)
+ call error (1, err)
+ }
+
+ nx = TAB_LEN(sym,1)
+ ny = TAB_LEN(sym,2)
+ if (nx > 1 && ny > 1) { # 2D interpolation
+ xp = TAB_COORD(sym,1)
+ do i = 1, nx-2
+ if (x < Memr[xp+i])
+ break
+ xi = i + (x - Memr[xp+i-1]) / (Memr[xp+i] - Memr[xp+i-1])
+ xi = max (1., min (real(nx), xi))
+ yp = TAB_COORD(sym,2)
+ do i = 1, ny-2
+ if (y < Memr[yp+i])
+ break
+ yi = i + (y - Memr[yp+i-1]) / (Memr[yp+i] - Memr[yp+i-1])
+ yi = max (1., min (real(nx), yi))
+ z = msieval (TAB_INTERP(sym), xi, yi)
+ } else if (ny == 1) { # 1D interpolation in x
+ if (nx == 1)
+ xi = 1
+ else {
+ xp = TAB_COORD(sym,1)
+ do i = 1, nx-2
+ if (x < Memr[xp+i])
+ break
+
+ xi = i + (x - Memr[xp+i-1]) / (Memr[xp+i] - Memr[xp+i-1])
+ }
+ xi = max (1., min (real(nx), xi))
+ z = asieval (TAB_INTERP(sym), xi)
+ } else { # 1D interpolation in y
+ yp = TAB_COORD(sym,2)
+ do i = 1, ny-2
+ if (y < Memr[yp+i])
+ break
+
+ yi = i + (y - Memr[yp+i-1]) / (Memr[yp+i] - Memr[yp+i-1])
+ yi = max (1., min (real(ny), yi))
+ z = asieval (TAB_INTERP(sym), yi)
+ }
+
+ return (z)
+end
diff --git a/noao/obsutil/src/sptime/x_spectime.x b/noao/obsutil/src/sptime/x_spectime.x
new file mode 100644
index 00000000..216aaa66
--- /dev/null
+++ b/noao/obsutil/src/sptime/x_spectime.x
@@ -0,0 +1,2 @@
+task sptime = t_sptime,
+ cgiparse = t_cgiparse
generated by cgit (git 2.34.1) at 2025-09-20 09:04:43 -0400