diff --git a/docs/api.rst b/docs/api.rst index 604bc0553..c14d1eafe 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -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 ================================== diff --git a/docs/user-guide/gravitational-harmonics/index.rst b/docs/user-guide/gravitational-harmonics/index.rst index f4452666c..5f1505677 100644 --- a/docs/user-guide/gravitational-harmonics/index.rst +++ b/docs/user-guide/gravitational-harmonics/index.rst @@ -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 @@ -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 `) -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 =============================================== @@ -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 ` and :func:`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 `__ 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()