From 471538f358365990f85c6df65d3a392fc23cd4eb Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Sun, 15 Sep 2024 22:33:57 +0200
Subject: [PATCH 01/39] A version working for triple systems!

---
 phoebe/frontend/bundle.py | 60 +++++++++++++++++++++------------------
 1 file changed, 32 insertions(+), 28 deletions(-)

diff --git a/phoebe/frontend/bundle.py b/phoebe/frontend/bundle.py
index 149ed6168..56b9aeb3b 100644
--- a/phoebe/frontend/bundle.py
+++ b/phoebe/frontend/bundle.py
@@ -2777,25 +2777,26 @@ def set_hierarchy(self, *args, **kwargs):
                 self.add_constraint(constraint.potential_contact_max, component,
                                     constraint=self._default_label('pot_max', context='constraint'))
 
-        for component in self.hierarchy.get_orbits():
-            for constraint_func in ['teffratio', 'requivratio', 'requivsumfrac']:
-                logger.debug('re-creating {} constraint for {}'.format(constraint_func, component))
-                if len(self.filter(context='constraint',
-                                   constraint_func=constraint_func,
-                                   component=component,
-                                   **_skip_filter_checks)):
-                    constraint_param = self.get_constraint(constraint_func=constraint_func,
-                                                           component=component,
-                                                           **_skip_filter_checks)
-                    self.remove_constraint(constraint_func=constraint_func,
-                                           component=component,
-                                           **_skip_filter_checks)
-                    self.add_constraint(getattr(constraint, constraint_func), component,
-                                        solve_for=constraint_param.constrained_parameter.uniquetwig,
-                                        constraint=constraint_param.constraint)
-                else:
-                    self.add_constraint(getattr(constraint, constraint_func), component,
-                                        constraint=self._default_label(constraint_func, context='constraint'))
+# NOTE: COMMENTED DUE TO DEFAULT_TRIPLE()
+#        for component in self.hierarchy.get_orbits():
+#            for constraint_func in ['teffratio', 'requivratio', 'requivsumfrac']:
+#                logger.debug('re-creating {} constraint for {}'.format(constraint_func, component))
+#                if len(self.filter(context='constraint',
+#                                   constraint_func=constraint_func,
+#                                   component=component,
+#                                   **_skip_filter_checks)):
+#                    constraint_param = self.get_constraint(constraint_func=constraint_func,
+#                                                           component=component,
+#                                                           **_skip_filter_checks)
+#                    self.remove_constraint(constraint_func=constraint_func,
+#                                           component=component,
+#                                           **_skip_filter_checks)
+#                    self.add_constraint(getattr(constraint, constraint_func), component,
+#                                        solve_for=constraint_param.constrained_parameter.uniquetwig,
+#                                        constraint=constraint_param.constraint)
+#                else:
+#                    self.add_constraint(getattr(constraint, constraint_func), component,
+#                                        constraint=self._default_label(constraint_func, context='constraint'))
 
 
         for component in self.hierarchy.get_stars():
@@ -2834,8 +2835,10 @@ def set_hierarchy(self, *args, **kwargs):
                                         solve_for=constraint_param.constrained_parameter.uniquetwig,
                                         constraint=constraint_param.constraint)
                 else:
-                    self.add_constraint(constraint.mass, component,
-                                        constraint=self._default_label('mass', context='constraint'))
+# NOTE: COMMENTED DUE TO DEFAULT_TRIPLE()
+#                    self.add_constraint(constraint.mass, component,
+#                                        constraint=self._default_label('mass', context='constraint'))
+                    pass
 
 
                 logger.debug('re-creating comp_sma constraint for {}'.format(component))
@@ -4133,13 +4136,14 @@ def _get_surf_area(comp):
                                     [self.hierarchy
                                      ]+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')
+# NOTE: COMMENTED DUE TO DEFAULT_TRIPLE()
+#            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')
 
             # sample_from and solution checks
             # check if any parameter is in sample_from but is constrained

From bf5fd08f1973f8966b408f213a6d856162f21ee5 Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Thu, 19 Sep 2024 14:14:40 +0200
Subject: [PATCH 02/39] Triple systems with two if's.

---
 phoebe/frontend/bundle.py | 68 +++++++++++++++++++++------------------
 1 file changed, 36 insertions(+), 32 deletions(-)

diff --git a/phoebe/frontend/bundle.py b/phoebe/frontend/bundle.py
index 56b9aeb3b..da198287f 100644
--- a/phoebe/frontend/bundle.py
+++ b/phoebe/frontend/bundle.py
@@ -2777,26 +2777,31 @@ def set_hierarchy(self, *args, **kwargs):
                 self.add_constraint(constraint.potential_contact_max, component,
                                     constraint=self._default_label('pot_max', context='constraint'))
 
-# NOTE: COMMENTED DUE TO DEFAULT_TRIPLE()
-#        for component in self.hierarchy.get_orbits():
-#            for constraint_func in ['teffratio', 'requivratio', 'requivsumfrac']:
-#                logger.debug('re-creating {} constraint for {}'.format(constraint_func, component))
-#                if len(self.filter(context='constraint',
-#                                   constraint_func=constraint_func,
-#                                   component=component,
-#                                   **_skip_filter_checks)):
-#                    constraint_param = self.get_constraint(constraint_func=constraint_func,
-#                                                           component=component,
-#                                                           **_skip_filter_checks)
-#                    self.remove_constraint(constraint_func=constraint_func,
-#                                           component=component,
-#                                           **_skip_filter_checks)
-#                    self.add_constraint(getattr(constraint, constraint_func), component,
-#                                        solve_for=constraint_param.constrained_parameter.uniquetwig,
-#                                        constraint=constraint_param.constraint)
-#                else:
-#                    self.add_constraint(getattr(constraint, constraint_func), component,
-#                                        constraint=self._default_label(constraint_func, 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',
+                                   constraint_func=constraint_func,
+                                   component=component,
+                                   **_skip_filter_checks)):
+                    constraint_param = self.get_constraint(constraint_func=constraint_func,
+                                                           component=component,
+                                                           **_skip_filter_checks)
+                    self.remove_constraint(constraint_func=constraint_func,
+                                           component=component,
+                                           **_skip_filter_checks)
+                    self.add_constraint(getattr(constraint, constraint_func), component,
+                                        solve_for=constraint_param.constrained_parameter.uniquetwig,
+                                        constraint=constraint_param.constraint)
+                else:
+                    self.add_constraint(getattr(constraint, constraint_func), component,
+                                        constraint=self._default_label(constraint_func, context='constraint'))
 
 
         for component in self.hierarchy.get_stars():
@@ -2835,10 +2840,14 @@ def set_hierarchy(self, *args, **kwargs):
                                         solve_for=constraint_param.constrained_parameter.uniquetwig,
                                         constraint=constraint_param.constraint)
                 else:
-# NOTE: COMMENTED DUE TO DEFAULT_TRIPLE()
-#                    self.add_constraint(constraint.mass, component,
-#                                        constraint=self._default_label('mass', context='constraint'))
-                    pass
+                    # NOTE: IN ORDER TO DEFAULT_TRIPLE() WORKS
+                    sibling = self.hierarchy.get_sibling_of(component)
+                    kind = self.hierarchy.get_kind_of(sibling)
+                    if kind != 'star':
+                        logger.warning('constraint mass not working for multiple systems')
+                    else:
+                        self.add_constraint(constraint.mass, component,
+                                            constraint=self._default_label('mass', context='constraint'))
 
 
                 logger.debug('re-creating comp_sma constraint for {}'.format(component))
@@ -4136,14 +4145,9 @@ def _get_surf_area(comp):
                                     [self.hierarchy
                                      ]+addl_parameters,
                                      True, 'run_compute')
-# NOTE: COMMENTED DUE TO DEFAULT_TRIPLE()
-#            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')
+            elif len(self.hierarchy.get_stars()) > 2:
+                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

From 54cdb1c9bf20769446303873bf16f7eaa24e5b2b Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Thu, 19 Sep 2024 15:46:43 +0200
Subject: [PATCH 03/39] Correction of nbody.py (pip3 install rebound). Note:
 Remeshing does not work. :-(

---
 phoebe/backend/backends.py | 7 +++++--
 phoebe/dynamics/nbody.py   | 2 +-
 2 files changed, 6 insertions(+), 3 deletions(-)

diff --git a/phoebe/backend/backends.py b/phoebe/backend/backends.py
index 011ff74f7..9906d6a56 100644
--- a/phoebe/backend/backends.py
+++ b/phoebe/backend/backends.py
@@ -1094,10 +1094,13 @@ 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)
+#                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
+                # NOTE: IN ORDER TO DYNAMICS_METHOD='REBOUND' WORKS
+                di = None
+                Fi = None
             else:
                 # then allow d to be determined from orbit and original sma
                 # and F to remain fixed
diff --git a/phoebe/dynamics/nbody.py b/phoebe/dynamics/nbody.py
index ef05de081..49cbb9959 100644
--- a/phoebe/dynamics/nbody.py
+++ b/phoebe/dynamics/nbody.py
@@ -292,7 +292,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

From 460a5e121a57af204cc900f677a87e354828165c Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Thu, 19 Sep 2024 21:12:55 +0200
Subject: [PATCH 04/39] Correction of the Euler angle.

---
 phoebe/dynamics/nbody.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/phoebe/dynamics/nbody.py b/phoebe/dynamics/nbody.py
index 49cbb9959..f3343b5f5 100644
--- a/phoebe/dynamics/nbody.py
+++ b/phoebe/dynamics/nbody.py
@@ -302,9 +302,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

From beb696045326410b733da41cf89456f8bcbb111e Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Sun, 22 Sep 2024 23:34:57 +0200
Subject: [PATCH 05/39] Correction of constraints <- double precision!

---
 phoebe/parameters/parameters.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/phoebe/parameters/parameters.py b/phoebe/parameters/parameters.py
index 7a8917960..14e08825c 100644
--- a/phoebe/parameters/parameters.py
+++ b/phoebe/parameters/parameters.py
@@ -7687,7 +7687,7 @@ def __rmath__(self, other, symbol, mathfunc):
                 else:
                     # assume dimensionless
                     other = float(other)*u.dimensionless_unscaled
-                return ConstraintParameter(self._bundle, "%f %s {%s}" % (_value_for_constraint(other), symbol, self.uniquetwig), default_unit=(getattr(self.quantity, mathfunc)(other).unit))
+                return ConstraintParameter(self._bundle, "%0.30f %s {%s}" % (_value_for_constraint(other), symbol, self.uniquetwig), default_unit=(getattr(self.quantity, mathfunc)(other).unit))
             elif isinstance(other, u.Unit) and mathfunc=='__mul__':
                 return self.quantity*other
             else:
@@ -11817,7 +11817,7 @@ def __math__(self, other, symbol, mathfunc):
             else:
                 # assume dimensionless
                 other = float(other)*u.dimensionless_unscaled
-            return ConstraintParameter(self._bundle, "(%s) %s %f" % (self.expr, symbol, _value_for_constraint(other, self)), default_unit=(getattr(self.result, mathfunc)(other).unit))
+            return ConstraintParameter(self._bundle, "(%s) %s %0.30f" % (self.expr, symbol, _value_for_constraint(other, self)), default_unit=(getattr(self.result, mathfunc)(other).unit))
         elif isinstance(other, str):
             return ConstraintParameter(self._bundle, "(%s) %s %s" % (self.expr, symbol, other), default_unit=(getattr(self.result, mathfunc)(eval(other)).unit))
         elif _is_unit(other) and mathfunc=='__mul__':
@@ -11843,7 +11843,7 @@ def __rmath__(self, other, symbol, mathfunc):
             else:
                 # assume dimensionless
                 other = float(other)*u.dimensionless_unscaled
-            return ConstraintParameter(self._bundle, "%f %s (%s)" % (_value_for_constraint(other, self), symbol, self.expr), default_unit=(getattr(self.result, mathfunc)(other).unit))
+            return ConstraintParameter(self._bundle, "%0.30f %s (%s)" % (_value_for_constraint(other, self), symbol, self.expr), default_unit=(getattr(self.result, mathfunc)(other).unit))
         elif isinstance(other, str):
             return ConstraintParameter(self._bundle, "%s %s (%s)" % (other, symbol, self.expr), default_unit=(getattr(self.result, mathfunc)(eval(other)).unit))
         elif _is_unit(other) and mathfunc=='__mul__':

From 05f7a8c9d80866d31fc82e827896272f9b35442f Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Thu, 26 Sep 2024 20:15:22 +0200
Subject: [PATCH 06/39] Epsilon parameter for integrators.

---
 phoebe/dynamics/nbody.py        | 14 ++++++++++++--
 phoebe/parameters/compute.py    |  1 +
 phoebe/parameters/parameters.py |  4 ++--
 3 files changed, 15 insertions(+), 4 deletions(-)

diff --git a/phoebe/dynamics/nbody.py b/phoebe/dynamics/nbody.py
index f3343b5f5..25011db2d 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")
@@ -175,6 +179,12 @@ def residual(t):
     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
diff --git a/phoebe/parameters/compute.py b/phoebe/parameters/compute.py
index 2707cb91b..3fe5957c5 100644
--- a/phoebe/parameters/compute.py
+++ b/phoebe/parameters/compute.py
@@ -128,6 +128,7 @@ def phoebe(**kwargs):
         # 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 += [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='integrator', value=kwargs.get('integrator', 'ias15'), choices=['ias15', 'whfast', 'sei', 'leapfrog', 'hermes'], description='Which integrator to use within rebound')]
 
 
diff --git a/phoebe/parameters/parameters.py b/phoebe/parameters/parameters.py
index 14e08825c..27ab7e01c 100644
--- a/phoebe/parameters/parameters.py
+++ b/phoebe/parameters/parameters.py
@@ -242,8 +242,8 @@
 
 # from compute:
 _forbidden_labels += ['enabled', 'dynamics_method', 'ltte', 'comments',
-                      'gr', 'stepsize', 'integrator',
-                      'irrad_method', 'mesh_method', 'distortion_method',
+                      'gr', 'stepsize', 'integrator', 'epsilon',
+                      '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',

From 7e535e32ea57993084204ac1ddd637c36b791862 Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Thu, 26 Sep 2024 20:16:08 +0200
Subject: [PATCH 07/39] A re-implementation of keplerian.py for triples!
 Elements in Jacobi coordinates -> barycentric frame. A 1:1 correspondence to
 nbody.py (rebound). As of now, no ltte, Euler angles, dpdt, dperdt, ...

---
 phoebe/dynamics/coord_j2b.py   |  77 +++++
 phoebe/dynamics/keplerian.py   | 509 ++++++++-------------------------
 phoebe/dynamics/orbel_ehie.py  |  68 +++++
 phoebe/dynamics/orbel_el2xv.py |  82 ++++++
 4 files changed, 349 insertions(+), 387 deletions(-)
 create mode 100755 phoebe/dynamics/coord_j2b.py
 mode change 100644 => 100755 phoebe/dynamics/keplerian.py
 create mode 100755 phoebe/dynamics/orbel_ehie.py
 create mode 100755 phoebe/dynamics/orbel_el2xv.py

diff --git a/phoebe/dynamics/coord_j2b.py b/phoebe/dynamics/coord_j2b.py
new file mode 100755
index 000000000..f7d15c0c5
--- /dev/null
+++ b/phoebe/dynamics/coord_j2b.py
@@ -0,0 +1,77 @@
+#!/usr/bin/env python3
+
+import numpy as np
+
+def coord_j2b(mass, rj, vj):
+    """
+    Converts from Jacobi to barycentric coordinateas.
+
+    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 = np.array(nbod*[0.0])
+    rb = np.array((nbod*[[0.0, 0.0, 0.0]]))
+    vb = np.array((nbod*[[0.0, 0.0, 0.0]]))
+
+    # 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 = np.array([0.0, 0.0, 0.0])
+    vtmp = np.array([0.0, 0.0, 0.0])
+
+    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
+
+if __name__ == "__main__":
+
+    mass = np.array([1.0, 1.0, 0.5])
+    rj = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]])
+    vj = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]])
+
+    print("mass = ", mass)
+    print("rj[0] = ", rj[0])
+    print("rj[1] = ", rj[1])
+    print("rj[2] = ", rj[2])
+
+    rb, vb = coord_j2b(mass, rj, vj)
+
+    print("rb[0] = ", rb[0])
+    print("rb[1] = ", rb[1])
+    print("rb[2] = ", rb[2])
+
+    com = np.array([0.0, 0.0, 0.0])
+    for i in range(3):
+        com[i] = np.dot(mass, rb[:,i])/np.sum(mass)
+
+    print("com = ", com)
+
diff --git a/phoebe/dynamics/keplerian.py b/phoebe/dynamics/keplerian.py
old mode 100644
new mode 100755
index e6f6c199b..60142a4bd
--- a/phoebe/dynamics/keplerian.py
+++ b/phoebe/dynamics/keplerian.py
@@ -1,11 +1,14 @@
-
+#!/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_j2b
+from phoebe.dynamics import orbel_el2xv
+from phoebe.dynamics import orbel_ehie
 
 import logging
 logger = logging.getLogger("DYNAMICS.KEPLERIAN")
@@ -47,7 +50,6 @@ 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)
@@ -56,402 +58,135 @@ def dynamics_from_bundle(b, times, compute=None, return_euler=False, **kwargs):
     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)
-
-                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)
+    G = c.G.to('AU3 / (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.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]
+    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]
 
-        # 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)
-
-
-
-
-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):
-    """
-    Compute the positions and velocites of each star in their nested
-    Keplerian orbits at a given list of times.
-
-    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.
-
-    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
-
-    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 = [], [], []
-
-    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.
-
-        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()
-
-            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
+        rotperiods = [b.get_value(qualifier='period', unit=u.d, component=component, context='component', **_skip_filter_checks) for component in starrefs]
     else:
-        return times, \
-            xs, ys, zs, \
-            vxs, vys, vzs
-
-
+        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)
 
