''' Calisto 30 flyby on Galileo summary plot.
'''
import logging
logging.basicConfig()
logging.getLogger().setLevel(logging.INFO)
import datetime
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, Wedge
import matplotlib.cm
#import pyana.pep.pep_spice
import pyana.pep.galileo_spice
import pyana.pep.pep_spice
import pyana.pep.moon_image
def time_arange(t0, t1, tr, append=None):
tlist = []
t = t0
while t <= t1:
tlist.append(t)
t = t + tr
if append:
for t in append:
tlist.append(t)
tlist = list(set(tlist))
tlist.sort()
return tlist
def fov_at_surface(fov_at_sc, height, radius):
''' Return the FOV at the surface (half angle in degrees both for in- and output)
'''
fov_radian = fov_at_sc * np.pi / 180.
return (np.arcsin((height+radius) / radius * np.sin(fov_radian)) - fov_radian) * 180. / np.pi
def main():
rj = 71492.
# Center time from AAREADME.TXT in GALILEO PDS.
tc = datetime.datetime(2000, 5, 20, 10, 10, 10)
# Time info for flyby period
t0f = tc - datetime.timedelta(minutes=30)
t1f = tc + datetime.timedelta(minutes=30)
# trf = datetime.timedelta(seconds=4)
trf = datetime.timedelta(seconds=10)
tlistf = time_arange(t0f, t1f, trf, append=[tc])
# Time info for browse period
t0b = tc - datetime.timedelta(hours=3)
t1b = tc + datetime.timedelta(hours=3)
trb = datetime.timedelta(minutes=30)
tlistb = time_arange(t0b, t1b, trb, append=[tc])
# Time info for global period
t0g = tc - datetime.timedelta(days=3)
t1g = tc + datetime.timedelta(days=3)
trg = datetime.timedelta(minutes=30)
tlistg = time_arange(t0g, t1g, trg, append=[tc])
# Instance orbit information classes
gal = pyana.pep.galileo_spice.GalileoSpice.get_default_instance()
pepsp = pyana.pep.pep_spice.get_default_pepspice()
# Instnce figure information
fig = plt.figure(figsize=(18, 12))
# Position of Callisto/Galileo relative to Jupiter
gnm_jjse = []
gal_jjse = gal.get_positions(tlistg) / rj
for t in tlistg:
pos = gal.get_ganymede_position(t)
gnm_jjse.append(pos / rj)
gnm_jjse = np.array(gnm_jjse)
ax1 = fig.add_subplot(221)
ax1.plot(gnm_jjse[:, 0], gnm_jjse[:, 1], 'k', label='Ganymede')
ax1.plot(gal_jjse[:, 0], gal_jjse[:, 1], 'k:', label='Galileo')
# Positions at epoch
galileo_jjse_tc = gal.get_position(tc) / rj
ganymede_jjse_tc = gal.get_ganymede_position(tc) / rj
ax1.plot(galileo_jjse_tc[0], galileo_jjse_tc[1], 'ro')
ax1.plot(ganymede_jjse_tc[0], ganymede_jjse_tc[1], 'bo')
print('Galileo (TC)', galileo_jjse_tc)
print('ganymede (TC)', ganymede_jjse_tc)
print('G-G (TC)', galileo_jjse_tc - ganymede_jjse_tc)
galileo_jjse_t0g = gal.get_position(t0g) / rj
ganymede_jjse_t0g = gal.get_ganymede_position(t0g) / rj
galileo_jjse_t1g = gal.get_position(t1g) / rj
ganymede_jjse_t1g = gal.get_ganymede_position(t1g) / rj
ax1.text(galileo_jjse_t0g[0], galileo_jjse_t0g[1], 't0')
ax1.text(galileo_jjse_t1g[0], galileo_jjse_t1g[1], 't1')
ax1.text(ganymede_jjse_t0g[0], ganymede_jjse_t0g[1], 't0')
ax1.text(ganymede_jjse_t1g[0], ganymede_jjse_t1g[1], 't1')
# Display Jupiter
jupiter_body = Circle((0, 0), 1, edgecolor='k', fill=False, alpha=0.5)
ax1.add_patch(jupiter_body)
jupiter_dark = Wedge((0, 0), 1, 90., 270., color='k', alpha=0.5)
ax1.add_patch(jupiter_dark)
ax1.set_aspect('equal')
ax1.set_xlim(ax1.get_xlim()[::-1])
ax1.set_ylim(ax1.get_ylim()[::-1])
ax1.set_title('%s-%s' % (t0g, t1g))
ax1.set_xlabel('Xjse [Rj]; Jupiter-center')
ax1.set_ylabel('Yjse [Rj]; Jupiter-center')
ax1.legend()
# Position of galileo in calisto-centered JSE
galileo_posb = []
for t in tlistb:
pos = gal.get_ganymede_position(t)
ganymede_pos = pos / rj
galileo_pos = gal.get_position(t) / rj
galileo_pos = galileo_pos - ganymede_pos
galileo_posb.append(galileo_pos)
galileo_posb = np.array(galileo_posb)
galileo_posf = []
for t in tlistf:
pos = gal.get_ganymede_position(t)
ganymede_pos = pos / rj
galileo_pos = gal.get_position(t) / rj
galileo_pos = galileo_pos - ganymede_pos
galileo_posf.append(galileo_pos)
galileo_posf = np.array(galileo_posf)
# Plot.
ax2 = fig.add_subplot(222)
ax2.plot(galileo_posb[:, 0] * rj, galileo_posb[:, 1] * rj, 'k.:')
# ax2.plot(galileo_posf[:, 0], galileo_posf[:, 1], 'k.')
rgnm = 2634.1 / rj
ganymede_body = Circle((0, 0), rgnm * rj, edgecolor='k', fill=False, alpha=0.5)
ax2.add_patch(ganymede_body)
ganymede_dark = Wedge((0, 0), rgnm * rj, 90., 270., color='k', alpha=0.5)
ax2.add_patch(ganymede_dark)
ax2.set_xlabel('JSE_X (Callisto-center) [km]')
ax2.set_ylabel('JSE_Y (Callisto-center) [km]')
ax2.set_title('Galileo trajectory. CA=%s' % tc)
ax2.set_aspect('equal')
ax2.set_xlim(-rj, rj)
ax2.set_ylim(-rj, rj)
ax2.set_xlim(ax2.get_xlim()[::-1])
ax2.set_ylim(ax2.get_ylim()[::-1])
# Direction of Jupiter
jupiter_dir = pepsp.get_posvel('Jupiter', tc, 'JSE', 'Ganymede')[0]
# Plot of the Jupiter direction and velocity vector
xlim = ax2.get_xlim()
xlen = np.abs(xlim[1] - xlim[0])
arrow_len = xlen * 0.1
logging.info('Length of arrow %f' % arrow_len)
# Jupiter direction plot
jupiter_dir = jupiter_dir / np.sqrt((jupiter_dir ** 2).sum()) * arrow_len
logging.info('Juptier dir = %f, %f' % (jupiter_dir[0], jupiter_dir[1]))
ax2.arrow(0, 0, jupiter_dir[0], jupiter_dir[1], color='r')
ax2.text(jupiter_dir[0], jupiter_dir[1], 'J', color='r')
# Plot of height profile
ax3 = fig.add_subplot(223)
# x axis is time, tlistf
# y asis is height. Calculate from galileo_posf.
height = (np.sqrt((galileo_posf ** 2).sum(1)))
ax3.plot_date(tlistf, (height - rgnm) * rj, 'r.')
ax3.set_ylabel('Height of S/C from Ganymede [km]')
ax3.set_ylim(0, ax3.get_ylim()[1])
# y2 axis is Horizontal velocity relative to Ganymede.
ax3r = ax3.twinx()
vel_hol = []
for tt in tlistf:
vel_vec = gal.get_velocity(tt, 'Ganymede', 'IAU_GANYMEDE')
vtot2 = (vel_vec ** 2).sum()
pos = gal.get_position(tt, 'Ganymede', 'IAU_GANYMEDE')
dist = np.sqrt((pos ** 2).sum())
pos = pos / dist
rv = (vel_vec * pos).sum()
vhol = np.sqrt(vtot2 - rv ** 2)
vel_hol.append(vhol)
vel_hol = np.array(vel_hol)
ax3r.plot_date(tlistf, vel_hol, 'g.')
ax3r.set_ylim(0, ax3r.get_ylim()[1])
ax3r.set_ylabel('Relative horizontal velocity [km/s]')
ax3.set_xlabel('UT')
ax3.set_title('CA = %.1f km' % ((min(height) - rgnm) * rj))
# Plot over the map
fov_full = 5.0 # Full width
fov = fov_full / 2.0
# Ganymede map
map = pyana.pep.moon_image.get_ganymede_image()
ax4 = fig.add_subplot(224)
ax4.imshow(map, origin='bootom', extent=[-180, 180, -90, 90], aspect='equal', cmap=matplotlib.cm.gray)
hlist = []
lonlist = []
latlist = []
galileo_posf = []
fovlist = []
for t in tlistf:
pos, vel = pepsp.get_posvel('Ganymede', t, 'JSE', 'Jupiter')
ganymede_pos = pos / rj
galileo_pos = gal.get_position(t) / rj
galileo_pos = galileo_pos - ganymede_pos # Galileo in G-JSE
# JSE to IAU_CALISTO
galileo_iaugnm = pepsp.convert_vector(galileo_pos, 'JSE', 'IAU_GANYMEDE', t)
# print galileo_iaugnm
# H, Lon, Lat
h = np.sqrt((galileo_iaugnm ** 2).sum())
lon = np.arctan2(galileo_iaugnm[1], galileo_iaugnm[0]) * 180. / np.pi
lat = 90 - np.arccos(galileo_iaugnm[2] / h) * 180. / np.pi
# FoV mapping
# theta_radius = (np.arcsin((rgnm + h) / rcgnm* np.sin(fov)) - fov) * 180. / np.pi
theta_radius = fov_at_surface(fov, h - rgnm, rgnm)
print(theta_radius, fov, h * rj, rgnm * rj)
if np.isfinite(theta_radius):
hlist.append(h)
lonlist.append(lon)
latlist.append(lat)
fovlist.append(theta_radius)
# print lonlist, latlist
# ax4.plot(lonlist, latlist, 'r.')
for lon, lat, r in zip(lonlist, latlist, fovlist):
circ = Circle((lon, lat), radius=r, alpha=0.4, fc='r', color='k')
ax4.add_patch(circ)
ax4.set_xlim(-180, 180)
ax4.set_ylim(-90, 90)
ax4.set_xlabel('East longitude')
ax4.set_ylabel('Latitude')
ax4.set_title('FoV: Rough sketch: fov=%.1f %s-%s dt=%.1fs' % (fov_full, t0f.strftime('%T'), t1f.strftime('%T'), trf.seconds))
fig.savefig('plot_ga28_%.1f.png' % (fov_full), dpi=180)
if __name__ == '__main__':
val = main()