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 in the inertial frame are computed:

• R = Nothing : no inertial-frame rotation available; per-face forces in the body-fixed frame can still be obtained via save_face_forces.
• R = Vector{SMatrix{3,3,Float64,9}} : forces and torques are rotated to 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 inertial-frame force/torque output is 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, output_times_secondary;
    subsurface_face_ids_primary   = [1, 2],
    subsurface_face_ids_secondary = [3],
    save_surface_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=[1, 2], save_forces=true)
output_secondary = SingleAsteroidOutputSpec(output_times_secondary; subsurface_face_ids=[3])
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.
• Shapes with surface roughness are rejected with an ArgumentError: the eclipse shadowing and the mutual heating are applied to the global faces only, so the sub-faces of a roughness model would ignore both. Surface roughness is supported for single asteroids.

BinaryAsteroidThermoPhysicalProblem(shape, thermo_params, grid_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 : ThermoParams applied to both bodies, or Tuple (thermo_params1, thermo_params2) for per-body settings
• grid_params : GridParams applied to both bodies, or Tuple (grid_params1, grid_params2) for per-body settings

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

grid_params = GridParams(z_max, Δz, n_depth)
problem = BinaryAsteroidThermoPhysicalProblem(
    (shape1, shape2),
    (thermo_params1, thermo_params2),
    grid_params;
    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.

struct GridParams

Numerical grid settings for the 1D heat conduction equation, shared across all facets.

Fields

• z_max : Depth of the lower boundary [m]
• n_depth : Number of depth nodes
• Δz : Depth step width [m]

Notes

Currently assumes a uniform depth grid. Variable-spacing support (e.g., finer near the surface) is planned for a future version.

GridParams(; z_max, n_depth)

Construct GridParams from keyword arguments, with Δz computed automatically as z_max / (n_depth - 1), placing nodes uniformly at 0, Δz, 2Δz, …, z_max.

Keyword Arguments

• z_max : Depth of the lower boundary [m]
• n_depth : Number of depth nodes

Notes

To specify Δz explicitly (advanced use), use the positional constructor GridParams(z_max, n_depth, Δz).

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 RoughnessNeighbours

For one global face that carries a roughness model: which of its sub-faces are seen from the direction of each neighbouring global face, in the order of the face's list in the face visibility graph. This is the geometric part of the exchange of radiation between roughness models (see update_flux_rad_single!); it does not change during a run and is built once, at state construction, when with_self_heating is enabled.

Fields

• visible_sub_faces : visible_sub_faces[p][m] is true when sub-face m is seen from the direction of the p-th visible global face (in the local frame of the model, including self-shadowing by the model's own topography)
• index_in_neighbour : Index of this face in the list of the p-th visible face, so that the neighbour's mask towards this face can be read directly; 0 when the neighbour has no roughness model

struct SingleAsteroidEphemerides{R} <: AbstractSingleAsteroidEphemerides

Ephemerides for a single-asteroid thermophysical simulation.

The type parameter R controls whether force and torque in the inertial frame are computed:

• R = Nothing : no inertial-frame rotation available; per-face forces in the body-fixed frame can still be obtained via save_face_forces.
• R = Vector{SMatrix{3,3,Float64,9}} : forces and torques are rotated to 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 inertial-frame force/torque output is 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. Quantities recorded for every face are switched on and off by a save_* flag; quantities recorded for selected faces only are switched on by listing those faces — an empty list means "not saved".

Fields

• output_times : Timesteps at which data are saved [s]; must be a subset of ephem.times
• save_surface_temperature : Save the surface temperature of every face at output_times (default: true)
• subsurface_face_ids : Faces whose subsurface temperature profiles to save at output_times (default: none)
• roughness_face_ids : Faces whose roughness-model surface temperatures (every sub-face) to save at output_times (default: none); requires a shape with surface roughness
• save_face_forces : Save the thermal force on every face at output_times (default: false)
• save_forces : Save the net thermal force at output_times (default: false)
• save_torques : Save the net thermal torque at output_times (default: false)

Notes

• roughness_face_ids requires a problem built on a shape with surface roughness and a roughness model on every listed face; both are checked when the solution is allocated at solve time.
• save_face_forces stores per-face forces in the body-fixed frame; it works with both SingleAsteroidEphemerides{Nothing} and SingleAsteroidEphemerides{<:AbstractVector}.
• 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;
    save_surface_temperature = true,
    subsurface_face_ids      = [1, 2, 3],
    roughness_face_ids       = [1, 7],
    save_face_forces         = false,
    save_forces              = true,
    save_torques             = true,
)

SingleAsteroidOutputSpec(output_times; kwargs...)

Keyword constructor; see SingleAsteroidOutputSpec for the meaning of each keyword.

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, optionally carrying surface roughness)
• thermo_params : Thermophysical material parameters (all vectors expanded to length n_face)
• grid_params : Numerical grid settings
• 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

thermo_params = ThermoParams(k, ρ, Cₚ, R_vis, R_ir, ε)
grid_params   = GridParams(z_max, n_depth, Δz)
problem = SingleAsteroidThermoPhysicalProblem(shape, thermo_params, grid_params;
    with_self_shadowing = true,
    with_self_heating   = true,
    upper_boundary_condition = RadiationBoundaryCondition(),
    lower_boundary_condition = InsulationBoundaryCondition(),
)
solution = solve(problem, CrankNicolson(); ephem = ephem, initial_temperature = 200.0)

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

Construct a thermophysical problem for a single asteroid.

Arguments

• shape : Shape model of the asteroid
• thermo_params : Thermophysical material parameters (ThermoParams)
• grid_params : Numerical grid settings (GridParams)

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

• The geometric data required by the enabled flags is computed automatically when missing: with_self_shadowing = true needs face_visibility_graph and face_max_elevations, while with_self_heating = true needs face_visibility_graph alone. To avoid this, pass with_face_visibility=true when loading the shape.
• ThermoParams with length-1 vectors is expanded to n_face at construction time.

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
• roughness_surface_temperature : Surface temperature of the sub-faces of the roughness model, by global face ID, each entry (n_sub, n_save) [K], or nothing. Sub-face j is face j of the roughness model attached to that global face

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.
• surface_temperature and subsurface_temperature are those of the global faces, i.e. the smooth-surface baseline of a hierarchical run; roughness_surface_temperature is where the rough-surface temperatures live.

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]

