qscat.tuning

The automatic discretisation tuner: derive a minimal FEM-DVR-ECS grid at a target precision from the potential and energy range, instead of hand-picking element lengths. See docs/physics/discretisation-tuning.md, including the documented limits of the 1-D convergence probes.

GridCost	grid_cost's report: one 1-D grid's discretisation cost.
IncidentSpec	TD Gaussian-wavepacket placement + observation boundary.
PotentialProfile	The potential-analysis output consumed by the mesh/ECS generators.
ProbeResult	A convergence probe's verdict: the computed value, whether it is converged (change under one refinement / vs. a fine reference is < rtol), the discretisation cost (grid point count) it was computed at, and a detail dict with the raw numbers behind the verdict (refined value, relative error, reference value, ...) for diagnostics/reporting.
Refine2dReport	refine_to_2d_convergence's third return value.
Refine2dStep	One ADOPTED refinement step recorded by refine_to_2d_convergence.
RefinementCoordinate
TensorCost	tensor_cost's report: the tensor-product (2-D) problem's size, plus ROUGH anchored estimates of its sparse-LU cost.
analyze_potential	Analyze V on [x_min, x_max], returning its PotentialProfile.
equidistribution_elements	Return real-region element lengths equidistributing de Broglie phase.
grid_cost	A single grid's discretisation cost: just its DVR point count.
interaction_region	The (R_lo, R_hi) bracket over which V_int(r, R) is still transitioning between its low-R and high-R regimes.
max_stable_angle	Largest ECS rotation angle (deg) for which V stays bounded on the tail.
optimal_real_mesh	Sweep DVR orders, returning the (mesh, order) with fewest points.
order_for_wavenumber	Smallest DVR order in orders resolving wavenumber k at a fixed real element length element_len, in points-per-wavelength.
probe_channel_representation	Is the energy-normalized channel function F_{k,l} REPRESENTABLE on grid's real region?
probe_electronic	Convergence of the lowest anion bound electronic-state energy at internuclear distance R, under one h-refinement of elec_grid.
probe_nuclear	Convergence of the n_vib lowest vibrational eigenvalues under one h-refinement of nuclear_grid.
propose_grid	The one-shot a-priori FemDvrEcsGrid for model/coordinate over energy_range = (E_min, E_max).
refine	One h-refinement step: split every REAL element into two half-length elements; leave the ECS tail, quadrature order, and x_min untouched.
refine_elements_in_window	Span-preservingly subdivide every real element overlapping [R_lo, R_hi] until each piece is <= target_len; elements entirely outside the window are returned untouched.
refine_to_2d_convergence	Iteratively refine whichever of g_r (electronic) / g_R (nuclear) grid is under-resolved for observable, until it stops moving.
required_extent	Free-function form of IncidentSpec.required_extent, for callers who prefer a function over a bound method.
resonance_curve_arrays	Efficient adiabatic resonance curve (R, V_d(R), Gamma(R)).
tensor_cost	ROUGH, ANCHORED cost ESTIMATES for the tensor-product 2-D problem.
tune_ecs_tail	Exp-growth ECS-tail element lengths absorbing wavenumber K at angle.
tw_analysis	Auto-place the TD Gaussian wavepacket for energy_range = (E_min, E_max) -- the "TW analysis".

Automatic FEM-DVR-ECS discretisation tuner.

Public API:

• analyze_potential – pure potential analysis (no models): samples a plain callable V(x) on a dense real grid and returns the local- wavenumber / forbidden-region-decay profile, classical turning points, and boundary singularities (e.g. the -1/r origin). This is the sole input the mesh/ECS generators consume.

• PotentialProfile – the frozen dataclass analyze_potential returns.

• equidistribution_elements – adaptive real-region element lengths equidistributing de Broglie phase per element (capped by kappa-decay length in classically forbidden runs, refined near turning points / singularities).

• optimal_real_mesh – h/p sweep over DVR orders, returning the (mesh, order) combination with the fewest DVR points.

• order_for_wavenumber – the smallest DVR order (from a fixed candidate set) resolving a given wavenumber at a fixed element length, in points-per-wavelength – see qscat.tuning.mesh.

