The starry Python API¶
Contents¶
Introduction¶
This page documents the starry
API, which is coded
in C++ with a pybind11
Python interface. The API consists
of a Map
class, which houses all of the surface map photometry
stuff, and a kepler
module, which connects photometry to
dynamics via a simple (multi)Keplerian solver.
There are two broad ways in which users can access
the core starry
functionality:
 Users can instantiate a
Map
class to compute phase curves and occultation light curves by directly specifying the rotational state of the object and (optionally) the position and size of an occultor. Users can specify spherical harmonic coefficients, limb darkening coefficients, or both. There is also support for both monochromatic and wavelengthdependent maps and light curves. This class may be particularly useful for users who wish to integratestarry
with their own dynamical code or for users wishing to compute simple light curves without any orbital solutions. Users can instantiate a
kepler.Primary
and one or morekepler.Secondary
objects and feed them into akepler.System
instance for integration with the Keplerian solver. Thekepler.Primary
andkepler.Secondary
classes inherit fromMap
, so users have access to all of the functionality of the previous method.
At present, starry
uses a simple Keplerian solver to compute orbits,
so the second approach listed above is limited to twobody systems or systems
where the secondary masses are negligible. Stay tuned for the next release
of the code, which will incorporate an Nbody solver.
The Map class¶

class
starry.
Map
¶ Instantiate a
starry
surface map. The map is described as an expansion in spherical harmonics, with optional arbitrary order limb darkening. Users can set the spherical harmonic coefficients by direct assignment to thel, m
index of the map instance. Ifmap
is an instance of this class,map[1, 0] = 0.5
sets the coefficient of the \(Y_{1,0}\) harmonic to \(\frac{1}{2}\). Users can set limb darkening coefficients by direct assignment to the
l
index of the map:map[1] = 0.4
sets the first order limb darkening coefficient \(u_1\) to \(0.4\).
Note
Map instances are normalized such that the average diskintegrated intensity is equal to the coefficient of the \(Y_{0,0}\) term, which defaults to unity. The total luminosity over all \(4\pi\) steradians is therefore four times the \(Y_{0,0}\) coefficient. This normalization is particularly convenient for constant or purely limbdarkened maps, whose diskintegrated intensity is always equal to unity.
Warning
Note that the total degree of the map can never exceed
lmax
. For example, iflmax=5
and you set spherical harmonic coefficients up tolmax=3
, you may only set limb darkening coefficients up to second order.As of version 0.2.1,
starry
supports multiwavelength maps via thenwav
keyword argument. If this argument is set to a value greater than one, light curves will be computed for each wavelength bin of the map. All coefficients at a given value ofl, m
and all limb darkening coefficients at a given value ofl
are now vectors, with one value per wavelength bin. Functions likeflux()
and__call__()
will return one value per wavelength bin, and gradients will be computed at every wavelength. Note that this feature is still experimental, so please raise an issue on GitHub if you run into any trouble.Parameters:  lmax (int) – Largest spherical harmonic degree in the surface map. Default 2.
 nwav (int) – Number of map wavelength bins. Default 1.
 multi (bool) – Use multiprecision to perform all calculations? Default
False
. IfTrue
, defaults to 32digit (approximately 128bit) floating point precision. This can be adjusted by changing theSTARRY_NMULTI
compiler macro.

__call__
(theta=0, x=0, y=0)¶ Return the specific intensity at a point
(x, y)
on the map. Users may optionally provide a rotation state. Note that this does not rotate the base map.Parameters: Returns: The specific intensity at
(x, y)
.

flux
(theta=0, xo=0, yo=0, ro=0, gradient=False)¶ Return the total flux received by the observer from the map. Computes the total flux received by the observer during or outside of an occultation.
Parameters:  theta (float or ndarray) – Angle of rotation. Default 0.
 xo (float or ndarray) – The
x
position of the occultor (if any). Default 0.  yo (float or ndarray) – The
y
position of the occultor (if any). Default 0.  ro (float) – The radius of the occultor in units of this body’s radius. Default 0 (no occultation).
 gradient (bool) – Compute and return the gradient of the flux as well? Default
False
.
Returns: The flux received by the observer (a scalar or a vector). If
gradient
isTrue
, returns the tuple(F, dF)
, whereF
is the flux anddF
is a dictionary containing the derivatives with respect to each of the input parameters and each of the map coefficients.

