588 lines
23 KiB
Python
588 lines
23 KiB
Python
"""
|
|
Low-level mode-matching helpers for waveguide / EME workflows.
|
|
|
|
These helpers operate on already-solved and already-normalized port fields.
|
|
They do not build geometries or solve modes themselves; downstream users are
|
|
expected to supply compatible `(E, H)` modal field pairs from
|
|
`waveguide_2d`, `waveguide_3d`, or `waveguide_cyl`.
|
|
|
|
The returned matrices follow the usual port ordering:
|
|
|
|
- `get_tr(...)` returns `(T, R)` for left-incident modes.
|
|
- `get_abcd(...)` returns the 2-port block transfer matrix built from the two
|
|
directional `T/R` solves.
|
|
- `get_s(...)` returns the full block scattering matrix
|
|
`[[R12, T12], [T21, R21]]`.
|
|
- `get_taper_abcd(...)` and `get_taper_s(...)` build the same transfer /
|
|
scattering objects for sampled continuously varying sections using local-mode
|
|
CMT.
|
|
|
|
This module is intentionally a thin library layer rather than an integrated
|
|
simulation suite. It provides the overlap algebra that downstream users can
|
|
compose into larger workflows.
|
|
"""
|
|
|
|
from collections.abc import Sequence
|
|
import dataclasses
|
|
import numpy
|
|
from numpy.typing import NDArray
|
|
from scipy import linalg
|
|
from scipy import sparse
|
|
|
|
from ..fdmath import dx_lists2_t, vcfdfield2
|
|
from .waveguide_2d import inner_product
|
|
|
|
type wavenumber_seq = Sequence[complex] | NDArray[numpy.complexfloating] | NDArray[numpy.floating]
|
|
type mode_seq = Sequence[Sequence[vcfdfield2]]
|
|
|
|
|
|
@dataclasses.dataclass(frozen=True)
|
|
class ModeSection:
|
|
"""
|
|
Local modal basis at one longitudinal sample of a continuously varying guide.
|
|
|
|
Args:
|
|
z: Longitudinal coordinate of this section.
|
|
modes: Forward modes as `(E, H)` field pairs.
|
|
wavenumbers: Forward propagation constants for `modes`.
|
|
backward_modes: Optional explicit backward modes. If omitted, backward
|
|
modes are synthesized as `(E, -H)`.
|
|
backward_wavenumbers: Optional propagation constants for
|
|
`backward_modes`. If omitted, they are synthesized as `-wavenumbers`.
|
|
dual_modes: Optional forward dual / adjoint projection modes. If
|
|
omitted, `modes` are used as their own projection basis.
|
|
dual_backward_modes: Optional backward dual / adjoint projection modes.
|
|
If omitted, they are synthesized from `dual_modes` when available,
|
|
otherwise from `backward_modes`.
|
|
"""
|
|
|
|
z: float
|
|
modes: mode_seq
|
|
wavenumbers: wavenumber_seq
|
|
backward_modes: mode_seq | None = None
|
|
backward_wavenumbers: wavenumber_seq | None = None
|
|
dual_modes: mode_seq | None = None
|
|
dual_backward_modes: mode_seq | None = None
|
|
|
|
|
|
def _validate_port_modes(
|
|
name: str,
|
|
ehs: mode_seq,
|
|
wavenumbers: wavenumber_seq,
|
|
) -> tuple[tuple[int, ...], tuple[int, ...]]:
|
|
if len(ehs) != len(wavenumbers):
|
|
raise ValueError(f'{name} mode list and wavenumber list must have the same length')
|
|
if not ehs:
|
|
raise ValueError(f'{name} must contain at least one mode')
|
|
|
|
e_shape: tuple[int, ...] | None = None
|
|
h_shape: tuple[int, ...] | None = None
|
|
for index, mode in enumerate(ehs):
|
|
if len(mode) != 2:
|
|
raise ValueError(f'{name}[{index}] must be a 2-tuple of (E, H) modal fields')
|
|
e_field, h_field = mode
|
|
mode_e_shape = numpy.shape(e_field)
|
|
mode_h_shape = numpy.shape(h_field)
|
|
if mode_e_shape != mode_h_shape:
|
|
raise ValueError(f'{name}[{index}] has mismatched E/H field shapes')
|
|
if e_shape is None:
|
|
e_shape = mode_e_shape
|
|
h_shape = mode_h_shape
|
|
elif mode_e_shape != e_shape or mode_h_shape != h_shape:
|
|
raise ValueError(f'{name} modal fields must all share the same shape')
|
|
|
|
assert e_shape is not None
|
|
assert h_shape is not None
|
|
return e_shape, h_shape
|
|
|
|
|
|
def _validate_dual_modes(
|
|
name: str,
|
|
dual_ehs: mode_seq | None,
|
|
reference_shape: tuple[int, ...],
|
|
wavenumbers: wavenumber_seq,
|
|
) -> mode_seq | None:
|
|
if dual_ehs is None:
|
|
return None
|
|
|
|
dual_e_shape, dual_h_shape = _validate_port_modes(name, dual_ehs, wavenumbers)
|
|
if dual_e_shape != reference_shape or dual_h_shape != reference_shape:
|
|
raise ValueError(f'{name} modal fields must share the same E/H shapes as the corresponding modes')
|
|
return dual_ehs
|
|
|
|
|
|
def _as_wavenumber_array(
|
|
name: str,
|
|
wavenumbers: wavenumber_seq,
|
|
) -> NDArray[numpy.complex128]:
|
|
array = numpy.asarray(wavenumbers, dtype=complex)
|
|
if array.ndim != 1:
|
|
raise ValueError(f'{name} must be a one-dimensional sequence')
|
|
if not numpy.isfinite(array).all():
|
|
raise ValueError(f'{name} must contain only finite values')
|
|
return array
|
|
|
|
|
|
def _as_mode_arrays(
|
|
ehs: mode_seq,
|
|
) -> list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]]:
|
|
return [
|
|
(numpy.asarray(e_field, dtype=complex), numpy.asarray(h_field, dtype=complex))
|
|
for e_field, h_field in ehs
|
|
]
|
|
|
|
|
|
def _lorentz_overlap(
|
|
mode_a: tuple[vcfdfield2, vcfdfield2],
|
|
mode_b: tuple[vcfdfield2, vcfdfield2],
|
|
dxes: dx_lists2_t,
|
|
) -> complex:
|
|
e_a, h_a = mode_a
|
|
e_b, h_b = mode_b
|
|
return 0.5 * (
|
|
inner_product(e_a, h_b, dxes=dxes, conj_h=False)
|
|
+ inner_product(e_b, h_a, dxes=dxes, conj_h=False)
|
|
)
|
|
|
|
|
|
def _lorentz_derivative_overlap(
|
|
mode_a: tuple[vcfdfield2, vcfdfield2],
|
|
derivative_b: tuple[vcfdfield2, vcfdfield2],
|
|
dxes: dx_lists2_t,
|
|
) -> complex:
|
|
e_a, h_a = mode_a
|
|
de_b, dh_b = derivative_b
|
|
return 0.5 * (
|
|
inner_product(e_a, dh_b, dxes=dxes, conj_h=False)
|
|
+ inner_product(de_b, h_a, dxes=dxes, conj_h=False)
|
|
)
|
|
|
|
|
|
def _phase_align_modes(
|
|
previous: Sequence[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]],
|
|
current: Sequence[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]],
|
|
dxes: dx_lists2_t,
|
|
previous_dual: Sequence[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]] | None = None,
|
|
) -> list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]]:
|
|
aligned = []
|
|
test_modes = previous if previous_dual is None else previous_dual
|
|
for index, (previous_mode, current_mode, test_mode) in enumerate(zip(previous, current, test_modes, strict=True)):
|
|
overlap = _lorentz_overlap(test_mode, current_mode, dxes)
|
|
self_overlap = _lorentz_overlap(test_mode, previous_mode, dxes)
|
|
if overlap == 0:
|
|
raise ValueError(f'cannot phase-track mode {index}: adjacent section overlap is zero')
|
|
if self_overlap == 0:
|
|
raise ValueError(f'cannot phase-track mode {index}: mode dual-overlap is zero')
|
|
phase = (overlap / abs(overlap)) / (self_overlap / abs(self_overlap))
|
|
e_field, h_field = current_mode
|
|
aligned.append((e_field / phase, h_field / phase))
|
|
return aligned
|
|
|
|
|
|
def _section_branches(
|
|
section: ModeSection,
|
|
index: int,
|
|
expected_count: int | None,
|
|
expected_shape: tuple[int, ...] | None,
|
|
) -> tuple[
|
|
float,
|
|
list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]],
|
|
NDArray[numpy.complex128],
|
|
list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]],
|
|
NDArray[numpy.complex128],
|
|
list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]],
|
|
list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]],
|
|
tuple[int, ...],
|
|
]:
|
|
z_coord = float(section.z)
|
|
if not numpy.isfinite(z_coord):
|
|
raise ValueError(f'sections[{index}].z must be finite')
|
|
|
|
shape, _h_shape = _validate_port_modes(f'sections[{index}].modes', section.modes, section.wavenumbers)
|
|
wavenumbers = _as_wavenumber_array(f'sections[{index}].wavenumbers', section.wavenumbers)
|
|
if expected_count is not None and len(wavenumbers) != expected_count:
|
|
raise ValueError('all taper sections must contain the same number of modes')
|
|
if expected_shape is not None and shape != expected_shape:
|
|
raise ValueError('all taper section modal fields must share the same E/H shapes')
|
|
|
|
if (section.backward_modes is None) != (section.backward_wavenumbers is None):
|
|
raise ValueError('backward_modes and backward_wavenumbers must be supplied together')
|
|
|
|
forward_modes = _as_mode_arrays(section.modes)
|
|
if section.backward_modes is None:
|
|
backward_modes = [(e_field.copy(), -h_field.copy()) for e_field, h_field in forward_modes]
|
|
backward_wavenumbers = -wavenumbers
|
|
else:
|
|
backward_shape, _backward_h_shape = _validate_port_modes(
|
|
f'sections[{index}].backward_modes',
|
|
section.backward_modes,
|
|
section.backward_wavenumbers,
|
|
)
|
|
if backward_shape != shape:
|
|
raise ValueError('forward and backward modal fields must share the same E/H shapes')
|
|
backward_wavenumbers = _as_wavenumber_array(
|
|
f'sections[{index}].backward_wavenumbers',
|
|
section.backward_wavenumbers,
|
|
)
|
|
backward_modes = _as_mode_arrays(section.backward_modes)
|
|
|
|
if len(backward_wavenumbers) != len(wavenumbers):
|
|
raise ValueError('forward and backward mode counts must match')
|
|
|
|
if section.dual_modes is None:
|
|
dual_modes = forward_modes
|
|
else:
|
|
dual_shape, _dual_h_shape = _validate_port_modes(
|
|
f'sections[{index}].dual_modes',
|
|
section.dual_modes,
|
|
section.wavenumbers,
|
|
)
|
|
if dual_shape != shape:
|
|
raise ValueError('modal fields and dual modal fields must share the same E/H shapes')
|
|
dual_modes = _as_mode_arrays(section.dual_modes)
|
|
|
|
if section.dual_backward_modes is None:
|
|
if section.dual_modes is None and section.backward_modes is not None:
|
|
dual_backward_modes = backward_modes
|
|
else:
|
|
dual_backward_modes = [(e_field.copy(), -h_field.copy()) for e_field, h_field in dual_modes]
|
|
else:
|
|
dual_backward_shape, _dual_backward_h_shape = _validate_port_modes(
|
|
f'sections[{index}].dual_backward_modes',
|
|
section.dual_backward_modes,
|
|
section.backward_wavenumbers if section.backward_wavenumbers is not None else backward_wavenumbers,
|
|
)
|
|
if dual_backward_shape != shape:
|
|
raise ValueError('backward modal fields and dual backward modal fields must share the same E/H shapes')
|
|
dual_backward_modes = _as_mode_arrays(section.dual_backward_modes)
|
|
|
|
if len(dual_modes) != len(forward_modes) or len(dual_backward_modes) != len(backward_modes):
|
|
raise ValueError('dual mode counts must match the corresponding modal basis counts')
|
|
|
|
return z_coord, forward_modes, wavenumbers, backward_modes, backward_wavenumbers, dual_modes, dual_backward_modes, shape
|
|
|
|
|
|
def _validate_taper_sections(
|
|
sections: Sequence[ModeSection],
|
|
dxes: dx_lists2_t,
|
|
) -> tuple[
|
|
NDArray[numpy.float64],
|
|
list[list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]]],
|
|
list[list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]]],
|
|
list[NDArray[numpy.complex128]],
|
|
int,
|
|
]:
|
|
if len(sections) < 2:
|
|
raise ValueError('at least two taper sections are required')
|
|
|
|
z_coords: list[float] = []
|
|
branch_modes: list[list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]]] = []
|
|
branch_dual_modes: list[list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]]] = []
|
|
branch_wavenumbers: list[NDArray[numpy.complex128]] = []
|
|
explicit_duals: list[bool] = []
|
|
expected_count: int | None = None
|
|
expected_shape: tuple[int, ...] | None = None
|
|
|
|
for index, section in enumerate(sections):
|
|
z_coord, forward_modes, forward_wavenumbers, backward_modes, backward_wavenumbers, dual_modes, dual_backward_modes, shape = _section_branches(
|
|
section,
|
|
index,
|
|
expected_count,
|
|
expected_shape,
|
|
)
|
|
if expected_count is None:
|
|
expected_count = len(forward_wavenumbers)
|
|
expected_shape = shape
|
|
z_coords.append(z_coord)
|
|
branch_modes.append([*forward_modes, *backward_modes])
|
|
branch_dual_modes.append([*dual_modes, *dual_backward_modes])
|
|
branch_wavenumbers.append(numpy.concatenate((forward_wavenumbers, backward_wavenumbers)))
|
|
explicit_duals.append(section.dual_modes is not None or section.dual_backward_modes is not None)
|
|
|
|
z_array = numpy.asarray(z_coords, dtype=float)
|
|
if not (numpy.diff(z_array) > 0).all():
|
|
raise ValueError('taper section z coordinates must be strictly increasing')
|
|
|
|
for index in range(1, len(branch_modes)):
|
|
branch_modes[index] = _phase_align_modes(branch_modes[index - 1], branch_modes[index], dxes, branch_dual_modes[index - 1])
|
|
if not explicit_duals[index]:
|
|
branch_dual_modes[index] = branch_modes[index]
|
|
|
|
assert expected_count is not None
|
|
return z_array, branch_modes, branch_dual_modes, branch_wavenumbers, expected_count
|
|
|
|
|
|
def _taper_interval_generator(
|
|
left_modes: Sequence[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]],
|
|
left_dual_modes: Sequence[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]],
|
|
right_modes: Sequence[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]],
|
|
left_wavenumbers: NDArray[numpy.complex128],
|
|
right_wavenumbers: NDArray[numpy.complex128],
|
|
dz: float,
|
|
dxes: dx_lists2_t,
|
|
) -> NDArray[numpy.complex128]:
|
|
mode_count = len(left_modes)
|
|
gram = numpy.zeros((mode_count, mode_count), dtype=complex)
|
|
derivative_overlap = numpy.zeros((mode_count, mode_count), dtype=complex)
|
|
|
|
for row, left_row_mode in enumerate(left_dual_modes):
|
|
for col, left_col_mode in enumerate(left_modes):
|
|
gram[row, col] = _lorentz_overlap(left_row_mode, left_col_mode, dxes)
|
|
for col, (left_col_mode, right_col_mode) in enumerate(zip(left_modes, right_modes, strict=True)):
|
|
derivative = (
|
|
(right_col_mode[0] - left_col_mode[0]) / dz,
|
|
(right_col_mode[1] - left_col_mode[1]) / dz,
|
|
)
|
|
derivative_overlap[row, col] = _lorentz_derivative_overlap(left_row_mode, derivative, dxes)
|
|
|
|
coupling = numpy.linalg.pinv(gram) @ derivative_overlap
|
|
propagation = numpy.diag(-1j * 0.5 * (left_wavenumbers + right_wavenumbers))
|
|
return propagation - coupling
|
|
|
|
|
|
def _abcd_to_s(
|
|
abcd: NDArray[numpy.complex128],
|
|
n_modes: int,
|
|
) -> NDArray[numpy.complex128]:
|
|
A = abcd[:n_modes, :n_modes]
|
|
B = abcd[:n_modes, n_modes:]
|
|
C = abcd[n_modes:, :n_modes]
|
|
D = abcd[n_modes:, n_modes:]
|
|
D_inv = numpy.linalg.pinv(D)
|
|
r12 = -D_inv @ C
|
|
t21 = D_inv
|
|
t12 = A - B @ D_inv @ C
|
|
r21 = B @ D_inv
|
|
return numpy.block([[r12, t12],
|
|
[t21, r21]])
|
|
|
|
|
|
def get_tr(
|
|
ehLs: mode_seq,
|
|
wavenumbers_L: wavenumber_seq,
|
|
ehRs: mode_seq,
|
|
wavenumbers_R: wavenumber_seq,
|
|
dxes: dx_lists2_t,
|
|
dual_ehLs: mode_seq | None = None,
|
|
) -> tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]:
|
|
"""
|
|
Compute left-incident transmission and reflection matrices.
|
|
|
|
Args:
|
|
ehLs: Left-port modes as `(E, H)` field pairs.
|
|
wavenumbers_L: Propagation constants for `ehLs`.
|
|
ehRs: Right-port modes as `(E, H)` field pairs.
|
|
wavenumbers_R: Propagation constants for `ehRs`.
|
|
dxes: Two-dimensional Yee-cell edge lengths for the shared port plane.
|
|
dual_ehLs: Optional left-port dual / adjoint projection modes. If
|
|
omitted, `ehLs` are used as their own projection basis.
|
|
|
|
Returns:
|
|
`(T12, R12)` where columns index left-incident modes and rows index
|
|
outgoing right-going / left-going modes respectively.
|
|
|
|
Raises:
|
|
ValueError: If the port mode lists are empty, malformed, or defined on
|
|
incompatible field shapes.
|
|
"""
|
|
left_e_shape, left_h_shape = _validate_port_modes('ehLs', ehLs, wavenumbers_L)
|
|
right_e_shape, right_h_shape = _validate_port_modes('ehRs', ehRs, wavenumbers_R)
|
|
if left_e_shape != right_e_shape or left_h_shape != right_h_shape:
|
|
raise ValueError('left and right modal fields must share the same E/H shapes')
|
|
dual_projection_ehLs = _validate_dual_modes('dual_ehLs', dual_ehLs, left_e_shape, wavenumbers_L)
|
|
projection_ehLs = ehLs if dual_projection_ehLs is None else dual_projection_ehLs
|
|
|
|
nL = len(wavenumbers_L)
|
|
nR = len(wavenumbers_R)
|
|
A12 = numpy.zeros((nL, nR), dtype=complex)
|
|
A21 = numpy.zeros((nL, nR), dtype=complex)
|
|
B11 = numpy.zeros((nL,), dtype=complex)
|
|
for ll in range(nL):
|
|
eL, hL = ehLs[ll]
|
|
eP, hP = projection_ehLs[ll]
|
|
B11[ll] = inner_product(eL, hP, dxes=dxes, conj_h=False)
|
|
for rr in range(nR):
|
|
eR, hR = ehRs[rr]
|
|
A12[ll, rr] = inner_product(eP, hR, dxes=dxes, conj_h=False) # TODO optimize loop?
|
|
A21[ll, rr] = inner_product(eR, hP, dxes=dxes, conj_h=False)
|
|
|
|
# tt0 = 2 * numpy.linalg.pinv(A21 + numpy.conj(A12))
|
|
tt0, _resid, _rank, _sing = numpy.linalg.lstsq(A21 + A12, numpy.diag(2 * B11), rcond=None)
|
|
|
|
U, st, V = numpy.linalg.svd(tt0)
|
|
gain = st > 1
|
|
st[gain] = 1 / st[gain]
|
|
tt = U @ numpy.diag(st) @ V
|
|
|
|
# rr = 0.5 * (A21 - numpy.conj(A12)) @ tt
|
|
rr = numpy.diag(0.5 / B11) @ (A21 - A12) @ tt
|
|
|
|
return tt, rr
|
|
|
|
|
|
def get_abcd(
|
|
ehLs: mode_seq,
|
|
wavenumbers_L: wavenumber_seq,
|
|
ehRs: mode_seq,
|
|
wavenumbers_R: wavenumber_seq,
|
|
dual_ehLs: mode_seq | None = None,
|
|
dual_ehRs: mode_seq | None = None,
|
|
**kwargs,
|
|
) -> sparse.sparray:
|
|
"""
|
|
Build the 2-port block transfer matrix for an interface.
|
|
|
|
The blocks are assembled from the forward and reverse `get_tr(...)`
|
|
solutions using the standard
|
|
|
|
`[[A, B], [C, D]] = [[T12 - R21 T21^-1 R12, R21 T21^-1], [-T21^-1 R12, T21^-1]]`
|
|
|
|
convention.
|
|
"""
|
|
t12, r12 = get_tr(ehLs, wavenumbers_L, ehRs, wavenumbers_R, dual_ehLs=dual_ehLs, **kwargs)
|
|
t21, r21 = get_tr(ehRs, wavenumbers_R, ehLs, wavenumbers_L, dual_ehLs=dual_ehRs, **kwargs)
|
|
t21i = numpy.linalg.pinv(t21)
|
|
A = t12 - r21 @ t21i @ r12
|
|
B = r21 @ t21i
|
|
C = -t21i @ r12
|
|
D = t21i
|
|
return sparse.block_array(
|
|
[
|
|
[sparse.csr_array(A), sparse.csr_array(B)],
|
|
[sparse.csr_array(C), sparse.csr_array(D)],
|
|
],
|
|
format='csr',
|
|
)
|
|
|
|
|
|
def get_s(
|
|
ehLs: mode_seq,
|
|
wavenumbers_L: wavenumber_seq,
|
|
ehRs: mode_seq,
|
|
wavenumbers_R: wavenumber_seq,
|
|
force_nogain: bool = False,
|
|
force_reciprocal: bool = False,
|
|
dual_ehLs: mode_seq | None = None,
|
|
dual_ehRs: mode_seq | None = None,
|
|
**kwargs,
|
|
) -> NDArray[numpy.complex128]:
|
|
"""
|
|
Build the full block scattering matrix for a two-sided interface.
|
|
|
|
The returned matrix is ordered as `[[R12, T12], [T21, R21]]`, where the
|
|
first block-row/column corresponds to the left port and the second to the
|
|
right port.
|
|
|
|
Args:
|
|
force_nogain: If `True`, clamp singular values of the assembled
|
|
scattering matrix to at most one.
|
|
force_reciprocal: If `True`, symmetrize the assembled matrix as
|
|
`0.5 * (S + S.T)`.
|
|
dual_ehLs: Optional left-port dual / adjoint projection modes.
|
|
dual_ehRs: Optional right-port dual / adjoint projection modes.
|
|
"""
|
|
t12, r12 = get_tr(ehLs, wavenumbers_L, ehRs, wavenumbers_R, dual_ehLs=dual_ehLs, **kwargs)
|
|
t21, r21 = get_tr(ehRs, wavenumbers_R, ehLs, wavenumbers_L, dual_ehLs=dual_ehRs, **kwargs)
|
|
|
|
ss = numpy.block([[r12, t12],
|
|
[t21, r21]])
|
|
|
|
if force_nogain:
|
|
# force S @ S.H diagonal
|
|
U, sing, Vh = numpy.linalg.svd(ss)
|
|
ss = U @ numpy.diag(numpy.minimum(sing, 1.0)) @ Vh
|
|
|
|
if force_reciprocal:
|
|
ss = 0.5 * (ss + ss.T)
|
|
|
|
return ss
|
|
|
|
|
|
def get_taper_abcd(
|
|
sections: Sequence[ModeSection],
|
|
dxes: dx_lists2_t,
|
|
*,
|
|
rtol: float = 1e-7,
|
|
atol: float = 1e-9,
|
|
max_step: float | None = None,
|
|
) -> sparse.sparray:
|
|
"""
|
|
Build a bidirectional transfer matrix for a continuously varying taper.
|
|
|
|
The taper is represented by local modal bases sampled at increasing `z`
|
|
coordinates. Adjacent modal phases are tracked with the same non-conjugated
|
|
Lorentz/Poynting overlap used by the abrupt-interface helpers, then each
|
|
interval is propagated with a finite-difference local-mode CMT generator.
|
|
If a `ModeSection` supplies dual / adjoint modes, those modes are used for
|
|
the CMT projection. This supports leaky or radiative mode sets whose natural
|
|
projection basis is biorthogonal rather than self-projected.
|
|
|
|
Args:
|
|
sections: Local modal samples ordered by increasing `z`.
|
|
dxes: Two-dimensional Yee-cell edge lengths for the shared port plane.
|
|
rtol: Relative tolerance reserved for future adaptive CMT integrators.
|
|
Must be positive.
|
|
atol: Absolute tolerance reserved for future adaptive CMT integrators.
|
|
Must be positive.
|
|
max_step: Optional maximum matrix-exponential step inside each sampled
|
|
interval. This does not change the piecewise-constant interval
|
|
generator, but can improve conditioning for long intervals.
|
|
|
|
Returns:
|
|
Sparse block transfer matrix ordered as `[forward, backward]`.
|
|
"""
|
|
if rtol <= 0:
|
|
raise ValueError('rtol must be positive')
|
|
if atol <= 0:
|
|
raise ValueError('atol must be positive')
|
|
if max_step is not None and max_step <= 0:
|
|
raise ValueError('max_step must be positive')
|
|
|
|
z_coords, branch_modes, branch_dual_modes, branch_wavenumbers, n_modes = _validate_taper_sections(sections, dxes)
|
|
branch_count = 2 * n_modes
|
|
transfer = numpy.eye(branch_count, dtype=complex)
|
|
|
|
for index, dz in enumerate(numpy.diff(z_coords)):
|
|
generator = _taper_interval_generator(
|
|
branch_modes[index],
|
|
branch_dual_modes[index],
|
|
branch_modes[index + 1],
|
|
branch_wavenumbers[index],
|
|
branch_wavenumbers[index + 1],
|
|
float(dz),
|
|
dxes,
|
|
)
|
|
step_count = 1 if max_step is None else max(1, int(numpy.ceil(dz / max_step)))
|
|
interval_transfer = linalg.expm(generator * (dz / step_count))
|
|
for _step in range(step_count):
|
|
transfer = interval_transfer @ transfer
|
|
|
|
return sparse.csr_array(transfer)
|
|
|
|
|
|
def get_taper_s(
|
|
sections: Sequence[ModeSection],
|
|
dxes: dx_lists2_t,
|
|
*,
|
|
force_nogain: bool = False,
|
|
force_reciprocal: bool = False,
|
|
**kwargs,
|
|
) -> NDArray[numpy.complex128]:
|
|
"""
|
|
Build the full block scattering matrix for a continuously varying taper.
|
|
|
|
The returned matrix uses the same ordering as `get_s(...)`:
|
|
`[[R12, T12], [T21, R21]]`.
|
|
"""
|
|
_z_coords, _branch_modes, _branch_dual_modes, _branch_wavenumbers, n_modes = _validate_taper_sections(sections, dxes)
|
|
abcd = get_taper_abcd(sections, dxes, **kwargs).toarray()
|
|
ss = _abcd_to_s(abcd, n_modes)
|
|
|
|
if force_nogain:
|
|
U, sing, Vh = numpy.linalg.svd(ss)
|
|
ss = U @ numpy.diag(numpy.minimum(sing, 1.0)) @ Vh
|
|
|
|
if force_reciprocal:
|
|
ss = 0.5 * (ss + ss.T)
|
|
|
|
return ss
|