Note
This tutorial was generated from a Jupyter notebook that can be downloaded here.
Star spots¶
A major part of the philosophy of starry
is a certain amount of agnosticism about what the surface of a star or planet actually looks like. Many codes fit stellar light curves by solving for the number, location, size, and contrast of star spots. This is usually fine if you know the stellar surface consists of a certain number of discrete star spots of a given shape. In general, however, that’s a very strong prior to assume. And even in cases where it holds, the problem is still extremely
degenerate and lies in a space that is quite difficult to sample.
Instead, in starry
we assume the surface is described by a vector of spherical harmonic coefficients. The advantage of this is that (1) it automatically bakes in a Gaussianprocess smoothness prior over the surface map, in which the scale of features is dictated by the degree of the expansion; and (2) under gaussian priors and gaussian errors, the posterior over surface maps is analytic. If and only if the data and prior support the existence of discrete star spots on the surface, the
posterior will reflect that.
However, sometimes it’s convenient to restrict ourselves to the case of discrete star spots. In starry
, we therefore implement the spot
method, which adds a spotlike feature to the surface by expanding a circular tophat in terms of spherical harmonics.
Note
This method replaced the add_spot
method introduced in version 1.0.0
.
In this notebook, we’ll take a look at how this new method works. For reference, here is the docstring of starry.Map.spot
:

spot
( contrast=1.0, radius=None, lat=0.0, lon=0.0 ) ¶ 
Add the expansion of a circular spot to the map.
This function adds a spot whose functional form is a top hat in \(\Delta\theta\), the angular separation between the center of the spot and another point on the surface. The spot intensity is controlled by the parameter
contrast
, defined as the fractional change in the intensity at the center of the spot. Parameters


contrast (scalar or vector, optional) – The contrast of the spot. This is equal to the fractional change in the intensity of the map at the center of the spot relative to the baseline intensity of an unspotted map. If the map has more than one wavelength bin, this must be a vector of length equal to the number of wavelength bins. Positive values of the contrast result in dark spots; negative values result in bright spots. Default is
1.0
, corresponding to a spot with central intensity close to zero. 
radius (scalar, optional) – The angular radius of the spot in units of
angle_unit
Defaults to20.0
degrees. 
lat (scalar, optional) – The latitude of the spot in units of
angle_unit
Defaults to20.0
degrees. 
lon (scalar, optional) – The longitude of the spot in units of
angle_unit
Defaults to20.0
degrees.

Note
Keep in mind that things are normalized in starry
such that
the diskintegrated flux (not the intensity!)
of an unspotted body is unity. The default intensity of an
unspotted map is 1.0 / np.pi
everywhere (this ensures the
integral over the unit disk is unity).
So when you instantiate a map and add a spot of contrast c
,
you’ll see that the intensity at the center is actually
(1  c) / np.pi
. This is expected behavior, since that’s
a factor of 1  c
smaller than the baseline intensity.
Note
This function computes the spherical harmonic expansion of a
circular spot with uniform contrast. At finite spherical
harmonic degree, this will return an approximation that
may be subject to ringing. Users can control the amount of
ringing and the smoothness of the spot profile (see below).
In general, however, at a given spherical harmonic degree
ydeg
, there is always minimum spot radius that can be
modeled well. For ydeg = 15
, for instance, that radius
is about 10
degrees. Attempting to add a spot smaller
than this will in general result in a large amount of ringing and
a smaller contrast than desired.

There are a few additional underthehood keywords that control the behavior of the spot expansion. These are
 Parameters


spot_pts (int, optional) – The number of points in the expansion of the (1dimensional) spot profile. Default is
1000
. 
spot_eps (float, optional) – Regularization parameter in the expansion. Default is
1e9
. 
spot_smoothing (float, optional) – Standard deviation of the Gaussian smoothing applied to the spot to suppress ringing (unitless). Default is
2.0 / self.ydeg
. 
spot_fac (float, optional) – Parameter controlling the smoothness of the spot profile. Increasing this parameter increases the steepness of the profile (which approaches a top hat as
spot_fac > inf
). Decreasing it results in a smoother sigmoidal function. Default is300
. Changing this parameter is not recommended; change spot_smoothing instead.

