search

lbolla / EMpy

Electromagnetic Python
MIT License
194 stars 83 forks source link

Symetric Double Grating not acting appropriately #7

Closed strspeed closed 8 years ago

strspeed commented 8 years ago

class SymmetricDoubleGrating(object): in utils.py

Initial assumptions:

if..

Top = IsotropicMaterial('Top', n0=RefractiveIndex(n0_const=1.0))
Bottom = IsotropicMaterial('Bottom', n0=RefractiveIndex(n0_const=1.5))
..
EMpy.utils.SymmetricDoubleGrating(Top,Top,Bottom,.25,.25,580e-9,d*nanometers)
..

Should effectively act as a 50% DC binary grating. I tested this theory in G solver and in it this holds true.

EMpy's

EMpy.utils.AsymmetricDoubleGrating(Top,Top,Bottom,.25,.25,.0,580e-9,d*nanometers),
..
EMpy.utils.AsymmetricDoubleGrating(Top,Top,Bottom,.25,.25,1,580e-9,d*nanometers),

Both act effectively as a 50% DC binary solution

However, the EMpy solution does not follow this case, Results below with results of proposed fix...

Plot 1 is the results of the Binary grating. Under it, labeled figure 3, are the results from the Asymmetrical grating class, which are in agreement with our binary test. The Symmetrical grating class initial results are shown in figure two, and the results from the Symmetrical class fix are shown in the final plot labeled figure 4.

fix

Hot fix is as follows:

Find utils.py --> class

def getEPSFourierCoeffs(self, wl, n, anisotropic=True):
        """Return the Fourier coefficients of eps and eps**-1, orders [-n,n]."""
        nood = 2 * n + 1
        hmax = nood - 1
        if not anisotropic:
            # isotropic
            rix1 = self.mat1.n(wl)
            rix2 = self.mat2.n(wl)
            rix3 = self.mat3.n(wl)
            f1 = self.dc1
            f2 = self.dc2
            h = numpy.arange(-hmax, hmax + 1)
            N = len(h)
            A = -N / 4
            B = N  / 4
            print(A,B)
            EPS = (
                rix3 ** 2 * (h == 0) + (rix1 ** 2 - rix3 ** 2) * f1 *
                numpy.sinc(h * f1) * numpy.exp(2j * numpy.pi * h / N * A) +
                (rix2 ** 2 - rix3 ** 2) * f2 * numpy.sinc(h * f2) *
                numpy.exp(2j * numpy.pi * h / N * B)
            )

            EPS1 = (
                rix3 ** -2 * (h == 0) + (rix1 ** -2 - rix3 ** -2) * f1 *
                numpy.sinc(h * f1) * numpy.exp(2j * numpy.pi * h / N * A) +
                (rix2 ** -2 - rix3 ** -2) * f2 * numpy.sinc(h * f2) *
                numpy.exp(2j * numpy.pi * h / N * B)
            )
            return EPS, EPS1
        else:
            # anisotropic
            EPS = numpy.zeros((3, 3, 2 * hmax + 1), dtype=complex)
            EPS1 = numpy.zeros_like(EPS)
            eps1 = numpy.squeeze(
                self.mat1.epsilonTensor(wl)) / EMpy.constants.eps0
            eps2 = numpy.squeeze(
                self.mat2.epsilonTensor(wl)) / EMpy.constants.eps0
            eps3 = numpy.squeeze(
                self.mat3.epsilonTensor(wl)) / EMpy.constants.eps0
            f1 = self.dc1
            f2 = self.dc2
            h = numpy.arange(-hmax, hmax + 1)
            N = len(h)
            A = -N / 4.
            B = N / 4.
            for ih, hh in enumerate(h):
                EPS[:, :, ih] = (
                    (eps1 - eps3) * f1 * numpy.sinc(hh * f1) *
                    numpy.exp(2j * numpy.pi * hh / N * A) +
                    (eps2 - eps3) * f2 * numpy.sinc(hh * f2) *
                    numpy.exp(2j * numpy.pi * hh / N * B) +
                    eps3 * (hh == 0)
                )
                EPS1[:, :, ih] = (
                    (scipy.linalg.inv(eps1) - scipy.linalg.inv(eps3)) * f1 *
                    numpy.sinc(hh * f1) *
                    numpy.exp(2j * numpy.pi * hh / N * A) +
                    (scipy.linalg.inv(eps2) - scipy.linalg.inv(eps3)) * f2 *
                    numpy.sinc(hh * f2) *
                    numpy.exp(2j * numpy.pi * hh / N * B) +
                    scipy.linalg.inv(eps3) * (hh == 0)
                )
            return EPS, EPS1

Referencing the Asymmetric class I changed assignment statement of A,B vars to the following..

 def getEPSFourierCoeffs(self, wl, n, anisotropic=True):
        """Return the Fourier coefficients of eps and eps**-1, orders [-n,n]."""
        nood = 2 * n + 1
        hmax = nood - 1
        if not anisotropic:
            # isotropic
            rix1 = self.mat1.n(wl)
            rix2 = self.mat2.n(wl)
            rix3 = self.mat3.n(wl)
            f1 = self.dc1
            f2 = self.dc2
            h = numpy.arange(-hmax, hmax + 1)
            N = len(h)
            A = -N*f1 / 2
            B = N*f2 / 2
            EPS = (
                rix3 ** 2 * (h == 0) + (rix1 ** 2 - rix3 ** 2) * f1 *
                numpy.sinc(h * f1) * numpy.exp(2j * numpy.pi * h / N * A) +
                (rix2 ** 2 - rix3 ** 2) * f2 * numpy.sinc(h * f2) *
                numpy.exp(2j * numpy.pi * h / N * B)
            )

            EPS1 = (
                rix3 ** -2 * (h == 0) + (rix1 ** -2 - rix3 ** -2) * f1 *
                numpy.sinc(h * f1) * numpy.exp(2j * numpy.pi * h / N * A) +
                (rix2 ** -2 - rix3 ** -2) * f2 * numpy.sinc(h * f2) *
                numpy.exp(2j * numpy.pi * h / N * B)
            )
            return EPS, EPS1
        else:
            # anisotropic
            EPS = numpy.zeros((3, 3, 2 * hmax + 1), dtype=complex)
            EPS1 = numpy.zeros_like(EPS)
            eps1 = numpy.squeeze(
                self.mat1.epsilonTensor(wl)) / EMpy.constants.eps0
            eps2 = numpy.squeeze(
                self.mat2.epsilonTensor(wl)) / EMpy.constants.eps0
            eps3 = numpy.squeeze(
                self.mat3.epsilonTensor(wl)) / EMpy.constants.eps0
            f1 = self.dc1
            f2 = self.dc2
            h = numpy.arange(-hmax, hmax + 1)
            N = len(h)
            A = -N / 4.
            B = N / 4.
            for ih, hh in enumerate(h):
                EPS[:, :, ih] = (
                    (eps1 - eps3) * f1 * numpy.sinc(hh * f1) *
                    numpy.exp(2j * numpy.pi * hh / N * A) +
                    (eps2 - eps3) * f2 * numpy.sinc(hh * f2) *
                    numpy.exp(2j * numpy.pi * hh / N * B) +
                    eps3 * (hh == 0)
                )
                EPS1[:, :, ih] = (
                    (scipy.linalg.inv(eps1) - scipy.linalg.inv(eps3)) * f1 *
                    numpy.sinc(hh * f1) *
                    numpy.exp(2j * numpy.pi * hh / N * A) +
                    (scipy.linalg.inv(eps2) - scipy.linalg.inv(eps3)) * f2 *
                    numpy.sinc(hh * f2) *
                    numpy.exp(2j * numpy.pi * hh / N * B) +
                    scipy.linalg.inv(eps3) * (hh == 0)
                )
            return EPS, EPS1

I have not tested Anisotropic just yet, but am working on getting it up shortly. By chance can you reference any literature on getting the Fourier coefficients? I read the following work and see the similarities, but I'm at a loss of how the multiple indices come together with relation to duty cycle. https://www.osapublishing.org/josaa/abstract.cfm?uri=josaa-12-5-1068

lbolla commented 8 years ago

Thanks for the analysis! Can you please issue a pull request? It's very difficult to spot the code differences.

lbolla commented 8 years ago

Fixed.