RZIp vertical stability criterion

bluemira/equilibria/coils/_coil.py

@@ -127,6 +127,7 @@
         "_flag_sizefix",
         "_j_max",
         "_number",
+        "_resistance",
         "_x",
         "_x_boundary",
         "_z",
@@ -148,6 +149,7 @@
         b_max: float = np.nan,
         discretisation: float = np.nan,
         n_turns: int = 1,
+        resistance: float = 0,
         *,
         psi_analytic: bool = False,
         Bx_analytic: bool = True,
@@ -172,6 +174,7 @@
         self.ctype = ctype
         self.name = name
         self.n_turns = n_turns
+        self.resistance = resistance
 
         self._number = CoilNumber.generate(self.ctype)
         if self.name is None:
@@ -196,6 +199,7 @@
             f"{type(self).__name__}({self.name} ctype={self.ctype.name} x={self.x:.2g}"
             f" z={self.z:.2g} dx={self.dx:.2g} dz={self.dz:.2g}"
             f" current={self.current:.2g} j_max={self.j_max:.2g} b_max={self.b_max:.2g}"
+            f" resistance={self.resistance:.2g}"
             f" discretisation={self.discretisation:.2g})"
         )
 
@@ -415,6 +419,16 @@
         """Set coil max field"""
         self._b_max = floatify(value)
 
+    @property
+    def resistance(self):
+        """Get coil resistance"""
+        return self._resistance
+
+    @resistance.setter
+    def resistance(self, value: float):
+        """Set coil resistance"""
+        self._resistance = value
+
     @discretisation.setter
     def discretisation(self, value: float):
         """Set coil discretisation"""
@@ -425,6 +439,7 @@
         self,
         j_max: float = NBTI_J_MAX,
         b_max: float = NBTI_B_MAX,
+        resistance: float = 0,
     ) -> None:
         """
         Assigns EM material properties to coil
@@ -443,6 +458,7 @@
         """
         self.j_max = j_max
         self.b_max = b_max
+        self.resistance = resistance
 
     def get_max_current(self):
         """

bluemira/equilibria/coils/_field.py

@@ -16,7 +16,12 @@
 
 from bluemira.base.constants import MU_0
 from bluemira.equilibria.constants import X_TOLERANCE
-from bluemira.magnetostatics.greens import greens_Bx, greens_Bz, greens_psi
+from bluemira.magnetostatics.greens import (
+    greens_Bx,
+    greens_Bz,
+    greens_dbz_dx,
+    greens_psi,
+)
 from bluemira.magnetostatics.semianalytic_2d import (
     semianalytic_Bx,
     semianalytic_Bz,
@@ -59,6 +64,50 @@
         self._Bx_analytic = Bx_analytic
         self._Bz_analytic = Bz_analytic
 
+    def dB_d(self, x: float | np.ndarray, z: float | np.ndarray):
+        """
+        dB_d of Coilset
+
+        Parameters
+        ----------
+        x:
+            The x values at which to calculate the dB_d response
+        z:
+            The z values at which to calculate the dB_d response
+        sum_coils:
+            sum over coils
+        control:
+            operations on control coils only
+
+        Returns
+        -------
+        :
+            Differential of magnetic field
+        """
+        return self.dB_d_response(x, z) * self.current
+
+    def dB_d_response(self, x: float | np.ndarray, z: float | np.ndarray):
+        """
+        Unit dB_d of Coilset
+
+        Parameters
+        ----------
+        x:
+            The x values at which to calculate the dB_d response
+        z:
+            The z values at which to calculate the dB_d response
+        sum_coils:
+            sum over coils
+        control:
+            operations on control coils only
+
+        Returns
+        -------
+        :
+            Differential of magnetic field response
+        """
+        return self._mix_control_method(x, z, greens_dbz_dx, None, disable_analytic=True)
+
     def psi(self, x: float | np.ndarray, z: float | np.ndarray):
         """
         Calculate poloidal flux at (x, z)
@@ -615,6 +664,34 @@
         """
         return self._sum(super().Bz(x, z), sum_coils=sum_coils, control=control)
 
+    def dB_d(
+        self,
+        x: np.ndarray,
+        z: np.ndarray,
+        *,
+        sum_coils: bool = True,
+        control: bool = False,
+    ) -> np.ndarray:
+        """
+        dB_d of Coilset
+
+        Parameters
+        ----------
+        x:
+            The x values at which to calculate the dB_d response
+        z:
+            The z values at which to calculate the dB_d response
+        sum_coils:
+            sum over coils
+        control:
+            operations on control coils only
+
+        Returns
+        -------
+        Differential of magnetic field
+        """
+        return self._sum(super().dB_d(x, z), sum_coils=sum_coils, control=control)
+
     def psi_response(
         self,
         x: np.ndarray,
@@ -682,7 +759,7 @@
         control: bool = False,
     ) -> np.ndarray:
         """
-        Bz of Coilset
+        Unit Bz of Coilset
 
         Parameters
         ----------
