appsUnsorted.juice_orbit_quicklook

''' There are several juice orbit files. Compare them.
'''

import os
import sys
import logging
logging.basicConfig()
import datetime
import dateutil.parser
import math
from optparse import OptionParser

import matplotlib.pyplot as plt
import matplotlib.dates
import numpy as np
import scipy as sp

import spice

import pyana.pep.juice_spice


def makeall(kernels=None, t0=datetime.datetime(2025, 1, 1), t1=datetime.datetime(2030, 1, 1), dt=datetime.timedelta(seconds=3600), timemaskfunction=None):

    if timemaskfunction == None:  # Or use pyana.pep.juice_spice.is_valid_interval_mantra001
        timemaskfunction = lambda x: True

    if kernels == None:
        kernels = [os.path.join('..', 'pyana', 'pep', 'spice_juice', 'mantra.jgo_2020_001_ipc_cal_pso_res_e40_562_200.bsp')]

    kernsdir = os.path.join('..', 'pyana', 'pep', 'spice_juice')
    kerns = ['naif0010.tls', 'pck00009.tpc', 'de405.bsp', 'jup230.bsp', ]

    for k in kerns:
        kern = os.path.join(kernsdir, k)
        print(kern)
        spice.furnsh(kern)

    spice.furnsh(os.path.join('..', 'pyana', 'pep', 'spice', 'jse_111130.tf'))

    for k in kernels:
        print(k)
        spice.furnsh(k)

    frame = "JSE"
    center = "JUPITER"
    center_id = spice.bodn2c(center)

    moons = ["IO", 'EUROPA', 'GANYMEDE', "CALLISTO"]
    moons_id = [spice.bodn2c(m) for m in moons]
    print(moons, moons_id)

    spacecraft = -999

    xyz = []
#    vxyz = []
    etlist = []

    xyz_moon = [[], [], [], []]
        
    t = t0

    while t <= t1:
        if not timemaskfunction(t):
            t = t + dt
            continue

        et = spice.str2et(t.strftime('%FT%T'))

        try:
            posvel = spice.spkez(spacecraft, et, frame, 'LT+S', center_id)[0]
        except spice.SpiceException:
            posvel = [np.nan, np.nan, np.nan]

        xyz.append(posvel[:3])
#        vxyz.append(posvel[3:])

        for i in range(4):
            try:
                pos = spice.spkez(spacecraft, et, frame, 'LT+S', moons_id[i])[0][:3]
            except spice.SpiceException:
                pos = [np.nan, np.nan, np.nan]
            xyz_moon[i].append(pos)

        etlist.append(matplotlib.dates.date2num(t))

        logging.debug('et=%e / t=%f / pos_x=%f' % (et, matplotlib.dates.date2num(t), posvel[0]))

        t = t + dt

    xyz = np.array(xyz) / 71492.
    xyz = np.ma.masked_invalid(xyz)

    etlist = np.array(etlist)
    etlist = np.ma.masked_array(etlist, xyz[:, 0].mask)

    if xyz.mask.all():
        print('no data')
        return

    fig = plt.figure()
    ax = fig.add_subplot(221)
    ax.plot(xyz[:, 0], xyz[:, 1], '.-')
    for rkm in (421600, 670900, 1070000, 1833000):
        r = rkm / 71492.
        ax.plot(r * np.cos(np.linspace(0, 2*np.pi, 100)),
                r * np.sin(np.linspace(0, 2*np.pi, 100)), 'k--')
    ax.set_xlim(-100, 100)
    ax.set_ylim(-100, 100)
    ax.set_xlabel('x [Rj]')
    ax.set_ylabel('y [Rj]')
    ax.set_aspect('equal')

    dist = np.sqrt(xyz[:, 0] ** 2 + xyz[:, 1] ** 2 + xyz[:, 2] ** 2)
    xyzio = np.ma.masked_invalid(xyz_moon[0])
    distio = np.sqrt(xyzio[:, 0] ** 2 + xyzio[:, 1] ** 2 + xyzio[:, 2] ** 2) / 3632. * 2
    xyzeu = np.ma.masked_invalid(xyz_moon[1])
    disteu = np.sqrt(xyzeu[:, 0] ** 2 + xyzeu[:, 1] ** 2 + xyzeu[:, 2] ** 2) / 3138. * 2
    xyzga = np.ma.masked_invalid(xyz_moon[2])
    distga = np.sqrt(xyzga[:, 0] ** 2 + xyzga[:, 1] ** 2 + xyzga[:, 2] ** 2) / 5262. * 2
    xyzca = np.ma.masked_invalid(xyz_moon[3])
    distca = np.sqrt(xyzca[:, 0] ** 2 + xyzca[:, 1] ** 2 + xyzca[:, 2] ** 2) / 4820. * 2

#    vxyz = np.array(vxyz)  # in km/s
#    vxyz = np.ma.masked_invalid(vxyz)
#    vel = np.sqrt(vxyz[:, 0] ** 2 + vxyz[:, 1] ** 2 + vxyz[:, 2] ** 2)


    ax2 = fig.add_subplot(212)
    ax2.plot_date(etlist, dist, fmt='-', label='J')
    ax2.plot_date(etlist, distio, fmt='-', label='I')
    ax2.plot_date(etlist, disteu, fmt='-', label='E')
    ax2.plot_date(etlist, distga, fmt='-', label='G')
    ax2.plot_date(etlist, distca, fmt='-', label='C')
    ax2.set_ylabel('Distance [R_body]')
    ax2.set_yscale('log')
    ax2.legend(loc='best')

    fmt = matplotlib.dates.DateFormatter('%FT%T')
    ax2.xaxis.set_major_formatter(fmt)

    ax2.set_xlim(t0, t1)

    fig.autofmt_xdate()

    
    
def main():

    usage = "usage: %prog [options]"
    parser = OptionParser(usage)

    parser.add_option("-k", "--kernel", dest="kernels", default=[],
        help="Load non-default kernel(s) of FILENAME", action='append')
    parser.add_option("-s", dest="start_time", default='2025-01-01T00:00:00',
        help="Set start time")
    parser.add_option("-e", dest="stop_time", default='2030-01-01T00:00:00',
        help="Set stop time")

    parser.add_option("-d", dest="delta_time", default=3600,
        help="Set time resolution in seconds.", type='int')

    parser.add_option("-v", "--verbose",
            action="store_true", dest="verbose")
    parser.add_option("-q", "--quiet",
            action="store_false", dest="verbose")

    (options, args) = parser.parse_args()

    if len(args) != 0:
        parser.error("incorrect number of arguments")

    if options.verbose:
        logging.getLogger().setLevel(logging.DEBUG)

    start_time = dateutil.parser.parse(options.start_time)
    stop_time = dateutil.parser.parse(options.stop_time)
    delta_time = datetime.timedelta(seconds=options.delta_time)
    logging.debug('Time:  %s --(%s)--> %s' % (str(start_time),
                    repr(delta_time), str(stop_time)))

    print('K=', options.kernels)
    if len(options.kernels) == 0:
	options.kernels = None
	timemaskfunction = pyana.pep.juice_spice.is_valid_interval_mantra001
    else:
        timemaskfunction = None
    makeall(kernels=options.kernels, t0=start_time, t1=stop_time, dt=delta_time,
        timemaskfunction=timemaskfunction)

if __name__ == '__main__':
    main()