Spherical harmonic maps¶

class starry.Map

The default starry map class.

This class handles light curves and phase curves of objects in emitted light.

Note

Instantiate this class by calling starry.Map() with ydeg > 0 and both rv and reflected set to False.

property N

Total number of map coefficients. Read-only

This is equal to \(N_\mathrm{y} + N_\mathrm{u} + N_\mathrm{f}\).

property Nf

Number of spherical harmonic coefficients in the filter. Read-only

This is equal to \((f_\mathrm{deg} + 1)^2\).

property Nu

Number of limb darkening coefficients, including \(u_0\). Read-only

This is equal to \(u_\mathrm{deg} + 1\).

property Ny

Number of spherical harmonic coefficients. Read-only

This is equal to \((y_\mathrm{deg} + 1)^2\).

add_spot(amp=None, intensity=None, relative=True, sigma=0.1, lat=0.0, lon=0.0)

Add the expansion of a gaussian spot to the map.

This function adds a spot whose functional form is the spherical harmonic expansion of a gaussian in the quantity \(\cos\Delta\theta\), where \(\Delta\theta\) is the angular separation between the center of the spot and another point on the surface. The spot brightness is controlled by either the parameter amp, defined as the fractional change in the total luminosity of the object due to the spot, or the parameter intensity, defined as the fractional change in the intensity at the center of the spot.

Parameters

amp (scalar or vector, optional) – The amplitude of the spot. This is equal to the fractional change in the luminosity of the map due to the spot. If the map has more than one wavelength bin, this must be a vector of length equal to the number of wavelength bins. Default is None. Either amp or intensity must be given.
intensity (scalar or vector, optional) – The intensity of the spot. This is equal to the fractional change in the intensity of the map at the center of the spot. If the map has more than one wavelength bin, this must be a vector of length equal to the number of wavelength bins. Default is None. Either amp or intensity must be given.
relative (bool, optional) – If True, computes the spot expansion assuming the fractional amp or intensity change is relative to the current map amplitude/intensity. If False, computes the spot expansion assuming the fractional change is relative to the original map amplitude/intensity (i.e., that of a featureless map). Defaults to True. Note that if True, adding two spots with the same values of amp or intensity will generally result in different intensities at their centers, since the first spot will have changed the map intensity everywhere! Defaults to True.
sigma (scalar, optional) – The standard deviation of the gaussian. Defaults to 0.1.
lat (scalar, optional) – The latitude of the spot in units of angle_unit. Defaults to 0.0.
lon (scalar, optional) – The longitude of the spot in units of angle_unit. Defaults to 0.0.

property alpha

The rotational shear coefficient, a number in the range [0, 1].

The parameter \(\alpha\) is used to model linear differential rotation. The angular velocity at a given latitude \(\theta\) is

\(\omega = \omega_{eq}(1 - \alpha \sin^2\theta)\)

where \(\omega_{eq}\) is the equatorial angular velocity of the object.

Note

This parameter is only used if drorder is greater than zero and/or radial velocity mode is enabled.

property amp: The overall amplitude of the map in arbitrary units. This factor multiplies the intensity and the flux and is thus proportional to the luminosity of the object. For multi-wavelength maps, this is a vector corresponding to the amplitude of each wavelength bin.

property angle_unit: An astropy.units unit defining the angle metric for this map.

property deg

Total degree of the map. Read-only

This is equal to \(y_\mathrm{deg} + u_\mathrm{deg} + f_\mathrm{deg}\).

design_matrix(**kwargs)

Compute and return the light curve design matrix \(A\).

This matrix is used to compute the flux \(f\) from a vector of spherical harmonic coefficients \(y\) and the map amplitude \(a\): \(f = a A y\).

Parameters

xo (scalar or vector, optional) – x coordinate of the occultor relative to this body in units of this body’s radius.
yo (scalar or vector, optional) – y coordinate of the occultor relative to this body in units of this body’s radius.
zo (scalar or vector, optional) – z coordinate of the occultor relative to this body in units of this body’s radius.
ro (scalar, optional) – Radius of the occultor in units of this body’s radius.
theta (scalar or vector, optional) – Angular phase of the body in units of angle_unit.

draw()

Draw a map from the posterior distribution.

This method draws a random map from the posterior distribution and sets the y map vector and amp map amplitude accordingly. Users should call solve() to enable this attribute.

property drorder: Differential rotation order. Read-only

property fdeg: Degree of the multiplicative filter. Read-only

