Fix reflector phases and ghost path losses

docs/architecture.md

@@ -55,8 +55,8 @@ Do not retain completed acquisitions in a step: that would defeat streaming.
 
 | Stage | Main fields it supplies or updates |
 | --- | --- |
-| `FiberGenerator` | `z`, `dz`, `fiber_profile`, `attenuation_envelope` |
-| Strain and temperature | `fiber_profile`, applied perturbation fields |
+| `FiberGenerator` | `z`, `dz`, `fiber_profile`, `attenuation_envelope`, `ghost_profile` |
+| Strain and temperature | `fiber_profile`, `ghost_profile`, applied perturbation fields |
 | `SweptLaser` | `t`, `dt`, `n_samples`, `E_source`, frequency arrays |
 | `MachZehnder` | `photocurrent_main` (optical beat power, despite the name) |
 | `AuxMZI` | `aux_signal`, `aux_valid_start` when enabled |
@@ -72,7 +72,9 @@ The HDF5 writer in `output/` is a thin layer on top.
 
 Array fields start as `None` until the relevant stage supplies them. Main-channel
 arrays have shape `(n_cores, n_samples)`; `fiber_profile` has shape
-`(n_cores, n_z)`. The laser and auxiliary signal are shared across cores and
+`(n_cores, n_z)`, as does the optional `ghost_profile`. Ghosts already include
+their path losses and are added after weighting the direct profile by
+`attenuation_envelope`. The laser and auxiliary signal are shared across cores and
 use one-dimensional arrays. Call `acq.add_log(name, ...)` for per-sweep
 metadata. Adding a dataclass field alone does not save it to HDF5: datasets
 are written explicitly in `output/hdf5_writer.py`.

docs/config_reference.md

@@ -105,7 +105,13 @@ outside the table clamp to the edge with a warning).
 | --- | --- | --- | --- |
 | `z` | length | -- | position along the fiber |
 | `R` | float in [0, 1] | -- | power reflectivity |
