.clang-format

@@ -1,4 +1,8 @@
+UseTab: Never
 IndentWidth: 4
+TabWidth: 4
+AccessModifierOffset: -4
+IndentAccessModifiers: false
 AllowShortFunctionsOnASingleLine: false
 BreakBeforeBraces: Custom
 BraceWrapping:

src/PenningTrap.cpp

@@ -103,6 +103,65 @@ arma::vec PenningTrap::total_force(int i)
 
 void PenningTrap::evolve_RK4(double dt)
 {
+    std::vector<Particle> 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];
+
+    int size = this->particles.size();
+
+    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; i<size; i++) {
+        Particle *p = &this->particles.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];
+    }
+
+
+    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; i<size; i++) {
+        Particle *p = &this->particles.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];
+    }
+
+    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; i<size; i++) {
+        Particle *p = &this->particles.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];
+    }
+
+
+    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; i<size; i++) {
+        Particle *p = &this->particles.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;
+    }
+
+    delete [] k_v;
+    delete [] k_r;
 }
 
 void PenningTrap::evolve_forward_euler(double dt)

src/animate_100_particles.py

@@ -7,11 +7,11 @@ def get_data(files):
     res = []
     for file in files:
         arr = [[], [], []]
-        with open(file) as f:
+        with open(file, encoding="utf8") as f:
             lines = f.readlines()
 
             for line in lines:
-                xi,yi,zi = map(lambda x: float(x), line.strip().split(","))
+                xi,yi,zi = map(float, line.strip().split(","))
                 arr[0].append(xi)
                 arr[1].append(yi)
                 arr[2].append(zi)
@@ -27,11 +27,12 @@ def update(num, lines, arr):
 
 
 def animate():
+    plt.style.use("dark_background")
     fig = plt.figure()
     ax = fig.add_subplot(projection="3d")
 
 
-    arr = get_data([f"output/p{i}.txt" for i in range(100)])
+    arr = get_data([f"output/p{i}_RK4.txt" for i in range(100)])
 
     arr = arr[:,:,::10]
 
@@ -53,7 +54,7 @@ def animate():
                                   blit=False)
 
 
-    # ani.save("100_particles.gif", writer=animation.FFMpegFileWriter(fps=30))
+    ani.save("100_particles.gif", writer=animation.FFMpegFileWriter(fps=30))
     plt.show()
 
 

src/main.cpp

@@ -21,52 +21,54 @@
 #define CHARGE 1.
 #define MASS 40. // unit: amu
 
