Add chromatic dispersion with sign-verified validation

Fiber gains a `dispersion` parameter (D in ps/nm/km). Propagation solves the
linear NLSE exactly in the frequency domain: A(z,w) = A(0,w)*exp(i*b2*w^2*z/2),
an all-pass phase rotation, so it conserves energy and inverts exactly.

beta2 is computed per band from that band's own wavelength. Two channels
hundreds of nanometres apart really do disperse differently, and the multi-band
signal model gets that right with no special handling — a single-carrier model
could not express it at all.

New: kernels.py (frequency grid, D->beta2, dispersion propagator) as the narrow
array-in/array-out boundary a CuPy or native back-end can later replace;
analysis.py (RMS time width, peak/instantaneous power); GaussianPulse source,
which exists because a Gaussian through pure GVD has a closed-form solution and
is therefore the reference input for validating the fiber.

Validation (18 new tests, 85 total):
- Unchirped broadening T1/T0 = sqrt(1 + (z/L_D)^2) across z/L_D = 0.5 to 3
- L_D is where the pulse broadens by sqrt(2); broadening scales as T0^2
- Chirped broadening, including compression to ~45% for C=+2 in anomalous fiber
- Energy conservation (Parseval), exact reversibility with +D then -D
- beta2 = -21.7 ps^2/km for D = 17 ps/nm/km at 1550 nm, with the sign asserted
- Each band dispersed at its own wavelength

The chirped tests are not redundant with the unchirped ones: the unchirped
broadening factor is even in beta2 and passes under either sign. Verified by
flipping the sign in the propagator — 13 unchirped tests still passed while all
5 chirped ones failed. That blind spot is documented in kernels.py.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

Author: ehsun-sh

README.md

@@ -12,11 +12,11 @@
 > ### ⚠️ Project status: pre-alpha, Phase 0
 >
 > The engine core exists and is tested: simulation context, multi-band optical signal model,
-> typed ports, component/parameter system, and the block-mode scheduler. Four components are
-> implemented — CW laser, fiber (attenuation only), combiner/attenuator, power meter — and the
-> physics is validated against closed-form results in CI.
+> typed ports, component/parameter system, and the block-mode scheduler. Implemented components:
+> CW laser, Gaussian pulse source, fiber (attenuation + chromatic dispersion), combiner,
+> attenuator, power meter — with the physics validated against closed-form results in CI.
 >
