diff --git a/nom-eme.py b/nom-eme.py new file mode 100644 index 0000000..d69bdff --- /dev/null +++ b/nom-eme.py @@ -0,0 +1,350 @@ + +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 + 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. + + 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: + C: 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) + + +def innerconnect_s( + S: NDArray[numpy.complex128], + k: int, + l: int, + ) -> NDArray[numpy.complex128]: + """ + TODO + n x n x freq + 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. + + 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") + + l = [l] + k = [k] + + mkl = 1 - S[k, l] + mlk = 1 - S[l, k] + C = S + ( + S[k, :] * S[:, l] * mlk + + S[l, :] * S[:, k] * mkl + + S[k, :] * S[l, l] * S[:, k] + + S[l, :] * S[k, k] * S[:, l] + ) / ( + mlk * mkl - S[k, k] * S[l, l] + ) + + # remove connected ports + C = npy.delete(C, (k, l), 1) + C = npy.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[0]) + 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[0]) + 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) diff --git a/spar.py b/spar.py new file mode 100644 index 0000000..2efd573 --- /dev/null +++ b/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 diff --git a/spar_tests.py b/spar_tests.py new file mode 100644 index 0000000..73ec60f --- /dev/null +++ b/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') + diff --git a/spconv.py b/spconv.py new file mode 100644 index 0000000..f76e2db --- /dev/null +++ b/spconv.py @@ -0,0 +1,86 @@ +import numpy +from numpy import sqrt, real, abs +from numpy.linalg import pinv + + +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.eye(SL.shape[-1])[None, :, :] + 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.eye(SL.shape[-1])[None, :, :] + 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,) + ): + nf = zref0.shape[0] + s = numpy.empty((nf, 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 + +