""" Boundary model.
Magnetopause:
- mp_xyz_aMSO_winslow13
Bow shock:
- bs_xyz_aMSO_winslow13
"""
from collections import namedtuple as _namedtuple
import numpy as _np
BoundaryTuple = _namedtuple('BoundaryTuple',
[
'theta',
'phi',
'x',
'y',
'z',
])
""" The returned type for the boundary model.
:param theta: Theta array (radians)
:param phi: Phi array (radians)
:param x: x
:param y: y
:param z: z
"""
[docs]class Magnetopause_Winslow13:
""" The magnetopause model by Winslow 2013
The model is described by two angles, |theta| and |phi|.
In the aMSO frame,
.. math::
r = R_{ss} \left(\frac{2}{1+\cos\theta}\right)^\alpha
.. math::
x = r\cos\theta
.. math::
y = r\cos\phi
.. math::
z = r\sin\phi + z_d
Not well-tested.
>>> mp = Magnetopause_MSM_Winslow13()
>>> t, p, x, y, z = mp.xyz()
>>> # x, y = mp.xy() # This method is not yet implemented
>>> y, z = mp.yz()
>>> x, z = mp.xz()
"""
def __init__(self, rss=1.45, alpha=0.5, zd=0.196):
""" Initialize with (average) parameters
:keyword rss: Rss in the unit of Rm.
:keyword alpha: alpha
:keyword zd: zd, the offset of the dipole (defaut, 0.196 Rm, for aMSO frame). If zd=0, aMSM frame.
See :func:`mp_xyz_aMSO_winslow13` for more details.
"""
self.rss = rss
self.alpha = alpha
self.zd = zd
[docs] def xyz(self, theta_array=None, phi_array=None):
""" Getting x, y, and z
:keyword theta_array: Array of theta to calculate, in radians. (0 to pi)
:keyword phi_array: Array of phi to calculate, in radians. (0 to 2pi)
"""
return mp_xyz_aMSO_winslow13(theta_array=theta_array, phi_array=phi_array, rss=self.rss, alpha=self.alpha, zd=self.zd)
[docs] def xy(self, theta_array=None):
# z = 0, not simple...
raise NotImplementedError('')
equator = xy
[docs] def yz(self, phi_array=None):
""" Return the x=0 slice, i.e., terminator plane
"""
# x = 0 means theta=90 deg.
if phi_array is None:
phi_array = _np.linspace(0, 2 * _np.pi, 361, endpoint=True)
else:
phi_array = _np.array(phi_array)
y = self.rss * _np.power(2, self.alpha) * _np.cos(phi_array)
z = self.rss * _np.power(2, self.alpha) * _np.sin(phi_array) + self.zd
return (y, z)
terminator = yz
[docs] def xz(self, theta_array=None):
""" Return the y=0 slice, i.e., noon midnight meridian plane.
y=0 slice.
"""
# y = 0
# phi = 90 deg or 270 deg
t, p, x, y, z = self.xyz(theta_array, _np.deg2rad([90, 270]))
x = x.flatten()
z = z.flatten()
a = _np.argsort(z)
x = x[a]
z = z[a]
return x, z
noon_midnight = xz
[docs]class Magnetopause_MSM_Winslow13(Magnetopause_Winslow13):
""" Magnetopause in the MSM frame.
Not well-tested.
>>> mp = Magnetopause_MSM_Winslow13()
>>> t, p, x, y, z = mp.xyz()
>>> x, y = mp.xy()
>>> y, z = mp.yz()
>>> x, z = mp.xz()
"""
def __init__(self, rss=1.45, alpha=0.5):
Magnetopause_Winslow13.__init__(self, rss=rss, alpha=alpha, zd=0)
[docs] def xy(self, theta_array=None):
return self.xz(theta_array=theta_array)
[docs]def mp_xyz_aMSO_winslow13(theta_array=None, phi_array=None, rss=1.45, alpha=0.5, zd=0.196):
r""" Winslow [2013] Magnetopause model
:keyword theta_array: Array of theta to calculate, in radians. (0 to pi)
:keyword phi_array: Array of phi to calculate, in radians. (0 to 2pi)
:keyword rss: Rss in the unit of Rm.
:keyword alpha: alpha
:keyword zd: zd, the offset of the dipole (0.196 Rm). If zd=0, aMSM frame.
:returns: Array of boundary. (theta, phi, x, y, z). In Rm.
:rtype: :class:`BoundaryTuple`
>>> t, p, x, y, z = mp_xyz_aMSO_winslow13()
"""
if theta_array is None:
theta_array = _np.linspace(0, _np.pi, 180, endpoint=False)
else:
theta_array = _np.array(theta_array)
if phi_array is None:
phi_array = _np.linspace(0, 2 * _np.pi, 361, endpoint=True)
else:
phi_array = _np.array(phi_array)
theta, phi = _np.meshgrid(theta_array, phi_array, indexing='ij')
r = rss * _np.power(2 / (1 + _np.cos(theta)), alpha)
x = r * _np.cos(theta)
rho = r * _np.sin(theta)
y = rho * _np.cos(phi)
z = rho * _np.sin(phi) + zd
return BoundaryTuple(theta, phi, x, y, z)
[docs]def bs_xyz_aMSO_winslow13(theta_array=None, phi_array=None, zd=0.196, X0=0.5, e=1.04, p=2.75):
""" Winslow [2013] bow shock model
:keyword theta_array: Array of theta to calculate, in radians. (0 to pi)
:keyword phi_array: Array of phi to calculate, in radians. (0 to 2pi)
:keyword zd: zd, the offset of the dipole (0.196 Rm)
:keyword X0: X0 in the unit of Rm.
:keyword e: e
:keyword p: p
:returns: Array of boundary. (theta, phi, x, y, z). In Rm.
:rtype: :class:`BoundaryTuple`
>>> t, p, x, y, z = bs_xyz_aMSO_winslow13()
"""
if theta_array is None:
theta_array = _np.linspace(0, _np.pi, 180, endpoint=False)
else:
theta_array = _np.array(theta_array)
if phi_array is None:
phi_array = _np.linspace(0, 2 * _np.pi, 361, endpoint=True)
else:
phi_array = _np.array(phi_array)
theta, phi = _np.meshgrid(theta_array, phi_array, indexing='ij')
r2 = p * e / (1 + e * _np.cos(theta))
x = r2 * _np.cos(theta) + X0
rho = r2 * _np.sin(theta)
y = rho * _np.cos(phi)
z = rho * _np.sin(phi) + zd
return BoundaryTuple(theta, phi, x, y, z)