diff --git a/bfps/cpp/rFFTW_interpolator.cpp b/bfps/cpp/rFFTW_interpolator.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0968b20fd32b05c82afbc35267a8b60f290e0761
--- /dev/null
+++ b/bfps/cpp/rFFTW_interpolator.cpp
@@ -0,0 +1,144 @@
+/**********************************************************************
+* *
+* 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 *
+* *
+**********************************************************************/
+
+
+
+#include "rFFTW_interpolator.hpp"
+
+template <class rnumber, int interp_neighbours>
+rFFTW_interpolator<rnumber, interp_neighbours>::rFFTW_interpolator(
+ fluid_solver_base<rnumber> *fs,
+ base_polynomial_values BETA_POLYS)
+{
+ int tdims[4];
+ this->descriptor = fs->rd;
+ 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));
+ }
+ 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;
+}
+
+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));
+ }
+ else if (sizeof(rnumber) == 8)
+ {
+ fftw_free((double*)((void*)this->f0));
+ fftw_free((double*)((void*)this->f1));
+ }
+ delete this->descriptor;
+}
+
+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->unbuffered_descriptor->local_size,
+ dst + this->buffer_size);
+ return EXIT_SUCCESS;
+}
+
+template <class rnumber, int interp_neighbours>
+void rFFTW_interpolator<rnumber, interp_neighbours>::operator()(
+ double t,
+ int *xg,
+ double *xx,
+ double *dest,
+ int *deriv)
+{
+ double bx[interp_neighbours*2+2], by[interp_neighbours*2+2], bz[interp_neighbours*2+2];
+ if (deriv == NULL)
+ {
+ this->compute_beta(0, xx[0], bx);
+ this->compute_beta(0, xx[1], by);
+ this->compute_beta(0, xx[2], bz);
+ }
+ else
+ {
+ this->compute_beta(deriv[0], xx[0], bx);
+ this->compute_beta(deriv[1], xx[1], by);
+ this->compute_beta(deriv[2], xx[2], bz);
+ }
+ std::fill_n(dest, 3, 0);
+ ptrdiff_t bigiz, bigiy, bigix;
+ 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])
+ for (int iy = -interp_neighbours; iy <= interp_neighbours+1; iy++)
+ {
+ bigiy = ptrdiff_t(MOD(xg[1]+iy, this->descriptor->sizes[1]));
+ for (int ix = -interp_neighbours; ix <= interp_neighbours+1; ix++)
+ {
+ bigix = ptrdiff_t(MOD(xg[0]+ix, this->descriptor->sizes[2]));
+ for (int c=0; c<3; c++)
+ {
+ ptrdiff_t tindex = ((bigiz *this->descriptor->sizes[1] +
+ bigiy)*this->descriptor->sizes[2] +
+ bigix)*3+c + this->buffer_size;
+ dest[c] += (this->f0[tindex]*(1-t) + t*this->f1[tindex])*(bz[iz+interp_neighbours]*
+ by[iy+interp_neighbours]*
+ bx[ix+interp_neighbours]);
+ }
+ }
+ }
+ }
+}
+
+template class rFFTW_interpolator<float, 1>;
+template class rFFTW_interpolator<float, 2>;
+template class rFFTW_interpolator<float, 3>;
+template class rFFTW_interpolator<float, 4>;
+template class rFFTW_interpolator<float, 5>;
+template class rFFTW_interpolator<float, 6>;
+template class rFFTW_interpolator<double, 1>;
+template class rFFTW_interpolator<double, 2>;
+template class rFFTW_interpolator<double, 3>;
+template class rFFTW_interpolator<double, 4>;
+template class rFFTW_interpolator<double, 5>;
+template class rFFTW_interpolator<double, 6>;
+
diff --git a/bfps/cpp/rFFTW_interpolator.hpp b/bfps/cpp/rFFTW_interpolator.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..59b8b0c0b6895829cd889bce1a6b74197308451f
--- /dev/null
+++ b/bfps/cpp/rFFTW_interpolator.hpp
@@ -0,0 +1,62 @@
+/**********************************************************************
+* *
+* 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 *
+* *
+**********************************************************************/
+
+
+
+#include "field_descriptor.hpp"
+#include "fftw_tools.hpp"
+#include "fluid_solver_base.hpp"
+#include "spline_n1.hpp"
+#include "spline_n2.hpp"
+#include "spline_n3.hpp"
+#include "spline_n4.hpp"
+#include "spline_n5.hpp"
+#include "spline_n6.hpp"
+#include "Lagrange_polys.hpp"
+#include "interpolator.hpp"
+
+#ifndef RFFTW_INTERPOLATOR
+
+#define RFFTW_INTERPOLATOR
+
+template <class rnumber, int interp_neighbours>
+class rFFTW_interpolator
+{
+ public:
+ ptrdiff_t field_size;
+ base_polynomial_values compute_beta;
+ field_descriptor<rnumber> *descriptor;
+ rnumber *f0, *f1, *temp;
+
+ interpolator(
+ fluid_solver_base<rnumber> *FSOLVER,
+ base_polynomial_values BETA_POLYS);
+ ~interpolator();
+
+ void operator()(double t, int *__restrict__ xg, double *__restrict__ xx, double *__restrict__ dest, int *deriv = NULL);
+ int read_rFFTW(void *src);
+};
+
+#endif//RFFTW_INTERPOLATOR
+