diff --git a/.gitignore b/.gitignore
index 3717951f2c8409fbe95675bc9975a689154f8c77..290ce8e46600e7d9e9b6e497fb0709d1c47070d1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -6,3 +6,7 @@ meta/.ipynb_checkpoints/*
 dist
 build
 bfps.egg-info
+MANIFEST.in
+bfps/install_info.pickle
+bfps/lib*.a
+obj
diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index 20e1093181a2cd018cb4a1f53a7a89979d93ade2..f3fa39d04c76ac5b400b481ecf423095b4135bcd 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -408,8 +408,7 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
             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}->temp);\n' +
-                          'acc_{0}->read_rFFTW(acc_{0}->temp);\n').format(name)
+                          'fs->compute_Lagrangian_acceleration(acc_{0}->field);\n').format(name)
         self.fluid_start += update_fields
         self.fluid_loop += update_fields
         return None
@@ -450,7 +449,7 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
                               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}->fs->rd->myrank == 0)
+                              if (ps{0}->myrank == 0)
                               {{
                                   //VELOCITY begin
                                   std::string temp_string = (std::string("/particles/") +
@@ -512,7 +511,7 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
                     neighbours,
                     self.particle_species)
             self.particle_start += ('ps{0} = new rFFTW_particles<VELOCITY_TRACER, {1}, {2}>(\n' +
-                                    'fname, fs, vel_{3},\n' +
+                                    'fname, vel_{3},\n' +
                                     'nparticles,\n' +
                                     'niter_part, tracers{0}_integration_steps);\n').format(
                                             self.particle_species, self.C_dtype, neighbours, fields_name)
@@ -536,10 +535,11 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
                 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_loop += (('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_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
         return None
     def get_data_file(self):
diff --git a/bfps/cpp/rFFTW_interpolator.cpp b/bfps/cpp/rFFTW_interpolator.cpp
index 6f1a49bc7b538b94286c397ab4613540a1b72f9f..a858eb104abf19837f7ccd30ae673caaf841a493 100644
--- a/bfps/cpp/rFFTW_interpolator.cpp
+++ b/bfps/cpp/rFFTW_interpolator.cpp
@@ -26,6 +26,7 @@
 
 #define NDEBUG
 
+#include <cmath>
 #include "rFFTW_interpolator.hpp"
 
 template <class rnumber, int interp_neighbours>
@@ -37,58 +38,104 @@ rFFTW_interpolator<rnumber, interp_neighbours>::rFFTW_interpolator(
     this->field_size = 2*fs->cd->local_size;
     this->compute_beta = BETA_POLYS;
     if (sizeof(rnumber) == 4)
-    {
-        this->f0 = (rnumber*)((void*)fftwf_alloc_real(this->field_size));
-        this->f1 = (rnumber*)((void*)fftwf_alloc_real(this->field_size));
-    }
+        this->field = (rnumber*)((void*)fftwf_alloc_real(this->field_size));
     else if (sizeof(rnumber) == 8)
-    {
-        this->f0 = (rnumber*)((void*)fftw_alloc_real(this->field_size));
-        this->f1 = (rnumber*)((void*)fftw_alloc_real(this->field_size));
-    }
-    this->temp = this->f1;
+        this->field = (rnumber*)((void*)fftw_alloc_real(this->field_size));
+
+    // compute dx, dy, dz;
+    this->dx = 4*acos(0) / (fs->dkx*this->descriptor->sizes[2]);
+    this->dy = 4*acos(0) / (fs->dky*this->descriptor->sizes[1]);
+    this->dz = 4*acos(0) / (fs->dkz*this->descriptor->sizes[0]);
+
+    // generate compute array
+    this->compute = new bool[this->descriptor->sizes[0]];
+    std::fill_n(this->compute, this->descriptor->sizes[0], false);
+    for (int iz = this->descriptor->starts[0]-interp_neighbours-1;
+            iz <= this->descriptor->starts[0]+this->descriptor->subsizes[0]+interp_neighbours;
+            iz++)
+        this->compute[((iz + this->descriptor->sizes[0]) % this->descriptor->sizes[0])] = true;
 }
 
 template <class rnumber, int interp_neighbours>
 rFFTW_interpolator<rnumber, interp_neighbours>::~rFFTW_interpolator()
 {
     if (sizeof(rnumber) == 4)
-    {
-        fftwf_free((float*)((void*)this->f0));
-        fftwf_free((float*)((void*)this->f1));
-    }
+        fftwf_free((float*)((void*)this->field));
     else if (sizeof(rnumber) == 8)
-    {
-        fftw_free((double*)((void*)this->f0));
-        fftw_free((double*)((void*)this->f1));
-    }
+        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)
 {
-    /* first, roll fields */
-    rnumber *tmp = this->f0;
-    this->f0 = this->f1;
-    this->f1 = tmp;
-    this->temp = this->f0;
-    /* now do regular things */
     rnumber *src = (rnumber*)void_src;
-    rnumber *dst = this->f1;
     /* do big copy of middle stuff */
     std::copy(src,
               src + this->field_size,
-              dst);
+              this->field);
     return EXIT_SUCCESS;
 }
 
 template <class rnumber, int interp_neighbours>
