# 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.
"""
Plane Surface (GPU-Capable)
Implements an infinite plane surface with GPU acceleration support.
Can be used as DETECTOR, OPTICAL, or ABSORBER.
"""
from __future__ import annotations
from dataclasses import dataclass, field
from typing import TYPE_CHECKING, Any
import numpy as np
import numpy.typing as npt
from ..protocol import Surface, SurfaceRole
from ..registry import register_surface_type
if TYPE_CHECKING:
from ...propagation.kernels.registry import IntersectionKernelID
[docs]
@dataclass
class PlaneSurface(Surface):
"""
Infinite plane surface with GPU acceleration.
The plane is defined by a point and a normal vector.
Signed distance from point p to plane: dot(p - point, normal)
Parameters
----------
point : tuple of float
A point (px, py, pz) on the plane.
normal : tuple of float
Unit normal vector (nx, ny, nz) pointing "outward" (front side).
role : SurfaceRole
What happens when a ray hits (DETECTOR, OPTICAL, or ABSORBER).
name : str
Human-readable name.
material_front : MaterialField, optional
Material on front side (where normal points). Required for OPTICAL.
material_back : MaterialField, optional
Material on back side. Required for OPTICAL.
Examples
--------
>>> # Detection plane at 35 km altitude
>>> detector = PlaneSurface(
... point=(0, 0, 35000),
... normal=(0, 0, 1),
... role=SurfaceRole.DETECTOR,
... name="altitude_detector",
... )
>>>
>>> # Glass-air interface
>>> interface = PlaneSurface(
... point=(0, 0, 0),
... normal=(0, 0, -1),
... role=SurfaceRole.OPTICAL,
... material_front=glass,
... material_back=air,
... name="glass_exit",
... )
"""
point: tuple[float, float, float]
normal: tuple[float, float, float]
role: SurfaceRole
name: str = "plane"
material_front: Any = None
material_back: Any = None
# GPU capability
_gpu_capable: bool = field(default=True, init=False, repr=False)
_geometry_id: int = field(default=1, init=False, repr=False) # plane = 1
# Internal normalized normal (set in __post_init__)
_normal: tuple[float, float, float] = field(
default=(0, 0, 1), init=False, repr=False
)
# Kernel ID for this instance (set in __post_init__)
_kernel_id: IntersectionKernelID | None = field(
default=None, init=False, repr=False
)
@classmethod
def _get_supported_kernels(cls) -> list[IntersectionKernelID]:
"""Get supported intersection kernels (lazy initialization)."""
from ...propagation.kernels.registry import IntersectionKernelID
return [IntersectionKernelID.PLANE_ANALYTICAL]
@classmethod
def _get_default_kernel(cls) -> IntersectionKernelID:
"""Get default intersection kernel."""
from ...propagation.kernels.registry import IntersectionKernelID
return IntersectionKernelID.PLANE_ANALYTICAL
[docs]
@classmethod
def supported_kernels(cls) -> list[IntersectionKernelID]:
"""Return list of intersection kernels supported by this surface type."""
return cls._get_supported_kernels()
[docs]
@classmethod
def default_kernel(cls) -> IntersectionKernelID:
"""Return the default intersection kernel for this surface type."""
return cls._get_default_kernel()
def __post_init__(self) -> None:
# Normalize the normal vector
n = np.array(self.normal, dtype=np.float64)
norm = np.linalg.norm(n)
if norm < 1e-12:
raise ValueError("Normal vector cannot be zero")
self._normal = tuple((n / norm).tolist())
# Set default kernel
self._kernel_id = self._get_default_kernel()
@property
def gpu_capable(self) -> bool:
"""This surface supports GPU acceleration."""
return True
@property
def geometry_id(self) -> int:
"""GPU geometry type ID (plane = 1)."""
return 1
[docs]
def get_gpu_parameters(self) -> tuple:
"""
Return parameters for GPU kernel.
Returns
-------
tuple
(normal_x, normal_y, normal_z, point_x, point_y, point_z)
"""
return (
self._normal[0],
self._normal[1],
self._normal[2],
self.point[0],
self.point[1],
self.point[2],
)
[docs]
def get_materials(self) -> tuple | None:
"""
Return (material_front, material_back) for Fresnel calculation.
Returns
-------
tuple or None
(material_front, material_back) or None if not OPTICAL
"""
if self.role == SurfaceRole.OPTICAL:
return (self.material_front, self.material_back)
return None
[docs]
def signed_distance(
self,
positions: npt.NDArray[np.float32],
) -> npt.NDArray[np.float32]:
"""
Compute signed distance from positions to plane.
Parameters
----------
positions : ndarray, shape (N, 3)
Points to compute distance for
Returns
-------
ndarray, shape (N,)
Signed distance (positive on normal side, negative on back side)
"""
normal = np.array(self._normal, dtype=np.float32)
point = np.array(self.point, dtype=np.float32)
# d = dot(p - point, normal)
diff = positions - point
return np.dot(diff, normal)
[docs]
def intersect(
self,
origins: npt.NDArray[np.float32],
directions: npt.NDArray[np.float32],
min_distance: float = 1e-6,
) -> tuple[npt.NDArray[np.float32], npt.NDArray[np.bool_]]:
"""
Compute ray-plane intersection.
Parameters
----------
origins : ndarray, shape (N, 3)
Ray origins
directions : ndarray, shape (N, 3)
Ray directions (normalized)
min_distance : float
Minimum valid intersection distance
Returns
-------
distances : ndarray, shape (N,)
Distance to intersection (inf if no hit)
hit_mask : ndarray, shape (N,)
Boolean mask of valid intersections
"""
normal = np.array(self._normal, dtype=np.float32)
point = np.array(self.point, dtype=np.float32)
# t = dot(point - origin, normal) / dot(direction, normal)
denom = np.dot(directions, normal)
# Parallel rays don't intersect
parallel_mask = np.abs(denom) < 1e-10
diff = point - origins
t = np.dot(diff, normal) / np.where(parallel_mask, 1.0, denom)
# Valid intersection: not parallel, t >= min_distance
hit_mask = (~parallel_mask) & (t >= min_distance)
distances = np.where(hit_mask, t, np.inf)
return distances.astype(np.float32), hit_mask
[docs]
def normal_at(
self,
positions: npt.NDArray[np.float32],
incoming_directions: npt.NDArray[np.float32] | None = None,
) -> npt.NDArray[np.float32]:
"""
Compute surface normal at positions.
For a plane, normal is constant everywhere.
Parameters
----------
positions : ndarray, shape (N, 3)
Points on the surface
incoming_directions : ndarray, shape (N, 3), optional
Ray directions (used to flip normal if needed)
Returns
-------
ndarray, shape (N, 3)
Normal vectors at each position
"""
n = len(positions)
normals = np.tile(np.array(self._normal, dtype=np.float32), (n, 1))
# Optionally flip normals to face incoming rays
if incoming_directions is not None:
dot = np.sum(normals * incoming_directions, axis=1)
flip_mask = dot > 0 # Normal facing same direction as ray
normals[flip_mask] = -normals[flip_mask]
return normals
# Register class with registry
register_surface_type("plane", "gpu", 1, PlaneSurface)