Keplerian system

class starry.System

A system of bodies in Keplerian orbits about a central primary body.

Parameters

  • primary (Primary) – The central body.

  • secondaries (Secondary) – One or more secondary bodies in orbit about the primary.

  • time_unit (optional) – An astropy.units unit defining the time metric for this object. Defaults to astropy.units.day.

  • light_delay (bool, optional) – Account for the light travel time delay to the barycenter of the system? Default is False.

  • texp (scalar) – The exposure time of each observation. This can be a scalar or a tensor with the same shape as t. If texp is provided, t is assumed to indicate the timestamp at the middle of an exposure of length texp.

  • oversample (int) – The number of function evaluations to use when numerically integrating the exposure time.

  • order (int) – The order of the numerical integration scheme. This must be one of the following: 0 for a centered Riemann sum (equivalent to the “resampling” procedure suggested by Kipping 2010), 1 for the trapezoid rule, or 2 for Simpson’s rule.

property bodies

A list of all objects in the Keplerian system.

design_matrix(t)

Compute the system flux design matrix at times t.

Note

This is the unweighted design matrix, i.e., it does not include the scaling by the amplitude of each body’s map. To perform this weighting, do

X = sys.design_matrix(**kwargs)
for i, body in zip(sys.map_indices, sys.bodies):
    X[:, i] *= body.map.amp
Parameters

t (scalar or vector) – An array of times at which to evaluate the design matrix in units of time_unit.

draw()

Draw a map from the posterior distribution and set the y map vector of each body.

Users should call solve() to enable this attribute.

flux(t, total=True)

Compute the system flux at times t.

Parameters

  • t (scalar or vector) – An array of times at which to evaluate the flux in units of time_unit.

  • total (bool, optional) – Return the total system flux? Defaults to True. If False, returns arrays corresponding to the flux from each body.

property light_delay

Account for the light travel time delay? Read-only

lnlike(*, design_matrix=None, t=None, woodbury=True)

Returns the log marginal likelihood of the data given a design matrix.

This method computes the marginal likelihood (marginalized over the spherical harmonic coefficients of all bodies) given a system light curve and its covariance (set via the set_data() method) and a Gaussian prior on the spherical harmonic coefficients (set via the set_prior() method).

Parameters

  • design_matrix (matrix, optional) – The flux design matrix, the quantity returned by design_matrix(). Default is None, in which case this is computed based on kwargs.

  • t (vector, optional) – The vector of times at which to evaluate design_matrix(), if a design matrix is not provided. Default is None.

  • woodbury (bool, optional) – Solve the linear problem using the Woodbury identity? Default is True. The Woodbury identity is used to speed up matrix operations in the case that the number of data points is much larger than the number of spherical harmonic coefficients. In this limit, it can speed up the code by more than an order of magnitude. Keep in mind that the numerical stability of the Woodbury identity is not great, so if you’re getting strange results try disabling this. It’s also a good idea to disable this in the limit of few data points and large spherical harmonic degree.

Returns

The log marginal likelihood.

Return type

lnlike

property map_indices

A list of the indices corresponding to each body in the design matrix.

property order

The order of the numerical integration scheme. Read-only

  • 0: a centered Riemann sum

  • 1: trapezoid rule

  • 2: Simpson’s rule

property oversample

Oversample factor when integrating over exposure time. Read-only

position(t)

Compute the Cartesian positions of all bodies at times t.

Parameters

t (scalar or vector) – An array of times at which to evaluate the position in units of time_unit.

property primary

The primary (central) object in the Keplerian system.

rv(t, keplerian=True, total=True)

Compute the observed radial velocity of the system at times t.

Parameters

  • t (scalar or vector) – An array of times at which to evaluate the radial velocity in units of time_unit.

  • keplerian (bool) – Include the Keplerian component of the radial velocity of the primary? Default is True. If False, this method returns a model for only the radial velocity anomaly due to transits (the Rossiter-McLaughlin effect) and time-variable surface features (Doppler tomography) for all bodies in the system.

  • total (bool, optional) – Return the total system RV? Defaults to True. If False, returns arrays corresponding to the RV contribution from each body.

property secondaries

A list of the secondary (orbiting) object(s) in the Keplerian system.

set_data(flux, C=None, cho_C=None)

Set the data vector and covariance matrix.

This method is required by the solve() method, which analytically computes the posterior over surface maps for all bodies in the system given a dataset and a prior, provided both are described as multivariate Gaussians.

Parameters

  • flux (vector) – The observed system light curve.

  • C (scalar, vector, or matrix) – The data covariance. This may be a scalar, in which case the noise is assumed to be homoscedastic, a vector, in which case the covariance is assumed to be diagonal, or a matrix specifying the full covariance of the dataset. Default is None. Either C or cho_C must be provided.

  • cho_C (matrix) – The lower Cholesky factorization of the data covariance matrix. Defaults to None. Either C or cho_C must be provided.

show(t, cmap='plasma', res=300, interval=75, file=None, figsize=(3, 3), html5_video=True, window_pad=1.0)

Visualize the Keplerian system.

Parameters

  • t (scalar or vector) – The time(s) at which to evaluate the orbit and the map in units of time_unit.

  • cmap (string or colormap instance, optional) – The matplotlib colormap to use. Defaults to plasma.

  • res (int, optional) – The resolution of the map in pixels on a side. Defaults to 300.

  • figsize (tuple, optional) – Figure size in inches. Default is (3, 3) for orthographic maps and (7, 3.5) for rectangular maps.

  • interval (int, optional) – Interval between frames in milliseconds (animated maps only). Defaults to 75.

  • file (string, optional) – The file name (including the extension) to save the animation to (animated maps only). Defaults to None.

  • html5_video (bool, optional) – If rendering in a Jupyter notebook, display as an HTML5 video? Default is True. If False, displays the animation using Javascript (file size will be larger.)

  • window_pad (float, optional) – Padding around the primary in units of the primary radius. Bodies outside of this window will be cropped. Default is 1.0.

property solution

The posterior probability distribution for the maps in the system.

This is a tuple containing the mean and lower Cholesky factorization of the covariance of the amplitude-weighted spherical harmonic coefficient vectors, obtained by solving the regularized least-squares problem via the solve() method.

Note that to obtain the actual covariance matrix from the lower Cholesky factorization \(L\), simply compute \(L L^\top\).

Note also that this is the posterior for the amplitude-weighted map vectors. Under this convention, the map amplitude is equal to the first term of the vector of each body and the spherical harmonic coefficients are equal to the vector normalized by the first term.

solve(*, design_matrix=None, t=None)

Solve the least-squares problem for the posterior over maps for all bodies.

This method solves the generalized least squares problem given a system light curve and its covariance (set via the set_data() method) and a Gaussian prior on the spherical harmonic coefficients (set via the set_prior() method). The map amplitudes and coefficients of each of the bodies in the system are then set to the maximum a posteriori (MAP) solution.

Parameters

  • design_matrix (matrix, optional) – The flux design matrix, the quantity returned by design_matrix(). Default is None, in which case this is computed based on kwargs.

  • t (vector, optional) – The vector of times at which to evaluate design_matrix(), if a design matrix is not provided. Default is None.

Returns

The posterior mean for the spherical harmonic coefficients l > 0 and the Cholesky factorization of the posterior covariance of all of the bodies in the system, stacked in order (primary, followed by each of the secondaries in the order they were provided.)

Note

Users may call the draw() method of this class to draw from the posterior after calling solve().

property texp

The exposure time in units of time_unit. Read-only

property time_unit

An astropy.units unit defining the time metric for the system.