diff --git a/desc/geometry/curve.py b/desc/geometry/curve.py index e2c583ea6f..a6896e38ae 100644 --- a/desc/geometry/curve.py +++ b/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,10 +809,10 @@ 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 @@ -817,17 +820,43 @@ def from_values(cls, coords, N=10, basis="xyz", name=""): 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] diff --git a/tests/test_coils.py b/tests/test_coils.py index 704ad5f761..ccf07b2f7f 100644 --- a/tests/test_coils.py +++ b/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.""" diff --git a/tests/test_curves.py b/tests/test_curves.py index 74fd41aa98..37ab5f18c9 100644 --- a/tests/test_curves.py +++ b/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."""