doc/user-tutorial.html

<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 3.2//EN">
<HTML><HEAD>
<TITLE>User Tutorial</TITLE>
<LINK rel="Contents" href="index.html">
<LINK rel="Copyright" href="license.html">
<LINK rel="Start" href="index.html">
</HEAD>
<BODY TEXT="#000000" BGCOLOR="#FFFFFF">

Go to the <a href="analysis-tutorial.html">next</a>, <a href="installation.html">previous</a>, or <a href="index.html">main</a> section.
<hr>

<h1>User Tutorial</h1>

<p>In this section, we'll walk through the process of computing the
band structure and outputting some fields for a two-dimensional
photonic crystal using the MIT Photonic-Bands package.  This should
give you the basic idea of how it works and some of the things that
are possible.  Here, we tell you the truth, but not the whole truth.
The <a href="user-ref.html">User Reference</a> section gives a more
complete, but less colloquial, description of the features supported
by Photonic-Bands.  See also the <a href="analysis-tutorial.html">data
analysis tutorial</a> in the next section for more examples, focused
on analyzing and visualizing the results of MPB calculations.

<h2><a name="ctl">The ctl File</a></h2>

<p>The use of the Photonic-Bands package revolves around the control
file, abbreviated "ctl" and typically called something like
<code>foo.ctl</code> (although you can use any filename extension you
wish).  The ctl file specifies the geometry you wish to study, the
number of eigenvectors to compute, what to output, and everything else
specific to your calculation.  Rather than a flat, inflexible file
format, however, the ctl file is actually written in a scripting
language.  This means that it can be everything from a simple sequence
of commands setting the geometry, etcetera, to a full-fledged program
with user input, loops, and anything else that you might need.

<p>Don't worry, though--simple things are simple (you don't need to be
a <a
href="http://www.tuxedo.org/~esr/jargon/html/entry/Real-Programmer.html">Real
Programmer</a>), and even there you will appreciate the flexibility
that a scripting language gives you.  (e.g. you can input things in
any order, without regard for whitespace, insert comments where you
please, omit things when reasonable defaults are available...)

<p>The ctl file is actually implemented on top of the libctl library,
a set of utilities that are in turn built on top of the Scheme
language.  Thus, there are three sources of possible commands and
syntax for a ctl file:

<ul>

<li>Scheme, a powerful and beautiful programming language
developed at MIT, which has a particularly simple syntax: all
statements are of the form <code>(<i>function
arguments...</i>)</code>.  We run Scheme under the GNU Guile
interpreter (designed to be plugged into programs as a scripting and
extension language).  You don't need to learn much Scheme for a basic
ctl file, but it is always there if you need it; you can learn more
about it from these <a
href="http://ab-initio.mit.edu/libctl/doc/guile-links.html">Guile and
Scheme links</a>.

<p><li>libctl, a library that we built on top of Guile to simplify
communication between Scheme and scientific computation software.
libctl sets the basic tone of the user interface and defines a number
of useful functions.  See the <a
href="http://ab-initio.mit.edu/libctl/">libctl home page</a>.

<p><li>The MIT Photonic-Bands package, which defines all the interface
features that are specific to photonic band structure calculations.
This manual is primarily focused on documenting these features.

</ul>

<p>It would be an excellent idea at this point for you to go read the
<a href="http://ab-initio.mit.edu/libctl/doc/">libctl manual</a>,
particularly the <a
href="http://ab-initio.mit.edu/libctl/doc/basic-user.html">Basic User
Experience</a> section, which will give you an overview of what the
user interface is like, provide a crash course in the Scheme features
that are most useful here, and describe some useful general features.
We're not going to repeat this material (much), so learn it now!

<p>Okay, let's continue with our tutorial.  The MIT Photonic-Bands
(MPB) program is normally invoked by running something like:

<pre>
unix% mpb <i>foo.ctl</i> &gt;&amp; <i>foo.out</i>
</pre>

<p>which reads the ctl file <code>foo.ctl</code> and executes it,
saving the output to the file <code>foo.out</code>.  (Some sample ctl
files are provided for you in the <code>mpb-ctl/examples/</code>
directory.)  However, as you already know (since you obediently read
the libctl manual, right?), if you invoke <code>mpb</code> with no
arguments, you are dropped into an <em>interactive</em> mode in which
you can type commands and see their results immediately.  Why don't
you do that right now, in your terminal window?  Then, you can paste
in the commands from the tutorial as you follow it and see what they
do.  Isn't this fun?

<h2><a name="sq-rods">Our First Band Structure</a></h2>

<p>As our beginning example, we'll compute the band structure of a
two-dimensional square lattice of dielectric rods in air.  In our
control file, we'll first specify the parameters and geometry of the
simulation, and then tell it to run and give us the output.

<p>All of the parameters (each of which corresponds to a Scheme
variable) have default setting, so we only need to specify the ones we
need to change.  (For a complete listing of the parameter variables
and their current values, along with some other info, type
<code>(help)</code> at the <code>guile&gt;</code> prompt.)  One of the
parameters, <code>num-bands</code>, controls how many bands
(eigenstates) are computed at each k point.  If you type
<code>num-bands</code> at the prompt, it will return the current
value, <code>1</code>--this is too small; let's set it to a larger
value:

<pre>
(set! num-bands 8)
</pre>

<p>This is how we change the value of variables in Scheme (if you type
<code>num-bands</code> now, it will return <code>8</code>).  The next
thing we want to set (although the order really doesn't matter), is
the set of k points (Bloch wavevectors) we want to compute the bands
at.  This is controlled by the variable <code>k-points</code>, a list
of 3-vectors (which is initially empty).  We'll set it to the corners
of the irreducible Brillouin zone of a square lattice, Gamma, X, M,
and Gamma again:

<pre>
(set! k-points (list (vector3 0 0 0)     ; Gamma
                     (vector3 0.5 0 0)   ; X
                     (vector3 0.5 0.5 0) ; M
                     (vector3 0 0 0)))   ; Gamma
</pre>

<p>Notice how we construct a list, and how we make 3-vectors; notice
also how we can break things into multiple lines if we want, and that
a semicolon ('<code>;</code>') marks the start of a comment.
Typically, we'll want to also compute the bands at a lot of
intermediate k points, so that we see the continuous band structure.
Instead of manually specifying all of these intermediate points,
however, we can just call one of the functions provided by libctl to
interpolate them for us:

<pre>
(set! k-points (interpolate 4 k-points))
</pre>

<p>This takes the <code>k-points</code> and linearly interpolates four
new points between each pair of consecutive points; if we type
<code>k-points</code> now at the prompt, it will show us all 16 points
in the new list:

<p><code>(#(0 0 0) #(0.1 0.0 0.0) #(0.2 0.0 0.0) #(0.3 0.0 0.0) #(0.4
0.0 0.0) #(0.5 0 0) #(0.5 0.1 0.0) #(0.5 0.2 0.0) #(0.5 0.3 0.0) #(0.5
0.4 0.0) #(0.5 0.5 0) #(0.4 0.4 0.0) #(0.3 0.3 0.0) #(0.2 0.2 0.0)
#(0.1 0.1 0.0) #(0 0 0))</code>

<p>As is <a href="#units">described below</a>, all spatial vectors in
the program are specified in the basis of the lattice directions
normalized to <code>basis-size</code> lengths (default is
unit-normalized), while the k points are specified in the basis of the
(unnormalized) reciprocal lattice vectors.  In this case, we don't
have to specify the lattice directions, because we are happy with the
defaults--the lattice vectors default to the cartesian unit axes
(i.e. the most common case, a square/cubic lattice).  (The reciprocal
lattice vectors in this case are also the unit axes.)  We'll see how
to change the lattice vectors in later subsections.

<p>Now, we want to set the geometry of the system--we need to specify
which objects are in the primitive cell of the lattice, centered on
the origin.  This is controlled by the variable <code>geometry</code>,
which is a list of geometric objects.  As you know from reading the
libctl documentation, objects (more complicated, structured data
types), are created by statements of the form <code>(make <i>type</i>
(<i>property1 value1</i>) (<i>property2 value2</i>) ...)</code>.
There are various kinds (sub-classes) of geometric object: cylinders,
spheres, blocks, ellipsoids, and perhaps others in the future.  Right
now, we want a square lattice of rods, so we put a single dielectric
cylinder at the origin:

<pre>
(set! geometry (list (make cylinder 
                       (center 0 0 0) (radius 0.2) (height infinity)
                       (material (make dielectric (epsilon 12))))))
</pre>

<p>Here, we've set several properties of the cylinder: the
<code>center</code> is the origin, its <code>radius</code> is 0.2, and
its <code>height</code> (the length along its axis) is
<code>infinity</code>.  Another property, the <code>material</code>,
it itself an object--we made it a dielectric with the property that
its <code>epsilon</code> is 12.  There is another property of the
cylinder that we can set, the direction of its axis, but we're happy
with the default value of pointing in the z direction.

