# 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_like = PlanetRelAstromLikelihood(
(epoch = 50000, ra = -505.7637580573554, dec = -66.92982418533026, σ_ra = 10, σ_dec = 10, cor=0),
(epoch = 50120, ra = -502.570356287689, dec = -37.47217527025044, σ_ra = 10, σ_dec = 10, cor=0),
(epoch = 50240, ra = -498.2089148883798, dec = -7.927548139010479, σ_ra = 10, σ_dec = 10, cor=0),
(epoch = 50360, ra = -492.67768482682357, dec = 21.63557115669823, σ_ra = 10, σ_dec = 10, cor=0),
(epoch = 50480, ra = -485.9770335870402, dec = 51.147204404903704, σ_ra = 10, σ_dec = 10, cor=0),
(epoch = 50600, ra = -478.1095526888573, dec = 80.53589069730698, σ_ra = 10, σ_dec = 10, cor=0),
(epoch = 50720, ra = -469.0801731788123, dec = 109.72870493064629, σ_ra = 10, σ_dec = 10, cor=0),
(epoch = 50840, ra = -458.89628893460525, dec = 138.65128697876773, σ_ra = 10, σ_dec = 10, cor=0),
)
@planet b ThieleInnesOrbit begin
e ~ Uniform(0.0, 0.5)
A ~ Normal(0, 10000) # milliarcseconds
B ~ Normal(0, 10000) # milliarcseconds
F ~ Normal(0, 10000) # milliarcseconds
G ~ Normal(0, 10000) # milliarcseconds

θ ~ UniformCircular()
tp = θ_at_epoch_to_tperi(system,b,50000.0)  # reference epoch for θ. Choose an MJD date near your data.
end astrom_like

@system TutoriaPrime begin
M ~ truncated(Normal(1.2, 0.1), lower=0)
plx ~ truncated(Normal(50.0, 0.02), lower=0)
end b

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


We now sample from the model as usual:

using Random
Random.seed!(0)

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

Iterations        = 1:1:1000
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 4.61 seconds
Compute duration  = 4.61 seconds
parameters        = M, plx, b_e, b_A, b_B, b_F, b_G, b_θy, b_θx, b_θ, b_tp
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   Flo ⋯

M       1.2183       0.0931      0.0042   484.0294   337.4978    1. ⋯
plx      50.0017       0.0202      0.0007   878.8417   521.0876    0. ⋯
b_e       0.3877       0.1021      0.0052   353.9388   391.8210    1. ⋯
b_A     157.5875     709.6447    125.4662    35.1322   394.1444    1. ⋯
b_B    -671.3603     257.4412     44.8047    43.0665    28.2628    1. ⋯
b_F    1182.5071     583.2842     71.8534    67.8065    35.4212    1. ⋯
b_G     308.3946     343.3486     51.8785    35.3775   209.8055    1. ⋯
b_θy      -0.1275       0.0172      0.0006   747.5913   654.3178    1. ⋯
b_θx      -1.0002       0.0994      0.0039   667.2774   511.5619    1. ⋯
b_θ      -1.6977       0.0127      0.0005   643.6210   425.0584    0. ⋯
b_tp   15892.7442   32994.4514   2512.8620   216.1465   585.2850    1. ⋯
2 columns omitted

Quantiles
parameters          2.5%         25.0%        50.0%        75.0%        97.5 ⋯
Symbol       Float64       Float64      Float64      Float64      Float6 ⋯

M        1.0280        1.1636       1.2127       1.2796       1.398 ⋯
plx       49.9627       49.9883      50.0023      50.0151      50.040 ⋯
b_e        0.1242        0.3337       0.4186       0.4656       0.497 ⋯
b_A    -1032.8015     -457.7115     157.0750     779.9485    1338.311 ⋯
b_B    -1069.0533     -865.5326    -711.2323    -522.2641     -72.616 ⋯
b_F      -23.5728      732.3522    1194.2301    1614.1392    2235.778 ⋯
b_G     -395.7447       56.8582     401.8087     594.4841     798.568 ⋯
b_θy       -0.1651       -0.1382      -0.1259      -0.1152      -0.097 ⋯
b_θx       -1.2021       -1.0628      -0.9945      -0.9338      -0.818 ⋯
b_θ       -1.7242       -1.7058      -1.6971      -1.6893      -1.672 ⋯
b_tp   -48997.4305   -12109.2660   23030.0212   48413.5625   49821.231 ⋯
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)
octocorner(model, results, small=false)

