Fitting in different orbital bases

In this tutorial, we show how one can perform orbit-fits in different coordinate bases amongst the ones supported by orbitize. Currently fitting in different bases is only supported in MCMC, so we will use MCMC to perform an orbit-fit in an orbital basis distinct from the default one. For a general introduction to MCMC, be sure to check out the MCMC Introduction tutorial first!

The “standard” and “XYZ” bases

The default way to define an orbit in orbitize is through what we call the ‘standard basis’, which consists of eight parameters: semi-major axis (sma), eccentricity (ecc), inclination (inc), argument of periastron (aop), position angle of the nodes (pan), epoch of periastron expressed as a fraction of the period past a reference epoch (tau), parallax (plx) and total system mass (mtot). Each orbital element has an associated default prior; to see how to explore and modify these priors check out the Modifying priors tutorial.

An alternative way to define an orbit is through its position and velocity components in XYZ space for a given epoch; we will call this the ‘XYZ basis’. The orbit is thus defined with the array (\(x\), \(y\), \(z\), \(\dot{x}\), \(\dot{y}\),\(\dot{z}\), plx, mtot), with position coordinates measured in AU and velocity components in \(\text{km s}^{-1}\). In this basis, the sky-plane coordinates (\(x,y\)) are the separations of the planet relative to the primary, with the positive \(x\) and \(y\) directions coinciding with the positive RA and Dec directions, respectively. The \(z\) direction is the line-of-sight coordinate, such that movement in the positive \(z\) direction causes a redshift. The default priors are uniform all uniform.

Setting up Sampler in the XYZ basis

The easiest way to run an orbit-fit in an alternative orbital basis in orbitize is through the orbitize.driver.Driver interface. The process is exactly like initializing a regular Driver object, but setting the fitting_basis keyword to ‘XYZ’:

import numpy as np

import orbitize
from orbitize import driver
import multiprocessing as mp

filename = "{}xyz_test_data.csv".format(orbitize.DATADIR) # a file with input in radec since rn it only works for that

# system parameters
num_secondary_bodies = 1
system_mass = 1.75 # [Msol]
plx = 51.44 # [mas]
mass_err = 0.05 # [Msol]
plx_err = 0.12 # [mas]

# MCMC parameters
num_temps = 5
num_walkers = 20
num_threads = mp.cpu_count() # or a different number if you prefer


my_driver = driver.Driver(
    filename, 'MCMC', num_secondary_bodies, system_mass, plx, mass_err=mass_err, plx_err=plx_err,
    mcmc_kwargs={'num_temps': num_temps, 'num_walkers': num_walkers, 'num_threads': num_threads},
    system_kwargs={'fitting_basis': 'XYZ'}
)

s = my_driver.sampler

Converting ra/dec data points in data_table to sep/pa. Original data are stored in input_table.

(Properly) initializing walkers in the XYZ basis

In the standard basis at this point we would be ready to use the s.run_sampler method to start the sampling, but with the XYZ basis we have to make sure that all our walkers are initialized in a valid region of parameter space. This is because randomly generated values of (\(x\), \(y\), \(z\), \(\dot{x}\), \(\dot{y}\), \(\dot{z}\)) can result in unbound, invalid orbits with, for example, negative eccentricities (which is not cool). This can be easily corrected with the s.validate_xyz_positions method:

s.validate_xyz_positions()

All walker positions validated.

/home/sblunt/Projects/orbitize/orbitize/basis.py:944: RuntimeWarning: invalid value encountered in arccos
  eanom = np.arccos(cos_eanom)
/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/units/quantity.py:481: RuntimeWarning: invalid value encountered in sqrt
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)

---

After this is done, the sampler can be run and the results saved normally:

total_orbits = 600 # number of steps x number of walkers (at lowest temperature)
burn_steps = 10 # steps to burn in per walker
thin = 2 # only save every 2nd step

s.run_sampler(total_orbits, burn_steps=burn_steps, thin=thin)
s.results.save_results('my_posterior.hdf5')

/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))

Starting Burn in

/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))
/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))
/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))
/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))
/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))
/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))
/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))

/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))
/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))

10/10 steps of burn-in complete
Burn in complete. Sampling posterior now.

/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))
/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))

/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))

/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))
/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))

/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))
/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))
/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))

/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))
/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))
/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))
/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))
/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))
/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))
/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))
/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))
/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))
/data/user/sblunt/miniconda3/envs/python3.7/lib/python3.7/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)
/home/sblunt/Projects/orbitize/orbitize/kepler.py:112: RuntimeWarning: invalid value encountered in sqrt
  tanom = 2.*np.arctan(np.sqrt((1.0 + ecc)/(1.0 - ecc))*np.tan(0.5*eanom))

30/30 steps completed
Run complete

Loading and converting results

You can load the results as you normally would. The orbit posteriors are saved in the results.post attribute, and the basis you used for the fit in the results.fitting_basis attribute:

myResults = orbitize.results.Results() # create empty Results object
myResults.load_results('my_posterior.hdf5')
print('The used basis for the fit was ', myResults.fitting_basis)
print('The posteriors are ', myResults.post)

Converting ra/dec data points in data_table to sep/pa. Original data are stored in input_table.
The used basis for the fit was  XYZ
The posteriors are  [[-1.55895340e+01 -3.20352269e+01  9.38119986e+00 ... -7.46195838e-03
   5.15299823e+01  1.72400897e+00]
 [-1.56081092e+01 -3.20562773e+01  1.49810951e+00 ...  8.53387015e-02
   5.13883909e+01  1.73579787e+00]
 [-1.55612462e+01 -3.20801251e+01 -2.96308303e-01 ...  7.18261775e-01
   5.14206555e+01  1.76608956e+00]
 ...
 [-1.55671475e+01 -3.20089823e+01  3.62094122e+01 ...  3.72598271e-01
   5.15487665e+01  1.74432532e+00]
 [-1.55794884e+01 -3.20303979e+01  2.75912018e+01 ... -2.40407388e-01
   5.14895599e+01  1.78876637e+00]
 [-1.55837449e+01 -3.20532712e+01  1.93382300e+01 ... -1.24185021e-01
   5.14598433e+01  1.80483891e+00]]

Let’s convert back to the good old standard basis:

xyz_posterior = myResults.post

standard_posterior = myResults.system.basis.to_standard_basis(xyz_posterior)

print('My posterior in standard basis is ', standard_posterior)

My posterior in standard basis is  [[  7.89543529   6.92670011   4.84646336 ...   0.36153132  -6.60018059
    1.52254522]
 [  1.00000429   1.00000659   0.96238341 ...   1.00053929   1.00014326
    0.9999977 ]
 [  1.43499876   1.22274162   0.79516331 ...   1.46937244   2.46248534
    1.23409917]
 ...
 [-15.56714749 -32.00898233  36.20941224 ...   0.37259827  51.54876654
    1.74432532]
 [-15.57948843 -32.03039793  27.59120179 ...  -0.24040739  51.48955988
    1.78876637]
 [-15.58374489 -32.0532712   19.33822997 ...  -0.12418502  51.45984325
    1.80483891]]

And we’re done! Enjoy the XYZ basis.