1 files changed, 82 insertions, 0 deletions
diff --git a/noao/astutil/asttools/asthjd.x b/noao/astutil/asttools/asthjd.x
new file mode 100644
index 00000000..29efd0b5
--- /dev/null
+++ b/noao/astutil/asttools/asthjd.x
@@ -0,0 +1,82 @@
+include	<math.h>
+
+# AST_HJD -- Helocentric Julian Day from Epoch
+
+procedure ast_hjd (ra, dec, epoch, lt, hjd)
+
+double	ra		# Right ascension of observation (hours)
+double	dec		# Declination of observation (degrees)
+double	epoch		# Julian epoch of observation
+double	lt		# Light travel time in seconds
+double	hjd		# Helocentric Julian Day
+
+double	ast_julday()
+
+begin
+	call ast_jd_to_hjd (ra, dec, ast_julday(epoch), lt, hjd)
+end
+
+
+# AST_JD_TO_HJD -- Helocentric Julian Day from UT Julian date
+
+procedure ast_jd_to_hjd (ra, dec, jd, lt, hjd)
+
+double	ra		# Right ascension of observation (hours)
+double	dec		# Declination of observation (degrees)
+double	jd		# Geocentric Julian date of observation
+double	lt		# Light travel time in seconds
+double	hjd		# Helocentric Julian Day
+
+double	t, manom, lperi, oblq, eccen, tanom, slong, r, d, l, b, rsun
+
+begin
+	# JD is the geocentric Julian date.
+	# T is the number of Julian centuries since J1900.
+
+	t = (jd - 2415020d0) / 36525d0
+
+	# MANOM is the mean anomaly of the Earth's orbit (degrees)
+	# LPERI is the mean longitude of perihelion (degrees)
+	# OBLQ is the mean obliquity of the ecliptic (degrees)
+	# ECCEN is the eccentricity of the Earth's orbit (dimensionless)
+
+	manom = 358.47583d0 +
+	    t * (35999.04975d0 - t * (0.000150d0 + t * 0.000003d0))
+	lperi = 101.22083d0 +
+	    t * (1.7191733d0 + t * (0.000453d0 + t * 0.000003d0))
+	oblq = 23.452294d0 -
+	    t * (0.0130125d0 + t * (0.00000164d0 - t * 0.000000503d0))
+	eccen = 0.01675104d0 - t * (0.00004180d0 + t * 0.000000126d0)
+
+	# Convert to principle angles
+	manom = mod (manom, 360.0D0)
+	lperi = mod (lperi, 360.0D0)
+
+	# Convert to radians
+	r = DEGTORAD(ra * 15)
+	d = DEGTORAD(dec)
+	manom = DEGTORAD(manom)
+	lperi = DEGTORAD(lperi)
+	oblq = DEGTORAD(oblq)
+
+	# TANOM is the true anomaly (approximate formula) (radians)
+	tanom = manom + (2 * eccen - 0.25 * eccen**3) * sin (manom) +
+	    1.25 * eccen**2 * sin (2 * manom) +
+	    13./12. * eccen**3 * sin (3 * manom)
+
+	# SLONG is the true longitude of the Sun seen from the Earth (radians)
+	slong = lperi + tanom + PI
+
+	# L and B are the longitude and latitude of the star in the orbital
+	# plane of the Earth (radians)
+
+	call ast_coord (double (0.), double (0.), double (-HALFPI),
+	    HALFPI - oblq, r, d, l, b)
+
+	# R is the distance to the Sun.
+	rsun = (1. - eccen**2) / (1. + eccen * cos (tanom))
+
+	# LT is the light travel difference to the Sun.
+	lt = -0.005770d0 * rsun * cos (b) * cos (l - slong)
+	hjd = jd + lt
+end