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 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.

When to consider Thiele-Innes

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
Different syntax for `θ_at_epoch_to_tperi`

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)
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.14 seconds
Compute duration  = 5.14 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

Use `describe(chains)` for summary statistics and quantiles.
Note

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)
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.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)
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 ] = 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 
──────────────────────────