from skyfield.api import load
import math

# using this guy: https://rhodesmill.org/skyfield/
# which is installed via pip: https://pypi.org/project/skyfield/
# for all the computational heavy-lifting


# Time object
# Note: Four ways to specify 2014 January 18 01:35:37.5
# t1 = ts.utc(2014, 1, 18.06640625)
# t2 = ts.utc(2014, 1, 18, 1.59375)
# t3 = ts.utc(2014, 1, 18, 1, 35.625)
# t4 = ts.utc(2014, 1, 18, 1, 35, 37.5)

# we can also specify the Julian date directly:
#  t = ts.tt_jd(2456675.56640625)

ts = load.timescale()
# t = ts.now()
t = ts.utc(1900, 1, 1)
print("Time: ", t)
# extract the Julian Date, which is internal to the time object
jd = float(t.tt)
print("Julian Date: ", jd)
# Also see: https://en.wikipedia.org/wiki/Julian_day


# Load the planetary ephemeris catalog
planets = load('de421.bsp')
# print(planets)
# Note: this catalog looks like it has a date range:   
#   JD 2414864.50 - JD 2471184.50  (1899-07-28 through 2053-10-08)


# Declare objects for all planets and sun
sun = planets['sun']
mercury = planets['mercury']
venus = planets['venus']
earth = planets['earth']
mars = planets['mars']
jupiter = planets['jupiter barycenter']
saturn = planets['saturn barycenter']
uranus = planets['uranus barycenter']
neptune = planets['neptune barycenter']
pluto = planets['pluto barycenter']


# From the center of the Sun (Heliocentric)
hc_sun = sun.at(t).observe(sun)
hc_mercury = sun.at(t).observe(mercury)
hc_venus = sun.at(t).observe(venus)
hc_earth = sun.at(t).observe(earth)
hc_mars = sun.at(t).observe(mars)
hc_jupiter = sun.at(t).observe(jupiter)
hc_saturn = sun.at(t).observe(saturn)
hc_uranus = sun.at(t).observe(uranus)
hc_neptune = sun.at(t).observe(neptune)
hc_pluto = sun.at(t).observe(pluto)

# calculate orbital right-ascention position
#   note: if you want fractional hours, do hours = ra.hours

# sun-from-the-sun causes a warning, but gives zero, which is
#  what we'd expect.
# ra, dec, distance = hc_sun.radec()
# sun_degrees = ra._degrees
ra, dec, distance = hc_mercury.radec()
mercury_degrees = ra._degrees
ra, dec, distance = hc_venus.radec()
venus_degrees = ra._degrees
ra, dec, distance = hc_earth.radec()
earth_degrees = ra._degrees
ra, dec, distance = hc_mars.radec()
mars_degrees = ra._degrees
ra, dec, distance = hc_jupiter.radec()
jupiter_degrees = ra._degrees
ra, dec, distance = hc_saturn.radec()
saturn_degrees = ra._degrees
ra, dec, distance = hc_uranus.radec()
uranus_degrees = ra._degrees
ra, dec, distance = hc_neptune.radec()
neptune_degrees = ra._degrees
ra, dec, distance = hc_pluto.radec()
pluto_degrees = ra._degrees

ppos = []

# print("Sun degrees: ", sun_degrees)
print("Mercury Degrees: ", mercury_degrees)
ppos.append(mercury_degrees)
print("Venus Degrees: ", venus_degrees)
ppos.append(venus_degrees)
print("Earth Degrees: ", earth_degrees)
ppos.append(earth_degrees)
print("Mars Degrees: ", mars_degrees)
ppos.append(mars_degrees)
print("Jupiter Degrees: ", jupiter_degrees)
ppos.append(jupiter_degrees)
print("Saturn Degrees: ", saturn_degrees)
ppos.append(saturn_degrees)
print("Uranus Degrees: ", uranus_degrees)
ppos.append(uranus_degrees)
print("Neptune Degrees: ", neptune_degrees)
ppos.append(neptune_degrees)
print("Pluto Degrees: ", pluto_degrees)
ppos.append(pluto_degrees)

goodness = 0.0

# calculate closeness to 90 for each pair of planetary positions
#  using absolute value of cosine of 2 * degrees, which should be
#  1.0 at multiples of 90 degrees

for i in range(0, len(ppos) - 1):
  for j in range(i + 1, len(ppos)):
    goodness += math.fabs(math.cos(2 * math.radians(math.fabs(ppos[i] - ppos[j]))))

print("Goodness function: ", goodness)

# to do: loop through entire ephemeris table dates, calculating 
#  starting date and then just increment t.tt by Julian Date 
#  equivalent of an hour (0.7417666666).  Cacl planetary positions with goodness
#  function to maximize at multiples of 90 degrees as per
#  J. H. Nelson's paper.  Save resultant output with JD.
#  Correlate with Kp index or other ionospheric reflectance
#  data, or number of sunspots, or reported dates of poor/good propagation