Algorithms
==========

This section describes the computational algorithms implemented in L-SURF.

Ray-Surface Intersection
-------------------------

Overview
~~~~~~~~

Ray-surface intersection is the most critical operation in raytracing. L-SURF implements
specialized algorithms for each surface type.

Ray Parametrization
~~~~~~~~~~~~~~~~~~~

A ray is defined by:

.. math::

   \vec{r}(t) = \vec{r}_0 + t \vec{d}

where:

- :math:`\vec{r}_0` = ray origin
- :math:`\vec{d}` = ray direction (normalized)
- :math:`t` = distance parameter (:math:`t \geq 0`)

Planar Surface Intersection
~~~~~~~~~~~~~~~~~~~~~~~~~~~~

For a plane defined by point :math:`\vec{p}` and normal :math:`\vec{n}`:

**Algorithm:**

1. Compute ray-plane angle:

   .. math::

      \cos\alpha = \vec{d} \cdot \vec{n}

2. If :math:`|\cos\alpha| < \epsilon`, ray is parallel (no intersection)

3. Compute intersection distance:

   .. math::

      t = \frac{(\vec{p} - \vec{r}_0) \cdot \vec{n}}{\vec{d} \cdot \vec{n}}

4. If :math:`t < 0`, intersection is behind ray origin

5. Intersection point: :math:`\vec{r}_{\text{hit}} = \vec{r}_0 + t \vec{d}`

**Complexity:** :math:`O(1)` - constant time

**Numerical Stability:**

- Use double precision for :math:`t` calculation
- Threshold :math:`\epsilon \approx 10^{-10}` for parallelism test
- Offset intersection point by :math:`\epsilon \vec{n}` to avoid self-intersection

Spherical Surface Intersection
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

For a sphere with center :math:`\vec{c}` and radius :math:`R`:

**Implicit equation:**

.. math::

   |\vec{r} - \vec{c}|^2 = R^2

**Algorithm:**

1. Transform to sphere-local coordinates: :math:`\vec{o} = \vec{r}_0 - \vec{c}`

2. Solve quadratic equation :math:`at^2 + bt + c = 0`:

   .. math::

      a &= \vec{d} \cdot \vec{d} = 1 \\
      b &= 2 \vec{o} \cdot \vec{d} \\
      c &= \vec{o} \cdot \vec{o} - R^2

3. Compute discriminant: :math:`\Delta = b^2 - 4ac`

4. If :math:`\Delta < 0`, no intersection

5. Solutions:

   .. math::

      t_{\pm} = \frac{-b \pm \sqrt{\Delta}}{2a}

6. Choose appropriate root:

   - For exterior ray: :math:`t_1 = \min(t_+, t_-)` (first intersection)
   - For interior ray: :math:`t_2 = \max(t_+, t_-)` (exit point)

**Complexity:** :math:`O(1)` - constant time

**Numerical Stability:**

- Use numerically stable quadratic solver:

  .. math::

     q = -\frac{1}{2}(b + \text{sign}(b)\sqrt{\Delta})

     t_1 = q/a, \quad t_2 = c/q

Wave Surface Intersection
~~~~~~~~~~~~~~~~~~~~~~~~~~

For parametric wave surfaces (e.g., Gerstner waves), no closed-form solution exists.
We use iterative root-finding.

**Algorithm (Newton-Raphson):**

1. Define implicit function:

   .. math::

      f(t) = z(\vec{r}_0 + t\vec{d}) - (z_0 + t d_z)

   where :math:`z(\vec{r})` is the wave surface height.

2. Initial guess from planar approximation:

   .. math::

      t_0 = \frac{-z_0}{d_z}

