add links to pypi and github

README.md

@@ -41,14 +41,15 @@ solver will need the ability to solve complex symmetric (non-Hermitian)
 linear systems, ideally with double precision.
 
 - [Source repository](https://mpxd.net/code/jan/meanas)
-- PyPI *TBD*
+- [PyPI](https://pypi.org/project/meanas)
+- [Github mirror](https://github.com/anewusername/meanas)
 
 
 ## Installation
 
 **Requirements:**
 
-* python >=3.8
+* python >=3.11
 * numpy
 * scipy
 

nom-eme.py

@@ -0,0 +1,391 @@
+import scipy
+import numpy
+from numpy.typing import ArrayLike, NDArray
+
+
+#from simphony.elements import Model
+#from simphony.netlist import Subcircuit
+#from simphony.simulation import SweepSimulation
+#
+#from matplotlib import pyplot as plt
+#
+#
+#class PeriodicLayer(Model):
+#    def __init__(self, left_modes, right_modes, s_params):
+#        self.left_modes = left_modes
+#        self.right_modes = right_modes
+#        self.left_ports = len(self.left_modes)
+#        self.right_ports = len(self.right_modes)
+#        self.normalize_fields()
+#        self.s_params = s_params
+#
+#    def normalize_fields(self):
+#        for mode in range(len(self.left_modes)):
+#            self.left_modes[mode].normalize()
+#        for mode in range(len(self.right_modes)):
+#            self.right_modes[mode].normalize()
+#
+#
+#class PeriodicEME:
+#    def __init__(self, layers=[], num_periods=1):
+#        self.layers = layers
+#        self.num_periods = num_periods
+#        self.wavelength = wavelength
+#
+#    def propagate(self):
+#        wl = self.wavelength
+#        if not len(self.layers):
+#            raise Exception("Must place layers before propagating")
+#
+#        num_modes = max([l.num_modes for l in self.layers])
+#        iface = InterfaceSingleMode if num_modes == 1 else InterfaceMultiMode
+#
+#        eme = EME(layers=self.layers)
+#        left, right = eme.propagate()
+#        self.single_period = eme.s_matrix
+#
+#        period_layer = PeriodicLayer(left.modes, right.modes, self.single_period)
+#        current_layer = PeriodicLayer(left.modes, right.modes, self.single_period)
+#        interface = iface(right, left)
+#
+#        for _ in range(self.num_periods - 1):
+#            current_layer.s_params = cascade(current_layer, interface, wl)
+#            current_layer.s_params = cascade(current_layer, period_layer, wl)
+#
+#        self.s_params = current_layer.s_params
+#
+#
+#class EME:
+#    def __init__(self, layers=[]):
+#        self.layers = layers
+#        self.wavelength = None
+#
+#    def propagate(self):
+#        layers = self.layers
+#        wl = layers[0].wavelength if self.wavelength is None else self.wavelength
+#        if not len(layers):
+#            raise Exception("Must place layers before propagating")
+#
+#        num_modes = max([l.num_modes for l in layers])
+#        iface = InterfaceSingleMode if num_modes == 1 else InterfaceMultiMode
+#
+#        first_layer = layers[0]
+#        current = Current(wl, first_layer)
+#        interface = iface(first_layer, layers[1])
+#
+#        current.s = cascade(current, interface, wl)
+#        current.right_pins = interface.right_pins
+#
+#        for index in range(1, len(layers) - 1):
+#            layer1 = layers[index]
+#            layer2 = layers[index + 1]
+#            interface = iface(layer1, layer2)
+#
+#            current.s = cascade(current, layer1, wl)
+#            current.right_pins = layer1.right_pins
+#
+#            current.s = cascade(current, interface, wl)
+#            current.right_pins = interface.right_pins
+#
+#        last_layer = layers[-1]
+#        current.s = cascade(current, last_layer, wl)
+#        current.right_pins = last_layer.right_pins
+#
+#        self.s_matrix = current.s
+#        return first_layer, last_layer
+#
+#
+#def stack(sa, sb):
+#    qab = numpy.eye() - sa.r11 @ sb.r11
+#    qba = numpy.eye() - sa.r11 @ sb.r11
+#    #s.t12 = sa.t12 @ numpy.pinv(qab) @ sb.t12
+#    #s.r21 = sa.t12 @ numpy.pinv(qab) @ sb.r22 @ sa.t21 + sa.r22
+#    #s.r12 = sb.t21 @ numpy.pinv(qba) @ sa.r11 @ sb.t12 + sb.r11
+#    #s.t21 = sb.t21 @ numpy.pinv(qba) @ sa.t21
+#    s.t12 = sa.t12 @ numpy.linalg.solve(qab, sb.t12)
+#    s.r21 = sa.t12 @ numpy.linalg.solve(qab, sb.r22 @ sa.t21) + sa.r22
+#    s.r12 = sb.t21 @ numpy.linalg.solve(qba, sa.r11 @ sb.t12) + sb.r11
+#    s.t21 = sb.t21 @ numpy.linalg.solve(qba, sa.t21)
+#    return s
+#
+#
+#def cascade(first, second, wavelength):
+#    circuit = Subcircuit("Device")
+#
+#    circuit.add([(first, "first"), (second, "second")])
+#    for port in range(first.right_ports):
+#        circuit.connect("first", "right" + str(port), "second", "left" + str(port))
+#
+#    simulation = SweepSimulation(circuit, wavelength, wavelength, num=1)
+#    result = simulation.simulate()
+#    return result.s
+#
+#
+#class InterfaceSingleMode(Model):
+#    def __init__(self, layer1, layer2, num_modes=1):
+#        self.num_modes = num_modes
+#        self.num_ports = 2 * num_modes
+#        self.s = self.solve(layer1, layer2, num_modes)
+#
+#    def solve(self, layer1, layer2, num_modes):
+#        nm = num_modes
+#        s = numpy.zeros((2 * nm, 2 * nm), dtype=complex)
+#
+#        for ii, left_mode in enumerate(layer1.modes):
+#            for oo, right_mode in enumerate(layer2.modes):
+#                r, t = get_rt(left_mode, right_mode)
+#                s[     oo, ii] = r
+#                s[nm + oo, ii] = t
+#
+#        for ii, right_mode in enumerate(layer2.modes):
+#            for oo, left_mode in enumerate(layer1.modes):
+#                r, t = get_rt(right_mode, left_mode)
+#                s[     oo, nm + ii] = t
+#                s[nm + oo, nm + ii] = r
+#        return s
+#
+#
+#class InterfaceMultiMode(Model):
+#    def __init__(self, layer1, layer2):
+#        self.s = self.solve(layer1, layer2)
+#
+#    def solve(self, layer1, layer2):
+#        n1p = layer1.num_modes
+#        n2p = layer2.num_modes
+#        num_ports = n1p + n2p
+#        s = numpy.zeros((num_ports, num_ports), dtype=complex)
+#
+#        for l1p in range(n1p):
+#            ts = get_t(l1p, layer1, layer2)
+#            rs = get_r(l1p, layer1, layer2, ts)
+#            s[n1p:, l1p] = ts
+#            s[:n1p, l1p] = rs
+#
+#        for l2p in range(n2p):
+#            ts = get_t(l2p, layer2, layer1)
+#            rs = get_r(l2p, layer2, layer1, ts)
+#            s[:n1p, n1p + l2p] = ts
+#            s[n1p:, n1p + l2p] = rs
+#
+#        return s
+
+
+def get_t(p, left, right):
+    A = numpy.empty(left.num_modes, right.num_modes, dtype=complex)
+    for i in range(left.num_modes):
+        for k in range(right.num_modes):
+            # TODO optimize loop
+            A[i, k] = inner_product(right[k], left[i]) + inner_product(left[i], right[k])
+
+    b = numpy.zeros(left.num_modes)
+    b[p] = 2 * inner_product(left[p], left[p])
+
+    x = numpy.linalg.solve(A, b)
+    # NOTE: `A` does not depend on `p`, so it might make sense to partially precompute
+    # the solution (pinv(A), or LU decomposition?)
+    # Actually solve() can take multiple vectors, so just pass it something with the full diagonal?
+
+    xx = numpy.matmul(numpy.linalg.pinv(A), b)        #TODO verify
+    assert numpy.allclose(xx, x)
+    return x
+
+
+def get_r(p, left, right, t):
+    r = numpy.empty(left.num_modes, dtype=complex)
+    for ii in range(left.num_modes):
+        r[ii] = sum((inner_product(right[kk], left[ii]) - inner_product(left[ii], right[kk])) * t[kk]
+                    for kk in range(right.num_modes)
+                   ) / (2 * inner_product(left[ii], left[ii]))
+    return r
+
+
+def get_rt(left, right):
+    a = 0.5 * (inner_product(left, right) + inner_product(right, left))
+    b = 0.5 * (inner_product(left, right) - inner_product(right, left))
+    t = (a ** 2 - b ** 2) / a
+    r = 1 - t / (a + b)
+    return -r, t
+
+
+def inner_product(left_E, right_H, dxes):
+    #  ExHy' - EyHx'
+    cross_z = left_E[0] * right_H[1].conj() - left_E[1] * right_H[0].conj()
+#    cross_z = numpy.cross(left_E, numpy.conj(right_H), axisa=0, axisb=0, axisc=0)[2]
+    return numpy.trapz(numpy.trapz(cross_z, dxes[0][0]), dxes[0][1]) / 2       # TODO might need cumsum on dxes
+
+
+def propagation_matrix(mode_neffs: ArrayLike, wavelength: float, distance: float):
+    eigenv = numpy.array(mode_neffs, copy=False) * 2 * numpy.pi / wavelength
+    prop_diag = numpy.diag(numpy.exp(distance * 1j * numpy.hstack((eigenv, eigenv))))
+    prop_matrix = numpy.roll(prop_diag, len(eigenv), axis=0)
+    return prop_matrix
+
+
+def connect_s(
+        A: NDArray[numpy.complex128],
+        k: int,
+        B: NDArray[numpy.complex128],
+        l: int,
+        ) -> NDArray[numpy.complex128]:
+    """
+    TODO
+    freq x ...  x n x n
+
+    Based on skrf implementation
+
+    Connect two n-port networks' s-matrices together.
+
+    Specifically, connect port `k` on network `A` to port `l` on network
+    `B`. The resultant network has nports = (A.rank + B.rank-2); first
+    (A.rank - 1) ports are from `A`, remainder are from B.
+
+    Assumes same reference impedance for both `k` and `l`; may need to
+    connect an "impedance mismatch" thru element first!
+
+    Args:
+        A: S-parameter matrix of `A`, shape is fxnxn
+        k: port index on `A` (port indices start from 0)
+        B: S-parameter matrix of `B`, shape is fxnxn
+        l: port index on `B`
+
+    Returns:
+        new S-parameter matrix
+    """
+    if k > A.shape[-1] - 1 or l > B.shape[-1] - 1:
+        raise ValueError("port indices are out of range")
+
+    #C = scipy.sparse.block_diag((A, B), dtype=complex)
+    #return innerconnect_s(C, k, A.shape[0] + l)
+
+    nA = A.shape[-1]
+    nB = B.shape[-1]
+    nC = nA + nB - 2
+    assert numpy.array_equal(A.shape[:-2], B.shape[:-2])
+
+    ll = slice(l, l + 1)
+    kk = slice(k, k + 1)
+
+    denom = 1 - A[..., kk, kk] * B[..., ll, ll]
+    Anew = A + A[..., kk, :] * B[..., ll, ll] * A[..., :, kk] / denom
+    Bnew = A[..., kk, :] * B[..., :, ll] / denom
+    Anew = numpy.delete(Anew, (k,), 1)
+    Anew = numpy.delete(Anew, (k,), 2)
+    Bnew = numpy.delete(Bnew, (l,), 1)
+    Bnew = numpy.delete(Bnew, (l,), 2)
+
+    dtype = (A[0, 0] * B[0, 0]).dtype
+    C = numpy.zeros(tuple(A.shape[:-2]) + (nC, nC), dtype=dtype)
+    C[..., :nA - 1, :nA - 1] = Anew
+    C[..., nA - 1:, nA - 1:] = Bnew
+    return C
+
+
+def innerconnect_s(
+        S: NDArray[numpy.complex128],
+        k: int,
+        l: int,
+        ) -> NDArray[numpy.complex128]:
+    """
+    TODO
+    freq x ...  x n x n
+
+    Based on skrf implementation
+
+
+    Connect two ports of a single n-port network's s-matrix.
+    Specifically, connect port `k`  to port `l` on `S`. This results in
+    a (n-2)-port network.
+
+    Assumes same reference impedance for both `k` and `l`; may need to
+    connect an "impedance mismatch" thru element first!
+
+    Args:
+        S: S-parameter matrix of `S`, shape is fxnxn
+        k: port index on `S` (port indices start from 0)
+        l: port index on `S`
+
+    Returns:
+        new S-parameter matrix
+
+    Notes:
+        - Compton, R.C., "Perspectives in microwave circuit analysis",
+            doi:10.1109/MWSCAS.1989.101955
+        - Filipsson, G., "A New General Computer Algorithm for S-Matrix Calculation
+            of Interconnected Multiports",
+            doi:10.1109/EUMA.1981.332972
+    """
+    if k > S.shape[-1] - 1 or l > S.shape[-1] - 1:
+        raise ValueError("port indices are out of range")
+
+    ll = slice(l, l + 1)
+    kk = slice(k, k + 1)
+
+    mkl = 1 - S[..., kk, ll]
+    mlk = 1 - S[..., ll, kk]
+    C = S + (
+              S[..., kk, :] * S[..., :, l] * mlk
+            + S[..., ll, :] * S[..., :, k] * mkl
+            + S[..., kk, :] * S[..., l, l] * S[..., :, kk]
+            + S[..., ll, :] * S[..., k, k] * S[..., :, ll]
+        ) / (
+            mlk * mkl - S[..., kk, kk] * S[..., ll, ll]
+        )
+
+    # remove connected ports
+    C = numpy.delete(C, (k, l), 1)
+    C = numpy.delete(C, (k, l), 2)
+
+    return C
+
+
+def s2abcd(
+        S: NDArray[numpy.complex128],
+        z0: NDArray[numpy.complex128],
+        ) -> NDArray[numpy.complex128]:
+    assert numpy.array_equal(S.shape[:2] == (2, 2))
+    Z1, Z2 = z0
+    cross = S[0, 1] * S[1, 0]
+
+    T = numpy.empty_like(S, dtype=complex)
+    T[0, 0, :] = (Z1.conj() + S[0, 0] * Z1) * (1 - S[1, 1]) + cross * Z1    # A numerator
+    T[0, 1, :] = (Z1.conj() + S[0, 0] * Z1) * (Z1.conj() + S[1, 1] * Z2) - cross * Z1 * Z2  # B numerator
+    T[1, 0, :] = (1 - S[0, 0]) * (1 - S[1, 1]) - cross                      # C numerator
+    T[1, 1, :] = (1 - S[0, 0]) * (Z2.conj() + S[1, 1] * Z2) + cross * Z2    # D numerator
+    det = 2 * S[1, 0] * numpy.sqrt(Z1.real * Z2.real)
+    T /= det
+    return T
+
+
+def generalize_S(
+        S: NDArray[numpy.complex128],
+        r0: float,
+        z0: NDArray[numpy.complex128],
+        ) -> NDArray[numpy.complex128]:
+    g = (z0 - r0) / (z0 + r0)
+    D = numpy.diag((1 - g) / numpy.abs(1 - g.conj()) * numpy.sqrt(1 - numpy.abs(g * g.conj())))
+    G = numpy.diag(g)
+    U = numpy.eye(S.shape[-1]).reshape((S.ndim - 2) * (1,) + (S.shape[-2], S.shape[-1]))
+    S_gen = pinv(D.conj()) @ (S - G.conj()) @ pinv(U - G @ S) @ D
+    return S_gen
+
+
+def change_R0(
+        S: NDArray[numpy.complex128],
+        r1: float,
+        r2: float,
+        ) -> NDArray[numpy.complex128]:
+    g = (r2 - r1) / (r2 + r1)
+    U = numpy.eye(S.shape[-1]).reshape((S.ndim - 2) * (1,) + (S.shape[-2], S.shape[-1]))
+    G = U * g
+    S_r2 = (S - G) @ pinv(U - G @ S)
+    return S_r2
+
+# Zc = numpy.sqrt(B / C)
+# gamma = numpy.arccosh(A) / L_TL
+# n_eff = -1j * gamma * c_light / (2 * pi * f)
+# n_eff_grp = n_eff + f * diff(n_eff) / diff(f)
+# attenuation = (1 - S[0, 0] * S[0, 0].conj()) / (S[1, 0] * S[1, 0].conj())
+# R = numpy.real(gamma * Zc)
+# C = numpy.real(gamma / Zc)
+# L = numpy.imag(gamma * Zc) / (-1j * 2 * pi * f)
+# G = numpy.imag(gamma / Zc) / (-1j * 2 * pi * f)

pyproject.toml

@@ -33,7 +33,7 @@ classifiers = [
     "License :: OSI Approved :: GNU Affero General Public License v3",
     "Topic :: Scientific/Engineering :: Physics",
     ]
-requires-python = ">=3.8"
+requires-python = ">=3.11"
 include = [
     "LICENSE.md"
     ]

spar.py

@@ -0,0 +1,387 @@
+# Based on scripts from Andy H. va rfcafe
+
+# IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES. VOL 42, NO 2. FEBRUARY 1994
+# Conversions Between S, Z, Y, h, ABCD, and T Parameters which are Valid for Complex Source and Load Impedances
+# Dean A. Frickey, Member, EEE
+# Tables I and II
+
+import numpy
+
+
+def s_to_z(s, z0):
+    """
+    Scattering (S) to Impedance (Z)
+
+    Args:
+        s: The scattering matrix.
+        z0: The port impedances (Ohms).
+
+    Returns:
+        The impedance matrix.
+    """
+    z0c = numpy.conj(z0)
+
+    z = numpy.empty([2, 2], dtype=complex)
+    z[0, 0] = (z0c[0] + s[0, 0] * z0[0]) * (1 - s[1, 1]) + s[0, 1] * s[1, 0] * z0[0]
+    z[0, 1] = 2 * s[0, 1] * numpy.sqrt(z0[0].real * z0[1].real)
+    z[1, 0] = 2 * s[1, 0] * numpy.sqrt(z0[0].real * z0[1].real)
+    z[1, 1] = (1 - s[0, 0]) * (z0c[1] + s[1, 1] * z0[1]) + s[0, 1] * s[1, 0] * z0[1]
+
+    z /= (1 - s[0, 0]) * (1 - s[1, 1]) - s[0, 1] * s[1, 0]
+    return z
+
+
+def z_to_s(z, z0):
+    """
+    Impedance (Z) to Scattering (S)
+
+    Args:
+        z: The impedance matrix.
+        z0: The port impedances (Ohms).
+
+    Returns:
+        The scattering matrix.
+    """
+    z0c = numpy.conj(z0)
+
+    s = numpy.empty([2, 2], dtype=complex)
+    s[0, 0] = (z[0, 0] - z0c[0]) * (z[1, 1] + z0[1]) - z[0, 1] * z[1, 0]
+    s[0, 1] = 2 * z[0, 1] * numpy.sqrt(z0[0].real * z0[1].real)
+    s[1, 0] = 2 * z[1, 0] * numpy.sqrt(z0[0].real * z0[1].real)
+    s[1, 1] = (z[0, 0] + z0[0]) * (z[1, 1] - z0c[1]) - z[0, 1] * z[1, 0]
+
+    s /= (z[0, 0] + z0[0]) * (z[1, 1] + z0[1]) - z[0, 1] * z[1, 0]
+    return s
+
+
+def s_to_y(s, z0):
+    """
+    Scattering (S) to Admittance (Y)
+
+    Args:
+        s: The scattering matrix.
+        z0: The port impedances (Ohms).
+
+    Returns:
+        The admittance matrix.
+    """
+    z0c = numpy.conj(z0)
+
+    y = numpy.empty([2, 2], dtype=complex)
+    y[0, 0] = (1 - s[0, 0]) * (z0c[1] + s[1, 1] * z0[1]) + s[0,1] * s[1, 0] * z0[1]
+    y[0, 1] = -2 * s[0,1] * numpy.sqrt(z0[0].real * z0[1].real)
+    y[1, 0] = -2 * s[1, 0] * numpy.sqrt(z0[0].real * z0[1].real)
+    y[1, 1] = (z0c[0] + s[0, 0] * z0[0]) * (1 - s[1,1]) + s[0,1] * s[1, 0] * z0[0]
+
+    y /= (z0c[0] + s[0, 0] * z0[0]) * (z0c[1] + s[1, 1] * z0[1]) - s[0,1] * s[1, 0] * z0[0] * z0[1]
+    return y
+
+
+def y_to_s(y, z0):
+    """
+    Admittance (Y) to Scattering (S)
+
+    Args:
+        y: The admittance matrix.
+        z0: The port impedances (Ohms).
+
+    Returns:
+        The scattering matrix.
+    """
+    z0c = numpy.conj(z0)
+
+    s = numpy.empty([2, 2], dtype=complex)
+    s[0, 0] = (1 - y[0, 0] * z0c[0]) * (1 + y[1, 1] * z0[1]) + y[0,1] * y[1, 0] * z0c[0] * z0[1]
+    s[0, 1] = -2 * y[0,1] * numpy.sqrt(z0[0].real * z0[1].real)
+    s[1, 0] = -2 * y[1, 0] * numpy.sqrt(z0[0].real * z0[1].real)
+    s[1, 1] = (1 + y[0, 0] * z0[0]) * (1 - y[1,1] * z0c[1]) + y[0,1] * y[1, 0] * z0[0] * z0c[1]
+
+    s /= (1 + y[0, 0] * z0[0]) * (1 + y[1, 1] * z0[1]) - y[0,1] * y[1, 0] * z0[0] * z0[1]
+    return s
+
+
+def s_to_h(s, z0):
+    """
+    Scattering (S) to Hybrid (H)
+
+    Args:
+        s: The scattering matrix.
+        z0: The port impedances (Ohms).
+
+    Returns:
+        The hybrid matrix.
+    """
+    z0c = numpy.conj(z0)
+
+    h = numpy.empty([2, 2], dtype=complex)
+    h[0, 0] = (z0c[0] + s[0, 0] * z0[0]) * (z0c[1] + s[1, 1] * z0[1]) - s[0,1] * s[1, 0] * z0[0] * z0[1]
+    h[0, 1] = 2 * s[0,1] * numpy.sqrt(z0[0].real * z0[1].real)
+    h[1, 0] = -2 * s[1, 0] * numpy.sqrt(z0[0].real * z0[1].real)
+    h[1, 1] = (1 - s[0, 0]) * (1 - s[1,1]) - s[0,1] * s[1, 0]
+
+    h /= (1 - s[0, 0]) * (z0c[1] + s[1, 1] * z0[1]) + s[0,1] * s[1, 0] * z0[1]
+    return h
+
+
+def h_to_s(h, z0):
+    """
+    Hybrid (H) to Scattering (S)
+
+    Args:
+        h: The hybrid matrix.
+        z0: The port impedances (Ohms).
+
+    Returns:
+        The scattering matrix.
+    """
+    z0c = numpy.conj(z0)
+
+    s = numpy.empty([2, 2], dtype=complex)
+    s[0, 0] = (h[0, 0] - z0c[0]) * (1 + h[1, 1] * z0[1]) - h[0,1] * h[1, 0] * z0[1]
+    s[0, 1] = 2 * h[0,1] * numpy.sqrt(z0[0].real * z0[1].real)
+    s[1, 0] = -2 * h[1, 0] * numpy.sqrt(z0[0].real * z0[1].real)
+    s[1, 1] = (z0[0] + h[0, 0]) * (1 - h[1,1] * z0c[1]) + h[0,1] * h[1, 0] * z0c[1]
+
+    s /= (z0[0] + h[0, 0]) * (1 + h[1, 1] * z0[1]) - h[0,1] * h[1, 0] * z0[1]
+    return s
+
+
+def s_to_abcd(s, z0):
+    """
+    Scattering to Chain (ABCD)
+
+    Args:
+        s: The scattering matrix.
+        z0: The port impedances (Ohms).
+
+    Returns:
+        The chain matrix.
+    """
+    z0c = numpy.conj(z0)
+
+    ans = numpy.empty([2, 2], dtype=complex)
+    ans[0, 0] = (z0c[0] + s[0, 0] * z0[0]) * (1 - s[1, 1]) + s[0,1] * s[1, 0] * z0[0]
+    ans[0, 1] = (z0c[0] + s[0, 0] * z0[0]) * (z0c[1] + s[1,1] * z0[1]) - s[0,1] * s[1, 0] * z0[0] * z0[1]
+    ans[1, 0] = (1 - s[0, 0]) * (1 - s[1, 1]) - s[0,1] * s[1, 0]
+    ans[1, 1] = (1 - s[0, 0]) * (z0c[1] + s[1,1] * z0[1]) + s[0,1] * s[1, 0] * z0[1]
+
+    ans /= 2 * s[1, 0] * numpy.sqrt(z0[0].real * z0[1].real)
+    return ans
+
+
+def abcd_to_s(abcd, z0):
+    """
+    Chain (ABCD) to Scattering (S)
+
+    Args:
+        abcd: The chain matrix.
+        z0: The port impedances (Ohms).
+
+    Return:
+        The scattering matrix.
+    """
+    A = abcd[0, 0]
+    B = abcd[0, 1]
+    C = abcd[1, 0]
+    D = abcd[1, 1]
+
+    z0c = numpy.conj(z0)
+
+    s = numpy.empty([2, 2], dtype=complex)
+    s[0, 0] = A * z0[1] + B - C * z0c[0] * z0[1] - D * z0c[0]
+    s[0, 1] = 2 * (A * D - B * C) * numpy.sqrt(z0[0].real * z0[1].real)
+    s[1, 0] = 2 * numpy.sqrt(z0[0].real * z0[1].real)
+    s[1, 1] = -A * z0c[1] + B - C * z0[0] * z0c[1] + D * z0[0]
+
+    s /= A * z0[1] + B + C * z0[0] * z0[1] + D * z0[0]
+    return s
+
+
+def t_to_z(t, z0):
+    """
+    Chain Transfer (T) to Impedance (Z)
+
+    Args:
+        t: The chain transfer matrix.
+        z0: The port impedances (Ohms).
+
+    Returns:
+        The impedance matrix.
+    """
+    z0c = numpy.conj(z0)
+
+    z = numpy.empty([2, 2], dtype=complex)
+    z[0, 0] = z0c[0] * (t[0, 0] + t[0, 1]) + z0[0] * (t[1, 0] + t[1,1])
+    z[0, 1] = 2 * numpy.sqrt(z0[0].real * z0[1].real) * (t[0, 0] * t[1,1] - t[0,1] * t[1, 0])
+    z[1, 0] = 2 * numpy.sqrt(z0[0].real * z0[1].real)
+    z[1, 1] = z0c[1] * (t[0, 0] - t[1, 0]) - z0[1] * (t[0,1] - t[1,1])
+
+    z /= t[0, 0] + t[0, 1] - t[1, 0] - t[1,1]
+    return z
+
+
+def z_to_t(z, z0):
+    """
+    Impedance (Z) to Chain Transfer (T)
+
+    Args:
+        z: The impedance matrix.
+        z0: The port impedances (Ohms).
+
+    Returns:
+        The chain transfer matrix.
+    """
+    z0c = numpy.conj(z0)
+
+    t = numpy.empty([2, 2], dtype=complex)
+    t[0, 0] = (z[0, 0] + z0[0]) * (z[1, 1] + z0[1]) - z[0,1] * z[1, 0]
+    t[0, 1] = (z[0, 0] + z0[0]) * (z0c[1] - z[1,1]) + z[0,1] * z[1, 0]
+    t[1, 0] = (z[0, 0] - z0c[0]) * (z[1, 1] + z0[1]) - z[0,1] * z[1, 0]
+    t[1, 1] = (z0c[0] - z[0, 0]) * (z[1,1] - z0c[1]) + z[0,1] * z[1, 0]
+
+    t /= 2 * z[1, 0] * numpy.sqrt(z0[0].real * z0[1].real)
+    return t
+
+
+def t_to_y(t, z0):
+    """
+    Chain Transfer (T) to Admittance (Y)
+
+    Args:
+        t: The chain transfer matrix.
+        z0: The port impedances (Ohms).
+
+    Returns:
+        The admittance matrix.
+    """
+    z0c = numpy.conj(z0)
+
+    y = numpy.empty([2, 2], dtype=complex)
+    y[0, 0] = z0c[1] * (t[0, 0] - t[1, 0]) - z0[1] * (t[0, 1] - t[1,1])
+    y[0, 1] = -2 * numpy.sqrt(z0[0].real * z0[1].real) * (t[0, 0] * t[1,1] - t[0,1] * t[1, 0])
+    y[1, 0] = -2 * numpy.sqrt(z0[0].real * z0[1].real)
+    y[1, 1] = z0c[0] * (t[0, 0] + t[0,1]) + z0[0] * (t[1, 0] + t[1,1])
+
+    y /= t[0, 0] * z0c[0] * z0c[1] - t[0, 1] * z0c[0] * z0[1] + t[1, 0] * z0[0] * z0c[1] - t[1,1] * z0[0] * z0[1]
+    return y
+
+
+def y_to_t(y, z0):
+    """
+    Admittance (Y) to Chain Transfer (T)
+
+    Args:
+        y: The admittance matrix.
+        z0: The port impedances (Ohms).
+
+    Returns:
+        The chain transfer matrix.
+    """
+    z0c = numpy.conj(z0)
+
+    t = numpy.empty([2, 2], dtype=complex)
+    t[0, 0] = (-1 - y[0, 0] * z0[0]) * (1 + y[1, 1] * z0[1]) + y[0,1] * y[1, 0] * z0[0] * z0[1]
+    t[0, 1] = (1 + y[0, 0] * z0[0]) * (1 - y[1,1] * z0c[1]) + y[0,1] * y[1, 0] * z0[0] * z0c[1]
+    t[1, 0] = (y[0, 0] * z0c[0] - 1) * (1 + y[1, 1] * z0[1]) - y[0,1] * y[1, 0] * z0c[0] * z0[1]
+    t[1, 1] = (1 - y[0, 0] * z0c[0]) * (1 - y[1,1] * z0c[1]) - y[0,1] * y[1, 0] * z0c[0] * z0c[1]
+
+    t /= 2 * y[1, 0] * numpy.sqrt(z0[0].real * z0[1].real)
+    return t
+
+
+def t_to_h(t, z0):
+    """
+    Chain Transfer (T) to Hybrid (H)
+
+    Args:
+        t: The chain transfer matrix.
+        z0: The port impedances (Ohms).
+
+    Returns:
+        The hybrid matrix.
+    """
+    z0c = numpy.conj(z0)
+
+
+    h = numpy.empty([2, 2], dtype=complex)
+    h[0, 0] = z0c[1]*(t[0, 0] * z0c[0] + t[1, 0] * z0[0]) - z0[1] * (t[0, 1] * z0c[0] + t[1,1] * z0[0])
+    h[0, 1] = 2 * numpy.sqrt(z0[0].real * z0[1].real) * (t[0, 0] * t[1,1] - t[0,1] * t[1, 0])
+    h[1, 0] = -2 * numpy.sqrt(z0[0].real * z0[1].real)
+    h[1, 1] = t[0, 0] + t[0,1] - t[1, 0] - t[1,1]
+
+    h /= z0c[1] * (t[0, 0] - t[1, 0]) - z0[1] * (t[0, 1] - t[1,1])
+    return h
+
+
+def h_to_t(h, z0):
+    """
+    Hybrid (H) to Chain Transfer (T)
+
+    Args:
+        t: The hybrid matrix.
+        z0: The port impedances (Ohms).
+
+    Returns:
+        The chain transfer matrix.
+    """
+    z0c = numpy.conj(z0)
+
+    t = numpy.empty([2, 2], dtype=complex)
+    t[0, 0] = (-h[0, 0] - z0[0]) * (1 + h[1, 1] * z0[1]) + h[0,1] * h[1, 0] * z0[1]
+    t[0, 1] = (h[0, 0] + z0[0]) * (1 - h[1,1] * z0c[1]) + h[0,1] * h[1, 0] * z0c[1]
+    t[1, 0] = (z0c[0] - h[0, 0]) * (1 + h[1, 1] * z0[1]) + h[0,1] * h[1, 0] * z0[1]
+    t[1, 1] = (h[0, 0] - z0c[0]) * (1 - h[1,1] * z0c[1]) + h[0,1] * h[1, 0] * z0c[1]
+
+    t /= 2 * h[1, 0] * numpy.sqrt(z0[0].real * z0[1].real)
+    return t
+
+
+def t_to_abcd(t, z0):
+    """
+    Chain Transfer (T) to Chain (ABCD)
+
+    Args:
+        t: The chain transfer matrix.
+        z0: The port impedances (Ohms).
+
+    Returns:
+        The chain matrix.
+    """
+    z0c = numpy.conj(z0)
+    ans = numpy.empty([2, 2], dtype=complex)
+    ans[0, 0] = z0c[0] * (t[0, 0] + t[0, 1]) + z0[0] * (t[1, 0] + t[1, 1])
+    ans[0, 1] = z0c[1] * (t[0, 0] * z0c[0] + t[1, 0] * z0[0]) - z0[1] * (t[0, 1] * z0c[0] + t[1, 1] * z0[0])
+    ans[1, 0] = t[0, 0] + t[0, 1] - t[1, 0] - t[1, 1]
+    ans[1, 1] = z0c[1] * (t[0, 0] - t[1, 0]) - z0[1] * (t[0, 1] - t[1, 1])
+
+    ans /= 2 * numpy.sqrt(z0[0].real * z0[1].real)
+    return ans
+
+
+def abcd_to_t(abcd, z0):
+    """
+    Chain (ABCD) to Chain Transfer (T)
+
+    Args:
+        abcd: The chain matrix.
+        z0: The port impedances (Ohms).
+
+    Returns:
+        The chain transfer matrix.
+    """
+    # Break out the components
+    A = abcd[0, 0]
+    B = abcd[0, 1]
+    C = abcd[1, 0]
+    D = abcd[1, 1]
+
+    z0c = numpy.conj(z0)
+
+    t = numpy.empty([2, 2], dtype=complex)
+    t[0, 0] = A * z0[1] + B + C * z0[0] * z0[1] + D * z0[0]
+    t[0, 1] = A * z0c[1] - B + C * z0[0] * z0c[1] - D * z0[0]
+    t[1, 0] = A * z0[1] + B - C * z0c[0] * z0[1] - D * z0c[0]
+    t[1, 1] = A * z0c[1] - B - C * z0c[0] * z0c[1] + D * z0c[0]
+
+    t /= 2 * numpy.sqrt(z0[0].real * z0[1].real)
+    return t

spar_tests.py

@@ -0,0 +1,58 @@
+# IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES. VOL 42, NO 2. FEBRUARY 1994
+# Conversions Between S, Z, Y, h, ABCD, and T Parameters which are Valid for Complex Source and Load Impedances
+# Dean A. Frickey, Member, EEE
+
+# Tables I and II
+
+import numpy as np
+from two_port_conversions import *
+
+""" Testing """
+# Analog Devices - HMC455LP3 - S parameters
+# High output IP3 GaAs InGaP Heterojunction Bipolar Transistor
+
+# MHz S (Magntidue and Angle (deg))
+# 1487.273 0.409 160.117 4.367 163.864 0.063 115.967 0.254 -132.654
+
+s11 = 0.409 * np.exp(1j * np.radians(160.117))
+s12 = 4.367 * np.exp(1j * np.radians(163.864))
+s21 = 0.063 * np.exp(1j * np.radians(115.967))
+s22 = 0.254 * np.exp(1j * np.radians(-132.654))
+
+s_orig = np.array([[s11, s12], [s21, s22]])
+
+# Data specified at 50 Ohms (adding small complex component to test conversions)
+z1, z2 = 50 + 0.01j, 50 - 0.02j
+z0 = np.array([z1, z2])
+
+""" Conversions """
+print(f'Original S: \n{s_orig}\n')
+
+# S --> Z --> T --> Z --> S
+z = s_to_z(s_orig, z0)
+t = z_to_t(z, z0)
+z = t_to_z(t, z0)
+s = z_to_s(z, z0)
+print(f'Test (S --> Z --> T --> Z --> S): \n{s}\n')
+
+# S --> Y --> T --> Y --> S
+y = s_to_y(s_orig, z0)
+t = y_to_t(y, z0)
+y = t_to_y(t, z0)
+s = y_to_s(y, z0)
+print(f'Test (S --> Y --> T --> Y --> S): \n{s}\n')
+
+# S --> H --> T --> H --> S
+h = s_to_h(s_orig, z0)
+t = h_to_t(h, z0)
+h = t_to_h(t, z0)
+s = h_to_s(h, z0)
+print(f'Test (S --> H --> T --> H --> S): \n{s}\n')
+
+# S --> ABCD --> T --> ABCD --> S
+abcd = s_to_abcd(s_orig, z0)
+t = abcd_to_t(abcd, z0)
+abcd = t_to_abcd(t, z0)
+s = abcd_to_s(abcd, z0)
+print(f'Test (S --> ABCD --> T --> ABCD --> S): \n{s}\n')
+

spconv.py

@@ -0,0 +1,268 @@
+import scipy
+import numpy
+from numpy.typing import ArrayLike, NDArray
+from numpy.linalg import pinv
+from numpy import sqrt, real, abs, pi
+
+
+def diag(twod):
+    # numpy.diag() but makes a stack of diagonal matrices
+    return numpy.einsum('ij,jk->ijk', twod, numpy.eye(twod.shape[-1], dtype=twod.dtype))
+
+
+def s2z(s, zref):
+    # G0_inv @ inv(I - S) @ (S Z0 + Z0*) @ G0
+    #  Where Z0 is diag(zref) and G0 = diag(1/sqrt(abs(real(zref))))
+    nf = s.shape[-1]
+    I = numpy.eye(nf)[None, :, :]
+    zref = numpy.array(zref, copy=False)
+    gref = 1 / sqrt(abs(zref.real))
+    z = diag(1 / gref) @ pinv(I - s) @ ( S @ diag(zref) + diag(zref).conj()) @ diag(gref)
+    return z
+
+
+def change_of_zref(
+        s,  # (nf, np, np)
+        zref_old,  # (nf, np)
+        zref_new,  # (nf, np)
+        ):
+    # Change-of-Z0 to Z0'
+    # S' = inv(A) @ (S - rho*) @ inv(I - rho @ S) @ A*
+    #  A = inv(G0') @ G0 @ inv(I - rho*)   (diagonal)
+    #  rho = (Z0' - Z0) @ inv(Z0' + Z0)    (diagonal)
+
+    I = numpy.zeros_like(SL)
+    numpy.einsum('...jj->...j', I)[...] = 1
+
+    zref_old = numpy.array(zref_old, copy=False)
+    zref_new = numpy.array(zref_new, copy=False)
+
+    g_old = 1 / sqrt(abs(zref_old.real))
+    r_new = sqrt(abs(zref_new.real))
+
+    rhov = (zref_new - zref_old) / (zref_new + zref_old)
+    av = r_new * g_old / (1 - rhov.conj())
+
+    s_new = diag(1 / av) @ (s - diag(rhov.conj())) @ pinv(I[None, :] - diag(rhov) @ S) @ diag(av.conj())
+    return s_new
+
+
+def embedding(
+        See,  # (nf, ne, ne)
+        Sei,  # (nf, ne, ni)
+        Sie,  # (nf, ni, ne)
+        Sii,  # (nf, ni, ni)
+        SL,   # (nf, ni, ni)
+        ):
+    # Reveyrand, doi:10.1109/INMMIC.2018.8430023
+
+    I = numpy.zeros_like(SL)
+    numpy.einsum('...jj->...j', I)[...] = 1
+
+    S_tot = See + Sei @ pinv(I - SL @ Sii) @ SL @ Sie
+    return S_tot
+
+
+def deembedding(
+        Sei,  # (nf, ne, ni)
+        Sie,  # (nf, ni, ne)
+        Sext, # (nf, ne, ne)
+        See,  # (nf, ne, ne)
+        Si,   # (nf, ni, ni)
+        ):
+    # Reveyrand, doi:10.1109/INMMIC.2018.8430023
+    # Requires balanced number of ports, similar to VNA calibration
+    Sei_inv = pinv(Sei)
+    Sdif = Sext - See
+    return Sei_inv @ Sdif @ pinv(Sie + Sii @ Sei_inv @ Sdif)
+
+
+def thru_with_Zref_change(
+        zref0,  # (nf,)
+        zref1,  # (nf,)
+        ):
+    s = numpy.empty(tuple(zref0.shape) + (2, 2), dtype=complex)
+    s[..., 0, 0] = zref1 - zref0
+    s[..., 0, 1] = 2 * sqrt(zref0 * zref1)
+    s[..., 1, 0] = s[..., 0, 1]
+    s[..., 1, 1] = -s[..., 0, 0]
+
+    s /= zref0 + zref1
+    return s
+
+
+def propagation_matrix(mode_neffs: ArrayLike, wavelength: float, distance: float):
+    eigenv = numpy.array(mode_neffs, copy=False) * 2 * pi / wavelength
+    prop_diag = numpy.diag(numpy.exp(distance * 1j * numpy.hstack((eigenv, eigenv))))
+    prop_matrix = numpy.roll(prop_diag, len(eigenv), axis=0)
+    return prop_matrix
+
+
+def connect_s(
+        A: NDArray[numpy.complex128],
+        k: int,
+        B: NDArray[numpy.complex128],
+        l: int,
+        ) -> NDArray[numpy.complex128]:
+    """
+    TODO
+    freq x ...  x n x n
+
+    Based on skrf implementation
+
+    Connect two n-port networks' s-matrices together.
+
+    Specifically, connect port `k` on network `A` to port `l` on network
+    `B`. The resultant network has nports = (A.rank + B.rank-2); first
+    (A.rank - 1) ports are from `A`, remainder are from B.
+
+    Assumes same reference impedance for both `k` and `l`; may need to
+    connect an "impedance mismatch" thru element first!
+
+    Args:
+        A: S-parameter matrix of `A`, shape is fxnxn
+        k: port index on `A` (port indices start from 0)
+        B: S-parameter matrix of `B`, shape is fxnxn
+        l: port index on `B`
+
+    Returns:
+        new S-parameter matrix
+    """
+    if k > A.shape[-1] - 1 or l > B.shape[-1] - 1:
+        raise ValueError("port indices are out of range")
+
+    #C = scipy.sparse.block_diag((A, B), dtype=complex)
+    #return innerconnect_s(C, k, A.shape[0] + l)
+
+    nA = A.shape[-1]
+    nB = B.shape[-1]
+    nC = nA + nB - 2
+    assert numpy.array_equal(A.shape[:-2], B.shape[:-2])
+
+    ll = slice(l, l + 1)
+    kk = slice(k, k + 1)
+
+    denom = 1 - A[..., kk, kk] * B[..., ll, ll]
+    Anew = A + A[..., kk, :] * B[..., ll, ll] * A[..., :, kk] / denom
+    Bnew = A[..., kk, :] * B[..., :, ll] / denom
+    Anew = numpy.delete(Anew, (k,), 1)
+    Anew = numpy.delete(Anew, (k,), 2)
+    Bnew = numpy.delete(Bnew, (l,), 1)
+    Bnew = numpy.delete(Bnew, (l,), 2)
+
+    dtype = (A[0, 0] * B[0, 0]).dtype
+    C = numpy.zeros(tuple(A.shape[:-2]) + (nC, nC), dtype=dtype)
+    C[..., :nA - 1, :nA - 1] = Anew
+    C[..., nA - 1:, nA - 1:] = Bnew
+    return C
+
+
+def innerconnect_s(
+        S: NDArray[numpy.complex128],
+        k: int,
+        l: int,
+        ) -> NDArray[numpy.complex128]:
+    """
+    TODO
+    freq x ...  x n x n
+
+    Based on skrf implementation
+
+
+    Connect two ports of a single n-port network's s-matrix.
+    Specifically, connect port `k`  to port `l` on `S`. This results in
+    a (n-2)-port network.
+
+    Assumes same reference impedance for both `k` and `l`; may need to
+    connect an "impedance mismatch" thru element first!
+
+    Args:
+        S: S-parameter matrix of `S`, shape is fxnxn
+        k: port index on `S` (port indices start from 0)
+        l: port index on `S`
+
+    Returns:
+        new S-parameter matrix
+
+    Notes:
+        - Compton, R.C., "Perspectives in microwave circuit analysis",
+            doi:10.1109/MWSCAS.1989.101955
+        - Filipsson, G., "A New General Computer Algorithm for S-Matrix Calculation
+            of Interconnected Multiports",
+            doi:10.1109/EUMA.1981.332972
+    """
+    if k > S.shape[-1] - 1 or l > S.shape[-1] - 1:
+        raise ValueError("port indices are out of range")
+
+    ll = slice(l, l + 1)
+    kk = slice(k, k + 1)
+
+    mkl = 1 - S[..., kk, ll]
+    mlk = 1 - S[..., ll, kk]
+    C = S + (
+              S[..., kk, :] * S[..., :, l] * mlk
+            + S[..., ll, :] * S[..., :, k] * mkl
+            + S[..., kk, :] * S[..., l, l] * S[..., :, kk]
+            + S[..., ll, :] * S[..., k, k] * S[..., :, ll]
+        ) / (
+            mlk * mkl - S[..., kk, kk] * S[..., ll, ll]
+        )
+
+    # remove connected ports
+    C = numpy.delete(C, (k, l), 1)
+    C = numpy.delete(C, (k, l), 2)
+
+    return C
+
+
+def s2abcd(
+        S: NDArray[numpy.complex128],
+        z0: NDArray[numpy.complex128],
+        ) -> NDArray[numpy.complex128]:
+    assert numpy.array_equal(S.shape[:2] == (2, 2))
+    Z1, Z2 = z0
+    cross = S[0, 1] * S[1, 0]
+
+    T = numpy.empty_like(S, dtype=complex)
+    T[0, 0, :] = (Z1.conj() + S[0, 0] * Z1) * (1 - S[1, 1]) + cross * Z1    # A numerator
+    T[0, 1, :] = (Z1.conj() + S[0, 0] * Z1) * (Z1.conj() + S[1, 1] * Z2) - cross * Z1 * Z2  # B numerator
+    T[1, 0, :] = (1 - S[0, 0]) * (1 - S[1, 1]) - cross                      # C numerator
+    T[1, 1, :] = (1 - S[0, 0]) * (Z2.conj() + S[1, 1] * Z2) + cross * Z2    # D numerator
+    det = 2 * S[1, 0] * numpy.sqrt(Z1.real * Z2.real)
+    T /= det
+    return T
+
+
+def generalize_S(
+        S: NDArray[numpy.complex128],
+        r0: float,
+        z0: NDArray[numpy.complex128],
+        ) -> NDArray[numpy.complex128]:
+    g = (z0 - r0) / (z0 + r0)
+    D = numpy.diag((1 - g) / numpy.abs(1 - g.conj()) * numpy.sqrt(1 - numpy.abs(g * g.conj())))
+    G = numpy.diag(g)
+    U = numpy.eye(S.shape[-1]).reshape((S.ndim - 2) * (1,) + (S.shape[-2], S.shape[-1]))
+    S_gen = pinv(D.conj()) @ (S - G.conj()) @ pinv(U - G @ S) @ D
+    return S_gen
+
+
+def change_R0(
+        S: NDArray[numpy.complex128],
+        r1: float,
+        r2: float,
+        ) -> NDArray[numpy.complex128]:
+    g = (r2 - r1) / (r2 + r1)
+    U = numpy.eye(S.shape[-1]).reshape((S.ndim - 2) * (1,) + (S.shape[-2], S.shape[-1]))
+    G = U * g
+    S_r2 = (S - G) @ pinv(U - G @ S)
+    return S_r2
+
+# Zc = numpy.sqrt(B / C)
+# gamma = numpy.arccosh(A) / L_TL
+# n_eff = -1j * gamma * c_light / (2 * pi * f)
+# n_eff_grp = n_eff + f * diff(n_eff) / diff(f)
+# attenuation = (1 - S[0, 0] * S[0, 0].conj()) / (S[1, 0] * S[1, 0].conj())
+# R = numpy.real(gamma * Zc)
+# C = numpy.real(gamma / Zc)
+# L = numpy.imag(gamma * Zc) / (-1j * 2 * pi * f)
+# G = numpy.imag(gamma / Zc) / (-1j * 2 * pi * f)