From f1b1efdb39214f734f293fc1ff5d496a8e6d1fa1 Mon Sep 17 00:00:00 2001
From: Forgejo Actions <forgejo-actions@localhost>
Date: Thu, 25 Jun 2026 11:52:42 -0700
Subject: [PATCH] [waveguide_cyl] enable mu

---
 meanas/fdfd/waveguide_cyl.py | 22 +++++++++++++++++-----
 1 file changed, 17 insertions(+), 5 deletions(-)

diff --git a/meanas/fdfd/waveguide_cyl.py b/meanas/fdfd/waveguide_cyl.py
index f2cb5c3..e1f8fd0 100644
--- a/meanas/fdfd/waveguide_cyl.py
+++ b/meanas/fdfd/waveguide_cyl.py
@@ -143,6 +143,7 @@ def cylindrical_operator(
         dxes: dx_lists2_t,
         epsilon: vfdslice,
         rmin: float,
+        mu: vfdslice | None = None,
         ) -> sparse.sparray:
     r"""
     Cylindrical coordinate waveguide operator of the form
@@ -176,10 +177,13 @@ def cylindrical_operator(
         dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D)
         epsilon: Vectorized dielectric constant grid
         rmin: Radius at the left edge of the simulation domain (at minimum 'x')
+        mu: Vectorized magnetic permeability grid (default 1 everywhere)
 
     Returns:
         Sparse matrix representation of the operator
     """
+    if mu is None:
+        mu = numpy.ones_like(epsilon)
 
     Dfx, Dfy = deriv_forward(dxes[0])
     Dbx, Dby = deriv_back(dxes[1])
@@ -191,12 +195,17 @@ def cylindrical_operator(
     eps_y = sparse.diags_array(eps_parts[1])
     eps_z_inv = sparse.diags_array(1 / eps_parts[2])
 
+    mu_parts = numpy.split(mu, 3)
+    mu_y = sparse.diags_array(mu_parts[1])
+    mu_x = sparse.diags_array(mu_parts[0])
+    mu_z_inv = sparse.diags_array(1 / mu_parts[2])
+
     omega2 = omega * omega
     diag = sparse.block_diag
 
-    sq0 = omega2 * diag((Tb @ Tb @ eps_x,
-                         Ta @ Ta @ eps_y))
-    lin0 = sparse.vstack((-Tb @ Dby, Ta @ Dbx)) @ Tb @ sparse.hstack((-Dfy, Dfx))
+    sq0 = omega2 * diag((Tb @ Tb @ mu_y @ eps_x,
+                         Ta @ Ta @ mu_x @ eps_y))
+    lin0 = sparse.vstack((-Tb @ mu_y @ Dby, Ta @ mu_x @ Dbx)) @ Tb @ mu_z_inv @ sparse.hstack((-Dfy, Dfx))
     lin1 = sparse.vstack((Dfx, Dfy)) @ Ta @ eps_z_inv @ sparse.hstack((Dbx @ Tb @ eps_x,
                                                                        Dby @ Ta @ eps_y))
     op = sq0 + lin0 + lin1
@@ -209,6 +218,7 @@ def solve_modes(
         dxes: dx_lists2_t,
         epsilon: vfdslice,
         rmin: float,
+        mu: vfdslice | None = None,
         mode_margin: int = 2,
         ) -> tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]:
     """
@@ -223,6 +233,7 @@ def solve_modes(
         epsilon: Dielectric constant
         rmin: Radius of curvature for the simulation. This should be the minimum value of
             r within the simulation domain.
+        mu: Magnetic permeability (default 1 everywhere)
 
     Returns:
         e_xys: NDArray of vfdfield_t specifying fields. First dimension is mode number.
@@ -233,8 +244,9 @@ def solve_modes(
     # Solve for the largest-magnitude eigenvalue of the real operator
     #
     dxes_real = [[numpy.real(dx) for dx in dxi] for dxi in dxes]
+    mu_real = None if mu is None else numpy.real(mu)
 
-    A_r = cylindrical_operator(numpy.real(omega), dxes_real, numpy.real(epsilon), rmin=rmin)
+    A_r = cylindrical_operator(numpy.real(omega), dxes_real, numpy.real(epsilon), rmin=rmin, mu=mu_real)
     eigvals, eigvecs = signed_eigensolve(A_r, max(mode_numbers) + mode_margin)
     keep_inds = -(numpy.array(mode_numbers) + 1)
     e_xys = eigvecs[:, keep_inds].T
@@ -244,7 +256,7 @@ def solve_modes(
     # Now solve for the eigenvector of the full operator, using the real operator's
     #  eigenvector as an initial guess for Rayleigh quotient iteration.
     #
-    A = cylindrical_operator(omega, dxes, epsilon, rmin=rmin)
+    A = cylindrical_operator(omega, dxes, epsilon, rmin=rmin, mu=mu)
     for nn in range(len(mode_numbers)):
         eigvals[nn], e_xys[nn, :] = rayleigh_quotient_iteration(A, e_xys[nn, :])