providing support for both pseudo and power wave formulations

Author: dmarek-flex

tests/test_plugins/smatrix/test_terminal_component_modeler.py

@@ -40,7 +40,7 @@ def run_component_modeler(monkeypatch, modeler: TerminalComponentModeler):
     monkeypatch.setattr(
         TerminalComponentModeler,
         "_compute_F",
-        lambda matrix: 1.0 / (2.0 * np.sqrt(np.abs(matrix) + 1e-4)),
+        lambda matrix, wave_type: 1.0 / (2.0 * np.sqrt(np.abs(matrix) + 1e-4)),
     )
     monkeypatch.setattr(
         TerminalComponentModeler,

tidy3d/plugins/smatrix/component_modelers/terminal.py

@@ -2,7 +2,7 @@
 
 from __future__ import annotations
 
-from typing import Optional, Union
+from typing import Literal, Optional, Union
 
 import numpy as np
 import pydantic.v1 as pd
@@ -33,7 +33,17 @@
 
 class TerminalComponentModeler(AbstractComponentModeler):
     """Tool for modeling two-terminal multiport devices and computing port parameters
-    with lumped and wave ports."""
+    with lumped and wave ports.
+
+    Notes
+    -----
+    **References**
+
+        [1] R. B. Marks and D. F. Williams, "A general waveguide circuit theory,"
+            J. Res. Natl. Inst. Stand. Technol., vol. 97, pp. 533, 1992.
+
+        [2] D. M. Pozar, Microwave Engineering, 4th ed. Hoboken, NJ, USA: John Wiley & Sons, 2012.
+    """
 
     ports: tuple[TerminalPortType, ...] = pd.Field(
         (),
@@ -220,7 +230,9 @@ def _internal_construct_smatrix(self, batch_data: BatchData) -> TerminalPortData
         # loop through source ports
         for port_in in self.ports:
             sim_data = batch_data[self._task_name(port=port_in)]
-            a, b = self.compute_power_wave_amplitudes_at_each_port(port_impedances, sim_data)
+            a, b = self.compute_wave_amplitudes_at_each_port(
+                port_impedances, sim_data, wave_type="pseudo"
+            )
             indexer = {"f": a.f, "port_in": port_in.name, "port_out": a.port}
             a_matrix.loc[indexer] = a
             b_matrix.loc[indexer] = b
@@ -280,10 +292,13 @@ def _check_grid_size_at_wave_ports(simulation: Simulation, ports: list[WavePort]
                         "for the simulation passed to the 'TerminalComponentModeler'."
                     )
 
-    def compute_power_wave_amplitudes_at_each_port(
-        self, port_reference_impedances: PortDataArray, sim_data: SimulationData
+    def compute_wave_amplitudes_at_each_port(
+        self,
+        port_reference_impedances: PortDataArray,
+        sim_data: SimulationData,
+        wave_type: Literal["pseudo", "power"] = "pseudo",
     ) -> tuple[PortDataArray, PortDataArray]:
-        """Compute the incident and reflected power wave amplitudes at each port.
+        """Compute the incident and reflected amplitudes at each port.
         The computed amplitudes have not been normalized.
 
         Parameters
@@ -292,11 +307,14 @@ def compute_power_wave_amplitudes_at_each_port(
             Reference impedance at each port.
         sim_data : :class:`.SimulationData`
             Results from the simulation.
+        wave_type : Literal['pseudo', 'power']
+            The type of waves computed, either pseudo waves defined by Equation 53 and Equation 54 in [1],
+            or power waves defined by Equation 4.67 in [2].
 
         Returns
         -------
         tuple[:class:`.PortDataArray`, :class:`.PortDataArray`]
-            Incident (a) and reflected (b) power wave amplitudes at each port.
+            Incident (a) and reflected (b) wave amplitudes at each port.
         """
         port_names = [port.name for port in self.ports]
         values = np.zeros(
@@ -331,14 +349,41 @@ def compute_power_wave_amplitudes_at_each_port(
         V_numpy = np.where(negative_real_Z, -V_numpy, V_numpy)
         Z_numpy = np.where(negative_real_Z, -Z_numpy, Z_numpy)
 
-        F_numpy = TerminalComponentModeler._compute_F(Z_numpy)
+        F_numpy = TerminalComponentModeler._compute_F(Z_numpy, wave_type)
+
+        b_Zref = Z_numpy
+        if wave_type == "power":
+            b_Zref = np.conj(Z_numpy)
 
+        # Equations 53 and 54 from [1]
         # Equation 4.67 - Pozar - Microwave Engineering 4ed
         a.values = F_numpy * (V_numpy + Z_numpy * I_numpy)
-        b.values = F_numpy * (V_numpy - np.conj(Z_numpy) * I_numpy)
+        b.values = F_numpy * (V_numpy - b_Zref * I_numpy)
 
         return a, b
 
+    def compute_power_wave_amplitudes_at_each_port(
+        self, port_reference_impedances: PortDataArray, sim_data: SimulationData
+    ) -> tuple[PortDataArray, PortDataArray]:
+        """DEPRECATED: Compute the incident and reflected power wave amplitudes at each port.
+        The computed amplitudes have not been normalized.
+
+        Parameters
+        ----------
+        port_reference_impedances : :class:`.PortDataArray`
+            Reference impedance at each port.
+        sim_data : :class:`.SimulationData`
+            Results from the simulation.
+
+        Returns
+        -------
+        tuple[:class:`.PortDataArray`, :class:`.PortDataArray`]
+            Incident (a) and reflected (b) power wave amplitudes at each port.
+        """
+        return self.compute_wave_amplitudes_at_each_port(
+            port_reference_impedances, sim_data, wave_type="power"
+        )
+
     @staticmethod
     def compute_port_VI(
         port_out: TerminalPortType, sim_data: SimulationData
@@ -365,7 +410,7 @@ def compute_port_VI(
     def compute_power_wave_amplitudes(
         port: Union[LumpedPort, CoaxialLumpedPort], sim_data: SimulationData
     ) -> tuple[FreqDataArray, FreqDataArray]:
-        """Compute the incident and reflected power wave amplitudes at a lumped port.
+        """DEPRECATED: Compute the incident and reflected power wave amplitudes at a lumped port.
         The computed amplitudes have not been normalized.
 
         Parameters
@@ -388,9 +433,10 @@ def compute_power_wave_amplitudes(
 
     @staticmethod
     def compute_power_delivered_by_port(
-        port: Union[LumpedPort, CoaxialLumpedPort], sim_data: SimulationData
+        port: Union[LumpedPort, CoaxialLumpedPort],
+        sim_data: SimulationData,
     ) -> FreqDataArray:
-        """Compute the power delivered to the network by a lumped port.
+        """DEPRECATED: Compute the power delivered to the network by a lumped port.
 
         Parameters
         ----------
@@ -412,7 +458,7 @@ def compute_power_delivered_by_port(
     def ab_to_s(
         a_matrix: TerminalPortDataArray, b_matrix: TerminalPortDataArray
     ) -> TerminalPortDataArray:
-        """Get the scattering matrix given the power wave matrices."""
+        """Get the scattering matrix given the wave amplitude matrices."""
         # Ensure dimensions are ordered properly
         a_matrix = a_matrix.transpose(*TerminalPortDataArray._dims)
         b_matrix = b_matrix.transpose(*TerminalPortDataArray._dims)
@@ -428,44 +474,63 @@ def ab_to_s(
 
     @staticmethod
     def s_to_z(
-        s_matrix: TerminalPortDataArray, reference: Union[complex, PortDataArray]
+        s_matrix: TerminalPortDataArray,
+        reference: Union[complex, PortDataArray],
+        wave_type: Literal["pseudo", "power"] = "pseudo",
     ) -> DataArray:
-        """Get the impedance matrix given the scattering matrix and a reference impedance."""
+        """Get the impedance matrix given the scattering matrix and a reference impedance.
+
+        Parameters
+        ----------
+        s_matrix : :class:`.TerminalPortDataArray`
+            Scattering matrix computed using either the pseudo or power wave formulation.
+        reference : Union[complex, :class:`.PortDataArray`]
+            The reference impedance used at each port.
+        wave_type : Literal['pseudo', 'power']
+            The type of wave amplitudes used for computing the scattering matrix, either pseudo waves
+            defined by Equation 53 and Equation 54 in [1] or power waves defined by Equation 4.67 in [2].
+        """
 
         # Ensure dimensions are ordered properly
         z_matrix = s_matrix.transpose(*TerminalPortDataArray._dims).copy(deep=True)
         s_vals = z_matrix.values
-        eye = np.eye(len(s_matrix.port_out.values), len(s_matrix.port_in.values))
+        eye = np.eye(len(s_matrix.port_out.values), len(s_matrix.port_in.values))[np.newaxis, :, :]
+        # Ensure that Zport, F, and Finv act as diagonal matrices when multiplying by left or right
+        shape_left = (len(s_matrix.f), len(s_matrix.port_out), 1)
+        shape_right = (len(s_matrix.f), 1, len(s_matrix.port_in))
+        # Setup the port reference impedance array (scalar)
         if isinstance(reference, PortDataArray):
-            # From Equation 4.68 - Pozar - Microwave Engineering 4ed
-            # Ensure that Zport, F, and Finv act as diagonal matrices when multiplying by left or right
-            shape_left = (len(s_matrix.f), len(s_matrix.port_out), 1)
-            shape_right = (len(s_matrix.f), 1, len(s_matrix.port_in))
             Zport = reference.values.reshape(shape_right)
-            F = TerminalComponentModeler._compute_F(Zport).reshape(shape_right)
+            F = TerminalComponentModeler._compute_F(Zport, wave_type).reshape(shape_right)
             Finv = (1.0 / F).reshape(shape_left)
-            FinvSF = Finv * s_vals * F
-            RHS = eye * np.conj(Zport) + FinvSF * Zport
-            LHS = eye - FinvSF
-            z_vals = np.matmul(AbstractComponentModeler.inv(LHS), RHS)
         else:
-            # Simpler case when all port impedances are the same
-            z_vals = (
-                np.matmul(AbstractComponentModeler.inv(eye - s_vals), (eye + s_vals)) * reference
-            )
+            Zport = reference
+            F = TerminalComponentModeler._compute_F(Zport, wave_type)
+            Finv = 1.0 / F
+        # Use conjugfate when S matrix is power-wave based
+        if wave_type == "power":
+            Zport_mod = np.conj(Zport)
+        else:
+            Zport_mod = Zport
+
+        # From equation 74 from [1] for pseudo waves
+        # From Equation 4.68 - Pozar - Microwave Engineering 4ed for power waves
+        FinvSF = Finv * s_vals * F
+        RHS = eye * Zport_mod + FinvSF * Zport
+        LHS = eye - FinvSF
+        z_vals = np.linalg.solve(LHS, RHS)
 
         z_matrix.data = z_vals
         return z_matrix
 
     @cached_property
     def port_reference_impedances(self) -> PortDataArray:
-        """The reference impedance used at each port for definining power wave amplitudes.
+        """The reference impedance used at each port for definining wave amplitudes.
 
         Note
         ----
-        By default, we choose reference impedances to be the conjugate of the load impedance.
-        For wave ports, this corresponds with the conjugate of the characteristic impedance.
-        This choice corresponds with Equation 4.64 - Pozar - Microwave Engineering 4ed.
+        By default, we choose reference impedances to be the load impedance.
+        For wave ports, this corresponds with the the characteristic impedance.
         """
         return self._port_reference_impedances(self.batch_data)
 
@@ -493,16 +558,19 @@ def _port_reference_impedances(self, batch_data: BatchData) -> PortDataArray:
                 # LumpedPorts have a constant reference impedance
                 port_impedances.loc[{"port": port.name}] = np.full(len(self.freqs), port.impedance)
 
-        port_impedances = TerminalComponentModeler._set_port_data_array_attributes(
-            port_impedances.conj()
-        )
+        port_impedances = TerminalComponentModeler._set_port_data_array_attributes(port_impedances)
         return port_impedances
 
     @staticmethod
-    def _compute_F(Z_numpy: np.array):
+    def _compute_F(Z_numpy: np.array, wave_type: Literal["pseudo", "power"] = "pseudo"):
         """Helper to convert port impedance matrix to F, which is used for
-        computing generalized scattering parameters."""
-        return 1.0 / (2.0 * np.sqrt(np.real(Z_numpy)))
+        computing scattering parameters
+        """
+        # Defined in [2] after equation 4.67
+        if wave_type == "power":
+            return 1.0 / (2.0 * np.sqrt(np.real(Z_numpy)))
+        # Equation 75 from [1]
+        return np.sqrt(np.real(Z_numpy)) / (2.0 * np.abs(Z_numpy))
 
     @cached_property
     def _lumped_ports(self) -> list[AbstractLumpedPort]:
@@ -585,7 +653,7 @@ def get_antenna_metrics_data(
     ) -> AntennaMetricsData:
         """Calculate antenna parameters using superposition of fields from multiple port excitations.
 
-        The method computes the radiated far fields and port excitation power wave amplitudes
+        The method computes the radiated far fields and port excitation wave amplitudes
         for a superposition of port excitations, which can be used to analyze antenna radiation
         characteristics.
 
@@ -621,7 +689,7 @@ def get_antenna_metrics_data(
         else:
             rad_mon = self.get_radiation_monitor_by_name(monitor_name)
 
-        # Create data arrays for holding the superposition of all port power wave amplitudes
+        # Create data arrays for holding the superposition of all port wave amplitudes
         f = list(rad_mon.freqs)
         coords = {"f": f, "port": port_names}
         a_sum = PortDataArray(np.zeros((len(f), len(port_names)), dtype=complex), coords=coords)
@@ -632,8 +700,8 @@ def get_antenna_metrics_data(
             sim_data_port = self.batch_data[self._task_name(port=port)]
             radiation_data = sim_data_port[rad_mon.name]
 
-            a, b = self.compute_power_wave_amplitudes_at_each_port(
-                self.port_reference_impedances, sim_data_port
+            a, b = self.compute_wave_amplitudes_at_each_port(
+                self.port_reference_impedances, sim_data_port, wave_type="pseudo"
             )
             # Select a possible subset of frequencies
             a = a.sel(f=f)
@@ -652,7 +720,7 @@ def get_antenna_metrics_data(
             a = scale_factor * a
             b = scale_factor * b
 
-            # Combine the possibly scaled directivity data and the power wave amplitudes
+            # Combine the possibly scaled directivity data and the wave amplitudes
             if combined_directivity_data is None:
                 combined_directivity_data = scaled_directivity_data
             else: