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:

.. code-block:: python

   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:

.. code-block:: python

   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
~~~~~~~~~~~~~

.. list-table::
   :header-rows: 1

   * - 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:

.. code-block:: python

   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:

.. code-block:: python

   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:

.. code-block:: python

   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:

.. code-block:: python

   from lsurf.simulation import Simulation

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

Enable logging to monitor progress:

.. code-block:: python

   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:

.. code-block:: python

   # 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
----------------

.. code-block:: python

   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:

.. code-block:: python

   # 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:

.. code-block:: python

   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:

.. code-block:: python

   # 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:

.. code-block:: python

   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")