diff --git a/meanas/fdfd/operators.py b/meanas/fdfd/operators.py
index afa5fbd..8c16ef7 100644
--- a/meanas/fdfd/operators.py
+++ b/meanas/fdfd/operators.py
@@ -40,7 +40,7 @@ __author__ = 'Jan Petykiewicz'
 def e_full(
         omega: complex,
         dxes: dx_lists_t,
-        epsilon: vfdfield_t,
+        epsilon: vfdfield_t | vcfdfield_t,
         mu: vfdfield_t | None = None,
         pec: vfdfield_t | None = None,
         pmc: vfdfield_t | None = None,
diff --git a/meanas/fdfd/solvers.py b/meanas/fdfd/solvers.py
index 517ecab..5b48493 100644
--- a/meanas/fdfd/solvers.py
+++ b/meanas/fdfd/solvers.py
@@ -35,9 +35,9 @@ def _scipy_qmr(
         Guess for solution (returned even if didn't converge)
     """
 
-    '''
-    Report on our progress
-    '''
+    #
+    #Report on our progress
+    #
     ii = 0
 
     def log_residual(xk: ArrayLike) -> None:
@@ -56,10 +56,9 @@ def _scipy_qmr(
     else:
         kwargs['callback'] = log_residual
 
-    '''
-    Run the actual solve
-    '''
-
+    #
+    # Run the actual solve
+    #
     x, _ = scipy.sparse.linalg.qmr(A, b, **kwargs)
     return x
 
diff --git a/meanas/fdfd/waveguide_2d.py b/meanas/fdfd/waveguide_2d.py
index 32e65bc..05215c4 100644
--- a/meanas/fdfd/waveguide_2d.py
+++ b/meanas/fdfd/waveguide_2d.py
@@ -845,13 +845,13 @@ def solve_modes(
             ability to find the correct mode. Default 2.
 
     Returns:
-        e_xys: list of vfdfield_t specifying fields
+        e_xys: NDArray of vfdfield_t specifying fields. First dimension is mode number.
         wavenumbers: list of wavenumbers
     """
 
-    '''
-    Solve for the largest-magnitude eigenvalue of the real operator
-    '''
+    #
+    # 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 = operator_e(numpy.real(omega), dxes_real, numpy.real(epsilon), mu_real)
@@ -859,10 +859,10 @@ def solve_modes(
     eigvals, eigvecs = signed_eigensolve(A_r, max(mode_numbers) + mode_margin)
     e_xys = eigvecs[:, -(numpy.array(mode_numbers) + 1)]
 
-    '''
-    Now solve for the eigenvector of the full operator, using the real operator's
-     eigenvector as an initial guess for Rayleigh quotient iteration.
-    '''
+    #
+    # Now solve for the eigenvector of the full operator, using the real operator's
+    #  eigenvector as an initial guess for Rayleigh quotient iteration.
+    #
     A = operator_e(omega, dxes, epsilon, mu)
     for nn in range(len(mode_numbers)):
         eigvals[nn], e_xys[:, nn] = rayleigh_quotient_iteration(A, e_xys[:, nn])
diff --git a/meanas/fdfd/waveguide_3d.py b/meanas/fdfd/waveguide_3d.py
index 3cffa94..8bb0513 100644
--- a/meanas/fdfd/waveguide_3d.py
+++ b/meanas/fdfd/waveguide_3d.py
@@ -53,9 +53,9 @@ def solve_mode(
 
     slices = tuple(slices)
 
-    '''
-    Solve the 2D problem in the specified plane
-    '''
+    #
+    # Solve the 2D problem in the specified plane
+    #
     # Define rotation to set z as propagation direction
     order = numpy.roll(range(3), 2 - axis)
     reverse_order = numpy.roll(range(3), axis - 2)
@@ -73,9 +73,10 @@ def solve_mode(
     }
     e_xy, wavenumber_2d = waveguide_2d.solve_mode(mode_number, **args_2d)
 
-    '''
-    Apply corrections and expand to 3D
-    '''
+    #
+    # Apply corrections and expand to 3D
+    #
+
     # Correct wavenumber to account for numerical dispersion.
     wavenumber = 2 / dx_prop * numpy.arcsin(wavenumber_2d * dx_prop / 2)
 
diff --git a/meanas/fdmath/types.py b/meanas/fdmath/types.py
index bc678ea..d44b30a 100644
--- a/meanas/fdmath/types.py
+++ b/meanas/fdmath/types.py
@@ -20,7 +20,7 @@ vcfdfield_t = NDArray[complexfloating]
 """Linearized complex vector field (single vector of length 3*X*Y*Z)"""
 
 
-dx_lists_t = Sequence[Sequence[NDArray[floating]]]
+dx_lists_t = Sequence[Sequence[NDArray[floating | complexfloating]]]
 """
  'dxes' datastructure which contains grid cell width information in the following format:
 
@@ -31,7 +31,7 @@ dx_lists_t = Sequence[Sequence[NDArray[floating]]]
    and `dy_h[0]` is the y-width of the `y=0` cells, as used when calculating dH/dy, etc.
 """
 
-dx_lists_mut = MutableSequence[MutableSequence[NDArray[floating]]]
+dx_lists_mut = MutableSequence[MutableSequence[NDArray[floating | complexfloating]]]
 """Mutable version of `dx_lists_t`"""