import os
import sys
import logging
logging.basicConfig()
import datetime
import math
import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
from optparse import OptionParser
from pyana.util import dipole
import pyana.pep.iotorus
from . import observer
from . import parameters as p
rj = p.rj
from . import plasmamodel
from . import jupiter
from . import frames
def plotscreen(xc, yc, val, ax=None, *args, **kwds):
if ax == None:
ax = plt.gca()
img = ax.pcolor(xc, yc, val[1:-1, 1:-1].T, **kwds) # , vmin=0, vmax=3)
ax.set_aspect(1)
return ax, img
def addiopos(ax, iopos_km, scpos_km, sc):
### Io position in x-y screen
iorelpos = iopos_km - scpos_km
print('iorelpos', iorelpos / rj)
th, ph = sc.get_angles(iorelpos)
th = np.rad2deg(th)
ximg = th * np.cos(ph)
yimg = th * np.sin(ph)
rio = np.rad2deg(np.arcsin(p.rio / np.linalg.norm(iorelpos))) * 10
# ax.plot([ximg], [yimg], marker='o', markerfacecolor='brown',
# markeredgecolor='none')
from matplotlib.patches import Circle
io = Circle((ximg, yimg), rio, color='#c89664')
ax.add_patch(io)
print('Io at projection', ximg, yimg)
# 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)) + ximg,
# r * np.sin(np.linspace(0, 2 * np.pi, 45)) + yimg, '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')
def addmagdipole(ax, sunsys3lon, scpos_km, iopos_km, sc):
rjup = np.rad2deg(np.arcsin(rj / np.linalg.norm(scpos_km)))
### 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')
### Same fo Io.
pointjeqs = iopos_km
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', lw=2)
def scangle2xyimg(th, ph):
''' th, ph in radian. xy in degrees.
'''
th2 = np.rad2deg(th)
return th2 * np.cos(ph), th2 * np.sin(ph)
def xyimg2scangle(ximg, yimg):
''' ximg, yimg in degrees, th, ph in radian.
Not ximg, yimg, but th, ph can go to sc.get_vector()
'''
th = np.deg2rad(np.sqrt(ximg **2 + yimg **2))
ph = np.arctan2(yimg, ximg)
return th, ph
def addjuppos(ax, scpos_km, sc):
### 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)
# print ximg, yimg
# ax.plot([ximg], [yimg], 'w+')
rjup = np.rad2deg(np.arcsin(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')
from matplotlib.patches import Circle
jup = Circle((0, 0), rjup, color='#c89664')
ax.add_patch(jup)
def main(output=None, iophase=96.7, sunsys3lon=-65.6):
scpos = np.array([11.36, -9.38, 4.10])
scpos_km = scpos * rj
print('scpos', scpos, scpos_km)
scvel = np.array([6.73, 10.89, -0.43])
print('scvel', scvel)
iophase = iophase
iopos = np.array([np.cos(np.deg2rad(iophase)), np.sin(np.deg2rad(iophase)), 0]) * 5.9
iopos_km = iopos * rj
print('iopos', iopos, iopos_km)
sunsys3 = sunsys3lon
# Neutral model
nm = pyana.pep.iotorus.TwoPeakModel()
nm = pyana.pep.iotorus.RotatingNeutralCloud(nm)
nm.set_iophase(np.deg2rad(iophase))
print('NM at Io', nm.get_density(iopos[0], iopos[1], iopos[2]))
# Plasma model
# pm = pyana.pep.iotorus.PlasmaTorusEnvdocLSech6(10)
pm = plasmamodel.CorotationalFlowMJconv(jupiter_phase=jupiter.JupiterPhase(sunsys3))
sc = observer.VirtualSpacecraft(scpos_km, scvel)
print(sc.primary_axis)
print(sc.secondary_axis)
print(sc.third_axis)
xcoords = np.linspace(-40, 40, 81)
ycoords = np.linspace(-30, 30, 61)
# xcoords = np.array([11, 11.5, 11.88, 12.5, 13])
# ycoords = np.array([2, 2.5, 3.27, 3.5, 4])
grid = np.zeros([len(xcoords), len(ycoords)])
xval, yval = np.meshgrid(xcoords, ycoords)
nval = np.zeros_like(grid)
pval = np.zeros_like(grid)
segment = np.linspace(0, 30, 200)
for ix, xx in enumerate(xcoords):
for iy, yy in enumerate(ycoords):
# nval[ix, iy] = np.abs(xx) + 1
# pval[ix, iy] = np.abs(yy) + 1
# Calculate the maximum of density along FoV.
th, ph = xyimg2scangle(xx, yy)
viewvec = sc.get_vector(th, ph)
losx = segment * viewvec[0]
losy = segment * viewvec[1]
losz = segment * viewvec[2]
nlos = nm.get_density(scpos[0] + losx, scpos[1] + losy, scpos[2] + losz)
nval[ix, iy] = nlos.max()
plos = pm.get_density(rj * (scpos[0] + losx), rj * (scpos[1] + losy), rj * (scpos[2] + losz)) / 1e6
pval[ix, iy] = plos.max()
# if ix == 2 and iy == 2:
# print 'Viewvec to io=', viewvec, th, ph
# for xxx, yyy, zzz, nnn in zip(scpos[0] + losx, scpos[1] + losy, scpos[2] + losz, nlos):
# print xxx, yyy, zzz, nnn
# print np.array(sc.get_angles(viewvec)) * 180 / np.pi
# 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.
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111)
ax, img = plotscreen(xc, yc, np.log10(nval), ax=ax, vmin=-1, vmax=3)
cb = fig.colorbar(img)
cb.set_label('Nn [/cc]')
addiopos(ax, iopos_km, scpos_km, sc)
addjuppos(ax, scpos_km, sc)
addmagdipole(ax, sunsys3, scpos_km, iopos_km, sc)
ax.set_xlim(-40, 40)
ax.set_ylim(-30, 30)
if output != None:
fig.savefig(output + '_n.png')
fig.savefig(output + '_n.eps')
fig2 = plt.figure(figsize=(8, 6))
ax2 = fig2.add_subplot(111)
ax2, img2 = plotscreen(xc, yc, np.log10(pval), ax=ax2, vmin=1, vmax=3.5)
cb2 = fig2.colorbar(img2)
cb2.set_label('Np [/cc]')
addiopos(ax2, iopos_km, scpos_km, sc)
addjuppos(ax2, scpos_km, sc)
addmagdipole(ax2, sunsys3, scpos_km, iopos_km, sc)
ax2.set_xlim(-40, 40)
ax2.set_ylim(-30, 30)
if output != None:
fig2.savefig(output + '_p.png')
fig2.savefig(output + '_p.eps')
def main_movie(output=None):
if output == None: output='plasmamodel_plot3d_movie'
if not os.path.exists(output):
try:
os.mkdir(output)
except:
pass
if not os.path.isdir(output):
raise ValueError('Path %s should be directory. Use -o option.' % output)
time_io = 40. # 40 hours
time_jm = 10. # 10 hours
w_io = 360 / time_io # deg/hr
w_jm = 360 / time_jm # deg/hr
dt = 0.5 # hr
iophase = 96.7
jmphase = -65.6
tlist = np.arange(0, 40, dt)
iophases = iophase + tlist * w_io
jmphases = jmphase - tlist * w_jm
for idx, (i, j) in enumerate(zip(iophases, jmphases)):
outputfile = os.path.join(output, 'plasmamodel_plot3d_%02d' % idx)
print(idx, i, j, outputfile)
main(output=outputfile, iophase=i, sunsys3lon=j)
plt.close('all')
if __name__ == "__main__":
usage = "usage: %prog [options]"
parser = OptionParser(usage)
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('--movie',
action='store_true', dest='movie', default=False)
(options, args) = parser.parse_args()
if options.movie:
main_movie(output=options.output)
else:
main(output=options.output)