# Path: skyfield/api.py
# N = E = +1.0
# S = W = -1.0
# 
# Path: skyfield/magnitudelib.py
# def planetary_magnitude(position):
#     """Given the position of a planet, return its visual magnitude.
# 
#     >>> from skyfield.api import load
#     >>> from skyfield.magnitudelib import planetary_magnitude
#     >>> ts = load.timescale()
#     >>> t = ts.utc(2020, 7, 31)
#     >>> eph = load('de421.bsp')
#     >>> astrometric = eph['earth'].at(t).observe(eph['jupiter barycenter'])
#     >>> print('%.2f' % planetary_magnitude(astrometric))
#     -2.73
# 
#     The formulae are from `Mallama and Hilton “Computing Apparent
#     Planetary Magnitude for the Astronomical Almanac” (2018)
#     <https://arxiv.org/pdf/1808.01973.pdf>`_.  Two of the formulae have
#     inherent limits:
# 
#     * Saturn’s magnitude is unknown and the function will return ``nan``
#       (the floating-point value “Not a Number”) if the “illumination
#       phase angle” — the angle of the vertex observer-Saturn-Sun —
#       exceeds 6.5°.
# 
#     * Neptune’s magnitude is unknown and will return ``nan`` if the
#       illumination phase angle exceeds 1.9° and the position's date is
#       before the year 2000.
# 
#     And one formula is not fully implemented (though contributions are
#     welcome!):
# 
#     * Skyfield does not compute which features on Mars are facing the
#       observer, which can introduce an error of ±0.06 magnitude.
# 
#     """
#     target = position.target
#     function = _FUNCTIONS.get(target)
#     if function is None:
#         name = _target_name(target)
#         raise ValueError('cannot compute the magnitude of target %s' % name)
# 
#     # Shamelessly treat the Sun as sitting at the Solar System Barycenter.
#     sun_to_observer = position.center_barycentric.xyz.au
#     observer_to_planet = position.xyz.au
#     sun_to_planet = sun_to_observer + observer_to_planet
# 
#     r = length_of(sun_to_planet)
#     delta = length_of(observer_to_planet)
#     ph_ang = angle_between(sun_to_planet, observer_to_planet) * RAD2DEG
# 
#     if function is _saturn_magnitude:
#         a = angle_between(_SATURN_POLE, sun_to_planet)
#         sun_sub_lat = a * RAD2DEG - 90.0
# 
#         a = angle_between(_SATURN_POLE, observer_to_planet)
#         observer_sub_lat = a * RAD2DEG - 90.0
# 
#         return function(r, delta, ph_ang, sun_sub_lat, observer_sub_lat)
# 
#     if function is _uranus_magnitude:
#         a = angle_between(_URANUS_POLE, sun_to_planet)
#         sun_sub_lat = a * RAD2DEG - 90.0
# 
#         a = angle_between(_URANUS_POLE, observer_to_planet)
#         observer_sub_lat = a * RAD2DEG - 90.0
# 
#         return function(r, delta, ph_ang, sun_sub_lat, observer_sub_lat)
# 
#     if function is _neptune_magnitude:
#         year = position.t.J
#         return function(r, delta, ph_ang, year)
# 
#     return function(r, delta, ph_ang)


from skyfield.api import load
from skyfield.magnitudelib import planetary_magnitude

def test_magnitudes():
    ts = load.timescale()
    #t = ts.utc(1995, 5, 22)  # Rings edge-on from Earth.
    t = ts.utc(2021, 10, 4)
    eph = load('de421.bsp')
    names = [
        'mercury', 'venus', 'mars', 'jupiter barycenter',
        'saturn barycenter', 'uranus barycenter', 'neptune barycenter',
    ]
    e = eph['earth'].at(t)
    positions = [e.observe(eph[name]) for name in names]

# Path: skyfield/api.py
# N = E = +1.0
# S = W = -1.0
# 
# Path: skyfield/magnitudelib.py
# def planetary_magnitude(position):
#     """Given the position of a planet, return its visual magnitude.
# 
#     >>> from skyfield.api import load
#     >>> from skyfield.magnitudelib import planetary_magnitude
#     >>> ts = load.timescale()
#     >>> t = ts.utc(2020, 7, 31)
#     >>> eph = load('de421.bsp')
#     >>> astrometric = eph['earth'].at(t).observe(eph['jupiter barycenter'])
#     >>> print('%.2f' % planetary_magnitude(astrometric))
#     -2.73
# 
#     The formulae are from `Mallama and Hilton “Computing Apparent
#     Planetary Magnitude for the Astronomical Almanac” (2018)
#     <https://arxiv.org/pdf/1808.01973.pdf>`_.  Two of the formulae have
#     inherent limits:
# 
#     * Saturn’s magnitude is unknown and the function will return ``nan``
#       (the floating-point value “Not a Number”) if the “illumination
#       phase angle” — the angle of the vertex observer-Saturn-Sun —
#       exceeds 6.5°.
# 
#     * Neptune’s magnitude is unknown and will return ``nan`` if the
#       illumination phase angle exceeds 1.9° and the position's date is
#       before the year 2000.
# 
#     And one formula is not fully implemented (though contributions are
#     welcome!):
# 
#     * Skyfield does not compute which features on Mars are facing the
#       observer, which can introduce an error of ±0.06 magnitude.
# 
#     """
#     target = position.target
#     function = _FUNCTIONS.get(target)
#     if function is None:
#         name = _target_name(target)
#         raise ValueError('cannot compute the magnitude of target %s' % name)
# 
#     # Shamelessly treat the Sun as sitting at the Solar System Barycenter.
#     sun_to_observer = position.center_barycentric.xyz.au
#     observer_to_planet = position.xyz.au
#     sun_to_planet = sun_to_observer + observer_to_planet
# 
#     r = length_of(sun_to_planet)
#     delta = length_of(observer_to_planet)
#     ph_ang = angle_between(sun_to_planet, observer_to_planet) * RAD2DEG
# 
#     if function is _saturn_magnitude:
#         a = angle_between(_SATURN_POLE, sun_to_planet)
#         sun_sub_lat = a * RAD2DEG - 90.0
# 
#         a = angle_between(_SATURN_POLE, observer_to_planet)
#         observer_sub_lat = a * RAD2DEG - 90.0
# 
#         return function(r, delta, ph_ang, sun_sub_lat, observer_sub_lat)
# 
#     if function is _uranus_magnitude:
#         a = angle_between(_URANUS_POLE, sun_to_planet)
#         sun_sub_lat = a * RAD2DEG - 90.0
# 
#         a = angle_between(_URANUS_POLE, observer_to_planet)
#         observer_sub_lat = a * RAD2DEG - 90.0
# 
#         return function(r, delta, ph_ang, sun_sub_lat, observer_sub_lat)
# 
#     if function is _neptune_magnitude:
#         year = position.t.J
#         return function(r, delta, ph_ang, year)
# 
#     return function(r, delta, ph_ang)


# Path: skyfield/tests/test_magnitudes.py
from skyfield.api import load
from skyfield.magnitudelib import planetary_magnitude


def test_magnitudes():
    ts = load.timescale()
    #t = ts.utc(1995, 5, 22)  # Rings edge-on from Earth.
    t = ts.utc(2021, 10, 4)
    eph = load('de421.bsp')
    names = [
        'mercury', 'venus', 'mars', 'jupiter barycenter',
        'saturn barycenter', 'uranus barycenter', 'neptune barycenter',
    ]
    e = eph['earth'].at(t)
    positions = [e.observe(eph[name]) for name in names]

    magnitudes = [planetary_magnitude(position) for position in positions]