import os
import sys
import logging
logging.basicConfig()
import datetime
import dateutil.parser
import math
from optparse import OptionParser
import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
import pyana.util.datafile
from pyana.util import dipole
from . import parameters as p
from . import observer
from . import jupiter
from . import frames
def main_versionized(filename, outname):
rj = p.rj
fp = open(filename)
effective_crosssection = 1.88 # Reference Xsection 1e-15. For O, 18.8e-16.
plasma_fraction = 0.41 # Density fraction of O+ and O++
efficiency = effective_crosssection * plasma_fraction
dfile = pyana.util.datafile.Datafile()
dfile.readheader(fp) # Only header is used. Data part is in the following
drfile = pyana.util.datafile.DatafileReader(dfile, missing_return="Unknown")
scpos_km = np.array(drfile.get_header('SCPOS').split()[-3:], dtype=np.float)
scpos = scpos_km / rj
print(scpos)
scvel = np.array(drfile.get_header('SCVEL').split()[-3:], dtype=np.float)
print(scvel)
ll = drfile.get_header('MASS, ENAVEL(SCf)').split()
velena = float(ll[-1])
try:
mass = float(ll[-2])
except:
mass = 16.
print(velena, mass)
ll = np.array(drfile.get_header('SUN_SYS3, IO').split(), dtype=np.float)
iopos = ll[-3:] / rj
try:
sunsys3lon = ll[-4]
except:
sunsys3lon = np.nan
print(sunsys3lon, iopos)
pax = np.array(drfile.get_header('PAX').split()[-3:], dtype=np.float)
sax = np.array(drfile.get_header('SAX').split()[-3:], dtype=np.float)
tax = np.array(drfile.get_header('TAX').split()[-3:], dtype=np.float)
print(pax, sax, tax)
maxangle = float(drfile.get_header('MAXANGLE').split()[-1])
print(maxangle)
sc = observer.VirtualSpacecraft(scpos, scvel, primary_axis=pax, secondary_axis=sax)
vals = np.loadtxt(fp) # (N, 9) array
print(vals.shape)
xcoords = np.unique(vals[:, 0]).tolist()
ycoords = np.unique(vals[:, 1]).tolist()
print(ycoords)
grid = np.zeros([len(xcoords), len(ycoords)])
xval = np.zeros_like(grid)
yval = np.zeros_like(grid)
fval = np.zeros_like(grid)
jval = np.zeros_like(grid)
for i in range(len(vals[:, 0])):
ximg, yimg, theta, phi, nx, ny, nz, f, j = vals[i, :]
idxx = xcoords.index(ximg)
idxy = ycoords.index(yimg)
# print idxx, idxy
xval[idxx, idxy] = ximg
yval[idxx, idxy] = yimg
fval[idxx, idxy] = f
jval[idxx, idxy] = j
jval = jval * efficiency
fig = plt.figure(figsize=(6, 6))
fig2 = plt.figure(figsize=(6, 6))
# A orbit information
ax1 = fig.add_subplot(111)
ax1.plot([0], [0], 'k+') # Jupiter
ax1.text(0, 0, 'J') # Jupiter
ax1.plot([scpos[0]], [scpos[1]], 'k*')
ax1.text(scpos[0], scpos[1], 'SC')
ax1.plot([iopos[0]], [iopos[1]], 'k+')
ax1.text(iopos[0], iopos[1], 'Io')
ax1.set_xlabel('X [Rj]')
ax1.set_ylabel('Y [Rj]')
ax1.set_xlim(-30, 30)
ax1.set_ylim(-30, 30)
ax = fig2.add_subplot(111)
# Convert "on-grid" value to "center of grids" value for pcolor.
xcoords = np.array(xcoords)
ycoords = np.array(ycoords)
xc = (xcoords[:-1] + xcoords[1:]) / 2. # Center of grids.
yc = (ycoords[:-1] + ycoords[1:]) / 2. # Center of grids.
img = ax.pcolor(xc, yc, np.log10(jval[1:-1, 1:-1]).T, vmin=0, vmax=6)
# cb = fig2.colorbar(img)
# cb.set_label('Log10 differential flux [/cm2 eV sr s]')
### Io position in x-y screen
iorelpos = iopos - scpos
th, ph = sc.get_angles(iorelpos)
th = np.rad2deg(th)
ximg = th * np.cos(ph)
yimg = th * np.sin(ph)
ax.plot([ximg], [yimg], 'ow')
r = np.rad2deg(np.arcsin(p.rio / p.rj / np.linalg.norm(iorelpos)))
print('Radius = ', r)
ax.plot(r * np.cos(np.linspace(0, 2 * np.pi, 45)),
r * np.sin(np.linspace(0, 2 * np.pi, 45)), 'w') # Wrong!!! To be revised!!! It provides plot in the center! Add ximg and yimg? Distortion also need to be considered!
### Jupiter position in x-y screen
juppos = - scpos_km
th, ph = sc.get_angles(juppos)
th = np.rad2deg(th)
ximg = th * np.cos(ph)
yimg = th * np.sin(ph)
# ax.plot([ximg], [yimg], '+w')
rjup = np.rad2deg(np.arcsin(p.rj / np.linalg.norm(juppos)))
print('Radius = ', rjup)
ax.plot(rjup * np.cos(np.linspace(0, 2 * np.pi, 45)),
rjup * np.sin(np.linspace(0, 2 * np.pi, 45)), 'w')
### Io orbit in x-y screen
ximg = []
yimg = []
for ang in np.linspace(0, np.pi * 2, 50):
iorelpos = 5.9 * rj * np.array([np.cos(ang), np.sin(ang), 0]) - scpos_km
th, ph = sc.get_angles(iorelpos)
th = np.rad2deg(th)
ximg.append(th * np.cos(ph))
yimg.append(th * np.sin(ph))
ax.plot(ximg, yimg, ':w')
### Magnetic field
if np.isfinite(sunsys3lon):
### System3 to Jeqs
dips3 = frames.magnetic_dipole_in_system3()
s3tojeqs = frames.system3_to_jeqs(sunsys3lon)
dipjeqs = s3tojeqs.dot(dips3)
print(dips3, dipjeqs)
dipedge = dipjeqs * 1.5 * p.rj - scpos_km
th, ph = sc.get_angles(dipedge)
th = np.rad2deg(th)
print(th, np.rad2deg(ph))
ximg = th * np.cos(ph)
yimg = th * np.sin(ph)
#### ax.plot([-ximg, ximg], [-yimg, yimg], '0.5')
### Fieldline to plot
ll = [10, ]
maglons = [0, 90, 180, 270]
magtojeqs = frames.jeqs_to_magnetic(sunsys3lon).T
for l in ll:
for maglon in maglons:
maglonr = np.deg2rad(maglon)
pointmag = l * np.array([np.cos(maglonr), np.sin(maglonr), 0]) * p.rj
pointjeqs = magtojeqs.dot(pointmag)
fl = dipole.fieldline(pointjeqs, ndiv=80, dipole_vector=dipjeqs) # [3,:]
ximg = []
yimg = []
for i in range(fl.shape[1]):
th, ph = sc.get_angles(fl[:, i] - scpos_km)
th = np.rad2deg(th)
ximg.append(th * np.cos(ph))
yimg.append(th * np.sin(ph))
r = (fl ** 2).sum(0)
print(r)
ximg = np.ma.masked_array(ximg, r < rj * rj)
yimg = np.ma.masked_array(yimg, r < rj * rj)
print(ximg)
behind = ((ximg * ximg + yimg * yimg < rjup * rjup) &
(fl.T.dot(scpos_km).T < 0))
print(behind)
ximg = np.ma.masked_where(behind, ximg)
yimg = np.ma.masked_where(behind, yimg)
ax.plot(ximg, yimg, '0.7', alpha=0.5)
# Same for Io field
pointjeqs = iopos * p.rj
fl = dipole.fieldline(pointjeqs, ndiv=80, dipole_vector=dipjeqs)
ximg = []
yimg = []
for i in range(fl.shape[1]):
th, ph = sc.get_angles(fl[:, i] - scpos_km)
th = np.rad2deg(th)
ximg.append(th * np.cos(ph))
yimg.append(th * np.sin(ph))
r = (fl ** 2).sum(0)
ximg = np.ma.masked_array(ximg, r < rj * rj)
yimg = np.ma.masked_array(yimg, r < rj * rj)
behind = ((ximg * ximg + yimg * yimg < rjup * rjup) &
(fl.T.dot(scpos_km).T < 0))
ximg = np.ma.masked_where(behind, ximg)
yimg = np.ma.masked_where(behind, yimg)
ax.plot(ximg, yimg, '0.7')
ax.set_xlim(xval.min(), xval.max())
ax.set_ylim(yval.min(), yval.max())
ax.set_aspect('equal')
ax.set_xlabel('Angle (Azimuth)')
ax.set_ylabel('Angle (Elevation)')
# ax.set_title('Pos: %.1f %.1f %.1f' % (scpos[0], scpos[1], scpos[2]))
ax.set_title('%s' % os.path.basename(filename)[:-4])
# ax.set_title('Oxygen ENA: 160 eV')
# fig.savefig('figure4a.eps')
fig2.savefig(outname)
if __name__ == "__main__":
# usage = "usage: %prog [options] arg"
# parser = OptionParser(usage)
#
# # Not used.
# parser.add_option("-f", "--file", dest="filename",
# help="read data from FILENAME")
# parser.add_option("-v", "--verbose",
# action="store_true", dest="verbose")
# parser.add_option("-q", "--quiet",
# action="store_false", dest="verbose")
# parser.add_option('-o', '--output',
# action='store', dest='output', default=None)
# parser.add_option('--use-classic-format',
# action='store_true', dest='use_classic_format', default=False)
#
# (options, args) = parser.parse_args()
# if len(args) != 1:
# parser.error("incorrect number of arguments")
#
# if options.use_classic_format:
# main(args[0], output=options.output)
#
# else:
t0 = dateutil.parser.parse('2031-05-19T02:00:00')
dt = datetime.timedelta(hours=4)
for i in range(14):
t = t0 + i * dt
infile = 'main_image_singlepoint_juice_orbit_O/2031-05-%02dT%02d:00:00.txt' % (t.day, t.hour)
outfile = 'figure5_%02d.png' % i
print(infile, outfile)
main_versionized(infile, outfile)