API Reference

High-fidelity forward modelling of PNPS (Parametrized Nonlinear Process Spectrum) pulse-characterisation traces by full spatially-resolved nonlinear propagation, using Luna.jl.

ModelPNPS is a standalone package for generating synthetic pulse-characterisation traces directly from the underlying experimental physics — spatial beam overlap, mode shape, mask edges, chromatic vignetting, material dispersion, phase-matching and full χ⁽ⁿ⁾ nonlinear propagation. Given an analytic input pulse and an experimental geometry, it produces the trace a real apparatus would record. These ground-truth traces are intended for testing advanced retrieval algorithms and developing new characterisation techniques.

The long-term ambition is a complete PNPS trace modeller spanning the full Geib et al. (2019) taxonomy of nonlinear process × parametrization (FROG, d-scan, time-domain ptychography, …). The currently implemented process is TG-FROG (Transient-Grating FROG); see the documentation roadmap for the planned methods.

Physical model

The TG-FROG signal in the boxcar geometry is a degenerate four-wave mixing process

    k_signal = k_g2 - k_g1 + k_test

Three input beams (gates g1, g2 and a delayed test pulse t) are crossed inside a thin nonlinear substrate (e.g. UV fused silica). Their interference produces a transient grating; the test pulse diffracts off this grating into the fourth corner of the boxcar. Scanning the test-pulse delay τ and spectrally resolving the diffracted signal yields a 2-D I(ω, τ) spectrogram.

Two beam models

Two AbstractInputBeam subtypes are provided:

• HE11Beam — the master experimental model. The HE₁₁ mode of a hollow capillary fibre is collimated by a long lens, clipped by a four-hole apodised mask in the collimated beam, then focused into the substrate. Each hole selects one of the four boxcar arms.
• GaussianBeam — a simplified Gaussian-beam model that places three Gaussian beams directly at the correct k-space angles (no mask, no fibre mode). Useful as a sanity-check baseline.

Three signal-extraction window types

Three AbstractSignalWindow subtypes are provided:

• PhysicalMaskWindow — the master experimental signal extraction: a frequency-dependent apodised hole in the mask plane (chromatic vignetting is captured exactly).
• PlanckWindow — a smooth, frequency-independent radial Planck taper in k-space (no chromatic vignetting; baseline for the Gaussian model).
• PlanckOmegaWindow — a smooth, frequency-dependent Planck taper that mimics the chromatic vignetting of PhysicalMaskWindow while keeping the smooth-edge advantage. Used to isolate the two effects (smooth-edge vs ω-scaling) within the Gaussian model.

Two field representations

By default the propagation is Luna's complex envelope (Grid.EnvGrid): an analytic field about a carrier. build_setup(field_mode=true) instead propagates the real, carrier-resolved field on a Grid.RealGrid, which has no envelope/carrier split, no dropped third-harmonic term and no negative-frequency wrap. That matters when the pulse is only a cycle or two long — at 260 nm a 1 fs pulse is 1.15 optical cycles — where the envelope approximation is marginal by construction and the two representations can be compared directly. It costs roughly twice the memory and three times the time per delay point, so it is a diagnostic, not the production default. See build_setup's field_mode, response and ffac keywords.

High-level usage

using ModelPNPS
import Luna.Scans

beam   = HE11Beam(125e-6, 5.0, 0.1)
window = PhysicalMaskWindow(
    holex=-0.75e-3, holey=-0.75e-3,
    holediam=0.5e-3, zmask=0.1)

setup = build_setup(; λ0=260e-9, τfwhm=2e-15, energy=0.2e-6,
                      thickness=10e-6, material=:SiO2,
                      mask_diam=1.0e-3, mask_spacing=0.5e-3,
                      beam, window)

τ = collect(range(-10e-15, 10e-15, 80))
exec = Scans.SlurmExec(@__FILE__, length(τ); memory="18G",
                       arraymode=:batch)
run_scan(setup, τ; scan_name="my_trace", exec)

The full simulation requires SLURM (it is hours of CPU per delay scan with typical grid sizes); the unit tests exercise everything except the actual Luna.run call by passing skip_propagation=true to simulate_delay_point.

AbstractInputBeam

Abstract supertype for the input-beam models from which the three TG-FROG beamlets are constructed. Concrete subtypes: HE11Beam, GaussianBeam.

AbstractSignalWindow

Abstract supertype for k-space windows used to extract the FWM signal beam from the propagated field. Concrete subtypes: PhysicalMaskWindow, PlanckWindow, PlanckOmegaWindow.

Return the same array for every requested propagation slice.

Callable wrapper for an already-built setup.

FrozenRamanPolarEnv(t, r)

Envelope Raman polarisation response with a frozen response kernel.

Wraps Luna.Nonlinear.RamanPolarEnv, precomputing the frequency-domain response function hω once at construction (at unit density). Luna's own response recomputes the time-domain kernel and its FFT on every call, which is negligible for modal simulations (one call per step) but dominant for free-space grids, where the response runs once per transverse point — ~10⁶ calls per RK stage on a 1024² grid, each re-evaluating the 13-mode Hollenbeck–Cantrell sum over the doubled time grid. Freezing is exact here: the density is constant (densityfun = z -> 1) and the :intermediate response ignores its density argument entirely.

The per-call convolution below reproduces Luna.Nonlinear.(::RamanPolar) line for line, minus the kernel update and with preallocated plan applications; the test suite verifies agreement with Luna's response to machine precision, which guards against drift in Luna's internals.

Fused error metric for SignalQuadrantNorm: the DP5 error estimate is computed element-by-element on the fly from the stepper's stage arrays instead of materialising a field-sized yerr array (Luna.RK45 allocates it lazily only for norms without a fused version). Same per-element expression and accumulation order as the materialised path, so the result is bit-identical.

GaussianBeam(w0, f_foc)
GaussianBeam(; w0, f_foc)

A Gaussian beam with 1/e² intensity radius w0 at the focus. f_foc is retained only so the crossing angle (and Δk) can be derived from the same mask geometry as the HE₁₁ model.

Fields

• w0::Float64 — 1/e² intensity radius at the focus [m]
• f_foc::Float64 — focusing-lens focal length [m] (geometry only)

HE11Beam(a, f_coll, f_foc)
HE11Beam(; a, f_coll, f_foc)

The HE₁₁ capillary mode imaged from the fibre output through a collimating lens (f_coll) onto a beam mask, then focused (f_foc) into the substrate. The Hankel transform of the mode has a closed form

    Ẽ(k_⊥) ∝ -a² u₁₁ J₁(u₁₁) J₀(a k_⊥) / (a² k_⊥² - u₁₁²)