rotate
(theta=0)¶ Rotate the base map an angle
theta
aboutaxis
. This performs a permanent rotation to the base map. Subsequent rotations and calculations will be performed relative to this rotational state.Parameters: theta (float or ndarray) – Angle of rotation in degrees. Default 0.

show
(cmap='plasma', res=300)¶ Convenience routine to quickly display the body’s surface map.
Parameters: Note
For maps with
nwav > 1
, this method displays an animated sequence of frames, one per wavelength bin. Users can save this to disk by specifying agif
keyword with the name of the GIF image to save to.

animate
(cmap='plasma', res=150, frames=50, interval=75, gif='')¶ Convenience routine to animate the body’s surface map as it rotates.
Parameters:  cmap (str) – The
matplotlib
colormap name. Defaultplasma
.  res (int) – The resolution of the map in pixels on a side. Default 150.
 frames (int) – The number of frames in the animation. Default 50.
 interval (int) – Interval in milliseconds between frames. Default 75.
 gif (str) – The name of the .gif file to save the animation to. If set, does not show the animation. Default
None
.
 cmap (str) – The

load_image
(image, lmax=None)¶ Load an array or an image from file. This routine loads an image file, computes its spherical harmonic expansion up to degree
lmax
, and sets the map vector. Alternatively, users may provide a twodimensionalnumpy
array, structured such that the pixel at index(0, 0)
is at latitude \(+90^\circ\) and longitude \(180^\circ W\).Parameters: Note
For maps with
nwav > 1
, users may specify anwav
keyword argument indicating the wavelength bin into which the image or array will be loaded.Load an array or an image from file. This routine loads an image file, computes its spherical harmonic expansion up to degree
lmax
, and sets the map vector. Alternatively, users may provide a twodimensionalnumpy
array, structured such that the pixel at index(0, 0)
is at latitude \(+90^\circ\) and longitude \(180^\circ W\).Parameters: Note
For maps with
nwav > 1
, users may specify anwav
keyword argument indicating the wavelength bin into which the image or array will be loaded.

load_healpix
(image, lmax=None)¶ Load a healpix image array. This routine loads a
healpix
array, computes its spherical harmonic expansion up to degreelmax
, and sets the map vector.Parameters: Note
For maps with
nwav > 1
, users may specify anwav
keyword argument indicating the wavelength bin into which the image or array will be loaded.

add_gaussian
(sigma=0.1, amp=1, lat=0, lon=0, lmax=None)¶ Add the spherical harmonic expansion of a gaussian to the current map. This routine adds a gaussianlike feature to the surface map by computing the spherical harmonic expansion of a 3D gaussian constrained to the surface of the sphere. This is useful for, say, modeling star spots or other discrete, localized features on a body’s surface.
Note
Because this routine wraps a Python function, it is slow and should probably not be used repeatedly when fitting a map to data!
Warning
This routine is still in beta and may change in the near future.
Parameters:  sigma (float) – The standard deviation of the gaussian. Default 0.1
 amp (float) – The amplitude. Default 1.0, resulting in a gaussian whose integral over the sphere is unity.
 lat (float) – The latitude of the center of the gaussian in degrees. Default 0.
 lon (float) – The longitude of the center of the gaussian in degrees. Default 0.

reset
()¶ Set all of the map coefficients to zero, except for \(Y_{0,0}\), which is set to unity.

lmax
¶ The highest spherical harmonic degree of the map. Readonly.

nwav
¶ The number of wavelength bins. Readonly.

multi
¶ Are calculations done using multiprecision? Readonly.

N
¶ The number of map coefficients, equal to \((l + 1)^2\). Readonly.

y
¶ The spherical harmonic map vector. This is a vector of the coefficients of the spherical harmonics \(\{Y_{0,0}, Y_{1,1}, Y_{1,0}, Y_{1,1}, ...\}\). Readonly.

u
¶ The limb darkening map vector. This is a vector of the limb darkening coefficients \(\{u_1, u_2, u_3, ...\}\). Readonly.

p
¶ The polynomial map vector. This is a vector of the coefficients of the polynomial basis \(\{1, x, z, y, x^2, xz, ...\}\). Readonly.

