Development manual

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)

Ynm(n, m, theta, phi)

Compute the spherical harmonic Ynm for given n, m, theta, and phi.

Arguments

• n::Int : Tidal degree.
• m::Int : Tidal order.
• theta::Array{Float64,1} : Array of colatitudes in radians.
• phi::Array{Float64,1} : Array of longitudes in radians.

Returns

• Ynm::Array{ComplexF64,2} : 2D array of spherical harmonic values for each combination of theta and phi.

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

Compute the M matrix, which is used to propagate the solution across the entire interior. This is used in the compute_y function.

Arguments

• r::Array{prec,2} : 2D array of layer boundaries.
• ρ::Array{prec,1} : 1D array of layer densities.
• g::Array{prec,2} : 2D array of gravity values at the layer boundaries.
• μ::Array{prec,1} : 1D array of layer shear moduli.
• K::Array{prec,1} : 1D array of layer bulk moduli.
• n::Int : Tidal degree.

Keyword Arguments

• core::String="liquid" : Type of core, either "liquid" or "solid". This is used to compute the starting vector for the numerical integration across the interior.

Returns

• M::Array{precc,2} : 3x3 M matrix, which is used to propagate the solution across the entire interior.
• y1_4::Array{precc,4} : 4D array of the y solutions across each layer, which is used in the compute_y function to compute the solution vector across the interior.

compute_strain_ten!(ϵ, y, n, rr, ρr, gr, μr, Ksr)

Calculate the strain tensor ϵ at a particular radial level.

Arguments

• ϵ::Array{ComplexF64,3} : 3D array to store the strain tensor at a particular radial level, with dimensions corresponding to latitude, longitude, and the 6 independent components of the strain tensor.
• y::Array{precc,1} : 1D array of the solution vector y at a particular radial level, with 6 components.
• n::Int : Tidal degree.
• rr::prec : Radius at which to compute the strain tensor.
• ρr::prec : Density at radius rr.
• gr::prec : Gravity at radius rr.
• μr::prec : Shear modulus at radius rr.
• Ksr::prec : Bulk modulus at radius rr.

compute_y(r, g, M, R, y1_4, n; load=false)

Compute the solution vector y across the entire interior, given the M matrix and the y1_4 solutions across each layer. This is used to compute the strain tensor and heating profile.

Arguments

• r::Array{prec,2} : 2D array of layer boundaries.
• g::Array{prec,2} : 2D array of gravity values at the layer boundaries.
• M::Array{precc,2} : 3x3 M matrix, which is used to propagate the solution across the entire interior.
• R::prec : Surface radius of the body.
• y1_4::Array{precc,4} : 4D array of the y solutions across each layer, which is used in the compute_y function to compute the solution vector across the interior.
• n::Int : Tidal degree.

Keyword Arguments

• load::Bool=false : If true, compute the solution for a loaded problem.

Returns

• y::Array{ComplexF64,3} : 3D array of the solution vector y across the interior.

convert_params_to_prec(r, ρ, g, μ, κs)

Convert input parameters into the required precision.

Arguments

• r::Array{Float64,2} : 2D array of layer boundaries.
• ρ::Array{Float64,1} : 1D array of layer densities.
• g::Array{Float64,2} : 2D array of gravity values at the layer boundaries.
• μ::Array{Float64,1} : 1D array of layer shear moduli.
• κs::Array{Float64,1} : 1D array of layer bulk moduli.

Returns

Tuple of converted parameters in 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.

Arguments

• res::Float64 : Desired angular resolution in degrees.
• n::Int : Tidal degree.
• m::Int : Tidal order.

Notes

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

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

expand_layers(r; nr::Int=80)

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

Arguments

• r::Array{Float64,2} : 2D array of primary layer boundaries.

Keyword Arguments

• nr::Int=80 : Number of secondary layers to discretize.

Returns

• rs::Array{Float64,2} : 2D array of secondary layer boundaries/

get_A!(A, r, ρ, g, μ, K, n; λ=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.

Arguments

• A::Array{precc,2} : 6x6 A matrix at radius r, which is used in the ODE for the solid-body problem.
• r::prec : Radius at which to compute the A matrix.
• ρ::prec : Density at radius r.
• g::prec : Gravity at radius r.
• μ::prec : Shear modulus at radius r.
• K::prec : Bulk modulus at radius r.
• n::Int : Tidal degree.

Keyword Arguments

• λ::prec=nothing : Lamé's first parameter at radius r. If not provided, it is computed as λ = K - 2μ/3.

Notes

See also get_A

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

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

Arguments

• r::prec : Radius at which to compute the A matrix.
• ρ::prec : Density at radius r.
• g::prec : Gravity at radius r.
• μ::prec : Shear modulus at radius r.
• K::prec : Bulk modulus at radius r.
• n::Int : Tidal degree.

Returns

• A::Array{precc,2} : 6x6 A matrix at radius r, which is used in the ODE for the solid-body problem.

Notes

See also get_A!

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).

Arguments

• B::Array{precc,2} : 6x6 numerical integrator matrix for integrating dy/dr from r1 to r2 for the solid-body problem.
• r1::prec : Starting radius for integration.
• r2::prec : Ending radius for integration.
• g1::prec : Gravity at radius r1.
• g2::prec : Gravity at radius r2.
• ρ::prec : Density at radius r.
• μ::prec : Shear modulus at radius r.
• K::prec : Bulk modulus at radius r.
• n::Int : Tidal degree.

Notes

See also get_B

get_B(r1, r2, g1, g2, ρ, μ, K, n)

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

Arguments

• r1::prec : Starting radius for integration.
• r2::prec : Ending radius for integration.
• g1::prec : Gravity at radius r1.
• g2::prec : Gravity at radius r2.
• ρ::prec : Density at radius r.
• μ::prec : Shear modulus at radius r.
• K::prec : Bulk modulus at radius r.
• n::Int : Tidal degree.

Returns

• B::Array{precc,2} : 6x6 numerical integrator matrix for integrating dy/dr from r1 to r2 for the solid-body problem.

Notes

See 'get_B!' for definition.

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

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).

Arguments

• Bprod2::Array{precc,4} : 6x6x(nr-1)x(nlayers-1) array to store the B products across each secondary layer within each primary layer.
• r::Array{prec,2} : 2D array of layer boundaries.
• ρ::Array{prec,1} : 1D array of layer densities.
• g::Array{prec,2} : 2D array of gravity values at the layer boundaries.
• μ::Array{prec,1} : 1D array of layer shear moduli.
• K::Array{prec,1} : 1D array of layer bulk moduli.
• n::Int : Tidal degree.

get_Ic(r, ρ, g, μ, type, n; M=6, N=3)

Get the core solution vector.

Arguments

• r::prec : Radius of the core boundary.
• ρ::prec : Density of the core.
• g::prec : Gravity at the core boundary.
• μ::prec : Shear modulus of the core.
• type::String : Type of core, either "liquid" or "solid".
• n::Int : Tidal degree.

Keyword Arguments

• M::Int=6 : Number of rows in the Ic matrix. This should be 6 for the solid-body problem.
• N::Int=3 : Number of linearly independent solutions to compute. This should be 3 for the solid-body problem.

Returns

• Ic::Array{precc,2} : MxN array of linearly independent solutions at the core boundary. These are used as starting vectors for the numerical integration across the interior.

get_g(r, ρ)

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

Arguments

• r::Array{Float64,2} : 2D array of layer boundaries.
• ρ::Array{Float64,1} : 1D array of layer densities. The length of ρ must be equal to the number of columns in r.

Returns

• g::Array{Float64,2} : 2D array of gravity values at the layer boundaries. The dimensions of g are the same as r.

