apps120803_jdcfly.app02_sw

A very simple exercise to observe Maxwell distributed SW by JCD.

The second lesson is to implement JDC virtual observation of Maxwell distribution. Maxwell should have a bulk velocity of 400 km/s in the -x direction 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 distributed SW by JCD.

The second lesson is to implement JDC virtual observation of Maxwell distribution.
Maxwell should have a bulk velocity of 400 km/s in
the -x direction 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 = (-400.e3, 0, 0)
    vth = 100e3

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

    # S/C attitude
    nsc = att.NadirLookingSc()
    nsc.set_posvel([2500., 0, 0], [0, 0, -1.5])

    ## Now setup JDC.
    # Direction
    enlist =np.arange(128)
    enevlist = energy.getEnergy()
    vellist = np.sqrt(enevlist) * 13800  # Assume protons, in m/s
    ellist = np.arange(32)
    azlist = np.arange(16)

    fval = np.zeros([len(enlist), len(ellist), len(azlist)])
    vveclist = np.zeros([3, len(enlist), len(ellist), len(azlist)])

    for iel in ellist:
        for iaz in azlist:
            el = fov.elev_pix_center(iel)
            az = fov.azim_pix_center(iaz)
            ## Looking direction.
            lookingdir_jdc = frame.angles2jdc(el, az)  # (3, 512)
            ## Looking direction in SC
            lookingdir_nsc = frame.jdc2nsc(lookingdir_jdc)
            # conversion to physical reference.
            lookingdir_ref = nsc.convert_to_ref(lookingdir_nsc)
            # The looking dir and observing velocity vector is opposite.
            velocitydir_ref = -lookingdir_ref

            for iene in enlist:

                vel = vellist[iene]
                vvec = velocitydir_ref * vel

                f = fmaxwell(vvec)

                fval[iene, iel, iaz] = f
                vveclist[:, iene, iel, iaz] = vvec

    vveclist = np.array(vveclist)   # (3, 128, 32, 16), m/s
    fval = np.array(fval)  # (128, 32, 16)
    enelist = ((vveclist / 13800.) ** 2).sum(0)   # (128, 32, 16), eV
    print(enelist.shape)

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

    plt.xlim(1, 10000)
    plt.ylim(1e-4, 1e5)

#    # Choose 1000 point to plot.
    try:
        from mayavi import mlab
        randomchoice = (np.random.random_sample(1000) * 128 * 32 * 16).astype(np.int)
        vveclist.shape = (3, 128 * 32 * 16)
        vveclist = vveclist[:, randomchoice] / 1e3
        dflux = np.log10(dflux.flatten()[randomchoice])
        print(vveclist.shape, dflux.shape)
        mlab.points3d(vveclist[0, :], vveclist[1, :], vveclist[2, :], dflux, vmin=-5)

        plt.figure()
        plt.plot(vveclist[0, :], dflux)
    except ImportError:
        print("NO 3D OUTPUT")

if __name__ == "__main__":
    main()