Implementation of Barycentric Coordinates

rebound/integrators/whfast.py

@@ -2,7 +2,7 @@
 from ..particle import Particle
 
 WHFAST_KERNELS = {"default": 0, "modifiedkick": 1, "composition": 2, "lazy": 3}
-WHFAST_COORDINATES = {"jacobi": 0, "democraticheliocentric": 1, "whds": 2}
+WHFAST_COORDINATES = {"jacobi": 0, "democraticheliocentric": 1, "whds": 2, "barycentric": 3}
 
 class IntegratorWHFast(ctypes.Structure):
     """
@@ -79,6 +79,7 @@ def coordinates(self):
         - ``'jacobi'`` (default)
         - ``'democraticheliocentric'``
         - ``'whds'``
+        - ``'barycentric'``
         """
         i = self._coordinates
         for name, _i in WHFAST_COORDINATES.items():

rebound/tests/test_transformations.py

@@ -17,6 +17,38 @@ def getc(sim):
     return c
 
 class TestTransformations(unittest.TestCase):
+
+    def test_barycentric(self):
+        sim = rebound.Simulation()
+        sim.add(m=1.2354)
+        sim.add(m=0.1,a=1.24,e=0.123,inc=0.14,omega=0.12,Omega=0.64,l=0.632)
+        sim.add(m=0.01,a=5.24,e=0.2123,inc=0.014,omega=0.012,Omega=0.0164,l=10.18632)
+        sim.add(m=1e-7,a=7.24,e=0.22123,inc=0.3014,omega=0.4012,Omega=0.110164,l=2.18632)
+
+        elems = (rebound.Particle * sim.N)()
+        p = ctypes.cast(elems,ctypes.POINTER(rebound.Particle))
+
+        c0 = getc(sim)
+        cl = rebound.clibrebound
+        cl.reb_particles_transform_inertial_to_barycentric_posvel(sim._particles,p,sim.N,sim.N)
+
+        for i in range(sim.N):
+            sim.particles[i].x = 1234.
+            sim.particles[i].vx = 1234.
+        cl.reb_particles_transform_barycentric_to_inertial_posvel(sim._particles,p,sim.N,sim.N)
+
+        c1 = getc(sim)
+
+        for i in range(len(c0)):
+            self.assertAlmostEqual(c0[i],c1[i],delta=1e-16)
+
+        for i in range(sim.N):
+            sim.particles[i].x = 1234.
+        cl.reb_particles_transform_barycentric_to_inertial_pos(sim._particles,p,sim.N,sim.N)
+
+        c1 = getc(sim)
+        for i in range(len(c0)):
+            self.assertAlmostEqual(c0[i],c1[i],delta=1e-16)
 
     def test_democratichelio(self):
         sim = rebound.Simulation()

rebound/tests/test_whfast_testparticles.py

@@ -3,7 +3,7 @@
 import rebound.data
 import warnings
 
-coordinatelist = ["democraticheliocentric","whds","jacobi"]
+coordinatelist = ["democraticheliocentric","whds","jacobi","barycentric"]
 class TestIntegratorWHFastTestParticle(unittest.TestCase):
     pass
 

src/integrator_whfast.c

@@ -427,6 +427,14 @@ void reb_whfast_interaction_step(struct reb_simulation* const r, const double _d
                 p_j[i].vz += _dt*particles[i].az;
             }
             break;
+        case REB_WHFAST_COORDINATES_BARYCENTRIC:
+#pragma omp parallel for
+            for (unsigned int i = 0; i < N_real; i++) {
+                p_j[i].vx += _dt * particles[i].ax;
+                p_j[i].vy += _dt * particles[i].ay;
+                p_j[i].vz += _dt * particles[i].az;
+            }
+            break;
     };
 }
 void reb_whfast_jump_step(const struct reb_simulation* const r, const double _dt){
@@ -482,6 +490,9 @@ void reb_whfast_jump_step(const struct reb_simulation* const r, const double _dt
             }
             }
             break;
+        case REB_WHFAST_COORDINATES_BARYCENTRIC:
+            // Nothing to be done.
+            break;
     };
 }
 
@@ -522,7 +533,15 @@ void reb_whfast_kepler_step(const struct reb_simulation* const r, const double _
                 }
                 reb_whfast_kepler_solver(r, p_j, eta*G, i, _dt);
             }
-            break;
+        break;
+        case REB_WHFAST_COORDINATES_BARYCENTRIC:
+            for (unsigned int i=1;i<N_active;i++) {
+                eta+=p_j[i].m;
+            }
+            for (unsigned int i=1;i<N_real;i++) {
+                reb_whfast_kepler_solver(r, p_j,eta*G, i, _dt);
+            }
+        break;
     };
 }
 
