✨ Inductance additions

bluemira/equilibria/coils/_tools.py

@@ -14,13 +14,19 @@
 
 import numpy as np
 
-from bluemira.magnetostatics.greens import circular_coil_inductance_elliptic, greens_psi
+from bluemira.magnetostatics.greens import (
+    circular_coil_inductance_elliptic,
+    greens_psi,
+    square_coil_inductance_kirchhoff,
+)
 
 if TYPE_CHECKING:
     from bluemira.equilibria.coils import CoilSet
 
 
-def make_mutual_inductance_matrix(coilset: CoilSet) -> np.ndarray:
+def make_mutual_inductance_matrix(
+    coilset: CoilSet, *, square_coil: bool = False, with_quadratures: bool = False
+) -> np.ndarray:
     """
     Calculate the mutual inductance matrix of a coilset.
 
@@ -37,28 +43,37 @@ def make_mutual_inductance_matrix(coilset: CoilSet) -> np.ndarray:
     -----
     Single-filament coil formulation; serves as a useful approximation.
     """
-    n_coils = coilset.n_coils()
+    if with_quadratures:
+        n_coils = coilset._quad_dx.ravel().size
+        xcoord = coilset._quad_x
+        zcoord = coilset._quad_z
+        dx = coilset._quad_dx
+        dz = coilset._quad_dz
+    else:
+        n_coils = coilset.n_coils()
+        xcoord = coilset.x
+        zcoord = coilset.z
+        dx = coilset.dx
+        dz = coilset.dz
+
     M = np.zeros((n_coils, n_coils))  # noqa: N806
-    xcoord = coilset.x
-    zcoord = coilset.z
-    dx = coilset.dx
-    dz = coilset.dz
     n_turns = coilset.n_turns
 
     itri, jtri = np.triu_indices(n_coils, k=1)
-
     M[itri, jtri] = (
         n_turns[itri]
         * n_turns[jtri]
-        * greens_psi(xcoord[itri], zcoord[itri], xcoord[jtri], zcoord[jtri])
+        * greens_psi(xcoord[itri], zcoord[itri], xcoord[jtri], zcoord[jtri]).ravel()
     )
     M[jtri, itri] = M[itri, jtri]
 
-    radius = np.hypot(dx, dz)
-    for i in range(n_coils):
-        M[i, i] = n_turns[i] ** 2 * circular_coil_inductance_elliptic(
-            xcoord[i], radius[i]
-        )
+    ind = np.diag_indices_from(M)
+
+    if square_coil:
+        M[ind] = square_coil_inductance_kirchhoff(xcoord, dx, dz)
+    else:
+        radius = np.hypot(dx, dz)
+        M[ind] = n_turns**2 * circular_coil_inductance_elliptic(xcoord, radius)
 
     return M
 

bluemira/magnetostatics/greens.py

@@ -215,7 +215,9 @@ def circular_coil_inductance_elliptic(
     return MU_0 * (2 * radius - rc) * ((1 - k**2 / 2) * ellipk_nb(k) - ellipe_nb(k))
 
 
-def circular_coil_inductance_kirchhoff(radius: float, rc: float) -> float:
+def circular_coil_inductance_kirchhoff(
+    radius: float | np.ndarray, rc: float | np.ndarray
+) -> float | np.ndarray:
     """
     Calculate the inductance of a circular coil by Kirchhoff's approximation.
 
@@ -239,6 +241,34 @@ def circular_coil_inductance_kirchhoff(radius: float, rc: float) -> float:
     return MU_0 * radius * (np.log(8 * radius / rc) - 2 + 0.25)
 
 
+def square_coil_inductance_kirchhoff(
+    radius: float | np.ndarray, width: float | np.ndarray, height: float | np.ndarray
+) -> float | np.ndarray:
+    """
+    Calculate the inductance of a square coil by Kirchhoff's approximation.
+
+    radius:
+        The radius of the square coil
+    width:
+        The width of the coil cross-section
+    height
+        The height of the coil cross-section
+
+    Returns
+    -------
+    The self-inductance of the square coil [H]
+
+    Notes
+    -----
+    .. math::
+
+        Inductance = \\mu_0 radius (ln(8\\frac{radius}{width + height}) - 0.5)
+
+    where :math:`\\mu_{0}` is the vacuum permeability
+    """
+    return MU_0 * radius * (np.log(8 * radius / (width + height)) - 0.5)
+
+
 @nb.jit(nopython=True)
 def greens_psi(
     xc: float | np.ndarray,