# 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.
"""
Material Field Base Class
Defines the abstract base class for all material types. Materials provide
spatially-varying optical properties including refractive index, absorption,
and scattering coefficients.
Design Notes
------------
- Materials are position-dependent by default (field concept)
- Homogeneous materials can optimize for constant properties
- GPU-compatible interface for raytracing kernels
- Materials declare supported kernels and propagators for compatibility checking
"""
from __future__ import annotations
import math
from abc import ABC, abstractmethod
from typing import TYPE_CHECKING
import numpy as np
from numpy.typing import NDArray
if TYPE_CHECKING:
from ...propagation.kernels.registry import PropagationKernelID, PropagatorID
[docs]
class MaterialField(ABC):
"""
Abstract base class for spatially-varying material properties.
All material properties are functions of position (x, y, z) and wavelength.
Subclasses must implement methods that can be called from both CPU and GPU.
Parameters
----------
name : str, optional
Descriptive name for this material. Default is "Material".
Attributes
----------
name : str
Material identifier.
Notes
-----
Material properties modeled:
- Refractive index n: Real part of complex refractive index
- Absorption coefficient α: Controls intensity decay
- Scattering coefficient μ_s: Controls scattering mean free path
- Anisotropy factor g: Henyey-Greenstein scattering parameter
References
----------
.. [1] Born, M., & Wolf, E. (1999). Principles of Optics (7th ed.).
Cambridge University Press.
"""
# =========================================================================
# CLASS-LEVEL COMPATIBILITY DECLARATIONS
# =========================================================================
# Subclasses should override these to declare supported kernels/propagators
_supported_kernels: list[PropagationKernelID] = []
_default_kernel: PropagationKernelID | None = None
_supported_propagators: list[PropagatorID] = []
_default_propagator: PropagatorID | None = None
[docs]
def __init__(
self,
name: str = "Material",
kernel: PropagationKernelID | None = None,
propagator: PropagatorID | None = None,
):
"""
Initialize material field.
Parameters
----------
name : str, optional
Descriptive name for this material. Default is "Material".
kernel : PropagationKernelID, optional
Override the default propagation kernel. Must be in supported_kernels().
If None, uses the class default.
propagator : PropagatorID, optional
Override the default propagator. Must be in supported_propagators().
If None, uses the class default.
Raises
------
ValueError
If kernel or propagator is not supported by this material type.
"""
self.name = name
self._is_homogeneous = False
# Resolve kernel preference
if kernel is None:
self._kernel_id = self._default_kernel
else:
if self._supported_kernels and kernel not in self._supported_kernels:
raise ValueError(
f"{self.__class__.__name__} does not support kernel {kernel}. "
f"Supported: {self._supported_kernels}"
)
self._kernel_id = kernel
# Resolve propagator preference
if propagator is None:
self._propagator_id = self._default_propagator
else:
if (
self._supported_propagators
and propagator not in self._supported_propagators
):
raise ValueError(
f"{self.__class__.__name__} does not support propagator {propagator}. "
f"Supported: {self._supported_propagators}"
)
self._propagator_id = propagator
[docs]
@abstractmethod
def get_refractive_index(
self,
x: float | NDArray[np.float64],
y: float | NDArray[np.float64],
z: float | NDArray[np.float64],
wavelength: float | NDArray[np.float64],
) -> float | NDArray[np.float64]:
"""
Get refractive index at position (x, y, z) for given wavelength.
Parameters
----------
x, y, z : float or ndarray
Position coordinates in meters.
wavelength : float or ndarray
Wavelength in meters.
Returns
-------
n : float or ndarray
Real part of refractive index (dimensionless).
Notes
-----
For absorbing materials, this returns only Re(n). Use
get_absorption_coefficient() for the imaginary part.
"""
pass
[docs]
@abstractmethod
def get_refractive_index_gradient(
self,
x: float | NDArray[np.float64],
y: float | NDArray[np.float64],
z: float | NDArray[np.float64],
wavelength: float | NDArray[np.float64],
) -> tuple[
float | NDArray[np.float64],
float | NDArray[np.float64],
float | NDArray[np.float64],
]:
"""
Get gradient of refractive index ∇n at position (x, y, z).
Parameters
----------
x, y, z : float or ndarray
Position coordinates in meters.
wavelength : float or ndarray
Wavelength in meters.
Returns
-------
grad_n : tuple of (float or ndarray)
(∂n/∂x, ∂n/∂y, ∂n/∂z) in units of m⁻¹.
Notes
-----
For homogeneous materials, this returns (0, 0, 0).
Gradient drives ray curvature via the eikonal equation:
d²r/ds² = ∇n(r)/n(r)
References
----------
.. [1] Born & Wolf (1999), Section 3.1.
"""
pass
[docs]
@abstractmethod
def get_absorption_coefficient(
self,
x: float | NDArray[np.float64],
y: float | NDArray[np.float64],
z: float | NDArray[np.float64],
wavelength: float | NDArray[np.float64],
) -> float | NDArray[np.float64]:
"""
Get absorption coefficient α at position (x, y, z).
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⁻¹.
Notes
-----
Intensity decays as I(d) = I₀ exp(-αd) (Beer-Lambert law).
References
----------
.. [1] https://en.wikipedia.org/wiki/Beer%E2%80%93Lambert_law
"""
pass
[docs]
@abstractmethod
def get_scattering_coefficient(
self,
x: float | NDArray[np.float64],
y: float | NDArray[np.float64],
z: float | NDArray[np.float64],
wavelength: float | NDArray[np.float64],
) -> float | NDArray[np.float64]:
"""
Get scattering coefficient μ_s at position (x, y, z).
Parameters
----------
x, y, z : float or ndarray
Position coordinates in meters.
wavelength : float or ndarray
Wavelength in meters.
Returns
-------
mu_s : float or ndarray
Scattering coefficient in m⁻¹.
Notes
-----
Mean free path between scattering events: ℓ_s = 1/μ_s.
"""
pass
[docs]
def get_extinction_coefficient(
self,
x: float | NDArray[np.float64],
y: float | NDArray[np.float64],
z: float | NDArray[np.float64],
wavelength: float | NDArray[np.float64],
) -> float | NDArray[np.float64]:
"""
Get total extinction coefficient (absorption + scattering).
Parameters
----------
x, y, z : float or ndarray
Position coordinates in meters.
wavelength : float or ndarray
Wavelength in meters.
Returns
-------
mu_t : float or ndarray
Total extinction coefficient in m⁻¹.
"""
alpha = self.get_absorption_coefficient(x, y, z, wavelength)
mu_s = self.get_scattering_coefficient(x, y, z, wavelength)
return alpha + mu_s
[docs]
def get_anisotropy_factor(
self,
x: float | NDArray[np.float64],
y: float | NDArray[np.float64],
z: float | NDArray[np.float64],
wavelength: float | NDArray[np.float64],
) -> float | NDArray[np.float64]:
"""
Get scattering anisotropy factor g (Henyey-Greenstein parameter).
Parameters
----------
x, y, z : float or ndarray
Position coordinates in meters.
wavelength : float or ndarray
Wavelength in meters.
Returns
-------
g : float or ndarray
Anisotropy factor, range [-1, 1].
- g = 0: isotropic scattering
- g > 0: forward scattering
- g < 0: backward scattering
Notes
-----
Default implementation returns 0 (isotropic).
Used in Henyey-Greenstein phase function:
p(cos θ) = (1-g²) / (1 + g² - 2g cos θ)^(3/2)
References
----------
.. [1] Henyey & Greenstein (1941). Astrophysical Journal, 93, 70-83.
"""
if isinstance(x, np.ndarray):
return np.zeros_like(x)
return 0.0
[docs]
def is_homogeneous(self) -> bool:
"""
Check if material has uniform properties (no spatial variation).
Returns
-------
bool
True if material is homogeneous, False if inhomogeneous.
Notes
-----
Homogeneous materials enable optimizations (straight-line propagation).
"""
return self._is_homogeneous
# =========================================================================
# COMPATIBILITY QUERY METHODS
# =========================================================================
[docs]
@classmethod
def supported_kernels(cls) -> list[PropagationKernelID]:
"""
Return list of propagation kernels supported by this material type.
Returns
-------
list[PropagationKernelID]
List of kernel IDs this material supports.
Examples
--------
>>> from lsurf.materials import ExponentialAtmosphere
>>> ExponentialAtmosphere.supported_kernels()
[<PropagationKernelID.SIMPLE_EULER: ...>, <PropagationKernelID.SIMPLE_RK4: ...>]
"""
# Return a copy to prevent modification
return list(cls._supported_kernels)
[docs]
@classmethod
def default_kernel(cls) -> PropagationKernelID | None:
"""
Return the default propagation kernel for this material type.
Returns
-------
PropagationKernelID or None
The default kernel ID, or None if no GPU kernels are supported.
Examples
--------
>>> from lsurf.materials import ExponentialAtmosphere
>>> ExponentialAtmosphere.default_kernel()
<PropagationKernelID.SIMPLE_RK4: ...>
"""
return cls._default_kernel
[docs]
@classmethod
def supported_propagators(cls) -> list[PropagatorID]:
"""
Return list of propagators supported by this material type.
Returns
-------
list[PropagatorID]
List of propagator IDs this material supports.
Examples
--------
>>> from lsurf.materials import ExponentialAtmosphere
>>> ExponentialAtmosphere.supported_propagators()
[<PropagatorID.GPU_GRADIENT: ...>, <PropagatorID.CPU_GRADIENT: ...>]
"""
return list(cls._supported_propagators)
[docs]
@classmethod
def default_propagator(cls) -> PropagatorID | None:
"""
Return the default propagator for this material type.
Returns
-------
PropagatorID or None
The default propagator ID, or None if not specified.
Examples
--------
>>> from lsurf.materials import ExponentialAtmosphere
>>> ExponentialAtmosphere.default_propagator()
<PropagatorID.GPU_GRADIENT: ...>
"""
return cls._default_propagator
@property
def kernel_id(self) -> PropagationKernelID | None:
"""
Return the kernel ID configured for this material instance.
Returns
-------
PropagationKernelID or None
The kernel ID selected for this instance, or None if not applicable.
"""
return self._kernel_id
@property
def propagator_id(self) -> PropagatorID | None:
"""
Return the propagator ID configured for this material instance.
Returns
-------
PropagatorID or None
The propagator ID selected for this instance, or None if not applicable.
"""
return self._propagator_id
[docs]
def get_refractive_index_gradient_magnitude(
self,
x: float | NDArray[np.float64],
y: float | NDArray[np.float64],
z: float | NDArray[np.float64],
wavelength: float | NDArray[np.float64],
) -> float | NDArray[np.float64]:
"""
Get magnitude of refractive index gradient |∇n|.
Parameters
----------
x, y, z : float or ndarray
Position coordinates in meters.
wavelength : float or ndarray
Wavelength in meters.
Returns
-------
grad_mag : float or ndarray
Magnitude of gradient |∇n| in m⁻¹.
"""
grad_x, grad_y, grad_z = self.get_refractive_index_gradient(x, y, z, wavelength)
if isinstance(grad_x, np.ndarray):
return np.sqrt(grad_x**2 + grad_y**2 + grad_z**2)
else:
return math.sqrt(grad_x**2 + grad_y**2 + grad_z**2)
[docs]
def compute_phase_velocity(
self,
x: float | NDArray[np.float64],
y: float | NDArray[np.float64],
z: float | NDArray[np.float64],
wavelength: float | NDArray[np.float64],
) -> float | NDArray[np.float64]:
"""
Compute phase velocity v_p = c/n at position.
Parameters
----------
x, y, z : float or ndarray
Position coordinates in meters.
wavelength : float or ndarray
Wavelength in meters.
Returns
-------
v_p : float or ndarray
Phase velocity in m/s.
"""
c = 299792458.0 # Speed of light in vacuum (m/s)
n = self.get_refractive_index(x, y, z, wavelength)
return c / n
[docs]
def __repr__(self) -> str:
"""Return string representation."""
homo_str = "homogeneous" if self._is_homogeneous else "inhomogeneous"
return f"<{self.__class__.__name__}(name='{self.name}', {homo_str})>"