diff --git a/phoebe/backend/backends.py b/phoebe/backend/backends.py
index bbbc3883e..1b0b2fadb 100644
--- a/phoebe/backend/backends.py
+++ b/phoebe/backend/backends.py
@@ -961,6 +961,9 @@ def _compute_intrinsic_system_at_t0(self, b, compute,
             if dynamics_method in ['nbody', 'rebound']:
                 t0, xs0, ys0, zs0, vxs0, vys0, vzs0, inst_ds0, inst_Fs0, ethetas0, elongans0, eincls0 = dynamics.nbody.dynamics_from_bundle(b, [t0], compute, return_roche_euler=True, **kwargs)
 
+            elif dynamics_method in ['xyz']:
+                t0, xs0, ys0, zs0, vxs0, vys0, vzs0, inst_ds0, inst_Fs0, ethetas0, elongans0, eincls0 = dynamics.xyz.dynamics_from_bundle(b, [t0], compute, return_roche_euler=True, **kwargs)
+
             elif dynamics_method == 'bs':
                 # TODO: pass stepsize
                 # TODO: pass orbiterror
@@ -1027,6 +1030,9 @@ def _worker_setup(self, b, compute, times, infolists, **kwargs):
             if dynamics_method in ['nbody', 'rebound']:
                 ts, xs, ys, zs, vxs, vys, vzs, inst_ds, inst_Fs, ethetas, elongans, eincls = dynamics.nbody.dynamics_from_bundle(b, times, compute, return_roche_euler=True, **kwargs)
 
+            elif dynamics_method in ['xyz']:
+                ts, xs, ys, zs, vxs, vys, vzs, inst_ds, inst_Fs, ethetas, elongans, eincls = dynamics.xyz.dynamics_from_bundle(b, times, compute, return_roche_euler=True, **kwargs)
+
             elif dynamics_method == 'bs':
                 # if distortion_method == 'roche':
                     # raise ValueError("distortion_method '{}' not compatible with dynamics_method '{}'".format(distortion_method, dynamics_method))
@@ -1039,6 +1045,8 @@ def _worker_setup(self, b, compute, times, infolists, **kwargs):
             elif dynamics_method=='keplerian':
                 # TODO: make sure that this takes systemic velocity and corrects positions and velocities (including ltte effects if enabled)
                 ts, xs, ys, zs, vxs, vys, vzs, ethetas, elongans, eincls = dynamics.keplerian.dynamics_from_bundle(b, times, compute, return_euler=True, **kwargs)
+                inst_ds = None
+                inst_Fs = None
 
             else:
                 raise NotImplementedError
@@ -1054,6 +1062,8 @@ def _worker_setup(self, b, compute, times, infolists, **kwargs):
             vxs, vys, vzs = [np.zeros(len(times))], [np.zeros(len(times))], [np.zeros(len(times))]
             xs, ys, zs = [np.zeros(len(times))], [np.zeros(len(times))], [np.full(len(times), vgamma)]
             ethetas, elongans, eincls = [np.zeros(len(times))], [np.full(len(times), long_an)], [np.full(len(times), incl)]
+            inst_ds = None
+            inst_Fs = None
 
             for i,t in enumerate(times):
                 zs[0][i] = vgamma*(t-t0)
@@ -1065,7 +1075,8 @@ def _worker_setup(self, b, compute, times, infolists, **kwargs):
                     dynamics_method=dynamics_method,
                     ts=ts, xs=xs, ys=ys, zs=zs,
                     vxs=vxs, vys=vys, vzs=vzs,
-                    ethetas=ethetas, elongans=elongans, eincls=eincls)
+                    ethetas=ethetas, elongans=elongans, eincls=eincls,
+                    inst_ds=inst_ds, inst_Fs=inst_Fs)
 
     def _run_single_time(self, b, i, time, infolist, **kwargs):
         logger.debug("rank:{}/{} PhoebeBackend._run_single_time(i={}, time={}, infolist={}, **kwargs.keys={})".format(mpi.myrank, mpi.nprocs, i, time, infolist, kwargs.keys()))
@@ -1096,10 +1107,12 @@ def _run_single_time(self, b, i, time, infolist, **kwargs):
         if True in [info['needs_mesh'] for info in infolist]:
 
             if dynamics_method in ['nbody', 'rebound']:
-                di = dynamics.at_i(inst_ds, i)
-                Fi = dynamics.at_i(inst_Fs, i)
                 # by passing these along to update_positions, volume conservation will
                 # handle remeshing the stars
+                inst_ds = kwargs.get('inst_ds')
+                inst_Fs = kwargs.get('inst_Fs')
+                di = dynamics.at_i(inst_ds, i)
+                Fi = dynamics.at_i(inst_Fs, i)
             else:
                 # then allow d to be determined from orbit and original sma
                 # and F to remain fixed
diff --git a/phoebe/dynamics/__init__.py b/phoebe/dynamics/__init__.py
index 89ee986fe..b712b362f 100644
--- a/phoebe/dynamics/__init__.py
+++ b/phoebe/dynamics/__init__.py
@@ -1,6 +1,7 @@
 
 from . import nbody as nbody
 from . import keplerian as keplerian
+from . import xyz as xyz
 from phoebe import u
 
 