where u₁₁ is the first zero of J₁. The "image" of the fibre core at the substrate has demagnified radius a_scaled = a · f_foc / f_coll.

Fields

• a::Float64 — fibre core radius [m]
• f_coll::Float64 — collimating-lens focal length [m]
• f_foc::Float64 — focusing-lens focal length [m]

InputPulseData(ω, Eω)

A measured or simulated input pulse: a complex spectrum Eω on an ABSOLUTE, ascending, (approximately) uniform angular-frequency axis ω [rad/s]. Inject it through build_setup's input_pulse keyword — HE₁₁ beam model only, because there the 1-D reference spectrum IS the pulse and the chromatic mask vignetting is applied to it downstream, so an arbitrary field composes exactly (see build_beamlets). Amplitude units are irrelevant: the beamlet builder rescales the assembled beam to the requested energy.

Companion utilities, typically chained in this order: load_input_pulse → spectral_window! → center_pulse! → build_setup(...; input_pulse=p).

Return a view of a propagation slice held in memory.

Read one propagation slice from an output backend.

PhysicalMaskWindow(holex, holey, holediam, zmask, apod, apod_param)
PhysicalMaskWindow(; holex, holey, holediam, zmask,
                   apod = :supergauss, apod_param = nothing)

A frequency-dependent mask hole: physical position (holex, holey) and diameter holediam in the mask plane, sitting zmask (= focal length) upstream of the substrate. The mask plane ↔ k-space mapping is

    (x_mask, y_mask) = (k_x, k_y) · zmask · c / ω

so the same physical hole transmits a wavelength-dependent k-space region (chromatic vignetting). Apodisation choices:

• :hard — binary (1 inside the hole, 0 outside).
• :supergauss — exp(-(2 r/d)^n) with n = apod_param (default 16).
• :tanh — smooth 0.5(1 - tanh((r - d/2)/Δ)) with Δ = apod_param in mask-plane metres (default 3 × Δx_mask evaluated at the carrier wavelength).

Fields

• holex, holey — hole centre in the mask plane [m]
• holediam — hole diameter [m]
• zmask — focal length / mask-to-focus distance [m]
• apod — apodisation type :hard | :supergauss | :tanh
• apod_param — apodisation parameter (nothing → defaults)

PlanckOmegaWindow(xc, yc, holediam, f_foc, pad)
PlanckOmegaWindow(; xc, yc, holediam, f_foc, pad = 1.25)

A frequency-dependent Planck-taper window. The hole is specified in the mask plane by its centre (xc, yc) and diameter holediam; at frequency ω the window centre and half-width in k-space are

    k_c(ω)    = (ω/c) · (xc, yc) / f_foc
    k_hole(ω) = (ω/c) · (holediam/2) / f_foc

This restores the chromatic vignetting of PhysicalMaskWindow while keeping the smooth-edge advantage of PlanckWindow.

Fields

• xc, yc — hole centre in the mask plane [m]
• holediam — hole diameter in the mask plane [m]
• f_foc — focusing-lens focal length [m]
• pad — outer roll-off multiplier (typically 1.25)

PlanckWindow(kxc, kyc, kwidth, pad)
PlanckWindow(; kxc, kyc, kwidth, pad = 1.25)

A radial Planck-taper window centred at (kxc, kyc) in k-space with flat half-width kwidth and an outer roll-off radius pad·kwidth. The window is frequency-independent: the same mask shape is applied to every spectral component, so chromatic vignetting is removed.

Fields

• kxc, kyc — k-space centre of the window [rad/m]
• kwidth — flat half-width of the window [rad/m]
• pad — multiplier setting the outer roll-off (typically 1.25)

Callable wrapper that resolves and builds setup arguments lazily.

Callable RK45 error norm built by signal_quadrant_norm.

Holds the signal-quadrant mask, the relative floor, and — on a device — a cached 0/1 indicator of the quadrant on the solver's own array type.

TGFROGSetup

Container holding everything that is built once (independent of the FROG delay τ): grids, propagation operators, FFT plan, the three pre-built input beamlets, the signal window(s) and the metadata dictionary.

Use build_setup to construct one and simulate_delay_point or run_scan to use it.

Fields

The struct is a passive bundle; fields are not part of the public API and may evolve. Use the constructors and methods provided.

TraceExtractOutput(setup, zvec, arraytype)

A Luna output handler that reduces each saved z-slice to the trace spectra immediately and keeps only the results, so the full field is never stored, streamed or transferred.

Satisfies the Output interface Luna.run uses (the save call, willsave, metadata calls, and the generic check_cache fallback). Metadata is discarded: ModelPNPS builds its own from setup.combined_grid, and the temp file this replaces was thrown away too.

Return the squared deviation of abs(value) from mean_amplitude.

_beamlet_profile(grid, xygrid, Eωk, holex, holey, zmask; nr, rmax, nθ=64)
    -> (r, Eωr, asym)

The spatially resolved complex focal field of one beamlet, reduced to a radial profile Eωr[ω, r] about its own centre, with the radius axis r in metres and a per-ω measure asym of how well the radial reduction describes it.

Diagnostic only: nothing here feeds the propagation.

Where the beamlet actually is

In this representation the beamlets do not sit at BOXCARS corners in real space — they all cross at the focus, and the corners are in k-space. build_he11_kspace builds the field in (ω, ky, kx) with the transverse amplitude the Hankel transform of the HE₁₁ mode, so k-space is the COLLIMATED (mask) plane — makemask maps x = kx·zmask·c/ω — and real space, after ifft over dims 2 and 3, is the FOCAL plane. A hole at mask position (holex, holey) therefore selects k around k₀ = hole·ω/(c·zmask), and in the focal plane that offset is a tilt, not a displacement: measured on the production geometry, the gate beamlet peaks at the real-space grid centre to the pixel, and carries a phase slope of 2.402e5 rad/m against the predicted k₀ = 2.417e5.

So the centre of this profile is the grid centre. Centring it on the mask-hole position mapped through the focus — 1 mm out, against a 26 µm spot — would sample nothing.

Why the tilt is removed first

k₀·r reaches 37 rad across the default sampling radius, so an azimuthal average of the raw complex field would annihilate it. The field is demodulated by exp(-i k₀(ω)·r) before sampling, leaving the beamlet's own envelope. k₀ ∝ ω, so the coefficient is a constant hole/(c·zmask) in s/m; it is stored, and multiplying the profile by exp(+iω(cₓx + c_yy)) restores the full field. Note that the removed tilt is physical — a linear delay across the beamlet, i.e. the pulse-front tilt of the crossing geometry — not an artefact.

Accuracy

