3.1 Creating a simple geometry

To use PSlab you must first define the structure for simulation. In order to do this you need to load the module Geometry. Start Python interpreter (usually by entering command python in UN*X shell) and type


 >>> from pslab import Geometry
 PSlab, version 0.1.1, Copyright (C) 2005 by Maciej Dems and Tomasz Czyszanowski
 Developed at Laboratory of Computer Physics, Technical University of Lodz
 
 Please note that PSlab comes with ABSOLUTELY NO WARRANTY. This is free
 software, and you are welcome to redistribute it under the conditions
 of GNU General Public Licence version 2, or any higher version.
 
 >>>
 
The appearance of the banner means that program has loaded successfully. Now you can create the Geometry3D object. Its constructor take two parameters, which are either the two-element tuples of the lattice vectors components3.1, or the $x$ and $y$ dimensions of the orthogonal (rectangular) unit cell. In the latter case the lattice vectors are assumed to be $[x,0]$ and $[0,y]$. Because we use the plane-wave expansion here, the boundary condtitions in $xy$ plane are always periodic.

We will run computation for the rectangular cell with all the distances expressed in its lattice constant, so we type

>>> geometry = Geometry.Geometry3D(1, 1)
>>>

Now, having geometry defined, let's create some layers. Assume we want to run a simulation of photonic crystal slab having refractive index equal to 3.5 with air as the cladding. As the structure have an inverse symmetry in the plane paralel to the layers, we need to define only two layers - half of the slab core and cladding.

>>> core = Geometry.Layer(3.5**2)
>>> cladding = Geometry.Layer(epsilon=1)
>>>
The Layer class has the constructor which require at least one parameter - the permittivity of the bulk of the layer. You can also specify permeability of the layer giving it as a second parameter. When you skip it the default value 1 is assumed. If you wish you can give the parameters with thier names as in the following example, in which case the order does not need to be preserved.

As every layer must have the finite thickness (specified at the later stage) we want to add also the absorbing PML as the boundary conditions. For this reason we must specify both $\epsilon$ and $\mu$ and what more they must have form of the diagonal tensor. With PSlab this is possible - you simple give the diagonal tensor components in the list in order $[x,y,z]$

>>> sz = 1-2j
>>> PML = Geometry.Layer(epsilon=[sz,sz,1/sz], mu=[sz,sz,1/sz])
>>>

To make our structure a photonic crystal let's add an air cylinder to the core

>>> core.addObject( Geometry.Cylinder(center=(0,0), radius=0.3, epsilon=1) )
>>>
The center parameter is always a tuple with two-dimensional coordinates.

Now we are ready to arrange our layers in a stack. We do this by using method addLayer of our geometry object (mind the order of the layers -- for symmetric position the symmetry plane is before the first layer)

>>> geometry.addLayer( core(0.3) )
>>> geometry.addLayer( cladding(4.0) )
>>> geometry.addLayer( PML(width=0.5) )
>>>

Having our structure defined lets run some simulation. For this purpose we must use the class Simulation from the pslab module

>>> from pslab import Simulation
>>> simulation = Simulation(geometry, size=3)
Creating Plane Wave grid...
  number of distinct layers: 3
  number of inverse vectors: [7,7], total matrix size : 98 (49 planewaves)
  whole mesh is [104x104]

Setting material coefficients for layer 0...

Setting material coefficients for layer 1...
  layer has diagonal Q-matrix :)

Setting material coefficients for layer 2...
  layer has diagonal Q-matrix :)

Creating simple diagonalizer...

Creating admittance-matrix solver...
  the stack consists of 3 layers, interface is after 0 layers
  the structure has symmetry in z direction

>>>
where the size parameter decides about the number of planewaves used. For three-dimenstional simulation $N = (2 \cdot size + 1)^2$. Now we can for example find the eigenmode around some frequency for the default wavevector (equal to zero).

>>> simulation.getMode(2.5)

Searching directly for the mode starting from (2.5+0j)
  searching for the mode with Broyden method starting from (2.5+0j)...
  found mode at (2.53406515-6.0354957e-16j)

>>>
It must be noted that in whole program frequencies are given as the $k_0 = \omega/c$ where $\omega$ is the angular frequency and $c$ the speed of light.

The short example presented here can be found in examples/example1.py. Below there is the full file.

#!/usr/bin/python

## Import the program modules
from pslab import Geometry, Simulation

## Declare the geometry. The unit cell is a 1x1 square
geometry = Geometry.Geometry3D(1, 1)

## Fill it with some layers

# PML has anisotropic electric and magnetic tensor
sz = 1-2j
PML = Geometry.Layer(epsilon=[sz,sz,1/sz], mu=[sz,sz,1/sz])

# Cladding is a uniform layer of air
cladding = Geometry.Layer(1)

# Core has isotropic epsilon equal to 12.25 and a circular air rod
core = Geometry.Layer(12.25)
core.addObject( Geometry.Cylinder(center=(0,0), radius=0.3, epsilon=1) )

## Now construct a stack
geometry.addLayer( core(0.3) )
geometry.addLayer( cladding(4.0) )
geometry.addLayer( PML(width=0.5) )

## Define the simulation
simulation = Simulation(geometry, size=3)

## Find some mode around angular frequency 0.5 c/a
mode = simulation.getMode(2.5)

print mode



Footnotes

... components3.1
The computational cell is a parallelogram with the edges defined by the lattice vectors.