-def _true_anomaly(M,ecc,itermax=8):
-    r"""
-    Calculation of true and eccentric anomaly in Kepler orbits.
+    return dynamics(times, masses, smas, eccs, incls, per0s, long_ans, mean_anoms, \
+        rotperiods, t0, vgamma, ltte, \
+        return_euler=return_euler)
 
-    ``M`` is the phase of the star, ``ecc`` is the eccentricity
 
-    See p.39 of Hilditch, 'An Introduction To Close Binary Stars':
+def dynamics(times, masses, smas, eccs, incls, per0s, long_ans, mean_anoms, \
+    rotperiods=None, t0=0.0, vgamma=0.0, ltte=False, \
+    return_euler=False):
 
-    Kepler's equation:
+    nbod = len(masses)
+    rj = np.array(nbod*[[0.0, 0.0, 0.0]])
+    vj = np.array(nbod*[[0.0, 0.0, 0.0]])
 
-    .. math::
+    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)))
 
-        E - e\sin E = \frac{2\pi}{P}(t-T)
+    if return_euler:
+        ethetas  = np.zeros((nbod, len(times)))
+        elongans = np.zeros((nbod, len(times)))
+        eincls   = np.zeros((nbod, len(times)))
+
+    fac = u.au.to('m')/u.solRad.to('m')
+
+    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 = sqrt(msum/a**3)
+            P = 2.0*np.pi/n
+            M = mean_anoms[j-1] + n*(t-t0)
+            elmts = [a, eccs[j-1], incls[j-1], long_ans[j-1], per0s[j-1], M]
+ 
+            rj[j], vj[j] = orbel_el2xv.orbel_el2xv(msum, ialpha, elmts)
+
+        # convert to barycentric frame
+        rb, vb = coord_j2b.coord_j2b(masses, rj, vj)
+
+        rb *= fac
+        vb *= fac
+
+        xs[:,i] = -rb[:,0]
+        ys[:,i] = -rb[:,1]
+        zs[:,i] = rb[:,2]
+        vxs[:,i] = -vb[:,0]
+        vys[:,i] = -vb[:,1]
+        vzs[:,i] = vb[:,2]
 
-    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:
+        if return_euler:
+            ethetas[:,i] = nbod*[0.0]
+            elongans[:,i] = nbod*[0.0]
+            eincls[:,i]  = nbod*[0.0]
 
-    .. math::
+    if return_euler:
+        return times, xs, ys, zs, vxs, vys, vzs, ethetas, elongans, eincls
+    else:
+        return times, xs, ys, zs, vxs, vys, vzs
+
+if __name__ == "__main__":
+
+    # Sun-Earth
+    times = np.array([0.0])
+    masses = np.array([1.0, 0.0])
+    smas = np.array([1.0])
+    eccs = np.array([0.0])
+    incls = np.array([0.0])
+    per0s = np.array([0.0])
+    long_ans = np.array([0.0])
+    mean_anoms = np.array([0.0])
+
+    G = c.G.to('AU3 / (Msun d2)').value
+    masses *= G
+
+    times, xs, ys, zs, vxs, vys, vzs = dynamics(times, masses, smas, eccs, incls, per0s, long_ans, mean_anoms)
+
+    print("xs = ", xs)
+    print("ys = ", ys)
+    print("zs = ", zs)
+    print("vxs = ", vxs)
+    print("vys = ", vys)
+    print("vzs = ", vzs)
+
+    print("v_of_Earth = ", vys[1][0]*(u.au/u.d).to('m s^-1'), " m/s")
+
+    # 3-body problem
+    times = np.array([0.0])
+    masses = np.array([1.0, 1.0, 1.0])
+    smas = np.array([1.0, 10.0])
+    eccs = np.array([0.0, 0.0])
+    incls = np.array([0.0, 0.0])
+    per0s = np.array([0.0, 0.0])
+    long_ans = np.array([0.0, 0.0])
+    mean_anoms = np.array([0.0, 0.0])
+
+    times, xs, ys, zs, vxs, vys, vzs = dynamics(times, masses, smas, eccs, incls, per0s, long_ans, mean_anoms)
+
+    print("")
+    print("xs = ", xs)
+    print("ys = ", ys)
+    print("zs = ", zs)
+    print("vxs = ", vxs)
+    print("vys = ", vys)
+    print("vzs = ", 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/orbel_ehie.py b/phoebe/dynamics/orbel_ehie.py
new file mode 100755
index 000000000..362d32c2e
--- /dev/null
+++ b/phoebe/dynamics/orbel_ehie.py
@@ -0,0 +1,68 @@
+#!/usr/bin/env python3
+
+from numpy import sin, cos, pi
+
+NMAX = 3
+TWOPI = 2.0*pi
+
+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
+    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
+
+if __name__ == "__main__":
+
+    m = 0.25
+    e = 0.6
+
+    print("m = ", m)
+    print("e = ", e)
+
+    cape = orbel_ehie(e, m)
+
+    print("cape = ", cape)
+
+    print("")
+    print("M = E - e*sin(E)")
+    print(m, " = ", cape - e*sin(cape))
+
diff --git a/phoebe/dynamics/orbel_el2xv.py b/phoebe/dynamics/orbel_el2xv.py
new file mode 100755
index 000000000..511a16149
--- /dev/null
+++ b/phoebe/dynamics/orbel_el2xv.py
@@ -0,0 +1,82 @@
+#!/usr/bin/env python3
+
+from numpy import sin, cos, sqrt, array
+
+from phoebe.dynamics import orbel_ehie
+
+TINY = 4.0e-15
+
+def orbel_el2xv(gm, ialpha, elmts):
+    """
+    Computes cartesian positions and velocities given orbital elements.
+
+    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.
+
+    """
+
+    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.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
+
+if __name__ == "__main__":
+
+    day = 86400.
+    au = 1.496e11
+    M_S = 2.e30
+    G = 6.67e-11
+    gm = G*M_S / (au**3 * day**-2)
+
+    elmts = [1.0, 0.1, 0.0, 0.0, 0.0, 0.0]
+
+    print("gm = ", gm)
+    print("elmts = ", elmts)
+
+    r, v = orbel_el2xv(gm, -1, elmts)
+
+    print("r = ", r)
+    print("v = ", v)
+
+

From 93ab3e7c4941eda0ef1464966a4406930b0bb368 Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Fri, 27 Sep 2024 12:33:31 +0200
Subject: [PATCH 08/39] Correction of constraints <- "%.17e" format.

---
 phoebe/parameters/parameters.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/phoebe/parameters/parameters.py b/phoebe/parameters/parameters.py
index 27ab7e01c..242f210e2 100644
--- a/phoebe/parameters/parameters.py
+++ b/phoebe/parameters/parameters.py
@@ -7654,7 +7654,7 @@ def __math__(self, other, symbol, mathfunc):
                 default_unit = getattr(self_quantity, mathfunc)(other_quantity).unit
                 return ConstraintParameter(self._bundle, "{%s} %s {%s}" % (self.uniquetwig, symbol, other.uniquetwig), default_unit=default_unit)
             elif isinstance(other, u.Quantity):
-                return ConstraintParameter(self._bundle, "{%s} %s %0.30f" % (self.uniquetwig, symbol, _value_for_constraint(other)), default_unit=(getattr(self.quantity, mathfunc)(other).unit))
+                return ConstraintParameter(self._bundle, "{%s} %s %0.17e" % (self.uniquetwig, symbol, _value_for_constraint(other)), default_unit=(getattr(self.quantity, mathfunc)(other).unit))
             elif isinstance(other, float) or isinstance(other, int):
                 if symbol in ['+', '-'] and hasattr(self, 'default_unit'):
                     # assume same units as self (NOTE: NOT NECESSARILY SI) if addition or subtraction
@@ -7679,7 +7679,7 @@ def __rmath__(self, other, symbol, mathfunc):
             elif isinstance(other, Parameter):
                 return ConstraintParameter(self._bundle, "{%s} %s {%s}" % (other.uniquetwig, symbol, self.uniquetwig), default_unit=(getattr(self.quantity, mathfunc)(other.quantity).unit))
             elif isinstance(other, u.Quantity):
-                return ConstraintParameter(self._bundle, "%0.30f %s {%s}" % (_value_for_constraint(other), symbol, self.uniquetwig), default_unit=(getattr(self.quantity, mathfunc)(other).unit))
+                return ConstraintParameter(self._bundle, "%0.17e %s {%s}" % (_value_for_constraint(other), symbol, self.uniquetwig), default_unit=(getattr(self.quantity, mathfunc)(other).unit))
             elif isinstance(other, float) or isinstance(other, int):
                 if symbol in ['+', '-'] and hasattr(self, 'default_unit'):
                     # assume same units as self if addition or subtraction
@@ -7687,7 +7687,7 @@ def __rmath__(self, other, symbol, mathfunc):
                 else:
                     # assume dimensionless
                     other = float(other)*u.dimensionless_unscaled
-                return ConstraintParameter(self._bundle, "%0.30f %s {%s}" % (_value_for_constraint(other), symbol, self.uniquetwig), default_unit=(getattr(self.quantity, mathfunc)(other).unit))
+                return ConstraintParameter(self._bundle, "%0.17e %s {%s}" % (_value_for_constraint(other), symbol, self.uniquetwig), default_unit=(getattr(self.quantity, mathfunc)(other).unit))
             elif isinstance(other, u.Unit) and mathfunc=='__mul__':
                 return self.quantity*other
             else:
@@ -11809,7 +11809,7 @@ def __math__(self, other, symbol, mathfunc):
             return ConstraintParameter(self._bundle, "(%s) %s {%s}" % (self.expr, symbol, other.uniquetwig), default_unit=(getattr(self.result, mathfunc)(other.quantity).unit))
         elif isinstance(other, u.Quantity):
             #print "***", other, type(other), isinstance(other, ConstraintParameter)
-            return ConstraintParameter(self._bundle, "(%s) %s %0.30f" % (self.expr, symbol, _value_for_constraint(other, self)), default_unit=(getattr(self.result, mathfunc)(other).unit))
+            return ConstraintParameter(self._bundle, "(%s) %s %0.17e" % (self.expr, symbol, _value_for_constraint(other, self)), default_unit=(getattr(self.result, mathfunc)(other).unit))
         elif isinstance(other, float) or isinstance(other, int):
             if symbol in ['+', '-']:
                 # assume same units as self (NOTE: NOT NECESSARILY SI) if addition or subtraction
@@ -11817,7 +11817,7 @@ def __math__(self, other, symbol, mathfunc):
             else:
                 # assume dimensionless
                 other = float(other)*u.dimensionless_unscaled
-            return ConstraintParameter(self._bundle, "(%s) %s %0.30f" % (self.expr, symbol, _value_for_constraint(other, self)), default_unit=(getattr(self.result, mathfunc)(other).unit))
+            return ConstraintParameter(self._bundle, "(%s) %s %0.17e" % (self.expr, symbol, _value_for_constraint(other, self)), default_unit=(getattr(self.result, mathfunc)(other).unit))
         elif isinstance(other, str):
             return ConstraintParameter(self._bundle, "(%s) %s %s" % (self.expr, symbol, other), default_unit=(getattr(self.result, mathfunc)(eval(other)).unit))
         elif _is_unit(other) and mathfunc=='__mul__':
@@ -11835,7 +11835,7 @@ def __rmath__(self, other, symbol, mathfunc):
             return ConstraintParameter(self._bundle, "{%s} %s (%s)" % (other.uniquetwig, symbol, self.expr), default_unit=(getattr(self.result, mathfunc)(other.quantity).unit))
         elif isinstance(other, u.Quantity):
             #~ print "*** rmath", other, type(other)
-            return ConstraintParameter(self._bundle, "%0.30f %s (%s)" % (_value_for_constraint(other, self), symbol, self.expr), default_unit=(getattr(self.result, mathfunc)(other).unit))
+            return ConstraintParameter(self._bundle, "%0.17e %s (%s)" % (_value_for_constraint(other, self), symbol, self.expr), default_unit=(getattr(self.result, mathfunc)(other).unit))
         elif isinstance(other, float) or isinstance(other, int):
             if symbol in ['+', '-']:
                 # assume same units as self if addition or subtraction
@@ -11843,7 +11843,7 @@ def __rmath__(self, other, symbol, mathfunc):
             else:
                 # assume dimensionless
                 other = float(other)*u.dimensionless_unscaled
-            return ConstraintParameter(self._bundle, "%0.30f %s (%s)" % (_value_for_constraint(other, self), symbol, self.expr), default_unit=(getattr(self.result, mathfunc)(other).unit))
+            return ConstraintParameter(self._bundle, "%0.17e %s (%s)" % (_value_for_constraint(other, self), symbol, self.expr), default_unit=(getattr(self.result, mathfunc)(other).unit))
         elif isinstance(other, str):
             return ConstraintParameter(self._bundle, "%s %s (%s)" % (other, symbol, self.expr), default_unit=(getattr(self.result, mathfunc)(eval(other)).unit))
         elif _is_unit(other) and mathfunc=='__mul__':

From 77e64d0848e208daa72a9cefeb2b3949d1fcb852 Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Sun, 29 Sep 2024 14:03:25 +0200
Subject: [PATCH 09/39] solRad units in Keplerian dynamics.

---
 phoebe/dynamics/keplerian.py | 9 ++-------
 1 file changed, 2 insertions(+), 7 deletions(-)

diff --git a/phoebe/dynamics/keplerian.py b/phoebe/dynamics/keplerian.py
index 60142a4bd..18c84eac2 100755
--- a/phoebe/dynamics/keplerian.py
+++ b/phoebe/dynamics/keplerian.py
@@ -59,9 +59,9 @@ def dynamics_from_bundle(b, times, compute=None, return_euler=False, **kwargs):
     starrefs = hier.get_stars()
     orbitrefs = hier.get_orbits()
 
-    G = c.G.to('AU3 / (Msun d2)').value
+    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.au, component=component, context='component', **_skip_filter_checks) for component in orbitrefs]
+    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]
@@ -103,8 +103,6 @@ def dynamics(times, masses, smas, eccs, incls, per0s, long_ans, mean_anoms, \
         elongans = np.zeros((nbod, len(times)))
         eincls   = np.zeros((nbod, len(times)))
 
-    fac = u.au.to('m')/u.solRad.to('m')
-
     for i, t in enumerate(times):
 
         # compute Jacobi coordinates
@@ -123,9 +121,6 @@ def dynamics(times, masses, smas, eccs, incls, per0s, long_ans, mean_anoms, \
         # convert to barycentric frame
         rb, vb = coord_j2b.coord_j2b(masses, rj, vj)
 
-        rb *= fac
-        vb *= fac
-
         xs[:,i] = -rb[:,0]
         ys[:,i] = -rb[:,1]
         zs[:,i] = rb[:,2]

From a9d0399e4eedada77f0688484d7d61dbe6ea623b Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Sun, 29 Sep 2024 18:42:58 +0200
Subject: [PATCH 10/39] Computation of the Euler angles.

---
 phoebe/dynamics/keplerian.py   | 21 ++++++++++++++++-----
 phoebe/dynamics/orbel_el2xv.py | 25 ++++++++++++++++++++++++-
 2 files changed, 40 insertions(+), 6 deletions(-)

diff --git a/phoebe/dynamics/keplerian.py b/phoebe/dynamics/keplerian.py
index 18c84eac2..28b5ea4a8 100755
--- a/phoebe/dynamics/keplerian.py
+++ b/phoebe/dynamics/keplerian.py
@@ -118,6 +118,22 @@ def dynamics(times, masses, smas, eccs, incls, per0s, long_ans, mean_anoms, \
  
             rj[j], vj[j] = orbel_el2xv.orbel_el2xv(msum, ialpha, elmts)
 
+            if return_euler:
+                euler = orbel_el2xv.get_euler(elmts) 
+
+                ethetas[j][i] = euler[0]
+                elongans[j][i] = euler[1]
+                eincls[j][i] = euler[2]
+
+                # need to copy the primary component
+                # need to add np.pi for secondary component
+                if j==1:
+                    ethetas[0][i] = euler[0]
+                    elongans[0][i] = euler[1]
+                    eincls[0][i] = euler[2]
+
+                    ethetas[j][i] += np.pi
+
         # convert to barycentric frame
         rb, vb = coord_j2b.coord_j2b(masses, rj, vj)
 
@@ -128,11 +144,6 @@ def dynamics(times, masses, smas, eccs, incls, per0s, long_ans, mean_anoms, \
         vys[:,i] = -vb[:,1]
         vzs[:,i] = vb[:,2]
 
-        if return_euler:
-            ethetas[:,i] = nbod*[0.0]
-            elongans[:,i] = nbod*[0.0]
-            eincls[:,i]  = nbod*[0.0]
-
     if return_euler:
         return times, xs, ys, zs, vxs, vys, vzs, ethetas, elongans, eincls
     else:
diff --git a/phoebe/dynamics/orbel_el2xv.py b/phoebe/dynamics/orbel_el2xv.py
index 511a16149..3a85e4bff 100755
--- a/phoebe/dynamics/orbel_el2xv.py
+++ b/phoebe/dynamics/orbel_el2xv.py
@@ -1,11 +1,13 @@
 #!/usr/bin/env python3
 
-from numpy import sin, cos, sqrt, array
+from numpy import sin, cos, sqrt, tan, arctan, array
 
 from phoebe.dynamics import orbel_ehie
 
 TINY = 4.0e-15
 
+cape = None
+
 def orbel_el2xv(gm, ialpha, elmts):
     """
     Computes cartesian positions and velocities given orbital elements.
