from skyfield.api import Loader
import numpy as np
import matplotlib.pyplot as plt

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
to_degs, to_rads = 180/pi, pi/180
R_sun =  695700. # kilometers

load = Loader('~/Documents/fishing/SkyData')  # avoids multiple copies of large files

ts = load.timescale() # include builtin=True if you want to use older files (you may miss some leap-seconds)
eph = load('de421.bsp') # Use de405 or another if you want to go past year 2053
earth, sun, mars = [eph[x] for x in ('earth', 'sun', 'mars')]

days = np.arange(0, 53*365.2564, 0.1)
years = 2000 + days/365.2564

times = ts.utc(2000, 1, days)
m = mars.at(times)
e, s = [m.observe(thing) for thing in (earth, sun)]
ems = e.separation_from(s).degrees
dsun = s.distance().km
half_angle_sun = to_degs*np.arctan2(R_sun, dsun)

plt.figure()
plt.rcParams['axes.formatter.useoffset'] = False # https://stackoverflow.com/a/29535832/3904031
plt.subplot(2, 1, 1)
plt.plot(years, ems)
plt.plot(years, half_angle_sun, '-k')
plt.ylim(0, None)
plt.xlim(2000, 2053)
plt.subplot(2, 1, 2)
plt.plot(years, ems)
plt.plot(years, half_angle_sun, '-k')
plt.ylim(0, 1)
plt.xlim(2000, 2053)
plt.suptitle('Sun-Mars-Earth angle and Sun half-angle', fontsize=16)
plt.show()