from skyfield.api import load 
from skyfield.api import Topos
from datetime import datetime
from pytz import timezone
import matplotlib.pyplot as plt
import math
import numpy as np

ts = load.timescale()
eph = load('de421.bsp')
jupiter, saturn, earth = eph['jupiter barycenter'], eph['saturn barycenter'], eph['earth']

morden = earth + Topos('49.2 N', '98.1 W')

def jupiter_saturn_separation(t):
    e = earth.at(t)
    j = e.observe(jupiter).apparent()
    s = e.observe(saturn).apparent()
    return j.separation_from(s).degrees

def jupiter_saturn_horizon_coordinates(t):
    tc = t.astimezone(central)
    jp = earth.at(t).observe(jupiter)
    jra, jdec, jdist = jp.radec()
    astrometric = morden.at(t).observe(jupiter)
    jalt, jaz, jd = astrometric.apparent().altaz()
    sp = earth.at(t).observe(saturn)
    sra, sdec, sdist = sp.radec()
    astrometric = morden.at(t).observe(saturn)
    salt, saz, sd = astrometric.apparent().altaz()
    return tc, jalt.degrees, jaz.degrees, salt.degrees, saz.degrees

central = timezone('US/Central')

t = ts.utc(2020, 12, range(1,31))

ei = jupiter_saturn_separation(t)
tcdt, jalt, jaz, salt, saz = jupiter_saturn_horizon_coordinates(t)

def midpoint(x1, y1, x2, y2):
    return ((x1 + x2)/2, (y1 + y2)/2)

plt.figure(figsize=(9,9))

plt.style.use("dark_background")
plt.title('Great Conjunction 2020\nJupiter/Saturn Angular Separation\n18:00 CDT Morden, Manitoba, Canada')
plt.xlabel('Azimuth (Degrees)')
plt.ylabel('Altitude (Degrees)')
plt.grid(True, linewidth=0.5, linestyle='dotted')
plt.plot(jaz, jalt, 'o', markersize=6.6, color=(.855, .765, .592), label='Jupiter')
plt.plot(saz, salt, 'o', markersize=3.0, color=(.914, .776, .376), label='Saturn')
plt.legend(loc="upper right")
plt.axis('square')
xmax = max(int(max(jaz)+1),int(max(saz)+1))
xmin = min(int(min(jaz)),int(min(saz)))
ymax = max(int(max(jalt)+1),int(max(salt)+1))
ymin = min(int(min(jalt)),int(min(salt)))
plt.xlim(221,232)
plt.ylim(0,11)
xint = range(221, xmax) 
plt.xticks(xint)
yint = range(0, 11) 
plt.yticks(yint)

degree_sign= u'\N{DEGREE SIGN}'

for i in range(len(jaz)):
    x = (jaz[i], saz[i])
    y = (jalt[i], salt[i]) 
    plt.plot(x, y, color='cyan', linewidth=0.5)
    md = midpoint(x[0], y[0], x[1], y[1])
    plt.annotate(str(round(ei[i],4)) + degree_sign + " " + tcdt[i].strftime('%Y-%m-%d'), # this is the text
                 (md[0], md[1]), # this is the point to label
                 size=7,
                 textcoords="offset points", # how to position the text
                 xytext=(5,-2), # distance from text to points (x,y)
                 ha='left') # horizontal alignment can be left, right or center

plt.savefig('great_conjunction_2020_skyfield_matplotlib_graph_800x800.svg',bbox_inches='tight')