"""
Define MagneticStatics class to calculate common static magnetic fields
as first raised in issue #100.
"""
__all__ = [
"CircularWire",
"FiniteStraightWire",
"GeneralWire",
"InfiniteStraightWire",
"MagneticDipole",
"MagnetoStatics",
"Wire",
]
import abc
import astropy.constants as const
import astropy.units as u
import numpy as np
import scipy.special
from plasmapy.utils.decorators import validate_quantities
[docs]
class MagnetoStatics(abc.ABC):
"""Abstract class for magnetostatic fields."""
[docs]
@abc.abstractmethod
def magnetic_field(self, p: u.Quantity[u.m]) -> u.Quantity[u.T]:
"""
Calculate magnetic field generated by this wire at position ``p``.
Parameters
----------
p : `~astropy.units.Quantity`
Three-dimensional position vector.
Returns
-------
B : `~astropy.units.Quantity`
Magnetic field at the specified position.
"""
[docs]
class MagneticDipole(MagnetoStatics):
r"""
Simple magnetic dipole — two nearby opposite point charges.
Parameters
----------
moment : `~astropy.units.Quantity`
Magnetic moment vector, in units of A m\ :sup:`2`\ .
p0 : `~astropy.units.Quantity`
Position of the dipole.
"""
@validate_quantities
def __init__(self, moment: u.Quantity[u.A * u.m**2], p0: u.Quantity[u.m]) -> None:
self.moment = moment.value
self._moment_u = moment.unit
self.p0 = p0.value
self._p0_u = p0.unit
def __repr__(self) -> str:
name = self.__class__.__name__
moment = self.moment
p0 = self.p0
moment_u = self._moment_u
p0_u = self._p0_u
return f"{name}(moment={moment}{moment_u}, p0={p0}{p0_u})"
[docs]
def magnetic_field(self, p: u.Quantity[u.m]) -> u.Quantity[u.T]:
r"""
Calculate magnetic field from this magnetic dipole at position ``p``.
Parameters
----------
p : `~astropy.units.Quantity`
Three-dimensional position vector.
Returns
-------
B : `~astropy.units.Quantity`
Magnetic field at the specified position.
Notes
-----
The magnetic field generated by a magnetic dipole is calculated
at a point in 3D space by taking the limit of a current loop as
its radius shrinks to a point and keeping its magnetic moment
constant.
Let the point where the magnetic field will be calculated be
represented by the point :math:`p` and the location of the
dipole by the point :math:`p_0` (with associated position vectors
:math:`\vec{p}` and :math:`\vec{p}_0`, respectively). Further,
let the displacement vector from the dipole at point
:math:`p_0` to the point :math:`p` be written as
:math:`\vec{r} = \vec{p} - \vec{p}_0`.
The magnetic field :math:`\vec{B}` from a magnetic dipole with
a magnetic dipole moment :math:`\vec{m}` is then found
at the point :math:`p` using the expression
.. math::
\vec{B} = \frac{\mu_0}{4\pi}
\left( \frac{3(\vec{m} \cdot \vec{r})\vec{r}}{|\vec{r}|^5}
- \frac{\vec{m}}{|\vec{r}|^3} \right)
where :math:`\mu_0` is vacuum permeability.
"""
r = p - self.p0
m = self.moment
B = (
const.mu0.value
/ 4
/ np.pi
* (
3 * r * np.dot(m, r) / np.linalg.norm(r) ** 5
- m / np.linalg.norm(r) ** 3
)
)
return B * u.T
[docs]
class Wire(MagnetoStatics):
"""Abstract wire class for concrete wires to be inherited from."""
[docs]
class GeneralWire(Wire):
r"""
General wire class described by its parametric vector equation.
Parameters
----------
parametric_eq : Callable
A vector-valued (with units of position) function of a single real
parameter.
t1 : `float`
Lower bound of the parameter, smaller than ``t2``.
t2 : `float`
Upper bound of the parameter, larger than ``t1``.
current: `~astropy.units.Quantity`
Electric current.
"""
@validate_quantities
def __init__(self, parametric_eq, t1, t2, current: u.Quantity[u.A]) -> None:
if callable(parametric_eq):
self.parametric_eq = parametric_eq
else:
raise TypeError("Argument parametric_eq should be a callable")
if t1 >= t2:
raise ValueError(f"t1={t1} is not smaller than t2={t2}")
self.t1 = t1
self.t2 = t2
self.current = current.value
self._current_u = current.unit
def __repr__(self) -> str:
name = self.__class__.__name__
parametric_eq = self.parametric_eq.__name__
t1 = self.t1
t2 = self.t2
current = self.current
current_u = self._current_u
return (
f"{name}(parametric_eq={parametric_eq}, t1={t1}, t2={t2}, "
f"current={current}{current_u})"
)
[docs]
def magnetic_field(self, p: u.Quantity[u.m], n: int = 1000) -> u.Quantity[u.T]:
r"""
Calculate magnetic field generated by this wire at position ``p``.
Parameters
----------
p : `~astropy.units.Quantity`
Three-dimensional position vector.
n : `int`, optional
Number of segments for Wire calculation (defaults to 1000).
Returns
-------
B : `astropy.units.Quantity`
Magnetic field at the specified position.
Notes
-----
The magnetic field generated by a wire with constant electric
current is found at a point in 3D space by approximating the
Biot–Savart law.
Let the point where the magnetic field will be calculated be
represented by the point :math:`p` (with associated position vector
:math:`\vec{p}`) and the curve :math:`C` defining the wire by the
parametric vector equation :math:`\vec{l}(t)` with
:math:`t_{\min} \leq t \leq t_{\max}`. Further, let the
displacement vector from the wire to the point :math:`p` be written
as :math:`\vec{r}(t) = \vec{p} - \vec{l}(t)`.
The magnetic field :math:`\vec{B}` generated by the wire with constant
current :math:`I` at point :math:`p` can then be expressed using
the Biot–Savart law, which takes the form
.. math::
\vec B = \frac{\mu_0 I}{4\pi}
\int_C \frac{d\vec{l} \times \vec{r}}{|\vec{r}|^3}
where :math:`\mu_0` is the permeability of free space.
This line integral is approximated by segmenting the wire into
:math:`n` straight pieces of equal length. The :math:`i\text{th}`
wire element can be written as
:math:`\Delta\vec{l}_i = \vec{l}(t_i) - \vec{l}(t_{i - 1})`
where :math:`t_i = t_{\min} + i(t_{\max}-t_{\min})/n`. Further,
the displacement vector from the center of the :math:`i\text{th}`
wire element to the position :math:`\vec{p}` is
:math:`\vec{r}_i = \vec{p} - \frac{\vec{l}(t_i) + \vec{l}(t_{i-1})}{2}`.
The integral is then approximated as
.. math::
\vec B \approx \frac{\mu_0 I}{4\pi}
\sum_{i=1}^{n}\frac{\vec{\Delta l}_i \times
\vec{r}_i}{\left| \vec{r}_i \right|^3}.
"""
p1 = self.parametric_eq(self.t1)
step = (self.t2 - self.t1) / n
t = self.t1
B = 0
for _ in range(n): # noqa: B007
t = t + step
p2 = self.parametric_eq(t)
dl = p2 - p1
p1 = p2
R = p - (p2 + p1) / 2
B += np.cross(dl, R) / np.linalg.norm(R) ** 3
B = B * const.mu0.value / 4 / np.pi * self.current
return B * u.T
[docs]
class FiniteStraightWire(Wire):
"""
Finite length straight wire class.
The ``p1`` to ``p2`` direction is the positive current direction.
Parameters
----------
p1 : `~astropy.units.Quantity`
Three-dimensional Cartesian coordinate of one end of the straight wire.
p2 : `~astropy.units.Quantity`
Three-dimensional Cartesian coordinate of another end of the straight wire.
current : `astropy.units.Quantity`
Electric current.
"""
@validate_quantities
def __init__(
self, p1: u.Quantity[u.m], p2: u.Quantity[u.m], current: u.Quantity[u.A]
) -> None:
self.p1 = p1.value
self.p2 = p2.value
self._p1_u = p1.unit
self._p2_u = p2.unit
if np.all(p1 == p2):
raise ValueError("p1, p2 should not be the same point.")
self.current = current.value
self._current_u = current.unit
def __repr__(self) -> str:
name = self.__class__.__name__
p1 = self.p1
p2 = self.p2
current = self.current
p1_u = self._p1_u
p2_u = self._p2_u
current_u = self._current_u
return f"{name}(p1={p1}{p1_u}, p2={p2}{p2_u}, current={current}{current_u})"
[docs]
def magnetic_field(self, p) -> u.Quantity[u.T]:
r"""
Calculate magnetic field generated by this wire at position ``p``.
Parameters
----------
p : `astropy.units.Quantity`
Three-dimensional position vector
Returns
-------
B : `astropy.units.Quantity`
Magnetic field at the specified position
Notes
-----
The magnetic field generated by a straight, finite wire with
constant electric current can be found at a point in 3D space
using the Biot–Savart law.
Let the point where the magnetic field will be calculated be
represented by the point :math:`p_0` (or ``p``) and the wire's
beginning and end as :math:`p_1` and :math:`p_2`, respectively
(with corresponding position vectors :math:`\vec{p}_0`,
:math:`\vec{p}_1`, and :math:`\vec{p}_2`, respectively).
Further, the vector from points
:math:`p_i` to :math:`p_j` can be written as
:math:`\vec{p}_{ij} = \vec{p}_j - \vec{p}_i`.
Next, consider the triangle with the points :math:`p_0`,
:math:`p_1`, and :math:`p_2` as vertices. The vector from the
vertex :math:`p_0` to the perpendicular foot opposite the
vertex :math:`p_0`, which will be used to find the unit vector
in the direction of the magnetic field, can be expressed as
.. math::
\vec{p}_f = \vec{p}_1 + \vec{p}_{12}
\frac{\vec{p}_{10} \cdot \vec{p}_{12}}
{|\vec{p}_{12}|^2}.
The magnetic field :math:`\vec{B}` generated by the wire with
current :math:`I` can be found at the point :math:`p_0` using
the Biot–Savart law which in this case simplifies to
.. math::
\vec{B} = \frac{\mu_0 I}{4π} (\cos θ_1 - \cos θ_2)
\hat{B}
where :math:`\mu_0` is the permeability of free space,
:math:`\theta_1` (:math:`\theta_2`) is the angle between
:math:`\vec{p}_{10}` (:math:`\vec{p}_{20}`) and
:math:`\vec{p}_{12}` with
.. math::
\cos\theta_1 = \frac{\vec{p}_{10} \cdot \vec{p}_{12}}
{|\vec{p}_{10}| |\vec{p}_{12}|}, \quad
\cos\theta_2 = \frac{\vec{p}_{20} \cdot \vec{p}_{12}}
{|\vec{p}_{20}| |\vec{p}_{12}|},
and
.. math::
\hat{B} = \frac{\vec{p}_{12} \times \vec{p}_{f0}}
{|\vec{p}_{12} \times \vec{p}_{f0}|}
is the unit vector in the direction of the magnetic field
at the point :math:`p_0`.
"""
# foot of perpendicular
p1, p2 = self.p1, self.p2
p2_p1 = p2 - p1
ratio = np.dot(p - p1, p2_p1) / np.dot(p2_p1, p2_p1)
pf = p1 + p2_p1 * ratio
# angles: theta_1 = <p - p1, p2 - p1>, theta_2 = <p - p2, p2 - p1>
cos_theta_1 = (
np.dot(p - p1, p2_p1) / np.linalg.norm(p - p1) / np.linalg.norm(p2_p1)
)
cos_theta_2 = (
np.dot(p - p2, p2_p1) / np.linalg.norm(p - p2) / np.linalg.norm(p2_p1)
)
B_unit = np.cross(p2_p1, p - pf)
B_unit = B_unit / np.linalg.norm(B_unit)
B = (
B_unit
/ np.linalg.norm(p - pf)
* (cos_theta_1 - cos_theta_2)
* const.mu0.value
/ 4
/ np.pi
* self.current
)
return B * u.T
[docs]
def to_GeneralWire(self):
"""Convert this `Wire` into a `GeneralWire`."""
p1, p2 = self.p1, self.p2
return GeneralWire(lambda t: p1 + (p2 - p1) * t, 0, 1, self.current * u.A)
[docs]
class InfiniteStraightWire(Wire):
"""
Infinite straight wire class.
Parameters
----------
direction:
Three-dimensional direction vector of the wire, also the
positive current direction.
p0 : `~astropy.units.Quantity`
One point on the wire.
current : `~astropy.units.Quantity`
Electric current.
"""
@validate_quantities
def __init__(
self, direction, p0: u.Quantity[u.m], current: u.Quantity[u.A]
) -> None:
self.direction = direction / np.linalg.norm(direction)
self.p0 = p0.value
self._p0_u = p0.unit
self.current = current.value
self._current_u = current.unit
def __repr__(self) -> str:
name = self.__class__.__name__
direction = self.direction
p0 = self.p0
current = self.current
p0_u = self._p0_u
current_u = self._current_u
return (
f"{name}(direction={direction}, p0={p0}{p0_u}, "
f"current={current}{current_u})"
)
[docs]
def magnetic_field(self, p) -> u.Quantity[u.T]:
r"""
Calculate magnetic field generated by this wire at position ``p``.
Parameters
----------
p : `astropy.units.Quantity`
Three-dimensional position vector.
Returns
-------
B : `astropy.units.Quantity`
Magnetic field at the specified position.
Notes
-----
The magnetic field generated by a straight wire with infinite
length and constant electric current is found at a point in 3D
space using the Biot–Savart law.
Let the point where the magnetic field will be calculated be
represented by the point :math:`p`, a point on the wire by
:math:`p_0`, and the direction of the wire as the vector
:math:`\vec{l}`. The magnetic field :math:`\vec{B}` generated by the
wire with constant current :math:`I` at point :math:`p` is
then expressed using the Biot–Savart law which takes the form
.. math::
\vec{B} = \frac{\mu_0 I}{2\pi |\vec{r}|} \hat{B}
where :math:`\mu_0` is the permeability of free space,
:math:`|\vec{r}| = |\vec{l} \times (\vec{p} - \vec{p}_0)|`
is the perpendicular distance between the wire and the point
:math:`p`, and
.. math::
\hat{B} = \frac{\vec{l} \times (\vec{p} - \vec{p}_0)}
{|\vec{l} \times (\vec{p} - \vec{p}_0)|}
is the unit vector in the direction of the magnetic field
at point :math:`p`.
"""
r = np.cross(self.direction, p - self.p0)
B_unit = r / np.linalg.norm(r)
r = np.linalg.norm(r)
return B_unit / r * const.mu0.value / 2 / np.pi * self.current * u.T
[docs]
class CircularWire(Wire):
"""
Circular wire (coil) class.
Parameters
----------
normal :
Three-dimensional normal vector of the circular coil.
center : `~astropy.units.Quantity`
Three-dimensional position vector of the circular coil's center.
radius: `~astropy.units.Quantity`
Radius of the circular coil.
current: `~astropy.units.Quantity`
Electric current.
"""
def __repr__(self) -> str:
name = self.__class__.__name__
normal = self.normal
center = self.center
radius = self.radius
current = self.current
center_u = self._center_u
radius_u = self._radius_u
current_u = self._current_u
return (
f"{name}(normal={normal}, center={center}{center_u}, "
f"radius={radius}{radius_u}, current={current}{current_u})"
)
@validate_quantities
def __init__(
self,
normal,
center: u.Quantity[u.m],
radius: u.Quantity[u.m],
current: u.Quantity[u.A],
n=300,
) -> None:
self.normal = normal / np.linalg.norm(normal)
self.center = center.value
self._center_u = center.unit
if radius > 0:
self.radius = radius.value
self._radius_u = radius.unit
else:
raise ValueError("Radius should be larger than 0")
self.current = current.value
self._current_u = current.unit
# parametric equation
# find other two axes in the disc plane
z = np.array([0, 0, 1])
axis_x = np.cross(z, self.normal)
axis_y = np.cross(self.normal, axis_x)
if np.linalg.norm(axis_x) == 0:
axis_x = np.array([1, 0, 0])
axis_y = np.array([0, 1, 0])
else:
axis_x = axis_x / np.linalg.norm(axis_x)
axis_y = axis_y / np.linalg.norm(axis_y)
self.axis_x = axis_x
self.axis_y = axis_y
def curve(t):
if not isinstance(t, np.ndarray):
return (
self.radius * (np.cos(t) * axis_x + np.sin(t) * axis_y)
+ self.center
)
t = np.expand_dims(t, 0)
axis_x_mat = np.expand_dims(axis_x, 1)
axis_y_mat = np.expand_dims(axis_y, 1)
return self.radius * (
np.matmul(axis_x_mat, np.cos(t)) + np.matmul(axis_y_mat, np.sin(t))
) + np.expand_dims(self.center, 1)
self.curve = curve
self.roots_legendre = scipy.special.roots_legendre(n)
self.n = n
[docs]
def magnetic_field(self, p) -> u.Quantity[u.T]:
r"""
Calculate magnetic field generated by this wire at position ``p``.
Parameters
----------
p : `~astropy.units.Quantity`
Three-dimensional position vector.
Returns
-------
B : `~astropy.units.Quantity`
Magnetic field at the specified position.
Notes
-----
The magnetic field generated by a circular wire with constant
electric current is found at a point in 3D space using the
Biot–Savart law. The integral in the Biot–Savart law is
approximated using the Gauss–Legendre quadrature.
Let the point where the magnetic field will be calculated be
represented by the point :math:`p` and the wire be represented by
the parametric vector equation
.. math::
\vec{l}(\theta) =
R\cos\theta \hat{x} + R\sin\theta \hat{y},\quad
-\pi \leq \theta \leq \pi
where :math:`R` is the radius of the circular wire.
Further, let the displacement vector from a point on the wire
to the point :math:`p` be written as
:math:`\vec{r}(\theta) = \vec{p} - \vec{l}(\theta)`.
The magnetic field :math:`B` due to a current :math:`I` is
then found at the point :math:`p` using the Biot–Savart law,
which takes the form
.. math::
\vec{B} = \frac{\mu_0 I}{4\pi}
\int_C \frac{d\vec{l} \times \vec{r}}{|\vec{r}|^3}.
This line integral is approximated using the Gauss–Legendre
quadrature with :math:`n` sample points:
.. math::
\hat{B} \approx \frac{\mu_0 I}{4\pi} \sum_{i=1}^n w_i
\frac{\Delta\vec{l}(\pi x_i) \times \vec{r}(\pi x_i)}
{|\vec{r}(\pi x_i)|^3}
where :math:`w_i` is the :math:`i\text{th}` quadrature weight and
:math:`x_i` is the :math:`i\text{th}` root of the :math:`n\text{th}`
Legendre polynomial.
"""
x, w = self.roots_legendre
t = x * np.pi
pt = self.curve(t)
dl = self.radius * (
-np.matmul(np.expand_dims(self.axis_x, 1), np.expand_dims(np.sin(t), 0))
+ np.matmul(np.expand_dims(self.axis_y, 1), np.expand_dims(np.cos(t), 0))
) # (3, n)
r = np.expand_dims(p, 1) - pt # (3, n)
r_norm_3 = np.linalg.norm(r, axis=0) ** 3
ft = np.cross(dl, r, axisa=0, axisb=0) / np.expand_dims(r_norm_3, 1) # (n, 3)
return (
np.pi
* np.matmul(np.expand_dims(w, 0), ft).squeeze(0)
* const.mu0.value
/ 4
/ np.pi
* self.current
* u.T
)
[docs]
def to_GeneralWire(self):
"""Convert this `Wire` into a `GeneralWire`."""
return GeneralWire(self.curve, -np.pi, np.pi, self.current * u.A)