Fun with Electromagnetic Python: http://empy.sourceforge.net

moharam1995.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright: 2009 Elson Liu
# License: BSD
"""
Reproduces Fig. 2 in Moharam, Grann, Pommet, and Gaylord, "Formulation for
stable and efficient implementation of the rigorous coupled-wave analysis
of binary gratings", J. Opt. Soc. Am. A 12(5), pp. 1068–1076, 1995
(doi:10.1364/JOSAA.12.001068)
"""
import numpy, EMpy, pylab
from EMpy.materials import IsotropicMaterial, RefractiveIndex


def moharam(theta_deg,psi_deg,phi_deg,wavelength,ambient,pitch,fillfactor,thicknesses,substrate,norders):
    """
    Rigorous coupled-wave simulation of a binary grating as a function of
    grating thickness
    """
    theta = EMpy.utils.deg2rad(theta_deg)   # angle of incidence
    delta = EMpy.utils.deg2rad(0.)
    psi = EMpy.utils.deg2rad(psi_deg)       # polarization (0 => TE, 90 => TM)
    phi = EMpy.utils.deg2rad(phi_deg)       # azimuth (90 => planar diffraction)
    
    solutions = []
    for d in thicknesses:
        multilayer = EMpy.utils.Multilayer([ \
        EMpy.utils.Layer(ambient, numpy.inf), \
        EMpy.utils.BinaryGrating(ambient, substrate, fillfactor, pitch, d), \
        EMpy.utils.Layer(substrate, numpy.inf), \
        ])
        
        solution = EMpy.RCWA.IsotropicRCWA(multilayer, theta, delta, psi, phi, norders).solve(wavelength)
        solutions.append(solution)
    return solutions


def main():
    wavelength = numpy.array([1.55e-6])
    thicknesses = numpy.linspace(0.,5.,100)*wavelength
    multiplier = 1
    pitch = multiplier*wavelength
    
    norders = 15    # orders of diffraction
    
    ambient = IsotropicMaterial('ambient', n0=RefractiveIndex(n0_const=1.))
    substrate = IsotropicMaterial('substrate', n0=RefractiveIndex(n0_const=2.04))
    
    # TE (psi = 0)
    print "TE"
    TEsolutions = moharam(10.,0.,90.,wavelength,ambient,pitch,0.5,thicknesses,substrate,norders)
    TE = numpy.zeros(len(TEsolutions))
    for ss,s in enumerate(TEsolutions):
		# DE3 is transmission efficiency, n+1 gives the +1 diffracted order
        TE[ss] = s.DE3[norders+1,0]
    
    # TM (psi = 90)
    print "TM"
    TMsolutions = moharam(10.,90.,90.,wavelength,ambient,pitch,0.5,thicknesses,substrate,norders)
    TM = numpy.zeros(len(TMsolutions))
    for ss,s in enumerate(TMsolutions):
		# DE3 is transmission efficiency, n+1 gives the +1 diffracted order
        TM[ss] = s.DE3[norders+1,0]

    pylab.plot(thicknesses/wavelength, TE[:], 'b.-', \
                thicknesses/wavelength, TM[:], 'r.-')
    pylab.xlabel('Normalized groove depth $d/\lambda$')
    pylab.ylabel('Diffraction efficiency $DE_1$')
    
    pylab.legend(('TE', 'TM'))
    pylab.axis('tight')
    pylab.ylim([0,1])
    if multiplier == 1:
        pylab.title("First-order transmitted diffraction efficiency\n\n" + \
                    "$n_I = 1$, $n_{II} = 2.04$, $\\theta = 10^\circ$, $\\Lambda = \lambda_0$")
        pylab.savefig("moharam1995_lambda.png")
        pylab.savefig("moharam1995_lambda.svg")
    else:
        pylab.title("First-order transmitted diffraction efficiency\n\n" + \
                    "$n_I = 1$, $n_{II} = 2.04$, $\\theta = 10^\circ$, $\\Lambda =$ " + \
                    str(multiplier) + "$\lambda_0$")
        pylab.savefig("moharam1995_"+str(multiplier)+"lambda.png")
        pylab.savefig("moharam1995_"+str(multiplier)+"lambda.svg")
    pylab.show()


if __name__ == "__main__":
    main()

