API Documentation

@planet [planet_name] [Orbit_Type] begin
    [prior_1] ~ [UnivariateDistribution]
    [prior_2] ~ [UnivariateDistribution]
    calculation_3 = [planet_name].[prior_1] + [planet_name].[prior_2] + [system.variable_x]
end [likelihood_objects...]

Generate a Planet model named planet_name. A variable will be created with the name [planet_name] in the current scope. Orbit_Type specifies the orbit parameterization from PlanetOrbits.jl. You must provide all input variables needed for the selected orbit type (see PlanetOrbits.jl documentation). Following that is a block of variable assignments. Variables with a ~ will be free variables with a prior distribution given by the right-hand-side (a UnivariateDistribution from Distributions.jl or a KDEDist). Calculated quantities are also allowed. These may reference other variables using the planet name followed by a dot and the name of the variable. Variables from other planets in a single system are not accessible. You can access other variables in the current local scope, but these bindings are only guaranteed to be resolved a single time. Note that using non-constant global variables in calculated expressions can lead to poor performance. Finally, the planet model can be conditioned on data by supplying zero or more likelihood objects.

@system [system_name] begin
    [prior_1] ~ [UnivariateDistribution]
    [prior_2] ~ [UnivariateDistribution]
    calculation_3 = system.[prior_1] + system.[prior_2]
end [likelihood_objects...] [planet_models...]

Generate a System model named system_name. A variable will be created with the name [system_name] in the current scope. Following that is a block of variable assignments. Variables with a ~ will be free variables with a prior distribution given by the right-hand-side (a UnivariateDistribution from Distributions.jl or a KDEDist). Calculated quantities are also allowed. These may reference other variables using the planet name followed by a dot and the name of the variable. Variables from other planets in a single system are not accessible. You can access other variables in the current local scope, but these bindings are only guaranteed to be resolved a single time. Note that using non-constant global variables in calculated expressions can lead to poor performance. After the end of the variable block, the system model can be conditioned on data by supplying zero or more likelihood objects. Finally, zero or more planet models can be attached to the system, potentially conditioned on likelihood objects of their own.

PlanetRelAstromLikelihood(
    (epoch = 5000, ra = -505.7637580573554, dec = -66.92982418533026, σ_ra = 10, σ_dec = 10, cor=0),
    (epoch = 5050, ra = -505.7637580573554, dec = -66.92982418533026, σ_ra = 10, σ_dec = 10, cor=0),
    (epoch = 5100, ra = -505.7637580573554, dec = -66.92982418533026, σ_ra = 10, σ_dec = 10, cor=0),
)

Represents a likelihood function of relative astometry between a host star and a secondary body. :epoch is a required column, in addition to either :ra, :dec, :σ_ra, :σ_dec or :pa, :sep, :σ_pa, :σ_sep. All units are in milliarcseconds or radians as appropriate.

In addition to the example above, any Tables.jl compatible source can be provided.

StarAbsoluteRVLikelihood(
    (;inst_idx=1, epoch=5000.0,  rv=−6.54, σ_rv=1.30),
    (;inst_idx=1, epoch=5050.1,  rv=−3.33, σ_rv=1.09),
    (;inst_idx=1, epoch=5100.2,  rv=7.90,  σ_rv=.11),
)

Represents a likelihood function of relative astometry between a host star and a secondary body. :epoch (mjd), :rv (m/s), and :σ_rv (m/s), and :inst_idx are all required.

:inst_idx is used to distinguish RV time series between insturments so that they may optionally be fit with different zero points and jitters. In addition to the example above, any Tables.jl compatible source can be provided.

PlanetRelativeRVLikelihood(
    (;inst_idx=1, epoch=5000.0,  rv=−6.54, σ_rv=1.30),
    (;inst_idx=1, epoch=5050.1,  rv=−3.33, σ_rv=1.09),
    (;inst_idx=1, epoch=5100.2,  rv=7.90,  σ_rv=.11),
)

Represents a likelihood function of relative astometry between a host star and a secondary body. :epoch (mjd), :rv (m/s), and :σ_rv (m/s), and :inst_idx are all required.

:inst_idx is used to distinguish RV time series between insturments so that they may optionally be fit with different zero points and jitters. In addition to the example above, any Tables.jl compatible source can be provided.

PhotometryLikelihood(
    (band = :Z, phot=15.0, σ_phot=3.),
    (band = :J, phot=13.5, σ_phot=0.5),
    (band = :L, phot=11.0, σ_phot=1.0)
)

A likelihood for comparing measured photometry points in one or more filter bands to data (provided here). Requires the :band, :phot', and:σ_phot` columns. Can be provided with any Tables.jl compatible data source.

HGCALikelihood(;gaia_id=1234)

Load proper motion anomaly data from the HIPPARCOS-GAIA Catalog of Accelerations (Brandt et al) for a star with catalog id gaia_id. The resulting velocities are in mas/yr and have the long term trend between HIPPARCOS and GAIA already subtracted out. e.g. we would expect 0 pma if there is no companion.

ObsPriorAstromONeil2019(astrometry_likelihood, period_prior)

Given a an astrometry likelihood (PlanetRelAstromLikelihood), apply the "observable based priors" of K. O'Neil 2019 "Improving Orbit Estimates for Incomplete Orbits with a New Approach to Priors: with Applications from Black Holes to Planets".

This prior correction is only correct if you supply Uniform priors on all Campbell orbital parameters and a Uniform prior on Period (not semi-major axis). This period prior has a significant impact in the fit and recommendations for its range were not published in the original paper.

Examples