@@ -15,6 +17,7 @@ def orbel_el2xv(gm, ialpha, elmts):
     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
 
@@ -61,6 +64,26 @@ def orbel_el2xv(gm, ialpha, elmts):
 
     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
+
+    return euler
+
 if __name__ == "__main__":
 
     day = 86400.

From 035ee7a93c1990669ca684d179d460254302a876 Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Sun, 29 Sep 2024 19:54:34 +0200
Subject: [PATCH 11/39] vgamma in Keplerian dynamics.

---
 phoebe/dynamics/keplerian.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/phoebe/dynamics/keplerian.py b/phoebe/dynamics/keplerian.py
index 28b5ea4a8..55b7c1cca 100755
--- a/phoebe/dynamics/keplerian.py
+++ b/phoebe/dynamics/keplerian.py
@@ -126,7 +126,7 @@ def dynamics(times, masses, smas, eccs, incls, per0s, long_ans, mean_anoms, \
                 eincls[j][i] = euler[2]
 
                 # need to copy the primary component
-                # need to add np.pi for secondary component
+                # need to add np.pi for the secondary component
                 if j==1:
                     ethetas[0][i] = euler[0]
                     elongans[0][i] = euler[1]
@@ -137,6 +137,9 @@ def dynamics(times, masses, smas, eccs, incls, per0s, long_ans, mean_anoms, \
         # convert to barycentric frame
         rb, vb = coord_j2b.coord_j2b(masses, rj, vj)
 
+        rb[:,2] -= vgamma*(t-t0)
+        vb[:,2] -= vgamma
+
         xs[:,i] = -rb[:,0]
         ys[:,i] = -rb[:,1]
         zs[:,i] = rb[:,2]

From 39327b18b212b710b614179ea8ca54e8fff54b3d Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Sun, 29 Sep 2024 20:10:31 +0200
Subject: [PATCH 12/39] Correction of vgamma unit.

---
 phoebe/dynamics/keplerian.py | 14 +++-----------
 1 file changed, 3 insertions(+), 11 deletions(-)

diff --git a/phoebe/dynamics/keplerian.py b/phoebe/dynamics/keplerian.py
index 55b7c1cca..41669afef 100755
--- a/phoebe/dynamics/keplerian.py
+++ b/phoebe/dynamics/keplerian.py
@@ -52,9 +52,6 @@ def dynamics_from_bundle(b, times, compute=None, return_euler=False, **kwargs):
 
     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()
@@ -70,21 +67,16 @@ def dynamics_from_bundle(b, times, compute=None, return_euler=False, **kwargs):
     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]
 
-    if return_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)
+    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)
 
     return dynamics(times, masses, smas, eccs, incls, per0s, long_ans, mean_anoms, \
-        rotperiods, t0, vgamma, ltte, \
+        t0, vgamma, ltte, \
         return_euler=return_euler)
 
 
 def dynamics(times, masses, smas, eccs, incls, per0s, long_ans, mean_anoms, \
-    rotperiods=None, t0=0.0, vgamma=0.0, ltte=False, \
+    t0=0.0, vgamma=0.0, ltte=False, \
     return_euler=False):
 
     nbod = len(masses)

From bab6f5886b0318df2a8ffa3c95d5708ed96aa9e3 Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Sun, 29 Sep 2024 23:07:27 +0200
Subject: [PATCH 13/39] dpdt, dperdt in Keplerian dynamics.

---
 phoebe/dynamics/keplerian.py | 20 ++++++++++++++++++--
 1 file changed, 18 insertions(+), 2 deletions(-)

diff --git a/phoebe/dynamics/keplerian.py b/phoebe/dynamics/keplerian.py
index 41669afef..5b3d922da 100755
--- a/phoebe/dynamics/keplerian.py
+++ b/phoebe/dynamics/keplerian.py
@@ -70,12 +70,17 @@ def dynamics_from_bundle(b, times, compute=None, return_euler=False, **kwargs):
     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, masses, smas, eccs, incls, per0s, long_ans, mean_anoms, \
+    dpdts, dperdts, \
     t0=0.0, vgamma=0.0, ltte=False, \
     return_euler=False):
 
@@ -104,9 +109,20 @@ def dynamics(times, masses, smas, eccs, incls, per0s, long_ans, mean_anoms, \
             ialpha = -1
             a = smas[j-1]
             n = sqrt(msum/a**3)
-            P = 2.0*np.pi/n
             M = mean_anoms[j-1] + n*(t-t0)
-            elmts = [a, eccs[j-1], incls[j-1], long_ans[j-1], per0s[j-1], M]
+            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
+
+            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_el2xv.orbel_el2xv(msum, ialpha, elmts)
 

From 0be34a1b30ac07fd1a346a6e9a06abd9cce0a1a1 Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Sun, 29 Sep 2024 23:36:10 +0200
Subject: [PATCH 14/39] Documentation of Keplerian dynamics.

---
 phoebe/dynamics/keplerian.py | 40 +++++++++++++++++++++++++++++++++---
 1 file changed, 37 insertions(+), 3 deletions(-)

diff --git a/phoebe/dynamics/keplerian.py b/phoebe/dynamics/keplerian.py
index 5b3d922da..a65905744 100755
--- a/phoebe/dynamics/keplerian.py
+++ b/phoebe/dynamics/keplerian.py
@@ -1,7 +1,7 @@
 #!/usr/bin/env python3
 
 import numpy as np
-from numpy import pi,sqrt,cos,sin,tan,arctan
+from numpy import pi, sqrt, cos, sin, tan, arctan
 from scipy.optimize import newton
 
 from phoebe import u, c
@@ -83,6 +83,38 @@ 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):
+    """
+    Computes Keplerian dynamics, including binaries, triples, quadruples...
+
+    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]
+        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: 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]
+
+    Note: Euler angles are only returned if return_euler is True.
+
+    """
 
     nbod = len(masses)
     rj = np.array(nbod*[[0.0, 0.0, 0.0]])
