user-tutorial.html

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

(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>