-| `loss_dB` | float >= 0 | 0.0 | one-way insertion loss at the reflector |
+| `loss_dB` | float >= 0 | 0.0 | extra one-way loss on the input side of the reflecting interface |
+
+Each reflector gets a seeded phase per core, fixed across sweeps. Direct
+returns beyond it include the round-trip field transmission `1-R` as well as
+`loss_dB`. Reflectors must occupy distinct bins after spatial rounding.
+`loss_dB` excludes the power removed by reflection; see the phase and loss
+conventions in [the theory guide](theory.md#discrete-reflectors-and-connectors).
 
 ### fiber.fbg_arrays[]
 
@@ -158,7 +164,7 @@ Full coupled-mode equations in #69.
 ### fiber.multiple_scattering
 
 Cascading multi-bounce ghosts between discrete reflectors (#35). Order N
-enumerates paths with 2N+1 reflection events; N=2 is the dominant case
+enumerates paths with 2N-1 reflection events; N=2 is the dominant case
 (e.g. input connector + end face -> ghost at `2*z_b - z_a`).
 
 | field | type | default | description |

docs/theory.md

@@ -12,7 +12,7 @@ For a perfectly linear sweep the beat frequency of a scatterer at distance
 index, `gamma = delta_nu / T_sweep` the chirp rate, factor 2 for the round
 trip. The whole time-domain beat is then the inverse transform of the
 attenuation weighted fiber profile:
-    s(t) = IFFT[ profile(z) * attenuation(z) ] * n_samples
+    s(t) = IFFT[ profile(z) * attenuation(z) + ghosts(z) ] * n_samples
 
 (`optics/mach_zehnder.py`). O(N log N) instead of the O(N^2) of a
 per-scatterer sum, and exact while the chirp stays linear. Mode hops and
@@ -91,15 +91,73 @@ Ref: Marcuse, "Curvature loss formula for optical fibers", JOSA 66, 216
 (1976).
 
 ### Discrete reflectors and connectors
-A point reflector just adds `sqrt(R)` into its bin of the profile (`fiber/reflectors.py`). A connector insertion loss is a multiplicative step on the envelope, `envelope[idx:] *= 10^(-loss_dB/10)`. Exponent /10 and not /20 because of the round-trip convention above.
+A point reflector adds `r = sqrt(R) exp(j phi)` at its nearest spatial bin
+(`fiber/reflectors.py`). Each core has its own seeded reflector phases, drawn
+once and reused across sweeps. This stream is separate from the Rayleigh RNG.
+Reflectors are sorted by position before drawing phases. Two entries that
+round to the same bin are rejected: their internal ordering is not resolved.
+
+We choose real field transmission `t = sqrt(1-R)`. Reflection from the right
+is then `-conj(r)`, so the lossless interface has the reciprocal, unitary
+matrix `[[r, t], [t, -conj(r)]]`. This fixes the relative phases of direct
+returns and ghosts; the phases of separate bounces are not drawn independently.
+For the underlying two-port constraints, see the scattering-matrix discussion
+in [Karlsson's electrodynamics notes, section 3.7](https://www.eit.lth.se/fileadmin/eit/courses/eitn80/root.pdf).
+
+`loss_dB` is an extra matched loss just before the interface, on its input
+(smaller-z) side. Writing `q = 10^(-loss_dB/20)`, the left reflection is
+`q^2 r`, right reflection is `-conj(r)`, and transmission is `q t` in either
+direction. This preserves the earlier convention that a direct reflection
+sees its own insertion loss twice. It also avoids counting the reflected
+fraction twice: `loss_dB` is additional to the loss of transmitted power `R`.
+
+The direct-return envelope includes `1-R` beyond each reflector and
+`10^(-loss_dB/10)` from its bin onwards. The reflector's own return does not
+cross its interface. Rayleigh and weak-FBG returns at that same bin use the
+upstream-side convention. With `R=1`, all later direct returns vanish.
 
 ### Multiple-scattering ghosts
-Light reflected by one strong reflector can bounce off another and come back via a longer path, showing up as a ghost peak
-(`fiber/multiple_scattering.py`). For an order-N path (2N-1 reflection events) the apparent depth and amplitude are
+Light can bounce between discrete reflectors and return at a later delay
+(`fiber/multiple_scattering.py`). An order-N path has `2N-1` reflections and
+apparent depth
+
     z_app = z_i0 - z_i1 + z_i2 - ... + z_i(2N-2)
-    a     = prod sqrt(R_im)
-with the constraint that every odd "valley" reflector sits strictly below both neighbours. Ghosts whose apparent bin falls outside `[0, n_z)`
-are dropped, extending the array to 2L would touch too much downstream code. The continuum Rayleigh-Rayleigh double-scatter floor is much weaker in standard SMF and tracked separately in #79.
+
+Every odd "valley" must sit strictly below both neighbours. Positions are
+rounded to bins before forming this sum. Each path multiplies the directional
+reflection coefficients and the transmission of every intervening interface,
+including the initial trip from the input and the final trip back.
+
+Distributed loss is integrated from the input to each bin. Its value at zero
+is zero; the loss assigned to bin k applies to `[z[k], z[k+1])`. For each path
+leg, the difference of these cumulative losses gives the one-way loss of that
+leg. Summing the legs and converting with `10^(-path_loss_dB/20)` gives the
+field attenuation. Segments and bends therefore count on every crossing.
+For uniform loss this agrees with a total distance `2*z_app`; for nonuniform
+loss, the visited sections matter.
+
+For example, the first ghost of two lossless interfaces at `a < b` is
+`-(1-R_a) conj(r_a) r_b^2` at `2*b-a`, before distributed attenuation.
+The two paths `(b,a,c)` and `(c,a,b)` share a delay and add coherently.
+
+`Acquisition.ghost_profile` stores these already attenuated returns separately.
+The optics step adds them to `fiber_profile * attenuation_envelope`, then
+applies the usual coherence roll-off and beat transform. Dividing by the
+direct envelope would fail behind an opaque reflector or an underflowed loss,
+even when a shorter physical path still produces a ghost at that delay.
+
+The expansion still stops at `max_order` (2 or 3), and out-of-range ghosts
+are dropped. It is a finite-bounce model, not a converged high-finesse cavity
+solver. Index, strain and temperature perturbations still use the phase at
+the apparent bin, so nonuniform perturbations along a ghost path remain an
+approximation. Rayleigh-mediated multiple scattering is outside this model.
+
+`tests/test_reflector_paths.py` checks complex two-mirror amplitudes, repeated
+crossings, localized loss, opaque interfaces, sweep reproducibility and the
+optics sum. A separate frequency-domain recursion sums the infinite series
+for three weak mirrors; its impulse response must agree within `1e-10` field
+amplitude, where omitted orders are below that tolerance. These are numerical
+checks of the scalar model, not experimental validation.
 
 ### Strain
 A scatterer at z picks up an extra round-trip phase from the integrated strain (`fiber/strain.py`):

src/pyofdr/core/acquisition.py

@@ -16,6 +16,7 @@ class Acquisition:
     # per-core arrays have shape (n_cores, ...). single-core runs use (1, ...).
     fiber_profile: np.ndarray | None = None       # (n_cores, n_z), complex
     attenuation_envelope: np.ndarray | None = None  # (n_z,), float
+    ghost_profile: np.ndarray | None = None       # (n_cores, n_z), already attenuated
 
     t: np.ndarray | None = None       # (n_t,), time samples [s]
     dt: float = 0.0                   # sample interval [s]

src/pyofdr/core/config_models.py

@@ -48,7 +48,7 @@ class SimulationConfig(_Strict):
 class ReflectorEntry(_Strict):
     z: Length
     R: float = Field(..., ge=0, le=1)
-    loss_dB: float = Field(0.0, ge=0)   # one-way insertion loss [dB]
+    loss_dB: float = Field(0.0, ge=0)   # extra one-way loss on the input side [dB]
 
 
 class FBGEntry(_Strict):
@@ -161,7 +161,7 @@ class CrosstalkConfig(_Strict):
 class MultipleScatteringConfig(_Strict):
     """Cascading multi-bounce ghosts between discrete reflectors (#35).
 
-    Order N enumerates paths with 2N+1 reflection events; the dominant
+    Order N enumerates paths with 2N-1 reflection events; the dominant
     case is N=2 (e.g. input connector + end face -> ghost at 2*z_b-z_a).
     Continuum Rayleigh-Rayleigh background is tracked in #79 -- it sits
     well below shot/RIN/thermal in standard SMF.

src/pyofdr/fiber/attenuation.py

@@ -37,8 +37,14 @@ def round_trip_attenuation_varying(z, alpha_dB_km_z, dz, xp=np):
     per spatial bin. Left-Riemann cumulative integral, so the per-bin
     discretization error is O(dz * d(alpha)/dz). For typical SMF
     (alpha ~ 0.2 dB/km, dz ~ mm) this is ~1e-6 dB, well below any
-    realistic loss budget. Step discontinuities in alpha land on the
-    left of the jumping bin rather than straddling it.
+    realistic loss budget. The loss at bin k applies to [z[k], z[k+1]).
     """
-    alpha_np = dB_per_km_to_neper_per_m(alpha_dB_km_z)
-    return xp.exp(-xp.cumsum(alpha_np) * dz)
+    loss_dB = cumulative_loss_dB(alpha_dB_km_z, dz, xp=xp)
+    return xp.exp(-loss_dB * _LN10 / 10.0)
+
+
+def cumulative_loss_dB(alpha_dB_km_z, dz, xp=np):
+    """One-way distributed loss from the input to each bin, in dB."""
+    loss = xp.zeros_like(alpha_dB_km_z, dtype=float)
+    loss[1:] = xp.cumsum(alpha_dB_km_z[:-1]) * (dz / 1000.0)
+    return loss

src/pyofdr/fiber/multiple_scattering.py

@@ -1,104 +1,84 @@
-"""Cascading multi-bounce ghost reflections from discrete reflectors (#35).
-
-Light reflected by one strong reflector can be re-scattered by another and
-end up back at the detector via a longer optical path. Each such path
-appears in the OFDR signal as a ghost peak at an apparent depth set by
-the round-trip geometry.
-
-For an order-N path (2N-1 reflection events) with reflector indices
-(i_0, i_1, ..., i_{2N-2}) the apparent depth and amplitude are
-
-    z_app = z_{i_0} - z_{i_1} + z_{i_2} - ... + z_{i_{2N-2}}
-    a     = prod sqrt(R_{i_m})
-
-with the trajectory constraint that every odd-indexed reflector (the
-"valleys" where the light bounces back forward) sits strictly below both
-of its neighbours: z_{i_{2m+1}} < z_{i_{2m}} and z_{i_{2m+1}} < z_{i_{2m+2}}.
-
-The continuum Rayleigh-Rayleigh background is several orders of magnitude
-weaker in standard SMF -- tracked separately in #79.
-"""
+"""Multi-bounce returns from discrete reflectors."""
 
 from __future__ import annotations
 
 import itertools
 import math
-from typing import Any
 
 import numpy as np
 
+from pyofdr.fiber.reflectors import reflector_bin
 
-def add_ghost_reflections(
-    profile: np.ndarray,
-    dz: float,
-    reflectors: list[dict[str, Any]],
-    *,
-    max_order: int = 2,
-    xp=np,
-) -> np.ndarray:
-    """Inject multi-bounce ghost peaks on top of the existing profile.
-
-    Ghosts whose apparent bin falls outside [0, n_z) are dropped --
-    extending the array to 2L would change a lot of downstream code,
-    so for now we only model ghosts that land within the fiber range.
-    For the canonical input-connector + end-face case the user can pad
-    the fiber length so 2L_eff stays inside the array.
-
-    Bin convention: each reflector z is rounded to its bin first, then
-    the apparent bin is the signed sum of those bins (NOT round(z_app/dz)
-    of the real-valued sum -- the two differ by at most 1).
 
-    Caveat: if an interior "valley" reflector carries a non-zero loss_dB,
-    the per-bin attenuation envelope undercounts the actual loss along
-    a ghost path (the ghost crosses the valley twice). For typical OFDR
-    setups -- where the loss reflector is the input connector at z=0 --
-    this is exact.
+def add_ghost_reflections(
+    profile, dz, reflectors, *, max_order=2, reflection=None,
+    loss_to_z=None, xp=np,
+):
+    """Add ghosts with all path losses already included, in-place.
+
+    Order N has 2N-1 reflections. Positions are snapped before forming
+    the apparent bin; ghosts outside the array are dropped.
+    loss_to_z is the cumulative one-way distributed loss in dB, without
+    connector losses. reflection contains the left-side coefficients,
+    with shape (n_cores, n_reflectors).
     """
     if max_order < 2 or len(reflectors) < 2:
         return profile
 
-    n_z = profile.shape[-1]
-
-    def _bin(r):
-        b = int(round(r["z"] / dz))
-        if b == n_z and r["z"] <= n_z * dz:
-            b = n_z - 1   # end-face clamp, must match inject_reflectors (#84)
-        return b
-
-    # drop reflectors whose bin sits outside [0, n_z): they were skipped
-    # by inject_reflectors so they don't actually reflect, can't seed ghosts
-    in_range = [r for r in reflectors if 0 <= _bin(r) < n_z]
-    if len(in_range) < 2:
-        return profile
-    bins = [_bin(r) for r in in_range]
-    amps = [math.sqrt(r["R"]) for r in in_range]
-    n_r = len(in_range)
-
-    # enumerate orders 2..max_order. each order N has 2N-1 reflection
-    # events; itertools.product is in C so the loop body is the only
-    # thing that runs in Python.
-    for N in range(2, max_order + 1):
-        path_len = 2 * N - 1
-        for path in itertools.product(range(n_r), repeat=path_len):
-            # valley constraint -- every odd-position reflector must sit
-            # strictly below both neighbours
-            valid = True
-            for m in range(1, path_len, 2):
-                v = bins[path[m]]
-                if v >= bins[path[m - 1]] or v >= bins[path[m + 1]]:
-                    valid = False
-                    break
-            if not valid:
+    n_cores, n_z = profile.shape
+    selected = sorted(
+        (i for i, r in enumerate(reflectors)
+         if 0 <= reflector_bin(r, dz, n_z) < n_z),
+        key=lambda i: reflectors[i]["z"],
+    )
+    refs = [reflectors[i] for i in selected]
+    bins = [reflector_bin(r, dz, n_z) for r in refs]
+    if len(set(bins)) != len(bins):
+        raise ValueError("discrete reflectors must occupy distinct spatial bins")
+    if reflection is None:
+        reflection = xp.tile(xp.asarray([math.sqrt(r["R"]) for r in refs]),
+                             (n_cores, 1))
+    else:
+        reflection = reflection[:, selected]
+
+    # Node 0 is the input, before a possible reflector at z=0.
+    losses = [0.0] + [0.0 if loss_to_z is None else float(loss_to_z[b])
+                      for b in bins]
+    extra = [10.0 ** (-r.get("loss_dB", 0.0) / 20.0) for r in refs]
+    transmission = [math.sqrt(1.0 - r["R"]) * q for r, q in zip(refs, extra)]
+    n_r = len(refs)
+    travel = [[1.0] * (n_r + 1) for _ in range(n_r + 1)]
+    for i in range(n_r + 1):
+        for j in range(i + 1, n_r + 1):
+            factor = 10.0 ** (-abs(losses[j] - losses[i]) / 20.0)
+            for k in range(i + 1, j):
+                factor *= transmission[k - 1]
+            travel[i][j] = travel[j][i] = factor
+
+    active = [i for i, r in enumerate(refs) if r["R"] > 0]
+    for order in range(2, max_order + 1):
+        for path in itertools.product(active, repeat=2 * order - 1):
+            if any(bins[path[m]] >= bins[path[m - 1]]
+                   or bins[path[m]] >= bins[path[m + 1]]
+                   for m in range(1, len(path), 2)):
                 continue
-            # signed sum -> apparent depth bin
-            b = 0
-            for i, idx in enumerate(path):
-                b += bins[idx] if i % 2 == 0 else -bins[idx]
-            if b < 0 or b >= n_z:
+            b = sum(bins[idx] if m % 2 == 0 else -bins[idx]
+                    for m, idx in enumerate(path))
+            if not 0 <= b < n_z:
                 continue
-            amp = 1.0
-            for idx in path:
-                amp *= amps[idx]
+
+            amp = xp.ones(n_cores, dtype=xp.complex128)
+            previous = 0
+            for m, idx in enumerate(path):
+                amp *= travel[previous][idx + 1]
+                if m % 2 == 0:
+                    # The extra loss sits on the left of the interface.
+                    amp *= reflection[:, idx] * extra[idx] ** 2
+                else:
+                    # Real transmission fixes the reverse reflection phase.
+                    amp *= -xp.conj(reflection[:, idx])
+                previous = idx + 1
+            amp *= travel[previous][0]
             profile[:, b] += amp
 
     return profile