How to interconvert celestial coordinate systems?

Having defined the prerequisites for celestial coordinate conversion, we can now discuss how to convert celestial coordinates from one system to the other.

Horizontal and Equatorial coordinate systems

1. The horizontal coordinate system uses azimuth and altitude/zenith angle pairs.
2. The equatorial coordinate system uses declination/right ascension or declination/hour angle pairs.

Convert equatorial to horizontal coordinates

def equatorial_to_horizontal(observer_latitude, declination, right_ascension=None, hour_angle=None, local_time=None):
    declination, latitude = dms_to_dd([declination, observer_latitude])
    if right_ascension is not None:
        hour_angle = rightasension2hourangle(right_ascension, local_time)
        hour_angle = hms_to_dd(hour_angle)
    elif hour_angle is not None:
        hour_angle = hms_to_dd(hour_angle)
    elif right_ascension is not None and hour_angle is not None:
        print('Both right_ascension and hour_angle parameters are provided.\nUsing hour_angle for calculations.')
        hour_angle = hms_to_dd(hour_angle)
    else:
        print('Either right_ascension or hour_angle must be provided.')

    latitude, hour_angle, declination = np.radians([latitude, hour_angle, declination])

    zenith_angle = np.arccos(np.sin(latitude) * np.sin(declination) + np.cos(latitude) * np.cos(declination) * np.cos(hour_angle))

    altitude = zenithangle2altitude(zenith_angle, deg=False)

    _num = np.sin(declination) - np.sin(latitude) * np.cos(zenith_angle)
    _den = np.cos(latitude) * np.sin(zenith_angle)
    azimuth = np.arccos(_num / _den)

    if latitude < 0:
        azimuth = np.pi - azimuth
    altitude, azimuth = np.degrees([altitude, azimuth])

    return dd_to_dms(altitude), dd_to_dms(azimuth)

Convert horizontal to equatorial coordinates

def horizontal_to_equatorial(observer_latitude, altitude, azimuth):
    altitude, azimuth, latitude = np.radians([dms_to_dd([altitude, azimuth, observer_latitude])])
    zenith_angle = zenithangle2altitude(altitude)

    zenith_angle = [-zenith_angle if latitude < 0 else zenith_angle][0]

    declination = np.sin(latitude) * np.cos(zenith_angle)
    declination = declination + (np.cos(latitude) * np.sin(zenith_angle) * np.cos(azimuth))
    declination = np.arcsin(declination)

    _num = np.cos(zenith_angle) - np.sin(latitude) * np.sin(declination)
    _den = np.cos(latitude) * np.cos(declination)
    hour_angle = np.arccos(_num / _den)

    if (latitude > 0 > declination) or (latitude < 0 < declination):
        hour_angle = 2 * np.pi - hour_angle

    declination, hour_angle = np.degrees([declination, hour_angle])

    return dd_to_dms(declination), dd_to_hms(hour_angle)

Equatorial <=> Ecliptic

The equatorial coordinate system uses declination/right ascension as the pair of coordinates.

The ecliptic coordinate system uses ecliptic latitude/ecliptic longitude as the pair of coordinates.

Since the ecliptic coordinate system uses the projection of the path traced by the Sun on Earth as the Earth moves around it, it makes an approximate angle of with the equator. This value is incorporated in the code as a constant _ecliptic.

_ecliptic = np.radians(23.43927944)

Convert equatorial to ecliptic coordinates

def equatorial_to_ecliptic(right_ascension, declination):
    ra, dec = hms_to_dd(right_ascension), dms_to_dd(declination)

    ec_latitude = np.arcsin(np.sin(dec) * np.cos(_ecliptic)
    ec_latitude = ec_latitude - (np.cos(dec) * np.sin(_ecliptic) * np.sin(ra)))

    _num, _den = np.sin(ra) * np.cos(_ecliptic) + np.tan(dec) * np.sin(_ecliptic), np.cos(ra)
    ec_longitude = np.arcsin(_num / _den)

    ec_latitude, ec_longitude = np.degrees([ec_latitude, ec_longitude])

    return dd_to_dms(ec_latitude), dd_to_dms(ec_longitude)

Convert ecliptic to equatorial coordinates

def ecliptic_to_equatorial(ecliptic_latitude, ecliptic_longitude):
    ec_latitude, ec_longitude = dms_to_dd([ecliptic_latitude, ecliptic_longitude])

    dec = np.sin(ec_latitude) * np.cos(_ecliptic)
    dec = dec + (np.cos(ec_latitude) * np.sin(_ecliptic) * np.sin(ec_longitude))
    dec = np.arcsin(dec)

    _num = np.cos(ec_latitude) * np.sin(ec_longitude) * np.cos(_ecliptic)
    _num = _num - (np.sin(ec_latitude) * np.sin(_ecliptic))
    _den = np.cos(ec_latitude) * np.cos(ec_longitude)
    ra = np.arctan2(_num, _den)

    if ra < 0:
        ra = ra + 2 * np.pi

    ra, dec = np.degrees([ra, dec])

    return dd_to_hms(ra), dd_to_dms(dec)

Equatorial <=> Galactic

The equatorial coordinate system uses declination/right ascension as the pair of coordinates.

The galactic coordinate system uses galactic latitude/galactic longitude as the pair of coordinates.

As the north galactic poles and the celestial north pole are used to reference these coordinate systems, the RA/Dec and longitude of these references are given below

1. Right ascension of north galactic pole: .
2. Declination of north galactic pole: .
3. Longitude of north celestial pole: .

Before performing the equatorial to galactic coordinate conversions, these constants (_ra_gal and _dec_gal) are introduced in the code.

_ra_gal, _dec_gal, _long_ncp = np.radians([192.85948, 27.12825, 122.93192])

Convert equatorial to galactic coordinates

def equatorial_to_galactic(right_ascension, declination):
    ra, dec = hms_to_dd(right_ascension), dms_to_dd(declination)

    gal_long = np.cos(dec) * np.sin(ra - _ra_gal), np.sin(dec) * np.cos(_dec_gal)
    gal_long = gal_long - (np.cos(dec) * np.sin(_dec_gal) * np.cos(ra - _ra_gal))
    gal_long = np.arctan2(gal_long) - _long_npc

    gal_lat = np.arcsin(np.sin(dec) * np.sin(_dec_gal)
    gal_lat = gal_lat + np.cos(dec) * np.cos(_dec_gal) * np.cos(ra - _ra_gal))

    gal_long, gal_lat = np.degrees([gal_long, gal_lat])

    return dd_to_dms(gal_long), dd_to_dms(gal_lat)

Convert galactic to equatorial coordinates

def galactic_to_equatorial(galactic_latitude, galactic_longitude):
    gal_lat, gal_long = dms_to_dd([galactic_latitude, galactic_longitude])

    dec = np.sin(gal_lat) * np.sin(_dec_gal)
    dec = dec + (np.cos(gal_lat) * np.cos(_dec_gal) * np.cos(_long_ncp - gal_long))
    dec = np.arcsin(dec)

    _ra = np.sin(gal_lat) * np.cos(_dec_gal)
    _ra = _ra - (np.cos(gal_lat) * np.sin(_dec_gal) * np.cos(_long_ncp - gal_long))
    ra = np.arctan2(np.cos(gal_lat) * np.sin(_long_ncp - gal_long), _ra)
    ra = ra + _ra_gal

    ra, dec = np.degrees([ra, dec])

    return dd_to_hms(ra), dd_to_dms(dec)

These codes are also available in the GitLab repository.