Notes

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_profile(y, r, ρ, g, μ, κ, n, ω; 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 heating profile for a specific layer is specified with lay, otherwise all layers will be caclulated.

Arguments

• y::Array{ComplexF64,4} : 4D array of the solution vector y across the interior, returned by compute_y.
• r::Array{Float64,2} : 2D array of layer boundaries.
• ρ::Array{Float64,1} : 1D array of layer densities.
• g::Array{Float64,2} : 2D array of gravity values at the layer boundaries.
• μ::Array{Float64,1} : 1D array of layer shear moduli.
• κ::Array{Float64,1} : 1D array of layer bulk moduli.
• n::Int : Tidal degree.
• ω::Float64 : Tidal frequency in radians per second.

Keyword Arguments

• lay::Int=nothing : If specified, compute the heating profile for only the layer corresponding to this index. Otherwise, compute for all layers.

Returns

• Eμ_tot::Array{Float64,1} : 1D array of total power dissipated in each primary layer due to shear, in W.
• Eμ_vol::Array{Float64,2} : 2D array of angular averaged volumetric heating profiles in W/m^3 for dissipation due to shear, with dimensions corresponding to the secondary layer and primary layer indices.
• Eκ_tot::Array{Float64,1} : 1D array of total power dissipated in each primary layer due to compaction, in W.
• Eκ_vol::Array{Float64,2} : 2D array of angular averaged volumetric heating profiles in W/m^3 for dissipation due to compaction, with dimensions corresponding to the secondary layer and primary layer indices.

Ynm(n, m, theta, phi)

Compute the spherical harmonic Ynm for given n, m, theta, and phi.

Arguments

• n::Int : Tidal degree.
• m::Int : Tidal order.
• theta::Array{Float64,1} : Array of colatitudes in radians.
• phi::Array{Float64,1} : Array of longitudes in radians.

Returns

• Ynm::Array{ComplexF64,2} : 2D array of spherical harmonic values for each combination of theta and phi.

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

Compute the 4x4 M matrix, which relates the solution at the surface and porous layer interface to the core solution.

Arguments

• r::Array{prec,2} : 2D array of layer boundaries.
• ρ::Array{prec,1} : 1D array of layer densities.
• g::Array{prec,2} : 2D array of gravity values at the layer boundaries.
• μ::Array{prec,1} : 1D array of layer shear moduli.
• K::Array{prec,1} : 1D array of layer bulk moduli.
• ω::prec : Forcing frequency.
• ρₗ::Array{prec,1} : 1D array of liquid densities at layer boundaries.
• Kl::Array{prec,1} : 1D array of liquid bulk moduli at layer boundaries.
• Kd::Array{prec,1} : 1D array of drained bulk moduli at layer boundaries.
• α::Array{prec,1} : 1D array of Biot coefficients at layer boundaries.
• ηₗ::Array{prec,1} : 1D array of liquid viscosities at layer boundaries.
• ϕ::Array{prec,1} : 1D array of porosities at layer boundaries.
• k::Array{prec,1} : 1D array of permeabilities at layer boundaries.
• n::Int : Tidal degree.

Keyword Arguments

• core::String="liquid" : Type of core, either "liquid" or "solid". This is used to compute the starting vector for the numerical integration across the interior.

Returns

• M::Array{precc,2} : 4x4 M matrix, which is used to propagate the solution across the entire interior.
• y1_4::Array{precc,4} : 4D array of the y solutions across each layer, which is used in the compute_y function to compute the solution vector across the interior.

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

Calculate the Darcy displacement vector at a particular radial level.

Arguments

• dis_rel::Array{ComplexF64,4} : 4D array to store the Darcy displacement vector at a particular radial level, with dimensions corresponding to latitude, longitude, and the 3 components of the Darcy displacement vector.
• y::Array{precc,1} : 1D array of the solution vector y at a particular radial level, with 8 components.
• n::Int : Tidal degree.
• r::prec : Radius at which to compute the Darcy displacement vector.
• ω::prec : Forcing frequency.
• ϕ::prec : Porosity at radius r.
• ηl::prec : Liquid viscosity at radius r.
• k::prec : Permeability at radius r.
• g::prec : Gravity at radius r.
• ρl::prec : Liquid density at radius r.

compute_pore_pressure!(p, y, n)

Calculate the fluid pore pressur at a particular radial level.

Arguments

• p::Array{ComplexF64,4} : 4D array to store the pore pressure at a particular radial level, with dimensions corresponding to latitude and longitude.
• y::Array{precc,1} : 1D array of the solution vector y at a particular radial level, with 8 components.
• n::Int : Tidal degree.

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

Calculate the strain tensor ϵ at a particular radial level.

Arguments

• ϵ::Array{ComplexF64,4} : 4D array to store the strain tensor at a particular radial level, with dimensions corresponding to latitude, longitude, and the 6 independent components of the strain tensor.
• y::Array{precc,1} : 1D array of the solution vector y at a particular radial level, with 6 components.
• n::Int : Tidal degree.
• rr::prec : Radius at which to compute the strain tensor.
• ρr::prec : Density at radius rr.
• gr::prec : Gravity at radius rr.
• μr::prec : Shear modulus at radius rr.
• Ksr::prec : Bulk modulus at radius rr.
• ω::prec : Forcing frequency.
• ρlr::prec : Liquid density at radius rr.
• Klr::prec : Liquid bulk modulus at radius rr.
• Kdr::prec : Drained bulk modulus at radius rr.
• αr::prec : Biot coefficient at radius rr.
• ηlr::prec : Liquid viscosity at radius rr.
• ϕr::prec : Porosity at radius rr.
• kr::prec : Permeability at radius rr.

compute_y(r, g, M, R, y1_4, n; load=false)

Compute the 8x1 solution vector at the surface and porous layer interface, which is used to compute the strain, Darcy flux, and pore pressure at a particular radial level. This is given by Eq. S5.20 in Hay et al., (2025).

Arguments

• r::Array{prec,2} : 2D array of layer boundaries.
• g::Array{prec,2} : 2D array of gravity values at the layer boundaries.
• M::Array{precc,2} : 4x4 M matrix, which is used to propagate the solution across the entire interior.
• R::prec : Surface radius of the body.
• y1_4::Array{precc,4} : 4D array of the y solutions across each layer, which is used in the compute_y function to compute the solution vector across the interior.
• n::Int : Tidal degree.

Keyword Arguments

• load::Bool=false : If true, compute the solution for a loaded problem.

Returns

• y::Array{ComplexF64,4} : 4D array of the solution vector y across the interior.

convert_params_to_prec(r, ρ, g, μ, κs, ω, ρl, κl, κd, α, ηl, ϕ, k)

Convert input parameters into the required precision.

Arguments

• r::Array{Float64,2} : 2D array of layer boundaries.
• ρ::Array{Float64,1} : 1D array of layer densities.
• g::Array{Float64,2} : 2D array of gravity values at the layer boundaries.
• μ::Array{Float64,1} : 1D array of layer shear moduli.
• κs::Array{Float64,1} : 1D array of layer bulk moduli.
• ω::Float64 : Forcing frequency.
• ρl::Array{Float64,1} : 1D array of layer liquid densities.
• κl::Array{Float64,1} : 1D array of layer liquid bulk moduli.
• κd::Array{Float64,1} : 1D array of layer drained bulk moduli.
• α::Array{Float64,1} : 1D array of layer Biot coefficients.
• ηl::Array{Float64,1} : 1D array of layer liquid viscosities.
• ϕ::Array{Float64,1} : 1D array of layer porosities.
• k::Array{Float64,1} : 1D array of layer permeabilities.

Returns

Tuple of converted parameters in 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.

Arguments

• res::Float64 : Desired angular resolution in degrees.
• n::Int : Tidal degree.
• m::Int : Tidal order.

Notes

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

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

expand_layers(r; nr::Int=80)

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

Arguments

• r::Array{Float64,2} : 2D array of primary layer boundaries.

Keyword Arguments

• nr::Int=80 : Number of secondary layers to discretize.

Returns

• rs::Array{Float64,2} : 2D array of secondary layer boundaries/

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

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).

