Repository: PlasmaControl/DESC

Issue:
Add `FourierPlanarCurve.from_values` Method

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 75d36322a..cd3e219fc 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,6 +1,11 @@
 Changelog
 =========
 
+New Features
+
+- adds ``from_values`` method that was present in ``FourierRZCurve`` but missing in ``FourierRZCoil``
+- also adds new ``from_values`` method for ``FourierPlanarCurve`` and ``FourierPlanarCoil``
+
 
 v0.12.0
 -------
diff --git a/desc/coils.py b/desc/coils.py
index efea3b0ee..49c56d596 100644
--- a/desc/coils.py
+++ b/desc/coils.py
@@ -346,6 +346,75 @@ class _Coil(_MagneticField, Optimizable, ABC):
             self.current, coords, knots=knots, method=method, name=name, basis="xyz"
         )
 
+    def to_FourierRZ(self, N=10, grid=None, NFP=None, sym=False, name=""):
+        """Convert Coil to FourierRZCoil representation.
+
+        Note that some types of coils may not be representable in this basis.
+
+        Parameters
+        ----------
+        N : int
+            Fourier resolution of the new R,Z representation.
+        grid : Grid, int or None
+            Grid used to evaluate curve coordinates on to fit with FourierRZCoil.
+            If an integer, uses that many equally spaced points.
+        NFP : int
+            Number of field periods, the coil will have a discrete toroidal symmetry
+            according to NFP.
+        sym : bool, optional
+            whether the curve is stellarator-symmetric or not. Default
+            is False.
+        name : str
+            name for this coil
+
+        Returns
+        -------
+        curve : FourierRZCoil
+            New representation of the coil parameterized by Fourier series for R,Z.
+
+        """
+        NFP = 1 or NFP
+        if grid is None:
+            grid = LinearGrid(N=2 * N + 1)
+        coords = self.compute("x", grid=grid, basis="xyz")["x"]
+        return FourierRZCoil.from_values(
+            self.current, coords, N=N, NFP=NFP, basis="xyz", sym=sym, name=name
+        )
+
+    def to_FourierPlanar(self, N=10, grid=None, basis="xyz", name=""):
+        """Convert Coil to FourierPlanarCoil representation.
+
+        Note that some types of coils may not be representable in this basis.
+        In this case, a least-squares fit will be done to find the
+        planar coil that best represents the coil.
+
+        Parameters
+        ----------
+        N : int
+            Fourier resolution of the new FourierPlanarCoil representation.
+        grid : Grid, int or None
+            Grid used to evaluate curve coordinates on to fit with FourierPlanarCoil.
+            If an integer, uses that many equally spaced points.
+        basis : {'xyz', 'rpz'}
+            Coordinate system for center and normal vectors. Default = 'xyz'.
+        name : str
+            name for this coil
+
+        Returns
+        -------
+        coil : FourierPlanarCoil
+            New representation of the coil parameterized by Fourier series for
+            minor radius r in a plane specified by a center position and normal
+            vector.
+
+        """
+        if grid is None:
+            grid = LinearGrid(N=2 * N + 1)
+        coords = self.compute("x", grid=grid, basis=basis)["x"]
+        return FourierPlanarCoil.from_values(
+            self.current, coords, N=N, basis=basis, name=name
+        )
+
 
 class FourierRZCoil(_Coil, FourierRZCurve):
     """Coil parameterized by fourier series for R,Z in terms of toroidal angle phi.
@@ -353,7 +422,7 @@ class FourierRZCoil(_Coil, FourierRZCurve):
     Parameters
     ----------
     current : float
-        current through coil, in Amperes
+        Current through the coil, in Amperes.
     R_n, Z_n: array-like
         fourier coefficients for R, Z
     modes_R : array-like
@@ -414,6 +483,50 @@ class FourierRZCoil(_Coil, FourierRZCurve):
     ):
         super().__init__(current, R_n, Z_n, modes_R, modes_Z, NFP, sym, name)
 
+    @classmethod
+    def from_values(cls, current, coords, N=10, NFP=1, basis="rpz", sym=False, name=""):
+        """Fit coordinates to FourierRZCoil representation.
+
+        Parameters
+        ----------
+        current : float
+            Current through the coil, in Amperes.
+        coords: ndarray, shape (num_coords,3)
+            coordinates to fit a FourierRZCurve object with each column
+            corresponding to xyz or rpz depending on the basis argument.
+        N : int
+            Fourier resolution of the new R,Z representation.
+        NFP : int
+            Number of field periods, the curve will have a discrete toroidal symmetry
+            according to NFP.
+        basis : {"rpz", "xyz"}
+            basis for input coordinates. Defaults to "rpz"
+        sym : bool
+            Whether to enforce stellarator symmetry.
+        name : str
+            name for this coil
+
+
+        Returns
+        -------
+        coil : FourierRZCoil
+            New representation of the coil parameterized by Fourier series for R,Z.
+
+        """
+        curve = super().from_values(
+            coords=coords, N=N, NFP=NFP, basis=basis, sym=sym, name=name
+        )
+        return FourierRZCoil(
+            current=current,
+            R_n=curve.R_n,
+            Z_n=curve.Z_n,
+            modes_R=curve.R_basis.modes[:, 2],
+            modes_Z=curve.Z_basis.modes[:, 2],
+            NFP=NFP,
+            sym=curve.sym,
+            name=name,
+        )
+
 
 class FourierXYZCoil(_Coil, FourierXYZCurve):
     """Coil parameterized by fourier series for X,Y,Z in terms of arbitrary angle s.
