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 toastropy.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
. Iftexp
is provided,t
is assumed to indicate the timestamp at the middle of an exposure of lengthtexp
.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, or2
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 theset_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 onkwargs
.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 sum1
: trapezoid rule2
: 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.
- Returns
A vector or matrix of radial velocities in units of meters per second. If
total
isTrue
, this returns a vector of the total stellar radial velocity at every pointt
. This is the sum of all effects: the Keplerian motion of the star due to the planets, the radial velocity anomaly due to spots or features on its surface rotating in and out of view, as well as the Rossiter-McLaughlin effect due to transits by the secondaries. Iftotal
isFalse
, this returns a matrix whose rows are the radial velocity contributions to the measured stellar radial velocity from each of the bodies in the system. As a specific example, if there are three bodies in the system (onePrimary
and twoSecondary
bodies), the first row is the radial velocity anomaly corresponding to the star’s own Doppler signals (due to spots rotating in and out of view, or due to the Rossiter-McLaughlin effect of the other two bodies transiting it). The other two rows are the Keplerian contribution from each of theSecondary
bodies, plus any radial velocity anomaly due to spots or occultations of their surfaces. The sum across all rows is equal to the total radial velocity and is the same as what this method would return iftotal=True
.
Note, importantly, that when
total=False
, this method does not return the radial velocities of each of the bodies; instead, it returns the contribution to the stellar radial velocity due to each of the bodies. If you require knowing the radial velocity of the secondary objects, you can compute this from conservation of momentum:Note that this method implicitly assumes multi-Keplerian orbits; i.e., the
Secondary
bodies are treated as massive only when computing their gravitational effect on thePrimary
(as opposed to each other). This therefore ignores allSecondary
-Secondary
(i.e., planet-planet) interactions.
- 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, **kwargs)¶
Visualize the Keplerian system.
Note that the body surface intensities are not normalized.
- 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.
Note
If calling this method on an instance of
System
created within apymc3.Model()
, you may specify apoint
keyword with the model point at which to evaluate the map. This method also accepts amodel
keyword, although this is inferred automatically if called from within apymc3.Model()
context. If no point is provided, attempts to evaluate the map atmodel.test_point
and raises a warning.
- 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 theset_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 onkwargs
.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.)
- property time_unit¶
An
astropy.units
unit defining the time metric for the system.