OblateSystem#

This class is intended to be the primary user interface for generating transit, eclipse, and phase curve models for oblate planets. Users specify parameters describing a system, then vary those parameters through built-in functions to generate different models and compare them to data.

class oblate_system.OblateSystem(times=None, t_peri=None, period=None, a=None, tidally_locked=None, e=0.0, i=1.5707963267948966, Omega=3.141592653589793, omega=0.0, obliq=0.0, prec=0.0, r=None, f1=0.0, f2=0.0, ld_u_coeffs=Array([0., 0.], dtype=float64), hotspot_latitude=0.0, hotspot_longitude=0.0, hotspot_concentration=0.2, albedo=1.0, emitted_scale=1e-05, stellar_ellipsoidal_alpha=1e-06, stellar_doppler_alpha=1e-06, systematic_trend_coeffs=Array([0., 0.], dtype=float64), log_jitter=-inf, projected_effective_r=0.0, projected_f=0.0, projected_theta=0.0, extended_illumination_npts=1, compute_reflected_phase_curve=False, compute_emitted_phase_curve=False, compute_stellar_ellipsoidal_variations=False, compute_stellar_doppler_variations=False, parameterize_with_projected_ellipse=False, phase_curve_nsamples=50000, random_seed=0, data=Array([1.], dtype=float64), uncertainties=Array([0.01], dtype=float64), exposure_time=0.0, oversample=1, oversample_correction_order=2)[source]#

Bases: object

The core user interface for squishyplanet, used to model potentially-triaxial exoplanet transits/phase curves.

Note, all instances will have values associated with phase curve calculations, such as albedo and hotspot location. However, if inputs such are “compute_reflected_phase_curve” are set to False, these values will not be used. The lightcurve() method will return only the transit light curve in this case, and should be used if computing a transit, reflected, or emitted phase curve.

All arguments will be internally converted to jax.numpy dtypes, and all methods will similarly return jax.numpy arrays. These can be treated similarly to numpy arrays in most cases, but if passing outputs to external inference libraries expecting numpy, you may need to explicitly convert them.

Properties:

state (dict):: A dictionary of all the parameters of the system, including those specified by the user, default values, and those calculated by combinations of the two. Immutable, but can be accessed to see the current state of the system.

Parameters:

times (array-like, [Days], default=None) – The times at which to calculate the light curve. The gap between times is assumed to be in units of days, but any zero-point/standard system (e.g. BJD) will work. A required parameter, will raise an error if not provided.
t_peri (float, [Days], default=None) – The time of periastron passage. A required parameter, will raise an error if not provided.
period (float, [Days], default=None) – The period of the orbit. A required parameter, will raise an error if not provided.
a (float, [Rstar], default=None) – The semi-major axis of the orbit in units of the radius of the star. A required parameter, will raise an error if not provided.
tidally_locked (bool, default=None) – Whether the planet is tidally locked to the star. If True, then prec will always be set equal to the true anomaly, meaning the same face of the planet will always face the star. A required parameter, will raise an error if not provided.
e (float, default=0.0) – The eccentricity of the orbit.
i (float, [Radian], default=jnp.pi / 2) – The inclination of the orbit.
Omega (float, [Radian], default=jnp.pi) – The longitude of the ascending node. Changing this will not affect the transit light curve (more accurately, changes can always be compensated for via rotations in obliq or prec). It is included only because it naturally arises in the orbit rotations, and I guess could come into play if anyone ever wants to do a joint astrometry model.
omega (float, [Radian], default=0.0) – The argument of periapsis. Set to 0.0 for a circular orbit, otherwise there will be degeneracies with t_peri.
obliq (float, [Radian], default=0.0) – The obliquity of the planet. This is the angle between the planet’s rotation axis and the normal to the orbital plane. It is defined such that a planet on a circular orbit with \(\Omega = 0\) and \(\nu = 0\) (i.e., when it’s along the positive \(x\) axis) will have its north pole tipped away from the star.
prec (float, [Radian], default=0.0) – The “precession angle” of the planet. This defined as a rotation of the planet about an axis that’s aligned with its orbit normal and runs through the center of the planet (e.g., if obliq=0, it would set the planet’s instantaneous rotational phase, and if obliq \(\neq\) 0, it would set the “season” of the northern hemisphere at periastron passage.)
r (float, [Rstar], default=None) – The equatorial radius of the planet. This will always be the largest of the 3 axes of the triaxial ellipsoid. Either this or the entire set of projected_effective_r, projected_f, and projected_theta must be provided.
f1 (float, [Dimensionless], default=0.0) – The fractional difference between the (longest) equatorial and polar radii of the planet. This is defined as \((R_{eq} - R_{pol}) / R_{eq}\).
f2 (float, [Dimensionless], default=0.0) – The fractional difference between long and short radii of the ellipse that defines the equator of the planet. Defined similarly to f1.
ld_u_coeffs (array-like, default=jnp.array([0.0, 0.0])) –
The coefficients that determine the limb darkening profile of the star. The star is assumed to be azimuthally symmetric and have a radial profile described by:

\[\frac{I(\mu)}{I_0} = - \Sigma_{i=0}^N u_i (1 - \mu)^i\]

for some order polynomial \(N\). See Agol, Luger, and Foreman-Mackey 2020 for more.
hotspot_latitude (float, [Radian], default=0.0) – The latitude of a potential hotspot on the planet. This is defined according to the “physics” convention of spherical coordinates, not in the geography sense: 0 is the north pole, \(\pi/2\) is the equator, and \(\pi\) is the south pole.
hotspot_longitude (float, [Radian], default=0.0) – The longitude of a potential hotspot on the planet.
hotspot_concentration (float, default=0.2) – The “concentration” of the hotspot. This is the \(\kappa\) parameter in the von Mises-Fisher distribution that describes the hotspot.
albedo (float, default=1.0) – The (spatialy uniform) albedo of the planet. This is the fraction of light that is reflected, though the directional-dependent scattering is dictated by Lambert’s cosine law.
emitted_scale (float, default=1e-5) – A scale factor that sets the amplitude of the the emitted flux of the planet.
systematic_trend_coeffs (array-like, default=jnp.array([0.0,0.0])) – The coefficients that determine the polynomial trend in time added to the lightcurves. Used to optionally model long-term drifts in observed data.
log_jitter (float, default=-jnp.inf) – The log of the “jitter” term included in likelihood calculations. The jitter is added in quadrature to the provided uncertainties to account for any unmodeled noise in the data. This value is the standard deviation of the jitter, not the variance. If set to -jnp.inf, the jitter term will not affect the likelihood.
projected_effective_r (float, [Rstar], default=0.0) – The radius of a circle with the same area is the projected ellipse. This is only relevant if parameterize_with_projected_ellipse is set to True, which will override r, f1, f2, obliq, and prec.
projected_f (float, [Dimensionless], default=0.0) – The flattening value of the projected ellipse. This is only relevant if parameterize_with_projected_ellipse is set to True, which will override r, f1, f2, obliq, and prec.
projected_theta (float, [Radian], default=0.0) – The angle of the semi-major axis of the projected ellipse. This is only relevant if parameterize_with_projected_ellipse is set to True, which will override r, f1, f2, obliq, and prec.
extended_illumination_npts (int, default=1) – The number of points used to sample the star’s projected disk as seen by the planet when calculating the reflected flux and accounting for the star’s extended size. Closely follows the implementation in starry, specifically Sec. 4.1 of Luger et al. 2022. IMPLEMENTATION IS INCOMPLETE, WILL RAISE A NOTIMPLEMENTEDERROR IF SET TO ANYTHING OTHER THAN 1.
compute_reflected_phase_curve (bool, default=False) – Whether to include flux reflected by the planet when calling lightcurve().
compute_emitted_phase_curve (bool, default=False) – Whether to include flux emitted by the planet when calling lightcurve().
compute_stellar_ellipsoidal_variations (bool, default=False) – Whether to include stellar ellipsoidal variations in the light curve. This is the effect of the star’s shape changing due to the gravitational pull of the planet, and here is modeled as a simple sinusoidal variation with 4 peaks per orbit.
compute_stellar_doppler_variations (bool, default=False) – Whether to include stellar doppler variations in the light curve. This captures the effects of the star’s radial velocity changing and boosting the total flux/pushing some flux into/out of the bandpass of the observation. Here, it is modeled as a simple sinusoidal variation with 2 peaks per orbit.
parameterize_with_projected_ellipse (bool, default=False) – Whether to parameterize the planet as a projected ellipse rather than a triaxial ellipsoid. If True, then projected_effective_r, projected_f, and projected_theta will be used.
phase_curve_nsamples (int, default=50_000) – The number of random samples of the planet’s surface to draw when performing Monte Carlo estimates of the emitted/reflected flux. A larger number will increase the resolution/shrink the error of the estimate but result in longer computation times.
random_seed (int, default=0) – A random seed used for the Monte Carlo integrals in the phase curve. This feeds into jax.random.PRNGKey. Runs with the same random_seed will always return identical outputs, so if checking the affect of altering phase_curve_nsamples, you should change this as well.
data (array-like, default=jnp.array([1.0])) – The observed data to compare to the light curve. Must be the same length as times. Only needed if calling loglike().
uncertainties (array-like, default=jnp.array([0.01])) – The uncertainties on the observed data. Must be the same length as data, even if the errors are homoskedastic. Only needed if calling loglike().
exposure_time (float, [Days], default=0.0) – The length of each exposure in the light curve, used to correct for finite integration times if oversample is set to a value greater than 1. Important: the finite exposure time correction procedure assumes that the given times correspond to the midpoints of each exposure, not the start or end. No checks are made to ensure that it is shorter than the minimum time difference between the provided times.
oversample (int, default=1) – The factor by which to oversample the light curve to partially compensate for finite-time integrations. The overdense lightcurve is then binned down to the original provided times. See e.g. Kipping 2010, “Binning is Sinning” for more. Must be a positive integer. Will be rounded up to nearest odd number.
oversample_correction_order (int, default=2) – After oversampling the light curve, how do you want to integrate over the exposure time to get the final binned light curve? This follows starry’s treatment very closely: 0 is a centered Riemann sum like in Kipping 2010, 1 is a trapezoidal rule, and 2 is Simpson’s rule. Must be one of those values.

static fit_limb_darkening_profile(intensities, order=None, mus=None, rs=None)[source]#

Convert a stellar limb darkening profile to a polynomial representation.

Given a set of stellar parameters, one can use a grid of stellar models to compute the limb darkening profile as a function of projected r or of mu = sqrt(1 - r**2). These profiles are often then approximated with one of a few common limb darkening “laws”, such as the quadratic or 4-parameter non-linear laws. Since squishyplanet only supports polynomial limb darkening profiles, but can support nearly arbitrary orders, we can approximate the profile with a polynomial. This is a convenience function for converting between a grid-derived profile and its best-fit polynomial representation in the correct basis for squishyplanet.

Parameters:

intensities (array-like) – The relative intensities of the star at a given mu or r.
order (int) – The order of the polynomial to fit to the limb darkening profile. Note that in the squishyplanet basis, the polynomial is defined as 1 - sum_{i=1}^{order} u_i (1 - mu)^i, so the number of coefficients is order, not order+1.
mus (array-like, default=None) – The mu values at which the intensities were computed. If rs is not provided, this is required.
rs (array-like, default=None) – The r values at which the intensities were computed. If mus is not provided, this is required.

Returns:

The coefficients of the polynomial representation of the limb darkening profile. These can then be used as the ld_u_coeffs parameter in OblateSystem.

Return type:

array-like

illustrate(times=None, true_anomalies=None, orbit=True, reflected=False, emitted=False, star_fill=True, window_size=0.4, star_centered=False, nsamples=50000, figsize=(8, 8))[source]#

Visualize the layout of the system at one or more times.

This method, if run in a jupyter notebook, will display a plot of some combination of the star, planet, and its orbit. It can color in the planet according to its reflected or emission profile, and the star according to its limb darkening profile. Helpful for checking the orientation of planet hotspots and/or its orientation after deformation.

Parameters:

times (array-like, [Days], default=None) – The times at which to illustrate the system. The gap between times is assumed to be in units of days, but any zero-point/standard system (e.g. BJD) will work. Provide either this or true_anomalies but not both.
true_anomalies (array-like, [Radian], default=None) – The true anomalies at which to illustrate the system. Provide either this or times but not both.
orbit (bool, default=True) – Whether to plot a trace of the planet’s orbital path
reflected (bool, default=False) – Whether to color in the planet according to its reflected flux profile. Can optionally include this or emitted but not both.
emitted (bool, default=False) – Whether to color in the planet according to its emitted flux profile. Can optionally include this or reflected but not both.
star_fill (bool, default=True) – Whether to color in the star according to its limb darkening profile. Note that the lowest color contour is bounded at zero, so if you have an unphysical limb darkening law where some radii are negative, those will appear as gaps (most often the contours will not reach the black outline of the star, which is always drawn).
window_size (float, [Rstar], default=0.4) – The size of the plotting window. The window will be centered on the mean position of the planet across all of the suggested times, unless star_centered is set to True.
star_centered (bool, default=False) – Whether to center the plot on the star rather than the planet.
nsamples (int, default=50_000) – The number of random samples of the planet’s surface to draw when illustrating the system. A larger number will increase the resolution of the plot but result in longer computation times.
figsize (tuple, default=(8, 8)) – The size of the figure to display. Passed directly to matplotlib.pyplot.subplots.

Returns:

This method is used for its side effects of displaying a plot, not for its return value.

Return type:

None

lightcurve(params={})[source]#

Compute the light curve of the system.

This method will return the light curve of the system at the times specified when the system was initialized. If you want to compute the light curve at different times, or with different orbital parameters, you can pass those parameters as a dictionary to this method.

The first time this is run for a given system, JAX will jit-compile the function, which can take some time. Subsequent calls will be much faster unless you change the shape of any of the input arrays (e.g., changing the number of times or the order of the polynomial limb darkening law). In those cases, or if changing any of boolean flags, JAX will need to re-compile the function again.

Parameters:: params (dict, default={}) – A dictionary of parameters to update in the system state. Any keys not provided will be pulled from the current state of the system.
Returns:: The timeseries lightcurve of the system. The length will be equal to state[“times”], and each index corresponds to a time in that array.
Return type:: Array

Examples

>>> state = {
        "t_peri" : 0.0,
        "times" : jnp.linspace(-jnp.pi, 2*jnp.pi, 3504),
        "a" : 2.0,
        "period" : 2*jnp.pi,
        "r" : 0.1,
        "compute_reflected_phase_curve" : True,
        "compute_emitted_phase_curve" : True,
        "emitted_scale" : 1e-5,
    }
>>> system = OblateSystem(**state)
>>> system.lightcurve()

static limb_darkening_profile(ld_u_coeffs=None, r=None, mu=None)[source]#

Compute the limb darkening profile of the star at a given radius.

Meant as a helper function for sanity checks and plotting, especially if you’re using higher-order limb darkening laws and are concerned if the profile is positive/monotonic.

Parameters:

ld_u_coeffs (array-like, default=self.state["ld_u_coeffs"]) – The coefficients of the polynomial limb darkening law.
r (float or array-like, default=None) – The radius at which to compute the limb darkening profile. Must be between 0 and 1. If provided, mu should be None.
mu (float or array-like, default=None) – The cosine of the angle between the line of sight and the normal to the surface of the star. Must be between 0 and 1. If provided, r should be None.

Returns:

The limb darkening profile of the star at the given r or mu values.

Return type:

Array

loglike(params={})[source]#

Compute the log likelihood of the system given the observed data and some set of parameters.

This method will call lightcurve() with the provided parameters and compare the output to the observed data. The likelihood is assumed to be Gaussian with no correlation between times. The jitter term is added in quadrature to the provided uncertainties.

Parameters:: params (dict, default={}) – A dictionary of parameters to update in the system state. Any keys not provided will be pulled from the current state of the system.
Returns:: The log likelihood of the system given the observed data and the provided parameters.
Return type:: float

property state#

A dictionary that includes all of the parameters of the system.

This is an immutable property, and will raise an error if you try to set it. If altering parameters that would affect a lightcurve, pass those as a dictionary to the lightcurve() method. If altering the data or times at which to generate the lightcurve, just define a new system with those values.

Returns:: A dictionary of all the parameters of the system, including those specified by the user, default values, and those calculated by combinations of the two.
Return type:: dict

OblateSystem

Contents

OblateSystem#