API Reference

abstract type AbstractAsteroidEphemerides

Supertype for all ephemerides types used in thermophysical simulations.

Abstract type for asteroid thermophysical simulation state.

abstract type AbstractBinaryAsteroidEphemerides <: AbstractAsteroidEphemerides

Supertype for binary-asteroid ephemerides.

Abstract type of a boundary condition for a heat conduction equation

abstract type AbstractSingleAsteroidEphemerides <: AbstractAsteroidEphemerides

Supertype for single-asteroid ephemerides.

Abstract type for thermophysical model solving algorithms.

Concrete subtypes are passed as the second argument to solve:

solve(problem, CrankNicolson(); kwargs...)

Abstract type for thermophysical problem definitions.

struct BinaryAsteroidEphemerides{R} <: AbstractBinaryAsteroidEphemerides

Ephemerides for a binary-asteroid thermophysical simulation.

The type parameter R controls whether force and torque are computed:

• R = Nothing : temperature only; R_primary_to_inertial = nothing
• R = Vector{SMatrix{3,3,Float64,9}} : force and torque are computed in the inertial frame

The secondary-to-inertial rotation is not stored but can be derived as:

\[R_{s2i} = R_{p2i} \cdot R_{p2s}^{\top}\]

Fields

• times : Simulation timesteps [s]
• r_sun : Sun position vector in the primary body-fixed frame at each timestep [m]
• r_secondary : Secondary position vector in the primary body-fixed frame at each timestep [m]
• R_primary_to_secondary : Passive rotation matrices from primary body-fixed to secondary body-fixed frame
• R_primary_to_inertial : Passive rotation matrices from primary body-fixed to inertial frame, or nothing when force/torque are not needed.

Constructors

BinaryAsteroidEphemerides(times, r_sun, r_secondary, R_primary_to_secondary)
    -> BinaryAsteroidEphemerides{Nothing}
BinaryAsteroidEphemerides(times, r_sun, r_secondary, R_primary_to_secondary, nothing)
    -> BinaryAsteroidEphemerides{Nothing}
BinaryAsteroidEphemerides(times, r_sun, r_secondary, R_primary_to_secondary, R_primary_to_inertial)
    -> BinaryAsteroidEphemerides{Vector{SMatrix{3,3,Float64,9}}}

An AbstractRange (e.g. range(et_begin, et_end; length=n)) may be passed as times and is automatically collected to Vector{Float64}. Plain AbstractVector / AbstractMatrix elements are automatically converted to SVector{3,Float64} / SMatrix{3,3,Float64,9} types in all fields.

struct BinaryAsteroidOutputSpec

Output specification for a binary-asteroid thermophysical simulation. Wraps two SingleAsteroidOutputSpec instances, one for each body.

Fields

• primary : Output spec for the primary body
• secondary : Output spec for the secondary body

Example

# Convenience constructor: shared Bool flags, separate output_times and subsurface_face_ids
output = BinaryAsteroidOutputSpec(
    output_times_primary,   subsurface_face_ids_primary,
    output_times_secondary, subsurface_face_ids_secondary;
    save_surface_temperature    = true,
    save_subsurface_temperature = true,
    save_face_forces            = false,
    save_forces                 = true,
    save_torques                = true,
)

# Base constructor: independent spec per body
output_primary   = SingleAsteroidOutputSpec(output_times_primary,   subsurface_face_ids_primary;   save_forces=true)
output_secondary = SingleAsteroidOutputSpec(output_times_secondary, subsurface_face_ids_secondary; save_forces=false)
output = BinaryAsteroidOutputSpec(output_primary, output_secondary)

struct BinaryAsteroidThermoPhysicalProblem

Defines the thermophysical problem for a binary asteroid system.

Fields

• primary : Problem definition for the primary body
• secondary : Problem definition for the secondary body
• with_mutual_shadowing : Whether to include mutual shadowing (eclipses)
• with_mutual_heating : Whether to include mutual heating

Usage

prob1 = SingleAsteroidThermoPhysicalProblem(shape1, thermo_params1; ...)
prob2 = SingleAsteroidThermoPhysicalProblem(shape2, thermo_params2; ...)
problem = BinaryAsteroidThermoPhysicalProblem(prob1, prob2;
    with_mutual_shadowing = true,
    with_mutual_heating   = true,
)
solution = solve(problem, CrankNicolson(); ephem = ephem, T₀ = 200.0)

BinaryAsteroidThermoPhysicalProblem(primary, secondary; kwargs...) -> problem

Construct a thermophysical problem for a binary asteroid system.

Arguments

• primary : Problem definition for the primary body
• secondary : Problem definition for the secondary body

Keyword Arguments

• with_mutual_shadowing = true : Whether to include mutual shadowing (eclipses)
• with_mutual_heating = true : Whether to include mutual heating

Notes

• If with_mutual_shadowing = true and BVH is not yet built for either shape, it is built automatically. To avoid this, pre-build with build_bvh! or pass with_bvh=true when loading shapes.

BinaryAsteroidThermoPhysicalProblem(shape, thermo_params; kwargs...) -> problem

Convenience constructor: build both single-body problems from raw shapes and parameters.

Arguments

• shape : Tuple (shape1, shape2) of shape models
• thermo_params : Tuple (thermo_params1, thermo_params2) of thermophysical parameters

Keyword Arguments

Single-body kwargs (applied identically to both bodies):

• with_self_shadowing = true
• with_self_heating = true
• upper_boundary_condition = RadiationBoundaryCondition()
• lower_boundary_condition = InsulationBoundaryCondition()

Binary-system kwargs:

• with_mutual_shadowing = true
• with_mutual_heating = true

Notes

For per-body control of single-body kwargs, use SingleAsteroidThermoPhysicalProblem separately and pass the results to BinaryAsteroidThermoPhysicalProblem(primary, secondary; ...).

Example

problem = BinaryAsteroidThermoPhysicalProblem(
    (shape1, shape2),
    (thermo_params1, thermo_params2);
    with_mutual_shadowing = true,
    with_mutual_heating   = true,
)

struct BinaryAsteroidThermoPhysicalSolution

Solution data for a binary asteroid thermophysical simulation.

Fields

• primary : Solution for the primary body
• secondary : Solution for the secondary body

BinaryAsteroidThermoPhysicalSolution(state, ephem, output)

Allocate a binary solution from the given state, ephemerides, and output specification.

struct BinaryAsteroidThermoPhysicalState <: AbstractAsteroidThermoPhysicalState

Internal simulation state for a binary-asteroid thermophysical model.

Fields

• problem : Binary problem definition (mutual shadowing/heating flags)
• primary : Simulation state for the primary body
• secondary : Simulation state for the secondary body

Invariant

The inner constructor enforces:

state.primary.problem   === state.problem.primary
state.secondary.problem === state.problem.secondary

Use _build_binary_state rather than constructing directly to guarantee consistency.

CrankNicolson()

Crank-Nicolson method for heat conduction:

• Second order in time
• Unconditionally stable

References

• https://en.wikipedia.org/wiki/Crank–Nicolson_method

Internal cache for the Crank-Nicolson method. Holds pre-allocated vectors for the tridiagonal matrix algorithm.

References

• https://en.wikipedia.org/wiki/Crank–Nicolson_method

ExplicitEuler()

Explicit (forward) Euler method for heat conduction:

• First order in time
• Conditionally stable (requires Fourier number λ < 0.5)

Internal cache for the explicit (forward) Euler method. Holds a pre-allocated vector for the temperature at the next time step.

Abstract type for internal pre-allocated caches used in heat conduction computations. These types are not part of the public API.

ImplicitEuler()

Implicit (backward) Euler method for heat conduction:

• First order in time
• Unconditionally stable

Internal cache for the implicit (backward) Euler method. Holds pre-allocated vectors for the tridiagonal matrix algorithm.

Singleton type of insulation boundary condition

Type of isothermal boundary condition

Singleton type of radiation boundary condition

struct SingleAsteroidEphemerides{R} <: AbstractSingleAsteroidEphemerides

Ephemerides for a single-asteroid thermophysical simulation.

The type parameter R controls whether force and torque are computed:

• R = Nothing : temperature only; R_body_to_inertial = nothing
• R = Vector{SMatrix{3,3,Float64,9}} : force and torque are computed in the inertial frame

Fields

• times : Simulation timesteps [s]
• r_sun : Sun position vector in the body-fixed frame at each timestep [m]
• R_body_to_inertial : Passive rotation matrices from body-fixed to inertial frame, or nothing when force/torque are not needed.

Constructors

SingleAsteroidEphemerides(times, r_sun)
    -> SingleAsteroidEphemerides{Nothing}
SingleAsteroidEphemerides(times, r_sun, nothing)
    -> SingleAsteroidEphemerides{Nothing}
SingleAsteroidEphemerides(times, r_sun, R_body_to_inertial)
    -> SingleAsteroidEphemerides{Vector{SMatrix{3,3,Float64,9}}}

An AbstractRange (e.g. range(et_begin, et_end; length=n)) may be passed as times and is automatically collected to Vector{Float64}. Plain AbstractVector / AbstractMatrix elements in r_sun and R_body_to_inertial are automatically converted to the corresponding SVector{3,Float64} / SMatrix{3,3,Float64,9} types, so importing StaticArrays in user code is not required.

struct SingleAsteroidOutputSpec

Output specification for a single-asteroid thermophysical simulation.

Encapsulates which timesteps, face indices, and physical quantities to record.

Fields

• output_times : Timesteps at which data are saved [s]; must be a subset of ephem.times
• subsurface_face_ids : Face indices for which to save subsurface temperature profiles
• save_surface_temperature : Save surface temperature at output_times (default: true)
• save_subsurface_temperature : Save subsurface temperature profiles at output_times (default: true)
• save_face_forces : Save per-face thermal forces at output_times (default: false)
• save_forces : Save net thermal force at output_times (default: false)
• save_torques : Save net thermal torque at output_times (default: false)

Notes

• save_subsurface_temperature = true requires a non-empty subsurface_face_ids.
• save_forces and save_torques require ephemerides with R_body_to_inertial (i.e., SingleAsteroidEphemerides{<:AbstractVector}); using them with rotation-free ephemerides raises an ArgumentError at solve time.

Example

output = SingleAsteroidOutputSpec(output_times, subsurface_face_ids;
    save_surface_temperature    = true,
    save_subsurface_temperature = true,
    save_face_forces            = false,
    save_forces                 = true,
    save_torques                = true,
)

struct SingleAsteroidThermoPhysicalProblem

Defines the thermophysical problem for a single asteroid. Encapsulates all information needed to describe the physical problem, separate from the numerical method used to solve it.

Fields

• shape : Shape model of the asteroid (ShapeModel or HierarchicalShapeModel)
• thermo_params : Thermophysical parameters
• with_self_shadowing : Whether to include self-shadowing
• with_self_heating : Whether to include self-heating (re-absorption of thermal emission from other faces)
• upper_boundary_condition : Boundary condition at the surface (upper boundary)
• lower_boundary_condition : Boundary condition at depth (lower boundary)

Usage

problem = SingleAsteroidThermoPhysicalProblem(shape, thermo_params;
    with_self_shadowing = true,
    with_self_heating   = true,
    upper_boundary_condition = RadiationBoundaryCondition(),
    lower_boundary_condition = InsulationBoundaryCondition(),
)
solution = solve(problem, CrankNicolson(); ephem = ephem, T₀ = 200.0)

SingleAsteroidThermoPhysicalProblem(shape, thermo_params; kwargs...) -> problem

Construct a thermophysical problem for a single asteroid.

Arguments

• shape : Shape model of the asteroid
• thermo_params : Thermophysical parameters

Keyword Arguments

• with_self_shadowing = true : Whether to include self-shadowing
• with_self_heating = true : Whether to include self-heating
• upper_boundary_condition = RadiationBoundaryCondition() : Boundary condition at the surface
• lower_boundary_condition = InsulationBoundaryCondition() : Boundary condition at depth

Notes

