solprop.processing.fields#

Read fields, spectra and spectral densities out of a saved run.

The field, spectrum and density accessors, in solprop’s conventions. Each one takes an Output – a memory handler or an HDF5 handler, through the same protocol – selects a saved plane and a set of modes, and returns an abscissa and an array whose trivial mode and polarisation axes have been squeezed away.

Three rules hold throughout.

Convention. Arrays come back in the optics convention by default and in the FFT convention on request. For a spectrum the difference is the conjugation and the phase origin, both of which the output layer applies; for a time-domain array it is the conjugation alone. A real field, an energy and a density are convention-free.

Normalisation. get_spectrum() carries the spectral normalisation, so that abs(spectrum) ** 2 is the spectral energy density in J/(rad/s) – solprop.fourier.spectral_density() of the same plane, to round-off. The stored array itself, the raw forward DFT, is output.fields().

Band. The default is the grid’s active band, on an ascending frequency axis. Pass band="all" for every stored bin, remembering that a field grid then includes the zero-frequency bin, where a wavelength is undefined.

No window is ever applied implicitly. bandpass is a keyword, it is a wavelength range, and nothing happens without it.

A spatial file’s columns are not modes. get_power() and get_spectral_density() read each column of the mode axis as a modal amplitude in sqrt(W), so on a radial, Cartesian, line or planar file they refuse by name and name the functional of solprop.spatial that measures the quantity. get_spectrum(), get_field() and get_phase() return the stored columns as they are, in the file’s own units, and read every file.

The noise floor. A run with input shot noise saves a field that carries its seed, one photon per active bin and channel. noise_floor="subtract" on get_spectral_density() removes exactly that seed’s energy, bin by bin and channel by channel, and returns the signed remainder (SEP 37, section 5.5).

Module Attributes

`FieldKind`	Which time-domain representation `get_field()` returns.
`FIELD_KINDS`	The analytic signal, the real field, and the analytic signal without its carrier.
`NoiseFloor`	Whether a spectral density keeps or removes the input noise's one-photon seed.
`NOISE_FLOORS`	The saved density as it stands, and the density less the seed's energy.
`DEFAULT_GATE_POSITIONS`	How many gate positions a spectrogram uses when it is not told.

Functions

`get_field`(output, *[, kind, plane, z, ...])	Return the time-domain field of one plane, or of every plane.
`get_phase`(output, *[, plane, z, modes, band])	Return the unwrapped spectral phase of one plane, or of every plane.
`get_power`(output, *[, plane, z, modes, ...])	Return the instantaneous power of one plane, or of every plane.
`get_spectral_density`(output[, axis, plane, ...])	Return the spectral energy density of one plane, or of every plane.
`get_spectrogram`(output[, axis, plane, z, ...])	Return the time-frequency map of one saved plane.
`get_spectrum`(output, *[, plane, z, modes, ...])	Return the complex spectrum of one plane, or of every plane.
`nearest_plane`(output, z)	Return the index of the saved plane nearest a position.

solprop.processing.fields.DEFAULT_GATE_POSITIONS = 256#: How many gate positions a spectrogram uses when it is not told.

solprop.processing.fields.FIELD_KINDS: tuple[Literal['analytic', 'real', 'envelope'], ...] = ('analytic', 'real', 'envelope')#: The analytic signal, the real field, and the analytic signal without its carrier.

solprop.processing.fields.FieldKind#

Which time-domain representation get_field() returns.

alias of Literal[‘analytic’, ‘real’, ‘envelope’]

solprop.processing.fields.NOISE_FLOORS: tuple[Literal['keep', 'subtract'], ...] = ('keep', 'subtract')#: The saved density as it stands, and the density less the seed’s energy.

solprop.processing.fields.NoiseFloor#

Whether a spectral density keeps or removes the input noise’s one-photon seed.

alias of Literal[‘keep’, ‘subtract’]

solprop.processing.fields.get_field(output, *, kind='analytic', plane=None, z=None, modes=None, oversampling=1, bandpass=None, convention='optics', transfer=None, transform_limit=False)[source]#

Return the time-domain field of one plane, or of every plane.

On either grid and in three representations: the analytic signal E_a(t), whose modulus squared is the instantaneous power in watts; the real field Re E_a; and the envelope A(t) = E_a(t) exp(-i omega0 t). On an envelope grid the stored state already is the envelope, and on a field grid the analytic signal is formed by solprop.fourier.analytic_signal(); one definition serves both.

Oversampling zero-pads the spectrum, so the returned time axis is the grid’s own fine axis: centred on zero, of length oversampling * n_t.

Parameters:

outputsolprop.output.schema.Output: An opened output handler.
kind{‘analytic’, ‘real’, ‘envelope’}: Which representation to return.
planeint or None: Index of the plane wanted, negative from the end; see get_spectrum().
zfloat or None: Position of the plane wanted, m.
modesint, sequence of int, ‘sum’ or None: Which modes to return. 'sum' is refused: modal amplitudes are complex.
oversamplingint: Factor by which to resample the time axis, at least 1.
bandpasstuple of float, Passband, numpy.ndarray or None: Wavelength limits, m, or a Passband on any axis, of an explicit Planck-tapered window applied to the spectrum before the transform, or a ready-made window on the stored axis.
convention{‘optics’, ‘fft’}: Sign convention of the complex representations; 'real' is convention-free.
transferTransferFunction or None: A linear optical element of solprop.optics applied to the spectrum before the transform – a compressor, a window, a stretch of air – in the optics convention.
transform_limitbool: Whether to discard the spectral phase first, giving the shortest pulse the spectrum supports, centred on t = 0. Applied after bandpass and before transfer.

Returns:

tnumpy.ndarray: Time axis, shape (n_t,) or (oversampling n_t,), s, centred on zero.
fieldnumpy.ndarray: The field, in the square root of watts, with trivial axes squeezed; real for kind='real' and complex otherwise.

Raises:

solprop.exceptions.ProcessingError: If kind is unknown or oversampling is below one.
solprop.exceptions.ModeError: If modes='sum'.

Examples

>>> import numpy as np
>>> from solprop.processing._examples import example_output
>>> from solprop.processing.fields import get_field
>>> t, field = get_field(example_output(), plane=0, kind="real")
>>> bool(np.isrealobj(field))
True

solprop.processing.fields.get_phase(output, *, plane=None, z=None, modes=None, band='active')[source]#

Return the unwrapped spectral phase of one plane, or of every plane.

The phase of the optics-convention spectrum, unwrapped along the frequency axis. The optics convention has already moved the phase origin from the first time sample to the centre of the window (conventions.md section 2), so no linear ramp omega T / 2 is subtracted and a flat phase means a transform-limited pulse.

Parameters:

outputsolprop.output.schema.Output: An opened output handler.
planeint or None: Index of the plane wanted, negative from the end; see get_spectrum().
zfloat or None: Position of the plane wanted, m.
modesint, sequence of int, ‘sum’ or None: Which modes to return. 'sum' is refused: a phase is not additive.
band{‘active’, ‘all’}: Which part of the stored axis to use.

Returns:

omeganumpy.ndarray: Angular frequency, shape (n_x,), rad/s, ascending.
phasenumpy.ndarray: Unwrapped spectral phase, rad, with trivial axes squeezed.

Examples

>>> import numpy as np
>>> from solprop.processing._examples import example_output
>>> from solprop.processing.fields import get_phase
>>> omega, phase = get_phase(example_output(), plane=0)
>>> bool(np.all(np.isfinite(phase)))
True

solprop.processing.fields.get_power(output, *, plane=None, z=None, modes=None, oversampling=1, bandpass=None, transfer=None, transform_limit=False)[source]#

Return the instantaneous power of one plane, or of every plane.

P(t) = sum_pol abs(E_a(t)) ** 2 in watts, the quantity every temporal metric and every temporal plot is built on. Polarisation components are summed because they are two components of one field whose powers add at each instant, which is the reduction solprop.stats.peak_power() uses.

Parameters:

outputsolprop.output.schema.Output: An opened output handler.
planeint or None: Index of the plane wanted, negative from the end; see get_spectrum().
zfloat or None: Position of the plane wanted, m.
modesint, sequence of int, ‘sum’ or None: Which modes to return; 'sum' adds their powers.
oversamplingint: Factor by which to resample the time axis, at least 1.
bandpasstuple of float, Passband, numpy.ndarray or None: Wavelength limits, m, or a Passband on any axis, of an explicit window applied before the transform, or a ready-made window on the stored axis.
transferTransferFunction or None: A linear optical element of solprop.optics applied to the spectrum before the transform – a compressor, a window, a stretch of air – in the optics convention.
transform_limitbool: Whether to discard the spectral phase first, giving the shortest pulse the spectrum supports, centred on t = 0. Applied after bandpass and before transfer.

Returns:

tnumpy.ndarray: Time axis, shape (n_t,) or (oversampling n_t,), s.
powernumpy.ndarray: Instantaneous power, W, with trivial axes squeezed.

Raises:

solprop.exceptions.OutputError: If the stored run is a spatial one, whose columns are not modal amplitudes: solprop.spatial.spatial_power() of a saved plane is the beam’s power.

Examples

>>> import numpy as np
>>> from solprop.processing._examples import example_output
>>> from solprop.processing.fields import get_power
>>> t, power = get_power(example_output(), plane=0)
>>> bool(np.all(power >= 0.0))
True

solprop.processing.fields.get_spectral_density(output, axis='wavelength', *, plane=None, z=None, modes=None, band='active', spec_range=None, resolution=None, noise_floor='keep')[source]#

Return the spectral energy density of one plane, or of every plane.

The density is per unit of the abscissa asked for: J/(rad/s), J/Hz or J/m, with the Jacobian of solprop.processing.axes.spectral_axis() applied. The wavelength abscissa is returned ascending, so the density is reversed with it. Polarisation components are summed, because they are two components of one field and their densities add.

Parameters:

outputsolprop.output.schema.Output: An opened output handler.
axis{‘wavelength’, ‘f’, ‘omega’}: The abscissa, in m, Hz or rad/s.
planeint or None: Index of the plane wanted, negative from the end; see get_spectrum().
zfloat or None: Position of the plane wanted, m.
modesint, sequence of int, ‘sum’ or None: Which modes to return; 'sum' adds their densities.
band{‘active’, ‘all’}: Which part of the stored axis to use.
spec_rangetuple of float or None: Limits on the abscissa, in its own units, applied after the change of variable.
resolutionfloat or None: Instrument resolution, in the units of the abscissa, as a full width at half maximum. When given, the density is resampled onto a uniform axis and smoothed by solprop.processing.spectra.specres().
noise_floor{‘keep’, ‘subtract’}: 'keep' returns the saved density as it stands. 'subtract' removes the input noise’s seed: the energy of the stored realisation, which is one photon per active bin and channel in the route’s own photon measure (the one solprop.noise.photons_per_bin() inverts), before the polarisation components are summed and any mode selection or smoothing. The draw has a fixed modulus, so the subtraction is exact for the one realisation; where the linear operator only changes the phase, it equals the modified model’s corrected readout |psi + N|^2 - |N|^2 of the same draw. It is not a photon estimator under loss, where the seed has decayed and more is removed than remains. Values near the floor can be negative, by the interference of the pulse with the seed and by its decay, and are returned signed. Refused for a run without noise and for the 'modified' model, whose saved field holds no seed.

Returns:

xnumpy.ndarray: The abscissa, shape (n_x,), ascending.
densitynumpy.ndarray: The spectral energy density, J per unit of x; signed under noise_floor='subtract'.

Raises:

solprop.exceptions.ProcessingError: If noise_floor is not one of NOISE_FLOORS, or is 'subtract' for a run without input noise, or if spec_range selects no sample of the abscissa (a range in another axis’s units, say).
solprop.exceptions.OutputError: If the stored run is a spatial one, whose columns are not modal amplitudes.

Examples

>>> import numpy as np
>>> from solprop.processing._examples import example_output
>>> from solprop.processing.fields import get_spectral_density
>>> wavelength, density = get_spectral_density(example_output(), plane=0)
>>> bool(np.all(np.diff(wavelength) > 0))
True

solprop.processing.fields.get_spectrogram(output, axis='f', *, plane=None, z=None, modes=None, time_range=None, gate_fwhm=None, positions=256, bandpass=None, transfer=None)[source]#

Return the time-frequency map of one saved plane.

The real field (get_field()) is gated by a Gaussian of the given full width at half maximum and transformed at evenly spaced positions (scipy.signal.ShortTimeFFT), which trades time resolution against frequency resolution in the usual way: a gate of duration tau resolves frequencies no better than 1/tau. The map is abs(stft) ** 2, in arbitrary units, per unit of the abscissa asked for: the change of variable carries the Jacobian of solprop.processing.axes.spectral_axis(), as get_spectral_density() does, so a map on wavelength is a density per metre.

Parameters:

outputsolprop.output.schema.Output: An opened output handler.
axis{‘f’, ‘omega’, ‘wavelength’}: The spectral abscissa, in Hz, rad/s or m.
planeint or None: Index of the plane, negative from the end; the last plane when neither it nor z is given.
zfloat or None: Position of the plane, m: the nearest saved plane.
modesint, sequence of int, ‘sum’ or None: Which mode. A multimode run must name one: a spectrogram is of a field, and modal fields do not add away from a point in space.
time_rangetuple of float or None: The time range the gate positions cover, s, in either order and clipped to the time window; the whole window by default.
gate_fwhmfloat or None: Width of the gate, s; a twentieth of the time range by default.
positionsint: How many gate positions to place across the time range.
bandpasstuple of float, Passband, numpy.ndarray or None: A window applied to the spectrum before the transform, as get_field() takes it: the map of the light in that band.
transferTransferFunction or None: A linear optical element of solprop.optics applied first, so that the map is of the pulse after it.

Returns:

tnumpy.ndarray: The gate positions within the time range, s, shape (n_slices,).
xnumpy.ndarray: The positive frequencies of the transform on axis, ascending, shape (n_x,).
valuesnumpy.ndarray: The map, shape (n_x, n_slices), arbitrary units per unit of x.

Raises:

solprop.exceptions.ModeError: If a multimode run does not name a single mode.
solprop.exceptions.ProcessingError: If the time range misses the time window or selects no gate position, the gate is not positive or is wider than the window, the number of positions is not positive, or axis is unknown.

Examples

>>> import numpy as np
>>> from solprop.processing._examples import example_output
>>> from solprop.processing.fields import get_spectrogram
>>> t, f, values = get_spectrogram(example_output(), positions=32)
>>> values.shape == (f.size, t.size), bool(np.all(np.diff(f) > 0))
(True, True)

solprop.processing.fields.get_spectrum(output, *, plane=None, z=None, modes=None, band='active', convention='optics', noise_floor='keep')[source]#

Return the complex spectrum of one plane, or of every plane.

The stored raw forward DFT is scaled by the spectral normalisation – 2 dt / sqrt(2 pi) on a field grid and dt / sqrt(2 pi) on an envelope grid – so that the modulus squared is the spectral energy density in J/(rad/s) (conventions.md section 4). The axis is ascending and, by default, the active band.

Parameters:

outputsolprop.output.schema.Output: An opened output handler.
planeint or None: Index of the plane wanted, negative from the end. With z, at most one of the two; with neither, every plane is returned and the result carries a leading n_z axis.
zfloat or None: Position of the plane wanted, m; the nearest saved plane is used.
modesint, sequence of int, ‘sum’ or None: Which modes to return. 'sum' is refused here: modal amplitudes are complex and adding them is meaningless away from a point in space.
band{‘active’, ‘all’}: Which part of the stored axis to return.
convention{‘optics’, ‘fft’}: Sign convention of the returned spectrum.
noise_floor{‘keep’}: Only 'keep'. The noise floor is a power, one photon per bin, and an amplitude has no part that is that power: the saved amplitude is the pulse and the seed added coherently, and their interference is as large as the seed itself, so no amplitude has the floor-subtracted density as its modulus squared, which can be negative. 'subtract' is refused and named; get_spectral_density() subtracts it.

Returns:

omeganumpy.ndarray: Angular frequency, shape (n_x,), rad/s, ascending.
spectrumnumpy.ndarray: The spectrum, in the square root of J/(rad/s), with the leading n_z axis when no plane was selected and trivial mode and polarisation axes squeezed.

Raises:

solprop.exceptions.ModeError: If modes='sum'.
solprop.exceptions.ProcessingError: If noise_floor is not 'keep'.

Examples

>>> import numpy as np
>>> from solprop.processing._examples import example_output
>>> from solprop.processing.fields import get_spectrum
>>> omega, spectrum = get_spectrum(example_output(), plane=-1)
>>> spectrum.shape == omega.shape
True

solprop.processing.fields.nearest_plane(output, z)[source]#

Return the index of the saved plane nearest a position.

Parameters:

outputsolprop.output.schema.Output: An opened output handler holding at least one plane.
zfloat: Position along the propagation, m.

Returns:

int: Index into the saved planes.

Raises:

solprop.exceptions.OutputError: If no plane has been saved.

Examples

>>> from solprop.processing._examples import example_output
>>> from solprop.processing.fields import nearest_plane
>>> nearest_plane(example_output(), 0.049)
2