Sampled by bilinear interpolation on nr × nθ polar points (the focal spot is ~27 grid cells across, so the interpolation is not the limiting error) and averaged azimuthally. asym is the azimuthal RMS of |E| over its mean, restricted to radii carrying signal: 1.2–3.8 % on the production geometry, i.e. the radial reduction is a good description but not an exact one, and a consumer can see how good.

Integrating |Eωr|² with the 2πr dr Jacobian and dividing by the cell area reproduces the stored Iω_beamlet to ~1.5 % at the default nr = 64 — the shortfall is the Airy wings beyond rmax plus that asymmetry, and it converges (0.980/0.987/0.988 at nr = 128 against 0.975/0.984/0.986 at 64, at 200/260/350 nm). The test suite asserts this closure, which is the check that the centre, the Jacobian and the normalisation are all right.

_build_setup_resolved(setup_args) -> TGFROGSetup

Build the setup, resolving arraytype FIRST and then calling build_setup through Base.invokelatest.

Luna.resolve_arraytype(:cuda) loads the GPU package at run time, and methods defined by a package loaded during a call are not visible to that same call — Julia rejects them as "too new to be called from this world context". Resolving first and invoking afterwards puts the construction in a world where the array type's constructors exist.

This is why a scan script should pass arraytype=:cuda inside setup_args and let this happen on the compute node, rather than loading the GPU package itself.

_completed_scanidcs(scan_name) -> Set{Int}

Scan indices already present in <scan_name>_collected.h5, i.e. those whose trace data is not all zero. An empty set if the file does not exist yet.

Reads one point at a time: the file may be large and this runs before any propagation.

Return no additional values for Output.scansave.

_envelope_intensity(grid, Et)

|A(t)|² — the ENVELOPE intensity — from whatever time-domain field grid produces. On an EnvGrid that is Et itself; on a RealGrid the field is carrier-resolved and the envelope is recovered through its analytic signal.

Both conventions coincide numerically: Luna builds a real-grid pulse as √I·cos(ω₀t) and an envelope-grid pulse as √I·exp(iΔωt), so |A|² = I either way. Keeping the envelope intensity in the output metadata means a consumer of a field-mode file sees the same physical quantity in It/Ito as in every envelope file, rather than a carrier-modulated one it would have to demodulate.

_extract_slice_device!(Iint, Ireim, Ifull, Ez, wsgn, quadrng)

Reduce one (Nω, Nky, Nkx) device slice into the three per-ω spectra, writing into column views of the host result arrays. Three reductions over dims (2, 3) and one small copy back per slice; the field itself never moves.

Mathematically identical to extract_signal_spectra plus _quadrant_spectrum!, which is what a slice arriving on the host uses instead. The sums are formed in a different order here, so results agree to rounding rather than bitwise — the standard everywhere else on the device path.

Both operands of the two windowed reductions have the SAME shape: mapreduce over several arrays does not broadcast (Base throws DimensionMismatch, and a GPU backend may silently compute something else), which is why the re-imaging sign lives inside wsgn rather than in a (1, Nky, Nkx) array of its own. The quadrant sum is a single-array reduction over a strided view, so it reads only the quadrant instead of masking the whole field.

_field_responses(response, χ3) -> Tuple

The nonlinear response for a field-resolved (RealGrid) run.

• :nothg (and :auto) — (3/4) ε₀ χ³ |E_a|² E via Luna.Nonlinear.Kerr_field_nothg. This is the SAME physics content as the envelope Kerr_env, evaluated on a carrier-resolved field, so an envelope-versus-field comparison made with it isolates representation error with nothing else changed. It is the default for exactly that reason.
• :thg — ε₀ χ³ E³ via Luna.Nonlinear.Kerr_field, which adds the third-harmonic and counter-rotating terms the envelope drops. The difference between the two runs is precisely what the envelope omits.

Return the HDF5 dataset stored at key, or throw for a malformed file.

Return the HDF5 group stored at key, or throw for a malformed file.

Whether grid is field-resolved (real) rather than an envelope grid.

Return whether an HDF5 key names a scan result dataset.

Provenance label for the error norm used by a scan.

_plan_1d(grid)

Forward transform plan for the 1-D reference pulse: complex for an EnvGrid, real-to-complex for a RealGrid. Fields.GaussField dispatches its time-domain shape on the grid type but takes the plan from the caller, so the two have to be chosen together.

_profile_meta(r, Eωr, asym, coef, rmax_req, holex, holey, zmask,
              rmax_units, which) -> Dict

Package _beamlet_profile's output for the output file. Complex data is split into two real datasets, matching the Eω_beamlet_re/_im convention (h5py reads HDF5.jl's native complex compound awkwardly), and enough geometry is recorded for the file to be self-describing without the script that made it.

_quadrant_spectrum!(out, Ez, quad)

Accumulate |E|² over the transverse points selected by the (Nky, Nkx) mask quad into the length-Nω vector out, without materialising any field-sized temporary (the previous broadcast allocated two per z-slice).

Read key when it exists, otherwise return nothing.

_reduce_slice!(o::TraceExtractOutput, Ez, iz)

Reduce one saved slice into column iz, routing on where the slice actually is.

A device propagation delivers every slice on the device — including z = 0, which is the step start rather than an endpoint and so comes through the interpolant, because Luna.needs_host_save declines the copy HostOutput would otherwise make.

The host branch therefore serves a host propagation (extract_on_save=true on the CPU), and any save that is genuinely interpolated. It goes through the original extract_signal_spectra/_quadrant_spectrum! kernels against the plain host window the setup already holds: no extra memory, no transfer, and the host result is bit-identical to the save-the-stack route by construction rather than by a parallel implementation that could drift.

_resolve_zsave(zsave, zmax) -> Vector{Float64}

Resolve the zsave propagation-snapshot specification into a validated, sorted vector of z positions [m] at which the field is saved during propagation.

• zsave::Integer — a uniform grid of zsave points over [0, zmax] (range(0, zmax, zsave)), reproducing the legacy nz behaviour exactly (including the entrance slice at z=0 and the exit slice at z=zmax).
• zsave::AbstractVector — explicit material thicknesses [m]. Must be strictly increasing, all >= 0, and all <= zmax. zmax is appended if not already present (within rtol=1e-12) so the full-thickness (":end") slice always exists.

Because the propagation is a forward-marching integrator with z-independent dynamics, the field saved at an intermediate z is identical to a dedicated run of thickness z, so a single zmax run yields every shorter thickness for free.

The function is idempotent: re-resolving an already-resolved vector (which the integer path produces with an entrance slice at z=0) returns it unchanged, so it is safe to call more than once on the same grid.

Canonical name of the response actually used, for the output metadata.

_scan_peak(dset) -> Float64

Largest absolute value over every computed delay point of a collected trace dataset. Read one point at a time rather than whole: this runs against a file a scan may still be writing, and the datasets grow with the delay count.

