diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 9b619a7c66a011400466213591acf8e60dcff74e..e89c26197a74b20d0abec37bb8cb3a9bfc218b7a 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -32,8 +32,8 @@ image: gitlab-registry.mpcdf.mpg.de/mpcdf/module-image:latest
     export CXX=${GCC_HOME}/bin/g++
 
 .export_INTEL_compilers: &export_INTEL_compilers |
-    export CC=icc
-    export CXX=icpc
+    export CC=icx
+    export CXX=icpx
 
 #.allow_docker_run_as_root: &allow_docker_run_as_root |
 #    export OMPI_ALLOW_RUN_AS_ROOT=1
@@ -49,7 +49,7 @@ build-gcc-openmpi:
       COMPILER: "gcc/11"
       MPI: "openmpi/4"
     tags:
-      - docker
+      - docker, 1core
     artifacts:
         paths:
             - build/
@@ -68,7 +68,7 @@ build-gcc-impi:
       COMPILER: "gcc/11"
       MPI: "impi/2021.7"
     tags:
-      - docker
+      - docker, 1core
     artifacts:
         paths:
             - build/
@@ -84,10 +84,10 @@ build-intel:
         export MPI_HOME=${I_MPI_ROOT}
       - *build
     variables:
-      COMPILER: "intel/21.7.1"
-      MPI: "impi/2021.7"
+      COMPILER: "intel/2023.1.0.x"
+      MPI: "impi/2021.9"
     tags:
-      - docker
+      - docker, 1core
     artifacts:
         paths:
             - build/
@@ -138,8 +138,8 @@ test-intel:
     tags:
       - docker
     variables:
-      COMPILER: "intel/21.7.1"
-      MPI: "impi/2021.7"
+      COMPILER: "intel/2023.1.0.x"
+      MPI: "impi/2021.9"
     needs:
         - job: build-intel
           artifacts: true
@@ -158,7 +158,7 @@ build-doc:
       - make VERBOSE=1 doc_html_full
         #- make VERBOSE=1 doc_latex
     tags:
-      - docker
+      - docker, 1core
     variables:
       COMPILER: "gcc/11"
       MPI: "impi/2021.7"
@@ -202,7 +202,7 @@ pages:
 #    only:
 #      - tags
     tags:
-      - docker
+      - docker, 1core
     needs:
         - job: build-doc
           artifacts: true
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9a16ec039f9b6b3569b4c693296bcd3e6c0ba037..bd7846bc24d58d6723cc4dc8a31849001b4ad652 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -96,29 +96,31 @@ target_link_libraries(TurTLE
 
 
 #####################################################################################
+## CMake stuff
 if (DEFINED ENV{CMAKE_INSTALL_PREFIX})
     set(CMAKE_INSTALL_PREFIX $ENV{CMAKE_INSTALL_PREFIX})
 endif()
 
-
 set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake ${CMAKE_MODULE_PATH})
-#####################################################################################
-
 
-#####################################################################################
 if(NOT CMAKE_BUILD_TYPE)
     set(CMAKE_BUILD_TYPE Release CACHE STRING "Build type" FORCE)
 endif()
 #####################################################################################
 
 #####################################################################################
-## Miscellaneous options
+## TurTLE options
 option(NDEBUG "Define NDEBUG macro" ON)
 if(NDEBUG)
     add_definitions(-DNDEBUG)
 endif()
 
 option(USE_TIMING_OUTPUT "Toggle timing output. WARNING: memory usage is proportional to `niter_todo`" OFF)
+
+set(MIN_INTERPOLATION_NEIGHBOURS 1 CACHE STRING "Minimum interpolation neighbours available to `particles_system`")
+set(MIN_INTERPOLATION_SMOOTHNESS 1 CACHE STRING "Minimum interpolation smoothness available to `particles_system`")
+set(MAX_INTERPOLATION_NEIGHBOURS 2 CACHE STRING "Maximum interpolation neighbours available to `particles_system`")
+set(MAX_INTERPOLATION_SMOOTHNESS 2 CACHE STRING "Maximum interpolation smoothness available to `particles_system`")
 #####################################################################################
 
 
diff --git a/TurTLE/DNS.py b/TurTLE/DNS.py
index 40d1da539a798c0add8fbef311b1b2d4de5765c0..896c9ba8954a9b3ad6b72e98b573c5f053c4ffad 100644
--- a/TurTLE/DNS.py
+++ b/TurTLE/DNS.py
@@ -88,19 +88,7 @@ class DNS(_code):
                 '* version {0}\n'.format(TurTLE.__version__) +
                 '***********************************************************************/\n\n\n')
         self.include_list = [
-                '"base.hpp"',
-                '"scope_timer.hpp"',
-                '"fftw_interface.hpp"',
                 '"full_code/main_code.hpp"',
-                '<cmath>',
-                '<iostream>',
-                '<hdf5.h>',
-                '<string>',
-                '<cstring>',
-                '<fftw3-mpi.h>',
-                '<omp.h>',
-                '<cfenv>',
-                '<cstdlib>',
                 '"full_code/{0}.hpp"\n'.format(self.dns_type)]
         self.main = """
             int main(int argc, char *argv[])
diff --git a/TurTLE/PP.py b/TurTLE/PP.py
index d77d550fab534e8337f1a28da3cf170bb747957b..6b86bb41930f8bb8597f0cb4df9a6328e0314e51 100644
--- a/TurTLE/PP.py
+++ b/TurTLE/PP.py
@@ -80,19 +80,8 @@ class PP(_code):
                 '* version {0}\n'.format(TurTLE.__version__) +
                 '***********************************************************************/\n\n\n')
         self.include_list = [
-                '"base.hpp"',
-                '"scope_timer.hpp"',
-                '"fftw_interface.hpp"',
                 '"full_code/main_code.hpp"',
                 '<cmath>',
-                '<iostream>',
-                '<hdf5.h>',
-                '<string>',
-                '<cstring>',
-                '<fftw3-mpi.h>',
-                '<omp.h>',
-                '<cfenv>',
-                '<cstdlib>',
                 '"full_code/{0}.hpp"\n'.format(self.dns_type)]
         self.main = """
             int main(int argc, char *argv[])
diff --git a/TurTLE/TEST.py b/TurTLE/TEST.py
index 269cdef17e0902c696ffd7b1713f995c7665a613..42522472751b2871624989d9ab69e2ea686653c5 100644
--- a/TurTLE/TEST.py
+++ b/TurTLE/TEST.py
@@ -80,19 +80,8 @@ class TEST(_code):
                 '* version {0}\n'.format(TurTLE.__version__) +
                 '***********************************************************************/\n\n\n')
         self.include_list = [
-                '"base.hpp"',
-                '"scope_timer.hpp"',
-                '"fftw_interface.hpp"',
                 '"full_code/main_code.hpp"',
                 '<cmath>',
-                '<iostream>',
-                '<hdf5.h>',
-                '<string>',
-                '<cstring>',
-                '<fftw3-mpi.h>',
-                '<omp.h>',
-                '<cfenv>',
-                '<cstdlib>',
                 '"full_code/{0}.hpp"\n'.format(self.dns_type)]
         self.main = """
             int main(int argc, char *argv[])
diff --git a/cpp/base.cpp b/cpp/base.cpp
index 1ad16d3794e7632a3e107f04b65320bf8283c85b..94d53d94bb866dde70d5af7d0e626d09f9ea5ab6 100644
--- a/cpp/base.cpp
+++ b/cpp/base.cpp
@@ -1,5 +1,7 @@
 #include "base.hpp"
 
+#include <iostream>
+
 int maximum_MPI_tag_value;
 
 int read_maximum_MPI_tag_value()
@@ -23,3 +25,6 @@ int read_maximum_MPI_tag_value()
     }
     return EXIT_SUCCESS;
 }
+
+/////////////////////////////////////////////////////////////
+
diff --git a/cpp/base.hpp b/cpp/base.hpp
index ca33bf6876048864a0e148d8f7f14a18fa908e69..d098f797782bcba04b7efbeb883af49d3cc9504c 100644
--- a/cpp/base.hpp
+++ b/cpp/base.hpp
@@ -28,13 +28,14 @@
 #define BASE_HPP
 
 
+#include "turtle_config.hpp"
+
 #include <cassert>
 #include <mpi.h>
 #include <stdarg.h>
-#include <iostream>
 #include <typeinfo>
+#include <cstdlib>
 
-#include "turtle_config.hpp"
 
 /////////////////////////////////////////////////////////////
 /////////////////////////////////////////////////////////////
@@ -167,12 +168,13 @@ class mpi_real_type<long long int>
 
 #ifndef NDEBUG
 
-static char debug_message_buffer[message_buffer_length];
+#include <iostream>
 
 inline void DEBUG_MSG(const char * format, ...)
 {
     va_list argptr;
     va_start(argptr, format);
+    static char debug_message_buffer[message_buffer_length];
     sprintf(
             debug_message_buffer,
             "MPIrank%.4d ",
@@ -190,6 +192,7 @@ inline void DEBUG_MSG_WAIT(MPI_Comm communicator, const char * format, ...)
 {
     va_list argptr;
     va_start(argptr, format);
+    static char debug_message_buffer[message_buffer_length];
     sprintf(
             debug_message_buffer,
             "MPIrank%.4d ",
diff --git a/cpp/env_utils.hpp b/cpp/env_utils.hpp
index f030c7c8a4fec908fff0b2f7514ead49d151c416..e3cd73052d72109b62a92a86841cff363c0f45bb 100644
--- a/cpp/env_utils.hpp
+++ b/cpp/env_utils.hpp
@@ -29,10 +29,7 @@
 
 #include <cstdlib>
 #include <sstream>
-#include <iostream>
 #include <cstring>
-#include <cstring>
-#include <array>
 
 /** \class env_utils
  * \brief utilities
diff --git a/cpp/fftw_interface.hpp b/cpp/fftw_interface.hpp
index 8528aba39dbe56855181ae26be9fc1b984171aa5..3dd4c74d506cdd676aa677e2d670921b74f9c595 100644
--- a/cpp/fftw_interface.hpp
+++ b/cpp/fftw_interface.hpp
@@ -26,7 +26,6 @@
 #define FFTW_INTERFACE_HPP
 
 #include <fftw3-mpi.h>
-#include <map>
 #include <string>
 
 #ifdef USE_FFTWESTIMATE
@@ -42,7 +41,6 @@
 #ifdef SPLIT_FFTW_MANY
 #include <vector>
 #include <memory>
-#include <algorithm>
 #include <cassert>
 #include <cstring>
 #include <type_traits>
diff --git a/cpp/fftw_tools.cpp b/cpp/fftw_tools.cpp
index ae03966d686baab544d865cf1eab253faea6584d..6ee6338ffa519b70bc27dca2ac6050fb83638db5 100644
--- a/cpp/fftw_tools.cpp
+++ b/cpp/fftw_tools.cpp
@@ -22,9 +22,6 @@
 *                                                                     *
 **********************************************************************/
 
-#include <stdlib.h>
-#include <algorithm>
-#include <iostream>
 #include "base.hpp"
 #include "fftw_tools.hpp"
 #include "fftw_interface.hpp"
diff --git a/cpp/fftw_tools.hpp b/cpp/fftw_tools.hpp
index 5fbcdf81bb3e3e82def8e897ed80a8c803584282..064202ccf2a0758b6535ea42b755f4d91f1295fb 100644
--- a/cpp/fftw_tools.hpp
+++ b/cpp/fftw_tools.hpp
@@ -28,12 +28,15 @@
 #define FFTW_TOOLS
 
 
-#include <mpi.h>
 #include <fftw3-mpi.h>
 #include <string>
 #include <map>
 #include <cmath>
 
+extern int myrank, nprocs;
+
+extern std::map<std::string, unsigned> fftw_planner_string_to_flag;
+
 
 // helper function
 inline void complex_number_phase_shifter(double &xx, double &yy, const double cos_val, const double sin_val)
@@ -46,8 +49,8 @@ inline void complex_number_phase_shifter(double &xx, double &yy, const double co
 
 inline void complex_number_phase_shifter(double &xx, double &yy, const double phase)
 {
-    const double cval = cos(phase);
-    const double sval = sin(phase);
+    const double cval = std::cos(phase);
+    const double sval = std::sin(phase);
     complex_number_phase_shifter(xx, yy, cval, sval);
 }
 
@@ -61,14 +64,10 @@ inline void complex_number_phase_shifter(float &xx, float &yy, const double cos_
 
 inline void complex_number_phase_shifter(float &xx, float &yy, const double phase)
 {
-    const double cval = cos(phase);
-    const double sval = sin(phase);
+    const double cval = std::cos(phase);
+    const double sval = std::sin(phase);
     complex_number_phase_shifter(xx, yy, cval, sval);
 }
 
-extern int myrank, nprocs;
-
-extern std::map<std::string, unsigned> fftw_planner_string_to_flag;
-
 #endif//FFTW_TOOLS
 
diff --git a/cpp/field.cpp b/cpp/field.cpp
index 91905c9457ac81427240bbe69c061e064a1a0249..e1147209781ed7cebd5d76e2fe78ef2ba576eedb 100644
--- a/cpp/field.cpp
+++ b/cpp/field.cpp
@@ -24,17 +24,13 @@
 
 
 
-#include "field.hpp"
 #include "fftw_tools.hpp"
 #include "scope_timer.hpp"
 #include "shared_array.hpp"
 #include "field_binary_IO.hpp"
 
 #include <sys/stat.h>
-#include <cmath>
-#include <cstdlib>
 #include <algorithm>
-#include <cassert>
 #include <random>
 
 template <typename rnumber,
@@ -156,6 +152,39 @@ field<rnumber, be, fc>::~field()
     }
 }
 
+template <typename rnumber,
+          field_backend be,
+          field_components fc>
+int field<rnumber, be, fc>::print_plan(const std::string preamble)
+{
+    char *c2r_plan_information = fftw_interface<rnumber>::sprint(this->c2r_plan);
+    char *r2c_plan_information = fftw_interface<rnumber>::sprint(this->r2c_plan);
+    if (this->myrank == 0)
+    {
+        std::cout << preamble <<
+                     std::endl <<
+                     "----c2r plan is:\n" <<
+                     c2r_plan_information <<
+                     std::endl <<
+                     "----r2c plan is:\n" <<
+                     r2c_plan_information <<
+                     std::endl;
+    }
+    std::string err_message = (
+            std::string("MPI rank ") +
+            preamble +
+            std::to_string(this->myrank) +
+            std::string("\n----c2r plan is:\n") +
+            std::string(c2r_plan_information) +
+            std::string("\n----r2c plan is:\n") +
+            std::string(r2c_plan_information) +
+            std::string("\n"));
+    std::cerr << err_message;
+    free(c2r_plan_information);
+    free(r2c_plan_information);
+    return EXIT_SUCCESS;
+}
+
 template <typename rnumber,
           field_backend be,
           field_components fc>
@@ -1170,7 +1199,7 @@ void field<rnumber, be, fc>::compute_rspace_stats(
             }
             if (nvals == int(4))
             {
-                val_tmp[3] = sqrt(val_tmp[3]);
+                val_tmp[3] = std::sqrt(val_tmp[3]);
                 if (val_tmp[3] < local_moments[0*nvals+3])
                     local_moments[0*nvals+3] = val_tmp[3];
                 if (val_tmp[3] > local_moments[9*nvals+3])
@@ -1882,22 +1911,22 @@ void field<rnumber, be, fc>::symmetrize_alternate()
                     ptrdiff_t cindex = this->get_cindex(0, iyy, iz);
                     for (int cc = 0; cc < int(ncomp(fc)); cc++)
                     {
-                        const double ampp = sqrt(
+                        const double ampp = std::sqrt(
                                 (*(bufferp + ncomp(fc)*iz+cc))[0]*(*(bufferp + ncomp(fc)*iz+cc))[0] +
                                 (*(bufferp + ncomp(fc)*iz+cc))[1]*(*(bufferp + ncomp(fc)*iz+cc))[1]);
                         const double phip = atan2(
                                 (*(bufferp + ncomp(fc)*iz+cc))[1],
                                 (*(bufferp + ncomp(fc)*iz+cc))[0]);
-                        const double ampm = sqrt(
+                        const double ampm = std::sqrt(
                                 (*(bufferm + ncomp(fc)*izz+cc))[0]*(*(bufferm + ncomp(fc)*izz+cc))[0] +
                                 (*(bufferm + ncomp(fc)*izz+cc))[1]*(*(bufferm + ncomp(fc)*izz+cc))[1]);
                         const double phim = atan2(
                                 (*(bufferm + ncomp(fc)*izz+cc))[1],
                                 (*(bufferm + ncomp(fc)*izz+cc))[0]);
-                        const double amp = sqrt(ampp*ampm);
+                        const double amp = std::sqrt(ampp*ampm);
                         const double phi = (phip - phim)/2;
-                        (*(cdata + ncomp(fc)*cindex + cc))[0] =  amp*cos(phi);
-                        (*(cdata + ncomp(fc)*cindex + cc))[1] =  amp*sin(phi);
+                        (*(cdata + ncomp(fc)*cindex + cc))[0] =  amp*std::cos(phi);
+                        (*(cdata + ncomp(fc)*cindex + cc))[1] =  amp*std::sin(phi);
                     }
                 }
             }
@@ -1912,22 +1941,22 @@ void field<rnumber, be, fc>::symmetrize_alternate()
                     ptrdiff_t cindex = this->get_cindex(0, iyy, izz);
                     for (int cc = 0; cc < int(ncomp(fc)); cc++)
                     {
-                        const double ampp = sqrt(
+                        const double ampp = std::sqrt(
                                 (*(bufferp + ncomp(fc)*iz+cc))[0]*(*(bufferp + ncomp(fc)*iz+cc))[0] +
                                 (*(bufferp + ncomp(fc)*iz+cc))[1]*(*(bufferp + ncomp(fc)*iz+cc))[1]);
                         const double phip = atan2(
                                 (*(bufferp + ncomp(fc)*iz+cc))[1],
                                 (*(bufferp + ncomp(fc)*iz+cc))[0]);
-                        const double ampm = sqrt(
+                        const double ampm = std::sqrt(
                                 (*(bufferm + ncomp(fc)*izz+cc))[0]*(*(bufferm + ncomp(fc)*izz+cc))[0] +
                                 (*(bufferm + ncomp(fc)*izz+cc))[1]*(*(bufferm + ncomp(fc)*izz+cc))[1]);
                         const double phim = atan2(
                                 (*(bufferm + ncomp(fc)*izz+cc))[1],
                                 (*(bufferm + ncomp(fc)*izz+cc))[0]);
-                        const double amp = sqrt(ampp*ampm);
+                        const double amp = std::sqrt(ampp*ampm);
                         const double phi = (phip - phim)/2;
-                        (*(cdata + ncomp(fc)*cindex + cc))[0] =  amp*cos(phi);
-                        (*(cdata + ncomp(fc)*cindex + cc))[1] = -amp*sin(phi);
+                        (*(cdata + ncomp(fc)*cindex + cc))[0] =  amp*std::cos(phi);
+                        (*(cdata + ncomp(fc)*cindex + cc))[1] = -amp*std::sin(phi);
                     }
                 }
             }
@@ -1947,24 +1976,24 @@ void field<rnumber, be, fc>::symmetrize_alternate()
             ptrdiff_t cindex1 = this->get_cindex(0, 0, this->clayout->sizes[1] - iz);
             for (int cc = 0; cc < int(ncomp(fc)); cc++)
             {
-                const double ampp = sqrt(
+                const double ampp = std::sqrt(
                         (*(cdata + cc + ncomp(fc)*cindex0))[0]*(*(cdata + cc + ncomp(fc)*cindex0))[0] +
                         (*(cdata + cc + ncomp(fc)*cindex0))[1]*(*(cdata + cc + ncomp(fc)*cindex0))[1]);
                 const double phip = atan2(
                         (*(cdata + cc + ncomp(fc)*cindex0))[1],
                         (*(cdata + cc + ncomp(fc)*cindex0))[0]);
-                const double ampm = sqrt(
+                const double ampm = std::sqrt(
                         (*(cdata + cc + ncomp(fc)*cindex1))[0]*(*(cdata + cc + ncomp(fc)*cindex1))[0] +
                         (*(cdata + cc + ncomp(fc)*cindex1))[1]*(*(cdata + cc + ncomp(fc)*cindex1))[1]);
                 const double phim = atan2(
                         (*(cdata + cc + ncomp(fc)*cindex1))[1],
                         (*(cdata + cc + ncomp(fc)*cindex1))[0]);
-                const double amp = sqrt(ampp*ampm);
+                const double amp = std::sqrt(ampp*ampm);
                 const double phi = (phip - phim)/2;
-                (*(cdata + cc + ncomp(fc)*cindex0))[0] =  amp*cos(phi);
-                (*(cdata + cc + ncomp(fc)*cindex0))[1] =  amp*sin(phi);
-                (*(cdata + cc + ncomp(fc)*cindex1))[0] =  amp*cos(phi);
-                (*(cdata + cc + ncomp(fc)*cindex1))[1] = -amp*sin(phi);
+                (*(cdata + cc + ncomp(fc)*cindex0))[0] =  amp*std::cos(phi);
+                (*(cdata + cc + ncomp(fc)*cindex0))[1] =  amp*std::sin(phi);
+                (*(cdata + cc + ncomp(fc)*cindex1))[0] =  amp*std::cos(phi);
+                (*(cdata + cc + ncomp(fc)*cindex1))[1] = -amp*std::sin(phi);
             }
         }
     }
@@ -2228,7 +2257,7 @@ void field<rnumber, be, fc>::compute_stats(
             break;
         case THREE:
             max_estimate_vector.resize(4, max_estimate);
-            max_estimate_vector[3] *= sqrt(3);
+            max_estimate_vector[3] *= std::sqrt(3);
             break;
         case THREExTHREE:
             max_estimate_vector.resize(9, max_estimate);
@@ -2303,7 +2332,7 @@ double field<rnumber, be, fc>::L2norm(
                 &m2,
                 1,
                 MPI_DOUBLE, MPI_SUM, this->comm);
-        return sqrt(m2 / this->npoints);
+        return std::sqrt(m2 / this->npoints);
     }
 }
 
@@ -2622,8 +2651,8 @@ int joint_rspace_PDF(
                             local_histc[(bin1*nbins + bin2)*9 + i*3 + j]++;
                     }
                 }
-                bin1 = int(floor(sqrt(mag1)/bin1size[3]));
-                bin2 = int(floor(sqrt(mag2)/bin2size[3]));
+                bin1 = int(floor(std::sqrt(mag1)/bin1size[3]));
+                bin2 = int(floor(std::sqrt(mag2)/bin2size[3]));
             }
             else if (fc == ONE)
             {
@@ -2912,7 +2941,6 @@ field<rnumber, be, fc> &field<rnumber, be, fc>::operator=(
         // use CLOOP, such that memory caching per thread stuff is not messed up
         #pragma omp parallel
         {
-            hsize_t npoints = 0;
             const hsize_t start = OmpUtils::ForIntervalStart(this->clayout->subsizes[1]);
             const hsize_t end = OmpUtils::ForIntervalEnd(this->clayout->subsizes[1]);
 
@@ -2923,7 +2951,6 @@ field<rnumber, be, fc> &field<rnumber, be, fc>::operator=(
                             zindex*this->clayout->subsizes[2]);
                     for (hsize_t xindex = 0; xindex < this->clayout->subsizes[2]; xindex++)
                     {
-                        npoints++;
                         for (unsigned int cc=0; cc < ncomp(fc); cc++)
                         {
                             *(this->data + 2*((cindex+xindex)*ncomp(fc) + cc)) = value;
@@ -3118,24 +3145,24 @@ int make_gaussian_random_field(
                         // Pope spectrum, Pope (2000), ''Turbulent flows'', p. 232
                         double C = 1.5;
                         double beta = 5.2;
-                        double f_L = pow(sqrt(k2)*Lint/sqrt(k2*pow(Lint, 2) + c_L), 11./3.);
-                        double f_eta = exp(-beta*(pow(pow(k2, 2)*pow(etaK, 4) + pow(c_eta, 4), 1./4.) - c_eta));
-                        double spectrum = coefficient*C*pow(dissipation, 2./3.)*pow(k2, -5./6.)*f_L*f_eta;
+                        double f_L = std::pow(std::sqrt(k2)*Lint/std::sqrt(k2*std::pow(Lint, 2) + c_L), 11./3.);
+                        double f_eta = exp(-beta*(std::pow(std::pow(k2, 2)*std::pow(etaK, 4) + std::pow(c_eta, 4), 1./4.) - c_eta));
+                        double spectrum = coefficient*C*std::pow(dissipation, 2./3.)*std::pow(k2, -5./6.)*f_L*f_eta;
                         // TODO: What about Delta k?
                         switch(fc)
                         {
                         case ONE:
                         {
-                            output_field->cval(cindex,0) = rdist(rgen[omp_get_thread_num()]) * sqrt(spectrum * kk->dkx*kk->dky*kk->dkz / (4.*M_PI*k2));
-                            output_field->cval(cindex,1) = rdist(rgen[omp_get_thread_num()]) * sqrt(spectrum * kk->dkx*kk->dky*kk->dkz / (4.*M_PI*k2));
+                            output_field->cval(cindex,0) = rdist(rgen[omp_get_thread_num()]) * std::sqrt(spectrum * kk->dkx*kk->dky*kk->dkz / (4.*M_PI*k2));
+                            output_field->cval(cindex,1) = rdist(rgen[omp_get_thread_num()]) * std::sqrt(spectrum * kk->dkx*kk->dky*kk->dkz / (4.*M_PI*k2));
                             break;
                         }
                         case THREE:
                         for (int cc = 0; cc<3; cc++)
                         {
                             // factor 3 compensates the trace between the spectral tensor and the energy spectrum
-                            output_field->cval(cindex,cc,0) = rdist(rgen[omp_get_thread_num()]) * sqrt(spectrum * kk->dkx*kk->dky*kk->dkz / 3. / (4.*M_PI*k2));
-                            output_field->cval(cindex,cc,1) = rdist(rgen[omp_get_thread_num()]) * sqrt(spectrum * kk->dkx*kk->dky*kk->dkz / 3. / (4.*M_PI*k2));
+                            output_field->cval(cindex,cc,0) = rdist(rgen[omp_get_thread_num()]) * std::sqrt(spectrum * kk->dkx*kk->dky*kk->dkz / 3. / (4.*M_PI*k2));
+                            output_field->cval(cindex,cc,1) = rdist(rgen[omp_get_thread_num()]) * std::sqrt(spectrum * kk->dkx*kk->dky*kk->dkz / 3. / (4.*M_PI*k2));
                         }
                             break;
                         case THREExTHREE:
@@ -3143,8 +3170,8 @@ int make_gaussian_random_field(
                         for (int ccc = 0; ccc<3; ccc++)
                         {
                             // factor 9 compensates the trace between the spectral tensor and the energy spectrum
-                            output_field->cval(cindex,cc,ccc,0) = rdist(rgen[omp_get_thread_num()]) * sqrt(spectrum * kk->dkx*kk->dky*kk->dkz / 9. / (4.*M_PI*k2));
-                            output_field->cval(cindex,cc,ccc,1) = rdist(rgen[omp_get_thread_num()]) * sqrt(spectrum * kk->dkx*kk->dky*kk->dkz / 9. / (4.*M_PI*k2));
+                            output_field->cval(cindex,cc,ccc,0) = rdist(rgen[omp_get_thread_num()]) * std::sqrt(spectrum * kk->dkx*kk->dky*kk->dkz / 9. / (4.*M_PI*k2));
+                            output_field->cval(cindex,cc,ccc,1) = rdist(rgen[omp_get_thread_num()]) * std::sqrt(spectrum * kk->dkx*kk->dky*kk->dkz / 9. / (4.*M_PI*k2));
                         }
                             break;
                         }
@@ -3165,7 +3192,7 @@ int generate_random_phase_field(
 {
     TIMEZONE("generate_random_phase_field");
     assert(output_field->real_space_representation == false);
-    const double twopi = 4*acos(0);
+    const double twopi = 4*std::acos(0.0);
 
     // initialize a separate random number generator for each thread
     std::vector<std::mt19937_64> rgen;
@@ -3195,8 +3222,8 @@ int generate_random_phase_field(
                     {
                         // single phase shift for all components
                         const double phi = rdist(rgen[omp_get_thread_num()]);
-                        output_field->cval(cindex, 0) = cos(phi);
-                        output_field->cval(cindex, 1) = sin(phi);
+                        output_field->cval(cindex, 0) = std::cos(phi);
+                        output_field->cval(cindex, 1) = std::sin(phi);
                     }
                     else
                     {
@@ -3217,7 +3244,7 @@ int generate_random_phase_field(
                     if (k2 > 0 && k2 < kk->kM2)
                     {
                         {
-                            const double amplitude = sqrt(
+                            const double amplitude = std::sqrt(
                                 output_field->cval(cindex, 0)*output_field->cval(cindex, 0)+
                                 output_field->cval(cindex, 1)*output_field->cval(cindex, 1));
                             if (amplitude > 0)
diff --git a/cpp/field.hpp b/cpp/field.hpp
index 1e558461bc5e943a4ecdd3a0a153ef20cee5b250..71f04444c5f16a541393b988afe4c3866915a473 100644
--- a/cpp/field.hpp
+++ b/cpp/field.hpp
@@ -28,12 +28,7 @@
 
 #define FIELD_HPP
 
-#include <hdf5.h>
-#include <unordered_map>
-#include <vector>
-#include <string>
 #include "kspace.hpp"
-#include "omputils.hpp"
 
 /** \class field
  *  \brief Holds field data, performs FFTs and HDF5 I/O operations.
@@ -341,35 +336,7 @@ class field
                               in_global_z - this->rlayout->starts[0]);
         }
 
-        int print_plan(const std::string preamble)
-        {
-            char *c2r_plan_information = fftw_interface<rnumber>::sprint(this->c2r_plan);
-            char *r2c_plan_information = fftw_interface<rnumber>::sprint(this->r2c_plan);
-            if (this->myrank == 0)
-            {
-                std::cout << preamble <<
-                             std::endl <<
-                             "----c2r plan is:\n" <<
-                             c2r_plan_information <<
-                             std::endl <<
-                             "----r2c plan is:\n" <<
-                             r2c_plan_information <<
-                             std::endl;
-            }
-            std::string err_message = (
-                    std::string("MPI rank ") +
-                    preamble +
-                    std::to_string(this->myrank) +
-                    std::string("\n----c2r plan is:\n") +
-                    std::string(c2r_plan_information) +
-                    std::string("\n----r2c plan is:\n") +
-                    std::string(r2c_plan_information) +
-                    std::string("\n"));
-            std::cerr << err_message;
-            free(c2r_plan_information);
-            free(r2c_plan_information);
-            return EXIT_SUCCESS;
-        }
+        int print_plan(const std::string preamble);
 };
 
 /** \brief Compute gradient of a field
diff --git a/cpp/field_binary_IO.cpp b/cpp/field_binary_IO.cpp
index 0d11599d5180c73ee590d156f1d82e0200ae0f48..1768c9e8cb467442e6b1c9cdc3c13f26f3c6d25b 100644
--- a/cpp/field_binary_IO.cpp
+++ b/cpp/field_binary_IO.cpp
@@ -24,12 +24,11 @@
 
 
 
-#include <vector>
-#include <array>
-#include "base.hpp"
 #include "scope_timer.hpp"
 #include "field_binary_IO.hpp"
 
+#include <array>
+
 template <typename rnumber, field_representation fr, field_components fc>
 field_binary_IO<rnumber, fr, fc>::field_binary_IO(
                 const hsize_t *SIZES,
diff --git a/cpp/field_binary_IO.hpp b/cpp/field_binary_IO.hpp
index c615f9e266fa7de5816c0cd98b5a830dc81d6e48..468e142bd6176f89c2e38b10a108e567325a789a 100644
--- a/cpp/field_binary_IO.hpp
+++ b/cpp/field_binary_IO.hpp
@@ -28,11 +28,6 @@
 
 #define FIELD_BINARY_IO_HPP
 
-#include <vector>
-#include <string>
-#include "base.hpp"
-#include "fftw_interface.hpp"
-#include "field_layout.hpp"
 #include "field.hpp"
 
 /* could this be a boolean somehow?*/
diff --git a/cpp/field_layout.cpp b/cpp/field_layout.cpp
index 0f19a6cd98a5dbe6d2be54d80f887d77a33f67ea..d918b722b0f4562d1a5fdc13994dd11e5af21f24 100644
--- a/cpp/field_layout.cpp
+++ b/cpp/field_layout.cpp
@@ -24,12 +24,10 @@
 
 
 
-#include <cassert>
 #include "field_layout.hpp"
 #include "scope_timer.hpp"
 
 
-
 template <field_components fc>
 field_layout<fc>::field_layout(
         const hsize_t *SIZES,
diff --git a/cpp/field_layout.hpp b/cpp/field_layout.hpp
index e7ba578cbeb5f7e990e2c21f95e0cb9cc1054985..deb80014c92c80f828a8da6a551a265180f830b9 100644
--- a/cpp/field_layout.hpp
+++ b/cpp/field_layout.hpp
@@ -28,9 +28,10 @@
 
 #define FIELD_LAYOUT_HPP
 
+#include "base.hpp"
+
 #include <vector>
 #include <hdf5.h>
-#include "base.hpp"
 
 enum field_components {ONE, THREE, THREExTHREE};
 
diff --git a/cpp/full_code/Gauss_field_test.cpp b/cpp/full_code/Gauss_field_test.cpp
index 6f7cdfd16877f3d207dd8f1a30c5058bf2a03dfb..3b40895a24d53f0e3cdba114fcbfd16afd10d0f2 100644
--- a/cpp/full_code/Gauss_field_test.cpp
+++ b/cpp/full_code/Gauss_field_test.cpp
@@ -24,16 +24,12 @@
 
 
 
-#include <string>
-#include <cmath>
-#include <random>
 #include "Gauss_field_test.hpp"
 #include "scope_timer.hpp"
 #include "hdf5_tools.hpp"
 
-
-
-
+#include <cmath>
+#include <random>
 
 template <typename rnumber>
 int Gauss_field_test<rnumber>::initialize(void)
diff --git a/cpp/full_code/Gauss_field_test.hpp b/cpp/full_code/Gauss_field_test.hpp
index 311e7a8c9f31f8172114d85242fdb519155f7856..3c74b2883c5927b5be68dee6d498ab59660de601 100644
--- a/cpp/full_code/Gauss_field_test.hpp
+++ b/cpp/full_code/Gauss_field_test.hpp
@@ -29,9 +29,6 @@
 
 
 
-#include <cstdlib>
-#include "base.hpp"
-#include "kspace.hpp"
 #include "field.hpp"
 #include "full_code/test.hpp"
 
diff --git a/cpp/full_code/NSE.cpp b/cpp/full_code/NSE.cpp
index 7a84475eabd194466f894ca740b4e89ff187686b..1ed1339fa0d089f8a05ace5693d23920e07b4c3f 100644
--- a/cpp/full_code/NSE.cpp
+++ b/cpp/full_code/NSE.cpp
@@ -23,15 +23,11 @@
 
 
 
-#include <string>
-#include <cmath>
 #include "NSE.hpp"
 #include "scope_timer.hpp"
 #include "fftw_tools.hpp"
-#include "hdf5_tools.hpp"
 #include "shared_array.hpp"
 
-
 template <typename rnumber>
 int NSE<rnumber>::initialize(void)
 {
@@ -139,9 +135,10 @@ int NSE<rnumber>::add_field_band(
                     const ptrdiff_t xindex,
                     const ptrdiff_t yindex,
                     const ptrdiff_t zindex){
-        const double knorm = sqrt(this->kk->kx[xindex]*this->kk->kx[xindex] +
-                                  this->kk->ky[yindex]*this->kk->ky[yindex] +
-                                  this->kk->kz[zindex]*this->kk->kz[zindex]);
+        const double knorm = std::sqrt(
+                this->kk->kx[xindex]*this->kk->kx[xindex] +
+                this->kk->ky[yindex]*this->kk->ky[yindex] +
+                this->kk->kz[zindex]*this->kk->kz[zindex]);
         if ((k0 <= knorm) &&
             (k1 >= knorm))
             for (int c=0; c<3; c++)
@@ -180,7 +177,7 @@ int NSE<rnumber>::add_forcing_term(
                         const ptrdiff_t zindex,
                         const double k2,
                         const int nxmodes){
-            const double knorm = sqrt(k2);
+            const double knorm = std::sqrt(k2);
             if ((k2 > 0) &&
                 (this->fk0 <= knorm) &&
                 (this->fk1 >= knorm))
@@ -558,7 +555,7 @@ int NSE<rnumber>::do_stats()
             stat_group,
             "velocity",
             this->iteration / niter_stat,
-            this->max_velocity_estimate/sqrt(3));
+            this->max_velocity_estimate/std::sqrt(3));
 
     compute_curl(this->kk,
                  this->velocity,
@@ -569,7 +566,7 @@ int NSE<rnumber>::do_stats()
             stat_group,
             "vorticity",
             this->iteration / niter_stat,
-            this->max_vorticity_estimate/sqrt(3));
+            this->max_vorticity_estimate/std::sqrt(3));
 
     if (this->myrank == 0)
         H5Gclose(stat_group);
diff --git a/cpp/full_code/NSE.hpp b/cpp/full_code/NSE.hpp
index 0ab0d3d6c6ea72ac4d92b41d648f14b861003bab..ef0813e74aa3a68980468d0d0a377015eda0edb2 100644
--- a/cpp/full_code/NSE.hpp
+++ b/cpp/full_code/NSE.hpp
@@ -29,8 +29,6 @@
 
 
 
-#include <cstdlib>
-#include "base.hpp"
 #include "vorticity_equation.hpp"
 #include "full_code/direct_numerical_simulation.hpp"
 
diff --git a/cpp/full_code/NSVE.cpp b/cpp/full_code/NSVE.cpp
index 866633265b4e551486f1f16c37e7c8dc812e0245..24fc890fef55d50d07d2b2e0439f276b86e54913 100644
--- a/cpp/full_code/NSVE.cpp
+++ b/cpp/full_code/NSVE.cpp
@@ -23,8 +23,6 @@
 
 
 
-#include <string>
-#include <cmath>
 #include "NSVE.hpp"
 #include "scope_timer.hpp"
 #include "fftw_tools.hpp"
@@ -182,7 +180,7 @@ int NSVE<rnumber>::do_stats()
             stat_group,
             "vorticity",
             fs->iteration / niter_stat,
-            max_vorticity_estimate/sqrt(3));
+            max_vorticity_estimate/std::sqrt(3));
 
     fs->compute_velocity(fs->cvorticity);
     *tmp_vec_field = fs->cvelocity->get_cdata();
@@ -191,7 +189,7 @@ int NSVE<rnumber>::do_stats()
             stat_group,
             "velocity",
             fs->iteration / niter_stat,
-            max_velocity_estimate/sqrt(3));
+            max_velocity_estimate/std::sqrt(3));
 
     double CFL_velocity = tmp_vec_field->compute_CFL_velocity();
     if (this->myrank == 0)
diff --git a/cpp/full_code/NSVE.hpp b/cpp/full_code/NSVE.hpp
index d7c7ab3ad9c1df8c9e8c9b210df03a0a2f760dca..338919d1f0fa8c77563d7df0c3e03cdbb4327da7 100644
--- a/cpp/full_code/NSVE.hpp
+++ b/cpp/full_code/NSVE.hpp
@@ -28,13 +28,9 @@
 #define NSVE_HPP
 
 
-#include "base.hpp"
 #include "vorticity_equation.hpp"
 #include "full_code/direct_numerical_simulation.hpp"
 
-#include <cstdlib>
-#include <string>
-
 /** \brief Navier-Stokes solver.
  *
  * Solves the vorticity equation for a Navier-Stokes fluid.
diff --git a/cpp/full_code/NSVE_Stokes_particles.cpp b/cpp/full_code/NSVE_Stokes_particles.cpp
index 8da3bcd0ac5dd8091d05b7c938856858dd8ffef8..19c193e3d4a80b67f6bea61b4af05a549e96be4a 100644
--- a/cpp/full_code/NSVE_Stokes_particles.cpp
+++ b/cpp/full_code/NSVE_Stokes_particles.cpp
@@ -24,14 +24,13 @@
 
 
 
-#include <string>
-#include <cmath>
 #include "full_code/NSVE_Stokes_particles.hpp"
 #include "scope_timer.hpp"
-#include "particles/particles_sampling.hpp"
 #include "particles/p2p/p2p_ghost_collisions.hpp"
 #include "particles/inner/particles_inner_computer_2nd_order.hpp"
 
+#include <cmath>
+
 template <typename rnumber>
 int NSVE_Stokes_particles<rnumber>::initialize(void)
 {
diff --git a/cpp/full_code/NSVE_Stokes_particles.hpp b/cpp/full_code/NSVE_Stokes_particles.hpp
index fb8497c99074e52c246ac83ba9a071e53d3023c6..f8ac865a76b84d25466e19b88b573ecf3f3356a7 100644
--- a/cpp/full_code/NSVE_Stokes_particles.hpp
+++ b/cpp/full_code/NSVE_Stokes_particles.hpp
@@ -29,12 +29,8 @@
 
 
 
-#include <cstdlib>
-#include "base.hpp"
-#include "vorticity_equation.hpp"
 #include "full_code/NSVE.hpp"
 #include "particles/particles_system_builder.hpp"
-#include "particles/particles_output_hdf5.hpp"
 #include "particles/particles_sampling.hpp"
 
 /** \brief Navier-Stokes solver that includes simple Lagrangian tracers.
diff --git a/cpp/full_code/NSVE_field_stats.cpp b/cpp/full_code/NSVE_field_stats.cpp
index 4ee109b88463436a9448a421dad31b5833af76b0..aa24666ae489ce2286e3d4a59345a1c41b3614ec 100644
--- a/cpp/full_code/NSVE_field_stats.cpp
+++ b/cpp/full_code/NSVE_field_stats.cpp
@@ -23,14 +23,11 @@
 
 
 
-#include <string>
-#include <cmath>
 #include "NSVE_field_stats.hpp"
 #include "fftw_tools.hpp"
 #include "scope_timer.hpp"
 #include "hdf5_tools.hpp"
 
-
 template <typename rnumber>
 int NSVE_field_stats<rnumber>::initialize(void)
 {
diff --git a/cpp/full_code/NSVE_field_stats.hpp b/cpp/full_code/NSVE_field_stats.hpp
index a04215389a7773f8ecbdf052d4a2a6fd4119b1b0..d260926152a5dc6e43330414f2417703c197dd39 100644
--- a/cpp/full_code/NSVE_field_stats.hpp
+++ b/cpp/full_code/NSVE_field_stats.hpp
@@ -27,12 +27,6 @@
 #ifndef NSVE_FIELD_STATS_HPP
 #define NSVE_FIELD_STATS_HPP
 
-#include <cstdlib>
-#include <sys/types.h>
-#include <sys/stat.h>
-#include <vector>
-#include "base.hpp"
-#include "field.hpp"
 #include "field_binary_IO.hpp"
 #include "full_code/postprocess.hpp"
 
diff --git a/cpp/full_code/NSVEcomplex_particles.cpp b/cpp/full_code/NSVEcomplex_particles.cpp
index 0c71cf723ca18e982a9b684a76b9f9c9f8834ae9..ef2980181689557e8d6c8eb53d23f50b5067338c 100644
--- a/cpp/full_code/NSVEcomplex_particles.cpp
+++ b/cpp/full_code/NSVEcomplex_particles.cpp
@@ -24,14 +24,12 @@
 
 
 
-#include <string>
-#include <cmath>
 #include "full_code/NSVEcomplex_particles.hpp"
-#include "scope_timer.hpp"
-#include "particles/particles_sampling.hpp"
 #include "particles/p2p/p2p_computer.hpp"
 #include "particles/inner/particles_inner_computer.hpp"
 
+#include <cmath>
+
 template <typename rnumber>
 int NSVEcomplex_particles<rnumber>::initialize(void)
 {
diff --git a/cpp/full_code/NSVEcomplex_particles.hpp b/cpp/full_code/NSVEcomplex_particles.hpp
index 3f6adaac45e79594e1e15a1e1fde737e4a1ea27f..75a77241e6623ea1579fd51945a4832aa8f879fe 100644
--- a/cpp/full_code/NSVEcomplex_particles.hpp
+++ b/cpp/full_code/NSVEcomplex_particles.hpp
@@ -29,12 +29,8 @@
 
 
 
-#include <cstdlib>
-#include "base.hpp"
-#include "vorticity_equation.hpp"
 #include "full_code/NSVE.hpp"
 #include "particles/particles_system_builder.hpp"
-#include "particles/particles_output_hdf5.hpp"
 #include "particles/particles_sampling.hpp"
 
 /** \brief Navier-Stokes solver that includes complex particles.
diff --git a/cpp/full_code/NSVEp_extra_sampling.cpp b/cpp/full_code/NSVEp_extra_sampling.cpp
index 0b96e259add4a1b6eca00eecd1ee430052ea8755..824c824455fa212e2e50347f415bd66343f4622d 100644
--- a/cpp/full_code/NSVEp_extra_sampling.cpp
+++ b/cpp/full_code/NSVEp_extra_sampling.cpp
@@ -25,8 +25,6 @@
 
 #include "full_code/NSVEp_extra_sampling.hpp"
 
-
-
 template <typename rnumber>
 int NSVEp_extra_sampling<rnumber>::initialize(void)
 {
diff --git a/cpp/full_code/NSVEp_extra_sampling.hpp b/cpp/full_code/NSVEp_extra_sampling.hpp
index a1f5a75a07e83c75291b7b699e4a518f36b18047..1e7b95a39e39a15f489c1b3b49f596a958cc539d 100644
--- a/cpp/full_code/NSVEp_extra_sampling.hpp
+++ b/cpp/full_code/NSVEp_extra_sampling.hpp
@@ -28,13 +28,7 @@
 
 
 
-#include <cstdlib>
-#include "base.hpp"
-#include "vorticity_equation.hpp"
 #include "full_code/NSVEparticles.hpp"
-#include "particles/particles_system_builder.hpp"
-#include "particles/particles_output_hdf5.hpp"
-#include "particles/particles_sampling.hpp"
 
 /** \brief Navier-Stokes solver with tracers that sample velocity gradient
  *  and pressure Hessian.
diff --git a/cpp/full_code/NSVEparticles.cpp b/cpp/full_code/NSVEparticles.cpp
index 237cd23ffba670a71d11a1c040417f2b5807c0d1..2fae3155002504575a833b4adae4476a6eb28af2 100644
--- a/cpp/full_code/NSVEparticles.cpp
+++ b/cpp/full_code/NSVEparticles.cpp
@@ -24,13 +24,12 @@
 
 
 
-#include <string>
-#include <cmath>
 #include "NSVEparticles.hpp"
-#include "scope_timer.hpp"
 #include "particles/interpolation/particle_set.hpp"
 #include "particles/particles_input_random.hpp"
 
+#include <cmath>
+
 template <typename rnumber>
 int NSVEparticles<rnumber>::initialize(void)
 {
diff --git a/cpp/full_code/NSVEparticles.hpp b/cpp/full_code/NSVEparticles.hpp
index dfa5598a911c3a6f7071312a8c22a43e9d0095e7..718d1d557a1716dde82ea38038be8fd0f36385bc 100644
--- a/cpp/full_code/NSVEparticles.hpp
+++ b/cpp/full_code/NSVEparticles.hpp
@@ -29,12 +29,8 @@
 
 
 
-#include <cstdlib>
-#include "base.hpp"
-#include "vorticity_equation.hpp"
 #include "full_code/NSVE.hpp"
 #include "particles/particles_system_builder.hpp"
-#include "particles/particles_output_hdf5.hpp"
 #include "particles/particles_sampling.hpp"
 
 /** \brief Navier-Stokes solver that includes simple Lagrangian tracers.
diff --git a/cpp/full_code/bandpass_stats.cpp b/cpp/full_code/bandpass_stats.cpp
index 7985a26af9be291915a0fcc088d4ccc084a6280f..cc765889eec4158dd3cc390afdb9da467496bb80 100644
--- a/cpp/full_code/bandpass_stats.cpp
+++ b/cpp/full_code/bandpass_stats.cpp
@@ -23,12 +23,11 @@
 
 
 
-#include <string>
-#include <cmath>
 #include "bandpass_stats.hpp"
 #include "scope_timer.hpp"
 #include "hdf5_tools.hpp"
 
+#include <cmath>
 
 template <typename rnumber>
 int bandpass_stats<rnumber>::initialize(void)
@@ -136,7 +135,7 @@ int bandpass_stats<rnumber>::work_on_current_iteration(void)
 
     /// make max vel estimate std::vector, required by compute_rspace_stats
     std::vector<double> max_vel_estimate;
-    max_vel_estimate.resize(4, max_velocity_estimate / sqrt(3));
+    max_vel_estimate.resize(4, max_velocity_estimate / std::sqrt(3));
     max_vel_estimate[3] = max_velocity_estimate;
 
     /// initialize `stat_group`.
diff --git a/cpp/full_code/bandpass_stats.hpp b/cpp/full_code/bandpass_stats.hpp
index 2210fd3efe5c80fd5410da34decdc00b13a931d3..d40ee891f5fe093f73b978a57de8b49ecdc300a1 100644
--- a/cpp/full_code/bandpass_stats.hpp
+++ b/cpp/full_code/bandpass_stats.hpp
@@ -27,13 +27,6 @@
 #ifndef BANDPASS_STATS_HPP
 #define BANDPASS_STATS_HPP
 
-#include <cstdlib>
-#include <sys/types.h>
-#include <sys/stat.h>
-#include <vector>
-#include "base.hpp"
-#include "field.hpp"
-#include "field_binary_IO.hpp"
 #include "full_code/NSVE_field_stats.hpp"
 
 /** \brief A class for running bandpass statistics.
diff --git a/cpp/full_code/code_base.hpp b/cpp/full_code/code_base.hpp
index adf21f6b2dde29f2655ec6af573a9d8f30c2d105..2886447dc7b3a5d8658d1a8f08ec8d994318cc59 100644
--- a/cpp/full_code/code_base.hpp
+++ b/cpp/full_code/code_base.hpp
@@ -27,10 +27,11 @@
 #ifndef CODE_BASE_HPP
 #define CODE_BASE_HPP
 
-#include <cstdlib>
+#include "base.hpp"
+
 #include <sys/types.h>
 #include <sys/stat.h>
-#include "base.hpp"
+#include <iostream>
 
 /** \class code_base
  *  \brief Defines basic timer and method to check stopping condition.
diff --git a/cpp/full_code/dealias_test.cpp b/cpp/full_code/dealias_test.cpp
index 56980b54a9bbb8c7e35574e7c26908cd948a831f..b2f6b6f77d6d4d0043a63ae7f4d0afb91e0bbe68 100644
--- a/cpp/full_code/dealias_test.cpp
+++ b/cpp/full_code/dealias_test.cpp
@@ -24,16 +24,12 @@
 
 
 
-#include <string>
-#include <cmath>
-#include <random>
 #include "dealias_test.hpp"
 #include "scope_timer.hpp"
 #include "hdf5_tools.hpp"
 
-
-
-
+#include <cmath>
+#include <random>
 
 template <typename rnumber>
 int dealias_test<rnumber>::initialize(void)
diff --git a/cpp/full_code/dealias_test.hpp b/cpp/full_code/dealias_test.hpp
index ba765c45086760910cd8bbc7270b55189e58b48e..1bdc23a6bf41457fdd3957741c4910aa3850e999 100644
--- a/cpp/full_code/dealias_test.hpp
+++ b/cpp/full_code/dealias_test.hpp
@@ -29,9 +29,6 @@
 
 
 
-#include <cstdlib>
-#include "base.hpp"
-#include "kspace.hpp"
 #include "field.hpp"
 #include "full_code/test.hpp"
 
diff --git a/cpp/full_code/direct_numerical_simulation.cpp b/cpp/full_code/direct_numerical_simulation.cpp
index e7c8a226d419a9ce74b603a3d41cf2d73f5e0d52..050c8f42ea4fe8a2a02004512f04fda934f23750 100644
--- a/cpp/full_code/direct_numerical_simulation.cpp
+++ b/cpp/full_code/direct_numerical_simulation.cpp
@@ -23,12 +23,8 @@
 
 
 
-#include <cstdlib>
-#include <sys/types.h>
-#include <sys/stat.h>
 #include "direct_numerical_simulation.hpp"
 #include "scope_timer.hpp"
-#include "hdf5_tools.hpp"
 
 int direct_numerical_simulation::grow_file_datasets()
 {
diff --git a/cpp/full_code/direct_numerical_simulation.hpp b/cpp/full_code/direct_numerical_simulation.hpp
index 3c71e45324422cb743e5cf7fb26ecbabda70290b..a5132099e7b3178a8b3f61f220189fb46d1ce3a1 100644
--- a/cpp/full_code/direct_numerical_simulation.hpp
+++ b/cpp/full_code/direct_numerical_simulation.hpp
@@ -27,10 +27,6 @@
 #ifndef DIRECT_NUMERICAL_SIMULATION_HPP
 #define DIRECT_NUMERICAL_SIMULATION_HPP
 
-#include <cstdlib>
-#include <sys/types.h>
-#include <sys/stat.h>
-#include "base.hpp"
 #include "hdf5_tools.hpp"
 #include "full_code/code_base.hpp"
 
diff --git a/cpp/full_code/field_output_test.cpp b/cpp/full_code/field_output_test.cpp
index ec90948dab64280ad99cee472a297ba80fd1d97e..fdea1f4cb33fe4c565ad4d4b4dc8b49ac0e06ee8 100644
--- a/cpp/full_code/field_output_test.cpp
+++ b/cpp/full_code/field_output_test.cpp
@@ -28,7 +28,6 @@
 #include "hdf5_tools.hpp"
 #include "shared_array.hpp"
 
-#include <string>
 #include <cmath>
 #include <random>
 
diff --git a/cpp/full_code/field_output_test.hpp b/cpp/full_code/field_output_test.hpp
index 75ecafa391395c228d634844cd433ece99a27ed8..d8e0ecab7c7f029a6598abdcbb2e05e50a912923 100644
--- a/cpp/full_code/field_output_test.hpp
+++ b/cpp/full_code/field_output_test.hpp
@@ -29,9 +29,6 @@
 
 
 
-#include <cstdlib>
-#include "base.hpp"
-#include "kspace.hpp"
 #include "field.hpp"
 #include "full_code/test.hpp"
 
diff --git a/cpp/full_code/field_single_to_double.cpp b/cpp/full_code/field_single_to_double.cpp
index 9f2082e63613e7c2d64943c01d8c782f57e4a1cf..df9bacd5acde8cb4fb6aae91a370a5a493d8d75c 100644
--- a/cpp/full_code/field_single_to_double.cpp
+++ b/cpp/full_code/field_single_to_double.cpp
@@ -23,12 +23,11 @@
 
 
 
-#include <string>
-#include <cmath>
 #include "field_single_to_double.hpp"
 #include "scope_timer.hpp"
 #include "hdf5_tools.hpp"
 
+#include <cmath>
 
 template <typename rnumber>
 int field_single_to_double<rnumber>::initialize(void)
diff --git a/cpp/full_code/field_single_to_double.hpp b/cpp/full_code/field_single_to_double.hpp
index de98e63ed5fda3fa9e6ffadfbb076c38436d414f..a044e87b6807047c6db758455d5a52d6d62c8fb5 100644
--- a/cpp/full_code/field_single_to_double.hpp
+++ b/cpp/full_code/field_single_to_double.hpp
@@ -27,13 +27,6 @@
 #ifndef FIELD_SINGLE_TO_DOUBLE_HPP
 #define FIELD_SINGLE_TO_DOUBLE_HPP
 
-#include <cstdlib>
-#include <sys/types.h>
-#include <sys/stat.h>
-#include <vector>
-#include "base.hpp"
-#include "field.hpp"
-#include "field_binary_IO.hpp"
 #include "full_code/NSVE_field_stats.hpp"
 
 template <typename rnumber>
diff --git a/cpp/full_code/field_test.cpp b/cpp/full_code/field_test.cpp
index 88ebea24946aada4cb3030b69412f579b8461188..0ae1609438807e3fe6002211b994c3533de93e33 100644
--- a/cpp/full_code/field_test.cpp
+++ b/cpp/full_code/field_test.cpp
@@ -23,13 +23,12 @@
 
 
 
-#include <string>
-#include <cmath>
-#include <random>
 #include "field_test.hpp"
 #include "scope_timer.hpp"
 #include "hdf5_tools.hpp"
 
+#include <cmath>
+#include <random>
 
 template <typename rnumber>
 int field_test<rnumber>::initialize(void)
diff --git a/cpp/full_code/field_test.hpp b/cpp/full_code/field_test.hpp
index 5a57c1d2e6f5b0f92ddeb8a01f238521cf3e5bfe..85fa2ba7dfa377ab72632100a53d0381d523d1aa 100644
--- a/cpp/full_code/field_test.hpp
+++ b/cpp/full_code/field_test.hpp
@@ -29,9 +29,6 @@
 
 
 
-#include <cstdlib>
-#include "base.hpp"
-#include "kspace.hpp"
 #include "field.hpp"
 #include "full_code/test.hpp"
 
diff --git a/cpp/full_code/filter_test.cpp b/cpp/full_code/filter_test.cpp
index e70cc61c7ce8d40b606262a451e99f3829958832..860a1d16d7dc20326c622276f5288dd75b36b7b9 100644
--- a/cpp/full_code/filter_test.cpp
+++ b/cpp/full_code/filter_test.cpp
@@ -23,12 +23,11 @@
 
 
 
-#include <string>
-#include <cmath>
 #include "filter_test.hpp"
 #include "scope_timer.hpp"
 #include "hdf5_tools.hpp"
 
+#include <cmath>
 
 template <typename rnumber>
 int filter_test<rnumber>::initialize(void)
diff --git a/cpp/full_code/filter_test.hpp b/cpp/full_code/filter_test.hpp
index c7df600f7b638c34617be03e68691f2305f5087a..ebe5fa80b6a25f9657d3fbf8c3b478eb9e4136c2 100644
--- a/cpp/full_code/filter_test.hpp
+++ b/cpp/full_code/filter_test.hpp
@@ -29,9 +29,6 @@
 
 
 
-#include <cstdlib>
-#include "base.hpp"
-#include "kspace.hpp"
 #include "field.hpp"
 #include "full_code/test.hpp"
 
diff --git a/cpp/full_code/get_3D_correlations.cpp b/cpp/full_code/get_3D_correlations.cpp
index 6970e7af4a10e2f97d0a9fca542a0ee65989d800..d3cff01ddbd52c7ffdd0dbda19baf4e20b8312db 100644
--- a/cpp/full_code/get_3D_correlations.cpp
+++ b/cpp/full_code/get_3D_correlations.cpp
@@ -23,12 +23,12 @@
 
 
 
-#include <string>
-#include <cmath>
 #include "get_3D_correlations.hpp"
 #include "scope_timer.hpp"
 #include "hdf5_tools.hpp"
 
+#include <cmath>
+
 template <typename rnumber>
 int get_3D_correlations<rnumber>::initialize(void)
 {
diff --git a/cpp/full_code/get_3D_correlations.hpp b/cpp/full_code/get_3D_correlations.hpp
index e328d83cea5921de74a83474baf363a37540c2f4..229310963bfc2ace592623c0af42386e9efce83a 100644
--- a/cpp/full_code/get_3D_correlations.hpp
+++ b/cpp/full_code/get_3D_correlations.hpp
@@ -27,13 +27,6 @@
 #ifndef GET_3D_CORRELATIONS_HPP
 #define GET_3D_CORRELATIONS_HPP
 
-#include <cstdlib>
-#include <sys/types.h>
-#include <sys/stat.h>
-#include <vector>
-#include "base.hpp"
-#include "field.hpp"
-#include "field_binary_IO.hpp"
 #include "full_code/NSVE_field_stats.hpp"
 
 template <typename rnumber>
diff --git a/cpp/full_code/get_rfields.cpp b/cpp/full_code/get_rfields.cpp
index 85864ee4d8f4fda468539a94cd43512a6aea5573..30184913ffad591d0f0c6f60ccd3eec3d53387a6 100644
--- a/cpp/full_code/get_rfields.cpp
+++ b/cpp/full_code/get_rfields.cpp
@@ -23,12 +23,11 @@
 
 
 
-#include <string>
-#include <cmath>
 #include "get_rfields.hpp"
 #include "scope_timer.hpp"
 #include "hdf5_tools.hpp"
 
+#include <cmath>
 
 template <typename rnumber>
 int get_rfields<rnumber>::initialize(void)
diff --git a/cpp/full_code/get_rfields.hpp b/cpp/full_code/get_rfields.hpp
index 2d266d390322fbf3c2ecb631c95e074aef298491..25fe453f3e69331c72541f5cd53f40ad2f732096 100644
--- a/cpp/full_code/get_rfields.hpp
+++ b/cpp/full_code/get_rfields.hpp
@@ -27,13 +27,6 @@
 #ifndef GET_RFIELDS_HPP
 #define GET_RFIELDS_HPP
 
-#include <cstdlib>
-#include <sys/types.h>
-#include <sys/stat.h>
-#include <vector>
-#include "base.hpp"
-#include "field.hpp"
-#include "field_binary_IO.hpp"
 #include "full_code/NSVE_field_stats.hpp"
 
 template <typename rnumber>
diff --git a/cpp/full_code/get_velocity.cpp b/cpp/full_code/get_velocity.cpp
index a29cc660c72fd6e07cfcd41ad9f03482be7f1c1b..37863153c32c792c6abbd4f12d97f7354b98ef08 100644
--- a/cpp/full_code/get_velocity.cpp
+++ b/cpp/full_code/get_velocity.cpp
@@ -23,12 +23,11 @@
 
 
 
-#include <string>
-#include <cmath>
 #include "get_velocity.hpp"
 #include "scope_timer.hpp"
 #include "hdf5_tools.hpp"
 
+#include <cmath>
 
 template <typename rnumber>
 int get_velocity<rnumber>::initialize(void)
diff --git a/cpp/full_code/get_velocity.hpp b/cpp/full_code/get_velocity.hpp
index dababf2c16c290b2e6aa60d76f151026bc49c167..48e754d13adbea5f3adee3cda69d9d12260edd89 100644
--- a/cpp/full_code/get_velocity.hpp
+++ b/cpp/full_code/get_velocity.hpp
@@ -27,13 +27,6 @@
 #ifndef GET_VELOCITY_HPP
 #define GET_VELOCITY_HPP
 
-#include <cstdlib>
-#include <sys/types.h>
-#include <sys/stat.h>
-#include <vector>
-#include "base.hpp"
-#include "field.hpp"
-#include "field_binary_IO.hpp"
 #include "full_code/NSVE_field_stats.hpp"
 
 template <typename rnumber>
diff --git a/cpp/full_code/joint_acc_vel_stats.cpp b/cpp/full_code/joint_acc_vel_stats.cpp
index a685c1f82f01ba53fa849e5bbe88792f977bcb62..619b316478add62cb404843764394468e653cb94 100644
--- a/cpp/full_code/joint_acc_vel_stats.cpp
+++ b/cpp/full_code/joint_acc_vel_stats.cpp
@@ -23,12 +23,11 @@
 
 
 
-#include <string>
-#include <cmath>
 #include "joint_acc_vel_stats.hpp"
 #include "scope_timer.hpp"
 #include "hdf5_tools.hpp"
 
+#include <cmath>
 
 template <typename rnumber>
 int joint_acc_vel_stats<rnumber>::initialize(void)
@@ -149,8 +148,8 @@ int joint_acc_vel_stats<rnumber>::work_on_current_iteration(void)
     std::vector<double> max_acc_estimate;
     std::vector<double> max_vel_estimate;
 
-    max_acc_estimate.resize(4, max_acceleration_estimate / sqrt(3));
-    max_vel_estimate.resize(4, max_velocity_estimate / sqrt(3));
+    max_acc_estimate.resize(4, max_acceleration_estimate / std::sqrt(3));
+    max_vel_estimate.resize(4, max_velocity_estimate / std::sqrt(3));
     max_acc_estimate[3] = max_acceleration_estimate;
     max_vel_estimate[3] = max_velocity_estimate;
 
@@ -169,13 +168,13 @@ int joint_acc_vel_stats<rnumber>::work_on_current_iteration(void)
             stat_group,
             "acceleration",
             this->iteration / this->niter_out,
-            max_acceleration_estimate / sqrt(3));
+            max_acceleration_estimate / std::sqrt(3));
     vel->compute_stats(
             kk,
             stat_group,
             "velocity",
             this->iteration / this->niter_out,
-            max_velocity_estimate / sqrt(3));
+            max_velocity_estimate / std::sqrt(3));
 
     /// close stat group
     if (this->myrank == 0)
diff --git a/cpp/full_code/joint_acc_vel_stats.hpp b/cpp/full_code/joint_acc_vel_stats.hpp
index 814ddfd1b0b0c24af6e41c2397c862c35d4c17b1..cf293d38d242bd125fb28aa6b29929306de3f9e7 100644
--- a/cpp/full_code/joint_acc_vel_stats.hpp
+++ b/cpp/full_code/joint_acc_vel_stats.hpp
@@ -27,14 +27,6 @@
 #ifndef JOINT_ACC_VEL_STATS_HPP
 #define JOINT_ACC_VEL_STATS_HPP
 
-#include <cstdlib>
-#include <sys/types.h>
-#include <sys/stat.h>
-#include <vector>
-#include "base.hpp"
-#include "field.hpp"
-#include "field_binary_IO.hpp"
-#include "vorticity_equation.hpp"
 #include "full_code/NSVE_field_stats.hpp"
 
 template <typename rnumber>
diff --git a/cpp/full_code/kraichnan_field.cpp b/cpp/full_code/kraichnan_field.cpp
index bb63f5def0bef1b470801aff19bbd2d916e16235..213b730ff37a10efb3141d3b30788e6b138ea9bd 100644
--- a/cpp/full_code/kraichnan_field.cpp
+++ b/cpp/full_code/kraichnan_field.cpp
@@ -23,8 +23,6 @@
 
 
 
-#include <string>
-#include <cmath>
 #include "kraichnan_field.hpp"
 #include "scope_timer.hpp"
 #include "fftw_tools.hpp"
@@ -155,7 +153,7 @@ template <typename rnumber>
 int kraichnan_field<rnumber>::step(void)
 {
     TIMEZONE("kraichnan_file::step");
-    this->ps->completeLoop(sqrt(this->dt));
+    this->ps->completeLoop(std::sqrt(this->dt));
     this->generate_random_velocity();
     this->iteration++;
     return EXIT_SUCCESS;
@@ -234,7 +232,7 @@ int kraichnan_field<rnumber>::do_stats()
             stat_group,
             "velocity",
             this->iteration / this->niter_stat,
-            this->max_velocity_estimate/sqrt(3));
+            this->max_velocity_estimate/std::sqrt(3));
         if (this->myrank == 0)
             H5Gclose(stat_group);
         delete tvelocity;
diff --git a/cpp/full_code/kraichnan_field.hpp b/cpp/full_code/kraichnan_field.hpp
index 020cd0cdc7b2c5b8e35114bc7d78106e9a8c1fb2..f64de4e2982fec2aeefaa9ae35b3815de29e1838 100644
--- a/cpp/full_code/kraichnan_field.hpp
+++ b/cpp/full_code/kraichnan_field.hpp
@@ -29,13 +29,8 @@
 
 
 
-#include <cstdlib>
-#include "base.hpp"
-#include "kspace.hpp"
-#include "field.hpp"
 #include "full_code/direct_numerical_simulation.hpp"
 #include "particles/particles_system_builder.hpp"
-#include "particles/particles_output_hdf5.hpp"
 #include "particles/particles_sampling.hpp"
 
 template <typename rnumber>
diff --git a/cpp/full_code/main_code.hpp b/cpp/full_code/main_code.hpp
index 6715f26eadf1089fea53147ccf7afae38abd5286..8b3336460020a0bb985cbf7c8ff0408db7c91459 100644
--- a/cpp/full_code/main_code.hpp
+++ b/cpp/full_code/main_code.hpp
@@ -28,14 +28,12 @@
 #define MAIN_CODE_HPP
 
 
-
-#include <cfenv>
-#include <string>
-#include <iostream>
-#include "base.hpp"
 #include "field.hpp"
 #include "scope_timer.hpp"
 
+
+#include <cfenv>
+
 int myrank, nprocs;
 
 #ifdef PINCHECK_FOUND
diff --git a/cpp/full_code/native_binary_to_hdf5.cpp b/cpp/full_code/native_binary_to_hdf5.cpp
index ac0973d61144e0e7bdd5bb08533ece0779d7d757..5451e9e2abb1a8f439b60bad4237723f57e6c367 100644
--- a/cpp/full_code/native_binary_to_hdf5.cpp
+++ b/cpp/full_code/native_binary_to_hdf5.cpp
@@ -23,7 +23,6 @@
 
 
 
-#include <string>
 #include <cmath>
 #include "native_binary_to_hdf5.hpp"
 #include "scope_timer.hpp"
diff --git a/cpp/full_code/native_binary_to_hdf5.hpp b/cpp/full_code/native_binary_to_hdf5.hpp
index 3baaf609bfa4bba093dd097dce3be3d3aa590a4e..2884fbe96959989798555fea1ad677516cdd6bd1 100644
--- a/cpp/full_code/native_binary_to_hdf5.hpp
+++ b/cpp/full_code/native_binary_to_hdf5.hpp
@@ -27,12 +27,6 @@
 #ifndef NATIVE_BINARY_TO_HDF5_HPP
 #define NATIVE_BINARY_TO_HDF5_HPP
 
-#include <cstdlib>
-#include <sys/types.h>
-#include <sys/stat.h>
-#include <vector>
-#include "base.hpp"
-#include "field.hpp"
 #include "field_binary_IO.hpp"
 #include "full_code/postprocess.hpp"
 
diff --git a/cpp/full_code/ornstein_uhlenbeck_process.cpp b/cpp/full_code/ornstein_uhlenbeck_process.cpp
index 653c2540f01ee18661346b04945c6cade09450d2..23b2172708c23b68240f46f31262c9e9e5bd6c0e 100644
--- a/cpp/full_code/ornstein_uhlenbeck_process.cpp
+++ b/cpp/full_code/ornstein_uhlenbeck_process.cpp
@@ -96,7 +96,7 @@ void ornstein_uhlenbeck_process<rnumber,be>::step(
                     double k2){
         if (k2 <= this->ou_kmax_squ && k2 >= this->ou_kmin_squ)
         {
-            double g = this->gamma(sqrt(k2));
+            double g = this->gamma(std::sqrt(k2));
             int tr = omp_get_thread_num();
             double random[6] =
             {
@@ -106,7 +106,7 @@ void ornstein_uhlenbeck_process<rnumber,be>::step(
             for (int cc=0; cc<3; cc++) for (int i=0; i<2; i++)
                 this->ou_field->cval(cindex,cc,i) =
                     (1-dt*g) * this->ou_field->cval(cindex,cc,i) +
-                    sqrt(dt) * (
+                    std::sqrt(dt) * (
                         this->B->cval(cindex,0,cc,0) * random[0+3*i] +
                         this->B->cval(cindex,1,cc,0) * random[1+3*i] +
                         this->B->cval(cindex,2,cc,0) * random[2+3*i] );
@@ -135,7 +135,7 @@ void ornstein_uhlenbeck_process<rnumber,be>::initialize_B()
                           pow(this->kk->kz[zindex],2);
             if (ksqu <= this->ou_kmax_squ && ksqu >= this->ou_kmin_squ)
             {
-                double kabs = sqrt(ksqu);
+                double kabs = std::sqrt(ksqu);
                 double sigma;
                 if (ksqu == 0 || ksqu == this->ou_kmax_squ)
                 {
@@ -144,7 +144,7 @@ void ornstein_uhlenbeck_process<rnumber,be>::initialize_B()
                 else
                 {
                     sigma =
-                        sqrt(4*this->gamma(kabs)*this->energy(kabs)
+                        std::sqrt(4*this->gamma(kabs)*this->energy(kabs)
                         /(4*M_PI*ksqu));
                 }
 
@@ -188,7 +188,7 @@ void ornstein_uhlenbeck_process<rnumber,be>::initialize_B()
 template <class rnumber, field_backend be>
 void ornstein_uhlenbeck_process<rnumber, be>::let_converge(void)
 {
-  double ou_kmin = sqrt(this->ou_kmin_squ);
+  double ou_kmin = std::sqrt(this->ou_kmin_squ);
   double tau = 1.0/this->gamma(ou_kmin);
   double dt = tau/1000.;
 
diff --git a/cpp/full_code/ornstein_uhlenbeck_process.hpp b/cpp/full_code/ornstein_uhlenbeck_process.hpp
index 5a1288c929afac4cad4a3a6f36ff5b2a0436b1ca..3da3f0b92afb80fd27e24dafff5db73dd2d4c72d 100644
--- a/cpp/full_code/ornstein_uhlenbeck_process.hpp
+++ b/cpp/full_code/ornstein_uhlenbeck_process.hpp
@@ -1,11 +1,12 @@
 #ifndef ORNSTEIN_UHLENBECK_PROCESS_HPP
 #define ORNSTEIN_UHLENBECK_PROCESS_HPP
 
+#include "field.hpp"
+
 #include <sys/stat.h>
-#include <stdio.h>
-#include <stdlib.h>
+#include <cstdio>
+#include <cstdlib>
 #include <iostream>
-#include "field.hpp"
 #include <random>
 
 extern int myrank, nprocs;
diff --git a/cpp/full_code/ou_vorticity_equation.hpp b/cpp/full_code/ou_vorticity_equation.hpp
index 0f25356787a2c725af4e462b44305ec16a22521a..efe25122c681ece51bd4b49ec5a3bddcd6977207 100644
--- a/cpp/full_code/ou_vorticity_equation.hpp
+++ b/cpp/full_code/ou_vorticity_equation.hpp
@@ -1,7 +1,6 @@
 #ifndef OU_VORTICITY_EQUATION_HPP
 #define OU_VORTICITY_EQUATION_HPP
 
-#include <cstdlib>
 #include "vorticity_equation.hpp"
 #include "ornstein_uhlenbeck_process.hpp"
 
diff --git a/cpp/full_code/phase_shift_test.cpp b/cpp/full_code/phase_shift_test.cpp
index 7e6c2ad402ca3becf99c497eceefd511ae61f113..02cafeb21c09671df968eb6a269fe765c91217fb 100644
--- a/cpp/full_code/phase_shift_test.cpp
+++ b/cpp/full_code/phase_shift_test.cpp
@@ -23,11 +23,10 @@
 
 
 
-#include <string>
-#include <cmath>
-#include <random>
 #include "phase_shift_test.hpp"
 #include "scope_timer.hpp"
+#include <cmath>
+#include <random>
 
 
 template <typename rnumber>
diff --git a/cpp/full_code/phase_shift_test.hpp b/cpp/full_code/phase_shift_test.hpp
index e98f76cd568239a88f15da18c4fa72bc255fae9b..577b413e54d44e05ad59cdc3693d6a04200992d1 100644
--- a/cpp/full_code/phase_shift_test.hpp
+++ b/cpp/full_code/phase_shift_test.hpp
@@ -28,9 +28,6 @@
 
 
 
-#include <cstdlib>
-#include "base.hpp"
-#include "kspace.hpp"
 #include "field.hpp"
 #include "full_code/test.hpp"
 
diff --git a/cpp/full_code/postprocess.cpp b/cpp/full_code/postprocess.cpp
index 2f2feaae197efad44038858e90b732d708625cf5..b280034d0c6dd5338012e24673c2c6c4e8e8c3fb 100644
--- a/cpp/full_code/postprocess.cpp
+++ b/cpp/full_code/postprocess.cpp
@@ -22,10 +22,6 @@
 ******************************************************************************/
 
 
-
-#include <cstdlib>
-#include <sys/types.h>
-#include <sys/stat.h>
 #include "scope_timer.hpp"
 #include "hdf5_tools.hpp"
 #include "full_code/postprocess.hpp"
diff --git a/cpp/full_code/postprocess.hpp b/cpp/full_code/postprocess.hpp
index 52826cf7acf436823ad200e28ca29534bcf5e355..df41b3928b9948f9c098f2ceb0e0a477e596f542 100644
--- a/cpp/full_code/postprocess.hpp
+++ b/cpp/full_code/postprocess.hpp
@@ -29,12 +29,6 @@
 
 #include "full_code/code_base.hpp"
 #include "vorticity_equation.hpp"
-#include "base.hpp"
-
-#include <cstdlib>
-#include <sys/types.h>
-#include <sys/stat.h>
-#include <vector>
 
 class postprocess: public code_base
 {
diff --git a/cpp/full_code/pressure_stats.cpp b/cpp/full_code/pressure_stats.cpp
index c652f346d939888ba57cd7063cffffabdeb83c11..99eb9777929e3bba527871d6fe98a7dea5153d0a 100644
--- a/cpp/full_code/pressure_stats.cpp
+++ b/cpp/full_code/pressure_stats.cpp
@@ -23,12 +23,11 @@
 
 
 
-#include <string>
-#include <cmath>
 #include "pressure_stats.hpp"
 #include "scope_timer.hpp"
 #include "hdf5_tools.hpp"
 
+#include <cmath>
 
 template <typename rnumber>
 int pressure_stats<rnumber>::initialize(void)
diff --git a/cpp/full_code/pressure_stats.hpp b/cpp/full_code/pressure_stats.hpp
index f20ba08ee767fcd44bf6017670f3a3a41d7ba96e..c60b56ffed70894c85be988664d0e6b0f96c66de 100644
--- a/cpp/full_code/pressure_stats.hpp
+++ b/cpp/full_code/pressure_stats.hpp
@@ -27,14 +27,6 @@
 #ifndef PRESSURE_STATS_HPP
 #define PRESSURE_STATS_HPP
 
-#include <cstdlib>
-#include <sys/types.h>
-#include <sys/stat.h>
-#include <vector>
-#include "base.hpp"
-#include "field.hpp"
-#include "field_binary_IO.hpp"
-#include "vorticity_equation.hpp"
 #include "full_code/NSVE_field_stats.hpp"
 
 template <typename rnumber>
diff --git a/cpp/full_code/resize.cpp b/cpp/full_code/resize.cpp
index 12f125e026e9faf864b5d510d59371390d5874e9..8c1553b34a7a958b93c6f7aa509f7aebb9910ebf 100644
--- a/cpp/full_code/resize.cpp
+++ b/cpp/full_code/resize.cpp
@@ -23,12 +23,11 @@
 
 
 
-#include <string>
-#include <cmath>
 #include "resize.hpp"
 #include "scope_timer.hpp"
 #include "hdf5_tools.hpp"
 
+#include <cmath>
 
 template <typename rnumber>
 int resize<rnumber>::initialize(void)
diff --git a/cpp/full_code/resize.hpp b/cpp/full_code/resize.hpp
index 9123dd9c85dbc767222b03ef9c9291b3ff7d9896..b10b5101869327e3e8dc1e71fc6d90e33123f085 100644
--- a/cpp/full_code/resize.hpp
+++ b/cpp/full_code/resize.hpp
@@ -27,13 +27,6 @@
 #ifndef RESIZE_HPP
 #define RESIZE_HPP
 
-#include <cstdlib>
-#include <sys/types.h>
-#include <sys/stat.h>
-#include <vector>
-#include "base.hpp"
-#include "field.hpp"
-#include "field_binary_IO.hpp"
 #include "full_code/NSVE_field_stats.hpp"
 
 template <typename rnumber>
diff --git a/cpp/full_code/shared_array_test.hpp b/cpp/full_code/shared_array_test.hpp
index b81098381072520210bac5120dccff1600ece9db..5a4b95eba72fe186411ea1e4f368bf2aa7501abc 100644
--- a/cpp/full_code/shared_array_test.hpp
+++ b/cpp/full_code/shared_array_test.hpp
@@ -28,7 +28,6 @@
 
 
 
-#include "base.hpp"
 #include "test.hpp"
 #include "shared_array.hpp"
 
diff --git a/cpp/full_code/static_field.cpp b/cpp/full_code/static_field.cpp
index 1fa52b646deb48f53b2066440221e10b01c3ab55..faccf7f3434424d3428978869ec0343fb08d24b6 100644
--- a/cpp/full_code/static_field.cpp
+++ b/cpp/full_code/static_field.cpp
@@ -23,8 +23,6 @@
 
 
 
-#include <string>
-#include <cmath>
 #include "static_field.hpp"
 #include "scope_timer.hpp"
 #include "fftw_tools.hpp"
diff --git a/cpp/full_code/static_field.hpp b/cpp/full_code/static_field.hpp
index 90b563d1cd777768a63591e78ef121c9a594fc27..ccc1d4c45bef7bb8566a718ae9c6676e1f3287f5 100644
--- a/cpp/full_code/static_field.hpp
+++ b/cpp/full_code/static_field.hpp
@@ -29,13 +29,8 @@
 
 
 
-#include <cstdlib>
-#include "base.hpp"
-#include "kspace.hpp"
-#include "field.hpp"
 #include "full_code/direct_numerical_simulation.hpp"
 #include "particles/particles_system_builder.hpp"
-#include "particles/particles_output_hdf5.hpp"
 #include "particles/particles_sampling.hpp"
 
 template <typename rnumber>
diff --git a/cpp/full_code/symmetrize_test.cpp b/cpp/full_code/symmetrize_test.cpp
index 904b57ae2e32004c97f0ca1046e450822f853809..943014ee5f5b9753f9031b36ce17a6bc006bbdb4 100644
--- a/cpp/full_code/symmetrize_test.cpp
+++ b/cpp/full_code/symmetrize_test.cpp
@@ -23,14 +23,11 @@
 
 
 
-#include <string>
-#include <cmath>
 #include "symmetrize_test.hpp"
 #include "fftw_tools.hpp"
 #include "scope_timer.hpp"
 #include "hdf5_tools.hpp"
 
-
 template <typename rnumber>
 int symmetrize_test<rnumber>::initialize(void)
 {
@@ -140,8 +137,14 @@ int symmetrize_test<rnumber>::do_work(void)
     // temporary variables for field comparison
     double max_diff;
     ptrdiff_t ix, iy, iz;
+    variable_used_only_in_assert(ix);
+    variable_used_only_in_assert(iy);
+    variable_used_only_in_assert(iz);
     double k_at_max_diff;
+    variable_used_only_in_assert(k_at_max_diff);
     double a0, a1;
+    variable_used_only_in_assert(a0);
+    variable_used_only_in_assert(a1);
 
     // now compare the two fields
     max_diff = 0;
@@ -158,8 +161,9 @@ int symmetrize_test<rnumber>::do_work(void)
             double diff_im0 = test_field0->cval(cindex, 0, 1) - test_field1->cval(cindex, 0, 1);
             double diff_im1 = test_field0->cval(cindex, 1, 1) - test_field1->cval(cindex, 1, 1);
             double diff_im2 = test_field0->cval(cindex, 2, 1) - test_field1->cval(cindex, 2, 1);
-            double diff = sqrt(diff_re0*diff_re0 + diff_re1*diff_re1 + diff_re2*diff_re2 +
-                               diff_im0*diff_im0 + diff_im1*diff_im1 + diff_im2*diff_im2);
+            double diff = std::sqrt(
+                    diff_re0*diff_re0 + diff_re1*diff_re1 + diff_re2*diff_re2 +
+                    diff_im0*diff_im0 + diff_im1*diff_im1 + diff_im2*diff_im2);
             double amplitude0 = (test_field0->cval(cindex, 0, 0)*test_field0->cval(cindex, 0, 0) +
                                  test_field0->cval(cindex, 1, 0)*test_field0->cval(cindex, 1, 0) +
                                  test_field0->cval(cindex, 2, 0)*test_field0->cval(cindex, 2, 0) +
@@ -172,7 +176,7 @@ int symmetrize_test<rnumber>::do_work(void)
                                  test_field1->cval(cindex, 0, 1)*test_field1->cval(cindex, 0, 1) +
                                  test_field1->cval(cindex, 1, 1)*test_field1->cval(cindex, 1, 1) +
                                  test_field1->cval(cindex, 2, 1)*test_field1->cval(cindex, 2, 1));
-            double amplitude = sqrt((amplitude0 + amplitude1)/2);
+            double amplitude = std::sqrt((amplitude0 + amplitude1)/2);
             if (amplitude > 0)
             if (diff/amplitude > max_diff)
             {
@@ -180,9 +184,9 @@ int symmetrize_test<rnumber>::do_work(void)
                 ix = xindex;
                 iy = yindex + test_field0->clayout->starts[0];
                 iz = zindex;
-                k_at_max_diff = sqrt(k2);
-                a0 = sqrt(amplitude0);
-                a1 = sqrt(amplitude1);
+                k_at_max_diff = std::sqrt(k2);
+                a0 = std::sqrt(amplitude0);
+                a1 = std::sqrt(amplitude1);
             }
             });
     DEBUG_MSG("rseed1 found maximum relative difference %g at ix = %ld, iy = %ld, iz = %ld, wavenumber = %g, amplitudes %g %g\n",
@@ -243,8 +247,9 @@ int symmetrize_test<rnumber>::do_work(void)
             double diff_im0 = test_field0->cval(cindex, 0, 1) - test_field1->cval(cindex, 0, 1);
             double diff_im1 = test_field0->cval(cindex, 1, 1) - test_field1->cval(cindex, 1, 1);
             double diff_im2 = test_field0->cval(cindex, 2, 1) - test_field1->cval(cindex, 2, 1);
-            double diff = sqrt(diff_re0*diff_re0 + diff_re1*diff_re1 + diff_re2*diff_re2 +
-                               diff_im0*diff_im0 + diff_im1*diff_im1 + diff_im2*diff_im2);
+            double diff = std::sqrt(
+                    diff_re0*diff_re0 + diff_re1*diff_re1 + diff_re2*diff_re2 +
+                    diff_im0*diff_im0 + diff_im1*diff_im1 + diff_im2*diff_im2);
             double amplitude0 = (test_field0->cval(cindex, 0, 0)*test_field0->cval(cindex, 0, 0) +
                                  test_field0->cval(cindex, 1, 0)*test_field0->cval(cindex, 1, 0) +
                                  test_field0->cval(cindex, 2, 0)*test_field0->cval(cindex, 2, 0) +
@@ -257,7 +262,7 @@ int symmetrize_test<rnumber>::do_work(void)
                                  test_field1->cval(cindex, 0, 1)*test_field1->cval(cindex, 0, 1) +
                                  test_field1->cval(cindex, 1, 1)*test_field1->cval(cindex, 1, 1) +
                                  test_field1->cval(cindex, 2, 1)*test_field1->cval(cindex, 2, 1));
-            double amplitude = sqrt((amplitude0 + amplitude1)/2);
+            double amplitude = std::sqrt((amplitude0 + amplitude1)/2);
             if (amplitude > 0)
             if (diff/amplitude > max_diff)
             {
@@ -265,9 +270,9 @@ int symmetrize_test<rnumber>::do_work(void)
                 ix = xindex;
                 iy = yindex + test_field0->clayout->starts[0];
                 iz = zindex;
-                k_at_max_diff = sqrt(k2);
-                a0 = sqrt(amplitude0);
-                a1 = sqrt(amplitude1);
+                k_at_max_diff = std::sqrt(k2);
+                a0 = std::sqrt(amplitude0);
+                a1 = std::sqrt(amplitude1);
             }
             });
     DEBUG_MSG("rseed2 found maximum relative difference %g at ix = %ld, iy = %ld, iz = %ld, wavenumber = %g, amplitudes %g %g\n",
diff --git a/cpp/full_code/symmetrize_test.hpp b/cpp/full_code/symmetrize_test.hpp
index 69443d1a8dbeee28d035c5c2364f2d79dce8f8af..43b74e1b00690fffdc89e2f368e9a6c86a383867 100644
--- a/cpp/full_code/symmetrize_test.hpp
+++ b/cpp/full_code/symmetrize_test.hpp
@@ -29,9 +29,6 @@
 
 
 
-#include <cstdlib>
-#include "base.hpp"
-#include "kspace.hpp"
 #include "field.hpp"
 #include "full_code/test.hpp"
 
diff --git a/cpp/full_code/test.hpp b/cpp/full_code/test.hpp
index f7f9d8eeae30246d993edc75d5d5196b9ee5dabd..69aef319c3d0671704a0cf113001ea7fe61ebc9b 100644
--- a/cpp/full_code/test.hpp
+++ b/cpp/full_code/test.hpp
@@ -27,11 +27,6 @@
 #ifndef TEST_HPP
 #define TEST_HPP
 
-#include <cstdlib>
-#include <sys/types.h>
-#include <sys/stat.h>
-#include <vector>
-#include "base.hpp"
 #include "full_code/code_base.hpp"
 
 /** \brief base class for miscellaneous tests.
diff --git a/cpp/full_code/test_interpolation.hpp b/cpp/full_code/test_interpolation.hpp
index 7c9e0fff1f1d286865fd1f29b62005063f141e75..8881f307478719a5736db753259e282d1c157c0b 100644
--- a/cpp/full_code/test_interpolation.hpp
+++ b/cpp/full_code/test_interpolation.hpp
@@ -28,12 +28,8 @@
 
 
 
-#include <cstdlib>
-#include "base.hpp"
-#include "kspace.hpp"
 #include "full_code/test.hpp"
 #include "particles/particles_system_builder.hpp"
-#include "particles/particles_output_hdf5.hpp"
 #include "particles/particles_sampling.hpp"
 
 /** \brief Interpolation tester.
diff --git a/cpp/full_code/test_interpolation_methods.hpp b/cpp/full_code/test_interpolation_methods.hpp
index 62f661729cfbec51a25e6bb73819c8daac6dec66..93cee40529b9af2f7b25526a081a6ccc901ddbf5 100644
--- a/cpp/full_code/test_interpolation_methods.hpp
+++ b/cpp/full_code/test_interpolation_methods.hpp
@@ -28,9 +28,6 @@
 
 
 
-#include <cstdlib>
-#include "base.hpp"
-#include "kspace.hpp"
 #include "full_code/test.hpp"
 #include "particles/interpolation/field_tinterpolator.hpp"
 #include "particles/interpolation/particle_set.hpp"
diff --git a/cpp/full_code/test_particle_integration.hpp b/cpp/full_code/test_particle_integration.hpp
index 88a32dc19216887d619832fddcf6c5dc6ec18118..ed6ddee11fdf5d38a9841e29d72ac2754310cc71 100644
--- a/cpp/full_code/test_particle_integration.hpp
+++ b/cpp/full_code/test_particle_integration.hpp
@@ -28,13 +28,10 @@
 
 
 
-#include <cstdlib>
-#include "base.hpp"
-#include "kspace.hpp"
 #include "full_code/test.hpp"
-#include "particles/interpolation/field_tinterpolator.hpp"
-#include "particles/interpolation/particle_set.hpp"
 #include "particles/particle_solver.hpp"
+#include "particles/interpolation/particle_set.hpp"
+#include "particles/interpolation/field_tinterpolator.hpp"
 
 /** \brief Test of interpolation methods.
  *
diff --git a/cpp/full_code/test_tracer_set.cpp b/cpp/full_code/test_tracer_set.cpp
index ec8893577686cab0410ec49b396133535868b586..fdf2d5ab46b53a9c19701c64b3b90ddfe97c3154 100644
--- a/cpp/full_code/test_tracer_set.cpp
+++ b/cpp/full_code/test_tracer_set.cpp
@@ -23,17 +23,11 @@
 
 
 
-#include <string>
-#include <cmath>
-#include <random>
 #include "test_tracer_set.hpp"
-#include "particles/rhs/tracer_rhs.hpp"
 #include "particles/rhs/tracer_with_collision_counter_rhs.hpp"
 #include "particles/interpolation/particle_set.hpp"
 #include "particles/particle_solver.hpp"
 #include "particles/particles_input_random.hpp"
-#include "scope_timer.hpp"
-
 
 template <typename rnumber>
 int test_tracer_set<rnumber>::initialize(void)
diff --git a/cpp/full_code/test_tracer_set.hpp b/cpp/full_code/test_tracer_set.hpp
index b6d0ad4a0a68500639fea1f5b1e214a2fa61aad4..e78e12b4a3dbc13b6f18bbe8808708f251c53a2a 100644
--- a/cpp/full_code/test_tracer_set.hpp
+++ b/cpp/full_code/test_tracer_set.hpp
@@ -29,9 +29,6 @@
 
 
 
-#include <cstdlib>
-#include "base.hpp"
-#include "kspace.hpp"
 #include "field.hpp"
 #include "full_code/test.hpp"
 
diff --git a/cpp/full_code/write_filtered_test.cpp b/cpp/full_code/write_filtered_test.cpp
index ead8916e26d0b3b3c0be3e3e98985ca6dfc237b0..9ae70b2f33ff6cf4b835b296a94f9e7e89c138bc 100644
--- a/cpp/full_code/write_filtered_test.cpp
+++ b/cpp/full_code/write_filtered_test.cpp
@@ -22,13 +22,11 @@
 ******************************************************************************/
 
 
-
-#include <string>
-#include <cmath>
-#include <random>
 #include "write_filtered_test.hpp"
 #include "scope_timer.hpp"
 
+#include <cmath>
+#include <random>
 
 template <typename rnumber>
 int write_filtered_test<rnumber>::initialize(void)
diff --git a/cpp/full_code/write_filtered_test.hpp b/cpp/full_code/write_filtered_test.hpp
index 34cab4de7b1027b7520980365ea1f03a9df2f12e..5059731fdc762929e2f007f810ee833dcf0da30b 100644
--- a/cpp/full_code/write_filtered_test.hpp
+++ b/cpp/full_code/write_filtered_test.hpp
@@ -29,9 +29,6 @@
 
 
 
-#include <cstdlib>
-#include "base.hpp"
-#include "kspace.hpp"
 #include "field.hpp"
 #include "full_code/test.hpp"
 
diff --git a/cpp/full_code/write_rpressure.cpp b/cpp/full_code/write_rpressure.cpp
index 7833be2abea27bf8598f2cf93111e06fce504f47..4de5a2797a71703c24b1e20abc1c8c0775de2a25 100644
--- a/cpp/full_code/write_rpressure.cpp
+++ b/cpp/full_code/write_rpressure.cpp
@@ -1,10 +1,9 @@
-#include <string>
-#include <cmath>
-#include <hdf5.h>
 #include "write_rpressure.hpp"
 #include "scope_timer.hpp"
 #include "hdf5_tools.hpp"
 
+#include <cmath>
+
 template <typename rnumber>
 int write_rpressure<rnumber>::initialize(void)
 {
diff --git a/cpp/full_code/write_rpressure.hpp b/cpp/full_code/write_rpressure.hpp
index 2a14dcb6d58188cb5ad7b6d9333cff76ef963282..fdf54728b442b7e2de48ae26fbd33af023989a78 100644
--- a/cpp/full_code/write_rpressure.hpp
+++ b/cpp/full_code/write_rpressure.hpp
@@ -27,14 +27,6 @@
 #ifndef WRITE_RPRESSURE
 #define WRITE_RPRESSURE
 
-#include <cstdlib>
-#include <sys/types.h>
-#include <sys/stat.h>
-#include <vector>
-#include "base.hpp"
-#include "field.hpp"
-#include "field_binary_IO.hpp"
-#include "vorticity_equation.hpp"
 #include "full_code/NSVE_field_stats.hpp"
 
 /** generate and write out pressure field in real space **/
diff --git a/cpp/hdf5_tools.cpp b/cpp/hdf5_tools.cpp
index 22319198f4fbb15bcf5d3a048a51dce3c4c8a990..645c6c4af9c0e05154b24d32eb62ba90f546910d 100644
--- a/cpp/hdf5_tools.cpp
+++ b/cpp/hdf5_tools.cpp
@@ -1,8 +1,6 @@
 #include "hdf5_tools.hpp"
 #include "scope_timer.hpp"
 
-#include <limits>
-
 template <> hid_t hdf5_tools::hdf5_type_id<char>()
 {
     return H5T_NATIVE_CHAR;
@@ -341,6 +339,7 @@ int hdf5_tools::update_time_series_with_single_rank(
     H5Sget_simple_extent_dims(wspace, dims, NULL);
 
     int file_nvalues = -1;
+    variable_used_only_in_assert(file_nvalues);
     switch(ndims)
     {
         case 1:
diff --git a/cpp/hdf5_tools.hpp b/cpp/hdf5_tools.hpp
index a75c9a9585badfca264e94c8d45f67739e837cf6..8e8d870a39442cd2ff4b5e525eacb311d4eb225e 100644
--- a/cpp/hdf5_tools.hpp
+++ b/cpp/hdf5_tools.hpp
@@ -27,9 +27,11 @@
 #ifndef HDF5_TOOLS_HPP
 #define HDF5_TOOLS_HPP
 
+#include "base.hpp"
+
 #include <vector>
 #include <hdf5.h>
-#include "base.hpp"
+#include <string>
 
 namespace hdf5_tools
 {
diff --git a/cpp/kspace.cpp b/cpp/kspace.cpp
index dea114cfd7198ea842aa7d039c2af47e6ec9fd29..08f991edd55ec83a2a994b09846785404c1941f0 100644
--- a/cpp/kspace.cpp
+++ b/cpp/kspace.cpp
@@ -24,16 +24,13 @@
 
 
 
-#include <cmath>
-#include <cstdlib>
-#include <algorithm>
-#include <cassert>
-#include <stdexcept>
 #include "kspace.hpp"
 #include "scope_timer.hpp"
 #include "shared_array.hpp"
 #include "hdf5_tools.hpp"
 
+#include <cmath>
+
 template <field_backend be,
           kspace_dealias_type dt>
 template <field_components fc>
@@ -151,7 +148,7 @@ kspace<be, dt>::kspace(
                 const int nxmodes){
             if (k2 < this->kM2)
             {
-                const double knorm = sqrt(k2);
+                const double knorm = std::sqrt(k2);
                 kshell_local_thread.getMine()[int(knorm/this->dk)] += nxmodes*knorm;
                 nshell_local_thread.getMine()[int(knorm/this->dk)] += nxmodes;
             }
@@ -314,7 +311,7 @@ void kspace<be, dt>::ball_filter(
         typename fftw_interface<rnumber>::complex *__restrict__ a,
         const double ell)
 {
-    const double prefactor0 = double(3) / pow(ell/2, 3);
+    const double prefactor0 = double(3) / std::pow(ell/2, 3);
     this->CLOOP_K2(
             [&](const ptrdiff_t cindex,
                 const ptrdiff_t xindex,
@@ -323,12 +320,12 @@ void kspace<be, dt>::ball_filter(
                 const double k2){
                 if (k2 > 0)
                 {
-                    const double argument = sqrt(k2)*ell / 2;
-                    const double prefactor = prefactor0 / pow(k2, 1.5);
+                    const double argument = std::sqrt(k2)*ell / 2;
+                    const double prefactor = prefactor0 / std::pow(k2, 1.5);
                     for (unsigned int tcounter=0; tcounter<2*ncomp(fc); tcounter++)
                         ((rnumber*)a)[2*ncomp(fc)*cindex + tcounter] *= (
                             prefactor *
-                            (sin(argument) - argument * cos(argument)));
+                            (std::sin(argument) - argument * std::cos(argument)));
                 }
                 });
 }
@@ -359,12 +356,12 @@ void kspace<be, dt>::general_M_filter(
                 const double k2){
                 if (k2 > 0)
                 {
-                    const double argument = sqrt(k2)*ell;
+                    const double argument = std::sqrt(k2)*ell;
                     const double prefactor = prefactor0;
                     for (unsigned int tcounter=0; tcounter<2*ncomp(fc); tcounter++)
                         ((rnumber*)a)[2*ncomp(fc)*cindex + tcounter] *= (
-//                            prefactor * (exp(-0.5*pow((2.9*argument),(68.0*(pow(ell,0.74))))) * sqrt(1.0 + (pow((argument/0.06),3.8))/(1.0 + (pow((argument/0.057),3.8))))));
-                            prefactor * (exp(-1.7*argument)) * sqrt(1.0 + 0.68*(pow((argument/0.065),4.6))/(1.0 + (pow((argument/0.065),4.6)))));
+//                            prefactor * (std::exp(-0.5*std::pow((2.9*argument),(68.0*(std::pow(ell,0.74))))) * std::sqrt(1.0 + (std::pow((argument/0.06),3.8))/(1.0 + (std::pow((argument/0.057),3.8))))));
+                            prefactor * (std::exp(-1.7*argument)) * std::sqrt(1.0 + 0.68*(std::pow((argument/0.065),4.6))/(1.0 + (std::pow((argument/0.065),4.6)))));
                 }
                 });
 }
@@ -398,7 +395,7 @@ void kspace<be, dt>::Gauss_filter(
                 const double k2){
                 {
                     for (unsigned int tcounter=0; tcounter<2*ncomp(fc); tcounter++)
-                        ((rnumber*)a)[2*ncomp(fc)*cindex + tcounter] *= exp(prefactor*k2);
+                        ((rnumber*)a)[2*ncomp(fc)*cindex + tcounter] *= std::exp(prefactor*k2);
                 }
                 });
 }
@@ -594,10 +591,10 @@ void kspace<be, dt>::dealias(typename fftw_interface<rnumber>::complex *__restri
                     const ptrdiff_t xindex,
                     const ptrdiff_t yindex,
                     const ptrdiff_t zindex){
-                    const double kk2 = (pow(this->kx[xindex]/this->kMx, 2) +
-                                        pow(this->ky[yindex]/this->kMy, 2) +
-                                        pow(this->kz[zindex]/this->kMz, 2));
-                    const double tval = exp(-36.0 * (pow(kk2, 18)));
+                    const double kk2 = (std::pow(this->kx[xindex]/this->kMx, 2) +
+                                        std::pow(this->ky[yindex]/this->kMy, 2) +
+                                        std::pow(this->kz[zindex]/this->kMz, 2));
+                    const double tval = std::exp(-36.0 * (std::pow(kk2, 18)));
                     for (unsigned int tcounter=0; tcounter<2*ncomp(fc); tcounter++)
                         ((rnumber*)a)[2*ncomp(fc)*cindex + tcounter] *= tval;
                 });
@@ -632,7 +629,7 @@ void kspace<be, dt>::project_divfree(
                     double projected_size = 1;
                     if (maintain_energy)
                     {
-                        initial_size = sqrt(
+                        initial_size = std::sqrt(
                                 ((*(a + cindex*3  ))[0])*((*(a + cindex*3  ))[0]) +
                                 ((*(a + cindex*3+1))[0])*((*(a + cindex*3+1))[0]) +
                                 ((*(a + cindex*3+2))[0])*((*(a + cindex*3+2))[0]) +
@@ -648,7 +645,7 @@ void kspace<be, dt>::project_divfree(
                     }
                     if (maintain_energy)
                     {
-                        projected_size = sqrt(
+                        projected_size = std::sqrt(
                                 ((*(a + cindex*3  ))[0])*((*(a + cindex*3  ))[0]) +
                                 ((*(a + cindex*3+1))[0])*((*(a + cindex*3+1))[0]) +
                                 ((*(a + cindex*3+2))[0])*((*(a + cindex*3+2))[0]) +
@@ -698,24 +695,24 @@ void kspace<be, dt>::rotate_divfree(typename fftw_interface<rnumber>::complex *_
                             this->ky[yindex]*((*(a + cindex*3+1))[cc]) +
                             this->kz[zindex]*((*(a + cindex*3+2))[cc]) );
                 // now compute size of initial velocity vector
-                usize = sqrt(pow((*(a + cindex*3  ))[cc], 2) +
-                             pow((*(a + cindex*3+1))[cc], 2) +
-                             pow((*(a + cindex*3+2))[cc], 2));
+                usize = std::sqrt(std::pow((*(a + cindex*3  ))[cc], 2) +
+                                  std::pow((*(a + cindex*3+1))[cc], 2) +
+                                  std::pow((*(a + cindex*3+2))[cc], 2));
                 // finalize computation of cos Theta
                 if (usize != 0)
                 {
-                    cosTheta /= (usize * sqrt(k2));
+                    cosTheta /= (usize * std::sqrt(k2));
                     // now compute cross product
                     // cross product for complex vectors is complex conjugate of regular cross product
                     double cp[3];
                     cp[0] = (*(a + cindex*3+1))[cc]*this->kz[zindex] - (*(a + cindex*3+2))[cc]*this->ky[yindex];
                     cp[1] = (*(a + cindex*3+2))[cc]*this->kx[xindex] - (*(a + cindex*3+0))[cc]*this->kz[zindex];
                     cp[2] = (*(a + cindex*3+0))[cc]*this->ky[yindex] - (*(a + cindex*3+1))[cc]*this->kx[xindex];
-                    double cpsize = sqrt(cp[0]*cp[0] + cp[1]*cp[1] + cp[2]*cp[2]);
+                    double cpsize = std::sqrt(cp[0]*cp[0] + cp[1]*cp[1] + cp[2]*cp[2]);
                     cp[0] /= cpsize;
                     cp[1] /= cpsize;
                     cp[2] /= cpsize;
-                    double sinTheta = sqrt(1 - cosTheta*cosTheta);
+                    double sinTheta = std::sqrt(1 - cosTheta*cosTheta);
                     // we are actually interested in rotating the vector with pi/2 - Theta
                     double tmpdouble = cosTheta;
                     cosTheta = sinTheta;
@@ -765,10 +762,10 @@ void kspace<be, dt>::cospectrum(
             if (k2 <= this->kM2)
             {
                 double* spec_local = spec_local_thread.getMine();
-                const int tmp_int = int(sqrt(k2) / this->dk)*ncomp(fc)*ncomp(fc);
+                const int tmp_int = int(std::sqrt(k2) / this->dk)*ncomp(fc)*ncomp(fc);
                 for (hsize_t i=0; i<ncomp(fc); i++)
                 for (hsize_t j=0; j<ncomp(fc); j++){
-                    spec_local[tmp_int + i*ncomp(fc)+j] += pow(k2,wavenumber_exp/2.)*
+                    spec_local[tmp_int + i*ncomp(fc)+j] += std::pow(k2,wavenumber_exp/2.)*
 			nxmodes * (
                         (a[ncomp(fc)*cindex + i][0] * b[ncomp(fc)*cindex + j][0]) +
                         (a[ncomp(fc)*cindex + i][1] * b[ncomp(fc)*cindex + j][1]));
@@ -895,11 +892,11 @@ void kspace<be, dt>::cospectrum(
             if (k2 <= this->kM2)
             {
                 double* spec_local = spec_local_thread.getMine();
-                const int tmp_int = int(sqrt(k2) / this->dk)*ncomp(fc)*ncomp(fc);
+                const int tmp_int = int(std::sqrt(k2) / this->dk)*ncomp(fc)*ncomp(fc);
 
                 for (hsize_t i=0; i<ncomp(fc); i++)
                 for (hsize_t j=0; j<ncomp(fc); j++){
-                    spec_local[tmp_int + i*ncomp(fc)+j] += pow(k2, wavenumber_exp/2.)*
+                    spec_local[tmp_int + i*ncomp(fc)+j] += std::pow(k2, wavenumber_exp/2.)*
 			nxmodes * (
                         (a[ncomp(fc)*cindex + i][0] * a[ncomp(fc)*cindex + j][0]) +
                         (a[ncomp(fc)*cindex + i][1] * a[ncomp(fc)*cindex + j][1]));
@@ -954,7 +951,7 @@ double kspace<be, dt>::L2norm(
             &L2,
             1,
             MPI_DOUBLE, MPI_SUM, this->layout->comm);
-    return sqrt(L2 * this->dkx * this->dky * this->dkz);
+    return std::sqrt(L2 * this->dkx * this->dky * this->dkz);
 }
 
 
diff --git a/cpp/kspace.hpp b/cpp/kspace.hpp
index 6457c8bd0a65d229d154984c05341f69a7728f7b..dd401dd9ac97c3b2a1c87ffdc395f4a05f8db8d5 100644
--- a/cpp/kspace.hpp
+++ b/cpp/kspace.hpp
@@ -28,9 +28,6 @@
 
 #define KSPACE_HPP
 
-#include <hdf5.h>
-#include <vector>
-#include <string>
 #include "omputils.hpp"
 #include "fftw_interface.hpp"
 #include "field_layout.hpp"
diff --git a/cpp/particles/abstract_particles_output.hpp b/cpp/particles/abstract_particles_output.hpp
index af690dcbd72ca06cf756d1d639a6efaebfd6b677..23fb1e3f6590d2e36899761e996254527c0304ac 100644
--- a/cpp/particles/abstract_particles_output.hpp
+++ b/cpp/particles/abstract_particles_output.hpp
@@ -26,18 +26,12 @@
 #ifndef ABSTRACT_PARTICLES_OUTPUT
 #define ABSTRACT_PARTICLES_OUTPUT
 
-#include <memory>
-#include <vector>
-#include <cassert>
-#include <algorithm>
-#include <cstddef>
-
-#include "base.hpp"
-#include "particles_utils.hpp"
 #include "alltoall_exchanger.hpp"
-#include "scope_timer.hpp"
 #include "env_utils.hpp"
 
+#include <algorithm>
+#include <cstddef>
+
 /** \brief Class to handle distributed output of particle data
  *
  * \tparam partsize_t data type of particle index
diff --git a/cpp/particles/abstract_particles_system.hpp b/cpp/particles/abstract_particles_system.hpp
index 4edf20c7ca2d34fb023c6a0ea1e370c29e288df0..c8416621b950379ac8b4722fae685b62aa2113ce 100644
--- a/cpp/particles/abstract_particles_system.hpp
+++ b/cpp/particles/abstract_particles_system.hpp
@@ -26,14 +26,11 @@
 #ifndef ABSTRACT_PARTICLES_SYSTEM_HPP
 #define ABSTRACT_PARTICLES_SYSTEM_HPP
 
-#include <memory>
-#include <vector>
-
 //- Not generic to enable sampling begin
 #include "field.hpp"
-#include "kspace.hpp"
 //- Not generic to enable sampling end
 
+#include <memory>
 
 template <class partsize_t, class real_number>
 class abstract_particles_system {
diff --git a/cpp/particles/alltoall_exchanger.hpp b/cpp/particles/alltoall_exchanger.hpp
index 288ff47e3b68e9fc5c759ce3547ec530b76e74d8..2296c13ce46aa80fc270df7a9b157b3594e279cd 100644
--- a/cpp/particles/alltoall_exchanger.hpp
+++ b/cpp/particles/alltoall_exchanger.hpp
@@ -26,10 +26,6 @@
 #ifndef ALLTOALL_EXCHANGER_HPP
 #define ALLTOALL_EXCHANGER_HPP
 
-#include <mpi.h>
-#include <cassert>
-
-#include "base.hpp"
 #include "particles_utils.hpp"
 #include "scope_timer.hpp"
 
diff --git a/cpp/particles/inner/particles_inner_computer.cpp b/cpp/particles/inner/particles_inner_computer.cpp
index 2a8ac8da1ab7978a123569b9a052f9821423e9da..ecba7deb561581015dd94cab72651fcf719eeeb2 100644
--- a/cpp/particles/inner/particles_inner_computer.cpp
+++ b/cpp/particles/inner/particles_inner_computer.cpp
@@ -63,7 +63,7 @@ void particles_inner_computer<double, long long>::add_Lagrange_multipliers<6,6>(
             const long long idx0 = idx_part*6 + 3;
             const long long idx1 = idx_part*6 + 3;
             // check that orientation is unit vector:
-            double orientation_size = sqrt(
+            double orientation_size = std::sqrt(
                     pos_part[idx0+IDXC_X]*pos_part[idx0+IDXC_X] +
                     pos_part[idx0+IDXC_Y]*pos_part[idx0+IDXC_Y] +
                     pos_part[idx0+IDXC_Z]*pos_part[idx0+IDXC_Z]);
@@ -174,7 +174,7 @@ void particles_inner_computer<double, long long>::enforce_unit_orientation<6>(
     for(long long idx_part = 0 ; idx_part < nb_particles ; ++idx_part){
         const long long idx0 = idx_part*6 + 3;
         // compute orientation size:
-        double orientation_size = sqrt(
+        double orientation_size = std::sqrt(
                 pos_part[idx0+IDXC_X]*pos_part[idx0+IDXC_X] +
                 pos_part[idx0+IDXC_Y]*pos_part[idx0+IDXC_Y] +
                 pos_part[idx0+IDXC_Z]*pos_part[idx0+IDXC_Z]);
diff --git a/cpp/particles/inner/particles_inner_computer.hpp b/cpp/particles/inner/particles_inner_computer.hpp
index fd9dd678c27c2f1f878a40cecd7b2a1df4a0c727..65d3d6cdaee44b3989dba1daa4059b1b2692d743 100644
--- a/cpp/particles/inner/particles_inner_computer.hpp
+++ b/cpp/particles/inner/particles_inner_computer.hpp
@@ -28,7 +28,6 @@
 
 #include <cstring>
 #include <cassert>
-#include <iostream>
 
 template <class real_number, class partsize_t>
 class particles_inner_computer{
diff --git a/cpp/particles/inner/particles_inner_computer_2nd_order.hpp b/cpp/particles/inner/particles_inner_computer_2nd_order.hpp
index eebe162e627709bafd73b067d8660b6cdc2fd6fd..93d2d2df41ecd6eee792872e76acfd417e8f192f 100644
--- a/cpp/particles/inner/particles_inner_computer_2nd_order.hpp
+++ b/cpp/particles/inner/particles_inner_computer_2nd_order.hpp
@@ -26,8 +26,6 @@
 #ifndef PARTICLES_INNER_COMPUTER_2ND_ORDER_HPP
 #define PARTICLES_INNER_COMPUTER_2ND_ORDER_HPP
 
-#include <cstring>
-#include <cassert>
 
 /** \brief Computes action of Stokes drag on particles.
  *
diff --git a/cpp/particles/inner/particles_inner_computer_empty.hpp b/cpp/particles/inner/particles_inner_computer_empty.hpp
index aeff11c4e1c87256f216b3c70eda6fe78ea39f65..e09467bc2e9c604f1f8ace5f658afb3c1e242b88 100644
--- a/cpp/particles/inner/particles_inner_computer_empty.hpp
+++ b/cpp/particles/inner/particles_inner_computer_empty.hpp
@@ -26,9 +26,6 @@
 #ifndef PARTICLES_INNER_COMPUTER_EMPTY_HPP
 #define PARTICLES_INNER_COMPUTER_EMPTY_HPP
 
-#include <cstring>
-#include <cassert>
-
 template <class real_number, class partsize_t>
 class particles_inner_computer_empty{
 public:
diff --git a/cpp/particles/interpolation/abstract_particle_set.hpp b/cpp/particles/interpolation/abstract_particle_set.hpp
index af907584ddd55dedc7d8a6d5a6b59f6f01262928..2717975ff875ac6dced861f80bc3bf3bddd7efd8 100644
--- a/cpp/particles/interpolation/abstract_particle_set.hpp
+++ b/cpp/particles/interpolation/abstract_particle_set.hpp
@@ -28,18 +28,10 @@
 
 
 
-#include "field.hpp"
 #include "particles/p2p/p2p_distr_mpi.hpp"
 #include "particles/particles_sampling.hpp"
 #include "particles/abstract_particles_input.hpp"
 
-#include <array>
-#include <cassert>
-#include <vector>
-#include <memory>
-#include <string>
-
-
 
 /** \brief Brings together particle information with interpolation functionality.
  *
diff --git a/cpp/particles/interpolation/field_tinterpolator.hpp b/cpp/particles/interpolation/field_tinterpolator.hpp
index 2dc27b56a69d5462b131347b23c8fce1903864af..c308e635ce1f81e86904c2e3dc5d5f677f89ce77 100644
--- a/cpp/particles/interpolation/field_tinterpolator.hpp
+++ b/cpp/particles/interpolation/field_tinterpolator.hpp
@@ -26,14 +26,11 @@
 #ifndef FIELD_TINTERPOLATOR_HPP
 #define FIELD_TINTERPOLATOR_HPP
 
-#include <array>
-#include <cassert>
-#include <vector>
-#include "field.hpp"
 #include "particles/particles_distr_mpi.hpp"
 #include "particles/particles_output_hdf5.hpp"
 #include "particles/interpolation/abstract_particle_set.hpp"
 
+
 enum temporal_interpolation_type {NONE, LINEAR, CUBIC};
 
 template <typename rnumber,
diff --git a/cpp/particles/interpolation/particle_set.hpp b/cpp/particles/interpolation/particle_set.hpp
index 5f5eb7aa1bf30ead417c0409fbb6d5a768d294f6..52bb971630d5942dcf777f7928fe9e2325fc02b6 100644
--- a/cpp/particles/interpolation/particle_set.hpp
+++ b/cpp/particles/interpolation/particle_set.hpp
@@ -26,21 +26,11 @@
 #ifndef PARTICLE_SET_HPP
 #define PARTICLE_SET_HPP
 
-#include <array>
-#include <cassert>
-#include <vector>
-#include <set>
-#include <memory>
-#include "field_layout.hpp"
 #include "field.hpp"
 #include "particles/particles_distr_mpi.hpp"
-#include "particles/abstract_particles_output.hpp"
-#include "particles/abstract_particles_input.hpp"
 #include "particles/interpolation/particles_field_computer.hpp"
 #include "particles/interpolation/particles_generic_interp.hpp"
 #include "particles/interpolation/abstract_particle_set.hpp"
-#include "particles/p2p/p2p_distr_mpi.hpp"
-
 
 
 /** Practical implementation of particle sets.
diff --git a/cpp/particles/interpolation/particles_field_computer.hpp b/cpp/particles/interpolation/particles_field_computer.hpp
index f81f8e53a692904fe35ffb842d60eff07e7bb31d..bf504fae17018807ce2313e80746ce45d9c5d0f9 100644
--- a/cpp/particles/interpolation/particles_field_computer.hpp
+++ b/cpp/particles/interpolation/particles_field_computer.hpp
@@ -26,13 +26,13 @@
 #ifndef PARTICLES_FIELD_COMPUTER_HPP
 #define PARTICLES_FIELD_COMPUTER_HPP
 
+#include "scope_timer.hpp"
+#include "particles/particles_utils.hpp"
+
 #include <array>
 #include <utility>
 #include <cmath>
 
-#include "scope_timer.hpp"
-#include "particles/particles_utils.hpp"
-
 
 
 /** \brief Computes interpolation kernel sum operation.
diff --git a/cpp/particles/lock_free_bool_array.hpp b/cpp/particles/lock_free_bool_array.hpp
index a08abb808217c1952a8b9b6893e9b788534b8d25..61673c1be4829e86c1b32f64865de24a3c6ba211 100644
--- a/cpp/particles/lock_free_bool_array.hpp
+++ b/cpp/particles/lock_free_bool_array.hpp
@@ -29,7 +29,6 @@
 #include <vector>
 #include <atomic>
 #include <unistd.h>
-#include <cstdio>
 #include <omp.h>
 
 class lock_free_bool_array{
diff --git a/cpp/particles/p2p/p2p_computer.hpp b/cpp/particles/p2p/p2p_computer.hpp
index 5a399e515da2bf59bdfcb3874c7591daf0af1b4c..8e048cbba85508edd1c45a6e4cb9eecff0ef6ba1 100644
--- a/cpp/particles/p2p/p2p_computer.hpp
+++ b/cpp/particles/p2p/p2p_computer.hpp
@@ -26,7 +26,6 @@
 #ifndef P2P_COMPUTER_HPP
 #define P2P_COMPUTER_HPP
 
-#include <cstring>
 #include <cassert>
 
 template <class real_number, class partsize_t>
diff --git a/cpp/particles/p2p/p2p_computer_empty.hpp b/cpp/particles/p2p/p2p_computer_empty.hpp
index b13a95e1f2bb8fd4f51f51c4cb3b6f74e00b079b..f1f41448d14720dbed5f7c4d76f9b939b2e2e9b7 100644
--- a/cpp/particles/p2p/p2p_computer_empty.hpp
+++ b/cpp/particles/p2p/p2p_computer_empty.hpp
@@ -26,8 +26,6 @@
 #ifndef P2P_COMPUTER_EMPTY_HPP
 #define P2P_COMPUTER_EMPTY_HPP
 
-#include <cstring>
-
 template <class real_number, class partsize_t>
 class p2p_computer_empty{
 public:
diff --git a/cpp/particles/p2p/p2p_distr_mpi.hpp b/cpp/particles/p2p/p2p_distr_mpi.hpp
index 3051bcee05b5d6078abdbac571cc4500ff6af475..3a0a305f8cb79047b7f6caa1b20f4e45969fd41e 100644
--- a/cpp/particles/p2p/p2p_distr_mpi.hpp
+++ b/cpp/particles/p2p/p2p_distr_mpi.hpp
@@ -26,21 +26,14 @@
 #ifndef P2P_DISTR_MPI_HPP
 #define P2P_DISTR_MPI_HPP
 
-#include <mpi.h>
-
-#include <vector>
-#include <memory>
-#include <cassert>
-
-#include <type_traits>
-#include <omp.h>
-#include <algorithm>
-
 #include "scope_timer.hpp"
 #include "particles/particles_utils.hpp"
 #include "particles/p2p/p2p_tree.hpp"
 #include "particles/lock_free_bool_array.hpp"
 
+#include <type_traits>
+#include <algorithm>
+
 template <class partsize_t, class real_number>
 class p2p_distr_mpi {
 protected:
diff --git a/cpp/particles/p2p/p2p_ghost_collisions.hpp b/cpp/particles/p2p/p2p_ghost_collisions.hpp
index 4124acfacdde9fff3144c8a11ed219c0405a8a6f..82d3113d46eb3a24bd62d5a7d1685e5ea36f7e5a 100644
--- a/cpp/particles/p2p/p2p_ghost_collisions.hpp
+++ b/cpp/particles/p2p/p2p_ghost_collisions.hpp
@@ -116,7 +116,7 @@ public:
     double rod_distance(double dist_pow2, double xp, double xq, double pq, double t, double s){
         double min_distance;
         min_distance = dist_pow2 + 0.25*this->get_cylinder_length()*this->get_cylinder_length() * ( t*t + s*s - 2.0*t*s*pq ) + this->get_cylinder_length() * ( t*xp - s*xq );
-        min_distance = sqrt( min_distance );
+        min_distance = std::sqrt( min_distance );
         return min_distance;
     }
 
@@ -241,7 +241,7 @@ public:
                     const double x = xseparation;
                     const double y = yseparation;
                     const double z = zseparation;
-                    /* const double r = sqrt(dist_pow2); This variable is not needed. */
+                    /* const double r = std::sqrt(dist_pow2); This variable is not needed. */
 
                     /* p and q are the orientation vectors of the first and second particles. */
                     px = pos_part1[IDXC_X+3];
@@ -256,7 +256,7 @@ public:
                     xq = x * qx + y * qy + z * qz;
                     /* t and s parametrize the two rods. Find min distance: */
                     if( pq == 1.0 ){
-                        min_distance = sqrt(dist_pow2-xp*xp);
+                        min_distance = std::sqrt(dist_pow2-xp*xp);
                     }
                     else{
                     t = 2.0/(this->get_cylinder_length()*(pq*pq-1.0))*(xp-pq*xq);
@@ -338,7 +338,7 @@ public:
                     p0[IDXC_X] = p1[IDXC_Y]*p2[IDXC_Z] - p1[IDXC_Z]*p2[IDXC_Y];
                     p0[IDXC_Y] = p1[IDXC_Z]*p2[IDXC_X] - p1[IDXC_X]*p2[IDXC_Z];
                     p0[IDXC_Z] = p1[IDXC_X]*p2[IDXC_Y] - p1[IDXC_Y]*p2[IDXC_X];
-                    p0_norm = sqrt(p0[IDXC_X]*p0[IDXC_X]+p0[IDXC_Y]*p0[IDXC_Y]+p0[IDXC_Z]*p0[IDXC_Z]);
+                    p0_norm = std::sqrt(p0[IDXC_X]*p0[IDXC_X]+p0[IDXC_Y]*p0[IDXC_Y]+p0[IDXC_Z]*p0[IDXC_Z]);
                     p0[IDXC_X] = p0[IDXC_X] / p0_norm;
                     p0[IDXC_Y] = p0[IDXC_Y] / p0_norm;
                     p0[IDXC_Z] = p0[IDXC_Z] / p0_norm;
@@ -355,8 +355,8 @@ public:
                         t1[0] = x1p0;
                         t1[1] = x1p0;
                     } else if (det > 0.0) {
-                        t1[0] =  x1p0 - sqrt(det);
-                        t1[1] =  x1p0 + sqrt(det);
+                        t1[0] =  x1p0 - std::sqrt(det);
+                        t1[1] =  x1p0 + std::sqrt(det);
                     } else {return;}
 
                     det = x2p0 * x2p0 + 0.25 * this->get_disk_width() * this->get_disk_width() - x2_x0_2;
@@ -364,8 +364,8 @@ public:
                         t2[0] = x2p0;
                         t2[1] = x2p0;
                     } else if (det > 0.0) {
-                        t2[0] =  x2p0 - sqrt(det);
-                        t2[1] =  x2p0 + sqrt(det);
+                        t2[0] =  x2p0 - std::sqrt(det);
+                        t2[1] =  x2p0 + std::sqrt(det);
                     } else {return;}
 
                     if (( t1[1]>=t2[0] ) and
diff --git a/cpp/particles/p2p/p2p_merge_collisions.hpp b/cpp/particles/p2p/p2p_merge_collisions.hpp
index 4b728c8f3c3db1c26e41947fe1a23325fbe688e8..071c2e92f0cc6c751fb269e343e2c87b4cef2405 100644
--- a/cpp/particles/p2p/p2p_merge_collisions.hpp
+++ b/cpp/particles/p2p/p2p_merge_collisions.hpp
@@ -23,8 +23,6 @@
 #ifndef P2P_MERGE_COLLISIONS_HPP
 #define P2P_MERGE_COLLISIONS_HPP
 
-#include <cstring>
-
 template <class real_number, class partsize_t>
 class p2p_merge_collisions{
     long int collision_counter;
diff --git a/cpp/particles/p2p/p2p_tree.hpp b/cpp/particles/p2p/p2p_tree.hpp
index 9edf7a17050470d3c10d434a1a5816c1564771ca..2022f0abc9615b2f1c4118f5e5ae369e98cc1c59 100644
--- a/cpp/particles/p2p/p2p_tree.hpp
+++ b/cpp/particles/p2p/p2p_tree.hpp
@@ -27,7 +27,6 @@
 #define P2P_TREE_HPP
 
 #include <unordered_map>
-#include <vector>
 #include <array>
 
 template <class CellClass>
diff --git a/cpp/particles/particle_solver.hpp b/cpp/particles/particle_solver.hpp
index 54bbc3c846022c00f220d9f427e121241fd10223..55e753110c7d43e2de54ce43b4ff93aa59eb23a3 100644
--- a/cpp/particles/particle_solver.hpp
+++ b/cpp/particles/particle_solver.hpp
@@ -26,13 +26,10 @@
 #ifndef PARTICLE_SOLVER_HPP
 #define PARTICLE_SOLVER_HPP
 
-#include <cassert>
-#include "particles/interpolation/abstract_particle_set.hpp"
 #include "particles/abstract_particle_rhs.hpp"
 #include "particles/particles_output_hdf5.hpp"
 
 
-
 /** Time-stepping and checkpointing functionality for particle systems.
  *
  * This class implements the universal functionality of a particle solver:
diff --git a/cpp/particles/particles_adams_bashforth.hpp b/cpp/particles/particles_adams_bashforth.hpp
index f760d473890d8210a5b9927e524bad0a6d1df761..4db24b3892523abe5eb21bb70db81f524c5ce83a 100644
--- a/cpp/particles/particles_adams_bashforth.hpp
+++ b/cpp/particles/particles_adams_bashforth.hpp
@@ -26,9 +26,6 @@
 #ifndef PARTICLES_ADAMS_BASHFORTH_HPP
 #define PARTICLES_ADAMS_BASHFORTH_HPP
 
-#include <stdexcept>
-#include <omp.h>
-
 #include "scope_timer.hpp"
 #include "particles_utils.hpp"
 
diff --git a/cpp/particles/particles_distr_mpi.hpp b/cpp/particles/particles_distr_mpi.hpp
index 8b25297e61e09935ab2632a1b4b8f4e6343231a9..f330e638c039b03bb25e2ce7fed2dc7799b7627d 100644
--- a/cpp/particles/particles_distr_mpi.hpp
+++ b/cpp/particles/particles_distr_mpi.hpp
@@ -26,18 +26,11 @@
 #ifndef PARTICLES_DISTR_MPI_HPP
 #define PARTICLES_DISTR_MPI_HPP
 
-#include <mpi.h>
-
-#include <vector>
-#include <memory>
-#include <cassert>
-
-#include <type_traits>
-#include <omp.h>
-
 #include "scope_timer.hpp"
 #include "particles_utils.hpp"
 
+#include <array>
+
 
 /** Particle interpolation/redistribution functionality
  *
diff --git a/cpp/particles/particles_input_grid.hpp b/cpp/particles/particles_input_grid.hpp
index 784c08e0623b602dc6e81a7b4bd13ed3c8a8dfe0..b0b8f28792a5ef7fb92aa6b6ae46aa51925df07a 100644
--- a/cpp/particles/particles_input_grid.hpp
+++ b/cpp/particles/particles_input_grid.hpp
@@ -26,9 +26,10 @@
 #ifndef PARTICLES_INPUT_GRID_HPP
 #define PARTICLES_INPUT_GRID_HPP
 
+#include "particles/abstract_particles_input.hpp"
+
 #include <mpi.h>
 #include <cassert>
-#include "particles/abstract_particles_input.hpp"
 
 template <class partsize_t, class particle_rnumber, int state_size>
 class particles_input_grid: public abstract_particles_input<partsize_t, particle_rnumber>
diff --git a/cpp/particles/particles_input_hdf5.hpp b/cpp/particles/particles_input_hdf5.hpp
index b1ea662ed14eaf4e725c532429eb631ade4ef8b7..c21943270c2889b4dad486b678b6dac9aa6dabea 100644
--- a/cpp/particles/particles_input_hdf5.hpp
+++ b/cpp/particles/particles_input_hdf5.hpp
@@ -26,20 +26,10 @@
 #ifndef PARTICLES_INPUT_HDF5_HPP
 #define PARTICLES_INPUT_HDF5_HPP
 
-#include <tuple>
-#include <mpi.h>
-#include <hdf5.h>
-#include <cassert>
-#include <vector>
-
 #include "abstract_particles_input.hpp"
-#include "base.hpp"
 #include "alltoall_exchanger.hpp"
-#include "particles/particles_utils.hpp"
-#include "scope_timer.hpp"
 #include "hdf5_tools.hpp"
 
-
 template <class partsize_t, class real_number, int size_particle_positions, int size_particle_rhs>
 class particles_input_hdf5 : public abstract_particles_input<partsize_t, real_number> {
     const std::string filename;
diff --git a/cpp/particles/particles_input_random.hpp b/cpp/particles/particles_input_random.hpp
index 0b40ae57f70729b82984c4c5b4591568e4078485..35d6d4c7726bd49f39d8443e83a0c7f4f893cae3 100644
--- a/cpp/particles/particles_input_random.hpp
+++ b/cpp/particles/particles_input_random.hpp
@@ -26,10 +26,11 @@
 #ifndef PARTICLES_INPUT_RANDOM_HPP
 #define PARTICLES_INPUT_RANDOM_HPP
 
+#include "particles/abstract_particles_input.hpp"
+
 #include <mpi.h>
 #include <cassert>
 #include <random>
-#include "particles/abstract_particles_input.hpp"
 
 template <class partsize_t, class particle_rnumber, int state_size>
 class particles_input_random: public abstract_particles_input<partsize_t, particle_rnumber>
diff --git a/cpp/particles/particles_output_hdf5.hpp b/cpp/particles/particles_output_hdf5.hpp
index 0ef809e63114948509c96f5cbccc5cdd615b3599..48b953c7a2019d6440db4d70c3c44b98b99fc812 100644
--- a/cpp/particles/particles_output_hdf5.hpp
+++ b/cpp/particles/particles_output_hdf5.hpp
@@ -26,14 +26,11 @@
 #ifndef PARTICLES_OUTPUT_HDF5_HPP
 #define PARTICLES_OUTPUT_HDF5_HPP
 
-#include <memory>
-#include <vector>
+#include "abstract_particles_output.hpp"
+
 #include <hdf5.h>
 #include <sys/stat.h>
 
-#include "abstract_particles_output.hpp"
-#include "scope_timer.hpp"
-
 template <class partsize_t,
           class real_number,
           int size_particle_positions>
diff --git a/cpp/particles/particles_output_mpiio.hpp b/cpp/particles/particles_output_mpiio.hpp
index ee9bd9cee42d3973e1e9e1bb93eb67a6329836d3..7a3da6cf42a0814fe5452ab97095cd39b68d6dff 100644
--- a/cpp/particles/particles_output_mpiio.hpp
+++ b/cpp/particles/particles_output_mpiio.hpp
@@ -26,14 +26,7 @@
 #ifndef PARTICLES_OUTPUT_MPIIO
 #define PARTICLES_OUTPUT_MPIIO
 
-#include <memory>
-#include <vector>
-#include <string>
-#include <cassert>
-
 #include "abstract_particles_output.hpp"
-#include "scope_timer.hpp"
-#include "particles_utils.hpp"
 
 template <class partsize_t, class real_number, int size_particle_positions, int size_particle_rhs>
 class particles_output_mpiio : public abstract_particles_output<partsize_t, real_number, real_number, size_particle_positions>{
diff --git a/cpp/particles/particles_sampling.hpp b/cpp/particles/particles_sampling.hpp
index 8e83a5a53690b6cb3abaddaafe2245b1f448d164..a341929612e63e859f18534bbc7c8875a79fd2f6 100644
--- a/cpp/particles/particles_sampling.hpp
+++ b/cpp/particles/particles_sampling.hpp
@@ -26,15 +26,9 @@
 #ifndef PARTICLES_SAMPLING_HPP
 #define PARTICLES_SAMPLING_HPP
 
-#include <memory>
-#include <string>
-
 #include "abstract_particles_system.hpp"
 #include "particles_output_sampling_hdf5.hpp"
 
-#include "field.hpp"
-#include "kspace.hpp"
-
 
 template <class partsize_t, class particles_rnumber, class rnumber, field_backend be, field_components fc>
 void sample_from_particles_system(const field<rnumber, be, fc>& in_field, // a pointer to a field<rnumber, FFTW, fc>
diff --git a/cpp/particles/particles_system.hpp b/cpp/particles/particles_system.hpp
index 9589e4dcc2738c554f7ff364235b7bc92d8a155a..d58f92f8bc4fa2942a30a91625385781f5f4c138 100644
--- a/cpp/particles/particles_system.hpp
+++ b/cpp/particles/particles_system.hpp
@@ -26,12 +26,6 @@
 #ifndef PARTICLES_SYSTEM_HPP
 #define PARTICLES_SYSTEM_HPP
 
-#include <array>
-#include <set>
-#include <algorithm>
-#include <cmath>
-
-#include "particles/abstract_particles_system.hpp"
 #include "particles/abstract_particles_system_with_p2p.hpp"
 #include "particles/particles_distr_mpi.hpp"
 #include "particles/particles_output_hdf5.hpp"
@@ -39,10 +33,11 @@
 #include "particles/interpolation/particles_field_computer.hpp"
 #include "particles/abstract_particles_input.hpp"
 #include "particles/particles_adams_bashforth.hpp"
-#include "scope_timer.hpp"
 
 #include "particles/p2p/p2p_distr_mpi.hpp"
 
+#include <set>
+
 template <class partsize_t, class real_number, class field_rnumber, class field_class, class interpolator_class, int interp_neighbours,
           int size_particle_positions, int size_particle_rhs, class p2p_computer_class, class particles_inner_computer_class>
 class particles_system : public abstract_particles_system_with_p2p<partsize_t, real_number, p2p_computer_class> {
diff --git a/cpp/particles/particles_system_builder.hpp b/cpp/particles/particles_system_builder.hpp
index dd2b775032fa5ffe8e3439a500419e83d219de1e..fbbf4a41994f7f99e6b09b41d6ffe73d52cdec24 100644
--- a/cpp/particles/particles_system_builder.hpp
+++ b/cpp/particles/particles_system_builder.hpp
@@ -26,22 +26,12 @@
 #ifndef PARTICLES_SYSTEM_BUILDER_HPP
 #define PARTICLES_SYSTEM_BUILDER_HPP
 
-#include <string>
-
-#include <cmath>
-#include "particles/abstract_particles_system.hpp"
 #include "particles/particles_system.hpp"
 #include "particles/particles_input_hdf5.hpp"
 #include "particles/interpolation/particles_generic_interp.hpp"
 #include "particles/p2p/p2p_computer_empty.hpp"
 #include "particles/inner/particles_inner_computer_empty.hpp"
 
-#include "field.hpp"
-#include "kspace.hpp"
-
-const int MAX_INTERPOLATION_NEIGHBOURS=7;
-const int MAX_INTERPOLATION_SMOOTHNESS=4;
-
 //////////////////////////////////////////////////////////////////////////////
 ///
 /// Double template "for"
@@ -306,8 +296,8 @@ inline std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>>
         MPI_Comm mpi_comm,
         const int in_current_iteration){
     return Template_double_for_if::evaluate<std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>>,
-                       int, 1, MAX_INTERPOLATION_NEIGHBOURS, 1, // interpolation_size
-                       int, 0, MAX_INTERPOLATION_SMOOTHNESS, 1, // spline_mode
+                       int, MIN_INTERPOLATION_NEIGHBOURS, MAX_INTERPOLATION_NEIGHBOURS, 1, // interpolation_size
+                       int, MIN_INTERPOLATION_SMOOTHNESS, MAX_INTERPOLATION_SMOOTHNESS, 1, // spline_mode
                        particles_system_build_container<partsize_t, field_rnumber,be,fc,particles_rnumber,
                                                         p2p_computer_empty<particles_rnumber,partsize_t>,
                                                         particles_inner_computer_empty<particles_rnumber,partsize_t>,
@@ -336,8 +326,8 @@ inline std::unique_ptr<abstract_particles_system_with_p2p<partsize_t, particles_
         particles_inner_computer_class inner_computer,
         const particles_rnumber cutoff){
     return Template_double_for_if::evaluate<std::unique_ptr<abstract_particles_system_with_p2p<partsize_t, particles_rnumber, p2p_computer_class>>,
-                       int, 1, MAX_INTERPOLATION_NEIGHBOURS, 1, // interpolation_size
-                       int, 0, MAX_INTERPOLATION_SMOOTHNESS, 1, // spline_mode
+                       int, MIN_INTERPOLATION_NEIGHBOURS, MAX_INTERPOLATION_NEIGHBOURS, 1, // interpolation_size
+                       int, MIN_INTERPOLATION_SMOOTHNESS, MAX_INTERPOLATION_SMOOTHNESS, 1, // spline_mode
                        particles_system_build_container<partsize_t, field_rnumber,be,fc,particles_rnumber,
                                                         p2p_computer_class,
                                                         particles_inner_computer_class,
diff --git a/cpp/particles/rhs/deformation_tensor_rhs.hpp b/cpp/particles/rhs/deformation_tensor_rhs.hpp
index bac9624852db2482e0e79e73bf24e9986fc5b371..91c880d57d9ef62c177636ce5077f43e066ea2ff 100644
--- a/cpp/particles/rhs/deformation_tensor_rhs.hpp
+++ b/cpp/particles/rhs/deformation_tensor_rhs.hpp
@@ -26,7 +26,6 @@
 #ifndef DEFORMATION_TENSOR_RHS_HPP
 #define DEFORMATION_TENSOR_RHS_HPP
 
-#include "particles/particles_utils.hpp"
 #include "particles/interpolation/field_tinterpolator.hpp"
 #include "particles/abstract_particle_rhs.hpp"
 #include <gsl/gsl_linalg.h>
diff --git a/cpp/particles/rhs/tracer_with_collision_counter_rhs.hpp b/cpp/particles/rhs/tracer_with_collision_counter_rhs.hpp
index 4fdc0686586003b913f3ef09dd70bba69365b5e4..892954147fc121e3b6721e5081969a732da31ce0 100644
--- a/cpp/particles/rhs/tracer_with_collision_counter_rhs.hpp
+++ b/cpp/particles/rhs/tracer_with_collision_counter_rhs.hpp
@@ -26,8 +26,6 @@
 #ifndef TRACER_WITH_COLLISION_COUNTER_RHS_HPP
 #define TRACER_WITH_COLLISION_COUNTER_RHS_HPP
 
-#include "particles/interpolation/field_tinterpolator.hpp"
-#include "particles/abstract_particle_rhs.hpp"
 #include "particles/rhs/tracer_rhs.hpp"
 #include "particles/p2p/p2p_ghost_collisions.hpp"
 
diff --git a/cpp/scope_timer.hpp b/cpp/scope_timer.hpp
index d6085050ef233abf50d27f6c17d332c78e14ee27..c149be61cb873a60285c14432c696e92d188d3c8 100644
--- a/cpp/scope_timer.hpp
+++ b/cpp/scope_timer.hpp
@@ -25,25 +25,23 @@
 #ifndef SCOPE_TIMER_HPP
 #define SCOPE_TIMER_HPP
 
+#include "base.hpp"
+#include "turtle_timer.hpp"
+
 #include <memory>
 #include <iostream>
 #include <vector>
 #include <stack>
 #include <string>
 #include <limits>
-#include <cassert>
 #include <sstream>
 #include <unordered_map>
-#include <mpi.h>
 #include <cstring>
 #include <stdexcept>
 #include <omp.h>
 #include <iomanip>
 #include <fstream>
 
-#include "base.hpp"
-#include "turtle_timer.hpp"
-
 //< To add it as friend of EventManager
 class ScopeEvent;
 
diff --git a/cpp/shared_array.hpp b/cpp/shared_array.hpp
index 2a17eaf89d62e152d3dfb798dc0ff3137b7bf56c..68029157f0b7ed2d16a4303e309c2f5ed2568a89 100644
--- a/cpp/shared_array.hpp
+++ b/cpp/shared_array.hpp
@@ -28,7 +28,6 @@
 
 #include <omp.h>
 #include <functional>
-#include <iostream>
 
 template <class ValueType, size_t dim>
 void shared_array_zero_initializer(ValueType *data) {
diff --git a/cpp/spectrum_function.hpp b/cpp/spectrum_function.hpp
index 7be4474c0a6151e4eed81248140ad95fc2cdec51..57b6c42654441cedd7a2364ede8b7fda562902fc 100644
--- a/cpp/spectrum_function.hpp
+++ b/cpp/spectrum_function.hpp
@@ -1,10 +1,10 @@
 #ifndef SPECTRUM_FUNCTION_CPP
 #define SPECTRUM_FUNCTION_CPP
 
+#include "kspace.hpp"
+
 #include<cmath>
-#include<vector>
 
-#include "kspace.hpp"
 
 template <field_backend be,
           kspace_dealias_type dt>
@@ -28,7 +28,7 @@ class spectrum_function
         double operator()(double kvalue)
         {
             assert(kvalue >= double(0));
-            int index = floor(kvalue / this->kk->dk);
+            int index = std::floor(kvalue / this->kk->dk);
             assert(index < this->values.size());
             return this->values[index];
         }
diff --git a/cpp/turtle_config.hpp.in b/cpp/turtle_config.hpp.in
index 8085ce38e2652f45c172a6d2c75b2d54db7557f3..59f877c226bc9ffd558f3e1ef027efc7bf1aad1c 100644
--- a/cpp/turtle_config.hpp.in
+++ b/cpp/turtle_config.hpp.in
@@ -31,5 +31,10 @@
 #cmakedefine USE_FFTW_OMP
 #cmakedefine USE_TIMING_OUTPUT
 
+const int MIN_INTERPOLATION_NEIGHBOURS=@MIN_INTERPOLATION_NEIGHBOURS@;
+const int MIN_INTERPOLATION_SMOOTHNESS=@MIN_INTERPOLATION_SMOOTHNESS@;
+const int MAX_INTERPOLATION_NEIGHBOURS=@MAX_INTERPOLATION_NEIGHBOURS@;
+const int MAX_INTERPOLATION_SMOOTHNESS=@MAX_INTERPOLATION_SMOOTHNESS@;
+
 #endif//TURTLE_CONFIG_HPP
 
diff --git a/cpp/vorticity_equation.cpp b/cpp/vorticity_equation.cpp
index 9539b6167c3800721f6b6b56094e12ce9a446739..e13fe231c3902f091987ae89acb80432b441f298 100644
--- a/cpp/vorticity_equation.cpp
+++ b/cpp/vorticity_equation.cpp
@@ -24,17 +24,12 @@
 
 
 
-#include <limits>
-#include <cassert>
-#include <cmath>
-#include <cstring>
-#include <stdexcept>
 #include "fftw_tools.hpp"
 #include "vorticity_equation.hpp"
 #include "scope_timer.hpp"
 #include "shared_array.hpp"
 
-
+#include <sys/stat.h>
 
 template <class rnumber,
           field_backend be>
@@ -270,9 +265,10 @@ void vorticity_equation<rnumber, be>::add_field_band(
                     const ptrdiff_t xindex,
                     const ptrdiff_t yindex,
                     const ptrdiff_t zindex){
-        const double knorm = sqrt(this->kk->kx[xindex]*this->kk->kx[xindex] +
-                                  this->kk->ky[yindex]*this->kk->ky[yindex] +
-                                  this->kk->kz[zindex]*this->kk->kz[zindex]);
+        const double knorm = std::sqrt(
+                this->kk->kx[xindex]*this->kk->kx[xindex] +
+                this->kk->ky[yindex]*this->kk->ky[yindex] +
+                this->kk->kz[zindex]*this->kk->kz[zindex]);
         if ((k0 <= knorm) &&
             (k1 >= knorm))
             for (int c=0; c<3; c++)
@@ -320,7 +316,7 @@ void vorticity_equation<rnumber, be>::add_forcing(
     if (this->forcing_type == std::string("Kolmogorov_and_compensated_drag"))
     {
         double amplitude = this->famplitude * (
-                1 + this->friction_coefficient / sqrt(this->fmode  * this->famplitude));
+                1 + this->friction_coefficient / std::sqrt(this->fmode  * this->famplitude));
         this->add_Kolmogorov_forcing(dst, this->fmode, amplitude);
         this->add_field_band(
                 dst, vort_field,
@@ -352,7 +348,7 @@ void vorticity_equation<rnumber, be>::add_forcing(
                         const ptrdiff_t zindex,
                         const double k2,
                         const int nxmodes){
-            const double knorm = sqrt(k2);
+            const double knorm = std::sqrt(k2);
             if ((k2 > 0) &&
                 (this->fk0 <= knorm) &&
                 (this->fk1 >= knorm))
@@ -382,7 +378,7 @@ void vorticity_equation<rnumber, be>::add_forcing(
         double temp_famplitude = 0;
         if (this->forcing_type == std::string("sinusoidal_energy_injection_rate")){
             temp_famplitude = (this->injection_rate + this->variation_strength
-                                *sin(t/this->variation_time_scale))
+                                *std::sin(t/this->variation_time_scale))
                                 / energy_in_shell;
         }else{
             temp_famplitude = this->injection_rate / energy_in_shell;
@@ -440,7 +436,7 @@ void vorticity_equation<rnumber, be>::impose_forcing(
                             onew->cval(cindex, 2, 0)*onew->cval(cindex, 2, 0) + onew->cval(cindex, 2, 1)*onew->cval(cindex, 2, 1)
                             ) / k2;
                 *local_total_energy.getMine() += mode_energy;
-                double knorm = sqrt(k2);
+                double knorm = std::sqrt(k2);
                 if ((this->fk0 <= knorm) && (this->fk1 >= knorm))
                     *local_energy_in_shell.getMine() += mode_energy;
             }
@@ -467,14 +463,16 @@ void vorticity_equation<rnumber, be>::impose_forcing(
         energy_in_shell /= 2;
         // now, add forcing term
         // see Michael's thesis, page 38
-        double temp_famplitude = sqrt((this->energy - total_energy + energy_in_shell) / energy_in_shell);
+        double temp_famplitude = std::sqrt(
+                (this->energy - total_energy + energy_in_shell)
+               / energy_in_shell);
         this->kk->CLOOP_K2(
                     [&](const ptrdiff_t cindex,
                         const ptrdiff_t xindex,
                         const ptrdiff_t yindex,
                         const ptrdiff_t zindex,
                         const double k2){
-            const double knorm = sqrt(k2);
+            const double knorm = std::sqrt(k2);
             if ((this->fk0 <= knorm) &&
                 (this->fk1 >= knorm))
                 for (int c=0; c<3; c++)
@@ -811,7 +809,7 @@ void vorticity_equation<rnumber, be>::compute_Eulerian_acceleration(
                         - this->nu*k2*this->cvelocity->get_cdata()[tindex+cc][i];
             if (this->forcing_type == std::string("linear"))
             {
-                double knorm = sqrt(k2);
+                double knorm = std::sqrt(k2);
                 if ((this->fk0 <= knorm) &&
                         (this->fk1 >= knorm))
                 {
diff --git a/cpp/vorticity_equation.hpp b/cpp/vorticity_equation.hpp
index e790ae75c06e14979753affb6d114c2acd54caa3..b6ec1a9e03aba151b3cdd452da0ed1b3c230a460 100644
--- a/cpp/vorticity_equation.hpp
+++ b/cpp/vorticity_equation.hpp
@@ -29,15 +29,6 @@
 
 #include "field.hpp"
 
-#include <sys/stat.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <iostream>
-#include <string>
-
-
-extern int myrank, nprocs;
-
 
 /** container for field descriptor, fields themselves, parameters, etc
  *