From c249e78cd91294d86f975d533b163623c0ebf735 Mon Sep 17 00:00:00 2001
From: Forgejo Actions <forgejo-actions@localhost>
Date: Thu, 25 Jun 2026 12:01:36 -0700
Subject: [PATCH] [waveugide_cyl] adjust radius calculation

---
 meanas/fdfd/waveguide_cyl.py | 23 +++++++++++++++++++----
 1 file changed, 19 insertions(+), 4 deletions(-)

diff --git a/meanas/fdfd/waveguide_cyl.py b/meanas/fdfd/waveguide_cyl.py
index 96a2ff6..a00fbf9 100644
--- a/meanas/fdfd/waveguide_cyl.py
+++ b/meanas/fdfd/waveguide_cyl.py
@@ -321,12 +321,20 @@ def linear_wavenumbers(
 
     shape2d = (len(dxes[0][0]), len(dxes[0][1]))
     epsilon2d = unvec(epsilon, shape2d)[:2]
-    grid_radii = rmin + numpy.cumsum(dxes[0][0])
+    ra = rmin + numpy.cumsum(dxes[0][0])
+    rb = rmin + dxes[0][0] / 2.0 + numpy.concatenate((
+        numpy.zeros(1, dtype=dxes[1][0].dtype),
+        numpy.cumsum(dxes[1][0][:-1]),
+        ))
     for ii in range(angular_wavenumbers.size):
         efield = unvec(e_xys[ii], shape2d, 2)
         energy = numpy.real((efield * efield.conj()) * epsilon2d)
-        energy_vs_x = energy.sum(axis=(0, 2))
-        mode_radii[ii] = (grid_radii * energy_vs_x).sum() / energy_vs_x.sum()
+        er_energy_vs_r = energy[0].sum(axis=1)
+        ey_energy_vs_r = energy[1].sum(axis=1)
+        energy_vs_r = er_energy_vs_r + ey_energy_vs_r
+        mode_radii[ii] = (
+            (rb * er_energy_vs_r).sum() + (ra * ey_energy_vs_r).sum()
+            ) / energy_vs_r.sum()
 
     logger.info(f'{mode_radii=}')
     lin_wavenumbers = angular_wavenumbers / mode_radii
@@ -474,7 +482,14 @@ def dxes2T(
         Sparse diagonal matrices `(T_a, T_b)`.
     """
     ra = rmin + numpy.cumsum(dxes[0][0])                      # Radius at Ey points
-    rb = rmin + dxes[0][0] / 2.0 + numpy.cumsum(dxes[1][0])   # Radius at Ex points
+    rb = (
+        rmin
+        + dxes[0][0] / 2.0
+        + numpy.concatenate((
+            numpy.zeros(1, dtype=dxes[1][0].dtype),
+            numpy.cumsum(dxes[1][0][:-1]),
+            ))
+        )                                                     # Radius at Er points
     ta = ra / rmin
     tb = rb / rmin