12 changed files with 54 additions and 1083 deletions
--- 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.
--- a/examples/bloch.py
+++ b/examples/bloch.py
@ -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 ))
--- 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()
--- a/fdfd_tools/bloch.py
+++ b/fdfd_tools/bloch.py
@ -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
--- a/fdfd_tools/eigensolvers.py
+++ b/fdfd_tools/eigensolvers.py
@ -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]
--- a/fdfd_tools/farfield.py
+++ b/fdfd_tools/farfield.py
@ -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
--- 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
+    Utility operator for performing a circular shift along a specified axis by 1 element.
     specified number of elements.
    :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')]
+    dxag = [dx.flatten(order='C') for dx in numpy.meshgrid(*dxes[0], indexing='ij')]
-    dbgx, dbgy, dbgz = [sparse.diags(dx.ravel(order='C'))
+    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')]
+    dxbg = [dx.flatten(order='C') for dx in numpy.meshgrid(*dxes[1], indexing='ij')]
-    dagx, dagy, dagz = [sparse.diags(dx.ravel(order='C'))
+    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)]
--- 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.
--- 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:
--- 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
--- 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)
+    # Use power iteration for 20 steps to estimate the dominant eigenvector
-    v = eigvecs[:, -(mode_number + 1)]
+    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
--- 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',