appsUnsorted.plot_scpos_ganymede

''' Plot the spacecraft position at a time range given.

This will give a map of the coverage of JUICE with
0.33 degree gridded map.

To simulate the height difference, mainly for 200 km
and 500 km altitude, a simple trick is used.

200 km height will give 0.4 degrees surface area with
FOV=5 degrees.  500 km height is 1 degree at surface.
Thus, the gridding at the sphere is settled to be 0.33 deg,
and for 200 km height I set only 1 point in the spacecraft
footprint, and for 500 km height 5 points at and four neighbor
of the spacecraft footprint.
'''

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

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

import spice

import pyana.util.gridsphere

import pyana.pep.juice_spice
import pyana.pep.moon_image

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

def main():

    setup_kernel()
    timemaskfunction = pyana.pep.juice_spice.is_valid_interval_mantra001

    juice = pyana.pep.juice_spice.JuiceSpice()
    juice_onr = pyana.pep.juice_spice.OrbitNumber001()

    o0, o1 = juice_onr.orb0, juice_onr.orb1   # Full period
#    o0, o1 = 1, 3
    dt = datetime.timedelta(seconds=15)

    fig = plt.figure(figsize=(15, 12))
    ax = fig.add_subplot(111)
    gmap = pyana.pep.moon_image.get_ganymede_image()
    ax.imshow(gmap, origin='bootom', extent=[-180, 180, -90, 90],
                    aspect='equal', cmap=matplotlib.cm.gray)
    ax.imshow(gmap, origin='bootom', extent=[180, 540, -90, 90],
                    aspect='equal', cmap=matplotlib.cm.gray)

    rg = 2634.1

    grid = pyana.util.gridsphere.SimpleGridSphere(nlon=1080, nlat=540)
    dlon = 1/3.
    dlat = 1/3.

    import time
    realstart = time.time()
    
    for oo in range(o0, o1 + 1):
        current = time.time()
        print(oo, '%d/%d' % (oo, o1 - o0 + 1), '%f passed.' % (current - realstart), '%f sec expected.' % ((current - realstart) * (o1 - o0 + 1)/ oo))
        t0 = juice_onr.get_starttime(oo)
        t1 = juice_onr.get_stoptime(oo)

        t = t0

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

            pos = juice.get_position(t, relative_to='ganymede',
                    frame='iau_ganymede')

            h = np.linalg.norm(pos)
            lat = np.arcsin(pos[2] / h)
            lon = np.arctan2(pos[1], pos[0])

#            print t, h, lon, lat, pos

            h -= rg
            lat, lon = np.rad2deg((lat, lon))   # Spacecraft position

            if h >= 350:   # For 500 and 5000 km altitude data
                grid.append_value_lonlat(lon + dlon, lat, h)
                grid.append_value_lonlat(lon, lat, h)
                grid.append_value_lonlat(lon - dlon, lat, h)
                try:
                    grid.append_value_lonlat(lon, lat + dlat, h)
                except:
                    pass
                try:
                    grid.append_value_lonlat(lon, lat - dlat, h)
                except:
                    pass
                
            else:   # For 200 km altitude data
                grid.append_value_lonlat(lon, lat, h)

            t = t + dt

#        draw = ax.scatter(lonlist, latlist, s=1, alpha=1, c=-np.log10(hlist), vmin=-np.log10(5000), vmax=-np.log10(200), edgecolors='none')
    bgrid = grid.get_bgrid()
    min_height = grid.get_min()

    clr = -np.log10(min_height)

    draw = ax.pcolor(bgrid[0], bgrid[1], clr, vmin=-np.log10(5000), vmax=-np.log10(200))

#    clb = plt.colorbar(draw)
#    clb.set_label('log(height [km])')

    ax.set_xlim(0, 360)
    ax.set_ylim(-90, 90)
    ax.set_aspect('equal')
        

if __name__ == '__main__':
    main()