diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index 2f758ca2d0395e76634d01e5c0500a5374d14203..70a4426d6f81c0610aff559b7cb8620194f74a0a 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -119,80 +119,76 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
         self.fluid_includes += '#include "fftw_tools.hpp"\n'
         self.stat_src += """
                 //begincpp
+                double *velocity_moments  = new double[10*4];
+                double *vorticity_moments = new double[10*4];
+                ptrdiff_t *hist_velocity  = new ptrdiff_t[histogram_bins*4];
+                ptrdiff_t *hist_vorticity = new ptrdiff_t[histogram_bins*4];
+                double max_estimates[4];
+                fs->compute_velocity(fs->cvorticity);
+                double *spec_velocity  = new double[fs->nshells*9];
+                double *spec_vorticity = new double[fs->nshells*9];
+                fs->cospectrum(fs->cvelocity, fs->cvelocity, spec_velocity);
+                fs->cospectrum(fs->cvorticity, fs->cvorticity, spec_vorticity);
                 //endcpp
                 """
-        self.stat_src += """
-                //begincpp
-                    double *velocity_moments  = new double[10*4];
-                    double *vorticity_moments = new double[10*4];
-                    ptrdiff_t *hist_velocity  = new ptrdiff_t[histogram_bins*4];
-                    ptrdiff_t *hist_vorticity = new ptrdiff_t[histogram_bins*4];
-                    double max_estimates[4];
-                    fs->compute_velocity(fs->cvorticity);
-                    double *spec_velocity  = new double[fs->nshells*9];
-                    double *spec_vorticity = new double[fs->nshells*9];
-                    fs->cospectrum(fs->cvelocity, fs->cvelocity, spec_velocity);
-                    fs->cospectrum(fs->cvorticity, fs->cvorticity, spec_vorticity);
-                    //endcpp
-                    """
         if self.QR_stats_on:
             self.stat_src += """
                 //begincpp
-                    double *trS2_Q_R_moments  = new double[10*3];
-                    double *gradu_moments     = new double[10*9];
-                    ptrdiff_t *hist_trS2_Q_R  = new ptrdiff_t[histogram_bins*3];
-                    ptrdiff_t *hist_gradu     = new ptrdiff_t[histogram_bins*9];
-                    ptrdiff_t *hist_QR2D      = new ptrdiff_t[QR2D_histogram_bins*QR2D_histogram_bins];
-                    double trS2QR_max_estimates[3];
-                    double gradu_max_estimates[9];
-                    trS2QR_max_estimates[0] = max_trS2_estimate;
-                    trS2QR_max_estimates[1] = max_Q_estimate;
-                    trS2QR_max_estimates[2] = max_R_estimate;
-                    std::fill_n(gradu_max_estimates, 9, sqrt(3*max_trS2_estimate));
-                    fs->compute_gradient_statistics(
-                        fs->cvelocity,
-                        gradu_moments,
-                        trS2_Q_R_moments,
-                        hist_gradu,
-                        hist_trS2_Q_R,
-                        hist_QR2D,
-                        trS2QR_max_estimates,
-                        gradu_max_estimates,
-                        histogram_bins,
-                        QR2D_histogram_bins);
-                    //endcpp
-                    """
+                double *trS2_Q_R_moments  = new double[10*3];
+                double *gradu_moments     = new double[10*9];
+                ptrdiff_t *hist_trS2_Q_R  = new ptrdiff_t[histogram_bins*3];
+                ptrdiff_t *hist_gradu     = new ptrdiff_t[histogram_bins*9];
+                ptrdiff_t *hist_QR2D      = new ptrdiff_t[QR2D_histogram_bins*QR2D_histogram_bins];
+                double trS2QR_max_estimates[3];
+                double gradu_max_estimates[9];
+                trS2QR_max_estimates[0] = max_trS2_estimate;
+                trS2QR_max_estimates[1] = max_Q_estimate;
+                trS2QR_max_estimates[2] = max_R_estimate;
+                std::fill_n(gradu_max_estimates, 9, sqrt(3*max_trS2_estimate));
+                fs->compute_gradient_statistics(
+                    fs->cvelocity,
+                    gradu_moments,
+                    trS2_Q_R_moments,
+                    hist_gradu,
+                    hist_trS2_Q_R,
+                    hist_QR2D,
+                    trS2QR_max_estimates,
+                    gradu_max_estimates,
+                    histogram_bins,
+                    QR2D_histogram_bins);
+                //endcpp
+                """
         self.stat_src += """
                 //begincpp
-                    fs->ift_velocity();
-                    max_estimates[0] = max_velocity_estimate/sqrt(3);
-                    max_estimates[1] = max_estimates[0];
-                    max_estimates[2] = max_estimates[0];
-                    max_estimates[3] = max_velocity_estimate;
-                    fs->compute_rspace_stats4(fs->rvelocity,
-                                             velocity_moments,
-                                             hist_velocity,
-                                             max_estimates,
-                                             histogram_bins);
-                    fs->ift_vorticity();
-                    max_estimates[0] = max_vorticity_estimate/sqrt(3);
-                    max_estimates[1] = max_estimates[0];
-                    max_estimates[2] = max_estimates[0];
-                    max_estimates[3] = max_vorticity_estimate;
-                    fs->compute_rspace_stats4(fs->rvorticity,
-                                             vorticity_moments,
-                                             hist_vorticity,
-                                             max_estimates,
-                                             histogram_bins);
-                    if (fs->cd->myrank == 0)
-                    {{
-                        hid_t Cdset, wspace, mspace;
-                        int ndims;
-                        hsize_t count[4], offset[4], dims[4];
-                        offset[0] = fs->iteration/niter_stat;
-                        offset[1] = 0;
-                        offset[2] = 0;
-                        offset[3] = 0;
+                fs->ift_velocity();
+                max_estimates[0] = max_velocity_estimate/sqrt(3);
+                max_estimates[1] = max_estimates[0];
+                max_estimates[2] = max_estimates[0];
+                max_estimates[3] = max_velocity_estimate;
+                fs->compute_rspace_stats4(fs->rvelocity,
+                                         velocity_moments,
+                                         hist_velocity,
+                                         max_estimates,
+                                         histogram_bins);
+                fs->ift_vorticity();
+                max_estimates[0] = max_vorticity_estimate/sqrt(3);
+                max_estimates[1] = max_estimates[0];
+                max_estimates[2] = max_estimates[0];
+                max_estimates[3] = max_vorticity_estimate;
+                fs->compute_rspace_stats4(fs->rvorticity,
+                                         vorticity_moments,
+                                         hist_vorticity,
+                                         max_estimates,
+                                         histogram_bins);
+                if (fs->cd->myrank == 0)
+                {{
+                    hid_t Cdset, wspace, mspace;
+                    int ndims;
+                    hsize_t count[4], offset[4], dims[4];
+                    offset[0] = fs->iteration/niter_stat;
+                    offset[1] = 0;
+                    offset[2] = 0;
+                    offset[3] = 0;
                 //endcpp
                 """.format(self.C_dtype)
         if self.dtype == np.float32:
@@ -300,23 +296,23 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
                         """)
         self.stat_src += """
                 //begincpp
-                    }
-                    delete[] spec_velocity;
-                    delete[] spec_vorticity;
-                    delete[] velocity_moments;
-                    delete[] vorticity_moments;
-                    delete[] hist_velocity;
-                    delete[] hist_vorticity;
+                }
+                delete[] spec_velocity;
+                delete[] spec_vorticity;
+                delete[] velocity_moments;
+                delete[] vorticity_moments;
+                delete[] hist_velocity;
+                delete[] hist_vorticity;
                 //endcpp
                 """
         if self.QR_stats_on:
             self.stat_src += """
                 //begincpp
-                    delete[] trS2_Q_R_moments;
-                    delete[] gradu_moments;
-                    delete[] hist_trS2_Q_R;
-                    delete[] hist_gradu;
-                    delete[] hist_QR2D;
+                delete[] trS2_Q_R_moments;
+                delete[] gradu_moments;
+                delete[] hist_trS2_Q_R;
+                delete[] hist_gradu;
+                delete[] hist_QR2D;
                 //endcpp
                 """
         return None
@@ -378,16 +374,27 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
                           '}\n' +
                           'delete fs;\n')
         return None
-    def add_particle_fields(
+    def add_3D_rFFTW_field(
+            self,
+            name = 'rFFTW_acc'):
+        if self.dtype == np.float32:
+            FFTW = 'fftwf'
+        elif self.dtype == np.float64:
+            FFTW = 'fftw'
+        self.fluid_variables += '{0} *{1};\n'.format(self.C_dtype, name)
+        self.fluid_start += '{0} = {1}_alloc_real(2*fs->cd->local_size);\n'.format(name, FFTW)
+        self.fluid_end   += '{0}_free({1});\n'.format(FFTW, name)
+        return None
+    def add_interpolator(
             self,
             interp_type = 'spline',
-            kcut = None,
             neighbours = 1,
             smoothness = 1,
-            name = 'particle_field',
-            field_class = 'rFFTW_interpolator'):
-        self.fluid_includes += '#include "{0}.hpp"\n'.format(field_class)
-        self.fluid_variables += field_class + '<{0}, {1}> *vel_{2}, *acc_{2};\n'.format(self.C_dtype, neighbours, name)
+            name = 'field_interpolator',
+            field_name = 'fs->rvelocity'):
+        self.fluid_includes += '#include "rFFTW_interpolator.hpp"\n'
+        self.fluid_variables += 'rFFTW_interpolator <{0}, {1}> *{2};\n'.format(
+                self.C_dtype, neighbours, name)
         self.parameters[name + '_type'] = interp_type
         self.parameters[name + '_neighbours'] = neighbours
         if interp_type == 'spline':
@@ -395,41 +402,66 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
             beta_name = 'beta_n{0}_m{1}'.format(neighbours, smoothness)
         elif interp_type == 'Lagrange':
             beta_name = 'beta_Lagrange_n{0}'.format(neighbours)
-        self.fluid_start += ('vel_{0} = new {1}<{2}, {3}>(fs, {4});\n' +
-                             'acc_{0} = new {1}<{2}, {3}>(fs, {4});\n').format(name,
-                                                                               field_class,
-                                                                               self.C_dtype,
-                                                                               neighbours,
-                                                                               beta_name)
-        self.fluid_end += ('delete vel_{0};\n' +
-                           'delete acc_{0};\n').format(name)
-        update_fields = 'fs->compute_velocity(fs->cvorticity);\n'
-        if not type(kcut) == type(None):
-            update_fields += 'fs->low_pass_Fourier(fs->cvelocity, 3, {0});\n'.format(kcut)
-        update_fields += ('fs->ift_velocity();\n' +
-                          'vel_{0}->read_rFFTW(fs->rvelocity);\n' +
-                          'fs->compute_Lagrangian_acceleration(acc_{0}->field);\n').format(name)
-        self.fluid_start += update_fields
-        self.fluid_loop += update_fields
+        self.fluid_start += '{0} = new rFFTW_interpolator<{1}, {2}>(fs, {3}, {4});\n'.format(
+                name,
+                self.C_dtype,
+                neighbours,
+                beta_name,
+                field_name)
+        self.fluid_end += 'delete {0};\n'.format(name)
         return None
     def add_particles(
             self,
-            integration_method = 'AdamsBashforth',
             integration_steps = 2,
-            kcut = 'fs->kM',
+            kcut = None,
+            interpolator = 'field_interpolator',
             frozen_particles = False,
-            fields_name = None,
-            particle_class = 'rFFTW_particles'):
-        if integration_method == 'cRK4':
-            integration_steps = 4
-        elif integration_method == 'Heun':
-            integration_steps = 2
-        neighbours = self.parameters[fields_name + '_neighbours']
-        self.parameters['tracers{0}_field'.format(self.particle_species)] = fields_name
-        self.parameters['tracers{0}_integration_method'.format(self.particle_species)] = integration_method
-        self.parameters['tracers{0}_kcut'.format(self.particle_species)] = kcut
-        self.parameters['tracers{0}_integration_steps'.format(self.particle_species)] = integration_steps
-        self.file_datasets_grow += """
+            acc_name = None):
+        """Adds code for tracking a series of particle species, each
+        consisting of `nparticles` particles.
+
+        :type integration_steps: int, list of int
+        :type kcut: None (default), str, list of str
+        :type interpolator: str, list of str
+        :type frozen_particles: bool
+        :type acc_name: str
+
+        .. warning :: if not None, kcut must be a list of decreasing
+        wavenumbers, since filtering is done on the same field...
+        """
+        if self.dtype == np.float32:
+            FFTW = 'fftwf'
+        elif self.dtype == np.float64:
+            FFTW = 'fftw'
+        s0 = self.particle_species
+        if type(integration_steps) == int:
+            integration_steps = [integration_steps]
+        if type(kcut) == str:
+            kcut = [kcut]
+        if type(interpolator) == str:
+            interpolator = [interpolator]
+        nspecies = max(len(integration_steps), len(interpolator))
+        if type(kcut) == list:
+            nspecies = max(nspecies, len(kcut))
+        if len(integration_steps) == 1:
+            integration_steps = [integration_steps[0] for s in range(nspecies)]
+        if len(interpolator) == 1:
+            interpolator = [interpolator[0] for s in range(nspecies)]
+        if type(kcut) == list:
+            if len(kcut) == 1:
+                kcut = [kcut[0] for s in range(nspecies)]
+        assert(len(integration_steps) == nspecies)
+        assert(len(interpolator) == nspecies)
+        if type(kcut) == list:
+            assert(len(kcut) == nspecies)
+        for s in range(nspecies):
+            neighbours = self.parameters[interpolator[s] + '_neighbours']
+            if type(kcut) == list:
+                self.parameters['tracers{0}_kcut'.format(s0 + s)] = kcut[s]
+            self.parameters['tracers{0}_acc_on'.format(s0 + s)] = not type(acc_name) == type(None)
+            self.parameters['tracers{0}_interpolator'.format(s0 + s)] = interpolator[s]
+            self.parameters['tracers{0}_integration_steps'.format(s0 + s)] = integration_steps[s]
+            self.file_datasets_grow += """
                         //begincpp
                         temp_string = (std::string("/particles/") +
                                        std::string(ps{0}->name));
@@ -437,110 +469,118 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
                         grow_particle_datasets(group, temp_string.c_str(), NULL, NULL);
                         H5Gclose(group);
                         //endcpp
-                        """.format(self.particle_species)
-        if self.dtype == np.float32:
-            FFTW = 'fftwf'
-        elif self.dtype == np.float64:
-            FFTW = 'fftw'
-        output_vel_acc =  """
-                          //begincpp
-                          {{
-                              double *acceleration = new double[ps{0}->array_size];
-                              double *velocity     = new double[ps{0}->array_size];
-                              ps{0}->sample_vec_field(vel_{1}, velocity);
-                              ps{0}->sample_vec_field(acc_{1}, acceleration);
-                              if (ps{0}->myrank == 0)
-                              {{
-                                  //VELOCITY begin
-                                  std::string temp_string = (std::string("/particles/") +
-                                                             std::string(ps{0}->name) +
-                                                             std::string("/velocity"));
-                                  hid_t Cdset = H5Dopen(stat_file, temp_string.c_str(), H5P_DEFAULT);
-                                  hid_t mspace, wspace;
-                                  int ndims;
-                                  hsize_t count[3], offset[3];
-                                  wspace = H5Dget_space(Cdset);
-                                  ndims = H5Sget_simple_extent_dims(wspace, count, NULL);
-                                  count[0] = 1;
-                                  offset[0] = ps{0}->iteration / ps{0}->traj_skip;
-                                  offset[1] = 0;
-                                  offset[2] = 0;
-                                  mspace = H5Screate_simple(ndims, count, NULL);
-                                  H5Sselect_hyperslab(wspace, H5S_SELECT_SET, offset, NULL, count, NULL);
-                                  H5Dwrite(Cdset, H5T_NATIVE_DOUBLE, mspace, wspace, H5P_DEFAULT, velocity);
-                                  H5Dclose(Cdset);
-                                  //VELOCITY end
-                                  //ACCELERATION begin
-                                  temp_string = (std::string("/particles/") +
-                                                 std::string(ps{0}->name) +
-                                                 std::string("/acceleration"));
-                                  Cdset = H5Dopen(stat_file, temp_string.c_str(), H5P_DEFAULT);
-                                  H5Dwrite(Cdset, H5T_NATIVE_DOUBLE, mspace, wspace, H5P_DEFAULT, acceleration);
-                                  H5Sclose(mspace);
-                                  H5Sclose(wspace);
-                                  H5Dclose(Cdset);
-                                  //ACCELERATION end
-                              }}
-                              delete[] acceleration;
-                              delete[] velocity;
-                          }}
-                          //endcpp
-                          """.format(self.particle_species, fields_name)
-        self.particle_start += 'sprintf(fname, "tracers{0}");\n'.format(self.particle_species)
-        self.particle_end += ('ps{0}->write(stat_file);\n' +
-                              'delete ps{0};\n').format(self.particle_species)
-        self.particle_includes += '#include "{0}.hpp"\n'.format(particle_class)
-        if particle_class == 'particles':
-            if integration_method == 'AdamsBashforth':
-                multistep = 'true'
-            else:
-                multistep = 'false'
-            self.particle_variables += 'particles<VELOCITY_TRACER, {0}, {1}, {2}> *ps{3};\n'.format(
-                    self.C_dtype,
-                    multistep,
-                    neighbours,
-                    self.particle_species)
-            self.particle_start += ('ps{0} = new particles<VELOCITY_TRACER, {1}, {2},{3}>(\n' +
-                                    'fname, fs, vel_{4},\n' +
-                                    'nparticles,\n' +
-                                    'niter_part, tracers{0}_integration_steps);\n').format(
-                                            self.particle_species, self.C_dtype, multistep, neighbours, fields_name)
+                        """.format(s0 + s)
+
+        #### code that outputs statistics
+        output_vel_acc = '{\n'
+        # array for putting sampled velocity in
+        # must compute velocity, just in case it was messed up by some
+        # other particle species before the stats
+        output_vel_acc += ('double *velocity = new double[3*nparticles];\n' +
+                           'fs->compute_velocity(fs->cvorticity);\n')
+        if not type(kcut) == list:
+            output_vel_acc += 'fs->ift_velocity();\n'
+        if not type(acc_name) == type(None):
+            # array for putting sampled acceleration in
+            # must compute acceleration
+            output_vel_acc += 'double *acceleration = new double[3*nparticles];\n'
+            output_vel_acc += 'fs->compute_Lagrangian_acceleration({0});\n'.format(acc_name)
+        for s in range(nspecies):
+            if type(kcut) == list:
+                output_vel_acc += 'fs->low_pass_Fourier(fs->cvelocity, 3, {0});\n'.format(kcut[s])
+                output_vel_acc += 'fs->ift_velocity();\n'
+            output_vel_acc += """
+                {0}->field = fs->rvelocity;
+                ps{1}->sample_vec_field({0}, velocity);
+                """.format(interpolator[s], s0 + s)
+            if not type(acc_name) == type(None):
+                output_vel_acc += """
+                    {0}->field = {1};
+                    ps{2}->sample_vec_field({0}, acceleration);
+                    """.format(interpolator[s], acc_name, s0 + s)
+            output_vel_acc += """
+                if (myrank == 0)
+                {{
+                //VELOCITY begin
+                std::string temp_string = (std::string("/particles/") +
+                                           std::string(ps{0}->name) +
+                                           std::string("/velocity"));
+                hid_t Cdset = H5Dopen(stat_file, temp_string.c_str(), H5P_DEFAULT);
+                hid_t mspace, wspace;
+                int ndims;
+                hsize_t count[3], offset[3];
+                wspace = H5Dget_space(Cdset);
+                ndims = H5Sget_simple_extent_dims(wspace, count, NULL);
+                count[0] = 1;
+                offset[0] = ps{0}->iteration / ps{0}->traj_skip;
+                offset[1] = 0;
+                offset[2] = 0;
+                mspace = H5Screate_simple(ndims, count, NULL);
+                H5Sselect_hyperslab(wspace, H5S_SELECT_SET, offset, NULL, count, NULL);
+                H5Dwrite(Cdset, H5T_NATIVE_DOUBLE, mspace, wspace, H5P_DEFAULT, velocity);
+                H5Dclose(Cdset);
+                //VELOCITY end\n""".format(s0 + s)
+            if not type(acc_name) == type(None):
+                output_vel_acc += """
+                    //ACCELERATION begin
+                    temp_string = (std::string("/particles/") +
+                                   std::string(ps{0}->name) +
+                                   std::string("/acceleration"));
+                    Cdset = H5Dopen(stat_file, temp_string.c_str(), H5P_DEFAULT);
+                    H5Dwrite(Cdset, H5T_NATIVE_DOUBLE, mspace, wspace, H5P_DEFAULT, acceleration);
+                    H5Sclose(mspace);
+                    H5Sclose(wspace);
+                    H5Dclose(Cdset);
+                    //ACCELERATION end\n""".format(s0 + s)
+            output_vel_acc += '}\n'
+        output_vel_acc += 'delete[] velocity;\n'
+        if not type(acc_name) == type(None):
+            output_vel_acc += 'delete[] acceleration;\n'
+        output_vel_acc += '}\n'
+
+        #### initialize, stepping and finalize code
+        if not type(kcut) == list:
+            update_fields = ('fs->compute_velocity(fs->cvorticity);\n' +
+                             'fs->ift_velocity();\n')
+            self.particle_start += update_fields
+            self.particle_loop  += update_fields
         else:
+            self.particle_loop += 'fs->compute_velocity(fs->cvorticity);\n'
+        self.particle_includes += '#include "rFFTW_particles.hpp"\n'
+        self.particle_stat_src += (
+                'if (ps0->iteration % niter_part == 0)\n' +
+                '{\n')
+        for s in range(nspecies):
+            self.particle_start += 'sprintf(fname, "tracers{0}");\n'.format(s0 + s)
+            self.particle_end += ('ps{0}->write(stat_file);\n' +
+                                  'delete ps{0};\n').format(s0 + s)
             self.particle_variables += 'rFFTW_particles<VELOCITY_TRACER, {0}, {1}> *ps{2};\n'.format(
                     self.C_dtype,
                     neighbours,
-                    self.particle_species)
+                    s0 + s)
             self.particle_start += ('ps{0} = new rFFTW_particles<VELOCITY_TRACER, {1}, {2}>(\n' +
-                                    'fname, vel_{3},\n' +
+                                    'fname, {3},\n' +
                                     'nparticles,\n' +
                                     'niter_part, tracers{0}_integration_steps);\n').format(
-                                            self.particle_species, self.C_dtype, neighbours, fields_name)
-        self.particle_start += ('ps{0}->dt = dt;\n' +
-                                'ps{0}->iteration = iteration;\n' +
-                                'ps{0}->read(stat_file);\n').format(self.particle_species)
+                                            s0 + s,
+                                            self.C_dtype,
+                                            neighbours,
+                                            interpolator[s])
+            self.particle_start += ('ps{0}->dt = dt;\n' +
+                                    'ps{0}->iteration = iteration;\n' +
+                                    'ps{0}->read(stat_file);\n').format(s0 + s)
+            if not frozen_particles:
+                if type(kcut) == list:
+                    update_field = ('fs->low_pass_Fourier(fs->cvelocity, 3, {0});\n'.format(kcut[s]) +
+                                    'fs->ift_velocity();\n')
+                    self.particle_loop += update_field
+                self.particle_loop += '{0}->field = fs->rvelocity;\n'.format(interpolator[s])
+                self.particle_loop += 'ps{0}->step();\n'.format(s0 + s)
+            self.particle_stat_src += 'ps{0}->write(stat_file, false);\n'.format(s0 + s)
         self.particle_start += output_vel_acc
-        if not frozen_particles:
-            if particle_class == 'particles':
-                if integration_method == 'AdamsBashforth':
-                    self.particle_loop += 'ps{0}->AdamsBashforth((ps{0}->iteration < ps{0}->integration_steps) ? ps{0}->iteration+1 : ps{0}->integration_steps);\n'.format(self.particle_species)
-                elif integration_method == 'Euler':
-                    self.particle_loop += 'ps{0}->Euler();\n'.format(self.particle_species)
-                elif integration_method == 'Heun':
-                    assert(integration_steps == 2)
-                    self.particle_loop += 'ps{0}->Heun();\n'.format(self.particle_species)
-                elif integration_method == 'cRK4':
-                    assert(integration_steps == 4)
-                    self.particle_loop += 'ps{0}->cRK4();\n'.format(self.particle_species)
-                self.particle_loop += 'ps{0}->iteration++;\n'.format(self.particle_species)
-                self.particle_loop += 'ps{0}->synchronize();\n'.format(self.particle_species)
-            elif particle_class == 'rFFTW_particles':
-                self.particle_loop += 'ps{0}->step();\n'.format(self.particle_species)
-        self.particle_stat_src += (
-                ('if (ps{0}->iteration % niter_part == 0)\n' +
-                 '{{\n' +
-                 'ps{0}->write(stat_file, false);\n').format(self.particle_species) +
-                output_vel_acc + '}\n')
-        self.particle_species += 1
+        self.particle_stat_src += output_vel_acc
+        self.particle_stat_src += '}\n'
+        self.particle_species += nspecies
         return None
     def get_data_file(self):
         return h5py.File(os.path.join(self.work_dir, self.simname + '.h5'), 'r')
@@ -553,7 +593,8 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
         with self.get_data_file() as data_file:
             if 'moments' not in data_file['statistics'].keys():
                 return None
-            iter0 = min(data_file['statistics/moments/velocity'].shape[0]*self.parameters['niter_stat']-1,
+            iter0 = min((data_file['statistics/moments/velocity'].shape[0] *
+                         self.parameters['niter_stat']-1),
                         iter0)
             if type(iter1) == type(None):
                 iter1 = data_file['iteration'].value
@@ -565,9 +606,9 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
             self.statistics['kM'] = data_file['kspace/kM'].value
             self.statistics['dk'] = data_file['kspace/dk'].value
             if self.particle_species > 0:
-                self.trajectories = [
-                        data_file['particles/' + key + '/state'][
-                            iter0//self.parameters['niter_part']:iter1//self.parameters['niter_part']+1]
+                self.trajectories = [data_file['particles/' + key + '/state'][
+                                        iter0//self.parameters['niter_part'] :
+                                        iter1//self.parameters['niter_part']+1]
                                      for key in data_file['particles'].keys()]
             computation_needed = True
             pp_file = h5py.File(self.get_postprocess_file_name(), 'a')
@@ -612,7 +653,8 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
             self.statistics[key + '(t)'] = (self.statistics['dk'] *
                                             np.sum(self.statistics[key + '(t, k)'], axis = 1))
         self.statistics['Uint(t)'] = np.sqrt(2*self.statistics['energy(t)'] / 3)
-        self.statistics['Lint(t)'] = ((self.statistics['dk']*np.pi / (2*self.statistics['Uint(t)']**2)) *
+        self.statistics['Lint(t)'] = ((self.statistics['dk']*np.pi /
+                                       (2*self.statistics['Uint(t)']**2)) *
                                       np.nansum(self.statistics['energy(t, k)'] /
                                                 self.statistics['kshell'][None, :], axis = 1))
         for key in ['energy',
@@ -639,7 +681,8 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
             self.statistics['Rlambda' + suffix] = (self.statistics['Uint' + suffix] *
                                                    self.statistics['lambda' + suffix] /
                                                    self.parameters['nu'])
-            self.statistics['kMeta' + suffix] = self.statistics['kM']*self.statistics['etaK' + suffix]
+            self.statistics['kMeta' + suffix] = (self.statistics['kM'] *
+                                                 self.statistics['etaK' + suffix])
             if self.parameters['dealias_type'] == 1:
                 self.statistics['kMeta' + suffix] *= 0.8
         self.statistics['Tint'] = self.statistics['Lint'] / self.statistics['Uint']
@@ -731,6 +774,51 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
                                      dtype = np.int64,
                                      compression = 'gzip')
         return None
+    def add_particle_fields(
+            self,
+            interp_type = 'spline',
+            kcut = None,
+            neighbours = 1,
+            smoothness = 1,
+            name = 'particle_field',
+            field_class = 'rFFTW_interpolator',
+            acc_field_name = 'rFFTW_acc'):
+        self.fluid_includes += '#include "{0}.hpp"\n'.format(field_class)
+        self.fluid_variables += field_class + '<{0}, {1}> *vel_{2}, *acc_{2};\n'.format(
+                self.C_dtype, neighbours, name)
+        self.parameters[name + '_type'] = interp_type
+        self.parameters[name + '_neighbours'] = neighbours
+        if interp_type == 'spline':
+            self.parameters[name + '_smoothness'] = smoothness
+            beta_name = 'beta_n{0}_m{1}'.format(neighbours, smoothness)
+        elif interp_type == 'Lagrange':
+            beta_name = 'beta_Lagrange_n{0}'.format(neighbours)
+        if field_class == 'rFFTW_interpolator':
+            self.fluid_start += ('vel_{0} = new {1}<{2}, {3}>(fs, {4}, fs->rvelocity);\n' +
+                                 'acc_{0} = new {1}<{2}, {3}>(fs, {4}, {5});\n').format(name,
+                                                                                   field_class,
+                                                                                   self.C_dtype,
+                                                                                   neighbours,
+                                                                                   beta_name,
+                                                                                   acc_field_name)
+        elif field_class == 'interpolator':
+            self.fluid_start += ('vel_{0} = new {1}<{2}, {3}>(fs, {4});\n' +
+                                 'acc_{0} = new {1}<{2}, {3}>(fs, {4});\n').format(name,
+                                                                                   field_class,
+                                                                                   self.C_dtype,
+                                                                                   neighbours,
+                                                                                   beta_name,
+                                                                                   acc_field_name)
+        self.fluid_end += ('delete vel_{0};\n' +
+                           'delete acc_{0};\n').format(name)
+        update_fields = 'fs->compute_velocity(fs->cvorticity);\n'
+        if not type(kcut) == type(None):
+            update_fields += 'fs->low_pass_Fourier(fs->cvelocity, 3, {0});\n'.format(kcut)
+        update_fields += ('fs->ift_velocity();\n' +
+                          'fs->compute_Lagrangian_acceleration(acc_{0}->field);\n').format(name)
+        self.fluid_start += update_fields
+        self.fluid_loop += update_fields
+        return None
 
 def launch(
         opt,
diff --git a/bfps/cpp/rFFTW_interpolator.cpp b/bfps/cpp/rFFTW_interpolator.cpp
index a858eb104abf19837f7ccd30ae673caaf841a493..c701fbd64bba8e1c68ab2de953eae0795baead2f 100644
--- a/bfps/cpp/rFFTW_interpolator.cpp
+++ b/bfps/cpp/rFFTW_interpolator.cpp
@@ -32,15 +32,13 @@
 template <class rnumber, int interp_neighbours>
 rFFTW_interpolator<rnumber, interp_neighbours>::rFFTW_interpolator(
         fluid_solver_base<rnumber> *fs,
-        base_polynomial_values BETA_POLYS)
+        base_polynomial_values BETA_POLYS,
+        rnumber *FIELD)
 {
     this->descriptor = fs->rd;
     this->field_size = 2*fs->cd->local_size;
     this->compute_beta = BETA_POLYS;
-    if (sizeof(rnumber) == 4)
-        this->field = (rnumber*)((void*)fftwf_alloc_real(this->field_size));
-    else if (sizeof(rnumber) == 8)
-        this->field = (rnumber*)((void*)fftw_alloc_real(this->field_size));
+    this->field = FIELD;
 
     // compute dx, dy, dz;
     this->dx = 4*acos(0) / (fs->dkx*this->descriptor->sizes[2]);
@@ -59,24 +57,9 @@ rFFTW_interpolator<rnumber, interp_neighbours>::rFFTW_interpolator(
 template <class rnumber, int interp_neighbours>
 rFFTW_interpolator<rnumber, interp_neighbours>::~rFFTW_interpolator()
 {
-    if (sizeof(rnumber) == 4)
-        fftwf_free((float*)((void*)this->field));
-    else if (sizeof(rnumber) == 8)
-        fftw_free((double*)((void*)this->field));
     delete[] this->compute;
 }
 
-template <class rnumber, int interp_neighbours>
-int rFFTW_interpolator<rnumber, interp_neighbours>::read_rFFTW(void *void_src)
-{
-    rnumber *src = (rnumber*)void_src;
-    /* do big copy of middle stuff */
-    std::copy(src,
-              src + this->field_size,
-              this->field);
-    return EXIT_SUCCESS;
-}
-
 template <class rnumber, int interp_neighbours>
 void rFFTW_interpolator<rnumber, interp_neighbours>::get_grid_coordinates(
         const int nparticles,
diff --git a/bfps/cpp/rFFTW_interpolator.hpp b/bfps/cpp/rFFTW_interpolator.hpp
index 154296018926202bb53367e4313510079343c4b8..3dccd7150d17bea612128bf6d70b37033d281db8 100644
--- a/bfps/cpp/rFFTW_interpolator.hpp
+++ b/bfps/cpp/rFFTW_interpolator.hpp
@@ -60,7 +60,8 @@ class rFFTW_interpolator
 
         rFFTW_interpolator(
                 fluid_solver_base<rnumber> *FSOLVER,
-                base_polynomial_values BETA_POLYS);
+                base_polynomial_values BETA_POLYS,
+                rnumber *FIELD_DATA);
         ~rFFTW_interpolator();
 
         /* map real locations to grid coordinates */
@@ -83,7 +84,6 @@ class rFFTW_interpolator
                 const double *__restrict__ xx,
                 double *__restrict__ dest,
                 const int *deriv = NULL);
-        int read_rFFTW(void *src);
 };
 
 #endif//RFFTW_INTERPOLATOR
diff --git a/bfps/fluid_base.py b/bfps/fluid_base.py
index ba44f306bce344d73ef39c51618e51b48bf98a6c..ba0be9ed8e3660dc66c9fa0a89a2eb2b8979a99f 100644
--- a/bfps/fluid_base.py
+++ b/bfps/fluid_base.py
@@ -130,9 +130,10 @@ class fluid_particle_base(bfps.code):
                              'hsize_t dset;\n')
         for key in ['state', 'velocity', 'acceleration']:
             self.definitions += ('full_name = (std::string(name) + std::string("/{0}"));\n'.format(key) +
+                                 'if (H5Lexists(g_id, full_name.c_str(), H5P_DEFAULT))\n{\n' +
                                  'dset = H5Dopen(g_id, full_name.c_str(), H5P_DEFAULT);\n' +
                                  'grow_single_dataset(dset, niter_todo/niter_part);\n' +
-                                 'H5Dclose(dset);\n')
+                                 'H5Dclose(dset);\n}\n')
         self.definitions += ('full_name = (std::string(name) + std::string("/rhs"));\n' +
                              'if (H5Lexists(g_id, full_name.c_str(), H5P_DEFAULT))\n{\n' +
                              'dset = H5Dopen(g_id, full_name.c_str(), H5P_DEFAULT);\n' +
@@ -445,25 +446,24 @@ class fluid_particle_base(bfps.code):
                                      dtype = np.int64,
                                      compression = 'gzip')
             for s in range(self.particle_species):
-                if self.parameters['tracers{0}_integration_method'.format(s)] == 'AdamsBashforth':
-                    time_chunk = 2**20 // (8*3*
-                                           self.parameters['nparticles']*
-                                           self.parameters['tracers{0}_integration_steps'.format(s)])
-                    time_chunk = max(time_chunk, 1)
-                    ofile.create_dataset('particles/tracers{0}/rhs'.format(s),
-                                         (1,
-                                          self.parameters['tracers{0}_integration_steps'.format(s)],
-                                          self.parameters['nparticles'],
-                                          3),
-                                         maxshape = (None,
-                                                     self.parameters['tracers{0}_integration_steps'.format(s)],
-                                                     self.parameters['nparticles'],
-                                                     3),
-                                         chunks =  (time_chunk,
-                                                    self.parameters['tracers{0}_integration_steps'.format(s)],
-                                                    self.parameters['nparticles'],
-                                                    3),
-                                         dtype = np.float64)
+                time_chunk = 2**20 // (8*3*
+                                       self.parameters['nparticles']*
+                                       self.parameters['tracers{0}_integration_steps'.format(s)])
+                time_chunk = max(time_chunk, 1)
+                ofile.create_dataset('particles/tracers{0}/rhs'.format(s),
+                                     (1,
+                                      self.parameters['tracers{0}_integration_steps'.format(s)],
+                                      self.parameters['nparticles'],
+                                      3),
+                                     maxshape = (None,
+                                                 self.parameters['tracers{0}_integration_steps'.format(s)],
+                                                 self.parameters['nparticles'],
+                                                 3),
+                                     chunks =  (time_chunk,
+                                                self.parameters['tracers{0}_integration_steps'.format(s)],
+                                                self.parameters['nparticles'],
+                                                3),
+                                     dtype = np.float64)
                 time_chunk = 2**20 // (8*3*self.parameters['nparticles'])
                 time_chunk = max(time_chunk, 1)
                 ofile.create_dataset(
@@ -474,14 +474,15 @@ class fluid_particle_base(bfps.code):
                     chunks = (time_chunk, self.parameters['nparticles'], 3),
                     maxshape = (None, self.parameters['nparticles'], 3),
                     dtype = np.float64)
-                ofile.create_dataset(
-                    '/particles/tracers{0}/acceleration'.format(s),
-                    (1,
-                     self.parameters['nparticles'],
-                     3),
-                    chunks = (time_chunk, self.parameters['nparticles'], 3),
-                    maxshape = (None, self.parameters['nparticles'], 3),
-                    dtype = np.float64)
+                if self.parameters['tracers{0}_acc_on'.format(s)]:
+                    ofile.create_dataset(
+                        '/particles/tracers{0}/acceleration'.format(s),
+                        (1,
+                         self.parameters['nparticles'],
+                         3),
+                        chunks = (time_chunk, self.parameters['nparticles'], 3),
+                        maxshape = (None, self.parameters['nparticles'], 3),
+                        dtype = np.float64)
             ofile.close()
         return None
 
diff --git a/done.txt b/done.txt
index e150748784fced7730cfe899f0e01a50d8fe778c..79d3a56f93a3ca8227643655bae6e46bfa4b2ee0 100644
--- a/done.txt
+++ b/done.txt
@@ -6,3 +6,7 @@ x 2015-12-26 call interpolation only when needed                     @optimizati
 x 2015-12-26 clean up tox files, make sure all tests run             @tests        +v1.0
 x 2016-01-03 check divfree function
 x 2016-01-03 compute kMeta(t) as well
+x 2016-01-03 split library into core and extra                                       @optimization +v1.0
+x 2016-01-07 FFTW interpolator doesn't need its own field                            @optimization +v1.0 +particle_api
+x 2016-01-08 simplify tracer/field addition mechanism                                @design +v1.0 +particle_api
+x 2016-01-08 add stat choice parameter to add_particles                              @design +v1.0 +particle_api
diff --git a/tests/base.py b/tests/base.py
index c050245b85910bd35e2c62bf937945ba8d897e52..abe5cf861ecc35b9211b76bd32a7f43e84a82148 100644
--- a/tests/base.py
+++ b/tests/base.py
@@ -102,29 +102,19 @@ def launch(
     c.parameters['famplitude'] = 0.2
     c.fill_up_fluid_code()
     if c.parameters['nparticles'] > 0:
-        c.add_particle_fields(
-                name = 'regular',
+        c.add_3D_rFFTW_field(name = 'rFFTW_acc')
+        c.add_interpolator(
+                name = 'spline',
                 neighbours = opt.neighbours,
                 smoothness = opt.smoothness)
-        c.add_particle_fields(kcut = 'fs->kM/2', name = 'filtered', neighbours = opt.neighbours)
         c.add_particles(
-                kcut = 'fs->kM/2',
-                integration_steps = 1,
-                fields_name = 'filtered')
-        #for integr_steps in range(1, 7):
-        #    c.add_particles(
-        #            integration_steps = integr_steps,
-        #            neighbours = opt.neighbours,
-        #            smoothness = opt.smoothness,
-        #            fields_name = 'regular')
-        for info in [(2, 'AdamsBashforth'),
-                     (3, 'AdamsBashforth'),
-                     (4, 'AdamsBashforth'),
-                     (6, 'AdamsBashforth')]:
-            c.add_particles(
-                    integration_steps = info[0],
-                    integration_method = info[1],
-                    fields_name = 'regular')
+                kcut = ['fs->kM/2', 'fs->kM/3'],
+                integration_steps = 3,
+                interpolator = 'spline')
+        c.add_particles(
+                integration_steps = [2, 3, 4, 6],
+                interpolator = 'spline',
+                acc_name = 'rFFTW_acc')
     c.finalize_code()
     c.write_src()
     c.write_par()
@@ -158,6 +148,8 @@ def launch(
 
 
 def acceleration_test(c, m = 3, species = 0):
+    if not c.parameters['tracers{0}_acc_on'.format(species)]:
+        return None
     import numpy as np
     import matplotlib.pyplot as plt
     from bfps.tools import get_fornberg_coeffs
@@ -182,13 +174,12 @@ def acceleration_test(c, m = 3, species = 0):
     pid = np.argmin(SNR(num_acc1, acc[n+1:-n-1]))
     pars = d['parameters']
     to_print = (
-            'integration={0}, steps={1}, interp={2}, neighbours={3}, '.format(
-                pars['tracers{0}_integration_method'.format(species)].value,
+            'steps={0}, interp={1}, neighbours={2}, '.format(
                 pars['tracers{0}_integration_steps'.format(species)].value,
-                pars[str(pars['tracers{0}_field'.format(species)].value) + '_type'].value,
-                pars[str(pars['tracers{0}_field'.format(species)].value) + '_neighbours'].value))
-    if 'spline' in pars['tracers{0}_field'.format(species)].value:
-        to_print += 'smoothness = {0}, '.format(pars[str(pars['tracers{0}_field'.format(species)].value) + '_smoothness'].value)
+                pars[str(pars['tracers{0}_interpolator'.format(species)].value) + '_type'].value,
+                pars[str(pars['tracers{0}_interpolator'.format(species)].value) + '_neighbours'].value))
+    if 'spline' in pars['tracers{0}_interpolator'.format(species)].value:
+        to_print += 'smoothness = {0}, '.format(pars[str(pars['tracers{0}_interpolator'.format(species)].value) + '_smoothness'].value)
     to_print += (
             'SNR d1p-vel={0:.3f}, d1v-acc={1:.3f}, d2p-acc={2:.3f}'.format(
                 np.mean(SNR(num_vel1, vel[n+1:-n-1])),
diff --git a/todo.txt b/todo.txt
index a4dac0865e9dc62c2ce01076eadc69499eb3b2c2..41ab0f009a43dd95179363bdb499aab43624983c 100644
--- a/todo.txt
+++ b/todo.txt
@@ -1,19 +1,18 @@
-(B) FFTW interpolator doesn't need its own field            @optimization +v1.0
-(B) compute z polynomials only when needed                  @optimization +v1.0
+(B) compute z polynomials only when needed                                  @optimization +v1.0
 (B) read https://www.xsede.org/documents/271087/369161/ExtScale-Koziol.pdf  @optimization @HDF5 +I/O
-(B) set up mechanism for adding in new PDEs                 @design +v2.0 +alternate_algorithms
-(B) split library into core and extra                       @optimization +v1.0
-(B) tweak HDF5 settings                                     @optimization @HDF5 +I/O
-(B) use less memory                                         @optimization
-(C) clean up machine_settings mess                          @design @documentation +v2.0
-(C) code overview                                           @documentation
-(C) move stat I/O to cpp lib                                @design @HDF5
-(C) test involving hydrodynamic similarity                  @tests
-(C) use HDF5 io for fields                                  @design @HDF5 +I/O
-(D) generalize interpolation comparison test                @tests
-(D) test anisotropic grids                                  @tests
-(D) test non-cubic domains                                  @tests
-(D) tests should not overwrite other tests (tox_full)       @tests
-(E) add u-equation algorithm for testing purposes           @tests +alternate_algorithms
-(E) pure python DNS addon: pros and cons                    @tests +alternate_algorithms
+(B) set up mechanism for adding in new PDEs                                 @design +v2.0 +alternate_algorithms
+(B) tweak HDF5 settings                                                     @optimization @HDF5 +I/O
+(B) use less memory                                                         @optimization
+(C) clean up machine_settings mess                                          @design @documentation +v2.0
+(C) code overview                                                           @documentation
+(C) move stat I/O to cpp lib                                                @design @HDF5
+(C) test involving hydrodynamic similarity                                  @tests
+(C) use HDF5 io for fields                                                  @design @HDF5 +I/O
+(D) generalize interpolation comparison test                                @tests
+(D) generate separate lib(s) with extra classes                             @tests +alternate_algorithms
+(D) test anisotropic grids                                                  @tests
+(D) test non-cubic domains                                                  @tests
+(D) tests should not overwrite other tests (tox_full)                       @tests
+(E) add u-equation algorithm for testing purposes                           @tests +alternate_algorithms
+(E) pure python DNS addon: pros and cons                                    @tests +alternate_algorithms
 (F) add switch to turn off simulation