Simulation Module

The simulation module provides the main Simulation class for running ray tracing simulations, along with SimulationConfig for configuration and SimulationResult for results.

Overview

The simulation module contains:

  1. Simulation - Main class for running ray tracing

  2. SimulationConfig - Configuration parameters

  3. SimulationResult - Results container with detected rays and statistics

  4. SimulationStatistics - Detailed statistics about the simulation run

Simulation

class lsurf.simulation.Simulation(geometry, config=None)[source]

Ray tracing simulation with geometry-based configuration.

The simulation takes a pre-built Geometry object which defines: - All optical surfaces with their materials - All detector surfaces - The background propagation medium - Named media for material consistency

Parameters:
  • geometry (Geometry) – Pre-built geometry from GeometryBuilder containing all surfaces, detectors, and materials.

  • config (SimulationConfig, optional) – Simulation configuration. Uses defaults if not provided.

Examples

>>> from lsurf.geometry import GeometryBuilder
>>> from lsurf.materials import LinsleyAtmosphere, WATER
>>> from lsurf.surfaces import SphereSurface, PlaneSurface, SurfaceRole
>>> from lsurf.simulation import Simulation, SimulationConfig
>>>
>>> # Build geometry
>>> EARTH_RADIUS = 6.371e6
>>> atmosphere = LinsleyAtmosphere()
>>>
>>> ocean = SphereSurface(
...     center=(0, 0, -EARTH_RADIUS),
...     radius=EARTH_RADIUS,
...     role=SurfaceRole.OPTICAL,
...     name="ocean",
... )
>>> detector = PlaneSurface(
...     point=(0, 0, 35000),
...     normal=(0, 0, 1),
...     role=SurfaceRole.DETECTOR,
...     name="detector_35km",
... )
>>>
>>> geometry = (
...     GeometryBuilder()
...     .register_medium("atmosphere", atmosphere)
...     .register_medium("ocean", WATER)
...     .set_background("atmosphere")
...     .add_surface(ocean, front="atmosphere", back="ocean")
...     .add_detector(detector)
...     .build()
... )
>>>
>>> # Create simulation with geometry
>>> config = SimulationConfig(step_size=100.0, max_bounces=5)
>>> sim = Simulation(geometry, config)
>>> result = sim.run(rays)
>>> print(f"Detected: {result.statistics.rays_detected}")
__init__(geometry, config=None)[source]
property geometry: Geometry

The simulation geometry.

property config: SimulationConfig

The simulation configuration.

property num_surfaces: int

Total number of surfaces (optical + absorber + detector).

property detector_surfaces: list

List of detector surfaces.

property optical_surfaces: list

List of optical surfaces.

property absorber_surfaces: list

List of absorber surfaces.

run(rays)[source]

Run the ray tracing simulation.

Parameters:

rays (RayBatch) – Initial rays to trace.

Returns:

Complete simulation results including: - detected: RecordedRays from detector surfaces - remaining: RayBatch of rays still active - statistics: SimulationStatistics with counts - detections_per_surface: dict mapping detector names to hit counts

Return type:

SimulationResult

Examples

>>> result = sim.run(rays)
>>> print(f"Detected {result.statistics.rays_detected} rays")
>>> print(f"Absorbed {result.statistics.rays_absorbed} rays")
>>> for name, count in result.detections_per_surface.items():
...     print(f"  {name}: {count} hits")
run_single_bounce(rays)[source]

Run a single propagation + interaction cycle.

Useful for step-by-step debugging or custom simulation loops.

Parameters:

rays (RayBatch) – Rays to propagate.

Returns:

  • continuing_rays (RayBatch) – Rays that should continue (reflected, refracted, no-hit).

  • result (SimulationResult) – Results from this single bounce (detections, absorptions).

Return type:

tuple[‘RayBatch’, SimulationResult]

SimulationConfig

class lsurf.simulation.SimulationConfig(step_size=100.0, min_step_size=0.0003, adaptive_stepping=True, surface_proximity_factor=0.5, surface_proximity_threshold=10.0, max_steps_per_leg=10000, max_bounces=10, min_intensity=1e-10, bounding_radius=500000.0, bounding_center=(0.0, 0.0, -6371000.0), apply_absorption=True, polarization='unpolarized', track_polarization_vector=False, track_surface_hits=False, track_refracted_rays=False, use_gpu=True, gradient_adaptive_stepping=False, target_dn_frac=0.01)[source]

