appsUnsorted.juice_ganymede_pseudo_orbitnr

''' This script is for calculation of the pseudo-orbit numbering.

The pseudo-orbit numbering would be helpful in case you are interested in orbit numbering sorted pre-analysis.
The number is defined from north pole to north pole.
If there is a "data gap", the corresponding orbit number is ended at the time of the
start of the data gap, and the next orbit number starts immediately ending at the end of the data gap.  The next will only start the end of the data gap.
This means that 1 data gap corresponds to 1 orbit number, which includes no data.
The origin, i.e. orbit number 1, is defined in the script.

The script will print the results to standard output as a python dictionary form.


'''



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

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

import spice
import pyana.pep.juice_spice

def main():

    t0 = datetime.datetime(2028, 9, 21, 23)  # When Rg < Rj in .
    t1 = datetime.datetime(2030, 1, 1)
#    t1 = datetime.datetime(2028, 10, 19, 1)

    jsp = pyana.pep.juice_spice.JuiceSpice()
    jsp.load_default_kernels()

    kernels = [os.path.join('..', 'pyana', 'pep',
                'spice_juice', 'mantra.jgo_2020_001_ipc_cal_pso_res_e40_562_200.bsp')]
    for k in kernels:
        spice.furnsh(k)

    timemaskfunction = pyana.pep.juice_spice.is_valid_interval_mantra001

    # Just calculate the latitude in Ganymede every 1 minutes,
    # then, compare with the latitude 1 minutes before,
    # and output.

    dt = datetime.timedelta(seconds=60)
    t = t0

    pos_t_minus_dt = jsp.get_position(t0, relative_to="Ganymede", frame="IAU_GANYMEDE")
    pos_t = jsp.get_position(t0 + dt, relative_to='Ganymede', frame="IAU_GANYMEDE")
    lat0 = np.arcsin(pos_t_minus_dt[2] / np.linalg.norm(pos_t_minus_dt)) / np.pi * 180.
    lat1 = np.arcsin(pos_t[2] / np.linalg.norm(pos_t)) / np.pi * 180.
    northgoing = (lat0 < lat1)

    current_onr = 1

    nodata = False

    print("import dateutil.parser")
    print("pseudo_onr = {")
    while t < t1:
        t = t + dt

        if not timemaskfunction(t):
            if not nodata:
                print("%d %s [nan nan nan] %f nan" % (current_onr, t, lat0), file=sys.stderr)
                print('%d: dateutil.parser.parse("%s"), # datagap' % (current_onr, t))
                current_onr += 1
            nodata = True
            continue
        nodata = False

        pos_t = jsp.get_position(t, relative_to='Ganymede', frame='IAU_GANYMEDE')

        height = np.linalg.norm(pos_t)
        lat = np.arcsin(pos_t[2] / height) / np.pi * 180.
        lon = np.arctan2(pos_t[1], pos_t[0]) / np.pi * 180.

        if northgoing and not (lat0 < lat):
            print("%d %s %s %f %f" % (current_onr, t, pos_t, lat0, lat), file=sys.stderr)
            print('%d: dateutil.parser.parse("%s"), ' % (current_onr, t))
            current_onr += 1

        northgoing = (lat0 < lat)
        lat0 = lat  # Replace the present position to the previous one.

        sys.stdout.flush()

    print("}")

    



if __name__ == '__main__':
    main()