Cleaned up the Chambers2013 example and added a new run_simulation sc…

…ript that can generate a short test case.

examples/Chambers2013/.gitignore

@@ -1,5 +1,6 @@
 *
 !.gitignore
-!init_cond.py
+!initial_conditions.py
 !scattermovie.py
+!run_simulation.py
 !README.txt

examples/Chambers2013/README.txt

@@ -15,9 +15,10 @@ Date             : December 6, 2022
 
 Included in the Chambers2013 example directory are the following files:
 
-	- README.txt      : This file
-	- init_cond.py    : A Python Script that generates a set of initial conditions.
-	- scattermovie.py : A Python Script that processes data.nc and generates Chambers2013-aescatter.mp4 or Chambers2013-aiscatter.mp4
+	- README.txt        : This file
+	- init_cond.py      : A Python script that generates a set of initial conditions.
+	- run_simulation.py : A Python script that will run the simulation.
+	- scattermovie.py   : A Python script that processes data.nc and generates Chambers2013-aescatter.mp4 or Chambers2013-aiscatter.mp4
 
 This example is intended to be run with Swiftest SyMBA. For details on how to generate, run, and analyze this example, 
 see the Swiftest User Manual.

examples/Chambers2013/init_cond.py → examples/Chambers2013/initial_conditions.py

@@ -32,14 +32,16 @@
 import matplotlib.pyplot as plt
 
 # Initialize simulation object
-sim = swiftest.Simulation(compute_conservation_values=True,  rotation=True, init_cond_format="EL",collision_model="fraggle",encounter_save="none")
+sim = swiftest.Simulation(compute_conservation_values=True, rotation=True, init_cond_format="EL",collision_model="fraggle",encounter_save="none")
+sim.clean()
 
 # Add bodies described in Chambers (2013) Sec. 2.1, with the uniform spatial distribution and two bodies sizes (big and small)
 Nb = 14
 Ns = 140
 Mb = 2.8e-7 * 14 / Nb
 Ms = 2.8e-8 * 140 / Ns
 dens = 3000.0 * sim.KG2MU / sim.M2DU**3
+rot = 1e-6 / sim.TU2S # Use a small but non-zero value for the initial rotation state to prevent divide-by-zero errors in analysis
 
 mtiny = 1e-2 * Ms
 minimum_fragment_mass = 1e-5 * Ms