-void rFFTW_interpolator<rnumber, interp_neighbours>::operator()(
-        double t,
+void rFFTW_interpolator<rnumber, interp_neighbours>::get_grid_coordinates(
+        const int nparticles,
+        const int pdimension,
+        const double *x,
         int *xg,
-        double *xx,
+        double *xx)
+{
+    static double grid_size[] = {this->dx, this->dy, this->dz};
+    double tval;
+    std::fill_n(xg, nparticles*3, 0);
+    std::fill_n(xx, nparticles*3, 0.0);
+    for (int p=0; p<nparticles; p++)
+    {
+        for (int c=0; c<3; c++)
+        {
+            tval = floor(x[p*pdimension+c]/grid_size[c]);
+            xg[p*3+c] = MOD(int(tval), this->descriptor->sizes[2-c]);
+            xx[p*3+c] = (x[p*pdimension+c] - tval*grid_size[c]) / grid_size[c];
+        }
+    }
+}
+
+template <class rnumber, int interp_neighbours>
+void rFFTW_interpolator<rnumber, interp_neighbours>::sample(
+        const int nparticles,
+        const int pdimension,
+        const double *__restrict__ x,
+        double *__restrict__ y,
+        const int *deriv)
+{
+    /* get grid coordinates */
+    int *xg = new int[3*nparticles];
+    double *xx = new double[3*nparticles];
+    double *yy =  new double[3*nparticles];
+    std::fill_n(yy, 3*nparticles, 0.0);
+    this->get_grid_coordinates(nparticles, pdimension, x, xg, xx);
+    /* perform interpolation */
+    for (int p=0; p<nparticles; p++)
+        if (this->compute[xg[p*3+2]])
+            this->operator()(xg + p*3, xx + p*3, yy + p*3, deriv);
+    MPI_Allreduce(
+            yy,
+            y,
+            3*nparticles,
+            MPI_DOUBLE,
+            MPI_SUM,
+            this->descriptor->comm);
+    delete[] yy;
+    delete[] xg;
+    delete[] xx;
+}
+
+template <class rnumber, int interp_neighbours>
+void rFFTW_interpolator<rnumber, interp_neighbours>::operator()(
+        const int *xg,
+        const double *xx,
         double *dest,
-        int *deriv)
+        const int *deriv)
 {
     double bx[interp_neighbours*2+2], by[interp_neighbours*2+2], bz[interp_neighbours*2+2];
     if (deriv == NULL)
@@ -105,13 +152,11 @@ void rFFTW_interpolator<rnumber, interp_neighbours>::operator()(
     }
     std::fill_n(dest, 3, 0);
     ptrdiff_t bigiz, bigiy, bigix;
-    double tval[3];
     for (int iz = -interp_neighbours; iz <= interp_neighbours+1; iz++)
     {
         bigiz = ptrdiff_t(((xg[2]+iz) + this->descriptor->sizes[0]) % this->descriptor->sizes[0]);
         if (this->descriptor->myrank == this->descriptor->rank[bigiz])
         {
-            std::fill_n(tval, 3, 0);
             for (int iy = -interp_neighbours; iy <= interp_neighbours+1; iy++)
             {
                 bigiy = ptrdiff_t(MOD(xg[1]+iy, this->descriptor->sizes[1]));
@@ -122,17 +167,11 @@ void rFFTW_interpolator<rnumber, interp_neighbours>::operator()(
                                          bigiy)*(this->descriptor->sizes[2]+2) +
                                          bigix)*3;
                     for (int c=0; c<3; c++)
-                    {
-                        dest[c] += (this->f0[tindex+c]*(1-t) + t*this->f1[tindex+c])*(bz[iz+interp_neighbours]*
-                                                                                      by[iy+interp_neighbours]*
-                                                                                      bx[ix+interp_neighbours]);
-                        tval[c] += (this->f0[tindex+c]*(1-t) + t*this->f1[tindex+c])*(bz[iz+interp_neighbours]*
-                                                                                      by[iy+interp_neighbours]*
-                                                                                      bx[ix+interp_neighbours]);
-                    }
+                        dest[c] += this->field[tindex+c]*(bz[iz+interp_neighbours]*
+                                                          by[iy+interp_neighbours]*
+                                                          bx[ix+interp_neighbours]);
                 }
             }
-            DEBUG_MSG("%ld %d %d %g %g %g\n", bigiz, xg[1], xg[0], tval[0], tval[1], tval[2]);
         }
     }
 }
