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.1-green.svg?style=flat
+.. image:: http://img.shields.io/badge/rebound-v2.16.5-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/conf.py

@@ -135,7 +135,7 @@
 # The short X.Y version.
 version = '2.16'
 # The full version, including alpha/beta/rc tags.
-release = '2.16.1'
+release = '2.16.5'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

rebound/particle.py

@@ -421,85 +421,35 @@ def calculate_orbit(self, primary=None, G=None):
 
     def __sub__(self, other):
         if isinstance(other, Particle):
-            np = Particle()
-            np.x = self.x - other.x
-            np.y = self.y - other.y
-            np.z = self.z - other.z
-            np.vx = self.vx - other.vx
-            np.vy = self.vy - other.vy
-            np.vz = self.vz - other.vz
-            np.ax = self.ax - other.ax
-            np.ay = self.ay - other.ay
-            np.az = self.az - other.az
-            np.m = self.m - other.m
-            return np
+            clibrebound.reb_particle_minus.restype = Particle
+            return clibrebound.reb_particle_minus(self, other)
         return NotImplemented 
 
     def __add__(self, other):
         if isinstance(other, Particle):
-            np = Particle()
-            np.x = self.x + other.x
-            np.y = self.y + other.y
-            np.z = self.z + other.z
-            np.vx = self.vx + other.vx
-            np.vy = self.vy + other.vy
-            np.vz = self.vz + other.vz
-            np.ax = self.ax + other.ax
-            np.ay = self.ay + other.ay
-            np.az = self.az + other.az
-            np.m = self.m + other.m
-            return np
+            clibrebound.reb_particle_plus.restype = Particle
+            return clibrebound.reb_particle_plus(self, other)
         return NotImplemented 
 
     def __mul__(self, other):
         if isinstance(other, float):
-            np = Particle()
-            np.x  = other*self.x
-            np.y  = other*self.y
-            np.z  = other*self.z
-            np.vx = other*self.vx
-            np.vy = other*self.vy
-            np.vz = other*self.vz 
-            np.ax = other*self.ax
-            np.ay = other*self.ay
-            np.az = other*self.az 
-            np.m  = other*self.m 
-            return np
+            clibrebound.reb_particle_multiply.restype = Particle
+            return clibrebound.reb_particle_multiply(self, c_double(other))
         return NotImplemented 
 
     def __rmul__(self, other):
         if isinstance(other, float):
-            np = Particle()
-            np.x  = other*self.x
-            np.y  = other*self.y
-            np.z  = other*self.z
-            np.vx = other*self.vx
-            np.vy = other*self.vy
-            np.vz = other*self.vz 
-            np.ax = other*self.ax
-            np.ay = other*self.ay
-            np.az = other*self.az 
-            np.m  = other*self.m 
-            return np
+            clibrebound.reb_particle_multiply.restype = Particle
+            return clibrebound.reb_particle_multiply(self, c_double(other))
         return NotImplemented 
 
     def __truediv__(self, other):
         return self.__div__(other)
 
     def __div__(self, other):
         if isinstance(other, float):
-            np = Particle()
-            np.x  = self.x / other
-            np.y  = self.y / other
-            np.z  = self.z / other
-            np.vx = self.vx/ other
-            np.vy = self.vy/ other
-            np.vz = self.vz/ other 
-            np.ax = self.ax/ other
-            np.ay = self.ay/ other
-            np.az = self.az/ other 
-            np.m  = self.m / other 
-            return np
+            clibrebound.reb_particle_divide.restype = Particle
+            return clibrebound.reb_particle_divide(self, c_double(other))
         return NotImplemented 
 
     @property
@@ -587,3 +537,7 @@ def T(self):
     def orbit(self):
         return self.calculate_orbit()
 
+    @property
+    def jacobi_com(self):
+        clibrebound.reb_get_jacobi_com.restype = Particle
+        return clibrebound.reb_get_jacobi_com(byref(self))

rebound/simulation.py

@@ -868,37 +868,39 @@ def calculate_orbits(self, heliocentric=False, barycentric=False):
         return orbits
 
 # COM calculation 
-    def calculate_com(self, last=None):
+    def calculate_com(self, first=0, last=None):
         """
         Returns the center of momentum for all particles in the simulation.
 
         Parameters
         ----------
+        first: int, optional
+            If ``first`` is specified, only calculate the center of momentum starting
+            from index=``first``.
         last : int or None, optional
-            If ``last`` is specified only calculate the center of momentum for the
-            first ``last`` particles in the array (i.e., indices up to i-1, as used 
-            in Jacobi coordinates).
+            If ``last`` is specified only calculate the center of momentum up to 
+            (but excluding) index=``last``.  Same behavior as Python's range function.
 
         Examples
         --------
         >>> sim = rebound.Simulation()
+        >>> sim.add(m=1, x=-20)
+        >>> sim.add(m=1, x=-10)
         >>> sim.add(m=1, x=0)
-        >>> sim.add(m=1, x=1)
+        >>> sim.add(m=1, x=10)
+        >>> sim.add(m=1, x=20)
         >>> com = sim.calculate_com()
         >>> com.x
-        0.5
+        0.0 
+        >>> com = sim.calculate_com(first=2,last=4) # Considers indices 2,3
+        >>> com.x
+        5.0
 
         """
-        if last is not None:
-            last = min(last, self.N_real-1)
-            clibrebound.reb_get_jacobi_com.restype = Particle
-            com = clibrebound.reb_get_jacobi_com(byref(self.particles[last]))
-            return com
-        else:
-            clibrebound.reb_get_com.restype = Particle
-            com = clibrebound.reb_get_com(byref(self))
-            return com
-
+        if last is None:
+            last = self.N_real
+        clibrebound.reb_get_com_range.restype = Particle
+        return clibrebound.reb_get_com_range(byref(self), c_int(first), c_int(last))
 
 # Tools
     def move_to_com(self):
@@ -916,7 +918,15 @@ def calculate_energy(self):
         """
         clibrebound.reb_tools_energy.restype = c_double
         return clibrebound.reb_tools_energy(byref(self))
-
+
+    def calculate_angular_momentum(self):
+        """
+        Returns a list of the three (x,y,z) components of the total angular momentum of all particles in the simulation.
+        """
+        clibrebound.reb_tools_angular_momentum.restype = reb_vec3d
+        L = clibrebound.reb_tools_angular_momentum(byref(self))
+        return [L.x, L.y, L.z]
+
     def configure_box(self, boxsize, root_nx=1, root_ny=1, root_nz=1):
         """
         Initialize the simulation box.

rebound/tests/test_simulation.py

@@ -89,7 +89,33 @@ def test_com(self):
         self.sim.move_to_com()
         com = self.sim.calculate_com()
         self.assertAlmostEqual(com.x, 0., delta=1e-15)
-
+
+    def test_com_range(self):
+        sim = rebound.Simulation()
+        sim.add(m=1.)
+        sim.add(m=1., x=2.)
+        sim.add(m=2., x=5.)
+        com = sim.calculate_com(first=1)
+        self.assertAlmostEqual(com.x, 4., delta=1e-15)
+        com = sim.calculate_com(last=2)
+        self.assertAlmostEqual(com.x, 1., delta=1e-15)
+        com = sim.calculate_com(first=1,last=2)
+        self.assertAlmostEqual(com.x, 2., delta=1e-15)
+        com = sim.calculate_com(first=4, last=-3)
+        self.assertAlmostEqual(com.x, 0., delta=1e-15)
+
+    def test_jacobi_com(self):
+        sim = rebound.Simulation()
+        sim.add(m=1., x=1.)
+        sim.add(m=1., x=3.)
+        sim.add(m=2., x=5.)
+        com = sim.particles[1].jacobi_com
+        self.assertAlmostEqual(com.x, 1., delta=1e-15)
+        com = sim.particles[2].jacobi_com
+        self.assertAlmostEqual(com.x, 2., delta=1e-15)
+        com = sim.particles[0].jacobi_com
+        self.assertAlmostEqual(com.x, 0., delta=1e-15)
+
     def test_init_megno(self):
         self.sim.init_megno()
         self.assertEqual(self.sim.N,4)
@@ -102,6 +128,16 @@ def test_calculate_energy(self):
         energy = self.sim.calculate_energy()
         self.assertAlmostEqual(energy, -0.5e-3, delta=1e-14)
 
+    def test_calculate_angular_momentum(self):
+        sim = rebound.Simulation()
+        sim.add(m=1.)
+        sim.add(m=1.e-3, a=1., inc=0.3, Omega=0.5)
+        sim.add(m=1.e-3, a=3., inc=0.2, Omega = -0.8)
+        L0 = sim.calculate_angular_momentum()
+        sim.integrate(1.)
+        Lf = sim.calculate_angular_momentum()
+        for i in range(3):
+            self.assertAlmostEqual(abs((Lf[i]-L0[i])/L0[i]), 0., delta=1e-15)
 
     def test_additional_forces(self):
         def af(sim):

setup.py

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

src/particle.c

@@ -182,3 +182,63 @@ int reb_remove_by_id(struct reb_simulation* const r, int id, int keepSorted){
 	}
 	return success;
 }
+
+struct reb_particle reb_particle_minus(struct reb_particle p1, struct reb_particle p2){
+    struct reb_particle p = {0};
+    p.x = p1.x - p2.x;
+    p.y = p1.y - p2.y;
+    p.z = p1.z - p2.z;
+    p.vx = p1.vx - p2.vx;
+    p.vy = p1.vy - p2.vy;
+    p.vz = p1.vz - p2.vz;
+    p.ax = p1.ax - p2.ax;
+    p.ay = p1.ay - p2.ay;
+    p.az = p1.az - p2.az;
+    p.m = p1.m - p2.m;
+    return p;
+}
+
+struct reb_particle reb_particle_plus(struct reb_particle p1, struct reb_particle p2){
+    struct reb_particle p = {0};
+    p.x = p1.x + p2.x;
+    p.y = p1.y + p2.y;
+    p.z = p1.z + p2.z;
+    p.vx = p1.vx + p2.vx;
+    p.vy = p1.vy + p2.vy;
+    p.vz = p1.vz + p2.vz;
+    p.ax = p1.ax + p2.ax;
+    p.ay = p1.ay + p2.ay;
+    p.az = p1.az + p2.az;
+    p.m = p1.m + p2.m;
+    return p;
+}
+
+struct reb_particle reb_particle_multiply(struct reb_particle p1, double value){
+    struct reb_particle p = {0};
+    p.x = p1.x * value;
+    p.y = p1.y * value;
+    p.z = p1.z * value;
+    p.vx = p1.vx * value;
+    p.vy = p1.vy * value;
+    p.vz = p1.vz * value;
+    p.ax = p1.ax * value;
+    p.ay = p1.ay * value;
+    p.az = p1.az * value;
+    p.m = p1.m * value;
+    return p;
+}
+
+struct reb_particle reb_particle_divide(struct reb_particle p1, double value){
+    struct reb_particle p = {0};
+    p.x = p1.x / value;
+    p.y = p1.y / value;
+    p.z = p1.z / value;
+    p.vx = p1.vx / value;
+    p.vy = p1.vy / value;
+    p.vz = p1.vz / value;
+    p.ax = p1.ax / value;
+    p.ay = p1.ay / value;
+    p.az = p1.az / value;
+    p.m = p1.m / value;
+    return p;
+}

src/rebound.c

@@ -63,7 +63,7 @@
 static const char* logo[];              /**< Logo of rebound. */
 #endif // LIBREBOUND
 const char* reb_build_str = __DATE__ " " __TIME__;  // Date and time build string. 
-const char* reb_version_str = "2.16.1";         // **VERSIONLINE** This line gets updated automatically. Do not edit manually.
+const char* reb_version_str = "2.16.5";         // **VERSIONLINE** This line gets updated automatically. Do not edit manually.
 
 
 void reb_step(struct reb_simulation* const r){