Merge branch 'Simulation_API_improvements' into debug

examples/Basic_Simulation/initial_conditions.py

@@ -30,48 +30,14 @@
 from numpy.random import default_rng
 
 # Initialize the simulation object as a variable
-sim = swiftest.Simulation()
+sim = swiftest.Simulation(init_cond_file_type="ASCII")
+
+sim.set_simulation_time(tstart=0.0, tstop=10.0, dt=0.005, tstep_out=1.0)
 
 # Add parameter attributes to the simulation object
-sim.param['T0']                 = 0.0
-sim.param['TSTOP']              = 10
-sim.param['DT']                 = 0.005
-sim.param['ISTEP_OUT']          = 200
-sim.param['ISTEP_DUMP']         = 200
-sim.param['OUT_FORM']           = 'XVEL'
-sim.param['OUT_TYPE']           = 'NETCDF_DOUBLE'
-sim.param['OUT_STAT']           = 'REPLACE'
-sim.param['IN_FORM']            = 'EL'
-sim.param['IN_TYPE']            = 'ASCII'
-sim.param['PL_IN']              = 'pl.in'
-sim.param['TP_IN']              = 'tp.in'
-sim.param['CB_IN']              = 'cb.in'
-sim.param['BIN_OUT']            = 'out.nc'
-sim.param['CHK_QMIN']           = swiftest.RSun / swiftest.AU2M
-sim.param['CHK_RMIN']           = swiftest.RSun / swiftest.AU2M
-sim.param['CHK_RMAX']           = 1000.0
-sim.param['CHK_EJECT']          = 1000.0
-sim.param['CHK_QMIN_COORD']     = 'HELIO'
-sim.param['CHK_QMIN_RANGE']     = f'{swiftest.RSun / swiftest.AU2M} 1000.0'
-sim.param['MU2KG']              = swiftest.MSun
-sim.param['TU2S']               = swiftest.YR2S
-sim.param['DU2M']               = swiftest.AU2M
-sim.param['EXTRA_FORCE']        = 'NO'
-sim.param['BIG_DISCARD']        = 'NO'
-sim.param['CHK_CLOSE']          = 'YES'
-sim.param['GR']                 = 'YES'
-sim.param['INTERACTION_LOOPS']  = 'ADAPTIVE'
-sim.param['ENCOUNTER_CHECK']    = 'ADAPTIVE'
-sim.param['RHILL_PRESENT']      = 'YES'
-sim.param['FRAGMENTATION']      = 'YES'
-sim.param['ROTATION']           = 'YES'
-sim.param['ENERGY']             = 'YES'
 sim.param['GMTINY']             = 1e-6
 sim.param['MIN_GMFRAG']         = 1e-9
 
-# Set gravitational units of the system
-GU = swiftest.GC / (sim.param['DU2M'] ** 3 / (sim.param['MU2KG'] * sim.param['TU2S'] ** 2))
-
 # Add the modern planets and the Sun using the JPL Horizons Database
 sim.add("Sun", idval=0, date="2022-08-08")
 sim.add("Mercury", idval=1, date="2022-08-08")
@@ -95,9 +61,9 @@
 capom_pl    = default_rng().uniform(0.0, 360.0, npl)
 omega_pl    = default_rng().uniform(0.0, 360.0, npl)
 capm_pl     = default_rng().uniform(0.0, 360.0, npl)
-GM_pl       = (np.array([6e23, 8e23, 1e24, 3e24, 5e24]) / sim.param['MU2KG']) * GU
-R_pl        = np.full(npl, (3 * (GM_pl / GU) / (4 * np.pi * density_pl)) ** (1.0 / 3.0))
-Rh_pl       = a_pl * ((GM_pl) / (3 * GU)) ** (1.0 / 3.0)
+GM_pl       = (np.array([6e23, 8e23, 1e24, 3e24, 5e24]) / sim.param['MU2KG']) * sim.GU
+R_pl        = np.full(npl, (3 * (GM_pl / sim.GU) / (4 * np.pi * density_pl)) ** (1.0 / 3.0))
+Rh_pl       = a_pl * ((GM_pl) / (3 * sim.GU)) ** (1.0 / 3.0)
 Ip1_pl      = np.array([0.4, 0.4, 0.4, 0.4, 0.4])
 Ip2_pl      = np.array([0.4, 0.4, 0.4, 0.4, 0.4])
 Ip3_pl      = np.array([0.4, 0.4, 0.4, 0.4, 0.4])

