1 files changed, 70 insertions, 0 deletions
diff --git a/noao/astutil/asttools/astvorbit.x b/noao/astutil/asttools/astvorbit.x
new file mode 100644
index 00000000..61c5f0c5
--- /dev/null
+++ b/noao/astutil/asttools/astvorbit.x
@@ -0,0 +1,70 @@
+include	<math.h>
+
+# AST_VORBIT -- Radial velocity component of the Earth-Moon barycenter
+# relative to the Sun.
+
+procedure ast_vorbit (ra, dec, epoch, v)
+
+double	ra		# Right ascension of observation (hours)
+double	dec		# Declination of observation (degrees)
+double	epoch		# Julian epoch of observation
+double	v		# Component of orbital velocity (km/s)
+
+double	t, manom, lperi, oblq, eccen, tanom, slong, r, d, l, b, vorb
+double	ast_julday()
+
+begin
+	# T is the number of Julian centuries since J1900.
+	t = (ast_julday (epoch) - 2415020d0) / 36525.
+
+	# 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)
+
+	# VORB is the component of the Earth's orbital velocity perpendicular
+	# to the radius vector (km/s) where the Earth's semi-major axis is
+	# 149598500 km and the year is 365.2564 days.
+
+	vorb = ((TWOPI / 365.2564d0) *
+	    149598500.d0 / sqrt (1. - eccen**2)) / 86400.d0
+
+	# V is the projection onto the line of sight to the observation of
+	# the velocity of the Earth-Moon barycenter with respect to the
+	# Sun (km/s).
+
+	v = vorb * cos (b) * (sin (slong - l) - eccen * sin (lperi - l))
+end