Meep Reference

From AbInitio

Revision as of 15:47, 2 August 2006; Stevenj (Talk | contribs)
(diff) ←Older revision | Current revision | Newer revision→ (diff)
Jump to: navigation, search
Release notes
Meep manual
C++ Tutorial
C++ Reference
License and Copyright

Here, we document the features exposed to the user by the Meep package. We do not document the Scheme language or the functions provided by libctl (see also the libctl User Reference section of the libctl manual).

This page is simply a compact listing of the functions exposed by the interface; for a gentler introduction, see the Meep tutorial. Also, we note that this page is not, and probably never will be, a complete listing of all functions. In particular, because of the SWIG wrappers, every function in the C++ interface is accessible from Scheme, but not all of these functions are documented or intended for end users.

See also our parallel Meep instructions for parallel (MPI) machines.


Input Variables

These are global variables that you can set to control various parameters of the Meep computation. In brackets after each variable is the type of value that it should hold. (The classes, complex datatypes like geometric-object, are described in a later subsection. The basic datatypes, like integer, boolean, cnumber, and vector3, are defined by libctl.)

geometry [list of geometric-object class]
Specifies the geometric objects making up the structure being simulated. When objects overlap, later objects in the list take precedence. Defaults to no objects (empty list).
sources [list of source class]
Specifies the current sources to be present in the simulation; defaults to none.
symmetries [list of symmetry class]
Specifies the spatial (mirror/rotation) symmetries to exploit in the simulation (defaults to none). The symmetries must be obeyed by both the structure and by the sources. See also: Exploiting symmetry in Meep.
pml-layers [list of pml class]
Specifies the absorbing PML boundary layers to use; defaults to none.
default-material [material-type class]
Holds the default material that is used for points not in any object of the geometry list. Defaults to air (ε of 1).
geometry-lattice [lattice class]
Specifies the the size of the unit cell (which is centered on the origin of the coordinate system). Any sizes of no-size imply (effectively) a reduced-dimensionality calculation (but a 2d xy calculation is especially optimized); see dimensions below. Defaults to a cubic cell of unit size.
dimensions [integer]
Explicitly specifies the dimensionality of the simulation, if the value is less than 3. If the value is 3 (the default), then the dimensions are automatically reduced to 2 if possible when geometry-lattice size in the z direction is no-size. If dimensions is the special value of CYLINDRICAL, then cylindrical coordinates are used and the x and z dimensions are interpreted as r and z, respectively. If dimensions is 1, then the cell must be along the z direction and only Ex and Hy field components are permitted. If dimensions is 2, then the cell must be in the xy plane.
m [number]
For CYLINDRICAL simulations, specifies that the angular φ dependence of the fields is of the form eimφ (default is m=0). If the simulation cell includes the origin r = 0, then m must be an integer.
resolution [number or vector3]
Specifies the computational grid resolution, in pixels per distance unit. Defaults to 10.
k-point [false or vector3]
If false (the default), then the boundaries are perfect metallic (zero electric field). If a vector, then the boundaries are Bloch-periodic: the fields at one side are \exp(i\mathbf{k}\cdot\mathbf{R}) times the fields at the other side, separated by the lattice vector \mathbf{R}. The k-point vector is specified in Cartesian coordinates, in units of 2π/distance. (This is different from MPB, equivalent to taking MPB's k-points through the function reciprocal->cartesian.)
ensure-periodicity [boolean]
If true (the default), and if the boundary conditions are periodic (k-point is not false), then the geometric objects are automatically repeated periodically according to the lattice vectors (the size of the computational cell).
eps-averaging? [boolean]
If true, then sub-pixel averaging is used when initializing the dielectric function. Still somewhat slow and experimental, so the default is currently false.
force-complex-fields? [boolean]
By default, Meep runs its simulations with purely real fields whenever possible. It uses complex fields (which require twice the memory and computation) if the k-point is non-zero or if m is non-zero. However, by setting force-complex-fields? to true, Meep will always use complex fields. See also: Complex fields in Meep.
filename-prefix [string]
A string prepended to all output filenames. If "" (the default), then Meep uses the name of the current ctl file, with ".ctl" replaced by "-" (e.g. foo.ctl uses a "foo-" prefix). See also Output file names.
Courant [number]
Specify the Courant factor S which relates the time step size to the spatial discretization: cΔt = SΔx. Default is 0.5. For numerical stability, the Courant factor must be at most n_\textrm{min}/\sqrt{\textrm{\# dimensions}}, where nmin is the minimum refractive index (usually 1), and in practice S should be slightly smaller.
output-volume [meep::geometric_volume*]
Specifies the default region of space that is output by the HDF5 output functions (below); see also the (volume ...) function to create meep::geometric_volume* objects. The default is '() (null), which means that the whole computational cell is output. Normally, you should use the (in-volume ...) function to modify the output volume instead of setting output-volume directly.
output-single-precision? [boolean]
Meep performs its computations in double precision, and by default its output HDF5 files are in the same format. However, by setting this variable to true (default is false) you can instead output in single precision which saves a factor of two in space.
progress-interval [number]
Time interval (seconds) after which Meep prints a progress message; default is 4 seconds.

The require a bit more understanding of the inner workings of Meep to use (see also the SWIG wrappers).

structure [meep::structure*]
Pointer to the current structure being simulated; initialized by (init-structure) which is called automatically by (init-fields) which is called automatically by any of the (run) functions.
fields [meep::fields*]
Pointer to the current fields being simulated; initialized by (init-fields) which is called automatically by any of the (run) functions.
num-chunks [integer]
Minimum number of "chunks" (subarrays) to divide the structure/fields into (default 0); actual number is determined by number of processors, PML layers, etcetera. (Mainly useful for debugging.)

Predefined Variables

Variables predefined for your convenience and amusement.

air, vacuum [material-type class]
Two aliases for a predefined material type with a dielectric constant of 1.
metal [material-type class]
A predefined material type corresponding to a perfectly conducting metal (in which the electric field is zero).
nothing [material-type class]
A material that, effectively, punches a hole through other objects to the background (default-material).
infinity [number]
A big number (1.0e20) to use for "infinite" dimensions of objects.
pi [number]
π (3.14159...).

Constants (enumerated types)

Several of the functions/classes in Meep ask you to specify e.g. a field component or a direction in the grid. These should be one of the following constants:

direction constants
Specify a direction in the grid. One of X, Y, Z, R, P for: x, y, z, r, φ, respectively.
boundary-side constants
Specify particular boundary in some direction (e.g. + x or x). One of Low or High.
component constants
Specify a particular field (or other) component. One of Ex, Ey, Ez, Er, Ep, Hx, Hy, Hz, Hy, Hp, Hz, Dx, Dy, Dz, Dr, Dp, Dielectric, for: Ex, Ey, Ez, Er, Eφ, Hx, Hy, Hz, Hr, Hφ, Dx, Dy, Dz, Dr, Dφ, \varepsilon, respectively.
derived-component constants
These are additional components, which are not actually stored by Meep but are computed as needed, mainly for use in output functions. One of Sx, Sy, Sz, Sr, Sp, EnergyDensity, D-EnergyDensity, H-EnergyDensity for: Sx, Sy, Sz, Sr, Sφ (components of the Poynting vector \mathrm{Re}\,\mathbf{E}^* \times \mathbf{H}), (\mathbf{E}^* \cdot \mathbf{D} + |\mathbf{H}|^2)/2, \mathbf{E}^* \cdot \mathbf{D}/2, |\mathbf{H}|^2/2, respectively.


Classes are complex datatypes with various "properties" which may have default values. Classes can be "subclasses" of other classes; subclasses inherit all the properties of their superclass, and can be used any place the superclass is expected. An object of a class is constructed with:

(make class (prop1 val1) (prop2 val2) ...)

See also the libctl manual.

Meep defines several types of classes, the most numerous of which are the various geometric object classes (which are the same as those used in MPB. You can also get a list of the available classes, along with their property types and default values, at runtime with the (help) command.


The lattice class is normally used only for the geometry-lattice variable, which sets the size of the computational cell. In MPB, you can use this to specify a variety of affine lattice structures. In Meep, only rectangular Cartesian computational cells are supported, so the only property of lattice that you should normally use is its size.

size [vector3]
The size of the computational cell. Defaults to unit lengths.

If any dimension has the special size no-size, then the dimensionality of the problem is (essentially) reduced by one; strictly speaking, the dielectric function is taken to be uniform along that dimension.

Because Maxwell's equations are scale-invariant, you can use any units of distance you want to specify the cell size: nanometers, inches, parsecs, whatever. However, it is usually convenient to pick some characteristic lengthscale of your problem and set that length to 1. See also Units in Meep.


This class is used to specify the materials that geometric objects are made of. Currently, there are three subclasses, dielectric, perfect-metal, and material-function.

A uniform, isotropic, possibly nonlinear or dispersive, dielectric material; see also Dielectric materials in Meep. It has three properties:
epsilon [number]
The dielectric constant (must be positive). No default value. You can also use (index n) as a synonym for (epsilon (* n n)).
chi3 [number]
The Kerr susceptibility χ(3). (See e.g. Meep Tutorial/Third harmonic generation.)
polarizations [list of polarizability class]
List of dispersive polarizabilities (see below) added to the dielectric constant, in order to model material dispersion; defaults to none. (See e.g. Meep Tutorial/Material dispersion.)
A perfectly conducting metal; this class has no properties and you normally just use the predefined metal object, above. (To model imperfect conductors, use a dispersive dielectric material.)
This material type allows you to specify the material as an arbitrary function of position. It has one property:
material-func [function]
A function of one argument, the position vector3, that returns the material at that point. Note that the function you supply can return any material; wild and crazy users could even return another material-function object (which would then have its function invoked in turn).
Instead of material-func, you can use epsilon-func: for epsilon-func, you give it a function of position that returns the dielectric constant at that point.

Dispersive dielectric materials, above, are specified via a list of objects of type polarizability, which is another class with four properties:

Specifies a single dispersive polarizability of damped harmonic form (see Material dispersion), with the parameters:
omega [number]
The resonance frequency ωn.
gamma [number]
The resonance frequency γn.
delta-epsilon [number]
The amplitude \Delta\varepsilon_n.
sigma [number]
The scale factor σn. (Default value is 1.0.)
energy-saturation [number]
See Saturable gain in Meep. (Default value is 0.0.)


This class, and its descendants, are used to specify the solid geometric objects that form the dielectric structure being simulated. The base class is:

material [material-type class]
The material that the object is made of (usually some sort of dielectric). No default value (must be specified).
center [vector3]
Center point of the object. No default value.

One normally does not create objects of type geometric-object directly, however; instead, you use one of the following subclasses. Recall that subclasses inherit the properties of their superclass, so these subclasses automatically have the material and center properties (which must be specified, since they have no default values).

In a two-dimensional calculation, only the intersections of the objects with the x-y plane are considered.

A sphere. Properties:
radius [number]
Radius of the sphere. No default value.
A cylinder, with circular cross-section and finite height. Properties:
radius [number]
Radius of the cylinder's cross-section. No default value.
height [number]
Length of the cylinder along its axis. No default value.
axis [vector3]
Direction of the cylinder's axis; the length of this vector is ignored. Defaults to point parallel to the z axis.
A cone, or possibly a truncated cone. This is actually a subclass of cylinder, and inherits all of the same properties, with one additional property. The radius of the base of the cone is given by the radius property inherited from cylinder, while the radius of the tip is given by the new property:
radius2 [number]
Radius of the tip of the cone (i.e. the end of the cone pointed to by the axis vector). Defaults to zero (a "sharp" cone).
A parallelepiped (i.e., a brick, possibly with non-orthogonal axes). Properties:
size [vector3]
The lengths of the block edges along each of its three axes. Not really a 3-vector, but it has three components, each of which should be nonzero. No default value.
e1, e2, e3 [vector3]
The directions of the axes of the block; the lengths of these vectors are ignored. Must be linearly independent. They default to the three lattice directions.
An ellipsoid. This is actually a subclass of block, and inherits all the same properties, but defines an ellipsoid inscribed inside the block.

Here are some examples of geometric objects created using the above classes, assuming mat is some material we have defined:

; A cylinder of infinite radius and height 0.25 pointing along the x axis,
; centered at the origin:
(make cylinder (center 0 0 0) (material mat) 
               (radius infinity) (height 0.25) (axis 1 0 0))
; An ellipsoid with its long axis pointing along (1,1,1), centered on
; the origin (the other two axes are orthogonal and have equal
; semi-axis lengths):
(make ellipsoid (center 0 0 0) (material mat)
                (size 0.8 0.2 0.2)
               (e1 1 1 1)
               (e2 0 1 -1)
               (e3 -2 1 1))
; A unit cube of material m with a spherical air hole of radius 0.2 at
; its center, the whole thing centered at (1,2,3):
(set! geometry (list
               (make block (center 1 2 3) (material mat) (size 1 1 1))
               (make sphere (center 1 2 3) (material air) (radius 0.2))))


This class is used for the symmetries input variable, to specify symmetries (which must preserve both the structure and the sources) for Meep to exploit. Any number of symmetries can be exploited simultaneously, but there is no point in specifying redundant symmetries: the computational cell can be reduced by at most a factor of 4 in 2d and 8 in 3d. See also Exploiting symmetry in Meep.

A single symmetry to exploit. This is the base class of the specific symmetries below, so normally you don't create it directly. However, it has two properties which are shared by all symmetries:
direction [direction constant]
The direction of the symmetry (the normal to a mirror plane or the axis for a rotational symmetry). e.g. X, Y, ... (only Cartesian/grid directions are allowed). No default value.
phase [cnumber]
An additional phase to multiply the fields by when operating the symmetry on them; defaults to 1.0. e.g. a phase of -1 for a mirror plane corresponds to an odd mirror. (Technically, you are essentially specifying the representation of the symmetry group that your fields and sources transform under.)

The specific symmetry sub-classes are:

A mirror symmetry plane. Here, the direction is the direction normal to the mirror plane.
A 180° (twofold) rotational symmetry (a.k.a. C2). Here, the direction is the axis of the rotation.
A 90° (fourfold) rotational symmetry (a.k.a. C4). Here, the direction is the axis of the rotation.


This class is used for specifying the PML absorbing boundary layers around the cell, if any, via the pml-layers input variable. See also Perfectly matched layers. pml-layers can be zero or more pml objects, with multiple objects allowing you to specify different PML layers on different boundaries.

A single PML layer specification, which sets up one or more PML layers around the boundaries according to the following four properties.
thickness [number]
The spatial thickness of the PML layer (which extends from the boundary towards the inside of the computational cell). The thinner it is, the more numerical reflections become a problem. No default value.
direction [direction constant]
Specify the direction of the boundaries to put the PML layers next to. e.g. if X, then specifies PML on the \pm x boundaries (depending on the value of side, below). Default is the special value ALL, which puts PML layers on the boundaries in all directions.
side [boundary-side constant]
Specify which side, Low or High of the boundary or boundaries to put PML on. e.g. if side is Low and direction is X, then a PML layer is added to the x boundary. Default is the special value ALL, which puts PML layers on both sides.
strength [number]
A strength (default is 1.0) to multiply the PML absorption coefficient by. A strength of 2.0 will square the theoretical reflection coefficient of the PML (making it smaller), but will also increase numerical reflections.


The source class is used to specify the current sources (via the sources input variable). Note that all sources in Meep are separable in time and space, i.e. of the form \mathbf{J}(\mathbf{x},t) = \mathbf{A}(\mathbf{x}) \cdot f(t) for some functions \mathbf{A} and f. (Non-separable sources can be simulated, however, by modifying the sources after each time step.) When real fields are being used (which is the default in many cases...see the force-complex-fields? input variable), only the real part of the current source is used by Meep.

The source class has the following properties:
src [src-time class]
Specify the time-dependence of the source (see below). No default.
component [component constant]
Specify the direction and type of the current component: e.g. Ex, Ey, etcetera for an electric current, and Hx, Hy, etcetera for a magnetic current. Note that currents pointing in an arbitrary direction are specified simply as multiple current sources with the appropriate amplitudes for each component. No default.
center [vector3]
The location of the center of the current source in the computational cell; no default.
size [vector3]
The size of the current distribution along each direction of the computational cell. The default is (0,0,0): a point source.
amplitude [cnumber]
An overall (complex) amplitude multiplying the the current source. Default is 1.0.
amp-func [function]
A Scheme function of a single argument, that takes a vector3 giving a position and returns a (complex) current amplitude for that point. The position argument is relative to the center of the current source, so that you can move your current around without changing your function. The default is '() (null), meaning that a constant amplitude of 1.0 is used. Note that your amplitude function (if any) is multiplied by the amplitude property, so both properties can be used simultaneously.

The src-time object, which specifies the time dependence of the source, can be one of the following three classes.

A continuous-wave source proportional to exp( − iωt), possibly with a smooth (exponential/tanh) turn-on/turn-off. It has the properties:
frequency [number]
The frequency ω, in units of 2πc (see Units in Meep). No default value. You can instead specify (wavelength x) or (period x), which are both a synonym for (frequency (/ 1 x)); i.e. 1/ω in these units is the vacuum wavelength or the temporal period.
start-time [number]
The starting time for the source; default is 0 (turn on at t = 0).
end-time [number]
The end time for the source; default is infinity (never turn off).
width [number]
Roughly, the temporal width of the smoothing (technically, the inverse of the exponential rate at which the current turns off and on). Default is 0 (no smoothing). You can instead specify (fwidth x), which is a synonym for (width (/ 1 x)) (i.e. the frequency width is proportional to the inverse of the temporal width).
cutoff [number]
How many widths the current decays for before we cut it off and set it to zero; the default is 3.0. A larger value of cutoff will reduce the amount of high-frequency components that are introduced by the start/stop of the source, but will of course lead to longer simulation times.
A Gaussian-pulse source roughly proportional to exp( − iωt − (tt0)2 / 2w2), It has the properties:
frequency [number]
The center frequency ω, in units of 2πc (see Units in Meep). No default value. You can instead specify (wavelength x) or (period x), which are both a synonym for (frequency (/ 1 x)); i.e. 1/ω in these units is the vacuum wavelength or the temporal period.
width [number]
The width w used in the Gaussian. No default value. You can instead specify (fwidth x), which is a synonym for (width (/ 1 x)) (i.e. the frequency width is proportional to the inverse of the temporal width).
start-time [number]
The starting time for the source; default is 0 (turn on at t = 0).
cutoff [number]
How many widths the current decays for before we cut it off and set it to zero—this applies for both turn-on and turn-off of the pulse. The default is 5.0. A larger value of cutoff will reduce the amount of high-frequency components that are introduced by the start/stop of the source, but will of course lead to longer simulation times.
A user-specified source function f(t). You can also specify start/end times (at which point your current is set to zero whether or not your function is actually zero). These are optional, but you must specify an end-time explicitly if you want functions like run-sources to work, since they need to know when your source turns off.
src-func [function]
The function f(t) specifying the time-dependence of the source. It should take one argument (the time in Meep units) and return a complex number.
start-time [number]
The starting time for the source; default is (- infinity) (turn on at t=-\infty). Note, however, that the simulation normally starts at t = 0 with zero fields as the initial condition, so there is implicitly a sharp turn-on at t = 0 whether you specify it or not.
end-time [number]
The end time for the source; default is infinity (never turn off).


A flux-region object is used with add-flux below to specify a region in which Meep should accumulate the appropriate Fourier-transformed fields in order to compute a flux spectrum.

A region (volume, plane, line, or point) in which to compute the integral of the Poynting vector of the Fourier-transformed fields. Its properties are:
center [vector3]
The center of the flux region (no default).
size [vector3]
The size of the flux region along each of the coordinate axes; default is (0,0,0) (a single point).
direction [direction constant]
The direction in which to compute the flux (e.g. X, Y, etcetera). The default is AUTOMATIC, in which the direction is determined by taking the normal direction if the flux region is a plane (or a line, in 2d). If the normal direction is ambiguous (e.g. for a point or volume), then you must specify the direction explicitly (not doing so will lead to an error).
weight [cnumber]
A weight factor to multiply the flux by when it is computed; default is 1.0.

Note that the flux is always computed in the positive coordinate direction, although this can effectively be flipped by using a weight of -1.0. (This is useful, for example, if you want to compute the outward flux through a box, so that the sides of the box add instead of subtract!)

Miscellaneous functions

Here, we describe a number of miscellaneous useful functions provided by Meep.

See also the reference section of the libctl manual, which describes a number of useful functions defined by libctl.

Geometry utilities

Some utility functions are provided to help you manipulate geometric objects:

(shift-geometric-object obj shift-vector)
Translate obj by the 3-vector shift-vector.
(geometric-object-duplicates shift-vector min-multiple max-multiple obj)
Return a list of duplicates of obj, shifted by various multiples of shift-vector (from min-multiple to max-multiple, inclusive, in steps of 1).
(geometric-objects-duplicates shift-vector min-multiple max-multiple obj-list)
Same as geometric-object-duplicates, except operates on a list of objects, obj-list. If A appears before B in the input list, then all the duplicates of A appear before all the duplicates of B in the output list.
(geometric-objects-lattice-duplicates obj-list [ ux uy uz ])
Duplicates the objects in obj-list by multiples of the Cartesian basis vectors, making all possible shifts of the "primitive cell" (see below) that fit inside the lattice cell. (This is useful for supercell calculations; see the [user-tutorial.html tutorial].) The primitive cell to duplicate is ux by uy by uz, in units of the Cartesian basis vectors. These three parameters are optional; any that you do not specify are assumed to be 1.
(point-in-object? point obj)
Returns whether or not the given 3-vector point is inside the geometric object obj.
(point-in-periodic-object? point obj)
As point-in-object?, but also checks translations of the given object by the lattice vectors.
(display-geometric-object-info indent-by obj)
Outputs some information about the given obj, indented by indent-by spaces.

Output file names

The output file names used by Meep, e.g. for HDF5 files, are automatically prefixed by the input variable filename-prefix. If filename-prefix is "" (the default), however, then Meep constructs a default prefix based on the current ctl file name with ".ctl" replaced by "-": e.g. tst.ctl implies a prefix of "tst-". You can get this prefix by running:

Return the current prefix string that is prepended, by default, to all file names.

If you don't want to use any prefix, then you should set filename-prefix to false.

In addition to the filename prefix, you can also specify that all the output files be written into a newly-created directory (if it does not yet exist). This is done by running:

(use-output-directory [dirname])
Put output in a subdirectory, which is created if necessary. If the optional argument dirname is specified, that that is the name of the directory. Otherwise, the directory name is the current ctl file name with ".ctl" replaced by "-out": e.g. tst.ctl implies a directory of "tst-out".


(volume (center ...) (size ...))
Many Meep functions require you to specify a volume in space, corresponding to the C++ type meep::geometric_volume. This function creates such a volume object, given the center and size properties (just like e.g. a block object). If the size is not specified, it defaults to (0,0,0), i.e. a single point.
Return the current simulation time (in simulation time units, not wall-clock time!). (e.g. during a run function.)

Field computations

Meep supports a large number of functions to perform computations on the fields. Currently, most of them are accessed via the lower-level C++/SWIG interface, but we are slowly adding simpler, higher-level versions of them here.

(get-field-point c pt)
Given a component or derived-component constant c and a vector3 pt, returns the value of that component at that point.
(get-epsilon-point pt)
Equivalent to (get-field-point Dielectric pt).
(meep-fields-flux-in-box fields dir box)
Given a meep::fields*, a direction constant, and a meep::volume*, returns the flux (the integral of \Re [\mathbf{E}^* \times \mathbf{H}]) in that volume. Most commonly, you will use the global fields object, a volume that is a plane or a line, and a direction perpendicular to it, e.g. (meep-fields-flux-in-box fields X (volume (center 0) (size 0 1 1))).
(meep-fields-electric-energy-in-box field box)
Given a meep::fields* (e.g. the global fields object) and a meep::volume*, returns the integral of the electric-field energy \mathbf{E}^* \cdot \mathbf{D}/2 in the given volume. (If the volume has zero size along a dimension, a lower-dimensional integral is used.)
(meep-fields-magnetic-energy-in-box field box)
Given a meep::fields* (e.g. the global fields object) and a meep::volume*, returns the integral of the magnetic-field energy |\mathbf{H}|^2/2in the given volume. (If the volume has zero size along a dimension, a lower-dimensional integral is used.)
(meep-fields-field-energy-in-box field box)
Given a meep::fields* (e.g. the global fields object) and a meep::volume*, returns the integral of the electric+magnetic-field energy \mathbf{E}^* \cdot \mathbf{D}/2 + |\mathbf{H}|^2/2in the given volume. (If the volume has zero size along a dimension, a lower-dimensional integral is used.)

Note that if you are at a fixed frequency and you use complex fields (Bloch-periodic boundary conditions or fields-complex?=true), then one half of the flux or energy integrals above corresponds to the time-average of the flux or energy for a simulation with real fields.

One powerful feature is that you can supply an arbitrary function f(\mathbf{x},c_1,c_2,\ldots) of position \mathbf{x} and various field components c_1,\ldots and ask Meep to integrate it over a given volume, find its maximum, or output it (via output-field-function, described later). This is done via the functions:

(integrate-field-function cs func [where] [fields-var])
Returns the integral of the complex-valued function func over the meep::geometric_volume specified by where (defaults to entire computational cell) for the meep::fields specified by fields-var (defaults to fields). func is a function of position (a vector3, its first argument) and zero or more field components specified by cs: a list of component constants. func can be real- or complex-valued.
(If any dimension of where is zero, that dimension is not integrated over. In this way you can specify one-, two-, or three-dimensional integrals.)
(max-abs-field-function cs func [where] [fields-var])
As integrate-field-function, but returns the maximum absolute value of func in the volume where instead of its integral.

(The integration is performed by summing over the grid points with a simple trapezoidal rule, and the maximum is similarly over the grid points.)

Reloading parameters

Once the fields/simulation have been initialized, you can change the values of various parameters by using the following functions:

Reset all of Meep's parameters, deleting the fields, structures, etcetera, from memory as if you had not run any computations.
Restart the fields at time zero, with zero fields. (Does not reset the Fourier transforms of the flux planes, which continue to be accumulated.)
(change-k-point! k)
Change the k-point (the Bloch periodicity).
(change-sources! new-sources)
Change the sources input variable to new-sources, and changes the sources used for the current simulation.

(More to come...)

Flux spectra

Given a bunch of flux-region objects (see above), you can tell Meep to accumulate the Fourier transforms of the fields in those regions in order to compute flux spectra. See also the transmission/reflection spectra introduction and the Meep tutorial. The most important function is:

(add-flux fcen df nfreq flux-regions...)
Add a bunch of flux-regions to the current simulation (initializing the fields if they have not yet been initialized), telling Meep to accumulate the appropriate field Fourier transforms for nfreq equally spaced frequencies covering the frequency range fcen-df/2 to fcen+df/2. Return a flux object, which you can pass to the functions below to get the flux spectrum, etcetera.

As described in the tutorial, you normally use add-flux via statements like:

(define transmission (add-flux ...))

to store the flux object in a variable. add-flux initializes the fields if necessary, just like calling run, so you should only call it after setting up your geometry, sources, pml-layers, etcetera. You can create as many flux objects as you want, e.g. to look at powers flowing in different regions or in different frequency ranges. Note, however, that Meep has to store (and update at every time step) a number of Fourier components equal to the number of grid points intersecting the flux region multiplied by the number of electric and magnetic field components required to get the Poynting vector multiplied by nfreq, so this can get quite expensive (in both memory and time) if you want a lot of frequency points over large regions of space.

Once you have called add-flux, the Fourier transforms of the fields are accumulated automatically during time-stepping by the run functions. At any time, you can ask for Meep to print out the current flux spectrum via:

(display-fluxes fluxes...)
Given a number of flux objects, this displays a comma-separated table of frequencies and flux spectra, prefixed by "flux1:" or similar (where the number is incremented after each run). All of the fluxes should be for the same fcen/df/nfreq. The first column are the frequencies, and subsequent columns are the flux spectra.

You might have to do something lower-level if you have multiple flux regions corresponding to different frequency ranges, or have other special needs. (display-fluxes f1 f2 f3) is actually equivalent to (display-csv "flux" (get-flux-freqs f1) (get-fluxes f1) (get-fluxes f2) (get-fluxes f3), where display-csv takes a bunch of lists of numbers and prints them as a comma-separated table, and we are calling two lower-level functions:

(get-flux-freqs flux)
Given a flux object, returns a list of the frequencies that it is computing the spectrum for.
(get-fluxes flux)
Given a flux object, returns a list of the current flux spectrum that it has accumulated.

As described in the Meep tutorial, for a reflection spectrum you often want to save the Fourier-transformed fields from a "normalization" run and then load them into another run to be subtracted. This can be done via:

(save-flux filename flux)
Save the Fourier-transformed fields corresponding to the given flux object in an HDF5 file of the given name (without the ".h5" suffix) (the current filename-prefix is prepended automatically).
(load-flux filename flux)
Load the Fourier-transformed fields into the given flux object (replacing any values currently there) from an HDF5 file of the given name (without the ".h5" suffix) (the current filename-prefix is prepended automatically). You must load from a file that was saved by save-flux in a simulation of the same dimensions (for both the computational cell and the flux regions) with the same number of processors.
(load-minus-flux filename flux)
As load-flux, but negates the Fourier-transformed fields after they are loaded. This means that they will be subtracted from any future field Fourier transforms that are accumulated.
(scale-flux-fields s flux)
Scale the Fourier-transformed fields in flux by the complex number s. e.g. load-minus-flux is equivalent to load-flux followed by scale-flux-fields with s=-1.

Run and step functions

The actual work in Meep is performed by run functions, which time-step the simulation for a given amount of time or until a given condition is satisfied.

The run functions, in turn, can be modified by use of step functions: these are called at every time step and can perform any arbitrary computation on the fields, do outputs and I/O, or even modify the simulation. The step functions can be transformed by many modifier functions, like at-beginning, during-sources, etcetera which cause them to only be called at certain times, etcetera, instead of at every time step.

Run functions

The following run functions are available. (You can also write your own, using the lower-level C++/SWIG functions, but these should suffice for most needs.)

(run-until cond?/time step-functions...)
Run the simulation until a certain time or condition, calling the given step functions (if any) at each timestep. The first argument is either a number, in which case it is an additional time (in Meep units) to run for, or it is a function (of no arguments) which returns true when the simulation should stop.
(run-sources step-functions...)
Run the simulation until all sources have turned off, calling the given step functions (if any) at each timestep. Note that this does not mean that the fields will be zero at the end: in general, some fields will still be bouncing around that were excited by the sources.
(run-sources+ cond?/time step-functions...)
As run-sources, but with an additional first argument: either a number, in which case it is an additional time (in Meep units) to run for after the sources are off, or it is a function (of no arguments). In the latter case, the simulation runs until the sources are off and (cond?) returns true.

In particular, a useful first argument to run-sources+ or run-until is often given by (as in the Meep tutorial):

(stop-when-fields-decayed dT c pt decay-by)
Return a cond? function, suitable for passing to run-until/run-sources+, that examines the component c (e.g. Ex, etc.) at the point pt (a vector3) and keeps running until its absolute value squared has decayed by at least decay-by from its maximum previous value. In particular, it keeps incrementing the run time by dT (in Meep units) and checks the maximum value over that time period—in this way, it won't be fooled just because the field happens to go through 0 at some instant.
Note that, if you make decay-by very small, you may need to increase the cutoff property of your source(s), to decrease the amplitude of the small high-frequency components that are excited when the source turns off. (High frequencies near the Nyquist frequency of the grid have slow group velocities and are absorbed poorly by PML.)

Finally, another two run functions, useful for computing ω(k) band diagrams, are

(run-k-point T k)
Given a vector3 k, runs a simulation for each k point (i.e. specifying Bloch-periodic boundary conditions) and extracts the eigen-frequencies, and returns a list of the (complex) frequencies. In particular, you should have specified one or more Gaussian sources. It will run the simulation until the sources are turned off plus an additional T time units. It will run harminv (see below) at the same point/component as the first Gaussian source and look for modes in the union of the frequency ranges for all sources.
(run-k-points T k-points)
Given a list k-points of k vectors, runs run-k-point for each one, and returns a list of lists of frequencies (one list of frequencies for each k). Also prints out a comma-delimited list of frequencies, prefixed by freqs:, and their imaginary parts, prefixed by freqs-im:. (See e.g. this band diagram tutorial.)

Predefined step functions

Several useful step functions are predefined for you by Meep.

Output functions

The most common step function is an output function, which outputs some field component to an HDF5 file. Normally, you will want to modify this by one of the at- functions, below, as outputting a field at every time step can get quite time- and storage-consuming.

Note that although the various field components are stored at different places in the Yee lattice, when they are outputted they are all linearly interpolated to the same grid: to the points at the centers of the Yee cells, i.e. (i+0.5,j+0.5,k+0.5)\cdot\Delta in 3d.

The predefined output functions are:

Output the dielectric function ε.
Output the magnetic-field energy density |\mathbf{H}|^2 / 2
Output the electric-field energy density \mathbf{E}^* \cdot \mathbf{D} / 2
Output the total electric and magnetic energy density.
output-Xfield-x, output-Xfield-y, output-Xfield-z, output-Xfield-r, output-Xfield-p
Output the x, y, z, r, or φ component respectively, of the field X, where X is either h, e, d, or s for the magnetic, electric, displacement, or Poynting field, respectively. If the field is complex, outputs two datasets, e.g. ex.r and ex.i, within the same HDF5 file for the real and imaginary parts, respectively.
Outputs all the components of the field X, where X is either h, e, d, or s as above, to an HDF5 file. (That is, the different components are stored as different datasets within the same file.)
(output-png component h5topng-options)
Output the given field component (e.g. Ex, etc.) as a PNG image, by first outputting the HDF5 file, then converting to PNG via h5topng, then deleting the HDF5 file. The second argument is a string giving options to pass to h5topng (e.g. "-Zc bluered"). See also the Meep tutorial.

More generally, it is possible to output an arbitrary function of position and zero or more field components, similar to the integrate-field-function described above. This is done by:

(output-field-function name cs func)
Output the field function func to an HDF5 file in the datasets named name.r and name.i (for the real and imaginary parts). Similar to integrate-field-function, func is a function of position (a vector3) and the field components corresponding to cs: a list of component constants.
(output-real-field-function name cs func)
As output-field-function, but only outputs the real part of func to the dataset given by the string name.


The following step function collects field data from a given point and runs Harminv on that data to extract the frequencies, decay rates, and other information.

(harminv c pt fcen df [maxbands])
Returns a step function that collects data from the field component c (e.g. Ex, etc.) at the given point pt (a vector3). Then, at the end of the run, it uses Harminv to look for modes in the given frequency range (center fcen and width df), printing the results to standard output (prefixed by harminv:) as comma-delimited text, and also storing them to the variable harminv-results. The optional argument maxbands is the maximum number of modes to search for; defaults to 100.

Important: normally, you should only use harminv to analyze data after the sources are off. Wrapping it in (after-sources (harminv ...)) is sufficient.

In particular, Harminv takes the time series f(t) corresponding to the given field component as a function of time and decomposes it (within the specified bandwidth) as:

f(t) = \sum_n a_n e^{-i\omega_n t}

The results are stored in the list harminv-results, which is a list of tuples holding the frequency, amplitude, and error of the modes. Given one of these tuples, you can extract its various components with one of the accessor functions:

(harminv-freq result)
Return the complex frequency ω (in the usual Meep 2πc units).
(harminv-freq-re result)
Return the real part of the frequency ω.
(harminv-freq-im result)
Return the imaginary part of the frequency ω.
(harminv-Q result)
Return dimensionless lifetime, or "quality factor", Q, defined as -\mathrm{Re}\,\omega / 2 \mathrm{Im}\,\omega.
(harminv-amp result)
Return the complex amplitude a.
(harminv-err result)
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!

For example, (map harminv-freq-re harminv-results) gives you a list of the real parts of the frequencies, using the Scheme built-in map.

Step-function modifiers

Rather than writing a brand-new step function every time we want to do something a bit different, the following "modifier" functions take a bunch of step functions and produce new step functions with modified behavior. (See also the Meep tutorial for examples.)

Controlling when a step function executes

(when-true cond? step-functions...)
Given zero or more step functions and a condition function cond? ( a function of no arguments), evaluate the step functions whenever (cond?) returns true.
(when-false cond? step-functions...)
Given zero or more step functions and a condition function cond? ( a function of no arguments), evaluate the step functions whenever (cond?) returns false.
(at-every dT step-functions...)
Given zero or more step functions, evaluates them at every time interval of dT units (rounded up to the next time step).
(after-time T step-functions...)
Given zero or more step functions, evaluates them only for times after a T time units have elapsed from the start of the run.
(before-time T step-functions...)
Given zero or more step functions, evaluates them only for times before a T time units have elapsed from the start of the run.
(at-time T step-functions...)
Given zero or more step functions, evaluates them only once, after a T time units have elapsed from the start of the run.
(after-sources step-functions...)
Given zero or more step functions, evaluates them only for times after all of the sources have turned off.
(after-sources+ T step-functions...)
Given zero or more step functions, evaluates them only for times after all of the sources have turned off, plus an additional T time units have elapsed.
(during-sources step-functions...)
Given zero or more step functions, evaluates them only for times before all of the sources have turned off.
(at-beginning step-functions...)
Given zero or more step functions, evaluates them only once, at the beginning of the run.
(at-end step-functions...)
Given zero or more step functions, evaluates them only once, at the end of the run.

Modifying HDF5 output

(in-volume v step-functions...)
Given zero or more step functions, modifies any output functions among them to only output a subset (or a superset) of the computational cell, corresponding to the meep::geometric_volume* v (created by the volume function).
(in-point pt step-functions...)
Given zero or more step functions, modifies any output functions among them to only output a single point of data, at pt (a vector3).
(to-appended filename step-functions...)
Given zero or more step functions, modifies any output functions among them to append their data to datasets in a single newly-created file named filename (plus an .h5 suffix and the current filename prefix). They append by adding an extra dimension to their datasets, corresponding to time.

Writing your own step functions

A step function can take two forms. The simplest is just a function of no arguments, which is called at every time step (unless modified by one of the modifier functions above). e.g.

(define (my-step) (print "Hello world!\n"))

If one then does (run-until 100 my-step), Meep will run for 100 time units and print "Hello world!" at every time step.

This suffices for most purposes. However, sometimes you need a step function that opens a file, or accumulates some computation, and you need to clean up (e.g. close the file or print the results) at the end of the run. For this case, you can write a step function of one argument: that argument will either be 'step when it is called during time-stepping, or 'finish when it is called at the end of the run.

Low-level functions

By default, Meep reads input functions like sources and geometry and creates global variables structure and fields to store the corresponding C++ objects. Given these, you can then call essentially any function in the C++ interface, because all of the C++ functions are automatically made accessible to Scheme by a wrapper-generator program called SWIG

Initializing the structure and fields

The structure and fields variables are automatically initialized when any of the run functions is called, or by various other functions such as add-flux. To initialize them separately, you can call (init-fields) manually, or (init-structure k-point) to just initialize the structure.

If you want to time step more than one field simultaneously, the easiest way is probably to do something like:

(define my-fields fields)
(set! fields '())

and then change the geometry etc. and re-run (init-fields). Then you'll have two field objects in memory.

SWIG wrappers

If you look at a function in the C++ interface, then there are a few simple rules to infer the name of the corresponding Scheme function.

  • First, all Meep functions (in the meep:: namespace) are prefixed with meep- in the Scheme interface.
  • Second, any method of a class is prefixed with the name of the class and a hyphen. For example, meep::fields::step, which is the function that performs a time-step, is exposed to Scheme as meep-fields-step. Moreover, you pass the object as the first argument in the Scheme wrapper. e.g. f.step() becomes (meep-fields-step f).
  • To call the C++ constructor for a type, you use new-type. e.g. (new-meep-fields ...arguments...) returns a new meep::fields object. Conversely, to call the destructor and deallocate an object, you use delete-type; most of the time, this is not necessary because objects are automatically garbage-collected.

Some argument type conversion is performed automatically, e.g. types like complex numbers are converted to complex<double>, etcetera. vector3 vectors are converted to meep::vec, but to do this we need to know the dimensionality of the problem in C++. The problem dimensions are automatically initialized by init-structure, but if you want to pass vector arguments to C++ before that time you should call (require-dimensions!), which infers the dimensions from the geometry-lattice, k-point, and dimensions variables.

Personal tools