Simulation Structure

Simulations with EMUstack are generally carried out using a python script file. This file is kept in its own directory which is placed in the EMUstack directory. All results of the simulation are automatically created within this directory. This directory then serves as a complete record of the calculation. Often, we will also save the simulation objects (scattering matrices, propagation constants etc.) within this folder for future inspection, manipulation, plotting, etc. Traditionally the name of the python script file begins with simo_. This is convenient for setting terminal alias’ for running the script. Throughout the tutorial the script file will be called simo.py.

To start a simulation open a terminal and change into the directory containing the simo.py file. To run this script:

$ python simo.py

To have direct access to the simulation objects upon the completion of the script use,:

$ python -i simo.py

This will return you into an interactive python session in which all simulation objects are accessible. In this session you can access the docstrings of objects, classes and methods. For example:

>>> from pydoc import help
>>> help(objects.Light)

where we have accessed the docstring of the Light class from objects.py

In the remainder of the guide we go through a number of example simo.py files. These cover a wide range (though non-exhaustive) of established applications of EMUstack. The source files for these examples are in EMUstack/examples/ The first 8 examples are pretty essential for using EMUstack, while those thereafter show EMUstack applied to a number of (IMHO) interesting situations.

Another tip to mention before diving into the examples is running simulations within Screen Sessions. These allow you to disconnect from the terminal instance and are discusses in Screen Sessions.

Single Interface

"""
Simulating an interface between 2 homogeneous, non-dispersive media.
"""

import time
import datetime
import numpy as np
import sys
from multiprocessing import Pool
sys.path.append("../backend/")

import objects
import materials
import plotting
from stack import *

start = time.time()

################ Light parameters #####################
wl_1     = 500
wl_2     = 600
no_wl_1  = 4
# Set up light objects, starting with the wavelengths,
wavelengths = np.linspace(wl_1, wl_2, no_wl_1)
# and also specifying angles of incidence and refractive medium of semi-infinite
# layer that the light is incident upon (default value is n_inc = 1.0).
# Fields in homogeneous layers are expressed in a Fourier series of diffraction
# orders,where all orders within a radius of max_order_PWs in k-space are included.
light_list  = [objects.Light(wl, max_order_PWs = 1, theta = 0.0, phi = 0.0, \
    n_inc=1.5) for wl in wavelengths]

# Our structure must have a period, even if this is artificially imposed
# on a homogeneous thin film. What's more,
# it is critical that the period be consistent throughout a simulation!
period = 300

# Define each layer of the structure.
superstrate = objects.ThinFilm(period, height_nm = 'semi_inf',
    material = materials.Material(1.5 + 0.0j))
substrate   = objects.ThinFilm(period, height_nm = 'semi_inf',
    material = materials.Material(3.0 + 0.0j))

def simulate_stack(light):
    ################ Evaluate each layer individually ##############
    sim_superstrate = superstrate.calc_modes(light)
    sim_substrate   = substrate.calc_modes(light)
    ###################### Evaluate structure ######################
    """ Now define full structure. Here order is critical and
        stack list MUST be ordered from bottom to top!
    """

    stack = Stack((sim_substrate, sim_superstrate))
    # Calculate scattering matrices of the stack (for all polarisations).
    stack.calc_scat(pol = 'TE') # Incident light has TE polarisation,
    # which only effects the net transmission etc, not the matrices.

    return stack

stacks_list = map(simulate_stack, light_list)
# Save full simo data to .npz file for safe keeping!
np.savez('Simo_results', stacks_list=stacks_list)

# Calculation of the modes and scattering matrices of each layer
# as well as the scattering matrices of the interfaces of the stack
# is complete.
# From here on we can print, plot or manipulate the results.

# Alternatively, you may wish to finish the simo file here,
# and be output into an interactive python instance were you
# have access to all simulation objects and results for further
# manipulation. In this case you run this file as
# $ python -i simo_010-single_interface.py
# In this session the docstrings of objects/classes/methods
# can be accessed by typing

# >>> from pydoc import help
# >>> help(objects.Light)

# where we have accessed the docstring of the Light class from objects.py


######################## Post Processing ########################
# We can retrieve the propagation constants (k_z) of each layer.
# Let's print the values at the short wavelength in the superstrate,
wl_num = 0
lay = 1
betas = stacks_list[wl_num].layers[lay].k_z
print 'k_z of superstrate \n', betas
# and save the values for the longest wavelength for the substrate.
wl_num = -1
lay = 0
betas = stacks_list[wl_num].layers[lay].k_z
np.savetxt('Substrate_k_zs.txt', betas.view(float).reshape(-1, 2))
# Note that saving to txt files is slower than saving data as .npz
# However txt files may be easily read by other programs...


# We can also access the scattering matrices of individual layers,
# and of interfaces of the stack.
# For instance the reflection scattering matrix off the top
# of the substrate when considered as an isolated layer.
wl_num = -1
lay = 0
R12_sub = stacks_list[wl_num].layers[lay].R12
print 'R12 of substrate \n', R12_sub

# The reflection matrix for the reflection off the top of the
# superstrate-substrate interface meanwhile is a property of the stack.
R_interface = stacks_list[wl_num].R_net
# Let us plot this matrix in greyscale.
plotting.vis_scat_mats(R_interface)
# Since all layers are homogeneous this matrix should only have non-zero
# entries on the diagonal.

# Lastly, we can also plot the transmission, reflection, absorption
# of each layer and of the stack as a whole.
plotting.t_r_a_plots(stacks_list)

# p.s. we'll keep an eye on the time...
######################## Wrapping up ########################
print '\n*******************************************'
# Calculate and record the (real) time taken for simulation,
elapsed = (time.time() - start)
hms     = str(datetime.timedelta(seconds=elapsed))
hms_string = 'Total time for simulation was \n \
    %(hms)s (%(elapsed)12.3f seconds)'% {
            'hms'       : hms,
            'elapsed'   : elapsed, }
print hms_string
print '*******************************************'
print ''

# and store this info.
python_log = open("python_log.log", "w")
python_log.write(hms_string)
python_log.close()

Dispersion & Parallel Computation

"""
Simulating an interface between 2 homogeneous, dispersive media.
We use multiple CPUs.
"""

import time
import datetime
import numpy as np
import sys
from multiprocessing import Pool
sys.path.append("../backend/")

import objects
import materials
import plotting
from stack import *

start = time.time()

# We begin by remove all results of previous simulations.
plotting.clear_previous()

################ Simulation parameters ################
# Select the number of CPUs to use in simulation.
num_cores = 2

################ Light parameters #####################
wl_1     = 400
wl_2     = 800
no_wl_1  = 4
# Set up light objects (no need to specifiy n_inc as light incident from
# Air with n_inc = 1.0).
wavelengths = np.linspace(wl_1, wl_2, no_wl_1)
light_list  = [objects.Light(wl, max_order_PWs = 1, theta = 0.0, phi = 0.0) \
    for wl in wavelengths]

# The period must be consistent throughout a simulation!
period = 300

# Define each layer of the structure, now with dispersive media.
# The refractive indices are interpolated from tabulated data.
superstrate = objects.ThinFilm(period, height_nm = 'semi_inf',
    material = materials.Air)
substrate   = objects.ThinFilm(period, height_nm = 'semi_inf',
    material = materials.SiO2_a) # Amorphous silica

def simulate_stack(light):
    ################ Evaluate each layer individually ##############
    sim_superstrate = superstrate.calc_modes(light)
    sim_substrate   = substrate.calc_modes(light)
    ###################### Evaluate structure ######################
    """ Now define full structure. Here order is critical and
        stack list MUST be ordered from bottom to top!
    """

    stack = Stack((sim_substrate, sim_superstrate))
    stack.calc_scat(pol = 'TM') # This time TM polarised light is incident.

    return stack

# Run wavelengths in parallel across num_cores CPUs using multiprocessing package.
pool = Pool(num_cores)
stacks_list = pool.map(simulate_stack, light_list)
# Save full simo data to .npz file for safe keeping!
np.savez('Simo_results', stacks_list=stacks_list)

######################## Post Processing ########################
# This time let's visualise the net Transmission scattering matrix,
# which describes the propagation of light all the way from the superstrate into
# the substrate. When studying diffractive layers it is useful to know how many
# of theplane waves of the substrate are propagating, so lets include this.
wl_num = -1
T_net = stacks_list[wl_num].T_net
nu_prop = stacks_list[wl_num].layers[0].num_prop_pw_per_pol
plotting.vis_scat_mats(T_net, nu_prop_PWs=nu_prop)

# Let's just plot the spectra and see the effect of changing refractive indices.
plotting.t_r_a_plots(stacks_list)

######################## Wrapping up ########################
print '\n*******************************************'
# Calculate and record the (real) time taken for simulation,
elapsed = (time.time() - start)
hms     = str(datetime.timedelta(seconds=elapsed))
hms_string = 'Total time for simulation was \n \
    %(hms)s (%(elapsed)12.3f seconds)'% {
            'hms'       : hms,
            'elapsed'   : elapsed, }
print hms_string
print '*******************************************'
print ''

# and store this info.
python_log = open("python_log.log", "w")
python_log.write(hms_string)
python_log.close()

Thin Film Stack

"""
Simulating a stack of homogeneous, dispersive media.
"""

import time
import datetime
import numpy as np
import sys
from multiprocessing import Pool
sys.path.append("../backend/")

import objects
import materials
import plotting
from stack import *

start = time.time()

# Remove results of previous simulations.
plotting.clear_previous()

################ Simulation parameters ################
# Select the number of CPUs to use in simulation.
num_cores = 2

################ Light parameters #####################
wl_1     = 400
wl_2     = 800
no_wl_1  = 4
wavelengths = np.linspace(wl_1, wl_2, no_wl_1)
light_list  = [objects.Light(wl, max_order_PWs = 1, theta = 0.0, phi = 0.0) \
for wl in wavelengths]

# The period must be consistent throughout a simulation!
period = 300

# Define each layer of the structure.
superstrate = objects.ThinFilm(period, height_nm = 'semi_inf',
    material = materials.Air)
# Define a thin film with (finite) thickness in nm and constant refractive index
TF_1 = objects.ThinFilm(period, height_nm = 100,
    material = materials.Material(1.0 + 0.05j))
# EMUstack calculation time is independent dispersion and thickness of layer!
# This layer is made of Indium Phosphide, the tabulated refractive index of which
# is stored in EMUstack/data/
# We artificially set the imaginary part of the layer to zero for all wavelengths.
TF_2 = objects.ThinFilm(period, height_nm = 5e6,
    material = materials.InP, loss=False)
# By default loss = True
TF_3 = objects.ThinFilm(period, height_nm = 52,
    material = materials.Si_a)
# Note that the semi-inf substrate must be lossess so that EMUstack can distinguish
# propagating plane waves that carry energy from evanescent waves which do not.
# This layer is therefore crystalline silicon with Im(n) == 0.
substrate   = objects.ThinFilm(period, height_nm = 'semi_inf',
    material = materials.Si_c, loss=False)

def simulate_stack(light):
    ################ Evaluate each layer individually ##############
    sim_superstrate = superstrate.calc_modes(light)
    sim_TF_1 = TF_1.calc_modes(light)
    sim_TF_2 = TF_2.calc_modes(light)
    sim_TF_3 = TF_3.calc_modes(light)
    sim_substrate   = substrate.calc_modes(light)
    ###################### Evaluate structure ######################
    """ Now define full structure. Here order is critical and
        stack list MUST be ordered from bottom to top!
    """
# We can now stack these layers of finite thickness however we wish.
    stack = Stack((sim_substrate, sim_TF_1, sim_TF_3, sim_TF_2, sim_TF_1, \
        sim_superstrate))
    stack.calc_scat(pol = 'TM')

    return stack

# Run wavelengths in parallel across num_cores CPUs using multiprocessing package.
pool = Pool(num_cores)
stacks_list = pool.map(simulate_stack, light_list)
# Save full simo data to .npz file for safe keeping!
np.savez('Simo_results', stacks_list=stacks_list)

######################## Post Processing ########################
# We will now see the absorption in each individual layer as well as of the stack.
plotting.t_r_a_plots(stacks_list)

######################## Wrapping up ########################
print '\n*******************************************'
# Calculate and record the (real) time taken for simulation,
elapsed = (time.time() - start)
hms     = str(datetime.timedelta(seconds=elapsed))
hms_string = 'Total time for simulation was \n \
    %(hms)s (%(elapsed)12.3f seconds)'% {
            'hms'       : hms,
            'elapsed'   : elapsed, }
print hms_string
print '*******************************************'
print ''

# and store this info.
python_log = open("python_log.log", "w")
python_log.write(hms_string)
python_log.close()

Including Metals

"""
EUMstack loves metal \m/
However, as we saw in the previous example the substrate layer must be lossless,
so that we can distinguish propagating waves from evanescent ones.
To terminate the stack with a metalic mirror we must make it finite, but very thick.
"""

import time
import datetime
import numpy as np
import sys
from multiprocessing import Pool
sys.path.append("../backend/")

import objects
import materials
import plotting
from stack import *

start = time.time()

# Remove results of previous simulations.
plotting.clear_previous()

################ Simulation parameters ################
# Select the number of CPUs to use in simulation.
num_cores = 2

################ Light parameters #####################
wl_1     = 400
wl_2     = 800
no_wl_1  = 4
wavelengths = np.linspace(wl_1, wl_2, no_wl_1)
light_list  = [objects.Light(wl, max_order_PWs = 1, theta = 0.0, phi = 0.0)\
    for wl in wavelengths]

# The period must be consistent throughout a simulation!
period = 300

# Define each layer of the structure, as in last example.
superstrate = objects.ThinFilm(period, height_nm = 'semi_inf',
    material = materials.Air)
TF_2 = objects.ThinFilm(period, height_nm = 5e6,
    material = materials.InP, loss=False)
TF_3 = objects.ThinFilm(period, height_nm = 52,
    material = materials.Si_a)
# Realistically a few micron thick mirror would do the trick,
# but EMUstack is height agnostic.... so what the hell.
mirror = objects.ThinFilm(period, height_nm = 1e5,
    material = materials.Ag)
substrate   = objects.ThinFilm(period, height_nm = 'semi_inf',
    material = materials.Air)

def simulate_stack(light):
    ################ Evaluate each layer individually ##############
    sim_superstrate = superstrate.calc_modes(light)
    sim_mirror = mirror.calc_modes(light)
    sim_TF_2 = TF_2.calc_modes(light)
    sim_TF_3 = TF_3.calc_modes(light)
    sim_substrate   = substrate.calc_modes(light)
    ###################### Evaluate structure ######################
    """ Now define full structure. Here order is critical and
        stack list MUST be ordered from bottom to top!
    """