Arguments

• A::Array{precc,2} : 6x6 A matrix at radius r, which is used in the ODE for the solid-body problem.
• r::prec : Radius at which to compute the A matrix.
• ρ::prec : Density at radius r.
• g::prec : Gravity at radius r.
• μ::prec : Shear modulus at radius r.
• K::prec : Bulk modulus at radius r.
• ω::prec : Forcing frequency.
• ρₗ::prec : Liquid density at radius r.
• Kl::prec : Liquid bulk modulus at radius r.
• Kd::prec : Drained bulk modulus at radius r.
• α::prec : Biot coefficient at radius r.
• ηₗ::prec : Liquid viscosity at radius r.
• ϕ::prec : Porosity at radius r.
• k::prec : Permeability at radius r.
• n::Int : Tidal degree.

Notes

See also get_A

get_A!(A, r, ρ, g, μ, K, n; λ=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.

Arguments

• A::Array{precc,2} : 6x6 A matrix at radius r, which is used in the ODE for the solid-body problem.
• r::prec : Radius at which to compute the A matrix.
• ρ::prec : Density at radius r.
• g::prec : Gravity at radius r.
• μ::prec : Shear modulus at radius r.
• K::prec : Bulk modulus at radius r.
• n::Int : Tidal degree.

Keyword Arguments

• λ::prec=nothing : Lamé's first parameter at radius r. If not provided, it is computed as λ = K - 2μ/3.

Notes

See also get_A

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

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).

