Meep is a versatile, free, and open-source software designed for electromagnetics simulations using the finite-difference time-domain (FDTD) method. Despite its extensive range of applications, broad features, and comprehensive documentation, new users might find it quite technical and challenging to learn.
In this beginner’s guide, we’ll recreate the famous double slit experiment in two dimensions, explaining the code at each step.
First we must quickly setup a Python environment from which we can import and utilize Meep. While many methods exist for installing Meep, here we use Conda, which is probably the simplest way. For alternative installation methods, please refer to Meep’s official documentation.
To start, download and install Miniconda and choose an installation directory, for instance, ~/.miniconda3/
. During installation, you’ll be prompted to initialize Miniconda. It’s recommended to select “yes”. This action will modify your .bashrc
file (or similar), ensuring that Conda’s base
environment is activated each time you open a new terminal.
If prefer not to have the base
environment automatically activated, disable this feature issuing the following command from a fresh terminal:
conda config --set auto_activate_base false
You can manually activate the base
environment using:
conda activate
We’ll use a Conda configuration file to create the environment. Copy and paste the following content into a file called environment.yml
:
name: meepchannels:- conda-forge- defaultsdependencies:- python=3.11- matplotlib=3.7- numpy=1.25- opencv=4.7- pymeep=1.27.0- pymeep-extras=1.27.0
Then, to create the environment, run:
conda env create -f environment.yml
Finally, activate the newly created environment using:
conda activate meep
You are now ready to use Meep from Python.
Launch Python and import the packages we’ll be using. Let me briefly explain what each package is for. The central focus of this post is, of course, meep
, which will be performing all of the heavy lifting in our electromagnetics simulation. To manipulate and work with the data produced by these simulations, we’ll use numpy
. Since our simulations can generate substantial data, we leverage h5py
for efficient storage. We use pyplot
for data visualization and plotting. Finally, os
is used for handling various filesystem tasks.
import h5pyimport matplotlib.pyplot as pltimport meep as mpimport numpy as npimport os
For convenience — and to avoid cluttering up the more important code — let’s also define some plotting functions now. We’ll use them later.
from matplotlib.colors import LinearSegmentedColormap# Increase the default resolution for imagesplt.rcParams['figure.dpi'] = 300plt.rcParams['savefig.dpi'] = 300# Plot everything on a dark backgroundplt.style.use('dark_background')# Some custom colormapscmap_alpha = LinearSegmentedColormap.from_list('custom_alpha', [[1, 1, 1, 0], [1, 1, 1, 1]])cmap_blue = LinearSegmentedColormap.from_list('custom_blue', [[0, 0, 0], [0, 0.66, 1], [1, 1, 1]])def label_plot(ax, title=None, xlabel=None, ylabel=None, elapsed=None):"""Add a title and x/y labels to the plot."""if title:ax.set_title(title)elif elapsed is not None:ax.set_title(f'{elapsed:0.1f} fs')if xlabel is not False:ax.set_xlabel('x (μm)'if xlabel is None else xlabel)if ylabel is not False:ax.set_ylabel('y (μm)'if ylabel is None else ylabel)def plot_eps_data(eps_data, domain, ax=None, **kwargs):"""Plot the wall geometry (dielectric data) within the domain."""ax = ax or plt.gca()ax.imshow(eps_data.T, cmap=cmap_alpha, extent=domain, origin='lower')label_plot(ax, **kwargs)def plot_ez_data(ez_data, domain, ax=None, vmax=None, aspect=None, **kwargs):"""Plot the amplitude of the complex-valued electric fielddata within the domain."""ax = ax or plt.gca()ax.imshow(np.abs(ez_data.T),interpolation='spline36',cmap=cmap_blue,extent=domain,vmax=vmax,aspect=aspect,origin='lower',)label_plot(ax, **kwargs)def plot_pml(pml_thickness, domain, ax=None):ax = ax or plt.gca()x_start = domain[0] + pml_thicknessx_end = domain[1] - pml_thicknessy_start = domain[2] + pml_thicknessy_end = domain[3] - pml_thicknessrect = plt.Rectangle((x_start, y_start),x_end - x_start,y_end - y_start,fill=None,color='#fff',linestyle='dashed',)ax.add_patch(rect)
Before going any further with the code, there’s an essential high-level concept that should be well-understood. Meep’s functionality stems from Maxwell’s equations, which are scale-invariant. Consequently, Meep employs dimensionless units.
Practically, this flexibility allows you to define the spatial and temporal dimensions of your simulation in any way you prefer. However, it also necessitates a careful conversion of any dimension-specific values into “Meep units” before using them as arguments in Meep’s various classes.
For the purposes of this post, let’s define the dimensions as follows:
Meep defines the speed of light internally as unity. Consequently, the units of distance and time in Meep both equate to 1 μm. A simple way to think about it is that one Meep time unit is how long it takes for light to traverse one Meep unit of distance.
For converting time from fs to Meep units, defining the speed of light in μm/fs is beneficial. So,
Therefore, time \(t\) measured in fs is equivalent to approximately \(0.3t\) in Meep time units. This conversion will be useful later.
The following piece of code first defines the speed of light (SOL
) in units of μm per fs. Then, it sets up the 2-dimensional spatial domain for the simulation, specified in μm. The domain
list represents the extent of the domain, and the center
and cell_size
variables define the center point and size of the cell in this domain, respectively, using Meep’s Vector3
object. These parameters will be used in the subsequent steps of the simulation.
# The speed of light in μm/fsSOL = 299792458e-9# 2D spatial domain measured in μmdomain = [0, 30, -10, 10]center = mp.Vector3((domain[1] + domain[0]) / 2,(domain[3] + domain[2]) / 2,)cell_size = mp.Vector3(domain[1] - domain[0],domain[3] - domain[2],)
The next code segment defines a wall with two apertures (holes). The wall’s dimensions and position, along with the size of the apertures, are all defined. The wall is constituted by three blocks arranged vertically and is created using a material with an exceptionally high dielectric constant, effectively blocking and reflecting light instead of allowing it to pass through.
The blocks are aligned to create the double-slit configuration with two outer blocks representing the edges of the wall and an inner block separating the two apertures. The geometry of the wall is created using these block definitions.
Moreover, the code sets up perfectly matched layers (PML) for the simulation boundaries. The PML is a layer that surrounds the simulation cell and is designed to absorb all incident waves without any reflection, thereby mimicking an infinite space.
# Dimensions of wall with two apertureswall_position = 10wall_thickness = 0.5aperture_width = 1inner_wall_len = 4 # wall separating the aperturesouter_wall_len = (cell_size[1]- 2*aperture_width- inner_wall_len) / 2# Define a wall material with high dielectric constant,# effectively blocking light and reflecting it insteadmaterial = mp.Medium(epsilon=1e6)# Define the wall as an array of 3 blocks arranged verticallygeometry = [mp.Block(mp.Vector3(wall_thickness, outer_wall_len, mp.inf),center=mp.Vector3(wall_position - center.x,domain[3] - outer_wall_len / 2),material=material),mp.Block(mp.Vector3(wall_thickness, outer_wall_len, mp.inf),center=mp.Vector3(wall_position - center.x,domain[2] + outer_wall_len / 2),material=material),mp.Block(mp.Vector3(wall_thickness, inner_wall_len, mp.inf),center=mp.Vector3(wall_position - center.x, 0),material=material),]# Perfectly matched layer of thickness 1pml_thickness = 1pml_layers = [mp.PML(pml_thickness)]
To quickly visualize the wall geometry, initialize a simulation (without running it) and extract the dielectric data.
# Extract and visualize the dielectric data (wall geometry)sim = mp.Simulation(cell_size=cell_size, geometry=geometry, resolution=10)sim.init_sim()eps_data = sim.get_array(center=mp.Vector3(), size=cell_size, component=mp.Dielectric)ax = plt.gca()plot_pml(pml_thickness, domain, ax=ax)plot_eps_data(eps_data, domain, ax=ax)
Next we define a light source with properties similar to a laser. We choose a wavelength of 470 nanometers, which corresponds approximately to light blue visible light, and a beam width of 6 units. The frequency of the light source is calculated as the reciprocal of the wavelength, given that the speed of light in Meep units is one.
# Light wavelength, frequency, and beam widthsource_lambda = 0.47 # in μmsource_frequency = 1 / source_lambdasource_beam_width = 6
To approximate a laser, we want to design a continuous Gaussian beam that originates from the left-hand side and propagates towards the right. Ideally, this beam should be collimated, which means that the light rays are parallel to each other, resulting in a beam that doesn’t spread out and maintains its narrow shape over distance. Unfortunately, achieving this setup is a somewhat involved, as Meep does not offer a built-in source that fulfills these requirements directly. At first glance, the GaussianBeam2DSource
may seem like what we want, but it requires a focal point, thus it isn’t truly collimated.
To overcome this, we will combine Meep’s various source classes with our own amplitude function, which simply returns a complex current amplitude for any given position. Let’s start the process of building our light source by constructing this function.
First, we define a function that returns a complex-valued plane wave propagating in the \(x\)-direction. This function uses the np.exp
function to generate the complex exponential, which is representative of a wave. The plot we create here shows the real part of the wave amplitude as a function of position, illustrating how it propagates along the \(x\)-axis, regardless of \(y\). The wave amplitude varies sinusoidally with \(x\), as expected.
# A method to return a complex-valued plane wave in the x-directiondef plane_wave(x):return np.exp(2j * np.pi / source_lambda * x)# Plot the plane wavexarr = np.linspace(0, 10*source_lambda, 1000)wave = plane_wave(xarr)plt.plot(xarr, wave.real)plt.xlabel('x (μm)')plt.ylabel('Amplitude')plt.grid(True)
Next, we define another function to compute the Gaussian profile along the \(y\)-direction. This profile decreases from a maximum at \(y=0\) (the beam center) to near-zero values as we move away from the center in either direction. The width of the Gaussian profile is determined by source_beam_width / 2
, representing the standard deviation of the Gaussian distribution. This profile is essentially a cross-sectional distribution of the beam amplitude. When we plot the profile, we see that it has a peak at \(y=0\) and falls off on either side, illustrating how the intensity of the laser light decreases away from the center of the beam.
# A method to compute the Gaussian profile in the y-directiondef gaussian_profile(y):return np.exp(-y**2 / (2 * (source_beam_width / 2)**2))# Plot the Guassian profileyarr = np.linspace(domain[2], domain[3], 200)prof = gaussian_profile(yarr)plt.plot(yarr, prof)plt.xlabel('y (μm)')plt.ylabel('Amplitude')plt.grid(True)
Now we can combine the two functions defined above to generate a complex amplitude distribution for a 2D Gaussian beam. This is achieved by multiplying the plane wave function by the Gaussian profile function. The resulting function gives the complex amplitude of the light wave at any position \((x, y)\) in the domain. We can plot it as a color-map, showing the spatial variation of the beam intensity in the 2D plane.
# A 2D grid of pointsX, Y = np.meshgrid(xarr, yarr)# Plot the combined termscombined = plane_wave(X) * gaussian_profile(Y)plt.imshow(np.real(combined),cmap='coolwarm',aspect='auto',extent=[xarr[0], xarr[-1], yarr[0], yarr[-1]],origin='lower',)plt.xlabel('x (μm)')plt.ylabel('y (μm)')
We use the combined function from the previous step to define a light source in Meep. This is achieved using Meep’s Source
object, to which we provide a ContinuousSource
configured to our specific frequency. We specify the direction and type of the current component with mp.Ez
. Using center
and size
we position the source to span the entire length of the left-hand edge of the domain, leaving some extra room for the PML on the left. Importantly, the amp_func
parameter is given the combined function we defined above, which ensures that the light source has the desired spatial amplitude variation.
Since our light source extends into the PML regions on the top and bottom, it’s useful to set is_integrated=True
. This helps to correctly handle the source’s overlap with the PML, minimizing any undesirable interference effects at the boundary.
Lastly, a note about the frequency
supplied to the ContinuousSource
object and the source_lambda
supplied to the amp_func
method: while they are related (since wavelength and frequency are inversely related in light waves), they control different aspects of the source. The frequency
sets the temporal oscillation of the source, i.e., how fast it’s oscillating in time. On the other hand, source_lambda
in amp_func
helps set up the initial spatial phase distribution across the source plane. Both are necessary to properly define the source’s temporal and spatial characteristics.
def amp_func(pos):return plane_wave(pos[0]) * gaussian_profile(pos[1])source = mp.Source(src=mp.ContinuousSource(frequency=source_frequency,is_integrated=True,),component=mp.Ez,center= mp.Vector3(1, 0, 0) - center, # positioned far-left, excluding PMLsize=mp.Vector3(y=cell_size[1]), # span entire height, including PMLamp_func=amp_func,)
Now we can finally construct and run the simulation. We’ll start by defining an appropriate resolution
for the simulation.
In Meep, the resolution
parameter in the Simulation
class determines the spatial discretization of the simulation domain. Higher resolution means more grid points per unit length and thus more detailed representation of the geometry and fields, but it also means a larger problem size, leading to longer simulation time and greater memory usage.
Here are some guidelines to choose a suitable resolution:
Wavelength of the light: The grid size should be small compared to the smallest wavelength in the simulation. A good rule of thumb is to have at least 8-10 pixels per smallest wavelength. For example, if the smallest wavelength of interest is 0.5 μm, a resolution of 16-20 pixels/μm could be a good starting point.
Feature size of the structures: The resolution should be high enough to accurately describe the smallest geometric features of the structure in your simulation. If you have small, intricate details in your geometry, you’ll need a higher resolution to accurately capture these.
Simulation time and resources: The resolution you choose also depends on how much computation resources you have and how long you’re willing to wait. Higher resolution will increase both the computation time and memory requirements.
Accuracy vs efficiency trade-off: Start with a moderate resolution, run the simulation, and analyze the results. Then increase the resolution and repeat. If the results change significantly with increasing resolution, it might be an indication that the initial resolution was too low. Once the results converge (i.e., further increases in resolution result in negligible changes in the results), you’ve found a resolution that’s high enough for your needs.
For our simple scenario we’ll limit ourselves to the first two guidelines. First identify the smallest feature length or wavelength. Then assign 10 pixels to that length to determine the simulation resolution. Rounding up to the nearest integer ensures that the grid volume has an integer number of pixels.
# Define resolution in terms of smallest componentsmallest_length = min(source_lambda,wall_thickness,aperture_width,inner_wall_len,)pixel_count = 10resolution = int(np.ceil(pixel_count / smallest_length))print('Simulation resolution:', resolution)
Output:
Simulation resolution: 22
Now create a Simulation
object, pulling together the variables we defined earlier, such as cell size, light source, boundary layers, geometry, and resolution. Also specify force_complex_fields=True
, forcing Meep to use complex-valued fields rather than real-valued.
Once the setup is complete, run the simulation for a duration slightly exceeding the time it would take for light to travel across the full length of the domain. This additional time is necessary to allow any initial wavefront disturbances to exit the simulation area, thereby cleaning up the final result.
Finally, extract the data from the completed simulation, and use our plotting functions to visualize it. We can now see how the electric field interacts with the geometry objects.
# Simulation objectsim = mp.Simulation(cell_size=cell_size,sources=[source],boundary_layers=pml_layers,geometry=geometry,resolution=resolution,force_complex_fields=True,)# Convenience method to extract Ez and dielectric datadef get_data(sim, cell_size):ez_data = sim.get_array(center=mp.Vector3(), size=cell_size, component=mp.Ez)eps_data = sim.get_array(center=mp.Vector3(), size=cell_size, component=mp.Dielectric)return ez_data, eps_data# Run simulation until light travels full length (plus some)sim.run(until=cell_size[0] + 10)ez_data, eps_data = get_data(sim, cell_size)# Plot simulationax = plt.gca()plot_ez_data(ez_data, domain, ax=ax)plot_eps_data(eps_data, domain, ax=ax)plot_pml(pml_thickness, domain, ax=ax)
In the next code block, we organize and run a simulation over a defined period, given in fs, capturing periodic snapshots of the process. These frames provide a time-resolved picture of how the light travels through the simulated domain. Keep in mind that we must convert into Meep time units by multiplying by SOL
, ensuring that our simulation aligns with the dimensional units used by the Meep framework.
The data generated during this simulation is saved in an efficient binary format, ideal for handling large datasets. This allows us to execute the simulation in one go, capturing the entire process without interruption, while also optimizing our computational resources.
For convenience and cleaner code, we’re creating a reusable wrapper function, simulate
, which handles the entire operation from start to finish. It prepares the data storage, runs the simulation in intervals, and saves the electric field data for each snapshot.
# Sim duration and number of snapshotssim_time = 120 # in fsn_frames = 6# Where to save the resultssim_path = 'simulation.h5'# Simulation objectsim = mp.Simulation(cell_size=cell_size,sources=[source],boundary_layers=pml_layers,geometry=geometry,resolution=resolution,force_complex_fields=True,)def simulate(sim, sim_path, sim_time, n_frames):# Remove previous sim file, if anyif os.path.exists(sim_path):os.remove(sim_path)# Time delta (in fs) between snapshots. Note that# we subtract 1 because we include the initial state# as the first frame.fs_delta = sim_time / (n_frames - 1)# Save data to an HDF5 binary filewith h5py.File(sim_path, 'a') as f:# Save simulation params for future referencef.attrs['sim_time'] = sim_timef.attrs['n_frames'] = n_framesf.attrs['fs_delta'] = fs_deltaf.attrs['resolution'] = sim.resolution# Save initial state as first framesim.init_sim()ez_data, eps_data = get_data(sim, cell_size)f.create_dataset('ez_data',shape=(n_frames, *ez_data.shape),dtype=ez_data.dtype,)f.create_dataset('eps_data',shape=eps_data.shape,dtype=eps_data.dtype,)f['ez_data'][0] = ez_dataf['eps_data'][:] = eps_data# Simulate and capture remaining snapshotsfor i in range(1, n_frames):# Run until the next frame timesim.run(until=SOL * fs_delta)# Capture electral field dataez_data, _ = get_data(sim, cell_size)f['ez_data'][i] = ez_datasimulate(sim, sim_path, sim_time, n_frames)
When the simulation completes, visualize the snapshots side-by-side.
fig_rows = 3fig_cols = 2n_subplots = fig_rows * fig_colsfig, ax = plt.subplots(fig_rows,fig_cols,figsize=(9, 12),sharex=False,sharey=True,)with h5py.File(sim_path, 'r') as f:for k in range(n_subplots):i, j = int(k / fig_cols), (k % fig_cols)# i, j = (k % fig_rows), int(k / fig_rows)_ax = ax[i][j]ez_data = f['ez_data'][k]eps_data = f['eps_data'][:]elapsed = k * f.attrs['fs_delta']vmax = 0.6 # force consistent brightnessplot_ez_data(ez_data, domain, ax=_ax, vmax=vmax, elapsed=elapsed)plot_eps_data(eps_data, domain, ax=_ax)plot_pml(pml_thickness, domain, ax=_ax)
For the final part of the tutorial, we are focusing on the visualization of the interference pattern created when light exits the double slits. By doing so, we can better understand how the interference pattern evolves as the light propagates further away from the slits. Keep in mind that the intensity of a wave is proportional to the square of its amplitude
Utilizing the final snapshot from our simulation — which offers the most “complete” view of the propagated light and has fewer residual ripples — we compute the light intensity at various slices along the \(x\)-axis, following the light’s exit from the double slits. Each slice provides us with an instantaneous snapshot of the light intensity at that distance, illustrating the evolving interference pattern.
# Grab final simulation snapshot without time-averagingwith h5py.File(sim_path, 'r') as f:final_snap = f['ez_data'][-1]# Compute intensity as square of the complex amplitudefinal_snap = np.abs(final_snap)**2# Pick slices at different distances from the double slitslice_dists = [10, 11, 12, 15, 20, 25]slices = [final_snap[i * resolution] for i in slice_dists]yarr = np.linspace(domain[2], domain[3], final_snap.shape[1])# A rather involved plotting functiondef plot_intensity(slice, yarr, ax1, ax2, vmax=None, xval=None, xlabel=False, ylabel=False):ax1.plot(yarr, slice)ax1.tick_params(axis='x', labelbottom=False)if ylabel:ax1.set_ylabel('$|E|^2$')else:ax1.tick_params('y', labelleft=False)if xval:ax1.annotate(f'x={xval}',xy=(1, 1),xytext=(-4, -4),xycoords='axes fraction',textcoords='offset pixels',horizontalalignment='right',verticalalignment='top',)ax2.imshow(np.vstack(slice).T,cmap='inferno',aspect='auto',vmax=vmax,extent=[yarr[0], yarr[-1], 0, 1],)ax2.set_xlim([yarr[0], yarr[-1]])ax2.tick_params('y', labelleft=False)ax2.set_yticks([])if xlabel:ax2.set_xlabel('y (μm)')else:ax2.tick_params(axis='x', labelbottom=False)fig, ax = plt.subplots(4, 3,figsize=(9, 6),gridspec_kw=dict(width_ratios=(4, 4, 4),height_ratios=(4, 1, 4, 1),wspace=0.12,hspace=0.1,),sharex='col',sharey='row',)for k, slice in enumerate(slices):i = 2 * int(k / 3)j = k % 3plot_intensity(slice, yarr, ax[i][j], ax[i+1][j],vmax=np.max(slices[:3]) if k < 3 else np.max(slices[3:]),xval=slice_dists[k],xlabel=(i==2),ylabel=(j==0))
That concludes the tutorial — I hope it has proven useful to you! The aim was to create an accessible, beginner’s guide to electromagnetics simulation using Meep, a free, open-source software built on the finite-difference time-domain (FDTD) method. We began by setting up the Python environment and importing the required packages. Then, we delved into the heart of the tutorial, recreating the iconic double-slit experiment. We covered key concepts such as Meep units, domain and geometry setup, light source creation, and simulation execution, using sample code for better understanding. We finished by exploring interference patterns, thereby showcasing Meep’s capabilities in visualizing and comprehending complex electromagnetics phenomena.
GET IN TOUCH
Have a project in mind?
Reach out directly to hello@humaticlabs.com or use the contact form.