@@ -701,6 +778,37 @@
         """
         return self._sum(super().Bz_response(x, z), sum_coils=sum_coils, control=control)
 
+    def dB_d_response(
+        self,
+        x: np.ndarray,
+        z: np.ndarray,
+        *,
+        sum_coils: bool = False,
+        control: bool = False,
+    ) -> np.ndarray:
+        """
+        Unit dB_d of Coilset
+
+        Parameters
+        ----------
+        x:
+            The x values at which to calculate the dB_d response
+        z:
+            The z values at which to calculate the dB_d response
+        sum_coils:
+            sum over coils
+        control:
+            operations on control coils only
+
+        Returns
+        -------
+        :
+            Differential of magnetic field response
+        """
+        return self._sum(
+            super().dB_d_response(x, z), sum_coils=sum_coils, control=control
+        )
+
     def control_F(self, coil_grp: CoilGroup, *, control: bool = False) -> np.ndarray:
         """
         Calculate the coil mutual force

bluemira/equilibria/coils/_grouping.py

@@ -115,7 +115,7 @@
         Bx_analytic: bool = True,
         Bz_analytic: bool = True,
     ):
-        if any(not isinstance(c, Coil | CoilGroup) for c in coils):
+        if any(not isinstance(c, Coil | CoilGroup) for c in coils) or not coils:
             raise TypeError("Not all arguments are a Coil or CoilGroup.")
         self._coils = coils
         self._pad_discretisation(self.__list_getter("_quad_x"))
@@ -337,7 +337,7 @@
                 pfcoils.append(coil)
             coils = pfcoils
             bluemira_warn(
-                "EQDSK coilset empty - dummy coilset in use."
+                "EQDSK coilset empty - dummy coilset in use. "
                 "Please replace with an appropriate coilset."
             )
             return cls(*coils)
@@ -563,14 +563,14 @@
 
         raise KeyError(f"Coil '{name}' not found in Group")
 
-    def _get_coiltype(self, ctype: CoilType | str) -> list[Coil]:
+    def _get_coiltype(self, *ctype: CoilType | str) -> list[Coil]:
         """Find coil by type"""  # noqa: DOC201
         coils = []
-        ctype = CoilType(ctype)
+        ctype = tuple(CoilType(ct) for ct in ctype)
         for c in self._coils:
             if isinstance(c, CoilGroup):
-                coils.extend(c._get_coiltype(ctype))
-            elif c.ctype == ctype:
+                coils.extend(c._get_coiltype(*ctype))
+            elif c.ctype in ctype:
                 coils.append(c)
         return coils
 
@@ -583,14 +583,22 @@
         """
         return [self[n] for n in self.name]
 
-    def get_coiltype(self, ctype: str | CoilType) -> CoilGroup | None:
+    def _get_type_index(self, *ctype: CoilType | str) -> npt.NDArray[int]:
+        coil_type = tuple(CoilType(ct) for ct in ctype) if ctype else CoilType
+
+        return np.asarray(
+            [no for no, coil in enumerate(self.all_coils()) if coil.ctype in coil_type],
+            dtype=int,
+        )
+
+    def get_coiltype(self, *ctype: str | CoilType) -> CoilGroup | None:
         """
         Returns
         -------
         :
             Coils matching coil type
         """
-        if coiltype := self._get_coiltype(ctype):
+        if coiltype := self._get_coiltype(*ctype):
             return CoilGroup(*coiltype)
         return None
 
@@ -724,6 +732,11 @@
         """Get coil max field"""
         return self.__getter("b_max")
 
+    @property
+    def resistance(self) -> np.ndarray:
+        """Get coil resistance"""
+        return self.__getter("resistance")
+
     @property
     def discretisation(self) -> np.ndarray:
         """Get coil discretisations"""
@@ -855,6 +868,11 @@
         """Set coil max fields"""
         self.__setter("b_max", values)
 
+    @resistance.setter
+    def resistance(self, values: float | Iterable[float]):
+        """Set coil resistance"""
+        self.__setter("b_max", values)
+
     @discretisation.setter
     def discretisation(self, values: float | Iterable[float]):
         """Set coil discretisations"""
@@ -1247,6 +1265,18 @@
                 coils.append(c)
         return CoilSet(*coils)
 
+    def get_uncontrolled_coils(self) -> CoilSet:
+        """Get uncontrolled coils"""
+        coils = []
+        for c in self._coils:
+            if isinstance(c, CoilSet):
+                coils.extend(c.get_uncontrolled_coils()._coils)
+            elif (isinstance(c, Coil) and c.name not in self.control) or (
+                isinstance(c, CoilGroup) and not any(n in self.control for n in c.name)
+            ):
+                coils.append(c)
+        return CoilSet(*coils)
+
     @property
     def area(self) -> float:
         """
@@ -1285,10 +1315,10 @@
 
         return np.sum(output[..., inds], axis=-1) if sum_coils else output[..., inds]
 
-    def get_coiltype(self, ctype: str | CoilType) -> CoilSet | None:
+    def get_coiltype(self, *ctype: str | CoilType) -> CoilSet | None:
         """Get coils by coils type"""  # noqa: DOC201
-        if coiltype := self._get_coiltype(ctype):
-            return CoilSet(*coiltype)
+        if coiltype := self._get_coiltype(*ctype):
+            return CoilSet(*coiltype)  # noqa: DOC201
         return None
 
     @classmethod
@@ -1301,11 +1331,7 @@
         """  # noqa: DOC201
         self = super().from_group_vecs(eqdsk)
 
-        self.control = [
-            coil.name
-            for ctype in control_coiltypes
-            for coil in self._get_coiltype(ctype)
-        ]
+        self.control = [coil.name for coil in self._get_coiltype(*control_coiltypes)]
         return self
 
     def get_optimisation_state(