Merge remote-tracking branch 'origin/debug' into debug

examples/Basic_Simulation/initial_conditions.py

@@ -11,24 +11,37 @@
 
 #!/usr/bin/env python3
 """
-Generates a set of Swiftest input files from initial conditions.
+Generates and runs a set of Swiftest input files from initial conditions with the SyMBA integrator. All simulation 
+outputs are stored in the /simdata subdirectory.
 
-Returns
--------
-Updates sim.data with the simulation data
+Input
+------
+None.
+
+Output
+------
+bin.nc         : A NetCDF file containing the simulation output.
+dump_bin1.nc   : A NetCDF file containing the necessary inputs to restart a simulation from t!=0.
+dump_bin2.nc   : A NetCDF file containing the necessary inputs to restart a simulation from t!=0.
+dump_param1.in : An ASCII file containing the necessary parameters to restart a simulation.
+dump_param2.in : An ASCII file containing the necessary parameters to restart a simulation.
+fraggle.log    : An ASCII file containing the information of any collisional events that occured.
+init_cond.nc   : A NetCDF file containing the initial conditions for the simulation.
+param.in       : An ASCII file containing the parameters for the simulation.
+swiftest.log   : An ASCII file containing the information on the status of the simulation as it runs.
 """
 
 import swiftest
 import numpy as np
 from numpy.random import default_rng
 
-# Initialize the simulation object as a variable
+# Initialize the simulation object as a variable. Arguments may be defined here or through the sim.run() method.
 sim = swiftest.Simulation(fragmentation=True, minimum_fragment_mass = 2.5e-11, mtiny=2.5e-8)
 
-# Add the modern planets and the Sun using the JPL Horizons Database
+# Add the modern planets and the Sun using the JPL Horizons Database.
 sim.add_solar_system_body(["Sun","Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune","Pluto"])
 
-# Add 5 user-defined massive bodies
+# Add 5 user-defined massive bodies.
 npl         = 5
 density_pl  = 3000.0 / (sim.param['MU2KG'] / sim.param['DU2M'] ** 3)
 
@@ -47,7 +60,7 @@
 
 sim.add_body(name=name_pl, a=a_pl, e=e_pl, inc=inc_pl, capom=capom_pl, omega=omega_pl, capm=capm_pl, Gmass=GM_pl, radius=R_pl, rhill=Rh_pl, Ip=Ip_pl, rot=rot_pl)
 
-# Add 10 user-defined test particles
+# Add 10 user-defined test particles.
 ntp = 10
 
 name_tp     = ["TestParticle_01", "TestParticle_02", "TestParticle_03", "TestParticle_04", "TestParticle_05", "TestParticle_06", "TestParticle_07", "TestParticle_08", "TestParticle_09", "TestParticle_10"]
@@ -59,8 +72,8 @@
 capm_tp     = default_rng().uniform(0.0, 360.0, ntp)
 
 sim.add_body(name=name_tp, a=a_tp, e=e_tp, inc=inc_tp, capom=capom_tp, omega=omega_tp, capm=capm_tp)
-# Display the run configuration parameters
+# Display the run configuration parameters.
 sim.get_parameter()
 
-# Run the simulation
+# Run the simulation. Arguments may be defined here or thorugh the swiftest.Simulation() method.
 sim.run(tstart=0.0, tstop=1.0e3, dt=0.01, tstep_out=1.0e0, dump_cadence=0)

examples/Basic_Simulation/output_reader.py

@@ -14,24 +14,23 @@
 Reads and processes a Swiftest output file.
 
 Input
--------
-out.nc     : NetCDF file
-    Swiftest output file.
+------
+bin.nc     : A NetCDF file containing the simulation output.
 
