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 tranformation used to convert 2D surface integrals to 1D line integrals.

---

polynomial_limb_darkened_transit.lightcurve(state, parameterize_with_projected_ellipse)

The main function for computing a transit light curve.

This function will return a 1-D array representing the flux recieved 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 simulataneously 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 braches 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.parameterize_2d_helper(projected_r, projected_f, projected_theta, xc, yc)

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, b, g_coeffs, c_x1, c_x2, c_x3, c_y1, c_y2, c_y3)

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 uses the quadax package to numerically solve the required integrals. 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_solution_vec(a, b, g_coeffs, c_x1, c_x2, c_x3, c_y1, c_y2, c_y3)

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

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, though we do still have to numerically integrate 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(). Technically we probably could use the analytic solutions for these terms, but so far we have not.

Parameters:

• a (float) – The starting parameter for the path along the star’s outline, \(\alpha_0\). Note: here \(\alpha\) is the angle parameterizing the path on the planet’s outline, not the star’s, even though the path we will integrate over is on the star. We convert to the relevant parameters internally.

• b (float) – The ending parameter for the path along the star’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

planet_3d.extended_illumination_offsets(a, e, f, Omega, i, omega, extended_illumination_points, **kwargs)

Generate a set of points uniformly distributed on the portion of the star visible from the planet.

This takes a spherical cap of points centered on the north pole of the star, squishes them together so they’d all be visible from the planet, and then rotates them to be centered on the sub-planet point on the star.

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

• extended_illumination_points (Array) – A set of points lying on a unit hemisphere centered at the origin that are evenly distributed across the projected disk when viewed from above. Created in during initialization of an OblateSystem object.

• **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 set of points on the star that are visible from the planet, evenly distributed across the projected disk of the star as seen by the planet, and centered on the sub-planet point on the star.

Return type:

Array

planet_3d.planet_3d_coeffs(a, e, f, Omega, i, omega, r, obliq, prec, f1, f2, **kwargs)

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 coeffients 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_3d.planet_3d_coeffs_extended_illumination(a, e, f, Omega, i, omega, r, obliq, prec, f1, f2, offsets, **kwargs)

Generate many sets of p coefficients that describe same planet offset from its true position by different amounts.

Since the star is not actually a point source, we slightly underestimate the area of the illuminated portion of the planet. The limb of the star can “see around the horizon”, and this extra illumination will affect the reflected portion of a phase curve. To crudely account for this, we can break the star into many point sources distributed over the portion of the star that is visible from the planet, then add their resulting lightcurves. This isn’t perfect for a few reasons: how should we distribute this point sources, and how should we weight them? Also, for a non-spherical planet, what do we mean by “the portion of the star that is visible from the planet”? For now, we avoid those questions by assigning equal intensities to a set of points distributed uniformly over the portion of the hemisphere of the star that would be visible to an observer at the center of the planet.

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.

• offsets (Array [Rstar]) – An array of offsets from the planet’s true position. Each offset is a 3-element array representing the x, y, and z offsets from the planet’s true position in units of R_star. Used when splitting the star into many point sources to account for extended illumination.

• **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 similar to the one returned by planet_3d_coeffs(), but with now describing the planet after being translated by the provided offsets.

Return type:

dict

planet_2d.planet_2d_coeffs(p_xx, p_xy, p_xz, p_x0, p_yy, p_yz, p_y0, p_zz, p_z0, p_00, **kwargs)

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

phase_curve_utils.corrected_emission_profile(x, y, z, transform, r, f1, f2, hotspot_latitude, hotspot_longitude, hotspot_concentration, **kwargs)

A helper function to emission_at_timestep(), broken out only to be used for illustrations in OblateSystem.illustrate().

phase_curve_utils.emission_at_timestep(x, y, z, transform, r, f1, f2, hotspot_latitude, hotspot_longitude, hotspot_concentration)

Compute the emitted intensity at a given point on the planet’s surface.

Parameters:

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

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

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

• transform (Array) – The rotation matrix to transform the sky frame to the planet’s frame, calculated with pre_squish_transform().

• r (Array) – The equatorial radius of the planet.

• f1 (Array) – The planet’s \(z\) flattening coefficient.

• f2 (Array) – The planet’s \(y\) flattening coefficient.

• hotspot_latitude (Array) – The “latitude” of the hotspot on the planet. Defined the physics way for \(\theta\) though, not the geography way: 0 is the north pole, \(\pi/2\) is the equator, and \(\pi\) is the south pole.

• hotspot_longitude (Array) – The longitude of the hotspot on the planet.

• hotspot_concentration (Array) – The concentration of the hotspot on the planet. \(\kappa\) in the von Mises-Fisher distribution.

Returns:

The intensity of the emitted light at each point.

Return type:

Array

phase_curve_utils.emission_phase_curve(sample_radii, sample_thetas, two, three, state, **kwargs)

Compute the timeseries of the emitted light from the planet.

This function does a Monte Carlo estimation of the visible flux emitted by the planet at each time step assuming that a) the surface intensity is modeled by a von Mises-Fisher distribution and b) the planet is a Lambertian emitter. To save on computation, it takes one set of randomly generated samples of a unit disk, then coverts these to points on the planet’s visible disk at each timestep. This introduces some bias- be sure to use a large enough sample it is below an appropriate threshold. Also, to compute secondary eclipses, samples are zeroed out when they fall behind the star.

Parameters:

• sample_radii (Array) – Randomly sampled radii from a unit sphere, used to generate points on the visible disk of the planet at each timestep. Create with generate_sample_radii_thetas().

• sample_thetas (Array) – Randomly sampled thetas from a unit sphere, used to generate points on the visible disk of the planet at each timestep. Create with generate_sample_radii_thetas().

• two (dict) – A dictionary containing the rho coefficients of the planet’s implicit 2D representation, as seen from the observer and calculated with planet_2d.planet_2d_coeffs().

• three (dict) – A dictionary containing the p coefficients of the planet’s 3D shape, as seen from the observer and calculated with planet_3d.planet_3d_coeffs().

• state (dict) – A dictionary containing all of the parameters needed to compute the phase curve. This includes the planet’s orbital parameters, the observer’s parameters, and the hotspot parameters.

Returns:

The timeseries of observed emitted light from the planet. Each element of the array is the total observed emitted flux at that corresponding time in the state[“times”] array.

Return type:

Array

phase_curve_utils.extended_illumination_reflected_phase_curve(sample_radii, sample_thetas, two, three, state, x_c, y_c, z_c, offsets)

WIP, not yet implemented. Hiding behind a NotImplementedError when setting extended_illumination_npts to anything greater than 1 when initializing an OblateSystem object.

phase_curve_utils.generate_sample_radii_thetas(key, num_points)

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 re-used 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

phase_curve_utils.lambertian_reflection(surface_star_cos_angle, x, y, z)

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

phase_curve_utils.phase_curve(sample_radii, sample_thetas, two, three, state, x_c, y_c, z_c)

Compute the reflected and emitted phase curves of the planet.

This is essentially a wrapper for reflected_phase_curve() and emission_phase_curve(), except is reuses computations where it can, and also applies the appriate scalings to each (albedo/distance from star/area seen by star for reflection, and the emitted scale for emission).

Parameters:

• sample_radii (Array) – Randomly sampled radii from a unit sphere, used to generate points on the visible disk of the planet at each timestep. Create with generate_sample_radii_thetas().

• sample_thetas (Array) – Randomly sampled thetas from a unit sphere, used to generate points on the visible disk of the planet at each timestep. Create with generate_sample_radii_thetas().

• two (dict) – A dictionary containing the rho coefficients of the planet’s implicit 2D representation, as seen from the observer and calculated with planet_2d.planet_2d_coeffs().

• three (dict) – A dictionary containing the p coefficients of the planet’s 3D shape, as seen from the observer and calculated with planet_3d.planet_3d_coeffs().

• state (dict) – A dictionary containing all of the parameters needed to compute the phase curve. This includes the planet’s orbital parameters, the observer’s parameters, and the hotspot parameters.

• 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.

Returns:

The correctly scaled reflected and emitted contributions to the phase curve.

Return type:

Tuple

phase_curve_utils.planet_from_star(p_xx, p_xy, p_xz, p_x0, p_yy, p_yz, p_y0, p_zz, p_z0, p_00, x_c, y_c, z_c, **kwargs)

Compute the coefficients of the planet’s 3D shape from the star’s perspective, as if it were aligned the the \(z\) axis.

When computing the reflected flux from the planet, we need to know how much flux initial reaches it from the star. To do that, we need to know the planet’s projected area as seen from the star, which importantly, could be different than the projected area as seen from the observer. To compute this area, we first use this function to get a 3D representation of the planet as seen from the star, then will use those coefficients to compute an implicit 2D representation, then will use those to get the area.

The x_c, y_c, and z_c inputs are all technically encoded in the p inputs as well, but it was easier just to carry them around explicitly.

Parameters:

• 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.

• 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.

Returns:

A dictionary containing the coefficients of the planet’s shape as seen from the star. Will look identical to the output of planet_3d.planet_3d_coeffs().

Return type:

dict

phase_curve_utils.planet_surface_normal(x, y, z, p_xx, p_xy, p_xz, p_x0, p_yy, p_yz, p_y0, p_zz, p_z0, p_00)

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

phase_curve_utils.pre_squish_transform(a, e, f, Omega, i, omega, r, obliq, prec, **kwargs)

Compute the rotation matrix to go from the sky frame to the planet’s.

This is the underlying transformation behind everything in planet_3d.planet_3d_coeffs(), except that module never actually uses it in this form since it applies it then goes ahead and gathers terms.

Parameters:

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

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

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

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

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

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

• r (Array) – The equatorial radius of the planet.

• obliq (Array) – The obliquity of the planet.

• prec (Array) – The precession 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:

A matrix that can be used to rotate vectors from the sky frame to the planet’s frame.

Return type:

Array

phase_curve_utils.reflected_normalization(two, three, x_c, y_c, z_c, xo=0.0, yo=0.0, zo=0.0, **kwargs)

Compute the time-dependent normalization factor for the reflected light.

The reflected light computations are almost entirely carried out assuming the star is a point source 1 R_star from the center of the planet emitting plane-parallel rays. To convert these to actual reflected flux, we need to a) correct for the distance between the planet and the star and b) account for how much area the planet actually subtends as seen from the star. a) is easy and common across all implementations, it’s just the inverse square law. b) is more complicated for oblate planets than spherical planets however, since even on circular orbits, the subtended area (and consequently area that is recieves flux and and is able to reflect it) can change as a function of orbital phase. Note however that it will not vary with phase if the planet is tidally locked and always shows the same face to the star.

Parameters:

• two (dict) – A dictionary containing the rho coefficients of the planet’s implicit 2D representation, as seen from the observer and calculated with planet_2d.planet_2d_coeffs().

• three (dict) – A dictionary containing the p coefficients of the planet’s 3D shape, as seen from the observer and calculated with planet_3d.planet_3d_coeffs().

• 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.

• xo (float or Array) – An offset to add to the x coordinate of the center of the planet, used when correcting for extended source illuminations. Default is 0.0.

• yo (float or Array) – An offset to add to the y coordinate of the center of the planet, used when correcting for extended source illuminations. Default is 0.0.

• zo (float or Array) – An offset to add to the z coordinate of the center of the planet, used when correcting for extended source illuminations. Default is 0.0.

• **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 normalization factor for the reflected light.

Return type:

Array

phase_curve_utils.reflected_phase_curve(sample_radii, sample_thetas, two, three, state, x_c, y_c, z_c, xo=Array([0.], dtype=float64), yo=Array([0.], dtype=float64), zo=Array([0.], dtype=float64))

Compute the timeseries of light reflected from the planet.

This function computes the reflected light from the planet at each time step. It assume the planet is a) a Lambertian reflector, b) that the star is a point source sending out parallel rays, c) that the \(a/R_s >> R_p\) (i.e., the distance between each point on the surface to the star is essentially constant), and d) that the planet has a spatially uniform albedo of unity (the entire curve can be scaled by an actual albedo later, but the assumption of uniformity is baked-in). However, it does take into account the planet’s oblateness and orientation.

Parameters:

• sample_radii (Array) – Randomly sampled radii from a unit sphere, used to generate points on the visible disk of the planet at each timestep. Create with generate_sample_radii_thetas().

• sample_thetas (Array) – Randomly sampled thetas from a unit sphere, used to generate points on the visible disk of the planet at each timestep. Create with generate_sample_radii_thetas().

• two (dict) – A dictionary containing the rho coefficients of the planet’s implicit 2D representation, as seen from the observer and calculated with planet_2d.planet_2d_coeffs().

• three (dict) – A dictionary containing the p coefficients of the planet’s 3D shape, as seen from the observer and calculated with planet_3d.planet_3d_coeffs().

• state (dict) – A dictionary containing all of the parameters needed to compute the phase curve. This includes the planet’s orbital parameters, the observer’s parameters, and the hotspot parameters.

• 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.

• xo (Array) – An offset to add to the x coordinate of the center of the planet, used when correcting for extended source illuminations. Default is 0.0.

• yo (Array) – An offset to add to the y coordinate of the center of the planet, used when correcting for extended source illuminations. Default is 0.0.

• zo (Array) – An offset to add to the z coordinate of the center of the planet, used when correcting for extended source illuminations. Default is 0.0.

Returns:

The timeseries of reflected light from the planet. Each element of the array corresponds to the time of the corresponding element in state[“times”].

Return type:

Array

phase_curve_utils.sample_surface(sample_radii, sample_thetas, rho_xx, rho_xy, rho_x0, rho_yy, rho_y0, rho_00, p_xx, p_xy, p_xz, p_x0, p_yy, p_yz, p_y0, p_zz, p_z0, p_00, **kwargs)

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

phase_curve_utils.stellar_doppler_variations(true_anomalies, stellar_doppler_alpha, period)

Compute the contributions to a phase curve for a star with Doppler variations.

A simple sinusoid model with a phase of 90 degrees at primary transit meant to capture Doppler boosting/flux falling in and out of the bandpass.

Parameters:

• true_anomalies (Array) – The true anomaly of the planet at each time step.

• stellar_doppler_alpha (float) – The amplitude of the Doppler variations.

• period (float) – The orbital period of the planet.

Returns:

The contribution to the phase curve from the star’s Doppler variations.

Return type:

Array

phase_curve_utils.stellar_ellipsoidal_variations(true_anomalies, stellar_ellipsoidal_alpha, period)

Compute the contributions to a phase curve for a star with ellipsoidal variations.

A simple sinusoid model with minima at primary and secondary eclipse, meant to capture gravitational Works only for a circular orbit and assumes that \(\Omega=\pi\). Uses the model in Shporer et al. 2014.

Technically the amplitude in Morris, Heng, and Kitzmann 2024 is given by

\[A_{ellip} = \frac{\alpha}{0.077} \frac{M_p}{M_J} \left(\frac{R_s}{R_\odot}\right)^3 \left(\frac{P}{1 \text{day}}\right)^{-2}\]

But we instead roll everything into the alpha parameter.

Parameters:

• true_anomalies (Array) – The true anomaly of the planet at each time step.

• stellar_ellipsoidal_alpha (float) – The amplitude of the ellipsoidal variations.

• period (float) – The orbital period of the planet.

Returns:

The contribution to the phase curve from the star’s ellipsoidal variations.

Return type:

Array

phase_curve_utils.surface_star_cos_angle(planet_surface_normal, x_c, y_c, z_c, **kwargs)

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, ys, c_x1, c_x2, c_x3, c_y1, c_y2, c_y3)

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.poly_to_parametric(rho_xx, rho_xy, rho_x0, rho_yy, rho_y0, rho_00)

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, rho_xy, rho_x0, rho_yy, rho_y0, rho_00, **kwargs)

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

kepler.kepler(M, ecc)

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, e, f, Omega, i, omega, **kwargs)

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

greens_basis_transform.generate_change_of_basis_matrix(N)

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