View on GitHub

EMpy

Electromagnetic Python

Download this project as a .zip file Download this project as a tar.gz file

EMpy

EMpy (Electromagnetic Python) (not to be confused with empy) is a suite of numerical algorithms widely used in electromagnetism.

The package contains:

This list, very short by now, will hopefully enlarge to include an FDTD and an interface to the many very good software used in electromagnetism out there.

Examples

More (and more up-to-date) examples are available in the source code.

Here are some simple examples of the EMpy's functionalities. More will come as soon as possible. Impatients can look at the examples in the distribution source.

Single Layer Anti-Reflection Coating

The isotropic transfer matrix algorithm can be used to study anti-reflection coatings. A very simple single layer anti-reflection coating, described here, is computed by the following script.

#!/usr/bin/python

# taken from http://hyperphysics.phy-astr.gsu.edu/hbase/phyopt/antiref.html#c1

import numpy, EMpy, pylab

# define multilayer
n = numpy.array([1.,1.38,1.9044])
d = numpy.array([numpy.inf,387.5e-9/1.38,numpy.inf])
iso_layers = EMpy.utils.Multilayer()
for i in xrange(n.size):
    n0 = EMpy.materials.RefractiveIndex(n[i])
    iso_layers.append(EMpy.utils.Layer(EMpy.materials.IsotropicMaterial('mat', n0=n0).setRixFromConst(n[i]),d[i]))

# define incident wave plane
theta_inc = EMpy.utils.deg2rad(10.)
wls = numpy.linspace(0.85e-6,2.25e-6,300)

# solve
solution_iso   = EMpy.transfer_matrix.IsotropicTransferMatrix(iso_layers, theta_inc).solve(wls)

# plot
pylab.figure()
pylab.plot(wls, 10*numpy.log10(solution_iso.Rs), 'rx-', wls, 10*numpy.log10(solution_iso.Rp), 'g.-')
pylab.legend(('Rs', 'Rp'))
pylab.title('Single Layer Anti-Reflection Coating')
pylab.xlabel('wavelength /m')
pylab.ylabel('Power /dB')
pylab.grid()
pylab.xlim(wls.min(), wls.max())
pylab.savefig(__file__ + '.png')
pylab.show()

In the lines 8-12 the multilayer is defined. Indices and thicknesses are chosen so that a single layer of refractive index 1.38 and thickness about one quarter of the central wavelength 1.55um is positioned between two semi-infinite layers of air and glass. This gives theoretically no reflection for the central wavelength. The following graph is the result, which confirms the theoretical predictions.

Single Layer Anti-Reflection Coating

The graph shows two curves, one for each polarization. They differ because the incident plane wave is tilted by 10 degrees (line 15).

back to examples

Anisotropic Multilayer

The anisotropic transfer-matrix can be used to study the following example.

#!/usr/bin/python

import numpy, EMpy, pylab

# define the multilayer
epsilon = [1.0**2 * EMpy.constants.eps0 * numpy.eye(3), \
           EMpy.constants.eps0 * numpy.diag([2.1, 2.0, 1.9]), \
           2.3**2 * EMpy.constants.eps0 * numpy.eye(3), \
           4.3**2 * EMpy.constants.eps0 * numpy.eye(3), \
           3.0**2 * EMpy.constants.eps0 * numpy.eye(3)]

d = numpy.array([numpy.inf,1e-6,2.3e-6,0.1e-6,numpy.inf])

aniso_layers = EMpy.utils.Multilayer()
for i in xrange(len(epsilon)):
    aniso_layers.append(EMpy.utils.Layer(EMpy.materials.AnisotropicMaterial('Air').
                                         setEpsilonTensorFromConst(epsilon[i]), d[i]))

# define the planewave
theta_inc_x = EMpy.utils.deg2rad(0.)
theta_inc_y = 0.
wls = numpy.linspace(1.4e-6,1.7e-6,100)

# solve
solution_aniso = EMpy.transfer_matrix.AnisotropicTransferMatrix(aniso_layers, 
                                                                theta_inc_x, 
                                                                theta_inc_y).solve(wls)

# plot
pylab.figure()
pylab.plot(wls, solution_aniso.R[0,0,:], \
           wls, solution_aniso.R[1,0,:], \
           wls, solution_aniso.R[0,1,:], \
           wls, solution_aniso.R[1,1,:], \
           wls, solution_aniso.T[0,0,:], \
           wls, solution_aniso.T[1,0,:], \
           wls, solution_aniso.T[0,1,:], \
           wls, solution_aniso.T[1,1,:])
pylab.legend(('Rss', 'Rps','Rsp', 'Rpp','Tss', 'Tps','Tsp', 'Tpp'))
pylab.title('Anisotropic Multilayer')
pylab.xlabel('wavelength /m')
pylab.ylabel('Power /dB')
pylab.xlim(wls.min(), wls.max())
pylab.show()
The result is shown in the following graph. Anisotropic Multilayer

Four curves are clearly visible because the two polarizations behave in a different way, even if the incidence is normal (line 16-17), because of the anisotropic material (line 5).

back to examples

Single Ring Resonator

This is the script to study a single ring resonator:

import EMpy
import numpy
import pylab

wls = numpy.linspace(1.53e-6, 1.56e-6, 1000)
K1 = EMpy.devices.Coupler(wls, numpy.sqrt(0.08), 1.)
K2 = EMpy.devices.Coupler(wls, numpy.sqrt(0.08), 1.)
l1 = numpy.pi * 5e-6
l2 = numpy.pi * 5e-6
SWG = EMpy.devices.SWG(488, 220, 25).solve(wls)
SRR = EMpy.devices.SRR(K1, K2, SWG.neff, l1, l2).solve()

pylab.plot(wls, numpy.absolute(SRR.THRU), 'r.-',
           wls, numpy.absolute(SRR.DROP), 'g.-')
pylab.axis('tight')
pylab.ylim([0, 1])
pylab.xlabel('wavelength /m')
pylab.ylabel('power')
pylab.legend(('THRU', 'DROP'))
pylab.show()

In lines 4 and 5 the two coupling coefficients and the two coupler losses are defined. The length of the ring is defined in line 6, while its effective index is supposed to be the one of a silicon rib waveguide of 488nm x 220nm at room temperature.

This is its spectrum: SingleRing Resonator

The maximum drop power is 0 dB because no losses are considered, and periodic with the period being the FSR of the ring.

back to examples

N-Rings Resonator

This is the script to study a ring resonator made of 3 rings:

import EMpy
import numpy
import pylab

wls = numpy.linspace(1.53e-6, 1.57e-6, 1000)

Ks = [EMpy.devices.Coupler(wls, numpy.sqrt(0.08), 1.),
      EMpy.devices.Coupler(wls, numpy.sqrt(0.008), 1.),
      EMpy.devices.Coupler(wls, numpy.sqrt(0.006), 1.),
      EMpy.devices.Coupler(wls, numpy.sqrt(0.09), 1.)]

R = 5e-6
l1s = [numpy.pi * R, numpy.pi * R, numpy.pi * R]
l2s = [numpy.pi * R, numpy.pi * R, numpy.pi * R]

SWG = EMpy.devices.SWG(400, 220, 125).solve(wls)
neffs = [SWG.neff, SWG.neff, SWG.neff]

NRR = EMpy.devices.NRR(Ks, neffs, l1s, l2s).solve()

pylab.plot(wls, 20 * numpy.log10(numpy.absolute(NRR.THRU)), 'r.-',
           wls, 20 * numpy.log10(numpy.absolute(NRR.DROP)), 'g.-')
pylab.axis('tight')
pylab.ylim([-30, 0])
pylab.xlabel('wavelength /m')
pylab.ylabel('power /dB')
pylab.legend(('THRU', 'DROP'))
pylab.show()

After the "imports", a list of couplers is defined: four, in this case, because the resonator is made of 3 rings. Then, the arcs connecting the resonators are defined: the rings are supposed to be all equal, and the couplers are on the diameters of the rings. Finally, the rings are supposed to have an effective index equal to a 400nm x 200nm straight waveguide, at 125 degrees.

This is its spectrum: N-Ring Resonator

The drop is always lower than 0 dB because the coupling coefficients are not optimized.

back to examples

Finite Difference Fully Vectorial Modesolver

This is the script to find the first two modes of a channel SOI waveguide:

import numpy
import EMpy

def epsfunc(x, y):
    [X, Y] = numpy.meshgrid(x, y)
    return numpy.where((numpy.abs(X.T - 1.24e-6) <= .24e-6) *
                       (numpy.abs(Y.T - 1.11e-6) <= .11e-6),
                       3.4757**2,
                       1.446**2)

wl = 1.55e-6
x = numpy.linspace(0, 2.48e-6, 125)
y = numpy.linspace(0, 2.22e-6, 112)

neigs = 2
tol = 1e-8
boundary = '0000' 

