Fix v X B

include/PenningTrap.hpp

@@ -75,6 +75,10 @@ public:
     /** @brief Go forward one timestep using the forward Euler method
      * */
     void evolve_forward_euler(double dt);
+
+    arma::vec get_particle(int i);
+
+    double get_d();
 };
 
 #endif

src/PenningTrap.cpp

@@ -9,21 +9,19 @@
  *
  *  @bug No known bugs
  *
- *  @todo Implement add_particle
- *  @todo Implement external_E_field
- *  @todo Implement external_B_field
- *  @todo Implement force_on_particle
- *  @todo Implement total_force_external
- *  @todo Implement total_force_particles
- *  @todo Implement total_force
  *  @todo Implement evolve_RK4
  *  @todo Implement evolve_forward_euler
  * */
 
+#include "utils.hpp"
 #include "PenningTrap.hpp"
 #include "constants.hpp"
 #include <algorithm>
 #include <stdexcept>
+#include <omp.h>
+
+#pragma omp declare reduction( + : arma::vec : omp_out += omp_in ) \
+  initializer( omp_priv = omp_orig )
 
 PenningTrap::PenningTrap(double B_0, double V_0, double d)
 {
@@ -40,13 +38,12 @@ void PenningTrap::add_particle(Particle particle)
 arma::vec PenningTrap::external_E_field(arma::vec r)
 {
     arma::vec::fixed<3> res;
-    double x = r(0), y = r(1), z = r(2);
-    double f = this->V_0/2*this->d*this->d;
-    res(0) = f*2*x;
-    res(1) = f*2*y;
-    res(2) = -f*4*z;
+    double f = this->V_0/(this->d*this->d);
+    res(0) = r(0);
+    res(1) = r(1);
+    res(2) = -2.*r(2);
 
-    return res;
+    return f*res;
 }
 
 arma::vec PenningTrap::external_B_field(arma::vec r)
@@ -82,9 +79,9 @@ arma::vec PenningTrap::total_force_external(int i)
 
     arma::vec::fixed<3> B = this->external_B_field(p.r_vec);
 
-    v_cross_B(0) = p.v_vec(2)*B(3) - p.v_vec(3)*B(2);
-    v_cross_B(1) = p.v_vec(1)*B(3) - p.v_vec(3)*B(1);
-    v_cross_B(2) = p.v_vec(1)*B(2) - p.v_vec(2)*B(1);
+    v_cross_B(0) = p.v_vec(1)*B(2) - p.v_vec(2)*B(1);
+    v_cross_B(1) = p.v_vec(2)*B(0) - p.v_vec(0)*B(2);
+    v_cross_B(2) = p.v_vec(0)*B(1) - p.v_vec(1)*B(0);
 
     arma::vec force = p.q
         *(this->external_E_field(p.r_vec) + v_cross_B);
@@ -96,7 +93,7 @@ arma::vec PenningTrap::total_force_particles(int i)
 {
     Particle p = this->particles.at(i);
 
-    arma::vec::fixed<3> res;
+    arma::vec res(3);
 
     for (int j=0; j < this->particles.size(); j++) {
         if (i == j) {
@@ -106,7 +103,7 @@ arma::vec PenningTrap::total_force_particles(int i)
         res += this->force_on_particle(i, j);
     }
 
-    res *= K_E*p.q/p.m;
+    res *= K_E*(p.q/p.m);
 
     return res;
 }
@@ -123,5 +120,26 @@ void PenningTrap::evolve_RK4(double dt)
 
 void PenningTrap::evolve_forward_euler(double dt)
 {
+    std::vector<Particle> new_state = this->particles;
 
+    Particle *p;
+
+    #pragma omp parallel for private(p)
+    for (int i=0; i < this->particles.size(); i++) {
+        p = &new_state.at(i);
+        p->v_vec += dt*this->total_force(i)/new_state.at(i).m;
+        p->r_vec += dt*this->particles.at(i).v_vec;
+    }
+
+    this->particles = new_state;
+}
+
+arma::vec PenningTrap::get_particle(int i)
+{
+    return this->particles.at(i).r_vec;
+}
+
+double PenningTrap::get_d()
+{
+    return this->d;
 }