@@ -421,7 +534,7 @@ class FourierXYZCoil(_Coil, FourierXYZCurve):
     Parameters
     ----------
     current : float
-        current through coil, in Amperes
+        Current through the coil, in Amperes.
     X_n, Y_n, Z_n: array-like
         fourier coefficients for X, Y, Z
     modes : array-like
@@ -484,7 +597,7 @@ class FourierXYZCoil(_Coil, FourierXYZCurve):
         Parameters
         ----------
         current : float
-            Current through the coil, in Amps.
+            Current through the coil, in Amperes.
         coords: ndarray
             Coordinates to fit a FourierXYZCoil object with.
         N : int
@@ -503,12 +616,13 @@ class FourierXYZCoil(_Coil, FourierXYZCurve):
             New representation of the coil parameterized by Fourier series for X,Y,Z.
 
         """
-        curve = super().from_values(coords, N, s, basis)
-        return cls(
-            current,
+        curve = super().from_values(coords=coords, N=N, s=s, basis=basis, name=name)
+        return FourierXYZCoil(
+            current=current,
             X_n=curve.X_n,
             Y_n=curve.Y_n,
             Z_n=curve.Z_n,
+            modes=curve.X_basis.modes[:, 2],
             name=name,
         )
 
@@ -586,6 +700,41 @@ class FourierPlanarCoil(_Coil, FourierPlanarCurve):
     ):
         super().__init__(current, center, normal, r_n, modes, basis, name)
 
+    @classmethod
+    def from_values(cls, current, coords, N=10, basis="xyz", name=""):
+        """Fit coordinates to FourierPlanarCoil representation.
+
+        Parameters
+        ----------
+        current : float
+            Current through the coil, in Amperes.
+        coords: ndarray, shape (num_coords,3)
+            Coordinates to fit a FourierPlanarCurve object with each column
+            corresponding to xyz or rpz depending on the basis argument.
+        N : int
+            Fourier resolution of the new r representation.
+        basis : {"rpz", "xyz"}
+            Basis for input coordinates. Defaults to "xyz".
+        name : str
+            Name for this curve.
+
+        Returns
+        -------
+        curve : FourierPlanarCoil
+            New representation of the coil parameterized by a Fourier series for r.
+
+        """
+        curve = super().from_values(coords=coords, N=N, basis=basis, name=name)
+        return FourierPlanarCoil(
+            current=current,
+            center=curve.center,
+            normal=curve.normal,
+            r_n=curve.r_n,
+            modes=curve.r_basis.modes[:, 2],
+            basis="xyz",
+            name=name,
+        )
+
 
 class SplineXYZCoil(_Coil, SplineXYZCurve):
     """Coil parameterized by spline points in X,Y,Z.
@@ -593,7 +742,7 @@ class SplineXYZCoil(_Coil, SplineXYZCurve):
     Parameters
     ----------
     current : float
-        current through coil, in Amperes
+        Current through the coil, in Amperes.
     X, Y, Z: array-like
         Points for X, Y, Z describing the curve. If the endpoint is included
         (ie, X[0] == X[-1]), then the final point will be dropped.
