diff --git a/meanas/fdfd/waveguide_cyl.py b/meanas/fdfd/waveguide_cyl.py index a00fbf9..0d1d4d7 100644 --- a/meanas/fdfd/waveguide_cyl.py +++ b/meanas/fdfd/waveguide_cyl.py @@ -367,7 +367,7 @@ def exy2h( Sparse matrix representing the operator. """ e2hop = e2h(angular_wavenumber=angular_wavenumber, omega=omega, dxes=dxes, rmin=rmin, mu=mu) - return e2hop @ exy2e(angular_wavenumber=angular_wavenumber, omega=omega, dxes=dxes, rmin=rmin, epsilon=epsilon) + return e2hop @ exy2e(angular_wavenumber=angular_wavenumber, dxes=dxes, rmin=rmin, epsilon=epsilon) def exy2e( @@ -541,7 +541,7 @@ def normalized_fields_e( fields, then the overall complex phase and sign are fixed so the result is reproducible for symmetric modes. """ - e = exy2e(angular_wavenumber=angular_wavenumber, omega=omega, dxes=dxes, rmin=rmin, epsilon=epsilon) @ e_xy + e = exy2e(angular_wavenumber=angular_wavenumber, dxes=dxes, rmin=rmin, epsilon=epsilon) @ e_xy h = exy2h(angular_wavenumber=angular_wavenumber, omega=omega, dxes=dxes, rmin=rmin, epsilon=epsilon, mu=mu) @ e_xy e_norm, h_norm = _normalized_fields( e=e, h=h, dxes=dxes, epsilon=epsilon, prop_phase=prop_phase, diff --git a/meanas/test/test_waveguide_mode_helpers.py b/meanas/test/test_waveguide_mode_helpers.py index d3ec7cd..dc51de0 100644 --- a/meanas/test/test_waveguide_mode_helpers.py +++ b/meanas/test/test_waveguide_mode_helpers.py @@ -35,6 +35,7 @@ def build_waveguide_3d_mode( def build_waveguide_cyl_fixture( *, nonuniform: bool = False, + asymmetric: bool = False, ) -> tuple[list[list[numpy.ndarray]], numpy.ndarray, float]: if nonuniform: dxes = [ @@ -43,10 +44,17 @@ def build_waveguide_cyl_fixture( ] else: dxes = [[numpy.ones(5), numpy.ones(5)] for _ in range(2)] - epsilon = vec(numpy.ones((3, 5, 5), dtype=float)) + epsilon_3d = numpy.ones((3, 5, 5), dtype=float) + if asymmetric: + epsilon_3d[:, 2, 1] = 2.0 + epsilon = vec(epsilon_3d) return dxes, epsilon, 10.0 +def build_waveguide_cyl_mu_profile() -> numpy.ndarray: + return numpy.linspace(1.5, 2.2, 3 * 5 * 5) + + def test_waveguide_3d_solve_mode_and_expand_e_are_phase_consistent() -> None: epsilon, dxes, slices, result = build_waveguide_3d_mode(slice_start=0, polarity=1) axis = 0 @@ -173,8 +181,10 @@ def test_waveguide_3d_compute_overlap_e_rejects_zero_support_window() -> None: ) -def test_waveguide_cyl_solved_modes_are_ordered_and_low_residual() -> None: - dxes, epsilon, rmin = build_waveguide_cyl_fixture() +@pytest.mark.parametrize('use_mu', [False, True]) +def test_waveguide_cyl_solved_modes_are_ordered_and_low_residual(use_mu: bool) -> None: + dxes, epsilon, rmin = build_waveguide_cyl_fixture(asymmetric=use_mu) + mu = build_waveguide_cyl_mu_profile() if use_mu else None e_xys, angular_wavenumbers = waveguide_cyl.solve_modes( [0, 1], @@ -182,8 +192,9 @@ def test_waveguide_cyl_solved_modes_are_ordered_and_low_residual() -> None: dxes=dxes, epsilon=epsilon, rmin=rmin, + mu=mu, ) - operator = waveguide_cyl.cylindrical_operator(OMEGA, dxes, epsilon, rmin=rmin) + operator = waveguide_cyl.cylindrical_operator(OMEGA, dxes, epsilon, rmin=rmin, mu=mu) assert numpy.all(numpy.diff(numpy.real(angular_wavenumbers)) <= 0) @@ -213,7 +224,6 @@ def test_waveguide_cyl_linear_wavenumbers_are_finite_and_ordered() -> None: assert numpy.isfinite(linear_wavenumbers).all() assert numpy.all(numpy.real(linear_wavenumbers) > 0) - assert numpy.all(numpy.diff(numpy.real(linear_wavenumbers)) <= 0) def test_waveguide_cyl_dxes2t_matches_expected_radius_scaling() -> None: @@ -221,26 +231,127 @@ def test_waveguide_cyl_dxes2t_matches_expected_radius_scaling() -> None: Ta, Tb = waveguide_cyl.dxes2T(dxes, rmin) ta = (rmin + numpy.cumsum(dxes[0][0])) / rmin - tb = (rmin + dxes[0][0] / 2 + numpy.cumsum(dxes[1][0])) / rmin + tb = ( + rmin + + dxes[0][0] / 2 + + numpy.concatenate((numpy.zeros(1), numpy.cumsum(dxes[1][0][:-1]))) + ) / rmin numpy.testing.assert_allclose(Ta.diagonal(), numpy.repeat(ta, dxes[0][1].size)) numpy.testing.assert_allclose(Tb.diagonal(), numpy.repeat(tb, dxes[1][1].size)) +@pytest.mark.parametrize('use_mu', [False, True]) +def test_waveguide_cyl_operator_matches_2d_limit(use_mu: bool) -> None: + dxes, epsilon, _rmin = build_waveguide_cyl_fixture(asymmetric=True) + mu = build_waveguide_cyl_mu_profile() if use_mu else None + rmin = 1e15 + + cyl_operator = waveguide_cyl.cylindrical_operator(OMEGA, dxes, epsilon, rmin=rmin, mu=mu) + straight_operator = waveguide_2d.operator_e(OMEGA, dxes, epsilon, mu=mu) + + numpy.testing.assert_allclose( + cyl_operator.toarray(), + straight_operator.toarray(), + rtol=1e-9, + atol=1e-10, + ) + + +@pytest.mark.parametrize('use_mu', [False, True]) +def test_waveguide_cyl_reconstruction_matches_2d_limit(use_mu: bool) -> None: + dxes, epsilon, _rmin = build_waveguide_cyl_fixture(asymmetric=True) + mu = build_waveguide_cyl_mu_profile() if use_mu else None + rmin = 1e15 + e_xy, wavenumber = waveguide_2d.solve_mode( + 0, + omega=OMEGA, + dxes=dxes, + epsilon=epsilon, + mu=mu, + ) + angular_wavenumber = wavenumber * rmin + + cyl_e = waveguide_cyl.exy2e( + angular_wavenumber=angular_wavenumber, + dxes=dxes, + rmin=rmin, + epsilon=epsilon, + ) @ e_xy + straight_e = waveguide_2d.exy2e( + wavenumber=wavenumber, + dxes=dxes, + epsilon=epsilon, + ) @ e_xy + cyl_h = waveguide_cyl.e2h( + angular_wavenumber=angular_wavenumber, + omega=OMEGA, + dxes=dxes, + rmin=rmin, + mu=mu, + ) @ cyl_e + straight_h = waveguide_2d.e2h( + wavenumber=wavenumber, + omega=OMEGA, + dxes=dxes, + mu=mu, + ) @ straight_e + + numpy.testing.assert_allclose(cyl_e, straight_e, rtol=1e-8, atol=1e-8) + numpy.testing.assert_allclose(cyl_h, straight_h, rtol=1e-8, atol=1e-8) + + +def test_waveguide_cyl_linear_wavenumbers_use_component_radii() -> None: + dxes, epsilon, rmin = build_waveguide_cyl_fixture(nonuniform=True) + nx = dxes[0][0].size + ny = dxes[0][1].size + angular_wavenumber = numpy.array([2.0]) + + ra = rmin + numpy.cumsum(dxes[0][0]) + rb = ( + rmin + + dxes[0][0] / 2 + + numpy.concatenate((numpy.zeros(1), numpy.cumsum(dxes[1][0][:-1]))) + ) + + er_only = numpy.zeros((1, 2 * nx * ny), dtype=complex) + er_only[0] = vec(numpy.array([numpy.ones((nx, ny)), numpy.zeros((nx, ny))])) + ey_only = numpy.zeros_like(er_only) + ey_only[0] = vec(numpy.array([numpy.zeros((nx, ny)), numpy.ones((nx, ny))])) + + er_linear = waveguide_cyl.linear_wavenumbers( + er_only, + angular_wavenumber, + epsilon=epsilon, + dxes=dxes, + rmin=rmin, + ) + ey_linear = waveguide_cyl.linear_wavenumbers( + ey_only, + angular_wavenumber, + epsilon=epsilon, + dxes=dxes, + rmin=rmin, + ) + + numpy.testing.assert_allclose(er_linear[0], angular_wavenumber[0] / rb.mean()) + numpy.testing.assert_allclose(ey_linear[0], angular_wavenumber[0] / ra.mean()) + + def test_waveguide_cyl_exy2e_and_exy2h_return_finite_full_fields() -> None: dxes, epsilon, rmin = build_waveguide_cyl_fixture() - mu = vec(2 * numpy.ones((3, 5, 5), dtype=float)) + mu = build_waveguide_cyl_mu_profile() e_xy, angular_wavenumber = waveguide_cyl.solve_mode( 0, omega=OMEGA, dxes=dxes, epsilon=epsilon, rmin=rmin, + mu=mu, ) e_field = waveguide_cyl.exy2e( angular_wavenumber=angular_wavenumber, - omega=OMEGA, dxes=dxes, rmin=rmin, epsilon=epsilon, @@ -265,13 +376,14 @@ def test_waveguide_cyl_exy2e_and_exy2h_return_finite_full_fields() -> None: @pytest.mark.parametrize('use_mu', [False, True]) def test_waveguide_cyl_normalized_fields_are_unit_norm_and_silent(use_mu: bool) -> None: dxes, epsilon, rmin = build_waveguide_cyl_fixture() - mu = vec(2 * numpy.ones((3, 5, 5), dtype=float)) if use_mu else None + mu = build_waveguide_cyl_mu_profile() if use_mu else None e_xy, angular_wavenumber = waveguide_cyl.solve_mode( 0, omega=OMEGA, dxes=dxes, epsilon=epsilon, rmin=rmin, + mu=mu, ) output = io.StringIO()