<p>All of the geometric objects are ostensibly three-dimensional, but
since we're doing a two-dimensional simulation the only thing that
matters is their intersection with the xy plane (z=0).  Speaking of
which, let us set the dimensionality of the system.  Normally, we do
this when we define the size of the computational cell, controlled by
the <code>geometry-lattice</code> variable, an object of the
<code>lattice</code> class: we can set some of the dimensions to have
a size <code>no-size</code>, which reduces the dimensionality of the
system.

<pre>
(set! geometry-lattice (make lattice (size 1 1 no-size)))
</pre>

<p>Here, we define a 1x1 two-dimensional cell (defaulting to square).
This cell is <i>discretized</i> according to the
<code>resolution</code> variable, which defaults to <code>10</code>
(pixels/lattice-unit).  That's on the small side, and this is only a
2d calculation, so let's increase the resolution:

<pre>
(set! resolution 32)
</pre>

<p>This results in a 32x32 computational grid.  For efficient
calculation, it is best to make the grid sizes a power of two, or
factorizable into powers of small primes (such as 2, 3, 5 and 7).  As
a rule of thumb, you should use a resolution of at least
<code>8</code> in order to obtain reasonable accuracy.

<p>Now, we're done setting the parameters--there are other parameters,
but we're happy with their default values for now.  At this point,
we're ready to go ahead and compute the band structure.  The simplest
way to do this is to type <code>(run)</code>. Since this is a
two-dimensional calculation, however, we would like to split the bands
up into TE- and TM-polarized modes, and we do this by invoking
<code>(run-te)</code> and <code>(run-tm)</code>.

<p>These produce a lot of output, showing you exactly what the code is
doing as it runs.  Most of this is self-explanatory, but we should
point out one output in particular.  Among the output, you should see
lines like:

<p><code>tefreqs:, 13, 0.3, 0.3, 0, 0.424264, 0.372604, 0.540287, 0.644083, 0.81406, 0.828135, 0.890673, 1.01328, 1.1124</code>

<p>These lines are designed to allow you to easily extract the
band-structure information and import it into a spreadsheet for
graphing.  They comprise the k point index, the k components and
magnitude, and the frequencies of the bands, in comma-delimited
format.  Each line is prefixed by "tefreqs" for TE bands, "tmfreqs"
for TM bands, and "freqs" for ordinary bands (produced by
<code>(run)</code>), and using this prefix you can extract the data
you want from the output by passing it through a program like
<code>grep</code>.  For example, if you had redirected the output to a
file <code>foo.out</code> (as described earlier), you could extract
the TM bands by running <code>grep tmfreqs foo.out</code> at the Unix
prompt.  Note that the output includes a header line, like:

<p><code>tefreqs:, k index, kx, ky, kz, kmag/2pi, band 1, band 2, band 3, band 4, band 5, band 6, band 7, band 8</code>

<p>explaining each column of data.  Another output of the
<code>run</code> is the list of band gaps detected in the computed
bands.  For example the <code>(run-tm)</code> output includes the
following gap output:

<pre>
Gap from band 1 (0.282623311147724) to band 2 (0.419334798706834), 38.9514660888911%
Gap from band 4 (0.715673834754345) to band 5 (0.743682920649084), 3.8385522650349%
</pre>

<p>This data is also stored in the variable <code>gap-list</code>,
which is a list of <code>(<i>gap-percent gap-min gap-max</i>)</code>
lists.  It is important to realize, however, that this band-gap data
may include "false positives," from two possible sources:

<ul>

<li>If two bands cross, a false gap may result because the code
computes the gap (connecting the dots in the band diagram) by assuming
that bands never cross.  Such false gaps are typically quite small
(&lt; 1%).  To be sure of what's going on, you should either look at
the symmetry of the modes involved or compute k points very close to
the crossing (although even if the crossing occurs precisely at one of
your k-points, there usually won't be an exact degeneracy for numerical
reasons).

<li>One typically computes band diagrams by considering k-points
around the boundary of the irreducible Brillouin zone.  It is possible
(though rare) that the band edges may occur at points in the interior
of the Brillouin zone.  To be absolutely sure you have a band gap (and
of its size), you should compute the frequencies for points inside the
Brillouin zone, too.

</ul>

