''' JUICE orbit around ganymede.
Making three figures.
- Orbit around Jupiter (2D).
- Orbit around ganymde (3D).
- Map above the 2D surface.
'''
import matplotlib
matplotlib.use('Agg')
debug = False
import random
runid = int(random.random() * 100000)
import os
import sys
import logging
logging.basicConfig()
logging.getLogger().setLevel(logging.DEBUG)
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 visibility_append(gridsphere, xyz_sc, half_fov_deg):
logger = logging.getLogger('visibility_append')
logger.setLevel(logging.DEBUG)
logger.debug(xyz_sc)
xyz = np.array(xyz_sc)
r = np.linalg.norm(xyz)
xyz_norm = xyz / r
lon = np.arctan2(xyz[1], xyz[0])
lat = np.arcsin(xyz_norm[2])
lon, lat = np.rad2deg((lon, lat))
rg = 2634.1
visible_angle_surface = np.arcsin(r / rg * np.sin(half_fov_deg * np.pi / 180.)) - half_fov_deg * np.pi / 180.
logger.debug('R=%f, FOV=%f, Vis angle = %f' % (r, half_fov_deg, np.rad2deg(visible_angle_surface)))
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)
footprint_index = gridsphere.lonlat2index(lon, lat)
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 or index == footprint_index:
gridsphere.append_value(index, np.rad2deg(visible_angle_surface))
# print xyz_norm, spherepos_n
# Examine really the footprint.
# If FOV is too small, the footprint is considered to be "visible"
juice = pyana.pep.juice_spice.JuiceSpice()
orbitnumber = pyana.pep.juice_spice.OrbitNumber001()
def map_on_grid_per_orbit(onr, fov, nlon=360, nlat=180):
logger = logging.getLogger('map_on_grid_per_orbit')
logger.setLevel(logging.DEBUG)
logger.debug('Onr= %d' % onr)
grid_sphere = pyana.util.gridsphere.SimpleGridSphere(nlon=nlon, nlat=nlat)
rg = 2634.1
# The grid data includes the "FoV angle" as a data.
t0 = orbitnumber.get_starttime(onr)
t1 = orbitnumber.get_stoptime(onr)
print(t0, t1, t1-t0)
dt = datetime.timedelta(seconds=60)
if debug:
dt = datetime.timedelta(seconds=3600)
t = t0
while t <= t1:
logger.debug('t= %s' % t)
if not timemaskfunction(t):
t = t + dt
continue
pos = juice.get_position(t, relative_to="GAnymede", frame="IAU_GANYMEDE")
h = np.linalg.norm(pos)
lat = np.arcsin(pos[2] / h)
lon = np.arctan2(pos[1], pos[0])
h -= rg
lat, lon = np.rad2deg((lat, lon)) # Spacecraft position
visibility_append(grid_sphere, pos, fov)
t = t + dt
return grid_sphere
def juice_map_on_grid(fov=2.5):
''' Make grid data
'''
nlon = 360
nlat = 180
if debug:
nlon = 180
nlat = 90
master_grid_sphere = pyana.util.gridsphere.SimpleGridSphere(nlon=nlon, nlat=nlat)
onr0 = 700
onr1 = 690
step = -1
if debug:
onr0 = 1700
onr1 = 1698
step = -1
for onr in range(onr0, onr1, step):
grid_sphere = map_on_grid_per_orbit(onr, fov, nlon=nlon, nlat=nlat)
# The grid_sphere includes the orbit data
cgrid = grid_sphere.get_cgrid()
val = grid_sphere.get_min() # want a minimum angular resolution
for ilat in range(nlat):
for ilon in range(nlon):
if not val.mask[ilat, ilon]:
print(cgrid[0][ilat, ilon], cgrid[1][ilat, ilon], val[ilat, ilon])
master_grid_sphere.append_value_lonlat(cgrid[0][ilat, ilon], cgrid[1][ilat, ilon], val[ilat, ilon])
# Just display the status.
statusfig = plt.figure()
statusax = statusfig.add_subplot(111)
statusbg = grid_sphere.get_bgrid()
statusvalue = grid_sphere.get_min()
statusim = statusax.pcolor(statusbg[0], statusbg[1], statusvalue, vmin=0)
statusax.set_title('Current orbit %d' % onr)
statusfig.colorbar(statusim)
statusfig.savefig('status1_%d.png' % runid)
plt.close(statusfig)
statusfig = plt.figure()
statusax = statusfig.add_subplot(111)
statusbg = master_grid_sphere.get_bgrid()
statusvalue = master_grid_sphere.get_min()
statusim = statusax.pcolor(statusbg[0], statusbg[1], statusvalue, vmin=0)
statusax.set_title('Minimum angle')
statusfig.colorbar(statusim)
statusfig.savefig('status2_%d.png' % runid)
plt.close(statusfig)
statusfig = plt.figure()
statusax = statusfig.add_subplot(111)
statusbg = master_grid_sphere.get_bgrid()
statusvalue = master_grid_sphere.get_counts()
statusim = statusax.pcolor(statusbg[0], statusbg[1], statusvalue, vmin=0)
statusax.set_title('Number obs')
statusfig.colorbar(statusim)
statusfig.savefig('status3_%d.png' % runid)
plt.close(statusfig)
fig = plt.figure()
ax = fig.add_subplot(111)
bgrid = master_grid_sphere.get_bgrid()
n_data = master_grid_sphere.get_counts()
min_resol = np.ma.masked_where(n_data <= 0, master_grid_sphere.get_min())
im = ax.pcolor(bgrid[0], bgrid[1], min_resol, vmin=0)
fig.colorbar(im)
fig2 = plt.figure()
ax2 = fig2.add_subplot(111)
n_data = np.ma.masked_where(n_data <=0, n_data)
im2 = ax2.pcolor(bgrid[0], bgrid[1], n_data, vmin=0)
fig2.colorbar(im2)
return (fig, ax, im), (fig2, ax2, im2), master_grid_sphere
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")
#
parser.add_option("-f", dest='fov', default="2.5")
(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 = 10. # Half FOV
fov = float(options.fov)
logging.info('FOV=%g' % fov)
# fov = 5. # Half FOV
# fov = 2.5 # Half FOV
title = 'Coverage over Ganymede mapping phase (Half-FOV=%g)' % fov
title2 = 'Minimum angle (Half-FOV=%g)' % fov
plt1, plt0, master_grid_sphere = juice_map_on_grid(fov=fov)
plt0[1].set_title(title)
plt1[1].set_title(title2)
plt0[0].savefig('juice_ganymede_coverage_%.1f.png' % (fov, ), dpi=300)
plt1[0].savefig('juice_ganymede_mininum_%.1f.png' % (fov, ), dpi=300)
# import cPickle as pickle
# fp = open('juice_ganymede_mininum_%.1f.png' % (fov,), 'w')
# pickle.dump(master_grid_sphere, fp)
# fp.close()
return master_grid_sphere
if __name__ == '__main__':
retval = main()