@@ -126,6 +158,7 @@ def dynamics(times, masses, smas, eccs, incls, per0s, long_ans, mean_anoms, \
  
             rj[j], vj[j] = orbel_el2xv.orbel_el2xv(msum, ialpha, elmts)
 
+            # compute Euler angles
             if return_euler:
                 euler = orbel_el2xv.get_euler(elmts) 
 
@@ -133,8 +166,8 @@ def dynamics(times, masses, smas, eccs, incls, per0s, long_ans, mean_anoms, \
                 elongans[j][i] = euler[1]
                 eincls[j][i] = euler[2]
 
-                # need to copy the primary component
-                # need to add np.pi for the secondary component
+                # need to copy the primary
+                # need to add np.pi for the secondary
                 if j==1:
                     ethetas[0][i] = euler[0]
                     elongans[0][i] = euler[1]
@@ -145,6 +178,7 @@ def dynamics(times, masses, smas, eccs, incls, per0s, long_ans, mean_anoms, \
         # convert to barycentric frame
         rb, vb = coord_j2b.coord_j2b(masses, rj, vj)
 
+        # gamma velocity
         rb[:,2] -= vgamma*(t-t0)
         vb[:,2] -= vgamma
 

From a7b8655cbd8858ed0dc068fb2cc588237de3bbb4 Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Mon, 30 Sep 2024 18:43:40 +0200
Subject: [PATCH 15/39] ltte in Keplerian dynamics.

---
 phoebe/dynamics/keplerian.py | 136 ++++++++++++++++++++++-------------
 1 file changed, 85 insertions(+), 51 deletions(-)

diff --git a/phoebe/dynamics/keplerian.py b/phoebe/dynamics/keplerian.py
index a65905744..63ab8b8d7 100755
--- a/phoebe/dynamics/keplerian.py
+++ b/phoebe/dynamics/keplerian.py
@@ -116,9 +116,81 @@ def dynamics(times, masses, smas, eccs, incls, per0s, long_ans, mean_anoms, \
 
     """
 
+    # 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 = 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
+
+                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_el2xv.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_el2xv.get_euler(elmts)
+
+                    if j==1:
+                        euler[0,:] = euler[1,:]
+                        euler[1,0] += np.pi
+
+            # convert to barycentric frame
+            rb, vb = coord_j2b.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.array(nbod*[[0.0, 0.0, 0.0]])
     vj = np.array(nbod*[[0.0, 0.0, 0.0]])
+    euler = np.array(nbod*[[0.0, 0.0, 0.0]])
 
     xs = np.zeros((nbod, len(times)))
     ys = np.zeros((nbod, len(times)))
@@ -132,62 +204,24 @@ def dynamics(times, masses, smas, eccs, incls, per0s, long_ans, mean_anoms, \
         elongans = np.zeros((nbod, len(times)))
         eincls   = np.zeros((nbod, len(times)))
 
-    for i, t in enumerate(times):
+    # unfortunately, ltte must be computed for each * (separately)!
+    if ltte:
+        scale_factor = (u.solRad/c.c).to(u.d).value
 
-        # compute Jacobi coordinates
-        msum = masses[0]
-        for j in range(1,nbod):
-            msum += masses[j]
-            ialpha = -1
-            a = smas[j-1]
-            n = 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]
+        for j in range(0,nbod):
 
-            if dpdt != 0.0:
-                P_ = P + dpdt*(t-t0)
-                a *= (P_/P)**2
+            def propertime_residual(t, time):
+                _keplerian([t], component=j)
+                z = zs[j][0]
+                return t - z*scale_factor - time
 
-            if dperdt != 0.0:
-                omega += dperdt*(t-t0)
+            propertimes = [newton(propertime_residual, time, args=(time,)) for time in times]
+            propertimes = np.array(propertimes).ravel()
 
-            elmts = [a, eccs[j-1], incls[j-1], long_ans[j-1], omega, M]
- 
-            rj[j], vj[j] = orbel_el2xv.orbel_el2xv(msum, ialpha, elmts)
+            _keplerian(propertimes, component=j)
 
-            # compute Euler angles
-            if return_euler:
-                euler = orbel_el2xv.get_euler(elmts) 
-
-                ethetas[j][i] = euler[0]
-                elongans[j][i] = euler[1]
-                eincls[j][i] = euler[2]
-
-                # need to copy the primary
-                # need to add np.pi for the secondary
-                if j==1:
-                    ethetas[0][i] = euler[0]
-                    elongans[0][i] = euler[1]
-                    eincls[0][i] = euler[2]
-
-                    ethetas[j][i] += np.pi
-
-        # convert to barycentric frame
-        rb, vb = coord_j2b.coord_j2b(masses, rj, vj)
-
-        # gamma velocity
-        rb[:,2] -= vgamma*(t-t0)
-        vb[:,2] -= vgamma
-
-        xs[:,i] = -rb[:,0]
-        ys[:,i] = -rb[:,1]
-        zs[:,i] = rb[:,2]
-        vxs[:,i] = -vb[:,0]
-        vys[:,i] = -vb[:,1]
-        vzs[:,i] = vb[:,2]
+    else:
+        _keplerian(times)
 
     if return_euler:
         return times, xs, ys, zs, vxs, vys, vzs, ethetas, elongans, eincls

From 4b365b41d6a2209abfcef45e778c98d46d07e152 Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Tue, 1 Oct 2024 08:16:37 +0200
Subject: [PATCH 16/39] Rebuilt default binaries.

---
 .../default_bundles/default_binary.bundle     | 76 +++++++++++--------
 .../default_contact_binary.bundle             | 76 +++++++++++--------
 .../default_bundles/default_star.bundle       | 42 ++++++----
 3 files changed, 118 insertions(+), 76 deletions(-)

diff --git a/phoebe/frontend/default_bundles/default_binary.bundle b/phoebe/frontend/default_bundles/default_binary.bundle
index 80135c198..d51012567 100644
--- a/phoebe/frontend/default_bundles/default_binary.bundle
+++ b/phoebe/frontend/default_bundles/default_binary.bundle
@@ -209,7 +209,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "logg at requiv",
-"value": 4.437551877570185,
+"value": 4.437551868254479,
 "default_unit": "",
 "limits": [
 null,
@@ -259,7 +259,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Rotation frequency (wrt the sky)",
-"value": 6.283185,
+"value": 6.283185307179586,
 "default_unit": "rad / d",
 "limits": [
 0.0,
@@ -460,7 +460,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Mass",
-"value": 0.9988131358058301,
+"value": 0.9988131257959815,
 "default_unit": "solMass",
 "limits": [
 1e-12,
@@ -557,7 +557,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "logg at requiv",
-"value": 4.437551877570185,
+"value": 4.437551868254479,
 "default_unit": "",
 "limits": [
 null,
@@ -607,7 +607,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Rotation frequency (wrt the sky)",
-"value": 6.283185,
+"value": 6.283185307179586,
 "default_unit": "rad / d",
 "limits": [
 0.0,
@@ -808,7 +808,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Mass",
-"value": 0.9988131358058301,
+"value": 0.9988131257959815,
 "default_unit": "solMass",
 "limits": [
 1e-12,
@@ -857,7 +857,7 @@ null
 "kind": "orbit",
 "context": "component",
 "description": "Orbital frequency (sidereal)",
-"value": 6.283185,
+"value": 6.283185307179586,
 "default_unit": "rad / d",
 "limits": [
 null,
@@ -988,7 +988,7 @@ null
 "kind": "orbit",
 "context": "component",
 "description": "Mean anomaly at t0@system",
-"value": 89.99999559997653,
+"value": 90.0,
 "default_unit": "deg",
 "limits": [
 null,
@@ -1245,7 +1245,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "6.283185 / {period@primary@component}",
+"value": "6.28318530717958623e+00 / {period@primary@component}",
 "default_unit": "rad / d",
 "constraint_func": "freq",
 "constraint_kwargs": {
@@ -1261,7 +1261,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "log10((({mass@primary@component} / ({requiv@primary@component} ** 2.000000)) * 2942.206218) * 9.319541)",
+"value": "log10((({mass@primary@component} / ({requiv@primary@component} ** 2.000000)) * 2.94220621750441933e+03) * 9.31954089506172778e+00)",
 "default_unit": "",
 "constraint_func": "logg",
 "constraint_kwargs": {
@@ -1277,7 +1277,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "1.000000 - {irrad_frac_refl_bol@primary@component}",
+"value": "1.00000000000000000e+00 - {irrad_frac_refl_bol@primary@component}",
 "default_unit": "",
 "constraint_func": "irrad_frac",
 "constraint_kwargs": {
@@ -1293,7 +1293,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "6.283185 / {period@secondary@component}",
+"value": "6.28318530717958623e+00 / {period@secondary@component}",
 "default_unit": "rad / d",
 "constraint_func": "freq",
 "constraint_kwargs": {
@@ -1309,7 +1309,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "log10((({mass@secondary@component} / ({requiv@secondary@component} ** 2.000000)) * 2942.206218) * 9.319541)",
+"value": "log10((({mass@secondary@component} / ({requiv@secondary@component} ** 2.000000)) * 2.94220621750441933e+03) * 9.31954089506172778e+00)",
 "default_unit": "",
 "constraint_func": "logg",
 "constraint_kwargs": {
@@ -1325,7 +1325,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "1.000000 - {irrad_frac_refl_bol@secondary@component}",
+"value": "1.00000000000000000e+00 - {irrad_frac_refl_bol@secondary@component}",
 "default_unit": "",
 "constraint_func": "irrad_frac",
 "constraint_kwargs": {
@@ -1389,7 +1389,7 @@ null
 "kind": "orbit",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "{period@binary@component} / ((((-1.000000 * {period@binary@component}) * {dperdt@binary@component}) / 6.283185307179586231995926937088) + 1.000000000000000000000000000000)",
+"value": "{period@binary@component} / ((((-1.00000000000000000e+00 * {period@binary@component}) * {dperdt@binary@component}) / 6.28318530717958623e+00) + 1.00000000000000000e+00)",
 "default_unit": "d",
 "constraint_func": "period_anom",
 "constraint_kwargs": {
@@ -1405,7 +1405,7 @@ null
 "kind": "orbit",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "(6.283185 * ({t0@system} - {t0_perpass@binary@component})) / {period@binary@component}",
+"value": "(6.28318530717958623e+00 * ({t0@system} - {t0_perpass@binary@component})) / {period@binary@component}",
 "default_unit": "deg",
 "constraint_func": "mean_anom",
 "constraint_kwargs": {
@@ -1463,7 +1463,7 @@ null
 "kind": "orbit",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "6.283185 / {period@binary@component}",
+"value": "6.28318530717958623e+00 / {period@binary@component}",
 "default_unit": "rad / d",
 "constraint_func": "freq",
 "constraint_kwargs": {
@@ -1539,7 +1539,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "(39.478418 * ({sma@binary@component} ** 3.000000)) / ((({period@binary@component} ** 2.000000) * ({q@binary@component} + 1.000000)) * 2942.206217504419328179210424423218)",
+"value": "(3.94784176043574320e+01 * ({sma@binary@component} ** 3.000000)) / ((({period@binary@component} ** 2.000000) * ({q@binary@component} + 1.000000)) * 2.94220621750441933e+03)",
 "default_unit": "solMass",
 "constraint_func": "mass",
 "constraint_kwargs": {
@@ -1561,7 +1561,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "{sma@binary@component} / ((1.000000 / {q@binary@component}) + 1.000000)",
+"value": "{sma@binary@component} / ((1.00000000000000000e+00 / {q@binary@component}) + 1.00000000000000000e+00)",
 "default_unit": "solRad",
 "constraint_func": "comp_sma",
 "constraint_kwargs": {
@@ -1577,7 +1577,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "({sma@binary@component} * (sin({incl@binary@component}))) / ((1.000000 / {q@binary@component}) + 1.000000)",
+"value": "({sma@binary@component} * (sin({incl@binary@component}))) / ((1.00000000000000000e+00 / {q@binary@component}) + 1.00000000000000000e+00)",
 "default_unit": "solRad",
 "constraint_func": "comp_asini",
 "constraint_kwargs": {
@@ -1657,7 +1657,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "(39.478418 * ({sma@binary@component} ** 3.000000)) / ((({period@binary@component} ** 2.000000) * ((1.000000 / {q@binary@component}) + 1.000000)) * 2942.206217504419328179210424423218)",
+"value": "(3.94784176043574320e+01 * ({sma@binary@component} ** 3.000000)) / ((({period@binary@component} ** 2.000000) * ((1.00000000000000000e+00 / {q@binary@component}) + 1.00000000000000000e+00)) * 2.94220621750441933e+03)",
 "default_unit": "solMass",
 "constraint_func": "mass",
 "constraint_kwargs": {
@@ -1982,6 +1982,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",
@@ -2085,11 +2099,11 @@ null
 "context": "compute",
 "description": "Atmosphere table",
 "choices": [
-"extern_atmx",
-"extern_planckint",
 "blackbody",
+"extern_atmx",
 "phoenix",
-"ck2004"
+"ck2004",
+"extern_planckint"
 ],
 "value": "ck2004",
 "copy_for": {
@@ -2290,11 +2304,11 @@ null
 "context": "compute",
 "description": "Atmosphere table",
 "choices": [
-"extern_atmx",
-"extern_planckint",
 "blackbody",
+"extern_atmx",
 "phoenix",
-"ck2004"
+"ck2004",
+"extern_planckint"
 ],
 "value": "ck2004",
 "copy_for": false,
@@ -2308,11 +2322,11 @@ null
 "context": "compute",
 "description": "Atmosphere table",
 "choices": [
-"extern_atmx",
-"extern_planckint",
 "blackbody",
+"extern_atmx",
 "phoenix",
-"ck2004"
+"ck2004",
+"extern_planckint"
 ],
 "value": "ck2004",
 "copy_for": false,
@@ -2508,7 +2522,7 @@ null
 "qualifier": "phoebe_version",
 "context": "setting",
 "description": "Version of PHOEBE",
-"value": "2.4.10",
+"value": "2.4.11.dev+feature-interferometry",
 "copy_for": false,
 "readonly": true,
 "advanced": true,
diff --git a/phoebe/frontend/default_bundles/default_contact_binary.bundle b/phoebe/frontend/default_bundles/default_contact_binary.bundle
index df14b5a04..b529e9d08 100644
--- a/phoebe/frontend/default_bundles/default_contact_binary.bundle
+++ b/phoebe/frontend/default_bundles/default_contact_binary.bundle
@@ -209,7 +209,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "logg at requiv",
-"value": 4.089736163094955,
+"value": 4.089736153779247,
 "default_unit": "",
 "limits": [
 null,
@@ -259,7 +259,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Rotation frequency (wrt the sky)",
-"value": 12.56637,
+"value": 12.566370614359172,
 "default_unit": "rad / d",
 "limits": [
 0.0,
@@ -460,7 +460,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Mass",
-"value": 1.0089067994531355,
+"value": 1.0089067893421308,
 "default_unit": "solMass",
 "limits": [
 1e-12,
@@ -557,7 +557,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "logg at requiv",
-"value": 4.089736163094955,
+"value": 4.089736153779248,
 "default_unit": "",
 "limits": [
 null,
@@ -607,7 +607,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Rotation frequency (wrt the sky)",
-"value": 12.56637,
+"value": 12.566370614359172,
 "default_unit": "rad / d",
 "limits": [
 0.0,
@@ -808,7 +808,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Mass",
-"value": 1.0089067994531355,
+"value": 1.0089067893421308,
 "default_unit": "solMass",
 "limits": [
 1e-12,
@@ -857,7 +857,7 @@ null
 "kind": "orbit",
 "context": "component",
 "description": "Orbital frequency (sidereal)",
-"value": 12.56637,
+"value": 12.566370614359172,
 "default_unit": "rad / d",
 "limits": [
 null,
@@ -988,7 +988,7 @@ null
 "kind": "orbit",
 "context": "component",
 "description": "Mean anomaly at t0@system",
-"value": 89.99999559997653,
+"value": 90.0,
 "default_unit": "deg",
 "limits": [
 null,
@@ -1324,7 +1324,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "6.283185 / {period@primary@component}",
+"value": "6.28318530717958623e+00 / {period@primary@component}",
 "default_unit": "rad / d",
 "constraint_func": "freq",
 "constraint_kwargs": {
@@ -1340,7 +1340,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "log10((({mass@primary@component} / ({requiv@primary@component} ** 2.000000)) * 2942.206218) * 9.319541)",
+"value": "log10((({mass@primary@component} / ({requiv@primary@component} ** 2.000000)) * 2.94220621750441933e+03) * 9.31954089506172778e+00)",
 "default_unit": "",
 "constraint_func": "logg",
 "constraint_kwargs": {
@@ -1356,7 +1356,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "1.000000 - {irrad_frac_refl_bol@primary@component}",
+"value": "1.00000000000000000e+00 - {irrad_frac_refl_bol@primary@component}",
 "default_unit": "",
 "constraint_func": "irrad_frac",
 "constraint_kwargs": {
@@ -1372,7 +1372,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "6.283185 / {period@secondary@component}",
+"value": "6.28318530717958623e+00 / {period@secondary@component}",
 "default_unit": "rad / d",
 "constraint_func": "freq",
 "constraint_kwargs": {
@@ -1388,7 +1388,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "log10((({mass@secondary@component} / ({requiv@secondary@component} ** 2.000000)) * 2942.206218) * 9.319541)",
+"value": "log10((({mass@secondary@component} / ({requiv@secondary@component} ** 2.000000)) * 2.94220621750441933e+03) * 9.31954089506172778e+00)",
 "default_unit": "",
 "constraint_func": "logg",
 "constraint_kwargs": {
@@ -1404,7 +1404,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "1.000000 - {irrad_frac_refl_bol@secondary@component}",
+"value": "1.00000000000000000e+00 - {irrad_frac_refl_bol@secondary@component}",
 "default_unit": "",
 "constraint_func": "irrad_frac",
 "constraint_kwargs": {
@@ -1468,7 +1468,7 @@ null
 "kind": "orbit",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "{period@binary@component} / ((((-1.000000 * {period@binary@component}) * {dperdt@binary@component}) / 6.283185307179586231995926937088) + 1.000000000000000000000000000000)",
+"value": "{period@binary@component} / ((((-1.00000000000000000e+00 * {period@binary@component}) * {dperdt@binary@component}) / 6.28318530717958623e+00) + 1.00000000000000000e+00)",
 "default_unit": "d",
 "constraint_func": "period_anom",
 "constraint_kwargs": {
@@ -1484,7 +1484,7 @@ null
 "kind": "orbit",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "(6.283185 * ({t0@system} - {t0_perpass@binary@component})) / {period@binary@component}",
+"value": "(6.28318530717958623e+00 * ({t0@system} - {t0_perpass@binary@component})) / {period@binary@component}",
 "default_unit": "deg",
 "constraint_func": "mean_anom",
 "constraint_kwargs": {
@@ -1542,7 +1542,7 @@ null
 "kind": "orbit",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "6.283185 / {period@binary@component}",
+"value": "6.28318530717958623e+00 / {period@binary@component}",
 "default_unit": "rad / d",
 "constraint_func": "freq",
 "constraint_kwargs": {
@@ -1698,7 +1698,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "(39.478418 * ({sma@binary@component} ** 3.000000)) / ((({period@binary@component} ** 2.000000) * ({q@binary@component} + 1.000000)) * 2942.206217504419328179210424423218)",
+"value": "(3.94784176043574320e+01 * ({sma@binary@component} ** 3.000000)) / ((({period@binary@component} ** 2.000000) * ({q@binary@component} + 1.000000)) * 2.94220621750441933e+03)",
 "default_unit": "solMass",
 "constraint_func": "mass",
 "constraint_kwargs": {
@@ -1720,7 +1720,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "{sma@binary@component} / ((1.000000 / {q@binary@component}) + 1.000000)",
+"value": "{sma@binary@component} / ((1.00000000000000000e+00 / {q@binary@component}) + 1.00000000000000000e+00)",
 "default_unit": "solRad",
 "constraint_func": "comp_sma",
 "constraint_kwargs": {
@@ -1736,7 +1736,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "({sma@binary@component} * (sin({incl@binary@component}))) / ((1.000000 / {q@binary@component}) + 1.000000)",
+"value": "({sma@binary@component} * (sin({incl@binary@component}))) / ((1.00000000000000000e+00 / {q@binary@component}) + 1.00000000000000000e+00)",
 "default_unit": "solRad",
 "constraint_func": "comp_asini",
 "constraint_kwargs": {
@@ -1832,7 +1832,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "(39.478418 * ({sma@binary@component} ** 3.000000)) / ((({period@binary@component} ** 2.000000) * ((1.000000 / {q@binary@component}) + 1.000000)) * 2942.206217504419328179210424423218)",
+"value": "(3.94784176043574320e+01 * ({sma@binary@component} ** 3.000000)) / ((({period@binary@component} ** 2.000000) * ((1.00000000000000000e+00 / {q@binary@component}) + 1.00000000000000000e+00)) * 2.94220621750441933e+03)",
 "default_unit": "solMass",
 "constraint_func": "mass",
 "constraint_kwargs": {
@@ -2173,6 +2173,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",
@@ -2276,11 +2290,11 @@ null
 "context": "compute",
 "description": "Atmosphere table",
 "choices": [
-"extern_atmx",
-"extern_planckint",
 "blackbody",
+"extern_atmx",
 "phoenix",
-"ck2004"
+"ck2004",
+"extern_planckint"
 ],
 "value": "ck2004",
 "copy_for": {
@@ -2512,11 +2526,11 @@ null
 "context": "compute",
 "description": "Atmosphere table",
 "choices": [
-"extern_atmx",
-"extern_planckint",
 "blackbody",
+"extern_atmx",
 "phoenix",
-"ck2004"
+"ck2004",
+"extern_planckint"
 ],
 "value": "ck2004",
 "copy_for": false,
@@ -2530,11 +2544,11 @@ null
 "context": "compute",
 "description": "Atmosphere table",
 "choices": [
-"extern_atmx",
-"extern_planckint",
 "blackbody",
+"extern_atmx",
 "phoenix",
-"ck2004"
+"ck2004",
+"extern_planckint"
 ],
 "value": "ck2004",
 "copy_for": false,
@@ -2740,7 +2754,7 @@ null
 "qualifier": "phoebe_version",
 "context": "setting",
 "description": "Version of PHOEBE",
-"value": "2.4.10",
+"value": "2.4.11.dev+feature-interferometry",
 "copy_for": false,
 "readonly": true,
 "advanced": true,
diff --git a/phoebe/frontend/default_bundles/default_star.bundle b/phoebe/frontend/default_bundles/default_star.bundle
index aefed6400..ed30b2d3b 100644
--- a/phoebe/frontend/default_bundles/default_star.bundle
+++ b/phoebe/frontend/default_bundles/default_star.bundle
@@ -144,7 +144,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Critical (maximum) value of the equivalent radius for the given morphology",
-"value": 3.4292622920441116,
+"value": 3.429265859647371,
 "default_unit": "solRad",
 "limits": [
 0.0,
@@ -209,7 +209,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "logg at requiv",
-"value": 4.438067632266453,
+"value": 4.438067627303133,
 "default_unit": "",
 "limits": [
 null,
@@ -259,7 +259,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Rotation frequency (wrt the sky)",
-"value": 6.283185,
+"value": 6.283185307179586,
 "default_unit": "rad / d",
 "limits": [
 0.0,
@@ -489,7 +489,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "6.283185 / {period@starA@component}",
+"value": "6.28318530717958623e+00 / {period@starA@component}",
 "default_unit": "rad / d",
 "constraint_func": "freq",
 "constraint_kwargs": {
@@ -505,7 +505,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "log10((({mass@starA@component} / ({requiv@starA@component} ** 2.000000)) * 2942.206218) * 9.319541)",
+"value": "log10((({mass@starA@component} / ({requiv@starA@component} ** 2.000000)) * 2.94220621750441933e+03) * 9.31954089506172778e+00)",
 "default_unit": "",
 "constraint_func": "logg",
 "constraint_kwargs": {
@@ -521,7 +521,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "1.000000 - {irrad_frac_refl_bol@starA@component}",
+"value": "1.00000000000000000e+00 - {irrad_frac_refl_bol@starA@component}",
 "default_unit": "",
 "constraint_func": "irrad_frac",
 "constraint_kwargs": {
@@ -537,7 +537,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "0.814886 * (((2942.206218 * {mass@starA@component}) * (({period@starA@component} / 6.283185) ** 2.000000)) ** 0.333333)",
+"value": "8.14885676770049971e-01 * (((2.94220621750441933e+03 * {mass@starA@component}) * (({period@starA@component} / 6.283185) ** 2.00000000000000000e+00)) ** 3.33333333333333315e-01)",
 "default_unit": "solRad",
 "constraint_func": "requiv_single_max",
 "constraint_kwargs": {
@@ -760,6 +760,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",
@@ -863,11 +877,11 @@ null
 "context": "compute",
 "description": "Atmosphere table",
 "choices": [
-"extern_atmx",
-"extern_planckint",
 "blackbody",
+"extern_atmx",
 "phoenix",
-"ck2004"
+"ck2004",
+"extern_planckint"
 ],
 "value": "ck2004",
 "copy_for": {
@@ -1019,11 +1033,11 @@ null
 "context": "compute",
 "description": "Atmosphere table",
 "choices": [
-"extern_atmx",
-"extern_planckint",
 "blackbody",
+"extern_atmx",
 "phoenix",
-"ck2004"
+"ck2004",
+"extern_planckint"
 ],
 "value": "ck2004",
 "copy_for": false,
@@ -1134,7 +1148,7 @@ null
 "qualifier": "phoebe_version",
 "context": "setting",
 "description": "Version of PHOEBE",
-"value": "2.4.10",
+"value": "2.4.11.dev+feature-interferometry",
 "copy_for": false,
 "readonly": true,
 "advanced": true,

From dea000ae4615e3c3484534049867914976d5d8f5 Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Wed, 2 Oct 2024 13:06:51 +0200
Subject: [PATCH 17/39] ASCII art.

---
 phoebe/dynamics/keplerian.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/phoebe/dynamics/keplerian.py b/phoebe/dynamics/keplerian.py
index 63ab8b8d7..289cf63ce 100755
--- a/phoebe/dynamics/keplerian.py
+++ b/phoebe/dynamics/keplerian.py
@@ -86,6 +86,12 @@ def dynamics(times, masses, smas, eccs, incls, per0s, long_ans, mean_anoms, \
     """
     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.

From b3debd7b8a887059c0ae829c1fa067e2875e7868 Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Thu, 3 Oct 2024 20:34:39 +0200
Subject: [PATCH 18/39] Pass roche **kwarg parameters.

---
 phoebe/backend/backends.py | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/phoebe/backend/backends.py b/phoebe/backend/backends.py
index 9906d6a56..165f96b4d 100644
--- a/phoebe/backend/backends.py
+++ b/phoebe/backend/backends.py
@@ -1037,6 +1037,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
@@ -1063,7 +1065,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()))
@@ -1094,13 +1097,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
-                # NOTE: IN ORDER TO DYNAMICS_METHOD='REBOUND' WORKS
-                di = None
-                Fi = None
+                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

From 4b31a86e76a7a10e36bb7bb325d4005a2c89c2dd Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Fri, 4 Oct 2024 08:34:09 +0200
Subject: [PATCH 19/39] Don't reuse module variable (cape).

---
 phoebe/dynamics/orbel_el2xv.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/phoebe/dynamics/orbel_el2xv.py b/phoebe/dynamics/orbel_el2xv.py
index 3a85e4bff..9c04bd710 100755
--- a/phoebe/dynamics/orbel_el2xv.py
+++ b/phoebe/dynamics/orbel_el2xv.py
@@ -82,6 +82,8 @@ def get_euler(elmts):
     euler[1] = capom
     euler[2] = inc
 
+    cape = None
+
     return euler
 
 if __name__ == "__main__":

From 5ba5500faa35a8bc758bae88159b55ce1a22e2ad Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Mon, 7 Oct 2024 09:08:06 +0200
Subject: [PATCH 20/39] Correction of Euler angles for triples.

---
 phoebe/dynamics/keplerian.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/phoebe/dynamics/keplerian.py b/phoebe/dynamics/keplerian.py
index 289cf63ce..2352ddfd0 100755
--- a/phoebe/dynamics/keplerian.py
+++ b/phoebe/dynamics/keplerian.py
@@ -159,7 +159,8 @@ def _keplerian(times, component=None):
 
                     if j==1:
                         euler[0,:] = euler[1,:]
-                        euler[1,0] += np.pi
+                    if j>=1:
+                        euler[j,0] += np.pi
 
             # convert to barycentric frame
             rb, vb = coord_j2b.coord_j2b(masses, rj, vj)

From 1dd03365d5ae1ed48360ff299b1e0aab93b69049 Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Thu, 3 Oct 2024 16:25:10 +0200
Subject: [PATCH 21/39] Introducing a "xyz" dynamics with a geometry layer.

---
 phoebe/backend/backends.py      |   6 ++
 phoebe/dynamics/__init__.py     |   1 +
 phoebe/dynamics/coord_h2j.py    |  61 +++++++++++
 phoebe/dynamics/coord_j2b.py    |   2 +-
 phoebe/dynamics/geometry.py     | 162 +++++++++++++++++++++++++++++
 phoebe/dynamics/xyz.py          | 175 ++++++++++++++++++++++++++++++++
 phoebe/parameters/compute.py    |   1 +
 phoebe/parameters/parameters.py |  17 +++-
 8 files changed, 420 insertions(+), 5 deletions(-)
 create mode 100755 phoebe/dynamics/coord_h2j.py
 create mode 100755 phoebe/dynamics/geometry.py
 create mode 100644 phoebe/dynamics/xyz.py

diff --git a/phoebe/backend/backends.py b/phoebe/backend/backends.py
index 165f96b4d..d325ab68e 100644
--- a/phoebe/backend/backends.py
+++ b/phoebe/backend/backends.py
@@ -959,6 +959,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
@@ -1025,6 +1028,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))
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_h2j.py b/phoebe/dynamics/coord_h2j.py
new file mode 100755
index 000000000..7d6a526f4
--- /dev/null
+++ b/phoebe/dynamics/coord_h2j.py
@@ -0,0 +1,61 @@
+#!/usr/bin/env python3
+
+import numpy as np
+
+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 = np.array(nbod*[0.0])
+    rj = np.array((nbod*[[0.0, 0.0, 0.0]]))
+    vj = np.array((nbod*[[0.0, 0.0, 0.0]]))
+
+    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
+
+if __name__ == "__main__":
+
+    mass = np.array([1.0, 1.0, 0.5])
+    rh = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]])
+    vh = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]])
+
+    print("mass = ", mass)
+    print("rh[0] = ", rh[0])
+    print("rh[1] = ", rh[1])
+    print("rh[2] = ", rh[2])
+
+    rj, vj = coord_h2j(mass, rh, vh)
+
+    print("rj[0] = ", rj[0])
+    print("rj[1] = ", rj[1])
+    print("rj[2] = ", rj[2])
+
+
diff --git a/phoebe/dynamics/coord_j2b.py b/phoebe/dynamics/coord_j2b.py
index f7d15c0c5..591054685 100755
--- a/phoebe/dynamics/coord_j2b.py
+++ b/phoebe/dynamics/coord_j2b.py
@@ -4,7 +4,7 @@
 
 def coord_j2b(mass, rj, vj):
     """
-    Converts from Jacobi to barycentric coordinateas.
+    Converts from Jacobi to barycentric coordinates.
 
     Rewritten from swift/coord/coord_j2b.f.
 
diff --git a/phoebe/dynamics/geometry.py b/phoebe/dynamics/geometry.py
new file mode 100755
index 000000000..c26adf965
--- /dev/null
+++ b/phoebe/dynamics/geometry.py
@@ -0,0 +1,162 @@
+#!/usr/bin/env python3
+
+import numpy as np
+from phoebe.dynamics import coord_h2j
+from phoebe.dynamics import coord_j2b
+from phoebe.dynamics import orbel_el2xv
+
+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.array((nbod*[[0.0, 0.0, 0.0]]))
+    vj = np.array((nbod*[[0.0, 0.0, 0.0]]))
+
+    # compute Jacobi coordinates
+    msum = m[0]
+    for j in range(1, nbod):
+        msum += m[j]
+        ialpha = -1
+
+        rj[j], vj[j] = orbel_el2xv.orbel_el2xv(msum, ialpha, elmts[j-1])
+
+    # convert to barycentric frame
+    rb, vb = coord_j2b.coord_j2b(m, rj, vj)
+
+    return rb, vb
+
+def geometry_twopairs(m, elmts):
+
+    """
+     _          _               \ 
+    / \        / \               |
+    1 2        3 4               5
+    \_/        \_/               |
+                                / 
+    """
+
+    nbod = len(m)
+    rh = np.array((nbod*[[0.0, 0.0, 0.0]]))
+    vh = np.array((nbod*[[0.0, 0.0, 0.0]]))
+    rj = np.array((nbod*[[0.0, 0.0, 0.0]]))
+    vj = np.array((nbod*[[0.0, 0.0, 0.0]]))
+
+    # (1+2) pair, 1-centric coordinates
+    msum = m[0]+m[1]
+    ialpha = -1
+    r2_1, v2_1 = orbel_el2xv.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_el2xv.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_el2xv.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_h2j.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_el2xv.orbel_el2xv(msum, ialpha, elmts[j-1])
+
+    # convert to barycentric frame
+    rb, vb = coord_j2b.coord_j2b(m, rj, vj)
+
+    return rb, vb
+
+if __name__ == "__main__":
+
+    m = [1.0, 0.0]
+    elmts = [[1.0, 0.0, 0.0, 0.0, 0.0, 0.0]]
+
+    rb, vb = geometry_hierarchical(m, elmts)
+
+    print("m = ", m)
+    print("rb[0] = ", rb[0])
+    print("rb[1] = ", rb[1])
+
+    m = [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_twopairs(m, elmts)
+
+    print("rb[0] = ", rb[0])
+    print("rb[1] = ", rb[1])
+    print("rb[2] = ", rb[2])
+    print("rb[3] = ", rb[3])
+
+
diff --git a/phoebe/dynamics/xyz.py b/phoebe/dynamics/xyz.py
new file mode 100644
index 000000000..9420fadf5
--- /dev/null
+++ b/phoebe/dynamics/xyz.py
@@ -0,0 +1,175 @@
+"""
+"""
+
+import numpy as np
+from scipy.optimize import newton
+
+from phoebe import u, c
+from phoebe import conf
+
+import rebound
+from phoebe.dynamics import geometry
+
+import logging
+logger = logging.getLogger("DYNAMICS.NBODY")
+logger.addHandler(logging.NullHandler())
+
+_skip_filter_checks = {'check_default': False, 'check_visible': False}
+
+def dynamics_from_bundle(b, times, compute=None, return_roche_euler=False, **kwargs):
+    """
+    """
+
+    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)
+
+    nbod = len(masses)
+    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, t0, vgamma, stepsize, ltte, gr, \
+                    integrator, return_roche_euler=return_roche_euler, \
+                    epsilon=epsilon)
+
+
+def dynamics(times, masses, xi, yi, zi, vxi, vyi, vzi,
+        rotperiods=None, t0=0.0, vgamma=0.0, stepsize=0.01, ltte=False, gr=False,
+        integrator='ias15', return_roche_euler=False,
+        epsilon=1.0e-9):
+
+    def particle_ltte(sim, j, time):
+        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]
+
+    times = np.asarray(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()
+
+    nbod = len(masses)
+    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
+
+    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)))
+
+    for i,time in enumerate(times):
+
+        sim.integrate(time, exact_finish_time=True)
+
+        for j in range(len(masses)):
+
+            if ltte:
+                particle = particle_ltte(sim, j, time)
+            else:
+                particle = sim.particles[j]
+
+            xs[j][i] = particle.x
+            ys[j][i] = particle.y
+            zs[j][i] = particle.z
+            vxs[j][i] = particle.vx
+            vys[j][i] = particle.vy
+            vzs[j][i] = particle.vz
+
+            if return_roche_euler:
+                if j==0:
+                    particle = sim.particles[j+1]
+                else:
+                    particle = sim.particles[j]
+
+                orbit = particle.orbit()
+
+                ds[j][i] = orbit.d / orbit.a
+                Fs[j][i] = orbit.P / rotperiods[j]
+
+                ethetas[j][i] = orbit.f + orbit.omega
+                elongans[j][i] = orbit.Omega
+                eincls[j][i] = orbit.inc
+
+                if j==1:
+                    ethetas[j][i] += np.pi
+
+    au_to_solrad = (1*u.AU).to(u.solRad).value
+
+    xs *= au_to_solrad
+    ys *= au_to_solrad
+    zs *= au_to_solrad
+    vxs *= au_to_solrad
+    vys *= au_to_solrad
+    vzs *= au_to_solrad
+
+    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
+
+
diff --git a/phoebe/parameters/compute.py b/phoebe/parameters/compute.py
index 3fe5957c5..d6a61d70d 100644
--- a/phoebe/parameters/compute.py
+++ b/phoebe/parameters/compute.py
@@ -129,6 +129,7 @@ def phoebe(**kwargs):
         params += [BoolParameter(visible_if='dynamics_method:bs', qualifier='gr', value=kwargs.get('gr', False), description='Whether to account for general relativity 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/parameters.py b/phoebe/parameters/parameters.py
index 242f210e2..63f722bea 100644
--- a/phoebe/parameters/parameters.py
+++ b/phoebe/parameters/parameters.py
@@ -242,7 +242,7 @@
 
 # from compute:
 _forbidden_labels += ['enabled', 'dynamics_method', 'ltte', 'comments',
-                      'gr', 'stepsize', 'integrator', 'epsilon',
+                      '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',
@@ -10754,14 +10754,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):

From f07858dce5b49219cb70a715bbb4e6d99e3ff01d Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Sun, 6 Oct 2024 10:23:28 +0200
Subject: [PATCH 22/39] The "xyz" dynamics with an inverse geometry layer.
 Euler angles. Roche parameters. Hierarchical systems. Twopairs systems.

---
 phoebe/dynamics/coord_b2j.py   |  60 +++++++++++
 phoebe/dynamics/coord_h2j.py   |   6 +-
 phoebe/dynamics/coord_j2b.py   |   6 +-
 phoebe/dynamics/geometry.py    |  44 +++++---
 phoebe/dynamics/invgeometry.py | 186 +++++++++++++++++++++++++++++++++
 phoebe/dynamics/keplerian.py   |   9 +-
 phoebe/dynamics/orbel_el2xv.py |   9 +-
 phoebe/dynamics/orbel_xv2el.py | 172 ++++++++++++++++++++++++++++++
 phoebe/dynamics/xyz.py         |  91 ++++++++--------
 9 files changed, 508 insertions(+), 75 deletions(-)
 create mode 100755 phoebe/dynamics/coord_b2j.py
 create mode 100755 phoebe/dynamics/invgeometry.py
 create mode 100755 phoebe/dynamics/orbel_xv2el.py

diff --git a/phoebe/dynamics/coord_b2j.py b/phoebe/dynamics/coord_b2j.py
new file mode 100755
index 000000000..4e0387ddb
--- /dev/null
+++ b/phoebe/dynamics/coord_b2j.py
@@ -0,0 +1,60 @@
+#!/usr/bin/env python3
+
+import numpy as np
+
+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 = np.zeros((nbod))
+    rj = np.zeros((nbod, 3))
+    vj = np.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
+
+if __name__ == "__main__":
+
+    mass = np.array([1.0, 1.0, 0.5])
+    rb = np.array([[-0.5, -0.2, 0.0], [0.5, -0.2, 0.0], [0.0, 0.8, 0.0]])
+    vb = np.array([[-0.5, -0.2, 0.0], [0.5, -0.2, 0.0], [0.0, 0.8, 0.0]])
+
+    print("mass = ", mass)
+    print("rb[0] = ", rb[0])
+    print("rb[1] = ", rb[1])
+    print("rb[2] = ", rb[2])
+
+    rj, vj = coord_b2j(mass, rb, vb)
+
+    print("rj[0] = ", rj[0])
+    print("rj[1] = ", rj[1])
+    print("rj[2] = ", rj[2])
+
diff --git a/phoebe/dynamics/coord_h2j.py b/phoebe/dynamics/coord_h2j.py
index 7d6a526f4..a8bba1369 100755
--- a/phoebe/dynamics/coord_h2j.py
+++ b/phoebe/dynamics/coord_h2j.py
@@ -13,9 +13,9 @@ def coord_h2j(mass, rh, vh):
     """
 
     nbod = len(mass)