Configuration for ray tracing simulation.

step_size

Maximum integration step size in meters (default 100.0).

Type:

float

min_step_size

Minimum step size in meters for adaptive stepping (default 3e-4 = 0.3mm). This provides ~1ps time resolution near surfaces.

Type:

float

adaptive_stepping

Whether to use adaptive step sizing near surfaces (default True). When enabled, steps decrease as rays approach surfaces for precise timing.

Type:

bool

surface_proximity_factor

Step size = distance * factor when within proximity threshold (default 0.5).

Type:

float

surface_proximity_threshold

Distance (in meters) within which adaptive stepping activates (default 10.0).

Type:

float

max_steps_per_leg

Maximum steps before forcing surface check (default 10000).

Type:

int

max_bounces

Maximum surface interactions before termination (default 10).

Type:

int

min_intensity

Intensity threshold below which rays are terminated (default 1e-10).

Type:

float

bounding_radius

Radius of bounding sphere in meters (default 500_000.0).

Type:

float

bounding_center

Center of bounding sphere (default (0.0, 0.0, -6.371e6) for Earth center).

Type:

tuple[float, float, float]

apply_absorption

Whether to apply Beer-Lambert absorption (default True).

Type:

bool

polarization

Polarization state: ‘s’, ‘p’, or ‘unpolarized’ (default ‘unpolarized’).

Type:

str

track_polarization_vector

Whether to track 3D polarization vectors through interactions (default False).

Type:

bool

track_surface_hits

Whether to store intermediate surface hit positions (default False). When enabled, SimulationResult.surface_hits will contain hit positions for all optical surfaces, useful for visualization.

Type:

bool

track_refracted_rays

Whether to continue propagating refracted rays from optical surfaces (default False). When False, only reflected rays continue; refracted rays are discarded. Set to False for simulations where only reflected light is of interest (e.g., ocean surface reflection where underwater propagation is not needed).

Type:

bool

use_gpu

Whether to use GPU acceleration if available (default True).

Type:

bool

Examples

>>> config = SimulationConfig(step_size=50.0, max_bounces=5)
>>> sim = Simulation(geometry, config)

Notes

Adaptive stepping provides sub-nanosecond timing precision near surfaces:

Distance to Surface | Step Size | Time Resolution |

|---------------------|———–|-----------------| | > 10m | 100m | ~333ns | | 5m | 2.5m | ~8ns | | 1m | 0.5m | ~1.7ns | | 0.1m | 0.05m | ~167ps | | < 0.6mm | 0.3mm | ~1ps (minimum) |

step_size: float = 100.0
min_step_size: float = 0.0003
adaptive_stepping: bool = True
surface_proximity_factor: float = 0.5
surface_proximity_threshold: float = 10.0
max_steps_per_leg: int = 10000
max_bounces: int = 10
min_intensity: float = 1e-10
bounding_radius: float = 500000.0
bounding_center: tuple[float, float, float] = (0.0, 0.0, -6371000.0)
apply_absorption: bool = True
polarization: str = 'unpolarized'
track_polarization_vector: bool = False
track_surface_hits: bool = False
track_refracted_rays: bool = False
use_gpu: bool = True
gradient_adaptive_stepping: bool = False
target_dn_frac: float = 0.01
__post_init__()[source]

Validate configuration.

__init__(step_size=100.0, min_step_size=0.0003, adaptive_stepping=True, surface_proximity_factor=0.5, surface_proximity_threshold=10.0, max_steps_per_leg=10000, max_bounces=10, min_intensity=1e-10, bounding_radius=500000.0, bounding_center=(0.0, 0.0, -6371000.0), apply_absorption=True, polarization='unpolarized', track_polarization_vector=False, track_surface_hits=False, track_refracted_rays=False, use_gpu=True, gradient_adaptive_stepping=False, target_dn_frac=0.01)

Configuration Parameters

The SimulationConfig dataclass accepts the following parameters:

Propagation Parameters

step_sizefloat, default=100.0

Maximum integration step size in meters. Rays advance by at most this distance per propagation step.

min_step_sizefloat, default=3e-4

Minimum step size in meters for adaptive stepping (0.3mm). This provides approximately 1 picosecond time resolution near surfaces.

adaptive_steppingbool, default=True

Whether to use adaptive step sizing near surfaces. When enabled, steps decrease as rays approach surfaces for more precise intersection timing.

surface_proximity_factorfloat, default=0.5

