solprop.modes.capillary#

Marcatili modes of a gas-filled hollow capillary, with analytic derivatives.

The model is Marcatili and Schmeltzer’s (Bell Syst. Tech. J. 43, 1783-1809, 1964) asymptotic expansion for a hollow dielectric waveguide whose core radius is large compared with the wavelength. With \(k = \omega/c\), \(u_{nm}\) the Bessel zero of the mode, \(\varepsilon_{\rm co}\) the core permittivity and \(v_n\) the cladding term,

\[n_{\rm eff} = \sqrt{\varepsilon_{\rm co} - \left(\frac{u_{nm}}{ka}\right)^2\left(1 - \frac{i v_n}{ka}\right)^2}\]

is the "full" model, and expanding it to leading order in \(1/(ka)\) gives the "reduced" one,

\[n_{\rm eff} = 1 + \frac{\varepsilon_{\rm co} - 1}{2} - \frac{c^2u_{nm}^2}{2\omega^2a^2} + i\,\frac{c^3u_{nm}^2v_n}{a^3\omega^3}.\]

A positive imaginary part is loss (docs/dev/conventions.md section 5), so the reduced model’s loss term enters with a plus sign. Getting that backwards turns loss into gain silently, so there is one arithmetic kernel here and a test that asserts loss is loss.

Three properties of this implementation matter.

Analytic frequency derivatives. MarcatiliMode.dneff_domega() differentiates both models in closed form, chaining through the analytic derivative of the gas susceptibility and of the cladding index. Recomputing beta_1 every stage with a seven-point finite-difference stencil instead would cost about seven full index evaluations carrying 1e-10 of noise, which floors the cumulative operator at 1e-8 (see docs/dev/decisions.md).

No round trip through unity. MarcatiliMode.neff_minus_one() forms n_eff - 1 directly – x/(sqrt(1 + x) + 1) for the full model, term by term for the reduced one – rather than building 1 + chi/2 and subtracting one, which loses about eight digits at the 1e-4 susceptibility of a gas fill.

A factorisation the translation layer can use. z enters only through the core radius a(z) and the number density rho(z); everything frequency-dependent is MarcatiliTerms, sampled once. MarcatiliMode.sampled_terms() and MarcatiliMode.neff_from_terms() expose that split, so the Chebyshev builder can evaluate the operator at many z for two scalars each.

The manual page is docs/guide/waveguides/capillary-model.md.

Module Attributes

`ModeKind`	Which family a capillary mode belongs to: hybrid, transverse electric or magnetic.
`MarcatiliModel`	Whether to use the square-root expression or its leading-order expansion.
`WallImpedance`	Whether a coherent wall condenses one mean impedance or separate TE and TM ones.

Functions

`bessel_zero`(order, index)	Return the `index`-th positive zero of the Bessel function of the first kind.
`taper_connection`(modes[, z])	Return the taper connection of an `HE_1m` block at one position.
`taper_coupling_block`(modes)	Return the mode set the taper coupling is carried on, or `None`.
`vector_wall_operator`(modes, omega, *[, z, ...])	Build the condensed vector Robin wall generator of one azimuthal order.
`wall_connection`(modes[, z, impedance])	Return the moving-basis connection of a wall block at one position.
`wall_operator`(modes, omega, *[, z, active, ...])	Build the condensed Robin wall generator of a contiguous `HE_1m` block.
`wall_square_root_metric`(generator, wavenumber)	Return the square-root generator and its candidate flux metric, as a diagnostic.
`worst_wall_limits`(modes, omega, *, positions)	Return the validity ratios at the position where they are worst.

Classes

`CoherentWall`([condition_limit, ...])	Declare the coherent wall block of a capillary mode set.
`MarcatiliMode`(radius[, gas, pressure, kind, ...])	A Marcatili mode of a gas-filled hollow capillary.
`MarcatiliTerms`(omega, u_nm, ...)	Everything in a Marcatili effective index that depends only on frequency.
`VectorWallOperator`(omega, active, generator, ...)	The condensed vector Robin wall generator of one azimuthal order.
`WallFamily`(label, azimuthal_order, ...)	One radial family of a vector coherent wall block.
`WallLimits`(transverse_ratio, vacuum_ratio, ...)	The dimensionless validity ratios of one coherent wall block.
`WallOperator`(omega, active, generator, ...)	The condensed Robin wall generator of a contiguous `HE_1m` block.

class solprop.modes.capillary.CoherentWall(condition_limit=10000.0, validity_limit=0.25, sections=1, impedance='mean', method='frozen_residual', growth_limit=None, fallback=None)[source]#

Bases: Displayed

Declare the coherent wall block of a capillary mode set.

The capillary route’s diagonal operator gives every mode its own Marcatili attenuation and no coupling at all. The coherent wall instead condenses the outgoing Robin boundary condition onto the retained basis, which makes the wall a rank-one term that couples every mode to every other one with a definite phase. Declaring one on solprop.driver.simulation.Simulation.from_capillary() replaces the diagonal Marcatili loss of the block’s channels by that operator, exactly as wall_operator() builds it.

Its physical content, in the optics convention of docs/dev/conventions.md section 2 and with the carrier removed, is the paraxial generator

\[G = D + H,\qquad D = -\frac{i}{2ka^2}\operatorname{diag}(u_m^2), \qquad H = g(a)\,hh^{T},\]

with \(h_m = -\sqrt2\,\operatorname{sgn}J_1(u_m)\,u_m\) the radial derivative of the orthonormal basis function at the wall, and \(g(a) = i/[2ka^2(2N - i k_{\rm eff}a)]\). The manual page is docs/guide/waveguides/coherent-wall.md.

A guide that varies along z – a pressure gradient, a taper – is carried by a frozen-residual split: the complete generator G_s = D(z_s) + H(z_s) is frozen at each section’s midpoint inside the exponential flow, and its variation joins the HE1m connection

\[C_{mn} = \frac{a'}{a}\frac{2s_ms_nu_mu_n}{u_m^2-u_n^2}\ (m\ne n), \qquad C_{nn} = 0\]

in the explicit residual R_s(z) = D(z) - D(z_s) + H(z) - H(z_s) - C(z), which the interaction-picture Runge-Kutta source carries. Freezing is a rearrangement of the same equation and not a second-order approximation, because the residual is evaluated exactly; sections is what trades the residual’s stiffness against the number of frozen blocks (solprop.translate.coupled.frozen_residual_flow()).

What it replaces, and what it therefore refuses. The block is the run’s loss model, so a capillary carrying one must declare loss=True: loss=False, a scale or a measured curve is refused by name rather than discarded. It is derived for the full capillary model, so model="reduced" is refused too. It needs the contiguous HE_11 to HE_1N set, whose every channel declares the same radius and density profiles and whose profiles carry analytic derivatives, a scalar run and no noise model; all three backends execute the linear.coupled_flow leaf it becomes (method="cf4" on numpy_ref alone).

The vector boundary response is impedance="te_tm". The scalar block above condenses one impedance-averaged Robin condition, which is what a linearly polarised HE_1m field sees on azimuthal average; the vector block condenses the two conditions the wall actually imposes – the radial (TM-like) field sees k_o/n_r^2 and the tangential (TE-like) field sees k_o – and is therefore a block per total azimuthal order, coupling the two radial families of that order rather than one. See vector_wall_operator(); the averaged impedance is recovered from it as a limit rather than assumed.

Attributes:

impedancestr: "mean" for the impedance-averaged block, "te_tm" for the coupled azimuthal block, which gives the wall’s TE (tangential) and TM (radial) fields their own impedances. It names the wall’s two impedances and not the state’s polarisation axis: a "te_tm" run still carries one transverse component per channel. On a guide that varies along z – a taper, a pressure gradient – the vector block takes the same frozen or CF4 route as the scalar one, with a connection that is each radial family’s own and nothing between the families (wall_connection()). Which families the run carries is the mode set’s business and not this field’s: a set retaining one hybrid family alone is a legitimate vector block with no partner to couple to, and is the scalar block to the last bit, because that family’s Bessel order is zero and its diagonal impedance is the average. Declare the partner family – the hybrids of azimuthal index -nu – to have a conversion at all. The transverse order nu = 0 is not like that: TE_0m and TM_0m carry their own impedances whether or not the other is retained.
condition_limitfloat: Largest admitted norm(U) norm(U^-1) of one step’s flow action, dimensionless and at least one. The inverse action is what the interaction-picture source is carried back by, so an unbounded inverse is an unbounded amplification of that source; the backend measures the product per step and refuses above this. The default 1e4 admits a step that attenuates the most lossy retained mode by about nine e-foldings.
validity_limitfloat: Largest admitted value of the three dimensionless validity ratios of WallLimits. Above it the block is refused by name; above a tenth of it the construction warns, naming the ratio. The default 0.25 keeps the retained transverse scale well inside the paraxial and quasi-static regimes.
sectionsint: How many equal sections each interval of the run’s own linear-operator table is divided into, a positive integer. The placement is not this number’s: the edges are the table’s breakpoints, which already land on every declared kink of the radius and the density, subdivided sections times. One, the default, freezes one generator per declared interval; raising it reduces the residual’s norm and the inverse amplification, and is not permission to discard the residual. It has no effect on a constant guide, whose residual is exactly zero and which takes the single-section constant-Schur flow whatever this says.
methodstr: How a guide that varies along z is integrated. "frozen_residual", the default, is the frozen-residual split: each section’s generator frozen inside the exponential and its variation in the explicit residual the Runge-Kutta source carries, fifth order. "cf4" integrates the whole section generator by the fourth-order commutator-free exponential, with no explicit residual and so no advisory step, but at most fourth order and with a flow error the embedded pair does not see (solprop.translate.coupled.cf4_flow()). It is admitted only behind two guards: the internal-growth bound of the CF4 factors against growth_limit and the pair conditioning against condition_limit, for either impedance. It is refused by name on a guide that does not vary at all, whose flow the constant-Schur action already integrates exactly.
growth_limitfloat or None: Largest admitted internal-growth bound of a CF4 section, dimensionless and at least one; None takes solprop.translate.coupled.CF4_GROWTH_LIMIT. Declaring it without method="cf4" is refused, because nothing would read it.
fallbackstr or None: What a section the CF4 guards reject becomes: None refuses it by name, "frozen_residual" integrates it with the frozen-residual split and warns naming it. Declaring it without method="cf4" is refused.

Raises:

ModeError: If a limit or the section count is out of range, or if a method, a growth limit or a fallback is declared where nothing would honour it.

Examples

>>> from solprop.modes import CoherentWall
>>> CoherentWall().impedance, CoherentWall(impedance="te_tm").impedance
('mean', 'te_tm')
>>> repr(CoherentWall())[:43]
'CoherentWall(impedance=mean, condition_limi'

condition_limit: float = 10000.0#

fallback: Literal['frozen_residual'] | None = None#

growth_limit: float | None = None#

impedance: Literal['mean', 'te_tm'] = 'mean'#

method: Literal['frozen_residual', 'cf4'] = 'frozen_residual'#

sections: int = 1#

validity_limit: float = 0.25#

class solprop.modes.capillary.MarcatiliMode(radius, gas=None, pressure=0.0, *, kind='HE', n=1, m=1, phi=0.0, model='full', loss=True, cladding=None, temperature=None, susceptibility=None, mixture=None)[source]#

Bases: GasFillMixin

A Marcatili mode of a gas-filled hollow capillary.

Departs from the omega >= 0 of Mode in one place: a hollow capillary has no mode at zero frequency – the waveguide term (u c / (omega a))^2 diverges – so omega must be strictly positive here and a zero raises. The linear operator evaluates modes on the active band, which excludes the DC bin, so nothing in the package asks.

Both radius and pressure accept a number or a ZProfile, which is what makes a taper and a pressure gradient ordinary cases rather than special ones.

Attributes:

radiusfloat or ZProfile: Core radius in m, > 0. A number becomes a Constant.
gasGas or str or None: The fill species, or None for an evacuated capillary.
pressurefloat or ZProfile: Absolute pressure in bar, >= 0 (docs/dev/conventions.md section 1).
kindstr: "HE" (hybrid), "TE" or "TM". The transverse families require n == 0; hybrids require n != 0.
nint: Azimuthal index. May be negative for a hybrid mode.
mint: Radial index, >= 1: which zero of the Bessel function.
phifloat: Rotation of the mode about the axis in rad, positive from x towards y. It turns the whole transverse pattern: the Cartesian components of a hybrid mode are (cos((n - 1) theta - n phi), -sin((n - 1) theta - n phi)) times J_{n-1}(u r / a). For HE_1m the pattern is uniform, so phi is the polarisation angle measured from the x axis: phi = 0 is x-polarised and phi = pi / 2 is y-polarised. For n >= 2 an oblique phi is a rotated pattern and not a polarisation angle – the field points in different directions at different theta – and the degenerate partner of phi = 0 is the quarter turn phi = pi / (2 n), which orthogonal_polarisation returns.
modelstr: "full" for the square root, "reduced" for its leading-order expansion.
lossbool or float or callable: True for the model’s own loss, False for none, a float to scale it, or a callable alpha(omega) giving an absolute power loss coefficient in 1/m that replaces it. The float and the callable therefore mean different things, deliberately: the float is a fitting knob on the model, the callable is measured fibre loss. The model’s own loss is never negative: where the square root lands on the branch that reads as gain – a cladding index below unity, outside Marcatili’s regime – it is clamped at zero and the mode logs the fact once. A callable is the caller’s and is not clamped. A callable that carries a derivative attribute gets an analytic d(Im n_eff)/d(omega); otherwise that one term is differenced, and the constructor says so.
claddingGlassIndex or complex or None: The cladding index. None uses fused_silica(); a plain complex number is a constant index, which is what a comparison against a published or symbolic result needs.
temperaturefloat, ZProfile or None: Temperature of the fill in K, > 0, or a path of it along z in m, which the density follows (DensityPath); a constant path is stored as its number. None (the default) is 293.15 K for a gas and the mixture’s own beside mixture, which carries its temperature: a different one given beside it is refused by name. Resolved at construction, so the field is never None afterwards.
susceptibilityGasSusceptibility or None: The gas susceptibility per unit number density. None takes the built-in gas’s, optmat’s, which is what a run normally wants. Passing one is the gas-side counterpart of a constant or tabulated cladding: it replaces the refractive data while leaving the equation of state alone, so gas is still what sets the density along z. A comparison against a published result or another code whose gas model differs can give the mode the other model’s susceptibility and isolate what it is actually comparing.
mixtureMixture or None: Labelled populations with their own density paths. Supplies temperature and optical data; do not also pass gas, pressure or susceptibility.

Examples

>>> import numpy as np
>>> from solprop.modes import MarcatiliMode
>>> mode = MarcatiliMode(125e-6, "He", 1.0, model="reduced")
>>> omega = np.array([2.356e15])
>>> f"{complex(mode.neff(omega)[0]):.9f}"
'1.000028936+0.000000009j'
>>> f"{float(mode.aeff()):.4e}"
'2.3393e-08'

aeff(z=0.0)[source]#

Return the effective area of the mode.

Parameters:

zfloat, optional: Position along the waveguide, m.

Returns:

float: A_eff = (int |e|^2 dA)^2 / int |e|^4 dA in m^2, in closed form: the transverse profile is a Bessel function of r / a, so the whole radius dependence is the prefactor a^2.

Examples

>>> from solprop.modes import MarcatiliMode
>>> f"{MarcatiliMode(75e-6, 'He', 1.0).aeff():.6e}"
'8.421575e-09'

at_parameters(values)[source]#

Return the same mode with the named paths frozen at the given values.

The effective index reaches z only through the two paths, so replacing the two paths by constants gives a mode that is exactly this one at every position where the paths take those values. That is what lets solprop.linear.dispersion_operator() be evaluated at parameter values, and so what lets the cumulative operator be expanded in the parameters rather than tabulated along z (docs/guide/numerics/linear-operator.md).

The copy shares the parent’s derived frequency-only quantities – the Bessel zero, the cladding index, the effective-area factor – so freezing costs nothing per evaluation. Its radius field follows core_radius; its pressure field still records the parent’s declaration, because the density is given directly here and the equation of state is not inverted. Nothing reads pressure after construction; density is the path that is evaluated.

Parameters:

valuescollections.abc.Mapping: Any subset of {"core_radius": a in m, "density": rho in m^-3}. A name left out keeps its path.

Returns:

MarcatiliMode: A copy reporting z_constant.

Raises:

ModeError: If a name is not a declared path, the core radius is not positive, or the density is negative.

Examples

>>> import numpy as np
>>> from solprop.modes import MarcatiliMode
>>> from solprop.profiles import PressureGradient
>>> fill = PressureGradient.two_point(1.0, 5.0, 1.0)
>>> mode = MarcatiliMode(125e-6, "He", fill)
>>> frozen = mode.at_parameters(
...     {
...         name: float(path(np.array([0.3]))[0])
...         for name, path in mode.parameter_paths.items()
...     }
... )
>>> frozen.z_constant
True
>>> omega = np.array([2.356e15])
>>> bool(np.array_equal(frozen.neff(omega), mode.neff(omega, 0.3)))
True

property azimuthal_order: int[source]#

Return the highest azimuthal harmonic in the Cartesian field components.

Returns:

int: |n - 1| for a hybrid mode and 1 for a transverse one. The signed n - 1 would be negative for HE_{-1,m} and would then understate the theta rule a cubic product needs.

property channel_label: ChannelLabel[source]#

Return the family, the two indices and which polarisation partner this is.

Returns:

ChannelLabel

parity is None for a transverse family, which has no degenerate polarisation partner. For a hybrid mode it follows the one definition in ChannelLabel: the even member is the cos((n - 1) theta) pattern about the x axis, which here is the unrotated mode phi = 0 – at that angle the transverse components are (cos((n - 1) theta), -sin((n - 1) theta)), the same pattern the exact step-index even member carries, and for HE_1m it is the x-polarised state. The quarter turn phi = pi / (2 n) is therefore the odd member, y-polarised for HE_1m. Any other rotation is a gauge rather than a parity: its parity is None and its display label carries the angle, so that two rotated copies of one mode are still told apart.

A negative azimuthal index is the EH family, and the label says so: the label’s own azimuthal_order is the non-negative total azimuthal order |n| and the family is "EH", which is how an exact StepIndexMode of the same physical channel labels itself, so an EH capillary mode can be put in a ModeCollection.

Examples

>>> from solprop.modes import MarcatiliMode
>>> mode = MarcatiliMode(125e-6, "He", 1.0)
>>> mode.channel_label.display_label, mode.principal_component()
('HE1,1:even', 'x')
>>> partner = mode.orthogonal_polarisation
>>> partner.channel_label.display_label, partner.principal_component()
('HE1,1:odd', 'y')
>>> MarcatiliMode(125e-6, "He", 1.0, n=-1).channel_label.display_label
'EH1,1:even'

cladding: GlassIndex | complex | None = None#

cladding_index: GlassIndex | complex#

core_index_minus_one(omega, z=0.0)[source]#

Return the filled core’s n_core - 1 formed without cancellation.

n_core = sqrt(1 + chi) with chi the fill’s susceptibility at the position’s density, and n_core - 1 is formed as chi / (sqrt(1 + chi) + 1), whose denominator is near two, so it carries the relative rounding of chi itself rather than the absolute rounding of unity. Its difference with neff_minus_one() is n_core - n_eff to the rounding of two numbers of order chi, where forming k - beta from the two wavenumbers would carry the rounding of k.

Parameters:

omeganumpy.ndarray: Angular frequencies, rad/s, shape (n_omega,), strictly positive.
zfloat, optional: Position along the waveguide, m.

Returns:

numpy.ndarray: n_core - 1, complex, shape (n_omega,); zero for an evacuated core.

Examples

>>> import numpy as np
>>> from solprop.modes import MarcatiliMode
>>> mode = MarcatiliMode(125e-6, "He", 1.0)
>>> omega = np.array([2.356e15])
>>> core = mode.core_index_minus_one(omega)
>>> bool(core.real[0] > mode.neff_minus_one(omega).real[0] > 0.0)
True

core_radius: ZProfile#

density: ZProfile#

dimlimits(z=0.0)[source]#

Return the cross-section parametrisation and the integration limits.

Parameters:

zfloat, optional: Position along the waveguide, m.

Returns:

tuple: ("polar", (0, 0), (a(z), 2 pi)) in m and rad.

Examples

>>> from solprop.modes import MarcatiliMode
>>> MarcatiliMode(125e-6, 'He', 1.0).dimlimits()[0]
'polar'

dneff_domega(omega, z=0.0)[source]#

Return the analytic frequency derivative of the effective index.

Both models are differentiated in closed form, chaining through the analytic derivative of the gas susceptibility and of the cladding index (ADR-006). The only term that may be differenced is the imaginary part when the loss is given as a callable without its own derivative.

Parameters:

omeganumpy.ndarray: Angular frequencies, rad/s, shape (n_omega,), strictly positive.
zfloat, optional: Position along the waveguide, m.

Returns:

numpy.ndarray: dn_eff/domega in s/rad, complex, shape (n_omega,).

Examples

>>> import numpy as np
>>> from solprop.modes import MarcatiliMode
>>> mode = MarcatiliMode(125e-6, "He", 1.0)
>>> f"{mode.dneff_domega(np.array([2.356e15]))[0].real:.4e}"
'2.6423e-21'

dneff_domega_from_terms(terms, radius, density)[source]#

Return the analytic dn_eff/domega from pre-sampled terms.

Parameters:

termsMarcatiliTerms: The sampled frequency-dependent quantities.
radiusfloat: Core radius in m, > 0.
densityfloat: Number density of the fill in m^-3, >= 0.

Returns:

numpy.ndarray: dn_eff/domega in s/rad, complex.

Examples

>>> import numpy as np
>>> from solprop.modes import MarcatiliMode
>>> mode = MarcatiliMode(125e-6, "He", 1.0, model="reduced")
>>> terms = mode.sampled_terms(np.array([2.356e15]))
>>> f"{mode.dneff_domega_from_terms(terms, 125e-6, 2.46e25).real[0]:.4e}"
'2.6420e-21'

dneff_dz(omega, z=0.0)[source]#

Return the analytic derivative of the effective index along the capillary.

z reaches the index through exactly two scalars, the core radius and the fill’s number density, so dn/dz = (dn/da) a' + sum_s (dn/drho_s) rho_s' with the declared profiles’ own analytic slopes. For the full model n^2 = 1 + chi - W with chi = sum_s gamma_s rho_s and W = (u s)^2 (1 - i v s)^2, s = c/(omega a), so that dW/da = -(2 u^2 s^2/a)(1 - i v s)(1 - 2 i v s) and dn/dz = (sum_s gamma_s rho_s' - (dW/da) a')/(2 n); the reduced model differentiates its own polynomial. The loss specification acts on the imaginary part exactly as it does on the index: False drops it, a float scales it, a measured curve alpha(omega) does not depend on z and contributes nothing, where the no-gain clamp binds the imaginary part’s slope is zero, and where the evanescent floor binds the whole slope is zero, as dneff_domega()’s is. What the coherent wall route reads is the real part, the slope of the power weight w = Re n_eff.

Parameters:

omeganumpy.ndarray: Angular frequencies, rad/s, shape (n_omega,), strictly positive.
zfloat, optional: Position along the waveguide, m.

Returns:

numpy.ndarray: dn_eff/dz in 1/m, complex, shape (n_omega,).

Raises:

ModeError: If the mode is a subclass that computes its index another way – a structured cladding or a tabulated index, whose z derivative is not this one – or a declared profile reports no analytic derivative.

Examples

>>> import numpy as np
>>> from solprop.modes import MarcatiliMode
>>> from solprop.profiles import Linear
>>> mode = MarcatiliMode(Linear(50e-6, 40e-6, 1e-3), "He", 1.0)
>>> bool(mode.dneff_dz(np.array([2.356e15]), 5e-4)[0].real < 0.0)
True

field(points, z=0.0)[source]#

Return the transverse field components at points of the cross-section.

Parameters:

pointsnumpy.ndarray: Shape (n_points, 2): (r, theta) in m and rad.
zfloat, optional: Position along the waveguide, m.

Returns:

numpy.ndarray: (E_x, E_y) of the unnormalised profile, shape (n_points, 2), real. A hybrid mode has magnitude J_{n-1}(u r/a) independent of theta and azimuthal order |n - 1| in its Cartesian components; at phi = 0 those components are (cos((n-1) theta), -sin((n-1) theta)), the even pattern of every family, so HE_{11} is uniformly polarised along x there.

Examples

>>> import numpy as np
>>> from solprop.modes import MarcatiliMode
>>> mode = MarcatiliMode(125e-6, 'He', 1.0)
>>> mode.field(np.array([[0.0, 0.0]])).round(6)
array([[1., 0.]])

gas: Gas | str | None = None#

kind: Literal['HE', 'TE', 'TM'] = 'HE'#

loss: bool | float | Callable[[Float[ndarray, 'n_omega']], Float[ndarray, 'n_omega']] = True#

m: int = 1#

mixture: Mixture | None = None#

model: Literal['full', 'reduced'] = 'full'#

n: int = 1#

neff(omega, z=0.0)[source]#

Return the complex effective index.

Parameters:

omeganumpy.ndarray: Angular frequencies, rad/s, shape (n_omega,), strictly positive.
zfloat, optional: Position along the waveguide, m.

Returns:

numpy.ndarray: n_eff(omega, z), complex, shape (n_omega,); Im n_eff > 0 is loss.

Examples

>>> import numpy as np
>>> from solprop.modes import MarcatiliMode
>>> mode = MarcatiliMode(125e-6, "He", 1.0)
>>> bool(mode.neff(np.array([2.356e15]))[0].imag > 0.0)
True

neff_from_terms(terms, radius, density)[source]#

Return the effective index from pre-sampled terms and two scalars.

Parameters:

termsMarcatiliTerms: The sampled frequency-dependent quantities, from sampled_terms().
radiusfloat: Core radius in m at the position of interest, > 0.
densityfloat: Number density of the fill in m^-3, >= 0.

Returns:

numpy.ndarray: n_eff, complex, shape (n_omega,); Im n_eff > 0 is loss.

Examples

>>> import numpy as np
>>> from solprop.modes import MarcatiliMode
>>> mode = MarcatiliMode(125e-6, "He", 1.0)
>>> terms = mode.sampled_terms(np.array([2.356e15]))
>>> value = mode.neff_from_terms(terms, 125e-6, 2.46e25)
>>> bool(value.imag[0] > 0.0)
True

neff_minus_one(omega, z=0.0)[source]#

Return n_eff - 1 formed without cancellation.

Parameters:

omeganumpy.ndarray: Angular frequencies, rad/s, shape (n_omega,), strictly positive.
zfloat, optional: Position along the waveguide, m.

Returns:

numpy.ndarray: n_eff - 1, complex, exact to rounding for a gas-filled capillary.

Examples

>>> import numpy as np
>>> from solprop.modes import MarcatiliMode
>>> mode = MarcatiliMode(125e-6, "He", 1.0)
>>> f"{mode.neff_minus_one(np.array([2.356e15]))[0].real:.4e}"
'2.8935e-05'

normalisation(z=0.0)[source]#

Return the power carried by the unnormalised transverse profile.

Parameters:

zfloat, optional: Position along the waveguide, m.

Returns:

float: N = (1/2) sqrt(eps_0/mu_0) int |e|^2 dA in W, in closed form (pi/2) a^2 J_{n'}(u)^2 sqrt(eps_0/mu_0) with n' = n for a hybrid mode and 2 for a transverse one. The vacuum-normalised basis gives monochromatic power Re(n_eff) * abs(a)**2 for internal amplitude a (docs/dev/conventions.md section 4).

Examples

>>> from solprop.modes import MarcatiliMode
>>> f"{MarcatiliMode(125e-6, 'He', 1.0).normalisation():.6e}"
'1.755864e-11'

property orthogonal_polarisation: MarcatiliMode | None[source]#

Return the degenerate partner of a hybrid mode, or None.

Returns:

MarcatiliMode or None: The same mode turned by pi / (2 n), which turns the pattern a further quarter turn from x towards y and so exchanges the two Cartesian components: for HE_11 that is the y-polarised partner of the x-polarised default, and so the odd member of the pair (channel_label). None for a transverse family, whose members are not degenerate in polarisation – TE_0m and TM_0m are different modes, not two states of one.

Examples

>>> from solprop.modes import MarcatiliMode
>>> partner = MarcatiliMode(125e-6, "He", 1.0).orthogonal_polarisation
>>> f"{partner.phi:.6f}"
'1.570796'
>>> MarcatiliMode(125e-6, kind="TE", n=0).orthogonal_polarisation is None
True

property parameter_paths: Mapping[str, ZProfile][source]#

Return the scalar paths that carry the whole z dependence.

Returns:

collections.abc.Mapping: Radius and the material paths, including constant ones. One component retains density; a shared fill uses total pressure in bar; independent populations use density:<label> in m^-3. See solprop.modes.protocol.ParameterisedMode.

Examples

>>> from solprop.modes import MarcatiliMode
>>> sorted(MarcatiliMode(125e-6, "He", 1.0).parameter_paths)
['core_radius', 'density']

phi: float = 0.0#

pressure: float | ZProfile = 0.0#

principal_component(z=0.0)[source]#

Return the sole nonzero Cartesian component, or None for a vector profile.

Parameters:

zfloat, optional: Position along the waveguide, m; the polarisation is independent of it.

Returns:

{‘x’, ‘y’} or None: An axis-aligned HE_1m component, or None when both components are nonzero. phi = 0 is x (ADR-031) and the quarter turn is y. Angles are classified exactly modulo pi, without a tolerance that would erase a small physical component.

radius: float | ZProfile#

property radius_is_constant: bool[source]#

Return whether the transverse profile is the same at every z.

Returns:

bool: True when the core radius does not vary, so the synthesis and projection matrices of the multimode transform may be built once. A pressure gradient changes the effective index but not the field profile.

refractivity_from_terms(terms, radius, density)[source]#

Return n_eff - 1 from pre-sampled terms and two scalars.

Parameters:

termsMarcatiliTerms: The sampled frequency-dependent quantities.
radiusfloat: Core radius in m, > 0.
densityfloat: Number density of the fill in m^-3, >= 0.

Returns:

numpy.ndarray: n_eff - 1, complex, formed without cancellation.

Examples

>>> import numpy as np
>>> from solprop.modes import MarcatiliMode
>>> mode = MarcatiliMode(125e-6, "He", 1.0, model="reduced")
>>> terms = mode.sampled_terms(np.array([2.356e15]))
>>> f"{complex(mode.refractivity_from_terms(terms, 125e-6, 2.46e25)[0]):.3e}"
'2.881e-05+9.000e-09j'

rotated(phi)[source]#

Return the same mode turned about the axis.

Parameters:

phifloat: The new rotation angle in rad, in the sense of phi: positive from x towards y, so for HE_1m it is the polarisation angle measured from the x axis.

Returns:

MarcatiliMode: A copy with the new angle; everything else is unchanged.

Examples

>>> import math
>>> from solprop.modes import MarcatiliMode
>>> MarcatiliMode(125e-6, "He", 1.0).rotated(math.pi / 2).phi
1.5707963267948966

sampled_terms(omega)[source]#

Sample everything in the effective index that does not depend on z.

Parameters:

omeganumpy.ndarray: Angular frequencies in rad/s, shape (n_omega,), > 0.

Returns:

MarcatiliTerms: The gas susceptibility, the cladding term and their derivatives, ready to be combined with a radius and a density at any position.

Raises:

ModeError: If any frequency is not strictly positive or not finite.

Examples

>>> import numpy as np
>>> from solprop.modes import MarcatiliMode
>>> terms = MarcatiliMode(125e-6, "He", 1.0).sampled_terms(np.array([2.4e15]))
>>> terms.cladding_term.shape
(1,)

property scale_invariant: bool[source]#

Return whether the profile scales exactly with the core radius.

Returns:

bool: Always True: the field depends on r only through r / a and the normalisation goes as a^2, so a taper rescales the transform matrices by two scalars.

susceptibility: GasSusceptibility | None = None#

temperature: float | ZProfile | None = None#

u_nm: float#

property z_constant: bool[source]#

Return whether nothing about the mode varies along z.

Returns:

bool: True when neither the core radius nor the density varies. The multimode transform cares only about the radius, which radius_is_constant reports separately.

solprop.modes.capillary.MarcatiliModel#

Whether to use the square-root expression or its leading-order expansion.

alias of Literal[‘full’, ‘reduced’]

class solprop.modes.capillary.MarcatiliTerms(omega, u_nm, core_susceptibility, dcore_susceptibility_domega, cladding_term, dcladding_term_domega)[source]#

Bases: object

Everything in a Marcatili effective index that depends only on frequency.

z enters a capillary mode through exactly two scalars, the core radius and the number density, so the frequency-dependent pieces can be sampled once and reused at every position (docs/dev/architecture.md section 6.5). Building this is the expensive part: it calls the gas model and the cladding model.

The frequency is typed as inexact rather than real on purpose. The effective index is a holomorphic function of omega – powers, a square root and a rational cladding term – so sampling these terms at a complex frequency lets the analytic derivative be checked against Cauchy’s integral formula, which is free of the subtractive cancellation that limits any finite difference. Nothing in the public API does that; MarcatiliMode.sampled_terms() requires a real, positive frequency.

Attributes:

omeganumpy.ndarray: Angular frequencies in rad/s, shape (n_omega,), > 0 in every use except the contour check described above.
u_nmfloat: The Bessel zero of the mode, dimensionless.
core_susceptibilitynumpy.ndarray: Per-particle susceptibility gamma(omega) in m^3, so that chi = gamma * rho. Mixtures carry shape (n_species, n_omega); the existing single-gas spelling retains its frequency vector. Values are real for transparent built-in gases and complex for absorbing custom species.
dcore_susceptibility_domeganumpy.ndarray: Its analytic frequency derivative in m^3 s/rad.
cladding_termnumpy.ndarray: The complex v_n(omega) of the mode’s family.
dcladding_term_domeganumpy.ndarray: Its analytic frequency derivative in s/rad.

cladding_term: Complex[ndarray, 'n_clad']#

core_susceptibility: Inexact[ndarray, '*species n_terms']#

dcladding_term_domega: Complex[ndarray, 'n_clad']#

dcore_susceptibility_domega: Inexact[ndarray, '*species n_terms']#

omega: Inexact[ndarray, 'n_omega']#

u_nm: float#

solprop.modes.capillary.ModeKind#

Which family a capillary mode belongs to: hybrid, transverse electric or magnetic.

alias of Literal[‘HE’, ‘TE’, ‘TM’]

class solprop.modes.capillary.VectorWallOperator(omega, active, generator, families, sector, family_traces, weight, wavenumber, impedance, limits)[source]#

Bases: Displayed

The condensed vector Robin wall generator of one azimuthal order.

The output of vector_wall_operator(), and the vector counterpart of WallOperator: the dense paraxial generator in the optics convention and in 1/m, the pieces it was factored from and the validity ratios it was checked against. Crossing to the internal FFT convention is the complex conjugate of the whole equation, exactly as for the scalar block.

The wall term is rank two, not rank one: one boundary trace vector per family, weighted by the condensed 2 x 2 impedance of the order.

Attributes:

omeganumpy.ndarray: The angular frequencies in rad/s the block is tabulated on, shape (n_omega,); the whole grid, with the inactive bins left exactly zero.
activenumpy.ndarray: Where the operator was evaluated, bool, shape (n_omega,).
generatornumpy.ndarray: G = D + H in 1/m, complex, shape (n_modes, n_modes, n_omega), in the optics convention with the carrier k removed.
familiestuple of WallFamily: The one or two radial families of this order, in the state’s channel order.
sectorint: The total azimuthal order nu >= 0 every channel of the block belongs to.
family_tracesnumpy.ndarray: The boundary traces h_m arranged one column per family, dimensionless, shape (n_modes, n_families), exactly zero outside each family’s own channels. The wall term is V W V^T with V this matrix.
weightnumpy.ndarray: W = i S/(2 k a^2) in 1/m, complex, shape (n_families, n_families, n_omega): the condensed impedance of the order, scaled by the paraxial reduction. Zero off the active band.
wavenumbernumpy.ndarray: The in-core wavenumber k = n_core omega / c in 1/m, complex, shape (n_omega,), zero off the active band.
impedancenumpy.ndarray: The wall’s own two impedances in m, complex, shape (2, n_omega): the radial (TM-like) v_TM/k_0 and the tangential (TE-like) v_TE/k_0, built from this package’s own Marcatili cladding terms. Their mean is the scalar block’s 1/k_eff and their half-difference is the anisotropy.
limitsWallLimits: The validity ratios over the active band.

Examples

>>> import numpy as np
>>> from solprop.modes import MarcatiliMode, vector_wall_operator
>>> modes = [MarcatiliMode(125e-6, "He", 1.0, kind=k, n=0) for k in ("TM", "TE")]
>>> block = vector_wall_operator(modes, np.array([2.356e15]))
>>> block.sector, block.n_families, block.generator.shape
(0, 2, (2, 2, 1))
>>> bool(abs(block.generator[0, 1, 0]) == 0.0)
True

active: Bool[ndarray, 'n_omega_vec_a']#

families: tuple[WallFamily, ...]#

family_traces: Float[ndarray, 'n_modes_vec_h n_families_vec_h']#

generator: Complex[ndarray, 'n_modes_vec_out n_modes_vec_in n_omega_vec_g']#

impedance: Complex[ndarray, 'n_wall_components n_omega_vec_z']#

limits: WallLimits#

property n_families: int[source]#

Return how many radial families the block carries.

Returns:

int: One or two.

property n_modes: int[source]#

Return how many channels the block couples.

Returns:

int: The total over its families.

omega: Float[ndarray, 'n_omega_vec']#

sector: int#

wavenumber: Complex[ndarray, 'n_omega_vec_k']#

weight: Complex[ndarray, 'n_families_w_out n_families_w_in n_omega_vec_w']#

class solprop.modes.capillary.WallFamily(label, azimuthal_order, bessel_order, n_channels, angular, radial_sign)[source]#

Bases: object

One radial family of a vector coherent wall block.

A vector block is one total azimuthal order \(\nu\) of a capillary, and that order has at most two radial families, because the transverse field of an order-\(\nu\) mode is a combination of the two circular components \(\hat e_\pm\), whose Cartesian angular harmonics are \(\nu\mp1\) and whose radial functions are therefore \(J_{\nu-1}\) and \(J_{\nu+1}\). Each family is a contiguous m = 1 .. N_f set of this package’s own Marcatili channels.

Attributes:

labelstr: How the family prints: "HE_1m", "EH_1m", "TE_0m", "TM_0m" and the higher orders in the same spelling.
azimuthal_orderint: The channels’ own MarcatiliMode.n. Within one sector the two families carry +nu and -nu; nu = 0 carries the two transverse kinds.
bessel_orderint: p, the order of the Bessel function the family’s radial basis is built from, >= 0: |n - 1| for a hybrid family and 1 for a transverse one.
n_channelsint: How many radial orders of this family the block retains, >= 1.
angulartuple of float: The family’s angular content at the wall, (radial, tangential), a unit two-vector: the fraction of its boundary field that is normal to the wall and the fraction tangential to it. (1, 0) for TM_0m and (0, 1) for TE_0m, which is why those two see the wall’s own two impedances; and (1, -1)/sqrt 2 and (1, 1)/sqrt 2 for the two hybrid families of an order nu >= 1, which is why those see the impedance average on the diagonal and the impedance difference between them.
radial_signfloat: The sign of the family’s implemented radial profile against J_p, +1 or -1. MarcatiliMode.field() builds a hybrid from jv(n - 1, .), and J_{-p} = (-1)^p J_p, so the n = -nu family of an even order is implemented with the opposite sign to J_{nu+1}. It is a gauge of the mode and cancels out of every diagonal quantity, but it is the sign of the coupling between the two families, so the wall’s boundary traces carry it.

Examples

>>> from solprop.modes.capillary import WallFamily
>>> family = WallFamily("TE_0m", 0, 1, 4, (0.0, 1.0), 1.0)
>>> family.truncation
9

angular: tuple[float, float]#

azimuthal_order: int#

bessel_order: int#

label: str#

n_channels: int#

radial_sign: float#

property truncation: int[source]#

Return the family’s boundary truncation constant 2 N_f + p_f.

Returns:

int: The constant the condensation’s denominator carries for this family: the truncated boundary self term sum_m h_m^2/u_m^2 = 2 N_f of the retained orders plus the p_f of the exact recurrence q J_{p+1}(q) = (p - i eta) J_p(q). Dropping the p_f leaves a residue that does not fall with N.

solprop.modes.capillary.WallImpedance#

Whether a coherent wall condenses one mean impedance or separate TE and TM ones.

alias of Literal[‘mean’, ‘te_tm’]

class solprop.modes.capillary.WallLimits(transverse_ratio, vacuum_ratio, quasi_static, band, n_modes, radius)[source]#

Bases: Displayed

The dimensionless validity ratios of one coherent wall block.

Each is evaluated with the highest retained mode over the band the run scores (band below), which is where every approximation in it is worst there. None of them is a proof of Maxwell accuracy: they are the stated conditions under which the condensed paraxial model is the model it claims to be. In particular a converged paraxial Robin equation is not by itself Maxwell accuracy, and many retained basis components can exceed the in-core wavenumber without being propagating channels.

Attributes:

transverse_ratiofloat: u_N / (|k| a) with k the in-core wavenumber: how far the highest retained transverse scale is from paraxial. The paraxial generator’s own error is a quarter of its fourth power relative to the carrier.
vacuum_ratiofloat: u_N / (k_0 a) with k_0 the vacuum wavenumber: the grazing-incidence ratio the Marcatili expansion is asymptotic in, which the cladding term is referred to.
quasi_staticfloat: 2N / |a k_eff|: the fraction of the eliminated boundary coordinate the condensation discards in the quasi-static wall limit. Raising N raises this, so it cannot by itself validate the condensation. A vector block takes the same ratio family by family and keeps the worst: (2N_f + p_f) sum_g |Z_fg| / a, with p_f the family’s Bessel order and Z its impedance matrix in m, which is the row sum the condensed S = (C - i a Z^-1)^-1 departs from i Z/a by (vector_wall_operator()).
bandtuple of float: Lowest and highest angular frequency in rad/s the ratios were taken over. On a run this is the error norm’s scored band (solprop.grid.TimeGrid.inband_mask, the requested wavelength limits), not the active band the operator is built on: a field grid’s active band reaches about twice the longest requested wavelength, where every ratio (each proportional to the wavelength) is about twice its in-band value, and those collar bins are absorbed rather than scored. Without a scored mask, or with one that holds no evaluated bin, it is the whole evaluated band.
n_modesint: How many channels the block carries, over every family in it.
radiusfloat: Core radius in m the ratios were taken at.

Examples

>>> from solprop.modes.capillary import WallLimits
>>> limits = WallLimits(0.01, 0.01, 0.02, (1e15, 3e15), 4, 125e-6)
>>> f"{limits.worst:.3f}"
'0.020'

band: tuple[float, float]#

n_modes: int#

named()[source]#

Return the name and the value of the ratio that is worst.

Returns:

tuple: The ratio’s name as it is written here, and its value.

Examples

>>> from solprop.modes.capillary import WallLimits
>>> WallLimits(0.3, 0.1, 0.2, (1e15, 3e15), 4, 125e-6).named()[0]
'transverse_ratio'

quasi_static: float#

radius: float#

transverse_ratio: float#

vacuum_ratio: float#

property worst: float[source]#

Return the largest of the three ratios.

Returns:

float: The one a refusal or a warning is taken on.

class solprop.modes.capillary.WallOperator(omega, active, generator, trace_vector, coupling, wavenumber, impedance, limits)[source]#

Bases: Displayed

The condensed Robin wall generator of a contiguous HE_1m block.

The output of wall_operator(): the dense paraxial generator in the optics convention and in 1/m, together with the factored pieces it was built from and the validity ratios it was checked against. Crossing to the internal FFT convention is the complex conjugate of the whole equation, c' = G c becoming c' = conj(G) c, which is what the translation layer stores; it is not -G, and a passive mode propagated with the wrong one grows (docs/dev/conventions.md section 2, ADR-017).

Attributes:

omeganumpy.ndarray: The angular frequencies in rad/s the block is tabulated on, shape (n_omega,); the whole grid, with the inactive bins left exactly zero.
activenumpy.ndarray: Where the operator was evaluated, bool, shape (n_omega,).
generatornumpy.ndarray: G = D + H in 1/m, complex, shape (n_modes, n_modes, n_omega), in the optics convention with the carrier k removed.
trace_vectornumpy.ndarray: h_m, dimensionless, shape (n_modes,): the radial derivative of the m-th orthonormal basis function at the wall.
couplingnumpy.ndarray: g(a) in 1/m, complex, shape (n_omega,), zero off the active band.
wavenumbernumpy.ndarray: The in-core wavenumber k = n_core omega / c in 1/m, complex, shape (n_omega,), zero off the active band.
impedancenumpy.ndarray: k_eff = k_0 / v_n in 1/m, complex, shape (n_omega,), zero off the active band: the vacuum wavenumber over the current filled-core cladding term, which is what makes the isolated block reproduce this package’s own capillary loss.
limitsWallLimits: The validity ratios over the active band.
n_modesint: Return how many channels the block couples.

Examples

>>> import numpy as np
>>> from solprop.modes import MarcatiliMode, wall_operator
>>> modes = [MarcatiliMode(125e-6, "He", 1.0, m=index) for index in (1, 2)]
>>> block = wall_operator(modes, np.array([2.356e15]))
>>> block.generator.shape
(2, 2, 1)
>>> bool(block.generator[0, 0, 0].real < 0.0)
True

active: Bool[ndarray, 'n_omega_wall_a']#

coupling: Complex[ndarray, 'n_omega_wall_g2']#

generator: Complex[ndarray, 'n_modes_wall_out n_modes_wall_in n_omega_wall_g']#

impedance: Complex[ndarray, 'n_omega_wall_e']#

limits: WallLimits#

property n_modes: int[source]#

Return how many channels the block couples.

Returns:

int: The length of the trace vector.

omega: Float[ndarray, 'n_omega_wall']#

trace_vector: Float[ndarray, 'n_modes_wall_h']#

wavenumber: Complex[ndarray, 'n_omega_wall_k']#

solprop.modes.capillary.bessel_zero(order, index)[source]#

Return the index-th positive zero of the Bessel function of the first kind.

Parameters:

orderint: Order of the Bessel function. Negative orders are folded onto positive ones, since J_-p = (-1)^p J_p for integer p and the zeros coincide.
indexint: Which zero, counting from 1.

Returns:

float: The zero j_{order, index}, dimensionless.

Raises:

ModeError: If index is not a positive integer.

Examples

>>> from solprop.modes.capillary import bessel_zero
>>> f"{bessel_zero(0, 1):.10f}"
'2.4048255577'
>>> bessel_zero(-2, 1) == bessel_zero(2, 1)
True

solprop.modes.capillary.taper_connection(modes, z=0.0)[source]#

Return the taper connection of an HE_1m block at one position.

The connection \(C_{mn} = \langle\psi_m|\partial_z\psi_n\rangle\) of the field-normalised J_0(u_m r/a) basis, exactly wall_connection()’s scalar block (impedance="mean") – real, antisymmetric, proportional to a'/a – without its two conditions that belong to the wall: it admits any loss and either capillary model, because the basis is the Dirichlet one whatever the modes’ effective index says. It is what the taper coupling carries on a multimode capillary without a wall.

Parameters:

modessequence of Mode: The contiguous HE_11 to HE_1N set in state order, every channel declaring the same radius and density profiles, each with an analytic derivative.
zfloat, optional: Position along the capillary in m.

Returns:

numpy.ndarray: C in 1/m, real, shape (n_modes, n_modes), with a zero diagonal.

Raises:

ModeError: If the set is not such a block.

Examples

>>> import numpy as np
>>> from solprop.modes.capillary import MarcatiliMode, taper_connection
>>> from solprop.profiles import Linear
>>> taper = Linear(125e-6, 62.5e-6, 1e-2)
>>> modes = [MarcatiliMode(taper, "He", 1.0, m=m, loss=False) for m in (1, 2)]
>>> link = taper_connection(modes, 5e-3)
>>> bool(np.all(link == -link.T)), float(link[0, 0])
(True, 0.0)

solprop.modes.capillary.taper_coupling_block(modes)[source]#

Return the mode set the taper coupling is carried on, or None.

The taper coupling – the moving-basis connection <psi_m|d_z psi_n> of taper_connection() in a coupled flow, with each channel’s own propagation constant on the flow’s diagonal – is carried by default on a multimode capillary whose radius varies along z, lossless or not. That is a set of two or more plain MarcatiliMode channels forming the contiguous HE_11 to HE_1N block, whose core radius is not constant. Every other set is None: one channel needs no connection; a constant bore – a pressure gradient included, which changes the effective index and not the J_0(u_m r/a) field – has none; and a structured cladding (a Vincetti mode), a tabulated mode or a set of several azimuthal orders is not the block the connection is derived for, and keeps the adiabatic route with its AdiabaticTaperWarning.

Parameters:

modessequence of Mode: The run’s mode set, in the state’s channel order.

Returns:

tuple of MarcatiliMode or None: The block, or None when the coupling does not apply.

Examples

>>> from solprop.modes.capillary import MarcatiliMode, taper_coupling_block
>>> from solprop.profiles import Linear
>>> taper = Linear(125e-6, 62.5e-6, 1e-2)
>>> modes = [MarcatiliMode(taper, "He", 1.0, m=m) for m in (1, 2)]
>>> len(taper_coupling_block(modes))
2
>>> taper_coupling_block([MarcatiliMode(125e-6, "He", 1.0, m=m) for m in (1, 2)])
>>> taper_coupling_block([MarcatiliMode(taper, "He", 1.0, m=m) for m in (1, 3)])

solprop.modes.capillary.vector_wall_operator(modes, omega, *, z=0.0, active=None, scored=None, validity_limit=0.25)[source]#

Build the condensed vector Robin wall generator of one azimuthal order.

In outline (the manual page is docs/guide/waveguides/vector-wall.md), three steps.

The boundary response is a matrix, not a scalar. At grazing incidence the wall imposes \(\partial_rE = ik_oE\) on the field tangential to it and \(\partial_rE = i(k_o/n_r^2)E\) on the field normal to it. In the near-Dirichlet regime the condensation is linear in the impedance, so what enters is \(Z = \operatorname{diag}(1/k_{\rm TM}, 1/k_{\rm TE})\) in the \((\hat r,\hat\theta)\) components. Rotating it into the radial families of a total azimuthal order \(\nu\) – the angular quadrature \(Z_{fg} = \int_0^{2\pi}[Z_{rr}P_fP_g + Z_{\theta\theta}Q_fQ_g]\,d\theta\) over the unit-normalised angular patterns of WallFamily.angular – gives, for \(\nu \ge 1\),

\[Z = \bar z\,I + \delta\,\sigma_x,\qquad \bar z = \tfrac12(1/k_{\rm TM} + 1/k_{\rm TE}),\qquad \delta = \tfrac12(1/k_{\rm TM} - 1/k_{\rm TE}),\]

and, for \(\nu = 0\), the same matrix diagonalised: TM sees \(\bar z + \delta\) and TE sees \(\bar z - \delta\), with no coupling at all. The mean \(\bar z\) is exactly the scalar block’s \(1/k_{\rm eff}\), because v_HE is the mean of v_TE and v_TM; the anisotropy \(\delta\) is what the scalar average discards, and it is (n_r^2 - 1)/(n_r^2 + 1) of the mean – 35 % for silica, not a small correction.

The condensation is rank two. Truncating the Mittag-Leffler expansion of the exact recurrence \(qJ_{p+1}(q) = (p - i\eta)J_p(q)\) family by family gives \(T_p(\lambda) \simeq 2N_f + p_f + h_f^{T}(\lambda - D_f)^{-1}h_f\), and the exact two-family secular equation \(\det[\operatorname{diag}(T_{p_A}, T_{p_B}) - iaZ^{-1}] = 0\) is exactly the secular equation of

\[L_N = -\operatorname{diag}(u_m^2) + VSV^{T},\qquad S = (C - i\,aZ^{-1})^{-1},\qquad C = \operatorname{diag}(2N_f + p_f),\]

with V one boundary trace column per family. The truncation constant carries the family’s Bessel order p_f: without it the condensed eigenvalues stop converging to the exact roots. The paraxial reduction of \(M = k^2I + L_N/a^2\) is \(G = i(M - k^2)/(2k)\) in 1/m, as for the scalar block.

What is new beyond two copies of the scalar block. The scalar block is one family with one averaged impedance, and two copies of it for two polarisations are two copies of the same average. This block is per azimuthal order: it gives TE_0m and TM_0m their own impedances, which differ by n_r^2, and it couples the J_{nu-1} and J_{nu+1} families of every order nu >= 1 through \(\delta\), which converts power between them. Each family in isolation still reduces to this package’s own Marcatili loss for its own v_n.

Parameters:

modessequence of Mode: One total azimuthal order of one capillary, in the state’s channel order: either TM_0m and TE_0m sets, or the hybrid n = +nu set followed by the n = -nu one, each contiguous from m = 1. Either family may be retained alone, in which case the block carries no conversion and the report’s single-family loss limits are what it is. Every channel declares the same guide, fill and polarisation angle; anything else is refused by name.
omeganumpy.ndarray: Angular frequencies in rad/s, shape (n_omega,); the whole grid.
zfloat, optional: Position along the capillary in m the block is frozen at.
activenumpy.ndarray or None, optional: Bool mask of the bins to evaluate, shape (n_omega,). None takes every strictly positive frequency; an out-of-band bin is never evaluated.
scorednumpy.ndarray or None, optional: Bool mask of the bins the validity ratios are taken over, as for wall_operator(); None takes every evaluated bin.
validity_limitfloat, optional: Largest admitted validity ratio; above it the block is refused, above a tenth of it it warns. See CoherentWall.

Returns:

VectorWallOperator: The dense optics-convention generator, its factored pieces and its limits.

Raises:

ModeError: If the set is not one azimuthal order of a bare capillary, if a channel declares a loss setting or the reduced model, if the active band is empty, if the block is not passive on the branch the cladding terms land on, or if a validity ratio exceeds validity_limit.

Warns:

GrazingIncidenceWarning: When a validity ratio exceeds a tenth of validity_limit, naming it.

Examples

>>> import numpy as np
>>> from solprop.modes import MarcatiliMode, vector_wall_operator
>>> modes = [MarcatiliMode(125e-6, "He", 1.0, n=n, m=1) for n in (1, -1)]
>>> block = vector_wall_operator(modes, np.array([2.356e15]))
>>> block.sector, [family.label for family in block.families]
(1, ['HE_1m', 'EH_1m'])
>>> bool(abs(block.generator[0, 1, 0]) > 0.0)
True

solprop.modes.capillary.wall_connection(modes, z=0.0, *, impedance='mean')[source]#

Return the moving-basis connection of a wall block at one position.

The basis of the condensed wall block is attached to the core radius, so a taper moves it. Writing the field as \(E = \sum_n c_n \psi_n(r; a(z))\) with \(\psi_n = \sqrt2 J_p(u_n r/a)/(a|J_{p+1}(u_n)|)\) orthonormal under \(\int_0^a \cdot\, r\,dr\) at every z, \(u_n\) the zeros of \(J_p\), the coefficients pick up \(c' = (iB - C)c\) with \(C_{mn} = \langle\psi_m|\partial_z\psi_n\rangle\):

\[C_{mn} = \frac{a'}{a}\,\frac{2s_ms_nu_mu_n}{u_m^2-u_n^2}\quad(m\ne n), \qquad C_{nn} = 0,\]

with \(s_m = \operatorname{sgn}J_{p+1}(u_m)\) the sign of the basis slope at the wall, \(-J_p'(u_m) = J_{p+1}(u_m)\). For the HE_1m block p = 0 and this is the scalar block’s formula. It is real – so the optics/internal conjugation of docs/dev/conventions.md section 2 leaves it alone – and antisymmetric, because swapping the indices flips the denominator alone, so the connection on its own is norm preserving and moves no power out of the retained set. The sign is checked by differentiating the implemented profile (tests/modes/test_wall_operator.py), not asserted.

The vector block (impedance="te_tm") is one total azimuthal order with at most two radial families, and its connection is block diagonal: each family carries the formula above with its own Bessel order p_f (WallFamily.bessel_order), zeros and signs, and the two families do not connect at all. The reason is exact rather than small: a taper moves the radius and nothing else, and every channel’s field is its radial function of r/a times an angular pattern that does not depend on a, so \(\langle\psi^A_m|\partial_z\psi^B_n\rangle\) is a radial integral times the angular overlap of the two families’ patterns, which is zero – the two hybrid families of an order nu >= 1 overlap as \(\cos 2\nu\theta\), and TE_0m and TM_0m are orthogonal at every point. A family’s implemented sign (WallFamily.radial_sign) enters \(s_m s_n\) squared and cancels. Both statements are measured on the implemented fields, not assumed.

The whole matrix scales with a'(z)/a(z), and nothing in it depends on frequency: the basis is the same set of Bessel orders at every bin.

Parameters:

modessequence of Mode: The block in the state’s channel order: the contiguous HE_11 to HE_1N set for impedance="mean", one azimuthal order of a capillary for "te_tm". Refused by name on exactly the conditions wall_operator() and vector_wall_operator() refuse.
zfloat, optional: Position along the capillary in m.
impedancestr, optional: Which block the connection is of, "mean" or "te_tm", as CoherentWall.impedance.

Returns:

numpy.ndarray: C in 1/m, real, shape (n_modes, n_modes), with a zero diagonal.

Examples

>>> import numpy as np
>>> from solprop.modes.capillary import MarcatiliMode, wall_connection
>>> from solprop.profiles import Linear
>>> taper = Linear(125e-6, 62.5e-6, 1e-2)
>>> modes = [MarcatiliMode(taper, "He", 1.0, m=index) for index in (1, 2)]
>>> link = wall_connection(modes, 5e-3)
>>> bool(np.all(link == -link.T)), float(link[0, 0])
(True, 0.0)
>>> vector = [MarcatiliMode(taper, "He", 1.0, n=n, m=1) for n in (1, -1)]
>>> float(wall_connection(vector, 5e-3, impedance="te_tm")[0, 1])
0.0

solprop.modes.capillary.wall_operator(modes, omega, *, z=0.0, active=None, scored=None, validity_limit=0.25)[source]#

Build the condensed Robin wall generator of a contiguous HE_1m block.

In outline (the manual page is docs/guide/waveguides/coherent-wall.md), the orthonormal radial basis is \(\phi_m(x) = \sqrt2 J_0(u_m x)/|J_1(u_m)|\) on \(x = r/a\), whose derivative at the wall is \(h_m = -\sqrt2\,\operatorname{sgn}J_1(u_m)\,u_m\). Truncating the Mittag-Leffler expansion of the exact outgoing root equation \(qJ_1(q) + i\eta J_0(q) = 0\) at N orders and eliminating the boundary coordinate gives the condensed transverse operator \(L_N = -\operatorname{diag}(u_m^2) + hh^{T}/(2N - i\eta)\), and the paraxial reduction of \(M = k^2I + L_N/a^2\) about the in-core carrier is \(G = i(M - k^2)/(2k) = D + H\) in 1/m.

Two normalisations are re-derived here rather than assumed. The amplitude map is the identity up to one mode-independent positive scale: this package’s internal amplitude is a_m = E_m sqrt(N_m) with N_m = (pi/2) a^2 J_1(u_m)^2 eps_0 c (docs/dev/conventions.md section 4), and both that normalisation and the orthonormal basis above carry the same |J_1(u_m)|, so the two sets of coordinates differ by a sqrt(pi eps_0 c) for every channel alike and the generator is unchanged by the map. The impedance is k_eff = k_0 / v_n, the vacuum wavenumber over the cladding term this package’s own capillary index already uses, because that is the combination whose N = 1 limit is this package’s Marcatili transverse eigenvalue; the carrier k of the paraxial reduction is the in-core wavenumber. Referring the impedance to the core instead is identical for a vacuum core and differs at the order of half the gas susceptibility for a filled one.

Parameters:

modessequence of Mode: The contiguous block HE_11 to HE_1N, in the state’s channel order, sharing one polarisation angle, one capillary and one constant section. Anything else is refused by name.
omeganumpy.ndarray: Angular frequencies in rad/s, shape (n_omega,); the whole grid.
zfloat, optional: Position along the capillary in m the block is frozen at.
activenumpy.ndarray or None, optional: Bool mask of the bins to evaluate, shape (n_omega,). None takes every strictly positive frequency. An empty out-of-band bin is never evaluated: the wall formula divides by omega and its validity ratios are meaningless there, so the generator is left exactly zero.
scorednumpy.ndarray or None, optional: Bool mask of the bins the validity ratios are taken over, shape (n_omega,), intersected with the evaluated bins; a run passes the error norm’s scored band (solprop.grid.TimeGrid.inband_mask), so the absorbed collar a field grid adds below the requested band, at up to twice the longest requested wavelength, neither refuses nor warns. The generator is still built on every evaluated bin. None, or a mask that holds no evaluated bin, takes the ratios over every evaluated bin.
validity_limitfloat, optional: Largest admitted validity ratio; above it the block is refused, above a tenth of it it warns. See CoherentWall.

Returns:

WallOperator: The dense optics-convention generator, its factored pieces and its limits.

Raises:

ModeError: If the set is not one contiguous scalar block, if a channel declares a loss setting or the reduced model, if the active band is empty or holds a non-positive frequency, if the block is not passive on the branch the cladding term lands on, or if a validity ratio exceeds validity_limit.

Warns:

GrazingIncidenceWarning: When a validity ratio exceeds a tenth of validity_limit, naming it.

Examples

>>> import numpy as np
>>> from solprop.modes import MarcatiliMode, wall_operator
>>> modes = [MarcatiliMode(125e-6, "He", 1.0, m=index) for index in (1, 2, 3)]
>>> block = wall_operator(modes, np.array([2.356e15]))
>>> block.n_modes, f"{block.limits.worst:.4f}"
(3, '0.0090')

solprop.modes.capillary.wall_square_root_metric(generator, wavenumber)[source]#

Return the square-root generator and its candidate flux metric, as a diagnostic.

For a constant, real, positive carrier k, the square-root model c' = iBc with B^2 = M = k^2 I - 2ikG satisfies

\[\frac{d}{dz}\left(c^\dagger W c\right) = -c^\dagger\frac{M - M^\dagger}{2ik}c,\qquad W = \frac{B + B^\dagger}{2k},\]

which is the flux balance of this scalar model. Its domain is exactly that statement. Production semantics are the paraxial generator, not this square root; W is positive semidefinite only when the Hermitian part of B is, which positive eigenvalues of B do not establish, and which the principal branch does not establish either; a complex k, an evanescent sector, a varying cross-section – which adds W' – and a varying k all fall outside it. Use it to measure the identity, not to define an error norm or to transport noise.

Parameters:

generatornumpy.ndarray: The optics-convention paraxial generator at one frequency in 1/m, complex, shape (n_modes, n_modes).
wavenumberfloat: The carrier k in 1/m, real and positive.

Returns:

tuple of numpy.ndarray: B in 1/m on the principal branch, and the dimensionless W.

Raises:

ModeError: If the carrier is not real, finite and positive.

Examples

>>> import numpy as np
>>> from solprop.modes import MarcatiliMode, wall_operator
>>> from solprop.modes.capillary import wall_square_root_metric
>>> modes = [MarcatiliMode(125e-6, "He", 1.0, m=index) for index in (1, 2)]
>>> block = wall_operator(modes, np.array([2.356e15]))
>>> k = float(block.wavenumber[0].real)
>>> root, metric = wall_square_root_metric(block.generator[:, :, 0], k)
>>> bool(np.max(abs(root @ root - (k**2 * np.eye(2)
...      - 2j * k * block.generator[:, :, 0]))) < 1e-6 * k**2)
True

solprop.modes.capillary.worst_wall_limits(modes, omega, *, positions, active=None, scored=None, validity_limit=0.25, impedance='mean')[source]#

Return the validity ratios at the position where they are worst.

A guide that varies along z has a different block at every section, and the three ratios of WallLimits are not monotone in the declared parameters in general – the radius, the fill’s index and the cladding term all enter them. The scan is therefore over the positions the caller actually builds the flow from, and the one call that carries the caller’s validity_limit is the one at the worst of them, so a run is refused or warned about once and names the geometry that did it rather than the entrance’s.

Parameters:

modessequence of Mode: The contiguous block HE_11 to HE_1N, in the state’s channel order.
omeganumpy.ndarray: Angular frequencies in rad/s, shape (n_omega,); the whole grid.
positionssequence of float: Positions in m to scan, non-empty.
activenumpy.ndarray or None, optional: Bool mask of the bins to evaluate, shape (n_omega,).
scorednumpy.ndarray or None, optional: Bool mask of the bins the ratios are taken over, shape (n_omega,) (wall_operator()); None takes every evaluated bin.
validity_limitfloat, optional: Largest admitted validity ratio, applied at the worst position alone.
impedancestr, optional: Which block to build, "mean" (wall_operator()) or "te_tm" (vector_wall_operator()). The scan is the same either way.

Returns:

WallLimits: The ratios at that position.

Raises:

ModeError: If positions is empty, or on everything the block’s own constructor refuses.

Warns:

GrazingIncidenceWarning: Once, from the worst position, when a ratio exceeds a tenth of the limit.

Examples

>>> import numpy as np
>>> from solprop.modes.capillary import MarcatiliMode, worst_wall_limits
>>> from solprop.profiles import Linear
>>> taper = Linear(125e-6, 100e-6, 1e-2)
>>> modes = [MarcatiliMode(taper, "He", 1.0, m=index) for index in (1, 2)]
>>> limits = worst_wall_limits(
...     modes, np.array([2.356e15]), positions=[0.0, 5e-3, 1e-2]
... )
>>> f"{limits.radius * 1e6:.1f} um"
'100.0 um'