diff --git a/phoebe/dynamics/coord.py b/phoebe/dynamics/coord.py
new file mode 100644
index 000000000..659d68034
--- /dev/null
+++ b/phoebe/dynamics/coord.py
@@ -0,0 +1,134 @@
+
+from numpy import zeros
+
+def coord_j2b(mass, rj, vj):
+    """
+    Converts from Jacobi to barycentric coordinates.
+
+    Rewritten from swift/coord/coord_j2b.f.
+
+    Reference: Levison, H.F., Duncan, M.J., The long-term dynamical behavior of short-period comets. Icarus 108, 18-36, 1994.
+
+    """
+
+    nbod = len(mass)
+    eta = zeros((nbod))
+    rb = zeros((nbod, 3))
+    vb = zeros((nbod, 3))
+
+    # First compute auxiliary variables, then the barycentric positions
+    eta[0] = mass[0]
+    for i in range(1,nbod):
+       eta[i] = eta[i-1] + mass[i]
+
+    i = nbod-1
+    mtot = eta[i]
+    rb[i] = eta[i-1]*rj[i]/mtot
+    vb[i] = eta[i-1]*vj[i]/mtot
+
+    capr = mass[i]*rj[i]/mtot
+    capv = mass[i]*vj[i]/mtot
+
+    for i in range(nbod-2,0,-1):
+        rat = eta[i-1]/eta[i]
+        rb[i] = rat*rj[i] - capr
+        vb[i] = rat*vj[i] - capv
+       
+        rat2 = mass[i]/eta[i]
+        capr += rat2*rj[i]
+        capv += rat2*vj[i]
+	
+    # Now compute the Sun's barycentric position
+    rtmp = zeros((3))
+    vtmp = zeros((3))
+
+    for i in range(1,nbod):
+        rtmp += mass[i]*rb[i]
+        vtmp += mass[i]*vb[i]
+
+    rb[0] = -rtmp/mass[0]
+    vb[0] = -vtmp/mass[0]
+
+    return rb, vb
+
+
+def coord_b2j(mass, rb, vb):
+    """
+    Converts from barycentric to Jacobi coordinates.
+
+    Rewritten from swift/coord/coord_b2j.f.
+
+    Reference: Levison, H.F., Duncan, M.J., The long-term dynamical behavior of short-period comets. Icarus 108, 18-36, 1994.
+
+    """
+
+    nbod = len(mass)
+    eta = zeros((nbod))
+    rj = zeros((nbod, 3))
+    vj = zeros((nbod, 3))
+
+    # First compute auxiliary variables, then the Jacobi positions
+    eta[0] = mass[0]
+    for i in range(1,nbod):
+       eta[i] = eta[i-1] + mass[i]
+
+    sumr = mass[0]*rb[0]
+    sumv = mass[0]*vb[0]
+    capr = rb[0]
+    capv = vb[0]
+
+    for i in range(1, nbod-1):
+        rj[i] = rb[i] - capr
+        vj[i] = vb[i] - capv
+
+        sumr += mass[i]*rb[i]
+        sumv += mass[i]*vb[i]
+        capr = sumr/eta[i]
+        capv = sumv/eta[i]
+
+    rj[nbod-1] = rb[nbod-1] - capr
+    vj[nbod-1] = vb[nbod-1] - capv
+ 
+    return rj, vj
+
+
+def coord_h2j(mass, rh, vh):
+    """
+    Converts from heliocentric to Jacobi coordinates.
+
+    Rewritten from swift/coord/coord_h2j.f.
+
+    Reference: Levison, H.F., Duncan, M.J., The long-term dynamical behavior of short-period comets. Icarus 108, 18-36, 1994.
+
+    """
+
+    nbod = len(mass)
+    eta = zeros((nbod))
+    rj = zeros((nbod, 3))
+    vj = zeros((nbod, 3))
+
+    eta[0] = mass[0]
+    for i in range(1,nbod):
+       eta[i] = eta[i-1] + mass[i]
+
+    rj[1] = rh[1]
+    vj[1] = vh[1]
+
+    sumr = mass[1]*rh[1]
+    sumv = mass[1]*vh[1]
+    capr = sumr/eta[1]
+    capv = sumv/eta[1]
+
+    for i in range(2,nbod):
+        rj[i] = rh[i] - capr
+        vj[i] = vh[i] - capv
+
+        if i < nbod-1:
+             sumr += mass[i]*rh[i]
+             sumv += mass[i]*vh[i]
+             capr = sumr/eta[i]
+             capv = sumv/eta[i]
+
+    return rj, vj
+
+
diff --git a/phoebe/dynamics/geometry.py b/phoebe/dynamics/geometry.py
new file mode 100644
index 000000000..cfa5894a7
--- /dev/null
+++ b/phoebe/dynamics/geometry.py
@@ -0,0 +1,271 @@
+#!/usr/bin/env python3
+
+import numpy as np
+
+from phoebe.dynamics import coord
+from phoebe.dynamics import orbel
+
+def geometry(m, elmts, geometry='hierarchical'):
+    """
+    Convert elements to barycentric coordinates for various geometries.
+   
+    """
+
+    if geometry == 'hierarchical':
+
+        _geometry = geometry_hierarchical
+
+    elif geometry == 'twopairs':
+
+        _geometry = geometry_twopairs
+
+    else:
+        raise NotImplementedError
+
+    nbod = len(m)
+    xs = np.zeros((nbod))
+    ys = np.zeros((nbod))
+    zs = np.zeros((nbod))
+    vxs = np.zeros((nbod))
+    vys = np.zeros((nbod))
+    vzs = np.zeros((nbod))
+
+    rb, vb = _geometry(m, elmts)
+
+    # orientation
+    rb[:,0] *= -1.0
+    rb[:,1] *= -1.0
+    vb[:,0] *= -1.0
+    vb[:,1] *= -1.0
+
+    xs[:] = rb[:,0]
+    ys[:] = rb[:,1]
+    zs[:] = rb[:,2]
+    vxs[:] = vb[:,0]
+    vys[:] = vb[:,1]
+    vzs[:] = vb[:,2]
+
+    return xs, ys, zs, vxs, vys, vzs
+
+def geometry_hierarchical(m, elmts):
+    """
+     _          \                |
+    / \          |               |
+    1 2          3               4
+    \_/          |               |
+                /                |
+    """
+
+    nbod = len(m)
+    rj = np.zeros((nbod, 3))
+    vj = np.zeros((nbod, 3))
+
+    # compute Jacobi coordinates
+    msum = m[0]
+    for j in range(1, nbod):
+        msum += m[j]
+        ialpha = -1
+
+        rj[j], vj[j] = orbel.orbel_el2xv(msum, ialpha, elmts[j-1])
+
+    # convert to barycentric frame
+    rb, vb = coord.coord_j2b(m, rj, vj)
+
+    return rb, vb
+
+def geometry_twopairs(m, elmts):
+    """
+     _          _               \ 
+    / \        / \               |
+    1 2        3 4               5
+    \_/        \_/               |
+                                / 
+    """
+
+    nbod = len(m)
+    rh = np.zeros((nbod, 3))
+    vh = np.zeros((nbod, 3))
+    rj = np.zeros((nbod, 3))
+    vj = np.zeros((nbod, 3))
+
+    # (1+2) pair, 1-centric coordinates
+    msum = m[0]+m[1]
+    ialpha = -1
+    r2_1, v2_1 = orbel.orbel_el2xv(msum, ialpha, elmts[0])
+
+    # barycenter
+    r12_1 = m[1]/msum * r2_1
+    v12_1 = m[1]/msum * v2_1
+
+    # (3+4) pair, 3-centric
+    msum = m[2]+m[3]
+    r4_3, v4_3 = orbel.orbel_el2xv(msum, ialpha, elmts[1])
+
+    # barycenter
+    r34_3 = m[3]/msum * r4_3
+    v34_3 = m[3]/msum * v4_3
+
+    # (1+2)+(3+4) mutual orbit, (1+2)-centric
+    msum = m[0]+m[1]+m[2]+m[3]
+    r34_12, v34_12 = orbel.orbel_el2xv(msum, ialpha, elmts[2])
+
+    # everything to 1-centric
+    rh[0,:] = 0.0
+    vh[0,:] = 0.0
+    rh[1,:] = r2_1
+    vh[1,:] = v2_1
+    rh[2,:] = r12_1 + r34_12 - r34_3
+    vh[2,:] = v12_1 + v34_12 - v34_3
+    rh[3,:] = rh[2,:] + r4_3
+    vh[3,:] = vh[2,:] + v4_3
+
+    # everything to Jacobian
+    rj[0:4], vj[0:4] = coord.coord_h2j(m[0:4], rh[0:4], vh[0:4])
+
+    # other bodies (also Jacobian)
+    for j in range(4, nbod):
+        msum += m[j]
+
+        rj[j], vj[j] = orbel.orbel_el2xv(msum, ialpha, elmts[j-1])
+
+    # convert to barycentric frame
+    rb, vb = coord.coord_j2b(m, rj, vj)
+
+    return rb, vb
+
+########################################################################
+
+def invgeometry(m, rb, vb, geometry='hierarchical'):
+    """
+    Convert barycentric coordinates to elements for various geometries.
+   
+    """
+
+    if geometry == 'hierarchical':
+
+        _invgeometry = invgeometry_hierarchical
+
+    elif geometry == 'twopairs':
+
+        _invgeometry = invgeometry_twopairs
+
+    else:
+        raise NotImplementedError
+
+    # cf. mutable numpy arrays
+    rb = np.copy(rb)
+    vb = np.copy(vb)
+
+    # orientation
+    rb[:,0] *= -1.0
+    rb[:,1] *= -1.0
+    vb[:,0] *= -1.0
+    vb[:,1] *= -1.0
+
+    return _invgeometry(m, rb, vb)
+
+def invgeometry_hierarchical(m, rb, vb):
+    """
+     _          \                |
+    / \          |               |
+    1 2          3               4
+    \_/          |               |
+                /                |
+    """
+
+    nbod = len(m)
+    elmts = np.zeros((nbod-1, 6))
+    euler = np.zeros((nbod, 3))
+    roche = np.zeros((nbod, 2))
+
+    # convert to Jacobian coordinates
+    rj, vj = coord.coord_b2j(m, rb, vb)
+
+    # compute osculating elements
+    msum = m[0]
+    for j in range(1, nbod):
+        msum += m[j]
+        ialpha = -1
+
+        elmts[j-1,:] = orbel.orbel_xv2el(msum, rj[j], vj[j])
+        euler[j,:] = orbel.get_euler(elmts[j-1])
+        roche[j,:] = orbel.get_roche(msum, elmts[j-1])
+
+        if j==1:
+            euler[0,:] = euler[1,:]
+            roche[0,:] = roche[1,:]
+        if j>=1:
+            euler[j,0] += np.pi
+
+    return elmts, euler, roche
+
+def invgeometry_twopairs(m, rb, vb):
+    """
+     _          _               \ 
+    / \        / \               |
+    1 2        3 4               5
+    \_/        \_/               |
+                                / 
+    """
+
+    nbod = len(m)
+    elmts = np.zeros((nbod-1, 6))
+    euler = np.zeros((nbod, 3))
+    roche = np.zeros((nbod, 2))
+
+    # (1+2) pair, 1-centric coordinates
+    msum = m[0]+m[1]
+    r2_1 = rb[1] - rb[0]
+    v2_1 = vb[1] - vb[0]
+
+    elmts[0,:] = orbel.orbel_xv2el(msum, r2_1, v2_1)
+    euler[0,:] = orbel.get_euler(elmts[0])
+    roche[0,:] = orbel.get_roche(msum, elmts[0])
+
+    euler[1,:] = euler[0,:]
+    roche[1,:] = roche[0,:]
+    euler[1,0] += np.pi
+
+    # barycenter
+    r12 = (m[0]*rb[0] + m[1]*rb[1])/msum
+    v12 = (m[0]*vb[0] + m[1]*vb[1])/msum
+
+    # (3+4) pair, 3-centric
+    msum = m[2]+m[3]
+    r4_3 = rb[3] - rb[2]
+    v4_3 = vb[3] - vb[2]
+
+    elmts[1,:] = orbel.orbel_xv2el(msum, r4_3, v4_3)
+    euler[2,:] = orbel.get_euler(elmts[1])
+    roche[2,:] = orbel.get_roche(msum, elmts[1])
+    
+    euler[3,:] = euler[2,:]
+    roche[3,:] = roche[2,:]
+    euler[3,0] += np.pi
+
+    # barycenter
+    r34 = (m[2]*rb[2] + m[3]*rb[3])/msum
+    v34 = (m[2]*vb[2] + m[3]*vb[3])/msum
+
+    # (1+2)+(3+4) mutual orbit, (1+2)-centric
+    msum = m[0]+m[1]+m[2]+m[3]
+    r34_12 = r34 - r12
+    v34_12 = v34 - v12
+    elmts[2,:] = orbel.orbel_xv2el(msum, r34_12, v34_12)
+
+    # everything to Jacobian
+    rj, vj = coord.coord_b2j(m, rb, vb)
+    
+    # other bodies (also Jacobian)
+    for j in range(4, nbod):
+        msum += m[j]
+
+        elmts[j-1,:] = orbel.orbel_xv2el(msum, rj[j], vj[j])
+        euler[j,:] = orbel.get_euler(elmts[j-1])
+        roche[j,:] = orbel.get_roche(msum, elmts[j-1])
+
+        euler[j,0] += np.pi
+
+    return elmts, euler, roche
+
+
diff --git a/phoebe/dynamics/keplerian.py b/phoebe/dynamics/keplerian.py
index e6f6c199b..3d5e7fff7 100644
--- a/phoebe/dynamics/keplerian.py
+++ b/phoebe/dynamics/keplerian.py
@@ -1,11 +1,12 @@
-
+#!/usr/bin/env python3
 
 import numpy as np
-from numpy import pi,sqrt,cos,sin,tan,arctan
 from scipy.optimize import newton
 
 from phoebe import u, c
-from phoebe import conf
+
+from phoebe.dynamics import coord
+from phoebe.dynamics import orbel
 
 import logging
 logger = logging.getLogger("DYNAMICS.KEPLERIAN")
@@ -47,411 +48,190 @@ def dynamics_from_bundle(b, times, compute=None, return_euler=False, **kwargs):
     else:
         ltte = False
 
-    # make sure times is an array and not a list
     times = np.array(times)
 
-    vgamma = b.get_value(qualifier='vgamma', context='system', unit=u.solRad/u.d, **_skip_filter_checks)
-    t0 = b.get_value(qualifier='t0', context='system', unit=u.d, **_skip_filter_checks)
-
     hier = b.hierarchy
     starrefs = hier.get_stars()
     orbitrefs = hier.get_orbits()
-    s = b.filter(context='component', **_skip_filter_checks)
-
-    periods, eccs, smas, t0_perpasses, per0s, long_ans, incls, dpdts, \
-    deccdts, dperdts, components = [],[],[],[],[],[],[],[],[],[],[]
-
-
-    for component in starrefs:
-
-        # we need to build a list of all orbitlabels underwhich this component
-        # belongs.  For a simple binary this is just the parent, but for hierarchical
-        # systems we need to get the labels of the outer-orbits as well
-        ancestororbits = []
-        comp = component
-        while hier.get_parent_of(comp) in orbitrefs:
-            comp = hier.get_parent_of(comp)
-            ancestororbits.append(comp)
-
-        #print "***", component, ancestororbits
-
-        periods.append([s.get_value(qualifier='period_anom', unit=u.d, component=orbit, **_skip_filter_checks) for orbit in ancestororbits])
-        eccs.append([s.get_value(qualifier='ecc', component=orbit, **_skip_filter_checks) for orbit in ancestororbits])
-        t0_perpasses.append([s.get_value(qualifier='t0_perpass', unit=u.d, component=orbit, **_skip_filter_checks) for orbit in ancestororbits])
-        per0s.append([s.get_value(qualifier='per0', unit=u.rad, component=orbit, **_skip_filter_checks) for orbit in ancestororbits])
-        long_ans.append([s.get_value(qualifier='long_an', unit=u.rad, component=orbit, **_skip_filter_checks) for orbit in ancestororbits])
-        incls.append([s.get_value(qualifier='incl', unit=u.rad, component=orbit, **_skip_filter_checks) for orbit in ancestororbits])
-        # TODO: do we need to convert from anomalistic to sidereal?
-        dpdts.append([s.get_value(qualifier='dpdt', unit=u.d/u.d, component=orbit, **_skip_filter_checks) for orbit in ancestororbits])
-        if conf.devel:
-            try:
-                deccdts.append([s.get_value('deccdt', u.dimensionless_unscaled/u.d, component=orbit, **_skip_filter_checks) for orbit in ancestororbits])
-            except ValueError:
-                deccdts.append([0.0 for orbit in ancestororbits])
-        else:
-            deccdts.append([0.0 for orbit in ancestororbits])
-        dperdts.append([s.get_value(qualifier='dperdt', unit=u.rad/u.d, component=orbit, **_skip_filter_checks) for orbit in ancestororbits])
-
-        # sma needs to be the COMPONENT sma.  This is stored in the bundle for stars, but is NOT
-        # for orbits in orbits, so we'll need to recompute those from the mass-ratio and sma of
-        # the parent orbit.
-
-        smas_this = []
-        for comp in [component]+ancestororbits[:-1]:
-            if comp in starrefs:
-                smas_this.append(s.get_value(qualifier='sma', unit=u.solRad, component=comp, **_skip_filter_checks))
-            else:
-                q = s.get_value(qualifier='q', component=hier.get_parent_of(comp), **_skip_filter_checks)
-                comp_comp = hier.get_primary_or_secondary(comp)
 
