src/Particle.cpp

@@ -13,8 +13,8 @@
 #include "Particle.hpp"
 
 Particle::Particle(double q, double m,
-                   arma::vec::fixed<3> r_vec,
-                   arma::vec::fixed<3> v_vec)
+                   vec_3d r_vec,
+                   vec_3d v_vec)
 {
     // Giving the particle its properties
     this->q = q;

src/PenningTrap.cpp

@@ -11,9 +11,6 @@
  * */
 
 #include "PenningTrap.hpp"
-#include "constants.hpp"
-#include "typedefs.hpp"
-#include "utils.hpp"
 
 PenningTrap::PenningTrap(double B_0, std::function<double(double)> V_0,
                          double d, double t)
@@ -99,14 +96,13 @@ vec_3d PenningTrap::external_B_field(vec_3d r)
 
 vec_3d PenningTrap::force_on_particle(unsigned int i, unsigned int j)
 {
-    Particle p_j = this->particles[j];
     // Calculate the difference between the particles' position
-    vec_3d res = this->particles[i].r_vec - p_j.r_vec;
+    vec_3d res = this->particles[i].r_vec - this->particles[j].r_vec;
 
     // Get the distance between the particles
     double norm = arma::norm(res, 2);
 
-    return vec_3d(res * p_j.q / (norm * norm * norm));
+    return vec_3d(res * this->particles[j].q / (norm * norm * norm));
 }
 
 vec_3d PenningTrap::total_force_external(unsigned int i)
@@ -143,6 +139,9 @@ vec_3d PenningTrap::total_force_particles(unsigned int i)
 
 vec_3d PenningTrap::total_force(unsigned int i)
 {
+    if (arma::norm(this->particles[i].r_vec) > this->d) {
+        return vec_3d{0., 0., 0.};
+    }
     return this->total_force_external(i) - this->total_force_particles(i);
 }
 
@@ -217,7 +216,7 @@ simulation_t PenningTrap::simulate(double time, unsigned int steps,
     double dt = time / (double)steps;
 
     unsigned int size = this->particles.size();
-    // sim_arr res(this->particles.size(), sim_cols(steps));
+
     simulation_t res{sim_arr(size, sim_cols(steps)),
                      sim_arr(size, sim_cols(steps))};