- class starry.System¶
A system of bodies in Keplerian orbits about a central primary body.
Primary) – The central body.
Secondary) – One or more secondary bodies in orbit about the primary.
time_unit (optional) – An
astropy.unitsunit defining the time metric for this object. Defaults to
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
tis assumed to indicate the timestamp at the middle of an exposure of length
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:
0for a centered Riemann sum (equivalent to the “resampling” procedure suggested by Kipping 2010),
1for the trapezoid rule, or
2for Simpson’s rule.
- property bodies¶
A list of all objects in the Keplerian system.
Compute the system flux design matrix at times
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
t (scalar or vector) – An array of times at which to evaluate the design matrix in units of
Draw a map from the posterior distribution and set the
ymap vector of each body.
Users should call
solve()to enable this attribute.
- flux(t, total=True)¶
Compute the system flux at times
t (scalar or vector) – An array of times at which to evaluate the flux in units of
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
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
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.
The log marginal likelihood.
- Return type
- 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
Compute the Cartesian positions of all bodies at times
t (scalar or vector) – An array of times at which to evaluate the position in units of
- 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 (scalar or vector) – An array of times at which to evaluate the radial velocity in units of
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.
A vector or matrix of radial velocities in units of meters per second. If
True, this returns a vector of the total stellar radial velocity at every point
t. 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. If
False, 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 (one
Secondarybodies), 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 the
Secondarybodies, 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 if
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
Secondarybodies are treated as massive only when computing their gravitational effect on the
Primary(as opposed to each other). This therefore ignores all
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.
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.
t (scalar or vector) – The time(s) at which to evaluate the orbit and the map in units of
cmap (string or colormap instance, optional) – The matplotlib colormap to use. Defaults to
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.
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.
If calling this method on an instance of
Systemcreated within a
pymc3.Model(), you may specify a
pointkeyword with the model point at which to evaluate the map. This method also accepts a
modelkeyword, although this is inferred automatically if called from within a
pymc3.Model()context. If no point is provided, attempts to evaluate the map at
model.test_pointand 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
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.
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¶
astropy.unitsunit defining the time metric for the system.