Merge pull request #80 from mfschubert/recon2

improve longitudinal field calculation

Author: mfschubert

CHANGELOG.md

@@ -3,6 +3,7 @@
 ## Unreleased
 - Restructure codebase to integrate code originally in the `_orig` submodule within the main `src/fmmax` modules, separating the code available under the two licenses to the top (AGPL) and bottom (MIT) of each file.
 - Include `transverse_permittivity_matrix` in the `LayerSolveResult`.
+- Update the manner in which the `ez` and `hz` field components are computed. Instead of relying on the Toeplitz matrix obtained from the inverse of permittivity (permeability), they now make use of the inverse of the Toeplitz matrix obtained from the permittiivity (permeability).
 
 ## 1.4.2 (July 10, 2025)
 - Avoid use of `jax.pure_callback` in the eigendecomposition. This requires jax newer than 0.4.36.

src/fmmax/fields.py

@@ -529,21 +529,16 @@ def fields_from_wave_amplitudes(
     )
     angular_frequency = angular_frequency[..., jnp.newaxis, jnp.newaxis]
 
-    # We use the the Fourier convolution matrix for the inverse of permittivity,
-    # rather than inverting the Fourier convolution matrix of permittivity itself.
-    # This improves convergence of the computed z-component of electric field.
-    ez = (
-        -layer_solve_result.inverse_z_permittivity_matrix
-        @ (1j * kx * hy - 1j * ky * hx)
-        / (1j * angular_frequency)
+    ez = jnp.linalg.solve(
+        layer_solve_result.z_permittivity_matrix,
+        -(1j * kx * hy - 1j * ky * hx) / (1j * angular_frequency),
     )
 
     # Compute the z-directed magnetic field. The expression is similar to
     # equation 14 from [2012 Liu], but modified to allow for magnetic materials.
-    hz = (
-        layer_solve_result.inverse_z_permeability_matrix
-        @ (1j * kx * ey - 1j * ky * ex)
-        / (1j * angular_frequency)
+    hz = jnp.linalg.solve(
+        layer_solve_result.z_permeability_matrix,
+        (1j * kx * ey - 1j * ky * ex) / (1j * angular_frequency),
     )
 
     assert ex.shape == ey.shape == ez.shape == hx.shape == hy.shape == hz.shape

tests/examples/test_crystal.py

@@ -36,7 +36,7 @@ def test_regression(self):
         with self.subTest("efield"):
             onp.testing.assert_allclose(
                 [onp.mean(onp.abs(ex)), onp.mean(onp.abs(ey)), onp.mean(onp.abs(ez))],
-                [8.311728e00, 1.293643e-06, 5.980783e00],
+                [8.311728e00, 1.293643e-06, 5.886467e00],
                 rtol=1e-4,
             )
         with self.subTest("hfield"):
@@ -70,7 +70,7 @@ def test_CW_regression(self):
         with self.subTest("efield"):
             onp.testing.assert_allclose(
                 [onp.mean(onp.abs(ex)), onp.mean(onp.abs(ey)), onp.mean(onp.abs(ez))],
-                [4.830812e-01, 4.060432e-04, 3.205514e-01],
+                [4.830812e-01, 4.060432e-04, 3.208762e-01],
                 rtol=1e-4,
                 atol=1e-6,
             )
@@ -111,7 +111,7 @@ def test_broadband_regression(self):
                     onp.mean(onp.abs(ey[wavelength_idx, ...])),
                     onp.mean(onp.abs(ez[wavelength_idx, ...])),
                 ],
-                [4.830811e-01, 4.060428e-04, 3.205514e-01],
+                [4.830811e-01, 4.060428e-04, 3.208762e-01],
                 rtol=1e-4,
             )
         with self.subTest("hfield"):

tests/examples/test_metal_dipole.py

@@ -39,4 +39,4 @@ def test_regression(self):
         with self.subTest("hy"):
             onp.testing.assert_allclose(onp.mean(onp.abs(hy)), 0.0, atol=1e-4)
         with self.subTest("hz"):
-            onp.testing.assert_allclose(onp.mean(onp.abs(hz)), 9.952436, atol=1e-4)
+            onp.testing.assert_allclose(onp.mean(onp.abs(hz)), 9.893226, atol=1e-4)

tests/examples/test_metal_pillars.py

@@ -52,7 +52,7 @@ def test_compute_fields_regression(self):
             onp.testing.assert_allclose(
                 jnp.mean(jnp.abs(efields) ** 2, axis=(2, 3, 4, 5)),
                 onp.asarray(
-                    [[0.046316, 0.052124], [1.047591, 1.859426], [0.145127, 0.077935]]
+                    [[0.046316, 0.052124], [1.047591, 1.859426], [0.322618, 0.095657]]
                 ),
                 rtol=0.01,
             )

tests/examples/test_uled.py

@@ -43,7 +43,7 @@ def test_regression(self):
                 [
                     [64.48044, 8.389313, 2.740909],
                     [8.386348, 64.499756, 2.740864],
-                    [15.16424, 15.166618, 8.501956],
+                    [15.824419, 15.819216, 8.312174],
                 ],
                 rtol=1e-3,
             )