Fix bug in step-db update equations in cylindrical coordinates for `m…

…=±1` (#2382)

* fix bug in step-db update equations in cylindrical coordinates for m=±1

* verify correctness of m=±1 simulation using flux rather than LDOS

python/tests/test_ldos.py

@@ -1,7 +1,5 @@
 import unittest
-
 import numpy as np
-
 import meep as mp
 
 
@@ -17,11 +15,6 @@ def setUp(cls):
 
         cls.fcen = 1 / cls.wvl
 
-        # source properties (cylindrical)
-        cls.df = 0.05 * cls.fcen
-        cls.cutoff = 10.0
-        cls.src = mp.GaussianSource(cls.fcen, fwidth=cls.df, cutoff=cls.cutoff)
-
         # termination criteria
         cls.tol = 1e-8
 
@@ -37,7 +30,7 @@ def bulk_ldos_cyl(self):
 
         sources = [
             mp.Source(
-                src=self.src,
+                src=mp.GaussianSource(self.fcen, fwidth=0.1 * self.fcen),
                 component=mp.Er,
                 center=mp.Vector3(),
             )
@@ -76,7 +69,7 @@ def cavity_ldos_cyl(self, sz):
 
         sources = [
             mp.Source(
-                src=self.src,
+                src=mp.GaussianSource(self.fcen, fwidth=0.1 * self.fcen),
                 component=mp.Er,
                 center=mp.Vector3(),
             )
@@ -201,14 +194,15 @@ def purcell_enh_theory(self, c):
             4 * np.power(np.fix(c + 0.5), 3) - np.fix(c + 0.5)
         ) / (16 * np.power(c, 3))
 
-    def ext_eff_cyl(self, dmat, h):
+    def ext_eff_cyl(self, dmat, h, m):
         """Computes the extraction efficiency of a point dipole embedded
         within a dielectric layer above a lossless ground plane in
         cylindrical coordinates.
 
         Args:
           dmat: thickness of dielectric layer.
           h: height of dipole above ground plane as fraction of dmat.
+          m: angular dependence of the fields, exp(imφ).
         """
         sr = self.L + self.dpml
         sz = dmat + self.dair + self.dpml
@@ -221,7 +215,13 @@ def ext_eff_cyl(self, dmat, h):
 
         src_cmpt = mp.Er
         src_pt = mp.Vector3(0, 0, -0.5 * sz + h * dmat)
-        sources = [mp.Source(src=self.src, component=src_cmpt, center=src_pt)]
+        sources = [
+            mp.Source(
+                src=mp.GaussianSource(self.fcen, fwidth=0.1 * self.fcen),
+                component=src_cmpt,
+                center=src_pt,
+            ),
+        ]
 
         geometry = [
             mp.Block(
@@ -235,7 +235,7 @@ def ext_eff_cyl(self, dmat, h):
             resolution=self.resolution,
             cell_size=cell_size,
             dimensions=mp.CYLINDRICAL,
-            m=-1,
+            m=m,
             boundary_layers=boundary_layers,
             sources=sources,
             geometry=geometry,
@@ -428,15 +428,22 @@ def test_ldos_3D(self):
     def test_ldos_ext_eff(self):
         """Verifies that the extraction efficiency of a point dipole in a
         dielectric layer above a lossless ground plane computed in cylindrical
-        and 3D Cartesian coordinates agree.
+        coordinates (for m=±1, separately) and 3D Cartesian agree.
         """
         layer_thickness = 0.5 * self.wvl / self.n
         dipole_height = 0.5
 
-        ext_eff_cyl = self.ext_eff_cyl(layer_thickness, dipole_height)
+        ext_eff_cyl = self.ext_eff_cyl(layer_thickness, dipole_height, -1.0)
         ext_eff_3D = self.ext_eff_3D(layer_thickness, dipole_height)
 
-        self.assertAlmostEqual(ext_eff_cyl, ext_eff_3D, places=2)
+        self.assertAlmostEqual(ext_eff_cyl, ext_eff_3D, delta=0.01)
+
+        ext_eff_cyl_m_plus = self.ext_eff_cyl(
+            layer_thickness,
+            dipole_height,
+            +1.0,
+        )
+        self.assertEqual(ext_eff_cyl, ext_eff_cyl_m_plus)
 
 
 if __name__ == "__main__":

src/step_db.cpp

@@ -329,7 +329,7 @@ bool fields_chunk::step_db(field_type ft) {
         realnum *fu = siginvu && f_u[cc][cmp] ? f[cc][cmp] : 0;
         realnum *the_f = fu ? f_u[cc][cmp] : f[cc][cmp];
         int sd = ft == D_stuff ? +1 : -1;
-        realnum f_m_mult = ft == D_stuff ? 2 : (1 - 2 * cmp);
+        realnum f_m_mult = ft == D_stuff ? 2 : (1 - 2 * cmp) * m;
 
         for (int iz = (ft == D_stuff); iz < nz + (ft == D_stuff); ++iz) {
           realnum df;