Implementing Custom Materials

This guide explains how to create custom materials for the L-SURF ray tracing system. The materials module supports both homogeneous materials (constant properties) and inhomogeneous materials (spatially-varying properties with GPU acceleration).

Overview

The materials system uses a tiered architecture that balances ease of implementation with GPU performance:

Material Types

Type	Use Case	GPU Support	Implementation Effort
HomogeneousMaterial	Constant n (glass, water, air)	N/A (trivial)	None (use directly)
SimpleInhomogeneousModel	Radially symmetric (atmospheres)	Auto via LUT	Implement 1 method
GridInhomogeneousModel	3D fields (turbulence, thermal)	Auto via interpolation	Provide 3D array
FullInhomogeneousModel	Arbitrary Python	CPU only	Implement 1-2 methods

Choose the simplest tier that meets your needs. Higher tiers provide GPU acceleration with minimal implementation effort.

Module Structure

The materials module is organized as follows:

lsurf/materials/
    __init__.py              # Public API
    base/                    # Base classes
        material_field.py    # MaterialField ABC
        homogeneous.py       # HomogeneousMaterial
        simple_inhomogeneous.py   # SimpleInhomogeneousModel
        grid_inhomogeneous.py     # GridInhomogeneousModel
        full_inhomogeneous.py     # FullInhomogeneousModel
    implementations/         # Concrete implementations
        standard_materials.py     # VACUUM, AIR_STP, WATER, BK7_GLASS
        exponential_atmosphere.py # ExponentialAtmosphere
        ...
    utils/                   # Utilities
        constants.py         # Physical constants
        dispersion.py        # Dispersion equations
        factories.py         # Factory functions
    gpu/                     # CUDA kernels
        kernels_simple.py    # SimpleInhomogeneousModel kernels
        kernels_grid.py      # GridInhomogeneousModel kernels

Homogeneous Materials

For materials with constant optical properties, use HomogeneousMaterial directly. No subclassing required.

Basic Usage

from lsurf.materials import HomogeneousMaterial

# Simple material with constant refractive index
glass = HomogeneousMaterial(
    name="Crown Glass",
    refractive_index=1.52,
)

# Material with absorption
colored_glass = HomogeneousMaterial(
    name="Colored Glass",
    refractive_index=1.52,
    absorption_coef=0.1,  # m^-1
)

# Material with scattering
turbid_water = HomogeneousMaterial(
    name="Turbid Water",
    refractive_index=1.333,
    absorption_coef=0.01,
    scattering_coef=0.05,
    anisotropy=0.9,  # Forward scattering
)

Wavelength-Dependent Materials

For dispersion, provide a callable for the refractive index:

from lsurf.materials import HomogeneousMaterial

def water_dispersion(wavelength):
    """Cauchy-like dispersion for water."""
    wl_um = wavelength * 1e6  # Convert to micrometers
    return 1.324 + 0.00306 / (wl_um ** 2)

water = HomogeneousMaterial(
    name="Water (dispersive)",
    refractive_index=water_dispersion,
)

# Test at different wavelengths
n_blue = water.get_refractive_index(0, 0, 0, 450e-9)   # ~1.339
n_red = water.get_refractive_index(0, 0, 0, 650e-9)    # ~1.331

Using Factory Functions

For common dispersion models, use the provided factories:

from lsurf.materials import create_sellmeier_material, create_cauchy_material

# Sellmeier dispersion (accurate for glasses)
bk7 = create_sellmeier_material(
    name="N-BK7",
    B1=1.03961212, C1=0.00600069867,
    B2=0.231792344, C2=0.0200179144,
    B3=1.01046945, C3=103.560653,
)

# Cauchy dispersion (simpler, less accurate)
simple_glass = create_cauchy_material(
    name="Simple Glass",
    A=1.5,      # Base index
    B=0.004,    # First dispersion term
    C=0.0,      # Second dispersion term (usually 0)
)

SimpleInhomogeneousModel

Use SimpleInhomogeneousModel for materials where the refractive index depends only on the distance from a center point (altitude above a reference surface). This is ideal for atmospheric models.

GPU Support: Automatic via precomputed lookup tables (LUT).

When to Use

