# SPDX-FileCopyrightText: 2026 John C. Travers and contributors
#
# SPDX-License-Identifier: GPL-3.0-or-later
r"""Design calculators: soliton numbers and dispersive-wave phase matching.
The numbers a person works out before a run. For a pulse of energy :math:`E` and
full width at half maximum :math:`\tau_{\rm fwhm}` in a mode with group-velocity
dispersion :math:`\beta_2`, effective area :math:`A_{\rm eff}` and a fill of nonlinear
index :math:`n_2`,
.. math::
\gamma = \frac{n_2\omega_0}{cA_{\rm eff}},\quad
L_D = \frac{\tau_0^2}{|\beta_2|},\quad
L_{NL} = \frac{1}{\gamma P_0},\quad
N = \sqrt{L_D/L_{NL}},\quad
L_{\rm fiss} = L_D/N,
with the natural width :math:`\tau_0` and the peak power :math:`P_0` of the chosen
pulse shape (:func:`pulse_width`, :func:`pulse_peak_power`). One shape runs through the
whole chain, so the Gaussian numbers are never computed with a ``sech`` width.
A dispersive wave is phase matched to a soliton at the frequencies where
.. math::
\Delta\beta(\omega) = \beta_{\rm dw}(\omega) - \beta_{\rm s}(\omega_0)
- \beta_{1,\rm s}(\omega_0)(\omega-\omega_0) - k_{\rm NL}
vanishes (:func:`phase_mismatch`). The last term is the soliton's nonlinear wavenumber
and uses the largest unchirped constituent when pulse data are given
(:func:`soliton_content`). Without pulse data the condition is linear.
``docs/guide/workflows/design-tools.md`` derives both and says what the nonlinear term
assumes.
Everything is SI except pressure, which is in bar as at the rest of the user interface
(``docs/dev/conventions.md`` section 1).
"""
import logging
import math
import warnings
from collections.abc import Callable, Iterator, Mapping
from contextlib import contextmanager
from dataclasses import dataclass
from itertools import pairwise
from typing import Literal, cast
import numpy as np
from jaxtyping import Float
from scipy.integrate import solve_ivp
from scipy.optimize import brentq
from scipy.special import erfcinv
from solprop import spatial
from solprop._display import Displayed, line
from solprop.constants import EPS0, ROOM_TEMPERATURE, C
from solprop.exceptions import (
DiluteScalingWarning,
MaterialError,
MaterialRangeWarning,
ModeError,
ToolsError,
)
from solprop.modes import base
from solprop.modes.antiresonant import TubeLattice, VincettiMode
from solprop.modes.capillary import LossSpec, MarcatiliMode, MarcatiliModel, ModeKind
from solprop.modes.protocol import Mode
from solprop.modes.quadrature import BASIS_NORMALISATION, mode_matrix, on_axis_rule
from solprop.physdata import Gas, Mixture, chi3_from_n2, n2_from_chi3, species
from solprop.physdata.glasses import GlassIndex, fused_silica
from solprop.physdata.nonlinear_data import carriers_agree
from solprop.physdata.selection import MaterialSelection
from solprop.profiles import ZProfile
__all__ = [
"Fill",
"GasRatio",
"NonlinearPhase",
"PulseShape",
"SolitonContent",
"SolitonParams",
"capillary_design",
"capillary_params",
"dispersion_length",
"field_to_intensity",
"fission_length",
"gas_ratio",
"intensity_to_field",
"lambda_rdw",
"lambda_rdw_all",
"nonlinear_coefficient",
"nonlinear_length",
"phase_mismatch",
"pressure_rdw",
"pressure_zdw",
"pulse_peak_power",
"pulse_width",
"self_focusing_power",
"soliton_content",
"soliton_order",
"soliton_params",
]
logger = logging.getLogger(__name__)
PulseShape = Literal["sech", "gauss"]
"""The two pulse shapes the design numbers are defined for."""
NonlinearPhase = Literal["largest", "supplied_power", "linear"]
"""Choose the constituent, supplied-power or linear RDW phase model."""
Fill = Gas | str | Mapping[str, float] | Mixture
"""A fill as the run interfaces accept it: a species, mole fractions or a mixture."""
_SECH_WIDTH = 2.0 * math.log(1.0 + math.sqrt(2.0))
"""``fwhm / tau0`` of ``sech^2(t / tau0)``: twice the inverse cosh of the root of 2."""
_GAUSS_WIDTH = 2.0 * math.sqrt(math.log(2.0))
"""``fwhm / tau0`` of ``exp(-t^2 / tau0^2)``."""
_MISMATCH_SAMPLES = 2048
"""Frequency samples the phase-matching search brackets its roots on.
The mismatch of a hollow waveguide is smooth on the scale of the search window (its
structure is set by the zero-dispersion points, of which there are a handful), so two
roots closer than a two-thousandth of the window are not separated. The cost is one
vectorised index evaluation.
"""
_PRESSURE_FLOOR = 1e-6
"""Lower end of the pressure brackets in bar; a true vacuum has no gas term."""
_ZDW_RANGE = (100e-9, 3e-6)
"""Window of the zero-dispersion search in m: :func:`modes.base.zdw`'s."""
_ZDW_SAMPLES = 512
"""Samples, evenly spaced in frequency, on which the dispersion's sign is read."""
_POLE_FACTOR = 2.0
"""How far the group delay at a refined sign change of ``beta2`` may sit from its value
at the bracket's ends, in units of the change a bounded ``beta2`` could produce across
the bracket, before the sign change is read as a pole rather than a zero. At a zero
``beta1`` is extremal and moves by at most half that unit; at a gas resonance it
diverges. The differenced ``beta2`` itself cannot tell them apart: ``beta1`` is even
about a pole, so its central difference vanishes there too."""
@contextmanager
def _probing() -> Iterator[None]:
"""Silence material-range and dilute-scaling warnings while a search samples.
A root search evaluates the index across a window far wider than its answer --
in frequency, or in pressure up to a hundred bar, where an optmat fill's linear
density scaling departs from its finite-density law -- and the ends are probes,
not results. The answer itself is re-evaluated outside this context, so a root
that lies outside a model's range, or at a density beyond the dilute law,
still warns.
Yields
------
None
Control, with the warning categories ignored.
"""
with warnings.catch_warnings():
warnings.simplefilter("ignore", MaterialRangeWarning)
warnings.simplefilter("ignore", DiluteScalingWarning)
yield
def _positive(name: str, value: float) -> float:
"""Return a finite positive number or refuse it by name.
Parameters
----------
name : str
The argument's name, for the message.
value : float
The value.
Returns
-------
float
The value as a float.
Raises
------
ToolsError
If it is not finite and positive.
"""
number = float(value)
if not math.isfinite(number) or number <= 0.0:
raise ToolsError(f"{name} must be finite and positive, got {value!r}")
return number
def _omega(wavelength: float) -> float:
"""Return the angular frequency of a vacuum wavelength.
Parameters
----------
wavelength : float
Vacuum wavelength in m, > 0.
Returns
-------
float
Angular frequency in rad/s.
"""
return 2.0 * math.pi * C / _positive("wavelength", wavelength)
def _check_shape(shape: str) -> None:
"""Refuse a pulse shape the design numbers are not defined for.
Parameters
----------
shape : str
The shape keyword.
Raises
------
ToolsError
If it is neither ``"sech"`` nor ``"gauss"``.
"""
if shape not in ("sech", "gauss"):
raise ToolsError(f"shape must be 'sech' or 'gauss', got {shape!r}")
[docs]
def pulse_width(fwhm: float, shape: PulseShape = "sech") -> float:
r"""Return the natural width of a pulse from its full width at half maximum.
Parameters
----------
fwhm : float
Full width at half maximum of the **power** in s, > 0.
shape : {'sech', 'gauss'}, optional
``"sech"`` for :math:`P\propto{\rm sech}^2(t/\tau_0)`, ``"gauss"`` for
:math:`P\propto\exp(-t^2/\tau_0^2)`.
Returns
-------
float
:math:`\tau_0` in s: ``fwhm / (2 ln(1 + sqrt 2))`` or ``fwhm / (2 sqrt(ln 2))``.
Raises
------
ToolsError
If the shape is unknown or the width is not positive.
Examples
--------
>>> from solprop.tools import pulse_width
>>> f"{pulse_width(10e-15) * 1e15:.4f}"
'5.6730'
>>> f"{pulse_width(10e-15, 'gauss') * 1e15:.4f}"
'6.0056'
"""
_check_shape(shape)
return _positive("fwhm", fwhm) / (_SECH_WIDTH if shape == "sech" else _GAUSS_WIDTH)
[docs]
def pulse_peak_power(energy: float, fwhm: float, shape: PulseShape = "sech") -> float:
r"""Return the peak power of a pulse of given energy, duration and shape.
Parameters
----------
energy : float
Pulse energy in J, > 0.
fwhm : float
Full width at half maximum of the power in s, > 0.
shape : {'sech', 'gauss'}, optional
The pulse shape; see :func:`pulse_width`.
Returns
-------
float
Peak power in W: :math:`E/(2\tau_0)` for ``sech`` and
:math:`E/(\sqrt\pi\,\tau_0)` for a Gaussian, the integrals of the two shapes.
Raises
------
ToolsError
If the shape is unknown or a quantity is not positive.
Examples
--------
>>> from solprop.tools import pulse_peak_power
>>> f"{pulse_peak_power(1e-3, 10e-15) / 1e9:.3f}"
'88.137'
>>> f"{pulse_peak_power(1e-3, 10e-15, 'gauss') / 1e9:.3f}"
'93.944'
"""
width = pulse_width(fwhm, shape)
integral = 2.0 if shape == "sech" else math.sqrt(math.pi)
return _positive("energy", energy) / (integral * width)
[docs]
def dispersion_length(fwhm: float, beta2: float, shape: PulseShape = "sech") -> float:
r"""Return the dispersion length :math:`L_D = \tau_0^2/|\beta_2|`.
Parameters
----------
fwhm : float
Full width at half maximum of the power in s, > 0.
beta2 : float
Group-velocity dispersion in s^2/m, non-zero.
shape : {'sech', 'gauss'}, optional
The pulse shape; see :func:`pulse_width`.
Returns
-------
float
Dispersion length in m.
Raises
------
ToolsError
If the dispersion is zero or not finite, or the width is not positive.
Examples
--------
>>> from solprop.tools import dispersion_length
>>> f"{dispersion_length(10e-15, -8e-30):.3f}"
'4.023'
"""
if not math.isfinite(beta2) or beta2 == 0.0:
raise ToolsError(f"beta2 must be finite and non-zero, got {beta2!r}")
return pulse_width(fwhm, shape) ** 2 / abs(beta2)
[docs]
def nonlinear_length(peak_power: float, gamma: float) -> float:
r"""Return the nonlinear length :math:`L_{NL} = 1/(\gamma P_0)`.
Parameters
----------
peak_power : float
Peak power in W, > 0.
gamma : float
Nonlinear coefficient in rad/(W m), > 0.
Returns
-------
float
Nonlinear length in m.
Raises
------
ToolsError
If either quantity is not positive.
Examples
--------
>>> from solprop.tools import nonlinear_length
>>> nonlinear_length(1e9, 1e-9)
1.0
"""
return 1.0 / (_positive("gamma", gamma) * _positive("peak_power", peak_power))
[docs]
def soliton_order(
peak_power: float,
fwhm: float,
gamma: float,
beta2: float,
shape: PulseShape = "sech",
) -> float:
r"""Return the soliton order :math:`N = \sqrt{L_D/L_{NL}}`.
The number is defined for either sign of the dispersion; it is a
soliton order in the physical sense only where ``beta2 < 0``.
Parameters
----------
peak_power : float
Peak power in W, > 0.
fwhm : float
Full width at half maximum of the power in s, > 0.
gamma : float
Nonlinear coefficient in rad/(W m), > 0.
beta2 : float
Group-velocity dispersion in s^2/m, non-zero.
shape : {'sech', 'gauss'}, optional
The pulse shape, used for the width **throughout**.
Returns
-------
float
The soliton order, dimensionless.
Examples
--------
>>> from solprop.tools import soliton_order
>>> f"{soliton_order(1e9, 10e-15, 1e-9, -8e-30):.4f}"
'2.0057'
"""
return math.sqrt(
dispersion_length(fwhm, beta2, shape) / nonlinear_length(peak_power, gamma)
)
[docs]
@dataclass(frozen=True)
class SolitonContent:
"""Describe the discrete soliton content of an unchirped single-lobe pulse.
Attributes
----------
area_order : float
Dimensionless field area divided by pi, not the conventional soliton order.
count : int
Number of strictly positive discrete eigenvalues; zero at first threshold.
largest_eigenvalue : float or None
Imaginary part of the largest Zakharov-Shabat eigenvalue, dimensionless
with time scaled by the input shape's natural width. None if count is zero.
largest_peak_power : float or None
Peak power of the largest constituent in W, or None if count is zero.
largest_width : float or None
Natural sech width of that constituent in s, or None if count is zero.
largest_nonlinear_wavenumber : float or None
Positive optics-convention phase rate in 1/m, or None if count is zero.
"""
area_order: float
count: int
largest_eigenvalue: float | None
largest_peak_power: float | None
largest_width: float | None
largest_nonlinear_wavenumber: float | None
def _scattering_angle(
x: float, theta: Float[np.ndarray, " n_components"], order: float, eta: float
) -> list[float]:
"""Evaluate the real Zakharov-Shabat phase equation for a Gaussian.
Parameters
----------
x : float
Time divided by the input width.
theta : numpy.ndarray
Unwrapped scattering angle, one element.
order : float
Conventional soliton order, the field's peak.
eta : float
Trial imaginary eigenvalue.
Returns
-------
list of float
Angle derivative.
"""
return [-order * math.exp(-x * x / 2.0) - eta * math.sin(2 * theta[0])]
def _gaussian_eigenvalue(order: float) -> float:
"""Solve the half-line scattering problem for the largest Gaussian eigenvalue.
Parameters
----------
order : float
Conventional soliton order, above the first Gaussian threshold.
Returns
-------
float
Positive imaginary part of the dimensionless eigenvalue.
Raises
------
ToolsError
If integration or eigenvalue refinement fails.
"""
# The omitted field area is <= 1e-13. Moving the boundary with the amplitude
# avoids assuming that the Gaussian tail is negligible for every input N.
extent = math.sqrt(2.0) * float(erfcinv(1e-13 / (order * math.sqrt(math.pi / 2.0))))
def residual(eta: float) -> float:
"""Return the ground-state symmetry residual at the pulse centre.
Parameters
----------
eta : float
Trial imaginary eigenvalue.
Returns
-------
float
Phase angle plus pi/4.
"""
if eta == 0.0:
# The exact area protects the first threshold from tail truncation.
return math.pi / 4.0 - order * math.sqrt(math.pi / 2.0)
solution = solve_ivp(
_scattering_angle,
(-extent, 0.0),
[0.0],
method="DOP853",
args=(order, eta),
rtol=1e-11,
atol=1e-13,
)
if not solution.success:
raise ToolsError(
f"Gaussian scattering integration failed: {solution.message}"
)
return float(solution.y[0, -1]) + math.pi / 4.0
# The ground-state angle crosses -pi/4 once for 0 < eta < max(q).
# Unlike shooting the amplitudes, this does not overflow in the tails.
try:
return float(brentq(residual, 0.0, order, xtol=1e-13, rtol=1e-11))
except (ValueError, RuntimeError) as error:
raise ToolsError(f"Gaussian eigenvalue refinement failed: {error}") from error
[docs]
def soliton_content(
peak_power: float,
fwhm: float,
gamma: float,
beta2: float,
shape: PulseShape = "sech",
) -> SolitonContent:
r"""Calculate the soliton count and largest constituent of an unchirped pulse.
This is the scattering content of the lossless cubic focusing NLSE, not a
prediction of fission dynamics with higher-order dispersion, Raman or plasma.
The conventional :func:`soliton_order` remains a separate design parameter.
Parameters
----------
peak_power : float
Input peak power in W, finite and positive.
fwhm : float
Input power FWHM in s, finite and positive.
gamma : float
Nonlinear coefficient in rad/(W m), finite and positive.
beta2 : float
Anomalous group-velocity dispersion in s^2/m, finite and negative.
shape : {'sech', 'gauss'}, optional
Unchirped input field shape; the Gaussian power is exp(-t^2/tau0^2).
Returns
-------
SolitonContent
Analytic count and largest constituent. The sech eigenvalue is analytic;
the Gaussian eigenvalue is numerical, with dimensionless absolute accuracy
about 1e-10 away from thresholds. Relative accuracy degrades at threshold.
Raises
------
ToolsError
If inputs are invalid, dispersion is not anomalous, or scattering fails.
References
----------
.. [1] Klaus and Shaw, SIAM J. Math. Anal. 34, 759 (2003),
doi:10.1137/S0036141002403067.
.. [2] Kim, Lee and Lyng, arXiv:1310.4145, sections 1-2.
"""
if not math.isfinite(beta2) or beta2 >= 0.0:
raise ToolsError(
"soliton content requires finite anomalous dispersion (beta2 < 0)"
)
order = soliton_order(peak_power, fwhm, gamma, beta2, shape)
area_order = order if shape == "sech" else math.sqrt(2.0 / math.pi) * order
if not math.isfinite(area_order):
raise ToolsError("the normalised pulse area must be finite")
# Strictly positive eigenvalues only: at a half-integer area the new state is
# at zero and is not a square-integrable soliton (Klaus-Shaw threshold).
count = max(0, math.ceil(area_order - 0.5))
if count == 0:
return SolitonContent(area_order, 0, None, None, None, None)
eta = order - 0.5 if shape == "sech" else _gaussian_eigenvalue(order)
if eta <= 0.0:
raise ToolsError(
"the positive Gaussian eigenvalue is too close to threshold to resolve"
)
power = peak_power * (2.0 * eta / order) ** 2
width = pulse_width(fwhm, shape) / (2.0 * eta)
return SolitonContent(area_order, count, eta, power, width, 0.5 * gamma * power)
[docs]
def fission_length(
peak_power: float,
fwhm: float,
gamma: float,
beta2: float,
shape: PulseShape = "sech",
) -> float:
r"""Return the soliton fission length :math:`L_{\rm fiss} = L_D/N`.
Parameters
----------
peak_power : float
Peak power in W, > 0.
fwhm : float
Full width at half maximum of the power in s, > 0.
gamma : float
Nonlinear coefficient in rad/(W m), > 0.
beta2 : float
Group-velocity dispersion in s^2/m, non-zero.
shape : {'sech', 'gauss'}, optional
The pulse shape; see :func:`pulse_width`.
Returns
-------
float
Fission length in m, the geometric mean of the two lengths.
Examples
--------
>>> from solprop.tools import fission_length
>>> f"{fission_length(1e9, 10e-15, 1e-9, -8e-30):.4f}"
'2.0057'
"""
order = soliton_order(peak_power, fwhm, gamma, beta2, shape)
return dispersion_length(fwhm, beta2, shape) / order
[docs]
def nonlinear_coefficient(
mode: Mode, n2: float, wavelength: float, z: float = 0.0
) -> float:
r"""Return the nonlinear coefficient :math:`\gamma = n_2\omega/(cA_{\rm eff})`.
Parameters
----------
mode : Mode
The mode, for its effective area.
n2 : float
Nonlinear index of the fill in m^2/W, >= 0.
wavelength : float
Vacuum wavelength in m, > 0.
z : float, optional
Position along the waveguide, m.
Returns
-------
float
Nonlinear coefficient in rad/(W m).
Raises
------
ToolsError
If the nonlinear index is negative or not finite.
Examples
--------
>>> from solprop.modes import MarcatiliMode
>>> from solprop.tools import nonlinear_coefficient
>>> mode = MarcatiliMode(125e-6, "He", 1.0)
>>> f"{nonlinear_coefficient(mode, 3.55e-25, 800e-9):.4e}"
'1.1919e-10'
"""
if not math.isfinite(n2) or n2 < 0.0:
raise ToolsError(f"n2 must be finite and not negative, got {n2!r}")
return n2 * _omega(wavelength) / (C * base.aeff(mode, z))
[docs]
def self_focusing_power(wavelength: float, index: float, n2: float) -> float:
r"""Return the critical power for self-focusing, as physical power.
:math:`P_{\rm cr} = 1.86226\,\lambda^2/(4\pi n_0 n_2)`, the Townes power of Fibich
and Gaeta (Opt. Lett. 25, 335 (2000)). It is
:func:`solprop.spatial.critical_power` -- the package's one Townes constant --
at the linear index, both measuring the physical axial power in the medium
(SEP 34). The paper prints the constant as 1.86225; the
computed value is 1.8622555, of which 1.86226 is the rounding.
Parameters
----------
wavelength : float
Vacuum wavelength in m, > 0.
index : float
Linear refractive index :math:`n_0`, > 0.
n2 : float
Nonlinear index in m^2/W, > 0.
Returns
-------
float
Critical power in W.
Examples
--------
>>> from solprop.tools import self_focusing_power
>>> f"{self_focusing_power(800e-9, 1.0, 3.55e-25) / 1e9:.2f}"
'267.17'
"""
linear = _positive("index", index)
chi3 = chi3_from_n2(_positive("n2", n2), linear)
return spatial.critical_power(
_positive("wavelength", wavelength), chi3, index=linear
)
[docs]
def field_to_intensity(field: float) -> float:
r"""Return the intensity of a field amplitude, :math:`I = \epsilon_0 c E^2/2`.
Parameters
----------
field : float
Field amplitude in V/m.
Returns
-------
float
Intensity in W/m^2.
Examples
--------
>>> from solprop.tools import field_to_intensity
>>> f"{field_to_intensity(1e10):.4e}"
'1.3272e+17'
"""
return 0.5 * EPS0 * C * float(field) ** 2
[docs]
def intensity_to_field(intensity: float) -> float:
r"""Return the field amplitude of an intensity, :math:`E = \sqrt{2I/\epsilon_0 c}`.
Parameters
----------
intensity : float
Intensity in W/m^2, >= 0.
Returns
-------
float
Field amplitude in V/m.
Raises
------
ToolsError
If the intensity is negative.
Examples
--------
>>> from solprop.tools import field_to_intensity, intensity_to_field
>>> f"{intensity_to_field(field_to_intensity(1e10)):.6e}"
'1.000000e+10'
"""
if not intensity >= 0.0:
raise ToolsError(f"an intensity must not be negative, got {intensity!r}")
return math.sqrt(2.0 * float(intensity) / (EPS0 * C))
# ---------------------------------------------------------------- soliton parameters
[docs]
@dataclass(frozen=True)
class GasRatio:
"""Ratios of the bulk dispersion and the Kerr response of two gases.
Attributes
----------
beta2 : float
Ratio of the bulk group-velocity dispersions, first gas over second.
chi3 : float
Ratio of the third-order susceptibilities, first gas over second.
"""
beta2: float
chi3: float
[docs]
@dataclass(frozen=True, eq=False)
class SolitonParams(Displayed):
"""The design numbers of one pulse in one mode and one fill.
Attributes
----------
energy : float
Pulse energy in J.
fwhm : float
Full width at half maximum of the power in s.
width : float
Natural width ``tau0`` of the shape in s (:func:`pulse_width`).
omega : float
Angular frequency of the pump in rad/s.
wavelength : float
Vacuum wavelength of the pump in m.
gas : str
The fill's species labels, joined by ``+``.
pressure : float
Total pressure at the position in bar; NaN for a mixture of independent
populations, which has no single pressure.
temperature : float
Fill temperature in K.
shape : {'sech', 'gauss'}
``"sech"`` or ``"gauss"``.
peak_power : float
Peak power in W.
beta2 : float
Group-velocity dispersion of the mode at the pump in s^2/m.
density : float
Total number density of the fill in m^-3.
index : float
Linear refractive index of the fill at the pump.
n2 : float
Nonlinear index of the fill in m^2/W.
gamma : float
Nonlinear coefficient in rad/(W m).
soliton_order : float
``sqrt(L_D / L_NL)``.
soliton_content : SolitonContent or None
Unchirped NLSE constituent count and largest soliton; None for normal GVD.
intensity : float
Mode-averaged intensity ``peak_power / effective_area`` in W/m^2.
peak_intensity : float
On-axis intensity of the mode at the peak power in W/m^2: what a run's
``peak_intensity`` statistic reports at its first plane. For the fundamental
of a capillary it is 1.77 times `intensity`. NaN for a mode without a
transverse profile (a dispersion-specified one).
critical_power : float
Self-focusing power in W (:func:`self_focusing_power`).
dispersion_length : float
In m.
nonlinear_length : float
In m.
fission_length : float
In m.
zdw : float
Zero-dispersion wavelength in m within 100 nm to 3 um, or NaN when the
dispersion does not change sign there. No second, guess-seeded search is
tried.
loss_length : float
The 1/e **power** loss length of the mode at the pump in m.
effective_area : float
In m^2.
mode : Mode
The mode.
"""
energy: float
fwhm: float
width: float
omega: float
wavelength: float
gas: str
pressure: float
temperature: float
shape: PulseShape
peak_power: float
beta2: float
density: float
index: float
n2: float
gamma: float
soliton_order: float
soliton_content: SolitonContent | None
intensity: float
peak_intensity: float
critical_power: float
dispersion_length: float
nonlinear_length: float
fission_length: float
zdw: float
loss_length: float
effective_area: float
mode: Mode
def __repr__(self) -> str:
"""Return a one-line description: the pulse, the fill and the soliton order.
Returns
-------
str
For example ``SolitonParams(300 uJ, 10 fs at 800 nm, He 0.4 bar, N=2.23)``.
"""
return line(
"SolitonParams",
[
("", f"{self.energy * 1e6:.4g} uJ"),
("", f"{self.fwhm * 1e15:.4g} fs at {self.wavelength * 1e9:.4g} nm"),
("", f"{self.gas} {self.pressure:.4g} bar"),
("N", f"{self.soliton_order:.3g}"),
],
)
def __str__(self) -> str:
"""Return the parameter summary, section by section.
Returns
-------
str
MODE, WAVEGUIDE, FILL, PULSE, DISPERSION, POWER, INTENSITY and SOLITON,
in the customary units: ``n2`` in cm^2/W, intensities in W/cm^2.
"""
zdw = (
"none in 100-3000 nm"
if math.isnan(self.zdw)
else f"{self.zdw * 1e9:.1f} nm"
)
ratio = self.peak_power / self.critical_power
sections = (
("MODE", repr(self.mode)),
(
"WAVEGUIDE",
f"Aeff = {self.effective_area:.2e} m^2, "
f"Lloss = {self.loss_length:.2e} m (1/e power)",
),
(
"FILL",
f"{self.pressure:.3g} bar {self.gas} at {self.temperature:g} K, "
f"Pcr = {self.critical_power:.2e} W,\n "
f"gamma = {self.gamma:.2e} /(W m), n2 = {self.n2 * 1e4:.2e} cm^2/W",
),
(
"PULSE",
f"{self.energy:.3e} J, {self.fwhm:.3e} s @ "
f"{self.wavelength * 1e9:.1f} nm (shape: {self.shape})",
),
(
"DISPERSION",
f"{self.beta2:.3e} s^2/m @ {self.wavelength * 1e9:.1f} nm, ZDW = {zdw}",
),
("POWER", f"{self.peak_power:.2e} W ({ratio:.4f} of Pcr)"),
(
"INTENSITY",
f"{self.intensity * 1e-4:.2e} W/cm^2 mode-averaged, "
f"{self.peak_intensity * 1e-4:.2e} W/cm^2 on axis",
),
(
"SOLITON",
f"Ld = {self.dispersion_length:.2e} m, Lnl = "
f"{self.nonlinear_length:.2e} m, Lfiss = {self.fission_length:.2e} m, "
f"N = {self.soliton_order:.2f}",
),
)
content = self.soliton_content
if content is None:
constituents = "unavailable (normal dispersion)"
elif content.count == 0:
constituents = "0 discrete solitons"
else:
constituents = (
f"{content.count} discrete solitons, "
f"area/pi = {content.area_order:.3g}\n "
f"largest P = {content.largest_peak_power:.3e} W, "
f"T = {content.largest_width:.3e} s, "
f"kNL = {content.largest_nonlinear_wavenumber:.3e} 1/m"
)
return "\n".join(f"{name}:\n {text}" for name, text in sections) + (
f"\nCONSTITUENTS (unchirped NLSE):\n {constituents}"
)
def _resolve_fill(
gas: Fill,
pressure: float | None,
temperature: float | None,
materials: MaterialSelection | None = None,
) -> Mixture:
"""Return the mixture a fill declaration means, as the run interfaces read it.
Parameters
----------
gas : Gas, str, Mapping or Mixture
A built-in gas, ``"air"`` (nitrogen and oxygen), mole fractions, or a
complete mixture.
pressure : float or None
Total pressure in bar; must be ``None`` beside a `Mixture`, which carries its
own, and defaults to one bar otherwise.
temperature : float or None
Fill temperature in K; ``None`` is room temperature, or the mixture's own
beside a `Mixture`, which refuses any other.
materials : MaterialSelection or None, optional
Resolve a name, a fraction mapping's names and ``"air"`` from optmat under
this selection, as the drivers' ``capillary_fill`` does; ``None`` is
``MaterialSelection()``. A mixture is taken as it is.
Returns
-------
Mixture
The ordered populations.
Raises
------
ToolsError
If a pressure, or a temperature other than its own, is given beside a
mixture, or the mixture's temperature varies along z.
"""
if isinstance(gas, Mixture):
if pressure is not None:
raise ToolsError("a Mixture carries its own pressure; omit pressure")
if isinstance(gas.temperature, ZProfile):
raise ToolsError(
"the design calculators take one temperature in K, and this mixture's "
"varies along z; give a Mixture at one temperature"
)
if temperature is not None and float(temperature) != gas.temperature:
# As the run interfaces: the mixture's densities were set at its own
# temperature, so another one cannot be honoured.
raise ToolsError(
f"the Mixture carries its own temperature ({gas.temperature:g} K), "
f"and temperature={temperature!r} was given beside it: omit "
"temperature=, or build the Mixture at the temperature you mean"
)
return gas
temperature = ROOM_TEMPERATURE if temperature is None else temperature
total = 1.0 if pressure is None else float(pressure)
if isinstance(gas, Mapping):
fractions = {
species(name, selection=materials): fraction
for name, fraction in gas.items()
}
return Mixture.fill(fractions, total, temperature=temperature)
if gas is Gas.Air or (isinstance(gas, str) and gas.strip().lower() == "air"):
return Mixture.air(
total, temperature=temperature, argon=False, selection=materials
)
return Mixture.fill(
{species(gas, selection=materials): 1.0}, total, temperature=temperature
)
def _check_carrier(materials: MaterialSelection | None, wavelength: float) -> None:
"""Refuse a selection whose carrier is not the tool's pump wavelength.
A run under a selection refuses a carrier other than its grid's reference
wavelength, so a design tool refuses one other than the wavelength it evaluates
the Kerr coefficient at: the two then evaluate it at the same carrier.
Parameters
----------
materials : MaterialSelection or None
The design's material selection.
wavelength : float
The pump wavelength in m.
Raises
------
ToolsError
If the selection names another carrier.
"""
carrier = None if materials is None else materials.carrier_wavelength
if carrier is not None and not carriers_agree(float(carrier), wavelength):
msg = (
f"the selection's carrier_wavelength is {carrier:g} m and the design's "
f"wavelength {wavelength:g} m: one carrier sets the Kerr coefficient of "
"every optmat species, as in a run. Leave the selection's "
"carrier_wavelength unset, or design at it"
)
raise ToolsError(msg)
def _default_cladding(
cladding: GlassIndex | complex | None,
) -> GlassIndex | complex:
"""Return the cladding a design's capillary takes, as the capillary driver does.
Parameters
----------
cladding : GlassIndex, complex or None
The cladding given; ``None`` is the package default.
Returns
-------
GlassIndex or complex
A given cladding as it is; the default as optmat's recommended fused
silica.
"""
return fused_silica() if cladding is None else cladding
def _fill_chi3(mixture: Mixture, wavelength: float, z: float = 0.0) -> float:
"""Return a fill's density-weighted third-order susceptibility.
Parameters
----------
mixture : Mixture
The fill.
wavelength : float
Vacuum wavelength in m, the carrier at which an optmat Kerr model is taken.
z : float, optional
Position in m at which the densities are evaluated.
Returns
-------
float
``chi3`` in m^2/V^2.
Raises
------
ToolsError
If a component of the fill has no Kerr data.
"""
try:
return float(mixture.gamma3_path(wavelength)(np.array([float(z)]))[0])
except MaterialError as error:
raise ToolsError(f"the fill has no nonlinear index: {error}") from error
def _fill_n2(
mixture: Mixture, wavelength: float, z: float = 0.0
) -> tuple[float, float]:
"""Return a fill's nonlinear index and its linear index at a wavelength.
Parameters
----------
mixture : Mixture
The fill.
wavelength : float
Vacuum wavelength in m, the carrier at which an optmat Kerr model is taken.
z : float, optional
Position in m at which the densities are evaluated.
Returns
-------
tuple of float
``n2`` in m^2/W, from the density-weighted ``chi3`` and the exact ``n0^2`` of
the fill, and that linear index ``n0``.
Raises
------
ToolsError
If a component of the fill has no Kerr data.
"""
chi3 = _fill_chi3(mixture, wavelength, z)
here = np.array([float(z)])
susceptibility = mixture.susceptibility(np.array([_omega(wavelength)]), here)[0, 0]
index = float(np.real(np.sqrt(1.0 + susceptibility)))
return n2_from_chi3(chi3, index), index
def _on_axis_intensity(mode: Mode, peak_power: float, z: float) -> float:
"""Return the on-axis intensity of a mode carrying a given power.
Parameters
----------
mode : Mode
The mode.
peak_power : float
Power in W.
z : float
Position along the waveguide, m.
Returns
-------
float
Intensity on the axis in W/m^2, both transverse components included; NaN for
a mode that declares no transverse profile.
"""
# The same synthesis the run's peak_intensity statistic uses: the normalised mode
# field at the origin in V/m per sqrt(W), so that eps0 c |e|^2 / 2 is the
# intensity per watt.
try:
on_axis = mode_matrix(
(mode,), on_axis_rule(mode.dimlimits(z)[0]), z, components=("x", "y")
)
except ModeError:
# A dispersion-specified mode has an area and no transverse profile, so it has
# a mode-averaged intensity and no on-axis one.
return math.nan
per_watt = BASIS_NORMALISATION * float(np.sum(np.abs(on_axis.synthesis) ** 2))
return per_watt * peak_power
def _zdw_near(mode: Mode, omega: float, z: float) -> float:
"""Return the zero-dispersion wavelength nearest a pump frequency.
:func:`solprop.modes.base.zdw` brackets one sign change across its whole window,
which a gas resonance in the vacuum ultraviolet defeats: krypton flips the sign of
the dispersion twice near 120 nm, the ends then agree, and the zero at 769 nm of a
14 um capillary at 15 bar is refused. Here the sign is sampled, every change is
refined, the changes through a pole are dropped, and the zero nearest the pump is
the one reported. Two limits of a sampled search remain: a zero and a pole inside
one sample interval (3.6e13 rad/s) restore the sign and are both missed, and a
sample that is not finite cannot bracket anything.
Parameters
----------
mode : Mode
The mode.
omega : float
Pump angular frequency in rad/s.
z : float
Position along the waveguide, m.
Returns
-------
float
Wavelength in m, or NaN when the dispersion has no zero in 100 nm to 3 um.
"""
grid = np.linspace(
2.0 * math.pi * C / _ZDW_RANGE[1],
2.0 * math.pi * C / _ZDW_RANGE[0],
_ZDW_SAMPLES,
)
zeros = []
# Brent's method converges onto a pole as readily as onto a zero, and a model may
# divide by zero exactly there; the group-delay test below tells the two apart.
with _probing(), np.errstate(divide="ignore", invalid="ignore"):
sampled = base.beta2(mode, grid, z)
for i in np.flatnonzero(np.sign(sampled[:-1]) * np.sign(sampled[1:]) < 0.0):
low, high = float(grid[i]), float(grid[i + 1])
root = float(
brentq(
lambda value: float(base.beta2(mode, np.array([value]), z)[0]),
low,
high,
)
)
delay = base.beta1(mode, np.array([low, root, high]), z)
bounded = max(abs(sampled[i]), abs(sampled[i + 1])) * (high - low)
if abs(delay[1] - delay[0]) <= _POLE_FACTOR * bounded:
zeros.append(root)
if not zeros:
return math.nan
nearest = min(zeros, key=lambda value: abs(value - omega))
# Outside the probing context, so a zero beyond the index model's range warns.
base.beta(mode, np.array([nearest]), z)
return 2.0 * math.pi * C / nearest
[docs]
def soliton_params( # noqa: PLR0913 - a pulse (three numbers and a shape) and a fill
mode: Mode,
gas: Fill,
pressure: float | None = None,
*,
energy: float,
fwhm: float,
wavelength: float,
shape: PulseShape = "sech",
temperature: float | None = None,
z: float = 0.0,
materials: MaterialSelection | None = None,
) -> SolitonParams:
"""Return the design numbers of a pulse in any mode.
The mode supplies the dispersion, the area and the loss; the fill supplies the
nonlinear index. They are separate arguments because a mode need not know its
fill's Kerr response (a tabulated mode does not), and the caller is trusted to
give the fill the mode was built with -- :func:`capillary_params` does both.
Parameters
----------
mode : Mode
Any mode: a capillary, an antiresonant or tabulated fibre, a rectangular
guide.
gas : Gas, str, Mapping or Mixture
The fill: a registry species, ``"air"``, mole fractions, or a mixture. A
mixture of species selected from optmat takes their Kerr coefficients at
`wavelength`, the carrier.
pressure : float or None, optional
Total pressure in **bar**; ``None`` is one bar, and is required beside a
`Mixture`.
energy : float
Pulse energy in J, > 0.
fwhm : float
Full width at half maximum of the power in s, > 0.
wavelength : float
Pump wavelength in m, > 0.
shape : {'sech', 'gauss'}, optional
The pulse shape, used consistently for the width and the peak power.
temperature : float or None, optional
Fill temperature in K; ``None`` is room temperature, or the mixture's own
beside a `Mixture`, which refuses any other.
z : float, optional
Position along the waveguide at which everything is evaluated, m.
materials : MaterialSelection or None, optional
Resolve a named fill from optmat under this selection, as a run's
``materials=`` does: its index, density and Kerr coefficient at
`wavelength`. ``None`` is ``MaterialSelection()``; a `Mixture` is taken as
it is. The mode is the caller's either way.
Returns
-------
SolitonParams
The numbers, with the mode.
Raises
------
ToolsError
If a quantity is not positive, the shape is unknown, the mode has no
dispersion at the pump, a component of the fill has no Kerr data, the
selection names a carrier other than `wavelength`, or a pressure or a
temperature other than its own is given beside a `Mixture`.
Examples
--------
>>> from solprop.modes import MarcatiliMode
>>> from solprop.tools import soliton_params
>>> mode = MarcatiliMode(125e-6, "He", 0.4, model="reduced", cladding=1.45)
>>> p = soliton_params(mode, "He", 0.4, energy=300e-6, fwhm=10e-15,
... wavelength=800e-9)
>>> f"{p.soliton_order:.2f} {p.fission_length:.2f}"
'2.26 1.75'
"""
materials = MaterialSelection() if materials is None else materials
_check_carrier(materials, wavelength)
mixture = _resolve_fill(gas, pressure, temperature, materials)
omega = _omega(wavelength)
here = np.array([float(z)])
n2, index = _fill_n2(mixture, wavelength, float(z))
peak_power = pulse_peak_power(energy, fwhm, shape)
beta2 = float(base.beta2(mode, np.array([omega]), z)[0])
gamma = nonlinear_coefficient(mode, n2, wavelength, z)
area = base.aeff(mode, z)
zdw = _zdw_near(mode, omega, float(z))
shared = mixture.shared
return SolitonParams(
energy=float(energy),
fwhm=float(fwhm),
width=pulse_width(fwhm, shape),
omega=omega,
wavelength=float(wavelength),
gas="+".join(mixture.labels),
pressure=math.nan if shared is None else float(shared(here)[0]),
# _resolve_fill refused a temperature that varies along z.
temperature=cast("float", mixture.temperature),
shape=shape,
peak_power=peak_power,
beta2=beta2,
density=float(mixture.total_density(here)[0]),
index=index,
n2=n2,
gamma=gamma,
soliton_order=soliton_order(peak_power, fwhm, gamma, beta2, shape),
soliton_content=(
soliton_content(peak_power, fwhm, gamma, beta2, shape)
if beta2 < 0
else None
),
intensity=peak_power / area,
peak_intensity=_on_axis_intensity(mode, peak_power, float(z)),
critical_power=self_focusing_power(wavelength, index, n2),
dispersion_length=dispersion_length(fwhm, beta2, shape),
nonlinear_length=nonlinear_length(peak_power, gamma),
fission_length=fission_length(peak_power, fwhm, gamma, beta2, shape),
zdw=zdw,
loss_length=float(base.loss_length(mode, np.array([omega]), z)[0]),
effective_area=area,
mode=mode,
)
def _capillary( # noqa: PLR0913 - the capillary's keywords
radius: float,
mixture: Mixture,
*,
n: int,
m: int,
kind: ModeKind,
model: MarcatiliModel,
loss: LossSpec,
cladding: GlassIndex | complex | None,
) -> MarcatiliMode:
"""Build the HE, TE or TM capillary the capillary design tools report on.
Parameters
----------
radius : float
Core radius in m, > 0.
mixture : Mixture
The resolved fill.
n : int
Azimuthal mode index.
m : int
Radial mode index, >= 1.
kind : {'HE', 'TE', 'TM'}
Mode family.
model : {'full', 'reduced'}
Marcatili model.
loss : bool, float or callable
Loss policy of the mode.
cladding : GlassIndex, complex or None
Cladding index; ``None`` is fused silica.
Returns
-------
MarcatiliMode
The mode, as `prop_capillary` builds it for the same arguments.
"""
return MarcatiliMode(
_positive("radius", radius),
mixture=mixture,
n=n,
m=m,
kind=kind,
model=model,
loss=loss,
cladding=_default_cladding(cladding),
)
[docs]
def capillary_params( # noqa: PLR0913 - the capillary's keywords beside the pulse's
radius: float,
gas: Fill,
pressure: float | None = None,
*,
energy: float,
fwhm: float,
wavelength: float,
shape: PulseShape = "sech",
temperature: float | None = None,
n: int = 1,
m: int = 1,
kind: ModeKind = "HE",
model: MarcatiliModel = "full",
loss: LossSpec = True,
cladding: GlassIndex | complex | None = None,
materials: MaterialSelection | None = None,
) -> SolitonParams:
"""Return the design numbers of a pulse in a gas-filled capillary.
Builds the mode as :func:`solprop.driver.interface.prop_capillary` does for the same
arguments -- a :class:`~solprop.modes.capillary.MarcatiliMode` over the fill's
mixture -- and hands it to :func:`soliton_params`.
Parameters
----------
radius : float
Core radius in m, > 0.
gas : Gas, str, Mapping or Mixture
The fill, as `prop_capillary` takes it.
pressure : float or None, optional
Total pressure in **bar**; ``None`` is one bar. For a gradient, evaluate at
the pressure of the position of interest.
energy : float
Pulse energy in J, > 0.
fwhm : float
Full width at half maximum of the power in s, > 0.
wavelength : float
Pump wavelength in m, > 0.
shape : {'sech', 'gauss'}, optional
The pulse shape.
temperature : float or None, optional
Fill temperature in K; ``None`` is room temperature, or the mixture's own
beside a `Mixture`, which refuses any other.
n : int, optional
Azimuthal mode index.
m : int, optional
Radial mode index, >= 1.
kind : {'HE', 'TE', 'TM'}, optional
Mode family.
model : {'full', 'reduced'}, optional
Marcatili model.
loss : bool, float or callable, optional
Loss policy of the mode.
cladding : GlassIndex, complex or None, optional
Cladding index; ``None`` is fused silica.
materials : MaterialSelection or None, optional
Build the capillary as `prop_capillary` does under the same ``materials=``:
the fill from optmat (index, density and Kerr coefficient at `wavelength`)
and ``cladding=None`` as optmat's recommended fused silica. ``None`` is
``MaterialSelection()``.
Returns
-------
SolitonParams
The numbers, with the mode that was built.
Examples
--------
>>> from solprop.tools import capillary_params
>>> p = capillary_params(125e-6, "He", 0.4, energy=300e-6, fwhm=10e-15,
... wavelength=800e-9, model="reduced", cladding=1.45)
>>> f"{p.gamma:.3e} {p.intensity:.3e}"
'4.884e-11 1.130e+18'
"""
materials = MaterialSelection() if materials is None else materials
_check_carrier(materials, wavelength)
mixture = _resolve_fill(gas, pressure, temperature, materials)
mode = _capillary(
radius,
mixture,
n=n,
m=m,
kind=kind,
model=model,
loss=loss,
cladding=cladding,
)
return soliton_params(
mode, mixture, energy=energy, fwhm=fwhm, wavelength=wavelength, shape=shape
)
def _area_per_order(shape: PulseShape) -> float:
"""Return the pulse area over pi per unit conventional soliton order.
Parameters
----------
shape : {'sech', 'gauss'}
The pulse shape.
Returns
-------
float
One for sech, ``sqrt(2 / pi)`` for a Gaussian (:func:`soliton_content`).
"""
return 1.0 if shape == "sech" else math.sqrt(2.0 / math.pi)
def _gaussian_order(eigenvalue: float) -> float:
"""Return the conventional order whose Gaussian has a given largest eigenvalue.
Parameters
----------
eigenvalue : float
Target imaginary part of the largest Zakharov-Shabat eigenvalue, > 0.
Returns
-------
float
The conventional soliton order.
Raises
------
ToolsError
If the eigenvalue solve fails or the bracket does not hold the root.
"""
threshold = 0.5 / _area_per_order("gauss")
def residual(order: float) -> float:
"""Return the largest eigenvalue at an order less the target.
Parameters
----------
order : float
Conventional soliton order.
Returns
-------
float
The difference, dimensionless.
"""
# At and below the first threshold there is no bound state: the eigenvalue
# that appears there starts from zero.
if order <= threshold:
return -eigenvalue
return _gaussian_eigenvalue(order) - eigenvalue
# N - eta of the Gaussian falls from sqrt(pi/8) = 0.627 at the first threshold
# towards 1/2 as N grows (0.579 at N = 1, 0.506 at N = 10, measured with
# _gaussian_eigenvalue), so the root lies between eta + 1/2 and eta + 0.627. The
# upper end carries a margin; brentq refuses a bracket that does not hold the root.
low = max(threshold, eigenvalue + 0.5)
high = eigenvalue + 0.7
try:
return float(brentq(residual, low, high, xtol=1e-12, rtol=1e-11))
except ValueError as error:
raise ToolsError(
f"no Gaussian soliton order between {low:.6g} and {high:.6g} has a largest "
f"eigenvalue of {eigenvalue:.6g}: {error}"
) from error
def _target_order(
order: float | None,
count: int | None,
eigenvalue: float | None,
shape: PulseShape,
) -> float:
"""Return the conventional soliton order that meets exactly one design target.
Parameters
----------
order : float or None
A conventional soliton order, > 0.
count : int or None
A number of solitons, >= 1.
eigenvalue : float or None
The largest Zakharov-Shabat eigenvalue, > 0.
shape : {'sech', 'gauss'}
The pulse shape.
Returns
-------
float
The conventional soliton order.
Raises
------
ToolsError
If not exactly one target is given, or it is out of range.
"""
given = [
name
for name, value in (
("soliton_order", order),
("soliton_count", count),
("largest_eigenvalue", eigenvalue),
)
if value is not None
]
if len(given) != 1:
raise ToolsError(
"give exactly one of soliton_order, soliton_count and largest_eigenvalue, "
f"got {', '.join(given) or 'none'}"
)
if order is not None:
return _positive("soliton_order", order)
if count is not None:
if count < 1:
raise ToolsError(f"soliton_count must be at least 1, got {count}")
# The count is M = ceil(S - 1/2) for an area S (soliton_content), so the
# pulses with M solitons have S in (M - 1/2, M + 1/2]. The design takes the
# middle, S = M: the pure M-soliton for sech, the farthest from either
# threshold for a Gaussian.
return count / _area_per_order(shape)
eta = _positive("largest_eigenvalue", cast("float", eigenvalue))
return eta + 0.5 if shape == "sech" else _gaussian_order(eta)
[docs]
def capillary_design( # noqa: PLR0913 - the capillary's keywords beside the pulse's and targets
radius: float,
gas: Gas | str,
*,
zdw_wavelength: float,
fwhm: float,
wavelength: float,
soliton_order: float | None = None,
soliton_count: int | None = None,
largest_eigenvalue: float | None = None,
shape: PulseShape = "sech",
temperature: float = ROOM_TEMPERATURE,
m: int = 1,
model: MarcatiliModel = "full",
loss: LossSpec = True,
cladding: GlassIndex | complex | None = None,
materials: MaterialSelection | None = None,
pressure_max: float = 100.0,
) -> SolitonParams:
r"""Return the pressure and energy that give a capillary a soliton design.
The inverse of :func:`capillary_params`: the pressure puts the zero-dispersion
point at `zdw_wavelength` (:func:`pressure_zdw`), and the energy then gives the
pulse one soliton target. At that pressure :math:`\beta_2` and :math:`\gamma` do
not depend on the energy and :math:`N^2 = \gamma P_0\tau_0^2/|\beta_2|`, so the
peak power is :math:`P_0 = N^2|\beta_2|/(\gamma\tau_0^2)`. Give exactly one
target:
- `soliton_order`, the conventional :math:`N`;
- `soliton_count`, a number :math:`M` of solitons in the scattering content
(:func:`soliton_content`): the pulse area is set to :math:`M\pi`, the middle
of the band of areas with :math:`M` solitons, which is :math:`N = M` for
sech and :math:`N = M\sqrt{\pi/2}` for a Gaussian;
- `largest_eigenvalue`, the :math:`\eta` of the largest constituent, which
fixes its width :math:`\tau_0/(2\eta)`: :math:`N = \eta + 1/2` for sech, and
a root of the Gaussian scattering problem otherwise.
Parameters
----------
radius : float
Core radius in m, > 0.
gas : Gas or str
A single built-in gas, as :func:`pressure_zdw` takes it.
zdw_wavelength : float
Wavelength at which the dispersion should vanish, m.
fwhm : float
Full width at half maximum of the power in s, > 0.
wavelength : float
Pump wavelength in m, > 0.
soliton_order : float or None, optional
Target conventional soliton order, > 0.
soliton_count : int or None, optional
Target number of solitons, >= 1.
largest_eigenvalue : float or None, optional
Target imaginary part of the largest Zakharov-Shabat eigenvalue,
dimensionless with time in units of the shape's natural width, > 0.
shape : {'sech', 'gauss'}, optional
The pulse shape, used for the width, the peak power and the content.
temperature : float, optional
Fill temperature in K.
m : int, optional
Radial index of the HE_1m mode.
model : {'full', 'reduced'}, optional
Marcatili model.
loss : bool, float or callable, optional
Loss policy of the mode in the returned numbers.
cladding : GlassIndex, complex or None, optional
Cladding index; ``None`` is fused silica.
materials : MaterialSelection or None, optional
Design under this selection, as :func:`capillary_params` and
:func:`pressure_zdw` do. ``None`` is ``MaterialSelection()``.
pressure_max : float, optional
Upper end of the pressure bracket in **bar**.
Returns
-------
SolitonParams
The numbers at the designed pressure and energy, with the mode; the
pressure is ``.pressure`` (bar) and the energy ``.energy`` (J).
Raises
------
ToolsError
If not exactly one target is given or it is out of range, no pressure in the
bracket gives the zero-dispersion point, or a content target is asked of a
pump in normal dispersion.
Examples
--------
>>> from solprop.tools import capillary_design
>>> p = capillary_design(125e-6, "Ar", zdw_wavelength=500e-9, fwhm=10e-15,
... wavelength=800e-9, soliton_order=3.0, cladding=1.45)
>>> f"{p.pressure:.4f} bar {p.energy * 1e6:.1f} uJ N = {p.soliton_order:.3f}"
'0.0664 bar 100.5 uJ N = 3.000'
>>> g = capillary_design(125e-6, "Ar", zdw_wavelength=500e-9, fwhm=10e-15,
... wavelength=800e-9, soliton_count=2, shape="gauss",
... cladding=1.45)
>>> f"{g.soliton_order:.4f} {g.soliton_content.count}"
'2.5066 2'
"""
materials = MaterialSelection() if materials is None else materials
_check_carrier(materials, wavelength)
_check_shape(shape)
order = _target_order(soliton_order, soliton_count, largest_eigenvalue, shape)
pressure = pressure_zdw(
radius,
gas,
zdw_wavelength,
pressure_max=pressure_max,
temperature=temperature,
m=m,
model=model,
cladding=cladding,
materials=materials,
)
mixture = _resolve_fill(gas, pressure, temperature, materials)
mode = _capillary(
radius,
mixture,
n=1,
m=m,
kind="HE",
model=model,
loss=loss,
cladding=cladding,
)
beta2 = float(base.beta2(mode, np.array([_omega(wavelength)]))[0])
if soliton_order is None and beta2 >= 0.0:
raise ToolsError(
f"a soliton count or eigenvalue needs anomalous dispersion at the pump, "
f"but beta2 = {beta2:.4g} s^2/m at {wavelength * 1e9:.4g} nm with the "
f"zero-dispersion point at {zdw_wavelength * 1e9:.4g} nm"
)
n2, _ = _fill_n2(mixture, wavelength, 0.0)
gamma = nonlinear_coefficient(mode, n2, wavelength)
peak_power = order**2 / (gamma * dispersion_length(fwhm, beta2, shape))
# The peak power is linear in the energy at a fixed width and shape.
energy = peak_power / pulse_peak_power(1.0, fwhm, shape)
return soliton_params(
mode, mixture, energy=energy, fwhm=fwhm, wavelength=wavelength, shape=shape
)
[docs]
def gas_ratio(
gas_a: Gas | str,
gas_b: Gas | str,
wavelength: float,
*,
pressure: float = 1.0,
temperature: float = ROOM_TEMPERATURE,
materials: MaterialSelection | None = None,
) -> GasRatio:
"""Return the ratios of bulk dispersion and Kerr response of two gases.
The bulk (unguided) dispersion of each gas, not a mode's: it answers "how much
more pressure of B matches A", which is a property of the gases alone.
Parameters
----------
gas_a : Gas or str
The numerator's species.
gas_b : Gas or str
The denominator's species.
wavelength : float
Wavelength in m at which the dispersions are compared, > 0.
pressure : float, optional
Pressure of both gases in **bar**, > 0.
temperature : float, optional
Temperature of both in K.
materials : MaterialSelection or None, optional
Take both gases from optmat under this selection, as a run's ``materials=``
does: their index, their density and their Kerr coefficient at
`wavelength`; ``"air"`` is nitrogen and oxygen. ``None`` is
``MaterialSelection()``.
Returns
-------
GasRatio
Named ``beta2`` and ``chi3`` ratios.
Raises
------
ToolsError
If the selection names a carrier other than `wavelength`, or a gas has no
Kerr data under it.
Examples
--------
>>> from solprop.tools import gas_ratio
>>> ratio = gas_ratio("Ar", "He", 800e-9)
>>> f"{ratio.beta2:.1f} {ratio.chi3:.1f}"
'18.6 28.4'
"""
omega = _omega(wavelength)
materials = MaterialSelection() if materials is None else materials
_check_carrier(materials, wavelength)
dispersions = []
susceptibilities = []
for gas in (gas_a, gas_b):
fill = _resolve_fill(
gas, _positive("pressure", pressure), temperature, materials
)
dispersions.append(_bulk_beta2(_BulkGas(fill), omega))
susceptibilities.append(_fill_chi3(fill, wavelength))
return GasRatio(
beta2=dispersions[0] / dispersions[1],
chi3=susceptibilities[0] / susceptibilities[1],
)
@dataclass(frozen=True)
class _BulkGas:
"""A uniform gas as an index model, for bulk dispersion.
Attributes
----------
fill : Mixture
The gas resolved under a material selection, whose index is taken.
"""
fill: Mixture
def excess_beta(
self, omega: Float[np.ndarray, " n_omega"]
) -> Float[np.ndarray, " n_omega"]:
"""Return the gas's part of the bulk propagation constant.
The vacuum's ``omega / c`` has no second derivative, and leaving it out of a
difference quotient removes its rounding, which at a gas's ``n - 1`` of 1e-5 is
the larger part of the quotient.
Parameters
----------
omega : numpy.ndarray
Angular frequencies in rad/s.
Returns
-------
numpy.ndarray
``omega (n - 1) / c`` in 1/m, with ``n - 1`` formed without cancellation.
"""
chi = self.fill.susceptibility(omega, np.zeros(1))[0]
excess = chi / (1.0 + np.sqrt(1.0 + chi))
return np.asarray(omega * np.real(excess) / C)
def _bulk_beta2(bulk: _BulkGas, omega: float) -> float:
"""Return the bulk group-velocity dispersion by a Richardson central difference.
Parameters
----------
bulk : _BulkGas
The gas.
omega : float
Angular frequency in rad/s.
Returns
-------
float
``d^2 beta / d omega^2`` in s^2/m.
"""
# A second difference of the gas's part of beta: only a ratio of two such numbers is
# used, so the few parts in 1e5 this reaches at the step modes.base chose are ample.
step = base.BETA2_RELATIVE_STEP * omega
def second(h: float) -> float:
"""Return the three-point second difference at step ``h``.
Parameters
----------
h : float
Step in rad/s.
Returns
-------
float
The difference quotient in s^2/m.
"""
values = bulk.excess_beta(np.array([omega - h, omega, omega + h]))
return float(values[0] - 2.0 * values[1] + values[2]) / h**2
return (4.0 * second(step / 2.0) - second(step)) / 3.0
# -------------------------------------------------------------------- phase matching
def _check_phase(nonlinear_phase: str) -> None:
"""Refuse a nonlinear phase model the phase matching does not define.
Parameters
----------
nonlinear_phase : str
The phase keyword.
Raises
------
ToolsError
If it is not ``"largest"``, ``"supplied_power"`` or ``"linear"``.
"""
if nonlinear_phase not in ("largest", "supplied_power", "linear"):
raise ToolsError(
"nonlinear_phase must be 'largest', 'supplied_power' or 'linear', "
f"got {nonlinear_phase!r}"
)
@dataclass(frozen=True)
class _PhasePulse:
"""Validate the pulse declaration before a phase-matching search."""
peak_power: float | None
gamma: float | None
fwhm: float | None
shape: PulseShape
nonlinear_phase: NonlinearPhase
def __post_init__(self) -> None:
"""Refuse incomplete or contradictory phase declarations.
Raises
------
ToolsError
If the phase choice or pulse declaration is invalid.
"""
_check_shape(self.shape)
_check_phase(self.nonlinear_phase)
if (self.peak_power is None) != (self.gamma is None):
raise ToolsError(
"the nonlinear term needs both peak_power and gamma; "
"give both or neither"
)
if self.nonlinear_phase == "linear" and (
self.peak_power is not None or self.fwhm is not None
):
raise ToolsError(
"nonlinear_phase='linear' requires omitting pulse power, gamma and fwhm"
)
if self.peak_power is None:
if self.fwhm is not None:
raise ToolsError("fwhm requires peak_power and gamma")
return
_positive("peak_power", self.peak_power)
assert self.gamma is not None # The paired declaration was validated above.
_positive("gamma", self.gamma)
if self.nonlinear_phase == "largest":
if self.fwhm is None:
raise ToolsError(
"the largest constituent needs fwhm; supply it or select "
"nonlinear_phase='supplied_power'"
)
_positive("fwhm", self.fwhm)
elif self.fwhm is not None:
raise ToolsError(
"nonlinear_phase='supplied_power' does not use fwhm; omit it"
)
def wavenumber(self, mode: Mode, omega: float, z: float) -> float:
"""Calculate the nonlinear phase in the soliton mode.
Parameters
----------
mode : Mode
Mode carrying the soliton.
omega : float
Pump angular frequency in rad/s.
z : float
Position in m.
Returns
-------
float
Nonlinear wavenumber in 1/m.
Raises
------
ToolsError
If the largest-constituent model has no discrete soliton.
"""
if self.peak_power is None or self.gamma is None:
return 0.0
if self.fwhm is None:
return 0.5 * self.gamma * self.peak_power
dispersion = float(base.beta2(mode, np.array([omega]), z)[0])
content = soliton_content(
self.peak_power, self.fwhm, self.gamma, dispersion, self.shape
)
if content.largest_nonlinear_wavenumber is None:
raise ToolsError(
"the pulse has no discrete soliton for "
"largest-constituent phase matching"
)
return content.largest_nonlinear_wavenumber
@dataclass(frozen=True)
class _PhaseMismatch:
"""Evaluate a mismatch with the pump and nonlinear phase prepared once."""
mode: Mode
omega0: float
beta0: float
slope: float
nonlinear: float
z: float
def __call__(
self, omega: Float[np.ndarray, " n_omega"]
) -> Float[np.ndarray, " n_omega"]:
"""Evaluate the mismatch at trial frequencies.
Parameters
----------
omega : numpy.ndarray
Angular frequencies in rad/s.
Returns
-------
numpy.ndarray
Mismatch in 1/m.
"""
return np.asarray(
base.beta(self.mode, omega, self.z)
- self.beta0
- self.slope * (omega - self.omega0)
- self.nonlinear
)
def _prepare_mismatch(
mode: Mode,
wavelength: float,
soliton_mode: Mode | None,
z: float,
pulse: _PhasePulse,
) -> _PhaseMismatch:
"""Prepare a phase-matching calculation at one position.
Parameters
----------
mode : Mode
Dispersive-wave mode.
wavelength : float
Pump wavelength in m.
soliton_mode : Mode or None
Pump mode, or None to use mode.
z : float
Position in m.
pulse : _PhasePulse
Validated pulse and phase choice.
Returns
-------
_PhaseMismatch
Callable with the nonlinear phase evaluated once.
"""
soliton = mode if soliton_mode is None else soliton_mode
omega0 = _omega(wavelength)
pump = np.array([omega0])
return _PhaseMismatch(
mode,
omega0,
float(base.beta(soliton, pump, z)[0]),
float(base.beta1(soliton, pump, z)[0]),
pulse.wavenumber(soliton, omega0, z),
z,
)
[docs]
def phase_mismatch( # noqa: PLR0913 - two modes, a position and the two Kerr numbers
mode: Mode,
omega: Float[np.ndarray, " n_omega"],
wavelength: float,
*,
soliton_mode: Mode | None = None,
z: float = 0.0,
peak_power: float | None = None,
gamma: float | None = None,
fwhm: float | None = None,
shape: PulseShape = "sech",
nonlinear_phase: NonlinearPhase = "largest",
) -> Float[np.ndarray, " n_omega"]:
r"""Return the phase mismatch between a soliton and a dispersive wave.
.. math::
\Delta\beta(\omega) = \beta_{\rm dw}(\omega) - \beta_{\rm s}(\omega_0)
- \beta_{1,\rm s}(\omega_0)(\omega-\omega_0) - k_{\rm NL} .
Parameters
----------
mode : Mode
The mode the dispersive wave is emitted into.
omega : numpy.ndarray
Angular frequencies in rad/s, shape ``(n_omega,)``.
wavelength : float
Wavelength of the soliton in m, > 0.
soliton_mode : Mode or None, optional
The mode the soliton travels in; ``None`` is `mode`. A different one gives
intermodal phase matching.
z : float, optional
Position along the waveguide, m.
peak_power : float or None, optional
Input peak power in W. With `gamma` and `fwhm`, adds the soliton's nonlinear
wavenumber; both ``None`` (the default) gives the linear condition.
gamma : float or None, optional
Nonlinear coefficient of the soliton's mode in rad/(W m)
(:func:`nonlinear_coefficient`).
fwhm : float or None, optional
Input power FWHM in s, positive; required with power for the default phase.
shape : {'sech', 'gauss'}, optional
Unchirped input shape for the largest constituent.
nonlinear_phase : {'largest', 'supplied_power', 'linear'}, optional
Default uses the largest constituent. Supplied power uses gamma P0/2 and
requires omitting fwhm. Linear requires omitting all pulse numbers.
Without pulse numbers the default is also linear.
Returns
-------
numpy.ndarray
Mismatch in 1/m, shape ``(n_omega,)``. It vanishes at the pump in the linear
case, with zero slope there.
Raises
------
ToolsError
If pulse arguments are incomplete or contradictory, or the largest
constituent does not exist in anomalous dispersion.
Notes
-----
The default nonlinear term is gamma times the largest constituent's power / 2
(:func:`soliton_content`). The supplied-power option treats the quoted power as
a fundamental soliton, advancing as exp(i gamma P0 z / 2) (Akhmediev and
Karlsson, Phys. Rev. A 51, 2602 (1995)). Both are local estimates; loss,
higher-order dispersion, Raman, chirp and plasma evolution are not included.
Examples
--------
>>> import numpy as np
>>> from solprop.modes import MarcatiliMode
>>> from solprop.tools import phase_mismatch
>>> mode = MarcatiliMode(125e-6, "He", 1.0)
>>> pump = 2 * np.pi * 299792458.0 / 800e-9
>>> abs(float(phase_mismatch(mode, np.array([pump]), 800e-9)[0])) < 1e-6
True
"""
pulse = _PhasePulse(peak_power, gamma, fwhm, shape, nonlinear_phase)
prepared = _prepare_mismatch(mode, wavelength, soliton_mode, z, pulse)
return prepared(np.asarray(omega, dtype=float))
[docs]
def lambda_rdw_all( # noqa: PLR0913 - phase_mismatch's keywords and the search window
mode: Mode,
wavelength: float,
*,
soliton_mode: Mode | None = None,
z: float = 0.0,
peak_power: float | None = None,
gamma: float | None = None,
fwhm: float | None = None,
shape: PulseShape = "sech",
nonlinear_phase: NonlinearPhase = "largest",
wavelength_limits: tuple[float, float] | None = None,
) -> tuple[float, ...]:
"""Return every phase-matched wavelength in a window, shortest first.
Parameters
----------
mode : Mode
The mode the dispersive wave is emitted into.
wavelength : float
Wavelength of the soliton in m, > 0.
soliton_mode : Mode or None, optional
The soliton's mode; ``None`` is `mode`.
z : float, optional
Position along the waveguide, m.
peak_power : float or None, optional
Peak power in W for the nonlinear term; see :func:`phase_mismatch`.
gamma : float or None, optional
Nonlinear coefficient in rad/(W m) for the nonlinear term.
fwhm : float or None, optional
Input power FWHM in s, positive; required with power for the default phase.
shape : {'sech', 'gauss'}, optional
Unchirped input shape for the largest constituent.
nonlinear_phase : {'largest', 'supplied_power', 'linear'}, optional
Default uses the largest constituent. Supplied power uses gamma P0/2 and
requires omitting fwhm. Linear requires omitting all pulse numbers.
Without pulse numbers the default is also linear.
wavelength_limits : tuple of float or None, optional
Shortest and longest wavelength searched, m. ``None`` is 100 nm to 0.9 of the
pump wavelength, which keeps the double root at the pump out.
Returns
-------
tuple of float
The roots in m, ascending; empty when there is none. The mismatch is sampled
on 2048 frequencies across the window, each sign change is refined with
Brent's method, and roots closer than one sample are not separated.
Raises
------
ToolsError
If the window is not positive and ordered.
Examples
--------
>>> from solprop.modes import MarcatiliMode
>>> from solprop.tools import lambda_rdw_all
>>> mode = MarcatiliMode(125e-6, "He", 1.0, cladding=1.45)
>>> [f"{root * 1e9:.0f}" for root in lambda_rdw_all(mode, 800e-9)]
['200']
"""
pulse = _PhasePulse(peak_power, gamma, fwhm, shape, nonlinear_phase)
prepared = _prepare_mismatch(mode, wavelength, soliton_mode, z, pulse)
limits = (
(100e-9, 0.9 * wavelength) if wavelength_limits is None else wavelength_limits
)
shortest, longest = limits
if not 0.0 < shortest < longest:
raise ToolsError(
f"the wavelength limits must be positive and ordered: {limits}"
)
grid = np.linspace(
2.0 * math.pi * C / longest, 2.0 * math.pi * C / shortest, _MISMATCH_SAMPLES
)
def mismatch(value: float) -> float:
"""Return the mismatch at one frequency, for the root finder.
Parameters
----------
value : float
Angular frequency in rad/s.
Returns
-------
float
Mismatch in 1/m.
"""
return float(prepared(np.array([value]))[0])
with _probing():
sampled = prepared(grid)
if not np.all(np.isfinite(sampled)):
# A sample that is not a number compares false with everything, so it
# would delete any root beside it without a word.
missing = int(np.sum(~np.isfinite(sampled)))
raise ToolsError(
f"the phase mismatch is not finite at {missing} of {sampled.size} "
f"samples in {shortest * 1e9:.4g}-{longest * 1e9:.4g} nm: the index "
"model does not cover the window"
)
signs = np.sign(sampled)
found = [float(grid[i]) for i in np.flatnonzero(signs == 0.0)]
found.extend(
float(brentq(mismatch, float(grid[i]), float(grid[i + 1])))
for i in np.flatnonzero(signs[:-1] * signs[1:] < 0.0)
)
if found:
# Outside the probing context, so a root beyond the index model's range warns.
base.beta(mode, np.array(sorted(found)), z)
return tuple(sorted(2.0 * math.pi * C / root for root in found))
[docs]
def lambda_rdw( # noqa: PLR0913 - phase_mismatch's keywords and the search window
mode: Mode,
wavelength: float,
*,
soliton_mode: Mode | None = None,
z: float = 0.0,
peak_power: float | None = None,
gamma: float | None = None,
fwhm: float | None = None,
shape: PulseShape = "sech",
nonlinear_phase: NonlinearPhase = "largest",
wavelength_limits: tuple[float, float] | None = None,
) -> float:
"""Return the phase-matched dispersive-wave wavelength nearest the pump.
Parameters
----------
mode : Mode
The mode the dispersive wave is emitted into.
wavelength : float
Wavelength of the soliton in m, > 0.
soliton_mode : Mode or None, optional
The soliton's mode; ``None`` is `mode`.
z : float, optional
Position along the waveguide, m.
peak_power : float or None, optional
Peak power in W for the nonlinear term; see :func:`phase_mismatch`.
gamma : float or None, optional
Nonlinear coefficient in rad/(W m) for the nonlinear term.
fwhm : float or None, optional
Input power FWHM in s, positive; required with power for the default phase.
shape : {'sech', 'gauss'}, optional
Unchirped input shape for the largest constituent.
nonlinear_phase : {'largest', 'supplied_power', 'linear'}, optional
Default uses the largest constituent. Supplied power uses gamma P0/2 and
requires omitting fwhm. Linear requires omitting all pulse numbers.
Without pulse numbers the default is also linear.
wavelength_limits : tuple of float or None, optional
Window searched, m; ``None`` is 100 nm to 0.9 of the pump wavelength.
Returns
-------
float
The longest root inside the window in m, which with the default window is the
short-wavelength dispersive wave nearest the pump. :func:`lambda_rdw_all`
returns all of them.
Raises
------
ToolsError
If the window holds no root, rather than returning a missing value that would
propagate silently.
Examples
--------
>>> from solprop.modes import MarcatiliMode
>>> from solprop.tools import lambda_rdw
>>> mode = MarcatiliMode(125e-6, "He", 1.0, cladding=1.45)
>>> f"{lambda_rdw(mode, 800e-9) * 1e9:.1f}"
'200.1'
"""
roots = lambda_rdw_all(
mode,
wavelength,
soliton_mode=soliton_mode,
z=z,
peak_power=peak_power,
gamma=gamma,
fwhm=fwhm,
shape=shape,
nonlinear_phase=nonlinear_phase,
wavelength_limits=wavelength_limits,
)
if not roots:
limits = (
(100e-9, 0.9 * wavelength)
if wavelength_limits is None
else wavelength_limits
)
ends = phase_mismatch(
mode,
2.0 * math.pi * C / np.array(limits[::-1]),
wavelength,
soliton_mode=soliton_mode,
z=z,
peak_power=peak_power,
gamma=gamma,
fwhm=fwhm,
shape=shape,
nonlinear_phase=nonlinear_phase,
)
raise ToolsError(
f"no phase-matched wavelength in {limits[0] * 1e9:.4g}-"
f"{limits[1] * 1e9:.4g} nm: the mismatch is {ends[1]:.4g} 1/m at the short "
f"end and {ends[0]:.4g} 1/m at the long end and never changes sign"
)
return roots[-1]
def _pressure_root(
residual: Callable[[float], float], pressure_max: float, what: str
) -> float:
"""Return the pressure at which a function of pressure vanishes.
Parameters
----------
residual : callable
A function of the pressure in bar.
pressure_max : float
Upper end of the bracket in bar.
what : str
What is being matched, for the refusal.
Returns
-------
float
The root in bar.
Raises
------
ToolsError
If the function has the same sign at both ends of the bracket.
"""
top = _positive("pressure_max", pressure_max)
if top <= _PRESSURE_FLOOR:
raise ToolsError(f"pressure_max must exceed {_PRESSURE_FLOOR} bar, got {top}")
with _probing():
low, high = residual(_PRESSURE_FLOOR), residual(top)
if not (math.isfinite(low) and math.isfinite(high)):
raise ToolsError(
f"the residual for {what} is not finite at the ends of the pressure "
f"bracket ({low!r}, {high!r}): the index model does not reach that far"
)
if low * high > 0.0:
raise ToolsError(
f"no pressure between {_PRESSURE_FLOOR} and {top:.4g} bar gives {what}: "
f"the residual is {low:.4g} at the low end and {high:.4g} at the high end"
)
with _probing():
root = float(brentq(residual, _PRESSURE_FLOOR, top))
# Outside the probing context, so the answer's own warnings are raised.
residual(root)
return root
def _design_capillary( # noqa: PLR0913 - the capillary's keywords
radius: float,
gas: Gas | str,
*,
temperature: float,
m: int,
model: MarcatiliModel,
cladding: GlassIndex | complex | None,
lattice: TubeLattice | None,
materials: MaterialSelection | None,
) -> MarcatiliMode:
"""Build the single-fill capillary a pressure search re-evaluates.
The mode `prop_capillary` builds with the same ``materials=``: the gas (or
``"air"``) as an optmat fill at one bar and ``cladding=None`` as optmat's
recommended silica.
Parameters
----------
radius : float
Core radius in m, > 0.
gas : Gas or str
The species.
temperature : float
Fill temperature in K.
m : int
Radial index of the HE_1m mode.
model : {'full', 'reduced'}
Marcatili model.
cladding : GlassIndex, complex or None
Cladding index; ``None`` is fused silica.
lattice : TubeLattice or None
Tube geometry for a VincettiMode; None builds a Marcatili capillary.
materials : MaterialSelection or None
The design's material selection.
Returns
-------
MarcatiliMode
The capillary at one bar.
"""
bore = _positive("radius", radius)
wall = _default_cladding(cladding)
# The fill is an optmat mixture, which carries its own pressure and temperature.
mixture = _resolve_fill(gas, 1.0, temperature, materials)
if lattice is None:
return MarcatiliMode(
bore, None, 0.0, m=m, model=model, cladding=wall, mixture=mixture
)
return VincettiMode(
bore,
None,
0.0,
m=m,
model=model,
cladding=wall,
mixture=mixture,
lattice=lattice,
)
def _capillary_at(mode: MarcatiliMode, pressure: float) -> MarcatiliMode:
"""Return a single-fill capillary re-evaluated at another pressure.
The mode is built once and frozen at each density through its parameterised form,
so a root finder costs an equation-of-state call and an index evaluation per step
rather than a mode construction. A selected fill takes its density from its
species' own equation of state (optmat's); a selected air, two populations at
one shared pressure, is frozen at that pressure.
Parameters
----------
mode : MarcatiliMode
The capillary, filled at any pressure.
pressure : float
Pressure in bar.
Returns
-------
MarcatiliMode
The same capillary at that pressure.
"""
# pressure_rdw and pressure_zdw build this mode at their one scalar temperature.
temperature = cast("float", mode.temperature)
# `_design_capillary` always builds the mode on an optmat fill.
fill = cast("Mixture", mode.mixture)
if len(fill.components) > 1:
return mode.at_parameters({"pressure": float(pressure)})
state = next(iter(fill.components)).species.equation_of_state
rho = float(state.density(float(pressure), temperature))
return mode.at_parameters({"density": rho})
def _rdw_pressure_root(
residual: Callable[[float], float | None], limits: tuple[float, float]
) -> float:
"""Find a unique pressure root without crossing a sampled inadmissible region.
Parameters
----------
residual : callable
Mismatch in 1/m, or None only where no focusing discrete soliton exists.
limits : tuple of float
Positive ordered pressure bracket in bar.
Returns
-------
float
Unique pressure root in bar.
Raises
------
ToolsError
If the bracket is invalid, the mismatch is nonfinite, a refinement enters
an inadmissible region, or there is no unique sampled root.
"""
low, high = limits
if not (math.isfinite(low) and math.isfinite(high) and 0 < low < high):
raise ToolsError("pressure_limits must be finite, positive and ordered")
# The logarithmic grid resolves the low-pressure end; the linear grid resolves
# high-pressure structure. As with the wavelength search, sub-grid pairs of
# roots or admissible intervals are not guaranteed; narrow the explicit bracket.
grid = np.unique(
np.concatenate((np.geomspace(low, high, 129), np.linspace(low, high, 129)))
)
def checked(pressure: float) -> float | None:
"""Evaluate a pressure without disguising numerical failures as invalid physics.
Parameters
----------
pressure : float
Trial pressure in bar.
Returns
-------
float or None
Finite mismatch, or an explicitly inadmissible state.
Raises
------
ToolsError
If the mismatch is nonfinite.
"""
value = residual(pressure)
if value is not None and not math.isfinite(value):
raise ToolsError(
f"the RDW pressure residual is not finite at {pressure:g} bar"
)
return value
def refine(pressure: float) -> float:
"""Evaluate a root bracket inside one admissible interval.
Parameters
----------
pressure : float
Trial pressure in bar.
Returns
-------
float
Mismatch in 1/m.
Raises
------
ToolsError
If a bracket crossed an unresolved inadmissible region.
"""
value = checked(pressure)
if value is None:
raise ToolsError(
"pressure refinement crossed a no-soliton region; "
"narrow pressure_limits"
)
return value
with _probing():
sampled = [checked(float(pressure)) for pressure in grid]
roots = [
float(pressure)
for pressure, value in zip(grid, sampled, strict=True)
if value == 0.0
]
for i, (left, right) in enumerate(pairwise(sampled)):
if left is not None and right is not None and left * right < 0:
roots.append(float(brentq(refine, float(grid[i]), float(grid[i + 1]))))
if not roots:
raise ToolsError(
f"no pressure between {low:g} and {high:g} bar gives an admissible "
"RDW root; narrow pressure_limits if needed"
)
if len(roots) > 1:
found = ", ".join(f"{root:.4g}" for root in sorted(roots))
raise ToolsError(
f"ambiguous RDW pressure: {len(roots)} roots, at {found} bar; "
"choose one with pressure_limits"
)
# Re-evaluate the accepted answer outside the probing warning filter.
refine(roots[0])
return roots[0]
[docs]
def pressure_rdw( # noqa: PLR0913 - the capillary's keywords beside the two wavelengths
radius: float,
gas: Gas | str,
target_wavelength: float,
wavelength: float,
*,
pressure_max: float = 100.0,
peak_power: float | None = None,
fwhm: float | None = None,
shape: PulseShape = "sech",
nonlinear_phase: NonlinearPhase = "largest",
lattice: TubeLattice | None = None,
pressure_limits: tuple[float, float] | None = None,
temperature: float = ROOM_TEMPERATURE,
m: int = 1,
model: MarcatiliModel = "full",
cladding: GlassIndex | complex | None = None,
materials: MaterialSelection | None = None,
) -> float:
"""Return the pressure that phase matches a dispersive wave at a target wavelength.
Parameters
----------
radius : float
Core radius in m, > 0.
gas : Gas or str
A single built-in gas.
target_wavelength : float
Wavelength the dispersive wave should have, m.
wavelength : float
Pump wavelength in m.
pressure_max : float, optional
Upper end of the pressure bracket in **bar**; the lower end is 1e-6 bar.
peak_power : float or None, optional
Peak power in W. When given, the soliton's nonlinear wavenumber is included
and its nonlinear coefficient follows the pressure.
fwhm : float or None, optional
Positive input power FWHM in s, required with power for the default phase.
shape : {'sech', 'gauss'}, optional
Unchirped input shape.
nonlinear_phase : {'largest', 'supplied_power', 'linear'}, optional
Largest constituent by default; supplied power preserves gamma P0/2
without fwhm. Linear requires omitting pulse numbers.
lattice : TubeLattice or None, optional
Tube geometry for VincettiMode; None builds a Marcatili capillary.
pressure_limits : tuple of float or None, optional
Positive ordered search range in bar, within (1e-6, pressure_max).
None uses that whole range. Searches 129 logarithmic and 129 linear
samples, refining admissible sign changes; sub-grid root pairs or valid
intervals can be missed. Narrow the range to resolve them.
temperature : float, optional
Fill temperature in K.
m : int, optional
Radial index of the HE_1m mode.
model : {'full', 'reduced'}, optional
Marcatili model.
cladding : GlassIndex, complex or None, optional
Cladding index; ``None`` is fused silica.
materials : MaterialSelection or None, optional
Design the capillary `prop_capillary` builds under the same
``materials=``: the gas from optmat (index, density, and the Kerr
coefficient at `wavelength` for the nonlinear phase) and ``cladding=None``
as optmat's recommended fused silica. ``None`` is
``MaterialSelection()``.
Returns
-------
float
Pressure in bar.
Raises
------
ToolsError
If no pressure in the bracket matches, or the selection names a carrier
other than `wavelength`.
Examples
--------
>>> from solprop.tools import pressure_rdw
>>> f"{pressure_rdw(125e-6, 'He', 199.75e-9, 800e-9, cladding=1.45):.2f}"
'1.00'
"""
# Validate before scanning, using a placeholder gamma; the real one follows p.
_PhasePulse(
peak_power, None if peak_power is None else 1.0, fwhm, shape, nonlinear_phase
)
top = _positive("pressure_max", pressure_max)
if top <= _PRESSURE_FLOOR:
raise ToolsError(f"pressure_max must exceed {_PRESSURE_FLOOR} bar, got {top}")
limits = (_PRESSURE_FLOOR, top) if pressure_limits is None else pressure_limits
if limits[0] < _PRESSURE_FLOOR or limits[1] > top:
raise ToolsError("pressure_limits must lie within (1e-6, pressure_max) bar")
materials = MaterialSelection() if materials is None else materials
_check_carrier(materials, wavelength)
base_mode = _design_capillary(
radius,
gas,
temperature=temperature,
m=m,
model=model,
cladding=cladding,
lattice=lattice,
materials=materials,
)
target = np.array([_omega(target_wavelength)])
def residual(pressure: float) -> float | None:
"""Return the mismatch at the target for one pressure.
Parameters
----------
pressure : float
Pressure in bar.
Returns
-------
float or None
Mismatch in 1/m, or None when the constituent does not exist.
"""
mode = _capillary_at(base_mode, pressure)
gamma = None
if peak_power is not None:
assert mode.mixture is not None
index = _fill_n2(mode.mixture, wavelength)[0]
gamma = nonlinear_coefficient(mode, index, wavelength)
if peak_power is not None and fwhm is not None and gamma is not None:
dispersion = float(base.beta2(mode, np.array([_omega(wavelength)]))[0])
if not math.isfinite(dispersion):
raise ToolsError(
"the pump dispersion is not finite during pressure search"
)
if dispersion >= 0:
return None
order = soliton_order(peak_power, fwhm, gamma, dispersion, shape)
area = order if shape == "sech" else math.sqrt(2 / math.pi) * order
if area <= 0.5:
return None
return float(
phase_mismatch(
mode,
target,
wavelength,
peak_power=peak_power,
gamma=gamma,
fwhm=fwhm,
shape=shape,
nonlinear_phase=nonlinear_phase,
)[0]
)
return _rdw_pressure_root(residual, limits)
[docs]
def pressure_zdw( # noqa: PLR0913 - the capillary's keywords
radius: float,
gas: Gas | str,
zdw_wavelength: float,
*,
pressure_max: float = 100.0,
temperature: float = ROOM_TEMPERATURE,
m: int = 1,
model: MarcatiliModel = "full",
cladding: GlassIndex | complex | None = None,
materials: MaterialSelection | None = None,
) -> float:
"""Return the pressure that puts the zero-dispersion point at a wavelength.
Parameters
----------
radius : float
Core radius in m, > 0.
gas : Gas or str
A single built-in gas.
zdw_wavelength : float
Wavelength at which the dispersion should vanish, m.
pressure_max : float, optional
Upper end of the pressure bracket in **bar**; the lower end is 1e-6 bar.
temperature : float, optional
Fill temperature in K.
m : int, optional
Radial index of the HE_1m mode.
model : {'full', 'reduced'}, optional
Marcatili model.
cladding : GlassIndex, complex or None, optional
Cladding index; ``None`` is fused silica.
materials : MaterialSelection or None, optional
Design the capillary `prop_capillary` builds under the same
``materials=``: the gas's index and density from optmat and
``cladding=None`` as optmat's recommended fused silica. ``None`` is
``MaterialSelection()``.
Returns
-------
float
Pressure in bar.
Raises
------
ToolsError
If no pressure in the bracket does.
Examples
--------
>>> from solprop.tools import pressure_zdw
>>> f"{pressure_zdw(125e-6, 'Ar', 500e-9, cladding=1.45):.3f}"
'0.066'
"""
materials = MaterialSelection() if materials is None else materials
base_mode = _design_capillary(
radius,
gas,
temperature=temperature,
m=m,
model=model,
cladding=cladding,
lattice=None,
materials=materials,
)
omega = np.array([_omega(zdw_wavelength)])
def residual(pressure: float) -> float:
"""Return the dispersion at the wavelength for one pressure.
Parameters
----------
pressure : float
Pressure in bar.
Returns
-------
float
``beta_2`` in s^2/m, scaled to order one.
"""
mode = _capillary_at(base_mode, pressure)
# Scaled so that brentq's absolute tolerance on the residual is meaningless
# and its tolerance on the pressure decides, as intended.
return float(base.beta2(mode, omega)[0]) * 1e30
return _pressure_root(
residual,
pressure_max,
f"a zero-dispersion wavelength of {zdw_wavelength * 1e9:.4g} nm",
)