apps120712_spice1205.dump_juice1205_2

''' Dump the 1205 version of spice kernel

Command line options:
'''

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

import math

import numpy as np

import pyana.pep.juice_spice1205 as jspice
from pyana.spice.pyspice import spice

def caldt(x, v, r):

    xx = math.sqrt((x * x).sum()) - r
    vv = math.sqrt((v * v).sum())

    t = xx / vv / 10.0

    return t   # in seconds.
    

def tfloor(t):
    ''' Return as follows.

    t > 86400  ->  1 day
    t > 3600   ->  1 hour
    t > 60     ->  1 min
    otherwize  ->  10 sec
    '''

    if t > 86400:
        return 86400
    elif t > 3600:
        return 3600
    elif t > 60:
        return 60
    else:
        return 10

def tupdate(t, tmin):
    ''' Update the time.
    '''
    tmin2 = tfloor(tmin)

    tnext = t + datetime.timedelta(seconds=tmin2)

    # Replace the time in case not on the hour, munutes or day.
    if tmin2 == 60:
        tnext = tnext.replace(second=0)
    elif tmin2 == 3600:
        tnext = tnext.replace(second=0, minute=0)
    elif tmin2 == 86400:
        tnext = tnext.replace(second=0, minute=0, hour=0)
        
    return tnext


def main(t0=datetime.datetime(2022, 6, 1, 7),
                t1=datetime.datetime(2033, 7, 4, 18)):

    print('#yymmdd hhmmss x_jse y_jse z_jse r_jup x_io y_io z_io r_io x_europa y_europa z_europa r_europa x_ganymede y_ganymede z_ganymede r_ganymede x_callisto y_callisto z_callisto r_callisto')

    js = jspice.JuiceSpice.get_default_instance()

    t = t0

    dt = datetime.timedelta(seconds=1)

    logging.getLogger('JuiceSpice').setLevel(logging.WARN)
    while t <= t1:
        # Convert to et.
#        et = spice.utc2et(t.strftime('%FT%T'))

        # Get position in the proper reference frame.
        ju = js.get_position(t, relative_to='Jupiter', frame='JSE')
        io = js.get_position(t, relative_to='Io', frame='IAU_IO')
        eu = js.get_position(t, relative_to='Europa', frame='IAU_EUROPA')
        ga = js.get_position(t, relative_to='Ganymede', frame='IAU_GANYMEDE')
        ca = js.get_position(t, relative_to='Callisto', frame='IAU_CALLISTO')

        if np.isnan(np.array([ju[0], io[0], eu[0], ga[0], ca[0]])).any():
            t = t + dt
            continue

        print(t, end=' ')
        print(ju[0], ju[1], ju[2], np.sqrt((ju * ju).sum()), end=' ')
        print(io[0], io[1], io[2], np.sqrt((io * io).sum()), end=' ')
        print(eu[0], eu[1], eu[2], np.sqrt((eu * eu).sum()), end=' ')
        print(ga[0], ga[1], ga[2], np.sqrt((ga * ga).sum()), end=' ')
        print(ca[0], ca[1], ca[2], np.sqrt((ca * ca).sum()), end=' ')
        print()

        ### Calculate next time.

        # Calculate velocity in JSE.
        vju = js.get_velocity(t, relative_to='Jupiter', frame='JSE')
        vio = js.get_velocity(t, relative_to='Io', frame='JSE')
        veu = js.get_velocity(t, relative_to='Europa', frame='JSE')
        vga = js.get_velocity(t, relative_to='Ganymede', frame='JSE')
        vca = js.get_velocity(t, relative_to='Callisto', frame='JSE')

        # Calculate minimum update time.
        tju = caldt(ju, vju, 65000.)
        tio = caldt(io, vio, 1821.3)
        teu = caldt(eu, veu, 1569.)
        tga = caldt(ga, vga, 2634.1)
        tca = caldt(ca, vca, 2410.3)

        t_min = np.min(np.array([tju, tio, teu, tga, tca]))

        # next time.
        t = tupdate(t, t_min)

if __name__ == "__main__":
    optlist, args = getopt.getopt(sys.argv[1:], 'c:f:l:s:e:d:')

    # Default values
    start = datetime.datetime(2022, 6, 1, 7)
    end = datetime.datetime(2033, 7, 4, 18)

    main(t0=start, t1=end)