Surface roughness (empty for a smooth surface)

• face_roughness_indices : Maps face index → roughness_states index (0 = no roughness); length = n_face when the shape has roughness, empty otherwise. Mirrors shape.roughness.face_roughness_indices, but maps to independent per-face states rather than shared models.
• roughness_states : Independent sub-face state per roughness-carrying face; empty when the shape has no roughness. Each sub-state is itself a SingleAsteroidThermoPhysicalState with empty roughness (the roughness models are smooth ShapeModels)
• roughness_neighbours : Per roughness-carrying face (same order as roughness_states), the visibility of its sub-faces from the direction of each neighbouring global face (RoughnessNeighbours). Built only when with_self_heating is enabled; empty otherwise

Notes

When the shape carries surface roughness (has_roughness(problem.shape)), the fields above describe the global faces, and each global face with a roughness model additionally has its own sub-face state in roughness_states. Code that iterates the roughness states can simply loop over them: for a smooth surface the vectors are empty and the loop does nothing.

struct ThermoParams <: AbstractThermoParams

Material thermal properties per facet.

Each field is a Vector{Float64} of length n_face. The outer constructor accepts Float64 (uniform) or Vector{Float64} (non-uniform) per field, and scalar arguments are automatically broadcast to match any vector arguments.

Fields

• conductivity : Thermal conductivity for each facet [W/m/K]
• density : Density for each facet [kg/m³]
• heat_capacity : Heat capacity for each facet [J/kg/K]
• reflectance_vis : Reflectance in visible light for each facet [-]
• reflectance_ir : Reflectance in thermal infrared for each facet [-]
• emissivity : Emissivity for each facet [-]

ThermoParams(conductivity, density, heat_capacity, reflectance_vis, reflectance_ir, emissivity)

Construct ThermoParams from scalar or vector arguments, which may be freely mixed. Float64 arguments are automatically broadcast to match the length of any Vector{Float64} arguments. All vector arguments must have the same length.

Arguments

• 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 [-]

Examples

# Uniform surface (all scalars)
ThermoParams(0.1, 1500.0, 800.0, 0.05, 0.0, 0.9)

