fix impedance calculation for lossy transmission lines

Author: dmarek-flex

CHANGELOG.md

@@ -28,6 +28,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Solver error for EME simulations with bends, introduced when support for 2D EME simulations was added.
 - Internal interpolation errors with some versions of `xarray` and `numpy`.
 - If `ModeSpec.angle_rotation=True` for a mode object, validate that the structure rotation can be successfully done. Also, error if the medium cannot be rotated (e.g. anisotropic or custom medium), which would previously have just produced wrong results.
+- Characteristic impedance calculations in the `ImpedanceCalculator` using definitions that rely on flux, which were giving incorrect results for lossy transmission lines.
 
 ### Changed
 - Relaxed bounds checking of path integrals during `WavePort` validation.

tidy3d/components/data/monitor_data.py

@@ -660,7 +660,7 @@ def intensity(self) -> ScalarFieldDataArray:
         return intensity.squeeze(dim=normal_dim, drop=True)
 
     @property
-    def poynting(self) -> ScalarFieldDataArray:
+    def complex_poynting(self) -> ScalarFieldDataArray:
         """Time-averaged Poynting vector for frequency-domain data associated to a 2D monitor,
         projected to the direction normal to the monitor plane."""
 
@@ -677,27 +677,36 @@ def poynting(self) -> ScalarFieldDataArray:
         e2_h1 = e2 * h1.conj()
 
         e_x_h_star = e1_h2 - e2_h1
-        poynting = 0.5 * np.real(e_x_h_star)
+        return 0.5 * e_x_h_star
 
-        return poynting
+    @property
+    def poynting(self) -> ScalarFieldDataArray:
+        """Time-averaged Poynting vector for frequency-domain data associated to a 2D monitor,
+        projected to the direction normal to the monitor plane."""
+        return self.complex_poynting.real
 
     def package_flux_results(self, flux_values: DataArray) -> Any:
         """How to package flux"""
         return FluxDataArray(flux_values)
 
     @cached_property
-    def flux(self) -> FluxDataArray:
+    def complex_flux(self) -> FluxDataArray:
         """Flux for data corresponding to a 2D monitor."""
 
         # Compute flux by integrating Poynting vector in-plane
         d_area = self._diff_area
-        poynting = self.poynting
+        poynting = self.complex_poynting
 
         flux_values = poynting * d_area
         flux_values = flux_values.sum(dim=d_area.dims)
 
         return self.package_flux_results(flux_values)
 
+    @cached_property
+    def flux(self) -> FluxDataArray:
+        """Flux for data corresponding to a 2D monitor."""
+        return self.complex_flux.real
+
     @cached_property
     def mode_area(self) -> FreqModeDataArray:
         r"""Effective mode area corresponding to a 2D monitor.

tidy3d/plugins/microwave/impedance_calculator.py

@@ -8,7 +8,9 @@
 import pydantic.v1 as pd
 
 from tidy3d.components.base import Tidy3dBaseModel
+from tidy3d.components.data.data_array import FreqDataArray, FreqModeDataArray, TimeDataArray
 from tidy3d.components.data.monitor_data import FieldTimeData
+from tidy3d.components.monitor import ModeMonitor, ModeSolverMonitor
 from tidy3d.constants import OHM
 from tidy3d.exceptions import ValidationError
 from tidy3d.log import log
@@ -60,9 +62,9 @@ def compute_impedance(self, em_field: MonitorDataTypes) -> IntegralResultTypes:
         AxisAlignedPathIntegral._check_monitor_data_supported(em_field=em_field)
 
         # If both voltage and current integrals have been defined then impedance is computed directly
-        if self.voltage_integral:
+        if self.voltage_integral is not None:
             voltage = self.voltage_integral.compute_voltage(em_field)
-        if self.current_integral:
+        if self.current_integral is not None:
             current = self.current_integral.compute_current(em_field)
 
         # If only one of the integrals has been provided, then the computation falls back to using
@@ -71,20 +73,31 @@ def compute_impedance(self, em_field: MonitorDataTypes) -> IntegralResultTypes:
         # a time signal, then it is real and flux corresponds to the instantaneous power. Otherwise
         # the input field is in frequency domain, where flux indicates the time-averaged power
         # 0.5*Re(V*conj(I)).
-        if not self.voltage_integral:
-            flux = em_field.flux
+        # We explicitly take the real part, in case Bloch BCs were used in the simulation.
+        flux_sign = 1.0
+        # Determine flux sign
+        if isinstance(em_field.monitor, ModeSolverMonitor):
+            flux_sign = 1 if em_field.monitor.direction == "+" else -1
+        if isinstance(em_field.monitor, ModeMonitor):
+            flux_sign = 1 if em_field.monitor.store_fields_direction == "+" else -1
+
+        if self.voltage_integral is None:
+            flux = flux_sign * em_field.complex_flux
             if isinstance(em_field, FieldTimeData):
-                voltage = flux.abs / current
+                impedance = flux / np.real(current) ** 2
             else:
-                voltage = 2 * flux.abs / np.conj(current)
-        if not self.current_integral:
-            flux = em_field.flux
+                impedance = 2 * flux / (current * np.conj(current))
+        elif self.current_integral is None:
+            flux = flux_sign * em_field.complex_flux
             if isinstance(em_field, FieldTimeData):
-                current = flux.abs / voltage
+                impedance = np.real(voltage) ** 2 / flux
             else:
-                current = np.conj(2 * flux.abs / voltage)
-
-        impedance = voltage / current
+                impedance = (voltage * np.conj(voltage)) / (2 * np.conj(flux))
+        else:
+            if isinstance(em_field, FieldTimeData):
+                impedance = np.real(voltage) / np.real(current)
+            else:
+                impedance = voltage / current
         impedance = ImpedanceCalculator._set_data_array_attributes(impedance)
         return impedance
 
@@ -101,6 +114,13 @@ def check_voltage_or_current(cls, val, values):
     @staticmethod
     def _set_data_array_attributes(data_array: IntegralResultTypes) -> IntegralResultTypes:
         """Helper to set additional metadata for ``IntegralResultTypes``."""
+        # Determine type based on coords present
+        if "mode_index" in data_array.coords:
+            data_array = FreqModeDataArray(data_array)
+        elif "f" in data_array.coords:
+            data_array = FreqDataArray(data_array)
+        else:
+            data_array = TimeDataArray(data_array)
         data_array.name = "Z0"
         return data_array.assign_attrs(units=OHM, long_name="characteristic impedance")