Shift selected envelope-spectrum dimensions into natural order.

Shift an envelope spectrum into natural order; field spectra are already ordered.

_signed_window(w, arraytype, Nω) -> array

The signal window with the re-imaging sign pattern (-1)^((iky-1)+(ikx-1)) folded in, on arraytype.

Folding the sign into the window lets one array serve both reductions of extract_signal_spectra: the signed sum needs it, and the intensity sum is unaffected because |±w·E|² == |w·E|². Carrying a separate sign array instead would cost a second field-sized device array, and could not be broadcast into the reduction anyway (see the shape note in _sqn_fused).

A 2-D window is expanded to 3-D: the reduction takes mapreduce over two arrays, which does not broadcast shapes.

The quadrant indicator as a (1, Nky, Nkx) array on y's array type, built once and cached on the norm. Broadcasting it into the reduction costs one small array (8 MB even at the largest campaign shape) and keeps the reduction over whole, contiguous arrays.

The mask is built from the same BitMatrix the host path uses, via quadrant_ranges — so its rectangle assertion still guards the device path.

Device version of _sqn_fused: the same six sums, computed as reductions along ω into one partial per transverse point, which are then split by quadrant.

Every operand of a reduction here has the SAME shape. That is deliberate: mapreduce over several arrays does not broadcast shapes (Base throws DimensionMismatch, and a GPU backend may quietly compute something else instead), so the (1, Nky, Nkx) quadrant mask cannot be folded into the reduction and is applied afterwards, to the small per-transverse-point partials. Reducing over strided views of the quadrant would also work in principle, but whole-array reductions are the shape the rest of Luna's device code uses and the one best supported across backends.

Both halves are summed directly rather than one being total - other, so no cancellation is involved. The error estimate is never materialised — it is formed inside the reduction kernel.

_to_time(grid, Eω)

The 1-D time-domain field for a spectrum on grid's own frequency axis: an inverse FFT for an EnvGrid (whose spectrum is FFT-ordered about the carrier) and an inverse real FFT for a RealGrid (whose spectrum is the monotonic rfft half-spectrum of a real field).

_trace_results(setup, o::TraceExtractOutput) -> NamedTuple

The same NamedTuple the save-the-stack route returns, assembled from an extraction handler. Kept next to that route's assembly block so the two cannot drift.

a_scaled = fibre core radius imaged onto the focal plane.

apply_delay(Eωk, grid, τ) -> Array{ComplexF64,3}

Apply a time delay τ (seconds) to a frequency-domain field by multiplying each spectral component by exp(-i ω τ). τ = 0 returns a copy equal to the input.

apply_tilt(Eωxy, xygrid, Δkx, Δky) -> Array{ComplexF64,3}

Multiply a real-space field E(ω, y, x) by the phase ramp exp(i Δkx · x) · exp(i Δky · y), which shifts its centre by (Δky, Δkx) in k-space (after FFT). Δkx = Δky = 0 is the identity.

build_beamlets(beam, grid, xygrid, geom, Eω, energy, energyfun_ω;
               apod=:supergauss, apod_param=nothing, ϕ=nothing,
               profile=true, profile_nr=64, profile_rmax_units=6)
    -> (Eωk_g1, Eωk_g2, Eωk_t_base, Iω_beamlet, beam_metadata::Dict)

Construct the input beamlets at the substrate, in k-space. The geometry geom is a NamedTuple(mask_diam, mask_spacing, f_foc, λ0, τfwhm, geometry) shared by both beam models. geom.geometry is :tg for the three-beam boxcar layout (g1, g2, t-base), or :sd for the two-beam self-diffraction layout — HE11Beam only, and rejected for any other beam by build_setup — which returns nothing in place of Eωk_g2 (a zero array of that size is half a gigabyte of pure waste) and puts the probe in Eωk_g1 and the delayed gate in Eωk_t_base.

ϕ is accepted for a uniform interface across beam models and ignored by both: the spectral phase is already carried by the 1-D reference Eω that build_setup passes in. profile, profile_nr and profile_rmax_units control the diagnostic radial focal profile added to beam_metadata; see _beamlet_profile.

For HE11Beam: builds the full HE₁₁ k-space field, rescales to the requested energy, then applies three apodised hole masks (g1, g2, t). Each beamlet sits at one of the boxcar corners. Iω_beamlet is the spatially-integrated spectrum of g1 (used as a chromatic-vignetting diagnostic in the output file).

For GaussianBeam: builds a Gaussian-Gaussian field with energy energy/3 per beam, ifft's to real space, then applies real-space tilts to position the three beams at the boxcar corners. Iω_beamlet here is just the (unvignetted) input spectrum scaled to energy/3; it is returned for uniformity with the HE₁₁ model so downstream code never special-cases the beam type.

build_gaussian_kspace(grid, xygrid, beam::GaussianBeam,
                      λ0, τfwhm, energy) -> Array{ComplexF64,3}

Construct the 3-D field E(ω, ky, kx) for a Gaussian-Gaussian spatio-temporal pulse: temporal Gaussian envelope (FWHM = τfwhm) at carrier λ0, spatial Gaussian (1/e² radius = beam.w0) centred on the grid, with total spectral energy normalised to energy. Internally uses Luna.Fields.GaussGaussField and Luna.setup (with no nonlinearity) to construct the field, then discards the throw-away transform/FT.

build_he11_kspace(grid, xygrid, beam::HE11Beam, Eω) -> Array{ComplexF64,3}

Construct the 3-D field E(ω, ky, kx) for the HE₁₁ capillary mode imaged onto the focal plane, multiplied by the 1-D spectral pulse Eω. Phase ramps shift the beam from the FFTW corner to the centre of the spatial grid.

The closed-form Hankel transform of the J₀ mode profile is used; the a²k² - u₁₁² denominator is finite at all (kx, ky) sample points for reasonable grid sizes (the singular ring is at radius u₁₁/a, well outside typical Nyquist limits at the focal-plane scale).

build_setup(; λ0, τfwhm, energy, thickness, material,
              mask_diam, mask_spacing, beam, window,
              kwargs...) -> TGFROGSetup

Build the once-per-simulation setup: temporal/spatial grids, propagation operators, FFT plans, the three input beamlets and the signal window(s). The defaults reproduce the master script context/tgfrog_DUV_mask_apod6.jl.

Required keyword arguments

• λ0, τfwhm, energy — pulse carrier wavelength [m], intensity FWHM [s], total pulse energy [J]
• thickness, material — substrate thickness [m] and Luna PhysData material symbol (e.g. :SiO2)
• mask_diam, mask_spacing — mask hole diameter [m] and edge-to-edge gap [m]
• beam::AbstractInputBeam — input-beam model (HE11Beam or GaussianBeam)
• window — signal-extraction window: a single AbstractSignalWindow or a vector of them (the latter is used by the Gaussian example to save both the ω-independent and ω-dependent windows in one run)

