From c18fa7a4776fd9d8485861b15f8153a73a1ee78e Mon Sep 17 00:00:00 2001
From: Ardavan Oskooi <ardavano@google.com>
Date: Fri, 12 Nov 2021 15:58:39 -0800
Subject: [PATCH 1/8] plot geometry without initializing structure class

---
 python/simulation.py               |  32 ++--
 python/tests/test_visualization.py |   2 +-
 python/visualization.py            | 231 +++++++++++++++--------------
 3 files changed, 141 insertions(+), 124 deletions(-)

diff --git a/python/simulation.py b/python/simulation.py
index 83742fc71..be93866cf 100644
--- a/python/simulation.py
+++ b/python/simulation.py
@@ -2755,8 +2755,8 @@ def get_farfield(self, near2far, x):
         (Fourier-transformed) "far" fields at `x` as list of length 6`nfreq`, 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,...)$
         in Cartesian coordinates and
-        $(E_r^1,E_\phi^1,E_z^1,H_r^1,H_\phi^1,H_z^1,E_r^2,E_\phi^2,E_z^2,H_r^2,H_\phi^2,H_z^2,...)$
-        in cylindrical coordinates for the frequencies 1,2,…,`nfreq`.
+        $(E_r^1,E_\\phi^1,E_z^1,H_r^1,H_\\phi^1,H_z^1,E_r^2,E_\\phi^2,E_z^2,H_r^2,H_\\phi^2,H_z^2,...)$
+        in cylindrical coordinates for the frequencies 1,2,...,`nfreq`.
         """
         return mp._get_farfield(near2far.swigobj, py_v3_to_vec(self.dimensions, x, is_cylindrical=self.is_cylindrical))
 
@@ -4149,20 +4149,22 @@ def plot2D(self, ax=None, output_plane=None, fields=None, labels=False,
             - `alpha=0.6`: transparency of fields
             - `post_process=np.real`: post processing function to apply to fields (must be
               a function object)
-        * `frequency`: for materials with a [frequency-dependent
-          permittivity](Materials.md#material-dispersion) $\\varepsilon(f)$, specifies the
-          frequency $f$ (in Meep units) of the real part of the permittivity to use in the
-          plot. Defaults to the `frequency` parameter of the [Source](#source) object.
-
         """
-        return vis.plot2D(self, ax=ax, output_plane=output_plane, fields=fields, labels=labels,
-            eps_parameters=eps_parameters, boundary_parameters=boundary_parameters,
-            source_parameters=source_parameters,monitor_parameters=monitor_parameters,
-            field_parameters=field_parameters, frequency=frequency,
-            plot_eps_flag=plot_eps_flag, plot_sources_flag=plot_sources_flag,
-            plot_monitors_flag=plot_monitors_flag, plot_boundaries_flag=plot_boundaries_flag,
-            **kwargs)
-
+        return vis.plot2D(self,
+                          ax=ax,
+                          output_plane=output_plane,
+                          fields=fields,
+                          labels=labels,
+                          eps_parameters=eps_parameters,
+                          boundary_parameters=boundary_parameters,
+                          source_parameters=source_parameters,
+                          monitor_parameters=monitor_parameters,
+                          field_parameters=field_parameters,
+                          plot_eps_flag=plot_eps_flag,
+                          plot_sources_flag=plot_sources_flag,
+                          plot_monitors_flag=plot_monitors_flag,
+                          plot_boundaries_flag=plot_boundaries_flag,
+                          **kwargs)
 
     def plot_fields(self,**kwargs):
         return vis.plot_fields(self,**kwargs)
diff --git a/python/tests/test_visualization.py b/python/tests/test_visualization.py
index e2b9d2846..1b030716e 100644
--- a/python/tests/test_visualization.py
+++ b/python/tests/test_visualization.py
@@ -161,7 +161,7 @@ def test_animation_output(self):
         sim = setup_sim() # generate 2D simulation
 
         Animate = mp.Animate2D(sim,fields=mp.Ez, realtime=False, normalize=False) # Check without normalization
-        Animate_norm = mp.Animate2D(sim,mp.Ez,realtime=False,normalize=True) # Check with normalization
+        Animate_norm = mp.Animate2D(sim, mp.Ez, realtime=False, normalize=True) # Check with normalization
 
         # test both animation objects during same run
         sim.run(
diff --git a/python/visualization.py b/python/visualization.py
index b42a3dcf7..e2bf3de62 100644
--- a/python/visualization.py
+++ b/python/visualization.py
@@ -95,7 +95,7 @@ def filter_dict(dict_to_filter, func_with_kwargs):
 # ------------------------------------------------------- #
 # Routines to add legends to plot
 
-def place_label(ax,label_text,x,y,centerx,centery,label_parameters=None):
+def place_label(ax, label_text, x, y, centerx, centery, label_parameters=None):
 
     label_parameters = default_label_parameters if label_parameters is None else dict(default_label_parameters, **label_parameters)
 
@@ -112,29 +112,37 @@ def place_label(ax,label_text,x,y,centerx,centery,label_parameters=None):
     else:
         ytext = offset
 
-    ax.annotate(label_text, xy=(x, y), xytext=(xtext,ytext),
-            textcoords='offset points', ha='center', va='bottom',
-            bbox=dict(boxstyle='round,pad=0.2', fc=color, alpha=alpha),
-            arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=0.5',
-                            color=color))
+    ax.annotate(label_text, xy=(x, y), xytext=(xtext, ytext),
+                textcoords='offset points', ha='center', va='bottom',
+                bbox=dict(boxstyle='round,pad=0.2', fc=color, alpha=alpha),
+                arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=0.5',
+                                color=color))
     return ax
 
 # ------------------------------------------------------- #
 # Helper functions used to plot volumes on a 2D plane
 
-# Returns the intersection points of 2 Volumes.
+# Returns the intersection points of two Volumes.
 # Volumes must be a line, plane, or rectangular prism
 # (since they are volume objects)
-def intersect_volume_volume(volume1,volume2):
+def intersect_volume_volume(volume1, volume2):
     # volume1 ............... [volume]
     # volume2 ............... [volume]
 
     # Represent the volumes by an "upper" and "lower" coordinate
-    U1 = [volume1.center.x+volume1.size.x/2,volume1.center.y+volume1.size.y/2,volume1.center.z+volume1.size.z/2]
-    L1 = [volume1.center.x-volume1.size.x/2,volume1.center.y-volume1.size.y/2,volume1.center.z-volume1.size.z/2]
-
-    U2 = [volume2.center.x+volume2.size.x/2,volume2.center.y+volume2.size.y/2,volume2.center.z+volume2.size.z/2]
-    L2 = [volume2.center.x-volume2.size.x/2,volume2.center.y-volume2.size.y/2,volume2.center.z-volume2.size.z/2]
+    U1 = [volume1.center.x+volume1.size.x/2,
+          volume1.center.y+volume1.size.y/2,
+          volume1.center.z+volume1.size.z/2]
+    L1 = [volume1.center.x-volume1.size.x/2,
+          volume1.center.y-volume1.size.y/2,
+          volume1.center.z-volume1.size.z/2]
+
+    U2 = [volume2.center.x+volume2.size.x/2,
+          volume2.center.y+volume2.size.y/2,
+          volume2.center.z+volume2.size.z/2]
+    L2 = [volume2.center.x-volume2.size.x/2,
+          volume2.center.y-volume2.size.y/2,
+          volume2.center.z-volume2.size.z/2]
 
     # Evaluate intersection
     U = np.min([U1,U2],axis=0)
@@ -164,13 +172,13 @@ def intersect_volume_volume(volume1,volume2):
 # All of the 2D plotting routines need an output plane over which to plot.
 # The user has many options to specify this output plane. They can pass
 # the output_plane parameter, which is a 2D volume object. They can specify
-# a volume using in_volume, which stores the volume as a C volume, not a python
+# a volume using in_volume, which stores the volume as a C volume, not a Python
 # volume. They can also do nothing and plot the XY plane through Z=0.
 #
 # Not only do we need to check for all of these possibilities, but we also need
 # to check if the user accidentally specifies a plane that stretches beyond the
 # simulation domain.
-def get_2D_dimensions(sim,output_plane):
+def get_2D_dimensions(sim, output_plane):
     from meep.simulation import Volume
 
     # Pull correct plane from user
@@ -182,12 +190,12 @@ def get_2D_dimensions(sim,output_plane):
         plane_center, plane_size = (sim.geometry_center, sim.cell_size)
     plane_volume = Volume(center=plane_center,size=plane_size)
 
-    if plane_size.x!=0 and plane_size.y!=0 and plane_size.z!=0:
+    if plane_size.x != 0 and plane_size.y != 0 and plane_size.z != 0:
         raise ValueError("Plane volume must be 2D (a plane).")
 
-    check_volume = Volume(center=sim.geometry_center,size=sim.cell_size)
+    check_volume = Volume(center=sim.geometry_center, size=sim.cell_size)
 
-    vertices = intersect_volume_volume(check_volume,plane_volume)
+    vertices = intersect_volume_volume(check_volume, plane_volume)
 
     if len(vertices) == 0:
         raise ValueError("The specified user volume is completely outside of the simulation domain.")
@@ -203,10 +211,7 @@ def get_2D_dimensions(sim,output_plane):
 # ------------------------------------------------------- #
 # actual plotting routines
 
-def plot_volume(sim,ax,volume,output_plane=None,plotting_parameters=None,label=None):
-    if not sim._is_initialized:
-        sim.init_sim()
-
+def plot_volume(sim, ax, volume, output_plane=None, plotting_parameters=None, label=None):
     import matplotlib.patches as patches
     from matplotlib import pyplot as plt
     from meep.simulation import Volume
@@ -215,39 +220,57 @@ def plot_volume(sim,ax,volume,output_plane=None,plotting_parameters=None,label=N
     plotting_parameters = default_volume_parameters if plotting_parameters is None else dict(default_volume_parameters, **plotting_parameters)
 
     # Get domain measurements
-    sim_center, sim_size = get_2D_dimensions(sim,output_plane)
+    sim_center, sim_size = get_2D_dimensions(sim, output_plane)
 
-    plane = Volume(center=sim_center,size=sim_size)
+    plane = Volume(center=sim_center, size=sim_size)
 
     # Pull volume parameters
     size = volume.size
     center = volume.center
 
-    xmax = center.x+size.x/2
-    xmin = center.x-size.x/2
-    ymax = center.y+size.y/2
-    ymin = center.y-size.y/2
-    zmax = center.z+size.z/2
-    zmin = center.z-size.z/2
+    xmax = center.x + size.x/2
+    xmin = center.x - size.x/2
+    ymax = center.y + size.y/2
+    ymin = center.y - size.y/2
+    zmax = center.z + size.z/2
+    zmin = center.z - size.z/2
 
     # Add labels if requested
     if label is not None and mp.am_master():
         if sim_size.x == 0:
-            ax = place_label(ax,label,center.y,center.z,sim_center.y,sim_center.z,label_parameters=plotting_parameters)
+            ax = place_label(ax,
+                             label,
+                             center.y,
+                             center.z,
+                             sim_center.y,
+                             sim_center.z,
+                             label_parameters=plotting_parameters)
         elif sim_size.y == 0:
-            ax = place_label(ax,label,center.x,center.z,sim_center.x,sim_center.z,label_parameters=plotting_parameters)
+            ax = place_label(ax,
+                             label,
+                             center.x,
+                             center.z,
+                             sim_center.x,
+                             sim_center.z,
+                             label_parameters=plotting_parameters)
         elif sim_size.z == 0:
-            ax = place_label(ax,label,center.x,center.y,sim_center.x,sim_center.y,label_parameters=plotting_parameters)
+            ax = place_label(ax,
+                             label,
+                             center.x,
+                             center.y,
+                             sim_center.x,
+                             sim_center.y,
+                             label_parameters=plotting_parameters)
 
     # Intersect plane with volume
-    intersection = intersect_volume_volume(volume,plane)
+    intersection = intersect_volume_volume(volume, plane)
 
     # Sort the points in a counter clockwise manner to ensure convex polygon is formed
     def sort_points(xy):
         xy = np.squeeze(xy)
-        xy_mean = np.mean(xy,axis=0)
-        theta = np.arctan2(xy[:,1]-xy_mean[1],xy[:,0]-xy_mean[0])
-        return xy[np.argsort(theta,axis=0),:]
+        xy_mean = np.mean(xy, axis=0)
+        theta = np.arctan2(xy[:,1] - xy_mean[1], xy[:,0] - xy_mean[0])
+        return xy[np.argsort(theta, axis=0), :]
 
     if mp.am_master():
         # Point volume
@@ -269,16 +292,16 @@ def sort_points(xy):
         elif len(intersection) == 2:
             line_args = {key:value for key, value in plotting_parameters.items() if key in ['color','linestyle','linewidth','alpha']}
             # Plot YZ
-            if sim_size.x==0:
-                ax.plot([a.y for a in intersection],[a.z for a in intersection], **line_args)
+            if sim_size.x == 0:
+                ax.plot([a.y for a in intersection], [a.z for a in intersection], **line_args)
                 return ax
-            #Plot XZ
-            elif sim_size.y==0:
-                ax.plot([a.x for a in intersection],[a.z for a in intersection], **line_args)
+            # Plot XZ
+            elif sim_size.y == 0:
+                ax.plot([a.x for a in intersection], [a.z for a in intersection], **line_args)
                 return ax
             # Plot XY
-            elif sim_size.z==0:
-                ax.plot([a.x for a in intersection],[a.y for a in intersection], **line_args)
+            elif sim_size.z == 0:
+                ax.plot([a.x for a in intersection], [a.y for a in intersection], **line_args)
                 return ax
             else:
                 return ax
@@ -287,15 +310,15 @@ def sort_points(xy):
         elif len(intersection) > 2:
             planar_args = {key:value for key, value in plotting_parameters.items() if key in ['edgecolor','linewidth','facecolor','hatch','alpha']}
             # Plot YZ
-            if sim_size.x==0:
+            if sim_size.x == 0:
                 ax.add_patch(patches.Polygon(sort_points([[a.y,a.z] for a in intersection]), **planar_args))
                 return ax
-            #Plot XZ
+            # Plot XZ
             elif sim_size.y==0:
                 ax.add_patch(patches.Polygon(sort_points([[a.x,a.z] for a in intersection]), **planar_args))
                 return ax
             # Plot XY
-            elif sim_size.z==0:
+            elif sim_size.z == 0:
                 ax.add_patch(patches.Polygon(sort_points([[a.x,a.y] for a in intersection]), **planar_args))
                 return ax
             else:
@@ -304,16 +327,12 @@ def sort_points(xy):
             return ax
     return ax
 
-def plot_eps(sim,ax,output_plane=None,eps_parameters=None,frequency=0):
-    if sim.structure is None:
-        sim.init_sim()
-
-
+def plot_eps(sim, ax, output_plane=None, eps_parameters=None, resolution=0):
     # consolidate plotting parameters
     eps_parameters = default_eps_parameters if eps_parameters is None else dict(default_eps_parameters, **eps_parameters)
 
     # Get domain measurements
-    sim_center, sim_size = get_2D_dimensions(sim,output_plane)
+    sim_center, sim_size = get_2D_dimensions(sim, output_plane)
 
     xmin = sim_center.x - sim_size.x/2
     xmax = sim_center.x + sim_size.x/2
@@ -322,28 +341,40 @@ def plot_eps(sim,ax,output_plane=None,eps_parameters=None,frequency=0):
     zmin = sim_center.z - sim_size.z/2
     zmax = sim_center.z + sim_size.z/2
 
-    center = Vector3(sim_center.x,sim_center.y,sim_center.z)
-    cell_size = Vector3(sim_size.x,sim_size.y,sim_size.z)
+    center = Vector3(sim_center.x, sim_center.y, sim_center.z)
+    cell_size = Vector3(sim_size.x, sim_size.y, sim_size.z)
+
+    grid_resolution = resolution if resolution else sim.resolution + 1
 
     if sim_size.x == 0:
         # Plot y on x axis, z on y axis (YZ plane)
         extent = [ymin,ymax,zmin,zmax]
         xlabel = 'Y'
         ylabel = 'Z'
+        xtics = np.array([0])
+        ytics = np.linspace(ymin,ymax, grid_resolution)
+        ztics = np.linspace(zmin,zmax, grid_resolution)
     elif sim_size.y == 0:
         # Plot x on x axis, z on y axis (XZ plane)
         extent = [xmin,xmax,zmin,zmax]
         xlabel = 'X'
         ylabel = 'Z'
+        xtics = np.linspace(xmin, xmax, grid_resolution)
+        ytics = np.array([0])
+        ztics = np.linspace(zmin, zmax, grid_resolution)
     elif sim_size.z == 0:
         # Plot x on x axis, y on y axis (XY plane)
         extent = [xmin,xmax,ymin,ymax]
         xlabel = 'X'
         ylabel = 'Y'
+        xtics = np.linspace(xmin, xmax, grid_resolution)
+        ytics = np.linspace(ymin, ymax, grid_resolution)
+        ztics = np.array([0])
     else:
         raise ValueError("A 2D plane has not been specified...")
 
-    eps_data = np.rot90(np.real(sim.get_array(center=center, size=cell_size, component=mp.Dielectric, frequency=frequency)))
+    eps_data = sim.get_epsilon_grid(xtics, ytics, ztics)
+
     if mp.am_master():
         if eps_parameters['contour']:
             ax.contour(eps_data, 0, colors='black', origin='upper', extent=extent, linewidths=eps_parameters['contour_linewidth'])
@@ -354,14 +385,11 @@ def plot_eps(sim,ax,output_plane=None,eps_parameters=None,frequency=0):
 
     return ax
 
-def plot_boundaries(sim,ax,output_plane=None,boundary_parameters=None):
-    if not sim._is_initialized:
-        sim.init_sim()
-
+def plot_boundaries(sim, ax, output_plane=None, boundary_parameters=None):
     # consolidate plotting parameters
     boundary_parameters = default_boundary_parameters if boundary_parameters is None else dict(default_boundary_parameters, **boundary_parameters)
 
-    def get_boundary_volumes(thickness,direction,side):
+    def get_boundary_volumes(thickness, direction, side):
         from meep.simulation import Volume
 
         thickness = boundary.thickness
@@ -427,10 +455,7 @@ def get_boundary_volumes(thickness,direction,side):
             ax = plot_volume(sim,ax,vol,output_plane,plotting_parameters=boundary_parameters)
     return ax
 
-def plot_sources(sim,ax,output_plane=None,labels=False,source_parameters=None):
-    if not sim._is_initialized:
-        sim.init_sim()
-
+def plot_sources(sim, ax, output_plane=None, labels=False, source_parameters=None):
     from meep.simulation import Volume
 
     # consolidate plotting parameters
@@ -443,13 +468,10 @@ def plot_sources(sim,ax,output_plane=None,labels=False,source_parameters=None):
         ax = plot_volume(sim,ax,vol,output_plane,plotting_parameters=source_parameters,label=label)
     return ax
 
-def plot_monitors(sim,ax,output_plane=None,labels=False,monitor_parameters=None):
-    if not sim._is_initialized:
-        sim.init_sim()
-
+def plot_monitors(sim, ax, output_plane=None, labels=False, monitor_parameters=None):
     from meep.simulation import Volume
 
-     # consolidate plotting parameters
+    # consolidate plotting parameters
     monitor_parameters = default_monitor_parameters if monitor_parameters is None else dict(default_monitor_parameters, **monitor_parameters)
 
     label = 'monitor' if labels else None
@@ -460,7 +482,7 @@ def plot_monitors(sim,ax,output_plane=None,labels=False,monitor_parameters=None)
             ax = plot_volume(sim,ax,vol,output_plane,plotting_parameters=monitor_parameters,label=label)
     return ax
 
-def plot_fields(sim,ax=None,fields=None,output_plane=None,field_parameters=None):
+def plot_fields(sim, ax=None, fields=None, output_plane=None, field_parameters=None):
     if not sim._is_initialized:
         sim.init_sim()
 
@@ -515,36 +537,17 @@ def plot_fields(sim,ax=None,fields=None,output_plane=None,field_parameters=None)
         return np.rot90(fields)
     return ax
 
-def plot2D(sim,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,
+def plot2D(sim, ax=None, output_plane=None, fields=None, labels=False,
+           eps_parameters=None, boundary_parameters=None,
+           source_parameters=None, monitor_parameters=None,
+           field_parameters=None, resolution=None,
            plot_eps_flag=True, plot_sources_flag=True,
            plot_monitors_flag=True, plot_boundaries_flag=True):
 
-    # Initialize the simulation
-    if sim.structure is None:
-        sim.init_sim()
     # Ensure a figure axis exists
     if ax is None and mp.am_master():
         from matplotlib import pyplot as plt
         ax = plt.gca()
-    # Determine a frequency to plot all epsilon
-    if frequency is None:
-        try:
-            # Continuous sources
-            frequency = sim.sources[0].frequency
-        except:
-            try:
-                # Gaussian sources
-                frequency = sim.sources[0].src.frequency
-            except:
-                try:
-                    # Custom sources
-                    frequency = sim.sources[0].src.center_frequency
-                except:
-                    # No sources
-                    frequency = 0
 
     # validate the output plane to ensure proper 2D coordinates
     from meep.simulation import Volume
@@ -553,35 +556,48 @@ def plot2D(sim,ax=None, output_plane=None, fields=None, labels=False,
 
     # Plot geometry
     if plot_eps_flag:
-        ax = plot_eps(sim,ax,output_plane=output_plane,eps_parameters=eps_parameters,frequency=frequency)
+        ax = plot_eps(sim, ax, output_plane=output_plane,
+                      eps_parameters=eps_parameters, resolution=resolution)
 
     # Plot boundaries
     if plot_boundaries_flag:
-        ax = plot_boundaries(sim,ax,output_plane=output_plane,boundary_parameters=boundary_parameters)
+        ax = plot_boundaries(sim, ax, output_plane=output_plane,
+                             boundary_parameters=boundary_parameters)
 
     # Plot sources
     if plot_sources_flag:
-        ax = plot_sources(sim,ax,output_plane=output_plane,labels=labels,source_parameters=source_parameters)
+        ax = plot_sources(sim, ax, output_plane=output_plane,
+                          labels=labels, source_parameters=source_parameters)
 
     # Plot monitors
     if plot_monitors_flag:
-        ax = plot_monitors(sim,ax,output_plane=output_plane,labels=labels,monitor_parameters=monitor_parameters)
+        ax = plot_monitors(sim, ax, output_plane=output_plane,
+                           labels=labels, monitor_parameters=monitor_parameters)
 
     # Plot fields
-    ax = plot_fields(sim,ax,fields,output_plane=output_plane,field_parameters=field_parameters)
+    if fields:
+        ax = plot_fields(sim, ax, fields, output_plane=output_plane,
+                         field_parameters=field_parameters)
 
     return ax
 
 def plot3D(sim):
     from mayavi import mlab
 
-    if not sim._is_initialized:
-        sim.init_sim()
-
     if sim.dimensions < 3:
         raise ValueError("Simulation must have 3 dimensions to visualize 3D")
 
-    eps_data = sim.get_epsilon()
+    xtics = np.linspace(sim.geometry_center.x - 0.5*sim.cell_size.x,
+                        sim.geometry_center.x + 0.5*sim.cell_size.x,
+                        sim.cell_size.x*sim.resolution + 1)
+    ytics = np.linspace(sim.geometry_center.y - 0.5*sim.cell_size.y,
+                        sim.geometry_center.y + 0.5*sim.cell_size.y,
+                        sim.cell_size.y*sim.resolution + 1)
+    ztics = np.linspace(sim.geometry_center.z - 0.5*sim.cell_size.z,
+                        sim.geometry_center.z + 0.5*sim.cell_size.z,
+                        sim.cell_size.z*sim.resolution + 1)
+
+    eps_data = sim.get_epsilon_grid(xtics, ytics, ztics)
     s = mlab.contour3d(eps_data, colormap="YlGnBu")
     return s
 
@@ -711,7 +727,7 @@ class Animate2D(object):
     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
     [step-function modifier](#step-function-modifiers). For example, one can record the
-    E<sub>z</sub> fields at every one time unit using:
+    $E_z$ fields at every one time unit using:
 
     ```py
     animate = mp.Animate2D(sim,
@@ -792,7 +808,6 @@ def mod1(ax):
         self.plot_modifiers = plot_modifiers
         self.customization_args = customization_args
 
-
         self.cumulative_fields = []
         self._saved_frames = []
 
@@ -813,7 +828,7 @@ def __call__(self,sim,todo):
             # Initialize the plot
             if not self.init:
                 filtered_plot2D = filter_dict(self.customization_args, plot2D)
-                ax = sim.plot2D(ax=self.ax,fields=self.fields,**filtered_plot2D)
+                ax = sim.plot2D(ax=self.ax, fields=self.fields, **filtered_plot2D)
                 # Run the plot modifier functions
                 if self.plot_modifiers:
                     for k in range(len(self.plot_modifiers)):
@@ -827,7 +842,7 @@ def __call__(self,sim,todo):
             else:
                 # Update the plot
                 filtered_plot_fields= filter_dict(self.customization_args, plot_fields)
-                fields = sim.plot_fields(fields=self.fields,**filtered_plot_fields)
+                fields = sim.plot_fields(fields=self.fields, **filtered_plot_fields)
                 if mp.am_master():
                     self.ax.images[-1].set_data(fields)
                     self.ax.images[-1].set_clim(vmin=0.8*np.min(fields), vmax=0.8*np.max(fields))
@@ -851,7 +866,7 @@ def __call__(self,sim,todo):
             if self.normalize and mp.am_master():
                 if mp.verbosity.meep > 0:
                     print("Normalizing field data...")
-                fields = np.array(self.cumulative_fields) / np.max(np.abs(self.cumulative_fields),axis=(0,1,2))
+                fields = np.array(self.cumulative_fields) / np.max(np.abs(self.cumulative_fields), axis=(0,1,2))
                 for k in range(len(self.cumulative_fields)):
                     self.ax.images[-1].set_data(fields[k,:,:])
                     self.ax.images[-1].set_clim(vmin=-0.8, vmax=0.8)

From 771560e9efef217399856a7c0a2914b49582229e Mon Sep 17 00:00:00 2001
From: Ardavan Oskooi <ardavano@google.com>
Date: Sat, 13 Nov 2021 08:47:20 -0800
Subject: [PATCH 2/8] update docstrings

---
 doc/docs/Python_User_Interface.md    | 18 +++++++--------
 doc/docs/Python_User_Interface.md.in |  2 +-
 python/simulation.py                 | 12 ++++++----
 python/visualization.py              | 34 ++++++++++++++--------------
 4 files changed, 35 insertions(+), 31 deletions(-)

diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md
index 1904b9fc7..488cf8dcf 100644
--- a/doc/docs/Python_User_Interface.md
+++ b/doc/docs/Python_User_Interface.md
@@ -2067,7 +2067,7 @@ including outside the cell and a `near2far` object, returns the computed
 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,...)$
 in Cartesian coordinates and
 $(E_r^1,E_\phi^1,E_z^1,H_r^1,H_\phi^1,H_z^1,E_r^2,E_\phi^2,E_z^2,H_r^2,H_\phi^2,H_z^2,...)$
-in cylindrical coordinates for the frequencies 1,2,…,`nfreq`.
+in cylindrical coordinates for the frequencies 1,2,...,`nfreq`.
 
 </div>
 
@@ -2609,7 +2609,7 @@ 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:
+This module provides basic visualization functionality for the simulation domain. The intent 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:
 
 
 <a id="Simulation.plot2D"></a>
@@ -2648,8 +2648,8 @@ 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})
+           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')
 ```
@@ -2657,6 +2657,10 @@ plt.savefig('sim_domain.png')
 Note: When running a [parallel simulation](Parallel_Meep.md), the `plot2D` function expects to be called
 on all processes, but only generates a plot on the master process.
 
+Note: The geometry function is evaluated at infinite frequency (i.e., the instantaneous response).
+Dispersive materials are ignored. To visualize dispersive materials, you can use the `get_array` or
+`output_epsilon` functions for which a `frequency` parameter can be specified.
+
 **Parameters:**
 
 * `ax`: a `matplotlib` axis object. `plot2D()` will add plot objects, like lines,
@@ -2713,10 +2717,6 @@ on all processes, but only generates a plot on the master process.
     - `alpha=0.6`: transparency of fields
     - `post_process=np.real`: post processing function to apply to fields (must be
       a function object)
-* `frequency`: for materials with a [frequency-dependent
-  permittivity](Materials.md#material-dispersion) $\varepsilon(f)$, specifies the
-  frequency $f$ (in Meep units) of the real part of the permittivity to use in the
-  plot. Defaults to the `frequency` parameter of the [Source](#source) object.
 
 </div>
 
@@ -6893,7 +6893,7 @@ A class used to record the fields during timestepping (i.e., a [`run`](#run-func
 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
 [step-function modifier](#step-function-modifiers). For example, one can record the
-E<sub>z</sub> fields at every one time unit using:
+$E_z$ fields at every one time unit using:
 
 ```py
 animate = mp.Animate2D(sim,
diff --git a/doc/docs/Python_User_Interface.md.in b/doc/docs/Python_User_Interface.md.in
index 362f652ec..2ed09ca3d 100644
--- a/doc/docs/Python_User_Interface.md.in
+++ b/doc/docs/Python_User_Interface.md.in
@@ -480,7 +480,7 @@ This feature is only available if Meep is built with [libGDSII](Build_From_Sourc
 
 ### 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:
+This module provides basic visualization functionality for the simulation domain. The intent 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:
 
 @@ Simulation.plot2D @@
 @@ Simulation.plot3D @@
diff --git a/python/simulation.py b/python/simulation.py
index be93866cf..4bfa69d6b 100644
--- a/python/simulation.py
+++ b/python/simulation.py
@@ -4084,8 +4084,8 @@ def plot2D(self, ax=None, output_plane=None, fields=None, labels=False,
         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})
+                   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')
         ```
@@ -4093,6 +4093,10 @@ def plot2D(self, ax=None, output_plane=None, fields=None, labels=False,
         Note: When running a [parallel simulation](Parallel_Meep.md), the `plot2D` function expects to be called
         on all processes, but only generates a plot on the master process.
 
+        Note: The geometry function is evaluated at infinite frequency (i.e., the instantaneous response).
+        Dispersive materials are ignored. To visualize dispersive materials, you can use the `get_array` or
+        `output_epsilon` functions for which a `frequency` parameter can be specified.
+
         **Parameters:**
 
         * `ax`: a `matplotlib` axis object. `plot2D()` will add plot objects, like lines,
@@ -4166,8 +4170,8 @@ def plot2D(self, ax=None, output_plane=None, fields=None, labels=False,
                           plot_boundaries_flag=plot_boundaries_flag,
                           **kwargs)
 
-    def plot_fields(self,**kwargs):
-        return vis.plot_fields(self,**kwargs)
+    def plot_fields(self, **kwargs):
+        return vis.plot_fields(self, **kwargs)
 
     def plot3D(self):
         """
diff --git a/python/visualization.py b/python/visualization.py
index e2bf3de62..db2a1d968 100644
--- a/python/visualization.py
+++ b/python/visualization.py
@@ -348,15 +348,15 @@ def plot_eps(sim, ax, output_plane=None, eps_parameters=None, resolution=0):
 
     if sim_size.x == 0:
         # Plot y on x axis, z on y axis (YZ plane)
-        extent = [ymin,ymax,zmin,zmax]
+        extent = [ymin, ymax, zmin, zmax]
         xlabel = 'Y'
         ylabel = 'Z'
         xtics = np.array([0])
-        ytics = np.linspace(ymin,ymax, grid_resolution)
-        ztics = np.linspace(zmin,zmax, grid_resolution)
+        ytics = np.linspace(ymin, ymax, grid_resolution)
+        ztics = np.linspace(zmin, zmax, grid_resolution)
     elif sim_size.y == 0:
         # Plot x on x axis, z on y axis (XZ plane)
-        extent = [xmin,xmax,zmin,zmax]
+        extent = [xmin, xmax, zmin, zmax]
         xlabel = 'X'
         ylabel = 'Z'
         xtics = np.linspace(xmin, xmax, grid_resolution)
@@ -364,7 +364,7 @@ def plot_eps(sim, ax, output_plane=None, eps_parameters=None, resolution=0):
         ztics = np.linspace(zmin, zmax, grid_resolution)
     elif sim_size.z == 0:
         # Plot x on x axis, y on y axis (XY plane)
-        extent = [xmin,xmax,ymin,ymax]
+        extent = [xmin, xmax, ymin, ymax]
         xlabel = 'X'
         ylabel = 'Y'
         xtics = np.linspace(xmin, xmax, grid_resolution)
@@ -431,20 +431,20 @@ def get_boundary_volumes(thickness, direction, side):
 
     import itertools
     for boundary in sim.boundary_layers:
-        # All 4 side are the same
+        # All four sides are the same
         if boundary.direction == mp.ALL and boundary.side == mp.ALL:
             if sim.dimensions == 1:
                 dims = [mp.X]
             elif sim.dimensions == 2:
-                dims = [mp.X,mp.Y]
+                dims = [mp.X, mp.Y]
             elif sim.dimensions == 3:
-                dims = [mp.X,mp.Y,mp.Z]
+                dims = [mp.X, mp.Y, mp.Z]
             else:
                 raise ValueError("Invalid simulation dimensions")
             for permutation in itertools.product(dims, [mp.Low, mp.High]):
                 vol = get_boundary_volumes(boundary.thickness,*permutation)
                 ax = plot_volume(sim,ax,vol,output_plane,plotting_parameters=boundary_parameters)
-        # 2 sides are the same
+        # two sides are the same
         elif boundary.side == mp.ALL:
             for side in [mp.Low, mp.High]:
                 vol = get_boundary_volumes(boundary.thickness,boundary.direction,side)
@@ -494,7 +494,7 @@ def plot_fields(sim, ax=None, fields=None, output_plane=None, field_parameters=N
     # user specifies a field component
     if fields in [mp.Ex, mp.Ey, mp.Ez, mp.Hx, mp.Hy, mp.Hz]:
         # Get domain measurements
-        sim_center, sim_size = get_2D_dimensions(sim,output_plane)
+        sim_center, sim_size = get_2D_dimensions(sim, output_plane)
 
         xmin = sim_center.x - sim_size.x/2
         xmax = sim_center.x + sim_size.x/2
@@ -503,22 +503,22 @@ def plot_fields(sim, ax=None, fields=None, output_plane=None, field_parameters=N
         zmin = sim_center.z - sim_size.z/2
         zmax = sim_center.z + sim_size.z/2
 
-        center = Vector3(sim_center.x,sim_center.y,sim_center.z)
-        cell_size = Vector3(sim_size.x,sim_size.y,sim_size.z)
+        center = Vector3(sim_center.x, sim_center.y, sim_center.z)
+        cell_size = Vector3(sim_size.x, sim_size.y, sim_size.z)
 
         if sim_size.x == 0:
             # Plot y on x axis, z on y axis (YZ plane)
-            extent = [ymin,ymax,zmin,zmax]
+            extent = [ymin, ymax, zmin, zmax]
             xlabel = 'Y'
             ylabel = 'Z'
         elif sim_size.y == 0:
             # Plot x on x axis, z on y axis (XZ plane)
-            extent = [xmin,xmax,zmin,zmax]
+            extent = [xmin, xmax, zmin, zmax]
             xlabel = 'X'
             ylabel = 'Z'
         elif sim_size.z == 0:
             # Plot x on x axis, y on y axis (XY plane)
-            extent = [xmin,xmax,ymin,ymax]
+            extent = [xmin, xmax, ymin, ymax]
             xlabel = 'X'
             ylabel = 'Y'
         fields = sim.get_array(center=center, size=cell_size, component=fields)
@@ -551,8 +551,8 @@ def plot2D(sim, ax=None, output_plane=None, fields=None, labels=False,
 
     # validate the output plane to ensure proper 2D coordinates
     from meep.simulation import Volume
-    sim_center, sim_size = get_2D_dimensions(sim,output_plane)
-    output_plane = Volume(center=sim_center,size=sim_size)
+    sim_center, sim_size = get_2D_dimensions(sim, output_plane)
+    output_plane = Volume(center=sim_center, size=sim_size)
 
     # Plot geometry
     if plot_eps_flag:

From d748b34c10b62f097d5fee8dfbcc093b4e69909e Mon Sep 17 00:00:00 2001
From: Ardavan Oskooi <ardavano@google.com>
Date: Sat, 13 Nov 2021 14:09:10 -0800
Subject: [PATCH 3/8] rotate epsilon grid by 90 degrees to properly orient axes

---
 doc/docs/Python_User_Interface.md | 17 +++++++----
 python/simulation.py              | 18 +++++++----
 python/visualization.py           | 50 ++++++++++++++++++-------------
 3 files changed, 54 insertions(+), 31 deletions(-)

diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md
index 488cf8dcf..67e8aa69e 100644
--- a/doc/docs/Python_User_Interface.md
+++ b/doc/docs/Python_User_Interface.md
@@ -2627,7 +2627,7 @@ def plot2D(self,
            source_parameters=None,
            monitor_parameters=None,
            field_parameters=None,
-           frequency=None,
+           resolution=None,
            plot_eps_flag=True,
            plot_sources_flag=True,
            plot_monitors_flag=True,
@@ -2653,13 +2653,18 @@ sim.plot2D(fields=mp.Ez,
 plt.show()
 plt.savefig('sim_domain.png')
 ```
+If you just want to quickly visualize the simulation domain without the fields (i.e., when
+debugging your simulation), there is no need to invoke the `run` function prior to calling
+`plot2D`. Just define the `Simulation` object followed by any DFT monitors and then
+invoke `plot2D`.
 
-Note: When running a [parallel simulation](Parallel_Meep.md), the `plot2D` function expects to be called
-on all processes, but only generates a plot on the master process.
+Note: When running a [parallel simulation](Parallel_Meep.md), the `plot2D` function expects
+to be called on all processes, but only generates a plot on the master process.
 
 Note: The geometry function is evaluated at infinite frequency (i.e., the instantaneous response).
-Dispersive materials are ignored. To visualize dispersive materials, you can use the `get_array` or
-`output_epsilon` functions for which a `frequency` parameter can be specified.
+Dispersive materials containing complex, wavelength-dependent $\varepsilon$ are ignored. To
+visualize dispersive materials, you can use the `get_array` or `output_epsilon` functions and
+specify a `frequency` parameter for evaluating the complex $\varepsilon$.
 
 **Parameters:**
 
@@ -2681,6 +2686,8 @@ Dispersive materials are ignored. To visualize dispersive materials, you can use
     - `alpha=1.0`: transparency of geometry
     - `contour=False`: if `True`, plot a contour of the geometry rather than its image
     - `contour_linewidth=1`: line width of the contour lines if `contour=True`
+    - `resolution=None`: the resolution to sample the $\varepsilon$ grid. Defaults to the
+       `resolution` of the `Simulation` object.
 * `boundary_parameters`: a `dict` of optional plotting parameters that override
   the default parameters for the boundary layers.
     - `alpha=1.0`: transparency of boundary layers
diff --git a/python/simulation.py b/python/simulation.py
index 4bfa69d6b..bb967574f 100644
--- a/python/simulation.py
+++ b/python/simulation.py
@@ -4068,7 +4068,7 @@ def get_sfield_p(self,snap=False):
     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,
+               field_parameters=None, resolution=None,
                plot_eps_flag=True, plot_sources_flag=True,
                plot_monitors_flag=True, plot_boundaries_flag=True,
                **kwargs):
@@ -4089,13 +4089,18 @@ def plot2D(self, ax=None, output_plane=None, fields=None, labels=False,
         plt.show()
         plt.savefig('sim_domain.png')
         ```
+        If you just want to quickly visualize the simulation domain without the fields (i.e., when
+        debugging your simulation), there is no need to invoke the `run` function prior to calling
+        `plot2D`. Just define the `Simulation` object followed by any DFT monitors and then
+        invoke `plot2D`.
 
-        Note: When running a [parallel simulation](Parallel_Meep.md), the `plot2D` function expects to be called
-        on all processes, but only generates a plot on the master process.
+        Note: When running a [parallel simulation](Parallel_Meep.md), the `plot2D` function expects
+        to be called on all processes, but only generates a plot on the master process.
 
         Note: The geometry function is evaluated at infinite frequency (i.e., the instantaneous response).
-        Dispersive materials are ignored. To visualize dispersive materials, you can use the `get_array` or
-        `output_epsilon` functions for which a `frequency` parameter can be specified.
+        Dispersive materials containing complex, wavelength-dependent $\\varepsilon$ are ignored. To
+        visualize dispersive materials, you can use the `get_array` or `output_epsilon` functions and
+        specify a `frequency` parameter for evaluating the complex $\\varepsilon$.
 
         **Parameters:**
 
@@ -4117,6 +4122,8 @@ def plot2D(self, ax=None, output_plane=None, fields=None, labels=False,
             - `alpha=1.0`: transparency of geometry
             - `contour=False`: if `True`, plot a contour of the geometry rather than its image
             - `contour_linewidth=1`: line width of the contour lines if `contour=True`
+            - `resolution=None`: the resolution to sample the $\\varepsilon$ grid. Defaults to the
+               `resolution` of the `Simulation` object.
         * `boundary_parameters`: a `dict` of optional plotting parameters that override
           the default parameters for the boundary layers.
             - `alpha=1.0`: transparency of boundary layers
@@ -4164,6 +4171,7 @@ def plot2D(self, ax=None, output_plane=None, fields=None, labels=False,
                           source_parameters=source_parameters,
                           monitor_parameters=monitor_parameters,
                           field_parameters=field_parameters,
+                          resolution=resolution,
                           plot_eps_flag=plot_eps_flag,
                           plot_sources_flag=plot_sources_flag,
                           plot_monitors_flag=plot_monitors_flag,
diff --git a/python/visualization.py b/python/visualization.py
index db2a1d968..02d668ce0 100644
--- a/python/visualization.py
+++ b/python/visualization.py
@@ -327,7 +327,7 @@ def sort_points(xy):
             return ax
     return ax
 
-def plot_eps(sim, ax, output_plane=None, eps_parameters=None, resolution=0):
+def plot_eps(sim, ax, output_plane=None, eps_parameters=None, resolution=None):
     # consolidate plotting parameters
     eps_parameters = default_eps_parameters if eps_parameters is None else dict(default_eps_parameters, **eps_parameters)
 
@@ -344,36 +344,39 @@ def plot_eps(sim, ax, output_plane=None, eps_parameters=None, resolution=0):
     center = Vector3(sim_center.x, sim_center.y, sim_center.z)
     cell_size = Vector3(sim_size.x, sim_size.y, sim_size.z)
 
-    grid_resolution = resolution if resolution else sim.resolution + 1
+    grid_resolution = resolution if resolution else sim.resolution
+    Nx = int((xmax - xmin) * grid_resolution + 1)
+    Ny = int((ymax - ymin) * grid_resolution + 1)
+    Nz = int((zmax - zmin) * grid_resolution + 1)
 
     if sim_size.x == 0:
         # Plot y on x axis, z on y axis (YZ plane)
         extent = [ymin, ymax, zmin, zmax]
         xlabel = 'Y'
         ylabel = 'Z'
-        xtics = np.array([0])
-        ytics = np.linspace(ymin, ymax, grid_resolution)
-        ztics = np.linspace(zmin, zmax, grid_resolution)
+        xtics = np.array([sim_center.x])
+        ytics = np.linspace(ymin, ymax, Ny)
+        ztics = np.linspace(zmin, zmax, Nz)
     elif sim_size.y == 0:
         # Plot x on x axis, z on y axis (XZ plane)
         extent = [xmin, xmax, zmin, zmax]
         xlabel = 'X'
         ylabel = 'Z'
-        xtics = np.linspace(xmin, xmax, grid_resolution)
-        ytics = np.array([0])
-        ztics = np.linspace(zmin, zmax, grid_resolution)
+        xtics = np.linspace(xmin, xmax, Nx)
+        ytics = np.array([sim_center.y])
+        ztics = np.linspace(zmin, zmax, Nz)
     elif sim_size.z == 0:
         # Plot x on x axis, y on y axis (XY plane)
         extent = [xmin, xmax, ymin, ymax]
         xlabel = 'X'
         ylabel = 'Y'
-        xtics = np.linspace(xmin, xmax, grid_resolution)
-        ytics = np.linspace(ymin, ymax, grid_resolution)
-        ztics = np.array([0])
+        xtics = np.linspace(xmin, xmax, Nx)
+        ytics = np.linspace(ymin, ymax, Ny)
+        ztics = np.array([sim_center.z])
     else:
         raise ValueError("A 2D plane has not been specified...")
 
-    eps_data = sim.get_epsilon_grid(xtics, ytics, ztics)
+    eps_data = np.rot90(sim.get_epsilon_grid(xtics, ytics, ztics))
 
     if mp.am_master():
         if eps_parameters['contour']:
@@ -587,15 +590,20 @@ def plot3D(sim):
     if sim.dimensions < 3:
         raise ValueError("Simulation must have 3 dimensions to visualize 3D")
 
-    xtics = np.linspace(sim.geometry_center.x - 0.5*sim.cell_size.x,
-                        sim.geometry_center.x + 0.5*sim.cell_size.x,
-                        sim.cell_size.x*sim.resolution + 1)
-    ytics = np.linspace(sim.geometry_center.y - 0.5*sim.cell_size.y,
-                        sim.geometry_center.y + 0.5*sim.cell_size.y,
-                        sim.cell_size.y*sim.resolution + 1)
-    ztics = np.linspace(sim.geometry_center.z - 0.5*sim.cell_size.z,
-                        sim.geometry_center.z + 0.5*sim.cell_size.z,
-                        sim.cell_size.z*sim.resolution + 1)
+    xmin = sim.geometry_center.x - 0.5*sim.cell_size.x
+    xmax = sim.geometry_center.x + 0.5*sim.cell_size.x
+    ymin = sim.geometry_center.y - 0.5*sim.cell_size.y
+    ymax = sim.geometry_center.y + 0.5*sim.cell_size.y
+    zmin = sim.geometry_center.z - 0.5*sim.cell_size.z
+    zmax = sim.geometry_center.z + 0.5*sim.cell_size.z
+
+    Nx = int(sim.cell_size.x * sim.resolution) + 1
+    Ny = int(sim.cell_size.y * sim.resolution) + 1
+    Nz = int(sim.cell_size.z * sim.resolution) + 1
+
+    xtics = np.linspace(xmin, xmax, Nx)
+    ytics = np.linspace(ymin, ymax, Ny)
+    ztics = np.linspace(zmin, zmax, Nz)
 
     eps_data = sim.get_epsilon_grid(xtics, ytics, ztics)
     s = mlab.contour3d(eps_data, colormap="YlGnBu")

From a7fb4c3bde5a36fb81eeb189f0ad5f73c9b066ee Mon Sep 17 00:00:00 2001
From: Ardavan Oskooi <ardavano@google.com>
Date: Wed, 17 Nov 2021 17:47:05 -0800
Subject: [PATCH 4/8] =?UTF-8?q?add=20support=20for=20dispersive=20=CE=B5?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 python/meep.i           | 12 +++++++-----
 python/simulation.py    | 23 +++++++++++++++--------
 python/visualization.py | 26 ++++++++++++++++++++++----
 src/meepgeom.cpp        | 27 +++++++++++++++++++--------
 src/meepgeom.hpp        |  3 ++-
 5 files changed, 65 insertions(+), 26 deletions(-)

diff --git a/python/meep.i b/python/meep.i
index fe301cf58..1f9c5eb03 100644
--- a/python/meep.i
+++ b/python/meep.i
@@ -1097,12 +1097,12 @@ void _get_gradient(PyObject *grad, double scalegrad, PyObject *fields_a, PyObjec
     $1 = (double *)array_data($input);
 }
 
-%typecheck(SWIG_TYPECHECK_POINTER, fragment="NumPy_Fragments") double* grid_vals {
+%typecheck(SWIG_TYPECHECK_POINTER, fragment="NumPy_Fragments") std::complex<double>* grid_vals {
     $1 = is_array($input);
 }
 
-%typemap(in, fragment="NumPy_Macros") double* grid_vals {
-    $1 = (double *)array_data($input);
+%typemap(in, fragment="NumPy_Macros") std::complex<double>* grid_vals {
+    $1 = (std::complex<double> *)array_data($input);
 }
 
 // typemap for solve_cw:
@@ -2041,7 +2041,8 @@ void _get_epsilon_grid(geometric_object_list gobj_list,
                        int nx, double *xtics,
                        int ny, double *ytics,
                        int nz, double *ztics,
-                       double *grid_vals) {
+                       std::complex<double> *grid_vals,
+                       double frequency) {
      meep_geom::get_epsilon_grid(gobj_list,
                                  mlist,
                                  _default_material,
@@ -2052,7 +2053,8 @@ void _get_epsilon_grid(geometric_object_list gobj_list,
                                  nx, xtics,
                                  ny, ytics,
                                  nz, ztics,
-                                 grid_vals);
+                                 grid_vals,
+                                 frequency);
 }
 
 %}
diff --git a/python/simulation.py b/python/simulation.py
index bb967574f..2feab2033 100644
--- a/python/simulation.py
+++ b/python/simulation.py
@@ -2225,18 +2225,19 @@ def get_mu_point(self, pt, frequency=0):
         v3 = py_v3_to_vec(self.dimensions, pt, self.is_cylindrical)
         return self.fields.get_mu(v3,frequency)
 
-    def get_epsilon_grid(self, xtics=None, ytics=None, ztics=None):
+    def get_epsilon_grid(self, xtics=None, ytics=None, ztics=None, frequency=0):
         """
         Given three 1d NumPy arrays (`xtics`,`ytics`,`ztics`) which define the coordinates of a Cartesian
-        grid anywhere within the cell volume, compute the trace of the $\\varepsilon$ tensor from the `geometry`
-        exactly at each grid point. (For [`MaterialGrid`](#materialgrid)s, the $\\varepsilon$ at each grid
-        point is computed using bilinear interpolation from the nearest `MaterialGrid` points and possibly also
-        projected to form a level set.) Note that this is different from `get_epsilon_point` which computes
+        grid anywhere within the cell volume, compute the trace of the $\\varepsilon(f)$ tensor at frequency
+        $f$ (in Meep units) from the `geometry` exactly at each grid point. `frequency` defaults to 0 which is
+        the instantaneous $\\varepsilon$. (For [`MaterialGrid`](#materialgrid)s, the $\\varepsilon$ at each
+        grid point is computed using bilinear interpolation from the nearest `MaterialGrid` points and possibly
+        also projected to form a level set.) Note that this is different from `get_epsilon_point` which computes
         $\\varepsilon$ by bilinearly interpolating from the nearest Yee grid points. This function is useful for
         sampling the material geometry to any arbitrary resolution. The return value is a NumPy array with shape
         equivalent to `numpy.meshgrid(xtics,ytics,ztics)`. Empty dimensions are collapsed.
         """
-        grid_vals = np.squeeze(np.empty((len(xtics), len(ytics), len(ztics))))
+        grid_vals = np.squeeze(np.empty((len(xtics), len(ytics), len(ztics)), dtype=np.complex128))
         gv = self._create_grid_volume(False)
         mp._get_epsilon_grid(self.geometry,
                              self.extra_materials,
@@ -2247,7 +2248,8 @@ def get_epsilon_grid(self, xtics=None, ytics=None, ztics=None):
                              len(xtics), xtics,
                              len(ytics), ytics,
                              len(ztics), ztics,
-                             grid_vals)
+                             grid_vals,
+                             frequency)
         return grid_vals
 
     def get_filename_prefix(self):
@@ -4068,7 +4070,7 @@ def get_sfield_p(self,snap=False):
     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, resolution=None,
+               field_parameters=None, frequency=None, resolution=None,
                plot_eps_flag=True, plot_sources_flag=True,
                plot_monitors_flag=True, plot_boundaries_flag=True,
                **kwargs):
@@ -4160,6 +4162,10 @@ def plot2D(self, ax=None, output_plane=None, fields=None, labels=False,
             - `alpha=0.6`: transparency of fields
             - `post_process=np.real`: post processing function to apply to fields (must be
               a function object)
+        * `frequency`: for materials with a [frequency-dependent
+          permittivity](Materials.md#material-dispersion) $\\varepsilon(f)$, specifies the
+          frequency $f$ (in Meep units) of the real part of the permittivity to use in the
+          plot. Defaults to the `frequency` parameter of the [Source](#source) object.
         """
         return vis.plot2D(self,
                           ax=ax,
@@ -4171,6 +4177,7 @@ def plot2D(self, ax=None, output_plane=None, fields=None, labels=False,
                           source_parameters=source_parameters,
                           monitor_parameters=monitor_parameters,
                           field_parameters=field_parameters,
+                          frequency=frequency,
                           resolution=resolution,
                           plot_eps_flag=plot_eps_flag,
                           plot_sources_flag=plot_sources_flag,
diff --git a/python/visualization.py b/python/visualization.py
index 02d668ce0..68beb1ed8 100644
--- a/python/visualization.py
+++ b/python/visualization.py
@@ -327,7 +327,7 @@ def sort_points(xy):
             return ax
     return ax
 
-def plot_eps(sim, ax, output_plane=None, eps_parameters=None, resolution=None):
+def plot_eps(sim, ax, output_plane=None, eps_parameters=None, resolution=None, frequency=0):
     # consolidate plotting parameters
     eps_parameters = default_eps_parameters if eps_parameters is None else dict(default_eps_parameters, **eps_parameters)
 
@@ -376,7 +376,7 @@ def plot_eps(sim, ax, output_plane=None, eps_parameters=None, resolution=None):
     else:
         raise ValueError("A 2D plane has not been specified...")
 
-    eps_data = np.rot90(sim.get_epsilon_grid(xtics, ytics, ztics))
+    eps_data = np.rot90(np.real(sim.get_epsilon_grid(xtics, ytics, ztics, frequency)))
 
     if mp.am_master():
         if eps_parameters['contour']:
@@ -543,7 +543,7 @@ def plot_fields(sim, ax=None, fields=None, output_plane=None, field_parameters=N
 def plot2D(sim, ax=None, output_plane=None, fields=None, labels=False,
            eps_parameters=None, boundary_parameters=None,
            source_parameters=None, monitor_parameters=None,
-           field_parameters=None, resolution=None,
+           field_parameters=None, frequency=None, resolution=None,
            plot_eps_flag=True, plot_sources_flag=True,
            plot_monitors_flag=True, plot_boundaries_flag=True):
 
@@ -552,6 +552,23 @@ def plot2D(sim, ax=None, output_plane=None, fields=None, labels=False,
         from matplotlib import pyplot as plt
         ax = plt.gca()
 
+    # Determine a frequency to plot all epsilon
+    if frequency is None:
+        try:
+            # Continuous sources
+            frequency = sim.sources[0].frequency
+        except:
+            try:
+                # Gaussian sources
+                frequency = sim.sources[0].src.frequency
+            except:
+                try:
+                    # Custom sources
+                    frequency = sim.sources[0].src.center_frequency
+                except:
+                    # No sources
+                    frequency = 0
+
     # validate the output plane to ensure proper 2D coordinates
     from meep.simulation import Volume
     sim_center, sim_size = get_2D_dimensions(sim, output_plane)
@@ -560,7 +577,8 @@ def plot2D(sim, ax=None, output_plane=None, fields=None, labels=False,
     # Plot geometry
     if plot_eps_flag:
         ax = plot_eps(sim, ax, output_plane=output_plane,
-                      eps_parameters=eps_parameters, resolution=resolution)
+                      eps_parameters=eps_parameters, resolution=resolution,
+                      frequency=frequency)
 
     # Plot boundaries
     if plot_boundaries_flag:
diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp
index 077f2f52b..46928be78 100644
--- a/src/meepgeom.cpp
+++ b/src/meepgeom.cpp
@@ -2557,13 +2557,11 @@ void invert_tensor(std::complex<double> t_inv[9], std::complex<double> t[9]) {
 #undef minv(x,y)
 }
 
-void eff_chi1inv_row_disp(meep::component c, std::complex<double> chi1inv_row[3],
-  const meep::vec &r, double freq, geom_epsilon *geps) {
+void get_chi1_tensor_disp(std::complex<double> tensor[9], const meep::vec &r, double freq, geom_epsilon *geps) {
   // locate the proper material
   material_type md;
   geps->get_material_pt(md, r);
   const medium_struct *mm = &(md->medium);
-  std::complex<double> tensor[9], tensor_inv[9];
 
   // loop over all the tensor components
   for (int i = 0; i < 9; i++) {
@@ -2574,7 +2572,7 @@ void eff_chi1inv_row_disp(meep::component c, std::complex<double> chi1inv_row[3]
     double conductivityCur = vec_to_value(mm->D_conductivity_diag, dummy, i);
     a = std::complex<double>(1.0, conductivityCur / (freq));
 
-    // compute lorentzian component
+    // compute lorentzian component including the instantaneous ε
     b = cvec_to_value(mm->epsilon_diag, mm->epsilon_offdiag, i);
     for (const auto &mm_susc: mm->E_susceptibilities) {
       meep::lorentzian_susceptibility sus = meep::lorentzian_susceptibility(
@@ -2586,7 +2584,12 @@ void eff_chi1inv_row_disp(meep::component c, std::complex<double> chi1inv_row[3]
     // elementwise multiply
     tensor[i] = a * b;
   }
+}
 
+void eff_chi1inv_row_disp(meep::component c, std::complex<double> chi1inv_row[3],
+  const meep::vec &r, double freq, geom_epsilon *geps) {
+  std::complex<double> tensor[9], tensor_inv[9];
+  get_chi1_tensor_disp(tensor, r, freq, geps);
   // invert the matrix
   invert_tensor(tensor_inv, tensor);
 
@@ -3028,7 +3031,8 @@ void get_epsilon_grid(geometric_object_list gobj_list,
                       int nx, const double *x,
                       int ny, const double *y,
                       int nz, const double *z,
-                      double *grid_vals) {
+                      std::complex<double> *grid_vals,
+                      double frequency) {
   double min_val[3], max_val[3];
   for (int n = 0; n < 3; ++n) {
     int ndir = (n == 0) ? nx : ((n == 1) ? ny : nz);
@@ -3044,10 +3048,17 @@ void get_epsilon_grid(geometric_object_list gobj_list,
   geom_epsilon geps(gobj_list, mlist, vol);
   for (int i = 0; i < nx; ++i)
     for (int j = 0; j < ny; ++j)
-      for (int k = 0; k < nz; ++k)
-        /* obtain the trace of the \varepsilon tensor for each
+      for (int k = 0; k < nz; ++k) {
+        /* obtain the trace of the ε tensor (dispersive or non) for each
            grid point in row-major order (the order used by NumPy) */
-        grid_vals[k + nz*(j + ny*i)] = geps.chi1p1(meep::E_stuff, meep::vec(x[i],y[j],z[k]));
+        if (frequency == 0)
+          grid_vals[k + nz*(j + ny*i)] = geps.chi1p1(meep::E_stuff, meep::vec(x[i],y[j],z[k]));
+        else {
+          std::complex<double> tensor[9];
+          get_chi1_tensor_disp(tensor, meep::vec(x[i],y[j],z[k]), frequency, &geps);
+          grid_vals[k + nz*(j + ny*i)] = (tensor[0] + tensor[4] + tensor[8]) / 3.0;
+        }
+      }
 }
 
 } // namespace meep_geom
diff --git a/src/meepgeom.hpp b/src/meepgeom.hpp
index d052d5b4a..f66c33e3b 100644
--- a/src/meepgeom.hpp
+++ b/src/meepgeom.hpp
@@ -281,7 +281,8 @@ void get_epsilon_grid(geometric_object_list gobj_list,
                       int nx, const double *x,
                       int ny, const double *y,
                       int nz, const double *z,
-                      double *grid_vals);
+                      std::complex<double> *grid_vals,
+                      double frequency = 0);
 void init_libctl(material_type default_mat, bool ensure_per,
                  meep::grid_volume *gv, vector3 cell_size, vector3 cell_center,
                  geometric_object_list *geom_list);

From b47d272ca94fc515834c695b803386d69b7a2148 Mon Sep 17 00:00:00 2001
From: Ardavan Oskooi <ardavano@google.com>
Date: Wed, 17 Nov 2021 17:57:43 -0800
Subject: [PATCH 5/8] update markdown pages from docstrings

---
 doc/docs/Python_User_Interface.md | 22 ++++++++++++++++------
 1 file changed, 16 insertions(+), 6 deletions(-)

diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md
index 67e8aa69e..20293337c 100644
--- a/doc/docs/Python_User_Interface.md
+++ b/doc/docs/Python_User_Interface.md
@@ -676,16 +676,21 @@ frequency-independent part of $\mu$ (the $\omega\to\infty$ limit).
 <div class="class_members" markdown="1">
 
 ```python
-def get_epsilon_grid(self, xtics=None, ytics=None, ztics=None):
+def get_epsilon_grid(self,
+                     xtics=None,
+                     ytics=None,
+                     ztics=None,
+                     frequency=0):
 ```
 
 <div class="method_docstring" markdown="1">
 
 Given three 1d NumPy arrays (`xtics`,`ytics`,`ztics`) which define the coordinates of a Cartesian
-grid anywhere within the cell volume, compute the trace of the $\varepsilon$ tensor from the `geometry`
-exactly at each grid point. (For [`MaterialGrid`](#materialgrid)s, the $\varepsilon$ at each grid
-point is computed using bilinear interpolation from the nearest `MaterialGrid` points and possibly also
-projected to form a level set.) Note that this is different from `get_epsilon_point` which computes
+grid anywhere within the cell volume, compute the trace of the $\varepsilon(f)$ tensor at frequency
+$f$ (in Meep units) from the `geometry` exactly at each grid point. `frequency` defaults to 0 which is
+the instantaneous $\varepsilon$. (For [`MaterialGrid`](#materialgrid)s, the $\varepsilon$ at each
+grid point is computed using bilinear interpolation from the nearest `MaterialGrid` points and possibly
+also projected to form a level set.) Note that this is different from `get_epsilon_point` which computes
 $\varepsilon$ by bilinearly interpolating from the nearest Yee grid points. This function is useful for
 sampling the material geometry to any arbitrary resolution. The return value is a NumPy array with shape
 equivalent to `numpy.meshgrid(xtics,ytics,ztics)`. Empty dimensions are collapsed.
@@ -2627,6 +2632,7 @@ def plot2D(self,
            source_parameters=None,
            monitor_parameters=None,
            field_parameters=None,
+           frequency=None,
            resolution=None,
            plot_eps_flag=True,
            plot_sources_flag=True,
@@ -2724,6 +2730,10 @@ specify a `frequency` parameter for evaluating the complex $\varepsilon$.
     - `alpha=0.6`: transparency of fields
     - `post_process=np.real`: post processing function to apply to fields (must be
       a function object)
+* `frequency`: for materials with a [frequency-dependent
+  permittivity](Materials.md#material-dispersion) $\varepsilon(f)$, specifies the
+  frequency $f$ (in Meep units) of the real part of the permittivity to use in the
+  plot. Defaults to the `frequency` parameter of the [Source](#source) object.
 
 </div>
 
@@ -4408,7 +4418,7 @@ def __init__(self,
              medium2,
              weights=None,
              grid_type='U_DEFAULT',
-             do_averaging=False,
+             do_averaging=True,
              beta=0,
              eta=0.5,
              damping=0):

From 00f7b480ee920523ef3eadd1d349d9600f061626 Mon Sep 17 00:00:00 2001
From: Ardavan Oskooi <ardavano@google.com>
Date: Tue, 23 Nov 2021 17:45:24 -0800
Subject: [PATCH 6/8] make frequency and resolution parameters of plot2D
 dictionary keys of eps_parameters

---
 doc/docs/Python_User_Interface.md | 21 ++++--------
 python/simulation.py              | 28 ++++++----------
 python/visualization.py           | 53 ++++++++++++++++---------------
 3 files changed, 44 insertions(+), 58 deletions(-)

diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md
index 20293337c..87064f006 100644
--- a/doc/docs/Python_User_Interface.md
+++ b/doc/docs/Python_User_Interface.md
@@ -2632,8 +2632,6 @@ def plot2D(self,
            source_parameters=None,
            monitor_parameters=None,
            field_parameters=None,
-           frequency=None,
-           resolution=None,
            plot_eps_flag=True,
            plot_sources_flag=True,
            plot_monitors_flag=True,
@@ -2660,18 +2658,13 @@ plt.show()
 plt.savefig('sim_domain.png')
 ```
 If you just want to quickly visualize the simulation domain without the fields (i.e., when
-debugging your simulation), there is no need to invoke the `run` function prior to calling
+setting up your simulation), there is no need to invoke the `run` function prior to calling
 `plot2D`. Just define the `Simulation` object followed by any DFT monitors and then
 invoke `plot2D`.
 
 Note: When running a [parallel simulation](Parallel_Meep.md), the `plot2D` function expects
 to be called on all processes, but only generates a plot on the master process.
 
-Note: The geometry function is evaluated at infinite frequency (i.e., the instantaneous response).
-Dispersive materials containing complex, wavelength-dependent $\varepsilon$ are ignored. To
-visualize dispersive materials, you can use the `get_array` or `output_epsilon` functions and
-specify a `frequency` parameter for evaluating the complex $\varepsilon$.
-
 **Parameters:**
 
 * `ax`: a `matplotlib` axis object. `plot2D()` will add plot objects, like lines,
@@ -2692,8 +2685,12 @@ specify a `frequency` parameter for evaluating the complex $\varepsilon$.
     - `alpha=1.0`: transparency of geometry
     - `contour=False`: if `True`, plot a contour of the geometry rather than its image
     - `contour_linewidth=1`: line width of the contour lines if `contour=True`
-    - `resolution=None`: the resolution to sample the $\varepsilon$ grid. Defaults to the
-       `resolution` of the `Simulation` object.
+    - `frequency=None`: for materials with a [frequency-dependent
+      permittivity](Materials.md#material-dispersion) $\varepsilon(f)$, specifies the
+      frequency $f$ (in Meep units) of the real part of the permittivity to use in the
+      plot. Defaults to the `frequency` parameter of the [Source](#source) object.
+    - `resolution=None`: the resolution of the $\varepsilon$ grid. Defaults to the
+      `resolution` of the `Simulation` object.
 * `boundary_parameters`: a `dict` of optional plotting parameters that override
   the default parameters for the boundary layers.
     - `alpha=1.0`: transparency of boundary layers
@@ -2730,10 +2727,6 @@ specify a `frequency` parameter for evaluating the complex $\varepsilon$.
     - `alpha=0.6`: transparency of fields
     - `post_process=np.real`: post processing function to apply to fields (must be
       a function object)
-* `frequency`: for materials with a [frequency-dependent
-  permittivity](Materials.md#material-dispersion) $\varepsilon(f)$, specifies the
-  frequency $f$ (in Meep units) of the real part of the permittivity to use in the
-  plot. Defaults to the `frequency` parameter of the [Source](#source) object.
 
 </div>
 
diff --git a/python/simulation.py b/python/simulation.py
index 2feab2033..09c1eed61 100644
--- a/python/simulation.py
+++ b/python/simulation.py
@@ -4070,10 +4070,9 @@ def get_sfield_p(self,snap=False):
     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, resolution=None,
-               plot_eps_flag=True, plot_sources_flag=True,
-               plot_monitors_flag=True, plot_boundaries_flag=True,
-               **kwargs):
+               field_parameters=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
@@ -4092,18 +4091,13 @@ def plot2D(self, ax=None, output_plane=None, fields=None, labels=False,
         plt.savefig('sim_domain.png')
         ```
         If you just want to quickly visualize the simulation domain without the fields (i.e., when
-        debugging your simulation), there is no need to invoke the `run` function prior to calling
+        setting up your simulation), there is no need to invoke the `run` function prior to calling
         `plot2D`. Just define the `Simulation` object followed by any DFT monitors and then
         invoke `plot2D`.
 
         Note: When running a [parallel simulation](Parallel_Meep.md), the `plot2D` function expects
         to be called on all processes, but only generates a plot on the master process.
 
-        Note: The geometry function is evaluated at infinite frequency (i.e., the instantaneous response).
-        Dispersive materials containing complex, wavelength-dependent $\\varepsilon$ are ignored. To
-        visualize dispersive materials, you can use the `get_array` or `output_epsilon` functions and
-        specify a `frequency` parameter for evaluating the complex $\\varepsilon$.
-
         **Parameters:**
 
         * `ax`: a `matplotlib` axis object. `plot2D()` will add plot objects, like lines,
@@ -4124,8 +4118,12 @@ def plot2D(self, ax=None, output_plane=None, fields=None, labels=False,
             - `alpha=1.0`: transparency of geometry
             - `contour=False`: if `True`, plot a contour of the geometry rather than its image
             - `contour_linewidth=1`: line width of the contour lines if `contour=True`
-            - `resolution=None`: the resolution to sample the $\\varepsilon$ grid. Defaults to the
-               `resolution` of the `Simulation` object.
+            - `frequency=None`: for materials with a [frequency-dependent
+              permittivity](Materials.md#material-dispersion) $\\varepsilon(f)$, specifies the
+              frequency $f$ (in Meep units) of the real part of the permittivity to use in the
+              plot. Defaults to the `frequency` parameter of the [Source](#source) object.
+            - `resolution=None`: the resolution of the $\\varepsilon$ grid. Defaults to the
+              `resolution` of the `Simulation` object.
         * `boundary_parameters`: a `dict` of optional plotting parameters that override
           the default parameters for the boundary layers.
             - `alpha=1.0`: transparency of boundary layers
@@ -4162,10 +4160,6 @@ def plot2D(self, ax=None, output_plane=None, fields=None, labels=False,
             - `alpha=0.6`: transparency of fields
             - `post_process=np.real`: post processing function to apply to fields (must be
               a function object)
-        * `frequency`: for materials with a [frequency-dependent
-          permittivity](Materials.md#material-dispersion) $\\varepsilon(f)$, specifies the
-          frequency $f$ (in Meep units) of the real part of the permittivity to use in the
-          plot. Defaults to the `frequency` parameter of the [Source](#source) object.
         """
         return vis.plot2D(self,
                           ax=ax,
@@ -4177,8 +4171,6 @@ def plot2D(self, ax=None, output_plane=None, fields=None, labels=False,
                           source_parameters=source_parameters,
                           monitor_parameters=monitor_parameters,
                           field_parameters=field_parameters,
-                          frequency=frequency,
-                          resolution=resolution,
                           plot_eps_flag=plot_eps_flag,
                           plot_sources_flag=plot_sources_flag,
                           plot_monitors_flag=plot_monitors_flag,
diff --git a/python/visualization.py b/python/visualization.py
index 68beb1ed8..b66f3729b 100644
--- a/python/visualization.py
+++ b/python/visualization.py
@@ -47,7 +47,9 @@
         'cmap':'binary',
         'alpha':1.0,
         'contour':False,
-        'contour_linewidth':1
+        'contour_linewidth':1,
+        'frequency':None,
+        'resolution':None
     }
 
 default_boundary_parameters = {
@@ -327,10 +329,27 @@ def sort_points(xy):
             return ax
     return ax
 
-def plot_eps(sim, ax, output_plane=None, eps_parameters=None, resolution=None, frequency=0):
+def plot_eps(sim, ax, output_plane=None, eps_parameters=None):
     # consolidate plotting parameters
     eps_parameters = default_eps_parameters if eps_parameters is None else dict(default_eps_parameters, **eps_parameters)
 
+    # Determine a frequency to plot all epsilon
+    if eps_parameters['frequency'] is None:
+        try:
+            # Continuous sources
+            eps_parameters['frequency'] = sim.sources[0].frequency
+        except:
+            try:
+                # Gaussian sources
+                eps_parameters['frequency'] = sim.sources[0].src.frequency
+            except:
+                try:
+                    # Custom sources
+                    eps_parameters['frequency'] = sim.sources[0].src.center_frequency
+                except:
+                    # No sources
+                    eps_parameters['frequency'] = 0
+
     # Get domain measurements
     sim_center, sim_size = get_2D_dimensions(sim, output_plane)
 
@@ -344,7 +363,7 @@ def plot_eps(sim, ax, output_plane=None, eps_parameters=None, resolution=None, f
     center = Vector3(sim_center.x, sim_center.y, sim_center.z)
     cell_size = Vector3(sim_size.x, sim_size.y, sim_size.z)
 
-    grid_resolution = resolution if resolution else sim.resolution
+    grid_resolution = eps_parameters['resolution'] if eps_parameters['resolution'] else sim.resolution
     Nx = int((xmax - xmin) * grid_resolution + 1)
     Ny = int((ymax - ymin) * grid_resolution + 1)
     Nz = int((zmax - zmin) * grid_resolution + 1)
@@ -376,7 +395,7 @@ def plot_eps(sim, ax, output_plane=None, eps_parameters=None, resolution=None, f
     else:
         raise ValueError("A 2D plane has not been specified...")
 
-    eps_data = np.rot90(np.real(sim.get_epsilon_grid(xtics, ytics, ztics, frequency)))
+    eps_data = np.rot90(np.real(sim.get_epsilon_grid(xtics, ytics, ztics, eps_parameters['frequency'])))
 
     if mp.am_master():
         if eps_parameters['contour']:
@@ -543,32 +562,15 @@ def plot_fields(sim, ax=None, fields=None, output_plane=None, field_parameters=N
 def plot2D(sim, 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, resolution=None,
-           plot_eps_flag=True, plot_sources_flag=True,
-           plot_monitors_flag=True, plot_boundaries_flag=True):
+           field_parameters=None, plot_eps_flag=True,
+           plot_sources_flag=True, plot_monitors_flag=True,
+           plot_boundaries_flag=True):
 
     # Ensure a figure axis exists
     if ax is None and mp.am_master():
         from matplotlib import pyplot as plt
         ax = plt.gca()
 
-    # Determine a frequency to plot all epsilon
-    if frequency is None:
-        try:
-            # Continuous sources
-            frequency = sim.sources[0].frequency
-        except:
-            try:
-                # Gaussian sources
-                frequency = sim.sources[0].src.frequency
-            except:
-                try:
-                    # Custom sources
-                    frequency = sim.sources[0].src.center_frequency
-                except:
-                    # No sources
-                    frequency = 0
-
     # validate the output plane to ensure proper 2D coordinates
     from meep.simulation import Volume
     sim_center, sim_size = get_2D_dimensions(sim, output_plane)
@@ -577,8 +579,7 @@ def plot2D(sim, ax=None, output_plane=None, fields=None, labels=False,
     # Plot geometry
     if plot_eps_flag:
         ax = plot_eps(sim, ax, output_plane=output_plane,
-                      eps_parameters=eps_parameters, resolution=resolution,
-                      frequency=frequency)
+                      eps_parameters=eps_parameters)
 
     # Plot boundaries
     if plot_boundaries_flag:

From 00c0c3bb5a6e9e0b7889bfac87735e59a5face49 Mon Sep 17 00:00:00 2001
From: Ardavan Oskooi <ardavano@google.com>
Date: Tue, 23 Nov 2021 18:48:58 -0800
Subject: [PATCH 7/8] reinstate frequency parameter of plot2D and add warning
 that it is deprecated

---
 doc/docs/Python_User_Interface.md |  1 +
 python/simulation.py              |  8 +++++---
 python/visualization.py           | 34 +++++++++++++++++--------------
 3 files changed, 25 insertions(+), 18 deletions(-)

diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md
index 87064f006..581dd3382 100644
--- a/doc/docs/Python_User_Interface.md
+++ b/doc/docs/Python_User_Interface.md
@@ -2632,6 +2632,7 @@ def plot2D(self,
            source_parameters=None,
            monitor_parameters=None,
            field_parameters=None,
+           frequency=None,
            plot_eps_flag=True,
            plot_sources_flag=True,
            plot_monitors_flag=True,
diff --git a/python/simulation.py b/python/simulation.py
index 09c1eed61..5fbf66155 100644
--- a/python/simulation.py
+++ b/python/simulation.py
@@ -4070,9 +4070,10 @@ def get_sfield_p(self,snap=False):
     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, plot_eps_flag=True,
-               plot_sources_flag=True, plot_monitors_flag=True,
-               plot_boundaries_flag=True, **kwargs):
+               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
@@ -4171,6 +4172,7 @@ def plot2D(self, ax=None, output_plane=None, fields=None, labels=False,
                           source_parameters=source_parameters,
                           monitor_parameters=monitor_parameters,
                           field_parameters=field_parameters,
+                          frequency=frequency,
                           plot_eps_flag=plot_eps_flag,
                           plot_sources_flag=plot_sources_flag,
                           plot_monitors_flag=plot_monitors_flag,
diff --git a/python/visualization.py b/python/visualization.py
index b66f3729b..29e4167b9 100644
--- a/python/visualization.py
+++ b/python/visualization.py
@@ -329,26 +329,30 @@ def sort_points(xy):
             return ax
     return ax
 
-def plot_eps(sim, ax, output_plane=None, eps_parameters=None):
+def plot_eps(sim, ax, output_plane=None, eps_parameters=None, frequency=None):
     # consolidate plotting parameters
     eps_parameters = default_eps_parameters if eps_parameters is None else dict(default_eps_parameters, **eps_parameters)
 
     # Determine a frequency to plot all epsilon
     if eps_parameters['frequency'] is None:
-        try:
-            # Continuous sources
-            eps_parameters['frequency'] = sim.sources[0].frequency
-        except:
+        if frequency is None:
             try:
-                # Gaussian sources
-                eps_parameters['frequency'] = sim.sources[0].src.frequency
+                # Continuous sources
+                eps_parameters['frequency'] = sim.sources[0].frequency
             except:
                 try:
-                    # Custom sources
-                    eps_parameters['frequency'] = sim.sources[0].src.center_frequency
+                    # Gaussian sources
+                    eps_parameters['frequency'] = sim.sources[0].src.frequency
                 except:
-                    # No sources
-                    eps_parameters['frequency'] = 0
+                    try:
+                        # Custom sources
+                        eps_parameters['frequency'] = sim.sources[0].src.center_frequency
+                    except:
+                        # No sources
+                        eps_parameters['frequency'] = 0
+        else:
+            warnings.warn('The frequency parameter of plot2D has been deprecated. Use the frequency key of the eps_parameters dictionary instead.')
+            eps_parameters['frequency'] = frequency
 
     # Get domain measurements
     sim_center, sim_size = get_2D_dimensions(sim, output_plane)
@@ -562,9 +566,9 @@ def plot_fields(sim, ax=None, fields=None, output_plane=None, field_parameters=N
 def plot2D(sim, ax=None, output_plane=None, fields=None, labels=False,
            eps_parameters=None, boundary_parameters=None,
            source_parameters=None, monitor_parameters=None,
-           field_parameters=None, plot_eps_flag=True,
-           plot_sources_flag=True, plot_monitors_flag=True,
-           plot_boundaries_flag=True):
+           field_parameters=None, frequency=None,
+           plot_eps_flag=True, plot_sources_flag=True,
+           plot_monitors_flag=True, plot_boundaries_flag=True):
 
     # Ensure a figure axis exists
     if ax is None and mp.am_master():
@@ -579,7 +583,7 @@ def plot2D(sim, ax=None, output_plane=None, fields=None, labels=False,
     # Plot geometry
     if plot_eps_flag:
         ax = plot_eps(sim, ax, output_plane=output_plane,
-                      eps_parameters=eps_parameters)
+                      eps_parameters=eps_parameters, frequency=frequency)
 
     # Plot boundaries
     if plot_boundaries_flag:

From 6d1fa5cd7b50a3d55cc0f1cc8df7898e9893c38e Mon Sep 17 00:00:00 2001
From: Ardavan Oskooi <ardavano@google.com>
Date: Tue, 23 Nov 2021 18:53:57 -0800
Subject: [PATCH 8/8] fix order of frequency check

---
 python/visualization.py | 27 +++++++++++++--------------
 1 file changed, 13 insertions(+), 14 deletions(-)

diff --git a/python/visualization.py b/python/visualization.py
index 29e4167b9..1e17f6a08 100644
--- a/python/visualization.py
+++ b/python/visualization.py
@@ -334,25 +334,24 @@ def plot_eps(sim, ax, output_plane=None, eps_parameters=None, frequency=None):
     eps_parameters = default_eps_parameters if eps_parameters is None else dict(default_eps_parameters, **eps_parameters)
 
     # Determine a frequency to plot all epsilon
+    if frequency is not None:
+        warnings.warn('The frequency parameter of plot2D has been deprecated. Use the frequency key of the eps_parameters dictionary instead.')
+        eps_parameters['frequency'] = frequency
     if eps_parameters['frequency'] is None:
-        if frequency is None:
+        try:
+            # Continuous sources
+            eps_parameters['frequency'] = sim.sources[0].frequency
+        except:
             try:
-                # Continuous sources
-                eps_parameters['frequency'] = sim.sources[0].frequency
+                # Gaussian sources
+                eps_parameters['frequency'] = sim.sources[0].src.frequency
             except:
                 try:
-                    # Gaussian sources
-                    eps_parameters['frequency'] = sim.sources[0].src.frequency
+                    # Custom sources
+                    eps_parameters['frequency'] = sim.sources[0].src.center_frequency
                 except:
-                    try:
-                        # Custom sources
-                        eps_parameters['frequency'] = sim.sources[0].src.center_frequency
-                    except:
-                        # No sources
-                        eps_parameters['frequency'] = 0
-        else:
-            warnings.warn('The frequency parameter of plot2D has been deprecated. Use the frequency key of the eps_parameters dictionary instead.')
-            eps_parameters['frequency'] = frequency
+                    # No sources
+                    eps_parameters['frequency'] = 0
 
     # Get domain measurements
     sim_center, sim_size = get_2D_dimensions(sim, output_plane)