# We need to import some toolboxes 
import pykep as pk
from matplotlib import pyplot as plt
import numpy as np
%matplotlib notebook

"""
Explanation: Using SPICE kernels to study HERA trajectory
In this tutorial we will use the spice kernels available for download at ftp://spiftp.esac.esa.int/data/SPICE/ to study the interplanetary trajectory of the HERA mission
End of explanation
"""


# Assuming the kernels have been downloaded and
# placed in the same directory of this notebook
# 1 - We load the solar system ephemerides
pk.util.load_spice_kernel("de432s.bsp")
# 2 - We load the ephemerides for HERA trajectory
pk.util.load_spice_kernel("HERA_sc_LPC_EMA_2024c.bsp")
# 3 - We load the ephemerides for the Didymos position
pk.util.load_spice_kernel("didymos_hor_200101_300101_v01.bsp")

"""
Explanation: Loading the SPICE kernels
End of explanation
"""


# We now need to define an epoch (date) to query for positions of the various bodies
# Keep in mind that the kernels do not contain the information for all epochs, only for some
# ranges. You will get an error if you query outside of these bounds
ep = pk.epoch_from_string("2026-Jan-04 12:00:00")
print(ep)

# The code assigned to HERA seems to be -999, so we stick to it. Here we create a fake "planet"
# describing the position of HERA w.r.t. the SUN in the J2000 ecliptic frame and no aberration.
# (ignore the parameters 100 .. 100 and 100 they describe [in pykep syntax] the mass, dimensions and safe
# distance of planets which in this case makes no sense)
hera = pk.planet.spice('-999', 'SUN', 'ECLIPJ2000', 'NONE', pk.MU_SUN, 100, 100, 100) 

r, v = hera.eph(ep)
print("Position vector at epoch: ", r)
print("Velocity vector at epoch: ", v)

# We want also the position with respect to didymos? The code assigned to Didymos seems 
# to be 2065803 (apparently using names is not fashionable). So we define a new fake planet object
hera_wrt_didymos = pk.planet.spice('-999', '2065803', 'ECLIPJ2000', 'NONE', pk.MU_SUN, 100, 100, 100) 

r, v = hera_wrt_didymos.eph(ep)
print("Position vector at epoch: ", r)
print("Velocity vector at epoch: ", v)

# We want also the position with respect to the Earth? We define a new fake planet object
hera_wrt_earth = pk.planet.spice('-999', 'EARTH', 'ECLIPJ2000', 'NONE', pk.MU_SUN, 100, 100, 100) 

r, v = hera_wrt_earth.eph(ep)
print("Position vector at epoch: ", r)
print("Velocity vector at epoch: ", v)

"""
Explanation: A simple use
End of explanation
"""


# The documentation of HERA_sc_LPC_EMA_2024c.bsp kernel reveals that the kernel contains the following data:
#
# HERA_sc_LPC_EMA_2024c.bsp
# Bodies                                    Start of Interval (UTC)         End of Interval (UTC)
# -------                                   -----------------------------   -----------------------------
# -999 HERA w.r.t. 4 MARS BARYCENTER        2025-MAR-17 10:05:00.284        2025-MAR-18 22:29:36.525
# -999 HERA w.r.t. 10 SUN                   2024-OCT-27 11:26:16.873        2025-MAR-17 10:05:00.284
#                                           2025-MAR-18 22:29:36.525        2026-DEC-28 11:58:50.816
# -999 HERA w.r.t. 399 EARTH                2024-OCT-25 11:58:50.817        2024-OCT-27 11:26:16.873
#
# In reality it seems the data starts at 2024-OCT-25 12:00:00. So we use that as launch date (anything before seems
# to be absent in the data)

launch = pk.epoch_from_string("2024-OCT-25 12:00:00")
arrival = pk.epoch_from_string("2026-DEC-28 11:58:50.816")
flyby = pk.epoch_from_string("2025-MAR-17 10:05:00.284")


# We reinstantiate the spacecraft orbit for convenience.
hera = pk.planet.spice('-999', 'SUN', 'ECLIPJ2000', 'NONE', pk.MU_SUN, 100, 100, 100) 
# These are the solar system plants anorbiting objects we are intrerested in
earth = pk.planet.spice('EARTH', 'SUN', 'ECLIPJ2000', 'NONE', pk.MU_SUN, pk.MU_EARTH, pk.EARTH_RADIUS, pk.EARTH_RADIUS*1.05) 
mars = pk.planet.spice('4', 'SUN', 'ECLIPJ2000', 'NONE', pk.MU_SUN, 100, 100,100) 
didymos = pk.planet.spice('2065803', 'SUN', 'ECLIPJ2000', 'NONE', pk.MU_SUN, 100, 100, 100) 

# And these are the plotting commands
ax = pk.orbit_plots.plot_planet(hera, t0 = launch, tf = arrival, N=500, units=pk.AU)
pk.orbit_plots.plot_planet(earth, t0 = launch, axes=ax, color = 'b', units=pk.AU)
pk.orbit_plots.plot_planet(mars, t0 = flyby, axes=ax, color='r', units=pk.AU)
pk.orbit_plots.plot_planet(didymos, t0 = launch, axes=ax, color='g', units=pk.AU)
ax.set_xlim(-4,4)
ax.set_ylim(-4,4)

"""
Explanation: Visualizing the whole trajectory
End of explanation
"""


# Say we are interested in the HERA-EARTH distance.
tspan = np.linspace(launch.jd, arrival.jd, 100)
Reh = []
Red = []
for jd in tspan:
    r1, _ = hera_wrt_earth.eph(pk.epoch(jd, 'jd'))
    r2, _ = hera_wrt_didymos.eph(pk.epoch(jd, 'jd'))
    Reh.append(np.linalg.norm(r1) / pk.AU)
    Red.append(np.linalg.norm(r2) / pk.AU)



days_from_launch = tspan - launch.jd
plt.figure()
plt.plot(days_from_launch, Reh, label='SC / Earth distance (AU)')
plt.plot(days_from_launch, Red, label='SC / Didymos distance (AU)')
plt.legend()

"""
Explanation: Plotting other relevant quantities
End of explanation
"""


