Development manual

andrade_mu_complex(ω, μ, η, α)

Return the complex shear modulus μ̃(ω) for Andrade rheology, using Eq 82b from Efroimsky, M. 2012. μ, η, τA can be scalar or arrays.

Note: see Efroimsky, M. 2012, Section 5.3 for detailed description. Note: J is the unrelaxed compliance which is inverse to the unrelaxed rigidity μ, The same is through in the frequency domain (i.e. with ar{J} = 1/ar{μ}) Appendix B1.

compute_darcy_displacement!(dis_rel, y, m, r, ω, ϕ, ηl, k, g, ρl)

Calculate the Darcy displacement vector at a particular radial level.

compute_displacement!(dis, y, m)

Calculate the solid displacement vector at a particular radial level.

compute_pore_pressure!(p, y, m)

Calculate the fluid pore pressur at a particular radial level.

compute_strain_ten!(ϵ, y, m, rr, ρr, gr, μr, Ksr, ω, ρlr, Klr, Kdr, αr, ηlr, ϕr, kr)

Calculate the strain tensor ϵ at a particular radial level.

compute_strain_ten!(ϵ, y, m, rr, ρr, gr, μr, Ksr, ω, ρlr, Klr, Kdr, αr, ηlr, ϕr, kr)

Calculate the strain tensor ϵ at a particular radial level.

compute_stress_ten!(σ, ϵ, y, m, rr, ρr, gr, μr, Ksr, ω, ρlr, Klr, Kdr, αr, ηlr, ϕr, kr)

Calculate the stress tensor σ at a particular radial level.

compute_stress_ten!(σ, ϵ, y, m, rr, ρr, gr, μr, Ksr, ω, ρlr, Klr, Kdr, αr, ηlr, ϕr, kr)

Calculate the stress tensor σ at a particular radial level.

   compute_y(r, ρ, g, μ, K, ω, ρₗ, Kl, Kd, α, ηₗ, ϕ, k, core="liquid")

Compute the two-phase tidal solution y with the propagator matrix method.

Returns y, an array with size (8, nr, nlay), where nr is the number of secondary layers, and nlay is the number of primary layers, e.g., core, mantle, crust. The first index corresponds to y1 to y8, where

y1 = radial displacement
y2 = tangential displacement
y3 = radial stress
y4 = tangential stress
y5 = gravitational potential
y6 = potential stress
y7 = pore pressure
y8 = radial Darcy displacement

Example

   # Define four layer body
   ρs = [7640, 3300, 3300, 3000.]  # Solid density
   μ =  [60, 60, 60, 60] .* 1e9    # Shear modulus
   κs = [100, 200, 200, 200].*1e9  # Solid bulk modulus
   η =  [1e25, 1e25, 1e14, 1e25]   # Shear viscosity
   ζ =  [1e25, 1e25, 1e15, 1e25]   # Compaction viscosity

   ρl = [0, 0, 3300, 0]            # Liquid density
   κl = [0, 0, 1e9, 0]             # Liquid bulk modulus
   κd = 0.01κs                     # Drained compaction modulus
   k =  [0, 0, 1e-7, 0]            # Permeability
   α =  1.0 .- κd ./ κs            # Biot's modulus
   ηl = [0, 0, 1.0, 0]             # Liquid viscosity
   ϕ =  [0, 0, 0.1, 0]             # Porosity

   ρ = (1 .- ϕ) .* ρs + ϕ .* ρl    # Bulk density

   # Discretize primary layers into 100 secondary layers
   rr = expand_layers(r, nr=100)   
   g = get_g(rr, ρ);               # get gravity profile 

   y1 = compute_y(rr, ρ, g, μ, κs, ω, ρl, κl, κd, α, ηl, ϕ, k);

compute_y(r, ρ, g, μ, K; core="liquid")

Compute the single-phase tidal solution y with the propagator matrix method.

Returns y, an array with size (6, nr, nlay), where nr is the number of secondary layers, and nlay is the number of primary layers, e.g., core, mantle, crust. The first index corresponds to y1 to y6, where

y1 = radial displacement
y2 = tangential displacement
y3 = radial stress
y4 = tangential stress
y5 = gravitational potential
y6 = potential stress

Example

