Initial commit

1 files changed, 1107 insertions, 0 deletions

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