-Returns
--------
-output.eps : Encapsulated PostScript file.
-    A figure containing the eccentricity vs. semi-major axis for all bodies at the start of the simulation.
+Output
+------
+output.eps : Encapsulated PostScript file depicting the eccentricity vs. semi-major axis for all bodies at the start 
+             of the simulation.
 """
 
 import swiftest
 import xarray as xr
 import matplotlib.pyplot as plt
 
-# Read in the simulation output and store it as an Xarray dataset
+# Read in the simulation output and store it as an Xarray dataset.
 sim = swiftest.Simulation(read_old_output_file=True)
 
-# Plot of the data and save the output plot
+# Plot of the data and save the output plot.
 colors = ['white' if x == 'Massive Body' else 'black' for x in sim.data['particle_type']]
 sizes = [100 if x == 'Massive Body' else 10 for x in sim.data['particle_type']]
 

examples/Fragmentation/swiftest_fragmentation.py

@@ -10,32 +10,74 @@
 """
 #!/usr/bin/env python3
 """
-Generates a set of Swiftest input files from initial conditions, runs Swiftest, and generates fragmentation movies.
-Returns
--------
-Updates sim.data with the simulation data, produces fragmentation movies.
+Generates and runs a set of Swiftest input files from initial conditions with the SyMBA integrator. All simulation 
+outputs for the disruption case are stored in the /disruption subdirectory. All simulation outputs for the hit and run 
+case are stored in the /hitandrun subdirectory. All simulation outputs for the super-catastrophic disruption case are 
+stored in the /supercat subdirectory.
+
+Input
+------
+None.
+
+Output
+------
+disruption/bin.nc         : A NetCDF file containing the simulation output.
+disruption/dump_bin1.nc   : A NetCDF file containing the necessary inputs to restart a simulation from t!=0.
+disruption/dump_bin2.nc   : A NetCDF file containing the necessary inputs to restart a simulation from t!=0.
+disruption/dump_param1.in : An ASCII file containing the necessary parameters to restart a simulation.
+disruption/dump_param2.in : An ASCII file containing the necessary parameters to restart a simulation.
+disruption/fraggle.log    : An ASCII file containing the information of any collisional events that occured.
+disruption/init_cond.nc   : A NetCDF file containing the initial conditions for the simulation.
+disruption/param.in       : An ASCII file containing the parameters for the simulation.
+disruption/swiftest.log   : An ASCII file containing the information on the status of the simulation as it runs.
+hitandrun/bin.nc          : A NetCDF file containing the simulation output.
+hitandrun/dump_bin1.nc    : A NetCDF file containing the necessary inputs to restart a simulation from t!=0.
+hitandrun/dump_bin2.nc    : A NetCDF file containing the necessary inputs to restart a simulation from t!=0.
+hitandrun/dump_param1.in  : An ASCII file containing the necessary parameters to restart a simulation.
+hitandrun/dump_param2.in  : An ASCII file containing the necessary parameters to restart a simulation.
+hitandrun/fraggle.log     : An ASCII file containing the information of any collisional events that occured.
+hitandrun/init_cond.nc    : A NetCDF file containing the initial conditions for the simulation.
+hitandrun/param.in        : An ASCII file containing the parameters for the simulation.
+hitandrun/swiftest.log    : An ASCII file containing the information on the status of the simulation as it runs.
+supercat/bin.nc           : A NetCDF file containing the simulation output.
+supercat/dump_bin1.nc     : A NetCDF file containing the necessary inputs to restart a simulation from t!=0.
+supercat/dump_bin2.nc     : A NetCDF file containing the necessary inputs to restart a simulation from t!=0.
+supercat/dump_param1.in   : An ASCII file containing the necessary parameters to restart a simulation.
+supercat/dump_param2.in   : An ASCII file containing the necessary parameters to restart a simulation.
+supercat/fraggle.log      : An ASCII file containing the information of any collisional events that occured.
+supercat/init_cond.nc     : A NetCDF file containing the initial conditions for the simulation.
+supercat/param.in         : An ASCII file containing the parameters for the simulation.
+supercat/swiftest.log     : An ASCII file containing the information on the status of the simulation as it runs.
 """
 import swiftest
 import numpy as np
 from numpy.random import default_rng
 
