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;
         }