appsUnsorted.juice_orbit_file

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

import os
import sys
import logging
logging.basicConfig()
import datetime
import math

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

import spice

import pyana.pep.juice_spice

def main():
    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'))

    spice.furnsh(sys.argv[1])

    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
    years = list(range(2020, 2031))

    dt = datetime.timedelta(hours=3)

    for year in years:

        xyz = []
        etlist = []

        xyz_moon = [[], [], [], []]
        
        t0 = datetime.datetime(year, 1, 1)
        t1 = datetime.datetime(year+1, 1, 1)

        t = t0

        while t < t1:
            if sys.argv[1].endswith('mantra.jgo_2020_001_ipc_cal_pso_res_e40_562_200.bsp') and (not pyana.pep.juice_spice.is_valid_interval_mantra001(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])

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

            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 for %d' % year)
            continue

        output = "juice_orbit_%s_%04d.png" % (os.path.basename(sys.argv[1]), year)
        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, 10)),
                    r * np.sin(np.linspace(0, 2*np.pi, 10)), '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')
        ax.set_title(output[12:-4])

        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

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

        ax2.set_xlim(t0, t1)


        fig.savefig(output)
    
    

if __name__ == '__main__':
    main()