''' A plasma model.
Final goal is to get distribution function
at position (x,y,z) at time t (independent, so far).
'''
import logging
logging.basicConfig()
import numpy as np
import pyana.util.maxwell
import pyana.pep.iotorus
from . import frames
from . import jupiter
class CorotationalFlowMJconv:
''' Support coordinate system conversion.
The existing :class:`SimpleCorotationalFlow` or
:class:`SimpleCorotationalFlow2` support the plasma torus
model in JEqS frame implincitly.
This class support the frame conversion.
The interface of this class is for JEqS system.
The model used inside is for Magnetic frame.
'''
rj = 71492.
def __init__(self, jupiter_phase=None, m_amu=16, plasmatorusmodel=None):
if jupiter_phase == None:
self.jupiter_phase = jupiter.JupiterPhase()
else:
self.jupiter_phase = jupiter_phase
if plasmatorusmodel == None:
self.plasmatorus = pyana.pep.iotorus.PlasmaTorusEnvdocLSech6(10)
else:
self.plasmatorus = plasmatorusmodel
self.m = m_amu
@classmethod
def _km2rj(cls, v):
return np.array(v) / cls.rj
def get_density(self, x, y, z, tarr=None):
''' Get the density.
Given coordinate system is in JEqS.
Internally it is converted t
:keyword tarr: An array of sun direction (longitude) in system 3 frame, degrees.
'''
xrj, yrj, zrj = self._km2rj([x, y, z])
if np.isscalar(xrj):
xrj = np.array([xrj])
yrj = np.array([yrj])
zrj = np.array([zrj])
if tarr == None:
tarr = np.zeros_like(xrj)
if np.isscalar(tarr):
tarr = np.array([tarr])
sun_long_sys3 = self.jupiter_phase.get_subsolar_longitude_sys3(tarr)
xlist = []
ylist = []
zlist = []
for xx, yy, zz, ss in zip(xrj, yrj, zrj, sun_long_sys3):
jeqs2sys3 = frames.system3_to_jeqs(ss).T
sys32mag = frames.system3_to_magnetic()
jeqs2mag = sys32mag.dot(jeqs2sys3)
xm, ym, zm = jeqs2mag.dot(np.array([xx, yy, zz]))
xlist.append(xm)
ylist.append(ym)
zlist.append(zm)
density_list = self.plasmatorus.get_density(xlist, ylist, zlist)
return np.array(density_list) * 1e6 # in m^3
def get_distribution_function(self, x, y, z, vx, vy, vz, tarr=None):
xrj, yrj, zrj = self._km2rj([x, y, z])
if np.isscalar(xrj):
xrj = np.array([xrj])
yrj = np.array([yrj])
zrj = np.array([zrj])
if tarr == None:
tarr = np.zeros_like(xrj)
if np.isscalar(tarr):
tarr = np.array([tarr])
sun_long_sys3 = self.jupiter_phase.get_subsolar_longitude_sys3(tarr)
xlist = []
ylist = []
zlist = []
for xx, yy, zz, ss in zip(xrj, yrj, zrj, sun_long_sys3):
# jeqs2sys3 = frames.system3_to_jeqs(ss).T
# sys32mag = frames.system3_to_magnetic()
# jeqs2mag = sys32mag.dot(jeqs2sys3)
jeqs2mag = frames.jeqs_to_magnetic(ss)
xm, ym, zm = jeqs2mag.dot(np.array([xx, yy, zz]))
xlist.append(xm)
ylist.append(ym)
zlist.append(zm)
n = self.plasmatorus.get_density(xlist, ylist, zlist) * 1e6 # in m^3 (N,)
v = self.plasmatorus.get_velocity(xrj, yrj, zrj) * 1e3 # in m/s (3,N), it is in JEqS, OK?
kT = 70. # in eV. Kivelson 2004.
vth = np.sqrt(kT / self.m / 1.67e-27 * 1.60e-19) # in m/s
vx = np.array(vx) * 1e3 # in m/s
vy = np.array(vy) * 1e3 # in m/s
vz = np.array(vz) * 1e3 # in m/s
coeff = (2 * np.pi * vth ** 2) ** (-1.5)
vals = n * coeff * np.exp(-((v[0, :] - vx) ** 2 +
(v[1, :] - vy) ** 2 + (v[2, :] - vz) ** 2) / (2 * vth ** 2))
### Distribution function, MKSA.
return vals
class SimpleCorotationalFlow2:
''' A simple corotational flow model. Another high performance version.
Similar to :class:`SimpleCorotationalFlow`, but array support.
'''
rj = 71492.
logger = logging.getLogger('SimpleCorotationalFlow2')
#logger.setLevel(logging.DEBUG)
def __init__(self, m_amu=16., plasmatorusmodel=None):
if plasmatorusmodel == None:
self.plasmatorus = pyana.pep.iotorus.PlasmaTorusEnvdocLSech6(10)
else:
self.plasmatorus = plasmatorusmodel
self.m = m_amu
@classmethod
def _km2rj(cls, v):
return np.array(v) / cls.rj
def get_density(self, x, y, z):
''' Return the density in m^3 with inpput in km'''
xrj, yrj, zrj = self._km2rj([x, y, z])
return self.plasmatorus.get_density(xrj, yrj, zrj) * 1e6 # in m^3
def get_distribution_function(self, x, y, z, vx, vy, vz):
''' Return a distribution function at x, y, z
:param x: x coord in km.
:param vx: vx in km/s
:returns: VDF in s3/m6
>>> cf = SimpleCorotationalFlow()
>>> f = cf.get_distribution_function(420000., 0, 0)
'''
rj = self.rj
xrj = np.array(x) / rj
yrj = np.array(y) / rj
zrj = np.array(z) / rj
n = self.plasmatorus.get_density(xrj, yrj, zrj) * 1e6 # in m^3 (N,)
v = self.plasmatorus.get_velocity(xrj, yrj, zrj) * 1e3 # in m/s (3,N)
kT = 70. # in eV. Kivelson 2004.
vth = np.sqrt(kT / self.m / 1.67e-27 * 1.60e-19) # in m/s
vx = np.array(vx) * 1e3 # in m/s
vy = np.array(vy) * 1e3 # in m/s
vz = np.array(vz) * 1e3 # in m/s
coeff = (2 * np.pi * vth ** 2) ** (-1.5)
vals = n * coeff * np.exp(-((v[0, :] - vx) ** 2 +
(v[1, :] - vy) ** 2 + (v[2, :] - vz) ** 2) / (2 * vth ** 2))
### Distribution function, MKSA.
return vals
class SimpleCorotationalFlow:
''' A simple corotational flow model.
'''
rj = 71492.
logger = logging.getLogger('SimpleCorotationalFlow')
#logger.setLevel(logging.DEBUG)
def __init__(self, m_amu=16., rjlimit=10, plasmatorusmodel=None):
if plasmatorusmodel == None:
self.plasmatorus = pyana.pep.iotorus.PlasmaTorusEnvdocLSech6(10)
else:
self.plasmatorus = plasmatorusmodel
self.m = m_amu
self.set_rjlimit(rjlimit)
def set_rjlimit(self, rj):
self.rjlimit2 = rj * rj
@classmethod
def _km2rj(cls, v):
return np.array(v) / cls.rj
def get_density(self, x, y, z):
''' Return the density in m^3 with inpput in km'''
xrj, yrj, zrj = self._km2rj([x, y, z])
if xrj * xrj + yrj * yrj + zrj * zrj >= self.rjlimit2: # beyond 10 Rj
return 0
return self.plasmatorus.get_density(xrj, yrj, zrj) * 1e6 # in m^3
def get_distribution_function(self, x, y, z):
''' Return a distribution function at x, y, z
:param x: x coord in km.
>>> cf = SimpleCorotationalFlow()
>>> f = cf.get_distribution_function(420000., 0, 0)
'''
rj = self.rj
xrj = x / rj
yrj = y / rj
zrj = z / rj
if xrj * xrj + yrj * yrj + zrj * zrj >= self.rjlimit2: # beyond 10 Rj
return lambda x: 0
n = self.plasmatorus.get_density(xrj, yrj, zrj) * 1e6 # in m^3
v = self.plasmatorus.get_velocity(xrj, yrj, zrj) * 1e3 # in m/s
kT = 70. # in eV. Kivelson 2004.
vth = np.sqrt(kT / self.m / 1.67e-27 * 1.60e-19) # in m/s
self.logger.debug('N=%g' % n)
self.logger.debug('V=%g %g %g' % (v[0], v[1], v[2]))
self.logger.debug('Vth=%g' % vth)
### Distribution function, MKSA.
return pyana.util.maxwell.mkfunc(n, v, vth)
import unittest
import doctest
def doctests():
return unittest.TestSuite((
doctest.DocTestSuite(),
))
if __name__ == '__main__':
unittest.main(defaultTest='doctests')