How to Make a Meep Have a Baby

From AbInitio

(Redirected from Meep tutorial)

Jump to: navigation, search

Meep

Download

Release notes

FAQ

Meep manual

Introduction

Installation

Tutorial

Reference

C++ Tutorial

C++ Reference

Acknowledgements

License and Copyright

In this page, nosotros'll go through a couple of simple examples that illustrate the process of computing fields, manual/reflection spectra, and resonant modes in Meep. All of the examples here are 2-dimensional calculations, simply because they are quicker than 3d computations and they illustrate most of the essential features. For more advanced functionality involving 3d computations, encounter the Simpetus projects page.

This tutorial uses the libctl/Scheme scripting interface to Meep, which is what we look nearly users to apply most of the time. There is as well a C++ interface that may give additional flexibility in some situations; that is described in the C++ tutorial.

In social club to convert the HDF5 output files of Meep into images of the fields and and then on, this tutorial uses our free h5utils programs. Y'all could as well use whatsoever other programme, such as Matlab, that supports reading HDF5 files.

Contents

• i The ctl file
• 2 Fields in a waveguide
  • ii.1 A directly waveguide
  • 2.two A 90° bend
    • 2.2.1 Output tips and tricks
• 3 Transmission spectrum effectually a waveguide curve
• iv Modes of a ring resonator
  • 4.1 Exploiting symmetry
• v More examples
• half-dozen Editors and ctl

The ctl file

The use of Meep revolves around the control file, abbreviated "ctl" and typically chosen something like foo.ctl (although you can use whatever file proper name you wish). The ctl file specifies the geometry you wish to report, the current sources, the outputs computed, and everything else specific to your calculation. Rather than a flat, inflexible file format, all the same, the ctl file is actually written in a scripting language. This means that it can be everything from a elementary sequence of commands setting the geometry, etcetera, to a full-fledged program with user input, loops, and annihilation else that yous might need.

Don't worry, though—simple things are simple (you don't need to be a Real Programmer), and even there yous volition appreciate the flexibility that a scripting language gives yous. (east.grand. you can input things in whatsoever order, without regard for whitespace, insert comments where you delight, omit things when reasonable defaults are available...)

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

