solprop.grid#

Time and frequency grids for field-resolved and envelope propagation.

Two grids implement one protocol. FieldGrid carries a real field E(t) and stores its one-sided spectrum; EnvelopeGrid carries a complex envelope A(t) and stores the full spectrum in FFT order. Both carry an oversampled fine grid on which nonlinear responses are evaluated, the masks that say which bins are physical, and the Planck-taper collars that absorb what leaves the window.

The construction arithmetic, the invariants, the layout rules and the sign convention are fixed by docs/dev/conventions.md §2-§3, and the manual page docs/guide/setup/grids.md derives them. Four hazards are guarded against, each noted where it is handled: overlapping coarse spectral bins when the transform size is not a power of two, a degenerate temporal collar when the realised window equals the requested one, an envelope fine spacing that ignores the low-frequency edge, and an envelope crop that is derived from the upper edge alone and then applied symmetrically.

Functions

`planck_taper`(x, zero_low, one_low, one_high, ...)	Evaluate a Planck taper that rises then falls between four breakpoints.
`planck_taper_derivative`(x, zero_low, ...)	Evaluate the analytic derivative of `planck_taper()` with respect to `x`.
`planck_taper_derivative_error`(x, zero_low, ...)	Bound the absolute error of `planck_taper_derivative()` at each point.

Classes

`Collar`(zero_low, one_low, one_high, zero_high)	Geometry of one Planck-taper collar, as four breakpoints.
`EnvelopeGrid`(length, reference_wavelength, ...)	Grid for propagation of a complex envelope about a carrier.
`FieldGrid`(length, reference_wavelength, t, ...)	Grid for field-resolved propagation of a real electric field.
`FreeGrid`(radius_x, n_x, radius_y, n_y[, ...])	Construct a periodic Cartesian grid with x fastest in flattened columns.
`LineGrid`(radius_x, n_x[, window_factor])	Construct a periodic line grid in x for a beam uniform in y, per metre of y.
`RadialGrid`(radius, n_r[, window_factor])	Construct an axisymmetric spatial grid with a declared outer collar.
`TimeGrid`(length, reference_wavelength, t, ...)	Shared state of the field and envelope grids.

class solprop.grid.Collar(zero_low, one_low, one_high, zero_high)[source]#

Bases: object

Geometry of one Planck-taper collar, as four breakpoints.

A collar is where the propagation is absorbed rather than modelled: the spectral collar outside the requested wavelength limits, and the temporal collar outside the requested time window. solprop.boundaries turns these profiles into absorber rates; here they define the profile itself.

Attributes:

zero_lowfloat: Lower point at which the profile reaches 0.
one_lowfloat: Lower point at which the profile reaches 1.
one_highfloat: Upper point at which the profile leaves 1.
zero_highfloat: Upper point at which the profile returns to 0.

Examples

>>> import numpy as np
>>> from solprop.grid import Collar
>>> collar = Collar(0.0, 1.0, 3.0, 4.0)
>>> collar.window(np.array([0.5, 2.0])).round(3).tolist()
[0.5, 1.0]

one_high: float#

one_low: float#

property width_high: float[source]#

Return the width of the falling edge, in the collar’s own units.

Returns:

float: The distance from one_high to zero_high.

property width_low: float[source]#

Return the width of the rising edge, in the collar’s own units.

Returns:

float: The distance from zero_low to one_low.

window(x)[source]#

Evaluate the collar’s Planck taper at the given points.

Parameters:

xnumpy.ndarray: Points at which to evaluate the profile, shape (n,).

Returns:

numpy.ndarray: Window values in [0, 1], shape (n,).

window_derivative(x)[source]#

Evaluate the analytic derivative of the collar’s profile.

Parameters:

xnumpy.ndarray: Points at which to evaluate it, shape (n,).

Returns:

numpy.ndarray: d window / dx in reciprocal units of x, shape (n,); see planck_taper_derivative().