Optional keyword arguments

• trange = 40e-15 — temporal window [s]
• λlims = (160e-9, 500e-9) — wavelength window [m]
• R, N — spatial half-width [m] and grid size; if either is nothing, both are computed via optimal_spatial_grid
• apod, apod_param — apodisation for the input-beamlet masks (only relevant for HE11Beam)
• geometry = :tg — beam layout. :tg is the four-hole boxcar TG-FROG geometry (three inputs, signal in the fourth corner); :sd places two collinear holes for self-diffraction, whose 2k_E - k_G signal sits one slot further out on the same axis. It selects both the beamlet layout and the k-space bound used by optimal_spatial_grid. :sd is implemented for HE11Beam only; any other value throws an ArgumentError
• fftsize = :pow2 — how the temporal sample count is rounded up: :pow2 to the next power of two, :smooth to the next even 2,3,5-smooth size (a smaller grid for the same resolution). Envelope mode only — Grid.RealGrid has no such control
• GDD = 0.0, TOD = 0.0 — group-delay and third-order dispersion [s², s³] applied to the input pulse
• input_pulse = nothing — an InputPulseData: use this measured/simulated complex spectrum as the source instead of the analytic Gaussian (HE11Beam only). λ0/τfwhm then serve only as nominal values (mask apodisation defaults, diagnostics, metadata); energy still sets the beam energy (the data's amplitude scale is irrelevant); GDD/TOD compose on top if nonzero. See load_input_pulse, spectral_window!, center_pulse!
• raman = false — include the delayed (Raman) part of the nonlinear response via FrozenRamanPolarEnv; requires a material with an :intermediate Raman model in Luna.PhysData.raman_parameters (for :SiO2 the multimode Hollenbeck–Cantrell response). The total polarisation is (3/4)ε₀χ³[(1-f_R)|E|²E + f_R E(h_R⊛|E|²)] — equal prefactors on both terms, the envelope-defined f_R convention of Luna's prop_gnlse, so the quasi-static limit reproduces the Kerr-only response exactly
• raman_fraction = 0.18 — envelope-defined nuclear fraction f_R of χ³ (the Blow–Wood silica value)
• raman_impl = :batched — Raman implementation: :batched computes the convolution for all transverse points at once (two batched FFTs per RHS evaluation); :frozen is the legacy per-column FrozenRamanPolarEnv. Results agree to rounding accuracy
• field_mode = false — propagate the real, carrier-resolved field on a Luna.Grid.RealGrid instead of the complex envelope on an EnvGrid. There is then no carrier/envelope split, no dropped third-harmonic term and no negative-frequency wrap; the cost is roughly 2× the memory and 3× the time per delay point (measured 3.0× at N = 64 and 3.3× at N = 128, at matched step counts). The envelope path is untouched and remains the default
• response = :auto — field-mode nonlinearity: :nothg (= :auto) for (3/4) ε₀ χ³ |E_a|² E, the same physics content as the envelope Kerr_env and hence the response for an envelope-versus-field comparison; :thg for ε₀ χ³ E³, which adds what the envelope drops. Ignored unless field_mode = true
• ffac = 6 — field-mode nonlinear-grid sampling factor, forwarded to Grid.RealGrid. 6 (the default) sizes the fine grid for E³; 4 is enough for :nothg alone and typically removes the oversampling entirely, halving memory and per-step cost. It changes the grid, so use it only with a convergence check against the default
• raman — not implemented in field mode (see the error message there for why)
• beamlet_profile = true — store the gate beamlet's spatially resolved complex focal field as a radial profile Eω_beamlet_r_re/_im (Nω, nr) plus the radius axis beamlet_r in metres, so the pulse that actually drives the signal can be computed rather than assumed. Diagnostic only — no propagation result depends on it — costing one 2-D inverse transform per ω ONCE at setup and ~130 kB in the file. See _beamlet_profile for where the beamlet is (the focus, not a BOXCARS corner) and why the geometric tilt is removed first
• beamlet_profile_nr = 64 — radial samples. Measured radial closure against Iω_beamlet on the production geometry, at 200 / 260 / 350 nm: 0.955/0.972/0.980 at nr = 32, 0.975/0.984/0.986 at 64, 0.980/0.987/0.988 at 128. 64 is where it has essentially converged, for 262 kB at Nω = 256; the residual ~1.5 % is truncation at rmax plus the beamlet's real azimuthal asymmetry
• beamlet_profile_rmax_units = 6 — outer radius, in units of λ0·f_foc/mask_diam (w0 for GaussianBeam)
• factored_linop = true — use Luna's lazy (factored) linear operator and normalisation, saving two field-sized arrays; bit-identical to the materialised versions
• frozen_transverse = false — ABLATION, not physics: build the linear operator with k_z(ω, k⊥) replaced by k_z(ω, 0), so every k⊥ component gets the same ω-dependent phase and the transverse field pattern (beamlet profiles, crossing interference, tilt phases) is frozen exactly at its entrance-face form, while temporal dispersion, the nonlinearity, the apodisation and the k-space collection all run unchanged. Note the pulse-front tilts live in the initial condition and remain. Recorded in the output metadata as frozen_transverse (absent = 0 = normal propagation, which is what every pre-existing file is)
• store_window = true — store the materialised window array(s) in the output metadata (≈1 GiB at production size); the window parameters (window_def) are always stored and reconstruct the array via build_window
• arraytype = Array — array type the propagation runs on. Pass :cuda to build the beamlets, operators and window on the GPU; it is resolved lazily, so a scan script passes it inside setup_args rather than as a run_scan keyword and the GPU package is then loaded on the compute node, never on the submitting host
• beamlets_on_host = false — on a device run, keep the pre-built beamlets in host memory and upload the delayed sum once per delay point — two fewer resident device fields in exchange for one transfer per point. Use it when the card is memory-bound; see memory_budget
• optimal_grid_kwargs — extra kwargs forwarded to optimal_spatial_grid
• extra_grid_metadata — additional entries merged into the output combined_grid dict

build_window(w::AbstractSignalWindow, grid, xygrid; λ0=nothing)
    -> Array{Float64, N}

Materialise the precomputed signal-extraction window. Returns a (Nky, Nkx) 2-D array for PlanckWindow and a (Nω, Nky, Nkx) 3-D array for PhysicalMaskWindow and PlanckOmegaWindow. λ0 is forwarded to makemask for default :tanh apodisation widths only.

center_pulse!(p::InputPulseData; oversample=8) -> (p, tshift)

Remove the linear spectral-phase component so the temporal intensity envelope peaks at the data FFT's natural origin (array index 1), returning the applied shift tshift [s] (positive = the pulse arrived late and was advanced). A pure linear phase is physically irrelevant; numerically, centring minimises the trange the simulation needs to hold the pulse plus the delay scan, and — more importantly — it is what makes the spectrum interpolatable: a pulse far from its grid's natural time origin has a spectral phase rotating by up to π per sample, which no Re/Im interpolation can resample (interp_input_pulse warns if it sees this). interp_input_pulse then re-anchors the interpolated field at t = 0, the middle sample of Luna's centred target time grid. Requires an (approximately) uniform ω grid. The returned shift is reported modulo the data grid's time period (the on-grid phase is identical for any branch).

delayed_input(setup, τ) -> Array{ComplexF64,3}

Coherent input field for scan delay τ, in the gate-delay convention: the stored trace $T(ω, τ)$ has the GATE pair delayed by +τ relative to the probe. Physically the probe arm carries the delay stage, so the probe is delayed by -τ, which equals gating at +τ up to a global time shift that the time-integrating measurement cannot see. The same-τ gate pair stays untouched, so the geometrical smearing structure of the crossed-beam layout (a same-τ gate pair) is preserved. Files written with this convention carry /grid/delay_convention = "gate" and need NO delay-axis reversal on loading (retrieval loaders can detect the marker; legacy marker-less files are reversed as before).

extract_signal_spectra(Eωk, window_array, xygrid)
    -> (Iω_integrated, Iω_reimaged)

Apply a precomputed signal window to a propagated field and extract two spectral diagnostics:

1. Iω_integrated — |E|² summed over all (ky, kx). Models a spectrometer collecting all the signal light.
2. Iω_reimaged — |E|² at the centre pixel of the IFFT'd field. Models a spectrometer fed only by the on-axis re-collimated signal.

Eωk is either a single (Nω, Nky, Nkx) slice, for which both spectra are length-Nω vectors, or the (Nω, Nky, Nkx, Nz) stack Luna.run produces, for which both are (Nω, Nz). The 4-D method loops over z and calls the 3-D one, so its peak extra memory is one windowed slice rather than a second copy of the whole stack (tens of GB at production size).

window_array is broadcast over ω (if 2-D) or matched directly (if 3-D), and over the Nz z-slices in either case.

interp_input_pulse(grid, p::InputPulseData) -> Vector{ComplexF64}

The pulse's complex spectrum on grid.ω (the grid's ABSOLUTE frequency axis), zero outside the data's range. Real and imaginary parts are interpolated separately with cubic B-splines, which is accurate when the data grid is finer than the simulation grid — a warning is emitted if it is not (then spectral detail is being invented between samples; supply denser data instead). The input is expected to have been moved to its data FFT's natural origin with center_pulse!. After interpolation, the field is shifted to the middle sample of Luna's centred target time grid, matching Luna's native Fields.DataField convention.

