# SPDX-FileCopyrightText: 2026 John C. Travers and contributors
#
# SPDX-License-Identifier: GPL-3.0-or-later
"""Translate constant Cartesian propagation into the schema-5 separable seam."""
import warnings
import numpy as np
from jaxtyping import Float
from solprop.boundaries import TEMPORAL_ALPHA_LENGTH, Absorber
from solprop.exceptions import GrazingIncidenceWarning, TranslationError
from solprop.grid import FreeGrid, TimeGrid
from solprop.spatial import SpatialDispersion
from solprop.translate.linop import DEFAULT_MARGIN, LinearOperator, tabulate
from solprop.translate.modal import DEFAULT_SPATIAL_NODE_BUDGET, Transform
from solprop.translate.spatial import SeparableOperator, SpatialGeometry
__all__ = ["free_operators", "free_transform"]
def _validate_span(n_pol: int, z_max: float, margin: float) -> float:
"""Validate spatial polarisation and finite positive translation extent.
Parameters
----------
n_pol : int
One or two independent Cartesian components.
z_max : float
Run length in m.
margin : float
Fractional nonnegative extension.
Returns
-------
float
Extended positive span in m.
"""
if isinstance(n_pol, bool) or n_pol not in (1, 2):
msg = "Cartesian polarisation count must be one or two"
raise TranslationError(msg)
span = z_max * (1 + margin)
if not np.isfinite(span) or z_max <= 0 or margin < 0:
msg = "Cartesian free operator requires a positive finite span"
raise TranslationError(msg)
return span
def _warn_grazing(space: FreeGrid, grid: TimeGrid, data: SeparableOperator) -> None:
"""Report the first active cutoff inside the split-square uncertainty bound.
Parameters
----------
space : solprop.grid.FreeGrid
Cartesian shape for the reported column coordinates.
grid : solprop.grid.TimeGrid
Frequency bins in rad/s and their active mask.
data : solprop.translate.spatial.SeparableOperator
Validated square high/low parts in m^-2.
"""
order = np.argsort(data.kperp_sq, kind="stable")
transverse = data.kperp_sq[order]
epsilon = np.finfo(float).eps
# Search every high-part candidate, including ties with different residuals.
# The first-order interval encloses all second-order-uncertain pairs without
# constructing a spatial-frequency volume; see the predeclared contract.
for omega_bin in np.flatnonzero(grid.active_mask):
material = data.k_squared[omega_bin]
radius = 64 * epsilon * abs(material) + 64 * np.nextafter(0.0, 1.0)
first = np.searchsorted(transverse, material.real - radius, side="left")
last = np.searchsorted(transverse, material.real + radius, side="right")
columns = order[first:last]
candidates = data.kperp_sq[columns]
difference = (material.real - candidates) + (
data.k_squared_low[omega_bin].real - data.kperp_sq_low[columns]
)
uncertainty = (
8 * epsilon * np.abs(difference)
+ 64 * epsilon**2 * (abs(material) + candidates)
+ 32 * np.nextafter(0.0, 1.0)
)
unresolved = np.flatnonzero(np.abs(difference) <= uncertainty)
if unresolved.size == 0:
continue
selected = int(unresolved[0])
column = int(columns[selected])
msg = (
"Cartesian grazing-incidence rounding: active omega bin "
f"{omega_bin} (omega={grid.omega[omega_bin]:.17g} rad/s), transverse "
f"column {column} (ix={column % space.n_x}, iy={column // space.n_x}) "
f"has Re(k²)-kperp²={difference[selected]:.17g} m^-2 "
f"within arithmetic uncertainty {uncertainty[selected]:.17g} "
"m^-2 after split-square correction. Roundoff can still change the "
"clamp branch and nonlinear 1/Re(k_z) prefactor. Continuing with "
"the corrected coefficients and existing nonpositive clamp; "
"this grid diagnostic does not measure the affected field strength."
)
warnings.warn(msg, GrazingIncidenceWarning, stacklevel=3)
return
[docs]
def free_operators( # noqa: PLR0913 - geometry, operator settings and controls
space: FreeGrid,
dispersion: SpatialDispersion,
*,
n_pol: int,
z_max: float,
absorber: Absorber | None = None,
margin: float = DEFAULT_MARGIN,
transverse_rate: Float[np.ndarray, " n_modes"] | None = None,
) -> tuple[LinearOperator, SeparableOperator]:
"""Keep constant Cartesian material coefficients in separable form.
Parameters
----------
space : solprop.grid.FreeGrid
Cartesian transverse wave-number squares in m^-2.
dispersion : solprop.spatial.SpatialDispersion
Constant bulk material, temporal grid and frame.
n_pol : int
State polarisation count, 1 or 2.
z_max : float
Positive propagation length in m.
absorber : solprop.boundaries.Absorber or None, optional
Shared spectral absorption rate in m^-1.
margin : float, optional
Relative extension of the valid interval, nonnegative.
transverse_rate : numpy.ndarray or None, optional
A transverse collar's power rate per column in m^-1
(:func:`solprop.boundaries.transverse_collar_rate`); ``None`` is zero.
Returns
-------
tuple
Linear operator and the same separable record for the nonlinear prefactor.
Raises
------
TranslationError
If the medium varies along z or the interval is invalid.
Warns
-----
solprop.exceptions.GrazingIncidenceWarning
If an active transverse/frequency bin lies within the elementary-arithmetic
uncertainty of the grazing clamp branch. Translation continues unchanged;
this diagnostic does not assess the field or nonlinear source in that bin.
"""
if not all(path.is_constant for path in dispersion.medium.parameter_paths.values()):
msg = (
"Cartesian free propagation requires constant media; "
"varying density profiles are unsupported"
)
raise TranslationError(msg)
span = _validate_span(n_pol, z_max, margin)
grid = dispersion.grid
_, frame, numerator, clamped = dispersion.shared_values()
square, square_low, transverse_low = dispersion.cartesian_squares(space)
data = SeparableOperator(
k_squared=square,
kperp_sq=space.kperp_sq,
k_squared_low=square_low,
kperp_sq_low=transverse_low,
frame=frame,
spectral_rate=np.zeros(grid.n_omega)
if absorber is None
else absorber.spectral_rate,
prefactor_numerator=numerator * grid.crop_scale,
prefactor_clamped=clamped * grid.crop_scale,
active=grid.active_mask,
evanescent_rate=dispersion.evanescent_rate,
n_pol=n_pol,
z_span=(0.0, span),
kperp_rate=transverse_rate,
)
_warn_grazing(space, grid, data)
beta0, beta1 = dispersion.frame_quantities(0.0)
return LinearOperator(
phi_representation="separable",
phi_breakpoints=np.array([0.0, span]),
phi_coefficients=None,
phi_degrees=np.ones(1, dtype=np.intp),
hard_kinks=np.empty(0),
phi_constant=None,
phi_expansion=None,
beta0=tabulate(lambda z: np.full_like(z, beta0), 0.0, span),
beta1=tabulate(lambda z: np.full_like(z, beta1), 0.0, span),
separable=data,
), data