""" Simple implementation of Volland-Stern potential
This is an empirical model of electric potential of the Earth's magnetosphere.
The Volland-Stern model is a parametric model. We prepared several presets in the
original paper. (``V78_Model1``, ``V78_Model2``).
A typical parameter for Mercury is also available (``Mercury_Model``).
To obtain the electric potential, one can use methods starting with "from".
For example, one can get the electric potential for L-value 10 and Local time 18:00 as follows.
The unit is kV.
>>> potential = V78_Model1.fromLLT(10, 18)
>>> print('{:.2f}'.format(potential))
-6.57
To produce a model using own parameter, one can start from :class:`ModelParameter` as an
argument of :class:`VollandPotential`, as seen below.
>>> V78_Model1 = VollandPotential(ModelParams(15, 10, 15, 10, -1, 4, 4, 1.5 * (1 + 1j), 14.8, 13.0, 14.0, 16.5, 270))
"""
from dataclasses import dataclass as _dataclass
import math as _math
import numpy as _np
[docs]
class ModelParams:
""" Model parameter for Volland-Stern model. See the paper Volland 78 for details.
"""
def __init__(self,
sigmaPbar: float,
sigmaHbar: float,
sigmaP: float,
sigmaH: float,
P1: complex,
P2: complex,
n1: complex,
n2: complex,
A1: complex,
theta1: float, # Degrees
theta0: float,
theta2: float,
tauMax: float,
):
self.sigmaPbar = sigmaPbar
self.sigmaHbar = sigmaHbar
self.sigmaP = sigmaP
self.sigmaH = sigmaH
self.P1 = P1
self.P2 = P2
self.n1 = n1
self.n2 = n2
self.A1 = A1
self.theta1 = theta1
self.theta0 = theta0
self.theta2 = theta2
self.tauMax = tauMax
# Calculation from the given parameter
self.tauMax_r = _np.radians(tauMax)
self.L1 = lshell_deg(1, self.theta1)
self.L0 = lshell_deg(1, self.theta0)
self.L2 = lshell_deg(1, self.theta2)
self.p1bar = self._calc_p1bar(self.sigmaPbar, self.sigmaHbar, self.sigmaP, self.sigmaH,
self.P1, _np.radians(self.theta1))
self.p2bar = self._calc_p2bar(self.sigmaPbar, self.sigmaHbar, self.sigmaP, self.sigmaH,
self.P2, _np.radians(self.theta2))
self.C1 = self.p1bar * self.A1 / 2
self.C2 = self._calc_C2(self.p1bar, self.p2bar, self.L0, self.L1, self.L2, self.n1, self.n2, self.C1)
self.B1 = self.A1
self.B2 = self.A2 = 2 * self.C2 / self.p2bar
def _calc_p1bar(self, sigmaPbar, sigmaHbar, sigmaP, sigmaH, P1, theta1) -> complex:
t1 = sigmaPbar * P1
t2 = 1j * (sigmaH - sigmaHbar) / _np.cos(theta1)
p1bar = (t1 + t2) / sigmaP
return p1bar
def _calc_p2bar(self, sigmaPbar, sigmaHbar, sigmaP, sigmaH, P2, theta2) -> complex:
t1 = sigmaPbar * P2
t2 = 1j * (sigmaH - sigmaHbar) / _np.cos(theta2)
p2bar = (t1 + t2) / sigmaP
return p2bar
def _calc_C2(self, p1bar, p2bar, L0, L1, L2, n1, n2, C1) -> complex:
t1 = 2 / p1bar
t2 = abs(L0 / L1 - 1)
t3 = n1
t4 = 2 + n1
t_up = t1 + t2 * t3 / t4
t1 = 2 / p2bar
t2 = abs(L0 / L2 - 1)
t3 = n2
t4 = 2 + n2
t_dn = t1 + t2 * t3 / t4
return C1 * t_up / t_dn
[docs]
def lt2deg(tau_LT: float):
""" Convert the local time (0 to 24) to the angle (degrees)
:param tau_LT: Local time in hour
:return: Angle in degrees.
>>> print(lt2deg(18))
270.0
"""
return tau_LT / 24.0 * 360
[docs]
def lt2rad(tau_LT: float):
""" Convert the local time (0 to 24) to radians.
:param tau_LT: Local time in hour.
:return: Angle in radians
>>> print('{:.2f}'.format(lt2rad(12)))
3.14
"""
return tau_LT / 24.0 * (2 * _math.pi)
[docs]
def lshell(r: float, theta_rad: float):
r""" L-shell obtained from the colatitude in radians
:param r: Distance from the center, unit Re
:param theta_rad: Colatitude in radians
:return: The L-shell, $r/\sin\theta^2$
"""
return r / (_np.sin(theta_rad) ** 2)
[docs]
def lshell_deg(r: float, theta_deg: float):
r""" L-shell obtaine from the colatitude in degrees
:param r: Distance form the center. Rm.
:param theta_deg: Colatitude in degrees
:return: The L-shell, $r/\sin\theta^2$
>>> print('{:.2f}'.format(lshell_deg(2.5, 90)))
2.50
"""
return lshell(r, _np.radians(theta_deg))
[docs]
class VollandPotential:
""" V78 potential model implementation
"""
def __init__(self, model_parameters: ModelParams):
""" Make potential calculation class from the given parameter (Model)
:param model_parameters: Model parameters.
Model shall be given as in the paper by Volland-1978.
>>> model1 = VollandPotential(ModelParams(15, 10, 15, 10, -1, 4, 4, 1.5 * (1 + 1j), 14.8, 13.0, 14.0, 16.5, 270))
>>> model1.fromPosition(5, 15, 14) # Potential (in kV) at R=5, coLat=15deg, Localtime at 14:00
-3.807631535691775
>>> model1.fromLLT(17.5, 15.5) # L=17.5 and LT=15:30
-11.229792294255114
"""
self.model = model_parameters
[docs]
def fromLLT(self, lvalue: float, localtime: float) -> float:
""" Get the potential from the L-value and the local time
:params lvalue: The L-value
:params localtime: The local time in the unit of hours.
"""
tau = lt2rad(localtime)
if lvalue > self.model.L1:
phi_c = self._phi_pol(lvalue, tau)
elif lvalue > self.model.L0:
phi_c = self._phi_highopen(lvalue, tau)
elif lvalue > self.model.L2:
phi_c = self._phi_lowopen(lvalue, tau)
else:
phi_c = self._phi_midlat(lvalue, tau)
phi = abs(phi_c)
phi = -phi * _np.cos(tau - self.model.tauMax_r)
return float(phi)
[docs]
def fromLLTs(self, lvalues, localtimes):
""" Obtain the potentials from an array of L-value and LT.
:param lvalues: A numpy array of L-value
:param localtimes: A numpy array of local time in the unit of Hours.
"""
taus = lt2rad(localtimes)
phi_regA = self._phi_pol(lvalues, taus)
phi_regB = self._phi_highopen(lvalues, taus)
phi_regC = self._phi_lowopen(lvalues, taus)
phi_regD = self._phi_midlat(lvalues, taus)
# A bit weired but should work
phi_c = phi_regD
phi_c = _np.where(lvalues > self.model.L2, phi_regC, phi_c)
phi_c = _np.where(lvalues > self.model.L0, phi_regB, phi_c)
phi_c = _np.where(lvalues > self.model.L1, phi_regA, phi_c)
phi = _np.absolute(phi_c)
phi = -phi * _np.cos(taus - self.model.tauMax_r)
return phi
[docs]
def fromPositions(self, rs, coLats_deg, localtimes):
phi_c = self.fromLLTs(lshell_deg(rs, coLats_deg), localtimes)
return phi_c
[docs]
def fromPosition(self, r: float, coLat_deg: float, localtime: float) -> float:
""" Obtain the potential at the given position
:params r: The radial distance in the unit of Re
:params coLat_deg: The colatitude in degrees.
:params localtime: The local time in hours.
"""
phi_c = self.fromLLT(lshell_deg(r, coLat_deg), localtime)
return float(phi_c)
def _phi_pol(self, L, tau):
""" Calculating potential for the polar cap, colatitude below theta1
"""
return -1j * self.model.A1 * _np.power(L / self.model.L1, self.model.P1 / 2.) * _np.exp(1j * tau)
def _phi_midlat(self, L, tau):
""" Calculate the potential for the mid latitude, co-latitude higher than theta2
"""
return -1j * self.model.A2 * _np.power(L / self.model.L2, self.model.P2 / 2.) * _np.exp(1j * tau)
def _phi_highopen(self, L, tau):
""" Calculate the potential for the high-latitude cusp. Between theta1 and theta0
"""
m = self.model
t1 = -1j
t2 = m.B1
t3 = 1 - L / m.L1
t4 = 2 / (2 + m.n1)
t5 = 1 - L / m.L1
t6 = 1 - m.L0 / m.L1
t7 = _np.exp(1j * tau)
phi = t1 * (t2 + abs(t3) * (1 - t4 * _np.power(abs(t5 / t6), m.n1 / 2)) * m.C1) * t7
return phi
def _phi_lowopen(self, L, tau):
""" Calculate the potential for the low-latitude cusp. Between theta0 and theta2
"""
m = self.model
t1 = -1j
t2 = m.B2
t3 = 1 - L / m.L2
t4 = 2 / (2 + m.n2)
t5 = 1 - L / m.L2
t6 = 1 - m.L0 / m.L2
t7 = _np.exp(1j * tau)
phi = t1 * (t2 + abs(t3) * (1 - t4 * _np.power(abs(t5 / t6), m.n2 / 2)) * m.C2) * t7
return phi
V78_Model1 = VollandPotential(ModelParams(15, 10, 15, 10, -1, 4, 4, 1.5 * (1 + 1j), 14.8, 13.0, 14.0, 16.5, 270))
""" A parameter set for the Model1 in V78 paper
"""
V78_Model2 = VollandPotential(ModelParams(15, 10, 22, 24, -1, 4, 4, 1.5 * (1 + 1j), 44.1, 13.5, 15.5, 20.0, 19 / 24 * 360))
""" A parameter set for the Model2 in V78 paper
"""
Mercury_Model = VollandPotential(ModelParams(1, 1, 1, 1, -1, 4, 4, 1.5 * (1 + 1j), 10, 40, 42.5, 45, 270))
""" A parameter set for typical Mercury.
Experimental implementation. Please use carefully. Mercury is different from the Earth.
Note that the ionospheric parameters will be disappered regardless of the numbers (conductivity) by formulation.
"""