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

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

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

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

Arguments

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

Returns

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

Mathematical Formula

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

Physical Interpretation

The function accounts for different reflectance properties at different wavelengths:

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

Example

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

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

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

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

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

Arguments

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

Returns

• T : Temperature [K]

blackbody_radiance(λ, T) -> L_λ

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

Arguments

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

Return

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

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

blackbody_radiance(T) -> L

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

Arguments

• T : Temperature [K]

Return

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

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

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

Arguments

• thermo_params : Thermophysical parameters
• shape : Shape model

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

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

Arguments

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

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

Calculate the z-coordinate (depth) of a point on a concave spherical segment (crater).

Arguments

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

Returns

• z::Real: Depth below the reference surface (negative value) [m]

Mathematical Model

The crater is modeled as a spherical cap. For a point (x,y) within the crater:

\[z = R - h - \sqrt{R^2 - d^2}\]

where:

• R is the curvature radius calculated from crater_curvature_radius(r, h)
• d is the distance from the crater center: d = √((x-xc)² + (y-yc)²)
• z = 0 outside the crater (when d > r)

Example

# A crater at (0.5, 0.5) with radius 0.2 and depth 0.05
z = concave_spherical_segment(0.2, 0.05, 0.5, 0.5, 0.6, 0.5)  # Point at crater rim

Notes

• Returns z ≤ 0 (crater is a depression)
• Returns z = 0 for points outside the crater
• The reference surface is at z = 0

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

Generate a grid representation of a concave spherical segment (crater) on a unit square domain.

Arguments

• r::Real : Crater radius (in normalized coordinates, 0-1)
• h::Real : Crater depth (in same units as radius)

Keyword Arguments

• Nx::Integer=32: Number of grid points in the x-direction
• Ny::Integer=32: Number of grid points in the y-direction
• xc::Real=0.5: x-coordinate of crater center (normalized, 0-1)
• yc::Real=0.5: y-coordinate of crater center (normalized, 0-1)

Returns

• xs::LinRange: x-coordinates of grid points (0 to 1)
• ys::LinRange: y-coordinates of grid points (0 to 1)
• zs::Matrix: z-coordinates (depths) at each grid point

Example

# Create a crater grid with radius 0.3 and depth 0.1
xs, ys, zs = concave_spherical_segment(0.3, 0.1; Nx=64, Ny=64)

# Off-center crater
xs, ys, zs = concave_spherical_segment(0.2, 0.05; xc=0.3, yc=0.7)

Notes

• The domain is normalized to [0,1] × [0,1]
• The crater parameters should be scaled accordingly
• Grid has (Nx+1) × (Ny+1) total points
• Useful for generating synthetic rough surfaces

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

Calculate the curvature radius of a concave spherical segment (crater shape).

Arguments

• r::Real: Crater radius [m]
• h::Real: Crater depth [m]

Returns

• R::Real: Curvature radius of the spherical cap [m]

Mathematical Formula

The curvature radius R of a spherical cap is given by:

\[R = \frac{r^2 + h^2}{2h}\]

where r is the crater rim radius and h is the crater depth.

Example

# A crater with 10m radius and 2m depth
R = crater_curvature_radius(10.0, 2.0)  # Returns 26.0

Notes

• The crater is modeled as a spherical cap (portion of a sphere)
• Larger R values correspond to shallower craters
• As h → 0, R → ∞ (flat surface)

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

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

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]

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

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☉, R_vis, ε) -> Tₛₛ
subsolar_temperature(r☉, params::AbstractThermoParams) -> Tₛₛ

Calculate the subsolar temperature on an asteroid at a given heliocentric distance.

Arguments

• r☉::Vector : Sun's position vector in the asteroid's fixed frame [m]
• R_vis::Real : Visible light reflectance (albedo) [-]
• ε::Real : Emissivity [-]
• params::AbstractThermoParams : Thermal parameters (alternative input)

Returns

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

Mathematical Formula

Assuming radiative equilibrium with zero thermal conductivity (zero thermal inertia):

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

where:

• $\Phi_\odot = \Phi_0 / r^2$ is the solar flux at distance r [au]
• $\Phi_0 = 1366$ W/m² is the solar constant at 1 au
• $\sigma$ is the Stefan-Boltzmann constant

Notes

• This gives the maximum temperature for a non-rotating asteroid
• Actual subsolar temperature may be lower due to thermal inertia
• Valid for airless bodies with negligible heat conduction

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

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!(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})

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

Arguments

• btpm::BinaryAsteroidTPM : Thermophysical model for a binary asteroid
• r☉₁::StaticVector{3} : Sun's position vector in the primary's body-fixed frame [m]
• r☉₂::StaticVector{3} : Sun's position vector in the secondary's body-fixed frame [m]

Notes

• Each component's solar flux is calculated independently using its own Sun vector.
• The vectors should account for the different orientations of the two bodies.
• Mutual shadowing between components is handled separately by mutual_shadowing!.

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

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

Arguments

• stpm::SingleAsteroidTPM : Thermophysical model for a single asteroid
• r̂☉::StaticVector{3} : Unit vector pointing toward the Sun in body-fixed frame
• F☉::Real : Solar flux at the asteroid's location [W/m²]

Algorithm

For each face, the solar flux is calculated as:

F_sun = F☉ × max(0, n̂ · r̂☉)

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

Notes

• Faces with negative dot product (facing away from Sun) receive zero flux
• Shadowed faces (when SELF_SHADOWING = true) also receive zero flux
• The input solar direction r̂☉ is normalized internally for safety

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

Update solar irradiation flux on every face using the Sun's position vector.

Arguments

• stpm::SingleAsteroidTPM : Thermophysical model for a single asteroid
• r☉::StaticVector{3} : Position vector from asteroid to Sun in body-fixed frame [m]

Algorithm

1. Calculates the solar flux using the inverse square law:
2. Normalizes the Sun direction vector
3. Calls the main update_flux_sun! function with computed values

Notes

• The input vector r☉ should be in meters
• Solar flux is automatically computed from the solar constant and distance
• This is a convenience function that handles flux calculation

update_lower_temperature!(stpm::SingleAsteroidTPM)

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

Arguments

• stpm::SingleAsteroidTPM : Thermophysical model 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] = stpm.BC_LOWER.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!(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