diff --git a/README.md b/README.md
index 14e93f3204e46076b3245f4675bbd8a3f302c789..502f39c4a1decab2d668a7d53267316f3b1fe86f 100644
--- a/README.md
+++ b/README.md
@@ -8,6 +8,10 @@ This is a submodule of SWIFT containing only the reader.
 
 For details on the coverage, please see https://gitlab.cosma.dur.ac.uk/jenkins/job/CSDS-Pipeline/.
 
+To compile the code, run `./autogen.sh` and then `./configure; make -j 10`.
+To ensure its quality, the code is compiled by default with `-Wall`.
+If you need to disable it, configure the code with `./configure --enable-compiler-warnings=yes`.
+
 
 # File Description
 
diff --git a/src/csds_cosmology.cpp b/src/csds_cosmology.cpp
index 9a3f7845d8916ee1d8bf0999bd35b36f0f8186b8..c6da174eb373ab1b1bed12b233e4bb84aa0e287d 100644
--- a/src/csds_cosmology.cpp
+++ b/src/csds_cosmology.cpp
@@ -25,45 +25,38 @@
 #include "csds_tools.hpp"
 
 /**
- * @brief Compute the first time derivative of the scale factor
- *
- * @param cosmo The cosmological model
- * @param scale The scale factor
- *
- * @return The first derivative of the scale factor.
+ * @brief Initialize the #Cosmology from a file.
  */
-void csds_cosmology_init(struct csds_cosmology *cosmo,
-                         struct csds_parser *params) {
+Cosmology::Cosmology(struct csds_parser &params) {
 
   /* Read Omega_cdm */
-  cosmo->Omega_cdm = parser_get_param_double(params, "Cosmology:Omega_cdm");
+  mOmegaCdm = parser_get_param_double(params, "Cosmology:Omega_cdm");
 
   /* Read Omega_lambda */
-  cosmo->Omega_lambda =
-      parser_get_param_double(params, "Cosmology:Omega_lambda");
+  mOmegaLambda = parser_get_param_double(params, "Cosmology:Omega_lambda");
 
   /* Read Omega_b */
-  cosmo->Omega_b = parser_get_param_double(params, "Cosmology:Omega_b");
+  mOmegaBaryon = parser_get_param_double(params, "Cosmology:Omega_b");
 
   /* Read Omega_r */
-  cosmo->Omega_r = parser_get_param_double(params, "Cosmology:Omega_r");
+  mOmegaRadiation = parser_get_param_double(params, "Cosmology:Omega_r");
 
   /* Read Omega_k */
-  cosmo->Omega_k = parser_get_param_double(params, "Cosmology:Omega_k");
+  mOmegaCurvature = parser_get_param_double(params, "Cosmology:Omega_k");
 
   /* Read Omega_nu_0 */
-  double Omega_nu_0 =
+  double omega_nu_0 =
       parser_get_opt_param_double(params, "Cosmology:Omega_nu_0", 0);
 
-  if (Omega_nu_0 != 0)
+  if (omega_nu_0 != 0)
     csds_error("Neutrinos are not implemented in the cosmology.");
 
   /* Read w_0 */
-  cosmo->w_0 = parser_get_param_double(params, "Cosmology:w_0");
+  mW0 = parser_get_param_double(params, "Cosmology:w_0");
 
   /* Read w_a */
-  cosmo->w_a = parser_get_param_double(params, "Cosmology:w_a");
+  mWa = parser_get_param_double(params, "Cosmology:w_a");
 
   /* Read the Hubble constant at z=0 */
-  cosmo->H0 = parser_get_param_double(params, "Cosmology:Hubble0");
+  mH0 = parser_get_param_double(params, "Cosmology:Hubble0");
 }
diff --git a/src/csds_cosmology.hpp b/src/csds_cosmology.hpp
index b865aa620c3289924544ebae409589d78b3bb3be..f7c3adf5753cf446bca1e2cf4a6629a8988c7f30 100644
--- a/src/csds_cosmology.hpp
+++ b/src/csds_cosmology.hpp
@@ -36,158 +36,145 @@ struct csds_parser;
 /**
  * @brief Structure containing the cosmological parameters.
  */
-struct csds_cosmology {
+class Cosmology {
+
+ public:
+  /**
+   * @brief Compute E(a) for the derivatives of the scale factor.
+   * @param scale The scale factor
+   * @return The first derivative of the scale factor.
+   */
+  template <class Type>
+  INLINE Type GetE(const Type scale) const {
+
+    const Type a_inv = 1. / scale;
+    const Type a_inv2 = a_inv * a_inv;
+    const Type omega_matter = mOmegaBaryon + mOmegaCdm;
+    const Type dark_eos = GetWTilde<Type>(scale);
+    const Type E = sqrt(mOmegaRadiation * a_inv2 * a_inv2 + /* Radiation */
+                        omega_matter * a_inv * a_inv2 +     /* Matter */
+                        mOmegaCurvature * a_inv2 +          /* Curvature */
+                        mOmegaLambda * dark_eos);
+    return E;
+  }
+
+  /**
+   * @brief Compute the first time derivative of the scale factor
+   * @param scale The scale factor
+   * @return The first derivative of the scale factor.
+   */
+  template <class Type>
+  INLINE Type GetScaleFactorDerivative(const Type scale) const {
+    const Type E = GetE<Type>(scale);
+    return mH0 * scale * E;
+  }
+
+  /**
+   * @brief Compute the first time derivative of the scale factor
+   * @param scale The scale factor
+   * @return The first derivative of the scale factor.
+   */
+  template <class Type>
+  INLINE Type GetScaleFactorSecondDerivative(const Type scale) const {
+
+    /* Compute some powers of a */
+    const Type a_inv = 1. / scale;
+    const Type a_inv2 = a_inv * a_inv;
+    const Type a_inv4 = a_inv2 * a_inv2;
+
+    /* Compute E */
+    const Type E = GetE<Type>(scale);
+
+    /* Compute the derivative of E */
+    Type dark_eos_der = (mWa) - (1 + mW0 + mWa) * a_inv;
+    dark_eos_der *= GetWTilde<Type>(scale);
+    const Type omega_matter = mOmegaBaryon + mOmegaCdm;
+    const Type E_der =
+        (-4 * mOmegaRadiation * a_inv4 * a_inv - 3 * omega_matter * a_inv4 -
+         2 * mOmegaCurvature * a_inv2 * a_inv + dark_eos_der);
+
+    /* Compute the full expression */
+    const Type a_dot = GetScaleFactorDerivative<Type>(scale);
+    return mH0 * (a_dot * E + 0.5 * scale * E_der / E);
+  }
+
+  /**
+   * @brief Add the cosmological factors for the position.
+   * See Hausammann et al. 2021 for details.
+   *
+   * @param scale The scale factor.
+   * @param vel The velocity to modify.
+   * @param acc The acceleration to modify.
+   */
+  template <class TypeVel, class TypeAcc>
+  INLINE void AddFactorPosition(const TypeVel scale, TypeVel *vel,
+                                TypeAcc *acc) const {
+
+    const TypeVel a_dot = GetScaleFactorDerivative<TypeVel>(scale);
+    const TypeVel scale2 = scale * scale;
+    if (acc) {
+      const TypeAcc a_dot_inv = 1. / a_dot;
+      const TypeAcc scale2_inv = 1. / scale2;
+      const TypeAcc a_ddot = GetScaleFactorSecondDerivative<TypeAcc>(scale);
+      *acc = (*acc - 2 * a_dot * *vel) * scale;
+      *acc -= scale2 * a_ddot * *vel * a_dot_inv;
+      *acc *= scale2_inv * scale2_inv * a_dot_inv * a_dot_inv;
+    }
+    *vel /= scale2 * a_dot;
+  }
+
+  /**
+   * @brief Add the cosmological factors for the velocities.
+   * \f$ \frac{dv'}{da} = \frac{g'}{\dot{a}a} \f$.
+   *
+   * @param scale The scale factor.
+   * @param acc The acceleration to modify.
+   */
+  template <class TypeAcc>
+  INLINE void AddFactorVelocity(const TypeAcc scale, TypeAcc *acc) const {
+
+    const TypeAcc a_dot = GetScaleFactorDerivative<TypeAcc>(scale);
+    *acc /= scale * a_dot;
+  }
+
+  Cosmology(struct csds_parser &params);
+  Cosmology() {}
+
+ protected:
+  /**
+   * @brief Compute the dark energy equation of state.
+   * @param scale The scale factor
+   * @return The equation of state.
+   */
+  template <class Type>
+  INLINE Type GetWTilde(const Type scale) const {
+    const Type w_tilde = (scale - 1.) * mWa - (1. + mW0 + mWa) * log(scale);
+    return exp(3. * w_tilde);
+  }
 
   /*! Cold dark matter density parameter */
-  double Omega_cdm;
+  double mOmegaCdm;
 
   /*! Baryon density parameter */
-  double Omega_b;
+  double mOmegaBaryon;
 
   /*! Cosmological constant density parameter */
-  double Omega_lambda;
+  double mOmegaLambda;
 
   /*! Total radiation density parameter (photons and other relics) */
-  double Omega_r;
+  double mOmegaRadiation;
 
   /*! Curvature density parameter */
-  double Omega_k;
+  double mOmegaCurvature;
 
   /*! Dark-energy equation of state at z=0 */
-  double w_0;
+  double mW0;
 
   /*! Dark-energy evolution parameter */
-  double w_a;
+  double mWa;
 
   /*! Hubble constant at z = 0 (in internal units) */
-  double H0;
+  double mH0;
 };
 
-/**
- * @brief Compute the dark energy equation of state.
- *
- * @param cosmo The cosmological model
- * @param scale The scale factor
- *
- * @return The equation of state.
- */
-INLINE static double csds_cosmology_get_w_tilde(
-    const struct csds_cosmology *cosmo, const double scale) {
-  const double w_tilde =
-      (scale - 1.) * cosmo->w_a - (1. + cosmo->w_0 + cosmo->w_a) * log(scale);
-  return exp(3. * w_tilde);
-}
-
-/**
- * @brief Compute E(a) for the derivatives of the scale factor.
- *
- * @param cosmo The cosmological model
- * @param scale The scale factor
- *
- * @return The first derivative of the scale factor.
- */
-INLINE static double csds_cosmology_get_E(const struct csds_cosmology *cosmo,
-                                          const double scale) {
-
-  const double a_inv = 1. / scale;
-  const double a_inv2 = a_inv * a_inv;
-  const double Omega_m = cosmo->Omega_b + cosmo->Omega_cdm;
-  const double dark_eos = csds_cosmology_get_w_tilde(cosmo, scale);
-  const double E = sqrt(cosmo->Omega_r * a_inv2 * a_inv2 + /* Radiation */
-                        Omega_m * a_inv * a_inv2 +         /* Matter */
-                        cosmo->Omega_k * a_inv2 +          /* Curvature */
-                        cosmo->Omega_lambda * dark_eos);
-  return E;
-}
-
-/**
- * @brief Compute the first time derivative of the scale factor
- *
- * @param cosmo The cosmological model
- * @param scale The scale factor
- *
- * @return The first derivative of the scale factor.
- */
-INLINE static double csds_cosmology_get_scale_factor_derivative(
-    const struct csds_cosmology *cosmo, const double scale) {
-  const double E = csds_cosmology_get_E(cosmo, scale);
-  return cosmo->H0 * scale * E;
-}
-
-/**
- * @brief Compute the first time derivative of the scale factor
- *
- * @param cosmo The cosmological model
- * @param scale The scale factor
- *
- * @return The first derivative of the scale factor.
- */
-INLINE static double csds_cosmology_get_scale_factor_second_derivative(
-    const struct csds_cosmology *cosmo, const double scale) {
-
-  /* Compute some powers of a */
-  const double a_inv = 1. / scale;
-  const double a_inv2 = a_inv * a_inv;
-  const double a_inv4 = a_inv2 * a_inv2;
-
-  /* Compute E */
-  const double E = csds_cosmology_get_E(cosmo, scale);
-
-  /* Compute the derivative of E */
-  double dark_eos_der = (cosmo->w_a) - (1 + cosmo->w_0 + cosmo->w_a) * a_inv;
-  dark_eos_der *= csds_cosmology_get_w_tilde(cosmo, scale);
-  const double Omega_m = cosmo->Omega_b + cosmo->Omega_cdm;
-  const double E_der =
-      (-4 * cosmo->Omega_r * a_inv4 * a_inv - 3 * Omega_m * a_inv4 -
-       2 * cosmo->Omega_k * a_inv2 * a_inv + dark_eos_der);
-
-  /* Compute the full expression */
-  const double a_dot = csds_cosmology_get_scale_factor_derivative(cosmo, scale);
-  return cosmo->H0 * (a_dot * E + 0.5 * scale * E_der / E);
-}
-
-/**
- * @brief Add the cosmological factors for the position.
- * See Hausammann et al. 2021 for details.
- *
- * @param cosmo The cosmological model.
- * @param scale The scale factor.
- * @param vel The velocity to modify.
- * @param acc The acceleration to modify.
- */
-INLINE static void csds_cosmology_add_factor_position(
-    const struct csds_cosmology *cosmo, const double scale, float *vel,
-    float *acc) {
-
-  const double a_dot = csds_cosmology_get_scale_factor_derivative(cosmo, scale);
-  const double scale2 = scale * scale;
-  if (acc) {
-    const double a_dot_inv = 1. / a_dot;
-    const double scale2_inv = 1. / scale2;
-    const double a_ddot =
-        csds_cosmology_get_scale_factor_second_derivative(cosmo, scale);
-    *acc = (*acc - 2 * a_dot * *vel) * scale;
-    *acc -= scale2 * a_ddot * *vel * a_dot_inv;
-    *acc *= scale2_inv * scale2_inv * a_dot_inv * a_dot_inv;
-  }
-  *vel /= scale2 * a_dot;
-}
-
-/**
- * @brief Add the cosmological factors for the velocities.
- * \f$ \frac{dv'}{da} = \frac{g'}{\dot{a}a} \f$.
- *
- * @param cosmo The cosmological model.
- * @param scale The scale factor.
- * @param acc The acceleration to modify.
- */
-INLINE static void csds_cosmology_add_factor_velocity(
-    const struct csds_cosmology *cosmo, const double scale, float *acc) {
-
-  const double a_dot = csds_cosmology_get_scale_factor_derivative(cosmo, scale);
-  *acc /= scale * a_dot;
-}
-
-void csds_cosmology_init(struct csds_cosmology *cosmo,
-                         struct csds_parser *params);
 #endif  // CSDS_CSDS_PARAMETERS_H
diff --git a/src/csds_interpolation.hpp b/src/csds_interpolation.hpp
index 9c807c8c19ce7a80e5cf50c3606e22c51f96301d..d9dd74a5e3107d30bf91c3dee37b750cb29702d3 100644
--- a/src/csds_interpolation.hpp
+++ b/src/csds_interpolation.hpp
@@ -188,14 +188,16 @@ INLINE static TypeField interpolate_linear_spline(const TypeField t0,
  * @param t_after Time of field_after (> t).
  * @param t Requested time.
  * @param periodic Should the periodic boundary be applied?
- * @param params The simulation's #csds_parameters.
+ * @param params The simulation's #Parameters.
  */
 template <class TypeField, class TypeFirst, class TypeSecond>
-INLINE static void interpolate_quintic(
-    const TypeField t_before, const struct csds_reader_field &before,
-    const TypeField t_after, const struct csds_reader_field &after,
-    void *__restrict__ output, const TypeField t, const int dimension,
-    int periodic, const struct csds_parameters &params) {
+INLINE static void interpolate_quintic(const TypeField t_before,
+                                       const struct csds_reader_field &before,
+                                       const TypeField t_after,
+                                       const struct csds_reader_field &after,
+                                       void *__restrict__ output,
+                                       const TypeField t, const int dimension,
+                                       int periodic, const Parameters &params) {
 
   /* Get the arrays */
   const TypeField *x_bef = (TypeField *)before.field;
@@ -219,21 +221,21 @@ INLINE static void interpolate_quintic(
     TypeSecond a1 = a_aft ? a_aft[i] : -1;
 
     /* Add cosmological factor for velocity and acceleration */
-    if (params.with_cosmology && v_bef && v_aft) {
-      csds_cosmology_add_factor_position(&params.cosmo, t_before, &v0,
-                                         a_bef ? &a0 : NULL);
-      csds_cosmology_add_factor_position(&params.cosmo, t_after, &v1,
-                                         a_aft ? &a1 : NULL);
+    if (params.IsWithCosmology() && v_bef && v_aft) {
+      params.GetCosmology().AddFactorPosition<TypeFirst, TypeSecond>(
+          t_before, &v0, a_bef ? &a0 : NULL);
+      params.GetCosmology().AddFactorPosition<TypeFirst, TypeSecond>(
+          t_after, &v1, a_aft ? &a1 : NULL);
     }
 
     /* Apply the periodic boundaries */
     if (periodic) {
       /* Move towards the origin for more accuracy */
-      const TypeField half = 0.5 * params.box_size[i];
+      const TypeField half = 0.5 * params.GetBoxSize()[i];
       if (x0 > x1 + half)
-        x0 -= params.box_size[i];
+        x0 -= params.GetBoxSize()[i];
       else if (x1 > x0 + half)
-        x1 -= params.box_size[i];
+        x1 -= params.GetBoxSize()[i];
     }
 
     /* Use quintic hermite spline. */
@@ -258,13 +260,13 @@ INLINE static void interpolate_quintic(
     if (periodic) {
       /* We suppose that only 1 wrapping is enough
          otherwise something really bad happened */
-      x[i] = box_wrap(x[i], 0., params.box_size[i]);
+      x[i] = box_wrap(x[i], 0., params.GetBoxSize()[i]);
 
 #ifdef CSDS_DEBUG_CHECKS
       /* Check if the periodic box is correctly applied */
-      if (x[i] >= params.box_size[i] || x[i] < 0) {
+      if (x[i] >= params.GetBoxSize()[i] || x[i] < 0) {
         csds_error("A particle is outside the periodic box "
-                   << params.box_size[i] << ": before=" << x0
+                   << params.GetBoxSize()[i] << ": before=" << x0
                    << ", after=" << x1 << ", interpolation=" << x[i]);
       }
 #endif
@@ -292,7 +294,7 @@ INLINE static void interpolate_cubic(
     const TypeField t_before, const struct csds_reader_field &before,
     const TypeField t_after, const struct csds_reader_field &after,
     void *__restrict__ output, const TypeField t, const int dimension,
-    ATTR_UNUSED int periodic, const struct csds_parameters &params) {
+    ATTR_UNUSED int periodic, const Parameters &params) {
 
 #ifdef CSDS_DEBUG_CHECKS
   /* The position is expected to be a double => error with periodic */
@@ -316,9 +318,9 @@ INLINE static void interpolate_cubic(
     TypeFirst a1 = a_aft ? a_aft[i] : -1;
 
     /* Add cosmological factor for acceleration */
-    if (params.with_cosmology && a_bef && a_aft) {
-      csds_cosmology_add_factor_velocity(&params.cosmo, t_before, &a0);
-      csds_cosmology_add_factor_velocity(&params.cosmo, t_after, &a0);
+    if (params.IsWithCosmology() && a_bef && a_aft) {
+      params.GetCosmology().AddFactorVelocity<TypeFirst>(t_before, &a0);
+      params.GetCosmology().AddFactorVelocity<TypeFirst>(t_after, &a0);
     }
 
     /* Use a cubic hermite spline. */
@@ -352,8 +354,7 @@ INLINE static void interpolate_linear(
     const TypeField t_before, const struct csds_reader_field &before,
     const TypeField t_after, const struct csds_reader_field &after,
     void *__restrict__ output, const TypeField t, const int dimension,
-    ATTR_UNUSED int periodic,
-    ATTR_UNUSED const struct csds_parameters &params) {
+    ATTR_UNUSED int periodic, ATTR_UNUSED const Parameters &params) {
 
 #ifdef CSDS_DEBUG_CHECKS
   /* The position is expected to be a double => error with periodic */
diff --git a/src/csds_parameters.cpp b/src/csds_parameters.cpp
index 749f871ab3af911ddfcc92456b65c78466d36244..ab317e2c7404208c79ae56a506015c332d3c32fc 100644
--- a/src/csds_parameters.cpp
+++ b/src/csds_parameters.cpp
@@ -30,60 +30,61 @@ using namespace std;
 /**
  * @brief Initialize the parameters from the yaml file.
  *
- * @param reader The #csds_reader.
  * @param filename The filename to use.
  * @param verbose the verbose level
  *
  */
-void csds_parameters_init(struct csds_parameters *params, const string filename,
-                          int verbose) {
+Parameters::Parameters(const string filename, int verbose) {
 
   /* Initialize the parser */
   struct csds_parser csds_parser;
-  parser_read_file(filename, &csds_parser);
+  parser_read_file(filename, csds_parser);
 
   /* Print the parameters */
-  if (verbose > 1) parser_print_params(&csds_parser);
+  if (verbose > 1) parser_print_params(csds_parser);
 
   /* Read the number of dimension */
-  params->dimension = parser_get_param_int(&csds_parser, "Header:Dimension");
+  mDimension = parser_get_param_int(csds_parser, "Header:Dimension");
 
   /* Read the box size */
-  parser_get_param_double_array(&csds_parser, "Header:BoxSize", 3,
-                                params->box_size);
+  double box_size[3];
+  parser_get_param_double_array(csds_parser, "Header:BoxSize", 3, box_size);
+  for (int i = 0; i < 3; i++) {
+    mBoxSize[i] = box_size[i];
+  }
 
   /* Read the periodicity */
-  params->periodic = parser_get_param_int(&csds_parser, "Header:Periodic");
+  mPeriodic = parser_get_param_int(csds_parser, "Header:Periodic");
 
   /* Read if we are running with cosmology */
-  params->with_cosmology =
-      parser_get_param_int(&csds_parser, "Policy:cosmological integration");
+  mWithCosmology =
+      parser_get_param_int(csds_parser, "Policy:cosmological integration");
 
   /* Initialize the cosmology */
-  if (params->with_cosmology) csds_cosmology_init(&params->cosmo, &csds_parser);
+  if (mWithCosmology) mCosmo = Cosmology(csds_parser);
 
   /* Read the initial number of particles. */
   for (int i = 0; i < csds_type_count; i++) {
-    params->approximate_number_particles[i] = 0;
+    mApproximateNumberParticles[i] = 0;
   }
-  params->approximate_number_particles[csds_type_gas] =
-      parser_get_param_longlong(&csds_parser, "Header:NumberParts");
-  params->approximate_number_particles[csds_type_stars] =
-      parser_get_param_longlong(&csds_parser, "Header:NumberSParts");
-  params->approximate_number_particles[csds_type_dark_matter] =
-      parser_get_param_longlong(&csds_parser, "Header:NumberGParts");
+  mApproximateNumberParticles[csds_type_gas] =
+      parser_get_param_longlong(csds_parser, "Header:NumberParts");
+  mApproximateNumberParticles[csds_type_stars] =
+      parser_get_param_longlong(csds_parser, "Header:NumberSParts");
+  mApproximateNumberParticles[csds_type_dark_matter] =
+      parser_get_param_longlong(csds_parser, "Header:NumberGParts");
 
   /* Read the cache parameters */
-  params->cache_over_alloc =
-      parser_get_opt_param_float(&csds_parser, "Cache:OverAllocation", 1.2);
+  mCacheOverAllocation =
+      parser_get_opt_param_float(csds_parser, "Cache:OverAllocation", 1.2);
 
   long long init_alloc[csds_type_count];
   for (int i = 0; i < csds_type_count; i++) {
     init_alloc[i] = 0;
   }
-  parser_get_opt_param_longlong_array(&csds_parser, "Cache:InitialAllocation",
+  parser_get_opt_param_longlong_array(csds_parser, "Cache:InitialAllocation",
                                       csds_type_count, init_alloc);
   for (int i = 0; i < csds_type_count; i++) {
-    params->cache_init_size[i] = init_alloc[i];
+    mCacheInitSize[i] = init_alloc[i];
   }
 }
diff --git a/src/csds_parameters.hpp b/src/csds_parameters.hpp
index 7fe4acfd0bb6cb0ec6d11298eeac8254ccad0c14..345a3d5c15875857b7dc8e2a0095d44a600f5782 100644
--- a/src/csds_parameters.hpp
+++ b/src/csds_parameters.hpp
@@ -30,42 +30,57 @@
 #include "csds_part_type.hpp"
 
 /* Standard headers */
+#include <array>
 #include <stdint.h>
 #include <string>
 
-class Reader;
-
 /**
  * @brief Structure containing the simulation's parameters.
  */
-struct csds_parameters {
-
+class Parameters {
+
+ public:
+  Parameters() {}
+  Parameters(const std::string filename, int verbose);
+
+  int GetDimension() const { return mDimension; }
+  std::array<double, 3> const& GetBoxSize() const { return mBoxSize; }
+  bool IsPeriodic() const { return mPeriodic; }
+  bool IsWithCosmology() const { return mWithCosmology; }
+  std::array<long long, csds_type_count> const& GetApproximateNumberParticles()
+      const {
+    return mApproximateNumberParticles;
+  }
+  Cosmology const& GetCosmology() const { return mCosmo; }
+  float GetCacheOverAllocation() const { return mCacheOverAllocation; }
+  std::array<int64_t, csds_type_count> const& GetCacheInitialSize() const {
+    return mCacheInitSize;
+  }
+
+ protected:
   /* Number of dimension */
-  int dimension;
+  int mDimension;
 
   /* Simulation volume */
-  double box_size[3];
+  std::array<double, 3> mBoxSize;
 
   /* Is the volume periodic? */
-  int periodic;
+  bool mPeriodic;
 
   /* Are we doing a cosmological simulation? */
-  int with_cosmology;
+  bool mWithCosmology;
 
   /* Number of particles initially present. */
-  long long approximate_number_particles[csds_type_count];
+  std::array<long long, csds_type_count> mApproximateNumberParticles;
 
   /* The cosmological model */
-  struct csds_cosmology cosmo;
+  Cosmology mCosmo;
 
   /* Over allocation for the cache. */
-  float cache_over_alloc;
+  float mCacheOverAllocation;
 
   /* Initial cache for in-existing particles */
-  int64_t cache_init_size[csds_type_count];
+  std::array<int64_t, csds_type_count> mCacheInitSize;
 };
 
-void csds_parameters_init(struct csds_parameters *params,
-                          const std::string filename, int verbose);
-
 #endif  // CSDS_PARAMETERS_H
diff --git a/src/csds_parser.cpp b/src/csds_parser.cpp
index aa5c546a85cdb13378f35c9052dc0b3abccca526..5d4795611bf202d7d32cdca792ad2d9c86e9d367 100644
--- a/src/csds_parser.cpp
+++ b/src/csds_parser.cpp
@@ -49,13 +49,13 @@ typedef long long longlong;
 /* Private functions. */
 static int is_empty(const char *str);
 static int count_indentation(const char *str);
-static void parse_line(char *line, struct csds_parser *params);
-static void parse_value(char *line, struct csds_parser *params);
+static void parse_line(char *line, struct csds_parser &params);
+static void parse_value(char *line, struct csds_parser &params);
 static void parse_section_param(char *line, int *isFirstParam,
-                                char *sectionName, struct csds_parser *params);
-static void find_duplicate_params(const struct csds_parser *params,
+                                char *sectionName, struct csds_parser &params);
+static void find_duplicate_params(const struct csds_parser &params,
                                   const char *param_name);
-static void find_duplicate_section(const struct csds_parser *params,
+static void find_duplicate_section(const struct csds_parser &params,
                                    const char *section_name);
 static int lineNumber = 0;
 
@@ -108,10 +108,10 @@ char *trim_both(char *s) {
  * @param file_name Name of file to be read
  * @param params Structure to be populated from file
  */
-void parser_init(const std::string filename, struct csds_parser *params) {
-  params->paramCount = 0;
-  params->sectionCount = 0;
-  params->fileName = filename;
+void parser_init(const std::string filename, struct csds_parser &params) {
+  params.paramCount = 0;
+  params.sectionCount = 0;
+  params.fileName = filename;
 }
 
 /**
@@ -120,7 +120,7 @@ void parser_init(const std::string filename, struct csds_parser *params) {
  * @param file_name Name of file to be read
  * @param params Structure to be populated from file
  */
-void parser_read_file(const std::string file_name, struct csds_parser *params) {
+void parser_read_file(const std::string file_name, struct csds_parser &params) {
   /* Open file for reading */
   FILE *file = fopen(file_name.data(), "r");
 
@@ -192,10 +192,10 @@ static int is_empty(const char *str) {
  * @param param_name Name of parameter to be searched for
  */
 
-static void find_duplicate_params(const struct csds_parser *params,
+static void find_duplicate_params(const struct csds_parser &params,
                                   const char *param_name) {
-  for (int i = 0; i < params->paramCount; i++) {
-    if (!strcmp(param_name, params->data[i].name)) {
+  for (int i = 0; i < params.paramCount; i++) {
+    if (!strcmp(param_name, params.data[i].name)) {
       csds_error("Invalid line:" << lineNumber << " '" << param_name
                                  << "', parameter is a duplicate.");
     }
@@ -209,10 +209,10 @@ static void find_duplicate_params(const struct csds_parser *params,
  * @param section_name Name of section to be searched for
  */
 
-static void find_duplicate_section(const struct csds_parser *params,
+static void find_duplicate_section(const struct csds_parser &params,
                                    const char *section_name) {
-  for (int i = 0; i < params->sectionCount; i++) {
-    if (!strcmp(section_name, params->section[i].name)) {
+  for (int i = 0; i < params.sectionCount; i++) {
+    if (!strcmp(section_name, params.section[i].name)) {
       csds_error("Invalid line:" << lineNumber << " '" << section_name
                                  << "', parameter is a duplicate.");
     }
@@ -225,7 +225,7 @@ static void find_duplicate_section(const struct csds_parser *params,
  * @param params Structure to be populated from file.
  */
 
-static void parse_line(char *line, struct csds_parser *params) {
+static void parse_line(char *line, struct csds_parser &params) {
   /* Parse line if it doesn't begin with a comment. */
   if (line[0] != PARSER_COMMENT_CHAR) {
     char trim_line[PARSER_MAX_LINE_SIZE];
@@ -267,7 +267,7 @@ static void parse_line(char *line, struct csds_parser *params) {
  * @param params Structure to be written to
  *
  */
-static void parse_value(char *line, struct csds_parser *params) {
+static void parse_value(char *line, struct csds_parser &params) {
   static int inSection = 0;
   static char section[PARSER_MAX_LINE_SIZE]; /* Keeps track of current section
                                                 name. */
@@ -310,12 +310,12 @@ static void parse_value(char *line, struct csds_parser *params) {
       find_duplicate_params(params, tmpStr);
 
       strcpy(section, tmpSectionName);
-      strcpy(params->section[params->sectionCount].name, tmpSectionName);
-      if (params->sectionCount == PARSER_MAX_NO_OF_SECTIONS - 1) {
+      strcpy(params.section[params.sectionCount].name, tmpSectionName);
+      if (params.sectionCount == PARSER_MAX_NO_OF_SECTIONS - 1) {
         csds_error(
             "Maximal number of sections in parameter file reached. Aborting !");
       } else {
-        params->sectionCount++;
+        params.sectionCount++;
       }
       inSection = 1;
       isFirstParam = 1;
@@ -333,14 +333,14 @@ static void parse_value(char *line, struct csds_parser *params) {
 
       /* Must be a standalone parameter so no need to prefix name with a
        * section. */
-      strcpy(params->data[params->paramCount].name, tmpStr);
-      strcpy(params->data[params->paramCount].value, token);
-      if (params->paramCount == PARSER_MAX_NO_OF_PARAMS - 1) {
+      strcpy(params.data[params.paramCount].name, tmpStr);
+      strcpy(params.data[params.paramCount].value, token);
+      if (params.paramCount == PARSER_MAX_NO_OF_PARAMS - 1) {
         csds_error(
             "Maximal number of parameters in parameter file reached. Aborting "
             "!");
       } else {
-        params->paramCount++;
+        params.paramCount++;
       }
       inSection = 0;
       isFirstParam = 1;
@@ -360,7 +360,7 @@ static void parse_value(char *line, struct csds_parser *params) {
  */
 
 static void parse_section_param(char *line, int *isFirstParam,
-                                char *sectionName, struct csds_parser *params) {
+                                char *sectionName, struct csds_parser &params) {
   static int sectionIndent = 0;
   char tmpStr[PARSER_MAX_LINE_SIZE];
   char paramName[PARSER_MAX_LINE_SIZE];
@@ -391,40 +391,40 @@ static void parse_section_param(char *line, int *isFirstParam,
   /* Check for duplicate parameter name. */
   find_duplicate_params(params, paramName);
 
-  strcpy(params->data[params->paramCount].name, paramName);
-  strcpy(params->data[params->paramCount].value, token);
-  if (params->paramCount == PARSER_MAX_NO_OF_PARAMS - 1) {
+  strcpy(params.data[params.paramCount].name, paramName);
+  strcpy(params.data[params.paramCount].value, token);
+  if (params.paramCount == PARSER_MAX_NO_OF_PARAMS - 1) {
     csds_error(
         "Maximal number of parameters in parameter file reached. Aborting !");
   } else {
-    params->paramCount++;
+    params.paramCount++;
   }
 }
 
 // Retrieve parameter value from structure. TYPE is the data type, float, int
 // etc. FMT the format required for that data type, i.e. %f, %d etc. and DESC
 // a one word description of the type, "float", "int" etc.
-#define PARSER_GET_VALUE(TYPE, FMT, DESC)                                     \
-  static int get_param_##TYPE(struct csds_parser *params, const char *name,   \
-                              TYPE *def, TYPE *result) {                      \
-    char str[PARSER_MAX_LINE_SIZE];                                           \
-    for (int i = 0; i < params->paramCount; i++) {                            \
-      if (strcmp(name, params->data[i].name) == 0) {                          \
-        /* Check that exactly one number is parsed, capture junk. */          \
-        if (sscanf(params->data[i].value, " " FMT "%s ", result, str) != 1) { \
-          csds_error("Tried parsing "                                         \
-                     << DESC << " '" << params->data[i].name                  \
-                     << "' but found '" << params->data[i].value              \
-                     << "' with illegal trailing characters '" << str         \
-                     << "'.");                                                \
-        }                                                                     \
-        return 1;                                                             \
-      }                                                                       \
-    }                                                                         \
-    if (def == NULL)                                                          \
-      csds_error("Cannot find '" << name << "' in the structure, in file: "   \
-                                 << params->fileName);                        \
-    return 0;                                                                 \
+#define PARSER_GET_VALUE(TYPE, FMT, DESC)                                      \
+  static int get_param_##TYPE(struct csds_parser &params, const char *name,    \
+                              TYPE *def, TYPE *result) {                       \
+    char str[PARSER_MAX_LINE_SIZE];                                            \
+    for (int i = 0; i < params.paramCount; i++) {                              \
+      if (strcmp(name, params.data[i].name) == 0) {                            \
+        /* Check that exactly one number is parsed, capture junk. */           \
+        if (sscanf(params.data[i].value, " " FMT "%s ", result, str) != 1) {   \
+          csds_error("Tried parsing "                                          \
+                     << DESC << " '" << params.data[i].name << "' but found '" \
+                     << params.data[i].value                                   \
+                     << "' with illegal trailing characters '" << str          \
+                     << "'.");                                                 \
+        }                                                                      \
+        return 1;                                                              \
+      }                                                                        \
+    }                                                                          \
+    if (def == NULL)                                                           \
+      csds_error("Cannot find '" << name << "' in the structure, in file: "    \
+                                 << params.fileName);                          \
+    return 0;                                                                  \
   }
 
 /* Instantiations. */
@@ -441,7 +441,7 @@ PARSER_GET_VALUE(longlong, "%lld", "long long");
  * @param name Name of the parameter to be found
  * @return Value of the parameter found
  */
-int parser_get_param_int(struct csds_parser *params, const char *name) {
+int parser_get_param_int(struct csds_parser &params, const char *name) {
   int result = 0;
   get_param_int(params, name, NULL, &result);
   return result;
@@ -454,7 +454,7 @@ int parser_get_param_int(struct csds_parser *params, const char *name) {
  * @param name Name of the parameter to be found
  * @return Value of the parameter found
  */
-char parser_get_param_char(struct csds_parser *params, const char *name) {
+char parser_get_param_char(struct csds_parser &params, const char *name) {
   char result = 0;
   get_param_char(params, name, NULL, &result);
   return result;
@@ -467,7 +467,7 @@ char parser_get_param_char(struct csds_parser *params, const char *name) {
  * @param name Name of the parameter to be found
  * @return Value of the parameter found
  */
-float parser_get_param_float(struct csds_parser *params, const char *name) {
+float parser_get_param_float(struct csds_parser &params, const char *name) {
   float result = 0;
   get_param_float(params, name, NULL, &result);
   return result;
@@ -480,7 +480,7 @@ float parser_get_param_float(struct csds_parser *params, const char *name) {
  * @param name Name of the parameter to be found
  * @return Value of the parameter found
  */
-double parser_get_param_double(struct csds_parser *params, const char *name) {
+double parser_get_param_double(struct csds_parser &params, const char *name) {
   double result = 0;
   get_param_double(params, name, NULL, &result);
   return result;
@@ -493,7 +493,7 @@ double parser_get_param_double(struct csds_parser *params, const char *name) {
  * @param name Name of the parameter to be found
  * @return Value of the parameter found
  */
-long long parser_get_param_longlong(struct csds_parser *params,
+long long parser_get_param_longlong(struct csds_parser &params,
                                     const char *name) {
   long long result = 0;
   get_param_longlong(params, name, NULL, &result);
@@ -507,12 +507,12 @@ long long parser_get_param_longlong(struct csds_parser *params,
  * @param name Name of the parameter to be found
  * @param retParam (return) Value of the parameter found
  */
-void parser_get_param_string(struct csds_parser *params, const char *name,
+void parser_get_param_string(struct csds_parser &params, const char *name,
                              char *retParam) {
 
-  for (int i = 0; i < params->paramCount; i++) {
-    if (!strcmp(name, params->data[i].name)) {
-      strcpy(retParam, params->data[i].value);
+  for (int i = 0; i < params.paramCount; i++) {
+    if (!strcmp(name, params.data[i].name)) {
+      strcpy(retParam, params.data[i].value);
       return;
     }
   }
@@ -528,7 +528,7 @@ void parser_get_param_string(struct csds_parser *params, const char *name,
  * @param def Default value of the parameter of not found.
  * @return Value of the parameter found
  */
-int parser_get_opt_param_int(struct csds_parser *params, const char *name,
+int parser_get_opt_param_int(struct csds_parser &params, const char *name,
                              int def) {
   int result = 0;
   if (get_param_int(params, name, &def, &result)) return result;
@@ -543,7 +543,7 @@ int parser_get_opt_param_int(struct csds_parser *params, const char *name,
  * @param def Default value of the parameter of not found.
  * @return Value of the parameter found
  */
-char parser_get_opt_param_char(struct csds_parser *params, const char *name,
+char parser_get_opt_param_char(struct csds_parser &params, const char *name,
                                char def) {
   char result = 0;
   if (get_param_char(params, name, &def, &result)) return result;
@@ -558,7 +558,7 @@ char parser_get_opt_param_char(struct csds_parser *params, const char *name,
  * @param def Default value of the parameter of not found.
  * @return Value of the parameter found
  */
-float parser_get_opt_param_float(struct csds_parser *params, const char *name,
+float parser_get_opt_param_float(struct csds_parser &params, const char *name,
                                  float def) {
   float result = 0;
   if (get_param_float(params, name, &def, &result)) return result;
@@ -573,7 +573,7 @@ float parser_get_opt_param_float(struct csds_parser *params, const char *name,
  * @param def Default value of the parameter of not found.
  * @return Value of the parameter found
  */
-double parser_get_opt_param_double(struct csds_parser *params, const char *name,
+double parser_get_opt_param_double(struct csds_parser &params, const char *name,
                                    double def) {
   double result = 0;
   if (get_param_double(params, name, &def, &result)) return result;
@@ -588,7 +588,7 @@ double parser_get_opt_param_double(struct csds_parser *params, const char *name,
  * @param def Default value of the parameter of not found.
  * @return Value of the parameter found
  */
-long long parser_get_opt_param_longlong(struct csds_parser *params,
+long long parser_get_opt_param_longlong(struct csds_parser &params,
                                         const char *name, long long def) {
   long long result = 0;
   if (get_param_longlong(params, name, &def, &result)) return result;
@@ -603,12 +603,12 @@ long long parser_get_opt_param_longlong(struct csds_parser *params,
  * @param def Default value of the parameter of not found.
  * @param retParam (return) Value of the parameter found
  */
-void parser_get_opt_param_string(struct csds_parser *params, const char *name,
+void parser_get_opt_param_string(struct csds_parser &params, const char *name,
                                  char *retParam, const char *def) {
 
-  for (int i = 0; i < params->paramCount; i++) {
-    if (!strcmp(name, params->data[i].name)) {
-      strcpy(retParam, params->data[i].value);
+  for (int i = 0; i < params.paramCount; i++) {
+    if (!strcmp(name, params.data[i].name)) {
+      strcpy(retParam, params.data[i].value);
       return;
     }
   }
@@ -622,16 +622,16 @@ void parser_get_opt_param_string(struct csds_parser *params, const char *name,
  * i.e. "float".
  */
 #define PARSER_GET_ARRAY(TYPE, FMT, DESC)                                   \
-  static int get_param_##TYPE##_array(struct csds_parser *params,           \
+  static int get_param_##TYPE##_array(struct csds_parser &params,           \
                                       const char *name, int required,       \
                                       int nval, TYPE *values) {             \
     char str[PARSER_MAX_LINE_SIZE];                                         \
     char cpy[PARSER_MAX_LINE_SIZE];                                         \
                                                                             \
-    for (int i = 0; i < params->paramCount; i++) {                          \
-      if (!strcmp(name, params->data[i].name)) {                            \
+    for (int i = 0; i < params.paramCount; i++) {                           \
+      if (!strcmp(name, params.data[i].name)) {                             \
         char *cp = cpy;                                                     \
-        strcpy(cp, params->data[i].value);                                  \
+        strcpy(cp, params.data[i].value);                                   \
         cp = trim_both(cp);                                                 \
                                                                             \
         /* Strip off [], if present. */                                     \
@@ -652,15 +652,15 @@ void parser_get_opt_param_string(struct csds_parser *params, const char *name,
             TYPE tmp_value;                                                 \
             if (sscanf(p, fmt, &tmp_value, str) != 1) {                     \
               csds_error("Tried parsing "                                   \
-                         << DESC << " '" << params->data[i].name            \
-                         << "' but found '" << params->data[i].value        \
+                         << DESC << " '" << params.data[i].name             \
+                         << "' but found '" << params.data[i].value         \
                          << "' with illegal trailing characters '" << str   \
                          << "'.");                                          \
             }                                                               \
             values[k] = tmp_value;                                          \
           } else {                                                          \
             csds_error("Array '"                                            \
-                       << name << "' with value '" << params->data[i].value \
+                       << name << "' with value '" << params.data[i].value  \
                        << "' has too few values, expected " << nval);       \
           }                                                                 \
           if (k < nval - 1) p = strtok(NULL, ",");                          \
@@ -670,7 +670,7 @@ void parser_get_opt_param_string(struct csds_parser *params, const char *name,
     }                                                                       \
     if (required)                                                           \
       csds_error("Cannot find '" << name << "' in the structure, in file '" \
-                                 << params->fileName << "'.");              \
+                                 << params.fileName << "'.");               \
     return 0;                                                               \
   }
 
@@ -689,7 +689,7 @@ PARSER_GET_ARRAY(longlong, "%lld", "long long");
  * @param nval number of values expected.
  * @param values Values of the parameter found, of size at least nvals.
  */
-void parser_get_param_char_array(struct csds_parser *params, const char *name,
+void parser_get_param_char_array(struct csds_parser &params, const char *name,
                                  int nval, char *values) {
   get_param_char_array(params, name, 1, nval, values);
 }
@@ -705,7 +705,7 @@ void parser_get_param_char_array(struct csds_parser *params, const char *name,
  *               unmodified, so should be set to the default values.
  * @return whether the parameter has been found.
  */
-int parser_get_opt_param_char_array(struct csds_parser *params,
+int parser_get_opt_param_char_array(struct csds_parser &params,
                                     const char *name, int nval, char *values) {
   if (get_param_char_array(params, name, 0, nval, values) != 1) {
     return 0;
@@ -721,7 +721,7 @@ int parser_get_opt_param_char_array(struct csds_parser *params,
  * @param nval number of values expected.
  * @param values Values of the parameter found, of size at least nvals.
  */
-void parser_get_param_int_array(struct csds_parser *params, const char *name,
+void parser_get_param_int_array(struct csds_parser &params, const char *name,
                                 int nval, int *values) {
   get_param_int_array(params, name, 1, nval, values);
 }
@@ -737,7 +737,7 @@ void parser_get_param_int_array(struct csds_parser *params, const char *name,
  *               unmodified, so should be set to the default values.
  * @return whether the parameter has been found.
  */
-int parser_get_opt_param_int_array(struct csds_parser *params, const char *name,
+int parser_get_opt_param_int_array(struct csds_parser &params, const char *name,
                                    int nval, int *values) {
   if (get_param_int_array(params, name, 0, nval, values) != 1) {
     return 0;
@@ -753,7 +753,7 @@ int parser_get_opt_param_int_array(struct csds_parser *params, const char *name,
  * @param nval number of values expected.
  * @param values Values of the parameter found, of size at least nvals.
  */
-void parser_get_param_float_array(struct csds_parser *params, const char *name,
+void parser_get_param_float_array(struct csds_parser &params, const char *name,
                                   int nval, float *values) {
   get_param_float_array(params, name, 1, nval, values);
 }
@@ -769,7 +769,7 @@ void parser_get_param_float_array(struct csds_parser *params, const char *name,
  *               unmodified, so should be set to the default values.
  * @return whether the parameter has been found.
  */
-int parser_get_opt_param_float_array(struct csds_parser *params,
+int parser_get_opt_param_float_array(struct csds_parser &params,
                                      const char *name, int nval,
                                      float *values) {
   if (get_param_float_array(params, name, 0, nval, values) != 1) {
@@ -786,7 +786,7 @@ int parser_get_opt_param_float_array(struct csds_parser *params,
  * @param nval number of values expected.
  * @param values Values of the parameter found, of size at least nvals.
  */
-void parser_get_param_double_array(struct csds_parser *params, const char *name,
+void parser_get_param_double_array(struct csds_parser &params, const char *name,
                                    int nval, double *values) {
   get_param_double_array(params, name, 1, nval, values);
 }
@@ -802,7 +802,7 @@ void parser_get_param_double_array(struct csds_parser *params, const char *name,
  *               unmodified, so should be set to the default values.
  * @return whether the parameter has been found.
  */
-int parser_get_opt_param_double_array(struct csds_parser *params,
+int parser_get_opt_param_double_array(struct csds_parser &params,
                                       const char *name, int nval,
                                       double *values) {
   if (get_param_double_array(params, name, 0, nval, values) != 1) {
@@ -819,7 +819,7 @@ int parser_get_opt_param_double_array(struct csds_parser *params,
  * @param nval number of values expected.
  * @param values Values of the parameter found, of size at least nvals.
  */
-void parser_get_param_longlong_array(struct csds_parser *params,
+void parser_get_param_longlong_array(struct csds_parser &params,
                                      const char *name, int nval,
                                      long long *values) {
   get_param_longlong_array(params, name, 1, nval, values);
@@ -836,7 +836,7 @@ void parser_get_param_longlong_array(struct csds_parser *params,
  *               unmodified, so should be set to the default values.
  * @return whether the parameter has been found.
  */
-int parser_get_opt_param_longlong_array(struct csds_parser *params,
+int parser_get_opt_param_longlong_array(struct csds_parser &params,
                                         const char *name, int nval,
                                         long long *values) {
   if (get_param_longlong_array(params, name, 0, nval, values) != 1) {
@@ -850,13 +850,13 @@ int parser_get_opt_param_longlong_array(struct csds_parser *params,
  *
  * @param params Structure that holds the parameters
  */
-void parser_print_params(const struct csds_parser *params) {
+void parser_print_params(const struct csds_parser &params) {
   printf("\n--------------------------\n");
   printf("|   CSDS Parameter File  |\n");
   printf("--------------------------\n");
 
-  for (int i = 0; i < params->paramCount; i++) {
-    printf("Parameter name: %s\n", params->data[i].name);
-    printf("Parameter value: %s\n", params->data[i].value);
+  for (int i = 0; i < params.paramCount; i++) {
+    printf("Parameter name: %s\n", params.data[i].name);
+    printf("Parameter value: %s\n", params.data[i].value);
   }
 }
diff --git a/src/csds_parser.hpp b/src/csds_parser.hpp
index ada8a23a63170452bb86c46df6376583b64acf20..58a6ecb495655d6159f79d9d1666438c6866b57c 100644
--- a/src/csds_parser.hpp
+++ b/src/csds_parser.hpp
@@ -56,52 +56,52 @@ struct csds_parser {
 };
 
 /* Public API. */
-void parser_init(const std::string file_name, struct csds_parser *params);
-void parser_read_file(const std::string file_name, struct csds_parser *params);
-void parser_print_params(const struct csds_parser *params);
+void parser_init(const std::string file_name, struct csds_parser &params);
+void parser_read_file(const std::string file_name, struct csds_parser &params);
+void parser_print_params(const struct csds_parser &params);
 
-char parser_get_param_char(struct csds_parser *params, const char *name);
-int parser_get_param_int(struct csds_parser *params, const char *name);
-float parser_get_param_float(struct csds_parser *params, const char *name);
-double parser_get_param_double(struct csds_parser *params, const char *name);
-long long parser_get_param_longlong(struct csds_parser *params,
+char parser_get_param_char(struct csds_parser &params, const char *name);
+int parser_get_param_int(struct csds_parser &params, const char *name);
+float parser_get_param_float(struct csds_parser &params, const char *name);
+double parser_get_param_double(struct csds_parser &params, const char *name);
+long long parser_get_param_longlong(struct csds_parser &params,
                                     const char *name);
-void parser_get_param_string(struct csds_parser *params, const char *name,
+void parser_get_param_string(struct csds_parser &params, const char *name,
                              char *retParam);
 
-char parser_get_opt_param_char(struct csds_parser *params, const char *name,
+char parser_get_opt_param_char(struct csds_parser &params, const char *name,
                                char def);
-int parser_get_opt_param_int(struct csds_parser *params, const char *name,
+int parser_get_opt_param_int(struct csds_parser &params, const char *name,
                              int def);
-float parser_get_opt_param_float(struct csds_parser *params, const char *name,
+float parser_get_opt_param_float(struct csds_parser &params, const char *name,
                                  float def);
-double parser_get_opt_param_double(struct csds_parser *params, const char *name,
+double parser_get_opt_param_double(struct csds_parser &params, const char *name,
                                    double def);
-long long parser_get_opt_param_longlong(struct csds_parser *params,
+long long parser_get_opt_param_longlong(struct csds_parser &params,
                                         const char *name, long long def);
-void parser_get_opt_param_string(struct csds_parser *params, const char *name,
+void parser_get_opt_param_string(struct csds_parser &params, const char *name,
                                  char *retParam, const char *def);
-void parser_get_param_char_array(struct csds_parser *params, const char *name,
+void parser_get_param_char_array(struct csds_parser &params, const char *name,
                                  int nval, char *values);
-void parser_get_param_int_array(struct csds_parser *params, const char *name,
+void parser_get_param_int_array(struct csds_parser &params, const char *name,
                                 int nval, int *values);
-void parser_get_param_float_array(struct csds_parser *params, const char *name,
+void parser_get_param_float_array(struct csds_parser &params, const char *name,
                                   int nval, float *values);
-void parser_get_param_double_array(struct csds_parser *params, const char *name,
+void parser_get_param_double_array(struct csds_parser &params, const char *name,
                                    int nval, double *values);
-void parser_get_param_longlong_array(struct csds_parser *params,
+void parser_get_param_longlong_array(struct csds_parser &params,
                                      const char *name, int nval,
                                      long long *values);
-int parser_get_opt_param_char_array(struct csds_parser *params,
+int parser_get_opt_param_char_array(struct csds_parser &params,
                                     const char *name, int nval, char *values);
-int parser_get_opt_param_int_array(struct csds_parser *params, const char *name,
+int parser_get_opt_param_int_array(struct csds_parser &params, const char *name,
                                    int nval, int *values);
-int parser_get_opt_param_float_array(struct csds_parser *params,
+int parser_get_opt_param_float_array(struct csds_parser &params,
                                      const char *name, int nval, float *values);
-int parser_get_opt_param_double_array(struct csds_parser *params,
+int parser_get_opt_param_double_array(struct csds_parser &params,
                                       const char *name, int nval,
                                       double *values);
-int parser_get_opt_param_longlong_array(struct csds_parser *params,
+int parser_get_opt_param_longlong_array(struct csds_parser &params,
                                         const char *name, int nval,
                                         long long *values);
 
diff --git a/src/csds_particle.hpp b/src/csds_particle.hpp
index 68c84daf342fb931b7f6d82db95e82a76cd762fb..2be50411ea13cf346393226e38f0f8c4a9fee0af 100644
--- a/src/csds_particle.hpp
+++ b/src/csds_particle.hpp
@@ -167,13 +167,13 @@ INLINE static enum csds_special_flags csds_particle_read_special_flag(
  * @param time_after Time of field_after (> t).
  * @param time Requested time.
  * @param field The field to reconstruct.
- * @param params The simulation's #csds_parameters.
+ * @param params The simulation's #Parameters.
  */
 INLINE static void csds_particle_interpolate_field(
     const double time_before, const struct csds_reader_field &before,
     const double time_after, const struct csds_reader_field &after,
     void *__restrict__ output, const double time, const Field &field,
-    const struct csds_parameters &params) {
+    const Parameters &params) {
 
 #ifdef CSDS_DEBUG_CHECKS
   if (time_before > time || time_after < time) {
@@ -187,7 +187,7 @@ INLINE static void csds_particle_interpolate_field(
     case field_enum_coordinates:
       interpolate_quintic<double, float, float>(
           time_before, before, time_after, after, output, time,
-          /* dimension */ 3, params.periodic, params);
+          /* dimension */ 3, params.IsPeriodic(), params);
       break;
 
       /* Velocities */
diff --git a/src/csds_reader.cpp b/src/csds_reader.cpp
index b29de40d3cd04d2a22d5411ffe02030d981056a4..57810a1b7030f625bb226d0b8e3be26d430e2e81 100644
--- a/src/csds_reader.cpp
+++ b/src/csds_reader.cpp
@@ -63,7 +63,7 @@ Reader::Reader(const string basename, int verbose, int number_index,
   /* Remove the rank number */
   params_filename.erase(params_filename.size() - 5, params_filename.size());
   params_filename += ".yml";
-  csds_parameters_init(&this->mParams, params_filename, this->mVerbose);
+  mParams = Parameters(params_filename, this->mVerbose);
 
   /* Initialize the log file. */
   this->mLog = new LogFile(basename, /* only_header */ 0, this->mVerbose);
@@ -711,8 +711,8 @@ void Reader::ReadAllParticles(
   if (this->mUseCache) {
     for (int i = 0; i < csds_type_count; i++) {
       mCache[i] =
-          Cache(n_part[i], output, this->mParams.cache_over_alloc,
-                this->mParams.cache_init_size[i], this->mTime.time_offset);
+          Cache(n_part[i], output, mParams.GetCacheOverAllocation(),
+                mParams.GetCacheInitialSize()[i], this->mTime.time_offset);
     }
   }
 
@@ -814,8 +814,8 @@ void Reader::ReadParticlesFromIds(
   if (this->mUseCache) {
     for (int i = 0; i < csds_type_count; i++) {
       this->mCache[i] =
-          Cache(n_part[i], output, this->mParams.cache_over_alloc,
-                this->mParams.cache_init_size[i], this->mTime.time_offset);
+          Cache(n_part[i], output, mParams.GetCacheOverAllocation(),
+                mParams.GetCacheInitialSize()[i], this->mTime.time_offset);
     }
   }
 
diff --git a/src/csds_reader.hpp b/src/csds_reader.hpp
index c6efb84174e7b3d5da78d811f077a76db2e7edb1..15b6310fde2857e11404874e5b9a2740daa5abd8 100644
--- a/src/csds_reader.hpp
+++ b/src/csds_reader.hpp
@@ -128,7 +128,7 @@ class Reader {
   } mTime;
 
   /* Information from the yaml file */
-  struct csds_parameters mParams;
+  Parameters mParams;
 
   /* Level of verbosity. */
   int mVerbose;
diff --git a/src/csds_reader_generate_index.cpp b/src/csds_reader_generate_index.cpp
index 81a5a6968de9151626be3cf2953e698cce1168d1..448f786954aa4b0e89a1ab9d36bb55aaf3b3ecd1 100644
--- a/src/csds_reader_generate_index.cpp
+++ b/src/csds_reader_generate_index.cpp
@@ -674,7 +674,7 @@ void Reader::GenerateIndexFiles(int number_index, int current_index) {
 
     /* Allocate the arrays for the current state */
     for (int i = 0; i < csds_type_count; i++) {
-      current_state[i].reserve(this->mParams.approximate_number_particles[i] *
+      current_state[i].reserve(mParams.GetApproximateNumberParticles()[i] *
                                hashmap_overallocation);
     }
 
diff --git a/tests/testInterpolation.cpp b/tests/testInterpolation.cpp
index 79d002ee987ca4c4886f70f7a17f2ad4771c21a4..45887aee71b0f1a92b164bf06ade37fa8c53ad82 100644
--- a/tests/testInterpolation.cpp
+++ b/tests/testInterpolation.cpp
@@ -57,8 +57,7 @@ long long id_func(void) { return 1698241543; }
 
 int main(void) {
   /* Initialize the parameters */
-  struct csds_parameters params;
-  csds_parameters_init(&params, "test_interpolation.yml", 2);
+  Parameters params("test_interpolation.yml", 2);
 
   /* Construct the fields */
   const int n_fields = 4;
diff --git a/tests/testParser.cpp b/tests/testParser.cpp
index 6b8469e273f2c29711fc991d41e6c6fca2cef545..58429d3dd461f342f9bdcb1667e99d511180900a 100644
--- a/tests/testParser.cpp
+++ b/tests/testParser.cpp
@@ -29,42 +29,42 @@
 int main(void) {
 
   struct csds_parser parser;
-  parser_init("test_parser.yml", &parser);
-  parser_read_file("test_parser.yml", &parser);
-  parser_print_params(&parser);
+  parser_init("test_parser.yml", parser);
+  parser_read_file("test_parser.yml", parser);
+  parser_print_params(parser);
 
   /* scalar reader */
-  assert(parser_get_param_char(&parser, "Test:char") == 'c');
-  assert(parser_get_param_int(&parser, "Test:int") == 10);
-  float ftest = parser_get_param_float(&parser, "Test:float");
+  assert(parser_get_param_char(parser, "Test:char") == 'c');
+  assert(parser_get_param_int(parser, "Test:int") == 10);
+  float ftest = parser_get_param_float(parser, "Test:float");
   assert(fabsf(ftest - 0.2f) < 1e-6);
-  assert(parser_get_param_double(&parser, "Test:double") == 0.2);
-  assert(parser_get_param_longlong(&parser, "Test:longlong") == 1254108);
+  assert(parser_get_param_double(parser, "Test:double") == 0.2);
+  assert(parser_get_param_longlong(parser, "Test:longlong") == 1254108);
 
   char ret_param[CSDS_STRING_SIZE];
-  parser_get_param_string(&parser, "Test:string", ret_param);
+  parser_get_param_string(parser, "Test:string", ret_param);
   assert(strcmp(ret_param, "test") == 0);
 
   /* Optional scalar reader (failed) */
-  assert(parser_get_opt_param_char(&parser, "Test:test", 'd') == 'd');
-  assert(parser_get_opt_param_int(&parser, "Test:test", 2) == 2);
-  ftest = parser_get_opt_param_float(&parser, "Test:test", 0.1);
+  assert(parser_get_opt_param_char(parser, "Test:test", 'd') == 'd');
+  assert(parser_get_opt_param_int(parser, "Test:test", 2) == 2);
+  ftest = parser_get_opt_param_float(parser, "Test:test", 0.1);
   assert(fabsf(ftest - 0.1f) < 1e-6);
-  assert(parser_get_opt_param_double(&parser, "Test:test", 0.1) == 0.1);
-  assert(parser_get_opt_param_longlong(&parser, "Test:test", 1) == 1);
+  assert(parser_get_opt_param_double(parser, "Test:test", 0.1) == 0.1);
+  assert(parser_get_opt_param_longlong(parser, "Test:test", 1) == 1);
 
-  parser_get_opt_param_string(&parser, "Test:test", ret_param, "test1");
+  parser_get_opt_param_string(parser, "Test:test", ret_param, "test1");
   assert(strcmp(ret_param, "test1") == 0);
 
   /* Optional scalar reader (pass) */
-  assert(parser_get_opt_param_char(&parser, "Test:char", 'd') == 'c');
-  assert(parser_get_opt_param_int(&parser, "Test:int", 1) == 10);
-  ftest = parser_get_opt_param_float(&parser, "Test:float", 0.1);
+  assert(parser_get_opt_param_char(parser, "Test:char", 'd') == 'c');
+  assert(parser_get_opt_param_int(parser, "Test:int", 1) == 10);
+  ftest = parser_get_opt_param_float(parser, "Test:float", 0.1);
   assert(fabsf(ftest - 0.2f) < 1e-6);
-  assert(parser_get_opt_param_double(&parser, "Test:double", 0.1) == 0.2);
-  assert(parser_get_opt_param_longlong(&parser, "Test:longlong", 1) == 1254108);
+  assert(parser_get_opt_param_double(parser, "Test:double", 0.1) == 0.2);
+  assert(parser_get_opt_param_longlong(parser, "Test:longlong", 1) == 1254108);
 
-  parser_get_opt_param_string(&parser, "Test:string", ret_param, "test1");
+  parser_get_opt_param_string(parser, "Test:string", ret_param, "test1");
   assert(strcmp(ret_param, "test") == 0);
 
   /* Array reader */
@@ -75,11 +75,11 @@ int main(void) {
   double d[nval];
   long long l[nval];
 
-  parser_get_param_char_array(&parser, "Test:char_array", nval, c);
-  parser_get_param_int_array(&parser, "Test:int_array", nval, i);
-  parser_get_param_float_array(&parser, "Test:float_array", nval, f);
-  parser_get_param_double_array(&parser, "Test:double_array", nval, d);
-  parser_get_param_longlong_array(&parser, "Test:longlong_array", nval, l);
+  parser_get_param_char_array(parser, "Test:char_array", nval, c);
+  parser_get_param_int_array(parser, "Test:int_array", nval, i);
+  parser_get_param_float_array(parser, "Test:float_array", nval, f);
+  parser_get_param_double_array(parser, "Test:double_array", nval, d);
+  parser_get_param_longlong_array(parser, "Test:longlong_array", nval, l);
 
   /* Check the values */
   assert(c[0] == 'a');
@@ -115,11 +115,11 @@ int main(void) {
   l[1] = 4;
   l[2] = 5;
 
-  parser_get_opt_param_char_array(&parser, "Test:test", nval, c);
-  parser_get_opt_param_int_array(&parser, "Test:test", nval, i);
-  parser_get_opt_param_float_array(&parser, "Test:test", nval, f);
-  parser_get_opt_param_double_array(&parser, "Test:test", nval, d);
-  parser_get_opt_param_longlong_array(&parser, "Test:test", nval, l);
+  parser_get_opt_param_char_array(parser, "Test:test", nval, c);
+  parser_get_opt_param_int_array(parser, "Test:test", nval, i);
+  parser_get_opt_param_float_array(parser, "Test:test", nval, f);
+  parser_get_opt_param_double_array(parser, "Test:test", nval, d);
+  parser_get_opt_param_longlong_array(parser, "Test:test", nval, l);
 
   /* Check the values */
   assert(c[0] == 'e');
@@ -159,11 +159,11 @@ int main(void) {
   l[1] = 4;
   l[2] = 5;
 
-  parser_get_opt_param_char_array(&parser, "Test:char_array", nval, c);
-  parser_get_opt_param_int_array(&parser, "Test:int_array", nval, i);
-  parser_get_opt_param_float_array(&parser, "Test:float_array", nval, f);
-  parser_get_opt_param_double_array(&parser, "Test:double_array", nval, d);
-  parser_get_opt_param_longlong_array(&parser, "Test:longlong_array", nval, l);
+  parser_get_opt_param_char_array(parser, "Test:char_array", nval, c);
+  parser_get_opt_param_int_array(parser, "Test:int_array", nval, i);
+  parser_get_opt_param_float_array(parser, "Test:float_array", nval, f);
+  parser_get_opt_param_double_array(parser, "Test:double_array", nval, d);
+  parser_get_opt_param_longlong_array(parser, "Test:longlong_array", nval, l);
 
   /* Check the values */
   assert(c[0] == 'a');