1 files changed, 54 insertions, 0 deletions
diff --git a/noao/onedspec/getairm.x b/noao/onedspec/getairm.x
new file mode 100644
index 00000000..5226702d
--- /dev/null
+++ b/noao/onedspec/getairm.x
@@ -0,0 +1,54 @@
+# GET_AIRM -- Derive airmass value from telescope data if possible
+#             Otherwise return INDEF
+#
+# Note that HA must be negative to the East.
+# If HA is not reasonable, then ST-RA is used
+
+procedure get_airm (ra, dec, ha, st, latitude, airm)
+
+real	ra, dec, ha, st, latitude, airm
+
+begin
+	# Verify realistic value for HA
+
+	if (IS_INDEF (ha)) {
+	    if (IS_INDEF (st) || IS_INDEF (ra))
+		call error (0, "Can't determine airmass")
+
+	    ha = st - ra
+	}
+
+	# Now verify DEC
+	if (IS_INDEF (dec))
+	    call error (0, "Can't determine airmass")
+
+	# Everything should be just fine now
+	# Compute airmass using method of John Ball
+
+	call airmass (dec, ha, latitude, airm)
+end
+
+# AIRMASS -- Compute airmass from RA, DEC, and HA
+#
+# Airmass formulation from Allen "Astrophysical Quantities" 1973 p.125,133.
+# and John Ball's book on Algorithms for the HP-45
+
+procedure airmass (dec, ha, latitude, airm)
+
+real	dec, ha, latitude, airm
+
+real	scale, rads, cos_zd, sin_elev
+real	x
+
+data	rads /57.29577951/		# Degrees per radian
+data	scale/750.0      /		# Atmospheric scale height approx
+
+begin
+	cos_zd = sin (latitude/rads) * sin (dec/rads) +
+	    cos (latitude/rads) * cos (dec/rads) * cos (ha*15/rads)
+
+	sin_elev = cos_zd		# SIN of elev = cos of Zenith dist
+
+	x  = scale * sin_elev
+	airm = sqrt (x**2 + 2*scale + 1) - x
+end