Engine

squishyplanet was designed to expose OblateSystem() to the user, while keeping most of the actual computation in an “engine” directory. Modules included here are meant to rely on as few libraries as possible aside from basic jax.numpy, and are meant to contain functions that are entirely jit-able. Actually calling these functions may be helpful for building a more complex likelihood function (e.g., if you want to jointly fit multiple spectral channels, or are building a more sophisticated emission model), but in general we assume most interaction will be through the OblateSystem() class.

The code itself is based heavily on Agol, Luger, and Foreman-Mackey 2020, and the documentation makes references to specific equations in that paper. Those going through this code are encouraged to review that paper and be familiar with the Green’s basis transformation used to convert 2D surface integrals to 1D line integrals.

---

polynomial_limb_darkened_transit.ellipse_bound(c_x1: Array, c_x2: Array, c_x3: Array, c_y1: Array, c_y2: Array, c_y3: Array) → Array

Decide whether a step straddles the star edge (and so needs the quartic solve).

The outline center is (c_x3, c_y3) and a safe upper bound on the distance of any outline point from the center is bound = sqrt(c_x1^2 + c_x2^2 + c_y1^2 + c_y2^2) (it equals sqrt(r1^2 + r2^2) >= max(r1, r2)). With center distance d = sqrt(c_x3^2 + c_y3^2): d + bound < 1 is provably fully inside, d - bound > 1 provably fully outside; otherwise the step straddles and the quartic must be solved. A conservative bound is safe – a straddling-classified step with no real root is harmless (the center-inside test resolves it).

Parameters:

• c_x1 (float) – Parametric outline coefficient (x, cos term).

• c_x2 (float) – Parametric outline coefficient (x, sin term).

• c_x3 (float) – Parametric outline coefficient (x, center).

• c_y1 (float) – Parametric outline coefficient (y, cos term).

• c_y2 (float) – Parametric outline coefficient (y, sin term).

• c_y3 (float) – Parametric outline coefficient (y, center).

Returns:

True if straddling (solve the quartic), False if provably fully inside or fully outside.

Return type:

bool

polynomial_limb_darkened_transit.ellipse_star_term(para: dict, two: dict, g_coeffs: Array) → Array

The Green’s-theorem boundary integral of the blocked flux over one ellipse.

Computes \(F(E \cap S)\), the flux blocked by the intersection of a single opaque ellipse \(E\) (a planet outline or a ring edge) with the stellar disk \(S\), dotted with the Green’s-basis limb-darkening coefficients. The caller is responsible for multiplying by the overall normalization constant.

Rather than assuming a fixed number of ellipse-star intersections, this takes an “exploratory” approach that handles 0 through 4 crossings (and containment in either direction) with one code path: find all crossings, split both the ellipse outline and the stellar limb into arcs at those crossings, keep each arc if its midpoint lies inside the other region (its inside/outside status is constant along the arc, since arcs are split at every crossing), and sum the kept arcs’ contributions. Because the boundary of the intersection of convex regions is exactly the set of arcs of each curve that lie inside the other, and both curves are traversed counterclockwise, Green’s theorem makes the summation order irrelevant – no case enumeration or boundary stitching is needed.

Parameters:

• para (dict) – The parametric coefficients c_x1 ... c_y3 of the ellipse outline for a single timestep. Must trace the outline counterclockwise (true of everything produced by poly_to_parametric(), outline_prelude(), and rings.ring_para_coeffs()).

• two (dict) – The implicit-conic coefficients rho_xx ... rho_00 of the same ellipse.

• g_coeffs (Array) – The system-specific limb darkening coefficients in the Green’s basis.

Returns:

Scalar: the blocked flux from this ellipse, before normalization.

Return type:

Array

polynomial_limb_darkened_transit.lightcurve(state: dict, parameterize_with_projected_ellipse: bool) → Array

The main function for computing a transit light curve.

This function will return a 1-D array representing the flux received from the star, where each entry corresponds to a time in the input state dictionary. It first transforms the state into the implicit 3D surface of the planet, the implicit 2D sky-projected outline of the planet, and a parametric form of that outline for each time step. These are vectorized operations that are computed simultaneously across all times. It then solves for the intersection points of the planet and star, and if the planet is either partially or fully transiting, numerically solves the required 1D integrals that leverage Green’s Theorem to compute the blocked flux. The flux-blocking calculations are done sequentially for each timestep using jax.lax.scan, which seemed to be more efficient than vectorizing again while switching between branches with something like jax.lax.cond. Keep these different behaviors in mind when computing dense lightcurves with ~100s of thousands of time steps: the first part will require enough memory to compute and store ~30 values for each step, but then the actual 1D integrals will be computed sequentially.

Parameters:

• state (dict) – A dictionary containing all of the keys that are included in an OblateSystem() state attribute.

• parameterize_with_projected_ellipse (bool) – If True, the planet’s outline will be parameterized by the projected ellipse as seen by the observer. If False, the planet’s outline will be set by the full 3D parameterization of the planet. When dealing with planets that are not tidally locked and/or far from their host star and/or very close to spherical, you won’t be able to tell the difference between these two parameterizations since the projected area won’t be changing. In that case, it’s better to use the simpler 2D parameterization to avoid the degeneracies and extra computation that can arise from the 3D parameterization. This argument is static for the JIT-compiled function.

Returns:

The flux received from the star at each time step for the times included as state["times"].

Return type:

Array

polynomial_limb_darkened_transit.outline_prelude(state: dict) → tuple[dict, dict]

Direct construction of the planet’s projected outline from orbital elements.