3. Iterate until convergence:

   .. math::

      t_{n+1} = t_n - \frac{f(t_n)}{f'(t_n)}

4. Derivative:

   .. math::

      f'(t) = \nabla z \cdot \vec{d} - d_z

5. Convergence criterion: :math:`|f(t_n)| < \epsilon` (typically :math:`\epsilon = 10^{-6}` m)

6. Maximum iterations: 20 (typically converges in 3-5 iterations)

**Complexity:** :math:`O(N_{\text{iter}} \cdot N_{\text{waves}})` where :math:`N_{\text{waves}}` is number of wave components

**Robustness:**

- Fallback to bisection if Newton fails to converge
- Check for multiple intersections in rough seas
- Use adaptive step size for steep waves

GPU Implementation
^^^^^^^^^^^^^^^^^^

For GPU execution, wave surface intersection is parallelized:

.. code-block:: python

   @cuda.jit
   def intersect_wave_surface_kernel(ray_positions, ray_directions,
                                     wave_params, intersection_points):
       idx = cuda.grid(1)
       if idx < ray_positions.shape[0]:
           # Newton-Raphson iteration in parallel
           t = initial_guess(ray_positions[idx], ray_directions[idx])
           for iteration in range(MAX_ITERATIONS):
               pos = ray_positions[idx] + t * ray_directions[idx]
               z_surf = evaluate_wave_height(pos, wave_params)
               f = z_surf - pos[2]
               if abs(f) < EPSILON:
                   break
               df_dt = evaluate_wave_gradient(pos, wave_params).dot(
                   ray_directions[idx]
               ) - ray_directions[idx][2]
               t -= f / df_dt
           intersection_points[idx] = ray_positions[idx] + t * ray_directions[idx]

Reflection Vector Calculation
------------------------------

Law of Reflection
~~~~~~~~~~~~~~~~~

Given incident direction :math:`\vec{d}_i` and surface normal :math:`\vec{n}`:

.. math::

   \vec{d}_r = \vec{d}_i - 2(\vec{d}_i \cdot \vec{n})\vec{n}

**Algorithm:**

1. Ensure normal points toward incident medium:

   .. code-block:: python

      if dot(incident_dir, normal) > 0:
          normal = -normal

2. Compute projection: :math:`p = \vec{d}_i \cdot \vec{n}`

3. Reflect: :math:`\vec{d}_r = \vec{d}_i - 2p\vec{n}`

4. Normalize (if needed): :math:`\vec{d}_r = \vec{d}_r / |\vec{d}_r|`

**Verification:**

Check angle equality:

.. math::

   \theta_r = \arccos(-\vec{d}_r \cdot \vec{n}) = \arccos(-\vec{d}_i \cdot \vec{n}) = \theta_i

Refraction Vector Calculation
------------------------------

Snell's Law in Vector Form
~~~~~~~~~~~~~~~~~~~~~~~~~~~

Given incident direction :math:`\vec{d}_i`, surface normal :math:`\vec{n}`, and refractive
indices :math:`n_1`, :math:`n_2`:

.. math::

   \vec{d}_t = \frac{n_1}{n_2}\vec{d}_i + \left(\frac{n_1}{n_2}\cos\theta_i - \cos\theta_t\right)\vec{n}

where:

.. math::

   \cos\theta_i &= -\vec{d}_i \cdot \vec{n} \\
   \cos\theta_t &= \sqrt{1 - \left(\frac{n_1}{n_2}\right)^2(1 - \cos^2\theta_i)}

**Algorithm:**

1. Compute refractive index ratio: :math:`\eta = n_1/n_2`

2. Compute incident angle: :math:`\cos\theta_i = -\vec{d}_i \cdot \vec{n}`

3. Compute discriminant: :math:`k = 1 - \eta^2(1 - \cos^2\theta_i)`

4. Check for total internal reflection:

   .. code-block:: python

      if k < 0:
          return None  # Total internal reflection

5. Compute refracted direction:

   .. code-block:: python

      cos_theta_t = sqrt(k)
      d_t = eta * d_i + (eta * cos_theta_i - cos_theta_t) * normal

**Numerical Stability:**

- Use :math:`k = 1 - \eta^2 + \eta^2\cos^2\theta_i` form to avoid cancellation
- Threshold for TIR: :math:`k < -\epsilon` with :math:`\epsilon = 10^{-10}`

Intensity and Amplitude Tracking
---------------------------------

Ray Intensity Update
~~~~~~~~~~~~~~~~~~~~

At each surface interaction:

**Reflection:**

.. math::

   I_{\text{refl}} = I_{\text{inc}} \cdot R(\theta_i, n_1, n_2)

**Refraction:**

.. math::

   I_{\text{trans}} = I_{\text{inc}} \cdot T(\theta_i, n_1, n_2)

where :math:`R` and :math:`T` are Fresnel coefficients.

**Propagation through absorbing medium:**

.. math::

   I(d) = I_0 e^{-\alpha d}

where :math:`\alpha` is the absorption coefficient and :math:`d` is distance traveled.

Polarization State Update
~~~~~~~~~~~~~~~~~~~~~~~~~~

For polarized rays, track s and p components separately:

**Algorithm:**

1. Compute local s and p reference frames at intersection:

   .. code-block:: python

      # Plane of incidence
      plane_normal = cross(incident_dir, surface_normal)
      plane_normal /= norm(plane_normal)

      # s-direction (perpendicular to plane of incidence)
      s_dir = plane_normal

      # p-direction (in plane of incidence)
      p_dir = cross(s_dir, incident_dir)

2. Compute Fresnel coefficients: :math:`R_s, R_p, T_s, T_p`

3. Update intensities:

   .. code-block:: python

      I_s_refl = I_s_inc * R_s
      I_p_refl = I_p_inc * R_p
      I_s_trans = I_s_inc * T_s
      I_p_trans = I_p_inc * T_p

4. Rotate polarization frame for refracted ray (if needed)

Multi-Bounce Ray Tracing
-------------------------

Algorithm Overview
~~~~~~~~~~~~~~~~~~

**Input:** Initial rays, list of surfaces, max bounces

**Output:** Final ray states, bounce history

**Algorithm:**

.. code-block:: python

   function trace_rays_multi_bounce(rays, surfaces, max_bounces):
       active_rays = rays
       all_bounces = []

       for bounce in range(max_bounces):
           if no active_rays:
               break

           # Find next intersection for each ray
           intersections = []
           for ray in active_rays:
               nearest_t = infinity
               nearest_surface = None
               for surface in surfaces:
                   t, hit = surface.intersect(ray)
                   if hit and t < nearest_t:
                       nearest_t = t
                       nearest_surface = surface
               intersections.append((nearest_t, nearest_surface))

           # Update ray states
           new_rays = []
           for ray, (t, surface) in zip(active_rays, intersections):
               if surface is None:
                   continue  # Ray escaped

               # Compute interaction point
               hit_point = ray.position + t * ray.direction
               normal = surface.normal_at(hit_point)

               # Compute Fresnel coefficients
               n1 = surface.material_front.n(ray.wavelength, hit_point)
               n2 = surface.material_back.n(ray.wavelength, hit_point)
               R, T = fresnel_coefficients(ray.direction, normal, n1, n2)

               # Generate reflected ray
               if ray.intensity * R > threshold:
                   refl_dir = reflect(ray.direction, normal)
                   refl_ray = Ray(
                       position=hit_point + epsilon * normal,
                       direction=refl_dir,
                       intensity=ray.intensity * R,
                       time=ray.time + t * n1 / c
                   )
                   new_rays.append(refl_ray)

               # Generate refracted ray
               if ray.intensity * T > threshold:
                   refr_dir = refract(ray.direction, normal, n1, n2)
                   if refr_dir is not None:  # Check TIR
                       refr_ray = Ray(
                           position=hit_point - epsilon * normal,
                           direction=refr_dir,
                           intensity=ray.intensity * T,
                           time=ray.time + t * n1 / c
                       )
                       new_rays.append(refr_ray)

               all_bounces.append(BounceRecord(ray, surface, hit_point))

           active_rays = new_rays

       return active_rays, all_bounces

Computational Complexity
~~~~~~~~~~~~~~~~~~~~~~~~

For :math:`N` rays, :math:`S` surfaces, :math:`B` bounces:

- **Naïve intersection:** :math:`O(N \cdot S \cdot B)`
- **With spatial acceleration:** :math:`O(N \cdot \log S \cdot B)`

GPU parallelization reduces to:

- **Per-ray cost:** :math:`O(S \cdot B / N_{\text{threads}})`

Accuracy and Error Analysis
----------------------------

Floating Point Precision
~~~~~~~~~~~~~~~~~~~~~~~~

L-SURF uses double precision (float64) for:

- Intersection distance calculations
- Time-of-flight accumulation
- Fresnel coefficient computation

Single precision (float32) for:

- Ray positions and directions (after intersection)
- Intensity values
- Surface coordinates (GPU arrays)

Intersection Epsilon
~~~~~~~~~~~~~~~~~~~~

To avoid self-intersection, rays are offset from surfaces by :math:`\epsilon`:

- **Reflected rays:** offset in +normal direction
- **Refracted rays:** offset in -normal direction
- **Typical value:** :math:`\epsilon = 10^{-6}` m

Convergence Criteria
~~~~~~~~~~~~~~~~~~~~

**Wave surface intersection:**

- Position error: :math:`|f(t)| < 10^{-6}` m
- Angle error: :math:`< 10^{-4}` rad

**Fresnel coefficients:**

- Energy conservation: :math:`|R + T - 1| < 10^{-10}`

**Intensity tracking:**

- Termination threshold: :math:`I < 10^{-6} I_0`