-                # NOTE: similar logic is also in constraints.comp_sma
-                # If changing any of the logic here, it should be changed there as well.
-                if comp_comp == 'primary':
-                    qthing = (1. + 1./q)
-                else:
-                    qthing = (1. + q)
+    G = c.G.to('solRad3 / (Msun d2)').value
+    masses = [b.get_value(qualifier='mass', unit=u.solMass, component=component, context='component', **_skip_filter_checks) * G for component in starrefs]
+    smas = [b.get_value(qualifier='sma', unit=u.solRad, component=component, context='component', **_skip_filter_checks) for component in orbitrefs]
+    eccs = [b.get_value(qualifier='ecc', component=component, context='component', **_skip_filter_checks) for component in orbitrefs]
+    incls = [b.get_value(qualifier='incl', unit=u.rad, component=component, context='component', **_skip_filter_checks) for component in orbitrefs]
+    per0s = [b.get_value(qualifier='per0', unit=u.rad, component=component, context='component', **_skip_filter_checks) for component in orbitrefs]
+    long_ans = [b.get_value(qualifier='long_an', unit=u.rad, component=component, context='component', **_skip_filter_checks) for component in orbitrefs]
+    t0_perpasses = [b.get_value(qualifier='t0_perpass', unit=u.d, component=component, context='component', **_skip_filter_checks) for component in orbitrefs]
+    periods = [b.get_value(qualifier='period', unit=u.d, component=component, context='component', **_skip_filter_checks) for component in orbitrefs]
+    mean_anoms = [b.get_value(qualifier='mean_anom', unit=u.rad, component=component, context='component', **_skip_filter_checks) for component in orbitrefs]
 
-                smas_this.append(s.get_value(qualifier='sma', unit=u.solRad, component=hier.get_parent_of(comp), **_skip_filter_checks) / qthing)
-
-        smas.append(smas_this)
-
-        # components is whether an entry is the primary or secondary in its parent orbit, so here we want
-        # to start with component and end one level short of the top-level orbit
-        components.append([hier.get_primary_or_secondary(component=comp) for comp in [component]+ancestororbits[:-1]])
-
-
-    return  dynamics(times, periods, eccs, smas, t0_perpasses, per0s, \
-                    long_ans, incls, dpdts, deccdts, dperdts, \
-                    components, t0, vgamma, \
-                    mass_conservation=True, ltte=ltte, return_euler=return_euler)
+    vgamma = b.get_value(qualifier='vgamma', context='system', unit=u.solRad/u.d, **_skip_filter_checks)
+    t0 = b.get_value(qualifier='t0', context='system', unit=u.d, **_skip_filter_checks)
 
+    dpdts = [b.get_value(qualifier='dpdt', unit=u.d/u.d, component=component, context='component', **_skip_filter_checks) for component in orbitrefs]
+    dperdts = [b.get_value(qualifier='dperdt', unit=u.rad/u.d, component=component, context='component', **_skip_filter_checks) for component in orbitrefs]
 
+    return dynamics(times, masses, smas, eccs, incls, per0s, long_ans, mean_anoms, \
+        dpdts, dperdts, \
+        t0, vgamma, ltte, \
+        return_euler=return_euler)
 
 
-def dynamics(times, periods, eccs, smas, t0_perpasses, per0s, long_ans, incls,
-            dpdts, deccdts, dperdts, components, t0=0.0, vgamma=0.0,
-            mass_conservation=True, ltte=False, return_euler=False):
+def dynamics(times, masses, smas, eccs, incls, per0s, long_ans, mean_anoms, \
+    dpdts, dperdts, \
+    t0=0.0, vgamma=0.0, ltte=False, \
+    return_euler=False):
     """
-    Compute the positions and velocites of each star in their nested
-    Keplerian orbits at a given list of times.
+    Computes Keplerian dynamics, including binaries, triples, quadruples...
+
+     _          \                |
+    / \          |               |
+    1 2          3               4
+    \_/          |               |
+                /                |
 
     See :func:`dynamics_from_bundle` for a wrapper around this function
     which automatically handles passing everything in the correct order
     and in the correct units.
 
+    Note: orbits = stars - 1
+
     Args:
-        times: (iterable) times at which to compute positions and
-            velocities for each star
-        periods: (iterable) anomalistic period of the parent orbit for each star
-            [days]
-        eccs: (iterable) eccentricity of the parent orbit for each star
-        smas: (iterable) semi-major axis of the parent orbit for each
-            star [solRad]
-        t0_perpasses: (iterable) t0_perpass of the parent orbit for each
-            star [days]
-        per0s: (iterable) longitudes of periastron of the parent orbit
-            for each star [rad]
-        long_ans: (iterable) longitudes of the ascending node of the
-            parent orbit for each star [rad]
-        incls: (iterable) inclination of the parent orbit for each
-            star [rad]
-        dpdts: (iterable) change in period with respect to time of the
-            parent orbit for each star [days/day]
-        deccdts: (iterable) change in eccentricity with respect to time
-            of the parent orbit for each star [1/day]
-        dperdts: (iterable) change in periastron with respect to time
-            of the parent orbit for each star [rad/d]
-        components: (iterable) component ('primary' or 'secondary') of
-            each star within its parent orbit [string]
-        t0: (float, default=0) time at which all initial values (ie period, per0)
-            are given [days]
-        mass_conservation: (bool, optional) whether to require mass
-            conservation if any of the derivatives (dpdt, dperdt, etc)
-            are non-zero [default: True]
-        return_euler: (bool, default=False) whether to include euler angles
-            in the return
+        times: (iterable) times when to compute positions and velocities [days]
+        masses: (iterable) mass for each star [solMass]
+        smas: (iterable) semi-major axis for each orbit [solRad]
+        eccs: (iterable) eccentricity for each orbit [1]
+        incls: (iterable) inclination for each orbit [rad]
+        per0s: (iterable) longitudes of periastron for each orbit [rad]
+        long_ans: (iterable) longitude of the ascending node for each orbit [rad]
+        mean_anoms: (iterable) mean anomaly for each orbit [rad]
+        dpdts: (iterable) change in period with respect to time [days/day]
+        dperdts: (iterable) change in periastron with respect to time [deg/day]
+        t0: (float, default=0) epoch at which elements were given [days]
+        return_euler: (bool, default=False) whether to return euler angles
 
     Returns:
-        t, xs, ys, zs, vxs, vys, vzs [, theta, longan, incl].
-        t is a numpy array of all times,
-        the remaining are a list of numpy arrays (a numpy array per
-        star - in order given by b.hierarchy.get_stars()) for the cartesian
-        positions and velocities of each star at those same times.
-        Euler angles (theta, longan, incl) are only returned if return_euler is
-        set to True.
-    """
-    # TODO: NOTE: smas must be per-component, not per-orbit
-    # TODO: steal some documentation from 2.0a:keplerorbit.py:get_orbit
-    # TODO: deal with component number more smartly
-
-    def binary_dynamics(times, period, ecc, sma, t0_perpass, per0, long_an,
-                        incl, dpdt, deccdt, dperdt, component='primary',
-                        t0=0.0, vgamma=0.0, mass_conservation=True,
-                        com_pos=(0.,0.,0.), com_vel=(0.,0.,0.), com_euler=(0.,0.,0.)):
-        """
-        """
-        # TODO: steal some documentation from 2.0a:keplerorbit.py:get_orbit
-
-
-        #-- if dpdt is non-zero, the period is actually an array, and the semi-
-        #   major axis changes to match Kepler's third law (unless
-        #   `mass_conservation` is set to False)
-        # p0 = period
-        # sma0 = sma
-
-        if dpdt!=0:
-            p0 = period
-            period = dpdt*(times-t0) + p0
-            if mass_conservation:
-                # Pieter used to have a if not np.isscale(period) to use period[0] instead of p0,
-                # effectively starting mass conservation at the first dataset time instead of t0
-                sma = sma/p0**2*period**2
-
-        #-- if dperdt is non-zero, the argument of periastron is actually an
-        #   array
-        if dperdt!=0.:
-            per0 = dperdt*(times-t0) + per0
-        #-- if deccdt is non-zero, the eccentricity is actually an array
-        if deccdt!=0:
-            ecc = deccdt*(times-t0) + ecc
-
-
-        #-- compute orbit
-        n = 2*pi/period
-        ma = n*(times-t0_perpass)
-        E,theta = _true_anomaly(ma,ecc)
-        r = sma*(1-ecc*cos(E))
-        PR = r*sin(theta)
-        #-- compute rdot and thetadot
-        l = r*(1+ecc*cos(theta))#-omega))
-        L = 2*pi*sma**2/period*sqrt(1-ecc**2)
-        rdot = L/l*ecc*sin(theta)#-omega)
-        thetadot = L/r**2
-        #-- the secondary is half an orbit further than the primary
-        if 'sec' in component.lower():
-            theta += pi
-        #-- take care of Euler angles
-        theta_ = theta+per0
-        #-- convert to the right coordinate frame
-        #-- create some shortcuts
-        sin_theta_ = sin(theta_)
-        cos_theta_ = cos(theta_)
-        sin_longan = sin(long_an)
-        cos_longan = cos(long_an)
-        #-- spherical coordinates to cartesian coordinates. Note that we actually
-        #   set incl=-incl (but of course it doesn't show up in the cosine). We
-        #   do this to match the convention that superior conjunction (primary
-        #   eclipsed) happens at periastron passage when per0=90 deg.
-        x = r*(cos_longan*cos_theta_ - sin_longan*sin_theta_*cos(incl))
-        y = r*(sin_longan*cos_theta_ + cos_longan*sin_theta_*cos(incl))
-        z = r*(sin_theta_*sin(-incl))
-        #-- spherical vectors to cartesian vectors, and then rotated for
-        #   the Euler angles Omega and i.
-        vx_ = cos_theta_*rdot - sin_theta_*r*thetadot
-        vy_ = sin_theta_*rdot + cos_theta_*r*thetadot
-        vx = cos_longan*vx_ - sin_longan*vy_*cos(incl)
-        vy = sin_longan*vx_ + cos_longan*vy_*cos(incl)
-        vz = sin(-incl)*vy_
-        #-- that's it!
-
-        # correct by vgamma (only z-direction)
-        # NOTE: vgamma is in the direction of positive RV or negative vz
-        vz -= vgamma
-        z -= vgamma * (times-t0)
-
-        return (x+com_pos[0],y+com_pos[1],z+com_pos[2]),\
-                (vx+com_vel[0],vy+com_vel[1],vz+com_vel[2]),\
-                (theta_,long_an,incl)
-                # (theta_+com_euler[0],long_an+com_euler[1],incl+com_euler[2])
-
-    def binary_dynamics_nested(times, periods, eccs, smas, \
-                            t0_perpasses, per0s, long_ans, incls, dpdts, deccdts, \
-                            dperdts, components, t0, vgamma, \
-                            mass_conservation):
-
-        """
-        compute the (possibly nested) positions of a single component (ltte should be
-        handle externally)
-        """
-
-
-        if not hasattr(periods, '__iter__'):
-            # then we don't have to worry about nesting, and each item should
-            # be a single value ready to pass to binary_dynamics
-            pos, vel, euler = binary_dynamics(times, periods, eccs, smas, t0_perpasses, \
-                                per0s, long_ans, incls, dpdts, deccdts, dperdts, components, \
-                                t0, vgamma, mass_conservation)
-
-        else:
-            # Handle nesting - if this star is not in the top-level orbit, then
-            # all values should actually be lists.  We want to sort by period to handle
-            # the outer orbits first and then apply those offsets to the inner-orbit(s)
-
-            # let's convert to arrays so we can use argsort easily
-            periods = np.array(periods)
-            eccs = np.array(eccs)
-            smas = np.array(smas)
-            t0_perpasses = np.array(t0_perpasses)
-            per0s = np.array(per0s)
-            long_ans = np.array(long_ans)
-            incls = np.array(incls)
-            dpdts = np.array(dpdts)
-            deccdts = np.array(deccdts)
-            dperdts = np.array(dperdts)
-            components = np.array(components)
-
-            si = periods.argsort()[::-1]
-
-            #print "***", periods, si
-
-            pos = (0.0, 0.0, 0.0)
-            vel = (0.0, 0.0, 0.0)
-            euler = (0.0, 0.0, 0.0)
-
-            for period, ecc, sma, t0_perpass, per0, long_an, incl, dpdt, \
-                deccdt, dperdt, component in zip(periods[si], eccs[si], \
-                smas[si], t0_perpasses[si], per0s[si], long_ans[si], \
-                incls[si], dpdts[si], deccdts[si], dperdts[si], components[si]):
-
-                pos, vel, euler = binary_dynamics(times, period, ecc, sma, t0_perpass, \
-                                    per0, long_an, incl, dpdt, deccdt, dperdt, component, \
-                                    t0, vgamma, mass_conservation,
-                                    com_pos=pos, com_vel=vel, com_euler=euler)
-
-
-
-        return pos, vel, euler
-
-
-
-
-    xs, ys, zs = [], [], []
-    vxs, vys, vzs = [], [], []
+        t: is a numpy array of all times,
+        xs, ys, zs: Cartesian positions [solRad]
+        vxs, vys, vzs: Cartesian velocities [solRad/day]
+        theta, longan, incl: Euler angles [rad]
 