python/swiftest/swiftest/init_cond.py

@@ -146,7 +146,7 @@ def solar_system_horizons(plname, idval, param, ephemerides_start_date, ds):
         Rpl = Rcb
         J2 = J2RP2
         J4 = J4RP4
-        if param['ROTATION'] == 'YES':
+        if param['ROTATION']:
             Ip1 = [Ipsun[0]]
             Ip2 = [Ipsun[1]]
             Ip3 = [Ipsun[2]]
@@ -202,19 +202,19 @@ def solar_system_horizons(plname, idval, param, ephemerides_start_date, ds):
         if ispl:
             GMpl = []
             GMpl.append(GMcb[0] / MSun_over_Mpl[plname])
-            if param['CHK_CLOSE'] == 'YES':
+            if param['CHK_CLOSE']:
                 Rpl = []
                 Rpl.append(planetradius[plname] * DCONV)
             else:
                 Rpl = None
 
             # Generate planet value vectors
-            if (param['RHILL_PRESENT'] == 'YES'):
+            if (param['RHILL_PRESENT']):
                 rhill = []
                 rhill.append(pldata[plname].elements()['a'][0] * DCONV * (3 * MSun_over_Mpl[plname]) ** (-THIRDLONG))
             else:
                 rhill = None
-            if (param['ROTATION'] == 'YES'):
+            if (param['ROTATION']):
                 Ip1 = []
                 Ip2 = []
                 Ip3 = []
@@ -304,7 +304,7 @@ def vec2xr(param, idvals, namevals, v1, v2, v3, v4, v5, v6, GMpl=None, Rpl=None,
     else:
         iscb = False
 
-    if param['ROTATION'] == 'YES':
+    if param['ROTATION']:
         if Ip1 is None:
             Ip1 = np.full_like(v1, 0.4)
         if Ip2 is None:
@@ -325,10 +325,10 @@ def vec2xr(param, idvals, namevals, v1, v2, v3, v4, v5, v6, GMpl=None, Rpl=None,
     else:
         ispl = False
 
-    if ispl and param['CHK_CLOSE'] == 'YES' and Rpl is None:
+    if ispl and param['CHK_CLOSE'] and Rpl is None:
         print("Massive bodies need a radius value.")
         return None
-    if ispl and rhill is None and param['RHILL_PRESENT'] == 'YES':
+    if ispl and rhill is None and param['RHILL_PRESENT']:
         print("rhill is required.")
         return None
 
@@ -342,19 +342,19 @@ def vec2xr(param, idvals, namevals, v1, v2, v3, v4, v5, v6, GMpl=None, Rpl=None,
     if iscb:
         label_float = clab.copy()
         vec_float = np.vstack([GMpl,Rpl,J2,J4])
-        if param['ROTATION'] == 'YES':
+        if param['ROTATION']:
             vec_float = np.vstack([vec_float, Ip1, Ip2, Ip3, rotx, roty, rotz])
         particle_type = "Central Body"
     else:
         vec_float = np.vstack([v1, v2, v3, v4, v5, v6])
         if ispl:
             label_float = plab.copy()
             vec_float = np.vstack([vec_float, GMpl])
-            if param['CHK_CLOSE'] == 'YES':
+            if param['CHK_CLOSE']:
                 vec_float = np.vstack([vec_float, Rpl])
-            if param['RHILL_PRESENT'] == 'YES':
+            if param['RHILL_PRESENT']:
                 vec_float = np.vstack([vec_float, rhill])
-            if param['ROTATION'] == 'YES':
+            if param['ROTATION']:
                 vec_float = np.vstack([vec_float, Ip1, Ip2, Ip3, rotx, roty, rotz])
             particle_type = np.repeat("Massive Body",idvals.size)
         else: