From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- noao/astutil/asttools/precessmgb.x | 154 +++++++++++++++++++++++++++++++++++++ 1 file changed, 154 insertions(+) create mode 100644 noao/astutil/asttools/precessmgb.x (limited to 'noao/astutil/asttools/precessmgb.x') 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 + +# 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 -- cgit