add wave ports

added path integral plotting functionality

added convenience function for voltage integral

docs/api/plugins/smatrix.rst

@@ -11,6 +11,7 @@ Scattering Matrix Calculator
    tidy3d.plugins.smatrix.Port
    tidy3d.plugins.smatrix.ModalPortDataArray
    tidy3d.plugins.smatrix.TerminalComponentModeler
+   tidy3d.plugins.smatrix.TerminalPortDataArray
    tidy3d.plugins.smatrix.LumpedPort
-   tidy3d.plugins.smatrix.LumpedPortDataArray
-   tidy3d.plugins.smatrix.CoaxialLumpedPort
+   tidy3d.plugins.smatrix.CoaxialLumpedPort
+   tidy3d.plugins.smatrix.WavePort

tests/test_plugins/terminal_component_modeler_def.py

@@ -1,9 +1,13 @@
+from typing import Union
+
 import numpy as np
 import tidy3d as td
+import tidy3d.plugins.microwave as microwave
 from tidy3d.plugins.smatrix import (
     CoaxialLumpedPort,
     LumpedPort,
     TerminalComponentModeler,
+    WavePort,
 )
 
 # Microstrip dimensions
@@ -25,8 +29,8 @@
 freq_stop = 10e9
 
 # Coaxial dimensions
-Rinner = 0.2768 * 1e-3
-Router = 1.0 * 1e-3
+Rinner = 0.2768 * 1e3
+Router = 1.0 * 1e3
 
 
 def make_simulation(planar_pec: bool, length: float = None, auto_grid: bool = True):
@@ -176,7 +180,13 @@ def make_coaxial_simulation(length: float = None, auto_grid: bool = True):
 
     # Spatial grid specification
     if auto_grid:
-        grid_spec = td.GridSpec.auto(min_steps_per_wvl=10, wavelength=td.C_0 / freq_stop)
+        mesh_override = td.MeshOverrideStructure(
+            geometry=td.Box(center=(0, 0, 0), size=(2 * Router, 2 * Router, coax_length)),
+            dl=(50, 50, 10e3),
+        )
+        grid_spec = td.GridSpec.auto(
+            min_steps_per_wvl=10, wavelength=td.C_0 / freq_stop, override_structures=[mesh_override]
+        )
     else:
         grid_spec = td.GridSpec.uniform(wavelength0 / 11)
 
@@ -246,6 +256,10 @@ def make_coaxial_component_modeler(
     length: float = None,
     port_refinement: bool = True,
     auto_grid: bool = True,
+    port_types: tuple[Union[CoaxialLumpedPort, WavePort], Union[CoaxialLumpedPort, WavePort]] = (
+        CoaxialLumpedPort,
+        CoaxialLumpedPort,
+    ),
     **kwargs,
 ):
     if length:
@@ -255,24 +269,55 @@ def make_coaxial_component_modeler(
 
     sim = make_coaxial_simulation(length=coax_length, auto_grid=auto_grid)
 
-    center_src1 = [0, 0, -coax_length / 2]
-
-    port_cells = None
-    if port_refinement:
-        port_cells = 21
+    def make_port(center, direction, type, name) -> Union[CoaxialLumpedPort, WavePort]:
+        if type is CoaxialLumpedPort:
+            port_cells = None
+            if port_refinement:
+                port_cells = 21
+            port = CoaxialLumpedPort(
+                center=center,
+                outer_diameter=2 * Router,
+                inner_diameter=2 * Rinner,
+                normal_axis=2,
+                direction=direction,
+                name=name,
+                num_grid_cells=port_cells,
+                impedance=reference_impedance,
+            )
+        else:
+            mean_radius = (Router + Rinner) / 2
+            voltage_center = list(center)
+            voltage_center[0] += mean_radius
+            voltage_size = [Router - Rinner, 0, 0]
+
+            port = WavePort(
+                center=center,
+                size=[2 * Router, 2 * Router, 0],
+                direction=direction,
+                name=name,
+                mode_spec=td.ModeSpec(num_modes=1),
+                mode_index=0,
+                voltage_integral=microwave.VoltageIntegralAxisAligned(
+                    center=voltage_center,
+                    size=voltage_size,
+                    extrapolate_to_endpoints=True,
+                    snap_path_to_grid=True,
+                    sign="+",
+                ),
+                current_integral=microwave.CustomCurrentIntegral2D.from_circular_path(
+                    center=center,
+                    radius=mean_radius,
+                    num_points=41,
+                    normal_axis=2,
+                    clockwise=False,
+                ),
+            )
+        return port
 
-    port_1 = CoaxialLumpedPort(
-        center=center_src1,
-        outer_diameter=2 * Router,
-        inner_diameter=2 * Rinner,
-        normal_axis=2,
-        direction="+",
-        name="coax_port_1",
-        num_grid_cells=port_cells,
-        impedance=reference_impedance,
-    )
+    center_src1 = [0, 0, -coax_length / 2]
+    port_1 = make_port(center_src1, direction="+", type=port_types[0], name="coax_port_1")
     center_src2 = [0, 0, coax_length / 2]
-    port_2 = port_1.updated_copy(name="coax_port_2", center=center_src2, direction="-")
+    port_2 = make_port(center_src2, direction="-", type=port_types[1], name="coax_port_2")
     ports = [port_1, port_2]
     freqs = np.linspace(freq_start, freq_stop, 100)
 

tests/test_plugins/test_microwave.py

@@ -1,5 +1,6 @@
 """Test the microwave plugin."""
 
+import matplotlib.pyplot as plt
 import numpy as np
 import pydantic.v1 as pydantic
 import pytest
@@ -33,6 +34,8 @@
 )
 STRIP_WIDTH = 1.5
 STRIP_HEIGHT = 0.5
+COAX_R1 = 0.04
+COAX_R2 = 0.5
 
 SIM_Z = td.Simulation(
     size=(2, 1, 1),
@@ -103,6 +106,51 @@ def make_stripline_scalar_field_data_array(grid_key: str):
     return td.ScalarFieldDataArray(values, coords=dict(x=XS, y=YS, z=ZS, f=FS))
 
 
+def make_coaxial_field_data_array(grid_key: str):
+    """Populate FIELD_MONITOR with a coaxial transmission line mode."""
+
+    # Get a normalized electric field that represents the electric field within a coaxial cable transmission line.
+    def compute_coax_radial_electric(rin, rout, x, y, is_x):
+        # Radial distance
+        r = np.sqrt((x) ** 2 + (y) ** 2)
+        # Remove division by 0
+        r_valid = np.where(r == 0.0, 1, r)
+        # Compute current density so that the total current
+        # is 1 flowing through surfaces of constant r
+        # Extra r is for changing to Cartesian coordinates
+        denominator = 2 * np.pi * r_valid**2
+        if is_x:
+            Exy = np.where(r <= rin, 0, (x / denominator))
+            Exy = np.where(r >= rout, 0, Exy)
+        else:
+            Exy = np.where(r <= rin, 0, (y / denominator))
+            Exy = np.where(r >= rout, 0, Exy)
+        return Exy
+
+    XS, YS, ZS = get_xyz(FIELD_MONITOR, grid_key)
+    XGRID, YGRID = np.meshgrid(XS, YS, indexing="ij")
+    XGRID = XGRID.reshape((len(XS), len(YS), 1, 1))
+    YGRID = YGRID.reshape((len(XS), len(YS), 1, 1))
+    values = np.zeros((len(XS), len(YS), len(ZS), len(FS)))
+
+    XGRID = np.broadcast_to(XGRID, values.shape)
+    YGRID = np.broadcast_to(YGRID, values.shape)
+
+    is_x = grid_key[1] == "x"
+    if grid_key[0] == "E":
+        field = compute_coax_radial_electric(COAX_R1, COAX_R2, XGRID, YGRID, is_x)
+    else:
+        field = compute_coax_radial_electric(COAX_R1, COAX_R2, XGRID, YGRID, not is_x)
+        # H field is perpendicular and oriented for positive z propagation
+        # We want to compute Hx which is -Ey/eta0, Hy which is Ex/eta0
+        if is_x:
+            field /= -ETA_0
+        else:
+            field /= ETA_0
+
+    return td.ScalarFieldDataArray(field, coords=dict(x=XS, y=YS, z=ZS, f=FS))
+
+
 def make_field_data():
     return FieldData(
         monitor=FIELD_MONITOR,
@@ -116,6 +164,19 @@ def make_field_data():
     )
 
 
+def make_coax_field_data():
+    return FieldData(
+        monitor=FIELD_MONITOR,
+        Ex=make_coaxial_field_data_array("Ex"),
+        Ey=make_coaxial_field_data_array("Ey"),
+        Hx=make_coaxial_field_data_array("Hx"),
+        Hy=make_coaxial_field_data_array("Hy"),
+        symmetry=SIM_Z.symmetry,
+        symmetry_center=SIM_Z.center,
+        grid_expanded=SIM_Z.discretize_monitor(FIELD_MONITOR),
+    )
+
+
 @pytest.mark.parametrize("axis", [0, 1, 2])
 def test_voltage_integral_axes(axis):
     """Check VoltageIntegralAxisAligned runs."""
@@ -210,6 +271,15 @@ def test_current_missing_fields():
     with pytest.raises(DataError):
         _ = current_integral.compute_current(SIM_Z_DATA["ExHx"])
 
+    current_integral = CurrentIntegralAxisAligned(
+        center=center,
+        size=[length, length, 0],
+        sign="+",
+    )
+
+    with pytest.raises(DataError):
+        _ = current_integral.compute_current(SIM_Z_DATA["ExHx"])
+
 
 def test_time_monitor_voltage_integral():
     """Check VoltageIntegralAxisAligned runs on time domain data."""
@@ -417,6 +487,20 @@ def test_fields_missing_custom_current_integral():
         current_integral.compute_current(SIM_Z_DATA["ExHx"])
 
 
+def test_wrong_monitor_data():
+    """Check that the function arguments to integrals are correctly validated."""
+    length = 0.5
+    size = [0, 0, 0]
+    size[1] = length
+    # Make box
+    vertices = [(0.2, -0.2), (0.2, 0.2), (-0.2, 0.2), (-0.2, -0.2), (0.2, -0.2)]
+    current_integral = CustomCurrentIntegral2D(axis=2, position=0, vertices=vertices)
+    with pytest.raises(DataError):
+        current_integral.compute_current(SIM_Z.sources[0].source_time)
+    with pytest.raises(DataError):
+        current_integral.compute_current(em_field=SIM_Z.sources[0].source_time)
+
+
 def test_time_monitor_custom_current_integral():
     length = 0.5
     size = [0, 0, 0]
@@ -435,3 +519,109 @@ def test_mode_solver_custom_current_integral():
     vertices = [(0.2, -0.2), (0.2, 0.2), (-0.2, 0.2), (-0.2, -0.2), (0.2, -0.2)]
     current_integral = CustomCurrentIntegral2D(axis=2, position=0, vertices=vertices)
     current_integral.compute_current(SIM_Z_DATA["mode"])
+
+
+def test_custom_current_integral_normal_y():
+    current_integral = CustomCurrentIntegral2D.from_circular_path(
+        center=(0, 0, 0), radius=0.4, num_points=31, normal_axis=1, clockwise=False
+    )
+    current_integral.compute_current(SIM_Z_DATA["field"])
+
+
+def test_custom_path_integral_accuracy():
+    """Test the accuracy of the custom path integral."""
+    field_data = make_coax_field_data()
+
+    current_integral = CustomCurrentIntegral2D.from_circular_path(
+        center=(0, 0, 0), radius=0.4, num_points=31, normal_axis=2, clockwise=False
+    )
+    current = current_integral.compute_current(field_data)
+    assert np.allclose(current.values, 1.0 / ETA_0, rtol=0.01)
+
+    current_integral = CustomCurrentIntegral2D.from_circular_path(
+        center=(0, 0, 0), radius=0.4, num_points=31, normal_axis=2, clockwise=True
+    )
+    current = current_integral.compute_current(field_data)
+    assert np.allclose(current.values, -1.0 / ETA_0, rtol=0.01)
+
+
+def test_impedance_accuracy_on_coaxial():
+    """Test the accuracy of the ImpedanceCalculator."""
+    field_data = make_coax_field_data()
+    # Setup path integrals
+    mean_radius = (COAX_R2 + COAX_R1) * 0.5
+    size = [COAX_R2 - COAX_R1, 0, 0]
+    center = [mean_radius, 0, 0]
+
+    voltage_integral = VoltageIntegralAxisAligned(
+        center=center, size=size, sign="-", extrapolate_to_endpoints=True, snap_path_to_grid=True
+    )
+    current_integral = CustomCurrentIntegral2D.from_circular_path(
+        center=(0, 0, 0), radius=mean_radius, num_points=31, normal_axis=2, clockwise=False
+    )
+
+    def impedance_of_coaxial_cable(r1, r2, wave_impedance=td.ETA_0):
+        return np.log(r2 / r1) * wave_impedance / 2 / np.pi
+
+    Z_analytic = impedance_of_coaxial_cable(COAX_R1, COAX_R2)
+
+    # Compute impedance using the tool
+    Z_calculator = ImpedanceCalculator(
+        voltage_integral=voltage_integral, current_integral=current_integral
+    )
+    Z_calc = Z_calculator.compute_impedance(field_data)
+    assert np.allclose(Z_calc, Z_analytic, rtol=0.04)
+
+
+def test_path_integral_plotting():
+    """Test that all types of path integrals correctly plot themselves."""
+
+    mean_radius = (COAX_R2 + COAX_R1) * 0.5
+    size = [COAX_R2 - COAX_R1, 0, 0]
+    center = [mean_radius, 0, 0]
+
+    voltage_integral = VoltageIntegralAxisAligned(
+        center=center, size=size, sign="-", extrapolate_to_endpoints=True, snap_path_to_grid=True
+    )
+
+    current_integral = CustomCurrentIntegral2D.from_circular_path(
+        center=(0, 0, 0), radius=0.4, num_points=31, normal_axis=2, clockwise=False
+    )
+
+    ax = voltage_integral.plot(z=0)
+    current_integral.plot(z=0, ax=ax)
+    plt.close()
+
+    # Test off center plotting
+    ax = voltage_integral.plot(z=2)
+    current_integral.plot(z=2, ax=ax)
+    plt.close()
+
+    # Plot
+    voltage_integral = CustomVoltageIntegral2D(
+        axis=1, position=0, vertices=[(-1, -1), (0, 0), (1, 1)]
+    )
+
+    current_integral = CurrentIntegralAxisAligned(
+        center=(0, 0, 0),
+        size=(2, 0, 1),
+        sign="-",
+        extrapolate_to_endpoints=False,
+        snap_contour_to_grid=False,
+    )
+
+    ax = voltage_integral.plot(y=0)
+    current_integral.plot(y=0, ax=ax)
+    plt.close()
+
+    # Test off center plotting
+    ax = voltage_integral.plot(y=2)
+    current_integral.plot(y=2, ax=ax)
+    plt.close()
+
+
+def test_creation_from_terminal_positions():
+    """Test creating an VoltageIntegralAxisAligned using terminal positions."""
+    _ = VoltageIntegralAxisAligned.from_terminal_positions(
+        plus_terminal=2, minus_terminal=1, y=2.2, z=1
+    )