How to set up an Astropy SkyCoord object to match the True Equinox and Equator Frame

Or how to make RA match LST

The True Equinox and Equator frame is also known as Current Epoch or E_gamma, this is the frame where zenith ra, dec = apparent lst, latitude

Bryna Hazelton

June 14, 2018

import numpy as np
from astropy.coordinates import SkyCoord, Angle, EarthLocation
from astropy.time import Time
from astropy import units

mjd = 55780.1
latitude = Angle('-26d42m11.94986s')
longitude = Angle('116d40m14.93485s')

obs_time = Time(mjd, format='mjd', location = (longitude, latitude))
lst_apparent = obs_time.sidereal_time('apparent')
print(lst_apparent)

7h11m46.2716s

# The frame radio astronomers call the apparent or current epoch is the
# "true equator & equinox" frame, notated E_gamma in the USNO circular 179
# (http://aa.usno.navy.mil/publications/docs/Circular_179.php)
# astropy doesn't have this frame but it's pretty easy to adapt the CIRS frame
# by modifying the ra to reflect the difference between
# GAST (Grenwich Apparent Sidereal Time) and the earth rotation angle (theta).
# The earth rotation angle is most clearly explained in the 2nd paragraph of Chapter 6
# of the USNO circular (page 50)
# The relationship between E_gamma and CIRS is shown in the figure on page 57 of the circular
def egamma_to_cirs_ra(egamma_ra, time):
    from astropy import _erfa as erfa
    from astropy.coordinates.builtin_frames.utils import get_jd12
    era = erfa.era00(*get_jd12(Time(mjd, format='mjd'), 'ut1'))
    theta_earth = Angle(era, unit='rad')

    assert(isinstance(time, Time))
    gast = time.sidereal_time('apparent', longitude=0)
    cirs_ra = egamma_ra - (gast-theta_earth)
    return cirs_ra

# If we define a source at RA = apparent LST and Dec = latitude it should be at zenith

loc_obj = EarthLocation.from_geodetic(lon=longitude, lat=latitude)

# use CIRS but with a different ra to account for the difference between LST & earth's rotation angle
cirs_ra = egamma_to_cirs_ra(lst_apparent, Time(mjd, format='mjd'))

egamma_zenith_coord = SkyCoord(ra=cirs_ra, dec=latitude, frame='cirs',
                               obstime=Time(mjd, format='mjd'), location = loc_obj)
egamma_zenith_coord

<SkyCoord (CIRS: obstime=55780.1): (ra, dec) in deg
    (107.78961213, -26.70331941)>

# check where it is in altaz (should be near zenith):
egamma_zenith_altaz = egamma_zenith_coord.transform_to('altaz')
egamma_zenith_altaz

<SkyCoord (AltAz: obstime=55780.1, location=(-2559302.5737783727, 5095070.526830904, -2848887.400942108) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0, obswl=1.0 micron): (az, alt) in deg
    (28.60732437, 89.99985797)>

(egamma_zenith_altaz.alt - Angle('90d')).to_string(unit=units.degree, sep=('deg', 'm', 's'))

'-0deg00m00.5113s'