Fitting Planets

Although Starlord is designed for fitting stellar models to observations, it is perfectly capable of fitting planet evolution models instead. If anything, planet models are simpler as the mass and/or radius are usually measured directly, rather than inferred from a collection of observed bands. Existing model grids can be readily input into Starlord by the user, allowing modellers to avoid the laborious process of retrieval setup.

For illustrative purposes we’ll use the hot Jupiter evolution model grids of Thorngren & Fortney 2018; these will be made publicly available when the Starlord built-in grid repository is setup, but for now may be obtained by contacting the developers. This grid is 5 dimensional, in mass, metallicity, incident flux, internal heating power (as a fraction of flux), and age. The outputs are the specific entropy, luminosity, and radius.

Defining the Model Grid

We’ll start by showing how to convert your planet grid into the Starlord format (required to fit models to it). This is also covered in Grid Management, but will be so common for planet modelling as to merit a tailored example. This mostly consists of naming the axes, defining any derived parameters, and handing the data to starlord.GridGenerator.create_grid() for processing. For this example, we’ll assume the data was stored in a csv file and open it with Pandas (in reality it wasn’t but this is an important case to cover).

from collections import OrderedDict
import numpy as np
import pandas as pd
import starlord

grid = pd.read_csv("hotJupiters.csv", skipinitialspace=True, sep=",")
print(grid.columns)

# Rearrange the grid into 3 outputs in 5 dimensions
inputs, outputs = starlord.GridGenerator.restructure_grid(grid.values, (0, 1, 2, 3, 4), (5, 6, 7))
(mass, zpl, flux, heating, age) = inputs
radius, entropy, luminosity = outputs

# The axes of the grids -- order matters!
# Mass, flux, and age span multiple order of magnitude so it's best to interpolate in logspace
inputs = OrderedDict(
    log_mass=np.log10(mass),
    zpl=zpl,
    log_flux=np.log10(flux),
    heating=heating,
    log_age=np.log10(age),
)

# The output values to interpolate on, order does not matter anymore.
outputs = dict(
    log_radius=np.log10(radius),
    entropy=entropy,
    log_luminosity=np.log10(luminosity),
)

# Optionally, derived values that can be calculated from the outputs or inputs
derived = dict(
    age="10**d.hotJupiters.log_age",
    flux="10**d.hotJupiters.log_flux",
    luminosity="10**d.hotJupiters.log_luminosity",
    mass="10**d.hotJupiters.log_mass",
    radius="10**d.hotJupiters.log_radius",
    tint="math.pow(d.hotJupiters.luminosity / (7.125593e-4 * (d.hotJupiters.radius * 6.9911e9)**2), 0.25)",
    typical_heating="0.0237 * math.exp(-(d.hotJupiters.log_flux - 9.14)**2 / (2 * .37**2))",
)

# Input everything into the
starlord.GridGenerator.create_grid(
    "hotJupiters",
    inputs=inputs,
    outputs=outputs,
    derived=derived,
)

The steps here are:

Load your data from whatever format it is saved in.
Rearrange it so that the input axes are 1d and sorted and the outputs are nd-arrays arranged so that each axis has the same length as the corresponding input axis.
Define the input names with an OrderedDict, optionally transforming axes (e.g. taking the log) for better interpolation.
Name the output axes with a dict, again optionally transforming them.
If desired, define formulas for derived parameters as Cython code, using math for math operations and d.gridname.input_or_output to refer to other grid data.
Dump all of this along with a grid name into starlord.GridGenerator.create_grid().

If all went well, you should see your new grid listed when you run starlord -g. You would be wise to run starlord -g gridname and double-check that the axes, minima, maxima, and lengths all make sense to you – it’s very easy to make unit errors here. For this example, you would see:

Grid hotJupiters
    Input                       Min        Max     Length     Default Mapping
log_mass                 -1.326      1.301         30     p.log_mass--i
zpl                           0          1         10     p.zpl--i
log_flux                      6         11         30     p.log_flux--i
heating                       0       0.05         30     p.heating--i
log_age                      -1      1.176         12     p.log_age--i
Outputs
    Output                      Min        Max
entropy                   6.148         14
log_luminosity           -42.73      32.65
log_radius              -0.6131      2.109
Derived
    Derived              Code
