diff --git a/examples/Chambers2013/.gitignore b/examples/Chambers2013/.gitignore index 64af50e94..0ea510e82 100644 --- a/examples/Chambers2013/.gitignore +++ b/examples/Chambers2013/.gitignore @@ -1,5 +1,6 @@ * !.gitignore -!init_cond.py +!initial_conditions.py !scattermovie.py +!run_simulation.py !README.txt diff --git a/examples/Chambers2013/README.txt b/examples/Chambers2013/README.txt index aeb95fc1f..c281be7fc 100644 --- a/examples/Chambers2013/README.txt +++ b/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. \ No newline at end of file diff --git a/examples/Chambers2013/init_cond.py b/examples/Chambers2013/initial_conditions.py similarity index 77% rename from examples/Chambers2013/init_cond.py rename to examples/Chambers2013/initial_conditions.py index d5c9f9f19..596d6829c 100755 --- a/examples/Chambers2013/init_cond.py +++ b/examples/Chambers2013/initial_conditions.py @@ -32,7 +32,8 @@ 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 @@ -40,6 +41,7 @@ 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)) diff --git a/examples/Chambers2013/run_simulation.py b/examples/Chambers2013/run_simulation.py new file mode 100755 index 000000000..5ae2fb2a3 --- /dev/null +++ b/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") diff --git a/examples/Chambers2013/scattermovie.py b/examples/Chambers2013/scattermovie.py index 98b0645ef..e1eda643a 100755 --- a/examples/Chambers2013/scattermovie.py +++ b/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')