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)
t = ts.tt_jd(2414865) # set time to start of catalog Julian day
print("Time: ", t)
# extract the Julian Date, which is internal to the time object
jd = float(t.tt)
print("Julian Date start: ", 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)

jdays_begin = 2414865
jdays_end = 2471183
jdays_diff = jdays_end - jdays_begin

# 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']

# open output file
out_file = open("output_daily.txt", "w")

# main loop

while (jd < 2471183): # go through end of catalog

  # 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 gives weird results, ignore for now
  # 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)


  # 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

  goodness = 0.0
  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)

  if ((jd % 1000) == 0):
    diff = jdays_end - jd
    percent_done = float(diff) / float(jdays_diff)
    percent_done = int(100.0 - (percent_done * 100.0))
    print(jd, " : ", percent_done, "%")

  # print data to file
  out_file.write(str(jd) + "," + str(goodness) + "\n")

  # increment Julian day
  jd = jd + 1
  t = ts.tt_jd(jd)

# end of main loop

# close output file
out_file.close()

# 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