solprop.modes.quadrature#

Fixed transverse quadrature rules and the mode matrices built on them.

The multimode transform integrates over the cross-section with one fixed rule chosen once per problem, never an adaptive one (ADR-010). That is not a compromise: the Kerr polarisation is cubic in the field, so its projection is a quartic in the transverse profiles, and for Marcatili modes those are Bessel functions of \(r/a\) times \(\cos h\theta\) – analytic, and integrated to rounding by a Gauss rule of modest order. An adaptive rule buys nothing and costs a data-dependent time per step and an error estimate that can be fooled.

The rules:

Radial, \(n_r\) Gauss-Legendre nodes on \(\xi = r/a \in (0, 1)\), with the Jacobian and the azimuthal \(2\pi\) folded into the weights, \(w_i = 2\pi a^2 \xi_i w^{\rm GL}_i / 2\). With radial_rule="gauss_kronrod" the count is rounded up to \(2n+1\) and the rule becomes the Kronrod extension of the \(n\)-point Gauss rule, whose Gauss subset is an embedded coarse rule for the error estimate. Nodes are strictly interior, so the wall is never sampled.

Field scale (SEP 19), for one Marcatili HE_11 mode in a transversely uniform medium, where every node sees the axial waveform scaled by \(\eta = J_0(u_{01}\xi)\): "field_scale" is the Gauss rule in \(\eta\) for the HE_11 projection measure and "field_scale_radau" its Radau variant with a node on the axis (field_scale_rule()). They exist for responses so expensive that the number of atoms per call matters, carry area weights, which are finite at the axis and integrate the source rather than the area, and have no embedded coarse rule.

Azimuthal, \(n_\theta\) equispaced nodes offset by half a step, \(\theta_j = (j - \tfrac12) 2\pi/n_\theta\): the periodic trapezoid rule, which is exact for trigonometric polynomials of degree at most \(n_\theta - 1\). The offset matters. On-grid nodes at \(0, \pi/2, \pi, 3\pi/2\) annihilate products such as \(E_xE_y\) for a superposition of hybrid modes, which can make an adaptive cubature such as pcubature report convergence on a wrong value; the offset rule cannot be fooled that way and there is no convergence decision in it to fool. A quartic product of modes with azimuthal order up to \(h\) reaches degree \(4h\), so transverse_quadrature() warns below \(n_\theta = 4h + 1\).

Cartesian, Gauss-Legendre in both directions, for rectangular waveguides. A product of trigonometric eigenfunctions has no finite exactness degree there, so the bound is the measured cartesian_nodes(), \(6h + 6\) per axis, and transverse_quadrature() warns below it.

Slab, Gauss-Legendre in x alone on \([-a, a]\), for the hollow planar guide uniform in y (solprop.modes.slab). The rule emits (x, 0) points with the weights \(a\,w_i\) times a unit y weight, so everything downstream of the rule is the Cartesian pipeline with one y node: an integral over it is a per-unit-width one. The same cartesian_nodes() bound holds – re-measured on the one-dimensional quartic products – and with radial_rule="gauss_kronrod" the x rule is the Kronrod extension, so the slab has the quadrature alarm a Cartesian product rule lacks.

The Kronrod extension is computed, not tabulated, by Laurie’s algorithm [1]: build the Jacobi-Kronrod matrix from the Legendre recurrence coefficients and take its eigenvalues and first eigenvector components. That gives any order rather than the two that scipy.integrate happens to ship, and it is checked node by node against QuadGK’s values in tests/data/luna_modes.json.

Tapers. For a scale-invariant mode the profile depends on r only through r / a and the normalisation goes as \(a^2\), so at nodes placed at \(r = a\xi\) the normalised field is exactly proportional to \(1/a\) and the weights to \(a^2\). Moving along a taper therefore rescales the synthesis matrix by \(a(0)/a(z)\) and the projection matrix by \(a(z)/a(0)\), exactly, instead of re-evaluating thousands of Bessel functions: ModeMatrices.rescaled().

The dual pair. The propagated basis is \(e_m/\sqrt{N_m}\) with \(N_m = (\varepsilon_0 c/2)\int|e_m|^2\,dA\) (conventions.md section 4), so its overlaps are \(2/(\varepsilon_0 c)\) rather than one, and the coefficient of a field in that basis is \((\varepsilon_0 c/2)\int P\,e_m/\sqrt{N_m}\,dA\). The projection matrix therefore carries BASIS_NORMALISATION beside the quadrature weights, so that projection @ synthesis is the identity on an orthonormal set and a backend applies the pair as it is.

The manual page is docs/guide/waveguides/mode-sets.md.

References#

[1]

D. P. Laurie, “Calculation of Gauss-Kronrod quadrature rules”, Math. Comp. 66, 1133-1145 (1997), doi:10.1090/S0025-5718-97-00861-2.

Module Attributes

`BASIS_NORMALISATION`	The factor that makes the projection the dual of the synthesis, in F/(m s).
`DEFAULT_N_THETA`	Default azimuthal node count, exact for cubic products up to azimuthal order 3.
`RadialRule`	The first-coordinate rule families (SEP 19).
`FIELD_SCALE_RULES`	The rule families that are Gaussian in the HE_11 field scale.

Functions

`cartesian_nodes`(order)	Return the nodes per axis a Cartesian quartic product needs.
`check_field_scale_modes`(modes)	Refuse a field-scale rule for anything but one Marcatili HE_11 profile.
`field_scale_rule`(n, *[, radau])	Return a Gauss or Gauss-Radau rule in the HE_11 field scale.
`gauss_kronrod`(n)	Return the `2 n + 1` point Kronrod extension of the n-point Gauss rule.
`mode_matrix`(modes, quadrature[, z, components])	Build the synthesis and projection matrices of a mode set on a rule.
`on_axis_rule`(kind)	Return a one-node rule at the transverse origin of a waveguide.
`quadrature_rule`(kind, *, full, n_r[, ...])	Build a transverse rule from its geometry alone.
`taper_scale`(mode, z[, reference])	Return the factor a scale-invariant mode's matrices rescale by.
`transverse_quadrature`(modes, *, n_r[, ...])	Build the transverse rule a mode set needs.

Classes

`ModeMatrices`(synthesis, projection, ...)	Synthesis and projection matrices in the layout the kernel contract expects.
`TransverseQuadrature`(kind, full, n_r, ...[, ...])	A fixed product rule over the cross-section of a waveguide.

solprop.modes.quadrature.BASIS_NORMALISATION: Final = 0.0013272093648958553#

The factor that makes the projection the dual of the synthesis, in F/(m s).

The propagated basis e_m / sqrt(N_m) has overlaps 2 / (eps0 c) and not one, so the coefficient of a field in it is (eps0 c / 2) int P e_m/sqrt(N_m) dA: mode_matrix() folds this into both projections beside the quadrature weights. Without it a multimode run is wrong by the constant 753.5.

solprop.modes.quadrature.DEFAULT_N_THETA: Final = 16#: Default azimuthal node count, exact for cubic products up to azimuthal order 3.

solprop.modes.quadrature.FIELD_SCALE_RULES: Final = ('field_scale', 'field_scale_radau')#: The rule families that are Gaussian in the HE_11 field scale.

class solprop.modes.quadrature.ModeMatrices(synthesis, projection, projection_coarse, n_nodes, n_pol)[source]#

Bases: object

Synthesis and projection matrices in the layout the kernel contract expects.

Synthesis takes modal amplitudes to the field at the quadrature nodes; projection takes a polarisation at the nodes back to modal amplitudes, with the quadrature weights folded in. Column i * n_pol + p of the projection – row of the synthesis – is polarisation component p at node i, so that a (n_to, n_pol, n_nodes) view of the real-space arrays is free and each node’s column is contiguous (docs/dev/kernel-contract.md section 4).

Attributes:

synthesisnumpy.ndarray: Shape (n_nodes * n_pol, n_modes): the normalised transverse fields.
projectionnumpy.ndarray: Shape (n_modes, n_nodes * n_pol): the dual basis – the same fields with the quadrature weights and BASIS_NORMALISATION folded in, so that projection @ synthesis is the identity on an orthonormal set.
projection_coarsenumpy.ndarray: The projection under the embedded coarse rule; equal to projection when the rule has no embedded pair, in which case the error estimate is identically zero and must not be read as convergence.
n_nodesint: Quadrature nodes.
n_polint: Polarisation components carried, 1 or 2.

n_nodes: int#

n_pol: int#

projection: Float[ndarray, 'n_modes n_columns']#

projection_coarse: Float[ndarray, 'n_modes_coarse n_columns_coarse']#

rescaled(scale)[source]#

Return the matrices for a scale-invariant mode set at a new core radius.

Parameters:

scalefloat: The ratio a(z_reference) / a(z), > 0.

Returns:

ModeMatrices: Synthesis multiplied by scale and projection divided by it, which is exact: the normalised field goes as 1 / a at nodes placed at r = a xi, and the weights as a^2.

Raises:

ModeError: If scale is not positive.

Examples

>>> from solprop.modes import MarcatiliMode
>>> from solprop.modes.quadrature import mode_matrix, transverse_quadrature
>>> modes = [MarcatiliMode(100e-6, "He", 1.0)]
>>> rule = transverse_quadrature(modes, n_r=8)
>>> matrices = mode_matrix(modes, rule)
>>> matrices.rescaled(2.0).synthesis.shape
(16, 1)

synthesis: Float[ndarray, 'n_rows n_modes']#

solprop.modes.quadrature.RadialRule#

The first-coordinate rule families (SEP 19).

"gauss_kronrod" and "gauss" are Gauss-Legendre on xi = r / a, the former with the Kronrod extension and its embedded coarse rule. "field_scale" and "field_scale_radau" are Gauss and Gauss-Radau rules in the HE_11 field scale eta = J_0(u_01 xi) (field_scale_rule()), valid for one Marcatili HE_11 mode in a transversely uniform medium only.

alias of Literal[‘gauss_kronrod’, ‘gauss’, ‘field_scale’, ‘field_scale_radau’]

class solprop.modes.quadrature.TransverseQuadrature(kind, full, n_r, n_theta, kronrod, xi_r, w_r, w_r_coarse, xi_theta, w_theta, w_theta_coarse, breaks=(), radial_rule='gauss', area_r=None, area_r_coarse=None, requested_n_r=None)[source]#

Bases: object

A fixed product rule over the cross-section of a waveguide.

Node p = i + j * n_r is coordinate-1 node i and coordinate-2 node j, so the radial index varies fastest.

The nodes and weights here are on the reference cross-section; nodes() and weights() map them onto the geometry a mode declares, folding the Jacobian into the weights.

Attributes:

kindstr: "polar", "cartesian" or "slab".
fullbool: Whether the second coordinate is resolved. A polar rule with full=False is radial only, valid when every product of the mode fields is azimuthally symmetric – any set of HE_1m modes. A slab rule is never full: its one y node sits at y = 0 with unit weight.
n_rint: Nodes in the first coordinate, after any rounding up for a Kronrod rule. With breaks it is the TOTAL count across the panels, each of which carries the requested count, so a rule with one break holds twice what was asked for.
n_thetaint: Nodes in the second coordinate; 1 for a radial-only rule.
kronrodbool: Whether the first coordinate carries an embedded coarse rule.
xi_rnumpy.ndarray: First-coordinate nodes on the reference interval.
w_r, w_r_coarsenumpy.ndarray or None: First-coordinate weights for d xi (polar) or dx on the reference interval; None when the rule carries area weights instead.
xi_theta, w_theta, w_theta_coarsenumpy.ndarray: Second-coordinate nodes and weights.
breakstuple of float: Interior radial panel boundaries on the unit interval, in ascending order.
radial_rulestr: The first-coordinate family, one of RadialRule.
area_r, area_r_coarsenumpy.ndarray or None: Dimensionless radial area weights q of int f xi d xi, for a field-scale rule (SEP 19 section 4.2), whose Radau member has a node on the axis where the d xi form would divide by zero. Exactly one of w_r and area_r is present. They are source-adapted: they integrate the linear and Kerr source exactly, not a constant area density, and must not be renormalised to the geometrical area.
requested_n_rint or None: The first-coordinate count asked of quadrature_rule(), before any Kronrod rounding or panel multiplication; None for a rule built directly.

Raises:

ModeError: If both or neither weight representation is present, or area weights are given for anything but a radial-only polar rule without panels.

area_r: Float[ndarray, 'n_first_area'] | None = None#

area_r_coarse: Float[ndarray, 'n_first_area_coarse'] | None = None#

breaks: tuple[float, ...] = ()#

full: bool#

property has_embedded_rule: bool[source]#

Return whether the coarse weights differ from the fine ones.

Returns:

bool: True when there is a genuine embedded rule, so that the difference of the two projections is an error estimate rather than identically zero.

kind: Literal['polar', 'cartesian', 'slab']#

kronrod: bool#

property n_nodes: int[source]#

Return the number of quadrature nodes.

Returns:

int: n_r * n_theta.

n_r: int#

n_theta: int#

nodes(dimlimits)[source]#

Return the quadrature nodes on a mode’s cross-section.

Parameters:

dimlimitstuple: As returned by solprop.modes.protocol.Mode.dimlimits().

Returns:

numpy.ndarray: Shape (n_nodes, 2): (r, theta) in m and rad, or (x, y) in m; a slab rule’s y is zero at every node.

Raises:

ModeError: If the cross-section kind does not match the rule.

Examples

>>> import math
>>> from solprop.modes.quadrature import quadrature_rule
>>> rule = quadrature_rule("polar", full=False, n_r=4)
>>> rule.nodes(("polar", (0.0, 0.0), (100e-6, 2 * math.pi))).shape
(4, 2)

radial_rule: Literal['gauss_kronrod', 'gauss', 'field_scale', 'field_scale_radau'] = 'gauss'#

record(dimlimits)[source]#

Return the rule’s identity for the run metadata and its resume key.

Parameters:

dimlimitstuple: As returned by solprop.modes.protocol.Mode.dimlimits().

Returns:

dict: The family, the requested and realised counts, the total node count, whether there is an embedded rule, the panels, whether this is the one-node effective-area rule of SEP 19, and the nodes and physical weights on that cross-section (m and m^2, or m for a slab).

Examples

>>> import math
>>> from solprop.modes.quadrature import quadrature_rule
>>> rule = quadrature_rule(
...     "polar", full=False, n_r=32, radial_rule="gauss_kronrod"
... )
>>> record = rule.record(("polar", (0.0, 0.0), (100e-6, 2 * math.pi)))
>>> record["radial_rule"], record["requested_n_r"], record["n_nodes"]
('gauss_kronrod', 32, 33)

requested_n_r: int | None = None#

w_r: Float[ndarray, 'n_first_weights'] | None#

w_r_coarse: Float[ndarray, 'n_first_coarse'] | None#

w_theta: Float[ndarray, 'n_second_weights']#

w_theta_coarse: Float[ndarray, 'n_second_coarse']#

weights(dimlimits, *, coarse=False)[source]#

Return the quadrature weights, with the area element folded in.

Parameters:

dimlimitstuple: As returned by solprop.modes.protocol.Mode.dimlimits().
coarsebool, optional: Return the embedded coarse rule’s weights instead of the fine ones.

Returns:

numpy.ndarray: Shape (n_nodes,). Polar: a^2 xi_i w_i w_j, the Jacobian r of r dr dtheta included, or a^2 q_i w_j for a rule carrying radial area weights q, which already include it. Cartesian: the product of the half-widths and the two weights. Slab: the half-gap times the x weight, times the unit y weight, in m – a per-unit-width measure.

Raises:

ModeError: If the cross-section kind does not match the rule.

Examples

>>> import math
>>> from solprop.modes.quadrature import quadrature_rule
>>> rule = quadrature_rule("polar", full=False, n_r=8)
>>> limits = ("polar", (0.0, 0.0), (100e-6, 2 * math.pi))
>>> area = float(rule.weights(limits).sum())
>>> bool(abs(area / (math.pi * 100e-6**2) - 1.0) < 1e-15)
True

xi_r: Float[ndarray, 'n_first']#

xi_theta: Float[ndarray, 'n_second']#

solprop.modes.quadrature.cartesian_nodes(order)[source]#

Return the nodes per axis a Cartesian quartic product needs.

A product of trigonometric eigenfunctions has no finite polynomial exactness degree, so this is a measured bound rather than a derived one: 6 * order + 6 covers the counts that reach a relative 1e-12 on the quartic field product of a rectangular mode – 12 nodes at order 1, 18 at 2, 22 at 3, 26 at 4 and 35 at 6. Higher orders need their own resolution check.

The one-dimensional quartic products of a slab set (solprop.modes.slab) were measured separately against 30-digit mpmath integrals: every product of orders up to h reaches 1e-12 of the int V^4 scale at 12, 17, 22, 26, 30 and 34 Gauss nodes for h = 1 to 6, and stays there at every count up to 6 h + 6, so the same bound serves the slab’s x rule.

Parameters:

orderint: Highest transverse order of the set, max(n, m), >= 0.

Returns:

int: Nodes per axis for quartic products of that order.

Examples

>>> from solprop.modes.quadrature import cartesian_nodes
>>> cartesian_nodes(1), cartesian_nodes(3)
(12, 24)

solprop.modes.quadrature.check_field_scale_modes(modes)[source]#

Refuse a field-scale rule for anything but one Marcatili HE_11 profile.

Parameters:

modessequence of Mode: The mode set the rule would integrate.

Raises:

ModeError: Naming the mode set, if it has more than one mode, or a mode whose transverse profile is not Marcatili’s J_0(u_01 r / a): a mode merely labelled HE_11 with another radial profile does not qualify.

solprop.modes.quadrature.field_scale_rule(n, *, radau=False)[source]#

Return a Gauss or Gauss-Radau rule in the HE_11 field scale.

For one HE_11 mode in a transversely uniform medium every node sees the same waveform scaled by \(\eta = J_0(u_{01}\xi)\), so the radial projection is \(\int_0^1 F[\eta A]\,d\mu(\eta)\) with \(d\mu = J_0(u_{01}\xi)\,\xi\,d\xi\) (SEP 19 section 4). The Gauss rule for \(\mu\) is computed from its recurrence coefficients, by the Lanczos process on a refined Gauss-Legendre discretisation in \(\xi\) [1], and the Golub-Welsch eigenvalue method [2]; it is exact for polynomials in \(\eta\) of degree \(2n - 1\), which from two nodes includes the linear and Kerr terms. The Radau rule fixes one node at \(\eta = 1\), the axis, by replacing the last diagonal entry with \(1 - \beta_{n-1} e^T(I - J_{n-1})^{-1}e\) [3] and is exact to degree \(2n - 2\). One Gauss node is the effective-area special rule, not the one-point Gauss rule: \(\eta_e^2 = M_3/M_1\) and \(\lambda_e = M_3/\eta_e^3\), with \(M_k = \int\eta^k d\mu\), which matches moments 1 and 3.

Every rule is checked before it is returned: finite positive weights, distinct nodes in \((0, 1]\), and every moment it should integrate against a measure refined once more, to a relative 1e-12.

Parameters:

nint: Nodes, >= 1 for the Gauss family and >= 2 for Radau.
radaubool, optional: Fix one node on the axis.

Returns:

tuple of numpy.ndarray: eta in decreasing order (increasing radius) and the weights lambda of \(\mu\), read-only.

Raises:

ModeError: If the order is too small or the rule fails its own check.

References

[1]

W. Gautschi, Orthogonal Polynomials: Computation and Approximation (Oxford University Press, 2004), doi:10.1093/oso/9780198506720.001.0001.

[2]

G. H. Golub and J. H. Welsch, Math. Comp. 23, 221 (1969), doi:10.1090/S0025-5718-69-99647-1.

[3]

1. Golub, SIAM Rev. 15, 318 (1973), doi:10.1137/1015032.

Examples

>>> from solprop.modes.quadrature import field_scale_rule
>>> eta, weights = field_scale_rule(5)
>>> eta.size, bool(eta[0] < 1.0)
(5, True)
>>> radau, _ = field_scale_rule(5, radau=True)
>>> float(radau[0])
1.0

solprop.modes.quadrature.gauss_kronrod(n)[source]#

Return the 2 n + 1 point Kronrod extension of the n-point Gauss rule.

Parameters:

nint: Order of the underlying Gauss rule, >= 1.

Returns:

tuple of numpy.ndarray: Nodes on [-1, 1] in increasing order, their weights, and the embedded coarse weights: the Gauss weights at the Gauss nodes and zero elsewhere, so that the two rules share every node they have in common and their difference is an error estimate.

Raises:

ModeError: If n is less than 1.

Examples

>>> from solprop.modes.quadrature import gauss_kronrod
>>> nodes, weights, coarse = gauss_kronrod(7)
>>> nodes.size, int((coarse > 0.0).sum())
(15, 7)
>>> bool(abs(weights.sum() - 2.0) < 16 * np.spacing(2.0))
True

solprop.modes.quadrature.mode_matrix(modes, quadrature, z=0.0, *, components=('x', 'y'))[source]#

Build the synthesis and projection matrices of a mode set on a rule.

Takes the rule rather than bare nodes: the projection needs the weights, and passing the rule keeps the two from ever being paired wrongly.

Every mode is evaluated on the cross-section of the first, which is what makes one product rule serve the whole set; a set whose limits disagree is rejected.

Parameters:

modessequence of Mode: The mode set, non-empty.
quadratureTransverseQuadrature: The rule, whose cross-section kind must match the modes’.
zfloat, optional: Position along the waveguide, m.
componentssequence of str, optional: Which transverse components to carry: ("x", "y") for a vector run, or one of them for a scalar run.

Returns:

ModeMatrices: The matrices.

Raises:

ModeError: If modes is empty, the components are not a subset of ("x", "y"), the modes do not share a cross-section, or a field-scale rule is paired with anything but one Marcatili HE_11 mode – checked here as well as in transverse_quadrature(), so that a rule built by quadrature_rule() cannot reach a translation with the wrong profile.

Examples

>>> from solprop.modes import MarcatiliMode
>>> from solprop.modes.quadrature import mode_matrix, transverse_quadrature
>>> modes = [MarcatiliMode(125e-6, "He", 1.0, m=m) for m in (1, 2, 3)]
>>> rule = transverse_quadrature(modes, n_r=24)
>>> matrices = mode_matrix(modes, rule, components=("x",))
>>> matrices.synthesis.shape, matrices.projection.shape
((24, 3), (3, 24))

solprop.modes.quadrature.on_axis_rule(kind)[source]#

Return a one-node rule at the transverse origin of a waveguide.

The axis is not a quadrature node – a Gauss rule on (0, 1) never places one there – so the on-axis intensity needs the mode fields evaluated at r = 0. The rule is the vehicle for that: solprop.modes.quadrature.mode_matrix() is what knows how to normalise a mode’s profile, and giving it a rule keeps that knowledge in the physics layer where it belongs. Its weights are never used.

Parameters:

kindstr: The polar, Cartesian or slab cross-section the mode set declares; a slab’s origin is its mid-plane x = 0.

Returns:

solprop.modes.quadrature.TransverseQuadrature: A rule with one node at the origin.

Raises:

ModeError: If the cross-section kind is unknown.

solprop.modes.quadrature.quadrature_rule(kind, *, full, n_r, n_theta=16, radial_rule='gauss', breaks=())[source]#

Build a transverse rule from its geometry alone.

transverse_quadrature() is the entry point for a mode set; this is the one for a bare cross-section, which the tests and the rectangular waveguides use.

Parameters:

kindstr: "polar", "cartesian" or "slab".
fullbool: Whether to resolve the second coordinate. A Cartesian rule must; a slab rule must not, having no y to resolve.
n_rint: Nodes in the first coordinate: >= 2, rounded up to the smallest odd count (at least 3) for "gauss_kronrod"; >= 1 for "field_scale", whose one-node member is the effective-area rule; >= 4 for "field_scale_radau".
n_thetaint, optional: Nodes in the second coordinate, >= 2 when full; ignored by a slab rule.
radial_rulestr, optional: The first-coordinate family, one of RadialRule. "gauss_kronrod" uses the Kronrod extension, giving an embedded coarse rule for the transform’s error estimate; a Cartesian product rule has none and reads it as "gauss". The field-scale families need a radial-only polar rule without breaks; which modes they may integrate is checked where the modes are known (check_field_scale_modes()).
breakstuple of float, optional: Strictly increasing radial panel boundaries inside (0, 1), for polar rules. Each panel receives n_r nodes before any Kronrod rounding.

Returns:

TransverseQuadrature: The rule.

Raises:

ModeError: If the kind or the rule family is unknown, a Cartesian rule is asked not to be full, a slab rule is asked to be full, a node count is not an integer or is too small, or a field-scale rule is asked of anything but a radial-only polar cross-section without panels.

Examples

>>> from solprop.modes.quadrature import quadrature_rule
>>> rule = quadrature_rule("polar", full=True, n_r=16, n_theta=8)
>>> rule.n_nodes
128
>>> slab = quadrature_rule("slab", full=False, n_r=12)
>>> limits = ("slab", (-50e-6,), (50e-6,))
>>> slab.n_nodes, bool(slab.nodes(limits)[:, 1].any())
(12, False)
>>> f"{slab.weights(limits).sum():.1e}"
'1.0e-04'

solprop.modes.quadrature.taper_scale(mode, z, reference=0.0)[source]#

Return the factor a scale-invariant mode’s matrices rescale by.

Parameters:

modeMode: The mode, which must report scale_invariant.
zfloat: Position along the waveguide, m.
referencefloat, optional: Position the matrices were built at, m.

Returns:

float: a(reference) / a(z), to be passed to ModeMatrices.rescaled(). For a slab it is the square root of that: the normalised field goes as 1/sqrt(a) there, because N goes as a and not a^2, and the weights as a.

Raises:

ModeError: If the mode is not scale invariant.

Examples

>>> from solprop.modes import MarcatiliMode
>>> from solprop.modes.quadrature import taper_scale
>>> from solprop.profiles import Linear
>>> mode = MarcatiliMode(Linear(100e-6, 50e-6, 1.0), "He", 1.0)
>>> f"{taper_scale(mode, 1.0):.1f}"
'2.0'

solprop.modes.quadrature.transverse_quadrature(modes, *, n_r, n_theta=16, radial_rule='gauss', full=None, z=0.0)[source]#

Build the transverse rule a mode set needs.

The geometry comes from the first mode’s cross-section, and whether the azimuthal coordinate has to be resolved from the modes’ azimuthal orders: a set of HE_1m modes has azimuthally symmetric products and needs the radial rule only.

Parameters:

modessequence of Mode: The mode set the rule will integrate, non-empty.
n_rint: Nodes in the first coordinate; the minimum depends on the family (quadrature_rule()).
n_thetaint, optional: Nodes in the azimuthal coordinate, >= 2. Ignored by a radial-only rule.
radial_rulestr, optional: The first-coordinate family, one of RadialRule: "gauss_kronrod" gives an embedded coarse rule; "field_scale" and "field_scale_radau" are Gaussian in the HE_11 field scale and need exactly one Marcatili HE_11 mode (SEP 19).
fullbool or None, optional: Force the azimuthal coordinate on or off. None decides from the modes; a slab set is never full.
zfloat, optional: Position along the waveguide the geometry is taken at, m.

Returns:

TransverseQuadrature: The rule.

Raises:

ModeError: If modes is empty, a field-scale rule is asked of anything but one Marcatili HE_11 mode, or the rule cannot be built.

Warns:

UserWarning: If the azimuthal rule is too coarse for a quartic product of these modes, which needs n_theta >= 4 h + 1 with h the highest azimuthal order, or if either axis of a Cartesian rule, or the x rule of a slab, is below cartesian_nodes() of the highest transverse order. The rule is still built, because a caller may know the products it will actually integrate; the warning is the only thing standing between an aliased Kerr term and a plausible-looking answer.

Examples

>>> from solprop.modes import MarcatiliMode
>>> from solprop.modes.quadrature import transverse_quadrature
>>> modes = [MarcatiliMode(125e-6, "He", 1.0, m=m) for m in (1, 2)]
>>> rule = transverse_quadrature(modes, n_r=24)
>>> rule.full, rule.n_nodes
(False, 24)