• refine_elements_in_window – span-preservingly subdivide every real element overlapping an [R_lo, R_hi] window until each piece is below a target length; elements outside the window are untouched – the LOCAL, min_len-overriding refinement the resonance-aware nuclear mesh uses to super-refine a narrow resonance crossing – see qscat.tuning.mesh.

• max_stable_angle – the largest ECS rotation angle (capped at the double-ECS bound, ~35 deg) for which a potential V stays bounded on the rotated tail contour.

• tune_ecs_tail – exp-growth ECS-tail element lengths sized to absorb a given outgoing wavenumber K down to a target decay.

• ProbeResult, refine, probe_nuclear, probe_electronic, probe_channel_representation – the decoupled 1-D convergence probes: empirical validators that tell the tuner whether a candidate grid resolves the physics (nuclear vibrational levels, the electronic bound-state energy, and the cheap/diagnostic channel-representation check that catches an unresolved fast outgoing wave) – see qscat.tuning.probes.

• grid_cost, tensor_cost, GridCost, TensorCost – the cost model: exact DVR point counts, plus ROUGH anchored estimates (nnz, factor memory/time) of the tensor-product problem’s sparse-LU cost, for RELATIVE ranking of candidate grids. The two TypedDict`s are the fixed shape of the reports, for annotating a report a caller keeps or passes on – see `qscat.tuning.metrics.

• propose_grid – the one-shot a-priori grid assembler: model adapter -> analyze_potential -> optimal_real_mesh -> max_stable_angle + tune_ecs_tail -> a complete FemDvrEcsGrid, the a-priori half of the hybrid tuner – see qscat.tuning.propose.

• IncidentSpec, required_extent, tw_analysis – incident-state / test-function placement: the TD Gaussian-wavepacket spec, the real- region extent it forces, and the (best-effort) Tannor-Weeks auto-tune that places it for a target energy range – see qscat.tuning.incident.

• interaction_region – the R-window where the electron-molecule interaction V_int(r, R) is non-negligible, from model.v_int alone – see qscat.tuning.resonance.

• resonance_curve_arrays – the efficient adiabatic resonance-curve sampler (R, V_d(R), Gamma(R)): dense inside interaction_region, a single far point at the asymptote – see qscat.tuning.resonance.

• refine_to_2d_convergence – the general, model-agnostic 2-D-convergence FALLBACK: given any scalar observable closing over a real cross-section (or a synthetic test function), iteratively refine`s whichever of the electronic/nuclear grids moves the observable more, until the larger relative move is under `rtol or max_iter adopted steps is hit – the step-6 spot-check’s supervised loop, generalized – see qscat.tuning.refine2d.

• Refine2dReport, Refine2dStep, RefinementCoordinate – the fixed shape of that loop’s report dict and of one adopted step within it, and the two coordinate names a step can carry.

class qscat.tuning.GridCost

grid_cost’s report: one 1-D grid’s discretisation cost.

• n_points: the grid’s DVR point count (FemDvrEcsGrid.n). EXACT – the number of unknowns a 1-D problem on this grid carries, and the quantity a grid-vs-deck size comparison is made on.

class qscat.tuning.IncidentSpec(position, impulse, sigma, observation=None)

TD Gaussian-wavepacket placement + observation boundary.

For the TI route this dataclass is not used – see the module docstring.

Parameters:

• position (float)

• impulse (float)

• sigma (float)

• observation (float | None)

position

The wavepacket centre r0 (bohr) at t=0.

Type:

float

impulse

The mean momentum p0 (bohr^-1); negative launches the packet INWARD (toward the interaction region), matching the n2_2d_td_cross_section convention (wp_in’s p0 < 0 – see projects/n2_2d_td_cross_section/convergence.py).

Type:

float

sigma

The Gaussian width (bohr).

Type:

float

observation

The flux-extraction / test-function boundary (bohr); None if the caller has not set one (required_extent then falls back to just the wavepacket’s own tail).

Type:

float or None

required_extent()

How far the REAL region must reach to CONTAIN this wavepacket and its observation boundary: max(position + N*sigma, observation), N = _TAIL_SIGMAS.

A METHOD on the spec (propose_grid calls it directly on its typed incident: IncidentSpec | None parameter); the module-level required_extent(spec) below remains for callers who prefer a free function.

Return type:

float

incident_energy()

The wavepacket’s mean kinetic energy, impulse**2 / 2 (mass 1, the electronic-coordinate convention every IncidentSpec is defined against).

