diff --git a/python/tests/test_cyl_weights.py b/python/tests/test_cyl_weights.py
new file mode 100644
index 000000000..2f0f10354
--- /dev/null
+++ b/python/tests/test_cyl_weights.py
@@ -0,0 +1,145 @@
+"""Tests that the integration weights are computed properly in cylindrical coordinates.
+
+"""
+
+import unittest
+from typing import Tuple
+import meep as mp
+import numpy as np
+
+
+resolution = 80  # pixels/μm
+n = 2.4  # refractive index of dielectric layer
+wvl = 1.0  # wavelength (in vacuum)
+fcen = 1 / wvl  # center frequency of source/monitor
+
+
+def led_flux(dmat: float, h: float, rpos: float, m: int) -> Tuple[float, float]:
+    """Computes the radiated and total flux of a point source embedded
+       within a dielectric layer above a lossless ground plane.
+
+    Args:
+       dmat: thickness of dielectric layer.
+       h: height of dipole above ground plane as a fraction of dmat.
+       rpos: position of source in radial direction.
+       m: angular φ dependence of the fields exp(imφ).
+
+    Returns:
+       The radiated and total flux as a 2-Tuple.
+    """
+    L = 20  # length of non-PML region in radial direction
+    dair = 1.0  # thickness of air padding
+    dpml = 1.0  # PML thickness
+    sr = L + dpml
+    sz = dmat + dair + dpml
+    cell_size = mp.Vector3(sr, 0, sz)
+
+    boundary_layers = [
+        mp.PML(dpml, direction=mp.R),
+        mp.PML(dpml, direction=mp.Z, side=mp.High),
+    ]
+
+    src_pt = mp.Vector3(rpos, 0, -0.5 * sz + h * dmat)
+    sources = [
+        mp.Source(
+            src=mp.GaussianSource(fcen, fwidth=0.1 * fcen),
+            component=mp.Er,
+            center=src_pt,
+        ),
+    ]
+
+    geometry = [
+        mp.Block(
+            material=mp.Medium(index=n),
+            center=mp.Vector3(0, 0, -0.5 * sz + 0.5 * dmat),
+            size=mp.Vector3(mp.inf, mp.inf, dmat),
+        )
+    ]
+
+    sim = mp.Simulation(
+        resolution=resolution,
+        cell_size=cell_size,
+        dimensions=mp.CYLINDRICAL,
+        m=m,
+        boundary_layers=boundary_layers,
+        sources=sources,
+        geometry=geometry,
+    )
+
+    flux_air_mon = sim.add_flux(
+        fcen,
+        0,
+        1,
+        mp.FluxRegion(
+            center=mp.Vector3(0.5 * L, 0, 0.5 * sz - dpml),
+            size=mp.Vector3(L, 0, 0),
+        ),
+        mp.FluxRegion(
+            center=mp.Vector3(L, 0, 0.5 * sz - dpml - 0.5 * dair),
+            size=mp.Vector3(0, 0, dair),
+        ),
+    )
+
+    sim.run(
+        mp.dft_ldos(fcen, 0, 1),
+        until_after_sources=mp.stop_when_fields_decayed(
+            50.0,
+            mp.Er,
+            src_pt,
+            1e-8,
+        ),
+    )
+
+    flux_air = mp.get_fluxes(flux_air_mon)[0]
+
+    if rpos == 0:
+        dV = np.pi / (resolution**3)
+    else:
+        dV = 2 * np.pi * rpos / (resolution**2)
+
+    # total flux from point source via LDOS
+    flux_src = -np.real(sim.ldos_Fdata[0] * np.conj(sim.ldos_Jdata[0])) * dV
+
+    print(f"flux-cyl:, {rpos:.2f}, {m:3d}, {flux_src:.6f}, {flux_air:.6f}")
+
+    return flux_air, flux_src
+
+class TestCylWeights(unittest.TestCase):
+    def test_weights(self):
+        layer_thickness = 0.7 * wvl / n
+        dipole_height = 0.5
+
+        rpos = [0.5, 1.5]
+        flux_tol = 1e-5  # threshold flux to determine when to truncate expansion
+        normalized_flux = []
+        for rp in rpos:
+            # analytic upper bound on m based on coupling to free-space modes
+            # in light cone of source medium
+            cutoff_M = int(rp * 2 * np.pi * fcen * n)
+            ms = range(cutoff_M + 1)
+            flux_src_tot = 0
+            flux_air_tot = 0
+            flux_air_max = 0
+            for m in ms:
+                flux_air, flux_src = led_flux(
+                    layer_thickness,
+                    dipole_height,
+                    rp,
+                    m,
+                )
+                flux_air_tot += flux_air if m == 0 else 2 * flux_air
+                flux_src_tot += flux_src if m == 0 else 2 * flux_src
+                if flux_air > flux_air_max:
+                    flux_air_max = flux_air
+                if m > 0 and (flux_air / flux_air_max) < flux_tol:
+                    break
+            
+            normalized_flux.append(flux_air_tot / rp**2)
+        
+        print("--------------------------")
+        print(normalized_flux)
+        # TODO
+        # check for equivalence
+
+if __name__ == "__main__":
+    unittest.main()
\ No newline at end of file
diff --git a/src/loop_in_chunks.cpp b/src/loop_in_chunks.cpp
index b8456e585..572ea80a8 100644
--- a/src/loop_in_chunks.cpp
+++ b/src/loop_in_chunks.cpp
@@ -257,39 +257,100 @@ static inline int iabs(int i) { return (i < 0 ? -i : i); }
 void compute_boundary_weights(grid_volume gv, const volume &where, ivec &is, ivec &ie,
                               bool snap_empty_dimensions, vec &s0, vec &e0, vec &s1, vec &e1) {
   LOOP_OVER_DIRECTIONS(gv.dim, d) {
-    double w0, w1;
-    w0 = 1. - where.in_direction_min(d) * gv.a + 0.5 * is.in_direction(d);
-    w1 = 1. + where.in_direction_max(d) * gv.a - 0.5 * ie.in_direction(d);
+      // w0 = (i+1) - a (from above)
+      double w0 = 1. - where.in_direction_min(d) * gv.a + 0.5 * is.in_direction(d);
+      // w1 = b - (i+3) (from above)
+      double w1 = 1. + where.in_direction_max(d) * gv.a - 0.5 * ie.in_direction(d);
+
+      // only needed for d==R cases
+      double idx_i = is.in_direction(d);
+      double idx_in = ie.in_direction(d);
+
+    // Case 1 (a and b separated by at least 2 grid points)
     if (ie.in_direction(d) >= is.in_direction(d) + 3 * 2) {
-      s0.set_direction(d, w0 * w0 / 2);
-      s1.set_direction(d, 1 - (1 - w0) * (1 - w0) / 2);
-      e0.set_direction(d, w1 * w1 / 2);
-      e1.set_direction(d, 1 - (1 - w1) * (1 - w1) / 2);
+      // cylindrical coordinates in R direction require extra care
+      if (d == R) {
+        // include dV in the weight to prevent division by zero
+        if (is.in_direction(R) == 0)
+          s0.set_direction(d, -w0 * w0 * w0 / 3 + w0 * w0 / 2);
+        else
+          s0.set_direction(d, (-w0 * w0 * w0 / 3 + (idx_i+2) * w0 * w0 / 2) / idx_i );
+        
+        s1.set_direction(d, (w0 * w0 * w0 / 3 - (idx_i+4) * w0 * w0 / 2 + (idx_i+2) * w0 + (idx_i+1)/2 + 1.0 / 6.0) / (idx_i+2));
+        e0.set_direction(d, (w1 * w1 * w1 / 3 + (idx_in-2) * w1 * w1 / 2) / (idx_in));
+        e1.set_direction(d, (-w1 * w1 * w1 / 3 - (idx_in-4) * w1 * w1 / 2 + (idx_in-2) * w1 + (idx_in-2)/2 - 1/6) / (idx_in-2));    
+      } else {
+        s0.set_direction(d, w0 * w0 / 2);
+        s1.set_direction(d, 1 - (1 - w0) * (1 - w0) / 2);
+        e0.set_direction(d, w1 * w1 / 2);
+        e1.set_direction(d, 1 - (1 - w1) * (1 - w1) / 2);
+      }
+
     }
+    // Case 2 (one grid point between a and b)
     else if (ie.in_direction(d) == is.in_direction(d) + 2 * 2) {
-      s0.set_direction(d, w0 * w0 / 2);
-      s1.set_direction(d, 1 - (1 - w0) * (1 - w0) / 2 - (1 - w1) * (1 - w1) / 2);
-      e0.set_direction(d, w1 * w1 / 2);
-      e1.set_direction(d, s1.in_direction(d));
+      // cylindrical coordinates in R direction require extra care
+      if (d == R) {
+        // include dV in the weight to prevent division by zero
+        if (is.in_direction(R) == 0)
+          s0.set_direction(d, -w0 * w0 * w0 / 3 + w0 * w0 / 2);
+        else
+          s0.set_direction(d, -(w0 * w0 * w0 / 3 + (idx_i+2) * w0 * w0 / 2) / idx_i );
+        
+        s1.set_direction(d, (w0 * w0 * w0 / 3 - (idx_i+4) * w0 * w0 / 2 + (idx_i+2) * w0 
+          - w1 * w1 * w1 / 3 - (idx_i) * w1 * w1 / 2 + (idx_i+2) * w1) / (idx_i+2));
+        e0.set_direction(d, (w1 * w1 * w1 / 3 + (idx_i+2) * w1 * w1 / 2) / (idx_i+2));
+        e1.set_direction(d, s1.in_direction(d));
+      } else {
+        s0.set_direction(d, w0 * w0 / 2);
+        s1.set_direction(d, 1 - (1 - w0) * (1 - w0) / 2 - (1 - w1) * (1 - w1) / 2);
+        e0.set_direction(d, w1 * w1 / 2);
+        e1.set_direction(d, s1.in_direction(d));
+      }
     }
+    // Case 4 (as (3), but a = b: interpolation, not integration)
+    // to make the logic easier, we check for 4 before 3...
     else if (where.in_direction_min(d) == where.in_direction_max(d)) {
-      if (snap_empty_dimensions) {
-        if (w0 > w1)
-          ie.set_direction(d, is.in_direction(d));
+      // cylindrical coordinates should be the same as cartesian, with the
+      // exception at r=0 (due to our convention when defining dv1).
+        if (snap_empty_dimensions) {
+          if (w0 > w1)
+            ie.set_direction(d, is.in_direction(d));
+          else
+            is.set_direction(d, ie.in_direction(d));
+          w0 = w1 = 1.0;
+        }
+        // include dV in the weight to prevent division by zero
+        if ((d == R) && (is.in_direction(R) == 0))
+          s0.set_direction(d, w0);
         else
-          is.set_direction(d, ie.in_direction(d));
-        w0 = w1 = 1.0;
+          s0.set_direction(d, w0);
+        s1.set_direction(d, w1);
+        e0.set_direction(d, w1);
+        e1.set_direction(d, w0);
       }
-      s0.set_direction(d, w0);
-      s1.set_direction(d, w1);
-      e0.set_direction(d, w1);
-      e1.set_direction(d, w0);
-    }
+    // Case 3 (no grid points between a and b)
     else if (ie.in_direction(d) == is.in_direction(d) + 1 * 2) {
-      s0.set_direction(d, w0 * w0 / 2 - (1 - w1) * (1 - w1) / 2);
-      e0.set_direction(d, w1 * w1 / 2 - (1 - w0) * (1 - w0) / 2);
-      s1.set_direction(d, e0.in_direction(d));
-      e1.set_direction(d, s0.in_direction(d));
+      // cylindrical coordinates in R direction require extra care
+      if (d == R) {
+        // include dV in the weight to prevent division by zero
+        if (is.in_direction(R) == 0)
+          s0.set_direction(d, -w0 * w0 * w0 / 3 + w0 * w0 / 2 - w1 * w1 * w1 / 3 + w1 *w1 /2 - 1/6);
+        else
+          s0.set_direction(d, (-w0 * w0 * w0 / 3 + (idx_i+2) * w0 * w0 / 2 
+          - w1 * w1 * w1 / 3 - (idx_i-2) * w1 * w1 / 2 + idx_i * w1 - idx_i/2 - 1/6) / idx_i );
+        
+        e0.set_direction(d, (w0 * w0 * w0 / 3 - (idx_i+4) * w0 * w0 / 2 + (idx_i+2) * w0
+         + w1 * w1 * w1 /3 + (idx_i) * w1 * w1 / 2 - (idx_i+2) / 2 + 1/6) / (idx_i+2));
+        s1.set_direction(d, e0.in_direction(d));
+        e1.set_direction(d, s0.in_direction(d)); 
+
+      } else {      
+        s0.set_direction(d, w0 * w0 / 2 - (1 - w1) * (1 - w1) / 2);
+        e0.set_direction(d, w1 * w1 / 2 - (1 - w0) * (1 - w0) / 2);
+        s1.set_direction(d, e0.in_direction(d));
+        e1.set_direction(d, s0.in_direction(d));
+      }
     }
     else
       meep::abort("bug: impossible(?) looping boundaries");
@@ -352,9 +413,19 @@ void fields::loop_in_chunks(field_chunkloop chunkloop, void *chunkloop_data, con
   vec yee_c(gv.yee_shift(Centered) - gv.yee_shift(cgrid));
   ivec iyee_c(gv.iyee_shift(Centered) - gv.iyee_shift(cgrid));
   volume wherec(where + yee_c);
+  
+  /* Integrating zero crossings in cylindrical coordinates is really hairy...
+  let's just punt for now...*/
+  if ((gv.dim == Dcyl) && (wherec.get_min_corner().in_direction(R) < 0))
+    meep::abort("Monitors and sources must not extend below r=0.");
+
   ivec is(vec2diel_floor(wherec.get_min_corner(), gv.a, zero_ivec(gv.dim)) - iyee_c);
   ivec ie(vec2diel_ceil(wherec.get_max_corner(), gv.a, zero_ivec(gv.dim)) - iyee_c);
 
+  //In cylindrical coordinates, we don't need the extra pad point beyond r=0
+  if ((gv.dim == Dcyl) && (is.in_direction(R) < 0))
+    is.set_direction(R,is.in_direction(R) + 2); // add a yee voxel width
+
   vec s0(gv.dim), e0(gv.dim), s1(gv.dim), e1(gv.dim);
   compute_boundary_weights(gv, where, is, ie, snap_empty_dimensions, s0, e0, s1, e1);
 
@@ -509,7 +580,12 @@ void fields::loop_in_chunks(field_chunkloop chunkloop, void *chunkloop_data, con
           }
           if (gv.dim == Dcyl) {
             dV1 = dV0 * 2 * pi * gv.inva;
-            dV0 *= 2 * pi * fabs((S.transform(chunks[i]->gv[isc], sn) + shift).in_direction(R));
+            // we include the dV term in the computation of the s0 weight when
+            // r=0 so that we don't have to worry about dividing by zero. 
+            if (is.in_direction(R) == 0)
+              dV0 = 1;
+            else
+              dV0 *= 2 * pi * fabs((S.transform(chunks[i]->gv[isc], sn) + shift).in_direction(R));
           }
 
           chunkloop(chunks[i], i, cS, isc, iec, s0c, s1c, e0c, e1c, dV0, dV1, shifti, ph, S, sn,