Notebook

Orbit visualization¶

In this notebook, we'll visualize a simple star-planet system and its light curve.

In [ ]:

%matplotlib inline

In [ ]:

%run notebook_setup.py

In [ ]:

from matplotlib.animation import FuncAnimation
from IPython.display import HTML

First we import stuff. Note that we disable lazy evaluation in starry: all map attributes and method return values will be actual numpy floats and arrays.

In [ ]:

import matplotlib.pyplot as plt
import numpy as np
import starry

starry.config.lazy = False
starry.config.quiet = True

Let's instantiate a star. We'll give it a bit of quadratic limb darkening.

In [ ]:

A = starry.Primary(starry.Map(udeg=2, amp=1.0), r=1.0, m=1.0)
A.map[1:] = [0.5, 0.25]

Now we instantiate a planet and give it the 10th degree spherical harmonic expansion of the map of the Earth. The values below aren't at all realistic, but they'll make for a cool visualization below. By default, mass and radius units are solar, angles are in degrees, and times are in days. Finally, we set the map inclination and the orbital inclination to the same value, and likewise the map obliquity and longitude of ascending node. This causes the axis of rotation to be perpendicular to the orbital plane (i.e., the orbital and rotational angular momentum vectors are parallel to each other, as we'd expect for a tidally-locked planet).

In [ ]:

b = starry.Secondary(
    starry.Map(ydeg=10, inc=80.0, obl=30.0, amp=0.1),
    r=0.5,
    m=0.5,
    porb=1.0,
    prot=1.0,
    t0=0.0,
    inc=80.0,
    Omega=30.0,
)
b.map.load("earth")

Finally, we instantiate a Keplerian system:

In [ ]:

sys = starry.System(A, b)

Compute the positions of the two bodies over the course of one orbit of the planet:

In [ ]:

npts = 200
t = np.linspace(-0.5, 0.5, npts)
x, y, z = sys.position(t)

Render the maps over the same time period:

In [ ]:

res = 300
theta_sec = [360.0 / sec.prot * (t - sec.t0) - sec.theta0 for sec in sys.secondaries]
img = np.array(
    [np.tile(sys.primary.map.render(res=300), (npts, 1, 1))]
    + [
        sec.map.render(theta=theta_sec[i], res=res)
        for i, sec in enumerate(sys.secondaries)
    ]
)

Compute the full system light curve:

In [ ]:

flux = sys.flux(t)

Let's visualize everything. We can normally visualize the orbit by calling

sys.show(t)

but here's a more souped-up version just for fun:

In [ ]:

# Set up the plot
fig, ax = plt.subplots(1, figsize=(6.5, 7))
ax_xz = fig.add_axes([0.275, 0.8, 0.2, 0.2])
ax_xz.annotate(
    "Top", fontsize=12, xy=(0, 0), xycoords="axes fraction", ha="left", va="bottom"
)
ax_zy = fig.add_axes([0.525, 0.8, 0.2, 0.2])
ax_zy.annotate(
    "Side", fontsize=12, xy=(0, 0), xycoords="axes fraction", ha="left", va="bottom"
)
ax_lc = fig.add_axes([0.125, 0.05, 0.775, 0.2])

xz = [None] + [None for sec in sys.secondaries]
xy = [None] + [None for sec in sys.secondaries]
zy = [None] + [None for sec in sys.secondaries]
circ = [None] + [None for sec in sys.secondaries]
maps = [sys.primary.map] + [sec.map for sec in sys.secondaries]
radii = np.array([sys.primary.r] + [sec.r for sec in sys.secondaries])

for axis, arrs in zip([ax, ax_xz, ax_zy], [(x, y), (x, z), (z, y)]):
    axis.axis("off")
    R = 1.2 * max(-np.min(arrs), np.max(arrs))
    axis.set_xlim(-R, R)
    axis.set_ylim(-R, R)

# Plot the light curve
ax_lc.plot(t, flux, "k-")
(lc,) = ax_lc.plot(t[0], flux[0], "o", color="k")
ax_lc.axis("off")

# Plot the first frame
for i, xi, yi, zi, map, r in zip(range(1 + len(sys.secondaries)), x, y, z, maps, radii):

    # Orbit outlines
    ax_xz.plot(xi, zi)
    ax_zy.plot(zi, yi)

    # Body positions
    xz[i] = ax_xz.scatter(xi[0], zi[0])
    zy[i] = ax_zy.scatter(zi[0], yi[0])

    # Maps
    extent = np.array([xi[0], xi[0], yi[0], yi[0]]) + np.array([-1, 1, -1, 1]) * r
    xy[i] = ax.imshow(
        img[i, 0],
        origin="lower",
        cmap="plasma",
        extent=extent,
        clip_on=False,
        zorder=zi[0],
    )
    circ[i] = plt.Circle(
        (xi[0], yi[0]), r, color="k", fill=False, zorder=zi[0] + 1e-3, lw=3
    )
    ax.add_artist(circ[i])

# Animation
def updatefig(k):
    for i, xi, yi, zi, map, r in zip(
        range(1 + len(sys.secondaries)), x, y, z, maps, radii
    ):
        xz[i].set_offsets((xi[k], zi[k]))
        zy[i].set_offsets((zi[k], yi[k]))
        xy[i].set_extent(
            np.array([xi[k], xi[k], yi[k], yi[k]]) + np.array([-1, 1, -1, 1]) * r
        )
        xy[i].set_zorder(zi[k])
        xy[i].set_data(img[i, k])
        circ[i].center = (xi[k], yi[k])
        circ[i].set_zorder(zi[k] + 1e-3)
    lc.set_xdata(t[k])
    lc.set_ydata(flux[k])
    return xz + xy + zy + circ + [lc]


ani = FuncAnimation(fig, updatefig, interval=30, blit=False, frames=len(t))
plt.close()
display(HTML(ani.to_html5_video()))