-> **Not implemented yet:** modulators, photodetectors, dispersion, SSFM/nonlinearity, DSP,
+> **Not implemented yet:** modulators, photodetectors, SSFM/nonlinearity, PMD, amplifiers, DSP,
 > BER/eye analysis, and the GUI. Those are Phases 1–2; see the [roadmap](#roadmap).
 >
 > This is not yet a useful simulator. It is a foundation with the expensive decisions made and
@@ -222,8 +222,11 @@ Every physics block ships with a test against a closed-form result, run in CI
 | Source power | Independent of the simulated time window | ✅ |
 | Phase noise | Broadens the line, conserves average power | ✅ |
 | Multi-carrier | Channels stay separate bands; spacing does not drive `Fs` | ✅ |
-| Gaussian pulse, CD only | `T(z) = T₀·√(1 + (z/L_D)²)`, `L_D = T₀²/\|β₂\|` | ⬜ |
-| Lossless SSFM | Energy conserved (Parseval) | ⬜ |
+| Gaussian pulse, CD only | `T₁/T₀ = √(1 + (z/L_D)²)`, `L_D = T₀²/\|β₂\|` | ✅ |
+| Chirped Gaussian | `T₁/T₀ = √((1 + Cβ₂z/T₀²)² + (β₂z/T₀²)²)` — pins the sign of β₂ | ✅ |
+| Dispersion compensation | `+D` then `−D` restores the input sample-for-sample | ✅ |
+| GVD | Energy conserved (Parseval); β₂ = −Dλ²/2πc per band | ✅ |
+| Lossless SSFM | Energy conserved with nonlinearity | ⬜ |
 | Fundamental soliton (N=1) | Envelope magnitude invariant along propagation | ⬜ |
 | Ideal push-pull MZM | `P_out/P_in = cos²(πV / 2V_π)` | ⬜ |
 | PIN detector | `I = R·P`; shot `σ² = 2qIB`; thermal `σ² = 4kTB/R_L` | ⬜ |

src/oosim/analysis.py

@@ -0,0 +1,50 @@
+"""Measurements computed from signals.
+
+These are plain functions rather than components on purpose: the same reductions
+are needed by the test suite, by scripts, and later by measurement blocks and the
+GUI. Keeping them here means one implementation and one set of tests.
+
+They are also where result data gets *reduced*. A browser must never receive a
+million raw samples; it receives what these functions return.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from .signals import Band
+
+
+def instantaneous_power(band: Band) -> np.ndarray:
+    """Instantaneous power [W] per sample, summed over both polarizations."""
+    return (np.abs(band.Ex.astype(np.complex128)) ** 2) + (
+        np.abs(band.Ey.astype(np.complex128)) ** 2
+    )
+
+
+def rms_time_width(band: Band) -> float:
+    """RMS width of the intensity envelope in time [s].
+
+    The second central moment of ``|A(t)|**2``. For a Gaussian intensity profile
+    ``exp(-(T/T0)**2)`` this equals ``T0 / sqrt(2)``, so *ratios* of this quantity
+    track ``T1/T0`` exactly — which is what the dispersion validation compares
+    against the analytical broadening factor.
+
+    The time window is periodic, so this is only meaningful while the pulse stays
+    well inside it. A pulse that has spread far enough to wrap around will report
+    a width that is wrong rather than merely imprecise.
+    """
+    power = instantaneous_power(band)
+    total = float(power.sum())
+    if total <= 0.0:
+        raise ValueError("cannot measure the width of a signal carrying no power")
+
+    t = np.arange(band.num_samples, dtype=np.float64) / band.fs
+    mean = float((power * t).sum() / total)
+    variance = float((power * (t - mean) ** 2).sum() / total)
+    return float(np.sqrt(variance))
+
+
+def peak_power(band: Band) -> float:
+    """Highest instantaneous power in the window [W]."""
+    return float(instantaneous_power(band).max())

src/oosim/components/__init__.py

@@ -10,6 +10,6 @@
 from .fiber import Fiber
 from .meters import PowerMeter
 from .passive import Attenuator, Combiner
-from .sources import CWLaser
+from .sources import CWLaser, GaussianPulse
 
-__all__ = ["Attenuator", "CWLaser", "Combiner", "Fiber", "PowerMeter"]
+__all__ = ["Attenuator", "CWLaser", "Combiner", "Fiber", "GaussianPulse", "PowerMeter"]

src/oosim/components/fiber.py

@@ -1,32 +1,45 @@
 """Optical fiber.
 
-Currently the linear, lossy model only. Chromatic dispersion is the next slice;
-Kerr nonlinearity and PMD follow in Phase 1.5 with an adaptive-step SSFM.
+Currently the linear model: attenuation and chromatic dispersion. Kerr
+nonlinearity and PMD follow in Phase 1.5, with an adaptive-step SSFM built on
+the same frequency-domain kernel this module already uses.
+
+Model references: G. P. Agrawal, *Nonlinear Fiber Optics*, ch. 2-3
+(NLSE, GVD-induced pulse broadening); ITU-T G.652 for typical parameter values.
 """
 
 from __future__ import annotations
 
 from ..component import Component, Param, PortType
 from ..context import SimulationContext
-from ..signals import OpticalSignal, Signal
+from ..kernels import dispersion_to_beta2, propagate_dispersion
+from ..signals import Band, OpticalSignal, Signal
 from ..units import db_to_linear
 
 
 class Fiber(Component):
-    """Single-mode fiber, attenuation only.
+    """Single-mode fiber: attenuation and chromatic dispersion.
+
+    Attenuation follows ``P_out = P_in * 10**(-alpha_dB_per_km * L_km / 10)``, so
+    the field amplitude is scaled by the square root of that factor.
+
+    Dispersion is applied per band, using each band's *own* centre wavelength to
+    compute β₂. Two channels a few nanometres apart really do see different
+    dispersion, and because every band carries its own centre frequency the model
+    gets that right for free — a single-carrier signal model could not express it.
 
-    Power obeys ``P_out = P_in * 10**(-alpha_dB_per_km * L_km / 10)``, so the
-    field amplitude is scaled by the square root of that factor. Attenuation is
-    wavelength-independent in this model: every band and every noise bin is
-    scaled identically. Making loss wavelength-dependent is a later refinement
-    and does not change the interface.
+    Attenuation is wavelength-independent here, and the dispersion slope is not
+    modelled (β₃ = 0); both are refinements that do not change the interface.
     """
 
     display_name = "Optical Fiber"
     category = "Fiber"
 
     length = Param(80.0, unit="km", min=0.0, doc="Fiber span length")
     attenuation = Param(0.2, unit="dB/km", min=0.0, doc="Attenuation coefficient")
+    dispersion = Param(
+        0.0, unit="ps/nm/km", doc="Dispersion parameter D at the band wavelength (0 disables)"
+    )
 
     inputs = {"in": PortType.OPTICAL}
     outputs = {"out": PortType.OPTICAL}
@@ -36,14 +49,33 @@ def loss_db(self) -> float:
         # si() gives dB/m and metres, so their product is dB.
         return self.si("attenuation") * self.si("length")
 
+    def beta2_at(self, wavelength: float) -> float:
+        """Group-velocity dispersion β₂ [s²/m] at ``wavelength`` [m]."""
+        return dispersion_to_beta2(self.si("dispersion"), wavelength)
+
     def run(self, ctx: SimulationContext, inputs: dict[str, Signal]) -> dict[str, Signal]:
         signal: OpticalSignal = inputs["in"]
+        distance = self.si("length")
         power_factor = db_to_linear(-self.loss_db())
         amplitude_factor = power_factor**0.5
 
+        bands = []
+        for band in signal.bands:
+            beta2 = self.beta2_at(band.wavelength)
+            Ex = propagate_dispersion(band.Ex, band.fs, beta2, distance) * amplitude_factor
+            Ey = propagate_dispersion(band.Ey, band.fs, beta2, distance) * amplitude_factor
+            bands.append(
+                Band(
+                    Ex=Ex.astype(ctx.complex_dtype),
+                    Ey=Ey.astype(ctx.complex_dtype),
+                    f0=band.f0,
+                    fs=band.fs,
+                )
+            )
+
         return {
             "out": OpticalSignal(
-                bands=tuple(b.scale_amplitude(amplitude_factor) for b in signal.bands),
+                bands=tuple(bands),
                 noise=tuple(n.scale_power(power_factor) for n in signal.noise),
             )
         }

src/oosim/components/sources.py

@@ -54,3 +54,45 @@ def run(self, ctx: SimulationContext, inputs: dict[str, Signal]) -> dict[str, Si
             fs=ctx.sample_rate,
         )
         return {"out": OpticalSignal(bands=(band,))}
+
+
+class GaussianPulse(Component):
+    """A single chirped Gaussian pulse, centred in the time window.
+
+    Follows the standard form (Agrawal, *Nonlinear Fiber Optics*, eq. 3.2.1)::
+
+        A(0, T) = sqrt(P0) * exp(-(1 + i*C) / 2 * (T / T0)**2)
+
+    so the intensity envelope is ``P0 * exp(-(T/T0)**2)`` and ``width`` is T0, the
+    1/e half-width of the *intensity* (``T_FWHM = 1.665 * T0``).
+
+    This exists because dispersion has an exact analytical solution for a Gaussian
+    input, which makes it the reference input for validating the fiber model.
+    """
+
+    display_name = "Gaussian Pulse"
+    category = "Optical Sources"
+
+    peak_power = Param(0.0, unit="dBm", doc="Peak power P0 (not average power)")
+    width = Param(10.0, unit="ps", min=0.0, doc="T0, the 1/e intensity half-width")
+    chirp = Param(0.0, doc="Linear chirp parameter C; sign matters against beta2")
+    wavelength = Param(1550.0, unit="nm", min=1200.0, max=1700.0, doc="Vacuum wavelength")
+
+    outputs = {"out": PortType.OPTICAL}
+
+    def run(self, ctx: SimulationContext, inputs: dict[str, Signal]) -> dict[str, Signal]:
+        t0 = self.si("width")
+        if t0 <= 0.0:
+            raise ValueError(f"{self.label}: width must be positive, got {self.width}")
+
+        tau = (ctx.time_axis() - ctx.time_window / 2.0) / t0
+        amplitude = np.sqrt(self.si("peak_power"))
+        Ex = amplitude * np.exp(-(1.0 + 1j * self.chirp) * tau**2 / 2.0)
+
+        band = Band(
+            Ex=Ex.astype(ctx.complex_dtype),
+            Ey=np.zeros(ctx.num_samples, dtype=ctx.complex_dtype),
+            f0=C_LIGHT / self.si("wavelength"),
+            fs=ctx.sample_rate,
+        )
+        return {"out": OpticalSignal(bands=(band,))}

src/oosim/kernels.py

@@ -0,0 +1,73 @@
+"""Numerical kernels.
+
+Everything computationally heavy goes through this module. It is deliberately
+narrow — a handful of array-in/array-out functions with no knowledge of
+components, graphs, or units — so that the back-end can change without touching
+any physics code above it. Today that back-end is NumPy; CuPy (`cupy.fft` is a
+drop-in for `numpy.fft`) and a native module are the intended next options.
+
+**FFT library.** `numpy.fft` uses pocketfft (BSD). FFTW is *not* used and must
+not be introduced: it is GPL-2.0-or-later, and linking it — directly or through
+`pyFFTW` — would relicense the whole project.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from .units import C_LIGHT
+
+
+def angular_frequency_grid(num_samples: int, sample_rate: float) -> np.ndarray:
+    """Angular frequency offsets from the band centre [rad/s], in FFT order.
+
+    Returned in `numpy.fft` output order (positive frequencies first, then
+    negative), so it multiplies an un-shifted spectrum directly.
+    """
+    return 2.0 * np.pi * np.fft.fftfreq(num_samples, d=1.0 / sample_rate)
+
+
+def dispersion_to_beta2(dispersion: float, wavelength: float) -> float:
+    """Convert the dispersion parameter D [s/m²] to the GVD parameter β₂ [s²/m].
+
+    ``beta2 = -D * lambda**2 / (2*pi*c)``
+
+    The sign matters: standard single-mode fiber has D > 0 at 1550 nm and
+    therefore β₂ < 0 (anomalous dispersion).
+    """
+    return -dispersion * wavelength**2 / (2.0 * np.pi * C_LIGHT)
+
+
+def propagate_dispersion(
+    field: np.ndarray, sample_rate: float, beta2: float, distance: float
+) -> np.ndarray:
+    """Propagate a complex envelope through pure group-velocity dispersion.
+
+    Solves the linear part of the NLSE, ``dA/dz = -(i*beta2/2) * d2A/dT2``, which
+    in the frequency domain is an exact all-pass phase rotation::
+
+        A(z, w) = A(0, w) * exp(i * beta2 * w**2 * z / 2)
+
+    Because the transfer function has unit magnitude, this conserves energy
+    exactly (up to floating-point) and is exactly invertible by propagating
+    ``-distance`` — both of which are asserted in the test suite.
+
+    Only β₂ is modelled, so the result is insensitive to the sign convention of
+    the Fourier transform (ω appears squared). That stops being true as soon as
+    β₃ or a group-delay term is added, so this note is worth keeping.
+
+    The sign of β₂ itself, however, is *not* free — and the unchirped broadening
+    formula cannot detect an error in it, being even in β₂. Only the chirped case
+    can, which is why ``test_chirped_pulse_compresses_before_broadening`` exists.
+
+    The phase argument reaches thousands of radians over a realistic span, so the
+    transform runs in double precision regardless of the storage precision and
+    the caller casts the result back. Correctness first; if profiling later shows
+    this matters, the precision policy belongs here, in one place.
+    """
+    if distance == 0.0 or beta2 == 0.0:
+        return field.astype(np.complex128, copy=True)
+
+    omega = angular_frequency_grid(field.shape[0], sample_rate)
+    transfer = np.exp(0.5j * beta2 * omega**2 * distance)
+    return np.fft.ifft(np.fft.fft(field.astype(np.complex128)) * transfer)

src/oosim/units.py

@@ -74,6 +74,8 @@ def frequency_to_wavelength(f_hz: float) -> float:
     # Ratios
     "dB": db_to_linear,
     "dB/km": lambda x: x * 1e-3,  # -> dB/m (still logarithmic, per metre)
+    # Dispersion: ps/(nm*km) -> s/m^2
+    "ps/nm/km": lambda x: x * 1e-6,
     # Length
     "m": lambda x: x,
     "km": lambda x: x * 1e3,
@@ -98,6 +100,7 @@ def frequency_to_wavelength(f_hz: float) -> float:
     "dBm": w_to_dbm,
     "dB": linear_to_db,
     "dB/km": lambda x: x * 1e3,
+    "ps/nm/km": lambda x: x * 1e6,
     "m": lambda x: x,
     "km": lambda x: x * 1e-3,
     "nm": lambda x: x * 1e9,