-void euler_100_particles()
+void simulate_100_particles()
 {
     PenningTrap trap;
 
     // Add particles inside trap
-    for (int i=0; i < PARTICLES; i++) {
-        arma::vec r = arma::vec(3).randn() * 0.1 * trap.get_d();  // random initial position
-        arma::vec v = arma::vec(3).randn() * 0.1 * trap.get_d();  // random initial velocity
+    for (int i = 0; i < PARTICLES; i++) {
+        arma::vec r = arma::vec(3).randn() * 0.1 *
+                      trap.get_d(); // random initial position
+        arma::vec v = arma::vec(3).randn() * 0.1 *
+                      trap.get_d(); // random initial velocity
         trap.add_particle(Particle(CHARGE, MASS, r, v));
     }
 
     double time = 50.; // microseconds
-    double dt = time / (double) N;
+    double dt = time / (double)N;
 
     auto res = new arma::vec::fixed<3>[PARTICLES][N];
 
     int counter = 0;
 
     // Get the path of all particles
-    for (int j=0; j < N; j++) {
-        #pragma omp parallel for 
-        for (int i=0; i < PARTICLES; i++) {
+    for (int j = 0; j < N; j++) {
+#pragma omp parallel for
+        for (int i = 0; i < PARTICLES; i++) {
             res[i][j] = trap.get_particle(i);
         }
-        trap.evolve_forward_euler(dt);
+        trap.evolve_RK4(dt);
     }
 
     std::cout << counter << std::endl;
 
-    arma::vec::fixed<3>* cur_row;
+    arma::vec::fixed<3> *cur_row;
     arma::vec::fixed<3> cur_elem;
 
     mkdir("output", 0777);
+    mkdir("output/simulate_100_particles", 0777);
 
     std::ofstream ofile;
 
-    // Write particle paths to file
-    #pragma omp parallel for private(cur_row, cur_elem, ofile)
-    for (int i=0; i < PARTICLES; i++) {
+// Write particle paths to file
+#pragma omp parallel for private(cur_row, cur_elem, ofile)
+    for (int i = 0; i < PARTICLES; i++) {
         cur_row = res[i];
-        ofile.open("output/p" + std::to_string(i) + ".txt");
-        for (int j=0; j < N; j++) {
+        ofile.open("output/simulate_100_particles/p" + std::to_string(i) + ".txt");
+        for (int j = 0; j < N; j++) {
             cur_elem = cur_row[j];
-            ofile << cur_elem(0) << "," 
-                << cur_elem(1) << "," 
-                << cur_elem(2) << "\n";
+            ofile << cur_elem(0) << "," << cur_elem(1) << "," << cur_elem(2)
+                  << "\n";
         }
         ofile.close();
     }
@@ -76,7 +78,7 @@ int main()
 {
     double start = omp_get_wtime();
 
-    euler_100_particles();
+    simulate_100_particles();
 
     double end = omp_get_wtime();
 

src/test_suite.cpp

@@ -26,71 +26,68 @@ public:
         // Vector containing inputs and expected results
         std::vector<std::pair<arma::vec, arma::vec>> tests;
 
-        tests.push_back(std::make_pair(arma::vec{0.,0.,0.}, 
-                                       arma::vec{0.,0.,0.}));
+        tests.push_back(
+            std::make_pair(arma::vec{0., 0., 0.}, arma::vec{0., 0., 0.}));
 
-        tests.push_back(std::make_pair(arma::vec{10.,0.,0.}, 
-                                       arma::vec{96.4852558,0.,0.}));
+        tests.push_back(std::make_pair(arma::vec{10., 0., 0.},
+                                       arma::vec{96.4852558, 0., 0.}));
 
-        tests.push_back(std::make_pair(arma::vec{10.,0.,0.}, 
-                                       arma::vec{96.4852558,0.,0.}));
+        tests.push_back(std::make_pair(arma::vec{10., 0., 0.},
+                                       arma::vec{96.4852558, 0., 0.}));
 
-        tests.push_back(std::make_pair(arma::vec{0.,10.,0.}, 
-                                       arma::vec{0.,96.4852558,0.}));
-
-        tests.push_back(std::make_pair(arma::vec{0.,0.,10.}, 
-                                       arma::vec{0.,0.,-192.9705116}));
+        tests.push_back(std::make_pair(arma::vec{0., 10., 0.},
+                                       arma::vec{0., 96.4852558, 0.}));
 
+        tests.push_back(std::make_pair(arma::vec{0., 0., 10.},
+                                       arma::vec{0., 0., -192.9705116}));
 
         arma::vec result;
         arma::vec v;
         std::stringstream msg;
-        for (int i=0; i < tests.size(); i++) {
+        for (int i = 0; i < tests.size(); i++) {
             v = tests.at(i).first;
             result = trap.external_E_field(v);
 
             msg.str("");
-            msg << "Testing the external E field at (" << std::setprecision(2) 
+            msg << "Testing the external E field at (" << std::setprecision(2)
                 << v(0) << "," << v(1) << "," << v(2) << ").";
 
             ASSERT(arma_vector_close_to(result, tests.at(i).second), msg.str());
         }
     }
 
-
     static void test_external_B_field()
     {
-        // No point in testing at different points since it's not dependent on 
-        // position.
+        // No point in testing at different points since it's not dependent
+        // on position.
         PenningTrap trap;
-        arma::vec expected{0.,0.,T};
-        arma::vec result = trap.external_B_field(arma::vec{0.,0.,0.});
+        arma::vec expected{0., 0., T};
+        arma::vec result = trap.external_B_field(arma::vec{0., 0., 0.});
         ASSERT(arma_vector_close_to(expected, result),
                "Testing the external B field at (0,0,0)");
     }
 
-
     static void test_force_on_particle()
     {
         PenningTrap trap;
-        arma::vec v{0.,0.,0.};
+        arma::vec v{0., 0., 0.};
 
         // Add particles to test
-        trap.add_particle(Particle(1.,40.,arma::vec{0.,0.,0.},v));
-        trap.add_particle(Particle(1.,40.,arma::vec{1.,0.,0.},v));
-        trap.add_particle(Particle(1.,40.,arma::vec{0.,3.,4.},v));
+        trap.add_particle(Particle(1., 40., arma::vec{0., 0., 0.}, v));
+        trap.add_particle(Particle(1., 40., arma::vec{1., 0., 0.}, v));
+        trap.add_particle(Particle(1., 40., arma::vec{0., 3., 4.}, v));
 
         // Test p0 and p1
-        arma::vec expected{-1.,0.,0.};
+        arma::vec expected{-1., 0., 0.};
         arma::vec result = trap.force_on_particle(0, 1);
         ASSERT(arma_vector_close_to(expected, result),
-               "Testing the force on a particle at (0,0,0) from a " 
-               "particle at (1,0,0)."); 
+               "Testing the force on a particle at (0,0,0) from a "
+               "particle at (1,0,0).");
 
         // Test p0 and p2
         expected = arma::vec{0, -.024, -.032};
         result = trap.force_on_particle(0, 2);
-        ASSERT(arma_vector_close_to(expected, result), 
+        ASSERT(arma_vector_close_to(expected, result),
                "Testing the force on a particle at (0,0,0) from a "
                "particle at (0,3,4).");
     }
@@ -98,10 +95,10 @@ public:
     static void test_total_force_external()
     {
         PenningTrap trap;
-        trap.add_particle(Particle(1.,40.,arma::vec{1.,2.,3.},
-                                   arma::vec{3.,4.,5.}));
+        trap.add_particle(
+            Particle(1., 40., arma::vec{1., 2., 3.}, arma::vec{3., 4., 5.}));
 
-        arma::vec expected{395.58954878,-270.15871624,-57.89115348};
+        arma::vec expected{395.58954878, -270.15871624, -57.89115348};
         arma::vec result = trap.total_force_external(0);
         ASSERT(arma_vector_close_to(expected, result),
                "Testing the total external force on a particle at "
@@ -111,23 +108,23 @@ public:
     static void test_total_force_particles()
     {
         PenningTrap trap;
-        trap.add_particle(Particle(1.,40.,arma::vec{0.,0.,0.},
-                                   arma::vec{0.,0.,0.}));
+        trap.add_particle(
+            Particle(1., 40., arma::vec{0., 0., 0.}, arma::vec{0., 0., 0.}));
 
-        arma::vec expected{0.,0.,0.};
+        arma::vec expected{0., 0., 0.};
         arma::vec result = trap.total_force_particles(0);
         ASSERT(arma_vector_close_to(expected, result),
                "Testing the total force of all particles on particle 0 "
                "with only a single particle");
 
-        trap.add_particle(Particle(1.,40.,arma::vec{1.,0.,0.},
-                                   arma::vec{0.,0.,0.}));
-        trap.add_particle(Particle(1.,40.,arma::vec{0.,1.,0.},
-                                   arma::vec{0.,0.,0.}));
-        trap.add_particle(Particle(1.,40.,arma::vec{0.,0.,1.},
-                                   arma::vec{0.,0.,0.}));
+        trap.add_particle(
+            Particle(1., 40., arma::vec{1., 0., 0.}, arma::vec{0., 0., 0.}));
+        trap.add_particle(
+            Particle(1., 40., arma::vec{0., 1., 0.}, arma::vec{0., 0., 0.}));
+        trap.add_particle(
+            Particle(1., 40., arma::vec{0., 0., 1.}, arma::vec{0., 0., 0.}));
 
-        expected = arma::vec(3,arma::fill::value(-3473.383325));
+        expected = arma::vec(3, arma::fill::value(-3473.383325));
         result = trap.total_force_particles(0);
         ASSERT(arma_vector_close_to(expected, result),
                "Testing the total force of all particles on particle 0 "
@@ -135,8 +132,7 @@ public:
     }
 };
 
-
-int main() 
+int main()
 {
     PenningTrapTest::test_external_E_field();
     PenningTrapTest::test_external_B_field();

src/utils.cpp

@@ -18,10 +18,10 @@ std::string scientific_format(double d, int width, int prec)
     return ss.str();
 }
 
-std::string scientific_format(const std::vector<double>& v, int width, int prec)
+std::string scientific_format(const std::vector<double> &v, int width, int prec)
 {
     std::stringstream ss;
-    for(double elem : v) {
+    for (double elem : v) {
         ss << scientific_format(elem, width, prec);
     }
     return ss.str();
@@ -37,12 +37,8 @@ static void print_message(std::string msg)
     }
 }
 
-void m_assert(bool expr, 
-              std::string expr_str, 
-              std::string f, 
-              std::string file, 
-              int line,
-              std::string msg)
+void m_assert(bool expr, std::string expr_str, std::string f, std::string file,
+              int line, std::string msg)
 {
     std::string new_assert(f.size() + (expr ? 4 : 6), '-');
     std::cout << "\x1B[36m" << new_assert << "\033[0m\n";
@@ -54,8 +50,8 @@ void m_assert(bool expr,
     else {
         std::cout << "\x1B[31mFAIL\033[0m\n";
         print_message(msg);
-        std::cout << file <<  " " <<  line << ": Assertion \"" 
-            << expr_str << "\" Failed\n\n";
+        std::cout << file << " " << line << ": Assertion \"" << expr_str
+                  << "\" Failed\n\n";
         abort();
     }
 }
@@ -66,7 +62,7 @@ bool arma_vector_close_to(arma::vec &a, arma::vec &b, double tol)
         return false;
     }
 
-    for (int i=0; i < a.n_elem; i++) {
+    for (int i = 0; i < a.n_elem; i++) {
         if (std::abs(a(i) - b(i)) >= tol) {
             return false;
         }