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 wavelength-dependent maps and light curves. This class may be particularly useful for users who wish to integrate starry 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 more kepler.Secondary objects and feed them into a kepler.System instance for integration with the Keplerian solver. The kepler.Primary and kepler.Secondary classes inherit from Map, 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 two-body systems or systems where the secondary masses are negligible. Stay tuned for the next release of the code, which will incorporate an N-body 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 the l, m index of the map instance. If map 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 disk-integrated 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 limb-darkened maps, whose disk-integrated intensity is always equal to unity.

Warning

Note that the total degree of the map can never exceed lmax. For example, if lmax=5 and you set spherical harmonic coefficients up to lmax=3, you may only set limb darkening coefficients up to second order.

As of version 0.2.1, starry supports multi-wavelength maps via the nwav 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 of l, m and all limb darkening coefficients at a given value of l are now vectors, with one value per wavelength bin. Functions like flux() 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 multi-precision to perform all calculations? Default False. If True, defaults to 32-digit (approximately 128-bit) floating point precision. This can be adjusted by changing the STARRY_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:
  • theta (float or ndarray) – Angle of rotation in degrees. Default 0.
  • x (float or ndarray) – Position scalar or vector.
  • y (float or ndarray) – Position scalar or vector.
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 is True, returns the tuple (F, dF), where F is the flux and dF 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 about axis. 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:
  • cmap (str) – The matplotlib colormap name. Default plasma.
  • res (int) – The resolution of the map in pixels on a side. Default 300.

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 a gif 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. Default plasma.
  • 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.
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 two-dimensional numpy array, structured such that the pixel at index (0, 0) is at latitude \(+90^\circ\) and longitude \(180^\circ W\).

Parameters:
  • image (str or array) – The full path to the image file, or a 2D numpy array of floats.
  • lmax (int) – The maximum degree of the spherical harmonic expansion of the image. Default lmax.

Note

For maps with nwav > 1, users may specify a nwav 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 two-dimensional numpy array, structured such that the pixel at index (0, 0) is at latitude \(+90^\circ\) and longitude \(180^\circ W\).

Parameters:
  • image (str or array) – The full path to the image file, or a 2D numpy array of floats.
  • lmax (int) – The maximum degree of the spherical harmonic expansion of the image. Default lmax.

Note

For maps with nwav > 1, users may specify a nwav 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 degree lmax, and sets the map vector.

Parameters:
  • image (ndarray) – The ring-ordered healpix array.
  • lmax (int) – The maximum degree of the spherical harmonic expansion of the image. Default lmax.

Note

For maps with nwav > 1, users may specify a nwav 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 gaussian-like 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. Read-only.

nwav

The number of wavelength bins. Read-only.

multi

Are calculations done using multi-precision? Read-only.

N

The number of map coefficients, equal to \((l + 1)^2\). Read-only.

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}, ...\}\). Read-only.

u

The limb darkening map vector. This is a vector of the limb darkening coefficients \(\{u_1, u_2, u_3, ...\}\). Read-only.

p

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

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, ...\}\). Read-only.

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. Read-only.

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. Read-only.

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 of Map 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 non-rotating 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. For nwav > 1, this is a matrix whose columns are the timeseries in each wavelength bin.

Note

Users must call the compute method of the System 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 current System object, formatted as body.parameter, where body is A for the primary and b, c, d, etc. for the secondaries. The parameter label is the name of the parameter.

Note

Users must call the compute method of the System object with gradient=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 limb-darkened 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 limb-darkening, 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 second-order 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 of Map 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 edge-on 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 non-rotating 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. For nwav > 1, this is a matrix whose columns are the timeseries in each wavelength bin.

Note

Users must call the compute method of the System 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 current System object, formatted as body.parameter, where body is A for the primary and b, c, d, etc. for the secondaries. The parameter label is the name of the parameter.

Note

Users must call the compute method of the System object with gradient=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 limb-darkened 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 limb-darkening, for which the derivatives with respect to the limb darkening coefficients are not computed by default.

a

The semi-major 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, edge-on orbit, a transit occurs when \(\lambda_0 = 90^\circ\).

X

The vector of x positions of the body in units of the primary radius. Read-only.

Note

Users must call the compute method of the System object before accessing this attribute.

Y

The vector of y positions of the body in units of the primary radius. Read-only.

Note

Users must call the compute method of the System object before accessing this attribute.

Z

The vector of z positions of the body in units of the primary radius. Read-only.

Note

Users must call the compute method of the System 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 (Primary) – The primary body. This body has unit radius and luminosity and is fixed at the origin.
  • secondaries (Secondary) – The secondary body, or a sequence of secondaries.
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 in lightcurve. The light curve for each body in the system is stored in the body’s lightcurve attribute. Optionally, also compute the gradient of the light curve and store it in the gradient attribute of the system and each of the body instances.

Parameters:
  • time (ndarray) – Time array, measured in days.
  • gradient (bool) – Compute the gradient of the light curve with respect to all body parameters? Default False
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. For nwav > 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 of Primary.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).