Quick Start (@id quick-start)
This guide introduces the key concepts in Octofitter:
- Likelihood objects to hold your data
@planet
and@system
models to specify variables, priors, and system architecture- Sampling from the posterior using MCMC
- Plotting the results
- Saving the chain
For installation instructions, see Installation.
Example: Fit a Single Planet Orbit
Load the required packages:
using Octofitter, Distributions, CairoMakie, PairPlots
Create a PlanetRelAstromLikelihood
object containing your observational data. In this case its the position of the planet relative to the star, but many other kinds of data are supported:
astrom = PlanetRelAstromLikelihood(Table(
epoch = [50000, 50120, 50240], # Dates in MJD
ra = [-505.7, -502.5, -498.2], # [mas] East positive
dec = [-66.9, -37.4, -7.9], # [mas] North positive
σ_ra = [10.0, 10.0, 10.0], # [mas] Uncertainties
σ_dec = [10.0, 10.0, 10.0], # [mas] Uncertainties
cor = [0.0, 0.0, 0.0] # RA/Dec correlations
))
Define a planet model with orbital elements and their prior distributions:
@planet b Visual{KepOrbit} begin
a ~ Uniform(0, 100) # Semi-major axis [AU]
e ~ Uniform(0.0, 0.5) # Eccentricity
i ~ Sine() # Inclination [rad]
ω ~ UniformCircular() # Argument of periastron [rad]
Ω ~ UniformCircular() # Longitude of ascending node [rad]
θ ~ UniformCircular() # Position angle at reference epoch [rad]
# Epoch of periastron passage
# We calculate it from the position angle above
tp = θ_at_epoch_to_tperi(system,b,50000)
end astrom
Make sure to adjust the epoch 50000
above to match your most constraining data epoch.
Define the system with its mass and distance - see System Construction for more options:
@system HD1234 begin
M ~ truncated(Normal(1.2, 0.1), lower=0.1) # Total mass (solar masses)
plx ~ truncated(Normal(50.0, 0.02), lower=0.1) # Parallax (mas)
end b
That there are many different orbit parameterizations, each requiring different of parameters names. The KepOrbit
is a full 3D keplerian orbit with Campbell parameters. Visual
means that we have a defined parallax distance plx
that can map separations in AU to arcseconds.
Compile the model into efficient sampling code:
model = Octofitter.LogDensityModel(HD1234)
Sample from the posterior using Hamiltonian Monte Carlo (see Samplers for other options):
chain = octofit(model)
Visualize the results with orbit plots and a corner plot:
octoplot(model, chain) # Plot orbits and data
octocorner(model, chain) # Corner plot of posterior
Save the results to a FITS file (see Loading and Saving Data for other formats):
savechain("output.fits", chain)
Working with Dates
Convert dates to and from Modified Julian Days using these helper functions:
mjd("2020-01-01") # Date string to MJD
years2mjd(2020.0) # Decimal year to MJD
mjd2date(50000) # MJD to date
Next Steps
See the Tutorials section for complete examples.