DE Solutions don't follow other applications solutions

[strspeed]

I've found the EMpy solutions break down when changing the wave angle. I've attached my binary grating test below.

The following example outputs a solution which is the same as the Gsolver solution for the same binary grating. The Gsolver parameters are as follows: Theta: 0 Phi: 0 Alpha: 90 Beta: 0

import numpy
import pylab

import EMpy
from EMpy.materials import IsotropicMaterial, RefractiveIndex

alpha = EMpy.utils.deg2rad(0.) #Angle of incidence
delta = EMpy.utils.deg2rad(0.) 
psi = EMpy.utils.deg2rad(0.)  # TE
phi = EMpy.utils.deg2rad(90.)  #azimuth

LAMBDA = 900e-9

ns = 1 # orders of diffraction

Top = IsotropicMaterial('Top', n0=RefractiveIndex(n0_const=1.0))
Bottom = IsotropicMaterial('Bottom', n0=RefractiveIndex(n0_const=1.45032))

solutions = []

wls = numpy.linspace(.909e-6, 1.500e-6, 30)

multilayer = EMpy.utils.Multilayer([
        EMpy.utils.Layer(Top, numpy.inf),
        #EMpy.utils.SymmetricDoubleGrating(Top,Top,Bottom,.25,.25, LAMBDA, 1200e-9),
        EMpy.utils.BinaryGrating(Top,Bottom, .5, LAMBDA, 1200e-9),
        EMpy.utils.Layer(Bottom, numpy.inf),
    ])

solution = EMpy.RCWA.AnisotropicRCWA(
    multilayer, alpha, delta, psi, phi, ns).solve(wls)

#Ns +1 == Gsolver -1?
#Ns == Gsolver 0
#Ns -1 == Gsolver +1?

n=0
pylab.plot(wls*1e9, solution.DEE1[ns+n, :], 'ko-',
           wls*1e9, solution.DEE3[ns+n, :], 'ro-',
           wls*1e9, solution.DEO1[ns+n, :], 'ko-',
           wls*1e9, solution.DEO3[ns+n, :], 'ro-',
           #wls*1e9, solution.DEE1[2, :], 'kx-',
           #wls*1e9, solution.DEE3[2, :], 'rx-',
           #wls*1e9, solution.DE1[ns + 1, 0], 'k.-',
           #wls*1e9, solution.DE3[ns + 1, 0], 'r.-',

           )
pylab.xlabel('wavelength /m')
pylab.ylabel('diffraction efficiency')
pylab.legend(('DE1:0', 'DE3:0', 'DE1:1', 'DE3:1', 'DE1:+2', 'DE3:+2'))
pylab.axis('tight')
pylab.ylim([0, 1])
pylab.show()

If we change our phi and Gsolvers Alpha to 0 degrees our solutions no longer follow.

We update the above code to have..

alpha = EMpy.utils.deg2rad(0.) #Angle of incidence
delta = EMpy.utils.deg2rad(0.) 
psi = EMpy.utils.deg2rad(0.)  # TE
phi = EMpy.utils.deg2rad(0.)  #azimuth

Gsolver parameters as follows: Theta: 0 Phi: 0 Alpha: 0 Beta: 0

What are your thoughts, Is my code incorrect?

[lbolla]

What result does GSolver give you? Are you sure the angles in GSolver match the definitions of the angles in EMpy's RCWA solver (Glytsis article)?

Given that you are using isotropic materials, have you tried comparing the results with IsotropicRCWA?

[strspeed]

I'll attach the results of both shortly. Regarding the angles, attached is an image of the Gsolver dialog annotated with my understanding of the EMpy RCWA angels. Let the green line be the incoming wave, and the red line coming off it is the polarization. As shown the wave is in S-Polarization or TE. My understanding is Alpha is the angle of incidence (0-180 deg) and Phi is the azimuth (0-360 deg).

Edit: Actually I'll hold off posting my results. I feel I might not understand the angle relationship here.

I see the source explanation of the angles. I don't have access to the article however.

[strspeed]

TL:DC: When n>0 EMpy solutions do not match Gsolvers solutions for sub-wavelength periods. EMpy solutions are always 0 when n>0 for sub-wavelength periods.

I figured out my angle issues thankfully and got exact matches! However, I'm seeing some major issues with sub wavelength periods, and/or the diffraction orders retained in memory.

I'm operating EMpy following that LAMBDA is the width of one period of a grating. Such that LAMBDA = 500e-9 means in a binary grating that a single period is 500 nanometers wide.

The Following code produces a solution which is exact to the one Gsolver produces should we let both N retained be 0 | n=0. The following code sweeps from a period of 200nm to 1000nm.

"""Rigorous Coupled Wave Analysis example.
Inspired by Moharam, "Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of
binary gratings", JOSA A, 12(5), 1995
"""

import numpy
import pylab

import EMpy
from EMpy.materials import IsotropicMaterial, RefractiveIndex

alpha = numpy.deg2rad(45)  #EMpy.utils.deg2rad(40.)
delta = EMpy.utils.deg2rad(0.) #Gsolver psi
psi = EMpy.utils.deg2rad(0.)  # TE 0  == Gsolver alpha=90
phi = EMpy.utils.deg2rad(90.) #Remains Unknown

wl = numpy.array([1000.0e-9])
ds = numpy.linspace(200e-9, 1000e-9, 1000)

#LAMBDA = numpy.array([250.0e-9])

VARLAM = numpy.linspace(200e-9, 1000e-9, 1000)
n = 0  # orders of diffraction

Top = IsotropicMaterial('Top', n0=RefractiveIndex(n0_const=1.))
Bottom = IsotropicMaterial('Bottom', n0=RefractiveIndex(n0_const=1.5))

solutions = []

for d in VARLAM:

    multilayer = EMpy.utils.Multilayer([
        EMpy.utils.Layer(Top, numpy.inf),
        EMpy.utils.BinaryGrating(Top, Bottom,   .5, d, 1000e-9),
        EMpy.utils.Layer(Bottom, numpy.inf),
    ])

    solution = EMpy.RCWA.IsotropicRCWA(multilayer, alpha, delta, psi, phi, n).solve(wl)
    solutions.append(solution)

DE1 = numpy.zeros(len(solutions))
DE3 = numpy.zeros(len(solutions))

offset = 0

for x in solutions:
    print(x.DE3)
for ss, s in enumerate(solutions):
    DE1[ss] = s.DE1[n+offset, 0]
    DE3[ss] = s.DE3[n+offset, 0]

pylab.plot(VARLAM*1.0e9 , DE1[:], 'k.-',
           VARLAM*1.0e9 , DE3[:], 'r.-')
pylab.xlabel('normalized groove depth')
pylab.ylabel('diffraction efficiency')
pylab.legend(('DE1', 'DE3'))
pylab.axis('tight')
pylab.ylim([0, 1])
pylab.show()

The result of said code is shown:

Next we will alter our Gsolver N to be 1. or in other words we retain 0R 0T, -1R,-1T We do the same with our EMpy code. Shown additionally is the results from our first run where n=0. No other changes are made to our code.

We see our EMpy N=1 has flatlines to absolute zero. However, our Gsolver results show a decrease in power as we approach LAMBDA = 1000.

Working.. The next test changes the total height of the grating. We check that EMpy's solutions are the same as Gsolvers for both n=0 and n=1. Which they are as shown!

"""Rigorous Coupled Wave Analysis example.
Inspired by Moharam, "Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of
binary gratings", JOSA A, 12(5), 1995
"""

import numpy
import pylab

import EMpy
from EMpy.materials import IsotropicMaterial, RefractiveIndex

alpha = numpy.deg2rad(45)  #EMpy.utils.deg2rad(40.)
delta = EMpy.utils.deg2rad(0.) #Gsolver psi
psi = EMpy.utils.deg2rad(90.)  # TE  == Gsolver alpha=90
phi = EMpy.utils.deg2rad(90.) #Gsolver Beta

wl = numpy.array([1000e-9])
ds = numpy.linspace(500e-9, 1000e-9, 500)
LAMBDA = numpy.array([1e-5])

n = 0  # orders of diffraction

Top = IsotropicMaterial('Top', n0=RefractiveIndex(n0_const=1.))
Bottom = IsotropicMaterial('Bottom', n0=RefractiveIndex(n0_const=1.5))

solutions = []
for d in ds:
    multilayer = EMpy.utils.Multilayer([
        EMpy.utils.Layer(Top, numpy.inf),
        EMpy.utils.BinaryGrating(Top, Bottom, .5, LAMBDA, d),
        EMpy.utils.Layer(Bottom, numpy.inf),
    ])

    solution = EMpy.RCWA.IsotropicRCWA(multilayer, alpha, delta, psi, phi, n).solve(wl)
    solutions.append(solution)

