''' 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_nadir_1.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)
### Note that the longitude here is
m = Basemap(projection='cyl', lon_0=90, resolution='c', ax=ax)
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)
# LINEAR scale
# pc = m.contourf(x, y, flxlist, [1e8, 1.25e8, 1.5e8, 1.75e8, 2e8, 3e8])
# fig.colorbar(pc)
### 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.arange(10) + 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')
fig.savefig('flyover_%s_1.png' % type)
fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)
jnatime = t0 + datetime.timedelta(minutes=2) * np.arange(90)
jnacnt = np.zeros((len(jnatime), 7)) + np.nan
import pyana.jdc.fluxsimulator
jdcsim = pyana.jdc.fluxsimulator.SimpleSimulator()
jdcsim.set_species(16.0, 1.0)
pp = pyana.pep.mhddata.PlasmaParameter1205()
jdcflx = np.zeros((len(jnatime), 128)) + np.nan
for ilat, t in enumerate(jnatime):
print(ilat, t)
scpos_gany = juice.get_position(t, relative_to='GANYMEDE', frame='IAU_GANYMEDE')
scvel_gany = juice.get_velocity(t, relative_to='GANYMEDE', frame='IAU_GANYMEDE')
# Convert to plasma frame
scpos_gany = geo2pla(scpos_gany) / rg
scvel_gany = geo2pla(scvel_gany)
sc.set_posvel(scpos_gany, scvel_gany)
for ch in (1, 2, 3, 4, 5): # Only center 5 channels.
# Center direction
az = jnafov.azim_pix_center(ch)
el = jnafov.elev_pix_center(ch)
vecjna = jnaframe.angles2jna(el, az)
vecnsc = jnaframe.jna2nsc(vecjna)
vecgan = sc.convert_to_ref(vecnsc) # Looking vector in Ganymede coordinates.
v2 = pyana.util.intersection.los_sphere(scpos_gany, vecgan, [0, 0, 0], 1)
if v2 == None:
continue
longit = np.arctan2(v2[1], v2[0]) * 180. / np.pi
if longit < 0: longit += 360
latit = np.arcsin(v2[2]) * 180. / np.pi
print('LONLAT', longit, latit)
try:
idxlon, idxlat = fmap.nearest_neighbor(longit, latit)
except IndexError as e:
from pudb import set_trace; set_trace()
print('Strange, index error. Anyway, skip it.', longit, latit)
continue
jnacnt[ilat, ch] = flxlist[idxlon, idxlat] * 0.2 / (2 * np.pi) * 2e-6 # in /s. Indeed , count rate
# 0.2 comes from 20% from CENA result.
# Sensor g-factor is taken from proposal v0.3, 2e-6 cm2 sr/pixel.
# print jnacnt[ilat, ch]
### JDC
n, vx, vy, vz, temp, pth = pp.interpolate3d(scpos_gany[0], scpos_gany[1], scpos_gany[2])
vth = pyana.pep.mhddata.t2vth(temp, mass=16.) # in m/s
n *= 1e6 # in /m3
vx *= 1e3; vy *= 1e3; vz *= 1e3 # in m/s
jdcsim.set_spacecraft_pvat(scpos_gany, scvel_gany)
jdcsim.set_plasmaparameter(n, vx, vy, vz, vth)
for ie in range(128):
jdcflx[ilat, ie] = jdcsim.get_countrate(ie, 31, 0)
print(jdcflx[ilat, :])
# xax = [(j - jnatime[0]).seconds for j in jnatime]
import matplotlib.dates as mpldates
xax = [mpldates.date2num(j) for j in jnatime]
yax = np.arange(8)
x, y = np.meshgrid(xax, yax)
# img = ax2.pcolor(x, y, np.ma.log10(jnacnt).T)
img = ax2.pcolor(x, y, np.ma.masked_invalid(jnacnt).T, vmin=0, vmax=10)
ax2.set_ylabel('JNA: Direction')
cb = plt.colorbar(img, ax=ax2)
cb.set_label('JNA count rate (log)')
cb.set_label('JNA count rate')
ax2.xaxis.set_major_formatter(mpldates.DateFormatter('%H:%M'))
# locator = mpldates.AutoDateLocator()
# ax2.xaxis.set_major_locator(locator)
ax2.xaxis.set_ticks([mpldates.date2num(j) for j in tick_time])
ax2.set_xlim(xax[0], xax[-1])
import pyana.jdc.energy0 as energy
yax = energy.getEnergy()
x, y = np.meshgrid(xax, yax)
img = ax1.pcolor(x, y, np.ma.log10(jdcflx).T, vmin=0.1, vmax=4)
ax1.set_ylabel('JDC: Energy [eV]')
ax1.set_yscale('log')
ax1.set_ylim(yax.min(), yax.max())
cb = plt.colorbar(img, ax=ax1)
cb.set_label('JDC count rate (log)')
ax1.xaxis.set_ticks([mpldates.date2num(j) for j in tick_time])
ax1.xaxis.set_ticklabels('')
# ax1.xaxis.set_major_locator(locator)
ax1.set_xlim(xax[0], xax[-1])
fig.autofmt_xdate()
fig.savefig('flyover_%s_2.png' % type)
if __name__ == "__main__":
main()
main(type='footprint')
#import unittest
#import doctest
#
#
#def doctests():
# return unittest.TestSuite((
# doctest.DocTestSuite(),
# ))
#if __name__ == '__main__':
# unittest.main(defaultTest='doctests')