''' Plot flyby informations
'''
import datetime
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mpldates
import matplotlib.ticker
import pyana.util.timeseries
import pyana.juice.jspice as js
js.init()
import pyana.juice.mission_summary as jms
radius = {'Ganymede': 2634, 'Callisto': 2410, 'Europa': 1561, 'Jupiter': 71490}
def plotall(eventname, tcenter, catarget):
''' Plot.
'''
t0 = tcenter - datetime.timedelta(hours=12)
t1 = tcenter + datetime.timedelta(hours=12)
t0_1hr = tcenter - datetime.timedelta(hours=1)
t1_1hr = tcenter + datetime.timedelta(hours=1)
d1min = datetime.timedelta(minutes=1)
d10min = datetime.timedelta(minutes=10)
d60min = datetime.timedelta(minutes=60)
tl1min = pyana.util.timeseries.dtrange(t0, t1, d1min)
tl10min = pyana.util.timeseries.dtrange(t0, t1, d10min)
tl60min = pyana.util.timeseries.dtrange(t0, t1, d60min)
fig = plt.figure(figsize=(16, 16))
### JSE frame
ax = fig.add_subplot(231)
cent = len(tl1min) / 2
io_pos = np.array([js.get_position(t, frame='JSE', origin='Jupiter', target='Io') for t in tl1min]) / radius['Jupiter']
ax.plot(io_pos[:, 0], io_pos[:, 1], 'k:')
ax.plot(io_pos[cent, 0], io_pos[cent, 1], 'ko')
europa_pos = np.array([js.get_position(t, frame='JSE', origin='Jupiter', target='Europa') for t in tl1min]) / radius['Jupiter']
ax.plot(europa_pos[:, 0], europa_pos[:, 1], 'k:')
ax.plot(europa_pos[cent, 0], europa_pos[cent, 1], 'ko')
ganymede_pos = np.array([js.get_position(t, frame='JSE', origin='Jupiter', target='Ganymede') for t in tl1min]) / radius['Jupiter']
ax.plot(ganymede_pos[:, 0], ganymede_pos[:, 1], 'k:')
ax.plot(ganymede_pos[cent, 0], ganymede_pos[cent, 1], 'ko')
callisto_pos = np.array([js.get_position(t, frame='JSE', origin='Jupiter', target='Callisto') for t in tl1min]) / radius['Jupiter']
ax.plot(callisto_pos[:, 0], callisto_pos[:, 1], 'k:')
ax.plot(callisto_pos[cent, 0], callisto_pos[cent, 1], 'ko')
juice_pos = np.array([js.get_position(t, frame='JSE', origin='Jupiter', target='JUICE') for t in tl1min]) / radius['Jupiter']
ax.plot(juice_pos[:, 0], juice_pos[:, 1], 'b:')
# ax.plot(juice_pos[cent, 0], juice_pos[cent, 1], 'bo')
ax.set_aspect(1)
from matplotlib.patches import Circle
jup = Circle((0, 0), 1, color='k')
ax.add_patch(jup)
ax.set_aspect(1)
ax.set_xlabel('JSE x [Rj]')
ax.set_ylabel('JSE y [Rj]')
xrng = ax.get_xlim()
ax.set_xlim(xrng[1], xrng[0])
yrng = ax.get_ylim()
ax.set_ylim(yrng[1], yrng[0])
### Height profile
from . import height_profile
ax = fig.add_subplot(232)
ax = height_profile.main(t0, t1, d10min, origin=catarget, ax=ax)
ax.axhline(radius[catarget], color='k')
# ax.set_yscale('log')
ax.set_title(eventname + " " + tcenter.strftime('%FT%T'))
axtwin = ax.twinx()
axylim = ax.get_ylim()
axtwin.set_ylim(axylim[0] / radius[catarget], axylim[1] / radius[catarget])
locator=mpldates.AutoDateLocator(maxticks=8)
# locator.intervald[mpldates.MINUTELY] = [20]
ax.xaxis.set_major_locator(locator)
ax.xaxis.set_major_formatter(mpldates.DateFormatter('%H%M'))
ax.set_ylabel('Height [km]')
ax1 = fig.add_subplot(233)
ax1 = height_profile.main(t0_1hr, t1_1hr, d1min, origin=catarget, ax=ax1)
ax1.axhline(radius[catarget], color='k')
ax1.set_title(eventname + " " + tcenter.strftime('%FT%T'))
ax1twin = ax1.twinx()
ax1ylim = ax1.get_ylim()
ax1twin.set_ylim(ax1ylim[0] / radius[catarget], ax1ylim[1] / radius[catarget])
locator=mpldates.AutoDateLocator(maxticks=8)
locator.intervald[mpldates.MINUTELY] = [20]
ax1.xaxis.set_major_locator(locator)
ax1.xaxis.set_major_formatter(mpldates.DateFormatter('%H%M'))
ax1twin.set_ylabel('Height [R%s]' % catarget[0])
cavec = np.array(js.get_position(tcenter, origin=catarget, frame='IAU_%s' % catarget))
cadist = np.sqrt((cavec * cavec).sum())
ax1.set_title('%.0fkm %.2fR%s from surface' % (cadist - radius[catarget], cadist / radius[catarget] - 1, catarget[0]))
### Slices
from . import slices
tcenter2 = datetime.datetime(tcenter.year, tcenter.month, tcenter.day, tcenter.hour, tcenter.minute / 10 * 10, 0)
t0 = tcenter2 - datetime.timedelta(hours=12)
t1 = tcenter2 + datetime.timedelta(hours=12)
t0_1hr = tcenter2 - datetime.timedelta(hours=1)
t1_1hr = tcenter2 + datetime.timedelta(hours=1)
d1min = datetime.timedelta(minutes=1)
d10min = datetime.timedelta(minutes=10)
d30min = datetime.timedelta(minutes=30)
tl1min = pyana.util.timeseries.dtrange(t0, t1, d1min)
tl10min = pyana.util.timeseries.dtrange(t0, t1, d10min)
tl30min = pyana.util.timeseries.dtrange(t0, t1, d30min)
ax2 = fig.add_subplot(2,3,4)
slices.main(t0_1hr, t1_1hr, d1min, tickdt=d10min, labeldt=d30min, origin=catarget, ax=ax2, projection='yx')
ax3 = fig.add_subplot(2,3,5)
slices.main(t0_1hr, t1_1hr, d1min, tickdt=d10min, labeldt=d30min, origin=catarget, ax=ax3, projection='yz')
ax4 = fig.add_subplot(2,3,6)
slices.main(t0_1hr, t1_1hr, d1min, tickdt=d10min, labeldt=d30min, origin=catarget, ax=ax4, projection='xz')
xrng = ax2.get_xlim()
yrng = ax2.get_ylim()
zrng = ax3.get_ylim()
rng0 = min([xrng[0], -xrng[1], yrng[0], -yrng[1], zrng[0], -zrng[1]])
rng1 = -rng0
print(rng0, rng1)
ax2.set_xlim(rng1, rng0)
ax2.set_ylim(rng0, rng1)
ax2.set_title('See from north')
ax3.set_xlim(rng1, rng0)
ax3.set_ylim(rng0, rng1)
ax3.set_title('See toward Jupiter')
ax4.set_xlim(rng1, rng0)
ax4.set_ylim(rng0, rng1)
ax4.set_title('See from plasma upstream')
# Corotational flow
for cent in np.linspace(rng0, rng1, 6):
ax2.arrow(rng1, cent, (rng0 - rng1) * .1, 0, color='b', width=0.3)
ax3.arrow(rng1, cent, (rng0 - rng1) * .1, 0, color='b', width=0.3)
ax4.add_patch(Circle((0, .9 * rng0), 0.3, color='b'))
ax4.add_patch(Circle((0, .9 * rng1), 0.3, color='b'))
# Bfield
for cent in np.linspace(rng1 * 0.95 + rng0 * 0.05, rng1 * 0.9 + rng0 * 0.1, 6):
ax3.axvline(cent, color='r', alpha=0.5)
for cent in np.linspace(-0.5, 0.5, 6):
ax4.axvline(cent, color='r', alpha=0.5)
ax2.add_patch(Circle((.9 * rng1, .9 * rng0), 0.3, color='r'))
ax2.add_patch(Circle((.9 * rng1, .9 * rng1), 0.3, color='r'))
##### Disabled, because negative and positive (visible and invisible) is not straightforward.
# ### Sun direction
# sunvec = np.array(js.get_position(tcenter2, origin=catarget, target='Sun', frame='IAU_%s' % catarget))
# sunvec = sunvec / np.sqrt((sunvec * sunvec).sum()) * radius[catarget]
# ax2.plot([0, sunvec[0]], [0, sunvec[1]], 'r')
# ax3.plot([0, sunvec[1]], [0, sunvec[2]], 'r')
# ax4.plot([0, sunvec[0]], [0, sunvec[2]], 'r')
# ax2.plot([sunvec[0]], [sunvec[1]], 'rD')
# ax3.plot([sunvec[1]], [sunvec[2]], 'rD')
# ax4.plot([sunvec[0]], [sunvec[2]], 'rD')
# print sunvec
#
# ### Jupiter direction
# jupvec = np.array(js.get_position(tcenter2, origin=catarget, target='Jupiter', frame='IAU_%s' % catarget))
# jupvec = jupvec / np.sqrt((jupvec * jupvec).sum()) * radius[catarget]
# ax2.plot([0, jupvec[0]], [0, jupvec[1]], 'g')
# ax3.plot([0, jupvec[1]], [0, jupvec[2]], 'g')
# ax4.plot([0, jupvec[0]], [0, jupvec[2]], 'g')
# ax2.plot([jupvec[0]], [jupvec[1]], 'gD')
# ax3.plot([jupvec[1]], [jupvec[2]], 'gD')
# ax4.plot([jupvec[0]], [jupvec[2]], 'gD')
# print jupvec
fig.savefig('html/png/flyby-%s.png' % eventname)
def main():
'''Main script'''
evlist = jms.get_label()[:-4] # Only flyby, not circular orbit.
# evlist = jms.get_label()[:-35] # Only several for debuggin
print(evlist)
t0list = [jms.get_datetime(e) for e in evlist]
print(t0list)
flybytarget = {'G': 'Ganymede', 'E': 'Europa', 'C': 'Callisto'}
calist = [flybytarget[e[0]] for e in evlist]
for e, t, c in zip(evlist, t0list, calist):
plotall(e, t, c)
if __name__ == "__main__":
main()