-    if return_euler:
-        ethetas, elongans, eincls = [], [], []
-
-    for period, ecc, sma, t0_perpass, per0, long_an, incl, dpdt, deccdt, \
-            dperdt, component in zip(periods, eccs, smas, t0_perpasses, per0s, long_ans, \
-            incls, dpdts, deccdts, dperdts, components):
-
-        # We now have the orbital parameters for a single star/component.
+    Note: Euler angles are only returned if return_euler is True.
 
-        if ltte:
-            #scale_factor = 1.0/c.c.value * c.R_sun.value/(24*3600.)
-            scale_factor = (c.R_sun/c.c).to(u.d).value
-
-            def propertime_barytime_residual(t, time):
-                pos, vel, euler = binary_dynamics_nested(t, period, ecc, sma, \
-                                t0_perpass, per0, long_an, incl, dpdt, deccdt, \
-                                dperdt, components=component, t0=t0, vgamma=vgamma, \
-                                mass_conservation=mass_conservation)
-                z = pos[2]
-                return t - z*scale_factor - time
+    """
 
-            # Finding that right time is easy with a Newton optimizer:
-            propertimes = [newton(propertime_barytime_residual, time, args=(time,)) for \
-               time in times]
-            propertimes = np.array(propertimes).ravel()
+    # a part called repeatedly (cf. ltte)
+    def _keplerian(times, component=None):
+
+        for i, t in enumerate(times):
+
+            # compute Jacobi coordinates
+            msum = masses[0]
+            for j in range(1, nbod):
+                msum += masses[j]
+                ialpha = -1
+                a = smas[j-1]
+                n = np.sqrt(msum/a**3)
+                M = mean_anoms[j-1] + n*(t-t0)
+                P = 2.0*np.pi/n
+                dpdt = dpdts[j-1]
+                omega = per0s[j-1]
+                dperdt = dperdts[j-1]
+
+                if dpdt != 0.0:
+                    P_ = P + dpdt*(t-t0)
+                    a *= (P_/P)**2
+                    M -= n*dpdt/(2.*P)*(t-t0)**2
+
+                if dperdt != 0.0:
+                    omega += dperdt*(t-t0)
+
+                elmts = [a, eccs[j-1], incls[j-1], long_ans[j-1], omega, M]
+
+                rj[j], vj[j] = orbel.orbel_el2xv(msum, ialpha, elmts)
+
+                # compute Euler angles
+                # need to copy the primary
+                # need to add np.pi for the secondary
+                if return_euler:
+                    euler[j] = orbel.get_euler(elmts)
+
+                    if j==1:
+                        euler[0,:] = euler[1,:]
+                    if j>=1:
+                        euler[j,0] += np.pi
+
+            # convert to barycentric frame
+            rb, vb = coord.coord_j2b(masses, rj, vj)
+
+            # gamma velocity
+            rb[:,2] -= vgamma*(t-t0)
+            vb[:,2] -= vgamma
+
+            # orientation
+            rb[:,0] *= -1.0
+            rb[:,1] *= -1.0
+            vb[:,0] *= -1.0
+            vb[:,1] *= -1.0
+
+            # update only selected components
+            k = range(0, nbod)
+            if component != None:
+                k = component
+
+            xs[k,i] = rb[k,0]
+            ys[k,i] = rb[k,1]
+            zs[k,i] = rb[k,2]
+            vxs[k,i] = vb[k,0]
+            vys[k,i] = vb[k,1]
+            vzs[k,i] = vb[k,2]
+            ethetas[k,i] = euler[k,0]
+            elongans[k,i] = euler[k,1]
+            eincls[k,i] = euler[k,2]
+
+        return
+
+
+    # ...
+    nbod = len(masses)
+    rj = np.zeros((nbod, 3))
+    vj = np.zeros((nbod, 3))
+    euler = np.zeros((nbod, 3))
+
+    xs = np.zeros((nbod, len(times)))
+    ys = np.zeros((nbod, len(times)))
+    zs = np.zeros((nbod, len(times)))
+    vxs = np.zeros((nbod, len(times)))
+    vys = np.zeros((nbod, len(times)))
+    vzs = np.zeros((nbod, len(times)))
 
-            pos, vel, euler = binary_dynamics_nested(propertimes, period, ecc, sma, \
-                            t0_perpass, per0, long_an, incl, dpdt, deccdt, \
-                            dperdt, components=component, t0=t0, vgamma=vgamma, \
-                            mass_conservation=mass_conservation)
-
-
-        else:
-            pos, vel, euler = binary_dynamics_nested(times, period, ecc, sma, \
-                    t0_perpass, per0, long_an, incl, dpdt, deccdt, \
-                    dperdt, components=component, t0=t0, vgamma=vgamma, \
-                    mass_conservation=mass_conservation)
-
-
-        xs.append(pos[0])
-        ys.append(pos[1])
-        zs.append(pos[2])
-        vxs.append(vel[0])
-        vys.append(vel[1])
-        vzs.append(vel[2])
-        if return_euler:
-            ethetas.append(euler[0])
-            elongans.append([euler[1]]*len(euler[0]))
-            eincls.append([euler[2]]*len(euler[0]))
-
-
-    # if return_euler:
-    #     return times, \
-    #         xs*u.solRad, ys*u.solRad, zs*u.solRad, \
-    #         vxs*u.solRad/u.d, vys*u.solRad/u.d, vzs*u.solRad/u.d, \
-    #         ethetas*u.rad, elongans*u.rad, eincls*u.rad
-    # else:
-    #     return times, \
-    #         xs*u.solRad, ys*u.solRad, zs*u.solRad, \
-    #         vxs*u.solRad/u.d, vys*u.solRad/u.d, vzs*u.solRad/u.d    if return_euler:
-
-    # d, solRad, solRad/d, rad
     if return_euler:
-        return times, \
-            xs, ys, zs, \
-            vxs, vys, vzs, \
-            ethetas, elongans, eincls
-    else:
-        return times, \
-            xs, ys, zs, \
-            vxs, vys, vzs
-
-
-
-
-def _true_anomaly(M,ecc,itermax=8):
-    r"""
-    Calculation of true and eccentric anomaly in Kepler orbits.
+        ethetas  = np.zeros((nbod, len(times)))
+        elongans = np.zeros((nbod, len(times)))
+        eincls   = np.zeros((nbod, len(times)))
 
-    ``M`` is the phase of the star, ``ecc`` is the eccentricity
+    # unfortunately, ltte must be computed for each * (separately)!
+    if ltte:
+        scale_factor = (u.solRad/c.c).to(u.d).value
 
-    See p.39 of Hilditch, 'An Introduction To Close Binary Stars':
+        for j in range(0,nbod):
 
-    Kepler's equation:
+            def propertime_residual(t, time):
+                _keplerian([t], component=j)
+                z = zs[j][0]
+                return t - z*scale_factor - time
 
-    .. math::
+            propertimes = [newton(propertime_residual, time, args=(time,)) for time in times]
+            propertimes = np.array(propertimes).ravel()
 
-        E - e\sin E = \frac{2\pi}{P}(t-T)
+            _keplerian(propertimes, component=j)
 
-    with :math:`E` the eccentric anomaly. The right hand size denotes the
-    observed phase :math:`M`. This function returns the true anomaly, which is
-    the position angle of the star in the orbit (:math:`\theta` in Hilditch'
-    book). The relationship between the eccentric and true anomaly is as
-    follows:
+    else:
+        _keplerian(times)
 
-    .. math::
+    if return_euler:
+        return times, xs, ys, zs, vxs, vys, vzs, ethetas, elongans, eincls
+    else:
+        return times, xs, ys, zs, vxs, vys, vzs
 
-        \tan(\theta/2) = \sqrt{\frac{1+e}{1-e}} \tan(E/2)
 
