Dynamics for multis

phoebe/backend/backends.py

@@ -961,6 +961,9 @@ def _compute_intrinsic_system_at_t0(self, b, compute,
             if dynamics_method in ['nbody', 'rebound']:
                 t0, xs0, ys0, zs0, vxs0, vys0, vzs0, inst_ds0, inst_Fs0, ethetas0, elongans0, eincls0 = dynamics.nbody.dynamics_from_bundle(b, [t0], compute, return_roche_euler=True, **kwargs)
 
+            elif dynamics_method in ['xyz']:
+                t0, xs0, ys0, zs0, vxs0, vys0, vzs0, inst_ds0, inst_Fs0, ethetas0, elongans0, eincls0 = dynamics.xyz.dynamics_from_bundle(b, [t0], compute, return_roche_euler=True, **kwargs)
+
             elif dynamics_method == 'bs':
                 # TODO: pass stepsize
                 # TODO: pass orbiterror
@@ -1027,6 +1030,9 @@ def _worker_setup(self, b, compute, times, infolists, **kwargs):
             if dynamics_method in ['nbody', 'rebound']:
                 ts, xs, ys, zs, vxs, vys, vzs, inst_ds, inst_Fs, ethetas, elongans, eincls = dynamics.nbody.dynamics_from_bundle(b, times, compute, return_roche_euler=True, **kwargs)
 
+            elif dynamics_method in ['xyz']:
+                ts, xs, ys, zs, vxs, vys, vzs, inst_ds, inst_Fs, ethetas, elongans, eincls = dynamics.xyz.dynamics_from_bundle(b, times, compute, return_roche_euler=True, **kwargs)
+
             elif dynamics_method == 'bs':
                 # if distortion_method == 'roche':
                     # raise ValueError("distortion_method '{}' not compatible with dynamics_method '{}'".format(distortion_method, dynamics_method))
@@ -1039,6 +1045,8 @@ def _worker_setup(self, b, compute, times, infolists, **kwargs):
             elif dynamics_method=='keplerian':
                 # TODO: make sure that this takes systemic velocity and corrects positions and velocities (including ltte effects if enabled)
                 ts, xs, ys, zs, vxs, vys, vzs, ethetas, elongans, eincls = dynamics.keplerian.dynamics_from_bundle(b, times, compute, return_euler=True, **kwargs)
+                inst_ds = None
+                inst_Fs = None
 
             else:
                 raise NotImplementedError
@@ -1054,6 +1062,8 @@ def _worker_setup(self, b, compute, times, infolists, **kwargs):
             vxs, vys, vzs = [np.zeros(len(times))], [np.zeros(len(times))], [np.zeros(len(times))]
             xs, ys, zs = [np.zeros(len(times))], [np.zeros(len(times))], [np.full(len(times), vgamma)]
             ethetas, elongans, eincls = [np.zeros(len(times))], [np.full(len(times), long_an)], [np.full(len(times), incl)]
+            inst_ds = None
+            inst_Fs = None
 
             for i,t in enumerate(times):
                 zs[0][i] = vgamma*(t-t0)
@@ -1065,7 +1075,8 @@ def _worker_setup(self, b, compute, times, infolists, **kwargs):
                     dynamics_method=dynamics_method,
                     ts=ts, xs=xs, ys=ys, zs=zs,
                     vxs=vxs, vys=vys, vzs=vzs,
-                    ethetas=ethetas, elongans=elongans, eincls=eincls)
+                    ethetas=ethetas, elongans=elongans, eincls=eincls,
+                    inst_ds=inst_ds, inst_Fs=inst_Fs)
 
     def _run_single_time(self, b, i, time, infolist, **kwargs):
         logger.debug("rank:{}/{} PhoebeBackend._run_single_time(i={}, time={}, infolist={}, **kwargs.keys={})".format(mpi.myrank, mpi.nprocs, i, time, infolist, kwargs.keys()))
@@ -1096,10 +1107,12 @@ def _run_single_time(self, b, i, time, infolist, **kwargs):
         if True in [info['needs_mesh'] for info in infolist]:
 
             if dynamics_method in ['nbody', 'rebound']:
-                di = dynamics.at_i(inst_ds, i)
-                Fi = dynamics.at_i(inst_Fs, i)
                 # by passing these along to update_positions, volume conservation will
                 # handle remeshing the stars
+                inst_ds = kwargs.get('inst_ds')
+                inst_Fs = kwargs.get('inst_Fs')
+                di = dynamics.at_i(inst_ds, i)
+                Fi = dynamics.at_i(inst_Fs, i)
             else:
                 # then allow d to be determined from orbit and original sma
                 # and F to remain fixed

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
 
 

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