API Reference

Insolation

A Julia package to calculate top-of-atmosphere (TOA) insolation (incoming solar radiation) based on Earth's (or another planetary body's) orbital parameters.

The calculations follow fundamental principles of celestial mechanics and solar geometry, as described in "Physics of Earth's Climate" by Tapio Schneider and Lenka Novak.

The package provides functions to:

• Calculate instantaneous insolation for a specific time and location.
• Calculate diurnally averaged insolation.
• Fetch and use orbital parameters (eccentricity, obliquity, and longitude of perihelion) from Laskar et al. (2004) to compute insolation for paleoclimate studies.

InsolationParameters{FT}

The orbital, solar, and epoch parameters needed for insolation calculations.

This is the main concrete subtype of AbstractInsolationParams. It can be constructed directly with the keyword constructor (see the fields below), or, when ClimaParams.jl is loaded, with the convenience constructor InsolationParameters(FT), which fills in default values for modern Earth.

The values stored here describe the planet's orbit at a fixed reference epoch (typically J2000); time-varying (Milankovitch) orbital parameters are supplied separately through Insolation.OrbitalDataSplines.

Fields

• year_anom: Anomalistic year (perihelion to perihelion) [seconds].
• day: Length of a solar day [seconds].
• eccentricity_epoch: Eccentricity at epoch [unitless].
• obliq_epoch: Obliquity (axial tilt) at epoch [radians].
• lon_perihelion_epoch: Longitude of perihelion at epoch [radians].
• tot_solar_irrad: Total solar irradiance at the mean orbital distance (semi-major axis) [W m⁻²].
• epoch: Reference epoch time [DateTime].
• mean_anom_epoch: Mean anomaly at epoch [radians].

orbital_params(od::OrbitalDataSplines, dt::FT) where {FT <: Real}

Interpolate the time-varying orbital parameters (Milankovitch cycles).

Interpolate orbital parameters from the Laskar et al. (2004) dataset for paleoclimate studies. The parameters vary over geological timescales.

Arguments

• od::OrbitalDataSplines: The struct containing the orbital parameter splines.
• dt::FT: The time for interpolation [Julian years since the J2000 epoch].

Returns

• (ϖ, γ, e): A tuple containing:
  • ϖ: Longitude of perihelion [radians]
  • γ: Obliquity (axial tilt) [radians]
  • e: Orbital eccentricity [unitless]

See also the orbital_params(param_set) method for the fixed epoch parameters.

orbital_params(param_set::AIP)

Return the fixed orbital parameters at epoch.

Uses the constant orbital parameter values at the reference epoch (typically J2000) stored in the parameter set. Suitable for modern climate simulations where orbital variations are negligible.

Arguments

• param_set::AIP: Parameter struct containing epoch orbital parameters

Returns

• (ϖ, γ, e): A tuple containing:
  • ϖ: Longitude of perihelion at epoch [radians]
  • γ: Obliquity (axial tilt) at epoch [radians]
  • e: Orbital eccentricity at epoch [unitless]

OrbitalDataSplines

A container struct that holds cubic spline interpolators for Earth's orbital parameters, based on the Laskar et al. (2004) dataset.

The orbital parameter time series is loaded from the laskar2004 artifact, which is downloaded lazily on first use. The splines are functions of time in Julian years since the J2000 epoch (see julian_years_since_epoch).

Fields

• e_spline: Spline for eccentricity [unitless].
• γ_spline: Spline for obliquity [radians].
• ϖ_spline: Spline for longitude of perihelion [radians].

GPU Support

This struct is GPU-compatible via Adapt.jl. To transfer to GPU memory:

using CUDA, Adapt
cpu_od = OrbitalDataSplines()  # Create on CPU
gpu_od = adapt(CuArray, cpu_od)  # Transfer to GPU

TSIDataSpline

A container struct that holds a cubic interpolator for the monthly mean total solar irradiance, from the CMIP record of the Sun's irradiance at 1 au.

The spline is a function of date between 1850-01-15T12:00 and 2229-12-15T12:00. Dates outside this range will result in NaN.

Fields

• tsi_spline: Spline of total solar irradiance in months since ref_date [W m⁻²].
• dates: The monthly dates of the underlying data points.
• ref_date: The reference date corresponding to the first data point.

GPU Support

This struct is GPU-compatible via Adapt.jl. To transfer to GPU memory:

using CUDA, Adapt
cpu_tsi = TSIDataSpline(Float32)  # Create on CPU
gpu_tsi = adapt(CuArray, cpu_tsi)  # Transfer to GPU