# Define four layer body
ρs = [7640, 3300, 3300, 3000.]  # Solid density
μ =  [60, 60, 60, 60] .* 1e9    # Shear modulus
κs = [100, 200, 200, 200].*1e9  # Solid bulk modulus

# Discretize primary layers into 100 secondary layers
rr = expand_layers(r, nr=100)   
g = get_g(rr, ρs);               # get gravity profile 

y1 = compute_y(rr, ρ, g, μ, κs);

Convert input parameters into the required precision.

Convert input parameters into the required precision.

define_spherical_grid(res)

Create the spherical grid of angular resolution res in degrees. This is used for numerical integrations over solid angle. A new grid can easily be defined by recalling the function with a new res.

The grid is internal to Love, but can be accessed with

Love.clats[:] # colatitude grid
Love.lons[:]  # longitude grid

expand_layers(r; nr::Int=80)

Discretize the primary layers given by r into nr discrete secondary layers.

find_mush_index(ϕ)

Find the mush region using ϕ (melt fraction) profile.:

1. Identify all indices where ϕ ≠ 0.
2. Partition these into connected segments.
3. Pick the connected segment with the smallest starting index

(the shallowest mush layer).

1. Return the largest index in that segment (bottom of mush).

Returns nothing if ϕ is zero everywhere.

get_A!(A::Matrix, r, ρ, g, μ, K, ω, ρₗ, Kl, Kd, α, ηₗ, ϕ, k)

Compute the 8x8 A matrix in the ODE for the two-phase problem. These correspond to the coefficients given in Equation S4.6 in Hay et al., (2025).

See also get_A

get_A!(r, ρ, g, μ, K; λ=nothing)

Compute the 6x6 A matrix in the ODE for the solid-body problem. These correspond to the coefficients given in Equation S4.6 in Hay et al., (2025) when α=φ=0, as well as Sabadini and Vermeersen (2016) Eq. 1.95.

See also get_A

get_A(A::Matrix, r, ρ, g, μ, K, ω, ρₗ, Kl, Kd, α, ηₗ, ϕ, k)

Compute the 8x8 A matrix in the ODE for the two-phase problem. These correspond to the coefficients given in Equation S4.6 in Hay et al., (2025).

See also get_A!

get_A(r, ρ, g, μ, K)

Compute the 6x6 A matrix in the ODE for the solid-body problem.

See also get_A!

get_B!(B, r1, r2, g1, g2, ρ, μ, K, ω, ρₗ, Kl, Kd, α, ηₗ, ϕ, k)

Compute the 8x8 numerical integrator matrix, which integrates dy/dr from r1 to r2 for the two-phase problem.

B here represnts the RK4 integrator, given by Eq. S5.5 in Hay et al., (2025).

See also get_B

get_B!(B, r1, r2, g1, g2, ρ, μ, K)

Compute the 6x6 numerical integrator matrix, which integrates dy/dr from r1 to r2 for the solid-body problem.

B here represnts the RK4 integrator, given by Eq. S5.5 in Hay et al., (2025).

See also get_B

See 'get_B!' for definition.

See 'get_B!' for definition.

get_B_product!(Brod, r, ρ, g, μ, K, ω, ρₗ, Kl, Kd, α, ηₗ, ϕ, k)

Compute the product of the 8x8 B matrices within a primary layer. This is used to propgate the y solution across a single two-phase primary layer.

Bprod is denoted by D in Eq. S5.14 in Hay et al., (2025).

get_B_product!(Brod, r, ρ, g, μ, K)

Compute the product of the 6x6 B matrices within a primary layer. This is used to propgate the y solution across one single-phase (solid) primary layer.

Bprod is denoted by D in Eq. S5.14 in Hay et al., (2025).

Get the core solution vector

get_g(r, ρ)

Compute the radial gravity structure associated with a density profile r at intervals given by r.

r must be be a 2D array, with index 1 representing the top radius of secondary layers, and index 2 representing the top radius of primary layers.

function get_heating_map(y, r, ρ, g, μ, Ks, ω, ρl, Kl, Kd, α, ηl, ϕ, k, ecc; vol=false)

Get a surface heating map for two-phase tides and eccentricity forcing, assuming synchronous rotation. Heating rate is computed with numerical integration using the solution y returned by compute_y, using Eq. 2.39a/b/c integrated in radius. The forcing magnitude is determined by orbital frequency ω and eccentricity ecc.

Returns the radially integrated heating maps in W/m^2 for dissipation due to shear, compaction, and Darcy flow. If vol=true, then additional volumentric heating arrays will also be returned in W/m^3.

get_heating_map(y, r, ρ, g, μ, κ, ω, ecc; vol=false)

Get a surface heating map for solid-body tides and eccentricity forcing, assuming synchronous rotation. Heating rate is computed with numerical integration using the solution y returned by compute_y, using Eq. 2.39a/b integrated in radius. The forcing magnitude is determined by orbital frequency ω and eccentricity ecc.

Returns the radially integrated heating maps in W/m^2 for dissipation due to shear and compaction. If vol=true, then additional volumentric heating arrays will also be returned in W/m^3.

get_heating_profile(y, r, ρ, g, μ, Ks, ω, ρl, Kl, Kd, α, ηl, ϕ, k, ecc; lay=nothing)

Get the radial volumetric heating for two-phase tides and eccentricity forcing, assuming synchronous rotation. Heating rate is computed with numerical integration using the solution y returned by compute_y, using Eq. 2.39a/b/c integrated over solid angle. The forcing magnitude is determined by orbital frequency ω and eccentricity ecc. The heating profile for a specific layer is specified with lay, otherwise all layers will be caclulated.

Returns the angular averaged volumetric heating profiles in W/m^3 for dissipation due to shear, compaction and Darcy flow, as well as the total power dissipated in each primary layer.

function get_heating_profile(y, r, ρ, g, μ, κ, ω, ecc; lay=nothing)

Get the radial volumetric heating for solid-body tides and eccentricity forcing, assuming synchronous rotation. Heating rate is computed with numerical integration using the solution y returned by compute_y, using Eq. 2.39a/b integrated over solid angle. The forcing magnitude is determined by orbital frequency ω and eccentricity ecc. The heating profile for a specific layer is specified with lay, otherwise all layers will be caclulated.

Returns the angular averaged volumetric heating profiles in W/m^3 for dissipation due to shear and compaction, as well as the power dissipated in each primary layer in W/m^3.

get_solution(y, n, m, r, ρ, g, μ, K, ω)

Expand the tidal solution from y at degree n and order m.

Prerequisite calls

Note that define_spherical_grid() must be called first.

Returns

A tuple of the complex fourier coefficients on a lat-lon grid for

1. displacement, u
2. solid strain, εs
3. total stress, σ

For tensors, the array dimensions represent (nlat, nlon, 6xtensor components, secondary lay., primary lay.) For vectors, the array dimensions represent (nlat, nlon, 3xvector components, secondary lay., primary lay.) For scalars, the array dimensions represent (nlat, nlon, secondary lay., primary lay.)

get_solution(y, n, m, r, ρ, g, μ, K, ω, ρₗ, Kl, Kd, α, ηₗ, ϕ, k)

Expand the tidal solution from y at degree n and order m.

Prerequisite calls

Note that define_spherical_grid() must be called first.

Returns

A tuple of the complex fourier coefficients on a lat-lon grid for

1. displacement, u
2. solid strain, εs
3. total stress, σ
4. pore pressure, pl
5. darcy (relative) displacement, u_rel

For tensors, the array dimensions represent (nlat, nlon, 6xtensor components, secondary lay., primary lay.) For vectors, the array dimensions represent (nlat, nlon, 3xvector components, secondary lay., primary lay.) For scalars, the array dimensions represent (nlat, nlon, secondary lay., primary lay.)

To recover the real solution at time t, the forcing magnitude and direction must be known.

Example

sol_22 = get_solution(y1, 2, 2, rr, ρ, g, μ, κs, ω, ρl, κl, κd, α, ηl, ϕ, k);
sol_20  = get_solution(y1, 2,  0, rr, ρ, g, μ, κs, ω, ρl, κl, κd, α, ηl, ϕ, k);

# eccentricity forcing 
U22E =  7/8 * ω^2*R^2*ecc 
U22W = -1/8 * ω^2*R^2*ecc 
U20  = -3/2 * ω^2*R^2*ecc 

