''' A main script to plot the FOV map, trajectory, and E-t (JDC) and C-t (JNA).
A related scripts are:
- :mod:`apps120718_jna_fov.traj_on_surface`.
- :mod:`apps120803_jdcfly.app08_jdc_jna_flyover`.
Note that the plasma coordinate and the geologic coordinate has 90 degrees difference.
For the final map, geologic coordinate is used, where 0 degree longitude is
Jupiter facing meridian. However, 0 degree longitude in plasma coordinate is
the facing the bodies orbital motion.
To convert the plasma coordinate to the geological coordinate for plotting,
90 degrees should be subtracted in the longitude.
.. image:: ../../../src/scripts/apps120816_flyover/flyover_mini.png
:width: 90%
'''
import datetime
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import pyana.pep.juice_spice1205 as jspice
import pyana.pep.pep_attitude
import pyana.jna.frame0 as frame
import pyana.jna.fov0 as jnafov
import pyana.jna.frame0 as jnaframe
import pyana.util.intersection
def pla2geo(vecs_pla):
vecs_geo = np.ma.zeros(vecs_pla.shape)
vecs_geo[0] = vecs_pla[1]
vecs_geo[1] = -vecs_pla[0]
vecs_geo[2] = vecs_pla[2]
return vecs_geo
def geo2pla(vecs_geo):
vecs_pla = np.ma.zeros(vecs_geo.shape)
vecs_pla[0] = -vecs_geo[1]
vecs_pla[1] = vecs_geo[0]
vecs_pla[2] = vecs_geo[2]
return vecs_pla
def vec2lonlat(vecs_geo):
lons = np.arctan2(vecs_geo[1], vecs_geo[0]) * 180. / np.pi
# lons is -180 to 180. Lons should range -90 to 270.
lons = np.array(lons)
lons = np.where(lons < -90, lons + 360, lons)
r = np.sqrt((vecs_geo ** 2).sum(0))
lats = np.arcsin(vecs_geo[2] / r) * 180. / np.pi
return lons, lats
def segmentize_lons(lons):
''' Return index boundary that separates the segment
Segment is separated by a boundary that has difference of
longitudal angle more than 180 degrees.
For example,
>>> lons = np.array([-175, -178, -179, 179, 177, 173, 178, 179, -178, -175])
>>> print(lons.shape)
(10,)
will be separated into three segments, namely first three negatives,
following five positives, and the last two negatives.
Returned is boundary index. In this example, [0, 3, 8, 10] will be returned.
The number of the segment is the length of the returned minus 1.
>>> idx = segmentize_lons(lons)
>>> print(idx)
[0, 3, 8, 10]
You can use the returned as
>>> nseg = len(idx) - 1
>>> segmented_lons = []
>>> for i in range(nseg):
... segmented_lons.append(lons[idx[i]:idx[i+1]])
>>> print(segmented_lons)
[array([-175, -178, -179]), array([179, 177, 173, 178, 179]), array([-178, -175])]
Obviously, the array as
>>> lons = np.arange(1000)
will return the boundary as [0, 1000], while it exceeds 360 degrees.
>>> print(segmentize_lons(lons))
[0, 1000]
'''
nlons = len(lons)
diff = np.abs(lons[:-1] - lons[1:])
segb = (np.where(diff > 180)[0] + 1).tolist()
segb.insert(0, 0)
segb.append(nlons)
return segb
def main(type="nadir"):
'''Main script'''
t0 = datetime.datetime(2033, 4, 4, 18, 50)
### Plot on map
fig = plt.figure(figsize=(24, 10))
ax = fig.add_subplot(111)
m = Basemap(projection='cyl', lon_0=90, ax=ax, llcrnrlat=25, llcrnrlon=40, urcrnrlat=70, urcrnrlon=140)
m.drawparallels(np.arange(-90, 120, 30), labels=[True, False, False, False])
m.drawmeridians(np.arange(-180, 210, 30), labels=[False, False, False, True])
m.warpimage('../../pyana/pep/images/ganymede.jpg')
### Map
import pyana.pep.mhddata
fmap = pyana.pep.mhddata.PrecipitationFlux(type=type)
lonlist = fmap.reshape(fmap.lonlist()) # In plasma frame
lonlist = lonlist - 90 # In geo-frame
latlist = fmap.reshape(fmap.latlist())
flxlist = fmap.reshape(fmap.flist())
lflxlist = np.ma.log10(flxlist)
x, y = m(lonlist, latlist)
# LOG scale
pc = m.contourf(x, y, lflxlist, np.linspace(7.9, 8.5, 7))
cb = fig.colorbar(pc)
cb.set_label('Log of flux [/cm2 s]')
from matplotlib import cm
pc2 = m.contour(x, y, lflxlist, np.linspace(7.5, 7.9, 5), cmap=cm.Blues)
plt.clabel(pc2)
### Open-close line
import pyana.pep.mhddata
lval = pyana.pep.mhddata.LValue1201()
ocb = lval.opencloseboundaries()
lon = np.arange(-90, 270.001, 2)
x, y = m(lon, ocb[0])
m.plot(x, y, 'y', lw=3)
x, y = m(lon, ocb[1])
m.plot(x, y, 'y', lw=3)
### Trajectory using spice
juice = jspice.JuiceSpice.get_default_instance()
### Track.
for tlen, style in zip((1080, 3600), ('w-', 'w--')):
track_time = t0 + datetime.timedelta(seconds=10) * np.arange(tlen) # 3 or 10 hours
poslist = juice.get_positions(track_time, relative_to='GANYMEDE', frame='IAU_GANYMEDE')
poslist = poslist.T
lon, lat = vec2lonlat(poslist)
segidx = segmentize_lons(lon)
seglon = [lon[segidx[i]:segidx[i+1]] for i in range(len(segidx) - 1)]
seglat = [lat[segidx[i]:segidx[i+1]] for i in range(len(segidx) - 1)]
for slon, slat in zip(seglon, seglat):
print(slon, slat)
x, y = m(slon, slat)
m.plot(x, y, style)
### Track every 10 minutes
tick_time = t0 + datetime.timedelta(minutes=10) * np.arange(60)
poslist = juice.get_positions(tick_time, relative_to='GANYMEDE', frame='IAU_GANYMEDE')
poslist = poslist.T
lon, lat = vec2lonlat(poslist)
x, y = m(lon, lat)
m.plot(x, y, 'wo')
### JNA FoV on map
fovtlist = t0 + datetime.timedelta(minutes=2) * np.array([4]) + datetime.timedelta(minutes=60)
### Need instance of sc
sc = pyana.pep.pep_attitude.NadirLookingSc()
rg = 2634.1
for fovt in fovtlist:
scpos_gany = juice.get_position(fovt, relative_to='GANYMEDE', frame='IAU_GANYMEDE')
scvel_gany = juice.get_velocity(fovt, relative_to='GANYMEDE', frame='IAU_GANYMEDE')
sc.set_posvel(scpos_gany, scvel_gany)
for ch in range(0, 7):
thetalist = np.array([-3., 0, 3, 3, 3, 3, 3, 3, 3, 3, 3, 0, -3, -3, -3, -3, -3, -3, -3, -3, -3], dtype=np.float)
philist = -70 + 20 * ch + np.array([0., 0, 0, 2.5, 5, 7.5, 10, 12.5, 15, 17.5, 20, 20, 20, 17.5, 15, 12.5, 10, 7.5, 5, 2.5, 0], dtype=np.float)
jnacoords = np.array([frame.angles2jna(t, p) for t, p in zip(thetalist, philist)])
pepcoords = frame.jna2nsc(jnacoords.T).T
ganycoords = sc.convert_to_ref(pepcoords.T).T
xsect = np.ma.ones(ganycoords.shape) * np.nan
for idx, vec in enumerate(ganycoords):
v2 = pyana.util.intersection.los_limb_sphere(scpos_gany, vec, [0, 0, 0], rg)
if v2 != None:
xsect[idx, :] = np.array([v2.x, v2.y, v2.z])
xsect = np.ma.masked_invalid(xsect)
xsectn = xsect / rg
lat = np.arcsin(xsectn[:, 2].compressed()) * 180. / np.pi
lon = np.unwrap(np.arctan2(xsectn[:, 1].compressed(), xsectn[:, 0].compressed()))
print(lon)
lon = lon * 180. / np.pi
x, y = m(lon, lat)
m.plot(x, y, 'r')
# For 5 deg FOV -- Copy and paste...
thetalist = np.array([-2.5, 0, 2.5, 2.5, 2.5, 0, -2.5, -2.5, -2.5])
philist = np.array([-2.5, -2.5, -2.5, 0, 2.5, 2.5, 2.5, 0, -2.5])
jnacoords = np.array([frame.angles2jna(t, p) for t, p in zip(thetalist, philist)])
pepcoords = frame.jna2nsc(jnacoords.T).T
ganycoords = sc.convert_to_ref(pepcoords.T).T
xsect = np.ma.ones(ganycoords.shape) * np.nan
for idx, vec in enumerate(ganycoords):
v2 = pyana.util.intersection.los_limb_sphere(scpos_gany, vec, [0, 0, 0], rg)
if v2 != None:
xsect[idx, :] = np.array([v2.x, v2.y, v2.z])
xsect = np.ma.masked_invalid(xsect)
xsectn = xsect / rg
lat = np.arcsin(xsectn[:, 2].compressed()) * 180. / np.pi
lon = np.unwrap(np.arctan2(xsectn[:, 1].compressed(), xsectn[:, 0].compressed()))
print(lon)
lon = lon * 180. / np.pi
x, y = m(lon, lat)
m.plot(x, y, 'y')
#### End of CaP
fig.savefig('flyover_mini.png')
if __name__ == "__main__":
main()
#import unittest
#import doctest
#
#
#def doctests():
# return unittest.TestSuite((
# doctest.DocTestSuite(),
# ))
#if __name__ == '__main__':
# unittest.main(defaultTest='doctests')