# Put semi-inf substrate below thick mirror so that propagating energy is defined.
    stack = Stack((sim_substrate, sim_mirror, sim_TF_3, sim_TF_2, sim_superstrate))
    stack.calc_scat(pol = 'TM')

    return stack

pool = Pool(num_cores)
stacks_list = pool.map(simulate_stack, light_list)
# Save full simo data to .npz file for safe keeping!
np.savez('Simo_results', stacks_list=stacks_list)

######################## Post Processing ########################
# The total transmission should be zero.
plotting.t_r_a_plots(stacks_list)

######################## Wrapping up ########################
print '\n*******************************************'
# Calculate and record the (real) time taken for simulation,
elapsed = (time.time() - start)
hms     = str(datetime.timedelta(seconds=elapsed))
hms_string = 'Total time for simulation was \n \
    %(hms)s (%(elapsed)12.3f seconds)'% {
            'hms'       : hms,
            'elapsed'   : elapsed, }
print hms_string
print '*******************************************'
print ''

# and store this info.
python_log = open("python_log.log", "w")
python_log.write(hms_string)
python_log.close()

1D Grating

"""
Simulating a lamellar grating that is periodic in x only.
For this simulation EMUstack uses the 1D diffraction orders for the basis
of the plane waves and carries out a 1D FEM calculation for the modes of
the grating.
"""

import time
import datetime
import numpy as np
import sys
from multiprocessing import Pool
sys.path.append("../backend/")

import objects
import materials
import plotting
from stack import *

start = time.time()

# Remove results of previous simulations.
plotting.clear_previous()

################ Simulation parameters ################
# Select the number of CPUs to use in simulation.
num_cores = 2

################ Light parameters #####################
wl_1     = 400
wl_2     = 800
no_wl_1  = 2
wavelengths = np.linspace(wl_1, wl_2, no_wl_1)
light_list  = [objects.Light(wl, max_order_PWs = 5, theta = 0.0, phi = 0.0) for wl in wavelengths]

# The period must be consistent throughout a simulation!
period = 300

# Define each layer of the structure
# We need to inform EMUstack at this point that all layers in the stack will
# be at most be periodic in one dimension (i.e. there are no '2D_arrays's).
# This is done with the Keyword Arg 'world_1d' and all homogenous layers are
# calculated using the PW basis of 1D diffraction orders.
superstrate = objects.ThinFilm(period, height_nm = 'semi_inf', world_1d=True,
    material = materials.Air)

substrate   = objects.ThinFilm(period, height_nm = 'semi_inf', world_1d=True,
    material = materials.Air)
# Define 1D grating that is periodic in x.
# The mesh for this is always made 'live' in objects.py the number of
# FEM elements used is given by 1/lc_bkg.
# See Fortran Backends section of tutorial for more details.
grating = objects.NanoStruct('1D_array', period, int(round(0.75*period)), height_nm = 2900,
    background = materials.Material(1.46 + 0.0j), inclusion_a = materials.Material(5.0 + 0.0j),
    loss = True, lc_bkg = 0.0051)

def simulate_stack(light):
    ################ Evaluate each layer individually ##############
    sim_superstrate = superstrate.calc_modes(light)
    sim_grating     = grating.calc_modes(light)
    sim_substrate   = substrate.calc_modes(light)
    ###################### Evaluate structure ######################
    """ Now define full structure. Here order is critical and
        stack list MUST be ordered from bottom to top!
    """

    stack = Stack((sim_substrate, sim_grating, sim_superstrate))
    stack.calc_scat(pol = 'TE')

    return stack

pool = Pool(num_cores)
# stacks_list = pool.map(simulate_stack, light_list)
stacks_list = map(simulate_stack, light_list)
# Save full simo data to .npz file for safe keeping!
np.savez('Simo_results', stacks_list=stacks_list)

######################## Post Processing ########################
# The total transmission should be zero.
plotting.t_r_a_plots(stacks_list)

######################## Wrapping up ########################
print '\n*******************************************'
# Calculate and record the (real) time taken for simulation,
elapsed = (time.time() - start)
hms     = str(datetime.timedelta(seconds=elapsed))
hms_string = 'Total time for simulation was \n \
    %(hms)s (%(elapsed)12.3f seconds)'% {
            'hms'       : hms,
            'elapsed'   : elapsed, }
print hms_string
print '*******************************************'
print ''

# and store this info.
python_log = open("python_log.log", "w")
python_log.write(hms_string)
python_log.close()

2D Grating

"""
Simulating a nanowire array with period 600 nm and NW diameter 120 nm.
"""

import time
import datetime
import numpy as np
import sys
from multiprocessing import Pool
sys.path.append("../backend/")

import objects
import materials
import plotting
from stack import *

start = time.time()
################ Simulation parameters ################

# Number of CPUs to use in simulation
num_cores = 7

# Remove results of previous simulations
plotting.clear_previous()

################ Light parameters #####################
wl_1     = 310
wl_2     = 1127
no_wl_1  = 3
# Set up light objects
wavelengths = np.linspace(wl_1, wl_2, no_wl_1)
light_list  = [objects.Light(wl, max_order_PWs = 2, theta = 0.0, phi = 0.0) \
    for wl in wavelengths]

# Period must be consistent throughout simulation!!!
period = 600

# In this example we set the number of Bloch modes to use in the simulation
# Be default it is set to be slightly greater than the number of PWs.
num_BMs = 200

superstrate = objects.ThinFilm(period, height_nm = 'semi_inf',
    material = materials.Air, loss = False)

substrate  = objects.ThinFilm(period, height_nm = 'semi_inf',
    material = materials.SiO2, loss = False)

NW_diameter = 120
NW_array = objects.NanoStruct('2D_array', period, NW_diameter, height_nm = 2330,
    inclusion_a = materials.Si_c, background = materials.Air, loss = True,
    make_mesh_now = True, force_mesh = True, lc_bkg = 0.1, lc2= 2.0)
# Here we get EMUstack to make the FEM mesh automagically using our input parameters.
# the lc_bkg parameter sets the baseline distance between points on the FEM mesh,
# lc_bkg/lc2 is the distance between mesh points that lie on the inclusion boundary.
# There are higher lc parameters which are used when including multiple inclusions.

# Alternatively we can specify a pre-made mesh as follows.
NW_array2 = objects.NanoStruct('2D_array', period, NW_diameter, height_nm = 2330,
    inclusion_a = materials.Si_c, background = materials.Air, loss = True,
    make_mesh_now = False, mesh_file='4testing-600_120.mail')


def simulate_stack(light):

    ################ Evaluate each layer individually ##############
    sim_superstrate = superstrate.calc_modes(light)
    sim_substrate   = substrate.calc_modes(light)
    sim_NWs         = NW_array.calc_modes(light, num_BMs=num_BMs)

    ###################### Evaluate structure ######################
    """ Now define full structure. Here order is critical and
        stack list MUST be ordered from bottom to top!
    """

    stack = Stack((sim_substrate, sim_NWs, sim_superstrate))
    stack.calc_scat(pol = 'TE')

    return stack


# Run in parallel across wavelengths.
pool = Pool(num_cores)
stacks_list = pool.map(simulate_stack, light_list)
# Save full simo data to .npz file for safe keeping!
np.savez('Simo_results', stacks_list=stacks_list)


######################## Plotting ########################

# We here wish to know the photovoltaic performance of the structure,
# where all light absorbed in the NW layer is considered to produce exactly
# one electron-hole pair.
# To do this we specify which layer of the stack is the PV active layer
# (default active_layer_nu=1), and indicate that we want to calculate
# the ideal short circuit current (J_sc) of the cell.
# We could also calculate the 'ultimate efficiency' by setting ult_eta=True.
plotting.t_r_a_plots(stacks_list, active_layer_nu=1, J_sc=True)

# We also plot the dispersion relation for each layer.
plotting.omega_plot(stacks_list, wavelengths)

######################## Wrapping up ########################
# Calculate and record the (real) time taken for simulation
elapsed = (time.time() - start)
hms     = str(datetime.timedelta(seconds=elapsed))
hms_string = 'Total time for simulation was \n \
    %(hms)s (%(elapsed)12.3f seconds)'% {
            'hms'       : hms,
            'elapsed'   : elapsed, }

python_log = open("python_log.log", "w")
python_log.write(hms_string)
python_log.close()

print hms_string
print '*******************************************'
print ''

Angles of Incidence & Eliptical Inclusions

"""
Simulating circular dichroism effect in elliptic nano hole arrays
as in T Cao1 and Martin J Cryan doi:10.1088/2040-8978/14/8/085101.
"""

import time
import datetime
import numpy as np
import sys
from multiprocessing import Pool
sys.path.append("../backend/")

import objects
import materials
import plotting
from stack import *

start = time.time()
################ Simulation parameters ################

# Number of CPUs to use in simulation
num_cores = 4

# Remove results of previous simulations
plotting.clear_previous()

################ Light parameters #####################
wl_1     = 300
wl_2     = 1000
no_wl_1  = 21
# Set up light objects
wavelengths = np.linspace(wl_1, wl_2, no_wl_1)
light_list  = [objects.Light(wl, theta = 45, phi = 45, max_order_PWs = 2) \
    for wl in wavelengths]



# Period must be consistent throughout simulation!!!
period = 165
diam1  = 140
diam2  = 60
ellipticity = (float(diam1-diam2))/float(diam1)

# Replicating the geometry of the paper we set up a gold layer with elliptical air
# holes. To get good agreement with the published work we use the Drude model for Au.
# Note that better physical results are obtained using the tabulated data for Au!
Au_NHs = objects.NanoStruct('2D_array', period, diam1, inc_shape = 'ellipse',
    ellipticity = ellipticity, height_nm = 60,
    inclusion_a = materials.Air, background = materials.Au_drude, loss = True,
    make_mesh_now = True, force_mesh = True, lc_bkg = 0.2, lc2= 5.0)

superstrate = objects.ThinFilm(period = period, height_nm = 'semi_inf',
    material = materials.Air, loss = True)
substrate = objects.ThinFilm(period = period, height_nm = 'semi_inf',
    material = materials.Air, loss = False)

# Again for this example we fix the number of BMs.
num_BMs = 50

def simulate_stack(light):
    ################ Evaluate each layer individually ##############
    sim_superstrate = superstrate.calc_modes(light)
    sim_Au   = Au_NHs.calc_modes(light, num_BMs = num_BMs)
    sim_substrate = substrate.calc_modes(light)

    stackSub = Stack((sim_substrate, sim_Au, sim_superstrate))
    stackSub.calc_scat(pol = 'R Circ')
    stackSub2 = Stack((sim_substrate, sim_Au, sim_superstrate))
    stackSub2.calc_scat(pol = 'L Circ')
    saveStack = Stack((sim_substrate, sim_Au, sim_superstrate))

    a_CD = []
    t_CD = []
    r_CD = []
    for i in range(len(stackSub.a_list)):
        a_CD.append(stackSub.a_list.pop() - stackSub2.a_list.pop())
    for i in range(len(stackSub.t_list)):
        t_CD.append(stackSub.t_list.pop() - stackSub2.t_list.pop())
    for i in range(len(stackSub.r_list)):
        r_CD.append(stackSub.r_list.pop() - stackSub2.r_list.pop())
    saveStack.a_list = a_CD
    saveStack.t_list = t_CD
    saveStack.r_list = r_CD

    return saveStack


# Run in parallel across wavelengths.
pool = Pool(num_cores)
stacks_list = pool.map(simulate_stack, light_list)
# Save full simo data to .npz file for safe keeping!
np.savez('Simo_results', stacks_list=stacks_list)

######################## Plotting ########################
# Just to show how it's done we can add the height of the layer and some extra
# details to the file names and plot titles.
title = 'what_a_lovely_day-'

plotting.t_r_a_plots(stacks_list, add_height=Au_NHs.height_nm, add_name=title)


# Calculate and record the (real) time taken for simulation
elapsed = (time.time() - start)
hms     = str(datetime.timedelta(seconds=elapsed))
hms_string = 'Total time for simulation was \n \
    %(hms)s (%(elapsed)12.3f seconds)'% {
            'hms'       : hms,
            'elapsed'   : elapsed, }

python_log = open("python_log.log", "w")
python_log.write(hms_string)
python_log.close()

print '*******************************************'
print hms_string
print '*******************************************'
print ''

Plotting Fields 1D

"""
Show how to plot electric fields.
"""

import time
import datetime
import numpy as np
import sys
from multiprocessing import Pool
sys.path.append("../backend/")

import objects
import materials
import plotting
from stack import *

start = time.time()
################ Simulation parameters ################

# Number of CPUs to use in simulation
num_cores = 7

# Remove results of previous simulations
plotting.clear_previous()

################ Light parameters #####################
wl     = 615
light_list  = [objects.Light(wl, max_order_PWs = 10, theta = 0.0, phi = 0.0)]

# Period must be consistent throughout simulation!!!
period = 600

superstrate = objects.ThinFilm(period, height_nm = 'semi_inf', world_1d = True,
    material = materials.Air, loss = False)

substrate  = objects.ThinFilm(period, height_nm = 'semi_inf', world_1d = True,
    material = materials.Air, loss = False)

spacer  = objects.ThinFilm(period, height_nm = 200, world_1d = True,
    material = materials.SiO2_a, loss = True)

grating = objects.NanoStruct('1D_array', period, int(round(0.7*period)), height_nm = 400,
    background = materials.Material(1.45 + 0.0j),
    inclusion_a = materials.Material(3.77 + 0.01j),
    loss = True, lc_bkg = 0.005, plotting_fields = True)


def simulate_stack(light):

    ################ Evaluate each layer individually ##############
    sim_superstrate = superstrate.calc_modes(light)
    sim_substrate   = substrate.calc_modes(light)
    sim_grating     = grating.calc_modes(light)
    sim_spacer      = spacer.calc_modes(light)

    ###################### Evaluate structure ######################
    """ Now define full structure. Here order is critical and
        stack list MUST be ordered from bottom to top!
    """

    stack = Stack((sim_substrate, sim_spacer, sim_grating, sim_superstrate))
    stack.calc_scat(pol = 'TE')

    return stack


# Run in parallel across wavelengths.
pool = Pool(num_cores)
stacks_list = pool.map(simulate_stack, light_list)
# Save full simo data to .npz file for safe keeping!
np.savez('Simo_results', stacks_list = stacks_list)


######################## Plotting ########################

# Plot fields on slices through stack.
#
# Note that all field plots of previous simulations are deleted! Move any
# results that you wish to keep into a different folder, ideally copying the
# whole simo directory to future reference to simo parameters.
#
# plotting.fields_vertically(stacks_list)
# # We can also plot only the scattered field (disregarding the incident field)
# plotting.fields_vertically(stacks_list, no_incoming = True, add_name = '-no_incoming')
#
# The above fields are the total fields, we can also look at the fields of
# each individual Bloch mode, which for a 1D array is done like so,
plotting.Bloch_fields_1d(stacks_list)

######################## Wrapping up ########################
# Calculate and record the (real) time taken for simulation
elapsed = (time.time() - start)
hms     = str(datetime.timedelta(seconds=elapsed))
hms_string = 'Total time for simulation was \n \
    %(hms)s (%(elapsed)12.3f seconds)'% {
            'hms'       : hms,
            'elapsed'   : elapsed, }

python_log = open("python_log.log", "w")
python_log.write(hms_string)
python_log.close()

print hms_string
print '*******************************************'
print ''

Plotting Fields 2D

"""
Show how to plot electric fields.
"""

import time
import datetime
import numpy as np
import sys
from multiprocessing import Pool
sys.path.append("../backend/")

import objects
import materials
import plotting
from stack import *

start = time.time()
################ Simulation parameters ################

# Number of CPUs to use in simulation
num_cores = 7

# Remove results of previous simulations
plotting.clear_previous()

################ Light parameters #####################
wl     = 615
light_list  = [objects.Light(wl, max_order_PWs = 15, theta = 0.0, phi = 0.0)]

# Period must be consistent throughout simulation!!!
period = 600

superstrate = objects.ThinFilm(period, height_nm = 'semi_inf',
    material = materials.Air, loss = False)

substrate  = objects.ThinFilm(period, height_nm = 'semi_inf',
    material = materials.Air, loss = False)

spacer  = objects.ThinFilm(period, height_nm = 200,
    material = materials.SiO2_a, loss = True)

NW_diameter = 120
NW_array = objects.NanoStruct('2D_array', period, NW_diameter,
    height_nm = 2330, inclusion_a = materials.Si_c, background = materials.Air,
    loss = True, make_mesh_now = True, force_mesh = True, lc_bkg = 0.1,
    lc2= 2.0, plotting_fields = True)


def simulate_stack(light):

    ################ Evaluate each layer individually ##############
    sim_superstrate = superstrate.calc_modes(light)
    sim_substrate   = substrate.calc_modes(light)
    sim_NWs         = NW_array.calc_modes(light)
    sim_spacer      = spacer.calc_modes(light)

    ###################### Evaluate structure ######################
    """ Now define full structure. Here order is critical and
        stack list MUST be ordered from bottom to top!
    """

    stack = Stack((sim_substrate, sim_spacer, sim_NWs, sim_superstrate))
    stack.calc_scat(pol = 'TE')

    return stack


# Run in parallel across wavelengths.
pool = Pool(num_cores)
stacks_list = pool.map(simulate_stack, light_list)
# Save full simo data to .npz file for safe keeping!
np.savez('Simo_results', stacks_list = stacks_list)


######################## Plotting ########################

# Plot fields on slices through stack along the x & y axis,
# and along the diagonals.
# This is done through all layers of the stack and saved as png files.
#
# Note that all field plots of previous simulations are deleted! Move any
# results that you wish to keep into a different folder, ideally copying the
# whole simo directory to future reference to simo parameters.
#
plotting.fields_vertically(stacks_list)

# Plot fields in the x-y plane at a list of specified heights.
plotting.fields_in_plane(stacks_list, lay_interest = 2, z_values = [0.0, 2.0])
plotting.fields_in_plane(stacks_list, lay_interest = 1, z_values = [1.0, 3.2])

# Plot fields inside nanostructures in 3D which are viewed using gmsh.
plotting.fields_3d(stacks_list, lay_interest = 2)

# Save electric field values (all components) at a list of selected point.
plotting.field_values(stacks_list, lay_interest = 0, xyz_values = [(4.0, 2.5, 7.0), (1.0, 1.5, 3.0)])

######################## Wrapping up ########################
# Calculate and record the (real) time taken for simulation
elapsed = (time.time() - start)
hms     = str(datetime.timedelta(seconds=elapsed))
hms_string = 'Total time for simulation was \n \
    %(hms)s (%(elapsed)12.3f seconds)'% {
            'hms'       : hms,
            'elapsed'   : elapsed, }

python_log = open("python_log.log", "w")
python_log.write(hms_string)
python_log.close()

print hms_string
print '*******************************************'
print ''

Plotting Amplitudes

"""
Here we investigate how efficiently a stack of 1D gratings excite diffraction orders.
"""

import time
import datetime
import numpy as np
import sys
from multiprocessing import Pool
sys.path.append("../backend/")

import objects
import materials
import plotting
from stack import *

start = time.time()
################ Simulation parameters ################

# Number of CPUs to use in simulation
num_cores = 5

# Remove results of previous simulations
plotting.clear_previous()
################ Light parameters #####################
wavelengths = np.linspace(1500,1600,10)
light_list  = [objects.Light(wl, max_order_PWs = 6, theta = 0.0, phi = 0.0) \
    for wl in wavelengths]


################ Grating parameters #####################
# The period must be consistent throughout a simulation!
period = 700

superstrate = objects.ThinFilm(period, height_nm = 'semi_inf', world_1d = True,
    material = materials.Air, loss = False)

substrate  = objects.ThinFilm(period, height_nm = 'semi_inf', world_1d = True,
    material = materials.Air, loss = False)

absorber    = objects.ThinFilm(period, height_nm = 10, world_1d = True,
    material = materials.Material(1.0 + 0.05j), loss = True)

grating_1 = objects.NanoStruct('1D_array', period, int(round(0.75*period)),
    height_nm = 2900, background = materials.Material(1.46 + 0.0j),
    inclusion_a = materials.Material(3.61 + 0.0j), loss = True,
    lc_bkg = 0.005)


def simulate_stack(light):

    ################ Evaluate each layer individually ##############
    sim_superstrate = superstrate.calc_modes(light)
    sim_substrate   = substrate.calc_modes(light)
    sim_absorber    = absorber.calc_modes(light)
    sim_grating_1   = grating_1.calc_modes(light)

    ###################### Evaluate structure ######################
    """ Now define full structure. Here order is critical and
        stack list MUST be ordered from bottom to top!
    """

    stack = Stack((sim_substrate, sim_absorber, sim_grating_1, sim_superstrate))
    stack.calc_scat(pol = 'TE')

    return stack


# Run in parallel across wavelengths.
pool = Pool(num_cores)
stacks_list = pool.map(simulate_stack, light_list)
# Save full simo data to .npz file for safe keeping!
np.savez('Simo_results', stacks_list = stacks_list)

######################## Post Processing ########################
# We can plot the amplitudes of each transmitted plane wave order as a
# function of angle.
plotting.PW_amplitudes(stacks_list, add_name = '-default_substrate')
# By default this will plot the amplitudes in the substrate, however we can also give
# the index in the stack of a different homogeneous layer and calculate them here.
# We here chose a subset of orders to plot.
plotting.PW_amplitudes(stacks_list, chosen_PWs = [-1,0,2], \
    lay_interest = 1)

# When many plane wave orders are included these last plots can become confusing,
# so instead one may wish to sum together the amplitudes of all propagating orders,
# of all evanescent orders, and all far-evanescent orders
# (which have in plane k>n_H * k0).
plotting.evanescent_merit(stacks_list, lay_interest = 0)


plotting.BM_amplitudes(stacks_list, lay_interest = 2, chosen_BMs = [0,1,2,3,4,5])

# Lastly we also plot the transmission, reflection and absorption of each
# layer and the stack.
plotting.t_r_a_plots(stacks_list, xvalues = wavelengths)


######################## Wrapping up ########################
# Calculate and record the (real) time taken for simulation
elapsed = (time.time() - start)
hms     = str(datetime.timedelta(seconds=elapsed))
hms_string = 'Total time for simulation was \n \
    %(hms)s (%(elapsed)12.3f seconds)'% {
            'hms'       : hms,
            'elapsed'   : elapsed, }

python_log = open("python_log.log", "w")
python_log.write(hms_string)
python_log.close()

print hms_string
print '*******************************************'
print ''

Shear Transformations

"""
Here we introduce a shear transformation to shift layers relative to one
another in the plane.
"""

import time
import datetime
import numpy as np
import sys
from multiprocessing import Pool
sys.path.append("../backend/")

import objects
import materials
import plotting
from stack import *

start = time.time()
################ Simulation parameters ################

# Number of CPUs to use in simulation
num_cores = 5

# Remove results of previous simulations
plotting.clear_previous()

################ Light parameters #####################
azi_angles = np.linspace(0,20,5)
wl = 1600
light_list  = [objects.Light(wl, max_order_PWs = 2, theta = p, phi = 0.0) \
    for p in azi_angles]

################ Grating parameters #####################
period = 760

superstrate = objects.ThinFilm(period, height_nm = 'semi_inf',
    material = materials.Air, loss = False)

substrate  = objects.ThinFilm(period, height_nm = 'semi_inf',
    material = materials.Air, loss = False)

grating_1 = objects.NanoStruct('1D_array', period, small_d=period/2,
    diameter1=int(round(0.25*period)), diameter2=int(round(0.25*period)),
    height_nm = 150, inclusion_a = materials.Material(3.61 + 0.0j),
    inclusion_b = materials.Material(3.61 + 0.0j),
    background = materials.Material(1.46 + 0.0j),
    loss = True, make_mesh_now = True, force_mesh = False, lc_bkg = 0.1, lc2= 3.0)

grating_2 = objects.NanoStruct('1D_array', period, int(round(0.75*period)),
    height_nm = 2900, background = materials.Material(1.46 + 0.0j),
    inclusion_a = materials.Material(3.61 + 0.0j),
    loss = True, make_mesh_now = True, force_mesh = False, lc_bkg = 0.1, lc2= 3.0)

num_BMs = 60


def simulate_stack(light):
    ################ Evaluate each layer individually ##############
    sim_superstrate = superstrate.calc_modes(light)
    sim_substrate   = substrate.calc_modes(light)
    sim_grating_1   = grating_1.calc_modes(light, num_BMs = num_BMs)
    sim_grating_2   = grating_2.calc_modes(light, num_BMs = num_BMs)

    ###################### Evaluate structure ######################
    """ Now define full structure. Here order is critical and
        stack list MUST be ordered from bottom to top!
    """

    # Shear is relative to top layer (ie incident light) and in units of d.
    stack = Stack((sim_substrate, sim_grating_1, sim_grating_2, sim_superstrate), \
        shears = ([(0.1,0.0),(-0.3,0.1),(0.2,0.5)]) )
    stack.calc_scat(pol = 'TE')

    return stack


# Run in parallel across wavelengths.
pool = Pool(num_cores)
stacks_list = pool.map(simulate_stack, light_list)
# Save full simo data to .npz file for safe keeping!
np.savez('Simo_results', stacks_list=stacks_list)

plotting.t_r_a_plots(stacks_list)

######################## Wrapping up ########################
# Calculate and record the (real) time taken for simulation
elapsed = (time.time() - start)
hms     = str(datetime.timedelta(seconds=elapsed))
hms_string = 'Total time for simulation was \n \
    %(hms)s (%(elapsed)12.3f seconds)'% {
            'hms'       : hms,
            'elapsed'   : elapsed, }

python_log = open("python_log.log", "w")
python_log.write(hms_string)
python_log.close()

print hms_string
print '*******************************************'
print ''

Ultrathin Absorption Limit - Varying n

"""
Simulating an ultrathin film with a range of real and imaginary refractive
indices. Can we reach the theoretical limit of 0.5 absorption?
"""

import time
import datetime
import numpy as np
import sys
from multiprocessing import Pool
sys.path.append("../backend/")

import objects
import materials
import plotting
from stack import *

start = time.time()

# Remove results of previous simulations.
plotting.clear_previous()

################ Simulation parameters ################
# Select the number of CPUs to use in simulation.
num_cores = 8

################ Light parameters #####################
wl     = 700
light  = objects.Light(wl, max_order_PWs = 0, theta = 0.0, phi = 0.0)

# The period must be consistent throughout a simulation!
period = 660

# Define each layer of the structure.
superstrate = objects.ThinFilm(period, height_nm = 'semi_inf',
    material = materials.Air, world_1d=True)
substrate   = objects.ThinFilm(period, height_nm = 'semi_inf',
    material = materials.Air, loss=False, world_1d=True)

n_min = 1
n_max = 10
num_n_re = 51
num_n_im = num_n_re
Re_n = np.linspace(n_min,n_max,num_n_re)
Im_n = np.linspace(n_max,n_min,num_n_im)
# Having lists run this way will ease plotting, as matshow plots from top left

def simulate_stack(Re):
    ################ Evaluate each layer individually ##############
    sim_superstrate = superstrate.calc_modes(light)
    sim_substrate   = substrate.calc_modes(light)

    # Re_stack = []
    # for Re in Re_n:
    Im_stack = []
    for Im in Im_n:
        TF_1 = objects.ThinFilm(period, height_nm = 10,
            material = materials.Material(Re + Im*1j))
        sim_TF_1 = TF_1.calc_modes(light)

        stack = Stack((sim_substrate, sim_TF_1, sim_superstrate))
        stack.calc_scat(pol = 'TM')

        Im_stack.append(stack)
        # Re_stack.append(Im_stack)

    return Im_stack

# Run wavelengths in parallel across num_cores CPUs using multiprocessing package.
pool = Pool(num_cores)
stacks_list = pool.map(simulate_stack, Re_n)
# # Save full simo data to .npz file for safe keeping!
# np.savez('Simo_results', stacks_list=stacks_list)

######################## Post Processing ########################

abs_mat = np.zeros((num_n_im,num_n_re))
for i in range(num_n_re):
    for j in range(num_n_im):
        abs_mat[j,i] = stacks_list[i][j].a_list[-1]

# Now plot as a function of Real and Imaginary refractive index.
# Requires a bit of manipulation of axis...
import matplotlib
matplotlib.use('pdf')
import matplotlib.pyplot as plt
fig = plt.figure()
linesstrength = 3
font = 18
ax1 = fig.add_subplot(1,1,1)
mat = ax1.matshow(abs_mat,cmap=plt.cm.hot)
cbar1 = plt.colorbar(mat, extend='neither',alpha=1)
ax1.xaxis.set_ticks_position('bottom')
ax1.set_xticks(np.linspace(n_min,(num_n-1),n_max))
ax1.set_yticks(np.linspace(n_min,(num_n-1),n_max))
ax1.set_xticklabels([str(i) for i in np.linspace(n_min,n_max,n_max-n_min+1)])
ax1.set_yticklabels([str(i) for i in np.linspace(n_max,n_min,n_max-n_min+1)])
ax1.set_xlabel('Re(n)',fontsize=font)
ax1.set_ylabel('Im(n)',fontsize=font)
plt.title('Absorption of %(h)5.1f nm thick film @ wl = %(wl)5.1f'% \
    {'h' : stacks_list[0][0].heights_nm()[0],'wl' : wl})
plt.savefig('ultrathin_limit')



######################## Wrapping up ########################
print '\n*******************************************'
# Calculate and record the (real) time taken for simulation,
elapsed = (time.time() - start)
hms     = str(datetime.timedelta(seconds=elapsed))
hms_string = 'Total time for simulation was \n \
    %(hms)s (%(elapsed)12.3f seconds)'% {
            'hms'       : hms,
            'elapsed'   : elapsed, }
print hms_string
print '*******************************************'
print ''

# and store this info.
python_log = open("python_log.log", "w")
python_log.write(hms_string)
python_log.close()

Varying a Layer of a Stack

"""
Simulating solar cell efficiency of nanohole array as a function of
substrate refractive indices (keeping geometry fixed).
We also average over a range of thicknesses to remove sharp Fabry-Perot resonances.
"""

import time
import datetime
import numpy as np
import sys
from multiprocessing import Pool
sys.path.append("../backend/")

import objects
import materials
import plotting
from stack import *

start = time.time()
################ Simulation parameters ################

# Number of CPUs to use in simulation
num_cores = 4

# Remove results of previous simulations
plotting.clear_previous()

################ Light parameters #####################
wl_1     = 310
wl_2     = 1127
no_wl_1  = 3
# Set up light objects
wavelengths = np.linspace(wl_1, wl_2, no_wl_1)
light_list  = [objects.Light(wl, max_order_PWs = 2, theta = 0.0, phi = 0.0) \
    for wl in wavelengths]


# Period must be consistent throughout simulation!!!
period = 550

cover  = objects.ThinFilm(period = period, height_nm = 'semi_inf',
    material = materials.Air, loss = True)

sub_ns = np.linspace(1.0,4.0,100)

NW_diameter = 480
NWs = objects.NanoStruct('1D_array', period, NW_diameter, height_nm = 2330,
    inclusion_a = materials.Si_c, background = materials.Air, loss = True,
    make_mesh_now = True, force_mesh = True, lc_bkg = 0.17, lc2= 2.5)

def simulate_stack(light):
    ################ Evaluate each layer individually ##############
    sim_cover = cover.calc_modes(light)
    sim_NWs   = NWs.calc_modes(light)

    # Loop over substrates
    stack_list = []
    for s in sub_ns:
        sub = objects.ThinFilm(period = period, height_nm = 'semi_inf',
        material = materials.Material(s + 0.0j), loss = False)
        sim_sub = sub.calc_modes(light)

        # Loop over heights to wash out sharp FP resonances
        average_t = 0
        average_r = 0
        average_a = 0

        num_h = 21
        for h in np.linspace(2180,2480,num_h):
            stackSub = Stack((sim_sub, sim_NWs, sim_cover), heights_nm = ([h]))
            stackSub.calc_scat(pol = 'TE')
            average_t += stackSub.t_list[-1]/num_h
            average_r += stackSub.r_list[-1]/num_h
            average_a += stackSub.a_list[-1]/num_h
        stackSub.t_list[-1] = average_t
        stackSub.r_list[-1] = average_r
        stackSub.a_list[-1] = average_a
        stack_list.append(stackSub)

    return stack_list


# Run in parallel across wavelengths.
pool = Pool(num_cores)
stacks_list = pool.map(simulate_stack, light_list)
# Save full simo data to .npz file for safe keeping!
np.savez('Simo_results', stacks_list=stacks_list)



######################## Plotting ########################
eta = []
for s in range(len(sub_ns)):
    stack_label = s # Specify which stack you are dealing with.
    stack1_wl_list = []
    for i in range(len(wavelengths)):
        stack1_wl_list.append(stacks_list[i][stack_label])
    sub_n = sub_ns[s]
    Efficiency = plotting.t_r_a_plots(stack1_wl_list, ult_eta=True,
        stack_label=stack_label, add_name = str(s))
    eta.append(100.0*Efficiency[0])
    # Dispersion of structured layer is the same for all cases.
    if s == 0:
        plotting.omega_plot(stack1_wl_list, wavelengths, stack_label=stack_label)

# Now plot as a function of substrate refractive index.
import matplotlib
matplotlib.use('pdf')
import matplotlib.pyplot as plt
fig = plt.figure()
linesstrength = 3
font = 18
ax1 = fig.add_subplot(1,1,1)
ax1.plot(sub_ns,eta, 'k-o', linewidth=linesstrength)
ax1.set_xlabel('Substrate refractive index',fontsize=font)
ax1.set_ylabel(r'$\eta$ (%)',fontsize=font)
plt.savefig('eta_substrates')

# Animate spectra as a function of substrates.
from os import system as ossys
delay = 30 # delay between images in gif in hundredths of a second
names = 'Total_Spectra_stack'
gif_cmd = 'convert -delay %(d)i +dither -layers Optimize -colors 16 \
%(n)s*.pdf %(n)s.gif'% {
'd' : delay, 'n' : names}
ossys(gif_cmd)
opt_cmd = 'gifsicle -O2 %(n)s.gif -o %(n)s-opt.gif'% {'n' : names}
ossys(opt_cmd)
rm_cmd = 'rm %(n)s.gif'% {'n' : names}
ossys(rm_cmd)



# Calculate and record the (real) time taken for simulation
elapsed = (time.time() - start)
hms     = str(datetime.timedelta(seconds=elapsed))
hms_string = 'Total time for simulation was \n \
    %(hms)s (%(elapsed)12.3f seconds)'% {
            'hms'       : hms,
            'elapsed'   : elapsed, }

python_log = open("python_log.log", "w")
python_log.write(hms_string)
python_log.close()

print '*******************************************'
print hms_string
print '*******************************************'

Convergence Testing

"""
    Replicate Fig 2a from Handmer Opt Lett 2010
"""


import time
import datetime
import numpy as np
import sys
from multiprocessing import Pool
sys.path.append("../backend/")

import objects
import materials
import plotting
from stack import *

start = time.time()
################ Simulation parameters ################

# Number of CPUs to use in simulation
num_cores = 8

# Remove results of previous simulations
plotting.clear_previous()

################ Light parameters #####################
wavelengths = np.linspace(1600,900,1)

BMs = [11,27,59,99,163,227,299,395,507,635,755,883,1059,1227,1419]
B = 0

for PWs in np.linspace(1,10,10):
    light_list  = [objects.Light(wl, max_order_PWs = PWs, theta = 28.0, phi = 0.0) for wl in wavelengths]

    ################ Grating parameters #####################
    period = 760

    superstrate = objects.ThinFilm(period, height_nm = 'semi_inf',
        material = materials.Air, loss = False)

    substrate  = objects.ThinFilm(period, height_nm = 'semi_inf',
        material = materials.Air, loss = False)

    grating_1 = objects.NanoStruct('1D_array', period, small_d=period/2,
        diameter1=int(round(0.25*period)), diameter2=int(round(0.25*period)), height_nm = 150,
        inclusion_a = materials.Material(3.61 + 0.0j), inclusion_b = materials.Material(3.61 + 0.0j),
        background = materials.Material(1.46 + 0.0j),
        loss = True, make_mesh_now = True, force_mesh = True, lc_bkg = 0.1, lc2= 3.0)

    grating_2 = objects.NanoStruct('1D_array', period, int(round(0.75*period)), height_nm = 2900,
        background = materials.Material(1.46 + 0.0j), inclusion_a = materials.Material(3.61 + 0.0j),
        loss = True, make_mesh_now = True, force_mesh = True, lc_bkg = 0.1, lc2= 3.0)

    num_BMs = BMs[B]+30
    B += 1


    def simulate_stack(light):

        ################ Evaluate each layer individually ##############
        sim_superstrate = superstrate.calc_modes(light)
        sim_substrate   = substrate.calc_modes(light)
        sim_grating_1   = grating_1.calc_modes(light, num_BMs = num_BMs)
        sim_grating_2   = grating_2.calc_modes(light, num_BMs = num_BMs)

        ###################### Evaluate structure ######################
        """ Now define full structure. Here order is critical and
            stack list MUST be ordered from bottom to top!
        """

        stack = Stack((sim_substrate, sim_grating_1, sim_grating_2, sim_superstrate))
        # stack = Stack((sim_substrate, sim_grating_2, sim_superstrate))
        stack.calc_scat(pol = 'TE')
        return stack



    # Run in parallel across wavelengths.
    pool = Pool(num_cores)
    stacks_list = pool.map(simulate_stack, light_list)
    # Save full simo data to .npz file for safe keeping!
    np.savez('Simo_results', stacks_list=stacks_list)


    additional_name = str(int(PWs))
    plotting.t_r_a_plots(stacks_list, add_name = additional_name)



######################## Wrapping up ########################
# Calculate and record the (real) time taken for simulation
elapsed = (time.time() - start)
hms     = str(datetime.timedelta(seconds=elapsed))
hms_string = 'Total time for simulation was \n \
    %(hms)s (%(elapsed)12.3f seconds)'% {
            'hms'       : hms,
            'elapsed'   : elapsed, }

# python_log = open("python_log.log", "w")
# python_log.write(hms_string)
# python_log.close()

print hms_string
print '*******************************************'
print ''

Extraordinary Optical Transmission

"""
Simulating Extraordinary Optical Transmission
as in H. Liu, P. Lalanne, Nature 452 2008 doi:10.1038/nature06762
"""

import time
import datetime
import numpy as np
import sys
from multiprocessing import Pool
sys.path.append("../backend/")

import objects
import materials
import plotting
from stack import *

start = time.time()
################ Simulation parameters ################

# Number of CPUs to use in simulation
num_cores = 16

# Remove results of previous simulations
plotting.clear_previous()

################ Light parameters #####################
wl_1     = 0.85*940
wl_2     = 1.15*940
no_wl_1  = 600
# Set up light objects
wavelengths = np.linspace(wl_1, wl_2, no_wl_1)
# wavelengths = np.array([785,788,790,792,795])
light_list  = [objects.Light(wl, max_order_PWs = 5, theta = 0.0, phi = 0.0) for wl in wavelengths]


#period must be consistent throughout simulation!!!
period = 940
diam1 = 266
NHs = objects.NanoStruct('2D_array', period, diam1, height_nm = 200,
    inclusion_a = materials.Air, background = materials.Au, loss = True,
    inc_shape = 'square',
    make_mesh_now = True, force_mesh = True, lc_bkg = 0.12, lc2= 5.0, lc3= 3.0)#lc_bkg = 0.08, lc2= 5.0)

strate  = objects.ThinFilm(period = period, height_nm = 'semi_inf',
    material = materials.Air, loss = False)

NH_heights = [200]
# num_h = 21
# NH_heights = np.linspace(50,3000,num_h)

def simulate_stack(light):
    ################ Evaluate each layer individually ##############
    sim_NHs    = NHs.calc_modes(light)
    sim_strate = strate.calc_modes(light)

# Loop over heights
    height_list = []
    for h in NH_heights:
        stackSub = Stack((sim_strate, sim_NHs, sim_strate), heights_nm = ([h]))
        stackSub.calc_scat(pol = 'TE')
        height_list.append(stackSub)

    return [height_list]


# Run in parallel across wavelengths.
pool = Pool(num_cores)
stacks_list = pool.map(simulate_stack, light_list)
# Save full simo data to .npz file for safe keeping!
np.savez('Simo_results', stacks_list=stacks_list)



######################## Plotting ########################
last_light_object = light_list.pop()

wls_normed = wavelengths/period

for h in range(len(NH_heights)):
    height = NH_heights[h]
    wl_list = []
    stack_label = 0
    for wl in range(len(wavelengths)):
        wl_list.append(stacks_list[wl][stack_label][h])
    mess_name = '_h%(h)i'% {'h'   : h, }
    plotting.EOT_plot(wl_list, wls_normed, add_name = mess_name, savetxt = True)
# Dispersion
plotting.omega_plot(wl_list, wavelengths)


# Calculate and record the (real) time taken for simulation
elapsed = (time.time() - start)
hms     = str(datetime.timedelta(seconds=elapsed))
hms_string = 'Total time for simulation was \n \
    %(hms)s (%(elapsed)12.3f seconds)'% {
            'hms'       : hms,
            'elapsed'   : elapsed, }

python_log = open("python_log.log", "w")
python_log.write(hms_string)
python_log.close()

print '*******************************************'
print hms_string
print '*******************************************'
print ''

Screen Sessions

screen

is an extremely useful little linux command. In the context of long-ish calculations it has two important applications; ensuring your calculation is unaffected if your connection to a remote machine breaks, and terminating calculations that have hung without closing the terminal. For more information see the manual:

$ man screen

or see online discussions here, and here.

The screen session or also called screen instance looks just like your regular terminal/putty, but you can disconnect from it (close putty, turn off your computer etc.) and later reconnect to the screen session and everything inside of this will have kept running. You can also reconnect to the session from a different computer via ssh.

Basic Usage

To install screen:

$ sudo apt-get install screen

To open a new screen session:

$ screen

We can start a new calculation here:

$ cd EMUstack/examples/
$ python simo_040-2D_array.py

We can then detach from the session (leaving everything in the screen running) by typing:

Ctrl +a
Ctrl +d

We can now monitor the processes in that session:

$ top

Where we note the numerous running python processes that EMUstack has started. Watching the number of processes is useful for checking if a long simulation is near completion (which is indicated by the number of processes dropping to less than the specified num_cores).

We could now start another screen and run some more calculations in this terminal (or do anything else). If we want to access the first session we ‘reattach’ by typing:

Ctrl +a +r

Or entering the following into the terminal:

$ screen -r

If there are multiple sessions use:

$ screen -ls

to get a listing of the sessions and their ID numbers. To reattach to a particular screen, with ID 1221:

$ screen -r 1221

To terminate a screen from within type:

Ctrl+d

Or, taking the session ID from the previous example:

screen -X -S 1221 kill

Terminating EMU stacks

If (for some estranged reason) a simulation hangs, we can kill all python instances upon the machine:

$ pkill python

If a calculation hangs from within a screen session one must first detach from that session then kill python. A more targeted way to kill processes is using their PID:

$ kill PID

Or if this does not suffice be a little more forceful:

$ kill -9 PID

The PID is found from one of two ways:

$ top
$ ps -fe | grep username