# The Clear BSD License
#
# Copyright (c) 2026 Tobias Heibges
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted (subject to the limitations in the disclaimer
# below) provided that the following conditions are met:
#
# * Redistributions of source code must retain the above copyright notice,
# this list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in the
# documentation and/or other materials provided with the distribution.
#
# * Neither the name of the copyright holder nor the names of its
# contributors may be used to endorse or promote products derived from this
# software without specific prior written permission.
#
# NO EXPRESS OR IMPLIED LICENSES TO ANY PARTY'S PATENT RIGHTS ARE GRANTED BY
# THIS LICENSE. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
# CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
# IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
"""
Exponential Atmosphere Material
Implements a spherically symmetric exponential atmosphere where the air density
(and thus refractive index) depends on the radial distance from Earth's center.
The refractive index follows:
n(r) = 1 + (n_0 - 1) * exp(-(r - R_E) / H)
where:
r = radial distance from Earth's center
R_E = Earth's radius
n_0 = refractive index at sea level (~1.000293)
H = scale height (~8.5 km for Earth's atmosphere)
This model is appropriate for studying atmospheric refraction effects
such as mirages, stellar aberration, and ray bending in grazing incidence
scenarios.
References
----------
.. [1] Garfinkel, B. (1967). "Astronomical Refraction in a Polytropic
Atmosphere". The Astronomical Journal, 72, 235-254.
.. [2] Smart, W.M. (1977). "Textbook on Spherical Astronomy", 6th ed.
Cambridge University Press, Chapter XI.
"""
from __future__ import annotations
import math
from typing import TYPE_CHECKING, Union
import numpy as np
from numpy.typing import NDArray
from ..base import SimpleInhomogeneousModel
from ..utils.constants import EARTH_RADIUS, SCALE_HEIGHT_DEFAULT, N_SEA_LEVEL
if TYPE_CHECKING:
from ...propagation.kernels.registry import PropagationKernelID, PropagatorID
[docs]
class ExponentialAtmosphere(SimpleInhomogeneousModel):
"""
Exponential atmosphere with radially-dependent refractive index.
Models an atmosphere where air density decreases exponentially with
altitude, causing the refractive index to approach 1 at high altitudes.
The model is spherically symmetric about Earth's center, which is
assumed to be at the origin (0, 0, 0) by default, or can be specified.
This class inherits from SimpleInhomogeneousModel, implementing the
`n_at_altitude()` method for the exponential profile. GPU support is
provided automatically via lookup table interpolation.
Parameters
----------
name : str, optional
Descriptive name for this material. Default is "Exponential Atmosphere".
n_sea_level : float, optional
Refractive index at sea level (Earth's surface). Default is 1.000293.
scale_height : float, optional
Atmospheric scale height H in meters. Default is 8500.0 m.
earth_radius : float, optional
Radius of Earth in meters. Default is 6,371,000 m.
earth_center : tuple of float, optional
Position of Earth's center in meters. Default is (0, 0, 0).
absorption_coef : float, optional
Absorption coefficient at sea level in m⁻¹. Default is 0.0.
absorption_scale_height : float, optional
Scale height for absorption (can differ from density). Default is same as scale_height.
Examples
--------
>>> atmosphere = ExponentialAtmosphere()
>>> # Get refractive index at 10 km altitude
>>> n = atmosphere.get_refractive_index(0, 0, EARTH_RADIUS + 10000, 532e-9)
>>> print(f"n at 10 km: {n:.6f}") # ~1.000089
>>> # Get gradient at sea level directly above Earth's center
>>> grad = atmosphere.get_refractive_index_gradient(0, 0, EARTH_RADIUS, 532e-9)
>>> print(f"dn/dz at sea level: {grad[2]:.2e}") # ~ -3.4e-8 m^-1
"""
[docs]
def __init__(
self,
name: str = "Exponential Atmosphere",
n_sea_level: float = N_SEA_LEVEL,
scale_height: float = SCALE_HEIGHT_DEFAULT,
earth_radius: float = EARTH_RADIUS,
earth_center: tuple[float, float, float] = (0.0, 0.0, 0.0),
absorption_coef: float = 0.0,
absorption_scale_height: float | None = None,
kernel: PropagationKernelID | None = None,
propagator: PropagatorID | None = None,
):
# Validate inputs
if scale_height <= 0:
raise ValueError(f"Scale height must be positive, got {scale_height}")
if n_sea_level < 1.0:
raise ValueError(f"Refractive index must be >= 1.0, got {n_sea_level}")
# Initialize base class
super().__init__(
name=name,
center=earth_center,
reference_radius=earth_radius,
altitude_range=(0.0, 15 * scale_height), # ~127 km for default
lut_resolution=10000,
kernel=kernel,
propagator=propagator,
)
# Store exponential atmosphere specific parameters
self.n_sea_level = n_sea_level
self.scale_height = scale_height
self.earth_radius = earth_radius
self.earth_center = earth_center
self.absorption_coef_sea_level = absorption_coef
self.absorption_scale_height = (
absorption_scale_height
if absorption_scale_height is not None
else scale_height
)
# Precompute refractivity
self.delta_n = n_sea_level - 1.0
# =========================================================================
# SimpleInhomogeneousModel Interface (Required)
# =========================================================================
[docs]
def n_at_altitude(self, altitude: float, wavelength: float | None = None) -> float:
"""
Return refractive index at given altitude.
Implements the exponential profile:
n(h) = 1 + delta_n * exp(-h / H)
Parameters
----------
altitude : float
Altitude above Earth's surface in meters (clamped to >= 0)
wavelength : float, optional
Wavelength in meters (not used - no dispersion)
Returns
-------
float
Refractive index at altitude
"""
altitude_clamped = max(altitude, 0.0)
return 1.0 + self.delta_n * math.exp(-altitude_clamped / self.scale_height)
[docs]
def dn_dh_at_altitude(
self, altitude: float, wavelength: float | None = None
) -> float:
"""
Return analytical dn/dh at given altitude.
Derivative of exponential profile:
dn/dh = -(delta_n / H) * exp(-h / H)
Parameters
----------
altitude : float
Altitude above Earth's surface in meters
wavelength : float, optional
Wavelength in meters (not used)
Returns
-------
float
Derivative dn/dh in m^-1
"""
if altitude < 0:
return 0.0
return -(self.delta_n / self.scale_height) * math.exp(
-altitude / self.scale_height
)
[docs]
def alpha_at_altitude(
self, altitude: float, wavelength: float | None = None
) -> float:
"""
Return absorption coefficient at given altitude.
Implements exponential absorption profile:
α(h) = α₀ * exp(-h / H_α)
Parameters
----------
altitude : float
Altitude above Earth's surface in meters
wavelength : float, optional
Wavelength in meters (not used)
Returns
-------
float
Absorption coefficient α in m⁻¹
"""
if self.absorption_coef_sea_level == 0.0:
return 0.0
altitude_clamped = max(altitude, 0.0)
return self.absorption_coef_sea_level * math.exp(
-altitude_clamped / self.absorption_scale_height
)
# =========================================================================
# Absorption (Override base class default)
# =========================================================================
[docs]
def get_absorption_coefficient(
self,
x: Union[float, NDArray[np.float64]],
y: Union[float, NDArray[np.float64]],
z: Union[float, NDArray[np.float64]],
wavelength: Union[float, NDArray[np.float64]],
) -> Union[float, NDArray[np.float64]]:
"""
Get absorption coefficient at position (x, y, z).
Absorption follows an exponential profile with altitude.
Parameters
----------
x, y, z : float or ndarray
Position coordinates in meters.
wavelength : float or ndarray
Wavelength in meters.
Returns
-------
alpha : float or ndarray
Absorption coefficient in m⁻¹.
"""
if self.absorption_coef_sea_level == 0.0:
if isinstance(x, np.ndarray):
return np.zeros_like(x)
return 0.0
altitude = self._compute_altitude(x, y, z)
if isinstance(altitude, np.ndarray):
alpha = self.absorption_coef_sea_level * np.exp(
-altitude / self.absorption_scale_height
)
else:
alpha = self.absorption_coef_sea_level * math.exp(
-altitude / self.absorption_scale_height
)
return alpha
# =========================================================================
# Atmosphere-Specific Utilities
# =========================================================================
[docs]
def get_density_ratio(
self,
x: Union[float, NDArray[np.float64]],
y: Union[float, NDArray[np.float64]],
z: Union[float, NDArray[np.float64]],
) -> Union[float, NDArray[np.float64]]:
"""
Get atmospheric density ratio relative to sea level.
Parameters
----------
x, y, z : float or ndarray
Position coordinates in meters.
Returns
-------
rho_ratio : float or ndarray
Density relative to sea level (dimensionless).
rho_ratio = rho(r) / rho_0 = exp(-altitude / H)
"""
altitude = self._compute_altitude(x, y, z)
if isinstance(altitude, np.ndarray):
return np.exp(-altitude / self.scale_height)
else:
return math.exp(-altitude / self.scale_height)
[docs]
def compute_curvature_radius(
self,
x: Union[float, NDArray[np.float64]],
y: Union[float, NDArray[np.float64]],
z: Union[float, NDArray[np.float64]],
wavelength: float = 532e-9,
) -> Union[float, NDArray[np.float64]]:
"""
Compute the radius of curvature for a ray at given position.
The radius of curvature is R_c = n / |∇n|.
Parameters
----------
x, y, z : float or ndarray
Position coordinates in meters.
wavelength : float, optional
Wavelength in meters. Default is 532 nm.
Returns
-------
R_c : float or ndarray
Radius of curvature in meters.
Very large values indicate nearly straight propagation.
"""
n = self.get_refractive_index(x, y, z, wavelength)
grad_mag = self.get_refractive_index_gradient_magnitude(x, y, z, wavelength)
if isinstance(grad_mag, np.ndarray):
grad_mag_safe = np.where(grad_mag < 1e-15, 1e-15, grad_mag)
return n / grad_mag_safe
else:
if grad_mag < 1e-15:
return float("inf")
return n / grad_mag
[docs]
def __repr__(self) -> str:
"""Return string representation."""
return (
f"<ExponentialAtmosphere("
f"n_sea_level={self.n_sea_level:.6f}, "
f"H={self.scale_height/1000:.1f} km, "
f"R_E={self.earth_radius/1000:.0f} km)>"
)
# Pre-configured atmosphere instances
STANDARD_ATMOSPHERE = ExponentialAtmosphere(
name="Standard Atmosphere",
n_sea_level=N_SEA_LEVEL,
scale_height=SCALE_HEIGHT_DEFAULT,
)
"""
Standard exponential atmosphere model.
Uses typical values:
- Sea level refractive index: 1.000293
- Scale height: 8.5 km
- Earth radius: 6371 km
- Earth center at origin
"""
[docs]
def create_exponential_atmosphere(
n_sea_level: float = N_SEA_LEVEL,
scale_height: float = SCALE_HEIGHT_DEFAULT,
earth_radius: float = EARTH_RADIUS,
earth_center: tuple[float, float, float] = (0.0, 0.0, 0.0),
name: str = "Custom Exponential Atmosphere",
) -> ExponentialAtmosphere:
"""
Factory function to create an exponential atmosphere.
Parameters
----------
n_sea_level : float, optional
Refractive index at sea level. Default is 1.000293.
scale_height : float, optional
Atmospheric scale height in meters. Default is 8500.0 m.
earth_radius : float, optional
Radius of Earth in meters. Default is 6,371,000 m.
earth_center : tuple of float, optional
Position of Earth's center. Default is (0, 0, 0).
name : str, optional
Descriptive name for the material.
Returns
-------
ExponentialAtmosphere
Configured atmosphere material.
Examples
--------
>>> # Mars-like atmosphere (thinner, smaller planet)
>>> mars_atmo = create_exponential_atmosphere(
... n_sea_level=1.0001,
... scale_height=11_100.0, # Mars scale height
... earth_radius=3_389_500.0, # Mars radius
... )
"""
return ExponentialAtmosphere(
name=name,
n_sea_level=n_sea_level,
scale_height=scale_height,
earth_radius=earth_radius,
earth_center=earth_center,
)
__all__ = [
"ExponentialAtmosphere",
"STANDARD_ATMOSPHERE",
"create_exponential_atmosphere",
]