GitHub

Fit with a Thiele-Innes Basis

This example shows how to fit relative astrometry using a Thiele-Innes orbital basis instead of the traditional Campbell basis used in other tutorials. The Thiele-Innes basis is more suitable then Campbell for fitting low-eccentricity orbits, because it does not have the issues where ω, Ω, and tp become degenerate as eccentricity and/or inclination fall to zero.

At the end, we will convert our results back into the Campbell basis to compare.

using Octofitter
using CairoMakie
using PairPlots
using Distributions

astrom_dat = Table(;
    epoch = [50000, 50120, 50240, 50360, 50480, 50600, 50720, 50840],
    ra    = [-505.7637580573554, -502.570356287689, -498.2089148883798, -492.67768482682357, -485.9770335870402, -478.1095526888573, -469.0801731788123, -458.89628893460525],
    dec   = [-66.92982418533026, -37.47217527025044, -7.927548139010479, 21.63557115669823, 51.147204404903704, 80.53589069730698, 109.72870493064629, 138.65128697876773],
    σ_ra  = [10, 10, 10, 10, 10, 10, 10, 10],
    σ_dec = [10, 10, 10, 10, 10, 10, 10, 10],
    cor   = [0, 0, 0, 0, 0, 0, 0, 0]
)

astrom_like = PlanetRelAstromLikelihood(
    astrom_dat,
    name = "GPI",
    variables = @variables begin
        # Fixed values for this example - could be free variables:
        jitter = 0        # mas [could use: jitter ~ Uniform(0, 10)]
        northangle = 0    # radians [could use: northangle ~ Normal(0, deg2rad(1))]
        platescale = 1    # relative [could use: platescale ~ truncated(Normal(1, 0.01), lower=0)]
    end
)

planet_b = Planet(
    name="b",
    basis=ThieleInnesOrbit,
    likelihoods=[astrom_like],
    variables=@variables begin
        e ~ Uniform(0.0, 0.5)
        A ~ Normal(0, 1000) # milliarcseconds
        B ~ Normal(0, 1000) # milliarcseconds
        F ~ Normal(0, 1000) # milliarcseconds
        G ~ Normal(0, 1000) # milliarcseconds

        M = system.M
        θ ~ UniformCircular()
        tp = θ_at_epoch_to_tperi(θ, 50000.0; system.plx, M, e, A, B, F, G)  # reference epoch for θ. Choose an MJD date near your data.
    end
)

sys = System(
    name="TutoriaPrime",
    companions=[planet_b],
    likelihoods=[],
    variables=@variables begin
        M ~ truncated(Normal(1.2, 0.1), lower=0.1)
        plx ~ truncated(Normal(50.0, 0.02), lower=0.1)
    end
)

model = Octofitter.LogDensityModel(sys)
LogDensityModel for System TutoriaPrime of dimension 9 and 9 epochs with fields .ℓπcallback and .∇ℓπcallback

Initialize the starting points, and confirm the data are entered correcly:

init_chain = initialize!(model)
octoplot(model, init_chain)
Example block output

We now sample from the model as usual:

results = octofit(model)
Chains MCMC chain (1000×28×1 Array{Float64, 3}):

Iterations        = 1:1:1000
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 5.93 seconds
Compute duration  = 5.93 seconds
parameters        = M, plx, b_e, b_A, b_B, b_F, b_G, b_θx, b_θy, b_θ, b_M, b_tp, b_GPI_jitter, b_GPI_northangle, b_GPI_platescale
internals         = n_steps, is_accept, acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size, is_adapt, loglike, logpost, tree_depth, numerical_error

Summary Statistics
        parameters         mean          std        mcse   ess_bulk   ess_tail ⋯
            Symbol      Float64      Float64     Float64    Float64    Float64 ⋯

                 M       1.2093       0.1006      0.0041   605.4687   503.4560 ⋯
               plx      49.9992       0.0196      0.0007   812.7774   490.3898 ⋯
               b_e       0.2910       0.1404      0.0062   519.6361   605.4305 ⋯
               b_A      44.0218     632.5368     47.8796   184.0012   308.8855 ⋯
               b_B    -516.5600     315.2948     24.6555   196.5038   232.9641 ⋯
               b_F     742.7807     561.4974     45.1048   151.1392   308.1562 ⋯
               b_G     163.6639     380.8540     28.4828   196.0964   420.1776 ⋯
              b_θx      -0.1297       0.0192      0.0008   620.9341   558.5930 ⋯
              b_θy      -0.9999       0.1000      0.0037   738.8254   546.0207 ⋯
               b_θ      -1.6996       0.0130      0.0005   606.8835   445.5038 ⋯
               b_M       1.2093       0.1006      0.0041   605.4687   503.4560 ⋯
              b_tp   30339.9234   23674.9524   1410.4102   341.4213   380.1178 ⋯
      b_GPI_jitter       0.0000       0.0000         NaN        NaN        NaN ⋯
  b_GPI_northangle       0.0000       0.0000         NaN        NaN        NaN ⋯
  b_GPI_platescale       1.0000       0.0000         NaN        NaN        NaN ⋯
                                                               2 columns omitted

