Epoch conversion of Equatorial Coordinates

Convert equatorial coordinates from one precession epoch to another. Algorithm is from the book «Astronomical Algorithms» by Jean Meeus from 1991. The current standard reference system used is the International Celestial Reference System (ICRS) which reference frame is on the barycenter of the Solar System and fixed to measured centers of extragalactic sources (mostly quasars). It aligns up to a few milli arc seconds to the FK5 (J2000.0) reference system.

epoch.py

from numpy import sin,cos,arcsin,arctan2,deg2rad,rad2deg

# Epochs for Besselian year (B) and Julian year (J)
JDE = {
    "B1850.0": 2396758.203,
    "B1900.0": 2415020.3135, # 1900 January 0.8135
    "B1950.0": 2433282.4235, # 1950 January 0.9235
    "J2000.0": 2451545.00,   # 2000 January 1 12:00 UT
}


def convert_epoch(to_epoch, coords0):
    """
    Converts source coordinates with epoch as array to destination epoch.
    Reads and returns coordinates with epoch, right ascension (hours, decimal)
    and declination (degrees, decimal). Example:
    """

    jd0 = JDE[coords0['epoch']]
    jd  = JDE[to_epoch]
    
    # Convert right ascension from hours to degrees.
    # Convert degrees to radians
    alpha0 = deg2rad(coords0['ra'] * 15.0)
    delta0 = deg2rad(coords0['dec'])
 
    # Let t be the time interval, in Julian centuries, between J2000.0
    # and the starting epoch, and let dt be the interval, in the same
    # units, between the starting epoch and the final epoch.
    # In other words, if jd0 and jd are the Julian Days correspon-
    # ding to the initial and the final epoch, respectively we have 

    t = (jd0 - 2451545.0) / 36525.0
    dt = (jd - jd0) / 36525.0

    # Then we have the following numerical expressions for the quanti-
    # ties zeta, z and theta which are needed for the accurate reduction of po-
    # sitions from one equinox to another. Values are in seconds.
    
    zeta = (2306.2181 + 1.39656 * t - 0.000139 * t**2) * dt \
        + (0.30188 - 0.000344 * t) * dt**2 + 0.017998 * dt**3
    z = (2306.2181 + 1.39656 * t - 0.000139 * t**2) * dt \
        + (1.09468 + 0.000066 * t) * dt**2 + 0.018203 * dt**3
    theta = (2004.3109 - 0.85330 * t - 0.000217 * t**2) * dt \
        - (0.42665 + 0.000217 * t) * dt**2 - 0.041833 * dt**3

    # Convert seconds to radians
    zeta = deg2rad(zeta / 3600.0)
    z = deg2rad(z / 3600.0)
    theta = deg2rad(theta / 3600.0)
 
    # Then, te rigorous formulae for the reduction of the given equa-
    # torial coordinates $alpha0 and $delta0 of the starting epoch to the coordi-
    # nates alpha and delta of the final epoch are:
    
    a = cos(delta0) * sin(alpha0 + zeta)
    b = cos(theta) * cos(delta0) * cos(alpha0 + zeta) - sin(theta) * sin(delta0)
    c = sin(theta) * cos(delta0) * cos(alpha0 + zeta) + cos(theta) * sin(delta0)
    
    # Get equatorial coordinates
    alpha = arctan2(a, b) + z
    delta = arcsin(c)
    
    # Convert from radians to hours or degrees.
    # Keep RA between 0...24
    coords1 = {
        "epoch": to_epoch,
        "ra": (rad2deg(alpha) / 15.0) % 24,
        "dec": rad2deg(delta)
    }
    return(coords1)

sample.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

from epoch import convert_epoch

coords0 = {
    "epoch": "B1950.0",
    "ra": 17 + 2/60.0 + 52.5/3600.0,
    "dec": -(10 + 4/60.0 + 31.0/3600.0)
}

coords1 = convert_epoch("J2000.0", coords0)

print(coords0)
print(coords1)