Core physics

Grids and transforms, pulse parameterisations, nonlinear interactions, the forward model, material dispersion and numerical helpers. Begin with measurement fundamentals for units and conventions, or modelling the measurement for the physical assumptions behind these interfaces. The first retrieval shows how the grid, pulse and trace functions fit together.

Grid and transforms

Time/frequency grid and Fourier transforms.

A Grid holds matched, centred time and angular-frequency axes and provides the physics-convention Fourier transform pair used throughout croak:

• forward (time → frequency): \(\tilde E(\omega) = \int E(t)\,e^{+i\omega t}\,dt\)

• inverse (frequency → time): \(E(t) = \tfrac{1}{2\pi}\int \tilde E(\omega)\,e^{-i\omega t}\,d\omega\)

discretised as \(\tilde E \approx \Delta t \sum E\,e^{+i\omega t}\) and \(E \approx \tfrac{\Delta\omega}{2\pi}\sum \tilde E\,e^{-i\omega t}\).

These Grid.fft()/Grid.ifft() operate on centred arrays (zero frequency in the middle). The forward model (see croak.forward) instead works in raw DFT-bin order; the unshifted, unnormalised transforms it uses are exposed as the module-level raw_fft()/raw_ifft().

This NumPy-facing module uses scipy.fft; the JAX forward models use the corresponding JAX transforms with the same conventions.

class croak.grid.Grid(n, *, dt=None, domega=None)

Bases: object

A matched, centred time/angular-frequency grid.

Parameters:

• n (int) – Number of grid points.

• dt (float, optional) – Time step (s). Provide exactly one of dt or domega.

• domega (float, optional) – Angular-frequency step (rad/s). The other step follows from \(\Delta\omega = 2\pi/(\Delta t\,N)\).

omega

Centred angular-frequency axis (rad/s), relative to the carrier.

Type:

numpy.ndarray

t

Centred time axis (s).

Type:

numpy.ndarray

dt, domega, df

Time step, angular-frequency step, and ordinary-frequency step domega / (2*pi).

Type:

float

n

Number of points.

Type:

int

classmethod from_omega(omega)

Construct a Grid matching a uniformly spaced omega axis.

Parameters:

omega (ArrayLike)

Return type:

Grid

fft(et)

Forward transform time → frequency (centred), scaled by dt.

Parameters:

et (ArrayLike)

Return type:

NDArray[complex128]

ifft(ew)

Inverse transform frequency → time (centred), scaled by df.

Parameters:

ew (ArrayLike)

Return type:

NDArray[complex128]

tshift(et, tau)

Shift a time-domain field by tau via a spectral phase ramp.

Applies \(E(t) \to E(t-\tau)\) using \(\tilde E \to \tilde E\,e^{+i\omega\tau}\) in the raw (DFT-bin) domain.

Parameters:

• et (ArrayLike)

• tau (float)

Return type:

NDArray[complex128]

croak.grid.raw_fft(x)

Unnormalised, unshifted forward DFT (scipy.fft.fft).

Used by the forward model as the spectrum→time transform in DFT-bin order.

Parameters:

x (NDArray[complexfloating])

Return type:

NDArray[complex128]

croak.grid.raw_ifft(x)

Inverse DFT with 1/N normalisation (scipy.fft.ifft).

Companion to raw_fft(); the time→frequency transform in DFT-bin order.

Parameters:

x (NDArray[complexfloating])

Return type:

NDArray[complex128]

croak.grid.gridparams_omega(trange, *, omega_min, omega_max, omega0=None)

Grid parameters spanning a time range and frequency range.

Picks the time step from the Nyquist limit of the largest frequency offset from the carrier and the number of points from the requested time window (rounded up to an FFT-friendly length).

Parameters:

• trange (float) – Required temporal window (s).

• omega_min (float) – Extremes of the absolute angular-frequency content (rad/s).

• omega_max (float) – Extremes of the absolute angular-frequency content (rad/s).

• omega0 (float, optional) – Carrier angular frequency; defaults to the midpoint.

Returns:

tuple – (n, dt, omega0).

Return type:

tuple[int, float, float]

croak.grid.gridparams_lambda(trange, *, lambda_min, lambda_max, omega0=None)

As gridparams_omega(), but specified by a wavelength range (m).

Parameters:

• trange (float)

• lambda_min (float)

• lambda_max (float)

• omega0 (float | None)

Return type:

tuple[int, float, float]

Pulses

Pulse representations and test-pulse generators.

Two groups of functionality live here:

• Generators — gaussian_pulse(), sech_pulse() and random_gaussian_spectrum() build complex spectra \(\tilde E(\omega)\) (centred order) for tests, examples and retrieval initial guesses. Chirp and higher-order phase are added as a Taylor expansion of the spectral phase.

• Parameterisations — ComplexPulse and ArrayPulse expose a real parameter vector u to the gradient-based (L-BFGS) solver, together with the hand-derived vector–Jacobian product Pulse.ew_vjp() that maps a spectral cotangent back to u, and optional spectral-smoothness penalties.

The carrier (centre wavelength) is not encoded in the baseband pulse shape: a pulse lives on the grid’s centred omega axis, and the carrier enters absolute-frequency calculations through omega0, including material dispersion, chromatic focusing, collection and wavelength conversion. “Different wavelengths” therefore means the same baseband pulse with a different omega0 in ForwardModel.

croak.pulses.gaussian_pulse(grid, fwhm, *, phases=())

Transform-limited Gaussian pulse with optional spectral-phase chirp.

Parameters:

• grid (Grid) – Time/frequency grid.

• fwhm (float) – Intensity FWHM of the (unchirped) temporal pulse (s).

• phases (sequence of float, optional) – Taylor spectral-phase coefficients [GDD, TOD, ...] (s², s³, …).

Returns:

numpy.ndarray – Complex spectrum (centred order), normalised to unit peak spectral intensity.

Return type:

NDArray[complex128]

croak.pulses.sech_pulse(grid, fwhm, *, phases=())

Transform-limited \(\operatorname{sech}^2\) pulse with optional chirp.

Parameters as in gaussian_pulse(); fwhm is the intensity FWHM.

Parameters:

• grid (Grid)

• fwhm (float)

• phases (ArrayLike)

Return type:

NDArray[complex128]

croak.pulses.random_gaussian_spectrum(grid, fwhm=2.5e-15, *, rng=None)

Random initial-guess spectrum: a Gaussian envelope with random phase.

Builds a Gaussian temporal envelope of the given FWHM, multiplies by a uniform random temporal phase, and transforms to the spectral domain.

Parameters:

• grid (Grid) – Time/frequency grid.

• fwhm (float, optional) – Temporal FWHM of the envelope (s).

• rng (numpy.random.Generator, optional) – Random generator (defaults to a fresh default generator).

Return type:

NDArray[complex128]

class croak.pulses.Pulse

Bases: object

Abstract real parameterisation of a complex spectrum for optimisation.

Subclasses map between a real parameter vector u and the complex spectrum, and provide the hand-derived adjoint ew_vjp().

property nparams: int

Number of real parameters in u.

get_params()

Return the current parameter vector.

Return type:

NDArray[float64]

set_params(u)

Set the pulse from a parameter vector.

Parameters:

u (ArrayLike)

Return type:

None

spectrum()

Return the complex spectrum (centred order).

Return type:

NDArray[complex128]

ew_vjp(ew_bar)

Map a spectral cotangent ew_bar to the parameter gradient dL/du.

ew_bar is the Wirtinger cotangent \(\partial L/\partial\tilde E^*\).

Parameters:

ew_bar (NDArray[complex128])

Return type:

NDArray[float64]

bounds()

Return (lower, upper) parameter bounds.

Return type:

tuple[NDArray[float64], NDArray[float64]]

amplitude_penalty(u)

Spectral-amplitude smoothness penalty (0 by default).

Parameters:

u (ArrayLike)

Return type:

float

amplitude_penalty_grad(u)

Gradient of amplitude_penalty() (zeros by default).

Parameters:

u (ArrayLike)

Return type:

NDArray[float64]

phase_penalty(u)

Spectral-phase smoothness penalty (0 by default).

Parameters:

u (ArrayLike)

Return type:

float

phase_penalty_grad(u)

Gradient of phase_penalty() (zeros by default).

Parameters:

u (ArrayLike)

Return type:

NDArray[float64]

class croak.pulses.ComplexPulse(grid, ew)

Bases: Pulse

Full complex spectrum parameterised by u = [Re E; Im E] (length 2N).

The natural representation for step-based algorithms and a fully general, unconstrained parameterisation for L-BFGS.

Parameters:

• grid (Grid)

• ew (ArrayLike)

property nparams: int

real and imaginary parts).

Type:

Number of real parameters (2N

get_params()

Return u = [Re E; Im E].

set_params(u)

Set the spectrum from u = [Re E; Im E].

spectrum()

Return the complex spectrum (centred order).

ew_vjp(ew_bar)

Map a spectral cotangent to dL/du (real then imaginary blocks).

amplitude_penalty(u)

Second-difference smoothness penalty on the spectral amplitude.

class croak.pulses.ArrayPulse(grid, amplitude, phase=None, *, support=None, phase_only=False)

Bases: Pulse

Amplitude + spectral phase parameterisation.

Parameters u are the spectral amplitude followed by the spectral phase over an optional support window: u = [A[support]; phi[support]] (or just the phase when phase_only=True). This is the standard L-BFGS FROG parameterisation; phase_only retrieves the phase with a fixed amplitude.

Parameters:

• grid (Grid) – Time/frequency grid.

• amplitude (array_like) – Spectral amplitude \(|\tilde E|\) (length N).

• phase (array_like, optional) – Spectral phase (length N); defaults to zero.

• support (tuple of int, optional) – (start, stop) half-open index range of the optimised window; defaults to the full grid.

• phase_only (bool, optional) – Optimise only the phase, keeping amplitude fixed.

classmethod from_spectrum(grid, ew, *, support=None, phase_only=False)

Build an ArrayPulse from a complex spectrum.

Parameters:

• grid (Grid)

• ew (ArrayLike)

• support (tuple[int, int] | None)

• phase_only (bool)

Return type:

ArrayPulse

property nparams: int

Number of real parameters (phase only, or amplitude + phase).

get_params()

Return the parameter vector over the support window.

set_params(u)

Set amplitude/phase over the support window from u.

spectrum()

Return the complex spectrum A exp(i phi).

ew_vjp(ew_bar)

Map a spectral cotangent to dL/du (amplitude/phase blocks).

amplitude_penalty(u)

Second-difference amplitude smoothness penalty (0 if phase-only).

amplitude_penalty_grad(u)

Gradient of amplitude_penalty() (phase block zero).

phase_penalty(u)

Second-difference phase smoothness penalty (mean squared curvature).

phase_penalty_grad(u)

Gradient of phase_penalty().

spectral_penalty(u, target)

Relative spectral-amplitude mismatch to a target spectrum.

Penalises the squared difference between the retrieved spectral amplitude \(|\tilde E(\omega)|\) (the amplitude parameters) and a target amplitude target over the optimised support:

\[P = \lVert a - s\rVert^2 / \max(\lVert s\rVert^2, \varepsilon),\]

a dimensionless relative mismatch. Pulls a full-mode retrieval towards a measured independent spectrum (and softly suppresses amplitude outside the measured band). Returns 0 in phase-only mode (no amplitude parameters). target is the full-grid amplitude \(s\) (length N), sliced here to the optimisation support.

Parameters:

target (ArrayLike)

Return type:

float

spectral_penalty_grad(u, target)

Gradient of spectral_penalty() (amplitude block; phase block zero).

Parameters:

target (ArrayLike)

Return type:

NDArray[float64]

Interactions

Nonlinear FROG interactions and their Wirtinger adjoints.

Each interaction defines the instantaneous nonlinear signal field produced by a (undelayed) test field E and a (delayed) gate field G in the time domain, together with the reverse-mode (Wirtinger) adjoint used to back-propagate a cotangent through that signal. The three standard geometries are supported:

Interaction	signal s	geometry
SHG	E * G	second-harmonic generation
SD	E**2 * conj(G)	self-diffraction
PG	E * |G|**2	polarisation gating

Adjoint convention: the cotangent of a complex variable v is v_bar = dL/dv* (Wirtinger). Interaction.signal_adjoint() returns the cotangents (E_bar, G_bar) for given output cotangent s_bar.

Geometric smearing (croak.smearing) needs the same signals written out as a product of three separately shifted replicas, because the two factors that coincide in the thin 1-D model acquire different arrival times across the focal spot. That form is exposed by Interaction.smeared_shifts() and Interaction.smeared_signal(); at zero smearing they reproduce signal() exactly.

class croak.interactions.Interaction

Bases: object

Base class for FROG nonlinear interactions.

Subclasses set name and implement signal() and signal_adjoint().

omega0_scale: ClassVar[float] = 1.0

Factor relating the trace carrier to the pulse carrier (SHG doubles it).

time_reversal_ambiguous: ClassVar[bool] = False

Whether the trace is invariant under time reversal of the pulse, \(E(t) \to E^*(-t)\) — equivalently conjugating the spectrum. When it is, the retrieval cannot distinguish a pulse from its mirror image and the solver returns whichever branch it happened to land in. True only for SHG, whose signal E*G is symmetric in the two arms; the SD and PG kernels treat their arms differently and so fix the direction of time. See croak.processing.resolve_time_direction().

signal(test, gate)

Return the time-domain nonlinear signal field.

Parameters:

• test (NDArray[complex128])

• gate (NDArray[complex128])

Return type:

NDArray[complex128]

signal_adjoint(test, gate, s_bar)

Return (test_bar, gate_bar) for output cotangent s_bar.

Parameters:

• test (NDArray[complex128])

• gate (NDArray[complex128])

• s_bar (NDArray[complex128])

Return type:

tuple[NDArray[complex128], NDArray[complex128]]

smeared_shifts(tau)

Time shifts of the three replicas under geometric smearing.

The smeared signal is a product of three replicas of E. One of them does not depend on the smearing parameter p; the other two are shifted by varying_base -/+ p/2 and therefore swap under p -> -p. With symmetric quadrature nodes that lets the forward model build each shifted field once and use it twice.

Parameters:

tau (float) – Delay (s).

Returns:

tuple of float – (fixed_shift, varying_base), both in seconds. Replica time shifts are fixed_shift and varying_base -/+ p/2, where a shift s means the replica is E(t - s).

Return type:

tuple[float, float]

smeared_signal(fixed, varying, varying_swapped)

Combine the three smeared replicas into the nonlinear signal.

Parameters:

• fixed (numpy.ndarray) – The replica at smeared_shifts()’ fixed_shift.

• varying (numpy.ndarray) – The replica at varying_base - p/2.

• varying_swapped (numpy.ndarray) – The replica at varying_base + p/2 — i.e. varying evaluated at the mirrored node -p.

Return type:

NDArray[complex128]

class croak.interactions.SHG

Bases: Interaction

Second-harmonic generation: s = E * G.

omega0_scale: ClassVar[float] = 2.0

Factor relating the trace carrier to the pulse carrier (SHG doubles it).

time_reversal_ambiguous: ClassVar[bool] = True

Whether the trace is invariant under time reversal of the pulse, \(E(t) \to E^*(-t)\) — equivalently conjugating the spectrum. When it is, the retrieval cannot distinguish a pulse from its mirror image and the solver returns whichever branch it happened to land in. True only for SHG, whose signal E*G is symmetric in the two arms; the SD and PG kernels treat their arms differently and so fix the direction of time. See croak.processing.resolve_time_direction().

signal(test, gate)

Return the SHG signal E * G.

Parameters:

• test (NDArray[complex128])

• gate (NDArray[complex128])

Return type:

NDArray[complex128]

signal_adjoint(test, gate, s_bar)

Return (test_bar, gate_bar) for the SHG signal E * G.

class croak.interactions.SD

Bases: Interaction

Self-diffraction: s = E**2 * conj(G).

signal(test, gate)

Return the SD signal E**2 * conj(G).

signal_adjoint(test, gate, s_bar)

Return (test_bar, gate_bar) for the SD signal E**2 * conj(G).

smeared_shifts(tau)

Return (tau, 0.0): the conjugated replica is the delayed, fixed one.

The smeared SD signal is E(t + p/2) E(t - p/2) conj(E(t - tau)) — the two squared replicas separate by p while the conjugated one only ever moves with the delay.

smeared_signal(fixed, varying, varying_swapped)

Return E(t+p/2) E(t-p/2) conj(E(t-tau)).

class croak.interactions.PG

Bases: Interaction

Polarisation gating: s = E * |G|**2.

signal(test, gate)

Return the PG signal E * |G|**2.

signal_adjoint(test, gate, s_bar)

Return (test_bar, gate_bar) for the PG signal E * |G|**2.

smeared_shifts(tau)

Return (0.0, tau): the undelayed probe is the fixed replica.

The smeared PG signal is E(t) E(t - tau + p/2) conj(E(t - tau - p/2)) — the gate pair, which coincide in |G|**2 when unsmeared, separate by p.

smeared_signal(fixed, varying, varying_swapped)

Return E(t) E(t-tau+p/2) conj(E(t-tau-p/2)).

croak.interactions.get_interaction(interaction)

Resolve an interaction from a name (case-insensitive) or instance.

Parameters:

interaction (str or Interaction) – "shg", "sd", "pg" (any case) or an Interaction.

Returns:

Interaction

Return type:

Interaction

croak.interactions.INTERACTIONS: dict[str, Interaction] = {'pg': PG(), 'sd': SD(), 'shg': SHG()}

Registry of interaction instances keyed by lowercase name.

Forward model

The unified nonlinear-process forward model and its reverse-mode adjoint.

ForwardModel is the NumPy computational core and its analytic adjoint. The JAX solvers use croak.forward_jax, which follows the same maths and adds the explicit focal and collection models. It maps a complex spectrum \(\tilde E(\omega)\) to a delay-scanned nonlinear-process spectrum (a “PNPS” trace, e.g. FROG or TDP), optionally through a dispersive PG/SD slab, and provides the exact reverse-mode (Wirtinger) adjoint of that map.

Thin medium is simply the special case of the dispersive model with \(\beta = 0\), zero thickness and a single quadrature node — so one code path and one adjoint kernel serve COPRA and L-BFGS.

Dispersive model (depth integral over a slab of thickness L), per delay \(\tau\):

\[\psi(\omega) = \sum_q w_q\, e^{i\beta(L-z_q)} \odot \mathcal{F}^{-1}\!\big[\, \mathrm{sig}( \mathcal{F}[\tilde E e^{i\beta z_q}],\; \mathcal{F}[\tilde E e^{i\omega\tau} e^{i\beta z_q}]) \big]\]

with Gauss–Legendre nodes \(z_q\) and weights \(w_q\) on [0, L]. \(\mathcal{F}\) here is the raw (unshifted) DFT; the model works internally in DFT-bin order and exposes centred arrays at its boundary. The trace is \(T(\omega,\tau)=|\psi(\omega,\tau)|^2\).

class croak.forward.ForwardModel(omega, delays, interaction, *, material=None, thickness=0.0, npoints=1, omega0=0.0, quadrature='gausslegendre', smearing=None)

Bases: object

Nonlinear-process forward map and adjoint for fixed grid/delays/interaction.

Parameters:

• omega (array_like) – Centred angular-frequency grid (rad/s), relative to omega0.

• delays (array_like) – Delay axis (s).

• interaction (str or Interaction) – Nonlinear interaction ("shg", "sd", "pg").

• material (str or None, optional) – Dispersive material for the slab (e.g. "SiO2"); None (default) is a thin, dispersionless medium.

• thickness (float, optional) – Slab thickness (m). Default 0.

• npoints (int, optional) – Number of depth quadrature points (default 1).

• omega0 (float, optional) – Carrier angular frequency (rad/s), used to build beta.

• quadrature (str, optional) – Quadrature method (see quadrature_nodes_weights()).

• smearing (SmearingKernel or None, optional) – Geometric time-smearing kernel of a non-collinear geometry (see croak.smearing). None (default) switches the correction off and selects the local unsmeared model.

Notes

All public array arguments/returns (spectra, cotangents) are in centred frequency order; the DFT-bin ordering is an internal detail.

The hand-derived adjoint_single() is not available with smearing; use the autodiff solvers (see croak.forward_jax) for that case.

signal_single(ew, tau, *, record=False)

Frequency-domain signal \(\psi(\omega)\) for one delay (centred).

Parameters:

• ew (numpy.ndarray) – Complex spectrum (centred order).

• tau (float) – Delay (s).

• record (bool, optional) – If True, store the forward time-domain fields so a following adjoint_single() call (same ew, tau) can reuse them.

Return type:

NDArray[complex128]

signal_smeared(ew, tau)

Per-node signals \(\psi_k(\omega)\) for one delay, shape (N, K).

One column per Gauss–Hermite node of the shape parameter p; the delay offset delta is not applied here because it acts across the whole delay axis (see trace()).

The three replicas are built as exact spectral phases, so no interpolation and no delay-grid alignment is involved. Because the nodes are antisymmetric, the replica at varying_base + p_k/2 is the one already computed for the mirrored node -p_k: the [::-1] below is what keeps the cost at K + 1 field constructions rather than 2K + 1.

Parameters:

• ew (NDArray[complex128])

• tau (float)

Return type:

NDArray[complex128]

trace(ew, *, normalize=True)

Full FROG trace (Nomega, Ndelay) for spectrum ew (centred).

Parameters:

• ew (numpy.ndarray) – Complex spectrum (centred order).

• normalize (bool, optional) – Scale the trace to unit maximum (default True).

Return type:

NDArray[float64]

adjoint_single(psi_bar, tau)

Apply the signal adjoint to a centred field cotangent.

Adjoint of signal_single(): map psi_bar → ew_bar (centred).

Requires the forward tape for the same tau (call signal_single(ew, tau, record=True) first). psi_bar is the Wirtinger cotangent \(\partial L/\partial\psi^*\); the return is \(\partial L/\partial \tilde E^*\) (centred order).

Parameters:

• psi_bar (NDArray[complex128])

• tau (float)

Return type:

NDArray[complex128]

croak.forward.quadrature_nodes_weights(thickness, npoints, method='gausslegendre')

Nodes and weights for the depth integral over [0, thickness].

Parameters:

• thickness (float) – Slab thickness (m). 0 returns the single thin-medium node.

• npoints (int) – Number of quadrature points.

• method ({"gausslegendre", "uniform"}) – "gausslegendre" (default) uses Gauss–Legendre quadrature mapped to [0, L]; "uniform" uses equally spaced nodes with unit weights (a plain sum, useful for convergence checks).

Returns:

tuple of numpy.ndarray – (nodes, weights).

Return type:

tuple[NDArray[float64], NDArray[float64]]

croak.forward.maketrace(omega, delays, ew, interaction, *, material=None, thickness=0.0, npoints=1, omega0=0.0, quadrature='gausslegendre', smearing=None, normalize=True)

Generate a (normalised) FROG trace for spectrum ew.

Convenience wrapper around ForwardModel for tests and examples, and the forward-only entry point for synthesising a geometrically smeared trace without running a retrieval. Returns an (Nomega, Ndelay) array. See ForwardModel for the parameter meanings.

Parameters:

• omega (ArrayLike)

• delays (ArrayLike)

• ew (NDArray[complex128])

• interaction (str | Interaction)

• material (str | None)

• thickness (float)

• npoints (int)

• omega0 (float)

• quadrature (str)

• smearing (SmearingKernel | None)

• normalize (bool)

Return type:

NDArray[float64]

Geometric smearing

Mask geometry to smearing-kernel statistics for a non-collinear BOXCARS setup; consumed by the forward model.

Geometric (pulse-front-tilt) time smearing in a non-collinear BOXCARS geometry.

In a non-collinear geometry the input beams reach the interaction plane with tilted pulse fronts, so the relative arrival time between arms varies linearly across the transverse coordinate \(\mathbf r\) of the focal spot. The detector integrates over that coordinate, so the measured trace is an average over a distribution of internal relative delays rather than a single one. This module turns a mask geometry into the statistics of that distribution; croak.forward and croak.forward_jax consume it.

The reduction

Write the third-order signal at focal-plane position \(\mathbf r\) as the product of three arm fields, one of them phase-conjugated. A beam apertured at mask position \(\mathbf r_j\) and focused by a lens of focal length \(f\) crosses the focus at angle \(\mathbf r_j/f\), so its pulse front makes it arrive at focal position \(\mathbf r\) at time \(t - \boldsymbol\alpha_j\cdot\mathbf r\) with the tilt coefficient

\[\boldsymbol\alpha_j = -\mathbf r_j / (f c).\]

Because \(|\mathcal F\{\cdot\}|^2\) is blind to a global time origin, and the origin may be chosen per \(\mathbf r\), only two scalar combinations of the three arrival times survive. In croak’s own operator conventions they are

\[\begin{split}\mathrm{PG:}\quad S &= E(t)\,E(t-\tau+\delta+p/2)\,E^*(t-\tau+\delta-p/2) \\ \mathrm{SD:}\quad S &= E(t+p/2)\,E(t-p/2)\,E^*(t-\tau+\delta)\end{split}\]

with, for arms given in the interaction’s role order (\(\mathbf r_1,\mathbf r_2,\mathbf r_3\); see kernel_from_arms()),

\[\begin{split}\mathrm{PG:}\quad p &= (\alpha_3-\alpha_2)\cdot\mathbf r, & \delta &= \bigl(\alpha_1-\tfrac12(\alpha_2+\alpha_3)\bigr)\cdot\mathbf r \\ \mathrm{SD:}\quad p &= (\alpha_2-\alpha_1)\cdot\mathbf r, & \delta &= \bigl(\tfrac12(\alpha_1+\alpha_2)-\alpha_3\bigr)\cdot\mathbf r.\end{split}\]

Setting \(p=\delta=0\) recovers test * |gate|**2 and test**2 * conj(gate) exactly, so the correction vanishes identically when the kernel is switched off.

p changes the shape of the signal and needs an explicit quadrature. delta is a pure translation of the delay axis (\(\tau\to\tau-\delta\)), so integrating over it is a convolution along the measured delay axis, applied once to the finished trace.

The focal weight

The local signal amplitude is the product of the three focal-plane amplitudes, so the weight with which position \(\mathbf r\) contributes to the incoherent sum is \(|A(\mathbf r)|^6\). For a uniformly illuminated circular hole of diameter D, \(A\) is the Airy amplitude \(2J_1(u)/u\) with \(u=\pi D r/(\lambda f)\), and the rms of one Cartesian component of \(\mathbf r\) under that weight is

\[\sigma_r = C\,\lambda f / D, \qquad C = 0.245673\ldots\]

(see AIRY6_RMS_COEFF). Since \(\boldsymbol\alpha_j\propto 1/f\), the focal length cancels: the kernel widths depend only on the mask ratio and the wavelength.

Since p and delta are both linear functionals of \(\mathbf r\), replacing the focal weight by an isotropic Gaussian of the same second moment makes (p, delta) jointly Gaussian, fully described by SmearingKernel. The leading trace correction is quadratic in (p, delta), so only those second moments matter to first order; the true \(|A|^6\) weight has heavier tails, which tests/test_smearing.py quantifies.

Validity limits

• Full-beam collection. The incoherent sum over \(\mathbf r\) follows from Parseval and holds only if the detector collects the whole signal beam. Under strong spatial filtering the average is partly coherent, which changes the trace in ways a narrower blur width does not capture. croak.collection lifts this: it transforms the focal field to transverse wavevector before squaring it, so a finite hole is modelled. The reduction to \((p,\vartheta)\) survives that intact (the phase the incoherent sum discards is exactly \(\omega p\), with \(\omega\) the centred frequency), but the kernel in this module remains the full-collection form.

• Achromatic kernel. \(\sigma\propto\lambda\) because the focal spot scales with wavelength while the tilt does not; the kernel is evaluated at the carrier only. Across a 1 fs deep-UV pulse this is a −16 %/+24 % variation of width; such broadband pulses need the chromatic filters of croak.focal.

• Uniform hole illumination, beams collimated at the focusing lens, and a medium thin compared with the Rayleigh range (so the delay-vs-position map is linear).

• No spatial chirp in the input beam.

References

Derivation, the arrangement-dependent widths and the discriminating experiments: docs/howto/geometric_smearing.md. Note that the closed forms differ according to which arm carries the scanned delay; square_boxcars_kernel() covers both cases via its interaction argument.

croak.smearing.collinear_sd_kernel(*, hole_diameter, hole_spacing, wavelength, npoints=5)

Smearing kernel for a two-beam self-diffraction geometry.

For two apertures whose centres lie on one line in the mask plane, the probe entering the nonlinearity twice and the gate once with conjugation, signal at \(2k_E - k_G\). It is NOT what square_boxcars_kernel() with "sd" returns — that closed form assumes the two unconjugated arms come from separate holes, a three-hole layout, and gives a spurious sigma_p here.

The p channel vanishes identically. p is proportional to \(\alpha_A - \alpha_B\), the difference of the two unconjugated arms’ tilts, and in a two-beam SD those arms are the same beam. So

\[\sigma_p = 0, \qquad \sigma_\vartheta = 0.49135\,(d/D)\,\lambda/c, \qquad \rho = 0,\]

verified against kernel_from_arms() to six figures over d/D = 0.75–1.25. Two consequences matter for retrieval. The smearing reduces to a pure one-dimensional convolution along the delay axis, so it can be deconvolved directly — unlike the TG kernel, whose rho = 1/sqrt(5) makes it a correlated two-dimensional shear. And the p channel is the one that shifts the frequency axis in proportion to the gate chirp, so its absence removes the mechanism behind the chirped-arm difficulties of the TG geometry.

The width coefficient is numerically the same 0.49135 that appears as sigma_p for PG, which is a coincidence of the geometry (both reduce to one arm-separation times the focal weight), not a shared mechanism.

Parameters:

• hole_diameter (float) – Aperture diameter D (m).

• hole_spacing (float) – Edge-to-edge gap between the two holes (m); the beam offset from the axis is d = (spacing + D) / 2 and their centre-to-centre separation is 2d.

• wavelength (float) – Carrier wavelength (m) — use the measured trace carrier.

• npoints (int, optional) – Gauss–Hermite nodes. Retained for interface symmetry; with sigma_p = 0 the p quadrature is a single node regardless.

Returns:

SmearingKernel

Return type:

SmearingKernel

croak.smearing.AIRY6_RMS_COEFF = 0.24567310004243084

RMS of one Cartesian component of the focal position under the \(|A|^6\) weight, in units of \(\lambda f / D\). Obtained from the moments of the Airy amplitude \(g(u)=(2J_1(u)/u)^6\):

<|u|^2> = int_0^inf u^3 g du / int_0^inf u g du = 1.19137
sigma_u = sqrt(<|u|^2> / 2) = 0.771805        (one Cartesian component)
C       = sigma_u / pi = 0.245673             (u = pi D r / (lambda f))

Both integrals converge: \(g\sim u^{-9}\) at large u. tests/test_smearing.py recomputes this from scipy’s Bessel functions.

croak.smearing.deconvolve_delay(trace, delays, sigma, *, reg=0.001)

Remove a Gaussian blur of width sigma from the delay axis of a trace.

The \(\vartheta\) channel of the smearing kernel enters the trace as a convolution along delay, and it factorises out of the \(p\) mixture exactly. Writing the joint kernel as \(p\) times \(\vartheta\) conditioned on \(p\) (SmearingKernel.conditional_delay()), the conditional mean \(\mu(p)\) is a per-node delay shift while the conditional width \(\sigma_c = \sigma_\vartheta\sqrt{1-\rho^2}\) is the same for every node. It therefore commutes with the sum:

\[\tilde T(\omega,\tau) = \mathcal N(0,\sigma_c) \ast \Big[\textstyle\sum_k w_k\, \tilde T_{p_k}(\omega, \tau - \mu_k)\Big].\]

Deconvolving a measurement by \(\sigma_c\) thus leaves exactly the \(p\) mixture with its delay shifts — no approximation. Note the width to remove is the conditional one, not the marginal \(\sigma_\vartheta\); using the marginal over-corrects by \(1/\sqrt{1-\rho^2}\) (12% for the square BOXCARS \(\rho = 1/\sqrt5\)).

This matters because it moves the delay-axis half of the smearing out of the forward model and into preprocessing, where any algorithm can use it — including projection-based ones such as COPRA, whose local step is per-delay and cannot represent a convolution that couples delays.

Deconvolution is ill-posed: dividing by a Gaussian amplifies high delay frequencies without bound, so a Wiener filter is used, trading exactness for noise control via reg.

Parameters:

• trace (array_like) – Trace (Nomega, Ndelay).

• delays (array_like) – Delay axis (s). Must be uniformly spaced.

• sigma (float) – Gaussian width to remove (s). Zero or negative returns the trace unchanged.

• reg (float, optional) – Wiener regularisation as a fraction of the filter’s peak power. Larger leaves more blur but suppresses noise amplification. Default 1e-3.

Returns:

numpy.ndarray – Deconvolved trace, same shape.

Return type:

NDArray[float64]

class croak.smearing.SmearingKernel(sigma_p, sigma_delta, rho=0.0, npoints=5)

Bases: object

Bivariate-Gaussian geometric-smearing kernel in (p, delta).

p separates the gate pair for PG (one conjugated) or the two unconjugated replicas for SD. delta shifts their relative delay with respect to the remaining replica, a pure translation of the delay axis. p changes the shape of the nonlinear signal. Both are zero-mean.

Parameters:

• sigma_p (float) – RMS of the shape parameter p (s). Must be >= 0.

• sigma_delta (float) – RMS of the delay offset delta (s). Must be >= 0.

• rho (float, optional) – Correlation coefficient between p and delta, in (-1, 1). It is zero only for particular mask layouts (see square_boxcars_kernel()).

• npoints (int, optional) – Number of Gauss–Hermite nodes along p (default 5). The delta integral is exact and needs no nodes. Cost scales as (npoints + 1) / 2 times the unsmeared model.

Examples

>>> k = SmearingKernel(sigma_p=0.32e-15, sigma_delta=0.36e-15, rho=0.447)
>>> p, w = k.nodes_weights()
>>> float(w.sum())
1.0
>>> bool(np.allclose(p, -p[::-1]))
True

nodes_weights()

Gauss–Hermite nodes and normalised weights for the p integral.

Uses the probabilists’ Hermite rule, whose weight function is \(e^{-x^2/2}\), so the nodes are simply sigma_p * x. The nodes are forced to be exactly antisymmetric (p[k] == -p[-1-k]) and the weights exactly symmetric, which the forward model relies on to reuse each shifted field twice.

Returns:

tuple of numpy.ndarray – (p, w), each of length npoints, with w.sum() == 1.

Return type:

tuple[NDArray[float64], NDArray[float64]]

conditional_delay(p)

Mean and standard deviation of delta conditioned on p.

For a bivariate normal, delta | p is normal with mean \(\rho\,(\sigma_\delta/\sigma_p)\,p\) and standard deviation \(\sigma_\delta\sqrt{1-\rho^2}\). Splitting the joint integral this way lets the delta direction stay an exact convolution even when the two are correlated.

Parameters:

p (numpy.ndarray) – Quadrature nodes along p (s).

Returns:

tuple – (mu, sigma) — the per-node conditional means (s, same shape as p) and the single conditional standard deviation (s).

Return type:

tuple[NDArray[float64], float]

scaled(factor)

Return a copy with both widths multiplied by factor.

Scaling the kernel is equivalent to scaling the mask ratio d/D, which is the single physical lever on the effect. This is the parameterisation used when the smearing width is fitted.

Parameters:

factor (float) – Non-negative multiplier on sigma_p and sigma_delta.

Return type:

SmearingKernel

croak.smearing.delay_frequency_grid(delays, sigma_delta)

Angular-frequency grid conjugate to a uniform delay axis.

The delta integral is a convolution along tau, and it is applied by multiplying the delay-axis DFT of the trace by the analytic Gaussian transform \(e^{-i\Omega\mu-\sigma^2\Omega^2/2}\). Doing it in the transform domain rather than sampling the kernel on the delay grid makes the operation exact for the periodic band-limited representation on the delay grid (which a Nyquist-sampled trace has), so there is no requirement that the kernel span several delay steps.

The convolution is circular. This function raises if the kernel is not comfortably narrower than the scan; the trace must also decay before the delay boundaries, or its ends wrap into each other.

Parameters:

• delays (array_like) – Delay axis (s). Must be uniformly spaced and have at least two points.

• sigma_delta (float) – RMS delay offset (s), used only for the sanity check above.

Returns:

numpy.ndarray – Omega (rad/s) in DFT-bin order, the same length as delays.

Raises:

ValueError – If the delay axis is too short, is not uniformly spaced, or is not wide enough to contain the kernel.

Return type:

NDArray[float64]

croak.smearing.kernel_from_arms(r1, r2, r3, *, interaction, hole_diameter, wavelength, npoints=5)

Build a SmearingKernel from three mask hole positions.

Fully general: any three-arm layout, including unequal spacings. Positions are given in the role order of the interaction, which is what fixes the meaning of p and delta:

• "pg" — (probe, gate unconjugated, gate conjugated). The probe is the arm that appears undelayed in E(t) |E(t - tau)|**2.

• "sd" — (unconjugated A, unconjugated B, conjugated). The conjugated arm is the delayed one in E(t)**2 conj(E(t - tau)).

Parameters:

• r1 (array_like) – Mask hole centres (x, y) in metres, in the role order above. Only differences matter, so the origin is arbitrary.

• r2 (array_like) – Mask hole centres (x, y) in metres, in the role order above. Only differences matter, so the origin is arbitrary.

• r3 (array_like) – Mask hole centres (x, y) in metres, in the role order above. Only differences matter, so the origin is arbitrary.

• interaction (str or Interaction) – "pg" (also used for TG) or "sd". SHG raises.

• hole_diameter (float) – Mask hole diameter (m); sets the focal spot size.

• wavelength (float) – Carrier wavelength (m).

• npoints (int, optional) – Gauss–Hermite nodes along p (default 5).

Returns:

SmearingKernel

Return type:

SmearingKernel

Notes

The focal length never appears: the pulse-front tilt scales as \(1/f\) and the focal spot as \(f\), so the product is focal-length independent.

croak.smearing.square_boxcars_kernel(interaction, *, hole_diameter, hole_spacing, wavelength, npoints=5, reverse_delay=False)

Kernel for the standard folded square BOXCARS mask.

Four holes of diameter D at the corners (±d, ±d), with d = (hole_spacing + D) / 2, uniformly illuminated. The closed forms are

interaction	sigma_p	sigma_delta	rho
"pg"	0.49135 (d/D) lambda/c	0.54933 (d/D) lambda/c	1/√5
"sd"	0.69487 (d/D) lambda/c	0.34743 (d/D) lambda/c	0

so the only physical lever is the mask ratio d/D (and the wavelength).

Parameters:

• interaction (str or Interaction) – "pg" (also used for TG) or "sd"; this also selects which arm carries the scanned delay — see kernel_from_arms().

• hole_diameter (float) – Mask hole diameter D (m).

• hole_spacing (float) – Edge-to-edge gap between adjacent holes (m), so d = (spacing + D) / 2.

• wavelength (float) – Carrier wavelength (m).

• npoints (int, optional) – Gauss–Hermite nodes along p (default 5).

• reverse_delay (bool, optional) – Negate rho (default False). The sign of the p–delta correlation is tied to the orientation of the delay axis: a trace whose tau runs the other way relative to croak’s operator needs the opposite sign. It affects only a small tau-asymmetry of the smeared trace, not the bulk broadening.

Returns:

SmearingKernel

Return type:

SmearingKernel

Examples

An example DUV geometry (1 mm holes, 0.5 mm edge-to-edge, 260 nm):

>>> k = square_boxcars_kernel("pg", hole_diameter=1e-3, hole_spacing=0.5e-3,
...                           wavelength=260e-9)
>>> round(k.sigma_delta * 1e15, 3)
0.357
>>> round(k.sigma_p * 1e15, 3)
0.32

croak.smearing.square_boxcars_delay_width(interaction, *, hole_diameter, hole_spacing, wavelength)

Delay-offset width sigma_delta (s) of a square BOXCARS mask.

Thin wrapper over square_boxcars_kernel() for callers (the GUI) that need the geometry → width map on its own. See that function for the parameters.

Parameters:

• interaction (str | Interaction)

• hole_diameter (float)

• hole_spacing (float)

• wavelength (float)

Return type:

float

croak.smearing.square_boxcars_spacing(interaction, *, hole_diameter, sigma_delta, wavelength)

Invert square_boxcars_delay_width() for the edge-to-edge hole spacing.

sigma_delta is proportional to d = (spacing + D) / 2, so a single reference evaluation fixes the slope and no constants are duplicated. The result may be negative when the requested width is below the geometric floor reached with the holes touching (spacing = 0); callers should clamp.

Parameters:

• interaction (str or Interaction) – "pg" or "sd".

• hole_diameter (float) – Mask hole diameter D (m).

• sigma_delta (float) – Target delay-offset width (s).

• wavelength (float) – Carrier wavelength (m).

Returns:

float – Edge-to-edge hole spacing (m).

Return type:

float

Chromatic focal mixture

An explicit quadrature over the focal plane, with wavelength-dependent amplitude filters and arrival offsets at each node. It resolves effects approximated by the reduced Gaussian kernel; matching second moments does not make the models identical for every field.

Chromatic focal mixture: the smearing kernel without its achromatic reduction.

croak.smearing reduces the incoherent sum over focal position \(\mathbf r\) to two arrival-time parameters \((p,\vartheta)\) carried by a bivariate Gaussian, with the focal weight \(|A(\mathbf r)|^6\) evaluated at the carrier wavelength alone. That reduction keeps the phases the focal average imposes and discards its amplitude structure. Because the mask maps to transverse wavevector as \(x = k_\perp z_{\text{mask}} c/\omega\), the focal amplitude \(A(r,\omega_{\rm abs})\) depends on the product \(r\omega_{\rm abs}\), and the local field at radius \(r\) is the retrieved field passed through a frequency-dependent filter. Different focal positions see different spectra, hence different local pulses, and the measured trace is a generation-weighted mixture of them.

For a transform-limited pulse that costs a few percent of duration. For a chirped one it is much larger, because reshaping the spectrum changes how far the local pulse stretches: in our earlier full 3D simulations, an apertured instrument under-reported a \(+2\,\mathrm{fs}^2\) pulse by 23 %, against 8 % for an otherwise identical instrument whose beams carry no chromatic transverse structure.

This module keeps the mixture explicit. The focal plane is sampled by a quadrature and each node carries its area weight and

• the chromatic amplitude filter \(A(r,\omega_{\rm abs})\), applied to the pulse spectrum before transforming to time and forming the nonlinear product; and

• the definite arrival-time offsets \((p,\vartheta)\) of that node, in place of the reduced model’s distribution over them.

For an achromatic filter the three-arm product supplies \(A^3\) and the trace \(|A|^6\), the weight of the reduced model; a chromatic filter cannot be factored outside the nonlinear time-domain product, which is the point of the mixture. When \(A\) is held at the carrier, the second moments of \((p,\vartheta)\) under the \(|A|^6\) weight are exactly the closed-form \(\sigma_p\), \(\sigma_\vartheta\), \(\rho\) of croak.smearing, which reduced_moments() checks (matching moments does not make the Gaussian kernel identical to the full Airy-weighted distribution).

Without a collection aperture the model sums node intensities, which is exact for full-beam collection. With croak.collection it restores their relative phase, transforms to the aperture and forms the selected detector observable.

The filters have unit on-axis amplitude, so the pulse must be supplied in the on-axis field reference; see airy_amplitude(), Conventions and Modelling the collection aperture.

Cost is one forward evaluation per node instead of one per Gauss–Hermite \(p\)-node: roughly 40 times the reduced model per iteration at the default 24 x 16 node grid.

class croak.focal.EvolveSpec(material, thickness, npoints, quadrature='gausslegendre')

Bases: object

Slab parameters for the depth-resolved mixture (fc-z).

The evolved per-arm filters need the depth-quadrature nodes, which are set by the slab; storing the parameters here (rather than a prebuilt table) lets the table be built at model-build time on the model’s own frequency grid, and lets croak.forward_jax.make_param_trace_fn() verify that the mixture and the trace map agree about the slab — a silent mismatch between the filter depths and the propagation depths would be a wrong model with no error message.

Parameters:

• material (str)

• thickness (float)

• npoints (int)

• quadrature (str)

material

Slab material name (must match the trace map’s material).

Type:

str

thickness

Slab thickness in metres (must match the trace map’s thickness).

Type:

float

npoints

Depth-quadrature node count (must match the trace map’s npoints).

Type:

int

quadrature

Depth-quadrature rule (must match the trace map’s quadrature).

Type:

str

class croak.focal.FocalMixture(r, phi, area, p, theta, hole_diameter, f_foc, profile='airy', waist=0.0, arms=None, collection=None, evolve=None)

Bases: object

Quadrature over the focal plane, carrying the chromatic amplitude filter.

Parameters:

• r (NDArray[float64])

• phi (NDArray[float64])

• area (NDArray[float64])

• p (NDArray[float64])

• theta (NDArray[float64])

• hole_diameter (float)

• f_foc (float)

• profile (str)

• waist (float)

• arms (NDArray[float64] | None)

• collection (CollectionAperture | None)

• evolve (EvolveSpec | None)

r

Radial coordinate of each node (m).

Type:

ndarray, shape (K,)

phi

Azimuth of each node (rad).

Type:

ndarray, shape (K,)

area

Quadrature weight including the r dr dphi Jacobian (m^2). The \(|A|^6\) weight is NOT folded in here: it arrives through spectral_filter(), which is the whole point of the model.

Type:

ndarray, shape (K,)

p

Gate-splitting arrival offset at each node (s).

Type:

ndarray, shape (K,)

theta

Delay-axis arrival offset at each node (s).

Type:

ndarray, shape (K,)

hole_diameter

Aperture diameter (m).

Type:

float

f_foc

Focal length (m).

Type:

float

arms

Mask hole centres (m) in the interaction’s role order, as passed to croak.smearing.kernel_from_arms(). The incoherent sum never needs them – only the differences that make p and theta – but a finite collection aperture does, because the absolute tilts fix where the signal goes (see croak.collection). None for a mixture built without them.

Type:

ndarray, shape (3, 2), or None

collection

Finite collection aperture. None (default) is the full-beam incoherent sum for complete collection; anything else transforms the focal field to k, windows it and integrates there.

Type:

CollectionAperture or None

evolve

Depth-resolved mixture (fc-z). None (default) is the entrance-face model: one filter per node, shared by all three arms and all depth nodes. An EvolveSpec makes each arm at each depth node carry its own linearly evolved complex profile — see evolved_arm_filters(). The (p, theta) offsets stay depth-independent either way (exact: transverse wavevector is conserved).

Type:

EvolveSpec or None

property nodes: int

Number of quadrature nodes.

property signal_position: NDArray[float64]

Mask-plane position of the phase-matched signal, \(r_1+r_2-r_3\) (m).

The degenerate carrier construction places the signal along \(\mathbf k_1+\mathbf k_2-\mathbf k_3\), which maps back to this mask position. For the square BOXCARS layout it is the fourth corner – the one the collection hole sits on.

Raises:

ValueError – If the mixture was built without arms.

spectral_filter(omega)

Amplitude filter \(A(r_k, \omega)\) for every node.

Parameters:

omega (ndarray, shape (N,)) – ABSOLUTE angular frequency (rad/s), not the centred grid: the chromatic scaling is in \(r\omega\), so an offset grid would mis-scale it.

Returns:

ndarray, shape (K, N) – One filter per node. Apply it to each arm spectrum before transforming to time and multiplying the three fields. Only an achromatic common filter reduces to an external \(A^3\) amplitude factor.

Return type:

NDArray[float64]

reduced_weight(wavelength)

\(|A|^6\) weight at one wavelength, times the area element.

This is what the reduced kernel of croak.smearing integrates against; it exists here so the two models can be compared on identical quadratures.

Parameters:

wavelength (float)

Return type:

NDArray[float64]

croak.focal.airy_amplitude(u)

Evaluate the unit-peak Airy amplitude \(2J_1(u)/u\).

Parameters:

u (ndarray) – Reduced coordinate \(u=\pi D r/(\lambda f)=D r\omega_{\rm abs}/(2fc)\).

Returns:

ndarray – Real field amplitude, with the removable singularity set to one at zero.

Return type:

NDArray[float64]

Notes

The normalisation is unit on-axis amplitude, not unit power. Because the focal spot shrinks as \(1/\omega_{\rm abs}\), the power integral carries a frequency dependence, \(\int |A|^2\,\mathrm d^2r\propto\omega_{\rm abs}^{-2}\), whereas an unclipped aperture physically passes the same power at every frequency. The field that multiplies these filters is therefore the ON-AXIS field, not the spatially integrated one. An integrated beamlet spectrum converts to it with the field factor \(\omega_{\rm abs}/\omega_0\) (intensity \((\omega_{\rm abs}/\omega_0)^2\)), where \(\omega_0\) is the pulse’s carrier frequency: this is spectrum_frame_p = 1 in croak.session.pipeline.apply_spectrum_frame(). Weighting a truth by \(\sum_k a_k |A_k(\omega)|^2\) instead imposes this convention’s \(\omega^{-2}\) as if it were a chromatic transmission and shortens the apparent pulse. See Conventions.

At fixed frequency the relative weighting between radii, the chromatic focal structure the mixture exists to represent, is exactly \(|A(r,\omega)|^2\). The local filter changes the pulse inside the nonlinear product, so it is not equivalent to a single output spectral response that R_omega can absorb. The Airy shape describes focusing of a uniformly illuminated, unclipped hole: it carries chromatic focal structure but no chromatic vignetting, which an instrument whose mask clips a guided mode also has and which needs its own treatment (e.g. a beamlet reference spectrum).

croak.focal.evolved_arm_filters(mixture, omega, omega0)

Spec-driven fc-z filter table on the mixture’s own depth quadrature.

Thin wrapper over evolved_profile_table() using the depth nodes of the mixture’s EvolveSpec — the same croak.forward.quadrature_nodes_weights() nodes the trace map propagates on, which is what makes the per-depth filters line up with input_prop column for column.

Parameters:

• mixture (FocalMixture)

• omega (NDArray[float64])

• omega0 (float)

Return type:

NDArray[complex128]

croak.focal.evolved_profile_table(mixture, omega, omega0, z_nodes, n_kperp=256)

Per-arm, per-depth complex focal profiles at explicit depths (fc-z).

Between generation events the beamlet evolution is linear, so the profile of arm \(j\) at depth \(z\) is the entrance profile propagated in transverse wavevector and evaluated at the arm’s walked radius:

\[A_j(r_k, \omega; z) \;=\; N(\omega) \int_0^{k_{\max}(\omega)} S(k_\perp, \omega)\, J_0\!\bigl(k_\perp \rho_{jk}(z)\bigr)\, e^{\,i\,[k_z(\omega,k_\perp) - k_z(\omega,0)]\, z}\, k_\perp\, \mathrm{d}k_\perp , \qquad \rho_{jk}(z) = \lvert \mathbf r_k - \bm\delta_j(z)\rvert ,\]

with the in-glass walk \(\bm\delta_j(z) = -z\,\mathbf r_j/(f\,n_0)\) (Snell: transverse wavevector conserved, so the ray angle divides by the index; the sign follows the +r_j/(f c) tilt convention of _arm_offsets) and \(k_z = \sqrt{(n_0\,\omega/c)^2 - k_\perp^2}\) in the glass — the full square root, computed cancellation-safely as \(k_z - k_{z0} = -k_\perp^2/(k_{z0} + \sqrt{k_{z0}^2 - k_\perp^2})\). Only the TRANSVERSE part \(k_z - k_z(\omega,0)\) appears: the on-axis dispersion stays in the trace map’s input_prop, so there is no double counting and \(z \to 0\) reduces exactly to the entrance filter.

The entrance k-space \(S\) is the aperture’s: a unit disc of radius \(k_R(\omega) = (\omega/c)\,(D/2)/f\) for profile="airy" (whose closed form at \(z=0\) is airy_amplitude()), or the Gaussian \(e^{-k_\perp^2 w_0^2/4}\) for profile="gaussian" (closed form \(e^{-\rho^2/w_0^2}\)). \(N\) normalises to unit peak at \((\rho, z) = (0, 0)\), preserving airy_amplitude()’s unit-peak convention and its documented \(\omega^{-2}\) caveat.

Numerically the table is built in the DIFFERENCE form

\[A_j = A_0(\rho_{jk}(z), \omega) \;+\; N \int S J_0 k_\perp \bigl(e^{i[\cdot]z} - 1\bigr)\, \mathrm{d}k_\perp ,\]

with \(A_0\) the closed-form entrance profile: the quadrature error of the Gauss–Legendre \(k_\perp\) integral cancels between the two terms as \(z \to 0\), so a zero-phase/zero-walk table is BIT-identical to the entrance filter — the reduction the tests assert.

Two approximations are used: the glass index is evaluated at the carrier (\(n_0 = n(\omega_0)\)) in both the walk and \(k_z\) (the chromatic correction is second order in the small quantities), and bins with \(\omega_{\rm abs} \le 0\) or a non-finite index keep the unevolved entrance profile (they carry no signal; a NaN would poison the sums — the croak.materials.beta() policy).

Parameters:

• mixture (FocalMixture) – Must carry arms (the walk needs the hole centres).

• omega (ndarray, shape (N,)) – CENTRED angular-frequency grid (rad/s), as the forward model uses.

• omega0 (float) – Carrier angular frequency (rad/s).

• z_nodes (ndarray, shape (Q,)) – Depths (m) at which to evaluate the evolved profiles.

• n_kperp (int, optional) – Gauss–Legendre nodes of the transverse-wavevector integral. The evolved profiles oscillate more than the entrance Airy; 256 holds the propagator-truth test at ~1e-6 of peak over 40 µm. Refine it for thicker slabs.

Returns:

ndarray, shape (3, K, Q, N), complex – One profile per (arm, node, depth) on the centred frequency axis, in the PG arm order of _ARMS_PG (probe, gate, conjugated gate).

Return type:

NDArray[complex128]

croak.focal.focal_mixture(*, hole_diameter, hole_spacing, f_foc, wavelength, n_radial=24, n_azimuth=16, r_max_units=6.0, profile='airy', waist=None, collection=None, evolve_profiles=False, material=None, thickness=0.0, npoints=0, quadrature='gausslegendre')

Build a focal-plane quadrature for the square BOXCARS PG/TG geometry.

Parameters:

• hole_diameter (float) – Aperture diameter and edge-to-edge gap (m). The hole-centre offset is d = (hole_spacing + hole_diameter)/2, matching croak.smearing.

• hole_spacing (float) – Aperture diameter and edge-to-edge gap (m). The hole-centre offset is d = (hole_spacing + hole_diameter)/2, matching croak.smearing.

• f_foc (float) – Focal length (m).

• wavelength (float) – Carrier wavelength (m), used ONLY to set the radial extent of the quadrature in physical units. The model itself is chromatic.

• n_radial (int, optional) – Quadrature nodes. The azimuthal integrand is smooth and periodic, so uniform sampling converges spectrally and 16 nodes is generous for the incoherent sum; the radial integrand carries the Airy oscillations and wants more. A finite collection aperture can need finer sampling of both.

• n_azimuth (int, optional) – Quadrature nodes. The azimuthal integrand is smooth and periodic, so uniform sampling converges spectrally and 16 nodes is generous for the incoherent sum; the radial integrand carries the Airy oscillations and wants more. A finite collection aperture can need finer sampling of both.

• r_max_units (float, optional) – Radial cut in units of \(\lambda f/D\). \(|A|^6\) falls as \(u^{-9}\), so 6 captures the weight to ~1e-4. Note that a coherent collection integrand is not \(|A|^6\), so a mixture carrying a collection may want more radial nodes – see croak.collection.

• collection (CollectionAperture or None, optional) – Finite collection aperture (croak.collection). None (default) keeps the full-beam incoherent sum.

• evolve_profiles (bool, optional) – Depth-resolved mixture (fc-z): each arm at each depth node carries its own linearly evolved complex profile, evaluated at the arm’s walked radius, instead of the single shared entrance-face filter. Requires the slab parameters below, which must MATCH the ones later passed to croak.forward_jax.make_param_trace_fn() (verified there). See evolved_arm_filters() for the physics and conventions.

• material (str or None, optional) – Slab material for evolve_profiles (in-glass evolution and walk).

• thickness (float, optional) – Slab thickness in metres for evolve_profiles.

• npoints (int, optional) – Depth-quadrature node count for evolve_profiles.

• quadrature (str, optional) – Depth-quadrature rule for evolve_profiles.

• profile (str)

• waist (float | None)

Returns:

FocalMixture

Return type:

FocalMixture

croak.focal.mixture_from_nodes(*, x, y, area, arms, hole_diameter, f_foc, profile='airy', waist=None, collection=None)

Build a mixture on arbitrary focal-plane nodes from three arm positions.

focal_mixture() is this with a polar Gauss–Legendre quadrature and the square BOXCARS layout filled in. Splitting them out keeps the (p, \vartheta) algebra in one place: a Cartesian node set (which the collection tests need, to compare against a plain 2-D FFT and to check Parseval exactly) then cannot drift from it.

Parameters:

• x (ndarray, shape (K,)) – Focal-plane node positions (m).

• y (ndarray, shape (K,)) – Focal-plane node positions (m).

• area (ndarray, shape (K,)) – Quadrature weights including any Jacobian (m^2).

• arms (ndarray, shape (3, 2)) – Mask hole centres (m) in the PG role order (probe, gate unconjugated, gate conjugated). Only differences enter (p, theta), but the absolute positions fix the signal direction and so the aperture centring.

• hole_diameter (float) – Mask hole diameter and focal length (m).

• f_foc (float) – Mask hole diameter and focal length (m).

• profile ({"airy", "gaussian"}, optional) – Focal amplitude model; see FocalMixture.spectral_filter().

• waist (float or None, optional) – 1/e^2 radius (m); required for profile="gaussian".

• collection (CollectionAperture or None, optional) – Finite collection aperture (croak.collection).

Returns:

FocalMixture

Return type:

FocalMixture

croak.focal.reduced_moments(mix, wavelength)

Second moments of \((p,\vartheta)\) under the achromatic weight.

Self-test of the geometry and quadrature: these must reproduce the closed forms of croak.smearing — \(\sigma_p = 0.49135\,(d/D)\lambda/c\), \(\sigma_\vartheta = 0.54933\,(d/D)\lambda/c\), \(\rho = 1/\sqrt5\) — to the quadrature’s accuracy. Agreement means the mixture samples the same focal average the reduced kernel models. It checks second moments only: the full spatial distribution is not Gaussian, and a finite collection aperture adds coherent effects the reduced kernel lacks.

Returns:

dict – sigma_p, sigma_theta (s) and rho.

Parameters:

• mix (FocalMixture)

• wavelength (float)

Return type:

dict[str, float]

Collection aperture

The focal mixture without its full-beam-collection assumption. See Modelling the collection aperture.

Finite collection aperture: the focal mixture without full-beam collection.

croak.focal sums incoherently over focal-plane nodes. By Parseval that is exact only if the detector collects the entire signal beam. For a physical amplitude transmission \(W\), the integrated spectral intensity is

\[T(\omega,\tau)=\frac{1}{(2\pi)^2}\int W^2(\mathbf k,\omega)\,\bigl|\widetilde S(\mathbf k,\omega,\tau)\bigr|^2\, \mathrm d^2k ,\]

with \(\widetilde S\) the transverse Fourier transform of the focal-plane signal field, and only complete collection turns that into \(\int|S|^2\,\mathrm d^2r\). A finite aperture must transform before squaring. This module supplies the aperture; the transform itself lives in croak.forward_jax.make_param_trace_fn().

The phase the incoherent sum throws away

Write the arrival-time offset of arm \(j\) at focal position \(\mathbf r\) as \(\delta t_j=\boldsymbol\alpha_j\cdot\mathbf r\) with \(\boldsymbol\alpha_j=\mathbf r_j/(fc)\), as croak.focal.focal_mixture() does (the opposite overall sign to croak.smearing — a pure point inversion of the focal plane, unobservable for the centred circular aperture modelled here). Each local arm is the pulse spectrum filtered by \(A(\mathbf r,\omega_{\rm abs})\) and transformed to time, \(E_{\mathbf r}\), and the exact local PG signal is

\[S(t,\mathbf r) = E_{\mathbf r}(t-\delta t_1)\, E_{\mathbf r}(t-\tau-\delta t_2)\,E_{\mathbf r}^*(t-\tau-\delta t_3).\]

croak’s per-node build is \(E_{\mathbf r}(t)E_{\mathbf r}(t-\tau+\vartheta+p/2) E^*_{\mathbf r}(t-\tau+\vartheta-p/2)\) with \(\vartheta=\delta t_1-\tfrac12(\delta t_2+\delta t_3)\) and \(p=\delta t_3-\delta t_2\). The two are identical up to a rigid time shift, \(S(t,\mathbf r)=S_{\text{croak}}(t-\delta t_1,\mathbf r)\), so in croak’s convention (a delay \(s\) is a spectral factor \(e^{+i\omega s}\))

\[\Psi(\omega,\mathbf r) = \Psi_{\text{croak}}(\omega,\mathbf r)\, e^{+i\omega\,\boldsymbol\alpha_1\cdot\mathbf r}.\]

The \((p,\vartheta)\) reduction therefore loses nothing: the incoherent assumption enters only in discarding that one phase. Removing the signal’s own phase-matched carrier direction \(\boldsymbol\alpha_s=\boldsymbol\alpha_1+\boldsymbol\alpha_2-\boldsymbol\alpha_3\) leaves \(\boldsymbol\alpha_1-\boldsymbol\alpha_s=\boldsymbol\alpha_3- \boldsymbol\alpha_2\), which is precisely the p vector, so

\[\widetilde S(\boldsymbol\kappa,\Omega,\tau) =\sum_k a_k\,\Psi_{\text{croak}}(\Omega,\tau,\mathbf r_k)\, e^{+i\Omega p_k}\,e^{-i\boldsymbol\kappa\cdot\mathbf r_k}, \qquad \boldsymbol\kappa=\mathbf k-\omega_{\rm abs}\boldsymbol\alpha_s .\]

The discarded phase is exactly the frequency times one of the two reduced smearing parameters: p for PG and, by the same algebra, \(-\vartheta\) for SD, which is why only PG is implemented here.

Warning

That frequency is the centred envelope frequency \(\Omega=\omega_{\rm abs}-\omega_0\), while the aperture’s own position \(\boldsymbol\kappa\) scales with the absolute one. croak’s forward model applies every time shift as \(e^{+i\Omega s}\) on a centred spectrum, so each replica drops a constant \(e^{+i\omega_0 s}\); the three together drop \(e^{-i\omega_0 p}\), which \(|\Psi|^2\) never sees but a coherent sum does. Written against the absolute frequency the transform phase is \(\omega_{\rm abs}(p_k-\mathbf u_j\cdot\mathbf r_k)-\omega_0 p_k\), and that trailing term is the whole difference between a model that reproduces a dense 2-D FFT and one that is wrong by four orders of magnitude. See Modelling the collection aperture for the arm conventions.

Where the aperture sits

A hole at mask-plane position \(\mathbf x\) a distance \(z_{\rm mask}\) from the focus selects transverse wavevectors about \(\mathbf k=(\omega_{\rm abs}/c)\mathbf x/z_{\rm mask}\), so in the carrier-removed frame the aperture is centred at \((\omega_{\rm abs}/c)(\mathbf x_{\text{hole}}-\mathbf r_s)/z_{\rm mask}\) with half-width \((\omega_{\rm abs}/c)(D_{\text{hole}}/2)/z_{\rm mask}\): both scale with frequency, i.e. the aperture is chromatic. For a square BOXCARS mask collected at the signal corner (the degenerate carrier direction) the centre is identically zero, which tests/test_collection.py asserts rather than assumes.

Two collection observables are supported, both functionals of the same \(\widetilde S\):

"integrated"

\(\sum_j \mathrm dk_j\,W_j^2\,|\widetilde S_j|^2/(2\pi)^2\) — a spectrometer fed by everything the hole passes (ModelPNPS.jl’s Iω_win).

"reimaged"

\(\bigl|\sum_j \mathrm dk_j\,W_j\,\widetilde S_j/(2\pi)^2\bigr|^2\) — the on-axis re-imaged pixel, i.e. the windowed field back at the focal origin, and the fully coherent limit (Iω_win_reimaged).

Transmission and quadrature

With weighting="amplitude" the physical field transmission enters the integrated intensity as \(W^2\). With weighting="quadrature" a softened numerical weight approximates integration over a hard hole and enters once in both observables. Use the convention the simulator recorded; see WEIGHTINGS.

Regime

For the reference instrument at 260 nm the aperture passes a k half-width of 6.04e4 rad/m against a signal footprint of 3.63e5 rad/m: six times narrower, about 3 % of the solid angle. The corresponding real-space coherence length (~63 µm) exceeds the region that actually emits (~15–30 µm, where \(|p|,|\vartheta|\) stay inside the pulse duration), so this collection is much closer to the coherent limit than to the incoherent one, and a coherent sum sharpens where an incoherent one blurs.

Validity

The collection model is exact given croak’s focal-plane field, but that field comes from a one-dimensional depth quadrature: signal generated at every depth is added with the same transverse phase, i.e. fully coherently in k. A real slab decorrelates depth from depth (diffraction, Gouy phase, walk-through), so the entrance-face model is too coherent, increasingly so the thicker the slab. In our earlier scans at the known truth against the 3-D instrument traces, it improved the residual by 24–30 % at 2–9.5 µm and over-corrected beyond about 20 µm. The leading missing term, the per-depth transverse phase \([k_z(\omega,k_\perp)-k_z(\omega,0)](L-z)\), has a parameter-free model: depth_transverse=True on croak.forward_jax.make_param_trace_fn() moves the depth sum inside the aperture transform with that phase applied per depth node, and evolved arm profiles additionally model incident transverse propagation. In the papers’ current validation protocol the full model with collection retrieves 1.047–1.079 fs (against a 1.053 fs truth) at all 15 thicknesses from 1 to 40 µm through a 0.5 mm hole; see Validation and Modelling the collection aperture. Check radial, azimuthal, aperture and depth convergence for new geometries.

croak.collection.APODISATIONS = ('hard', 'supergauss', 'tanh', 'rcos')

Mask apodisations, ported from ModelPNPS.makemask so a simulated window can be reproduced exactly rather than approximated by a top hat. The reference instrument’s tanh edge is 96.9 um on a 500 um hole – 19 % of the diameter – so the difference is not a detail.

croak.collection.COLLECTION_MODES = ('integrated', 'reimaged')

Which functional of the transformed field the detector measures.

class croak.collection.CollectionAperture(x, y, weight, transmission, z_mask, chromatic=True, mode='integrated', weighting='amplitude')

Bases: object

Quadrature over the collecting aperture, in the mask (far-field) plane.

Parameters:

• x (NDArray[float64])

• y (NDArray[float64])

• weight (NDArray[float64])

• transmission (NDArray[float64])

• z_mask (float)

• chromatic (bool)

• mode (str)

• weighting (str)

x, y

Node positions. Mask-plane metres when chromatic (a physical hole, whose k-space footprint scales with frequency), else absolute transverse wavevectors in rad/m (a frequency-independent k window).

Type:

ndarray, shape (J,)

weight

Quadrature weight including the r dr dphi Jacobian: m^2 when chromatic, else rad^2/m^2. The transmission is NOT folded in – "integrated" needs its square and "reimaged" its first power, so it is kept separate.

Type:

ndarray, shape (J,)

transmission

Amplitude transmission \(W\) at each node (see mask_transmission()).

Type:

ndarray, shape (J,)

z_mask

Distance from the focus to the mask plane (m). Unused when not chromatic.

Type:

float

chromatic

Whether the aperture is a physical hole (True, the usual case) or a fixed window in k. The latter exists because that is the frame a signal-quadrant or Planck k-window lives in, and it is what makes the exact Parseval test possible.

Type:

bool

mode

Which functional of the transformed field the detector measures; see the module docstring.

Type:

{“integrated”, “reimaged”}

weighting

How transmission enters the collected energy: squared (a field filter, the ModelPNPS convention) or once (a quadrature weight over a hard hole, the pnps convention). See WEIGHTINGS. The re-imaged observable carries it once either way.

Type:

{“amplitude”, “quadrature”}

property nodes: int

Number of quadrature nodes.

croak.collection.DEFAULT_N_RADIAL = 6

nodes per radial panel, azimuthal nodes, and the outer radial cut in hole radii. Chosen by the convergence measurement in tests/test_collection.py – on the reference geometry these reproduce a 6x denser quadrature to better than 1e-5 of peak in both collection modes. Named constants rather than repeated literals so mask_hole_aperture() and aperture_from_scan() cannot drift apart.

Type:

Default aperture quadrature

class croak.collection.TransformPhases(chromatic_delay, static_phase, weight_sq, weight_amp)

Bases: object

Precomputed geometry of the focal-plane-to-aperture transform.

The transform matrix is

\[M_{jk}(\omega) = a_k\, \exp\bigl[i(\omega\,\Phi_{jk} - \Theta_{jk})\bigr],\]

with \(\omega\) the absolute optical angular frequency and \(a_k\) the mixture’s area weights, so that \(\tilde S_j = \sum_k M_{jk}\Psi_k\) is the transform of the carrier-removed field at aperture node \(j\). Splitting the phase into a part linear in \(\omega\) and a constant part is what lets a chromatic hole and a fixed k window share one implementation.

Parameters:

• chromatic_delay (NDArray[float64])

• static_phase (NDArray[float64])

• weight_sq (NDArray[float64])

• weight_amp (NDArray[float64])

chromatic_delay

\(\Phi_{jk}\) in seconds: the node’s own arrival offset minus the aperture node’s phase ramp, both of which scale with frequency.

Type:

ndarray, shape (J, K)

static_phase

\(\Theta_{jk}\) in radians, including the carrier-reference term omega0 * p for a chromatic hole. A fixed transverse-wavevector window also contributes its spatial phase.

Type:

ndarray, shape (J, K)

weight_sq

\(\mathrm dk_j W_j^2/(2\pi)^2\) for physical amplitude transmission, or \(\mathrm dk_j W_j/(2\pi)^2\) for quadrature weighting.

Type:

ndarray, shape (J, Nomega)

weight_amp

\(\mathrm dk_j W_j/(2\pi)^2\) — the "reimaged" reduction weight.

Type:

ndarray, shape (J, Nomega)

croak.collection.aperture_from_scan(path, *, mode='integrated', window=None, n_radial=6, n_azimuth=12, pad=2.5)

Build the aperture the reference simulation actually collected through.

Reads the flattened window_def_* scalars a scansave file records (see croak.io.read_simulated_mask_window()) and rebuilds the same hole, resolving a "default" tanh width from the file’s own transverse k-grid.

Parameters:

• path (str) – Path to the simulated scan file.

• mode ({"integrated", "reimaged"}, optional) – Which collection model to build; the two correspond to the file’s Iω_win and Iω_win_reimaged datasets.

• window (str or int, optional) – Which collection hole of a multi-aperture scan (one propagation reduced through several holes, Iω_win_2 …): the trace-window dataset name (_reimaged accepted) or the bare index. None is the first window — the default selection.

• n_radial (optional) – Quadrature controls; see mask_hole_aperture().

• n_azimuth (optional) – Quadrature controls; see mask_hole_aperture().

• pad (optional) – Quadrature controls; see mask_hole_aperture().

Returns:

CollectionAperture

Raises:

KeyError – If the file records no mask window.

Return type:

CollectionAperture

croak.collection.mask_hole_aperture(*, hole_x, hole_y, hole_diameter, z_mask, apod='tanh', apod_param=None, n_radial=6, n_azimuth=12, pad=2.5, mode='integrated', weighting='amplitude')

Build the quadrature for a single apodised collection hole.

The radial quadrature is split at the hole edge into [0, D/2] and [D/2, pad\,D/2]. The integrand is smooth on each panel but has a knee at the edge, so two panels of n_radial nodes resolve an apodised rim that a single panel of 2 n_radial would smear. A compact "rcos" edge of half-width \(\Delta\) is instead split into the flat core [0, D/2-\Delta] and the edge band [D/2-\Delta, D/2+\Delta], outside which it transmits nothing.

Parameters:

• hole_x (float) – Hole centre in the mask plane (m).

• hole_y (float) – Hole centre in the mask plane (m).

• hole_diameter (float) – Hole diameter (m).

• z_mask (float) – Distance from the focus to the mask plane (m).

• apod ({"hard", "supergauss", "tanh", "rcos"}, optional) – Apodisation form; see mask_transmission().

• apod_param (float or None, optional) – Apodisation parameter; see mask_transmission() and resolve_tanh_width().

• n_radial (int, optional) – Quadrature nodes per radial panel and in azimuth, so the total is 2 n_radial n_azimuth (n_radial n_azimuth for a hard edge, which has no outer panel to integrate).

• n_azimuth (int, optional) – Quadrature nodes per radial panel and in azimuth, so the total is 2 n_radial n_azimuth (n_radial n_azimuth for a hard edge, which has no outer panel to integrate).

• pad (float, optional) – Outer radial cut in units of the hole radius. A soft edge leaks well past the nominal rim – the reference tanh window still passes 0.6 % of amplitude at r = D – so the default carries the quadrature to 2.5 times the radius. The "reimaged" mode sums the amplitude rather than its square and so feels that tail ten times more strongly: measured on the 1 fs geometry, dropping to pad = 2 costs it 4e-5 of peak against 2e-7 for "integrated". Ignored (clamped to 1) for a hard edge, which transmits nothing outside the rim.

• mode ({"integrated", "reimaged"}, optional) – Collection model; see the module docstring.

• weighting ({"amplitude", "quadrature"}, optional) – How the edge profile enters the collected energy (WEIGHTINGS). With "quadrature" the edge is built at quadrature_hole_diameter() so the weight encloses exactly the nominal hole area, as the pnps simulator does.

Returns:

CollectionAperture

Return type:

CollectionAperture

croak.collection.mask_transmission(radius, *, hole_diameter, apod='tanh', apod_param=None)

Amplitude transmission of an apodised mask hole at a mask-plane radius.

A direct port of ModelPNPS.makemask; the supported forms are

\[\begin{split}\text{hard} &: \; \mathbb 1[r \le D/2] \\ \text{supergauss} &: \; \exp\bigl[-(2r/D)^n\bigr] \\ \text{tanh} &: \; \tfrac12\bigl[1-\tanh\bigl((r-D/2)/\Delta\bigr)\bigr] \\ \text{rcos} &: \; \tfrac12\bigl[1-\sin\bigl(\tfrac{\pi}{2}\, \mathrm{clip}((r-D/2)/\Delta,-1,1)\bigr)\bigr]\end{split}\]

the last (pnps.beams.chromatic_mask) a raised cosine of compact support: exactly 1 inside \(D/2-\Delta\) and exactly 0 outside \(D/2+\Delta\).

Parameters:

• radius (ndarray) – Distance from the hole centre in the mask plane (m).

• hole_diameter (float) – Hole diameter \(D\) (m).

• apod ({"hard", "supergauss", "tanh", "rcos"}, optional) – Apodisation form.

• apod_param (float or None, optional) – The super-Gaussian exponent \(n\) (default 16), the tanh smoothing width \(\Delta\) or the rcos edge half-width \(\Delta\), both in mask-plane metres. Required for "tanh" and "rcos": their defaults in the simulators depend on the simulation’s own k-grid, so they must be read from the file (or resolved with resolve_tanh_width()) rather than guessed here.

Returns:

ndarray – Amplitude (not intensity) transmission, same shape as radius.

Raises:

ValueError – If apod is unknown, hole_diameter is not positive, or a tanh window is requested without a width.

Return type:

NDArray[float64]

croak.collection.resolve_tanh_width(*, delta_k, z_mask, omega_reference)

Reproduce ModelPNPS’ default tanh smoothing width.

The simulator sets \(\Delta = 3\,\Delta x_{\text{mask}}\) with \(\Delta x_{\text{mask}} = \Delta k\, z_{\text{mask}} c/\omega_0\). Files record the resolved window as the string "default", so it has to be re-derived from the stored k-grid.

Parameters:

• delta_k (float) – Transverse k-grid spacing of the simulation (rad/m).

• z_mask (float) – Mask-plane distance from the focus (m).

• omega_reference (float) – The frequency the simulator evaluated the default at (rad/s): its own grid point nearest \(2\pi c/\lambda_0\), not that value itself. On a coarse spectral grid the two differ, which is why croak.io.MaskWindowSpec resolves it at read time.

Returns:

float – Smoothing width in mask-plane metres.

Return type:

float

croak.collection.transform_phases(mixture, omega, omega0)

Build the transform geometry for a mixture and its collection aperture.

Parameters:

• mixture (FocalMixture) – The focal quadrature. Must carry arms and a collection.

• omega (ndarray, shape (Nomega,)) – Centred angular-frequency grid (rad/s), as the forward model uses.

• omega0 (float) – Carrier angular frequency (rad/s). Both are needed and they are not interchangeable: the aperture’s k footprint scales with the absolute frequency, while the phase croak’s envelope bookkeeping dropped scales with the centred one. See the frequency convention in the module docstring.

Returns:

TransformPhases

Return type:

TransformPhases

Notes

Memory is 16 J K Nomega bytes for the complex matrix the caller builds from this, so the aperture and mixture node counts multiply.

croak.collection.WEIGHTINGS = ('amplitude', 'quadrature')

How the edge profile enters the collected energy. "amplitude" (ModelPNPS/Luna files): the profile is a transmission that filters the FIELD, so the collected energy carries it squared and the re-imaged field once. "quadrature" (pnps files, window_def_*_weighting = "quadrature"): the profile is the integration WEIGHT over a hard hole and enters once in both observables; its edge is shrunk by its own second moment so the weight encloses exactly the nominal hole area (quadrature_hole_diameter()). Use the convention recorded in the file: squaring a tanh edge that is not antisymmetric about the rim under-counts a 0.5 mm aperture by a quarter, and the two recorded observables of such a file describe different apertures.

croak.collection.EDGE_SECOND_MOMENT = {'hard': 0.0, 'rcos': 0.04735763271532445}

Second moment of the antisymmetric edge profiles, \(\int t\,\phi'(t)\,dt\) over the unit-width edge (ported from pnps.beams.EDGE_SECOND_MOMENT); it is what quadrature_hole_diameter() cancels.

croak.collection.aperture_offset(mixture)

Mask-plane offset of the aperture from the phase-matched signal direction.

The transmission-weighted centroid of the aperture nodes minus \(\mathbf r_s=\mathbf r_1+\mathbf r_2-\mathbf r_3\), in metres. Divided by \(z_{\text{mask}}\) and multiplied by \(\omega/c\) it is where the aperture sits in the carrier-removed k frame, so a value of zero means the instrument collects exactly on the phase-matched direction.

Returns:

ndarray, shape (2,) – (dx, dy) in metres.

Raises:

ValueError – If the mixture carries no collection aperture, or the aperture is not chromatic (a fixed k window has no mask-plane position).

Parameters:

mixture (FocalMixture)

Return type:

NDArray[np.float64]

croak.collection.quadrature_hole_diameter(hole_diameter, apod, apod_param)

Edge diameter whose quadrature weight encloses exactly the nominal hole area.

Port of pnps.beams.quadrature_holediam. An antisymmetric edge (\(\phi(t)+\phi(-t)=1\)) is unbiased in one dimension, but in two the ring just outside the nominal radius \(R\) is longer than the ring just inside, so a profile of half-width \(w\) encloses \(\pi R^2 + 4\pi m_2 w^2\) with \(m_2\) from EDGE_SECOND_MOMENT. Centring the edge on \(\sqrt{R^2 - 4 m_2 w^2}\) cancels that term exactly (residual \(O((w/R)^4)\)). This is what a pnps weighting = "quadrature" window does before rendering, so the aperture croak rebuilds must do the same.

Parameters:

• hole_diameter (float) – Nominal hole diameter \(D = 2R\) (m).

• apod ({"hard", "rcos"}) – Edge profile; must have an antisymmetric expansion.

• apod_param (float or None) – Edge half-width \(w\) (m); ignored for "hard".

Returns:

float – The corrected diameter to build the edge at.

Raises:

ValueError – If the profile has no antisymmetric expansion ("supergauss", "tanh"), or the correction would consume the whole aperture.

Return type:

float

Materials

Optical material dispersion: the refractive index and the propagation constant.

The dispersive-propagation model needs a registry of refractive-index data and the derived propagation constant \(\beta(\omega)\).

Each material maps to a callable n(wavelength_m) returning the (real) refractive index. beta() builds the slab propagation constant on a frequency grid, referenced to the group velocity at omega0 so that a pulse propagates without an overall group delay — matching the dispersive forward model.

Materials are resolved in this order: runtime-registered callables (register_material()), then the built-in materials of MATERIALS. The built-in materials come from optmat, a curated library of optical constants: each name uses optmat’s recommended model of that material, whose sources, validity range and reason for the choice optmat records (material_model() returns it). Its indices are absolute (relative to vacuum), as the propagation model needs.

Experimental data of your own can be registered through make_tabulated_material(), which represents the index by a generalised-cross-validation smoothing spline rather than a raw interpolant: the dispersion quantities (beta1(), croak.dispersion.material_gdd(), croak.dispersion.material_tod()) differentiate \(n(\lambda)\), and finite differences of an interpolant through noisy samples amplify that noise catastrophically in the second and third derivatives. The smoothing spline regularises the noise once, at the representation layer, so the finite-difference machinery is unchanged.

croak.materials.refractive_index(material)

Return a callable n(wavelength_m) for material.

Parameters:

material (str) – Material name (case-insensitive): one added via register_material(), or a built-in material (see MATERIALS).

Returns:

callable – Function mapping vacuum wavelength in metres to the real refractive index. A built-in material gives nan at a non-positive or non-finite wavelength and is extrapolated beyond its model’s stated range; a registered material returns what its callable returns. Callers working on a frequency grid should mask non-finite values.

Raises:

ValueError – If material is neither registered nor built in.

Return type:

Callable[[ArrayLike], NDArray[float64]]

croak.materials.beta(material, omega, omega0)

Propagation constant on a frequency grid, referenced to omega0.

Returns \(\beta(\omega) = (\omega+\omega_0)\,n/c - \beta_1(\omega_0)\,\omega\), where omega is the centred grid (relative to omega0) and the group delay term keeps the pulse centred in the time window. Non-finite entries (a non-positive frequency, or outside a registered material’s range) are set to zero.

Parameters:

• material (str or None) – Material name, or None for a dispersionless medium (returns zeros).

• omega (array_like) – Centred angular-frequency grid (rad/s), relative to omega0.

• omega0 (float) – Reference (carrier) angular frequency (rad/s).

Returns:

numpy.ndarray – Real propagation constant (rad/m) of the same length as omega.

Return type:

NDArray[float64]

croak.materials.beta1(material, omega0)

Group delay per unit length \(\beta_1 = d\beta/d\omega\) at omega0.

Computed by numerically differentiating \(\beta(\omega) = \omega\, n(\lambda(\omega))/c\).

Parameters:

• material (str)

• omega0 (float)

Return type:

float

croak.materials.register_material(name, n_func)

Register a custom n(wavelength_m) callable under name.

The callable must accept a vacuum wavelength array in metres and return the real refractive index; values outside its validity range should be returned as nan so beta() masks them. Registered materials are resolved case-insensitively and take precedence over the built-in materials.

Parameters:

• name (str)

• n_func (Callable[[ArrayLike], NDArray[float64]])

Return type:

None

croak.materials.unregister_material(name)

Remove a previously registered custom material (no-op if absent).

Parameters:

name (str)

Return type:

None

croak.materials.make_tabulated_material(lam_m, n)

Build a smooth n(wavelength_m) callable from tabulated samples.

Intended for experimental refractive-index data. The index is represented by a smoothing spline (scipy.interpolate.make_smoothing_spline, with the smoothing penalty chosen automatically by generalised cross-validation) so that the derived dispersion quantities — which differentiate \(n(\lambda)\) up to third order — are not corrupted by measurement noise (calibration errors in the data remain). The fit is performed in log-wavelength, because broadband optical data is typically (near-)logarithmically sampled; this gives uniform fractional resolution across many decades of wavelength.

Parameters:

• lam_m (array_like) – Vacuum wavelength samples in metres (any order; duplicates are dropped).

• n (array_like) – Real refractive index at each wavelength.

Returns:

(n_of_lambda_m, wl_range_m) – n_of_lambda_m(lam_m) returns the smoothed refractive index, with nan outside the sampled range so beta() masks it. wl_range_m is (min_m, max_m).

Raises:

ValueError – If fewer than four valid samples remain after cleaning.

Return type:

tuple[Callable[[ArrayLike], NDArray[float64]], tuple[float, float]]

croak.materials.material_model(material)

Return the optmat model behind a built-in material.

Parameters:

material (str) – A name in MATERIALS (or an earlier name of one), any case.

Returns:

optmat.solid.models.base.Model – optmat’s recommended model of the material, with its name, citation, wavelength_range (m) and identity for provenance.

Raises:

ValueError – If material is not a built-in material.

Return type:

Model

Examples

>>> material_model("BK7").name
'schott_2023/sellmeier'

croak.materials.canonical_name(material)

Return the current name of a built-in material given by an earlier name.

Parameters:

material (str) – A material name, any case.

Returns:

str – The MATERIALS name that material refers to ("SiO2" for "sio2-franta"), or material unchanged if it names no built-in material (a registered one, say).

Return type:

str

Examples

>>> canonical_name("SiO2-Franta")
'SiO2'
>>> canonical_name("MyGlass")
'MyGlass'

croak.materials.read_nk_csv(path)

Read the refractive index from a stacked wl,n / wl,k CSV.

The file format (as exported by some optical-constants databases, e.g. Franta et al.) stacks two blank-line-separated sections: a wl,n refractive-index table followed by a wl,k extinction-coefficient table, each with a wavelength column in micrometres. Only the real refractive index n is read; the extinction k is ignored (the propagation model is lossless).

Parameters:

path (pathlib.Path or str) – Path to the CSV file.

Returns:

(lam_m, n) – Vacuum wavelength in metres and the real refractive index, in file order (ascending wavelength).

Raises:

ValueError – If the file contains no wl,n section.

Return type:

tuple[NDArray[float64], NDArray[float64]]

croak.materials.MATERIALS: tuple[str, ...] = ('SiO2', 'BK7', 'CaF2', 'BaF2', 'MgF2', 'Al2O3', 'Si')

Names of the built-in materials, each optmat’s recommended model of it.

croak.materials.DESCRIPTIONS: dict[str, str] = {'Al2O3': 'sapphire (ordinary ray)', 'BK7': 'SCHOTT N-BK7', 'BaF2': 'BaF₂', 'CaF2': 'CaF₂', 'MgF2': 'MgF₂ (ordinary ray)', 'Si': 'crystalline silicon (opaque below 1.1 µm)', 'SiO2': 'fused silica'}

A short description of each built-in material, by name (for labels and help).

Gases

Gas refractive index at a pressure and temperature, from optmat.

A gas medium is dispersive like a solid (see croak.materials), but its refractive index depends on the gas density, and hence on pressure and temperature, as well as on wavelength. This module provides gas_refractive_index(), which builds an n(wavelength_m) callable for a gas at a chosen pressure and temperature. That callable has the same shape as croak.materials.refractive_index()’s result, so a gas plugs into the whole dispersion pipeline by registering it with croak.materials.register_material() and treating the optical path length as the “thickness”.

Physics

The dispersion of each gas is optmat’s recommended model of its per-molecule polarisability \(\alpha(\lambda)\) (optmat, optmat.gas.models.registry.default_model); air is optmat’s dry air, the mole-fraction sum of its constituents’ models. The index at a number density \(N\) follows from the Lorentz–Lorenz relation,

\[\frac{n^2 - 1}{n^2 + 2} = \frac{N\alpha}{3\varepsilon_0},\]

and the density from the real gas, \(N = P/(Z k_\mathrm{B}T)\), with the compressibility factor \(Z\) of each gas’s reference equation of state (CoolProp, through optmat; air is Lemmon et al.’s pseudo-pure dry air, and ozone has none, so it is taken as ideal). \(n - 1\) is then nearly proportional to \(P/T\), but not exactly: at 20 °C the real gas differs from the ideal by 0.05 % for helium and air at 1 bar, 0.5 % at 10 bar, and up to 15 % for SF₆ at 10 bar. A state where the gas would condense (SF₆ above about 21 bar at 20 °C, for example) is refused.

croak.gases.GASES: tuple[str, ...] = ('Air', 'He', 'Ne', 'Ar', 'Kr', 'Xe', 'H2', 'D2', 'N2', 'O2', 'CH4', 'SF6', 'N2O', 'CO2', 'O3', 'CO')

air, then optmat’s species.

Type:

Names of the gases available (case-insensitive)

croak.gases.ROOM_TEMPERATURE_K: float = 293.15

Default gas temperature (K) when none is given (ca. 20 °C).

croak.gases.canonical_name(gas)

Return the current name of a gas given by any case or an earlier name.

Parameters:

gas (str) – A gas name.

Returns:

str – The GASES name gas refers to ("He" for "Helium").

Raises:

ValueError – If gas is not a known gas.

Return type:

str

Examples

>>> canonical_name("helium")
'He'

croak.gases.gas_model(gas)

Return the optmat model of a gas’s polarisability.

Parameters:

gas (str) – A name in GASES (or an earlier name of one), any case.

Returns:

optmat.gas.models.base.Model or optmat.gas.air.AirModel – optmat’s recommended model of the gas, or its dry-air model, with its name, wavelength_range (m) and identity for provenance.

Raises:

ValueError – If gas is not a known gas.

Return type:

Model | AirModel

Examples

>>> gas_model("Helium").name
'fit/He'

croak.gases.gas_refractive_index(gas, pressure_bar, temperature_k=293.15)

Build an n(wavelength_m) callable for a gas at given pressure/temperature.

The index follows from optmat’s polarisability model of the gas through the Lorentz–Lorenz relation at the real-gas density \(P/(Z k_\mathrm{B}T)\) (see the module docstring). The returned callable has the same contract as croak.materials.refractive_index()’s result, so it can be registered with croak.materials.register_material() and used wherever a material is.

Parameters:

• gas (str) – Gas name (case-insensitive); one of GASES.

• pressure_bar (float) – Gas pressure in bar, >= 0. 0.0 gives vacuum (n = 1).

• temperature_k (float, optional) – Gas temperature in kelvin, > 0 (default ROOM_TEMPERATURE_K, ca. 20 °C).

Returns:

callable – Function mapping vacuum wavelength in metres to the real refractive index, nan at a non-positive or non-finite wavelength (masked by croak.materials.beta()) and extrapolated beyond the model’s stated range.

Raises:

ValueError – If gas is not a known gas, the pressure is negative or not finite, the temperature is not positive and finite, or the gas is not a gas at that state (condensed, or beyond its equation of state’s range; optmat.exceptions.StateError, a ValueError).

Return type:

Callable[[ArrayLike], NDArray[float64]]

Examples

>>> n = gas_refractive_index("Air", 1.0, 273.15)
>>> float(n(800e-9)) - 1.0
0.000286...

Maths utilities

Numerical helpers.

This module provides the small set of mathematical utilities the retrieval core needs to be self-contained: pulse envelope shapes (gauss(), sech()), distribution statistics (fwhm(), moment()), a numerical derivative, and the angular-frequency/wavelength conversion wlfreq().

All quantities are SI unless stated otherwise.

croak.maths.wlfreq(x)

Convert between angular frequency and wavelength.

Implements \(2\pi c / x\), which maps angular frequency \(\omega\) (rad/s) to vacuum wavelength \(\lambda\) (m) and vice versa. The function is its own inverse.

Parameters:

x (array_like) – Angular frequency (rad/s) or wavelength (m).

Returns:

numpy.ndarray – Wavelength (m) or angular frequency (rad/s).

Return type:

NDArray[float64]

croak.maths.fwhm_to_sigma(fwhm, power=2.0)

Convert a full width at half maximum to a Gaussian sigma.

For a (hyper-)Gaussian \(\exp[-\tfrac12 (x/\sigma)^{p}]\), the FWHM and width parameter are related by \(\sigma = \mathrm{FWHM} / (2\,(2\ln 2)^{1/p})\).

Parameters:

• fwhm (float) – Full width at half maximum.

• power (float, optional) – Exponent p of the (hyper-)Gaussian. 2 (default) is an ordinary Gaussian.

Return type:

float

croak.maths.gauss(x, sigma=None, *, fwhm=None, x0=0.0, power=2.0)

Evaluate a (hyper-)Gaussian.

Computes \(\exp[-\tfrac12 ((x-x_0)/\sigma)^{p}]\). Exactly one of sigma or fwhm must be given.

Parameters:

• x (array_like) – Sample points.

• sigma (float, optional) – Standard-deviation-like width parameter.

• fwhm (float, optional) – Full width at half maximum (converted to sigma internally).

• x0 (float, optional) – Centre of the distribution.

• power (float, optional) – Exponent of the (hyper-)Gaussian. 2 (default) is an ordinary Gaussian.

Return type:

NDArray[float64]

croak.maths.sech(x, *, fwhm=None, tau=None, x0=0.0)

Evaluate a hyperbolic-secant amplitude profile.

Computes \(\operatorname{sech}((x-x_0)/\tau)\). Exactly one of tau or fwhm must be given. The intensity \(\operatorname{sech}^2\) has full width at half maximum \(\mathrm{FWHM} = 2\,\tau\,\operatorname{arcsech} (1/\sqrt2) = 2\,\tau\,\ln(1+\sqrt2)\), so tau = fwhm / (2 ln(1+sqrt 2)).

Parameters:

• x (array_like) – Sample points.

• fwhm (float, optional) – Intensity FWHM.

• tau (float, optional) – Characteristic width.

• x0 (float, optional) – Centre of the distribution.

Return type:

NDArray[float64]

croak.maths.fwhm(x, y, *, level=0.5)

Estimate the full width at level of the maximum by linear interpolation.

Finds the outermost crossings of level * max(y) on either side of the peak and returns their separation. Assumes a single dominant peak.

Parameters:

• x (array_like) – Monotonically increasing coordinate.

• y (array_like) – Signal values.

• level (float, optional) – Fraction of the maximum at which the width is measured (0.5 for the conventional FWHM).

Returns:

float – The full width. Returns nan if the threshold is not crossed on both sides of the peak.

Return type:

float

croak.maths.moment(x, y, n=1)

Return the n-th moment of a distribution.

Computes \(\sum_i x_i^n\, y_i / \sum_i y_i\), the weighted average of \(x^n\) with weights y. The first moment (n=1) is the centroid.

Parameters:

• x (array_like) – Coordinate values.

• y (array_like) – Non-negative weights (e.g. an intensity distribution).

• n (int, optional) – Moment order (default 1).

Return type:

float

croak.maths.derivative(f, x, order=1, *, dx=None)

Numerical derivative of a scalar function via central finite differences.

Used for low-order dispersion quantities (e.g. the group delay \(\beta_1 = d\beta/d\omega\)). For order 1 and 2 this uses the standard second-order-accurate central stencils.

Parameters:

• f (callable) – Scalar function of a scalar argument.

• x (float) – Point at which to evaluate the derivative.

• order (int, optional) – Derivative order (1 or 2).

• dx (float, optional) – Step size. Defaults to a scale-aware value based on x.

Returns:

float – The estimated derivative.

Return type:

float

croak.maths.planck_taper(x, left0, left1, right1, right0)

Planck-taper window defined by four corner points.

The window is 0 for x <= left0, rises smoothly to 1 across [left0, left1], stays 1 over [left1, right1], falls back to 0 across [right1, right0], and is 0 for x >= right0. The smooth transition is the Planck-taper of arXiv:1003.2939 eq. (7).

Parameters:

• x (array_like) – Sample points.

• left0 (float) – Corner positions (must satisfy left0 <= left1 <= right1 <= right0).

• left1 (float) – Corner positions (must satisfy left0 <= left1 <= right1 <= right0).

• right1 (float) – Corner positions (must satisfy left0 <= left1 <= right1 <= right0).

• right0 (float) – Corner positions (must satisfy left0 <= left1 <= right1 <= right0).

Returns:

numpy.ndarray – Window values in [0, 1].

Return type:

NDArray[float64]

Constants

Physical constants (SI units).

Values follow the CODATA recommended set.

croak.constants.C: float = 299792458.0

Speed of light in vacuum (m/s).