-sim_disruption = swiftest.Simulation(param_file="param.disruption.in", output_file_name="disruption.nc", tstart=0.0, tstop=1.0e-5, dt=1.0e-8, istep_out=1.0, fragmentation=True, minimum_fragment_gmass=1.0e-11, gmtiny=1.0e-11, output_file_format="XVEL", init_cond_format="XV")
+# Initialize the simulation object as a variable with arguments.
+sim_disruption = swiftest.Simulation(simdir="disruption", tstart=0.0, tstop=1.0e-5, dt=1.0e-8, istep_out=1.0, fragmentation=True, minimum_fragment_gmass=1.0e-11, gmtiny=1.0e-11, output_format="XVEL", init_cond_format="XV")
+# Add the Sun using the JPL Horizons Database.
 sim_disruption.add_solar_system_body(["Sun"])
-sim_disruption.add_body(name="Target", v1=1.0, v2=-1.807993e-05, v3=0.0, v4=-2.562596e-04, v5=6.280005, v6=0.0, idvals=1, Gmass=1e-7, radius=7e-6, rhill=9e-4, Ip1=0.4, Ip2=0.4, Ip3=0.4, rotx=0.0, roty=0.0, rotz=0.0)
-sim_disruption.add_body(name="Projectile", v1=1.0, v2=1.807993e-05, v3=0.0, v4=-2.562596e-04, v5=-6.280005, v6=0.0, idvals=2, Gmass=7e-10, radius=3.25e-6, rhill=4e-4, Ip1=0.4, Ip2=0.4, Ip3=0.4, rotx=0.0, roty=0.0, rotz=0.0)
+# Add a user-defined target body.
+sim_disruption.add_body(name="Target", rh=[[1.0, -1.807993e-05, 0.0]], vh=[[-2.562596e-04, 6.280005, 0.0]], Gmass=1e-7, radius=7e-6, rhill=9e-4, Ip=[[0.4, 0.4, 0.4]], rot=[[0.0, 0.0, 0.0]])
+# Add a user-defined projectile body.
+sim_disruption.add_body(name="Projectile", rh=[[1.0, 1.807993e-05, 0.0]], vh=[[-2.562596e-04, -6.280005, 0.0]], Gmass=7e-10, radius=3.25e-6, rhill=4e-4, Ip=[[0.4, 0.4, 0.4]], rot=[[0.0, 0.0, 0.0]])
+# Display the run configuration parameters.
 sim_disruption.get_parameter()
+# Run the simulation.
 sim_disruption.run()
 
-sim_hitandrun = swiftest.Simulation(param_file="param.hitandrun.in", output_file_name="hitandrun.nc", tstart=0.0, tstop=1.0e-5, dt=1.0e-8, istep_out=1.0, fragmentation=True, minimum_fragment_gmass=1.0e-11, gmtiny=1.0e-11, output_file_format="XVEL", init_cond_format="XV")
+# Do the same as above for the hit and run case.
+sim_hitandrun = swiftest.Simulation(simdir="hitandrun", tstart=0.0, tstop=1.0e-5, dt=1.0e-8, istep_out=1.0, fragmentation=True, minimum_fragment_gmass=1.0e-11, gmtiny=1.0e-11, output_format="XVEL", init_cond_format="XV")
 sim_hitandrun.add_solar_system_body(["Sun"])
-sim_hitandrun.add_body(name="Target", v1=1.0, v2=-4.2e-05, v3=0.0, v4=0.0, v5=6.28, v6=0.0, idvals=1, Gmass=1e-7, radius=7e-6, rhill=9e-4, Ip1=0.4, Ip2=0.4, Ip3=0.4, rotx=0.0, roty=0.0, rotz=6.0e4)
-sim_hitandrun.add_body(name="Projectile", v1=1.0, v2=4.2e-05, v3=0.0, v4=-1.5, v5=-6.28, v6=0.0, idvals=2, Gmass=7e-10, radius=3.25e-6, rhill=4e-4, Ip1=0.4, Ip2=0.4, Ip3=0.4, rotx=0.0, roty=0.0, rotz=1.0e5)
+sim_hitandrun.add_body(name="Target", rh=[[1.0, -4.2e-05, 0.0]], vh=[[0.0, 6.28, 0.0]], Gmass=1e-7, radius=7e-6, rhill=9e-4, Ip=[[0.4, 0.4, 0.4]], rot=[[0.0, 0.0, 6.0e4]])
+sim_hitandrun.add_body(name="Projectile", rh=[[1.0, 4.2e-05, 0.0]], vh=[[-1.5, -6.28, 0.0]], Gmass=7e-10, radius=3.25e-6, rhill=4e-4, Ip=[[0.4, 0.4, 0.4]], rot=[[0.0, 0.0, 1.0e5]])
 sim_hitandrun.get_parameter()
 sim_hitandrun.run()
 
