From 5870bfdf9005f3b1447716b9c4ecda543bc91e7d Mon Sep 17 00:00:00 2001 From: anand43 Date: Fri, 1 Mar 2024 13:30:38 -0500 Subject: [PATCH] Adding coefficients computation --- .../gravitational-harmonics/index.rst | 70 ++++++++++++++++--- 1 file changed, 60 insertions(+), 10 deletions(-) diff --git a/docs/user-guide/gravitational-harmonics/index.rst b/docs/user-guide/gravitational-harmonics/index.rst index ac926de91..285b26869 100644 --- a/docs/user-guide/gravitational-harmonics/index.rst +++ b/docs/user-guide/gravitational-harmonics/index.rst @@ -4,34 +4,84 @@ 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 `__ to compute the gravitational -harmonics coefficients for a given body and calculate it's associated acceleration kick. The additional acceleration -kick is based on the gravitaional potential described `here `__. +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 `__. 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. +: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. +Swiftest adopts the :math:`4\pi` or geodesy normalization for the gravitational harmonics coefficients described +in `Weiczorek et al. (2015) `__. + The coefficients can be computed in a number of ways: -- Using the axes measurements of the body ( :func:`clm_from_ellipsoid `) - -- Using a surface relief grid ( :func:`clm_from_relief `). +- Using the axes measurements of the body. (:func:`clm_from_ellipsoid `) +- Using a surface relief grid (:func:`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 `) - +- Manually entering the coefficients when adding the central body. (:func:`add_body `) +.. - Ensure to correctly normalize the coefficients if manually entering them. Computing coefficients from axes measurements =============================================== +Given the axes measurements of a body, the gravitational harmonics coefficients can be computed in a straightforward +manner. 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() + +Define the central body parameters. Here we use Chariklo as an example body and refer to Jacobi Ellipsoid model from +`Leiva et al. (2017) `__ for the axes measurements. :: + + 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. :: + + sim.add_body(name = 'Chariklo', mass = cb_mass, rot = cb_rot, radius = cb_radius, c_lm = c_lm) + +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: \ No newline at end of file