Updated examples and documentation using new synatx without "calculate"

docs/simulationdiagnostics.md

@@ -6,12 +6,12 @@ You can calculate the total energy (kinetic plus potential energy) of a simulati
 === "C"
     ```c
     struct reb_simulation* r = reb_create_simulation();
-    double energy = reb_calculate_energy(r);
+    double energy = reb_tools_energy(r);
     ```
 === "Python"
     ```python
     sim = rebound.Simulation()
-    energy = sim.calculate_energy()
+    energy = sim.energy()
     ```
 
 ## Angular momentum
@@ -20,12 +20,12 @@ You can calculate the angular momentum of a simulation using the following funct
 === "C"
     ```c
     struct reb_simulation* r = reb_create_simulation();
-    double angular_momentum = reb_calculate_angular_momentum(r);
+    struct reb_vec3d angular_momentum = reb_tools_angular_momentum(r);
     ```
 === "Python"
     ```python
     sim = rebound.Simulation()
-    Lx, Ly, Lz = sim.calculate_angular_momentum()
+    Lx, Ly, Lz = sim.angular_momentum()
     ```
 
 ## Center-of-mass

ipython_examples/EccentricComets.ipynb

@@ -119,7 +119,7 @@
    "outputs": [],
    "source": [
     "sim.move_to_com()\n",
-    "E0 = sim.calculate_energy()"
+    "E0 = sim.energy()"
    ]
   },
   {
@@ -172,7 +172,7 @@
    ],
    "source": [
     "sim.integrate(tmax)\n",
-    "dE = abs((sim.calculate_energy() - E0)/E0)\n",
+    "dE = abs((sim.energy() - E0)/E0)\n",
     "print(dE)"
    ]
   },

ipython_examples/EmbeddedOperatorSplittingMethods.ipynb

@@ -51,10 +51,10 @@
     "    \n",
     "    # First run to measure energy error\n",
     "    Emax = 0\n",
-    "    E0 = sim.calculate_energy()\n",
+    "    E0 = sim.energy()\n",
     "    while sim.t<tmax:\n",
     "        sim.integrate(sim.t+1.23456, exact_finish_time=0)\n",
-    "        E1 = sim.calculate_energy()\n",
+    "        E1 = sim.energy()\n",
     "        Emax = np.max([Emax,np.abs((E0-E1)/E0)])\n",
     "        \n",
     "    # Second run to measuring run time\n",

ipython_examples/HighOrderSymplectic.ipynb

@@ -82,13 +82,13 @@
     "    Nsamples = 1000\n",
     "    tmax = 2.*np.pi*1e3 # 1000 years\n",
     "    t_samples = tmax*np.sort(np.random.random(Nsamples))\n",
-    "    E0 = sim.calculate_energy() # initial energy\n",
+    "    E0 = sim.energy() # initial energy\n",
     "    Emax = 0. # maximum energy error\n",
     "    for t in t_samples:\n",
     "        # we do not want to change the timestep to reach t exactly, thus \n",
     "        # we need to set exact_finish_time=False and slighlty overshoot.\n",
     "        sim.integrate(t,exact_finish_time=False) \n",
-    "        E = sim.calculate_energy()\n",
+    "        E = sim.energy()\n",
     "        Emax = max(Emax, np.abs((E-E0)/E0))\n",
     "    return Emax"
    ]

ipython_examples/HybridIntegrationsWithMercurius.ipynb

@@ -52,7 +52,7 @@
     "sim.integrator = \"WHFast\" # This will end badly!\n",
     "sim.dt = sim.particles[1].P * 0.002 # Timestep a small fraction of innermost planet's period\n",
     "sim.move_to_com()\n",
-    "E0 = sim.calculate_energy() # Calculate initial energy \n",
+    "E0 = sim.energy() # Calculate initial energy \n",
     "rebound.OrbitPlot(sim)"
    ]
   },
@@ -78,7 +78,7 @@
    ],
    "source": [
     "sim.integrate(600*2.*math.pi)\n",
-    "E1 = sim.calculate_energy()\n",
+    "E1 = sim.energy()\n",
     "print(\"Relative energy error with WHFast: %f\"%((E0-E1)/E0))"
    ]
   },