g
¶ The Green’s polynomial map vector. This is a vector of the coefficients of the polynomial basis \(\{1, 2x, z, y, 3x^2, 3xz, ...\}\). Readonly.

r
¶ The current rotation solution vector r. Each term in this vector corresponds to the total flux observed from each of the terms in the polynomial basis. Readonly.

s
¶ The current occultation solution vector s. Each term in this vector corresponds to the total flux observeed from each of the terms in the Green’s basis. Readonly.
Note
For pure linear and quadratic limb darkening, the full solution vector is not computed, so this vector will not necessarily reflect the true solution coefficients.

axis
¶ A normalized unit vector specifying the default axis of rotation for the map. Default \(\hat{y} = (0, 1, 0)\).
The kepler module¶
The starry.kepler
module allows users to incorporate the
features of the starry.Map
class into a simple Keplerian
solver to model systems of stars, planets, or moons. The bodies in the
starry.kepler
module inherit from starry.Map
, so
all of the functionality outlined above is also available here.

class
starry.kepler.
Primary
¶ Instantiate a primary body. This body is assumed to be fixed at the origin. This class inherits from
Map
, so users can assign to and retrieve spherical harmonic and limb darkening coefficients in the same way. Refer to the documentation ofMap
for all options.
r
¶ The radius of the primary body, fixed at unity.

L
¶ The luminosity of the primary body, fixed at unity.

tref
¶ A reference time in days. Several of the orbital elements are specified at this time.

prot
¶ The rotation period of the body in days. For nonrotating bodies, set this to
np.inf
(or zero).

lightcurve
¶ The computed light curve for this body. If
nwav = 1
, this is a timeseries vector of fluxes. Fornwav > 1
, this is a matrix whose columns are the timeseries in each wavelength bin.Note
Users must call the
compute
method of theSystem
object before accessing this attribute.

gradient
¶ The gradient of the body’s light curve. This is a dictionary of vectors (
nwav = 1
) or matrices (nwav > 1
). The dictionary keys are the names of all parameters of all bodies in the currentSystem
object, formatted asbody.parameter
, wherebody
isA
for the primary andb
,c
,d
, etc. for the secondaries. Theparameter
label is the name of the parameter.Note
Users must call the
compute
method of theSystem
object withgradient=True
before accessing this attribute.Note
The gradients with respect to map coefficients are stored as vectors in the y and u attributes of the gradient dictionary, starting with the first degree coefficients. The derivatives with respect to the constant terms \(Y_{0,0}\) and \(u_0\) are not computed, since those values are constant!
Note
Depending on the properties of a body’s map, not all map coefficient gradients may be computed. For instance, for purely limbdarkened maps (whose \(Y_{l,m}\) coefficients are zero for \(l > 1\)), the derivatives of the flux with respect to the spherical harmonic coefficients are not computed. This is entirely in the interest of speed. To force the code to compute these, set one of the spherical harmonic coefficients to a very small (i.e., \(10^{15}\)) value. A similar caveat applies to maps with no limbdarkening, for which the derivatives with respect to the limb darkening coefficients are not computed by default.

r_m
¶ The radius of the primary body in meters. This is used exclusively for calculating the light travel time delay in the system. The default value is zero, in which case the speed of light is effectively infinite and there is no time delay. When this parameter is set, the time delay is computed with a secondorder Taylor expansion, which should be accurate enough for most applications. Note that the reference point for the time delay (where \(\Delta t = 0\)) is the barycenter of the system. Transits of bodies across the primary will therefore occur earlier, while occultations will occur later, than if there were no delay.


class
starry.kepler.
Secondary
¶ Instantiate a secondary body. This body is assumed to be massless and orbits the primary in a pure Keplerian orbit. This class inherits from
Map
, so users can assign to and retrieve spherical harmonic and limb darkening coefficients in the same way. Refer to the documentation ofMap
for all options.Note
When specifying the map coefficients of a
Secondary
body, we need to be careful about the frame in which the map is defined.starry
was developed with exoplanet secondary eclipses in mind, so the map coefficients correspond to the map that would be visible at opposition, i.e., at “secondary eclipse” for an exoplanet in an edgeon orbit. For bodies in inclined/rotated orbits, the map coefficients correspond to the map that would be visible to an observer located down the \(z\) axis, along the plane of the orbit. This means the map coefficients are independent of the inclination or longitude of ascending node of the orbit. If this is all very confusing, check out the tutorial on orbital geometry.
r
¶ The radius of the body in units of the primary’s radius. Default 0.1.

