""" JUICE senario 1.
Jupiter pointing senario.
Jupiter pointing phase, the spacecraft
"""
import datetime
import irfpy.juice.jspice as js
import matplotlib.dates as mpldates
import numpy as np
import os
# spicepath = '../apps160221_juice/'
# js.spice.furnsh(os.path.join(spicepath, 'de421.bsp'))
# js.spice.furnsh(os.path.join(spicepath, 'mantra_juice_jup_a5d_141a_lau_c5e_016.bsp'))
# js.spice.furnsh(os.path.join(spicepath, 'naif0009.tls'))
# js.spice.furnsh(os.path.join(spicepath, 'pck00010.tpc'))
# js.spice.furnsh(os.path.join(spicepath, 'jup230l.bsp'))
# js.spice.furnsh(os.path.join(spicepath, 'jse_111130.tf'))
js.spice.furnsh('scenario.mk')
n_jp = js.spice.bodn2c('JUPITER')
n_sun = js.spice.bodn2c('SUN')
n_earth = js.spice.bodn2c('EARTH')
n_ganymede = js.spice.bodn2c('GANYMEDE')
n_io = js.spice.bodn2c('IO')
n_europa = js.spice.bodn2c('EUROPA')
n_callisto = js.spice.bodn2c('CALLISTO')
rj = 71492.
import irfpy.juice.attitude0 as jatt
import matplotlib.pyplot as plt
def length(v):
return np.sqrt(np.inner(v, v))
def unit_vector(v):
return v / np.sqrt(np.inner(v, v))
def longitude(v):
return np.rad2deg(np.arctan2(v[1], v[0]))
def latitude(v):
vv = unit_vector(v)
return np.rad2deg(np.arcsin(vv[2]))
def get_state(n_target, et, frame, n_origin):
posvel, lt = js.spice.spkez(n_target, et, frame, 'LT+S', n_origin)
return posvel[:3] / rj, posvel[3:]
def main():
fp = open('scenario1.txt', 'w')
print('#datetime mpldates et vx_jse vy_jse vz_jse vx_sc vy_sc vz_sc', file=fp)
print('# Velocity is measured in the spacecraft frame', file=fp)
pathname = 'scenario1'
try:
os.mkdir(pathname)
except:
pass
### Jupiter orbit from 2031-05-10T17:00:00 for 30 days
t0 = datetime.datetime(2031, 5, 10, 17)
t1 = datetime.datetime(2031, 6, 10, 17)
dt = datetime.timedelta(minutes=60*2)
t = t0 - dt
# t = t1 - dt
# juice_frame = JuiceNominal() # Convert between SC frame and JSE frame
while t <= t1:
t = t + dt
et = js.spice.str2et(t.strftime('%FT%T'))
### Ok, now many object positions can be calculated.
### The variable names are "target"_"frame"_"origin"
### Here SC centered JSE positions
# jupiter_jse_juice, v_jupiter_jse_juice = get_state(n_jp, et, "JSE", -907)
sun_jse_juice, v_sun_jse_juice = get_state(n_sun, et, "JSE", -907)
# earth_jse_juice, __vel = get_state(n_earth, et, "JSE", -907)
# io_jse_juice, __vel = get_state(n_io, et, "JSE", -907)
# europa_jse_juice, __vel = get_state(n_europa, et, "JSE", -907)
# ganymede_jse_juice, __vel = get_state(n_ganymede, et, "JSE", -907)
# callisto_jse_juice, __vel = get_state(n_callisto, et, "JSE", -907)
### Here JSE positions
juice_jse_jupiter, v_juice_jse_jupiter = get_state(-907, et, "JSE", n_jp)
io_jse_jupiter, __vel = get_state(n_io, et, "JSE", n_jp)
europa_jse_jupiter, __vel = get_state(n_europa, et, "JSE", n_jp)
ganymede_jse_jupiter, __vel = get_state(n_ganymede, et, "JSE", n_jp)
callisto_jse_jupiter, __vel = get_state(n_callisto, et, "JSE", n_jp)
sun_jse_jupiter, __vel = get_state(n_sun, et, "JSE", n_jp)
sun_jse_jupiter_unitvector = unit_vector(sun_jse_jupiter)
earth_jse_jupiter, __vel = get_state(n_earth, et, "JSE", n_jp)
earth_jse_jupiter_unitvector = unit_vector(earth_jse_jupiter)
### Here System3 positions (many velocity not used)
juice_sys3_jupiter, v_juice_sys3_jupiter = get_state(-907, et, "IAU_JUPITER", n_jp)
io_sys3_jupiter, __vel = get_state(n_io, et, "IAU_JUPITER", n_jp)
europa_sys3_jupiter, __vel = get_state(n_europa, et, "IAU_JUPITER", n_jp)
ganymede_sys3_jupiter, __vel = get_state(n_ganymede, et, "IAU_JUPITER", n_jp)
callisto_sys3_jupiter, __vel = get_state(n_callisto, et, "IAU_JUPITER", n_jp)
sun_sys3_jupiter, __vel = get_state(n_sun, et, "IAU_JUPITER", n_jp)
sun_sys3_jupiter_unitvector = unit_vector(sun_sys3_jupiter)
earth_sys3_jupiter, __vel = get_state(n_earth, et, "IAU_JUPITER", n_jp)
earth_sys3_jupiter_unitvector = unit_vector(earth_sys3_jupiter)
### Conversion matrixes.
cmat_eclip_jse = js.spice.pxform('ECLIPJ2000', 'JSE', et) # Ecliptic -> JSE
cmat_sys3_jse = js.spice.pxform('IAU_JUPITER', 'JSE', et) ### System III to JSE frame conversion matrix
### Earth's orbital plane.
eclip_norm_vec = cmat_eclip_jse.dot(np.array([0, 0, 1])) # Almost along z by definition.
### Spacecraft frame definition
# juice_frame.set_configuration(juice_jse_jupiter,
# ecliptic_normal=eclip_norm_vec,
# sundir=sun_jse_juice,
# vel=v_juice_jse_jupiter)
juice_frame = jatt.get_juice_jupiter_pointing(juice_jse_jupiter,
eclip_norm_vec,
sun_jse_juice,
vel=v_juice_jse_jupiter)
### Objects in Spacecraft frame
# jupiter_sc_juice = juice_frame.convert_to_nsc(jupiter_jse_juice)
# sun_sc_juice = juice_frame.convert_to_nsc(sun_jse_juice)
# earth_sc_juice = juice_frame.convert_to_nsc(earth_jse_juice)
# io_sc_juice = juice_frame.convert_to_nsc(io_jse_juice)
# europa_sc_juice = juice_frame.convert_to_nsc(europa_jse_juice)
# ganymede_sc_juice = juice_frame.convert_to_nsc(ganymede_jse_juice)
# callisto_sc_juice = juice_frame.convert_to_nsc(callisto_jse_juice)
fig = plt.figure(figsize=(12, 12))
axes={}
### Orbit plot (xy)
ax = fig.add_axes([0.1, 0.4, 0.35, 0.35])
axes['orbit_jse_xy'] = ax
# ax = patt.axis_configuration_xy(juice_frame, axis=ax)
# Spacecraft position in JSE
ax.plot([juice_jse_jupiter[0]], [juice_jse_jupiter[1]], 'ko')
# SC->JSE conversion, plot
armlen = 3
vec_xaxis_jse = juice_frame.convert_to_ref([1, 0 ,0])
vec_yaxis_jse = juice_frame.convert_to_ref([0, 1 ,0])
vec_zaxis_jse = juice_frame.convert_to_ref([0, 0 ,1])
ax.plot(juice_jse_jupiter[0] + [0, vec_xaxis_jse[0] * armlen],
juice_jse_jupiter[1] + [0, vec_xaxis_jse[1] * armlen],
'r')
ax.text(juice_jse_jupiter[0] + vec_xaxis_jse[0] * armlen,
juice_jse_jupiter[1] + vec_xaxis_jse[1] * armlen,
'x', color='r')
ax.plot(juice_jse_jupiter[0] + [0, vec_yaxis_jse[0] * armlen],
juice_jse_jupiter[1] + [0, vec_yaxis_jse[1] * armlen],
'g')
ax.text(juice_jse_jupiter[0] + vec_yaxis_jse[0] * armlen,
juice_jse_jupiter[1] + vec_yaxis_jse[1] * armlen,
'y', color='g')
ax.plot(juice_jse_jupiter[0] + [0, vec_zaxis_jse[0] * armlen],
juice_jse_jupiter[1] + [0, vec_zaxis_jse[1] * armlen],
'b')
ax.text(juice_jse_jupiter[0] + vec_zaxis_jse[0] * armlen,
juice_jse_jupiter[1] + vec_zaxis_jse[1] * armlen,
'z', color='b')
# Decoration
ax.plot([0, 0], [-100, 100], 'k:')
ax.plot([-100, 100], [0, 0], 'k:')
ax.plot([-100, 100], [100, -100], 'k:')
ax.plot([-100, 100], [-100, 100], 'k:')
x = np.cos(np.linspace(0, np.pi * 2, 100))
y = np.sin(np.linspace(0, np.pi * 2, 100))
for r in np.arange(5, 100, 5): ax.plot(r * x, r * y, 'k:')
ax.plot(x, y, 'k')
ax.set_xlim(-20, 20)
ax.set_ylim(-20, 20)
ax.set_xlabel('X [JSE]')
ax.set_ylabel('Y [JSE]')
ax.set_title('{:%FT%T}'.format(t))
ax.set_aspect('equal')
ax.text(io_jse_jupiter[0], io_jse_jupiter[1], 'I', color='g', ha='center', va='center')
ax.text(europa_jse_jupiter[0], europa_jse_jupiter[1], 'E', color='g', ha='center', va='center')
ax.text(ganymede_jse_jupiter[0], ganymede_jse_jupiter[1], 'G', color='g', ha='center', va='center')
ax.text(callisto_jse_jupiter[0], callisto_jse_jupiter[1], 'C', color='g', ha='center', va='center')
# Sun,earth direction
ax.plot([0, sun_jse_jupiter_unitvector[0] * 10],
[0, sun_jse_jupiter_unitvector[1] * 10],
'r-')
ax.text(sun_jse_jupiter_unitvector[0]* 10, sun_jse_jupiter_unitvector[1] * 10,
'S', color='r')
ax.plot([0, earth_jse_jupiter_unitvector[0] * 10],
[0, earth_jse_jupiter_unitvector[1] * 10],
'b-')
ax.text(earth_jse_jupiter_unitvector[0] * 10, earth_jse_jupiter_unitvector[1] * 10,
'E', color='b')
### Orbit plot (xz)
ax = fig.add_axes([0.1, 0.78, 0.35, 0.20])
axes['orbit_jse_xz'] = ax
# ax = patt.axis_configuration_xz(juice_frame, axis=ax)
# Spacecraft position in JSE
ax.plot([juice_jse_jupiter[0]], [juice_jse_jupiter[2]], 'ko')
# SC->JSE conversion, plot
ax.plot(juice_jse_jupiter[0] + [0, vec_xaxis_jse[0] * armlen],
juice_jse_jupiter[2] + [0, vec_xaxis_jse[2] * armlen],
'r')
ax.text(juice_jse_jupiter[0] + vec_xaxis_jse[0] * armlen,
juice_jse_jupiter[2] + vec_xaxis_jse[2] * armlen,
'x', color='r')
ax.plot(juice_jse_jupiter[0] + [0, vec_yaxis_jse[0] * armlen],
juice_jse_jupiter[2] + [0, vec_yaxis_jse[2] * armlen],
'g')
ax.text(juice_jse_jupiter[0] + vec_yaxis_jse[0] * armlen,
juice_jse_jupiter[2] + vec_yaxis_jse[2] * armlen,
'y', color='g')
ax.plot(juice_jse_jupiter[0] + [0, vec_zaxis_jse[0] * armlen],
juice_jse_jupiter[2] + [0, vec_zaxis_jse[2] * armlen],
'b')
ax.text(juice_jse_jupiter[0] + vec_zaxis_jse[0] * armlen,
juice_jse_jupiter[2] + vec_zaxis_jse[2] * armlen,
'z', color='b')
# Decoration
ax.plot([0, 0], [-1.1, 1.1], 'k-')
ax.plot([-100, 100], [0, 0], 'k:')
x = np.cos(np.linspace(0, np.pi * 2, 100))
y = np.sin(np.linspace(0, np.pi * 2, 100))
ax.plot(x, y, 'k')
ax.set_xlim(-20, 20)
ax.set_ylim(-6, 6)
ax.set_xlabel('X [JSE]')
ax.set_ylabel('Z [JSE]')
ax.set_title('{:%FT%T}'.format(t))
ax.set_aspect('equal')
ax.text(io_jse_jupiter[0], io_jse_jupiter[2], 'I', color='g', ha='center', va='center')
ax.text(europa_jse_jupiter[0], europa_jse_jupiter[2], 'E', color='g', ha='center', va='center')
ax.text(ganymede_jse_jupiter[0], ganymede_jse_jupiter[2], 'G', color='g', ha='center', va='center')
ax.text(callisto_jse_jupiter[0], callisto_jse_jupiter[2], 'C', color='g', ha='center', va='center')
### Sun direction
ax.plot([0, sun_jse_jupiter_unitvector[0] * 10],
[0, sun_jse_jupiter_unitvector[2] * 10],
'r-')
ax.text(sun_jse_jupiter_unitvector[0] * 19, sun_jse_jupiter_unitvector[2] * 10,
'S', color='r')
ax.plot([0, earth_jse_jupiter_unitvector[0] * 10],
[0, earth_jse_jupiter_unitvector[2] * 10],
'b-')
ax.text(earth_jse_jupiter_unitvector[0] * 10, earth_jse_jupiter_unitvector[2] * 10,
'E', color='b')
### S/C at Jupiter's System III frame
### Orbit plot (xy)
ax = fig.add_axes([0.5, 0.4, 0.35, 0.35])
axes['or bit_sys3_xy'] = ax
ax.plot([juice_sys3_jupiter[0]], [juice_sys3_jupiter[1]], 'ko')
ax.plot([0, 0], [-100, 100], 'k:')
ax.plot([-100, 100], [0, 0], 'k:')
ax.plot([-100, 100], [100, -100], 'k:')
ax.plot([-100, 100], [-100, 100], 'k:')
x = np.cos(np.linspace(0, np.pi * 2, 100))
y = np.sin(np.linspace(0, np.pi * 2, 100))
for r in np.arange(5, 100, 5): ax.plot(r * x, r * y, 'k:')
ax.plot(x, y, 'k')
ax.set_xlim(-20, 20)
ax.set_ylim(-20, 20)
ax.set_xlabel('X [SYS3]')
ax.set_ylabel('Y [SYS3]')
ax.set_title('{:%FT%T}'.format(t))
ax.text(io_sys3_jupiter[0], io_sys3_jupiter[1], 'I', color='g', ha='center', va='center')
ax.text(europa_sys3_jupiter[0], europa_sys3_jupiter[1], 'E', color='g', ha='center', va='center')
ax.text(ganymede_sys3_jupiter[0], ganymede_sys3_jupiter[1], 'G', color='g', ha='center', va='center')
ax.text(callisto_sys3_jupiter[0], callisto_sys3_jupiter[1], 'C', color='g', ha='center', va='center')
# Sun,earth direction
ax.plot([0, sun_sys3_jupiter_unitvector[0] * 10],
[0, sun_sys3_jupiter_unitvector[1] * 10],
'r-')
ax.text(sun_sys3_jupiter_unitvector[0] * 10, sun_sys3_jupiter_unitvector[1] * 10,
'S', color='r')
ax.plot([0, earth_sys3_jupiter_unitvector[0] * 10],
[0, earth_sys3_jupiter_unitvector[1] * 10],
'b-')
ax.text(earth_sys3_jupiter_unitvector[0] * 10, earth_sys3_jupiter_unitvector[1] * 10,
'E', color='b')
### Plot of SC in x-z plane.
ax = fig.add_axes([0.5, 0.78, 0.35, 0.20])
axes['orbit_sys3_xz'] = ax
ax.plot([juice_sys3_jupiter[0]], [juice_sys3_jupiter[2]], 'ko')
ax.plot([0, 0], [-1.1, 1.1], 'k-')
ax.plot([-100, 100], [0, 0], 'k:')
x = np.cos(np.linspace(0, np.pi * 2, 100))
y = np.sin(np.linspace(0, np.pi * 2, 100))
ax.plot(x, y, 'k')
ax.set_xlim(-20, 20)
ax.set_ylim(-6, 6)
ax.set_xlabel('X [SYS3]')
ax.set_ylabel('Z [SYS3]')
ax.set_aspect('equal')
### Other objects
ax.text(io_sys3_jupiter[0], io_sys3_jupiter[2], 'I', color='g', ha='center', va='center')
ax.text(europa_sys3_jupiter[0], europa_sys3_jupiter[2], 'E', color='g', ha='center', va='center')
ax.text(ganymede_sys3_jupiter[0], ganymede_sys3_jupiter[2], 'G', color='g', ha='center', va='center')
ax.text(callisto_sys3_jupiter[0], callisto_sys3_jupiter[2], 'C', color='g', ha='center', va='center')
# Sun,earth direction
ax.plot([0, sun_sys3_jupiter_unitvector[0] * 10],
[0, sun_sys3_jupiter_unitvector[2] * 10],
'r-')
ax.text(sun_sys3_jupiter_unitvector[0] * 10, sun_sys3_jupiter_unitvector[2] * 10,
'S', color='r')
ax.plot([0, earth_sys3_jupiter_unitvector[0] * 10],
[0, earth_sys3_jupiter_unitvector[2] * 10],
'b-')
ax.text(earth_sys3_jupiter_unitvector[0] * 10, earth_sys3_jupiter_unitvector[2] * 10,
'E', color='b')
### Calculate the corotation flow.
# Assume the corotation flow is axisymmetric flow wrt the rotation axis (System 3)
# Position of the spacecraft in sys3 is
# The flow vector is perpendicular to the position vector in x-y plane, without z component.
speed_solid_corotation_sys3 = length(juice_sys3_jupiter) * rj * 2 * np.pi / 36000. # 10 hour per rotation assumed
vel_solid_corotation_sys3 = [-juice_sys3_jupiter[1], juice_sys3_jupiter[0], 0]
vel_solid_corotation_sys3_unitvector = unit_vector(vel_solid_corotation_sys3)
vel_solid_corotation_sys3 = vel_solid_corotation_sys3_unitvector * speed_solid_corotation_sys3
print('Corotation/1st order', vel_solid_corotation_sys3, speed_solid_corotation_sys3)
### Geometric conversion from SystemIII to JSE. (Velocity measured at spacecraft -> geometric/position conversion,
### not velocity conversion)
vel_solid_corotation_jse = cmat_sys3_jse.dot(vel_solid_corotation_sys3)
### Plasma flow at the spacecraft. Expressed in JSE frame.
vel_solid_corotation_sys3 = vel_solid_corotation_jse - v_juice_jse_jupiter
### Plasma flow at the spacecraft. Convert it to SC frame.
vel_solid_corotation_sc = juice_frame.convert_to_nsc(vel_solid_corotation_sys3)
print('{0:%FT%T} {4:.5f} {1:.3f} {2[0]:.2f} {2[1]:.2f} {2[2]:.2f} {3[0]:.2f} {3[1]:.2f} {3[2]:.2f}'.format(
t, et, vel_solid_corotation_sys3, vel_solid_corotation_sc, mpldates.date2num(t)), file=fp)
print('{0:%FT%T} {4:.5f} {1:.3f} {2[0]:.2f} {2[1]:.2f} {2[2]:.2f} {3[0]:.2f} {3[1]:.2f} {3[2]:.2f}'.format(
t, et, vel_solid_corotation_sys3, vel_solid_corotation_sc, mpldates.date2num(t)))
### Viewing direction axis
ax = fig.add_axes([0.1, 0.08, 0.35, 0.20])
## Boundary of the viewing direction
import fovplot
fovplot.frame_plot_limit(ax=ax, alpha=0.1)
ax.text(0, 0, 'Tilt=0', color='r')
fovplot.frame_plot_limit(ax=ax, tilt=15, alpha=0.1)
ax.text(0, 15, 'Tilt=15', color='r')
fovplot.frame_plot_limit(ax=ax, tilt=30, alpha=0.1)
ax.text(0, 30, 'Tilt=30', color='r')
fovplot.frame_plot_limit(ax=ax, tilt=45, alpha=0.1)
ax.text(0, 45, 'Tilt=45', color='r')
## Direction of the corotation flow in SC
looking_vco_long = longitude(-vel_solid_corotation_sc)
looking_vco_lat = latitude(-vel_solid_corotation_sc)
ax.plot([looking_vco_long], [looking_vco_lat], 'bo')
ax.text(looking_vco_long, looking_vco_lat, 'Corot', color='b')
## Galileans / Earth / Sun
for num, nam, clr in zip((n_jp, n_sun, n_earth, n_io, n_europa, n_ganymede, n_callisto),
('Jupiter', 'Sun', 'Earth', 'Io', 'Europa', 'Ganymede', 'Callisto'),
('g', 'r', 'b', 'g', 'g', 'g', 'g')):
looking_object_jse, __vel = get_state(num, et, 'JSE', -907)
looking_object_sc = juice_frame.convert_to_nsc(looking_object_jse)
looking_object_long = longitude(looking_object_sc)
looking_object_lat = latitude(looking_object_sc)
if np.abs(looking_object_lat) > 89.5:
looking_object_long = 0.
ax.plot([looking_object_long], [looking_object_lat], clr + 'o')
ax.text(looking_object_long, looking_object_lat, nam, color=clr)
print(nam, looking_object_jse, looking_object_sc, looking_object_long, looking_object_lat)
plt.gcf().savefig(os.path.join(pathname, 'scenario1_{:d}.png'.format(int(et))))
plt.close('all')
fp.flush()
if __name__ == '__main__':
main()