Merge branch '4-write-a-basic-test-simulation-and-user-guide-for-sphe…

…rical-harmonics-use' into oblate_coord_rot

docs/user-guide/basic-simulation/index.rst

@@ -79,6 +79,12 @@ file. The default value is 10, which means that Swiftest will store 10 outputs i
 Setting ``dump_cadence`` to 0 is a a special case that tells Swiftest to store *all* output in memory until the end of the 
 simulation. This is useful for short simulations, but can be memory intensive for long simulations. 
 
+.. note::
+  Changing the value of ``dump_cadence`` does not change the amount of data that is output by the end of the simulation. It only
+  changes how often the data is written to disk. Changing the value of ``tstep_out`` (or ``istep_out``) *does* change the amount of
+  data that is output by the end of the simulation. Intermediate steps between output steps are not saved to disk, and cannot be
+  recovered after the simulation has finished.
+
 The choice of what values to set for ``tstep_out`` (or ``istep_out``), ``nstep_out``, and ``dump_cadence`` depends on the particular
 simulation. Higher values of ``dump_cadence`` are typically useful for simulations with small numbers of bodies and small values
 of ```tstep_out`` where frequent writing to disk can severely impact performance. For simulations with large numbers of bodies and 
@@ -87,10 +93,10 @@ so that the memory usage does not become too large. The default value of ``dump_
 caes.
 
 We can set these simulation parameters using the :func:`set_parameter <swiftest.Simulation.set_parameter>` method. 
-Here we have a simulation that runs for 1 My a step size of 0.01 y. We will also save the system every 1000 y and wait until the end
+Here we have a simulation that runs for 100,000 y a step size of 0.01 y. We will also save the system every 1000 y and wait until the end
 of the simulation to write the simulation data to file using the ``dump_cadence=0`` argument::
 
-    sim.set_parameter(tstop=1.0e6, tstep_out=1e3, dt=0.01, dump_cadence=0)
+    sim.set_parameter(tstop=1.0e5, tstep_out=1e3, dt=0.01, dump_cadence=0)
 
 Once everything is set up, we call the :func:`run <swiftest.Simulation.run>` method to integrate the system forward in time::
 
@@ -99,7 +105,7 @@ Once everything is set up, we call the :func:`run <swiftest.Simulation.run>` met
 Swiftest is relatively flexible with arguments. You can pass the parameters in when initializing the simulation object, or even later when running.
 So the following are all equivalent::
 
-    sim = swiftest.Simulation(tstop=1.0e6, tstep_out=1e3, dt=0.01, dump_cadence=0)
+    sim = swiftest.Simulation(tstop=1.0e5, tstep_out=1e3, dt=0.01, dump_cadence=0)
     sim.add_solar_system_body(["Sun","Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune"])
     sim.run()
 
@@ -110,7 +116,7 @@ So the following are all equivalent::
 
     sim = swiftest.Simulation()
     sim.add_solar_system_body(["Sun","Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune"])
-    sim.run(tstop=1.0e6, tstep_out=1e3, dt=0.01, dump_cadence=0)
+    sim.run(tstop=1.0e5, tstep_out=1e3, dt=0.01, dump_cadence=0)
 
 .. note::
     Swiftest uses OpenMP parallelization to help speed up the integration, however the parallelization is most effective when there
@@ -131,7 +137,69 @@ Once a simulation has been run, its output data is stored in the ``./simdata`` d
 default name of ``data.nc``, which is a netCDF file. It is read in and stored as an 
 `Xarray Dataset <https://docs.xarray.dev/en/stable/>`__ object in the ``sim.data`` attribute.
 
+Here is an example of what the dataset looks like after the above simulation has been run::
+
+  In [5]: sim.data
+  Out[5]: 
+  <xarray.Dataset> Size: 229kB
+
+    Dimensions:          (time: 101, space: 3, name: 9)
+      Coordinates:
+      * time             (time) float64 808B 0.0 1e+03 2e+03 ... 9.9e+04 1e+05
+      * space            (space) <U1 12B 'x' 'y' 'z'
+      * name             (name) <U32 1kB 'Sun' 'Mercury' ... 'Uranus' 'Neptune'
+    Data variables:
+        id               (name) int64 72B 0 1 2 3 4 5 6 7 8
+        status           (time, name) int64 7kB 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0
+        npl              (time) int64 808B 8 8 8 8 8 8 8 8 8 8 ... 8 8 8 8 8 8 8 8 8
+        ntp              (time) int64 808B 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0
+        nplm             (time) int64 808B 8 8 8 8 8 8 8 8 8 8 ... 8 8 8 8 8 8 8 8 8
+        particle_type    (name) <U32 1kB 'Central Body' ... 'Massive Body'
+        rh               (time, name, space) float64 22kB nan nan ... -15.76 -0.4033
+        vh               (time, name, space) float64 22kB nan nan ... -0.03418
+        gr_pseudo_vh     (time, name, space) float64 22kB nan nan ... -0.03418
+        a                (time, name) float64 7kB nan 0.3871 0.7233 ... 19.16 30.17
+        e                (time, name) float64 7kB nan 0.2056 ... 0.03379 0.008614
+        inc              (time, name) float64 7kB nan 7.003 3.394 ... 1.306 1.884
+        capom            (time, name) float64 7kB nan 48.3 76.6 ... 150.3 124.7
+        omega            (time, name) float64 7kB nan 29.2 54.96 ... 135.7 310.1
+        capm             (time, name) float64 7kB nan 338.3 200.5 ... 212.3 254.7
+        varpi            (time, name) float64 7kB nan 77.5 131.6 ... 286.0 74.82
+        lam              (time, name) float64 7kB nan 55.84 332.0 ... 138.3 329.5
+        f                (time, name) float64 7kB nan 327.0 200.2 ... 210.3 253.8
+        cape             (time, name) float64 7kB nan 333.0 200.3 ... 211.3 254.3
+        Gmass            (time, name) float64 7kB 39.48 6.554e-06 ... 0.002033
+        mass             (time, name) float64 7kB 1.0 1.66e-07 ... 5.15e-05
+        rhill            (time, name) float64 7kB nan nan nan nan ... nan nan nan
+        radius           (time, name) float64 7kB 0.00465 1.631e-05 ... 0.0001646
+        origin_time      (name) float64 72B 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+        origin_type      (name) <U32 1kB 'Initial conditions' ... 'Initial condit...
+        origin_rh        (name, space) float64 216B 0.0 0.0 0.0 ... 2.045 -0.7287
+        origin_vh        (name, space) float64 216B 0.0 0.0 0.0 ... 1.149 -0.02168
+        collision_id     (name) int64 72B 0 0 0 0 0 0 0 0 0
+        discard_time     (name) float64 72B 0.0 1.798e+308 ... 1.798e+308 1.798e+308
+        discard_rh       (name, space) float64 216B 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
+        discard_vh       (name, space) float64 216B 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
+        discard_body_id  (name) int64 72B -2147483647 -2147483647 ... -2147483647
+        Ip               (time, name, space) float64 22kB 0.0 0.0 0.07 ... 0.0 0.23
+        rot              (time, name, space) float64 22kB 642.2 ... 1.721e+05
+        rotphase         (time) float64 808B 0.0 281.4 197.1 ... 176.4 126.8 344.4
+        j2rp2            (time) float64 808B 4.754e-12 4.754e-12 ... 4.754e-12
+        j4rp4            (time) float64 808B -2.247e-18 -2.247e-18 ... -2.247e-18
+
+
+As you can see, even in this very simple example, the dataset contains a large amount of information about the simulated system. 
+For details about the definitions of *variables*, *dimensions*, and *coordinates*, see the 
+`Terminology <https://docs.xarray.dev/en/stable/user-guide/terminology.html>`__. section of the Xarray documentation. Xarray 
+Datasets are very powerful and flexible, and can be used to analyze and visualize the simulation data in a variety of ways. 
+Here is an example where we can generate a simple plot of the semimajor axis vs. time history of all the planets in the system::
+
+  sim.data['a'].where(sim.data.particle_type != 'Central Body', drop=True).plot(x='time',hue='name')
+
+.. image:: ../../_static/basic_simulation_a_vs_t_plot.png
+
+This is just a simple example of what you can do with the simulation data. Xarray has a large number of built-in plotting and 
+data processing functions. For more information, see the `Xarray documentation <https://docs.xarray.dev/en/stable/>`__.
 
