Set preferred normal direction for conversion from points to FourierPlanarCurve

desc/geometry/curve.py

@@ -612,6 +612,9 @@ def __init__(
         assert issubclass(modes.dtype.type, np.integer)
         assert r_n.size == modes.size, "r_n size and modes must be the same size"
         assert basis.lower() in ["xyz", "rpz"]
+        assert not np.allclose(
+            (normal / np.linalg.norm(normal)), [0, 0, -1]
+        ), "normal vector cannot be antiparallel to z-axis"
 
         N = np.max(abs(modes))
         self._r_basis = FourierSeries(N, NFP=1, sym=False)
@@ -806,28 +809,54 @@ def from_values(cls, coords, N=10, basis="xyz", name=""):
 
         # center
         center = np.mean(coords, axis=0)
-        coords = coords - center  # shift to origin
+        coords_centered = coords - center  # shift to origin
 
         # normal
-        U, _, _ = np.linalg.svd(coords.T)
+        U, _, _ = np.linalg.svd(coords_centered.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
+        coords_rotated = coords_centered @ rotmat  # rotate to X-Y plane
 
         warnif(
-            np.max(np.abs(coords[:, 2])) > 1e-14,  # check that Z=0 for all points
+            np.max(np.abs(coords_rotated[:, 2]))
+            > 1e-14,  # check that Z=0 for all points
             UserWarning,
             "Curve values are not planar! Using the projection onto a plane.",
         )
 
+        # polar angle
+        s = np.arctan2(coords_rotated[:, 1], coords_rotated[:, 0])
+        unwrapped_s = np.unwrap(s)
+        # Determine if the sequence is monotonically increasing or decreasing
+        if np.all(np.diff(unwrapped_s) > 0):
+            curve_sign = 1
+        elif np.all(np.diff(unwrapped_s) < 0):
+            curve_sign = -1
+        else:
+            warnings.warn(
+                "The curve parameter s is not strictly increasing or decreasing. "
+                "Assuming default direction.",
+                UserWarning,
+            )
+            curve_sign = 1
+
+        if curve_sign == -1:
+            # original curve was parameterized "backwards" (clockwise),
+            # compared to FourierPlanarCurve assumption
+            normal = -normal  # flip normal vector direction
+            rotmat = rotation_matrix(axis, np.pi)
+            coords_rotated = (
+                coords_rotated @ rotmat
+            )  # flip on X-Y plane to match normal vector
+
         # polar radius and angle
-        r = np.sqrt(coords[:, 0] ** 2 + coords[:, 1] ** 2)
-        s = np.arctan2(coords[:, 1], coords[:, 0])
+        r = np.sqrt(coords_rotated[:, 0] ** 2 + coords_rotated[:, 1] ** 2)
+        s = np.arctan2(coords_rotated[:, 1], coords_rotated[:, 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]

tests/test_coils.py

@@ -252,6 +252,63 @@ def test_convert_type(self):
         np.testing.assert_allclose(B1, B4, rtol=1e-8, atol=1e-8)
         np.testing.assert_allclose(B1, B5, rtol=1e-6, atol=1e-7)
 
+    @pytest.mark.unit
+    def test_SplineXYZ_to_FourierPlanar(self):
+        """Test converting SplineXYZCoil to FourierPlanarCoil object."""
+        # Necessary as SplineXYZ can be either clockwise or counterclockwise
+        npts = 1000
+        N = 50
+
+        # Create a SplineXYZCoil of a planar ellipse
+        s = np.linspace(0, 2 * np.pi, npts)
+        X = 2 * np.cos(s)
+        Y = np.ones(npts)
+        Z = 1 * np.sin(s)
+        c = SplineXYZCoil(X=X, Y=Y, Z=Z, current=1e6)
+
+        # Create a backwards SplineXYZCoil by flipping the coordinates
+        c_backwards = SplineXYZCoil(
+            X=np.flip(X),
+            Y=np.flip(Y),
+            Z=np.flip(Z),
+            current=1e6,
+        )
+
+        # Convert to FourierPlanarCoil
+        c_planar = c.to_FourierPlanar(N=N, grid=npts, basis="xyz")
+        c_backwards_planar = c_backwards.to_FourierPlanar(N=N, grid=npts, basis="xyz")
+
+        grid = LinearGrid(zeta=100)
+
+        field_spline = c.compute_magnetic_field(
+            np.zeros((1, 3)), source_grid=grid, basis="xyz"
+        )
+        field_planar = c_planar.compute_magnetic_field(
+            np.zeros((1, 3)), source_grid=grid, basis="xyz"
+        )
+        field_backwards_spline = c_backwards.compute_magnetic_field(
+            np.zeros((1, 3)), source_grid=grid, basis="xyz"
+        )
+        field_backwards_planar = c_backwards_planar.compute_magnetic_field(
+            np.zeros((1, 3)), source_grid=grid, basis="xyz"
+        )
+
+        np.testing.assert_allclose(
+            field_spline,
+            field_planar,
+            atol=1e-5,
+        )
+        np.testing.assert_allclose(
+            field_backwards_spline,
+            field_backwards_planar,
+            atol=1e-5,
+        )
+        np.testing.assert_allclose(
+            field_planar,
+            -field_backwards_planar,
+            atol=1e-5,
+        )
+
 
 class TestCoilSet:
     """Tests for sets of multiple coils."""

tests/test_curves.py

@@ -657,6 +657,49 @@ def test_asserts(self):
         with pytest.raises(AssertionError):
             _ = FourierPlanarCurve(r_n=[1], modes=[1, 2])
 
+    @pytest.mark.unit
+    def test_to_FourierPlanarCurve(self):
+        """Test converting SplineXYZCurve to FourierPlanarCurve object."""
+        npts = 1000
+        N = 5
+
+        # Create a SplineXYZCurve of a planar circle
+        s = np.linspace(0, 2 * np.pi, npts)
+        X = 2 * np.cos(s)
+        Y = np.ones(npts)
+        Z = 2 * np.sin(s)
+        c = SplineXYZCurve(X=X, Y=Y, Z=Z)
+
+        # Create a backwards SplineXYZCurve by flipping the coordinates
+        c_backwards = SplineXYZCurve(
+            X=np.flip(X),
+            Y=np.flip(Y),
+            Z=np.flip(Z),
+        )
+
+        # Convert to FourierPlanarCurve
+        c_planar = c.to_FourierPlanar(N=N, grid=npts, basis="xyz")
+        c_backwards_planar = c_backwards.to_FourierPlanar(N=N, grid=npts, basis="xyz")
+
+        grid = LinearGrid(N=20, endpoint=True)
+
+        coords_spline = c.compute("x", grid=grid)["x"]
+        coords_planar = c_planar.compute("x", grid=grid)["x"]
+        coords_backwards_spline = c_backwards.compute("x", grid=grid)["x"]
+        coords_backwards_planar = c_backwards_planar.compute("x", grid=grid)["x"]
+
+        # Assertions for point positions
+        np.testing.assert_allclose(coords_spline, coords_planar, atol=1e-10)
+        np.testing.assert_allclose(
+            coords_backwards_spline,
+            coords_backwards_planar,
+            atol=1e-10,
+        )
+        # backwards coordinate order is important for Biot-Savart current
+        np.testing.assert_allclose(
+            coords_planar, np.flip(coords_backwards_planar, axis=0), atol=1e-10
+        )
+
 
 class TestSplineXYZCurve:
     """Tests for SplineXYZCurve class."""