-    eta = np.array(nbod*[0.0])
-    rj = np.array((nbod*[[0.0, 0.0, 0.0]]))
-    vj = np.array((nbod*[[0.0, 0.0, 0.0]]))
+    eta = np.zeros((nbod))
+    rj = np.zeros((nbod, 3))
+    vj = np.zeros((nbod, 3))
 
     eta[0] = mass[0]
     for i in range(1,nbod):
diff --git a/phoebe/dynamics/coord_j2b.py b/phoebe/dynamics/coord_j2b.py
index 591054685..d2f683b06 100755
--- a/phoebe/dynamics/coord_j2b.py
+++ b/phoebe/dynamics/coord_j2b.py
@@ -13,9 +13,9 @@ def coord_j2b(mass, rj, vj):
     """
 
     nbod = len(mass)
-    eta = np.array(nbod*[0.0])
-    rb = np.array((nbod*[[0.0, 0.0, 0.0]]))
-    vb = np.array((nbod*[[0.0, 0.0, 0.0]]))
+    eta = np.zeros((nbod))
+    rb = np.zeros((nbod, 3))
+    vb = np.zeros((nbod, 3))
 
     # First compute auxiliary variables, then the barycentric positions
     eta[0] = mass[0]
diff --git a/phoebe/dynamics/geometry.py b/phoebe/dynamics/geometry.py
index c26adf965..827ba2cc5 100755
--- a/phoebe/dynamics/geometry.py
+++ b/phoebe/dynamics/geometry.py
@@ -13,11 +13,11 @@ def geometry(m, elmts, geometry='hierarchical'):
 
     if geometry == 'hierarchical':
 
-        _geometry = geometry_hierarchical
+        _geometry = hierarchical
 
     elif geometry == 'twopairs':
 
-        _geometry = geometry_twopairs
+        _geometry = twopairs
 
     else:
         raise NotImplementedError
@@ -47,8 +47,7 @@ def geometry(m, elmts, geometry='hierarchical'):
 
     return xs, ys, zs, vxs, vys, vzs
 
-def geometry_hierarchical(m, elmts):
-
+def hierarchical(m, elmts):
     """
      _          \                |
     / \          |               |
@@ -58,8 +57,8 @@ def geometry_hierarchical(m, elmts):
     """
 
     nbod = len(m)
-    rj = np.array((nbod*[[0.0, 0.0, 0.0]]))
-    vj = np.array((nbod*[[0.0, 0.0, 0.0]]))
+    rj = np.zeros((nbod, 3))
+    vj = np.zeros((nbod, 3))
 
     # compute Jacobi coordinates
     msum = m[0]
@@ -74,8 +73,7 @@ def geometry_hierarchical(m, elmts):
 
     return rb, vb
 
-def geometry_twopairs(m, elmts):
-
+def twopairs(m, elmts):
     """
      _          _               \ 
     / \        / \               |
