import datetime
import pandas as pd
from skyfield.api import load
from skyfield import almanac
from skyfield.searchlib import find_discrete

ts = load.timescale()

# set analysis start and end times (Year and month)
t0 = ts.utc(2010,3)
t1 = ts.utc(2060,6)

# set observing season start and end elongations (degrees)
target_elong_start = 60
target_elong_end = 90

# load ephemeris kernels
main_ephemeris = load('de430.bsp')
jupiter_ephemeris = load('jup310.bsp')
jupiter_ephemeris.segments.extend(main_ephemeris.segments)
earth = jupiter_ephemeris['earth']
sun = jupiter_ephemeris['sun']
jup = jupiter_ephemeris['jupiter']

#----------------------------------------------------------------------
# get conjunctions and oppositions
#----------------------------------------------------------------------

# find the dates & times of conjunction and opposition
# t list: dates & times when Jupiter is at conjunction or opposition
# values list:
#   TRUE = opposition
#   FALSE = conjunction
f = almanac.oppositions_conjunctions(jupiter_ephemeris, jup)
t, values = almanac.find_discrete(t0, t1, f)

# created list of event labels
w=['Conjunction'] * len(values)

# Modify list of labels so that w[n]='Opposition' when y[n]>0
for index, e in enumerate(values):
    if e>0:
        w[index]='Opposition'

# Created a dataframe with the dates of opposition and conjunction as a single column
df = pd.DataFrame(t.utc_strftime('%Y-%b-%d %H%M UTC'),
                  columns=['Date-Time'])

# append a column containing the event name
df['Event'] = w
df['JD'] = t.tt

#----------------------------------------------------------------------
# get season starts
#----------------------------------------------------------------------

target_elong = target_elong_start

def jup_quadrature(t):
    e = earth.at(t)
    s = e.observe(sun).apparent()
    j = e.observe(jup).apparent()
    return s.separation_from(j).degrees >= target_elong

# approximate jupiter elongation cycle in days;
# required to make the find_discrete function work
jup_quadrature.rough_period = 380.0

# t list: dates & times when Jupiter is at target elongation
# values list:
#   TRUE = positive elongation (Jupiter leading; observing season start)
#   FALSE = negative elongation (Jupiter following; observing season end)
t, values = find_discrete(t0, t1, jup_quadrature)

# append rows to dataframe, only those dates at +60 elongation, skip -60 elongation results
for ti, vi in zip(t, values):
    if vi!=0:
        df = df.append(pd.Series(data=[ti.utc_strftime('%Y-%b-%d %H%M UTC'),\
             'Observing season starts at elongation +60°', ti.tt], index=df.columns), ignore_index = True)

#----------------------------------------------------------------------
# get season ends
#----------------------------------------------------------------------

target_elong = target_elong_end

# t list: dates & times when Jupiter is at target elongation
# values list:
#   TRUE = positive elongation (Jupiter leading; observing season start)
#   FALSE = negative elongation (Jupiter following; observing season end)
t, values = find_discrete(t0, t1, jup_quadrature)

# append rows to dataframe, only those dates at -90 elongation, skip +90 elongation results
for ti, vi in zip(t, values):
    if vi==0:
        df = df.append(pd.Series(data=[ti.utc_strftime('%Y-%b-%d %H%M UTC'),\
             'Observing season ends at elongation –90°', ti.tt], index=df.columns), ignore_index = True)

#----------------------------------------------------------------------
# sort and print
#----------------------------------------------------------------------

# sort data frame by date & time
df = df.sort_values(by='JD')

# print a nicely formatted list of elongations
# print header
print('Table of Jupiter elongations 2010-2060, from JPL DE430 & JUP310')
print('Positive elongation indicates Jupiter is west of the Sun')
print('Printed ',datetime.datetime.now().strftime('%d %b %Y'),' by AJ4CO Observatory', sep='') #Date is UTC
print()

for row in df.loc[:,:].itertuples():
    if row[2]=='Conjunction':
        print()
        print(row[1],'->',row[2])
        print()
    else:
        print(row[1],'->',row[2])