Note
These last four parameters are cached. That means that
changing their value in a call to spot
will result in
all future calls to spot
“remembering” those settings,
unless you change them back!
Adding a simple spot¶
Let’s begin by importing stuff as usual:
[3]:
import numpy as np
import starry
starry.config.lazy = False
starry.config.quiet = True
The first thing we’ll do is create a dummy featureless map, which we’ll use for comparisons below.
[4]:
map0 = starry.Map(ydeg=1)
map0.show()
Now let’s instantiate a very high degree map and add a spot with a contrast of \(25\%\) and a radius of \(15^\circ\) at latitude/longitude \((0, 0)\):
[5]:
contrast = 0.25
radius = 15
map = starry.Map(ydeg=30)
map.spot(contrast=contrast, radius=radius)
Here’s what we get:
[6]:
map.show(theta=np.linspace(0, 360, 50))
The spot contrast¶
Let’s look at what a contrast of \(25\%\) implies for the actual intensity of points on the surface of our map. Fore definiteness, let’s look at the intensity along the equator:
[7]:
lon = np.linspace(90, 90, 100)
plt.plot(lon, map.intensity(lon=lon))
plt.xlabel("longitude [deg]")
plt.ylabel("intensity");
It’s worth recalling how intensities are normalized in starry
(see the note in the docstring above). The baseline intensity of an unspotted map is 1.0 / np.pi
, so the spot intensity is a \(25\%\) reduction relative to that value. Let’s normalize the function above to the continuum level to see that:
[8]:
plt.plot(lon, np.pi * map.intensity(lon=lon))
plt.xlabel("longitude [deg]")
plt.ylabel("normalized intensity");
As expected, the intensity of the spot is \(25\%\) lower than the baseline.
The spot expansion¶
As mentioned in the docstring, the spot is modeled as a top hat in \(\Delta\theta\), the angular distance along the surface of the body. The profile in the figure above doesn’t really look like a top hat\(\)that’s because we’re actually expanding the functional form of a top hat in the spherical harmonic basis. Since we’re at finite spherical harmonic degree, the expansion is not exact. In particular, the spherical harmonic basis enforces smoothing at small scales. You get the same behavior as if you tried to fit a 1d top hat with a loworder polynomial.
Because of this, there are a few things we must keep in mind.
First, the function we’re actually expanding is a sigmoid
with smoothing parameter \(c\). In the limit that \(c \rightarrow \infty\), \(f(\Delta\theta)\) approaches an actual top hat. In practice, this can be problematic because of the discontinuity at \(\Delta\theta = r\), where the function is not differentiable. Setting \(c\) to a large but finite value makes the function continuous and differentiable everywhere.
Here’s the function we’re expanding for different values of \(c\):
[9]:
radius = 15
dtheta = np.abs(lon)  radius
for c in [0.25, 0.5, 1, 10, 300]:
plt.plot(lon, 1 / (1 + np.exp(c * dtheta)), label=c)
plt.legend(title="c", fontsize=10)
plt.xlabel("angular distance [deg]")
plt.ylabel("intensity");
In starry
, the smoothing parameter \(c\) is specified as the keyword spot_fac
(see docs above). By default it’s set to 300
(purple line), so it’s actually almost indistinguishable from a true top hat.
OK, but even though the function we’re expanding is almost exactly a top hat, the actual intensity on the surface of the sphere is much smoother. That’s because of the finite spherical harmonic expansion, as mentioned above. Let’s compare the input function and the actual intensity profile after the expansion:
[10]:
# Plot the starry intensity along the equator
map.reset()
map.spot(contrast=1, radius=radius)
plt.plot(lon, np.pi * map.intensity(lon=lon), label="measured")
# Plot the actual function we're expanding
c = 300
dtheta = np.abs(lon)  radius
plt.plot(lon, 1 / (1 + np.exp(c * dtheta)), label="exact")
plt.legend()
plt.xlabel("longitude [deg]")
plt.ylabel("intensity");
There’s actually another level of smoothing happening in starry
that is controlled by the spot_smoothing
parameter. Let’s look at what the measured intensity profile would be if we got rid of smoothing altogether:
[11]:
# Plot the starry intensity along the equator
map.reset()
map.spot(contrast=1, radius=radius, spot_smoothing=0)
plt.plot(lon, np.pi * map.intensity(lon=lon), label="measured (no smoothing)")
# Plot the actual function we're expanding
c = 300
dtheta = np.abs(lon)  radius
plt.plot(lon, 1 / (1 + np.exp(c * dtheta)), label="exact")
plt.legend()
plt.xlabel("longitude [deg]")
plt.ylabel("intensity");
The measured intensity is steeper and closer to the exact profile at \(\Delta\theta = r\), but we’re paying a high price for this with all the ringing! This is the classical Gibbs phenomenon, which you’ve seen if you’ve ever tried to compute the Fourier expansion (or transform) of a square wave. At finite degree, the spherical harmonic basis simply can’t capture sharp discontinuities. The spot_smoothing
parameter convolves the surface
map with a function of the form
where \(\sigma =\) spot_smoothing
. By default, this is equal to 2 / self.ydeg
, which results in the smooth (but gentler) spot profile in the previous figure. Increase this to suppress ringing, at the expense of a smoother profile; decrease it for the opposite effect.
The final point we should address concerns the minimum radius
we can model at a given spherical harmonic degree. Let’s visualize the measured spot profile for different radii ranging from \(3^\circ\) to \(25^\circ\) (using the original, smoothed map
):
[12]:
lon = np.linspace(45, 45, 300)
radii = np.arange(3, 26)
cmap = plt.get_cmap("plasma_r")
for k, radius in enumerate(radii):
map.reset()
map.spot(contrast=1, radius=radius, spot_smoothing=2 / map.ydeg)
plt.plot(
lon, np.pi * map.intensity(lon=lon), color=cmap(0.25 + 0.5 * k / len(radii))
)
if radius <= 10:
plt.text(
0, 0.01 + np.pi * map.intensity(lon=0), radius, ha="center", fontsize=6
)
if radius == 25:
plt.text(
0, 0.05 + np.pi * map.intensity(lon=0), radius, ha="center", fontsize=6
)
plt.xlabel("longitude [deg]")
plt.ylabel("intensity");
For \(r \ge 10^\circ\), the spherical harmonic expansion has no trouble in capturing the spot profile. But for smaller radii, the profile doesn’t actually get any smaller\(\)instead, its amplitude starts to decrease! That’s because a spherical harmonic expansion of degree ydeg
has a minimum angular scale it can model\(\)below that scale, things get smoothed away (and their amplitude gets washed out). The angular scale is proportional to 1/ydeg
. As a very rough rule
of thumb, it scales as 180 / ydeg
. That would mean our minimum diameter is \(6^\circ\) (or a radius of \(3^\circ\), but that doesn’t account for the fact that we are adding additional smoothing to suppress ringing. That smoothing makes the minimum radius closer to \(10^\circ\) for ydeg=30
.
If we weren’t adding any smoothing, then we can see that the limiting radius is in fact very close to \(3^\circ\):
[13]:
lon = np.linspace(45, 45, 300)
radii = np.arange(3, 26)
cmap = plt.get_cmap("plasma_r")
for k, radius in enumerate(radii):
map.reset()
map.spot(contrast=1, radius=radius, spot_smoothing=0)
plt.plot(
lon, np.pi * map.intensity(lon=lon), color=cmap(0.25 + 0.5 * k / len(radii))
)
if radius == 3:
plt.text(
0, 0.01 + np.pi * map.intensity(lon=0), radius, ha="center", fontsize=6
)
plt.xlabel("longitude [deg]")
plt.ylabel("intensity");
But in many cases, that may be an unacceptable level of ringing.
Long story short: there’s always going to be a tradeoff between the amount of ringing and the maximum resolution of features you can model at a given spherical harmonic degree.
The last thing we’ll do in this tutorial is compute the minimum spot radius we can model at a given spherical harmonic degree and with a given smoothing strength. We’ll define this minimum value as the radius below which we get more than \(10\%\) error in the mean contrast of the spot. You can change this tolerance by tweaking the tol
parameter below:
[14]:
smstrs = [0.0, 1.0, 2.0]
ydegs = [10, 15, 20, 25, 30]
radii = np.arange(1, 90)
tol = 0.1
rmin = np.zeros((len(smstrs), len(ydegs)))
for i, smstr in enumerate(smstrs):
for j, ydeg in enumerate(ydegs):
map = starry.Map(ydeg)
error = np.zeros(len(radii))
for radius in radii:
map.reset()
map.spot(contrast=1, radius=radius, spot_smoothing=smstr / ydeg)
if (
np.abs(np.mean(map.intensity(lon=np.linspace(0, 0.75 * radius, 50))))
< tol
):
rmin[i, j] = radius
break
plt.figure(figsize=(8, 8))
plt.plot(ydegs, rmin[0], "o", label=r"$0$")
plt.plot(ydegs, rmin[1], "o", label=r"$1 / l_\mathrm{max}$")
plt.plot(ydegs, rmin[2], "o", label=r"$2 / l_\mathrm{max}$ (default)")
plt.legend(title="smoothing", fontsize=10)
plt.xticks([10, 15, 20, 25, 30])
plt.xlabel(r"spherical harmonic degree $l_\mathrm{max}$")
plt.ylabel(r"minimum spot radius [degrees]");
The green line shows the default smoothing strength (the spot_smoothing
parameter); the orange line shows half that amount of smoothing, and the blue line shows no smoothing. In all cases, the minimum radius we can safely model decreases roughly as \(1/l_\mathrm{max}\).
Use this figure as a rule of thumb when modeling discrete spots with starry
.