zero_high: float#

zero_low: float#

class solprop.grid.EnvelopeGrid(length, reference_wavelength, t, omega, t_fine, omega_fine, dt, dt_fine, ffac, omega_collar, t_collar)[source]#

Bases: TimeGrid

Grid for propagation of a complex envelope about a carrier.

The state is the full spectrum on omega = omega0 + v with v in FFT order, so n_omega = n_t and the time samples are complex. The kernels never shift; only processing and plotting do. Build one with from_limits(); see TimeGrid for the stored attributes.

Examples

>>> from solprop.grid import EnvelopeGrid
>>> grid = EnvelopeGrid.from_limits(1.0, 835e-9, (400e-9, 2e-6), 20e-12)
>>> grid.n_omega == grid.n_t, grid.n_to == grid.n_t
(True, True)

classmethod from_limits(length, reference_wavelength, wavelength_limits, trange, *, dt=None, thg=False, ffac=None, fft_size='pow2', collar_fraction=0.05)[source]#

Build an envelope grid from the wavelength limits and the time window.

Without third harmonic the grid is sampled at twice the widest relative frequency it must carry, so the fine and coarse grids coincide and there is no oversampling unless a finer dt is asked for. With thg=True it is sampled at six times that, and the coarse grid is cropped back to the apodised window.

The relative-frequency span is measured symmetrically, as the larger of 1.1 omega_max - omega0 (omega_max - omega0 with thg=True, where the factor 6 already leaves room for the collar) and omega0 - omega_min/2. Measuring only the upper side and applying the result to both would silently drop the whole low-frequency half of the band whenever the reference wavelength sits near the blue edge of the window.

Parameters:

lengthfloat: Propagation distance, m.
reference_wavelengthfloat: Reference wavelength, m; the carrier the envelope is written about.
wavelength_limitstuple of float: Shortest and longest wavelength of the frequency window, m, strictly increasing and positive.
trangefloat: Requested extent of the time window, s.
dtfloat or None, optional: Requested fine sample spacing, s. The value used is this or the value the frequency window demands, whichever is smaller; a finer value here turns on oversampling.
thgbool, optional: Leave room for third-harmonic generation. Sets the default ffac to 6 instead of 2 and turns on oversampling.
ffacfloat or None, optional: Fine-grid oversampling factor, at least 2. Defaults to 6 with thg and 2 without.
fft_size{‘pow2’, ‘smooth’}, optional: Round the transform sizes up to the next power of two, or to the next even 2-3-5-smooth size.
collar_fractionfloat, optional: Minimum width of each temporal collar as a fraction of the full window, strictly between 0 and 0.5. The default is 0.05.

Returns:

EnvelopeGrid: A validated grid.

Raises:

GridError: If any argument is outside its valid range, or an invariant fails.

Examples

>>> from solprop.grid import EnvelopeGrid
>>> grid = EnvelopeGrid.from_limits(1.0, 835e-9, (400e-9, 2e-6), 20e-12)
>>> grid.n_t == grid.n_to
True
>>> limits = (400e-9, 2e-6)
>>> thg = EnvelopeGrid.from_limits(1.0, 835e-9, limits, 20e-12, thg=True)
>>> thg.n_to // thg.n_t
2

kind: ClassVar[Literal['field', 'envelope']] = 'envelope'#

class solprop.grid.FieldGrid(length, reference_wavelength, t, omega, t_fine, omega_fine, dt, dt_fine, ffac, omega_collar, t_collar)[source]#

Bases: TimeGrid

Grid for field-resolved propagation of a real electric field.

The state is the one-sided spectrum on omega_k = k d_omega, k = 0 ... n_t/2, so n_omega = n_t/2 + 1 and the time samples are real. Build one with from_limits(); see TimeGrid for the stored attributes.

Examples

>>> from solprop.grid import FieldGrid
>>> grid = FieldGrid.from_limits(1.5, 800e-9, (90e-9, 4e-6), 400e-15)
>>> grid.n_t, grid.n_to, grid.n_omega
(4096, 8192, 2049)

classmethod from_limits(length, reference_wavelength, wavelength_limits, trange, *, dt=None, ffac=6.0, fft_size='pow2', collar_fraction=0.05)[source]#

Build a field grid from the wavelength limits and the requested time window.

The fine grid is sampled at ffac times the highest frequency of the window, its length is rounded up to a transform size, and the coarse grid keeps the fine bins out to just past 1.1 omega_max. Note that ffac changes the whole grid — the bin spacing, the realised window and n_omega — not only the fine one, so two values of ffac give two different grids that both satisfy the request; compare them as physical spectra on their common band, never array to array.

Parameters:

lengthfloat: Propagation distance, m.
reference_wavelengthfloat: Reference wavelength, m, for example the pump’s centre wavelength.
wavelength_limitstuple of float: Shortest and longest wavelength of the frequency window, m, strictly increasing and positive.
trangefloat: Requested extent of the time window, s. The realised window is at least this long, because the transform size is rounded up.
dtfloat or None, optional: Requested fine sample spacing, s. The value used is this or the value the frequency window demands, whichever is smaller.
ffacfloat, optional: Fine-grid oversampling factor, at least 2. The default 6 resolves a cubic response, whose third harmonic reaches 3 omega_max; 4 suffices for a response that generates no third harmonic, such as |E_a|^2 E. At exactly 2 the fine Nyquist frequency is omega_max itself, so the apodisation collar above it is truncated rather than represented.
fft_size{‘pow2’, ‘smooth’}, optional: Round the transform sizes up to the next power of two, or to the next even 2-3-5-smooth size, which tracks the requirement far more closely.
collar_fractionfloat, optional: Minimum width of each temporal collar as a fraction of the full window, strictly between 0 and 0.5. The default is 0.05.

Returns:

FieldGrid: A validated grid.

Raises:

GridError: If any argument is outside its valid range, or an invariant fails.

Examples

>>> from solprop.grid import FieldGrid
>>> grid = FieldGrid.from_limits(1.5, 800e-9, (90e-9, 4e-6), 400e-15)
>>> grid.n_t, grid.n_to
(4096, 8192)
>>> float(grid.dt / grid.dt_fine)
2.0

kind: ClassVar[Literal['field', 'envelope']] = 'field'#

class solprop.grid.FreeGrid(radius_x, n_x, radius_y, n_y, window_factor=0.1)[source]#

Bases: Displayed

Construct a periodic Cartesian grid with x fastest in flattened columns.

Attributes:

radius_x, radius_yfloat: Requested physical half-widths in m, finite and positive.
n_x, n_yint: Independent sample counts, each an integer of at least 2.
window_factorfloat: Added collar width divided by each requested half-width, in [0, 1). The periodic half-width is radius * (1 + window_factor). Zero disables the collar. The default is 0.1.
x, ynumpy.ndarray: Centred physical axes in m, shapes (n_x,) and (n_y,). Odd sample counts place zero halfway between samples.
k_x, k_ynumpy.ndarray: Angular transverse frequencies in rad/m, in unshifted FFT order.
rnumpy.ndarray: Distance from the origin in m, shape (n_nodes,), x fastest.
spatial_weights, spectral_weightsnumpy.ndarray: Area and raw-DFT Parseval weights in m^2, shapes (n_nodes,).
windownumpy.ndarray: Product of the two Planck collars, shape (n_nodes,).
on_axisnumpy.ndarray: Complex reconstruction weights at x=y=0, shape (n_nodes,).
kperp_sqnumpy.ndarray: Transverse wave-number squares in m^-2, shape (n_nodes,).

Raises:

GridError: If dimensions, the collar, or derived floating-point arrays are invalid.

Notes

Public axes are x before y. Column ix + n_x*iy represents (x[ix], y[iy]). A C-order reshape for an FFT therefore has shape (n_y, n_x). Forward transforms are unscaled; inverse transforms divide by n_x*n_y.

k_x: Float[ndarray, 'n_x']#

k_y: Float[ndarray, 'n_y']#

kperp_sq: Float[ndarray, 'n_nodes']#

property n_columns: int[source]#

Return the number of flattened transverse spectral columns.

Returns:

int: n_x*n_y.

n_x: int#

n_y: int#

on_axis: Complex[ndarray, 'n_nodes']#

r: Float[ndarray, 'n_nodes']#

radius_x: float#

radius_y: float#

spatial_weights: Float[ndarray, 'n_nodes']#

spectral_weights: Float[ndarray, 'n_nodes']#

validate()[source]#

Validate the immutable derived axes and area measures.

Raises:

GridError: If a derived array has invalid shape, dtype, finiteness or mutability.

window: Float[ndarray, 'n_nodes']#

window_factor: float = 0.1#

x: Float[ndarray, 'n_x']#

y: Float[ndarray, 'n_y']#

class solprop.grid.LineGrid(radius_x, n_x, window_factor=0.1)[source]#

Bases: FreeGrid

Construct a periodic line grid in x for a beam uniform in y, per metre of y.

The line is the x axis of a FreeGrid – the same centred samples, the same unshifted FFT frequencies and the same Planck collar – with y reduced to one sample at y = 0 carrying k_y = 0, no collar and a unit measure. It therefore runs as the "free" transform with n_y = 1, and every integral over the transverse plane becomes an integral over x per metre of y: energy in J/m and power in W/m, while intensity and fluence stay local, in W/m^2 and J/m^2 (transverse_measure). This is x-only free space; docs/guide/free-space/model.md derives how its physics differs from the plane’s.

Attributes:

radius_xfloat

Requested physical half-width in m, finite and positive.

n_xint

Number of samples in x, an integer of at least 2.

window_factorfloat

Added collar width divided by the half-width, in [0, 1). The periodic half-width is radius_x * (1 + window_factor); zero disables the collar. The default is 0.1.

radius_yfloat

math.inf: the field is uniform and unbounded in y.

n_yint

x, ynumpy.ndarray

The centred x axis in m, shape (n_x,), and y = [0.0].

k_x, k_ynumpy.ndarray

Angular frequencies in rad/m in unshifted FFT order, and k_y = [0.0].

rnumpy.ndarray

|x| in m, shape (n_x,).

spatial_weights, spectral_weightsnumpy.ndarray

Line and raw-DFT Parseval weights in m (the spacing times one metre of y per metre of y), shapes (n_x,).

windownumpy.ndarray

The x collar, shape (n_x,).

on_axisnumpy.ndarray

Complex reconstruction weights at x = 0, shape (n_x,).

kperp_sqnumpy.ndarray

k_x**2 in m^-2, shape (n_x,).

Raises:

GridError: If the half-width, the sample count or the collar is invalid.

Notes

It is a FreeGrid so that everything that runs a Cartesian grid runs it unchanged, the backends included; FreeGrid itself keeps two samples per axis. Refusals that are specific to the line, such as solprop.spatial.critical_power(), test for this class by name.

Examples

>>> from solprop.grid import LineGrid
>>> line = LineGrid(1e-3, 64)
>>> line.n_columns, line.y.tolist(), line.transverse_measure
(64, [0.0], 'per_unit_width')
>>> f"{line.spatial_weights[0]:.4e} m"
'3.4375e-05 m'

n_y: int#

radius_y: float#

property transverse_measure: Literal['per_unit_width'][source]#

Return what the grid’s transverse integrals measure.

Returns:

str: "per_unit_width": energies are J/m and powers W/m, per metre of y (solprop.modes.protocol.TransverseMeasure).

class solprop.grid.RadialGrid(radius, n_r, window_factor=0.1)[source]#

Bases: Displayed

Construct an axisymmetric spatial grid with a declared outer collar.

Attributes:

radiusfloat: Outer Dirichlet boundary in m, finite and positive.
n_rint: Number of interior radial and spectral samples, at least 2.
window_factorfloat: Fraction of the outer radius occupied by the absorbing collar, in [0, 1). Zero disables the spatial collar; the default 0.1 uses its outer tenth.
hankelsolprop.hankel.QDHT: The raw transform, its physical scalings and integration weights.
windownumpy.ndarray: Planck window on the radial samples, shape (n_r,), dimensionless.

Raises:

GridError: If a dimension or the collar fraction is invalid.

Examples

>>> from solprop.grid import RadialGrid
>>> grid = RadialGrid(1e-3, 64)
>>> grid.r.shape, grid.k_r.shape, bool(grid.r[0] > 0)
((64,), (64,), True)

hankel: QDHT#

property k_r: Float[ndarray, 'n_r'][source]#

Return the transverse spectral samples.

Returns:

numpy.ndarray: Positive angular frequencies in rad/m, shape (n_r,).

property n_columns: int[source]#

Return the number of transverse spectral columns.

Returns:

int: The radial sample count.

n_r: int#

property r: Float[ndarray, 'n_r'][source]#

Return the positive interior radial samples.

Returns:

numpy.ndarray: Radii in m, shape (n_r,); zero and the outer boundary are excluded.

property radial_weights: Float[ndarray, 'n_r'][source]#

Return the quadrature weights of the transverse area integral.

Returns:

numpy.ndarray: Weights of 2*pi*r*dr, in m^2, shape (n_r,).

radius: float#

property spectral_weights: Float[ndarray, 'n_r'][source]#

Return the quadrature weights of the transverse spectral integral.

Returns:

numpy.ndarray: Weights of 2*pi*k_r*dk_r, in m^-2, shape (n_r,).

validate()[source]#

Check the radial axes, weights and spatial window.

Raises:

GridError: If an axis is unordered or an array is mutable, nonfinite or inconsistent.

window: Float[ndarray, 'n_r']#

window_factor: float = 0.1#

class solprop.grid.TimeGrid(length, reference_wavelength, t, omega, t_fine, omega_fine, dt, dt_fine, ffac, omega_collar, t_collar)[source]#

Bases: Displayed

Shared state of the field and envelope grids.

The eleven stored attributes below are the whole of a grid’s state; every other quantity named in docs/dev/kernel-contract.md §2 is a property derived from them, so each has exactly one definition. Instances are frozen and hold private read-only copies of their axes, so a grid can be shared without defensive copying. Equality is by identity, because comparing whole axes elementwise is never what a caller means.

Build one with FieldGrid.from_limits() or EnvelopeGrid.from_limits() rather than by hand; the constructor is the low-level path and only validate() stands behind it.

Attributes:

lengthfloat: Propagation distance the grid was built for, m.
reference_wavelengthfloat: Reference wavelength, m; fixes omega0 and the co-moving frame.
tnumpy.ndarray: Coarse time axis, s, shape (n_t,), centred on zero.
omeganumpy.ndarray: Coarse angular frequency axis, rad/s, shape (n_omega,). Non-negative and increasing on a field grid; in FFT order about omega0 on an envelope grid.
t_finenumpy.ndarray: Fine (oversampled) time axis, s, shape (n_to,), centred on zero.
omega_finenumpy.ndarray: Fine angular frequency axis, rad/s, shape (n_omega_fine,).
dtfloat: Coarse sample spacing, s.
dt_finefloat: Fine sample spacing, s.
ffacfloat: Oversampling factor: the fine grid is sampled at ffac times the highest frequency of the window. At least 2.
omega_collarCollar: Spectral collar geometry, rad/s: rises over [omega_min/2, omega_min] and falls over [omega_max, 1.1 omega_max].
t_collarCollar: Temporal collar geometry, s, symmetric about zero and guaranteed non-degenerate.
kind{‘field’, ‘envelope’}: Spectral layout, fixed by the subclass and passed to solprop.fourier.

Examples

>>> from solprop.grid import FieldGrid, TimeGrid
>>> isinstance(FieldGrid.from_limits(1.5, 800e-9, (90e-9, 4e-6), 400e-15), TimeGrid)
True

property active_mask: Bool[ndarray, 'n_omega'][source]#

Return the bins where the physics lives, inside the apodised window.

The dispersion relation, the noise and the error norm are evaluated here and nowhere else.

Returns:

numpy.ndarray: Boolean mask over the coarse frequency axis, shape (n_omega,).

property crop_index: ndarray[source]#

Return the positions of the coarse spectral bins inside the fine spectrum.

Returns:

numpy.ndarray: Index array of numpy.intp, shape (n_omega,).

property crop_scale: float[source]#

Return the factor applied when cropping the fine spectrum back.

Returns:

float: n_t / n_to.

property d_omega: float[source]#

Return the spectral bin spacing, rad/s.

Returns:

float: 2 pi / (n_to dt_fine), the same on the coarse and the fine grid.

dt: float#

dt_fine: float#

ffac: float#

property inband_mask: Bool[ndarray, 'n_omega'][source]#

Return the bins inside the requested wavelength limits.

Returns:

numpy.ndarray: Boolean mask over the coarse frequency axis, shape (n_omega,), true where omega_min <= omega <= omega_max.

kind: ClassVar[Literal['field', 'envelope']] = 'field'#

length: float#

property n_omega: int[source]#

Return the number of coarse spectral bins.

Returns:

int: The length of omega.

property n_omega_fine: int[source]#

Return the number of fine spectral bins.

Returns:

int: The length of omega_fine.

property n_t: int[source]#

Return the number of coarse time samples.

Returns:

int: The length of t.

property n_to: int[source]#

Return the number of fine time samples.

Returns:

int: The length of t_fine.

omega: Float[ndarray, 'n_omega']#

property omega0: float[source]#

Return the reference angular frequency, rad/s.

Returns:

float: 2 pi c / reference_wavelength.

omega_collar: Collar#

omega_fine: Float[ndarray, 'n_omega_fine']#

property omega_max: float[source]#

Return the highest in-band angular frequency, rad/s.

Returns:

float: The frequency of the shortest requested wavelength.

property omega_max_window: float[source]#

Return the upper edge of the apodised frequency window, rad/s.

Returns:

float: 1.1 omega_max, where the spectral collar reaches zero.

property omega_min: float[source]#

Return the lowest in-band angular frequency, rad/s.

Returns:

float: The frequency of the longest requested wavelength.

property omega_window: Float[ndarray, 'n_omega'][source]#

Return the spectral collar profile on the coarse frequency axis.

Returns:

numpy.ndarray: Window values in [0, 1], shape (n_omega,), equal to 1 in band.

property pad_scale: float[source]#

Return the factor applied when padding the coarse spectrum out.

Returns:

float: n_to / n_t.

reference_wavelength: float#

t: Float[ndarray, 'n_t']#

t_collar: Collar#

t_fine: Float[ndarray, 'n_to']#

property t_window: Float[ndarray, 'n_t'][source]#

Return the temporal collar profile on the coarse time axis.

Returns:

numpy.ndarray: Window values in [0, 1], shape (n_t,), equal to 1 inside the window.

property t_window_fine: Float[ndarray, 'n_to'][source]#

Return the temporal collar profile on the fine time axis.

Both temporal windows come from the one collar geometry, so the coarse and fine apodisations agree exactly wherever their samples coincide, rather than each window being built from its own axis’s extremes, which differ at the upper edge.

Returns:

numpy.ndarray: Window values in [0, 1], shape (n_to,), 1 inside the window.

validate()[source]#

Check every grid invariant, raising on the first that fails.

The invariants are those of docs/dev/conventions.md §3. Both from_limits classmethods call this, so a grid they return is always valid; call it again after building or loading one by hand.

Raises:

GridError: If any invariant fails, naming the offending quantity.

Examples

>>> from solprop.grid import FieldGrid
>>> FieldGrid.from_limits(1.5, 800e-9, (90e-9, 4e-6), 400e-15).validate()

solprop.grid.planck_taper(x, zero_low, one_low, one_high, zero_high)[source]#

Evaluate a Planck taper that rises then falls between four breakpoints.

The window is 0 at or below zero_low, rises smoothly to 1 at one_low, stays 1 up to one_high, and falls back to 0 at zero_high. Every derivative is continuous, which is what distinguishes a Planck taper from the Tukey window that scipy.signal.windows provides, and why it is written out here.

Parameters:

xnumpy.ndarray: Points at which to evaluate the window, shape (n,), in the same units as the breakpoints: rad/s for a spectral collar, s for a temporal one.
zero_lowfloat: Lower point at which the window reaches 0.
one_lowfloat: Lower point at which the window reaches 1. Must exceed zero_low.
one_highfloat: Upper point at which the window leaves 1. Must be at least one_low.
zero_highfloat: Upper point at which the window returns to 0. Must exceed one_high.

Returns:

numpy.ndarray: Window values in [0, 1], shape (n,).

Raises:

GridError: If the four breakpoints are not ordered, so that a taper edge would have zero width and the window would step rather than taper.

References

D. J. A. McKechan, C. Robinson and B. S. Sathyaprakash, “A tapering window for time-domain templates and simulated signals in the detection of gravitational waves from coalescing compact binaries”, Class. Quantum Grav. 27, 084020 (2010), eq. (7).

Examples

>>> import numpy as np
>>> from solprop.grid import planck_taper
>>> window = planck_taper(np.array([0.0, 1.0, 2.0, 3.0, 4.0]), 0.0, 1.0, 3.0, 4.0)
>>> window.tolist()
[0.0, 1.0, 1.0, 1.0, 0.0]

solprop.grid.planck_taper_derivative(x, zero_low, one_low, one_high, zero_high)[source]#

Evaluate the analytic derivative of planck_taper() with respect to x.

On an edge with outer breakpoint a and inner b, width = b - a (negative on the falling edge), s = (x - a) / width and r = (b - x) / width, both in (0, 1). The gate is g = 1 / (1 + exp(zeta)) with zeta = 1/s - 1/r, and

\[g'(x) = \frac{h(s)}{\mathrm{width}},\qquad h = \frac{1}{4\cosh^2(\zeta/2)}\left(\frac{1}{s^2}+\frac{1}{r^2}\right) > 0,\]

so g' is positive on the rising edge, negative on the falling edge, integrates to +1 and -1 over them, and is zero outside the tapers. The factor 1/(4 cosh^2(zeta/2)) is evaluated as e/(1 + e)^2 with e = exp(-|zeta|), which never overflows. s and r are formed separately, so neither is obtained by cancellation. Where min(s, r) <= 1/740 the asymptotic zero is returned and no reciprocal is formed: there the analytic value is below 1e-300 / |width|. Points at or beyond the breakpoints return exactly zero. The evaluation never forms an infinite intermediate (SEP 26 section 4.2).

Parameters:

xnumpy.ndarray: Points at which to evaluate the derivative, shape (n,), in the units of the breakpoints.
zero_low, one_low, one_high, zero_highfloat: Breakpoints of planck_taper(), in the same units.

Returns:

numpy.ndarray: dg/dx in reciprocal units of x, shape (n,).

Raises:

GridError: If the breakpoints are not ordered, or an edge is so narrow that its peak derivative 2 / width would overflow.