''' JUICE orbit around ganymede.
Making three figures.
- Orbit around Jupiter (2D).
- Orbit around ganymde (3D).
- Map above the 2D surface.
'''
import os
import sys
import logging
logging.basicConfig()
import datetime
import dateutil.parser
import math
from optparse import OptionParser
import matplotlib.pyplot as plt
import matplotlib.dates
import matplotlib.patches
import numpy as np
import scipy as sp
from mpl_toolkits.mplot3d import Axes3D
import spice
import pyana.pep.juice_spice
import pyana.pep.moon_image
import pyana.util.gridsphere
import pyana.util.cone
timemaskfunction = lambda x: True
def setup_kernel(kernels=None, timemaskfunc=None):
global timemaskfunction
if timemaskfunc != None:
timemaskfunction = timemaskfunc
if kernels == None:
kernels = [os.path.join('..', 'pyana', 'pep', 'spice_juice', 'mantra.jgo_2020_001_ipc_cal_pso_res_e40_562_200.bsp')]
timemaskfunction = pyana.pep.juice_spice.is_valid_interval_mantra001
kernsdir = os.path.join('..', 'pyana', 'pep', 'spice_juice')
kerns = ['naif0010.tls', 'pck00009.tpc', 'de405.bsp', 'jup230.bsp', ]
for k in kerns:
kern = os.path.join(kernsdir, k)
print(kern)
spice.furnsh(kern)
spice.furnsh(os.path.join('..', 'pyana', 'pep', 'spice', 'jse_111130.tf'))
for k in kernels:
print(k)
spice.furnsh(k)
def ganymede_around_jupiter(t0=datetime.datetime(2028, 11, 28, 0, 0, 0), t1=datetime.datetime(2028, 11, 29, 0, 0, 0),
dt=datetime.timedelta(minutes=4)):
frame = "JSE"
center = "JUPITER"
center_id = spice.bodn2c(center)
moon = 'Ganymede'
moon_id = spice.bodn2c(moon)
# print moon, moon_id
spacecraft = -999
t = t0
rj = 71492.
ganymedepos = []
spacecraftpos = []
while t <= t1:
if not timemaskfunction(t):
t = t + dt
continue
et = spice.str2et(t.strftime('%FT%T'))
posvel = spice.spkez(spacecraft, et, frame, 'LT+S', center_id)[0]
spacecraftpos.append(posvel[:3])
posvel = spice.spkez(moon_id, et, frame, 'LT+S', center_id)[0]
ganymedepos.append(posvel[:3])
# print t, posvel
t = t + dt
ganymedepos = np.array(ganymedepos) / rj
spacecraftpos = np.array(spacecraftpos) / rj
jupiter_body = matplotlib.patches.Circle((0, 0), 1, edgecolor='none')
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(ganymedepos[:, 0], ganymedepos[:, 1], '-')
ax.plot(spacecraftpos[:, 0], spacecraftpos[:, 1], '-')
ax.add_patch(jupiter_body)
ax.set_aspect('equal')
ax.set_xlim(20, -20)
ax.set_ylim(20, -20)
ax.grid()
ax.set_xlabel('X(JSE) [km]')
ax.set_ylabel('Y(JSE) [km]')
return fig
def juice_around_ganymede(t0=datetime.datetime(2028, 11, 28, 0, 0, 0), t1=datetime.datetime(2028, 11, 29, 0, 0, 0),
dt=datetime.timedelta(minutes=10)):
frame = 'IAU_GANYMEDE'
center = 'GANYMEDE'
center_id = spice.bodn2c(center)
spacecraft = -999
t = t0
rg = 2634.1
spacecraftpos = []
while t <= t1:
if not timemaskfunction(t):
t = t + dt
continue
et = spice.str2et(t.strftime('%FT%T'))
posvel = spice.spkez(spacecraft, et, frame, 'LT+S', center_id)[0]
spacecraftpos.append(posvel[:3])
t = t + dt
spacecraftpos = np.array(spacecraftpos) / rg
# print spacecraftpos
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x = np.outer(np.cos(u), np.sin(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.ones(np.size(u)), np.cos(v))
ax.plot_surface(x, y, z, rstride=4, cstride=4, color='b')
ax.plot(spacecraftpos[:, 0], spacecraftpos[:, 1], spacecraftpos[:, 2], color='r')
ax.plot(spacecraftpos[:, 0], spacecraftpos[:, 1], spacecraftpos[:, 2], zdir='z')
ax.set_aspect('equal')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
return ax
def visibility_append(gridsphere, xyz_sc, half_fov_deg):
xyz = np.array(xyz_sc)
r = np.linalg.norm(xyz)
xyz_norm = xyz / r
rg = 2634.1
# horizon_angle = np.arcsin((rg + r) / rg * np.sin(half_fov_deg * np.pi / 180.)) - half_fov_deg * np.pi / 180.
# cos_hor_angle = np.cos(horizon_angle)
cos_hor_angle = rg / r
horizon_angle = np.arccos(cos_hor_angle)
n = len(gridsphere)
cosfov = np.cos(half_fov_deg * np.pi / 180.0)
for index in range(n):
spherepos_n = gridsphere.index2xyz(index) # Postion vector at surface, norm.
# Invisible surface, ignore.
if np.dot(spherepos_n, xyz_norm) < cos_hor_angle:
continue
spherepos = spherepos_n * rg
# Sphere position vector (normalized) seen from spacecraft.
spherepos_from_sc = spherepos - xyz
spherepos_from_sc = spherepos_from_sc / np.linalg.norm(spherepos_from_sc)
cos_spherepos = np.dot(spherepos_from_sc, -xyz_norm)
if cos_spherepos >= cosfov:
gridsphere.append_value(index, 1.0)
# print xyz_norm, spherepos_n
def addfovcircle(ax, tlist, fov):
'''
'''
frame = 'IAU_GANYMEDE'
center = 'GANYMEDE'
center_id = spice.bodn2c(center)
spacecraft = -999
rg = 2634.1
# print tlist
for t in tlist:
print('Prcessing', t)
et = spice.str2et(t.strftime('%FT%T'))
# Spacecraft position.
posvel = spice.spkez(spacecraft, et, frame, 'LT+S', center_id)[0]
x, y, z, vx, vy, vz = posvel
# Spacecraft distance
r = np.sqrt(x ** 2 + y ** 2 + z ** 2)
# Horizon angle seen from spacecraft
horizon_angle = np.arcsin(rg / r)
# print 'S/C distance=', r
#
# print 'Horizon angle =', np.degrees(horizon_angle), 90 - np.degrees(horizon_angle)
if fov >= np.degrees(horizon_angle):
# Horizon is the visible angle from the s/c
# In this case, the "intersection" is defined by the cone
# centered at the Ganymede with direction of spacecraft as center
# with the angle of 90-horizon_angle.
visible_angles = 90 - np.degrees(horizon_angle)
else:
# Fov is the visible angle from the s/c
# Converting to the visible angle of the sphere.
visible_angles = np.degrees(np.arcsin(r / rg * np.sin(np.radians(fov)))) - fov
# print 'Visible angle =', visible_angles
vecs = pyana.util.cone.get_surface_vectors((x, y, z), visible_angles)
# (3, 180), normalized.
lats = 90 - np.degrees(np.arccos(vecs[2, :]))
lons = np.degrees(np.arctan2(vecs[1, :], vecs[0, :]))
while lons.max() >= 360:
lons = np.where(lons >= 360, lons - 360, lons)
while lons.min() < 0:
lons = np.where(lons < 0, lons + 360, lons)
lons = np.degrees(np.unwrap(np.radians(lons)))
ax.plot(lons, lats, 'r-', alpha=0.7)
ax.plot(lons - 360, lats, 'r-', alpha=0.7)
ax.plot(lons + 360, lats, 'r-', alpha=0.7)
return lons, lats
def juice_on_ganymede(t0 = datetime.datetime(2028, 11, 28, 0, 0, 0), t1=datetime.datetime(2028, 11, 28, 12, 0, 0),
dt = datetime.timedelta(seconds=600), fov=None, fov_circle_time=None, fill='counts'):
'''
:keyword fov: Field of view angle (half angle)
:keyword fov_circle_time: If "start", the start time FoV is shown as circle.
If "stop", the stop time FoV is shown as circle.
If a sequence of ``datetime.datetime`` objects is given, the FoVs at the specific times are shown.
:keyword fill: Fill type. 'counts' or 'coverage'
Note that if ``fov`` is not given (or given as ``None``), ``fov_circle_time`` is neglected.
'''
frame = 'IAU_GANYMEDE'
center = 'GANYMEDE'
center_id = spice.bodn2c(center)
spacecraft = -999
t = t0
spacecraftpos = []
sphere = pyana.util.gridsphere.SimpleGridSphere()
while t <= t1:
if not timemaskfunction(t):
t = t + dt
continue
et = spice.str2et(t.strftime('%FT%T'))
posvel = spice.spkez(spacecraft, et, frame, 'LT+S', center_id)[0]
x, y, z, vx, vy, vz = posvel
r = np.sqrt(x ** 2 + y ** 2 + z ** 2)
vel = np.sqrt(vx**2 + vy ** 2 + vz ** 2)
nx = x/r
ny = y/r
nz = z/r
lat = np.arcsin(nz) * 180 / np.pi
lon = np.arctan2(ny, nx) * 180 / np.pi
if lat < 60:
print('Lat = %f, skip' % lat)
t = t + dt
continue
print('Processing at %s (lon=%.1f, lat=%.1f) d=%.3f vel=%.3f' % (str(t), lon, lat, r, vel))
spacecraftpos.append([lon, lat])
if fov:
visibility_append(sphere, (x, y, z), fov)
t = t + dt
spacecraftpos = np.array(spacecraftpos) # N x 2 array
fig = plt.figure()
ax = fig.add_subplot(111)
gmap = pyana.pep.moon_image.get_ganymede_image()
ax.imshow(gmap, origin='bootom', extent=[-180, 180, -90, 90], aspect='equal', cmap=matplotlib.cm.gray)
ax.imshow(gmap, origin='bootom', extent=[180, 540, -90, 90], aspect='equal', cmap=matplotlib.cm.gray)
if fov:
lons, lats = sphere.get_bgrid()
if fill == 'coverage':
cnts = sphere.get_average() # Average is 1 where data is there.
else:
cnts = np.ma.masked_less_equal(0, sphere.get_counts())
ax.pcolor(lons, lats, cnts, alpha=0.4, edgecolor='none')
if fov_circle_time != None:
if fov_circle_time == "start":
addfovcircle(ax, [t0], fov)
elif fov_circle_time == "stop":
addfovcircle(ax, [t1], fov)
elif fov_circle_time == 'all':
tlist = []
t = t0
while t <= t1:
tlist.append(t)
t = t + dt
addfovcircle(ax, tlist, fov)
else:
addfovcircle(ax, fov_circle_time, fov)
spacecraftpos[:, 0] = np.where(spacecraftpos[:, 0] < 0, spacecraftpos[:, 0] + 360, spacecraftpos[:, 0])
ax.plot(spacecraftpos[:, 0], spacecraftpos[:, 1], 'k+', markersize=0.5, alpha=0.5)
ax.set_xlim(0, 360)
ax.set_ylim(-90, 90)
return fig, ax
def main():
usage = "usage: %prog [options]"
parser = OptionParser(usage)
# parser.add_option("-k", "--kernel", dest="kernels", default=[],
# help="Load non-default kernel(s) of FILENAME", action='append')
# parser.add_option("-s", dest="start_time", default='2025-01-01T00:00:00',
# help="Set start time")
# parser.add_option("-e", dest="stop_time", default='2030-01-01T00:00:00',
# help="Set stop time")
#
# parser.add_option("-d", dest="delta_time", default=3600,
# help="Set time resolution in seconds.", type='int')
#
# parser.add_option("-v", "--verbose",
# action="store_true", dest="verbose")
# parser.add_option("-q", "--quiet",
# action="store_false", dest="verbose")
#
(options, args) = parser.parse_args()
#
# if len(args) != 0:
# parser.error("incorrect number of arguments")
#
# if options.verbose:
# logging.getLogger().setLevel(logging.DEBUG)
#
# start_time = dateutil.parser.parse(options.start_time)
# stop_time = dateutil.parser.parse(options.stop_time)
# delta_time = datetime.timedelta(seconds=options.delta_time)
# logging.debug('Time: %s --(%s)--> %s' % (str(start_time),
# repr(delta_time), str(stop_time)))
#
# print 'K=', options.kernels
# makeall(kernels=options.kernels, t0=start_time, t1=stop_time, dt=delta_time)
setup_kernel()
fov = 2.5
# prmsetname = 'p1'
# t0 = datetime.datetime(2028, 11, 28, 5, 0, 0)
# dt_duration_sec = 42000
# dt_sec = 240
# dt_fov_sec = 660
# linewidth = None
# title = '~2 Rm, FOV=%gdeg, dt=%gs (total %gs)' % (fov, dt_fov_sec, dt_duration_sec)
prmsetname = 'p2'
t0 = datetime.datetime(2029, 3, 1, 0, 30, 0)
dt_duration_sec = 10 * 3600
dt_sec = 20
dt_fov_sec = 15
linewidth = 0.5
title = '500 km, %g deg, dt=%gs (total %gs)' % (fov, dt_fov_sec, dt_duration_sec)
# prmsetname = 'p3'
# t0 = datetime.datetime(2029, 7, 1, 0, 0, 0)
# dt_duration_sec = 10000
# dt_sec = 10
# dt_fov_sec = 5
# title = '200 km, %g deg, dt=%gs (total %gs)' % (fov, dt_fov_sec, dt_duration_sec)
ndat = int(dt_duration_sec / dt_fov_sec)
dt_duration = datetime.timedelta(seconds=dt_duration_sec)
t1 = t0 + dt_duration
dt = datetime.timedelta(seconds=dt_sec)
dt_fov = datetime.timedelta(seconds=dt_fov_sec)
fig, ax = juice_on_ganymede(fov=fov, fov_circle_time=[t0, t0 + dt_fov, t0 + 2 * dt_fov,
t0 + dt_fov * (ndat / 2 - 1), t0 + dt_fov * (ndat / 2),
t0 + dt_fov * (ndat / 2 + 1)],
fill='coverage', t0=t0, t1=t1, dt=dt)
ax.set_title(title)
# fig.savefig('juice_ganymede_quicklook_%s.png' % (prmsetname, ), dpi=300)
if __name__ == '__main__':
main()