From 0dfb102cef2372753c8c57779928b6b3c388cb9f Mon Sep 17 00:00:00 2001 From: Cory Date: Mon, 23 Oct 2023 12:41:03 +0200 Subject: [PATCH] Update versions and some other things - Update version number for Doxygen - Use vec3 instead of vec_3d --- src/Particle.cpp | 4 +- src/PenningTrap.cpp | 91 ++++++++++++++-------------- src/main.cpp | 19 +++--- src/scripts/animate_100_particles.py | 2 +- src/scripts/plot_particles_left.py | 82 ++++++++++++++++++------- src/test_suite.cpp | 87 +++++++++++++------------- src/utils.cpp | 3 +- 7 files changed, 160 insertions(+), 128 deletions(-) diff --git a/src/Particle.cpp b/src/Particle.cpp index e1f837b..9d55274 100644 --- a/src/Particle.cpp +++ b/src/Particle.cpp @@ -3,7 +3,7 @@ * @author Cory Alexander Balaton (coryab) * @author Janita Ovidie Sandtrøen Willumsen (janitaws) * - * @version 0.1 + * @version 1.0 * * @brief The implementation of the Particle class. * @@ -12,7 +12,7 @@ #include "Particle.hpp" -Particle::Particle(vec_3d r_vec, vec_3d v_vec, double q, double m) +Particle::Particle(vec3 r_vec, vec3 v_vec, double q, double m) { // Giving the particle its properties this->r_vec = r_vec; diff --git a/src/PenningTrap.cpp b/src/PenningTrap.cpp index dbf7991..f045bae 100644 --- a/src/PenningTrap.cpp +++ b/src/PenningTrap.cpp @@ -3,7 +3,7 @@ * @author Cory Alexander Balaton (coryab) * @author Janita Ovidie Sandtrøen Willumsen (janitaws) * - * @version 0.1 + * @version 1.0 * * @brief The implementation of the PenningTrap class. * @@ -31,8 +31,8 @@ PenningTrap::PenningTrap(uint i, double B_0, double V_0, double d, double t) { for (size_t j = 0; j < i; j++) { this->particles.push_back( - Particle(vec_3d(vec_3d().randn() * .1 * this->d), - vec_3d(vec_3d().randn() * .1 * this->d))); + Particle(vec3(vec3().randn() * .1 * this->d), + vec3(vec3().randn() * .1 * this->d))); } } @@ -46,7 +46,7 @@ PenningTrap::PenningTrap(std::vector particles, double B_0, void PenningTrap::set_pertubation(double f, double omega_V) { this->perturbation = [f, omega_V](double t) { - return f * std::cos(omega_V * t); + return 1 + f * std::cos(omega_V * t); }; } @@ -56,11 +56,11 @@ void PenningTrap::reinitialize(double f, double omega_V, double t) this->set_pertubation(f, omega_V); for (size_t i = 0; i < this->particles.size(); i++) { - this->particles[i].r_vec = vec_3d().randn() * .1 * this->d; + this->particles[i].r_vec = vec3().randn() * .1 * this->d; } } -vec_3d PenningTrap::v_func(uint i, uint j, double dt) +vec3 PenningTrap::v_func(uint i, uint j, double dt) { switch (i) { case 0: @@ -70,16 +70,16 @@ vec_3d PenningTrap::v_func(uint i, uint j, double dt) case 2: return dt * this->k_v[2][j]; case 3: - return vec_3d((dt / 6.) - * (this->k_v[0][j] + 2. * this->k_v[1][j] - + 2. * this->k_v[2][j] + this->k_v[3][j])); + return vec3((dt / 6.) + * (this->k_v[0][j] + 2. * this->k_v[1][j] + + 2. * this->k_v[2][j] + this->k_v[3][j])); default: std::cout << "Not valid!" << std::endl; abort(); } } -vec_3d PenningTrap::r_func(uint i, uint j, double dt) +vec3 PenningTrap::r_func(uint i, uint j, double dt) { switch (i) { case 0: @@ -89,9 +89,9 @@ vec_3d PenningTrap::r_func(uint i, uint j, double dt) case 2: return dt * this->k_r[2][j]; case 3: - return vec_3d((dt / 6.) - * (this->k_r[0][j] + 2. * this->k_r[1][j] - + 2. * this->k_r[2][j] + this->k_r[3][j])); + return vec3((dt / 6.) + * (this->k_r[0][j] + 2. * this->k_r[1][j] + + 2. * this->k_r[2][j] + this->k_r[3][j])); default: std::cout << "Not valid!" << std::endl; abort(); @@ -103,63 +103,61 @@ void PenningTrap::add_particle(Particle particle) this->particles.push_back(particle); } -vec_3d PenningTrap::external_E_field(vec_3d r) +vec3 PenningTrap::external_E_field(vec3 r) { r(2) *= -2.; - return vec_3d( - (this->V_0 * this->perturbation(this->t) / (this->d * this->d)) * r); + return vec3((this->V_0 * this->perturbation(this->t) / (this->d * this->d)) + * r); } -vec_3d PenningTrap::external_B_field(vec_3d r) +vec3 PenningTrap::external_B_field(vec3 r) { - return vec_3d{0., 0., this->B_0}; + return vec3{0., 0., this->B_0}; } -vec_3d PenningTrap::force_on_particle(uint i, uint j) +vec3 PenningTrap::force_on_particle(uint i, uint j) { // Calculate the difference between the particles' position - vec_3d res = this->particles[i].r_vec - this->particles[j].r_vec; + vec3 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((this->particles[j].q / (norm * norm * norm)) * res); + return vec3((this->particles[j].q / (norm * norm * norm)) * res); } -vec_3d PenningTrap::total_force_external(uint i) +vec3 PenningTrap::total_force_external(uint i) { Particle *p = &this->particles[i]; if (arma::norm(p->r_vec) > this->d) { - return vec_3d{0., 0., 0.}; + return vec3{0., 0., 0.}; } - return vec_3d( - p->q - * (this->external_E_field(p->r_vec) - + arma::cross(p->v_vec, this->external_B_field(p->r_vec)))); + return vec3(p->q + * (this->external_E_field(p->r_vec) + + arma::cross(p->v_vec, this->external_B_field(p->r_vec)))); } -vec_3d PenningTrap::total_force_particles(uint i) +vec3 PenningTrap::total_force_particles(uint i) { - vec_3d res; + vec3 res; for (size_t j = 0; j < this->particles.size(); j++) { if (i != j) res += this->force_on_particle(i, j); } - return vec_3d(res * (K_E * this->particles[i].q)); + return vec3(res * (K_E * this->particles[i].q)); } -vec_3d PenningTrap::total_force(uint i) +vec3 PenningTrap::total_force(uint i) { if (arma::norm(this->particles[i].r_vec) > this->d) { - return vec_3d{0., 0., 0.}; + return vec3{0., 0., 0.}; } - return vec_3d(this->total_force_external(i) - - this->total_force_particles(i)); + return vec3(this->total_force_external(i) - this->total_force_particles(i)); } void PenningTrap::evolve_RK4(double dt, bool particle_interaction) @@ -168,9 +166,9 @@ void PenningTrap::evolve_RK4(double dt, bool particle_interaction) std::vector original_particles = this->particles; std::vector tmp_particles = this->particles; - vec_3d (PenningTrap::*force)(uint) = - particle_interaction ? &PenningTrap::total_force - : &PenningTrap::total_force_external; + vec3 (PenningTrap::*force)(uint) = particle_interaction + ? &PenningTrap::total_force + : &PenningTrap::total_force_external; size_t size = this->particles.size(); @@ -203,12 +201,11 @@ void PenningTrap::evolve_RK4(double dt, bool particle_interaction) void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction) { size_t size = this->particles.size(); - vec_3d force_res[size]; - Particle *p; + vec3 force_res[size]; - vec_3d (PenningTrap::*force)(uint) = - particle_interaction ? &PenningTrap::total_force - : &PenningTrap::total_force_external; + vec3 (PenningTrap::*force)(uint) = particle_interaction + ? &PenningTrap::total_force + : &PenningTrap::total_force_external; // Calculating the force for each particle is independent and therefore // a good candidate for parallel execution @@ -219,10 +216,10 @@ void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction) // Updating the particles is also independent, so we can parallelize // this as well -#pragma omp parallel for private(p) +#pragma omp parallel for for (size_t i = 0; i < size; i++) { - p->r_vec += dt * this->particles[i].v_vec; - p->v_vec += dt * force_res[i] / this->particles[i].m; + this->particles[i].r_vec += dt * this->particles[i].v_vec; + this->particles[i].v_vec += dt * force_res[i] / this->particles[i].m; } this->t += dt; @@ -281,14 +278,14 @@ void PenningTrap::write_simulation_to_dir(std::string path, double time, #pragma omp parallel for private(ofile) for (size_t i = 0; i < this->particles.size(); i++) { ofile.open(path + "particle_" + std::to_string(i) + "_r.txt"); - for (vec_3d &vec : res.r_vecs[i]) { + for (vec3 &vec : res.r_vecs[i]) { ofile << scientific_format(vec(0), 10, 8) << ',' << scientific_format(vec(1), 10, 8) << ',' << scientific_format(vec(2), 10, 8) << '\n'; } ofile.close(); ofile.open(path + "particle_" + std::to_string(i) + "_v.txt"); - for (vec_3d &vec : res.v_vecs[i]) { + for (vec3 &vec : res.v_vecs[i]) { ofile << scientific_format(vec(0), 10, 8) << ',' << scientific_format(vec(1), 10, 8) << ',' << scientific_format(vec(2), 10, 8) << '\n'; diff --git a/src/main.cpp b/src/main.cpp index 8c5b300..f6a707c 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -3,7 +3,7 @@ * @author Cory Alexander Balaton (coryab) * @author Janita Ovidie Sandtrøen Willumsen (janitaws) * - * @version 0.1 + * @version 1.0 * * @brief The main program for this project * @@ -25,16 +25,16 @@ #define N 40000 // Particles used for testing -Particle p1(vec_3d{20., 0., 20.}, vec_3d{0., 25., 0.}); ///< Particle 1 -Particle p2(vec_3d{25., 25., 0.}, vec_3d{0., 40., 5.}); ///< Particle 2 +Particle p1(vec3{20., 0., 20.}, vec3{0., 25., 0.}); ///< Particle 1 +Particle p2(vec3{25., 25., 0.}, vec3{0., 40., 5.}); ///< Particle 2 /** @brief The analytical solution for particle p1 * * @param t Time * - * @return vec_3d + * @return vec3 * */ -vec_3d analytical_solution_particle_1(double t) +vec3 analytical_solution_particle_1(double t) { double w_0 = T / CA_MASS; double w_z2 = (50. * V / 1000.) / (CA_MASS * 500. * 500.); @@ -45,8 +45,7 @@ vec_3d analytical_solution_particle_1(double t) std::complex f = A_p * std::exp(std::complex(0., -w_p * t)) + A_n * std::exp(std::complex(0., -w_n * t)); - vec_3d res{std::real(f), std::imag(f), - 20. * std::cos(std::sqrt(w_z2) * t)}; + vec3 res{std::real(f), std::imag(f), 20. * std::cos(std::sqrt(w_z2) * t)}; return res; } @@ -213,8 +212,7 @@ void potential_resonance_narrow_sweep() double freq_start = 1.; double freq_end = 1.7; double freq_increment = .002; - size_t freq_iterations = - (size_t)((freq_end - freq_start) / freq_increment); + size_t freq_iterations = (size_t)((freq_end - freq_start) / freq_increment); double res[4][freq_iterations]; @@ -270,8 +268,7 @@ void potential_resonance_narrow_sweep_interaction() double freq_start = 1.; double freq_end = 1.7; double freq_increment = .002; - size_t freq_iterations = - (size_t)((freq_end - freq_start) / freq_increment); + size_t freq_iterations = (size_t)((freq_end - freq_start) / freq_increment); double res[4][freq_iterations]; diff --git a/src/scripts/animate_100_particles.py b/src/scripts/animate_100_particles.py index 0fc0083..076b186 100644 --- a/src/scripts/animate_100_particles.py +++ b/src/scripts/animate_100_particles.py @@ -32,7 +32,7 @@ def animate(): fig = plt.figure() ax = fig.add_subplot(projection="3d") - arr = get_data([f"output/simulate_100_particles/particle_{i}.txt" for i in range(100)]) + arr = get_data([f"output/simulate_100_particles/particle_{i}_r.txt" for i in range(100)]) arr = arr[:, :, ::10] diff --git a/src/scripts/plot_particles_left.py b/src/scripts/plot_particles_left.py index 85bfd78..61185d9 100644 --- a/src/scripts/plot_particles_left.py +++ b/src/scripts/plot_particles_left.py @@ -1,30 +1,68 @@ import matplotlib.pyplot as plt +import seaborn as sns + +sns.set_theme() +params = { + "font.family": "Serif", + "font.serif": "Roman", + "text.usetex": True, + "axes.titlesize": "large", + "axes.labelsize": "large", + "xtick.labelsize": "large", + "ytick.labelsize": "large", + "legend.fontsize": "medium" +} +plt.rcParams.update(params) def main(): - with open("output/time_dependent_potential/res.txt") as f: - lines = f.readlines() - x = [] - y1 = [] - y2 = [] - y3 = [] - for line in lines: - l = line.strip().split(",") - x.append(float(l[0])) - y1.append(float(l[1])) - y2.append(float(l[2])) - y3.append(float(l[3])) + colors = [ + "lightskyblue", + "deepskyblue", + "salmon", + "tomato", + "mediumaquamarine", + "mediumseagreen" + ] + files = [ + "output/time_dependent_potential/wide_sweep.txt", + "output/time_dependent_potential/narrow_sweep.txt", + "output/time_dependent_potential/narrow_sweep_interactions.txt", + ] + outputs = [ + "../latex/images/wide_sweep.pdf", + "../latex/images/narrow_sweep.pdf", + "../latex/images/narrow_sweep_interactions.pdf", + ] + for file, output in zip(files, outputs): + with open(file) as f: + lines = f.readlines() + x = [] + y1 = [] + y2 = [] + y3 = [] + for line in lines: + l = line.strip().split(",") + x.append(float(l[0])) + y1.append(float(l[1])) + y2.append(float(l[2])) + y3.append(float(l[3])) - plt.plot(x,y1,label=f"amplitude: 0.1") - plt.plot(x,y2,label=f"amplitude: 0.4") - plt.plot(x,y3,label=f"amplitude: 0.7") + fig, ax = plt.subplots() + ax.plot(x, y1, label=r"$f_{1} = 0.1$", color=colors[0]) + ax.plot(x, y2, label=r"$f_{2} = 0.4$", color=colors[2]) + ax.plot(x, y3, label=r"$f_{3} = 0.7$", color=colors[4]) - plt.xlabel(r"$\omega_V$ (MHz)") - plt.ylabel(r"Fraction of particles left") - plt.title(r"The fraction of particles left in the Penning trap " - "after 500 microseconds for different amplitudes and frequencies") - plt.legend() - # plt.show() - plt.savefig("../latex/images/particles_left.pdf") + ax.set_xlabel(r"Frequency $\omega_V$ (MHz)") + # ax.set_xlim((0, 2.8)) + ax.set_ylabel(r"Fraction of particles left") + # ax.set_ylim((-0.1, 1.1)) + # plt.title(r"The fraction of particles left in the Penning trap " + # "after 500 microseconds for different amplitudes and frequencies") + + ax.legend(loc="upper right") + + # plt.show() + fig.savefig(output, bbox_inches="tight") if __name__ == "__main__": main() diff --git a/src/test_suite.cpp b/src/test_suite.cpp index 3828541..7a1e2ad 100644 --- a/src/test_suite.cpp +++ b/src/test_suite.cpp @@ -3,7 +3,7 @@ * @author Cory Alexander Balaton (coryab) * @author Janita Ovidie Sandtrøen Willumsen (janitaws) * - * @version 0.1 + * @version 1.0 * * @brief The test suite for the project * @@ -25,30 +25,29 @@ class PenningTrapTest public: /** @brief Test that the external E field gives correct values. * */ - static void test_external_E_field() + void test_external_E_field() { PenningTrap trap; // Vector containing inputs and expected results - std::vector> tests; + std::vector> tests; + + tests.push_back(std::make_pair(vec3{0., 0., 0.}, vec3{0., 0., 0.})); tests.push_back( - std::make_pair(vec_3d{0., 0., 0.}, vec_3d{0., 0., 0.})); + std::make_pair(vec3{10., 0., 0.}, vec3{96.4852558, 0., 0.})); tests.push_back( - std::make_pair(vec_3d{10., 0., 0.}, vec_3d{96.4852558, 0., 0.})); + std::make_pair(vec3{10., 0., 0.}, vec3{96.4852558, 0., 0.})); tests.push_back( - std::make_pair(vec_3d{10., 0., 0.}, vec_3d{96.4852558, 0., 0.})); + std::make_pair(vec3{0., 10., 0.}, vec3{0., 96.4852558, 0.})); tests.push_back( - std::make_pair(vec_3d{0., 10., 0.}, vec_3d{0., 96.4852558, 0.})); + std::make_pair(vec3{0., 0., 10.}, vec3{0., 0., -192.9705116})); - tests.push_back( - std::make_pair(vec_3d{0., 0., 10.}, vec_3d{0., 0., -192.9705116})); - - vec_3d result; - vec_3d v; + vec3 result; + vec3 v; std::stringstream msg; for (size_t i = 0; i < tests.size(); i++) { v = tests.at(i).first; @@ -64,38 +63,39 @@ public: /** @brief Test that the external B field gives correct values. * */ - static void test_external_B_field() + void test_external_B_field() { // No point in testing at different points since it's not dependent // on position. PenningTrap trap; - vec_3d expected{0., 0., T}; - vec_3d result = trap.external_B_field(vec_3d{0., 0., 0.}); + vec3 expected{0., 0., T}; + vec3 result = trap.external_B_field(vec3{0., 0., 0.}); ASSERT(close_to(expected, result), "Testing the external B field at (0,0,0)"); } - /** @brief Test that the force between particles gives expected results. + /** @brief Test that the force between particles gives expected + * results. * */ - static void test_force_on_particle() + void test_force_on_particle() { PenningTrap trap; - vec_3d v{0., 0., 0.}; + vec3 v{0., 0., 0.}; // Add particles to test - trap.add_particle(Particle(vec_3d{0., 0., 0.}, v)); - trap.add_particle(Particle(vec_3d{1., 0., 0.}, v)); - trap.add_particle(Particle(vec_3d{0., 3., 4.}, v)); + trap.add_particle(Particle(vec3{0., 0., 0.}, v)); + trap.add_particle(Particle(vec3{1., 0., 0.}, v)); + trap.add_particle(Particle(vec3{0., 3., 4.}, v)); // Test p0 and p1 - vec_3d expected{-1., 0., 0.}; - vec_3d result = trap.force_on_particle(0, 1); + vec3 expected{-1., 0., 0.}; + vec3 result = trap.force_on_particle(0, 1); ASSERT(close_to(expected, result), "Testing the force on a particle at (0,0,0) from a " "particle at (1,0,0)."); // Test p0 and p2 - expected = vec_3d{0, -.024, -.032}; + expected = vec3{0, -.024, -.032}; result = trap.force_on_particle(0, 2); ASSERT(close_to(expected, result), "Testing the force on a particle at (0,0,0) from a " @@ -104,37 +104,37 @@ public: /** @brief Test that the total external force returns expected results * */ - static void test_total_force_external() + void test_total_force_external() { PenningTrap trap; - trap.add_particle(Particle(vec_3d{1., 2., 3.}, vec_3d{3., 4., 5.})); + trap.add_particle(Particle(vec3{1., 2., 3.}, vec3{3., 4., 5.})); - vec_3d expected{395.58954878, -270.15871624, -57.89115348}; - vec_3d result = trap.total_force_external(0); + vec3 expected{395.58954878, -270.15871624, -57.89115348}; + vec3 result = trap.total_force_external(0); ASSERT(close_to(expected, result), "Testing the total external force on a particle at " "(1,2,3) with velocity (3,4,5)"); } - /** @brief Test that the total force of all particles on a single particle - * returns expected results. + /** @brief Test that the total force of all particles on a single + * particle returns expected results. * */ - static void test_total_force_particles() + void test_total_force_particles() { PenningTrap trap; - trap.add_particle(Particle(vec_3d{0., 0., 0.}, vec_3d{0., 0., 0.})); + trap.add_particle(Particle(vec3{0., 0., 0.}, vec3{0., 0., 0.})); - vec_3d expected{0., 0., 0.}; - vec_3d result = trap.total_force_particles(0); + vec3 expected{0., 0., 0.}; + vec3 result = trap.total_force_particles(0); ASSERT(close_to(expected, result), "Testing the total force of all particles on particle 0 " "with only a single particle"); - trap.add_particle(Particle(vec_3d{1., 0., 0.}, vec_3d{0., 0., 0.})); - trap.add_particle(Particle(vec_3d{0., 1., 0.}, vec_3d{0., 0., 0.})); - trap.add_particle(Particle(vec_3d{0., 0., 1.}, vec_3d{0., 0., 0.})); + trap.add_particle(Particle(vec3{1., 0., 0.}, vec3{0., 0., 0.})); + trap.add_particle(Particle(vec3{0., 1., 0.}, vec3{0., 0., 0.})); + trap.add_particle(Particle(vec3{0., 0., 1.}, vec3{0., 0., 0.})); - expected = vec_3d().fill(-138935.333); + expected = vec3().fill(-138935.333); result = trap.total_force_particles(0); ASSERT(close_to(expected, result), "Testing the total force of all particles on particle 0 " @@ -144,10 +144,11 @@ public: int main() { - PenningTrapTest::test_external_E_field(); - PenningTrapTest::test_external_B_field(); - PenningTrapTest::test_force_on_particle(); - PenningTrapTest::test_total_force_external(); - PenningTrapTest::test_total_force_particles(); + PenningTrapTest test; + test.test_external_E_field(); + test.test_external_B_field(); + test.test_force_on_particle(); + test.test_total_force_external(); + test.test_total_force_particles(); return 0; } diff --git a/src/utils.cpp b/src/utils.cpp index a7a1507..ca2d50d 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -19,8 +19,7 @@ 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) {