## 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(results, :b, :) # colon means all rows
1000-element Vector{ThieleInnesOrbit{Float64}}:
ThieleInnesOrbit{Float64}(0.36354368256064273, 48195.97013158614, 1.3559047603538406, 49.99227342286215, -746.1625377193155, -832.2954540915464, 1572.4464751657574, -31.52871629775228, -27.660669779608366, -16.592715218472797, 66537.48296668445, 0.03449083630945921, 1.4636937789323032)
ThieleInnesOrbit{Float64}(0.4720111935425532, 49668.25741937575, 1.12493024638176, 50.007837325674636, -228.98326906728744, -984.5315134001428, 2337.777649418175, 467.5216241379114, -44.11001727289982, -9.025520835177861, 116392.88103726132, 0.019717128856984456, 1.6697182673132045)
ThieleInnesOrbit{Float64}(0.4980592121065302, 49380.88534970021, 1.25546835481278, 50.01275144562243, -437.1567841673255, -1046.6264241813763, 1989.8340566137606, 184.45027173659855, -35.05906701247158, -12.12100757678117, 87954.88082789727, 0.026092166936566907, 1.7275803108249494)
ThieleInnesOrbit{Float64}(0.4983152972280249, 48743.6888054888, 1.2718189697025561, 50.006179357063466, -733.3267727500815, -1096.0950590020104, 1724.4650910413409, 106.71469414429252, -29.24374391131641, -18.892580222819557, 80035.37137836032, 0.02867398993625261, 1.7281688681673697)
ThieleInnesOrbit{Float64}(0.482894804646395, 48298.23706322935, 1.4995717311579722, 50.00482211941061, -959.4778053482489, -1049.8172714377051, 1415.241918455757, -67.40332008997015, -22.62102879848148, -22.755531051627095, 65343.6306194151, 0.035120996670874705, 1.6934240418492894)
ThieleInnesOrbit{Float64}(0.43711809435253846, 47790.21475954833, 1.3901593032314394, 50.006980261498754, -1045.6712350463536, -861.0495389867634, 1116.4886842673309, -269.08594298050383, -16.88245576539567, -22.16555603825093, 55861.895976966785, 0.041082268929676154, 1.5978558046783493)
ThieleInnesOrbit{Float64}(0.45792318671258564, 47794.31347337137, 1.280689493353286, 50.00849862381803, -1048.9695518844553, -883.1307634170082, 1062.00772563123, -279.90704311260833, -15.398791067143975, -22.508900806273587, 56919.86176818299, 0.04031867545276021, 1.639973784946812)
ThieleInnesOrbit{Float64}(0.33685295655193714, 48942.832581962364, 1.1889530133569115, 50.005099595336176, -485.1889312294952, -791.2301432491122, 817.6131007887748, -86.3783139233569, -9.967926416160756, -13.173641896589874, 32393.539640449202, 0.0708454049455498, 1.419831808903814)
ThieleInnesOrbit{Float64}(0.32117288523960974, 47330.99450706643, 1.1650499738171631, 50.00406430139065, -886.308624556389, -715.367681488608, 880.1294699983442, -228.71115750012444, -13.03139927109989, -18.919029361425906, 45490.984466731985, 0.05044809340465366, 1.3950836064107217)
ThieleInnesOrbit{Float64}(0.10858833937852247, 27842.026733726365, 1.2485333979083482, 50.02388594657734, -21.581380298094366, -559.5273638010947, 805.4440421179258, 153.01474280580953, -12.421358172907569, -3.3136534730462164, 22350.322963272545, 0.1026801016351541, 1.1151826336162194)
⋮
ThieleInnesOrbit{Float64}(0.49386073026902694, -54142.56263180979, 1.382345999835243, 49.99529766073284, 242.48029887413668, -919.314881579546, 2334.0019423482195, 656.1887716032575, -44.61133322400853, -0.3344587796688264, 104916.12272970915, 0.021873982508481382, 1.7179876534822074)
ThieleInnesOrbit{Float64}(0.4977796479467208, 49426.37971464023, 1.1976492134024725, 49.94748071257903, -447.45383294035895, -1047.7232639040787, 1634.7115979647408, 144.7397196044447, -27.032638279120388, -13.0947377283813, 70207.70822307114, 0.03268777021115183, 1.7269382485352798)
ThieleInnesOrbit{Float64}(0.4982715632568989, 49676.46815721267, 1.2330235248395707, 49.96860883950446, -315.8326804741043, -1007.2201203549718, 1801.0189783154683, 117.91827502113541, -30.64568630387775, -8.985988930347482, 74696.06675959166, 0.030723618163638496, 1.7280683265668149)
ThieleInnesOrbit{Float64}(0.4691908574862223, 46205.11937448415, 1.3190717388824915, 50.00694981024282, -1518.9559256434197, -907.5219726221013, 1260.4976918011903, -310.52542934884815, -20.616314157336664, -31.671577282907958, 83342.52709519042, 0.027536163270237352, 1.6636804428540914)
ThieleInnesOrbit{Float64}(0.3610579147195968, 48509.908478704834, 1.4637296947267235, 49.99390507342605, -575.3711351663536, -877.6243868855549, 2214.7244805529317, 331.83456391559724, -42.25672610483789, -14.82269476292227, 97847.17161557954, 0.023454264395741994, 1.4595116142214921)
ThieleInnesOrbit{Float64}(0.46601658787267486, 49128.109278355296, 1.0586612834073497, 50.003907766293025, -476.6522484845627, -1063.3438647327082, 1459.8094209345604, 308.72361033124406, -24.870405785362472, -16.46842033104126, 70634.27747464809, 0.03249036467133166, 1.6569353331122758)
ThieleInnesOrbit{Float64}(0.48146286286609474, 48757.752181599004, 1.2697013237493116, 49.99879100220449, -697.1049295712749, -988.3559165559998, 1935.1219642639132, 5.325750070633453, -34.12866672129228, -15.873017271883471, 87702.18357269332, 0.02616734658088818, 1.6902675416677857)
ThieleInnesOrbit{Float64}(0.4351070127472574, 48698.54115022885, 1.190332802613484, 50.015970698251415, -607.9393007904579, -943.3782601803766, 1952.1223538158251, 181.2293507580277, -35.57676447037571, -15.255718747689183, 91326.73851889443, 0.025128822847128706, 1.5938925902595544)
ThieleInnesOrbit{Float64}(0.4204790250745677, -11403.47666888393, 1.0535984785854713, 49.96898521878589, 534.9253015879914, -685.9540095423714, 1438.671824549482, 572.6082528511167, -26.271679204799778, 5.74408538766902, 62958.47313397388, 0.036451542091305, 1.565607283281273)

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),
)
pairplot(table)

We can also convert the orbit objects into Campbell parameters:

orbits_campbell = Visual{KepOrbit}.(orbits_ti)
orbits_campbell[1]
KepOrbit{Float64}
─────────────────────────
a   [au ] = 35.6
e         = 0.36354368
i   [°  ] = 65.1
ω   [°  ] = 239.0
Ω   [°  ] = 13.0
tp  [day] = 48200.0
M   [M⊙ ] = 1.36
period      [yrs ] : 182.2
mean motion [°/yr] : 1.98
──────────────────────────
plx [mas] = 50.0
distance    [pc  ] : 20.0
──────────────────────────