add support for dispersive ε

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);
 }
 
 %}

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,

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:

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

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);