Merge branch 'master' into ari-clean

README.rst

@@ -1,7 +1,7 @@
 REBOUND - An open-source multi-purpose N-body code
 ==================================================
 
-.. image:: http://img.shields.io/badge/rebound-v2.16.0-green.svg?style=flat
+.. image:: http://img.shields.io/badge/rebound-v2.16.1-green.svg?style=flat
     :target: http://rebound.readthedocs.org
 .. image:: https://badge.fury.io/py/rebound.svg
     :target: https://badge.fury.io/py/rebound

doc/Makefile

@@ -16,6 +16,7 @@ I18NSPHINXOPTS  = $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) .
 
 .PHONY: help clean html dirhtml singlehtml pickle json htmlhelp qthelp devhelp epub latex latexpdf text man changes linkcheck doctest coverage gettext
 
+
 help:
 	@echo "Please use \`make <target>' where <target> is one of"
 	@echo "  html       to make standalone HTML files"

doc/conf.py

@@ -66,7 +66,10 @@
     try:
         os.chdir("ipython")
         for example in glob.glob("../../ipython_examples/*.ipynb"):
-            subprocess.check_output(["jupyter", "nbconvert", example, "--to", "rst"])
+            print("Trying file: ", example)
+            outp = subprocess.check_output(["cp", example, example[23:]])
+            outp = subprocess.check_output(["jupyter", "nbconvert", example[23:], "--to", "rst"])
+            print(outp)
     except:
         with open("ipython.rst","w") as fd:
             fd.write("Examples can be found on github\n")
@@ -132,7 +135,7 @@
 # The short X.Y version.
 version = '2.16'
 # The full version, including alpha/beta/rc tags.
-release = '2.16.0'
+release = '2.16.1'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

doc/requirements.txt

@@ -1,6 +1,6 @@
 breathe==4.0.0
 jupyter==1.0.0
-jupyter-client==4.0.0
+jupyter-client==4.1.1
 jupyter-console==4.0.0
-jupyter-core==4.0.3
+jupyter-core==4.0.6
 numpydoc==0.5

rebound/tests/test_collisions.py

@@ -0,0 +1,96 @@
+import rebound
+import unittest
+import math
+import numpy as np
+
+class TestCollisions(unittest.TestCase):
+
+    def test_tree_remove_both(self):
+        sim = rebound.Simulation()
+        boxsize = 50000.           
+        sim.configure_box(boxsize)
+        sim.integrator = "leapfrog"
+        sim.boundary   = "open"
+        sim.gravity    = "tree"
+        sim.collision  = "tree"
+        def cor_remove_both(r, c):
+            r.contents.collisions_Nlog += 1
+            return 3
+        sim.collision_resolve = cor_remove_both
+
+        while sim.N< 10:
+            sim.add(m=1., r=100., x=np.random.uniform(-20,20),
+                                y=np.random.uniform(-20,20),
+                                z=np.random.uniform(-20,20))
+        sim.dt = 0.001
+        with self.assertRaises(rebound.NoParticles):
+            sim.integrate(1000.)
+        self.assertEqual(sim.collisions_Nlog,5)
+
+    def test_direct_remove_both(self):
+        sim = rebound.Simulation()
+        boxsize = 50000.           
+        sim.configure_box(boxsize)
+        sim.integrator = "leapfrog"
+        sim.boundary   = "open"
+        sim.collision  = "direct"
+        def cor_remove_both(r, c):
+            r.contents.collisions_Nlog += 1
+            return 3
+        sim.collision_resolve = cor_remove_both
+
+        while sim.N< 10:
+            sim.add(m=1., r=100., x=np.random.uniform(-20,20),
+                                y=np.random.uniform(-20,20),
+                                z=np.random.uniform(-20,20))
+        sim.dt = 0.001
+        with self.assertRaises(rebound.NoParticles):
+            sim.integrate(1000.)
+        self.assertEqual(sim.collisions_Nlog,5)
+
+    def test_tree_remove_one(self):
+        sim = rebound.Simulation()
+        boxsize = 50000.           
+        sim.configure_box(boxsize)
+        sim.integrator = "leapfrog"
+        sim.boundary   = "open"
+        sim.gravity    = "tree"
+        sim.collision  = "tree"
+        def cor_remove_both(r, c):
+            r.contents.collisions_Nlog += 1
+            return np.random.randint(1,3)
+        sim.collision_resolve = cor_remove_both
+
+        while sim.N< 50:
+            sim.add(m=1., r=100., x=np.random.uniform(-2,2),
+                                y=np.random.uniform(-2,2),
+                                z=np.random.uniform(-2,2))
+        sim.dt = 0.001
+        sim.integrate(sim.dt)
+        sim.integrate(2.*sim.dt)
+        self.assertLess(sim.N,25)
+
+    def test_direct_remove_one(self):
+        sim = rebound.Simulation()
+        boxsize = 50000.           
+        sim.configure_box(boxsize)
+        sim.integrator = "leapfrog"
+        sim.boundary   = "open"
+        sim.collision  = "direct"
+        def cor_remove_both(r, c):
+            r.contents.collisions_Nlog += 1
+            return np.random.randint(1,3)
+        sim.collision_resolve = cor_remove_both
+
+        while sim.N< 50:
+            sim.add(m=1., r=100., x=np.random.uniform(-2,2),
+                                y=np.random.uniform(-2,2),
+                                z=np.random.uniform(-2,2))
+        sim.dt = 0.001
+        sim.integrate(sim.dt)
+        sim.integrate(2.*sim.dt)
+        self.assertLess(sim.N,25)
+
+
+if __name__ == "__main__":
+    unittest.main()

rebound/tests/test_integrator.py

@@ -63,6 +63,28 @@ def setUp(self):
     def tearDown(self):
         self.sim = None
 
+    def test_ias15_globaloff(self):
+        self.sim.integrator = "ias15"
+        self.sim.ri_ias15.epsilon_global = 0
+        jupyr = 11.86*2.*math.pi
+        e0 = self.sim.calculate_energy()
+        self.assertNotEqual(e0,0.)
+        self.sim.integrate(1e3*jupyr)
+        e1 = self.sim.calculate_energy()
+        self.assertLess(math.fabs((e0-e1)/e1),1e-14)
+
+    def test_ias15_small_initial_dt(self):
+        self.sim.integrator = "ias15"
+        jupyr = 11.86*2.*math.pi
+        self.sim.dt = jupyr*1e-7
+        e0 = self.sim.calculate_energy()
+        self.assertNotEqual(e0,0.)
+        self.sim.integrate(self.sim.dt*1.001)
+        self.sim.dt = jupyr
+        self.sim.integrate(1e3*jupyr)
+        e1 = self.sim.calculate_energy()
+        self.assertLess(math.fabs((e0-e1)/e1),1e-14)
+
     def test_ias15(self):
         self.sim.integrator = "ias15"
         jupyr = 11.86*2.*math.pi

rebound/tests/test_interruble_pool.py

@@ -15,13 +15,15 @@ def test_nopool(self):
         pool = InterruptiblePool(2)
         params = [1.,1.1]
         res = [runsim(params[0]), runsim(params[1])]
-        self.assertAlmostEqual(res,[0.9950041652780258,1.095870355119381],delta=1e-15)
+        self.assertAlmostEqual(res[0],0.9950041652780258,delta=1e-15)
+        self.assertAlmostEqual(res[1],1.095870355119381,delta=1e-15)
 
     def test_pool(self):
         pool = InterruptiblePool(2)
         params = [1.,1.1]
         res = pool.map(runsim,params)
-        self.assertAlmostEqual(res,[0.9950041652780258,1.095870355119381],delta=1e-15)
+        self.assertAlmostEqual(res[0],0.9950041652780258,delta=1e-15)
+        self.assertAlmostEqual(res[1],1.095870355119381,delta=1e-15)
 
 if __name__ == "__main__":
     unittest.main()

rebound/tests/test_shearingsheet.py

@@ -5,7 +5,7 @@
 
 class TestShearingSheet(unittest.TestCase):
 
-    def test_whfast_nosafemode(self):
+    def test_saturnsrings(self):
         sim = rebound.Simulation()
         OMEGA = 0.00013143527     # [1/s]
         sim.ri_sei.OMEGA = OMEGA

rebound/tests/test_simulation.py

@@ -88,7 +88,7 @@ def test_calculate_orbits(self):
     def test_com(self):
         self.sim.move_to_com()
         com = self.sim.calculate_com()
-        self.assertEqual(com.x, 0.)
+        self.assertAlmostEqual(com.x, 0., delta=1e-15)
 
     def test_init_megno(self):
         self.sim.init_megno()

setup.py

@@ -49,7 +49,7 @@
     long_description = f.read()
 
 setup(name='rebound',
-    version='2.16.0',
+    version='2.16.1',
     description='An open-source multi-purpose N-body code',
     long_description=long_description,
     url='http://github.com/hannorein/rebound',

src/collision.c

@@ -190,28 +190,17 @@ void reb_collision_search(struct reb_simulation* const r){
 
         struct reb_collision c = r->collisions[i];
         if (c.p1 != -1 && c.p2 != -1){
-
             // Resolve collision
             int outcome = resolve(r, c);
 
             // Remove particles
-            int shift_pos = 0;
-            int keepsorted = 1;
-            if (r->tree_root){
-                keepsorted = 0;
-            }
             if (outcome & 1){
                 // Remove p1
-                if (keepsorted==0){
-                    if (c.p2==r->N-1){
-                        shift_pos = -c.p2+c.p1;
-                    }
-                }else{
-                    if (c.p1<c.p2){
-                        shift_pos = -1;
-                    }
+                if (c.p2==r->N-1 && !(r->tree_root)){
+                    // Particles swapped
+                    c.p2 = c.p1;
                 }
-                reb_remove(r,c.p1,keepsorted);
+                reb_remove(r,c.p1,0);
                 // Check for pair
                 for (int j=i+1;j<collisions_N;j++){
                     struct reb_collision cp = r->collisions[j];
@@ -220,11 +209,17 @@ void reb_collision_search(struct reb_simulation* const r){
                         r->collisions[j].p2 = -1;
                         // Will be skipped.
                     }
+                    if (cp.p1==r->N){
+                        r->collisions[j].p1 = c.p1;
+                    }
+                    if (cp.p2==r->N){
+                        r->collisions[j].p2 = c.p1;
+                    }
                 }
             }
             if (outcome & 2){
                 // Remove p2
-                reb_remove(r,c.p2+shift_pos,keepsorted);
+                reb_remove(r,c.p2,0);
                 // Check for pair
                 for (int j=i+1;j<collisions_N;j++){
                     struct reb_collision cp = r->collisions[j];
@@ -233,6 +228,12 @@ void reb_collision_search(struct reb_simulation* const r){
                         r->collisions[j].p2 = -1;
                         // Will be skipped.
                     }
+                    if (cp.p1==r->N){
+                        r->collisions[j].p1 = c.p2;
+                    }
+                    if (cp.p2==r->N){
+                        r->collisions[j].p2 = c.p2;
+                    }
                 }
             }
         }