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