moharam1995a.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright: 2009 Elson Liu
# License: BSD
"""
Reproduces Fig. 3 in Moharam, Grann, Pommet, and Gaylord, "Formulation for
stable and efficient implementation of the rigorous coupled-wave analysis
of binary gratings", J. Opt. Soc. Am. A 12(5), pp. 1068–1076, 1995
(doi:10.1364/JOSAA.12.001068)
"""
import numpy, EMpy, pylab
from EMpy.materials import IsotropicMaterial, RefractiveIndex

def moharam(theta_deg,psi_deg,phi_deg,wavelength,ambient,pitch,fillfactor,thickness,substrate,nlist):
    """
    Rigorous coupled-wave simulation of a binary grating as a function of
    the number of retained diffracted orders
    """
    theta = EMpy.utils.deg2rad(theta_deg)   # angle of incidence
    delta = EMpy.utils.deg2rad(0.)
    psi = EMpy.utils.deg2rad(psi_deg)       # polarization (0 => TE, 90 => TM)
    phi = EMpy.utils.deg2rad(phi_deg)       # azimuth (90 => planar diffraction)

    multilayer = EMpy.utils.Multilayer([ \
                 EMpy.utils.Layer(ambient, numpy.inf), \
                 EMpy.utils.BinaryGrating(ambient, substrate, fillfactor, pitch, thickness), \
                 EMpy.utils.Layer(substrate, numpy.inf), \
    ])
    
    solutions = []
    for norders in nlist:
        print "\t" + str(norders)
        solution = EMpy.RCWA.IsotropicRCWA(multilayer, theta, delta, psi, phi, norders).solve(wavelength)
        solutions.append(solution)
    return solutions


def main():
    wavelength = numpy.array([1.55e-6])
    thicknessmultiplier = 0.5
    thickness = thicknessmultiplier*wavelength
    pitchmultiplier = 10
    pitch = pitchmultiplier*wavelength
    
    nlist = range(5,50,1)       # orders of diffraction
    
    ambient = IsotropicMaterial('ambient', n0=RefractiveIndex(n0_const=1.))
    substrate = IsotropicMaterial('substrate', n0=RefractiveIndex(n0_const=2.04))
    
    # TE (psi = 0)
    print "TE"
    TEsolutions = moharam(10.,0.,90.,wavelength,ambient,pitch,0.5,thickness,substrate,nlist)
    TE = numpy.zeros(len(TEsolutions))
    for ss,s in enumerate(TEsolutions):
        # DE3 is transmission efficiency, n+1 gives the +1 diffracted order
        TE[ss] = s.DE3[nlist[ss]+1,0]
    
    # TM (psi = 90)
    print "TM"
    TMsolutions = moharam(10.,90.,90.,wavelength,ambient,pitch,0.5,thickness,substrate,nlist)
    TM = numpy.zeros(len(TMsolutions))
    for ss,s in enumerate(TMsolutions):
        # DE3 is transmission efficiency, n+1 gives the +1 diffracted order
        TM[ss] = s.DE3[nlist[ss]+1,0]

    pylab.plot(nlist, TE[:], 'b.-', \
                nlist, TM[:], 'r.-')
    pylab.xlabel('Number of harmonics')
    pylab.ylabel('Diffraction efficiency $DE_1$')
    
    pylab.legend(('TE', 'TM'))
    pylab.axis('tight')
    pylab.ylim([0,1])
    if pitchmultiplier == 1:
        pylab.title("First-order transmitted diffraction efficiency\n\n" + \
                    "$n_I = 1$, $n_{II} = 2.04$, $\\theta = 10^\circ$, $\\Lambda = \lambda_0$, $d =$ " + \
                    str(thicknessmultiplier) + "$\lambda_0$")
    else:
        pylab.title("First-order transmitted diffraction efficiency\n\n" + \
                    "$n_I = 1$, $n_{II} = 2.04$, $\\theta = 10^\circ$, $\\Lambda =$ " + \
                    str(pitchmultiplier) + "$\lambda_0$, $d =$ " + \
                    str(thicknessmultiplier) + "$\lambda_0$")
    pylab.savefig("moharam1995_p"+str(pitchmultiplier)+"_d" + str(thicknessmultiplier) + ".png")
    pylab.savefig("moharam1995_p"+str(pitchmultiplier)+"_d" + str(thicknessmultiplier) + ".svg")
    pylab.show()

if __name__ == "__main__":
    main()