• If with_self_shadowing = true and face_visibility_graph is not yet built, it is built automatically. To avoid this, pre-build with build_face_visibility_graph! or pass with_face_visibility=true when loading the shape.
• Thermophysical parameters are broadcast to all faces via broadcast_thermo_params!.

struct SingleAsteroidThermoPhysicalSolution

Solution data for a single asteroid thermophysical simulation.

Fields

Saved at all timesteps

• times : All simulation timesteps [s]
• absorbed_power : Total absorbed power on the whole surface [W]
• emitted_power : Total emitted thermal radiation power from the whole surface [W]

Metadata

• output : Output specification (controls which data are saved and when)
• depth_nodes : Depth of each subsurface calculation node [m], size (n_depth,)

Saved only at output.output_times (nothing when the corresponding flag is false)

• surface_temperature : Surface temperature [K], size (n_face, n_save), or nothing
• subsurface_temperature : Subsurface temperature [K] by face ID, each entry (n_depth, n_save), or nothing
• face_forces : Per-face thermal force in the body-fixed frame [N], size (n_face, n_save), or nothing
• forces : Net thermal force in the inertial frame [N], size (n_save,), or nothing
• torques : Net thermal torque in the inertial frame [N⋅m], size (n_save,), or nothing

Notes

• forces and torques are non-nothing only when the ephemerides include R_body_to_inertial (i.e., SingleAsteroidEphemerides{<:AbstractVector}) and the corresponding flag in output is true.

SingleAsteroidThermoPhysicalSolution(state, ephem, output)

Allocate a solution from the given state, ephemerides, and output specification.

struct SingleAsteroidThermoPhysicalState <: AbstractAsteroidThermoPhysicalState

Internal simulation state for a single-asteroid thermophysical model. Holds the mutable arrays that evolve during a solve call. The problem definition (shape, parameters, flags, boundary conditions) is accessed via the problem field to avoid duplication.

Fields

• problem : Problem definition (shape, thermo_params, flags, BCs)
• solver_cache : Pre-allocated cache for the heat-conduction solver
• illuminated_faces : Illumination flag for each face
• flux_sun : Direct solar flux on each face [W/m²]
• flux_scat : Scattered-light flux on each face [W/m²]
• flux_rad : Thermal-emission flux from surrounding faces [W/m²]
• temperature : Temperature matrix (n_depth, n_face) [K]
• face_forces : Thermal recoil force on each face [N]
• force : Net thermal recoil force in body-fixed frame [N]
• torque : Net thermal recoil torque in body-fixed frame [N⋅m]

struct ThermoParams

Fields

• thermal_conductivity : Vector of thermal conductivity for each facet [W/m/K]

• density : Vector of density for each facet [kg/m³]

• heat_capacity : Vector of heat capacity for each facet [J/kg/K]

• reflectance_vis : Vector of reflectance in visible light for each facet [-]

• reflectance_ir : Vector of reflectance in thermal infrared for each facet [-]

• emissivity : Vector of emissivity for each facet [-]

• z_max : Depth of the lower boundary of a heat conduction equation [m]

• Δz : Depth step width [m]

• n_depth : Number of depth steps

ThermoParams(
    thermal_conductivity ::Float64,
    density              ::Float64,
    heat_capacity        ::Float64,
    reflectance_vis      ::Float64,
    reflectance_ir       ::Float64,
    emissivity           ::Float64,
    z_max                ::Float64,
    Δz                   ::Float64,
    n_depth              ::Int
)

Outer constructor for ThermoParams. You can give the same parameters to all facets by Float64.

Arguments

• thermal_conductivity : Thermal conductivity [W/m/K]

• density : Density [kg/m³]

• heat_capacity : Heat capacity [J/kg/K]

• reflectance_vis : Reflectance in visible light [-]

• reflectance_ir : Reflectance in thermal infrared [-]

• emissivity : Emissivity [-]

• z_max : Depth of the lower boundary of a heat conduction equation [m]

• Δz : Depth step width [m]

• n_depth : Number of depth steps

absorbed_energy_flux(R_vis, R_ir, F_sun, F_scat, F_rad) -> F_abs

Calculate the total energy flux absorbed by a surface element, accounting for wavelength-dependent reflectance properties.

Arguments

• R_vis::Real : Reflectance for visible light [-], valid between 0 and 1.
• R_ir::Real : Reflectance for thermal infrared [-], valid between 0 and 1.
• F_sun::Real : Direct solar radiation flux [W/m²]
• F_scat::Real : Scattered sunlight flux from other surfaces [W/m²]
• F_rad::Real : Thermal radiation flux from surrounding surfaces [W/m²]

Returns

• F_abs::Real : Total absorbed energy flux [W/m²]

Mathematical Formula

F_abs = (1 - R_vis) × F_sun + (1 - R_vis) × F_scat + (1 - R_ir) × F_rad

Physical Interpretation

The function accounts for different reflectance properties at different wavelengths:

• Solar radiation (Fsun) and scattered light (Fscat) are in the visible spectrum
• Thermal radiation (F_rad) is in the infrared spectrum
• The absorbed fraction is (1 - reflectance) for each component

Example

R_vis = 0.1   # 10% reflectance in visible
R_ir = 0.05   # 5% reflectance in IR
F_sun = 1000.0   # Direct solar flux
F_scat = 50.0    # Scattered light
F_rad = 100.0    # Thermal radiation
F_abs = absorbed_energy_flux(R_vis, R_ir, F_sun, F_scat, F_rad)
# Returns: 0.9 × 1000 + 0.9 × 50 + 0.95 × 100 = 1040.0 W/m²

analytical_solution_isothermal(x, t, L, α; n_max=100) -> T

Calculate the analytical solution of the 1D heat equation with isothermal boundary conditions.

• Equation: ∂T/∂t = α ∂²T/∂x²
• Domain: 0 ≤ x ≤ L
• Boundary conditions: T(0,t) = T(L,t) = 0
• Initial condition: T₀(x) = x < 0.5L ? 2x/L : 2(1 - x/L) # Triangular profile as follows: T₀ ^ 1 | ・ | ・ ・ | ・ ・ |・ ・ 0 +–-+–-+–> x 0 L/2 L

