'''Io torus model
'''
import logging
logging.basicConfig()
import numpy as np
import scipy.interpolate
import irfpy.util.dipole
[docs]class SimpleNeutralCloudModel:
''' A simple neutral cloud model.
The coordinate system is Io fixed frame, meaning that the Io is
at (5.9, 0, 0) regardless of time. Jupiter is obviously the origin
of the frame. The Io rotation axis around the Jupiter is in z direction.
The unit of the length is normalized by the Jupiter radius.
The density is retrieved by the :meth:`get_density` method.
The returned is the in the unit of /cm^3.
The density model is assumed very simple.
* At Io, the maximum density of 500 /cc
* The center of the neutral cloud is r=5.9. Note that people knows
that this assumption is too crude.
* The radius (HWHM) of the neutral cloud is 0.1 Rj (~7000 km). Exponential distribution assumed. (Smyth+ 2005)
* In the azimuthal HWHM is 45 deg for the heading (plus) direction and 10 deg for trailing (minus) direction. Exponential distribution assumed.
>>> ncm = SimpleNeutralCloudModel()
>>> print(ncm.get_density(5.9, 0, 0))
500.0
'''
logger = None
def __init__(self, nmax=500., radius=0.1, head=45.0, trail=10.0):
self.logger = logging.getLogger(str(self.__class__))
self.logger.setLevel(logging.DEBUG)
self.logger.debug('Class instanced')
self.nmax=nmax
self.radius=radius
self.head=head
self.trail=trail
[docs] def get_density(self, x, y, z):
''' Return the density in the unit of cm-3
:param x: In Jupiter radii
>>> ncm = SimpleNeutralCloudModel()
>>> print(ncm.get_density(5.9, 0, 0))
500.0
>>> print(ncm.get_density(0, 5.9, 0))
125.0
>>> print(ncm.get_density(0, -5.9, 0))
7.8125
>>> print(ncm.get_density([0, 0], [5.9, -5.9], [0, 0]))
[ 125. 7.8125]
'''
nmax, radius, head, trail = self.nmax, self.radius, self.head, self.trail
head2 = np.deg2rad(head)
trail2 = np.deg2rad(trail)
# First, calculate the azimuthal angle.
azi = np.array(np.arctan2(y, x))
while (azi < 0).any():
azi = np.where(azi < 0, azi + 2 * np.pi, azi)
azi_head = azi / head2
azi_trail = (2 * np.pi - azi ) / trail2
n_azi_head = 2 ** (-azi_head ** 1)
n_azi_trail = 2 ** (-azi_trail ** 1)
n_azi = np.where(n_azi_head > n_azi_trail, n_azi_head, n_azi_trail)
n_azi = n_azi * nmax
# Then calculate the radius from the torus center
xt0 = np.cos(azi) * 5.9
yt0 = np.sin(azi) * 5.9
zt0 = np.zeros_like(xt0)
r2 = (xt0 - x) ** 2 + (yt0 - y) ** 2 + (zt0 - z) ** 2
r = np.sqrt(r2)
n = n_azi * (2 ** (-(r / radius) ** 1))
return n
[docs]class TwoPeakModel:
'''Two peak model
Two peak model of Io Neutral cloud.
Io is assumed to be at (5.9, 0, 0).
'''
def __init__(self):
self.first = SimpleNeutralCloudModel(nmax=15000, radius=0.05, head=10.0, trail=2.5)
self.second = SimpleNeutralCloudModel(nmax=50, radius=0.5, head=60.0, trail=10.)
[docs] def get_density(self, x, y, z):
return self.first.get_density(x, y, z) + self.second.get_density(x, y, z)
[docs]class RotatingNeutralCloud2:
''' Rotating neutral cloud model with array support.
Similar to :class:`RotatingNeutralCloud`, but support multiple
phases and positions.
For example, :class:`RotatingNeutralCloud` only support single
phase each time.
>>> incr = RotatingNeutralCloud(TwoPeakModel())
>>> incr.set_iophase(np.deg2rad(0.)) # Io phase 0.
>>> print('%.1f' % incr.get_density(0, 6.0, 0))
22.7
>>> print('%.1f' % incr.get_density(0, 6.0, 0))
22.7
>>> incr.set_iophase(np.deg2rad(30.)) # Io phase 30 deg.
>>> print('%.1f' % incr.get_density(0, 6.0, 0))
80.4
It may also accept the array of positions, if the
Io phase is fixed.
>>> incr.set_iophase(np.deg2rad(45.)) # Io phase 45 deg.
>>> xyz = incr.get_density([0., 0, 0, 0], [4, 5, 6, 7], [0, 0, 0, 0])
>>> print('%.1f %.1f %.1f %.1f' % (xyz[0], xyz[1], xyz[2], xyz[3]))
2.1 8.5 191.6 6.5
On the other hand, this :class:`RotatingNeutralCloud2` can give
phase for each data.
The usage is as follows. No ``set_iophase`` method allowed, while
you may use iophase keyword in calling :meth:`get_density`.
>>> incr2 = RotatingNeutralCloud2(TwoPeakModel())
>>> print('%.1f' % incr2.get_density(0, 6, 0, iophase=np.deg2rad(30)))
80.4
>>> n = incr2.get_density([0., 0, 0, 0, 0, 0], [6., 6, 4, 5, 6, 7],
... [0., 0, 0, 0, 0, 0], iophase=np.deg2rad([0., 30., 45., 45., 45., 45.]))
>>> print('%.1f %.1f %.1f %.1f %.1f %.1f' % (n[0], n[1], n[2], n[3], n[4], n[5]))
22.7 80.4 2.1 8.5 191.6 6.5
Of course you may still use :class:`RotatingNeutralCloud`, but this
:class:`RotatingNeutralCloud2` may provide better performance.
From a simple test in :mod:`apps130106_iotorus_performance.ioncrot_getdensity.py`,
about 10 times better performance.
'''
def __init__(self, model):
self.model = model
[docs] def get_density(self, x, y, z, iophase=None):
if np.isscalar(x): # If x, y, z is scalar, delegate to the RotatingNeutralCloud.
rc = RotatingNeutralCloud(self.model)
if iophase == None:
iophase = 0
rc.set_iophase(iophase)
return rc.get_density(x, y, z)
# pos = np.array([x, y, z]) # (3, N)
npos = len(x)
pos_converted = np.empty([3, npos])
pos_converted[2, :] = z
# Decide not to use matrix shape as for speed.
cosp, sinp = np.cos(iophase), np.sin(iophase)
xi, yi = cosp * x + sinp * y, -sinp * x + cosp * y
pos_converted[0, :] = xi
pos_converted[1, :] = yi
xarr, yarr, zarr = pos_converted
n = self.model.get_density(xarr, yarr, zarr)
return n
[docs]class RotatingNeutralCloud:
''' Rotating neutral cloud model
You can make object of the "non-rotating" model (:class:`TwoPeakModel`
and :class:`SimpleNeutralCloudModel`)
first, and wrap by this class.
>>> inc = TwoPeakModel()
>>> print('%.1f' % inc.get_density(0, 6.0, 0))
22.7
This is the density at 90 degrees ahead of the Io.
>>> incr = RotatingNeutralCloud(inc)
>>> incr.set_iophase(np.deg2rad(30))
>>> print('%.1f' % incr.get_density(-3.0, 5.196, 0))
22.7
To get the same value for Io phase 30 degrees,
the position is 120 degrees from the x axis, i.e. (-3.0, 5.196, 0).
>>> print('%.5f' %inc.get_density(5, 0, 0))
14.41595
>>> incr.set_iophase(np.deg2rad(90))
>>> print('%.5f' %incr.get_density(0, 5, 0))
14.41595
'''
def __init__(self, model):
'''
:param model: Neutral cloud model. Should have
``get_density`` method.
'''
self.phase = 0.
self.model = model
[docs] def set_iophase(self, phase_r):
''' Set the Io phase
:param phase_r: Phase in radians.
'''
self.phase = phase_r
self.matrix = np.array([[np.cos(phase_r), np.sin(phase_r), 0],
[-np.sin(phase_r), np.cos(phase_r), 0],
[0, 0, 1]])
[docs] def get_density(self, x, y, z):
pos = np.array([x, y, z])
xi, yi, zi = self.matrix.dot(pos)
return self.model.get_density(xi, yi, zi)
[docs]class PlasmaTorusEnvdoc:
'''Plasma torus model employed from environment document with extrapolation.
'''
def __init__(self):
r_rj = np.array([0, 1, 2, 4, 5., 6, 7, 9, 11, 13, 15, 20, 25])
n_cp_rj = np.log10(np.array([0.1, 0.4, 1.6, 25.6, 1e2, 1.7e3, 7.2e2, 4.3e1, 9.9, 3.9, 1.7, .23, .13])) # in cm-3
self.n_cp = scipy.interpolate.interp1d(r_rj, n_cp_rj)
[docs] def get_density(self, x, y, z):
''' Return the density
:param x: Scalar or (N,) shaped np.array. Unit of Rj.
:param y: Scalar or (N,) shaped np.array. Unit of Rj.
:param z: Scalar or (N,) shaped np.array. Unit of Rj.
:returns: Density in cm^3. Scalar or (N,) shaped np.array
>>> p = PlasmaTorusEnvdoc()
>>> print(p.get_density(3, 4, 0)) # r=5.
100.0
>>> print(p.get_density([2, 0, 0, 13], [0, 0, 6, 0], [0, 5, 0, 0]))
[ 1.60000000e+00 1.00000000e+02 1.70000000e+03 3.90000000e+00]
>>> print(p.get_density(3, 5, 100)) # r>25.
0.0
>>> print(p.get_density([3, 0], [5, 0], [100, 0])) # r>25.
[ 0. 0.1]
'''
x = np.array(x)
y = np.array(y)
z = np.array(z)
r = np.sqrt(x * x + y * y + z * z)
# Here, extrapolation to be 0.
r = np.ma.masked_where(r > 25, r)
rval = np.array(r.tolist(fill_value=np.nan))
n = 10 ** self.n_cp(rval)
n = np.ma.masked_invalid(n)
return np.array(n.tolist(fill_value=0.))
[docs] def get_velocity(self, x, y, z):
''' Return the velocity, in km/s
Corotational flow assumed, but no tilt of magnetic field assumed.
So in this method, the magnetic field axis is parallel to the
rotation axis.
>>> p = PlasmaTorusEnvdoc()
>>> print('%.2f' % p.get_velocity(5.9, 0, 0)[1])
74.17
>>> print('%.2f' % p.get_velocity(0, 0, 5.9)[1])
0.00
>>> print('%.2f' % p.get_velocity(0, 5.9 / np.sqrt(2), 5.9 / np.sqrt(2))[0])
-52.45
You can get multiple data at once.
>>> vels = p.get_velocity([0, 1, 2, 3], [0, 1, 2, 3], [1, 1, 2, 3])
>>> print(vels.shape)
(3, 4)
'''
x = np.array(x)
y = np.array(y)
z = np.array(z)
r = np.sqrt(x * x + y * y + z * z)
x = x / r # Normalize
y = y / r
z = z / r
tj = 9 * 3600 + 55.5 * 60 # Jupiter SYSTEM-III
rj = 71492.
v = 2 * np.pi * r / tj * rj
vx = -y * v
vy = x * v
vz = np.zeros_like(vx)
return np.array([vx, vy, vz])
[docs]class PlasmaTorusEnvdocL(PlasmaTorusEnvdoc):
''' Plasma torus model from environment document with extrapolation.
Constant density along same L-value.
Use :class:PlasmaTorusEnvdocLSech6` for more realistic analysis.
'''
def __init__(self):
PlasmaTorusEnvdoc.__init__(self)
[docs] def get_density(self, x, y, z):
x = np.array(x)
y = np.array(y)
z = np.array(z)
xy = np.sqrt(x * x + y * y)
l = irfpy.util.dipole.lvalue(xy, z)
l = np.ma.masked_where(l > 25, l)
lval = np.array(l.tolist(fill_value=np.nan))
n = 10 ** self.n_cp(lval)
n = np.ma.masked_invalid(n)
return np.array(n.tolist(fill_value=0.))
[docs]class PlasmaTorusEnvdocLScaleHeight(PlasmaTorusEnvdocL):
r''' Plasma torus model extended toward latitude.
Starting from Envdoc.
The "equatorial" plasma density is from the environment document.
See also (:class:`PlasmaTorusEnvdocL`).
``exp (-(z/H)^2)`` dependence is assumed, where z is the distance to the
"equator" from the observation point along B field.
H is the scale height, which can be retrieved assuming
.. math::
H = 0.64 \times \sqrt{Ti/Ai}
where Ti (=100 eV as default) is the temperature and Ai is the average mass (20, as deafult).
Under the dipole assumption, z is taken from the
:func:`util:irfpy.util.dipole.distance_along` function.
*More for developer*
See same-level subclass :class:`PlasmaTorusEnvdocLSech6` how to implement.
'''
def __init__(self, Ti=100.0, Ap=20.0):
PlasmaTorusEnvdocL.__init__(self)
self.scaleheight = 0.64 * np.sqrt(Ti / Ap)
[docs] def get_density(self, xx, yy, zz):
x = np.array(xx)
y = np.array(yy)
z = np.array(zz)
n = PlasmaTorusEnvdocL.get_density(self, x, y, z) # At equator, i.e. z=0
longitude = np.arctan2(y, x)
xy = np.sqrt(x ** 2 + y ** 2)
r = np.sqrt(x ** 2 + y ** 2 + z ** 2)
colatitude = np.rad2deg(np.arccos(z / r))
lvalue = irfpy.util.dipole.lvalue(xy, z)
distance = irfpy.util.dipole.distance_along(colatitude, 90, lvalue)
exppart = np.exp(-((distance / self.scaleheight) ** 2))
return n * exppart
[docs]class PlasmaTorusEnvdocLSech6(PlasmaTorusEnvdocL):
''' Plasma torus model extended from :class:`PlasmaTorusEnvdoc`.
The equatorial plasma density is from the environment document
(:class:`PlasmaTorusEnvdoc`).
In addition, ``sech(6t)`` dependence is added along the
L-value to simulate the
latidudal dependence. Note that ``t`` is tha latitude.
The "equator" here is slightly problematic.
The maximum density plane should be the best term.
As the coordinate system used here is indeed arbitrary
(or precisely spealing maximum density plane is THE plane
for z=0).
The maximum density plane is, indeed, between the Jupiter equatorial
plane (system III) and the magnetic equator.
'''
def __init__(self, idx=6.):
PlasmaTorusEnvdoc.__init__(self)
self.idx = idx
[docs] def get_density(self, x, y, z):
''' Return the density.
:param x: Scalar or (N,) shaped np.array. Unit of Rj.
:param y: Scalar or (N,) shaped np.array. Unit of Rj.
:param z: Scalar or (N,) shaped np.array. Unit of Rj.
:returns: Density in cm^3. Scalar or (N,) shaped np.array
>>> p = PlasmaTorusEnvdocLSech6()
>>> print(p.get_density(3, 4, 0)) # r=5. at equator.
100.0
>>> print(p.get_density([2, 0, 0, 13], [0, 0, 6, 0], [0, 5, 0, 0]))
[ 1.60000000e+00 0.00000000e+00 1.70000000e+03 3.90000000e+00]
>>> print(p.get_density(3, 5, 100)) # r>25.
0.0
>>> print(p.get_density([3, 0], [5, 0], [100, 1])) # r>25.
[ 0. 0.]
'''
n = PlasmaTorusEnvdocL.get_density(self, x, y, z)
x = np.array(x)
y = np.array(y)
z = np.array(z)
r = x * x + y * y + z * z
xy = x * x + y * y
t = np.arccos(np.sqrt(xy / r))
latdep = 1 / np.cosh(t * self.idx)
return n * latdep
[docs]class PlasmaTorusEnvdocCos2lat():
''' Plasma torus model extended from :class:`PlasmaTorusEnvdoc`.
The equatorial plasma density is from the environment document
(:class:`PlasmaTorusEnvdoc`).
In addition, the ``cos^2(lat)`` dependence is added.
The dependence is totally conceptional, so that it does not
provide any real configuration.
Rather, use :class:`PlasmaTorusEnvdocSech6`.
'''
def __init__(self):
self.envmod = PlasmaTorusEnvdoc()
[docs] def get_density(self, x, y, z):
''' Return the density
:param x: Scalar or (N,) shaped np.array. Unit of Rj.
:param y: Scalar or (N,) shaped np.array. Unit of Rj.
:param z: Scalar or (N,) shaped np.array. Unit of Rj.
:returns: Density in cm^3. Scalar or (N,) shaped np.array
>>> tor0 = PlasmaTorusEnvdoc()
>>> tor1 = PlasmaTorusEnvdocCos2lat()
In the equator, these models returns the same.
>>> print('%.2f' % tor0.get_density(5.9, 0, 0))
1280.57
>>> print('%.2f' % tor1.get_density(5.9, 0, 0))
1280.57
In the pole, first model returns high density
while cosine model returns 0.
>>> print('%.2f' % tor0.get_density(0, 0, 5.9))
1280.57
>>> print('%.2f' % tor1.get_density(0, 0, 5.9))
0.00
'''
n = self.envmod.get_density(x, y, z)
r = x * x + y * y + z * z
xy = x * x + y * y
t = np.arccos(np.sqrt(xy / r))
cos2t = np.cos(t) ** 2
# cos2t = 1 / np.cosh(t * 6)
return n * cos2t
import unittest
import doctest
[docs]def doctests():
return unittest.TestSuite((
doctest.DocTestSuite(),
))
if __name__ == '__main__':
unittest.main(defaultTest='doctests')