[waveguide_cyl] adjust operators to match waveguide_2d in the large-radius limit

meanas/fdfd/waveguide_cyl.py

@@ -43,39 +43,9 @@ T_b &= \operatorname{diag}(r_b / r_{\min}).
 $$
 
 With the same forward/backward derivative notation used in `waveguide_2d`, the
-coordinate-transformed discrete curl equations used here are
+implementation treats the transverse electric eigenproblem as the canonical
-
+cylindrical discretization. It reduces to `waveguide_2d.operator_e(...)` in the
-$$
+large-radius limit `T_a, T_b \to I`. The eigenproblem implemented by
-\begin{aligned}
--\imath \omega \mu_{rr} H_r &= \tilde{\partial}_y E_z + \imath \beta T_a^{-1} E_y, \\
--\imath \omega \mu_{yy} H_y &= -\imath \beta T_b^{-1} E_r
-  - T_b^{-1} \tilde{\partial}_r (T_a E_z), \\
--\imath \omega \mu_{zz} H_z &= \tilde{\partial}_r E_y - \tilde{\partial}_y E_r, \\
-\imath \beta H_y &= -\imath \omega T_b \epsilon_{rr} E_r - T_b \hat{\partial}_y H_z, \\
-\imath \beta H_r &= \imath \omega T_a \epsilon_{yy} E_y
-  - T_b T_a^{-1} \hat{\partial}_r (T_b H_z), \\
-\imath \omega E_z &= T_a \epsilon_{zz}^{-1}
-  \left(\hat{\partial}_r H_y - \hat{\partial}_y H_r\right).
-\end{aligned}
-$$
-
-The first three equations are the cylindrical analogue of the straight-guide
-relations for `H_r`, `H_y`, and `H_z`. The next two are the metric-weighted
-versions of the straight-guide identities for `\imath \beta H_y` and
-`\imath \beta H_r`, and the last equation plays the same role as the
-longitudinal `E_z` reconstruction in `waveguide_2d`.
-
-Following the same elimination steps as in `waveguide_2d`, apply
-`\imath \beta \tilde{\partial}_r` and `\imath \beta \tilde{\partial}_y` to the
-equation for `E_z`, substitute for `\imath \beta H_r` and `\imath \beta H_y`,
-and then eliminate `H_z` with
-
-$$
-H_z = \frac{1}{-\imath \omega \mu_{zz}}
-\left(\tilde{\partial}_r E_y - \tilde{\partial}_y E_r\right).
-$$
-
-This yields the transverse electric eigenproblem implemented by
 `cylindrical_operator(...)`:
 
 $$
@@ -111,6 +81,33 @@ T_a \epsilon_{zz}^{-1}
 \begin{bmatrix} E_r \\ E_y \end{bmatrix}.
 $$
 
+The full fields reconstructed by `exy2e(...)` and `e2h(...)` use the matching
+large-radius-compatible identities
+
+$$
+E_z =
+\frac{1}{\imath \beta} T_a \epsilon_{zz}^{-1}
+\begin{bmatrix}
+\hat{\partial}_r T_b \epsilon_{rr} &
+\hat{\partial}_y T_a \epsilon_{yy}
+\end{bmatrix}
+\begin{bmatrix} E_r \\ E_y \end{bmatrix},
+$$
+
+and
+
+$$
+\begin{bmatrix} H_r \\ H_y \\ H_z \end{bmatrix}
+=
+\frac{1}{-\imath \omega}\mu^{-1}
+\begin{bmatrix}
+0 & \imath\beta T_a^{-1} & \tilde{\partial}_y \\
+-\imath\beta T_b^{-1} & 0 & -T_b^{-1}\tilde{\partial}_r T_a \\
+-\tilde{\partial}_y & \tilde{\partial}_r & 0
+\end{bmatrix}
+\begin{bmatrix} E_r \\ E_y \\ E_z \end{bmatrix}.
+$$
+
 Since `\beta = m / r_{\min}`, the solver implemented in this file returns the
 angular wavenumber `m`, while the operator itself is most naturally written in
 terms of the linear quantity `\beta`. The helpers below reconstruct the full
@@ -367,7 +364,6 @@ def exy2h(
 
 def exy2e(
         angular_wavenumber: complex,
-        omega: float,
         dxes: dx_lists2_t,
         rmin: float,
         epsilon: vfdslice,
@@ -383,7 +379,6 @@ def exy2e(
         angular_wavenumber: Wavenumber assuming fields have theta-dependence of
             `exp(-i * angular_wavenumber * theta)`. It should satisfy
             `operator_e() @ e_xy == (angular_wavenumber / rmin) ** 2 * e_xy`
-        omega: The angular frequency of the system
         dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D)
         rmin: Radius at the left edge of the simulation domain (at minimum 'x')
         epsilon: Vectorized dielectric constant grid
@@ -391,30 +386,22 @@ def exy2e(
     Returns:
         Sparse matrix representing the operator.
     """
-    Dfx, Dfy = deriv_forward(dxes[0])
     Dbx, Dby = deriv_back(dxes[1])
     wavenumber = angular_wavenumber / rmin
 
     Ta, Tb = dxes2T(dxes=dxes, rmin=rmin)
-    Tai = sparse.diags_array(1 / Ta.diagonal())
-    #Tbi = sparse.diags_array(1 / Tb.diagonal())
 
     epsilon_parts = numpy.split(epsilon, 3)
     epsilon_x, epsilon_y = (sparse.diags_array(epsi) for epsi in epsilon_parts[:2])
     epsilon_z_inv = sparse.diags_array(1 / epsilon_parts[2])
 
     n_pts = dxes[0][0].size * dxes[0][1].size
-    zeros = sparse.coo_array((n_pts, n_pts))
+    exy2ez = (
-
+        Ta @ epsilon_z_inv
-    mu_z = numpy.ones(n_pts)
+        @ sparse.hstack((Dbx @ Tb @ epsilon_x,
-    mu_z_inv = sparse.diags_array(1 / mu_z)
+                         Dby @ Ta @ epsilon_y))
-    exy2hz = 1 / (-1j * omega) * mu_z_inv @ sparse.hstack((Dfy, -Dfx))
+        / (1j * wavenumber)
-    hxy2ez = 1 / (1j * omega) * epsilon_z_inv @ sparse.hstack((Dby, -Dbx))
+        )
-
-    exy2hy = Tb / (1j * wavenumber) @ (-1j * omega * sparse.hstack((epsilon_x, zeros)) - Dby @ exy2hz)
-    exy2hx = Tb / (1j * wavenumber) @ ( 1j * omega * sparse.hstack((zeros, epsilon_y)) - Tai @ Dbx @ Tb @ exy2hz)
-
-    exy2ez = hxy2ez @ sparse.vstack((exy2hx, exy2hy))
 
     op = sparse.vstack((sparse.eye_array(2 * n_pts),
                         exy2ez))
@@ -460,9 +447,9 @@ def e2h(
     Tbi = sparse.diags_array(1 / Tb.diagonal())
 
     jB = 1j * angular_wavenumber / rmin
-    op = sparse.block_array([[    None, -jB * Tai,           -Dfy],
+    op = sparse.block_array([[     None,  jB * Tai,            Dfy],
-                             [jB * Tbi,      None, Tbi @ Dfx @ Ta],
+                             [-jB * Tbi,       None, -Tbi @ Dfx @ Ta],
-                             [     Dfy,      -Dfx,           None]]) / (-1j * omega)
+                             [     -Dfy,        Dfx,           None]]) / (-1j * omega)
     if mu is not None:
         op = sparse.diags_array(1 / mu) @ op
     return op
@@ -565,19 +552,16 @@ def _normalized_fields(
 
     The normalization procedure is:
 
-    1. Flip the magnetic field sign so the reconstructed `(e, h)` pair follows
+    1. Compute the time-averaged forward power with
-       the same forward-power convention as `waveguide_2d`.
-    2. Compute the time-averaged forward power with
        `waveguide_2d.inner_product(..., conj_h=True)`.
-    3. Scale by `1 / sqrt(S_z)` so the resulting mode has unit forward power.
+    2. Scale by `1 / sqrt(S_z)` so the resulting mode has unit forward power.
-    4. Remove the arbitrary complex phase and apply a quadrant-sum sign heuristic
+    3. Remove the arbitrary complex phase and apply a quadrant-sum sign heuristic
        so symmetric modes choose a stable representative.
 
     `prop_phase` has the same meaning as in `waveguide_2d`: it compensates for
     the half-cell longitudinal staggering between the E and H fields when the
     propagation direction is itself discretized.
     """
-    h *= -1
     shape = [s.size for s in dxes[0]]
 
     # Find time-averaged Sz and normalize to it