solver = EMpy.modesolvers.VFDModeSolver(wl, x, y, epsfunc, boundary).solve(neigs, tol)
print [m.neff for m in solver.modes]

Effective indices are printed on stdout.

This is a screenshot of the first (TE) mode: Mode Channel WG

back to examples

Multilayer with grating

This is the script to study the reflectivity of a grating embedded in an anisotropic multilayer. The multilayer is made of 7 layers, the 4th being a binary grating made with an anisotropic (SiN, rotated with the optical axis not aligned with the coordinate axis) and an isotropic (BPTEOS) material. Incidence is not normal. A comparison of the reflectivity is made with a similar multilayer, where the grating layer has been substituted by an "effective" isotropic layer.

#!/usr/bin/python

import numpy, EMpy, pylab

alpha = 0.
delta = 0.
# psi = EMpy.utils.deg2rad(0.)                       # TM
# psi = EMpy.utils.deg2rad(90.)                      # TE
psi = EMpy.utils.deg2rad(70.)                      # hybrid
phi = EMpy.utils.deg2rad(90.)

LAMBDA = 1016e-9                        # grating periodicity
n = 2                                   # orders of diffraction

UV6    = EMpy.materials.IsotropicMaterial('UV6').setRixFromConst(1.560)
SiN    = EMpy.materials.AnisotropicMaterial('SiN').setEpsilonTensorFromConst(EMpy.constants.eps0 * \
         EMpy.utils.euler_rotate(numpy.diag(numpy.asarray([1.8550, 1.8750, 1.9130])**2), \
                                 EMpy.utils.deg2rad(14), \
                                 EMpy.utils.deg2rad(25), \
                                 EMpy.utils.deg2rad(32)))
BPTEOS = EMpy.materials.IsotropicMaterial('BPTEOS').setRixFromConst(1.448)
ARC1   = EMpy.materials.IsotropicMaterial('ARC1').setRixFromConst(1.448)

EFF    = EMpy.materials.IsotropicMaterial('EFF').setRixFromConst(1.6)

multilayer1 = EMpy.utils.Multilayer([ \
    EMpy.utils.Layer(EMpy.materials.Air, numpy.inf), \
    EMpy.utils.Layer(SiN, 226e-9), \
    EMpy.utils.Layer(BPTEOS, 226e-9), \
    EMpy.utils.BinaryGrating(SiN, BPTEOS, .659, LAMBDA, 123e-9), \
    EMpy.utils.Layer(SiN, 219e-9), \
    EMpy.utils.Layer(EMpy.materials.SiO2, 2188e-9), \
    EMpy.utils.Layer(EMpy.materials.Si, numpy.inf), \
    ])

multilayer2 = EMpy.utils.Multilayer([ \
    EMpy.utils.Layer(EMpy.materials.Air, numpy.inf), \
    EMpy.utils.Layer(SiN, 226e-9), \
    EMpy.utils.Layer(BPTEOS, 226e-9), \
    EMpy.utils.Layer(EMpy.materials.IsotropicMaterial().setRixFromConst(1.6), 123e-9), \
    EMpy.utils.Layer(SiN, 219e-9), \
    EMpy.utils.Layer(EMpy.materials.SiO2, 2188e-9), \
    EMpy.utils.Layer(EMpy.materials.Si, numpy.inf), \
    ])

wls = numpy.linspace(1.50e-6, 1.60e-6, 301)

solution1 = EMpy.RCWA.AnisotropicRCWA(multilayer1, alpha, delta, psi, phi, n).solve(wls)
solution2 = EMpy.RCWA.AnisotropicRCWA(multilayer2, alpha, delta, psi, phi, n).solve(wls)

pylab.plot(wls, solution1.DEO1[n,:], 'k.-', \
           wls, solution1.DEO3[n,:], 'r.-', \
           wls, solution1.DEE1[n,:], 'b.-', \
           wls, solution1.DEE3[n,:], 'g.-', \
           wls, solution2.DEO1[n,:], 'k--', \
           wls, solution2.DEO3[n,:], 'r--', \
           wls, solution2.DEE1[n,:], 'b--', \
           wls, solution2.DEE3[n,:], 'g--', \
           )
pylab.xlabel('wavelength')
pylab.ylabel('diffraction efficiency')
pylab.legend(('DEO1', 'DEO3', 'DEE1', 'DEE3'))
pylab.axis('tight')
pylab.ylim([0,1])
pylab.show()

The reflectivity are shown in the picture below (in dashed lines is the reflectivity withouth the grating): the effect of the grating is clearly visible as a peak in the reflectivity. RCWA