apps120718_jna_fov.traj_on_surface

''' Plot trajectory on the surface

.. image:: ../../../src/scripts/apps120718_jna_fov/traj_on_surface_500km.png
    :width: 80%

.. image:: ../../../src/scripts/apps120718_jna_fov/traj_on_surface_200km.png
    :width: 80%
'''

import numpy as np
#from pyana.spice.pyspice import spice
import datetime
import matplotlib.pyplot as plt
import matplotlib.dates as mpldates

from mpl_toolkits.basemap import Basemap

import pyana.util.intersection

import pyana.pep.juice_spice1205 as jspice
import pyana.pep.pep_attitude

import pyana.jna.fov0 as fov
import pyana.jna.frame0 as frame

from pyana.pep import mhddata


def draw_enaflux(m):
    f = mhddata.Flux1205()

    lon = f.reshape(f.lonlist())
    lat = f.reshape(f.latlist())
    flx = f.reshape(f.flist())
    flx *= 1.44  # 1.2 Rm map to 1.0 Rm, very rough estimation.

    enaflx = flx * 0.2   # 20% fraction backscattered
    enaflx = enaflx / (2 * np.pi)   # 2 pi steradian, ie isotropic
    enaflx = enaflx / 100   # 100eV width assumed

    enaflx = np.log10(enaflx)  # LOG scale

    x, y = m(lon, lat)
    enaflx = np.ma.masked_less(enaflx, -3)
    p = m.pcolor(x, y, enaflx, vmin=0)

    clb = plt.colorbar(p)
    clb.set_label('LOG of Backscattered ENA flux [/cm2 sr s]')

def plotonmap(t0, t1, dt, fovtlist=None, ticktlist=None, onmap=False):
    '''
    '''
    if fovtlist == None: fovtlist = []

    juice = jspice.JuiceSpice.get_default_instance()

    fig = plt.figure(figsize=(20, 15))
    ax = fig.add_subplot(111)

    m = Basemap(projection='moll', lon_0=0, resolution='c', ax=ax)
#    m = Basemap(projection='merc', lon_0=0, resolution='c', ax=ax)

    draw_enaflux(m)

    tlist = [mpldates.num2date(t).replace(tzinfo=None) for t in mpldates.drange(t0, t1, dt)]
    print('TLIST', tlist)

    poslist = juice.get_positions(tlist, relative_to='GANYMEDE', frame='IAU_GANYMEDE')   # (N, 3) array
    print(poslist)

    rg = 2634.1

    r = np.sqrt((poslist ** 2).sum(axis=1))
    h = r - rg

    posnorm = (poslist.T / r)   # Now (3, N) array
    lat = np.arcsin(posnorm[2, :])
    lon = np.unwrap(np.arctan2(posnorm[1, :], posnorm[0, :]))
    
    x, y = m(np.rad2deg(lon) % 360, np.rad2deg(lat))

    m.plot(x, y, '.')

    # Same for ticktlist
    if ticktlist != None:
        poslist = juice.get_positions(ticktlist, relative_to='GANYMEDE', frame='IAU_GANYMEDE')
        r = np.sqrt((poslist ** 2).sum(axis=1))
        posnorm = (poslist.T / r)
        lat = np.arcsin(posnorm[2, :])
        lon = np.unwrap(np.arctan2(posnorm[1, :], posnorm[0, :]))
        x, y = m(np.rad2deg(lon) % 360, np.rad2deg(lat))

        m.plot(x, y, 'o')
        strlist = [t.strftime('%T') for t in ticktlist]
#        for x_, y_, str_ in zip(x, y, strlist):
#            plt.text(x_, y_, str_)
    

    m.drawparallels(np.arange(-90, 120, 30))
    m.drawmeridians(np.arange(-180, 210, 30))

    sc = pyana.pep.pep_attitude.NadirLookingSc()

    for fovt in fovtlist:
        ### Instance the nadir-looking spacecraft.
        ### See also sample map_on_surface.py.

        scpos_gany = juice.get_position(fovt, relative_to='GANYMEDE', frame='IAU_GANYMEDE')
        scvel_gany = juice.get_velocity(fovt, relative_to='GANYMEDE', frame='IAU_GANYMEDE')
        sc.set_posvel(scpos_gany, scvel_gany)
        print(scpos_gany, scvel_gany)

        for ch in range(0, 7):
#            thetalist, philist = fov.pixel_shape(ch, azimresolution=10, elevresolution=3)
            thetalist = np.array([-3., 0, 3, 3, 3, 3, 3, 3, 3, 3, 3, 0, -3, -3, -3, -3, -3, -3, -3, -3, -3], dtype=np.float)
            philist = -70 + 20 * ch + np.array([0., 0, 0, 2.5, 5, 7.5, 10, 12.5, 15, 17.5, 20, 20, 20, 17.5, 15, 12.5, 10, 7.5, 5, 2.5, 0], dtype=np.float)
            print('T', thetalist)
            print('P', philist)
            jnacoords = np.array([frame.angles2jna(t, p) for t, p in zip(thetalist, philist)])
#            pepcoords = np.array([frame.jna2nsc(vec) for vec in jnacoords])
            pepcoords = frame.jna2nsc(jnacoords.T).T
            ganycoords = sc.convert_to_ref(pepcoords.T).T
            xsect = np.ma.ones(ganycoords.shape) * np.nan
            for idx, vec in enumerate(ganycoords):
                v2 = pyana.util.intersection.los_limb_sphere(scpos_gany, vec, [0, 0, 0], rg)
                if v2 != None:
                    xsect[idx, :] = np.array([v2.x, v2.y, v2.z])

            xsect = np.ma.masked_invalid(xsect)

            xsectn = xsect / rg

            lat = np.arcsin(xsectn[:, 2].compressed()) * 180. / np.pi
            lon = np.unwrap(np.arctan2(xsectn[:, 1].compressed(), xsectn[:, 0].compressed()))
            print(lon)
            lon = lon * 180. / np.pi

            x, y = m(lon, lat)
            m.plot(x, y, 'r')
        if onmap:
            m.warpimage('../../pyana/pep/images/ganymede.jpg')



    return fig


def main():

    tbl = jspice.JuiceSummary.table

    dt = datetime.timedelta(seconds=10)
    dt2 = datetime.timedelta(minutes=2)
    dt3 = datetime.timedelta(minutes=10)

    # 500 km orbit.
    t0_500 = tbl['GCO-500'][1]  # Start
    t0_200 = tbl['GCO-200'][1]  # End

#    t500 = t0_500 + (t0_200 - t0_500) / 3  # Average is in 500 km orbit.
    t500 = datetime.datetime(2033, 4, 7, 3)


    # 200 km orbit
    t0_end = tbl['END'][1]   # Mission end
    t200 = t0_200 + (t0_end - t0_200) / 2  # Average is in 500 km orbit.
    print(t200)

    t200 = datetime.datetime(2033, 6, 17, 14,  0)

    fig = plotonmap(t500, t500 + datetime.timedelta(hours=10), dt, ticktlist=t500 + dt3 * np.arange(60), fovtlist=t500 + dt2 * np.arange(0, 10),
                onmap=True)
    fig.savefig('traj_on_surface_500km.png')
    fig = plotonmap(t200, t200 + datetime.timedelta(hours=10), dt, ticktlist=t200 + dt3 * np.arange(60), fovtlist=t200 + dt2 * np.arange(0, 10),
                onmap=True)
    fig.savefig('traj_on_surface_200km.png')


if __name__ == "__main__":
    main()