''' A very simple exercise to observe Maxwell distributed SW by JCD over half the orbit.
The next lesson is to have the data in count rate.
Now it is time to plot using JIA's data.
:mod:`irfpy.pep.mhddata` support reading and handling the data.
For simplicity, no interface is prepared.
.. image:: ../../../src/scripts/apps120803_jdcfly/app07_jdc_flyover.py_1.2_0.png
:width: 80%
.. image:: ../../../src/scripts/apps120803_jdcfly/app07_jdc_flyover.py_1.2_180.png
:width: 80%
'''
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
import irfpy.pep.mhddata
from app05_op_longitude import get_dflux2
def main(pp=None, longitude=180., r=1.2):
## Orbit around Ganymede.
ndiv = 19
posthetas = np.linspace(-np.pi / 2., np.pi / 2., ndiv) # Only dayside
posphi = longitude * 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
if pp == None:
pp = irfpy.pep.mhddata.PlasmaParameter1205()
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)])
pos = pos * r
n, vx, vy, vz, temp, pth = pp.interpolate3d(pos[0], pos[1], pos[2])
vth = irfpy.pep.mhddata.t2vth(temp, mass=16.) # in m/s
n *= 1e6 # in /m3
vx *= 1e3; vy *= 1e3; vz *= 1e3 # in m/s
fmaxwell = mkfunc(n, [vx, vy, vz], vth) # maxwell distribution function in physical coord. No mass dependence.
vel = [-np.sin(postheta) * np.cos(posphi), -np.sin(postheta) * np.sin(posphi), np.cos(postheta)]
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' % longitude)
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('app07_jdc_flyover.py_%g_%g.png' % (r, longitude))
return pp
if __name__ == "__main__":
pp = main()
main(pp=pp, longitude=0)