# Path: EMpy/utils.py
# def snell(theta_inc, n):
#     """Snell law.
# 
#     INPUT
#     theta_inc = angle of incidence.
#     n = 1D numpy.array of refractive indices.
# 
#     OUTPUT
#     theta = 1D numpy.array.
#     """
# 
#     theta = numpy.zeros_like(n)
#     theta[0] = theta_inc
#     for i in range(1, n.size):
#         theta[i] = numpy.arcsin(n[i - 1] / n[i] * numpy.sin(theta[i - 1]))
#     return theta
# 
# def norm(x):
#     """Return the norm of a 1D vector."""
#     return numpy.sqrt(numpy.vdot(x, x))
# 
# Path: EMpy/constants.py



from builtins import zip
from builtins import range
from builtins import object
from scipy.linalg import inv
from EMpy.utils import snell, norm
from EMpy.constants import c, mu0
import scipy as S

            gamma = S.roots(coeff)
            tmp = S.sort_complex(gamma)
            gamma = tmp[[3, 0, 2, 1]]  # convention

            k = (
                k0
                * S.array(
                    [alpha * S.ones(gamma.shape), beta * S.ones(gamma.shape), gamma]
                ).T
            )
            v = S.zeros((4, 3), dtype=complex)

            for i, g in enumerate(gamma):

                # H = K + [
                #     [-beta2 - g ** 2, alpha * beta, alpha * g],
                #     [alpha * beta, -alpha2 - g ** 2, beta * g],
                #     [alpha * g, beta * g, -alpha2 - beta2],
                # ]
                v[i, :] = [
                    (K[1, 1] - alpha2 - g ** 2) * (K[2, 2] - alpha2 - beta2)
                    - (K[1, 2] + beta * g) ** 2,
                    (K[1, 2] + beta * g) * (K[2, 0] + alpha * g)
                    - (K[0, 1] + alpha * beta) * (K[2, 2] - alpha2 - beta2),
                    (K[0, 1] + alpha * beta) * (K[1, 2] + beta * g)
                    - (K[0, 2] + alpha * g) * (K[1, 1] - alpha2 - g ** 2),
                ]

            p3 = v[0, :]

# Path: EMpy/utils.py
# def snell(theta_inc, n):
#     """Snell law.
# 
#     INPUT
#     theta_inc = angle of incidence.
#     n = 1D numpy.array of refractive indices.
# 
#     OUTPUT
#     theta = 1D numpy.array.
#     """
# 
#     theta = numpy.zeros_like(n)
#     theta[0] = theta_inc
#     for i in range(1, n.size):
#         theta[i] = numpy.arcsin(n[i - 1] / n[i] * numpy.sin(theta[i - 1]))
#     return theta
# 
# def norm(x):
#     """Return the norm of a 1D vector."""
#     return numpy.sqrt(numpy.vdot(x, x))
# 
# Path: EMpy/constants.py



# Path: EMpy/transfer_matrix.py
from builtins import zip
from builtins import range
from builtins import object
from scipy.linalg import inv
from EMpy.utils import snell, norm
from EMpy.constants import c, mu0
import scipy as S


            gamma = S.roots(coeff)
            tmp = S.sort_complex(gamma)
            gamma = tmp[[3, 0, 2, 1]]  # convention

            k = (
                k0
                * S.array(
                    [alpha * S.ones(gamma.shape), beta * S.ones(gamma.shape), gamma]
                ).T
            )
            v = S.zeros((4, 3), dtype=complex)

            for i, g in enumerate(gamma):

                # H = K + [
                #     [-beta2 - g ** 2, alpha * beta, alpha * g],
                #     [alpha * beta, -alpha2 - g ** 2, beta * g],
                #     [alpha * g, beta * g, -alpha2 - beta2],
                # ]
                v[i, :] = [
                    (K[1, 1] - alpha2 - g ** 2) * (K[2, 2] - alpha2 - beta2)
                    - (K[1, 2] + beta * g) ** 2,
                    (K[1, 2] + beta * g) * (K[2, 0] + alpha * g)
                    - (K[0, 1] + alpha * beta) * (K[2, 2] - alpha2 - beta2),
                    (K[0, 1] + alpha * beta) * (K[1, 2] + beta * g)
                    - (K[0, 2] + alpha * g) * (K[1, 1] - alpha2 - g ** 2),
                ]

            p3 = v[0, :]

            p3 /= norm(p3)