appsUnsorted.plot_moons_orbit

''' Plot Gallilean Moon position
'''
import sys
import datetime
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

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, unit='km', frame='JSE', output=None, with_longitude=False):

    # Argument to get from commandline, but not yet supported.
    moons = ['Io', 'Europa', 'Ganymede', 'Callisto']
    rmoon = {'Io': 1816., 'Europa': 1569., 'Ganymede': 2631., 'Callisto': 2410.}
    center = "Jupiter"

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

    # Unit conversion
    if unit == 'km':
        norm = 1.
    elif unit == 'Rj':
        norm = 71492.
    else:
        raise ValueError('No unit definition for %s' % unit)

    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,
                    ]

    # 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--')


    # x axis in IAU_IO to be plotted on JSE.
    if with_longitude:
        # IAU_XXX frame
        xaxis = np.array([100000, 0, 0])
        for moon_pos, moon in zip(moons_positions, moons):
            xaxis_jse = pepsp.convert_vector(xaxis, 'IAU_%s' % moon, 'JSE', tc) / norm
            ax.arrow(moon_pos[0], moon_pos[1], xaxis_jse[0], xaxis_jse[1], alpha=0.5)


    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_%s [%s]' % (frame, unit))
    ax.set_ylabel('Y_%s [%s]' % (frame, unit))

    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('-u', dest='unit', default='km', help='Unit. "km" or "Rj"')
    parser.add_option('-f', dest='frame', default='JSE', help='Frame.  Should be defined in SPICE')

    parser.add_option('-l', dest='with_longitude', default=False, action='store_true', help='Indicate longitude')

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

    fig = plotmain(datetime.datetime(2025, 1, 1, 0, 0, 0), unit=options.unit, frame=options.frame, output=options.outputfigname, with_longitude=options.with_longitude)

    
    return fig

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