DE1 = numpy.zeros(len(solutions))
DE3 = numpy.zeros(len(solutions))
for ss, s in enumerate(solutions):
    DE1[ss] = s.DE1[n, 0]
    DE3[ss] = s.DE3[n, 0]

pylab.plot(ds*1e9, DE1[:], 'k.-',
           ds*1e9, DE3[:], 'r.-')
pylab.xlabel('normalized groove depth')
pylab.ylabel('diffraction efficiency')
pylab.legend(('DE1', 'DE3'))
pylab.axis('tight')
pylab.ylim([0, 1])
pylab.show()

|

For the next sets of tests we vary the period from sub wavelength height of 500nm to 1000nm. Our first test of n=0 produces the same result in both programs. However, when we change into n=1 EMpy's solutions flatline to zero. Note: I tested up to n=15 in both programs. EMpy gave the same results shown for all, Gsolver did the same.

"""Rigorous Coupled Wave Analysis example.
Inspired by Moharam, "Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of
binary gratings", JOSA A, 12(5), 1995
"""

import numpy
import pylab

import EMpy
from EMpy.materials import IsotropicMaterial, RefractiveIndex

alpha = numpy.deg2rad(45)  #EMpy.utils.deg2rad(40.)
delta = EMpy.utils.deg2rad(0.) #Gsolver psi
psi = EMpy.utils.deg2rad(90.)  # TE  == Gsolver alpha=90
phi = EMpy.utils.deg2rad(90.) #Gsolver Beta

wl = numpy.array([1000e-9])
ds = numpy.linspace(500e-9, 1000e-9, 500)
LAMBDA = numpy.array([1000e-9])

n = 1  # orders of diffraction

Top = IsotropicMaterial('Top', n0=RefractiveIndex(n0_const=1.))
Bottom = IsotropicMaterial('Bottom', n0=RefractiveIndex(n0_const=1.5))

solutions = []
for d in ds:
    multilayer = EMpy.utils.Multilayer([
        EMpy.utils.Layer(Top, numpy.inf),
        EMpy.utils.BinaryGrating(Top, Bottom, .5, LAMBDA, d),
        EMpy.utils.Layer(Bottom, numpy.inf),
    ])

    solution = EMpy.RCWA.IsotropicRCWA(multilayer, alpha, delta, psi, phi, n).solve(wl)
    solutions.append(solution)

DE1 = numpy.zeros(len(solutions))
DE3 = numpy.zeros(len(solutions))
for ss, s in enumerate(solutions):
    DE1[ss] = s.DE1[n, 0]
    DE3[ss] = s.DE3[n, 0]

pylab.plot(ds*1e9, DE1[:], 'k.-',
           ds*1e9, DE3[:], 'r.-')
pylab.xlabel('normalized groove depth')
pylab.ylabel('diffraction efficiency')
pylab.legend(('DE1', 'DE3'))
pylab.axis('tight')
pylab.ylim([0, 1])
pylab.show()

A final test i ran was sweeping the period from 200nm to 10,000nm with n=1

[lbolla]

Interesting. Using the AnisotropicRCWA solver I think I get the same results as GSolver. Try to use solution = EMpy.RCWA.AnisotropicRCWA(multilayer, alpha, delta, psi, phi, n).solve(wl) and let me know if you agree.

[strspeed]

Hi,

That fixed it! I've included now my initial testing of TM polarization which exhibits the same issues in EMpy both solvers (Iso, Ano). Attached is my TM code. The first image is a TE solution in both applications, the second is TM solutions in both applications. I'll do a better analysis later.

"""Rigorous Coupled Wave Analysis example.
Inspired by Moharam, "Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of
binary gratings", JOSA A, 12(5), 1995
"""

import numpy
import pylab

import EMpy
from EMpy.materials import IsotropicMaterial, RefractiveIndex

alpha = numpy.deg2rad(45)  #EMpy.utils.deg2rad(40.)
delta = EMpy.utils.deg2rad(0.) #Gsolver psi
psi = EMpy.utils.deg2rad(90.)  # TE 0  == Gsolver alpha=90

##First image TE EMpy psi = 0rad / Gsolver Alpha 90##
##Second Image psi = EMpy.utils.deg2rad(90.) / Gsolver Alpha = 0 ##

phi = EMpy.utils.deg2rad(90.) #Remains Unknown