@@ -55,37 +57,33 @@ def a_profile(n_bodies):
     The region with a < 0.7 AU deviates from this law, declining linearly with decreasing distance until R = 0 at 0.3 AU. 
     The outer edge of the disk is 2 AU in all cases.
     """
-    def sample(r_inner, r_break, r_outer, slope1, slope2):
-       """ 
-       Define the functions to pull random semi-major axes from a distribution using a rejection sampling method
-       This defines a 2-slope model with a break at r_break
-       Based on (https://stackoverflow.com/questions/66874819/random-numbers-with-user-defined-continuous-probability-distribution)
-       """
-       while True:
-          x = rng.uniform(r_inner, r_outer, size=1)
+    def sample(r_inner, r_break, r_outer, slope1, slope2, size):
+        """ 
+        Define the functions to pull random semi-major axes from a distribution using a rejection sampling method
+        This defines a 2-slope model with a break at r_break
+        Based on (https://stackoverflow.com/questions/66874819/random-numbers-with-user-defined-continuous-probability-distribution)
+        """
+        y=np.ones([size])
+        pdf=np.zeros([size])
+        x=np.empty_like(y)
+        while np.any(y>pdf): 
+            x = np.where(y>pdf,rng.uniform(r_inner, r_outer, size=size),x)
 
-          # The proportionality factor A ensures that the PDF approaches the same value on either side of the break point
-          # Assumes the break point is the max of the PDF
-          if x < r_break:
-              slope = slope1 + 1
-              A = 1.0
-          else:
-              slope = slope2 + 1
-              A = r_break**(slope1-slope2) 
-          y = rng.uniform(0, A*r_break**slope, size=1)
-          pdf = A*x**(slope)
-          if (y < pdf):
-                return x
-
+            # The proportionality factor A ensures that the PDF approaches the same value on either side of the break point
+            # Assumes the break point is the max of the PDF
+            A=np.where(x < r_break,1.0,r_break**(slope1-slope2))
+            slope=np.where(x < r_break,slope1+1,slope2+1)
+            y = np.where(y>pdf,rng.uniform(0, A*r_break**slope, size=size),y)
+            pdf = np.where(y>pdf,A*x**(slope),pdf)
+        return x
+
     a_inner = 0.3
     a_break = 0.7
     a_outer = 2.0
     slope1 = 1.0
     slope2 = -1.5
-
-    a_vals = np.zeros(n_bodies)
-    for k in range(n_bodies):
-        a_vals[k] = sample(a_inner, a_break, a_outer, slope1, slope2)
+
+    a_vals = sample(a_inner, a_break, a_outer, slope1, slope2, n_bodies)
     return a_vals
 
 # Define the initial orbital elements of the big and small bodies
@@ -107,8 +105,8 @@ def sample(r_inner, r_break, r_outer, slope1, slope2):
 capmvals = rng.uniform(0.0, 360.0, Ns)
 Ipvalb = np.full((Nb,3), 0.4)
 Ipvals = np.full((Ns,3), 0.4)
-rotvalb = np.zeros_like(Ipvalb)
-rotvals = np.zeros_like(Ipvals)
+rotvalb = np.full_like(Ipvalb,rot)
+rotvals = np.full_like(Ipvals,rot)
 
 noise_digits = 4 # Approximately the number of digits of precision to vary the mass values to avoid duplicate masses
 epsilon = np.finfo(float).eps
@@ -132,8 +130,7 @@ def sample(r_inner, r_break, r_outer, slope1, slope2):
 sim.set_parameter(mtiny=mtiny, minimum_fragment_mass=minimum_fragment_mass, nfrag_reduction=nfrag_reduction)
 
 sim.set_parameter(tstop=3e8, dt=6.0875/365.25, istep_out=60000, dump_cadence=10)
-sim.clean()
-sim.write_param()
+sim.save()
 
 # Plot the initial conditions
 fig = plt.figure(figsize=(8.5, 11))

examples/Chambers2013/run_simulation.py

@@ -0,0 +1,33 @@
+#!/usr/bin/env python3
+
+"""
+ Copyright 2023 - David Minton
+ This file is part of Swiftest.
+ Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License 
+ as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
+ Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty 
+ of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
+ You should have received a copy of the GNU General Public License along with Swiftest. 
+ If not, see: https://www.gnu.org/licenses. 
+"""
+
+"""
+This will run the simulation from a set of initial conditions. The simulation parameters in this file are set to generate
+a very short simulation for testing purposes. Edit the values passed to the run() function as necessary.
+
+Input
+------
+simdata/param.in    : ASCII Swiftest parameter input file.
+
+Output
+------
+Outputs are stored in the /simdata subdirectory.
+
+"""
+import swiftest
+sim = swiftest.Simulation(read_param=True)
+
+# Original run parameters
+# sim.run(tstop=3e8, dt=6.0875/365.25, istep_out=60000, dump_cadence=10,integreator="symba")
+# 
+sim.run(tstop=10000.0, dt=6.0875/365.25, istep_out=10000, dump_cadence=0, integrator="symba")

examples/Chambers2013/scattermovie.py

@@ -15,6 +15,15 @@
 Creates a movie from a set of Swiftest output files. All simulation 
 outputs are stored in the /simdata subdirectory.
 
+**NOTE: You must have ffmpeg installed on your system before running this script. For instance, on MacOS:
+
+```brew install ffmpeg```
+
+on Ubuntu:
+
+```sudo apt-get install ffmpeg```
+
+
 Input
 ------
 param.in    : ASCII Swiftest parameter input file.
@@ -33,6 +42,8 @@
 from matplotlib import animation
 import matplotlib.colors as mcolors
 from collections import namedtuple
+
+
 plt.switch_backend('agg')
 
 titletext = "Chambers (2013)"
@@ -108,8 +119,8 @@ def init_plot(self):
         self.ax.set_ylabel(ylabel[plot_style], fontsize='16', labelpad=1) 
 
         leg = plt.legend(loc="upper left", scatterpoints=1, fontsize=10)
-        for i,l in enumerate(leg.legendHandles):
-            leg.legendHandles[i]._sizes = [20]
+        for i,l in enumerate(leg.legend_handles):
+            leg.legend_handles[i]._sizes = [20]
 
         if plot_style == "arotscatter":
             self.ax.set_yscale('log')