r''' Dipole fields, such as electric or magnetic fields.
.. codeauthor:: Yoshifumi Futaana
See also :class:`irfpy.util.fields.DipoleField` to obtain the 3-D field
(vector at a specific position).
A dipole is expressed by the following equations.
(Look at any text book for derivation.)
The sample is for magnetic field.
.. math::
\vec{B} = \frac{\mu_0}{4\pi} \Bigl[\frac{3\vec{r}(\vec{r}\cdot\vec{m})-
r^2\vec{m}}{r^5} +
\frac{8\pi}{3}\vec{m}\delta^3(\vec{r})\Bigr]
If you change the first coefficient :math:`\frac{\mu_0}{4\pi}` to
:math:`\frac{1}{4\pi\epsilon_0}`, it expresses the electric dipole.
The above equation can be expressed in the polar coordinate system very simply.
.. math::
B_r = \frac{\mu_0 m}{4\pi}\frac{2\cos\theta}{r^3} \\
B_\theta = \frac{\mu_0 m}{4\pi}\frac{\sin\theta}{r^3}
For L-value calculation, use :func:`lvalue`.
'''
import numpy as np
[docs]def fieldline(crossing, dipole_vector=np.array([0, 0, 1]), ndiv=181):
r''' Calculate the field line in 3-D space.
Calculate the filed line that crosses the ``crossing`` point.
:param crossing: A numpy array shape=(3,) specifying the crossing point.
:keyword dipole_vector: A dipole vector (np.array with shape (3,)
placed at the origin.
:keyword ndiv: Resolution.
:returns: np.array in (3, ndiv) shaped coordinates of the field line.
See YF-Bnote 130304 for derivation.
'''
m = np.array(dipole_vector)
m = m / np.linalg.norm(m) # Dipole moment.
r = np.array(crossing) # ``r``, the reference point.
rr = np.linalg.norm(r) # Length of ``r``.
tr = np.arccos(np.inner(m, r) / rr) # Angle of ``r``.
lval = rr / (np.sin(tr) ** 2) # L-value.
fldlin = np.zeros([3, ndiv])
for i, tx in enumerate(np.linspace(0, np.pi, ndiv)[1:-1]): # ``x``, the field line element.
# ``tx``, the angle of ``x``.
ix = i + 1
rx = lval * (np.sin(tx) ** 2)
# ``rx``, the radius of ``x``.
alpha = rx * np.cos(tx) - rx * np.sin(tx) / np.tan(tr)
beta = rx * np.sin(tx) / (rr * np.sin(tr))
fldlin[:, ix] = alpha * m + beta * r
return fldlin
[docs]def fieldline_rt(l, ndiv=181):
r''' With a specific L-value, list of the coordinate system is returned.
:param l: The L-value. You may use :func:`lvalue` to derive the L-value
of the arbitrary position.
:keyword ndiv: Resolution. Returned array will have a shape of (``ndiv``)
:returns: Tuple of numpy array, (``r``, ``t``), where ``r`` is the radius
and ``t`` is the angle in radians.
For this function, a dipole moment is assumed to align z-axis.
The field line is (see :func:`lvalue`) expressed by
.. math::
r = L \sin^2\theta
>>> rlist, thlist = fieldline_rt(1.)
>>> print('%.2f %.2f' % (rlist[90], thlist[90])) # 90 deg.
1.00 1.57
>>> print('%.2f %.2f' % (rlist[45], thlist[45])) # 45 deg.
0.50 0.79
'''
anglist = np.linspace(0, np.pi, ndiv)
rlist = l * (np.sin(anglist) ** 2)
return rlist, anglist
[docs]def lvalue_rt(r, t):
r''' Return the L-value in the polar coordinates.
:param r: The (list of) distance.
:param t: The (list of) angle in degrees from the dipole vector direction.
>>> print(lvalue_rt([1, 2, 3], [45, 90, 135])) # doctest: +NORMALIZE_WHITESPACE
[2. 2. 6.]
'''
return r / (np.sin(np.deg2rad(t)) ** 2)
[docs]def lvalue(x, z):
r''' Return the L-value.
Return the L-value. Consider the magnetic dipole (or equivalent),
:math:`\vec{m}`, along z-axis.
The field line starting position of (x, z) hit the z=0 is
so-called L-value.
:param x: The (list of) position.
:param z: The (list of) position.
:returns l: The (list of) L-value.
>>> print(lvalue(5.0, 0.0))
5.0
>>> print(lvalue([5.0, 5.0], [0.0, 5.0])) # doctest: +NORMALIZE_WHITESPACE
[ 5. 14.14213562]
>>> print(lvalue(5.0, -5.0))
14.142135623730951
The field line can be traced by the equation
.. math::
\frac{B_r}{dr} = \frac{B_\theta}{rd\theta}
This can be rewritten as
.. math::
\frac{dr}{r} = 2\frac{\cos\theta}{\sin\theta}d\theta
Integration will give us
.. math::
\ln r = 2 \ln|\sin\theta| + C
or
.. math::
r = L \sin^2\theta
'''
xl2 = np.array(x) ** 2
zl2 = np.array(z) ** 2
rl2 = xl2 + zl2
return np.sqrt(rl2) * rl2 / xl2
def _distance_mother(theta, lvalue=1.0):
''' Mother function of distance calculator.
:param theta: in radian.
'''
t = np.cos(theta)
val = - np.array(lvalue) * ( 1 / 6. * (3 * np.sqrt(3 * (t ** 2) + 1) * t + np.sqrt(3)
* np.arcsinh(np.sqrt(3) * t)))
return val
[docs]def distance_along(theta0, theta1, lvalue):
r''' Calculate distance between the specified angles.
On the same l-value, the distance between the latitudes
:math:`\theta_0` and :math:`\theta_1` can be calculated,
and returned.
:param theta0: Angle in degrees.
:param theta1: Angle in degrees.
:param lvalue: The L-value.
For L=10, the distance between the pole and the equator is
>>> print('%.3f' % distance_along(0, 90, 10.0))
13.802
The distance between the poles are
>>> print('%.3f' % distance_along(0, 180, 10.0))
27.603
Numpy array can be given for angles.
>>> l0, l1, l2 = distance_along([0, 0, 90], [0, 90, 180], [5.0, 5.0, 20.0])
>>> print('%.3f' % l0)
0.000
>>> print('%.3f' % l1)
6.901
>>> print('%.3f' % l2)
27.603
*Derivation*
.. math::
x &=& r \sin\theta \\
z &=& r \cos\theta \\
l &=& \int d\theta\sqrt{\dot{x}^2+\dot{z}^2} \\
Here the dot specifies the differentiation by :math:`\theta`,
and therefore
.. math::
\dot{x} &=& L\sin^3\theta \\
\dot{z} &=& L\sin^2\theta\cos\theta \\
Substituting to :math:`l` and proceed the calculation, you will get
.. math::
l &=& L \int d\theta \sqrt{\sin^2\theta(1 + 3 \cos^2\theta)}
Using :math:`t=\cos\theta`,
.. math::
l &=& -L \int dt \sqrt{1+3t^2}
With a help of Mathematica :) you get
.. math::
l &=& -L \Bigl(\frac{1}{6}(3t\sqrt{3t^2+1}+\sqrt{3}\mathrm{arcsinh}(\sqrt{3}t))\Bigr)
Here, arcsinh is the inverse-hyperbolic-sin function.
This forumulation is valid for the angle between 0 and 180 degrees.
'''
l0 = _distance_mother(np.deg2rad(theta0), lvalue=lvalue)
l1 = _distance_mother(np.deg2rad(theta1), lvalue=lvalue)
return l1 - l0
import unittest
import doctest
[docs]def doctests():
return unittest.TestSuite((
doctest.DocTestSuite(),))
if __name__ == '__main__':
unittest.main(defaultTest='doctests')