poly_to_parametric capable of flipping the ellipse

[ben-cassese]

As discovered by @tigerchenlu98, the poly_to_parametric function in squishyplanet.engine.parametric_ellipse, which converts the $\rho$ coefficients that describe the ellipse in an implicit form to $c$ coefficients that describe it in a parametric form, can return an ellipse that's rotated by 90 degrees in the sky plane. Thankfully, it looks like this does not happen when using it within the usual squishyplanet workflow: I can't find a set of $\rho$ coefficients generated after creating an OblateSystem that do not map to the correct $c$ coefficients. But, it's possible to fool the function when providing coefficients that were generated externally:

import jax

jax.config.update("jax_enable_x64", True)
import jax.numpy as jnp
import numpy as np
import matplotlib.pyplot as plt
from squishyplanet.engine.parametric_ellipse import (
    poly_to_parametric,
    poly_to_parametric_helper,
)
from squishyplanet import OblateSystem

def ellipse_poly_eq(rho_xx, rho_xy, rho_x0, rho_yy, rho_y0, rho_00, x, y):
    return (
        rho_xx * x**2
        + rho_xy * x * y
        + rho_x0 * x
        + rho_yy * y**2
        + rho_y0 * y
        + rho_00
    )

def equator_poly_coeffs(theta, phi, r, xp, yp):
    q_x1 = np.sin(theta) ** 2
    q_x2 = 2 * np.cos(theta) * np.sin(theta) * np.sin(phi)
    q_x3 = -2 * xp * np.sin(theta) ** 2 - 2 * yp * np.cos(theta) * np.sin(
        theta
    ) * np.sin(phi)
    q_y1 = np.cos(phi) ** 2 + np.cos(theta) ** 2 * np.sin(phi) ** 2
    q_y2 = (
        -2 * yp * np.cos(phi) ** 2
        - 2 * xp * np.cos(theta) * np.sin(theta) * np.sin(phi)
        - 2 * yp * np.cos(theta) ** 2 * np.sin(phi) ** 2
    )
    q_y3 = 0

    norm = (
        r**2 * np.sin(theta) ** 2 * np.cos(phi) ** 2
        - yp**2 * np.cos(phi) ** 2
        - xp**2 * np.sin(theta) ** 2
        - 2 * xp * yp * np.cos(theta) * np.sin(theta) * np.sin(phi)
        - yp**2 * np.cos(theta) ** 2 * np.sin(phi) ** 2
    )

    return {
        "rho_xx": q_x1 / norm,
        "rho_xy": q_x2 / norm,
        "rho_x0": q_x3 / norm,
        "rho_yy": q_y1 / norm,
        "rho_y0": q_y2 / norm,
        "rho_00": q_y3 / norm,
    }

def parametric_equator(theta, phi, r, xp, yp):
    equator_2d_coeffs = equator_poly_coeffs(theta=theta, phi=phi, r=r, xp=xp, yp=yp)

    coeffs = poly_to_parametric(**equator_2d_coeffs)

    alpha = jnp.linspace(0, 2 * jnp.pi, 100)
    x_para = (
        coeffs["c_x1"] * np.cos(alpha) + coeffs["c_x2"] * np.sin(alpha) + coeffs["c_x3"]
    )
    y_para = (
        coeffs["c_y1"] * np.cos(alpha) + coeffs["c_y2"] * np.sin(alpha) + coeffs["c_y3"]
    )
    return x_para, y_para

equator_params = {
    "theta": jnp.pi / 6,
    "phi": 1.0,
    "r": 0.3,
    "xp": 0.16,
    "yp": 0.0,
}

xs = np.linspace(-1.5, 1.5, 100)
ys = np.linspace(-1.5, 1.5, 100)
XS, YS = np.meshgrid(xs, ys)

fig, ax = plt.subplots(figsize=(10, 5), ncols=2)
equator_2d_coeffs = equator_poly_coeffs(**equator_params)
ax[0].contour(
    XS, YS, ellipse_poly_eq(**equator_2d_coeffs, x=XS, y=YS), levels=[1], colors="teal"
)
x_para, y_para = parametric_equator(**equator_params)
ax[0].plot(x_para, y_para, color="red", label="Parametric Form", ls="--")
circle = plt.Circle((0, 0), 1, fill=False)
ax[0].add_patch(circle)

equator_params["xp"] = 0.17
equator_para_coeffs = equator_poly_coeffs(**equator_params)
ax[1].contour(
    XS, YS, ellipse_poly_eq(**equator_2d_coeffs, x=XS, y=YS), levels=[1], colors="teal"
)
x_para, y_para = parametric_equator(**equator_params)
ax[1].plot(x_para, y_para, color="red", label="Parametric Form", ls="--")
circle = plt.Circle((0, 0), 1, fill=False)
ax[1].add_patch(circle);

[ben-cassese]

Fixed via a re-write of the poly_to_parametric_helper function, which now converts the 2D implicit coefficients to the 2D parametric coefficients through some algebra instead of solving an eigenvalue system.