L
¶ The luminosity of the body in units of the primary’s luminosity. Default 0.0.

tref
¶ A reference time in days. Several of the orbital elements are specified at this time.

prot
¶ The rotation period of the body in days. For nonrotating bodies, set this to
np.inf
(or zero).

lightcurve
¶ The computed light curve for this body. If
nwav = 1
, this is a timeseries vector of fluxes. Fornwav > 1
, this is a matrix whose columns are the timeseries in each wavelength bin.Note
Users must call the
compute
method of theSystem
object before accessing this attribute.

gradient
¶ The gradient of the body’s light curve. This is a dictionary of vectors (
nwav = 1
) or matrices (nwav > 1
). The dictionary keys are the names of all parameters of all bodies in the currentSystem
object, formatted asbody.parameter
, wherebody
isA
for the primary andb
,c
,d
, etc. for the secondaries. Theparameter
label is the name of the parameter.Note
Users must call the
compute
method of theSystem
object withgradient=True
before accessing this attribute.Note
The gradients with respect to map coefficients are stored as vectors in the y and u attributes of the gradient dictionary, starting with the first degree coefficients. The derivatives with respect to the constant terms \(Y_{0,0}\) and \(u_0\) are not computed, since those values are constant!
Note
Depending on the properties of a body’s map, not all map coefficient gradients may be computed. For instance, for purely limbdarkened maps (whose \(Y_{l,m}\) coefficients are zero for \(l > 1\)), the derivatives of the flux with respect to the spherical harmonic coefficients are not computed. This is entirely in the interest of speed. To force the code to compute these, set one of the spherical harmonic coefficients to a very small (i.e., \(10^{15}\)) value. A similar caveat applies to maps with no limbdarkening, for which the derivatives with respect to the limb darkening coefficients are not computed by default.

a
¶ The semimajor axis of the body in units of the primary radius. Default 50.

porb
¶ The orbital period of the body in days. Default 1.

inc
¶ The inclination of the body in degrees. Default 90.

ecc
¶ The eccentricity of the body. Default 0.

w
¶ The longitude of pericenter for the body’s orbit in degrees. This parameter is typically denoted \(\varpi\). Default 90.

Omega
¶ The longitude of ascending node in degrees. Default 0.

lambda0
¶ The mean longitude of the body in degrees at the reference time. Default 90. Note that for a circular, edgeon orbit, a transit occurs when \(\lambda_0 = 90^\circ\).

X
¶ The vector of
x
positions of the body in units of the primary radius. Readonly.Note
Users must call the
compute
method of theSystem
object before accessing this attribute.


class
starry.kepler.
System
¶ Instantiate a Keplerian orbital system. The primary is fixed at the origin, and all secondary bodies are assumed to be massless.
Parameters: 
primary
¶ The primary body in the Keplerian system.

secondaries
¶ A list of the secondary bodies in the Keplerian system.

compute
(time, gradient=False)¶ Compute the system light curve analytically. Compute the full system light curve at the times given by the
time
array and store the result inlightcurve
. The light curve for each body in the system is stored in the body’slightcurve
attribute. Optionally, also compute the gradient of the light curve and store it in thegradient
attribute of the system and each of the body instances.Parameters:

lightcurve
¶ The computed light curve for the system, equal to the sum of the light curves of each of the bodies. If
nwav = 1
, this is a timeseries vector of fluxes. Fornwav > 1
, this is a matrix whose columns are the timeseries in each wavelength bin.Note
Users should call
compute()
first.

gradient
¶ The gradient of the systems’s light curve. This is a dictionary of vectors (
nwav = 1
) or matrices (nwav > 1
). See the docstring ofPrimary.gradient
for more details.Note
Users should call
compute()
first.

exposure_time
¶ Exposure time for each data point in the light curve. Default 0. If nonzero, integrates the light curve over the exposure time using a recursive technique to approximate the integral.

exposure_tol
¶ Tolerance of the recursive exposure time algorithm. Default is square root of machine epsilon.

exposure_max_depth
¶ Maximum number of recursions in the exposure time algorithm. Default 4. Increase this for higher accuracy (and longer run times).