astrom_like = PlanetRelAstromLikelihood(astrom_table)

# Apply observable based priors ontop of our uniform Campbell priors:
obs_prior = ObsPriorAstromONeil2019(astrom_like)

# The astrometry lieklihood object is passed as a first parameter
# since the obserable-based priors depend on the observation 
# epochs.

@planet b Visual{KepOrbit} begin
    # Instead of a prior on sma
    # a ~ Uniform(0.001, 10000)

    # Put a prior on period:
	P ~ Uniform(0.001, 2000) # yrs
    a = cbrt(system.M * b.P^2)

    # Keep sine prior on inclination
    i ~ Sine()

    # Rest are uniform
    e ~ Uniform(0.0, 1.0)
    ω ~ UniformCircular()
    Ω ~ UniformCircular()
    τ ~ UniformCircular(1.0)
end astrom_like obs_prior

Sine()

A custom univariate distribution. The pdf is a sine function defined between 0 and π. This is a common prior distribution used when fitting orbits to astrometry.

The full Distributions.jl interface is not yet defined for this distribution, but the following methods work: pdf, logpdf, minimum, maximum, insupport, mean, var, cdf, quantile

mjd("2020-01-01")

Get the modfied julian day of a date, or in general a UTC timestamp.

mjd(Date("2020-01-01"))

Get the modfied julian day of a Date or DateTime object.

mjd()

Get the current modified julian day of right now.

gaia_plx(gaia_id=12123)

Get a distribution (truncated Normal) of parallax distance in mas of a source with GAIA catalog id gaia_id.

The method signature of Octofitter.hmc is as follows:

advancedhmc(
    [rng::Random.AbstractRNG],
    model::Octofitter.LogDensityModel
    target_accept::Number=0.8,
    ensemble::AbstractMCMC.AbstractMCMCEnsemble=MCMCSerial();
    adaptation,
    iterations,
    drop_warmup=true,
    max_depth=12,
    initial_samples= pathfinder ? 500 : 250_000,  # deprecated
    initial_parameters=nothing, # deprecated
    step_size=nothing,
    verbosity=2,
)

The only required arguments are system, adaptation, and iterations. The two positional arguments are system, the model you wish to sample; and targetaccept, the acceptance rate that should be targeted during windowed adaptation. During this time, the step size and mass matrix will be adapted (see AdvancedHMC.jl for more information). The number of steps taken during adaptation is controlled by adaptation. You can prevent these samples from being dropped by pasing includeadaptation=false. The total number of posterior samples produced are given by iterations. These include the adaptation steps that may be discarded. treedepth controls the maximum tree depth of the sampler. initialparameters is an optional way to pass a starting point for the chain. If you don't pass a default position, one will be selected by drawing initial_samples from the priors. The sample with the highest posterior value will be used as the starting point.

Visual{OrbitType}(..., plx=...)

This wraps another orbit to add the parallax distance field plx, thus allowing projected quantities to be calculated. It forwards everything else to the parent orbit.

For example, the KepOrbit type supports calculating x and y positions in AU. A Visual{KepOrbit} additionally supports calculating projected right ascension and declination offsets.

sonora_photometry_interpolator(:Keck_L′, [metalicity="+0.0"])

Given a supported photometric band and [M/H] metalicity (default=solar), return a function of temperature (K) and mass (M_jup) that gives the absolute magnitude of the planet in that bandpass.

Supported bands: :MKOY, :MKOZ, :MKOJ, :MKOH, :MKOK, :MKOL′, :MKOM′, :TwoMASSJ, :TwoMASSH, :TwoMASSKs, :KeckKs, :KeckL′, :KeckMs, :SDSSg′, :SDSSr′, :SDSSi′, :SDSSz′, :IRAC36, :IRAC45, :IRAC57, :IRAC79, :WISEW1, :WISEW2, :WISEW3, :WISE_W4

Supported metalicities: "+0.0", "-0.5", "+0.5"

itp = sonora_cooling_interpolator()

Create a function mapping (ageMyr, massMjup) -> temp_K using Sonora Bobcat cooling model grids.

plotchains(
    chain, planet_key;
    N=1500,
    ii = rand(1:size(chain,1)*size(chain,3), N),
    color=length(model.system.planets) == 0 || !haskey(chain, string(planet_key)*"_a") ? nothing : string(planet_key)*"_a",
    colorbartitle=color,
    clims=nothing,
    cmap=:plasma,
    alpha=30/length(ii),
    attime=nothing,
    kwargs...,
)

Draw samples from a posterior chain for a given planet given by name planet_key and visualize them in some way. Use kind to control what plot is made. A few options: :astrometry, :radvel, :trueanom, :meananom, :eccanom, :x, :y, :z, (:x, :y), :raoff, :decoff, :pmra, :pmdec, :accra, :accdec, :radvel, :posangle, :projectedseparation. See PlanetOrbits documentation for more details.

Inputs:

• chain The chain to draw from
• planet_key Planet name in the model (symbol)
• N=1500 Number of samples to draw for the plot
• kind=nothing Specify what kind of plot to make.
• ii=... Specific row numbers to use, if you want to e.g. plot the same 100 samples in a few different plots
• color="planetkeya" Column name to to map colors to. Semi-major axis by default but can be any column or an arbitrary array.
• colorbartitle=color Name for colourbar
• clims=nothing Tuple of colour limits (min and max)
• cmap=:plasma Colormap
• alpha=... Transparency of the lines

projectpositions(chains.planets[1], mjd("2020-02-02"))

Given the posterior for a particular planet in the model and a modified julian date(s), return ra and dec offsets in mas for each sampling in the posterior.