include/Particle.hpp

@@ -14,14 +14,16 @@
 
 #include <armadillo>
 
+#include "typedefs.hpp"
+
 /** @brief A class that holds attributes of a particle
  * */
 class Particle {
 private:
-    double q; ///< Charge
+    double q;     ///< Charge
-    double m; ///< Mass
+    double m;     ///< Mass
-    arma::vec::fixed<3> r_vec; ///< position
+    vec_3d r_vec; ///< position
-    arma::vec::fixed<3> v_vec; ///< velocity
+    vec_3d v_vec; ///< velocity
 
 public:
     /** @brief Initialize the particle.
@@ -29,9 +31,7 @@ public:
      *  @details Initialize the particle with a charge, mass, position and
      *  velocity.
      * */
-    Particle(double q, double m, 
+    Particle(double q, double m, vec_3d r_vec, vec_3d v_vec);
-             arma::vec::fixed<3> r_vec, 
-             arma::vec::fixed<3> v_vec);
 
     /** @brief Make private attributes available for PenningTrap.
      * */

include/PenningTrap.hpp

@@ -17,12 +17,10 @@
 
 #include "Particle.hpp"
 #include "constants.hpp"
+#include "typedefs.hpp"
 
-typedef std::vector<arma::vec::fixed<3>> sim_cols;
-typedef std::vector<arma::vec::fixed<3>> sim_rows;
-typedef std::vector<sim_cols> sim_arr;
 
-#pragma omp declare reduction(+ : arma::vec : omp_out += omp_in)               \
+#pragma omp declare reduction(+ : vec_3d : omp_out += omp_in)               \
     initializer(omp_priv = omp_orig)
 
 /** @brief A class that simulates a Penning trap.
@@ -36,10 +34,10 @@ private:
     double V_0;                      ///< Applied potential
     double d;                        ///< Characteristic dimension
     std::vector<Particle> particles; ///< The particles in the Penning trap
-    arma::vec::fixed<3> *k_v;
+    sim_arr k_v;
-    arma::vec::fixed<3> *k_r;
+    sim_arr k_r;
-    std::function<arma::vec(int, double)> v_funcs[4];
+    std::function<vec_3d(int, double)> v_funcs[4];
-    std::function<arma::vec(int, double)> r_funcs[4];
+    std::function<vec_3d(int, double)> r_funcs[4];
 
 public:
     /** @brief Set B_0, V_0 and d.
@@ -58,33 +56,33 @@ public:
 
     /** @brief Calculate E at point r
      * */
-    arma::vec external_E_field(arma::vec r);
+    vec_3d external_E_field(vec_3d r);
 
     /** @brief Calculate B at point r
      * */
-    arma::vec external_B_field(arma::vec r);
+    vec_3d external_B_field(vec_3d r);
 
     /** @brief Calculate the force between 2 particles.
      *
      *  @details Calculate the force exhibited on particle p_i from
      *  particle p_j.
      * */
-    arma::vec force_on_particle(int i, int j);
+    vec_3d force_on_particle(int i, int j);
 
     /** @brief Calculate the total external force on a particle.
      *
      *  @details Calculate the total amount of force that E and B exhibits
      *  on particle p_i.
      * */
-    arma::vec total_force_external(int i);
+    vec_3d total_force_external(int i);
 
     /** @brief Calculate the total force on a particle from other particles.
      * */
-    arma::vec total_force_particles(int i);
+    vec_3d total_force_particles(unsigned int i);
 
     /** @brief calculate the total force on a particle.
      * */
-    arma::vec total_force(int i);
+    vec_3d total_force(int i);
 
     /** @brief Go forward one timestep using the RK4 method
      * */
@@ -94,8 +92,9 @@ public:
      * */
     void evolve_forward_euler(double dt, bool particle_interaction = true);
 
-    sim_arr simulate(double time, int steps, std::string method = "rk4",
+    sim_arr simulate(double time, unsigned int steps,
-                     bool particle_interaction = true);
+                                              std::string method = "rk4",
+                                              bool particle_interaction = true);
 
     void write_simulation_to_dir(std::string path, double time, int steps,
                                  std::string method = "rk4",

include/typedefs.hpp

@@ -0,0 +1,12 @@
+#ifndef __TYPEDEFS__
+#define __TYPEDEFS__
+
+#include <vector>
+#include <armadillo>
+
+typedef std::vector<arma::vec::fixed<3>> sim_cols;
+typedef std::vector<arma::vec::fixed<3>> sim_rows;
+typedef std::vector<sim_cols> sim_arr;
+typedef arma::vec::fixed<3> vec_3d;
+
+#endif

src/PenningTrap.cpp

@@ -15,6 +15,7 @@
 
 #include "PenningTrap.hpp"
 #include "constants.hpp"
+#include "typedefs.hpp"
 #include "utils.hpp"
 #include <cstdlib>
 #include <functional>
@@ -27,42 +28,31 @@ PenningTrap::PenningTrap(double B_0, double V_0, double d)
     this->V_0 = V_0;
     this->d = d;
 
-    v_funcs[0] = [this](int i, double dt) { return .5 * dt * this->k_v[i]; };
+    // Create lambda functions for each iteration in RK4
-    v_funcs[1] = [this](int i, double dt) {
+    v_funcs[0] = [this](int i, double dt) { return .5 * dt * this->k_v[0][i]; };
-        return .5 * dt * this->k_v[this->particles.size() + i];
+    v_funcs[1] = [this](int i, double dt) { return .5 * dt * this->k_v[1][i]; };
-    };
+    v_funcs[2] = [this](int i, double dt) { return dt * this->k_v[2][i]; };
-    v_funcs[2] = [this](int i, double dt) {
-        return dt * this->k_v[2 * this->particles.size() + i];
-    };
     v_funcs[3] = [this](int i, double dt) {
-        arma::vec res = this->k_v[i] + this->k_v[this->particles.size() + i] +
+        return (dt / 6) * (this->k_v[0][i].eval() + this->k_v[1][i].eval() +
-                        this->k_v[2 * this->particles.size() + i] +
+                           this->k_v[2][i].eval() + this->k_v[3][i].eval());
-                        this->k_v[3 * this->particles.size() + i];
-        return (dt / 6) * res;
     };
 
-    r_funcs[0] = [this](int i, double dt) { return .5 * dt * this->k_r[i]; };
+    r_funcs[0] = [this](int i, double dt) { return .5 * dt * this->k_r[0][i]; };
-    r_funcs[1] = [this](int i, double dt) {
+    r_funcs[1] = [this](int i, double dt) { return .5 * dt * this->k_r[0][i]; };
-        return .5 * dt * this->k_r[this->particles.size() + i];
+    r_funcs[2] = [this](int i, double dt) { return dt * this->k_r[2][i]; };
-    };
-    r_funcs[2] = [this](int i, double dt) {
-        return dt * this->k_r[2 * this->particles.size() + i];
-    };
     r_funcs[3] = [this](int i, double dt) {
-        arma::vec res = this->k_r[i] + this->k_r[this->particles.size() + i] +
+        return (dt / 6) * (this->k_r[0][i].eval() + this->k_r[1][i].eval() +
-                        this->k_r[2 * this->particles.size() + i] +
+                           this->k_r[2][i].eval() + this->k_r[3][i].eval());
-                        this->k_r[3 * this->particles.size() + i];
-        return (dt / 6) * res;
     };
 }
 
 PenningTrap::PenningTrap(int i, double B_0, double V_0, double d)
     : PenningTrap::PenningTrap(B_0, V_0, d)
 {
-    arma::vec r, v;
+    vec_3d r, v;
     for (int j = 0; j < i; j++) {
-        r = arma::vec::fixed<3>().randn() * .1 * this->d;
+        r = vec_3d().randn() * .1 * this->d;
-        v = arma::vec::fixed<3>().randn() * .1 * this->d;
+        v = vec_3d().randn() * .1 * this->d;
         this->add_particle(Particle(1., 40., r, v));
     }
 }
@@ -79,63 +69,60 @@ void PenningTrap::add_particle(Particle particle)
     this->particles.push_back(particle);
 }
 
-arma::vec PenningTrap::external_E_field(arma::vec r)
+vec_3d PenningTrap::external_E_field(vec_3d r)
 {
-    arma::vec::fixed<3> res{r(0), r(1), -2. * r(2)};
+    r(2) *= -2.;
     double f = this->V_0 / (this->d * this->d);
 
-    return f * res;
+    return f * r;
 }
 
-arma::vec PenningTrap::external_B_field(arma::vec r)
+vec_3d PenningTrap::external_B_field(vec_3d r)
 {
-    arma::vec::fixed<3> res{0., 0., this->B_0};
+    return vec_3d{0.,0.,this->B_0};
-
-    return res;
 }
 
-arma::vec PenningTrap::force_on_particle(int i, int j)
+vec_3d PenningTrap::force_on_particle(int i, int j)
 {
     // Calculate the difference between the particles' position
-    arma::vec::fixed<3> res =
+    vec_3d res = this->particles[i].r_vec - this->particles[j].r_vec;
-        this->particles.at(i).r_vec - this->particles.at(j).r_vec;
 
     // Get the distance between the particles
     double norm = arma::norm(res);
 
     // Multiply res with p_j's charge divided by the norm cubed
-    res *= this->particles.at(j).q / (norm * norm * norm);
+    res *= this->particles[j].q / (norm * norm * norm);
 
     return res;
 }
 
-arma::vec PenningTrap::total_force_external(int i)
+vec_3d PenningTrap::total_force_external(int i)
 {
-    if (arma::norm(this->particles.at(i).r_vec) > this->d) {
+    Particle p = this->particles[i];
-        return arma::vec::fixed<3>{0., 0., 0.};
+
+    if (arma::norm(p.r_vec) > this->d) {
+        return vec_3d{0., 0., 0.};
     }
 
-    Particle p = this->particles.at(i);
+    vec_3d B = this->external_B_field(p.r_vec);
 
-    arma::vec::fixed<3> B = this->external_B_field(p.r_vec);
+    vec_3d v_cross_B{p.v_vec(1) * B(2) - p.v_vec(2) * B(1),
+                     p.v_vec(2) * B(0) - p.v_vec(0) * B(2),
+                     p.v_vec(0) * B(1) - p.v_vec(1) * B(0)};
 
-    arma::vec::fixed<3> v_cross_B{p.v_vec(1) * B(2) - p.v_vec(2) * B(1),
+    vec_3d force = p.q * (this->external_E_field(p.r_vec) + v_cross_B);
-                                  p.v_vec(2) * B(0) - p.v_vec(0) * B(2),
-                                  p.v_vec(0) * B(1) - p.v_vec(1) * B(0)};
-
-    arma::vec::fixed<3> force =
-        p.q * (this->external_E_field(p.r_vec) + v_cross_B);
 
     return force;
 }
 
-arma::vec PenningTrap::total_force_particles(int i)
+vec_3d PenningTrap::total_force_particles(unsigned int i)
 {
-    Particle p = this->particles.at(i);
+    Particle p = this->particles[i];
 
-    arma::vec::fixed<3> res;
+    vec_3d res;
+    size_t size = this->particles.size();
 
-    for (int j = 0; j < this->particles.size(); j++) {
+    for (size_t j = 0; j < size; j++) {
         if (i == j) {
             continue;
         }
@@ -148,7 +135,7 @@ arma::vec PenningTrap::total_force_particles(int i)
     return res;
 }
 
-arma::vec PenningTrap::total_force(int i)
+vec_3d PenningTrap::total_force(int i)
 {
     return this->total_force_external(i) - this->total_force_particles(i);
 }
@@ -159,7 +146,7 @@ void PenningTrap::evolve_RK4(double dt, bool particle_interaction)
     std::vector<Particle> original_particles = this->particles;
     std::vector<Particle> tmp_particles = this->particles;
 
-    std::function<arma::vec(int)> force;
+    std::function<vec_3d(int)> force;
     if (particle_interaction) {
         force = [this](int i) { return this->total_force(i); };
     }
@@ -169,25 +156,22 @@ void PenningTrap::evolve_RK4(double dt, bool particle_interaction)
 
     int size = this->particles.size();
 
-    k_v = new arma::vec::fixed<3>[size * 4];
+    this->k_v = sim_arr(4, sim_cols(size));
-    k_r = new arma::vec::fixed<3>[size * 4];
+    this->k_r = sim_arr(4, sim_cols(size));
 
     for (int i = 0; i < 4; i++) {
 #pragma omp parallel for
         for (int j = 0; j < size; j++) {
-            k_v[i * size + j] = this->total_force(j) / this->particles.at(j).m;
+            this->k_v[i][j] = this->total_force(j) / this->particles[j].m;
-            k_r[i * size + j] = this->particles.at(j).v_vec;
+            this->k_r[i][j] = this->particles[j].v_vec;
 
-            Particle *p = &tmp_particles.at(j);
+            Particle *p = &tmp_particles[j];
 
-            p->v_vec = original_particles.at(j).v_vec + this->v_funcs[i](j, dt);
+            p->v_vec = original_particles[j].v_vec + this->v_funcs[i](j, dt);
-            p->r_vec = original_particles.at(j).r_vec + this->r_funcs[i](j, dt);
+            p->r_vec = original_particles[j].r_vec + this->r_funcs[i](j, dt);
         }
         this->particles = tmp_particles;
     }
-
-    delete[] k_v;
-    delete[] k_r;
 }
 
 void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction)
@@ -196,7 +180,7 @@ void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction)
 
     Particle *p;
 
-    std::function<arma::vec(int)> force;
+    std::function<vec_3d(int)> force;
     if (particle_interaction) {
         force = [this](int i) { return this->total_force(i); };
     }
@@ -205,19 +189,20 @@ void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction)
     }
 
 #pragma omp parallel for private(p)
-    for (int i = 0; i < this->particles.size(); i++) {
+    for (size_t i = 0; i < this->particles.size(); i++) {
-        p = &new_state.at(i);
+        p = &new_state[i];
-        p->v_vec += dt * force(i) / new_state.at(i).m;
+        p->v_vec += dt * force(i) / p->m;
-        p->r_vec += dt * this->particles.at(i).v_vec;
+        p->r_vec += dt * this->particles[i].v_vec;
     }
 
     this->particles = new_state;
 }
 
-sim_arr PenningTrap::simulate(double time, int steps, std::string method,
+sim_arr PenningTrap::simulate(double time, unsigned int steps,
-                              bool particle_interaction)
+                              std::string method, bool particle_interaction)
 {
     double dt = time / (double)steps;
+
     sim_arr res(this->particles.size(), sim_cols(steps));
 
     std::function<void(double, bool)> func;
@@ -236,9 +221,8 @@ sim_arr PenningTrap::simulate(double time, int steps, std::string method,
         abort();
     }
 
-    for (int j = 0; j < steps; j++) {
+    for (size_t j = 0; j < steps; j++) {
-        #pragma omp parallel for
+        for (size_t i = 0; i < this->particles.size(); i++) {
-        for (int i = 0; i < this->particles.size(); i++) {
             res[i][j] = this->particles[i].r_vec;
         }
         func(dt, particle_interaction);
@@ -261,17 +245,13 @@ void PenningTrap::write_simulation_to_dir(std::string path, double time,
 
     sim_arr res = this->simulate(time, steps, method, particle_interaction);
 
-    int number_of_particles = this->particles.size();
-
     std::ofstream ofile;
-    sim_rows row;
 
-#pragma omp parallel for private(ofile, row)
+#pragma omp parallel for private(ofile)
-    for (int i = 0; i < number_of_particles; i++) {
+    for (size_t i = 0; i < this->particles.size(); i++) {
-        row = res[i];
         ofile.open(path + "particle_" + std::to_string(i) + ".txt");
-        for (arma::vec vector : row) {
+        for (vec_3d &vec : res[i]) {
-            ofile << vector(0) << "," << vector(1) << "," << vector(2) << "\n";
+            ofile << vec(0) << "," << vec(1) << "," << vec(2) << "\n";
         }
         ofile.close();
     }

src/main.cpp

@@ -24,8 +24,8 @@
 #define CHARGE 1.
 #define MASS 40. // unit: amu
 
-Particle p1(CHARGE, MASS, arma::vec{20., 0., 20.}, arma::vec{0., 25., 0.});
+Particle p1(CHARGE, MASS, vec_3d{20., 0., 20.}, vec_3d{0., 25., 0.});
-Particle p2(CHARGE, MASS, arma::vec{25., 25., 0.}, arma::vec{0., 40., 5.});
+Particle p2(CHARGE, MASS, vec_3d{25., 25., 0.}, vec_3d{0., 40., 5.});
 
 void simulate_single_particle()
 {

src/test_suite.cpp

@@ -44,7 +44,7 @@ public:
         arma::vec result;
         arma::vec v;
         std::stringstream msg;
-        for (int i = 0; i < tests.size(); i++) {
+        for (size_t i = 0; i < tests.size(); i++) {
             v = tests.at(i).first;
             result = trap.external_E_field(v);
 

src/utils.cpp

@@ -65,7 +65,7 @@ bool arma_vector_close_to(arma::vec &a, arma::vec &b, double tol)
         return false;
     }
 
-    for (int i = 0; i < a.n_elem; i++) {
+    for (size_t i = 0; i < a.n_elem; i++) {
         if (std::abs(a(i) - b(i)) >= tol) {
             return false;
         }