'''CPEMv2 model implementation.
Reference is Futaana et al., 2017 (to be submitted).
CPEMv2 is composed of the following plasma parameters.
- Charge density (:attr:`plasma_density`)
- Ion velocity vector (:attr:`ion_velocity_corotation`, :attr:`ion_velocity_polar` and :attr:`ion_velocity_radial`)
- Ion temperature (:attr:`ion_temperature`)
- Hot (~keV) electron density (:attr:`hot_electron_density`)
- Hot (~keV) electron temperature (:attr:`hot_electron_temperature`)
There are four functions for each parameter:
- Reference model, only depending on the Jovi-center distance, :math:`R`.
- Mean model, depending on the Jovi-center distance and the magnetic local time, :math:`R, LT`.
- Percentile model, depending on the percentile (user can specify) together with above.
- Full model, depending on the magnetic latitude, but only for Charge density.
The new approach of CPEM against other older models are the introduction of percentile model.
For the percentile model, each parameter is expressed by a function of the position
in a magnetic coordinate system
(r, MLT) together with a percentile.
The unit for the interface are human-oriented, as follows:
- Distance, Rj = 71492 km.
- Angle, degrees = :math:`\pi` / 180 (=0.017453) radians.
- LT, hours, = 15 degrees.
- Percentile, fraction (0 to 1).
- Density, cm^-3, = 1e-6 m^-3
- Velocity, km/s, = 1e3 m/s.
The simplest usage of CPEM is
>>> from irfpy.cpem.cpemv2 import *
>>> n = plasma_density.percentile_model(10, 2.5, 0.5) # Density at 10 Rj, 2.5 hr, and p=0.5
>>> print('{:.2f}'.format(n))
12.69
'''
import abc
import numpy as np
import scipy.stats
def _hour2rad(hour):
""" A helper function to convert the time (Hour) to the angle (Radian)
"""
return (hour / 24.) * 2 * np.pi
class __ABCDMeta(abc.ABCMeta):
"""Metaclass that allows docstring 'inheritance' for abc.
Not at all useful for the code, or implementation. Only for docstring.
You do not need to worry about this class.
.. note::
(This is for maintainer)
In case you meet the error to import the following classes:
- Just remove this class, change the metaclass for CpemIEEE as ``abc.ABCMeta``.
- The docstring may not shown correctly, but the implementation should work as expected.
"""
def __new__(mcls, classname, bases, cls_dict):
cls = super().__new__(mcls, classname, bases, cls_dict)
for name, member in cls_dict.items():
if not getattr(member, '__doc__'):
member.__doc__ = getattr(bases[-1], name).__doc__
return cls
class _CpemIEEE(metaclass=__ABCDMeta):
""" This is a base class of CPEM/IEEE.
"""
@abc.abstractmethod
def reference_model(self, r):
""" Reference model implementation for CPEM/IEEE.
:param r: The Jovi-center radius in Rj.
"""
raise NotImplementedError('')
@abc.abstractmethod
def mean_model(self, r, mlt):
""" Mean model implementation for CPEM/IEEE.
:param r: The Jovi-center radius in Rj.
:param mlt: The magnetic local time in hours.
"""
raise NotImplementedError('')
@abc.abstractmethod
def percentile_model(self, r, mlt, p):
""" Percentile model implementation for CPEM/IEEE.
:param r: The Jovi-center radius in Rj.
:param mlt: The magnetic local time in hours.
:param p: The percentile. The number should be between 0 and 1. (Not percent...)
"""
raise NotImplementedError('')
def full_model(self, r, mlt, mlat, p):
""" Full model implementation for CPEM/IEEE.
:param r: The Jovi-center radius in Rj.
:param mlt: The magnetic local time in hours.
:param mlat: The magnetic latitude in degrees.
:param p: The percentile. The number should be between 0 and 1. (Not percent...)
"""
return self.percentile_model(r, mlt, p)
[docs]class CpemIonChargeDensity(_CpemIEEE):
""" CPEM/IEEE ion charge density models.
Let's see the models at R=20 [Rj] and MLT=10 [h].
>>> ne = CpemIonChargeDensity()
>>> print('{:.3f}'.format(ne.reference_model(20))) # Reference model at 20 Rj.
0.422
>>> print('{:.3f}'.format(ne.mean_model(20, 10))) # Mean model
0.238
>>> print('{:.3f}'.format(ne.percentile_model(20, 10, 0.5))) # P=50% (median) model
0.239
>>> print('{:.3f}'.format(ne.percentile_model(20, 10, 0.023))) # P=2.3% model (-2 sigma)
0.060
>>> print('{:.3f}'.format(ne.percentile_model(20, 10, 0.977))) # P=97.7% model (+2 sigma)
0.956
"""
def __init__(self):
_CpemIEEE.__init__(self)
self.parameter_a = [-0.6674, 2.160, 0.3509]
self.parameter_b = [0.02551, 1.929, -0.02353]
self.parameter_ppf = [0.3014, 0.002546]
[docs] def reference_model(self, r):
return _plasma_density_bagenal(r)
[docs] def mean_model(self, r, mlt):
ref = self.reference_model(r)
coeff_a = self.parameter_a[0] + np.sin(_hour2rad(mlt - self.parameter_a[1])) * self.parameter_a[2]
coeff_b = self.parameter_b[0] + np.sin(_hour2rad(mlt - self.parameter_b[1])) * self.parameter_b[2]
coeff = coeff_a + coeff_b * r
mean = (10 ** coeff) * ref
return mean
[docs] def percentile_model(self, r, mlt, p):
if not (0 <= p <= 1):
raise ValueError("Percentile should be given in the value between 0 and 1")
mean = self.mean_model(r, mlt)
ap = scipy.stats.norm.ppf(p, loc=self.parameter_ppf[1], scale=self.parameter_ppf[0])
perc = mean * np.power(10, ap)
return perc
[docs] def full_model(self, r, mlt, mlat, p):
""" Ion density full model.
"""
z = r * np.sin(np.radians(mlat))
Ai = 20
ti = CpemIonTemperature()
Ti = ti.percentile_model(r, mlt, 0.5)
H = 0.64 * np.sqrt(Ti / Ai)
LN = np.exp(-(z / H) ** 2)
return LN * self.percentile_model(r, mlt, p)
[docs]class CpemIonTemperature(_CpemIEEE):
""" CPEM/IEEE ion temperature models.
Let's see the models at R=15 [Rj] and MLT=20 [h].
>>> kTi = CpemIonTemperature()
>>> print('{:.1f}'.format(kTi.reference_model(15))) # Reference model at 20 Rj.
379.1
>>> print('{:.1f}'.format(kTi.mean_model(15, 20))) # Mean model
560.7
>>> print('{:.1f}'.format(kTi.percentile_model(15, 20, 0.5))) # P=50% (median) model
559.9
>>> print('{:.1f}'.format(kTi.percentile_model(15, 20, 0.023))) # P=2.3% model (-2 sigma)
166.0
>>> print('{:.0f}'.format(kTi.percentile_model(15, 20, 0.977))) # P=97.7% model (+2 sigma)
1888
"""
def __init__(self):
_CpemIEEE.__init__(self)
self.parameter_a = [0.2535, -1.849, 0.2711]
self.parameter_b = [-0.006874, -0.6528, -0.01425]
self.parameter_ppf = [0.2646, -0.0005726]
[docs] def reference_model(self, r):
return _ion_temperature_bagenal(r)
[docs] def mean_model(self, r, mlt):
ref = self.reference_model(r)
coeff_a = self.parameter_a[0] + np.sin(_hour2rad(mlt - self.parameter_a[1])) * self.parameter_a[2]
coeff_b = self.parameter_b[0] + np.sin(_hour2rad(mlt - self.parameter_b[1])) * self.parameter_b[2]
coeff = coeff_a + coeff_b * r
mean = (10 ** coeff) * ref
return mean
[docs] def percentile_model(self, r, mlt, p):
if not (0 <= p <= 1):
raise ValueError("Percentile should be given in the value between 0 and 1")
mean = self.mean_model(r, mlt)
ap = scipy.stats.norm.ppf(p, loc=self.parameter_ppf[1], scale=self.parameter_ppf[0])
perc = mean * np.power(10, ap)
return perc
[docs]class CpemIonVelCorotation(_CpemIEEE):
""" Corotational velocity model of CPEM.
Let's see the models at R=10 [Rj] and MLT=7.5 [h].
>>> vc = CpemIonVelCorotation()
>>> print('{:.1f}'.format(vc.reference_model(10))) # Reference model at 20 Rj.
116.8
>>> print('{:.1f}'.format(vc.mean_model(10, 7.5))) # Mean model
106.1
>>> print('{:.1f}'.format(vc.percentile_model(10, 7.5, 0.5))) # P=50% (median) model
106.8
>>> print('{:.2f}'.format(vc.percentile_model(10, 7.5, 0.023))) # P=2.3% model (-2 sigma)
77.66
>>> print('{:.1f}'.format(vc.percentile_model(10, 7.5, 0.977))) # P=97.7% model (+2 sigma)
147.0
The full model should return the results same sa the percentile model, disregarding the Mlat.
>>> print('{:.1f}'.format(vc.full_model(10, 7.5, -30, 0.977))) # -30 degrees mlat, but disregarded.
147.0
"""
def __init__(self):
_CpemIEEE.__init__(self)
self.parameter_a = [-0.1200, 1.415, -0.04894]
self.parameter_b = [0.007300, -0.7488, 0.006530]
self.parameter_ppf = [0.06941, 0.002892]
[docs] def reference_model(self, r):
return _ion_convection_speed_DG83(r)
[docs] def mean_model(self, r, mlt):
ref = self.reference_model(r)
coeff_a = self.parameter_a[0] + np.sin(_hour2rad(mlt - self.parameter_a[1])) * self.parameter_a[2]
coeff_b = self.parameter_b[0] + np.sin(_hour2rad(mlt - self.parameter_b[1])) * self.parameter_b[2]
coeff = coeff_a + coeff_b * r
mean = (10 ** coeff) * ref
return mean
[docs] def percentile_model(self, r, mlt, p):
if not (0 <= p <= 1):
raise ValueError("Percentile should be given in the value between 0 and 1")
mean = self.mean_model(r, mlt)
ap = scipy.stats.norm.ppf(p, loc=self.parameter_ppf[1], scale=self.parameter_ppf[0])
perc = mean * np.power(10, ap)
return perc
[docs]class CpemIonVelRadial(_CpemIEEE):
""" CPEM model for radial velocity.
>>> vr = CpemIonVelRadial()
>>> print('{:.2f}'.format(vr.reference_model(10))) # Reference model at 20 Rj.
0.00
>>> print('{:.2f}'.format(vr.mean_model(10, 7.5))) # Mean model
4.85
>>> print('{:.2f}'.format(vr.percentile_model(10, 7.5, 0.5))) # P=50% (median) model
3.14
>>> print('{:.2f}'.format(vr.percentile_model(10, 7.5, 0.023))) # P=2.3% model (-2 sigma)
-35.79
>>> print('{:.2f}'.format(vr.percentile_model(10, 7.5, 0.977))) # P=97.7% model (+2 sigma)
42.07
"""
def __init__(self):
_CpemIEEE.__init__(self)
self.parameter_a = [3.545e-3, 4.421, -9.322e-2]
self.parameter_b = [1.943e-3, 3.072, 9.787e-3]
self.parameter_ppf = [19.51, -1.708]
[docs] def reference_model(self, r):
if np.isscalar(r):
return 0
return np.zeros_like(r)
[docs] def mean_model(self, r, mlt):
vr_ref = self.reference_model(r)
vc_median = ion_velocity_corotation.percentile_model(r, mlt, 0.5) # Median model of corotation flow
coeff_a = self.parameter_a[0] + np.sin(_hour2rad(mlt - self.parameter_a[1])) * self.parameter_a[2]
coeff_b = self.parameter_b[0] + np.sin(_hour2rad(mlt - self.parameter_b[1])) * self.parameter_b[2]
coeff = coeff_a + coeff_b * r
vr_mean = coeff * vc_median
return vr_mean
[docs] def percentile_model(self, r, mlt, p):
vr_mean = self.mean_model(r, mlt)
ap = scipy.stats.norm.ppf(p, loc=self.parameter_ppf[1], scale=self.parameter_ppf[0])
vr_perc = vr_mean + ap
return vr_perc
[docs]class CpemIonVelPolar(_CpemIEEE):
""" CPEM model for polar velocity.
Let's see the models at R=10 [Rj] and MLT=7.5 [h].
>>> vp = CpemIonVelPolar()
>>> print('{:.2f}'.format(vp.reference_model(10))) # Reference model at 20 Rj.
0.00
>>> print('{:.2f}'.format(vp.mean_model(10, 7.5))) # Mean model
-17.67
>>> print('{:.2f}'.format(vp.percentile_model(10, 7.5, 0.5))) # P=50% (median) model
-15.29
>>> print('{:.2f}'.format(vp.percentile_model(10, 7.5, 0.023))) # P=2.3% model (-2 sigma)
-52.96
>>> print('{:.2f}'.format(vp.percentile_model(10, 7.5, 0.977))) # P=97.7% model (+2 sigma)
22.39
"""
def __init__(self):
_CpemIEEE.__init__(self)
self.parameter_a = [-0.2004, 4.466, -0.2248]
self.parameter_b = [0.009830, -1.535, 0.01385]
self.parameter_ppf = [18.88, 2.387]
[docs] def reference_model(self, r):
ref = np.zeros_like(r)
if np.isscalar(r):
ref = 0
return ref
[docs] def mean_model(self, r, mlt):
vr_ref = self.reference_model(r)
vc_median = ion_velocity_corotation.percentile_model(r, mlt, 0.5) # Median model of corotation flow
coeff_a = self.parameter_a[0] + np.sin(_hour2rad(mlt - self.parameter_a[1])) * self.parameter_a[2]
coeff_b = self.parameter_b[0] + np.sin(_hour2rad(mlt - self.parameter_b[1])) * self.parameter_b[2]
coeff = coeff_a + coeff_b * r
vr_mean = coeff * vc_median
return vr_mean
[docs] def percentile_model(self, r, mlt, p):
vr_mean = self.mean_model(r, mlt)
ap = scipy.stats.norm.ppf(p, loc=self.parameter_ppf[1], scale=self.parameter_ppf[0])
vr_perc = vr_mean + ap
return vr_perc
[docs]class CpemElectronDensity(_CpemIEEE):
"""CPEM/IEEE electron density models.
The reference model was produced by ourselves from GAL/PLS data analysis.
Let's see the models at R=20 [Rj] and MLT=10 [h].
>>> ne = CpemElectronDensity()
>>> print('{:.4f}'.format(ne.reference_model(20))) # Reference model at 20 Rj.
0.0235
>>> print('{:.4f}'.format(ne.mean_model(20, 10))) # Mean model
0.0180
>>> print('{:.4f}'.format(ne.percentile_model(20, 10, 0.5))) # P=50% (median) model
0.0180
>>> print('{:.5f}'.format(ne.percentile_model(20, 10, 0.023))) # P=2.3% model (-2 sigma)
0.00258
>>> print('{:.3f}'.format(ne.percentile_model(20, 10, 0.977))) # P=97.7% model (+2 sigma)
0.126
"""
def __init__(self):
_CpemIEEE.__init__(self)
self.parameter_a_ref = -0.6839
self.parameter_b_ref = -0.04727
self.parameter_a_mean = [0.009778, 3.785, -0.2286]
self.parameter_b_mean = [5.059e-4, -4.256e0, -8.443e-3]
self.parameter_ppf = [0.4228, -0.0004700]
[docs] def reference_model(self, r):
return np.power(10, self.parameter_a_ref + self.parameter_b_ref * r)
[docs] def mean_model(self, r, mlt):
ref = self.reference_model(r)
coeff_a = self.parameter_a_mean[0] + np.sin(_hour2rad(mlt - self.parameter_a_mean[1])) * self.parameter_a_mean[2]
coeff_b = self.parameter_b_mean[0] + np.sin(_hour2rad(mlt - self.parameter_b_mean[1])) * self.parameter_b_mean[2]
coeff = coeff_a + coeff_b * r
mean = (10 ** coeff) * ref
return mean
[docs] def percentile_model(self, r, mlt, p):
if not (0 <= p <= 1):
raise ValueError("Percentile should be given in the value between 0 and 1")
mean = self.mean_model(r, mlt)
ap = scipy.stats.norm.ppf(p, loc=self.parameter_ppf[1], scale=self.parameter_ppf[0])
perc = mean * np.power(10, ap)
return perc
[docs]class CpemElectronTemperature(_CpemIEEE):
"""CPEM/IEEE electron temperature models.
Reference model is produced by ourselves.
.. math::
\log_{10}{kT_e} = a + br + cr^2 + dr^3
Let's see the models at R=15 [Rj] and MLT=20 [h].
>>> kTe = CpemElectronTemperature()
>>> print('{:.0f}'.format(kTe.reference_model(15))) # Reference model at 15 Rj.
3385
>>> print('{:.0f}'.format(kTe.mean_model(15, 20))) # Mean model
2893
>>> print('{:.0f}'.format(kTe.percentile_model(15, 20, 0.5))) # P=50% (median) model
2804
>>> print('{:.0f}'.format(kTe.percentile_model(15, 20, 0.023))) # P=2.3% model (-2 sigma)
1038
>>> print('{:.0f}'.format(kTe.percentile_model(15, 20, 0.977))) # P=97.7% model (+2 sigma)
7571
"""
def __init__(self):
_CpemIEEE.__init__(self)
self.parameter_a_ref = 3.189
self.parameter_b_ref = 0.04540
self.parameter_c_ref = -0.001820
self.parameter_d_ref = 0.00002046
self.parameter_a_mean = [-0.08202, 0.5147, 0.1345]
self.parameter_b_mean = [0.002386, 0.9014, -0.007121]
self.parameter_ppf = [0.2162, -0.01354]
[docs] def reference_model(self, r):
return np.power(10, self.parameter_a_ref + r * (self.parameter_b_ref + r * (self.parameter_c_ref + r * self.parameter_d_ref)))
[docs] def mean_model(self, r, mlt):
ref = self.reference_model(r)
coeff_a = self.parameter_a_mean[0] + np.sin(_hour2rad(mlt - self.parameter_a_mean[1])) * self.parameter_a_mean[2]
coeff_b = self.parameter_b_mean[0] + np.sin(_hour2rad(mlt - self.parameter_b_mean[1])) * self.parameter_b_mean[2]
coeff = coeff_a + coeff_b * r
mean = (10 ** coeff) * ref
return mean
[docs] def percentile_model(self, r, mlt, p):
if not (0 <= p <= 1):
raise ValueError("Percentile should be given in the value between 0 and 1")
mean = self.mean_model(r, mlt)
ap = scipy.stats.norm.ppf(p, loc=self.parameter_ppf[1], scale=self.parameter_ppf[0])
perc = mean * np.power(10, ap)
return perc
ion_charge_density = plasma_density = CpemIonChargeDensity()
""" Ion density model.
"""
ion_temperature = CpemIonTemperature()
""" Ion temperature model.
"""
ion_velocity_corotation = CpemIonVelCorotation()
""" Ion corotation velocity model.
"""
ion_velocity_radial = CpemIonVelRadial()
""" Ion radial velocity model.
"""
ion_velocity_polar = CpemIonVelPolar()
""" Ion polar velocity model.
"""
hot_electron_density = CpemElectronDensity()
""" Electron density model.
"""
hot_electron_temperature = CpemElectronTemperature()
""" Electron temperature model.
"""
def _plasma_density_bagenal(r):
''' Plasma model by Bagenal and Delareme, 2011.
.. math::
n_\mathrm{Bagenal} = a_1 (R/6)^{b_1}+a_2 (R/6)^{b_2}+a_3 (R/6)^{b_3}
where a1 = 1987., b1 = -8.2, a2 = 14, b2 = -3.2, a3 = 0.05 and b3= -0.65.
:param r: Distance from Jupiter. In Rj. NOT L-shell value, although can be taken the same in practice.
:returns n: plasma density in /cm^3
.. note::
The distance is the actual distance of Galileo within the latitude range [-30, 30] degrees. In geographical coordinate system.
>>> print('%.3f' % _plasma_density_bagenal(10))
32.899
'''
a1 = 1987.
b1 = -8.2
a2 = 14
b2 = -3.2
a3 = 0.05
b3 = -0.65
r6 = r / 6.
n1 = a1 * r6 ** b1
n2 = a2 * r6 ** b2
n3 = a3 * r6 ** b3
return n1 + n2 + n3
def _ion_scaleheight_bagenal(R):
''' Ion scale height, in equation (6) in B&D, 2011.
.. math::
h = \log_{10}(H) = a_1 + a_2 r + a_3 r^2 + a_4 r^3 + a_5 r^4
and r=log_10(R).
a1 = -0.116, a2=2.14, a3=-2.05, a4=0.491 and a5=0.126. H is the distance from the centrifugal disk.
.. warning::
We wonder if R in the paper B&D 2011 is a typo of R/6.The division by 6 reproduces the plots almost perfectly.
Also, see equation (1) in B&D 2011 for density, where R/6 was used instead of just R.
:param R: Distance from Jupiter in Rj.
:returns H: Ion scale height in Rj
>>> print('%.3f' % _ion_scaleheight_bagenal(10))
1.835
'''
r = np.log10(R / 6.0) # Division by 6 creates almost perfect match to the plot. Also, the paper uses R/6 for density equations (1).
a1 = -0.116
a2 = 2.14
a3 = -2.05
a4 = 0.491
a5 = 0.126
h = a1 + r * a2 + r * r * a3 + r * r * r * a4 + r * r * r * r * a5
H = np.power(10, h)
return H
def _ion_temperature_bagenal(R):
''' Ion temperature model
From Bagenal and Delamere, 2011.
.. math::
T_i = A_i(H/0.64)^2
where Ai = 20 (ion mass). H is defined by
.. math::
h = \log_{10}(H) = a_1 + a_2 r + a_3 r^2 + a_4 r^3 + a_5 r^4
and r=log_10(R).
a1 = -0.116, a2=2.14, a3=-2.05, a4=0.491 and a5=0.126.
:param R: Distance from Jupiter in Rj.
:returns Ti: Ion temperature in eV
>>> print('%.3f' % _ion_temperature_bagenal(10))
164.334
'''
H = _ion_scaleheight_bagenal(R)
ai = 20.
Ti = ai * ((H / 0.64) ** 2)
return Ti
def _ion_convection_speed_DG83(rmeq):
'''Model to determine ion convection speed in the magnetic equatorial plane, from DG83
From Table 6 (DG83), the ion convections speed (azimuthal velocity) is:
.. math::
\mathrm{for}\quad r_meq < 7.9 R_j,
v = v_c
v = v_c , namely take the corotation velocity.
.. math::
v_c = r_meq*71492*\Omega_{\mathrm{jupiter}} \mathrm{km/s}
r_meq is in R_j, and should be converted to km by multiplying it with Jupiters radius (71492 km),
omega_jupiter is Jupiter's rotation rate, equal to: 2*pi/(9.925*60*60)
.. math::
\mathrm{for}\quad 7.9<= r_meq < 20 R_j,
v = 8.32 r_meq + 33.6 \mathrm{km/s}
.. math::
\mathrm{for}\quad 20 <= r_meq < 170 R_j,
v = 200 \mathrm{km/s}
.. warning::
Note: Jupiter's equatorial radius and rotation period (9.925 hours) are hard-coded!
Source of these constants: http://nssdc.gsfc.nasa.gov/planetary/factsheet/jupiterfact.html
:param r_meq: distance from planet center in Rj, in the magnetic equatorial plane
:returns v: ion convection speed in km/s
>>> print('%.3f' % _ion_convection_speed_DG83(10))
116.800
>>> print('%.3f' % _ion_convection_speed_DG83(20))
200.000
'''
r_meq = np.asarray(rmeq)
omega_jupiter = 2 * np.pi / (9.925 * 60 * 60)
r_jupiter = 71492
vc = r_meq * r_jupiter * omega_jupiter
v = np.asarray(vc)
flag_792 = np.logical_and(r_meq > 7.9, r_meq <= 20)
if flag_792.any():
v[flag_792] = 8.32 * r_meq[flag_792] + 33.6 # km/s
flag_20 = (r_meq > 20)
if flag_20.any():
v[flag_20] = 200 # km/s
if np.isscalar(rmeq):
v = float(v)
return(v)