load_input_pulse(path; ω_key="ω", Eω_key="Eω") -> InputPulseData

Read an InputPulseData from an HDF5 file: an absolute angular frequency axis under ω_key and a complex spectrum under Eω_key (a native complex dataset, e.g. as written by HDF5.jl or h5py).

load_simulated_scan(filename; window_key="Iω_win", z_index=:end,
                    z_thickness=nothing) -> NamedTuple

Read the raw HDF5 file produced by run_scan and return its contents as a NamedTuple, with all ω-dependent arrays fftshifted into natural (centred) order and the requested z slice(s) extracted from the propagated trace.

Arguments

• filename: path to the <scan_name>_collected.h5 file.

Keyword arguments

• window_key="Iω_win": which scansave dataset to use as the FROG trace. Common choices:
  • "Iω_win" — full-beam k-space integrated spectrum
  • "Iω_win_reimaged" — on-axis re-imaged spectrum
  • "Iω_win_ωdep" — ω-dependent window (Gaussian two-window setup)
  • "Iω_win_ωdep_reimaged" — ω-dependent re-imaged
• z_index=:end: which propagation z slice to use; the default :end picks the final (full-propagation) slice. Pass an Int for a specific slice index, or :all to return every z slice as a (Nω, nz, Nτ) stack (the equivalent of the trace at every saved material thickness).
• z_thickness=nothing: select the slice whose saved z position [m] is nearest this material thickness. Requires /grid/zsave in the file (written by recent run_scan runs); takes precedence over z_index.

Returned NamedTuple

The optional zsave, It_beamlet, Ito_beamlet, To, and Ito fields are returned only when their corresponding datasets are present.

To inspect the full signal-beam collection (and hence the exact collection / chromatic-vignetting efficiency Iω_win ./ Iω_full), load the signal-quadrant reference with window_key="Iω_full".

makemask(holex, holey, holediam, grid, xygrid;
         zmask, apod=:supergauss, apod_param=nothing,
         λ0_for_default=nothing) -> Array{Float64,3}

Build a 3-D (Nω, Nky, Nkx) apodised-hole mask. For each (ω, ky, kx) sample, the k-vector is mapped to the mask-plane position x = kx · zmask · c / ω (and likewise for y), and a hole of diameter holediam centred at (holex, holey) is evaluated.

λ0_for_default is only used when apod=:tanh and apod_param===nothing, in which case the smoothing width is set to 3·Δx_mask evaluated at the carrier wavelength.

memory_budget(setup_args::NamedTuple) -> NamedTuple

Resident device memory one delay point of setup_args will need, and the host peak build_setup will reach, broken down by buffer. setup_args is the same NamedTuple run_scan and verify_against_collected take; only the grid-determining entries are read, and building the 1-D time grid is the whole cost, so this is free to call.

This exists because guessing is expensive. The envelope path obeys a simple rule — 9 RK45 registers plus one transform buffer, i.e. 10× the field size, measured exactly on an A40 — and the field path does not: its state is twice as long in ω, its nonlinear evaluation runs on a grid twice as long again in time, and the no-THG response carries a complex analytic-signal buffer on that grid. At the 40 µm production shape (N = 768) that is 92 GiB against the envelope's 24. Finding this out by running is an hour of rented GPU and a dead process.

Fields: Nω, Nt, Nto, Nωo, field (one state array), the per-buffer terms state, et_win, eto, ewo, pto, analytic, window, input, and the totals device and host. All in GiB.

state … window are what the transform and solver hold; input is the per-delay field delayed_input produces. Only the first group is allocated by build_setup, so a measurement taken across build_setup alone will fall short of device by state, analytic and input — those appear when the first delay point runs.

