# 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.
"""
Atmospheric Duct Material
Implements an exponential atmosphere with a refractive index duct (inversion layer).
The duct creates a localized modification to the refractive index profile that can
trap or guide electromagnetic waves.
The refractive index follows:
n(h) = 1 + delta_n * exp(-h / H) * D(h)
where D(h) is a duct factor based on hyperbolic tangent functions that creates
a localized perturbation at a specified altitude.
This model is appropriate for studying:
- Atmospheric ducting phenomena
- Anomalous radio propagation
- Mirage formation in inversion layers
- Over-the-horizon radar propagation
References
----------
.. [1] Hitney, H.V. et al. (1985). "Tropospheric radio propagation assessment."
Proceedings of the IEEE, 73(2), 265-283.
.. [2] Turton, J.D. et al. (1988). "The structure of evaporation ducts."
Radio Science, 23(4), 519-528.
"""
import math
import numpy as np
from ..base import SimpleInhomogeneousModel
from ..utils.constants import EARTH_RADIUS, SCALE_HEIGHT_DEFAULT, N_SEA_LEVEL
[docs]
class DuctAtmosphere(SimpleInhomogeneousModel):
"""
Exponential atmosphere with a refractive index duct layer.
Models an atmosphere where air density decreases exponentially with
altitude, modified by a duct (inversion layer) at a specified altitude.
The duct creates a localized region where the refractive index gradient
is modified, which can trap or guide electromagnetic waves.
The duct is modeled using hyperbolic tangent functions to create smooth
transitions at the duct boundaries.
Parameters
----------
name : str, optional
Descriptive name for this material. Default is "Duct Atmosphere".
n_sea_level : float, optional
Refractive index at sea level. 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).
duct_center : float, optional
Center altitude of the duct layer in meters. Default is 0.0.
duct_width : float, optional
Width of the duct layer in meters. Default is 100.0.
duct_intensity : float, optional
Intensity of the duct effect (0 = no duct, 1 = full strength).
Default is 0.0.
duct_sharpness : float, optional
Sharpness of duct boundaries (0 = gradual, 1 = sharp).
Default is 0.5.
Examples
--------
>>> # Surface duct at 100m altitude
>>> atmosphere = DuctAtmosphere(
... duct_center=100.0,
... duct_width=50.0,
... duct_intensity=0.5,
... duct_sharpness=0.8,
... )
>>> # Elevated duct for radio propagation study
>>> atmosphere = DuctAtmosphere(
... duct_center=1000.0,
... duct_width=200.0,
... duct_intensity=0.8,
... duct_sharpness=0.9,
... earth_center=(0.0, 0.0, -EARTH_RADIUS),
... )
"""
[docs]
def __init__(
self,
name: str = "Duct 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),
duct_center: float = 0.0,
duct_width: float = 100.0,
duct_intensity: float = 0.0,
duct_sharpness: float = 0.5,
):
# 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}")
if duct_width <= 0:
raise ValueError(f"Duct width must be positive, got {duct_width}")
if not 0.0 <= duct_sharpness <= 1.0:
raise ValueError(f"Duct sharpness must be in [0, 1], got {duct_sharpness}")
# Initialize base class
super().__init__(
name=name,
center=earth_center,
reference_radius=earth_radius,
altitude_range=(0.0, 15 * scale_height),
lut_resolution=10000,
)
# Store atmosphere parameters
self.n_sea_level = n_sea_level
self.scale_height = scale_height
self.earth_radius = earth_radius
self.earth_center = earth_center
# Store duct parameters
self.duct_center = duct_center
self.duct_width = duct_width
self.duct_intensity = duct_intensity
self.duct_sharpness = duct_sharpness
# Precompute refractivity
self.delta_n = n_sea_level - 1.0
def _compute_duct_factor(self, altitude: float) -> float:
"""
Compute the duct modification factor at given altitude.
The duct factor modifies the exponential profile to create a
localized perturbation. Uses hyperbolic tangent functions for
smooth transitions.
Parameters
----------
altitude : float
Altitude above Earth's surface in meters.
Returns
-------
float
Duct factor (1.0 = no modification, <1.0 = reduced refractivity).
"""
if self.duct_intensity == 0.0:
return 1.0
# Compute transition width based on sharpness
alpha_min = self.duct_width / 50
alpha_max = self.duct_width / 2
alpha = alpha_max * (1 - self.duct_sharpness) + alpha_min * self.duct_sharpness
# Normalization factor
norm = 2 * np.tanh(self.duct_width / (2 * alpha))
# Duct profile using difference of tanh functions
lower_edge = (altitude - self.duct_center + self.duct_width / 2) / alpha
upper_edge = (altitude - self.duct_center - self.duct_width / 2) / alpha
duct_factor = 1 - self.duct_intensity / norm * (
np.tanh(lower_edge) - np.tanh(upper_edge)
)
return float(duct_factor)
def _compute_duct_factor_derivative(self, altitude: float) -> float:
"""
Compute the derivative of the duct factor with respect to altitude.
Parameters
----------
altitude : float
Altitude above Earth's surface in meters.
Returns
-------
float
Derivative dD/dh of the duct factor.
"""
if self.duct_intensity == 0.0:
return 0.0
# Compute transition width based on sharpness
alpha_min = self.duct_width / 50
alpha_max = self.duct_width / 2
alpha = alpha_max * (1 - self.duct_sharpness) + alpha_min * self.duct_sharpness
# Normalization factor
norm = 2 * np.tanh(self.duct_width / (2 * alpha))
# Derivative of tanh is sech^2
lower_edge = (altitude - self.duct_center + self.duct_width / 2) / alpha
upper_edge = (altitude - self.duct_center - self.duct_width / 2) / alpha
# sech²(x) → 0 for |x| >> 1; guard against cosh overflow
sech2_lower = 1.0 / np.cosh(lower_edge) ** 2 if abs(lower_edge) < 350 else 0.0
sech2_upper = 1.0 / np.cosh(upper_edge) ** 2 if abs(upper_edge) < 350 else 0.0
d_duct_dh = -(
self.duct_intensity / (norm * alpha) * (sech2_lower - sech2_upper)
)
return float(d_duct_dh)
# =========================================================================
# 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 with duct modification:
n(h) = 1 + delta_n * exp(-h / H) * D(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)
duct_factor = self._compute_duct_factor(altitude_clamped)
return float(
1.0
+ self.delta_n
* math.exp(-altitude_clamped / self.scale_height)
* duct_factor
)
[docs]
def dn_dh_at_altitude(
self, altitude: float, wavelength: float | None = None
) -> float:
"""
Return analytical dn/dh at given altitude.
Uses the product rule for the derivative:
dn/dh = delta_n * [exp' * D + exp * D']
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
# Exponential factor and its derivative
exp_factor = math.exp(-altitude / self.scale_height)
d_exp_dh = -exp_factor / self.scale_height
# Duct factor and its derivative
duct_factor = self._compute_duct_factor(altitude)
d_duct_dh = self._compute_duct_factor_derivative(altitude)
# Product rule: d(exp * D)/dh = exp' * D + exp * D'
return float(self.delta_n * (d_exp_dh * duct_factor + exp_factor * d_duct_dh))
# =========================================================================
# GPU Material Interface (Overrides)
# =========================================================================
@property
def gpu_material_id(self) -> int:
"""Return GPU material ID for analytical duct kernels."""
from ...propagation.propagator_protocol import GPUMaterialID
return GPUMaterialID.DUCT_ATMOSPHERE
[docs]
def get_gpu_kernels(self) -> dict:
"""Return GPU kernels for propagation (analytical, no LUT)."""
from ...propagation.kernels.propagation.duct import (
_kernel_duct_euler,
_kernel_duct_rk4,
_kernel_duct_euler_adaptive,
_kernel_duct_rk4_adaptive,
)
return {
"euler": _kernel_duct_euler,
"rk4": _kernel_duct_rk4,
"euler_adaptive": _kernel_duct_euler_adaptive,
"rk4_adaptive": _kernel_duct_rk4_adaptive,
}
[docs]
def get_gpu_parameters(self) -> tuple:
"""Return scalar parameters for GPU kernel."""
# Precompute alpha and norm from duct_sharpness and duct_width
alpha_min = self.duct_width / 50
alpha_max = self.duct_width / 2
alpha = alpha_max * (1 - self.duct_sharpness) + alpha_min * self.duct_sharpness
norm = 2 * np.tanh(self.duct_width / (2 * alpha))
return (
float(self.earth_center[0]),
float(self.earth_center[1]),
float(self.earth_center[2]),
float(self.earth_radius),
float(self.delta_n),
float(self.scale_height),
float(self.duct_center),
float(self.duct_width),
float(self.duct_intensity),
float(alpha),
float(norm),
)
[docs]
def get_duct_gpu_parameters(self) -> tuple:
"""Return scalar parameters for analytical duct GPU kernels.
Equivalent to get_gpu_parameters() for DuctAtmosphere, but provides
a consistent interface for propagators to detect analytical duct kernels.
"""
return self.get_gpu_parameters()
[docs]
def get_gpu_arrays(self) -> dict:
"""Return device arrays for GPU kernel (none needed for analytical)."""
return {}
[docs]
def __repr__(self) -> str:
"""Return string representation."""
return (
f"<DuctAtmosphere("
f"n_sea_level={self.n_sea_level:.6f}, "
f"H={self.scale_height / 1000:.1f} km, "
f"duct_center={self.duct_center:.0f} m, "
f"duct_width={self.duct_width:.0f} m, "
f"duct_intensity={self.duct_intensity:.2f})>"
)
[docs]
def create_duct_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),
duct_center: float = 0.0,
duct_width: float = 100.0,
duct_intensity: float = 0.0,
duct_sharpness: float = 0.5,
name: str = "Custom Duct Atmosphere",
) -> DuctAtmosphere:
"""
Factory function to create a duct 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).
duct_center : float, optional
Center altitude of the duct in meters. Default is 0.0.
duct_width : float, optional
Width of the duct in meters. Default is 100.0.
duct_intensity : float, optional
Intensity of the duct (0-1). Default is 0.0.
duct_sharpness : float, optional
Sharpness of duct edges (0-1). Default is 0.5.
name : str, optional
Descriptive name for the material.
Returns
-------
DuctAtmosphere
Configured atmosphere material with duct.
Examples
--------
>>> # Create a surface evaporation duct
>>> atmo = create_duct_atmosphere(
... duct_center=50.0,
... duct_width=30.0,
... duct_intensity=0.7,
... duct_sharpness=0.9,
... )
"""
return DuctAtmosphere(
name=name,
n_sea_level=n_sea_level,
scale_height=scale_height,
earth_radius=earth_radius,
earth_center=earth_center,
duct_center=duct_center,
duct_width=duct_width,
duct_intensity=duct_intensity,
duct_sharpness=duct_sharpness,
)
__all__ = [
"DuctAtmosphere",
"create_duct_atmosphere",
]