-    @parameter M: phase
-    @type M: float
-    @parameter ecc: eccentricity
-    @type ecc: float
-    @keyword itermax: maximum number of iterations
-    @type itermax: integer
-    @return: eccentric anomaly (E), true anomaly (theta)
-    @rtype: float,float
-    """
-    # Initial value
-    Fn = M + ecc*sin(M) + ecc**2/2.*sin(2*M)
-
-    # Iterative solving of the transcendent Kepler's equation
-    for i in range(itermax):
-        F = Fn
-        Mn = F-ecc*sin(F)
-        Fn = F+(M-Mn)/(1.-ecc*cos(F))
-        keep = F!=0 # take care of zerodivision
-        if hasattr(F,'__iter__'):
-            if np.all(abs((Fn-F)[keep]/F[keep])<0.00001):
-                break
-        elif (abs((Fn-F)/F)<0.00001):
-            break
-
-    # relationship between true anomaly (theta) and eccentric anomaly (Fn)
-    true_an = 2.*arctan(sqrt((1.+ecc)/(1.-ecc))*tan(Fn/2.))
-
-    return Fn,true_an
diff --git a/phoebe/dynamics/nbody.py b/phoebe/dynamics/nbody.py
index ef05de081..fc584bef3 100644
--- a/phoebe/dynamics/nbody.py
+++ b/phoebe/dynamics/nbody.py
@@ -6,6 +6,7 @@
 
 
 from phoebe import u, c
+from phoebe import conf
 
 try:
     import photodynam
@@ -94,6 +95,7 @@ def dynamics_from_bundle(b, times, compute=None, return_roche_euler=False, use_k
     ltte = computeps.get_value(qualifier='ltte', ltte=kwargs.get('ltte', None), **_skip_filter_checks)
     gr = computeps.get_value(qualifier='gr', gr=kwargs.get('gr', None), **_skip_filter_checks)
     integrator = computeps.get_value(qualifier='integrator', integrator=kwargs.get('integrator', None), **_skip_filter_checks)
+    epsilon = computeps.get_value(qualifier='epsilon', epsilon=kwargs.get('epsilon', None), **_skip_filter_checks)
 
     starrefs = hier.get_stars()
     orbitrefs = hier.get_orbits() if use_kepcart else [hier.get_parent_of(star) for star in starrefs]
@@ -127,12 +129,14 @@ def mean_anom(t0, t0_perpass, period):
 
     return dynamics(times, masses, smas, eccs, incls, per0s, long_ans, \
                     mean_anoms, rotperiods, t0, vgamma, stepsize, ltte, gr,
-                    integrator, use_kepcart=use_kepcart, return_roche_euler=return_roche_euler)
+                    integrator, use_kepcart=use_kepcart, return_roche_euler=return_roche_euler,
+                    epsilon=epsilon)
 
 
 def dynamics(times, masses, smas, eccs, incls, per0s, long_ans, mean_anoms,
         rotperiods=None, t0=0.0, vgamma=0.0, stepsize=0.01, ltte=False, gr=False,
-        integrator='ias15', return_roche_euler=False, use_kepcart=False):
+        integrator='ias15', return_roche_euler=False, use_kepcart=False,
+        epsilon=1.0e-9):
 
     if not _can_rebound:
         raise ImportError("rebound is not installed")
@@ -165,16 +169,24 @@ def residual(t):
 
     sim = rebound.Simulation()
 
+    # TODO: switch between different GR setups based on masses/hierarchy
+    # http://reboundx.readthedocs.io/en/latest/effects.html#general-relativity
     if gr:
-        logger.info("enabling 'gr_full' in reboundx")
+        logger.info("enabling 'gr' in reboundx")
         rebx = reboundx.Extras(sim)
-        # TODO: switch between different GR setups based on masses/hierarchy
-        # http://reboundx.readthedocs.io/en/latest/effects.html#general-relativity
-        params = rebx.add_gr_full()
+        gr = rebx.load_force("gr")
+        gr.params["c"] = c.c.to("AU/d").value
+        rebx.add_force(gr)
 
     sim.integrator = integrator
     # NOTE: according to rebound docs: "stepsize will change for adaptive integrators such as IAS15"
     sim.dt = stepsize
+    sim.ri_ias15.epsilon = epsilon
+    sim.ri_whfast.corrector = 17
+    sim.ri_whfast.safe_mode = 0;
+    sim.G = 1.0
+    if conf.devel:
+        sim.status()
 
     if use_kepcart:
         # print "*** bs.kep2cartesian", masses, smas, eccs, incls, per0s, long_ans, mean_anoms, t0
@@ -292,7 +304,7 @@ def residual(t):
 
                 # get the orbit based on the primary component defined already
                 # in the simulation.
-                orbit = particle.calculate_orbit()
+                orbit = particle.orbit()
 
                 # for instantaneous separation, we need the current separation
                 # from the sibling component in units of its instantaneous (?) sma
@@ -302,9 +314,12 @@ def residual(t):
                 # on the INSTANTANEOUS orbital PERIOD.
                 Fs[j][i] = orbit.P / rotperiods[j]
 
-                # TODO: need to add np.pi for secondary component
                 ethetas[j][i] = orbit.f + orbit.omega # true anomaly + periastron
 
+                # need to add np.pi for secondary component
+                if j==1:
+                    ethetas[j][i] += np.pi
+
                 elongans[j][i] = orbit.Omega
 
                 eincls[j][i] = orbit.inc
diff --git a/phoebe/dynamics/orbel.py b/phoebe/dynamics/orbel.py
new file mode 100644
index 000000000..a4fe5072e
--- /dev/null
+++ b/phoebe/dynamics/orbel.py
@@ -0,0 +1,257 @@
+
+from numpy import pi, sin, cos, tan, arccos, arctan, arctan2, sqrt, dot, cross, array
+
+TINY = 4.0e-15
+TWOPI = 2.0*pi
+
+cape = None
+absr = None
+
+def orbel_el2xv(gm, ialpha, elmts):
+    """
+    Computes cartesian positions and velocities given orbital elements.
+
+    Rewritten from swift/orbel/orbel_el2xv.f.
+
+    Reference: Levison, H.F., Duncan, M.J., The long-term dynamical behavior of short-period comets. Icarus 108, 18-36, 1994.
+
+    """
+    global cape
+
+    a, e, inc, capom, omega, capm = elmts
+
+    if e < 0.0:
+       print('WARNING in orbel_el2xv: e<0, setting e=0!')
+       e = 0.0
+
+    # check for inconsistencies between ialpha and e
+    em1 = e - 1.0
+    if (ialpha == 0 and abs(em1) > TINY) or \
+        (ialpha < 0 and e > 1.0) or \
+        (ialpha > 0 and e < 1.0):
+
+        print('ERROR in orbel_el2xv: ialpha and e inconsistent')
+        print('ialpha = ', ialpha)
+        print('e = ', e)
+        raise ValueError
+
+    # Generate rotation matrices (on p. 42 of Fitzpatrick)
+    sp = sin(omega); cp = cos(omega)
+    so = sin(capom); co = cos(capom)
+    si = sin(inc); ci = cos(inc)
+    d1 = array([ cp*co - sp*so*ci,  cp*so + sp*co*ci, sp*si])
+    d2 = array([-sp*co - cp*so*ci, -sp*so + cp*co*ci, cp*si])
+
+    # Get the other quantities depending on orbit type (i.e., ialpha)
+    if ialpha == -1:
+
+        cape = orbel_ehie(e, capm)
+        scap = sin(cape); ccap = cos(cape)
+        sqe = sqrt(1.0 - e*e)
+        sqgma = sqrt(gm*a)
+        xfac1 = a*(ccap - e)
+        xfac2 = a*sqe*scap
+        ri = 1.0/(a*(1.0 - e*ccap))
+        vfac1 = -ri * sqgma * scap
+        vfac2 = ri * sqgma * sqe * ccap
+
+    else:
+        raise NotImplementedError
+
+    r = d1*xfac1 + d2*xfac2
+    v = d1*vfac1 + d2*vfac2
+
+    return r, v
+
+def get_euler(elmts):
+    """
+    Get corresponding Euler angles.
+
+    Note: Module variable (cape) is used from previous computation in orbel_el2xv().
+
+    """
+    global cape
+
+    a, e, inc, capom, omega, capm = elmts
+
+    theta = 2.0*arctan(sqrt((1.0+e)/(1.0-e))*tan(cape/2.0))
+    
+    euler = array([0.0, 0.0, 0.0])
+    euler[0] = theta + omega
+    euler[1] = capom
+    euler[2] = inc
+
+    cape = None
+
+    return euler
+
+
+def orbel_xv2el(gm, r, v):
+    """
+    Computes orbial elements given cartesian positions and velocities.
+
+    Rewritten from swift/orbel/orbel_xv2el.f.
+
+    Reference: Levison, H.F., Duncan, M.J., The long-term dynamical behavior of short-period comets. Icarus 108, 18-36, 1994.
+
+    """
+    global cape
+    global absr
+
+    # Compute the angular momentum, and thereby the inclination.
+    h = cross(r, v)
+    h2 = dot(h, h)
+    absh = sqrt(h2)
+    inc = arccos(h[2]/absh)
+
+    fac = sqrt(h[0]*h[0] + h[1]*h[1])/absh
+
+    if fac < TINY or inc == 0.0 or inc == pi:
+        capom = 0.0
+        u = arctan2(r[1], r[0])
+        if abs(inc - pi) < 10.0*TINY:
+           u = -u
+    else:
+        capom = arctan2(h[0], -h[1])  
+        u = arctan2(r[2]/sin(inc), r[0]*cos(capom) + r[1]*sin(capom))
+
+    if capom < 0.0:
+        capom += 2.0*pi
+
+    if u < 0.0:
+        u += 2.0*pi
+
+    # Compute the radius R and velocity squared V2, and the dot product RDOTV, the energy per unit mass ENERGY .
+
+    absr = sqrt(dot(r, r))
+    v2 = dot(v, v)
+    absv = sqrt(v2)
+    vdotr = dot(r, v)
+    energy = 0.5*v2 - gm/absr
+
+    # Determine type of conic section and label it via IALPHA
+    if abs(energy*absr/gm) < sqrt(TINY):
+        ialpha = 0
+    else:
+        if energy < 0.0:
+            ialpha = -1
+        else:
+            ialpha = +1
+
+    # Depending on the conic type, determine the remaining elements
+    if ialpha == -1:
+        a = -0.5*gm/energy  
+        fac = 1.0 - h2/(gm*a)
+
+        if fac > TINY:
+            e = sqrt(fac)
+            face = (a-absr)/(a*e)
+
+            if face > 1.0:
+               cape = 0.0
+            else:
+                if face > -1.0:
+                   cape = arccos(face)
+                else:
+                   cape = pi
+
+            if vdotr < 0.0:
+                cape = 2.0*pi - cape
+
+            cw = (cos(cape) - e)/(1.0 - e*cos(cape))
+            sw = sqrt(1.0 - e*e)*sin(cape)/(1.0 - e*cos(cape))
+            w = arctan2(sw, cw)
+
+            if w < 0.0:
+                w += 2.0*pi
+
+        else:
+            e = 0.0
+            w = u
+            cape = u
+
+        capm = cape - e*sin(cape)
+        omega = u - w
+        if omega < 0.0:
+            omega += 2.0*pi
+        omega = omega - int(omega/(2.0*pi))*2.0*pi
+
+    else:
+        raise NotImplementedError
+
+    elmts = array([a, e, inc, capom, omega, capm])
+
+    return elmts
+
+def get_roche(gm, elmts):
+    """
+    Get corresponding Roche parameters.
+
+    Note: Module variable (absr) is used from previous computation in orbel_el2xv().
+
+    Note: The division by rotperiod is done elsewhere...
+
+    """
+    global absr
+
+    a, e, inc, capom, omega, capm = elmts
+
+    n = sqrt(gm/a**3)
+    P = 2.*pi/n
+
+    roche = array([0.0, 0.0])
+    roche[0] = absr/a
+    roche[1] = P
+
+    absr = None
+
+    return roche
+
+
+def orbel_ehie(e, m):
+    """
+    Solves Kepler's equation.
+
+    Rewritten from swift/orbel/orbel_ehie.f.
+
+    Reference: Levison, H.F., Duncan, M.J., The long-term dynamical behavior of short-period comets. Icarus 108, 18-36, 1994.
+
+    """
+
+    # In this section, bring M into the range (0,TWOPI) and if
+    # the result is greater than PI, solve for (TWOPI - M).
+    iflag = 0
+    nper = int(m/TWOPI)
+    m = m - nper*TWOPI
+    if m < 0.0:
+         m = m + TWOPI
+
+    if m > pi:
+        m = TWOPI - m
+        iflag = 1
+
+    # Make a first guess that works well for e near 1.
+    x = (6.0*m)**(1.0/3.0) - m
+
+    # Iteration loop
+    NMAX = 3
+    for niter in range(NMAX):
+        sa = sin(x + m); ca = cos(x + m)
+        esa = e*sa
+        eca = e*ca
+        f = x - esa
+        fp = 1.0 - eca
+        dx = -f/fp
+        dx = -f/(fp + 0.5*dx*esa)
+        dx = -f/(fp + 0.5*dx*(esa+0.3333333333333333*eca*dx))
+        x = x + dx
+
+    cape = m + x
+
+    if iflag == 1:
+        cape = TWOPI - cape
+        m = TWOPI - m
+
+    return cape
+
+
diff --git a/phoebe/dynamics/xyz.py b/phoebe/dynamics/xyz.py
new file mode 100644
index 000000000..5d7d2674f
--- /dev/null
+++ b/phoebe/dynamics/xyz.py
@@ -0,0 +1,327 @@
+"""
+"""
+
+import numpy as np
+from scipy.optimize import newton
+
+from phoebe import u, c
+from phoebe import conf
+from phoebe.dynamics import geometry
+
+try:
+    import rebound
+except ImportError:
+    _is_rebound = False
+else:
+    _is_rebound = True
+
+try:
+    import reboundx
+except ImportError:
+    _is_reboundx = False
+else:
+    _is_reboundx = True
+
+import logging
+logger = logging.getLogger("DYNAMICS.NBODY")
+logger.addHandler(logging.NullHandler())
+
+_skip_filter_checks = {'check_default': False, 'check_visible': False}
+
+_geometry = None
+
+# Note: From phoebe.backend.mesh, which can't be imported (circularly)!
+
+def Rx(x):
+    c = np.cos(x)
+    s = np.sin(x)
+    return np.array([[1., 0., 0.], [0., c, -s], [0., s, c]])
+
+def Rz(x):
+    c = np.cos(x)
+    s = np.sin(x)
+    return np.array([[c, -s, 0.], [s, c, 0.], [0., 0., 1.]])
+
+def spin_in_system(incl, long_an):
+    return np.dot(Rz(long_an), np.dot(Rx(-incl), np.array([0.,0.,1.])))
+
+
+def dynamics_from_bundle(b, times, compute=None, return_roche_euler=False, **kwargs):
+    """
+    Parse parameters in the bundle and call :func:`dynamics`.
+
+    See :func:`dynamics` for more detailed information.
+
+    NOTE: you must either provide compute (a label) OR all relevant options
+    as kwargs (integrator, stepsize, epsilon, ltte, gr, geometry).
+
+    Args:
+        b: (Bundle) the bundle with a set hierarchy
+        times: (list or array) times at which to run the dynamics
+        return_roche_euler: (bool) whether to return Roche parameters and Euler angles
+        geometry: (string) which geometry to use ('hierarchical', 'twopairs')
+
+    Returns:
+        t: (numpy array) all computed times
+        xs, ys, zs: (numpy arrays) time-dependent Cartesian positions [solRad]
+        vxs, vys, vzs: (numpy arrays) time-dependent Cartesian velocities [solRad/day]
+        ds, Fs: (numpy arrays) Roche parameters [1]
+        theta, longan, incl: (numpy arrays) Euler angles [rad]
+
+    """
+    global _geometry
+
+    b.run_delayed_constraints()
+
+    hier = b.hierarchy
+
+    computeps = b.get_compute(compute=compute, force_ps=True, **_skip_filter_checks)
+    stepsize = computeps.get_value(qualifier='stepsize', stepsize=kwargs.get('stepsize', None), **_skip_filter_checks)
+    ltte = computeps.get_value(qualifier='ltte', ltte=kwargs.get('ltte', None), **_skip_filter_checks)
+    gr = computeps.get_value(qualifier='gr', gr=kwargs.get('gr', None), **_skip_filter_checks)
+    integrator = computeps.get_value(qualifier='integrator', integrator=kwargs.get('integrator', None), **_skip_filter_checks)
+    epsilon = computeps.get_value(qualifier='epsilon', epsilon=kwargs.get('epsilon', None), **_skip_filter_checks)
+    _geometry = computeps.get_value(qualifier='geometry', geometry=kwargs.get('geometry', None), **_skip_filter_checks)
+
+    starrefs = hier.get_stars()
+    orbitrefs = hier.get_orbits()
+
+    masses = [b.get_value(qualifier='mass', unit=u.solMass, component=component, context='component', **_skip_filter_checks) * c.G.to('AU3 / (Msun d2)').value for component in starrefs]  # GM
+    smas = [b.get_value(qualifier='sma', unit=u.AU, component=component, context='component', **_skip_filter_checks) for component in orbitrefs]
+    eccs = [b.get_value(qualifier='ecc', component=component, context='component', **_skip_filter_checks) for component in orbitrefs]
+    incls = [b.get_value(qualifier='incl', unit=u.rad, component=component, context='component', **_skip_filter_checks) for component in orbitrefs]
+    per0s = [b.get_value(qualifier='per0', unit=u.rad, component=component, context='component', **_skip_filter_checks) for component in orbitrefs]
+    long_ans = [b.get_value(qualifier='long_an', unit=u.rad, component=component, context='component', **_skip_filter_checks) for component in orbitrefs]
+    mean_anoms = [b.get_value(qualifier='mean_anom', unit=u.rad, component=component, context='component', **_skip_filter_checks) for component in orbitrefs]
+
+    if return_roche_euler:
+        rotperiods = [b.get_value(qualifier='period', unit=u.d, component=component, context='component', **_skip_filter_checks) for component in starrefs]
+    else:
+        rotperiods = None
+
+    vgamma = b.get_value(qualifier='vgamma', context='system', unit=u.AU/u.d, **_skip_filter_checks)
+    t0 = b.get_value(qualifier='t0', context='system', unit=u.d, **_skip_filter_checks)
+
+    j2 = computeps.get_value(qualifier='j2', j2=kwargs.get('j2', None), **_skip_filter_checks)
+    j2s = [b.get_value(qualifier='j2', unit=u.dimensionless_unscaled, component=component, context='component', **_skip_filter_checks) for component in starrefs]
+    requivs = [b.get_value(qualifier='requiv', unit=u.AU, component=component, context='component', **_skip_filter_checks) for component in starrefs]
+    incls_ = [b.get_value(qualifier='incl', unit=u.rad, component=component, context='component', **_skip_filter_checks) for component in starrefs]
+    long_ans_ = [b.get_value(qualifier='long_an', unit=u.rad, component=component, context='component', **_skip_filter_checks) for component in starrefs]
+    
+    nbod = len(masses)
+    spins = []
+    for j in range(0, nbod):
+        spins.append(spin_in_system(incls_[j], long_ans_[j]))
+
+    elmts = []
+    for j in range(0, nbod-1):
+        elmts.append([smas[j], eccs[j], incls[j], long_ans[j], per0s[j], mean_anoms[j]])
+
+    xi, yi, zi, vxi, vyi, vzi = geometry.geometry(masses, elmts, geometry=_geometry)
+
+    return dynamics(times, masses, xi, yi, zi, vxi, vyi, vzi, \
+               rotperiods=rotperiods, \
+               t0=t0, \
+               vgamma=vgamma, \
+               integrator=integrator, \
+               stepsize=stepsize, \
+               epsilon=epsilon, \
+               ltte=ltte, \
+               gr=gr, \
+               j2=j2, \
+               j2s=j2s, \
+               requivs=requivs, \
+               spins=spins, \
+               return_roche_euler=return_roche_euler \
+               )
+
+
+def dynamics(times, masses, xi, yi, zi, vxi, vyi, vzi, \
+        rotperiods=None, \
+        t0=0.0, \
+        vgamma=0.0, \
+        integrator='ias15', \
+        epsilon=1.0e-9, \
+        stepsize=0.01, \
+        ltte=False, \
+        gr=False, \
+        j2=False, \
+        j2s=None, \
+        requivs=None, \
+        spins=None, \
+        return_roche_euler=False \
+        ):
+
+    """
+    Computes N-body dynamics for initial conditions in Cartesian coordinates ("xyz").
+
+    See :func:`dynamics_from_bundle` for a wrapper around this function
+    which automatically handles passing everything in the correct order
+    and in the correct units.
+
+    Note: orbits = stars - 1
+
+    Args:
+        times: (iterable) times when to compute positions and velocities [days]
+        masses: (iterable) mass for each star [solMass]
+        xi, yi, zi: (iterable) initial Cartesian positions for each star [AU]
+        vxi, vyi, vzi: (iterable) initial Cartesian velocities for each star [AU/day]
+        rotperiods: (iterable) rotation periods for each star [day]
+        t0: (float) epoch at which elements were given [days]
+        vgamma: (float) gamma velocity [AU/day]
+        integrator: (string) which integrator from the Rebound package
+        stepsize: (float) step size [day]
+        epsilon: (float) relative error controlling the step size [1]
+        ltte: (bool) whether to account for light travel time effects
+        gr: (bool) whether to account for general relativity effects
+        j2: (bool) whether to account for J2 = -C20 or oblateness
+        j2s: (iterable) J2 = -C20 or oblateness for each star [1]
+        requivs: (iterable) equatorial radius for each star [AU]
+        spins: (iterable) spin axes (x, y, z) for each star [1]
+        return_roche_euler: (bool) whether to return Roche parameters and Euler angles
+
+    Returns:
+        t: (numpy array) all computed times
+        xs, ys, zs: (numpy arrays) time-dependent Cartesian positions [solRad]
+        vxs, vys, vzs: (numpy arrays) time-dependent Cartesian velocities [solRad/day]
+        ds, Fs: (numpy arrays) Roche parameters [1]
+        theta, longan, incl: (numpy arrays) Euler angles [rad]
+
+    """
+
+    global _geometry
+
+    if not _is_rebound:
+        raise ImportError("rebound is not installed")
+
+    if gr and not _is_reboundx:
+        raise ImportError("reboundx is not installed (required for gr effects)")
+
+    if j2 and not _is_reboundx:
+        raise ImportError("reboundx is not installed (required for j2 effects)")
+
+    times = np.asarray(times)
+
+    nbod = len(masses)
+    xs = np.zeros((nbod, len(times)))
+    ys = np.zeros((nbod, len(times)))
+    zs = np.zeros((nbod, len(times)))
+    vxs = np.zeros((nbod, len(times)))
+    vys = np.zeros((nbod, len(times)))
+    vzs = np.zeros((nbod, len(times)))
+
+    if return_roche_euler:
+        ds = np.zeros((nbod, len(times)))
+        Fs = np.zeros((nbod, len(times)))
+        ethetas = np.zeros((nbod, len(times)))
+        elongans = np.zeros((nbod, len(times)))
+        eincls = np.zeros((nbod, len(times)))
+
+    sim = rebound.Simulation()
+
+    sim.integrator = integrator
+    sim.dt = stepsize
+    sim.ri_ias15.epsilon = epsilon
+    sim.ri_whfast.corrector = 17
+    sim.ri_whfast.safe_mode = 0;
+    sim.G = 1.0
+
+    if conf.devel:
+        sim.status()
+
+    for j in range(0, nbod):
+        sim.add(primary=None, m=masses[j], x=xi[j], y=yi[j], z=zi[j], vx=vxi[j], vy=vyi[j], vz=vzi[j])
+
+    sim.move_to_com()
+
+    for particle in sim.particles:
+        particle.vz -= vgamma
+
+    if _is_reboundx:
+        rebx = reboundx.Extras(sim)
+
+    if gr:
+        logger.info("enabling 'gr_full' in reboundx")
+        gr = rebx.load_force("gr_full")
+        gr.params["c"] = c.c.to("AU/d").value
+        rebx.add_force(gr)
+
+    if j2:
+        logger.info("enabling 'gravitational_harmonics' in reboundx")
+        rebx = reboundx.Extras(sim)
+        gh = rebx.load_force("gravitational_harmonics")
+        rebx.add_force(gh)
+
+        for j in range(0, nbod):
+            sim.particles[j].params["J2"] = j2s[j]
+            sim.particles[j].params["R_eq"] = requivs[j]
+            sim.particles[j].params["Omega"] = spins[j]
+
+    rb = np.zeros((nbod, 3))
+    vb = np.zeros((nbod, 3))
+
+    for i,time in enumerate(times):
+
+        sim.integrate(time, exact_finish_time=True)
+
+        for j in range(0, nbod):
+
+            if ltte:
+                particle = _ltte(sim, j, time)
+            else:
+                particle = sim.particles[j]
+
+            rb[j][0] = particle.x
+            rb[j][1] = particle.y
+            rb[j][2] = particle.z
+            vb[j][0] = particle.vx
+            vb[j][1] = particle.vy
+            vb[j][2] = particle.vz
+
+        if return_roche_euler:
+
+            elmts, euler, roche = geometry.invgeometry(masses, rb, vb, geometry=_geometry)
+
+        fac = (1*u.AU).to(u.solRad).value
+
+        xs[:,i] = fac * rb[:,0]
+        ys[:,i] = fac * rb[:,1]
+        zs[:,i] = fac * rb[:,2]
+        vxs[:,i] = fac * vb[:,0]
+        vys[:,i] = fac * vb[:,1]
+        vzs[:,i] = fac * vb[:,2]
+
+        if return_roche_euler:
+            ds[:,i] = roche[:,0]
+            Fs[:,i] = roche[:,1]/rotperiods[:]
+            ethetas[:,i] = euler[:,0]
+            elongans[:,i] = euler[:,1]
+            eincls[:,i] = euler[:,2]
+
+    if return_roche_euler:
+        return times, xs, ys, zs, vxs, vys, vzs, ds, Fs, ethetas, elongans, eincls
+    else:
+        return times, xs, ys, zs, vxs, vys, vzs
+
+
+def _ltte(sim, j, time):
+    """
+    Light-time effect.
+
+    """
+    scale_factor = (u.AU/c.c).to(u.d).value
+
+    def residual(t):
+        if sim.t != t:
+            sim.integrate(t, exact_finish_time=True)
+        z = sim.particles[j].z
+        return t - z*scale_factor - time
+
+    propertime = newton(residual, time)
+
+    if sim.t != propertime:
+        sim.integrate(propertime)
+
+    return sim.particles[j]
+
+
diff --git a/phoebe/frontend/bundle.py b/phoebe/frontend/bundle.py
index 4727ba902..d7a9164dc 100644
--- a/phoebe/frontend/bundle.py
+++ b/phoebe/frontend/bundle.py
@@ -2778,6 +2778,12 @@ def set_hierarchy(self, *args, **kwargs):
                                     constraint=self._default_label('pot_max', context='constraint'))
 
         for component in self.hierarchy.get_orbits():
+            # NOTE: IN ORDER TO DEFAULT_TRIPLE() WORKS
+            children = self.hierarchy.get_stars_of_children_of(component)
+            if len(children) > 2:
+                logger.warning('constraint {} not working for multiple systems'.format(constraint_func))
+                continue
+
             for constraint_func in ['teffratio', 'requivratio', 'requivsumfrac']:
                 logger.debug('re-creating {} constraint for {}'.format(constraint_func, component))
                 if len(self.filter(context='constraint',
@@ -4134,12 +4140,8 @@ def _get_surf_area(comp):
                                      ]+addl_parameters,
                                      True, 'run_compute')
             elif len(self.hierarchy.get_stars()) > 2:
-                if compute_kind not in []:
-                    report.add_item(self,
-                                    "{} (compute='{}') does not support multiple systems".format(compute_kind, compute),
-                                    [self.hierarchy
-                                     ]+addl_parameters,
-                                     True, 'run_compute')
+                if not conf.devel:
+                     raise NotImplementedError("{} (compute='{}') does not support multiple systems.  Enable developer mode to test.".format(compute_kind, compute))
 
             # sample_from and solution checks
             # check if any parameter is in sample_from but is constrained
diff --git a/phoebe/frontend/default_bundles/default_binary.bundle b/phoebe/frontend/default_bundles/default_binary.bundle
index 6682f0f9f..1be196720 100644
--- a/phoebe/frontend/default_bundles/default_binary.bundle
+++ b/phoebe/frontend/default_bundles/default_binary.bundle
@@ -340,6 +340,22 @@ null
 "Class": "FloatParameter"
 },
 {
+"qualifier": "j2",
+"component": "primary",
+"kind": "star",
+"context": "component",
+"description": "J2 = -C20, oblateness, gravitational quadrupole moment",
+"value": 0.0,
+"default_unit": "",
+"limits": [
+null,
+null
+],
+"copy_for": false,
+"advanced": true,
+"Class": "FloatParameter"
+},
+{
 "qualifier": "gravb_bol",
 "component": "primary",
 "kind": "star",
@@ -426,8 +442,8 @@ null
 "description": "Source for bolometric limb darkening coefficients (used only for irradiation; 'auto' to interpolate from the applicable table according to the 'atm' parameter, or the name of a specific atmosphere table)",
 "choices": [
 "auto",
-"phoenix",
-"ck2004"
+"ck2004",
+"phoenix"
 ],
 "value": "auto",
 "visible_if": "ld_mode_bol:lookup",
@@ -688,6 +704,22 @@ null
 "Class": "FloatParameter"
 },
 {
+"qualifier": "j2",
+"component": "secondary",
+"kind": "star",
+"context": "component",
+"description": "J2 = -C20, oblateness, gravitational quadrupole moment",
+"value": 0.0,
+"default_unit": "",
+"limits": [
+null,
+null
+],
+"copy_for": false,
+"advanced": true,
+"Class": "FloatParameter"
+},
+{
 "qualifier": "gravb_bol",
 "component": "secondary",
 "kind": "star",
@@ -774,8 +806,8 @@ null
 "description": "Source for bolometric limb darkening coefficients (used only for irradiation; 'auto' to interpolate from the applicable table according to the 'atm' parameter, or the name of a specific atmosphere table)",
 "choices": [
 "auto",
-"phoenix",
-"ck2004"
+"ck2004",
+"phoenix"
 ],
 "value": "auto",
 "visible_if": "ld_mode_bol:lookup",
@@ -1982,6 +2014,20 @@ null,
 "Class": "ChoiceParameter"
 },
 {
+"qualifier": "boosting_method",
+"compute": "phoebe01",
+"kind": "phoebe",
+"context": "compute",
+"description": "Type of boosting method",
+"choices": [
+"none"
+],
+"value": "none",
+"copy_for": false,
+"advanced": true,
+"Class": "ChoiceParameter"
+},
+{
 "qualifier": "mesh_method",
 "component": "_default",
 "compute": "phoebe01",
@@ -2087,10 +2133,10 @@ null
 "choices": [
 "phoenix",
 "extern_atmx",
+"ck2004",
 "extern_planckint",
-"blackbody",
 "phoenix",
-"ck2004"
+"blackbody"
 ],
 "value": "ck2004",
 "copy_for": {
@@ -2293,10 +2339,10 @@ null
 "choices": [
 "phoenix",
 "extern_atmx",
+"ck2004",
 "extern_planckint",
-"blackbody",
 "phoenix",
-"ck2004"
+"blackbody"
 ],
 "value": "ck2004",
 "copy_for": false,
@@ -2312,10 +2358,10 @@ null
 "choices": [
 "phoenix",
 "extern_atmx",
+"ck2004",
 "extern_planckint",
-"blackbody",
 "phoenix",
-"ck2004"
+"blackbody"
 ],
 "value": "ck2004",
 "copy_for": false,
diff --git a/phoebe/frontend/default_bundles/default_contact_binary.bundle b/phoebe/frontend/default_bundles/default_contact_binary.bundle
index e58fd0f87..80be3bb06 100644
--- a/phoebe/frontend/default_bundles/default_contact_binary.bundle
+++ b/phoebe/frontend/default_bundles/default_contact_binary.bundle
@@ -340,6 +340,22 @@ null
 "Class": "FloatParameter"
 },
 {
+"qualifier": "j2",
+"component": "primary",
+"kind": "star",
+"context": "component",
+"description": "J2 = -C20, oblateness, gravitational quadrupole moment",
+"value": 0.0,
+"default_unit": "",
+"limits": [
+null,
+null
+],
+"copy_for": false,
+"advanced": true,
+"Class": "FloatParameter"
+},
+{
 "qualifier": "gravb_bol",
 "component": "primary",
 "kind": "star",
@@ -426,8 +442,8 @@ null
 "description": "Source for bolometric limb darkening coefficients (used only for irradiation; 'auto' to interpolate from the applicable table according to the 'atm' parameter, or the name of a specific atmosphere table)",
 "choices": [
 "auto",
-"phoenix",
-"ck2004"
+"ck2004",
+"phoenix"
 ],
 "value": "auto",
 "visible_if": "ld_mode_bol:lookup",
@@ -688,6 +704,22 @@ null
 "Class": "FloatParameter"
 },
 {
+"qualifier": "j2",
+"component": "secondary",
+"kind": "star",
+"context": "component",
+"description": "J2 = -C20, oblateness, gravitational quadrupole moment",
+"value": 0.0,
+"default_unit": "",
+"limits": [
+null,
+null
+],
+"copy_for": false,
+"advanced": true,
+"Class": "FloatParameter"
+},
+{
 "qualifier": "gravb_bol",
 "component": "secondary",
 "kind": "star",
@@ -774,8 +806,8 @@ null
 "description": "Source for bolometric limb darkening coefficients (used only for irradiation; 'auto' to interpolate from the applicable table according to the 'atm' parameter, or the name of a specific atmosphere table)",
 "choices": [
 "auto",
-"phoenix",
-"ck2004"
+"ck2004",
+"phoenix"
 ],
 "value": "auto",
 "visible_if": "ld_mode_bol:lookup",
@@ -2173,6 +2205,20 @@ null,
 "Class": "ChoiceParameter"
 },
 {
+"qualifier": "boosting_method",
+"compute": "phoebe01",
+"kind": "phoebe",
+"context": "compute",
+"description": "Type of boosting method",
+"choices": [
+"none"
+],
+"value": "none",
+"copy_for": false,
+"advanced": true,
+"Class": "ChoiceParameter"
+},
+{
 "qualifier": "mesh_method",
 "component": "_default",
 "compute": "phoebe01",
@@ -2278,10 +2324,10 @@ null
 "choices": [
 "phoenix",
 "extern_atmx",
+"ck2004",
 "extern_planckint",
-"blackbody",
 "phoenix",
-"ck2004"
+"blackbody"
 ],
 "value": "ck2004",
 "copy_for": {
@@ -2515,10 +2561,10 @@ null
 "choices": [
 "phoenix",
 "extern_atmx",
+"ck2004",
 "extern_planckint",
-"blackbody",
 "phoenix",
-"ck2004"
+"blackbody"
 ],
 "value": "ck2004",
 "copy_for": false,
@@ -2534,10 +2580,10 @@ null
 "choices": [
 "phoenix",
 "extern_atmx",
+"ck2004",
 "extern_planckint",
-"blackbody",
 "phoenix",
-"ck2004"
+"blackbody"
 ],
 "value": "ck2004",
 "copy_for": false,
diff --git a/phoebe/frontend/default_bundles/default_star.bundle b/phoebe/frontend/default_bundles/default_star.bundle
index a93ac8444..cc042cd63 100644
--- a/phoebe/frontend/default_bundles/default_star.bundle
+++ b/phoebe/frontend/default_bundles/default_star.bundle
@@ -340,6 +340,22 @@ null
 "Class": "FloatParameter"
 },
 {
+"qualifier": "j2",
+"component": "starA",
+"kind": "star",
+"context": "component",
+"description": "J2 = -C20, oblateness, gravitational quadrupole moment",
+"value": 0.0,
+"default_unit": "",
+"limits": [
+null,
+null
+],
+"copy_for": false,
+"advanced": true,
+"Class": "FloatParameter"
+},
+{
 "qualifier": "gravb_bol",
 "component": "starA",
 "kind": "star",
@@ -426,8 +442,8 @@ null
 "description": "Source for bolometric limb darkening coefficients (used only for irradiation; 'auto' to interpolate from the applicable table according to the 'atm' parameter, or the name of a specific atmosphere table)",
 "choices": [
 "auto",
-"phoenix",
-"ck2004"
+"ck2004",
+"phoenix"
 ],
 "value": "auto",
 "visible_if": "ld_mode_bol:lookup",
@@ -760,6 +776,20 @@ null,
 "Class": "ChoiceParameter"
 },
 {
+"qualifier": "boosting_method",
+"compute": "phoebe01",
+"kind": "phoebe",
+"context": "compute",
+"description": "Type of boosting method",
+"choices": [
+"none"
+],
+"value": "none",
+"copy_for": false,
+"advanced": true,
+"Class": "ChoiceParameter"
+},
+{
 "qualifier": "mesh_method",
 "component": "_default",
 "compute": "phoebe01",
@@ -865,10 +895,10 @@ null
 "choices": [
 "phoenix",
 "extern_atmx",
+"ck2004",
 "extern_planckint",
-"blackbody",
 "phoenix",
-"ck2004"
+"blackbody"
 ],
 "value": "ck2004",
 "copy_for": {
@@ -1022,10 +1052,10 @@ null
 "choices": [
 "phoenix",
 "extern_atmx",
+"ck2004",
 "extern_planckint",
-"blackbody",
 "phoenix",
-"ck2004"
+"blackbody"
 ],
 "value": "ck2004",
 "copy_for": false,
diff --git a/phoebe/parameters/component.py b/phoebe/parameters/component.py
index 257208974..2fb6a9d26 100644
--- a/phoebe/parameters/component.py
+++ b/phoebe/parameters/component.py
@@ -202,6 +202,7 @@ def star(component, **kwargs):
 
     params += [FloatParameter(qualifier='incl', latexfmt=r'i_\mathrm{{ {component} }}', visible_if='hierarchy.is_contact_binary:False', value=kwargs.get('incl', 90), default_unit=u.deg, advanced=True, description='Inclination of the stellar rotation axis')]
     params += [FloatParameter(qualifier='long_an', visible_if='hierarchy.is_contact_binary:False', value=kwargs.get('long_an', 0.0), default_unit=u.deg, advanced=True, description='Longitude of the ascending node (ie. equator) of the star')]
+    params += [FloatParameter(qualifier='j2', value=kwargs.get('j2', 0.0), default_unit=u.dimensionless_unscaled, advanced=True, description='J2 = -C20, oblateness, gravitational quadrupole moment')]
 
     # params += [ChoiceParameter(qualifier='gravblaw_bol', value=kwargs.get('gravblaw_bol', 'zeipel'), choices=['zeipel', 'espinosa', 'claret'], description='Gravity brightening law')]
 
diff --git a/phoebe/parameters/compute.py b/phoebe/parameters/compute.py
index 2707cb91b..7c72d4450 100644
--- a/phoebe/parameters/compute.py
+++ b/phoebe/parameters/compute.py
@@ -127,7 +127,10 @@ def phoebe(**kwargs):
     if conf.devel:
         # note: even though bs isn't an option, its manually added as an option in test_dynamics and test_dynamics_grid
         params += [BoolParameter(visible_if='dynamics_method:bs', qualifier='gr', value=kwargs.get('gr', False), description='Whether to account for general relativity effects')]
+        params += [BoolParameter(visible_if='dynamics_method:bs', qualifier='j2', value=kwargs.get('j2', False), description='Whether to account for J2 = -C20 effects')]
         params += [FloatParameter(visible_if='dynamics_method:bs', qualifier='stepsize', value=kwargs.get('stepsize', 0.01), default_unit=None, description='stepsize for the N-body integrator')]         # TODO: improve description (and units??)
+        params += [FloatParameter(visible_if='dynamics_method:bs', qualifier='epsilon', value=kwargs.get('epsilon', 1.0e-9), default_unit=None, description='epsilon for the N-body integrator')]         # TODO: improve description (and units??)
+        params += [ChoiceParameter(visible_if='dynamics_method:bs', qualifier='geometry', value=kwargs.get('geometry', 'hierarchical'), choices=['hierarchical', 'twopairs'], description='geometry for N-body-xyz integrator')]
         params += [ChoiceParameter(visible_if='dynamics_method:bs', qualifier='integrator', value=kwargs.get('integrator', 'ias15'), choices=['ias15', 'whfast', 'sei', 'leapfrog', 'hermes'], description='Which integrator to use within rebound')]
 
 
diff --git a/phoebe/parameters/constraint.py b/phoebe/parameters/constraint.py
index 077d05219..fb8216168 100644
--- a/phoebe/parameters/constraint.py
+++ b/phoebe/parameters/constraint.py
@@ -1553,14 +1553,27 @@ def mass(b, component, solve_for=None, **kwargs):
 
     component_ps = _get_system_ps(b, component)
 
-    sibling = hier.get_sibling_of(component)
-    sibling_ps = _get_system_ps(b, sibling)
-
     parentorbit = hier.get_parent_of(component)
     parentorbit_ps = _get_system_ps(b, parentorbit)
 
-    mass = component_ps.get_parameter(qualifier='mass', **_skip_filter_checks)
-    mass_sibling = sibling_ps.get_parameter(qualifier='mass', **_skip_filter_checks)
+    m1 = component_ps.get_parameter(qualifier='mass', **_skip_filter_checks)
+    masses = []
+    masses.append(m1)
+
+    siblings = hier.get_stars_of_sibling_of(component)
+    siblings = siblings if isinstance(siblings, list) else [siblings]
+
+    for i, sibling in enumerate(siblings):
+        sibling_ps = _get_system_ps(b, sibling)
+        m2 = sibling_ps.get_parameter(qualifier='mass', **_skip_filter_checks)
+        if i==0:
+            msum = m2
+        else:
+            msum += m2
+        masses.append(m2)
+
+    mass = m1
+    mass_sibling = msum
 
     # we need to find the constraint attached to the other component... but we
     # don't know who is constrained, or whether it belongs to the sibling or parent
@@ -1619,7 +1632,7 @@ def mass(b, component, solve_for=None, **kwargs):
     else:
         raise NotImplementedError
 
-    return lhs, rhs, [mass, mass_sibling, period, sma, q], {'component': component}
+    return lhs, rhs, [period, sma, q] + masses, {'component': component}
 
 
     # ecosw_def = FloatParameter(qualifier='ecosw', value=0.0, default_unit=u.dimensionless_unscaled, limits=(-1.0,1.0), description='Eccentricity times cos of argument of periastron')
diff --git a/phoebe/parameters/parameters.py b/phoebe/parameters/parameters.py
index 0cb8e83a4..a623c325e 100644
--- a/phoebe/parameters/parameters.py
+++ b/phoebe/parameters/parameters.py
@@ -208,7 +208,8 @@
                       'mass', 'dpdt', 'per0',
                       'dperdt', 'ecc', 'deccdt', 't0_perpass', 't0_supconj',
                       't0_ref', 'mean_anom', 'q', 'sma', 'asini', 'ecosw', 'esinw',
-                      'teffratio', 'requivratio', 'requivsumfrac'
+                      'teffratio', 'requivratio', 'requivsumfrac',
+                      'j2'
                       ]
 
 # from dataset:
@@ -242,8 +243,8 @@
 
 # from compute:
 _forbidden_labels += ['enabled', 'dynamics_method', 'ltte', 'comments',
-                      'gr', 'stepsize', 'integrator',
-                      'irrad_method', 'mesh_method', 'distortion_method',
+                      'gr', 'stepsize', 'integrator', 'epsilon', 'geometry',
+                      'irrad_method', 'boosting_method', 'mesh_method', 'distortion_method',
                       'ntriangles', 'rv_grav',
                       'mesh_offset', 'mesh_init_phi', 'horizon_method', 'eclipse_method',
                       'atm', 'lc_method', 'rv_method', 'fti_method', 'fti_oversample',
@@ -10759,14 +10760,23 @@ def get_orbits(self):
         -------
         * (list of strings)
         """
-        #~ l = re.findall(r"[\w']+", self.get_value())
-        # now search for indices of orbit and take the next entry from this flat list
-        #~ return [l[i+1] for i,s in enumerate(l) if s=='orbit']
+        # old version, listing orbits for 2+1 systems in incorrect order :-(
+        # l = re.findall(r"[\w']+", self.get_value())
+        # return [l[i+1] for i,s in enumerate(l) if s=='orbit']
+
+        # adding parents of stars, in the correct order
         orbits = []
         for star in self.get_stars():
             parent = self.get_parent_of(star)
             if parent not in orbits and parent!='component' and parent is not None:
                 orbits.append(parent)
+
+        # adding also parents of orbits for 2+2 systems!
+        for orbit in orbits:
+            parent = self.get_parent_of(orbit)
+            if parent not in orbits and parent!='component' and parent is not None:
+                orbits.append(parent)
+
         return orbits
 
     def _compute_meshables(self):
diff --git a/tests/tests/test_geometry/test_geometry.py b/tests/tests/test_geometry/test_geometry.py
new file mode 100755
index 000000000..e25d10210
--- /dev/null
+++ b/tests/tests/test_geometry/test_geometry.py
@@ -0,0 +1,101 @@
+#!/usr/bin/env python3
+
+import numpy as np
+
+from phoebe.dynamics import geometry
+from phoebe import u
+from phoebe import c
+
+TINY = 4.0e-16
+
+def test_geometry(verbose=False):
+
+    day = u.d.to('s')
+    au = u.au.to('m')
+    M_S = u.solMass.to('kg')
+    G = c.G.to('kg^-1 m^3 s^-2').value
+    gms = G*M_S / (au**3 * day**-2)
+
+    m = gms*np.array([1.0, 0.0])
+    elmts = [[1.0, 0.0, 0.0, 0.0, 0.0, 0.0]]
+
+    rb, vb = geometry.geometry_hierarchical(m, elmts)
+
+    if verbose:
+        print("m = ", m)
+        print("rb[0] = ", rb[0])
+        print("rb[1] = ", rb[1])
+        print("vb[0] = ", vb[0])
+        print("vb[1] = ", vb[1])
+
+    assert(abs(rb[1,0] - 1.0) < TINY)
+    assert(abs(rb[1,1] - 0.0) < TINY)
+    assert(abs(rb[1,2] - 0.0) < TINY)
+    assert(abs(vb[1,0] - 0.0) < TINY)
+    assert(abs(vb[1,1] - 1.7202098947281922e-02) < TINY)
+    assert(abs(vb[1,2] - 0.0) < TINY)
+
+    elmts, euler, roche = geometry.invgeometry_hierarchical(m, rb, vb)
+
+    if verbose:
+        print("")
+        print("m = ", m)
+        print("elmts[0] = ", elmts[0])
+        print("euler[0] = ", euler[0])
+        print("euler[1] = ", euler[1])
+        print("roche[0] = ", roche[0])
+        print("roche[1] = ", roche[1])
+
+    assert(abs(elmts[0,0] - 1.0) < TINY)
+    assert(abs(elmts[0,1] - 0.0) < TINY)
+    assert(abs(elmts[0,2] - 0.0) < TINY)
+    assert(abs(elmts[0,3] - 0.0) < TINY)
+    assert(abs(elmts[0,4] - 0.0) < TINY)
+    assert(abs(elmts[0,5] - 0.0) < TINY)
+
+    m = gms*np.array([1.0, 1.0, 0.5, 0.5])
+    elmts = []
+    elmts.append([0.1, 0.0, 0.0, 0.0, 0.0, 0.0])
+    elmts.append([0.2, 0.0, 0.0, 0.0, 0.0, 0.0])
+    elmts.append([1.0, 0.0, 0.0, 0.0, 0.0, 0.0])
+
+    rb, vb = geometry.geometry_twopairs(m, elmts)
+
+    if verbose:
+        print("")
+        print("m = ", m)
+        print("rb[0] = ", rb[0])
+        print("rb[1] = ", rb[1])
+        print("rb[2] = ", rb[2])
+        print("rb[3] = ", rb[3])
+        print("vb[0] = ", vb[0])
+        print("vb[1] = ", vb[1])
+        print("vb[2] = ", vb[2])
+        print("vb[3] = ", vb[3])
+
+    elmts, euler, roche = geometry.invgeometry_twopairs(m, rb, vb)
+
+    if verbose:
+        print("")
+        print("m = ", m)
+        print("elmts[0] = ", elmts[0])
+        print("elmts[1] = ", elmts[1])
+        print("elmts[2] = ", elmts[2])
+        print("euler[0] = ", euler[0])
+        print("euler[1] = ", euler[1])
+        print("euler[2] = ", euler[2])
+        print("euler[3] = ", euler[3])
+        print("roche[0] = ", roche[0])
+        print("roche[1] = ", roche[1])
+        print("roche[2] = ", roche[2])
+        print("roche[3] = ", roche[3])
+
+    assert(abs(elmts[0,0] - 0.1) < TINY)
+    assert(abs(elmts[1,0] - 0.2) < TINY)
+    assert(abs(elmts[2,0] - 1.0) < TINY)
+
+if __name__ == "__main__":
+
+    test_geometry(verbose=True)
+
+