# Path: ceviche/fdfd.py
# class fdfd():
# class fdfd_ez(fdfd):
# class fdfd_hz(fdfd):
# class fdfd_mf_ez(fdfd):
# class fdfd_3d(fdfd):
#     def __init__(self, omega, dL, eps_r, npml, bloch_phases=None):
#     def eps_r(self):
#     def eps_r(self, new_eps):
#     def _make_A(self, eps_r):
#     def _solve_fn(self, entries_a, indices_a, source_vec):
#     def solve(self, source_z):
#     def _setup_derivatives(self):
#     def _setup_bloch_phases(self, bloch_phases):
#     def _vec_to_grid(self, vec):
#     def _grid_to_vec(self, grid):
#     def _save_shape(self, grid):
#     def _default_val(val, default_val=None):
#     def _Ex_Ey_to_Hz(self, Ex_vec, Ey_vec):
#     def _Ez_to_Hx(self, Ez_vec):
#     def _Ez_to_Hy(self, Ez_vec):
#     def _Ez_to_Hx_Hy(self, Ez_vec):
#     def _Hz_to_Ex(self, Hz_vec, eps_vec_xx):
#     def _Hz_to_Ey(self, Hz_vec, eps_vec_yy):
#     def _Hx_Hy_to_Ez(self, Hx_vec, Hy_vec, eps_vec_zz):
#     def _Hz_to_Ex_Ey(self, Hz_vec, eps_vec_xx, eps_vec_yy):
#     def __init__(self, omega, dL, eps_r, npml, bloch_phases=None):
#     def _make_A(self, eps_vec):
#     def _solve_fn(self, eps_vec, entries_a, indices_a, Jz_vec):
#     def __init__(self, omega, dL, eps_r, npml, bloch_phases=None):
#     def _grid_average_2d(self, eps_vec):
#     def _make_A(self, eps_vec):
#     def _solve_fn(self, eps_vec, entries_a, indices_a, Mz_vec):
#     def __init__(self, omega, dL, eps_r, omega_mod, delta, phi, Nsb, npml, bloch_phases=None):
#     def solve(self, source_z):
#     def _make_A(self, eps_vec, delta_matrix, phi_matrix):
#     def _solve_fn(self, eps_vec, entries_a, indices_a, Jz_vec):
#     def _Ez_to_Hx(self, Ez_vec):
#     def _Ez_to_Hy(self, Ez_vec):
#     def _Ez_to_Hx_Hy(self, Ez_vec):
#     def _vec_to_grid(self, vec):
#     def __init__(self, omega, dL, eps_r, npml, bloch_phases=None):
#     def _grid_average_3d(self, eps_vec):
#     def _make_A(self, eps_vec):
#     def _solve_fn(self, eps_vec, entries_a, indices_a, Mz_vec):
#         C = - 1 / MU_0 * self.Dxf.dot(self.Dxb) \
#             - 1 / MU_0 * self.Dyf.dot(self.Dyb)
#         M = 2*self.Nsb + 1
#         N = self.Nx * self.Ny
#         W = self.omega + npa.arange(-self.Nsb,self.Nsb+1)*self.omega_mod
#         C = sp.kron(sp.eye(M), - 1 / MU_0 * self.Dxf.dot(self.Dxb) - 1 / MU_0 * self.Dyf.dot(self.Dyb))
#         M = 2*self.Nsb + 1
#         N = self.Nx * self.Ny
#         W = self.omega + npa.arange(-self.Nsb,self.Nsb+1)*self.omega_mod 
#         P = sp.kron(sp.spdiags(W,[0],M,M), sp.eye(N))
#         M = 2*self.Nsb + 1
#         M = 2*self.Nsb + 1
#         C11 = -1 / MU_0 * self.Dyf.dot(self.Dyb)
#         C22 = -1 / MU_0 * self.Dxf.dot(self.Dxb)
#         C12 =  1 / MU_0 * self.Dyf.dot(self.Dxb)
#         C21 =  1 / MU_0 * self.Dxf.dot(self.Dyb)


import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spl
    import matplotlib.pylab as plt
from copy import deepcopy
from ceviche.constants import *
from ceviche.fdfd import compute_derivative_matrices#, Ez_to_Hx_Hy
    from ceviche.fdfd import fdfd_ez as fdfd

    # elements to keep, all for starts
    keep_elements = np.ones(values.shape)

    for f in filters:
        keep_f = f(values)
        keep_elements = np.logical_and(keep_elements, keep_f)

    # get the indeces you want to keep
    keep_indeces = np.where(keep_elements)[0]

    # filter and return arrays
    return values[keep_indeces], vectors[:, keep_indeces]


def normalize_modes(vectors):
    """ Normalize each `vec` in `vectors` such that `sum(|vec|^2)=1`
            vectors: array with shape (n_points, n_vectors)
        NOTE: eigs already normalizes for you, so you technically dont need this function
    """

    powers = np.sum(np.square(np.abs(vectors)), axis=0)

    return vectors / np.sqrt(powers)

def Ez_to_H(Ez, omega, dL, npml):
    """ Converts the Ez output of mode solver to Hx and Hy components
    """

    N = Ez.size

# Path: ceviche/fdfd.py
# class fdfd():
# class fdfd_ez(fdfd):
# class fdfd_hz(fdfd):
# class fdfd_mf_ez(fdfd):
# class fdfd_3d(fdfd):
#     def __init__(self, omega, dL, eps_r, npml, bloch_phases=None):
#     def eps_r(self):
#     def eps_r(self, new_eps):
#     def _make_A(self, eps_r):
#     def _solve_fn(self, entries_a, indices_a, source_vec):
#     def solve(self, source_z):
#     def _setup_derivatives(self):
#     def _setup_bloch_phases(self, bloch_phases):
#     def _vec_to_grid(self, vec):
#     def _grid_to_vec(self, grid):
#     def _save_shape(self, grid):
#     def _default_val(val, default_val=None):
#     def _Ex_Ey_to_Hz(self, Ex_vec, Ey_vec):
#     def _Ez_to_Hx(self, Ez_vec):
#     def _Ez_to_Hy(self, Ez_vec):
#     def _Ez_to_Hx_Hy(self, Ez_vec):
#     def _Hz_to_Ex(self, Hz_vec, eps_vec_xx):
#     def _Hz_to_Ey(self, Hz_vec, eps_vec_yy):
#     def _Hx_Hy_to_Ez(self, Hx_vec, Hy_vec, eps_vec_zz):
#     def _Hz_to_Ex_Ey(self, Hz_vec, eps_vec_xx, eps_vec_yy):
#     def __init__(self, omega, dL, eps_r, npml, bloch_phases=None):
#     def _make_A(self, eps_vec):
#     def _solve_fn(self, eps_vec, entries_a, indices_a, Jz_vec):
#     def __init__(self, omega, dL, eps_r, npml, bloch_phases=None):
#     def _grid_average_2d(self, eps_vec):
#     def _make_A(self, eps_vec):
#     def _solve_fn(self, eps_vec, entries_a, indices_a, Mz_vec):
#     def __init__(self, omega, dL, eps_r, omega_mod, delta, phi, Nsb, npml, bloch_phases=None):
#     def solve(self, source_z):
#     def _make_A(self, eps_vec, delta_matrix, phi_matrix):
#     def _solve_fn(self, eps_vec, entries_a, indices_a, Jz_vec):
#     def _Ez_to_Hx(self, Ez_vec):
#     def _Ez_to_Hy(self, Ez_vec):
#     def _Ez_to_Hx_Hy(self, Ez_vec):
#     def _vec_to_grid(self, vec):
#     def __init__(self, omega, dL, eps_r, npml, bloch_phases=None):
#     def _grid_average_3d(self, eps_vec):
#     def _make_A(self, eps_vec):
#     def _solve_fn(self, eps_vec, entries_a, indices_a, Mz_vec):
#         C = - 1 / MU_0 * self.Dxf.dot(self.Dxb) \
#             - 1 / MU_0 * self.Dyf.dot(self.Dyb)
#         M = 2*self.Nsb + 1
#         N = self.Nx * self.Ny
#         W = self.omega + npa.arange(-self.Nsb,self.Nsb+1)*self.omega_mod
#         C = sp.kron(sp.eye(M), - 1 / MU_0 * self.Dxf.dot(self.Dxb) - 1 / MU_0 * self.Dyf.dot(self.Dyb))
#         M = 2*self.Nsb + 1
#         N = self.Nx * self.Ny
#         W = self.omega + npa.arange(-self.Nsb,self.Nsb+1)*self.omega_mod 
#         P = sp.kron(sp.spdiags(W,[0],M,M), sp.eye(N))
#         M = 2*self.Nsb + 1
#         M = 2*self.Nsb + 1
#         C11 = -1 / MU_0 * self.Dyf.dot(self.Dyb)
#         C22 = -1 / MU_0 * self.Dxf.dot(self.Dxb)
#         C12 =  1 / MU_0 * self.Dyf.dot(self.Dxb)
#         C21 =  1 / MU_0 * self.Dxf.dot(self.Dyb)


# Path: ceviche/modes.py
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spl
    import matplotlib.pylab as plt
from copy import deepcopy
from ceviche.constants import *
from ceviche.fdfd import compute_derivative_matrices#, Ez_to_Hx_Hy
    from ceviche.fdfd import fdfd_ez as fdfd


    # elements to keep, all for starts
    keep_elements = np.ones(values.shape)

    for f in filters:
        keep_f = f(values)
        keep_elements = np.logical_and(keep_elements, keep_f)

    # get the indeces you want to keep
    keep_indeces = np.where(keep_elements)[0]

    # filter and return arrays
    return values[keep_indeces], vectors[:, keep_indeces]


def normalize_modes(vectors):
    """ Normalize each `vec` in `vectors` such that `sum(|vec|^2)=1`
            vectors: array with shape (n_points, n_vectors)
        NOTE: eigs already normalizes for you, so you technically dont need this function
    """

    powers = np.sum(np.square(np.abs(vectors)), axis=0)

    return vectors / np.sqrt(powers)

def Ez_to_H(Ez, omega, dL, npml):
    """ Converts the Ez output of mode solver to Hx and Hy components
    """

    N = Ez.size

    matrices = compute_derivative_matrices(omega, (N, 1), [npml, 0], dL=dL)