When within proximity threshold, step_size = distance × factor.

surface_proximity_thresholdfloat, default=10.0

Distance (meters) within which adaptive stepping activates.

max_steps_per_legint, default=10000

Maximum propagation steps before forcing a surface check.

Termination Parameters

max_bouncesint, default=10

Maximum surface interactions before ray termination.

min_intensityfloat, default=1e-10

Intensity threshold below which rays are terminated.

bounding_radiusfloat, default=500000.0

Radius of bounding sphere in meters. Rays outside are terminated.

bounding_centertuple, default=(0.0, 0.0, -6.371e6)

Center of bounding sphere (default is Earth center).

Physics Parameters

apply_absorptionbool, default=True

Whether to apply Beer-Lambert absorption during propagation.

polarizationstr, default=”unpolarized”

Polarization state: "s", "p", or "unpolarized".

track_polarization_vectorbool, default=False

Whether to track 3D polarization vectors through interactions.

Tracking Parameters

track_surface_hitsbool, default=False

Whether to store intermediate surface hit positions. Useful for visualization of ray paths.

track_refracted_raysbool, default=False

Whether to continue propagating refracted rays from optical surfaces. When False, only reflected rays continue.

Performance Parameters

use_gpubool, default=True

Whether to use GPU acceleration if available.

SimulationResult

class lsurf.simulation.SimulationResult(detected, remaining, statistics, detections_per_surface=<factory>, surface_hits=None)[source]

Complete result from a ray tracing simulation.

detected

Rays that were recorded by detector surfaces.

Type:

DetectorResult

remaining

Rays that are still active after simulation completed (either didn’t hit any surface or exceeded max bounces).

Type:

RayBatch

statistics

Detailed simulation statistics.

Type:

SimulationStatistics

detections_per_surface

Number of detections per detector surface name.

Type:

dict

surface_hits

If track_surface_hits was enabled, contains SurfaceHitRecord objects for each optical surface, keyed by surface name. Each surface may have multiple records (one per bounce). If disabled, this is None.

Type:

dict or None

Examples

>>> result = sim.run(rays)
>>> print(f"Detected {result.statistics.rays_detected} rays")
>>> print(f"Completed in {result.statistics.bounces_completed} bounces")
>>>
>>> # Access detection results
>>> print(f"Total detected intensity: {result.detected.total_intensity:.3e}")
>>> times = result.detected.times
>>> positions = result.detected.positions
>>>
>>> # Access surface hits (if track_surface_hits=True)
>>> if result.surface_hits:
...     for name, records in result.surface_hits.items():
...         for rec in records:
...             print(f"{name} bounce {rec.bounce}: {rec.num_hits} hits")
detected: DetectorResult
remaining: RayBatch
statistics: SimulationStatistics
detections_per_surface: dict[str, int]
surface_hits: dict[str, list[SurfaceHitRecord]] | None = None
property num_detected: int

Number of rays that hit detector surfaces.

property num_remaining: int

Number of rays still active.

property bounces: int

Number of bounce iterations performed (backwards compatibility).

property total_rays_processed: int

Total rays processed including splits (backwards compatibility).

get_surface_hit_positions(surface_name)[source]

Get all hit positions for a specific surface across all bounces.

Parameters:

surface_name (str) – Name of the surface.

Returns:

Concatenated hit positions from all bounces. Empty array if no hits or if track_surface_hits was not enabled.

Return type:

ndarray, shape (N, 3)

__init__(detected, remaining, statistics, detections_per_surface=<factory>, surface_hits=None)

Result Attributes

detectedRecordedRays

All rays detected by detector surfaces. Contains positions, directions, times, intensities, and wavelengths.

remainingRayBatch

Rays still active after simulation (not detected/absorbed/terminated).

statisticsSimulationStatistics

Detailed statistics about the simulation run.

detections_per_surfacedict[str, int]

Mapping from detector name to number of hits.

surface_hitsdict[str, list] or None

Intermediate surface hit positions (if track_surface_hits=True).

SimulationStatistics

class lsurf.simulation.SimulationStatistics(total_rays_initial=0, total_rays_created=0, rays_detected=0, rays_absorbed=0, rays_terminated_intensity=0, rays_terminated_bounds=0, rays_terminated_max_bounces=0, bounces_completed=0, max_depth_reached=0, adaptive_step_stats=None)[source]

Statistics from a simulation run.

total_rays_initial

Number of rays at simulation start.

Type:

int

total_rays_created

