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 <: AbstractBinaryAsteroidEphemerides

Ephemerides for a binary-asteroid thermophysical simulation. Only temperature is computed; force and torque are not.

Fields

• times : Simulation timesteps [s]
• r_sun : Sun position vector in the primary body-fixed frame [m]
• r_secondary : Secondary position vector in the primary body-fixed frame [m]
• R_primary_to_secondary : Passive rotation matrix from the primary body-fixed frame to the secondary body-fixed frame

struct BinaryAsteroidEphemeridesWithDynamics <: AbstractBinaryAsteroidEphemerides

Ephemerides for a binary-asteroid thermophysical simulation with force and torque output in the inertial frame. Requires the primary body orientation at each timestep.

The secondary-to-inertial rotation can be derived as R_primary_to_inertial[i] * R_primary_to_secondary[i]' when needed.

Fields

• times : Simulation timesteps [s]
• r_sun : Sun position vector in the primary body-fixed frame [m]
• r_secondary : Secondary position vector in the primary body-fixed frame [m]
• R_primary_to_secondary : Passive rotation matrix from the primary body-fixed frame to the secondary body-fixed frame
• R_primary_to_inertial : Passive rotation matrix from the primary body-fixed frame to the inertial frame

BinaryAsteroidEphemeridesWithDynamics(ephem, R_primary_to_inertial)

Upgrade constructor: extend a BinaryAsteroidEphemerides with primary body orientation data.

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
• secondary : Solution for the secondary

Outer constructor of BinaryAsteroidThermoPhysicalSolution

Arguments

• state : Thermophysical simulation state for a binary asteroid
• ephem : Ephemerides
• times_to_save : Timesteps to save temperature (Common to both the primary and the secondary)
• face_ID_pri : Face indices to save subsurface temperature of the primary
• face_ID_sec : Face indices to save subsurface temperature of the secondary

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 <: AbstractSingleAsteroidEphemerides

Ephemerides for a single-asteroid thermophysical simulation. Only temperature is computed; force and torque are not.

Fields

• times : Simulation timesteps [s]
• r_sun : Sun position vector in the body-fixed frame [m]

struct SingleAsteroidEphemeridesWithDynamics <: AbstractSingleAsteroidEphemerides

Ephemerides for a single-asteroid thermophysical simulation with force and torque output in the inertial frame. Requires the body orientation at each timestep.

Fields

• times : Simulation timesteps [s]
• r_sun : Sun position vector in the body-fixed frame [m]
• R_body_to_inertial : Passive rotation matrix from the body-fixed frame to the inertial frame

SingleAsteroidEphemeridesWithDynamics(ephem, R_body_to_inertial)

Upgrade constructor: extend a SingleAsteroidEphemerides with body orientation data.

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 time steps

• times : Timesteps, given the same vector as ephem.times [s]
• E_in : Input energy per second on the whole surface [W]
• E_out : Output energy per second from the whole surface [W]
• force : Thermal force on the asteroid [N]
• torque : Thermal torque on the asteroid [N ⋅ m]

Saved only at the time steps desired by the user

• times_to_save : Timesteps to save temperature and thermal force on every face [s]
• depth_nodes : Depths of the calculation nodes for 1-D heat conduction [m], a vector of size n_depth
• surface_temperature : Surface temperature [K], a matrix in size of (n_face, n_time).
  • n_face : Number of faces
  • n_time : Number of time steps to save surface temperature
• subsurface_temperature : Temperature [K] as a function of depth [m] and time [s], Dict with face ID as key and a matrix (n_depth, n_time) as an entry.
  • n_depth : The number of the depth nodes
  • n_time : The number of time steps to save temperature
• face_forces : Thermal force on every face of the shape model [N], a matrix in size of (n_face, n_time).
  • n_face : Number of faces
  • n_time : Number of time steps to save surface temperature

Outer constructor of SingleAsteroidThermoPhysicalSolution

Arguments

• state : Thermophysical simulation state for a single asteroid
• ephem : Ephemerides
• times_to_save : Timesteps to save temperature
• face_ID : Face indices to save subsurface temperature

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

energy_in(state::SingleAsteroidThermoPhysicalState) -> E_in

Calculate the total energy input rate (power) absorbed by the entire asteroid surface.

Arguments

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

Returns

• E_in::Float64 : Total absorbed power [W]

Calculation

Integrates the absorbed energy flux over all surface facets:

E_in = Σᵢ F_abs,ᵢ × Aᵢ

where:

• Fabs,ᵢ is the absorbed energy flux on facet i (calculated by `absorbedenergy_flux`)
• Aᵢ is the area of facet i

Components

The absorbed energy includes:

1. Direct solar radiation
2. Scattered light from other facets (if self-heating is enabled)
3. Thermal radiation from other facets (if self-heating is enabled)

Usage

This function is typically used to check energy conservation in the model by comparing with energy_out.

See Also

• energy_out for the total emitted power
• absorbed_energy_flux for the flux calculation

energy_out(state::SingleAsteroidThermoPhysicalState) -> E_out

Calculate the total energy output rate (power) emitted by the entire asteroid surface through thermal radiation.

Arguments

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

Returns

• E_out::Float64 : Total emitted power [W]

Calculation

Integrates the thermal emission over all surface facets using the Stefan-Boltzmann law:

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

where:

• εᵢ is the emissivity of facet i
• σ is the Stefan-Boltzmann constant (5.67×10⁻⁸ W/m²/K⁴)
• Tᵢ is the surface temperature of facet i [K]
• Aᵢ is the area of facet i [m²]

Energy Conservation

In thermal equilibrium, Eout should approximately equal Ein (from energy_in). The ratio Eout/Ein is often used as a convergence criterion in thermophysical models.

See Also

• energy_in for the total absorbed power
• update_thermal_force! for thermal recoil effects from this emission

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, result::BinaryAsteroidThermoPhysicalSolution)

Export the result of BinaryAsteroidThermoPhysicalState to CSV files. The output files are saved in the following directory structure:

dirpath
├── primary
│   ├── physical_quantities.csv
│   ├── subsurface_temperature.csv
│   ├── surface_temperature.csv
│   └── thermal_force.csv
└── secondary
    ├── physical_quantities.csv
    ├── subsurface_temperature.csv
    ├── surface_temperature.csv
    └── thermal_force.csv

Arguments

• dirpath : Path to the directory to save CSV files.
• solution : Solution object for a binary asteroid

export_solution(dirpath, result::SingleAsteroidThermoPhysicalSolution)

Export the result of SingleAsteroidThermoPhysicalState to CSV files. The output files are saved in the following directory structure:

dirpath
├── physical_quantities.csv
├── subsurface_temperature.csv
├── surface_temperature.csv
└── thermal_force.csv

Arguments

• dirpath : Path to the directory to save CSV files.
• solution : Solution object for a single asteroid

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]

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::BinaryAsteroidThermoPhysicalSolution, state::BinaryAsteroidThermoPhysicalState, i_time::Integer)

Record the simulation state at time step i_time into solution.

Arguments

• solution : Solution object for a binary asteroid
• state : Simulation state for a binary asteroid
• i_time : Time step

record_timestep!(solution::SingleAsteroidThermoPhysicalSolution, state::SingleAsteroidThermoPhysicalState, i_time::Integer)

Record the simulation state at time step i_time into solution.

Arguments

• solution : Solution object for a single asteroid
• state : Simulation state for a single asteroid
• i_time : Time step to save data

subsolar_temperature(r☉, R_vis, ε) -> Tₛₛ
subsolar_temperature(r☉, params::AbstractThermoParams) -> 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, times_to_save, face_ID_pri, face_ID_sec, T₀_primary, T₀_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)
• times_to_save : Time points at which to save detailed temperature data [s]
• face_ID_pri : Face indices for subsurface profiles of the primary
• face_ID_sec : Face indices for subsurface profiles of the secondary
• T₀_primary : Initial temperature for the primary; Real or AbstractMatrix of size (n_depth, n_face) [K]
• T₀_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₀_primary   = subsolar_temperature(ephem.r_sun[begin], params1)
T₀_secondary = subsolar_temperature(ephem.r_sun[begin], params2)
solution = solve(problem, CrankNicolson();
    ephem         = ephem,
    times_to_save = times_to_save,
    face_ID_pri   = face_ID_pri,
    face_ID_sec   = face_ID_sec,
    T₀_primary    = T₀_primary,
    T₀_secondary  = T₀_secondary,
)

solve(problem, algorithm; ephem, times_to_save, face_ID, T₀, 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)
• times_to_save : Time points at which to save detailed temperature data [s]
• face_ID : Face indices for which to save subsurface temperature profiles
• T₀ : 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,
)
solution = solve(problem, CrankNicolson();
    ephem         = ephem,
    times_to_save = times_to_save,
    face_ID       = face_ID,
    T₀            = 200.0,
)