from skyfield.api import Loader, Topos
from skyfield import almanac
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, r_moon =  695700., 1738. # 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')
earth, sun, moon = [eph[x] for x in ('earth', 'sun', 'moon')]

time = ts.utc(2017, 8, 21, 18, 50)  # puts the Sun, Earth and Moon in almost the same plane
epos, spos, mpos = [thing.at(time).ecliptic_position().km for thing in (earth, sun, moon)]

dpos = np.array([ 1.283504916395889E+08, -7.695767366817768E+07, -1.657242113438845E+05]) #  Horizons DSCOVR: 2457987.284722222, A.D. 2017-Aug-21 18:50:00.0000
# DSCOVR: 2457987.284722222, A.D. 2017-Aug-21 18:50:00.0000,  1.283504916395889E+08, -7.695767366817768E+07, -1.657242113438845E+05,  1.488346029647132E+01,  2.515585471355211E+01,  4.287686810002356E-02,

# move Sun to origin for shadow producing reasons
epos, spos, mpos, dpos = [pos - spos for pos in (epos, spos, mpos, dpos)]

# rotate around z to xy-plane
theta = np.arctan2(epos[1], epos[0])
# c, s, z, o = [f(-theta) for f in (np.cos, np.sin, np.zeros_like, np.ones_like)]
# R1 = np.array([[c, -s, z], [s, c, z], [z, z, o]])
# epos1, spos1, mpos1 = [(pos*R1).sum(axis=1) for pos in (epos, spos, mpos)]
c, s = [f(-theta) for f in (np.cos, np.sin)]
R1 = np.array([[c, -s, 0], [s, c, 0], [0, 0, 1]])
epos1, spos1, mpos1, dpos1 = [(pos*R1).sum(axis=1) for pos in (epos, spos, mpos, dpos)]

# rotate around y to x-axis
phi = np.arctan2(epos1[2], epos1[0])
# c, s, z, o = [f(-phi) for f in (np.cos, np.sin, np.zeros_like, np.ones_like)]
# R2 = np.array([[c, z, -s], [z, o, z], [s, z, c]])
# epos2, spos2, mpos2 = [(pos*R2).sum(axis=1) for pos in (epos1, spos1, mpos1)]
c, s = [f(-phi) for f in (np.cos, np.sin)]
R2 = np.array([[c, 0, -s], [0, 1, 0], [s, 0, c]])
epos2, spos2, mpos2, dpos2 = [(pos*R2).sum(axis=1) for pos in (epos1, spos1, mpos1, dpos1)]

# move everything to the geocenter

epos3, spos3, mpos3, dpos3 = [pos-epos2 for pos in (epos2, spos2, mpos2, dpos2)]

Re, Rm = 6378., 1737.
Rd = 500. # small but not too small

plt.figure()
colors = 'blue', 'green', 'black'
radii = Re, Rm, Rd
positions = epos3, mpos3, dpos3
th = np.linspace(0, twopi, 201)

plt.figure()
for pos, r, color in zip(positions, radii, colors):
    x = pos[1] + r * np.cos(th)
    y = pos[2] + r * np.sin(th)
    plt.plot(x, y, color)
    plt.gca().set_aspect('equal')
    plt.title('YZ projection; Sun is behind you')
plt.show()

plt.figure()
plt.subplot(2, 1, 1)
for pos, r, color in zip(positions, radii, colors):
    x = pos[0] + r * np.cos(th)
    y = pos[2] + r * np.sin(th)
    plt.plot(x, y, color)
    plt.gca().set_aspect('equal')
    plt.title('XZ projection; Sun is to the left')

plt.subplot(2, 1, 2)
for pos, r, color in zip(positions, radii, colors):
    x = pos[0] + r * np.cos(th)
    y = pos[1] + r * np.sin(th)
    plt.plot(x, y, color)
    plt.gca().set_aspect('equal')
    plt.title('XY projection (top-down); Sun is to the left')
plt.show()