apps120803_jdcfly.app01_maxwell

A very simple exercise to observe Maxwell distribution by JCD.

The first lesson is to implement JDC virtual observation of Maxwell distribution. Maxwell should have a bulk velocity of 0 with a certain density, say 5 part/cm3, and the thermal velocity of 100 km/s.

The Maxwell is implemented in irfpy.util.maxwell.

''' A very simple exercise to observe Maxwell distribution by JCD.

The first lesson is to implement JDC virtual observation of Maxwell distribution.
Maxwell should have a bulk velocity of 0 with a certain density, say 5 part/cm3,
and the thermal velocity of 100 km/s.

The Maxwell is implemented in :mod:`irfpy.util.maxwell`.
'''

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

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

from irfpy.util.maxwell import mkfunc

import irfpy.jdc.energy0 as energy
import irfpy.jdc.fov0 as fov
import irfpy.jdc.frame0 as frame
import irfpy.pep.pep_attitude as att

def main():
    '''Main script'''

    n = 5e6
    v = (0., 0, 0)
    vth = 100e3

    fmaxwell = mkfunc(n, v, vth)  # maxwell distribution function.

    ## Now setup JDC.
    # Direction
    ellist = np.arange(32)
    azlist = np.arange(16)
    azs2, els2 = np.meshgrid(azlist, ellist)  # els stores elevation value for (32, 16)
    azs = azs2.flatten()  # Now (512)
    els = els2.flatten()  # Now (512)

    ## Looking direction.
    azangles = fov.azim_pix_center(azs)
    elangles = fov.elev_pix_center(els)
    lookingdir_jdc = frame.angles2jdc(elangles, azangles)  # (3, 512)

    ## Looking direction in SC
    lookingdir_nsc = frame.jdc2nsc(lookingdir_jdc)

    # conversion to physical reference.
    nsc = att.NadirLookingSc()
    nsc.set_posvel([2500., 0, 0], [0, 0, -1.5])
    lookingdir_ref = nsc.convert_to_ref(lookingdir_nsc)

    velocitydir_ref = -lookingdir_ref   # The looking dir and observing velocity vector is opposite.

#    ## To display date points in 3D.
#    from mayavi import mlab
#    mlab.points3d(lookingdir_jdc[0, :], lookingdir_jdc[1, :], lookingdir_jdc[2, :], np.arange(512))
#    mlab.points3d(lookingdir_jdc[0, :], lookingdir_jdc[1, :], lookingdir_jdc[2, :], azs)
#    mlab.points3d(lookingdir_jdc[0, :], lookingdir_jdc[1, :], lookingdir_jdc[2, :], els)
    
    # For each looking direction, energy is calculated. 

    enelist = energy.getEnergy()
    vellist = np.sqrt(enelist) * 13800  # Assume protons

    flist = []
    vveclist = []

    for idir in range(len(azangles)):
        az = azs[idir]
        el = els[idir]
        vecref = velocitydir_ref[:, idir]
        for ene, vel in zip(enelist, vellist):
            vvec = vecref * vel
            f = fmaxwell(vvec)
            vveclist.append(vvec)
            flist.append(f)

    vveclist = np.array(vveclist)   # (65536, 3)
    flist = np.array(flist)  # (65536)
    enelist = ((vveclist / 13800.) ** 2).sum(1)   # (65536.) in eV

#    # Choose 1000 point to plot.
#    from mayavi import mlab
#    randomchoice = (np.random.random_sample(1000) * len(flist)).astype(np.int)
#    vveclist = vveclist[randomchoice, :]
#    flist = flist[randomchoice]
#    print vveclist.shape, flist.shape
#    mlab.points3d(vveclist[:, 0], vveclist[:, 1], vveclist[:, 2], np.log10(flist))

    dflux = 2 * enelist * (1.60e-19 / 1.67e-27) * flist / 1.67e-27
    dflux = dflux / 1e4 * 1.60e-19   # To convert to #/cm2 sr eV s
    
    import matplotlib.pyplot as plt
    plt.plot(enelist, dflux, '.')   # Regardless of direction, everything on 1 line
    plt.xscale('log')
    plt.yscale('log')

    plt.xlabel('Energy')
    plt.ylabel('Differential flux [/cm2 sr eV s]')


    plt.savefig('app01.png')

if __name__ == "__main__":
    main()