user-tutorial.html

麻省理工的计算光子晶体的程序

<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