Simulation Workflow

This guide explains how to set up and run ray tracing simulations with L-SURF.

Overview

A typical L-SURF simulation involves:

1. Define materials - Create or use built-in materials

2. Define surfaces - Create surfaces with roles (optical, detector, absorber)

3. Build geometry - Use GeometryBuilder to assemble the simulation

4. Configure simulation - Set propagation and termination parameters

5. Generate rays - Create initial ray batch from a source

6. Run simulation - Execute ray tracing

7. Analyze results - Process detected rays and statistics

Step 1: Define Materials

L-SURF provides built-in materials and supports custom materials:

import lsurf as sr

# Built-in homogeneous materials
air = sr.AIR_STP          # n ≈ 1.000293
water = sr.WATER          # n ≈ 1.33
glass = sr.BK7_GLASS      # n ≈ 1.52
vacuum = sr.VACUUM        # n = 1.0

# Gradient-index atmosphere (height-dependent n)
atmosphere = sr.ExponentialAtmosphere(
    n_sea_level=1.000293,
    scale_height=8500.0,  # meters
)

# Custom Sellmeier material
custom_glass = sr.create_sellmeier_material(
    B1=1.03961212, B2=0.23179234, B3=1.01046945,
    C1=0.00600069867, C2=0.02001791, C3=103.560653,
)

Step 2: Define Surfaces

Surfaces define the geometry and have roles that determine behavior:

from lsurf.surfaces import (
    PlaneSurface,
    SphereSurface,
    SurfaceRole,
    EARTH_RADIUS,
)

# Optical surface: reflects and refracts rays
ocean = SphereSurface(
    center=(0, 0, -EARTH_RADIUS),
    radius=EARTH_RADIUS,
    role=SurfaceRole.OPTICAL,
    name="ocean",
)

# Detector surface: records and terminates rays
detector = PlaneSurface(
    point=(0, 0, 35000),
    normal=(0, 0, 1),
    role=SurfaceRole.DETECTOR,
    name="detector_35km",
)

# Absorber surface: terminates rays without recording
boundary = SphereSurface(
    center=(0, 0, 0),
    radius=500000,
    role=SurfaceRole.ABSORBER,
    name="outer_boundary",
)

Surface Roles

Role	Behavior	Use Case
OPTICAL	Fresnel reflection/refraction	Physical interfaces (water, glass)
DETECTOR	Record and terminate rays	Sensors, measurement planes
ABSORBER	Terminate without recording	Boundaries, black surfaces

Step 3: Build Geometry

Use GeometryBuilder to assemble surfaces with consistent materials:

from lsurf.geometry import GeometryBuilder

geometry = (
    GeometryBuilder()
    # Register named media
    .register_medium("atmosphere", atmosphere)
    .register_medium("ocean", water)
    # Set background propagation medium
    .set_background("atmosphere")
    # Add optical surface with front/back materials
    .add_surface(ocean, front="atmosphere", back="ocean")
    # Add detector (no materials needed)
    .add_detector(detector)
    # Build immutable geometry
    .build()
)

The builder ensures:

• Material consistency across surfaces

• Validation of surface configurations

• Immutable geometry during simulation

Step 4: Configure Simulation

SimulationConfig controls simulation behavior:

from lsurf.simulation import SimulationConfig

config = SimulationConfig(
    # Propagation
    step_size=100.0,           # Max step (meters)
    min_step_size=3e-4,        # Min step for precision (0.3mm)
    adaptive_stepping=True,    # Reduce steps near surfaces
    max_steps_per_leg=10000,   # Steps before forced surface check

    # Termination
    max_bounces=5,             # Max surface interactions
    min_intensity=1e-10,       # Terminate weak rays
    bounding_radius=500000.0,  # Terminate outside this radius

    # Physics
    polarization="unpolarized",  # 's', 'p', or 'unpolarized'
    apply_absorption=True,       # Beer-Lambert absorption

    # Performance
    use_gpu=True,              # GPU acceleration if available
)

Key Parameters

step_size and min_step_size

Control the balance between speed and precision. Larger steps are faster but may miss fine surface details. Adaptive stepping automatically reduces step size near surfaces.

max_bounces

Limits simulation depth. Set high enough to capture all relevant reflections but low enough to avoid excessive computation.

bounding_radius

Defines the simulation domain. Rays escaping this radius are terminated.

Step 5: Generate Rays

Create initial rays from a source:

import lsurf as sr

# Collimated beam (parallel rays)
source = sr.CollimatedBeam(
    center=(0, 0, 1000),           # Beam center position
    direction=(0.17, 0, -0.98),    # Propagation direction
    radius=10.0,                   # Beam radius
    num_rays=50000,                # Number of rays
    wavelength=532e-9,             # Wavelength (meters)
)
rays = source.generate()

# Diverging beam (cone of rays)
source = sr.DivergingBeam(
    origin=(0, 0, 1000),
    direction=(0, 0, -1),
    half_angle=0.017,  # radians (~1°)
    num_rays=10000,
    wavelength=532e-9,
)

# Point source (spherical emission)
source = sr.PointSource(
    position=(0, 0, 1000),
    num_rays=10000,
    wavelength=532e-9,
)

Step 6: Run Simulation

Execute the simulation:

from lsurf.simulation import Simulation

# Create and run simulation
sim = Simulation(geometry, config)
result = sim.run(rays)

Enable logging to monitor progress:

import lsurf as sr

sr.configure_logging("INFO")  # Show simulation summaries
# sr.configure_logging("DEBUG")  # Show per-bounce details

result = sim.run(rays)

Step 7: Analyze Results

The SimulationResult contains all simulation outputs:

# Detected rays
detected = result.detected
print(f"Detected {detected.num_rays} rays")
print(f"Positions shape: {detected.positions.shape}")
print(f"Times: {detected.times[:5]}")
print(f"Intensities: {detected.intensities[:5]}")

# Statistics
stats = result.statistics
print(f"Initial rays: {stats.total_rays_initial}")
print(f"Created (with splits): {stats.total_rays_created}")
print(f"Detected: {stats.rays_detected}")
print(f"Absorbed: {stats.rays_absorbed}")
print(f"Out of bounds: {stats.rays_terminated_bounds}")
print(f"Below threshold: {stats.rays_terminated_intensity}")
print(f"Bounces completed: {stats.bounces_completed}")

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

Detected Ray Data

The detected object (type RecordedRays) contains:

• positions - (N, 3) array of hit positions

• directions - (N, 3) array of ray directions at hit

• times - (N,) array of arrival times (seconds)

• intensities - (N,) array of ray intensities

• wavelengths - (N,) array of wavelengths

• generations - (N,) array of ray generation (bounce count)

Complete Example

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

# Enable logging
sr.configure_logging("INFO")

# Materials
atmosphere = sr.ExponentialAtmosphere(n_sea_level=1.000293)

# Surfaces
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",
)

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

# Configure
config = SimulationConfig(
    step_size=100.0,
    max_bounces=3,
    adaptive_stepping=True,
)

# Run
sim = Simulation(geometry, config)
source = sr.CollimatedBeam(
    center=(0, 0, 2000),
    direction=(0.17, 0, -0.98),
    radius=10.0,
    num_rays=10000,
    wavelength=532e-9,
)
result = sim.run(source.generate())

# Results
print(f"Detected: {result.statistics.rays_detected}")
print(f"Detection times: {result.detected.times[:5]}")

Advanced Topics

GPU vs CPU Execution

L-SURF automatically uses GPU when available:

# Force CPU execution
config = SimulationConfig(use_gpu=False)

# Check if GPU is available
from numba import cuda
print(f"GPU available: {cuda.is_available()}")

Adaptive Stepping

Adaptive stepping provides high timing precision near surfaces while maintaining efficiency during long propagation paths:

config = SimulationConfig(
    step_size=100.0,              # Far from surfaces: 100m steps
    min_step_size=3e-4,           # Near surfaces: 0.3mm steps
    adaptive_stepping=True,
    surface_proximity_threshold=10.0,  # Start adapting within 10m
    surface_proximity_factor=0.5,      # step = distance × 0.5
)

This provides approximately 1 picosecond time resolution at surfaces.

Polarization

Enable polarization tracking for accurate Fresnel calculations:

# Single polarization state
config = SimulationConfig(polarization="s")  # or "p"

# Unpolarized (average of s and p)
config = SimulationConfig(polarization="unpolarized")

# Track 3D polarization vectors
config = SimulationConfig(
    polarization="s",
    track_polarization_vector=True,
)

Surface Hit Tracking

Track intermediate surface hits for visualization:

config = SimulationConfig(
    track_surface_hits=True,
)
result = sim.run(rays)

# Access intermediate hits
if result.surface_hits:
    for surface_name, hits in result.surface_hits.items():
        print(f"{surface_name}: {len(hits)} hit records")