Source code for irfpy.earth.volland

""" 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. """