The solution is given by the Fourier series: T(x, t) = Σ Bₙ * sin(nπx/L) * exp(-αn²π²t/L²) where Bₙ = (2/L) * ∫₀^L T₀(ξ) * sin(nπξ/L) dξ For the triangular initial condition, the coefficients can be calculated analytically: Bₙ = (8/n²π²) * sin(nπ/2) only for odd n. The sum of even-n terms is zero due to symmetry.

Arguments

• x : Position [m]
• t : Time [s]
• L : Length of the domain [m]
• α : Thermal diffusivity [m²/s]
• n_max : Number of terms in the Fourier series

Returns

• T : Temperature [K]

blackbody_radiance(λ, T) -> L_λ

Accroding to Planck's law, calculate the spectral intensity of blackbody radiation at wavelength λ and temperature T.

Arguments

• λ : Wavelength [m]
• T : Temperature [K]

Return

• L_λ : Spectral radiance [W/m²/m/steradian]

cf. https://github.com/JuliaAstro/Planck.jl/blob/main/src/Planck.jl

blackbody_radiance(T) -> L

According to Stefan-Boltzmann law, calculate the total radiance of blackbody radiation at temperature T, integrated over all wavelength.

Arguments

• T : Temperature [K]

Return

• L : Radiance [W/m²/steradian]

broadcast_thermo_params!(thermo_params::ThermoParams, shape::ShapeModel)

Broadcast the thermophysical parameters to all faces if the values are uniform globally.

Arguments

• thermo_params : Thermophysical parameters
• shape : Shape model

broadcast_thermo_params!(thermo_params::ThermoParams, n_face::Int)

Broadcast the thermophysical parameters to all faces if the values are uniform globally.

Arguments

• thermo_params : Thermophysical parameters
• n_face : Number of faces on the shape model

crank_nicolson!(state::SingleAsteroidThermoPhysicalState, Δt)

Solve the 1D heat conduction equation using the Crank-Nicolson method. This method combines the explicit and implicit Euler methods for improved accuracy.

Arguments

• state::SingleAsteroidThermoPhysicalState : Thermophysical simulation state for a single asteroid
• Δt::Real : Time step [s]

Method Properties

• Time discretization: Semi-implicit (average of forward and backward differences)
• Accuracy: Second-order in both time and space
• Stability: Unconditionally stable for any time step size

Discretization

The heat conduction equation ∂T/∂t = α∂²T/∂z² is discretized using the average of explicit and implicit schemes:

T[i,n+1] - T[i,n] = (α∆t)/(2∆z²) × 
    [(T[i+1,n+1] - 2T[i,n+1] + T[i-1,n+1]) + (T[i+1,n] - 2T[i,n] + T[i-1,n])]

This leads to a tridiagonal system:

-rT[i-1,n+1] + (1+2r)T[i,n+1] - rT[i+1,n+1] = 
    rT[i-1,n] + (1-2r)T[i,n] + rT[i+1,n]

where r = α∆t/(2∆z²)

Advantages

• Higher accuracy than both explicit and implicit Euler methods
• Unconditionally stable
• Optimal balance between accuracy and computational cost
• Second-order accuracy in both time and space

Implementation Details

The method requires solving a tridiagonal system at each time step, similar to the implicit Euler method but with a modified right-hand side that includes information from the current time step.

See Also

• tridiagonal_matrix_algorithm! for the solution algorithm
• implicit_euler!, explicit_euler! for comparison with other methods

explicit_euler!(state::SingleAsteroidThermoPhysicalState, Δt)

Solve the 1D heat conduction equation using the explicit (forward) Euler method. This method is conditionally stable and requires careful time step selection.

Arguments

• state::SingleAsteroidThermoPhysicalState : Thermophysical simulation state for a single asteroid
• Δt::Real : Time step [s]

Method Properties

• Time discretization: Explicit (forward difference)
• Accuracy: First-order in time, second-order in space
• Stability: Conditionally stable, requires λ = αΔt/Δz² < 0.5

Discretization

The heat conduction equation ∂T/∂t = α∂²T/∂z² is discretized as:

T[i,n+1] = T[i,n] + λ(T[i+1,n] - 2T[i,n] + T[i-1,n])

where:

• λ = αΔt/Δz² is the dimensionless time step
• α = k/(ρCₚ) is the thermal diffusivity
• n is the time index, i is the depth index

Stability Criterion

The method is stable only when λ < 0.5. If this condition is violated, an error is thrown.

Boundary Conditions

• Upper boundary: Determined by update_upper_temperature!
• Lower boundary: Determined by update_lower_temperature!

Performance Notes

• This method is simple and fast but requires small time steps for stability
• Consider using implicit methods for larger time steps

Errors

• Throws an error if λ ≥ 0.5 (stability violation)

export_solution(dirpath, solution::BinaryAsteroidThermoPhysicalSolution)

Export results for both bodies to dirpath/primary/ and dirpath/secondary/.

export_solution(dirpath, solution::SingleAsteroidThermoPhysicalSolution)

Export simulation results to CSV files in dirpath.

Files written depend on the output specification:

• diagnostics.csv : always (absorbed_power, emitted_power at all timesteps)
• surface_temperature.csv : when output.save_surface_temperature = true
• subsurface_temperature.csv : when output.save_subsurface_temperature = true
• thermal_face_forces.csv : when output.save_face_forces = true
• thermal_net_forces.csv : when output.save_forces = true or output.save_torques = true

implicit_euler!(state::SingleAsteroidThermoPhysicalState, Δt)

Solve the 1D heat conduction equation using the implicit (backward) Euler method. This method is unconditionally stable, allowing for larger time steps than explicit methods.

Arguments

• state::SingleAsteroidThermoPhysicalState : Thermophysical simulation state for a single asteroid
• Δt::Real : Time step [s]

Method Properties

• Time discretization: Implicit (backward difference)
• Accuracy: First-order in time, second-order in space
• Stability: Unconditionally stable for any time step size

Discretization

The heat conduction equation ∂T/∂t = α∂²T/∂z² is discretized as:

T[i,n+1] - T[i,n] = λ(T[i+1,n+1] - 2T[i,n+1] + T[i-1,n+1])

This leads to a tridiagonal system:

-λT[i-1,n+1] + (1+2λ)T[i,n+1] - λT[i+1,n+1] = T[i,n]

where λ = αΔt/Δz²

Solution Method

The resulting tridiagonal system is solved using the Thomas algorithm (tridiagonal matrix algorithm) for each face.

Boundary Conditions

Different boundary conditions modify the tridiagonal matrix:

• Radiation BC: Special treatment after solving the system
• Insulation BC: Modified coefficients at boundaries
• Isothermal BC: Direct temperature assignment

Advantages

• Unconditionally stable - no restriction on time step size
• Allows for larger time steps compared to explicit methods
• More computationally intensive per step but often faster overall

See Also

• tridiagonal_matrix_algorithm! for the solution algorithm
• update_upper_temperature!, update_lower_temperature! for boundary conditions

init_temperature!(state::BinaryAsteroidThermoPhysicalState, T₀_primary, T₀_secondary)

Initialize temperatures with separate values for the primary and secondary bodies.

Each argument can be a Real (uniform) or an AbstractMatrix of size (n_depth, n_face).

Arguments

• state : Thermophysical simulation state for a binary asteroid
• T₀_primary : Initial temperature for the primary body [K]
• T₀_secondary : Initial temperature for the secondary body [K]

init_temperature!(state::BinaryAsteroidThermoPhysicalState, T₀::Real)

Initialize all temperature cells in both bodies at the uniform temperature T₀.

Arguments

• state : Thermophysical simulation state for a binary asteroid
• T₀ : Initial temperature of all cells [K]

init_temperature!(state::SingleAsteroidThermoPhysicalState, T₀::AbstractMatrix)

Initialize temperatures from a full depth–face temperature matrix. The matrix must have size (n_depth, n_face), matching state.temperature.

Arguments

• state : Thermophysical simulation state for a single asteroid
• T₀ : Temperature matrix of size (n_depth, n_face) [K]

init_temperature!(state::SingleAsteroidThermoPhysicalState, T₀::Real)

Initialize all temperature cells at the uniform temperature T₀.

Arguments

• state : Thermophysical simulation state for a single asteroid
• T₀ : Initial temperature [K]

integrate_absorbed_power(state::SingleAsteroidThermoPhysicalState) -> Float64

Integrate the absorbed energy flux over all surface facets to obtain total absorbed power [W]:

P_abs = Σᵢ F_abs,ᵢ × Aᵢ

See Also

• integrate_emitted_power for the total emitted power
• absorbed_energy_flux for the per-facet flux calculation

integrate_emitted_power(state::SingleAsteroidThermoPhysicalState) -> Float64

Integrate the thermal emission over all surface facets to obtain total emitted power [W]:

P_emit = Σᵢ εᵢ × σ × Tᵢ⁴ × Aᵢ

In thermal equilibrium, integrate_emitted_power ≈ integrate_absorbed_power.

See Also

• integrate_absorbed_power for the total absorbed power

mutual_heating!(state::BinaryAsteroidThermoPhysicalState, r₁₂, R₂₁)

Calculate the mutual heating between the primary and secondary asteroids.

Arguments

• state::BinaryAsteroidThermoPhysicalState : Thermophysical simulation state for a binary asteroid
• r₁₂::StaticVector{3} : Position vector of secondary's center in primary's frame [m]
• R₂₁::StaticMatrix{3,3} : Rotation matrix from secondary to primary frame

TODO

• Need to consider local horizon?

record_timestep!(solution, state, i_time, R₁ᵢ, R₂ᵢ)

Record simulation data for both bodies of a binary system at timestep i_time, including force and torque rotated to the inertial frame for each body.

Arguments

• solution : Solution container (BinaryAsteroidThermoPhysicalSolution)
• state : Current simulation state (BinaryAsteroidThermoPhysicalState)
• i_time : Index into solution.primary.times for the current timestep
• R₁ᵢ : Rotation matrix from the primary body-fixed frame to the inertial frame
• R₂ᵢ : Rotation matrix from the secondary body-fixed frame to the inertial frame

record_timestep!(solution, state, i_time)

Record simulation data for both bodies of a binary system at timestep i_time. Delegates to the single-body form for each body independently.

Arguments

• solution : Solution container (BinaryAsteroidThermoPhysicalSolution)
• state : Current simulation state (BinaryAsteroidThermoPhysicalState)
• i_time : Index into solution.primary.times for the current timestep

record_timestep!(solution, state, i_time, R)

Record simulation data for a single asteroid at timestep i_time, including force and torque rotated to the inertial frame.

Extends the 3-argument form by additionally recording net thermal force and torque (when save_forces/save_torques are true) after rotating from the body-fixed frame to the inertial frame via R.

Arguments

• solution : Solution container (SingleAsteroidThermoPhysicalSolution)
• state : Current simulation state (SingleAsteroidThermoPhysicalState)
• i_time : Index into solution.times for the current timestep
• R : Rotation matrix from the body-fixed frame to the inertial frame

record_timestep!(solution, state, i_time)

Record simulation data for a single asteroid at timestep i_time.

Saves absorbed_power and emitted_power at every timestep. At timesteps that coincide with solution.output.output_times, also records snapshot data (surface temperature, subsurface temperature, face forces) according to the flags in solution.output.

Arguments

• solution : Solution container (SingleAsteroidThermoPhysicalSolution)
• state : Current simulation state (SingleAsteroidThermoPhysicalState)
• i_time : Index into solution.times for the current timestep

subsolar_temperature(r☉, R_vis, ε) -> Tₛₛ

Calculate the subsolar equilibrium temperature at a given heliocentric distance.

Arguments

• r☉ : Sun's position vector in the asteroid-fixed frame [m]
• R_vis : Visible-light reflectance (Bond albedo) [-]
• ε : Emissivity [-]

Returns

• Tₛₛ::Float64 : Subsolar point temperature [K]

Notes

Assumes instantaneous radiative equilibrium (zero thermal inertia). Useful as an upper bound for surface temperatures and as an initial guess for T₀.

Mathematical Formula

\[T_{ss} = \left[\frac{(1 - A) \Phi_\odot}{\varepsilon \sigma}\right]^{1/4}\]

where $\Phi_\odot = \Phi_0 / r^2$ is the solar flux at heliocentric distance $r$.

surface_temperature(state::SingleAsteroidThermoPhysicalState) -> T_surface

Extract the surface temperature (uppermost layer) for all faces.

Returns

• T_surface::Vector{Float64} : Surface temperature for each face [K]

thermal_diffusivity(k, ρ, Cp) -> α

Calculate the thermal diffusivity of a material.

Arguments

• k::Real : Thermal conductivity [W/m/K]
• ρ::Real : Material density [kg/m³]
• Cₚ::Real : Heat capacity [J/kg/K]

Returns

• α::Real : Thermal diffusivity [m²/s]

Mathematical Formula

\[\alpha = \frac{k}{\rho C_p}\]

Physical Meaning

• Measures how quickly temperature propagates through material
• Appears in the heat diffusion equation: ∂T/∂t = α∇²T
• High α: rapid heat diffusion
• Low α: slow heat diffusion

thermal_inertia(k, ρ, Cp) -> Γ

Calculate the thermal inertia of a material.

Arguments

• k::Real : Thermal conductivity [W/m/K]
• ρ::Real : Material density [kg/m³]
• Cₚ::Real : Heat capacity [J/kg/K]

Returns

• Γ::Real : Thermal inertia [J m⁻² K⁻¹ s⁻¹/²]

Mathematical Formula

\[\Gamma = \sqrt{k \rho C_p}\]

Physical Meaning

• Measures resistance to temperature change
• High Γ: slow temperature response (rock-like)
• Low Γ: rapid temperature response (dust-like)
• Typical values: 50-2500 J m⁻² K⁻¹ s⁻¹/² for a planetary surface

Note

The unit is sometimes called "tiu" (thermal inertia unit).

thermal_radiance(shape, emissivities, temperatures, obs) -> L

Calculate the radiance from the temperature distribution based on a shape model.

Arguments

• shape : Shape model of an asteroid
• emissivities : Emissivity of each facet of the shape model [-]
• temperatures : Temperature of each facet of the shape model [K]
• obs : Position vector of the observer in the same coordinate system as shape [m]

Return

• L : Radiance [W/m²]

thermal_skin_depth(P, k, ρ, Cp) -> l_2π

Calculate the thermal skin depth for a periodic temperature variation.

Arguments

• P::Real : Period of thermal cycle [s]
• k::Real : Thermal conductivity [W/m/K]
• ρ::Real : Material density [kg/m³]
• Cₚ::Real : Heat capacity [J/kg/K]

Returns

• l_2π::Real : Thermal skin depth [m]

Mathematical Formula

The thermal skin depth is defined as:

\[l_{2\pi} = \sqrt{\frac{4\pi P k}{\rho C_p}}\]

Physical Meaning

• Represents the e-folding depth of temperature variations
• Temperature amplitude decreases by factor e^(-2π) ≈ 0.0019 at this depth
• Useful for determining computational domain depth

Reference

• Rozitis & Green (2011), MNRAS 415, 2042-2062

tridiagonal_matrix_algorithm!(a, b, c, d, x)
tridiagonal_matrix_algorithm!(state::SingleAsteroidThermoPhysicalModel)

Tridiagonal matrix algorithm to solve the heat conduction equation by the implicit (backward) Euler and Crank-Nicolson methods.

| b₁ c₁ 0  ⋯  0   | | x₁ |   | d₁ |
| a₂ b₂ c₂ ⋯  0   | | x₂ |   | d₂ |
| 0  a₃ b₃ ⋯  0   | | x₃ | = | d₃ |
| ⋮  ⋮  ⋮  ⋱  cₙ₋₁| | ⋮  |   | ⋮  |
| 0  0  0  aₙ bₙ  | | xₙ |   | dₙ |

References

• https://en.wikipedia.org/wiki/Tridiagonalmatrixalgorithm

update_flux_all!(state::BinaryAsteroidThermoPhysicalState, r☉₁::StaticVector{3}, r₁₂::StaticVector{3}, R₁₂::StaticMatrix{3,3})

Update all energy fluxes (solar, scattered, thermal radiation) to the surface for a binary asteroid. This is a convenience function that computes necessary coordinate transformations and calls individual flux update functions.

Arguments

• state::BinaryAsteroidThermoPhysicalState : Thermophysical simulation state for a binary asteroid
• r☉₁::StaticVector{3} : Sun's position in the primary's body-fixed frame (NOT normalized) [m]
• r₁₂::StaticVector{3} : Position vector of secondary's center in primary's frame [m]
• R₁₂::StaticMatrix{3,3} : Rotation matrix from primary to secondary frame

Algorithm

1. Computes all necessary coordinate transformations
2. Updates solar flux considering eclipse (mutual shadowing)
3. Updates scattered light flux (self-heating)
4. Updates thermal radiation flux (self-heating)
5. Applies mutual heating between components

Notes

• This function internally handles all coordinate transformations
• Automatically respects SELFSHADOWING, SELFHEATING, MUTUALSHADOWING, and MUTUALHEATING flags

update_flux_all!(state::SingleAsteroidThermoPhysicalState, r☉::StaticVector{3})

Update all energy fluxes (solar, scattered, thermal radiation) to the surface for a single asteroid.

Arguments

• state::SingleAsteroidThermoPhysicalState : Thermophysical simulation state for a single asteroid
• r☉::StaticVector{3} : Sun's position in the asteroid-fixed frame (NOT normalized) [m]

Algorithm

1. Updates direct solar flux on all faces considering self-shadowing
2. Updates scattered sunlight flux from other faces (self-heating)
3. Updates thermal radiation flux from other faces (self-heating)

Notes

• This is a convenience function that calls all individual flux update functions
• Automatically respects SELFSHADOWING and SELFHEATING flags

update_flux_rad_single!(state::BinaryAsteroidThermoPhysicalState)

Update flux of absorption of thermal radiation from surrounding surface. Single radiation-absorption is only considered, assuming albedo is close to zero at thermal infrared wavelength.

Arguments

• state : Thermophysical simulation state for a binary asteroid

update_flux_rad_single!(state::SingleAsteroidThermoPhysicalState)

Update flux of absorption of thermal radiation from surrounding surface. Single radiation-absorption is only considered, assuming albedo is close to zero at thermal infrared wavelength.

Arguments

• state : Thermophysical simulation state for a single asteroid

update_flux_scat_single!(state::BinaryAsteroidThermoPhysicalState)

Update flux of scattered sunlight, only considering single scattering.

Arguments

• state : Thermophysical simulation state for a binary asteroid

update_flux_scat_single!(state::SingleAsteroidThermoPhysicalState)

Update flux of scattered sunlight, only considering single scattering.

Arguments

• state : Thermophysical simulation state for a single asteroid

update_flux_sun!(
    state::BinaryAsteroidThermoPhysicalState,
    r☉₁::StaticVector{3},   r☉₂::StaticVector{3}, 
    r₁₂::StaticVector{3},   r₂₁::StaticVector{3},
    R₁₂::StaticMatrix{3,3}, R₂₁::StaticMatrix{3,3},
)

Update solar irradiation flux on both components of a binary asteroid system with mutual shadowing.

Arguments

• state::BinaryAsteroidThermoPhysicalState : Thermophysical simulation state for a binary asteroid
• r☉₁::StaticVector{3} : Sun's position vector in the primary's body-fixed frame (NOT normalized) [m]
• r☉₂::StaticVector{3} : Sun's position vector in the secondary's body-fixed frame (NOT normalized) [m]
• r₁₂::StaticVector{3} : Position vector of secondary's center in primary's frame [m]
• r₂₁::StaticVector{3} : Position vector of primary's center in secondary's frame [m]
• R₁₂::StaticMatrix{3,3} : Rotation matrix from primary to secondary frame
• R₂₁::StaticMatrix{3,3} : Rotation matrix from secondary to primary frame

Notes

• All coordinate transformations should be pre-computed by the caller
• Uses the new apply_eclipse_shadowing! API from AsteroidShapeModels.jl v0.4.1
• Requires BVH to be built for both shapes (should be done when loading with with_bvh=true)
• Combines self-shadowing and mutual shadowing in a single call

update_flux_sun!(state::SingleAsteroidThermoPhysicalState, r☉::StaticVector{3})

Update the direct solar irradiation flux on every face of the asteroid.

Arguments

• state::SingleAsteroidThermoPhysicalState : Thermophysical simulation state for a single asteroid
• r☉::StaticVector{3} : Position vector from asteroid to Sun in body-fixed frame (NOT normalized) [m]

Algorithm

For each face, the solar flux is calculated as:

1. Solar flux at asteroid's location: F☉ = SOLAR_CONST / distance²
2. Normalize sun direction: r̂☉ = r☉ / |r☉|
3. Face flux: F_sun = F☉ × max(0, n̂ · r̂☉)

where n̂ is the face normal. If SELF_SHADOWING is enabled, the function also checks whether each face is shadowed by other parts of the asteroid.

Notes

• The input vector r☉ must not be normalized (used for distance calculation)
• Faces with negative dot product (facing away from Sun) receive zero flux
• Shadowed faces (when SELF_SHADOWING = true) also receive zero flux

update_lower_temperature!(state::SingleAsteroidThermoPhysicalState)

Update the temperature at the lower boundary (deepest layer) based on the boundary condition.

Arguments

• state::SingleAsteroidThermoPhysicalState : Thermophysical simulation state for a single asteroid

Boundary Conditions

The function applies one of the following boundary conditions at the bottom of the computational domain:

1. Insulation (Neumann): ∂T/∂z = 0

   • No heat flux through the lower boundary
   • Temperature gradient is zero: T[end] = T[end-1]
   • Most commonly used for asteroid modeling

2. Isothermal (Dirichlet): T = T_iso

   • Fixed temperature at the lower boundary
   • Used when deep interior temperature is known
   • T[end] = state.problem.lower_boundary_condition.T_iso

Notes

• This function is called after solving the heat conduction equation
• For explicit Euler method, it directly updates the temperature vector
• The lower boundary should be deep enough that the chosen condition doesn't affect surface temperatures

update_surface_temperature!(T::AbstractVector, F_abs::Real, k::Real, ρ::Real, Cₚ::Real, ε::Real, Δz::Real)

Newton's method to update the surface temperature under radiation boundary condition.

Arguments

• T : 1-D array of temperatures
• F_abs : Total energy flux absorbed by the facet
• k : Thermal conductivity [W/m/K]
• ρ : Density [kg/m³]
• Cₚ : Heat capacity [J/kg/K]
• ε : Emissivity [-]
• Δz : Depth step width [m]

update_temperature!(state::BinaryAsteroidThermoPhysicalState, Δt)

Calculate the temperature for the next time step based on 1D heat conductivity equation.

Arguments

• state : Thermophysical simulation state for a binary asteroid
• Δt : Time step [s]

update_temperature!(state::SingleAsteroidThermoPhysicalState, Δt)

Update the temperature distribution for the next time step by solving the 1D heat conduction equation. The solver method is determined by state.solver_cache, and special handling is applied for zero conductivity.

Arguments

• state::SingleAsteroidThermoPhysicalState : Thermophysical simulation state for a single asteroid
• Δt::Real : Time step [s]

Solver Selection

The function automatically selects the appropriate solver based on state.solver_cache:

• ExplicitEulerCache: Forward Euler method (conditionally stable, requires λ < 0.5)
• ImplicitEulerCache: Backward Euler method (unconditionally stable)
• CrankNicolsonCache: Crank-Nicolson method (unconditionally stable, second-order accurate)

Special Cases

• If thermal conductivity is zero, calls update_temperature_zero_conductivity! instead
• The zero-conductivity case uses instantaneous radiative equilibrium

Mathematical Background

Solves the 1D heat conduction equation:

∂T/∂t = α ∂²T/∂z²

where α = k/(ρCₚ) is the thermal diffusivity.

See Also

• explicit_euler!, implicit_euler!, crank_nicolson! for specific solver implementations
• update_temperature_zero_conductivity! for the zero-conductivity case

