# 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.
"""
Recording Sphere Surface (GPU-Capable)
Specialized spherical detector surface for recording rays at a specific
altitude above Earth. Follows the Surface protocol with GPU acceleration.
This is the primary detector setup for Earth-scale simulations.
"""
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
if TYPE_CHECKING:
from ...propagation.kernels.registry import IntersectionKernelID
# Earth's mean radius in meters
EARTH_RADIUS = 6.371e6
[docs]
@dataclass
class RecordingSphereSurface(Surface):
"""
Recording sphere surface at a specified altitude above Earth.
A specialized detector surface designed for Earth-scale simulations.
The sphere is centered at Earth's center with radius = earth_radius + altitude.
This is the recommended detector setup for atmospheric ray tracing.
Parameters
----------
altitude : float
Altitude above Earth's surface in meters (default 33 km).
earth_center : tuple of float
Center of Earth in simulation coordinates.
Default (0, 0, -EARTH_RADIUS) places Earth surface at z=0.
earth_radius : float
Earth radius in meters. Default is 6.371e6 m.
name : str
Human-readable name for the detector.
Examples
--------
>>> # Create a recording sphere at 33 km altitude
>>> detector = RecordingSphereSurface(
... altitude=33000.0,
... name="satellite_detector",
... )
>>>
>>> # Use in geometry builder
>>> geometry = (
... GeometryBuilder()
... .register_medium("atmosphere", atmosphere)
... .set_background("atmosphere")
... .add_detector(detector)
... .build()
... )
>>>
>>> # Check GPU capability
>>> print(detector.gpu_capable) # True
>>> print(detector.geometry_id) # 2 (sphere)
Notes
-----
This surface is GPU-capable and uses the same GPU kernels as SphereSurface.
The geometry_id is 2 (sphere), which uses analytical ray-sphere intersection.
The coordinate system places:
- Earth's center at (0, 0, -EARTH_RADIUS) by default
- Earth's surface at z=0
- Recording sphere at z = altitude at the pole
"""
altitude: float = 33000.0
earth_center: tuple[float, float, float] = (0, 0, -EARTH_RADIUS)
earth_radius: float = EARTH_RADIUS
name: str = "recording_sphere"
# Role is always DETECTOR
role: SurfaceRole = field(default=SurfaceRole.DETECTOR, init=False)
# No materials needed for detector
material_front: Any = field(default=None, init=False)
material_back: Any = field(default=None, init=False)
# GPU capability - same as SphereSurface
_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
_kernel_id: "IntersectionKernelID | None" = field(
default=None, init=False, repr=False
)
# Computed property
_sphere_radius: float = field(default=0.0, init=False, repr=False)
_center_array: npt.NDArray[np.float64] = field(
default_factory=lambda: np.zeros(3), 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()
[docs]
def __post_init__(self) -> None:
"""Initialize computed properties."""
if self.altitude < 0:
raise ValueError(f"Altitude must be non-negative, got {self.altitude}")
# Compute sphere radius
self._sphere_radius = self.earth_radius + self.altitude
# Convert center to array
self._center_array = np.array(self.earth_center, dtype=np.float64)
# 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
@property
def sphere_radius(self) -> float:
"""Total radius of the recording sphere (earth_radius + altitude)."""
return self._sphere_radius
@property
def center(self) -> tuple[float, float, float]:
"""Center of the sphere (same as earth_center)."""
return self.earth_center
@property
def center_array(self) -> npt.NDArray[np.float64]:
"""Center of the sphere as numpy array."""
return self._center_array
[docs]
def get_gpu_parameters(self) -> tuple:
"""
Return parameters for GPU kernel.
Returns
-------
tuple
(center_x, center_y, center_z, radius)
"""
return (
self.earth_center[0],
self.earth_center[1],
self.earth_center[2],
self._sphere_radius,
)
[docs]
def get_materials(self) -> tuple | None:
"""
Return materials for Fresnel calculation.
Returns None since DETECTOR surfaces don't use Fresnel equations.
"""
return None
[docs]
def signed_distance(
self,
positions: npt.NDArray[np.float32],
) -> npt.NDArray[np.float32]:
"""
Compute signed distance from positions to recording sphere surface.
Parameters
----------
positions : ndarray, shape (N, 3)
Points to compute distance for
Returns
-------
ndarray, shape (N,)
Signed distance (positive outside, negative inside)
"""
center = self._center_array.astype(np.float32)
diff = positions - center
dist = np.linalg.norm(diff, axis=1)
return (dist - self._sphere_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 = self._center_array.astype(np.float32)
r = self._sphere_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 = self._center_array.astype(np.float32)
diff = positions - center
# Normalize
norms = np.linalg.norm(diff, axis=1, keepdims=True)
normals = diff / np.maximum(norms, 1e-12)
# 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)
[docs]
def compute_angular_coordinates(
self,
positions: npt.NDArray[np.float32],
) -> dict[str, npt.NDArray[np.float32]]:
"""
Compute angular coordinates for intersection points on the sphere.
Computes spherical coordinates (latitude/longitude) of points
on the detection sphere relative to Earth's center.
Parameters
----------
positions : ndarray, shape (N, 3)
Intersection positions on the recording sphere
Returns
-------
dict
Dictionary with:
- 'latitude': Latitude angle (-pi/2 to pi/2)
- 'longitude': Longitude angle (-pi to pi)
- 'colatitude': Colatitude angle (0 to pi)
"""
# Vector from Earth center to position
to_pos = positions.astype(np.float64) - self._center_array
r = np.linalg.norm(to_pos, axis=1, keepdims=True)
# Latitude: angle above equatorial plane
latitude = np.arcsin(np.clip(to_pos[:, 2] / r.squeeze(), -1.0, 1.0))
# Longitude: angle in XY plane
longitude = np.arctan2(to_pos[:, 1], to_pos[:, 0])
# Colatitude: angle from +Z axis
colatitude = np.arccos(np.clip(to_pos[:, 2] / r.squeeze(), -1.0, 1.0))
return {
"latitude": latitude.astype(np.float32),
"longitude": longitude.astype(np.float32),
"colatitude": colatitude.astype(np.float32),
}
[docs]
def __repr__(self) -> str:
"""Return string representation."""
return (
f"RecordingSphereSurface("
f"altitude={self.altitude/1000:.1f}km, "
f"name='{self.name}', GPU)"
)
# Also export as LocalRecordingSphereSurface for local-scale simulations
[docs]
@dataclass
class LocalRecordingSphereSurface(Surface):
"""
Local recording sphere surface centered at an arbitrary position.
A simplified recording sphere for local-scale simulations without
Earth curvature considerations.
Parameters
----------
radius : float
Sphere radius in meters (default 33 km).
center : tuple of float
Center position (default (0, 0, 0)).
name : str
Human-readable name for the detector.
Examples
--------
>>> # Create a local recording sphere at origin
>>> detector = LocalRecordingSphereSurface(
... radius=1000.0,
... center=(0, 0, 0),
... name="local_detector",
... )
"""
radius: float = 33000.0
center: tuple[float, float, float] = (0, 0, 0)
name: str = "local_recording_sphere"
# Role is always DETECTOR
role: SurfaceRole = field(default=SurfaceRole.DETECTOR, init=False)
# No materials needed for detector
material_front: Any = field(default=None, init=False)
material_back: Any = field(default=None, init=False)
# 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
_kernel_id: "IntersectionKernelID | None" = field(
default=None, init=False, repr=False
)
# Computed property
_center_array: npt.NDArray[np.float64] = field(
default_factory=lambda: np.zeros(3), init=False, repr=False
)
@classmethod
def _get_supported_kernels(cls) -> list["IntersectionKernelID"]:
"""Get supported intersection kernels."""
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()
[docs]
def __post_init__(self) -> None:
"""Initialize computed properties."""
if self.radius <= 0:
raise ValueError(f"Radius must be positive, got {self.radius}")
self._center_array = np.array(self.center, dtype=np.float64)
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
@property
def sphere_radius(self) -> float:
"""Radius of the recording sphere."""
return self.radius
@property
def center_array(self) -> npt.NDArray[np.float64]:
"""Center as numpy array."""
return self._center_array
[docs]
def get_gpu_parameters(self) -> tuple:
"""Return parameters for GPU kernel."""
return (
self.center[0],
self.center[1],
self.center[2],
self.radius,
)
[docs]
def get_materials(self) -> tuple | None:
"""Return None (detectors don't use Fresnel)."""
return None
[docs]
def signed_distance(
self,
positions: npt.NDArray[np.float32],
) -> npt.NDArray[np.float32]:
"""Compute signed distance from positions to sphere surface."""
center = self._center_array.astype(np.float32)
diff = positions - center
dist = np.linalg.norm(diff, axis=1)
return (dist - 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."""
center = self._center_array.astype(np.float32)
r = self.radius
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_hit = discriminant < 0
sqrt_disc = np.sqrt(np.maximum(discriminant, 0))
t1 = (-b - sqrt_disc) / (2 * a)
t2 = (-b + sqrt_disc) / (2 * a)
t1_valid = t1 >= min_distance
t2_valid = t2 >= min_distance
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."""
center = self._center_array.astype(np.float32)
diff = positions - center
norms = np.linalg.norm(diff, axis=1, keepdims=True)
normals = diff / np.maximum(norms, 1e-12)
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)
[docs]
def __repr__(self) -> str:
"""Return string representation."""
return (
f"LocalRecordingSphereSurface("
f"radius={self.radius/1000:.1f}km, "
f"center={self.center}, "
f"name='{self.name}', GPU)"
)
# Note: We don't register these classes with the registry because they use
# the same geometry_id (2 = sphere) as SphereSurface. They share GPU kernels
# but provide a specialized interface for recording sphere functionality.