diff --git a/src/main.cpp b/src/main.cpp index 14e8fa7..1c55279 100644 --- a/src/main.cpp +++ b/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(); diff --git a/src/test_suite.cpp b/src/test_suite.cpp index 17f9c0e..81ea5c1 100644 --- a/src/test_suite.cpp +++ b/src/test_suite.cpp @@ -26,71 +26,68 @@ public: // Vector containing inputs and expected results std::vector> 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(); diff --git a/src/utils.cpp b/src/utils.cpp index 9694610..b55fa34 100644 --- a/src/utils.cpp +++ b/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& v, int width, int prec) +std::string scientific_format(const std::vector &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; }