Total rays including split rays.

Type:

int

rays_detected

Rays that hit detector surfaces.

Type:

int

rays_absorbed

Rays terminated by absorber surfaces.

Type:

int

rays_terminated_intensity

Rays terminated due to low intensity.

Type:

int

rays_terminated_bounds

Rays that exited bounding sphere.

Type:

int

rays_terminated_max_bounces

Rays that reached max bounce limit.

Type:

int

bounces_completed

Number of bounce iterations completed.

Type:

int

max_depth_reached

Maximum tree depth from ray splitting.

Type:

int

total_rays_initial: int = 0
total_rays_created: int = 0
rays_detected: int = 0
rays_absorbed: int = 0
rays_terminated_intensity: int = 0
rays_terminated_bounds: int = 0
rays_terminated_max_bounces: int = 0
bounces_completed: int = 0
max_depth_reached: int = 0
adaptive_step_stats: AdaptiveStepStats | None = None
__init__(total_rays_initial=0, total_rays_created=0, rays_detected=0, rays_absorbed=0, rays_terminated_intensity=0, rays_terminated_bounds=0, rays_terminated_max_bounces=0, bounces_completed=0, max_depth_reached=0, adaptive_step_stats=None)

Statistics Attributes

total_rays_initialint

Number of rays at simulation start.

total_rays_createdint

Total rays including those created by splitting.

rays_detectedint

Rays that hit detector surfaces.

rays_absorbedint

Rays absorbed by absorber surfaces.

rays_terminated_boundsint

Rays terminated for leaving bounding sphere.

rays_terminated_intensityint

Rays terminated for falling below intensity threshold.

rays_terminated_max_bouncesint

Rays terminated for exceeding max_bounces.

bounces_completedint

Number of bounce iterations completed.

max_depth_reachedint

Maximum ray generation depth reached.

Usage Example

from lsurf.geometry import GeometryBuilder
from lsurf.simulation import Simulation, SimulationConfig
from lsurf.surfaces import SphereSurface, PlaneSurface, SurfaceRole
import lsurf as sr

# Build geometry (see geometry module)
geometry = (
    GeometryBuilder()
    .register_medium("atmosphere", sr.ExponentialAtmosphere())
    .register_medium("ocean", sr.WATER)
    .set_background("atmosphere")
    .add_surface(ocean, front="atmosphere", back="ocean")
    .add_detector(detector)
    .build()
)

# Configure simulation
config = SimulationConfig(
    step_size=100.0,
    max_bounces=5,
    adaptive_stepping=True,
    min_step_size=3e-4,  # 0.3mm for ~1ps timing
    polarization="unpolarized",
)

# Create simulation
sim = Simulation(geometry, config)

# Generate rays
source = sr.CollimatedBeam(
    center=(0, 0, 1000),
    direction=(0.17, 0, -0.98),
    radius=10.0,
    num_rays=50000,
    wavelength=532e-9,
)
rays = source.generate()

# Run simulation
result = sim.run(rays)

# Process results
print(f"Detected: {result.statistics.rays_detected}")
print(f"Absorbed: {result.statistics.rays_absorbed}")
print(f"Bounces: {result.statistics.bounces_completed}")

# Access detected rays
detected = result.detected
print(f"Detection times: {detected.times[:10]}")
print(f"Detection positions: {detected.positions[:10]}")

# Per-detector breakdown
for name, count in result.detections_per_surface.items():
    print(f"  {name}: {count} hits")

Adaptive Stepping

The adaptive stepping feature provides high timing precision near surfaces while maintaining efficiency during long propagation paths:

Adaptive Step Size vs Time Resolution

Distance to Surface

Step Size

Time Resolution

> 10m

100m

~333ns

5m

2.5m

~8ns

1m

0.5m

~1.7ns

0.1m

0.05m

~167ps

< 0.6mm

0.3mm (min)

~1ps

Logging

Enable logging to monitor simulation progress:

import lsurf as sr

# INFO level: simulation summaries
sr.configure_logging("INFO")

# DEBUG level: per-bounce details
sr.configure_logging("DEBUG")

Example output at INFO level:

INFO - Simulation initialized: 1 optical, 1 detector, 0 absorber surfaces
INFO - Starting simulation: 50000 rays, max 5 bounces, step_size=100 m
INFO - Simulation complete: 2 bounces, 48523 detected, 0 absorbed, 1477 remaining
INFO -   Detector 'detector_35km': 48523 hits