@@ -85,10 +83,10 @@ def geometry_twopairs(m, elmts):
     """
 
     nbod = len(m)
-    rh = np.array((nbod*[[0.0, 0.0, 0.0]]))
-    vh = np.array((nbod*[[0.0, 0.0, 0.0]]))
-    rj = np.array((nbod*[[0.0, 0.0, 0.0]]))
-    vj = np.array((nbod*[[0.0, 0.0, 0.0]]))
+    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]
@@ -137,26 +135,40 @@ def geometry_twopairs(m, elmts):
 
 if __name__ == "__main__":
 
-    m = [1.0, 0.0]
+    day = 86400.
+    au = 1.496e11
+    M_S = 2.e30
+    G = 6.67e-11
+    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_hierarchical(m, elmts)
+    rb, vb = hierarchical(m, elmts)
 
     print("m = ", m)
     print("rb[0] = ", rb[0])
     print("rb[1] = ", rb[1])
+    print("vb[0] = ", vb[0])
+    print("vb[1] = ", vb[1])
 
-    m = [1.0, 1.0, 0.5, 0.5]
+    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_twopairs(m, elmts)
+    rb, vb = twopairs(m, elmts)
 
+    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])
 
 
diff --git a/phoebe/dynamics/invgeometry.py b/phoebe/dynamics/invgeometry.py
new file mode 100755
index 000000000..34b75a431
--- /dev/null
+++ b/phoebe/dynamics/invgeometry.py
@@ -0,0 +1,186 @@
+#!/usr/bin/env python3
+
+import numpy as np
+from phoebe.dynamics import coord_b2j
+from phoebe.dynamics import orbel_xv2el
+
+def invgeometry(m, rb, vb, geometry='hierarchical'):
+    """
+    Convert barycentric coordinates to elements for various geometries.
+   
+    """
+
+    if geometry == 'hierarchical':
+
+        _invgeometry = hierarchical
+
+    elif geometry == 'twopairs':
+
+        _invgeometry = twopairs
+
+    else:
+        raise NotImplementedError
+
+    # orientation
+    rb[:,0] *= -1.0
+    rb[:,1] *= -1.0
+    vb[:,0] *= -1.0
+    vb[:,1] *= -1.0
+
+    return _invgeometry(m, rb, vb)
+
+def 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_b2j.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_xv2el.orbel_xv2el(msum, rj[j], vj[j])
+        euler[j,:] = orbel_xv2el.get_euler(msum, elmts[j-1])
+        roche[j,:] = orbel_xv2el.get_roche(msum, elmts[j-1])
+
+        if j==1:
+            euler[0,:] = euler[1,:]
+            roche[0,:] = roche[1,:]
+            euler[0,0] += np.pi
+
+    return elmts, euler, roche
+
+
+def 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_xv2el.orbel_xv2el(msum, r2_1, v2_1)
+    euler[0,:] = orbel_xv2el.get_euler(msum, elmts[0])
+    roche[0,:] = orbel_xv2el.get_roche(msum, elmts[0])
+
+    euler[1,:] = euler[0,:]
+    roche[1,:] = roche[0,:]
+    euler[0,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_xv2el.orbel_xv2el(msum, r4_3, v4_3)
+    euler[2,:] = orbel_xv2el.get_euler(msum, elmts[1])
+    roche[2,:] = orbel_xv2el.get_roche(msum, elmts[1])
+    
+    euler[3,:] = euler[2,:]
+    roche[3,:] = roche[2,:]
+    euler[2,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_xv2el.orbel_xv2el(msum, r34_12, v34_12)
+
+    rj, vj = coord_b2j.coord_b2j(m, rb, vb)
+    
+    # compute osculating elements
+    for j in range(4, nbod):
+        msum += m[j]
+
+        elmts[j-1,:] = orbel_xv2el.orbel_xv2el(msum, rj[j], vj[j])
+        euler[j,:] = orbel_xv2el.get_euler(msum, elmts[j-1])
+        roche[j,:] = orbel_xv2el.get_roche(msum, elmts[j-1])
+
+
+    return elmts, euler, roche
+
+if __name__ == "__main__":
+
+    day = 86400.
+    au = 1.496e11
+    M_S = 2.e30
+    G = 6.67e-11
+    gms = G*M_S / (au**3 * day**-2)
+    vk = np.sqrt(gms/1.0)
+
+    m = gms*np.array([1.0, 0.0])
+    rb = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]])
+    vb = np.array([[0.0, 0.0, 0.0], [0.0, vk, 0.0]])
+
+    elmts, euler, roche = hierarchical(m, rb, vb)
+
+    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])
+
+    m = gms*np.array([1.0, 1.0, 0.5, 0.5])
+    rb = np.array([ \
+        [-0.38333333, 0.0, 0.0], \
+        [-0.28333333, 0.0, 0.0], \
+        [ 0.56666667, 0.0, 0.0], \
+        [ 0.76666667, 0.0, 0.0], \
+        ])
+    vb = np.array([ \
+        [0.0, -0.04852087, 0.0], \
+        [0.0,  0.02860663, 0.0], \
+        [0.0,  0.00063236, 0.0], \
+        [0.0,  0.03919611, 0.0], \
+        ])
+
+    elmts, euler, roche = twopairs(m, rb, vb)
+
+    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])
+
+
diff --git a/phoebe/dynamics/keplerian.py b/phoebe/dynamics/keplerian.py
index 2352ddfd0..69d6b21c9 100755
--- a/phoebe/dynamics/keplerian.py
+++ b/phoebe/dynamics/keplerian.py
@@ -1,7 +1,6 @@
 #!/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
@@ -133,7 +132,7 @@ def _keplerian(times, component=None):
                 msum += masses[j]
                 ialpha = -1
                 a = smas[j-1]
-                n = sqrt(msum/a**3)
+                n = np.sqrt(msum/a**3)
                 M = mean_anoms[j-1] + n*(t-t0)
                 P = 2.0*np.pi/n
                 dpdt = dpdts[j-1]
@@ -195,9 +194,9 @@ def _keplerian(times, component=None):
 
     # ...
     nbod = len(masses)
-    rj = np.array(nbod*[[0.0, 0.0, 0.0]])
-    vj = np.array(nbod*[[0.0, 0.0, 0.0]])
-    euler = np.array(nbod*[[0.0, 0.0, 0.0]])
+    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)))
diff --git a/phoebe/dynamics/orbel_el2xv.py b/phoebe/dynamics/orbel_el2xv.py
index 9c04bd710..c18d44415 100755
--- a/phoebe/dynamics/orbel_el2xv.py
+++ b/phoebe/dynamics/orbel_el2xv.py
@@ -12,7 +12,7 @@ def orbel_el2xv(gm, ialpha, elmts):
     """
     Computes cartesian positions and velocities given orbital elements.
 
-    Rewritten from swift/coord/coord_j2b.f.
+    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.
 
@@ -92,16 +92,17 @@ def get_euler(elmts):
     au = 1.496e11
     M_S = 2.e30
     G = 6.67e-11
-    gm = G*M_S / (au**3 * day**-2)
+    gms = G*M_S / (au**3 * day**-2)
 
     elmts = [1.0, 0.1, 0.0, 0.0, 0.0, 0.0]
 
-    print("gm = ", gm)
+    print("gms = ", gms)
     print("elmts = ", elmts)
 
-    r, v = orbel_el2xv(gm, -1, elmts)
+    r, v = orbel_el2xv(gms, -1, elmts)
 
     print("r = ", r)
     print("v = ", v)
+    print("v[1] = %.22e" % v[1])
 
 
diff --git a/phoebe/dynamics/orbel_xv2el.py b/phoebe/dynamics/orbel_xv2el.py
new file mode 100755
index 000000000..24bb21755
--- /dev/null
+++ b/phoebe/dynamics/orbel_xv2el.py
@@ -0,0 +1,172 @@
+#!/usr/bin/env python3
+
+from numpy import pi, sin, cos, tan, arccos, arctan, arctan2, sqrt, dot, cross, array
+
+TINY = 4.0e-15
+
+cape = None
+absr = None
+
+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:
+        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_euler(gm, 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 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
+
+if __name__ == "__main__":
+
+    day = 86400.
+    au = 1.496e11
+    M_S = 2.e30
+    G = 6.67e-11
+    gms = G*M_S / (au**3 * day**-2)
+
+    r = array([0.9, 0.0, 0.0])
+    v = array([0.0, 1.9066428749416830523700e-02, 0.0])
+
+    print("gms = ", gms)
+    print("r = ", r)
+    print("v = ", v)
+
+    elmts = orbel_xv2el(gms, r, v)
+    
+    print("elmts = ", elmts)
+
+
diff --git a/phoebe/dynamics/xyz.py b/phoebe/dynamics/xyz.py
index 9420fadf5..4342d4db1 100644
--- a/phoebe/dynamics/xyz.py
+++ b/phoebe/dynamics/xyz.py
@@ -9,6 +9,7 @@
 
 import rebound
 from phoebe.dynamics import geometry
+from phoebe.dynamics import invgeometry
 
 import logging
 logger = logging.getLogger("DYNAMICS.NBODY")
@@ -16,9 +17,12 @@
 
 _skip_filter_checks = {'check_default': False, 'check_visible': False}
 
+_geometry = None
+
 def dynamics_from_bundle(b, times, compute=None, return_roche_euler=False, **kwargs):
     """
     """
+    global _geometry
 
     b.run_delayed_constraints()
 
@@ -30,7 +34,7 @@ def dynamics_from_bundle(b, times, compute=None, return_roche_euler=False, **kwa
     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)
+    _geometry = computeps.get_value(qualifier='geometry', geometry=kwargs.get('geometry', None), **_skip_filter_checks)
 
     starrefs = hier.get_stars()
     orbitrefs = hier.get_orbits()
@@ -56,7 +60,7 @@ def dynamics_from_bundle(b, times, compute=None, return_roche_euler=False, **kwa
     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_)
+    xi, yi, zi, vxi, vyi, vzi = geometry.geometry(masses, elmts, geometry=_geometry)
 
     return dynamics(times, masses, xi, yi, zi, vxi, vyi, vzi, \
                     rotperiods, t0, vgamma, stepsize, ltte, gr, \
@@ -69,7 +73,10 @@ def dynamics(times, masses, xi, yi, zi, vxi, vyi, vzi,
         integrator='ias15', return_roche_euler=False,
         epsilon=1.0e-9):
 
+    global _geometry
+
     def particle_ltte(sim, j, time):
+
         scale_factor = (u.AU/c.c).to(u.d).value
 
         def residual(t):
@@ -87,6 +94,21 @@ def residual(t):
 
     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
@@ -98,7 +120,6 @@ def residual(t):
     if conf.devel:
         sim.status()
 
-    nbod = len(masses)
     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])
 
@@ -107,19 +128,8 @@ def residual(t):
     for particle in sim.particles:
         particle.vz -= vgamma
 
-    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)))
+    rb = np.zeros((nbod, 3))
+    vb = np.zeros((nbod, 3))
 
     for i,time in enumerate(times):
 
@@ -132,39 +142,32 @@ def residual(t):
             else:
                 particle = sim.particles[j]
 
-            xs[j][i] = particle.x
-            ys[j][i] = particle.y
-            zs[j][i] = particle.z
-            vxs[j][i] = particle.vx
-            vys[j][i] = particle.vy
-            vzs[j][i] = particle.vz
-
-            if return_roche_euler:
-                if j==0:
-                    particle = sim.particles[j+1]
-                else:
-                    particle = sim.particles[j]
-
-                orbit = particle.orbit()
+            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
 
-                ds[j][i] = orbit.d / orbit.a
-                Fs[j][i] = orbit.P / rotperiods[j]
+        if return_roche_euler:
 
-                ethetas[j][i] = orbit.f + orbit.omega
-                elongans[j][i] = orbit.Omega
-                eincls[j][i] = orbit.inc
+            elmts, euler, roche = invgeometry.invgeometry(masses, rb, vb, geometry=_geometry)
 
-                if j==1:
-                    ethetas[j][i] += np.pi
+        fac = (1*u.AU).to(u.solRad).value
 
-    au_to_solrad = (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]
 
-    xs *= au_to_solrad
-    ys *= au_to_solrad
-    zs *= au_to_solrad
-    vxs *= au_to_solrad
-    vys *= au_to_solrad
-    vzs *= au_to_solrad
+        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

From 8322b8692073f95105e045235120982ae6fbb77a Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Sun, 6 Oct 2024 17:48:36 +0200
Subject: [PATCH 23/39] Correction of orientation and Euler angles.

---
 phoebe/dynamics/invgeometry.py | 16 +++++++++++-----
 1 file changed, 11 insertions(+), 5 deletions(-)

diff --git a/phoebe/dynamics/invgeometry.py b/phoebe/dynamics/invgeometry.py
index 34b75a431..b318aa0a4 100755
--- a/phoebe/dynamics/invgeometry.py
+++ b/phoebe/dynamics/invgeometry.py
@@ -21,6 +21,10 @@ def invgeometry(m, rb, vb, geometry='hierarchical'):
     else:
         raise NotImplementedError
 
+    # cf. mutable numpy arrays
+    rb = np.copy(rb)
+    vb = np.copy(vb)
+
     # orientation
     rb[:,0] *= -1.0
     rb[:,1] *= -1.0
@@ -59,11 +63,11 @@ def hierarchical(m, rb, vb):
         if j==1:
             euler[0,:] = euler[1,:]
             roche[0,:] = roche[1,:]
-            euler[0,0] += np.pi
+        if j>=1:
+            euler[j,0] += np.pi
 
     return elmts, euler, roche
 
-
 def twopairs(m, rb, vb):
     """
      _          _               \ 
@@ -89,7 +93,7 @@ def twopairs(m, rb, vb):
 
     euler[1,:] = euler[0,:]
     roche[1,:] = roche[0,:]
-    euler[0,0] += np.pi
+    euler[1,0] += np.pi
 
     # barycenter
     r12 = (m[0]*rb[0] + m[1]*rb[1])/msum
@@ -106,7 +110,7 @@ def twopairs(m, rb, vb):
     
     euler[3,:] = euler[2,:]
     roche[3,:] = roche[2,:]
-    euler[2,0] += np.pi
+    euler[3,0] += np.pi
 
     # barycenter
     r34 = (m[2]*rb[2] + m[3]*rb[3])/msum
@@ -118,9 +122,10 @@ def twopairs(m, rb, vb):
     v34_12 = v34 - v12
     elmts[2,:] = orbel_xv2el.orbel_xv2el(msum, r34_12, v34_12)
 
+    # everything to Jacobian
     rj, vj = coord_b2j.coord_b2j(m, rb, vb)
     
-    # compute osculating elements
+    # other bodies (also Jacobian)
     for j in range(4, nbod):
         msum += m[j]
 
@@ -128,6 +133,7 @@ def twopairs(m, rb, vb):
         euler[j,:] = orbel_xv2el.get_euler(msum, elmts[j-1])
         roche[j,:] = orbel_xv2el.get_roche(msum, elmts[j-1])
 
+        euler[j,0] += np.pi
 
     return elmts, euler, roche
 

From 41f8024ca7e7d860f97db9f4c59666bf4e3f7fde Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Sun, 6 Oct 2024 21:51:03 +0200
Subject: [PATCH 24/39] Reboundx support. If import reboundx doesn't work, use
 workaround export LD_LIBRARY_PATH=~/.local/lib/python3.9/site-packages/

---
 phoebe/dynamics/xyz.py | 30 ++++++++++++++++++++++++++++--
 1 file changed, 28 insertions(+), 2 deletions(-)

diff --git a/phoebe/dynamics/xyz.py b/phoebe/dynamics/xyz.py
index 4342d4db1..e974fb7bf 100644
--- a/phoebe/dynamics/xyz.py
+++ b/phoebe/dynamics/xyz.py
@@ -6,11 +6,23 @@
 
 from phoebe import u, c
 from phoebe import conf
-
-import rebound
 from phoebe.dynamics import geometry
 from phoebe.dynamics import invgeometry
 
+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())
@@ -75,6 +87,12 @@ def dynamics(times, masses, xi, yi, zi, vxi, vyi, vzi,
 
     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)")
+
     def particle_ltte(sim, j, time):
 
         scale_factor = (u.AU/c.c).to(u.d).value
@@ -117,6 +135,14 @@ def residual(t):
     sim.ri_whfast.corrector = 17
     sim.ri_whfast.safe_mode = 0;
     sim.G = 1.0
+
+    if gr:
+        logger.info("enabling 'gr' in reboundx")
+        rebx = reboundx.Extras(sim)
+        gr = rebx.load_force("gr")
+        gr.params["c"] = c.c.to("AU/d").value
+        rebx.add_force(gr)
+
     if conf.devel:
         sim.status()
 

From 7be6a57732403229e3d094aef77913cdcaa92dc4 Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Sun, 6 Oct 2024 22:08:24 +0200
Subject: [PATCH 25/39] Reboundx update.

---
 phoebe/dynamics/nbody.py | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/phoebe/dynamics/nbody.py b/phoebe/dynamics/nbody.py
index 25011db2d..fc584bef3 100644
--- a/phoebe/dynamics/nbody.py
+++ b/phoebe/dynamics/nbody.py
@@ -169,12 +169,14 @@ 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"

From 5a9da3aaa8ec5552a9b8736b85e6f0f970e9ad53 Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Mon, 7 Oct 2024 09:41:25 +0200
Subject: [PATCH 26/39] Single star without ds, Fs.

---
 phoebe/backend/backends.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/phoebe/backend/backends.py b/phoebe/backend/backends.py
index d325ab68e..f9847bf71 100644
--- a/phoebe/backend/backends.py
+++ b/phoebe/backend/backends.py
@@ -1060,6 +1060,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)

From 74b5234ad23ce0ba877299c5d4a077cfa6cf1c86 Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Mon, 7 Oct 2024 10:30:16 +0200
Subject: [PATCH 27/39] The "xyz" dynamics documentation.

---
 phoebe/dynamics/xyz.py | 53 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 53 insertions(+)

diff --git a/phoebe/dynamics/xyz.py b/phoebe/dynamics/xyz.py
index e974fb7bf..385726663 100644
--- a/phoebe/dynamics/xyz.py
+++ b/phoebe/dynamics/xyz.py
@@ -33,6 +33,26 @@
 
 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
 
@@ -85,6 +105,39 @@ def dynamics(times, masses, xi, yi, zi, vxi, vyi, vzi,
         integrator='ias15', return_roche_euler=False,
         epsilon=1.0e-9):
 
+    """
+    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]
+        t0: (float) epoch at which elements were given [days]
+        rotperiods: (iterable) rotation periods for each star [day]
+        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
+        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:

From e9e02a19e5e7eb82dc9a0301f4fb7af07f6e861a Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Mon, 7 Oct 2024 22:34:16 +0200
Subject: [PATCH 28/39] J2 = -C20 or oblateness.

---
 phoebe/dynamics/xyz.py          | 117 ++++++++++++++++++++++----------
 phoebe/parameters/component.py  |   1 +
 phoebe/parameters/compute.py    |   1 +
 phoebe/parameters/parameters.py |   3 +-
 4 files changed, 85 insertions(+), 37 deletions(-)

diff --git a/phoebe/dynamics/xyz.py b/phoebe/dynamics/xyz.py
index 385726663..3c8980e7c 100644
--- a/phoebe/dynamics/xyz.py
+++ b/phoebe/dynamics/xyz.py
@@ -87,6 +87,10 @@ def dynamics_from_bundle(b, times, compute=None, return_roche_euler=False, **kwa
     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]
+
     nbod = len(masses)
     elmts = []
     for j in range(0, nbod-1):
