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.
''' 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()