import ephem
import numpy as np

bodies = [ephem.Sun(),
          ephem.Moon(),
          ephem.Mars(),
          ephem.Venus(),
          ephem.Mars(),
          ephem.Jupiter(),
          ephem.Saturn(),
          ephem.Uranus(),
          ephem.Neptune(),
          ephem.Pluto()]

mydate = '2026/1/1'  # can be changed to whatever

def sexi2deci(sexi):
# convert sexigesimal to decimal using the above formula
    splits = str(sexi).split(':')
    deci = abs(int(splits[0])) + int(splits[1])/60 + float(splits[2])/3600
    if str(sexi)[0]=='-':
        deci = -deci
    return deci

# sanity check using values from https://rhodesmill.org/pyephem/tutorial.html
body = bodies[7]
body.compute('1781/3/13')
ra, dec, mag = str(body.ra), str(body.dec), body.mag
assert body.name=='Uranus' and ra=='5:35:45.28' and dec=='23:32:54.1' and mag==5.6

# quick-and-dirty demo: elongation, magnitude, phase and constellation
print('{:8s} {:>8s} {:>10s} {:>7s}  {}'.format('', 'ELONGATION', 'MAGNITUDE', 'PHASE', 'CONSTELLATION'))
print('{:8s} {:>8s} {:>10s} {:>7s}  {}'.format('', ' (radians)', '', '(%)', 'constellation'))
for body in bodies:
    body.compute(mydate)
    print('{:10s} {:>8.4f} {:>10.2f} {:7.2f}  {}'.format(body.name, body.elong, body.mag, body.phase, ephem.constellation(body)))

# RA: astrometric geocentric, apparent geocentric, apparent topocentric
print('{:10s} {:>20s} {:>22s} {:>23s}'.format('', 'ASTROMETRIC_GEOCENTRIC', 'APPARENT_GEOCENTRIC', 'APPARENT_TOPOCENTRIC'))
print('{:10s} {:>12s} {:5s} {:>12s} {:5s} {:>13s} {:5s}'.format('', '(degrees)', '(radians)', '(degrees)', '(radians)', '(degrees)', '(radians)'))
for body in bodies:
    body.compute(mydate)
    a_ra, g_ra, ra = body.a_ra, body.g_ra, body.ra
    print('{:10s} {:>12s} {:9.3f} {:>12s} {:9.3f} {:>13s} {:9.3f}'.format(body.name, str(a_ra), a_ra, str(g_ra), g_ra, str(ra), ra))
# not much difference between the three, as expected

# "Angles become strings when printed or given to str(), but otherwise act like Python floating point numbers." [https://rhodesmill.org/pyephem/tutorial.html]

# RA: manually convert degrees to radians, just to get a feel
print('{:10s} {:>32s}'.format('', 'APPARENT_TOPOCENTRIC'))
print('{:10s} {:>12s} {:5s} {:5s}'.format('', '(degrees)', '(radians)', '(radians)'))
for body in bodies:
    body.compute(mydate)
    ra = body.ra
    manual = 15*sexi2deci(ra)
    print('{:10s} {:>12s} {:9.3f} {:9.3f}'.format(body.name, str(ra), ra, np.deg2rad(manual)))

# Dec: manually convert degrees to radians, just to get a feel
print('{:10s} {:>32s}'.format('', 'DECLINATION'))
print('{:10s} {:>12s} {:5s} {:5s}'.format('', '(degrees)', '(radians)', '(radians)'))
for body in bodies:
    body.compute(mydate)
    dec = body.dec
    manual = np.deg2rad(sexi2deci(dec))
    assert np.allclose(dec, manual, rtol=1e-4)
    print('{:10s} {:>12s} {:>9.3f} {:>9.3f}'.format(body.name, str(dec), dec, manual))

# heliocentric vs geocentric separations
# manually calculate geocentric separation, just to get a feel
for body in bodies:
    body.compute(mydate)
    xx, sun = body, bodies[0]   # bodies[0] is the sun
    xx_ra, xx_dec, sun_ra, sun_dec = xx.ra, xx.dec, sun.ra, sun.dec

    # Haversine formula
    manual = np.sin(xx_dec)*np.sin(sun_dec) + np.cos(xx_dec)*np.cos(sun_dec)*np.cos(xx_ra-sun_ra)
    manual = np.rad2deg(np.arccos(manual))
    
    print('{:8s} {:>12s} [{:6.2f}] \t {:>12s}'.format(body.name, str(ephem.separation(xx, sun)), manual,
                                               str(ephem.separation((xx.hlon, xx.hlat), (sun.hlon, sun.hlat)))))