Python User Interface
This page is a listing of the functions exposed by the Python interface. For a gentler introduction, see Tutorial/Basics. Note that this page is not a complete listing of all functions. In particular, because of the SWIG wrappers, every function in the C++ interface is accessible from the Python module, but not all of these functions are documented or intended for end users. See also the instructions for Parallel Meep.
The Python API functions and classes can be found in the meep
module, which should be installed in your Python system by Meep's make install
script. If you installed into a nonstandard location (e.g. your home directory), you may need to set the PYTHONPATH
environment variable as documented in Building From Source. You typically import the meep
module in Python via import meep as mp
.
 Python User Interface
 Predefined Variables
 Constants (Enumerated Types)
 The Simulation Class
 Simulation
 Output File Names
 Simulation Time
 Field Computations
 Reloading Parameters
 Flux Spectra
 Mode Decomposition
 DiffractedPlanewave
 Energy Density Spectra
 Force Spectra
 LDOS spectra
 NeartoFarField Spectra
 Load and Dump Structure
 FrequencyDomain Solver
 FrequencyDomain Eigensolver
 GDSII Support
 Data Visualization
 Run and Step Functions
 Predefined Step Functions
 StepFunction Modifiers
 Writing Your Own Step Functions
 LowLevel Functions
 Class Reference
 Medium
 Susceptibility
 LorentzianSusceptibility
 DrudeSusceptibility
 MultilevelAtom
 Transition
 NoisyLorentzianSusceptibility
 NoisyDrudeSusceptibility
 GyrotropicLorentzianSusceptibility
 GyrotropicDrudeSusceptibility
 GyrotropicSaturatedSusceptibility
 Vector3
 GeometricObject
 Sphere
 Cylinder
 Wedge
 Cone
 Block
 Ellipsoid
 Prism
 Matrix
 Symmetry
 Rotate2
 Rotate4
 Mirror
 Identity
 PML
 Absorber
 Source
 SourceTime
 EigenModeSource
 GaussianBeamSource
 ContinuousSource
 GaussianSource
 CustomSource
 FluxRegion
 EnergyRegion
 ForceRegion
 Volume
 DftObj
 DftFlux
 DftForce
 DftNear2Far
 DftEnergy
 DftFields
 Animate2D
 Harminv
 Verbosity
 Miscellaneous Functions Reference
Predefined Variables
These are available directly via the meep
package.
air
, vacuum
[Medium
class ]
—
Two aliases for a predefined material type with a dielectric constant of 1.
perfect_electric_conductor
or metal
[Medium
class ]
—
A predefined material type corresponding to a perfect electric conductor at the boundary of which the parallel electric field is zero. Technically, .
perfect_magnetic_conductor
[Medium
class ]
—
A predefined material type corresponding to a perfect magnetic conductor at the boundary of which the parallel magnetic field is zero. Technically, .
inf
[number
]
—
A big number (10^{20}) to use for "infinite" dimensions of objects.
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 (which are available directly via the meep
package):
direction
constants
—
Specify a direction in the grid. One of X
, Y
, Z
, R
, P
for , , , , , respectively.
side
constants
—
Specify particular boundary in the positive High
(e.g., +X
) or negative Low
(e.g., X
) direction.
boundary_condition
constants
—
Metallic
(i.e., zero electric field) or Magnetic
(i.e., zero magnetic field).
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 , , , , , , , , , , , , , , , , , , , , ε, μ, 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 , , , , (components of the Poynting vector ), , , , respectively.
The Simulation Class
The Simulation
class is the primary abstraction of the highlevel interface. Minimally, a simulation script amounts to passing the desired keyword arguments to the Simulation
constructor and calling the run
method on the resulting instance.
Simulation
class Simulation(object):
The Simulation
class contains all the attributes that you can set to
control various parameters of the Meep computation.
def __init__(self,
cell_size,
resolution,
geometry=[],
sources=[],
eps_averaging=True,
dimensions=3,
boundary_layers=[],
symmetries=[],
force_complex_fields=False,
default_material=Medium(),
m=0,
k_point=False,
kz_2d='complex',
extra_materials=[],
material_function=None,
epsilon_func=None,
epsilon_input_file='',
progress_interval=4,
subpixel_tol=0.0001,
subpixel_maxeval=100000,
ensure_periodicity=True,
num_chunks=0,
Courant=0.5,
accurate_fields_near_cylorigin=False,
filename_prefix=None,
output_volume=None,
output_single_precision=False,
load_structure='',
geometry_center=Vector3<0.0, 0.0, 0.0>,
force_all_components=False,
split_chunks_evenly=True,
chunk_layout=None,
collect_stats=False):
All Simulation
attributes are described in further detail below. 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
, complex
, and string
are defined by
Python. Vector3
is a meep
class.

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). 
geometry_center
[Vector3
class ] — Specifies the coordinates of the center of the cell. Defaults to (0, 0, 0), but changing this allows you to shift the coordinate system used in Meep (for example, to put the origin at the corner). Passinggeometry_center=c
is equivalent to adding thec
vector to the coordinates of every other object in the simulation, i.e.c
becomes the new origin that other objects are defined with respect to. 
sources
[ list ofSource
class ] — Specifies the current sources to be present in the simulation. Defaults to none (empty list). 
symmetries
[ list ofSymmetry
class ] — Specifies the spatial symmetries (mirror or rotation) to exploit in the simulation. Defaults to none (empty list). The symmetries must be obeyed by both the structure and the sources. See also Exploiting Symmetry. 
boundary_layers
[ list ofPML
class ] — Specifies the PML absorbing boundary layers to use. Defaults to none. 
cell_size
[Vector3
] — Specifies the size of the cell which is centered on the origin of the coordinate system. Any sizes of 0 imply a reduceddimensionality calculation. Strictly speaking, the dielectric function is taken to be uniform along that dimension. A 2d calculation is especially optimized. Seedimensions
below. Note: because Maxwell's equations are scale invariant, you can use any units of distance you want to specify the cell size: nanometers, microns, centimeters, etc. However, it is usually convenient to pick some characteristic lengthscale of your problem and set that length to 1. See also Units. Required argument (no default). 
default_material
[Medium
class ] — Holds the default material that is used for points not in any object of the geometry list. Defaults toair
(ε=1). This can also be a NumPy array that defines a dielectric function much likeepsilon_input_file
below (see below). If you want to use a material function as the default material, use thematerial_function
keyword argument (below). 
material_function
[ function ] — A Python function that takes aVector3
and returns aMedium
. See also Material Function. Defaults toNone
. 
epsilon_func
[ function ] — A Python function that takes aVector3
and returns the dielectric constant at that point. See also Material Function. Defaults toNone
. 
epsilon_input_file
[string
] — If this string is not empty (the default), then it should be the name of an HDF5 file whose first/only dataset defines a scalar, realvalued, frequencyindependent 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 ofdefault_material
(i.e. where there are nogeometry
objects). The grid of the epsilon file dataset need not match the computational grid; it is scaled and/or linearly interpolated as needed to map the file onto the cell. The structure is warped if the proportions of the grids do not match. Note: the file contents only override the ε property of thedefault_material
, whereas other properties (μ, susceptibilities, nonlinearities, etc.) ofdefault_material
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 whencell_size
in the direction is0
. Ifdimensions
is the special value ofCYLINDRICAL
, then cylindrical coordinates are used and the and dimensions are interpreted as and , respectively. Ifdimensions
is 1, then the cell must be along the direction and only and field components are permitted. Ifdimensions
is 2, then the cell must be in the plane. 
m
[number
] — ForCYLINDRICAL
simulations, specifies that the angular dependence of the fields is of the form (default ism=0
). If the simulation cell includes the origin , thenm
must be an integer. 
accurate_fields_near_cylorigin
[boolean
] — ForCYLINDRICAL
simulations with m > 1, compute more accurate fields near the origin at the expense of requiring a smaller Courant factor. Empirically, when this option is set toTrue
, a Courant factor of roughly or smaller seems to be needed. Default isFalse
, in which case the , , and 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 . 
resolution
[number
] — Specifies the computational grid resolution in pixels per distance unit. Required argument. No default. 
k_point
[False
orVector3
] — IfFalse
(the default), then the boundaries are perfect metallic (zero electric field). If aVector3
, then the boundaries are Blochperiodic: the fields at one side are times the fields at the other side, separated by the lattice vector . A nonzeroVector3
will produce complex fields. Thek_point
vector is specified in Cartesian coordinates in units of 2π/distance. Note: this is different from MPB, equivalent to taking MPB'sk_points
through its functionreciprocal>cartesian
. 
kz_2d
["complex"
,"real/imag"
, or"3d"
] — A 2d cell (i.e.,dimensions=2
) combined with ak_point
that has a nonzero component in would normally result in a 3d simulation with complex fields. However, by default (kz_2d="complex"
), Meep will use a 2d computational cell in which is incorporated as an additional term in Maxwell's equations, which still results in complex fields but greatly improved performance. Settingkz_2d="3d"
will instead use a 3d cell that is one pixel thick (with Blochperiodic boundary conditions), which is considerably more expensive. The third possibility,kz_2d="real/imag"
, saves an additional factor of two by storing some field components as purely real and some as purely imaginary in a "real" field, but this option requires some care to use. See 2d Cell with OutofPlane Wavevector. 
ensure_periodicity
[boolean
] — IfTrue
(the default) and if the boundary conditions are periodic (k_point
is notFalse
), then the geometric objects are automatically repeated periodically according to the lattice vectors which define the size of the cell. 
eps_averaging
[boolean
] — IfTrue
(the default), then subpixel averaging is used when initializing the dielectric function. For simulations involving a material function,eps_averaging
isFalse
(the default) and must be enabled in which case the input variablessubpixel_maxeval
(default 10^{4}) andsubpixel_tol
(default 10^{4}) specify the maximum number of function evaluations and the integration tolerance for the adaptive numerical integration. Increasing/decreasing these, respectively, will cause a more accurate but slower computation of the average ε with diminishing returns for the actual FDTD error. Disabling subpixel averaging will lead to staircasing effects and irregular convergence. 
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 thek_point
is nonzero or ifm
is nonzero. However, by settingforce_complex_fields
toTrue
, Meep will always use complex fields. 
force_all_components
[boolean
] — By default, in a 2d simulation Meep uses only the field components that might excited by your current sources: either the inplane or outofplane polarization, depending on the source. (Both polarizations are excited if you use multiple source polarizations, or if an anisotropic medium is present that couples the two polarizations.) In rare cases (primarily for combining results of multiple simulations with differing polarizations), you might want to force it to simulate all fields, even those that remain zero throughout the simulation, by settingforce_all_components
toTrue
. 
filename_prefix
[string
] — A string prepended to all output filenames (e.g., for HDF5 files). IfNone
(the default), then Meep constructs a default prefix based on the current Python filename ".py" replaced by "" (e.g.foo.py
uses a"foo"
prefix). You can get this prefix by callingget_filename_prefix
. 
Courant
[number
] — Specify the Courant factor which relates the time step size to the spatial discretization: . Default is 0.5. For numerical stability, the Courant factor must be at most , where is the minimum refractive index (usually 1), and in practice should be slightly smaller. 
output_volume
[Volume
class ] — Specifies the default region of space that is output by the HDF5 output functions (below); see also theVolume
class which managesmeep::volume*
objects. Default isNone
, which means that the whole cell is output. Normally, you should use thein_volume(...)
function to modify the output volume instead of settingoutput_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 toTrue
(default isFalse
) 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. 
extra_materials
[ list ofMedium
class ] — By default, Meep turns off support for material dispersion (susceptibilities or conductivity) or nonlinearities if none of the objects ingeometry
have materials with these properties — since they are not needed, it is faster to omit their calculation. This doesn't work, however, if you use amaterial_function
: materials via a userspecified function of position instead of just geometric objects. If your material function only returns a nonlinear material, for example, Meep won't notice this unless you tell it explicitly viaextra_materials
.extra_materials
is a list of materials that Meep should look for in the cell in addition to any materials that are specified by geometric objects. You should list any materials other than scalar dielectrics that are returned bymaterial_function
here. 
load_structure
[string
] — If not empty, Meep will load the structure file specified by this string. The file must have been created bymp.dump_structure
. Defaults to an empty string. See Load and Dump Structure for more information. 
chunk_layout
[string
orSimulation
instance] — This will cause theSimulation
to use the chunk layout described by either an h5 file (created bySimulation.dump_chunk_layout
) or anotherSimulation
. See Load and Dump Structure for more information.
The following require a bit more understanding of the inner workings of Meep to use. See also SWIG Wrappers.

structure
[meep::structure*
] — Pointer to the current structure being simulated; initialized by_init_structure
which is called automatically byinit_sim()
which is called automatically by any of the run functions. The structure initialization is handled by theSimulation
class, and most users will not need to call_init_structure
. 
fields
[meep::fields*
] — Pointer to the current fields being simulated; initialized byinit_sim()
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. 
split_chunks_evenly
[boolean
] — WhenTrue
(the default), the work per chunk is not taken into account when splitting chunks up for multiple processors. The cell is simply split up into equal chunks (with the exception of PML regions, which must be on their own chunk). WhenFalse
, Meep attempts to allocate an equal amount of work to each processor, which can increase the performance of parallel simulations.
def run(self, *step_funcs, **kwargs):
def run(step_functions..., until=condition/time):
def run(step_functions..., until_after_sources=condition/time):
run(step_functions..., until=condition/time)
Run the simulation until a certain time or condition, calling the given step
functions (if any) at each timestep. The keyword argument until
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. until
can also be a list of stopping conditions which may include a number
and additional functions.
run(step_functions..., until_after_sources=condition/time)
Run the simulation until all sources have turned off, calling the given step
functions (if any) at each timestep. The keyword argument until_after_sources
is
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 condition
returns True
. Like until
above, until_after_sources
can take a list of
stopping conditions.
Output File Names
The output filenames used by Meep, e.g. for HDF5 files, are automatically prefixed by the
input variable filename_prefix
. If filename_prefix
is None
(the default), however,
then Meep constructs a default prefix based on the current Python filename with ".py"
replaced by ""
: e.g. test.py
implies a prefix of "test"
. You can get this prefix,
or set the output folder, with these methods of the Simulation
class:
def get_filename_prefix(self):
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 the
empty string ''
.
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 calling Simulation.use_output_directory([dirname])
def use_output_directory(self, dname=''):
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 Python file name with ".py"
replaced by "out"
:
e.g. test.py
implies a directory of "testout"
.
Simulation Time
The Simulation
class provides the following timerelated methods:
def meep_time(self):
Return the current simulation time in simulation time units (e.g. during a run function). This is not the wallclock time.
Occasionally, e.g. for termination conditions of the form , 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 Simulation.round_time()
instead, which returns the time
rounded to single precision.
def print_times(self):
Call after running a simulation to print the times spent on various types of work. Example output:
Field time usage:
connecting chunks: 0.0819176 s +/ 0.000428381 s
time stepping: 0.198949 s +/ 0.0225551 s
communicating: 0.410577 s +/ 0.278853 s
outputting fields: 0.512352 s +/ 0.0238399 s
Fourier transforming: 0.0738274 s +/ 0.0967926 s
everything else: 0.324933 s +/ 0.377573 s
def time_spent_on(self, time_sink):
Return a list of times spent by each process for a type of work time_sink
which
can be one of nine integer values 0
8
: (0
) connecting chunks, (1
) time
stepping, (2
) boundaries, (3
) MPI/synchronization, (4
) field output, (5
)
Fourier transforming, (6
) MPB, (7
) near to far field transformation, and (8
)
other.
def mean_time_spent_on(self, time_sink):
Return the mean time spent by all processes for a type of work time_sink
which
can be one of nine integer values 0
8
: (0
) connecting chunks, (1
) time
stepping, (2
) boundaries, (3
) MPI/synchronization, (4
) field output, (5
)
Fourier transforming, (6
) MPB, (7
) near to far field transformation, and (8
)
other.
def output_times(self, fname):
Call after running a simulation to output to a file with filename fname
the
times spent on various types of work as CSV (comma separated values) with headers
for each column and one row per process.
Field Computations
Meep supports a large number of functions to perform computations on the fields. Most of them are accessed via the lowerlevel C++/SWIG interface. Some of them are based on the following simpler, higherlevel versions. They are accessible as methods of a Simulation
instance.
def set_boundary(self, side, direction, condition):
Sets the condition of the boundary on the specified side in the specified
direction. See the Constants (Enumerated Types)
section for valid side
, direction
, and boundary_condition
values.
def phase_in_material(self, structure, time):
newstructure
should be the structure
field of another
Simulation
object with the same cell size and resolution.
Over the next time period phasetime
(in the current
simulation's time units), the current structure
(, , and conductivity ) will be
gradually changed to newstructure
. In particular, at each
timestep it linearly interpolates between the old structure
and the new structure. After phasetime
has elapsed, the
structure will remain equal to newstructure
. This is
demonstrated in the following image for two
Cylinder objects (the simulation script is in
examples/phase_in_material.py).
def get_field_point(self, c, pt):
Given a component
or derived_component
constant c
and a Vector3
pt
,
returns the value of that component at that point.
def get_epsilon_point(self, pt, frequency=0, omega=0):
Given a frequency frequency
and a Vector3
pt
, returns the average eigenvalue
of the permittivity tensor at that location and frequency. If frequency
is
nonzero, the result is complex valued; otherwise it is the real,
frequencyindependent part of (the limit).
def get_mu_point(self, pt, frequency=0, omega=0):
Given a frequency frequency
and a Vector3
pt
, returns the average eigenvalue
of the permeability tensor at that location and frequency. If frequency
is
nonzero, the result is complex valued; otherwise it is the real,
frequencyindependent part of (the limit).
def initialize_field(self, cmpnt, amp_func):
Initialize the component c
fields using the function func
which has a single
argument, a Vector3
giving a position and returns a complex number for the value
of the field at that point.
def add_dft_fields(self, *args, **kwargs):
def add_dft_fields(cs, fcen, df, nfreq, freq, where=None, center=None, size=None, yee_grid=False):
Given a list of field components cs
, compute the Fourier transform of these
fields for nfreq
equally spaced frequencies covering the frequency range
fcendf/2
to fcen+df/2
or an array/list freq
for arbitrarily spaced
frequencies over the Volume
specified by where
(default to the entire cell).
The volume can also be specified via the center
and size
arguments. The
default routine interpolates the Fourier transformed fields at the center of each
voxel within the specified volume. Alternatively, the exact Fourier transformed
fields evaluated at each corresponding Yee grid point is available by setting
yee_grid
to True
.
def flux_in_box(self, d, box=None, center=None, size=None):
Given a direction
constant, and a mp.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.
flux_in_box(d=mp.X,mp.Volume(center=mp.Vector3(0,0,0),size=mp.Vector3(0,1,1)))
If the center
and size
arguments are provided instead of box
, Meep will
construct the appropriate volume for you.
def electric_energy_in_box(self, box=None, center=None, size=None):
Given a mp.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. If the center
and
size
arguments are provided instead of box
, Meep will construct the
appropriate volume for you. Note: in cylindrical coordinates , the
integrand is
multiplied
by the circumference , or equivalently the integral is over an annular
volume.
def magnetic_energy_in_box(self, box=None, center=None, size=None):
Given a mp.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. If the center
and
size
arguments are provided instead of box
, Meep will construct the
appropriate volume for you. Note: in cylindrical coordinates , the
integrand is
multiplied
by the circumference , or equivalently the integral is over an annular
volume.
def field_energy_in_box(self, box=None, center=None, size=None):
Given a mp.Volume
, returns the integral of the electric and magneticfield
energy in the
given volume. If the volume has zero size along a dimension, a lowerdimensional
integral is used. If the center
and size
arguments are provided instead of
box
, Meep will construct the appropriate volume for you. Note: in cylindrical
coordinates , the integrand is
multiplied
by the circumference , or equivalently the integral is over an annular
volume.
def modal_volume_in_box(self, box=None, center=None, size=None):
Given a mp.Volume
, returns the instantaneous modal volume
according to the Purcelleffect definition:
.
If no volume argument is provided, the entire cell is used by
default. If the center
and size
arguments are provided
instead of box
, Meep will construct the appropriate volume
for you.
Note that if you are at a fixed frequency and you use complex fields (via
Blochperiodic 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.
Often, you want the integration box to be the entire cell. A useful function to
return this box, which you can then use for the box
arguments above, is
Simulation.total_volume()
.
One versatile 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 output_field_function
, described later). This is done
via the functions:
def integrate_field_function(self,
cs,
func,
where=None,
center=None,
size=None):
Returns the integral of the complexvalued function func
over the Volume
specified by where
(defaults to entire cell) for the meep::fields
contained in
the Simulation
instance that calls this method. 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 complexvalued. The
volume can optionally be specified via the center
and size
arguments.
If any dimension of where
is zero, that dimension is not integrated over. In
this way you can specify 1d, 2d, or 3d integrals.
Note: in cylindrical coordinates , the integrand is multiplied by the circumference , or equivalently the integral is over an annular volume.
def max_abs_field_function(self,
cs,
func,
where=None,
center=None,
size=None):
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. See Field
Functions for examples of how to call
integrate_field_function
and max_abs_field_function
. See Synchronizing the
Magnetic and Electric Fields
if you want to do computations combining the electric and magnetic fields. The
volume can optionally be specified via the center
and size
arguments.
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 cell, the same resolution, the same boundary conditions and symmetries (if any), and the same PML layers (if any).
def integrate2_field_function(self,
fields2,
cs1,
cs2,
func,
where=None,
center=None,
size=None):
Similar to integrate_field_function
, but takes additional parameters fields2
and cs2
. fields2
is a meep::fields*
object similar to the global fields
variable (see below) specifying the fields from another simulation. cs1
is a
list of components to integrate with from the meep::fields
instance in
Simulation.fields
, as for integrate_field_function
, while cs2
is a list of
components to integrate from fields2
. Similar to integrate_field_function
,
func
is a function that returns an number given arguments consisting of: the
position vector, followed by the values of the components specified by cs1
(in
order), followed by the values of the components specified by cs2
(in order).
The volume can optionally be specified via the center
and size
arguments.
To get two fields in memory at once for integrate2_field_function
, the easiest
way is to run one simulation within a given Python 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...
fields2 = sim.fields # save the fields in a variable
sim.fields = None # prevent the fields from getting deallocated by resetmeep
sim.reset_meep()
...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 fields_step()
.
Reloading Parameters
Once the fields/simulation have been initialized, you can change the values of various parameters by using the following functions (which are members of the Simulation
class):
def reset_meep(self):
Reset all of Meep's parameters, deleting the fields, structures, etcetera, from memory as if you had not run any computations.
def restart_fields(self):
Restart the fields at time zero, with zero fields. Does not reset the Fourier transforms of the flux planes, which continue to be accumulated.
def change_k_point(self, k):
Change the k_point
(the Bloch periodicity).
def change_sources(self, new_sources):
Change the list of sources in Simulation.sources
to new_sources
, and changes
the sources used for the current simulation. new_sources
must be a list of
Source
objects.
def set_materials(self, geometry=None, default_material=None):
This can be called in a step function, and is useful for changing the geometry or default material as a function of time.
Flux Spectra
Given a bunch of FluxRegion
objects, you can tell Meep to accumulate the Fourier transforms of the fields in those regions in order to compute the Poynting flux spectra. (Note: as a matter of convention, the "intensity" of the electromagnetic fields refers to the Poynting flux, not to the energy density.) See also Introduction/Transmittance/Reflectance Spectra and Tutorial/Basics/Transmittance Spectrum of a Waveguide Bend. These are attributes of the Simulation
class. The most important function is:
def add_flux(self, *args):
def add_flux(fcen, df, nfreq, freq, 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 for nfreq
equally spaced frequencies covering the
frequency range fcendf/2
to fcen+df/2
or an array/list freq
for arbitrarily
spaced frequencies. 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:
transmission = sim.add_flux(...)
to store the flux object in a variable. 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 timestepping by the run functions. At any time, you can ask for Meep to print out the current flux spectrum via the display_fluxes
method.
def display_fluxes(self, *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
or freq
. 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. display_fluxes(f1, f2, f3)
is actually equivalent to meep.display_csv("flux", meep.get_flux_freqs(f1), meep.get_fluxes(f1), meep.get_fluxes(f2), meep.get_fluxes(f3))
, where display_csv
takes a bunch of lists of numbers and prints them as a commaseparated table; this involves calling two lowerlevel functions:
def get_flux_freqs(f):
Given a flux object, returns a list of the frequencies that it is computing the spectrum for.
def get_fluxes(f):
Given a flux object, returns a list of the current flux spectrum that it has accumulated.
As described in Introduction/Transmittance/Reflectance Spectra and Tutorial/Basics/Transmittance Spectrum of a Waveguide Bend, 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:
def save_flux(self, fname, flux):
Save the Fouriertransformed fields corresponding to the given flux object in an
HDF5 file of the given filename
without the .h5
suffix (the current
filenameprefix is prepended automatically).
def load_flux(self, fname, flux):
Load the Fouriertransformed fields into the given flux object (replacing any
values currently there) from an HDF5 file of the given filename
without the
.h5
suffix (the current filenameprefix 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 cell and the flux regions) with the same number of
processors.
def load_minus_flux(self, fname, flux):
As load_flux
, 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.
Sometimes it is more convenient to keep the Fouriertransformed fields in memory rather than writing them to a file and immediately loading them back again. To that end, the Simulation
class exposes the following three methods:
def get_flux_data(self, flux):
Get the Fouriertransformed fields corresponding to the given flux object as a
FluxData
, which is just a named tuple of NumPy arrays. Note that this object is
only useful for passing to load_flux_data
below and should be considered opaque.
def load_flux_data(self, flux, fdata):
Load the Fouriertransformed fields into the given flux object (replacing any
values currently there) from the FluxData
object fdata
. You must load from an
object that was created by get_flux_data
in a simulation of the same dimensions
(for both the cell and the flux regions) with the same number of processors.
def load_minus_flux_data(self, flux, fdata):
As load_flux_data
, 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.
The Simulation
class also provides some aliases for the corresponding "flux" methods.
save_mode
load_mode
load_minus_mode
get_mode_data
load_mode_data
load_minus_mode_data
Mode Decomposition
Given a structure, Meep can decompose the Fouriertransformed fields into a superposition of its harmonic modes. For a theoretical background, see Features/Mode Decomposition.
def get_eigenmode_coefficients(self,
flux,
bands,
eig_parity=mp.NO_PARITY,
eig_vol=None,
eig_resolution=0,
eig_tolerance=1e12,
kpoint_func=None,
direction=mp.AUTOMATIC):
Given a flux object and list of band indices bands
or DiffractedPlanewave
, return a namedtuple
with the
following fields:
alpha
: the complex eigenmode coefficients as a 3d NumPy array of size (len(bands)
,flux.nfreqs
,2
). The last/third dimension refers to modes propagating in the forward (+) or backward () directions.vgrp
: the group velocity as a NumPy array.kpoints
: a list ofmp.Vector3
s of thekpoint
used in the mode calculation.kdom
: a list ofmp.Vector3
s of the mode's dominant wavevector.cscale
: a NumPy array of each mode's scaling coefficient. Useful for adjoint calculations.
The flux object should be created using add_mode_monitor
. (You could also use add_flux
, but with add_flux
you need to be more careful about symmetries that bisect the flux plane: the add_flux
object should only be used with get_eigenmode_coefficients
for modes of the same symmetry, e.g. constrained via eig_parity
. On the other hand, the performance of add_flux
planes benefits more from symmetry.) eig_vol
is the volume passed to MPB for the eigenmode calculation (based on interpolating the discretized materials from the Yee grid); in most cases this will simply be the volume over which the frequencydomain fields are tabulated, which is the default (i.e. flux.where
). eig_parity
should be one of [mp.NO_PARITY
(default), mp.EVEN_Z
, mp.ODD_Z
, mp.EVEN_Y
, mp.ODD_Y
]. It is the parity (= polarization in 2d) of the mode to calculate, assuming the structure has and/or mirror symmetry in the source region, just as for EigenModeSource
above. If the structure has both and mirror symmetry, you can combine more than one of these, e.g. EVEN_Z+ODD_Y
. Default is NO_PARITY
, 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. eig_resolution
is the spatial resolution to use in MPB for the eigenmode calculations. This defaults to twice the Meep resolution
in which case the structure is linearly interpolated from the Meep pixels. eig_tolerance
is the tolerance to use in the MPB eigensolver. MPB terminates when the eigenvalues stop changing to less than this fractional tolerance. Defaults to 1e12
. (Note that this is the tolerance for the frequency eigenvalue ω; the tolerance for the mode profile is effectively the square root of this.) For examples, see Tutorial/Mode Decomposition.
Technically, MPB computes ωₙ(k)
and then inverts it with Newton's method to find the wavevector k
normal to eig_vol
and mode for a given frequency; in rare cases (primarily waveguides with nonmonotonic dispersion relations, which doesn't usually happen in simple dielectric waveguides), MPB may need you to supply an initial "guess" for k
in order for this Newton iteration to converge. You can supply this initial guess with kpoint_func
, which is a function kpoint_func(f, n)
that supplies a rough initial guess for the k
of band number n
at frequency f = ω/2π
. (By default, the k components in the plane of the eig_vol
region are zero. However, if this region spans the entire cell in some directions, and the cell has Blochperiodic boundary conditions via the k_point
parameter, then the mode's k components in those directions will match k_point
so that the mode satisfies the Meep boundary conditions, regardless of kpoint_func
.) If direction
is set to mp.NO_DIRECTION
, then kpoint_func
is not only the initial guess and the search direction of the k vectors, but is also taken to be the direction of the waveguide, allowing you to detect modes in oblique waveguides (not perpendicular to the flux plane).
Note: for planewaves in homogeneous media, the kpoints
may not necessarily be equivalent to the actual wavevector of the mode. This quantity is given by kdom
.
def add_mode_monitor(self, *args, **kwargs):
def add_mode_monitor(fcen, df, nfreq, freq, ModeRegions...):
Similar to add_flux
, but for use with get_eigenmode_coefficients
.
add_mode_monitor
works properly with arbitrary symmetries, but may be suboptimal because the Fouriertransformed region does not exploit the symmetry. As an optimization, if you have a mirror plane that bisects the mode monitor, you can instead use add_flux
to gain a factor of two, but in that case you must also pass the corresponding eig_parity
to get_eigenmode_coefficients
in order to only compute eigenmodes with the corresponding mirror symmetry.
def get_eigenmode(self,
frequency,
direction,
where,
band_num,
kpoint,
eig_vol=None,
match_frequency=True,
parity=mp.NO_PARITY,
resolution=0,
eigensolver_tol=1e12):
The parameters of this routine are the same as that of
get_eigenmode_coefficients
or EigenModeSource
, but this function returns an
object that can be used to inspect the computed mode. In particular, it returns
an EigenmodeData
instance with the following fields:
band_num
: same as a single element of thebands
parameterfreq
: the computed frequency, same as thefrequency
input parameter ifmatch_frequency=True
group_velocity
: the group velocity of the mode indirection
k
: the Bloch wavevector of the mode indirection
kdom
: the dominant planewave of modeband_num
amplitude(point, component)
: the (complex) value of the given E or H fieldcomponent
(Ex
,Hy
, etcetera) at a particularpoint
(aVector3
) in space (interpreted with Blochperiodic boundary conditions if you give a point outside the originaleig_vol
).
If match_frequency=False
or kpoint
is not zero in the given direction
, the
frequency
input parameter is ignored.
The following toplevel function is also available:
def get_eigenmode_freqs(f):
Given a flux object, returns a list of the frequencies that it is computing the spectrum for.
DiffractedPlanewave
class DiffractedPlanewave(object):
For mode decomposition or eigenmode source, specify a diffracted planewave in homogeneous media. Should be passed as the bands
argument of get_eigenmode_coefficients
, band_num
of get_eigenmode
, or eig_band
of EigenModeSource
.
def __init__(self, g=None, axis=None, s=None, p=None):
Construct a DiffractedPlanewave
.

g
[ list of 3integer
s ] — The diffraction order corresponding to the wavevector . The diffraction order should be nonzero only in the 1 periodic directions of a dimensional cell (e.g., a plane in 3d) in which the mode monitor or source extends the entire length of the cell. 
axis
[Vector3
] — The plane of incidence for each planewave (used to define the and polarizations below) is defined to be the plane that contains theaxis
vector and the planewave's wavevector. IfNone
,axis
defaults to the first direction that lies in the plane of the monitor or source (e.g., direction for a plane in 3d, either or in 2d). 
s
[complex
] — The complex amplitude of the polarziation (i.e., electric field perpendicular to the plane of incidence). 
p
[complex
] — The complex amplitude of the polarziation (i.e., electric field parallel to the plane of incidence).
Energy Density Spectra
Very similar to flux spectra, you can also compute energy density spectra: the energy density of the electromagnetic fields as a function of frequency, computed by Fourier transforming the fields and integrating the energy density:
The usage is similar to the flux spectra: you define a set of EnergyRegion
objects telling Meep where it should compute the Fouriertransformed fields and energy densities, and call add_energy
to add these regions to the current simulation over a specified frequency bandwidth, and then use display_electric_energy
, display_magnetic_energy
, or display_total_energy
to display the energy density spectra at the end. There are also save_energy
, load_energy
, and load_minus_energy
functions that you can use to subtract the fields from two simulation, e.g. in order to compute just the energy from scattered fields, similar to the flux spectra. The function used to add an EnergyRegion
is as follows:
def add_energy(self, *args):
def add_energy(fcen, df, nfreq, freq, EnergyRegions...):
Add a bunch of EnergyRegion
s 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 fcendf/2
to fcen+df/2
or an array/list freq
for arbitrarily
spaced frequencies. Return an energy object, which you can pass to the functions
below to get the energy spectrum, etcetera.
As for flux regions, you normally use add_energy
via statements like:
En = sim.add_energy(...)
to store the energy object in a variable. You can create as many energy objects as you want, e.g. to look at the energy densities in different objects 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 energy region 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_energy
, 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 energy density spectrum via:
def display_electric_energy(self, *energys):
Given a number of energy objects, this displays a commaseparated table of
frequencies and energy density spectra for the electric fields prefixed by
"electric_energy1:" or similar (where the number is incremented after each run).
All of the energy should be for the same fcen
/df
/nfreq
or freq
. The first
column are the frequencies, and subsequent columns are the energy density spectra.
def display_magnetic_energy(self, *energys):
Given a number of energy objects, this displays a commaseparated table of
frequencies and energy density spectra for the magnetic fields prefixed by
"magnetic_energy1:" or similar (where the number is incremented after each run).
All of the energy should be for the same fcen
/df
/nfreq
or freq
. The first
column are the frequencies, and subsequent columns are the energy density spectra.
def display_total_energy(self, *energys):
Given a number of energy objects, this displays a commaseparated table of
frequencies and energy density spectra for the total fields "total_energy1:" or
similar (where the number is incremented after each run). All of the energy should
be for the same fcen
/df
/nfreq
or freq
. The first column are the
frequencies, and subsequent columns are the energy density spectra.
You might have to do something lowerlevel if you have multiple energy regions corresponding to different frequency ranges, or have other special needs. display_electric_energy(e1, e2, e3)
is actually equivalent to meep.display_csv("electric_energy", meep.get_energy_freqs(e1), meep.get_electric_energy(e1), meep.get_electric_energy(e2), meep.get_electric_energy(e3))
, where display_csv
takes a bunch of lists of numbers and prints them as a commaseparated table; this involves calling lowerlevel functions:
def get_energy_freqs(f):
Given an energy object, returns a list of the frequencies that it is computing the spectrum for.
def get_electric_energy(f):
Given an energy object, returns a list of the current energy density spectrum for the electric fields that it has accumulated.
def get_magnetic_energy(f):
Given an energy object, returns a list of the current energy density spectrum for the magnetic fields that it has accumulated.
def get_total_energy(f):
Given an energy object, returns a list of the current energy density spectrum for the total fields that it has accumulated.
As described in Introduction/Transmittance/Reflectance Spectra and Tutorial/Basics/Transmittance Spectrum of a Waveguide Bend for flux computations, to compute the energy density from the scattered fields 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:
def save_energy(self, fname, energy):
Save the Fouriertransformed fields corresponding to the given energy object in an
HDF5 file of the given filename
without the .h5
suffix (the current
filenameprefix is prepended automatically).
def load_energy(self, fname, energy):
Load the Fouriertransformed fields into the given energy object (replacing any
values currently there) from an HDF5 file of the given filename
without the
.h5
suffix (the current filenameprefix is prepended automatically). You must
load from a file that was saved by save_energy
in a simulation of the same
dimensions for both the cell and the energy regions with the same number of
processors.
def load_minus_energy(self, fname, energy):
As load_energy
, 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.
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 via . You should 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 Tutorial/Optical Forces.
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 by the incident power on an object, in order to compute the useful dimensionless ratio / where in Meep units. For example, it is a simple exercise to show that the force on a perfectly reflecting mirror with normalincident power satisfies /, and for a perfectly absorbing (black) surface /.
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 add_force
to add these regions to the current simulation over a specified frequency bandwidth, and then use display_forces
to display the force spectra at the end. There are also save_force
, load_force
, and load_minus_force
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. The function used to add a ForceRegion
object is defined as follows:
def add_force(self, *args):
def add_force(fcen, df, nfreq, freq, 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 for nfreq
equally spaced frequencies covering the
frequency range fcendf/2
to fcen+df/2
or an array/list freq
for arbitrarily
spaced frequencies. 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 add_force
via statements like:
Fx = sim.add_force(...)
to store the force object in a variable. 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 add_force
, 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:
def display_forces(self, *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
or freq
. 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. display_forces(f1, f2, f3)
is actually equivalent to meep.display_csv("force", meep.get_force_freqs(f1), meep.get_forces(f1), meep.get_forces(f2), meep.get_forces(f3))
, where display_csv
takes a bunch of lists of numbers and prints them as a commaseparated table; this involves calling two lowerlevel functions:
def get_force_freqs(f):
Given a force object, returns a list of the frequencies that it is computing the spectrum for.
def get_forces(f):
Given a force object, returns a list of the current force spectrum that it has accumulated.
As described in Introduction/Transmittance/Reflectance Spectra and Tutorial/Basics/Transmittance Spectrum of a Waveguide Bend for flux computations, to compute the force from the scattered fields often requires saving the Fouriertransformed fields from a "normalization" run and then loading them into another run to be subtracted. This can be done via these Simulation
methods:
def save_force(self, fname, force):
Save the Fouriertransformed fields corresponding to the given force object in an
HDF5 file of the given filename
without the .h5
suffix (the current
filenameprefix is prepended automatically).
def load_force(self, fname, force):
Load the Fouriertransformed fields into the given force object (replacing any
values currently there) from an HDF5 file of the given filename
without the
.h5
suffix (the current filenameprefix is prepended automatically). You must
load from a file that was saved by save_force
in a simulation of the same
dimensions for both the cell and the force regions with the same number of
processors.
def load_minus_force(self, fname, force):
As load_force
, 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.
To keep the fields in memory and avoid writing to and reading from a file, use the following three Simulation
methods:
def get_force_data(self, force):
Get the Fouriertransformed fields corresponding to the given force object as a
ForceData
, which is just a named tuple of NumPy arrays. Note that this object is
only useful for passing to load_force_data
below and should be considered
opaque.
def load_force_data(self, force, fdata):
Load the Fouriertransformed fields into the given force object (replacing any
values currently there) from the ForceData
object fdata
. You must load from an
object that was created by get_force_data
in a simulation of the same dimensions
(for both the cell and the flux regions) with the same number of processors.
def load_minus_force_data(self, force, fdata):
As load_force_data
, 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 Tutorial/Local Density of States. To do this, you simply pass the following step function to your run
command:
def Ldos(*args):
def Ldos(fcen, df, nfreq, freq):
Create an LDOS object with either frequency bandwidth df
centered at fcen
and
nfreq
equally spaced frequency points or an array/list freq
for arbitrarily spaced
frequencies. This can be passed to the dft_ldos
step function below as a keyword
argument.
def get_ldos_freqs(l):
Given an LDOS object, returns a list of the frequencies that it is computing the spectrum for.
def dft_ldos(*args, **kwargs):
def dft_ldos(fcen=None, df=None, nfreq=None, freq=None, ldos=None):
Compute the power spectrum of the sources (usually a single point dipole source),
normalized to correspond to the LDOS, in either a frequency bandwidth df
centered at
fcen
and nfreq
equally spaced frequency points or an array/list freq
for
arbitrarily spaced frequencies. One can also pass in an Ldos
object as
dft_ldos(ldos=my_Ldos)
.
The resulting spectrum is outputted as commadelimited text, prefixed by ldos:,
, and
is also stored in the ldos_data
variable of the Simulation
object after the run
is complete.
Analytically, the perpolarization LDOS is exactly proportional to the power radiated by an oriented pointdipole current, , at a given position in space. For a more mathematical treatment of the theory behind the LDOS, refer to the relevant discussion in Section 4.4 ("Currents and Fields: The Local Density of States") in Chapter 4 ("Electromagnetic Wave Source Conditions") of the book Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology, but for now it is defined as:
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 dft_ldos
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 can compute a neartofarfield transformation in the frequency domain as described in Tutorial/NeartoFar Field Spectra: given the fields on a "near" bounding surface inside the 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, 3d, or cylindrical region. That is, in a simulation surrounded by PML that absorbs outgoing waves, the neartofarfield feature can compute the fields outside the 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. For details, see Section 4.2.1 ("The Principle of Equivalence") in Chapter 4 ("Electromagnetic Wave Source Conditions") of the book Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology.
Note: in order for the farfield results to be accurate, the far region must be separated from the near region by at least 2D^{2}/λ, the Fraunhofer distance, where D is the largest dimension of the radiator and λ is the vacuum wavelength.
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 this Simulation
method:
def add_near2far(self, *args, **kwargs):
def add_near2far(fcen, df, nfreq, freq, Near2FarRegions..., nperiods=1):
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 for nfreq
equally spaced frequencies
covering the frequency range fcendf/2
to fcen+df/2
or an array/list freq
for arbitrarily spaced frequencies. 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 is complete, you can compute the far fields. This is usually for a pulsed source so that the fields have decayed away and the Fourier transforms have finished accumulating.
If you have Blochperiodic boundary conditions, then the corresponding neartofar transformation actually needs to perform a "lattice sum" of infinitely many periodic copies of the near fields. This doesn't happen by default, which means the default near2far
calculation may not be what you want for periodic boundary conditions. However, if the Near2FarRegion
spans the entire cell along the periodic directions, you can turn on an approximate lattice sum by passing nperiods > 1
. In particular, it then sums 2*nperiods+1
Blochperiodic copies of the near fields whenever a far field is requested. You can repeatedly double nperiods
until the answer converges to your satisfaction; in general, if the far field is at a distance d, and the period is a, then you want nperiods
to be much larger than d/a. (Future versions of Meep may use fancier techniques like Ewald summation to compute the lattice sum more rapidly at large distances.)
def get_farfield(self, near2far, x):
Given a Vector3
point x
which can lie anywhere outside the nearfield surface,
including outside the cell and a near2far
object, returns the computed
(Fouriertransformed) "far" fields at x
as list of length 6nfreq
, consisting
of fields
(E_{x}^{1},E_{y}^{1},E_{z}^{1},H_{x}^{1},H_{y}^{1},H_{z}^{1},E_{x}^{2},E_{y}^{2},E_{z}^{2},H_{x}^{2},H_{y}^{2},H_{z}^{2},...)
for the frequencies 1,2,…,nfreq
.
def output_farfields(self,
near2far,
fname,
resolution,
where=None,
center=None,
size=None):
Given an HDF5 file name fname
(does not include the .h5
suffix), a Volume
given by where
(may be 0d, 1d, 2d, or 3d), and a resolution
(in grid points /
distance unit), outputs the far fields in where
(which may lie outside the
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 datasets ex.r
, ex.i
,
..., hz.r
, hz.i
, giving the real and imaginary parts of the
Fouriertransformed and fields on this grid. Each dataset is an
nx×ny×nz×nfreq 4d array of space×frequency although dimensions
that =1 are omitted. The volume can optionally be specified via center
and
size
.
def get_farfields(self,
near2far,
resolution,
where=None,
center=None,
size=None):
Like output_farfields
but returns a dictionary of NumPy arrays instead of
writing to a file. The dictionary keys are Ex
, Ey
, Ez
, Hx
, Hy
, Hz
.
Each array has the same shape as described in output_farfields
.
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.
This lowerlevel function is also available:
def get_near2far_freqs(f):
Given a near2far
object, returns a list of the frequencies that it is computing the
spectrum for.
(Multifrequency get_farfields
and output_farfields
can be accelerated by
compiling Meep with withopenmp
and using the
OMP_NUM_THREADS
environment variable to specify multiple threads.)
For a scatteredfield computation, you often want to separate the scattered and incident fields. As described in Introduction/Transmittance/Reflectance Spectra and Tutorial/Basics/Transmittance Spectrum of a Waveguide Bend 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 these Simulation
methods:
def save_near2far(self, fname, near2far):
Save the Fouriertransformed fields corresponding to the given near2far
object
in an HDF5 file of the given filename
(without the .h5
suffix). The current
filenameprefix is prepended automatically.
def load_near2far(self, fname, near2far):
Load the Fouriertransformed fields into the given near2far
object (replacing
any values currently there) from an HDF5 file of the given filename
without the
.h5
suffix (the current filenameprefix is prepended automatically). You must
load from a file that was saved by save_near2far
in a simulation of the same
dimensions for both the cell and the near2far regions with the same number of
processors.
def load_minus_near2far(self, fname, near2far):
As load_near2far
, 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.
To keep the fields in memory and avoid writing to and reading from a file, use the following three methods:
def get_near2far_data(self, near2far):
Get the Fouriertransformed fields corresponding to the given near2far
object as
a NearToFarData
, which is just a named tuple of NumPy arrays. Note that this
object is only useful for passing to load_near2far_data
below and should be
considered opaque.
def load_near2far_data(self, near2far, n2fdata):
Load the Fouriertransformed fields into the near2far
object (replacing any
values currently there) from the NearToFarData
object n2fdata
. You must load
from an object that was created by get_near2far_data
in a simulation of the same
dimensions (for both the cell and the flux regions) with the same number of
processors.
def load_minus_near2far_data(self, near2far, n2fdata):
As load_near2far_data
, 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.
See also this lowerlevel function:
def scale_near2far_fields(s, near2far):
Scale the Fouriertransformed fields in near2far
by the complex number s
. e.g.
load_minus_near2far
is equivalent to load_near2far
followed by
scale_near2far_fields
with s=1
.
And this DftNear2Far
method:
def flux(self, direction, where, resolution):
Given a Volume
where
(may be 0d, 1d, 2d, or 3d) and a resolution
(in grid
points / distance unit), compute the far fields in where
(which may lie
outside the cell) in a grid with the given resolution (which may differ from the
FDTD solution) and return its Poynting flux in direction
as a list. The dataset
is a 1d array of nfreq
dimensions.
Load and Dump Structure
These functions dump the raw ε and μ data to disk and load it back for doing multiple simulations with the same materials but different sources etc. The only prerequisite is that the dump/load simulations have the same chunks (i.e. the same grid, number of processors, symmetries, and PML). When using split_chunks_evenly=False
, you must also dump the original chunk layout using dump_chunk_layout
and load it into the new Simulation
using the chunk_layout
parameter. Currently only stores dispersive and nondispersive and but not nonlinearities. Note that loading data from a file in this way overwrites any geometry
data passed to the Simulation
constructor.
def dump_structure(self, fname):
Dumps the structure to the file fname
.
def load_structure(self, fname):
Loads a structure from the file fname
. A file name to load can also be passed to
the Simulation
constructor via the load_structure
keyword argument.
def dump_chunk_layout(self, fname):
Dumps the chunk layout to file fname
.
To load a chunk layout into a Simulation
, use the chunk_layout
argument to the constructor, passing either a file obtained from dump_chunk_layout
or another Simulation
instance. Note that when using split_chunks_evenly=False
this parameter is required when saving and loading flux spectra, force spectra, or neartofar spectra so that the two runs have the same chunk layout. Just pass the Simulation
object from the first run to the second run:
# Split chunks based on amount of work instead of size
sim1 = mp.Simulation(..., split_chunks_evenly=False)
norm_flux = sim1.add_flux(...)
sim1.run(...)
sim1.save_flux(...)
# Make sure the second run uses the same chunk layout as the first
sim2 = mp.Simulation(..., chunk_layout=sim1)
flux = sim2.add_flux(...)
sim2.load_minus_flux(...)
sim2.run(...)
FrequencyDomain Solver
Meep contains a frequencydomain solver that computes the fields produced in a geometry in response to a continuouswave (CW) source. This is based on an iterative linear solver instead of timestepping. For details, see Section 5.3 ("Frequencydomain solver") of Computer Physics Communications, Vol. 181, pp. 687702, 2010. Benchmarking results have shown that in many instances, such as cavities (e.g., ring resonators) with longlived resonant modes, this solver converges much faster than simply running an equivalent timedomain simulation with a CW source (using the default width
of zero for no transient turnon), timestepping until all transient effects from the source turnon have disappeared, especially if the fields are desired to a very high accuracy.
To use the frequencydomain solver, simply define a ContinuousSource
with the desired frequency and initialize the fields and geometry via init_sim()
:
sim = mp.Simulation(...)
sim.init_sim()
sim.solve_cw(tol, maxiters, L)
The first two parameters to the frequencydomain 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 that determines a tradeoff between memory and work per step and convergence rate of the iterative algorithm, biconjugate gradient stabilized (BiCGSTABL), that is used; larger values of will often lead to faster convergence at the expense of more memory and more work per iteration. Default is , and normally a value ≥ 2 should be used.
The frequencydomain solver supports arbitrary geometries, PML, boundary conditions, symmetries, parallelism, conductors, and arbitrary nondispersive materials. LorentzDrude 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 force_complex_fields=True
.
After solve_cw
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. For examples, see Tutorial/Frequency Domain Solver and Tutorial/Mode Decomposition/Reflectance and Transmittance Spectra for Planewave at Oblique Incidence.
Note: The convergence of the iterative solver can sometimes encounter difficulties. For example, increasing the diameter of a ring resonator relative to the wavelength increases the condition number, which worsens the convergence of iterative solvers. The general way to improve this is to implement a more sophisticated iterative solver that employs preconditioners. Preconditioning wave equations (Helmholtzlike equations) is notoriously difficult to do well, but some possible strategies are discussed in Issue #548. In the meantime, a simpler way improving convergence (at the expense of computational cost) is to increase the parameter and the number of iterations.
FrequencyDomain Eigensolver
Building on the frequencydomain solver above, Meep also includes a frequencydomain eigensolver that computes resonant frequencies and modes in the frequency domain. The usage is very similar to solve_cw
:
sim = mp.Simulation(...)
sim.init_sim()
eigfreq = sim.solve_eigfreq(tol, maxiters, guessfreq, cwtol, cwmaxiters, L)
The solve_eig
routine performs repeated calls to solve_cw
in a way that converges to the resonant mode whose frequency is closest to the source frequency. The complex resonantmode frequency is returned, and the mode Q can be computed from eigfreq.real / (2*eigfreq.imag)
. Upon return, the fields should be the corresponding resonant mode (with an arbitrary scaling).
The resonant mode is converged to a relative error of roughly tol
, which defaults to 1e7
. A maximum of maxiters
(defaults to 100
) calls to solve_cw
are performed. The tolerance for each solve_cw
call is cwtol
(defaults to tol*1e3
) and the maximum iterations is cwmaxiters
(10^{4}, by default); the L
parameter (defaults to 10
) is also passed through to solve_cw
.
The closer the input frequency is to the resonantmode frequency, the faster solve_eig
should converge. Instead of using the source frequency, you can instead pass a guessfreq
argument to solve_eigfreq
specifying an input frequency (which may even be complex).
Technically, solve_eig
is using a shiftandinvert power iteration to compute the resonant mode, as reviewed in FrequencyDomain Eigensolver.
As for solve_cw
above, you are required to set force_complex_fields=True
to use solve_eigfreq
.
GDSII Support
This feature is only available if Meep is built with libGDSII. It so, then the following functions are available:
def GDSII_layers(fname):
Returns a list of integervalued layer indices for the layers present in the specified GDSII file.
mp.GDSII_layers('python/examples/coupler.gds')
Out[2]: [0, 1, 2, 3, 4, 5, 31, 32]
def GDSII_prisms(material, fname, layer=1, zmin=0.0, zmax=0.0):
Returns a list of GeometricObject
s with material
(mp.Medium
) on layer number
layer
of a GDSII file fname
with zmin
and zmax
(default 0).
def GDSII_vol(fname, layer, zmin, zmax):
Returns a mp.Volume
read from a GDSII file fname
on layer number layer
with
zmin
and zmax
(default 0). This function is useful for creating a FluxRegion
from a GDSII file as follows:
fr = mp.FluxRegion(volume=mp.GDSII_vol(fname, layer, zmin, zmax))
Data Visualization
This module provides basic visualization functionality for the simulation domain. The spirit of the module is to provide functions that can be called with no customization options whatsoever and will do useful relevant things by default, but which can also be customized in cases where you do want to take the time to spruce up the output. The Simulation
class provides the following methods:
def plot2D(self,
ax=None,
output_plane=None,
fields=None,
labels=False,
eps_parameters=None,
boundary_parameters=None,
source_parameters=None,
monitor_parameters=None,
field_parameters=None,
frequency=None,
plot_eps_flag=True,
plot_sources_flag=True,
plot_monitors_flag=True,
plot_boundaries_flag=True,
**kwargs):
Plots a 2D cross section of the simulation domain using matplotlib
. The plot
includes the geometry, boundary layers, sources, and monitors. Fields can also be
superimposed on a 2D slice. Requires matplotlib. Calling
this function would look something like:
sim = mp.Simulation(...)
sim.run(...)
field_func = lambda x: 20*np.log10(np.abs(x))
import matplotlib.pyplot as plt
sim.plot2D(fields=mp.Ez,
field_parameters={'alpha':0.8, 'cmap':'RdBu', 'interpolation':'none', 'post_process':field_func},
boundary_parameters={'hatch':'o', 'linewidth':1.5, 'facecolor':'y', 'edgecolor':'b', 'alpha':0.3})
plt.show()
plt.savefig('sim_domain.png')
Note: When running a parallel simulation, the plot2D
function expects to be called
on all processes, but only generates a plot on the master process.
Parameters:
ax
: amatplotlib
axis object.plot2D()
will add plot objects, like lines, patches, and scatter plots, to this object. If noax
is supplied, then the routine will create a new figure and grab its axis.output_plane
: aVolume
object that specifies the plane over which to plot. Must be 2D and a subset of the grid volume (i.e., it should not extend beyond the cell).fields
: the field component (mp.Ex
,mp.Ey
,mp.Ez
,mp.Hx
,mp.Hy
,mp.Hz
) to superimpose over the simulation geometry. Default isNone
, where no fields are superimposed.labels
: ifTrue
, then labels will appear over each of the simulation elements.eps_parameters
: adict
of optional plotting parameters that override the default parameters for the geometry.interpolation='spline36'
: interpolation algorithm used to upsample the pixels.cmap='binary'
: the color map of the geometryalpha=1.0
: transparency of geometrycontour=False
: ifTrue
, plot a contour of the geometry rather than its imagecontour_linewidth=1
: line width of the contour lines ifcontour=True
boundary_parameters
: adict
of optional plotting parameters that override the default parameters for the boundary layers.alpha=1.0
: transparency of boundary layersfacecolor='g'
: color of polygon faceedgecolor='g'
: color of outline strokelinewidth=1
: line width of outline strokehatch='\'
: hatching pattern
source_parameters
: adict
of optional plotting parameters that override the default parameters for the sources.color='r'
: color of line and pt sourcesalpha=1.0
: transparency of sourcefacecolor='none'
: color of polygon face for planar sourcesedgecolor='r'
: color of outline stroke for planar sourceslinewidth=1
: line width of outline strokehatch='\'
: hatching patternlabel_color='r'
: color of source labelslabel_alpha=0.3
: transparency of source label boxoffset=20
: distance from source center and label box
monitor_parameters
: adict
of optional plotting parameters that override the default parameters for the monitors.color='g'
: color of line and point monitorsalpha=1.0
: transparency of monitorsfacecolor='none'
: color of polygon face for planar monitorsedgecolor='r'
: color of outline stroke for planar monitorslinewidth=1
: line width of outline strokehatch='\'
: hatching patternlabel_color='g'
: color of source labelslabel_alpha=0.3
: transparency of monitor label boxoffset=20
: distance from monitor center and label box
field_parameters
: adict
of optional plotting parameters that override the default parameters for the fields.interpolation='spline36'
: interpolation function used to upsample field pixelscmap='RdBu'
: color map for field pixelsalpha=0.6
: transparency of fieldspost_process=np.real
: post processing function to apply to fields (must be a function object)
frequency
: for materials with a frequencydependent permittivity , specifies the frequency (in Meep units) of the real part of the permittivity to use in the plot. Defaults to thefrequency
parameter of the Source object.
def plot3D(self):
Uses Mayavi to render a 3D simulation domain. The simulation object must be 3D. Can also be embedded in Jupyter notebooks.
def visualize_chunks(self):
Displays an interactive image of how the cell is divided into chunks. Each rectangular region is a chunk, and each color represents a different processor. Requires matplotlib.
An animated visualization is also possible via the Animate2D class.
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. These are attributes of the Simulation
class.
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.
A common point of confusion is described in The Run Function Is Not A Loop. 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.
def run(self, *step_funcs, **kwargs):
def run(step_functions..., until=condition/time):
def run(step_functions..., until_after_sources=condition/time):
run(step_functions..., until=condition/time)
Run the simulation until a certain time or condition, calling the given step
functions (if any) at each timestep. The keyword argument until
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. until
can also be a list of stopping conditions which may include a number
and additional functions.
run(step_functions..., until_after_sources=condition/time)
Run the simulation until all sources have turned off, calling the given step
functions (if any) at each timestep. The keyword argument until_after_sources
is
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 condition
returns True
. Like until
above, until_after_sources
can take a list of
stopping conditions.
In particular, a useful value for until_after_sources
or until
is often stop_when_field_decayed
, which is demonstrated in Tutorial/Basics. These toplevel functions are available:
def stop_when_fields_decayed(dt, c, pt, decay_by):
Return a condition
function, suitable for passing to Simulation.run
as the until
or until_after_sources
parameter, 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 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.
def stop_after_walltime(t):
Return a condition
function, suitable for passing to Simulation.run
as the until
parameter. Stops the simulation after t
seconds of wall time have passed.
def stop_on_interrupt():
Return a condition
function, suitable for passing to Simulation.run
as the until
parameter. Instead of terminating when receiving a SIGINT or SIGTERM signal from the
system, the simulation will abort time stepping and continue executing any code that
follows the run
function (e.g., outputting fields).
Finally, another run function, useful for computing ω(k) band diagrams, is available via these Simulation
methods:
def run_k_points(self, t, k_points):
Given a list of Vector3
, k_points
of k vectors, 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 additional time units. It will run
Harminv at the same point/component as the first Gaussian source and
look for modes in the union of the frequency ranges for all sources. Returns a
list of lists of frequencies (one list of frequencies for each k). Also prints
out a commadelimited list of frequencies, prefixed by freqs:
, and their
imaginary parts, prefixed by freqsim:
. See Tutorial/Resonant Modes and
Transmission in a Waveguide
Cavity.
def run_k_point(self, t, k):
Lower level function called by run_k_points
that runs a simulation for a single
k point k_point
and returns a Harminv
instance. Useful when you need to
access more Harminv
data than just the frequencies.
Predefined Step Functions
Several useful step functions are predefined by Meep. These are available directly via the meep
package but require a Simulation
instance as an argument.
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.
def output_dft(self, dft_fields, fname):
Output the Fouriertransformed fields in dft_fields
(created by
add_dft_fields
) to an HDF5 file with name fname
(does not include the .h5
suffix).
def output_epsilon(sim=None, *step_func_args, **kwargs):
Given a frequency frequency
, (provided as a keyword argument) output (relative
permittivity); for an anisotropic tensor the output is the harmonic
mean of the eigenvalues. If
frequency
is nonzero, the output is complex; otherwise it is the real,
frequencyindependent part of (the limit).
When called as part of a step function,
the sim
argument specifying the Simulation
object can be omitted, e.g.,
sim.run(mp.at_beginning(mp.output_epsilon(frequency=1/0.7)),until=10)
.
def output_mu(sim=None, *step_func_args, **kwargs):
Given a frequency frequency
, (provided as a keyword argument) output (relative
permeability); for an anisotropic tensor the output is the harmonic
mean of the eigenvalues. If
frequency
is nonzero, the output is complex; otherwise it is the real,
frequencyindependent part of (the limit).
When called as part of a step function,
the sim
argument specifying the Simulation
object can be omitted, e.g.,
sim.run(mp.at_beginning(mp.output_mu(frequency=1/0.7)),until=10)
.
def output_poynting(sim):
Output the Poynting flux . Note that you
might want to wrap this step function in synchronized_magnetic
to compute it more
accurately. See Synchronizing the Magnetic and Electric
Fields.
def output_hpwr(sim):
Output the magneticfield energy density
def output_dpwr(sim):
Output the electricfield energy density
def output_tot_pwr(sim):
Output the total electric and magnetic energy density. Note that you might want to
wrap this step function in synchronized_magnetic
to compute it more accurately. See
Synchronizing the Magnetic and Electric
Fields.
def output_png(compnt, options, rm_h5=True):
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 Tutorial/Basics/Output Tips and
Tricks.
It is often useful to use the h5topng C
or A
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
output_epsilon
. To make this easier, a builtin shell variable $EPS
is provided
which refers to the lastoutput dielectricfunction .h5
file. So, for example
output_png(mp.Ez,"C $EPS")
will output the field and overlay the dielectric
contours.
By default, output_png
deletes the .h5
file when it is done. To preserve the .h5
file requires output_png(component, h5topng_options, rm_h5=False)
.
def output_hfield(sim):
Outputs all the components of the field h, (magnetic) to an HDF5 file. That is, the different components are stored as different datasets within the same file.
def output_hfield_x(sim):
Output the component of the field h (magnetic). 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.
def output_hfield_y(sim):
Output the component of the field h (magnetic). 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.
def output_hfield_z(sim):
Output the component of the field h (magnetic). 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.
def output_hfield_r(sim):
Output the component of the field h (magnetic). 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.
def output_hfield_p(sim):
Output the component of the field h (magnetic). 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.
def output_bfield(sim):
Outputs all the components of the field b, (magnetic) to an HDF5 file. That is, the different components are stored as different datasets within the same file.
def output_bfield_x(sim):
Output the component of the field b (magnetic). 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.
def output_bfield_y(sim):
Output the component of the field b (magnetic). 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.
def output_bfield_z(sim):
Output the component of the field b (magnetic). 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.
def output_bfield_r(sim):
Output the component of the field b (magnetic). 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.
def output_bfield_p(sim):
Output the component of the field b (magnetic). 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. Note that for outputting the Poynting flux, you
might want to wrap the step function in synchronized_magnetic
to compute it more
accurately. See Synchronizing the Magnetic and Electric
Fields.
def output_efield(sim):
Outputs all the components of the field e, (electric) to an HDF5 file. That is, the different components are stored as different datasets within the same file.
def output_efield_x(sim):
Output the component of the field e (electric). 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.
def output_efield_y(sim):
Output the component of the field e (electric). 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.
def output_efield_z(sim):
Output the component of the field e (electric). 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.
def output_efield_r(sim):
Output the component of the field e (electric). 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.
def output_efield_p(sim):
Output the component of the field e (electric). 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. Note that for outputting the Poynting flux, you
might want to wrap the step function in synchronized_magnetic
to compute it more
accurately. See Synchronizing the Magnetic and Electric
Fields.
def output_dfield(sim):
Outputs all the components of the field d, (displacement) to an HDF5 file. That is, the different components are stored as different datasets within the same file.
def output_dfield_x(sim):
Output the component of the field d (displacement). 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.
def output_dfield_y(sim):
Output the component of the field d (displacement). 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.
def output_dfield_z(sim):
Output the component of the field d (displacement). 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.
def output_dfield_r(sim):
Output the component of the field d (displacement). 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.
def output_dfield_p(sim):
Output the component of the field d (displacement). 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. Note that for outputting the Poynting flux, you
might want to wrap the step function in synchronized_magnetic
to compute it more
accurately. See Synchronizing the Magnetic and Electric
Fields.
def output_sfield(sim):
Outputs all the components of the field s, (poynting flux) to an HDF5 file. That
is, the different components are stored as different datasets within the same file.
Note that you might want to wrap this step function in synchronized_magnetic
to
compute it more accurately. See Synchronizing the Magnetic and Electric
Fields.
def output_sfield_x(sim):
Output the component of the field s (poynting flux). 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.
def output_sfield_y(sim):
Output the component of the field s (poynting flux). 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.
def output_sfield_z(sim):
Output the component of the field s (poynting flux). 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.
def output_sfield_r(sim):
Output the component of the field s (poynting flux). 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.
def output_sfield_p(sim):
Output the component of the field s (poynting flux). 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. Note that for outputting the Poynting flux, you
might want to wrap the step function in synchronized_magnetic
to compute it more
accurately. See Synchronizing the Magnetic and Electric
Fields.
More generally, it is possible to output an arbitrary function of position and zero or more field components, similar to the Simulation.integrate_field_function
method, described above. This is done by:
def output_field_function(self,
name,
cs,
func,
real_only=False,
h5file=None):
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. If
real_only
is True, only outputs the real part of func
.
See also Field Functions, and Synchronizing the Magnetic and Electric Fields if you want to do computations combining the electric and magnetic fields.
Array Slices
The output functions described above write the data for the fields and materials for the entire cell to an HDF5 file. This is useful for postprocessing large datasets which may not fit into memory as you can later read in the HDF5 file to obtain field/material data as a NumPy array. However, in some cases it is convenient to bypass the disk altogether to obtain the data directly in the form of a NumPy array without writing/reading HDF5 files. Additionally, you may want the field/material data on just a subregion (or slice) of the entire volume. This functionality is provided by the get_array
method which takes as input a subregion of the cell and the field/material component. The method returns a NumPy array containing values of the field/material at the current simulation time.
def get_array(self,
component=None,
vol=None,
center=None,
size=None,
cmplx=None,
arr=None,
frequency=0,
snap=False):
Takes as input a subregion of the cell and the field/material component. The method returns a NumPy array containing values of the field/material at the current simulation time.
Parameters:

vol
:Volume
; the orthogonal subregion/slice of the computational volume. The return value ofget_array
has the same dimensions as theVolume
'ssize
attribute. IfNone
(default), then asize
andcenter
must be specified. 
center
,size
:Vector3
; if both are specified, the library will construct an appropriateVolume
. This is a convenience feature and alternative to supplying aVolume
. 
component
: field/material component (i.e.,mp.Ex
,mp.Hy
,mp.Sz
,mp.Dielectric
, etc). Defaults toNone
. 
cmplx
:boolean
; ifTrue
, return complexvalued data otherwise return realvalued data (default). 
arr
: optional field to pass a preallocated NumPy array of the correct size, which will be overwritten with the field/material data instead of allocating a new array. Normally, this will be the array returned from a previous call toget_array
for a similar slice, allowing one to reusearr
(e.g., when fetching the same slice repeatedly at different times). 
frequency
: optional frequency point over which the average eigenvalue of the and tensors are evaluated (defaults to 0). 
snap
: By default, the elements of the grid slice are obtained using a bilinear interpolation of the nearest Yee grid points. Empty dimensions of the grid slice are "collapsed" into a single element. However, ifsnap
is set toTrue
, this interpolation behavior is disabled and the grid slice is instead "snapped" everywhere to the nearest grid point. (Empty slice dimensions are still of size one.) This feature is mainly useful for comparing results with theoutput_
routines (e.g.,output_epsilon
,output_efield_z
, etc.).
For convenience, the following wrappers for get_array
over the entire cell are
available: get_epsilon()
, get_mu()
, get_hpwr()
, get_dpwr()
,
get_tot_pwr()
, get_Xfield()
, get_Xfield_x()
, get_Xfield_y()
,
get_Xfield_z()
, get_Xfield_r()
, get_Xfield_p()
where X
is one of h
, b
,
e
, d
, or s
. The routines get_Xfield_*
all return an array type consistent
with the fields (real or complex). The routines get_epsilon()
and get_mu()
accept the optional argument frequency
(defaults to 0) and all routines accept
snap
(defaults to False
).
Note on arrayslice dimensions: The routines get_epsilon
, get_Xfield_z
,
etc. use as default size=meep.Simulation.fields.total_volume()
which for
simulations involving Blochperiodic boundaries (via k_point
) will result in
arrays that have slightly different dimensions than e.g.
get_array(center=meep.Vector3(), size=cell_size, component=meep.Dielectric
, etc.
(i.e., the slice spans the entire cell volume cell_size
). Neither of these
approaches is "wrong", they are just slightly different methods of fetching the
boundaries. The key point is that if you pass the same value for the size
parameter, or use the default, the slicing routines always give you the samesize
array for all components. You should not try to predict the exact size of these
arrays; rather, you should simply rely on Meep's output.
def get_dft_array(self, dft_obj, component, num_freq):
Returns the Fouriertransformed fields as a NumPy array.
Parameters:

dft_obj
: adft_flux
,dft_force
,dft_fields
, ordft_near2far
object obtained from calling the appropriateadd
function (e.g.,mp.add_flux
). 
component
: a field component (e.g.,mp.Ez
). 
num_freq
: the index of the frequency. An integer in the range0...nfreq1
, wherenfreq
is the number of frequencies stored indft_obj
as set by thenfreq
parameter toadd_dft_fields
,add_dft_flux
, etc.
Array Metadata
def get_array_metadata(self,
vol=None,
center=None,
size=None,
dft_cell=None,
return_pw=False):
This routine provides geometric information useful for interpreting the arrays
returned by get_array
or get_dft_array
for the spatial region defined by vol
or center
/size
. In both cases, the return value is a tuple (x,y,z,w)
, where:
x,y,z
are 1d NumPy arrays storing the coordinates of the points in the grid slicew
is a NumPy array of the same dimensions as the array returned byget_array
/get_dft_array
, whose entries are the weights in a cubature rule for integrating over the spatial region (with the points in the cubature rule being just the grid points contained in the region). Thus, if is some spatiallyvarying quantity whose value at the th grid point is , the integral of over the region may be approximated by the sum:
This is a 1, 2, or 3dimensional integral depending on the number of dimensions
in which has zero extent. If the samples are stored in an
array Q
of the same dimensions as w
, then evaluating the sum on the RHS is
just one line: np.sum(w*Q).
A convenience parameter dft_cell
is provided as an alternative to vol
or
center
/size
. Set dft_cell
to a dft_flux
or dft_fields
object to define the
region covered by the array. If the dft_cell
argument is provided then all other
arguments related to the spatial region (vol
, center
, and size
) are ignored.
If no arguments are provided, then the entire cell is used.
For empty dimensions of the grid slice get_array_metadata
will collapse
the two elements corresponding to the nearest Yee grid points into a single
element using linear interpolation.
If return_pw=True
, the return value is a 2tuple (p,w)
where p
(points) is a
list of mp.Vector3
s with the same dimensions as w
(weights). Otherwise, by
default the return value is a 4tuple (x,y,z,w)
.
The following are some examples of how array metadata can be used.
Labeling Axes in Plots of Grid Quantities
# using the geometry from the bendflux tutorial example
import matplotlib.pyplot as plt
import numpy as np
eps_array=sim.get_epsilon()
(x,y,z,w)=sim.get_array_metadata()
plt.figure()
ax = plt.subplot(111)
plt.pcolormesh(x,y,np.transpose(eps_array),shading='gouraud')
ax.set_aspect('equal')
plt.show()
Computing Quantities Defined by Integrals of FieldDependent Functions Over Grid Regions

energy stored in the field in a region :

Poynting flux through a surface :
import numpy as np
# Efield modal volume in box from timedomain fields
box = mp.Volume(center=box_center, size=box_size)
(Ex,Ey,Ez) = [sim.get_array(vol=box, component=c, cmplx=True) for c in [mp.Ex, mp.Ey, mp.Ez]]
eps = sim.get_array(vol=box, component=mp.Dielectric)
(x,y,z,w) = sim.get_array_metadata(vol=box)
energy_density = np.real(eps*(np.conj(Ex)*Ex + np.conj(Ey)*Ey + np.conj(Ez)*Ez)) # array
energy = np.sum(w*energy_density) # scalar
# xdirected Poynting flux through monitor from frequencydomain fields
monitor = mp.FluxRegion(center=mon_center, size=mon_size)
dft_cell = sim.add_flux(freq, 0, 1, monitor)
sim.run(...) # timestep until DFTs converged
(Ey,Ez,Hy,Hz) = [sim.get_dft_array(dft_cell,c,0) for c in [mp.Ey, mp.Ez, mp.Hy, mp.Hz]]
(x,y,z,w) = sim.get_array_metadata(dft=dft_cell)
flux_density = np.real( np.conj(Ey)*Hz  np.conj(Ez)*Hy ) # array
flux = np.sum(w*flux_density) # scalar
Source Slices
def get_source(self,
component,
vol=None,
center=None,
size=None):
Return an array of complex values of the source amplitude for
component
over the given vol
or center
/size
. The array has the same
dimensions as that returned by get_array
.
Not supported for cylindrical coordinates.
Harminv Step Function
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.
StepFunction Modifiers
Rather than writing a brandnew step function every time something a bit different is required, the following "modifier" functions take a bunch of step functions and produce new step functions with modified behavior. See also Tutorial/Basics for examples.
Miscellaneous StepFunction Modifiers
def combine_step_funcs(*step_funcs):
Given zero or more step functions, return a new step function that on each step calls all of the passed step functions.
def synchronized_magnetic(*step_funcs):
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
def when_true(cond, *step_funcs):
Given zero or more step functions and a condition function condition
(a function of
no arguments), evaluate the step functions whenever condition
returns True
.
def when_false(cond, *step_funcs):
Given zero or more step functions and a condition function condition
(a function of
no arguments), evaluate the step functions whenever condition
returns False
.
def at_every(dt, *step_funcs):
Given zero or more step functions, evaluates them at every time interval of units (rounded up to the next time step).
def after_time(t, *step_funcs):
Given zero or more step functions, evaluates them only for times after a time units have elapsed from the start of the run.
def before_time(t, *step_funcs):
Given zero or more step functions, evaluates them only for times before a time units have elapsed from the start of the run.
def at_time(t, *step_funcs):
Given zero or more step functions, evaluates them only once, after a time units have elapsed from the start of the run.
def after_sources(*step_funcs):
Given zero or more step functions, evaluates them only for times after all of the sources have turned off.
def after_sources_and_time(t, *step_funcs):
Given zero or more step functions, evaluates them only for times after all of the sources have turned off, plus an additional time units have elapsed.
def during_sources(*step_funcs):
Given zero or more step functions, evaluates them only for times before all of the sources have turned off.
def at_beginning(*step_funcs):
Given zero or more step functions, evaluates them only once, at the beginning of the run.
def at_end(*step_funcs):
Given zero or more step functions, evaluates them only once, at the end of the run.
Modifying HDF5 Output
def in_volume(v, *step_funcs):
Given zero or more step functions, modifies any output functions among them to only
output a subset (or a superset) of the cell, corresponding to the meep::volume* v
(created by the Volume
function).
def in_point(pt, *step_funcs):
Given zero or more step functions, modifies any output functions among them to only
output a single point of data, at pt
(a Vector3
).
def to_appended(fname, *step_funcs):
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.
def with_prefix(pre, *step_funcs):
Given zero or more step functions, modifies any output functions among them to prepend
the string prefix
to the file names (much like filename_prefix
, above).
Writing Your Own Step Functions
A step function can take two forms. The simplest is just a function with one argument (the simulation instance), which is called at every time step unless modified by one of the modifier functions above. e.g.
def my_step(sim):
print("Hello world!")
If one then does sim.run(my_step, until=100)
, 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 two arguments: the second argument will either be step
when it is called during timestepping, or finish
when it is called at the end of the run:
def my_step(sim, todo):
if todo == 'step':
# do something
elif todo == 'finish':
# do something else
# access simulation attributes
sim.fields ...etc.
LowLevel Functions
By default, Meep initializes C++ objects like meep::structure
and meep::fields
in the Simulation
object based on attributes like sources
and geometry
. Theses objects are then accessible via simulation_instance.structure
and simulation_instance.fields
. Given these, you can then call essentially any function in the C++ interface, because all of the C++ functions are automatically made accessible to Python by the wrappergenerator program 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 Simulation.init_sim()
manually, or Simulation._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:
sim = Simulation(cell_size, resolution).init_sim()
my_fields = sim.fields
sim.fields = None
sim.reset_meep()
and then change the geometry etc. and rerun sim.init_sim()
. 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 Python function.
 First, all functions in the
meep::
namespace are available in the Meep Python module from the toplevelmeep
package.  Second, any method of a class is accessible via the standard Python class interface. For example,
meep::fields::step
, which is the function that performs a timestep, is exposed to Python asfields_instance.step()
where a fields instance is usually accessible from Simulation.fields.  C++ constructors are called using the normal Python class instantiation. E.g.,
fields = mp.fields(...)
returns a newmeep::fields
object. Calling destructors 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 it is necessary to know the dimensionality of the problem in C++. The problem dimensions are automatically initialized by Simulation._init_structure
, but if you want to pass vector arguments to C++ before that time you should call Simulation.require_dimensions()
, which infers the dimensions from the cell_size
, k_point
, and dimensions
variables.
Class Reference
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 in any place the superclass is expected.
The meep
package defines several types of classes. The most important of these is the Simulation
class. Classes which are available directly from the meep
package are constructed with:
mp.ClassName(prop1=val1, prop2=val2, ...)
The most numerous are the geometric object classes which are the same as those used in MPB. You can get a list of the available classes (and constants) in the Python interpreter with:
import meep
[x for x in dir(meep) if x[0].isupper()]
More information, including their property types and default values, is available with the standard python help
function: help(mp.ClassName)
.
The following classes are available directly via the meep
package.
Medium
class Medium(object):
This class is used to specify the materials that geometric objects are made of. It
represents an electromagnetic medium which is possibly nonlinear and/or dispersive.
See also Materials. To model a perfectlyconducting metal, use the
predefined metal
object, above. To model imperfect conductors, use a dispersive
dielectric material. See also the Predefined Variables:
metal
, perfect_electric_conductor
, and perfect_magnetic_conductor
.
Material Function
Any function that accepts a Medium
instance can also accept a userdefined Python
function. This allows you to specify the material as an arbitrary function of
position. The function must have one argument, the position Vector3
, and return the
material at that point, which should be a Python Medium
instance. This is
accomplished by passing a function to the material_function
keyword argument in the
Simulation
constructor, or the material
keyword argument in any GeometricObject
constructor. For an example, see Subpixel Smoothing/Enabling Averaging for Material
Function.
Instead of the material
or material_function
arguments, you can also use the
epsilon_func
keyword argument to Simulation
and GeometricObject
, which takes a
function of position that returns the dielectric constant at that point.
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 extra_materials
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 values of and , so
you can only have a finite number of these, whereas can vary
continuously and a matching need not be specified in
extra_materials
. For nonlinear or conductivity materials, your extra_materials
list need not match the actual values of or returned by your material function,
which can vary continuously.
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. Second, if
you care about a broad bandwidth, then all physical materials have
a frequencydependent complex 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
.
def __init__(self,
epsilon_diag=Vector3<1.0, 1.0, 1.0>,
epsilon_offdiag=Vector3<0.0, 0.0, 0.0>,
mu_diag=Vector3<1.0, 1.0, 1.0>,
mu_offdiag=Vector3<0.0, 0.0, 0.0>,
E_susceptibilities=[],
H_susceptibilities=[],
E_chi2_diag=Vector3<0.0, 0.0, 0.0>,
E_chi3_diag=Vector3<0.0, 0.0, 0.0>,
H_chi2_diag=Vector3<0.0, 0.0, 0.0>,
H_chi3_diag=Vector3<0.0, 0.0, 0.0>,
D_conductivity_diag=Vector3<0.0, 0.0, 0.0>,
D_conductivity_offdiag=Vector3<0.0, 0.0, 0.0>,
B_conductivity_diag=Vector3<0.0, 0.0, 0.0>,
B_conductivity_offdiag=Vector3<0.0, 0.0, 0.0>,
epsilon=None,
index=None,
mu=None,
chi2=None,
chi3=None,
D_conductivity=None,
B_conductivity=None,
E_chi2=None,
E_chi3=None,
H_chi2=None,
H_chi3=None,
valid_freq_range=FreqRange(min=1e+20, max=1e+20)):
Creates a Medium
object.

epsilon
[number
] The frequencyindependent isotropic relative permittivity or dielectric constant. Default is 1. You can also useindex=n
as a synonym forepsilon=n*n
; note that this is not really the refractive index if you also specify μ, since the true index is . Usingepsilon=ep
is actually a synonym forepsilon_diag=mp.Vector3(ep, ep, ep)
. 
epsilon_diag
andepsilon_offdiag
[Vector3
] — These properties allow you to specify ε as an arbitrary realsymmetric tensor by giving the diagonal and offdiagonal parts. Specifyingepsilon_diag=Vector3(a, b, c)
and/orepsilon_offdiag=Vector3(u, v, w)
corresponds to a relative permittivity ε tensor Default is the identity matrix ( and ). 
mu
[number
] — The frequencyindependent isotropic relative permeability μ. Default is 1. Usingmu=pm
is actually a synonym formu_diag=mp.Vector3(pm, pm, pm)
. 
mu_diag
andmu_offdiag
[Vector3
] — These properties allow you to specify μ as an arbitrary realsymmetric tensor by giving the diagonal and offdiagonal parts exactly as for ε above. Default is the identity matrix. 
D_conductivity
[number
] — The frequencyindependent electric conductivity . Default is 0. You can also specify a diagonal anisotropic conductivity tensor by using the propertyD_conductivity_diag
which takes aVector3
to give the tensor diagonal. See also Conductivity. 
B_conductivity
[number
] — The frequencyindependent magnetic conductivity . Default is 0. You can also specify a diagonal anisotropic conductivity tensor by using the propertyB_conductivity_diag
which takes aVector3
to give the tensor diagonal. See also Conductivity. 
chi2
[number
] — The nonlinear electric Pockels susceptibility (quadratic nonlinearity). Default is 0. See also Nonlinearity. This is equivalent to settingE_chi2
; alternatively, an analogous magnetic nonlinearity can be specified usingH_chi2
. These are isotropic nonlinearities, but diagonal anisotropic polarizations of the form can be specified withE_chi2_diag
(which defaults to[E_chi2,E_chi2,E_chi2]
). 
chi3
[number
] — The nonlinear electric Kerr susceptibility (cubic nonlinearity). Default is 0. See also Nonlinearity. This is equivalent to settingE_chi3
; alternatively, an analogous magnetic nonlinearity can be specified usingH_chi3
. These are isotropic nonlinearities, but diagonal anisotropic polarizations of the form can be specified withE_chi3_diag
(which defaults to[E_chi3,E_chi3,E_chi3]
). 
E_susceptibilities
[ list ofSusceptibility
class ] — List of dispersive susceptibilities (see below) added to the dielectric constant ε in order to model material dispersion. Defaults to none (empty list). See also Material Dispersion. 
H_susceptibilities
[ list ofSusceptibility
class ] — List of dispersive susceptibilities (see below) added to the permeability μ in order to model material dispersion. Defaults to none (empty list). See also Material Dispersion.
def __repr__(self):
Return repr(self).
def epsilon(self, freq):
Returns the medium's permittivity tensor as a 3x3 Numpy array at the specified
frequency freq
which can be either a scalar, list, or Numpy array. In the case
of a list/array of N frequency points, a Numpy array of size Nx3x3 is returned.
def mu(self, freq):
Returns the medium's permeability tensor as a 3x3 Numpy array at the specified
frequency freq
which can be either a scalar, list, or Numpy array. In the case
of a list/array of N frequency points, a Numpy array of size Nx3x3 is returned.
def transform(self, m):
Transforms epsilon
, mu
, and sigma
of any susceptibilities
by the 3×3 matrix m
. If m
is a rotation
matrix, then the principal axes of
the susceptibilities are rotated by m
. More generally, the susceptibilities χ
are transformed to MχMᵀ/det M, which corresponds to transformation
optics for an
arbitrary curvilinear coordinate transformation with Jacobian matrix M. The
absolute value of the determinant is to prevent inadvertent construction of
lefthanded materials, which are problematic in nondispersive
media.
Susceptibility
class Susceptibility(object):
Parent class for various dispersive susceptibility terms, parameterized by an anisotropic amplitude σ. See Material Dispersion.
def __init__(self,
sigma_diag=Vector3<0.0, 0.0, 0.0>,
sigma_offdiag=Vector3<0.0, 0.0, 0.0>,
sigma=None):
sigma
[number
] — The scale factor σ.
You can also specify an anisotropic σ tensor by using the property sigma_diag
which takes three numbers or a Vector3
to give the σ tensor diagonal, and
sigma_offdiag
which specifies the offdiagonal elements (defaults to 0). That is,
sigma_diag=mp.Vector3(a, b, c)
and sigma_offdiag=mp.Vector3(u, v, w)
corresponds to a σ tensor
LorentzianSusceptibility
class LorentzianSusceptibility(Susceptibility):
Specifies a single dispersive susceptibility of Lorentzian (damped harmonic oscillator) form. See Material Dispersion, with the parameters (in addition to σ):
def __init__(self, frequency=0.0, gamma=0.0, **kwargs):

frequency
[number
] — The resonance frequency . 
gamma
[number
] — The resonance loss rate .
Note: multiple objects with identical values for the frequency
and gamma
but
different sigma
will appear as a single Lorentzian susceptibility term in the
preliminary simulation info output.
DrudeSusceptibility
class DrudeSusceptibility(Susceptibility):
Specifies a single dispersive susceptibility of Drude form. See Material Dispersion, with the parameters (in addition to σ):
def __init__(self, frequency=0.0, gamma=0.0, **kwargs):

frequency
[number
] — The frequency scale factor which multiplies σ (not a resonance frequency). 
gamma
[number
] — The loss rate .
MultilevelAtom
class MultilevelAtom(Susceptibility):
Specifies a multievel atomic susceptibility for modeling saturable gain and
absorption. This is a subclass of E_susceptibilities
which contains two objects: (1)
transitions
: a list of atomic Transition
s (defined below), and (2)
initial_populations
: a list of numbers defining the initial population of each
atomic level. See Materials/Saturable Gain and
Absorption.
def __init__(self,
initial_populations=[],
transitions=[],
**kwargs):
sigma
[number
] — The scale factor σ.
You can also specify an anisotropic σ tensor by using the property sigma_diag
which takes three numbers or a Vector3
to give the σ tensor diagonal, and
sigma_offdiag
which specifies the offdiagonal elements (defaults to 0). That is,
sigma_diag=mp.Vector3(a, b, c)
and sigma_offdiag=mp.Vector3(u, v, w)
corresponds to a σ tensor
Transition
class Transition(object):
def __init__(self,
from_level,
to_level,
transition_rate=0,
frequency=0,
sigma_diag=Vector3<1.0, 1.0, 1.0>,
gamma=0,
pumping_rate=0):
Construct a Transition
.

frequency
[number
] — The radiative transition frequency . 
gamma
[number
] — The loss rate . 
sigma_diag
[Vector3
] — The perpolarization coupling strength . 
from_level
[number
] — The atomic level from which the transition occurs. 
to_level
[number
] — The atomic level to which the transition occurs. 
transition_rate
[number
] — The nonradiative transition rate . Default is 0. 
pumping_rate
[number
] — The pumping rate . Default is 0.
NoisyLorentzianSusceptibility
class NoisyLorentzianSusceptibility(LorentzianSusceptibility):
Specifies a single dispersive susceptibility of Lorentzian (damped harmonic
oscillator) or Drude form. See Material
Dispersion, with the same sigma
, frequency
, and
gamma
parameters, but with an additional Gaussian random noise term (uncorrelated in
space and time, zero mean) added to the P dampedoscillator equation.
def __init__(self, noise_amp=0.0, **kwargs):
noise_amp
[number
] — The noise has rootmean square amplitude σnoise_amp
.
This is 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 and the near field. 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 frequencydomain methods, e.g. via SCUFFEM, as described e.g. here or here.
NoisyDrudeSusceptibility
class NoisyDrudeSusceptibility(DrudeSusceptibility):
Specifies a single dispersive susceptibility of Lorentzian (damped harmonic
oscillator) or Drude form. See Material
Dispersion, with the same sigma
, frequency
, and
gamma
parameters, but with an additional Gaussian random noise term (uncorrelated in
space and time, zero mean) added to the P dampedoscillator equation.
def __init__(self, noise_amp=0.0, **kwargs):
noise_amp
[number
] — The noise has rootmean square amplitude σnoise_amp
.
This is 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 and the near field. 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 frequencydomain methods, e.g. via SCUFFEM, as described e.g. here or here.
GyrotropicLorentzianSusceptibility
class GyrotropicLorentzianSusceptibility(LorentzianSusceptibility):
(Experimental feature) Specifies a single dispersive gyrotropic
susceptibility of Lorentzian (damped harmonic
oscillator) or Drude form. Its
parameters are sigma
, frequency
, and gamma
, which have the usual
meanings, and an additional 3vector bias
:
def __init__(self, bias=Vector3<0.0, 0.0, 0.0>, **kwargs):
bias
[Vector3
] — The gyrotropy vector. Its direction determines the orientation of the gyrotropic response, and the magnitude is the precession frequency .
GyrotropicDrudeSusceptibility
class GyrotropicDrudeSusceptibility(DrudeSusceptibility):
(Experimental feature) Specifies a single dispersive gyrotropic
susceptibility of Lorentzian (damped harmonic
oscillator) or Drude form. Its
parameters are sigma
, frequency
, and gamma
, which have the usual
meanings, and an additional 3vector bias
:
def __init__(self, bias=Vector3<0.0, 0.0, 0.0>, **kwargs):
bias
[Vector3
] — The gyrotropy vector. Its direction determines the orientation of the gyrotropic response, and the magnitude is the precession frequency .
GyrotropicSaturatedSusceptibility
class GyrotropicSaturatedSusceptibility(Susceptibility):
(Experimental feature) Specifies a single dispersive gyrotropic
susceptibility governed by a linearized
LandauLifshitzGilbert
equation.
This class takes parameters sigma
, frequency
, and gamma
, whose meanings are
different from the Lorentzian and Drude case. It also takes a 3vector bias
parameter and an alpha
parameter:
def __init__(self,
bias=Vector3<0.0, 0.0, 0.0>,
frequency=0.0,
gamma=0.0,
alpha=0.0,
**kwargs):

sigma
[number
] — The coupling factor between the polarization and the driving field. In magnetic ferrites, this is the Larmor precession frequency at the saturation field. 
frequency
[number
] — The Larmor precession frequency, . 
gamma
[number
] — The loss rate in the offdiagonal response. 
alpha
[number
] — The loss factor in the diagonal response. Note that this parameter is dimensionless and contains no 2π factor. 
bias
[Vector3
] — Vector specifying the orientation of the gyrotropic response. Unlike the similarlynamedbias
parameter for the gyrotropic Lorentzian/Drude susceptibilities, the magnitude is ignored; instead, the relevant precession frequencies are determined by thesigma
andfrequency
parameters.
Vector3
class Vector3(object):
Properties:
x
, y
, z
[float
or complex
] — The x
, y
, and z
components of the
vector. Generally, functions that take a Vector3
as an argument will accept an
iterable (e.g., a tuple or list) and automatically convert to a Vector3
.
def __add__(self, other):
Return the sum of the two vectors.
v3 = v1 + v2
def __eq__(self, other):
Returns whether or not the two vectors are numerically equal. Beware of using this
function after operations that may have some error due to the finite precision of
floatingpoint numbers; use close
instead.
v1 == v2
def __init__(self, x=0.0, y=0.0, z=0.0):
Create a new Vector3
with the given components. All three components default to
zero. This can also be represented simply as (x,y,z)
or [x,y,z]
.
def __mul__(self, other):
If other
is a Vector3
, returns the dot product of v1
and other
. If other
is a number, then v1
is scaled by the number.
c = v1 * other
def __ne__(self, other):
Returns whether or not the two vectors are numerically unequal. Beware of using
this function after operations that may have some error due to the finite
precision of floatingpoint numbers; use close
instead.
v1 != v2
def __repr__(self):
Return repr(self).
def __rmul__(self, other):
If other
is a Vector3
, returns the dot product of v1
and other
. If other
is a number, then v1
is scaled by the number.
c = other * v1
def __sub__(self, other):
Return the difference of the two vectors.
v3 = v1  v2
def cdot(self, v):
Returns the conjugated dot product: conj(self)
dot v
.
def close(self, v, tol=1e07):
Returns whether or not the corresponding components of the self
and v
vectors
are within tol
of each other. Defaults to 1e7.
v1.close(v2, [tol])
def cross(self, v):
Return the cross product of self
and v
.
v3 = v1.cross(v2)
def dot(self, v):
Returns the dot product of self
and v
.
v3 = v1.dot(v2)
def norm(self):
Returns the length math.sqrt(abs(self.dot(self)))
of the given vector.
v2 = v1.norm()
def rotate(self, axis, theta):
Returns the vector rotated by an angle theta
(in radians) in the righthand
direction around the axis
vector (whose length is ignored). You may find the
python functions math.degrees
and math.radians
useful to convert angles
between degrees and radians.
v2 = v1.rotate(axis, theta)
def unit(self):
Returns a unit vector in the direction of the vector.
v2 = v1.unit()
GeometricObject
class GeometricObject(object):
This class, and its descendants, are used to specify the solid geometric objects that form the dielectric structure being simulated.
In a 2d calculation, only the intersections of the objects with the plane are considered.
Geometry Utilities
See the MPB documentation for utility functions to help manipulate geometric objects.
Examples
These are some examples of geometric objects created using some GeometricObject
subclasses:
# A cylinder of infinite radius and height 0.25 pointing along the x axis,
# centered at the origin:
cyl = mp.Cylinder(center=mp.Vector3(0,0,0), height=0.25, radius=mp.inf,
axis=mp.Vector3(1,0,0), material=mp.Medium(index=3.5))
# 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):
ell = mp.Ellipsoid(center=mp.Vector3(0,0,0), size=mp.Vector3(0.8,0.2,0.2),
e1=Vector3(1,1,1), e2=Vector3(0,1,1), e3=Vector3(2,1,1),
material=mp.Medium(epsilon=13))
# A unit cube of material metal with a spherical air hole of radius 0.2 at
# its center, the whole thing centered at (1,2,3):
geometry=[mp.Block(center=Vector3(1,2,3), size=Vector3(1,1,1), material=mp.metal),
mp.Sphere(center=Vector3(1,2,3), radius=0.2, material=mp.air)]
# A hexagonal prism defined by six vertices centered on the origin
# of material crystalline silicon (from the materials library)
vertices = [mp.Vector3(1,0),
mp.Vector3(0.5,math.sqrt(3)/2),
mp.Vector3(0.5,math.sqrt(3)/2),
mp.Vector3(1,0),
mp.Vector3(0.5,math.sqrt(3)/2),
mp.Vector3(0.5,math.sqrt(3)/2)]
geometry = [mp.Prism(vertices, height=1.5, center=mp.Vector3(), material=cSi)]
def __init__(self,
material=Medium(),
center=Vector3<0.0, 0.0, 0.0>,
epsilon_func=None):
Construct a GeometricObject
.

material
[Medium
class or function ] — The material that the object is made of (usually some sort of dielectric). Uses defaultMedium
. If a function is supplied, it must take one argument and return a PythonMedium
. 
epsilon_func
[ function ] — A function that takes one argument (aVector3
) and returns the dielectric constant at that point. Can be used instead ofmaterial
. Default isNone
. 
center
[Vector3
] — Center point of the object. Defaults to(0,0,0)
.
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 and can be specified in a subclass's
constructor via keyword arguments.
def info(self, indent_by=0):
Displays all properties and current values of a GeometricObject
, indented by
indent_by
spaces (default is 0).
def shift(self, vec):
Shifts the object's center
by vec
(Vector3
), returning a new object.
This can also be accomplished via the +
operator:
geometric_obj + Vector3(10,10,10)
Using +=
will shift the object in place.
Sphere
class Sphere(GeometricObject):
Represents a sphere.
Properties:
radius
[number
] — Radius of the sphere. No default value.
def __init__(self, radius, **kwargs):
Constructs a Sphere
Cylinder
class Cylinder(GeometricObject):
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 toVector3(x=0, y=0, z=1)
.
def __init__(self,
radius,
axis=Vector3<0.0, 0.0, 1.0>,
height=1e+20,
**kwargs):
Constructs a Cylinder
.
Wedge
class Wedge(Cylinder):
Represents a cylindrical wedge.
def __init__(self,
radius,
wedge_angle=6.283185307179586,
wedge_start=Vector3<1.0, 0.0, 0.0>,
**kwargs):
Constructs a Wedge
.
Cone
class Cone(Cylinder):
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
. The center
of a cone
is halfway between the two circular ends.
def __init__(self, radius, radius2=0, **kwargs):
Construct a Cone
.
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
class Block(GeometricObject):
A parallelepiped (i.e., a brick, possibly with nonorthogonal axes).
def __init__(self,
size,
e1=Vector3<1.0, 0.0, 0.0>,
e2=Vector3<0.0, 1.0, 0.0>,
e3=Vector3<0.0, 0.0, 1.0>,
**kwargs):
Construct a Block
.

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
class Ellipsoid(Block):
An ellipsoid. This is actually a subclass of Block
, and inherits all the same
properties, but defines an ellipsoid inscribed inside the block.
def __init__(self, **kwargs):
Construct an Ellipsoid
.
Prism
class Prism(GeometricObject):
Polygonal prism type.
def __init__(self,
vertices,
height,
axis=Vector3<0.0, 0.0, 1.0>,
center=None,
sidewall_angle=0,
**kwargs):
Construct a Prism
.

vertices
[list ofVector3
] — The vertices that make up the prism. They must lie in a plane that's perpendicular to theaxis
. Note that infinite lengths are not supported. To simulate infinite geometry, just extend the edge of the prism beyond the cell. 
height
[number
] — The prism thickness, extruded in the direction ofaxis
.mp.inf
can be used for infinite height. No default value. 
axis
[Vector3
] — The axis perpendicular to the prism. Defaults toVector3(0,0,1)
. 
center
[Vector3
] — Ifcenter
is not specified, then the coordinates of thevertices
define the bottom of the prism with the top of the prism being at the same coordinates shifted byheight*axis
. Ifcenter
is specified, thencenter
is the coordinates of the centroid of all the vertices (top and bottom) of the resulting 3d prism so that the coordinates of thevertices
are shifted accordingly.
Matrix
class Matrix(object):
The Matrix
class represents a 3x3 matrix with c1, c2, and c3 as its columns.
m.transpose()
m.getH() or m.H
m.determinant()
m.inverse()
Return the transpose, adjoint (conjugate transpose), determinant, or inverse of the given matrix.
m1 + m2
m1  m2
m1 * m2
Return the sum, difference, or product of the given matrices.
v * m
m * v
Returns the Vector3
product of the matrix m
by the vector v
, with the vector
multiplied on the left or the right respectively.
s * m
m * s
Scales the matrix m
by the number s
.
def __init__(self,
c1=Vector3<0.0, 0.0, 0.0>,
c2=Vector3<0.0, 0.0, 0.0>,
c3=Vector3<0.0, 0.0, 0.0>,
diag=Vector3<0.0, 0.0, 0.0>,
offdiag=Vector3<0.0, 0.0, 0.0>):
Constructs a Matrix
.
def __repr__(self):
Return repr(self).
Related function:
def get_rotation_matrix(axis, theta):
Like Vector3.rotate
, except returns the (unitary) rotation matrix that performs the
given rotation. i.e., get_rotation_matrix(axis, theta) * v
produces the same result
as v.rotate(axis, theta)
.

axis
[Vector3
] — The vector around which the rotation is applied in the righthand direction. 
theta
[number
] — The rotation angle (in radians).
Symmetry
class Symmetry(object):
This class is used for the symmetries
input variable to specify symmetries which
must preserve both the structure and the sources. Any number of symmetries can be
exploited simultaneously but there is no point in specifying redundant symmetries: the
cell can be reduced by at most a factor of 4 in 2d and 8 in 3d. See also Exploiting
Symmetry. 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:
The specific symmetry subclasses are:
Mirror
— A mirror symmetry plane. direction
is the direction normal to the
mirror plane.
Rotate2
— A 180° (twofold) rotational symmetry (a.k.a. ). direction
is
the axis of the rotation.
Rotate4
— A 90° (fourfold) rotational symmetry (a.k.a. ). direction
is
the axis of the rotation.
def __init__(self, direction, phase=1):
Construct a Symmetry
.

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
, orZ
(only Cartesian/grid directions are allowed). No default value. 
phase
[complex
] — An additional phase to multiply the fields by when operating the symmetry on them. Default is +1, 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.
Rotate2
class Rotate2(Symmetry):
A 180° (twofold) rotational symmetry (a.k.a. ). direction
is the axis of the
rotation.
Rotate4
class Rotate4(Symmetry):
A 90° (fourfold) rotational symmetry (a.k.a. ). direction
is the axis of the
rotation.
Mirror
class Mirror(Symmetry):
A mirror symmetry plane. direction
is the direction normal to the mirror plane.
Identity
class Identity(Symmetry):
PML
class PML(object):
This class is used for specifying the PML absorbing boundary layers around the cell,
if any, via the boundary_layers
input variable. See also Perfectly Matched
Layers. boundary_layers
can be zero or more PML
objects, with multiple objects allowing you to specify different PML layers on
different boundaries. The class represents a single PML layer specification, which
sets up one or more PML layers around the boundaries according to the following
properties.
def __init__(self,
thickness,
direction=mp.ALL,
side=mp.ALL,
R_asymptotic=1e15,
mean_stretch=1.0,
pml_profile=DefaultPMLProfile):

thickness
[number
] — The spatial thickness of the PML layer which extends from the boundary towards the inside of the 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. ifX
, 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
[side
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 boundary. Default is the special valueALL
, which puts PML layers on both sides. 
R_asymptotic
[number
] — The asymptotic reflection in the limit of infinite resolution or infinite PML thickness, for reflections 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. Default value is 10^{−15}, which should suffice for most purposes. You want to set this to be small enough so that waves propagating within the PML are attenuated sufficiently, but makingR_asymptotic
too small will increase the numerical reflection due to discretization. 
pml_profile
[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, withpml_profile
one can specify an arbitrary PML "profile" function 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 . 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 Optics Express, Vol. 16, pp. 1137692, 2008. For example, one can use a cubic profile by specifyingpml_profile=lambda u: u*u*u
.
Absorber
class Absorber(PML):
Instead of a PML
layer, there is an alternative class called Absorber
which is a
dropin replacement for PML
. For example, you can do
boundary_layers=[mp.Absorber(thickness=2)]
instead of
boundary_layers=[mp.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 pml_profile
(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 computationally less expensive and generally just as good. See Optics Express, Vol. 16, pp. 1137692, 2008.
 PML can lead to divergent fields for certain waveguides with "backwardwave" modes; this can readily occur in metals with surface plasmons, and a scalar absorber is your only choice. See Physical Review E, Vol. 79, 065601, 2009.
 PML can fail if you have a waveguide hitting the edge of your cell at an angle. See J. Computational Physics, Vol. 230, pp. 236977, 2011.
Source
class Source(object):
The Source
class is used to specify the current sources via the Simulation.sources
attribute. Note that all sources in Meep are separable in time and space, i.e. of the
form for some functions
and . 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 Simulation.force_complex_fields
), only the real part of
the current source is used.
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 FAQ and in Section 4.4 ("Currents and Fields: The Local Density of States") in Chapter 4 ("Electromagnetic Wave Source Conditions") of the book Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology.
def __init__(self,
src,
component,
center=None,
volume=None,
size=Vector3<0.0, 0.0, 0.0>,
amplitude=1.0,
amp_func=None,
amp_func_file='',
amp_data=None):
Construct a Source
.

src
[SourceTime
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.mp.Ex
,mp.Ey
, etcetera for an electriccharge current, andmp.Hx
,mp.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 cell. No default. 
size
[Vector3
] — The size of the current distribution along each direction of the cell. Default is(0,0,0)
: a pointdipole source. 
volume
[Volume
] — Ameep.Volume
can be used to specify the source region instead of acenter
and asize
. 
amplitude
[complex
] — An overall complex amplitude multiplying the current source. Default is 1.0. Note that specifying a complexamplitude
imparts a phase shift to the real part of the overall current and thus does not require using complex fields for the entire simulation (viaforce_complex_fields=True
). 
amp_func
[function
] — A Python function of a single argument, that takes aVector3
giving a position and returns a complex current amplitude for that point. The position argument is relative to thecenter
of the current source, so that you can move your current around without changing your function. Default isNone
, 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. 
amp_func_file
[string
] — String of the formpath_to_h5_file.h5:dataset
. The.h5
extension is optional. Meep will read the HDF5 file and create an amplitude function that interpolates into the grid specified by the file. Meep expects the data to be split into real and imaginary parts, so in the above example it will look fordataset.re
anddataset.im
in the filepath_to_h5_file.h5
. Defaults to the empty string. 
amp_data
[numpy.ndarray with dtype=numpy.complex128
] — Likeamp_func_file
above, but instead of interpolating into an HDF5 file, interpolates into a complex NumPy array. The array should be three dimensions. For a 2d simulation, just pass 1 for the third dimension, e.g.,arr = np.zeros((N, M, 1), dtype=np.complex128)
. Defaults toNone
.
As described in Section 4.2 ("Incident Fields and Equivalent Currents") in
Chapter 4 ("Electromagnetic Wave Source
Conditions") of the book Advances in FDTD Computational Electrodynamics:
Photonics and Nanotechnology,
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 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
SourceTime
class SourceTime(object):
This is the parent for classes describing the time dependence of sources; it should not be instantiated directly.
EigenModeSource
class EigenModeSource(Source):
This is a subclass of Source
and has all of the properties of Source
above.
However, you normally do not specify a component
. Instead of component
, 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 the size
and center
of the
source, with the modes computed as if the source region were repeated periodically in
all directions. If an amplitude
and/or amp_func
are supplied, they are
multiplied by this current profile. The desired eigenmode and other features are
specified by the properties shown in __init__
.
Eigenmode sources are normalized so that in the case of a timeharmonic simulation with all sources and fields having monochromatic time dependence where is the frequency of the eigenmode, the total timeaverage power of the fields — the integral of the normal Poynting vector over the entire crosssectional line or plane — is equal to 1. This convention has two use cases:

For frequencydomain calculations involving a
ContinuousSource
time dependence, the timeaverage power of the fields is 1. 
For timedomain calculations involving a time dependence which is typically a Gaussian, the amplitude of the fields at frequency will be multiplied by , the Fourier transform of , while fieldbilinear quantities like the Poynting flux and energy density are multiplied by . For the particular case of a Gaussian time dependence, the Fourier transform at can be obtained via the
fourier_transform
class method.
In either case, the eig_power
method returns the total power at frequency f
.
However, for a userdefined CustomSource
, eig_power
will not
include the factor since Meep does not know the Fourier
transform of your source function . You will have to multiply by this yourself
if you need it.
Note: Due to discretization effects, the normalization of eigenmode sources to
yield unit power transmission is only approximate: at any finite resolution, the power
of the fields as measured using DFT flux monitors will not precisely
match that of calling eig_power
but will rather include discretization errors that
decrease with resolution. Generally, the most reliable procedure is to normalize your
calculations by the power computed in a separate normalization run at the same
resolution, as shown in several of the tutorial examples.
Note that Meep's MPB interface 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 SourceTime
object (Source.src
), which specifies the time dependence of the
source, can be one of ContinuousSource
, GaussianSource
or CustomSource
.
def __init__(self,
src,
center=None,
volume=None,
eig_lattice_size=None,
eig_lattice_center=None,
component=mp.ALL_COMPONENTS,
direction=mp.AUTOMATIC,
eig_band=1,
eig_kpoint=Vector3<0.0, 0.0, 0.0>,
eig_match_freq=True,
eig_parity=mp.NO_PARITY,
eig_resolution=0,
eig_tolerance=1e12,
**kwargs):
Construct an EigenModeSource
.

eig_band
[integer
orDiffractedPlanewave
] — Either the index n (1,2,3,...) of the desired band ω_{n}(k) to compute in MPB where 1 denotes the lowestfrequency band at a given k point, and so on, or alternatively a diffracted planewave in homogeneous media. 
direction
[mp.X
,mp.Y
, ormp.Z;
defaultmp.AUTOMATIC
],eig_match_freq
[boolean;
defaultTrue
],eig_kpoint
[Vector3
] — By default (ifeig_match_freq
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'sfind_k
functionality. Alternatively, ifeig_kpoint
is supplied, it is used as an initial guess for k. By default,direction
is the direction normal to the source region, assumingsize
is –1 dimensional in a dimensional simulation (e.g. a plane in 3d). Ifdirection
is set tomp.NO_DIRECTION
, theneig_kpoint
is not only the initial guess and the search direction of the k vectors, but is also taken to be the direction of the waveguide, allowing you to launch modes in oblique ridge waveguides (not perpendicular to the source plane). Ifeig_match_freq
isFalse
, then the k vector of the desired mode is specified witheig_kpoint
(in Meep units of 2π/(unit length)). Also, the eigenmode frequency computed by MPB overwrites thefrequency
parameter of thesrc
property for aGaussianSource
andContinuousSource
but notCustomSource
(thewidth
or any other parameter ofsrc
is unchanged). By default, the k components in the plane of the source region are zero. However, if the source region spans the entire cell in some directions, and the cell has Blochperiodic boundary conditions via thek_point
parameter, then the mode's k components in those directions will matchk_point
so that the mode satisfies the Meep boundary conditions, regardless ofeig_kpoint
. Note that once k is either found by MPB, or specified byeig_kpoint
, the field profile used to create the current sources corresponds to the Bloch mode, , multiplied by the appropriate exponential factor, . 
eig_parity
[mp.NO_PARITY
(default),mp.EVEN_Z
,mp.ODD_Z
,mp.EVEN_Y
,mp.ODD_Y
] — The parity (= polarization in 2d) of the mode to calculate, assuming the structure has and/or mirror symmetry in the source region, with respect to thecenter
of the source region. (In particular, it does not matter if your simulation as a whole has that symmetry, only the cross section where you are introducing the source.) If the structure has both and mirror symmetry, you can combine more than one of these, e.g.EVEN_Z + ODD_Y
. Default isNO_PARITY
, 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. 
eig_resolution
[integer
, defaults to2*resolution
] — The spatial resolution to use in MPB for the eigenmode calculations. This defaults to twice the Meepresolution
in which case the structure is linearly interpolated from the Meep pixels. 
eig_tolerance
[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 toALL_COMPONENTS
] — 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 Section 4.2 ("Incident Fields and Equivalent Currents") in Chapter 4 ("Electromagnetic Wave Source Conditions") of the book Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology. If you specify acomponent
property, however, you can include only one component of these currents if you wish. Most users won't need this feature. 
eig_lattice_size
[Vector3
],eig_lattice_center
[Vector3
] — Normally, the MPB computational unit cell is the same as the source volume given by thesize
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, 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 optionaleig_lattice_size
andeig_lattice_center
, 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.
def eig_power(self, freq):
Returns the total power of the fields from the eigenmode source at frequency freq
.
GaussianBeamSource
class GaussianBeamSource(Source):
This is a subclass of Source
and has all of the properties of Source
above. However, the component
parameter of the Source
object is ignored. The Gaussian beam is a transverse electromagnetic mode for which the source region must be a line (in 2d) or plane (in 3d). For a beam polarized in the direction with propagation along , the electric field is defined by where is the radial distance from the center axis of the beam, is the axial distance from the beam's focus (or "waist"), is the wavenumber (for a freespace wavelength and refractive index of the homogeneous, lossless medium in which the beam propagates), is the electric field amplitude at the origin, is the radius at which the field amplitude decays by of its axial values, is the beam waist radius, and is the radius of curvature of the beam's wavefront at . The only independent parameters that need to be specified are , , , and the location of the beam focus (i.e., the origin: ).
(In 3d, we use a "complex pointsource" method to define a source that generates an exact Gaussianbeam solution. In 2d, we currently use the simple approximation of taking a crosssection of the 3d beam. In both cases, the beam is most accurate near the source's center frequency.)
The SourceTime
object (Source.src
), which specifies the time dependence of the source, should normally be a narrowband ContinuousSource
or GaussianSource
. (For a CustomSource
, the beam frequency is determined by the source's center_frequency
parameter.)
def __init__(self,
src,
center=None,
volume=None,
component=mp.ALL_COMPONENTS,
beam_x0=Vector3<0.0, 0.0, 0.0>,
beam_kdir=Vector3<0.0, 0.0, 0.0>,
beam_w0=None,
beam_E0=Vector3<0.0, 0.0, 0.0>,
**kwargs):
Construct a GaussianBeamSource
.

beam_x0
[Vector3
] — The location of the beam focus relative to the center of the source. The beam focus does not need to lie within the source region (i.e., the beam focus can be anywhere, inside or outside the cell, independent of the position of the source). 
beam_kdir
[Vector3
] — The propagation direction of the beam. The length is ignored. The wavelength of the beam is determined by the center frequency of theSource.src
object and the refractive index (real part only) at the center of the source region. 
beam_w0
[number
] — The beam waist radius. 
beam_E0
[Vector3
] — The polarization vector of the beam. Elements can be complex valued (i.e., for circular polarization). The polarization vector must be parallel to the source region in order to generate a transverse mode.
ContinuousSource
class ContinuousSource(SourceTime):
A continuouswave (CW) source is proportional to , possibly with a smooth (exponential/tanh) turnon/turnoff. In practice, the CW source never produces an exact singlefrequency response.
def __init__(self,
frequency=None,
start_time=0,
end_time=1e+20,
width=0,
fwidth=inf,
cutoff=3.0,
wavelength=None,
**kwargs):
Construct a ContinuousSource
.

frequency
[number
] — The frequency f in units of /distance or ω in units of 2π/distance. See Units. No default value. You can instead specifywavelength=x
orperiod=x
, which are both a synonym forfrequency=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 ). 
end_time
[number
] — The end time for the source. Default is 10^{20} (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 specifyfwidth=x
, which is a synonym forwidth=1/x
(i.e. the frequency width is proportional to the inverse of the temporal width). 
slowness
[number
] — Controls how far into the exponential tail of the tanh function the source turns on. Default is 3.0. A larger value means that the source turns on more gradually at the beginning. 
is_integrated
[boolean
] — IfTrue
, the source is the integral of the current (the dipole moment) which oscillates but does not increase for a sinusoidal current. In practice, there is little difference between integrated and nonintegrated sources except for planewaves extending into PML. Default isFalse
.
GaussianSource
class GaussianSource(SourceTime):
A Gaussianpulse source roughly proportional to . 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.
def __init__(self,
frequency=None,
width=0,
fwidth=inf,
start_time=0,
cutoff=5.0,
wavelength=None,
**kwargs):
Construct a GaussianSource
.

frequency
[number
] — The center frequency in units of /distance (or ω in units of 2π/distance). See Units. No default value. You can instead specifywavelength=x
orperiod=x
, which are both a synonym forfrequency=1/x
; i.e. 1/ω in these units is the vacuum wavelength or the temporal period. 
width
[number
] — The width used in the Gaussian. No default value. You can instead specifyfwidth=x
, which is a synonym forwidth=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 ). This is not the time of the peak. See below. 
cutoff
[number
] — How manywidth
s the current decays for before it is cut off and set to zero — this applies for both turnon and turnoff of the pulse. Default is 5.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 =start_time + cutoff*width
. 
is_integrated
[boolean
] — IfTrue
, the source is the integral of the current (the dipole moment) which is guaranteed to be zero after the current turns off. In practice, there is little difference between integrated and nonintegrated sources except for planewaves extending into PML. Default isFalse
. 
fourier_transform(f)
— Returns the Fourier transform of the current evaluated at frequencyf
(ω=2πf
) given by: where is the current (not the dipole moment). In this formula, is thefwidth
of the source, is times itsfrequency,
and is the peak time discussed above. Note that this does not include anyamplitude
oramp_func
factor that you specified for the source.
CustomSource
class CustomSource(SourceTime):
A userspecified source function . 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 run
functions like until_after_sources
to work, since they need to know when your source
turns off. To use a custom source within an EigenModeSource
, you must specify the
center_frequency
parameter, since Meep does not know the frequency content of the
CustomSource
. The resultant eigenmode is calculated at this frequency only. For a
demonstration of a linearchirped pulse, see
examples/chirped_pulse.py
.
def __init__(self,
src_func,
start_time=1e+20,
end_time=1e+20,
center_frequency=0,
**kwargs):
Construct a CustomSource
.

src_func
[function
] — The function specifying the timedependence 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 10^{20}: turn on at . Note, however, that the simulation normally starts at with zero fields as the initial condition, so there is implicitly a sharp turnon at whether you specify it or not. 
end_time
[number
] — The end time for the source. Default is 10^{20} (never turn off). 
is_integrated
[boolean
] — IfTrue
, the source is the integral of the current (the dipole moment) which is guaranteed to be zero after the current turns off. In practice, there is little difference between integrated and nonintegrated sources except for planewaves extending into PML. Default isFalse
. 
center_frequency
[number
] — Optional center frequency so that theCustomSource
can be used within anEigenModeSource
. Defaults to 0.
FluxRegion
class FluxRegion(object):
A FluxRegion
object is used with add_flux
to specify a region in
which Meep should accumulate the appropriate Fouriertransformed fields in order to
compute a flux spectrum. It represents a region (volume, plane, line, or point) in
which to compute the integral of the Poynting vector of the Fouriertransformed
fields. ModeRegion
is an alias for FluxRegion
for use with add_mode_monitor
.
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.
def __init__(self,
center=None,
size=Vector3<0.0, 0.0, 0.0>,
direction=mp.AUTOMATIC,
weight=1.0,
volume=None):
Construct a FluxRegion
object.

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.mp.X
,mp.Y
, etcetera). 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
[complex
] — A weight factor to multiply the flux by when it is computed. Default is 1.0. 
volume
[Volume
] — Ameep.Volume
can be used to specify the flux region instead of acenter
and asize
.
EnergyRegion
class EnergyRegion(FluxRegion):
A region (volume, plane, line, or point) in which to compute the integral of the energy density of the Fouriertransformed fields. Its properties are:

center
[Vector3
] — The center of the energy region (no default). 
size
[Vector3
] — The size of the energy region along each of the coordinate axes. Default is (0,0,0): a single point. 
weight
[complex
] — A weight factor to multiply the energy density by when it is computed. Default is 1.0.
ForceRegion
class ForceRegion(FluxRegion):
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
[complex
] — A weight factor to multiply the force by when it is computed. Default is 1.0. 
volume
[Volume
] — Ameep.Volume
can be used to specify the force region instead of acenter
and asize
.
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.
Volume
class Volume(object):
Many Meep functions require you to specify a volume in space, corresponding to the C++
type meep::volume
. This class 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. Any method that accepts such a volume
also accepts center
and size
keyword arguments. If these are specified instead of
the volume, the library will construct a volume for you. Alternatively, you can
specify a list of Vector3
vertices using the vertices
parameter. The center
and
size
will automatically be computed from this list.
def __init__(self,
center=Vector3<0.0, 0.0, 0.0>,
size=Vector3<0.0, 0.0, 0.0>,
dims=2,
is_cylindrical=False,
vertices=[]):
Construct a Volume.
Related function:
def get_center_and_size(vol):
Utility function that takes a meep::volume
vol
and returns the center and size of
the volume as a tuple of Vector3
.
DftObj
class DftObj(object):
Wrapper around dft objects that allows delayed initialization of the structure.
When splitting the structure into chunks for parallel simulations, we want to know all of the details of the simulation in order to ensure that each processor gets a similar amount of work. The problem with DFTs is that the 'add_flux' style methods immediately initialize the structure and fields. So, if the user adds multiple DFT objects to the simulation, the load balancing code only knows about the first one and can't split the work up nicely. To circumvent this, we delay the execution of the 'add_flux' methods as late as possible. When 'add_flux' (or add_near2far, etc.) is called, we:

Create an instance of the appropriate subclass of DftObj (DftForce, DftFlux, etc.). Set its args property to the list of arguments passed to add_flux, and set its func property to the 'real' add_flux, which is prefixed by an underscore.

Add this DftObj to the list Simulation.dft_objects. When we actually run the simulation, we call Simulation._evaluate_dft_objects, which calls dft.func(*args) for each dft in the list.
If the user tries to access a property or call a function on the DftObj before Simulation._evaluate_dft_objects is called, then we initialize the C++ object through swigobj_attr and return the property they requested.
def __init__(self, func, args):
Construct a DftObj
.
DftFlux
class DftFlux(DftObj):
def __init__(self, func, args):
Construct a DftFlux
.
DftForce
class DftForce(DftObj):
def __init__(self, func, args):
Construct a DftForce
.
DftNear2Far
class DftNear2Far(DftObj):
def __init__(self, func, args):
Construct a DftNear2Far
.
def flux(self, direction, where, resolution):
Given a Volume
where
(may be 0d, 1d, 2d, or 3d) and a resolution
(in grid
points / distance unit), compute the far fields in where
(which may lie
outside the cell) in a grid with the given resolution (which may differ from the
FDTD solution) and return its Poynting flux in direction
as a list. The dataset
is a 1d array of nfreq
dimensions.
DftEnergy
class DftEnergy(DftObj):
def __init__(self, func, args):
Construct a DftEnergy
.
DftFields
class DftFields(DftObj):
def __init__(self, func, args):
Construct a DftFields
.
Animate2D
class Animate2D(object):
A class used to record the fields during timestepping (i.e., a run
function). The object is initialized prior to timestepping by specifying the
simulation object and the field component. The object can then be passed to any
stepfunction modifier. For example, one can record the
E_{z} fields at every one time unit using:
animate = mp.Animate2D(sim,
fields=mp.Ez,
realtime=True,
field_parameters={'alpha':0.8, 'cmap':'RdBu', 'interpolation':'none'},
boundary_parameters={'hatch':'o', 'linewidth':1.5, 'facecolor':'y', 'edgecolor':'b', 'alpha':0.3})
sim.run(mp.at_every(1,animate),until=25)
By default, the object saves each frame as a PNG image into memory (not disk). This is
typically more memory efficient than storing the actual fields. If the user sets the
normalize
argument, then the object will save the actual field information as a
NumPy array to be normalized for post processing. The fields of a figure can also be
updated in realtime by setting the realtime
flag. This does not work for
IPython/Jupyter notebooks, however.
Once the simulation is run, the animation can be output as an interactive JSHTML object, an mp4, or a GIF.
Multiple Animate2D
objects can be initialized and passed to the run function to
track different volume locations (using mp.in_volume
) or field components.
def __call__(self, sim, todo):
Call self as a function.
def __init__(self,
sim,
fields,
f=None,
realtime=False,
normalize=False,
plot_modifiers=None,
**customization_args):
Construct an Animate2D
object.

sim
— Simulation object. 
fields
— Field component to record at each time instant. 
f=None
— Optionalmatplotlib
figure object that the routine will update on each call. If not supplied, then a new one will be created upon initialization. 
realtime=False
— Whether or not to update a figure window in realtime as the simulation progresses. Disabled by default. Not compatible with IPython/Jupyter notebooks. 
normalize=False
— Records fields at each time step in memory in a NumPy array and then normalizes the result by dividing by the maximum field value at a single point in the cell over all the time snapshots. 
plot_modifiers=None
— A list of functions that can modify the figure'saxis
object. Each function modifier accepts a single argument, anaxis
object, and must return that same axis object. The following modifier changes thexlabel
:
def mod1(ax):
ax.set_xlabel('Testing')
return ax
plot_modifiers = [mod1]
**customization_args
— Customization keyword arguments passed toplot2D()
(i.e.labels
,eps_parameters
,boundary_parameters
, etc.)
def to_gif(self, fps, filename):
Generates and outputs a GIF file of the animation with the filename, filename
,
and the frame rate, fps
. Note that GIFs are significantly larger than mp4 videos
since they don't use any compression. Artifacts are also common because the GIF
format only supports 256 colors from a predefined color palette. Requires
ffmpeg
.
def to_jshtml(self, fps):
Outputs an interactable JSHTML animation object that is embeddable in Jupyter
notebooks. The object is packaged with controls to manipulate the video's
playback. User must specify a frame rate fps
in frames per second.
def to_mp4(self, fps, filename):
Generates and outputs an mp4 video file of the animation with the filename,
filename
, and the frame rate, fps
. Default encoding is h264 with yuv420p
format. Requires ffmpeg
.
Harminv
class Harminv(object):
Harminv
is implemented as a class with a __call__
method,
which allows it to be used as a step function that collects field data from a given
point and runs Harminv on that data to extract
the frequencies, decay rates, and other information.
See __init__
for details about constructing a Harminv
.
Important: normally, you should only use Harminv to analyze data after the
sources are off. Wrapping it in after_sources(mp.Harminv(...))
is sufficient.
In particular, Harminv takes the time series 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 Harminv.modes
, which is a list of tuples holding
the frequency, amplitude, and error of the modes. Given one of these tuples (e.g.,
first_mode = harminv_instance.modes[0]
), you can extract its various components:

freq
— The real part of frequency ω (in the usual Meep 2πc units). 
decay
— The imaginary part of the frequency ω. 
Q
— The dimensionless lifetime, or quality factor defined as . 
amp
— The complex amplitude . 
err
— 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 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, [m.freq for m in harminv_instance.modes]
gives a list of the real parts
of the frequencies. Be sure to save a reference to the Harminv
instance if you wish
to use the results after the simulation:
sim = mp.Simulation(...)
h = mp.Harminv(...)
sim.run(mp.after_sources(h))
# do something with h.modes
def __call__(self, sim, todo):
Allows a Haminv instance to be used as a step function.
def __init__(self, c, pt, fcen, df, mxbands=None):
Construct a Harminv object.
A Harminv
is a step function that collects data from the field component c
(e.g. E_{x}, 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 commadelimited text, and also storing them to the variable
Harminv.modes
. The optional argument mxbands
is the maximum number of modes to
search for. Defaults to 100.
Verbosity
class Verbosity(object):
A class to help make accessing and setting the global verbosity level a bit more pythonic. It manages one or more verbosity flags that are located in the C/C++ libraries used by Meep.
The verbosity levels are:
 0: minimal output
 1: a little
 2: a lot
 3: debugging
An instance of Verbosity
is created when meep is imported, and is
accessible as meep.verbosity
. The meep.mpb
package also has a verbosity
flag in its C library, and it can also be managed via the Verbosity
class
after meep.mpb
is imported.
Note that this class is a Singleton, meaning that each call to create a new
Verbosity
actually gives you the same instance. The new C verbosity
flag will be added to a list of verbosity flags managed by this class.
The Verbosity
instance can be used as a global verbosity controller, and
assignments to any instance of Verbosity
will set the global verbosity
level for all library components. For example, this:
meep.verbosity(2)
# or meep.verbosity.set(2) if you prefer being more explicit
will set all of the managed verbosity flags to level 2.
Each managed verbosity flag can also be accessed individually if desired.
Each time a new C/C++ library verbosity flag is added to this Python class a
new property is added which can be used to access that individual flag.
Currently the properties that are available are named simply meep
and
mpb
. This means that you can set two different verbosity levels like this:
verbosity = meep.verbosity # not required, it's just to save some typing
verbosity.meep = 2
verbosity.mpb = 1
def __init__(self, cvar=None, name=None, initial_level=1):
See add_verbosity_var()
def __call__(self, level):
Convenience for setting the verbosity level. This lets you set the
global level by calling the instance like a function. For example, if
verbosity
is an instance of this class, then it's value can be changed
like this:
verbosity(0)
def add_verbosity_var(self, cvar=None, name=None, initial_level=1):
Add a new verbosity flag to be managed. cvar
should be some object
that has a verbosity
attribute, such as meep.cvar
or mpb.cvar
.
def get(self):
Returns the current global verbosity level.
def set(self, level):
Validates the range, and sets the global verbosity level. Returns the former value.
Miscellaneous Functions Reference
def quiet(quietval=True):
Meep ordinarily prints various diagnostic and progress information to standard output.
This output can be suppressed by calling this function with True
(the default). The
output can be enabled again by passing False
. This sets a global variable, so the
value will persist across runs within the same script.
This function is deprecated, please use the Verbosity class instead.
def interpolate(n, nums):
Given a list of numbers or Vector3
s nums
, linearly interpolates between them to
add n
new evenlyspaced values between each pair of consecutive values in the
original list.
Flux functions
def get_flux_freqs(f):
Given a flux object, returns a list of the frequencies that it is computing the spectrum for.
def get_fluxes(f):
Given a flux object, returns a list of the current flux spectrum that it has accumulated.
def scale_flux_fields(s, flux):
Scale the Fouriertransformed 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
.
def get_eigenmode_freqs(f):
Given a flux object, returns a list of the frequencies that it is computing the spectrum for.
Energy Functions
def get_energy_freqs(f):
Given an energy object, returns a list of the frequencies that it is computing the spectrum for.
def get_electric_energy(f):
Given an energy object, returns a list of the current energy density spectrum for the electric fields that it has accumulated.
def get_magnetic_energy(f):
Given an energy object, returns a list of the current energy density spectrum for the magnetic fields that it has accumulated.
def get_total_energy(f):
Given an energy object, returns a list of the current energy density spectrum for the total fields that it has accumulated.
Force Functions
def get_force_freqs(f):
Given a force object, returns a list of the frequencies that it is computing the spectrum for.
def get_forces(f):
Given a force object, returns a list of the current force spectrum that it has accumulated.
LDOS Functions
def Ldos(*args):
def Ldos(fcen, df, nfreq, freq):
Create an LDOS object with either frequency bandwidth df
centered at fcen
and
nfreq
equally spaced frequency points or an array/list freq
for arbitrarily spaced
frequencies. This can be passed to the dft_ldos
step function below as a keyword
argument.
def get_ldos_freqs(l):
Given an LDOS object, returns a list of the frequencies that it is computing the spectrum for.
def dft_ldos(*args, **kwargs):
def dft_ldos(fcen=None, df=None, nfreq=None, freq=None, ldos=None):
Compute the power spectrum of the sources (usually a single point dipole source),
normalized to correspond to the LDOS, in either a frequency bandwidth df
centered at
fcen
and nfreq
equally spaced frequency points or an array/list freq
for
arbitrarily spaced frequencies. One can also pass in an Ldos
object as
dft_ldos(ldos=my_Ldos)
.
The resulting spectrum is outputted as commadelimited text, prefixed by ldos:,
, and
is also stored in the ldos_data
variable of the Simulation
object after the run
is complete.
Near2Far Functions
def get_near2far_freqs(f):
Given a near2far
object, returns a list of the frequencies that it is computing the
spectrum for.
def scale_near2far_fields(s, near2far):
Scale the Fouriertransformed fields in near2far
by the complex number s
. e.g.
load_minus_near2far
is equivalent to load_near2far
followed by
scale_near2far_fields
with s=1
.
GDSII Functions
def GDSII_layers(fname):
Returns a list of integervalued layer indices for the layers present in the specified GDSII file.
mp.GDSII_layers('python/examples/coupler.gds')
Out[2]: [0, 1, 2, 3, 4, 5, 31, 32]
def GDSII_prisms(material, fname, layer=1, zmin=0.0, zmax=0.0):
Returns a list of GeometricObject
s with material
(mp.Medium
) on layer number
layer
of a GDSII file fname
with zmin
and zmax
(default 0).
def GDSII_vol(fname, layer, zmin, zmax):
Returns a mp.Volume
read from a GDSII file fname
on layer number layer
with
zmin
and zmax
(default 0). This function is useful for creating a FluxRegion
from a GDSII file as follows:
fr = mp.FluxRegion(volume=mp.GDSII_vol(fname, layer, zmin, zmax))
Run and Step Functions
def stop_when_fields_decayed(dt, c, pt, decay_by):
Return a condition
function, suitable for passing to Simulation.run
as the until
or until_after_sources
parameter, 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 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.
def stop_after_walltime(t):
Return a condition
function, suitable for passing to Simulation.run
as the until
parameter. Stops the simulation after t
seconds of wall time have passed.
def stop_on_interrupt():
Return a condition
function, suitable for passing to Simulation.run
as the until
parameter. Instead of terminating when receiving a SIGINT or SIGTERM signal from the
system, the simulation will abort time stepping and continue executing any code that
follows the run
function (e.g., outputting fields).
Output Functions
def output_epsilon(sim=None, *step_func_args, **kwargs):
Given a frequency frequency
, (provided as a keyword argument) output (relative
permittivity); for an anisotropic tensor the output is the harmonic
mean of the eigenvalues. If
frequency
is nonzero, the output is complex; otherwise it is the real,
frequencyindependent part of (the limit).
When called as part of a step function,
the sim
argument specifying the Simulation
object can be omitted, e.g.,
sim.run(mp.at_beginning(mp.output_epsilon(frequency=1/0.7)),until=10)
.
def output_mu(sim=None, *step_func_args, **kwargs):
Given a frequency frequency
, (provided as a keyword argument) output (relative
permeability); for an anisotropic tensor the output is the harmonic
mean of the eigenvalues. If
frequency
is nonzero, the output is complex; otherwise it is the real,
frequencyindependent part of (the limit).
When called as part of a step function,
the sim
argument specifying the Simulation
object can be omitted, e.g.,
sim.run(mp.at_beginning(mp.output_mu(frequency=1/0.7)),until=10)
.
def output_poynting(sim):
Output the Poynting flux . Note that you
might want to wrap this step function in synchronized_magnetic
to compute it more
accurately. See Synchronizing the Magnetic and Electric
Fields.
def output_hpwr(sim):
Output the magneticfield energy density
def output_dpwr(sim):
Output the electricfield energy density
def output_tot_pwr(sim):
Output the total electric and magnetic energy density. Note that you might want to
wrap this step function in synchronized_magnetic
to compute it more accurately. See
Synchronizing the Magnetic and Electric
Fields.
def output_png(compnt, options, rm_h5=True):
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 Tutorial/Basics/Output Tips and
Tricks.
It is often useful to use the h5topng C
or A
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
output_epsilon
. To make this easier, a builtin shell variable $EPS
is provided
which refers to the lastoutput dielectricfunction .h5
file. So, for example
output_png(mp.Ez,"C $EPS")
will output the field and overlay the dielectric
contours.
By default, output_png
deletes the .h5
file when it is done. To preserve the .h5
file requires output_png(component, h5topng_options, rm_h5=False)
.
def output_hfield(sim):
Outputs all the components of the field h, (magnetic) to an HDF5 file. That is, the different components are stored as different datasets within the same file.
def output_hfield_x(sim):
Output the component of the field h (magnetic). 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.
def output_hfield_y(sim):
Output the component of the field h (magnetic). 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.
def output_hfield_z(sim):
Output the component of the field h (magnetic). 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.
def output_hfield_r(sim):
Output the component of the field h (magnetic). 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.
def output_hfield_p(sim):
Output the component of the field h (magnetic). 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.
def output_bfield(sim):
Outputs all the components of the field b, (magnetic) to an HDF5 file. That is, the different components are stored as different datasets within the same file.
def output_bfield_x(sim):
Output the component of the field b (magnetic). 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.
def output_bfield_y(sim):
Output the component of the field b (magnetic). 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.
def output_bfield_z(sim):
Output the component of the field b (magnetic). 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.
def output_bfield_r(sim):
Output the component of the field b (magnetic). 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.
def output_bfield_p(sim):
Output the component of the field b (magnetic). 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. Note that for outputting the Poynting flux, you
might want to wrap the step function in synchronized_magnetic
to compute it more
accurately. See Synchronizing the Magnetic and Electric
Fields.
def output_efield(sim):
Outputs all the components of the field e, (electric) to an HDF5 file. That is, the different components are stored as different datasets within the same file.
def output_efield_x(sim):
Output the component of the field e (electric). 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.
def output_efield_y(sim):
Output the component of the field e (electric). 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.
def output_efield_z(sim):
Output the component of the field e (electric). 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.
def output_efield_r(sim):
Output the component of the field e (electric). 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.
def output_efield_p(sim):
Output the component of the field e (electric). 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. Note that for outputting the Poynting flux, you
might want to wrap the step function in synchronized_magnetic
to compute it more
accurately. See Synchronizing the Magnetic and Electric
Fields.
def output_dfield(sim):
Outputs all the components of the field d, (displacement) to an HDF5 file. That is, the different components are stored as different datasets within the same file.
def output_dfield_x(sim):
Output the component of the field d (displacement). 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.
def output_dfield_y(sim):
Output the component of the field d (displacement). 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.
def output_dfield_z(sim):
Output the component of the field d (displacement). 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.
def output_dfield_r(sim):
Output the component of the field d (displacement). 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.
def output_dfield_p(sim):
Output the component of the field d (displacement). 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. Note that for outputting the Poynting flux, you
might want to wrap the step function in synchronized_magnetic
to compute it more
accurately. See Synchronizing the Magnetic and Electric
Fields.
def output_sfield(sim):
Outputs all the components of the field s, (poynting flux) to an HDF5 file. That
is, the different components are stored as different datasets within the same file.
Note that you might want to wrap this step function in synchronized_magnetic
to
compute it more accurately. See Synchronizing the Magnetic and Electric
Fields.
def output_sfield_x(sim):
Output the component of the field s (poynting flux). 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.
def output_sfield_y(sim):
Output the component of the field s (poynting flux). 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.
def output_sfield_z(sim):
Output the component of the field s (poynting flux). 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.
def output_sfield_r(sim):
Output the component of the field s (poynting flux). 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.
def output_sfield_p(sim):
Output the component of the field s (poynting flux). 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. Note that for outputting the Poynting flux, you
might want to wrap the step function in synchronized_magnetic
to compute it more
accurately. See Synchronizing the Magnetic and Electric
Fields.