API Reference

Abstract type of an asteroid's thermophysical model. The AbstractAsteroidTPM type is an alias for AbstractAsteroidThermoPhysicalModel.

struct BinaryAsteroidThermoPhysicalModel{M1, M2} <: AbstractAsteroidThermoPhysicalModel

Fields

• pri : TPM for the primary
• sec : TPM for the secondary
• MUTUAL_SHADOWING : Flag to consider mutual shadowing
• MUTUAL_HEATING : Flag to consider mutual heating

BinaryAsteroidThermoPhysicalModel(pri, sec; MUTUAL_SHADOWING=true, MUTUAL_HEATING=true) -> btpm

Construct a thermophysical model for a binary asteroid (BinaryAsteroidThermoPhysicalModel).

struct BinaryAsteroidThermoPhysicalModelResult

Output data format for BinaryAsteroidThermoPhysicalModel

Fields

• pri : TPM result for the primary
• sec : TPM result for the secondary

Outer constructor of BinaryAsteroidThermoPhysicalModelResult

Arguments

• btpm : Thermophysical model 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

Abstract type of a boundary condition for a heat conduction equation

Type of the Crank-Nicolson method:

• Implicit in time (Unconditionally stable in the heat conduction equation)
• Second order in time

The CrankNicolsonSolver type has vectors for the tridiagonal matrix algorithm.

References

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

Type of the explicit (forward) Euler method:

• Explicit in time
• First order in time

The ExplicitEulerSolver type includes a vector for the temperature at the next time step.

Abstract type of a solver for a heat conduction equation

Type of the implicit (backward) Euler method:

• Implicit in time (Unconditionally stable in the heat conduction equation)
• First order in time

The ImplicitEulerSolver type has vectors for the tridiagonal matrix algorithm.

Singleton type of insulation boundary condition

Type of isothermal boundary condition

Singleton type of radiation boundary condition

ShapeModel

A polyhedral shape model of an asteroid.

Fields

• nodes : Vector of node positions
• faces : Vector of vertex indices of faces
• face_centers : Center position of each face
• face_normals : Normal vector of each face
• face_areas : Area of of each face
• visiblefacets : Vector of vector of VisibleFacet

struct SingleAsteroidThermoPhysicalModel <: AbstractAsteroidThermoPhysicalModel

Fields

• shape : Shape model

• thermo_params : Thermophysical parameters

• flux_sun : Flux of direct sunlight on each face [W/m²]

• flux_scat : Flux of scattered light on each face [W/m²]

• flux_rad : Flux of thermal emission from surrounding surface on each face [W/m²]

• temperature : Temperature matrix (n_depth, n_face) according to the number of depth cells n_depth and the number of faces n_face.

• face_forces : Thermal force on each face

• force : Thermal recoil force at body-fixed frame (Yarkovsky effect)

• torque : Thermal recoil torque at body-fixed frame (YORP effect)

• SELF_SHADOWING : Flag to consider self-shadowing

• SELF_HEATING : Flag to consider self-heating

• SOLVER : Solver of heat conduction equation

• BC_UPPER : Boundary condition at the upper boundary

• BC_LOWER : Boundary condition at the lower boundary

TODO:

• roughness_maps ::ShapeModel[]

SingleAsteroidThermoPhysicalModel(shape, thermo_params; SELF_SHADOWING=true, SELF_HEATING=true) -> stpm

Construct a thermophysical model for a single asteroid (SingleAsteroidThermoPhysicalModel).

Arguments

• shape : Shape model
• thermo_params : Thermophysical parameters

Keyword arguments

• SELF_SHADOWING : Flag to consider self-shadowing
• SELF_HEATING : Flag to consider self-heating
• SOLVER : Solver of heat conduction equation
• BC_UPPER : Boundary condition at the upper boundary
• BC_LOWER : Boundary condition at the lower boundary

struct SingleAsteroidThermoPhysicalModelResult

Output data format for SingleAsteroidThermoPhysicalModel

Fields

Saved at all time steps

• times : Timesteps, given the same vector as ephem.time [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 SingleAsteroidThermoPhysicalModelResult

Arguments

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

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

struct VisibleFacet

Index of an interfacing facet and its view factor

Fields

• id : Index of the interfacing facet
• f : View factor from facet i to j
• d : Distance from facet i to j
• d̂ : Normal vector from facet i to j

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, 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

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

concave_spherical_segment(r, h, xc, yc, x, y) -> z

Return the z-coordinate of a concave spherical segment.

Arguments

• r : Crater radius
• h : Crater depth
• xc: x-coordinate of crater center
• yc: y-coordinate of crater center
• x : x-coordinate where to calculate z
• y : y-coordinate where to calculate z

concave_spherical_segment(r, h; Nx=2^5, Ny=2^5, xc=0.5, yc=0.5) -> xs, ys, zs

Return (x, y, z) grid of a concave spherical segment.

Arguments

• r : Crater radius
• h : Crater depth
• Nx: Number of nodes in the x-direction
• Ny: Number of nodes in the y-direction
• xc: x-coordinate of crater center
• yc: y-coordinate of crater center

crank_nicolson!(stpm::SingleAsteroidTPM, Δ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

• stpm::SingleAsteroidTPM : Thermophysical model 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

crater_curvature_radius(r, h) -> R

Return the curvature radius of a concave spherical segment.

Arguments

• r: Crater radius
• h: Crater depth

energy_in(stpm::SingleAsteroidTPM) -> E_in

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

Arguments

• stpm::SingleAsteroidTPM : Thermophysical model 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(stpm::SingleAsteroidTPM) -> E_out

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

Arguments

• stpm::SingleAsteroidTPM : Thermophysical model 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

equivalent_radius(VOLUME::Real) -> R_eq
equivalent_radius(shape::ShapeModel) -> R_eq

Calculate the radius of a sphere with the same volume as the given volume or shape.

Arguments

• VOLUME::Real : Volume of the object [unit³]
• shape::ShapeModel : Shape model to compute equivalent radius for

Returns

• R_eq::Float64 : Equivalent spherical radius [unit]

Formula

R_eq = (3V/4π)^(1/3)

Example

# For a cube with side length 2
volume = 8.0  # 2³
R_eq = equivalent_radius(volume)  # Returns ≈ 1.24

# For a shape model
shape = load_shape_obj("asteroid.obj")
R_eq = equivalent_radius(shape)

See Also

• polyhedron_volume for volume calculation
• maximum_radius, minimum_radius for other size metrics

explicit_euler!(stpm::SingleAsteroidTPM, Δ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

• stpm::SingleAsteroidTPM : Thermophysical model 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_TPM_results(dirpath, result::BinaryAsteroidThermoPhysicalModelResult)

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

dirpath
├── pri
│   ├── physical_quantities.csv
│   ├── subsurface_temperature.csv
│   ├── surface_temperature.csv
│   └── thermal_force.csv
└── sec
    ├── physical_quantities.csv
    ├── subsurface_temperature.csv
    ├── surface_temperature.csv
    └── thermal_force.csv

Arguments

• dirpath : Path to the directory to save CSV files.
• result : Output data format for BinaryAsteroidThermoPhysicalModel

export_TPM_results(dirpath, result::SingleAsteroidThermoPhysicalModelResult)

Export the result of SingleAsteroidThermoPhysicalModel 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.
• result : Output data format for SingleAsteroidThermoPhysicalModel

face_area(vs::StaticVector{3, <:StaticVector{3}}) -> area
face_area(v1::StaticVector{3}, v2::StaticVector{3}, v3::StaticVector{3}) -> area

Calculate the area of a triangular face using the cross product method.

Arguments

• vs : A static vector containing three vertices of the triangle
• v1, v2, v3 : Three vertices of the triangle

Returns

• area::Float64 : The area of the triangle [unit²]

Mathematical Formula

The area is calculated as: A = ||(v2 - v1) × (v3 - v2)|| / 2

Example

v1 = SVector(0.0, 0.0, 0.0)
v2 = SVector(1.0, 0.0, 0.0)
v3 = SVector(0.0, 1.0, 0.0)
area = face_area(v1, v2, v3)  # Returns 0.5

face_center(vs::StaticVector{3, <:StaticVector{3}}) -> center
face_center(v1::StaticVector{3}, v2::StaticVector{3}, v3::StaticVector{3}) -> center

Calculate the centroid (center) of a triangular face.

Arguments

• vs : A static vector containing three vertices of the triangle
• v1, v2, v3 : Three vertices of the triangle

Returns

• center::StaticVector{3} : The centroid position vector of the triangle

Example

v1 = SVector(0.0, 0.0, 0.0)
v2 = SVector(1.0, 0.0, 0.0)
v3 = SVector(0.0, 1.0, 0.0)
center = face_center(v1, v2, v3)  # Returns SVector(1/3, 1/3, 0.0)

face_normal(vs::StaticVector{3, <:StaticVector{3}}) -> n̂
face_normal(v1::StaticVector{3}, v2::StaticVector{3}, v3::StaticVector{3}) -> n̂

Calculate the unit normal vector of a triangular face using the right-hand rule. The normal direction follows the ordering of vertices: v1 → v2 → v3.

Arguments

• vs : A static vector containing three vertices of the triangle
• v1, v2, v3 : Three vertices of the triangle in counter-clockwise order

Returns

• n̂::StaticVector{3} : The unit normal vector of the triangle

Mathematical Formula

The normal vector is calculated as: n̂ = normalize((v2 - v1) × (v3 - v2))

Example

v1 = SVector(0.0, 0.0, 0.0)
v2 = SVector(1.0, 0.0, 0.0)
v3 = SVector(0.0, 1.0, 0.0)
n̂ = face_normal(v1, v2, v3)  # Returns SVector(0.0, 0.0, 1.0)

find_visiblefacets!(shape::ShapeModel; show_progress=true)

Find facets that is visible from the facet where the observer is located.

Arguments

• shape : Shape model of an asteroid

Keyword arguments

• show_progress : Switch to show a progress meter

grid_to_faces(xs::AbstractVector, ys::AbstractVector, zs::AbstractMatrix) -> nodes, faces

Convert a regular grid (x, y) and corresponding z-coordinates to triangular facets

| ⧹| ⧹| ⧹|

j+1 ・–C–D–・ |⧹ |⧹ |⧹ | | ⧹| ⧹| ⧹| j ・–A–B–・ |⧹ |⧹ |⧹ | i i+1

Arguments

• xs::AbstractVector : x-coordinates of grid points (should be sorted)
• ys::AbstractVector : y-coordinates of grid points (should be sorted)
• zs::AbstractMatrix : z-coordinates of grid points

implicit_euler!(stpm::SingleAsteroidTPM, Δ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

• stpm::SingleAsteroidTPM : Thermophysical model 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!(btpm::BinaryTBinaryAsteroidThermoPhysicalModelPM, T₀::Real)

Initialize all temperature cells at the given temperature T₀

Arguments

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

init_temperature!(stpm::SingleAsteroidThermoPhysicalModel, T₀::Real)

Initialize all temperature cells at the given temperature T₀

Arguments

• stpm : Thermophysical model for a single asteroid
• T₀ : Initial temperature of all cells [K]

isilluminated(shape::ShapeModel, r☉::StaticVector{3}, i::Integer) -> Bool

Return if the i-th face of the shape model is illuminated by the direct sunlight or not

Arguments

• shape : Shape model of an asteroid
• r☉ : Sun's position in the asteroid-fixed frame, which doesn't have to be normalized.
• i : Index of the face to be checked

load_shape_grid(xs, ys, zs; scale=1.0, find_visible_facets=false) -> shape

Convert a regular grid (x, y) to a shape model

Arguments

• xs::AbstractVector : x-coordinates of grid points
• ys::AbstractVector : y-coordinates of grid points
• zs::AbstractMatrix : z-coordinates of grid points

load_shape_obj(shapepath; scale=1.0, find_visible_facets=false; show_progress=true)

Load a shape model from a Wavefront OBJ file.

Arguments

• shapepath : Path to a Wavefront OBJ file

Keyword arguments

• scale : Scale factor of the shape model
• find_visible_facets : Switch to find visible facets
• show_progress : Switch to show a progress meter

loadobj(shapepath::String; scale=1, message=true) -> nodes, faces

maximum_radius(nodes) -> R_max
maximum_radius(shape::ShapeModel) -> R_max

Find the maximum distance from the origin to any vertex in the shape.

Arguments

• nodes::Vector{SVector{3}} : Vector of vertex positions
• shape::ShapeModel : Shape model

Returns

• R_max::Float64 : Maximum radius [unit]

Notes

• Assumes the shape is centered at the origin
• Useful for determining bounding spheres
• For binary systems, ensure each component is centered appropriately

See Also

• minimum_radius for the minimum distance
• equivalent_radius for volume-based radius

minimum_radius(nodes) -> R_min
minimum_radius(shape::ShapeModel) -> R_min

Find the minimum distance from the origin to any vertex in the shape.

Arguments

• nodes::Vector{SVector{3}} : Vector of vertex positions
• shape::ShapeModel : Shape model

Returns

• R_min::Float64 : Minimum radius [unit]

Notes

• Assumes the shape is centered at the origin
• For highly irregular shapes, this may be much smaller than the equivalent radius
• Useful for determining the deepest surface features

See Also

• maximum_radius for the maximum distance
• equivalent_radius for volume-based radius

mutual_heating!(btpm::BinaryAsteroidTPM, rₛ, R₂₁)

Calculate the mutual heating between the primary and secondary asteroids.

Arguments

• btpm : Thermophysical model for a binary asteroid
• rₛ : Position of the secondary relative to the primary (NOT normalized)
• R₂₁ : Rotation matrix from secondary to primary

TODO

• Need to consider local horizon?

mutual_shadowing!(btpm::BinaryAsteroidTPM, r☉, rₛ, R₂₁)

Detect eclipse events between the primary and secondary, and update the solar fluxes of the faces.

Arguments

• btpm : Thermophysical model for a binary asteroid
• r☉ : Position of the sun relative to the primary (NOT normalized)
• rₛ : Position of the secondary relative to the primary (NOT normalized)
• R₂₁ : Rotation matrix from secondary to primary

polyhedron_volume(nodes, faces) -> vol
polyhedron_volume(shape::ShapeModel) -> vol

Calculate the volume of a closed polyhedral shape using the divergence theorem.

Arguments

• nodes::Vector{SVector{3}} : Vector of vertex positions
• faces::Vector{SVector{3,Int}} : Vector of triangular face definitions (vertex indices)
• shape::ShapeModel : Complete shape model containing nodes and faces

Returns

• vol::Float64 : Volume of the polyhedron [unit³]

Algorithm

Uses the divergence theorem to compute volume from surface triangles:

V = (1/6) Σᵢ (Aᵢ × Bᵢ) ⋅ Cᵢ

where Aᵢ, Bᵢ, Cᵢ are the three vertices of triangle i.

Requirements

• The polyhedron must be closed (watertight)
• Face normals should point outward for positive volume
• Faces must be consistently oriented (all clockwise or all counter-clockwise)

Example

# Create a tetrahedron
nodes = [SVector(0.0, 0.0, 0.0),
         SVector(1.0, 0.0, 0.0),
         SVector(0.0, 1.0, 0.0),
         SVector(0.0, 0.0, 1.0)]
faces = [SVector(1, 2, 3),
         SVector(1, 2, 4),
         SVector(1, 3, 4),
         SVector(2, 3, 4)]
vol = polyhedron_volume(nodes, faces)  # Returns 1/6

See Also

• equivalent_radius to convert volume to equivalent spherical radius

raycast(A, B, C, R) -> Bool

Detect intersection between a ray and a triangle using the Möller-Trumbore algorithm. The ray starts from the origin (0, 0, 0) and extends in direction R.

Arguments

• A::StaticVector{3} : First vertex of the triangle
• B::StaticVector{3} : Second vertex of the triangle
• C::StaticVector{3} : Third vertex of the triangle
• R::StaticVector{3} : Direction vector of the ray (does not need to be normalized)

Returns

• Bool : true if the ray intersects the triangle, false otherwise

Algorithm

Uses the Möller-Trumbore ray-triangle intersection algorithm which:

1. Computes the intersection point using barycentric coordinates
2. Checks if the intersection point lies within the triangle
3. Ensures the intersection occurs in the positive ray direction (t > 0)

Example

A = SVector(1.0, 0.0, 0.0)
B = SVector(0.0, 1.0, 0.0)
C = SVector(0.0, 0.0, 1.0)
R = SVector(1.0, 1.0, 1.0)  # Ray direction
intersects = raycast(A, B, C, R)

References

• Möller, T., & Trumbore, B. (1997). Fast, minimum storage ray-triangle intersection

raycast(A, B, C, R, O) -> Bool

Detect intersection between a ray and a triangle, with arbitrary ray origin. This is a convenience wrapper that translates the triangle to place the ray origin at (0,0,0).

Arguments

• A::StaticVector{3} : First vertex of the triangle
• B::StaticVector{3} : Second vertex of the triangle
• C::StaticVector{3} : Third vertex of the triangle
• R::StaticVector{3} : Direction vector of the ray (does not need to be normalized)
• O::StaticVector{3} : Origin point of the ray

Returns

• Bool : true if the ray intersects the triangle, false otherwise

Example

A = SVector(1.0, 0.0, 0.0)
B = SVector(0.0, 1.0, 0.0)
C = SVector(0.0, 0.0, 1.0)
O = SVector(0.5, 0.5, -1.0)  # Ray origin
R = SVector(0.0, 0.0, 1.0)   # Ray direction
intersects = raycast(A, B, C, R, O)

run_TPM!(btpm::BinaryAsteroidThermoPhysicalModel, ephem, savepath)

Run TPM for a binary asteroid.

Arguments

• btpm : Thermophysical model for a binary asteroid
• ephem : Ephemerides
  • time : Ephemeris times
  • sun1 : Sun's position in the primary's frame
  • sun2 : Sun's position in the secondary's frame
  • sec : Secondary's position in the primary's frame
  • P2S : Rotation matrix from primary to secondary frames
  • S2P : Rotation matrix from secondary to primary frames
• times_to_save : Timesteps to save temperature
• face_ID_pri : Face indices where to save subsurface termperature for the primary
• face_ID_sec : Face indices where to save subsurface termperature for the secondary

Keyword arguments

• show_progress : Flag to show the progress meter

run_TPM!(stpm::SingleAsteroidThermoPhysicalModel, ephem, times_to_save, face_ID; show_progress=true) -> result

Execute the thermophysical model simulation for a single asteroid over the specified time period.

Arguments

• stpm::SingleAsteroidThermoPhysicalModel : Thermophysical model containing shape, thermal parameters, and state
• ephem : Ephemerides data structure containing:
  • time::Vector{Float64} : Time points for the simulation [s]
  • sun::Vector{SVector{3}} : Sun's position vectors in the asteroid-fixed frame (not normalized) [m]
• times_to_save::Vector{Float64} : Specific time points at which to save detailed temperature data [s]
• face_ID::Vector{Int} : Face indices for which to save subsurface temperature profiles

Keyword Arguments

• show_progress::Bool=true : Display progress meter during simulation

Returns

• result::SingleAsteroidThermoPhysicalModelResult : Structure containing:
  • Time series of energy balance (Ein, Eout)
  • Thermal forces and torques at each time step
  • Surface temperatures at specified save times
  • Subsurface temperature profiles for selected faces

Algorithm

1. For each time step:
   • Update solar flux based on sun position
   • Calculate scattered light flux (if self-heating enabled)
   • Calculate thermal radiation flux (if self-heating enabled)
   • Compute thermal forces and torques
   • Save results if at a save point
   • Update temperature distribution for next step

Example

# Setup ephemerides
ephem = (
    time = collect(0:60:86400),  # One day, 1-minute steps
    sun = [SVector(au2m, 0.0, 0.0) for _ in 1:1441]  # Sun at 1 AU
)

# Run simulation
times_to_save = [0.0, 21600.0, 43200.0, 64800.0, 86400.0]  # Every 6 hours
face_ID = [1, 100, 500]  # Save subsurface data for these faces
result = run_TPM!(stpm, ephem, times_to_save, face_ID)

# Check energy balance
println("Final E_out/E_in ratio: ", result.E_out[end]/result.E_in[end])

Performance Notes

• Computation time scales with number of faces and time steps
• Self-shadowing and self-heating calculations add significant overhead
• Consider using larger time steps if thermal inertia is low

See Also

• init_temperature! to set initial conditions
• export_TPM_results to save results to files
• update_temperature! for the core temperature update algorithm

subsolar_temperature(r☉) -> Tₛₛ

Subsolar temperature [K] on an asteroid at a heliocentric distance r☉ [m], assuming radiative equilibrium with zero conductivity.

surface_temperature(stpm::SingleAsteroidThermoPhysicalModel) -> T_surface

Extract the surface temperature values for all faces of the asteroid.

Arguments

• stpm::SingleAsteroidThermoPhysicalModel : Thermophysical model for a single asteroid

Returns

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

Notes

• Returns temperatures from the uppermost layer (index 1) of the temperature matrix
• The length of the returned vector equals the number of faces in the shape model
• Useful for visualization, thermal emission calculations, and analysis

Example

T_surf = surface_temperature(stpm)
T_mean = mean(T_surf)  # Average surface temperature
T_max = maximum(T_surf)  # Hottest point
T_min = minimum(T_surf)  # Coldest point

See Also

• stpm.temperature for the full temperature matrix including subsurface

thermal_diffusivity(k, ρ, Cp) -> α

Arguments

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

Return

• α : Thermal diffusivity [m²/s]

thermal_inertia(k, ρ, Cp) -> Γ

Arguments

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

Return

• Γ : Thermal inertia [tiu (thermal intertia 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π

Arguments

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

Return

• l_2π : Thermal skin depth [m], as defined in Rozitis & Green (2011).

tridiagonal_matrix_algorithm!(a, b, c, d, x)
tridiagonal_matrix_algorithm!(stpm::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_TPM_result!(result::BinaryAsteroidThermoPhysicalModelResult, btpm::BinaryAsteroidThermoPhysicalModel, ephem, i_time::Integer)

Save the results of TPM at the time step i_time to result.

Arguments

• result : Output data format for BinaryAsteroidThermoPhysicalModel
• btpm : Thermophysical model for a binary asteroid
• ephem : Ephemerides
• i_time : Time step

update_TPM_result!(result::SingleAsteroidThermoPhysicalModelResult, stpm::SingleAsteroidThermoPhysicalModel, i_time::Integer)

Save the results of TPM at the time step i_time to result.

Arguments

• result : Output data format for SingleAsteroidThermoPhysicalModel
• stpm : Thermophysical model for a single asteroid
• i_time : Time step to save data

update_flux_rad_single!(btpm::BinaryAsteroidTPM)

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

• btpm : Thermophysical model for a binary asteroid

update_flux_rad_single!(stpm::SingleAsteroidTPM)

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

• stpm : Thermophysical model for a single asteroid

update_flux_scat_single!(btpm::BinaryAsteroidTPM)

Update flux of scattered sunlight, only considering single scattering.

Arguments

• btpm : Thermophysical model for a binary asteroid

update_flux_scat_single!(stpm::SingleAsteroidTPM)

Update flux of scattered sunlight, only considering single scattering.

Arguments

• stpm : Thermophysical model for a single asteroid

update_flux_sun!(btpm::BinaryAsteroidTPM, r☉₁::StaticVector{3}, r☉₂::StaticVector{3})

Arguments

• btpm : Thermophysical model for a binary asteroid
• r☉₁ : Sun's position in the body-fixed frame of the primary, which is not normalized.
• r☉₂ : Sun's position in the body-fixed frame of the secondary, which is not normalized.

update_flux_sun!(stpm::SingleAsteroidTPM, r̂☉::StaticVector{3}, F☉::Real)

Update solar irradiation flux on every face of a shape model.

• shape : Shape model
• r̂☉ : Normalized vector indicating the direction of the sun in the body-fixed frame
• F☉ : Solar radiation flux [W/m²]

update_flux_sun!(stpm::SingleAsteroidTPM, r☉::StaticVector{3})

Update solar irradiation flux on every face of a shape model.

Arguments

• stpm : Thermophysical model for a single asteroid
• r☉ : Position of the sun in the body-fixed frame (NOT normalized)

update_bottom_temperature!(shape::ShapeModel)

Update the temperature of the bottom surface based on the boundary condition stpm.BC_LOWER.

Arguments

• stpm : Thermophysical model for a single asteroid

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!(btpm::BinaryAsteroidTPM, Δt)

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

Arguments

• btpm : Thermophysical model for a binary asteroid
• Δt : Time step [sec]

update_temperature!(stpm::SingleAsteroidTPM, Δt)

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

Arguments

• stpm::SingleAsteroidTPM : Thermophysical model for a single asteroid
• Δt::Real : Time step [s]

Solver Selection

The function automatically selects the appropriate solver based on stpm.SOLVER:

• ExplicitEulerSolver: Forward Euler method (conditionally stable, requires λ < 0.5)
• ImplicitEulerSolver: Backward Euler method (unconditionally stable)
• CrankNicolsonSolver: 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!(stpm::SingleAsteroidTPM)

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

• stpm::SingleAsteroidTPM : Thermophysical model 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 stpm.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!(btpm::BinaryAsteroidTPM)

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

Arguments

• btpm : Thermophysical model for a binary asteroid

update_thermal_force!(stpm::SingleAsteroidTPM)

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

• stpm::SingleAsteroidTPM : Thermophysical model 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 stpm)

• stpm.face_forces : Thermal force vector on each facet [N]
• stpm.force : Total thermal force in body-fixed frame [N]
• stpm.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!(stpm::SingleAsteroidTPM, i::Integer)

Update the temperature of the upper surface based on the boundary condition stpm.BC_UPPER.

Arguments

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

view_factor(cᵢ, cⱼ, n̂ᵢ, n̂ⱼ, aⱼ) -> fᵢⱼ, dᵢⱼ, d̂ᵢⱼ

Calculate the view factor from facet i to facet j, assuming Lambertian emission. The view factor represents the fraction of radiation leaving surface i that directly reaches surface j.

Arguments

• cᵢ::StaticVector{3} : Center position of facet i
• cⱼ::StaticVector{3} : Center position of facet j
• n̂ᵢ::StaticVector{3} : Unit normal vector of facet i
• n̂ⱼ::StaticVector{3} : Unit normal vector of facet j
• aⱼ::Float64 : Area of facet j [m²]

Returns

• fᵢⱼ::Float64 : View factor from facet i to facet j [-]
• dᵢⱼ::Float64 : Distance between facet centers [m]
• d̂ᵢⱼ::StaticVector{3} : Unit direction vector from facet i to facet j

Mathematical Formula

The view factor for differential areas is calculated as:

fᵢⱼ = (cosθᵢ × cosθⱼ) / (π × dᵢⱼ²) × aⱼ

where:

• θᵢ is the angle between n̂ᵢ and the line connecting the centers
• θⱼ is the angle between n̂ⱼ and the line connecting the centers

Diagram

(i)   fᵢⱼ   (j)
 △    -->    △
 cᵢ          cⱼ  : Center of each face
 n̂ᵢ          n̂ⱼ  : Normal vector of each face

References

• Howell, J. R., et al. (2010). Thermal Radiation Heat Transfer