apps120718_jna_fov.map_on_surface

''' A simple script to create a map on the surface.

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

import math
import numpy as np

import matplotlib as mpl
import matplotlib.pyplot as plt

import pyana.util.intersection
import pyana.util.triangle2d

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

import pyana.pep.pep_attitude

if __name__ == '__main__':

    # Define the spacecraft height.
    h = 500.

    # First, calculate the FoVs and their FWHMs.

    azim_center = [fov.azim_pix_center(i) for i in range(7)]
    azim_fwhm = [fov.azim_pix_fwhm(i) for i in range(7)]

    # In the elevation direction
    elev_center = [fov.elev_pix_center(i) for i in range(7)]
    elev_fwhm = [fov.elev_pix_fwhm(i) for i in range(7)]

    # Prepare the canvas
    fig = plt.figure()
    ax = fig.add_subplot(111)

    # Some inertial ganymede frame
    scpos_gany = np.array([2800, 0, 0])
    scvel_gany = [0, 0, -2.4]
    scpos_norm = scpos_gany / np.sqrt((scpos_gany ** 2).sum())

    rg = 2634.1

    sc = pyana.pep.pep_attitude.NadirLookingSc()
    sc.set_posvel(scpos_gany, scvel_gany)

    colors = [None, 'b', 'r', 'k', 'g', 'm', None]
    for ch in range(0, 7):
        thetalist, philist = fov.pixel_shape(ch, azimresolution=10, elevresolution=3)
        print(thetalist[0], philist[0])

        jnacoords = np.array([frame.angles2jna(t, p) for t, p in zip(thetalist, philist)])
        print(jnacoords[0])
        
        pepcoords = np.array([frame.jna2nsc(vec) for vec in jnacoords])
        print(pepcoords[0])

        ganycoords = sc.convert_to_ref(pepcoords.T).T
        print(ganycoords[0])

        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

        ax.plot(lon, lat, 'o')

        polyarea = pyana.util.triangle2d.Polygon2D(lon, lat)
        triangles = polyarea.triangulate()

        for trian in triangles:
            ax.fill([trian.p0[0], trian.p1[0], trian.p2[0], trian.p0[0]],
                [trian.p0[1], trian.p1[1], trian.p2[1], trian.p0[1]], alpha=0.1, color='b')


        print(len(triangles), len(lon))

    ax.set_aspect('equal')
    fig.savefig('./map_on_surface.png')