@@ -705,7 +854,7 @@ class SplineXYZCoil(_Coil, SplineXYZCurve):
         Parameters
         ----------
         current : float
-            Current through the coil, in Amps.
+            Current through the coil, in Amperes.
         coords: ndarray
             Points for X, Y, Z describing the curve. If the endpoint is included
             (ie, X[0] == X[-1]), then the final point will be dropped.
@@ -737,9 +886,11 @@ class SplineXYZCoil(_Coil, SplineXYZCurve):
             New representation of the coil parameterized by splines in X,Y,Z.
 
         """
-        curve = super().from_values(coords, knots, method, basis=basis)
-        return cls(
-            current,
+        curve = super().from_values(
+            coords=coords, knots=knots, method=method, basis=basis, name=name
+        )
+        return SplineXYZCoil(
+            current=current,
             X=curve.X,
             Y=curve.Y,
             Z=curve.Z,
diff --git a/desc/geometry/core.py b/desc/geometry/core.py
index fd754dfe8..a389eb444 100644
--- a/desc/geometry/core.py
+++ b/desc/geometry/core.py
@@ -209,7 +209,7 @@ class Curve(IOAble, Optimizable, ABC):
             + " (name={})".format(self.name)
         )
 
-    def to_FourierXYZ(self, N=None, grid=None, s=None, name=""):
+    def to_FourierXYZ(self, N=10, grid=None, s=None, name=""):
         """Convert Curve to FourierXYZCurve representation.
 
         Parameters
@@ -284,7 +284,7 @@ class Curve(IOAble, Optimizable, ABC):
             coords, knots=knots, method=method, name=name, basis="xyz"
         )
 
-    def to_FourierRZ(self, N=None, grid=None, NFP=None, name=""):
+    def to_FourierRZ(self, N=10, grid=None, NFP=None, sym=False, name=""):
         """Convert Curve to FourierRZCurve representation.
 
         Note that some types of curves may not be representable in this basis.
@@ -299,6 +299,8 @@ class Curve(IOAble, Optimizable, ABC):
         NFP : int
             Number of field periods, the curve will have a discrete toroidal symmetry
             according to NFP.
+        sym : bool, optional
+            Whether the curve is stellarator-symmetric or not. Default is False.
         name : str
             name for this curve
 
@@ -312,9 +314,45 @@ class Curve(IOAble, Optimizable, ABC):
 
         NFP = 1 or NFP
         if grid is None:
-            grid = LinearGrid(N=2 * N + 1, NFP=NFP)
+            grid = LinearGrid(N=2 * N + 1)
         coords = self.compute("x", grid=grid, basis="xyz")["x"]
-        return FourierRZCurve.from_values(coords, N=N, NFP=NFP, basis="xyz", name=name)
+        return FourierRZCurve.from_values(
+            coords, N=N, NFP=NFP, basis="xyz", name=name, sym=sym
+        )
+
+    def to_FourierPlanar(self, N=10, grid=None, basis="xyz", name=""):
+        """Convert Curve to FourierPlanarCurve representation.
+
+        Note that some types of curves may not be representable in this basis.
+        In this case, a least-squares fit will be done to find the
+        planar curve that best represents the curve.
+
+        Parameters
+        ----------
+        N : int
+            Fourier resolution of the new FourierPlanarCurve representation.
+        grid : Grid, int or None
+            Grid used to evaluate curve coordinates on to fit with FourierPlanarCurve.
+            If an integer, uses that many equally spaced points.
+        basis : {'xyz', 'rpz'}
+            Coordinate system for center and normal vectors. Default = 'xyz'.
+        name : str
+            name for this curve
+
+        Returns
+        -------
+        curve : FourierPlanarCurve
+            New representation of the curve parameterized by Fourier series for
+            minor radius r in a plane specified by a center position and normal
+            vector.
+
+        """
+        from .curve import FourierPlanarCurve
+
+        if grid is None:
+            grid = LinearGrid(N=2 * N + 1)
+        coords = self.compute("x", grid=grid, basis=basis)["x"]
+        return FourierPlanarCurve.from_values(coords, N=N, basis=basis, name=name)
 
 
 class Surface(IOAble, Optimizable, ABC):
diff --git a/desc/geometry/curve.py b/desc/geometry/curve.py
index a62ed8fda..e2c583ea6 100644
--- a/desc/geometry/curve.py
+++ b/desc/geometry/curve.py
@@ -6,12 +6,13 @@ import numpy as np
 
 from desc.backend import jnp, put
 from desc.basis import FourierSeries
-from desc.compute import rpz2xyz, xyz2rpz
+from desc.compute import rpz2xyz, rpz2xyz_vec, xyz2rpz, xyz2rpz_vec
+from desc.compute.geom_utils import rotation_matrix
 from desc.grid import LinearGrid
 from desc.io import InputReader
 from desc.optimizable import optimizable_parameter
 from desc.transform import Transform
-from desc.utils import check_nonnegint, check_posint, copy_coeffs, errorif
+from desc.utils import check_nonnegint, check_posint, copy_coeffs, errorif, warnif
 
 from .core import Curve
 
@@ -230,7 +231,7 @@ class FourierRZCurve(Curve):
         return curve
 
     @classmethod
-    def from_values(cls, coords, N=10, NFP=1, basis="rpz", name="", sym=False):
+    def from_values(cls, coords, N=10, NFP=1, sym=False, basis="rpz", name=""):
         """Fit coordinates to FourierRZCurve representation.
 
         Parameters
@@ -243,8 +244,12 @@ class FourierRZCurve(Curve):
         NFP : int
             Number of field periods, the curve will have a discrete toroidal symmetry
             according to NFP.
+        sym : bool
+            Whether to enforce stellarator symmetry.
         basis : {"rpz", "xyz"}
             basis for input coordinates. Defaults to "rpz"
+        name : str
+            Name for this curve.
 
         Returns
         -------
@@ -287,7 +292,15 @@ class FourierRZCurve(Curve):
         transform = Transform(grid, basis, build_pinv=True)
         R_n = transform.fit(R)
         Z_n = transform.fit(Z)
-        return FourierRZCurve(R_n=R_n, Z_n=Z_n, NFP=NFP, name=name, sym=sym)
+        return FourierRZCurve(
+            R_n=R_n,
+            Z_n=Z_n,
+            modes_R=basis.modes[:, 2],
+            modes_Z=basis.modes[:, 2],
+            NFP=NFP,
+            sym=sym,
+            name=name,
+        )
 
 
 def _unclose_curve(X, Y, Z):
@@ -489,17 +502,18 @@ class FourierXYZCurve(Curve):
         Parameters
         ----------
         coords: ndarray, shape (num_coords,3)
-            coordinates to fit a FourierXYZCurve object, with each column
+            Coordinates to fit a FourierXYZCurve object, with each column
             corresponding to xyz or rpz depending on the basis argument.
         N : int
             Fourier resolution of the new X,Y,Z representation.
         s : ndarray or "arclength"
-            arbitrary curve parameter to use for the fitting.
+            Arbitrary curve parameter to use for the fitting.
             Should be monotonic, 1D array of same length as
             coords. if None, defaults linearly spaced in [0,2pi)
             Alternative, can pass "arclength" to use normalized distance between points.
         basis : {"rpz", "xyz"}
-            basis for input coordinates. Defaults to "xyz"
+            Basis for input coordinates. Defaults to "xyz".
+
         Returns
         -------
         curve : FourierXYZCurve
@@ -547,7 +561,9 @@ class FourierXYZCurve(Curve):
         X_n = transform.fit(X)
         Y_n = transform.fit(Y)
         Z_n = transform.fit(Z)
-        return FourierXYZCurve(X_n=X_n, Y_n=Y_n, Z_n=Z_n, name=name)
+        return FourierXYZCurve(
+            X_n=X_n, Y_n=Y_n, Z_n=Z_n, modes=basis.modes[:, 2], name=name
+        )
 
 
 class FourierPlanarCurve(Curve):
@@ -637,7 +653,7 @@ class FourierPlanarCurve(Curve):
         else:
             raise ValueError(
                 "center should be a 3 element vector in "
-                + self._basis
+                + self.basis
                 + " coordinates, got {}".format(new)
             )
 
@@ -654,7 +670,7 @@ class FourierPlanarCurve(Curve):
         else:
             raise ValueError(
                 "normal should be a 3 element vector in "
-                + self._basis
+                + self.basis
                 + " coordinates, got {}".format(new)
             )
 
