Coryab/implement penning trap simulate #6

Merged
coryab merged 16 commits from coryab/implement-PenningTrap-simulate into develop 2023-10-14 01:13:37 +00:00
2 changed files with 36 additions and 10 deletions
Showing only changes of commit 0b52008d4d - Show all commits

View File

@ -42,6 +42,12 @@ public:
* */
PenningTrap(double B_0 = T, double V_0 = 25. * V / 1000., double d = 500.);
PenningTrap(int i, double B_0 = T, double V_0 = 25. * V / 1000.,
double d = 500.);
PenningTrap(std::vector<Particle> particles, double B_0 = T,
double V_0 = 25. * V / 1000., double d = 500.);
/** @brief Add a particle to the system
* */
void add_particle(Particle particle);

View File

@ -17,6 +17,7 @@
#include "constants.hpp"
#include "utils.hpp"
#include <cstdlib>
#include <vector>
PenningTrap::PenningTrap(double B_0, double V_0, double d)
{
@ -25,6 +26,29 @@ PenningTrap::PenningTrap(double B_0, double V_0, double d)
this->d = d;
}
PenningTrap::PenningTrap(int i, double B_0, double V_0, double d)
{
this->B_0 = B_0;
this->V_0 = V_0;
this->d = d;
arma::vec r, v;
for (int j = 0; j < i; j++) {
r = arma::vec(3).randn() * .1 * this->d;
v = arma::vec(3).randn() * .1 * this->d;
this->add_particle(Particle(1., 40., r, v));
}
}
PenningTrap::PenningTrap(std::vector<Particle> particles, double B_0,
double V_0, double d)
{
this->B_0 = B_0;
this->V_0 = V_0;
this->d = d;
this->particles = particles;
}
void PenningTrap::add_particle(Particle particle)
{
this->particles.push_back(particle);
@ -196,14 +220,10 @@ sim_arr PenningTrap::simulate(double time, int steps, std::string method)
std::function<void(double)> func;
if (method == "rk4") {
func = [this](double dt) {
this->evolve_RK4(dt);
};
func = [this](double dt) { this->evolve_RK4(dt); };
}
else if(method == "euler") {
func = [this](double dt) {
this->evolve_forward_euler(dt);
};
else if (method == "euler") {
func = [this](double dt) { this->evolve_forward_euler(dt); };
}
else {
std::cout << "Not a valid method!" << std::endl;
@ -212,8 +232,8 @@ sim_arr PenningTrap::simulate(double time, int steps, std::string method)
int size = this->particles.size();
for (int j=0; j<steps; j++) {
for (int i=0;i<size;i++) {
for (int j = 0; j < steps; j++) {
for (int i = 0; i < size; i++) {
res[i][j] = this->particles[i].r_vec;
}
func(dt);