# 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.
"""
Sphere Surface (GPU-Capable)
Implements a spherical 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 SphereSurface(Surface):
"""
Sphere surface with GPU acceleration.
The sphere is defined by center and radius.
Signed distance from point p: |p - center| - radius
(positive outside, negative inside)
Parameters
----------
center : tuple of float
Center point (cx, cy, cz).
radius : float
Sphere radius (positive for convex, negative for concave).
role : SurfaceRole
What happens when a ray hits (DETECTOR, OPTICAL, or ABSORBER).
name : str
Human-readable name.
material_front : MaterialField, optional
Material outside the sphere. Required for OPTICAL.
material_back : MaterialField, optional
Material inside the sphere. Required for OPTICAL.
Examples
--------
>>> # Earth's surface as ocean
>>> EARTH_RADIUS = 6.371e6
>>> ocean = SphereSurface(
... center=(0, 0, -EARTH_RADIUS),
... radius=EARTH_RADIUS,
... role=SurfaceRole.OPTICAL,
... material_front=atmosphere,
... material_back=seawater,
... name="ocean_surface",
... )
>>>
>>> # Spherical detector
>>> detector = SphereSurface(
... center=(0, 0, 0),
... radius=1000.0,
... role=SurfaceRole.DETECTOR,
... name="spherical_detector",
... )
"""
center: tuple[float, float, float]
radius: float
role: SurfaceRole
name: str = "sphere"
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=2, init=False, repr=False) # sphere = 2
# 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.SPHERE_ANALYTICAL]
@classmethod
def _get_default_kernel(cls) -> IntersectionKernelID:
"""Get default intersection kernel."""
from ...propagation.kernels.registry import IntersectionKernelID
return IntersectionKernelID.SPHERE_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:
if self.radius == 0:
raise ValueError("Radius cannot be zero")
# 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 (sphere = 2)."""
return 2
[docs]
def get_gpu_parameters(self) -> tuple:
"""
Return parameters for GPU kernel.
Returns
-------
tuple
(center_x, center_y, center_z, radius)
"""
return (
self.center[0],
self.center[1],
self.center[2],
self.radius,
)
[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 sphere surface.
Parameters
----------
positions : ndarray, shape (N, 3)
Points to compute distance for
Returns
-------
ndarray, shape (N,)
Signed distance (positive outside, negative inside)
"""
center = np.array(self.center, dtype=np.float32)
diff = positions - center
dist = np.linalg.norm(diff, axis=1)
return (dist - abs(self.radius)).astype(np.float32)
[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-sphere 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
"""
center = np.array(self.center, dtype=np.float32)
r = abs(self.radius)
# Ray-sphere intersection:
# |origin + t*direction - center|^2 = r^2
oc = origins - center
a = np.sum(directions * directions, axis=1)
b = 2.0 * np.sum(oc * directions, axis=1)
c = np.sum(oc * oc, axis=1) - r * r
discriminant = b * b - 4 * a * c
# No intersection if discriminant < 0
no_hit = discriminant < 0
# Compute both roots
sqrt_disc = np.sqrt(np.maximum(discriminant, 0))
t1 = (-b - sqrt_disc) / (2 * a)
t2 = (-b + sqrt_disc) / (2 * a)
# Choose closest positive intersection >= min_distance
t1_valid = t1 >= min_distance
t2_valid = t2 >= min_distance
# Prefer t1 (closer) if valid, else t2
t = np.where(t1_valid, t1, np.where(t2_valid, t2, np.inf))
hit_mask = (~no_hit) & (t1_valid | t2_valid)
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 sphere, normal points radially outward from center.
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
"""
center = np.array(self.center, dtype=np.float32)
diff = positions - center
# Normalize
norms = np.linalg.norm(diff, axis=1, keepdims=True)
normals = diff / np.maximum(norms, 1e-12)
# For negative radius (concave), flip normals
if self.radius < 0:
normals = -normals
# 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
normals[flip_mask] = -normals[flip_mask]
return normals.astype(np.float32)
# Register class with registry
register_surface_type("sphere", "gpu", 2, SphereSurface)