wl = numpy.array([1000.0e-9])
ds = numpy.linspace(200e-9, 1000e-9, 1000)

#LAMBDA = numpy.array([250.0e-9])

VARLAM = numpy.linspace(200e-9, 1000e-9, 1000)
n = 2  # orders of diffraction

Top = IsotropicMaterial('Top', n0=RefractiveIndex(n0_const=1.))
Bottom = IsotropicMaterial('Bottom', n0=RefractiveIndex(n0_const=1.5))

solutions = []

for d in VARLAM:

    multilayer = EMpy.utils.Multilayer([
        EMpy.utils.Layer(Top, numpy.inf),
        EMpy.utils.BinaryGrating(Top, Bottom,   .5, 500e-9, d ),
        EMpy.utils.Layer(Bottom, numpy.inf),
    ])

    solution = EMpy.RCWA.AnisotropicRCWA(multilayer, alpha, delta, psi, phi, n).solve(wl)
    solutions.append(solution)

DE1 = numpy.zeros(len(solutions))
DE3 = numpy.zeros(len(solutions))

offset = 0

for ss, s in enumerate(solutions):
    DE1[ss] = s.DEE1[n+offset, 0]
    DE3[ss] = s.DEE3[n+offset, 0]

pylab.plot(VARLAM*1.0e9 , DE1[:], 'k.-',
           VARLAM*1.0e9 , DE3[:], 'r.-')
pylab.xlabel('normalized groove depth')
pylab.ylabel('diffraction efficiency')
pylab.legend(('DE1', 'DE3'))
pylab.axis('tight')
pylab.ylim([0, 1])
pylab.show()

[lbolla]

Right, so now the issue is matching results for n=1 in Iso and Aniso. Maybe you are able to create a very very simple example (and unittest) to expose the issue? Are you willing to contribute fixing this issue?

[strspeed]

Sure, I'll give it a shot and see what I can turn up!

[strspeed]

I believe I found the issue!

RCWA.py --> dispersion_relation_ordinary()

def dispersion_relation_ordinary(kx, ky, k, nO):
    """Dispersion relation for the ordinary wave.

    NOTE
    See eq. 15 in Glytsis, "Three-dimensional (vector) rigorous
    coupled-wave analysis of anisotropic grating diffraction",
    JOSA A, 7(8), 1990 **Always give positive real or negative
    imaginary.**
    """

    if kx.shape != ky.shape:
        raise ValueError('kx and ky must have the same length')

    delta = (k * nO) ** 2 - (kx ** 2 + ky ** 2)
    if S.all(delta > 0):
        kz = S.sqrt(delta)
    else:
        kz = -1j * S.sqrt(-delta)
    return kz

Code comment states solutions must be positive real or negative imaginary. However, the return of this function will return a negative real number. If we fix this our solutions for IsotropicRCWA() n>0 follow Gsolver perfectly in my initial testing!

def dispersion_relation_ordinary(kx, ky, k, nO):
    """Dispersion relation for the ordinary wave.

    NOTE
    See eq. 15 in Glytsis, "Three-dimensional (vector) rigorous
    coupled-wave analysis of anisotropic grating diffraction",
    JOSA A, 7(8), 1990 Always give positive real or negative
    imaginary.
    """

    if kx.shape != ky.shape:
        raise ValueError('kx and ky must have the same length')

    delta = (k * nO) ** 2 - (kx ** 2 + ky ** 2)
    else:
        kz = -1j * S.sqrt(-delta)
        kz.real = abs(kz.real)
    return kz

[strspeed]

The above solution works for all my tests with abstract grating shapes of a simple refractive index. The above hotfix fails to produce similar computations for wire mesh polarizers. Multilayer (Air, Binary(gold,air), const 1.5) .

Do note that I have had experimental DE results to compare to up until this point. EMpy's solutions make more sense at first glance. I'll have to see if I can get access to some polarizer results in the long run. In the meantime I'll look for more anomalies.

[lbolla]

Thanks! I've pushed a fix for it. Can you pull the latest version and give it whirl? The problem, I think, was that delta can be a complex vector, but I was in fact assuming that it was a real scalar. Now, I just take the sqrt and adjust the real/imag part to fit the requirements from the referenced article. It should work for Gold, too, which should have a complex delta. Let me know your findings!

And thanks a lot for the debugging.

[strspeed]

No problem. Everything I tested looks fantastic and all solutions are lining up. Great work on this library, its fantastic!