• Scheme, a powerful and cute programming language adult at MIT, which has a especially simple syntax: all statements are of the form (function arguments...). Nosotros run Scheme nether the GNU Guile interpreter (designed to be plugged into programs as a scripting and extension language). Y'all don't need to know much Scheme for a basic ctl file, simply it is ever there if y'all need it; you can learn more about information technology from these Guile and Scheme links.
• libctl, a library that nosotros built on height of Guile to simplify communication between Scheme and scientific ciphering software. libctl sets the basic tone of the interface and defines a number of useful functions (such as multi-variable optimization, numeric integration, and then on). Come across the libctl manual pages.
• 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, delight have a moment to foliage through the libctl tutorial to get a feel for the basic style of the interface, before we get to the Meep-specific stuff below. (If you've used MPB, all of this stuff should already be familiar, although Meep is somewhat more circuitous because it can perform a wider variety of computations.)

Okay, let'south go on with our tutorial. The Meep program is normally invoked by running something like the post-obit at the Unix command-line (herein denoted past the unix% prompt):

unix% meep          foo.ctl          >&          foo.out        

which reads the ctl file foo.ctl and executes it, saving the output to the file foo.out. However, if you lot invoke meep with no arguments, you are dropped into an interactive mode in which you tin blazon commands and run across their results immediately. If you do that now, y'all can paste in the commands from the tutorial every bit you follow information technology and see what they do.

Fields in a waveguide

For our first case, let'due south examine the field pattern excited by a localized CW source in a waveguide— offset directly, then bent. Our waveguide volition have (not-dispersive) and width 1. That is, we selection units of length then that the width is 1, and ascertain everything in terms of that (see also units in meep).

A straight waveguide

Before nosotros define the structure, however, nosotros have to define the computational cell. We're going to put a source at ane terminate and watch it propagate downwardly the waveguide in the x direction, then let'southward utilize a cell of length 16 in the ten direction to give information technology some distance to propagate. In the y management, we merely need plenty room so that the boundaries (beneath) don't affect the waveguide style; allow's requite information technology a size of viii. We at present specify these sizes in our ctl file via the geometry-lattice variable:

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

(The proper noun geometry-lattice comes from MPB, where information technology tin exist used to define a more general periodic lattice. Although Meep supports periodic structures, it is less full general than MPB in that affine grids are not supported.) ready! is a Scheme command to set the value of an input variable. The terminal no-size parameter says that the computational cell has no size in the z direction, i.e. it is two-dimensional.

Now, nosotros can add the waveguide. About commonly, the construction is specified by a listing of geometric objects, stored in the geometry variable. Here, nosotros do:

(set up! geometry (list                 (make cake (center 0 0) (size infinity 1 infinity)                       (fabric (make dielectric (epsilon 12))))))        

Dielectric function (black = high, white = air), for directly waveguide simulation.

The waveguide is specified by a block (parallelepiped) of size , with ε=12, centered at (0,0) (the centre of the computational cell). Past default, any place where there are no objects in that location is air (ε=1), although this can be inverse by setting the default-material variable. The resulting structure is shown at right.

At present that we have the structure, nosotros need to specify the current sources, which is specified as a listing called sources of source objects. The simplest affair is to add together a point source J _z :

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

Here, we gave the source a frequency of 0.15, and specified a continuous-src which is just a fixed-frequency sinusoid exp( − iωt) that (past 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 vacuum wavelength. Thus, 0.15 corresponds to a vacuum wavelength of about one / 0.xv = 6.67, or a wavelength of nigh 2 in the cloth—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 one/2√11 or roughly 0.15076.) Notation also that to specify a J _z , we specify a component Due east z (e.g. if nosotros wanted a magnetic current, we would specify Hx, Hy, or Hz). The current is located at ( − 7,0), which is one unit to the right of the left edge of the cell—we ever desire to get out a piffling infinite betwixt sources and the prison cell boundaries, to go along the boundary weather from interfering with them.

Speaking of boundary weather, we want to add absorbing boundaries around our prison cell. Arresting boundaries in Meep are handled past perfectly matched layers (PML)— which aren't really a boundary condition at all, but rather a fictitious absorbing textile added effectually the edges of the cell. To add together an absorbing layer of thickness 1 around all sides of the cell, we exercise:

(set! pml-layers (list (brand pml (thickness ane.0))))        

pml-layers is a list of pml objects—you may take more than one pml object if you desire PML layers only on certain sides of the cell, e.g. (make pml (thickness i.0) (management X) (side High)) specifies a PML layer on but the + x side. Now, we note an important signal: the PML layer is inside the jail cell, overlapping whatever objects y'all have at that place. So, in this example 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; see perfectly matched layers for more information.

Meep volition discretize this construction in space and time, and that is specified by a single variable, resolution, that gives the number of pixels per distance unit of measurement. We'll set this resolution to 10, which corresponds to effectually 67 pixels/wavelength, or effectually 20 pixels/wavelength in the high-dielectric fabric. (In general, at least 8 pixels/wavelength in the highest dielectric is a skillful idea.) This will give us a jail cell.

(set! resolution 10)        

At present, we are ready to run the simulation! Nosotros do this by calling the run-until role. 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-showtime output-epsilon)            (at-end output-efield-z))        

Hither, we are outputting the dielectric role ε 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-stop, which practice simply what they say. There are several other such functions to modify the output beliefs—and you can, of class, write your own, and in fact yous can do whatsoever ciphering or output you desire at any time during the time evolution (and fifty-fifty change the simulation while information technology is running).

Information technology should consummate in a few seconds. If you lot are running interactively, the two output files will be called eps-000000.00.h5 and ez-000200.00.h5 (observe that the file names include the time at which they were output). If nosotros were running a tutorial.ctl file, and then the outputs will be tutorial-eps-000000.00.h5 and tutorial-ez-000200.00.h5. In any instance, we can at present analyze and visualize these files with a wide variety of programs 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 (then that it is around 450 pixels broad, in this case). In fact, precisely this command is what created the dielectric image above. Much more than 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 colour scale become from night bluish (negative) to white (zero) to dark red (positive), and the -a/-A options overlay the dielectric part as light gray contours. This results in the image:

Here, we see that the the source has excited the waveguide mode, but has as well excited radiating fields propagating away from the waveguide. At the boundaries, the field quickly goes to zip due to the PML layers. If we look carefully (click on the image to see a larger view), we run into somethinge else—the epitome is "speckled" towards the right side. This is because, by turning on the current abruptly at t = 0, we take excited high-frequency components (very loftier club modes), and nosotros take non waited long enough for them to die away; we'll eliminate these in the side by side department by turning on the source more smoothly.

A 90° bend

Now, nosotros'll get-go a new simulation where we await at the fields in a aptitude waveguide, and we'll exercise a couple of other things differently besides. If you are running Meep interactively, you will want to get rid of the one-time structure and fields then that Meep volition re-initialize them:

(reset-meep)        

So permit'due south set up the bent waveguide, in a slightly bigger computational prison cell, via:

(set up! geometry-lattice (make lattice (size sixteen 16 no-size)))  (fix! geometry (listing                 (make block (center -2 -3.five) (size 12 i infinity)                       (material (make dielectric (epsilon 12))))                 (brand block (center 3.5 two) (size 1 12 infinity)                       (textile (make dielectric (epsilon 12))))))  (set! pml-layers (list (make pml (thickness one.0)))) (ready! resolution 10)        

Bent waveguide dielectric function and coordinate system.

Notation that we at present have ii blocks, both off-heart to produce the bent waveguide structure pictured at correct. Equally illustrated in the effigy, the origin (0,0) of the coordinate system is at the center of the computational jail cell, with positive y being downwards in h5topng, and thus the block of size 12×ane is centered at ( − 2, − iii.5). Also shown in green is the source plane at x = − seven (see below).

We also need to shift our source to y = − 3.5 so that it is even so inside the waveguide. While we're at it, we'll brand a couple of other changes. Beginning, 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 past adding a size property to the source (Meep also has an eigenmode source feature which can be used here and is covered in a separate tutorial). Second, instead of turning the source on all of a sudden at t = 0 (which excites many other frequencies because of the discontinuity), we volition ramp it on slowly (technically, Meep uses a tanh turn-on role) over a time proportional to the width of xx time units (a little over 3 periods). Finally, just for variety, nosotros'll specify the (vacuum) wavelength instead of the frequency; again, we'll use a wavelength such that the waveguide is half a wavelength broad.

(set! sources (list                (make source                  (src (brand continuous-src                         (wavelength (* two (sqrt 12))) (width 20)))                  (component Ez)                  (center -7 -3.v) (size 0 1))))        

Finally, nosotros'll run the simulation. Instead of running output-efield-z only at the cease of the simulation, still, we'll run it at every 0.6 fourth dimension units (about ten times per period) via (at-every 0.six output-efield-z). Past itself, this would output a split file for every different output fourth dimension, simply instead nosotros'll use another feature of Meep to output to a unmarried 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)))        

Hither, "ez" determines the proper 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 proper noun for a ctl file (eastward.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 unmarried dataset ez that is a 162×162×330 array, where the final dimension is time. (This is rather a large file, 69MB; later, we'll meet ways to reduce this size if we just want images.) Now, we accept a number of choices of how to output the fields. To output a single fourth dimension slice, we tin can employ the aforementioned h5topng command equally earlier, but with an additional -t option to specify the time alphabetize: due east.chiliad. h5topng -t 229 will output the last fourth dimension slice, similar to before. Instead, let's create an animation of the fields equally a function of fourth dimension. First, we accept to create images for all of the time slices:

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

This is similar to the command before, with ii new options: -t 0:329 outputs images for all time indices from 0 to 329, i.e. all of the times, and the the -R flag tells h5topng to use a consistent colour scale for every image (instead of scaling each epitome independently). Then, we have to convert these images into an blitheness in some format. For this, nosotros'll use the free ImageMagick catechumen program (although there is other software that will practise the trick as well).

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

Hither, we are using an animated GIF format for the output, which is not the most efficient animation format (eastward.g. ez.mpg, for MPEG format, would be better), but information technology is unfortunately the only format supported past this Wiki software. This results in the post-obit animation :

x past fourth dimension piece of bent waveguide (vertical = time).

It is clear that the transmission effectually the bend is rather low for this frequency and structure—both large reflection and large radiations loss are conspicuously visible. Moreover, since we operating are just barely below the cutoff for single-manner beliefs, we are able to excite a second leaky mode after the waveguide bend, whose 2d-guild mode pattern (superimposed with the primal mode) is apparent in the blitheness. At right, we show a field snapshot from a simulation with a larger cell along the y management, in which y'all can see that the second-guild leaky mode decays abroad, leaving us with the fundamental way propagating downward.

Instead of doing an blitheness, another interesting possibility is to make an paradigm from a slice. Hither is the y = − 3.five slice, which gives the states an image of the fields in the first waveguide branch as a office of time.

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

Here, the -0y -35 specifies the y = − 3.5 slice, where we accept multiplied by x (our resolution) to get the pixel coordinate.

Output tips and tricks

Above, we outputted the full 2d data slice at every 0.half dozen fourth dimension units, resulting in a 69MB file. This is non likewise bad by today's standards, just you lot can imagine how large 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 y'all know what yous want to look at.

To create the movie to a higher place, all we really need are the images corresponding to each time. Images tin be stored much more than efficiently than raw arrays of numbers—to exploit this fact, Meep allows you to output PNG images instead of HDF5 files. In item, 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 paradigm files). That is:

(run-until 200 (at-every 0.half-dozen (output-png Ez "-Zc bluered")))        

will output a PNG file file every 0.6 fourth dimension units, which can then exist combined with convert equally above to create a moving picture. The movie will exist like to the one before, but non identical because of how the color scale is determined. Before, we used the -R choice to make h5topng apply a uniform color scale for all images, based on the minimum/maximum field values over all fourth dimension steps. That is not possible, hither, because we output an image before knowing the field values at future time steps. Thus, what output-png does is to set its color calibration based on the minimum/maximum field values from all past times—therefore, the color scale will slowly "ramp upwards" 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 apply 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 past default filename-out/ if our ctl file is filename.ctl.

What if we want to output an slice, equally above? To do this, nosotros simply really wanted the values at y = − 3.5, and therefore we tin exploit another powerful Meep output feature—Meep allows us to output but a subset of the computational prison cell. This is washed using the in-volume office, which (like at-every and to-appended) is another function that modifies the behavior of other output functions. In particular, nosotros can do:

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

The first argument to in-volume is a volume, specified by (volume (heart ...) (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 162×330 corresponding to the desired slice.

Manual spectrum effectually a waveguide curve

In a higher place, nosotros computed the field patterns for light propagating effectually a waveguide curve. While this is pretty, the results are not quantitatively satisfying. We'd like to know exactly how much ability makes it around the bend, how much is reflected, and how much is radiated abroad. How tin can nosotros practise this?

The basic principles were described in the Meep introduction; please re-read that section if yous have forgotten. Basically, we'll tell Meep to keep track of the fields and their Fourier transforms in a certain region, and from this compute the flux of electromagnetic energy as a function of ω. Moreover, we'll become an unabridged spectrum of the transmission in a single run, by Fourier-transforming the response to a short pulse. Withal, in order to normalize the transmission (to get manual as a fraction of incident power), we'll have to do two runs, one with and one without a curve.

This command file volition be more complicated than before, so you'll definitely want it as a divide file rather than typing it interactively. Run into the bend-flux.ctl file included with Meep in its examples/ directory.

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

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

(Notice that a semicolon ";" begins a comment, which is ignored past Meep.) ascertain-param is a libctl feature to define variables that tin can be overridden from the command line. Nosotros could now practise meep sx=17 tut-wvg-bend-trans.ctl to change the X size to 17, without editing the ctl file, for case. 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 computational cell:

(define-param pad 4) ; padding distance betwixt waveguide and cell edge          (define-param due west 1) ; width of waveguide        

In order to define the waveguide positions, etcetera, we will now take to use arithmetic. For example, the y heart 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 (+ one two), and so on, so we will define the vertical and horizontal waveguide centers as:

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

Now, we have to make the geometry, as before. This time, withal, we really want 2 geometries: the bend, and also a directly waveguide for normalization. We could do this with two separate ctl files, only that is annoying. Instead, we'll define a parameter no-curve? which is true for the straight-waveguide case and faux for the bend.

(define-param no-bend? faux) ; if true, have straight waveguide, non bend        

At present, we define the geometry via 2 cases, with an if statement—the Scheme syntax is (if predicate? if-true if-false).

(set! geometry       (if no-bend?           (list            (brand block              (middle 0 wvg-ycen)              (size infinity w infinity)              (cloth (make dielectric (epsilon 12)))))           (list            (brand block              (eye (* -0.5 pad) wvg-ycen)              (size (- sx pad) w infinity)              (cloth (brand dielectric (epsilon 12))))            (make block              (center wvg-xcen (* 0.five pad))              (size west (- sy pad) infinity)              (material (make dielectric (epsilon 12)))))))        

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

The source is at present 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 every bit usual.

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

Find how nosotros're using our parameters similar wvg-ycen and w: if we alter the dimensions, everything volition now shift automatically. The boundary weather and resolution are gear up as earlier, except that now we'll use set-param! so that we can override the resolution from the control line.:

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

Finally, nosotros take to specify where nosotros want Meep to compute the flux spectra, and at what frequencies. (This must be washed subsequently specifying the geometry, sources, resolution, etcetera, because all of the field parameters are initialized when flux planes are created.)

(define-param nfreq 100) ; number of frequencies at which to compute flux              (define trans ; transmitted flux                                                       (add together-flux fcen df nfreq                 (if no-bend?                     (make flux-region                      (heart (- (/ sx 2) i.5) wvg-ycen) (size 0 (* w two)))                     (brand flux-region                      (center wvg-xcen (- (/ sy ii) 1.5)) (size (* w ii) 0))))) (define refl ; reflected flux                                                          (add-flux fcen df nfreq                  (brand flux-region                     (eye (+ (* -0.v sx) 1.v) 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 from the purlieus of the cell, so that they do not lie inside the absorbing PML regions.) Again, in that location are two cases: the transmitted flux is either computed at the correct or the lesser of the computational jail cell, depending on whether the waveguide is straight or bent.

Here, the fluxes will be computed for 100 (nfreq) 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, to far exterior the pulse bandwidth, the spectral ability is so low that numerical errors make the computed fluxes useless.

At present, every bit described in the Meep introduction, calculating reflection spectra is a bit catchy considering nosotros need to separate the incident and reflected fields. We exercise 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, nosotros might decrease these after the run, simply it turns out to be more convenient to decrease the incident fields start and then accrue 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). Nosotros can call them every bit follows:

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

This uses a file chosen refl-flux.h5, or actually bend-flux-refl-flux.h5 (the ctl file name is used as a prefix) to shop/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 washed automatically once information technology has decayed for a few standard deviations), plus an additional 500 time units.

Why practice we proceed running after the source has turned off? Because we must give the pulse fourth dimension to propagate completely across the prison cell. Moreover, the time required is a bit tricky to predict when you have circuitous structures, because there might exist 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 | E _z | ² at the finish of the waveguide must have rust-covered past a given corporeality (east.g. 1/one thousand) from its peak value. We tin can do this via:

(run-sources+   (terminate-when-fields-rust-covered l Ez                            (if no-bend?                                 (vector3 (- (/ sx 2) one.5) wvg-ycen)                                (vector3 wvg-xcen (- (/ sy 2) 1.5)))                            1e-3))        

terminate-when-fields-decayed takes four arguments: (finish-when-fields-decayed dT component pt decay-by). What information technology does is, after the sources have turned off, information technology keeps running for an additional dT fourth dimension units every time the given |component|ⁱⁱ at the given bespeak has not rust-covered by at least decay-by from its height value for all times within the previous dT. In this case, dT=50, the component is East _z , the point is at the center of the flux airplane at the stop of the waveguide, and decay-by=0.001. So, it keeps running for an additional 50 fourth dimension units until the square amplitude has decayed by 1/1000 from its pinnacle: 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 similar:

flux1:, 0.1, seven.91772317108475e-seven, -iii.16449591437196e-7 flux1:, 0.101010101010101, 1.18410865137737e-six, -4.85527604203706e-vii flux1:, 0.102020202020202, 1.77218779386503e-half-dozen, -7.37944901819701e-7 flux1:, 0.103030303030303, 2.63090852112034e-vi, -one.11118350510327e-6 flux1:, ...        

This is comma-delimited data, which tin can hands be imported into any spreadsheet or plotting programme (due east.g. Matlab): the first column is the frequency, the second is the transmitted ability, and the 3rd is the reflected power.

Now, nosotros 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 curve-flux.ctl | tee bend.out        

(The tee control is a useful Unix command that saves the output to a file and displays it on the screen, so that we can run into what is going on equally it runs.) So, 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 > curve.dat        

Now, we import them to Matlab (using its dlmread command), and plot the results:

What are we plotting hither? The transmission is the transmitted flux (2d column of bend.dat) divided by the incident flux (second column of bend0.dat), to give the states the fraction of power transmitted. The reflection is the reflected flux (third cavalcade of bend.dat) divided past the incident flux (second column of bend0.dat); we too take 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 management. Finally, the loss is merely 1 - transmission - reflection.

We should also check whether our information is converged, by increasing the resolution and cell size and seeing by how much the numbers change. In this instance, nosotros'll just try doubling the prison cell size:

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

Over again, nosotros must run both simulations in order to become the normalization right. The results are included in the plot above as dotted lines—you can encounter that the numbers accept changed slightly for manual and loss, probably stemming from interference between light radiated directly from the source and light propagating effectually the waveguide. To be actually confident, we should probably run the simulation again with an even bigger jail cell, merely we'll telephone call information technology plenty for this tutorial.

Modes of a ring resonator

Every bit described in the Meep introduction, another common task for FDTD simulation is to observe the resonant modes—frequencies and disuse rates—of some electromagnetic cavity structure. (You might desire to read that introduction again to recall the basic computational strategy.) Here, we volition testify how this works for perhaps the simplest example of a dielectric crenel: a band resonator, which is just a waveguide bent into a circle. (This can exist likewise institute in the examples/ring.ctl file included with Meep.) In fact, since this structure has cylindrical symmetry, we tin simulate information technology much more efficiently by using cylindrical coordinates, but for illustration here we'll only use an ordinary 2d simulation.

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

(ascertain-param n 3.4) ; index of waveguide                                        (define-param westward one) ; width of waveguide                                          (define-param r ane) ; inner radius of ring                                         (define-param pad iv) ; padding between waveguide and edge of PML                 (define-param dpml two) ; thickness of PML                                          (define sxy (* ii (+ r w pad dpml))) ; cell size                                  (set! geometry-lattice (make lattice (size sxy sxy no-size)))        

How do we make a round waveguide? And so far, nosotros've only seen cake objects, only Meep likewise lets you specify cylinders, spheres, ellipsoids, and cones, likewise every bit user-specified dielectric functions. In this case, nosotros'll use ii cylinder objects, i within the other:

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

Afterward objects in the geometry list take precedence over (lie "on top of") earlier objects, then the second air () cylinder cuts a round pigsty out of the larger cylinder, leaving a ring of width west.

Now, we don't know the frequency of the way(s) alee of time, and so nosotros'll just hit the structure with a broad Gaussian pulse to excite all of the (TM polarized) modes in a chosen bandwidth:

(define-param fcen 0.15) ; pulse centre frequency                                (define-param df 0.ane)  ; pulse width (in frequency)                              (ready! sources (list                (make source                  (src (brand gaussian-src (frequency fcen) (fwidth df)))                  (component Ez) (eye (+ r 0.ane) 0))))        

Finally, nosotros 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 fourth dimension. In that additional period, we'll perform some signal-processing on the fields at some indicate with harminv to identify the frequencies and decay rates of the modes that were excited:

(run-sources+ 300               (at-outset output-epsilon)               (after-sources (harminv Ez (vector3 (+ r 0.one)) fcen df)))        

The betoken-processing is performed past the harminv function, which takes four arguments: the field component (here Ez) and position (here (r + 0.1,0)) to analyze, and a frequency range given by a heart frequency and bandwidth (here, the same every bit the source pulse). Note that we wrap harminv in (after-sources ...), since we simply want to analyze the frequencies in the source-free system (the presence of a source volition distort the analysis). At the finish of the run, harminv prints a serial of lines (beginning with harminv0:, to make it easy to grep for) listing the frequencies it constitute:

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

There are six columns (in add-on to the label), comma-delimited for like shooting fish in a barrel import into other programs. The significant of these columns is as follows. Harminv analyzes the fields f(t) at the given signal, and expresses this as a sum of modes (in the specified bandwidth):

for complex amplitudes a _n and complex frequencies ω _n . The half dozen 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 role—a negative imaginary part corresponds to an exponential decay. This disuse charge per unit, for a crenel, is more often expressed equally a dimensionless "lifetime" Q , defined by:

( Q is the number of optical periods for the energy to decay by exp( − 2π), and 1 / Q is the fractional bandwidth at half-maximum of the resonance peak in Fourier domain.) This Q is the third cavalcade of the output. The fourth and fifth columns are the absolute value | a _{due north} | and complex amplitudes a _n . The last column is a crude mensurate 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 exist accurate. Note: this error is simply the doubt in the signal processing, and tells you nothing about the errors from finite resolution, finite cell size, so on!

An interesting question is how long should we run the simulation, afterwards 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. Hither, for example, there are three modes. The last has a Q of 1677, which means that the mode decays for most 2000 periods or almost 2000/0.175 = ten^four time units. We take simply analyzed it for about 300 fourth dimension units, however, and the estimated uncertainty in the frequency is 10 ^{− 7} (with an actual error of about x ^{− 6} , from below)! In general, you need to increment the run time to get more accurateness, and to find very high Q values, simply not past much—in our ain work, nosotros have successfully institute Q = 10⁹ modes past analyzing only 200 periods.

In this example, we found 3 modes in the specified bandwith, at frequencies of 0.118, 0.147, and 0.175, with corresponding Q values of 81, 316, and 1677. (Equally was shown by Marcatilli in 1969, the Q of a ring resonator increases exponentially with the product of ω and ring radius.) Now, suppose that we desire to actually see the field patterns of these modes. No problem: we just re-run the simulation with a narrow-band source effectually each mode and output the field at the terminate.

In item, to output the field at the stop nosotros might add an (at-end output-efield-z) argument to our run-sources+ function, but this is problematic: we might exist unlucky and output at a time when the E _z field is almost cypher (i.eastward. when all of the energy is in the magnetic field), in which case the motion-picture show will be deceptive. Instead, at the end of the run we'll output twenty field snapshots over a whole period 1/fcen by appending the control:

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

Now, nosotros can get our modes only by running east.thousand.:

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

After each one of these commands, we'll catechumen the fields into PNG images and thence into an animated GIF (as with the bend motion picture, higher up), via:

unix% h5topng -RZc dkbluered -C ring-eps-000000.00.h5 ring-ez-*.h5 unix% catechumen 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 meet the radiating fields that produce the losses:

(Each of these modes is, of course, doubly-degenerate co-ordinate to the representations of the symmetry group. The other mode is only a slight rotation of this way to get in 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 fifty-fifty/odd modes, corresponding an athwart φ dependence for thou = iii, 4, and 5 in this case.)

You may take noticed, by the fashion, that when yous run with the narrow-bandwidth source, harminv gives you slightly unlike frequency and Q estimates, with a much smaller error judge—this is not too foreign, since by heady a single mode you generate a cleaner betoken that can be analyzed more accurately. For example, the narrow-bandwidth source for the ω = 0.175 mode gives:

harminv0:, 0.175247426698716, -5.20844416909221e-v, 1682.33949533974, 0.185515412838043, 0.127625313330642-0.13463932485617i, vii.35320734698267e-12        

which differs by about 0.000001 (10 ^{− vi} ) from the earlier judge; the difference in Q is, of course, larger considering a small-scale absolute error in ω gives a larger relative error in the minor imaginary frequency.

Exploiting symmetry

In this case, considering we accept a mirror symmetry airplane (the x axis) that preserves both the structure and the sources, we tin exploit this mirror symmetry to speed up the computation. (See too exploiting symmetry in Meep.) In particular, everything about the input file is the aforementioned except that nosotros add a unmarried line, correct after nosotros 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—y'all should only specify symmetries that actually preserve your structure and your sources.

Everything else near your simulation is the aforementioned: yous 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 (YMMV).

In general, the symmetry of the sources may crave some stage. For example, if our source was in the y management instead of the z direction, then the source would exist odd nether mirror flips through the ten centrality. We would specify this by (make mirror-sym (management Y) (phase -one)). See the Meep reference for more than symmetry possibilities.

In this instance, we really have a lot more symmetry that we could potentially exploit, if we are willing to restrict the symmetry of our source/fields to a item angular momentum (i.east. athwart dependence e ^{i mφ} ). See likewise: Band resonator in cylindrical coordinates for how to solve for modes of this cylindrical geometry much more than efficiently.

More than examples

The examples higher up suffice to illustrate the nearly bones features of Meep. Yet, at that place are many more advanced features that have not been demonstrated here. So, we hope to add, over time, a sequence of examples that exhibit more than complicated structures and computational techniques. The ones we accept so far are listed below.

If you have a good case you would like to share, feel free to register for a user proper name and log in, and then add a link to a new folio for your example beneath. (Delight prefix the link with "Meep Tutorial/". Simple examples illustrating a specific concept are preferred.)

• Meep Tutorial/Ring resonator in cylindrical coordinates
• Meep Tutorial/Band diagram, resonant modes, and manual in a holey waveguide
• Meep Tutorial/Third harmonic generation (Kerr nonlinearity)
• Meep Tutorial/Material dispersion
• Meep Tutorial/Casimir forces
• Meep Tutorial/Local density of states (Purcell enhancement)
• Meep Tutorial/Optical forces (Maxwell stress tensor; besides includes an eigenmode-source case)
• Meep Tutorial/Near-to-far-field spectra (radiation pattern of an antenna)
• Meep Tutorial/Multilevel-diminutive susceptibility (saturable absorption/proceeds)

Editors and ctl

It is useful to accept emacs use its scheme-fashion for editing ctl files, so that hitting tab indents nicely, and so on. emacs does this automatically for files ending with ".scm"; to do information technology for files ending with ".ctl" too, add the following lines to your ~/.emacs file:

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

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

(if (assoc "\\.ctl" auto-mode-alist)     nil   (add-to-listing 'automobile-mode-alist '("\\.ctl\\'" . scheme-way))))        

(Incidentally, emacs scripts are written in "elisp," a linguistic communication closely related to Scheme.)

If you lot don't employ emacs (or derivatives such every bit Aquamacs), it would be skilful to find another editor that supports a Scheme mode. For example, jEdit is a costless/open-source cantankerous-platform editor with Scheme-syntax support. Another option is GNU gedit (for GNU/Linux and Unix); in fact, S. Hessam Moosavi Mehr has donated a hilighting fashion for Meep/MPB that specially highlights the Meep/MPB keywords.

floreshopil1968.blogspot.com

Source: http://ab-initio.mit.edu/wiki/index.php/Meep_tutorial

Share this post