WIP: IIR Filter Susceptibilites

python/geom.py

@@ -306,6 +306,24 @@ def eval_susceptibility(self,freq):
         else:
             return self.frequency*self.frequency / (self.frequency*self.frequency - freq*freq - 1j*self.gamma*freq) * sigma
 
+class IIR_Susceptibility(Susceptibility):
+    def __init__(self, num, den, **kwargs):
+        super(IIR_Susceptibility, self).__init__(**kwargs)
+        self.num = num
+        self.den = den
+
+    def eval_susceptibility(self,freq):
+        sigma = np.expand_dims(Matrix(diag=self.sigma_diag,offdiag=self.sigma_offdiag),axis=0)
+
+        freq = np.squeeze(freq)
+        eval_num = np.zeros(freq.shape,dtype=np.complex128)
+        eval_den = np.zeros(freq.shape,dtype=np.complex128)
+        N = len(self.num)
+        D = len(self.den)
+        for k in range(N): eval_num += pow(1j*freq,N-k)*self.num[k]
+        for k in range(D): eval_den += pow(1j*freq,D-k)*self.den[k]
+
+        return sigma * eval_num[:,np.newaxis,np.newaxis] / eval_den[:,np.newaxis,np.newaxis]
 
 class DrudeSusceptibility(Susceptibility):
 

python/meep.i

@@ -1169,6 +1169,12 @@ meep::volume_list *make_volume_list(const meep::volume &v, int c,
     }
 }
 
+//--------------------------------------------------
+// IIR transfer function typemaps
+//--------------------------------------------------
+%apply (double *INPLACE_ARRAY1, int DIM1) {(double *r, int N)};
+%apply (double *INPLACE_ARRAY1, int DIM1) {(double *vec_numS, int N), (double *vec_denS, int D), (double *vec_numZ, int N_z), (double *vec_denZ, int D_z)};
+
 // For some reason SWIG needs the namespaced version too
 %apply material_type_list { meep_geom::material_type_list };
 
@@ -1398,6 +1404,7 @@ PyObject *_get_array_slice_dimensions(meep::fields *f, const meep::volume &where
         GyrotropicDrudeSusceptibility,
         GyrotropicLorentzianSusceptibility,
         GyrotropicSaturatedSusceptibility,
+        IIR_Susceptibility,
         Lattice,
         LorentzianSusceptibility,
         Matrix,

python/typemap_utils.cpp

@@ -336,6 +336,9 @@ static int py_susceptibility_to_susceptibility(PyObject *po, susceptibility_stru
   s->saturated_gyrotropy = false;
   s->transitions.resize(0);
   s->initial_populations.resize(0);
+  s->iir = false;
+  s->numS.resize(0);
+  s->denS.resize(0);
 
   if (PyObject_HasAttrString(po, "frequency")) {
     if (!get_attr_dbl(po, &s->frequency, "frequency")) { return 0; }
@@ -390,6 +393,33 @@ static int py_susceptibility_to_susceptibility(PyObject *po, susceptibility_stru
     s->drude = false;
   }
 
+  // ---------- IIR filter sus. typemaps --------------------//
+  if (PyObject_HasAttrString(po, "num")) {
+    s->iir = true;
+
+    // numerator python list to C++ vector
+    PyObject *py_pop = PyObject_GetAttrString(po, "num");
+    if (!py_pop) { return 0; }
+    int length = PyList_Size(py_pop);
+    s->numS.resize(length);
+
+    for (int i = 0; i < length; ++i) {
+      s->numS[i] = PyFloat_AsDouble(PyList_GetItem(py_pop, i));
+    }
+
+    // denominator python list to C++ vector
+    py_pop = PyObject_GetAttrString(po, "den");
+    if (!py_pop) { return 0; }
+    length = PyList_Size(py_pop);
+    s->denS.resize(length);
+
+    for (int i = 0; i < length; ++i) {
+      s->denS[i] = PyFloat_AsDouble(PyList_GetItem(py_pop, i));
+    }
+    Py_DECREF(py_pop); 
+
+  }
+
   s->is_file = false;
 
   return 1;

src/Makefile.am

@@ -12,7 +12,7 @@ bicgstab.hpp meepgeom.hpp material_data.hpp
 libmeep_la_SOURCES = array_slice.cpp anisotropic_averaging.cpp 		\
 bands.cpp boundaries.cpp bicgstab.cpp casimir.cpp 	\
 control_c.cpp cw_fields.cpp dft.cpp dft_ldos.cpp energy_and_flux.cpp 	\
-fields.cpp loop_in_chunks.cpp h5fields.cpp h5file.cpp 	\
+fields.cpp loop_in_chunks.cpp h5fields.cpp h5file.cpp iir-susceptibility.cpp 	\
 initialize.cpp integrate.cpp integrate2.cpp monitor.cpp mympi.cpp 	\
 multilevel-atom.cpp near2far.cpp output_directory.cpp random.cpp 	\
 sources.cpp step.cpp step_db.cpp stress.cpp structure.cpp structure_dump.cpp		\