apps120820_joee.joee_pitchangle_bound

''' Time series of pitch angle coverage


.. image:: ../../../src/scripts/apps120820_joee/joee_pitchangle_bound_mhd.png
    :width: 90%

'''
import datetime

import numpy as np
import matplotlib.pyplot as plt

import pyana.pep.juice_spice1205 as spice
import pyana.joee.fov0 as fov
import pyana.joee.frame0 as frame
from pyana.pep.mhddata import MagneticField1201 as MF
import pyana.pep.pep_attitude

import pyana.pep.ganymede_magfield

def geo2pla(vecs_geo):
    vecs_pla = np.ma.zeros(vecs_geo.shape)
    vecs_pla[0] = -vecs_geo[1]
    vecs_pla[1] = vecs_geo[0]
    vecs_pla[2] = vecs_geo[2]
    return vecs_pla


def main(model='mhd'):
    t0 = datetime.datetime(2033, 4, 4, 18, 50)
    t1 = t0 + datetime.timedelta(minutes=180)
    dt = datetime.timedelta(seconds=30)
#    dt = datetime.timedelta(seconds=180)

    tlist = [t0]
    while tlist[-1] < t1:
        tlist.append(tlist[-1] + dt)
    print(len(tlist))

    if model == 'dipole':
        mf = pyana.pep.ganymede_magfield.TiltedDipole()
    else:
        mf = MF(interpolate='trilin')

    juice = spice.JuiceSpice.get_default_instance()
    rg = 2634.1
    sc = pyana.pep.pep_attitude.NadirLookingSc()

    pitchangles = np.zeros([16, len(tlist)]) + np.nan
    nadirangles = np.zeros_like(pitchangles) + np.nan
    lons = np.zeros_like(tlist) + np.nan
    lats = np.zeros_like(tlist) + np.nan
    bs = np.zeros([3, len(tlist)]) + np.nan

    pitchbound = np.zeros([2, len(tlist)]) + np.nan

    for idx, t in enumerate(tlist):
        print(idx)
        # Position and velocity in plasma frame.
        position = geo2pla(juice.get_position(t, relative_to="Ganymede", frame="IAU_GANYMEDE"))
        velocity = geo2pla(juice.get_velocity(t, relative_to="Ganymede", frame="IAU_GANYMEDE"))
        print('POS,VEL=', position, velocity)

        h = np.sqrt((position ** 2).sum())
        lon = np.arctan2(position[1] / h, position[0] / h)  # In rad
        lat = np.arcsin(position[2] / h)  # In rad

        lons[idx] = np.rad2deg(lon)
        lats[idx] = np.rad2deg(lat)

        print('LON,LAT', lons[idx], lats[idx])

        pos = position / rg   # In Rg
        posnorm = position / h

        mag = mf.interpolate3d(pos[0], pos[1], pos[2])   # In nT, (3,)
        bs[:, idx] = mag
        magmag = np.sqrt((mag ** 2).sum())
        magnorm = mag / magmag

        sc.set_posvel(position, velocity)

        # Now start conversions
        for ch in np.arange(16):    # Each channel
            sc.set_posvel(position, velocity)
            lookdir_el = pyana.joee.fov0.elev_pix_center(ch)
            lookdir_az = pyana.joee.fov0.azim_pix_center(ch)
            vecjoee = pyana.joee.frame0.angles2joee(lookdir_el, lookdir_az)
            vecnsc = pyana.joee.frame0.joee2nsc(vecjoee)
            vecganpla = sc.convert_to_ref(vecnsc)   # Looking direciton of the certain channel.

            angle_from_nadir = np.arccos((-vecganpla * posnorm).sum()) * 180. / np.pi
            angle_from_b = np.arccos((vecganpla * magnorm).sum()) * 180. / np.pi

            pitchangles[ch, idx] = 180 - angle_from_b   # Now the velocity vector with 180-
            nadirangles[ch, idx] = angle_from_nadir
        
        # Ok, same for -0.5 and 7.5
        for ich, ch in enumerate([-0.5, 7.5]):
            sc.set_posvel(position, velocity)
            lookdir_el = pyana.joee.fov0.elev_pix_center(ch)
            lookdir_az = pyana.joee.fov0.azim_pix_center(ch)
            vecjoee = pyana.joee.frame0.angles2joee(lookdir_el, lookdir_az)
            vecnsc = pyana.joee.frame0.joee2nsc(vecjoee)
            vecganpla = sc.convert_to_ref(vecnsc)   # Looking direciton of the certain channel.
            angle_from_nadir = np.arccos((-vecganpla * posnorm).sum()) * 180. / np.pi
            angle_from_b = np.arccos((vecganpla * magnorm).sum()) * 180. / np.pi            
            pitchbound[ich, idx] = 180 - angle_from_b   # Now the velocity vector with 180-             nadirangles[ch, idx] = angle_from_nadir
            

        
    fig = plt.figure(figsize=(8, 14))
    ax1 = fig.add_subplot(411)
    for ch in np.arange(16):
        ax1.plot(tlist, pitchangles[ch, :], 'b:', label='Ch%02d' % ch)
    for ch in range(2):
        ax1.plot(tlist, pitchbound[ch, :], 'b')

    ax1.set_title('JoEE Pitch angle coverage')
    ax1.set_ylabel('Pitch angle')
    ax1.yaxis.set_ticks([0, 45, 90, 135, 180])
#    ax1.legend(loc='best', prop={'size': 7})
    ax1.grid()

if __name__ == "__main__":
    retval = main('mhd')
#    plt.savefig('joee_pitchangle_bound_mhd.png')
    plt.savefig('joee_pitchangle_bound_mhd.eps')