diff --git a/include/PenningTrap.hpp b/include/PenningTrap.hpp index 74759e2..eb7bff7 100644 --- a/include/PenningTrap.hpp +++ b/include/PenningTrap.hpp @@ -36,8 +36,9 @@ private: std::vector particles; ///< The particles in the Penning trap sim_arr k_v; sim_arr k_r; - std::function* v_funcs; - std::function* 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. diff --git a/src/PenningTrap.cpp b/src/PenningTrap.cpp index 4ab38e9..0c50836 100644 --- a/src/PenningTrap.cpp +++ b/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[] { - [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[] { - [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 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->r_vec = original_particles[j].r_vec + this->r_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_func(i, j, dt); } DEBUG("After inner loop"); this->particles = tmp_particles;