TSIDataSpline(::Type{FT}) where {FT <: AbstractFloat}

Construct a TSIDataSpline that interpolates monthly mean total solar irradiance as a function of the date and time.

evaluate(tsi::TSIDataSpline, date::Dates.DateTime)

Interpolate the monthly mean total solar irradiance at date (the Sun's irradiance at 1 au) via a cubic interpolation.

The spline is a function of date between 1850-01-15T12:00 and 2229-12-15T12:00. Dates outside this range will result in NaN.

insolation(
    θ::FT,
    d::FT,
    param_set::IP.AIP,
    date::Union{DateTime,Nothing} = nothing,
    solar_variability_spline::Union{TSIDataSpline,Nothing} = nothing,
) where {FT <: Real}

Calculate top-of-atmosphere (TOA) insolation and the cosine of the solar zenith angle.

Implements $F = S \cos(\theta)$ where S is the solar flux at the given planet-star distance. Insolation is set to 0 at night (when $\cos(\theta) < 0$).

Arguments

• θ::FT: Solar zenith angle [radians]
• d::FT: Planet-star distance, in units of the semi-major axis [unitless]
• param_set::IP.AIP: Parameter struct
• date::Union{DateTime,Nothing}: (default nothing) Current date and time used to evaluate the solar flux if solar_variability_spline is available.
• solar_variability_spline::Union{TSIDataSpline,Nothing}: (default nothing) Use time-varying solar luminosity if a TSIDataSpline is passed as an argument.

Returns

A NamedTuple with fields:

• F: TOA insolation [W m⁻²]
• S: Solar flux at the given planet-star distance [W m⁻²]
• μ: Cosine of solar zenith angle [unitless], clamped to [0, 1]

insolation(
    date::DateTime,
    latitude::FT1,
    longitude::FT2,
    param_set::IP.AIP,
    orbital_data::Union{OrbitalDataSplines, Nothing} = nothing,
    milankovitch::Bool = false,
    solar_variability_spline::Union{TSIDataSpline, Nothing} = nothing,
    eot_correction::Bool = true,
) where {FT1 <: Real, FT2 <: Real}

Calculate instantaneous TOA insolation with optional long-term variations in Earth's orbital parameters (Milankovitch cycles) and solar luminosity.

Arguments

• date::DateTime: Current date and time
• latitude::FT1: Latitude [degrees]
• longitude::FT2: Longitude [degrees]
• param_set::IP.AIP: Parameter struct
• orbital_data::Union{OrbitalDataSplines, Nothing}: (default nothing) Orbital parameter splines. Required when milankovitch=true for GPU compatibility.
• milankovitch::Bool: (default false) Use time-varying orbital parameters (Milankovitch cycles). The OrbitalDataSplines are Earth-specific (Laskar et al., 2004), so this is meaningful for Earth only.
• solar_variability_spline::Union{TSIDataSpline, Nothing}: (default nothing) Use time-varying solar luminosity if TSIDataSpline is passed as an argument. The TSIDataSpline gives the Sun's irradiance at 1 au, so this is meaningful for Earth only.
• eot_correction::Bool: (default true) Apply equation of time correction

Returns

A NamedTuple with fields:

• F: TOA insolation [W m⁻²]
• S: Solar flux [W m⁻²]
• μ: Cosine of solar zenith angle [unitless]
• ζ: Solar azimuth angle [radians], 0 = due East, increasing counterclockwise

Examples

# Modern climate (fixed epoch parameters)
(; F, S, μ, ζ) = insolation(date, lat, lon, param_set)

# Paleoclimate with Milankovitch cycles
od = OrbitalDataSplines()  # Load once
milankovitch = true
(; F, S, μ, ζ) = insolation(date, lat, lon, param_set, od, milankovitch)

# Without equation of time correction
orbital_data = nothing
milankovitch = false
solar_variability_spline = nothing
eot_correction = false
result = insolation(date, lat, lon, param_set, orbital_data, milankovitch, solar_variability_spline, eot_correction)

See also daily_insolation and solar_geometry.

GPU Usage

For GPU execution, create orbital and solar variability data on CPU and transfer to GPU using Adapt.jl:

using CUDA, Adapt
cpu_od = OrbitalDataSplines()  # Create on CPU
gpu_od = adapt(CuArray, cpu_od)  # Transfer to GPU
cpu_solar = TSIDataSpline(Float32) # Create on CPU
gpu_solar = adapt(CuArray, cpu_solar) # Transfer to GPU
# In GPU kernel:
milankovitch=true
result = insolation(date, lat, lon, param_set, gpu_od, milankovitch, gpu_solar)

daily_insolation(
    date::DateTime,
    latitude::Real,
    param_set::IP.AIP,
    orbital_data::Union{OrbitalDataSplines, Nothing} = nothing,
    milankovitch::Bool = false,
    solar_variability_spline::Union{TSIDataSpline, Nothing} = nothing,
)

Calculate diurnally averaged TOA insolation with optional long-term variations in orbital parameters (Milankovitch cycles) and solar luminosity. The insolation is averaged over a full day.

Arguments

• date::DateTime: Current date
• latitude::Real: Latitude [degrees]
• param_set::IP.AIP: Parameter struct
• orbital_data::Union{OrbitalDataSplines, Nothing}: (default nothing) Orbital parameter splines. Required when milankovitch=true for GPU compatibility.
• milankovitch::Bool: (default false) Use time-varying orbital parameters (Milankovitch cycles). The OrbitalDataSplines are Earth-specific (Laskar et al., 2004), so this is meaningful for Earth only.
• solar_variability_spline::Union{TSIDataSpline, Nothing}: (default nothing) Use time-varying solar luminosity if TSIDataSpline is passed as an argument. The TSIDataSpline gives the Sun's irradiance at 1 au, so this is meaningful for Earth only.

Returns

A NamedTuple with fields:

• F: Daily averaged TOA insolation [W m⁻²]
• S: Solar flux [W m⁻²]
• μ: Daily averaged cosine of solar zenith angle [unitless]

Examples

# Modern climate (fixed epoch parameters)
result = daily_insolation(date, lat, param_set)
# Access fields: result.F, result.S, result.μ

# Paleoclimate with Milankovitch cycles
od = OrbitalDataSplines()  # Load once
milankovitch = true
result = daily_insolation(date, lat, param_set, od, milankovitch)

See also insolation.

GPU Usage

For GPU execution, create orbital and solar variability data on CPU and transfer to GPU using Adapt.jl:

using CUDA, Adapt
cpu_od = OrbitalDataSplines()  # Create on CPU
gpu_od = adapt(CuArray, cpu_od)  # Transfer to GPU
cpu_solar = TSIDataSpline(Float32) # Create on CPU
gpu_solar = adapt(CuArray, cpu_solar)
milankovitch = true
# In GPU kernel:
result = daily_insolation(date, lat, param_set, gpu_od, milankovitch, gpu_solar)

solar_flux(
    d::FT,
    param_set::IP.AIP,
    date::Union{DateTime,Nothing} = nothing,
    solar_variability_spline::Union{TSIDataSpline,Nothing} = nothing,
) where {FT<:Real}

Calculate the solar radiative energy flux at the top of the atmosphere (TOA) based on the planet-star distance and the inverse square law.

The total solar irradiance S0 (either tot_solar_irrad or, if a solar_variability_spline is supplied, the interpolated value) is defined at the mean orbital distance (the semi-major axis). Since d is the planet-star distance in units of the semi-major axis, the flux is simply $S = S_0 / d^2$. The absolute size of the orbit does not enter: the flux depends only on S0 and the orbital shape (eccentricity and true anomaly).

Arguments

• d::FT: Planet-star distance, in units of the semi-major axis [unitless]
• param_set::IP.AIP: Struct containing tot_solar_irrad, the total solar irradiance at the mean orbital distance [W m⁻²]
• date::Union{DateTime,Nothing}: (default nothing) Current date and time used to evaluate the solar flux if solar_variability_spline is available.
• solar_variability_spline::Union{TSIDataSpline,Nothing}: (default nothing) Use time-varying solar luminosity if a TSIDataSpline is passed as an argument.

Returns

• S: Solar flux at the given planet-star distance [W m⁻²]

solar_geometry(
    date::DateTime,
    latitude::Real,
    longitude::Real,
    orb_params::Tuple{<:Real, <:Real, <:Real},
    param_set::AIP;
    eot_correction::Bool = true,
)

Calculate the planet-star distance, solar zenith angle, and azimuth angle.

This is a high-level function that combines all necessary astronomical calculations to determine the position of the star in the sky (zenith and azimuth angles) and the planet-star distance at a specific time and location.

All real-valued inputs are converted to eltype(param_set) internally, so the computation is carried out consistently in the parameter set's floating-point type.

Arguments

• date::DateTime: Current date and time
• latitude::Real: Latitude [degrees]
• longitude::Real: Longitude [degrees]
• orb_params::Tuple{<:Real, <:Real, <:Real}: Orbital parameters tuple (ϖ, γ, e):
  • ϖ: Longitude of perihelion [radians]
  • γ: Obliquity (axial tilt) [radians]
  • e: Orbital eccentricity [unitless]
• param_set::AIP: Parameter struct
• eot_correction::Bool: (default true) Apply equation of time correction

Returns

A NamedTuple with fields:

• d: Planet-Sun distance, in units of the semi-major axis [unitless]
• θ: Solar zenith angle [radians]
• ζ: Solar azimuth angle [radians], 0 = due East, increasing CCW

daily_distance_zenith_angle(
    date::DateTime,
    latitude::Real,
    orb_params::Tuple{<:Real, <:Real, <:Real},
    param_set::IP.AIP,
)

Calculate the effective zenith angle for diurnally averaged insolation and the planet-star distance.

Return the effective zenith angle corresponding to the diurnally averaged insolation (averaging the cosine of the zenith angle over 24 hours) and the planet-star distance for a given date and latitude.

All real-valued inputs are converted to eltype(param_set) internally, so the computation is carried out consistently in the parameter set's floating-point type.

Arguments

• date::DateTime: Current date
• latitude::Real: Latitude [degrees]
• orb_params::Tuple{<:Real, <:Real, <:Real}: Orbital parameters tuple (ϖ, γ, e):
  • ϖ: Longitude of perihelion [radians]
  • γ: Obliquity (axial tilt) [radians]
  • e: Orbital eccentricity [unitless]
• param_set::IP.AIP: Parameter struct

Returns

A NamedTuple with fields:

• daily_θ: Effective solar zenith angle [radians]
• d: Planet-star distance, in units of the semi-major axis [unitless]

mean_anomaly(Δt_years::FT, param_set::IP.AIP) where {FT}

Calculate the mean anomaly at a given time since epoch.

The mean anomaly is the angle the planet would have traveled from perihelion if it moved in a circular orbit at constant angular velocity.

Arguments

• Δt_years::FT: Time since epoch [anomalistic years]
• param_set::IP.AIP: Parameter struct containing mean_anom_epoch

Returns

• MA: Mean anomaly [radians]

true_anomaly(MA::FT, e::FT) where {FT <: Real}

Calculate the true anomaly from the mean anomaly.

The true anomaly is the actual angular distance of the planet from perihelion along its orbital path. This function uses a series expansion (the "equation of the center") that is accurate to third order in the eccentricity e, with an error of O(e⁴) (see Fitzpatrick (2012), Appendix A.10).

Arguments

• MA::FT: Mean anomaly [radians]
• e::FT: Orbital eccentricity [unitless]

Returns

• TA: True anomaly [radians]

solar_longitude(TA::FT, ϖ::FT) where {FT <: Real}

Calculate the solar longitude (ecliptic longitude of the Sun).

The solar longitude is the angular distance of the planet along its orbital path, measured from the vernal equinox. It is the sum of the true anomaly (angle from perihelion) and the longitude of perihelion.

Arguments

• TA::FT: True anomaly [radians]
• ϖ::FT: Longitude of perihelion [radians]

Returns

• SL: Solar longitude [radians]

hour_angle(
    date::DateTime,
    λ::FT,
    MA::FT,
    (ϖ, γ, e)::Tuple{FT, FT, FT};
    eot_correction::Bool = true,
) where {FT}

Calculate the hour angle at a given time and longitude.

The hour angle is zero at local solar noon and increases with time. Optionally applies the equation of time correction to account for the difference between apparent and mean solar time.

Arguments

• date::DateTime: Current date and time
• λ::FT: Longitude [radians]
• MA::FT: Mean anomaly [radians]
• (ϖ, γ, e)::Tuple{FT, FT, FT}: Orbital parameters tuple containing:
  • ϖ: Longitude of perihelion [radians]
  • γ: Obliquity (axial tilt) [radians]
  • e: Orbital eccentricity [unitless]
• eot_correction::Bool: (default true) Apply equation of time correction

Returns

• η: Hour angle [radians]

equation_of_time(MA::FT, (ϖ, γ, e)::Tuple{FT, FT, FT}) where {FT <: Real}

Calculate the equation of time correction for the hour angle.

The equation of time accounts for the difference between apparent solar time (based on the actual Sun's position in the sky) and mean solar time (based on a fictitious mean Sun moving at constant speed). This difference arises from two effects:

1. The elliptical orbit (eccentricity e)
2. The axial tilt (obliquity γ)

It is computed as the difference between the mean longitude $L = M + \varpi$ and the right ascension $\alpha$ of the true Sun. The right ascension is obtained from the exact projection of the ecliptic longitude onto the equatorial plane, $\alpha = \mathrm{atan2}(\cos\gamma \sin L_s, \cos L_s)$ with $L_s$ the solar longitude, so the obliquity contribution is exact for any $\gamma$ (it does not rely on a small-tilt expansion). The eccentricity contribution enters through the true anomaly and is therefore accurate to the same order in e as true_anomaly.

Arguments

• MA::FT: Mean anomaly [radians]
• (ϖ, γ, e)::Tuple{FT, FT, FT}: Orbital parameters tuple containing:
  • ϖ: Longitude of perihelion [radians]
  • γ: Obliquity (axial tilt) [radians]
  • e: Orbital eccentricity [unitless]

Returns

• Δη: Hour angle correction [radians]

planet_star_distance(TA::FT, e::FT) where {FT <: Real}

Calculate the planet-star distance at a given true anomaly, in units of the semi-major axis (the mean orbital distance).

The distance varies due to the planet's elliptical orbit, being shortest at perihelion and longest at aphelion. The calculation uses the orbit equation for an ellipse. The result is normalized by the semi-major axis $a$, so it is the dimensionless ratio $d/a = (1-e^2)/(1+e\cos A)$, which ranges from 1-e (perihelion) to 1+e (aphelion). The absolute orbit size is irrelevant to the insolation, so it is not carried; multiply by the semi-major axis to recover a physical distance.

Arguments

• TA::FT: True anomaly [radians]
• e::FT: Orbital eccentricity [unitless]

Returns

• d: Planet-star distance, in units of the semi-major axis [unitless]

years_since_epoch(
    param_set::IP.InsolationParameters{FT},
    date::DateTime,
) where {FT}

Calculate the time elapsed since epoch (typically J2000) in anomalistic years (the time from perihelion to perihelion).

Converts the time difference between two dates from Julian days to anomalistic years, which is the natural time unit for orbital calculations.

Arguments

• param_set::IP.InsolationParameters{FT}: Parameter struct
• date::DateTime: Current date and time

Returns

• Δt_years: Time since epoch [anomalistic years]

See also julian_years_since_epoch, which uses Julian years (the time unit of the Laskar et al. (2004) orbital tables).

julian_years_since_epoch(
    param_set::IP.InsolationParameters{FT},
    date::DateTime,
) where {FT}

Calculate the time elapsed since epoch (typically J2000) in Julian years (a Julian year is exactly 365.25 days of 86400 s each).

This is the time unit of the Laskar et al. (2004) orbital-parameter tables, so it is used to index OrbitalDataSplines. It differs slightly from years_since_epoch, which returns anomalistic years (the natural unit for the mean anomaly); the two differ by the ratio of the anomalistic year to the Julian year (~0.0026% for Earth).

Arguments

• param_set::IP.InsolationParameters{FT}: Parameter struct
• date::DateTime: Current date and time

Returns

• Δt_years: Time since epoch [Julian years]

distance_declination_mean_anomaly(
    Δt_years::FT,
    (ϖ, γ, e)::Tuple{FT, FT, FT},
    param_set::IP.AIP,
) where {FT}

Compute the planet-star distance, solar declination angle, and mean anomaly.

This function calculates key astronomical parameters from orbital elements. The declination determines the subsolar latitude, while the planet-star distance varies due to orbital eccentricity. The mean anomaly is returned for use in hour angle calculations.

Arguments

• Δt_years::FT: Time since epoch [anomalistic years]
• (ϖ, γ, e)::Tuple{FT, FT, FT}: Orbital parameters tuple containing:
  • ϖ: Longitude of perihelion [radians]
  • γ: Obliquity (axial tilt) [radians]
  • e: Orbital eccentricity [unitless]
• param_set::IP.AIP: Parameter struct

Returns

• d: Planet-star distance, in units of the semi-major axis [unitless]
• δ: Solar declination angle [radians]
• MA: Mean anomaly [radians]

AbstractInsolationParams

Abstract base type for insolation parameter sets.

All parameter structs used in Insolation.jl should inherit from this type. The main concrete implementation is InsolationParameters.

This type hierarchy allows for flexible parameter management and enables package extensions to provide alternative parameter implementations while maintaining API compatibility.