# Scheme Tutorial

We'll go through several examples using the Scheme interface that demonstrate the process of computing fields, transmittance/reflectance spectra, and resonant modes. The examples are 1d or 2d calculations, simply because they are quicker than 3d and they illustrate most of the essential features. For more advanced functionality involving 3d computations with a focus on technology applications, see the Simpetus projects page.

In order to convert the HDF5 output files of Meep into images of the fields, this tutorial uses the h5utils package. You could also use any other package (i.e., Octave or Matlab) that supports reading HDF5 files.

## The Scheme Script File

The use of Meep revolves around the script (or control) file, abbreviated "ctl" and typically called something like foo.ctl. The script file specifies the geometry, the current sources, the outputs computed, and everything else specific to your calculation. Rather than a flat, inflexible file format, however, the script 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.

Don't worry, though — simple things are simple and you don't need to be an experienced programmer. 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, etc.

The script 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 script file:

• Scheme: a programming language developed at MIT. The syntax is particularly simple: all statements are of the form (function arguments...). We run Scheme under the Guile interpreter which is designed to be plugged into programs as a scripting and extension language. You don't need to know much Scheme for a basic script file, but it is always there if you need it. For more details, see Guile and Scheme Information.
• libctl: a library built on top of Guile to simplify protocols for scientific computation. libctl sets the basic tone of the interface and defines a number of useful functions (such as multi-variable optimization, numeric integration, and so on). See the libctl documentation.
• Meep itself, which defines all the interface features that are specific to FDTD calculations. This manual is primarily focused on documenting these features.

At this point, please take a moment to leaf through the libctl tutorial to get a feel for the basic style of the interface, before we get to the Meep-specific stuff below. MPB has a similar interface.

Let's continue with our tutorial. The Meep program is normally invoked by running something like the following at the Unix command line:

 unix% meep foo.ctl >& foo.out


which reads foo.ctl and executes it, saving the output to the file foo.out. However, if you invoke meep with no arguments, you are dropped into an interactive mode in which you can type commands and see their results immediately. You can paste in the commands from the tutorial as you follow it and see what they do.

## Fields in a Waveguide

For our first example, let's examine the field pattern excited by a localized CW source in a waveguide — first straight, then bent. The waveguide will have frequency-independent ε=12 and width 1 μm. The unit length in this example is 1 μm. See also Units.

### A Straight Waveguide

The simulation script is in examples/straight-waveguide.ctl.

Before we define the structure, however, we have to define the computational cell. We're going to put a source at one end and watch it propagate down the waveguide in the x direction, so let's use a cell of length 16 μm in the x direction to give it some distance to propagate. In the y direction, we just need enough room so that the boundaries don't affect the waveguide mode; let's give it a size of 8 μm. We specify these sizes via the geometry-lattice variable:

(set! geometry-lattice (make lattice (size 16 8 no-size)))


The name geometry-lattice comes from MPB, where it can be used to define a more general periodic lattice. Although Meep supports periodic structures, it is less general than MPB in that affine grids are not supported. set! is a Scheme command to set the value of an input variable. The last no-size parameter indicates that the cell has no size in the z direction, i.e. it is two-dimensional.

We can add the waveguide. Most commonly, the structure is specified by a list of geometric-objects, stored in the geometry variable.

(set! geometry (list
(make block (center 0 0) (size infinity 1 infinity)
(material (make medium (epsilon 12))))))


The waveguide is specified by a block (parallelepiped) of size $\infty \times 1 \times \infty$, with ε=12, centered at (0,0) which is the center of the cell. By default, any place where there are no objects there is air (ε=1), although this can be changed by setting the default-material variable. The resulting structure is shown below.

We have the structure and need to specify the current sources using the sources object. The simplest thing is to add a single point source $J_z$:

(set! sources (list
(make source
(src (make continuous-src (frequency 0.15)))
(component Ez)
(center -7 0))))


We gave the source a frequency of 0.15, and specified a continuous-src which is just a fixed-frequency sinusoid $\exp(-i \omega t)$ that by default is turned on at $t$=0. Recall that, in Meep units, frequency is specified in units of 2πc, which is equivalent to the inverse of the vacuum wavelength. Thus, 0.15 corresponds to a vacuum wavelength of about 1/0.15=6.67 μm, or a wavelength of about 2 μm in the ε=12 material — thus, our waveguide is half a wavelength wide, which should hopefully make it single mode. In fact, the cutoff for single-mode behavior in this waveguide is analytically solvable, and corresponds to a frequency of 1/2√11 or roughly 0.15076. Note also that to specify a Jz, we specify a component $Ez$ (e.g. if we wanted a magnetic current, we would specify Hx, Hy, or Hz). The current is located at (-7,0), which is 1 μm to the right of the left edge of the cell — we always want to leave a little space between sources and the cell boundaries, to keep the boundary conditions from interfering with them.

As for boundary conditions, we want to add absorbing boundaries around our cell. Absorbing boundaries in Meep are handled by perfectly matched layers (PML) — which aren't really a boundary condition at all, but rather a fictitious absorbing material added around the edges of the cell. To add an absorbing layer of thickness 1 μm around all sides of the cell, we do:

(set! pml-layers (list (make pml (thickness 1.0))))


pml-layers is a list of pml objects — you may have more than one pml object if you want PML layers only on certain sides of the cell, e.g. (make pml (thickness 1.0) (direction X) (side High)) specifies a PML layer on only the $+x$ side. We note an important point: the PML layer is inside the cell, overlapping whatever objects you have there. So, in this case our PML overlaps our waveguide, which is what we want so that it will properly absorb waveguide modes. The finite thickness of the PML is important to reduce numerical reflections. For more information, see Perfectly Matched Layer.

Meep will discretize this structure in space and time, and that is specified by a single variable, resolution, that gives the number of pixels per distance unit. We'll set this resolution to 10 pixels/μm, which corresponds to around 67 pixels/wavelength, or around 20 pixels/wavelength in the high-index material. In general, at least 8 pixels/wavelength in the highest dielectric is a good idea. This will give us a 160×80 cell.

(set! resolution 10)


We are ready to run the simulation. We do this by calling the run-until function. The first argument to run-until is the time to run for, and the subsequent arguments specify fields to output or other kinds of analyses at each time step:

(run-until 200
(at-beginning output-epsilon)
(at-end output-efield-z))


We are outputting the dielectric function $\epsilon$ and the electric-field component $E_z$, but have wrapped the output functions which would otherwise run at every time step in at-beginning and at-end, which do just what they say. There are several other such functions to modify the output behavior — and you can, of course, write your own, and in fact you can do any computation or output you want at any time during the time evolution and even modify the simulation while it is running.

It should complete in a few seconds. If you are running interactively, the two output files will be called eps-000000.00.h5 and ez-000200.00.h5 (notice that the file names include the time at which they were output). If we were running a tutorial.ctl file, then the outputs will be tutorial-eps-000000.00.h5 and tutorial-ez-000200.00.h5. In any case, we can analyze and visualize these files with a wide variety of packages that support the HDF5 format, including our own h5utils, and in particular the h5topng program to convert them to PNG images.

unix% h5topng -S3 eps-000000.00.h5


This will create eps-000000.00.png, where the -S3 increases the image scale by 3 so that it is around 450 pixels wide, in this case. In fact, precisely this command is what created the dielectric image above. Much more interesting, however, are the fields:

unix% h5topng -S3 -Zc dkbluered -a yarg -A eps-000000.00.h5 ez-000200.00.h5


Briefly, the -Zc dkbluered makes the color scale go from dark blue (negative) to white (zero) to dark red (positive), and the -a/-A options overlay the dielectric function as light gray contours. This results in the image:

We see that the the source has excited the waveguide mode, but has also excited radiating fields propagating away from the waveguide. At the boundaries, the field quickly goes to zero due to the PML layers. If we look carefully, we see something else — the image is "speckled" towards the right side. This is because, by turning on the current abruptly at $t=0$, we have excited high-frequency components (very high order modes), and we have not waited long enough for them to die away; we'll eliminate these in the next section by turning on the source more smoothly.

### A 90° Bend

We'll start a new simulation where we look at the fields in a bent waveguide, and we'll do a couple of other things differently as well. The simulation script is in examples/bent-waveguide.ctl. If you are running Meep interactively, you will want to get rid of the old structure and fields so that Meep will re-initialize them:

(reset-meep)


Then let's set up the bent waveguide, in a slightly bigger cell, via:

(set! geometry-lattice (make lattice (size 16 16 no-size)))

(set! geometry (list
(make block (center -2 -3.5) (size 12 1 infinity)
(material (make medium (epsilon 12))))
(make block (center 3.5 2) (size 1 12 infinity)
(material (make medium (epsilon 12))))))

(set! pml-layers (list (make pml (thickness 1.0))))

(set! resolution 10)


Note that we have two blocks, both off-center to produce the bent waveguide structure pictured at right. As illustrated in the figure, the origin (0,0) of the coordinate system is at the center of the cell, with positive $y$ being downwards in h5topng, and thus the block of size 12$\times$1 is centered at (-2,-3.5). Also shown in green is the source plane at $x=-7$ which is shifted to $y=-3.5$ so that it is still inside the waveguide.

There are a couple of items to note. First, a point source does not couple very efficiently to the waveguide mode, so we'll expand this into a line source the same width as the waveguide by adding a size property to the source. An eigenmode source can also be used which is described in Tutorial/Optical Forces. Second, instead of turning the source on suddenly at t=0 which excites many other frequencies because of the discontinuity, we will ramp it on slowly. Meep uses a hyperbolic tangent (tanh) turn-on function over a time proportional to the width of 20 time units which is a little over three periods. Finally, just for variety, we'll specify the vacuum wavelength instead of the frequency; again, we'll use a wavelength such that the waveguide is half a wavelength wide.

(set! sources (list
(make source
(src (make continuous-src (wavelength (* 2 (sqrt 11))) (width 20)))
(component Ez)
(center -7 -3.5)
(size 0 1))))


Finally, we'll run the simulation. Instead of running output-efield-z only at the end of the simulation, however, we'll run it at every 0.6 time units (about 10 times per period) via (at-every 0.6 output-efield-z). By itself, this would output a separate file for every different output time, but instead we'll use another feature of Meep to output to a single three-dimensional HDF5 file, where the third dimension is time:

(run-until 200
(at-beginning output-epsilon)
(to-appended "ez" (at-every 0.6 output-efield-z)))


"ez" determines the name of the output file, which will be called ez.h5 if you are running interactively or will be prefixed with the name of the file name for a script file (e.g. tutorial-ez.h5 for tutorial.ctl). If we run h5ls on this file (a standard utility, included with HDF5, that lists the contents of the HDF5 file), we get:

unix% h5ls ez.h5
ez                       Dataset {161, 161, 330/Inf}


That is, the file contains a single dataset ez that is a 160$\times$160$\times$333 array, where the last dimension is time. This is rather a large file, 66MB; later, we'll see ways to reduce this size if we only want images. We have a number of choices of how to output the fields. To output a single time slice, we can use the same h5topng command as before, but with an additional -t option to specify the time index: e.g. h5topng -t 332 will output the last time slice, similar to before. Instead, let's create an animation of the fields as a function of time. First, we have to create images for all of the time slices:

unix% h5topng -t 0:332 -R -Zc dkbluered -a yarg -A eps-000000.00.h5 ez.h5


This is similar to the command before, with two new options: -t 0:332 outputs images for all time indices from 0 to 332, i.e. all of the times, and the the -R flag tells h5topng to use a consistent color scale for every image instead of scaling each image independently. Then, we have to convert these images into an animation in some format. For this, we'll use the free ImageMagick convert program and there are other tools that work as well.

unix% convert ez.t*.png ez.gif


We are using an animated GIF format for the output. This results in the following animation:

It is clear that the transmission around the bend is rather low for this frequency and structure — both large reflection and large radiation loss are clearly visible. Moreover, since we are operating just barely below the cutoff for single-mode behavior, we are able to excite a second leaky mode after the waveguide bend, whose second-order mode pattern (superimposed with the fundamental mode) is apparent in the animation. At right, we show a field snapshot from a simulation with a larger cell along the $y$ direction, in which you can see that the second-order leaky mode decays away, leaving us with the fundamental mode propagating downward.

Instead of doing an animation, another interesting possibility is to make an image from a $x \times t$ slice. Here is the $y=-3.5$ slice, which gives us an image of the fields in the first waveguide branch as a function of time.

unix% h5topng -0y -35 -Zc dkbluered ez.h5


The -0y -35 specifies the $y=-3.5$ slice, where we have multiplied by 10 (our resolution) to get the pixel coordinate.

#### Output Tips and Tricks

Above, we outputted the full 2d data slice at every 0.6 time units, resulting in a 69MB file. This is not large but you can imagine how big the output file would get if we were doing a 3d simulation, or even a larger 2d simulation — one can easily generate gigabytes of files, which is not only wasteful but is also slow. Instead, it is possible to output more efficiently if you know what you want to look at.

To create the movie above, all we really need are the images corresponding to each time. Images can be stored much more efficiently than raw arrays of numbers — to exploit this fact, Meep allows you to output PNG images instead of HDF5 files. In particular, instead of output-efield-z as above, we can use (output-png Ez "-Zc dkbluered"), where Ez is the component to output and the "-Zc dkbluered" are options for h5topng which is the program that is actually used to create the image files. That is:

(run-until 200 (at-every 0.6 (output-png Ez "-Zc bluered")))


will output a PNG file file every 0.6 time units, which can then be combined with convert as above to create a movie. The movie will be similar to the one before, but not identical because of how the color scale is determined. Before, we used the -R option to make h5topng use a uniform color scale for all images, based on the minimum/maximum field values over all time steps. That is not possible because we output an image before knowing the field values at future time steps. Thus, what output-png does is to set its color scale based on the minimum/maximum field values from all past times — therefore, the color scale will slowly "ramp up" as the source turns on.

The above command outputs zillions of PNG files, and it is somewhat annoying to have them clutter up our directory. Instead, we can use the following command before run-until:

(use-output-directory)


This will put all of the output files (.h5, .png, etcetera) into a newly-created subdirectory, called by default filename-out/ if our script file is filename.ctl.

What if we want to output an $x \times t$ slice, as above? To do this, we only really wanted the values at $y=-3.5$, and therefore we can exploit another powerful Meep output feature — Meep allows us to output only a subset of the cell. This is done using the in-volume function, which similar to at-every and to-appended is another function that modifies the behavior of other output functions. In particular, we can do:

 (run-until 200
(to-appended "ez-slice"
(at-every 0.6
(in-volume (volume (center 0 -3.5) (size 16 0))
output-efield-z))))


The first argument to in-volume is a volume, specified by (volume (center ...) (size ...)), which applies to all of the nested output functions. Note that to-appended, at-every, and in-volume are cumulative regardless of what order you put them in. This creates the output file ez-slice.h5 which contains a dataset of size 160x333 corresponding to the desired $x \times t$ slice.

## Transmittance Spectrum of a Waveguide Bend

We have computed the field patterns for light propagating around a waveguide bend. While this can be visually informative, the results are not quantitatively satisfying. We'd like to know exactly how much power makes it around the bend (transmittance), how much is reflected (reflectance), and how much is radiated away (scattered loss). How can we do this?

The basic principles are described in Introduction. The computation involves keeping track of the fields and their Fourier transform in a certain region, and from this computing the flux of electromagnetic energy as a function of ω. Moreover, we'll get an entire spectrum of the transmittance in a single run, by Fourier-transforming the response to a short pulse. However, in order to normalize the transmitted flux by the incident power to obtain the transmittance, we'll have to do two runs, one with and one without the bend (i.e., a straight waveguide).

The simulation script is in examples/bend-flux.ctl.

Above, we hard-coded all of the parameters like the cell size, the waveguide width, etcetera. For serious work, however, this is inefficient — we often want to explore many different values of such parameters. For example, we may want to change the size of the cell, so we'll define it as:

(define-param sx 16) ; size of cell in X direction
(define-param sy 32) ; size of cell in Y direction
(set! geometry-lattice (make lattice (size sx sy no-size)))


Notice that a semicolon ";" begins a comment, which is ignored by Meep. define-param is a libctl feature to define variables that can be overridden from the command line. We could do meep sx=17 bend-flux.ctl to change the $x$ size to 17, without editing the script file, for example. We'll also define a couple of parameters to set the width of the waveguide and the "padding" between it and the edge of the cell:

(define-param pad 4) ; padding distance between waveguide and cell edge
(define-param w 1)   ; width of waveguide


In order to define the waveguide positions, we will have to use arithmetic. For example, the $y$ center of the horizontal waveguide will be given by -0.5*(sy-w-2*pad). At least, that is what the expression would look like in C; in Scheme, the syntax for 1+2 is (+ 1 2), and so on, so we will define the horizontal and vertical waveguide centers as:

(define wvg-xcen (*  0.5 (- sx w (* 2 pad)))) ; x center of vert. wvg
(define wvg-ycen (* -0.5 (- sy w (* 2 pad)))) ; y center of horiz. wvg


We proceed to define the geometry, as before. This time, however, we really want two geometries: the bend, and also a straight waveguide for normalization. We could do this with two separate script files, but that is annoying. Instead, we'll define a parameter no-bend? which is true for the straight-waveguide case and false for the bend.

(define-param no-bend? false) ; if true, have straight waveguide, not bend


We define the geometry via two cases, with an if statement — the Scheme syntax is (if predicate? if-true if-false).

(set! geometry
(if no-bend?
(list
(make block
(center 0 wvg-ycen)
(size infinity w infinity)
(material (make medium (epsilon 12)))))
(list
(make block
(size (- sx pad) w infinity)
(material (make medium (epsilon 12))))
(make block
(size w (- sy pad) infinity)
(material (make medium (epsilon 12)))))))


Thus, if no-bend? is true we make a single block for a straight waveguide, and otherwise we make two blocks for a bent waveguide.

The source is a gaussian-src instead of a continuous-src, parameterized by a center frequency and a frequency width (the width of the Gaussian spectrum), which we'll define via define-param as usual.

(define-param fcen 0.15) ; pulse center frequency
(define-param df 0.1)    ; pulse width (in frequency)
(set! sources (list
(make source
(src (make gaussian-src (frequency fcen) (fwidth df)))
(component Ez)
(center (+ 1 (* -0.5 sx)) wvg-ycen)
(size 0 w))))


Notice how we're using our parameters like wvg-ycen and w: if we change the dimensions, everything will shift automatically. The boundary conditions and resolution are set as before, except that we'll use set-param! so that we can override the resolution from the command line.:

(set! pml-layers (list (make pml (thickness 1.0))))
(set-param! resolution 10)


Finally, we have to specify where we want Meep to compute the flux spectra, and at what frequencies. This must be done after specifying the geometry, sources, resolution, etcetera, because all of the field parameters are initialized when flux planes are created. As described in Introduction, the flux is the integral of the Poynting vector over the specified flux-region. It only integrates one component of the Poynting vector and the direction property specifies which component. In this example, since the flux-region is a line, the direction is its normal by default which therefore does not need to be explicitly defined.

(define-param nfreq 100) ; number of frequencies at which to compute flux
(define trans ; transmitted flux
(if no-bend?
(make flux-region
(center (- (/ sx 2) 1.5) wvg-ycen) (size 0 (* w 2)))
(make flux-region
(center wvg-xcen (- (/ sy 2) 1.5)) (size (* w 2) 0)))))
(define refl ; reflected flux
(make flux-region
(center (+ (* -0.5 sx) 1.5) wvg-ycen) (size 0 (* w 2)))))


We compute the fluxes through a line segment twice the width of the waveguide, located at the beginning or end of the waveguide. Note that the flux lines are separated by 1 μm from the boundary of the cell, so that they do not lie within the absorbing PML regions. Again, there are two cases: the transmitted flux is either computed at the right or the bottom of the cell, depending on whether the waveguide is straight or bent.

The fluxes will be computed for nfreq=100 frequencies centered on fcen, from fcen-df/2 to fcen+df/2. That is, we only compute fluxes for frequencies within our pulse bandwidth. This is important because, far outside the pulse bandwidth, the spectral power is so low that numerical errors make the computed fluxes useless.

As described in Introduction, computing the reflection spectra requires some care because we need to separate the incident and reflected fields. We do this in Meep by saving the Fourier-transformed fields from the normalization run (no-bend?=true), and loading them, negated, before the other runs. The latter subtracts the Fourier-transformed incident fields from the Fourier transforms of the scattered fields. Logically, we might subtract these after the run, but it turns out to be more convenient to subtract the incident fields first and then accumulate the Fourier transform. All of this is accomplished with two commands, save-flux (after the normalization run) and load-minus-flux (before the other runs). We can call them as follows:

(if (not no-bend?) (load-minus-flux "refl-flux" refl))
(run-sources+ 500 (at-beginning output-epsilon))
(if no-bend? (save-flux "refl-flux" refl))


This uses a file called refl-flux.h5, or actually bend-flux-refl-flux.h5 (the script file name is used as a prefix) to store/load the Fourier transformed fields in the flux planes. The (run-sources+ 500) runs the simulation until the Gaussian source has turned off (which is done automatically once it has decayed for a few standard deviations), plus an additional 500 time units.

Why do we keep running after the source has turned off? Because we must give the pulse time to propagate completely across the cell. Moreover, the time required is a bit tricky to predict when you have complex structures, because there might be resonant phenomena that allow the source to bounce around for a long time. Therefore, it is convenient to specify the run time in a different way: instead of using a fixed time, we require that the |Ez|2 at the end of the waveguide must have decayed by a given amount (e.g. 1/1000) from its peak value. We can do this via:

(run-sources+
(stop-when-fields-decayed 50 Ez
(if no-bend? (vector3 (- (/ sx 2) 1.5) wvg-ycen)
(vector3 wvg-xcen (- (/ sy 2) 1.5)))
1e-3))


stop-when-fields-decayed takes four arguments: (stop-when-fields-decayed dT component pt decay-by). What it does is, after the sources have turned off, it keeps running for an additional dT time units every time the given |component|2 at the given point has not decayed by at least decay-by from its peak value for all times within the previous dT. In this case, dT=50, the component is Ez, the point is at the center of the flux plane at the end of the waveguide, and decay-by=0.001. So, it keeps running for an additional 50 time units until the square amplitude has decayed by 1/1000 from its peak. This should be sufficient to ensure that the Fourier transforms have converged.

Finally, we have to output the flux values:

(display-fluxes trans refl)


This prints a series of outputs like:

flux1:, 0.1, 7.91772317108475e-7, -3.16449591437196e-7
flux1:, 0.101010101010101, 1.18410865137737e-6, -4.85527604203706e-7
flux1:, 0.102020202020202, 1.77218779386503e-6, -7.37944901819701e-7
flux1:, 0.103030303030303, 2.63090852112034e-6, -1.11118350510327e-6
flux1:, ...


This is comma-delimited data, which can easily be imported into any spreadsheet or plotting program (e.g. Octave/Matlab): the first column is the frequency, the second is the transmitted power, and the third is the reflected power.

We need to run the simulation twice, once with no-bend?=true and once with no-bend?=false (the default):

unix% meep no-bend?=true bend-flux.ctl | tee bend0.out
unix% meep bend-flux.ctl | tee bend.out


The tee command is a useful Unix command that saves the output to a file and displays it on the screen, so that we can see what is going on as it runs. Then, we should pull out the flux1 lines into a separate file to import them into our plotting program:

unix% grep flux1: bend0.out > bend0.dat
unix% grep flux1: bend.out > bend.dat


We import them to Octave/Matlab (using its dlmread command), and plot the results:

What are we plotting here? The transmittance is the transmitted flux (second column of bend.dat) divided by the incident flux (second column of bend0.dat), to give us the fraction of power transmitted. The reflectance is the reflected flux (third column of bend.dat) divided by the incident flux (second column of bend0.dat). We also have to multiply by -1 because all fluxes in Meep are computed in the positive-coordinate direction by default, and we want the flux in the $-x$ direction. Finally, the scattered loss is simply $1-transmittance-reflectance$.

We should also check whether our data is converged. We can do this by increasing the resolution and cell size and seeing by how much the numbers change. In this case, we'll try doubling the cell size:

unix% meep sx=32 sy=64 no-bend?=true bend-flux.ctl |tee bend0-big.out
unix% meep sx=32 sy=64 bend-flux.ctl |tee bend-big.out


Again, we must run both simulations in order to get the normalization right. The results are included in the plot above as dotted lines — you can see that the numbers have changed slightly for transmittance and loss, probably stemming from interference between light radiated directly from the source and light propagating around the waveguide.

## Angular Reflectance Spectrum of a Planar Interface

We turn to a similar but slightly different example for which there exists an analytic solution via the Fresnel equations: computing the broadband reflectance spectrum of a planar air-dielectric interface for an incident planewave over a range of angles. Similar to the previous example, we will need to run two simulations: (1) an empty cell with air/vacuum (n=1) everywhere to obtain the incident flux, and (2) with the dielectric (n=3.5) interface to obtain the reflected flux. For each angle of the incident planewave, a separate simulation is necessary.

A 1d cell must be used since a higher-dimensional cell will introduce artificial modes due to band folding. We will use a Gaussian source spanning visible wavelengths of 0.4 to 0.8 μm. Unlike a continuous-wave (CW) source, a pulsed source turns off. This enables a termination condition of when there are no fields remaining in the cell (due to absorption by the PMLs) via the run function stop-when-fields-decayed, similar to the previous example.

Creating an oblique planewave source typically requires specifying two parameters: (1) for periodic structures, the Bloch-periodic wavevector $\vec{k}$ via k-point, and (2) the source amplitude function amp-func for setting the $e^{i\vec{k} \cdot \vec{r}}$ spatial dependence ($\vec{r}$ is the position vector). Since we have a 1d cell and the source is at a single point, it is not necessary to specify the source amplitude (see this 2d example for how this is done). The magnitude of the Bloch-periodic wavevector is specified according to the dispersion relation formula for a planewave in homogeneous media with index n: $\omega=c|\vec{k}|/n$. As the source in this example is incident from air, $|\vec{k}|$ is simply equal to the frequency ω (the minimum frequency of the pulse which excludes the 2π factor). Note that a fixed wavevector only applies to a single frequency. Any broadband source is therefore incident at a specified angle for only a single frequency. This is described in more detail in Section 4.5 ("Efficient Frequency-Angle Coverage") in Chapter 4 ("Electromagnetic Wave Source Conditions") of the book Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology.

In this example, the plane of incidence which contains $\vec{k}$ and the surface normal vector is $xz$. The source angle θ is defined in degrees in the counterclockwise (CCW) direction around the $y$ axis with 0 degrees along the +$z$ axis. In Meep, a 1d cell is defined along the $z$ direction. When $\vec{k}$ is not set, only the Ex and Hy field components are permitted. A non-zero $\vec{k}$ results in a 3d simulation where all field components are allowed and are complex (the fields are real, by default). A current source with Ex polarization lies in the plane of incidence and corresponds to the convention of $\mathcal{P}$-polarization. In order to model the $\mathcal{S}$-polarization, we must use an Ey source. This example involves just the $\mathcal{P}$-polarization.

The simulation script is examples/refl-angular.ctl

(set-param! resolution 200)         ; pixels/um

(define-param dpml 1)               ; PML thickness
(define-param sz 10)                ; size of cell (without PMLs)
(set! sz (+ sz (* 2 dpml)))
(set! pml-layers (list (make pml (thickness dpml))))

(set! geometry-lattice (make lattice (size no-size no-size sz)))

(define-param wvl-min 0.4)          ; minimum wavelength of source
(define-param wvl-max 0.8)          ; maximum wavelength of source
(define fmin (/ wvl-max))           ; minimum frequency of source
(define fmax (/ wvl-min))           ; maximum frequency of source
(define fcen (* 0.5 (+ fmin fmax))) ; center frequency of source
(define df (- fmax fmin))           ; frequency width of source
(define-param nfreq 50)             ; number of frequency bins

; rotation angle (in degrees) of source: CCW around Y axis, 0 degrees along +Z axis
(define-param theta 0)

; if source is at normal incidence, force number of dimensions to be 1
(set! dimensions (if (= theta-r 0) 1 3))

; plane of incidence is xz
(set! k-point (vector3* fmin (vector3 (sin theta-r) 0 (cos theta-r))))

(set! sources (list (make source (src (make gaussian-src (frequency fcen) (fwidth df)))
(component Ex) (center 0 0 (+ (* -0.5 sz) dpml)))))

(define-param empty? true)

; add a block with n=3.5 for the air-dielectric interface
(if (not empty?)
(set! geometry (list (make block (size infinity infinity (* 0.5 sz)) (center 0 0 (* 0.25 sz)) (material (make medium (index 3.5)))))))

(define refl (add-flux fcen df nfreq (make flux-region (center 0 0 (* -0.25 sz)))))

(if (not empty?) (load-minus-flux "refl-flux" refl))

(run-sources+ (stop-when-fields-decayed 50 Ex (vector3 0 0 (+ (* -0.5 sz) dpml)) 1e-9))

(if empty? (save-flux "refl-flux" refl))

(display-fluxes refl)


The simulation script above computes and prints to standard output the reflectance at each frequency. Also included in the output is the wavevector component $k_x$ and the corresponding angle for the ($k_x$, ω) pair. For those frequencies not equal to the minimum frequency of the source, this is not the same as the specified angle of the incident planewave, but rather sin-1(kx/ω).

The following Bash shell script runs the simulation for the wavelength range of 0$^\circ$ to 80$^\circ$ in increments of 5$^\circ$. For each run, the script pipes the output to one file and extracts the reflectance data to a different file.

#!/bin/bash

for i in seq 0 5 80; do
meep empty?=true theta=${i} refl-angular.ctl |tee -a flux0_t${i}.out;
grep flux1: flux0_t${i}.out |cut -d , -f2- > flux0_t${i}.dat
meep empty?=false theta=${i} refl-angular.ctl |tee -a flux_t${i}.out;
grep flux1: flux_t${i}.out |cut -d , -f2- > flux_t${i}.dat
done


Two-dimensional plots of the angular reflectance spectrum based on the simulated data and the analytic Fresnel equations are generated using the Octave/Matlab script below. The plots are shown in the accompanying figure with four insets. The top left inset shows the simulated and analytic reflectance spectra at a wavelength of 0.6 μm. The top right inset shows the simulated reflectance spectrum as a function of the wavelength λ and wavevector $k_x$: $R(\lambda, k_x)$. The lower left inset is a transformation of $R(\lambda, k_x)$ into $R(\lambda, \theta)$. Note how the range of angles depends on the wavelength. For a particular angle, the reflectance is a constant for all wavelengths due to the dispersionless dielectric. The lower right inset is the analytic reflectance spectrum computed using the Fresnel equations. There is agreement between the simulated and analytic results. The Brewster's angle, where the transmittance is 1 and the reflectance 0, is tan-1(3.5/1)=74.1°. This is also verified by the simulated results.

In order to generate results for the missing portion of the reflectance spectrum (i.e., the white region), we will need to rerun the simulations for different wavelength spectra.

theta_in = [0:5:80];
Rmeep = [];
for j = 1:length(theta_in)
Rmeep = [Rmeep -f(:,2)./f0(:,2)];
endfor

freqs = f(:,1);
% convert frequency to wavelength
wvl = 1./freqs;
% create a 2d matrix for the wavelength by repeating the column vector for each angle
wvls = repmat(wvl,1,length(theta_in));

wvl_min = 0.4;
wvl_max = 0.8;
fcen = (1/wvl_min+1/wvl_max)/2;
kx = fcen*sind(theta_in);
kxs = repmat(kx,length(wvl),1);
thetas = asind(kxs./freqs);

figure;
pcolor(kxs,wvls,Rmeep);
shading interp; c = colormap("hot"); colormap(c); colorbar;
eval(sprintf("axis([%0.2g %0.2g %0.2g %0.2g])",kx(1),kx(end),min(wvl),max(wvl)));
eval(sprintf("set(gca, 'xtick', [%0.2g:0.2:%0.2g])",kx(1),kx(end)));
eval(sprintf("set(gca, 'ytick', [%0.1g:0.1:%0.1g])",wvl(end),wvl(1)));
xlabel("wavevector of Bloch-Periodic boundary condition (k_x/2π)");
ylabel("wavelength (μm)");
title("reflectance (meep)");

figure;
pcolor(thetas,wvls,Rmeep);
shading interp; c = colormap("hot"); colormap(c); colorbar;
eval(sprintf("axis([%0.2g %0.2g %0.2g %0.2g])",min(min(thetas)),max(max(thetas)),min(wvl),max(wvl)));
eval(sprintf("set(gca, 'xtick', [%d:20:%d])",min(min(thetas)),max(max(thetas))));
eval(sprintf("set(gca, 'ytick', [%0.1g:0.1:%0.1g])",wvl(end),wvl(1)));
xlabel("angle of incident planewave (degrees)");
ylabel("wavelength (μm)");
title("reflectance (meep)");

n1 = 1;
n2 = 3.5;

% compute angle of refracted planewave in medium n2
% for incident planewave in medium n1 at angle theta_in
theta_out = @(theta_in) asin(n1*sin(theta_in)/n2);

% compute Fresnel reflectance for P-polarization in medium n2
% for incident planewave in medium n1 at angle theta_in
R_fresnel = @(theta_in) abs((n1*cos(theta_out(theta_in))-n2*cos(theta_in))./(n1*cos(theta_out(theta_in))+n2*cos(theta_in))).^2;

Ranalytic = R_fresnel(thetas*pi/180);

figure;
pcolor(thetas,wvls,Ranalytic);
shading interp; c = colormap("hot"); colormap(c); colorbar;
eval(sprintf("axis([%0.2g %0.2g %0.2g %0.2g])",min(min(thetas)),max(max(thetas)),min(wvl),max(wvl)));
eval(sprintf("set(gca, 'xtick', [%d:20:%d])",min(min(thetas)),max(max(thetas))));
eval(sprintf("set(gca, 'ytick', [%0.1g:0.1:%0.1g])",wvl(end),wvl(1)));
xlabel("angle of incident planewave (degrees)");
ylabel("wavelength (μm)");
title("reflectance (analytic)");


## Modes of a Ring Resonator

As described in Introduction, another common task for FDTD simulation is to find the resonant modes — frequencies and decay rates — of some cavity structure. You might want to read that again to recall the basic simulation strategy. We will show how this works for a ring resonator, which is simply a waveguide bent into a circle. This script can be also found in examples/ring.ctl. In fact, since this structure has cylindrical symmetry, we can simulate it much more efficiently by using cylindrical coordinates, but for illustration here we'll just use an ordinary 2d simulation.

As before, we'll define some parameters to describe the geometry, so that we can easily change the structure:

(define-param n 3.4)                     ; index of waveguide
(define-param w 1)                       ; width of waveguide
(define-param r 1)                       ; inner radius of ring
(define-param dpml 2)                    ; thickness of PML
(define sxy (* 2 (+ r w pad dpml)))      ; cell size
(set! geometry-lattice (make lattice (size sxy sxy no-size)))


How do we make a circular waveguide? So far, we've only seen block objects, but Meep also lets you specify cylinders, spheres, ellipsoids, and cones, as well as user-specified dielectric functions. In this case, we'll use two cylinder objects, one inside the other:

(set! geometry (list
(make cylinder (center 0 0) (height infinity)
(radius (+ r w)) (material (make medium (index n))))
(make cylinder (center 0 0) (height infinity)
(set! pml-layers (list (make pml (thickness dpml))))
(set-param! resolution 10)


Later objects in the geometry list take precedence over (lie "on top of") earlier objects, so the second air (ε=1) cylinder cuts a circular hole out of the larger cylinder, leaving a ring of width w.

We don't know the frequency of the mode(s) ahead of time, so we'll just hit the structure with a broad Gaussian pulse to excite all of the Ez-polarized modes in a chosen bandwidth:

(define-param fcen 0.15)   ; pulse center frequency
(define-param df 0.1)      ; pulse frequency width
(set! sources (list
(make source
(src (make gaussian-src (frequency fcen) (fwidth df)))
(component Ez) (center (+ r 0.1) 0))))


Finally, we are ready to run the simulation. The basic idea is to run until the sources are finished, and then to run for some additional period of time. In that additional period, we'll perform some signal processing on the fields at some point with Harminv to identify the frequencies and decay rates of the modes that were excited:

(run-sources+ 300
(at-beginning output-epsilon)
(after-sources (harminv Ez (vector3 (+ r 0.1)) fcen df)))


The signal processing is performed by the harminv function, which takes four arguments: the field component Ez and position ($r$+0.1,0) to analyze, and a frequency range given by a center frequency and bandwidth (same as the source pulse). Note that we wrap harminv in (after-sources ...), since we only want to analyze the frequencies in the source-free system (the presence of a source will distort the analysis). At the end of the run, harminv prints a series of lines (beginning with harminv0:, to make it easy to grep for) listing the frequencies it found:

harminv0:, frequency, imag. freq., Q, |amp|, amplitude, error
harminv0:, 0.118101575043663, -7.31885828253851e-4, 80.683059081382, 0.00341388964904578, -0.00305022905294175-0.00153321402956404i, 1.02581433904604e-5
harminv0:, 0.147162555528154, -2.32636643253225e-4, 316.29272471914, 0.0286457663908165, 0.0193127882016469-0.0211564681361413i, 7.32532621851082e-7
harminv0:, 0.175246750722663, -5.22349801171605e-5, 1677.48461212767, 0.00721133215656089, -8.12770506086109e-4-0.00716538314235085i, 1.82066436470489e-7


There are six columns (in addition to the label), comma-delimited for easy import into other programs. The meaning of these columns is as follows. Harminv analyzes the fields $f(t)$ at the given point, and expresses this as a sum of modes (in the specified bandwidth):

for complex amplitudes $a_n$ and complex frequencies ω$_n$. The six columns relate to these quantities. The first column is the real part of ω$_n$, expressed in our usual 2πc units, and the second column is the imaginary part — a negative imaginary part corresponds to an exponential decay. This decay rate, for a cavity, is more often expressed as a dimensionless "lifetime" $Q$, defined by:

$Q$ is the number of optical periods for the energy to decay by $\exp(-2\pi)$, and 1/$Q$ is the fractional bandwidth at half-maximum of the resonance peak in Fourier domain. This $Q$ is the third column of the output. The fourth and fifth columns are the absolute value $|a_n|$ and complex amplitudes $a_n$. The last column is a crude measure of the error in the frequency (both real and imaginary). If the error is much larger than the imaginary part, for example, then you can't trust the $Q$ to be accurate. Note: this error is only the uncertainty in the signal processing, and tells you nothing about the errors from finite resolution, finite cell size, and so on.

An interesting question is how long should we run the simulation, after the sources are turned off, in order to analyze the frequencies. With traditional Fourier analysis, the time would be proportional to the frequency resolution required, but with harminv the time is much shorter. For example, there are three modes. The last has a $Q$ of 1677, which means that the mode decays for about 2000 periods or about 2000/0.175 = 104 time units. We have only analyzed it for about 300 time units, however, and the estimated uncertainty in the frequency is 10-7 (with an actual error of about 10-6, from below). In general, you need to increase the run time to get more accuracy, and to find very high $Q$ values, but not by much. In some cases, modes with $Q$ of around 109 can be found with only 200 periods.

In this case, we found three modes in the specified bandwidth, at frequencies of 0.118, 0.147, and 0.175, with corresponding $Q$ values of 81, 316, and 1677. As was shown by Marcatilli in 1969, the $Q$ of a ring resonator increases exponentially with the product of ω and ring radius. Suppose that we want to actually see the field patterns of these modes. No problem: we just re-run the simulation with a narrow-band source around each mode and output the field at the end.

In particular, to output the field at the end we might add an (at-end output-efield-z) argument to our run-sources+ function, but this is problematic: we might be unlucky and output at a time when the Ez field is almost zero (i.e. when all of the energy is in the magnetic field), in which case the picture will be deceptive. Instead, at the end of the run we'll output 20 field snapshots over a whole period 1/fcen by appending the command:

(run-until (/ 1 fcen) (at-every (/ 1 fcen 20) output-efield-z))


We can get our modes just by running e.g.:

unix% meep fcen=0.118 df=0.01 ring.ctl


After each one of these commands, we'll convert the fields into PNG images and thence into an animated GIF (as with the bend movie, above), via:

unix% h5topng -RZc dkbluered -C ring-eps-000000.00.h5 ring-ez-*.h5
unix% convert ring-ez-*.png ring-ez-0.118.gif


The resulting animations for (from left to right) 0.118, 0.147, and 0.175, are below, in which you can clearly see the radiating fields that produce the losses:

Each of these modes is, of course, doubly-degenerate according to the representations of the $C_{\infty\mathrm{v}}$ symmetry group. The other mode is simply a slight rotation of this mode to make it odd through the $x$ axis, whereas we excited only the even modes due to our source symmetry. Equivalently, one can form clockwise and counter-clockwise propagating modes by taking linear combinations of the even/odd modes, corresponding to an angular $\phi$ dependence $e^{\pm i m\phi}$ for m=3, 4, and 5 in this case.

You may have noticed, by the way, that when you run with the narrow-bandwidth source, harminv gives you slightly different frequency and $Q$ estimates, with a much smaller error estimate — this is not too strange, since by exciting a single mode you generate a cleaner signal that can be analyzed more accurately. For example, the narrow-bandwidth source for the ω=0.175 mode gives:

harminv0:, 0.175247426698716, -5.20844416909221e-5, 1682.33949533974, 0.185515412838043, 0.127625313330642-0.13463932485617i, 7.35320734698267e-12


which differs by about 10-6 from the earlier estimate; the difference in $Q$ is, of course, larger because a small absolute error in ω gives a larger relative error in the small imaginary frequency.

### Exploiting Symmetry

In this case, because we have a mirror symmetry plane (the $x$ axis) that preserves both the structure and the sources, we can exploit this mirror symmetry to speed up the computation. See also Exploiting Symmetry. In particular, everything about the input file is the same except that we add a single line, right after we specify the sources:

(set! symmetries (list (make mirror-sym (direction Y))))


This tells Meep to exploit a mirror-symmetry plane through the origin perpendicular to the $y$ direction. Meep does not check whether your system really has this symmetry — you should only specify symmetries that really preserve your structure and your sources.

Everything else about your simulation is the same: you can still get the fields at any point, the output file still covers the whole ring, and the harminv outputs are exactly the same. Internally, however, Meep is only doing computations with half of the structure, and the simulation is around twice as fast.

In general, the symmetry of the sources may require some phase. For example, if our source was in the $y$ direction instead of the $z$ direction, then the source would be odd under mirror flips through the $x$ axis. We would specify this by (make mirror-sym (direction Y) (phase -1)). See User Interface for more symmetry possibilities.

In this case, we actually have a lot more symmetry that we could potentially exploit, if we are willing to restrict the symmetry of our source/fields to a particular angular momentum (i.e. angular dependence $e^{im\phi}$). See also Tutorial/Ring Resonator in Cylindrical Coordinates for how to solve for modes of this cylindrical geometry much more efficiently.

## Visualizing 3d Structures

The previous examples were based on a 1d or 2d cell in which the structures and fields can be visualized using h5topng of the h5utils package. In order to visualize 3d structures, you can use Mayavi. The following example, which includes a simulation script and shell commands, involves a sphere with index 3.5 perforated by a conical hole. There are no other simulation parameters specified. The permittivity data is written to an HDF5 file using output-epsilon. The HDF5 data is then converted to VTK using h5tovtk. VTK data can be visualized using Mayavi or Paraview via the IsoSurface module.

(set-param! resolution 50)

(set! geometry-lattice (make lattice (size 3 3 3)))

(set! geometry (list (make sphere (radius 1) (material (make medium (index 3.5))) (center 0 0 0))
(make cone (radius 0.8) (radius2 0.1) (height 2) (material air) (center 0 0 0))))

(init-fields)

(output-epsilon)

(exit)

#!/bin/bash

meep sphere-cone.ctl;

h5tovtk -o epsilon.vtk structure_demo-eps-000000.00.h5;

mayavi2 -d epsilon.vtk -m IsoSurface &> /dev/null &



## Editors and ctl

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

 (push '("\\.ctl\\'" . scheme-mode) auto-mode-alist)


or if your emacs version is 24.3 or earlier and you have other ".ctl" files which are not Scheme:

 (if (assoc "\\.ctl" auto-mode-alist)
nil

Incidentally, emacs scripts are written in "elisp," a language closely related to Scheme.