Repository: sblunt/orbitize (language: python)

Issue:
Changebasis
(Hopefully this is right this time around) Added the conversions.py module for converting between bases and a Jupyter notebook tutorial draft. Also edited sampler.py, results.py and system.py to support the XYZ basis. (Still some corrections need to be made e.g. currently only works for data provided in RA/Dec).


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index 1248101e..c7a7b695 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,8 @@
+# Added by orbitize! Developers
+orbitize/example_data/rebound*.csv
+orbitize/example_data/*test.hdf5
+*.vscode
+
 # Byte-compiled / optimized / DLL files
 __pycache__/
 *.py[cod]
diff --git a/.vscode/launch.json b/.vscode/launch.json
new file mode 100644
index 00000000..17e15f27
--- /dev/null
+++ b/.vscode/launch.json
@@ -0,0 +1,15 @@
+{
+    // Use IntelliSense to learn about possible attributes.
+    // Hover to view descriptions of existing attributes.
+    // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387
+    "version": "0.2.0",
+    "configurations": [
+        {
+            "name": "Python: Current File",
+            "type": "python",
+            "request": "launch",
+            "program": "${file}",
+            "console": "integratedTerminal"
+        }
+    ]
+}
\ No newline at end of file
diff --git a/docs/tutorials/Changing_bases_tutorial.ipynb b/docs/tutorials/Changing_bases_tutorial.ipynb
new file mode 100644
index 00000000..466ee4d4
--- /dev/null
+++ b/docs/tutorials/Changing_bases_tutorial.ipynb
@@ -0,0 +1,228 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Fitting in different orbital bases\n",
+    "\n",
+    "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](https://orbitize.readthedocs.io/en/latest/tutorials/MCMC_tutorial.html) first!"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The \"standard\" and \"XYZ\" bases"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "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](http://orbitize.info/en/latest/tutorials/Modifying_Priors.html).\n",
+    "\n",
+    "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."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Setting up Sampler in the XYZ basis"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "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':"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "import orbitize\n",
+    "from orbitize import driver\n",
+    "import multiprocessing as mp\n",
+    "\n",
+    "filename = \"{}/ \".format(orbitize.DATADIR) # a file with input in radec since rn it only works for that\n",
+    "\n",
+    "# system parameters\n",
+    "num_secondary_bodies = 1\n",
+    "system_mass = 1.75 # [Msol]\n",
+    "plx = 51.44 # [mas]\n",
+    "mass_err = 0.05 # [Msol]\n",
+    "plx_err = 0.12 # [mas]\n",
+    "\n",
+    "# MCMC parameters\n",
+    "num_temps = 5\n",
+    "num_walkers = 20\n",
+    "num_threads = mp.cpu_count() # or a different number if you prefer\n",
+    "\n",
+    "\n",
+    "my_driver = driver.Driver(\n",
+    "    filename, 'MCMC', num_secondary_bodies, system_mass, plx, mass_err=mass_err, plx_err=plx_err, fitting_basis='XYZ',\n",
+    "    mcmc_kwargs={'num_temps': num_temps, 'num_walkers': num_walkers, 'num_threads': num_threads}\n",
+    ")\n",
+    "\n",
+    "s = my_driver.sampler"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### (Properly) initializing walkers in the XYZ basis"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "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:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "s.validate_xyz_positions()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "___________________"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "After this is done, the sampler can be run and the results saved normally:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "total_orbits = 6000 # number of steps x number of walkers (at lowest temperature)\n",
+    "burn_steps = 10 # steps to burn in per walker\n",
+    "thin = 2 # only save every 2nd step\n",
+    "\n",
+    "s.run_sampler(total_orbits, burn_steps=burn_steps, thin=thin)\n",
+    "s.results.save_results('my_posterior.hdf5')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Loading and converting results"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "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:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from orbitize import results \n",
+    "\n",
+    "myResults = orbitize.results.Results() # create empty Results object\n",
+    "myResults.load_results('my_posterior.hdf5') \n",
+    "print('The used basis for the fit was ', myResults.fitting_basis)\n",
+    "print('The posteriors are ', myResults.post)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "It is important to note that, as you can probably see, the results will be saved in the XYZ basis, so we might want to convert those back to the standard basis to work with them. For this we can use the `orbitize.conversions` module, which can convert from the standard basis to the XYZ basis and back. To do this, remember that the XYZ basis gives the position and velocity of the planet at a given epoch; so before converting back to standard we need to check the epoch that was used for the fit through the `results.xyz_epochs` attribute:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(myResults.xyz_epochs)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Since we only fit for one planet, we only have one epoch. Now, we can use that epoch to convert back to the good old standard basis:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from orbitize import conversions\n",
+    "xyz_posterior = myResults.post\n",
+    "xyz_epoch = myResults.xyz_epochs\n",
+    "\n",
+    "standard_posterior = conversions.xyz_to_standard(xyz_epoch, xyz_posterior)\n",
+    "\n",
+    "print('My posterior in standard basis is ', standard_posterior)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "And we're done! Enjoy the XYZ basis."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": []
+  }
+ ],
+ "metadata": {
+  "interpreter": {
+   "hash": "cbafa42b6ff5e39c327df57f5a3bbaee7fd2ec257998acb6b9d49a010b26573d"
+  },
+  "kernelspec": {
+   "display_name": "Python 3.8.3  ('base': conda)",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": ""
+  },
+  "orig_nbformat": 4
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
\ No newline at end of file
diff --git a/orbitize/conversions.py b/orbitize/conversions.py
new file mode 100644
index 00000000..e27ba73d
--- /dev/null
+++ b/orbitize/conversions.py
@@ -0,0 +1,249 @@
+import numpy as np
+from astropy import units as u
+from astropy import constants as consts
+from orbitize import kepler
+
+"""
+This module converts orbital elements between different coordinate bases.
+"""
+
+def ecosomega_to_standard(elems):
+    """
+    Converts array of orbital elements with ecc*cos aop and ecc*sin aop to the usual one with ecc and aop.
+
+    Args:
+        elems (np.array of floats): Orbital elements with e cosine omega and e sin omega (semi-major axis [au], ecc*cos aop,
+            inclination [radians], ecc*sin aop, longitude of the ascending node [radians], epoch of 
+            periastron passage in fraction of orbital period past MJD=0 [0,1], parallax [mas], total mass of the two-body orbit
+            (M_* + M_planet) [Solar masses]). If more than 1 set of parameters is passed, the first dimension must be
+            the number of orbital parameter sets, and the second the orbital elements.
+
+    Return:
+        np.array: Orbital elements with ecc and aop
+    """
+    if elems.ndim == 1:
+        elems = elems[np.newaxis, :]
+    ecos, esin = elems[:,1], elems[:,3]
+    # Solve for aop, in the interval [-pi, pi]
+    aop = np.arctan2(esin, ecos)
+    # Fix interval to [0, 2 pi]
+    aop[aop < 0.0] = 2*np.pi + aop[aop < 0.0]
+    #Solve for ecc, taking care of possible negative cosines
+    ecc = np.sqrt(ecos**2 / (np.cos(aop)**2))
+
+    elems[:,1], elems[:,3] = ecc, aop
+
+    return np.squeeze(elems)
+
+def standard_to_ecosomega(elems):
+    """
+    Converts array of orbital elements from the usual basis to ecc*cos aop and ecc*sin aop.
+
+    Args:
+        elems (np.array of floats): Orbital elements with e cosine omega and e sin omega (semi-major axis [au], ecc,
+            inclination [radians], aop, longitude of the ascending node [radians], epoch of 
+            periastron passage in fraction of orbital period past MJD=0 [0,1], parallax [mas], total mass of the two-body orbit
+            (M_* + M_planet) [Solar masses]). If more than 1 set of parameters is passed, the first dimension must be
+            the number of orbital parameter sets, and the second the orbital elements.
+
+    Return:
+        np.array: Orbital elements with ecc*cos aop and ecc*sin aop at indices 1 and 3.
+    """
+    if elems.ndim == 1:
+        elems = elems[np.newaxis, :]
+    ecc, aop = elems[:,1], elems[:,3]
+    e_sin, e_cos = ecc*np.sin(aop), ecc*np.cos(aop)
+
+    elems[:,1], elems[:,3] = e_cos, e_sin
+
+    return np.squeeze(elems)
+
+def xyz_to_standard(epoch, elems, tau_ref_epoch=58849):
+    """
+    Converts array of orbital elements in terms of position and velocity in xyz to the regular one
+
+    Args:
+        epoch (float): Date in MJD of observation to calculate time of periastron passage (tau).
+        elems (np.array of floats): Orbital elements in xyz basis (x-coordinate [au], y-coordinate [au],
+            z-coordinate [au], velocity in x [km/s], velocity in y [km/s], velocity in z [km/s], parallax [mas], total mass of the two-body orbit
+            (M_* + M_planet) [Solar masses]). If more than 1 set of parameters is passed, the first dimension must be
+            the number of orbital parameter sets, and the second the orbital elements.
+
+    Return:
+        np.array: Orbital elements in the usual basis (sma, ecc, inc, aop, pan, tau, plx, mtot)
+    """
+
+    if elems.ndim == 1:
+        elems = elems[np.newaxis, :]
+    # Velocities and positions, with units
+    vel = elems[:,3:6] * u.km / u.s # velocities in km / s ?
+    pos = elems[:,0:3] * u.AU # positions in AU ?
+    vel_magnitude = np.linalg.norm(vel, axis=1)
+    pos_magnitude = np.linalg.norm(pos, axis=1)
+
+    # Mass
+    mtot = elems[:,7]*u.Msun
+    mu = consts.G * mtot # G in m3 kg-1 s-2, mtot in msun
+
+    # Angular momentum, making sure nodal vector is not exactly zero
+    h = (np.cross(pos, vel, axis=1)).si
+    # if h[0].value == 0.0 and h[1].value == 0.0:
+    #     pos[2] = 1e-8*u.AU
+    #     h = (np.cross(pos, vel)).si
+    h_magnitude = np.linalg.norm(h, axis=1)
+
+    sma = 1 / (2.0 / pos_magnitude - (vel_magnitude**2)/mu)
+    sma = sma.to(u.AU)
+
+    ecc = (np.sqrt(1 - h_magnitude**2 / (sma * mu))).value
+    e_vector = (np.cross(vel, h, axis=1) / mu[:, None] - pos / pos_magnitude[:, None]).si
+    e_vec_magnitude = np.linalg.norm(e_vector, axis=1)
+
+    unit_k = np.array((0,0,1))
+    cos_inc = (np.dot(h, unit_k) / h_magnitude).value
+    inc = np.arccos(-cos_inc)
+
+    #Nodal vector
+    n = np.cross(unit_k, h)
+    n_magnitude = np.linalg.norm(n, axis=1)
+
+    # Position angle of the nodes, checking for the right quadrant
+    # np.arccos yields angles in [0, pi]
+    unit_i = np.array((1,0,0))
+    unit_j = np.array((0,1,0))
+    cos_pan = (np.dot(n, unit_j) / n_magnitude).value
+    pan = np.arccos(cos_pan)
+    n_x = np.dot(n, unit_i)
+    pan[n_x < 0.0] = 2*np.pi - pan[n_x < 0.0]
+
+    # Argument of periastron, checking for the right quadrant
+    cos_aop = (np.sum(n*e_vector, axis=1) / (n_magnitude * e_vec_magnitude)).value
+    aop = np.arccos(cos_aop)
+    e_vector_z = np.dot(e_vector, unit_k)
+    aop[e_vector_z < 0.0] = 2.0*np.pi - aop[e_vector_z < 0.0]
+
+    # True anomaly, checking for the right quadrant
+    cos_tanom = (np.sum(pos*e_vector, axis=1) / (pos_magnitude*e_vec_magnitude)).value
+    tanom = np.arccos(cos_tanom)
+    # Check for places where tanom is nan, due to cos_tanom=1. (for some reason that was a problem)
+    tanom = np.where((0.9999<cos_tanom ) & (cos_tanom<1.001), 0.0, tanom)
+    rdotv = np.sum(pos*vel, axis=1)
+    tanom[rdotv < 0.0] = 2*np.pi - tanom[rdotv < 0.0]
+
+    # Eccentric anomaly to get tau, checking for the right quadrant
+    cos_eanom = ((1 - pos_magnitude / sma) / ecc).value
+    eanom = np.arccos(cos_eanom)
+    # Check for places where eanom is nan, due to cos_eanom = 1.(same problem as above)
+    eanom = np.where((0.9999<cos_eanom ) & (cos_eanom<1.001), 0.0, eanom)
+    eanom[tanom > np.pi] =  2*np.pi - eanom[tanom > np.pi]
+
+    # Time of periastron passage, using Kepler's equation, in MJD:
+    time_tau = epoch - ((np.sqrt(sma**3 / mu)).to(u.day)).value * (eanom - ecc*np.sin(eanom))
+
+    # Calculate period from Kepler's third law, in days:
+    period = np.sqrt(4*np.pi**2.0*(sma)**3/mu)
+    period = period.to(u.day).value
+
+    # Finally, tau
+    tau = (time_tau - tau_ref_epoch) / period
+    tau = tau%1.0
+
+    mtot = mtot.value
+    sma = sma.value
+
+    results = np.array((sma, ecc, inc, aop, pan, tau, elems[:,6], mtot)).T
+
+    return np.squeeze(results)
+
+def standard_to_xyz(epoch, elems, tau_ref_epoch=58849, tolerance=1e-9, max_iter=100):
+    """
+    Converts array of orbital elements from the regular base of Keplerian orbits to positions and velocities in xyz
+    Uses code from orbitize.kepler
+
+    Args:
+        epoch (float): Date in MJD of observation to calculate time of periastron passage (tau).
+        elems (np.array of floats): Orbital elements (sma, ecc, inc, aop, pan, tau, plx, mtot).
+                If more than 1 set of parameters is passed, the first dimension must be
+                the number of orbital parameter sets, and the second the orbital elements.
+
+    Return:
+        np.array: Orbital elements in xyz (x-coordinate [au], y-coordinate [au], z-coordinate [au], 
+        velocity in x [km/s], velocity in y [km/s], velocity in z [km/s], parallax [mas], total mass of the two-body orbit
+            (M_* + M_planet) [Solar masses])
+    """
+
+    # Use classical elements to obtain position and velocity in the perifocal coordinate system
+    # Then transform coordinates using matrix multiplication
+
+    if elems.ndim == 1:
+        elems = elems[np.newaxis, :]
+    # Define variables
+    sma = elems[:,0] # AU
+    ecc = elems[:,1] # [0.0, 1.0]
+    inc = elems[:,2] # rad [0, pi]
+    aop = elems[:,3] # rad [0, 2 pi]
+    pan = elems[:,4] # rad [0, 2 pi]
+    tau = elems[:,5] # [0.0, 1.0]
+    mtot = elems[:,7] # Msun
+
+    # Just in case so nothing breaks
+    ecc = np.where(ecc == 0.0, 1e-8, ecc)
+    inc = np.where(inc == 0.0, 1e-8, inc)
+
+    # Begin by calculating the eccentric anomaly
+    period = np.sqrt(4*np.pi**2.0*(sma*u.AU)**3/(consts.G*(mtot*u.Msun)))
+    period = period.to(u.day).value # Period in days
+    mean_motion = 2*np.pi/(period)
+
+    # Mean anomaly
+    manom = (mean_motion*(epoch - tau_ref_epoch) - 2*np.pi*tau) % (2.0*np.pi)
+    # Eccentric anomaly
+    eanom = kepler._calc_ecc_anom(manom, ecc, tolerance=tolerance, max_iter=max_iter)
+    # if eanom.ndim == 1:
+    #     eanom = eanom[np.newaxis, :]
+    # Magnitude of angular momentum:
+    h = np.sqrt(consts.G*(mtot*u.Msun)*(sma*u.AU)*(1 - ecc**2))
+
+    # Position vector in the perifocal system in AU
+    pos_peri_x = (sma*(np.cos(eanom) - ecc))
+    pos_peri_y = (sma*np.sqrt(1 - ecc**2)*np.sin(eanom))
+    pos_peri_z = np.zeros(len(elems))
+
+    pos = np.stack((pos_peri_x, pos_peri_y, pos_peri_z)).T
+    pos_magnitude = np.linalg.norm(pos, axis=1)
+
+    # Velocity vector in the perifocal system in km/s
+    vel_peri_x = - ((np.sqrt(consts.G*(mtot*u.Msun)*(sma*u.AU))*np.sin(eanom) / (pos_magnitude*u.AU)).to(u.km / u.s)).value 
+    vel_peri_y = ((h* np.cos(eanom) / (pos_magnitude*u.AU)).to(u.km / u.s)).value
+    vel_peri_z = np.zeros(len(elems))
+
+    vel = np.stack((vel_peri_x, vel_peri_y, vel_peri_z)).T
+
+    # Transformation matrix to inertial xyz system, component by component
+    pan = pan +np.pi / 2.0
+    T_11 = np.cos(pan)*np.cos(aop) - np.sin(pan)*np.sin(aop)*np.cos(inc)
+    T_12 = - np.cos(pan)*np.sin(aop) - np.sin(pan)*np.cos(aop)*np.cos(inc)
+    T_13 = np.sin(pan)*np.sin(inc)
+
+    T_21 = np.sin(pan)*np.cos(aop) + np.cos(pan)*np.sin(aop)*np.cos(inc)
+    T_22 = - np.sin(pan)*np.sin(aop) + np.cos(pan)*np.cos(aop)*np.cos(inc)
+    T_23 = - np.cos(pan)*np.sin(inc)
+
+    T_31 = np.sin(aop)*np.sin(inc)
+    T_32 = np.cos(aop)*np.sin(inc)
+    T_33 = np.cos(inc)
+
+    T = np.array([[T_11, T_12, T_13],
+                  [T_21, T_22, T_23],
+                  [T_31, T_32, T_33]])
+
+    pos_xyz = np.zeros((len(elems), 3))
+    vel_xyz = np.zeros((len(elems), 3))
+    for k in range(len(elems)):
+        pos_xyz[k,:] =  np.matmul(T[:,:,k], pos[k])
+        vel_xyz[k,:] =  np.matmul(T[:,:,k], vel[k])
+
+    # Flipping x-axis sign to increase X as RA increases
+    result = np.stack([-pos_xyz[:,0], pos_xyz[:,1], pos_xyz[:,2], -vel_xyz[:,0], vel_xyz[:,1], vel_xyz[:,2], elems[:,6], mtot]).T
+
+    return np.squeeze(result)
diff --git a/orbitize/driver.py b/orbitize/driver.py
index d9f80fc2..0cc91405 100644
--- a/orbitize/driver.py
+++ b/orbitize/driver.py
@@ -36,7 +36,7 @@ class Driver(object):
     def __init__(self, input_data, sampler_str,
                  num_secondary_bodies, system_mass, plx,
                  mass_err=0, plx_err=0, lnlike='chi2_lnlike',
-                 system_kwargs=None, mcmc_kwargs=None):
+                 system_kwargs=None, mcmc_kwargs=None, xyz_kwargs=None, fitting_basis='standard'):
 
         # Read in data
         # Try to interpret input as a filename first
@@ -53,12 +53,15 @@ def __init__(self, input_data, sampler_str,
         if system_kwargs is None:
             system_kwargs = {}
 
+        if fitting_basis == 'XYZ' and xyz_kwargs is not None:
+            system_kwargs.update(xyz_kwargs)
+
         if sampler_str == 'OFTI' and ('fit_secondary_mass' in system_kwargs and True == system_kwargs['fit_secondary_mass']):
             raise Exception('Run Astrometry+RV in MCMC for now.')
         # Initialize System object which stores data & sets priors
         self.system = orbitize.system.System(
             num_secondary_bodies, data_table, system_mass,
-            plx, mass_err=mass_err, plx_err=plx_err, **system_kwargs
+            plx, mass_err=mass_err, plx_err=plx_err, fitting_basis=fitting_basis, **system_kwargs
         )
 
         # Initialize Sampler object, which has System object as an attribute.
diff --git a/orbitize/example_data/HIP027321.d b/orbitize/example_data/HIP027321.d
new file mode 100644
index 00000000..0cadfe60
--- /dev/null
+++ b/orbitize/example_data/HIP027321.d
@@ -0,0 +1,112 @@
+ 27321  27251 111 1   5    0  -1.81  0           
+ 133 -1.245  0.624 -0.9065 -0.4222   -1.00   0.81
+ 133 -1.245  0.626 -0.9050 -0.4254   -0.31   0.80
+ 133 -1.245  0.622 -0.9076 -0.4198   -1.66   0.78
+ 194 -1.170 -0.651 -0.0680  0.9977    0.39   0.78
+ 194 -1.170 -0.649 -0.0716  0.9974   -0.79   0.86
+ 194 -1.170 -0.648 -0.0723  0.9974   -0.64   0.87
+ 257 -1.094  0.671 -0.8760  0.4823    0.01   0.80
+ 257 -1.094  0.674 -0.8778  0.4790   -0.23   0.78
+ 257 -1.094  0.671 -0.8756  0.4830   -0.20   0.85
+ 327 -1.009 -0.679  0.8136  0.5814    0.19   0.75
+ 327 -1.009 -0.680  0.8147  0.5799    0.02   0.82
+ 327 -1.009 -0.679  0.8137  0.5813   -1.26   0.81
+ 458 -0.849 -0.641  0.9457 -0.3250    0.51   0.91
+ 458 -0.849 -0.641  0.9457 -0.3250   -1.01   0.82
+ 458 -0.849 -0.640  0.9461 -0.3240    0.25   1.06
+ 458 -0.849 -0.641  0.9457 -0.3250   -1.17   0.84
+ 458 -0.849 -0.640  0.9458 -0.3247   -1.04   0.80
+ 503 -0.794  0.637  0.7353  0.6777    0.75   0.95
+ 503 -0.794  0.636  0.7367  0.6762   -0.01   0.82
+ 503 -0.794  0.636  0.7365  0.6764   -0.51   1.10
+ 503 -0.794  0.636  0.7365  0.6765   -0.34   0.96
+ 585 -0.695 -0.656  0.2592 -0.9658   -0.22   0.70
+ 585 -0.695 -0.655  0.2602 -0.9656   -0.52   0.86
+ 585 -0.695 -0.655  0.2609 -0.9654    0.73   0.81
+ 634 -0.636  0.674  0.9777 -0.2102   -0.19   0.80
+ 709 -0.544 -0.687 -0.6238 -0.7816    0.02   0.67
+ 709 -0.544 -0.685 -0.6225 -0.7826    0.09   0.73
+ 709 -0.544 -0.685 -0.6216 -0.7833   -0.41   0.72
+ 709 -0.544 -0.686 -0.6226 -0.7826   -0.21   0.76
+ 770 -0.470  0.673  0.3654 -0.9308   -1.15   0.84
+ 770 -0.470  0.672  0.3640 -0.9314   -0.04   0.84
+ 834 -0.392 -0.640 -0.9999  0.0141    0.66   0.78
+ 834 -0.392 -0.643 -0.9998  0.0182   -0.21   0.81
+ 834 -0.392 -0.640 -0.9999  0.0145   -1.21   0.84
+ 834 -0.392 -0.642 -0.9998  0.0177    0.20   0.76
+ 904 -0.307  0.617 -0.6578 -0.7532    0.20   0.71
+ 904 -0.307  0.617 -0.6584 -0.7526    0.05   0.74
+ 964 -0.234 -0.625 -0.4897  0.8719   -0.86   0.78
+ 964 -0.234 -0.623 -0.4919  0.8706    0.36   1.24
+ 964 -0.234 -0.623 -0.4930  0.8700   -1.03   0.86
+ 964 -0.234 -0.628 -0.4859  0.8740   -0.30   0.86
+ 964 -0.234 -0.625 -0.4904  0.8715    0.82   0.83
+1031 -0.153  0.655 -0.9905  0.1376    0.10   0.73
+1031 -0.153  0.651 -0.9899  0.1420   -0.37   0.88
+1031 -0.153  0.654 -0.9903  0.1389   -0.49   0.84
+1031 -0.153  0.652 -0.9901  0.1404   -0.85   0.84
+1097 -0.073 -0.675  0.5283  0.8491   -2.09   1.11
+1097 -0.073 -0.679  0.5330  0.8461   -0.90   1.08
+1097 -0.073 -0.677  0.5301  0.8479   -0.25   0.89
+1097 -0.073 -0.677  0.5305  0.8477    0.14   1.16
+1097 -0.073 -0.678  0.5318  0.8469   -0.05   0.80
+1153 -0.005  0.675 -0.4994  0.8664    0.04   0.83
+1153 -0.005  0.676 -0.5006  0.8657    0.21   0.91
+1153 -0.005  0.675 -0.4995  0.8663    0.73   1.15
+1153 -0.005  0.675 -0.4995  0.8663   -1.65   0.98
+1229  0.088 -0.653  0.9986  0.0528    0.30   0.86
+1229  0.088 -0.654  0.9987  0.0518   -0.53   1.04
+1229  0.088 -0.654  0.9987  0.0513   -0.11   1.00
+1229  0.088 -0.654  0.9987  0.0516   -0.84   0.83
+1275  0.144  0.636  0.4213  0.9069    0.34   0.85
+1276  0.145  0.644  0.4193  0.9078   -0.45   0.77
+1404  0.300  0.657  0.9862  0.1654   -0.43   1.06
+1404  0.300  0.656  0.9864  0.1643    1.55   0.96
+1404  0.300  0.657  0.9861  0.1661    0.71   0.90
+1404  0.300  0.657  0.9861  0.1659   -0.75   0.79
+1404  0.300  0.658  0.9860  0.1665   -2.69   1.19
+1404  0.300  0.658  0.9861  0.1664   -1.40   1.09
+1482  0.395 -0.683 -0.3021 -0.9533    0.38   0.85
+1482  0.395 -0.683 -0.3023 -0.9532   -1.30   1.06
+1482  0.395 -0.682 -0.3009 -0.9536    0.11   0.83
+1482  0.395 -0.683 -0.3023 -0.9532    1.19   1.15
+1539  0.464  0.690  0.7038 -0.7104   -0.50   0.85
+1539  0.464  0.691  0.7044 -0.7098    1.08   1.22
+1539  0.464  0.689  0.7033 -0.7109   -0.49   0.88
+1539  0.464  0.689  0.7027 -0.7115    1.45   1.07
+1607  0.547 -0.664 -0.9372 -0.3487   -1.27   1.23
+1607  0.547 -0.665 -0.9374 -0.3483   -0.21   0.93
+1674  0.628  0.631 -0.2787 -0.9604   -1.03   1.07
+1674  0.628  0.634 -0.2749 -0.9615   -1.85   1.17
+1674  0.628  0.633 -0.2756 -0.9613    0.30   1.21
+1674  0.628  0.631 -0.2785 -0.9604   -0.98   0.83
+1734  0.701 -0.623 -0.8109  0.5853    1.93   0.77
+1734  0.701 -0.620 -0.8130  0.5823   -0.85   1.18
+1734  0.701 -0.622 -0.8121  0.5835   -0.05   1.19
+1735  0.702 -0.621 -0.8083  0.5888   -1.82   1.18
+1804  0.786  0.630 -0.9717 -0.2363    0.09   0.77
+1804  0.786  0.630 -0.9715 -0.2369   -0.34   0.77
+1867  0.863 -0.664  0.1487  0.9889   -0.58   0.76
+2000  1.024 -0.669  0.9088  0.4172   -0.23   0.91
+2049  1.084  0.650  0.0446  0.9990    0.04   0.73
+2049  1.084  0.650  0.0441  0.9990   -0.02   0.77
+2049  1.084  0.650  0.0440  0.9990   -0.24   0.72
+2049  1.084  0.650  0.0441  0.9990    0.38   0.72
+2130  1.182 -0.633  0.8661 -0.4999   -0.27   0.70
+2130  1.182 -0.633  0.8659 -0.5002   -0.17   0.91
+2130  1.182 -0.632  0.8661 -0.4998   -0.16   0.90
+2130  1.182 -0.632  0.8661 -0.4999   -1.14   0.98
+2175  1.237  0.640  0.8524  0.5229   -0.59   1.30
+2175  1.237  0.640  0.8523  0.5231   -1.31   0.97
+2175  1.237  0.638  0.8536  0.5210    0.17   1.05
+2175  1.237  0.638  0.8537  0.5207    0.88   1.18
+2256  1.335 -0.663  0.0709 -0.9975   -0.22   1.26
+2256  1.335 -0.664  0.0702 -0.9975    0.26   0.97
+2506  1.639 -0.632 -0.9775  0.2107   -0.51   0.76
+2506  1.639 -0.628 -0.9787  0.2055   -0.17   0.92
+2506  1.639 -0.630 -0.9781  0.2081    0.04   0.84
+2506  1.639 -0.633 -0.9771  0.2128    0.39   0.89
+2506  1.639 -0.628 -0.9786  0.2058   -0.28   0.81
+2677  1.846  0.652 -0.9908  0.1352   -0.03   0.85
+2677  1.846  0.651 -0.9904  0.1380   -2.97   1.63
+2677  1.846  0.651 -0.9906  0.1364   -0.69   0.83
diff --git a/orbitize/example_data/betaPic.csv b/orbitize/example_data/betaPic.csv
new file mode 100644
index 00000000..918ce3d1
--- /dev/null
+++ b/orbitize/example_data/betaPic.csv
@@ -0,0 +1,47 @@
+epoch,object,sep,sep_err,pa,pa_err
+54781,1, 0.210,0.027,211.49,1.9
+52953,1,0.413,0.022,34,4
+55129,1, 0.299,0.014,211,3
+55194,1, 0.306,0.009,212.1,1.7
+55296,1, 0.346,0.007,209.9,1.2
+55467,1, 0.383,0.011,210.3,1.7
+55516,1, 0.387,0.008,212.4,1.4
+55517,1,0.390,0.013,212,2
+55593,1, 0.408,0.009,211.1,1.5
+55646,1, 0.426,0.013,210.1,1.8
+55168,1, 0.339,0.010,209.2,1.7
+55168,1, 0.323,0.010,209.3,1.8
+55555,1, 0.407,0.005,212.8,1.4
+55854,1, 0.452,0.003,211.6,0.4
+55854,1, 0.455,0.005,211.9,0.6
+56015,1,0.447,0.003,210.8,0.4
+56015,1, 0.448,0.005,211.8,0.6
+56263,1, 0.461,0.014,211.9,1.2
+56265,1, 0.470,0.010,212.0,1.2
+56612,1, 0.4308,0.0015,212.43,0.17
+56612,1, 0.4291,0.0010,212.58,0.15
+56614,1, 0.4302,0.0010,212.46,0.15
+56636,1, 0.4255,0.0010,212.51,0.15
+56636,1,0.4244,0.0010,212.85,0.15
+56637,1, 0.4253,0.0010,212.47,0.16
+56969,1, 0.3562,0.0010,213.02,0.19
+57114,1, 0.3173,0.0009,213.13,0.20
+57332,1, 0.2505,0.0015,214.14,0.34
+57361,1, 0.2402,0.0011,213.58,0.34
+57378,1, 0.2345,0.0010,213.81,0.30
+57408,1, 0.2226,0.0021,214.84,0.44
+56999,1,0.35051,0.00320,212.60,0.66
+57147,1, 0.33242,0.00170,212.58,0.35
+57296,1, 0.26202,0.00178,213.02,0.48
+57356,1, 0.24205,0.00251,213.30,0.74
+57382,1, 0.23484,0.00180,213.79,0.51
+57407,1, 0.22723,0.00155,213.15,0.46
+57473,1, 0.20366,0.00142,213.90,0.46
+57494,1, 0.19749,0.00236,213.88,0.83
+57647,1,0.14236,0.00234,214.62,1.10
+57675,1, 0.13450,0.00246,215.50,1.22
+57710,1, 0.12712,0.00644,215.80,3.37
+58378,1, 0.14046,0.00312,29.71,1.67
+57046,1, 0.3355,0.0009,212.88,0.20	
+58382,1, 0.1419,0.0053,28.16,1.82	
+58440,1,0.1645,0.0018,28.64,0.70
\ No newline at end of file
diff --git a/orbitize/hipparcos.py b/orbitize/hipparcos.py
new file mode 100644
index 00000000..01bdd383
--- /dev/null
+++ b/orbitize/hipparcos.py
@@ -0,0 +1,146 @@
+import numpy as np
+import pandas as pd
+
+from orbitize.kepler import calc_orbit
+from orbitize.radvel_utils.compute_sep import compute_sep
+from astroquery.vizier import Vizier
+from astropy.time import Time
+from astropy.coordinates import get_body_barycentric_posvel
+
+class HipparcosLogProb(object):
+    """
+    Class to compute the log probability of an orbit with respect to the 
+    Hipparcos Intermediate Astrometric Data (IAD). Queries Vizier for 
+    all metadata relevant to the IAD, and reads in the IAD themselves from
+    a specified location. Follows Nielsen+ 2020
+    (studying the orbit of beta Pic b).
+
+    Args:
+        iad_file (str): 
+        hip_num (str): the Hipparcos number of your target. Accessible on Simbad.
+    """
+
+    def __init__(self, iad_file, hip_num, num_secondary_bodies):
+
+        self.hip_num = hip_num
+        self.num_secondary_bodies = num_secondary_bodies
+
+        # load best-fit astrometric solution from van Leeuwen catalog
+        Vizier.ROW_LIMIT = -1
+        hip_cat = Vizier(
+            catalog='I/311/hip2', 
+            columns=[
+                'RArad', 'e_RArad', 'DErad', 'e_DErad', 'Plx', 'e_Plx', 'pmRA', 
+                'e_pmRA', 'pmDE', 'e_pmDE', 'F2'
+            ]
+        ).query_constraints(HIP=self.hip_num)[0]
+
+        self.plx0 = hip_cat['Plx'][0] # [mas]
+        self.pm_ra0 = hip_cat['pmRA'][0] # [mas/yr]
+        self.pm_dec0 = hip_cat['pmDE'][0] # [mas/yr]
+        self.alpha0 = hip_cat['RArad'][0] # [deg]
+        self.delta0 = hip_cat['DErad'][0] # [deg]
+        self.plx0_err = hip_cat['e_Plx'][0] # [mas]
+        self.pm_ra0_err = hip_cat['e_pmRA'][0] # [mas/yr]
+        self.pm_dec0_err = hip_cat['e_pmDE'][0] # [mas/yr]
+        self.alpha0_err = hip_cat['e_RArad'][0] # [mas]
+        self.delta0_err = hip_cat['e_DErad'][0] # [mas]
+
+        # read in IAD
+        iad = np.transpose(np.loadtxt(iad_file, skiprows=1))
+
+        times = iad[1] + 1991.25
+        epochs = Time(times, format='decimalyear')
+        self.epochs = epochs.decimalyear
+        self.epochs_mjd = epochs.mjd
+        self.cos_phi = iad[3] # scan direction
+        self.sin_phi = iad[4]
+        self.R = iad[5] # abscissa residual [mas]
+        self.eps = iad[6] # error on abscissa residual [mas]
+
+        # compute Earth XYZ position in barycentric coordinates
+        bary_pos, _ = get_body_barycentric_posvel('earth', epochs)
+        self.X = bary_pos.x.value # [au]
+        self.Y = bary_pos.y.value # [au]
+        self.Z = bary_pos.z.value # [au]
+
+        # reconstruct ephemeris of star given van Leeuwen best-fit (Nielsen+ 2020 Eqs 1-2) [mas]
+        changein_alpha_st = (
+            self.plx0 * (
+                self.X * np.sin(np.radians(self.alpha0)) - 
+                self.Y * np.cos(np.radians(self.alpha0))
+            ) + (self.epochs - 1991.25) * self.pm_ra0
+        )
+
+        changein_delta = (
+            self.plx0 * (
+                self.X * np.cos(np.radians(self.alpha0)) * np.sin(np.radians(self.delta0)) + 
+                self.Y * np.sin(np.radians(self.alpha0)) * np.sin(np.radians(self.delta0)) - 
+                self.Z * np.cos(np.radians(self.delta0))
+            ) + (self.epochs - 1991.25) * self.pm_dec0
+        )
+
+        # compute abcissa point (Nielsen+ Eq 3)
+        self.alpha_abs_st = self.R * self.cos_phi + changein_alpha_st
+        self.delta_abs = self.R * self.sin_phi + changein_delta
+
+    def compute_lnlike(self, raoff_model, deoff_model, samples, negative=False):
+        """
+        Computes the log probability of an orbit model with respect to the Hipparcos 
+        IAD. 
+
+        raoff/deoff: ra/deoffsets for star, computed for an arbitrary orbit model
+        astr_samples: pm_ra, pm_dec, alpha_H0, delta_H0, plx (fitted astrometric params)
+
+        Returns:
+            np.array of length M, where M is the number of input orbits (same as def'n
+                in description of `samples` arg above) representing the log probability
+                for each orbit with respect to the Hipparcos IAD
+        """
+        n_params = len(samples)
+
+        # variables for each of the astrometric fitting parameters
+        plx = samples[6 * self.num_secondary_bodies]
+        pm_ra = samples[6 * self.num_secondary_bodies + 1]
+        pm_dec = samples[6 * self.num_secondary_bodies + 2]
+        alpha_H0 = samples[6 * self.num_secondary_bodies + 3]
+        delta_H0 = samples[6 * self.num_secondary_bodies + 4]
+
+
+        try:
+            n_samples = len(pm_ra)
+        except TypeError:
+            n_samples = 1
+
+        n_epochs = len(self.epochs)
+        dist = np.empty((n_epochs, n_samples))
+
+        # add parallactic ellipse & proper motion to position (Nielsen+ 2020 Eq 8)
+        for i in np.arange(n_epochs):
+
+            # this is the expected offset from the Hipparcos photocenter in 1991.25
+            alpha_C_st = alpha_H0 + plx * (
+                self.X[i] * np.sin(np.radians(self.alpha0)) - 
+                self.Y[i] * np.cos(np.radians(self.alpha0))
+            ) + (self.epochs[i] - 1991.25) * pm_ra
+            delta_C = delta_H0 + plx * (
+                self.X[i] * np.cos(np.radians(self.alpha0)) * np.sin(np.radians(self.delta0)) + 
+                self.Y[i] * np.sin(np.radians(self.alpha0)) * np.sin(np.radians(self.delta0)) -
+                self.Z[i] * np.cos(np.radians(self.delta0))
+            ) + (self.epochs[i] - 1991.25) * pm_dec
+
+            # add in pre-computed secondary perturbations
+            alpha_C_st += raoff_model[i]
+            delta_C += deoff_model[i]
+
+            # calculate distance between line and expected measurement (Nielsen+ 2020 Eq 6) [mas]
+            dist[i, :] = np.abs(
+                (self.alpha_abs_st[i] - alpha_C_st) * self.cos_phi[i] + 
+                (self.delta_abs[i] - delta_C) * self.sin_phi[i]
+            )
+
+        # compute chi2 (Nielsen+ 2020 Eq 7)
+        chi2 = np.sum([(dist[:,i] / self.eps)**2 for i in np.arange(n_samples)], axis=1)
+        lnlike = -0.5 * chi2
+
+        return lnlike
\ No newline at end of file
diff --git a/orbitize/results.py b/orbitize/results.py
index ec12bfb1..ef9899ee 100644
--- a/orbitize/results.py
+++ b/orbitize/results.py
@@ -1,3 +1,4 @@
+from plistlib import Data
 import numpy as np
 import warnings
 import h5py
@@ -65,7 +66,8 @@ class Results(object):
     """
 
     def __init__(self, sampler_name=None, post=None, lnlike=None, tau_ref_epoch=None, labels=None,
-                 data=None, num_secondary_bodies=None, version_number=None, curr_pos=None):
+                 data=None, num_secondary_bodies=None, version_number=None, curr_pos=None, fitting_basis='standard', xyz_epochs=None):
+
 
         self.sampler_name = sampler_name
         self.post = post
@@ -80,6 +82,9 @@ def __init__(self, sampler_name=None, post=None, lnlike=None, tau_ref_epoch=None
         self.num_secondary_bodies=num_secondary_bodies
         self.curr_pos = curr_pos
         self.version_number = version_number
+        self.fitting_basis = fitting_basis
+        self.xyz_epochs = xyz_epochs
+
 
     def add_samples(self, orbital_params, lnlikes, labels, curr_pos=None):
         """
@@ -145,6 +150,10 @@ def save_results(self, filename):
         hf.attrs['sampler_name'] = self.sampler_name
         hf.attrs['tau_ref_epoch'] = self.tau_ref_epoch
         hf.attrs['version_number'] = self.version_number
+        hf.attrs['fitting_basis'] = self.fitting_basis
+        if self.fitting_basis == 'XYZ':
+            hf.create_dataset('xyz_epochs', data= self.xyz_epochs)
+            
         # Now add post and lnlike from the results object as datasets
         hf.create_dataset('post', data=self.post)
         hf.create_dataset('data', data=self.data)
@@ -212,12 +221,24 @@ def load_results(self, filename, append=False):
         except KeyError:
             curr_pos = None
 
+        try:
+            fitting_basis = np.str(hf.attrs['fitting_basis'])
+        except KeyError:
+            # if key does not exist, then it was fit in the standard basis
+            fitting_basis == 'standard'
+        try:
+            xyz_epochs = np.array(hf.get('xyz_epochs'))
+        except KeyError:
+            # if KeyError, this was not fit in xyz
+            xyz_epochs = None
+
         hf.close()  # Closes file object
 
         # doesn't matter if append or not. Overwrite curr_pos if not None
         if curr_pos is not None:
             self.curr_pos = curr_pos
 
+        # TODO: Check if this part is consistent with xyz_epoch
         # Adds loaded data to object as per append keyword
         if append:
             # if no sampler_name set, use the input file's value
@@ -263,6 +284,8 @@ def load_results(self, filename, append=False):
                 self.add_samples(post, lnlike, self.labels)
                 self.tau_ref_epoch = tau_ref_epoch
                 self.num_secondary_bodies = num_secondary_bodies
+                self.fitting_basis = fitting_basis
+                self.xyz_epochs = xyz_epochs
             else:
                 raise Exception(
                     'Unable to load file {} to Results object. append is set to False but object is not empty'.format(filename))
@@ -316,6 +339,10 @@ def plot_corner(self, param_list=None, **corner_kwargs):
             'mtot': '$M_T$ [M$_{{\\odot}}$]',
             'm0': '$M_0$ [M$_{{\\odot}}$]',
             'm': '$M_{0}$ [M$_\{{Jup\}}$]',
+            'pm_ra': 'PM RA',
+            'pm_dec': 'PM Dec',
+            'alpha0': 'alpha0',
+            'delta0': 'delta0',
         }
 
         if param_list is None:
@@ -336,11 +363,11 @@ def plot_corner(self, param_list=None, **corner_kwargs):
                 if label_key.startswith('m') and label_key != 'm0' and label_key != 'mtot':
                     secondary_mass_indices.append(i)
 
-
         samples = np.copy(self.post[:, param_indices])  # keep only chains for selected parameters
         samples[:, angle_indices] = np.degrees(
             samples[:, angle_indices])  # convert angles from rad to deg
         samples[:, secondary_mass_indices] *= u.solMass.to(u.jupiterMass) # convert to Jupiter masses for companions
+        samples = samples[:, param_indices]
 
         if 'labels' not in corner_kwargs:  # use default labels if user didn't already supply them
             reduced_labels_list = []
@@ -352,6 +379,8 @@ def plot_corner(self, param_list=None, **corner_kwargs):
                 elif label_key == 'm0' or label_key == 'mtot' or label_key.startswith('plx'):
                     body_num = ""
                     # maintain original label key
+                elif label_key in ['pm_ra', 'pm_dec', 'alpha0', 'delta0']:
+                    body_num = ""
                 else:
                     body_num = label_key[3]
                     label_key = label_key[0:3]
diff --git a/orbitize/sampler.py b/orbitize/sampler.py
index 2fdfdbbe..c593436d 100644
--- a/orbitize/sampler.py
+++ b/orbitize/sampler.py
@@ -5,10 +5,12 @@
 import abc
 import math
 import time
+import orbitize.conversions
 
 import emcee
 import ptemcee
 import multiprocessing as mp
+from multiprocessing import Pool
 
 import orbitize.lnlike
 import orbitize.priors
@@ -44,7 +46,7 @@ def __init__(self, system, like='chi2_lnlike', custom_lnlike=None):
     def run_sampler(self, total_orbits):
         pass
 
-    def _logl(self, params):
+    def _logl(self, params, hipparcos=False):
         """
         log likelihood function that interfaces with the orbitize objects
         Comptues the sum of the log likelihoods of the data given the input model
@@ -88,6 +90,19 @@ def _logl(self, params):
 
         if self.custom_lnlike is not None:
             lnlikes_sum += self.custom_lnlike(params)
+        
+
+        if self.system.hipparcos_number is not None:
+
+            # compute Ra/Dec predictions at the Hipparcos IAD epochs
+            raoff_model, deoff_model, _ = self.system.compute_all_orbits(
+                params, epochs=self.system.hipparcos_IAD.epochs_mjd
+            ) 
+
+            # select body 0 raoff/deoff predictions & feed into Hip IAD lnlike fn
+            lnlikes_sum += self.system.hipparcos_IAD.compute_lnlike(
+                raoff_model[:,0,:], deoff_model[:,0,:], params
+            )
 
         return lnlikes_sum
 
@@ -180,7 +195,9 @@ def __init__(self, system, like='chi2_lnlike', custom_lnlike=None):
             lnlike=None,
             tau_ref_epoch=self.system.tau_ref_epoch,
             data=self.system.data_table,
-            num_secondary_bodies=self.system.num_secondary_bodies
+            num_secondary_bodies=self.system.num_secondary_bodies,
+            fitting_basis=self.system.fitting_basis,
+            xyz_epochs= self.system.best_epochs
         )
 
     def prepare_samples(self, num_samples):
@@ -308,6 +325,7 @@ def reject(self, samples):
         # we just want the chi2 term for rejection, so compute the Gaussian normalization term and remove it
         errs = np.array([self.system.data_table['quant1_err'],
                          self.system.data_table['quant2_err']]).T
+
         if self.has_corr:
             corrs = self.system.data_table['quant12_corr']
         else:
@@ -528,7 +546,9 @@ def __init__(self, system, num_temps=20, num_walkers=1000, num_threads=1, like='
             lnlike=None,
             tau_ref_epoch=system.tau_ref_epoch,
             data=self.system.data_table,
-            num_secondary_bodies=system.num_secondary_bodies
+            num_secondary_bodies=system.num_secondary_bodies,
+            fitting_basis=self.system.fitting_basis,
+            xyz_epochs= self.system.best_epochs
         )
         
         if self.num_temps > 1:
@@ -689,6 +709,47 @@ def _update_chains_from_sampler(self, sampler, num_steps=None):
         # include fixed parameters in posterior
         self.post = self._fill_in_fixed_params(self.post)
 
+    def validate_xyz_positions(self):
+        """
+        If using the XYZ basis, walkers might be initialized in an invalid region of parameter space. This function fixes that
+        by replacing invalid positions by new randomly generated positions until all are valid.
+        """
+        if self.system.fitting_basis == 'XYZ':
+            if self.use_pt:
+                all_valid = False
+                while not all_valid:
+                    total_invalids = 0
+                    for temp in range(self.num_temps):
+                        to_stand = orbitize.conversions.xyz_to_standard(self.system.best_epochs[0], self.curr_pos[temp,:,:])
+                        invalids = np.where((to_stand[:, 1] < 0.) | (to_stand[:, 1] >=1.))[0]
+                        if len(invalids) > 0:
+                            newpos = []
+                            for prior in self.priors:
+                                randompos = prior.draw_samples(len(invalids))
+                                newpos.append(randompos)
+                            self.curr_pos[temp, invalids, :] = np.stack(newpos).T 
+                            total_invalids += len(invalids)
+                    if total_invalids == 0:
+                        all_valid = True
+                        print('All walker positions validated.')
+            else:
+                all_valid = False
+                while not all_valid:
+                    total_invalids = 0
+                    to_stand = orbitize.conversions.xyz_to_standard(self.system.best_epochs[0], self.curr_pos[:,:])
+                    invalids = np.where((to_stand[:, 1] < 0.) | (to_stand[:, 1] >=1.))[0]
+                    if len(invalids) > 0:
+                        newpos = []
+                        for prior in self.priors:
+                            randompos = prior.draw_samples(len(invalids))
+                            newpos.append(randompos)
+                        self.curr_pos[invalids, :] = np.stack(newpos).T 
+                        total_invalids += len(invalids)
+                    if total_invalids == 0:
+                        all_valid = True
+                        print('All walker positions validated.')
+
+
     def run_sampler(self, total_orbits, burn_steps=0, thin=1, examine_chains=False, output_filename=None, periodic_save_freq=None):
         """
         Runs PT MCMC sampler. Results are stored in ``self.chain`` and ``self.lnlikes``.
@@ -724,97 +785,98 @@ def run_sampler(self, total_orbits, burn_steps=0, thin=1, examine_chains=False,
         if nsteps <= 0:
             raise ValueError("Total_orbits must be greater than num_walkers.")
 
-        if self.use_pt:
-            sampler = ptemcee.Sampler(
-                self.num_walkers, self.num_params, self._logl, orbitize.priors.all_lnpriors,
-                ntemps=self.num_temps, threads=self.num_threads, logpargs=[self.priors, ]
-            )
-        else:
-            if self.num_threads != 1:
-                print('Setting num_threads=1. If you want parallel processing for emcee implemented in orbitize, let us know.')
-                self.num_threads = 1
-
-            sampler = emcee.EnsembleSampler(
-                self.num_walkers, self.num_params, self._logl,
-                kwargs={'include_logp': True}
-            )
-                
-        print("Starting Burn in")
-        for i, state in enumerate(sampler.sample(self.curr_pos, iterations=burn_steps, thin=thin)):
+        with Pool(processes=self.num_threads) as pool: 
             if self.use_pt:
-                self.curr_pos = state[0]
+                sampler = ptemcee.Sampler(
+                    self.num_walkers, self.num_params, self._logl, orbitize.priors.all_lnpriors,
+                    ntemps=self.num_temps, threads=self.num_threads, logpargs=[self.priors, ]
+                )
             else:
-                self.curr_pos = state.coords
+                # if self.num_threads != 1:
+                #     print('Setting num_threads=1. If you want parallel processing for emcee implemented in orbitize, let us know.')
+                #     self.num_threads = 1
+
+                sampler = emcee.EnsembleSampler(
+                    self.num_walkers, self.num_params, self._logl, pool=pool,
+                    kwargs={'include_logp': True}
+                )
+                    
+            print("Starting Burn in")
+            for i, state in enumerate(sampler.sample(self.curr_pos, iterations=burn_steps, thin=thin)):
+                if self.use_pt:
+                    self.curr_pos = state[0]
+                else:
+                    self.curr_pos = state.coords
 
-            if (i+1) % 5 == 0:
-                print(str(i+1)+'/'+str(burn_steps)+' steps of burn-in complete', end='\r')
+                if (i+1) % 5 == 0:
+                    print(str(i+1)+'/'+str(burn_steps)+' steps of burn-in complete', end='\r')
 
-            if periodic_save_freq is not None:
-                if (i+1) % periodic_save_freq == 0: # we've completed i+1 steps
-                    self.results.curr_pos = self.curr_pos
-                    self.results.save_results(output_filename)
+                if periodic_save_freq is not None:
+                    if (i+1) % periodic_save_freq == 0: # we've completed i+1 steps
+                        self.results.curr_pos = self.curr_pos
+                        self.results.save_results(output_filename)
 
-        sampler.reset()
-        print('')
-        print('Burn in complete. Sampling posterior now.')
-
-        saved_upto = 0 # keep track of how many steps of this chain we've saved. this is the next index that needs to be saved 
-        for i, state in enumerate(sampler.sample(self.curr_pos, iterations=nsteps, thin=thin)):
-            if self.use_pt:
-                self.curr_pos = state[0]
-            else:
-                self.curr_pos = state.coords
-                
-            # print progress statement
-            if (i+1) % 5 == 0:
-                print(str(i+1)+'/'+str(nsteps)+' steps completed', end='\r')
-
-            if periodic_save_freq is not None:
-                if (i+1) % periodic_save_freq == 0: # we've completed i+1 steps
-                    self._update_chains_from_sampler(sampler, num_steps=i+1)
-
-                    # figure out what is the new chunk of the chain and corresponding lnlikes that have been computed before last save
-                    # grab the current posterior and lnlikes and reshape them to have the Nwalkers x Nsteps dimension again
-                    post_shape = self.post.shape
-                    curr_chain_shape = (self.num_walkers, post_shape[0]//self.num_walkers, post_shape[-1])
-                    curr_chain = self.post.reshape(curr_chain_shape)
-                    curr_lnlike_chain = self.lnlikes.reshape(curr_chain_shape[:2])
-                    # use the reshaped arrays and find the new steps we computed
-                    curr_chunk = curr_chain[:, saved_upto:i+1]
-                    curr_chunk = curr_chunk.reshape(-1, curr_chunk.shape[-1]) # flatten nwalkers x nsteps dim
-                    curr_lnlike_chunk = curr_lnlike_chain[:, saved_upto:i+1].flatten()
-
-                    # add this current chunk to the results object (which already has all the previous chunks saved)
-                    self.results.add_samples(curr_chunk, curr_lnlike_chunk, 
-                                                labels=self.system.labels, curr_pos=self.curr_pos)
-                    self.results.save_results(output_filename)
-                    saved_upto = i+1
-
-        print('')
-        self._update_chains_from_sampler(sampler)
-
-        if periodic_save_freq is None:
-            # need to save everything
-            self.results.add_samples(self.post, self.lnlikes, labels=self.system.labels, curr_pos=self.curr_pos)
-        elif saved_upto < nsteps:
-            # just need to save the last few
-            # same code as above except we just need to grab the last few
-            post_shape = self.post.shape
-            curr_chain_shape = (self.num_walkers, post_shape[0]//self.num_walkers, post_shape[-1])
-            curr_chain = self.post.reshape(curr_chain_shape)
-            curr_lnlike_chain = self.lnlikes.reshape(curr_chain_shape[:2])
-            curr_chunk = curr_chain[:, saved_upto:]
-            curr_chunk = curr_chunk.reshape(-1, curr_chunk.shape[-1]) # flatten nwalkers x nsteps dim
-            curr_lnlike_chunk = curr_lnlike_chain[:, saved_upto:].flatten()
-
-            self.results.add_samples(curr_chunk, curr_lnlike_chunk, 
-                                                labels=self.system.labels, curr_pos=self.curr_pos)
-
-        if output_filename is not None:
-            self.results.save_results(output_filename)
-
-        print('Run complete')
+            sampler.reset()
+            print('')
+            print('Burn in complete. Sampling posterior now.')
 
+            saved_upto = 0 # keep track of how many steps of this chain we've saved. this is the next index that needs to be saved 
+            for i, state in enumerate(sampler.sample(self.curr_pos, iterations=nsteps, thin=thin)):
+                if self.use_pt:
+                    self.curr_pos = state[0]
+                else:
+                    self.curr_pos = state.coords
+                    
+                # print progress statement
+                if (i+1) % 5 == 0:
+                    print(str(i+1)+'/'+str(nsteps)+' steps completed', end='\r')
+
+                if periodic_save_freq is not None:
+                    if (i+1) % periodic_save_freq == 0: # we've completed i+1 steps
+                        self._update_chains_from_sampler(sampler, num_steps=i+1)
+
+                        # figure out what is the new chunk of the chain and corresponding lnlikes that have been computed before last save
+                        # grab the current posterior and lnlikes and reshape them to have the Nwalkers x Nsteps dimension again
+                        post_shape = self.post.shape
+                        curr_chain_shape = (self.num_walkers, post_shape[0]//self.num_walkers, post_shape[-1])
+                        curr_chain = self.post.reshape(curr_chain_shape)
+                        curr_lnlike_chain = self.lnlikes.reshape(curr_chain_shape[:2])
+                        # use the reshaped arrays and find the new steps we computed
+                        curr_chunk = curr_chain[:, saved_upto:i+1]
+                        curr_chunk = curr_chunk.reshape(-1, curr_chunk.shape[-1]) # flatten nwalkers x nsteps dim
+                        curr_lnlike_chunk = curr_lnlike_chain[:, saved_upto:i+1].flatten()
+
+                        # add this current chunk to the results object (which already has all the previous chunks saved)
+                        self.results.add_samples(curr_chunk, curr_lnlike_chunk, 
+                                                    labels=self.system.labels, curr_pos=self.curr_pos)
+                        self.results.save_results(output_filename)
+                        saved_upto = i+1
+
+            print('')
+            self._update_chains_from_sampler(sampler)
+
+            if periodic_save_freq is None:
+                # need to save everything
+                self.results.add_samples(self.post, self.lnlikes, labels=self.system.labels, curr_pos=self.curr_pos)
+            elif saved_upto < nsteps:
+                # just need to save the last few
+                # same code as above except we just need to grab the last few
+                post_shape = self.post.shape
+                curr_chain_shape = (self.num_walkers, post_shape[0]//self.num_walkers, post_shape[-1])
+                curr_chain = self.post.reshape(curr_chain_shape)
+                curr_lnlike_chain = self.lnlikes.reshape(curr_chain_shape[:2])
+                curr_chunk = curr_chain[:, saved_upto:]
+                curr_chunk = curr_chunk.reshape(-1, curr_chunk.shape[-1]) # flatten nwalkers x nsteps dim
+                curr_lnlike_chunk = curr_lnlike_chain[:, saved_upto:].flatten()
+
+                self.results.add_samples(curr_chunk, curr_lnlike_chunk, 
+                                                    labels=self.system.labels, curr_pos=self.curr_pos)
+
+            if output_filename is not None:
+                self.results.save_results(output_filename)
+
+            print('Run complete')
+        # Close pool
         if examine_chains:
             self.examine_chains()
 
@@ -936,7 +998,9 @@ def chop_chains(self, burn, trim=0):
             lnlike=flat_chopped_lnlikes,
             tau_ref_epoch=self.system.tau_ref_epoch,
             labels=self.system.labels,
-            num_secondary_bodies=self.system.num_secondary_bodies
+            num_secondary_bodies=self.system.num_secondary_bodies,
+            fitting_basis=self.system.fitting_basis,
+            xyz_epochs= self.system.best_epochs
         )
 
         # Print a confirmation
diff --git a/orbitize/system.py b/orbitize/system.py
index feab1196..d17aaf81 100644
--- a/orbitize/system.py
+++ b/orbitize/system.py
@@ -1,6 +1,7 @@
 import numpy as np
-from orbitize import priors, read_input, kepler
-
+from orbitize import priors, read_input, kepler, conversions, hipparcos
+import astropy.units as u
+import astropy.constants as consts
 
 class System(object):
     """
@@ -28,7 +29,6 @@ class System(object):
         results (list of orbitize.results.Results): results from an orbit-fit
             will be appended to this list as a Results class.
 
-
     Users should initialize an instance of this class, then overwrite
     priors they wish to customize.
 
@@ -52,7 +52,8 @@ class System(object):
 
     def __init__(self, num_secondary_bodies, data_table, stellar_mass,
                  plx, mass_err=0, plx_err=0, restrict_angle_ranges=None,
-                 tau_ref_epoch=58849, fit_secondary_mass=False, results=None):
+                 tau_ref_epoch=58849, fit_secondary_mass=False, results=None,
+                 hipparcos_number=None, fitting_basis='standard', hipparcos_filename=None):
 
         self.num_secondary_bodies = num_secondary_bodies
         self.sys_priors = []
@@ -61,6 +62,8 @@ def __init__(self, num_secondary_bodies, data_table, stellar_mass,
         self.fit_secondary_mass = fit_secondary_mass
         self.tau_ref_epoch = tau_ref_epoch
         self.restrict_angle_ranges = restrict_angle_ranges
+        self.hipparcos_number = hipparcos_number
+        self.fitting_basis = fitting_basis
 
         #
         # Group the data in some useful ways
@@ -144,6 +147,9 @@ def __init__(self, num_secondary_bodies, data_table, stellar_mass,
                                      ((len(self.radec[1]) + len(self.seppa[1]) + len(self.rv[1]) < len(data_table)) or \
                                       (self.num_secondary_bodies > 1))
 
+        if self.hipparcos_number is not None:
+            self.track_planet_perturbs = True
+
         if restrict_angle_ranges:
             angle_upperlim = np.pi
         else:
@@ -153,30 +159,154 @@ def __init__(self, num_secondary_bodies, data_table, stellar_mass,
         # Set priors for each orbital element
         #
 
-        for body in np.arange(num_secondary_bodies):
-            # Add semimajor axis prior
-            self.sys_priors.append(priors.LogUniformPrior(0.001, 1e7))
-            self.labels.append('sma{}'.format(body+1))
-
-            # Add eccentricity prior
-            self.sys_priors.append(priors.UniformPrior(0., 1.))
-            self.labels.append('ecc{}'.format(body+1))
-
-            # Add inclination angle prior
-            self.sys_priors.append(priors.SinPrior())
-            self.labels.append('inc{}'.format(body+1))
-
-            # Add argument of periastron prior
-            self.sys_priors.append(priors.UniformPrior(0., 2.*np.pi))
-            self.labels.append('aop{}'.format(body+1))
-
-            # Add position angle of nodes prior
-            self.sys_priors.append(priors.UniformPrior(0., angle_upperlim))
-            self.labels.append('pan{}'.format(body+1))
-
-            # Add epoch of periastron prior.
-            self.sys_priors.append(priors.UniformPrior(0., 1.))
-            self.labels.append('tau{}'.format(body+1))
+        if fitting_basis == 'standard':
+            self.best_epochs = None
+            for body in np.arange(num_secondary_bodies):
+                # Add semimajor axis prior
+                self.sys_priors.append(priors.LogUniformPrior(0.001, 1e7))
+                self.labels.append('sma{}'.format(body+1))
+
+                # Add eccentricity prior
+                self.sys_priors.append(priors.UniformPrior(0., 1.))
+                self.labels.append('ecc{}'.format(body+1))
+
+                # Add inclination angle prior
+                self.sys_priors.append(priors.SinPrior())
+                # self.sys_priors.append(priors.UniformPrior(0., np.pi))# TEST TO COMPARE, CHANGE LATER
+                self.labels.append('inc{}'.format(body+1))
+
+                # Add argument of periastron prior
+                self.sys_priors.append(priors.UniformPrior(0., 2.*np.pi))
+                self.labels.append('aop{}'.format(body+1))
+
+                # Add position angle of nodes prior
+                self.sys_priors.append(priors.UniformPrior(0., angle_upperlim))
+                self.labels.append('pan{}'.format(body+1))
+
+                # Add epoch of periastron prior.
+                self.sys_priors.append(priors.UniformPrior(0., 1.))
+                self.labels.append('tau{}'.format(body+1))
+
+        elif fitting_basis == 'XYZ':
+            epochs = self.data_table['epoch']
+            # Get epochs with least uncertainty, as is done in sampler.py
+            convert_warning_print = False
+            for body_num in np.arange(self.num_secondary_bodies) + 1:
+                if len(self.radec[body_num]) > 0:
+                    # only print the warning once. 
+                    if not convert_warning_print:
+                        print('Converting ra/dec data points in data_table to sep/pa. Original data are stored in input_table.')
+                        convert_warning_print = True
+                    self.convert_data_table_radec2seppa(body_num=body_num)
+
+            sep_err = self.data_table[np.where(
+                self.data_table['quant_type'] == 'seppa')]['quant1_err'].copy()
+            meas_object = self.data_table[np.where(
+                self.data_table['quant_type'] == 'seppa')]['object'].copy()
+
+            self.best_epochs = []
+            self.best_epoch_idx = []
+            min_sep_indices = np.argsort(sep_err) # indices of sep err sorted from smallest to higheset
+            min_sep_indices_body = meas_object[min_sep_indices] # the corresponding body_num that these sorted measurements correspond to
+            for i in range(self.num_secondary_bodies):
+                body_num = i + 1
+                this_object_meas = np.where(min_sep_indices_body == body_num)[0]
+                if np.size(this_object_meas) == 0:
+                    # no data, no scaling
+                    self.best_epochs.append(None)
+                    continue
+                # get the smallest measurement belonging to this body
+                this_best_epoch_idx = min_sep_indices[this_object_meas][0] # already sorted by argsort
+                self.best_epoch_idx.append(this_best_epoch_idx)
+                this_best_epoch = epochs[this_best_epoch_idx]
+                self.best_epochs.append(this_best_epoch)
+
+            for body in np.arange(num_secondary_bodies):
+                # Get the epoch with the least uncertainty for this body
+                # curr_idx = self.body_indices[body_num]
+                # radec_uncerts = self.data_table['quant1_err'][curr_idx] + self.data_table['quant2_err'][curr_idx]
+                # min_uncert = np.where(radec_uncerts == np.amin(radec_uncerts))[0]
+                # best_idx = curr_idx[0][min_uncert[0]]
+                datapoints_to_take = 3
+                best_idx = self.best_epoch_idx[body]
+                best_epochs = epochs[best_idx:(best_idx+datapoints_to_take)] # 0 is best, the others are for fitting velocity
+
+                # Get data near best epoch ASSUMING THE BEST IS NOT ONE OF THE LAST TWO EPOCHS OF A GIVEN BODY,
+                # also assuming this is in radec
+                best_ras = self.input_table['quant1'][best_idx:(best_idx+datapoints_to_take)].copy()
+                best_ras_err = self.input_table['quant1_err'][best_idx:(best_idx+datapoints_to_take)].copy()
+                best_decs =self.input_table['quant2'][best_idx:(best_idx+datapoints_to_take)].copy()
+                best_decs_err = self.input_table['quant2_err'][best_idx:(best_idx+datapoints_to_take)].copy()
+
+                # Convert to AU for prior limits
+                best_xs = best_ras / plx 
+                best_ys = best_decs / plx 
+                best_xs_err = np.sqrt((best_ras_err / best_ras)**2 + (plx_err / plx)**2)*np.absolute(best_xs)
+                best_ys_err = np.sqrt((best_decs_err / best_decs)**2 + (plx_err / plx)**2)*np.absolute(best_ys)
+
+                # Least-squares fit on velocity for prior limits
+                A = np.vander(best_epochs, 2)
+
+                ATA_x = np.dot(A.T, A / (best_xs_err ** 2)[:, None])
+                cov_x = np.linalg.inv(ATA_x)
+                w_x = np.linalg.solve(ATA_x, np.dot(A.T, best_xs / best_xs_err ** 2))
+
+                ATA_y = np.dot(A.T, A / (best_ys_err ** 2)[:, None])
+                cov_y = np.linalg.inv(ATA_y)
+                w_y = np.linalg.solve(ATA_y, np.dot(A.T, best_ys / best_ys_err ** 2))
+
+                x_vel = w_x[0]
+                x_vel_err = np.sqrt(cov_x[0, 0])
+                y_vel = w_y[0]
+                y_vel_err = np.sqrt(cov_y[0, 0])
+
+                x_vel = (( x_vel* u.AU / u.day).to(u.km / u.s)).value
+                x_vel_err = ((x_vel_err * u.AU / u.day).to(u.km / u.s)).value
+                y_vel = ((y_vel * u.AU / u.day).to(u.km / u.s)).value
+                y_vel_err = ((y_vel_err * u.AU / u.day).to(u.km / u.s)).value
+
+                # Propose bounds on absolute Z and Z dot given the energy equation
+                mu = consts.G * stellar_mass * u.Msun
+
+                mu_vel = 2 * mu / ((x_vel**2 + y_vel**2) * (u.km / u.s * u.km / u.s))
+                z_bound = (np.sqrt(mu_vel**2 - (best_xs[0]**2 + best_ys[0]**2)*u.AU *u.AU)).to(u.AU)
+                z_bound = z_bound.value
+
+                mu_pos = 2 * mu / np.sqrt((best_xs[0]**2 + best_ys[0]**2) * (u.AU *u.AU))
+                z_vel_bound = (np.sqrt(mu_pos - (x_vel**2 + y_vel**2)*(u.km / u.s * u.km / u.s))).to(u.km / u.s)
+                z_vel_bound = z_vel_bound.value
+
+                # Add x-coordinate prior
+                num_uncerts_x = 5
+                self.sys_priors.append(priors.UniformPrior(best_xs[0] - num_uncerts_x*best_xs_err[0], best_xs[0] + num_uncerts_x*best_xs_err[0]))
+                self.labels.append('x{}'.format(body+1))
+                
+                # Add y-coordinate prior
+                num_uncerts_y = 5
+                self.sys_priors.append(priors.UniformPrior(best_ys[0] - num_uncerts_y*best_ys_err[0], best_ys[0] + num_uncerts_y*best_ys_err[0]))
+                self.labels.append('y{}'.format(body+1))
+
+                # Add z-coordinate prior
+                # self.sys_priors.append(priors.UniformPrior(-z_bound,z_bound))
+                # self.sys_priors.append(priors.LogUniformPrior(0.0001,z_bound))
+                self.sys_priors.append(priors.GaussianPrior(0.,z_bound / 4, no_negatives=False))
+                self.labels.append('z{}'.format(body+1))
+
+                # Add x-velocity prior
+                num_uncerts_xvel = 5
+                self.sys_priors.append(priors.UniformPrior(x_vel - num_uncerts_xvel*x_vel_err, x_vel + num_uncerts_xvel*x_vel_err))
+                self.labels.append('xdot{}'.format(body+1))
+
+                # Add y-velocity prior
+                num_uncerts_yvel = 5
+                self.sys_priors.append(priors.UniformPrior(y_vel - num_uncerts_yvel*y_vel_err, y_vel + num_uncerts_yvel*y_vel_err))
+                self.labels.append('ydot{}'.format(body+1))
+
+                # Add z-velocity prior
+                # self.sys_priors.append(priors.UniformPrior(-z_vel_bound,z_vel_bound))
+                # self.sys_priors.append(priors.LogUniformPrior(0.0001,z_vel_bound))
+                self.sys_priors.append(priors.GaussianPrior(0.,z_vel_bound / 4, no_negatives=False))
+                self.labels.append('zdot{}'.format(body+1))
 
         #
         # Set priors on total mass and parallax
@@ -187,10 +317,40 @@ def __init__(self, num_secondary_bodies, data_table, stellar_mass,
         else:
             self.sys_priors.append(plx)
 
-        # checking for rv data to include appropriate rv priors:
+        # instantiate a HipparcosLogProb object to precompute & hold params relevant to IAD
+        if hipparcos_number is not None:
+            self.hipparcos_IAD = hipparcos.HipparcosLogProb(
+                hipparcos_filename, hipparcos_number, self.num_secondary_bodies
+            )
+
+            # for now, set broad uniform priors on astrometric params relevant for Hipparcos
+            self.sys_priors.append(priors.UniformPrior(
+                self.hipparcos_IAD.pm_ra0 - 10 * self.hipparcos_IAD.pm_ra0_err,
+                self.hipparcos_IAD.pm_ra0 + 10 * self.hipparcos_IAD.pm_ra0_err)
+            )
+            self.labels.append('pm_ra')
+
+            self.sys_priors.append(priors.UniformPrior(
+                self.hipparcos_IAD.pm_dec0 - 10 * self.hipparcos_IAD.pm_dec0_err,
+                self.hipparcos_IAD.pm_dec0 + 10 * self.hipparcos_IAD.pm_dec0_err)
+            )
+            self.labels.append('pm_dec')
 
+            self.sys_priors.append(priors.UniformPrior(
+                - 10 * self.hipparcos_IAD.alpha0_err,
+                10 * self.hipparcos_IAD.alpha0_err)
+            )
+            self.labels.append('alpha0')
+
+            self.sys_priors.append(priors.UniformPrior(
+                - 10 * self.hipparcos_IAD.delta0_err,
+                10 * self.hipparcos_IAD.delta0_err)
+            )
+            self.labels.append('delta0')
+
+        # checking for rv data to include appropriate rv priors:
         if len(self.rv[0]) > 0 and self.fit_secondary_mass:
-            # Rob and Lea:
+         
             # for instrument in rv_instruments:
                 # add gamma and sigma for each and label each unique gamma and sigma per instrument name (gamma+instrument1, ...)
             for instrument in self.rv_instruments:
@@ -217,9 +377,10 @@ def __init__(self, num_secondary_bodies, data_table, stellar_mass,
         # add labels dictionary for parameter indexing
         self.param_idx = dict(zip(self.labels, np.arange(len(self.labels))))
 
-    def compute_model(self, params_arr):
+    def compute_all_orbits(self, params_arr, epochs=None):
         """
-        Compute model predictions for an array of fitting parameters.
+        Calls orbitize.kepler.calc_orbit and optionally accounts for multi-body
+        interactions, as well as computes total quantities like RV (without jitter/gamma)
 
         Args:
             params_arr (np.array of float): RxM array
@@ -227,55 +388,81 @@ def compute_model(self, params_arr):
                 parameters being fit, and M is the number of orbits
                 we need model predictions for. Must be in the same order
                 documented in ``System()`` above. If M=1, this can be a 1d array.
-
+        
         Returns:
-            np.array of float: Nobsx2xM array model predictions. If M=1, this is
-            a 2d array, otherwise it is a 3d array.
+            tuple of:
+                raoff (np.array of float): N_epochs x N_bodies x N_orbits array of
+                    RA offsets from barycenter at each epoch.
+                decoff (np.array of float): N_epochs x N_bodies x N_orbits array of
+                    Dec offsets from barycenter at each epoch.
+                vz (np.array of float): N_epochs x N_bodies x N_orbits array of
+                    radial velocities at each epoch.
+
         """
 
+        if epochs is None:
+            epochs = self.data_table['epoch']
+
+        n_epochs = len(epochs)
+
         if len(params_arr.shape) == 1:
-            model = np.zeros((len(self.data_table), 2))
-            # Rob and Lea: adding gamma zeros
-            gamma = np.zeros((len(self.data_table), 2))
-            jitter = np.zeros((len(self.data_table), 2))
+            n_orbits = 1
         else:
-            model = np.zeros((len(self.data_table), 2, params_arr.shape[1]))
-            gamma = np.zeros((len(self.data_table), 2, params_arr.shape[1]))
-            jitter = np.zeros((len(self.data_table), 2, params_arr.shape[1]))
-            
-        if self.track_planet_perturbs:
-            radec_perturb = np.zeros(model.shape)
+            n_orbits = params_arr.shape[1]
 
-        if len(self.rv[0]) > 0 and self.fit_secondary_mass:
-            # looping through instruments to get the gammas
-            for rv_idx in range(len(self.rv_instruments)):
-                gamma[self.rv_inst_indices[rv_idx], 0] = params_arr[6 *
-                                                                    self.num_secondary_bodies+1+2*rv_idx]  # km/s
-                jitter[self.rv_inst_indices[rv_idx], 0] = params_arr[6 *
-                                                                     self.num_secondary_bodies+2+2*rv_idx]
-                gamma[self.rv_inst_indices[rv_idx], 1] = np.nan
-                jitter[self.rv_inst_indices[rv_idx], 1] = np.nan
-            # need to put planetary rv later
+        ra_kepler = np.zeros((n_epochs, self.num_secondary_bodies + 1, n_orbits)) # N_epochs x N_bodies x N_orbits
+        dec_kepler = np.zeros((n_epochs, self.num_secondary_bodies + 1, n_orbits))
+
+        ra_perturb = np.zeros((n_epochs, self.num_secondary_bodies + 1, n_orbits)) 
+        dec_perturb = np.zeros((n_epochs, self.num_secondary_bodies + 1, n_orbits))
+
+        vz = np.zeros((n_epochs, self.num_secondary_bodies + 1, n_orbits))
+
+        # mass/mtot used to compute each Keplerian orbit will be needed later to compute perturbations
+        if self.track_planet_perturbs:
+            masses = np.zeros((self.num_secondary_bodies + 1, n_orbits))
+            mtots = np.zeros((self.num_secondary_bodies + 1, n_orbits))
 
-            total_rv0 = 0  # If we're not fitting rv, then we don't regard the total rv and will not use this
+        total_rv0 = 0
 
         for body_num in np.arange(self.num_secondary_bodies)+1:
-            # we're going to compute at all epochs for convenience of indexing right now
-            # self.radec, and self.seppa index into the entire data table, not just the values for a particular body
-            epochs = self.data_table['epoch']#[self.body_indices[body_num]] 
+
             startindex = 6 * (body_num - 1)
-            sma = params_arr[startindex]
-            ecc = params_arr[startindex + 1]
-            inc = params_arr[startindex + 2]
-            argp = params_arr[startindex + 3]
-            lan = params_arr[startindex + 4]
-            tau = params_arr[startindex + 5]
+            if self.fitting_basis == 'standard':
+
+                sma = params_arr[startindex]
+                ecc = params_arr[startindex + 1]
+                inc = params_arr[startindex + 2]
+                argp = params_arr[startindex + 3]
+                lan = params_arr[startindex + 4]
+                tau = params_arr[startindex + 5]
+            
+            elif self.fitting_basis == 'XYZ':
+                
+                # curr_idx = self.body_indices[body_num]
+                # radec_uncerts = self.data_table['quant1_err'][curr_idx] + self.data_table['quant2_err'][curr_idx]
+                # min_uncert = np.where(radec_uncerts == np.amin(radec_uncerts))
+                best_idx = self.best_epoch_idx[body_num-1]
+                constrained_epoch = epochs[best_idx]
+
+                to_convert = np.array([*params_arr[startindex:(startindex+6)],params_arr[6 * self.num_secondary_bodies],params_arr[-1]])
+                standard_params = conversions.xyz_to_standard(constrained_epoch, to_convert)
+
+                sma = standard_params[0] 
+                ecc = standard_params[1]
+                inc = standard_params[2]
+                argp = standard_params[3]
+                lan = standard_params[4]
+                tau = standard_params[5]
+
             plx = params_arr[6 * self.num_secondary_bodies]
 
             if self.fit_secondary_mass:
+
                 # mass of secondary bodies are in order from -1-num_bodies until -2 in order.
                 mass = params_arr[-1-self.num_secondary_bodies+(body_num-1)]
                 m0 = params_arr[-1]
+
                 # For what mtot to use to calculate central potential, we should use the mass enclosed in a sphere with r <= distance of planet. 
                 # We need to select all planets with sma < this planet. 
                 all_smas = params_arr[0:6*self.num_secondary_bodies:6]
@@ -284,105 +471,174 @@ def compute_model(self, params_arr):
                 all_pl_masses = params_arr[-1-self.num_secondary_bodies:-1]
                 inside_masses = all_pl_masses[within_orbit]
                 mtot = np.sum(inside_masses) + m0
+
             else:
                 # if not fitting for secondary mass, then total mass must be stellar mass
                 mass = None
                 m0 = None
                 mtot = params_arr[-1]
+            
+            if self.track_planet_perturbs:
+                masses[body_num] = mass
+                mtots[body_num] = mtot
 
-            # i = 1,2,3... (companion index)
-
+            # solve Kepler's equation
             raoff, decoff, vz_i = kepler.calc_orbit(
                 epochs, sma, ecc, inc, argp, lan, tau, plx, mtot,
                 mass_for_Kamp=m0, tau_ref_epoch=self.tau_ref_epoch, tau_warning=False
             )
 
-            # raoff, decoff, vz are scalers if the length of epochs is 1. 
-            # Jason is too lazy to figure out how to make it return a one element array without breaking everything else
-            # so hard code it here to convert them into 1-element numpy arrays. 
+            # raoff, decoff, vz are scalers if the length of epochs is 1
             if len(epochs) == 1:
                 raoff = np.array([raoff])
                 decoff = np.array([decoff])
-                vz = np.array([vz_i])
+                vz_i = np.array([vz_i])
+
+            if n_orbits == 1:
+                raoff = raoff.reshape((n_epochs, 1))
+                decoff = decoff.reshape((n_epochs, 1))
+                vz_i = vz_i.reshape((n_epochs, 1))
 
+            # add Keplerian ra/deoff for this body to storage arrays
+            ra_kepler[:, body_num, :] = raoff 
+            dec_kepler[:, body_num, :] = decoff
+            vz[:, body_num, :] = vz_i
 
             # vz_i is the ith companion radial velocity
-            if len(self.rv[0]) > 0 and self.fit_secondary_mass:
-                vz0 = vz_i*-(mass/m0)  # calculating stellar velocity due to ith companion
-                total_rv0 = total_rv0 + vz0  # Adding stellar velocity and gamma
+            if self.fit_secondary_mass:
+                vz0 = vz_i * -(mass / m0)  # calculating stellar velocity due to ith companion
+                total_rv0 = total_rv0 + vz0  # adding stellar velocity and gamma
+
+        # if we are fitting for the mass of the planets, then they will perturb the star
+        # add the perturbation on the star due to this planet on the relative astrometry of the planet that was measured
+        # We are superimposing the Keplerian orbits, so we can add it linearly, scaled by the mass. 
+        # Because we are in Jacobi coordinates, for companions, we only should model the effect of planets interior to it. 
+        # (Jacobi coordinates mean that separation for a given companion is measured relative to the barycenter of all interior companions)
+        if self.track_planet_perturbs:
+            for body_num in np.arange(self.num_secondary_bodies + 1):
+                if body_num > 0:
+                    # for companions, only perturb companion orbits at larger SMAs than this one. 
+                    startindex = 6 * (body_num - 1) # subtract 1 because object 1 is 0th companion
+                    sma = params_arr[startindex]
+                    all_smas = params_arr[0:6*self.num_secondary_bodies:6]
+                    outside_orbit = np.where(all_smas > sma)[0]
+                    which_perturb_bodies = outside_orbit + 1
+
+                    # the planet will also perturb the star
+                    which_perturb_bodies = np.append([0], which_perturb_bodies)
+
+                else:
+                    # for the star, what we are measuring is its position relative to the system barycenter
+                    # so we want to account for all of the bodies.  
+                    which_perturb_bodies = np.arange(self.num_secondary_bodies+1)
+
+                for other_body_num in which_perturb_bodies:
+                    # skip itself since the the 2-body problem is measuring the planet-star separation already
+                    if (body_num == other_body_num) | (body_num == 0):
+                        continue
+
+                    ## NOTE: we are only handling astrometry right now (TODO: integrate RV into this)
+                    ra_perturb[:, other_body_num, :] += (masses[other_body_num]/mtots[other_body_num]) * ra_kepler[:, body_num, :]
+                    dec_perturb[:, other_body_num, :] += (masses[body_num]/mtots[body_num]) * dec_kepler[:, body_num, :] 
+                    
+                    # star is perturbed in opposite direction
+                    if other_body_num == 0:
+                        ra_perturb[:, other_body_num, :] *= -1
+                        dec_perturb[:, other_body_num, :] *= -1
+
+        raoff = ra_kepler + ra_perturb
+        deoff = dec_kepler + dec_perturb
+        vz[:, 0, :] = total_rv0
+        if self.fitting_basis == 'XYZ':
+            # To filter out unbound orbits
+            if ((ecc >= 1.) | (ecc < 0.)):
+                raoff[:,:,:] = np.inf
+                deoff[:,:,:] = np.inf 
+                vz[:,:,:] = np.inf
+                return raoff, deoff, vz
+            else: 
+                return raoff, deoff, vz 
+        elif self.fitting_basis == 'standard':
+            return raoff, deoff, vz
+
+
+    def compute_model(self, params_arr):
+        """
+        Compute model predictions for an array of fitting parameters. 
+        Calls the above compute_all_orbits() function, adds jitter/gamma to
+        RV measurements, and propagates these predictions to a model array that
+        can be subtracted from a data array to compute chi2. 
+        
+        Args:
+            params_arr (np.array of float): RxM array
+                of fitting parameters, where R is the number of
+                parameters being fit, and M is the number of orbits
+                we need model predictions for. Must be in the same order
+                documented in ``System()`` above. If M=1, this can be a 1d array.
+
+        Returns:
+            np.array of float: Nobsx2xM array model predictions. If M=1, this is
+            a 2d array, otherwise it is a 3d array.
+        """
+
+        raoff, decoff, vz = self.compute_all_orbits(params_arr)
+
+        if len(params_arr.shape) == 1:
+            n_orbits = 1
+        else:
+            n_orbits = params_arr.shape[1]
+
+        n_epochs = len(self.data_table)
+        model = np.zeros((n_epochs, 2, n_orbits))
+        jitter = np.zeros((n_epochs, 2, n_orbits))
+        gamma = np.zeros((n_epochs, 2, n_orbits))
+
+        if len(self.rv[0]) > 0 and self.fit_secondary_mass: 
+
+            # looping through instruments to get the gammas & jitters
+            for rv_idx in range(len(self.rv_instruments)):
+
+                jitter[self.rv_inst_indices[rv_idx], 0] = params_arr[ # [km/s]
+                    6 * self.num_secondary_bodies+2+2*rv_idx
+                ]
+                jitter[self.rv_inst_indices[rv_idx], 1] = np.nan
+
+
+                gamma[self.rv_inst_indices[rv_idx], 0] = params_arr[
+                    6 * self.num_secondary_bodies+1+2*rv_idx
+                ] 
+                gamma[self.rv_inst_indices[rv_idx], 1] = np.nan
+
+        for body_num in np.arange(self.num_secondary_bodies + 1):
 
             # for the model points that correspond to this planet's orbit, add the model prediction
             # RA/Dec
             if len(self.radec[body_num]) > 0: # (prevent empty array dimension errors)
-                model[self.radec[body_num], 0] = raoff[self.radec[body_num]]
-                model[self.radec[body_num], 1] = decoff[self.radec[body_num]]
+                model[self.radec[body_num], 0] = raoff[self.radec[body_num], body_num, :]  # N_epochs x N_bodies x N_orbits
+                model[self.radec[body_num], 1] = decoff[self.radec[body_num], body_num, :]
 
             # Sep/PA
             if len(self.seppa[body_num]) > 0:
                 sep, pa = radec2seppa(raoff, decoff)
 
-                model[self.seppa[body_num], 0] = sep[self.seppa[body_num]]
-                model[self.seppa[body_num], 1] = pa[self.seppa[body_num]]
+                model[self.seppa[body_num], 0] = sep[self.seppa[body_num], body_num, :]
+                model[self.seppa[body_num], 1] = pa[self.seppa[body_num], body_num, :]
 
             # RV
             if len(self.rv[body_num]) > 0:
-                model[self.rv[body_num], 0] = vz_i[self.rv[body_num]]
+                model[self.rv[body_num], 0] = vz[self.rv[body_num], body_num, :]
                 model[self.rv[body_num], 1] = np.nan
 
-            # for the other epochs, if we are fitting for the mass of the planets, then they will perturb the star
-            # add the perturbation on the star due to this planet on the relative astrometry of the planet that was measured
-            # We are superimposing the Keplerian orbits, so we can add it linearly, scaled by the mass. 
-            # Because we are in Jacobi coordinates, for companions, we only should model the effect of planets interior to it. 
-            # (Jacobi coordinates mean that separation for a given companion is measured relative to the barycenter of all interior companions)
-            if self.track_planet_perturbs:
-                if body_num > 0:
-                    # for companions, only perturb companion orbits at larger SMAs than this one. 
-                    # note the +1, since the 0th planet is body_num 1. 
-                    which_perturb_bodies = outside_orbit[0] + 1
-                else:
-                    # for the star, what we are measuring is it's position relative to the system barycenter
-                    # so we want to account for all of the bodies.  
-                    which_perturb_bodies = range(self.num_secondary_bodies+1)
-                    
-                for other_body_num in which_perturb_bodies:
-                    # skip itself since the the 2-body problem is measuring the planet-star separation already
-                    if (body_num == other_body_num) | (body_num == 0):
-                        continue
-                    ## NOTE: we are only handling ra/dec and sep/pa right now
-                    ## TODO: integrate RV into this
-                    if len(self.radec[other_body_num]) > 0:
-                        radec_perturb[self.radec[other_body_num], 0] += -(mass/mtot) * raoff[self.radec[other_body_num]]
-                        radec_perturb[self.radec[other_body_num], 1] += -(mass/mtot) * decoff[self.radec[other_body_num]] 
-                    if len(self.seppa[other_body_num]) > 0:
-                        radec_perturb[self.seppa[other_body_num], 0] += -(mass/mtot) * raoff[self.seppa[other_body_num]]
-                        radec_perturb[self.seppa[other_body_num], 1] += -(mass/mtot) * decoff[self.seppa[other_body_num]]
-
-        if len(self.rv[0]) > 0 and self.fit_secondary_mass:
-            if len(total_rv0[self.rv[0]]) > 0:
-                model[self.rv[0], 0] = total_rv0[self.rv[0]]
-                model[self.rv[0], 1] = np.nan  # nans only for rv indices
-
-        # add the effects of other planets on the measured astrometry
-        if self.track_planet_perturbs:
-            for body_num in range(self.num_secondary_bodies+1):
-                if len(self.radec[body_num]) > 0:
-                    model[self.radec[body_num]] -= radec_perturb[self.radec[body_num]]
-
-                if len(self.seppa[body_num]) > 0:
-                    # for seppa, add the perturbations in radec space and convert back
-                    ra_unperturb, dec_unperturb = seppa2radec(model[self.seppa[body_num], 0], model[self.seppa[body_num], 1])
-                    ra_perturb = ra_unperturb - radec_perturb[self.seppa[body_num], 0]
-                    dec_perturb = dec_unperturb - radec_perturb[self.seppa[body_num], 1]
-                    sep_perturb, pa_perturb = radec2seppa(ra_perturb, dec_perturb)
-
-                    model[self.seppa[body_num], 0] = sep_perturb
-                    model[self.seppa[body_num], 1] = pa_perturb
-
+        if n_orbits == 1:
+            model = model.reshape((n_epochs, 2))
+            jitter = jitter.reshape((n_epochs, 2))
+            gamma = gamma.reshape((n_epochs, 2))
 
         if self.fit_secondary_mass:
-            return model+gamma, jitter
+            return model + gamma, jitter
         else:
             return model, jitter
+
     def convert_data_table_radec2seppa(self, body_num=1):
         """
         Converts rows of self.data_table given in radec to seppa.
@@ -484,8 +740,8 @@ def seppa2radec(sep, pa):
 def transform_errors(x1, x2, x1_err, x2_err, x12_corr, transform_func, nsamps=100000):
     """
     Transform errors and covariances from one basis to another using a Monte Carlo apporach
-
-    Args:
+    
+   Args:
         x1 (float): planet location in first coordinate (e.g., RA, sep) before transformation
         x2 (float): planet location in the second coordinate (e.g., Dec, PA) before transformation)
         x1_err (float): error in x1
@@ -494,7 +750,6 @@ def transform_errors(x1, x2, x1_err, x2_err, x12_corr, transform_func, nsamps=10
         transform_func (function): function that transforms between (x1, x2) and (x1p, x2p) (the transformed coordinates)
                                     The function signature should look like: `x1p, x2p = transform_func(x1, x2)`
         nsamps (int): number of samples to draw more the Monte Carlo approach. More is slower but more accurate. 
-
     Returns:
         tuple (x1p_err, x2p_err, x12p_corr): the errors and correlations for x1p,x2p (the transformed coordinates)
     """
diff --git a/requirements.txt b/requirements.txt
index 5674bc85..3a7282e1 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -12,4 +12,5 @@ pytest
 pandas
 radvel
 pyerfa
+astroquery
 sphinx==3.5.3
\ No newline at end of file


Test patch:
diff --git a/tests/end-to-end-tests/betaPic_hipIAD.py b/tests/end-to-end-tests/betaPic_hipIAD.py
new file mode 100644
index 00000000..699de779
--- /dev/null
+++ b/tests/end-to-end-tests/betaPic_hipIAD.py
@@ -0,0 +1,46 @@
+import os
+import matplotlib.pyplot as plt
+
+import orbitize
+from orbitize import system, read_input, priors, sampler
+
+"""
+Attempts to reproduce case 3 (see table 3) of Nielsen+ 2020 (orbit fits of beta Pic b).
+"""
+
+fit_IAD = True
+
+input_file = os.path.join(orbitize.DATADIR, 'betaPic.csv')
+plx = 19.44
+
+num_secondary_bodies = 1
+data_table = read_input.read_file(input_file)
+
+if fit_IAD:
+    hipparcos_number='027321'
+    fit_secondary_mass=True
+    hipparcos_filename=os.path.join(orbitize.DATADIR, 'HIP027321.d')
+else:
+    hipparcos_number=None
+    fit_secondary_mass=False
+    hipparcos_filename=None
+
+betaPic_system = system.System(
+    num_secondary_bodies, data_table, 1, plx, hipparcos_number=hipparcos_number,
+    hipparcos_filename=hipparcos_filename,
+    fit_secondary_mass=fit_secondary_mass
+)
+
+# set uniform total mass prior
+betaPic_system.sys_priors[-1] = priors.UniformPrior(1.5, 2.0)
+
+# run MCMC
+betaPic_sampler = sampler.MCMC(betaPic_system, num_threads=20, num_temps=20, num_walkers=1000)
+betaPic_sampler.run_sampler(10000000, burn_steps=10000, thin=10)
+
+# save chains
+betaPic_sampler.results.save_results('betaPic_IAD{}.hdf5'.format(fit_IAD))
+
+# make corner plot
+fig = betaPic_sampler.results.plot_corner()
+plt.savefig('corner_IAD{}.png'.format(fit_IAD), dpi=250)
\ No newline at end of file
diff --git a/tests/test_OFTI.py b/tests/test_OFTI.py
index 3f245742..a82bf846 100644
--- a/tests/test_OFTI.py
+++ b/tests/test_OFTI.py
@@ -252,8 +252,8 @@ def test_OFTI_covariances():
 
 
 if __name__ == "__main__":
-    #test_scale_and_rotate()
-    test_run_sampler()
-    test_OFTI_covariances()
-    # test_OFTI_multiplanet()
-    # print("Done!")
+
+    # test_scale_and_rotate()
+    # test_run_sampler()
+    test_OFTI_multiplanet()
+    print("Done!")
diff --git a/tests/test_abs_astrometry.py b/tests/test_abs_astrometry.py
new file mode 100644
index 00000000..8d418665
--- /dev/null
+++ b/tests/test_abs_astrometry.py
@@ -0,0 +1,54 @@
+import numpy as np
+import os
+import astropy.table as table
+import astropy.units as u
+
+import orbitize
+from orbitize import kepler, read_input, system
+
+def test_1planet():
+    """
+    Check that for the 2-body case, the stellar orbit looks as expected
+    """
+
+    # generate a planet orbit
+    sma = 1
+    ecc = 0.1
+    inc = np.radians(45)
+    aop = np.radians(45)
+    pan = np.radians(45)
+    tau = 0.5
+    plx = 1
+    mtot = 1
+    tau_ref_epoch = 0
+    mjup = u.Mjup.to(u.Msun)
+    mass_b = 12 * mjup
+
+    epochs = np.linspace(0, 300, 100) + tau_ref_epoch # nearly the full period, MJD
+
+    ra_model, dec_model, vz_model = kepler.calc_orbit(epochs, sma, ecc, inc, aop, pan, tau, plx, mtot, tau_ref_epoch=tau_ref_epoch)
+
+    # generate some fake measurements just to feed into system.py to test bookkeeping
+    t = table.Table([epochs, np.ones(epochs.shape, dtype=int), ra_model, np.zeros(ra_model.shape), dec_model, np.zeros(dec_model.shape)], 
+                     names=["epoch", "object" ,"raoff", "raoff_err","decoff","decoff_err"])
+    filename = os.path.join(orbitize.DATADIR, "rebound_1planet.csv")
+    t.write(filename)
+
+    # create the orbitize system and generate model predictions using the ground truth
+    astrom_dat = read_input.read_file(filename)
+
+    sys = system.System(1, astrom_dat, mtot, plx, tau_ref_epoch=tau_ref_epoch, fit_secondary_mass=True)
+    sys.track_planet_perturbs = True
+
+    params = np.array([sma, ecc, inc, aop, pan, tau, plx, mass_b, mtot])
+    ra, dec, _ = sys.compute_all_orbits(params)
+
+    # the planet and stellar orbit should just be scaled versions of one another
+    planet_ra = ra[:,1,:]
+    planet_dec = dec[:,1,:]
+    star_ra = ra[:,0,:]
+    star_dec = dec[:,0,:]
+    assert np.all(star_ra + (mass_b / mtot) * planet_ra < 1e-3)
+
+if __name__ == '__main__':
+    test_1planet()
\ No newline at end of file
diff --git a/tests/test_api.py b/tests/test_api.py
index 43a495a4..8445fc91 100644
--- a/tests/test_api.py
+++ b/tests/test_api.py
@@ -28,6 +28,7 @@ def test_compute_model():
 
     params_arr = np.array([1., 0., 0., 0., 0., 245000., 10, 10])
     model, jitter = testSystem_parsing.compute_model(params_arr)
+
     assert model.shape == (4, 2)
 
 
@@ -126,4 +127,8 @@ def test_radec2seppa():
 
 
 if __name__ == "__main__":
+    test_compute_model()
+    test_systeminit()
+    test_chi2lnlike()
     test_custom_likelihood()
+    test_radec2seppa()
diff --git a/tests/test_multiplanet.py b/tests/test_multiplanet.py
index 36d5253b..3d29ca11 100644
--- a/tests/test_multiplanet.py
+++ b/tests/test_multiplanet.py
@@ -32,7 +32,7 @@ def test_compute_model():
     period_c = np.sqrt(c_params[0]**3/mtot)
     period_b = np.sqrt(b_params[0]**3/mtot)
 
-    epochs = np.linspace(0, period_c*365.25, 100) + tau_ref_epoch # the full period of c, MJD
+    epochs = np.linspace(0, period_c*365.25, 5) + tau_ref_epoch # the full period of c, MJD
 
     ra_model, dec_model, vz_model = kepler.calc_orbit(epochs, b_params[0], b_params[1], b_params[2], b_params[3], b_params[4], b_params[5], plx, mtot, tau_ref_epoch=tau_ref_epoch)
 
@@ -172,4 +172,4 @@ def test_fit_selfconsist():
 
 if __name__ == "__main__":
     test_compute_model()
-    #test_fit_selfconsist()
\ No newline at end of file
+    test_fit_selfconsist()
diff --git a/tests/test_multiplanet_rebound.py b/tests/test_multiplanet_rebound.py
index 7bf88319..68abe419 100644
--- a/tests/test_multiplanet_rebound.py
+++ b/tests/test_multiplanet_rebound.py
@@ -147,7 +147,7 @@ def test_2planet_massive():
     b_dec_orb_noc = radec_orb_noc[:,1]
 
     # check that planet c's perturbation is imprinted (nonzero))
-    #assert np.all(b_ra_orb_noc != b_ra_orb)
+    assert np.all(b_ra_orb_noc != b_ra_orb)
 
     # now project the orbit with rebound
     b_manom = basis.tau_to_manom(epochs[0], params[0], params[-1]+params[-3], params[5], tau_ref_epoch)
@@ -411,11 +411,10 @@ def test_2planet_nomass():
 
 
 if __name__ == "__main__":
-    #test_1planet()
+    test_1planet()
     test_2planet_massive()
     test_2planet_massive_reverse_order()
-    #test_2planet_nomass()
-    #test_OFTI_multiplanet()
+    test_2planet_nomass()