@@ -95,15 +99,35 @@ def dynamics_from_bundle(b, times, compute=None, return_roche_euler=False, **kwa
     xi, yi, zi, vxi, vyi, vzi = geometry.geometry(masses, elmts, geometry=_geometry)
 
     return dynamics(times, masses, xi, yi, zi, vxi, vyi, vzi, \
-                    rotperiods, t0, vgamma, stepsize, ltte, gr, \
-                    integrator, return_roche_euler=return_roche_euler, \
-                    epsilon=epsilon)
-
-
-def dynamics(times, masses, xi, yi, zi, vxi, vyi, vzi,
-        rotperiods=None, t0=0.0, vgamma=0.0, stepsize=0.01, ltte=False, gr=False,
-        integrator='ias15', return_roche_euler=False,
-        epsilon=1.0e-9):
+               rotperiods=rotperiods, \
+               t0=t0, \
+               vgamma=vgamma, \
+               integrator=integrator, \
+               stepsize=stepsize, \
+               epsilon=epsilon, \
+               ltte=ltte, \
+               gr=gr, \
+               j2=j2, \
+               j2s=j2s, \
+               requivs=j2s, \
+               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, \
+        return_roche_euler=False \
+        ):
 
     """
     Computes N-body dynamics for initial conditions in Cartesian coordinates ("xyz").
@@ -119,14 +143,17 @@ def dynamics(times, masses, xi, yi, zi, vxi, vyi, vzi,
         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]
-        t0: (float) epoch at which elements were given [days]
         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]
         return_roche_euler: (bool) whether to return Roche parameters and Euler angles
 
     Returns:
@@ -146,22 +173,8 @@ def dynamics(times, masses, xi, yi, zi, vxi, vyi, vzi,
     if gr and not _is_reboundx:
         raise ImportError("reboundx is not installed (required for gr effects)")
 
-    def particle_ltte(sim, j, time):
-
-        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]
+    if j2 and not _is_reboundx:
+        raise ImportError("reboundx is not installed (required for j2 effects)")
 
     times = np.asarray(times)
 
@@ -189,13 +202,6 @@ def residual(t):
     sim.ri_whfast.safe_mode = 0;
     sim.G = 1.0
 
-    if gr:
-        logger.info("enabling 'gr' in reboundx")
-        rebx = reboundx.Extras(sim)
-        gr = rebx.load_force("gr")
-        gr.params["c"] = c.c.to("AU/d").value
-        rebx.add_force(gr)
-
     if conf.devel:
         sim.status()
 
@@ -207,6 +213,25 @@ def residual(t):
     for particle in sim.particles:
         particle.vz -= vgamma
 
+    if _is_reboundx:
+        rebx = reboundx.Extras(sim)
+
+    if gr:
+        logger.info("enabling 'gr' in reboundx")
+        gr = rebx.load_force("gr")
+        gr.params["c"] = c.c.to("AU/d").value
+        rebx.add_force(gr)
+
+    if j2:
+        logger.info("enabling 'j2' 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]
+
     rb = np.zeros((nbod, 3))
     vb = np.zeros((nbod, 3))
 
@@ -214,10 +239,10 @@ def residual(t):
 
         sim.integrate(time, exact_finish_time=True)
 
-        for j in range(len(masses)):
+        for j in range(0, nbod):
 
             if ltte:
-                particle = particle_ltte(sim, j, time)
+                particle = _ltte(sim, j, time)
             else:
                 particle = sim.particles[j]
 
@@ -250,8 +275,28 @@ def residual(t):
 
     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/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 d6a61d70d..7c72d4450 100644
--- a/phoebe/parameters/compute.py
+++ b/phoebe/parameters/compute.py
@@ -127,6 +127,7 @@ 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')]
diff --git a/phoebe/parameters/parameters.py b/phoebe/parameters/parameters.py
index 63f722bea..216b68538 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:

From a8f894d6143d1fa66228d0545517d61bc5fa7cbd Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Mon, 7 Oct 2024 22:55:03 +0200
Subject: [PATCH 29/39] Rebuilt default binaries.

---
 .../default_bundles/default_binary.bundle     | 58 ++++++++++++++-----
 .../default_contact_binary.bundle             | 58 ++++++++++++++-----
 .../default_bundles/default_star.bundle       | 32 +++++++---
 3 files changed, 114 insertions(+), 34 deletions(-)

diff --git a/phoebe/frontend/default_bundles/default_binary.bundle b/phoebe/frontend/default_bundles/default_binary.bundle
index d51012567..52718cc3b 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",
@@ -2099,11 +2131,11 @@ null
 "context": "compute",
 "description": "Atmosphere table",
 "choices": [
-"blackbody",
 "extern_atmx",
-"phoenix",
 "ck2004",
-"extern_planckint"
+"extern_planckint",
+"phoenix",
+"blackbody"
 ],
 "value": "ck2004",
 "copy_for": {
@@ -2304,11 +2336,11 @@ null
 "context": "compute",
 "description": "Atmosphere table",
 "choices": [
-"blackbody",
 "extern_atmx",
-"phoenix",
 "ck2004",
-"extern_planckint"
+"extern_planckint",
+"phoenix",
+"blackbody"
 ],
 "value": "ck2004",
 "copy_for": false,
@@ -2322,11 +2354,11 @@ null
 "context": "compute",
 "description": "Atmosphere table",
 "choices": [
-"blackbody",
 "extern_atmx",
-"phoenix",
 "ck2004",
-"extern_planckint"
+"extern_planckint",
+"phoenix",
+"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 b529e9d08..06dfbb081 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",
@@ -2290,11 +2322,11 @@ null
 "context": "compute",
 "description": "Atmosphere table",
 "choices": [
-"blackbody",
 "extern_atmx",
-"phoenix",
 "ck2004",
-"extern_planckint"
+"extern_planckint",
+"phoenix",
+"blackbody"
 ],
 "value": "ck2004",
 "copy_for": {
@@ -2526,11 +2558,11 @@ null
 "context": "compute",
 "description": "Atmosphere table",
 "choices": [
-"blackbody",
 "extern_atmx",
-"phoenix",
 "ck2004",
-"extern_planckint"
+"extern_planckint",
+"phoenix",
+"blackbody"
 ],
 "value": "ck2004",
 "copy_for": false,
@@ -2544,11 +2576,11 @@ null
 "context": "compute",
 "description": "Atmosphere table",
 "choices": [
-"blackbody",
 "extern_atmx",
-"phoenix",
 "ck2004",
-"extern_planckint"
+"extern_planckint",
+"phoenix",
+"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 ed30b2d3b..a7bd67edc 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",
@@ -877,11 +893,11 @@ null
 "context": "compute",
 "description": "Atmosphere table",
 "choices": [
-"blackbody",
 "extern_atmx",
-"phoenix",
 "ck2004",
-"extern_planckint"
+"extern_planckint",
+"phoenix",
+"blackbody"
 ],
 "value": "ck2004",
 "copy_for": {
@@ -1033,11 +1049,11 @@ null
 "context": "compute",
 "description": "Atmosphere table",
 "choices": [
-"blackbody",
 "extern_atmx",
-"phoenix",
 "ck2004",
-"extern_planckint"
+"extern_planckint",
+"phoenix",
+"blackbody"
 ],
 "value": "ck2004",
 "copy_for": false,

From b4b607cd3bebcf42db53ef62da297e64922e0e46 Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Tue, 8 Oct 2024 13:01:50 +0200
Subject: [PATCH 30/39] Correction of J2 (requivs).

---
 phoebe/dynamics/xyz.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/phoebe/dynamics/xyz.py b/phoebe/dynamics/xyz.py
index 3c8980e7c..a457db3fc 100644
--- a/phoebe/dynamics/xyz.py
+++ b/phoebe/dynamics/xyz.py
@@ -109,7 +109,7 @@ def dynamics_from_bundle(b, times, compute=None, return_roche_euler=False, **kwa
                gr=gr, \
                j2=j2, \
                j2s=j2s, \
-               requivs=j2s, \
+               requivs=requivs, \
                return_roche_euler=return_roche_euler \
                )
 

From d91912739e8701ffd08410ac407f71e5614a152d Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Sun, 13 Oct 2024 18:52:03 +0200
Subject: [PATCH 31/39] Use 'gr_full' in reboundx.

---
 phoebe/dynamics/xyz.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/phoebe/dynamics/xyz.py b/phoebe/dynamics/xyz.py
index a457db3fc..d9e5c8257 100644
--- a/phoebe/dynamics/xyz.py
+++ b/phoebe/dynamics/xyz.py
@@ -217,8 +217,8 @@ def dynamics(times, masses, xi, yi, zi, vxi, vyi, vzi, \
         rebx = reboundx.Extras(sim)
 
     if gr:
-        logger.info("enabling 'gr' in reboundx")
-        gr = rebx.load_force("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)
 

From b5832f03550d91965c5913b29c794763ee944ab9 Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Mon, 14 Oct 2024 09:09:41 +0200
Subject: [PATCH 32/39] Correction of J2 = -C20 in reboundx. Using user-defined
 force 'j2' instead of 'gravitational harmonics', where the z-axis was assumed
 to be fixed, not perpendicular to the respective orbit(s). For the moment,
 'j2' must be added as follows:

https://sirrah.troja.mff.cuni.cz/~mira/tmp/phoebe2/reboundx/
---
 phoebe/dynamics/xyz.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/phoebe/dynamics/xyz.py b/phoebe/dynamics/xyz.py
index d9e5c8257..45d343e63 100644
--- a/phoebe/dynamics/xyz.py
+++ b/phoebe/dynamics/xyz.py
@@ -225,7 +225,7 @@ def dynamics(times, masses, xi, yi, zi, vxi, vyi, vzi, \
     if j2:
         logger.info("enabling 'j2' in reboundx")
         rebx = reboundx.Extras(sim)
-        gh = rebx.load_force("gravitational_harmonics")
+        gh = rebx.load_force("j2")
         rebx.add_force(gh)
 
         for j in range(0, nbod):

From d41abe65e652cb9fbf4936ff74725199c9b9d6b8 Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Mon, 14 Oct 2024 18:23:35 +0200
Subject: [PATCH 33/39] Consolidated to coord.py, orbel.py, geometry.py.

---
 phoebe/dynamics/coord.py                     | 134 +++++++++++++
 phoebe/dynamics/coord_b2j.py                 |  60 ------
 phoebe/dynamics/coord_h2j.py                 |  61 ------
 phoebe/dynamics/coord_j2b.py                 |  77 -------
 phoebe/dynamics/geometry.py                  | 201 ++++++++++++++-----
 phoebe/dynamics/invgeometry.py               | 192 ------------------
 phoebe/dynamics/keplerian.py                 |  57 +-----
 phoebe/dynamics/{orbel_xv2el.py => orbel.py} | 159 +++++++++++----
 phoebe/dynamics/orbel_ehie.py                |  68 -------
 phoebe/dynamics/orbel_el2xv.py               | 108 ----------
 10 files changed, 410 insertions(+), 707 deletions(-)
 create mode 100644 phoebe/dynamics/coord.py
 delete mode 100755 phoebe/dynamics/coord_b2j.py
 delete mode 100755 phoebe/dynamics/coord_h2j.py
 delete mode 100755 phoebe/dynamics/coord_j2b.py
 mode change 100755 => 100644 phoebe/dynamics/geometry.py
 delete mode 100755 phoebe/dynamics/invgeometry.py
 mode change 100755 => 100644 phoebe/dynamics/keplerian.py
 rename phoebe/dynamics/{orbel_xv2el.py => orbel.py} (55%)
 mode change 100755 => 100644
 delete mode 100755 phoebe/dynamics/orbel_ehie.py
 delete mode 100755 phoebe/dynamics/orbel_el2xv.py

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/coord_b2j.py b/phoebe/dynamics/coord_b2j.py
deleted file mode 100755
index 4e0387ddb..000000000
--- a/phoebe/dynamics/coord_b2j.py
+++ /dev/null
@@ -1,60 +0,0 @@
-#!/usr/bin/env python3
-
-import numpy as np
-
-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 = np.zeros((nbod))
-    rj = np.zeros((nbod, 3))
-    vj = np.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
-
-if __name__ == "__main__":
-
-    mass = np.array([1.0, 1.0, 0.5])
-    rb = np.array([[-0.5, -0.2, 0.0], [0.5, -0.2, 0.0], [0.0, 0.8, 0.0]])
-    vb = np.array([[-0.5, -0.2, 0.0], [0.5, -0.2, 0.0], [0.0, 0.8, 0.0]])
-
-    print("mass = ", mass)
-    print("rb[0] = ", rb[0])
-    print("rb[1] = ", rb[1])
-    print("rb[2] = ", rb[2])
-
-    rj, vj = coord_b2j(mass, rb, vb)
-
-    print("rj[0] = ", rj[0])
-    print("rj[1] = ", rj[1])
-    print("rj[2] = ", rj[2])
-
diff --git a/phoebe/dynamics/coord_h2j.py b/phoebe/dynamics/coord_h2j.py
deleted file mode 100755
index a8bba1369..000000000
--- a/phoebe/dynamics/coord_h2j.py
+++ /dev/null
@@ -1,61 +0,0 @@
-#!/usr/bin/env python3
-
-import numpy as np
-
-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 = np.zeros((nbod))
-    rj = np.zeros((nbod, 3))
-    vj = np.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
-
-if __name__ == "__main__":
-
-    mass = np.array([1.0, 1.0, 0.5])
-    rh = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]])
-    vh = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]])
-
-    print("mass = ", mass)
-    print("rh[0] = ", rh[0])
-    print("rh[1] = ", rh[1])
-    print("rh[2] = ", rh[2])
-
-    rj, vj = coord_h2j(mass, rh, vh)
-
-    print("rj[0] = ", rj[0])
-    print("rj[1] = ", rj[1])
-    print("rj[2] = ", rj[2])
-
-
diff --git a/phoebe/dynamics/coord_j2b.py b/phoebe/dynamics/coord_j2b.py
deleted file mode 100755
index d2f683b06..000000000
--- a/phoebe/dynamics/coord_j2b.py
+++ /dev/null
@@ -1,77 +0,0 @@
-#!/usr/bin/env python3
-
-import numpy as np
-
-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 = np.zeros((nbod))
-    rb = np.zeros((nbod, 3))
-    vb = np.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 = np.array([0.0, 0.0, 0.0])
-    vtmp = np.array([0.0, 0.0, 0.0])
-
-    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
-
-if __name__ == "__main__":
-
-    mass = np.array([1.0, 1.0, 0.5])
-    rj = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]])
-    vj = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]])
-
-    print("mass = ", mass)
-    print("rj[0] = ", rj[0])
-    print("rj[1] = ", rj[1])
-    print("rj[2] = ", rj[2])
-
-    rb, vb = coord_j2b(mass, rj, vj)
-
-    print("rb[0] = ", rb[0])
-    print("rb[1] = ", rb[1])
-    print("rb[2] = ", rb[2])
-
-    com = np.array([0.0, 0.0, 0.0])
-    for i in range(3):
-        com[i] = np.dot(mass, rb[:,i])/np.sum(mass)
-
-    print("com = ", com)
-
diff --git a/phoebe/dynamics/geometry.py b/phoebe/dynamics/geometry.py
old mode 100755
new mode 100644
index 827ba2cc5..cfa5894a7
--- a/phoebe/dynamics/geometry.py
+++ b/phoebe/dynamics/geometry.py
@@ -1,9 +1,9 @@
 #!/usr/bin/env python3
 
 import numpy as np
-from phoebe.dynamics import coord_h2j
-from phoebe.dynamics import coord_j2b
-from phoebe.dynamics import orbel_el2xv
+
+from phoebe.dynamics import coord
+from phoebe.dynamics import orbel
 
 def geometry(m, elmts, geometry='hierarchical'):
     """
@@ -13,11 +13,11 @@ def geometry(m, elmts, geometry='hierarchical'):
 
     if geometry == 'hierarchical':
 
-        _geometry = hierarchical
+        _geometry = geometry_hierarchical
 
     elif geometry == 'twopairs':
 
-        _geometry = twopairs
+        _geometry = geometry_twopairs
 
     else:
         raise NotImplementedError
@@ -47,7 +47,7 @@ def geometry(m, elmts, geometry='hierarchical'):
 
     return xs, ys, zs, vxs, vys, vzs
 
-def hierarchical(m, elmts):
+def geometry_hierarchical(m, elmts):
     """
      _          \                |
     / \          |               |
@@ -66,14 +66,14 @@ def hierarchical(m, elmts):
         msum += m[j]
         ialpha = -1
 
-        rj[j], vj[j] = orbel_el2xv.orbel_el2xv(msum, ialpha, elmts[j-1])
+        rj[j], vj[j] = orbel.orbel_el2xv(msum, ialpha, elmts[j-1])
 
     # convert to barycentric frame
-    rb, vb = coord_j2b.coord_j2b(m, rj, vj)
+    rb, vb = coord.coord_j2b(m, rj, vj)
 
     return rb, vb
 
-def twopairs(m, elmts):
+def geometry_twopairs(m, elmts):
     """
      _          _               \ 
     / \        / \               |
@@ -91,7 +91,7 @@ def twopairs(m, elmts):
     # (1+2) pair, 1-centric coordinates
     msum = m[0]+m[1]
     ialpha = -1
-    r2_1, v2_1 = orbel_el2xv.orbel_el2xv(msum, ialpha, elmts[0])
+    r2_1, v2_1 = orbel.orbel_el2xv(msum, ialpha, elmts[0])
 
     # barycenter
     r12_1 = m[1]/msum * r2_1
@@ -99,7 +99,7 @@ def twopairs(m, elmts):
 
     # (3+4) pair, 3-centric
     msum = m[2]+m[3]
-    r4_3, v4_3 = orbel_el2xv.orbel_el2xv(msum, ialpha, elmts[1])
+    r4_3, v4_3 = orbel.orbel_el2xv(msum, ialpha, elmts[1])
 
     # barycenter
     r34_3 = m[3]/msum * r4_3
@@ -107,7 +107,7 @@ def twopairs(m, elmts):
 
     # (1+2)+(3+4) mutual orbit, (1+2)-centric
     msum = m[0]+m[1]+m[2]+m[3]
-    r34_12, v34_12 = orbel_el2xv.orbel_el2xv(msum, ialpha, elmts[2])
+    r34_12, v34_12 = orbel.orbel_el2xv(msum, ialpha, elmts[2])
 
     # everything to 1-centric
     rh[0,:] = 0.0
@@ -120,55 +120,152 @@ def twopairs(m, elmts):
     vh[3,:] = vh[2,:] + v4_3
 
     # everything to Jacobian