Y = Love.clats;
X = Love.lons;
P20 = 0.5 .* (3 .* cos.(Y).^2 .- 1)
P22 = 3 .* (1. .- cos.(Y).^2)

# define phase of solution (0 to 2pi)
ωt = 0.0

# Obtain the radial tide
# Combine the solutions - note that the west component is a complex conjugate
tide = (U22E*sol_22[1] + U20*sol_20[1])[:,:,1,end,end] 
        .+ U22W*conj.(sol_22[1])[:,:,1,end,end];
tide = 0.5real.(tide .* exp(-1im * ωt) .+ conj.(tide) .* exp(1im * ωt));

get_total_heating(y, ω, R, ecc)

Get the total heating rate for solution y due to orbital eccentricity e for a body with synchronous rotation rate ω and radius R.

set_G(new_G)

Overwrite the value of the Universal Gravitational Constant. Only to be used for non-dimensional calculations!

compute_fluid_lovenumbers(omega, R, H_magma, g, ρ_ratio, n, σ_R)

Calculate k2 lovenumbers in fluid.

Arguments

• omega::Float64 : Forcing frequency range.
• R::prec : Outer radius of fluid segment in mantle.
• H_magma::prec : Height of fluid segment in mantle.
• g::prec : Surface gravity at top of fluid segment in mantle.
• ρ_ratio::prec : Density contrast between current (fluid) and lower (non-fluid) layer.
• n::Int=2 : Power of the radial factor (goes with (r/a)^{n}, since r<<a only n=2 contributes significantly).
• sigma_R::Float64=1e-3 : Rayleigh drag coefficient.

Returns

• k2_T::precc : Complex Tidal k2 Lovenumber.
• k2_L::precc : Complex Load k2 Lovenumber.

heat_profile(power, ρ, r)

Construct 1D uniform heating profile.

Arguments

• power::prec : Total dissipated power.
• ρ::Array{prec,1} : Density profile of the planet.
• r::Array{prec,1} : Radial positions of layers, from bottom to top of segment.

Returns

• power_prf::Array{prec,1} : Heating profile.

mean_rho(ρ, r)

Calculate the mean density in segment.

Arguments

• ρ::Array{prec,1} : Density profile of the planet.
• r::Array{prec,1} : Radial positions of layers, from bottom to top of segment.

Returns

• mean_density::prec : Mean density in segment.

compute_solid_lovenumbers(omega, R, H_magma, g, ρ_ratio, n, σ_R)

Calculate k2 lovenumbers in solid.

Arguments

• omega::Float64 : Forcing frequency range.
• R::prec : Outer radius of fluid segment in mantle.
• H_magma::prec : Height of fluid segment in mantle.
• g::prec : Surface gravity at top of fluid segment in mantle.
• ρ_ratio::prec : Density contrast between current (fluid) and lower (non-fluid) layer.
• n::Int=2 : Power of the radial factor (goes with (r/a)^{n}, since r<<a only n=2 contributes significantly).
• sigma_R::Float64=1e-3 : Rayleigh drag coefficient.

Returns

• k2_T::precc : Complex Tidal k2 Lovenumber.
• k2_L::precc : Complex Load k2 Lovenumber.

mean_cmu(μc, r)

Calculate the Hill-averaged complex shear modulus in a spherical segment (Hill+1952).

Arguments

• μc::Array{precc,1} : Complex shear modulus profile of the planet.
• r::Array{prec,1} : Radial positions of layers, from bottom to top of segment.

Returns

• mean_μc::precc : Mean complex shear modulus in segment.

Notes

• (DOI 10.1088/0370-1298/65/5/307)

get_hansen(ecc)

Wrapper used in your original code: computes Hansen coefficients for given eccentricity. Relies on global n, m, kmin, kmax.

hansen_fft(n, m, e, kmin, kmax; N=nothing)

Compute Hansen coefficients X_k^{n,m}(e) using FFT on mean anomaly.

Returns a tuple (k, Xkm), where:

• k is a vector of integer k-indices
• Xkm are the corresponding Hansen coefficients (real part)

kepler_newton(M, e)

Solve Kepler’s equation E - e*sin(E) = M for the eccentric anomaly E, using Newton iteration. M may be scalar or array.

nextpow2_int(x)

Return the integer p such that 2^p ≥ x.