apps120803_jdcfly.app06_opobs_alongorbit

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

The next lesson is to have the data in count rate.

Again, Maxwell should have a bulk velocity of 200 km/s in the -x direction with a certain density, say 5 part/cm3, and the thermal velocity of 50 km/s.

As in app04, we 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.

scripts/../../../src/scripts/apps120803_jdcfly/app06_opobs_alongorbit_0.png
''' A very simple exercise to observe Maxwell distributed SW by JCD over 1 orbit.

The next lesson is to have the data in count rate.

Again, Maxwell should have a bulk velocity of 200 km/s in
the -x direction with a certain density, say 5 part/cm3,
and the thermal velocity of 50 km/s.

As in app04, we 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/app06_opobs_alongorbit_0.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.jdc.flux0
import irfpy.pep.pep_attitude as att

from .app05_op_longitude import get_dflux2

def main():
    n = 5e6
    v = (-200.e3, 0, 0)   # corresponding to ~16 keV
    vth = 50e3

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

    ## Orbit around Ganymede.
    ndiv = 19
    posthetas = np.linspace(-np.pi / 2., np.pi / 2., ndiv)  # Only dayside
    phideg = 0.
    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 = 6e+4

    m = 16 * 1.67e-27
    q = 1.60e-19

    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_dflux2(fmaxwell, pos, vel, m=m)   # O+ assumede.
        enelist = ((vveclist / np.sqrt(2 * q / m)) ** 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)

    ## Convert to count rate
    j2c = irfpy.jdc.flux0.Flux2Count()
    c = np.zeros_like(dfluxall)
    for ie in range(128):
        c[:, :, ie] = j2c.getCounts(dfluxall[:, :, ie], ie)

    # 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('Count rate [c/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(c[iax]).T, vmin=np.log10(vmin), vmax=np.log10(vmax))

        print(c[iax, :, 97].max())

        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('app06_opobs_alongorbit_%g.png' % phideg)

if __name__ == "__main__":
    main()