-.. toctree::
 ..    :maxdepth: 2
 ..    :hidden:

docs/user-guide/gravitational-harmonics/index.rst

@@ -0,0 +1,116 @@
+##########################
+Gravitational Harmonics
+##########################
+
+Here, we show how to use Swiftest's Gravitational Harmonics capabilities. This is based on ``/spherical_harmonics_cb`` 
+in ``swiftest/examples``. Swiftest uses `SHTOOLS <https://shtools.github.io/SHTOOLS/>`__ to compute the gravitational 
+harmonics coefficients for a given body and calculate it's associated acceleration kick. The conventions and formulae used 
+to calculate the additional kick are described `here <https://sseh.uchicago.edu/doc/Weiczorek_2015.pdf>`__. The gravitational
+potential is given by the following equation:
+
+.. math::
+
+    U(r) = \frac{GM}{r} \sum_{l=0}^{\infty} \sum_{m=-l}^{l} \left( \frac{R_0}{r} \right)^l C_{lm} Y_{lm} (\theta, \phi) \label{grav_pot}
+
+Gravitational potential :math:`U` at a point :math:`\vec{r}`; :math:`\theta` is the polar angle; :math:`\phi` is the azimuthal angle; 
+:math:`R_0` is the central body radius; :math:`G` is the gravitational constant; :math:`Y_{lm}` is the spherical harmonic function at 
+degree :math:`l` and order :math:`m`; :math:`C_{lm}` is the corresponding coefficient.
+
+Gravitational Harmonics coefficients
+=====================================
+
+Swiftest adopts the  :math:`4\pi` or geodesy normalization for the gravitational harmonics coefficients described 
+in `Weiczorek et al. (2015) <https://sseh.uchicago.edu/doc/Weiczorek_2015.pdf>`__. 
+
+The coefficients can be computed in a number of ways: 
+
+- Using the axes measurements of the body. (:func:`clm_from_ellipsoid <swiftest.shgrav.clm_from_ellipsoid>`)
+
+- Using a surface relief grid (:func:`clm_from_relief <swiftest.shgrav.clm_from_relief>`). *Note: This function is still in development and may not work as expected.*
+
+- Manually entering the coefficients when adding the central body. (:func:`add_body <swiftest.Simulation.add_body>`)
+
+Set up a Simulation
+====================
+
+Let's start with setting up the simulation object with units of `km`, `days`, and `kg`. ::
+    
+    import swiftest
+
+    sim = swiftest.Simulation(DU2M = 1e3, TU = 'd', MU = 'kg', integrator = 'symba')
+    sim.clean() 
+
+Computing coefficients from axes measurements
+===============================================
+
+Given the axes measurements of a body, the gravitational harmonics coefficients can be computed in a straightforward 
+manner. Here we use Chariklo as an example body and refer to Jacobi Ellipsoid model from 
+`Leiva et al. (2017) <https://iopscience.iop.org/article/10.3847/1538-3881/aa8956>`__ for the axes measurements. ::
+
+    # Define the central body parameters. 
+    cb_mass = 6.1e18
+    cb_radius = 123
+    cb_a = 157 
+    cb_b = 139 
+    cb_c = 86 
+    cb_volume = 4.0 / 3 * np.pi * cb_radius**3 
+    cb_density = cb_mass / cb_volume 
+    cb_T_rotation = 7.004 # hours
+    cb_T_rotation/= 24.0 # converting to julian days (TU)
+    cb_rot = [[0, 0, 360.0 / cb_T_rotation]] # degrees/TU
+
+Once the central body parameters are defined, we can compute the gravitational harmonics coefficients (:math:`C_{lm}`). 
+The output coefficients are already correctly normalized. ::
+
+    c_lm, cb_radius = swiftest.clm_from_ellipsoid(mass = cb_mass, density = cb_density, a = cb_a, b = cb_b, c = cb_c, lmax = 6, lref_radius = True)
+
+Add the central body to the simulation along with the coefficients. ::
+
+    sim.add_body(name = 'Chariklo', mass = cb_mass, rot = cb_rot, radius = cb_radius, c_lm = c_lm)
+
+Now the user can set up the rest of the simulation as they prefer.
+
+Final Steps for Running the Simulation
+=======================================
+
+Add other bodies to the simulation. ::
+
+    # Add user-defined massive bodies
+    npl         = 5
+    density_pl  = cb_density
+
+    name_pl     = ["SemiBody_01", "SemiBody_02", "SemiBody_03", "SemiBody_04", "SemiBody_05"]
+    a_pl        = rng.uniform(250, 400, npl)
+    e_pl        = rng.uniform(0.0, 0.05, npl)
+    inc_pl      = rng.uniform(0.0, 10, npl)
+    capom_pl    = rng.uniform(0.0, 360.0, npl)
+    omega_pl    = rng.uniform(0.0, 360.0, npl)
+    capm_pl     = rng.uniform(0.0, 360.0, npl)
+    R_pl        = np.array([0.5, 1.0, 1.2, 0.75, 0.8])
+    M_pl        = 4.0 / 3 * np.pi * R_pl**3 * density_pl
+    Ip_pl       = np.full((npl,3),0.4,)
+    rot_pl      = np.zeros((npl,3))
+    mtiny       = 1.1 * np.max(M_pl)
+
+    sim.add_body(name=name_pl, a=a_pl, e=e_pl, inc=inc_pl, capom=capom_pl, omega=omega_pl, capm=capm_pl, mass=M_pl, radius=R_pl,  Ip=Ip_pl, rot=rot_pl)
+
+Set the parameters for the simulation and run. ::
+
+    sim.set_parameter(tstart=0.0, tstop=10.0, dt=0.01, istep_out=10, dump_cadence=0, compute_conservation_values=True, mtiny=mtiny)
+
+    # Run the simulation
+    sim.run()
+
+Setting a reference radius for the coefficients
+==================================================
+
+The coefficients can be computed with respect to a reference radius. This is useful when the user wants to explicitly set the reference radius. ::
+
+    c_lm, cb_radius = swiftest.clm_from_ellipsoid(mass = cb_mass, density = cb_density, a = cb_a, b = cb_b, c = cb_c, lmax = 6, lref_radius = True, ref_radius = cb_radius)
+
+
+
+
+.. .. toctree::
+..    :maxdepth: 2
+..    :hidden:

docs/user-guide/index.rst

@@ -8,12 +8,12 @@ In this user guide, you will find detailed descriptions and examples that descri
 
 - A more in-depth :doc:`Detailed Simulation Setup <detailed-simulation-setup/index>`
 
-- Understanding the :doc:`Spherical Harmonics capabilities <spherical-harmonics/index>`
+- Understanding the :doc:`Gravitational Harmonics capabilities <gravitational-harmonics/index>`
 
 .. toctree::
    :maxdepth: 2
    :hidden:
 
    Basic Simulation <basic-simulation/index>
    Detailed Simulation Setup <detailed-simulation-setup/index>
-   Spherical Harmonics <spherical-harmonics/index>
+   Gravitational Harmonics <gravitational-harmonics/index>