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 an alternative parameterization that replaces the angular elements (inclination i, longitude of ascending node Ω, and argument of periastron ω) with the Thiele-Innes constants (A, B, F, G). This avoids the coordinate singularities where ω, Ω, and tp become poorly defined as eccentricity and/or inclination approach zero.
The Thiele-Innes parameterization may be useful when:
- The orbit is nearly circular (e < 0.1) and/or nearly face-on (i near 0° or 180°)
- You want to avoid angular coordinate singularities
Both parameterizations should give consistent results for the physical orbital parameters. Choose based on your preference or specific analysis needs.
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_obs = PlanetRelAstromObs(
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,
observations=[astrom_obs],
variables=@variables begin
e ~ Uniform(0.0, 0.5)
# Thiele-Innes constants A, B, F, G are in milliarcseconds (not AU like semi-major axis).
# Set the prior width to encompass the expected angular separation of your target.
# A rough guide: if your astrometry spans ~500 mas, use Normal(0, 1000) or similar.
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],
observations=[],
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
When using ThieleInnesOrbit, the θ_at_epoch_to_tperi function requires the Thiele-Innes constants A, B, F, G as keyword arguments instead of the Campbell angular elements i, Ω, ω:
# Campbell (Visual{KepOrbit}):
tp = θ_at_epoch_to_tperi(θ, epoch; M, e, a, i, Ω, ω)
# Thiele-Innes:
tp = θ_at_epoch_to_tperi(θ, epoch; system.plx, M, e, A, B, F, G)Note that system.plx (parallax) is also required for Thiele-Innes since the constants are in angular units.
Initialize the starting points, and confirm the data are entered correcly:
init_chain = initialize!(model)
octoplot(model, init_chain)
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.25 seconds
Compute duration = 5.25 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.2182 0.0941 0.0047 392.6249 283.8885 ⋯
plx 49.9993 0.0183 0.0006 919.0183 593.7641 ⋯
b_e 0.2833 0.1500 0.0133 122.8939 201.3530 ⋯
b_A 96.9893 635.8191 64.0078 109.7296 356.4671 ⋯
b_B -440.8570 393.2580 57.0427 60.9548 47.6143 ⋯
b_F 629.5963 622.5680 79.3636 65.6041 43.2140 ⋯
b_G 178.2888 371.8124 40.2630 98.5434 240.9926 ⋯
b_θx -0.1306 0.0188 0.0008 557.7402 522.3495 ⋯
b_θy -0.9962 0.0969 0.0037 684.4661 568.7358 ⋯
b_θ -1.7011 0.0127 0.0006 469.1090 594.0827 ⋯
b_M 1.2182 0.0941 0.0047 392.6249 283.8885 ⋯
b_tp 31302.8220 20500.2219 1403.5477 282.3258 418.5125 ⋯
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.0340 1.1549 1.2156 1.2811 ⋯
plx 49.9635 49.9876 49.9992 50.0115 5 ⋯
b_e 0.0155 0.1600 0.3005 0.4157 ⋯
b_A -990.5591 -461.3443 117.3063 674.9287 115 ⋯
b_B -1015.0970 -728.9848 -515.3929 -217.4440 47 ⋯
b_F -619.1225 215.2350 632.3320 1060.4443 175 ⋯
b_G -495.4642 -137.1875 225.8036 522.1080 69 ⋯
b_θx -0.1690 -0.1431 -0.1302 -0.1175 - ⋯
b_θy -1.1895 -1.0634 -0.9961 -0.9254 - ⋯
b_θ -1.7238 -1.7098 -1.7010 -1.6918 - ⋯
b_M 1.0340 1.1549 1.2156 1.2811 ⋯
b_tp -18711.7850 21405.5934 38441.9607 48273.2104 4982 ⋯
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
The Thiele-Innes parameterization may reveal more complex posterior structure (e.g., multimodality) that Campbell masks through its angular parameterization. If your corner plot shows unexpected bimodality in the A, B, F, G parameters, this may reflect genuine orbital ambiguities rather than sampling issues.
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(model, results, :b, :) # colon means all rows1000-element Vector{ThieleInnesOrbit{Float64}}:
ThieleInnesOrbit{Float64}(0.41977044100683847, 49802.024200024796, 1.362857750709097, 50.01370995913912, -218.47939825247408, -901.9706471000799, 1561.7033605941408, 291.182829463059, -27.258627662805626, -8.856019829366323, 59245.104877234124, 0.038736254045002264, 1.5642607213030104)
ThieleInnesOrbit{Float64}(0.3719332054212324, 48976.16421513368, 1.3393441658597443, 50.01488652279577, -483.6944916541775, -889.4676453903868, 1477.8223061473193, 211.09946675106684, -25.965042625834425, -13.896275175646368, 59623.802128452335, 0.03849022288956321, 1.4779630182030248)
ThieleInnesOrbit{Float64}(0.2484980266254969, 9622.837831002107, 1.4272282144615207, 49.99767962524593, 278.4500272374752, -571.565029794722, 1235.089118265181, 462.1694738395639, -23.149056173019606, 1.3781657708196244, 41500.01590568711, 0.05529957961131405, 1.288928621090679)
ThieleInnesOrbit{Float64}(0.24205741465734046, 6489.472087829119, 1.3826978779184815, 50.01398940417366, 131.072146086494, -618.2645081143867, 1307.1421678839572, 382.069107488912, -24.1432298895384, -1.0744799986569975, 44186.70089472981, 0.05193719800251172, 1.2801258911832554)
ThieleInnesOrbit{Float64}(0.3960818748266758, 9437.724879785237, 1.0742380314349271, 49.99438770338508, 631.9661384371996, -605.3384818467999, 1026.6248192067014, 525.5326635440143, -16.935363686339784, 7.811848120407516, 42359.41374425026, 0.054177648616745816, 1.5204299373722747)
ThieleInnesOrbit{Float64}(0.21252514110468576, 21301.643684268347, 1.136160317505899, 49.996969897513395, 717.4551560327398, -326.9278151487115, 710.8451229150809, 500.19986632988923, -12.959424328263626, 10.695297500698969, 31600.112466928542, 0.07262421726660412, 1.2408721216045462)
ThieleInnesOrbit{Float64}(0.06245695241696649, 32558.05976513974, 1.2393300207605433, 50.01770360459174, 150.20999856216403, -479.0321935591564, 678.4080356279829, 263.9189684855658, -10.579323212493891, -0.9265113068225784, 18271.595365507354, 0.12560115236459646, 1.064535288595142)
ThieleInnesOrbit{Float64}(0.2641464181852836, 49288.90138238088, 1.1983302352779335, 49.99046085436754, -309.1517754227913, -686.8835727165989, 1248.0482985065316, 32.01727435362288, -21.33467272740086, -7.649222470172196, 44539.77087813213, 0.051525488079555815, 1.310699051510308)
ThieleInnesOrbit{Float64}(0.247153415789775, 49705.78763156389, 1.096959007873852, 50.008198083933195, -177.30087115288197, -684.0795516962748, 1144.7425631777103, 125.27071951022968, -19.154523526936824, -6.026030481716655, 40500.01347154029, 0.05666500419956683, 1.2870834267967388)
ThieleInnesOrbit{Float64}(0.38813439541018485, 23354.61089513812, 1.3635755190838434, 49.99877529076708, 1052.0088628521985, -135.68356485267597, 85.08399209204835, 587.8563478662782, 0.22180130909668797, -17.578158395122458, 30567.377766219222, 0.07507786408762654, 1.5062176675436811)
⋮
ThieleInnesOrbit{Float64}(0.3561527685422276, -2387.976914749517, 1.3796859872403215, 49.97105316259486, 1050.5131424083636, -294.5658222270212, 1130.9666858780981, 733.8040816235641, -23.120839931472936, 16.834477138763788, 55764.34944526468, 0.041154132636298905, 1.4513190343616242)
ThieleInnesOrbit{Float64}(0.1835652418746899, 19195.2277115923, 1.2966956829872762, 49.99384572630439, 599.9286504861303, -467.76837402211726, 922.4650064341998, 357.88761909846556, -15.938868807432753, 9.689513656485792, 33180.56945012991, 0.06916498033273985, 1.2040245949458965)
ThieleInnesOrbit{Float64}(0.13517044966625444, 24236.881094767457, 1.2384546161941417, 49.95495634455917, 94.80077020817187, -552.4090480895537, 919.2163915642483, 139.381923284966, -14.851575960193799, 0.2737727680918617, 26356.746895712895, 0.08707195324703104, 1.1456851351232096)
ThieleInnesOrbit{Float64}(0.4405080095646909, -12741.808718777269, 1.1732305475460045, 50.01126507316779, 874.4722916780488, -559.2219159827689, 1378.7827598583328, 715.6821630636338, -26.180209343721966, 12.301190241938832, 65119.02534069148, 0.03524213425248073, 1.6045782530257806)
ThieleInnesOrbit{Float64}(0.4787085515468824, -23583.323834162293, 1.3373823622279948, 50.02753751165003, 827.9496526591284, -640.1665729904777, 1679.4009137046987, 769.8026019300451, -32.383861273652016, 11.075511324939834, 75607.20103413027, 0.030353371134733253, 1.6842284218891277)
ThieleInnesOrbit{Float64}(0.3749702697759691, 48274.40984743167, 1.257236754965335, 50.02366354474863, -740.8737631118752, -843.2662383487615, 1188.1023321358205, -89.49583351392641, -18.84329074601381, -17.067163871240638, 51668.88938483146, 0.04441615565518759, 1.4831883860431936)
ThieleInnesOrbit{Float64}(0.2522447270157875, 49232.95275739171, 1.2410962856937933, 49.997133401036265, -333.8231112573156, -722.5919646812857, 961.9307959470615, 86.69662883672106, -14.987831469194278, -10.243110262095824, 33522.10743880381, 0.06846029706326946, 1.2940911396147654)
ThieleInnesOrbit{Float64}(0.08233825077710173, 49424.35908263738, 1.399820040394023, 49.99446365314282, -223.74786084303005, -567.8071228974902, 937.1768831190909, 51.38919853580408, -15.533433580973869, -6.152493771120915, 27109.263740675, 0.08465495246940269, 1.0860259155010676)
ThieleInnesOrbit{Float64}(0.22044839071113986, 24989.72850519949, 1.1660102384028084, 50.003384912933036, 589.1617070049595, -411.7275513865329, 712.064224542773, 461.36222801968, -11.876658502969825, 7.730603833475597, 27235.260727053465, 0.0842633179262253, 1.251230359143008)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)
We can also convert the orbit objects into Campbell parameters:
orbits_campbell = Visual{KepOrbit}.(orbits_ti)
orbits_campbell[1]KepOrbit{Float64}
─────────────────────────
a [au ] = 33.0
e = 0.41977044
i [° ] = 60.4
ω [° ] = 252.0
Ω [° ] = 19.7
tp [day] = 49800.0
M [M⊙ ] = 1.36
period [yrs ] : 162.2
mean motion [°/yr] : 2.22
──────────────────────────
plx [mas] = 50.0
distance [pc ] : 20.0
──────────────────────────