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}
To obtaint the L-value from cartesian coordinate, use :func:`lvalue`.
To obtaint the L-value from polar coordinate, use :func:`lvalue_rt`.
'''
import numpy as np
import numpy as _np
from irfpy.util import unitq as _u
[docs]def shape_function(latitude_rad):
r""" Shape function, see :func:`shape_function_u`.
:param latitude in radians.
:returns: The value of the shape function. It varies between 1 to infinity.
"""
s = np.sin(latitude_rad)
c = np.cos(latitude_rad)
up = np.sqrt(1 + 3 * (np.sin(s) ** 2))
dn = np.power(c, 6)
return up / dn
[docs]def shape_function_u(latitude):
r""" Shape function ($\frac{\sqrt{1+3\sin^2\lambda}}{\cos^6\lambda}$).
For dipole field, the function
.. math::
S_D(\lambda) = \frac{\sqrt{1+3\sin^2\lambda}}{\cos^6\lambda}
appears always. Not sure the name (shape function) is correct, but I implemented it
in a form of usual function.
:param latitude: The latitude. Unit shall be added (unitq.radian or unitq.degree).
:returns: The value of the shape function. It varies between 1 to infinity.
>>> import irfpy.util.unitq as u
>>> print('{:.2f}'.format(shape_function_u(30 * u.deg)))
3.08
>>> print(shape_function_u(0 * u.deg))
1.0
You can give array.
>>> s_d = shape_function_u(np.array([0, 10, 20, 30]) * u.deg)
"""
return shape_function(latitude.nr(_u.rad))
[docs]def bounce_time_approx_u(l_distance, equator_pitchangle, energy, mass):
r""" The bounce time, approximated, with unit
The bounce time is the time for a particle start from the equator, go to one pole,
then another pole, back to the equator.
The exact form is in an integral form, but it can be approximated by (see Baumjohann text book)
.. math::
\tau_b=\frac{LR_E}{\sqrt{W/m}}(3.7-1.6\sin\alpha_{eq})
It depends on is the l-distance (L Re, i.e., non-normalized l-value),
pitchangle at the equator (alpha_eq), energy (W) and mass (m).
:param l_distance: The distance. Length unit
:param equator_pitchangle: The angle.
:param energy: Energy of the particle.
:param mass: Mass of the particle
:returns: The approximate bounce time.
>>> import irfpy.util.unitq as u
>>> l = 10000 * u.km
>>> d = 45 * u.deg
>>> e = 1000 * u.eV
>>> m = u.k.m_p # Proton mass
>>> tau_b = bounce_time_approx_u(l, d, e, m)
>>> print('{:.1f}'.format(tau_b.nr('s')))
83.0
"""
taub = l_distance / np.sqrt(energy / mass) * (3.7 - 1.6 * np.sin(equator_pitchangle.nr(_u.radian)))
return taub.s
[docs]def bounce_time_approx(l_distance_km, equator_pitchangle_deg, energy_eV, mass_kg):
r""" The bounce time, see :func:`bounce_time_approx`.
Parameters shall be in the specific unit.
:param l_distance_km: The distance, km
:param equator_pitchangle_deg: The angle, deg
:param energy_eV: Energy of the particle, eV
:param mass_kg: Mass of the particle, kg
:returns: The approximate bounce time in second.
>>> l = 10000
>>> d = 45
>>> e = 1000
>>> m = 1.67e-27
>>> tau_b = bounce_time_approx(l, d, e, m)
>>> print('{:.1f}'.format(tau_b))
82.9
"""
taub = bounce_time_approx_u(l_distance_km * _u.km, equator_pitchangle_deg * _u.deg, energy_eV * _u.eV, mass_kg * _u.kg).nr('s')
return taub
[docs]def drift_time_approx_u(energy, ldistance, equator_pitch_angle, charge, surface_equator_magnetic_field=31100*_u.nT, planetary_radius=6378*_u.km):
""" Approxiated drift time.
The drift time (time to rotate around earth due to magnetic forces - curvature and grad B)
is approximated as,
.. math::
<\tau_d>\sim\frac{\pi q B_E R_E^2}{3 L W}(0.35 + 0.15 \sin\alpha_{eq})^{-1}
by Baumjohann text book (3.20).
All should be unitful.
:param energy: Energy
:param l_distance: L-value (unit with ditance)
:param equator_pitch_angle: Equatorial pitchangle
:param charge: Charge
:keyword surface_equator_magnetic_field: Default for Earth 31100 nT, but for other bodies, you have to specify the surface magnetic field strength at the equator.
:keyword planetary_radius: The planetary radius unless the Earth
:returns: The approximate time to rotate around planets.
>>> import irfpy.util.unitq as u
>>> e = 1 * u.keV
>>> l = 20000 * u.km
>>> p = 0 * u.deg
>>> q = 1.60e-19 * u.C
>>> td = drift_time_approx_u(e, l, p, q)
>>> print('{:.1f}'.format(td.nr('hour')))
334.9
"""
t1 = np.pi * charge * surface_equator_magnetic_field * (planetary_radius ** 2)
t2 = 3 * ldistance / planetary_radius * energy
t3 = 0.35 + 0.15 * np.sin(equator_pitch_angle.nr(_u.rad))
td = t1 / t2 / t3
td = td.s
return td
[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''' Calculate the L-value of specific locations given in the polar coordinates.
:param r: The (list of) distance of the locations.
:param t: The (list of) angle in degrees of the location(s) measured 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''' Calculate the L-value of the specific location given in the Cartsian..
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