propose_grid reads this to make sure the RESOLUTION (not just the extent) of the a-priori mesh covers the incident’s local wavenumber – see its docstring for why required_extent alone is not enough whenever a hand-built IncidentSpec implies an energy above energy_range’s own E_max.

Return type:

float

class qscat.tuning.PotentialProfile(x, V, k, kappa, turning_points, singularities)

The potential-analysis output consumed by the mesh/ECS generators.

Parameters:

• x (NDArray[float64])

• V (NDArray[float64])

• k (NDArray[float64])

• kappa (NDArray[float64])

• turning_points (NDArray[float64])

• singularities (NDArray[float64])

class qscat.tuning.ProbeResult(value, converged, cost, detail)

A convergence probe’s verdict: the computed value, whether it is converged (change under one refinement / vs. a fine reference is < rtol), the discretisation cost (grid point count) it was computed at, and a detail dict with the raw numbers behind the verdict (refined value, relative error, reference value, …) for diagnostics/reporting.

value and detail stay deliberately broad, unlike the fixed-shape tuner reports in qscat.tuning.metrics/refine2d. One ProbeResult type is shared by three probes that measure different things, so neither field has one shape: value is a vibrational-energy ARRAY from probe_nuclear but a scalar energy from probe_electronic and a scalar norm from probe_channel_representation, and the three detail payloads have DISJOINT key sets (eps_refined/max_rel_delta/refined_cost, e0_refined/rel_delta/refined_cost, reference/rel_error). A TypedDict here would either be a union a caller has to re-discriminate by probe anyway, or – with every key optional – a shape that promises nothing. detail is a diagnostic bag read by the supervising tuner loop and by reports, never a contract a solver depends on; the CONTRACT is converged and cost, which are precisely typed.

Parameters:

• value (Any)

• converged (bool)

• cost (int)

• detail (dict[str, Any])

value: Any

Alias for field number 0

converged: bool

Alias for field number 1

cost: int

Alias for field number 2

detail: dict[str, Any]

Alias for field number 3

class qscat.tuning.Refine2dReport

refine_to_2d_convergence’s third return value.

• converged: True when the larger of the two candidate relative moves fell below rtol; False when max_iter adopted steps were taken without ever meeting it – a real “did not settle” signal, not a silently accepted best effort.

• iterations: the adopted steps, in order. Empty when the input pair was already converged.

• final_value: the observable on the returned grid pair – the last adopted value, or the initial observable(g_r, g_R) if no step was adopted.

class qscat.tuning.Refine2dStep

One ADOPTED refinement step recorded by refine_to_2d_convergence.

• coordinate: which of the two grids was refined at this step.

• value: the observable on the newly adopted grid pair.

• rel_move: |value - previous| / max(|previous|, tiny), the relative move that made this coordinate the one adopted (it was the larger of the two candidate moves).

class qscat.tuning.TensorCost

tensor_cost’s report: the tensor-product (2-D) problem’s size, plus ROUGH anchored estimates of its sparse-LU cost. Only n_unknowns is exact; the three est_* fields are for RELATIVE ranking of candidate grids, never absolute prediction – see tensor_cost.

• n_unknowns: g_r.n * g_R.n, the tensor-grid unknown count. EXACT.

• est_nnz: estimated nonzeros in the assembled 2-D matrix, at a constant ~22.9 nonzeros per row.

• est_factor_gib: estimated peak resident memory of a SuperLU factorization, in GiB.

• est_factor_seconds: estimated SuperLU factorization wall-clock time, in seconds.

qscat.tuning.analyze_potential(V, x_min, x_max, m, e_max, *, n_sample=4000)

Analyze V on [x_min, x_max], returning its PotentialProfile.

V may return complex values (potentials do, on ECS tails); the analysis works on Re(V) and is only ever sampled at real x.

Parameters:

• V (Callable[[NDArray[float64]], NDArray[complexfloating | floating]])

• x_min (float)

• x_max (float)

• m (float)

• e_max (float)

• n_sample (int)

Return type:

PotentialProfile

qscat.tuning.equidistribution_elements(profile, order, *, phase_per_element, min_len, max_len)

Return real-region element lengths equidistributing de Broglie phase.

Boundaries are placed at Phi(x) = j * phase_per_element, where Phi(x) = cumulative_trapezoid(profile.k, profile.x). Where Phi fails to advance (classically forbidden: k ~ 0, kappa > 0), the local element length is instead capped by a kappa-based decay length. All lengths are then brought into [min_len, max_len] – oversized elements subdivided, undersized ones merged forward – and elements adjacent to a turning point or singularity are halved as a post-pass refinement. max_len is a HARD cap (never exceeded); min_len is a SOFT floor (met except where doing so would breach max_len). Every step only regroups/splits existing length, so the returned lengths always sum to the real-region domain span profile.x[-1] - profile.x[0].

Parameters:

• profile (PotentialProfile)

• order (int)

• phase_per_element (float)

• min_len (float)

• max_len (float)

Return type:

list[float]

qscat.tuning.grid_cost(grid)

A single grid’s discretisation cost: just its DVR point count.

Parameters:

grid (FemDvrEcsGrid)

Return type:

GridCost

qscat.tuning.interaction_region(model, *, r_probe=None, R_max=8.0, frac=0.02, n=400)

The (R_lo, R_hi) bracket over which V_int(r, R) is still transitioning between its low-R and high-R regimes.

Defines s(R) = max_r |Re(model.v_int(r_probe, R))| over a fixed electronic probe set r_probe, scans R on linspace(1e-3, R_max, n), and returns the first/last R whose s(R) falls within the middle (1 - 2*frac) fraction of s’s own [min, max] range on that scan – i.e. R_lo is the first R with s(R) - s.min() >= frac * (s.max() - s.min()) and R_hi is the last R with s(R) - s.min() <= (1 - frac) * (s.max() - s.min()). See the module docstring for why this is the range that localises (a raw threshold against s.max() alone does not, for these sigmoid-shaped interaction profiles).

R starts at 1e-3 rather than 0 to stay clear of any R=0 singularity in a model’s v_int (e.g. a Coulomb-tail ionic model).

Parameters:

• model (ResonanceModel)

• r_probe (FloatArray | None)

• R_max (float)

• frac (float)

• n (int)

Return type:

tuple[float, float]

qscat.tuning.max_stable_angle(V, R0, tail_extent, *, angle_cap=35.0, n_probe=40)

Largest ECS rotation angle (deg) for which V stays bounded on the tail.

Scans n_probe angles from angle_cap / n_probe up to angle_cap (ascending); at each, evaluates V on the rotated contour over _N_TAIL_SAMPLES points spanning [R0, R0 + tail_extent]. The first angle at which |V| grows anywhere along the tail (relative to its running minimum, beyond _GROWTH_RELTOL; a minimum below _GROWTH_ABS_FLOOR counts as that floor, so a potential already at round-off cannot “grow”) stops the scan; the last angle that did NOT diverge is returned. Never exceeds angle_cap.

Parameters:

• V (Callable[[NDArray[complexfloating]], ArrayLike])

• R0 (float)

• tail_extent (float)

• angle_cap (float)

• n_probe (int)

Return type:

float

qscat.tuning.optimal_real_mesh(profile, *, orders=(6, 8, 10, 14), phase_coeff=0.1, min_len, max_len)

Sweep DVR orders, returning the (mesh, order) with fewest points.

For each candidate order, phase_per_element = phase_coeff * (order - 1) sets the per-element phase budget; the resulting mesh’s DVR point count is estimated as len(elements) * (order - 1). The h/p optimum is the combination minimizing that count.

Parameters:

• profile (PotentialProfile)

• orders (tuple[int, ...])

• phase_coeff (float)

• min_len (float)

• max_len (float)

Return type:

tuple[list[float], int]

qscat.tuning.order_for_wavenumber(k, element_len, *, target_ppw=6.0, orders=(6, 8, 10, 14))

Smallest DVR order in orders resolving wavenumber k at a fixed real element length element_len, in points-per-wavelength.

Points-per-wavelength at a given order is order * (2*pi/k) / element_len – order DVR points span each element, and 2*pi/k is the de Broglie wavelength of the wave being resolved. Returns the smallest order (searched in the given, ascending order) with ppw >= target_ppw; falls back to max(orders) if none qualifies (a caller asking to resolve an implausibly fast wave at a fixed element length still gets the best available order rather than an exception).

k <= 0 (no wave to resolve) returns orders[0] – the cheapest order, since there is nothing to resolve.

Parameters:

• k (float)

• element_len (float)

• target_ppw (float)

• orders (tuple[int, ...])

Return type:

int

qscat.tuning.probe_channel_representation(grid, k, l, *, charge=0, mass=1.0, rtol=0.001)

Is the energy-normalized channel function F_{k,l} REPRESENTABLE on grid’s real region?

Compares the grid’s own composite Gauss-Lobatto quadrature estimate of integral |F|^2 dr over the real region [x_min, R0) – sum(Re(w_j) * F(r_j)^2) at the grid’s own real-region nodes/weights – against a fine-uniform-grid Simpson reference (_fine_reference_norm, integrated over that SAME [x_min, R0] span – grid.spec.x_min, not hardcoded to 0.0; every grid built by qscat.core.grids uses x_min=0.0 today, but GridSpec.x_min is a supported nonzero field). The node exactly at R0 (shared with the first ECS-tail element) is EXCLUDED: its global bridge-summed weight mixes in that neighbor’s complex Jacobian (see the inline comment at the mask below), which would corrupt a real-region-only quadrature at any nonzero ECS angle. charge=0 uses riccati_bessel_en_ mass(…, mass) (mass-mass free wave; reduces to the mass-1 electron channel at the default mass=1.0) and charge != 0 uses coulomb_f_en (which already takes mass) – both the DVR sum and the reference honor mass identically. No eigensolve is needed, so this is by far the cheapest of the three probes – and the most diagnostic: a grid whose element lengths are not small compared to the channel’s wavelength 2*pi/k cannot represent F and this probe catches it directly (this is what would have caught the K~58-on-1.0-bohr-elements coarse-grid DA failure).

Parameters:

• grid (FemDvrEcsGrid)

• k (float)

• l (int)

• charge (int)

• mass (float)

• rtol (float)

Return type:

ProbeResult

qscat.tuning.probe_electronic(model, elec_grid, R, *, window, rtol=0.001)

Convergence of the lowest anion bound electronic-state energy at internuclear distance R, under one h-refinement of elec_grid.

See the module docstring for why this bound-state energy (rather than a full two-angle resonance-pole match) is used as the electronic-grid convergence proxy; window is accepted for signature parity with a future pole-based probe but is UNUSED here.

Parameters:

• model (ResonanceModel)

• elec_grid (FemDvrEcsGrid)

• R (float)

• window (tuple[float, float] | None)

• rtol (float)

Return type:

ProbeResult

qscat.tuning.probe_nuclear(model, nuclear_grid, n_vib, *, rtol=0.001)

Convergence of the n_vib lowest vibrational eigenvalues under one h-refinement of nuclear_grid.

converged = max_v |eps_v(refined) - eps_v| / |eps_v| < rtol (an absolute floor is used for any eps_v too close to zero). cost is the (unrefined) grid’s point count nuclear_grid.n.

Parameters:

• model (ResonanceModel)

• nuclear_grid (FemDvrEcsGrid)

• n_vib (int)

• rtol (float)

Return type:

ProbeResult

qscat.tuning.propose_grid(model, coordinate, energy_range, *, incident=None, phase_coeff=None, channel='ve', elec_grids=None, resonance_n_dense=25)

The one-shot a-priori FemDvrEcsGrid for model/coordinate over energy_range = (E_min, E_max).

Pipeline: a small per-coordinate MODEL ADAPTER (_nuclear_adapter / _electronic_adapter) picks V/mass/extent/channel-k -> analyze_ potential profiles it at E_max -> optimal_real_mesh sweeps the h/p equidistribution mesh -> max_stable_angle finds the ECS rotation angle -> tune_ecs_tail sizes the absorbing tail for the channel wavenumber -> the resulting real + tail ElementSpec list becomes a GridSpec / FemDvrEcsGrid.

phase_coeff, if given, overrides optimal_real_mesh’s calibrated default de-Broglie phase-per-(order-1) constant C – the knob validation.tuning.calibrate sweeps to find that calibrated value against the eMoScat decks. None (the default) leaves optimal_real_mesh at its own calibrated default; ordinary callers never need to pass this.

incident (qscat.tuning.incident.IncidentSpec) is accepted here as BOTH an extent floor AND a resolution floor:

• EXTENT: incident.required_extent() extends the real-region cutoff to at least that value before the mesh/ECS steps run.

• RESOLUTION: incident.incident_energy() raises the effective E_max fed to analyze_potential/optimal_real_mesh to max(E_max, incident_energy). Without this, a HAND-BUILT IncidentSpec whose impulse implies an energy above energy_range’s own E_max would still get a mesh sized only for the (lower) energy_range’s local wavenumber – silently under-resolving the incident’s actual wave. (A tw_analysis-produced IncidentSpec never triggers this: its energy is bounded by energy_range by construction, so max(…) is then a no-op.) The ECS-tail channel_k is deliberately left keyed to energy_range’s own E_max alone, not this raised value: the incident wavepacket is a TD object CONTAINED in the real region by the extent floor above, not an outgoing wave the tail need absorb – widening channel_k too is a plausible future refinement, not required by this baseline.

The parameter is a real IncidentSpec (no duck typing – a wrong object now fails loudly at the call instead of silently contributing 0.0); the placement logic itself (impulse/width/observation boundary, tw_analysis) lives in qscat.tuning.incident – this is only the extent/resolution floor.

channel selects which physical channel the mesh targets:

• “ve” (the default): the VE (vibrational-excitation) path, v0-alone, exactly as before this parameter existed – BYTE-IDENTICAL to the pre-channel behavior; nothing below reads elec_grids or resonance_n_dense on this path.

• “dissociation”, coordinate=”nuclear” only: the resonance-aware nuclear path (see _resonant_nuclear_mesh). Builds the adiabatic resonance curve (R, V_d(R), Gamma(R)) (qscat.tuning.resonance.resonance_curve_arrays, via a two-angle ECS pole match – elec_grids, if given, overrides the default electronic grids used for that match, and resonance_n_dense overrides its dense-sampling point count; both exist so tests can inject small/cheap grids). A REDUCED real-region extent (_RESONANT_NUCLEAR_X_MAX_DEFAULT, not the VE path’s _NUCLEAR_X_MAX_DEFAULT) is used, since the ECS tail absorbs the outgoing wave. The DVR order is sized (order_for_wavenumber) to resolve the fast dissociation EXIT wave K_exit at the base min_len, and that same K_exit (not spec.channel_k) drives the ECS tail. The narrow resonance CROSSING R* (the outermost Re(V_d) - v0 sign change) is then LOCALLY super-refined (refine_elements_in_window, overriding min_len only inside a Gamma-closing-width window around R*). The local override is the point: merging the resonance’s wavenumber into the global profile instead would be inert, because the merged profile’s finer elements get floored right back up by the shared global min_len. “dissociation” with coordinate=”electronic” raises ValueError – the resonance-aware path is nuclear-only.

• any other value raises ValueError.

Parameters:

• model (ResonanceModel)

• coordinate (Coordinate)

• energy_range (tuple[float, float])

• incident (IncidentSpec | None)

• phase_coeff (float | None)

• channel (Channel)

• elec_grids (tuple[FemDvrEcsGrid, FemDvrEcsGrid] | None)

• resonance_n_dense (int)

Return type:

FemDvrEcsGrid

qscat.tuning.refine(grid)

One h-refinement step: split every REAL element into two half-length elements; leave the ECS tail, quadrature order, and x_min untouched.

grid.spec is the GridSpec the grid was built from (FemDvrEcsGrid stores it directly), so no separate spec-plumbing is needed – this is the “SIMPLEST” option the design brief calls out: reconstruct from the grid’s own stored spec rather than threading a GridSpec through every probe signature.

Parameters:

grid (FemDvrEcsGrid)

Return type:

FemDvrEcsGrid

qscat.tuning.refine_elements_in_window(real_lengths, x_min, R_lo, R_hi, target_len)

Span-preservingly subdivide every real element overlapping [R_lo, R_hi] until each piece is <= target_len; elements entirely outside the window are returned untouched.

This is a LOCAL override of min_len – unlike optimal_real_mesh’s global equidistribution sweep, this targets one narrow feature (e.g. the resonance-crossing region a resonant nuclear mesh must super-refine) without perturbing element lengths anywhere else. An element only counts as “overlapping” if it shares more than a single boundary point with [R_lo, R_hi] (hi > R_lo and lo < R_hi), so a neighbor that merely touches the window at one endpoint is left alone. Every input length is either kept whole or subdivided into equal pieces, so sum(result) == sum(real_lengths) exactly (up to float round-off) – nothing outside the window is coarsened, and nothing inside it is left above target_len.

Parameters:

• real_lengths (list[float])

• x_min (float)

• R_lo (float)

• R_hi (float)

• target_len (float)

Return type:

list[float]

qscat.tuning.refine_to_2d_convergence(observable, g_r, g_R, *, rtol=0.01, max_iter=4)

Iteratively refine whichever of g_r (electronic) / g_R (nuclear) grid is under-resolved for observable, until it stops moving.

observable(g_r, g_R) -> float is an arbitrary caller-supplied closure (e.g. a real cross-section at a fixed energy, or a synthetic test function) – this loop does not inspect what it computes.

Each of up to max_iter iterations:

1. Evaluates observable on a once-refine`d nuclear variant (`g_r, refine(g_R)) and a once-refine`d electronic variant (`refine(g_r), g_R), each against the CURRENT (adopted) grids.

2. Compares each candidate’s relative move |v - current| / max(|current|, tiny) off the current adopted value.

3. If the larger of the two relative moves is < rtol: STOP, converged=True, and do NOT record a step for this check.

4. Otherwise ADOPT the coordinate with the larger relative move (replace that grid with its refinement, update the current value), and record {“coordinate”: …, “value”: …, “rel_move”: …}.

If the loop exhausts max_iter adopted steps without ever satisfying step 3, converged=False – a genuine “did not settle” signal, not a silently accepted best-effort result.

Returns (g_r, g_R, detail): the (possibly refined) grid pair, and a Refine2dReport – a plain dict with exactly “converged”, “iterations” (a list of Refine2dStep, in order) and “final_value”; see those two definitions for each key’s meaning and type.

Parameters:

• observable (Callable[[FemDvrEcsGrid, FemDvrEcsGrid], float])

• g_r (FemDvrEcsGrid)

• g_R (FemDvrEcsGrid)

• rtol (float)

• max_iter (int)

Return type:

tuple[FemDvrEcsGrid, FemDvrEcsGrid, Refine2dReport]

qscat.tuning.required_extent(spec)

Free-function form of IncidentSpec.required_extent, for callers who prefer a function over a bound method.

Parameters:

spec (IncidentSpec)

Return type:

float

qscat.tuning.resonance_curve_arrays(model, elec_grid_a, elec_grid_b, *, R_max=22.0, n_dense=25, region=None)

Efficient adiabatic resonance curve (R, V_d(R), Gamma(R)).

Samples the pole continuation (qscat.core.lcp.resonance_pole_walk) DENSELY (n_dense points) only inside the interaction window region = region or interaction_region(model), plus a few inner points below R_lo and a SINGLE far point at R_max standing in for the saturated asymptote – the efficiency constraint that keeps this to O(n_dense) electronic solves rather than a scan over hundreds of R points. The walk runs over the samples DESCENDING (outer -> inner) and its freeze (see resonance_pole_walk) correctly holds V_d constant across the big R_max -> R_hi gap, since Gamma has already saturated to ~0 out there.

Renamed from resonance_curve (2026-08-25 API surface pass) to end the collision with qscat.core.bo.resonance_curve, which shares the same underlying pole walk but returns an ElectronicCurves carrying the eigenVECTORS, for building Born-Oppenheimer basis states. This one returns plain (R, V_d, Gamma) arrays – the name says so – and exists only to size a grid, so it discards the states and samples as sparsely as it can.

Seeded from the bound anion state at R_max (qscat.core.dissociation.anion_electronic_states). Returns arrays ascending in R; V_d(R) = Re(model.v0(R)) + shift.

Parameters:

• model (ResonanceModel)

• elec_grid_a (FemDvrEcsGrid)

• elec_grid_b (FemDvrEcsGrid)

• R_max (float)

• n_dense (int)

• region (tuple[float, float] | None)

Return type:

tuple[FloatArray, FloatArray, FloatArray]

qscat.tuning.tensor_cost(g_r, g_R)

ROUGH, ANCHORED cost ESTIMATES for the tensor-product 2-D problem.

n_unknowns = g_r.n * g_R.n is exact (the tensor-grid unknown count). est_nnz, est_factor_gib, and est_factor_seconds are simple monotone scalings of n_unknowns, calibrated so the CLAUDE.md/docs/physics/ mumps-sparse-backend.md ~143k-unknown N2 production-deck SuperLU numbers (260 s factor time, 7.4 GiB peak RSS, ~22.9 nonzeros/row) are ROUGHLY reproduced at that anchor size. They are NOT a validated performance model: no assembly-time-dependence on quadrature order, ECS-tail bandwidth, or ordering algorithm is captured, and a MUMPS- backed factorization of the same matrix would be far cheaper (see grid_cost’s module docstring). Use these ONLY to rank candidate grids against each other (bigger n_unknowns always costs more here by construction), never to predict a real wall-clock time or memory footprint.

Parameters:

• g_r (FemDvrEcsGrid)

• g_R (FemDvrEcsGrid)

Return type:

TensorCost

qscat.tuning.tune_ecs_tail(K, R0, *, angle, order, tail_alpha=0.2, tail_skip=2, decay_target=1e-12, base=None)

Exp-growth ECS-tail element lengths absorbing wavenumber K at angle.

Returns element LENGTHS (not a grid) spanning [R0, R0 + L] where L = -ln(decay_target) / (K * sin(angle)) is the distance over which the fastest outgoing wave exp(-K * (x - R0) * sin(angle)) decays to decay_target. The first tail_skip elements are flat at base; the rest grow as base * exp(tail_alpha * (i - tail_skip + 1)) (the same growth law as qscat.core.grids._ecs_tail), appended until the running sum reaches L.

order is the DVR quadrature order the caller will attach to each tail element when building the actual grid; it does not affect the LENGTHS returned here, but is validated (order >= 2, matching GridSpec) since tune_ecs_tail is meant to be a pure function of (K, R0, angle, order, growth) per the ECS-tail tuner design.

base, if not given, defaults to L / 8 – a small fraction of the target tail span, chosen so tail_skip flat elements plus a handful of growing ones reach L without either an excessively coarse first element or an excessively long list.

Parameters:

• K (float)

• R0 (float)

• angle (float)

• order (int)

• tail_alpha (float)

• tail_skip (int)

• decay_target (float)

• base (float | None)

Return type:

list[float]

qscat.tuning.tw_analysis(model, energy_range)

Auto-place the TD Gaussian wavepacket for energy_range = (E_min, E_max) – the “TW analysis”. BEST-EFFORT and not independently calibrated against the eMoScat decks: the calibration (validation.tuning.calibrate) targets the mesh’s de-Broglie phase constant C (qscat.tuning.mesh._PHASE_COEFF), not this placement heuristic. Calibrating the placement itself is a documented follow-on (docs/physics/discretisation-tuning.md).

Physics: g(r) = exp(-(r-position)^2 / (2 sigma^2)) exp(i impulse r) has mean energy impulse**2 / 2 and, near the mean, an energy spread delta_E ~ impulse * delta_p with delta_p ~ 1 / (2 sigma) – i.e. delta_E ~ impulse / (2 sigma). Inverting [E_min, E_max] -> (impulse, sigma, position):

• impulse: p = sqrt(2 * E_centre), E_centre = (E_min+E_max)/2, so the packet’s mean energy sits at the range’s centre; launched INWARD (impulse = -p), matching n2_2d_td_cross_section’s wp_in convention.

• sigma: the spread must COVER the range, delta_E >= (E_max - E_min) / 2 <=> sigma <= p / (E_max - E_min); sigma is set to _SIGMA_SAFETY_FACTOR times that marginal bound so the realized spread comfortably brackets [E_min, E_max] rather than just touching its edges (a smaller sigma only widens the momentum/ energy spread further, so this is the conservative direction).

• position: _interaction_extent(model) + _STANDOFF_SIGMAS * sigma – far enough out that the packet starts genuinely outside the interaction region, with room to spare.

• observation: set equal to position (a sane symmetric default – the outgoing test function is read at the same radius the packet launches from; construct IncidentSpec directly for a different boundary).

Parameters:

• model (ResonanceModel)

• energy_range (tuple[float, float])

Return type:

IncidentSpec