diff --git a/README.md b/README.md
index 35a190c..5a3f49c 100644
--- a/README.md
+++ b/README.md
@@ -1,11 +1,5 @@
 # fdfd_tools
 
-** DEPRECATED **
-
-The functionality in this module is now provided by [meanas](https://mpxd.net/code/jan/meanas).
-
------------------------
-
 **fdfd_tools** is a python package containing utilities for
 creating and analyzing 2D and 3D finite-difference frequency-domain (FDFD)
 electromagnetic simulations.
diff --git a/examples/bloch.py b/examples/bloch.py
deleted file mode 100644
index 8bbd30e..0000000
--- a/examples/bloch.py
+++ /dev/null
@@ -1,68 +0,0 @@
-import numpy, scipy, gridlock, fdfd_tools
-from fdfd_tools import bloch
-from numpy.linalg import norm
-import logging
-
-logging.basicConfig(level=logging.DEBUG)
-logger = logging.getLogger(__name__)
-
-
-dx = 40
-x_period = 400
-y_period = z_period = 2000
-g = gridlock.Grid([numpy.arange(-x_period/2, x_period/2, dx),
-                   numpy.arange(-1000, 1000, dx),
-                   numpy.arange(-1000, 1000, dx)],
-                  shifts=numpy.array([[0,0,0]]),
-                  initial=1.445**2,
-                  periodic=True)
-
-g.draw_cuboid([0,0,0], [200e8, 220, 220], eps=3.47**2)
-
-#x_period = y_period = z_period = 13000
-#g = gridlock.Grid([numpy.arange(3), ]*3,
-#                  shifts=numpy.array([[0, 0, 0]]),
-#                  initial=2.0**2,
-#                  periodic=True)
-
-g2 = g.copy()
-g2.shifts = numpy.zeros((6,3))
-g2.grids = [numpy.zeros(g.shape) for _ in range(6)]
-
-epsilon = [g.grids[0],] * 3
-reciprocal_lattice = numpy.diag(1e6/numpy.array([x_period, y_period, z_period])) #cols are vectors
-
-#print('Finding k at 1550nm')
-#k, f = bloch.find_k(frequency=1/1550,
-#                    tolerance=(1/1550 - 1/1551),
-#                    direction=[1, 0, 0],
-#                    G_matrix=reciprocal_lattice,
-#                    epsilon=epsilon,
-#                    band=0)
-#
-#print("k={}, f={}, 1/f={}, k/f={}".format(k, f, 1/f, norm(reciprocal_lattice @ k) / f ))
-
-print('Finding f at [0.25, 0, 0]')
-for k0x in [.25]:
-    k0 = numpy.array([k0x, 0, 0])
-
-    kmag = norm(reciprocal_lattice @ k0)
-    tolerance = (1e6/1550) * 1e-4/1.5  # df = f * dn_eff / n
-    logger.info('tolerance {}'.format(tolerance))
-
-    n, v = bloch.eigsolve(4, k0, G_matrix=reciprocal_lattice, epsilon=epsilon, tolerance=tolerance)
-    v2e = bloch.hmn_2_exyz(k0, G_matrix=reciprocal_lattice, epsilon=epsilon)
-    v2h = bloch.hmn_2_hxyz(k0, G_matrix=reciprocal_lattice, epsilon=epsilon)
-    ki = bloch.generate_kmn(k0, reciprocal_lattice, g.shape)
-
-    z = 0
-    e = v2e(v[0])
-    for i in range(3):
-        g2.grids[i] += numpy.real(e[i])
-        g2.grids[i+3] += numpy.imag(e[i])
-
-    f = numpy.sqrt(numpy.real(numpy.abs(n))) # TODO
-    print('k0x = {:3g}\n eigval = {}\n f = {}\n'.format(k0x, n, f))
-    n_eff = norm(reciprocal_lattice @ k0) / f
-    print('kmag/f = n_eff =  {} \n wl = {}\n'.format(n_eff, 1/f ))
-
diff --git a/examples/test.py b/examples/test.py
index a7e1746..a40309f 100644
--- a/examples/test.py
+++ b/examples/test.py
@@ -190,6 +190,7 @@ def test1(solver=generic_solver):
 
     s1x, s2x = poyntings(E)
     pyplot.plot(s1x[0].sum(axis=2).sum(axis=1))
+    pyplot.hold(True)
     pyplot.plot(s2x[0].sum(axis=2).sum(axis=1))
     pyplot.show()
 
diff --git a/fdfd_tools/bloch.py b/fdfd_tools/bloch.py
deleted file mode 100644
index b172e21..0000000
--- a/fdfd_tools/bloch.py
+++ /dev/null
@@ -1,513 +0,0 @@
-'''
-Bloch eigenmode solver/operators
-
-This module contains functions for generating and solving the
- 3D Bloch eigenproblem. The approach is to transform the problem
- into the (spatial) fourier domain, transforming the equation
-   1/mu * curl(1/eps * curl(H)) = (w/c)^2 H
- into
-   conv(1/mu_k, ik x conv(1/eps_k, ik x H_k)) = (w/c)^2 H_k
- where:
-  - the _k subscript denotes a 3D fourier transformed field
-  - each component of H_k corresponds to a plane wave with wavevector k
-  - x is the cross product
-  - conv denotes convolution
-
- Since k and H are orthogonal for each plane wave, we can use each
- k to create an orthogonal basis (k, m, n), with k x m = n, and
- |m| = |n| = 1. The cross products are then simplified with
-
- k @ h = kx hx + ky hy + kz hz = 0 = hk
- h = hk + hm + hn = hm + hn
- k = kk + km + kn = kk = |k|
-
- k x h = (ky hz - kz hy,
-          kz hx - kx hz,
-          kx hy - ky hx)
-       = ((k x h) @ k, (k x h) @ m, (k x h) @ n)_kmn
-       = (0, (m x k) @ h, (n x k) @ h)_kmn         # triple product ordering
-       = (0, kk (-n @ h), kk (m @ h))_kmn          # (m x k) = -|k| n, etc.
-       = |k| (0, -h @ n, h @ m)_kmn
-
- k x h = (km hn - kn hm,
-          kn hk - kk hn,
-          kk hm - km hk)_kmn
-       = (0, -kk hn, kk hm)_kmn
-       = (-kk hn)(mx, my, mz) + (kk hm)(nx, ny, nz)
-       = |k| (hm * (nx, ny, nz) - hn * (mx, my, mz))
-
- where h is shorthand for H_k, (...)_kmn deontes the (k, m, n) basis,
- and e.g. hm is the component of h in the m direction.
-
- We can also simplify conv(X_k, Y_k) as fftn(X * ifftn(Y_k)).
-
- Using these results and storing H_k as h = (hm, hn), we have
-  e_xyz = fftn(1/eps * ifftn(|k| (hm * n - hn * m)))
-  b_mn = |k| (-e_xyz @ n, e_xyz @ m)
-  h_mn = fftn(1/mu * ifftn(b_m * m + b_n * n))
- which forms the operator from the left side of the equation.
-
- We can then use a preconditioned block Rayleigh iteration algorithm, as in
-  SG Johnson and JD Joannopoulos, Block-iterative frequency-domain methods
-  for Maxwell's equations in a planewave basis, Optics Express 8, 3, 173-190 (2001)
- (similar to that used in MPB) to find the eigenvectors for this operator.
-
- ===
-
- Typically you will want to do something like
-
-    recip_lattice = numpy.diag(1/numpy.array(epsilon[0].shape * dx))
-    n, v = bloch.eigsolve(5, k0, recip_lattice, epsilon)
-    f = numpy.sqrt(-numpy.real(n[0]))
-    n_eff = norm(recip_lattice @ k0) / f
-
-    v2e = bloch.hmn_2_exyz(k0, recip_lattice, epsilon)
-    e_field = v2e(v[0])
-
-    k, f = find_k(frequency=1/1550,
-                  tolerance=(1/1550 - 1/1551),
-                  direction=[1, 0, 0],
-                  G_matrix=recip_lattice,
-                  epsilon=epsilon,
-                  band=0)
-
-'''
-
-from typing import List, Tuple, Callable, Dict
-import logging
-import numpy
-from numpy.fft import fftn, ifftn, fftfreq
-import scipy
-import scipy.optimize
-from scipy.linalg import norm
-import scipy.sparse.linalg as spalg
-
-from .eigensolvers import rayleigh_quotient_iteration
-from . import field_t
-
-logger = logging.getLogger(__name__)
-
-
-def generate_kmn(k0: numpy.ndarray,
-                 G_matrix: numpy.ndarray,
-                 shape: numpy.ndarray
-                 ) -> Tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray]:
-    """
-    Generate a (k, m, n) orthogonal basis for each k-vector in the simulation grid.
-
-    :param k0: [k0x, k0y, k0z], Bloch wavevector, in G basis.
-    :param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
-    :param shape: [nx, ny, nz] shape of the simulation grid.
-    :return: (|k|, m, n) where |k| has shape tuple(shape) + (1,)
-            and m, n have shape tuple(shape) + (3,).
-            All are given in the xyz basis (e.g. |k|[0,0,0] = norm(G_matrix @ k0)).
-    """
-    k0 = numpy.array(k0)
-
-    Gi_grids = numpy.meshgrid(*(fftfreq(n, 1/n) for n in shape[:3]), indexing='ij')
-    Gi = numpy.stack(Gi_grids, axis=3)
-
-    k_G = k0[None, None, None, :] - Gi
-    k_xyz = numpy.rollaxis(G_matrix @ numpy.rollaxis(k_G, 3, 2), 3, 2)
-
-    m = numpy.broadcast_to([0, 1, 0], tuple(shape[:3]) + (3,)).astype(float)
-    n = numpy.broadcast_to([0, 0, 1], tuple(shape[:3]) + (3,)).astype(float)
-
-    xy_non0 = numpy.any(k_xyz[:, :, :, 0:1] != 0, axis=3)
-    if numpy.any(xy_non0):
-        u = numpy.cross(k_xyz[xy_non0], [0, 0, 1])
-        m[xy_non0, :] = u / norm(u, axis=1)[:, None]
-
-    z_non0 = numpy.any(k_xyz != 0, axis=3)
-    if numpy.any(z_non0):
-        v = numpy.cross(k_xyz[z_non0], m[z_non0])
-        n[z_non0, :] = v / norm(v, axis=1)[:, None]
-
-    k_mag = norm(k_xyz, axis=3)[:, :, :, None]
-    return k_mag, m, n
-
-
-def maxwell_operator(k0: numpy.ndarray,
-                     G_matrix: numpy.ndarray,
-                     epsilon: field_t,
-                     mu: field_t = None
-                     ) -> Callable[[numpy.ndarray], numpy.ndarray]:
-    """
-    Generate the Maxwell operator
-        conv(1/mu_k, ik x conv(1/eps_k, ik x ___))
-    which is the spatial-frequency-space representation of
-        1/mu * curl(1/eps * curl(___))
-
-    The operator is a function that acts on a vector h_mn of size (2 * epsilon[0].size)
-
-    See the module-level docstring for more information.
-
-    :param k0: Bloch wavevector, [k0x, k0y, k0z].
-    :param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
-    :param epsilon: Dielectric constant distribution for the simulation.
-        All fields are sampled at cell centers (i.e., NOT Yee-gridded)
-    :param mu: Magnetic permability distribution for the simulation.
-        Default None (1 everywhere).
-    :return: Function which applies the maxwell operator to h_mn.
-    """
-
-    shape = epsilon[0].shape + (1,)
-    k_mag, m, n = generate_kmn(k0, G_matrix, shape)
-
-    epsilon = numpy.stack(epsilon, 3)
-    if mu is not None:
-        mu = numpy.stack(mu, 3)
-
-    def operator(h: numpy.ndarray):
-        """
-        Maxwell operator for Bloch eigenmode simulation.
-
-        h is complex 2-field in (m, n) basis, vectorized
-
-        :param h: Raveled h_mn; size (2 * epsilon[0].size).
-        :return: Raveled conv(1/mu_k, ik x conv(1/eps_k, ik x h_mn)).
-        """
-        hin_m, hin_n = [hi.reshape(shape) for hi in numpy.split(h, 2)]
-
-        #{d,e,h}_xyz fields are complex 3-fields in (1/x, 1/y, 1/z) basis
-
-        # cross product and transform into xyz basis
-        d_xyz = (n * hin_m -
-                 m * hin_n) * k_mag
-
-        # divide by epsilon
-        e_xyz = fftn(ifftn(d_xyz, axes=range(3)) / epsilon, axes=range(3))
-
-        # cross product and transform into mn basis
-        b_m = numpy.sum(e_xyz * n, axis=3)[:, :, :, None] * -k_mag
-        b_n = numpy.sum(e_xyz * m, axis=3)[:, :, :, None] * +k_mag
-
-        if mu is None:
-            h_m, h_n = b_m, b_n
-        else:
-            # transform from mn to xyz
-            b_xyz = (m * b_m[:, :, :, None] +
-                     n * b_n[:, :, :, None])
-
-            # divide by mu
-            h_xyz = fftn(ifftn(b_xyz, axes=range(3)) / mu, axes=range(3))
-
-            # transform back to mn
-            h_m = numpy.sum(h_xyz * m, axis=3)
-            h_n = numpy.sum(h_xyz * n, axis=3)
-        return numpy.hstack((h_m.ravel(), h_n.ravel()))
-
-    return operator
-
-
-def hmn_2_exyz(k0: numpy.ndarray,
-               G_matrix: numpy.ndarray,
-               epsilon: field_t,
-               ) -> Callable[[numpy.ndarray], field_t]:
-    """
-    Generate an operator which converts a vectorized spatial-frequency-space
-     h_mn into an E-field distribution, i.e.
-        ifft(conv(1/eps_k, ik x h_mn))
-
-    The operator is a function that acts on a vector h_mn of size (2 * epsilon[0].size)
-
-    See the module-level docstring for more information.
-
-    :param k0: Bloch wavevector, [k0x, k0y, k0z].
-    :param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
-    :param epsilon: Dielectric constant distribution for the simulation.
-        All fields are sampled at cell centers (i.e., NOT Yee-gridded)
-    :return: Function for converting h_mn into E_xyz
-    """
-    shape = epsilon[0].shape + (1,)
-    epsilon = numpy.stack(epsilon, 3)
-
-    k_mag, m, n = generate_kmn(k0, G_matrix, shape)
-
-    def operator(h: numpy.ndarray) -> field_t:
-        hin_m, hin_n = [hi.reshape(shape) for hi in numpy.split(h, 2)]
-        d_xyz = (n * hin_m -
-                 m * hin_n) * k_mag
-
-        # divide by epsilon
-        return [ei for ei in numpy.rollaxis(ifftn(d_xyz, axes=range(3)) / epsilon, 3)]
-
-    return operator
-
-
-def hmn_2_hxyz(k0: numpy.ndarray,
-               G_matrix: numpy.ndarray,
-               epsilon: field_t
-               ) -> Callable[[numpy.ndarray], field_t]:
-    """
-    Generate an operator which converts a vectorized spatial-frequency-space
-     h_mn into an H-field distribution, i.e.
-        ifft(h_mn)
-
-    The operator is a function that acts on a vector h_mn of size (2 * epsilon[0].size)
-
-    See the module-level docstring for more information.
-
-    :param k0: Bloch wavevector, [k0x, k0y, k0z].
-    :param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
-    :param epsilon: Dielectric constant distribution for the simulation.
-        Only epsilon[0].shape is used.
-    :return: Function for converting h_mn into H_xyz
-    """
-    shape = epsilon[0].shape + (1,)
-    k_mag, m, n = generate_kmn(k0, G_matrix, shape)
-
-    def operator(h: numpy.ndarray):
-        hin_m, hin_n = [hi.reshape(shape) for hi in numpy.split(h, 2)]
-        h_xyz = (m * hin_m +
-                 n * hin_n)
-        return [ifftn(hi) for hi in numpy.rollaxis(h_xyz, 3)]
-
-    return operator
-
-
-def inverse_maxwell_operator_approx(k0: numpy.ndarray,
-                                    G_matrix: numpy.ndarray,
-                                    epsilon: field_t,
-                                    mu: field_t = None
-                                    ) -> Callable[[numpy.ndarray], numpy.ndarray]:
-    """
-    Generate an approximate inverse of the Maxwell operator,
-        ik x conv(eps_k, ik x conv(mu_k, ___))
-     which can be used to improve the speed of ARPACK in shift-invert mode.
-
-    See the module-level docstring for more information.
-
-    :param k0: Bloch wavevector, [k0x, k0y, k0z].
-    :param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
-    :param epsilon: Dielectric constant distribution for the simulation.
-        All fields are sampled at cell centers (i.e., NOT Yee-gridded)
-    :param mu: Magnetic permability distribution for the simulation.
-        Default None (1 everywhere).
-    :return: Function which applies the approximate inverse of the maxwell operator to h_mn.
-    """
-    shape = epsilon[0].shape + (1,)
-    epsilon = numpy.stack(epsilon, 3)
-
-    k_mag, m, n = generate_kmn(k0, G_matrix, shape)
-
-    if mu is not None:
-        mu = numpy.stack(mu, 3)
-
-    def operator(h: numpy.ndarray):
-        """
-        Approximate inverse Maxwell operator for Bloch eigenmode simulation.
-
-        h is complex 2-field in (m, n) basis, vectorized
-
-        :param h: Raveled h_mn; size (2 * epsilon[0].size).
-        :return: Raveled ik x conv(eps_k, ik x conv(mu_k, h_mn))
-        """
-        hin_m, hin_n = [hi.reshape(shape) for hi in numpy.split(h, 2)]
-
-        #{d,e,h}_xyz fields are complex 3-fields in (1/x, 1/y, 1/z) basis
-
-        if mu is None:
-            b_m, b_n = hin_m, hin_n
-        else:
-            # transform from mn to xyz
-            h_xyz = (m * hin_m[:, :, :, None] +
-                     n * hin_n[:, :, :, None])
-
-            # multiply by mu
-            b_xyz = fftn(ifftn(h_xyz, axes=range(3)) * mu, axes=range(3))
-
-            # transform back to mn
-            b_m = numpy.sum(b_xyz * m, axis=3)
-            b_n = numpy.sum(b_xyz * n, axis=3)
-
-        # cross product and transform into xyz basis
-        e_xyz = (n * b_m -
-                 m * b_n) / k_mag
-
-        # multiply by epsilon
-        d_xyz = fftn(ifftn(e_xyz, axes=range(3)) * epsilon, axes=range(3))
-
-        # cross product and transform into mn basis   crossinv_t2c
-        h_m = numpy.sum(e_xyz * n, axis=3)[:, :, :, None] / +k_mag
-        h_n = numpy.sum(e_xyz * m, axis=3)[:, :, :, None] / -k_mag
-
-        return numpy.hstack((h_m.ravel(), h_n.ravel()))
-
-    return operator
-
-
-def eigsolve(num_modes: int,
-             k0: numpy.ndarray,
-             G_matrix: numpy.ndarray,
-             epsilon: field_t,
-             mu: field_t = None,
-             tolerance = 1e-8,
-             ) -> Tuple[numpy.ndarray, numpy.ndarray]:
-    """
-    Find the first (lowest-frequency) num_modes eigenmodes with Bloch wavevector
-     k0 of the specified structure.
-
-    :param k0: Bloch wavevector, [k0x, k0y, k0z].
-    :param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
-    :param epsilon: Dielectric constant distribution for the simulation.
-        All fields are sampled at cell centers (i.e., NOT Yee-gridded)
-    :param mu: Magnetic permability distribution for the simulation.
-        Default None (1 everywhere).
-    :return: (eigenvalues, eigenvectors) where eigenvalues[i] corresponds to the
-        vector eigenvectors[i, :]
-    """
-    h_size = 2 * epsilon[0].size
-
-    kmag = norm(G_matrix @ k0)
-
-    '''
-    Generate the operators
-    '''
-    mop = maxwell_operator(k0=k0, G_matrix=G_matrix, epsilon=epsilon, mu=mu)
-    imop = inverse_maxwell_operator_approx(k0=k0, G_matrix=G_matrix, epsilon=epsilon, mu=mu)
-
-    scipy_op = spalg.LinearOperator(dtype=complex, shape=(h_size, h_size), matvec=mop)
-    scipy_iop = spalg.LinearOperator(dtype=complex, shape=(h_size, h_size), matvec=imop)
-
-    y_shape = (h_size, num_modes)
-
-    def rayleigh_quotient(Z: numpy.ndarray, approx_grad: bool = True):
-        """
-        Absolute value of the block Rayleigh quotient, and the associated gradient.
-
-        See Johnson and Joannopoulos, Opt. Expr. 8, 3 (2001) for details (full
-         citation in module docstring).
-
-        ===
-
-        Notes on my understanding of the procedure:
-
-        Minimize f(Y) = |trace((Y.H @ A @ Y)|, making use of Y = Z @ inv(Z.H @ Z)^(1/2)
-         (a polar orthogonalization of Y). This gives f(Z) = |trace(Z.H @ A @ Z @ U)|,
-         where U = inv(Z.H @ Z). We minimize the absolute value to find the eigenvalues
-         with smallest magnitude.
-
-        The gradient is P @ (A @ Z @ U), where P = (1 - Z @ U @ Z.H) is a projection
-         onto the space orthonormal to Z. If approx_grad is True, the approximate
-         inverse of the maxwell operator is used to precondition the gradient.
-        """
-        z = Z.view(dtype=complex).reshape(y_shape)
-        U = numpy.linalg.inv(z.conj().T @ z)
-        zU = z @ U
-        AzU = scipy_op @ zU
-        zTAzU = z.conj().T @ AzU
-        f = numpy.real(numpy.trace(zTAzU))
-        if approx_grad:
-            df_dy = scipy_iop @ (AzU - zU @ zTAzU)
-        else:
-            df_dy = (AzU - zU @ zTAzU)
-        
-        df_dy_flat = df_dy.view(dtype=float).ravel()
-        return numpy.abs(f), numpy.sign(f) * df_dy_flat
-
-    '''
-    Use the conjugate gradient method and the approximate gradient calculation to
-     quickly find approximate eigenvectors.
-    '''
-    result = scipy.optimize.minimize(rayleigh_quotient,
-                                     numpy.random.rand(*y_shape, 2),
-                                     jac=True,
-                                     method='L-BFGS-B',
-                                     tol=1e-20,
-                                     options={'maxiter': 2000, 'gtol':0, 'ftol':1e-20 , 'disp':True})#, 'maxls':80, 'm':30})
-
-
-    result = scipy.optimize.minimize(lambda y: rayleigh_quotient(y, True),
-                                     result.x,
-                                     jac=True,
-                                     method='L-BFGS-B',
-                                     tol=1e-20,
-                                     options={'maxiter': 2000, 'gtol':0, 'disp':True})
-
-    result = scipy.optimize.minimize(lambda y: rayleigh_quotient(y, False),
-                                     result.x,
-                                     jac=True,
-                                     method='L-BFGS-B',
-                                     tol=1e-20,
-                                     options={'maxiter': 2000, 'gtol':0, 'disp':True})
-
-    for i in range(20):
-        result = scipy.optimize.minimize(lambda y: rayleigh_quotient(y, False),
-                                     result.x,
-                                     jac=True,
-                                     method='L-BFGS-B',
-                                     tol=1e-20,
-                                     options={'maxiter': 70, 'gtol':0, 'disp':True})
-        if result.nit == 0:
-            # We took 0 steps, so re-running won't help
-            break
-
-
-    z = result.x.view(dtype=complex).reshape(y_shape)
-
-    '''
-    Recover eigenvectors from Z
-    '''
-    U = numpy.linalg.inv(z.conj().T @ z)
-    y = z @ scipy.linalg.sqrtm(U)
-    w = y.conj().T @ (scipy_op @ y)
-
-    eigvals, w_eigvecs = numpy.linalg.eig(w)
-    eigvecs = y @ w_eigvecs
-
-    for i in range(len(eigvals)):
-        v = eigvecs[:, i]
-        n = eigvals[i]
-        v /= norm(v)
-        eigness = norm(scipy_op @ v - (v.conj() @ (scipy_op @ v)) * v )
-        f = numpy.sqrt(-numpy.real(n))
-        df = numpy.sqrt(-numpy.real(n + eigness))
-        neff_err = kmag * (1/df - 1/f)
-        logger.info('eigness {}: {}\n neff_err: {}'.format(i, eigness, neff_err))
-
-    order = numpy.argsort(numpy.abs(eigvals))
-    return eigvals[order], eigvecs.T[order]
-
-
-def find_k(frequency: float,
-           tolerance: float,
-           direction: numpy.ndarray,
-           G_matrix: numpy.ndarray,
-           epsilon: field_t,
-           mu: field_t = None,
-           band: int = 0,
-           k_min: float = 0,
-           k_max: float = 0.5,
-           ) -> Tuple[numpy.ndarray, float]:
-    """
-    Search for a bloch vector that has a given frequency.
-
-    :param frequency: Target frequency.
-    :param tolerance: Target frequency tolerance.
-    :param direction: k-vector direction to search along.
-    :param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
-    :param epsilon: Dielectric constant distribution for the simulation.
-        All fields are sampled at cell centers (i.e., NOT Yee-gridded)
-    :param mu: Magnetic permability distribution for the simulation.
-        Default None (1 everywhere).
-    :param band: Which band to search in. Default 0 (lowest frequency).
-    return: (k, actual_frequency) The found k-vector and its frequency
-    """
-
-    direction = numpy.array(direction) / norm(direction)
-
-    def get_f(k0_mag: float, band: int = 0):
-        k0 = direction * k0_mag
-        n, _v = eigsolve(band + 1, k0, G_matrix=G_matrix, epsilon=epsilon)
-        f = numpy.sqrt(numpy.abs(numpy.real(n[band])))
-        return f
-
-    res = scipy.optimize.minimize_scalar(lambda x: abs(get_f(x, band) - frequency),
-                                         (k_min + k_max) / 2,
-                                         method='Bounded',
-                                         bounds=(k_min, k_max),
-                                         options={'xatol': abs(tolerance)})
-    return res.x * direction, res.fun + frequency
-
-
diff --git a/fdfd_tools/eigensolvers.py b/fdfd_tools/eigensolvers.py
deleted file mode 100644
index d0ee541..0000000
--- a/fdfd_tools/eigensolvers.py
+++ /dev/null
@@ -1,120 +0,0 @@
-"""
-Solvers for eigenvalue / eigenvector problems
-"""
-from typing import Tuple, List
-import numpy
-from numpy.linalg import norm
-from scipy import sparse
-import scipy.sparse.linalg as spalg
-
-
-def power_iteration(operator: sparse.spmatrix,
-                    guess_vector: numpy.ndarray = None,
-                    iterations: int = 20,
-                    ) -> Tuple[complex, numpy.ndarray]:
-    """
-    Use power iteration to estimate the dominant eigenvector of a matrix.
-
-    :param operator: Matrix to analyze.
-    :param guess_vector: Starting point for the eigenvector. Default is a randomly chosen vector.
-    :param iterations: Number of iterations to perform. Default 20.
-    :return: (Largest-magnitude eigenvalue, Corresponding eigenvector estimate)
-    """
-    if numpy.any(numpy.equal(guess_vector, None)):
-        v = numpy.random.rand(operator.shape[0])
-    else:
-        v = guess_vector
-
-    for _ in range(iterations):
-        v = operator @ v
-        v /= norm(v)
-
-    lm_eigval = v.conj() @ (operator @ v)
-    return lm_eigval, v
-
-
-def rayleigh_quotient_iteration(operator: sparse.spmatrix or spalg.LinearOperator,
-                                guess_vector: numpy.ndarray,
-                                iterations: int = 40,
-                                tolerance: float = 1e-13,
-                                solver=None,
-                                ) -> Tuple[complex, numpy.ndarray]:
-    """
-    Use Rayleigh quotient iteration to refine an eigenvector guess.
-
-    :param operator: Matrix to analyze.
-    :param guess_vector: Eigenvector to refine.
-    :param iterations: Maximum number of iterations to perform. Default 40.
-    :param tolerance: Stop iteration if (A - I*eigenvalue) @ v < tolerance.
-        Default 1e-13.
-    :param solver: Solver function of the form x = solver(A, b).
-        By default, use scipy.sparse.spsolve for sparse matrices and
-        scipy.sparse.bicgstab for general LinearOperator instances.
-    :return: (eigenvalue, eigenvector)
-    """
-    try:
-        _test = operator - sparse.eye(operator.shape[0])
-        shift = lambda eigval: eigval * sparse.eye(operator.shape[0])
-        if solver is None:
-            solver = spalg.spsolve
-    except TypeError:
-        shift = lambda eigval: spalg.LinearOperator(shape=operator.shape,
-                                                    dtype=operator.dtype,
-                                                    matvec=lambda v: eigval * v)
-        if solver is None:
-            solver = lambda A, b: spalg.bicgstab(A, b)[0]
-
-    v = guess_vector
-    v /= norm(v)
-    for _ in range(iterations):
-        eigval = v.conj() @ (operator @ v)
-        if norm(operator @ v - eigval * v) < tolerance:
-            break
-
-        shifted_operator = operator - shift(eigval)
-        v = solver(shifted_operator, v)
-        v /= norm(v)
-    return eigval, v
-
-
-def signed_eigensolve(operator: sparse.spmatrix or spalg.LinearOperator,
-                      how_many: int,
-                      negative: bool = False,
-                      ) -> Tuple[numpy.ndarray, numpy.ndarray]:
-    """
-    Find the largest-magnitude positive-only (or negative-only) eigenvalues and
-     eigenvectors of the provided matrix.
-
-    :param operator: Matrix to analyze.
-    :param how_many: How many eigenvalues to find.
-    :param negative: Whether to find negative-only eigenvalues.
-        Default False (positive only).
-    :return: (sorted list of eigenvalues, 2D ndarray of corresponding eigenvectors)
-        eigenvectors[:, k] corresponds to the k-th eigenvalue
-    """
-    # Use power iteration to estimate the dominant eigenvector
-    lm_eigval, _ = power_iteration(operator)
-
-    '''
-    Shift by the absolute value of the largest eigenvalue, then find a few of the
-     largest-magnitude (shifted) eigenvalues. A positive shift ensures that we find the
-     largest _positive_ eigenvalues, since any negative eigenvalues will be shifted to the
-     range 0 >= neg_eigval + abs(lm_eigval) > abs(lm_eigval)
-    '''
-    shift = numpy.abs(lm_eigval)
-    if negative:
-        shift *= -1
-
-    # Try to combine, use general LinearOperator if we fail
-    try:
-        shifted_operator = operator + shift * sparse.eye(operator.shape[0])
-    except TypeError:
-        shifted_operator = operator + spalg.LinearOperator(shape=operator.shape,
-                                                           matvec=lambda v: shift * v)
-
-    shifted_eigenvalues, eigenvectors = spalg.eigs(shifted_operator, which='LM', k=how_many, ncv=50)
-    eigenvalues = shifted_eigenvalues - shift
-
-    k = eigenvalues.argsort()
-    return eigenvalues[k], eigenvectors[:, k]
-
diff --git a/fdfd_tools/farfield.py b/fdfd_tools/farfield.py
deleted file mode 100644
index 84a04ba..0000000
--- a/fdfd_tools/farfield.py
+++ /dev/null
@@ -1,220 +0,0 @@
-"""
-Functions for performing near-to-farfield transformation (and the reverse).
-"""
-from typing import Dict, List
-import numpy
-from numpy.fft import fft2, fftshift, fftfreq, ifft2, ifftshift
-from numpy import pi
-
-
-def near_to_farfield(E_near: List[numpy.ndarray],
-                     H_near: List[numpy.ndarray],
-                     dx: float,
-                     dy: float,
-                     padded_size: List[int] = None
-                     ) -> Dict[str]:
-    """
-    Compute the farfield, i.e. the distribution of the fields after propagation
-      through several wavelengths of uniform medium.
-
-    The input fields should be complex phasors.
-
-    :param E_near: List of 2 ndarrays containing the 2D phasor field slices for the transverse
-        E fields (e.g. [Ex, Ey] for calculating the farfield toward the z-direction).
-    :param H_near: List of 2 ndarrays containing the 2D phasor field slices for the transverse
-        H fields (e.g. [Hx, hy] for calculating the farfield towrad the z-direction).
-    :param dx: Cell size along x-dimension, in units of wavelength.
-    :param dy: Cell size along y-dimension, in units of wavelength.
-    :param padded_size: Shape of the output. A single integer `n` will be expanded to `(n, n)`.
-        Powers of 2 are most efficient for FFT computation.
-        Default is the smallest power of 2 larger than the input, for each axis.
-    :returns: Dict with keys
-        'E_far': Normalized E-field farfield; multiply by
-            (i k exp(-i k r) / (4 pi r)) to get the actual field value.
-        'H_far': Normalized H-field farfield; multiply by
-            (i k exp(-i k r) / (4 pi r)) to get the actual field value.
-        'kx', 'ky': Wavevector values corresponding to the x- and y- axes in E_far and H_far,
-            normalized to wavelength (dimensionless).
-        'dkx', 'dky': step size for kx and ky, normalized to wavelength.
-        'theta': arctan2(ky, kx) corresponding to each (kx, ky).
-            This is the angle in the x-y plane, counterclockwise from above, starting from +x.
-        'phi': arccos(kz / k) corresponding to each (kx, ky).
-            This is the angle away from +z.
-    """
-
-    if not len(E_near) == 2:
-        raise Exception('E_near must be a length-2 list of ndarrays')
-    if not len(H_near) == 2:
-        raise Exception('H_near must be a length-2 list of ndarrays')
-
-    s = E_near[0].shape
-    if not all(s == f.shape for f in E_near + H_near):
-        raise Exception('All fields must be the same shape!')
-
-    if padded_size is None:
-        padded_size = (2**numpy.ceil(numpy.log2(s))).astype(int)
-    if not hasattr(padded_size, '__len__'):
-        padded_size = (padded_size, padded_size)
-
-    En_fft = [fftshift(fft2(fftshift(Eni), s=padded_size)) for Eni in E_near]
-    Hn_fft = [fftshift(fft2(fftshift(Hni), s=padded_size)) for Hni in H_near]
-
-    # Propagation vectors kx, ky
-    k  = 2 * pi
-    kxx = 2 * pi * fftshift(fftfreq(padded_size[0], dx))
-    kyy = 2 * pi * fftshift(fftfreq(padded_size[1], dy))
-
-    kx, ky = numpy.meshgrid(kxx, kyy, indexing='ij')
-    kxy2 = kx * kx + ky * ky
-    kxy = numpy.sqrt(kxy2)
-    kz = numpy.sqrt(k * k - kxy2)
-
-    sin_th = ky / kxy
-    cos_th = kx / kxy
-    cos_phi = kz / k
-
-    sin_th[numpy.logical_and(kx == 0, ky == 0)] = 0
-    cos_th[numpy.logical_and(kx == 0, ky == 0)] = 1
-
-    # Normalized vector potentials N, L
-    N = [-Hn_fft[1] * cos_phi * cos_th + Hn_fft[0] * cos_phi * sin_th,
-          Hn_fft[1] * sin_th + Hn_fft[0] * cos_th]
-    L = [ En_fft[1] * cos_phi * cos_th - En_fft[0] * cos_phi * sin_th,
-         -En_fft[1] * sin_th - En_fft[0] * cos_th]
-
-    E_far = [-L[1] - N[0],
-              L[0] - N[1]]
-    H_far = [-E_far[1],
-              E_far[0]]
-
-    theta = numpy.arctan2(ky, kx)
-    phi = numpy.arccos(cos_phi)
-    theta[numpy.logical_and(kx == 0, ky == 0)] = 0
-    phi[numpy.logical_and(kx == 0, ky == 0)] = 0
-
-    # Zero fields beyond valid (phi, theta)
-    invalid_ind = kxy2 >= k * k
-    theta[invalid_ind] = 0
-    phi[invalid_ind] = 0
-    for i in range(2):
-        E_far[i][invalid_ind] = 0
-        H_far[i][invalid_ind] = 0
-
-    outputs = {
-        'E': E_far,
-        'H': H_far,
-        'dkx': kx[1]-kx[0],
-        'dky': ky[1]-ky[0],
-        'kx': kx,
-        'ky': ky,
-        'theta': theta,
-        'phi': phi,
-    }
-
-    return outputs
-
-
-
-def far_to_nearfield(E_far: List[numpy.ndarray],
-                     H_far: List[numpy.ndarray],
-                     dkx: float,
-                     dky: float,
-                     padded_size: List[int] = None
-                     ) -> Dict[str]:
-    """
-    Compute the farfield, i.e. the distribution of the fields after propagation
-      through several wavelengths of uniform medium.
-
-    The input fields should be complex phasors.
-
-    :param E_far: List of 2 ndarrays containing the 2D phasor field slices for the transverse
-        E fields (e.g. [Ex, Ey] for calculating the nearfield toward the z-direction).
-        Fields should be normalized so that
-            E_far = E_far_actual / (i k exp(-i k r) / (4 pi r))
-    :param H_far: List of 2 ndarrays containing the 2D phasor field slices for the transverse
-        H fields (e.g. [Hx, hy] for calculating the nearfield toward the z-direction).
-        Fields should be normalized so that
-            H_far = H_far_actual / (i k exp(-i k r) / (4 pi r))
-    :param dkx: kx discretization, in units of wavelength.
-    :param dky: ky discretization, in units of wavelength.
-    :param padded_size: Shape of the output. A single integer `n` will be expanded to `(n, n)`.
-        Powers of 2 are most efficient for FFT computation.
-        Default is the smallest power of 2 larger than the input, for each axis.
-    :returns: Dict with keys
-        'E': E-field nearfield
-        'H': H-field nearfield
-        'dx', 'dy': spatial discretization, normalized to wavelength (dimensionless)
-    """
-
-    if not len(E_far) == 2:
-        raise Exception('E_far must be a length-2 list of ndarrays')
-    if not len(H_far) == 2:
-        raise Exception('H_far must be a length-2 list of ndarrays')
-
-    s = E_far[0].shape
-    if not all(s == f.shape for f in E_far + H_far):
-        raise Exception('All fields must be the same shape!')
-
-    if padded_size is None:
-        padded_size = (2**numpy.ceil(numpy.log2(s))).astype(int)
-    if not hasattr(padded_size, '__len__'):
-        padded_size = (padded_size, padded_size)
-
-
-    k = 2 * pi
-    kxs = fftshift(fftfreq(s[0], 1/(s[0] * dkx)))
-    kys = fftshift(fftfreq(s[0], 1/(s[1] * dky)))
-
-    kx, ky = numpy.meshgrid(kxs, kys, indexing='ij')
-    kxy2 = kx * kx + ky * ky
-    kxy = numpy.sqrt(kxy2)
-
-    kz = numpy.sqrt(k * k - kxy2)
-
-    sin_th = ky / kxy
-    cos_th = kx / kxy
-    cos_phi = kz / k
-
-    sin_th[numpy.logical_and(kx == 0, ky == 0)] = 0
-    cos_th[numpy.logical_and(kx == 0, ky == 0)] = 1
-
-    # Zero fields beyond valid (phi, theta)
-    invalid_ind = kxy2 >= k * k
-    theta[invalid_ind] = 0
-    phi[invalid_ind] = 0
-    for i in range(2):
-        E_far[i][invalid_ind] = 0
-        H_far[i][invalid_ind] = 0
-
-
-    # Normalized vector potentials N, L
-    L = [0.5 * E_far[1],
-        -0.5 * E_far[0]]
-    N = [L[1],
-        -L[0]]
-
-    En_fft = [-( L[0] * sin_th + L[1] * cos_phi * cos_th)/cos_phi,
-              -(-L[0] * cos_th + L[1] * cos_phi * sin_th)/cos_phi]
-
-    Hn_fft = [( N[0] * sin_th + N[1] * cos_phi * cos_th)/cos_phi,
-              (-N[0] * cos_th + N[1] * cos_phi * sin_th)/cos_phi]
-
-    for i in range(2):
-        En_fft[i][cos_phi == 0] = 0
-        Hn_fft[i][cos_phi == 0] = 0
-
-    E_near = [ifftshift(ifft2(ifftshift(Ei), s=padded_size)) for Ei in En_fft]
-    H_near = [ifftshift(ifft2(ifftshift(Hi), s=padded_size)) for Hi in Hn_fft]
-
-    dx = 2 * pi / (s[0] * dkx)
-    dy = 2 * pi / (s[0] * dky)
-
-    outputs = {
-        'E': E_near,
-        'H': H_near,
-        'dx': dx,
-        'dy': dy,
-    }
-
-    return outputs
-
diff --git a/fdfd_tools/operators.py b/fdfd_tools/operators.py
index 02c1197..1691f36 100644
--- a/fdfd_tools/operators.py
+++ b/fdfd_tools/operators.py
@@ -182,10 +182,10 @@ def eh_full(omega: complex,
     :param mu: Vectorized magnetic permeability (default 1 everywhere)
     :param pec: Vectorized mask specifying PEC cells. Any cells where pec != 0 are interpreted
         as containing a perfect electrical conductor (PEC).
-        The PEC is applied per-field-component (i.e., pec.size == epsilon.size)
+        The PEC is applied per-field-component (ie, pec.size == epsilon.size)
     :param pmc: Vectorized mask specifying PMC cells. Any cells where pmc != 0 are interpreted
         as containing a perfect magnetic conductor (PMC).
-        The PMC is applied per-field-component (i.e., pmc.size == epsilon.size)
+        The PMC is applied per-field-component (ie, pmc.size == epsilon.size)
     :return: Sparse matrix containing the wave operator
     """
     if numpy.any(numpy.equal(pec, None)):
@@ -284,8 +284,7 @@ def m2j(omega: complex,
 
 def rotation(axis: int, shape: List[int], shift_distance: int=1) -> sparse.spmatrix:
     """
-    Utility operator for performing a circular shift along a specified axis by a
-     specified number of elements.
+    Utility operator for performing a circular shift along a specified axis by 1 element.
 
     :param axis: Axis to shift along. x=0, y=1, z=2
     :param shape: Shape of the grid being shifted
@@ -305,7 +304,7 @@ def rotation(axis: int, shape: List[int], shift_distance: int=1) -> sparse.spmat
     i_ind = numpy.arange(n)
     j_ind = numpy.ravel_multi_index(ijk, shape, order='C')
 
-    vij = (numpy.ones(n), (i_ind, j_ind.ravel(order='C')))
+    vij = (numpy.ones(n), (i_ind, j_ind.flatten(order='C')))
 
     d = sparse.csr_matrix(vij, shape=(n, n))
 
@@ -349,7 +348,7 @@ def shift_with_mirror(axis: int, shape: List[int], shift_distance: int=1) -> spa
     if len(shape) == 3:
         j_ind += ijk[2] * shape[0] * shape[1]
 
-    vij = (numpy.ones(n), (i_ind, j_ind.ravel(order='C')))
+    vij = (numpy.ones(n), (i_ind, j_ind.flatten(order='C')))
 
     d = sparse.csr_matrix(vij, shape=(n, n))
     return d
@@ -370,7 +369,7 @@ def deriv_forward(dx_e: List[numpy.ndarray]) -> List[sparse.spmatrix]:
     def deriv(axis):
         return rotation(axis, shape, 1) - sparse.eye(n)
 
-    Ds = [sparse.diags(+1 / dx.ravel(order='C')) @ deriv(a)
+    Ds = [sparse.diags(+1 / dx.flatten(order='C')) @ deriv(a)
           for a, dx in enumerate(dx_e_expanded)]
 
     return Ds
@@ -391,7 +390,7 @@ def deriv_back(dx_h: List[numpy.ndarray]) -> List[sparse.spmatrix]:
     def deriv(axis):
         return rotation(axis, shape, -1) - sparse.eye(n)
 
-    Ds = [sparse.diags(-1 / dx.ravel(order='C')) @ deriv(a)
+    Ds = [sparse.diags(-1 / dx.flatten(order='C')) @ deriv(a)
           for a, dx in enumerate(dx_h_expanded)]
 
     return Ds
@@ -462,8 +461,8 @@ def poynting_e_cross(e: vfield_t, dxes: dx_lists_t) -> sparse.spmatrix:
     fx, fy, fz = [avgf(i, shape) for i in range(3)]
     bx, by, bz = [avgb(i, shape) for i in range(3)]
 
-    dxag = [dx.ravel(order='C') for dx in numpy.meshgrid(*dxes[0], indexing='ij')]
-    dbgx, dbgy, dbgz = [sparse.diags(dx.ravel(order='C'))
+    dxag = [dx.flatten(order='C') for dx in numpy.meshgrid(*dxes[0], indexing='ij')]
+    dbgx, dbgy, dbgz = [sparse.diags(dx.flatten(order='C'))
                         for dx in numpy.meshgrid(*dxes[1], indexing='ij')]
 
     Ex, Ey, Ez = [sparse.diags(ei * da) for ei, da in zip(numpy.split(e, 3), dxag)]
@@ -491,8 +490,8 @@ def poynting_h_cross(h: vfield_t, dxes: dx_lists_t) -> sparse.spmatrix:
     fx, fy, fz = [avgf(i, shape) for i in range(3)]
     bx, by, bz = [avgb(i, shape) for i in range(3)]
 
-    dxbg = [dx.ravel(order='C') for dx in numpy.meshgrid(*dxes[1], indexing='ij')]
-    dagx, dagy, dagz = [sparse.diags(dx.ravel(order='C'))
+    dxbg = [dx.flatten(order='C') for dx in numpy.meshgrid(*dxes[1], indexing='ij')]
+    dagx, dagy, dagz = [sparse.diags(dx.flatten(order='C'))
                         for dx in numpy.meshgrid(*dxes[0], indexing='ij')]
 
     Hx, Hy, Hz = [sparse.diags(hi * db) for hi, db in zip(numpy.split(h, 3), dxbg)]
diff --git a/fdfd_tools/solvers.py b/fdfd_tools/solvers.py
index 066725c..bb230f2 100644
--- a/fdfd_tools/solvers.py
+++ b/fdfd_tools/solvers.py
@@ -3,7 +3,6 @@ Solvers for FDFD problems.
 """
 
 from typing import List, Callable, Dict, Any
-import logging
 
 import numpy
 from numpy.linalg import norm
@@ -12,9 +11,6 @@ import scipy.sparse.linalg
 from . import operators
 
 
-logger = logging.getLogger(__name__)
-
-
 def _scipy_qmr(A: scipy.sparse.csr_matrix,
                b: numpy.ndarray,
                **kwargs
@@ -33,20 +29,20 @@ def _scipy_qmr(A: scipy.sparse.csr_matrix,
     '''
     iter = 0
 
-    def log_residual(xk):
+    def print_residual(xk):
         nonlocal iter
         iter += 1
         if iter % 100 == 0:
-            logger.info('Solver residual at iteration {} : {}'.format(iter, norm(A @ xk - b)))
+            print('Solver residual at iteration', iter, ':', norm(A @ xk - b))
 
     if 'callback' in kwargs:
         def augmented_callback(xk):
-            log_residual(xk)
+            print_residual(xk)
             kwargs['callback'](xk)
 
         kwargs['callback'] = augmented_callback
     else:
-        kwargs['callback'] = log_residual
+        kwargs['callback'] = print_residual
 
     '''
     Run the actual solve
@@ -87,7 +83,7 @@ def generic(omega: complex,
               b: numpy.ndarray
               x: numpy.ndarray
         Default is a wrapped version of scipy.sparse.linalg.qmr()
-         which doesn't return convergence info and logs the residual
+         which doesn't return convergence info and prints the residual
          every 100 iterations.
     :param matrix_solver_opts: Passed as kwargs to matrix_solver(...)
     :return: E-field which solves the system.
diff --git a/fdfd_tools/vectorization.py b/fdfd_tools/vectorization.py
index e7b9645..2377d39 100644
--- a/fdfd_tools/vectorization.py
+++ b/fdfd_tools/vectorization.py
@@ -1,7 +1,7 @@
 """
 Functions for moving between a vector field (list of 3 ndarrays, [f_x, f_y, f_z])
 and a 1D array representation of that field [f_x0, f_x1, f_x2,... f_y0,... f_z0,...].
-Vectorized versions of the field use row-major (ie., C-style) ordering.
+Vectorized versions of the field use column-major (ie., Fortran, Matlab) ordering.
 """
 
 
@@ -27,7 +27,7 @@ def vec(f: field_t) -> vfield_t:
     """
     if numpy.any(numpy.equal(f, None)):
         return None
-    return numpy.hstack(tuple((fi.ravel(order='C') for fi in f)))
+    return numpy.hstack(tuple((fi.flatten(order='C') for fi in f)))
 
 
 def unvec(v: vfield_t, shape: numpy.ndarray) -> field_t:
diff --git a/fdfd_tools/waveguide.py b/fdfd_tools/waveguide.py
index 1566a74..a8ae1f2 100644
--- a/fdfd_tools/waveguide.py
+++ b/fdfd_tools/waveguide.py
@@ -23,7 +23,7 @@ import numpy
 from numpy.linalg import norm
 import scipy.sparse as sparse
 
-from . import vec, unvec, dx_lists_t, field_t, vfield_t
+from . import unvec, dx_lists_t, field_t, vfield_t
 from . import operators
 
 
@@ -307,62 +307,3 @@ def e_err(e: vfield_t,
         op = ch @ mu_inv @ ce @ e - omega ** 2 * (epsilon * e)
 
     return norm(op) / norm(e)
-
-
-def cylindrical_operator(omega: complex,
-                         dxes: dx_lists_t,
-                         epsilon: vfield_t,
-                         r0: float,
-                         ) -> sparse.spmatrix:
-    """
-    Cylindrical coordinate waveguide operator of the form
-
-    TODO
-
-    for use with a field vector of the form [E_r, E_y].
-
-    This operator can be used to form an eigenvalue problem of the form
-    A @ [E_r, E_y] = wavenumber**2 * [E_r, E_y]
-
-    which can then be solved for the eigenmodes of the system (an exp(-i * wavenumber * theta)
-     theta-dependence is assumed for the fields).
-
-    :param omega: The angular frequency of the system
-    :param dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header (2D)
-    :param epsilon: Vectorized dielectric constant grid
-    :param r0: Radius of curvature for the simulation. This should be the minimum value of
-        r within the simulation domain.
-    :return: Sparse matrix representation of the operator
-    """
-
-    Dfx, Dfy = operators.deriv_forward(dxes[0])
-    Dbx, Dby = operators.deriv_back(dxes[1])
-
-    rx = r0 + numpy.cumsum(dxes[0][0])
-    ry = r0 + dxes[0][0]/2.0 + numpy.cumsum(dxes[1][0])
-    tx = rx/r0
-    ty = ry/r0
-
-    Tx = sparse.diags(vec(tx[:, None].repeat(dxes[0][1].size, axis=1)))
-    Ty = sparse.diags(vec(ty[:, None].repeat(dxes[1][1].size, axis=1)))
-
-    eps_parts = numpy.split(epsilon, 3)
-    eps_x = sparse.diags(eps_parts[0])
-    eps_y = sparse.diags(eps_parts[1])
-    eps_z_inv = sparse.diags(1 / eps_parts[2])
-
-    pa = sparse.vstack((Dfx, Dfy)) @ Tx @ eps_z_inv @ sparse.hstack((Dbx, Dby))
-    pb = sparse.vstack((Dfx, Dfy)) @ Tx @ eps_z_inv @ sparse.hstack((Dby, Dbx))
-    a0 = Ty @ eps_x + omega**-2 * Dby @ Ty @ Dfy
-    a1 = Tx @ eps_y + omega**-2 * Dbx @ Ty @ Dfx
-    b0 = Dbx @ Ty @ Dfy
-    b1 = Dby @ Ty @ Dfx
-
-    diag = sparse.block_diag
-    op = (omega**2 * diag((Tx, Ty)) + pa) @ diag((a0, a1)) + \
-        - (sparse.bmat(((None, Ty), (Tx, None))) + omega**-2 * pb) @ diag((b0, b1))
-
-    return op
-
-
-
diff --git a/fdfd_tools/waveguide_mode.py b/fdfd_tools/waveguide_mode.py
index e50cc6a..df62143 100644
--- a/fdfd_tools/waveguide_mode.py
+++ b/fdfd_tools/waveguide_mode.py
@@ -1,10 +1,10 @@
 from typing import Dict, List
 import numpy
 import scipy.sparse as sparse
+import scipy.sparse.linalg as spalg
 
 from . import vec, unvec, dx_lists_t, vfield_t, field_t
 from . import operators, waveguide, functional
-from .eigensolvers import signed_eigensolve, rayleigh_quotient_iteration
 
 
 def solve_waveguide_mode_2d(mode_number: int,
@@ -12,12 +12,12 @@ def solve_waveguide_mode_2d(mode_number: int,
                             dxes: dx_lists_t,
                             epsilon: vfield_t,
                             mu: vfield_t = None,
-                            wavenumber_correction: bool = True,
+                            wavenumber_correction: bool = True
                             ) -> Dict[str, complex or field_t]:
     """
     Given a 2d region, attempts to solve for the eigenmode with the specified mode number.
 
-    :param mode_number: Number of the mode, 0-indexed.
+    :param mode_number: Number of the mode, 0-indexed
     :param omega: Angular frequency of the simulation
     :param dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header
     :param epsilon: Dielectric constant
@@ -29,19 +29,46 @@ def solve_waveguide_mode_2d(mode_number: int,
 
     '''
     Solve for the largest-magnitude eigenvalue of the real operator
+     by using power iteration.
     '''
     dxes_real = [[numpy.real(dx) for dx in dxi] for dxi in dxes]
+
     A_r = waveguide.operator(numpy.real(omega), dxes_real, numpy.real(epsilon), numpy.real(mu))
 
-    eigvals, eigvecs = signed_eigensolve(A_r, mode_number+3)
-    v = eigvecs[:, -(mode_number + 1)]
+    # Use power iteration for 20 steps to estimate the dominant eigenvector
+    v = numpy.random.rand(A_r.shape[0])
+    for _ in range(20):
+        v = A_r @ v
+        v /= numpy.linalg.norm(v)
+
+    lm_eigval = v @ A_r @ v
+
+    '''
+    Shift by the absolute value of the largest eigenvalue, then find a few of the
+     largest (shifted) eigenvalues. The shift ensures that we find the largest
+     _positive_ eigenvalues, since any negative eigenvalues will be shifted to the range
+     0 >= neg_eigval + abs(lm_eigval) > abs(lm_eigval)
+    '''
+    shifted_A_r = A_r + abs(lm_eigval) * sparse.eye(A_r.shape[0])
+    eigvals, eigvecs = spalg.eigs(shifted_A_r, which='LM', k=mode_number + 3, ncv=50)
+
+    # Pick the eigenvalue we want from the few we found
+    k = eigvals.argsort()[-(mode_number+1)]
+    v = eigvecs[:, k]
 
     '''
     Now solve for the eigenvector of the full operator, using the real operator's
      eigenvector as an initial guess for Rayleigh quotient iteration.
     '''
     A = waveguide.operator(omega, dxes, epsilon, mu)
-    eigval, v = rayleigh_quotient_iteration(A, v)
+
+    eigval = None
+    for _ in range(40):
+        eigval = v @ A @ v
+        if numpy.linalg.norm(A @ v - eigval * v) < 1e-13:
+            break
+        w = spalg.spsolve(A - eigval * sparse.eye(A.shape[0]), v)
+        v = w / numpy.linalg.norm(w)
 
     # Calculate the wave-vector (force the real part to be positive)
     wavenumber = numpy.sqrt(eigval)
@@ -272,69 +299,3 @@ def compute_overlap_e(E: field_t,
     overlap_e /= norm_factor * dx_forward
 
     return unvec(overlap_e, E[0].shape)
-
-
-def solve_waveguide_mode_cylindrical(mode_number: int,
-                                     omega: complex,
-                                     dxes: dx_lists_t,
-                                     epsilon: vfield_t,
-                                     r0: float,
-                                     wavenumber_correction: bool = True,
-                                     ) -> Dict[str, complex or field_t]:
-    """
-    Given a 2d (r, y) slice of epsilon, attempts to solve for the eigenmode
-     of the bent waveguide with the specified mode number.
-
-    :param mode_number: Number of the mode, 0-indexed
-    :param omega: Angular frequency of the simulation
-    :param dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header.
-        The first coordinate is assumed to be r, the second is y.
-    :param epsilon: Dielectric constant
-    :param r0: Radius of curvature for the simulation. This should be the minimum value of
-        r within the simulation domain.
-    :param wavenumber_correction: Whether to correct the wavenumber to
-        account for numerical dispersion (default True)
-    :return: {'E': List[numpy.ndarray], 'H': List[numpy.ndarray], 'wavenumber': complex}
-    """
-
-    '''
-    Solve for the largest-magnitude eigenvalue of the real operator
-    '''
-    dxes_real = [[numpy.real(dx) for dx in dxi] for dxi in dxes]
-
-    A_r = waveguide.cylindrical_operator(numpy.real(omega), dxes_real, numpy.real(epsilon), r0)
-    eigvals, eigvecs = signed_eigensolve(A_r, mode_number + 3)
-    v = eigvecs[:, -(mode_number+1)]
-
-    '''
-    Now solve for the eigenvector of the full operator, using the real operator's
-     eigenvector as an initial guess for Rayleigh quotient iteration.
-    '''
-    A = waveguide.cylindrical_operator(omega, dxes, epsilon, r0)
-    eigval, v = rayleigh_quotient_iteration(A, v)
-
-    # Calculate the wave-vector (force the real part to be positive)
-    wavenumber = numpy.sqrt(eigval)
-    wavenumber *= numpy.sign(numpy.real(wavenumber))
-
-    '''
-    Perform correction on wavenumber to account for numerical dispersion.
-
-     See Numerical Dispersion in Taflove's FDTD book.
-     This correction term reduces the error in emitted power, but additional
-      error is introduced into the E_err and H_err terms. This effect becomes
-      more pronounced as the wavenumber increases.
-    '''
-    if wavenumber_correction:
-        wavenumber -= 2 * numpy.sin(numpy.real(wavenumber / 2)) - numpy.real(wavenumber)
-
-    shape = [d.size for d in dxes[0]]
-    v = numpy.hstack((v, numpy.zeros(shape[0] * shape[1])))
-    fields = {
-        'wavenumber': wavenumber,
-        'E': unvec(v, shape),
-#        'E': unvec(e, shape),
-#        'H': unvec(h, shape),
-    }
-
-    return fields
diff --git a/setup.py b/setup.py
index ef1df08..4a3441a 100644
--- a/setup.py
+++ b/setup.py
@@ -3,7 +3,7 @@
 from setuptools import setup, find_packages
 
 setup(name='fdfd_tools',
-      version='0.4',
+      version='0.3',
       description='FDFD Electromagnetic simulation tools',
       author='Jan Petykiewicz',
       author_email='anewusername@gmail.com',