Some changes

include/PenningTrap.hpp

@@ -36,8 +36,8 @@ private:
     std::vector<Particle> particles; ///< The particles in the Penning trap
     sim_arr k_v;
     sim_arr k_r;
-    std::function<vec_3d(int, double)> v_funcs[4];
+    std::function<vec_3d(int, double)>* v_funcs;
-    std::function<vec_3d(int, double)> r_funcs[4];
+    std::function<vec_3d(int, double)>* r_funcs;
 
 public:
     /** @brief Set B_0, V_0 and d.

src/PenningTrap.cpp

@@ -29,21 +29,31 @@ PenningTrap::PenningTrap(double B_0, double V_0, double d)
     this->d = d;
 
     // Create lambda functions for each iteration in RK4
-    v_funcs[0] = [this](int i, double dt) { return .5 * dt * this->k_v[0][i]; };
+    this->v_funcs = new std::function<vec_3d(int, double)>[] {
-    v_funcs[1] = [this](int i, double dt) { return .5 * dt * this->k_v[1][i]; };
+        [this](int i, double dt) {
-    v_funcs[2] = [this](int i, double dt) { return dt * this->k_v[2][i]; };
+            DEBUG("Inside v_funcs");
-    v_funcs[3] = [this](int i, double dt) {
+            return .5 * dt * this->k_v[0][i];
+        },
+        [this](int i, double dt) { return .5 * dt * this->k_v[1][i]; },
+        [this](int i, double dt) { return dt * this->k_v[2][i]; },
+        [this](int i, double dt) {
             return (dt / 6) * (this->k_v[0][i].eval() + this->k_v[1][i].eval() +
                                this->k_v[2][i].eval() + this->k_v[3][i].eval());
+        }
     };
-
+    this->r_funcs = new std::function<vec_3d(int, double)>[] {
-    r_funcs[0] = [this](int i, double dt) { return .5 * dt * this->k_r[0][i]; };
+        [this](int i, double dt) {
-    r_funcs[1] = [this](int i, double dt) { return .5 * dt * this->k_r[0][i]; };
+            DEBUG("Inside v_funcs");
-    r_funcs[2] = [this](int i, double dt) { return dt * this->k_r[2][i]; };
+            return .5 * dt * this->k_r[0][i];
-    r_funcs[3] = [this](int i, double dt) {
+        },
+        [this](int i, double dt) { return .5 * dt * this->k_r[1][i]; },
+        [this](int i, double dt) { return dt * this->k_r[2][i]; },
+        [this](int i, double dt) {
             return (dt / 6) * (this->k_r[0][i].eval() + this->k_r[1][i].eval() +
                                this->k_r[2][i].eval() + this->k_r[3][i].eval());
+        }
     };
+
 }
 
 PenningTrap::PenningTrap(int i, double B_0, double V_0, double d)
@@ -143,9 +153,11 @@ vec_3d PenningTrap::total_force(int i)
 void PenningTrap::evolve_RK4(double dt, bool particle_interaction)
 {
 
+    DEBUG("Inside RK4");
     std::vector<Particle> original_particles = this->particles;
     std::vector<Particle> tmp_particles = this->particles;
 
+    DEBUG("Choose force func");
     std::function<vec_3d(int)> force;
     if (particle_interaction) {
         force = [this](int i) { return this->total_force(i); };
@@ -154,22 +166,28 @@ void PenningTrap::evolve_RK4(double dt, bool particle_interaction)
         force = [this](int i) { return this->total_force_external(i); };
     }
 
+    DEBUG("Create k_v and k_r");
     int size = this->particles.size();
 
     this->k_v = sim_arr(4, sim_cols(size));
     this->k_r = sim_arr(4, sim_cols(size));
 
+    DEBUG("Loop");
     for (int i = 0; i < 4; i++) {
+        DEBUG("Inside outer loop");
 #pragma omp parallel for
         for (int j = 0; j < size; j++) {
+            DEBUG("Set k_v and k_r");
             this->k_v[i][j] = this->total_force(j) / this->particles[j].m;
             this->k_r[i][j] = this->particles[j].v_vec;
 
             Particle *p = &tmp_particles[j];
 
+            DEBUG("Update v and r");
             p->v_vec = original_particles[j].v_vec + this->v_funcs[i](j, dt);
             p->r_vec = original_particles[j].r_vec + this->r_funcs[i](j, dt);
         }
+        DEBUG("After inner loop");
         this->particles = tmp_particles;
     }
 }
@@ -201,13 +219,16 @@ void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction)
 sim_arr PenningTrap::simulate(double time, unsigned int steps,
                               std::string method, bool particle_interaction)
 {
+    DEBUG("Inside simulate");
     double dt = time / (double)steps;
 
     sim_arr res(this->particles.size(), sim_cols(steps));
 
+    DEBUG("Choose function to use");
     std::function<void(double, bool)> func;
     if (method == "rk4") {
         func = [this](double dt, bool particle_interaction) {
+            DEBUG("Inside RK4 lambda");
             this->evolve_RK4(dt, particle_interaction);
         };
     }
@@ -221,10 +242,13 @@ sim_arr PenningTrap::simulate(double time, unsigned int steps,
         abort();
     }
 
+    DEBUG("Loop");
     for (size_t j = 0; j < steps; j++) {
+        DEBUG("Inside outer loop");
         for (size_t i = 0; i < this->particles.size(); i++) {
             res[i][j] = this->particles[i].r_vec;
         }
+        DEBUG("After inner loop");
         func(dt, particle_interaction);
     }
 
@@ -235,6 +259,7 @@ void PenningTrap::write_simulation_to_dir(std::string path, double time,
                                           int steps, std::string method,
                                           bool particle_interaction)
 {
+    DEBUG("Inside write to dir");
     if (path.back() != '/') {
         path += '/';
     }
@@ -243,10 +268,12 @@ void PenningTrap::write_simulation_to_dir(std::string path, double time,
         return;
     }
 
+    DEBUG("Create sim_arr");
     sim_arr res = this->simulate(time, steps, method, particle_interaction);
 
     std::ofstream ofile;
 
+    DEBUG("Loop");
 #pragma omp parallel for private(ofile)
     for (size_t i = 0; i < this->particles.size(); i++) {
         ofile.open(path + "particle_" + std::to_string(i) + ".txt");

src/main.cpp

@@ -29,10 +29,12 @@ Particle p2(CHARGE, MASS, vec_3d{25., 25., 0.}, vec_3d{0., 40., 5.});
 
 void simulate_single_particle()
 {
+    DEBUG("Inside single particle sim");
     PenningTrap trap(std::vector<Particle>{p1});
 
     double time = 50.; // microseconds
 
+    DEBUG("Write to dir");
     trap.write_simulation_to_dir("output/simulate_single_particle", time, N);
 }
 
@@ -82,14 +84,19 @@ void simulate_100_particles()
 
 int main()
 {
+    DEBUG("Start");
     double start = omp_get_wtime();
 
+    DEBUG("Single particle");
     simulate_single_particle();
 
+    DEBUG("Two particles");
     simulate_two_particles();
 
+    DEBUG("Single particle with different steps");
     simulate_single_particle_with_different_steps();
 
+    DEBUG("100 particles");
     simulate_100_particles();
 
     double end = omp_get_wtime();