Set preferred normal direction for conversion from points to FourierPlanarCurve

[ryanwu4]

Previously, the from_values() method in FourierPlanarCurve would convert the coords provided to an r(s) representation with s sorted to be monotonically increasing (counterclockwise). However, this means that if the original coordinates were parameterized "backwards" (clockwise), the new FourierPlanarCurve would return coordinates in a reversed order, which is problematic if you are trying to run something like a Biot Savart integral where the coordinate order is the implicit current direction in the loop. This PR fixes the issue by flipping the normal vector of the FourierPlanarCurve when the provided coordinates are in the clockwise direction.

[codecov]

Codecov Report

Attention: Patch coverage is 61.53846% with 5 lines in your changes missing coverage. Please review.

Project coverage is 95.58%. Comparing base (44bac68) to head (8fd19b2).
Report is 3067 commits behind head on master.

Files with missing lines	Patch %	Lines
desc/geometry/curve.py	61.53%	5 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1463      +/-   ##
==========================================
+ Coverage   95.31%   95.58%   +0.26%     
==========================================
  Files          87       98      +11     
  Lines       22262    25165    +2903     
==========================================
+ Hits        21220    24053    +2833     
- Misses       1042     1112      +70     

Files with missing lines	Coverage Δ
desc/geometry/curve.py	95.66% <61.53%> (-0.34%)	⬇️

... and 77 files with indirect coverage changes

[ddudt]

Why do we need the additional rotation and not just flipping the normal vector direction? I think we only need to flip the normal vector, since that would effectively reverse the parameterization in real space

[ddudt]

Alternatively, I think we could flip the sign of the Fourier coefficients r_n where $m < 0$, since $\cos(-\theta)=\cos(\theta)$ and $\sin(-\theta)=-\sin(\theta)$

[ryanwu4]

The weird thing is that the curve points test passes with just flipping the sign of the normal, but the fields are different between spline and normal-flipped-FourierPlanar unless I redo the rotation matrix with the flipped normal. I'll look into this a bit more.

[ryanwu4]

So I looked at the coil shapes with and without re-rotating, and it appears that if you only flip the sign of the normal without rotating as well, you mirror the curve as well -- this is because you swap the chirality of the curve. If you drew the curve on a piece of paper, flipping the normal would be looking at the sheet of paper from the back, rather than flipping the piece of paper itself.

[ddudt]

Can you try flipping the sign of the r_n coefficients? (Only the $\sin$ terms I think.) I think that would be a cleaner solution than flipping the normal + rotation, if it works.

[ryanwu4]

For some reason this doesn't work, no matter what it seems to need both an axis reversal and a rotation (pancake flip) to get the curves to line up again -- I clarified in the code what the rotation is now doing, but I couldn't find the best way to do it without another rotation.

[ryanwu4]

In the process of working on this , I noticed the following behavior:
If the normal vector of a FourierPlanarCurve is defined as [0,0,-1], the curve direction will be incorrect (clockwise about the normal instead of counterclockwise). This is due to the axis term in the corresponding x compute function being undefined for axes parallel to the Z axis, so the rotation matrix defaults to the unit diagonal matrix. I added this as an assertion check, though conceivably this can probably be fixed.

The problematic code in the compute function:

    Zaxis = jnp.array([0.0, 0.0, 1.0])  # 2D curve in X-Y plane has normal = +Z axis
    axis = cross(Zaxis, normal)
    angle = jnp.arccos(dot(Zaxis, safenormalize(normal)))
    A = rotation_matrix(axis=axis, angle=angle)

Minimal Example:

coil = FourierPlanarCoil(
            current=1e6,
            center=[0, 0, 0],
            normal=[0, 0, 1],
            r_n=2,
            modes=[0],
        )

coil2 = FourierPlanarCoil(
            current=1e6,
            center=[0, 0, 0],
            normal=[0, 0, -1],
            r_n=2,
            modes=[0],
        )

field1 = coil.compute_magnetic_field([0, 0, 0], source_grid=LinearGrid(N=100), basis = "xyz")
field2 = coil2.compute_magnetic_field([0, 0, 0], source_grid=LinearGrid(N=100), basis = "xyz")

np.testing.assert_allclose(field1, field2, rtol=1e-3, atol=1e-3) # these fields shouldn't be equal but they are

The expected behavior of the assertion should be to fail, as field1 and field2 should have equal and opposite magnitudes. Instead, they are actually the same number, and the assertion succeeds. Compare this to if you perturb the normal vectors just a tiny bit (eg. [0, 0.1, +-1]), and the field correctly has opposite signs.

[dpanici]

In the process of working on this , I noticed the following behavior: If the normal vector of a FourierPlanarCurve is defined as [0,0,-1], the curve direction will be incorrect (clockwise about the normal instead of counterclockwise). This is due to the axis term in the corresponding x compute function being undefined for axes parallel to the Z axis, so the rotation matrix defaults to the unit diagonal matrix. I added this as an assertion check, though conceivably this can probably be fixed.

The problematic code in the compute function:

    Zaxis = jnp.array([0.0, 0.0, 1.0])  # 2D curve in X-Y plane has normal = +Z axis
    axis = cross(Zaxis, normal)
    angle = jnp.arccos(dot(Zaxis, safenormalize(normal)))
    A = rotation_matrix(axis=axis, angle=angle)

Minimal Example:

coil = FourierPlanarCoil(
            current=1e6,
            center=[0, 0, 0],
            normal=[0, 0, 1],
            r_n=2,
            modes=[0],
        )

coil2 = FourierPlanarCoil(
            current=1e6,
            center=[0, 0, 0],
            normal=[0, 0, -1],
            r_n=2,
            modes=[0],
        )

field1 = coil.compute_magnetic_field([0, 0, 0], source_grid=LinearGrid(N=100), basis = "xyz")
field2 = coil2.compute_magnetic_field([0, 0, 0], source_grid=LinearGrid(N=100), basis = "xyz")

np.testing.assert_allclose(field1, field2, rtol=1e-3, atol=1e-3) # these fields shouldn't be equal but they are

The expected behavior of the assertion should be to fail, as field1 and field2 should have equal and opposite magnitudes. Instead, they are actually the same number, and the assertion succeeds. Compare this to if you perturb the normal vectors just a tiny bit (eg. [0, 0.1, +-1]), and the field correctly has opposite signs.

related to #1456