-sim_supercat = swiftest.Simulation(param_file="param.supercat.in", output_file_name="supercat.nc", tstart=0.0, tstop=1.0e-5, dt=1.0e-8, istep_out=1.0, fragmentation=True, minimum_fragment_gmass=1.0e-11, gmtiny=1.0e-11, output_file_format="XVEL", init_cond_format="XV")
+# Do the same as above for the super-catastrophic disruption case.
+sim_supercat = swiftest.Simulation(simdir="supercat", tstart=0.0, tstop=1.0e-5, dt=1.0e-8, istep_out=1.0, fragmentation=True, minimum_fragment_gmass=1.0e-11, gmtiny=1.0e-11, output_format="XVEL", init_cond_format="XV")
 sim_supercat.add_solar_system_body(["Sun"])
-sim_supercat.add_body(name="Target", v1=1.0, v2=-4.2e-05, v3=0.0, v4=0.0, v5=6.28, v6=0.0, idvals=1, Gmass=1e-7, radius=7e-6, rhill=9e-4, Ip1=0.4, Ip2=0.4, Ip3=0.4, rotx=0.0, roty=0.0, rotz=-6.0e4)
-sim_supercat.add_body(name="Projectile", v1=1.0, v2=4.2e-05, v3=0.0, v4=1.0, v5=-6.28, v6=0.0, idvals=2, Gmass=1e-8, radius=3.25e-6, rhill=4e-4, Ip1=0.4, Ip2=0.4, Ip3=0.4, rotx=0.0, roty=0.0, rotz=1.0e5)
+sim_supercat.add_body(name="Target", rh=[[1.0, -4.2e-05, 0.0]], vh=[[0.0, 6.28, 0.0]], Gmass=1e-7, radius=7e-6, rhill=9e-4, Ip=[[0.4, 0.4, 0.4]], rot=[[0.0, 0.0, -6.0e4]])
+sim_supercat.add_body(name="Projectile", rh=[[1.0, 4.2e-05, 0.0]], vh=[[1.0, -6.28, 0.0]], Gmass=1e-8, radius=3.25e-6, rhill=4e-4, Ip=[[0.4, 0.4, 0.4]], rot=[[0.0, 0.0, 1.0e5]])
 sim_supercat.get_parameter()
 sim_supercat.run()

examples/helio_gr_test/helio_gr_test.py

@@ -1,28 +1,77 @@
-#!/usr/bin/env python
+"""
+ Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh
+ 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. 
+"""
+
+#!/usr/bin/env python3
+"""
+Generates and runs two sets of Swiftest input files from initial conditions with the helio integrator. All simulation 
+outputs for the general relativity run are stored in the /gr subdirectory while all simulation outputs for the run 
+without general reelativity are stored in the /nogr subdirectory.
+
+Input
+------
+None.
+
+Output
+------
+helio_gr_mercury_precession.png : Portable Network Graphic file depicting the precession of Mercury's perihelion over time
+                                  with data sourced from the JPL Horizons database, Swiftest run with general relativity, 
+                                  and Swiftest run without general relativity.
+gr/bin.nc                       : A NetCDF file containing the simulation output.
+gr/dump_bin1.nc                 : A NetCDF file containing the necessary inputs to restart a simulation from t!=0.
+gr/dump_bin2.nc                 : A NetCDF file containing the necessary inputs to restart a simulation from t!=0.
+gr/dump_param1.in               : An ASCII file containing the necessary parameters to restart a simulation.
+gr/dump_param2.in               : An ASCII file containing the necessary parameters to restart a simulation.
+gr/init_cond.nc                 : A NetCDF file containing the initial conditions for the simulation.
+gr/param.in                     : An ASCII file containing the parameters for the simulation.
+gr/swiftest.log                 : An ASCII file containing the information on the status of the simulation as it runs.
+nogr/bin.nc                     : A NetCDF file containing the simulation output.
+nogr/dump_bin1.nc               : A NetCDF file containing the necessary inputs to restart a simulation from t!=0.
+nogr/dump_bin2.nc               : A NetCDF file containing the necessary inputs to restart a simulation from t!=0.
+nogr/dump_param1.in             : An ASCII file containing the necessary parameters to restart a simulation.
+nogr/dump_param2.in             : An ASCII file containing the necessary parameters to restart a simulation.
+nogr/init_cond.nc               : A NetCDF file containing the initial conditions for the simulation.
+nogr/param.in                   : An ASCII file containing the parameters for the simulation.
+nogr/swiftest.log               : An ASCII file containing the information on the status of the simulation as it runs.
+"""
+
 import swiftest
 from astroquery.jplhorizons import Horizons
 import datetime
 import numpy as np
 import matplotlib.pyplot as plt
 
+# Initialize the simulation object as a variable. Define the directory in which the output will be placed.
 sim_gr = swiftest.Simulation(simdir="gr")
+# Add the modern planets and the Sun using the JPL Horizons Database.
 sim_gr.add_solar_system_body(["Sun","Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune"])
 
+# Initialize the simulation object as a variable. Define the directory in which the output will be placed.
 sim_nogr = swiftest.Simulation(simdir="nogr")
+# Add the modern planets and the Sun using the JPL Horizons Database.
 sim_nogr.add_solar_system_body(["Sun","Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune"])
 
+# Define a set of arguments that apply to both runs. For a list of possible arguments, see the User Manual.
 run_args = {"tstop":1000.0, "dt":0.005, "tstep_out":10.0, "dump_cadence": 0,"integrator":"helio"}
 
+# Run both simulations.
 sim_gr.run(**run_args,general_relativity=True)
 sim_nogr.run(**run_args,general_relativity=False)
 
-# Get the start and end date of the simulation so we can compare with the real solar system
+# Get the start and end date of the simulation so we can compare with the real solar system.
 start_date = sim_gr.ephemeris_date
 tstop_d = sim_gr.param['TSTOP'] * sim_gr.param['TU2S'] / swiftest.JD2S
 
 stop_date = (datetime.datetime.fromisoformat(start_date) + datetime.timedelta(days=tstop_d)).isoformat()
 
-#Get the ephemerides of Mercury for the same timeframe as the simulation
+#Get the ephemerides of Mercury for the same timeframe as the simulation.
 obj = Horizons(id='1', location='@sun',
                epochs={'start':start_date, 'stop':stop_date,
                        'step':'10y'})
@@ -38,6 +87,7 @@
 dvarpi_nogr = np.diff(varpisim_nogr)  * 3600 * 100 / run_args['tstep_out']
 dvarpi_obs = np.diff(varpi_obs) / np.diff(t) * 3600 * 100
 
+# Plot of the data and save the output plot.
 fig, ax = plt.subplots()
 
 ax.plot(t, varpi_obs, label="JPL Horizons",linewidth=2.5)
@@ -48,6 +98,7 @@
 ax.legend()
 plt.savefig("helio_gr_mercury_precession.png",dpi=300)
 
+# Print the data to the terminal.
 print('Mean precession rate for Mercury long. peri. (arcsec/100 y)')
 print(f'JPL Horizons         : {np.mean(dvarpi_obs)}')
 print(f'Swiftest No GR       : {np.mean(dvarpi_nogr)}')