<p>You've computed the band structure, and extracted the
eigenfrequencies for each k point.  But what if you want to see what
the fields look like, or check that the dielectric function is what
you expect?  To do this, you need to output <a
href="http://hdf.ncsa.uiuc.edu/" title="HDF Home Page">HDF files</a>
for these functions (HDF is a binary format for multi-dimensional
scientific data, and can be read by many visualization programs).  (We
output files in HDF5 format, where the filenames are suffixed by
"<code>.h5</code>".)

<p>When you invoke one of the <code>run</code> functions, the
dielectric function in the unit cell is automatically written to the
file "<code>epsilon.h5</code>".  To output the fields or other
information, you need to pass one or more arguments to the
<code>run</code> function.  For example:

<pre>
(run-tm output-efield-z)
(run-te (output-at-kpoint (vector3 0.5 0 0) output-hfield-z output-dpwr))
</pre>

<p>This will output the electric (E) field z components for the TM
bands at all k-points; and the magnetic (H) field z components and
electric field energy densities (D power) for the TE bands at the X
point only.  The output filenames will be things like
"<code>e.k12.b03.z.te.h5</code>", which stands for the z component
(<code>.z</code>) of the TE (<code>.te</code>) electric field
(<code>e</code>) for the third band (<code>.b03</code>) of the twelfth
k point (<code>.k12</code>).  Each HDF5 file can contain multiple
datasets; in this case, it will contain the real and imaginary parts
of the field (in datasets "z.r" and "z.i"), and in general it may
include other data too (e.g. <code>output-efield</code> outputs all
the components in one file).  See also the <a
href="analysis-tutorial.html">Data Analysis</a> tutorial.

<p>There are several other output functions you can pass, described in
the <a href="user-ref.html">user reference section</a>, like
<code>output-dfield</code>, <code>output-hpwr</code>, and
<code>output-dpwr-in-objects</code>.  Actually, though, you can pass
in arbitrary functions that can do much more than just output the
fields--you can perform arbitrary analyses of the bands (using
functions that we will describe later).

<p>Instead of calling one of the <code>run</code> functions, it is
also possible to call lower-level functions of the code directly, to
have a finer control of the computation; such functions are described
in the reference section.

<h2><a name="units">A Few Words on Units</a></h2>

<p>In principle, you can use any units you want with the
Photonic-Bands package.  You can input all of your distances and
coordinates in furlongs, for example, as long as you set the size of
the lattice basis accordingly (the <code>basis-size</code> property of
<code>geometry-lattice</code>).  We do not recommend this
practice, however.

<p>Maxwell's equations possess an important property--they are
<em>scale-invariant</em>.  If you multiply all of your sizes by 10,
the solution scales are simply multiplied by 10 likewise (while the
frequencies are divided by 10).  So, you can solve a problem once and
apply that solution to all length-scales.  For this reason, we usually
pick some fundamental lengthscale <i>a</i> of a structure, such as its
lattice constant (unit of periodicity), and write all distances in
terms of that.  That is, we choose units so that <i>a</i> is unity.
Then, to apply to any physical system, one simply scales all distances
by <i>a</i>.  This is what we have done in the preceding and following
examples, and this is the default behavior of Photonic-Bands: the
lattice constant is one, and all coordinates are scaled accordingly.

<p>As has been mentioned already, (nearly) all 3-vectors in the
program are specified in the <em>basis</em> of the lattice vectors
<em>normalized</em> to lengths given by <code>basis-size</code>,
defaulting to the <em>unit-normalized</em> lattice vectors.  That is,
each component is multiplied by the corresponding basis vector and
summed to compute the corresponding cartesian vector.  It is worth
noting that a basis is not meaningful for scalar distances (such as
the cylinder radius), however: these are just the ordinary cartesian
distances in your chosen units of <i>a</i>.

<p>Note also that the <a href="user-ref.html#k-points">k-points</a>,
as mentioned above, are an exception: they are in the basis of the
reciprocal lattice vectors.  If a given dimension has size
<code>no-size</code>, its reciprocal lattice vector is taken to be
2*pi/<i>a</i>.

<p>We provide <a href="user-ref.html#coordconvert">conversion
functions</a> to transform vectors between the various bases.

<p>The frequency eigenvalues returned by the program are in units of
<i>c/a</i>, where <i>c</i> is the speed of light and <i>a</i> is the
unit of distance.  (Thus, the corresponding vacuum wavelength is
<i>a</i> over the frequency eigenvalue.)

