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.

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.

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_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.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.

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.

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_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