include/PenningTrap.hpp

@@ -36,8 +36,9 @@ 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;
+
-    std::function<vec_3d(int, double)>* r_funcs;
+    vec_3d v_func(int, int, double);
+    vec_3d r_func(int, int, double);
 
 public:
     /** @brief Set B_0, V_0 and d.

src/PenningTrap.cpp

@@ -28,32 +28,6 @@ PenningTrap::PenningTrap(double B_0, double V_0, double d)
     this->V_0 = V_0;
     this->d = d;
 
-    // Create lambda functions for each iteration in RK4
-    this->v_funcs = new std::function<vec_3d(int, double)>[] {
-        [this](int i, double dt) {
-            DEBUG("Inside v_funcs");
-            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)>[] {
-        [this](int i, double dt) {
-            DEBUG("Inside v_funcs");
-            return .5 * dt * this->k_r[0][i];
-        },
-        [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)
@@ -74,6 +48,42 @@ PenningTrap::PenningTrap(std::vector<Particle> particles, double B_0,
     this->particles = particles;
 }
 
+vec_3d PenningTrap::v_func(int i, int j, double dt)
+{
+    switch (i) {
+    case 0:
+        return .5 * dt * this->k_v[0][j];
+    case 1:
+        return .5 * dt * this->k_v[1][j];
+    case 2:
+        return dt * this->k_v[2][j];
+    case 3:
+        return (dt / 6.) * (this->k_v[0][j].eval() + this->k_v[1][j].eval() +
+                            this->k_v[2][j].eval() + this->k_v[3][j].eval());
+    default:
+        std::cout << "Not valid!" << std::endl;
+        abort();
+    }
+}
+
+vec_3d PenningTrap::r_func(int i, int j, double dt)
+{
+    switch (i) {
+    case 0:
+        return .5 * dt * this->k_r[0][j];
+    case 1:
+        return .5 * dt * this->k_r[1][j];
+    case 2:
+        return dt * this->k_r[2][j];
+    case 3:
+        return (dt / 6.) * (this->k_r[0][j].eval() + this->k_r[1][j].eval() +
+                            this->k_r[2][j].eval() + this->k_r[3][j].eval());
+    default:
+        std::cout << "Not valid!" << std::endl;
+        abort();
+    }
+}
+
 void PenningTrap::add_particle(Particle particle)
 {
     this->particles.push_back(particle);
@@ -184,8 +194,8 @@ void PenningTrap::evolve_RK4(double dt, bool particle_interaction)
             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->v_vec = original_particles[j].v_vec + this->v_func(i, j, dt);
-            p->r_vec = original_particles[j].r_vec + this->r_funcs[i](j, dt);
+            p->r_vec = original_particles[j].r_vec + this->r_func(i, j, dt);
         }
         DEBUG("After inner loop");
         this->particles = tmp_particles;