solprop.stats#

The catalogue of plane statistics, with one definition of each.

Statistics belong to saved planes, not to steps (see docs/dev/decisions.md). Evaluating every statistic on every accepted step, with a full inverse transform each time and a root-finding zero-dispersion search among them, would cost about a fifth of a step and dominate device runs; here they are evaluated where a field is saved, and an optional periodic subset (StatsSpec.periodic) covers step-level diagnostics.

Every definition here is written once, on internal FFT-convention arrays of shape (n_modes, n_pol, n_omega), and reduces in a fixed order. There is one definition per statistic in the run, not a host one and a device one that compute the same quantity differently, and a backend’s fused pass is validated against it. The analytic signal comes from solprop.fourier.analytic_signal() alone.

One quantity is measured differently after the run, and deliberately: a run forms its temporal power as the analytic Poynting flux from the internal amplitude and the plane’s index, while solprop.processing.time, given a saved plane and no mode, measures the squared true-power envelope |b(t)|^2. The two differ at second order in the spread of Re n_eff across the spectrum and have identical time integrals; see docs/guide/results/processing.md and the docstrings of solprop.processing.time.peak_power() and solprop.processing.time.fwhm_t().

Spatial statistics use field units and the geometry’s transverse integration measure. Their host definitions are in solprop.spatial; periodic backend reductions are qualified against those definitions. Plasma statistics need the compiled response, and spatial widths need the geometry declared at assembly.

Undefined is not an error. A mode that carries no light has no spectral centroid, and a trace that never falls to half its maximum – a continuous wave, an empty mode – has no full width. Those return nan, which the docstring of each states and the manual page repeats. A missing dependency, by contrast – a zero-dispersion wavelength with no mode, a windowed energy with no windows – raises solprop.exceptions.OutputError.

The manual page is docs/guide/results/statistics.md.

Module Attributes

`PROFILE_STATISTICS`	Statistics that read a `solprop.profiles.ZProfile` at the plane's position.
`ON_AXIS_PLASMA_STATISTICS`	The plasma statistics evaluated on the waveguide's axis (SEP 19, section 5.2).
`PLASMA_INTEGRATION_STATISTICS`	The field plasma column's run-time integration diagnostics (SEP 18, section 4.5).
`DEFAULT_ZDW_RANGE`	Search window of `zdw()` when no narrower one is given, m.
`CATALOGUE`	Every statistic solprop knows about, implemented or reserved, keyed by name.
`GENERAL_MODAL_ACCESSORS`	Every accessor a three-component general modal basis admits, by name.
`WITHHELD_WITH_A_FLOW`	The general-basis accessors a route withholds when its block couples, by name.
`ACCESSOR_ADVICE`	What to use instead of each legacy name on a general modal basis.
`METRIC_WEIGHTED_STATISTICS`	Statistics whose value depends on the basis's flux metric, by name.
`ApertureWeighting`	How an aperture's profile weights the collected spectrum (SEP 23 section 4.5).
`MAX_APERTURES`	The most apertures one `ApertureSet` carries.
`MAX_REGIONS`	The most rows `StatsSpec.spectral_regions` carries.
`WEIGHT_DEPENDENT_STATISTICS`	Names whose value depends on the modal power weight, so a backend evaluates it.

Functions

`analytic_field`(grid, field)	Return the time-domain analytic signal of a stored spectrum.
`analytic_flux`(grid, field, *[, metric])	Return the peak of the plane's same-frequency analytic flux, in watts.
`check_aperture_choices`(profile, weighting)	Refuse, by name, an aperture profile or weighting that does not exist.
`collect`(spec, grid, plane, *[, modes, ...])	Evaluate every statistic a spec asks for at one saved plane.
`default_spec`(kind, *[, with_profiles, with_zdw])	Return the statistics a run of the given kind records unless told otherwise.
`energy`(grid, field)	Return the pulse energy of each mode and polarisation, in joules.
`energy_window`(grid, field, windows)	Return the energy inside each spectral window, in joules.
`fwhm_t`(grid, field, *[, metric])	Return the width of the central feature at half maximum, in seconds.
`fwhm_t_outer`(grid, field, *[, metric])	Return the total extent above half maximum, in seconds.
`inband_wavelengths`(grid)	Return the shortest and longest in-band wavelength of a grid, in metres.
`omega0`(grid, field, *[, metric])	Return the centroid of the spectral energy density, per mode and polarisation.
`peak_power`(grid, field)	Return the largest instantaneous power each mode reaches, in watts.
`peak_power_total`(grid, field)	Return the largest instantaneous power of the whole field, in watts.
`peak_power_window`(grid, field, windows)	Return the peak power of each mode inside each spectral window, in watts.
`physical_energy`(grid, field, *[, metric])	Return the physical energy through the plane, in joules.
`profile_value`(profile, z)	Return the value of a z-profile at one position.
`species_statistics`(densities, fractions, ...)	Reduce species populations at the node of largest total electron density.
`spectral_flux`(grid, field, *[, metric])	Return the physical energy's spectral density, in J/(rad/s).
`weight_refusal_for`(names, weight)	Return why a general statistic cannot be formed under this weight, or `None`.
`zdw`(modes, z, *[, wavelength_range])	Return the zero-dispersion wavelength of each mode at one position.

Classes

`ApertureSet`(centres, diameters, distances[, ...])	Circular apertures in a far-field or Fourier plane behind a free-space run.
`StatDefinition`(name, summary, shape, units)	What one statistic is called, what it returns and what it needs.
`StatsSpec`([names, energy_windows, periodic, ...])	Which statistics a run asks for, and the windows and period they need.

solprop.stats.ACCESSOR_ADVICE: Mapping[str, str] = {'energy': 'physical_energy, which integrates the same-frequency metric product', 'energy_window': 'physical_energy with an explicit band, or the reconstructed E/H of a saved plane', 'peak_intensity': 'analytic_flux or reconstructed E/H with an explicit detector definition', 'peak_power': 'analytic_flux or reconstructed E/H with an explicit detector definition', 'peak_power_total': 'analytic_flux or reconstructed E/H with an explicit detector definition', 'peak_power_window': 'analytic_flux or reconstructed E/H with an explicit detector definition'}#

What to use instead of each legacy name on a general modal basis.

A name with no entry here gets the generic advice of solprop.output.capability: a general modal basis has physical_energy, analytic_flux and the reconstruction accessors, and nothing else of that family.

class solprop.stats.ApertureSet(centres, diameters, distances, edges=None, profile='hard', weighting='amplitude', labels=None)[source]#

Bases: object

Circular apertures in a far-field or Fourier plane behind a free-space run.

The tables of aperture_spectrum and aperture_on_axis_spectrum (SEP 23): each aperture collects the light whose paraxial image x = z c kappa / omega falls inside it, a chromatic weight on the stored transverse spectrum. Every position is physical: a beam travelling towards positive x lands at positive x. Equality is by identity; the arrays are owned read-only copies.

Attributes:

centresnumpy.ndarray: Aperture centres (x, y) in the mask plane, in m, shape (n_apertures, 2).
diametersnumpy.ndarray: Diameters in m, finite and positive, shape (n_apertures,). On a line the aperture is a slit of this full width along x.
distancesnumpy.ndarray: Distance from the exit face to the mask plane, or a lens’s focal length for a Fourier plane, in m, finite and positive, shape (n_apertures,).
edgesnumpy.ndarray or None: Edge widths w in m for tanh and rcos, orders n >= 2 for supergauss, shape (n_apertures,); None for hard, and only for hard.
profile{‘hard’, ‘tanh’, ‘rcos’, ‘supergauss’}: The edge profile the set shares (solprop.spatial.aperture_profile()). The default is hard.
weighting{‘amplitude’, ‘quadrature’}: amplitude (the default) models a mask of amplitude transmission W and weights the collected spectrum by W**2; quadrature uses a smooth profile once, at the corrected diameter, as a quadrature weight for a hard hole. Quadrature is refused with supergauss.
labelstuple of str: One unique non-empty label per aperture, recorded with the run’s configuration; None (the default) labels them "0", "1", …

Raises:

OutputError: If there are no apertures or more than MAX_APERTURES; an array is not finite or not of its stated shape; a diameter or distance is not positive; an edge is missing, given for hard, or out of its range (0 < w <= D/4 for rcos and for tanh under amplitude weighting, 0 < w <= D/8 for tanh under quadrature weighting, n >= 2 for supergauss); quadrature weighting is asked of supergauss; an aperture reaches the paraxial disc’s edge, |X| + D/2 (+ w for rcos) >= z; or the labels are not unique non-empty strings, one per aperture.

Warns:

solprop.exceptions.ApertureEdgeWarning: Once per set, for tanh under amplitude weighting with r/w < 4 (w > D/8): the profile slopes at its centre, so a bright beam there is collected with a lattice error near 1e-5. The set is built as asked; rcos is the recommended smooth profile.

Examples

>>> from solprop.stats import ApertureSet
>>> holes = ApertureSet(
...     centres=[[4e-3, 4e-3], [-4e-3, -4e-3]],
...     diameters=[1.5e-3, 1.5e-3],
...     distances=[0.25, 0.25],
...     edges=[40e-6, 40e-6],
...     profile="tanh",
...     weighting="quadrature",
...     labels=("probe", "signal"),
... )
>>> holes.n_apertures, holes.labels
(2, ('probe', 'signal'))

centres: Sequence[Sequence[float]] | ndarray#

diameters: Sequence[float] | ndarray#

distances: Sequence[float] | ndarray#

edges: Sequence[float] | ndarray | None = None#

labels: tuple[str, ...] | None = None#

property n_apertures: int[source]#

Return how many apertures the set holds.

Returns:

int: The number of apertures.

profile: Literal['hard', 'tanh', 'rcos', 'supergauss'] = 'hard'#

same_as(other)[source]#

Return whether another set describes the same apertures, entry for entry.

Parameters:

otherApertureSet or None: The set to compare.

Returns:

bool: True when every table, the profile, the weighting and the labels agree exactly: a structural comparison of two requests.

weighting: Literal['amplitude', 'quadrature'] = 'amplitude'#

solprop.stats.ApertureWeighting#

How an aperture’s profile weights the collected spectrum (SEP 23 section 4.5).

alias of Literal[‘amplitude’, ‘quadrature’]

solprop.stats.CATALOGUE: Mapping[str, StatDefinition] = {'absorbed_power': StatDefinition(name='absorbed_power', summary='Power each absorbing channel of the linear operator removes.', shape='(4,)', units='J/m', requires=('spatial_geometry', 'separable_operator'), implemented=True, backend=True, between_planes=True), 'analytic_flux': StatDefinition(name='analytic_flux', summary='Peak total analytic Poynting flux of a general modal basis.', shape='()', units='W', requires=(), implemented=True, backend=False, between_planes=True), 'aperture_on_axis_spectrum': StatDefinition(name='aperture_on_axis_spectrum', summary='On-axis spectral fluence behind each aperture in a Fourier plane.', shape='(n_apertures, n_omega)', units='J s/m^2', requires=('apertures', 'spatial_geometry'), implemented=True, backend=False, between_planes=False), 'aperture_spectrum': StatDefinition(name='aperture_spectrum', summary='Spectrum collected through each far-field aperture.', shape='(n_apertures, n_omega)', units='J/(rad/s)', requires=('apertures', 'spatial_geometry'), implemented=True, backend=False, between_planes=False), 'beam_radius': StatDefinition(name='beam_radius', summary='Fluence second-moment and 1/e^2 radii; NaN if undefined.', shape='(2,)', units='m', requires=('spatial_geometry',), implemented=True, backend=False, between_planes=False), 'beam_spectrum': StatDefinition(name='beam_spectrum', summary='Whole-beam spectral energy density.', shape='(n_omega,)', units='J/(rad/s)', requires=('spatial_geometry',), implemented=True, backend=False, between_planes=False), 'core_radius': StatDefinition(name='core_radius', summary='Core radius of the waveguide at the plane.', shape='()', units='m', requires=('core_radius',), implemented=True, backend=False, between_planes=True), 'density': StatDefinition(name='density', summary='Number density of the filling gas at the plane.', shape='()', units='1/m^3', requires=('density',), implemented=True, backend=False, between_planes=True), 'electron_density': StatDefinition(name='electron_density', summary='Free-electron density left behind by the pulse.', shape='()', units='1/m^3', requires=(), implemented=True, backend=True, between_planes=True), 'energy': StatDefinition(name='energy', summary='Modal column energies or whole spatial beam energy.', shape='(n_modes, n_pol); () for spatial beams', units='J', requires=(), implemented=True, backend=False, between_planes=True), 'energy_window': StatDefinition(name='energy_window', summary='Pulse energy inside each spectral window.', shape='(n_windows, n_modes, n_pol)', units='J', requires=('windows',), implemented=True, backend=False, between_planes=True), 'fwhm_r': StatDefinition(name='fwhm_r', summary='Central fluence FWHM; NaN if undefined.', shape='()', units='m', requires=('spatial_geometry',), implemented=True, backend=False, between_planes=False), 'fwhm_t': StatDefinition(name='fwhm_t', summary='Central half-maximum width of the total Poynting power.', shape='()', units='s', requires=(), implemented=True, backend=False, between_planes=True), 'fwhm_t_outer': StatDefinition(name='fwhm_t_outer', summary='Total extent of Poynting power above half its maximum.', shape='()', units='s', requires=(), implemented=True, backend=False, between_planes=True), 'omega0': StatDefinition(name='omega0', summary='Centroid of the spectral energy density.', shape='(n_modes, n_pol)', units='rad/s', requires=(), implemented=True, backend=False, between_planes=True), 'on_axis_electron_density': StatDefinition(name='on_axis_electron_density', summary='Free-electron density left behind on the axis.', shape='()', units='1/m^3', requires=(), implemented=True, backend=True, between_planes=False), 'on_axis_intensity': StatDefinition(name='on_axis_intensity', summary='Largest intensity in time on the spatial beam axis.', shape='()', units='W/m^2', requires=('spatial_geometry',), implemented=True, backend=False, between_planes=True), 'on_axis_peak_ionisation_rate': StatDefinition(name='on_axis_peak_ionisation_rate', summary='Peak rate on the axis; cycle-averaged for envelopes.', shape='()', units='1/s', requires=(), implemented=True, backend=True, between_planes=False), 'on_axis_species_electron_density': StatDefinition(name='on_axis_species_electron_density', summary="Each species' electron density left behind on the axis.", shape='(n_species,)', units='1/m^3', requires=(), implemented=True, backend=True, between_planes=False), 'on_axis_species_peak_ionisation_rate': StatDefinition(name='on_axis_species_peak_ionisation_rate', summary='Species peak rate on the axis; cycle-averaged for envelopes.', shape='(n_species,)', units='1/s', requires=(), implemented=True, backend=True, between_planes=False), 'on_axis_spectrum': StatDefinition(name='on_axis_spectrum', summary='On-axis spectral fluence density summed over polarisation.', shape='(n_omega,)', units='J s/m^2', requires=('spatial_geometry',), implemented=True, backend=False, between_planes=False), 'peak_intensity': StatDefinition(name='peak_intensity', summary='Spatial vacuum intensity; modal on-axis analytic Poynting flux.', shape='()', units='W/m^2', requires=(), implemented=True, backend=True, between_planes=True), 'peak_ionisation_rate': StatDefinition(name='peak_ionisation_rate', summary='Peak ionisation rate; cycle-averaged for envelopes.', shape='()', units='1/s', requires=(), implemented=True, backend=True, between_planes=True), 'peak_power': StatDefinition(name='peak_power', summary='Peak analytic Poynting power of each mode during a run.', shape='(n_modes,)', units='W', requires=(), implemented=True, backend=False, between_planes=True), 'peak_power_total': StatDefinition(name='peak_power_total', summary='Peak total analytic Poynting power during a run.', shape='()', units='W', requires=(), implemented=True, backend=False, between_planes=True), 'peak_power_window': StatDefinition(name='peak_power_window', summary='Peak analytic Poynting power per mode and spectral window.', shape='(n_windows, n_modes)', units='W', requires=('windows',), implemented=True, backend=False, between_planes=True), 'physical_energy': StatDefinition(name='physical_energy', summary='Physical energy of a general modal basis, metric-weighted.', shape='()', units='J', requires=(), implemented=True, backend=False, between_planes=True), 'plasma_boundary_level': StatDefinition(name='plasma_boundary_level', summary='Free-current integrand at the window start, relative to its peak.', shape='()', units='1', requires=(), implemented=True, backend=True, between_planes=False), 'plasma_undershoot': StatDefinition(name='plasma_undershoot', summary="Largest exposure undershoot of the plasma column's band integral.", shape='()', units='1', requires=(), implemented=True, backend=True, between_planes=False), 'plasma_work_defect': StatDefinition(name='plasma_work_defect', summary='Loss-current work against the ionisation it paid, relative.', shape='()', units='1', requires=(), implemented=True, backend=True, between_planes=False), 'pressure': StatDefinition(name='pressure', summary='Pressure of the filling gas at the plane.', shape='()', units='bar', requires=('pressure',), implemented=True, backend=False, between_planes=True), 'quadrature_error': StatDefinition(name='quadrature_error', summary='Relative error of the transverse quadrature on the polarisation.', shape='()', units='dimensionless', requires=(), implemented=True, backend=True, between_planes=True), 'region_spectrum': StatDefinition(name='region_spectrum', summary='Spectrum in each box of transverse wavevector.', shape='(n_regions, n_omega)', units='J/(rad/s)', requires=('regions', 'spatial_geometry'), implemented=True, backend=False, between_planes=False), 'species_density': StatDefinition(name='species_density', summary='Component number densities at this plane.', shape='(n_species,)', units='1/m^3', requires=(), implemented=True, backend=True, between_planes=False), 'species_electron_density': StatDefinition(name='species_electron_density', summary='Species electron densities at the densest total-electron node.', shape='(n_species,)', units='1/m^3', requires=(), implemented=True, backend=True, between_planes=True), 'species_ionisation': StatDefinition(name='species_ionisation', summary='Final ionised fraction of each species at every node.', shape='(n_species, n_nodes)', units='1', requires=(), implemented=True, backend=True, between_planes=True), 'species_ionisation_energy': StatDefinition(name='species_ionisation_energy', summary='Final ionised inventory energy per propagation length.', shape='(n_species,)', units='J/m', requires=(), implemented=True, backend=True, between_planes=False), 'species_peak_ionisation_rate': StatDefinition(name='species_peak_ionisation_rate', summary='Species peak rate over all nodes; cycle-averaged for envelopes.', shape='(n_species,)', units='1/s', requires=(), implemented=True, backend=True, between_planes=False), 'spectral_flux': StatDefinition(name='spectral_flux', summary="Spectral density of a general modal basis's physical energy.", shape='(n_omega,)', units='J/(rad/s)', requires=(), implemented=True, backend=False, between_planes=False), 'zdw': StatDefinition(name='zdw', summary='Zero-dispersion wavelength of each mode at the plane.', shape='(n_modes,)', units='m', requires=('modes',), implemented=True, backend=False, between_planes=False)}#: Every statistic solprop knows about, implemented or reserved, keyed by name.

solprop.stats.DEFAULT_ZDW_RANGE = (1e-07, 3e-06)#: Search window of zdw() when no narrower one is given, m.

solprop.stats.GENERAL_MODAL_ACCESSORS: frozenset[str] = frozenset({'analytic_flux', 'channel_labels', 'coordinate_intensity', 'electric_field', 'fwhm_t', 'fwhm_t_outer', 'instantaneous_flux', 'magnetic_field', 'omega0', 'physical_energy', 'plane_nodes', 'reconstructed_flux', 'spectral_flux', 'zdw'})#

Every accessor a three-component general modal basis admits, by name.

The sets live here, beside the catalogue whose names they select from, and solprop.output.capability imports them to build the gate: they are a statement about statistics, and keeping them here is what lets solprop.translate.general record the list a file publishes without a translation module importing an output one.

omega0, fwhm_t and fwhm_t_outer are functionals of the same-frequency flux a.H W a and take the basis metric (solprop.stats.METRIC_WEIGHTED_STATISTICS), so they are exact on a declared basis as well as on the identity one; zdw is a property of the channels’ propagation constants and of no basis at all – until a declared basis couples channels whose propagation constants differ, when the block executes through a coupled flow and each channel’s constant is the block’s mean; the route then withholds it (WITHHELD_WITH_A_FLOW).

spectral_flux is the one array of the family – the same flux as a density in frequency rather than integrated. It is admitted but not a default of the route: it is n_omega wide per plane where the other five are scalars, so a run asks for it by name.

Two names of solprop.processing.modal are deliberately absent. quadrature_error is refused separately, by the driver, because what disqualifies it is the absent embedded coarse rule and not the basis – a general guide’s node rule has no embedded pair to difference, so such a run has no transverse-resolution alarm, which is a known limitation – and it sits inside the schema-7 resume key, so publishing it here would promise a statistic the same record refuses. channel_amplitudes takes two arrays and a table rather than an opened handler, so no door ever asks the gate about it; it is in solprop.output.capability._GENERAL_ONLY all the same, so that the one direction that could be reached – an ordinary run asked for it – would be refused for the right reason.

solprop.stats.MAX_APERTURES = 16#: The most apertures one ApertureSet carries.

solprop.stats.MAX_REGIONS = 16#: The most rows StatsSpec.spectral_regions carries.

solprop.stats.METRIC_WEIGHTED_STATISTICS: frozenset[str] = frozenset({'analytic_flux', 'fwhm_t', 'fwhm_t_outer', 'omega0', 'physical_energy', 'spectral_flux'})#

Statistics whose value depends on the basis’s flux metric, by name.

These are the functionals of the physical same-frequency flux sum_mn conj(a_m) W_mn a_n – its integral (physical_energy), its density in frequency (spectral_flux), its peak (analytic_flux), its two widths (fwhm_t, fwhm_t_outer) and the spectral centroid of its per-channel share (omega0). On an orthonormal basis W is the identity and they are the unweighted quantities they have always been, bit for bit; on a general modal basis whose Gram is not the identity they are wrong without it, so a caller that has a metric must hand it over (collect()) and one that cannot must refuse rather than answer (solprop.processing.modal).

zdw is deliberately absent: it is a property of the channels’ propagation constants and of no basis at all, and a unitary change of basis of an exactly degenerate block leaves every one of them unchanged.

solprop.stats.ON_AXIS_PLASMA_STATISTICS = ('on_axis_electron_density', 'on_axis_peak_ionisation_rate', 'on_axis_species_electron_density', 'on_axis_species_peak_ionisation_rate')#

The plasma statistics evaluated on the waveguide’s axis (SEP 19, section 5.2).

Opt-in and saved-plane only. The axis field is synthesised from the mode profiles, as for peak_intensity, and its populations are scanned with exactly the rate tables, preionisation, densities and running-integral rule of the node scan, once per observation; a rule with a node on the axis supplies that row instead. They sit beside the unchanged node reductions – electron_density is the densest node’s – because a few-node field-scale rule may have no node near the axis, and a multimode field’s largest density need not lie on it. Defined on resolved modal routes; the spatial routes, general modal guides, the mode-averaged route and TDSE refuse them by name.

solprop.stats.PLASMA_INTEGRATION_STATISTICS = ('plasma_undershoot', 'plasma_work_defect', 'plasma_boundary_level')#

The field plasma column’s run-time integration diagnostics (SEP 18, section 4.5).

Evaluated with the other plasma statistics at each save plane. Each is a ratio of two maxima over every sample, node and population, so that a node the field barely ionises cannot set it: the undershoot max(-dS) / max(dS) of the exposure increments; the work defect max |paid - freed| / max freed, the loss current’s work per ionisation potential against the fraction it freed; and the boundary level, the free-current integrand’s largest magnitude in the 2K samples the filter reads past the window’s start against its largest anywhere.

solprop.stats.PROFILE_STATISTICS = ('density', 'pressure', 'core_radius')#: Statistics that read a solprop.profiles.ZProfile at the plane’s position.

class solprop.stats.StatDefinition(name, summary, shape, units, requires=(), implemented=True, backend=False, between_planes=False)[source]#

Bases: object

What one statistic is called, what it returns and what it needs.

Attributes:

namestr: Key the statistic is stored under.
summarystr: One line saying what the number is, used by the manual page.
shapestr: Shape of the value in the standard axis names, '()' for a scalar.
unitsstr: SI units of the value.
requirestuple: What collect() must be given for this statistic: 'windows', 'modes', the name of a z-profile, 'spatial_geometry', the far-field tables 'apertures' and 'regions' of the spec, or 'separable_operator', the run’s linear operator, which only a backend holds.
implementedbool: Whether the statistic exists yet. A reserved name is in the catalogue so that a spec asking for it is refused now instead of silently returning nothing.
backendbool: Whether only a backend can compute it. These are the statistics of quantities the saved field does not carry – the electron density, the ionisation rate, the on-axis intensity, the quadrature error – which a backend forms from its own buffers in one pass. collect() refuses them by name rather than pulling the state back across the seam to reconstruct them.
between_planesbool: Whether a backend can also evaluate it between save planes, from what an accepted step already leaves behind, for the periodic channel of StatsSpec.periodic. What disqualifies a statistic is cost per step, not difficulty: zdw needs a root search per mode, and paying that on every step would cost a large share of the step; peak_intensity is in the subset because its on-axis matrix is built once and rescaled rather than rebuilt. StatsSpec refuses a periodic name without it, by name.

Examples

>>> from solprop.stats import CATALOGUE
>>> CATALOGUE["energy"].units, CATALOGUE["energy"].shape
('J', '(n_modes, n_pol); () for spatial beams')
>>> CATALOGUE["electron_density"].backend
True
>>> CATALOGUE["energy"].between_planes, CATALOGUE["zdw"].between_planes
(True, False)

backend: bool = False#

between_planes: bool = False#

implemented: bool = True#

name: str#

requires: tuple[str, ...] = ()#

shape: str#

summary: str#

units: str#

class solprop.stats.StatsSpec(names=(), energy_windows=<factory>, periodic=(), period=1, apertures=None, spectral_regions=<factory>)[source]#

Bases: object

Which statistics a run asks for, and the windows and period they need.

This is the object of the kernel contract §8: the translation layer carries it into a backend, the driver hands it to collect(), and both read the same names. Equality is by identity (eq=False) because energy_windows is an array.

Attributes:

namestuple: Statistics to evaluate at every saved plane, in the order they are evaluated. Every name must be in CATALOGUE and must be implemented.
energy_windowsnumpy.ndarray: Spectral windows, shape (n_windows, 2), rad/s, each row (lo, hi) with lo < hi. Infinite edges are allowed, so that a pair of windows can tile the whole axis. Required by energy_window and peak_power_window.
periodictuple: The subset of names a backend also evaluates between planes, for step-level diagnostics. Must be a subset of names, and every entry must carry StatDefinition.between_planes: a statistic whose cost is only bearable once per plane is refused here rather than made the run’s per-step cost.
periodint: How many accepted steps apart the periodic subset is evaluated, at least 1.
aperturesApertureSet or None: Far-field apertures, required by aperture_spectrum and aperture_on_axis_spectrum (SEP 23). None by default, and then left out of the run’s configuration record, so every earlier record keeps its hash.
spectral_regionsnumpy.ndarray: Boxes of physical transverse wavevector in rad/m, shape (n_regions, 4), rows (kx_lo, kx_hi, ky_lo, ky_hi), each a half-open box [lo, hi) on both axes, at most MAX_REGIONS rows, infinite edges allowed. Required by region_spectrum. Empty by default, and then left out of the configuration record.

Raises:

OutputError: If a name is unknown, reserved or repeated; if a window or a region is malformed; if periodic is not a subset of names or names a statistic no backend evaluates between planes; if period is below 1; or if a name needs a far-field table the spec does not carry.

Examples

>>> import numpy as np
>>> from solprop.stats import StatsSpec
>>> spec = StatsSpec(
...     names=("energy", "energy_window"),
...     energy_windows=np.array([[2.0e15, 3.0e15]]),
... )
>>> spec.n_windows
1

apertures: ApertureSet | None = None#

energy_windows: Float[ndarray, 'n_windows 2']#

property n_apertures: int[source]#

Return how many far-field apertures the spec carries.

Returns:

int: The apertures of apertures, or zero without them.

property n_regions: int[source]#

Return how many regions of transverse wavevector the spec carries.

Returns:

int: The rows of spectral_regions.

property n_windows: int[source]#

Return how many spectral windows the spec carries.

Returns:

int: The number of windows.

names: tuple[str, ...] = ()#

period: int = 1#

periodic: tuple[str, ...] = ()#

same_tables(other)[source]#

Return whether another spec carries the same far-field tables.

Parameters:

otherStatsSpec: The spec to compare, for example a lockstep batch member’s.

Returns:

bool: True when both carry no apertures or the same apertures entry for entry, and identical regions: a structural comparison.

spectral_regions: Float[ndarray, 'n_regions 4']#

solprop.stats.WEIGHT_DEPENDENT_STATISTICS = frozenset({'analytic_flux', 'energy', 'energy_window', 'fwhm_t', 'fwhm_t_outer', 'omega0', 'peak_intensity', 'peak_power', 'peak_power_total', 'peak_power_window', 'physical_energy', 'spectral_flux'})#: Names whose value depends on the modal power weight, so a backend evaluates it.

solprop.stats.WITHHELD_WITH_A_FLOW: frozenset[str] = frozenset({'zdw'})#

The general-basis accessors a route withholds when its block couples, by name.

A declared change of basis over channels whose propagation constants differ executes through linear.coupled_flow, and the per-channel propagation constant the statistics would read is then the block’s mean: zdw from it is a number that propagates nothing. The route leaves these names out of the record’s published list and every accessor door refuses them with solprop.output.capability.refusal_for()’s withheld.

solprop.stats.analytic_field(grid, field)[source]#

Return the time-domain analytic signal of a stored spectrum.

The one construction of the analytic signal in the package, so that the statistics, the backends’ on-axis intensity and anything else that needs a power agree by construction. It serves both grids: an envelope grid’s state already is the analytic signal with the carrier removed, so it only has to be transformed back; a field grid’s is the analytic signal of the real field. solprop.fourier.analytic_signal() neither doubles the zero-frequency bin nor discards the Nyquist bin, because only that convention makes the real part an identity and the two energy functionals an exact Parseval pair.

Parameters:

gridsolprop.grid.TimeGrid: The grid the spectrum lives on.
fieldnumpy.ndarray: Stored spectrum, shape (..., n_omega), internal convention.

Returns:

numpy.ndarray: Complex analytic signal, shape (..., n_t), whose modulus squared is the instantaneous power in watts.

solprop.stats.analytic_flux(grid, field, *, metric=None)[source]#

Return the peak of the plane’s same-frequency analytic flux, in watts.

The general modal basis’s name for a peak power: max_t sum_m |a_m(t)|^2, the same-frequency diagonal of the Poynting flux, which is all a modal state can say about a power on its own.

It is not the cross-section flux, and the difference is physics. The physical int Re(E_a x conj(H_a))/2 . z dA of the reconstructed fields carries the guide’s cross-frequency flux kernel K_ij = sum_q w_q [e_x(w_i) conj(h_y(w_j)) - e_y(w_i) conj(h_x(w_j))]/2, whose diagonal a provider normalises to one and whose off-diagonal it says nothing about. The two peaks therefore differ, bounded term by term by sum_ij |B_i||B_j||K_ij - 1| and measured at 1.3e-5 of the peak on the general-basis test fixture, inside a bound of 3.4e-5: even a positive spectral metric does not make an instantaneous cross-frequency flux a sum of coordinate powers. The time integrals do agree, to 5.0e-7 there, because the cross-frequency terms are an integer number of cycles in the window. For the flux itself, reconstructed from the guide’s stored reference, use solprop.processing.modal.reconstructed_flux(), a trace in time whose peak is this statistic plus the cross-frequency terms; its real-field partner with the sum-frequency term is solprop.processing.modal.instantaneous_flux().

The arithmetic is exactly ``peak_power_total``’s on the identity metric, which is what metric=None evaluates, bit for bit. With a metric it is the same peak taken on the whitened amplitudes c = W^(1/2) a of _whitened, so that the quantity peaked is the physical same-frequency flux sum_mn conj(a_m) W_mn a_n and not a sum of coordinate powers. The old name is refused on this basis rather than redefined: peak_power_total promises that sum of coordinate powers, and on a basis that is not orthonormal there is no such sum to promise.

Parameters:

gridsolprop.grid.TimeGrid: The grid the spectrum lives on.
fieldnumpy.ndarray: Stored spectrum, shape (n_modes, n_pol, n_omega), internal convention, in the square root of watts.
metricnumpy.ndarray or None, optional: The basis flux metric W(omega), shape (n_modes, n_modes, n_omega), dimensionless; None, the default, is the identity.

Returns:

numpy.ndarray: Peak same-frequency analytic flux, a zero-dimensional array, W.

Raises:

solprop.exceptions.OutputError: If metric disagrees in shape with field, or is not positive definite at some frequency, so that no orthonormal basis exists to measure in.

Examples

>>> from solprop.grid import EnvelopeGrid
>>> from solprop.pulses import GaussianPulse
>>> from solprop.stats import analytic_flux, peak_power_total
>>> grid = EnvelopeGrid.from_limits(1.0, 800e-9, (400e-9, 2e-6), 1e-12)
>>> pulse = GaussianPulse(wavelength=800e-9, fwhm=30e-15, peak_power=1e9)
>>> built = pulse.build(grid)
>>> bool(analytic_flux(grid, built) == peak_power_total(grid, built))
True

solprop.stats.check_aperture_choices(profile, weighting)[source]#

Refuse, by name, an aperture profile or weighting that does not exist.

The check ApertureSet makes of its two choices. It takes any object, because a caller outside the test session’s runtime type checker can pass anything.

Parameters:

profileobject: The edge profile asked for.
weightingobject: The weighting asked for.

Raises:

OutputError: If either is not one of its choices.

Examples

>>> from solprop.stats import check_aperture_choices
>>> check_aperture_choices("tanh", "quadrature") is None
True

solprop.stats.collect(spec, grid, plane, *, modes=None, profiles=None, flux=None, metric=None)[source]#

Evaluate every statistic a spec asks for at one saved plane.

Each statistic uses its named public array function and deterministic reductions. For saved true-power modal amplitudes, temporal array functions use |b(t)|². Backend run statistics use analytic Poynting flux: their time integrals agree, but broadband peaks and widths can differ (docs/dev/conventions.md section 4).

Parameters:

specStatsSpec: Which statistics to evaluate, and the windows they need.
gridsolprop.grid.TimeGrid: The grid the plane’s field lives on.
planesolprop.output.schema.SavedPlane: The plane, whose field is in the internal FFT convention.
modescollections.abc.Sequence, optional: The propagated modes, in state order. Needed by zdw.
profilescollections.abc.Mapping, optional: The z-profiles, keyed by statistic name. Needed by density, pressure and core_radius.
fluxsolprop.spatial.SpatialFlux or None, optional: The plane’s measure for a radial or Cartesian spectrum, resolved at the plane’s own position (and jump side) in its medium: every spatial statistic is that medium’s axial Poynting flux (SEP 34), and widths use the signed time-integrated fluence. The far-field statistics aperture_spectrum, aperture_on_axis_spectrum, region_spectrum and beam_spectrum need a Cartesian grid or a line (SEP 23), and are formed in one pass by solprop.spatial.aperture_statistics(), the function a run uses.
metricnumpy.ndarray or None, optional: The basis flux metric W(omega) at this plane, shape (n_modes, n_modes, n_omega), for the names in METRIC_WEIGHTED_STATISTICS; None, the default, is the identity and every statistic is then bitwise what it was. A caller that holds a general modal problem gets it from solprop.translate.coupled.metric_at().

Returns:

dict: The value of each requested statistic, keyed by name, always as a NumPy array and zero-dimensional for a scalar. The shape of each is the shape field of its StatDefinition.

Raises:

OutputError: If a statistic’s dependency is missing – no windows for a windowed statistic, no modes for zdw, or no profile for a profile statistic – or if it is one a backend owns, which this function cannot form from a saved plane (absorbed_power among them, refused by name because a plane does not carry the linear operator’s rates); if a far-field statistic is asked of a radial grid, or of a line with an aperture off y = 0 or a region excluding kappa_y = 0; or if the plane kept no field (SEP 22).

Examples

>>> import numpy as np
>>> from solprop.grid import EnvelopeGrid
>>> from solprop.output.schema import SavedPlane
>>> from solprop.pulses import GaussianPulse
>>> from solprop.stats import collect, default_spec
>>> grid = EnvelopeGrid.from_limits(1.0, 800e-9, (400e-9, 2e-6), 1e-12)
>>> built = GaussianPulse(wavelength=800e-9, fwhm=30e-15, energy=1e-6).build(grid)
>>> plane = SavedPlane(z=0.0, field=built)
>>> values = collect(default_spec("mode_averaged"), grid, plane)
>>> sorted(values)
['energy', 'fwhm_t', 'fwhm_t_outer', 'omega0', 'peak_power']
>>> values["energy"].shape, values["fwhm_t"].shape
((1, 1), ())

solprop.stats.default_spec(kind, *, with_profiles=False, with_zdw=False)[source]#

Return the statistics a run of the given kind records unless told otherwise.

The field statistics every run wants, plus – for a multimode run, where the modes peak at different positions in time – the peak of the total power. The z-profile statistics and the zero-dispersion wavelength are opt-in: the first needs profiles the driver may not have, and the second is a root search per mode per plane.

Parameters:

kind{‘mode_averaged’, ‘multimode’, ‘radial’, ‘free’}: Which transform the run uses.
with_profilesbool, optional: Whether to record the gas density, the pressure and the core radius for modal runs, or only the gas density for spatial runs.
with_zdwbool, optional: Whether to record the zero-dispersion wavelength of each mode. Spatial runs reject this option because they carry no modal dispersion.

Returns:

StatsSpec: The spec, with no energy windows; add those for the run’s bands of interest.

Examples

>>> from solprop.stats import default_spec
>>> default_spec("mode_averaged").names
('omega0', 'energy', 'peak_power', 'fwhm_t', 'fwhm_t_outer')
>>> "peak_power_total" in default_spec("multimode").names
True

solprop.stats.energy(grid, field)[source]#

Return the pulse energy of each mode and polarisation, in joules.

This is solprop.fourier.energy_omega(), the exact discrete form of docs/dev/conventions.md §4, which agrees with the time-domain functional to round-off. No quadrature weights (such as an extended Simpson rule on a field grid) are derived: those would differ from this by their fourth-order quadrature error.

Parameters:

gridsolprop.grid.TimeGrid: The grid the spectrum lives on.
fieldnumpy.ndarray: Stored spectrum, shape (n_modes, n_pol, n_omega), internal convention, in the square root of watts.

Returns:

numpy.ndarray: Energy, shape (n_modes, n_pol), J.

Examples

>>> from solprop.grid import EnvelopeGrid
>>> from solprop.pulses import GaussianPulse
>>> from solprop.stats import energy
>>> grid = EnvelopeGrid.from_limits(1.0, 800e-9, (400e-9, 2e-6), 1e-12)
>>> pulse = GaussianPulse(wavelength=800e-9, fwhm=30e-15, energy=1e-6)
>>> f"{float(energy(grid, pulse.build(grid))[0, 0]):.6e}"
'1.000000e-06'

solprop.stats.energy_window(grid, field, windows)[source]#

Return the energy inside each spectral window, in joules.

The gate is hard and half-open, lo <= omega < hi, so windows that tile the band partition the bins and their energies sum to the total exactly.

Parameters:

gridsolprop.grid.TimeGrid: The grid the spectrum lives on.
fieldnumpy.ndarray: Stored spectrum, shape (n_modes, n_pol, n_omega), internal convention.
windowsnumpy.ndarray: Window edges, shape (n_windows, 2), rad/s, each row (lo, hi). Convert a wavelength pair with omega = 2 pi c / wavelength, remembering that it reverses the order.

Returns:

numpy.ndarray: Energy per window, shape (n_windows, n_modes, n_pol), J.

Examples

>>> import numpy as np
>>> from solprop.grid import EnvelopeGrid
>>> from solprop.pulses import GaussianPulse
>>> from solprop.stats import energy, energy_window
>>> grid = EnvelopeGrid.from_limits(1.0, 800e-9, (400e-9, 2e-6), 1e-12)
>>> built = GaussianPulse(wavelength=800e-9, fwhm=30e-15, energy=1e-6).build(grid)
>>> tiling = np.array([[-np.inf, 2.4e15], [2.4e15, np.inf]])
>>> total = energy_window(grid, built, tiling).sum()
>>> bool(np.isclose(total, energy(grid, built).sum(), rtol=1e-15))
True

solprop.stats.fwhm_t(grid, field, *, metric=None)[source]#

Return the width of the central feature at half maximum, in seconds.

The total power sum_modes sum_pol |E_a(t)|^2 is formed, the half-maximum level taken from its global peak with no baseline subtraction, and the crossings found by walking outwards from that peak to the first sample below the level on either side and interpolating linearly. Structure further out – a dispersive wave, a soliton that has walked off, a noise pedestal – does not widen the answer, which is why this is the number to read as a pulse duration. Use fwhm_t_outer() for the total extent above half maximum.

On a basis whose flux Gram is not the identity this takes the metric: the trace whose width is measured is then the physical same-frequency flux analytic_flux() peaks, evaluated on the whitened amplitudes of _whitened, so that a run does not report the width of a quantity whose peak it does not report.

Parameters:

gridsolprop.grid.TimeGrid: The grid the spectrum lives on.
fieldnumpy.ndarray: Stored spectrum, shape (n_modes, n_pol, n_omega), internal convention.
metricnumpy.ndarray or None, optional: The basis flux metric W(omega), shape (n_modes, n_modes, n_omega), dimensionless; None, the default, is the identity.

Returns:

numpy.ndarray: The width as a zero-dimensional array, s. A trace that never falls to half its maximum – a continuous wave, an empty plane – has no width and gives nan.

Raises:

solprop.exceptions.OutputError: If metric disagrees in shape with field, or is not positive definite at some frequency, so that no orthonormal basis exists to measure in.

Examples

>>> import numpy as np
>>> from solprop.grid import EnvelopeGrid
>>> from solprop.pulses import GaussianPulse
>>> from solprop.stats import fwhm_t
>>> grid = EnvelopeGrid.from_limits(1.0, 800e-9, (400e-9, 2e-6), 1e-12)
>>> built = GaussianPulse(wavelength=800e-9, fwhm=30e-15, energy=1e-6).build(grid)
>>> bool(np.isclose(fwhm_t(grid, built), 30e-15, rtol=1e-4))
True

solprop.stats.fwhm_t_outer(grid, field, *, metric=None)[source]#

Return the total extent above half maximum, in seconds.

The first and the last sample anywhere in the window that exceed half the peak power, interpolated outwards. On a clean single pulse this equals fwhm_t(); on a structured trace it spans every satellite, so the pair together says how much of the energy sits outside the central feature.

It takes the metric for the reason fwhm_t() does.

Parameters:

gridsolprop.grid.TimeGrid: The grid the spectrum lives on.
fieldnumpy.ndarray: Stored spectrum, shape (n_modes, n_pol, n_omega), internal convention.
metricnumpy.ndarray or None, optional: The basis flux metric W(omega), shape (n_modes, n_modes, n_omega), dimensionless; None, the default, is the identity.

Returns:

numpy.ndarray: The extent as a zero-dimensional array, s, or nan when the trace reaches the edge of the window still above half its maximum, or carries no power.

Raises:

solprop.exceptions.OutputError: If metric disagrees in shape with field, or is not positive definite at some frequency, so that no orthonormal basis exists to measure in.

Examples

>>> import numpy as np
>>> from solprop.grid import EnvelopeGrid
>>> from solprop.pulses import GaussianPulse
>>> from solprop.stats import fwhm_t, fwhm_t_outer
>>> grid = EnvelopeGrid.from_limits(1.0, 800e-9, (400e-9, 2e-6), 1e-12)
>>> built = GaussianPulse(wavelength=800e-9, fwhm=30e-15, energy=1e-6).build(grid)
>>> bool(np.isclose(fwhm_t_outer(grid, built), fwhm_t(grid, built), rtol=1e-12))
True

solprop.stats.inband_wavelengths(grid)[source]#

Return the shortest and longest in-band wavelength of a grid, in metres.

The search window zdw() uses by default when collect() calls it: a zero-dispersion point outside the simulated band says nothing about the run, and hunting for one there drives every material model into extrapolation.

Parameters:

gridsolprop.grid.TimeGrid: The grid.

Returns:

tuple of float: (shortest, longest) in m, from the grid’s in-band frequency limits.

Examples

>>> from solprop.grid import EnvelopeGrid
>>> from solprop.stats import inband_wavelengths
>>> grid = EnvelopeGrid.from_limits(1.0, 800e-9, (400e-9, 2e-6), 1e-12)
>>> [f"{value * 1e9:.1f}" for value in inband_wavelengths(grid)]
['400.0', '2000.0']

solprop.stats.omega0(grid, field, *, metric=None)[source]#

Return the centroid of the spectral energy density, per mode and polarisation.

\[\omega_0 = \frac{\sum_k \omega_k\,S(\omega_k)}{\sum_k S(\omega_k)}\]

with \(S\) from solprop.fourier.spectral_density(). The prefactor of the density cancels in the ratio, so this is the first moment of \(|\tilde E|^2\) over the frequency axis. The sum runs over every bin, including the zero-frequency one: outside the active band a propagated state is identically zero, so the bins there contribute nothing.

Because the centroid is an energy-weighted mean, the whole field’s centroid is the energy-weighted mean of these, which is why the statistic is per column rather than summed.

On a basis whose flux Gram is not the identity this takes the metric, and the columns are then the whitened ones of _whitened: a coordinate column of a nonorthogonal basis carries no separable spectral energy at all – the columns interfere – so a centroid of |a_m(omega)|^2 would be a promise the basis cannot keep. W^(1/2) is the identity plus half the basis’s Gram departure, so on a unitary change of basis of a certified orthonormal reference, which is the only kind this package executes, the whitened column is the declared channel to within that departure and the centroid moves by no more.

An empty column is empty relative to the plane. A column whose density sums to at most one machine epsilon of the largest column’s has no centroid and is nan. That floor is relative rather than exactly zero because whitening mixes: the unlaunched partner of a launched circular channel carries W^(1/2) times the other’s amplitude – 1e-32 of it on this route’s own fixture – and an exact-zero test would report a finite centroid of the whitening’s rounding noise for a channel nothing was launched into.

Parameters:

gridsolprop.grid.TimeGrid: The grid the spectrum lives on.
fieldnumpy.ndarray: Stored spectrum, shape (n_modes, n_pol, n_omega), internal convention.
metricnumpy.ndarray or None, optional: The basis flux metric W(omega), shape (n_modes, n_modes, n_omega), dimensionless; None, the default, is the identity and is the bitwise unchanged path.

Returns:

numpy.ndarray: Centroid frequency, shape (n_modes, n_pol), rad/s. A column carrying no energy has no centroid and is nan.

Raises:

solprop.exceptions.OutputError: If metric disagrees in shape with field, or is not positive definite at some frequency, so that no orthonormal basis exists to measure in.

Examples

>>> from solprop.grid import EnvelopeGrid
>>> from solprop.pulses import GaussianPulse
>>> from solprop.stats import omega0
>>> grid = EnvelopeGrid.from_limits(1.0, 800e-9, (400e-9, 2e-6), 1e-12)
>>> pulse = GaussianPulse(wavelength=800e-9, fwhm=30e-15, energy=1e-6)
>>> f"{float(omega0(grid, pulse.build(grid))[0, 0]):.6e}"
'2.354564e+15'

solprop.stats.peak_power(grid, field)[source]#

Return the largest instantaneous power each mode reaches, in watts.

max_t sum_pol |E_a(t)|^2: the polarisation components of one mode are summed before the maximum is taken, because they are two components of one field and their powers add at each instant. The analytic signal is formed once for the whole call. This array function measures squared true-power amplitude; run statistics use analytic Poynting flux, with the arithmetic-mean frequency kernel instead.

Parameters:

gridsolprop.grid.TimeGrid: The grid the spectrum lives on.
fieldnumpy.ndarray: Stored spectrum, shape (n_modes, n_pol, n_omega), internal convention.

Returns:

numpy.ndarray: Peak power, shape (n_modes,), W.

Examples

>>> from solprop.grid import EnvelopeGrid
>>> from solprop.pulses import GaussianPulse
>>> from solprop.stats import peak_power
>>> grid = EnvelopeGrid.from_limits(1.0, 800e-9, (400e-9, 2e-6), 1e-12)
>>> pulse = GaussianPulse(wavelength=800e-9, fwhm=30e-15, peak_power=1e9)
>>> f"{float(peak_power(grid, pulse.build(grid))[0]):.6e}"
'1.000000e+09'

solprop.stats.peak_power_total(grid, field)[source]#

Return the largest instantaneous power of the whole field, in watts.

max_t sum_modes sum_pol |E_a(t)|^2: the modes are summed before the maximum, so this is the peak of the total power and not the largest of the per-mode peaks, which arrive at different times.

Parameters:

gridsolprop.grid.TimeGrid: The grid the spectrum lives on.
fieldnumpy.ndarray: Stored spectrum, shape (n_modes, n_pol, n_omega), internal convention.

Returns:

numpy.ndarray: Peak total power, a zero-dimensional array, W.

Examples

>>> import numpy as np
>>> from solprop.grid import EnvelopeGrid
>>> from solprop.pulses import GaussianPulse
>>> from solprop.stats import peak_power, peak_power_total
>>> grid = EnvelopeGrid.from_limits(1.0, 800e-9, (400e-9, 2e-6), 1e-12)
>>> pulse = GaussianPulse(wavelength=800e-9, fwhm=30e-15, peak_power=1e9)
>>> built = pulse.build(grid)
>>> f"{float(peak_power_total(grid, built)):.6e}"
'1.000000e+09'

solprop.stats.peak_power_window(grid, field, windows)[source]#

Return the peak power of each mode inside each spectral window, in watts.

The spectrum is gated with the same hard half-open window as energy_window() and transformed back, so this is the peak power of the light in that band. One scalar per mode per window.

Parameters:

gridsolprop.grid.TimeGrid: The grid the spectrum lives on.
fieldnumpy.ndarray: Stored spectrum, shape (n_modes, n_pol, n_omega), internal convention.
windowsnumpy.ndarray: Window edges, shape (n_windows, 2), rad/s, each row (lo, hi).

Returns:

numpy.ndarray: Peak power, shape (n_windows, n_modes), W.

Examples

>>> import numpy as np
>>> from solprop.grid import EnvelopeGrid
>>> from solprop.pulses import GaussianPulse
>>> from solprop.stats import peak_power_window
>>> grid = EnvelopeGrid.from_limits(1.0, 800e-9, (400e-9, 2e-6), 1e-12)
>>> pulse = GaussianPulse(wavelength=800e-9, fwhm=30e-15, peak_power=1e9)
>>> built = pulse.build(grid)
>>> whole = np.array([[-np.inf, np.inf]])
>>> f"{float(peak_power_window(grid, built, whole)[0, 0]):.6e}"
'1.000000e+09'

solprop.stats.physical_energy(grid, field, *, metric=None)[source]#

Return the physical energy through the plane, in joules.

The general modal basis’s name for an energy, and the one number of that family that survives a change of coordinates: the same-frequency product of the state with the basis’s flux metric, sum_k sum_mn conj(a_m) W_mn a_n, integrated over frequency. A coordinate energy is not additive when W is not the identity, and this total is, which is why it is a scalar and energy is an array.

On the identity metric this is exactly ``sum(energy(grid, field))``, which is what metric=None evaluates, bit for bit. With a metric it is the same functional of the whitened amplitudes c = W^(1/2) a of _whitened, which is the definition above and not a second one: sum_m |c_m|^2 = sum_mn conj(a_m) W_mn a_n at every frequency, and the frequency integral is the Parseval weight of docs/dev/conventions.md section 4 either way.

Parameters:

gridsolprop.grid.TimeGrid: The grid the spectrum lives on.
fieldnumpy.ndarray: Stored spectrum, shape (n_modes, n_pol, n_omega), internal convention, in the square root of watts.
metricnumpy.ndarray or None, optional: The basis flux metric W(omega), shape (n_modes, n_modes, n_omega), dimensionless; None, the default, is the identity.

Returns:

numpy.ndarray: Physical energy, a zero-dimensional array, J.

Raises:

solprop.exceptions.OutputError: If metric disagrees in shape with field, or is not positive definite at some frequency, so that no orthonormal basis exists to measure in.

Examples

>>> import numpy as np
>>> from solprop.grid import EnvelopeGrid
>>> from solprop.pulses import GaussianPulse
>>> from solprop.stats import energy, physical_energy
>>> grid = EnvelopeGrid.from_limits(1.0, 800e-9, (400e-9, 2e-6), 1e-12)
>>> built = GaussianPulse(wavelength=800e-9, fwhm=30e-15, energy=1e-6).build(grid)
>>> bool(physical_energy(grid, built) == energy(grid, built).sum())
True

solprop.stats.profile_value(profile, z)[source]#

Return the value of a z-profile at one position.

The one implementation behind the density, pressure and core_radius statistics, which differ only in which profile they read; the catalogue gives each its own units.

Parameters:

profilesolprop.profiles.ZProfile: The path along the waveguide.
zfloat: Position, m.

Returns:

numpy.ndarray: The value as a zero-dimensional array, in the profile’s own units. Those are SI except for a pressure path, which is in bar: pressures are the one documented exception at the user interface (docs/dev/conventions.md §1) and a solprop.profiles.PressureGradient is a user-interface object.

Examples

>>> from solprop.profiles import PressureGradient
>>> from solprop.stats import profile_value
>>> gradient = PressureGradient.two_point(1.0, 2.0, 0.0)
>>> f"{float(profile_value(gradient, 0.75)):.6f}"
'1.000000'

solprop.stats.species_statistics(densities, fractions, peak_rates, potentials, weights)[source]#

Reduce species populations at the node of largest total electron density.

Parameters:

densitiesnumpy.ndarray: Component number densities in m^-3, in mixture order.
fractionsnumpy.ndarray: Final ionised fractions at each node, shape (n_species, n_nodes). Components without an ionisation channel have zero entries.
peak_ratesnumpy.ndarray: Largest rate over all nodes per species in s^-1; cycle-averaged on envelope grids.
potentialsnumpy.ndarray: Ionisation potentials in J, zero for components without channels.
weightsnumpy.ndarray: Physical transverse quadrature weights in m². The mode-averaged case uses its effective area as its single weight.

Returns:

dict: The five species statistics plus total electron density and peak rate. species_ionisation_energy is final ionised inventory in J/m, including preionisation. Pulse-deposited ionisation work counts only newly freed electrons and must subtract the initial fractions before this reduction.

Raises:

OutputError: If a population, rate, potential or quadrature weight is outside its domain.

solprop.stats.spectral_flux(grid, field, *, metric=None)[source]#

Return the physical energy’s spectral density, in J/(rad/s).

The spectrum analogue of physical_energy(), and the one array of the general modal family: the same same-frequency product of the state with the basis’s flux metric, sum_mn conj(a_m) W_mn a_n, published per bin instead of integrated over frequency. It is a density in the sense of docs/dev/conventions.md section 4 – dU/domega, the prefactor that pairs with the shot-noise amplitude so that one photon in a bin gives S(omega_k) domega = hbar omega_k – and it is solprop.fourier.spectral_density() summed over the channel axis, which is what makes it additive where a coordinate spectrum is not.

It is indexed by frequency and by nothing else. A per-channel spectral energy is no more additive than a per-channel energy on a basis whose Gram is not the identity, so the channel axis is summed after the whitening and not published.

The frequency sum is ``physical_energy``, exactly on every interior bin: domega times this density is the 4 dt / n_t weight of solprop.fourier.energy_omega() there. The zero-frequency and Nyquist bins are the documented exception – the analytic signal does not double them, so summing the density over every bin overcounts those two – and both are identically zero on any grid solprop.grid builds, whose spectral collar apodises them away. Ask physical_energy() when an exact energy is wanted.

Parameters:

gridsolprop.grid.TimeGrid: The grid the spectrum lives on.
fieldnumpy.ndarray: Stored spectrum, shape (n_modes, n_pol, n_omega), internal convention, in the square root of watts.
metricnumpy.ndarray or None, optional: The basis flux metric W(omega), shape (n_modes, n_modes, n_omega), dimensionless; None, the default, is the identity.

Returns:

numpy.ndarray: Spectral energy density, shape (n_omega,), J/(rad/s).

Raises:

solprop.exceptions.OutputError: If metric disagrees in shape with field, or is not positive definite at some frequency, so that no orthonormal basis exists to measure in.

Examples

>>> import numpy as np
>>> from solprop.grid import EnvelopeGrid
>>> from solprop.pulses import GaussianPulse
>>> from solprop.stats import physical_energy, spectral_flux
>>> grid = EnvelopeGrid.from_limits(1.0, 800e-9, (400e-9, 2e-6), 1e-12)
>>> built = GaussianPulse(wavelength=800e-9, fwhm=30e-15, energy=1e-6).build(grid)
>>> density = spectral_flux(grid, built)
>>> density.shape == (grid.n_omega,)
True
>>> total = float(np.sum(density) * grid.d_omega)
>>> bool(np.isclose(total, float(physical_energy(grid, built)), rtol=1e-14))
True

solprop.stats.weight_refusal_for(names, weight)[source]#

Return why a general statistic cannot be formed under this weight, or None.

physical_energy(), spectral_flux() and analytic_flux() evaluate the identity-metric, unit-weight arithmetic, which is what the general modal route certifies before it will run (linear.power_weight is the constant one there). This is unreachable today: no route carries both a modal table and a non-constant weight, and the translation layer refuses to build one. It is written because the failure would be silent rather than loud if a later change loosened either side – the three names are outside the set a backend routes through the weighted Poynting pass, so they would be answered from the unweighted field.

Parameters:

namescollections.abc.Sequence: The statistic names asked for.
weightnumpy.ndarray: The modal power weight at this plane, shape (n_modes, n_pol, n_omega).

Returns:

str or None: The refusal, or None when no general statistic is asked for or the weight is the constant one.

Examples

>>> import numpy as np
>>> from solprop.stats import weight_refusal_for
>>> weight_refusal_for(["physical_energy"], np.ones((1, 1, 4))) is None
True
>>> weight_refusal_for(["energy"], np.full((1, 1, 4), 2.0)) is None
True

solprop.stats.zdw(modes, z, *, wavelength_range=(1e-07, 3e-06))[source]#

Return the zero-dispersion wavelength of each mode at one position.

Off by default: it is a root search on beta_2 per mode, which is slow and does not belong on every accepted step. Each call stands alone: no previous answer is carried over as the next search’s starting guess, so one failure cannot poison every later plane with nan.

Parameters:

modescollections.abc.Sequence: The propagated modes, in state order.
zfloat: Position along the waveguide, m.
wavelength_rangetuple of float, optional: Shortest and longest wavelength in m to bracket the root in. Defaults to a 100 nm to 3 um window; collect() narrows it to the grid’s own band with inband_wavelengths(), which is where the answer means anything and where the material models are not extrapolating.

Returns:

numpy.ndarray: Zero-dispersion wavelength, shape (n_modes,), m.

Raises:

solprop.exceptions.ModeError: If a mode has no zero-dispersion point in the search window, rather than a missing value that would propagate silently.

Examples

>>> from solprop.modes import MarcatiliMode
>>> from solprop.stats import zdw
>>> mode = MarcatiliMode(125e-6, "He", 0.4, model="reduced", cladding=1.45)
>>> f"{zdw([mode], 0.0, wavelength_range=(150e-9, 900e-9))[0] * 1e9:.1f}"
'376.7'