diff --git a/python/_fast_interp.pyx b/python/_fast_interp.pyx
index 8e48a7b..b13b396 100644
--- a/python/_fast_interp.pyx
+++ b/python/_fast_interp.pyx
@@ -11,22 +11,33 @@ __all__=["fast_interp"]
 
 @cython.boundscheck(False)
 @cython.cdivision(True)
-def fast_interp(npx.float64_t xmin, npx.float64_t dx, npx.float64_t[:] A, npx.float64_t[:] y, npx.float64_t[:] out, npx.float64_t beyond_val=npx.nan):
+def fast_interp(xmin0, dx0, A0, y0, out0, beyond_val=np.nan):
   cdef double rq, q
   cdef int iq
   cdef long i, Asize, ysize
-  cdef npx.float64_t[:] out0
+  cdef npx.float64_t xmin, dx
+  cdef npx.float64_t[:] out
+  cdef npx.float64_t[:] A, y
+  cdef npx.float64_t beyond
 
+  beyond=beyond_val
+  xmin = xmin0
+  dx = dx0
+  A = A0
+  y = y0
   ysize = y.size
-  out0 = out
+  out = out0
   Asize = A.size
 
+  if out.size != ysize:
+    raise ValueError("out and y must have the same size")
+
   with nogil:
     for i in prange(ysize):
       q = (y[i] - xmin) / dx
       iq = int(floor(q))
       rq = (q-iq)
       if iq+1 >= Asize or iq < 0:
-        out0[i] = beyond_val
+        out[i] = beyond
       else:
-        out0[i] = rq * A[iq+1] + (1-rq)*A[iq]
+        out[i] = rq * A[iq+1] + (1-rq)*A[iq]
diff --git a/python_sample/gen_2lpt_density.py b/python_sample/gen_2lpt_density.py
index 85a23e7..76636db 100644
--- a/python_sample/gen_2lpt_density.py
+++ b/python_sample/gen_2lpt_density.py
@@ -1,3 +1,4 @@
+import numexpr as ne
 import numpy as np
 import cosmotool as ct
 import icgen as bic
@@ -25,6 +26,7 @@ parser.add_argument('--base', type=str, required=True)
 parser.add_argument('--N', type=int, default=256)
 parser.add_argument('--output', type=str, default="dcic_%d.npy")
 parser.add_argument('--supersample', type=int, default=1)
+parser.add_argument('--rsd', action='store_true') 
 args = parser.parse_args()
 
 
@@ -32,8 +34,19 @@ args = parser.parse_args()
 for i in xrange(args.start, args.end, args.step):
   print i
 #  pos,_,density,N,L,_ = bic.run_generation("/nethome/lavaux/remote/borg_2m++_128/initial_density_%d.dat" % i, 0.001, astart, cosmo, supersample=2, do_lpt2=True)
-  pos,_,density,N,L,_,_ = bic.run_generation("%s/initial_density_%d.dat" % (args.base,i), 0.001, astart, 
-                                           cosmo, supersample=args.supersample, do_lpt2=True, needvel=False)
+  pos,vel,density,N,L,_,_ = bic.run_generation("%s/initial_density_%d.dat" % (args.base,i), 0.001, astart, 
+                                           cosmo, supersample=args.supersample, do_lpt2=True, needvel=True)
+
+  if args.rsd:
+    inv_r2 = ne.evaluate('1/sqrt(x**2+y**2+z**2)', 
+        local_dict={'x':pos[0], 'y':pos[1], 'z':pos[2]})
+    rsd = lambda p,v: ne.evaluate('x + (x*vx)*inv_r2 / H', 
+        local_dict={'x':p, 'inv_r2':inv_r2,
+         'vx':v, 'H':100.0}, out=p, casting='unsafe')
+
+    rsd(pos[0], vel[0])
+    rsd(pos[1], vel[1])
+    rsd(pos[2], vel[2])
 
   dcic = ct.cicParticles(pos, L, args.N)
   dcic /= np.average(np.average(np.average(dcic, axis=0), axis=0), axis=0)
diff --git a/python_sample/icgen/borgicgen.py b/python_sample/icgen/borgicgen.py
index 706f2e4..78f1805 100644
--- a/python_sample/icgen/borgicgen.py
+++ b/python_sample/icgen/borgicgen.py
@@ -140,7 +140,7 @@ def run_generation(input_borg, a_borg, a_ic, cosmo, supersample=1, supergenerate
     lpt = LagrangianPerturbation(-density, L, fourier=True, supersample=supersample)
 
     # Generate grid
-    posq = gen_posgrid(N*supergenerate, L)
+    posq = gen_posgrid(N*supersample, L)
     vel= []
     posx = []
     
diff --git a/src/fourier/fft/fftw_calls_mpi.hpp b/src/fourier/fft/fftw_calls_mpi.hpp
index 410c930..fe3506e 100644
--- a/src/fourier/fft/fftw_calls_mpi.hpp
+++ b/src/fourier/fft/fftw_calls_mpi.hpp
@@ -46,7 +46,7 @@ public: \
   static void execute(plan_type p) { prefix ## _execute(p); } \
   static void execute_r2c(plan_type p, real_type *in, complex_type *out) { prefix ## _mpi_execute_dft_r2c(p, in, out); } \
   static void execute_c2r(plan_type p, std::complex<real_type> *in, real_type *out) { prefix ## _mpi_execute_dft_c2r(p, (complex_type*)in, out); } \
-  static void execute_r2c(plan_type p, real_type *in, std::compelex<real_type> *out) { prefix ## _mpi_execute_dft_r2c(p, in, (complex_type*)out); } \
+  static void execute_r2c(plan_type p, real_type *in, std::complex<real_type> *out) { prefix ## _mpi_execute_dft_r2c(p, in, (complex_type*)out); } \
 \
   static plan_type plan_dft_r2c_2d(int Nx, int Ny,  \
 		       real_type *in, complex_type *out, \
diff --git a/src/loadGadget.cpp b/src/loadGadget.cpp
index e307dbf..474cf9e 100644
--- a/src/loadGadget.cpp
+++ b/src/loadGadget.cpp
@@ -146,12 +146,12 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
       }
       
       if (h.flag_doubleprecision) {
-        cout << "Gadget format with double precision" << endl;
+        //cout << "Gadget format with double precision" << endl;
         readToDouble = boost::bind(myRead64<double>, f);
         readToSingle = boost::bind(myRead64<float>, f);
         float_size = sizeof(double);
       } else {
-        cout << "Gadget format with single precision" << endl;
+        //cout << "Gadget format with single precision" << endl;
         readToDouble = boost::bind(myRead32<double>, f);
         readToSingle = boost::bind(myRead32<float>, f);
         float_size = sizeof(float);
@@ -306,7 +306,6 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
   if (loadflags & NEED_MASS) {
     bool do_load = false;
 
-    cout << "=> Mass" << endl;
     for (int k = 0; k < 6; k++)
       {
          do_load = do_load || ((h.mass[k] == 0)&&(h.npart[k]>0));