@@ -674,6 +690,23 @@ class FourierPlanarCurve(Curve):
                 + f"basis with {self.r_basis.num_modes} modes."
             )
 
+    @property
+    def basis(self):
+        """Coordinate system for center and normal vectors."""
+        return self._basis
+
+    @basis.setter
+    def basis(self, new):
+        assert new.lower() in ["xyz", "rpz"]
+        if new != self.basis:
+            if new == "xyz":
+                self.normal = rpz2xyz_vec(self.normal, phi=self.center[1])
+                self.center = rpz2xyz(self.center)
+            else:
+                self.center = xyz2rpz(self.center)
+                self.normal = xyz2rpz_vec(self.normal, phi=self.center[1])
+            self._basis = new
+
     def get_coeffs(self, n):
         """Get Fourier coefficients for given mode number(s)."""
         n = np.atleast_1d(n).astype(int)
@@ -740,10 +773,81 @@ class FourierPlanarCurve(Curve):
             transforms=transforms,
             data=data,
             override_grid=override_grid,
-            basis_in=self._basis,
+            basis_in=self.basis,
             **kwargs,
         )
 
+    @classmethod
+    def from_values(cls, coords, N=10, basis="xyz", name=""):
+        """Fit coordinates to FourierPlanarCurve representation.
+
+        Parameters
+        ----------
+        coords: ndarray, shape (num_coords,3)
+            Coordinates to fit a FourierPlanarCurve object with each column
+            corresponding to xyz or rpz depending on the basis argument.
+        N : int
+            Fourier resolution of the new r representation.
+        basis : {"rpz", "xyz"}
+            Basis for input coordinates. Defaults to "xyz".
+        name : str
+            Name for this curve.
+
+        Returns
+        -------
+        curve : FourierPlanarCurve
+            New representation of the curve parameterized by a Fourier series for r.
+
+        """
+        # convert to xyz basis
+        if basis == "rpz":
+            coords = rpz2xyz(coords)
+        coords = np.atleast_2d(coords)
+
+        # center
+        center = np.mean(coords, axis=0)
+        coords = coords - center  # shift to origin
+
+        # normal
+        U, _, _ = np.linalg.svd(coords.T)
+        normal = U[:, -1].T  # left singular vector of the least singular value
+
+        # axis and angle of rotation
+        Z_axis = np.array([0, 0, 1])
+        axis = np.cross(Z_axis, normal)
+        angle = np.arccos(np.dot(Z_axis, normal))
+        rotmat = rotation_matrix(axis, angle)
+        coords = coords @ rotmat  # rotate to X-Y plane
+
+        warnif(
+            np.max(np.abs(coords[:, 2])) > 1e-14,  # check that Z=0 for all points
+            UserWarning,
+            "Curve values are not planar! Using the projection onto a plane.",
+        )
+
+        # polar radius and angle
+        r = np.sqrt(coords[:, 0] ** 2 + coords[:, 1] ** 2)
+        s = np.arctan2(coords[:, 1], coords[:, 0])
+        s = np.mod(s + 2 * np.pi, 2 * np.pi)  # mod angle to range [0, 2*pi)
+        idx = np.argsort(s)  # sort angle to be monotonically increasing
+        r = r[idx]
+        s = s[idx]
+
+        # Fourier transform
+        basis = FourierSeries(N, NFP=1, sym=False)
+        grid_fit = LinearGrid(zeta=s, NFP=1)
+        transform_fit = Transform(grid_fit, basis, build_pinv=True)
+        r_n = transform_fit.fit(r)
+
+        return FourierPlanarCurve(
+            center=center,
+            normal=normal,
+            r_n=r_n,
+            modes=basis.modes[:, 2],
+            basis="xyz",
+            name=name,
+        )
+
 
 class SplineXYZCurve(Curve):
     """Curve parameterized by spline knots in X,Y,Z.
@@ -969,7 +1073,7 @@ class SplineXYZCurve(Curve):
         )
 
     @classmethod
-    def from_values(cls, coords, knots=None, method="cubic", name="", basis="xyz"):
+    def from_values(cls, coords, knots=None, method="cubic", basis="xyz", name=""):
         """Create SplineXYZCurve from coordinate values.
 
         Parameters
