diff --git a/docs/architecture.md b/docs/architecture.md index 959afd0..e4facfe 100644 --- a/docs/architecture.md +++ b/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`. diff --git a/docs/config_reference.md b/docs/config_reference.md index ae94e2c..299b195 100644 --- a/docs/config_reference.md +++ b/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 | diff --git a/docs/theory.md b/docs/theory.md index 38c313d..ba02d64 100644 --- a/docs/theory.md +++ b/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`): diff --git a/src/pyofdr/core/acquisition.py b/src/pyofdr/core/acquisition.py index 158ee9a..ba91ba4 100644 --- a/src/pyofdr/core/acquisition.py +++ b/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] diff --git a/src/pyofdr/core/config_models.py b/src/pyofdr/core/config_models.py index 9f8de84..bf1a22e 100644 --- a/src/pyofdr/core/config_models.py +++ b/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. diff --git a/src/pyofdr/fiber/attenuation.py b/src/pyofdr/fiber/attenuation.py index 78c10e6..535d0ec 100644 --- a/src/pyofdr/fiber/attenuation.py +++ b/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 diff --git a/src/pyofdr/fiber/multiple_scattering.py b/src/pyofdr/fiber/multiple_scattering.py index f3fd62c..8ddef82 100644 --- a/src/pyofdr/fiber/multiple_scattering.py +++ b/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 diff --git a/src/pyofdr/fiber/profile.py b/src/pyofdr/fiber/profile.py index 79eb1d9..5545dac 100644 --- a/src/pyofdr/fiber/profile.py +++ b/src/pyofdr/fiber/profile.py @@ -7,11 +7,13 @@ from pyofdr.core.acquisition import Acquisition from pyofdr.core.pipeline import PipelineStep -from pyofdr.fiber.attenuation import round_trip_attenuation, round_trip_attenuation_varying +from pyofdr.fiber.attenuation import ( + cumulative_loss_dB, round_trip_attenuation, round_trip_attenuation_varying, +) from pyofdr.fiber.bends import bend_loss_total_dB from pyofdr.fiber.crosstalk import apply_crosstalk from pyofdr.fiber.multiple_scattering import add_ghost_reflections -from pyofdr.fiber.reflectors import apply_connector_losses, inject_reflectors +from pyofdr.fiber.reflectors import apply_connector_losses, inject_reflectors, reflector_bin from pyofdr.utils.constants import C from pyofdr.utils.seeding import derive_seed from pyofdr.utils.units import dB_to_linear, wavelength_range_to_freq_range @@ -34,7 +36,7 @@ def __init__(self, config: dict[str, Any]) -> None: self.attenuation_dB_km = fiber["attenuation_dB_per_km"] self.attenuation_segments = fiber.get("attenuation_segments", []) self.bends = fiber.get("bends", []) - self.reflectors = fiber.get("reflectors", []) + self.reflectors = sorted(fiber.get("reflectors", []), key=lambda r: r["z"]) self.index_segments = fiber.get("index_segments", []) self.index_fluctuations = fiber.get("index_fluctuations", None) self.crosstalk = fiber.get("crosstalk", None) @@ -51,6 +53,7 @@ def __init__(self, config: dict[str, Any]) -> None: self._cached_dz = None self._cached_profile = None self._cached_atten = None + self._cached_ghosts = None def process(self, acq: Acquisition) -> Acquisition: if self._cached_profile is not None: @@ -59,6 +62,7 @@ def process(self, acq: Acquisition) -> Acquisition: acq.dz = self._cached_dz acq.fiber_profile = self._cached_profile acq.attenuation_envelope = self._cached_atten + acq.ghost_profile = self._cached_ghosts return acq xp = self.bk.xp @@ -95,20 +99,20 @@ def process(self, acq: Acquisition) -> Acquisition: profile = xp.stack(parts) # discrete reflectors (connectors, splices, etc) + reflection = None if self.reflectors: - inject_reflectors(profile, z, dz, self.reflectors, xp=xp) - - # cascading multi-bounce ghosts between those reflectors (#35). - # added before the per-bin index phase so each ghost picks up - # the same phi_n at its apparent bin -- a small-perturbation - # approximation (ghost path is longer than 2*z_app, so strictly - # the phase is not phi_n[bin_app]; ok while delta_n stays small). - if self.multiple_scattering is not None and self.reflectors: - add_ghost_reflections( - profile, dz, self.reflectors, - max_order=self.multiple_scattering["max_order"], - xp=xp, - ) + bins = [reflector_bin(r, dz, n_z) for r in self.reflectors] + if len(set(bins)) != len(bins): + raise ValueError("discrete reflectors must occupy distinct spatial bins") + amps = xp.asarray([math.sqrt(r["R"]) for r in self.reflectors]) + reflection = xp.empty((self.n_cores, len(amps)), dtype=xp.complex128) + for c in range(self.n_cores): + rng = self.bk.random_generator( + derive_seed(self.seed, component="reflectors", core=c)) + phase = rng.uniform(0.0, 2.0 * math.pi, len(amps)) + reflection[c] = amps * xp.exp(1j * phase) + inject_reflectors(profile, z, dz, self.reflectors, + reflection=reflection, xp=xp) # small-signal n(z) perturbation -- each scatterer at z picks up # a round-trip phase from the integrated delta_n up to z. Two @@ -142,6 +146,7 @@ def process(self, acq: Acquisition) -> Acquisition: dn_rand[k] = a * dn_rand[k-1] + scale * eta[k] delta_n = dn_rand if delta_n is None else delta_n + dn_rand + phi_n = None if delta_n is not None: k0 = 2.0 * math.pi / self.center_wl phi_n = 2.0 * k0 * xp.cumsum(delta_n) * dz @@ -180,23 +185,46 @@ def process(self, acq: Acquisition) -> Acquisition: mask = (z >= bend["start"]) & (z < bend["end"]) alpha_z = xp.where(mask, alpha_z + extra_dB_km, alpha_z) attenuation = round_trip_attenuation_varying(z, alpha_z, dz, xp=xp) + loss_to_z = cumulative_loss_dB(alpha_z, dz, xp=xp) else: attenuation = round_trip_attenuation(z, self.attenuation_dB_km, xp=xp) + loss_to_z = z * (self.attenuation_dB_km / 1000.0) # connector insertion losses (step-downs in the envelope) if self.reflectors: apply_connector_losses(attenuation, dz, self.reflectors, xp=xp) + ghosts = None + if self.multiple_scattering is not None and self.reflectors: + ghosts = xp.zeros_like(profile) + add_ghost_reflections( + ghosts, dz, self.reflectors, reflection=reflection, + loss_to_z=loss_to_z, + max_order=self.multiple_scattering["max_order"], xp=xp, + ) + # Index/strain/temperature still use the apparent-bin phase. + if phi_n is not None: + ghosts *= xp.exp(1j * phi_n) + if self.crosstalk is not None and self.n_cores > 1: + rng_xt = self.bk.random_generator( + derive_seed(self.seed, component="crosstalk")) + ghosts = apply_crosstalk( + ghosts, z, self.crosstalk["xt_dB_per_km"], + self.crosstalk["topology"], rng_xt, xp=xp, + ) + acq.z = z acq.dz = dz acq.fiber_profile = profile acq.attenuation_envelope = attenuation + acq.ghost_profile = ghosts # cache so subsequent sweeps reuse the same fiber self._cached_z = z self._cached_dz = dz self._cached_profile = profile self._cached_atten = attenuation + self._cached_ghosts = ghosts acq.add_log("fiber", n_z=n_z, dz_mm=dz * 1e3, attenuation_dB_km=self.attenuation_dB_km, diff --git a/src/pyofdr/fiber/reflectors.py b/src/pyofdr/fiber/reflectors.py index ccc1b2f..e448846 100644 --- a/src/pyofdr/fiber/reflectors.py +++ b/src/pyofdr/fiber/reflectors.py @@ -8,19 +8,27 @@ import numpy as np +def reflector_bin(ref, dz, n_z): + idx = int(round(ref["z"] / dz)) + if idx == n_z and ref["z"] <= n_z * dz: + idx = n_z - 1 + return idx + + def inject_reflectors( profile: np.ndarray, z: np.ndarray, dz: float, reflectors: list[dict[str, Any]], *, + reflection=None, xp=np, ) -> np.ndarray: """Add discrete reflectors on top of the Rayleigh backscatter profile. Each reflector is a dict with keys 'z' (position in meters) and - 'R' (power reflectivity, 0-1). The reflector amplitude sqrt(R) - is added to every core at the nearest spatial bin. + 'R' (power reflectivity, 0-1). Optional complex reflection coefficients + have shape (n_cores, n_reflectors); otherwise sqrt(R) is used. Returns the modified profile (in-place). """ @@ -29,18 +37,13 @@ def inject_reflectors( n_z = profile.shape[-1] - for ref in reflectors: - pos = ref["z"] + for i, ref in enumerate(reflectors): R = ref["R"] - idx = int(round(pos / dz)) - # with n_z = ceil(length/dz) the end face can round to n_z -- - # clamp to the last bin instead of silently dropping it (#84) - if idx == n_z and pos <= n_z * dz: - idx = n_z - 1 + idx = reflector_bin(ref, dz, n_z) if idx < 0 or idx >= n_z: continue # amplitude for a point reflector with power reflectivity R - profile[:, idx] += math.sqrt(R) + profile[:, idx] += math.sqrt(R) if reflection is None else reflection[:, i] return profile @@ -52,30 +55,25 @@ def apply_connector_losses( *, xp=np, ) -> np.ndarray: - """Apply insertion-loss steps to the attenuation envelope. + """Apply direct-return transmission and insertion-loss steps. For each reflector with loss_dB > 0 the envelope beyond that point is scaled down by the round-trip amplitude factor: factor = 10^(-loss_dB / 10) (one-way loss_dB applied twice, converted to amplitude). - Reflectors are processed in order of position so that cascaded - connectors accumulate correctly. + The extra loss sits just before the reflecting interface. A return + from that interface sees the loss, but not its own (1-R) transmission. + Returns from later bins cross the interface twice. """ - lossy = [r for r in reflectors if r.get("loss_dB", 0) > 0] - if not lossy: - return attenuation - n_z = len(attenuation) - # sort by position so losses accumulate left-to-right - for ref in sorted(lossy, key=lambda r: r["z"]): - idx = int(round(ref["z"] / dz)) - if idx == n_z and ref["z"] <= n_z * dz: - idx = n_z - 1 # same end-face clamp as inject_reflectors + for ref in reflectors: + idx = reflector_bin(ref, dz, n_z) if idx < 0 or idx >= n_z: continue # round-trip amplitude factor for one-way loss_dB - factor = 10.0 ** (-ref["loss_dB"] / 10.0) + factor = 10.0 ** (-ref.get("loss_dB", 0.0) / 10.0) attenuation[idx:] *= factor + attenuation[idx + 1:] *= 1.0 - ref["R"] return attenuation diff --git a/src/pyofdr/fiber/strain.py b/src/pyofdr/fiber/strain.py index 40ec747..8c9312c 100644 --- a/src/pyofdr/fiber/strain.py +++ b/src/pyofdr/fiber/strain.py @@ -57,6 +57,8 @@ def __init__(self, config: dict[str, Any]) -> None: self._cached_profile = None self._cached_strain = None self._cached_base_profile = None # base (unstrained) profile reference + self._cached_base_ghosts = None + self._cached_ghosts = None # sweep-rate sampling sanity check. motion-specific thresholds # (Nyquist for harmonic, tau/width > 2*T_sweep for thermal/impulsive) @@ -73,6 +75,7 @@ def process(self, acq: Acquisition) -> Acquisition: # static case: cache strained profile, reuse it across sweeps if not self._has_motion and self._cached_profile is not None: acq.fiber_profile = self._cached_profile + acq.ghost_profile = self._cached_ghosts acq.strain_field = self._cached_strain return acq @@ -87,6 +90,7 @@ def process(self, acq: Acquisition) -> Acquisition: # the phase to the fresh profile, not to an already-strained one. if self._cached_base_profile is None: self._cached_base_profile = acq.fiber_profile + self._cached_base_ghosts = acq.ghost_profile base_profile = self._cached_base_profile # freeze motion at the sweep's lab time @@ -104,9 +108,12 @@ def process(self, acq: Acquisition) -> Acquisition: acq.strain_field = eps_fiber # uniform axial strain affects every core the same way acq.fiber_profile = base_profile * xp.exp(1j * phase) + if self._cached_base_ghosts is not None: + acq.ghost_profile = self._cached_base_ghosts * xp.exp(1j * phase) if not self._has_motion: self._cached_profile = acq.fiber_profile + self._cached_ghosts = acq.ghost_profile self._cached_strain = eps_fiber acq.add_log("strain", diff --git a/src/pyofdr/fiber/temperature.py b/src/pyofdr/fiber/temperature.py index 1f7ca84..532c939 100644 --- a/src/pyofdr/fiber/temperature.py +++ b/src/pyofdr/fiber/temperature.py @@ -72,6 +72,8 @@ def process(self, acq: Acquisition) -> Acquisition: # incoming profile so a dynamic upstream step still composes if not self._has_motion and self._cached_phase is not None: acq.fiber_profile = acq.fiber_profile * xp.exp(1j * self._cached_phase) + if acq.ghost_profile is not None: + acq.ghost_profile = acq.ghost_profile * xp.exp(1j * self._cached_phase) acq.temperature_field = self._cached_dT return acq @@ -93,6 +95,8 @@ def process(self, acq: Acquisition) -> Acquisition: acq.temperature_field = dT acq.fiber_profile = acq.fiber_profile * xp.exp(1j * phase) + if acq.ghost_profile is not None: + acq.ghost_profile = acq.ghost_profile * xp.exp(1j * phase) if not self._has_motion: self._cached_phase = phase diff --git a/src/pyofdr/optics/mach_zehnder.py b/src/pyofdr/optics/mach_zehnder.py index 5dcca8d..ac14eee 100644 --- a/src/pyofdr/optics/mach_zehnder.py +++ b/src/pyofdr/optics/mach_zehnder.py @@ -90,6 +90,8 @@ def process(self, acq: Acquisition) -> Acquisition: weighted = acq.fiber_profile if acq.attenuation_envelope is not None: weighted = weighted * acq.attenuation_envelope + if acq.ghost_profile is not None: + weighted = weighted + acq.ghost_profile # laser coherence roll-off: for a Lorentzian-linewidth laser each # scatterer at round-trip delay tau = 2*n*z/c sees its beat amplitude diff --git a/src/pyofdr/utils/seeding.py b/src/pyofdr/utils/seeding.py index 1652223..cf541b7 100644 --- a/src/pyofdr/utils/seeding.py +++ b/src/pyofdr/utils/seeding.py @@ -16,6 +16,7 @@ "adc": 3, "crosstalk": 4, "index_fluctuations": 5, + "reflectors": 6, } diff --git a/tests/test_fiber.py b/tests/test_fiber.py index e44ba4c..7796619 100644 --- a/tests/test_fiber.py +++ b/tests/test_fiber.py @@ -213,8 +213,7 @@ def test_default_disabled(self): np.testing.assert_array_equal(a.fiber_profile, b.fiber_profile) def test_double_bounce_ghost_position_and_amplitude(self): - """Two reflectors at z=0 and z=z_b -> ghost at 2*z_b with - amplitude sqrt(R_a) * R_b.""" + """The first ghost crosses the input reflector twice.""" z_b = 0.3 R_a, R_b = 0.01, 0.04 cfg = {**CFG, "fiber": {**CFG["fiber"], @@ -227,8 +226,8 @@ def test_double_bounce_ghost_position_and_amplitude(self): # we end up off-by-one when 2*round(x) != round(2*x)) bin_b = int(round(z_b / acq.dz)) ghost_bin = 2 * bin_b - amp = np.abs(acq.fiber_profile[0, ghost_bin]) - np.testing.assert_allclose(amp, np.sqrt(R_a) * R_b, atol=1e-10) + amp = np.abs(acq.ghost_profile[0, ghost_bin]) + np.testing.assert_allclose(amp, (1 - R_a) * np.sqrt(R_a) * R_b, atol=1e-10) def test_no_ghost_when_disabled(self): z_b = 0.3 @@ -261,8 +260,7 @@ def test_out_of_range_reflector_rejected(self): FiberGenerator(cfg) def test_three_reflector_amplitudes(self): - """End-to-end physics check: three known reflectors, verify all - five distinct order-2 ghost amplitudes against prod sqrt(R).""" + """Three reflectors give four bins, including two coincident paths.""" # carefully spaced so each ghost lands at its own distinct bin z_a, z_b, z_c = 0.0, 0.07, 0.19 R_a, R_b, R_c = 0.01, 0.04, 0.09 @@ -276,17 +274,21 @@ def test_three_reflector_amplitudes(self): dz = acq.dz Bb = int(round(z_b / dz)) Bc = int(round(z_c / dz)) - prof = np.abs(acq.fiber_profile[0]) + prof = np.abs(acq.ghost_profile[0]) # (b,a,b) -> 2*Bb, amp sqrt(R_a)*R_b - np.testing.assert_allclose(prof[2 * Bb], np.sqrt(R_a) * R_b, atol=1e-10) + np.testing.assert_allclose(prof[2 * Bb], + (1 - R_a) * np.sqrt(R_a) * R_b, atol=1e-10) # (c,a,c) -> 2*Bc, amp sqrt(R_a)*R_c - np.testing.assert_allclose(prof[2 * Bc], np.sqrt(R_a) * R_c, atol=1e-10) + np.testing.assert_allclose(prof[2 * Bc], + (1 - R_a) * (1 - R_b)**2 * np.sqrt(R_a) * R_c, + atol=1e-10) # (c,b,c) -> 2*Bc - Bb, amp sqrt(R_b)*R_c np.testing.assert_allclose(prof[2 * Bc - Bb], - np.sqrt(R_b) * R_c, atol=1e-10) + (1 - R_a) * (1 - R_b) * np.sqrt(R_b) * R_c, + atol=1e-10) # (b,a,c) and (c,a,b) both land at Bb+Bc and add coherently: # 2 * sqrt(R_a) * sqrt(R_b) * sqrt(R_c) - expected = 2.0 * np.sqrt(R_a * R_b * R_c) + expected = 2.0 * (1 - R_a) * (1 - R_b) * np.sqrt(R_a * R_b * R_c) np.testing.assert_allclose(prof[Bb + Bc], expected, atol=1e-10) def test_higher_order_adds_extra_paths(self): @@ -302,7 +304,7 @@ def test_higher_order_adds_extra_paths(self): "multiple_scattering": {"max_order": 3}}} a = FiberGenerator(cfg2).process(Acquisition()) b = FiberGenerator(cfg3).process(Acquisition()) - diff = np.abs(np.abs(b.fiber_profile) - np.abs(a.fiber_profile)).max() + diff = np.abs(b.ghost_profile - a.ghost_profile).max() assert diff > 0 @@ -550,10 +552,9 @@ def test_zero_loss_segment_matches_baseline(self): a = self._run(segs, base=0.5) b = FiberGenerator({**CFG, "fiber": {**CFG["fiber"], "attenuation_dB_per_km": 0.5}}).process(Acquisition()) - # numerical: varying path uses cumsum*dz, constant uses alpha*z, - # off by one bin -- allow 1e-6 tolerance + # A constant segment must integrate just like the scalar loss. np.testing.assert_allclose(a.attenuation_envelope, - b.attenuation_envelope, rtol=1e-4) + b.attenuation_envelope, rtol=1e-12) class TestIndexPerturbation: diff --git a/tests/test_reflector_paths.py b/tests/test_reflector_paths.py new file mode 100644 index 0000000..9868cce --- /dev/null +++ b/tests/test_reflector_paths.py @@ -0,0 +1,240 @@ +"""Small reflector networks with known complex returns.""" + +import numpy as np +import pytest + +from helpers import CFG +from pyofdr.core.acquisition import Acquisition +from pyofdr.fiber.attenuation import cumulative_loss_dB, round_trip_attenuation_varying +from pyofdr.fiber.crosstalk import apply_crosstalk +from pyofdr.fiber.multiple_scattering import add_ghost_reflections +from pyofdr.fiber.profile import FiberGenerator +from pyofdr.fiber.strain import StrainPerturbation +from pyofdr.fiber.temperature import TemperaturePerturbation +from pyofdr.optics.mach_zehnder import MachZehnder +from pyofdr.source.swept_laser import SweptLaser +from pyofdr.utils.seeding import derive_seed + + +def fiber_config(refs, **fiber_options): + return {**CFG, "fiber": {**CFG["fiber"], "reflectors": refs, + "rayleigh_coefficient_dB": -300, + "multiple_scattering": {"max_order": 3}, **fiber_options}} + + +def ghosts(refs, reflection, loss=None, order=2): + out = np.zeros((1, 100), complex) + add_ghost_reflections(out, 1.0, refs, reflection=np.array([reflection]), + loss_to_z=loss, max_order=order) + return out[0] + + +def test_two_reflector_complex_amplitudes(): + a, b = 5, 20 + ra, rb = 0.2 * np.exp(0.7j), 0.3 * np.exp(-1.1j) + refs = [{"z": a, "R": 0.04}, {"z": b, "R": 0.09}] + out = ghosts(refs, [ra, rb], order=3) + expected = np.zeros(100, complex) + expected[2*b-a] = -(1 - 0.04) * ra.conjugate() * rb**2 + expected[3*b-2*a] = (1 - 0.04) * ra.conjugate()**2 * rb**3 + np.testing.assert_allclose(out, expected, atol=1e-16) + + +@pytest.mark.parametrize("loss_bin, passes", [(2, 2), (12, 4), (25, 0)]) +def test_localized_loss_counts_real_crossings(loss_bin, passes): + refs = [{"z": 5, "R": 0.04}, {"z": 20, "R": 0.09}] + # A 3 dB step before, inside, or beyond the cavity. + loss = np.where(np.arange(100) > loss_bin, 3.0, 0.0) + actual = ghosts(refs, [0.2, 0.3], loss=loss)[35] + expected = -0.96 * 0.2 * 0.3**2 * 10**(-passes * 3 / 20) + np.testing.assert_allclose(actual, expected, atol=1e-16) + + +def test_uniform_loss_matches_apparent_distance(): + refs = [{"z": 5, "R": 0.04}, {"z": 20, "R": 0.09}] + out = ghosts(refs, [0.2, 0.3], loss=0.1 * np.arange(100)) + np.testing.assert_allclose(out[35], -0.96 * 0.2 * 0.3**2 * 10**(-7 / 20)) + + +@pytest.mark.parametrize("where, R, loss, passes", [ + (0, 0.16, 1.0, 2), (12, 0.16, 1.0, 4), (25, 1.0, 4000.0, 0), + (12, 0.0, 3.0, 4), (12, 1.0, 0.0, 4), +]) +def test_intermediate_interfaces(where, R, loss, passes): + refs = [{"z": 5, "R": 0.04}, {"z": 20, "R": 0.09}, + {"z": where, "R": R, "loss_dB": loss}] + out = ghosts(refs, [0.2, 0.3, np.sqrt(R)]) + t = np.sqrt(1 - R) * 10**(-loss / 20) + expected = -0.96 * 0.2 * 0.3**2 * t**passes + np.testing.assert_allclose(out[35], expected, atol=1e-16) + assert np.all(np.isfinite(out)) + + +def test_loss_at_bouncing_interfaces(): + refs = [{"z": 5, "R": 0.04, "loss_dB": 1.0}, + {"z": 20, "R": 0.09, "loss_dB": 2.0}] + # Twice through the loss before a, four times through that before b. + expected = -0.96 * 0.2 * 0.3**2 * 10**(-(2*1 + 4*2) / 20) + np.testing.assert_allclose(ghosts(refs, [0.2, 0.3])[35], expected) + + +@pytest.mark.parametrize("R_a, R_b", [(0.0, 0.09), (1.0, 0.09), (0.04, 0.0)]) +def test_no_accessible_cavity_gives_no_ghosts(R_a, R_b): + refs = [{"z": 0, "R": R_a}, {"z": 20, "R": R_b}] + out = ghosts(refs, [np.sqrt(R_a), np.sqrt(R_b)], order=3) + np.testing.assert_array_equal(out, 0) + + +def test_seeded_phases_and_direct_transmission(): + refs = [{"z": 0.1, "R": 0.04}, {"z": 0.3, "R": 0.09}] + cfg = fiber_config(refs, n_cores=2) + generator = FiberGenerator(cfg) + first = generator.process(Acquisition()) + again = generator.process(Acquisition(sweep_index=1)) + repeat = FiberGenerator(cfg).process(Acquisition()) + changed = FiberGenerator({**cfg, "simulation": {"seed": 43}}).process(Acquisition()) + for acq in (again, repeat): + np.testing.assert_array_equal(first.fiber_profile, acq.fiber_profile) + np.testing.assert_array_equal(first.ghost_profile, acq.ghost_profile) + a, b = [round(r["z"] / first.dz) for r in refs] + assert not np.allclose(first.fiber_profile[:, [a, b]].imag, 0) + assert not np.allclose(first.fiber_profile[0, [a, b]], first.fiber_profile[1, [a, b]]) + assert not np.allclose(first.fiber_profile[:, [a, b]], changed.fiber_profile[:, [a, b]]) + np.testing.assert_allclose(np.abs(first.fiber_profile[:, b]) * first.attenuation_envelope[b], + 0.96 * 0.3) + ra, rb = first.fiber_profile[:, a], first.fiber_profile[:, b] + np.testing.assert_allclose(first.ghost_profile[:, 2*b-a], -0.96 * ra.conj() * rb**2) + reversed_cfg = fiber_config(list(reversed(refs)), n_cores=2) + reversed_acq = FiberGenerator(reversed_cfg).process(Acquisition()) + np.testing.assert_array_equal(first.fiber_profile, reversed_acq.fiber_profile) + np.testing.assert_array_equal(first.ghost_profile, reversed_acq.ghost_profile) + bare = FiberGenerator(fiber_config([], n_cores=2)).process(Acquisition()) + mask = np.ones(len(first.z), bool) + mask[[a, b]] = False + np.testing.assert_array_equal(first.fiber_profile[:, mask], bare.fiber_profile[:, mask]) + + +@pytest.mark.parametrize("kind", ["segment", "bend"]) +@pytest.mark.parametrize("start,end,passes", [(0.02, 0.06, 2), (0.15, 0.20, 4), + (0.35, 0.40, 0)]) +def test_generated_ghost_uses_segment_and_bend_loss(kind, start, end, passes): + refs = [{"z": 0.1, "R": 0.04}, {"z": 0.3, "R": 0.09}] + options = {"multiple_scattering": {"max_order": 2}} + if kind == "segment": + options["attenuation_segments"] = [{"start": start, "end": end, + "attenuation_dB_per_km": 1000.0}] + loss_per_m = 1.0 + else: + options["bends"] = [{"start": start, "end": end, "radius": 0.005, + "turns": 1, "A_dB_per_turn": np.e, "R_c": 0.005}] + loss_per_m = 1.0 / (end - start) # 1 dB in total + acq = FiberGenerator(fiber_config(refs, **options)).process(Acquisition()) + a, b = [round(r["z"] / acq.dz) for r in refs] + covered_length = np.count_nonzero((acq.z >= start) & (acq.z < end)) * acq.dz + expected = 0.96 * 0.2 * 0.09 * 10**(-passes * loss_per_m * covered_length / 20) + np.testing.assert_allclose(abs(acq.ghost_profile[0, 2*b-a]), expected, rtol=1e-12) + + +def test_index_phase_and_crosstalk_also_apply_to_ghosts(): + refs = [{"z": 0.0, "R": 0.04}, {"z": 0.2, "R": 0.09}] + cfg = fiber_config(refs, n_cores=2) + bare = FiberGenerator(cfg).process(Acquisition()) + index = {"start": 0.1, "end": 0.5, "delta_n": 1e-7} + xt = {"xt_dB_per_km": -20.0, "topology": "linear"} + cfg["fiber"] = {**cfg["fiber"], "index_segments": [index], "crosstalk": xt} + actual = FiberGenerator(cfg).process(Acquisition()) + dn = np.where((bare.z >= 0.1) & (bare.z < 0.5), 1e-7, 0) + phase = 4 * np.pi / cfg["source"]["center_wavelength"] * np.cumsum(dn) * bare.dz + rng = np.random.default_rng(derive_seed(42, component="crosstalk")) + expected = apply_crosstalk(bare.ghost_profile * np.exp(1j * phase), bare.z, + -20.0, "linear", rng) + np.testing.assert_allclose(actual.ghost_profile, expected, atol=1e-15) + + +def test_full_mirror_blocks_direct_returns_but_not_earlier_cavity(): + refs = [{"z": 0.0, "R": 0.04}, {"z": 0.2, "R": 0.09}, + {"z": 0.317, "R": 1.0}] + cfg = fiber_config(refs, multiple_scattering={"max_order": 2}) + acq = FiberGenerator(cfg).process(Acquisition()) + b = round(0.2 / acq.dz) + mirror = round(0.317 / acq.dz) + np.testing.assert_array_equal(acq.attenuation_envelope[mirror + 1:], 0) + np.testing.assert_allclose(acq.attenuation_envelope[mirror], 0.96 * 0.91) + np.testing.assert_allclose(abs(acq.ghost_profile[0, 2*b]), 0.96 * 0.2 * 0.09) + + +def test_loss_integral_starts_at_input(): + alpha = np.array([1000., 2000., 0., 4000.]) + loss = cumulative_loss_dB(alpha, 1.0) + np.testing.assert_array_equal(loss, [0., 1., 3., 3.]) + envelope = round_trip_attenuation_varying(np.arange(4), alpha, 1.0) + np.testing.assert_allclose(envelope, 10**(-loss / 10)) + + +def test_unresolved_reflectors_are_rejected(): + refs = [{"z": 0.1, "R": 0.04}, {"z": 0.1, "R": 0.09}] + with pytest.raises(ValueError, match="distinct spatial bins"): + FiberGenerator(fiber_config(refs)).process(Acquisition()) + + +def test_frequency_domain_cavity_reference(): + # Sum the infinite reflection series recursively in frequency space. + # Weak mirrors make the omitted orders (>5 bounces) smaller than 1e-10. + bins = [5, 13, 27] + R = np.array([0.0004, 0.0009, 0.0006]) + r = np.sqrt(R) * np.exp(1j * np.array([0.7, -1.1, 2.0])) + extra_dB = np.array([0.2, 0.7, 0.4]) + refs = [{"z": b, "R": power, "loss_dB": loss} + for b, power, loss in zip(bins, R, extra_dB)] + loss = np.cumsum(np.where(np.arange(100) < 18, 0.01, 0.08)) + loss -= loss[0] + theta = 2 * np.pi * np.arange(512) / 512 + response = np.zeros(512, complex) + for i in range(2, -1, -1): + if i < 2: + response *= np.exp(1j * theta * (bins[i+1] - bins[i])) + response *= 10**(-(loss[bins[i+1]] - loss[bins[i]]) / 10) + q2 = 10**(-extra_dB[i] / 10) + response = q2 * (r[i] + (1 - R[i]) * response / (1 + r[i].conjugate() * response)) + response *= np.exp(1j * theta * bins[0]) * 10**(-loss[bins[0]] / 10) + expected = np.fft.fft(response)[:100] / 512 + actual = ghosts(refs, r, loss=loss, order=3) + upstream = 1.0 + for b, ri, power, extra in zip(bins, r, R, extra_dB): + actual[b] += upstream * ri * 10**(-(loss[b] + extra) / 10) + upstream *= (1 - power) * 10**(-extra / 10) + np.testing.assert_allclose(actual, expected, rtol=0, atol=1e-10) + + +@pytest.mark.parametrize("dynamic", [False, True]) +def test_perturbations_preserve_ghosts_across_sweeps(dynamic): + cfg = fiber_config([{"z": 0.0, "R": 0.04}, {"z": 0.2, "R": 0.09}]) + segment = {"start": 0.1, "end": 0.5, "epsilon": 1e-6} + if dynamic: + segment["motion"] = {"kind": "harmonic", "amplitude": 1e-6, "frequency": 3.0} + cfg["strain"] = {"segments": [segment]} + cfg["temperature"] = {"segments": [{"start": 0.1, "end": 0.5, "delta_T": 0.2}]} + fiber, strain, temp = FiberGenerator(cfg), StrainPerturbation(cfg), TemperaturePerturbation(cfg) + base = fiber.process(Acquisition()) + ghost_bin = 2 * round(0.2 / base.dz) + for sweep in range(3): + acq = temp.process(strain.process(fiber.process(Acquisition(sweep_index=sweep)))) + phase = acq.fiber_profile[0, ghost_bin] / base.fiber_profile[0, ghost_bin] + np.testing.assert_allclose(acq.ghost_profile[0, ghost_bin], + base.ghost_profile[0, ghost_bin] * phase) + + +def test_optics_adds_ghosts_without_another_attenuation(): + cfg = fiber_config([{"z": 0.0, "R": 0.04}, {"z": 0.2, "R": 0.09}, + {"z": 0.3, "R": 1.0}]) + cfg["source"] = {**cfg["source"], "sweep_duration": 1e-4, + "sweep_range": 1e-9, "linewidth": 0.0} + cfg["adc"] = {**cfg["adc"], "sample_rate": 200e6} + cfg["optics"] = {**cfg["optics"], "circulator": { + **cfg["optics"]["circulator"], "insertion_loss_dB": 0.0}} + acq = SweptLaser(cfg).process(FiberGenerator(cfg).process(Acquisition())) + weighted = acq.fiber_profile * acq.attenuation_envelope + acq.ghost_profile + expected = np.fft.ifft(weighted, n=acq.n_samples, axis=-1) * acq.n_samples + expected = np.abs(acq.E_source)**2 * expected.real + actual = MachZehnder(cfg).process(acq).photocurrent_main + np.testing.assert_allclose(actual, expected, atol=1e-16) diff --git a/tests/test_sampling.py b/tests/test_sampling.py index df8886b..ada90ce 100644 --- a/tests/test_sampling.py +++ b/tests/test_sampling.py @@ -67,7 +67,8 @@ def test_reflector_matches_physical_chirp(count, bin_index): position = cfg["fiber"]["reflectors"][0]["z"] frequency = 2 * cfg["fiber"]["n_core"] * position * span / (C * T) amplitude = 2 * 0.01 * np.sqrt(0.4 * 0.6 * 1e-4) - expected = amplitude * np.cos(2 * np.pi * frequency * acq.t) + phase = np.angle(acq.fiber_profile[0, bin_index]) + expected = amplitude * np.cos(2 * np.pi * frequency * acq.t + phase) np.testing.assert_allclose(acq.photocurrent_main[0], expected, rtol=0, atol=amplitude * 1e-9) np.testing.assert_allclose(acq.t, np.arange(acq.n_samples) / fs, rtol=5e-16, atol=0) assert acq.t[-1] < T @@ -86,9 +87,10 @@ def test_complex_multicore_profile_matches_direct_sum(count): profile = 0.01 * (rng.normal(size=acq.fiber_profile.shape) + 1j * rng.normal(size=acq.fiber_profile.shape)) acq.fiber_profile = profile.copy() + weighted = profile * acq.attenuation_envelope expected = np.zeros((2, acq.n_samples), dtype=complex) for k in range(profile.shape[1]): - expected += profile[:, k, None] * np.exp(2j * np.pi * k * acq.t / 1e-4) + expected += weighted[:, k, None] * np.exp(2j * np.pi * k * acq.t / 1e-4) prefactor = 2 * np.sqrt(0.4 * 0.6) * np.abs(acq.E_source)**2 expected = prefactor * expected.real acq = optics.process(acq) @@ -134,11 +136,13 @@ def test_nonlinear_warp_uses_physical_time_grid(a2): s = acq.t + a2 * acq.t**2 / (span / T) amplitude = 2 * 0.01 * np.sqrt(0.4 * 0.6 * 1e-4) # Separate the sampling correction from the existing linear interpolation error. - expected_interpolated = np.interp(s, acq.t, amplitude * np.cos(2 * np.pi * 16 * acq.t / T)) + phase = np.angle(acq.fiber_profile[0, 16]) + expected_interpolated = np.interp( + s, acq.t, amplitude * np.cos(2 * np.pi * 16 * acq.t / T + phase)) np.testing.assert_allclose(acq.photocurrent_main[0], expected_interpolated, rtol=0, atol=amplitude * 1e-9) inside = (s >= acq.t[0]) & (s <= acq.t[-1]) - exact = amplitude * np.cos(2 * np.pi * 16 * s[inside] / T) + exact = amplitude * np.cos(2 * np.pi * 16 * s[inside] / T + phase) rms = np.sqrt(np.mean((acq.photocurrent_main[0, inside] - exact)**2)) / amplitude assert rms < 1e-4 diff --git a/tests/test_seeding.py b/tests/test_seeding.py index 37be012..4b0ab3c 100644 --- a/tests/test_seeding.py +++ b/tests/test_seeding.py @@ -32,7 +32,7 @@ def test_neighboring_campaign_seeds_do_not_reuse_sweeps(): def test_long_campaign_seeds_are_distinct(): components = ["fiber", "laser", "detector", "adc", "crosstalk", - "index_fluctuations"] + "index_fluctuations", "reflectors"] seeds = set() for component in components: for core in range(2):