1 files changed, 154 insertions, 0 deletions
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