The buffer set and its aliasing are NonlinearRHS.TransFree's: Pωo always aliases Eωo (the inverse transform consumes it), Pto aliases Eto when every response is pointwise (the envelope Kerr, and field :thg, but not field :nothg), and Et_win exists only when the grid is oversampled. window is the extraction window, which is device-resident when save-time extraction is used — the default on a device.

optimal_spatial_grid(f, mask_diam, mask_spacing, λmin, λmax;
                     n_airy=5, pts_per_lobe=10, safety=1.5,
                     margin=1.1, geometry=:tg) -> (R, N)

Return (R, N) for a Luna FreeGrid(R, N) chosen so that the spatial grid

1. contains at least n_airy Airy diffraction patterns of the longest wavelength λmax from a mask hole of diameter mask_diam focused by a lens of focal length f (real-space containment), and
2. resolves the Airy pattern at the shortest wavelength λmin with at least pts_per_lobe points across the central lobe (real-space resolution), and
3. has a k-space half-extent that comfortably encloses the FWM nonlinear k-vectors generated at λmin from the outermost mask hole, with a safety headroom factor (k-space containment).

N is rounded up to the next power of 2 for FFT efficiency. Diagnostic information is printed via @info.

Arguments

• f: focal length of the focusing lens [m].
• mask_diam: diameter of each mask hole [m].
• mask_spacing: edge-to-edge spacing between adjacent mask holes [m].
• λmin, λmax: shortest and longest wavelengths the simulation must represent [m]. These should bracket the input spectrum and its FWM products.

Keyword arguments

• n_airy=5: number of Airy patterns the grid should contain at λmax.
• pts_per_lobe=10: real-space samples across the central Airy lobe at λmin.
• safety=1.5: multiplier on the required nonlinear k-vector envelope to guard against aliasing.
• margin=1.1: multiplier on the resolved grid size before rounding up to the next even 2,3,5-smooth FFT size (guards the containment against grid quantisation).
• geometry=:tg: the beam layout the k-space bound (3.) is computed for. :tg is the four-hole boxcar, whose χ⁽³⁾ combinations reach three times the hole offset; :sd is the two-hole self-diffraction layout, whose 2k₁ - k₂ signal sits one further slot out along the same axis. See the comment on x_max in the implementation for the two bounds.

quadrant_ranges(sig_quad) -> (ys, xs)

The signal quadrant as a pair of index ranges. In FFT ordering the negative half of each k axis is exactly the second half of its index range, so the mask is a dense rectangle — which lets the device norm use strided views instead of a boolean mask (a BitMatrix cannot enter a device kernel, and a masked reduction would need a gather).

Throws if the mask is not that rectangle, so a future change to the k-space layout cannot silently corrupt the solver's error control.

run_scan(setup, τs; scan_name, exec, kwargs...) -> Nothing

Build a Luna.Scans.Scan over the delay array τs and run simulate_delay_point at every τ, calling Output.scansave to write each result into the collected HDF5 file at "<scan_name>_collected.h5". The metadata block (combined_grid) is written once on the first scan point.

exec must be a Luna.Scans.AbstractExec instance (e.g. Scans.SlurmExec(...) or Scans.LocalExec()).

zsave selects the propagation snapshots saved at every delay (see simulate_delay_point): an Integer gives a uniform grid of that many points over [0, thickness] (default nz), or a Vector of explicit material thicknesses [m] (e.g. [1e-6, 10e-6, 20e-6, 40e-6]). thickness is appended to the vector if absent so the final slice is always the full-propagation output. The trace datasets become (Nω, nz, Nτ) and the realized z positions are stored once in /grid/zsave. Because the field at an intermediate z equals a dedicated thickness-z run, every shorter thickness comes free from one full-thickness run; note that peak memory scales with the number of z points.

extra_outputs(output_namedtuple) is an optional callable returning extra named tuples to splat into scansave. The default is empty.

Keywords

• scan_name: base name of the collected file, "<scan_name>_collected.h5".
• exec: the Luna.Scans.AbstractExec instance described above.
• nz = 2, zsave = nz: propagation snapshots, as above.
• init_dz = 5e-7, rtol = 1e-6, max_dz = 0.0: solver settings forwarded to simulate_delay_point; max_dz = 0.0 means thickness/2. They are recorded in the file's /grid block as provenance.
• norm = Luna.RK45.weaknorm: RK45 error norm.
• norm_builder = nothing: a callable setup -> norm, used instead of norm, for a norm that cannot exist before the setup does. Pass norm_builder = signal_quadrant_norm to get signal_quadrant_norm built lazily on the compute node.
• twin_period = 1: accepted steps between applications of the spectral/temporal windows. 1 applies them after every step, which makes the apodisation damping scale with the step count; a large value applies them only at saves, which with step_on sit at identical positions for any rtol.
• fftw_threads = 0: FFTW threads per process, set where the plans are created so that it reaches procs workers (a top-level Luna.set_fftw_threads does not). With procs workers sharing cpus cores, pass cpus ÷ procs. 0 leaves it alone.
• fftw_mode = :estimate: FFTW planning effort, set on the same path and for the same reason. MEASURE-class planning of production-size 3-D transforms costs tens of minutes per worker.
• stream = true: write the propagation slices to a node-local temp file rather than holding the whole (ω, ky, kx, z) stack in memory (~2.15 GB per slice at production size). Ignored when save-time extraction is active, which stores no slices at all.
• extract_on_save = nothing: reduce each slice as it is produced; see simulate_delay_point. nothing picks the per-device default.
• skip_existing = false: resume an interrupted scan by skipping delay points already present in the collected file. An all-zero slice is the "not yet computed" marker, the same test verify_against_collected uses.

Eager variant: wrap an already-built setup (costs nothing extra when the setup exists anyway, e.g. in interactive use or LocalExec runs).

run_scan(setup_args::NamedTuple, τs; kwargs...)

RECOMMENDED for scan scripts: pass the build_setup keyword arguments as a NamedTuple, e.g.

setup_args = (; λ0, τfwhm, energy, thickness, material,
                mask_diam, mask_spacing, λlims, beam, window,
                R=366.0e-6, N=1024)
run_scan(setup_args, τ; ...)

The setup is then built lazily on each process that executes scan points. This form is robust under EVERY execution mode, including multi-worker (procs > 0) queue scans: a NamedTuple of parameters serialises to the workers by value, whereas a NAMED function defined in a script (make_setup() = ...) serialises by reference and fails to deserialise on workers (Julia ships code only for anonymous closures). The wrapping closure here is defined inside ModelPNPS, which Luna loads on the workers.

signal_quadrant_norm(setup::TGFROGSetup; floor_rel=1e-6)

Region-relative RK45 error norm for weak-signal accuracy at moderate rtol.