@@ -823,6 +842,9 @@ void reb_integrator_whfast_from_inertial(struct reb_simulation* const r){
         case REB_WHFAST_COORDINATES_WHDS:
             reb_particles_transform_inertial_to_whds_posvel(particles, ri_whfast->p_jh, N_real, N_active);
             break;
+        case REB_WHFAST_COORDINATES_BARYCENTRIC:
+            reb_particles_transform_inertial_to_barycentric_posvel(particles, ri_whfast->p_jh, N_real, N_active);
+            break;
     };
 }
 
@@ -845,6 +867,9 @@ void reb_integrator_whfast_to_inertial(struct reb_simulation* const r){
             case REB_WHFAST_COORDINATES_WHDS:
                 reb_particles_transform_whds_to_inertial_posvel(particles, ri_whfast->p_jh, N_real, N_active);
                 break;
+            case REB_WHFAST_COORDINATES_BARYCENTRIC:
+                reb_particles_transform_barycentric_to_inertial_posvel(particles, ri_whfast->p_jh, N_real, N_active);
+                break;
         };
     }else{
         switch (ri_whfast->coordinates){
@@ -857,6 +882,9 @@ void reb_integrator_whfast_to_inertial(struct reb_simulation* const r){
             case REB_WHFAST_COORDINATES_WHDS:
                 reb_particles_transform_whds_to_inertial_posvel(particles, ri_whfast->p_jh, N_real, N_active);
                 break;
+            case REB_WHFAST_COORDINATES_BARYCENTRIC:
+                reb_particles_transform_barycentric_to_inertial_posvel(particles, ri_whfast->p_jh, N_real, N_active);
+                break;
         };
     }
 }
