1 files changed, 60 insertions, 0 deletions
diff --git a/noao/twodspec/longslit/airmass.x b/noao/twodspec/longslit/airmass.x
new file mode 100644
index 00000000..d47fab2d
--- /dev/null
+++ b/noao/twodspec/longslit/airmass.x
@@ -0,0 +1,60 @@
+include	<math.h>
+
+# IMG_AIRMASS -- Get or compute the image airmass from the image header.
+# If the airmass cannot be determined from header then INDEF is returned.
+#
+# Airmass formulation from Allen "Astrophysical Quantities" 1973 p.125,133
+# and John Ball's book on Algorithms for the HP-45.
+
+real procedure img_airmass (im)
+
+pointer	im			# IMIO pointer
+
+real	airmass, zd, ha, ra, dec, st, latitude, coszd, scale, x
+
+int	imaccf()
+real	imgetr()
+errchk	imgetr()
+
+data	scale/750.0/		# Atmospheric scale height approx
+
+begin
+	# If the airmass is in the header return its value.
+
+	if (imaccf (im, "airmass") == YES)
+	    return (imgetr (im, "airmass"))
+
+	# Compute zenith distance if not defined.
+
+	iferr (zd = imgetr (im, "zd")) {
+	    
+	    # Compute hour angle if not defined.
+
+	    iferr (ha = imgetr (im, "ha")) {
+		st = imgetr (im, "st")
+		ra = imgetr (im, "ra")
+		ha = st - ra
+		call imaddr (im, "ha", ha)
+	    }
+	
+	    dec = imgetr (im, "dec")
+	    latitude = imgetr (im, "latitude")
+
+	    ha = DEGTORAD (ha) * 15
+	    dec = DEGTORAD (dec)
+	    latitude = DEGTORAD (latitude)
+	    coszd = sin (latitude) * sin (dec) +
+		cos (latitude) * cos (dec) * cos (ha)
+	    zd = RADTODEG (acos (coszd))
+	    call imaddr (im, "zd", zd)
+	}
+
+	# Compute airmass from zenith distance.
+
+	zd = DEGTORAD (zd)
+	x = scale * cos (zd)
+	airmass = sqrt (x ** 2 + 2 * scale + 1) - x
+	call imaddr (im, "airmass", airmass)
+
+	return (airmass)
+end