     invert_invgammaN(invgammaNU,target,src); //=ne-1 = Gamma (ni-1)    
 }
 
-
+/*
+ * Interpolate physical values at the probe positions. Update probe
+ * vectors.
+ *
+ */
+template<class Matrix, class container, class P>
+void Feltor<Matrix, container, P> :: update_probes()
+{
+    static container probe_values(8); 
+    // Electron density
+    dg::blas2::gemv(probeinterp, npe[0], probe_values);
+    probevec[0] = probe_values;
+
+    // Electric potential
+    dg::blas2::gemv(probeinterp, phi[0], probe_values);
+    probevec[1] = probe_values;
+
+    // Radial particle transport
+    dg::blas2::gemv(poisson.dyrhs(), phi[0], omega);
+    dg::blas1::pointwiseDot(omega, binv, omega); // 1/B dy phi
+    dg::blas1::pointwiseDot(npe[0], omega, omega); // 1/B N dy phi
+    dg::blas2::gemv(probeinterp, omega, probe_values);
+    probevec[2] = probe_values;
+
+}
 
 
 template<class Matrix, class container, class P>

    initializene(profNi,profne); //ne = Gamma N_i
    dg::blas1::transform(profne,profne, dg::PLUS<>(+(p.bgprofamp + p.nprofileamp))); 
    dg::blas1::transform(profNi,profNi, dg::PLUS<>(+(p.bgprofamp + p.nprofileamp))); 
}

template<class Matrix, class container, class P>

#include "parameters.h"
// #include "geometry_circ.h"
#include "dg/average.h"
#include "dg/backend/linspace.cuh"

#ifdef DG_BENCHMARK
#include "dg/backend/timer.cuh"
