Merge branch 'debug'

examples/symba_chambers_2013/init_cond.py

@@ -0,0 +1,91 @@
+#!/usr/bin/env python3
+import swiftest
+import numpy as np
+from numpy.random import default_rng
+
+# Initialize simulation object
+sim = swiftest.Simulation()
+
+# Set unit conversion factors
+MU2KG = swiftest.MSun
+TU2S = swiftest.YR2S
+DU2M = swiftest.AU2M
+GU = swiftest.GC / (DU2M**3 / (MU2KG * TU2S**2))
+sim.param['MU2KG'] = MU2KG
+sim.param['TU2S'] = TU2S
+sim.param['DU2M'] = DU2M
+
+# Simulation time parameters
+sim.param['T0'] = 0.0
+sim.param['TSTOP'] = 300e6
+sim.param['DT'] = 6 * swiftest.JD2S / sim.param['TU2S']
+t_print = 1000.0
+iout = int(np.ceil(t_print / sim.param['DT']))
+sim.param['ISTEP_OUT']  = iout
+sim.param['ISTEP_DUMP'] = iout
+
+# Optional output file names
+sim.param['PARTICLE_OUT'] = "particle.dat"
+sim.param['ENERGY'] = "YES"
+sim.param['ENERGY_OUT'] = "energy.dat"
+sim.param['PL_IN'] = "pl_chambers_2013.in"
+sim.param['CB_IN'] = "sun_MsunAUYR.in"
+
+# Simulation parameters
+sim.param['FRAGMENTATION'] = "YES"
+sim.param['ROTATION'] = "YES"
+sim.param['CHK_RMAX'] = 1000.0
+sim.param['CHK_EJECT'] = 1000.0
+sim.param['IN_FORM'] = 'EL'
+sim.param['OUT_FORM'] = 'EL'
+
+# Add central body
+sim.add("Sun")
+GMcb = sim.ds['GMass'].values[0]
+sim.add("Jupiter")
+sim.add("Saturn")
+
+# 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
+Ms = 2.8e-8
+dens = 3000.0 / (MU2KG / DU2M**3)
+Rb = (3 * Mb / (4 * np.pi * dens) )**(1.0 / 3.0)
+Rs = (3 * Ms / (4 * np.pi * dens) )**(1.0 / 3.0)
+sim.param['GMTINY'] = 1e-2 * GU * Ms
+sim.param['MIN_GMFRAG'] = 1e-4 * GU * Ms
+
+# Define the initial orbital elements of the big and small bodies
+avalb = default_rng().uniform(0.3, 2.0, Nb)
+avals = default_rng().uniform(0.3, 2.0, Ns)
+evalb = default_rng().uniform(0.0, 0.001, Nb)
+evals = default_rng().uniform(0.0, 0.001, Ns)
+incvalb = default_rng().uniform(0.0, 0.0005 * 180 / np.pi, Nb)
+incvals = default_rng().uniform(0.0, 0.0005 * 180 / np.pi, Ns)
+capomvalb = default_rng().uniform(0.0, 360.0, Nb)
+capomvals = default_rng().uniform(0.0, 360.0, Ns)
+omegavalb = default_rng().uniform(0.0, 360.0, Nb)
+omegavals = default_rng().uniform(0.0, 360.0, Ns)
+capmvalb = default_rng().uniform(0.0, 360.0, Nb)
+capmvals = default_rng().uniform(0.0, 360.0, Ns)
+GMvalb = np.full(Nb, Mb * GU)
+GMvals = np.full(Ns, Ms * GU)
+Rvalb = np.full(Nb, Rb)
+Rvals = np.full(Ns, Rs)
+Rhb = avalb * (GMvalb / (3 * GMcb))**(1.0/3.0)
+Rhs = avals * (GMvals / (3 * GMcb))**(1.0/3.0)
+
+# Give the bodies unique ids
+idb = np.arange(100, 100 + Nb)
+ids = np.arange(100 + Nb, 100 + Nb + Ns)
+
+# Populate the simulation object with the two types of bodies
+sim.addp(idb, avalb, evalb, incvalb, capomvalb, omegavalb, capmvalb, GMpl=GMvalb, Rpl=Rvalb, Rhill=Rhb)
+sim.addp(ids, avals, evals, incvals, capomvals, omegavals, capmvals, GMpl=GMvals, Rpl=Rvals, Rhill=Rhs)
+
+# Save everything to a set of initial conditions files
+sim.save("param.in")
+
+
+

examples/symba_chambers_2013/param.in

@@ -1,32 +1,38 @@
-!
-! Parameter file for Chambers 2013 in units of Msun, AU, year
-!
-T0             0.0e0
-TSTOP          1e8        ! simulation length in years
-DT             0.016              ! stepsize in years
-CB_IN          sun_MsunAUYR.in
-PL_IN          pl_chambers_2013.in
-TP_IN          tp.in
-IN_TYPE        ASCII
-ISTEP_OUT      6250              ! output cadence 
-ISTEP_DUMP     6250            ! system dump cadence
-BIN_OUT        bin.dat
-PARTICLE_OUT   particle.dat
-OUT_TYPE       REAL8              ! double precision real output
-OUT_FORM       EL                 ! osculating element output
-OUT_STAT       REPLACE
-CHK_CLOSE      yes                ! check for planetary close encounters
-CHK_RMAX       100000.0            ! discard outside of
-EXTRA_FORCE    no                 ! no extra user-defined forces
-BIG_DISCARD    no                 ! output all planets if anything discarded
-RHILL_PRESENT  yes                ! Hill's sphere radii in input file
-MU2KG          1.98847e30         ! (M_sun-> kg)
-DU2M           1.495979e11        ! distance unit to meters (AU --> m)
-TU2S           3.15569259747e7    ! time unit to seconds (years --> seconds)
-GMTINY         1e-10		  
-ENERGY         yes
-ENERGY_OUT     energy.dat
-ROTATION       yes
-FRAGMENTATION  yes
-DISCARD_OUT    discard.out
-SEED           8 12261555 871132 92734722 21132443 36344777 4334443 219291656 3848566
+! VERSION        Swiftest parameter input
+T0               0.0
+TSTOP            300000000.0
+DT               0.01642710472279261
+ISTEP_OUT        60876
+ISTEP_DUMP       60876
+OUT_FORM         EL
+OUT_TYPE         REAL8
+OUT_STAT         REPLACE
+IN_TYPE          ASCII
+PL_IN            pl_chambers_2013.in
+TP_IN            tp.in
+CB_IN            sun_MsunAUYR.in
+BIN_OUT          bin.dat
+CHK_QMIN         -1.0
+CHK_RMIN         -1.0
+CHK_RMAX         1000.0
+CHK_EJECT        1000.0
+CHK_QMIN_COORD   HELIO
+CHK_QMIN_RANGE   -1.0 -1.0
+MU2KG            1.988409870698051e+30
+TU2S             31557600.0
+DU2M             149597870700.0
+IN_FORM          EL
+ENC_OUT          enc.dat
+EXTRA_FORCE      NO
+DISCARD_OUT      discard.out
+PARTICLE_OUT     particle.dat
+BIG_DISCARD      NO
+CHK_CLOSE        YES
+RHILL_PRESENT    YES
+FRAGMENTATION    YES
+ROTATION         YES
+TIDES            NO
+ENERGY           YES
+GR               NO
+GU               39.476926408897626
+ENERGY_OUT       energy.dat