Drop-in replacement for the planet_3d_coeffs -> planet_2d_coeffs -> poly_to_parametric chain (the full 3D-parameterization branch of lightcurve()). Instead of expanding the implicit polynomials, it builds the central 3D quadric M = R diag(d) R^T directly, eliminates z via a 2x2 Schur complement to get the projected conic, and uses the planet’s sky position skypos(...)[:2] as the outline center exactly. This avoids the near-cancellation divisions of the old chain and is machine-precision accurate (vs the chain’s ~1e-9 worst case) as well as slightly faster.

Expects state["f"] to already hold the per-timestep true anomalies. Applies the tidally-locked precession override internally (prec = where(tidally_locked, f, prec)).

Parameters:

state (dict) – An OblateSystem() state dictionary. Uses the keys a, e, f, Omega, i, omega, r, obliq, prec, f1, f2, tidally_locked.

Returns:

A tuple of two dictionaries. The first contains the implicit-conic coefficients rho_xx, rho_xy, rho_x0, rho_yy, rho_y0, rho_00 of the outline (matching planet_2d.planet_2d_coeffs()). The second contains the parametric coefficients c_x1, c_x2, c_x3, c_y1, c_y2, c_y3 with x = c_x1 cos(alpha) + c_x2 sin(alpha) + c_x3 (and similarly for y), describing the same curve as parametric_ellipse.poly_to_parametric() with a different outline-angle convention. Every leaf is broadcast to the time axis.

Return type:

tuple

polynomial_limb_darkened_transit.parameterize_2d_helper(projected_r: Array, projected_f: Array, projected_theta: Array, xc: Array, yc: Array) → tuple[dict, dict]

Convert from the alternative sky-projected parameterization to the same format used by the 3D parameterization.

A good chunk of the code assumes that the planet’s center is determined by the orbital elements and that it’s outline is derived from an equatorial radius r, a z-flattening f1, a y-flattening f2, and two body-centered rotations obliq and prec. This are useful to have when working with phase curves that are sensitive to the actual 3D orientation of the planet, but when dealing with transits only, this parameterization is overkill and allows a bunch of degeneracies. So, if only doing transits, it is more convenient to parameterize the planet by its projected radius in the x and y directions, and the angle of the projected ellipse. This function takes in those parameters and returns the same dictionaries you’d get if you fed a full 3D parameterization into planet_2d.planet_2d_coeffs().

Parameters:

• projected_r (float) – The projected “x” radius of the planet.

• projected_f (float) – The flattening of the projected ellipse.

• projected_theta (float) – The angle of the projected ellipse.

Returns:

A tuple of two dictionaries. The first dictionary contains the coefficients of the quadratic equation that describes the projected ellipse. The second dictionary contains coefficients that describe the parametric form of that same ellipse.

Return type:

tuple

polynomial_limb_darkened_transit.planet_solution_vec(a: Array, b: Array, g_coeffs: Array, c_x1: Array, c_x2: Array, c_x3: Array, c_y1: Array, c_y2: Array, c_y3: Array) → Array

Compute the “solution vector” for a 1D path across the star that lies on the outline of the planet.

This computes Eq. 21 of Agol, Luger, and Foreman-Mackey 2020. But, instead of doing it analytically, this numerically solves the required integrals with a fixed-grid Gauss-Legendre rule (see N_GL and the module-level notes). For terms s_2 and higher, this is straightforward to do based on the equations in the paper: we simply parameterize the outline of the planet by some angle \(\alpha\), then numerically integrate the dot product of Eq. 62 with that parameterization between the two endpoints of the path. For the first two lower-order terms however, Agol et al. do not provide an equivalent of Eq. 62 and instead provide only the analytic solutions. We therefore use the following as the equivalents for Eq. 62 for these terms:

\[G_0 = \{0, x\}\]

\[G_1 = \left\{0, \frac{1}{2} \left(x \sqrt{-x^2-y^2+1}-\left(y^2-1\right) \tan ^{-1}\left(\frac{x}{\sqrt{-x^2-y^2+1}}\right)\right)+\frac{\pi }{12} \right\}\]

These expressions were derived by solving the required PDE in Eq. 14 with the boundary conditions from Eq. 27. Finally, the C coefficients here describe the parametric form of the planet’s outline as seen by the observer, and they satisfy:

\[ \begin{align}\begin{aligned}x = c_{x1} \cos(\alpha) + c_{x2} \sin(\alpha) + c_{x3}\\y = c_{y1} \cos(\alpha) + c_{y2} \sin(\alpha) + c_{y3}\end{aligned}\end{align} \]

for some angle \(\alpha \in [0, 2\pi)\).

Parameters:

• a (float) – The starting parameter for the path along the planet’s outline, \(\alpha_0\).

• b (float) – The ending parameter for the path along the planet’s outline, \(\alpha_1\).

• g_coeffs (Array) – The system-specific limb darkening coefficients in the Green’s basis. Computed by multiplying the u coefficients with the change of basis matrix from greens_basis_transform.generate_change_of_basis_matrix().

• c_x1 (float) – The first coefficient of the parametric 2D outline of the planet.

• c_x2 (float) – The second coefficient of the parametric 2D outline of the planet.

• c_x3 (float) – The third coefficient of the parametric 2D outline of the planet.

• c_y1 (float) – The fourth coefficient of the parametric 2D outline of the planet.

• c_y2 (float) – The fifth coefficient of the parametric 2D outline of the planet.

• c_y3 (float) – The sixth coefficient of the parametric 2D outline of the planet.

Returns:

The solution vector for the path along the planet’s outline. The shape will match that of the input g_coeffs.

Return type:

Array

polynomial_limb_darkened_transit.star_arc_solution_vec(theta_lo: Array, theta_hi: Array, g_coeffs: Array) → Array

Compute the “solution vector” for a 1D path that lies on the edge of the star.

This is equivalent to planet_solution_vec(), but instead of integrating over paths that lie on the planet’s outline, we integrate over paths that lie on the edge of the star. As pointed out in the paragraph following Eq. 69 in Agol, Luger, and Foreman-Mackey 2020, the contribution of all terms higher than \(G_1\) will be zero in this case since we have limited ourselves to \(z=0\) by remaining on the star’s boundary. This simplifies things somewhat: the dot product of the parametric form of the star’s outline with the \(G_0\) and \(G_1\) terms written out in planet_solution_vec() reduces to elementary trigonometric polynomials, so we evaluate the required integrals in closed form via their exact antiderivatives rather than numerically.

Unlike its predecessor star_solution_vec, this function takes the polar angles of the arc endpoints on the star directly (not the parametric angles of another curve’s outline), and it makes no assumption about which of the arcs bounded by two intersection points is wanted: the caller chooses by passing the endpoints in counterclockwise order. This matters once more than one curve can cut the stellar limb into several arcs.

Parameters:

• theta_lo (float) – The starting polar angle of the arc on the star, in [0, 2 pi].

• theta_hi (float) – The ending polar angle of the arc on the star. Must satisfy theta_lo <= theta_hi <= 2 pi; the arc is traversed counterclockwise from theta_lo to theta_hi. Arcs crossing the 0 = 2 pi seam should be passed as two separate sub-arcs.

• g_coeffs (Array) – The system-specific limb darkening coefficients in the Green’s basis. Computed by multiplying the u coefficients with the change of basis matrix from greens_basis_transform.generate_change_of_basis_matrix().

Returns:

The solution vector for the path along the star’s edge. The shape will match that of the input g_coeffs.

Return type:

Array

planet_3d.planet_3d_coeffs(a: Array, e: Array, f: Array, Omega: Array, i: Array, omega: Array, r: Array, obliq: Array, prec: Array, f1: Array, f2: Array, **kwargs: object) → dict[str, Array]

Calculate and return the coefficients that describe the planet as an implicit surface in 3D space as a function of its orbital state.

This function computes a dictionary of coefficients related to the 3D position and orientation of a planet given its orbital and rotational characteristics. All inputs must be jnp.ndarrays. They can either be shape (1,) or (N,) with 1 unique N allowed per call (e.g., everything but f and prec are single-valued, but those two are length N).

Parameters:

• a (Array [Rstar]) – Semi-major axis of the orbit in units of R_star

• e (Array [Dimensionless]) – Eccentricity of the orbit

• f (Array [Radian]) – True anomaly, the angle between the direction of periapsis and the current position of the planet as seen from the star.

• Omega (Array [Radian]) – Longitude of the ascending node

• i (Array [Radian]) – Orbital inclination

• omega (Array [Radian]) – Argument of periapsis

• r (Array [Rstar]) – Equatorial radius of the planet, in units of R_star

• obliq (Array [Radian]) – Obliquity, the angle between the planet’s orbital plane and its equatorial plane. Defined when the planet is at periapsis with an Omega of zero as a rotation around the sky-frame y-axis, such that a positive obliquity tips the planet’s north pole away from the star.

• prec (Array [Radian]) – Precession angle, or equivalently the longitude of ascending node of the planet’s equatorial plane. This is defined at periapsis with an Omega of zero as a rotation about the sky-frame z-axis.

• f1 (Array [Dimensionless]) – The flattening coefficient of the planet that describes the compression along the planet’s polar axis. A value of 0.0 indicates no flattening.

• f2 (Array [Dimensionless]) – The flattening coefficient of the planet that describes the compression along the planet’s “y” axis, the vector in its equatorial plane that is perpendicular to the direction of motion at periapsis, assuming the 0.0 obliquity.

• **kwargs – Unused additional keyword arguments. These are included so that we can can take in a larger state dictionary that includes all of the required parameters along with other unnecessary ones.

Returns:

A dictionary with keys representing different coefficient names and their corresponding values. The coefficients satisfy the implicit equation:

\[p_{xx} x^2 + p_{xy} xy + p_{xz} xz + p_{x0} x + p_{yy} y^2 + p_{yz} yz + p_{y0} y + p_{zz} z^2 + p_{z0} z + p_{00} = 1\]

Return type:

dict

planet_2d.planet_2d_coeffs(p_xx: Array, p_xy: Array, p_xz: Array, p_x0: Array, p_yy: Array, p_yz: Array, p_y0: Array, p_zz: Array, p_z0: Array, p_00: Array, **kwargs: object) → dict[str, Array]

Compute the coefficients that describe the planet as an implicit 2D surface from the observer’s perspective.

This function transforms the coefficients that describe the planet’s 3D shape into a new set of coefficients that describe its projected outline as seen from z=infinity. The input p coefficients satisfy the equation:

\[p_{xx} x^2 + p_{xy} xy + p_{xz} xz + p_{x0} x + p_{yy} y^2 + p_{yz} yz + p_{y0} y + p_{zz} z^2 + p_{z0} z + p_{00} = 1\]

Parameters:

• p_xx (Array) – Coefficient representing xx interaction in the 3D model.

• p_xy (Array) – Coefficient representing xy interaction in the 3D model.

• p_xz (Array) – Coefficient representing xz interaction in the 3D model.

• p_x0 (Array) – Coefficient representing x0 interaction in the 3D model.

• p_yy (Array) – Coefficient representing yy interaction in the 3D model.

• p_yz (Array) – Coefficient representing yz interaction in the 3D model.

• p_y0 (Array) – Coefficient representing y0 interaction in the 3D model.

• p_zz (Array) – Coefficient representing zz interaction in the 3D model.

• p_z0 (Array) – Coefficient representing z0 interaction in the 3D model.

• p_00 (Array) – Coefficient representing 00 interaction in the 3D model.

• **kwargs – Unused additional keyword arguments. These are included to allow for flexibility in providing additional data that may be ignored during the computation but included for interface consistency.

Returns:

A dictionary with keys representing different transformed coefficient names (‘rho_xx’, ‘rho_xy’, ‘rho_x0’, ‘rho_yy’, ‘rho_y0’, ‘rho_00’) and their corresponding values. These coefficients describe the outline of the planet as an implicit curve that satisfies the equation:

\[\rho_{xx} x^2 + \rho_{xy} xy + \rho_{x0} x + \rho_{yy} y^2 + \rho_{y0} y + \rho_{00} = 1\]

Return type:

dict

illustrate_utils.generate_sample_radii_thetas(key: Array, num_points: Array) → tuple[Array, Array]

Create a random set of radii and thetas for sampling the planet’s surface.

These are uniformly distributed through a unit circle and will be scaled and rotated to match the planet’s shape and orientation at each timestep. However, they will be reused at every time step, which could introduce a bias but makes things much faster. Be sure you use sufficient samples to keep the bias small, then try multiple random keys to quantify it.

Parameters:

• key (Array) – A jax.random.PRNGKey for generating random numbers.

• num_points (int) – The number of points to generate.

Returns:

A tuple of two arrays, the first containing the radii and the second containing the thetas.

Return type:

Tuple

illustrate_utils.lambertian_reflection(surface_star_cos_angle: Array, x: Array, y: Array, z: Array) → Array

Compute the reflected intensity at a specific point on the planet’s surface assuming a simple Lambertian reflection model.

This is a simple model that assumes the planet reflects light according to Lambert’s cosine law, which states that the intensity of reflected light is proportional to the cosine of the angle between the surface normal and the illumination direction. That arrangement means it does not depend on the observer’s viewing angle, only the illumination angle. This helper function also assumes a uniform albedo of 1 across the planet’s surface (the final reflected flux will be scaled by the provided albedo, though is still always assumed to be uniform).

This function will also mask out any points on the planet’s surface that are on the wrong side of the terminator.

Parameters:

• surface_star_cos_angle (Array) – The cosine of the angle between the planet’s surface normal and the vector pointing from the planet’s center to the star.

• x (Array) – The x values of the points on the planet’s surface.

• y (Array) – The y values of the points on the planet’s surface.

• z (Array) – The z values of the points on the planet’s surface.

Returns:

The intensity of the reflected light at each point.

Return type:

Array

illustrate_utils.planet_surface_normal(x: Array, y: Array, z: Array, p_xx: Array, p_xy: Array, p_xz: Array, p_x0: Array, p_yy: Array, p_yz: Array, p_y0: Array, p_zz: Array, p_z0: Array, p_00: Array) → Array

Compute the unit normal vector to the planet’s surface at a given point.

The input \((x, y, z)\) points are assumed to lie on the planet’s surface.

Parameters:

• x (Array) – The x values of the points.

• y (Array) – The y values of the points.

• z (Array) – The z values of the points.

• p_xx (Array) – xx coefficient in the 3D implicit representation.

• p_xy (Array) – xy coefficient in the 3D implicit representation.

• p_xz (Array) – xz coefficient in the 3D implicit representation.

• p_x0 (Array) – x0 coefficient in the 3D implicit representation.

• p_yy (Array) – yy coefficient in the 3D implicit representation.

• p_yz (Array) – yz coefficient in the 3D implicit representation.

• p_y0 (Array) – y0 coefficient in the 3D implicit representation.

• p_zz (Array) – zz coefficient in the 3D implicit representation.

• p_z0 (Array) – z0 coefficient in the 3D implicit representation.

• p_00 (Array) – 00 coefficient in the 3D implicit representation.

Returns:

An array of shape (3, n) containing the unit normal vectors at each point.

Return type:

Array

illustrate_utils.sample_surface(sample_radii: Array, sample_thetas: Array, rho_xx: Array, rho_xy: Array, rho_x0: Array, rho_yy: Array, rho_y0: Array, rho_00: Array, p_xx: Array, p_xy: Array, p_xz: Array, p_x0: Array, p_yy: Array, p_yz: Array, p_y0: Array, p_zz: Array, p_z0: Array, p_00: Array, **kwargs: object) → tuple[Array, Array, Array]

Convert randomly sampled \((x, y)\) points on the projected planet to \((x, y, z)\) points on the planet’s surface.

The \(rho\) coefficients are calculated with planet_2d.planet_2d_coeffs(), the \(p\) coefficients are calculated with planet_3d.planet_3d_coeffs(), and the sample points are generated with generate_sample_radii_thetas().

Parameters:

• sample_radii (Array) – The radii of the sampled points.

• sample_thetas (Array) – The angles of the sampled points.

• rho_xx (Array) – xx coefficient in the 2D implicit representation.

• rho_xy (Array) – xy coefficient in the 2D implicit representation.

• rho_x0 (Array) – x0 coefficient in the 2D implicit representation.

• rho_yy (Array) – yy coefficient in the 2D implicit representation.

• rho_y0 (Array) – y0 coefficient in the 2D implicit representation.

• rho_00 (Array) – 00 coefficient in the 2D implicit representation.

• p_xx (Array) – xx coefficient in the 3D implicit representation.

• p_xy (Array) – xy coefficient in the 3D implicit representation.

• p_xz (Array) – xz coefficient in the 3D implicit representation.

• p_x0 (Array) – x0 coefficient in the 3D implicit representation.

• p_yy (Array) – yy coefficient in the 3D implicit representation.

• p_yz (Array) – yz coefficient in the 3D implicit representation.

• p_y0 (Array) – y0 coefficient in the 3D implicit representation.

• p_zz (Array) – zz coefficient in the 3D implicit representation.

• p_z0 (Array) – z0 coefficient in the 3D implicit representation.

• p_00 (Array) – 00 coefficient in the 3D implicit representation.

• **kwargs – Additional unused keyword arguments, included so that we can pass in a larger state dictionary that includes all of the required parameters along with other unnecessary ones.

Returns:

A tuple of three arrays, the first containing the x values, the second containing the y values, and the third containing the z values.

Return type:

Tuple

illustrate_utils.surface_star_cos_angle(planet_surface_normal: Array, x_c: Array, y_c: Array, z_c: Array, **kwargs: object) → Array

A helper function to compute the cosine of the angle between the planet’s surface normal vector and the vector linking the center of the planet to the star.

This is an approximation that the star is a) a point source and b) that all light coming from the star is parallel. Neither of these are strictly true. The former could be handled the same way starry does it, by distributing point sources across the surface of the star and averaging. I don’t know of any attempts to address the latter, though in principle it wouldn’t be hard to do here since we’re already doing so much numerically.

Parameters:

• planet_surface_normal (Array) – The unit normal vectors to the planet’s surface.

• x_c (Array) – The x coordinate of the center of the planet.

• y_c (Array) – The y coordinate of the center of the planet.

• z_c (Array) – The z coordinate of the center of the planet.

• **kwargs – Additional unused keyword arguments, included so that we can pass in a larger state dictionary that includes all of the required parameters along with other unnecessary ones.

Returns:

The cosine of the angle between the planet’s surface normal and the vector pointing from the planet’s center to the star.

Return type:

Array

parametric_ellipse.cartesian_intersection_to_parametric_angle(xs: Array, ys: Array, c_x1: Array, c_x2: Array, c_x3: Array, c_y1: Array, c_y2: Array, c_y3: Array) → Array

Given a set of x and y coordinates corresponding to the intersection of the planet and star, compute the angle \(\alpha\) that corresponds to each point.

Here, \(\alpha\) is the parameter in the parametric equations of the ellipse. See poly_to_parametric() for more details.

Parameters:

• xs (Array [Rstar]) – x-coordinates of the intersection points

• ys (Array [Rstar]) – y-coordinates of the intersection points

• c_x1 (Array [Dimensionless]) – Coefficient of x^2

• c_x2 (Array [Dimensionless]) – Coefficient of xy

• c_x3 (Array [Dimensionless]) – Coefficient of x

• c_y1 (Array [Dimensionless]) – Coefficient of y^2

• c_y2 (Array [Dimensionless]) – Coefficient of y

• c_y3 (Array [Dimensionless]) – Constant term

Returns:

The angle \(\alpha\) corresponding to each intersection point

Return type:

Array [Rstar]

parametric_ellipse.ellipse_angle(rho_xx: Array, rho_xy: Array, rho_yy: Array) → Array

Rotation angle, in [0, pi), of the centered conic rho_xx x^2 + rho_xy x y + rho_yy y^2 = 1.

The only genuinely degenerate input is a circle (rho_xx == rho_yy and rho_xy == 0), where any angle is valid and arctan2(0, 0) has a NaN gradient. A tilted ellipse with rho_xx == rho_yy (e.g. a major axis at exactly 45 degrees) is not degenerate: arctan2(rho_xy, 0) = +-pi/2 is well defined and must not be zeroed.

parametric_ellipse.parametric_to_poly(c_x1: Array, c_x2: Array, c_x3: Array, c_y1: Array, c_y2: Array, c_y3: Array) → dict[str, Array]

Convert the coefficients of a parametric ellipse to those of an implicit one.

This is the inverse of poly_to_parametric(). The input coefficients describe the ellipse as a parametric curve for parameter \(\alpha\):

\[x = c_{x1} * \cos(\alpha) + c_{x2} * \sin(\alpha) + c_{x3} y = c_{y1} * \cos(\alpha) + c_{y2} * \sin(\alpha) + c_{y3}\]

Note that the returned coefficients involve division by the squared determinant of the parametric coefficient matrix, so they blow up for nearly-degenerate (nearly zero-area) ellipses. Prefer keeping such curves in parametric form.

Parameters:

• c_x1 (Array [Rstar]) – Coefficient in the parametric equation

• c_x2 (Array [Rstar]) – Coefficient in the parametric equation

• c_x3 (Array [Rstar]) – Coefficient in the parametric equation

• c_y1 (Array [Rstar]) – Coefficient in the parametric equation

• c_y2 (Array [Rstar]) – Coefficient in the parametric equation

• c_y3 (Array [Rstar]) – Coefficient in the parametric equation

Returns:

Dictionary of coefficients for the implicit ellipse equation:

\[\rho_{xx} x^2 + \rho_{xy} xy + \rho_{x0} x + \rho_{yy} y^2 + \rho_{y0} y + \rho_{00} = 1\]

Return type:

dict

parametric_ellipse.point_in_ellipse(x: Array, y: Array, c_x1: Array, c_x2: Array, c_x3: Array, c_y1: Array, c_y2: Array, c_y3: Array) → Array

Test whether a point is strictly inside a parametrically-defined ellipse.

Solves the 2x2 linear system for \((\cos\alpha, \sin\alpha)\) via the adjugate, so the test involves no division and remains well-conditioned even for very thin ellipses (unlike evaluating the implicit form, whose coefficients scale as the inverse squared semi-minor axis). A degenerate (zero-area) ellipse has an empty interior and always returns False.

Parameters:

• x (Array [Rstar]) – x-coordinate(s) of the point(s) to test

• y (Array [Rstar]) – y-coordinate(s) of the point(s) to test

• c_x1 (Array [Rstar]) – Coefficient in the parametric equation

• c_x2 (Array [Rstar]) – Coefficient in the parametric equation

• c_x3 (Array [Rstar]) – Coefficient in the parametric equation

• c_y1 (Array [Rstar]) – Coefficient in the parametric equation

• c_y2 (Array [Rstar]) – Coefficient in the parametric equation

• c_y3 (Array [Rstar]) – Coefficient in the parametric equation

Returns:

Boolean array, True where the point is inside the ellipse.

Return type:

Array [Dimensionless]

parametric_ellipse.poly_to_parametric(rho_xx: Array, rho_xy: Array, rho_x0: Array, rho_yy: Array, rho_y0: Array, rho_00: Array) → dict[str, Array]

Convert between the coefficients that describe an implicit to those defining a parametric ellipse.

The input coefficients are those of the implicit ellipse equation:

\[\rho_{xx} x^2 + \rho_{xy} xy + \rho_{x0} x + \rho_{yy} y^2 + \rho_{y0} y + \rho_{00} = 1\]

Parameters:

• rho_xx (Array [Dimensionless]) – Coefficient of x^2

• rho_xy (Array [Dimensionless]) – Coefficient of xy

• rho_x0 (Array [Dimensionless]) – Coefficient of x

• rho_yy (Array [Dimensionless]) – Coefficient of y^2

• rho_y0 (Array [Dimensionless]) – Coefficient of y

• rho_00 (Array [Dimensionless]) – Constant term

Returns:

Dictionary of coefficients for the parametric ellipse. The ellipse can now be described by the following parametric equations for parameter \(\alpha\):

\[x = c_{x1} * \cos(\alpha) + c_{x2} * \sin(\alpha) + c_{x3} y = c_{y1} * \cos(\alpha) + c_{y2} * \sin(\alpha) + c_{y3}\]

Return type:

dict

parametric_ellipse.poly_to_parametric_helper(rho_xx: Array, rho_xy: Array, rho_x0: Array, rho_yy: Array, rho_y0: Array, rho_00: Array) → tuple[Array, Array, Array, Array, Array, Array]

A helper function for poly_to_parametric().

Parameters:

• rho_xx (Array [Dimensionless]) – Coefficient of x^2

• rho_xy (Array [Dimensionless]) – Coefficient of xy

• rho_x0 (Array [Dimensionless]) – Coefficient of x

• rho_yy (Array [Dimensionless]) – Coefficient of y^2

• rho_y0 (Array [Dimensionless]) – Coefficient of y

• rho_00 (Array [Dimensionless]) – Constant term

Returns:

• r1 (Array [Rstar]): Semi-major axis of the projected ellipse

• r2 (Array [Rstar]): Semi-minor axis of the projected ellipse

• xc (Array [Rstar]): x-coordinate of the center of the ellipse

• yc (Array [Rstar]): y-coordinate of the center of the ellipse

• cosa (Array [Dimensionless]): Cosine of the rotation angle

• sina (Array [Dimensionless]): Sine of the rotation angle

Return type:

Tuple

rings.para_eval(alpha: Array, coeffs: dict[str, Array]) → tuple

Evaluate a parametric ellipse at the given angle(s).

Parameters:

• alpha (Array [Radian]) – The parametric angle(s) at which to evaluate the curve.

• coeffs (dict) – Dictionary of the parametric coefficients c_x1 through c_y3.

Returns:

The x and y coordinates of the curve at each alpha.

Return type:

Tuple

rings.parametric_conic_intersections(c_x1: Array, c_x2: Array, c_x3: Array, c_y1: Array, c_y2: Array, c_y3: Array, rho_xx: Array, rho_xy: Array, rho_x0: Array, rho_yy: Array, rho_y0: Array, rho_00: Array) → tuple[Array, Array, Array]

Intersections between a parametric ellipse and an implicitly-defined conic.

This is a more general version of _single_intersection_points() in polynomial_limb_darkened_transit, which assumes the implicit curve is the unit circle (the star). Substituting the parametric curve into the implicit equation via the tangent half-angle substitution \(t = \tan(\alpha/2)\) yields a quartic in t.

The parametric curve should be the one at risk of degeneracy (e.g. a nearly edge-on ring): its coefficients stay bounded, whereas the implicit form of a thin ellipse has coefficients that diverge as the inverse squared semi-minor axis.

The implicit coefficients satisfy the following on the conic:

\[\rho_{xx} x^2 + \rho_{xy} xy + \rho_{x0} x + \rho_{yy} y^2 + \rho_{y0} y + \rho_{00} = 1\]

And the c coefficients describe the parametric curve:

\[x = c_{x1} \cos(\alpha) + c_{x2} \sin(\alpha) + c_{x3} y = c_{y1} \cos(\alpha) + c_{y2} \sin(\alpha) + c_{y3}\]

Parameters:

• c_x1 (Array [Rstar]) – Coefficient in the parametric form of the first curve.

• c_x2 (Array [Rstar]) – Coefficient in the parametric form of the first curve.

• c_x3 (Array [Rstar]) – Coefficient in the parametric form of the first curve.

• c_y1 (Array [Rstar]) – Coefficient in the parametric form of the first curve.

• c_y2 (Array [Rstar]) – Coefficient in the parametric form of the first curve.

• c_y3 (Array [Rstar]) – Coefficient in the parametric form of the first curve.

• rho_xx (Array [Dimensionless]) – Coefficient in the implicit form of the second curve.

• rho_xy (Array [Dimensionless]) – Coefficient in the implicit form of the second curve.

• rho_x0 (Array [Dimensionless]) – Coefficient in the implicit form of the second curve.

• rho_yy (Array [Dimensionless]) – Coefficient in the implicit form of the second curve.

• rho_y0 (Array [Dimensionless]) – Coefficient in the implicit form of the second curve.

• rho_00 (Array [Dimensionless]) – Coefficient in the implicit form of the second curve.

Returns:

Three arrays of length 4: the parametric angles alpha of the intersection points on the first curve (in (-pi, pi]), and the x and y coordinates of those points. Slots corresponding to imaginary quartic roots (no intersection) are filled with 999.

Return type:

Tuple

rings.ring_para_coeffs(a: Array, e: Array, f: Array, Omega: Array, i: Array, omega: Array, rRing: Array, ring_obliq: Array, ring_prec: Array, **kwargs: Array) → dict[str, Array]

Compute the coefficients describing the parametric form of a ring in the sky plane.

Remember that ring_obliq and ring_prec are defined in the planet’s orbital plane, at f=0. That implies that to get a face-on ring, you need ring_obliq=ring_prec=90 deg. Some other examples all at inc=90 deg:

• ring_obliq=0, ring_prec=0: the ring is an edge-on horizontal line

• ring_obliq=jnp.pi/4, ring_prec=0: the ring is still a line, but now it’s tilted 45 degrees in the sky frame. You’ve tipped the north pole away from the star.

• ring_obliq=0, ring_prec=anything: the ring is still a line

• ring_obliq=90, ring_prec=0: the ring is a face-on circle

Making inc != 90 deg will alter these: the angles are defined in the orbital plane, so if you tilt the orbit away from face-on, you’ll also tilt the ring away from face-on.

The returned coefficients are always oriented so that the curve is traversed counterclockwise on the sky as the parametric angle increases, as required for the Green’s theorem boundary integrals.

Parameters:

• a (Array [Rstar]) – The semi-major axis of the planet.

• e (Array [Dimensionless]) – The eccentricity of the planet.

• f (Array [Radian]) – The true anomaly of the planet.

• Omega (Array [Radian]) – The longitude of the ascending node of the planet.

• i (Array [Radian]) – The inclination of the planet.

• omega (Array [Radian]) – The argument of periapsis of the planet.

• rRing (Array [Rstar]) – The radius of the ring.

• ring_obliq (Array [Radian]) – The obliquity of the ring.

• ring_prec (Array [Radian]) – The precession angle of the ring.

• **kwargs – Additional (unused) keyword arguments.

Returns:

A dictionary with keys for each of the coefficients of the parametric form of the ring.

Return type:

dict

rings.ring_prelude(state: dict) → tuple[dict, dict]

Per-timestep parametric coefficients for both ring edges.

Expects state["f"] to already hold the per-timestep true anomalies (i.e. run polynomial_limb_darkened_transit._lightcurve_setup() first). When the ring orientation tracks the planet’s (ring_tracks_planet), a tidally locked planet drags the ring precession angle along with the true anomaly, mirroring the prec override in polynomial_limb_darkened_transit.outline_prelude().

Parameters:

state (dict) – A RingedSystem() state dictionary. Uses the orbital elements plus ring_inner_r, ring_outer_r, ring_obliq, ring_prec, ring_tracks_planet, and tidally_locked.

Returns:

Two dictionaries of parametric coefficients (each broadcast to the time axis): the outer ring edge and the inner ring edge.

Return type:

Tuple

Green’s-theorem machinery for transits of ringed planets.

The blocked flux of a planet P with an opaque annular ring (inner/outer edge ellipses I and O, with I inside O since the projected edges are concentric similar ellipses) in front of the star S decomposes exactly, for any surface brightness, as

\[F(B) = F(P \cap S) + F(O \cap S) - F(I \cap S) - F(P \cap O \cap S) + F(P \cap I \cap S)\]

Every term is an intersection of 2-3 convex regions, hence connected, and its boundary is exactly the set of arcs of each bounding curve that lie inside all the other regions. Each term is therefore computed with the same “exploratory” recipe as polynomial_limb_darkened_transit.ellipse_star_term(): find all pairwise crossings, split every curve into arcs at its crossings, keep the arcs whose midpoints lie inside the other regions, and sum their Green’s-theorem contributions (all curves are traversed counterclockwise, so the summation order is irrelevant).

Ring edges can be arbitrarily close to degenerate (edge-on), where their implicit conic coefficients diverge as the inverse squared semi-minor axis. All ring intersections are therefore solved with the ring on the parametric side (rings.parametric_conic_intersections() against the planet’s or star’s well-conditioned implicit form), and all membership tests use the division-free parametric_ellipse.point_in_ellipse().

ringed_transit.ring_star_term(ring_para: dict, rs_alphas: Array, rs_xs: Array, rs_ys: Array, g_coeffs: Array) → Array

The blocked-flux boundary integral F(E intersect S) for one ring edge.

The ring-edge analog of polynomial_limb_darkened_transit.ellipse_star_term(), except the ring-star crossings are precomputed by the caller (with rings.parametric_conic_intersections() against STAR_TWO, keeping the possibly-degenerate ring on the parametric side) and shared with the triple term.

Parameters:

• ring_para (dict) – Parametric coefficients of the ring edge (single timestep, counterclockwise).

• rs_alphas (Array) – Ring-side parametric angles of the ring-star crossings.

• rs_xs (Array) – Crossing x-coordinates (999 sentinels).

• rs_ys (Array) – Crossing y-coordinates.

• g_coeffs (Array) – Green’s-basis limb darkening coefficients.

Returns:

Scalar blocked flux from this ring edge, before normalization.

Return type:

Array

ringed_transit.ringed_lightcurve(state: dict) → Array

The main function for computing the transit light curve of a ringed planet.

The ringed analog of polynomial_limb_darkened_transit.lightcurve(): the same vectorized setup and sequential per-timestep jax.lax.scan, but each in-transit step evaluates the five-term inclusion-exclusion decomposition described in the module docstring instead of a single ellipse-star term. Two properties keep the common cases cheap:

• If the planet’s projected outline never pokes out of the inner ring edge (no planet-inner-edge crossings, the usual case for a mostly face-on ring), both triple terms equal F(P intersect S) exactly and cancel, so only the three two-region terms are evaluated.

• An exactly edge-on ring projects to a degenerate (zero-area) ellipse whose counterclockwise traversal doubles back on itself, making every ring term integrate to zero automatically – no special-casing is needed, and the ringed light curve collapses to the planet-only one continuously.

Requires the full 3D parameterization (parameterize_with_projected_ellipse is not meaningful for rings, which need the true orbital-frame orientation).

Parameters:

state (dict) – A dictionary containing all of the keys that are included in a RingedSystem() state attribute.

Returns:

The flux received from the star at each time step for the times included as state["times"].

Return type:

Array

ringed_transit.triple_term(planet_para: dict, ring_para: dict, ps_xs: Array, ps_ys: Array, rs_alphas: Array, rs_xs: Array, rs_ys: Array, pr_alphas: Array, pr_xs: Array, pr_ys: Array, g_coeffs: Array) → Array

The blocked-flux boundary integral F(P intersect E intersect S).

The three-region term of the ringed-planet inclusion-exclusion decomposition, for the planet P, one ring edge E, and the star S. All three pairwise crossing sets are precomputed by the caller (each pairwise quartic is solved once per timestep and shared across the decomposition’s terms). Each of the three curves is split into up to 9 arcs at its (up to 8) crossings with the other two curves, and an arc contributes exactly when its midpoint lies inside both other regions. Containment cases (e.g. P entirely inside E) fall out of the same code path via the bookended [0, 2 pi] arc, whose midpoint doubles as a whole-curve containment test.

Parameters:

• planet_para (dict) – Parametric coefficients of the planet outline (counterclockwise; single timestep).

• ring_para (dict) – Parametric coefficients of the ring edge (counterclockwise).

• ps_xs (Array) – Planet-star crossing x-coordinates (999 sentinels), from polynomial_limb_darkened_transit._single_intersection_points().

• ps_ys (Array) – Planet-star crossing y-coordinates.

• rs_alphas (Array) – Ring-side angles of the ring-star crossings, from rings.parametric_conic_intersections() against STAR_TWO.

• rs_xs (Array) – Ring-star crossing x-coordinates.

• rs_ys (Array) – Ring-star crossing y-coordinates.

• pr_alphas (Array) – Ring-side angles of the planet-ring crossings, from rings.parametric_conic_intersections() against the planet’s implicit conic.

• pr_xs (Array) – Planet-ring crossing x-coordinates.

• pr_ys (Array) – Planet-ring crossing y-coordinates.

• g_coeffs (Array) – Green’s-basis limb darkening coefficients.

Returns:

Scalar blocked flux of the triple-intersection region, before normalization.

Return type:

Array

kepler.kepler(M: Array, ecc: Array) → Array

Solve Kepler’s equation to compute the true anomaly.

This implementation is based on that within jaxoplanet, many thanks to the authors.

Parameters:

• M (Array [Radian]) – Mean anomaly

• ecc (Array [Dimensionless]) – Eccentricity

Returns:

True anomaly in radians

Return type:

Array

kepler.skypos(a: Array, e: Array, f: Array, Omega: Array, i: Array, omega: Array, **kwargs: object) → Array

Compute the cartesian coordinates of the center of the planet in the sky frame given its orbital elements.

Parameters:

• a (Array [Rstar]) – Semi-major axis of the orbit in units of R_star

• e (Array [Dimensionless]) – Eccentricity of the orbit

• f (Array [Radian]) – True anomaly, the angle between the direction of periapsis and the current position of the planet as seen from the star.

• Omega (Array [Radian]) – Longitude of the ascending node

• i (Array [Radian]) – Orbital inclination

• omega (Array [Radian]) – Argument of periapsis

• **kwargs – Unused additional keyword arguments. These are included so that we can can take in a larger state dictionary that includes all of the required parameters along with other unnecessary ones.

Returns:

The cartesian coordinates of the planet in the sky frame. Shape [3, N].

Return type:

Array

kepler.t0_to_t_peri(e: Array, i: Array, omega: Array, period: Array, t0: Array, **kwargs: object) → Array

kepler.true_anomaly_at_transit_center(e: Array, i: Array, omega: Array) → Array

Computes the true anomaly at the instant of minimum star/planet separation.

Uses equations 4.12-4.18 of Kipping 2011 to compute the true anomaly at the instant of minimum star/planet separation.

Parameters:

• e (Array [Dimensionless]) – Eccentricity of the orbit

• i (Array [Radian]) – Orbital inclination

• omega (Array [Radian]) – Argument of periapsis

Returns:

True anomaly at the instant of minimum star/planet separation in radians.

Return type:

Array

greens_basis_transform.generate_change_of_basis_matrix(N: int) → Array

Generate the change of basis matrix to convert limb darking u coefficients to Green’s basis coefficients.

This function is only run once per system, though the resulting matrix is used repeatedly in the light curve calculation. It implements Eq. 17 of Agol, Luger, and Foreman-Mackey 2020.

Parameters:

N (int) – The order of the polynomial limb darkening law.

Returns:

The change of basis matrix. When solving for the blocked flux, this will be multiplied by the u limb darkening coefficients to convert them to the Green’s basis.

Return type:

Array