Arguments

• r::prec : Radius at which to compute the A matrix.
• ρ::prec : Density at radius r.
• g::prec : Gravity at radius r.
• μ::prec : Shear modulus at radius r.
• K::prec : Bulk modulus at radius r.
• ω::prec : Forcing frequency.
• ρₗ::prec : Liquid density at radius r.
• Kl::prec : Liquid bulk modulus at radius r.
• Kd::prec : Drained bulk modulus at radius r.
• α::prec : Biot coefficient at radius r.
• ηₗ::prec : Liquid viscosity at radius r.
• ϕ::prec : Porosity at radius r.
• k::prec : Permeability at radius r.
• n::Int : Tidal degree.

Returns

• A::Array{precc,2} : 6x6 A matrix at radius r, which is used in the ODE for the solid-body problem.

See also get_A!

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

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

Arguments

• r::prec : Radius at which to compute the A matrix.
• ρ::prec : Density at radius r.
• g::prec : Gravity at radius r.
• μ::prec : Shear modulus at radius r.
• K::prec : Bulk modulus at radius r.
• n::Int : Tidal degree.

Returns

• A::Array{precc,2} : 6x6 A matrix at radius r, which is used in the ODE for the solid-body problem.

Notes

See also get_A!

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

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).

Arguments

• B::Array{precc,2} : 6x6 numerical integrator matrix for integrating dy/dr from r1 to r2 for the solid-body problem.
• r1::prec : Starting radius for integration.
• r2::prec : Ending radius for integration.
• g1::prec : Gravity at radius r1.
• g2::prec : Gravity at radius r2.
• ρ::prec : Density at radius r.
• μ::prec : Shear modulus at radius r.
• K::prec : Bulk modulus at radius r.
• ω::prec : Forcing frequency.
• ρₗ::prec : Liquid density at radius r.
• Kl::prec : Liquid bulk modulus at radius r.
• Kd::prec : Drained bulk modulus at radius r.
• α::prec : Biot coefficient at radius r.
• ηₗ::prec : Liquid viscosity at radius r.
• ϕ::prec : Porosity at radius r.
• k::prec : Permeability at radius r.
• n::Int : Tidal degree.

Notes

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).

Arguments

• B::Array{precc,2} : 6x6 numerical integrator matrix for integrating dy/dr from r1 to r2 for the solid-body problem.
• r1::prec : Starting radius for integration.
• r2::prec : Ending radius for integration.
• g1::prec : Gravity at radius r1.
• g2::prec : Gravity at radius r2.
• ρ::prec : Density at radius r.
• μ::prec : Shear modulus at radius r.
• K::prec : Bulk modulus at radius r.
• n::Int : Tidal degree.

Notes

See also get_B

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

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

Arguments

• r1::prec : Starting radius for integration.
• r2::prec : Ending radius for integration.
• g1::prec : Gravity at radius r1.
• g2::prec : Gravity at radius r2.
• ρ::prec : Density at radius r.
• μ::prec : Shear modulus at radius r.
• K::prec : Bulk modulus at radius r.
• ω::prec : Forcing frequency.
• ρₗ::prec : Liquid density at radius r.
• Kl::prec : Liquid bulk modulus at radius r.
• Kd::prec : Drained bulk modulus at radius r.
• α::prec : Biot coefficient at radius r.
• ηₗ::prec : Liquid viscosity at radius r.
• ϕ::prec : Porosity at radius r.
• k::prec : Permeability at radius r.
• n::Int : Tidal degree.

Returns

• B::Array{precc,2} : 8x8 numerical integrator matrix for integrating dy/dr from r1 to r2 for the two-phase problem.

Notes

See 'get_B!' for definition.

get_B(r1, r2, g1, g2, ρ, μ, K, n)

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

Arguments

• r1::prec : Starting radius for integration.
• r2::prec : Ending radius for integration.
• g1::prec : Gravity at radius r1.
• g2::prec : Gravity at radius r2.
• ρ::prec : Density at radius r.
• μ::prec : Shear modulus at radius r.
• K::prec : Bulk modulus at radius r.
• n::Int : Tidal degree.

Returns

• B::Array{precc,2} : 6x6 numerical integrator matrix for integrating dy/dr from r1 to r2 for the solid-body problem.

Notes

See 'get_B!' for definition.

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

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).

Arguments

• Bprod2::Array{precc,4} : 8x8x(nr-1)x(nlayers-1) array to store the B products across each secondary layer within each primary layer.
• r::Array{prec,2} : 2D array of layer boundaries.
• ρ::Array{prec,1} : 1D array of layer densities.
• g::Array{prec,2} : 2D array of gravity values at the layer boundaries.
• μ::Array{prec,1} : 1D array of layer shear moduli.
• K::Array{prec,1} : 1D array of layer bulk moduli.
• ω::prec : Forcing frequency.
• ρₗ::Array{prec,1} : 1D array of liquid densities at layer boundaries.
• Kl::Array{prec,1} : 1D array of liquid bulk moduli at layer boundaries.
• Kd::Array{prec,1} : 1D array of drained bulk moduli at layer boundaries.
• α::Array{prec,1} : 1D array of Biot coefficients at layer boundaries.
• ηₗ::Array{prec,1} : 1D array of liquid viscosities at layer boundaries.
• ϕ::Array{prec,1} : 1D array of porosities at layer boundaries.
• k::Array{prec,1} : 1D array of permeabilities at layer boundaries.
• n::Int : Tidal degree.

get_Ic(r, ρ, g, μ, type, n; M=8, N=4)

Get the core solution vector.

Arguments

• r::prec : Radius of the core boundary.
• ρ::prec : Density of the core.
• g::prec : Gravity at the core boundary.
• μ::prec : Shear modulus of the core.
• type::String : Type of core, either "liquid" or "solid".
• n::Int : Tidal degree.

Keyword Arguments