• Planetary atmospheres (exponential, US Standard, etc.)

• Radially symmetric lenses (GRIN lenses)

• Any spherically symmetric medium

Implementation

Inherit from SimpleInhomogeneousModel and implement n_at_altitude():

import math
from lsurf.materials.base import SimpleInhomogeneousModel

class MyExponentialAtmosphere(SimpleInhomogeneousModel):
    """Custom exponential atmosphere."""

    def __init__(
        self,
        n_surface: float = 1.000293,
        scale_height: float = 8500.0,
        **kwargs
    ):
        super().__init__(name="My Atmosphere", **kwargs)
        self.n_surface = n_surface
        self.scale_height = scale_height
        self.delta_n = n_surface - 1.0

    def n_at_altitude(self, altitude: float, wavelength: float = None) -> float:
        """
        Return refractive index at altitude.

        This is the ONLY method you must implement.

        Parameters
        ----------
        altitude : float
            Height above reference surface in meters (>= 0)
        wavelength : float, optional
            Wavelength in meters (for dispersion)

        Returns
        -------
        float
            Refractive index n >= 1.0
        """
        h = max(altitude, 0.0)
        return 1.0 + self.delta_n * math.exp(-h / self.scale_height)

That’s it! The base class automatically provides:

• get_refractive_index(x, y, z, wavelength) - converts position to altitude

• get_refractive_index_gradient(x, y, z, wavelength) - via radial direction

• GPU kernels via lookup table interpolation

Optional: Analytical Gradient

For better accuracy, override dn_dh_at_altitude() with an analytical derivative:

class MyExponentialAtmosphere(SimpleInhomogeneousModel):
    # ... __init__ and n_at_altitude as above ...

    def dn_dh_at_altitude(self, altitude: float, wavelength: float = None) -> float:
        """
        Return dn/d(altitude) analytically.

        Default uses numerical differentiation. Override for accuracy.
        """
        if altitude < 0:
            return 0.0
        return -(self.delta_n / self.scale_height) * math.exp(
            -altitude / self.scale_height
        )

Configuration Options

The base class accepts several configuration parameters:

atm = MyAtmosphere(
    name="Custom Atmosphere",
    center=(0.0, 0.0, 0.0),         # Center of spherical symmetry
    reference_radius=6_371_000.0,    # Radius where altitude = 0
    altitude_range=(0.0, 200_000.0), # Range for LUT
    lut_resolution=10000,            # LUT samples
)

Complete Example

Here’s a complete implementation of the US Standard Atmosphere 1976:

import math
from lsurf.materials.base import SimpleInhomogeneousModel

class USStandardAtmosphere(SimpleInhomogeneousModel):
    """
    US Standard Atmosphere 1976.

    Piecewise model with temperature lapse rates in different layers.
    """

    # Layer boundaries (altitude in meters)
    LAYERS = [
        (0, 11000, 288.15, -0.0065),      # Troposphere
        (11000, 20000, 216.65, 0.0),      # Tropopause
        (20000, 32000, 216.65, 0.001),    # Stratosphere 1
        (32000, 47000, 228.65, 0.0028),   # Stratosphere 2
    ]

    def __init__(self, **kwargs):
        super().__init__(
            name="US Standard Atmosphere 1976",
            altitude_range=(0.0, 100_000.0),
            **kwargs
        )

    def n_at_altitude(self, altitude: float, wavelength: float = None) -> float:
        """Compute n from pressure and temperature."""
        h = max(altitude, 0.0)

        # Get temperature and pressure at altitude
        T, P = self._get_T_P(h)

        # Refractivity from Edlen formula (simplified)
        N = 77.6 * (P / T)  # Parts per million
        return 1.0 + N * 1e-6

    def _get_T_P(self, h):
        """Get temperature (K) and pressure (Pa) at altitude."""
        # Implementation of piecewise atmosphere model
        # ... (full implementation in implementations/us_standard_atmosphere.py)
        pass

GridInhomogeneousModel

Use GridInhomogeneousModel for 3D inhomogeneous materials defined on a regular grid. The system automatically handles trilinear interpolation and gradient computation.

GPU Support: Automatic via 3D texture interpolation.

When to Use

• Atmospheric turbulence

• Thermal lensing / convection cells

• Pre-computed numerical simulations

• Weather model data

Implementation

For GridInhomogeneousModel, you provide a 3D array rather than implementing a method:

import numpy as np
from lsurf.materials.base import GridInhomogeneousModel

# Create a 3D refractive index field
nx, ny, nz = 64, 64, 32
n_grid = np.ones((nx, ny, nz), dtype=np.float32) * 1.000293

# Add some perturbation (e.g., turbulence)
perturbation = generate_turbulence(nx, ny, nz, strength=1e-6)
n_grid += perturbation

# Create the material
material = GridInhomogeneousModel(
    name="Turbulent Atmosphere",
    n_grid=n_grid,
    bounds=(
        (-5000.0, 5000.0),   # x range in meters
        (-5000.0, 5000.0),   # y range
        (0.0, 20000.0),      # z range (altitude)
    ),
)

Subclassing for Reusability

For reusable grid-based materials, create a subclass:

import numpy as np
from lsurf.materials.base import GridInhomogeneousModel

class GaussianLens(GridInhomogeneousModel):
    """
    Gaussian thermal lens / refractive index bubble.

    Creates a smooth 3D Gaussian perturbation in refractive index.
    """

    def __init__(
        self,
        bounds: tuple,
        grid_resolution: tuple = (64, 64, 32),
        center: tuple = (0.0, 0.0, 5000.0),
        sigma: tuple = (500.0, 500.0, 1000.0),
        delta_n: float = -0.00005,
        background_n: float = 1.000293,
        name: str = "Gaussian Lens",
    ):
        # Generate the grid
        n_grid = self._generate_gaussian_grid(
            bounds, grid_resolution, center, sigma, delta_n, background_n
        )

        super().__init__(
            name=name,
            n_grid=n_grid,
            bounds=bounds,
        )

    def _generate_gaussian_grid(self, bounds, resolution, center, sigma, delta_n, bg_n):
        """Generate 3D Gaussian refractive index field."""
        nx, ny, nz = resolution
        (x_min, x_max), (y_min, y_max), (z_min, z_max) = bounds

        x = np.linspace(x_min, x_max, nx)
        y = np.linspace(y_min, y_max, ny)
        z = np.linspace(z_min, z_max, nz)
        X, Y, Z = np.meshgrid(x, y, z, indexing='ij')

        cx, cy, cz = center
        sx, sy, sz = sigma

        gaussian = np.exp(
            -((X - cx)**2 / (2 * sx**2) +
              (Y - cy)**2 / (2 * sy**2) +
              (Z - cz)**2 / (2 * sz**2))
        )

        n_grid = bg_n + delta_n * gaussian
        return n_grid.astype(np.float32)

Pre-computed Gradients

For better performance, you can provide pre-computed gradients:

# Compute gradients externally
grad_x = np.gradient(n_grid, dx, axis=0)
grad_y = np.gradient(n_grid, dy, axis=1)
grad_z = np.gradient(n_grid, dz, axis=2)
gradient_grid = np.stack([grad_x, grad_y, grad_z], axis=-1)

material = GridInhomogeneousModel(
    name="With Precomputed Gradients",
    n_grid=n_grid,
    bounds=bounds,
    gradient_grid=gradient_grid,  # Shape: (nx, ny, nz, 3)
)

FullInhomogeneousModel

Use FullInhomogeneousModel when you need maximum flexibility and cannot express your model using the other tiers. Any Python code is allowed.

GPU Support: Not available. CPU propagation only.

When to Use

• Complex analytical models

• External data sources (databases, files, APIs)

• Models requiring Python libraries not available on GPU

• Prototyping before optimizing to a higher tier

Implementation

Inherit from FullInhomogeneousModel and implement get_refractive_index():

import numpy as np
from lsurf.materials.base import FullInhomogeneousModel

class LayeredAtmosphere(FullInhomogeneousModel):
    """
    Discrete horizontal layers with different refractive indices.
    """

    def __init__(
        self,
        layer_boundaries: list[float],  # [0, 10000, 20000, ...]
        layer_n_values: list[float],    # [1.0003, 1.0002, 1.0001]
        default_n: float = 1.0,
    ):
        super().__init__(name="Layered Atmosphere")
        self.boundaries = np.array(layer_boundaries)
        self.n_values = np.array(layer_n_values)
        self.default_n = default_n

    def get_refractive_index(
        self, x, y, z, wavelength
    ) -> float | np.ndarray:
        """
        Get refractive index at position.

        This is the ONLY method you must implement.
        Any Python code is allowed here.
        """
        z = np.asarray(z)

        # Find which layer each point is in
        result = np.full_like(z, self.default_n, dtype=float)

        for i in range(len(self.n_values)):
            lower = self.boundaries[i]
            upper = self.boundaries[i + 1]
            mask = (z >= lower) & (z < upper)
            result = np.where(mask, self.n_values[i], result)

        return float(result) if result.ndim == 0 else result

Optional: Analytical Gradient

The default uses numerical differentiation. Override for analytical gradients:

class LayeredAtmosphere(FullInhomogeneousModel):
    # ... __init__ and get_refractive_index as above ...

    def get_refractive_index_gradient(self, x, y, z, wavelength):
        """
        Gradient is zero everywhere except at boundaries.

        For discrete layers, the gradient is technically undefined
        at boundaries (discontinuity). We return zero.
        """
        if isinstance(x, np.ndarray):
            zeros = np.zeros_like(x)
            return zeros, zeros, zeros
        return 0.0, 0.0, 0.0

Using Materials with Propagators

Once you’ve created a material, use it with a propagator:

from lsurf.gpu import create_propagator
from lsurf.utilities.ray_data import create_ray_batch

# Create material (any type)
atm = MyExponentialAtmosphere()

# create_propagator auto-selects CPU or GPU based on material
propagator = create_propagator(atm)

# Create rays
rays = create_ray_batch(num_rays=10000)
# ... initialize ray positions and directions ...

# Propagate
propagator.propagate(
    rays,
    material=atm,
    total_distance=100_000.0,  # meters
    step_size=100.0,           # meters
    wavelength=532e-9,
)

For explicit control:

from lsurf.gpu import GPUGradientPropagator, GradientPropagator

# Force GPU (requires SimpleInhomogeneousModel or GridInhomogeneousModel)
gpu_prop = GPUGradientPropagator(method="rk4")

# Force CPU (works with any material)
cpu_prop = GradientPropagator(method="rk4")

Testing Materials

Comprehensive testing ensures your material works correctly with the propagation system.

Unit Tests

Test the basic interface:

import math
import pytest
import numpy as np

class TestMyAtmosphere:
    """Test suite for MyExponentialAtmosphere."""

    def test_n_at_sea_level(self):
        """Test refractive index at sea level."""
        atm = MyExponentialAtmosphere(n_surface=1.000293)
        assert atm.n_at_altitude(0) == pytest.approx(1.000293, rel=1e-6)

    def test_n_at_high_altitude(self):
        """Test that n approaches 1.0 at high altitude."""
        atm = MyExponentialAtmosphere()
        n_100km = atm.n_at_altitude(100_000)
        assert n_100km == pytest.approx(1.0, abs=1e-6)

    def test_n_decreases_with_altitude(self):
        """Test that n decreases monotonically."""
        atm = MyExponentialAtmosphere()
        altitudes = [0, 1000, 5000, 10000, 50000]
        n_values = [atm.n_at_altitude(h) for h in altitudes]
        assert all(n_values[i] > n_values[i+1] for i in range(len(n_values)-1))

    def test_get_refractive_index_at_position(self):
        """Test 3D position interface."""
        from lsurf.materials import EARTH_RADIUS
        atm = MyExponentialAtmosphere(
            center=(0, 0, 0),
            reference_radius=EARTH_RADIUS,
        )

        # At sea level directly above center
        n = atm.get_refractive_index(0, 0, EARTH_RADIUS, 532e-9)
        assert n == pytest.approx(1.000293, rel=1e-4)

    def test_gradient_direction(self):
        """Test gradient points in correct direction."""
        from lsurf.materials import EARTH_RADIUS
        atm = MyExponentialAtmosphere(
            center=(0, 0, 0),
            reference_radius=EARTH_RADIUS,
        )

        # At sea level directly above center
        grad = atm.get_refractive_index_gradient(0, 0, EARTH_RADIUS, 532e-9)

        # Gradient should point downward (toward higher n)
        assert grad[2] < 0  # Negative z-gradient
        assert abs(grad[0]) < 1e-12  # No x-gradient
        assert abs(grad[1]) < 1e-12  # No y-gradient

    def test_array_inputs(self):
        """Test vectorized evaluation."""
        from lsurf.materials import EARTH_RADIUS
        atm = MyExponentialAtmosphere(
            center=(0, 0, 0),
            reference_radius=EARTH_RADIUS,
        )

        z = EARTH_RADIUS + np.array([0, 1000, 10000])
        n = atm.get_refractive_index(
            np.zeros(3), np.zeros(3), z, 532e-9
        )

        assert n.shape == (3,)
        assert n[0] > n[1] > n[2]

