apps120803_jdcfly.app03_sw_scpos

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 get_dflux(fv, scpos, scvel):
    """Return the differential flux.  Replaced by app05_op_longitude.get_dflux2.

    Differential flux measured from the spacecraft will be returned.

    The spacecraft is assumed to have a nadir-pointing mode.

    :param scpos: Spacecraft position, used to convert to nadir pointing sensor
    :param scvel: Spacecraft velocity, used to convert to nadir pointing sensor
    :param fv: A function that returns the velocity distribution function.
    :returns: Differential flux for each energy, elevation and azimuthal channel.
        This is a tuple. (vveclist, dflux).
        vveclist is the velocity vector in the reference (physical) frame.
        dflux is the differential flux.
        Each has the shape of (nEne=128, nEl=32, nAz=16).
    """

    import warnings
    warnings.warn('The module app03_sw_scpos.get_dflux is deprected. Replaced by app05_op_longitude.get_dflux2', DeprecationWarning)

    # S/C attitude
    nsc = att.NadirLookingSc()
    nsc.set_posvel(scpos, scvel)

    ## 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 = fv(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)

    return vveclist, dflux   # Returns is a tuple 
    
def main():
    n = 5e6
    v = (-400.e3, 0, 0)
    vth = 100e3

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

    postheta = 30. * np.pi / 180.
    posphi = 120. * np.pi / 180.

    pos = [np.cos(postheta) * np.cos(posphi), np.cos(postheta) * np.sin(posphi), np.sin(postheta)]
    vel = [-np.sin(postheta) * np.cos(posphi), -np.sin(postheta) * np.sin(posphi), np.cos(postheta)]
    vveclist, dflux = get_dflux(fmaxwell, pos, vel)

    enelist = ((vveclist / 13800.) ** 2).sum(0)

#    # Choose 1000 point to plot in 3D.
    try:
        from mayavi import mlab
#        randomchoice = (np.random.random_sample(1000) * 128 * 32 * 16).astype(np.int)

        ### Choose only important data.
#        idx0 = (vveclist * vveclist).sum(0) < (1e6 ** 2)  # <1000 km/s
#        idx1 = (vveclist * vveclist).sum(0) > (1e4 ** 2)  # > 10 km/s
        vveclist.shape = (3, 128 * 32 * 16)
        dflux = dflux.flatten()

        randomchoice = np.where(dflux > 1e+3)[0]

        vveclist = vveclist[:, randomchoice] / 1e3
        dflux = np.log10(dflux[randomchoice])
        print(vveclist.shape, dflux.shape)
        mlab.points3d(vveclist[0, :], vveclist[1, :], vveclist[2, :], dflux, vmin=-3)
        mlab.axes(extent=[-500,500,-500,500,-500,500])
        mlab.outline(extent=[-1000, 1000, -1000, 1000, -1000, 1000])

    except ImportError:
        print("NO 3D OUTPUT")


if __name__ == "__main__":
    main()