apps120803_jdcfly.app05_op_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. But for this case, we try to use oxygen.
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.
''' 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.
But for this case, we try to use oxygen.
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/app05_op_longitude_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.pep.pep_attitude as att
#from app03_sw_scpos import get_dflux
def get_dflux2(fv, scpos, scvel, m=1.67e-27, q=1.60e-19):
"""Return the differential flux. Replaces app03_sw_scpos.get_dflux.
Differential flux measured from the spacecraft will be returned.
The spacecraft is assumed to have a nadir-pointing mode.
RAM velocity is not considered.
: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. MKSA unit. (s^3/m^6)
:keyword m: Mass. MKSA unit. (kg)
:keyword q: Charge. MKSA unit. (C)
: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).
The unit is #/cm2 sr s eV.
"""
# S/C attitude
nsc = att.NadirLookingSc()
nsc.set_posvel(scpos, scvel)
## Now setup JDC.
# Direction
enlist =np.arange(128)
enevlist = energy.getEnergy() # in eV
vellist = np.sqrt(enevlist) * np.sqrt(2 * q / m) # 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,)
## 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)
vlist2 = (vveclist ** 2).sum(0) # (128, 32, 16)
print('fmax', fval.max())
dflux = vlist2 * fval / m
dflux = dflux / 1e4 * q # To convert to #/cm2 sr eV s, (128, 32, 16)
print('jmax', dflux.max())
return vveclist, dflux # Returns is a tuple
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 = 6e+0
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)
print('Max dflux=', dflux.max())
#### 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))
print(dfluxall[iax].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('app05_op_longitude_%g.png' % phideg)
if __name__ == "__main__":
main()