AsteroidThermoPhysicalModels.jl

AsteroidThermoPhysicalModels.jl is a comprehensive Julia-based toolkit for thermophysical modeling (TPM) of asteroids. It allows you to simulate the temperature distribution of asteroids and predict non-gravitational perturbations on their dynamics (Yarkovsky and YORP effects).

📚 Documentation

For detailed documentation, please visit:

• Stable Documentation
• Development Documentation

Sample notebooks are available in Astroshaper-examples.

🚀 Installation

using Pkg
Pkg.add("AsteroidThermoPhysicalModels")
using AsteroidThermoPhysicalModels

Or in the Julia REPL package mode:

julia> ]  # Press ] to enter package mode
pkg> add AsteroidThermoPhysicalModels

🔍 Features

Thermophysical Processes

• Heat Conduction: 1-dimensional heat diffusion in depth direction
  • Multiple numerical solvers available (explicit Euler, implicit Euler, and Crank-Nicolson methods)
• Self-Shadowing: Local shadows cast by topography
• Self-Heating: Re-absorption of scattered and radiated photons by surrounding facets
• Binary Systems: Support for mutual shadowing (eclipses) and mutual heating between primary and secondary bodies
• Surface Roughness: Facets carry roughness models (e.g. spherical craters) solved as thermophysical models of their own, with self-shadowing and self-heating inside the roughness model — thermal-infrared beaming in the forces and in the radiance

Shape Models

• Supports Wavefront OBJ format (*.obj)

Non-Gravitational Effects

• Yarkovsky Effect: Orbital perturbation due to asymmetric thermal emission
• YORP Effect: Rotational perturbation due to asymmetric thermal emission

Observables

• Direction-dependent radiance and brightness temperature of every facet towards an observer, total or at a given wavelength, for comparison with thermal-infrared images

🆕 What's New in v0.3.x

v0.3.0

This release adds surface roughness on top of AsteroidShapeModels.jl v0.6: a facet can carry a small roughness model (e.g. a spherical crater) that is solved as a thermophysical model of its own, which produces thermal-infrared beaming in the recorded forces and in the new direction-dependent radiance.

shape  = load_shape_obj("shape.obj"; scale=1000)
crater = create_shape_crater(0.4, 0.1; Nx=8, Ny=8)   # radius 0.4, depth 0.1, on a 1 × 1 patch
add_roughness_models!(shape, crater)                  # every facet, or add_roughness_models!(shape, crater, face_idx)

problem  = SingleAsteroidThermoPhysicalProblem(shape, thermo_params, grid_params; with_self_shadowing=true, with_self_heating=true)
output   = SingleAsteroidOutputSpec(output_times; roughness_face_ids=1:length(shape.faces))
solution = solve(problem, CrankNicolson(); ephem=ephem, output=output, initial_temperature=200.0)

T_b = brightness_temperature(problem, solution, i_save, d̂_observer; λ=10e-6)   # per facet, towards the observer [K]

Breaking changes (see the Migration Guide):

• ThermoParams holds material properties only; the depth grid moves to the new GridParams, a third positional argument of the problem constructors
• SingleAsteroidOutputSpec(output_times; subsurface_face_ids, roughness_face_ids, save_*...) is keyword-only; face-specific outputs are switched on by listing faces
• Requires AsteroidShapeModels.jl v0.6 (HierarchicalShapeModel no longer exists)
• Input errors raise ArgumentError

Results that change: the net thermal force on non-spherical shapes was biased by a radial projection up to v0.2.1 (about 7 % in magnitude and 6° in direction on Ryugu); recompute archived net forces with v0.3.0.

🌟 Example

Temperature distribution of asteroid Didymos and its satellite Dimorphos:

📖 Basic Usage

The workflow follows a Problem → Solve → Export pattern.

using AsteroidShapeModels
using AsteroidThermoPhysicalModels

# --- Shape model ---
shape = load_shape_obj("path/to/shape.obj"; scale=1000, with_face_visibility=true, with_bvh=true)

# --- Ephemerides ---
# `times`  : epochs [s], any AbstractRange or Vector{Float64}
# `r_sun`  : Sun position in body-fixed frame [m], Vector of length-3 arrays
#   (typically computed from SPICE kernels — see integration test examples)
ephem = SingleAsteroidEphemerides(times, r_sun)

# --- Thermal parameters ---
k     = 0.1    # Thermal conductivity [W/m/K]
ρ     = 1270.0 # Density [kg/m³]
Cₚ    = 600.0  # Heat capacity [J/kg/K]
R_vis = 0.04   # Reflectance in visible light [-]
R_ir  = 0.0    # Reflectance in thermal infrared [-]
ε     = 1.0    # Emissivity [-]
z_max   = 0.6  # Lower boundary depth [m]
n_depth = 61   # Number of depth nodes
Δz      = z_max / (n_depth - 1)

thermo_params = ThermoParams(k, ρ, Cₚ, R_vis, R_ir, ε, z_max, Δz, n_depth)

# --- Problem definition ---
problem = SingleAsteroidThermoPhysicalProblem(shape, thermo_params;
    with_self_shadowing      = true,
    with_self_heating        = true,
    upper_boundary_condition = RadiationBoundaryCondition(),
    lower_boundary_condition = InsulationBoundaryCondition(),
)

# --- Output specification ---
output_times        = ephem.times[end-119:end]  # final rotation period
subsurface_face_ids = [1, 2, 3]                 # faces for saving subsurface temperature profiles
output = SingleAsteroidOutputSpec(output_times; subsurface_face_ids)

# --- Solve ---
solution = solve(problem, CrankNicolson();
    ephem               = ephem,
    output              = output,
    initial_temperature = 200.0,  # [K]
)

# --- Export results ---
export_solution("output/", solution)
# Writes: diagnostics.csv, surface_temperature.csv, subsurface_temperature.csv

For a complete end-to-end example with SPICE ephemerides, see test/TPM_Ryugu/TPM_Ryugu.jl.

📊 Output

The package produces detailed output files including:

• Surface and subsurface temperature distributions
• Thermal forces and torques
• Energy conservation metrics

🤝 Contributing

Contributions are welcome! Please feel free to submit a Pull Request.