GPU Tests

Test GPU compatibility:

import pytest

try:
    from numba import cuda
    CUDA_AVAILABLE = cuda.is_available()
except ImportError:
    CUDA_AVAILABLE = False


class TestMyAtmosphereGPU:
    """GPU-specific tests."""

    def test_gpu_material_id(self):
        """Test GPU material ID is correct."""
        from lsurf.gpu.interfaces import GPUMaterialID
        atm = MyExponentialAtmosphere()
        assert atm.gpu_material_id == GPUMaterialID.SIMPLE_INHOMOGENEOUS

    def test_lut_initialization(self):
        """Test LUT is created correctly."""
        atm = MyExponentialAtmosphere(lut_resolution=1000)

        # LUT should be lazy-initialized
        assert not atm._lut_initialized

        # Trigger initialization
        arrays = atm.get_gpu_arrays()

        assert atm._lut_initialized
        assert "lut_n" in arrays
        assert len(arrays["lut_n"]) == 1000

    def test_gpu_parameters(self):
        """Test GPU parameters are correct."""
        atm = MyExponentialAtmosphere(
            center=(100.0, 200.0, 300.0),
            reference_radius=6_371_000.0,
        )

        params = atm.get_gpu_parameters()

        assert params[0] == 100.0   # center_x
        assert params[1] == 200.0   # center_y
        assert params[2] == 300.0   # center_z
        assert params[3] == 6_371_000.0  # reference_radius

    @pytest.mark.skipif(not CUDA_AVAILABLE, reason="CUDA not available")
    def test_gpu_propagation(self):
        """Test actual GPU propagation."""
        from lsurf.gpu import GPUGradientPropagator
        from lsurf.utilities.ray_data import create_ray_batch
        from lsurf.materials import EARTH_RADIUS

        atm = MyExponentialAtmosphere()
        prop = GPUGradientPropagator(method="rk4")

        rays = create_ray_batch(num_rays=100)
        rays.positions[:] = [0, 0, EARTH_RADIUS + 10000]
        rays.directions[:] = [1, 0, 0]  # Horizontal
        rays.active[:] = True

        prop.propagate(
            rays,
            material=atm,
            total_distance=10000,
            step_size=100,
            wavelength=532e-9,
        )

        # Rays should have moved
        assert rays.positions[0, 0] > 0
        # Rays should have bent downward (toward higher n)
        assert rays.directions[0, 2] < 0

Integration Tests

Test with the full propagation pipeline:

class TestMyAtmosphereIntegration:
    """Integration tests with propagation system."""

    def test_ray_bending(self):
        """Test that rays bend toward higher refractive index."""
        from lsurf.gpu import create_propagator
        from lsurf.utilities.ray_data import create_ray_batch
        from lsurf.materials import EARTH_RADIUS

        atm = MyExponentialAtmosphere()
        prop = create_propagator(atm)

        # Create horizontal ray at 10 km altitude
        rays = create_ray_batch(num_rays=1)
        rays.positions[0] = [0, 0, EARTH_RADIUS + 10000]
        rays.directions[0] = [1, 0, 0]
        rays.active[0] = True

        # Propagate 100 km
        prop.propagate(
            rays,
            material=atm,
            total_distance=100_000,
            step_size=100,
            wavelength=532e-9,
        )

        # Ray should have bent downward
        final_dir = rays.directions[0]
        assert final_dir[2] < 0, "Ray should bend toward Earth"

    def test_optical_path_length(self):
        """Test optical path length accumulation."""
        from lsurf.gpu import create_propagator
        from lsurf.utilities.ray_data import create_ray_batch
        from lsurf.materials import EARTH_RADIUS

        atm = MyExponentialAtmosphere()
        prop = create_propagator(atm)

        rays = create_ray_batch(num_rays=1)
        rays.positions[0] = [0, 0, EARTH_RADIUS + 10000]
        rays.directions[0] = [1, 0, 0]
        rays.active[0] = True
        rays.optical_path_lengths[0] = 0

        prop.propagate(
            rays,
            material=atm,
            total_distance=1000,
            step_size=10,
            wavelength=532e-9,
        )

        # OPL should be approximately n * geometric_path
        n_avg = atm.n_at_altitude(10000)
        expected_opl = n_avg * 1000
        assert rays.optical_path_lengths[0] == pytest.approx(expected_opl, rel=0.01)

Running Tests

Run your tests with pytest:

# Run all material tests
pytest tests/test_materials.py -v

# Run specific test class
pytest tests/test_materials.py::TestMyAtmosphere -v

# Run with GPU tests (if CUDA available)
pytest tests/test_materials.py -v -m "not slow"

# Run with coverage
pytest tests/test_materials.py --cov=lsurf.materials --cov-report=html

Adding to the Package

To add your material to the L-SURF package:

1. Create the implementation file:

   src/lsurf/materials/implementations/my_atmosphere.py
   

2. Export from implementations/__init__.py:

   from .my_atmosphere import MyExponentialAtmosphere
   

3. Optionally export from main materials/__init__.py:

   from .implementations import MyExponentialAtmosphere
   

4. Add tests in tests/test_materials.py or a new test file.

5. Update documentation in docs/api/materials.rst.

Best Practices

1. Choose the right tier: Use the simplest tier that meets your needs. SimpleInhomogeneousModel provides GPU acceleration with minimal effort.

2. Provide analytical gradients: Override dn_dh_at_altitude() or get_refractive_index_gradient() when possible for better accuracy.

3. Validate inputs: Check that refractive index >= 1.0 and other physical constraints.

4. Document parameters: Include docstrings with units and physical meaning.

5. Test edge cases: Test at boundaries, extreme values, and with array inputs.

6. Consider performance: For GridInhomogeneousModel, balance grid resolution against memory usage and interpolation accuracy.

Example: Complete Custom Atmosphere

Here’s a complete, production-ready implementation:

"""
Custom Exponential Atmosphere

A configurable exponential atmosphere model with GPU support.
"""

import math
from lsurf.materials.base import SimpleInhomogeneousModel
from lsurf.materials.utils.constants import EARTH_RADIUS

class CustomExponentialAtmosphere(SimpleInhomogeneousModel):
    """
    Exponential atmosphere with configurable parameters.

    The refractive index follows:
        n(h) = 1 + (n_surface - 1) * exp(-h / H)

    Parameters
    ----------
    n_surface : float
        Refractive index at the reference surface. Default: 1.000293
    scale_height : float
        Atmospheric scale height in meters. Default: 8500.0
    earth_radius : float
        Reference radius in meters. Default: 6,371,000
    earth_center : tuple
        Center of spherical symmetry. Default: (0, 0, 0)

    Examples
    --------
    >>> atm = CustomExponentialAtmosphere()
    >>> n = atm.n_at_altitude(10000)  # At 10 km
    >>> print(f"n = {n:.6f}")
    n = 1.000089

    >>> # Mars atmosphere
    >>> mars = CustomExponentialAtmosphere(
    ...     n_surface=1.0001,
    ...     scale_height=11100,
    ...     earth_radius=3389500,
    ... )
    """

    def __init__(
        self,
        n_surface: float = 1.000293,
        scale_height: float = 8500.0,
        earth_radius: float = EARTH_RADIUS,
        earth_center: tuple[float, float, float] = (0.0, 0.0, 0.0),
        name: str = "Custom Exponential Atmosphere",
    ):
        # Validate inputs
        if n_surface < 1.0:
            raise ValueError(f"n_surface must be >= 1.0, got {n_surface}")
        if scale_height <= 0:
            raise ValueError(f"scale_height must be > 0, got {scale_height}")

        super().__init__(
            name=name,
            center=earth_center,
            reference_radius=earth_radius,
            altitude_range=(0.0, 15 * scale_height),
            lut_resolution=10000,
        )

        self.n_surface = n_surface
        self.scale_height = scale_height
        self.delta_n = n_surface - 1.0

    def n_at_altitude(self, altitude: float, wavelength: float = None) -> float:
        """Return refractive index at altitude."""
        h = max(altitude, 0.0)
        return 1.0 + self.delta_n * math.exp(-h / self.scale_height)

    def dn_dh_at_altitude(self, altitude: float, wavelength: float = None) -> float:
        """Return analytical dn/dh at altitude."""
        if altitude < 0:
            return 0.0
        return -(self.delta_n / self.scale_height) * math.exp(
            -altitude / self.scale_height
        )

    def __repr__(self) -> str:
        return (
            f"<CustomExponentialAtmosphere("
            f"n_surface={self.n_surface:.6f}, "
            f"H={self.scale_height/1000:.1f} km)>"
        )

Atmospheric Duct Example

The DuctAtmosphere class extends the exponential atmosphere with an inversion layer (duct) that can trap or guide electromagnetic waves. This is useful for modeling phenomena like:

• Surface evaporation ducts over oceans

• Elevated ducts causing anomalous radio propagation

• Mirage formation in temperature inversion layers

from lsurf.materials import DuctAtmosphere, EARTH_RADIUS

# Create an atmosphere with a surface duct at 3 km altitude
atmosphere = DuctAtmosphere(
    n_sea_level=1.000293,
    scale_height=8500.0,
    earth_center=(0.0, 0.0, -EARTH_RADIUS),
    duct_center=3000.0,      # Center of duct at 3 km
    duct_width=1000.0,       # 1 km wide duct
    duct_intensity=0.8,      # Strong duct effect
    duct_sharpness=0.9,      # Sharp duct boundaries
)

# The duct modifies the refractive index profile
# creating a localized perturbation that can trap rays
print(f"n at 2.5 km: {atmosphere.n_at_altitude(2500):.6f}")
print(f"n at 3.0 km: {atmosphere.n_at_altitude(3000):.6f}")
print(f"n at 3.5 km: {atmosphere.n_at_altitude(3500):.6f}")

The duct parameters control the shape of the perturbation:

• duct_center: Altitude of the duct center in meters

• duct_width: Vertical extent of the duct in meters

• duct_intensity: Strength of the duct effect (0 = none, 1 = maximum)

• duct_sharpness: How abrupt the duct boundaries are (0 = gradual, 1 = sharp)

Using with the source module:

from lsurf.materials import DuctAtmosphere, EARTH_RADIUS
from lsurf.sources import ParallelBeamFromPositions
from lsurf.gpu import create_propagator
import numpy as np

# Create duct atmosphere
atmosphere = DuctAtmosphere(
    earth_center=(0.0, 0.0, -EARTH_RADIUS),
    duct_center=3000.0,
    duct_width=1000.0,
    duct_intensity=0.8,
    duct_sharpness=0.9,
)

# Create rays at different impact parameters
impact_parameters = np.arange(0, 15000, 500)
prop_dist = np.sqrt((EARTH_RADIUS + 100000)**2 - (EARTH_RADIUS + impact_parameters)**2)

positions = np.column_stack([
    -prop_dist,
    np.zeros_like(impact_parameters),
    impact_parameters.astype(float),
])

source = ParallelBeamFromPositions(
    positions=positions,
    direction=(1.0, 0.0, 0.0),
)
rays = source.generate()

# Propagate through the duct
propagator = create_propagator(atmosphere)
propagator.propagate(
    rays,
    total_distance=500_000,
    step_size=500,
    wavelength=532e-9,
)

See Also

• Materials Module - API reference

• Algorithms - Propagation algorithms

• scripts/11_atmospheric_refraction.py - Example usage

• scripts/HAHA_refraction.py - Duct atmosphere example