diff --git a/bfps/cpp/distributed_particles.cpp b/bfps/cpp/distributed_particles.cpp
index bc19ea9417ca0ee943fb3db41928d07e66c2c04a..f486568160bdbd7ba6b8b9d7d0c66bc32ae678c9 100644
--- a/bfps/cpp/distributed_particles.cpp
+++ b/bfps/cpp/distributed_particles.cpp
@@ -39,7 +39,7 @@
 
 extern int myrank, nprocs;
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 distributed_particles<particle_type, rnumber, interp_neighbours>::distributed_particles(
         const char *NAME,
         const hid_t data_file_id,
@@ -61,12 +61,12 @@ distributed_particles<particle_type, rnumber, interp_neighbours>::distributed_pa
         this->rhs[i].reserve(2*this->nparticles / this->nprocs);
 }
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 distributed_particles<particle_type, rnumber, interp_neighbours>::~distributed_particles()
 {
 }
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void distributed_particles<particle_type, rnumber, interp_neighbours>::sample(
         interpolator<rnumber, interp_neighbours> *field,
         const std::unordered_map<int, single_particle_state<particle_type>> &x,
@@ -82,7 +82,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::sample(
     delete[] yy;
 }
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void distributed_particles<particle_type, rnumber, interp_neighbours>::get_rhs(
         const std::unordered_map<int, single_particle_state<particle_type>> &x,
         std::unordered_map<int, single_particle_state<particle_type>> &y)
@@ -99,7 +99,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::get_rhs(
     }
 }
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void distributed_particles<particle_type, rnumber, interp_neighbours>::sample(
         interpolator<rnumber, interp_neighbours> *field,
         const char *dset_name)
@@ -109,14 +109,14 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::sample(
     this->write(dset_name, y);
 }
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void distributed_particles<particle_type, rnumber, interp_neighbours>::roll_rhs()
 {
     for (int i=this->integration_steps-2; i>=0; i--)
         rhs[i+1] = rhs[i];
 }
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void distributed_particles<particle_type, rnumber, interp_neighbours>::redistribute(
         std::unordered_map<int, single_particle_state<particle_type>> &x,
         std::vector<std::unordered_map<int, single_particle_state<particle_type>>> &vals)
@@ -261,7 +261,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::redistrib
 
 
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void distributed_particles<particle_type, rnumber, interp_neighbours>::AdamsBashforth(
         const int nsteps)
 {
@@ -309,7 +309,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::AdamsBash
 }
 
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void distributed_particles<particle_type, rnumber, interp_neighbours>::step()
 {
     this->AdamsBashforth((this->iteration < this->integration_steps) ?
@@ -319,7 +319,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::step()
 }
 
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void distributed_particles<particle_type, rnumber, interp_neighbours>::read()
 {
     double *temp = new double[this->chunk_size*this->ncomponents];
@@ -363,7 +363,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::read()
     delete[] temp;
 }
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void distributed_particles<particle_type, rnumber, interp_neighbours>::write(
         const char *dset_name,
         std::unordered_map<int, single_particle_state<POINT3D>> &y)
@@ -395,7 +395,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::write(
     delete[] data;
 }
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void distributed_particles<particle_type, rnumber, interp_neighbours>::write(
         const bool write_rhs)
 {
diff --git a/bfps/cpp/distributed_particles.hpp b/bfps/cpp/distributed_particles.hpp
index 63a7ce35d68f7c01866008883c9443d27bc28241..cf6e124a7744c049b6fcf0c84c1618a0a214c30e 100644
--- a/bfps/cpp/distributed_particles.hpp
+++ b/bfps/cpp/distributed_particles.hpp
@@ -39,7 +39,7 @@
 
 #define DISTRIBUTED_PARTICLES
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 class distributed_particles: public particles_io_base<particle_type>
 {
     private:
diff --git a/bfps/cpp/particles.cpp b/bfps/cpp/particles.cpp
index 2182e769191d587aea2f5b174993ccae407bc109..9f55f65db12d05f12df9d7daa452c1d79f9b6fcf 100644
--- a/bfps/cpp/particles.cpp
+++ b/bfps/cpp/particles.cpp
@@ -39,7 +39,7 @@
 
 extern int myrank, nprocs;
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 particles<particle_type, rnumber, interp_neighbours>::particles(
         const char *NAME,
         const hid_t data_file_id,
@@ -65,7 +65,7 @@ particles<particle_type, rnumber, interp_neighbours>::particles(
     }
 }
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 particles<particle_type, rnumber, interp_neighbours>::~particles()
 {
     delete[] this->state;
@@ -75,7 +75,7 @@ particles<particle_type, rnumber, interp_neighbours>::~particles()
     }
 }
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void particles<particle_type, rnumber, interp_neighbours>::get_rhs(double *x, double *y)
 {
     switch(particle_type)
@@ -86,7 +86,7 @@ void particles<particle_type, rnumber, interp_neighbours>::get_rhs(double *x, do
     }
 }
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void particles<particle_type, rnumber, interp_neighbours>::roll_rhs()
 {
     for (int i=this->integration_steps-2; i>=0; i--)
@@ -97,7 +97,7 @@ void particles<particle_type, rnumber, interp_neighbours>::roll_rhs()
 
 
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void particles<particle_type, rnumber, interp_neighbours>::AdamsBashforth(
         const int nsteps)
 {
@@ -173,7 +173,7 @@ void particles<particle_type, rnumber, interp_neighbours>::AdamsBashforth(
 }
 
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void particles<particle_type, rnumber, interp_neighbours>::step()
 {
     this->AdamsBashforth((this->iteration < this->integration_steps) ?
@@ -183,7 +183,7 @@ void particles<particle_type, rnumber, interp_neighbours>::step()
 }
 
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void particles<particle_type, rnumber, interp_neighbours>::read()
 {
     if (this->myrank == 0)
@@ -210,7 +210,7 @@ void particles<particle_type, rnumber, interp_neighbours>::read()
                     this->comm);
 }
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void particles<particle_type, rnumber, interp_neighbours>::write(
         const bool write_rhs)
 {
@@ -224,7 +224,7 @@ void particles<particle_type, rnumber, interp_neighbours>::write(
         }
 }
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void particles<particle_type, rnumber, interp_neighbours>::sample(
         interpolator_base<rnumber, interp_neighbours> *field,
         const char *dset_name)
diff --git a/bfps/cpp/particles.hpp b/bfps/cpp/particles.hpp
index cf67b1fc32509b16d7a0774c0fc40a74759e39e0..06caefdd5e28881bf9a8e2dd80730645d5294583 100644
--- a/bfps/cpp/particles.hpp
+++ b/bfps/cpp/particles.hpp
@@ -37,7 +37,7 @@
 
 #define PARTICLES
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 class particles: public particles_io_base<particle_type>
 {
     private:
diff --git a/bfps/cpp/particles_base.cpp b/bfps/cpp/particles_base.cpp
index c0fc631799cf55e1fbd9b8c93db945f47aa9e980..5fd9e79c1dd535f552b88555c5836bd3b6f19b87 100644
--- a/bfps/cpp/particles_base.cpp
+++ b/bfps/cpp/particles_base.cpp
@@ -30,78 +30,56 @@
 #include <cassert>
 #include "particles_base.hpp"
 
-template <int particle_type>
+template <particle_types particle_type>
 single_particle_state<particle_type>::single_particle_state()
 {
-    switch(particle_type)
-    {
-        default:
-            this->data = new double[3];
-            std::fill_n(this->data, 3, 0);
-            break;
-    }
+    std::fill_n(this->data, state_dimension(particle_type), 0);
 }
 
-template <int particle_type>
+template <particle_types particle_type>
 single_particle_state<particle_type>::single_particle_state(
         const single_particle_state<particle_type> &src)
 {
-    switch(particle_type)
-    {
-        default:
-            this->data = new double[3];
-            std::copy(src.data, src.data + 3, this->data);
-            break;
-    }
+    std::copy(
+            src.data,
+            src.data + state_dimension(particle_type),
+            this->data);
 }
 
-template <int particle_type>
+template <particle_types particle_type>
 single_particle_state<particle_type>::single_particle_state(
         const double *src)
 {
-    switch(particle_type)
-    {
-        default:
-            this->data = new double[3];
-            std::copy(src, src + 3, this->data);
-            break;
-    }
+    std::copy(
+            src,
+            src + state_dimension(particle_type),
+            this->data);
 }
 
-template <int particle_type>
+template <particle_types particle_type>
 single_particle_state<particle_type>::~single_particle_state()
 {
-    switch(particle_type)
-    {
-        default:
-            delete[] this->data;
-            break;
-    }
 }
 
-template <int particle_type>
+template <particle_types particle_type>
 single_particle_state<particle_type> &single_particle_state<particle_type>::operator=(
         const single_particle_state &src)
 {
-    switch(particle_type)
-    {
-        default:
-            std::copy(src.data, src.data + 3, this->data);
-            break;
-    }
+    std::copy(
+            src.data,
+            src.data + state_dimension(particle_type),
+            this->data);
     return *this;
 }
 
-template <int particle_type>
+template <particle_types particle_type>
 single_particle_state<particle_type> &single_particle_state<particle_type>::operator=(
         const double *src)
 {
-    switch(particle_type)
-    {
-        default:
-            std::copy(src, src + 3, this->data);
-            break;
-    }
+    std::copy(
+            src,
+            src + state_dimension(particle_type),
+            this->data);
     return *this;
 }
 
@@ -134,19 +112,14 @@ int get_chunk_offsets(
     return EXIT_SUCCESS;
 }
 
-template <int particle_type>
+template <particle_types particle_type>
 particles_io_base<particle_type>::particles_io_base(
         const char *NAME,
         const int TRAJ_SKIP,
         const hid_t data_file_id,
         MPI_Comm COMM)
 {
-    switch(particle_type)
-    {
-        default:
-            this->ncomponents = 3;
-            break;
-    }
+    this->ncomponents = state_dimension(particle_type);
     this->name = std::string(NAME);
     this->traj_skip = TRAJ_SKIP;
     this->comm = COMM;
@@ -235,14 +208,14 @@ particles_io_base<particle_type>::particles_io_base(
     DEBUG_MSG("exiting particles_io_base constructor\n");
 }
 
-template <int particle_type>
+template <particle_types particle_type>
 particles_io_base<particle_type>::~particles_io_base()
 {
     if(this->myrank == 0)
         H5Gclose(this->hdf5_group_id);
 }
 
-template <int particle_type>
+template <particle_types particle_type>
 void particles_io_base<particle_type>::read_state_chunk(
         const int cindex,
         double *data)
@@ -276,7 +249,7 @@ void particles_io_base<particle_type>::read_state_chunk(
     DEBUG_MSG("exiting read_state_chunk\n");
 }
 
-template <int particle_type>
+template <particle_types particle_type>
 void particles_io_base<particle_type>::write_state_chunk(
         const int cindex,
         const double *data)
@@ -308,7 +281,7 @@ void particles_io_base<particle_type>::write_state_chunk(
     delete[] offset;
 }
 
-template <int particle_type>
+template <particle_types particle_type>
 void particles_io_base<particle_type>::read_rhs_chunk(
         const int cindex,
         const int rhsindex,
@@ -350,7 +323,7 @@ void particles_io_base<particle_type>::read_rhs_chunk(
     //DEBUG_MSG("exiting read_rhs_chunk\n");
 }
 
-template <int particle_type>
+template <particle_types particle_type>
 void particles_io_base<particle_type>::write_rhs_chunk(
         const int cindex,
         const int rhsindex,
@@ -387,7 +360,7 @@ void particles_io_base<particle_type>::write_rhs_chunk(
     delete[] offset;
 }
 
-template <int particle_type>
+template <particle_types particle_type>
 void particles_io_base<particle_type>::write_point3D_chunk(
         const std::string dset_name,
         const int cindex,
diff --git a/bfps/cpp/particles_base.hpp b/bfps/cpp/particles_base.hpp
index 923c29253ba1ddabff1791ca40d03f79109eaf20..a207e6b41c0670a7bf1db2946610ef9945bf5ea3 100644
--- a/bfps/cpp/particles_base.hpp
+++ b/bfps/cpp/particles_base.hpp
@@ -26,6 +26,7 @@
 
 #include <vector>
 #include <hdf5.h>
+#include <unordered_map>
 #include "interpolator_base.hpp"
 
 #ifndef PARTICLES_BASE
@@ -35,13 +36,21 @@
 /* particle types */
 enum particle_types {POINT3D, VELOCITY_TRACER};
 
+/* space dimension */
+constexpr unsigned int state_dimension(particle_types particle_type)
+{
+    return ((particle_type == POINT3D) ? 3 : (
+            (particle_type == VELOCITY_TRACER) ? 3 :
+            3));
+}
+
 /* 1 particle state type */
 
-template <int particle_type>
+template <particle_types particle_type>
 class single_particle_state
 {
     public:
-        double *data;
+        double data[state_dimension(particle_type)];
 
         single_particle_state();
         single_particle_state(const single_particle_state &src);
@@ -61,7 +70,7 @@ std::vector<std::vector<hsize_t>> get_chunk_offsets(
         std::vector<hsize_t> data_dims,
         std::vector<hsize_t> chnk_dims);
 
-template <int particle_type>
+template <particle_types particle_type>
 class particles_io_base
 {
     protected:
diff --git a/bfps/cpp/rFFTW_distributed_particles.cpp b/bfps/cpp/rFFTW_distributed_particles.cpp
index 0e982b491d9440cce27b96122066097bf31a0889..32aefaf3b7f40e1931d2ff5e66e7bb7d61655289 100644
--- a/bfps/cpp/rFFTW_distributed_particles.cpp
+++ b/bfps/cpp/rFFTW_distributed_particles.cpp
@@ -40,7 +40,7 @@
 
 extern int myrank, nprocs;
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::rFFTW_distributed_particles(
         const char *NAME,
         const hid_t data_file_id,
@@ -146,12 +146,12 @@ rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::rFFTW_di
     //}
 }
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::~rFFTW_distributed_particles()
 {
 }
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::sample(
         rFFTW_interpolator<rnumber, interp_neighbours> *field,
         const std::unordered_map<int, single_particle_state<particle_type>> &x,
@@ -213,7 +213,7 @@ void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::sam
     }
 }
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::get_rhs(
         const std::unordered_map<int, single_particle_state<particle_type>> &x,
         const std::unordered_map<int, std::unordered_set<int>> &dp,
@@ -231,7 +231,7 @@ void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::get
     }
 }
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::sample(
         rFFTW_interpolator<rnumber, interp_neighbours> *field,
         const char *dset_name)
@@ -241,14 +241,14 @@ void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::sam
     this->write(dset_name, y);
 }
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::roll_rhs()
 {
     for (int i=this->integration_steps-2; i>=0; i--)
         rhs[i+1] = rhs[i];
 }
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::redistribute(
         std::unordered_map<int, single_particle_state<particle_type>> &x,
         std::vector<std::unordered_map<int, single_particle_state<particle_type>>> &vals,
@@ -425,7 +425,7 @@ void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::red
 
 
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::AdamsBashforth(
         const int nsteps)
 {
@@ -473,7 +473,7 @@ void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::Ada
 }
 
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::step()
 {
     this->AdamsBashforth((this->iteration < this->integration_steps) ?
@@ -483,7 +483,7 @@ void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::ste
 }
 
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::sort_into_domains(
         const std::unordered_map<int, single_particle_state<particle_type>> &x,
         std::unordered_map<int, std::unordered_set<int>> &dp)
@@ -523,7 +523,7 @@ void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::sor
 }
 
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::read()
 {
     double *temp = new double[this->chunk_size*this->ncomponents];
@@ -575,7 +575,7 @@ void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::rea
     delete[] temp;
 }
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::write(
         const char *dset_name,
         std::unordered_map<int, single_particle_state<POINT3D>> &y)
@@ -611,7 +611,7 @@ void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::wri
     delete[] data;
 }
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::write(
         const bool write_rhs)
 {
diff --git a/bfps/cpp/rFFTW_distributed_particles.hpp b/bfps/cpp/rFFTW_distributed_particles.hpp
index cffb726baa35d3ab5020bf602ad9016bdde86371..e271bbfae56c0d49bf66cebcb5e8e8158f81940b 100644
--- a/bfps/cpp/rFFTW_distributed_particles.hpp
+++ b/bfps/cpp/rFFTW_distributed_particles.hpp
@@ -40,7 +40,7 @@
 
 #define RFFTW_DISTRIBUTED_PARTICLES
 
-template <int particle_type, class rnumber, int interp_neighbours>
+template <particle_types particle_type, class rnumber, int interp_neighbours>
 class rFFTW_distributed_particles: public particles_io_base<particle_type>
 {
     private: