irfpy.util.dipole

Dipole fields, such as electric or magnetic fields.

Code author: Yoshifumi Futaana

See also 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.

\[\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 \(\frac{\mu_0}{4\pi}\) to \(\frac{1}{4\pi\epsilon_0}\), it expresses the electric dipole.

The above equation can be expressed in the polar coordinate system very simply.

\[\begin{split}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}\end{split}\]

To obtaint the L-value from cartesian coordinate, use lvalue().

To obtaint the L-value from polar coordinate, use lvalue_rt().

irfpy.util.dipole.shape_function(latitude_rad)

Shape function, see shape_function_u().

:param latitude in radians. :returns: The value of the shape function. It varies between 1 to infinity.

irfpy.util.dipole.shape_function_u(latitude)

Shape function ($frac{sqrt{1+3sin^2lambda}}{cos^6lambda}$).

For dipole field, the function

\[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.

Parameters:

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)

irfpy.util.dipole.bounce_time_approx_u(l_distance, equator_pitchangle, energy, mass)

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)

\[\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).

Parameters:

• l_distance – The distance. Length unit

• equator_pitchangle – The angle.

• energy – Energy of the particle.

• 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

irfpy.util.dipole.bounce_time_approx(l_distance_km, equator_pitchangle_deg, energy_eV, mass_kg)

The bounce time, see bounce_time_approx().

Parameters shall be in the specific unit.

Parameters:

• l_distance_km – The distance, km

• equator_pitchangle_deg – The angle, deg

• energy_eV – Energy of the particle, eV

• 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

irfpy.util.dipole.drift_time_approx_u(energy, ldistance, equator_pitch_angle, charge, surface_equator_magnetic_field=array(31100.) * nT, planetary_radius=array(6378.) * km)

Approxiated drift time.

The drift time (time to rotate around earth due to magnetic forces - curvature and grad B) is approximated as,

\[<\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.

Parameters:

• energy – Energy

• l_distance – L-value (unit with ditance)

• equator_pitch_angle – Equatorial pitchangle

• charge – Charge

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

• 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

irfpy.util.dipole.fieldline(crossing, dipole_vector=array([0, 0, 1]), ndiv=181)

Calculate the field line in 3-D space.

Calculate the filed line that crosses the crossing point.

Parameters:

• crossing – A numpy array shape=(3,) specifying the crossing point.

• dipole_vector – A dipole vector (np.array with shape (3,) placed at the origin.

• ndiv – Resolution.

Returns:

np.array in (3, ndiv) shaped coordinates of the field line.

See YF-Bnote 130304 for derivation.

irfpy.util.dipole.fieldline_rt(l, ndiv=181)

With a specific L-value, list of the coordinate system is returned.

Parameters:

• l – The L-value. You may use lvalue() to derive the L-value of the arbitrary position.

• 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 lvalue()) expressed by

\[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

irfpy.util.dipole.lvalue_rt(r, t)

Calculate the L-value of specific locations given in the polar coordinates.

Parameters:

• r – The (list of) distance of the locations.

• 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]))
[2.  2.  6.]

irfpy.util.dipole.lvalue(x, z)

Calculate the L-value of the specific location given in the Cartsian..

Return the L-value. Consider the magnetic dipole (or equivalent), \(\vec{m}\), along z-axis. The field line starting position of (x, z) hit the z=0 is so-called L-value.

Parameters:

• x – The (list of) position.

• 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]))
[  5.          14.14213562]

>>> print(lvalue(5.0, -5.0))
14.142135623730951

The field line can be traced by the equation

\[\frac{B_r}{dr} = \frac{B_\theta}{rd\theta}\]

This can be rewritten as

\[\frac{dr}{r} = 2\frac{\cos\theta}{\sin\theta}d\theta\]

Integration will give us

\[\ln r = 2 \ln|\sin\theta| + C\]

or

\[r = L \sin^2\theta\]

irfpy.util.dipole.distance_along(theta0, theta1, lvalue)

Calculate distance between the specified angles.

On the same l-value, the distance between the latitudes \(\theta_0\) and \(\theta_1\) can be calculated, and returned.

Parameters:

• theta0 – Angle in degrees.

• theta1 – Angle in degrees.

• 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

\[\begin{split}x &=& r \sin\theta \\ z &=& r \cos\theta \\ l &=& \int d\theta\sqrt{\dot{x}^2+\dot{z}^2} \\\end{split}\]

Here the dot specifies the differentiation by \(\theta\), and therefore

\[\begin{split}\dot{x} &=& L\sin^3\theta \\ \dot{z} &=& L\sin^2\theta\cos\theta \\\end{split}\]

Substituting to \(l\) and proceed the calculation, you will get

\[l &=& L \int d\theta \sqrt{\sin^2\theta(1 + 3 \cos^2\theta)}\]

Using \(t=\cos\theta\),

\[l &=& -L \int dt \sqrt{1+3t^2}\]

With a help of Mathematica :) you get

\[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.