• M::Int=8 : Number of rows in the Ic matrix. This should be 8 for the two-phase problem.
• N::Int=4 : Number of linearly independent solutions to compute. This should be 4 for the two-phase problem.

Returns

• Ic::Array{precc,2} : MxN array of linearly independent solutions at the core boundary. These are used as starting vectors for the numerical integration across the interior.

get_g(r, ρ)

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

Arguments

• r::Array{Float64,2} : 2D array of layer boundaries.
• ρ::Array{Float64,1} : 1D array of layer densities. The length of ρ must be equal to the number of columns in r.

Returns

• g::Array{Float64,2} : 2D array of gravity values at the layer boundaries. The dimensions of g are the same as r.

Notes

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.

get_heating_profile(y, r, ρ, g, μ, Ks, ω, ρl, Kl, Kd, α, ηl, ϕ, k, n; 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 heating profile for a specific layer is specified with lay, otherwise all layers will be caclulated.

Arguments

• y::Array{ComplexF64,4} : 4D array of the solution vector y across the interior, returned by compute_y.
• r::Array{Float64,2} : 2D array of layer boundaries.
• ρ::Array{Float64,1} : 1D array of layer densities.
• g::Array{Float64,2} : 2D array of gravity values at the layer boundaries.
• μ::Array{Float64,1} : 1D array of layer shear moduli.
• Ks::Array{Float64,1} : 1D array of layer bulk moduli for shear dissipation.
• ω::Float64 : Tidal frequency in radians per second.
• ρl::Array{Float64,1} : 1D array of liquid densities at layer boundaries.
• Kl::Array{Float64,1} : 1D array of liquid bulk moduli at layer boundaries.
• Kd::Array{Float64,1} : 1D array of drained bulk moduli at layer boundaries.
• α::Array{Float64,1} : 1D array of Biot coefficients at layer boundaries.
• ηl::Array{Float64,1} : 1D array of liquid viscosities at layer boundaries.
• ϕ::Array{Float64,1} : 1D array of porosities at layer boundaries.
• k::Array{Float64,1} : 1D array of permeabilities at layer boundaries.
• n::Int : Tidal degree.

Keyword Arguments

• lay::Int=nothing : If specified, compute the heating profile for only the layer corresponding to this index. Otherwise, compute for all layers.

Returns

• Eμ_tot::Array{Float64,1} : 1D array of total power dissipated in each primary layer due to shear, in W.
• Eμ_vol::Array{Float64,2} : 2D array of angular averaged volumetric heating profiles in W/m^3 for dissipation due to shear, with dimensions corresponding to the secondary layer and primary layer indices.
• Eκ_tot::Array{Float64,1} : 1D array of total power dissipated in each primary layer due to compaction, in W.
• Eκ_vol::Array{Float64,2} : 2D array of angular averaged volumetric heating profiles in W/m^3 for dissipation due to compaction, with dimensions corresponding to the secondary layer and primary layer indices.
• El_tot::Array{Float64,1} : 1D array of total power dissipated in each primary layer due to Darcy flow, in W.
• El_vol::Array{Float64,2} : 2D array of angular averaged volumetric heating profiles in W/m^3 for dissipation due to Darcy flow, with dimensions corresponding to the secondary layer and primary layer indices.

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

Calculate k2 lovenumbers in fluid.

Arguments

• omega::prec : 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_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.

get_hansen(ecc)

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

get_k_range(ecc, n; tol=0.01)

Compute k-range for Hansen coefficients given eccentricity e and tidal degree n.

Arguments

• e::Float64 : Eccentricity of the orbit.
• n::Int : Tidal degree l.
• m::Int : Tidal order m.

Keyword Arguments

• tol::Float64=0.01 : Desired fractional contribution threshold for e^p.

Returns

• k_min::Int : Minimum k-index for Hansen coefficients.
• k_max::Int : Maximum k-index for Hansen coefficients.

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.