# Non-uniform conductivity only; other parameters are uniform
ThermoParams(k_vec, 1500.0, 800.0, 0.05, 0.0, 0.9)

# Fully non-uniform
ThermoParams(k_vec, ρ_vec, Cₚ_vec, R_vis_vec, R_ir_vec, ε_vec)

ThermoParams(; conductivity, density, heat_capacity, reflectance_vis, reflectance_ir, emissivity)

Construct ThermoParams from keyword arguments. Each argument can be a Float64 (uniform) or Vector{Float64} (non-uniform), and scalar/vector arguments may be freely mixed.

Keyword Arguments

• 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 [-]

_add_external_radiation!(state, k, i)
_add_external_scattering!(state, k, i)

Add to the sub-faces of the roughness model on global face i (sub-state k) the thermal radiation, respectively the reflected sunlight, that they receive from the other global faces.

For every face j visible from i (view factor f_ij, direction d̂_ij from the face visibility graph):

1. Emitter: the emission of j towards i. If j carries a roughness model, it is the directional emission of that model (_directional_emission) — the hot sunlit wall of a crater that faces i radiates more towards i than a smooth face would (thermal-infrared beaming), a shaded wall less. The sub-faces of j seen from i come from the mask that j precomputed towards i (RoughnessNeighbours). If j is smooth, it radiates as a Lambertian face, ε σ T_j⁴ or R_vis F_sun,j.
2. Far-field: the irradiance reaching face i is F_ij = E_j(d̂_ji) f_ij, the same form as the smooth-face term ε σ T_j⁴ f_ij — the view factor already carries the geometry of the pair, and both faces are taken as small compared to their distance, as the view factor does.
3. Receiver: F_ij is distributed over the sub-faces of i that see j (_add_directional_irradiance!), so the wall facing j is heated and the wall behind it is not.

Thermal emission uses the sub-face temperatures of the previous time step and reflected sunlight uses the direct solar flux only (single scattering, as on the global level), so the result does not depend on the order in which the roughness models are updated. The global-level fluxes of face i are not touched: they remain the smooth-surface baseline.

_directional_emission(patch::ShapeModel, visible, d̂_local, emission) -> E

Emission of a roughness model patch towards the direction d̂_local (unit vector in the local frame of the patch), per unit area of the patch's reference plane projected onto that direction:

E(d̂) = Σₙ Vₙ (n̂ₙ ⋅ d̂)⁺ Eₙ aₙ / (A_proj cos θ)

where visible[n] (Vₙ) tells whether sub-face n is seen from d̂, emission(n) (Eₙ) is the quantity emitted by sub-face n per unit area — ε σ T⁴ for thermal emission, R_vis F_sun for reflected sunlight, a radiance for roughness_radiance — and cos θ is the z component of d̂_local. For a uniform, unshadowed patch this reduces to E = Eₙ: the representative patch radiates like a smooth Lambertian face. Sub-faces facing away from d̂ do not contribute.

The caller guarantees cos θ > 0.

_expand_thermo_params(thermo_params::ThermoParams, n_face::Int) -> ThermoParams

Expand a ThermoParams with length-1 vectors to length n_face (uniform surface), or validate that all vectors already have length n_face (non-uniform surface).

_log_elapsed(f, what::AbstractString)

Run f, reporting what it is about to do and, once it returns, how long it took — both on a single line. Preparing the geometric data for a large shape model takes minutes, so the opening half tells the user what the wait is for and the closing half reports its cost.

The line is written to stdout rather than through @info, because a log record cannot be completed after the fact. This matches how solve already reports its progress.

_prepare_self_heating!(shape::ShapeModel)

Ensure that shape carries the geometric data required for self-heating, building it when missing. Self-heating needs the view factors of face_visibility_graph only; the maximum elevations used to accelerate self-shadowing are not involved.

_prepare_self_shadowing!(shape::ShapeModel)

Ensure that shape carries the geometric data required for self-shadowing, building whatever is missing. face_max_elevations depends on face_visibility_graph, so the visibility graph is built first.

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]

brightness_temperature(problem, solution, i_save, d̂; λ=nothing) -> T_b

Brightness temperature of every global facet towards the observer direction d̂: the temperature of a blackbody whose Lambertian radiance equals the facet's directional_radiance, $B(T_b)/\pi = L$. Emissivity is not divided out, so a smooth grey facet at temperature $T$ has $T_b = \varepsilon^{1/4} T$ for the total radiance.

With λ the spectral radiance at that wavelength is inverted through the Planck function. A facet seen from behind gives NaN. Returns a vector of length n_face [K].

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

directional_radiance(problem, solution, i_save, d̂; λ=nothing) -> L

Thermal radiance of every global facet towards the observer direction d̂ (body-fixed frame) at output time solution.output.output_times[i_save].

Facets whose roughness-model surface temperatures were recorded (output.roughness_face_ids) radiate anisotropically according to roughness_radiance; every other facet radiates as a smooth Lambertian surface, $\varepsilon B(T_i)/\pi$, from its recorded surface_temperature. A facet seen from behind gives NaN.

The result is one value per global facet, ready to be placed on an image by a ray-caster such as FOVSimulator.generate_image_radiance (pass $\varepsilon = 1$ with the corresponding brightness_temperature, or a radiance-taking variant).

Arguments

• problem : The problem that produced solution (shape and emissivities)
• solution : Solution with surface_temperature recorded (and roughness_surface_temperature for rough facets)
• i_save : Index into solution.output.output_times
• d̂ : Direction from the asteroid to the observer in the body-fixed frame

Keyword Arguments

• λ : Wavelength [m] for the spectral radiance; nothing (default) for the total radiance

Returns

• L::Vector{Float64} of length n_face [W/m²/sr], or [W/m²/m/sr] when λ is given

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 ArgumentError if λ ≥ 0.5 (stability violation)

export_solution(dirpath, solution::BinaryAsteroidThermoPhysicalSolution)

Export results for both bodies to dirpath/primary/ and dirpath/secondary/. All directories are created automatically if they do not exist.

export_solution(dirpath, solution::SingleAsteroidThermoPhysicalSolution)

Export simulation results to CSV files in dirpath. dirpath is created automatically if it does not exist.

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.subsurface_face_ids is non-empty
• thermal_face_forces.csv : when output.save_face_forces = true
• thermal_net_forces.csv : when output.save_forces = true or output.save_torques = true
• roughness_surface_temperature.csv : when output.roughness_face_ids is non-empty; long format with columns time, face_id, sub_face_id, temperature

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. When the shape carries surface roughness, each sub-face state is initialized to the surface temperature of its parent face.

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₀, including the sub-face states of any surface roughness (no-op for a smooth surface).

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ᵢ

When the shape carries surface roughness, a face with a roughness model is counted from its sub-faces and not from the global level: its roughness model is a patch that represents the face statistically, so the absorbed power of the patch, Σⱼ F_abs,ⱼ aⱼ in the units of the model, is scaled to the area of the face by Aᵢ / A_proj (see update_thermal_force!). A face without a roughness model contributes as usual.

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.

When the shape carries surface roughness, a face with a roughness model is counted from its sub-faces, scaled to the area of the face by Aᵢ / A_proj, exactly as in integrate_absorbed_power; the smooth-surface emission of that face is not added.

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

roughness_radiance(shape::ShapeModel, i, T_sub, ε, d̂; λ=nothing) -> L

Thermal radiance of facet i of shape, whose roughness model has the sub-facet surface temperatures T_sub, towards the observer direction d̂ given in the body-fixed frame.

The roughness model is a representative patch of the facet's surface, so the radiance is the emission of its visible sub-facets per unit projected area of the patch:

\[L_i(\hat{\mathbf d}) = \frac{1}{A_\mathrm{proj}\,(\hat{\mathbf z}\cdot\hat{\mathbf d}_\mathrm{local})} \sum_j V_j(\hat{\mathbf d}_\mathrm{local})\,(\hat{\mathbf n}_j\cdot\hat{\mathbf d}_\mathrm{local})^+\,a_j\,\frac{\varepsilon\,B(T_j)}{\pi}\]

where $\hat{\mathbf d}_\mathrm{local}$ is d̂ rotated into the local frame of the facet, $V_j$ is 1 when sub-facet j is visible from that direction (not hidden by the crater walls), $A_\mathrm{proj} = \sum_j a_j (\hat{\mathbf n}_j \cdot \hat{\mathbf z})$ is the projected area of the patch, and $B(T) = \sigma T^4$ (total) or the Planck function at wavelength λ (spectral). For an isothermal patch without shadowing this reduces to the Lambertian $\varepsilon B(T)/\pi$; a sunlit crater whose hot wall faces the observer radiates more than that — thermal-infrared beaming.

Arguments

• shape : Shape model with surface roughness; facet i must carry a roughness model
• i : Global facet index
• T_sub : Surface temperature of each sub-facet of the roughness model [K]
• ε : Emissivity of the facet (grey: independent of wavelength)
• d̂ : Direction from the facet to the observer in the body-fixed frame (normalised internally)

Keyword Arguments

• λ : Wavelength [m] for the spectral radiance; nothing (default) for the total radiance

Returns

• L : Radiance [W/m²/sr], or spectral radiance [W/m²/m/sr] when λ is given. NaN when the facet is seen from behind ($\hat{\mathbf z}\cdot\hat{\mathbf d}_\mathrm{local} \le 0$)

Notes

• Sub-facet visibility is evaluated with update_illumination! of the roughness model with the observer direction in place of the Sun: being lit from a direction and being visible from it are the same test. The roughness model's visibility graph and maximum elevations are built on demand if missing.
• The patch has no neighbours, so at grazing angles rays that would be blocked by the next patch are not; the representative-patch picture assumes patches much smaller than the facet.

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. For a shape with surface roughness, these are the global faces; the sub-face temperatures live in state.roughness_states.

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::SingleAsteroidThermoPhysicalState)

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 : Thermophysical simulation state for a single asteroid, with or without surface roughness
• 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 with_self_shadowing and with_self_heating
• When the shape carries surface roughness, each of the three updates handles the global faces first and then the sub-faces of every roughness model. The order of the three calls is fixed: the external irradiation of the sub-faces reads the scattered and thermal flux of the parent global face, which must therefore be complete before the sub-face update runs.

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 the flux of thermal radiation incident on each face from the surrounding surface.

Only the direct emission of the other faces is counted (single bounce); thermal radiation they reflect is neglected. flux_rad is an incident flux, like flux_sun and flux_scat: the thermal-infrared reflectance of the receiving face is applied where the flux is absorbed.

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 with_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 with_self_shadowing = true) also receive zero flux
• When the shape carries surface roughness, the sub-faces of every roughness model are then illuminated in the local frame of their parent face (dark when the parent is not illuminated)

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 : Thermophysical simulation state for a single asteroid, with or without surface roughness
• Δ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

Notes

• For a shape with surface roughness, the global faces are advanced first and then every sub-face state in roughness_states, each a full set of 1D columns driven by its own fluxes and using its own solver cache. The global faces are solved independently of their roughness models and serve as the smooth-surface baseline for the same run; nothing flows from the sub-faces back to them. For a smooth surface roughness_states is empty and only the global faces are advanced.

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.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: their momentum stays with the body, so it cancels part of the recoil. It is the momentum counterpart of the energy re-absorbed in self-heating, and is therefore applied only when with_self_heating is enabled; without it, every emitted photon counts as having left the body.

Outputs (stored in state)

• state.face_forces : Thermal force vector on each facet [N]
• state.force : Net thermal force Σᵢ F_i in the body-fixed frame [N]
• state.torque : Net thermal torque Σᵢ r_i × F_i about the body-fixed origin [N⋅m]; the origin is assumed to be the centre of mass

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

Surface roughness

When the shape carries surface roughness, a face with a roughness model gets its recoil from the sub-faces of its model instead of the smooth-surface baseline: the sub-face forces are summed in the local frame, scaled from the patch to the parent's area by Aᵢ / A_proj, and rotated into the body frame. With with_self_heating, the photons that escape the patch towards the sky and are intercepted by other faces are also accounted for (isotropic patch emission). The torque takes the centre of the parent face as the point of action; the torque of the patch about its own centre is neglected. The sub-face states' own force and torque are not used. Only the derived quantities of the faces are overwritten; their fluxes and temperatures remain the smooth-surface solution.

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=subsurface_face_ids_pri),
    SingleAsteroidOutputSpec(output_times; subsurface_face_ids=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,
)