Initial commit

20 files changed, 1344 insertions, 0 deletions

diff --git a/noao/astutil/asttools/PRECESS b/noao/astutil/asttools/PRECESS
new file mode 100644
index 00000000..48ca311a
--- /dev/null
+++ b/noao/astutil/asttools/PRECESS
@@ -0,0 +1,33 @@
+INPUT[1]           OUTPUT[2]		DEL[3]	     DEL[4]       DEL[5]
+ 1950.0              2000.0               SEC          SEC         SEC
+--------   -------------------------  ----------  ----------  ------------
+RA  DEC             RA           DEC    RA   DEC    RA   DEC     RA    DEC
+ 0    0     0:02:33.73    0:16:42.24  0.04  0.24  0.04  0.24   1.09   7.04
+ 6    0     6:02:33.72   -0:00:05.60  0.03  0.00  0.03  0.00  -0.45   0.89
+12    0    12:02:33.72   -0:16:42.24  0.03 -0.24  0.03 -0.24   0.71  -4.64
+18    0    18:02:33.72    0:00:05.60  0.03  0.00  0.03  0.00   2.25   1.51
+ 0   30     0:02:33.94   30:16:42.24  0.04  0.24  0.04  0.24   1.13  -3.24
+ 6   30     6:03:12.30   29:59:52.99  0.05 -0.01  0.05 -0.01  -0.43  -2.06
+12    0    12:02:33.72   -0:16:42.24  0.03 -0.24  0.03 -0.24   0.71  -4.64
+18    0    18:02:33.72    0:00:05.60  0.03  0.00  0.03  0.00   2.25   1.51
+ 0   60     0:02:34.38   60:16:42.24  0.04  0.24  0.04  0.24   1.32 -11.1
+ 6   60     6:04:29.44   59:59:50.18  0.06 -0.01  0.06 -0.01  -1.14  -6.33
+12   60    12:02:33.08   59:43:17.76  0.03 -0.24  0.03 -0.24   0.50  12.29
+18   60    18:00:37.99   60:00:01.38  0.01  0.00  0.01  0.00   2.92  -0.91
+ 0   90    12:01:16.87   89:43:17.75  0.02 -0.23 43200 -2004 -512529 -2019
+ 0  -30     0:02:33.51  -29:43:17.76  0.04  0.24  0.04  0.24   1.10  17.01
+ 6  -30     6:01:55.15  -30:00:04.20  0.03 -0.01  0.03 -0.01  -0.89   3.52
+12  -30    12:02:33.94  -30:16:42.24  0.04 -0.24  0.04 -0.24   0.70 -14.92
+18  -30    18:03:12.30  -29:59:52.99  0.05  0.01  0.05  0.01   2.68   1.37
+ 0  -60     0:02:33.08  -59:43:17.76  0.03  0.24  0.03  0.24   1.24  23.98
+ 6  -60     6:00:37.99  -60:00:01.38  0.01  0.00  0.01  0.00  -2.48   6.92
+12  -60    12:02:34.38  -60:16:42.24  0.04 -0.24  0.04 -0.24   0.58 -22.77
+18  -60    18:04:29.44  -59:59:50.18  0.06  0.01  0.06  0.01   4.26   2.73
+ 0  -90     0:01:16.87  -89:43:17.74  0.02  0.24  0.02  0.25 -733455 26.11
+
+ [1] Epoch 1950.0
+ [2] Epoch 2000.0 using ASTPRECESS.X based on 1984 Ephemeris
+ [3] Difference using PRECESSGJ.X (first IRAF procedure by G. Jacoby)
+ [4] Difference using routine originally for the Cyber by D. Wells
+ [5] Difference using PRECESSMGB.X (from DOPSET by Manchester, Gorder, and Ball)
+     It includes aberation and nutation.
diff --git a/noao/astutil/asttools/README b/noao/astutil/asttools/README
new file mode 100644
index 00000000..dbdb7d67
--- /dev/null
+++ b/noao/astutil/asttools/README
@@ -0,0 +1,137 @@
+Astronomical Tools:
+
+Precession:  These routines probably have some history even before the
+    authors quoted as original sources.  They have been collected and
+    some rewritten into SPP by F. Valdes (NOAO), March 1986.  PRECES.F
+    was used originally in V2.3-V2.5 IRAF by George Jacoby.  It has
+    been replaced by ASTPRECESS.X which is the only procedure in the
+    library.
+
+    astprecess.x -- Precession written by F. Valdes based on Astronomical
+	Almanac using new IAU system.
+    precessmgb.x -- Precession + aberration + nutation based on the work of
+	Manchester, Gordon, and Ball
+    precessgj.x -- Originally written by G. Jacoby in Fortran and distributed
+	with V2.3-V2.5 IRAF.  Transcribed to SPP by F. Valdes
+
+    Notes:
+	1. The differences between ASTPRECESS.X and PRECESSGJ.X (
+	   and a routine written by D. Wells for the Cyber and later
+	   used in other NOAO software) are on the order of a few
+	   tenths of a second of arc.  I believe the differences are
+	   due to using 1984 almanac methods in the former and much
+	   earlier methods for the latter.
+	3. PRECESSMGB.X differs considerably from the others to the order
+	   of many seconds of arc.  It does include other effects not
+	   present in the other routines.  It totally fails at DEC=+-90.
+	   It is based on roughly 1970 almanac methods.
+	4. See PRECESS.DOC for comparison.
+
+Radial Velocity:  These formulas for these routines were partly obtained
+    by inspection of the code for the subroutine DOP in the program DOPSET
+    written by R. N. Manchester and M. A. Gordon of NRAO dated January 1970.
+    They have been restructured, revised, and coded in SPP by F. Valdes.
+
+    astvr.x -- Project a velocity vector in radial velocity along line of sight.
+    astvbary.x -- Radial velocity component of center of the Earth relative to
+	to the barycenter of the Earth-Moon system.
+    astvrotate.x -- Radial velocity component of the observer relative to
+	the center of the Earth due to the Earth's rotation.
+    astvorbit.x -- Radial velocity component of the observer relative to
+	the center of the Earth due to the Earth's rotation.
+    astvsun.x -- Projection of the sun's velocity along the given direction.
+
+Coordinates:
+
+    astarcsep.x -- Arc distance (arcsec) between two spherical coordinates
+	(hours, degrees).
+    astcoord.x -- This procedure converts the longitude-latitude coordinates
+	(a1, b1) of a point on a sphere into corresponding coordinates
+	(a2, b2) in a different coordinate system that is specified by
+	the coordinates of its origin (ao, bo) and its north pole (ap,
+	bp) in the original coordinate system.  The range of a2 will be
+	from -pi to pi.
+    astgalactic.x -- Convert equatorial coordinates (1950) to galactic
+	coordinates.
+    astgaltoeq.x -- Convert galactic coordinates to equitorial (1950).
+
+Dates and times:
+
+    asttimes.x:
+	AST_DATE_TO_EPOCH -- Convert Gregorian date and solar mean time to
+	    a Julian epoch.  A Julian epoch has 365.25 days per year and 24
+	    hours per day.
+	AST_EPOCH_TO_DATE -- Convert a Julian epoch to year, month, day, and
+	    time.
+	AST_DAY_OF_YEAR -- The day number for the given year is returned.
+	AST_DAY_OF_WEEK -- Return the day of the week for the given Julian day.
+	    The integer day of the week is 0=Sunday - 6=Saturday.  The
+	    character string is the three character abbreviation for the day
+	    of the week.  Note that the day of the week is for Greenwich
+	    if the standard UT is used.
+	AST_JULDAY -- Convert epoch to Julian day.
+	AST_DATE_TO_JULDAY -- Convert date to Julian day.
+	AST_JULDAY_TO_DATE -- Convert Julian day to date.
+	AST_MST -- Mean sidereal time of the epoch at the given longitude.
+	    This procedure may be used to optain Greenwich Mean Sidereal Time
+	    (GMST) by setting the longitude to 0.
+    asthjd.x:
+	AST_HJD -- Helocentic Julian Day for a direction of observation.
+	AST_JD_TO_HJD -- Helocentic Julian Day for a direction of observation.
+
+Helocentric Parameters:
+
+    astdsun.x:
+	AST_DSUN - Distance to Sun in AU.
+
+Misc:
+
+    astlvac.x:
+	AST_LVAC - Convert air wavelengths to vacuum wavelengths (Angstroms)
+
+Y2K:
+    Most routines work in Julian days or epochs.  If they have an input
+    year it is converted to one of these forms by calling
+    ast_date_to_julday.  This is the only routine that has a Y2K
+    connection.  It assumes two digit years are 20th century.  These
+    routines are Y2K correct.
+
+
+The following are the comments and references from the DOPSET program
+noted above.  The HJD routine was also derived from this code.
+
+C     MODIFIED FOR IBM 360 BY R.N.MANCHESTER AND M.A.GORDON
+C     JANUARY 1970
+C
+C
+C  DOP CALCULATES THE VELOCITY COMPONENT OF THE OBSERVER WITH RESPECT
+C  TO THE LOCAL STANDARD OF REST AS PROJECTED ONTO A LINE SPECIFIED BY T
+C  ASCENSION AND DECLINATION (RAHRS, RAMIN, RASEC, DDEG, DMIN, DSEC) EPO
+C  DATE, FOR A TIME SPECIFIED AS FOLLOWS:  NYR = LAST TWO DIGITS OF THE
+C  (FOR 19XX A.D.), NDAY = DAY NUMBER (GMT), NHUT, NMUT, NSUT = HRS, MIN
+C  (GMT).  THE LOCATION OF THE OBSERVER IS SPECIFIED BY THE LATITUDE (AL
+C  LONGITUDE (OLONG) (GEODETIC) (IN DEGREES) AND ELEVATION (ELEV) (IN ME
+C  ABOVE MEAN SEA LEVEL.  THE SUBROUTINE OUTPUTS THE LOCAL MEAN SIDEREAL
+C  (XLST IN DAYS), THE COMPONENT OF THE SUN*S MOTION WITH RESPECT TO THE
+C  STANDARD OF REST AS PROJECTED ONTO THE LINE OF SIGHT TO THE SOURCE (V
+C  KM/SEC) AS WELL AS THE TOTAL VELOCITY COMPONENT V1 (KM/SEC).  POSITIV
+C  VELOCITY CORRESPONDS TO INCREASING DISTANCE BETWEEN SOURCE AND OBSERV
+C
+C  THIS VERSION OF DOP TAKES INTO ACCOUNT COMPONENTS OF THE OBSERVER*S
+C  MOTION DUE TO THE ROTATION OF THE EARTH, THE REVOLUTION OF THE EARTH-
+C  BARYCENTER ABOUT THE SUN, AND THE MOTION OF THE EARTH*S CENTER ABOUT
+C  EARTH-MOON BARYCENTER.  THE PERTURBATIONS OF THE EARTH*S ORBIT DUE TO
+C  PLANETS ARE NEGLECTED.  THE ABSOLUTE PRECISION OF THIS VERSION OF DOP
+C  ABOUT 0.004 KM/SEC, BUT SINCE THE DOMINANT ERROR TERM IS SLOWLY VARYI
+C  RELATIVE ERROR WILL BE CONSIDERABLY LESS FOR TIMES UP TO A WEEK OR SO
+C
+C  REFERENCES:  MCRAE, D. A., WESTERHOUT, G., TABLE FOR THE REDUCTION OF
+C                  VELOCITIES TO THE LOCAL STANDARD OF REST, THE OBSERVAY=2.
+C                  LUND, SWEDEN, 1956.
+C               SMART, W. M., TEXT-BOOK ON SPHERICAL ASTRONOMY, CAMBRIDG
+C                  UNIV. PRESS, 1962.
+C               THE AMERICAN EPHEMERIS AND NAUTICAL ALMANAC
+C               THE SUPPLEMENT TO THE ABOVE
+C
+C     VERSION OF JUNE 1969
+
diff --git a/noao/astutil/asttools/astarcsep.x b/noao/astutil/asttools/astarcsep.x
new file mode 100644
index 00000000..023fd20f
--- /dev/null
+++ b/noao/astutil/asttools/astarcsep.x
@@ -0,0 +1,33 @@
+include	<math.h>
+
+# AST_ARCSEP -- Arc distance (arcsec) between two spherical coordinates
+# (hours,deg).
+
+double procedure ast_arcsep (ra1, dec1, ra2, dec2) 
+
+double	ra1		# Right ascension (hours)
+double	dec1		# Declination (degrees)
+double	ra2		# Right ascension (hours)
+double	dec2		# Declination (degrees)
+
+double	a1, b1, c1, a2, b2, c2
+
+begin
+	# Direction cosines
+	a1 = cos (DEGTORAD (15D0 * ra1)) * cos (DEGTORAD (dec1))
+	b1 = sin (DEGTORAD (15D0 * ra1)) * cos (DEGTORAD (dec1))
+	c1 = sin (DEGTORAD (dec1))
+
+	a2 = cos (DEGTORAD (15D0 * ra2)) * cos (DEGTORAD (dec2))
+	b2 = sin (DEGTORAD (15D0 * ra2)) * cos (DEGTORAD (dec2))
+	c2 = sin (DEGTORAD (dec2))
+
+	# Modulus squared of half the difference vector.
+	a1 = ((a1-a2)**2 + (b1-b2)**2 + (c1-c2)**2) / 4
+
+	# Angle
+	a1 = 2D0 * atan2 (sqrt (a1), sqrt (max (0D0, 1D0-a1)))
+	a1 = RADTODEG (a1) * 3600D0
+
+	return (a1)
+end
diff --git a/noao/astutil/asttools/astcoord.x b/noao/astutil/asttools/astcoord.x
new file mode 100644
index 00000000..07a01ee3
--- /dev/null
+++ b/noao/astutil/asttools/astcoord.x
@@ -0,0 +1,56 @@
+# AST_COORD -- Convert spherical coordinates to new system.
+#
+# This procedure converts the longitude-latitude coordinates (a1, b1)
+# of a point on a sphere into corresponding coordinates (a2, b2) in a
+# different coordinate system that is specified by the coordinates of its
+# origin (ao, bo).  The range of a2 will be from -pi to pi.
+
+procedure ast_coord (ao, bo, ap, bp, a1, b1, a2, b2)
+
+double	ao, bo		# Origin of new coordinates (radians)
+double	ap, bp		# Pole of new coordinates (radians)
+double	a1, b1		# Coordinates to be converted (radians)
+double	a2, b2		# Converted coordinates (radians)
+
+double	sao, cao, sbo, cbo, sbp, cbp
+double	x, y, z, xp, yp, zp, temp
+
+begin
+	x = cos (a1) * cos (b1)
+	y = sin (a1) * cos (b1)
+	z = sin (b1)
+	xp = cos (ap) * cos (bp)
+	yp = sin (ap) * cos (bp)
+	zp = sin (bp)
+
+	# Rotate the origin about z.
+	sao = sin (ao)
+	cao = cos (ao)
+	sbo = sin (bo)
+	cbo = cos (bo)
+	temp = -xp * sao + yp * cao
+	xp = xp * cao + yp * sao
+	yp = temp
+	temp = -x * sao + y * cao
+	x = x * cao + y * sao
+	y = temp
+
+	# Rotate the origin about y.
+	temp = -xp * sbo + zp * cbo
+	xp = xp * cbo + zp * sbo
+	zp = temp
+	temp = -x * sbo + z * cbo
+	x = x * cbo + z * sbo
+	z = temp
+
+	# Rotate pole around x.
+	sbp = zp
+	cbp = yp
+	temp = y * cbp + z * sbp
+	y = y * sbp - z * cbp
+	z = temp
+
+	# Final angular coordinates.
+	a2 = atan2 (y, x)
+	b2 = asin (z)
+end
diff --git a/noao/astutil/asttools/astdsun.x b/noao/astutil/asttools/astdsun.x
new file mode 100644
index 00000000..89af95f7
--- /dev/null
+++ b/noao/astutil/asttools/astdsun.x
@@ -0,0 +1,19 @@
+include	<math.h>
+
+# AST_DSUN -- Distance to Sun in AU
+# Taken from Astronomical Almanac 1984, page C24.
+
+double procedure ast_dsun (epoch)
+
+double	epoch			# Epoch desired
+
+double	n, g, r
+double	ast_julday()
+
+begin
+	n = ast_julday (epoch) - 2451545d0
+	g = DEGTORAD (357.528d0 + 0.9856003d0 * n)
+	r = 1.00014d0 - 0.01671d0 * cos (g) - 0.00014d0 * cos (2 * g)
+
+	return (r)
+end
diff --git a/noao/astutil/asttools/astgalactic.x b/noao/astutil/asttools/astgalactic.x
new file mode 100644
index 00000000..ab51a918
--- /dev/null
+++ b/noao/astutil/asttools/astgalactic.x
@@ -0,0 +1,52 @@
+include <mach.h>
+include	<math.h>
+
+# Definition of system
+define	LONGNCP		123.00d0		# Longtitude of NCP
+define	RAGPOLE		192.25d0		# RA Galactic Pole (12:49)
+define	DECGPOLE	27.4d0			# DEC Galactic Pole (27:24)
+define	GEPOCH		1950.0d0		# Epoch of definition
+
+# AST_GALACTIC -- Convert equatorial coordinates to galactic coordinates.
+# Input coordinates and epoch are changed to epoch of galactic coordinates.
+
+procedure ast_galactic (ra, dec, epoch, lii, bii)
+
+double	ra		# Right ascension (hours)
+double	dec		# Declination (degrees)
+double	epoch		# Epoch of coordinates
+double	lii		# Galactic longitude (degrees)
+double	bii		# Galactic latitude (degrees)
+
+double	rar, decr, drar, cosdecg, sindecg, cosdecr, x, y, z, r, temp
+
+begin
+	# Precess the coordinates to 1950.0
+	call ast_precess (ra, dec, epoch, rar, decr, double (GEPOCH))
+
+	# Precompute the necessary constants.
+	drar = DEGTORAD (15.0d0 * rar - RAGPOLE) 
+	cosdecg = cos (DEGTORAD (DECGPOLE))
+	sindecg = sin (DEGTORAD(DECGPOLE))
+	cosdecr = cos (DEGTORAD (decr))
+
+	# Compute the tansformation equations
+	x = cosdecr * cos (drar)
+	y =  cosdecr * sin (drar)
+	z = sin (DEGTORAD (decr))
+	temp = z * cosdecg - x * sindecg 
+	z = z * sindecg + x * cosdecg 
+	x = temp
+	r = sqrt (x * x + y * y)
+
+	# Compute lii and bii and convert to degrees.
+	if (r < EPSILOND)
+	    lii = 0.0d0
+	else
+	    lii = DEGTORAD (LONGNCP)  + atan2 (-y, x)
+	if (lii < 0.0d0)
+	    lii = lii + TWOPI
+	bii = atan2 (z, r)
+	lii = RADTODEG (lii)
+	bii = RADTODEG (bii)
+end
diff --git a/noao/astutil/asttools/astgaltoeq.x b/noao/astutil/asttools/astgaltoeq.x
new file mode 100644
index 00000000..94f7861d
--- /dev/null
+++ b/noao/astutil/asttools/astgaltoeq.x
@@ -0,0 +1,37 @@
+include	<math.h>
+
+# Definition of system
+define	LP		123.00d0		# Longtitude of pole
+define	BP		27.40d0			# Latitude of pole
+define	LO		97.7422d0		# Longitude of origin
+define	BO		-60.1810d0		# Latitude of origin
+define	GEPOCH		1950.0d0		# Epoch of definition
+
+# AST_GALTOEQ -- Convert galactic coordinates (1950) to equatorial coordinates.
+
+procedure ast_galtoeq (lii, bii, ra, dec, epoch)
+
+double	lii		# Galactic longitude (degrees)
+double	bii		# Galactic latitude (degrees)
+double	ra		# Right ascension (hours)
+double	dec		# Declination (degrees)
+double	epoch		# Epoch of coordinates
+
+double	ao, bo, ap, bp, a1, b1, a2, b2
+
+begin
+	ao = DEGTORAD (LO)
+	bo = DEGTORAD (BO)
+	ap = DEGTORAD (LP)
+	bp = DEGTORAD (BP)
+	a1 = DEGTORAD (lii)
+	b1 = DEGTORAD (bii)
+
+	call ast_coord (ao, bo, ap, bp, a1, b1, a2, b2)
+
+	a2 = mod (24.0d0 + RADTODEG(a2) / 15.0d0, 24.0d0)
+	b2 = RADTODEG (b2)
+
+	# Precess the coordinates
+	call ast_precess (a2, b2, GEPOCH, ra, dec, epoch)
+end
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
diff --git a/noao/astutil/asttools/astlvac.x b/noao/astutil/asttools/astlvac.x
new file mode 100644
index 00000000..3ef4d895
--- /dev/null
+++ b/noao/astutil/asttools/astlvac.x
@@ -0,0 +1,29 @@
+# AST_LVAC -- Convert air wavelength to vacuum wavelength.
+#
+#     Convert LAMBDA(air) to LAMBDA(vacuum) using the formulae
+#
+#     (n-1)*1e8 = 8342.13 + 2406030/(130-s**2) + 15997/(38.9-s**2)
+#
+#                 where  s = 1/lambda(air)   (for lambda in MICRONS)
+#
+#     lambda(vac) = n * lambda(air)
+#
+#     NOTE: lambda is given in ANGSTROMS for this program.
+#
+
+procedure ast_lvac (lair, lvac, npts)
+
+double	lair[npts]		#I Air wavelength (Angstroms)
+double	lvac[npts]		#O Vacuum wavelength (Angstroms)
+int	npts			#I Number of points
+
+int	i
+double	s2, n
+
+begin
+	do i = 1, npts {
+	    s2 = 1D8 * (1. / lair[i]) ** 2
+	    n = 1 + ((8342.13 + 2406030 / (130-s2) + 15997 / (38.9-s2)) / 1D8)
+	    lvac[i] = n * lair[i]
+	}
+end
diff --git a/noao/astutil/asttools/astprecess.x b/noao/astutil/asttools/astprecess.x
new file mode 100644
index 00000000..6c69d792
--- /dev/null
+++ b/noao/astutil/asttools/astprecess.x
@@ -0,0 +1,117 @@
+include	<math.h>
+
+# AST_PRECESS -- Precess coordinates from epoch1 to epoch2.
+#
+# The method used here is based on the new IAU system described in the
+# supplement to the 1984 Astronomical Almanac.  The precession is
+# done in two steps; precess epoch1 to the standard epoch J2000.0 and then
+# precess from the standard epoch to epoch2.  The precession between
+# any two dates is done this way because the rotation matrix coefficients
+# are given relative to the standard epoch.
+
+procedure ast_precess (ra1, dec1, epoch1, ra2, dec2, epoch2)
+
+double	ra1, dec1, epoch1		# First coordinates
+double	ra2, dec2, epoch2		# Second coordinates
+
+double	r0[3], r1[3], p[3, 3]
+bool	fp_equald()
+
+begin
+	# If the input epoch is 0 or undefined then assume the input epoch
+	# is the same as the output epoch.  If the two epochs are the same
+	# then return the coordinates from epoch1.
+
+	if ((epoch1 == 0.) || IS_INDEFD (epoch1) || fp_equald(epoch1, epoch2)) {
+	    ra2 = ra1
+	    dec2 = dec1
+	    return
+	}
+
+	# Rectangular equitorial coordinates (direction cosines).
+	ra2 = DEGTORAD (ra1 * 15.)
+	dec2 = DEGTORAD (dec1)
+
+	r0[1] = cos (ra2) * cos (dec2)
+	r0[2] = sin (ra2) * cos (dec2)
+	r0[3] = sin (dec2)
+
+	# If epoch1 is not the standard epoch then precess to the standard
+	# epoch.
+
+	if (epoch1 != 2000.) {
+	    call ast_rotmatrix (epoch1, p)
+
+	    # Note that we multiply by the inverse of p which is the
+	    # transpose of p.
+
+	    r1[1] = p[1, 1] * r0[1] + p[1, 2] * r0[2] + p[1, 3] * r0[3]
+	    r1[2] = p[2, 1] * r0[1] + p[2, 2] * r0[2] + p[2, 3] * r0[3]
+	    r1[3] = p[3, 1] * r0[1] + p[3, 2] * r0[2] + p[3, 3] * r0[3]
+	    r0[1] = r1[1]
+	    r0[2] = r1[2]
+	    r0[3] = r1[3]
+	}
+
+	# If epoch2 is not the standard epoch then precess from the standard
+	# epoch to the desired epoch.
+
+	if (epoch2 != 2000.) {
+	    call ast_rotmatrix (epoch2, p)
+	    r1[1] = p[1, 1] * r0[1] + p[2, 1] * r0[2] + p[3, 1] * r0[3]
+	    r1[2] = p[1, 2] * r0[1] + p[2, 2] * r0[2] + p[3, 2] * r0[3]
+	    r1[3] = p[1, 3] * r0[1] + p[2, 3] * r0[2] + p[3, 3] * r0[3]
+	    r0[1] = r1[1]
+	    r0[2] = r1[2]
+	    r0[3] = r1[3]
+	}
+
+	# Convert from radians to hours and degrees.
+	ra2 = RADTODEG (atan2 (r0[2], r0[1]) / 15.)
+	dec2 = RADTODEG (asin (r0[3]))
+	if (ra2 < 0.)
+	    ra2 = ra2 + 24
+end
+
+
+# ROTMATRIX -- Compute the precession rotation matrix from the standard epoch
+# J2000.0 to the specified epoch.
+
+procedure ast_rotmatrix (epoch, p)
+
+double	epoch		# Epoch of date
+double	p[3, 3]		# Rotation matrix
+
+double	t, a, b, c, ca, cb, cc, sa, sb, sc
+double	ast_julday()
+
+begin
+	# The rotation matrix coefficients are polynomials in time measured
+	# in Julian centuries from the standard epoch.  The coefficients are
+	# in degrees.
+
+	t = (ast_julday (epoch) - 2451545.0d0) / 36525d0
+
+	a = t * (0.6406161d0 + t * (0.0000839d0 + t * 0.0000050d0))
+	b = t * (0.6406161d0 + t * (0.0003041d0 + t * 0.0000051d0))
+	c = t * (0.5567530d0 - t * (0.0001185d0 + t * 0.0000116d0))
+
+	# Compute the cosines and sines once for efficiency.
+	ca = cos (DEGTORAD (a))
+	sa = sin (DEGTORAD (a))
+	cb = cos (DEGTORAD (b))
+	sb = sin (DEGTORAD (b))
+	cc = cos (DEGTORAD (c))
+	sc = sin (DEGTORAD (c))
+
+	# Compute the rotation matrix from the sines and cosines.
+	p[1, 1] = ca * cb * cc - sa * sb
+	p[2, 1] = -sa * cb * cc - ca * sb
+	p[3, 1] = -cb * sc
+	p[1, 2] = ca * sb * cc + sa * cb
+	p[2, 2] = -sa * sb * cc + ca * cb
+	p[3, 2] = -sb * sc
+	p[1, 3] = ca * sc
+	p[2, 3] = -sa * sc
+	p[3, 3] = cc
+end
diff --git a/noao/astutil/asttools/asttimes.x b/noao/astutil/asttools/asttimes.x
new file mode 100644
index 00000000..f1a3674a
--- /dev/null
+++ b/noao/astutil/asttools/asttimes.x
@@ -0,0 +1,217 @@
+define	J2000		2000.0D0		# J2000
+define	JD2000		2451545.0D0		# J2000 Julian Date
+define	JYEAR		365.25D0		# Julian year
+
+
+# AST_DATE_TO_EPOCH -- Convert Gregorian date and solar mean time to
+# a Julian epoch.  A Julian epoch has 365.25 days per year and 24
+# hours per day.
+
+procedure ast_date_to_epoch (year, month, day, ut, epoch)
+
+int	year			# Year
+int	month			# Month (1-12)
+int	day			# Day of month
+double	ut			# Universal time for date (mean solar day)
+double	epoch			# Julian epoch
+
+double	jd, ast_date_to_julday()
+
+begin
+	jd = ast_date_to_julday (year, month, day, ut)
+	epoch = J2000 + (jd - JD2000) / JYEAR
+end
+
+
+# AST_EPOCH_TO_DATE -- Convert a Julian epoch to year, month, day, and time.
+
+procedure ast_epoch_to_date (epoch, year, month, day, ut)
+
+double	epoch			# Julian epoch
+int	year			# Year
+int	month			# Month (1-12)
+int	day			# Day of month
+double	ut			# Universal time for date
+
+double	jd
+
+begin
+	jd = JD2000 + (epoch - J2000) * JYEAR
+	call ast_julday_to_date (jd, year, month, day, ut)
+end
+
+
+# AST_DAY_OF_YEAR -- The day number for the given year is returned.
+
+int procedure ast_day_of_year (year, month, day)
+
+int	year			# Year
+int	month			# Month (1-12)
+int	day			# Day of month
+
+int	d
+int	bom[13]			# Beginning of month
+data	bom/1,32,60,91,121,152,182,213,244,274,305,335,366/
+
+begin
+	d = bom[month] + day - 1
+	if (month > 2 && mod (year, 4) == 0 &&
+	    (mod (year, 100) != 0 || mod (year, 400) == 0))
+	    d = d + 1
+	return (d)
+end
+
+
+# AST_DAY_OF_WEEK -- Return the day of the week for the given Julian day.
+# The integer day of the week is 0=Sunday - 6=Saturday.  The character string
+# is the three character abbreviation for the day of the week.  Note that
+# the day of the week is for Greenwich if the standard UT is used.
+
+procedure ast_day_of_week (jd, d, name, sz_name)
+
+double	jd		# Julian date
+int	d		# Day of the week (0=SUN)
+char	name[sz_name]	# Name for day of the week
+int	sz_name		# Size of name string
+
+begin
+	d = mod (int (jd - 0.5) + 2, 7)
+	switch (d) {
+	case 0:
+	    call strcpy ("SUN", name, sz_name)
+	case 1:
+	    call strcpy ("MON", name, sz_name)
+	case 2:
+	    call strcpy ("TUE", name, sz_name)
+	case 3:
+	    call strcpy ("WED", name, sz_name)
+	case 4:
+	    call strcpy ("THU", name, sz_name)
+	case 5:
+	    call strcpy ("FRI", name, sz_name)
+	case 6:
+	    call strcpy ("SAT", name, sz_name)
+	}
+end
+
+
+# AST_JULDAY -- Convert epoch to Julian day.
+
+double procedure ast_julday (epoch)
+
+double	epoch			# Epoch
+
+double	jd
+
+begin
+	jd = JD2000 + (epoch - J2000) * JYEAR
+	return (jd)
+end
+
+
+# AST_DATE_TO_JULDAY -- Convert date to Julian day.
+# This assumes dates after year 99.
+
+double procedure ast_date_to_julday (year, month, day, t)
+
+int	year			# Year
+int	month			# Month (1-12)
+int	day			# Day of month
+double	t			# Time for date (mean solar day)
+
+double	jd
+int	y, m, d
+
+begin
+	if (year < 100)
+	    y = 1900 + year
+	else
+	    y = year
+
+	if (month > 2)
+	    m = month + 1
+	else {
+	    m = month + 13
+	    y = y - 1
+	}
+
+	jd = int (JYEAR * y) + int (30.6001 * m) + day + 1720995
+	if (day + 31 * (m + 12 * y) >= 588829) {
+	    d = int (y / 100)
+	    m = int (y / 400)
+	    jd = jd + 2 - d + m
+	}
+	jd = jd - 0.5 + int (t * 360000. + 0.5) / 360000. / 24.
+	return (jd)
+end
+
+
+# AST_JULDAY_TO_DATE -- Convert Julian date to calendar date.
+# This is taken from Numerical Receipes by Press, Flannery, Teukolsy, and
+# Vetterling.
+
+procedure ast_julday_to_date (j, year, month, day, t)
+
+double	j			# Julian day
+int	year			# Year
+int	month			# Month (1-12)
+int	day			# Day of month
+double	t			# Time for date (mean solar day)
+
+int	ja, jb, jc, jd, je
+
+begin
+	ja = nint (j)
+	t = 24. * (j - ja + 0.5)
+
+	if (ja >= 2299161) {
+	    jb = int (((ja - 1867216) - 0.25) / 36524.25)
+	    ja = ja + 1 + jb - int (jb / 4)
+	}
+
+	jb = ja + 1524
+	jc = int (6680. + ((jb - 2439870) - 122.1) / JYEAR)
+	jd = 365 * jc + int (jc / 4)
+	je = int ((jb - jd) / 30.6001)
+	day = jb - jd - int (30.6001 * je)
+	month = je - 1
+	if (month > 12)
+	    month = month - 12
+	year = jc - 4715
+	if (month > 2)
+	    year = year - 1
+	if (year < 0)
+	    year = year - 1
+end
+
+
+# AST_MST -- Mean sidereal time of the epoch at the given longitude.
+# This procedure may be used to optain Greenwich Mean Sidereal Time (GMST)
+# by setting the longitude to 0.
+
+double procedure ast_mst (epoch, longitude)
+
+double	epoch		# Epoch
+double	longitude	# Longitude in degrees
+
+double	jd, ut, t, st
+double	ast_julday()
+
+begin
+	# Determine JD and UT, and T (JD in centuries from J2000.0).
+	jd = ast_julday (epoch)
+	ut = (jd - int (jd) - 0.5) * 24.
+	t = (jd - 2451545.d0) / 36525.d0
+
+	# The GMST at 0 UT in seconds is a power series in T.
+	st = 24110.54841d0 +
+	    t * (8640184.812866d0 + t * (0.093104d0 - t * 6.2d-6))
+
+	# Correct for longitude and convert to standard hours.
+	st = mod (st / 3600. + ut - longitude / 15., 24.0D0)
+
+	if (st < 0)
+	    st = st + 24
+
+	return (st)
+end
diff --git a/noao/astutil/asttools/astvbary.x b/noao/astutil/asttools/astvbary.x
new file mode 100644
index 00000000..826ed38d
--- /dev/null
+++ b/noao/astutil/asttools/astvbary.x
@@ -0,0 +1,76 @@
+include	<math.h>
+
+# AST_VBARY -- Radial velocity component of center of the Earth relative to
+# to the barycenter of the Earth-Moon system.
+
+procedure ast_vbary (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, oblq, omega, llong, lperi, inclin, em, anom, vmoon
+double	r, d, l, b, lm, bm, ast_julday()
+
+begin
+	# T is the number of Julian centuries since J1900.
+	t = (ast_julday (epoch) - 2415020) / 36525.
+
+	# OBLQ is the mean obliquity of the ecliptic
+	# OMEGA is the longitude of the mean ascending node
+	# LLONG is the mean lunar longitude (should be 13.1763965268)
+	# LPERI is the mean lunar longitude of perigee
+	# INCLIN is the inclination of the lunar orbit to the ecliptic
+	# EM is the eccentricity of the lunar orbit (dimensionless)
+	# All quantities except the eccentricity are in degrees.
+
+	oblq = 23.452294d0 -
+	    t * (0.0130125d0 + t * (0.00000164d0 - t * 0.000000503d0))
+	omega = 259.183275d0 -
+	    t * (1934.142008d0 + t * (0.002078d0 + t * 0.000002d0))
+	llong = 270.434164d0 +
+	    t * (481267.88315d0 + t * (-0.001133d0 + t * 0.0000019d0)) - omega
+	lperi = 334.329556d0 +
+	    t * (4069.034029d0 - t * (0.010325d0 + t * 0.000012d0)) - omega
+	em = 0.054900489d0
+	inclin = 5.1453964d0
+
+	# Determine true longitude.  Compute mean anomaly, convert to true
+	# anomaly (approximate formula), and convert back to longitude.
+	# The mean anomaly is only approximate because LPERI should
+	# be the true rather than the mean longitude of lunar perigee.
+
+	lperi = DEGTORAD (lperi)
+	llong = DEGTORAD (llong)
+	anom = llong - lperi
+	anom = anom + (2. * em - 0.25 * em**3) * sin (anom) + 1.25 * em**2 *
+	    sin (2. * anom) + 13./12. * em**3 * sin (3. * anom)
+	llong = anom + lperi
+
+	# L and B are the ecliptic longitude and latitude of the observation.
+	# LM and BM are the lunar longitude and latitude of the observation
+	# in the lunar orbital plane relative to the ascending node.
+
+	r = DEGTORAD (ra * 15)
+	d = DEGTORAD (dec)
+	omega = DEGTORAD (omega)
+	oblq = DEGTORAD (oblq)
+	inclin = DEGTORAD (inclin)
+
+	call ast_coord (double (0.), double (0.), double (-HALFPI),
+	    HALFPI - oblq, r, d, l, b)
+	call ast_coord (omega, double (0.), omega-HALFPI, HALFPI-inclin,
+	    l, b, lm, bm)
+
+	# VMOON is the component of the lunar velocity perpendicular to the
+	# radius vector.  V is the projection onto the line of sight to the
+	# observation of the velocity of the Earth's center with respect to
+	# the Earth-Moon barycenter.  The 81.53 is the ratio of the Earth's
+	# mass to the Moon's mass.
+
+	vmoon = (TWOPI / 27.321661d0) *
+	    384403.12040d0 / sqrt (1. - em**2) / 86400.
+	v = vmoon * cos (bm) * (sin (llong - lm) - em * sin (lperi - lm))
+	v = v / 81.53
+end
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
diff --git a/noao/astutil/asttools/astvr.x b/noao/astutil/asttools/astvr.x
new file mode 100644
index 00000000..f2139b59
--- /dev/null
+++ b/noao/astutil/asttools/astvr.x
@@ -0,0 +1,29 @@
+include	<math.h>
+
+# AST_VR -- Project a velocity vector in radial velocity along line of sight.
+
+procedure ast_vr (ra1, dec1, v1, ra2, dec2, v2) 
+
+double	ra1		# Right ascension of velocity vector (hours)
+double	dec1		# Declination of velocity vector (degrees)
+double	v1		# Magnitude of velocity vector
+double	ra2		# Right ascension of observation (hours)
+double	dec2		# Declination of observation (degrees)
+double	v2		# Radial velocity along direction of observation
+
+double	vx, vy, vz, cc, cs, s
+
+begin
+	# Cartisian velocity components of the velocity vector.
+	vx = v1 * cos (DEGTORAD (15. * ra1)) * cos (DEGTORAD (dec1))
+	vy = v1 * sin (DEGTORAD (15. * ra1)) * cos (DEGTORAD (dec1))
+	vz = v1 * sin (DEGTORAD (dec1))
+
+	# Direction cosines along the direction of observation.
+	cc = cos (DEGTORAD (dec2)) * cos (DEGTORAD (15. * ra2))
+	cs = cos (DEGTORAD (dec2)) * sin (DEGTORAD (15. * ra2))
+	s  = sin (DEGTORAD (dec2))
+
+	# Project velocity vector along the direction of observation.
+	v2 = (vx * cc + vy * cs + vz * s)
+end
diff --git a/noao/astutil/asttools/astvrotate.x b/noao/astutil/asttools/astvrotate.x
new file mode 100644
index 00000000..df5b9b1c
--- /dev/null
+++ b/noao/astutil/asttools/astvrotate.x
@@ -0,0 +1,43 @@
+include	<math.h>
+
+# AST_VROTATE -- Radial velocity component of the observer relative to
+# the center of the Earth due to the Earth's rotation.
+
+procedure ast_vrotate (ra, dec, epoch, latitude, longitude, altitude, v)
+
+double	ra		# Right Ascension of observation (hours)
+double	dec		# Declination of observation (degrees)
+double	epoch		# Epoch of observation (Julian epoch)
+double	latitude	# Latitude (degrees)
+double	longitude	# Latitude (degrees)
+double	altitude	# Altitude (meters)
+double	v		# Velocity (km / s)
+
+double	lat, dlat, r, vc, lmst, ast_mst()
+
+begin
+	# LAT is the latitude in radians.
+	lat = DEGTORAD (latitude) 
+
+	# Reduction of geodetic latitude to geocentric latitude (radians).
+	# Dlat is in arcseconds.
+
+	dlat = -(11. * 60. + 32.743000d0) * sin (2 * lat) +
+		1.163300d0 * sin (4 * lat) -0.002600d0 * sin (6 * lat)
+	lat = lat + DEGTORAD (dlat / 3600.)
+
+	# R is the radius vector from the Earth's center to the observer
+	# (meters).  Vc is the corresponding circular velocity
+	# (meters/sidereal day converted to km / sec).
+	# (sidereal day = 23.934469591229 hours (1986))
+
+	r = 6378160.0d0 * (0.998327073d0 + 0.00167643800d0 * cos (2 * lat) -
+	    0.00000351d0 * cos (4 * lat) + 0.000000008d0 * cos (6 * lat)) +
+	    altitude
+	vc = TWOPI * (r / 1000.)  / (23.934469591229d0 * 3600.)
+
+	# Project the velocity onto the line of sight to the star.
+	lmst = ast_mst (epoch, longitude)
+	v = vc * cos (lat) * cos (DEGTORAD (dec)) *
+	    sin (DEGTORAD ((ra - lmst) * 15.))
+end
diff --git a/noao/astutil/asttools/astvsun.x b/noao/astutil/asttools/astvsun.x
new file mode 100644
index 00000000..c466ab89
--- /dev/null
+++ b/noao/astutil/asttools/astvsun.x
@@ -0,0 +1,38 @@
+include	<math.h>
+
+# The sun's velocity with respect to the local standard of rest is:
+
+define	RAVSUN	18.	# RA of sun's velocity
+define	DECVSUN	30.	# DEC of sun's velocity
+define	VSUN	20.	# VLSR of sun
+define	EPOCH	1900.	# Epoch of sun's velocity
+
+# AST_VSUN -- Projection of the sun's velocity along the given direction.
+
+procedure ast_vsun (ra, dec, epoch, v) 
+
+double	ra		# Reference right ascension (hours)
+double	dec		# Reference declination (degrees)
+double	epoch		# Epoch (years)
+double	v		# VLSR of sun along reference direction
+
+double	ravsun, decvsun, vx, vy, vz, cc, cs, s
+
+begin
+	# Precess VLSR direction to current date.
+	call ast_precess (double (RAVSUN), double (DECVSUN), double (EPOCH),
+	    ravsun, decvsun, epoch)
+
+	# Cartisian velocity components of the sun's velocity.
+	vx = VSUN * cos (DEGTORAD (15. * ravsun)) * cos (DEGTORAD (decvsun))
+	vy = VSUN * sin (DEGTORAD (15. * ravsun)) * cos (DEGTORAD (decvsun))
+	vz = VSUN * sin (DEGTORAD (decvsun))
+
+	# Direction cosines along the reference direction.
+	cc = cos (DEGTORAD (dec)) * cos (DEGTORAD (15. * ra))
+	cs = cos (DEGTORAD (dec)) * sin (DEGTORAD (15. * ra))
+	s  = sin (DEGTORAD (dec))
+
+	# Project sun's motion along the reference direction.
+	v = -(vx * cc + vy * cs + vz * s)
+end
diff --git a/noao/astutil/asttools/mkpkg b/noao/astutil/asttools/mkpkg
new file mode 100644
index 00000000..4c59582c
--- /dev/null
+++ b/noao/astutil/asttools/mkpkg
@@ -0,0 +1,23 @@
+# NOAO astronomical toolslibasttools.a
+
+update:
+	$checkout libasttools.a noaolib$
+	$update   libasttools.a
+	$checkin  libasttools.a noaolib$
+	;
+
+libasttools.a:
+	astarcsep.x	<math.h>
+	astcoord.x	
+	astgalactic.x	<math.h> <mach.h>
+	astgaltoeq.x	<math.h>
+	asthjd.x	<math.h>
+	astlvac.x
+	astprecess.x	<math.h>
+	asttimes.x	
+	astvbary.x	<math.h>
+	astvorbit.x	<math.h>
+	astvr.x		<math.h>
+	astvrotate.x	<math.h>
+	astvsun.x	<math.h>
+	;
diff --git a/noao/astutil/asttools/precessgj.x b/noao/astutil/asttools/precessgj.x
new file mode 100644
index 00000000..4b836c35
--- /dev/null
+++ b/noao/astutil/asttools/precessgj.x
@@ -0,0 +1,48 @@
+include	<math.h>
+
+# PRECESSGJ -- Precess astronomical coordinates from epoch1 to epoch2.
+# Original IRAF/FORTRAN by G. Jacoby (NOAO).
+# Modified by by F. Valdes for IRAF/SPP (NOAO), March 1986
+
+procedure precessgj (ra1, dec1, epoch1, ra2, dec2, epoch2)
+
+double	ra1, dec1, epoch1	# Input coordinates
+double	ra2, dec2, epoch2	# Output coordinates
+
+double	t, tau, theta, zeta, z, ra, dec, a, ap, test
+bool	fp_equald()
+
+begin
+	if (fp_equald (epoch1, epoch2)) {
+	    ra2 = ra1
+	    dec2 = dec1
+	    return
+	}
+
+	t = (epoch2 - epoch1) / 100.
+	tau = (epoch1 - 1850.) / 100.
+
+	theta = t * ((2005.11d0 - 0.85 * tau) - t * (0.43 + t * 0.041)) / 3600.
+	zeta  = t * ((2303.55d0 + 1.40 * tau) + t * (0.30 + t * 0.017)) / 3600.
+	z = zeta + t * t * 0.79 / 3600.
+
+	ra = DEGTORAD (ra1 * 15.0)
+	dec = DEGTORAD (dec1)
+	theta = DEGTORAD (theta)
+	zeta = DEGTORAD (zeta)
+	z = DEGTORAD (z)
+
+	a = ra + zeta
+	dec2 = asin (cos(dec) * cos(a) * sin(theta) + sin(dec) * cos(theta))
+	ap = asin (cos(dec) * sin(a) / cos(dec2))
+	test = (cos(dec)*cos(a)*cos(theta) - sin(dec)*sin(theta)) / cos(dec2)
+
+	if (test < 0.)
+	    ap = PI - ap
+	ra2 = ap + z
+	if (ra2 < 0.)
+	    ra2 = ra2 + TWOPI
+
+	ra2 = RADTODEG (ra2) / 15.0
+	dec2 = RADTODEG (dec2)
+end
diff --git a/noao/astutil/asttools/precessmgb.x b/noao/astutil/asttools/precessmgb.x
new file mode 100644
index 00000000..6a61f9f5
--- /dev/null
+++ b/noao/astutil/asttools/precessmgb.x
@@ -0,0 +1,154 @@
+include	<math.h>
+
+# PRECESSMGB -- Calculate the apparent change in right ascension and
+# declination between two epochs.  Corrections are made for precession,
+# aberration, and some terms of nutation.
+#
+# Based on the work of R. N. Manchester, M. A. Gordon, J. A. Ball (NRAO)
+#     January 1970 (DOPSET, IBM360) 
+# Modified by E. Anderson (NOAO) April 1985 (DOPSET, VMS/VAX)
+# Converted to IRAF by F. Valdes (NOAO) March 1986.
+
+procedure precessmgb (ra1, dec1, epoch1, ra2, dec2, epoch2)
+
+double	ra1, dec1, epoch1	# First epoch coordinates
+double	ra2, dec2, epoch2	# Second epoch coordinates
+
+double	ra, dec, epoch, epoch3, dc
+bool	fp_equald()
+
+begin
+	# Return if the epochs are the same.
+
+	if (fp_equald (epoch1, epoch2)) {
+	    ra2 = ra1
+	    dec2 = dec1
+	    return
+	}
+
+	ra = ra1
+	dec = dec1
+	epoch = epoch1
+
+	# Check if epoch1 is the beginning of a year.  If not make correction
+	# to the beginning of the year.
+
+	if (epoch - int (epoch) > 1E-6) {
+	    epoch3 = epoch
+	    epoch = int (epoch)
+	    call mgb_precess (ra, dec, epoch, ra2, dec2, epoch3, dc)
+	    ra = ra - (ra2 - ra)
+	    dec = dec - (dec2 - dec)
+	}
+
+	# Precess to epoch2.
+
+	call mgb_precess (ra, dec, epoch, ra2, dec2, epoch2, dc)
+end
+
+
+# MGB_PRECESS -- Calculate the apparent change in right ascension and
+# declination between two epochs.  The first epoch is assumed to be the
+# beginning of a year (eg 1950.0).  Also calculate the equation of the
+# equinoxes (in minutes of time) which may be added to the mean siderial
+# time to give the apparent siderial time (AENA-469).  Corrections are
+# made for precession, aberration, and some terms of nutation.
+#    AENA - The American Ephemeris and Nautical Almanac (the blue book)
+#    ESE  - The explanatory suplement to AENA (the green book)
+
+procedure mgb_precess (ra1, dec1, epoch1, ra2, dec2, epoch2, dc)
+
+double	ra1, dec1, epoch1	# First epoch coordinates
+double	ra2, dec2, epoch2	# Second epoch coordinates
+double	dc			# Siderial time correction
+
+double	ra, dec, delr, deld
+double	t1, t2, nday2
+double	theta, zeta, z, am, an, al
+double	csr, snr, csd, snd, tnd, csl, snl
+double	omega, dlong, doblq
+bool	fp_equald()
+
+begin
+	# Check if epochs are the same.
+
+	if (fp_equald (epoch1, epoch2)) {
+	    ra2 = ra1
+	    dec2 = dec1
+	    return
+	}
+
+	# Convert input coordinates to radians.
+
+	ra = DEGTORAD (ra1 * 15.)
+	dec = DEGTORAD (dec1)
+
+	# Compute sines, cosines, and tangents.
+
+	csr = cos (ra)
+	snr = sin (ra)
+	csd = cos (dec)
+	snd = sin (dec)
+	tnd = snd / csd
+
+	# T1 is the time from 1900 to epoch1 (centuries),
+	# t2 is the time from epoch1 to epoch2 (centuries), and
+	# nday2 is the number of days since the beginning of the year
+	# for epoch2.  The number of ephemeris days in a tropical
+	# year is 365.2421988.
+
+	t1 = (epoch1 - 1900.) / 100.
+	t2 = (epoch2 - epoch1) / 100.
+	nday2 = (epoch2 - int (epoch2)) * 365.2421988
+
+	# Theta, zeta, and z are precessional angles from ESE-29 (arcseconds).
+
+	theta = t2 * ((2004.682 - 0.853 * t1) - t2 * (0.426 + t2 * 0.042))
+	zeta = t2 * ((2304.250 + 1.396 * t1) + t2 * (0.302 + t2 * 0.018))
+	z = zeta + 0.791 * t2 ** 2
+
+	# am and an are the M and N precessional numbers (see AENA-50, 474)
+	# (radians) and alam is an approximate mean longitude for the sun
+	# (AENA-50) (radians)
+
+	am = DEGTORAD ((zeta + z) / 3600.)
+	an = DEGTORAD (theta / 3600.)
+	al = DEGTORAD (0.985647 * nday2 + 278.5)
+
+	snl = sin (al)
+	csl = cos (al)
+
+	# Delr and deld are the annual aberation term in ra and dec (radians)
+	# (ESE-47,48)  (0.91745051   cos (obliquity of ecliptic))
+	# (0.39784993   sin (obliquity of ecliptic))
+	# (-9.92413605E-5   K   20.47 ARCSECONDS   constant of aberration)
+	# (ESE) plus precession terms (see AENA-50 and ESE-39).
+
+	delr = -9.92413605e-5 * (snl * snr + 0.91745051 * csl * csr) / csd +
+		am + an * snr * tnd
+	deld = -9.92413605e-5 * (snl * csr * snd - 0.91745051 * csl * snr *
+		snd + 0.39784993 * csl * csd) + an * csr
+
+	# The following calculates the nutation (approximately) (ESE-41,45)
+	# Omega is the angle of the first term of nutation (ESE-44)
+	# (approximate formula) (radians).
+	# Dlong is the nutation in longitude (delta-psi) (radians)
+	# Doblq is the nutation in obliquity (delta-epsilon) (radians)
+
+	omega = DEGTORAD (259.183275 - 1934.142 * (t1 + t2))
+	dlong = -8.3597e-5 * sin (omega)
+	doblq =  4.4678e-5 * cos (omega)
+
+	# Add nutation into delr and deld (ESE-43).
+
+	delr = delr + dlong * (0.91745051 + 0.39784993 * snr * tnd) -
+		csr * tnd * doblq
+	deld = deld + 0.39784993 * csr * dlong + snr * doblq
+
+	# Compute new position and the equation of the equinoxes
+	# (dc in minutes of tim, ESE-43)
+
+	ra2 = ra1 + RADTODEG (delr / 15.)
+	dec2 = dec1 + RADTODEG (deld)
+	dc = dlong * 210.264169
+end
diff --git a/noao/astutil/asttools/refrac.x b/noao/astutil/asttools/refrac.x
new file mode 100644
index 00000000..2834ff5f
--- /dev/null
+++ b/noao/astutil/asttools/refrac.x
@@ -0,0 +1,51 @@
+include	<math.h>
+
+
+# REFRAC -- Compute observed place from apparent place.
+#
+# This is a placeholder routine.  I am not completely sure this is
+# done correctly though the SLALIB routines are accurate.
+# This uses the quick (less precise) SLALIB routines for the
+# calculation.
+
+procedure refrac (ara, adec, aha, lat, w, t, p, h, ora, odec)
+
+double	ara, adec			# Apparent ra (hr) and dec (deg)
+double	aha				# Apparent hour angle (hr)
+double	lat				# Latitude (deg)
+double	w				# Effective wavelength (A)
+double	t				# Temperature (C)
+double	p				# Pressure (mbar)
+double	h				# Humidity (frac 0-1)
+double	ora, odec			# Observed ra (hr) and dec (deg)
+
+double	oha, refa, refb, az, el, vu[3], vr[3]
+
+begin
+	# Determine refraction coefficients.
+	call slRFCQ (t+273.15D0, p, h, w/10000D0, refa, refb)
+
+	# Convert (aha,adec) to (az,el).
+	call slDE2H (DEGTORAD(aha*15D0), DEGTORAD(adec), DEGTORAD(lat), az, el)
+
+	# Convert (az,el) to (x,y,z).
+	call slDS2C (az, el, vu)
+
+	# Apply refraction correction.
+	call slREFV (vu, refa, refb, vr)
+
+	# Convert (x,y,z) to (az,el)
+	call slDC2S (vr, az, el)
+
+	# Convert (az,el) to (ha,dec).
+	call slDH2E (az, el, DEGTORAD(lat), oha, odec)
+
+	# Convert (oha,odec) to (ora,odec).
+	oha = RADTODEG(oha) / 15D0
+	odec = RADTODEG(odec)
+	ora = ara + aha - oha
+	if (ara - ora < -12)
+	    ora = ora + 24
+	else if (ara - ora > 12)
+	    ora = ora - 24
+end