appsUnsorted.plot_galileo_moons_orbit

''' Plot Gallilean Moon position and Galileo satellite orbit.
'''
import sys
import datetime
import dateutil.parser
import calendar

from optparse import OptionParser

import matplotlib.pyplot as plt
import numpy as np

from matplotlib.patches import Circle, Wedge

import pyana.pep.pep_spice
import pyana.pep.galileo_orbit

def get_traj(moon, tc, dur_orbit, del_orbit, frame, center):
    ''' Return the trajectory of the moon

    :returns: np.array with (N, 3) shape where N is the data number
        calculated from ``tc``, ``dur_orbit`` and ``del_orbit``.
    '''
    # Plotting oribt, start and end time
    t0 = tc - dur_orbit
    t1 = tc + dur_orbit

    torblist = []
    torblist = []

    t = tc
    while t >= t0:
        torblist.append(t)
        t = t - del_orbit

    t = tc
    while t <= t1:
        torblist.append(t)
        t = t + del_orbit

    torblist = list(set(torblist))
    torblist.sort()

    pepsp = pyana.pep.pep_spice.get_default_pepspice()

    moons_orbit = []
    for t in torblist:
        pos, vel = pepsp.get_posvel(moon, t, frame, center)
        pos = pos
        moons_orbit.append(pos)

    return np.array(moons_orbit)


def plotmain(tc, output=None, duration=None):

    moons = ['Io', 'Europa', 'Ganymede', 'Callisto']
    rmoon = {'Io': 1816., 'Europa': 1569., 'Ganymede': 2631., 'Callisto': 2410.}
    center = "Jupiter"
    frame = "JSE"

    # For plotting orbit.  Duration and resolution
#    dur_orbit = datetime.timedelta(days=1)
    del_orbit = datetime.timedelta(hours=1)
#    label_orbit = datetime.timedelta(hours=6)

    norm = 71492.

    pepsp = pyana.pep.pep_spice.get_default_pepspice()

    moons_positions = [pepsp.get_posvel(moon, tc, frame, center)[0] / norm for moon in moons]

    moons_orbits = [get_traj('Io', tc, datetime.timedelta(days=1.769) / 2, del_orbit, frame, center) / norm,
                    get_traj('Europa', tc, datetime.timedelta(days=3.551) / 2, del_orbit, frame, center) / norm,
                    get_traj('Ganymede', tc, datetime.timedelta(days=7.155) / 2, del_orbit, frame, center) / norm,
                    get_traj('Callisto', tc, datetime.timedelta(days=16.69) / 2, del_orbit, frame, center) / norm,
                    ]

    # Galileo orbit

    gal = pyana.pep.galileo_orbit.GalileoOrbit()

    if tc < gal.t0() or tc > gal.t1():
        raise ValueError('Given time (%s) is not included in Galileo orbit data (%s - %s)' % (tc, t0, t1))
    

    del_gal = datetime.timedelta(minutes=120)  # Every 1 minutes
    if duration == None:
        t0 = gal.t0()
        t1 = gal.t1()
    else:
        t0 = tc - datetime.timedelta(seconds=duration)
        t1 = tc + datetime.timedelta(seconds=duration)

    tlist = [tc]
    tt = tc - del_gal
    while tt >= t0:
        print(tt, len(tlist))
        tlist.append(tt)
        tt = tt - del_gal
    tt = tc + del_gal
    while tt <= t1:
        print(tt, len(tlist))
        tlist.append(tt)
        tt = tt + del_gal

    tlist.sort()

    print('Time list length = %d' % len(tlist))

    poslist = gal.get_positions(tlist)

    gal_pos = gal.get_positions(tc)

    # Plotting

    fig = plt.figure()
    ax = fig.add_subplot(111)

#    ax.plot([0], [0], 'kx', label=center)
    jupiter_body = Circle((0, 0), 71492. / norm, edgecolor='k', fill=False, alpha=0.5)
    ax.add_patch(jupiter_body)
    jupiter_dark = Wedge((0, 0), 71492. / norm, 90., 270., color='k', alpha=0.5)
    ax.add_patch(jupiter_dark)
    ax.annotate('J', (0, 0))

    for moonspos, moon in zip(moons_positions, moons):
        ax.plot([moonspos[0]], [moonspos[1]], 'k.')
        moon_body = Circle((moonspos[0], moonspos[1]), rmoon[moon] / norm, edgecolor='k', fill=False, alpha=0.5)
        ax.add_patch(moon_body)
        ax.annotate(moon[0], (moonspos[0], moonspos[1]))

    for moons_orbit in moons_orbits:
        ax.plot(moons_orbit[:, 0], moons_orbit[:, 1], 'k--')

    ax.plot(poslist[:, 0], poslist[:, 1], 'k:')
    ax.plot([gal_pos[0]], [gal_pos[1]], 'ko')
    ax.annotate('SC', (gal_pos[0], gal_pos[1]))

    ax.set_xlim(-30, 30)
    ax.set_ylim(-30, 30)

    xlim = ax.get_xlim()
    ax.set_xlim(xlim[::-1])  # To swap. Left is the sun.
    ax.grid()
    ax.set_aspect('equal')

    ax.set_xlabel('X_JSE [Rj]')
    ax.set_ylabel('Y_JSE [Rj]')


    ax.set_title('%s' % tc.strftime('%F %T'))

    if output != None:
        fig.savefig(output)

    return fig

def main():
    # Parsing the argument

    usage = 'USAGE:  %prog  [options]  yyyy-mm-ddThh:mm:ss'
    parser = OptionParser(usage)

    parser.add_option('-o', dest='outputfigname', default=None, help='Output figure name')
    parser.add_option('-d', dest='duration', default=None,
            help='Duration of plotting Galileo trajectory in seconds w.r.t the given time.')

    (options, args) = parser.parse_args()
    if len(args) != 1:
        parser.error('Must specify the time.')

    fig = plotmain(dateutil.parser.parse(args[0]), output=options.outputfigname, duration=options.duration)
    
    return fig

if __name__ == "__main__":
    val = main()