diff --git a/bfps/cpp/rFFTW_interpolator.hpp b/bfps/cpp/rFFTW_interpolator.hpp
index 7a32e75808ce572b2e902d2ebd0a22b489a4b81e..0864bd964b77b5b3527a58898b8c22da63c33c2e 100644
--- a/bfps/cpp/rFFTW_interpolator.hpp
+++ b/bfps/cpp/rFFTW_interpolator.hpp
@@ -44,17 +44,52 @@ template <class rnumber, int interp_neighbours>
 class rFFTW_interpolator
 {
     public:
+        /* size of field to interpolate */
         ptrdiff_t field_size;
+
+        /* pointer to polynomial function */
         base_polynomial_values compute_beta;
+
+        /* descriptor of field to interpolate */
         field_descriptor<rnumber> *descriptor;
-        rnumber *f0, *f1, *temp;
+
+        /* pointers to fields that are to be interpolated
+         * */
+        rnumber *field;
+
+        /* physical parameters of field */
+        double dx, dy, dz;
+
+        /* compute[iz] is true if .
+         * local_zstart - neighbours <= iz <= local_zend + 1 + neighbours
+         * */
+        bool *compute;
 
         rFFTW_interpolator(
                 fluid_solver_base<rnumber> *FSOLVER,
                 base_polynomial_values BETA_POLYS);
         ~rFFTW_interpolator();
 
-        void operator()(double t, int *__restrict__ xg, double *__restrict__ xx, double *__restrict__ dest, int *deriv = NULL);
+        /* map real locations to grid coordinates */
+        void get_grid_coordinates(
+                const int nparticles,
+                const int pdimension,
+                const double *__restrict__ x,
+                int *__restrict__ xg,
+                double *__restrict__ xx);
+        /* interpolate field at an array of locations */
+        void sample(
+                const int nparticles,
+                const int pdimension,
+                const double *__restrict__ x,
+                double *__restrict__ y,
+                const int *deriv = NULL);
+        /* interpolate 1 point */
+        void operator()(
+                const int *__restrict__ xg,
+                const double *__restrict__ xx,
+                double *__restrict__ dest,
+                const int *deriv = NULL);
         int read_rFFTW(void *src);
 };
 
diff --git a/bfps/cpp/rFFTW_particles.cpp b/bfps/cpp/rFFTW_particles.cpp
index af225f9541d62adf93f1ff89aeadfb9a7cbe0cf1..4cd79c1259940861c1b0a8f07d147bc5a1356250 100644
--- a/bfps/cpp/rFFTW_particles.cpp
+++ b/bfps/cpp/rFFTW_particles.cpp
@@ -24,7 +24,7 @@
 
 
 
-//#define NDEBUG
+#define NDEBUG
 
 #include <cmath>
 #include <cassert>
@@ -42,7 +42,6 @@ extern int myrank, nprocs;
 template <int particle_type, class rnumber, int interp_neighbours>
 rFFTW_particles<particle_type, rnumber, interp_neighbours>::rFFTW_particles(
         const char *NAME,
-        fluid_solver_base<rnumber> *FSOLVER,
         rFFTW_interpolator<rnumber, interp_neighbours> *FIELD,
         const int NPARTICLES,
         const int TRAJ_SKIP,
@@ -60,7 +59,6 @@ rFFTW_particles<particle_type, rnumber, interp_neighbours>::rFFTW_particles(
     assert((INTEGRATION_STEPS <= 6) &&
            (INTEGRATION_STEPS >= 1));
     strncpy(this->name, NAME, 256);
-    this->fs = FSOLVER;
     this->nparticles = NPARTICLES;
     this->vel = FIELD;
     this->integration_steps = INTEGRATION_STEPS;
@@ -76,41 +74,6 @@ rFFTW_particles<particle_type, rnumber, interp_neighbours>::rFFTW_particles(
         this->rhs[i] = new double[this->array_size];
         std::fill_n(this->rhs[i], this->array_size, 0.0);
     }
-
-    // compute dx, dy, dz;
-    this->dx = 4*acos(0) / (this->fs->dkx*this->fs->rd->sizes[2]);
-    this->dy = 4*acos(0) / (this->fs->dky*this->fs->rd->sizes[1]);
-    this->dz = 4*acos(0) / (this->fs->dkz*this->fs->rd->sizes[0]);
-
-    // compute lower and upper bounds
-    this->lbound = new double[nprocs];
-    this->ubound = new double[nprocs];
-    double *tbound = new double[nprocs];
-    std::fill_n(tbound, nprocs, 0.0);
-    tbound[this->myrank] = this->fs->rd->starts[0]*this->dz;
-    MPI_Allreduce(
-            tbound,
-            this->lbound,
-            nprocs,
-            MPI_DOUBLE,
-            MPI_SUM,
-            this->comm);
-    std::fill_n(tbound, nprocs, 0.0);
-    tbound[this->myrank] = (this->fs->rd->starts[0] + this->fs->rd->subsizes[0])*this->dz;
-    MPI_Allreduce(
-            tbound,
-            this->ubound,
-            nprocs,
-            MPI_DOUBLE,
-            MPI_SUM,
-            this->comm);
-    delete[] tbound;
-    //for (int r = 0; r<nprocs; r++)
-    //    DEBUG_MSG(
-    //            "lbound[%d] = %lg, ubound[%d] = %lg\n",
-    //            r, this->lbound[r],
-    //            r, this->ubound[r]
-    //            );
 }
 
 template <int particle_type, class rnumber, int interp_neighbours>
@@ -121,59 +84,19 @@ rFFTW_particles<particle_type, rnumber, interp_neighbours>::~rFFTW_particles()
     {
         delete[] this->rhs[i];
     }
-    delete[] this->lbound;
-    delete[] this->ubound;
-}
-
-template <int particle_type, class rnumber,  int interp_neighbours>
-void rFFTW_particles<particle_type, rnumber, interp_neighbours>::sample_vec_field(
-        rFFTW_interpolator<rnumber, interp_neighbours> *vec,
-        double t,
-        double *x,
-        double *y,
-        int *deriv)
-{
-    /* get grid coordinates */
-    int *xg = new int[3*this->nparticles];
-    double *xx = new double[3*this->nparticles];
-    double *yy =  new double[3*this->nparticles];
-    std::fill_n(yy, 3*this->nparticles, 0.0);
-    this->get_grid_coordinates(x, xg, xx);
-    /* perform interpolation */
-    for (int p=0; p<this->nparticles; p++)
-        (*vec)(t, xg + p*3, xx + p*3, yy + p*3, deriv);
-    MPI_Allreduce(
-            yy,
-            y,
-            3*this->nparticles,
-            MPI_DOUBLE,
-            MPI_SUM,
-            this->comm);
-    delete[] yy;
-    delete[] xg;
-    delete[] xx;
 }
 
 template <int particle_type, class rnumber, int interp_neighbours>
