'''Models of ganymede magnetic field
A simple dipole is implemented in :class:`SimpleDipole`.
A tilted dipole is implemented in :class:`TiltedDipole`.
An MHD simulation provides more realisitic fields.
It is indeed implemented in :class:`irfpy.mhddata.MagneticField1201`.
'''
import numpy as np
[docs]class SimpleDipole:
r''' A simplest dipole approximation.
Centered at the center, no tilt.
It is obviously not valid, but for the very first approximation.
The formulation of the dipole field is
.. math::
\mathbf{B}(\mathbf{r}) = \frac{\mu_0}{4\pi}\frac{1}{r^3}
(3(\mathbf{m}\cdot\hat{\mathbf{r}})\hat{\mathbf{r}}
- \mathbf{m})
The usage is:
>>> sd = SimpleDipole()
>>> pos = np.array([2631.0, 0.0, 0.0]) # Position at equator
>>> bfield_equator = sd.get_field(pos)
>>> print('%.1f, %.1f, %.1f' % (bfield_equator[0], bfield_equator[1], bfield_equator[2]))
0.0, 0.0, 719.3
>>> pos = np.array([0.0, -2631.0, 0.0]) # Also at equator.
>>> bfield_equator2 = sd.get_field(pos)
>>> print('%.1f, %.1f, %.1f' % (bfield_equator2[0], bfield_equator2[1], bfield_equator2[2])) # The -0.0 will be a little bit ugly.
0.0, -0.0, 719.3
>>> pos = np.array([0.0, 0.0, 2631.0]) # At the pole.
>>> bfield_pole = sd.get_field(pos)
>>> print('%.1f, %.1f, %.1f' % (bfield_pole[0], bfield_pole[1], bfield_pole[2]))
-0.0, -0.0, -1438.6
The dipole_moment is set to the :attr:`dipole_moment`.
The default value is ``np.array([0, 0, -1.31e20])`` which produces ~719 nT
at the equatorial surface. This value is verified by Kivelson et al., 2002.
'''
mu0 = (4 * np.pi * 1e-7)
def __init__(self, dipole_moment=None):
self.dipole_moment = dipole_moment
''' The magnetic dipole moment. Preset to -1.31e20 only z direction
'''
if self.dipole_moment == None:
self.dipole_moment = np.array([0, 0, -1.31e20]) # in N m / T
[docs] def get_field(self, position_km):
''' Get the field at the position.
:param position_km: Position vector in km (x, y, z).
:type position_km: ``np.array`` with shape of (3,)
:returns: Mangetic field vector in nT.
:rtype: ``np.array`` with shape of (3,)
'''
pos_km = np.array(position_km)
if pos_km.shape != (3, ):
raise ValueError('Given position has diffent shape')
rvec = pos_km * 1000. # In meter
mvec = self.dipole_moment # For alias
r = np.sqrt((rvec ** 2).sum())
theta = np.arccos(rvec[2] / r) # Polar angle in radian
phi = np.arctan2(rvec[1], rvec[0]) # Azimuth angle in radian
rhat = rvec / r
field = mvec.dot(rhat) * 3. * rhat - mvec
coeff = self.mu0 / (4 * np.pi * (r ** 3))
return coeff * field * 1e9 # in nano-tesla
[docs] def interpolate3d(self, x, y, z):
''' Interface to get the magnetic field.
To make the same interface with :class:`irfpy.pep.mhddata.MagneticField1201`,
this method is prepared.
:param x: X position in Rg
:param y: Y position in Rg
:param z: Z position in Rg
:returns: A numpy array of the B-field with shape of (3,)
Since the dipole field is analytical and contiuous, no interpolation is
needed, but again ,this is prepared for the interface.
'''
pos = np.array([x, y, z]) * 2634.1
return self.get_field(pos)
[docs]class TiltedDipole(SimpleDipole):
r''' A second simplest model of the Ganymede field.
The second simplest model is dipole tilted a certain angle.
Indeed the implementation is quite the same to :class:`SimpleDipole`.
According to Kivelson et al., 2002, "it is tilted by 176 degrees
from the spin axis with the pole in the southern hemisphere
rotated by 24 degrees from the Jupiter-facing meridian plane
toward the trailing hemisphere."
The coordinate system plays into this game, indeed.
It is really annoying, but I will follow the SPICE for future consistency.
For SPICE, IAU_GANYMEDE coordinate system is defined.
This looks as "x points to jupiter; y opposite to ganymede motion".
The definition is indeed deferent from Figure 2 in Kivelson et al., 2002.
Nevertheless, the Kivelson's statement is quite complete.
Tilt of 176 deg is, as stated, in the southern hemisphere.
Rotate of 24 deg from Jupiter-facing meridian (=x-z plane in IAU)
toward the trailing hemisphere (=minus y in IAU).
This means that the magnetic dipole moment vector is
.. math::
m_x = m \sin(176.) * \cos(-24.) \\
m_y = m \sin(176.) * \sin(-24.) \\
m_z = m \cos(176.)
'''
def __init__(self):
m = 1.31e20
theta = 176. * np.pi / 180.
phi = -24. * np.pi / 180.
mvec = np.array([np.sin(theta) * np.cos(phi), np.sin(theta) * np.sin(phi), np.cos(theta)]) * m
SimpleDipole.__init__(self, dipole_moment=mvec)
import unittest
import doctest
[docs]def doctests():
return unittest.TestSuite((
doctest.DocTestSuite(),
))
if __name__ == '__main__':
unittest.main(defaultTest='doctests')