WIP: Allow unequal DFT frequencies

python/meep.i

@@ -965,6 +965,25 @@ meep::volume_list *make_volume_list(const meep::volume &v, int c,
 //--------------------------------------------------
 // typemaps needed for add_dft_fields
 //--------------------------------------------------
+/*
+%typecheck(SWIG_TYPECHECK_POINTER, fragment="NumPy_Fragments") double* freqs {
+    $1 = is_array($input);
+}
+
+%typemap(in, fragment="NumPy_Macros") double* freqs {
+    $1 = (double *)array_data($input);
+}
+
+%typecheck(SWIG_TYPECHECK_INTEGER) int Nfreqs {
+    $1 = PyInteger_Check($input);
+}
+
+%typemap(in) int Nfreqs {
+  int py_nfreqs = PyInteger_Check($input);
+  $1 = py_nfreqs;
+}
+*/
+%apply (double *IN_ARRAY1, int DIM1) {(double *freqs, int Nfreqs)};
 
 %typecheck(SWIG_TYPECHECK_POINTER) const volume where {
     int py_material = PyObject_IsInstance($input, py_volume_object());
@@ -1346,11 +1365,8 @@ void _get_eigenmode(meep::fields *f, double omega_src, meep::direction d, const
 // Make omega members of meep::dft_ldos available as 'freq' in python
 %extend meep::dft_ldos {
 
-    double get_omega_min() {
-        return $self->omega_min;
-    }
-    double get_domega() {
-        return $self->domega;
+    double* get_omegas() {
+        return $self->omegas;
     }
     int get_Nomega() {
         return $self->Nomega;
@@ -1363,12 +1379,10 @@ void _get_eigenmode(meep::fields *f, double omega_src, meep::direction d, const
             stop = start + (self.domega / (2 * math.pi)) * self.Nomega
             return np.linspace(start, stop, num=self.Nomega, endpoint=False).tolist()
 
-        __swig_getmethods__["freq_min"] = get_omega_min
+        __swig_getmethods__["freqs"] = get_omegas
         __swig_getmethods__["nfreq"] = get_Nomega
-        __swig_getmethods__["dfreq"] = get_domega
-        if _newclass: freq_min = property(get_omega_min)
+        if _newclass: freqs = property(get_omegas)
         if _newclass: nfreq = property(get_Nomega)
-        if _newclass: dfreq = property(get_domega)
     %}
 }
 

python/simulation.py

@@ -316,12 +316,8 @@ def remove(self):
         return self.swigobj_attr('remove')
 
     @property
-    def freq_min(self):
-        return self.swigobj_attr('freq_min')
-
-    @property
-    def dfreq(self):
-        return self.swigobj_attr('dfreq')
+    def freqs(self):
+        return self.swigobj_attr('freqs')
 
     @property
     def Nfreq(self):
@@ -336,8 +332,8 @@ class DftFlux(DftObj):
 
     def __init__(self, func, args):
         super(DftFlux, self).__init__(func, args)
-        self.nfreqs = args[2]
-        self.regions = args[3]
+        self.nfreqs = args[1]
+        self.regions = args[2]
         self.num_components = 4
 
     @property
@@ -721,8 +717,8 @@ def _check_material_frequencies(self):
 
         dft_freqs = []
         for dftf in self.dft_objects:
-            dft_freqs.append(dftf.freq_min)
-            dft_freqs.append(dftf.freq_min + dftf.Nfreq * dftf.dfreq)
+            dft_freqs.append(dftf.freqs[0])
+            dft_freqs.append(dftf.freqs[-1])
 
         warn_src = ('Note: your sources include frequencies outside the range of validity of the ' +
                     'material models. This is fine as long as you eventually only look at outputs ' +
@@ -1630,22 +1626,28 @@ def _evaluate_dft_objects(self):
             if dft.swigobj is None:
                 dft.swigobj = dft.func(*dft.args)
 
-    def add_dft_fields(self, components, freq_min, freq_max, nfreq, where=None, center=None, size=None, yee_grid=False):
+    # TODO: make add_dft_fields_uneven, wrap add_dft_fields around it, change parameters of _add_dft_fields
+    def add_dft_fields_uneven(self, components, freqs, nfreq, where=None, center=None, size=None, yee_grid=False):
+        freqs.sort()
         center_v3 = Vector3(*center) if center is not None else None
         size_v3 = Vector3(*size) if size is not None else None
         use_centered_grid = not yee_grid
-        dftf = DftFields(self._add_dft_fields, [components, where, center_v3, size_v3, freq_min, freq_max, nfreq, use_centered_grid])
+        dftf = DftFields(self._add_dft_fields, [components, where, center_v3, size_v3, freqs, nfreq, use_centered_grid])
         self.dft_objects.append(dftf)
         return dftf
 
-    def _add_dft_fields(self, components, where, center, size, freq_min, freq_max, nfreq, use_centered_grid):
+    def add_dft_fields(self, components, freq_min, freq_max, nfreq, where=None, center=None, size=None, yee_grid=False):
+        freqs = np.linspace(start=freq_min, end=freq_max, num=nfreq)
+        return self.add_dft_fields_uneven(components, freqs, nfreq, where, center, size, yee_grid)
+
+    def _add_dft_fields(self, components, where, center, size, freqs, nfreq, use_centered_grid):
         if self.fields is None:
             self.init_sim()
         try:
             where = self._volume_from_kwargs(where, center, size)
         except ValueError:
             where = self.fields.total_volume()
-        return self.fields.add_dft_fields(components, where, freq_min, freq_max, nfreq, use_centered_grid)
+        return self.fields.add_dft_fields(components, where, freqs, nfreq, use_centered_grid)
 
     def output_dft(self, dft_fields, fname):
         if self.fields is None:
@@ -1664,27 +1666,40 @@ def get_dft_data(self, dft_chunk):
         mp._get_dft_data(dft_chunk, arr)
         return arr
 
-    def add_near2far(self, fcen, df, nfreq, *near2fars, **kwargs):
+    # TODO: make add_near2far_uneven, wrap add_near2far around it, change parameters of _add_near2far
+    def add_near2far_uneven(self, freqs, nfreq, *near2fars, **kwargs):
+        freqs.sort()
         nperiods = kwargs.get('nperiods', 1)
-        n2f = DftNear2Far(self._add_near2far, [fcen, df, nfreq, nperiods, near2fars])
+        n2f = DftNear2Far(self._add_near2far, [freqs, nfreq, nperiods, near2fars])
         self.dft_objects.append(n2f)
         return n2f
 
-    def _add_near2far(self, fcen, df, nfreq, nperiods, near2fars):
+    def add_near2far(self, fcen, df, nfreq, *near2fars, **kwargs):
+        freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq)
+        return self.add_near2far_uneven(freqs, nfreq, *near2fars, **kwargs)
+
+    def _add_near2far(self, freqs, 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, nperiods)
+        return self._add_fluxish_stuff(self.fields.add_dft_near2far, freqs, nfreq, near2fars, nperiods)
 
-    def add_energy(self, fcen, df, nfreq, *energys):
-        en = DftEnergy(self._add_energy, [fcen, df, nfreq, energys])
+    # TODO: make add_energy_uneven, wrap add_energy around it, change parameters of _add_energy
+    def add_energy_uneven(self, freqs, nfreq, *energys):
+        freqs.sort()
+        en = DftEnergy(self._add_energy, [freqs, nfreq, energys])
         self.dft_objects.append(en)
         return en
 
-    def _add_energy(self, fcen, df, nfreq, energys):
+    def add_energy(self, fcen, df, nfreq, *energys):
+        freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq)
+        return self.add_energy_uneven(freqs, nfreq, *energys)
+
+    def _add_energy(self, freqs, nfreq, energys):
         if self.fields is None:
             self.init_sim()
-        return self._add_fluxish_stuff(self.fields.add_dft_energy, fcen, df, nfreq, energys)
+        return self._add_fluxish_stuff(self.fields.add_dft_energy, freqs, nfreq, energys)
 
+    # TODO: check if this needs to be changed
     def _display_energy(self, name, func, energys):
         if energys:
             freqs = get_energy_freqs(energys[0])
@@ -1770,16 +1785,23 @@ def load_minus_near2far_data(self, n2f, n2fdata):
         self.load_near2far_data(n2f, n2fdata)
         n2f.scale_dfts(complex(-1.0))
 
-    def add_force(self, fcen, df, nfreq, *forces):
-        force = DftForce(self._add_force, [fcen, df, nfreq, forces])
+    # TODO: make add_force_uneven, wrap add_force around it, change parameters of _add_force
+    def add_force_uneven(self, freqs, nfreq, *forces):
+        freqs.sort()
+        force = DftForce(self._add_force, [freqs, nfreq, forces])
         self.dft_objects.append(force)
         return force
 
-    def _add_force(self, fcen, df, nfreq, forces):
+    def add_force(self, fcen, df, nfreq, *forces):
+        freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq)
+        return self.add_force_uneven(freqs, nfreq, *forces)
+
+    def _add_force(self, freqs, nfreq, forces):
         if self.fields is None:
             self.init_sim()
-        return self._add_fluxish_stuff(self.fields.add_dft_force, fcen, df, nfreq, forces)
+        return self._add_fluxish_stuff(self.fields.add_dft_force, freqs, nfreq, forces)
 
+    # TODO: check if display_forces needs to be changed
     def display_forces(self, *forces):
         force_freqs = get_force_freqs(forces[0])
         display_csv(self, 'force', zip(force_freqs, *[get_forces(f) for f in forces]))
@@ -1812,22 +1834,34 @@ def load_minus_force_data(self, force, fdata):
         self.load_force_data(force, fdata)
         force.scale_dfts(complex(-1.0))
 
-    def add_flux(self, fcen, df, nfreq, *fluxes):
-        flux = DftFlux(self._add_flux, [fcen, df, nfreq, fluxes])
+    # TODO: make add_flux_uneven, wrap add_flux around it, change parameters of _add_flux
+    def add_flux_uneven(self, freqs, nfreq, *fluxes):
+        freqs.sort()
+        flux = DftFlux(self._add_flux, [freqs, nfreq, fluxes])
         self.dft_objects.append(flux)
         return flux
 
-    def _add_flux(self, fcen, df, nfreq, fluxes):
+    def add_flux(self, fcen, df, nfreq, *fluxes):
+        freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq)
+        return self.add_flux_uneven(freqs, nfreq, *fluxes)
+
+    def _add_flux(self, freqs, nfreq, fluxes):
         if self.fields is None:
             self.init_sim()
-        return self._add_fluxish_stuff(self.fields.add_dft_flux, fcen, df, nfreq, fluxes)
+        return self._add_fluxish_stuff(self.fields.add_dft_flux, freqs, nfreq, fluxes)
 
-    def add_mode_monitor(self, fcen, df, nfreq, *fluxes):
-        flux = DftFlux(self._add_mode_monitor, [fcen, df, nfreq, fluxes])
+    # TODO: make add_mode_monitor_uneven, wrap add_mode_monitor around it, change parameters of _add_mode_monitor
+    def add_mode_monitor_uneven(self, freqs, nfreq, *fluxes):
+        freqs.sort()
+        flux = DftFlux(self._add_mode_monitor, [freqs, nfreq, fluxes])
         self.dft_objects.append(flux)
         return flux
 
-    def _add_mode_monitor(self, fcen, df, nfreq, fluxes):
+    def add_mode_monitor(self, fcen, df, nfreq, *fluxes):
+        freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq)
+        return self.add_mode_monitor_uneven(freqs, nfreq, *fluxes)
+
+    def _add_mode_monitor(self, freqs, nfreq, fluxes):
         if self.fields is None:
             self.init_sim()
 
@@ -1839,7 +1873,7 @@ def _add_mode_monitor(self, fcen, df, nfreq, fluxes):
         d0 = region.direction
         d = self.fields.normal_direction(v.swigobj) if d0 < 0 else d0
 
-        return self.fields.add_mode_monitor(d, v.swigobj, fcen - df / 2, fcen + df / 2, nfreq)
+        return self.fields.add_mode_monitor(d, v.swigobj, freqs, nfreq)
 
     def add_eigenmode(self, fcen, df, nfreq, *fluxes):
         warnings.warn('add_eigenmode is deprecated. Please use add_mode_monitor instead.', DeprecationWarning)
@@ -1936,7 +1970,8 @@ 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, *args):
+    # TODO: change parameters of _add_fluxish_stuff
+    def _add_fluxish_stuff(self, add_dft_stuff, freqs, nfreq, stufflist, *args):
         vol_list = None
 
         for s in stufflist:
@@ -1949,7 +1984,7 @@ def _add_fluxish_stuff(self, add_dft_stuff, fcen, df, nfreq, stufflist, *args):
                         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, *args)
+        stuff = add_dft_stuff(vol_list, freqs, nfreq, *args)
         vol_list.__swig_destroy__(vol_list)
 
         return stuff
@@ -2945,11 +2980,13 @@ def Ldos(fcen, df, nfreq):
     return mp._dft_ldos(fcen - df / 2, fcen + df / 2, nfreq)
 
 
-def dft_ldos(fcen=None, df=None, nfreq=None, ldos=None):
+# TODO: make dft_ldos_uneven, wrap dft_ldos wrap around it
+def dft_ldos_uneven(freqs=None, nfreq=None, ldos=None):
     if ldos is None:
-        if fcen is None or df is None or nfreq is None:
-            raise ValueError("Either fcen, df, and nfreq, or an Ldos is required for dft_ldos")
-        ldos = mp._dft_ldos(fcen - df / 2, fcen + df / 2, nfreq)
+        if freqs is None or nfreq is None:
+            raise ValueError("Either freqs and nfreq or an Ldos is required for dft_ldos_uneven")
+        freqs.sort()
+        ldos = mp._dft_ldos(freqs, nfreq)
 
     def _ldos(sim, todo):
         if todo == 'step':
@@ -2962,12 +2999,21 @@ def _ldos(sim, todo):
     return _ldos
 
 
+def dft_ldos(fcen=None, df=None, nfreq=None, ldos=None):
+    if ldos is None:
+        if fcen is None or df is None or nfreq is None:
+            raise ValueError("Either fcen, df, and nfreq, or an Ldos is required for dft_ldos")
+        freqs = np.linspace(fcen - df / 2, fcen + df / 2, nfreq)
+        return dft_ldos_uneven(freqs, nfreq, ldos)
+
+
 def scale_flux_fields(s, flux):
     flux.scale_dfts(s)
 
 
+# TODO: change this
 def get_flux_freqs(f):
-    return np.linspace(f.freq_min, f.freq_min + f.dfreq * f.Nfreq, num=f.Nfreq, endpoint=False).tolist()
+    return f.freqs
 
 
 def get_fluxes(f):
@@ -2978,12 +3024,14 @@ def scale_force_fields(s, force):
     force.scale_dfts(s)
 
 
+# TODO: change this
 def get_eigenmode_freqs(f):
-    return np.linspace(f.freq_min, f.freq_min + f.dfreq * f.Nfreq, num=f.Nfreq, endpoint=False).tolist()
+    return f.freqs
 
 
+# TODO: change this
 def get_force_freqs(f):
-    return np.linspace(f.freq_min, f.freq_min + f.dfreq * f.Nfreq, num=f.Nfreq, endpoint=False).tolist()
+    return f.freqs
 
 
 def get_forces(f):
@@ -2994,16 +3042,19 @@ def scale_near2far_fields(s, n2f):
     n2f.scale_dfts(s)
 
 
+# TODO: change this
 def get_near2far_freqs(f):
-    return np.linspace(f.freq_min, f.freq_min + f.dfreq * f.Nfreq, num=f.Nfreq, endpoint=False).tolist()
+    return f.freqs
 
 
 def scale_energy_fields(s, ef):
-    df.scale_dfts(s)
+    # TODO: check if this is a typo
+    ef.scale_dfts(s)
 
 
+# TODO: change this
 def get_energy_freqs(f):
-    return np.linspace(f.freq_min, f.freq_min + f.dfreq * f.Nfreq, num=f.Nfreq, endpoint=False).tolist()
+    return f.freqs
 
 
 def get_electric_energy(f):