flux(**kwargs)

Compute and return the light curve.

Parameters

xo (scalar or vector, optional) – x coordinate of the occultor relative to this body in units of this body’s radius.
yo (scalar or vector, optional) – y coordinate of the occultor relative to this body in units of this body’s radius.
zo (scalar or vector, optional) – z coordinate of the occultor relative to this body in units of this body’s radius.
ro (scalar, optional) – Radius of the occultor in units of this body’s radius.
theta (scalar or vector, optional) – Angular phase of the body in units of angle_unit.

property inc: The inclination of the rotation axis in units of angle_unit.

intensity(lat=0, lon=0, **kwargs)

Compute and return the intensity of the map.

Parameters

lat (scalar or vector, optional) – latitude at which to evaluate the intensity in units of angle_unit.
lon (scalar or vector, optional) – longitude at which to evaluate the intensity in units of angle_unit`.
theta (scalar, optional) – For differentially rotating maps only, the angular phase at which to evaluate the intensity. Default 0.
limbdarken (bool, optional) – Apply limb darkening (only if udeg > 0)? Default True.

intensity_design_matrix(lat=0, lon=0)

Compute and return the pixelization matrix P.

This matrix is used to compute the intensity \(I\) from a vector of spherical harmonic coefficients \(y\) and the map amplitude \(a\): \(I = a P y\).

Parameters

lat (scalar or vector, optional) – latitude at which to evaluate the design matrix in units of angle_unit.
lon (scalar or vector, optional) – longitude at which to evaluate the design matrix in units of angle_unit.

Note

This method ignores any filters (such as limb darkening or velocity weighting) and illumination (for reflected light maps).

limbdark_is_physical()

Check whether the limb darkening profile (if any) is physical.

This method uses Sturm’s theorem to ensure that the limb darkening intensity is positive everywhere and decreases monotonically toward the limb.

Returns: Whether or not the limb darkening profile is physical.
Return type: bool

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

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

This method computes the marginal likelihood (marginalized over the spherical harmonic coefficients) given a 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.
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.
kwargs (optional) – Keyword arguments to be passed directly to design_matrix(), if a design matrix is not provided.

Returns

The log marginal likelihood, a scalar.

load(image, healpix=False, nside=32, max_iter=3, sigma=None, force_psd=False, **kwargs)

Load an image, array, or healpix map.

This routine uses various routines in healpix to compute the spherical harmonic expansion of the input image and sets the map’s y coefficients accordingly.

Parameters

image – A path to an image file, a two-dimensional numpy array, or a healpix map array (if healpix is True).
healpix (bool, optional) – Treat image as a healpix array? Default is False.
sigma (float, optional) – If not None, apply gaussian smoothing with standard deviation sigma to smooth over spurious ringing features. Smoothing is performed with the healpix.sphtfunc.smoothalm method. Default is None.
force_psd (bool, optional) – Force the map to be positive semi-definite? Default is False.
nside (int, optional) – The NSIDE argument to a Healpix map. This controls the angular resolution of the image; increase this for maps with higher fidelity, at the expense of extra compute time. Default is 32.
max_iter (int, optional) – Maximum number of iterations in trying to convert the map to a Healpix array. Each iteration, the input array is successively doubled in size in order to increase the angular coverage of the input. Default is 3. If the number of iterations exceeds this value, an error will be raised.
kwargs (optional) – Any other kwargs passed directly to minimize() (only if psd is True).

minimize(oversample=1, ntries=1, return_info=False)

Find the global minimum of the map intensity.

Parameters

oversample (int) – Factor by which to oversample the initial grid on which the brute force search is performed. Default 1.
ntries (int) – Number of times the nonlinear minimizer is called. Default 1.
return_info (bool) – Return the info from the minimization call? Default is False.

Returns

A tuple of the latitude, longitude, and the value of the intensity at the minimum. If return_info is True, also returns the detailed solver information.

property nw: Number of wavelength bins. Read-only

property obl: The obliquity of the rotation axis in units of angle_unit.

remove_prior(): Remove the prior on the map coefficients.

render(res=300, projection='ortho', theta=0.0)

Compute and return the intensity of the map on a grid.

Returns an image of shape (res, res), unless theta is a vector, in which case returns an array of shape (nframes, res, res), where nframes is the number of values of theta. However, if this is a spectral map, nframes is the number of wavelength bins and theta must be a scalar.

Parameters

res (int, optional) – The resolution of the map in pixels on a side. Defaults to 300.
projection (string, optional) – The map projection. Accepted values are ortho, corresponding to an orthographic projection (as seen on the sky), and rect, corresponding to an equirectangular latitude-longitude projection. Defaults to ortho.
theta (scalar or vector, optional) – The map rotation phase in units of angle_unit. If this is a vector, an animation is generated. Defaults to 0.0.

rotate(axis, theta)

Rotate the current map vector an angle theta about axis.

Parameters

axis (vector) – The axis about which to rotate the map.
theta (scalar) – The angle of (counter-clockwise) rotation.

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 given a dataset and a prior, provided both are described as multivariate Gaussians.

Parameters

flux (vector) – The observed 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.

set_prior(*, mu=None, L=None, cho_L=None)

Set the prior mean and covariance of the spherical harmonic coefficients.

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

Note that the prior is placed on the amplitude-weighted coefficients, i.e., the quantity x = map.amp * map.y. Because the first spherical harmonic coefficient is fixed at unity, x[0] is the amplitude of the map. The actual spherical harmonic coefficients are given by x / map.amp.

This convention allows one to linearly fit for an arbitrary map normalization at the same time as the spherical harmonic coefficients, while ensuring the starry requirement that the coefficient of the \(Y_{0,0}\) harmonic is always unity.

Parameters

mu (scalar or vector) – The prior mean on the amplitude-weighted spherical harmonic coefficients. Default is unity for the first term and zero for the remaining terms. If this is a vector, it must have length equal to Ny.
L (scalar, vector, or matrix) – The prior covariance. This may be a scalar, in which case the covariance is assumed to be homoscedastic, a vector, in which case the covariance is assumed to be diagonal, or a matrix specifying the full prior covariance. Default is None. Either L or cho_L must be provided.
cho_L (matrix) – The lower Cholesky factorization of the prior covariance matrix. Defaults to None. Either L or cho_L must be provided.

show(**kwargs)

Display an image of the map, with optional animation. See the docstring of render() for more details and additional keywords accepted by this method.

Parameters

ax (optional) – A matplotlib axis instance to use. Default is to create a new figure.
cmap (string or colormap instance, optional) – The matplotlib colormap to use. Defaults to plasma.
figsize (tuple, optional) – Figure size in inches. Default is (3, 3) for orthographic maps and (7, 3.5) for rectangular maps.
projection (string, optional) – The map projection. Accepted values are ortho, corresponding to an orthographic projection (as seen on the sky), and rect, corresponding to an equirectangular latitude-longitude projection. Defaults to ortho.
colorbar (bool, optional) – Display a colorbar? Default is False.
grid (bool, optional) – Show latitude/longitude grid lines? Defaults to True.
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 figure or animation to. 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.)
dpi (int, optional) – Image resolution in dots per square inch. Defaults to the value defined in matplotlib.rcParams.
bitrate (int, optional) – Bitrate in kbps (animations only). Defaults to the value defined in matplotlib.rcParams.
norm (optional) – The color normalization passed to matplotlib.pyplot.imshow. Default is None.

Note

Pure limb-darkened maps do not accept a projection keyword.

property solution

The posterior probability distribution for the map.

This is a tuple containing the mean and lower Cholesky factorization of the covariance of the amplitude-weighted spherical harmonic coefficient vector, 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 vector. Under this convention, the map amplitude is equal to the first term of the vector and the spherical harmonic coefficients are equal to the vector normalized by the first term.

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

Solve the linear least-squares problem for the posterior over maps.

This method solves the generalized least squares problem given a 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 amplitude and coefficients are 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.
kwargs (optional) – Keyword arguments to be passed directly to design_matrix(), if a design matrix is not provided.

Returns

A tuple containing the posterior mean for the amplitude-weighted spherical harmonic coefficients (a vector) and the Cholesky factorization of the posterior covariance (a lower triangular matrix).

Note

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

property u

The vector of limb darkening coefficients. Read-only

To set this vector, index the map directly using one index: map[n] = ... where n is the degree of the limb darkening coefficient. This may be an integer or an array of integers. Slice notation may also be used.

property udeg: Limb darkening degree. Read-only

property y

The spherical harmonic coefficient vector. Read-only

To set this vector, index the map directly using two indices: map[l, m] = ... where l is the spherical harmonic degree and m is the spherical harmonic order. These may be integers or arrays of integers. Slice notation may also be used.

property ydeg: Spherical harmonic degree of the map. Read-only