apps120803_jdcfly.app08_jdc_jna_flyover

A very simple exercise to observe Maxwell distributed SW by JCD/JNA over half the orbit.

''' A very simple exercise to observe Maxwell distributed SW by JCD/JNA over half the orbit.
'''

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

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

from irfpy.util.maxwell import mkfunc
import irfpy.util.intersection

import irfpy.jdc.energy0 as energy
import irfpy.jdc.fov0 as fov
import irfpy.jdc.frame0 as frame
import irfpy.jdc.flux0

import irfpy.jna.fov0 as jnafov
import irfpy.jna.frame0 as jnaframe
import irfpy.pep.pep_attitude as att
import irfpy.pep.mhddata

from .app05_op_longitude import get_dflux2

def main(pp=None, longitude=180., r=1.2, type="footprint"):

    ## Orbit around Ganymede.
    ndiv = 360
    posthetas = np.linspace(np.pi / 2., np.pi * 5 / 2., ndiv)  # Only dayside
    posphi = longitude * np.pi / 180.  # For the first trial, Noon-midnight meridian is taken.

    vmin = 4e-1
    vmax = 1e+1

    m = 16 * 1.67e-27
    q = 1.60e-19

    # For plotting
    fig = plt.figure(figsize=(15, 10))
    nullfmt = matplotlib.ticker.NullFormatter()

    # For JNA

    # Precipitation data
    precip = irfpy.pep.mhddata.PrecipitationFlux(type=type)
    lons = precip.reshape(precip.lonlist())
    lats = precip.reshape(precip.latlist())
    flx = precip.reshape(precip.flist())   # Flux at the surface.

    sc = att.NadirLookingSc()

    jnaflx = np.zeros((ndiv, 7)) + np.nan

    for ilat in range(ndiv):
        # For JNA
        postheta = posthetas[ilat]   # Position of S/C

        # S/C position
        pos = np.array([np.cos(postheta) * np.cos(posphi), np.cos(postheta) * np.sin(posphi), np.sin(postheta)])
        pos = pos * r
        # S/C Velocity direction.
        vel = np.array([-np.sin(postheta) * np.cos(posphi), -np.sin(postheta) * np.sin(posphi), np.cos(postheta)])
        # S/C instance
        sc.set_posvel(pos, vel)

        # Conversion.
        for ch in (1, 2, 3, 4, 5):  # Only center 5 channels.
            # Center direction
            az = jnafov.azim_pix_center(ch)
            el = jnafov.elev_pix_center(ch)
            vecjna = jnaframe.angles2jna(el, az)
            vecnsc = jnaframe.jna2nsc(vecjna)
            vecgan = sc.convert_to_ref(vecnsc)   # Looking vector in Ganymede coordinates.
            v2 = irfpy.util.intersection.los_sphere(pos, vecgan, [0, 0, 0], 1.)
            if v2 == None:
                continue
            longit = np.arctan2(v2[1], v2[0]) * 180. / np.pi
            latit = np.arcsin(v2[2]) * 180. / np.pi

            idxlon, idxlat = precip.nearest_neighbor(longit, latit)
            jnaflx[ilat, ch] = flx[idxlon, idxlat] * 0.2 / (2 * np.pi) * 2e-6   # in /s. Indeed, count rate
            # 0.2 comes from 20% from CENA result.
            # Sensor g-factor is taken from proposal v0.3, 2e-6 cm2 sr/pixel.
            print(jnaflx[ilat, ch])

    ax2 = fig.add_subplot(212)
    xax = np.linspace(0, 360, ndiv + 1)
    yax = np.arange(8)
    x, y = np.meshgrid(xax, yax)
    img = ax2.pcolor(x, y, np.ma.log10(jnaflx).T, vmin=np.log10(vmin), vmax=np.log10(vmax))
    ax2.set_xlim(0, 360)
    ax2.set_ylim(0, 7)
    ax2.set_xlabel('Angle (North pole=0)')
    ax2.set_ylabel('Direction')
    fig.colorbar(img)

    # FOR JDC
    vmin = 1e-1
    vmax = 6e+4

#    ndiv = 360   # High resolution
    ndiv = 8  # For debug

    # Now, this is to save the data.
    dfluxall = np.zeros([16, ndiv, 128]) - np.nan  # 16 is for number of panels


    if pp == None:
        pp = irfpy.pep.mhddata.PlasmaParameter1205()

    for ilat in range(ndiv):
        # For JDC
        postheta = posthetas[ilat]
        pos = np.array([np.cos(postheta) * np.cos(posphi), np.cos(postheta) * np.sin(posphi), np.sin(postheta)])
        pos = pos * r

        n, vx, vy, vz, temp, pth = pp.interpolate3d(pos[0], pos[1], pos[2])
        vth = irfpy.pep.mhddata.t2vth(temp, mass=16.)   # in m/s
        n *= 1e6   # in /m3
        vx *= 1e3; vy *= 1e3; vz *= 1e3  # in m/s

        fmaxwell = mkfunc(n, [vx, vy, vz], vth)  # maxwell distribution function in physical coord.  No mass dependence.

        vel = [-np.sin(postheta) * np.cos(posphi), -np.sin(postheta) * np.sin(posphi), np.cos(postheta)]
        vveclist, dflux = get_dflux2(fmaxwell, pos, vel, m=m)   # O+ assumede.
        enelist = ((vveclist / np.sqrt(2 * q / m)) ** 2).sum(0)

        #### dflux has (128, 32, 16) dimension.
        # panel0: el 24-31 (zenith), az 0-1,14-15 (RAM), max
        dfluxall[0, ilat, :] = dflux[:, 24:32, [0, 1, 14, 15]].max(2).max(1)


    ## Convert to count rate
    j2c = irfpy.jdc.flux0.Flux2Count()
    c = np.zeros_like(dfluxall)
    for ie in range(128):
        c[:, :, ie] = j2c.getCounts(dfluxall[:, :, ie], ie)


    ax = fig.add_subplot(2, 1, 1)

    xax = np.linspace(0, 360, ndiv + 1)
    yax = energy.getBound()

    x, y = np.meshgrid(xax, yax)

    ### Plot the data.
    img = ax.pcolor(x, y, np.ma.log10(c[0]).T, vmin=np.ma.log10(vmin), vmax=np.ma.log10(vmax))

    ax.set_yscale('log')
    ax.set_xlim(0, 360)
    ax.set_ylim(1, 41000)

    fig.colorbar(img)

    return pp


if __name__ == "__main__":
    pp = main(longitude=0)
    plt.savefig('app08_jdc_jna_flyover_1.png')
    pp = main(longitude=0, type='nadir')
    plt.savefig('app08_jdc_jna_flyover_2.png')