WIP: periodic green's functions (#769)

* first attempt to hack in periodic green's functions

* typo

* make nperiods an option for near2far

* only do lattice sum in non-empty dirs

* simplify kw handling in add_near2far

* fix symmetry again

python/simulation.py

@@ -321,7 +321,8 @@ class DftNear2Far(DftObj):
     def __init__(self, func, args):
         super(DftNear2Far, self).__init__(func, args)
         self.nfreqs = args[2]
-        self.regions = args[3]
+        self.nperiods = args[3]
+        self.regions = args[4]
         self.num_components = 4
 
     @property
@@ -1373,15 +1374,16 @@ def get_dft_data(self, dft_chunk):
         mp._get_dft_data(dft_chunk, arr)
         return arr
 
-    def add_near2far(self, fcen, df, nfreq, *near2fars):
-        n2f = DftNear2Far(self._add_near2far, [fcen, df, nfreq, near2fars])
+    def add_near2far(self, fcen, df, nfreq, *near2fars, **kwargs):
+        nperiods = kwargs.get('nperiods', 1)
+        n2f = DftNear2Far(self._add_near2far, [fcen, df, nfreq, nperiods, near2fars])
         self.dft_objects.append(n2f)
         return n2f
 
-    def _add_near2far(self, fcen, df, nfreq, near2fars):
+    def _add_near2far(self, fcen, df, nfreq, nperiods, near2fars):
         if self.fields is None:
             self.init_sim()
-        return self._add_fluxish_stuff(self.fields.add_dft_near2far, fcen, df, nfreq, near2fars)
+        return self._add_fluxish_stuff(self.fields.add_dft_near2far, fcen, df, nfreq, near2fars, nperiods)
 
     def add_energy(self, fcen, df, nfreq, *energys):
         en = DftEnergy(self._add_energy, [fcen, df, nfreq, energys])
@@ -1636,7 +1638,7 @@ def solve_cw(self, tol=1e-8, maxiters=10000, L=2):
         self._evaluate_dft_objects()
         return self.fields.solve_cw(tol, maxiters, L)
 
-    def _add_fluxish_stuff(self, add_dft_stuff, fcen, df, nfreq, stufflist):
+    def _add_fluxish_stuff(self, add_dft_stuff, fcen, df, nfreq, stufflist, *args):
         vol_list = None
 
         for s in stufflist:
@@ -1649,7 +1651,7 @@ def _add_fluxish_stuff(self, add_dft_stuff, fcen, df, nfreq, stufflist):
                         is_cylindrical=self.is_cylindrical).swigobj
             vol_list = mp.make_volume_list(v2, c, s.weight, vol_list)
 
-        stuff = add_dft_stuff(vol_list, fcen - df / 2, fcen + df / 2, nfreq)
+        stuff = add_dft_stuff(vol_list, fcen - df / 2, fcen + df / 2, nfreq, *args)
         vol_list.__swig_destroy__(vol_list)
 
         return stuff

src/meep.hpp

@@ -1117,7 +1117,8 @@ class dft_near2far {
   /* fourier tranforms of tangential E and H field components in a
      medium with the given scalar eps and mu */
   dft_near2far(dft_chunk *F, double fmin, double fmax, int Nf, double eps, double mu,
-               const volume &where_);
+               const volume &where_, const direction periodic_d_[2],
+               const int periodic_n_[2], const double periodic_k_[2], const double period_[2]);
   dft_near2far(const dft_near2far &f);
 
   /* return an array (Ex,Ey,Ez,Hx,Hy,Hz) x Nfreq of the far fields at x */
@@ -1157,6 +1158,9 @@ class dft_near2far {
   dft_chunk *F;
   double eps, mu;
   volume where;
+  direction periodic_d[2];
+  int periodic_n[2];
+  double periodic_k[2], period[2];
 };
 
 /* Class to compute local-density-of-states spectra: the power spectrum
@@ -1729,7 +1733,7 @@ class fields {
 
   // near2far.cpp
   dft_near2far add_dft_near2far(const volume_list *where, double freq_min, double freq_max,
-                                int Nfreq);
+                                int Nfreq, int Nperiods=1);
   // monitor.cpp
   double get_chi1inv(component, direction, const vec &loc) const;
   double get_inveps(component c, direction d, const vec &loc) const {

src/near2far.cpp

@@ -30,16 +30,30 @@ using namespace std;
 namespace meep {
 
 dft_near2far::dft_near2far(dft_chunk *F_, double fmin, double fmax, int Nf, double eps_, double mu_,
-                           const volume &where_)
+                           const volume &where_, const direction periodic_d_[2],
+                           const int periodic_n_[2], const double periodic_k_[2], const double period_[2])
     : Nfreq(Nf), F(F_), eps(eps_), mu(mu_), where(where_) {
   if (Nf <= 1) fmin = fmax = (fmin + fmax) * 0.5;
   freq_min = fmin;
   dfreq = Nf <= 1 ? 0.0 : (fmax - fmin) / (Nf - 1);
+  for (int i = 0; i < 2; ++i) {
+    periodic_d[i] = periodic_d_[i];
+    periodic_n[i] = periodic_n_[i];
+    periodic_k[i] = periodic_k_[i];
+    period[i] = period_[i];
+  }
 }
 
 dft_near2far::dft_near2far(const dft_near2far &f)
     : freq_min(f.freq_min), dfreq(f.dfreq), Nfreq(f.Nfreq), F(f.F), eps(f.eps), mu(f.mu),
-      where(f.where) {}
+      where(f.where) {
+  for (int i = 0; i < 2; ++i) {
+    periodic_d[i] = f.periodic_d[i];
+    periodic_n[i] = f.periodic_n[i];
+    periodic_k[i] = f.periodic_k[i];
+    period[i] = f.period[i];
+  }
+}
 
 void dft_near2far::remove() {
   while (F) {
@@ -241,9 +255,21 @@ void dft_near2far::farfield_lowlevel(std::complex<double> *EH, const vec &x) {
       x0 = f->S.transform(x0, f->sn) + rshift;
       for (int i = 0; i < Nfreq; ++i) {
         double freq = freq_min + i * dfreq;
-        green(EH6, x, freq, eps, mu, x0, c0, f->dft[Nfreq * idx_dft + i]);
-        for (int j = 0; j < 6; ++j)
-          EH[i * 6 + j] += EH6[j];
+        vec xs(x0);
+        for (int i0 = -periodic_n[0]; i0 <= periodic_n[0]; ++i0) {
+          if (periodic_d[0] != NO_DIRECTION)
+            xs.set_direction(periodic_d[0], x0.in_direction(periodic_d[0]) + i0*period[0]);
+          double phase0 = i0*periodic_k[0];
+          for (int i1 = -periodic_n[1]; i1 <= periodic_n[1]; ++i1) {
+            if (periodic_d[1] != NO_DIRECTION)
+              xs.set_direction(periodic_d[1], x0.in_direction(periodic_d[1]) + i1*period[1]);
+            double phase = phase0 + i1*periodic_k[1];
+            std::complex<double> cphase = std::polar(1.0, phase);
+            green(EH6, x, freq, eps, mu, xs, c0, f->dft[Nfreq * idx_dft + i]);
+            for (int j = 0; j < 6; ++j)
+              EH[i * 6 + j] += EH6[j] * cphase;
+          }
+        }
       }
       idx_dft++;
     }
@@ -418,11 +444,15 @@ double *dft_near2far::flux(direction df, const volume &where, double resolution)
 static double approxeq(double a, double b) { return fabs(a - b) < 0.5e-11 * (fabs(a) + fabs(b)); }
 
 dft_near2far fields::add_dft_near2far(const volume_list *where, double freq_min, double freq_max,
-                                      int Nfreq) {
+                                      int Nfreq, int Nperiods) {
   dft_chunk *F = 0; /* E and H chunks*/
   double eps = 0, mu = 0;
   volume everywhere = where->v;
 
+  direction periodic_d[2] = {NO_DIRECTION, NO_DIRECTION};
+  int periodic_n[2] = {0, 0};
+  double periodic_k[2] = {0, 0}, period[2] = {0, 0};
+
   for (const volume_list *w = where; w; w = w->next) {
     everywhere = everywhere | where->v;
     direction nd = component_direction(w->c);
@@ -465,6 +495,20 @@ dft_near2far fields::add_dft_near2far(const volume_list *where, double freq_min,
       default: abort("invalid normal direction in dft_near2far!");
     }
 
+    if (Nperiods > 1) {
+      for (int i = 0; i < 2; ++i) {
+        double user_width = user_volume.num_direction(fd[i]) / a;
+        if (has_direction(v.dim, fd[i]) &&
+            boundaries[High][fd[i]] == Periodic && boundaries[Low][fd[i]] == Periodic &&
+            float(w->v.in_direction(fd[i])) >= float(user_width)) {
+          periodic_d[i] = fd[i];
+          periodic_n[i] = Nperiods;
+          period[i] = user_width;
+          periodic_k[i] = 2*pi*real(k[fd[i]]) * period[i];
+        }
+      }
+    }
+
     for (int i = 0; i < 2; ++i) {   /* E or H */
       for (int j = 0; j < 2; ++j) { /* first or second component */
         component c = direction_component(i == 0 ? Ex : Hx, fd[j]);
@@ -480,7 +524,8 @@ dft_near2far fields::add_dft_near2far(const volume_list *where, double freq_min,
     }
   }
 
-  return dft_near2far(F, freq_min, freq_max, Nfreq, eps, mu, everywhere);
+  return dft_near2far(F, freq_min, freq_max, Nfreq, eps, mu, everywhere,
+                      periodic_d, periodic_n, periodic_k, period);
 }
 
 } // namespace meep