diff --git a/include/PenningTrap.hpp b/include/PenningTrap.hpp index 3601634..1c7d56c 100644 --- a/include/PenningTrap.hpp +++ b/include/PenningTrap.hpp @@ -18,8 +18,12 @@ #include "Particle.hpp" #include "constants.hpp" -#pragma omp declare reduction( + : arma::vec : omp_out += omp_in ) \ - initializer( omp_priv = omp_orig ) +typedef std::vector> sim_cols; +typedef std::vector> sim_rows; +typedef std::vector sim_arr; + +#pragma omp declare reduction(+ : arma::vec : omp_out += omp_in) \ + initializer(omp_priv = omp_orig) /** @brief A class that simulates a Penning trap. * @@ -28,15 +32,15 @@ * */ class PenningTrap { private: - double B_0; ///< Magnetic field strength - double V_0; ///< Applied potential - double d; ///< Characteristic dimension + double B_0; ///< Magnetic field strength + double V_0; ///< Applied potential + double d; ///< Characteristic dimension std::vector particles; ///< The particles in the Penning trap public: /** @brief Set B_0, V_0 and d. * */ - PenningTrap(double B_0 = T, double V_0 = 25.*V/1000., double d = 500.); + PenningTrap(double B_0 = T, double V_0 = 25. * V / 1000., double d = 500.); /** @brief Add a particle to the system * */ @@ -52,7 +56,7 @@ public: /** @brief Calculate the force between 2 particles. * - * @details Calculate the force exhibited on particle p_i from + * @details Calculate the force exhibited on particle p_i from * particle p_j. * */ arma::vec force_on_particle(int i, int j); @@ -69,7 +73,7 @@ public: arma::vec total_force_particles(int i); /** @brief calculate the total force on a particle. - * */ + * */ arma::vec total_force(int i); /** @brief Go forward one timestep using the RK4 method @@ -80,6 +84,8 @@ public: * */ void evolve_forward_euler(double dt); + sim_arr simulate(double time, int steps, std::string method = "rk4"); + arma::vec get_particle(int i); double get_d(); diff --git a/src/PenningTrap.cpp b/src/PenningTrap.cpp index 5d9a16b..c4fea93 100644 --- a/src/PenningTrap.cpp +++ b/src/PenningTrap.cpp @@ -16,6 +16,7 @@ #include "PenningTrap.hpp" #include "constants.hpp" #include "utils.hpp" +#include PenningTrap::PenningTrap(double B_0, double V_0, double d) { @@ -104,64 +105,72 @@ arma::vec PenningTrap::total_force(int i) void PenningTrap::evolve_RK4(double dt) { std::vector tmp_particles = this->particles; - - arma::vec::fixed<3> *k_v = new arma::vec::fixed<3>[this->particles.size()*4]; - arma::vec::fixed<3> *k_r = new arma::vec::fixed<3>[this->particles.size()*4]; + + arma::vec::fixed<3> *k_v = + new arma::vec::fixed<3>[this->particles.size() * 4]; + arma::vec::fixed<3> *k_r = + new arma::vec::fixed<3>[this->particles.size() * 4]; int size = this->particles.size(); - for (int i=0; itotal_force(i)/this->particles.at(i).m; + for (int i = 0; i < size; i++) { + k_v[i] = this->total_force(i) / this->particles.at(i).m; k_r[i] = this->particles.at(i).v_vec; } - for (int i=0; iparticles.at(i); - p->v_vec = tmp_particles.at(i).v_vec + (dt/2)*k_v[i]; - p->r_vec = tmp_particles.at(i).r_vec + (dt/2)*k_r[i]; + p->v_vec = tmp_particles.at(i).v_vec + (dt / 2) * k_v[i]; + p->r_vec = tmp_particles.at(i).r_vec + (dt / 2) * k_r[i]; } - - for (int i=0; itotal_force(i)/this->particles.at(i).m; - k_r[1*size + i] = this->particles.at(i).v_vec; + for (int i = 0; i < size; i++) { + k_v[1 * size + i] = this->total_force(i) / this->particles.at(i).m; + k_r[1 * size + i] = this->particles.at(i).v_vec; } - for (int i=0; iparticles.at(i); - p->v_vec = tmp_particles.at(i).v_vec + (dt/2)*k_v[1*size + i]; - p->r_vec = tmp_particles.at(i).r_vec + (dt/2)*k_r[1*size + i]; + p->v_vec = tmp_particles.at(i).v_vec + (dt / 2) * k_v[1 * size + i]; + p->r_vec = tmp_particles.at(i).r_vec + (dt / 2) * k_r[1 * size + i]; } - for (int i=0; itotal_force(i)/this->particles.at(i).m; - k_r[2*size + i] = this->particles.at(i).v_vec; + for (int i = 0; i < size; i++) { + k_v[2 * size + i] = this->total_force(i) / this->particles.at(i).m; + k_r[2 * size + i] = this->particles.at(i).v_vec; } - for (int i=0; iparticles.at(i); - p->v_vec = tmp_particles.at(i).v_vec + dt*k_v[2*size + i]; - p->r_vec = tmp_particles.at(i).r_vec + dt*k_r[2*size + i]; + p->v_vec = tmp_particles.at(i).v_vec + dt * k_v[2 * size + i]; + p->r_vec = tmp_particles.at(i).r_vec + dt * k_r[2 * size + i]; } - - for (int i=0; itotal_force(i)/this->particles.at(i).m; - k_r[3*size + i] = this->particles.at(i).v_vec; + for (int i = 0; i < size; i++) { + k_v[3 * size + i] = this->total_force(i) / this->particles.at(i).m; + k_r[3 * size + i] = this->particles.at(i).v_vec; } - for (int i=0; iparticles.at(i); - p->v_vec = tmp_particles.at(i).v_vec + dt*(k_v[i] + k_v[size + i] + k_v[2*size + i] + k_v[3*size + i])/6; - p->r_vec = tmp_particles.at(i).r_vec + dt*(k_r[i] + k_r[size + i] + k_r[2*size + i] + k_r[3*size + i])/6; + p->v_vec = tmp_particles.at(i).v_vec + + dt * + (k_v[i] + k_v[size + i] + k_v[2 * size + i] + + k_v[3 * size + i]) / + 6; + p->r_vec = tmp_particles.at(i).r_vec + + dt * + (k_r[i] + k_r[size + i] + k_r[2 * size + i] + + k_r[3 * size + i]) / + 6; } - delete [] k_v; - delete [] k_r; + delete[] k_v; + delete[] k_r; } void PenningTrap::evolve_forward_euler(double dt) @@ -180,6 +189,39 @@ void PenningTrap::evolve_forward_euler(double dt) this->particles = new_state; } +sim_arr PenningTrap::simulate(double time, int steps, std::string method) +{ + double dt = time / (double)steps; + sim_arr res(this->particles.size(), sim_cols(steps)); + + std::function func; + if (method == "rk4") { + func = [this](double dt) { + this->evolve_RK4(dt); + }; + } + else if(method == "euler") { + func = [this](double dt) { + this->evolve_forward_euler(dt); + }; + } + else { + std::cout << "Not a valid method!" << std::endl; + abort(); + } + + int size = this->particles.size(); + + for (int j=0; jparticles[i].r_vec; + } + func(dt); + } + + return res; +} + arma::vec PenningTrap::get_particle(int i) { return this->particles.at(i).r_vec;