Epochen-Umrechnung von äquatorialen Koordinaten

Konversion von äquatorialen Koordinaten einer Präzessionsepoche in eine andere. Der Algorithmus stammt aus dem Buch «Astronomical Algorithms» von Jean Meeus aus dem Jahr 1991. Das derzeit verwendete Standard-Referenzsystem ist das International Celestial Reference System (ICRS), dessen Bezugsrahmen auf dem Baryzentrum des Sonnensystems liegt und an den gemessenen Zentren extragalaktischer Quellen ( meist Quasare) ausgerichtet ist. Es stimmt bis auf wenige Millibogensekunden mit Bezugssystem FK5 (J2000.0) überein.

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)