-    rj[0:4], vj[0:4] = coord_h2j.coord_h2j(m[0:4], rh[0:4], vh[0:4])
+    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_el2xv.orbel_el2xv(msum, ialpha, elmts[j-1])
+        rj[j], vj[j] = orbel.orbel_el2xv(msum, ialpha, elmts[j-1])
 
     # convert to barycentric frame
-    rb, vb = coord_j2b.coord_j2b(m, rj, vj)
+    rb, vb = coord.coord_j2b(m, rj, vj)
 
     return rb, vb
 
-if __name__ == "__main__":
-
-    day = 86400.
-    au = 1.496e11
-    M_S = 2.e30
-    G = 6.67e-11
-    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 = hierarchical(m, elmts)
-
-    print("m = ", m)
-    print("rb[0] = ", rb[0])
-    print("rb[1] = ", rb[1])
-    print("vb[0] = ", vb[0])
-    print("vb[1] = ", vb[1])
-
-    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 = twopairs(m, elmts)
-
-    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])
+########################################################################
+
+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/invgeometry.py b/phoebe/dynamics/invgeometry.py
deleted file mode 100755
index b318aa0a4..000000000
--- a/phoebe/dynamics/invgeometry.py
+++ /dev/null
@@ -1,192 +0,0 @@
-#!/usr/bin/env python3
-
-import numpy as np
-from phoebe.dynamics import coord_b2j
-from phoebe.dynamics import orbel_xv2el
-
-def invgeometry(m, rb, vb, geometry='hierarchical'):
-    """
-    Convert barycentric coordinates to elements for various geometries.
-   
-    """
-
-    if geometry == 'hierarchical':
-
-        _invgeometry = hierarchical
-
-    elif geometry == 'twopairs':
-
-        _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 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_b2j.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_xv2el.orbel_xv2el(msum, rj[j], vj[j])
-        euler[j,:] = orbel_xv2el.get_euler(msum, elmts[j-1])
-        roche[j,:] = orbel_xv2el.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 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_xv2el.orbel_xv2el(msum, r2_1, v2_1)
-    euler[0,:] = orbel_xv2el.get_euler(msum, elmts[0])
-    roche[0,:] = orbel_xv2el.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_xv2el.orbel_xv2el(msum, r4_3, v4_3)
-    euler[2,:] = orbel_xv2el.get_euler(msum, elmts[1])
-    roche[2,:] = orbel_xv2el.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_xv2el.orbel_xv2el(msum, r34_12, v34_12)
-
-    # everything to Jacobian
-    rj, vj = coord_b2j.coord_b2j(m, rb, vb)
-    
-    # other bodies (also Jacobian)
-    for j in range(4, nbod):
-        msum += m[j]
-
-        elmts[j-1,:] = orbel_xv2el.orbel_xv2el(msum, rj[j], vj[j])
-        euler[j,:] = orbel_xv2el.get_euler(msum, elmts[j-1])
-        roche[j,:] = orbel_xv2el.get_roche(msum, elmts[j-1])
-
-        euler[j,0] += np.pi
-
-    return elmts, euler, roche
-
-if __name__ == "__main__":
-
-    day = 86400.
-    au = 1.496e11
-    M_S = 2.e30
-    G = 6.67e-11
-    gms = G*M_S / (au**3 * day**-2)
-    vk = np.sqrt(gms/1.0)
-
-    m = gms*np.array([1.0, 0.0])
-    rb = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]])
-    vb = np.array([[0.0, 0.0, 0.0], [0.0, vk, 0.0]])
-
-    elmts, euler, roche = hierarchical(m, rb, vb)
-
-    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])
-
-    m = gms*np.array([1.0, 1.0, 0.5, 0.5])
-    rb = np.array([ \
-        [-0.38333333, 0.0, 0.0], \
-        [-0.28333333, 0.0, 0.0], \
-        [ 0.56666667, 0.0, 0.0], \
-        [ 0.76666667, 0.0, 0.0], \
-        ])
-    vb = np.array([ \
-        [0.0, -0.04852087, 0.0], \
-        [0.0,  0.02860663, 0.0], \
-        [0.0,  0.00063236, 0.0], \
-        [0.0,  0.03919611, 0.0], \
-        ])
-
-    elmts, euler, roche = twopairs(m, rb, vb)
-
-    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])
-
-
diff --git a/phoebe/dynamics/keplerian.py b/phoebe/dynamics/keplerian.py
old mode 100755
new mode 100644
index 69d6b21c9..2921e73c9
--- a/phoebe/dynamics/keplerian.py
+++ b/phoebe/dynamics/keplerian.py
@@ -5,9 +5,8 @@
 
 from phoebe import u, c
 
-from phoebe.dynamics import coord_j2b
-from phoebe.dynamics import orbel_el2xv
-from phoebe.dynamics import orbel_ehie
+from phoebe.dynamics import coord
+from phoebe.dynamics import orbel
 
 import logging
 logger = logging.getLogger("DYNAMICS.KEPLERIAN")
@@ -148,13 +147,13 @@ def _keplerian(times, component=None):
 
                 elmts = [a, eccs[j-1], incls[j-1], long_ans[j-1], omega, M]
 
-                rj[j], vj[j] = orbel_el2xv.orbel_el2xv(msum, ialpha, elmts)
+                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_el2xv.get_euler(elmts)
+                    euler[j] = orbel.get_euler(elmts)
 
                     if j==1:
                         euler[0,:] = euler[1,:]
@@ -162,7 +161,7 @@ def _keplerian(times, component=None):
                         euler[j,0] += np.pi
 
             # convert to barycentric frame
-            rb, vb = coord_j2b.coord_j2b(masses, rj, vj)
+            rb, vb = coord.coord_j2b(masses, rj, vj)
 
             # gamma velocity
             rb[:,2] -= vgamma*(t-t0)
@@ -234,50 +233,4 @@ def propertime_residual(t, time):
     else:
         return times, xs, ys, zs, vxs, vys, vzs
 
-if __name__ == "__main__":
-
-    # Sun-Earth
-    times = np.array([0.0])
-    masses = np.array([1.0, 0.0])
-    smas = np.array([1.0])
-    eccs = np.array([0.0])
-    incls = np.array([0.0])
-    per0s = np.array([0.0])
-    long_ans = np.array([0.0])
-    mean_anoms = np.array([0.0])
-
-    G = c.G.to('AU3 / (Msun d2)').value
-    masses *= G
-
-    times, xs, ys, zs, vxs, vys, vzs = dynamics(times, masses, smas, eccs, incls, per0s, long_ans, mean_anoms)
-
-    print("xs = ", xs)
-    print("ys = ", ys)
-    print("zs = ", zs)
-    print("vxs = ", vxs)
-    print("vys = ", vys)
-    print("vzs = ", vzs)
-
-    print("v_of_Earth = ", vys[1][0]*(u.au/u.d).to('m s^-1'), " m/s")
-
-    # 3-body problem
-    times = np.array([0.0])
-    masses = np.array([1.0, 1.0, 1.0])
-    smas = np.array([1.0, 10.0])
-    eccs = np.array([0.0, 0.0])
-    incls = np.array([0.0, 0.0])
-    per0s = np.array([0.0, 0.0])
-    long_ans = np.array([0.0, 0.0])
-    mean_anoms = np.array([0.0, 0.0])
-
-    times, xs, ys, zs, vxs, vys, vzs = dynamics(times, masses, smas, eccs, incls, per0s, long_ans, mean_anoms)
-
-    print("")
-    print("xs = ", xs)
-    print("ys = ", ys)
-    print("zs = ", zs)
-    print("vxs = ", vxs)
-    print("vys = ", vys)
-    print("vzs = ", vzs)
-
 
diff --git a/phoebe/dynamics/orbel_xv2el.py b/phoebe/dynamics/orbel.py
old mode 100755
new mode 100644
similarity index 55%
rename from phoebe/dynamics/orbel_xv2el.py
rename to phoebe/dynamics/orbel.py
index 24bb21755..16ea69baf
--- a/phoebe/dynamics/orbel_xv2el.py
+++ b/phoebe/dynamics/orbel.py
@@ -1,12 +1,91 @@
-#!/usr/bin/env python3
 
 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.
@@ -104,28 +183,6 @@ def orbel_xv2el(gm, r, v):
 
     return elmts
 
-def get_euler(gm, 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 get_roche(gm, elmts):
     """
     Get corresponding Roche parameters.
@@ -150,23 +207,51 @@ def get_roche(gm, elmts):
 
     return roche
 
-if __name__ == "__main__":
 
-    day = 86400.
-    au = 1.496e11
-    M_S = 2.e30
-    G = 6.67e-11
-    gms = G*M_S / (au**3 * day**-2)
+def orbel_ehie(e, m):
+    """
+    Solves Kepler's equation.
+
+    Rewritten from swift/orbel/orbel_ehie.f.
 
-    r = array([0.9, 0.0, 0.0])
-    v = array([0.0, 1.9066428749416830523700e-02, 0.0])
+    Reference: Levison, H.F., Duncan, M.J., The long-term dynamical behavior of short-period comets. Icarus 108, 18-36, 1994.
 
-    print("gms = ", gms)
-    print("r = ", r)
-    print("v = ", v)
+    """
 
-    elmts = orbel_xv2el(gms, r, v)
-    
-    print("elmts = ", elmts)
+    # 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/orbel_ehie.py b/phoebe/dynamics/orbel_ehie.py
deleted file mode 100755
index 362d32c2e..000000000
--- a/phoebe/dynamics/orbel_ehie.py
+++ /dev/null
@@ -1,68 +0,0 @@
-#!/usr/bin/env python3
-
-from numpy import sin, cos, pi
-
-NMAX = 3
-TWOPI = 2.0*pi
-
-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
-    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
-
-if __name__ == "__main__":
-
-    m = 0.25
-    e = 0.6
-
-    print("m = ", m)
-    print("e = ", e)
-
-    cape = orbel_ehie(e, m)
-
-    print("cape = ", cape)
-
-    print("")
-    print("M = E - e*sin(E)")
-    print(m, " = ", cape - e*sin(cape))
-
diff --git a/phoebe/dynamics/orbel_el2xv.py b/phoebe/dynamics/orbel_el2xv.py
deleted file mode 100755
index c18d44415..000000000
--- a/phoebe/dynamics/orbel_el2xv.py
+++ /dev/null
@@ -1,108 +0,0 @@
-#!/usr/bin/env python3
-
-from numpy import sin, cos, sqrt, tan, arctan, array
-
-from phoebe.dynamics import orbel_ehie
-
-TINY = 4.0e-15
-
-cape = 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.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
-
-if __name__ == "__main__":
-
-    day = 86400.
-    au = 1.496e11
-    M_S = 2.e30
-    G = 6.67e-11
-    gms = G*M_S / (au**3 * day**-2)
-
-    elmts = [1.0, 0.1, 0.0, 0.0, 0.0, 0.0]
-
-    print("gms = ", gms)
-    print("elmts = ", elmts)
-
-    r, v = orbel_el2xv(gms, -1, elmts)
-
-    print("r = ", r)
-    print("v = ", v)
-    print("v[1] = %.22e" % v[1])
-
-

From c29971988cf4596f515e87afb3360c9ec2781650 Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Mon, 14 Oct 2024 18:35:38 +0200
Subject: [PATCH 34/39] Correction of import.

---
 phoebe/dynamics/xyz.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/phoebe/dynamics/xyz.py b/phoebe/dynamics/xyz.py
index 45d343e63..57e33c9a2 100644
--- a/phoebe/dynamics/xyz.py
+++ b/phoebe/dynamics/xyz.py
@@ -7,7 +7,6 @@
 from phoebe import u, c
 from phoebe import conf
 from phoebe.dynamics import geometry
-from phoebe.dynamics import invgeometry
 
 try:
     import rebound
@@ -255,7 +254,7 @@ def dynamics(times, masses, xi, yi, zi, vxi, vyi, vzi, \
 
         if return_roche_euler:
 
-            elmts, euler, roche = invgeometry.invgeometry(masses, rb, vb, geometry=_geometry)
+            elmts, euler, roche = geometry.invgeometry(masses, rb, vb, geometry=_geometry)
 
         fac = (1*u.AU).to(u.solRad).value
 

From de0daae070dbcc6694d9e5327be81dfcc4d8053e Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Mon, 14 Oct 2024 22:53:22 +0200
Subject: [PATCH 35/39] Test of geometry.

---
 tests/tests/test_geometry/test_geometry.py | 101 +++++++++++++++++++++
 1 file changed, 101 insertions(+)
 create mode 100755 tests/tests/test_geometry/test_geometry.py

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)
+
+

From b5e61b124fc2834d28435193aceccb672d9c181a Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Mon, 21 Oct 2024 21:24:31 +0200
Subject: [PATCH 36/39] Generalisation of J2 = -C20 in reboundx. Each component
 has its own rotation axis, with respect to which the acceleration must be
 computed (and transformed back to the inertial frame). For the moment,
 src/gravitational_harmonics.c must be modified:

https://sirrah.troja.mff.cuni.cz/~mira/tmp/phoebe2/reboundx_j2_/
---
 phoebe/dynamics/xyz.py | 32 +++++++++++++++++++++++++++++---
 1 file changed, 29 insertions(+), 3 deletions(-)

diff --git a/phoebe/dynamics/xyz.py b/phoebe/dynamics/xyz.py
index 57e33c9a2..5d7d2674f 100644
--- a/phoebe/dynamics/xyz.py
+++ b/phoebe/dynamics/xyz.py
@@ -30,6 +30,22 @@
 
 _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`.
@@ -89,8 +105,14 @@ def dynamics_from_bundle(b, times, compute=None, return_roche_euler=False, **kwa
     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]])
@@ -109,6 +131,7 @@ def dynamics_from_bundle(b, times, compute=None, return_roche_euler=False, **kwa
                j2=j2, \
                j2s=j2s, \
                requivs=requivs, \
+               spins=spins, \
                return_roche_euler=return_roche_euler \
                )
 
@@ -125,6 +148,7 @@ def dynamics(times, masses, xi, yi, zi, vxi, vyi, vzi, \
         j2=False, \
         j2s=None, \
         requivs=None, \
+        spins=None, \
         return_roche_euler=False \
         ):
 
@@ -153,6 +177,7 @@ def dynamics(times, masses, xi, yi, zi, vxi, vyi, vzi, \
         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:
@@ -222,14 +247,15 @@ def dynamics(times, masses, xi, yi, zi, vxi, vyi, vzi, \
         rebx.add_force(gr)
 
     if j2:
-        logger.info("enabling 'j2' in reboundx")
+        logger.info("enabling 'gravitational_harmonics' in reboundx")
         rebx = reboundx.Extras(sim)
-        gh = rebx.load_force("j2")
+        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))

From 214b0fb7c3622a4d11b0f5c8e4151f0f311f10bc Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Tue, 22 Oct 2024 14:57:08 +0200
Subject: [PATCH 37/39] Round-offs in elements.

---
 phoebe/dynamics/orbel.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/phoebe/dynamics/orbel.py b/phoebe/dynamics/orbel.py
index 16ea69baf..a4fe5072e 100644
--- a/phoebe/dynamics/orbel.py
+++ b/phoebe/dynamics/orbel.py
@@ -106,7 +106,7 @@ def orbel_xv2el(gm, r, v):
 
     fac = sqrt(h[0]*h[0] + h[1]*h[1])/absh
 
-    if fac < TINY:
+    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:

From c648a281f1eae45719f7cdc47db1fa13b8f52491 Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Thu, 7 Nov 2024 14:53:54 +0100
Subject: [PATCH 38/39] Mass constraint for multiples.

---
 phoebe/frontend/bundle.py       | 10 ++--------
 phoebe/parameters/constraint.py | 25 +++++++++++++++++++------
 2 files changed, 21 insertions(+), 14 deletions(-)

diff --git a/phoebe/frontend/bundle.py b/phoebe/frontend/bundle.py
index 9b082e37b..d7a9164dc 100644
--- a/phoebe/frontend/bundle.py
+++ b/phoebe/frontend/bundle.py
@@ -2840,14 +2840,8 @@ def set_hierarchy(self, *args, **kwargs):
                                         solve_for=constraint_param.constrained_parameter.uniquetwig,
                                         constraint=constraint_param.constraint)
                 else:
-                    # NOTE: IN ORDER TO DEFAULT_TRIPLE() WORKS
-                    sibling = self.hierarchy.get_sibling_of(component)
-                    kind = self.hierarchy.get_kind_of(sibling)
-                    if kind != 'star':
-                        logger.warning('constraint mass not working for multiple systems')
-                    else:
-                        self.add_constraint(constraint.mass, component,
-                                            constraint=self._default_label('mass', context='constraint'))
+                    self.add_constraint(constraint.mass, component,
+                                        constraint=self._default_label('mass', context='constraint'))
 
 
                 logger.debug('re-creating comp_sma constraint for {}'.format(component))
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')

From c28e09e0545a1357f9d9dcb798ab7ca6c000378d Mon Sep 17 00:00:00 2001
From: Miroslav Broz <miroslav.broz@email.cz>
Date: Tue, 12 Nov 2024 08:38:19 +0100
Subject: [PATCH 39/39] Quadratic term in the mean anomaly.

---
 phoebe/dynamics/keplerian.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/phoebe/dynamics/keplerian.py b/phoebe/dynamics/keplerian.py
index 2921e73c9..3d5e7fff7 100644
--- a/phoebe/dynamics/keplerian.py
+++ b/phoebe/dynamics/keplerian.py
@@ -141,6 +141,7 @@ def _keplerian(times, component=None):
                 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)