@@ -110,9 +110,9 @@
     "sim.integrator = \"mercurius\" \n",
     "sim.dt = sim.particles[1].P * 0.002 # Timestep a small fraction of innermost planet's period\n",
     "sim.move_to_com()\n",
-    "E0 = sim.calculate_energy() # Calculate initial energy \n",
+    "E0 = sim.energy() # Calculate initial energy \n",
     "sim.integrate(600*2.*math.pi)\n",
-    "E1 = sim.calculate_energy()\n",
+    "E1 = sim.energy()\n",
     "print(\"Relative energy error with MERCURIUS: %e\"%((E1-E0)/E0))"
    ]
   },

ipython_examples/IntegratingArbitraryODEs.ipynb

@@ -147,7 +147,7 @@
     "    \n",
     "    r_nbody[i] = sim.particles[1].d\n",
     "    x_ho[i] = ode_ho.y[0]\n",
-    "    energies_nbody[i] = sim.calculate_energy()\n",
+    "    energies_nbody[i] = sim.energy()\n",
     "    energies_ho[i] = energy_ho(ode_ho)"
    ]
   },
@@ -339,7 +339,7 @@
     "    \n",
     "    r_nbody[i] = sim.particles[1].d\n",
     "    x_ho[i] = ode_ho.y[0]\n",
-    "    energies_nbody[i] = sim.calculate_energy()\n",
+    "    energies_nbody[i] = sim.energy()\n",
     "    energies_ho[i] = energy_ho(ode_ho)"
    ]
   },

ipython_examples/PrimordialEarth.ipynb

