apps120911_iotorus.app05_enaflux_plot

Just plot the app04_eneflux_batch.dat

''' Just plot the app04_eneflux_batch.dat
'''

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

import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
from matplotlib.patches import Circle
import numpy as np
import scipy as sp

def thetaphi2xyscreen(theta, phi):
    x = -theta * np.cos(np.deg2rad(phi))
    y = theta * np.sin(np.deg2rad(phi))
    return x, y

def thetaphi2xyzvector(theta_deg, phi_deg, ax0, ax1, ax2):
    t = np.deg2rad(theta_deg)
    p = np.deg2rad(phi_deg)

    lookdir = ax0 * np.cos(t) + ax1 * np.sin(t) * np.cos(p) + ax2 * np.sin(t) * np.sin(p)

def xyzvector2thetaphi(xyz, ax0, ax1, ax2):
    r = xyz / np.sqrt((xyz * xyz).sum())
    cos_t = np.dot(r, ax0)
    if cos_t > 1: cos_t = 1
    if cos_t < -1: cos_t = -1

    t = np.arccos(cos_t)
    x = np.dot(r, ax1)
    y = np.dot(r, ax2)
    p = np.arctan2(y, x)

    t = np.rad2deg(t)
    p = np.rad2deg(p)

    return t, p

def main(filename='app04_eneflux_batch.dat'):
    '''Main script'''

    f = open(filename)

    pos = f.readline().split()
    print(pos)
    ax0 = f.readline().split()
    ax1 = f.readline().split()
    ax2 = f.readline().split()

    # Spacecraft position

    xsc = float(pos[2])
    ysc = float(pos[3])
    zsc = float(pos[4])

    # Axis information
    ax0 = np.array([float(v) for v in ax0[2:]])
    ax1 = np.array([float(v) for v in ax1[2:]])
    ax2 = np.array([float(v) for v in ax2[2:]])

    print('POS', xsc, ysc, zsc)
    print('Ax0', ax0)
    print('Ax1', ax1)
    print('Ax2', ax2)

    # Read data.
    dats = f.readlines()

    # Data into array.
    ths = []   # Theta list.
    phs = []   # Phi list.
    dflxs = [] # Dflux list.

    lookdir_xyz = []  # Looking direction in x, y, z

    for dat in dats:
        datsp = dat.split()

        th = float(datsp[0])  # Polar angle
        ph = float(datsp[1])  # Azimuth angle
        ldir = np.array([float(d) for d in datsp[2:5]])
        dflx = float(datsp[5]) # Differential flux

        ths.append(th)
        phs.append(ph)
        dflxs.append(dflx)
        lookdir_xyz.append(ldir)

    ths = np.array(ths)
    phs = np.array(phs)
    dflxs = np.array(dflxs)
    lookdir_xyz = np.array(lookdir_xyz)

    # Map to a screen
    x, y = thetaphi2xyscreen(ths, phs)

    # Irregular data converted into the regular grid data.
    ximg = np.linspace(-50, 50, 720*2)
    yimg = np.linspace(-50, 50, 720*2)
    fimg = mlab.griddata(x, y, dflxs, ximg, yimg, interp='nn')

#    plt.scatter(x, y, np.log10(dflxs))
    cs = plt.contourf(ximg, yimg, np.log10(fimg), [0, 1, 2, 3, 4, 5, 6, 7, ], cmap=plt.cm.jet)
    cb = plt.colorbar(cs)
    cb.set_label('Log of ENA flux [/cm2 sr s]')
#    plt.contour(ximg, yimg, np.log10(fimg), np.arange(-5, 8), cmap=plt.cm.gray)
    
    print(max(dflxs))

    # Map the Jupiter (center)
    sc = np.array([xsc, ysc, zsc])

    jcenter = np.array([0, 0, 0])

    lookdir = jcenter - sc
#    lookdir = lookdir / np.sqrt((lookdir ** 2).sum())  # Normalize
    print('Jupiter', lookdir)
#    theta = np.rad2deg(np.arccos(-lookdir[1]))
#    phi = np.arctan2(lookdir[2], lookdir[0]) 
    theta, phi = xyzvector2thetaphi(lookdir, ax0, ax1, ax2)
    print('Jupiter', theta, phi)
    phi = np.deg2rad(phi)

    jupxy = [-theta * np.cos(phi), theta * np.sin(phi)]
#    plt.plot(jupxy[0], jupxy[1], 'ro')

    # Jupiter size
    js = np.rad2deg(np.arcsin(1 / np.sqrt(np.dot(sc, sc))))
    circ = Circle(jupxy, js, edgecolor='none', facecolor='gray')
    print('J size', js)
    plt.gca().add_patch(circ)

    
    # Map the Io, Europa, Ganymede orbit?
    r = [5.9, 9.0, 14.9]
    for rm in r:
        phases = np.linspace(0, 2 * np.pi, 360)
        xyio = np.zeros([2, 360])
        for iphase, phase in enumerate(phases):
            iopos = np.array([np.cos(phase), np.sin(phase), 0]) * rm   # in Io frame
            lookdir = iopos - sc   # In Io frame
            len_lookdir = np.sqrt((lookdir ** 2).sum())
            lookdir = lookdir / len_lookdir   # Normalized.

            theta, phi = xyzvector2thetaphi(lookdir, ax0, ax1, ax2)
            phi = np.deg2rad(phi)
#            print theta, phi, '->', 
#
#            theta = np.rad2deg(np.arccos(-lookdir[1]))
#            phi = np.arctan2(lookdir[2], lookdir[0])
#            print theta, phi
            xyio[:, iphase] = [-theta * np.cos(phi), theta * np.sin(phi)]

            if phase == 0 and rm == 5.9:
                mpos = [-theta * np.cos(phi), theta * np.sin(phi)]
#                plt.plot(mpos[0], mpos[1], 'go')
                circm = np.rad2deg(np.arcsin(1821. / 71500 / len_lookdir))
                print(circm)
                circ = Circle(mpos, circm, facecolor='g', edgecolor='none')
                plt.gca().add_patch(circ)
                
        plt.plot(xyio[0, :], xyio[1, :], 'g:', lw=0.5)

    # FoV, 25x7
    t = np.array([3.5, 3.5, 3.5, 0, -3.5, -3.5, -3.5, 0, 3.5])
    p = np.array([12.5, 0, -12.5, -12.5, -12.5, 0, 12.5, 12.5, 12.5])

    plt.fill(t, p, edgecolor='none', alpha=0.5, color='purple')
    
    plt.gca().set_xticks([-90, -60, -30, 0, 30, 60, 90])
    plt.gca().set_yticks([-90, -60, -30, 0, 30, 60, 90])
    plt.gca().set_xticklabels(['', '', '', '', '', '', '', ''])
    plt.gca().set_yticklabels(['', '', '', '', '', '', '', ''])
    plt.gca().set_aspect('equal')
    plt.gca().set_xlim(-50, 50)
    plt.gca().set_ylim(-50, 50)

    plt.title('ENA flux (Io neutral cloud): 100-200eV')

    plt.savefig('app05_enaflux_plot.png')
    plt.savefig('app05_enaflux_plot.eps')

if __name__ == "__main__":
    try:
        filename = sys.argv[1]
    except:
        filename = 'app04_eneflux_batch.dat'
    
    main(filename=filename)