age                  10**d.hotJupiters.log_age
flux                 10**d.hotJupiters.log_flux
luminosity           10**d.hotJupiters.log_luminosity
mass                 10**d.hotJupiters.log_mass
radius               10**d.hotJupiters.log_radius
tint                 math.pow(d.hotJupiters.luminosity / (7.125593e-4 * (d.hotJupiters.radius * 6.991 ...
typical_heating      0.0237 * math.exp(-(d.hotJupiters.log_flux - 9.14)**2 / (2 * .37**2))

Specifying the Model

To actually fit a planet using this grid, you’ll need to define a model. This can be done with the Python API, but we’ll use the toml file approach here (see Defining a Model). You should start by writing down your likelihood terms as gridname.output = ["distribution", param1, param2,...]. For a transiting planet this will generally include the radius and mass.

[model]
# Likelihood terms (both normal distributions, in this case)
hotJupiters.radius = [1.5, 0.05]
hotJupiters.mass = [0.9, 0.1]

You can already run starlord -da planet.toml to get a sense of how Starlord will interpret your model. Pay close attention to the Params: p.log_age, p.log_mass, ... line; you need to either set priors for every parameter or change the model to not use it by e.g. fixing them. Below, we set log_flux to a constant value and fix heating to a formula defined in the grid (see above) which is itself a function of flux (you could also have directly input the formula on this line). Finally we set the priors in the same manner as the likelihood terms were set.

# Fixing the log_flux to a constant
override.hotJupiters.log_flux = "9.2"
# Fixing the heating to a function of flux defined by the grid
override.hotJupiters.heating = "d.hotJupiters.typical_heating"

# Mass and age priors set from observations
prior.log_mass = ['uniform', -0.5, 0.5]
prior.log_age = ['uniform', -1.0, 1.0]
# Uninformative Z_planet prior
prior.zpl = ['uniform', 0.0, 1.0]

It can be helpful for Starlord to output additional derived values like the intrinsic temperature even though they aren’t parameters of the model. This is set in the model outputs section as shown. Finally, sampling and output options are set in the [sampling] and [output] sections respectively – we’ll just use a simple setup here, see Sampling for more information.

# Additional outputs for the sampler to record
# (Mass and age only because the sampler records only their log10 otherwise)
outputs = [
    "d.hotJupiters.mass",
    'd.hotJupiters.age',
    'd.hotJupiters.entropy',
    'd.hotJupiters.tint',
]

[sampling]
sampler="emcee"
emcee_run.progress = true

[output]
terminal=true
file="hotJupiter.npz"

Running and Reading Outputs

Running the model with starlord planet.toml, we obtain:

100%|███████████████████████████████████████████████| 5500/5500 [00:01<00:00, 3761.72it/s]
Convergence: Tau = 50.18; N/Tau = 99.65
     Name                            Mean         Std         16%         50%         84%
   0 log_age                   -0.0003021      0.5785     -0.6835   0.0002408      0.6796
   1 log_mass                    -0.05069     0.04869    -0.09833    -0.04862    -0.00254
   2 zpl                            0.125     0.03176     0.09352      0.1246      0.1566
-----------------------------------------------------------------------------------------
   3 log_like                       2.471       0.997       1.644       2.779       3.292
   4 log_prior                    -0.6931   3.967e-13     -0.6931     -0.6931     -0.6931
   5 hotJupiters__mass             0.8954     0.09903      0.7974      0.8941      0.9942
   6 hotJupiters__age                2.15       2.495      0.2072       1.001       4.782
   7 hotJupiters__entropy           9.816     0.04114       9.776       9.814       9.857
   8 hotJupiters__tint              636.1     0.09286         636       636.1       636.1

The three model parameters are listed first, then the log_likelihood, log_prior, and output values requested in the toml file. We can see, for example, that this planet was inferred to have a metallicity of 0.125 +/- 0.032, and an intrinsic temperature of 636 K. Because we fixed the flux and heating, the uncertainties on the latter are very small – relaxing that assumption would get us more realistic uncertainties.

In the model file we also specified an output file of hotJupiter.npz – this was saved in the directory the model was run in. The data can be loaded in using np.load or with starlord.load_to_frame() to obtain a nicely-formatted Pandas data frame of the posterior.