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). .. contents:: Contents :local: :depth: 3 Overview -------- The materials system uses a tiered architecture that balances ease of implementation with GPU performance: .. list-table:: Material Types :header-rows: 1 :widths: 25 35 20 20 * - 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 ^^^^^^^^^^^ .. code-block:: python 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: .. code-block:: python 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: .. code-block:: python 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()``: .. code-block:: python 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: .. code-block:: python 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: .. code-block:: python 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: .. code-block:: python 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: .. code-block:: python 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: .. code-block:: python 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: .. code-block:: python # 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()``: .. code-block:: python 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: .. code-block:: python 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: .. code-block:: python 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: .. code-block:: python 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: .. code-block:: python 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: .. code-block:: python 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: .. code-block:: python 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: .. code-block:: bash # 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: .. code-block:: bash src/lsurf/materials/implementations/my_atmosphere.py 2. Export from ``implementations/__init__.py``: .. code-block:: python from .my_atmosphere import MyExponentialAtmosphere 3. Optionally export from main ``materials/__init__.py``: .. code-block:: python 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: .. code-block:: python """ 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"" ) 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 .. code-block:: python 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: .. code-block:: python 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 -------- - :doc:`../api/materials` - API reference - :doc:`algorithms` - Propagation algorithms - ``scripts/11_atmospheric_refraction.py`` - Example usage - ``scripts/HAHA_refraction.py`` - Duct atmosphere example