@@ -980,7 +1084,7 @@ class SplineXYZCurve(Curve):
             endpoint is included (ie, X[0] == X[-1]), then the final point will be
             dropped.
         knots : ndarray
-            arbitrary curve parameter values to use for spline knots,
+            Arbitrary curve parameter values to use for spline knots,
             should be an 1D ndarray of same length as the input.
             (input length in this case is determined by grid argument, since
             the input coordinates come from
@@ -996,10 +1100,10 @@ class SplineXYZCurve(Curve):
             - `'cubic2'`: C2 cubic splines (aka natural splines)
             - `'catmull-rom'`: C1 cubic centripetal "tension" splines
 
-        name : str
-            name for this curve
         basis : {"rpz", "xyz"}
-            basis for input coordinates. Defaults to "xyz"
+            Basis for input coordinates. Defaults to "xyz".
+        name : str
+            Name for this curve.
 
         Returns
         -------
@@ -1010,5 +1114,10 @@ class SplineXYZCurve(Curve):
         if basis == "rpz":
             coords = rpz2xyz(coords)
         return SplineXYZCurve(
-            coords[:, 0], coords[:, 1], coords[:, 2], knots, method, name
+            X=coords[:, 0],
+            Y=coords[:, 1],
+            Z=coords[:, 2],
+            knots=knots,
+            method=method,
+            name=name,
         )


Test patch:
diff --git a/tests/test_coils.py b/tests/test_coils.py
index defa278c3..8d8b9c7d4 100644
--- a/tests/test_coils.py
+++ b/tests/test_coils.py
@@ -209,13 +209,22 @@ class TestCoil:
     def test_converting_coil_types(self):
         """Test conversions between coil representations."""
         s = np.linspace(0, 2 * np.pi, 100, endpoint=False)
-        coil1 = FourierRZCoil(1e6, [0, 10, 2], [-2, 0, 0])
+        coil1 = FourierRZCoil(1e6, [0, 10, 1], [0, 0, 0])
         coil2 = coil1.to_FourierXYZ(s=s)
         coil3 = coil1.to_SplineXYZ(knots=s)
+        coil4 = coil1.to_FourierRZ(N=coil1.N)
+        coil5 = coil1.to_FourierPlanar(N=10, basis="rpz")
+
         grid = LinearGrid(zeta=s)
         x1 = coil1.compute("x", grid=grid, basis="xyz")["x"]
         x2 = coil2.compute("x", grid=grid, basis="xyz")["x"]
         x3 = coil3.compute("x", grid=grid, basis="xyz")["x"]
+        x4 = coil4.compute("x", grid=grid, basis="xyz")["x"]
+        grid = LinearGrid(
+            zeta=np.arctan2(x1[:, 1] - coil5.center[1], x1[:, 0] - coil5.center[0])
+        )  # zeta = arctan angle instead of s for the planar coil for the same points
+        x5 = coil5.compute("x", grid=grid, basis="xyz")["x"]
+
         B1 = coil1.compute_magnetic_field(
             np.zeros((1, 3)), source_grid=grid, basis="xyz"
         )
@@ -225,10 +234,21 @@ class TestCoil:
         B3 = coil3.compute_magnetic_field(
             np.zeros((1, 3)), source_grid=grid, basis="xyz"
         )
+        B4 = coil4.compute_magnetic_field(
+            np.zeros((1, 3)), source_grid=grid, basis="xyz"
+        )
+        B5 = coil5.compute_magnetic_field(
+            np.zeros((1, 3)), source_grid=grid, basis="xyz"
+        )
+
         np.testing.assert_allclose(x1, x2, atol=1e-12)
         np.testing.assert_allclose(x1, x3, atol=1e-12)
+        np.testing.assert_allclose(x1, x4, atol=1e-12)
+        np.testing.assert_allclose(x1, x5, atol=1e-12)
         np.testing.assert_allclose(B1, B2, rtol=1e-8, atol=1e-8)
         np.testing.assert_allclose(B1, B3, rtol=1e-3, atol=1e-8)
+        np.testing.assert_allclose(B1, B4, rtol=1e-8, atol=1e-8)
+        np.testing.assert_allclose(B1, B5, rtol=1e-6, atol=1e-7)
 
 
 class TestCoilSet: