diff --git a/tests/test_components/test_microwave.py b/tests/test_components/test_microwave.py
index 5fb08815e..62ca54718 100644
--- a/tests/test_components/test_microwave.py
+++ b/tests/test_components/test_microwave.py
@@ -3,13 +3,18 @@
 from __future__ import annotations
 
 from math import isclose
+from typing import Literal
 
+import matplotlib.pyplot as plt
 import numpy as np
+import pydantic.v1 as pd
 import pytest
 import xarray as xr
+from shapely.geometry import LineString, Polygon
+from shapely.plotting import plot_line, plot_polygon
 
+import tidy3d as td
 from tidy3d.components.data.monitor_data import FreqDataArray
-from tidy3d.components.microwave.data.monitor_data import AntennaMetricsData
 from tidy3d.components.microwave.formulas.circuit_parameters import (
     capacitance_colinear_cylindrical_wire_segments,
     capacitance_rectangular_sheets,
@@ -17,10 +22,207 @@
     mutual_inductance_colinear_wire_segments,
     total_inductance_colinear_rectangular_wire_segments,
 )
+from tidy3d.components.microwave.path_integral_factory import (
+    make_current_integral,
+    make_voltage_integral,
+)
+from tidy3d.components.microwave.path_spec_generator import PathSpecGenerator
+from tidy3d.components.types import Ax, Shapely
 from tidy3d.constants import EPSILON_0
+from tidy3d.exceptions import ValidationError
 
 from ..test_data.test_monitor_data import make_directivity_data
 
+mm = 1e3
+
+MAKE_PLOTS = False
+if MAKE_PLOTS:
+    # Interative plotting for debugging
+    from matplotlib import use
+
+    use("TkAgg")
+
+
+def make_mw_sim(
+    use_2D: bool = False,
+    colocate: bool = False,
+    transmission_line_type: Literal["microstrip", "cpw", "coax", "stripline"] = "microstrip",
+) -> td.Simulation:
+    """Helper to create a microwave simulation with a single type of transmission line present."""
+
+    freq_start = 1e9
+    freq_stop = 10e9
+
+    freq0 = (freq_start + freq_stop) / 2
+    fwidth = freq_stop - freq_start
+    freqs = np.arange(freq_start, freq_stop, 1e9)
+
+    run_time = 60 / fwidth
+
+    length = 40 * mm
+    width = 3 * mm
+    height = 1 * mm
+    thickness = 0.2 * mm
+    sim_width = length
+
+    pec = td.PEC
+    if use_2D:
+        thickness = 0.0
+        pec = td.PEC2D
+
+    epsr = 4.4
+    diel = td.Medium(permittivity=epsr)
+
+    metal_geos = []
+
+    if transmission_line_type == "microstrip":
+        metal_geos.append(
+            td.Box(
+                center=[0, 0, height + thickness / 2],
+                size=[td.inf, width, thickness],
+            )
+        )
+    elif transmission_line_type == "cpw":
+        metal_geos.append(
+            td.Box(
+                center=[0, 0, height + thickness / 2],
+                size=[td.inf, width, thickness],
+            )
+        )
+        gnd_width = 10 * width
+        gap = width / 5
+        gnd_shift = gnd_width / 2 + gap + width / 2
+        metal_geos.append(
+            td.Box(
+                center=[0, -gnd_shift, height + thickness / 2],
+                size=[td.inf, gnd_width, thickness],
+            )
+        )
+        metal_geos.append(
+            td.Box(
+                center=[0, gnd_shift, height + thickness / 2],
+                size=[td.inf, gnd_width, thickness],
+            )
+        )
+    elif transmission_line_type == "coax":
+        metal_geos.append(
+            td.GeometryGroup(
+                geometries=(
+                    td.ClipOperation(
+                        operation="difference",
+                        geometry_a=td.Cylinder(
+                            axis=0, radius=2 * mm, center=(0, 0, 5 * mm), length=td.inf
+                        ),
+                        geometry_b=td.Cylinder(
+                            axis=0, radius=1.8 * mm, center=(0, 0, 5 * mm), length=td.inf
+                        ),
+                    ),
+                    td.Cylinder(axis=0, radius=0.6 * mm, center=(0, 0, 5 * mm), length=td.inf),
+                )
+            )
+        )
+    elif transmission_line_type == "stripline":
+        metal_geos.append(
+            td.Box(
+                center=[0, 0, 0],
+                size=[td.inf, width, thickness],
+            )
+        )
+        gnd_width = 10 * width
+        metal_geos.append(
+            td.Box(
+                center=[0, 0, height + thickness],
+                size=[td.inf, gnd_width, thickness],
+            )
+        )
+        metal_geos.append(
+            td.Box(
+                center=[0, 0, -height - thickness],
+                size=[td.inf, gnd_width, thickness],
+            )
+        )
+    else:
+        raise AssertionError("Incorrect argument")
+
+    metal_structures = [td.Structure(geometry=geo, medium=pec) for geo in metal_geos]
+    substrate = td.Structure(
+        geometry=td.Box(
+            center=[0, 0, 0],
+            size=[td.inf, td.inf, 2 * height],
+        ),
+        medium=diel,
+    )
+
+    structures = [substrate, *metal_structures]
+    boundary_spec = td.BoundarySpec(
+        x=td.Boundary(plus=td.PML(), minus=td.PML()),
+        y=td.Boundary(plus=td.PML(), minus=td.PML()),
+        z=td.Boundary(plus=td.PML(), minus=td.PECBoundary()),
+    )
+
+    size_sim = [
+        length + 2 * width,
+        sim_width,
+        20 * mm + height + thickness,
+    ]
+    center_sim = [0, 0, size_sim[2] / 2]
+    # Slightly different setup for stripline substrate sandwiched between ground planes
+    if transmission_line_type == "stripline":
+        center_sim[2] = 0
+        boundary_spec = td.BoundarySpec(
+            x=td.Boundary(plus=td.PML(), minus=td.PML()),
+            y=td.Boundary(plus=td.PML(), minus=td.PML()),
+            z=td.Boundary(plus=td.PML(), minus=td.PML()),
+        )
+    size_port = [0, sim_width, size_sim[2]]
+    center_port = [0, 0, center_sim[2]]
+    mode_spec = td.ModeSpec(
+        num_modes=4, target_neff=1.8, microwave_mode_spec=td.MicrowaveModeSpec()
+    )
+
+    mode_monitor = td.ModeMonitor(
+        center=center_port, size=size_port, freqs=freqs, name="mode_1", colocate=colocate
+    )
+
+    gaussian = td.GaussianPulse(freq0=freq0, fwidth=fwidth)
+    mode_src = td.ModeSource(
+        center=(-length / 2, 0, center_sim[2]),
+        size=size_port,
+        direction="+",
+        mode_spec=mode_spec,
+        mode_index=0,
+        source_time=gaussian,
+    )
+    sim = td.Simulation(
+        center=center_sim,
+        size=size_sim,
+        grid_spec=td.GridSpec.uniform(dl=0.1 * mm),
+        structures=structures,
+        sources=[mode_src],
+        monitors=[mode_monitor],
+        run_time=run_time,
+        boundary_spec=boundary_spec,
+        plot_length_units="mm",
+        symmetry=(0, 0, 0),
+    )
+    return sim
+
+
+def plot_auto_path_spec(
+    path_spec: td.CompositeCurrentIntegralSpec, geoms: list[Shapely], ax: Ax = None
+) -> Ax:
+    """Helper to plot composite path specifications along with the Shapely geometries used to generate them."""
+    if ax is None:
+        _, ax = plt.subplots(1, 1, tight_layout=True, figsize=(15, 15))
+
+    for geom in geoms:
+        if isinstance(geom, Polygon):
+            plot_polygon(geom, ax=ax)
+        elif isinstance(geom, LineString):
+            plot_line(geom, ax=ax)
+    i_integral = make_current_integral(path_spec)
+    i_integral.plot(x=i_integral.center[0], ax=ax)
+
 
 def test_inductance_formulas():
     """Run the formulas for inductance and compare to precomputed results."""
@@ -68,7 +270,7 @@ def test_antenna_parameters():
     f = directivity_data.coords["f"]
     power_inc = FreqDataArray(0.8 * np.ones(len(f)), coords={"f": f})
     power_refl = 0.25 * power_inc
-    antenna_params = AntennaMetricsData.from_directivity_data(
+    antenna_params = td.AntennaMetricsData.from_directivity_data(
         directivity_data, power_inc, power_refl
     )
 
@@ -112,3 +314,208 @@ def test_antenna_parameters():
         antenna_params.partial_gain("invalid")
     with pytest.raises(ValueError):
         antenna_params.partial_realized_gain("invalid")
+
+
+@pytest.mark.parametrize("colocate", [False, True])
+@pytest.mark.parametrize("tline_type", ["microstrip", "cpw", "coax"])
+def test_auto_path_spec_canonical_shapes(colocate, tline_type):
+    """Test canonical transmission line types to make sure the correct path integrals are generated."""
+    # Interative plotting for debugging
+    # from matplotlib import use
+
+    # use("TkAgg")
+
+    sim = make_mw_sim(False, colocate, tline_type)
+    mode_monitor = sim.monitors[0]
+    modal_plane = td.Box(center=mode_monitor.center, size=mode_monitor.size)
+    comp_path_spec, geos = PathSpecGenerator.create_current_path_specs(
+        modal_plane,
+        sim.structures,
+        sim.grid,
+        sim.symmetry,
+        sim.bounding_box,
+        field_data_colocated=mode_monitor.colocate,
+    )
+
+    if tline_type == "coax":
+        assert len(comp_path_spec.path_specs) == 2
+        for path_spec in comp_path_spec.path_specs:
+            assert np.all(np.isclose(path_spec.center, (0, 0, 5 * mm)))
+    else:
+        assert len(comp_path_spec.path_specs) == 1
+        assert np.all(np.isclose(comp_path_spec.path_specs[0].center, (0, 0, 1.1 * mm)))
+
+    _, ax = plt.subplots(1, 1, tight_layout=True, figsize=(15, 15))
+    sim.plot(x=modal_plane.center[0], ax=ax, monitor_alpha=0)
+    sim.plot_grid(x=modal_plane.center[0], ax=ax, monitor_alpha=0)
+    plot_auto_path_spec(comp_path_spec, geos, ax)
+    ax.set_aspect("equal")
+    ax.set_xlim(modal_plane.bounds[0][1], modal_plane.bounds[1][1])
+    ax.set_ylim(modal_plane.bounds[0][2], modal_plane.bounds[1][2])
+    if MAKE_PLOTS:
+        plt.show()
+
+
+@pytest.mark.parametrize("use_2D", [False, True])
+@pytest.mark.parametrize("symmetry", [(0, 0, 1), (0, 1, 1), (0, 1, 0)])
+def test_auto_path_spec_advanced(use_2D, symmetry):
+    """The various symmetry permutations as well as with and without 2D structures."""
+    sim = make_mw_sim(use_2D, False, "stripline")
+
+    # Add shapes outside the portion considered for the symmetric simulation
+    bottom_left = td.Structure(
+        geometry=td.Box(
+            center=[0, -5 * mm, -5 * mm],
+            size=[td.inf, 1 * mm, 1 * mm],
+        ),
+        medium=td.PEC,
+    )
+    # Add shape only in the symmetric portion
+    top_right = td.Structure(
+        geometry=td.Box(
+            center=[0, 5 * mm, 5 * mm],
+            size=[td.inf, 1 * mm, 1 * mm],
+        ),
+        medium=td.PEC,
+    )
+
+    structures = [*list(sim.structures), bottom_left, top_right]
+    sim = sim.updated_copy(symmetry=symmetry, structures=structures)
+    mode_monitor = sim.monitors[0]
+
+    modal_plane = td.Box(center=mode_monitor.center, size=mode_monitor.size)
+    comp_path_spec, geos = PathSpecGenerator.create_current_path_specs(
+        modal_plane,
+        sim.structures,
+        sim.grid,
+        sim.symmetry,
+        sim.bounding_box,
+        field_data_colocated=mode_monitor.colocate,
+    )
+
+    if symmetry[1] == 1 and symmetry[2] == 1:
+        assert len(comp_path_spec.path_specs) == 7
+    else:
+        assert len(comp_path_spec.path_specs) == 5
+
+    _, ax = plt.subplots(1, 1, tight_layout=True, figsize=(15, 15))
+    sim.plot(x=modal_plane.center[0], ax=ax, monitor_alpha=0)
+    sim.plot_grid(x=modal_plane.center[0], ax=ax, monitor_alpha=0)
+    plot_auto_path_spec(comp_path_spec, geos, ax)
+    ax.set_aspect("equal")
+    ax.set_xlim(modal_plane.bounds[0][1], modal_plane.bounds[1][1])
+    ax.set_ylim(modal_plane.bounds[0][2], modal_plane.bounds[1][2])
+    if MAKE_PLOTS:
+        plt.show()
+
+
+def test_auto_path_spec_validation():
+    """Check that the auto path specs are validated properly."""
+
+    # First some quick sanity checks with the helper
+    test_path = td.Box(center=(0, 0, 0), size=(0, 0.9, 0.1))
+    test_shapely = [LineString([(-1, 0), (1, 0)])]
+    assert PathSpecGenerator._check_path_intersects_with_conductors(test_shapely, test_path)
+
+    test_path = td.Box(center=(0, 0, 0), size=(0, 2.1, 0.1))
+    test_shapely = [LineString([(-1, 0), (1, 0)])]
+    assert not PathSpecGenerator._check_path_intersects_with_conductors(test_shapely, test_path)
+
+    sim = make_mw_sim(False, False, "microstrip")
+    coax = td.GeometryGroup(
+        geometries=(
+            td.ClipOperation(
+                operation="difference",
+                geometry_a=td.Cylinder(axis=0, radius=2 * mm, center=(0, 0, 5 * mm), length=td.inf),
+                geometry_b=td.Cylinder(
+                    axis=0, radius=1.4 * mm, center=(0, 0, 5 * mm), length=td.inf
+                ),
+            ),
+            td.Cylinder(axis=0, radius=1 * mm, center=(0, 0, 5 * mm), length=td.inf),
+        )
+    )
+    coax_struct = td.Structure(geometry=coax, medium=td.PEC)
+    sim = sim.updated_copy(structures=[coax_struct])
+    mode_monitor = sim.monitors[0]
+    modal_plane = td.Box(center=mode_monitor.center, size=mode_monitor.size)
+    with pytest.raises(ValidationError):
+        PathSpecGenerator.create_current_path_specs(
+            modal_plane,
+            sim.structures,
+            sim.grid,
+            sim.symmetry,
+            sim.bounding_box,
+            field_data_colocated=mode_monitor.colocate,
+        )
+
+
+def test_path_integral_creation():
+    """Check that path integrals are correctly constructed from path specifications."""
+
+    path_spec = td.VoltageIntegralAxisAlignedSpec(center=(1, 2, 3), size=(0, 0, 1), sign="-")
+    voltage_integral = make_voltage_integral(path_spec)
+
+    path_spec = td.CurrentIntegralAxisAlignedSpec(center=(1, 2, 3), size=(0, 1, 1), sign="-")
+    current_integral = make_current_integral(path_spec)
+
+    path_spec = td.CustomVoltageIntegral2DSpec(vertices=[(0, 1), (0, 4)], axis=1, position=2)
+    voltage_integral = make_voltage_integral(path_spec)
+
+    path_spec = td.CustomCurrentIntegral2DSpec(
+        vertices=[
+            (0, 1),
+            (0, 4),
+            (3, 4),
+            (3, 1),
+        ],
+        axis=1,
+        position=2,
+    )
+    _ = make_current_integral(path_spec)
+
+    with pytest.raises(pd.ValidationError):
+        path_spec = td.CustomCurrentIntegral2DSpec(
+            vertices=[
+                (0, 1, 3),
+                (0, 4, 5),
+                (3, 4, 5),
+                (3, 1, 5),
+            ],
+            axis=1,
+            position=2,
+        )
+
+
+def test_microwave_mode_spec_validation():
+    """Check that the various allowed methods for supplying path specifications are validated."""
+
+    _ = td.MicrowaveModeSpec()
+    _ = td.MicrowaveModeSpec(voltage_spec=None, current_spec=None)
+
+    v_spec = td.VoltageIntegralAxisAlignedSpec(center=(1, 2, 3), size=(0, 0, 1), sign="-")
+    i_spec = td.CurrentIntegralAxisAlignedSpec(center=(1, 2, 3), size=(0, 1, 1), sign="-")
+
+    # All valid methods
+    _ = td.MicrowaveModeSpec(voltage_spec=(v_spec,), current_spec=(i_spec,))
+    _ = td.MicrowaveModeSpec(voltage_spec=(v_spec,))
+    _ = td.MicrowaveModeSpec(current_spec=(i_spec,))
+    _ = td.MicrowaveModeSpec(voltage_spec=(v_spec, v_spec), current_spec=(i_spec, i_spec))
+
+    # Different lengths is not valid
+    with pytest.raises(pd.ValidationError):
+        _ = td.MicrowaveModeSpec(voltage_spec=(v_spec, v_spec), current_spec=(i_spec,))
+    with pytest.raises(pd.ValidationError):
+        _ = td.MicrowaveModeSpec(voltage_spec=(v_spec,), current_spec=(i_spec, i_spec))
+
+    # Only one path spec missing is ok
+    _ = td.MicrowaveModeSpec(voltage_spec=(v_spec, None), current_spec=(None, i_spec))
+
+    # But at least one must be given for each pair
+    with pytest.raises(pd.ValidationError):
+        _ = td.MicrowaveModeSpec(voltage_spec=(v_spec, None), current_spec=(None, None))
+
+    with pytest.raises(pd.ValidationError):
+        _ = td.MicrowaveModeSpec(voltage_spec=(v_spec, None))
+
+    with pytest.raises(pd.ValidationError):
+        _ = td.MicrowaveModeSpec(current_spec=(None, i_spec))
diff --git a/tests/test_components/test_simulation.py b/tests/test_components/test_simulation.py
index 999ca4a00..6037f5524 100644
--- a/tests/test_components/test_simulation.py
+++ b/tests/test_components/test_simulation.py
@@ -3648,3 +3648,56 @@ def test_create_sim_multiphysics_with_incompatibilities():
                 ),
             ],
         )
+
+
+def test_validate_microwave_mode_spec_generation():
+    """Test that auto generation of path specs is correctly validated for currently unsupported structures."""
+    freq0 = 10e9
+    mm = 1e3
+    run_time_spec = td.RunTimeSpec(quality_factor=3.0)
+    size = (10 * mm, 10 * mm, 10 * mm)
+    size_mon = (0, 8 * mm, 8 * mm)
+
+    # Currently limited to generation of axis aligned boxes around conductors,
+    # so the path may intersect other nearby conductors, like in this coaxial cable
+    coaxial = td.Structure(
+        geometry=td.GeometryGroup(
+            geometries=(
+                td.ClipOperation(
+                    operation="difference",
+                    geometry_a=td.Cylinder(
+                        axis=0, radius=2.5 * mm, center=(0, 0, 0), length=td.inf
+                    ),
+                    geometry_b=td.Cylinder(
+                        axis=0, radius=1.3 * mm, center=(0, 0, 0), length=td.inf
+                    ),
+                ),
+                td.Cylinder(axis=0, radius=1 * mm, center=(0, 0, 0), length=td.inf),
+            )
+        ),
+        medium=td.PEC,
+    )
+    mode_spec = td.ModeSpec(
+        num_modes=2, target_neff=1.8, microwave_mode_spec=td.MicrowaveModeSpec()
+    )
+
+    mode_mon = td.ModeMonitor(
+        center=(0, 0, 0),
+        size=size_mon,
+        freqs=[freq0],
+        name="mode_1",
+        colocate=False,
+        mode_spec=mode_spec,
+    )
+    sim = td.Simulation(
+        run_time=run_time_spec,
+        size=size,
+        sources=[],
+        structures=[coaxial],
+        grid_spec=td.GridSpec.uniform(dl=0.1 * mm),
+        monitors=[mode_mon],
+    )
+
+    # check that validation error is caught
+    with pytest.raises(SetupError):
+        sim._validate_microwave_mode_specs()
diff --git a/tests/test_data/test_data_arrays.py b/tests/test_data/test_data_arrays.py
index f38f42572..8b00e367a 100644
--- a/tests/test_data/test_data_arrays.py
+++ b/tests/test_data/test_data_arrays.py
@@ -315,6 +315,16 @@ def test_abs():
     _ = data.abs
 
 
+def test_angle():
+    # Make sure works on real data and the type is correct
+    data = make_scalar_field_time_data_array("Ex")
+    angle_data = data.angle
+    assert type(data) is type(angle_data)
+    data = make_mode_amps_data_array()
+    angle_data = data.angle
+    assert type(data) is type(angle_data)
+
+
 def test_heat_data_array():
     T = [0, 1e-12, 2e-12]
     _ = td.HeatDataArray((1 + 1j) * np.random.random((3,)), coords={"T": T})
diff --git a/tests/test_plugins/test_microwave.py b/tests/test_plugins/test_microwave.py
index 0cea9b17f..7e72bfcc3 100644
--- a/tests/test_plugins/test_microwave.py
+++ b/tests/test_plugins/test_microwave.py
@@ -14,10 +14,18 @@
 import tidy3d as td
 import tidy3d.plugins.microwave as mw
 from tidy3d import FieldData
+from tidy3d.components.data.data_array import FreqModeDataArray
 from tidy3d.constants import ETA_0
 from tidy3d.exceptions import DataError
 
-from ..utils import get_spatial_coords_dict, run_emulated
+from ..utils import AssertLogLevel, get_spatial_coords_dict, run_emulated
+
+MAKE_PLOTS = False
+if MAKE_PLOTS:
+    # Interative plotting for debugging
+    from matplotlib import use
+
+    use("TkAgg")
 
 # Using similar code as "test_data/test_data_arrays.py"
 MON_SIZE = (2, 1, 0)
@@ -527,6 +535,41 @@ def test_custom_current_integral_normal_y():
     current_integral.compute_current(SIM_Z_DATA["field"])
 
 
+def test_composite_current_integral_warnings():
+    """Ensures that the checks function correctly on some test data."""
+    f = [2e9, 3e9, 4e9]
+    mode_index = list(np.arange(5))
+    coords = {"f": f, "mode_index": mode_index}
+    values = np.ones((3, 5))
+
+    path_spec = td.CurrentIntegralAxisAlignedSpec(center=(0, 0, 0), size=(2, 2, 0), sign="+")
+    composite_integral = mw.CompositeCurrentIntegral(
+        center=(0, 0, 0), size=(4, 4, 0), path_specs=[path_spec], sum_spec="split"
+    )
+
+    phase_diff = FreqModeDataArray(np.angle(values), coords=coords)
+    with AssertLogLevel(None):
+        assert composite_integral._check_phase_sign_consistency(phase_diff)
+
+    values[1, 2:] = -1
+    phase_diff = FreqModeDataArray(np.angle(values), coords=coords)
+    with AssertLogLevel("WARNING"):
+        assert not composite_integral._check_phase_sign_consistency(phase_diff)
+
+    values = np.ones((3, 5))
+    in_phase = FreqModeDataArray(values, coords=coords)
+    values = 0.5 * np.ones((3, 5))
+    out_phase = FreqModeDataArray(values, coords=coords)
+    with AssertLogLevel(None):
+        assert composite_integral._check_phase_amplitude_consistency(in_phase, out_phase)
+
+    values = 0.5 * np.ones((3, 5))
+    values[2, 4:] = 1.5
+    out_phase = FreqModeDataArray(values, coords=coords)
+    with AssertLogLevel("WARNING"):
+        assert not composite_integral._check_phase_amplitude_consistency(in_phase, out_phase)
+
+
 def test_custom_path_integral_accuracy():
     """Test the accuracy of the custom path integral."""
     field_data = make_coax_field_data()
@@ -785,13 +828,11 @@ def test_lobe_measurements(apply_cyclic_extension, include_endpoint):
 @pytest.mark.parametrize("min_value", [0.0, 1.0])
 def test_lobe_plots(min_value):
     """Run the lobe measurer on some test data and plot the results."""
-    # Interative plotting for debugging
-    # from matplotlib import use
-    # use("TkAgg")
     theta = np.linspace(0, 2 * np.pi, 301)
     Urad = np.cos(theta) ** 2 * np.cos(3 * theta) ** 2 + min_value
     lobe_measurer = mw.LobeMeasurer(angle=theta, radiation_pattern=Urad)
     _, ax = plt.subplots(1, 1, subplot_kw={"projection": "polar"})
     ax.plot(theta, Urad, "k")
     lobe_measurer.plot(0, ax)
-    plt.show()
+    if MAKE_PLOTS:
+        plt.show()
diff --git a/tidy3d/__init__.py b/tidy3d/__init__.py
index 6bc0bd31a..05fe4680d 100644
--- a/tidy3d/__init__.py
+++ b/tidy3d/__init__.py
@@ -17,6 +17,16 @@
 from tidy3d.components.microwave.data.monitor_data import (
     AntennaMetricsData,
 )
+from tidy3d.components.microwave.microwave_mode_spec import (
+    MicrowaveModeSpec,
+)
+from tidy3d.components.microwave.path_spec import (
+    CompositeCurrentIntegralSpec,
+    CurrentIntegralAxisAlignedSpec,
+    CustomCurrentIntegral2DSpec,
+    CustomVoltageIntegral2DSpec,
+    VoltageIntegralAxisAlignedSpec,
+)
 from tidy3d.components.spice.analysis.dc import (
     ChargeToleranceSpec,
     IsothermalSteadyChargeDCAnalysis,
@@ -439,6 +449,7 @@ def set_logging_level(level: str) -> None:
     "ChargeToleranceSpec",
     "ClipOperation",
     "CoaxialLumpedResistor",
+    "CompositeCurrentIntegralSpec",
     "ConstantDoping",
     "ConstantMobilityModel",
     "ContinuousWave",
@@ -449,8 +460,10 @@ def set_logging_level(level: str) -> None:
     "Coords1D",
     "CornerFinderSpec",
     "CurrentBC",
+    "CurrentIntegralAxisAlignedSpec",
     "CustomAnisotropicMedium",
     "CustomChargePerturbation",
+    "CustomCurrentIntegral2DSpec",
     "CustomCurrentSource",
     "CustomDebye",
     "CustomDrude",
@@ -463,6 +476,7 @@ def set_logging_level(level: str) -> None:
     "CustomPoleResidue",
     "CustomSellmeier",
     "CustomSourceTime",
+    "CustomVoltageIntegral2DSpec",
     "Cylinder",
     "DCCurrentSource",
     "DCVoltageSource",
@@ -580,6 +594,7 @@ def set_logging_level(level: str) -> None:
     "Medium2D",
     "MediumMediumInterface",
     "MeshOverrideStructure",
+    "MicrowaveModeSpec",
     "ModeAmpsDataArray",
     "ModeData",
     "ModeIndexDataArray",
@@ -682,6 +697,7 @@ def set_logging_level(level: str) -> None:
     "Updater",
     "VisualizationSpec",
     "VoltageBC",
+    "VoltageIntegralAxisAlignedSpec",
     "VoltageSourceType",
     "VolumetricAveraging",
     "YeeGrid",
diff --git a/tidy3d/components/data/data_array.py b/tidy3d/components/data/data_array.py
index caa053449..cb361578e 100644
--- a/tidy3d/components/data/data_array.py
+++ b/tidy3d/components/data/data_array.py
@@ -205,6 +205,13 @@ def abs(self):
         """Absolute value of data array."""
         return abs(self)
 
+    @property
+    def angle(self):
+        """Angle or phase value of data array."""
+        values = np.angle(self.values)
+        SelfType = type(self)
+        return SelfType(values, coords=self.coords)
+
     @property
     def is_uniform(self):
         """Whether each element is of equal value in the data array"""
diff --git a/tidy3d/components/data/monitor_data.py b/tidy3d/components/data/monitor_data.py
index b91f26ed2..0d45de72e 100644
--- a/tidy3d/components/data/monitor_data.py
+++ b/tidy3d/components/data/monitor_data.py
@@ -1602,11 +1602,19 @@ class ModeData(ModeSolverDataset, ElectromagneticFieldData):
 
     eps_spec: list[EpsSpecType] = pd.Field(
         None,
-        title="Permettivity Specification",
+        title="Permittivity Specification",
         description="Characterization of the permittivity profile on the plane where modes are "
         "computed. Possible values are 'diagonal', 'tensorial_real', 'tensorial_complex'.",
     )
 
+    Z0: Optional[FreqModeDataArray] = pd.Field(
+        None,
+        title="Characteristic Impedance",
+        description="Optional quantity calculated for transmission lines. "
+        "The characteristic impedance is only calculated when a :class:`MicrowaveModeSpec` "
+        "is provided to the :class:`ModeSpec` associated with this data.",
+    )
+
     @pd.validator("eps_spec", always=True)
     @skip_if_fields_missing(["monitor"])
     def eps_spec_match_mode_spec(cls, val, values):
@@ -2055,6 +2063,10 @@ def modes_info(self) -> xr.Dataset:
             info["wg TE fraction"] = self.pol_fraction_waveguide["te"]
             info["wg TM fraction"] = self.pol_fraction_waveguide["tm"]
 
+        if self.Z0 is not None:
+            info["Re(Z0)"] = self.Z0.real
+            info["Im(Z0)"] = self.Z0.imag
+
         return xr.Dataset(data_vars=info)
 
     def to_dataframe(self) -> DataFrame:
diff --git a/tidy3d/components/geometry/utils.py b/tidy3d/components/geometry/utils.py
index 8c166fdd2..4c82d04dd 100644
--- a/tidy3d/components/geometry/utils.py
+++ b/tidy3d/components/geometry/utils.py
@@ -2,12 +2,21 @@
 
 from __future__ import annotations
 
+from collections.abc import Iterable
 from enum import Enum
 from math import isclose
 from typing import Any, Optional, Union
 
 import numpy as np
-import pydantic
+import pydantic.v1 as pydantic
+from shapely.geometry import (
+    GeometryCollection,
+    MultiLineString,
+    MultiPoint,
+    MultiPolygon,
+    Polygon,
+)
+from shapely.geometry.base import BaseGeometry
 
 from tidy3d.components.base import Tidy3dBaseModel
 from tidy3d.components.geometry.base import Box
@@ -38,6 +47,46 @@
 ]
 
 
+def flatten_shapely_geometries(
+    geoms: Union[Shapely, Iterable[Shapely]], keep_types: tuple[type, ...] = (Polygon,)
+) -> list[Shapely]:
+    """
+    Flatten nested geometries into a flat list, while only keeping the specified types.
+
+    Recursively extracts and returns non-empty geometries of the given types from input geometries,
+    expanding any GeometryCollections or Multi* types.
+
+    Parameters
+    ----------
+    geoms : Union[Shapely, Iterable[Shapely]]
+        Input geometries to flatten.
+
+    keep_types : tuple[type, ...]
+        Geometry types to keep (e.g., (Polygon, LineString)). Default is
+        (Polygon).
+
+    Returns
+    -------
+    list[Shapely]
+        Flat list of non-empty geometries matching the specified types.
+    """
+    # Handle single Shapely object by wrapping it in a list
+    if isinstance(geoms, Shapely):
+        geoms = [geoms]
+
+    flat = []
+    for geom in geoms:
+        if geom.is_empty:
+            continue
+        if isinstance(geom, keep_types):
+            flat.append(geom)
+        elif isinstance(geom, (MultiPolygon, MultiLineString, MultiPoint, GeometryCollection)):
+            flat.extend(flatten_shapely_geometries(geom.geoms, keep_types))
+        elif isinstance(geom, BaseGeometry) and hasattr(geom, "geoms"):
+            flat.extend(flatten_shapely_geometries(geom.geoms, keep_types))
+    return flat
+
+
 def merging_geometries_on_plane(
     geometries: list[GeometryType],
     plane: Box,
@@ -373,6 +422,15 @@ class SnappingSpec(Tidy3dBaseModel):
         description="Describes how snapping positions will be chosen.",
     )
 
+    margin: Optional[
+        tuple[pydantic.NonNegativeInt, pydantic.NonNegativeInt, pydantic.NonNegativeInt]
+    ] = pydantic.Field(
+        (0, 0, 0),
+        title="Margin",
+        description="Number of additional grid points to consider when expanding or contracting "
+        "during snapping. Only applies when ``SnapBehavior`` is ``Expand`` or ``Contract``.",
+    )
+
 
 def get_closest_value(test: float, coords: np.ArrayLike, upper_bound_idx: int) -> float:
     """Helper to choose the closest value in an array to a given test value,
@@ -404,6 +462,8 @@ def get_lower_bound(
         using the index of the upper bound. If the test value is close to the upper
         bound, it assumes they are equal, and in that case the upper bound is returned.
         """
+        upper_bound_idx = min(upper_bound_idx, len(coords))
+        upper_bound_idx = max(upper_bound_idx, 0)
         if upper_bound_idx == len(coords):
             return coords[upper_bound_idx - 1]
         if upper_bound_idx == 0 or isclose(coords[upper_bound_idx], test, rel_tol=rel_tol):
@@ -417,6 +477,8 @@ def get_upper_bound(
         using the index of the upper bound. If the test value is close to the lower
         bound, it assumes they are equal, and in that case the lower bound is returned.
         """
+        upper_bound_idx = min(upper_bound_idx, len(coords))
+        upper_bound_idx = max(upper_bound_idx, 0)
         if upper_bound_idx == len(coords):
             return coords[upper_bound_idx - 1]
         if upper_bound_idx > 0 and isclose(coords[upper_bound_idx - 1], test, rel_tol=rel_tol):
@@ -424,7 +486,11 @@ def get_upper_bound(
         return coords[upper_bound_idx]
 
     def find_snapping_locations(
-        interval_min: float, interval_max: float, coords: np.ndarray, snap_type: SnapBehavior
+        interval_min: float,
+        interval_max: float,
+        coords: np.ndarray,
+        snap_type: SnapBehavior,
+        snap_margin: pydantic.NonNegativeInt,
     ) -> tuple[float, float]:
         """Helper that snaps a supplied interval [interval_min, interval_max] to a
         sorted array representing coordinate values.
@@ -436,9 +502,15 @@ def find_snapping_locations(
             min_snap = get_closest_value(interval_min, coords, min_upper_bound_idx)
             max_snap = get_closest_value(interval_max, coords, max_upper_bound_idx)
         elif snap_type == SnapBehavior.Expand:
+            min_upper_bound_idx -= snap_margin
+            max_upper_bound_idx += snap_margin
             min_snap = get_lower_bound(interval_min, coords, min_upper_bound_idx, rel_tol=rtol)
             max_snap = get_upper_bound(interval_max, coords, max_upper_bound_idx, rel_tol=rtol)
         else:  # SnapType.Contract
+            min_upper_bound_idx += snap_margin
+            max_upper_bound_idx -= snap_margin
+            if max_upper_bound_idx < min_upper_bound_idx:
+                raise SetupError("The supplied 'snap_buffer' is too large for this contraction.")
             min_snap = get_upper_bound(interval_min, coords, min_upper_bound_idx, rel_tol=rtol)
             max_snap = get_lower_bound(interval_max, coords, max_upper_bound_idx, rel_tol=rtol)
         return (min_snap, max_snap)
@@ -450,6 +522,7 @@ def find_snapping_locations(
     for axis in range(3):
         snap_location = snap_spec.location[axis]
         snap_type = snap_spec.behavior[axis]
+        snap_margin = snap_spec.margin[axis]
         if snap_type == SnapBehavior.Off:
             continue
         if snap_location == SnapLocation.Boundary:
@@ -460,7 +533,9 @@ def find_snapping_locations(
         box_min = min_b[axis]
         box_max = max_b[axis]
 
-        (new_min, new_max) = find_snapping_locations(box_min, box_max, snap_coords, snap_type)
+        (new_min, new_max) = find_snapping_locations(
+            box_min, box_max, snap_coords, snap_type, snap_margin
+        )
         min_b[axis] = new_min
         max_b[axis] = new_max
     return Box.from_bounds(min_b, max_b)
diff --git a/tidy3d/components/microwave/microwave_mode_spec.py b/tidy3d/components/microwave/microwave_mode_spec.py
new file mode 100644
index 000000000..f28774b0f
--- /dev/null
+++ b/tidy3d/components/microwave/microwave_mode_spec.py
@@ -0,0 +1,97 @@
+"""Specification for modes associated with transmission lines."""
+
+from __future__ import annotations
+
+from typing import Optional
+
+import pydantic.v1 as pd
+
+from tidy3d.components.base import Tidy3dBaseModel
+from tidy3d.exceptions import SetupError
+
+from .path_spec import CurrentPathSpecTypes, VoltagePathSpecTypes
+
+
+class MicrowaveModeSpec(Tidy3dBaseModel):
+    """Specification for computing transmission line voltages and currents in mode solvers.
+
+    The :class:`.MicrowaveModeSpec` class specifies how quantities related to transmission line
+    modes are computed. For example, it defines the paths for line integrals, which are used to
+    compute voltage, current, and characteristic impedance of the transmission line.
+
+    Users may supply their own voltage and current path specifications to control where these integrals
+    are evaluated. If neither voltage nor current specifications are provided, an automatic choice of
+    paths will be made based on the simulation geometry and context.
+    """
+
+    voltage_spec: Optional[tuple[Optional[VoltagePathSpecTypes], ...]] = pd.Field(
+        None,
+        title="Voltage Integration Path",
+        description="Path specification for computing the voltage associated with each mode. "
+        "The number of path specifications should equal the 'num_modes' field "
+        "in the 'ModeSpec'.",
+    )
+
+    current_spec: Optional[tuple[Optional[CurrentPathSpecTypes], ...]] = pd.Field(
+        None,
+        title="Current Integration Path",
+        description="Path specification for computing the current associated with each mode. "
+        "The number of path specifications should equal the 'num_modes' field "
+        "in the 'ModeSpec'.",
+    )
+
+    @property
+    def use_automatic_setup(self) -> bool:
+        """Whether to setup the :class:`.MicrowaveModeSpec` automatically or use the supplied
+        path specifications."""
+        return self.voltage_spec is None and self.current_spec is None
+
+    @property
+    def num_voltage_specs(self) -> Optional[int]:
+        """The number of voltage specifications supplied."""
+        if type(self.voltage_spec) is tuple:
+            return len(self.voltage_spec)
+        return None
+
+    @property
+    def num_current_specs(self) -> Optional[int]:
+        """The number of current specifications supplied."""
+        if type(self.current_spec) is tuple:
+            return len(self.current_spec)
+        return None
+
+    @pd.validator("current_spec", always=True)
+    def check_path_spec_combinations(cls, val, values):
+        """In order to define voltage/current/impedance, either a voltage or current path spec
+        must be provided for each associated mode index.
+        """
+
+        voltage_specs = values["voltage_spec"]
+        if val is None and voltage_specs is None:
+            return val
+        elif val is not None and voltage_specs is not None:
+            if len(val) != len(voltage_specs):
+                raise SetupError(
+                    f"The length of 'voltage_spec' is {len(voltage_specs)}, which is not equal to the length "
+                    f"of 'current_spec' {len(val)}. Please ensure that the same number of voltage and current "
+                    "specifications have been supplied."
+                )
+
+        if val is None:
+            current_specs = (None,) * len(voltage_specs)
+        else:
+            current_specs = val
+
+        if voltage_specs is None:
+            voltage_specs = (None,) * len(val)
+
+        for current_spec, voltage_spec, index in zip(
+            current_specs, voltage_specs, range(len(current_specs))
+        ):
+            if current_spec is None and voltage_spec is None:
+                raise SetupError(
+                    f"Both entries in 'voltage_spec' and 'current_spec' at position {index} are "
+                    "'None'. Please ensure at least one of them is a valid path specification."
+                )
+
+        return val
diff --git a/tidy3d/components/microwave/path_integral_factory.py b/tidy3d/components/microwave/path_integral_factory.py
new file mode 100644
index 000000000..a3ebf7128
--- /dev/null
+++ b/tidy3d/components/microwave/path_integral_factory.py
@@ -0,0 +1,148 @@
+"""Factory functions for creating current and voltage path integrals from path specifications."""
+
+from __future__ import annotations
+
+from typing import Union
+
+from tidy3d.components.microwave.microwave_mode_spec import MicrowaveModeSpec
+from tidy3d.components.microwave.path_spec import (
+    CompositeCurrentIntegralSpec,
+    CurrentIntegralAxisAlignedSpec,
+    CurrentPathSpecTypes,
+    CustomCurrentIntegral2DSpec,
+    CustomVoltageIntegral2DSpec,
+    VoltageIntegralAxisAlignedSpec,
+    VoltagePathSpecTypes,
+)
+from tidy3d.components.microwave.path_spec_generator import PathSpecGenerator
+from tidy3d.components.monitor import ModeMonitor, ModeSolverMonitor
+from tidy3d.components.simulation import Simulation
+from tidy3d.exceptions import SetupError, ValidationError
+from tidy3d.plugins.microwave import (
+    CompositeCurrentIntegral,
+    CurrentIntegralAxisAligned,
+    CurrentIntegralTypes,
+    CustomCurrentIntegral2D,
+    CustomVoltageIntegral2D,
+    VoltageIntegralAxisAligned,
+    VoltageIntegralTypes,
+)
+
+
+def make_voltage_integral(path_spec: VoltagePathSpecTypes) -> VoltageIntegralTypes:
+    """Create a voltage path integral from a path specification.
+
+    Parameters
+    ----------
+    path_spec : VoltagePathSpecTypes
+        Specification defining the path for voltage integration. Can be either an axis-aligned or
+        custom path specification.
+
+    Returns
+    -------
+    VoltageIntegralTypes
+        Voltage path integral instance corresponding to the provided specification type.
+    """
+    v_integral = None
+    if isinstance(path_spec, VoltageIntegralAxisAlignedSpec):
+        v_integral = VoltageIntegralAxisAligned(**path_spec.dict(exclude={"type"}))
+    elif isinstance(path_spec, CustomVoltageIntegral2DSpec):
+        v_integral = CustomVoltageIntegral2D(**path_spec.dict(exclude={"type"}))
+    else:
+        raise ValidationError(f"Unsupported voltage path specification type: {type(path_spec)}")
+    return v_integral
+
+
+def make_current_integral(path_spec: CurrentPathSpecTypes) -> CurrentIntegralTypes:
+    """Create a current path integral from a path specification.
+
+    Parameters
+    ----------
+    path_spec : CurrentPathSpecTypes
+        Specification defining the path for current integration. Can be either an axis-aligned,
+        custom, or composite path specification.
+
+    Returns
+    -------
+    CurrentIntegralTypes
+        Current path integral instance corresponding to the provided specification type.
+    """
+    i_integral = None
+    if isinstance(path_spec, CurrentIntegralAxisAlignedSpec):
+        i_integral = CurrentIntegralAxisAligned(**path_spec.dict(exclude={"type"}))
+    elif isinstance(path_spec, CustomCurrentIntegral2DSpec):
+        i_integral = CustomCurrentIntegral2D(**path_spec.dict(exclude={"type"}))
+    elif isinstance(path_spec, CompositeCurrentIntegralSpec):
+        i_integral = CompositeCurrentIntegral(**path_spec.dict(exclude={"type"}))
+    else:
+        raise ValidationError(f"Unsupported current path specification type: {type(path_spec)}")
+    return i_integral
+
+
+def make_path_integrals(
+    microwave_mode_spec: MicrowaveModeSpec,
+    monitor: Union[ModeMonitor, ModeSolverMonitor],
+    sim: Simulation,
+) -> tuple[tuple[VoltageIntegralTypes, CurrentIntegralTypes], ...]:
+    """
+    Given a ``MicrowaveModeSpec``, monitor, and simulation instance, create the voltage and current path integrals used for impedance computation.
+
+    Parameters
+    ----------
+    mw_mode_spec : MicrowaveModeSpec
+        Terminal specification containing voltage and current path specifications.
+    monitor : Union[ModeMonitor, ModeSolverMonitor]
+        The monitor for which the path integrals are being generated.
+    sim : Simulation
+        The simulation instance providing structures, grid, symmetry, and bounding box.
+
+    Returns
+    -------
+    tuple[tuple[VoltageIntegralTypes, CurrentIntegralTypes], ...]
+        Tuple containing the voltage and current path integral instances for each mode.
+
+    Raises
+    ------
+    SetupError
+        If path specifications cannot be auto-generated or path integrals cannot be constructed.
+    """
+    try:
+        v_specs = microwave_mode_spec.voltage_spec
+        i_specs = microwave_mode_spec.current_spec
+        if microwave_mode_spec.use_automatic_setup:
+            i_spec, _ = PathSpecGenerator.create_current_path_specs(
+                monitor.bounding_box,
+                sim.structures,
+                sim.grid,
+                sim.symmetry,
+                sim.bounding_box,
+                monitor.colocate,
+            )
+            i_specs = (i_spec,) * monitor.mode_spec.num_modes
+        if v_specs is None:
+            v_specs = (None,) * monitor.mode_spec.num_modes
+        if i_specs is None:
+            i_specs = (None,) * monitor.mode_spec.num_modes
+
+    except ValidationError as e:
+        raise SetupError(
+            f"Failed to auto-generate path specification for impedance calculation in monitor '{monitor.name}'."
+        ) from e
+
+    try:
+        path_integrals = []
+        for v_spec, i_spec in zip(v_specs, i_specs):
+            v_integral = None
+            i_integral = None
+            if v_spec is not None:
+                v_integral = make_voltage_integral(v_spec)
+            if i_spec is not None:
+                i_integral = make_current_integral(i_spec)
+            path_integrals.append((v_integral, i_integral))
+        path_integrals = tuple(path_integrals)
+    except Exception as e:
+        raise SetupError(
+            f"Failed to construct path integrals from the microwave mode specification for monitor '{monitor.name}'. "
+            "Please create a github issue so that the problem can be investigated."
+        ) from e
+    return path_integrals
diff --git a/tidy3d/components/microwave/path_spec.py b/tidy3d/components/microwave/path_spec.py
new file mode 100644
index 000000000..d74c6f313
--- /dev/null
+++ b/tidy3d/components/microwave/path_spec.py
@@ -0,0 +1,770 @@
+"""Module containing specifications for voltage and current path integrals."""
+
+from __future__ import annotations
+
+from abc import ABC, abstractmethod
+from typing import Literal, Optional, Union
+
+import numpy as np
+import pydantic.v1 as pd
+import shapely
+import xarray as xr
+
+from tidy3d.components.base import Tidy3dBaseModel, cached_property
+from tidy3d.components.geometry.base import Box, Geometry
+from tidy3d.components.types import (
+    ArrayFloat2D,
+    Ax,
+    Axis,
+    Bound,
+    Coordinate,
+    Coordinate2D,
+    Direction,
+)
+from tidy3d.components.validators import assert_line, assert_plane
+from tidy3d.components.viz import add_ax_if_none
+from tidy3d.constants import MICROMETER, fp_eps
+from tidy3d.exceptions import SetupError, Tidy3dError
+from tidy3d.log import log
+
+from .viz import (
+    ARROW_CURRENT,
+    plot_params_current_path,
+    plot_params_voltage_minus,
+    plot_params_voltage_path,
+    plot_params_voltage_plus,
+)
+
+
+class AbstractAxesRH(Tidy3dBaseModel, ABC):
+    """Represents an axis-aligned right-handed coordinate system with one axis preferred.
+    Typically `main_axis` would refer to the normal axis of a plane.
+    """
+
+    @cached_property
+    @abstractmethod
+    def main_axis(self) -> Axis:
+        """Get the preferred axis."""
+
+    @cached_property
+    def remaining_axes(self) -> tuple[Axis, Axis]:
+        """Get in-plane axes, ordered to maintain a right-handed coordinate system."""
+        axes: list[Axis] = [0, 1, 2]
+        axes.pop(self.main_axis)
+        if self.main_axis == 1:
+            return (axes[1], axes[0])
+        else:
+            return (axes[0], axes[1])
+
+    @cached_property
+    def remaining_dims(self) -> tuple[str, str]:
+        """Get in-plane dimensions, ordered to maintain a right-handed coordinate system."""
+        dim1 = "xyz"[self.remaining_axes[0]]
+        dim2 = "xyz"[self.remaining_axes[1]]
+        return (dim1, dim2)
+
+    @cached_property
+    def local_dims(self) -> tuple[str, str, str]:
+        """Get in-plane dimensions with in-plane dims first, followed by the `main_axis` dimension."""
+        dim3 = "xyz"[self.main_axis]
+        return self.remaining_dims + tuple(dim3)
+
+    @pd.root_validator(pre=False)
+    def _warn_rf_license(cls, values):
+        log.warning(
+            "ℹ️ ⚠️ RF simulations are subject to new license requirements in the future. You have instantiated at least one RF-specific component.",
+            log_once=True,
+        )
+        return values
+
+
+class AxisAlignedPathIntegralSpec(AbstractAxesRH, Box):
+    """Class for defining the simplest type of path integral, which is aligned with Cartesian axes."""
+
+    _line_validator = assert_line()
+
+    extrapolate_to_endpoints: bool = pd.Field(
+        False,
+        title="Extrapolate to Endpoints",
+        description="If the endpoints of the path integral terminate at or near a material interface, "
+        "the field is likely discontinuous. When this field is ``True``, fields that are outside and on the bounds "
+        "of the integral are ignored. Should be enabled when computing voltage between two conductors.",
+    )
+
+    snap_path_to_grid: bool = pd.Field(
+        False,
+        title="Snap Path to Grid",
+        description="It might be desirable to integrate exactly along the Yee grid associated with "
+        "a field. When this field is ``True``, the integration path will be snapped to the grid.",
+    )
+
+    @cached_property
+    def main_axis(self) -> Axis:
+        """Axis for performing integration."""
+        for index, value in enumerate(self.size):
+            if value != 0:
+                return index
+        raise Tidy3dError("Failed to identify axis.")
+
+    def _vertices_2D(self, axis: Axis) -> tuple[Coordinate2D, Coordinate2D]:
+        """Returns the two vertices of this path in the plane defined by ``axis``."""
+        min = self.bounds[0]
+        max = self.bounds[1]
+        _, min = Box.pop_axis(min, axis)
+        _, max = Box.pop_axis(max, axis)
+
+        u = [min[0], max[0]]
+        v = [min[1], max[1]]
+        return (u, v)
+
+
+class VoltageIntegralAxisAlignedSpec(AxisAlignedPathIntegralSpec):
+    """Class for specifying the voltage calculation between two points defined by an axis-aligned line."""
+
+    sign: Direction = pd.Field(
+        ...,
+        title="Direction of Path Integral",
+        description="Positive indicates V=Vb-Va where position b has a larger coordinate along the axis of integration.",
+    )
+
+    @staticmethod
+    def from_terminal_positions(
+        plus_terminal: float,
+        minus_terminal: float,
+        x: Optional[float] = None,
+        y: Optional[float] = None,
+        z: Optional[float] = None,
+        extrapolate_to_endpoints: bool = True,
+        snap_path_to_grid: bool = True,
+    ) -> VoltageIntegralAxisAlignedSpec:
+        """Helper to create a :class:`VoltageIntegralAxisAlignedSpec` from two coordinates that
+        define a line and two positions indicating the endpoints of the path integral.
+
+        Parameters
+        ----------
+        plus_terminal : float
+            Position along the voltage axis of the positive terminal.
+        minus_terminal : float
+            Position along the voltage axis of the negative terminal.
+        x : float = None
+            Position in x direction, only two of x,y,z can be specified to define line.
+        y : float = None
+            Position in y direction, only two of x,y,z can be specified to define line.
+        z : float = None
+            Position in z direction, only two of x,y,z can be specified to define line.
+        extrapolate_to_endpoints: bool = True
+            Passed directly to :class:`VoltageIntegralAxisAlignedSpec`
+        snap_path_to_grid: bool = True
+            Passed directly to :class:`VoltageIntegralAxisAlignedSpec`
+
+        Returns
+        -------
+        VoltageIntegralAxisAlignedSpec
+            The created path integral for computing voltage between the two terminals.
+        """
+        axis_positions = Geometry.parse_two_xyz_kwargs(x=x, y=y, z=z)
+        # Calculate center and size of the future box
+        midpoint = (plus_terminal + minus_terminal) / 2
+        length = np.abs(plus_terminal - minus_terminal)
+        center = [midpoint, midpoint, midpoint]
+        size = [length, length, length]
+        for axis, position in axis_positions:
+            size[axis] = 0
+            center[axis] = position
+
+        direction = "+"
+        if plus_terminal < minus_terminal:
+            direction = "-"
+
+        return VoltageIntegralAxisAlignedSpec(
+            center=center,
+            size=size,
+            extrapolate_to_endpoints=extrapolate_to_endpoints,
+            snap_path_to_grid=snap_path_to_grid,
+            sign=direction,
+        )
+
+    @add_ax_if_none
+    def plot(
+        self,
+        x: Optional[float] = None,
+        y: Optional[float] = None,
+        z: Optional[float] = None,
+        ax: Ax = None,
+        **path_kwargs,
+    ) -> Ax:
+        """Plot path integral at single (x,y,z) coordinate.
+
+        Parameters
+        ----------
+        x : float = None
+            Position of plane in x direction, only one of x,y,z can be specified to define plane.
+        y : float = None
+            Position of plane in y direction, only one of x,y,z can be specified to define plane.
+        z : float = None
+            Position of plane in z direction, only one of x,y,z can be specified to define plane.
+        ax : matplotlib.axes._subplots.Axes = None
+            Matplotlib axes to plot on, if not specified, one is created.
+        **path_kwargs
+            Optional keyword arguments passed to the matplotlib plotting of the line.
+            For details on accepted values, refer to
+            `Matplotlib's documentation <https://tinyurl.com/36marrat>`_.
+
+        Returns
+        -------
+        matplotlib.axes._subplots.Axes
+            The supplied or created matplotlib axes.
+        """
+        axis, position = self.parse_xyz_kwargs(x=x, y=y, z=z)
+        if axis == self.main_axis or not np.isclose(position, self.center[axis], rtol=fp_eps):
+            return ax
+
+        (xs, ys) = self._vertices_2D(axis)
+        # Plot the path
+        plot_params = plot_params_voltage_path.include_kwargs(**path_kwargs)
+        plot_kwargs = plot_params.to_kwargs()
+        ax.plot(xs, ys, markevery=[0, -1], **plot_kwargs)
+
+        # Plot special end points
+        end_kwargs = plot_params_voltage_plus.include_kwargs(**path_kwargs).to_kwargs()
+        start_kwargs = plot_params_voltage_minus.include_kwargs(**path_kwargs).to_kwargs()
+
+        if self.sign == "-":
+            start_kwargs, end_kwargs = end_kwargs, start_kwargs
+
+        ax.plot(xs[0], ys[0], **start_kwargs)
+        ax.plot(xs[1], ys[1], **end_kwargs)
+        return ax
+
+
+class CurrentIntegralAxisAlignedSpec(AbstractAxesRH, Box):
+    """Class for specifying the computation of conduction current via Ampère's circuital law on an axis-aligned loop."""
+
+    _plane_validator = assert_plane()
+
+    sign: Direction = pd.Field(
+        ...,
+        title="Direction of Contour Integral",
+        description="Positive indicates current flowing in the positive normal axis direction.",
+    )
+
+    extrapolate_to_endpoints: bool = pd.Field(
+        False,
+        title="Extrapolate to Endpoints",
+        description="This parameter is passed to :class:`AxisAlignedPathIntegral` objects when computing the contour integral.",
+    )
+
+    snap_contour_to_grid: bool = pd.Field(
+        False,
+        title="Snap Contour to Grid",
+        description="This parameter is passed to :class:`AxisAlignedPathIntegral` objects when computing the contour integral.",
+    )
+
+    @cached_property
+    def main_axis(self) -> Axis:
+        """Axis normal to loop"""
+        for index, value in enumerate(self.size):
+            if value == 0:
+                return index
+        raise Tidy3dError("Failed to identify axis.")
+
+    def _to_path_integral_specs(
+        self, h_horizontal=None, h_vertical=None
+    ) -> tuple[AxisAlignedPathIntegralSpec, ...]:
+        """Returns four ``AxisAlignedPathIntegralSpec`` instances, which represent a contour
+        integral around the surface defined by ``self.size``."""
+        ax1 = self.remaining_axes[0]
+        ax2 = self.remaining_axes[1]
+
+        horizontal_passed = h_horizontal is not None
+        vertical_passed = h_vertical is not None
+        if self.snap_contour_to_grid and horizontal_passed and vertical_passed:
+            (coord1, coord2) = self.remaining_dims
+
+            # Locations where horizontal paths will be snapped
+            v_bounds = [
+                self.center[ax2] - self.size[ax2] / 2,
+                self.center[ax2] + self.size[ax2] / 2,
+            ]
+            h_snaps = h_horizontal.sel({coord2: v_bounds}, method="nearest").coords[coord2].values
+            # Locations where vertical paths will be snapped
+            h_bounds = [
+                self.center[ax1] - self.size[ax1] / 2,
+                self.center[ax1] + self.size[ax1] / 2,
+            ]
+            v_snaps = h_vertical.sel({coord1: h_bounds}, method="nearest").coords[coord1].values
+
+            bottom_bound = h_snaps[0]
+            top_bound = h_snaps[1]
+            left_bound = v_snaps[0]
+            right_bound = v_snaps[1]
+        else:
+            bottom_bound = self.bounds[0][ax2]
+            top_bound = self.bounds[1][ax2]
+            left_bound = self.bounds[0][ax1]
+            right_bound = self.bounds[1][ax1]
+
+        # Horizontal paths
+        path_size = list(self.size)
+        path_size[ax1] = right_bound - left_bound
+        path_size[ax2] = 0
+        path_center = list(self.center)
+        path_center[ax2] = bottom_bound
+
+        bottom = AxisAlignedPathIntegralSpec(
+            center=path_center,
+            size=path_size,
+            extrapolate_to_endpoints=self.extrapolate_to_endpoints,
+            snap_path_to_grid=self.snap_contour_to_grid,
+        )
+        path_center[ax2] = top_bound
+        top = AxisAlignedPathIntegralSpec(
+            center=path_center,
+            size=path_size,
+            extrapolate_to_endpoints=self.extrapolate_to_endpoints,
+            snap_path_to_grid=self.snap_contour_to_grid,
+        )
+
+        # Vertical paths
+        path_size = list(self.size)
+        path_size[ax1] = 0
+        path_size[ax2] = top_bound - bottom_bound
+        path_center = list(self.center)
+
+        path_center[ax1] = left_bound
+        left = AxisAlignedPathIntegralSpec(
+            center=path_center,
+            size=path_size,
+            extrapolate_to_endpoints=self.extrapolate_to_endpoints,
+            snap_path_to_grid=self.snap_contour_to_grid,
+        )
+        path_center[ax1] = right_bound
+        right = AxisAlignedPathIntegralSpec(
+            center=path_center,
+            size=path_size,
+            extrapolate_to_endpoints=self.extrapolate_to_endpoints,
+            snap_path_to_grid=self.snap_contour_to_grid,
+        )
+
+        return (bottom, right, top, left)
+
+    @add_ax_if_none
+    def plot(
+        self,
+        x: Optional[float] = None,
+        y: Optional[float] = None,
+        z: Optional[float] = None,
+        ax: Ax = None,
+        **path_kwargs,
+    ) -> Ax:
+        """Plot path integral at single (x,y,z) coordinate.
+
+        Parameters
+        ----------
+        x : float = None
+            Position of plane in x direction, only one of x,y,z can be specified to define plane.
+        y : float = None
+            Position of plane in y direction, only one of x,y,z can be specified to define plane.
+        z : float = None
+            Position of plane in z direction, only one of x,y,z can be specified to define plane.
+        ax : matplotlib.axes._subplots.Axes = None
+            Matplotlib axes to plot on, if not specified, one is created.
+        **path_kwargs
+            Optional keyword arguments passed to the matplotlib plotting of the line.
+            For details on accepted values, refer to
+            `Matplotlib's documentation <https://tinyurl.com/36marrat>`_.
+
+        Returns
+        -------
+        matplotlib.axes._subplots.Axes
+            The supplied or created matplotlib axes.
+        """
+        axis, position = self.parse_xyz_kwargs(x=x, y=y, z=z)
+        if axis != self.main_axis or not np.isclose(position, self.center[axis], rtol=fp_eps):
+            return ax
+
+        plot_params = plot_params_current_path.include_kwargs(**path_kwargs)
+        plot_kwargs = plot_params.to_kwargs()
+        path_integrals = self._to_path_integral_specs()
+        # Plot the path
+        for path in path_integrals:
+            (xs, ys) = path._vertices_2D(axis)
+            ax.plot(xs, ys, **plot_kwargs)
+
+        (ax1, ax2) = self.remaining_axes
+
+        # Add arrow to bottom path, unless right path is longer
+        arrow_path = path_integrals[0]
+        if self.size[ax2] > self.size[ax1]:
+            arrow_path = path_integrals[1]
+
+        (xs, ys) = arrow_path._vertices_2D(axis)
+        X = (xs[0] + xs[1]) / 2
+        Y = (ys[0] + ys[1]) / 2
+        center = np.array([X, Y])
+        dx = xs[1] - xs[0]
+        dy = ys[1] - ys[0]
+        direction = np.array([dx, dy])
+        segment_length = np.linalg.norm(direction)
+        unit_dir = direction / segment_length
+
+        # Change direction of arrow depending on sign of current definition
+        if self.sign == "-":
+            unit_dir *= -1.0
+        # Change direction of arrow when the "y" axis is dropped,
+        # since the plotted coordinate system will be left-handed (x, z)
+        if self.main_axis == 1:
+            unit_dir *= -1.0
+
+        start = center - unit_dir * segment_length
+        end = center
+        ax.annotate(
+            "",
+            xytext=(start[0], start[1]),
+            xy=(end[0], end[1]),
+            arrowprops=ARROW_CURRENT,
+        )
+        return ax
+
+
+class CustomPathIntegral2DSpec(AbstractAxesRH):
+    """Class for specifying a custom path integral defined as a curve on an axis-aligned plane.
+
+    Notes
+    -----
+
+    Given a set of vertices :math:`\\vec{r}_i`, this class approximates path integrals over
+    vector fields of the form :math:`\\int{\\vec{F} \\cdot \\vec{dl}}`
+    as :math:`\\sum_i{\\vec{F}(\\vec{r}_i) \\cdot \\vec{dl}_i}`,
+    where the differential length :math:`\\vec{dl}` is approximated using central differences
+    :math:`\\vec{dl}_i = \\frac{\\vec{r}_{i+1} - \\vec{r}_{i-1}}{2}`.
+    If the path is not closed, forward and backward differences are used at the endpoints.
+    """
+
+    axis: Axis = pd.Field(
+        2, title="Axis", description="Specifies dimension of the planar axis (0,1,2) -> (x,y,z)."
+    )
+
+    position: float = pd.Field(
+        ...,
+        title="Position",
+        description="Position of the plane along the ``axis``.",
+    )
+
+    vertices: ArrayFloat2D = pd.Field(
+        ...,
+        title="Vertices",
+        description="List of (d1, d2) defining the 2 dimensional positions of the path. "
+        "The index of dimension should be in the ascending order, which means "
+        "if the axis corresponds with ``y``, the coordinates of the vertices should be (x, z). "
+        "If you wish to indicate a closed contour, the final vertex should be made "
+        "equal to the first vertex, i.e., ``vertices[-1] == vertices[0]``",
+        units=MICROMETER,
+    )
+
+    @staticmethod
+    def _compute_dl_component(coord_array: xr.DataArray, closed_contour=False) -> np.array:
+        """Computes the differential length element along the integration path."""
+        dl = np.gradient(coord_array)
+        if closed_contour:
+            # If the contour is closed, we can use central difference on the starting/end point
+            # which will be more accurate than the default forward/backward choice in np.gradient
+            grad_end = np.gradient([coord_array[-2], coord_array[0], coord_array[1]])
+            dl[0] = dl[-1] = grad_end[1]
+        return dl
+
+    @classmethod
+    def from_circular_path(
+        cls, center: Coordinate, radius: float, num_points: int, normal_axis: Axis, clockwise: bool
+    ) -> CustomPathIntegral2DSpec:
+        """Creates a ``CustomPathIntegral2DSpec`` from a circular path given a desired number of points
+        along the perimeter.
+
+        Parameters
+        ----------
+        center : Coordinate
+            The center of the circle.
+        radius : float
+            The radius of the circle.
+        num_points : int
+            The number of equidistant points to use along the perimeter of the circle.
+        normal_axis : Axis
+            The axis normal to the defined circle.
+        clockwise : bool
+            When ``True``, the points will be ordered clockwise with respect to the positive
+            direction of the ``normal_axis``.
+
+        Returns
+        -------
+        :class:`.CustomPathIntegral2DSpec`
+            A path integral defined on a circular path.
+        """
+
+        def generate_circle_coordinates(radius: float, num_points: int, clockwise: bool):
+            """Helper for generating x,y vertices around a circle in the local coordinate frame."""
+            sign = 1.0
+            if clockwise:
+                sign = -1.0
+            angles = np.linspace(0, sign * 2 * np.pi, num_points, endpoint=True)
+            xt = radius * np.cos(angles)
+            yt = radius * np.sin(angles)
+            return (xt, yt)
+
+        # Get transverse axes
+        normal_center, trans_center = Geometry.pop_axis(center, normal_axis)
+
+        # These x,y coordinates in the local coordinate frame
+        if normal_axis == 1:
+            # Handle special case when y is the axis that is popped
+            clockwise = not clockwise
+        xt, yt = generate_circle_coordinates(radius, num_points, clockwise)
+        xt += trans_center[0]
+        yt += trans_center[1]
+        circle_vertices = np.column_stack((xt, yt))
+        # Close the contour exactly
+        circle_vertices[-1, :] = circle_vertices[0, :]
+        return cls(axis=normal_axis, position=normal_center, vertices=circle_vertices)
+
+    @cached_property
+    def is_closed_contour(self) -> bool:
+        """Returns ``true`` when the first vertex equals the last vertex."""
+        return np.isclose(
+            self.vertices[0, :],
+            self.vertices[-1, :],
+            rtol=fp_eps,
+            atol=np.finfo(np.float32).smallest_normal,
+        ).all()
+
+    @cached_property
+    def main_axis(self) -> Axis:
+        """Axis for performing integration."""
+        return self.axis
+
+    @pd.validator("vertices", always=True)
+    def _correct_shape(cls, val):
+        """Makes sure vertices size is correct."""
+        # overall shape of vertices
+        if val.shape[1] != 2:
+            raise SetupError(
+                "'CustomPathIntegral2DSpec.vertices' must be a 2 dimensional array shaped (N, 2). "
+                f"Given array with shape of '{val.shape}'."
+            )
+        return val
+
+    @cached_property
+    def bounds(self) -> Bound:
+        """Helper to get the geometric bounding box of the path integral."""
+        path_min = np.amin(self.vertices, axis=0)
+        path_max = np.amax(self.vertices, axis=0)
+        min_bound = Geometry.unpop_axis(self.position, path_min, self.axis)
+        max_bound = Geometry.unpop_axis(self.position, path_max, self.axis)
+        return (min_bound, max_bound)
+
+    @cached_property
+    def sign(self) -> Direction:
+        """Uses the ordering of the vertices to determine the direction of the current flow."""
+        linestr = shapely.LineString(coordinates=self.vertices)
+        is_ccw = shapely.is_ccw(linestr)
+        # Invert statement when the vertices are given as (x, z)
+        if self.axis == 1:
+            is_ccw = not is_ccw
+        if is_ccw:
+            return "+"
+        else:
+            return "-"
+
+
+class CustomVoltageIntegral2DSpec(CustomPathIntegral2DSpec):
+    """Class for specfying the computation of voltage between two points defined by a custom path.
+    Computed voltage is :math:`V=V_b-V_a`, where position b is the final vertex in the supplied path.
+
+    Notes
+    -----
+
+    Use :class:`.VoltageIntegralAxisAlignedSpec` if possible, since interpolation
+    near conductors will not be accurate.
+
+    .. TODO Improve by including extrapolate_to_endpoints field, non-trivial extension."""
+
+    @add_ax_if_none
+    def plot(
+        self,
+        x: Optional[float] = None,
+        y: Optional[float] = None,
+        z: Optional[float] = None,
+        ax: Ax = None,
+        **path_kwargs,
+    ) -> Ax:
+        """Plot path integral at single (x,y,z) coordinate.
+
+        Parameters
+        ----------
+        x : float = None
+            Position of plane in x direction, only one of x,y,z can be specified to define plane.
+        y : float = None
+            Position of plane in y direction, only one of x,y,z can be specified to define plane.
+        z : float = None
+            Position of plane in z direction, only one of x,y,z can be specified to define plane.
+        ax : matplotlib.axes._subplots.Axes = None
+            Matplotlib axes to plot on, if not specified, one is created.
+        **path_kwargs
+            Optional keyword arguments passed to the matplotlib plotting of the line.
+            For details on accepted values, refer to
+            `Matplotlib's documentation <https://tinyurl.com/36marrat>`_.
+
+        Returns
+        -------
+        matplotlib.axes._subplots.Axes
+            The supplied or created matplotlib axes.
+        """
+        axis, position = Geometry.parse_xyz_kwargs(x=x, y=y, z=z)
+        if axis != self.main_axis or not np.isclose(position, self.position, rtol=fp_eps):
+            return ax
+
+        plot_params = plot_params_voltage_path.include_kwargs(**path_kwargs)
+        plot_kwargs = plot_params.to_kwargs()
+        xs = self.vertices[:, 0]
+        ys = self.vertices[:, 1]
+        ax.plot(xs, ys, markevery=[0, -1], **plot_kwargs)
+
+        # Plot special end points
+        end_kwargs = plot_params_voltage_plus.include_kwargs(**path_kwargs).to_kwargs()
+        start_kwargs = plot_params_voltage_minus.include_kwargs(**path_kwargs).to_kwargs()
+        ax.plot(xs[0], ys[0], **start_kwargs)
+        ax.plot(xs[-1], ys[-1], **end_kwargs)
+
+        return ax
+
+
+class CustomCurrentIntegral2DSpec(CustomPathIntegral2DSpec):
+    """Class for specifying the computation of conduction current via Ampère's circuital law on a custom path.
+    To compute the current flowing in the positive ``axis`` direction, the vertices should be
+    ordered in a counterclockwise direction."""
+
+    @add_ax_if_none
+    def plot(
+        self,
+        x: Optional[float] = None,
+        y: Optional[float] = None,
+        z: Optional[float] = None,
+        ax: Ax = None,
+        **path_kwargs,
+    ) -> Ax:
+        """Plot path integral at single (x,y,z) coordinate.
+
+        Parameters
+        ----------
+        x : float = None
+            Position of plane in x direction, only one of x,y,z can be specified to define plane.
+        y : float = None
+            Position of plane in y direction, only one of x,y,z can be specified to define plane.
+        z : float = None
+            Position of plane in z direction, only one of x,y,z can be specified to define plane.
+        ax : matplotlib.axes._subplots.Axes = None
+            Matplotlib axes to plot on, if not specified, one is created.
+        **path_kwargs
+            Optional keyword arguments passed to the matplotlib plotting of the line.
+            For details on accepted values, refer to
+            `Matplotlib's documentation <https://tinyurl.com/36marrat>`_.
+
+        Returns
+        -------
+        matplotlib.axes._subplots.Axes
+            The supplied or created matplotlib axes.
+        """
+        axis, position = Geometry.parse_xyz_kwargs(x=x, y=y, z=z)
+        if axis != self.main_axis or not np.isclose(position, self.position, rtol=fp_eps):
+            return ax
+
+        plot_params = plot_params_current_path.include_kwargs(**path_kwargs)
+        plot_kwargs = plot_params.to_kwargs()
+        xs = self.vertices[:, 0]
+        ys = self.vertices[:, 1]
+        ax.plot(xs, ys, **plot_kwargs)
+
+        # Add arrow at start of contour
+        ax.annotate(
+            "",
+            xytext=(xs[0], ys[0]),
+            xy=(xs[1], ys[1]),
+            arrowprops=ARROW_CURRENT,
+        )
+        return ax
+
+    @cached_property
+    def sign(self) -> Direction:
+        """Uses the ordering of the vertices to determine the direction of the current flow."""
+        linestr = shapely.LineString(coordinates=self.vertices)
+        is_ccw = shapely.is_ccw(linestr)
+        # Invert statement when the vertices are given as (x, z)
+        if self.axis == 1:
+            is_ccw = not is_ccw
+        if is_ccw:
+            return "+"
+        else:
+            return "-"
+
+
+class CompositeCurrentIntegralSpec(Box):
+    """Specification for a composite current integral.
+
+    This class is used to set up a CompositeCurrentIntegral, which combines
+    multiple current integrals. It does not perform any integration itself.
+    """
+
+    path_specs: list[Union[CurrentIntegralAxisAlignedSpec, CustomCurrentIntegral2DSpec]] = pd.Field(
+        ...,
+        title="Path Specifications",
+        description="Definition of the disjoint path specifications for each isolated contour integral.",
+    )
+
+    sum_spec: Literal["sum", "split"] = pd.Field(
+        ...,
+        title="Sum Specification",
+        description="Determines the method used to combine the currents calculated by the different "
+        "current integrals defined by ``path_specs``. ``sum`` simply adds all currents, while ``split`` "
+        "keeps contributions with opposite phase separate, which allows for isolating the current "
+        "flowing in opposite directions. In ``split`` version, the current returned is the maximum "
+        "of the two contributions.",
+    )
+
+    @add_ax_if_none
+    def plot(
+        self,
+        x: Optional[float] = None,
+        y: Optional[float] = None,
+        z: Optional[float] = None,
+        ax: Ax = None,
+        **path_kwargs,
+    ) -> Ax:
+        """Plot path integral at single (x,y,z) coordinate.
+
+        Parameters
+        ----------
+        x : float = None
+            Position of plane in x direction, only one of x,y,z can be specified to define plane.
+        y : float = None
+            Position of plane in y direction, only one of x,y,z can be specified to define plane.
+        z : float = None
+            Position of plane in z direction, only one of x,y,z can be specified to define plane.
+        ax : matplotlib.axes._subplots.Axes = None
+            Matplotlib axes to plot on, if not specified, one is created.
+        **path_kwargs
+            Optional keyword arguments passed to the matplotlib plotting of the line.
+            For details on accepted values, refer to
+            `Matplotlib's documentation <https://tinyurl.com/36marrat>`_.
+
+        Returns
+        -------
+        matplotlib.axes._subplots.Axes
+            The supplied or created matplotlib axes.
+        """
+        for path_spec in self.path_specs:
+            ax = path_spec.plot(x=x, y=y, z=z, ax=ax, **path_kwargs)
+        return ax
+
+
+VoltagePathSpecTypes = Union[VoltageIntegralAxisAlignedSpec, CustomVoltageIntegral2DSpec]
+CurrentPathSpecTypes = Union[
+    CurrentIntegralAxisAlignedSpec, CustomCurrentIntegral2DSpec, CompositeCurrentIntegralSpec
+]
diff --git a/tidy3d/components/microwave/path_spec_generator.py b/tidy3d/components/microwave/path_spec_generator.py
new file mode 100644
index 000000000..f286a4d13
--- /dev/null
+++ b/tidy3d/components/microwave/path_spec_generator.py
@@ -0,0 +1,326 @@
+"""Module for automatic determination of voltage and current integration specifications."""
+
+from __future__ import annotations
+
+from itertools import chain
+from math import isclose
+
+import shapely
+from shapely.geometry import (
+    LineString,
+    Polygon,
+)
+
+from tidy3d.components.base import Tidy3dBaseModel
+from tidy3d.components.geometry.base import Box, Geometry
+from tidy3d.components.geometry.utils import (
+    SnapBehavior,
+    SnapLocation,
+    SnappingSpec,
+    flatten_shapely_geometries,
+    merging_geometries_on_plane,
+    snap_box_to_grid,
+)
+from tidy3d.components.grid.grid import Grid
+from tidy3d.components.structure import Structure
+from tidy3d.components.types import Axis, Coordinate, Shapely, Symmetry
+from tidy3d.exceptions import ValidationError
+
+from .path_spec import CompositeCurrentIntegralSpec, CurrentIntegralAxisAlignedSpec
+
+
+class PathSpecGenerator(Tidy3dBaseModel):
+    """Automatically determines current integration paths based on structure geometry.
+
+    This class analyzes the geometry of conductors in a simulation cross-section to determine
+    appropriate paths for computing current line integrals. These paths are typically used
+    for setting up a :class:`.MicrowaveModeSpec` automatically.
+
+    The paths are chosen by:
+    1. Finding and merging conductor surfaces that intersect with the plane
+    2. Creating current paths that enclose isolated conductors
+    """
+
+    @staticmethod
+    def _create_snap_spec(normal_axis: Axis, field_data_colocated: bool) -> SnappingSpec:
+        """Creates snapping specification for grid alignment."""
+        behavior = [SnapBehavior.Expand] * 3
+        location = [SnapLocation.Center] * 3
+        behavior[normal_axis] = SnapBehavior.Off
+        # To avoid interpolated H field near metal surface
+        margin = (2, 2, 2) if field_data_colocated else (0, 0, 0)
+        return SnappingSpec(location=location, behavior=behavior, margin=margin)
+
+    @staticmethod
+    def _create_port_boundary(mode_plane: Box, normal_axis: Axis) -> shapely.LineString:
+        """Creates a Shapely LineString representing the port boundary."""
+        _, min_b = Geometry.pop_axis(mode_plane.bounds[0], normal_axis)
+        _, max_b = Geometry.pop_axis(mode_plane.bounds[1], normal_axis)
+        port = Geometry.make_shapely_box(min_b[0], min_b[1], max_b[0], max_b[1])
+        return shapely.LineString(port.exterior)
+
+    @staticmethod
+    def _create_path_spec(box_snapped: Box) -> CurrentIntegralAxisAlignedSpec:
+        """Creates a current integral specification from a snapped box."""
+        return CurrentIntegralAxisAlignedSpec(
+            center=box_snapped.center,
+            size=box_snapped.size,
+            sign="+",
+            extrapolate_to_endpoints=False,
+            snap_contour_to_grid=True,
+        )
+
+    @staticmethod
+    def _get_isolated_conductors_as_shapely(
+        plane: Box, structures: list[Structure]
+    ) -> list[Shapely]:
+        """Find and merge all conductor structures that intersect the given plane.
+
+        Parameters
+        ----------
+        plane : Box
+            The plane to check for conductor intersections
+        structures : list[Structure]
+            List of all simulation structures to analyze
+
+        Returns
+        -------
+        list[Shapely]
+            List of merged conductor geometries as Shapely Polygons and LineStrings
+            that intersect with the given plane
+        """
+        from tidy3d.components.medium import LossyMetalMedium, Medium
+
+        def is_conductor(med: Medium) -> bool:
+            return med.is_pec or isinstance(med, LossyMetalMedium)
+
+        geometry_list = [structure.geometry for structure in structures]
+        # For metal, we don't distinguish between LossyMetal and PEC,
+        # so they'll be merged to PEC. Other materials are considered as dielectric.
+        prop_list = [is_conductor(structure.medium) for structure in structures]
+        # merge geometries
+        geos = merging_geometries_on_plane(geometry_list, plane, prop_list)
+        conductor_geos = [item[1] for item in geos if item[0]]
+        shapely_list = flatten_shapely_geometries(conductor_geos, keep_types=(Polygon, LineString))
+        return shapely_list
+
+    @staticmethod
+    def create_current_path_specs(
+        mode_plane: Box,
+        structures: list[Structure],
+        grid: Grid,
+        symmetry: tuple[Symmetry, Symmetry, Symmetry],
+        sim_box: Box,
+        field_data_colocated: bool = False,
+    ) -> tuple[CompositeCurrentIntegralSpec, list[Shapely]]:
+        """Creates path specifications for path integrals that encompass each isolated conductor
+        in the mode plane.
+
+        This method identifies isolated conductor geometries in the given plane and creates
+        current paths that enclose each conductor. The paths are snapped to the simulation grid
+        to ensure alignment with field data.
+
+        Parameters
+        ----------
+        mode_plane : Box
+            The cross-sectional plane where current paths are determined.
+        structures : list[Structure]
+            List of structures in the simulation.
+        grid : Grid
+            Simulation grid for snapping paths.
+        symmetry : tuple[Symmetry, Symmetry, Symmetry]
+            Symmetry conditions for the simulation in (x, y, z) directions.
+        sim_box : Box
+            Simulation domain box used for boundary conditions.
+        field_data_colocated : bool
+            Whether field data is colocated with grid points.
+
+        Returns
+        -------
+            tuple[CompositeCurrentIntegralSpec, list[Shapely]]
+                Composite path specification and list of merged conductor geometries.
+        """
+
+        def bounding_box_from_shapely(
+            geom: Shapely, normal_axis: Axis, normal_center: float
+        ) -> Box:
+            """Helper to convert the shapely geometry bounds to a Box."""
+            bounds = geom.bounds
+            rmin = Geometry.unpop_axis(normal_center, (bounds[0], bounds[1]), normal_axis)
+            rmax = Geometry.unpop_axis(normal_center, (bounds[2], bounds[3]), normal_axis)
+            return Box.from_bounds(rmin, rmax)
+
+        # Simulation box restricted by symmetry
+        min_sym_bounds, max_sym_bounds = sim_box.bounds
+
+        min_sym_bounds = list(min_sym_bounds)
+        for dim in range(3):
+            if symmetry[dim] != 0:
+                min_sym_bounds[dim] = sim_box.center[dim]
+        sym_sim_bounds = (tuple(min_sym_bounds), max_sym_bounds)
+
+        conductor_polygons = PathSpecGenerator._get_isolated_conductors_as_shapely(
+            mode_plane, structures
+        )
+        normal_axis = mode_plane.size.index(0.0)
+        # Get desired snapping behavior of box enclosed conductors.
+        # Ideally, just large enough to coincide with the H field positions outside of the conductor.
+        # So a half grid cell, when the metal boundary is coincident with grid boundaries.
+        snap_spec = PathSpecGenerator._create_snap_spec(normal_axis, field_data_colocated)
+
+        shapely_port_boundary = PathSpecGenerator._create_port_boundary(mode_plane, normal_axis)
+
+        current_integral_specs = []
+        for shape in conductor_polygons:
+            if not shapely_port_boundary.touches(shape):
+                box = bounding_box_from_shapely(shape, normal_axis, mode_plane.center[normal_axis])
+                intersected_box = Geometry.bounds_intersection(box.bounds, sym_sim_bounds)
+                if any(
+                    rmax - rmin < 0 for rmin, rmax in zip(intersected_box[0], intersected_box[1])
+                ):
+                    continue
+                box = Box.from_bounds(rmin=intersected_box[0], rmax=intersected_box[1])
+                boxes = PathSpecGenerator._apply_symmetry(
+                    symmetry, sim_box.center, normal_axis, box
+                )
+                for box in boxes:
+                    box_snapped = snap_box_to_grid(grid, box, snap_spec)
+                    path_spec = PathSpecGenerator._create_path_spec(box_snapped)
+                    current_integral_specs.append(path_spec)
+
+        for path_spec in current_integral_specs:
+            if PathSpecGenerator._check_path_intersects_with_conductors(
+                conductor_polygons, path_spec
+            ):
+                raise ValidationError(
+                    "Failed to automatically generate path specification. "
+                    "Please create a github issue so that the problem can be investigated. "
+                    "In the meantime, please provide an explicit path specification for this structure."
+                )
+
+        path_spec = CompositeCurrentIntegralSpec(
+            center=mode_plane.center,
+            size=mode_plane.size,
+            path_specs=current_integral_specs,
+            sum_spec="split",
+        )
+        return path_spec, conductor_polygons
+
+    @staticmethod
+    def _check_path_intersects_with_conductors(shapely_list: Shapely, path: Box) -> bool:
+        """Makes sure the path of the integral does not intersect with conductor shapes.
+
+        Parameters
+        ----------
+        shapely_list : Shapely
+            Merged conductor geometries, expected to be Polygons or LineStrings.
+        path : Box
+            Box corresponding with path specification.
+
+        Returns
+        -------
+            bool: True if the path intersects with any conductor geometry, False otherwise.
+        """
+        normal_axis = path.size.index(0.0)
+        min_b, max_b = path.bounds
+        _, min_b = Geometry.pop_axis(min_b, normal_axis)
+        _, max_b = Geometry.pop_axis(max_b, normal_axis)
+        path_shapely = shapely.box(min_b[0], min_b[1], max_b[0], max_b[1])
+        for shapely_geo in shapely_list:
+            if path_shapely.intersects(shapely_geo) and not path_shapely.contains(shapely_geo):
+                return True
+        return False
+
+    @staticmethod
+    def _reflect_box(box: Box, axis: Axis, position: float) -> Box:
+        """Reflects a box across a plane perpendicular to the given axis at the specified position.
+
+        Parameters
+        ----------
+        box : Box
+            The box to reflect
+        axis : Axis
+            The axis perpendicular to the reflection plane (0,1,2) -> (x,y,z)
+        position : float
+            Position along the axis where the reflection plane is located
+
+        Returns
+        -------
+        Box
+            The reflected box
+        """
+        new_center = list(box.center)
+        new_center[axis] = 2 * position - box.center[axis]
+        return box.updated_copy(center=new_center)
+
+    @staticmethod
+    def _apply_symmetry_to_box(box: Box, axis: Axis, position: float) -> list[Box]:
+        """Applies symmetry operation to a box along specified axis.
+
+        If the box touches the symmetry plane, merges the box with its reflection.
+        Otherwise returns both the original and reflected box.
+
+        Parameters
+        ----------
+        box : Box
+            The box to apply symmetry to
+        axis : Axis
+            The axis along which to apply symmetry (0,1,2) -> (x,y,z)
+        position : float
+            Position of the symmetry plane along the axis
+
+        Returns
+        -------
+        list[Box]
+            List containing either merged box or original and reflected boxes
+        """
+        new_box = PathSpecGenerator._reflect_box(box, axis, position)
+        if isclose(new_box.bounds[0][axis], box.bounds[1][axis]) or isclose(
+            new_box.bounds[1][axis], box.bounds[0][axis]
+        ):
+            new_size = list(box.size)
+            new_size[axis] = 2 * box.size[axis]
+            new_center = list(box.center)
+            new_center[axis] = position
+            new_box = Box(size=new_size, center=new_center)
+            return [new_box]
+        return [box, new_box]
+
+    @staticmethod
+    def _apply_symmetry(
+        symmetry: tuple[Symmetry, Symmetry, Symmetry],
+        sim_center: Coordinate,
+        normal_axis: Axis,
+        box: Box,
+    ) -> list[Box]:
+        """Applies multiple symmetry operations to a box.
+
+        Parameters
+        ----------
+        symmetry : tuple[Symmetry, Symmetry, Symmetry]
+            Symmetry conditions for each axis
+        sim_center : Coordinate
+            Center coordinates where symmetry planes intersect
+        normal_axis : Axis
+            Axis normal to the plane of interest
+        box : Box
+            The box to apply symmetries to
+
+        Returns
+        -------
+        list[Box]
+            List of boxes after applying all symmetry operations
+        """
+        symmetry = list(symmetry)
+        dims = [0, 1, 2]
+        symmetry.pop(normal_axis)
+        dims.pop(normal_axis)
+        result = [box]
+        for dim, sym in zip(dims, symmetry):
+            if sym != 0:
+                tmp_list = [
+                    PathSpecGenerator._apply_symmetry_to_box(box, dim, sim_center[dim])
+                    for box in result
+                ]
+                result = list(chain.from_iterable(tmp_list))
+        return result
diff --git a/tidy3d/components/microwave/viz.py b/tidy3d/components/microwave/viz.py
new file mode 100644
index 000000000..3f2417450
--- /dev/null
+++ b/tidy3d/components/microwave/viz.py
@@ -0,0 +1,56 @@
+"""Utilities for plotting microwave components"""
+
+from __future__ import annotations
+
+from numpy import inf
+
+from tidy3d.components.viz import PathPlotParams
+
+""" Constants """
+VOLTAGE_COLOR = "red"
+CURRENT_COLOR = "blue"
+PATH_LINEWIDTH = 2
+ARROW_CURRENT = {
+    "arrowstyle": "-|>",
+    "mutation_scale": 32,
+    "linestyle": "",
+    "lw": PATH_LINEWIDTH,
+    "color": CURRENT_COLOR,
+}
+
+plot_params_voltage_path = PathPlotParams(
+    alpha=1.0,
+    zorder=inf,
+    color=VOLTAGE_COLOR,
+    linestyle="--",
+    linewidth=PATH_LINEWIDTH,
+    marker="o",
+    markersize=10,
+    markeredgecolor=VOLTAGE_COLOR,
+    markerfacecolor="white",
+)
+
+plot_params_voltage_plus = PathPlotParams(
+    alpha=1.0,
+    zorder=inf,
+    color=VOLTAGE_COLOR,
+    marker="+",
+    markersize=6,
+)
+
+plot_params_voltage_minus = PathPlotParams(
+    alpha=1.0,
+    zorder=inf,
+    color=VOLTAGE_COLOR,
+    marker="_",
+    markersize=6,
+)
+
+plot_params_current_path = PathPlotParams(
+    alpha=1.0,
+    zorder=inf,
+    color=CURRENT_COLOR,
+    linestyle="--",
+    linewidth=PATH_LINEWIDTH,
+    marker="",
+)
diff --git a/tidy3d/components/mode/mode_solver.py b/tidy3d/components/mode/mode_solver.py
index 9f5c41725..6b338e6bd 100644
--- a/tidy3d/components/mode/mode_solver.py
+++ b/tidy3d/components/mode/mode_solver.py
@@ -33,6 +33,7 @@
     IsotropicUniformMediumType,
     LossyMetalMedium,
 )
+from tidy3d.components.microwave.path_integral_factory import make_path_integrals
 from tidy3d.components.mode_spec import ModeSpec
 from tidy3d.components.monitor import ModeMonitor, ModeSolverMonitor
 from tidy3d.components.scene import Scene
@@ -66,6 +67,7 @@
 from tidy3d.exceptions import SetupError, ValidationError
 from tidy3d.log import log
 from tidy3d.packaging import supports_local_subpixel, tidy3d_extras
+from tidy3d.plugins.microwave.impedance_calculator import ImpedanceCalculator
 
 # Importing the local solver may not work if e.g. scipy is not installed
 IMPORT_ERROR_MSG = """Could not import local solver, 'ModeSolver' objects can still be constructed
@@ -498,6 +500,9 @@ def data_raw(self) -> ModeSolverData:
         self._field_decay_warning(mode_solver_data.symmetry_expanded)
 
         mode_solver_data = self._filter_components(mode_solver_data)
+        # Calculate and add the characteristic impedance
+        if self.mode_spec.microwave_mode_spec is not None:
+            mode_solver_data = self._add_characteristic_impedance(mode_solver_data)
         return mode_solver_data
 
     @cached_property
@@ -1309,6 +1314,30 @@ def _filter_polarization(self, mode_solver_data: ModeSolverData):
             ]:
                 data.values[..., ifreq, :] = data.values[..., ifreq, sort_inds]
 
+    def _add_characteristic_impedance(self, mode_solver_data: ModeSolverData):
+        """Calculate and add characteristic impedance to ``mode_solver_data`` with the path specifications.
+        If they were not supplied by the user, then create a specification automatically.
+        """
+        path_integrals = make_path_integrals(
+            self.mode_spec.microwave_mode_spec,
+            self.to_monitor(name=MODE_MONITOR_NAME),
+            self.simulation,
+        )
+        # Need to operate on the full symmetry expanded fields
+        mode_solver_data_expanded = mode_solver_data.symmetry_expanded_copy
+        Z0_list = []
+        for mode_index in range(self.mode_spec.num_modes):
+            integral_pair = path_integrals[mode_index]
+            impedance_calc = ImpedanceCalculator(
+                voltage_integral=integral_pair[0], current_integral=integral_pair[1]
+            )
+            single_mode_data = mode_solver_data_expanded._isel(mode_index=[mode_index])
+            Z0 = impedance_calc.compute_impedance(single_mode_data)
+            Z0_list.append(Z0)
+        all_mode_Z0 = xr.concat(Z0_list, dim="mode_index")
+        all_mode_Z0 = ImpedanceCalculator._set_data_array_attributes(all_mode_Z0)
+        return mode_solver_data.updated_copy(Z0=all_mode_Z0)
+
     @cached_property
     def data(self) -> ModeSolverData:
         """:class:`.ModeSolverData` containing the field and effective index data.
diff --git a/tidy3d/components/mode_spec.py b/tidy3d/components/mode_spec.py
index 0f85bb47c..c86a71653 100644
--- a/tidy3d/components/mode_spec.py
+++ b/tidy3d/components/mode_spec.py
@@ -3,7 +3,7 @@
 from __future__ import annotations
 
 from math import isclose
-from typing import Literal, Union
+from typing import Literal, Optional, Union
 
 import numpy as np
 import pydantic.v1 as pd
@@ -13,6 +13,7 @@
 from tidy3d.log import log
 
 from .base import Tidy3dBaseModel, skip_if_fields_missing
+from .microwave.microwave_mode_spec import MicrowaveModeSpec
 from .types import Axis2D, TrackFreq
 
 GROUP_INDEX_STEP = 0.005
@@ -20,7 +21,7 @@
 
 class ModeSpec(Tidy3dBaseModel):
     """
-    Stores specifications for the mode solver to find an electromagntic mode.
+    Stores specifications for the mode solver to find an electromagnetic mode.
 
     Notes
     -----
@@ -163,6 +164,14 @@ class ModeSpec(Tidy3dBaseModel):
         f"default of {GROUP_INDEX_STEP} is used.",
     )
 
+    microwave_mode_spec: Optional[MicrowaveModeSpec] = pd.Field(
+        None,
+        title="Microwave Mode Specification",
+        description="Additional specification for microwave specific mode properties. For example, "
+        "it is used for setting up the computation for the characteristic impedance of a transmission "
+        "line mode.",
+    )
+
     @pd.validator("bend_axis", always=True)
     @skip_if_fields_missing(["bend_radius"])
     def bend_axis_given(cls, val, values):
@@ -235,3 +244,32 @@ def angle_rotation_with_phi(cls, val, values):
                 "enabled."
             )
         return val
+
+    @pd.validator("microwave_mode_spec")
+    def check_microwave_mode_spec_consistent(cls, val, values):
+        """Either the user has requested an automatic impedance calculation by leaving the path specifications empty
+        or the number of path specifications is equal to the number of modes.
+        """
+        if val is None:
+            return val
+        num_modes = values["num_modes"]
+        valid_number_voltage_specs = (
+            val.num_voltage_specs is None or val.num_voltage_specs == num_modes
+        )
+        valid_number_current_specs = (
+            val.num_current_specs is None or val.num_current_specs == num_modes
+        )
+
+        if not valid_number_voltage_specs:
+            raise SetupError(
+                f"Given {val.num_voltage_specs} voltage specifications, but the number of modes requested is {num_modes}. "
+                "Please either ensure that the number of voltage specifications is equal to the "
+                "number of modes or leave this field as 'None' in the 'MicrowaveModeSpec'."
+            )
+        if not valid_number_current_specs:
+            raise SetupError(
+                f"Given {val.num_current_specs} current specifications, but the number of modes requested is {num_modes}. "
+                "Please either ensure that the number of voltage specifications is equal to the "
+                "number of modes or leave this field as 'None' in the 'MicrowaveModeSpec'."
+            )
+        return val
diff --git a/tidy3d/components/simulation.py b/tidy3d/components/simulation.py
index 852159f30..0965414d9 100644
--- a/tidy3d/components/simulation.py
+++ b/tidy3d/components/simulation.py
@@ -3931,7 +3931,9 @@ def validate_pre_upload(self, source_required: bool = True) -> None:
         self._warn_time_monitors_outside_run_time()
         self._validate_time_monitors_num_steps()
         self._validate_freq_monitors_freq_range()
+        self._validate_microwave_mode_specs()
         _ = self.volumetric_structures
+
         log.end_capture(self)
         if source_required and len(self.sources) == 0:
             raise SetupError("No sources in simulation.")
@@ -4118,6 +4120,17 @@ def _validate_freq_monitors_freq_range(self) -> None:
                     "(Hz) as defined by the sources."
                 )
 
+    def _validate_microwave_mode_specs(self) -> None:
+        """Raise error if any microwave mode specifications fail to instantiate path integrals."""
+        from .microwave.path_integral_factory import make_path_integrals
+
+        for monitor in self.monitors:
+            if not isinstance(monitor, AbstractModeMonitor):
+                continue
+            mw_mode_spec = monitor.mode_spec.microwave_mode_spec
+            if mw_mode_spec is not None:
+                _ = make_path_integrals(mw_mode_spec, monitor, self)
+
     @cached_property
     def monitors_data_size(self) -> dict[str, float]:
         """Dictionary mapping monitor names to their estimated storage size in bytes."""
diff --git a/tidy3d/plugins/microwave/__init__.py b/tidy3d/plugins/microwave/__init__.py
index bc47c58c2..2b057188f 100644
--- a/tidy3d/plugins/microwave/__init__.py
+++ b/tidy3d/plugins/microwave/__init__.py
@@ -8,6 +8,7 @@
 )
 from .auto_path_integrals import path_integrals_from_lumped_element
 from .custom_path_integrals import (
+    CompositeCurrentIntegral,
     CustomCurrentIntegral2D,
     CustomPathIntegral2D,
     CustomVoltageIntegral2D,
@@ -23,6 +24,7 @@
 
 __all__ = [
     "AxisAlignedPathIntegral",
+    "CompositeCurrentIntegral",
     "CurrentIntegralAxisAligned",
     "CurrentIntegralTypes",
     "CustomCurrentIntegral2D",
diff --git a/tidy3d/plugins/microwave/custom_path_integrals.py b/tidy3d/plugins/microwave/custom_path_integrals.py
index d1e2c89a1..9d8256935 100644
--- a/tidy3d/plugins/microwave/custom_path_integrals.py
+++ b/tidy3d/plugins/microwave/custom_path_integrals.py
@@ -2,40 +2,36 @@
 
 from __future__ import annotations
 
-from typing import Literal, Optional
+from typing import Literal, Union
 
 import numpy as np
-import pydantic.v1 as pd
-import shapely
 import xarray as xr
 
-from tidy3d.components.base import cached_property
+from tidy3d.components.data.data_array import FreqDataArray, FreqModeDataArray
+from tidy3d.components.data.monitor_data import FieldTimeData
 from tidy3d.components.geometry.base import Geometry
-from tidy3d.components.types import ArrayFloat2D, Ax, Axis, Bound, Coordinate, Direction
-from tidy3d.components.viz import add_ax_if_none
-from tidy3d.constants import MICROMETER, fp_eps
-from tidy3d.exceptions import SetupError
+from tidy3d.components.microwave.path_spec import (
+    CompositeCurrentIntegralSpec,
+    CustomCurrentIntegral2DSpec,
+    CustomPathIntegral2DSpec,
+    CustomVoltageIntegral2DSpec,
+)
+from tidy3d.components.types import Axis, Coordinate
+from tidy3d.exceptions import DataError
+from tidy3d.log import log
 
 from .path_integrals import (
-    AbstractAxesRH,
     AxisAlignedPathIntegral,
     CurrentIntegralAxisAligned,
     IntegralResultTypes,
     MonitorDataTypes,
     VoltageIntegralAxisAligned,
 )
-from .viz import (
-    ARROW_CURRENT,
-    plot_params_current_path,
-    plot_params_voltage_minus,
-    plot_params_voltage_path,
-    plot_params_voltage_plus,
-)
 
 FieldParameter = Literal["E", "H"]
 
 
-class CustomPathIntegral2D(AbstractAxesRH):
+class CustomPathIntegral2D(CustomPathIntegral2DSpec):
     """Class for defining a custom path integral defined as a curve on an axis-aligned plane.
 
     Notes
@@ -49,27 +45,6 @@ class CustomPathIntegral2D(AbstractAxesRH):
     If the path is not closed, forward and backward differences are used at the endpoints.
     """
 
-    axis: Axis = pd.Field(
-        2, title="Axis", description="Specifies dimension of the planar axis (0,1,2) -> (x,y,z)."
-    )
-
-    position: float = pd.Field(
-        ...,
-        title="Position",
-        description="Position of the plane along the ``axis``.",
-    )
-
-    vertices: ArrayFloat2D = pd.Field(
-        ...,
-        title="Vertices",
-        description="List of (d1, d2) defining the 2 dimensional positions of the path. "
-        "The index of dimension should be in the ascending order, which means "
-        "if the axis corresponds with ``y``, the coordinates of the vertices should be (x, z). "
-        "If you wish to indicate a closed contour, the final vertex should be made "
-        "equal to the first vertex, i.e., ``vertices[-1] == vertices[0]``",
-        units=MICROMETER,
-    )
-
     def compute_integral(
         self, field: FieldParameter, em_field: MonitorDataTypes
     ) -> IntegralResultTypes:
@@ -195,43 +170,8 @@ def generate_circle_coordinates(radius: float, num_points: int, clockwise: bool)
         circle_vertices[-1, :] = circle_vertices[0, :]
         return cls(axis=normal_axis, position=normal_center, vertices=circle_vertices)
 
-    @cached_property
-    def is_closed_contour(self) -> bool:
-        """Returns ``true`` when the first vertex equals the last vertex."""
-        return np.isclose(
-            self.vertices[0, :],
-            self.vertices[-1, :],
-            rtol=fp_eps,
-            atol=np.finfo(np.float32).smallest_normal,
-        ).all()
-
-    @cached_property
-    def main_axis(self) -> Axis:
-        """Axis for performing integration."""
-        return self.axis
-
-    @pd.validator("vertices", always=True)
-    def _correct_shape(cls, val):
-        """Makes sure vertices size is correct."""
-        # overall shape of vertices
-        if val.shape[1] != 2:
-            raise SetupError(
-                "'CustomPathIntegral2D.vertices' must be a 2 dimensional array shaped (N, 2). "
-                f"Given array with shape of '{val.shape}'."
-            )
-        return val
 
-    @cached_property
-    def bounds(self) -> Bound:
-        """Helper to get the geometric bounding box of the path integral."""
-        path_min = np.amin(self.vertices, axis=0)
-        path_max = np.amax(self.vertices, axis=0)
-        min_bound = Geometry.unpop_axis(self.position, path_min, self.axis)
-        max_bound = Geometry.unpop_axis(self.position, path_max, self.axis)
-        return (min_bound, max_bound)
-
-
-class CustomVoltageIntegral2D(CustomPathIntegral2D):
+class CustomVoltageIntegral2D(CustomPathIntegral2D, CustomVoltageIntegral2DSpec):
     """Class for computing the voltage between two points defined by a custom path.
     Computed voltage is :math:`V=V_b-V_a`, where position b is the final vertex in the supplied path.
 
@@ -261,57 +201,8 @@ def compute_voltage(self, em_field: MonitorDataTypes) -> IntegralResultTypes:
         voltage = VoltageIntegralAxisAligned._set_data_array_attributes(voltage)
         return voltage
 
-    @add_ax_if_none
-    def plot(
-        self,
-        x: Optional[float] = None,
-        y: Optional[float] = None,
-        z: Optional[float] = None,
-        ax: Ax = None,
-        **path_kwargs,
-    ) -> Ax:
-        """Plot path integral at single (x,y,z) coordinate.
-
-        Parameters
-        ----------
-        x : float = None
-            Position of plane in x direction, only one of x,y,z can be specified to define plane.
-        y : float = None
-            Position of plane in y direction, only one of x,y,z can be specified to define plane.
-        z : float = None
-            Position of plane in z direction, only one of x,y,z can be specified to define plane.
-        ax : matplotlib.axes._subplots.Axes = None
-            Matplotlib axes to plot on, if not specified, one is created.
-        **path_kwargs
-            Optional keyword arguments passed to the matplotlib plotting of the line.
-            For details on accepted values, refer to
-            `Matplotlib's documentation <https://tinyurl.com/36marrat>`_.
-
-        Returns
-        -------
-        matplotlib.axes._subplots.Axes
-            The supplied or created matplotlib axes.
-        """
-        axis, position = Geometry.parse_xyz_kwargs(x=x, y=y, z=z)
-        if axis != self.main_axis or not np.isclose(position, self.position, rtol=fp_eps):
-            return ax
-
-        plot_params = plot_params_voltage_path.include_kwargs(**path_kwargs)
-        plot_kwargs = plot_params.to_kwargs()
-        xs = self.vertices[:, 0]
-        ys = self.vertices[:, 1]
-        ax.plot(xs, ys, markevery=[0, -1], **plot_kwargs)
-
-        # Plot special end points
-        end_kwargs = plot_params_voltage_plus.include_kwargs(**path_kwargs).to_kwargs()
-        start_kwargs = plot_params_voltage_minus.include_kwargs(**path_kwargs).to_kwargs()
-        ax.plot(xs[0], ys[0], **start_kwargs)
-        ax.plot(xs[-1], ys[-1], **end_kwargs)
-
-        return ax
-
 
-class CustomCurrentIntegral2D(CustomPathIntegral2D):
+class CustomCurrentIntegral2D(CustomPathIntegral2D, CustomCurrentIntegral2DSpec):
     """Class for computing conduction current via Ampère's circuital law on a custom path.
     To compute the current flowing in the positive ``axis`` direction, the vertices should be
     ordered in a counterclockwise direction."""
@@ -334,64 +225,150 @@ def compute_current(self, em_field: MonitorDataTypes) -> IntegralResultTypes:
         current = CurrentIntegralAxisAligned._set_data_array_attributes(current)
         return current
 
-    @add_ax_if_none
-    def plot(
+
+class CompositeCurrentIntegral(CompositeCurrentIntegralSpec):
+    """Current integral comprising one or more disjoint paths"""
+
+    def compute_current(self, em_field: MonitorDataTypes) -> IntegralResultTypes:
+        """Compute current flowing in loop defined by the outer edge of a rectangle."""
+        if isinstance(em_field, FieldTimeData) and self.sum_spec == "split":
+            raise DataError(
+                "Only frequency domain field data is supported when using the 'split' sum_spec. "
+                "Either switch the sum_spec to 'sum' or supply frequency domain data."
+            )
+
+        from tidy3d.components.microwave.path_integral_factory import make_current_integral
+
+        current_integrals = [make_current_integral(path_spec) for path_spec in self.path_specs]
+
+        # Initialize arrays with first current term
+        first_term = current_integrals[0].compute_current(em_field)
+        current_in_phase = xr.zeros_like(first_term)
+        current_out_phase = xr.zeros_like(first_term)
+
+        # Get reference phase from first non-zero current
+        phase_reference = None
+        for path in current_integrals:
+            term = path.compute_current(em_field)
+            if np.any(term.abs > 0):
+                phase_reference = term.angle
+                break
+        if phase_reference is None:
+            raise DataError("Cannot complete calculation of current. No non-zero current found.")
+
+        # Accumulate currents based on phase comparison
+        for path in current_integrals:
+            term = path.compute_current(em_field)
+            if np.all(term.abs == 0):
+                continue
+
+            # Compare phase to reference
+            phase_diff = term.angle - phase_reference
+            # Wrap phase difference to [-pi, pi]
+            phase_diff.values = np.mod(phase_diff.values + np.pi, 2 * np.pi) - np.pi
+
+            # Check phase consistency across frequencies
+            self._check_phase_sign_consistency(phase_diff)
+
+            # Add to in-phase or out-of-phase current
+            is_in_phase = phase_diff <= np.pi / 2
+            current_in_phase += xr.where(is_in_phase, term, 0)
+            current_out_phase += xr.where(~is_in_phase, term, 0)
+
+        current_in_phase = CurrentIntegralAxisAligned._set_data_array_attributes(current_in_phase)
+        current_out_phase = CurrentIntegralAxisAligned._set_data_array_attributes(current_out_phase)
+
+        if self.sum_spec == "sum":
+            return current_in_phase + current_out_phase
+
+        # Check amplitude consistency across frequencies
+        self._check_phase_amplitude_consistency(current_in_phase, current_out_phase)
+
+        # For split mode, return the larger magnitude current
+        current = xr.where(
+            abs(current_in_phase) >= abs(current_out_phase), current_in_phase, current_out_phase
+        )
+        return CurrentIntegralAxisAligned._set_data_array_attributes(current)
+
+    def _check_phase_sign_consistency(
         self,
-        x: Optional[float] = None,
-        y: Optional[float] = None,
-        z: Optional[float] = None,
-        ax: Ax = None,
-        **path_kwargs,
-    ) -> Ax:
-        """Plot path integral at single (x,y,z) coordinate.
+        phase_difference: Union[FreqDataArray, FreqModeDataArray],
+    ) -> bool:
+        """
+        Check that the provided current data has a consistent phase with respect to the reference
+        phase. A consistent phase allows for the automatic identification of currents flowing in
+        opposite directions. However, when the provided data does not correspond with a transmission
+        line mode, this consistent phase condition will likely fail, so we emit a warning here to
+        notify the user.
+        """
 
-        Parameters
-        ----------
-        x : float = None
-            Position of plane in x direction, only one of x,y,z can be specified to define plane.
-        y : float = None
-            Position of plane in y direction, only one of x,y,z can be specified to define plane.
-        z : float = None
-            Position of plane in z direction, only one of x,y,z can be specified to define plane.
-        ax : matplotlib.axes._subplots.Axes = None
-            Matplotlib axes to plot on, if not specified, one is created.
-        **path_kwargs
-            Optional keyword arguments passed to the matplotlib plotting of the line.
-            For details on accepted values, refer to
-            `Matplotlib's documentation <https://tinyurl.com/36marrat>`_.
+        # Check phase consistency across frequencies
+        freq_axis = phase_difference.get_axis_num("f")
+        all_in_phase = np.all(abs(phase_difference) <= np.pi / 2, axis=freq_axis)
+        all_out_of_phase = np.all(abs(phase_difference) > np.pi / 2, axis=freq_axis)
+        consistent_phase = np.logical_or(all_in_phase, all_out_of_phase)
+
+        if not np.all(consistent_phase) and self.sum_spec == "split":
+            warning_msg = (
+                "Phase alignment of computed current is not consistent across frequencies. "
+                "The provided fields are not suitable for the 'split' method of computing current. "
+                "Please provide the current path specifications manually."
+            )
 
-        Returns
-        -------
-        matplotlib.axes._subplots.Axes
-            The supplied or created matplotlib axes.
+            if isinstance(phase_difference, FreqModeDataArray):
+                inconsistent_modes = []
+                mode_indices = phase_difference.mode_index.values
+                for mode_idx in range(len(mode_indices)):
+                    if not consistent_phase[mode_idx]:
+                        inconsistent_modes.append(mode_idx)
+
+                warning_msg += (
+                    f" Modes with indices {inconsistent_modes} violated the phase consistency "
+                    "requirement."
+                )
+
+            log.warning(warning_msg)
+
+            return False
+        return True
+
+    def _check_phase_amplitude_consistency(
+        self,
+        current_in_phase: Union[FreqDataArray, FreqModeDataArray],
+        current_out_phase: Union[FreqDataArray, FreqModeDataArray],
+    ) -> bool:
         """
-        axis, position = Geometry.parse_xyz_kwargs(x=x, y=y, z=z)
-        if axis != self.main_axis or not np.isclose(position, self.position, rtol=fp_eps):
-            return ax
-
-        plot_params = plot_params_current_path.include_kwargs(**path_kwargs)
-        plot_kwargs = plot_params.to_kwargs()
-        xs = self.vertices[:, 0]
-        ys = self.vertices[:, 1]
-        ax.plot(xs, ys, **plot_kwargs)
-
-        # Add arrow at start of contour
-        ax.annotate(
-            "",
-            xytext=(xs[0], ys[0]),
-            xy=(xs[1], ys[1]),
-            arrowprops=ARROW_CURRENT,
-        )
-        return ax
-
-    @cached_property
-    def sign(self) -> Direction:
-        """Uses the ordering of the vertices to determine the direction of the current flow."""
-        linestr = shapely.LineString(coordinates=self.vertices)
-        is_ccw = shapely.is_ccw(linestr)
-        # Invert statement when the vertices are given as (x, z)
-        if self.axis == 1:
-            is_ccw = not is_ccw
-        if is_ccw:
-            return "+"
-        return "-"
+        Check that the summed in phase and out of phase components of current have a consistent relative amplitude.
+        A consistent amplitude across frequencies allows for the automatic identification of the total conduction
+        current flowing in the transmission line. If the amplitudes are not consistent, we emit a warning.
+        """
+
+        # For split mode, return the larger magnitude current
+        freq_axis = current_in_phase.get_axis_num("f")
+        in_all_larger = np.all(abs(current_in_phase) >= abs(current_out_phase), axis=freq_axis)
+        in_all_smaller = np.all(abs(current_in_phase) < abs(current_out_phase), axis=freq_axis)
+        consistent_max_current = np.logical_or(in_all_larger, in_all_smaller)
+        if not np.all(consistent_max_current) and self.sum_spec == "split":
+            warning_msg = (
+                "There is not a consistently larger current across frequencies between the in-phase "
+                "and out-of-phase components. The provided fields are not suitable for the "
+                "'split' method of computing current. Please provide the current path "
+                "specifications manually."
+            )
+
+            if isinstance(current_in_phase, FreqModeDataArray):
+                inconsistent_modes = []
+                mode_indices = current_in_phase.mode_index.values
+                for mode_idx in range(len(mode_indices)):
+                    if not consistent_max_current[mode_idx]:
+                        inconsistent_modes.append(mode_idx)
+
+                warning_msg += (
+                    f" Modes with indices {inconsistent_modes} violated the amplitude consistency "
+                    "requirement."
+                )
+
+            log.warning(warning_msg)
+
+            return False
+        return True
diff --git a/tidy3d/plugins/microwave/impedance_calculator.py b/tidy3d/plugins/microwave/impedance_calculator.py
index a40081862..7800625b5 100644
--- a/tidy3d/plugins/microwave/impedance_calculator.py
+++ b/tidy3d/plugins/microwave/impedance_calculator.py
@@ -15,7 +15,11 @@
 from tidy3d.exceptions import ValidationError
 from tidy3d.log import log
 
-from .custom_path_integrals import CustomCurrentIntegral2D, CustomVoltageIntegral2D
+from .custom_path_integrals import (
+    CompositeCurrentIntegral,
+    CustomCurrentIntegral2D,
+    CustomVoltageIntegral2D,
+)
 from .path_integrals import (
     AxisAlignedPathIntegral,
     CurrentIntegralAxisAligned,
@@ -25,7 +29,9 @@
 )
 
 VoltageIntegralTypes = Union[VoltageIntegralAxisAligned, CustomVoltageIntegral2D]
-CurrentIntegralTypes = Union[CurrentIntegralAxisAligned, CustomCurrentIntegral2D]
+CurrentIntegralTypes = Union[
+    CurrentIntegralAxisAligned, CustomCurrentIntegral2D, CompositeCurrentIntegral
+]
 
 
 class ImpedanceCalculator(Tidy3dBaseModel):
@@ -98,6 +104,7 @@ def compute_impedance(self, em_field: MonitorDataTypes) -> IntegralResultTypes:
                 impedance = np.real(voltage) / np.real(current)
             else:
                 impedance = voltage / current
+                impedance = voltage / current
         impedance = ImpedanceCalculator._set_data_array_attributes(impedance)
         return impedance
 
@@ -115,6 +122,13 @@ def check_voltage_or_current(cls, val, values):
     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)
+        # Determine type based on coords present
         if "mode_index" in data_array.coords:
             data_array = FreqModeDataArray(data_array)
         elif "f" in data_array.coords:
diff --git a/tidy3d/plugins/microwave/path_integrals.py b/tidy3d/plugins/microwave/path_integrals.py
index 802f5b0d9..2bd9e4f66 100644
--- a/tidy3d/plugins/microwave/path_integrals.py
+++ b/tidy3d/plugins/microwave/path_integrals.py
@@ -2,14 +2,11 @@
 
 from __future__ import annotations
 
-from abc import ABC, abstractmethod
 from typing import Optional, Union
 
 import numpy as np
-import pydantic.v1 as pd
 import xarray as xr
 
-from tidy3d.components.base import Tidy3dBaseModel, cached_property
 from tidy3d.components.data.data_array import (
     FreqDataArray,
     FreqModeDataArray,
@@ -19,88 +16,23 @@
     TimeDataArray,
 )
 from tidy3d.components.data.monitor_data import FieldData, FieldTimeData, ModeData, ModeSolverData
-from tidy3d.components.geometry.base import Box, Geometry
-from tidy3d.components.types import Ax, Axis, Coordinate2D, Direction
-from tidy3d.components.validators import assert_line, assert_plane
-from tidy3d.components.viz import add_ax_if_none
-from tidy3d.constants import AMP, VOLT, fp_eps
-from tidy3d.exceptions import DataError, Tidy3dError
-from tidy3d.log import log
-
-from .viz import (
-    ARROW_CURRENT,
-    plot_params_current_path,
-    plot_params_voltage_minus,
-    plot_params_voltage_path,
-    plot_params_voltage_plus,
+from tidy3d.components.geometry.base import Geometry
+from tidy3d.components.microwave.path_spec import (
+    AxisAlignedPathIntegralSpec,
+    CurrentIntegralAxisAlignedSpec,
+    VoltageIntegralAxisAlignedSpec,
 )
+from tidy3d.constants import AMP, VOLT, fp_eps
+from tidy3d.exceptions import DataError
 
 MonitorDataTypes = Union[FieldData, FieldTimeData, ModeData, ModeSolverData]
 EMScalarFieldType = Union[ScalarFieldDataArray, ScalarFieldTimeDataArray, ScalarModeFieldDataArray]
 IntegralResultTypes = Union[FreqDataArray, FreqModeDataArray, TimeDataArray]
 
 
-class AbstractAxesRH(Tidy3dBaseModel, ABC):
-    """Represents an axis-aligned right-handed coordinate system with one axis preferred.
-    Typically `main_axis` would refer to the normal axis of a plane.
-    """
-
-    @cached_property
-    @abstractmethod
-    def main_axis(self) -> Axis:
-        """Get the preferred axis."""
-
-    @cached_property
-    def remaining_axes(self) -> tuple[Axis, Axis]:
-        """Get in-plane axes, ordered to maintain a right-handed coordinate system."""
-        axes: list[Axis] = [0, 1, 2]
-        axes.pop(self.main_axis)
-        if self.main_axis == 1:
-            return (axes[1], axes[0])
-        return (axes[0], axes[1])
-
-    @cached_property
-    def remaining_dims(self) -> tuple[str, str]:
-        """Get in-plane dimensions, ordered to maintain a right-handed coordinate system."""
-        dim1 = "xyz"[self.remaining_axes[0]]
-        dim2 = "xyz"[self.remaining_axes[1]]
-        return (dim1, dim2)
-
-    @cached_property
-    def local_dims(self) -> tuple[str, str, str]:
-        """Get in-plane dimensions with in-plane dims first, followed by the `main_axis` dimension."""
-        dim3 = "xyz"[self.main_axis]
-        return self.remaining_dims + tuple(dim3)
-
-    @pd.root_validator(pre=False)
-    def _warn_rf_license(cls, values):
-        log.warning(
-            "ℹ️ ⚠️ RF simulations are subject to new license requirements in the future. You have instantiated at least one RF-specific component.",
-            log_once=True,
-        )
-        return values
-
-
-class AxisAlignedPathIntegral(AbstractAxesRH, Box):
+class AxisAlignedPathIntegral(AxisAlignedPathIntegralSpec):
     """Class for defining the simplest type of path integral, which is aligned with Cartesian axes."""
 
-    _line_validator = assert_line()
-
-    extrapolate_to_endpoints: bool = pd.Field(
-        False,
-        title="Extrapolate to Endpoints",
-        description="If the endpoints of the path integral terminate at or near a material interface, "
-        "the field is likely discontinuous. When this field is ``True``, fields that are outside and on the bounds "
-        "of the integral are ignored. Should be enabled when computing voltage between two conductors.",
-    )
-
-    snap_path_to_grid: bool = pd.Field(
-        False,
-        title="Snap Path to Grid",
-        description="It might be desireable to integrate exactly along the Yee grid associated with "
-        "a field. When this field is ``True``, the integration path will be snapped to the grid.",
-    )
-
     def compute_integral(self, scalar_field: EMScalarFieldType) -> IntegralResultTypes:
         """Computes the defined integral given the input ``scalar_field``."""
 
@@ -176,25 +108,6 @@ def _get_field_along_path(self, scalar_field: EMScalarFieldType) -> EMScalarFiel
         scalar_field = scalar_field.reset_coords(drop=True)
         return scalar_field
 
-    @cached_property
-    def main_axis(self) -> Axis:
-        """Axis for performing integration."""
-        for index, value in enumerate(self.size):
-            if value != 0:
-                return index
-        raise Tidy3dError("Failed to identify axis.")
-
-    def _vertices_2D(self, axis: Axis) -> tuple[Coordinate2D, Coordinate2D]:
-        """Returns the two vertices of this path in the plane defined by ``axis``."""
-        min = self.bounds[0]
-        max = self.bounds[1]
-        _, min = Box.pop_axis(min, axis)
-        _, max = Box.pop_axis(max, axis)
-
-        u = [min[0], max[0]]
-        v = [min[1], max[1]]
-        return (u, v)
-
     @staticmethod
     def _check_monitor_data_supported(em_field: MonitorDataTypes):
         """Helper for validating that monitor data is supported."""
@@ -215,15 +128,9 @@ def _make_result_data_array(result: xr.DataArray) -> IntegralResultTypes:
         return FreqDataArray(data=result.data, coords=result.coords)
 
 
-class VoltageIntegralAxisAligned(AxisAlignedPathIntegral):
+class VoltageIntegralAxisAligned(AxisAlignedPathIntegral, VoltageIntegralAxisAlignedSpec):
     """Class for computing the voltage between two points defined by an axis-aligned line."""
 
-    sign: Direction = pd.Field(
-        ...,
-        title="Direction of Path Integral",
-        description="Positive indicates V=Vb-Va where position b has a larger coordinate along the axis of integration.",
-    )
-
     def compute_voltage(self, em_field: MonitorDataTypes) -> IntegralResultTypes:
         """Compute voltage along path defined by a line."""
         self._check_monitor_data_supported(em_field=em_field)
@@ -245,6 +152,13 @@ def compute_voltage(self, em_field: MonitorDataTypes) -> IntegralResultTypes:
     @staticmethod
     def _set_data_array_attributes(data_array: IntegralResultTypes) -> IntegralResultTypes:
         """Add explanatory attributes to the data array."""
+        # 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 = "V"
         return data_array.assign_attrs(units=VOLT, long_name="voltage")
 
@@ -305,82 +219,10 @@ def from_terminal_positions(
             sign=direction,
         )
 
-    @add_ax_if_none
-    def plot(
-        self,
-        x: Optional[float] = None,
-        y: Optional[float] = None,
-        z: Optional[float] = None,
-        ax: Ax = None,
-        **path_kwargs,
-    ) -> Ax:
-        """Plot path integral at single (x,y,z) coordinate.
 
-        Parameters
-        ----------
-        x : float = None
-            Position of plane in x direction, only one of x,y,z can be specified to define plane.
-        y : float = None
-            Position of plane in y direction, only one of x,y,z can be specified to define plane.
-        z : float = None
-            Position of plane in z direction, only one of x,y,z can be specified to define plane.
-        ax : matplotlib.axes._subplots.Axes = None
-            Matplotlib axes to plot on, if not specified, one is created.
-        **path_kwargs
-            Optional keyword arguments passed to the matplotlib plotting of the line.
-            For details on accepted values, refer to
-            `Matplotlib's documentation <https://tinyurl.com/36marrat>`_.
-
-        Returns
-        -------
-        matplotlib.axes._subplots.Axes
-            The supplied or created matplotlib axes.
-        """
-        axis, position = self.parse_xyz_kwargs(x=x, y=y, z=z)
-        if axis == self.main_axis or not np.isclose(position, self.center[axis], rtol=fp_eps):
-            return ax
-
-        (xs, ys) = self._vertices_2D(axis)
-        # Plot the path
-        plot_params = plot_params_voltage_path.include_kwargs(**path_kwargs)
-        plot_kwargs = plot_params.to_kwargs()
-        ax.plot(xs, ys, markevery=[0, -1], **plot_kwargs)
-
-        # Plot special end points
-        end_kwargs = plot_params_voltage_plus.include_kwargs(**path_kwargs).to_kwargs()
-        start_kwargs = plot_params_voltage_minus.include_kwargs(**path_kwargs).to_kwargs()
-
-        if self.sign == "-":
-            start_kwargs, end_kwargs = end_kwargs, start_kwargs
-
-        ax.plot(xs[0], ys[0], **start_kwargs)
-        ax.plot(xs[1], ys[1], **end_kwargs)
-        return ax
-
-
-class CurrentIntegralAxisAligned(AbstractAxesRH, Box):
+class CurrentIntegralAxisAligned(CurrentIntegralAxisAlignedSpec):
     """Class for computing conduction current via Ampère's circuital law on an axis-aligned loop."""
 
-    _plane_validator = assert_plane()
-
-    sign: Direction = pd.Field(
-        ...,
-        title="Direction of Contour Integral",
-        description="Positive indicates current flowing in the positive normal axis direction.",
-    )
-
-    extrapolate_to_endpoints: bool = pd.Field(
-        False,
-        title="Extrapolate to Endpoints",
-        description="This parameter is passed to :class:`AxisAlignedPathIntegral` objects when computing the contour integral.",
-    )
-
-    snap_contour_to_grid: bool = pd.Field(
-        False,
-        title="Snap Contour to Grid",
-        description="This parameter is passed to :class:`AxisAlignedPathIntegral` objects when computing the contour integral.",
-    )
-
     def compute_current(self, em_field: MonitorDataTypes) -> IntegralResultTypes:
         """Compute current flowing in loop defined by the outer edge of a rectangle."""
         AxisAlignedPathIntegral._check_monitor_data_supported(em_field=em_field)
@@ -410,174 +252,26 @@ def compute_current(self, em_field: MonitorDataTypes) -> IntegralResultTypes:
         current = CurrentIntegralAxisAligned._set_data_array_attributes(current)
         return current
 
-    @cached_property
-    def main_axis(self) -> Axis:
-        """Axis normal to loop"""
-        for index, value in enumerate(self.size):
-            if value == 0:
-                return index
-        raise Tidy3dError("Failed to identify axis.")
-
     def _to_path_integrals(
         self, h_horizontal=None, h_vertical=None
     ) -> tuple[AxisAlignedPathIntegral, ...]:
         """Returns four ``AxisAlignedPathIntegral`` instances, which represent a contour
         integral around the surface defined by ``self.size``."""
-        ax1 = self.remaining_axes[0]
-        ax2 = self.remaining_axes[1]
-
-        horizontal_passed = h_horizontal is not None
-        vertical_passed = h_vertical is not None
-        if self.snap_contour_to_grid and horizontal_passed and vertical_passed:
-            (coord1, coord2) = self.remaining_dims
-
-            # Locations where horizontal paths will be snapped
-            v_bounds = [
-                self.center[ax2] - self.size[ax2] / 2,
-                self.center[ax2] + self.size[ax2] / 2,
-            ]
-            h_snaps = h_horizontal.sel({coord2: v_bounds}, method="nearest").coords[coord2].values
-            # Locations where vertical paths will be snapped
-            h_bounds = [
-                self.center[ax1] - self.size[ax1] / 2,
-                self.center[ax1] + self.size[ax1] / 2,
-            ]
-            v_snaps = h_vertical.sel({coord1: h_bounds}, method="nearest").coords[coord1].values
-
-            bottom_bound = h_snaps[0]
-            top_bound = h_snaps[1]
-            left_bound = v_snaps[0]
-            right_bound = v_snaps[1]
-        else:
-            bottom_bound = self.bounds[0][ax2]
-            top_bound = self.bounds[1][ax2]
-            left_bound = self.bounds[0][ax1]
-            right_bound = self.bounds[1][ax1]
-
-        # Horizontal paths
-        path_size = list(self.size)
-        path_size[ax1] = right_bound - left_bound
-        path_size[ax2] = 0
-        path_center = list(self.center)
-        path_center[ax2] = bottom_bound
-
-        bottom = AxisAlignedPathIntegral(
-            center=path_center,
-            size=path_size,
-            extrapolate_to_endpoints=self.extrapolate_to_endpoints,
-            snap_path_to_grid=self.snap_contour_to_grid,
-        )
-        path_center[ax2] = top_bound
-        top = AxisAlignedPathIntegral(
-            center=path_center,
-            size=path_size,
-            extrapolate_to_endpoints=self.extrapolate_to_endpoints,
-            snap_path_to_grid=self.snap_contour_to_grid,
-        )
-
-        # Vertical paths
-        path_size = list(self.size)
-        path_size[ax1] = 0
-        path_size[ax2] = top_bound - bottom_bound
-        path_center = list(self.center)
-
-        path_center[ax1] = left_bound
-        left = AxisAlignedPathIntegral(
-            center=path_center,
-            size=path_size,
-            extrapolate_to_endpoints=self.extrapolate_to_endpoints,
-            snap_path_to_grid=self.snap_contour_to_grid,
-        )
-        path_center[ax1] = right_bound
-        right = AxisAlignedPathIntegral(
-            center=path_center,
-            size=path_size,
-            extrapolate_to_endpoints=self.extrapolate_to_endpoints,
-            snap_path_to_grid=self.snap_contour_to_grid,
+        path_specs = self._to_path_integral_specs(h_horizontal=h_horizontal, h_vertical=h_vertical)
+        path_integrals = tuple(
+            AxisAlignedPathIntegral(**path_spec.dict(exclude={"type"})) for path_spec in path_specs
         )
-
-        return (bottom, right, top, left)
+        return path_integrals
 
     @staticmethod
     def _set_data_array_attributes(data_array: IntegralResultTypes) -> IntegralResultTypes:
         """Add explanatory attributes to the data array."""
+        # 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 = "I"
         return data_array.assign_attrs(units=AMP, long_name="current")
-
-    @add_ax_if_none
-    def plot(
-        self,
-        x: Optional[float] = None,
-        y: Optional[float] = None,
-        z: Optional[float] = None,
-        ax: Ax = None,
-        **path_kwargs,
-    ) -> Ax:
-        """Plot path integral at single (x,y,z) coordinate.
-
-        Parameters
-        ----------
-        x : float = None
-            Position of plane in x direction, only one of x,y,z can be specified to define plane.
-        y : float = None
-            Position of plane in y direction, only one of x,y,z can be specified to define plane.
-        z : float = None
-            Position of plane in z direction, only one of x,y,z can be specified to define plane.
-        ax : matplotlib.axes._subplots.Axes = None
-            Matplotlib axes to plot on, if not specified, one is created.
-        **path_kwargs
-            Optional keyword arguments passed to the matplotlib plotting of the line.
-            For details on accepted values, refer to
-            `Matplotlib's documentation <https://tinyurl.com/36marrat>`_.
-
-        Returns
-        -------
-        matplotlib.axes._subplots.Axes
-            The supplied or created matplotlib axes.
-        """
-        axis, position = self.parse_xyz_kwargs(x=x, y=y, z=z)
-        if axis != self.main_axis or not np.isclose(position, self.center[axis], rtol=fp_eps):
-            return ax
-
-        plot_params = plot_params_current_path.include_kwargs(**path_kwargs)
-        plot_kwargs = plot_params.to_kwargs()
-        path_integrals = self._to_path_integrals()
-        # Plot the path
-        for path in path_integrals:
-            (xs, ys) = path._vertices_2D(axis)
-            ax.plot(xs, ys, **plot_kwargs)
-
-        (ax1, ax2) = self.remaining_axes
-
-        # Add arrow to bottom path, unless right path is longer
-        arrow_path = path_integrals[0]
-        if self.size[ax2] > self.size[ax1]:
-            arrow_path = path_integrals[1]
-
-        (xs, ys) = arrow_path._vertices_2D(axis)
-        X = (xs[0] + xs[1]) / 2
-        Y = (ys[0] + ys[1]) / 2
-        center = np.array([X, Y])
-        dx = xs[1] - xs[0]
-        dy = ys[1] - ys[0]
-        direction = np.array([dx, dy])
-        segment_length = np.linalg.norm(direction)
-        unit_dir = direction / segment_length
-
-        # Change direction of arrow depending on sign of current definition
-        if self.sign == "-":
-            unit_dir *= -1.0
-        # Change direction of arrow when the "y" axis is dropped,
-        # since the plotted coordinate system will be left-handed (x, z)
-        if self.main_axis == 1:
-            unit_dir *= -1.0
-
-        start = center - unit_dir * segment_length
-        end = center
-        ax.annotate(
-            "",
-            xytext=(start[0], start[1]),
-            xy=(end[0], end[1]),
-            arrowprops=ARROW_CURRENT,
-        )
-        return ax
diff --git a/tidy3d/plugins/microwave/viz.py b/tidy3d/plugins/microwave/viz.py
index 72f4570c8..8bdc86b3e 100644
--- a/tidy3d/plugins/microwave/viz.py
+++ b/tidy3d/plugins/microwave/viz.py
@@ -7,56 +7,9 @@
 from tidy3d.components.viz import PathPlotParams
 
 """ Constants """
-VOLTAGE_COLOR = "red"
-CURRENT_COLOR = "blue"
 LOBE_PEAK_COLOR = "tab:red"
 LOBE_WIDTH_COLOR = "tab:orange"
 LOBE_FNBW_COLOR = "tab:blue"
-PATH_LINEWIDTH = 2
-ARROW_CURRENT = {
-    "arrowstyle": "-|>",
-    "mutation_scale": 32,
-    "linestyle": "",
-    "lw": PATH_LINEWIDTH,
-    "color": CURRENT_COLOR,
-}
-
-plot_params_voltage_path = PathPlotParams(
-    alpha=1.0,
-    zorder=inf,
-    color=VOLTAGE_COLOR,
-    linestyle="--",
-    linewidth=PATH_LINEWIDTH,
-    marker="o",
-    markersize=10,
-    markeredgecolor=VOLTAGE_COLOR,
-    markerfacecolor="white",
-)
-
-plot_params_voltage_plus = PathPlotParams(
-    alpha=1.0,
-    zorder=inf,
-    color=VOLTAGE_COLOR,
-    marker="+",
-    markersize=6,
-)
-
-plot_params_voltage_minus = PathPlotParams(
-    alpha=1.0,
-    zorder=inf,
-    color=VOLTAGE_COLOR,
-    marker="_",
-    markersize=6,
-)
-
-plot_params_current_path = PathPlotParams(
-    alpha=1.0,
-    zorder=inf,
-    color=CURRENT_COLOR,
-    linestyle="--",
-    linewidth=PATH_LINEWIDTH,
-    marker="",
-)
 
 plot_params_lobe_peak = PathPlotParams(
     alpha=1.0,
diff --git a/tidy3d/plugins/smatrix/ports/coaxial_lumped.py b/tidy3d/plugins/smatrix/ports/coaxial_lumped.py
index 44655bffc..00820b5fc 100644
--- a/tidy3d/plugins/smatrix/ports/coaxial_lumped.py
+++ b/tidy3d/plugins/smatrix/ports/coaxial_lumped.py
@@ -15,6 +15,7 @@
 from tidy3d.components.geometry.utils_2d import increment_float
 from tidy3d.components.grid.grid import Grid, YeeGrid
 from tidy3d.components.lumped_element import CoaxialLumpedResistor
+from tidy3d.components.microwave.path_spec import AbstractAxesRH
 from tidy3d.components.monitor import FieldMonitor
 from tidy3d.components.source.current import CustomCurrentSource
 from tidy3d.components.source.time import GaussianPulse
@@ -23,7 +24,6 @@
 from tidy3d.constants import MICROMETER
 from tidy3d.exceptions import SetupError, ValidationError
 from tidy3d.plugins.microwave import CustomCurrentIntegral2D, VoltageIntegralAxisAligned
-from tidy3d.plugins.microwave.path_integrals import AbstractAxesRH
 
 from .base_lumped import AbstractLumpedPort