apps120803_jdcfly.app04_sw_longitude

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

The next 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.

Now it is time to rotate the spacecraft around the Ganymede. The half rotation can be simulated. You can specify the meridan by “phideg” parameter in this script. ndiv provides the number of observations, i.e. equivalent to the time resolution.

For simplicity, no interface is prepared.

''' A very simple exercise to observe Maxwell distributed SW by JCD over 1 orbit.

The next 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.

Now it is time to rotate the spacecraft around the Ganymede.
The half rotation can be simulated.
You can specify the meridan by "phideg" parameter in this script.
``ndiv`` provides the number of observations, i.e. equivalent to
the time resolution.

For simplicity, no interface is prepared.

.. image:: ../../../src/scripts/apps120803_jdcfly/app04_sw_longitude_30.png
    :width: 90%
'''

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.jdc.energy0 as energy
import irfpy.jdc.fov0 as fov
import irfpy.jdc.frame0 as frame
import irfpy.pep.pep_attitude as att

from .app03_sw_scpos import get_dflux


def main():
    n = 5e6
    v = (-400.e3, 0, 0)
    vth = 100e3

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

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

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

    vmin = 1e+1
    vmax = 3e+5

    for ilat in range(ndiv):
        postheta = posthetas[ilat]

        pos = np.array([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)]
        print(postheta * 180. / np.pi, pos, vel)
        vveclist, dflux = get_dflux(fmaxwell, pos, vel)
        enelist = ((vveclist / 13800.) ** 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)
        # panel1: el 24-31, az 2-5, max
        dfluxall[1, ilat, :] = dflux[:, 24:32, 2:6].max(2).max(1)
        # panel2: el 24-31, az 6-9, max
        dfluxall[2, ilat, :] = dflux[:, 24:32, 6:10].max(2).max(1)
        # panel3: el 24-31, az 10-13, max
        dfluxall[3, ilat, :] = dflux[:, 24:32, 10:14].max(2).max(1)
        # panel4: el 16-23, az 0-1,14-15, max
        dfluxall[4, ilat, :] = dflux[:, 16:24, [0, 1, 14, 15]].max(2).max(1)
        dfluxall[5, ilat, :] = dflux[:, 16:24, 2:6].max(2).max(1)
        dfluxall[6, ilat, :] = dflux[:, 16:24, 6:10].max(2).max(1)
        dfluxall[7, ilat, :] = dflux[:, 16:24, 10:14].max(2).max(1)
        # panel8: el 8-15, az 0-1,14-15, max
        dfluxall[8, ilat, :] = dflux[:, 8:16, [0, 1, 14, 15]].max(2).max(1)
        dfluxall[9, ilat, :] = dflux[:, 8:16, 2:6].max(2).max(1)
        dfluxall[10, ilat, :] = dflux[:, 8:16, 6:10].max(2).max(1)
        dfluxall[11, ilat, :] = dflux[:, 8:16, 10:14].max(2).max(1)
        # panel12: el 0-7, az 0-1,14-15, max
        dfluxall[12, ilat, :] = dflux[:, 0:8, [0, 1, 14, 15]].max(2).max(1)
        dfluxall[13, ilat, :] = dflux[:, 0:8, 2:6].max(2).max(1)
        dfluxall[14, ilat, :] = dflux[:, 0:8, 6:10].max(2).max(1)
        dfluxall[15, ilat, :] = dflux[:, 0:8, 10:14].max(2).max(1)

#    dfluxall = np.ma.masked_less(dfluxall, vmin)

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

    # Colorbar.
    ax = fig.add_axes([0.95, 0.1, 0.02, 0.79])
    xx = [0, 1]
    yy = np.logspace(np.log10(vmin), np.log10(vmax), 257)
    xX, yY = np.meshgrid(xx, yy)
    ax.pcolor(xX, yY, np.log10(yy[np.newaxis, :]).T)
    ax.xaxis.set_major_formatter(nullfmt)
    ax.set_ylim(vmin, vmax)
    ax.set_ylabel('Particle differential flux [#/cm2 sr eV s]')
    ax.set_yscale('log')
    ax.set_title('Phi = %g' % phideg)

    left = [0.1, 0.3, 0.5, 0.7]
    bottom = [0.7, 0.5, 0.3, 0.1]
    width = 0.19
    height = 0.19

    axs = [fig.add_axes([left[i%4], bottom[i/4], width, height]) for i in range(16)]

    xax = np.linspace(-90, 90, ndiv + 1)
    yax = energy.getBound()

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

    for iax in range(16):
        ### Plot the data.
        img = axs[iax].pcolor(x, y, np.log10(dfluxall[iax].T), vmin=np.log10(vmin), vmax=np.log10(vmax))

        axs[iax].set_yscale('log')
        axs[iax].set_xlim(-90, 90)
        axs[iax].set_ylim(1, 41000)

        if iax < 12:
            axs[iax].xaxis.set_major_formatter(nullfmt)
        else:
            axs[iax].set_xlabel('Latitude')

        if iax % 4 != 0:
            axs[iax].yaxis.set_major_formatter(nullfmt)
        else:
            axs[iax].set_ylabel('Energy')

        axs[0].set_title('RAM looking')
        axs[1].set_title('Left looking')
        axs[2].set_title('Anti-RAM looking')
        axs[3].set_title('Right looking')

        axs[0].text(-85, 1e4, 'Zenith looking', color='white')
        axs[12].text(-85, 1, 'Nadir looking', color='white')


    fig.savefig('app04_sw_longitude_%g.png' % phideg)

if __name__ == "__main__":
    main()