<h2><a name="tri-rods">Bands of a Triangular Lattice</a></h2>

<p>As a second example, we'll compute the TM band structure of a
<em>triangular</em> lattice of dielectric rods in air.  To do this, we
only need to change the lattice, controlled by the variable
<code>geometry-lattice</code>.  We'll set it so that the first two
basis vectors (the properties <code>basis1</code> and
<code>basis2</code>) point 30 degrees above and below the x axis,
instead of their default value of the x and y axes:

<pre>
(set! geometry-lattice (make lattice (size 1 1 no-size)
			 (basis1 (/ (sqrt 3) 2) 0.5)
			 (basis2 (/ (sqrt 3) 2) -0.5)))
</pre>

<p>We don't specify <code>basis3</code>, keeping its default value of
the z axis.  Notice that Scheme supplies us all the ordinary
arithmetic operators and functions, but they use prefix (Polish)
notation, in Scheme fashion.  The <code>basis</code> properties only
specify the directions of the lattice basis vectors, and not their
lengths--the lengths default to unity, which is fine here.

<p>The irreducible Brillouin zone of a triangular lattice is different from
that of a square lattice, so we'll need to modify the
<code>k-points</code> list accordingly:

<pre>
(set! k-points (list (vector3 0 0 0)          ; Gamma
                     (vector3 0 0.5 0)        ; M
                     (vector3 (/ -3) (/ 3) 0) ; K
                     (vector3 0 0 0)))        ; Gamma
(set! k-points (interpolate 4 k-points))
</pre>

<p>Note that these vectors are in the basis of the new reciprocal
lattice vectors, which are different from before.  (Notice also the
Scheme shorthand <code>(/ 3)</code>, which is the same as <code>(/ 1
3)</code> or 1/3.)

<p>All of the other parameters (<code>geometry</code>,
<code>num-bands</code>, and <code>grid-size</code>) can remain the
same as in the previous subsection, so we can now call
<code>(run-tm)</code> to compute the bands.  As it turns out, this
structure has an even larger TM gap than the square lattice:

<pre>
Gap from band 1 (0.275065617068082) to band 2 (0.446289918847647), 47.4729292989213%
</pre>

<h2><a name="max-gap">Maximizing the First TM Gap</a></h2>

<p>Just for fun, we will now show you a more sophisticated example,
utilizing the programming capabilities of Scheme: we will write a
script to choose the cylinder radius that maximizes the first TM gap
of the triangular lattice of rods from above.  All of the Scheme
syntax here won't be explained, but this should give you a flavor of
what is possible.

<p>First, we will write the function that want to maximize, a function
that takes a dielectric constant and returns the size of the first TM
gap.  This function will change the geometry to reflect the new
radius, run the calculation, and return the size of the first gap:

<pre>
(define (first-tm-gap r)
  (set! geometry (list (make cylinder
			 (center 0 0 0) (radius r) (height infinity)
			 (material (make dielectric (epsilon 12))))))
  (run-tm)
  (retrieve-gap 1)) ; return the gap from TM band 1 to TM band 2
</pre>

<p>We'll leave most of the other parameters the same as in the
previous example, but we'll also change <code>num-bands</code> to 2,
since we only need to compute the first two bands:

<pre>
(set! num-bands 2)
</pre>

<p>In order to distinguish small differences in radius during the
optimization, it might seem that we have to increase the grid
resolution, slowing down the computation.  Instead, we can simply
increase the <i>mesh</i> resolution.  This is the size of the mesh
over which the dielectric constant is averaged at each grid point, and
increasing the mesh size means that the average index better reflects
small changes in the structure.

<pre>
(set! mesh-size 7) ; increase from default value of 3
</pre>

<p>Now, we're ready to maximize our function
<code>first-tm-gap</code>.  We could write a loop to do this
ourselves, but libctl provides a built-in function <code>(maximize
function tolerance arg-min arg-max)</code> to do it for us (using
Brent's algorithm).  So, we just tell it to find the maximum,
searching in the range of radii from 0.1 to 0.5, with a tolerance of 0.1:

<pre>
(define result (maximize first-tm-gap 0.1 0.1 0.5))
(print "radius at maximum: " (max-arg result) "\n")
(print "gap size at maximum: " (max-val result) "\n")
</pre>

<p>(<code>print</code> is a function defined by libctl to apply
the built-in <code>display</code> function to zero or more arguments.)
After five iterations, the output is:

<pre>
radius at maximum: 0.176393202250021
gap size at maximum: 48.6252611051049
</pre>

<p>The tolerance of 0.1 that we specified means that the true maximum
is within 0.1 * 0.176393202250021, or about 0.02, of the radius found
here.  It doesn't make much sense here to specify a lower tolerance,
since the discretization of the grid means that the code can't
accurately distinguish small differences in radius.

<p>Before we continue, let's reset <code>mesh-size</code> to its
default value:

<pre>
(set! mesh-size 3) ; reset to default value of 3
</pre>

<h2><a name="aniso">A Complete 2D Gap with an Anisotropic Dielectric</a></h2>

<p>As another example, one which does not require so much Scheme
knowledge, let's construct a structure with a complete 2D gap (i.e.,
in both TE and TM polarizations), in a somewhat unusual way: using an
<a href="user-ref.html#dielectric-anisotropic">anisotropic dielectric</a>
structure.  An anisotropic dielectric presents a different
dielectric constant depending upon the direction of the electric
field, and can be used in this case to make the TE and TM
polarizations "see" different structures.

<p>We already know that the triangular lattice of rods has a gap for
TM light, but not for TE light.  The dual structure, a triangular
lattice of holes, has a gap for TE light but not for TM light (at
least for the small radii we will consider).  Using an anisotropic
dielectric, we can make both of these structures simultaneously, with
each polarization seeing the structure that gives it a gap.

<p>As before, our <code>geometry</code> will consist of a single
cylinder, this time with a radius of 0.3, but now it will have an
epsilon of 12 (dielectric rod) for TM light and 1 (air hole) for TE
light:

<pre>
(set! geometry (list (make cylinder
                       (center 0 0 0) (radius 0.3) (height infinity)
                       (material (make dielectric-anisotropic
                                   (epsilon-diag 1 1 12))))))
</pre>

<p>Here, <code>epsilon-diag</code> specifies the diagonal elements of
the dielectric tensor.  The off-diagonal elements (specified by
<code>epsilon-offdiag</code>) default to zero and are only needed when
the principal axes of the dielectric tensor are different from the
cartesian xyz axes.

<p>The background defaults to air, but we need to make it a dielectric
(epsilon of 12) for the TE light, so that the cylinder forms a hole.
This is controlled via the <code>default-material</code> variable:

<pre>
(set! default-material (make dielectric-anisotropic (epsilon-diag 12 12 1)))
</pre>

<p>Finally, we'll increase the number of bands back to eight and run
the computation:

<pre>
(set! num-bands 8)
(run) ; just use run, instead of run-te or run-tm, to find the complete gap
</pre>

<p>The result, as expected, is a complete band gap:

<pre>
Gap from band 2 (0.223977612336924) to band 3 (0.274704473679751), 20.3443687933601%
</pre>

<p>(If we had computed the TM and TE bands separately, we would have
found that the lower edge of the complete gap in this case comes from
the TM gap, and the upper edge comes from the TE gap.)

<h2><a name="point-defect">Finding a Point-defect State</a></h2>

<p>Here, we consider the problem of finding a point-defect state in
our square lattice of rods.  This is a state that is localized in a
small region by creating a point defect in the crystal--e.g., by
removing a single rod.  The resulting mode will have a frequency
within, and be confined by, the gap.

<p>To compute this, we need a supercell of bulk crystal, within which
to put the defect--we will use a 5x5 cell of rods.  To do this, we
must first increase the size of the lattice by five, and then add all
of the rods.  We create the lattice by:

<pre>
(set! geometry-lattice (make lattice (size 5 5 no-size)))
</pre>

<p>Here, we have used the default orthogonal basis, but have changed
the size of the cell.  To populate the cell, we could specify all 25
rods manually, but that would be tedious--and any time you find
yourself doing something tedious in a programming language, you're
making a mistake.  A better approach would be to write a loop, but in
fact this has already been done for you.  Photonic-Bands provides a
function, <code>geometric-objects-lattice-duplicates</code>, that
duplicates a list of objects over the lattice:

<pre>
(set! geometry (list (make cylinder
                       (center 0 0 0) (radius 0.2) (height infinity)
                       (material (make dielectric (epsilon 12))))))
(set! geometry (geometric-objects-lattice-duplicates geometry))
</pre>

<p>There, now the <code>geometry</code> list contains 25 rods--the
originaly <code>geometry</code> list (which contained one rod,
remember) duplicated over the 5x5 lattice.

<p>To remove a rod, we'll just add another rod in the center of the
cell with a dielectric constant of 1.  When objects overlap, the later
object in the list takes precedence, so we have to put the new rod at
the end of <code>geometry</code>:

<pre>
(set! geometry (append geometry 
		       (list (make cylinder (center 0 0 0) 
				   (radius 0.2) (height infinity)
				   (material air)))))
</pre>

<p>Here, we've used the Scheme <code>append</code> function to combine
two lists, and have also snuck in the predefined material type
<code>air</code> (which has an epsilon of 1).

<p>We'll be frugal and only use only 16 points per lattice unit
(resulting in an 80x80 grid) instead of the 32 from before:

<pre>
(set! resolution 16)
</pre>

<p>Only a single k point is needed for a point-defect calculation
(which, for an infinite supercell, would be independent of k):

<pre>
(set! k-points (list (vector3 0.5 0.5 0)))
</pre>

<p>Unfortunately, for a supercell the original bands are folded many
times over (in this case, 25 times), so we need to compute many more
bands to reach the same frequencies:

<pre>
(set! num-bands 50)
</pre>

<p>At this point, we can call <code>(run-tm)</code> to solve for the
TM bands.  It will take several minutes to compute, so be patient.
Recall that the gap for this structure was for the frequency range
0.2812 to 0.4174.  The bands of the solution include exactly one state
in this frequency range: band 25, with a frequency of 0.378166.  This
is exactly what we should expect--the lowest band was folded 25 times
into the supercell Brillouin zone, but one of these states was pushed
up into the gap by the defect.

<p>We haven't yet output any of the fields, but we don't have to
repeat the run to do so.  The fields from the last k-point computation
remain in memory and can continue to be accessed and analyzed.  For
example, to output the electric field z component of band 25, we just
do:

<pre>
(output-efield-z 25)
</pre>

<p>That's right, the output functions that we passed to
<code>(run)</code> in the first example are just functions of the band
index that are called on each band.  We can do other computations too,
like compute the fraction of the electric field energy near the defect
cylinder (within a radius 1.0 of the origin):

<pre>
(get-dfield 25)  ; compute the D field for band 25
(compute-field-energy)  ; compute the energy density from D
(print
 "energy in cylinder: "
 (compute-energy-in-objects (make cylinder (center 0 0 0)
                                  (radius 1.0) (height infinity)
                                  (material air)))
 "\n")
</pre>

<p>The result is <code>0.624794702341156</code>, or over 62% of the
field energy in this localized region; the field decays exponentially
into the bulk crystal.  The full range of available functions is
described in the <a href="user-ref.html">reference section</a>, but
the typical sequence is to first load a field with a <code>get-</code>
function and then to call other functions to perform computations and
transformations on it.

<p>Note that the <code>compute-energy-in-objects</code> returns the
energy fraction, but does not itself print this value.  This is fine
when you are running interactively, in which case Guile always
displays the result of the last expression, but when running as part
of a script you need to explicitly print the result as we have done
above with the <code>print</code> function.  The <code>"\n"</code>
string is newline character (like in C), to put subsequent output on a
separate line.

<p>Instead of computing all those bands, we can instead take advantage
of a special feature of the MIT Photonic-Bands software that allows you to
compute the bands closest to a "target" frequency, rather than the
bands with the lowest frequencies.  One uses this feature by setting
the <code>target-freq</code> variable to something other than zero
(e.g. the mid-gap frequency).  In order to get accurate results, it's
currently also recommended that you decrease the
<code>tolerance</code> variable, which controls when convergence is
judged to have occurred, from its default value of <code>1e-7</code>:

<pre>
(set! num-bands 1)  ; only need to compute a single band, now!
(set! target-freq (/ (+ 0.2812 0.4174) 2))
(set! tolerance 1e-8)
</pre>

<p>Now, we just call <code>(run-tm)</code> as before.  Convergence
requires more iterations this time, both because we've decreased the
tolerance and because of the nature of the eigenproblem that is now
being solved, but only by about 3-4 times in this case.  Since we now
have to compute only a single band, however, we arrive at an answer
much more quickly than before.  The result, of course, is again the
defect band, with a frequency of 0.378166.

<h2><a name="tune-defect">Tuning the Point-defect Mode</a></h2>

<p>As another example utilizing the programming power of Scheme, we
will write a script to "tune" the defect mode to a particular
frequency.  Instead of forming a defect by simply removing a rod, we
can decrease the radius or the dielectric constant of the defect rod,
thereby changing the corresponding mode frequency.  In this case,
we'll vary the dielectric constant, and try to find a mode with a
frequency of, say, 0.314159 (a random number).

<p>We could write a loop to search for this epsilon, but instead we'll
use a root-finding function provided by libctl, <code>(find-root
<i>function tolerance arg-min arg-max</i>)</code>, that will solve the
problem for us using a quadratically-convergent algorithm (Ridder's
method).  First, we need to define a function that takes an epsilon
for the center rod and returns the mode frequency minus 0.314159; this
is the function we'll be finding the root of:

<pre>
(define old-geometry geometry) ; save the 5x5 grid with a missing rod
(define (rootfun eps)
  ; add the cylinder of epsilon = eps to the old geometry:
  (set! geometry (append old-geometry
			 (list (make cylinder (center 0 0 0)
				     (radius 0.2) (height infinity)
				     (material (make dielectric
						 (epsilon eps)))))))
  (run-tm)  ; solve for the mode (using the targeted solver)
  (print "epsilon = " eps " gives freq. = " (list-ref freqs 0) "\n")
  (- (list-ref freqs 0) 0.314159))  ; return 1st band freq. - 0.314159
</pre>

<p>Now, we can solve for epsilon (searching in the range 1 to 12, with
a fractional tolerance of 0.01) by:

<pre>
(define rooteps (find-root rootfun 0.01 1 12))
(print "root (value of epsilon) is at: " rooteps "\n")
</pre>

<p>The sequence of dielectric constants that it tries, along with the
corresponding mode frequencies, is:

<table cellspacing=2 cellpadding=1>
<tr align=center> <th>epsilon</th> <th>frequency</th> </tr>
<tr align=center> <td>1</td> <td>0.378165893321125</td> </tr>
<tr align=center> <td>12</td> <td>0.283987088221692</td> </tr>
<tr align=center> <td>6.5</td> <td>0.302998920718043</td> </tr>
<tr align=center> <td>5.14623274327171</td> <td>0.317371748739314</td> </tr>
<tr align=center> <td>5.82311637163586</td> <td>0.309702408341706</td> </tr>
<tr align=center> <td>5.41898003340128</td> <td>0.314169110036439</td> </tr>
<tr align=center> <td>5.62104820251857</td> <td>0.311893530112625</td> </tr>
</table>

<p>The final answer that it returns is an epsilon of 5.41986120170136
Interestingly enough, the algorithm doesn't actually evaluate the
function at the final point; you have to do so yourself if you want to
find out how close it is to the root.  (Ridder's method successively
reduces the interval bracketing the root by alternating bisection and
interpolation steps.  At the end, it does one last interpolation to
give you its best guess for the root location within the current
interval.)  If we go ahead and evaluate the band frequency at this
dielectric constant (calling <code>(rootfun rooteps)</code>), we find
that it is 0.314159008193209, matching our desired frequency to nearly
eight decimal places after seven function evaluations!  (Of course,
the computation isn't really this accurate anyway, due to the finite
discretization.)

<p>A slight improvement can be made to the calculation above.
Ordinarily, each time you call the <code>(run-tm)</code> function, the
fields are initialized to random values.  It would speed convergence
somewhat to use the fields of the previous calculation as the starting
point for the next calculation.  We can do this by instead calling a
lower-level function, <code>(run-parity TM false)</code>.  (The
first parameter is the polarization to solve for, and the second tells
it not to reset the fields if possible.)

<h2><a name="emacs">emacs and ctl</a></h2>

<p>It is useful to have <code>emacs</code> use its
<code>scheme-mode</code> for editing ctl files, so that hitting tab
indents nicely, and so on.  <code>emacs</code> does this automatically
for files ending with "<code>.scm</code>"; to do it for files ending
with "<code>.ctl</code>" as well, add the following lines to your
<code>~/.emacs</code> file:

<pre>
(if (assoc "\\.ctl" auto-mode-alist)
    nil
  (nconc auto-mode-alist '(("\\.ctl" . scheme-mode))))
</pre>

<p>(Incidentally, <code>emacs</code> scripts are written in "elisp," a
language closely related to Scheme.)

<hr>
Go to the <a href="analysis-tutorial.html">next</a>, <a href="installation.html">previous</a>, or <a href="index.html">main</a> section.

</BODY>
</HTML>