The default Luna.RK45.weaknorm measures the step error relative to the norm of the WHOLE field, which the three pump beamlets dominate. The FWM signal is orders of magnitude weaker, so the stepper's error budget — concentrated on the fastest-evolving (signal) components — allows a per-step signal error of order rtol × ‖pump‖/‖signal‖ relative to the signal: with rtol = 1e-6 the collected signal carries measured solver errors of 0.1–1% mid-slab, growing to ~10% at 40 µm (see 90_solver_accuracy_test.jl). Brute-forcing rtol = 1e-8 fixes this at ~4× the step count.

This norm instead measures relative error separately in the signal k-space quadrant (kx < 0, ky < 0 — the same quadrant Iω_full integrates) and in the rest of the field, and returns the larger: rtol then controls the signal's OWN relative error directly, recovering weak-signal accuracy at close to the default-rtol step count.

While the signal quadrant is still (nearly) empty its error is measured against a floor of floor_rel × ‖rest‖ (never below atol), so early steps are not throttled by a 0/0 relative error; once the signal exceeds that fraction of the pump field the relative control takes over.

Pass the result to simulate_delay_point / run_scan via their norm keyword. Validate a new (rtol, floor_rel) choice against a tight-rtol reference before production use (the pass criterion used here: every z-slice of Iω_win within 1e-3 relative of an rtol = 1e-8 run).

simulate_delay_point(setup::TGFROGSetup, τi;
                     nz=2, zsave=nz, init_dz=5e-7, rtol=1e-6, max_dz=0.0,
                     norm=Luna.RK45.weaknorm, twin_period=1,
                     filename=nothing, extract_on_save=nothing,
                     skip_propagation=false)
    -> NamedTuple

Run the full per-delay computation: apply delay τi to the test beam, coherently superpose the three beamlets, propagate them through the substrate via Luna.run, apply each signal window and extract two spectra per window. The returned NamedTuple has, for a single window, fields (Iω_win, Iω_win_reimaged, Iω_full). For a vector of windows the suffixes recorded in setup.window_suffix are appended (e.g. Iω_win_ωdep, Iω_win_ωdep_reimaged), and the single Iω_full is shared. All extracted arrays have shape (Nω, nz). The returned NamedTuple also carries zsave, the vector of realized z save positions [m] (length nz) — this is metadata, not a per-delay trace, and is excluded from the scansave dataset splat by run_scan.

zsave selects the propagation snapshots. Pass an Integer for a uniform grid of that many points over [0, zmax] (default nz), or a Vector of explicit material thicknesses [m] (e.g. [1e-6, 10e-6, 20e-6, 40e-6]); zmax is appended to the vector if absent. Because the field at an intermediate z equals a dedicated thickness-z run, every shorter thickness comes free from one zmax run. Peak memory scales with nz (the in-memory 4-D field is held per slice).

Iω_full is the signal beam collected in full: |E|² integrated over the signal's k-space quadrant only. The propagated field holds the three strong pump beamlets (at the g1/g2/test boxcar corners) plus the weak FWM signal at the fourth corner; integrating over all of k-space would be dominated by the pumps, so we restrict to the quadrant the signal occupies (kx<0, ky<0), which captures the whole signal lobe without aperture vignetting while excluding the pumps. Iω_win ./ Iω_full is therefore the exact per-(ω, τ) collection / chromatic-vignetting efficiency of the signal aperture, so the trace can be corrected for collection vignetting exactly rather than via a power-law approximation. (This assumes the boxcar beams are well separated, so pump tails leaking into the signal quadrant are negligible vs. the signal.)

The solver keywords go straight to Luna.run: init_dz [m] is the first step, rtol the RK45 relative tolerance, max_dz [m] the step ceiling (0.0 means zmax/2), and norm the error norm — pass signal_quadrant_norm to control the weak signal's OWN relative error rather than the pump-dominated whole field's. twin_period is the number of accepted steps between applications of the spectral/temporal windows; the default 1 applies them on every step, and larger values change the result at the apodisation-leakage level.

extract_on_save reduces each z-slice to its spectra as it is produced, so the field is never stored, streamed or transferred. It defaults to true for a device propagation, where the saved stack costs 14–18 % of the delay point in temp-file traffic, and to false on the host; the two routes are bit-identical, so passing true on the host is safe.

Pass filename to persist the propagation to disk: the Luna.run then writes to an Output.HDF5Output at that path instead of an in-memory Output.MemoryOutput. Downstream extraction is identical either way (both outputs index as output["Eω"]/output["z"]); filename is ignored when skip_propagation=true.

Setting skip_propagation=true substitutes the input field for the Luna output, exercising every other code path. This is used by the unit tests to keep the suite fast and deterministic.

spectral_window!(p::InputPulseData, λmin, λmax;
                 wfrac_blue=0.05, wfrac_red=0.03) -> p

Smooth tanh band-pass in place: unity well inside (λmin, λmax) [m], rolling off with tanh edges of width wfrac_red · ω(λmax) on the red side and wfrac_blue · ω(λmin) on the blue side. Use before injection to remove content the simulation band does not (or should not) carry — e.g. a residual driver remnant that survived an imperfect spectral filter, which would otherwise dominate the χ³ interaction. The windowed field is the ground truth the retrieval is compared against, so keep the applied window with the run's provenance.

verify_against_collected(setup_args, collected, scanidcs;
                         zsave, init_dz=5e-7, rtol=1e-6, max_dz=0.0,
                         norm=Luna.RK45.weaknorm, twin_period=1,
                         stream=true, extract_on_save=nothing)
    -> Vector{Dict}

Recompute selected delay points of an existing scan and compare against the collected HDF5 file — the A/B harness for validating a new code path (or a changed grid) against reference data.

For each scan index in scanidcs, the delay τ is read from /scanvariables/τ in collected, the point is recomputed via simulate_delay_point with the given solver settings (pass the SAME settings the reference scan used, unless deliberately testing a change), and every returned trace dataset (Iω_win, Iω_full, ...) present in the file is compared. Reference points that are still all-zero (not yet computed by a running scan) are reported as NaN and skipped.

Returns one Dict per point with the delay, wall time, Sys.maxrss() [GiB], and for each dataset ks the global relative difference maximum(abs, new - ref)/maximum(abs, ref), plus three diagnostics:

Both normalisations matter. A delay-scan wing carries a signal orders of magnitude below the τ≈0 signal, so a difference that is irrelevant in the assembled trace can still be a large fraction of that point's own peak. relscan is what a FROG retrieval sees; the own-peak number is the stricter statement about the code path.

To test a grid change (e.g. N=640 against an N=1024 reference), pass the changed N inside setup_args — differences then reflect the grid, not the code.