@@ -1000,6 +1028,9 @@ void reb_integrator_whfast_synchronize(struct reb_simulation* const r){
             case REB_WHFAST_COORDINATES_WHDS:
                 reb_particles_transform_whds_to_inertial_posvel(r->particles, ri_whfast->p_jh, N_real, N_active);
                 break;
+            case REB_WHFAST_COORDINATES_BARYCENTRIC:
+                reb_particles_transform_barycentric_to_inertial_posvel(r->particles, ri_whfast->p_jh, N_real, N_active);
+                break;
         };
         for (int v=0;v<r->N_var_config;v++){
             struct reb_variational_configuration const vc = r->var_config[v];

src/rebound.h

@@ -342,6 +342,7 @@ struct reb_integrator_whfast {
         REB_WHFAST_COORDINATES_JACOBI = 0,                      // Jacobi coordinates (default)
         REB_WHFAST_COORDINATES_DEMOCRATICHELIOCENTRIC = 1,      // Democratic Heliocentric coordinates
         REB_WHFAST_COORDINATES_WHDS = 2,                        // WHDS coordinates (Hernandez and Dehnen, 2017)
+        REB_WHFAST_COORDINATES_BARYCENTRIC = 3,                 // Barycentric coordinates
     } coordinates;                                              // Coordinate system used in Hamiltonian splitting
     unsigned int recalculate_coordinates_this_timestep;         // 1: recalculate coordinates from inertial coordinates
     unsigned int safe_mode;                                     // 0: Drift Kick Drift scheme (default), 1: combine first and last sub-step.
@@ -1384,6 +1385,10 @@ DLLEXPORT void reb_particles_transform_inertial_to_whds_posvel(const struct reb_
 DLLEXPORT void reb_particles_transform_whds_to_inertial_pos(struct reb_particle* const particles, const struct reb_particle* const p_h, const unsigned int N, const unsigned int N_active);
 DLLEXPORT void reb_particles_transform_whds_to_inertial_posvel(struct reb_particle* const particles, const struct reb_particle* const p_h, const unsigned int N, const unsigned int N_active);
 
+// Barycentric coordinates
+DLLEXPORT void reb_particles_transform_inertial_to_barycentric_posvel(const struct reb_particle* const particles, struct reb_particle* const p_b, const unsigned int N, const unsigned int N_active);
+DLLEXPORT void reb_particles_transform_barycentric_to_inertial_pos(struct reb_particle* const particles, const struct reb_particle* const p_b, const unsigned int N, const unsigned int N_active);
+DLLEXPORT void reb_particles_transform_barycentric_to_inertial_posvel(struct reb_particle* const particles, const struct reb_particle* const p_b, const unsigned int N, const unsigned int N_active);
 
 // Temporary. Function declarations needed by REBOUNDx 
 DLLEXPORT void reb_integrator_ias15_reset(struct reb_simulation* r);         ///< Internal function used to call a specific integrator

src/transformations.c

@@ -472,3 +472,109 @@ void reb_particles_transform_democraticheliocentric_to_inertial_posvel(struct re
     particles[0].vz = p_h[0].vz -vz0;
 }
 
+/***************
+* Barycentric */
+
+void reb_particles_transform_barycentric_to_inertial_pos(struct reb_particle* const particles, const struct reb_particle* const p_b, const unsigned int N, const unsigned int N_active){
+    double mtot = 0;
+    double x0 = 0;
+    double y0 = 0;
+    double z0 = 0;
+#pragma omp parallel for reduction(+:mtot)
+    for (unsigned int i = 0; i<N_active; i++) {
+        const double m = p_b[i].m;
+        particles[i].m = m;
+        mtot += p_b[i].m;
+    }
+    const double mi = 1.0 / mtot;
+#pragma omp parallel for reduction(+:x0) reduction(+:y0) reduction(+:z0)
+    for (unsigned int i=1; i<N; i++) {
+        const double m = p_b[i].m;
+        particles[i].x = p_b[i].x + p_b[0].x*mi;
+        particles[i].y = p_b[i].y + p_b[0].y*mi;
+        particles[i].z = p_b[i].z + p_b[0].z*mi;
+        if(i < N_active){
+            x0 += particles[i].x*m;
+            y0 += particles[i].y*m;
+            z0 += particles[i].z*m;
+        }
+    }
+    particles[0].m = p_b[0].m;
+    particles[0].x = (p_b[0].x - x0)/p_b[0].m;
+    particles[0].y = (p_b[0].y - y0)/p_b[0].m;
+    particles[0].z = (p_b[0].z - z0)/p_b[0].m;
+}
+
+void reb_particles_transform_barycentric_to_inertial_posvel(struct reb_particle* const particles, const struct reb_particle* const p_b, const unsigned int N, const unsigned int N_active){
+    double mtot = 0;
+    double vx0 = 0;
+    double vy0 = 0;
+    double vz0 = 0;
+    double x0 = 0;
+    double y0 = 0;
+    double z0 = 0;
+#pragma omp parallel for reduction(+:mtot)
+    for (unsigned int i = 0; i<N_active;i++) {
+        const double m = p_b[i].m;
+        mtot += m;
+        particles[i].m = m;
+    }
+    const double mi = 1./mtot;
+#pragma omp parallel for reduction(+:x0) reduction(+:y0) reduction(+:z0) reduction(+:vx0) reduction(+:vy0) reduction(+:vz0)
+    for (unsigned int i=1; i<N; i++) {
+        particles[i].x = p_b[i].x + p_b[0].x*mi;
+        particles[i].y = p_b[i].y + p_b[0].y*mi;
+        particles[i].z = p_b[i].z + p_b[0].z*mi;
+        particles[i].vx = p_b[i].vx + p_b[0].vx*mi;
+        particles[i].vy = p_b[i].vy + p_b[0].vy*mi;
+        particles[i].vz = p_b[i].vz + p_b[0].vz*mi;
+        if(i < N_active) {
+            const double m = p_b[i].m;
+            x0 += particles[i].x*m;
+            y0 += particles[i].y*m;
+            z0 += particles[i].z*m;
+            vx0 += particles[i].vx*m;
+            vy0 += particles[i].vy*m;
+            vz0 += particles[i].vz*m;
+        }
+    }
+    particles[0].m = p_b[0].m;
+    particles[0].x = (p_b[0].x - x0)/p_b[0].m;
+    particles[0].y = (p_b[0].y - y0)/p_b[0].m;
+    particles[0].z = (p_b[0].z - z0)/p_b[0].m;
+    particles[0].vx = (p_b[0].vx - vx0)/p_b[0].m;
+    particles[0].vy = (p_b[0].vy - vy0)/p_b[0].m;
+    particles[0].vz = (p_b[0].vz - vz0)/p_b[0].m;
+}
+void reb_particles_transform_inertial_to_barycentric_posvel(const struct reb_particle* const particles, struct reb_particle* const p_b, const unsigned int N, const unsigned int N_active) {
+    p_b[0].x = 0;
+    p_b[0].y = 0;
+    p_b[0].z = 0;
+    p_b[0].vx = 0;
+    p_b[0].vy = 0;
+    p_b[0].vz = 0;
+    p_b[0].m = particles[0].m;
+    double mtot = 0;
+#pragma omp parallel for reduction(+:mtot)
+    for (unsigned int i = 0; i < N_active; i++) {
+        const double m = particles[i].m;
+        mtot += m;
+        p_b[i].m = m;
+        p_b[0].x += particles[i].x * m;
+        p_b[0].y += particles[i].y * m;
+        p_b[0].z += particles[i].z * m;
+        p_b[0].vx += particles[i].vx * m;
+        p_b[0].vy += particles[i].vy * m;
+        p_b[0].vz += particles[i].vz * m;
+    }
+    const double mi = 1.0 / mtot;
+#pragma omp parallel for
+    for (unsigned int i = 1;i<N;i++) {
+        p_b[i].x = particles[i].x - p_b[0].x*mi;
+        p_b[i].y = particles[i].y - p_b[0].y*mi;
+        p_b[i].z = particles[i].z - p_b[0].z*mi;
+        p_b[i].vx = particles[i].vx - p_b[0].vx*mi;
+        p_b[i].vy = particles[i].vy - p_b[0].vy*mi;
+        p_b[i].vz = particles[i].vz - p_b[0].vz*mi;
+    }
+}