-void rFFTW_particles<particle_type, rnumber, interp_neighbours>::get_rhs(double t, double *x, double *y)
+void rFFTW_particles<particle_type, rnumber, interp_neighbours>::get_rhs(double *x, double *y)
 {
     switch(particle_type)
     {
         case VELOCITY_TRACER:
-            this->sample_vec_field(this->vel, t, x, y);
+            this->vel->sample(this->nparticles, this->ncomponents, x, y);
             break;
     }
 }
 
-
-template <int particle_type, class rnumber, int interp_neighbours>
-int rFFTW_particles<particle_type, rnumber, interp_neighbours>::get_rank(double z)
-{
-    int tmp = this->fs->rd->rank[MOD(int(floor(z/this->dz)), this->fs->rd->sizes[0])];
-    assert(tmp >= 0 && tmp < this->fs->rd->nprocs);
-    return tmp;
-}
-
 template <int particle_type, class rnumber, int interp_neighbours>
 void rFFTW_particles<particle_type, rnumber, interp_neighbours>::roll_rhs()
 {
@@ -189,7 +112,7 @@ template <int particle_type, class rnumber, int interp_neighbours>
 void rFFTW_particles<particle_type, rnumber, interp_neighbours>::AdamsBashforth(int nsteps)
 {
     ptrdiff_t ii;
-    this->get_rhs(0, this->state, this->rhs[0]);
+    this->get_rhs(this->state, this->rhs[0]);
     switch(nsteps)
     {
         case 1:
@@ -270,29 +193,6 @@ void rFFTW_particles<particle_type, rnumber, interp_neighbours>::step()
 }
 
 
-
-template <int particle_type, class rnumber, int interp_neighbours>
-void rFFTW_particles<particle_type, rnumber, interp_neighbours>::get_grid_coordinates(double *x, int *xg, double *xx)
-{
-    static double grid_size[] = {this->dx, this->dy, this->dz};
-    double tval;
-    std::fill_n(xg, this->nparticles*3, 0);
-    std::fill_n(xx, this->nparticles*3, 0.0);
-    for (int p=0; p<this->nparticles; p++)
-    {
-        for (int c=0; c<3; c++)
-        {
-            tval = floor(x[p*this->ncomponents+c]/grid_size[c]);
-            xg[p*3+c] = MOD(int(tval), this->fs->rd->sizes[2-c]);
-            xx[p*3+c] = (x[p*this->ncomponents+c] - tval*grid_size[c]) / grid_size[c];
-        }
-        /*xg[p*3+2] -= this->fs->rd->starts[0];
-        if (this->myrank == this->fs->rd->rank[0] &&
-            xg[p*3+2] > this->fs->rd->subsizes[0])
-            xg[p*3+2] -= this->fs->rd->sizes[0];*/
-    }
-}
-
 template <int particle_type, class rnumber, int interp_neighbours>
 void rFFTW_particles<particle_type, rnumber, interp_neighbours>::read(hid_t data_file_id)
 {
diff --git a/bfps/cpp/rFFTW_particles.hpp b/bfps/cpp/rFFTW_particles.hpp
index ac44cc784d4a1f08bd50e7e3aa6070d3e399cbbe..68f8badea0359436d3aaf3f06b02543ee9ff0d89 100644
--- a/bfps/cpp/rFFTW_particles.hpp
+++ b/bfps/cpp/rFFTW_particles.hpp
@@ -43,20 +43,8 @@ class rFFTW_particles
     public:
         int myrank, nprocs;
         MPI_Comm comm;
-        fluid_solver_base<rnumber> *fs;
         rnumber *data;
 
-        /* watching is an array of shape [nparticles], with
-         * watching[p] being true if particle p is in the domain of myrank
-         * or in the buffer regions.
-         * only used if multistep is false.
-         * */
-        bool *watching;
-        /* computing is an array of shape [nparticles], with
-         * computing[p] being the rank that is currently working on particle p
-         * */
-        int *computing;
-
         /* state will generally hold all the information about the particles.
          * in the beginning, we will only need to solve 3D ODEs, but I figured
          * a general ncomponents is better, since we may change our minds.
@@ -68,8 +56,6 @@ class rFFTW_particles
         int array_size;
         int integration_steps;
         int traj_skip;
-        double *lbound;
-        double *ubound;
         rFFTW_interpolator<rnumber, interp_neighbours> *vel;
 
         /* simulation parameters */
@@ -77,23 +63,15 @@ class rFFTW_particles
         int iteration;
         double dt;
 
-        /* physical parameters of field */
-        double dx, dy, dz;
-
         /* methods */
 
         /* constructor and destructor.
          * allocate and deallocate:
          *  this->state
          *  this->rhs
-         *  this->lbound
-         *  this->ubound
-         *  this->computing
-         *  this->watching
          * */
         rFFTW_particles(
                 const char *NAME,
-                fluid_solver_base<rnumber> *FSOLVER,
                 rFFTW_interpolator<rnumber, interp_neighbours> *FIELD,
                 const int NPARTICLES,
                 const int TRAJ_SKIP,
@@ -101,19 +79,10 @@ class rFFTW_particles
         ~rFFTW_particles();
 
         void get_rhs(double *__restrict__ x, double *__restrict__ rhs);
-        void get_rhs(double t, double *__restrict__ x, double *__restrict__ rhs);
-
-        int get_rank(double z); // get rank for given value of z
-        void get_grid_coordinates(double *__restrict__ x, int *__restrict__ xg, double *__restrict__ xx);
-        void sample_vec_field(
-            rFFTW_interpolator<rnumber, interp_neighbours> *vec,
-            double t,
-            double *__restrict__ x,
-            double *__restrict__ y,
-            int *deriv = NULL);
+
         inline void sample_vec_field(rFFTW_interpolator<rnumber, interp_neighbours> *field, double *vec_values)
         {
-            this->sample_vec_field(field, 1.0, this->state, vec_values, NULL);
+            field->sample(this->nparticles, this->ncomponents, this->state, vec_values, NULL);
         }
 
         /* input/output */
diff --git a/bfps/fluid_base.py b/bfps/fluid_base.py
index 62dfcfd264e2b6f0987ac0cfe547873e840fb686..ba44f306bce344d73ef39c51618e51b48bf98a6c 100644
--- a/bfps/fluid_base.py
+++ b/bfps/fluid_base.py
@@ -90,13 +90,14 @@ class fluid_particle_base(bfps.code):
         self.fluid_loop = ''
         self.fluid_end  = ''
         self.fluid_output = ''
+        self.stat_src = ''
         self.particle_includes = ''
         self.particle_variables = ''
         self.particle_definitions = ''
         self.particle_start = ''
         self.particle_loop = ''
         self.particle_end  = ''
-        self.stat_src = ''
+        self.particle_stat_src = ''
         self.file_datasets_grow   = ''
         return None
     def finalize_code(self):
@@ -144,6 +145,7 @@ class fluid_particle_base(bfps.code):
                              'return file_problems;\n'
                              '}\n')
         self.definitions += 'void do_stats()\n{\n' + self.stat_src + '}\n'
+        self.definitions += 'void do_particle_stats()\n{\n' + self.particle_stat_src + '}\n'
         # take care of wisdom
         if self.use_fftw_wisdom:
             if self.dtype == np.float32:
@@ -191,6 +193,7 @@ class fluid_particle_base(bfps.code):
                                return EXIT_SUCCESS;
                            }
                            do_stats();
+                           do_particle_stats();
                            //endcpp
                            """
         output_time_difference = ('time1 = clock();\n' +
@@ -202,16 +205,21 @@ class fluid_particle_base(bfps.code):
                                       '<< iteration << " took " ' +
                                       '<< time_difference/nprocs << " seconds" << std::endl;\n' +
                                   'time0 = time1;\n')
-        self.main       += 'for (int max_iter = iteration+niter_todo; iteration < max_iter; iteration++)\n{\n'
+        self.main       += 'for (int max_iter = iteration+niter_todo; iteration < max_iter; iteration++)\n'
+        self.main       += '{\n'
         self.main       += output_time_difference
-        self.main       += self.fluid_loop
         if self.particle_species > 0:
             self.main   += self.particle_loop
-        self.main       += 'if (iteration % niter_stat == 0) do_stats();\n}\n'
+        self.main       += self.fluid_loop
+        self.main       += 'if (iteration % niter_stat == 0) do_stats();\n'
+        if self.particle_species > 0:
+            self.main       += 'if (iteration % niter_part == 0) do_particle_stats();\n'
+        self.main       += '}\n'
         self.main       += output_time_difference
+        self.main       += 'do_stats();\n'
+        self.main       += 'do_particle_stats();\n'
         if self.particle_species > 0:
             self.main   += self.particle_end
-        self.main       += 'do_stats();\n'
         self.main       += self.fluid_end
         return None
     def read_rfield(
diff --git a/done.txt b/done.txt
index 787392010cce3f4cb0f5921eba4a7fa8ba31e89a..344832ae21046d40f47702000ab8b31d0128114b 100644
--- a/done.txt
+++ b/done.txt
@@ -1,2 +1,5 @@
 x 2015-12-04 make code py3 compatible                                @python3
 x 2015-12-23 decide on versioning system                                           +merge0
+x 2015-12-24 move get_grid coords to interpolator                    @optimization +v1.0
+x 2015-12-25 get rid of temporal interpolation                       @optimization +v1.0
+x 2015-12-26 call interpolation only when needed                     @optimization +v1.0
diff --git a/tests/test_base.py b/tests/base.py
similarity index 69%
rename from tests/test_base.py
rename to tests/base.py
index 3be7d2f5dfdbcca79f0615fde141f340e8bd5be8..c050245b85910bd35e2c62bf937945ba8d897e52 100755
--- a/tests/test_base.py
+++ b/tests/base.py
@@ -25,9 +25,13 @@
 
 
 import os
+import sys
 import subprocess
 import pickle
 
+import numpy as np
+import matplotlib.pyplot as plt
+
 import bfps
 from bfps import fluid_resize
 
@@ -211,6 +215,93 @@ def acceleration_test(c, m = 3, species = 0):
     plt.close(fig)
     return pid
 
+def compare_stats(
+        opt,
+        c0, c1,
+        plots_on = False):
+    for key in ['energy', 'enstrophy', 'vel_max']:
+        print('maximum {0} difference is {1}'.format(
+            key,
+            np.max(np.abs(c0.statistics[key + '(t)'] - c0.statistics[key + '(t)']))))
+    for i in range(c0.particle_species):
+        print('maximum traj difference species {0} is {1}'.format(
+            i,
+            np.max(np.abs(c0.trajectories[i] - c1.trajectories[i]))))
+    if plots_on:
+        # plot energy and enstrophy
+        fig = plt.figure(figsize = (12, 12))
+        a = fig.add_subplot(221)
+        c0.set_plt_style({'label' : '1',
+                          'dashes' : (None, None),
+                          'color' : (1, 0, 0)})
+        c1.set_plt_style({'label' : '2',
+                          'dashes' : (2, 2),
+                          'color' : (0, 0, 1)})
+        for c in [c0, c1]:
+            a.plot(c.statistics['t'],
+                   c.statistics['energy(t)'],
+                   label = c.style['label'],
+                   dashes = c.style['dashes'],
+                   color = c.style['color'])
+        a.set_title('energy')
+        a.legend(loc = 'best')
+        a = fig.add_subplot(222)
+        for c in [c0, c1]:
+            a.plot(c.statistics['t'],
+                   c.statistics['enstrophy(t)'],
+                   dashes = c.style['dashes'],
+                   color = c.style['color'])
+        a.set_title('enstrophy')
+        a = fig.add_subplot(223)
+        for c in [c0, c1]:
+            a.plot(c.statistics['t'],
+                   c.statistics['kM']*c.statistics['etaK(t)'],
+                   dashes = c.style['dashes'],
+                   color = c.style['color'])
+        a.set_title('$k_M \\eta_K$')
+        a = fig.add_subplot(224)
+        for c in [c0, c1]:
+            a.plot(c.statistics['t'],
+                   c.statistics['vel_max(t)'] * (c.parameters['dt'] * c.parameters['dkx'] /
+                                                 (2*np.pi / c.parameters['nx'])),
+                   dashes = c.style['dashes'],
+                   color = c.style['color'])
+        a.set_title('$\\frac{\\Delta t \\| u \\|_\infty}{\\Delta x}$')
+        fig.savefig('plain_stats_{0}.pdf'.format(opt.precision), format = 'pdf')
+
+        fig = plt.figure(figsize = (12, 12))
+        a = fig.add_subplot(221)
+        a.plot(c0.statistics['t'],
+               c0.statistics['energy(t)'] - c1.statistics['energy(t)'])
+        a.set_title('energy')
+        a = fig.add_subplot(222)
+        a.plot(c0.statistics['t'],
+               c0.statistics['enstrophy(t)'] - c1.statistics['enstrophy(t)'])
+        a.set_title('enstrophy')
+        a = fig.add_subplot(223)
+        a.plot(c0.statistics['t'],
+               c0.statistics['kM']*c0.statistics['etaK(t)'] - c1.statistics['kM']*c1.statistics['etaK(t)'])
+        a.set_title('$k_M \\eta_K$')
+        a = fig.add_subplot(224)
+        data0 = c0.statistics['vel_max(t)'] * (c0.parameters['dt'] * c0.parameters['dkx'] /
+                                               (2*np.pi / c0.parameters['nx']))
+        data1 = c1.statistics['vel_max(t)'] * (c1.parameters['dt'] * c1.parameters['dkx'] /
+                                               (2*np.pi / c1.parameters['nx']))
+        a.plot(c0.statistics['t'],
+               data0 - data1)
+        a.set_title('$\\frac{\\Delta t \\| u \\|_\infty}{\\Delta x}$')
+        fig.savefig('plain_stat_diffs_{0}.pdf'.format(opt.precision), format = 'pdf')
+
+        # plot trajectory differences
+        for i in range(c0.particle_species):
+            fig = plt.figure(figsize=(12, 4))
+            for j in range(3):
+                a = fig.add_subplot(131 + j)
+                for t in range(c0.parameters['nparticles']):
+                    a.plot(c0.trajectories[i][:, t, j] - c1.trajectories[i][:, t, j])
+            fig.savefig('traj_s{0}_{1}.pdf'.format(i, opt.precision), format = 'pdf')
+    return None
+
 if __name__ == '__main__':
     print('this file doesn\'t do anything')
 
diff --git a/tests/test_plain.py b/tests/test_plain.py
index 7437f1f61c9694c0e8828b4946c54e840a22461e..8ada988ad741d9a2fe177c8e2ea09369a9e69d39 100755
--- a/tests/test_plain.py
+++ b/tests/test_plain.py
@@ -25,7 +25,7 @@
 
 
 
-from test_base import *
+from base import *
 
 import numpy as np
 import matplotlib.pyplot as plt
@@ -58,84 +58,19 @@ def plain(opt):
     opt.niter_todo = 2*opt.niter_todo//3
     c2 = launch(opt, dt = c0.parameters['dt'])
     c2.compute_statistics()
-    # plot energy and enstrophy
-    fig = plt.figure(figsize = (12, 12))
-    a = fig.add_subplot(221)
-    c0.set_plt_style({'label' : '1',
-                      'dashes' : (None, None),
-                      'color' : (1, 0, 0)})
-    c1.set_plt_style({'label' : '2',
-                      'dashes' : (2, 2),
-                      'color' : (0, 0, 1)})
-    c2.set_plt_style({'label' : '3',
-                      'dashes' : (3, 3),
-                      'color' : (0, 1, 0)})
-    for c in [c0, c1, c2]:
-        a.plot(c.statistics['t'],
-               c.statistics['energy(t)'],
-               label = c.style['label'],
-               dashes = c.style['dashes'],
-               color = c.style['color'])
-    a.set_title('energy')
-    a.legend(loc = 'best')
-    a = fig.add_subplot(222)
-    for c in [c0, c1, c2]:
-        a.plot(c.statistics['t'],
-               c.statistics['enstrophy(t)'],
-               dashes = c.style['dashes'],
-               color = c.style['color'])
-    a.set_title('enstrophy')
-    a = fig.add_subplot(223)
-    for c in [c0, c1, c2]:
-        a.plot(c.statistics['t'],
-               c.statistics['kM']*c.statistics['etaK(t)'],
-               dashes = c.style['dashes'],
-               color = c.style['color'])
-    a.set_title('$k_M \\eta_K$')
-    a = fig.add_subplot(224)
-    for c in [c0, c1, c2]:
-        a.plot(c.statistics['t'],
-               c.statistics['vel_max(t)'] * (c.parameters['dt'] * c.parameters['dkx'] /
-                                             (2*np.pi / c.parameters['nx'])),
-               dashes = c.style['dashes'],
-               color = c.style['color'])
-    a.set_title('$\\frac{\\Delta t \\| u \\|_\infty}{\\Delta x}$')
-    fig.savefig('plain_stats_{0}.pdf'.format(opt.precision), format = 'pdf')
-
-    fig = plt.figure(figsize = (12, 12))
-    a = fig.add_subplot(221)
-    a.plot(c0.statistics['t'],
-           c0.statistics['energy(t)'] - c1.statistics['energy(t)'])
-    a.set_title('energy')
-    a = fig.add_subplot(222)
-    a.plot(c0.statistics['t'],
-           c0.statistics['enstrophy(t)'] - c1.statistics['enstrophy(t)'])
-    a.set_title('enstrophy')
-    a = fig.add_subplot(223)
-    a.plot(c0.statistics['t'],
-           c0.statistics['kM']*c0.statistics['etaK(t)'] - c1.statistics['kM']*c1.statistics['etaK(t)'])
-    a.set_title('$k_M \\eta_K$')
-    a = fig.add_subplot(224)
-    data0 = c0.statistics['vel_max(t)'] * (c0.parameters['dt'] * c0.parameters['dkx'] /
-                                           (2*np.pi / c0.parameters['nx']))
-    data1 = c1.statistics['vel_max(t)'] * (c1.parameters['dt'] * c1.parameters['dkx'] /
-                                           (2*np.pi / c1.parameters['nx']))
-    a.plot(c0.statistics['t'],
-           data0 - data1)
-    a.set_title('$\\frac{\\Delta t \\| u \\|_\infty}{\\Delta x}$')
-    fig.savefig('plain_stat_diffs_{0}.pdf'.format(opt.precision), format = 'pdf')
-
-    # plot trajectory differences
-    for i in range(c0.particle_species):
-        fig = plt.figure(figsize=(12, 4))
-        for j in range(3):
-            a = fig.add_subplot(131 + j)
-            for t in range(c0.parameters['nparticles']):
-                a.plot(c0.trajectories[i][:, t, j] - c1.trajectories[i][:, t, j])
-        fig.savefig('traj_s{0}_{1}.pdf'.format(i, opt.precision), format = 'pdf')
+    compare_stats(opt, c0, c1)
     return None
 
 if __name__ == '__main__':
-    opt = parser.parse_args()
+    opt = parser.parse_args(
+            ['-n', '32',
+             '--run',
+             '--ncpu', '2',
+             '--nparticles', '1000',
+             '--niter_todo', '48',
+             '--precision', 'single',
+             '--multiplejob',
+             '--wd', 'data/single'] +
+            sys.argv[1:])
     plain(opt)
 
diff --git a/tests/test_scaling.py b/tests/test_scaling.py
new file mode 100755
index 0000000000000000000000000000000000000000..a5e675ccc5fe3dccf1e41b4270591d681dceb9d9
--- /dev/null
+++ b/tests/test_scaling.py
@@ -0,0 +1,64 @@
+#! /usr/bin/env python
+#######################################################################
+#                                                                     #
+#  Copyright 2015 Max Planck Institute                                #
+#                 for Dynamics and Self-Organization                  #
+#                                                                     #
+#  This file is part of bfps.                                         #
+#                                                                     #
+#  bfps is free software: you can redistribute it and/or modify       #
+#  it under the terms of the GNU General Public License as published  #
+#  by the Free Software Foundation, either version 3 of the License,  #
+#  or (at your option) any later version.                             #
+#                                                                     #
+#  bfps is distributed in the hope that it will be useful,            #
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of     #
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the      #
+#  GNU General Public License for more details.                       #
+#                                                                     #
+#  You should have received a copy of the GNU General Public License  #
+#  along with bfps.  If not, see <http://www.gnu.org/licenses/>       #
+#                                                                     #
+# Contact: Cristian.Lalescu@ds.mpg.de                                 #
+#                                                                     #
+#######################################################################
+
+
+
+from base import *
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+def scaling(opt):
+    wd = opt.work_dir
+    opt.work_dir = wd + '/N{0:0>3x}_1'.format(opt.n)
+    c0 = launch(opt, dt = 0.2/opt.n)
+    c0.compute_statistics()
+    print ('Re = {0:.0f}'.format(c0.statistics['Re']))
+    print ('Rlambda = {0:.0f}'.format(c0.statistics['Rlambda']))
+    print ('Lint = {0:.4e}, etaK = {1:.4e}'.format(c0.statistics['Lint'], c0.statistics['etaK']))
+    print ('Tint = {0:.4e}, tauK = {1:.4e}'.format(c0.statistics['Tint'], c0.statistics['tauK']))
+    print ('kMetaK = {0:.4e}'.format(c0.statistics['kMeta']))
+    for s in range(c0.particle_species):
+        acceleration_test(c0, species = s, m = 1)
+    opt.work_dir = wd + '/N{0:0>3x}_2'.format(opt.n)
+    opt.ncpu *= 2
+    c1 = launch(opt, dt = c0.parameters['dt'])
+    c1.compute_statistics()
+    compare_stats(opt, c0, c1)
+    return None
+
+if __name__ == '__main__':
+    opt = parser.parse_args(
+            ['-n', '32',
+             '--run',
+             '--initialize',
+             '--ncpu', '4',
+             '--nparticles', '10000',
+             '--niter_todo', '32',
+             '--precision', 'single',
+             '--wd', 'data/single'] +
+            sys.argv[1:])
+    scaling(opt)
+
diff --git a/todo.txt b/todo.txt
index 71c2677e5fa6af740a768f8fbee676c3357e1a75..2356510f3170e2841465fa2231af6b11411fcd5b 100644
--- a/todo.txt
+++ b/todo.txt
@@ -1,8 +1,8 @@
 (B) FFTW interpolator doesn't need its own field            @optimization +v1.0
+(B) clean up tox files, make sure all tests run             @tests        +v1.0
 (B) compute z polynomials only when needed                  @optimization +v1.0
-(B) get rid of temporal interpolation                       @optimization +v1.0
-(B) move get_grid coords to interpolator                    @optimization +v1.0
 (B) use less memory                                         @optimization
+(B) split library into core and extra                       @optimization +v1.0
 (C) clean up machine_settings mess                          @design @documentation +v2.0
 (C) code overview                                           @documentation
 (C) move stat I/O to cpp lib                                @design @HDF5
diff --git a/tox.ini b/tox.ini
index bc8543963e76e6b57db1e8558a8b7d8182bd8aff..8aca14f12a8a95079220ae5a70ebdcac32fd8737 100644
--- a/tox.ini
+++ b/tox.ini
@@ -1,6 +1,7 @@
 [tox]
-envlist = py27
+envlist = py34
 [testenv]
+sitepackages = True
 whitelist_externals =
     echo
     cp
@@ -9,43 +10,4 @@ changedir =
     {envtmpdir}
 commands =
     cp -r {toxinidir}/tests {envtmpdir}
-    python tests/test_io.py -n 32 --run --initialize --ncpu 2
-    python tests/test_plain.py -n 32 --run --initialize --multiplejob --ncpu 2 --nparticles 16 --niter_todo 64 --precision single --wd "data/single"
-    python tests/test_plain.py -n 32 --run --initialize --multiplejob --ncpu 2 --nparticles 16 --niter_todo 64 --precision double --wd "data/double"
-    python tests/test_particles.py \
-        -n 16 \
-        --run \
-        --initialize \
-        --ncpu 2 \
-        --nparticles 32 \
-        --niter_todo 16 \
-        --precision single \
-        --wd "data/particles_single"
-    python tests/test_particles.py \
-        -n 16 \
-        --run \
-        --initialize \
-        --ncpu 2 \
-        --nparticles 32 \
-        --niter_todo 16 \
-        --precision double \
-        --wd "data/particles_double"
-    python tests/test_convergence.py \
-        -n 32 \
-        --run \
-        --initialize \
-        --ncpu 2 \
-        --nparticles 32 \
-        --niter_todo 16 \
-        --precision single \
-        --wd "data/convergence_single"
-    python tests/test_convergence.py \
-        -n 32 \
-        --run \
-        --initialize \
-        --ncpu 2 \
-        --nparticles 32 \
-        --niter_todo 16 \
-        --precision double \
-        --wd "data/convergence_double"
-
+    python tests/{posargs:DEFAULTS}
diff --git a/tox_plain.ini b/tox_plain.ini
index 1aad555519ed22f2e3c678b96639089120dc7f8c..38f34b15fee4d479f34cbf805adebb61be8ec218 100644
--- a/tox_plain.ini
+++ b/tox_plain.ini
@@ -1,8 +1,7 @@
 [tox]
-envlist = py27
+envlist = py34
 [testenv]
 sitepackages = True
-#deps = matplotlib
 whitelist_externals =
     echo
     cp