Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Merge branch '4-write-a-basic-test-simulation-and-user-guide-for-sphe…
Browse files Browse the repository at this point in the history 
…rical-harmonics-use' into oblate_coord_rot
  • Loading branch information
@daminton
daminton committed Mar 1, 2024
2 parents c8c0784 + df993d5 commit 7e9a9d4
Show file tree
Hide file tree
Showing 2 changed files with 124 additions and 41 deletions.
9 changes: 9 additions & 0 deletions docs/api.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -100,6 +100,15 @@ Initial Conditions Generation Functions
swiftest.init_cond.horizons_get_physical_properties
swiftest.init_cond.vec2xr

Gravitional Harmonics Functions
===============================

.. autosummary::
:toctree: generated/

swiftest.shgrav.clm_from_ellipsoid
swiftest.shgrav.clm_from_relief


Input/Output Processing Functions
==================================
Expand Down
156 changes: 115 additions & 41 deletions docs/user-guide/gravitational-harmonics/index.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,8 +15,18 @@ potential is given by the following equation:
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.

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()
Gravitational Harmonics coefficients
Gravitational Harmonics Coefficients
=====================================

Swiftest adopts the :math:`4\pi` or geodesy normalization for the gravitational harmonics coefficients described
Expand All @@ -30,16 +40,6 @@ The coefficients can be computed in a number of ways:

- 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
===============================================

Expand All @@ -59,54 +59,128 @@ manner. Here we use Chariklo as an example body and refer to Jacobi Ellipsoid mo
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}`).
Once the central body parameters are defined, we can compute the gravitational harmonics coefficients (:math:`C_{lm}`). Here we set the
reference radius flag to `True` and ask the function to return the reference radius. More in the additional capabilities section below.
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)
c_lm, cb_radius = swiftest.clm_from_ellipsoid(mass = cb_mass, density = cb_density, a = cb_a, b = cb_b, c = cb_c, lref_radius = True)

*Note: The maximum degree is set to 6 by default to ensure computational efficiency.*

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.
If the :math:`J2` and :math:`J4` terms are passed as well, Swiftest ignores them and uses the :math:`C_{lm}` terms instead.
Now the user can set up the rest of the simulation as usual. ::

# add other bodies and set simulation parameters
.
.
.

Final Steps for Running the Simulation
=======================================
# run the simulation
sim.run()

Add other bodies to the simulation. ::
Additional Capabilities of Swiftest's Coefficient Generator Functions
===========================================================================================

# Add user-defined massive bodies
npl = 5
density_pl = cb_density
The output from :func:`clm_from_ellipsoid <swiftest.shgrav.clm_from_ellipsoid>` and :func:`clm_from_relief <swiftest.shgrav.clm_from_relief>`
can be customised to the user's needs. Here we show some of the additional capabilities of these functions.

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)
Setting a reference radius for the coefficients
-------------------------------------------------

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)
The coefficients are computed with respect to a reference radius. `SHTOOLS <https://shtools.github.io/SHTOOLS/>`__ calculates it's own radius from
the axes passed, but there are difference ways to calculate the reference radius for non-spherical bodies in the literature. As a result, Swiftest allows
the user to explicitly set a reference radius (``ref_radius``) which scales the coefficients accordingly. This is particularly useful when a
specific radius is desired.

Set the parameters for the simulation and run. ::
This is done by setting ``lref_radius = True`` and passing a ``ref_radius``. Here we pass the Central Body radius (``cb_radius``) manually set earlier as
the reference. ::

sim.set_parameter(tstart=0.0, tstop=10.0, dt=0.01, istep_out=10, dump_cadence=0, compute_conservation_values=True, mtiny=mtiny)
c_lm, ref_radius = swiftest.clm_from_ellipsoid(mass = cb_mass, density = cb_density, a = cb_a, b = cb_b, c = cb_c, lref_radius = True, ref_radius = cb_radius)

# Run the simulation
sim.run()
When ``lref_radius == True``, it tells the function to return the reference radius used to calculate the
coefficients and look for any reference radius (``ref_radius``) passed. If no reference radius is passed, the function returns the radius calculated
internally. ::
c_lm, ref_radius = swiftest.clm_from_ellipsoid(mass = cb_mass, density = cb_density, a = cb_a, b = cb_b, c = cb_c, lref_radius = True)

Setting a reference radius for the coefficients
==================================================
By default, ``lref_radius`` is set to ``False``. In this case, the function only returns the coefficients. ::

c_lm = swiftest.clm_from_ellipsoid(mass = cb_mass, density = cb_density, a = cb_a, b = cb_b, c = cb_c)

We recommend extracting the ``ref_radius`` from the function output and using it accordingly.

Combinations of Principal Axes
-------------------------------

The user can pass any combinations of the principal axes (``a``, ``b``, and ``c``) with ``a`` being the only required one. This is particularly
useful for cases like oblate spheroids (:math:`a = b \neq c`). For example, the following statements are equivalent: ::
c_lm, ref_radius = swiftest.clm_from_ellipsoid(mass = cb_mass, density = cb_density, a = cb_a, b = cb_b, c = cb_c, lref_radius = True)

c_lm, ref_radius = swiftest.clm_from_ellipsoid(mass = cb_mass, density = cb_density, a = cb_a, b = cb_a, c = cb_c, lref_radius = True)

c_lm, ref_radius = swiftest.clm_from_ellipsoid(mass = cb_mass, density = cb_density, a = cb_a, c = cb_c, lref_radius = True)

For bodies with :math:`a \neq b = c`, the following statements are equivalent: ::
c_lm, ref_radius = swiftest.clm_from_ellipsoid(mass = cb_mass, density = cb_density, a = cb_a, b = cb_b, c = cb_c, lref_radius = True)

c_lm, ref_radius = swiftest.clm_from_ellipsoid(mass = cb_mass, density = cb_density, a = cb_a, b = cb_b, c = cb_b, lref_radius = True)

c_lm, ref_radius = swiftest.clm_from_ellipsoid(mass = cb_mass, density = cb_density, a = cb_a, b = cb_b, lref_radius = True)


Setting the Maximum Degree :math:`l`
-------------------------------------

The user can set the maximum degree :math:`l` for the coefficients. ::

lmax = 4
c_lm, ref_radius = swiftest.clm_from_ellipsoid(mass = cb_mass, density = cb_density, a = cb_a, b = cb_b, c = cb_c, lmax = lmax, lref_radius = True)

``lmax`` is by currently capped to 6 to ensure computational efficiency. This is derived from Jean's law by setting the
characteristic wavelength (:math:`\lambda`) of a harmonic degree (:math:`l`) to the radius (:math:`R`) of the body.

.. math::
\lambda = \frac{2\pi R}{\sqrt{l(l+1)}}
\lambda = R \Rightarrow l = 6
.. 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. ::
The coefficients can be computed with respect to a reference radius. This is useful when the user wants to explicitly set the reference radius. ::
.. sim.set_parameter(tstart=0.0, tstop=10.0, dt=0.01, istep_out=10, dump_cadence=0, compute_conservation_values=True, mtiny=mtiny)
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)
.. # Run the simulation
.. sim.run()
Expand Down

0 comments on commit 7e9a9d4

Please sign in to comment.