@@ -157,7 +157,7 @@
    "outputs": [],
    "source": [
     "sim.move_to_com()\n",
-    "E0 = sim.calculate_energy()"
+    "E0 = sim.energy()"
    ]
   },
   {
@@ -214,7 +214,7 @@
     "    sim.integrate(t,exact_finish_time=0)\n",
     "    totalN[i] = sim.N\n",
     "    encounterN[i] = sim.ri_mercurius._encounterN\n",
-    "    errors[i] = abs((sim.calculate_energy() - E0)/E0)"
+    "    errors[i] = abs((sim.energy() - E0)/E0)"
    ]
   },
   {

ipython_examples/Rotations.ipynb

@@ -345,7 +345,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "rot = rebound.Rotation.to_new_axes(newz=sim.calculate_angular_momentum())"
+    "rot = rebound.Rotation.to_new_axes(newz=sim.angular_momentum())"
    ]
   },
   {
@@ -403,9 +403,9 @@
     }
    ],
    "source": [
-    "print(sim.calculate_angular_momentum())\n",
+    "print(sim.angular_momentum())\n",
     "sim.rotate(rot)\n",
-    "print(sim.calculate_angular_momentum())"
+    "print(sim.angular_momentum())"
    ]
   },
   {

ipython_examples/SimulationArchive.ipynb

@@ -355,11 +355,11 @@
     "sa = rebound.SimulationArchive(filename)\n",
     "sim0 = sa[0]\n",
     "P = sim0.particles[1].P\n",
-    "E0 = sim.calculate_energy()\n",
+    "E0 = sim.energy()\n",
     "\n",
     "Eerr = np.zeros(Nout)\n",
     "for i, sim in enumerate(sa):\n",
-    "    E = sim.calculate_energy()\n",
+    "    E = sim.energy()\n",
     "    Eerr[i] = np.abs((E-E0)/E0)\n",
     "\n",
     "%matplotlib inline\n",

rebound/simulation.py

@@ -197,9 +197,9 @@ def to_new_axes(cls, newz, newx=None):
         >>> sim.add(m=1)                                                 
         >>> sim.add(m=1e-3, a=1, inc=0.3, Omega=4.2)
         >>> sim.add(m=1e-3, a=2, inc=0.1, Omega=0.5)    # Get a rotation to new axes where z points along total ang. momentum L. Didn't specify newx, 
-        >>> rot = rebound.Rotation.to_new_axes(newz=sim.calculate_angular_momentum()) # so it points along line of nodes between our original xy plane
+        >>> rot = rebound.Rotation.to_new_axes(newz=sim.angular_momentum()) # so it points along line of nodes between our original xy plane
         >>> sim.rotate(rot)                             # and the new plane perp. to L. Rotate our simulation into this new reference system
-        >>> print(sim.calculate_angular_momentum())     # Now the total angular momentum points in the z direction as we expect.
+        >>> print(sim.angular_momentum())     # Now the total angular momentum points in the z direction as we expect.
         """
         if not newx: # newx not specified, newx will point along z cross newz (line of nodes)
             clibrebound.reb_vec3d_cross.restype = Vec3d

rebound/tests/test_bs.py

@@ -37,9 +37,9 @@ def test_bs_outersolarsystem(self):
             sim.integrator = "bs"
             sim.ri_bs.eps_rel = eps
             sim.ri_bs.eps_abs = eps
-            e0 = sim.calculate_energy()
+            e0 = sim.energy()
             sim.integrate(1000)
-            e1 = sim.calculate_energy()
+            e1 = sim.energy()
             self.assertLess(math.fabs((e0-e1)/e1),5*eps)
 
     def test_bs_high_e(self):
@@ -51,9 +51,9 @@ def test_bs_high_e(self):
             sim.integrator = "bs"
             sim.ri_bs.eps_rel = eps
             sim.ri_bs.eps_abs = eps
-            e0 = sim.calculate_energy()
+            e0 = sim.energy()
             sim.integrate(1000)
-            e1 = sim.calculate_energy()
+            e1 = sim.energy()
             self.assertLess(math.fabs((e0-e1)/e1),60*eps)
 
     def test_bs_inout(self):

rebound/tests/test_eos.py

@@ -16,10 +16,10 @@ def tearDown(self):
     def _run(self):
         tmax = 100.
         Emax = 0
-        E0 = self.sim.calculate_energy()
+        E0 = self.sim.energy()
         while self.sim.t<tmax:
             self.sim.integrate(self.sim.t+1.23456, exact_finish_time=0)
-            E1 = self.sim.calculate_energy()
+            E1 = self.sim.energy()
             Emax = max([Emax,abs((E0-E1)/E0)])
         return Emax
 

rebound/tests/test_integrator.py

@@ -48,23 +48,23 @@ def test_whfast_hyperbolic(self):
         sim.add(m=1.)
         sim.add(m=1e-3, a=-1.,e=2.5)
         sim.integrator = "whfast"
-        x0 = sim.calculate_energy()
+        x0 = sim.energy()
         yr = -sim.particles[1].P
         sim.dt = 0.12*yr
         sim.integrate(1e2*yr)
-        x1 = sim.calculate_energy()
+        x1 = sim.energy()
         self.assertAlmostEqual(x0, x1, delta=1e-14)
 
     def test_whfast_verylargedt_hyperbolic(self):
         sim = rebound.Simulation()
         sim.add(m=1.)
         sim.add(m=1e-3, a=-1.,e=2.5)
         sim.integrator = "whfast"
-        x0 = sim.calculate_energy()
+        x0 = sim.energy()
         yr = -sim.particles[1].P
         sim.dt = 4.56*yr
         sim.integrate(1e3*yr)
-        x1 = sim.calculate_energy()
+        x1 = sim.energy()
         self.assertAlmostEqual(x0, x1, delta=1e-14)
 
 
@@ -81,86 +81,86 @@ 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()
+        e0 = self.sim.energy()
         self.assertNotEqual(e0,0.)
         self.sim.integrate(1e3*jupyr)
-        e1 = self.sim.calculate_energy()
+        e1 = self.sim.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()
+        e0 = self.sim.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()
+        e1 = self.sim.energy()
         self.assertLess(math.fabs((e0-e1)/e1),1e-14)
 
     def test_ias15(self):
         self.sim.integrator = "ias15"
         jupyr = 11.86*2.*math.pi
-        e0 = self.sim.calculate_energy()
+        e0 = self.sim.energy()
         self.assertNotEqual(e0,0.)
         self.sim.integrate(1e3*jupyr)
-        e1 = self.sim.calculate_energy()
+        e1 = self.sim.energy()
         self.assertLess(math.fabs((e0-e1)/e1),1e-14)
         # Longer tests:
         #self.sim.integrate(1e4*jupyr)
-        #e1 = self.sim.calculate_energy()
+        #e1 = self.sim.energy()
         #self.assertLess(math.fabs((e0-e1)/e1),10**13.5)
 
     def test_ias15_compensated(self):
         self.sim.integrator = "ias15"
         self.sim.gravity = "compensated"
         jupyr = 11.86*2.*math.pi
-        e0 = self.sim.calculate_energy()
+        e0 = self.sim.energy()
         self.assertNotEqual(e0,0.)
         self.sim.integrate(1e3*jupyr)
-        e1 = self.sim.calculate_energy()
+        e1 = self.sim.energy()
         self.assertLess(math.fabs((e0-e1)/e1),1e-14)
 
     def test_whfast_largedt(self):
         self.sim.integrator = "whfast"
         jupyr = 11.86*2.*math.pi
         self.sim.dt = 0.123*jupyr
-        e0 = self.sim.calculate_energy()
+        e0 = self.sim.energy()
         self.sim.integrate(1e3*jupyr)
         self.assertNotEqual(e0,0.)
-        e1 = self.sim.calculate_energy()
+        e1 = self.sim.energy()
         self.assertLess(math.fabs((e0-e1)/e1),1e-4)
 
     def test_whfast_smalldt(self):
         self.sim.integrator = "whfast"
         jupyr = 11.86*2.*math.pi
         self.sim.dt = 0.0123*jupyr
-        e0 = self.sim.calculate_energy()
+        e0 = self.sim.energy()
         self.sim.integrate(1e3*jupyr)
         self.assertNotEqual(e0,0.)
-        e1 = self.sim.calculate_energy()
+        e1 = self.sim.energy()
         self.assertLess(math.fabs((e0-e1)/e1),1e-6)
 
     def test_whfast_smalldt_compensated(self):
         self.sim.integrator = "whfast"
         self.sim.gravity = "compensated"
         jupyr = 11.86*2.*math.pi
         self.sim.dt = 0.0123*jupyr
-        e0 = self.sim.calculate_energy()
+        e0 = self.sim.energy()
         self.sim.integrate(1e3*jupyr)
         self.assertNotEqual(e0,0.)
-        e1 = self.sim.calculate_energy()
+        e1 = self.sim.energy()
         self.assertLess(math.fabs((e0-e1)/e1),1e-6)
 
     def test_whfast_verysmalldt(self):
         self.sim.integrator = "whfast"
         jupyr = 11.86*2.*math.pi
         self.sim.dt = 0.00123*jupyr
-        e0 = self.sim.calculate_energy()
+        e0 = self.sim.energy()
         self.sim.integrate(1e3*jupyr)
         self.assertNotEqual(e0,0.)
-        e1 = self.sim.calculate_energy()
+        e1 = self.sim.energy()
         self.assertLess(math.fabs((e0-e1)/e1),1e-8)
 
     def test_whfast_nosafemode(self):
@@ -169,10 +169,10 @@ def test_whfast_nosafemode(self):
         self.sim.ri_whfast.corrector = 11
         jupyr = 11.86*2.*math.pi
         self.sim.dt = 0.0123*jupyr
-        e0 = self.sim.calculate_energy()
+        e0 = self.sim.energy()
         self.sim.integrate(1e3*jupyr)
         self.assertNotEqual(e0,0.)
-        e1 = self.sim.calculate_energy()
+        e1 = self.sim.energy()
         self.assertLess(math.fabs((e0-e1)/e1),1e-9)
 
 if __name__ == "__main__":

rebound/tests/test_janus.py

@@ -17,9 +17,9 @@ def test_janus_energy(self):
             sim.ri_janus.order = o
             sim.ri_janus.scale_pos = 1e-16
             sim.ri_janus.scale_vel = 1e-16
-            e0 = sim.calculate_energy()
+            e0 = sim.energy()
             sim.integrate(1e2)
-            e1 = sim.calculate_energy()
+            e1 = sim.energy()
             self.assertLess(math.fabs((e0-e1)/e1),eps)
 
     def test_janus_reverse(self):