update_temperature_zero_conductivity!(state::SingleAsteroidThermoPhysicalState)

Update surface temperature for the zero thermal conductivity case. When thermal conductivity is zero, there is no heat conduction into the subsurface, and the surface temperature is determined solely by instantaneous radiative equilibrium.

Arguments

• state::SingleAsteroidThermoPhysicalState : Thermophysical simulation state for a single asteroid

Mathematical Formula

For each face, the surface temperature T is calculated from:

εσT⁴ = (1-Rᵥᵢₛ)F_sun + (1-Rᵥᵢₛ)F_scat + (1-Rᵢᵣ)F_rad

where:

• ε : Emissivity
• σ : Stefan-Boltzmann constant
• Rᵥᵢₛ : Reflectance in visible light
• Rᵢᵣ : Reflectance in thermal infrared
• F_sun : Direct solar flux
• F_scat : Scattered light flux
• F_rad : Thermal radiation flux from surrounding surfaces

Notes

• This function is called when state.problem.thermo_params.thermal_conductivity is zero
• The temperature instantly adjusts to balance incoming and outgoing radiation
• No subsurface temperatures are updated (only surface layer)

update_thermal_force!(state::BinaryAsteroidThermoPhysicalState)

Calculate the thermal force and torque on every face and integrate them over all faces.

Arguments

• state : Thermophysical simulation state for a binary asteroid

update_thermal_force!(state::SingleAsteroidThermoPhysicalState)

Calculate the thermal recoil force (Yarkovsky effect) and torque (YORP effect) on the asteroid by integrating photon momentum from thermal emission and reflection over all surface facets.

Arguments

• state::SingleAsteroidThermoPhysicalState : Thermophysical simulation state for a single asteroid

Physics

The function calculates non-gravitational effects caused by anisotropic photon emission:

• Yarkovsky effect: Net force due to thermal lag causing asymmetric emission
• YORP effect: Net torque changing the asteroid's rotation state

Algorithm

For each facet i, the thermal force is computed as:

F_i = -(2/3) × (E_i × A_i)/c × n̂_i + Σⱼ (E_i × A_i)/c × f_ij × d̂_ij

where:

• E_i = total emittance from facet i (reflection + thermal emission) [W/m²]
• A_i = area of facet i [m²]
• c = speed of light [m/s]
• n̂_i = outward normal vector of facet i
• f_ij = view factor from facet i to j
• d̂_ij = unit vector from facet i to j

The first term represents direct photon recoil normal to the surface. The second term accounts for photons intercepted by other facets (self-heating contribution).

Outputs (stored in state)

• state.face_forces : Thermal force vector on each facet [N]
• state.force : Total thermal force in body-fixed frame [N]
• state.torque : Total thermal torque in body-fixed frame [N⋅m]

Physical Significance

• The force causes orbital drift (Yarkovsky effect)
• The torque changes rotation period and obliquity (YORP effect)
• Both effects are crucial for asteroid orbital evolution

References

• Bottke Jr, W. F., et al. (2006). The Yarkovsky and YORP effects
• Rozitis, B., & Green, S. F. (2012). The influence of rough surface thermal-infrared beaming

update_upper_temperature!(state::SingleAsteroidThermoPhysicalState, i::Integer)

Update the temperature of the upper surface based on the boundary condition state.problem.upper_boundary_condition.

Arguments

• state : Thermophysical simulation state for a single asteroid
• i : Index of the face of the shape model

Base.length(ephem::AbstractAsteroidEphemerides) -> Int

Return the number of timesteps in the ephemerides.

solve(problem, algorithm; ephem, output, initial_temperature_primary, initial_temperature_secondary, show_progress=true) -> solution

Run a thermophysical simulation for a binary asteroid system.

Arguments

• problem : Problem definition (BinaryAsteroidThermoPhysicalProblem)
• algorithm : Numerical method (ExplicitEuler(), ImplicitEuler(), or CrankNicolson())

Keyword Arguments

• ephem : Ephemerides (AbstractBinaryAsteroidEphemerides)
• output : Output specification (BinaryAsteroidOutputSpec); wraps a SingleAsteroidOutputSpec for each body
• initial_temperature_primary : Initial temperature for the primary; Real or AbstractMatrix of size (n_depth, n_face) [K]
• initial_temperature_secondary : Initial temperature for the secondary; Real or AbstractMatrix of size (n_depth, n_face) [K]
• show_progress = true : Display progress meter during simulation

Returns

• BinaryAsteroidThermoPhysicalSolution

Example

T_init = subsolar_temperature(ephem.r_sun[begin], R_vis, ε)
output = BinaryAsteroidOutputSpec(
    SingleAsteroidOutputSpec(output_times, subsurface_face_ids_pri),
    SingleAsteroidOutputSpec(output_times, subsurface_face_ids_sec),
)
solution = solve(problem, CrankNicolson();
    ephem                         = ephem,
    output                        = output,
    initial_temperature_primary   = T_init,
    initial_temperature_secondary = T_init,
)

solve(problem, algorithm; ephem, output, initial_temperature, show_progress=true) -> solution

Run a thermophysical simulation for a single asteroid.

Arguments

• problem : Problem definition (SingleAsteroidThermoPhysicalProblem)
• algorithm : Numerical method (ExplicitEuler(), ImplicitEuler(), or CrankNicolson())

Keyword Arguments

• ephem : Ephemerides (AbstractSingleAsteroidEphemerides)
• output : Output specification (SingleAsteroidOutputSpec); controls which timesteps, face indices, and physical quantities (temperatures, forces, torques) to record
• initial_temperature : Initial temperature; Real for uniform, or AbstractMatrix of size (n_depth, n_face) [K]
• show_progress = true : Display progress meter during simulation

Returns

• SingleAsteroidThermoPhysicalSolution

Example

problem = SingleAsteroidThermoPhysicalProblem(shape, thermo_params;
    with_self_shadowing = true,
    with_self_heating   = true,
)
output = SingleAsteroidOutputSpec(output_times, subsurface_face_ids)
solution = solve(problem, CrankNicolson();
    ephem               = ephem,
    output              = output,
    initial_temperature = 200.0,
)