Quantiles
        parameters          2.5%        25.0%        50.0%        75.0%        ⋯
            Symbol       Float64      Float64      Float64      Float64      F ⋯

                 M        1.0169       1.1383       1.2079       1.2822        ⋯
               plx       49.9602      49.9856      49.9993      50.0125      5 ⋯
               b_e        0.0264       0.1855       0.2985       0.4183        ⋯
               b_A     -923.5428    -509.1384     -40.9763     581.1510    116 ⋯
               b_B     -973.1430    -748.6135    -569.6612    -346.1234     24 ⋯
               b_F     -318.6783     384.9358     719.9969    1099.5475    190 ⋯
               b_G     -514.2538    -178.3602     211.0709     505.1457     70 ⋯
              b_θx       -0.1706      -0.1415      -0.1292      -0.1158      - ⋯
              b_θy       -1.2064      -1.0677      -0.9952      -0.9287      - ⋯
               b_θ       -1.7226      -1.7087      -1.6999      -1.6914      - ⋯
               b_M        1.0169       1.1383       1.2079       1.2822        ⋯
              b_tp   -31376.4192   18016.5528   40063.4993   48540.2505   4976 ⋯
      b_GPI_jitter        0.0000       0.0000       0.0000       0.0000        ⋯
  b_GPI_northangle        0.0000       0.0000       0.0000       0.0000        ⋯
  b_GPI_platescale        1.0000       1.0000       1.0000       1.0000        ⋯
                                                                1 column omitted

Notice that the fit was very very fast! The Thiele-Innes orbital paramterization is easier to explore than the default Campbell in many cases.

We now display the results:

octoplot(model,results)
Example block output
octocorner(model, results, small=false)
Example block output

Conversion back to Campbell Elements

To convert our chain into the more familiar Campbell parameterization, we have to do a few steps. We start by turning the chain table into a an array of orbit objects, and then convert their type:

orbits_ti = Octofitter.construct_elements(model, results, :b, :) # colon means all rows
1000-element Vector{ThieleInnesOrbit{Float64}}:
 ThieleInnesOrbit{Float64}(0.24114866491057388, 29403.72472838958, 1.2455147077757607, 49.95128688435789, 788.8336039429025, -262.0187643484763, 356.9593404230577, 508.0191994505008, 4.915176674820727, -12.10623280140728, 23654.68765804608, 0.09701812455201572, 1.278891058474485)
 ThieleInnesOrbit{Float64}(0.4568592605637595, 4727.849900030982, 1.332491484231123, 49.95236364447514, 1211.6556760610572, -341.90647545811873, 752.6019001539516, 714.5589880357832, -13.581555664296436, 19.69894768230192, 48472.60741956658, 0.047344955339063625, 1.6377688620041835)
 ThieleInnesOrbit{Float64}(0.4537402499069735, 49120.07855616792, 1.065095876874362, 50.029779264743766, -509.6395612365742, -933.8418270810007, 814.9910685331353, -165.55392515551566, -6.956604472380784, -14.975132801486142, 37454.54258567388, 0.0612724992755656, 1.6313374515235572)
 ThieleInnesOrbit{Float64}(0.4129988480614527, 49896.226642049805, 1.0627277607907768, 50.01120669701405, -154.55497925110402, -839.8031314174071, 592.4332133200999, -102.76833458665106, -0.1733955555075377, -12.124826223227267, 24999.745423348362, 0.09179827212576351, 1.551498739705961)
 ThieleInnesOrbit{Float64}(0.46284936808342236, 48815.161549291, 1.0799119326901607, 50.01114658443287, -630.9182324885733, -875.9775609595908, 703.3441669691937, -313.8256715382573, -4.29360819951474, -15.72320028618697, 36290.866538188595, 0.06323721785569389, 1.6502576654514964)
 ThieleInnesOrbit{Float64}(0.4137839479039797, 48682.33101311054, 1.087458578164918, 49.98895665429712, -619.7957321330287, -895.0968311173378, 1092.3473418883962, -76.08113676456769, -15.698097561715494, -15.522933496816743, 48724.021136599215, 0.0471006575383717, 1.5529685833310245)
 ThieleInnesOrbit{Float64}(0.22306594827936096, 48913.35772004048, 1.4304537258077072, 49.99404933504282, -468.8648355258894, -621.0228747751554, 897.2692500491453, -155.16476051657236, -13.490981098337008, -9.618702924300633, 28548.937786524486, 0.0803859481780995, 1.254679665235274)
 ThieleInnesOrbit{Float64}(0.38504255853303443, 48583.97441028858, 1.0823557927289056, 50.00318620201535, -565.1107167855997, -938.1474580600539, 1252.7344883056562, 166.9362828879431, -20.840869582021988, -16.59113803606729, 58364.34875237033, 0.039320809406858, 1.5007523377854524)
 ThieleInnesOrbit{Float64}(0.45217183358312557, 46943.802738684935, 1.3571934862096982, 49.999118825373586, -1242.9787013031612, -880.342740298963, 871.4029411984403, -325.58626455596055, -11.854278841282369, -26.87758159648409, 58597.03334973959, 0.03916466930586585, 1.628121553513892)
 ThieleInnesOrbit{Float64}(0.3147394784246087, 3931.5235071328643, 1.1726286712175615, 49.985050365734885, 65.85092738457612, -697.4369780909825, 1269.2987565820192, 379.8972217180148, -22.72403848939326, -3.1944739985506714, 46531.227138411516, 0.049320286065545796, 1.385134656947846)
 ⋮
 ThieleInnesOrbit{Float64}(0.4106450137151641, 49100.93751831243, 1.0911072171992084, 49.99409673480623, -445.4861128972526, -792.8406052792428, 569.8977621441844, -269.3327462314166, -1.2256565569569118, -13.169473644218101, 27221.524249270096, 0.08430583873380562, 1.5471071310123297)
 ThieleInnesOrbit{Float64}(0.41248885696109866, 48806.89197927719, 1.079166893557902, 50.03412947324056, -517.3987985323494, -786.2858251558375, 510.8281668015004, -318.4580266430267, -0.3838755595600431, -14.467070589258697, 28697.62616196478, 0.0799694518457767, 1.550545308041463)
 ThieleInnesOrbit{Float64}(0.04468143413825885, 44486.49030418944, 1.3433024196059613, 49.96902265338312, -902.3755160980788, -173.33939965422925, 24.010321794641435, -490.88559801551895, 1.625937192012877, -15.622288338935418, 24996.5621067155, 0.09180996265205582, 1.045725815142836)
 ThieleInnesOrbit{Float64}(0.48549601823886945, 49066.09096070463, 1.114772709467992, 50.02105779100882, -538.4937534835408, -901.5171990345531, 649.1408920682593, -317.2908458749223, -1.6565860041906457, -15.32348937288309, 33430.550480081634, 0.06864779073305106, 1.699187779412098)
 ThieleInnesOrbit{Float64}(0.19790804374152907, 45933.28111543249, 1.3714353694557306, 49.9937764766383, -837.5440069488465, -225.45600883928802, 0.025661078375257773, -518.82087137388, 3.2757990933465986, -14.28401515152966, 23139.002716806466, 0.09918030874254029, 1.22208006568898)
 ThieleInnesOrbit{Float64}(0.46482712022815714, 48821.562177002364, 1.0979255913866341, 50.00338172703536, -644.4593925100303, -926.7720422660256, 813.6897030767583, -218.56502991830902, -7.637196514060883, -16.85364427979508, 40555.097227505765, 0.056588039243827694, 1.6544213937601533)
 ThieleInnesOrbit{Float64}(0.41843517591562635, 48168.22820522976, 1.2153336981206821, 50.000034844054795, -832.7037485550725, -889.0854155117371, 1400.5594870632628, -98.11317174554465, -23.236610928816273, -18.57445811341379, 64723.64341052301, 0.035457420387960195, 1.5617290229331637)
 ThieleInnesOrbit{Float64}(0.4798516583800818, -5621.163928468872, 1.1365929424135295, 49.97030374213769, 499.6386269267906, -793.2917887390207, 1368.8962578904404, 613.897653915194, -23.6742714085268, 3.3316782621765686, 56880.0109291954, 0.0403469232153294, 1.686729663501245)
 ThieleInnesOrbit{Float64}(0.45134814574935855, -7051.51621200089, 1.1583096504941743, 50.0169278928235, 59.08681005033855, -902.4947304233606, 1510.8188832545864, 226.4346260344827, -24.686144770871408, -1.863534417262788, 57448.01319912853, 0.039948003519156644, 1.6264374868322444)

Here is one of those entries:

display(orbits_ti[1])

We can now make a table of results (and visualize them in a corner plot) by querying properties of these objects:

table = (;
    B_a = semimajoraxis.(orbits_ti),
    B_e = eccentricity.(orbits_ti),
    B_i = rad2deg.(inclination.(orbits_ti)),
)
pairplot(table)
Example block output

We can also convert the orbit objects into Campbell parameters:

orbits_campbell = Visual{KepOrbit}.(orbits_ti)
orbits_campbell[1]
KepOrbit{Float64}
─────────────────────────
a   [au ] = 17.4
e         = 0.24114866
i   [°  ] = 48.9
ω   [°  ] = 518.0
Ω   [°  ] = 177.0
tp  [day] = 29400.0
M   [M⊙ ] = 1.25 
period      [yrs ] : 64.8 
mean motion [°/yr] : 5.56 
──────────────────────────
plx [mas] = 50.0 
distance    [pc  ] : 20.0 
──────────────────────────