Meep Reference
From AbInitio
(diff) ←Older revision  Current revision  Newer revision→ (diff)
Meep 
Download 
Release notes 
FAQ 
Meep manual 
Introduction 
Installation 
Tutorial 
Reference 
C++ Tutorial 
C++ Reference 
Acknowledgements 
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.
Contents 
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 geometricobject
, are described in a later subsection. The basic datatypes, like integer
, boolean
, cnumber
, and vector3
, are defined by libctl.)

geometry
[list ofgeometricobject
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 ofsource
class]  Specifies the current sources to be present in the simulation; defaults to none.

symmetries
[list ofsymmetry
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.

pmllayers
[list ofpml
class]  Specifies the absorbing PML boundary layers to use; defaults to none.

geometrylattice
[lattice
class]  Specifies the size of the unit cell (which is centered on the origin of the coordinate system). Any sizes of
nosize
imply (effectively) a reduceddimensionality calculation (but a 2d xy calculation is especially optimized); seedimensions
below. Defaults to a cubic cell of unit size. 
defaultmaterial
[materialtype
class]  Holds the default material that is used for points not in any object of the geometry list. Defaults to
air
(ε of 1). See alsoepsiloninputfile
below. 
epsiloninputfile
[string
]  If this string is not
""
(the default), then it should be the name of an HDF5 file whose first/only dataset defines a scalar dielectric function (over some discrete grid); alternatively, the dataset name can be specified explicitly if the string is in the form "filename:dataset". This dielectric function is then used in place of the ε property ofdefaultmaterial
(i.e. where there are nogeometry
objects). The grid of the epsilon file dataset need not match Meep's computational grid; it is scaled and/or linearly interpolated as needed to map the file onto the computational cell (which warps the structure if the proportions of the grids do not match, however). Note: the file contents only override the ε property of thedefaultmaterial
, whereas other properties (μ, susceptibilities, nonlinearities, etc.) ofdefaultmaterial
are still used. 
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
geometrylattice
size in the z direction isnosize
. Ifdimensions
is the special value ofCYLINDRICAL
, then cylindrical coordinates are used and the x and z dimensions are interpreted as r and z, respectively. Ifdimensions
is 1, then the cell must be along the z direction and only E_{x} and H_{y} field components are permitted. Ifdimensions
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 e^{imφ} (default ism=0
). If the simulation cell includes the origin r = 0, thenm
must be an integer. 
accuratefieldsnearcylorigin?
[boolean
]  For
CYLINDRICAL
simulations with m > 1, compute more accurate fields near the origin r = 0 at the expense of requiring a smaller Courant factor. Empirically, when this option is set totrue
, a Courant factor of roughly min[0.5,1 / (  m  + 0.5)] (or smaller) seems to be needed. The default isfalse
, in which case the D_{r}, D_{z}, and B_{r} fields within m pixels of the origin are forced to zero, which usually ensures stability with the default Courant factor of 0.5, at the expense of slowing convergence of the fields near r = 0. 
resolution
[number
]  Specifies the computational grid resolution, in pixels per distance unit. Defaults to
10
. 
kpoint
[false
orvector3
]  If
false
(the default), then the boundaries are perfect metallic (zero electric field). If a vector, then the boundaries are Blochperiodic: the fields at one side are times the fields at the other side, separated by the lattice vector . Thekpoint
vector is specified in Cartesian coordinates, in units of 2π/distance. (This is different from MPB, equivalent to taking MPB'skpoints
through the functionreciprocal>cartesian
.) 
ensureperiodicity
[boolean
]  If
true
(the default), and if the boundary conditions are periodic (kpoint
is notfalse
), then the geometric objects are automatically repeated periodically according to the lattice vectors (the size of the computational cell). 
epsaveraging?
[boolean
]  If
true
(the default), then subpixel averaging is used when initializing the dielectric function (see the Farjadpour et al. reference in Citing Meep). The input variablessubpixelmaxeval
(default100000
) andsubpixeltol
(default1.0e4
) specify the maximum number of function evaluations and the integration tolerance for subpixel averaging. Increasing/decreasing these, respectively, will cause a more accurate (but slower) computation of the average ε (with diminishing returns for the actual FDTD error). 
forcecomplexfields?
[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
kpoint
is nonzero or ifm
is nonzero. However, by settingforcecomplexfields?
totrue
, Meep will always use complex fields. See also: Complex fields in Meep. 
filenameprefix
[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 , where n_{min} is the minimum refractive index (usually 1), and in practice S should be slightly smaller. 
outputvolume
[meep::geometric_volume*
]  Specifies the default region of space that is output by the HDF5 output functions (below); see also the
(volume ...)
function to createmeep::geometric_volume*
objects. The default is'()
(null), which means that the whole computational cell is output. Normally, you should use the(involume ...)
function to modify the output volume instead of settingoutputvolume
directly. 
outputsingleprecision?
[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 isfalse
) you can instead output in single precision which saves a factor of two in space. 
progressinterval
[number
]  Time interval (seconds) after which Meep prints a progress message; default is 4 seconds.

extramaterials
[list ofmaterialtype
class]  By default, Meep turns off support for material dispersion, nonlinearities, and similar properties if none of the objects in the
geometry
have materials with these properties—since they are not needed, it is faster to omit their calculation. This doesn't work if you use amaterialfunction
, though (materials via a userspecified function of position, instead of just geometric objects). If only your material function returns a nonlinear material, for example, Meep won't notice this unless you tell it explicitly viaextramaterials
.extramaterials
is a list of materials that Meep should look for in the computational cell, in addition to any materials that are specified by geometric objects. You should list any materials (other than scalar dielectrics) that are returned bymaterialfunction
s here.
The following 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
(initstructure)
which is called automatically by(initfields)
which is called automatically by any of the(run)
functions. 
fields
[meep::fields*
]  Pointer to the current fields being simulated; initialized by
(initfields)
which is called automatically by any of the(run)
functions. 
numchunks
[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
[materialtype
class]  Two aliases for a predefined material type with a dielectric constant of 1.

perfectelectricconductor
ormetal
[materialtype
class]  A predefined material type corresponding to a perfect electric conductor (at the boundary of which the parallel electric field is zero), technically .

perfectmagneticconductor
[materialtype
class]  A predefined material type corresponding to a perfect magnetic conductor (at the boundary of which the parallel magnetic field is zero), technically .

nothing
[materialtype
class]  A material that, effectively, punches a hole through other objects to the background (
defaultmaterial
). 
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. 
boundaryside
constants  Specify particular boundary in some direction (e.g. + x or − x). One of
Low
orHigh
. 
component
constants  Specify a particular field (or other) component. One of
Ex
,Ey
,Ez
,Er
,Ep
,Hx
,Hy
,Hz
,Hy
,Hp
,Hz
,Bx
,By
,Bz
,By
,Bp
,Bz
,Dx
,Dy
,Dz
,Dr
,Dp
,Dielectric
,Permeability
, for: E_{x}, E_{y}, E_{z}, E_{r}, E_{φ}, H_{x}, H_{y}, H_{z}, H_{r}, H_{φ}, B_{x}, B_{y}, B_{z}, B_{r}, B_{φ}, D_{x}, D_{y}, D_{z}, D_{r}, D_{φ}, , μ, respectively. 
derivedcomponent
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
,DEnergyDensity
,HEnergyDensity
for: S_{x}, S_{y}, S_{z}, S_{r}, S_{φ} (components of the Poynting vector ), , , , respectively.
Classes
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.
lattice
The lattice class is normally used only for the geometrylattice
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
.

lattice
 Properties:

size
[vector3
]  The size of the computational cell. Defaults to unit lengths.
If any dimension has the special size nosize
, 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 scaleinvariant, 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.
materialtype
This class is used to specify the materials that geometric objects are made of. Currently, there are three subclasses, dielectric
, perfectmetal
, and materialfunction
.

medium
 An electromagnetic medium (possibly nonlinear and/or dispersive); see also Materials in Meep. For backwards compatibility, a synonym for
medium
isdielectric
. It has several properties:
epsilon
[number
]  The (frequencyindependent) isotropic relative permittivity, or dielectric constant. Default value = 1. You can also use
(index n)
as a synonym for(epsilon (* n n))
; note that this is not really the refractive index if you also specify μ, since the true index is .  Using
(epsilon ε)
is actually a synonym for(epsilondiag ε ε ε)
. 
epsilondiag
andepsilonoffdiag
[vector3
]  These properties allow you to specify ε as an arbitrary realsymmetric tensor, by giving the diagonal and offdiagonal parts. Specifying
(epsilondiag a b c)
and/or(epsilonoffdiag u v w)
corresponds to a relative permittivity tensor  The default is the identity matrix (a = b = c = 1 and u = v = w = 0).

mu
[number
]  The (frequencyindependent) isotropic relative permeability μ. Default value = 1. Using
(mu μ)
is actually a synonym for(mudiag μ μ μ)
. 
mudiag
andmuoffdiag
[vector3
]  These properties allow you to specify μ as an arbitrary realsymmetric tensor, by giving the diagonal and offdiagonal parts exactly as for ε above, again defaulting to the identity matrix.

Dconductivity
[number
]  The (frequencyindependent) electric conductivity σ_{D}. Default value = 0. You can also specify an diagonal anisotropic conductivity tensor by using the property
Dconductivitydiag
[vector3], which takes three numbers or a vector3 to give the σ_{D} tensor diagonal. See also Conductivity in Meep. 
Bconductivity
[number
]  The (frequencyindependent) magnetic conductivity σ_{B}. Default value = 0. You can also specify an diagonal anisotropic conductivity tensor by using the property
Bconductivitydiag
[vector3], which takes three numbers or a vector3 to give the σ_{B} tensor diagonal. See also Conductivity in Meep. 
chi2
[number
]  The nonlinear (Pockels) susceptibility χ^{(2)}. Default value = 0. See also Nonlinearity in Meep.

chi3
[number
]  The nonlinear (Kerr) susceptibility χ^{(3)}. Default value = 0. (See e.g. Meep Tutorial/Third harmonic generation.) See also Nonlinearity in Meep.

Esusceptibilities
[list ofsusceptibility
class]  List of dispersive susceptibilities (see below) added to the dielectric constant ε, in order to model material dispersion; defaults to none. (See e.g. Meep Tutorial/Material dispersion.) See also Material dispersion in Meep. For backwards compatibility, synonyms of
Esusceptibilities
areEpolarizations
andpolarizations
. 
Hsusceptibilities
[list ofsusceptibility
class]  List of dispersive susceptibilities (see below) added to the permeability μ, in order to model material dispersion; defaults to none. (See e.g. Meep Tutorial/Material dispersion.) See also Material dispersion in Meep. For backwards compatibility, a synonym of
Hsusceptibilities
isHpolarizations
.


perfectmetal
 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.) See also the predefined variablesmetal
,perfectelectricconductor
, andperfectmagneticconductor
above. 
materialfunction
 This material type allows you to specify the material as an arbitrary function of position. It has one property:

materialfunc
[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 anothermaterialfunction
object (which would then have its function invoked in turn). Instead of
materialfunc
, you can useepsilonfunc
: forepsilonfunc
, you give it a function of position that returns the dielectric constant at that point.
 Instead of
 Important: If your material function returns nonlinear, dispersive (Lorentzian or conducting), or magnetic materials, you should also include a list of these materials in the
extramaterials
input variable (above) to let Meep know that it needs to support these material types in your simulation. (For dispersive materials, you need to include a material with the same γ_{n} and ω_{n} values, so you can only have a finite number of these, whereas σ_{n} can vary continuously if you want and a matching σ_{n} need not be specified inextramaterials
. For nonlinear or conductivity materials, yourextramaterials
list need not match the actual values of σ or χ returned by your material function, which can vary continuously if you want.)

Complex ε and μ: you cannot specify a frequencyindependent complex ε or μ in Meep (where the imaginary part is a frequencyindependent loss), but there is an alternative. That is because there are only two important physical situations. First, if you only care about the loss in a narrow bandwidth around some frequency, you can set the loss at that frequency via the conductivity (see Conductivity in Meep). Second, if you care about a broad bandwidth, then all physical materials have a frequencydependent imaginary (and real) ε (and/or μ), and you need to specify that frequency dependence by fitting to Lorentzian (and/or Drude) resonances via the lorentziansusceptibility
(or drudesusceptibility
) classes below.
Dispersive dielectric (and magnetic) materials, above, are specified via a list of objects that are subclasses of type susceptibility
.

susceptibility
 Parent class for various dispersive susceptibility terms, parameterized by an anisotropic amplitude σ (see Material dispersion in Meep):

sigma
[number
]  The scale factor σ. You can also specify an anisotropic σ tensor by using the property
sigmadiag
[vector3], which takes three numbers or a vector3 to give the σ_{n} tensor diagonal, andsigmaoffdiag
[vector3] which specifies the offdiagonal elements (defaults to 0). That is,(sigmadiag a b c)
and(sigmaoffdiag u v w)
corresponds to a σ tensor


lorentziansusceptibility
 Specifies a single dispersive susceptibility of Lorentzian (damped harmonic oscillator) form (see Material dispersion in Meep), with the parameters (in addition to σ):

frequency
[number
]  The resonance frequency f_{n} = ω_{n} / 2π.

gamma
[number
]  The resonance loss rate γ_{n} / 2π.


drudesusceptibility
 Specifies a single dispersive susceptibility of Drude form (see Material dispersion in Meep), with the parameters (in addition to σ):

frequency
[number
]  The frequency scale factor f_{n} = ω_{n} / 2π which multiplies σ (not a resonance frequency).

gamma
[number
]  The loss rate γ_{n} / 2π.

Meep also supports a somewhat unusual polarizable medium, a Lorentzian susceptibility with a random noise term added into the dampedoscillator equation at each point. This can be used to directly model thermal radiation in both the far field (as described in Luo et al., 2004) and the near field (as described in Rodriguez et al., 2011). (Note, however that it is more efficient to compute farfield thermal radiation using Kirchhoff's law of radiation, which states that emissivity equals absorptivity. Nearfield thermal radiation can usually be computed more efficiently using frequency domain techniques, e.g. via our SCUFFEM code, as described by Rodriguez et al., 2013.)

noisylorentziansusceptibility
ornoisydrudesusceptibility
 Specifies a single dispersive susceptibility of Lorentzian (damped harmonic oscillator) or Drude form (see Material dispersion in Meep), with the same σ,
frequency
, andgamma
) parameters, but with an additional Gaussian random noise term (uncorrelated in space and time, zero mean) added to the P dampedoscillator equation.
noiseamp
[number
]  The noise has rootmean square amplitude σ⋅
noiseamp
.

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

geometricobject
 Properties:

material
[materialtype
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 geometricobject
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 twodimensional calculation, only the intersections of the objects with the xy plane are considered.

sphere
 A sphere. Properties:

radius
[number
]  Radius of the sphere. No default value.


cylinder
 A cylinder, with circular crosssection and finite height. Properties:

radius
[number
]  Radius of the cylinder's crosssection. 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.


cone
 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 theradius
property inherited fromcylinder
, while the radius of the tip is given by the new property,radius2
. (Thecenter
of a cone is halfway between the two circular ends.)
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).


block
 A parallelepiped (i.e., a brick, possibly with nonorthogonal axes). Properties:

size
[vector3
]  The lengths of the block edges along each of its three axes. Not really a 3vector, 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.


ellipsoid
 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 ; semiaxis 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))))
symmetry
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.

symmetry
 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 of1
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 subclasses are:

mirrorsym
 A mirror symmetry plane. Here, the
direction
is the direction normal to the mirror plane. 
rotate2sym
 A 180° (twofold) rotational symmetry (a.k.a. C_{2}). Here, the
direction
is the axis of the rotation. 
rotate4sym
 A 90° (fourfold) rotational symmetry (a.k.a. C_{4}). Here, the
direction
is the axis of the rotation.
pml
This class is used for specifying the PML absorbing boundary layers around the cell, if any, via the pmllayers
input variable. See also Perfectly matched layers. pmllayers
can be zero or more pml
objects, with multiple objects allowing you to specify different PML layers on different boundaries.

pml
 A single PML layer specification, which sets up one or more PML layers around the boundaries according to the following 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 boundaries (depending on the value ofside
, below). Default is the special valueALL
, which puts PML layers on the boundaries in all directions. 
side
[boundaryside
constant]  Specify which side,
Low
orHigh
of the boundary or boundaries to put PML on. e.g. if side isLow
and direction isX
, then a PML layer is added to the − x boundary. Default is the special valueALL
, which puts PML layers on both sides. 
strength
[number
]  A strength (default is
1.0
) to multiply the PML absorption coefficient by. A strength of2.0
will square the theoretical asymptotic reflection coefficient of the PML (making it smaller), but will also increase numerical reflections. Alternatively, you can changeRasymptotic
, below. 
Rasymptotic
[number
]  The asymptotic reflection in the limit of infinite resolution or infinite PML thickness, for refections from air (an upper bound for other media with index > 1). (For a finite resolution or thickness, the reflection will be much larger, due to the discretization of Maxwell's equation.) The default value is 10^{−15}, which should be suffice for most purposes. (You want to set this to be small enough so that waves propagating within the PML are attenuated sufficiently, but making
Rasymptotic
too small will increase the numerical reflection due to discretization.) 
pmlprofile
[function
]  By default, Meep turns on the PML conductivity quadratically within the PML layer—one doesn't want to turn it on suddenly, because that exacerbates reflections due to the discretization. More generally, with
pmlprofile
one can specify an arbitrary PML "profile" function f(u) that determines the shape of the PML absorption profile up to an overall constant factor. u goes from 0 to 1 at the start and end of the PML, and the default is f(u)=u^{2}. In some cases where a very thick PML is required, such as in a periodic medium (where there is technically no such thing as a true PML, only a pseudoPML), it can be advantageous to turn on the PML absorption more smoothly (see Oskooi et al., 2008). For example, one can use a cubic profile f(u)=u^{3} by specifying(pmlprofile (lambda (u) (* u u u)))
.

absorber
Instead of a pml
layer, there is an alternative class called absorber
which is a dropin replacement for pml
. For example, you can do (set! pmllayers (list (make absorber (thickness 2))))
instead of (set! pmllayers (list (make pml (thickness 2))))
. All the parameters are the same as for pml
, above. You can have a mix of pml
on some boundaries and absorber
on others.
The absorber
class does not implement a perfectly matched layer (PML), however (except in 1d). Instead, it is simply a scalar electric and magnetic conductivity that turns on gradually within the layer according to the pmlprofile
(defaulting to quadratic). Such a scalar conductivity gradient is only reflectionless in the limit as the layer becomes sufficiently thick.
The main reason to use absorber
is if you have a case in which PML fails:
 No true PML exists for periodic media, and a scalar absorber is cheaper and is generally just as good. See Oskooi et al. (2008).
 PML can lead to divergent fields for certain waveguides with "backwardwave" modes; this can easily happen in metallic with surface plasmons, and a scalar absorber is your only choice. See Loh et al. (2009).
 PML can fail if you have a waveguide hitting the edge of your computational cell at an angle. See Oskooi et. al. (2011).
source
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 for some functions and f. (Nonseparable 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 forcecomplexfields?
input variable), only the real part of the current source is used by Meep.
Important note: These are current sources (J terms in Maxwell's equations), even though they are labelled by electric/magnetic field components. They do not specify a particular electric/magnetic field (which would be what is called a "hard" source in the FDTD literature). There is no fixed relationship between the current source and the resulting field amplitudes; it depends on the surrounding geometry, as described in the Meep FAQ and in our book chapter online.

source
 The source class has the following properties:

src
[srctime
class]  Specify the timedependence 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 electriccharge current, andHx
,Hy
, etcetera for a magneticcharge 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 pointdipole source.
amplitude
[cnumber
] An overall (complex) amplitude multiplying the the current source. Default is
1.0
. ampfunc
[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 theamplitude
property, so both properties can be used simultaneously.

As described in section 4.2 of our book chapter online, it is also possible to supply a source that is designed to couple exclusively into a single waveguide mode (or other mode of some cross section or periodic region) at a single frequency, and which couples primarily into that mode as long as the bandwidth is not too broad. This is possible if you have MPB (version 1.5 or later) installed: Meep will call MPB to compute the field profile of the desired mode, and uses the field profile to produce an equivalent current source. (Note: this feature does not work in cylindrical coordinates.) To do this, instead of a source
you should use an eigenmodesource
:

eigenmodesource
 This is a subclass of
source
and has all of the properties ofsource
above. However, you normally do not specify acomponent
. Instead ofcomponent
, the current source components and amplitude profile are computed by calling MPB to compute the modes of the dielectric profile in the region given by thesize
andcenter
of the source, with the modes computed as if the source region were repeated periodically in all directions. (If anamplitude
and/orampfunc
are supplied, they are multiplied by this current profile.) The desired eigenmode and other features are specified by the following properties:
eigband
[integer
]  The index n (1,2,3,...) of the desired band ω_{n}(k) to compute in MPB (1 denotes the lowestfrequency band at a given k point, and so on).

direction
[X
,Y
, orZ
; defaultAUTOMATIC
],eigmatchfreq?
[boolean
; defaulttrue
],eigkpoint
[vector3
]  By default (if
eigmatchfreq?
istrue
), Meep tries to find a mode with the same frequency ω_{n}(k) as thesrc
property (above), by scanning k vectors in the givendirection
(using MPB's findk functionality; alternatively, ifeigkpoint
is supplied, it is used as an initial guess and direction for k ). By default,direction
is the direction normal to the source region, assumingsize
is d–1 dimensional in a ddimensional simulation (e.g. a plane in 3d). Alternatively (ifeigmatchfreq?
isfalse
), you can specify a particular k vector of the desired mode witheigkpoint
(in Meep units of 2π/a). 
eigparity
[NOPARITY
(= default),EVENZ
(=TE
),ODDZ
(=TM
),EVENY
,ODDY
]  The parity (= polarization in 2d) of the mode to calculate, assuming the structure has z and/or y mirror symmetry in the source region. If the structure has both y and z mirror symmetry, you can combine more than one of these, e.g.
EVENZ + ODDY
. The default isNOPARITY
, in which case MPB computes all of the bands (which will still be even or odd if the structure has mirror symmetry, of course). This is especially useful in 2d simulations to restrict yourself to a desired polarization. 
eigresolution
[integer
, defaults to same as Meep resolution]  The spatial resolution to use in MPB for the eigenmode calculations. This defaults to the same resolution as Meep, but you can use a higher resolution (in which case the structure is linearly interpolated from the Meep pixels).

eigtolerance
[number
, defaults to 10^{–7}]  The tolerance to use in the MPB eigensolver. (MPB terminates when the eigenvalues stop changing to less than this fractional tolerance.)

component
[as above, but defaults toALLCOMPONENTS
]  Once the MPB modes are computed, equivalent electric and magnetic sources are created within Meep. By default, these sources include magnetic and electric currents in all transverse directions within the source region, corresponding to the mode fields as described in our book chapter. If you specify a
component
property, however, you can include only one component of these currents if you wish. Most people won't need this feature. 
eiglatticesize
[vector3
],eiglatticecenter
[vector3
]  Normally, the MPB computational unit cell is the same as the source volume (given by the
size
andcenter
parameters). However, occasionally you want the unit cell to be larger than the source volume. For example, to create an eigenmode source in a periodic medium (photonic crystal), you need to pass MPB the entire unit cell of the periodic medium, but once the mode is computed then the actual current sources need only lie on a cross section of that medium. To accomplish this, you can specify the optionaleiglatticesize
andeiglatticecenter
, which define a volume (which must enclosesize
andcenter
) that is used for the unit cell in MPB (with the dielectric function ε taken from the corresponding region in the Meep simulation).

Note that MPB only supports dispersionless nonmagnetic materials (but it does support anisotropic ε). Any nonlinearities, magnetic responses µ, conductivities σ, or dispersive polarizations in your materials will be ignored when computing the eigenmode source. PML will also be ignored.
The srctime
object, which specifies the time dependence of the source, can be one of the following three classes.

continuoussrc
 A continuouswave source proportional to exp( − iωt), possibly with a smooth (exponential/tanh) turnon/turnoff. It has the properties:
frequency
[number
] The frequency f in units of c/distance (or ω in units of 2πc/distance) (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. starttime
[number
] The starting time for the source; default is
0
(turn on at t = 0). endtime
[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
width
s the current decays for before we cut it off and set it to zero; the default is3.0
. A larger value ofcutoff
will reduce the amount of highfrequency components that are introduced by the start/stop of the source, but will of course lead to longer simulation times.

gaussiansrc
 A Gaussianpulse source roughly proportional to exp( − iωt − (t − t_{0})^{2} / 2w^{2}). Technically, the "Gaussian" sources in Meep are the (discretetime) derivative of a Gaussian, i.e. they are , but the difference between this and a true Gaussian is usually irrelevant. It has the properties:
frequency
[number
] The center frequency f in units of c/distance (or ω in units of 2πc/distance) (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). starttime
[number
] The starting time for the source; default is
0
(turn on at t = 0). (Not the time of the peak! See below.) cutoff
[number
] How many
width
s the current decays for before we cut it off and set it to zero—this applies for both turnon and turnoff of the pulse. The default is5.0
. A larger value ofcutoff
will reduce the amount of highfrequency components that are introduced by the start/stop of the source, but will of course lead to longer simulation times. The peak of the gaussian is reached at the time t_{0}=starttime
+cutoff
*width
.

customsrc
 A userspecified 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
endtime
explicitly if you want functions likerunsources
to work, since they need to know when your source turns off.
srcfunc
[function
]  The function f(t) specifying the timedependence of the source. It should take one argument (the time in Meep units) and return a complex number.
starttime
[number
] The starting time for the source; default is
( infinity)
(turn on at ). Note, however, that the simulation normally starts at t = 0 with zero fields as the initial condition, so there is implicitly a sharp turnon at t = 0 whether you specify it or not. endtime
[number
] The end time for the source; default is
infinity
(never turn off).

fluxregion
A fluxregion
object is used with addflux
below to specify a region in which Meep should accumulate the appropriate Fouriertransformed fields in order to compute a flux spectrum.
fluxregion
 A region (volume, plane, line, or point) in which to compute the integral of the Poynting vector of the Fouriertransformed 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 isAUTOMATIC
, 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 thedirection
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:

(shiftgeometricobject obj shiftvector)
 Translate
obj
by the 3vectorshiftvector
. 
(geometricobjectduplicates shiftvector minmultiple maxmultiple obj)
 Return a list of duplicates of
obj
, shifted by various multiples ofshiftvector
(fromminmultiple
tomaxmultiple
, inclusive, in steps of 1). 
(geometricobjectsduplicates shiftvector minmultiple maxmultiple objlist)
 Same as
geometricobjectduplicates
, except operates on a list of objects,objlist
. 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. 
(geometricobjectslatticeduplicates objlist [ ux uy uz ])
 Duplicates the objects in
objlist
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 [usertutorial.html tutorial].) The primitive cell to duplicate isux
byuy
byuz
, in units of the Cartesian basis vectors. These three parameters are optional; any that you do not specify are assumed to be1
. 
(pointinobject? point obj)
 Returns whether or not the given 3vector
point
is inside the geometric objectobj
. 
(pointinperiodicobject? point obj)
 As
pointinobject?
, but also checks translations of the given object by the lattice vectors. 
(displaygeometricobjectinfo indentby obj)
 Outputs some information about the given
obj
, indented byindentby
spaces.
Output file names
The output file names used by Meep, e.g. for HDF5 files, are automatically prefixed by the input variable filenameprefix
. If filenameprefix
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:
(getfilenameprefix)
 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 filenameprefix
to false
.
In addition to the filename prefix, you can also specify that all the output files be written into a newlycreated directory (if it does not yet exist). This is done by running:
(useoutputdirectory [dirname])
 Put output in a subdirectory, which is created if necessary. If the optional argument dirname is specified, 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"tstout"
.
Misc.
(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 thecenter
andsize
properties (just like e.g. ablock
object). If thesize
is not specified, it defaults to (0,0,0), i.e. a single point. (meeptime)
 Return the current simulation time (in simulation time units, not wallclock time!). (e.g. during a run function.)
 Occasionally, e.g. for termination conditions of the form time < T?, it is desirable to round the time to single precision in order to avoid small differences in roundoff error from making your results different by one timestep from machine to machine (a difference much bigger than roundoff error); in this case you can call
(meeproundtime)
instead, which returns the time rounded to single precision.
Field computations
Meep supports a large number of functions to perform computations on the fields. Currently, most of them are accessed via the lowerlevel C++/SWIG interface, but we are slowly adding simpler, higherlevel versions of them here.
(getfieldpoint c pt)
 Given a
component
orderivedcomponent
constantc
and avector3
pt
, returns the value of that component at that point. (getepsilonpoint pt)
 Equivalent to
(getfieldpoint Dielectric pt)
. (fluxinbox dir box)
 Given a
direction
constant, and ameep::volume*
, returns the flux (the integral of ) in that volume. Most commonly, you specify a volume that is a plane or a line, and a direction perpendicular to it, e.g.(fluxinbox X (volume (center 0) (size 0 1 1)))
. (electricenergyinbox box)
 Given a
meep::volume*
, returns the integral of the electricfield energy in the given volume. (If the volume has zero size along a dimension, a lowerdimensional integral is used.) (magneticenergyinbox box)
 Given a
meep::volume*
, returns the integral of the magneticfield energy in the given volume. (If the volume has zero size along a dimension, a lowerdimensional integral is used.) (fieldenergyinbox box)
 Given a
meep::volume*
, returns the integral of the electric+magneticfield energy in the given volume. (If the volume has zero size along a dimension, a lowerdimensional integral is used.)
Note that if you are at a fixed frequency and you use complex fields (Blochperiodic boundary conditions or fieldscomplex?=true
), then one half of the flux or energy integrals above corresponds to the timeaverage of the flux or energy for a simulation with real fields.
Often, you want the integration box to be the entire computational cell. A useful function to return this box, which you can then use for the box
arguments above, is (meepfieldstotalvolume fields)
, where fields
is the global variable (above) holding the current meep::fields
object.
One powerful feature is that you can supply an arbitrary function of position and various field components and ask Meep to integrate it over a given volume, find its maximum, or output it (via outputfieldfunction
, described later). This is done via the functions:
(integratefieldfunction cs func [where] [fieldsvar])
 Returns the integral of the complexvalued function
func
over themeep::geometric_volume
specified bywhere
(defaults to entire computational cell) for themeep::fields
specified byfieldsvar
(defaults tofields
).func
is a function of position (avector3
, its first argument) and zero or more field components specified bycs
: a list ofcomponent
constants.func
can be real or complexvalued.  (If any dimension of
where
is zero, that dimension is not integrated over. In this way you can specify one, two, or threedimensional integrals.) (maxabsfieldfunction cs func [where] [fieldsvar])
 As
integratefieldfunction
, but returns the maximum absolute value offunc
in the volumewhere
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.) See also Meep fieldfunction examples for illustrations of how to call integratefieldfunction
and maxabsfieldfunction
. See also Synchronizing the magnetic and electric fields if you want to do computations combining the electric and magnetic fields.
Occasionally, one wants to compute an integral that combines fields from two separate simulations (e.g. for nonlinear coupledmode calculations). This functionality is supported in Meep, as long as the two simulations have the same computational cell, the same resolution, the same boundary conditions and symmetries (if any), and the same PML layers (if any).
(integrate2fieldfunction fields2 cs1 cs2 func [where] [fieldsvar])
 Similar to
integratefieldfunction
, but takes additional parametersfields2
andcs2
.fields2
is ameep::fields*
object similar to the globalfields
variable (see below) specifying the fields from another simulation.cs1
is a list of components to integrate with fromfieldsvar
(defaults tofields
), as forintegratefieldfunction
, whilecs2
is a list of components to integrate fromfields2
. Similar tointegratefieldfunction
,func
is a function that returns an number given arguments consisting of: the position vector, followed by the values of the components specified bycs1
(in order), followed by the values of the components specified bycs2
(in order).
To get two fields in memory at once for integrate2fieldfunction
, the easiest way is to run one simulation within a given .ctl file, then save the results in another fields variable, then run a second simulation. This would look something like:
...set up and run first simulation... (define fields2 fields) ; save the fields in a variable (set! fields '()) ; prevent the fields from getting deallocated by resetmeep (resetmeep) ...set up and run second simulation...
It is also possible to timestep both fields simultaneously (e.g. doing one timestep of one simulation then one timestep of another simulation, and so on, but this requires you to call much lowerlevel functions like (meepfieldsstep fields)
.
Reloading parameters
Once the fields/simulation have been initialized, you can change the values of various parameters by using the following functions:
(resetmeep)
 Reset all of Meep's parameters, deleting the fields, structures, etcetera, from memory as if you had not run any computations.
(restartfields)
 Restart the fields at time zero, with zero fields. (Does not reset the Fourier transforms of the flux planes, which continue to be accumulated.)
(changekpoint! k)
 Change the
kpoint
(the Bloch periodicity). (changesources! newsources)
 Change the
sources
input variable tonewsources
, and changes the sources used for the current simulation.
(More to come...)
Flux spectra
Given a bunch of fluxregion
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:
(addflux fcen df nfreq fluxregions...)
 Add a bunch of
fluxregion
s to the current simulation (initializing the fields if they have not yet been initialized), telling Meep to accumulate the appropriate field Fourier transforms fornfreq
equally spaced frequencies covering the frequency rangefcendf/2
tofcen+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 addflux
via statements like:
(define transmission (addflux ...))
to store the flux object in a variable. addflux
initializes the fields if necessary, just like calling run
, so you should only call it after setting up your geometry
, sources
, pmllayers
, 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 addflux
, the Fourier transforms of the fields are accumulated automatically during timestepping by the run
functions. At any time, you can ask for Meep to print out the current flux spectrum via:
(displayfluxes fluxes...)
 Given a number of flux objects, this displays a commaseparated 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 lowerlevel if you have multiple flux regions corresponding to different frequency ranges, or have other special needs. (displayfluxes f1 f2 f3)
is actually equivalent to (displaycsv "flux" (getfluxfreqs f1) (getfluxes f1) (getfluxes f2) (getfluxes f3)
, where displaycsv
takes a bunch of lists of numbers and prints them as a commaseparated table, and we are calling two lowerlevel functions:
(getfluxfreqs flux)
 Given a flux object, returns a list of the frequencies that it is computing the spectrum for.
(getfluxes 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 Fouriertransformed fields from a "normalization" run and then load them into another run to be subtracted. This can be done via:
(saveflux filename flux)
 Save the Fouriertransformed fields corresponding to the given flux object in an HDF5 file of the given name (without the ".h5" suffix) (the current filenameprefix is prepended automatically).
(loadflux filename flux)
 Load the Fouriertransformed 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 filenameprefix is prepended automatically). You must load from a file that was saved by
saveflux
in a simulation of the same dimensions (for both the computational cell and the flux regions) with the same number of processors. (loadminusflux filename flux)
 As
loadflux
, but negates the Fouriertransformed fields after they are loaded. This means that they will be subtracted from any future field Fourier transforms that are accumulated. (scalefluxfields s flux)
 Scale the Fouriertransformed fields in
flux
by the complex numbers
. e.g.loadminusflux
is equivalent toloadflux
followed byscalefluxfields
withs=1
.
Force spectra
Very similar to flux spectra, you can also compute force spectra: forces on an object as a function of frequency, computed by Fourier transforming the fields and integrating the vacuum Maxwell stress tensor
over a surface S via . We recommend that you normally only evaluate the stress tensor over a surface lying in vacuum, as the interpretation and definition of the stress tensor in arbitrary media is often problematic (the subject of extensive and controversial literature). (It is fine if the surface encloses an object made of arbitrary materials, as long as the surface itself is in vacuum.)
See also the optical forces Meep tutorial.
Most commonly, you will want to normalize the force spectrum in some way, just as for flux spectra. Most simply, you could divide two different force spectra to compute the ratio of forces on two objects. Often, you will divide a force spectrum by a flux spectrum, to divide the force F by the incident power P on an object, in order to compute the useful dimensionless ratio Fc/P (where c=1 in Meep units). For example, it is a simple exercise to show that the force F on a perfectly reflecting mirror with normalincident power P satisfies Fc/P=2, and for a perfectly absorbing (black) surface Fc/P=1.
The usage is similar to the flux spectra: you define a set of forceregion
objects telling Meep where it should compute the Fouriertransformed fields and stress tensors, and call addforce
to add these regions to the current simulation over a specified frequency bandwidth, and then use displayforces
to display the force spectra at the end. There are also saveforce
, loadforce
, and loadminusforce
functions that you can use to subtract the fields from two simulation, e.g. in order to compute just the force from scattered fields, similar to the flux spectra. These types and functions are defined as follows:
forceregion
 A region (volume, plane, line, or point) in which to compute the integral of the stress tensor of the Fouriertransformed fields. Its properties are:
center
[vector3
] The center of the force region (no default).
size
[vector3
] The size of the force region along each of the coordinate axes; default is (0,0,0) (a single point).
direction
[direction
constant] The direction of the force that you wish to compute (e.g.
X
,Y
, etcetera). Unlikefluxregion
, you must specify this explicitly, because there is not generally any relationship between the direction of the force and the orientation of the force region. weight
[cnumber
] A weight factor to multiply the force by when it is computed; default is
1.0
.
In most circumstances, you should define a set of forceregion
s whose union is a closed surface (lying in vacuum and enclosing the object that is experiencing the force).
(addforce fcen df nfreq forceregions...)
 Add a bunch of
forceregion
s to the current simulation (initializing the fields if they have not yet been initialized), telling Meep to accumulate the appropriate field Fourier transforms fornfreq
equally spaced frequencies covering the frequency rangefcendf/2
tofcen+df/2
. Return a force object, which you can pass to the functions below to get the force spectrum, etcetera.
As for flux regions, you normally use addforce
via statements like:
(define Fx (addforce ...))
to store the flux object in a variable. addforce
initializes the fields if necessary, just like calling run
, so you should only call it after setting up your geometry
, sources
, pmllayers
, etcetera. You can create as many force objects as you want, e.g. to look at forces on different objects, in different directions, 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 force region, multiplied by the number of electric and magnetic field components required to get the stress 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 addforce
, the Fourier transforms of the fields are accumulated automatically during timestepping by the run
functions. At any time, you can ask for Meep to print out the current force spectrum via:
(displayforces forces...)
 Given a number of force objects, this displays a commaseparated table of frequencies and force spectra, prefixed by "force1:" or similar (where the number is incremented after each run). All of the forces should be for the same
fcen
/df
/nfreq
. The first column are the frequencies, and subsequent columns are the force spectra.
You might have to do something lowerlevel if you have multiple force regions corresponding to different frequency ranges, or have other special needs. (displayforces f1 f2 f3)
is actually equivalent to (displaycsv "force" (getforcefreqs f1) (getforces f1) (getforces f2) (getforces f3)
, where displaycsv
takes a bunch of lists of numbers and prints them as a commaseparated table, and we are calling two lowerlevel functions:
(getforcefreqs flux)
 Given a force object, returns a list of the frequencies that it is computing the spectrum for.
(getforces flux)
 Given a force object, returns a list of the current force spectrum that it has accumulated.
As described in the Meep tutorial, to compute the force from scattered fields often want to save the Fouriertransformed fields from a "normalization" run and then load them into another run to be subtracted. This can be done via:
(saveforce filename force)
 Save the Fouriertransformed fields corresponding to the given force object in an HDF5 file of the given name (without the ".h5" suffix) (the current filenameprefix is prepended automatically).
(loadforce filename force)
 Load the Fouriertransformed fields into the given force object (replacing any values currently there) from an HDF5 file of the given name (without the ".h5" suffix) (the current filenameprefix is prepended automatically). You must load from a file that was saved by
saveforce
in a simulation of the same dimensions (for both the computational cell and the force regions) with the same number of processors. (loadminusforce filename force)
 As
loadforce
, but negates the Fouriertransformed fields after they are loaded. This means that they will be subtracted from any future field Fourier transforms that are accumulated.
LDOS spectra
Meep can also calculate the LDOS (local density of states) spectrum, as described in the Meep LDOS tutorial. To do this, you simply pass the following step function to your run
command:
(dftldos fcen df nfreq)
 Compute the power spectrum of the sources (usually a single point dipole source), normalized to correspond to the LDOS, in a frequency bandwith
df
centered atfcen
, atnfreq
frequency points.
The resulting spectrum is outputted as commadelimited text, prefixed by ldos:,
, and is also stored in the dftldosdata
global variable (after the run
is complete).
Analytically, the perpolarization LDOS is exactly proportional to the power radiated by an oriented pointdipole current, p(t), at a given position in space. For a more mathematical treatment of the theory behind the LDOS, we refer you to the relevant discussion in chapter 4 of our book, but for now we simply give the definition:
where the normalization is necessary for obtaining the power exerted by a unitamplitude dipole (assuming linear materials), and hats denote Fourier transforms. It is this quantity that is computed by the dftldos
command for a single dipole source. (For a volumetric source, the numerator and denominator are both integrated over the current volume, but "LDOS" computation is less meaningful in this case.)
Neartofarfield spectra
Meep (1.3 or later) can compute a "neartofarfield transformation" in the frequency domain as described in this tutorial: given the fields on a "near" bounding surface inside the computational cell, it can compute the fields arbitrarily far away using an analytical transformation, assuming that the "near" surface and the "far" region lie in a single homogeneous nonperiodic 2d or 3d region. That is, in a simulation surrounded by PML that absorbs outgoing waves, the neartofarfield feature can compute the fields outside the computational cell as if the outgoing waves had not been absorbed (i.e. in the fictitious infinite open volume). Moreover, this operation is performed on the Fouriertransformed fields: like the flux and force spectra above, you specify a set of desired frequencies, Meep accumulates the Fourier transforms, and then Meep computes the fields at each frequency for the desired farfield points.
(This is based on the principle of equivalence — given the Fouriertransformed tangential fields on the "near" surface, Meep computes equivalent currents and convolves them with the analytical Green's functions in order to compute the fields at any desired point in the "far" region.)
There are three steps to using the neartofarfield feature: first, define the "near" surface(s) as a set of surfaces capturing all outgoing radiation in the desired direction(s); second, run the simulation, typically with a pulsed source, to allow Meep to accumulate the Fourier transforms on the near surface(s); third, tell Meep to compute the far fields at any desired points (optionally saving the far fields from a grid of points to an HDF5 file). To define the near surfaces, use:
(addnear2far fcen df nfreq near2farregions...)
 Add a bunch of
near2farregion
s to the current simulation (initializing the fields if they have not yet been initialized), telling Meep to accumulate the appropriate field Fourier transforms fornfreq
equally spaced frequencies covering the frequency rangefcendf/2
tofcen+df/2
. Return a near2far object, which you can pass to the functions below to get the far fields.
Each near2farregion
is identical to fluxregion
except for the name: in 3d, these give a set of planes (important: all these "near surfaces" must lie in a single homogeneous material with isotropic ε and μ — and they should not lie in the PML regions) surrounding the source(s) of outgoing radiation that you want to capture and convert to a far field. Ideally, these should form a closed surface, but in practice it is sufficient for the near2farregion
s to capture all of the radiation in the direction of the farfield points. Important: as for flux computations, each near2farregion
should be assigned a weight
of ±1 indicating the direction of the outward normal relative to the +coordinate direction. So, for example, if you have six regions defining the six faces of a cube, i.e. the faces in the +x, x, +y, y, +z, and z directions, then they should have weights +1, 1, +1, 1, +1, and 1 respectively. Note that, neglecting discretization errors, all nearfield surfaces that enclose the same outgoing fields are equivalent and will yield the same far fields (with a discretizationinduced difference that vanishes with increasing resolution etc.).
After the simulation run
(usually for a pulsed source) is complete (so that the fields have decayed away and the Fourier transforms have finished accumulating), you can compute the far fields.
(getfarfield near2far x)
 Given a
vector3
pointx
(which can lie anywhere outside the nearfield surface, including outside the computational cell) and a near2far object, returns the computed (Fouriertransformed) "far" fields at x as list of length 6nfreq
, consisting of fields (Ex1,Ey1,Ez1,Hx1,Hy1,Hz1,Ex2,Ey2,Ez2,Hx2,Hy2,Hz2,...) for the frequencies 1,2,…,nfreq
. (outputfarfields near2far fname where resolution)
 Given an HDF5 file name
fname
(does not include the.h5
suffix), avolume
given bywhere
(may be 0d, 1d, 2d, or 3d), and aresolution
(in grids points / distance a), outputs the far fields inwhere
(which may lie outside the computational cell) in a grid with the given resolution (which may differ from the FDTD grid resolution) to the HDF5 file as a set of twelve array datasetsex.r
,ex.i
, ...,hz.r
,hz.i
, giving the real and imaginary parts of the Fouriertransformed E and H fields on this grid. Each dataset is an nx×ny×nz×nfreq 4dimensional array of space×frequency (although dimensions that =1 are omitted).
Note that far fields have the same units and scaling as the Fourier transforms of the fields, and hence cannot be directly compared to timedomain fields. In practice, it is easiest to use the far fields in computations where overall scaling (units) cancel out or are irrelevant, e.g. to compute the fraction of the far fields in one region vs. another region.
For a scatteredfield computation, you often want to separate the scattered and incident fields. Just as is described in the Meep tutorial for flux computations, you can do this by saving the Fouriertransformed incident from a "normalization" run and then load them into another run to be subtracted. This can be done via:
(savenear2far filename near2far)
 Save the Fouriertransformed fields corresponding to the given near2far object in an HDF5 file of the given name (without the ".h5" suffix) (the current filenameprefix is prepended automatically).
(loadnear2far filename near2far)
 Load the Fouriertransformed fields into the given near2far object (replacing any values currently there) from an HDF5 file of the given name (without the ".h5" suffix) (the current filenameprefix is prepended automatically). You must load from a file that was saved by
savenear2far
in a simulation of the same dimensions (for both the computational cell and the near2far regions) with the same number of processors. (loadminusnear2far filename near2far)
 As
loadnear2far
, but negates the Fouriertransformed fields after they are loaded. This means that they will be subtracted from any future field Fourier transforms that are accumulated. (scalenear2farfields s near2far)
 Scale the Fouriertransformed fields in
near2far
by the complex numbers
. e.g.loadminusnear2far
is equivalent toloadnear2far
followed byscalenear2farfields
withs=1
.
Frequencydomain solver
Meep contains a frequencydomain solver that directly computes the fields produced in a geometry in response to a constantfrequency source, using an iterative linear solver instead of timestepping. Preliminary tests have shown that in many instances, this solver converges much faster than simply running an equivalent time domain simulation with a continuous wave source, timestepping until all transient effects from the source turnon have disappeared, especially if the fields are desired to a very high accuracy. To use it, simply define a continuoussrc
with the desired frequency, initialize the fields and geometry via (initfields)
, and then:
(meepfieldssolvecw fields tol maxiters L)
After the fields
variable (a global variable pointing to the meep::fields*
object initialized by initfields
, see above), the next two parameters to the frequencydomaine solver are the tolerance tol
for the iterative solver (10^{−8}, by default) and a maximum number of iterations maxiters
(10^{4}, by default). Finally, there is a parameter L
that determines a tradeoff between memory and work per step and convergence rate of the iterative algorithm biCGSTAB(L) that is used; larger values of L
of will often lead to faster convergence at the expense of more memory and more work per iteration. The default is L
=2, and normally a value ≥ 2 should be used.
The frequencydomain solver supports arbitrary geometries, PML, boundary conditions, symmetries, parallelism, conductors, and arbitrary nondispersive materials. Lorentz–Drude dispersive materials are not currently supported in the frequencydomain solver, but since you are solving at a known fixed frequency rather than timestepping, you should be able to pick conductivities etcetera in order to obtain any desired complex ε and μ at that frequency.
The frequencydomain solver requires you to use complexvalued fields, via (set! forcecomplexfields? true)
.
After meepfieldssolvecw
completes, it should be as if you had just run the simulation for an infinite time with the source at that frequency. You can call the various fieldoutput functions and so on as usual at this point.
Run and step functions
The actual work in Meep is performed by run functions, which timestep 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 atbeginning, duringsources, etcetera which cause them to only be called at certain times, etcetera, instead of at every time step.
A common point of confusion is described in the article: The run function is not a loop. Please read this article if you want to make Meep do some customized action on each time step, as many users make the same mistake. What you really want to in that case is to write a step function, as described below.
Run functions
The following run functions are available. (You can also write your own, using the lowerlevel C++/SWIG functions, but these should suffice for most needs.)
(rununtil cond?/time stepfunctions...)
 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. (runsources stepfunctions...)
 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.
(runsources+ cond?/time stepfunctions...)
 As
runsources
, 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?)
returnstrue
.
In particular, a useful first argument to runsources+
or rununtil
is often given by (as in the Meep tutorial):
(stopwhenfieldsdecayed dT c pt decayby)
 Return a
cond?
function, suitable for passing torununtil
/runsources+
, that examines the componentc
(e.g.Ex
, etc.) at the pointpt
(avector3
) and keeps running until its absolute value squared has decayed by at leastdecayby
from its maximum previous value. In particular, it keeps incrementing the run time bydT
(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
decayby
very small, you may need to increase thecutoff
property of your source(s), to decrease the amplitude of the small highfrequency 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
(runkpoint T k)
 Given a
vector3
k
, runs a simulation for each k point (i.e. specifying Blochperiodic boundary conditions) and extracts the eigenfrequencies, 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 additionalT
time units. It will runharminv
(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. (runkpoints T kpoints)
 Given a list
kpoints
of k vectors, runsrunkpoint
for each one, and returns a list of lists of frequencies (one list of frequencies for each k). Also prints out a commadelimited list of frequencies, prefixed byfreqs:
, and their imaginary parts, prefixed byfreqsim:
. (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 storageconsuming.
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. in 3d.
The predefined output functions are:
outputepsilon
 Output the dielectric function (relative permittivity) ε. Note that this only outputs the frequencyindependent part of ε (the limit).
outputmu
 Output the relative permeability function μ. Note that this only outputs the frequencyindependent part of μ (the limit).
outputhpwr
 Output the magneticfield energy density
outputdpwr
 Output the electricfield energy density
outputtotpwr
 Output the total electric and magnetic energy density. Note that you might want to wrap this step function in
synchronizedmagnetic
to compute it more accurately; see Synchronizing the magnetic and electric fields. outputXfieldx
,outputXfieldy
,outputXfieldz
,outputXfieldr
,outputXfieldp
 Output the x, y, z, r, or φ component respectively, of the field X, where X is either
h
,b
,e
,d
, ors
for the magnetic, electric, displacement, or Poynting field, respectively. If the field is complex, outputs two datasets, e.g.ex.r
andex.i
, within the same HDF5 file for the real and imaginary parts, respectively. Note that for outputting the Poynting field, you might want to wrap the step function insynchronizedmagnetic
to compute it more accurately; see Synchronizing the magnetic and electric fields. outputXfield
 Outputs all the components of the field X, where X is either
h
,b
,e
,d
, ors
as above, to an HDF5 file. (That is, the different components are stored as different datasets within the same file.) (outputpng component h5topngoptions)
 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.  It is often useful to use the
h5topng C
orA
options to overlay the dielectric function when outputting fields. To do this, you need to know the name of the dielectricfunction.h5
file (which must have been previously output by outputepsilon). To make this easier, a builtin shell variable$EPS
is provided which refers to the lastoutput dielectricfunction.h5
file. So, for example(outputpng Ez "C $EPS")
will output the E_{z} field and overlay the dielectric contours. (outputpng+h5 component h5topngoptions)
 Like
outputpng
, but also outputs the.h5
file for the component. (In contrast,outputpng
deletes the.h5
when it is done.)
More generally, it is possible to output an arbitrary function of position and zero or more field components, similar to the integratefieldfunction
described above. This is done by:
 (outputfieldfunction name cs func)
 Output the field function
func
to an HDF5 file in the datasets namedname.r
andname.i
(for the real and imaginary parts). Similar tointegratefieldfunction
,func
is a function of position (avector3
) and the field components corresponding tocs
: a list ofcomponent
constants.  (outputrealfieldfunction name cs func)
 As
outputfieldfunction
, but only outputs the real part offunc
to the dataset given by the stringname
.
See also Meep fieldfunction examples. See also Synchronizing the magnetic and electric fields if you want to do computations combining the electric and magnetic fields.
Harminv
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 pointpt
(avector3
). Then, at the end of the run, it uses Harminv to look for modes in the given frequency range (centerfcen
and widthdf
), printing the results to standard output (prefixed byharminv:
) as commadelimited text, and also storing them to the variableharminvresults
. The optional argumentmaxbands
is the maximum number of modes to search for; defaults to100
.
Important: normally, you should only use harminv
to analyze data after the sources are off. Wrapping it in (aftersources (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:
The results are stored in the list harminvresults
, 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:
(harminvfreq result)
 Return the complex frequency ω (in the usual Meep 2πc units).
(harminvfreqre result)
 Return the real part of the frequency ω.
(harminvfreqim result)
 Return the imaginary part of the frequency ω.
(harminvQ result)
 Return dimensionless lifetime, or "quality factor", Q, defined as .
(harminvamp result)
 Return the complex amplitude a.
(harminverr 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 harminvfreqre harminvresults)
gives you a list of the real parts of the frequencies, using the Scheme builtin map
.
Stepfunction modifiers
Rather than writing a brandnew 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.)
Miscellaneous stepfunction modifiers
(combinestepfuncs stepfunctions...)
 Given zero or more step functions, return a new step function that (on each step) calls all of the passed step functions.
(synchronizedmagnetic stepfunctions...)
 Given zero or more step functions, return a new step function that (on each step) calls all of the passed step functions with the magnetic field synchronized in time with the electric field. See Synchronizing the magnetic and electric fields.
Controlling when a step function executes
(whentrue cond? stepfunctions...)
 Given zero or more step functions and a condition function
cond?
( a function of no arguments), evaluate the step functions whenever(cond?)
returnstrue
. (whenfalse cond? stepfunctions...)
 Given zero or more step functions and a condition function
cond?
( a function of no arguments), evaluate the step functions whenever(cond?)
returnsfalse
. (atevery dT stepfunctions...)
 Given zero or more step functions, evaluates them at every time interval of
dT
units (rounded up to the next time step). (aftertime T stepfunctions...)
 Given zero or more step functions, evaluates them only for times after a
T
time units have elapsed from the start of the run. (beforetime T stepfunctions...)
 Given zero or more step functions, evaluates them only for times before a
T
time units have elapsed from the start of the run. (attime T stepfunctions...)
 Given zero or more step functions, evaluates them only once, after a
T
time units have elapsed from the start of the run. (aftersources stepfunctions...)
 Given zero or more step functions, evaluates them only for times after all of the sources have turned off.
(aftersources+ T stepfunctions...)
 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. (duringsources stepfunctions...)
 Given zero or more step functions, evaluates them only for times before all of the sources have turned off.
(atbeginning stepfunctions...)
 Given zero or more step functions, evaluates them only once, at the beginning of the run.
(atend stepfunctions...)
 Given zero or more step functions, evaluates them only once, at the end of the run.
Modifying HDF5 output
(involume v stepfunctions...)
 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 thevolume
function). (inpoint pt stepfunctions...)
 Given zero or more step functions, modifies any output functions among them to only output a single point of data, at
pt
(avector3
). (toappended filename stepfunctions...)
 Given zero or more step functions, modifies any output functions among them to append their data to datasets in a single newlycreated 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. (withprefix prefix stepfunctions...)
 Given zero or more step functions, modifies any output functions among them to prepend the string
prefix
to the file names (much likefilenameprefix
, above).
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 (mystep) (print "Hello world!\n"))
If one then does (rununtil 100 mystep)
, 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 timestepping, or 'finish
when it is called at the end of the run.
Lowlevel 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 wrappergenerator 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 addflux
. To initialize them separately, you can call (initfields)
manually, or (initstructure kpoint)
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:
(initfields) (define myfields fields) (set! fields '()) (resetmeep)
and then change the geometry etc. and rerun (initfields)
. 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 withmeep
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 timestep, is exposed to Scheme asmeepfieldsstep
. Moreover, you pass the object as the first argument in the Scheme wrapper. e.g.f.step()
becomes(meepfieldsstep f)
.  To call the C++ constructor for a type, you use
newtype
. e.g.(newmeepfields ...arguments...)
returns a newmeep::fields
object. Conversely, to call the destructor and deallocate an object, you usedeletetype
; most of the time, this is not necessary because objects are automatically garbagecollected.
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 initstructure
, but if you want to pass vector arguments to C++ before that time you should call (requiredimensions!)
, which infers the dimensions from the geometrylattice
, kpoint
, and dimensions
variables.