Janus test

rebound/simulation.py

@@ -1378,6 +1378,12 @@ def integrator_synchronize(self):
         """
         clibrebound.reb_integrator_synchronize(byref(self))
 
+    def integrator_janus_to_int(self):
+        """
+        Converts all floating point positions and velocities to integers manually. Mainly for testting.
+        """
+        clibrebound.reb_integrator_janus_to_int(byref(self))
+
     def tree_update(self):
         """
         Call this function to update the tree structure manually after removing particles.

rebound/tests/test_janus.py

@@ -0,0 +1,54 @@
+import rebound
+import unittest
+import math
+import rebound.data
+import warnings
+
+class TestIntegratorJanus(unittest.TestCase):
+    def test_janus_energy(self):
+        for o, eps in [ (2,1e-4), (4,1e-8), (6,1e-9), (8,1e-11), (10,1e-13)]:
+            sim = rebound.Simulation()
+            sim.add(m=1.)
+            sim.add(m=1e-3,a=1.12313)
+            sim.add(m=1e-3,a=2.32323)
+            sim.move_to_com()
+            sim.dt = 0.25
+            sim.integrator = "janus"
+            sim.ri_janus.order = o
+            sim.ri_janus.scale_pos = 1e16
+            sim.ri_janus.scale_vel = 1e16
+            e0 = sim.calculate_energy()
+            sim.integrate(1e2)
+            e1 = sim.calculate_energy()
+            self.assertLess(math.fabs((e0-e1)/e1),eps)
+
+    def test_janus_reverse(self):
+        for o in [2,4,6,8,10]:
+            sim = rebound.Simulation()
+            sim.add(m=1.)
+            sim.add(m=1e-3,a=1.12313,omega=0.32643,l=0.3788,e=0.012)
+            sim.add(m=1e-3,a=2.32323,omega=0.12314,l=0.1726,e=0.103)
+            sim.move_to_com()
+            sim.dt = 0.25
+            sim.integrator = "janus"
+            sim.ri_janus.order = o
+            sim.ri_janus.safe_mode = 0
+            sim.ri_janus.scale_pos = 1e16
+            sim.ri_janus.scale_vel = 1e16
+
+            sim.integrator_janus_to_int()
+            x1, x2 = sim.particles[1].x, sim.particles[2].x
+
+
+            sim.integrate(1e2,exact_finish_time=0)
+            sim.dt *= -1
+            sim.ri_janus.is_synchronized = 0
+            sim.integrate(0,exact_finish_time=0)
+
+            xf1, xf2 = sim.particles[1].x, sim.particles[2].x
+
+            self.assertEqual(x1,xf1)
+            self.assertEqual(x2,xf2)
+
+if __name__ == "__main__":
+    unittest.main()

src/integrator_janus.c

@@ -258,10 +258,22 @@ void reb_integrator_janus_part2(struct reb_simulation* r){
 }
 
 void reb_integrator_janus_synchronize(struct reb_simulation* r){
-    if(r->ri_janus.is_synchronized==0){
-        to_double(r->particles, r->ri_janus.p_int, r->N, r->ri_janus.scale_pos, r->ri_janus.scale_vel); 
-        r->ri_janus.is_synchronized = 1;
+    struct reb_simulation_integrator_janus* ri_janus = &(r->ri_janus);
+    if(ri_janus->is_synchronized==0){
+        to_double(r->particles, ri_janus->p_int, r->N, ri_janus->scale_pos, ri_janus->scale_vel); 
+        ri_janus->is_synchronized = 1;
+    }
+}
+
+void reb_integrator_janus_to_int(struct reb_simulation* r){
+    struct reb_simulation_integrator_janus* ri_janus = &(r->ri_janus);
+    if (ri_janus->allocated_N != r->N){
+        ri_janus->allocated_N = r->N;
+        ri_janus->p_int = realloc(ri_janus->p_int, sizeof(struct reb_particle_int)*r->N);
     }
+    to_int(ri_janus->p_int, r->particles, r->N, ri_janus->scale_pos, ri_janus->scale_vel); 
+    to_double(r->particles, ri_janus->p_int, r->N, ri_janus->scale_pos, ri_janus->scale_vel); 
+    ri_janus->is_synchronized = 0;
 }
 void reb_integrator_janus_reset(struct reb_simulation* r){
     struct reb_simulation_integrator_janus* const ri_janus = &(r->ri_janus);

src/rebound.h

@@ -936,6 +936,14 @@ void reb_integrator_synchronize(struct reb_simulation* r);
  **/
 void reb_integrator_reset(struct reb_simulation* r);
 
+
+/** 
+ * @brief Manually convert current floating point number to integers for JANUS integrator. Mainly for testing.
+ * @param r The rebound simulation to be considered
+ **/
+void reb_integrator_janus_to_int(struct reb_simulation* r);
+
+
 /**
  * @brief Configure the boundary/root box
  * @details This function helps to setup the variables for the simulation box.