From 8f315a095598cba6fc613d591aabe14461755d20 Mon Sep 17 00:00:00 2001 From: Cory Date: Mon, 23 Oct 2023 09:54:13 +0200 Subject: [PATCH] Change Particle - Give Particle default values for mass and charge - Modify PenningTrap to use the new Particle constructor - Modify main.cpp to use the new Particle constructor - Modify test_suite.cpp to use the new Particle constructor --- include/Particle.hpp | 7 ++++--- src/Particle.cpp | 8 +++----- src/PenningTrap.cpp | 2 +- src/main.cpp | 45 ++++++++++++++++++-------------------------- src/test_suite.cpp | 19 ++++++++++--------- 5 files changed, 36 insertions(+), 45 deletions(-) diff --git a/include/Particle.hpp b/include/Particle.hpp index 4b5e716..fb26a78 100644 --- a/include/Particle.hpp +++ b/include/Particle.hpp @@ -15,15 +15,16 @@ #include #include "typedefs.hpp" +#include "constants.hpp" /** @brief A class that holds attributes of a particle * */ class Particle { private: - double q; ///< Charge - double m; ///< Mass vec_3d r_vec; ///< position vec_3d v_vec; ///< velocity + double q; ///< Charge + double m; ///< Mass public: /** @brief Initialize the particle. @@ -36,7 +37,7 @@ public: * @param r_vec The initial position of the particle * @param v_vec The initial velocity of the particle * */ - Particle(double q, double m, vec_3d r_vec, vec_3d v_vec); + Particle(vec_3d r_vec, vec_3d v_vec, double q = CA_CHARGE, double m = CA_MASS); /** @brief Make private attributes available for PenningTrap. * */ diff --git a/src/Particle.cpp b/src/Particle.cpp index 374812f..e1f837b 100644 --- a/src/Particle.cpp +++ b/src/Particle.cpp @@ -12,13 +12,11 @@ #include "Particle.hpp" -Particle::Particle(double q, double m, - vec_3d r_vec, - vec_3d v_vec) +Particle::Particle(vec_3d r_vec, vec_3d v_vec, double q, double m) { // Giving the particle its properties - this->q = q; - this->m = m; this->r_vec = r_vec; this->v_vec = v_vec; + this->q = q; + this->m = m; } diff --git a/src/PenningTrap.cpp b/src/PenningTrap.cpp index 806901b..6ccf12d 100644 --- a/src/PenningTrap.cpp +++ b/src/PenningTrap.cpp @@ -31,7 +31,7 @@ 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(1., 40., vec_3d(vec_3d().randn() * .1 * this->d), + Particle(vec_3d(vec_3d().randn() * .1 * this->d), vec_3d(vec_3d().randn() * .1 * this->d))); } } diff --git a/src/main.cpp b/src/main.cpp index 641c66a..76286d1 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -18,18 +18,15 @@ #include #include "PenningTrap.hpp" +#include "constants.hpp" #include "utils.hpp" #define PARTICLES 100 #define N 40000 -#define CHARGE 1. // unit: e -#define MASS 40.078 // unit: amu // Particles used for testing -Particle p1(CHARGE, MASS, vec_3d{20., 0., 20.}, - vec_3d{0., 25., 0.}); ///< Particle 1 -Particle p2(CHARGE, MASS, vec_3d{25., 25., 0.}, - vec_3d{0., 40., 5.}); ///< Particle 2 +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 /** @brief The analytical solution for particle p1 * @@ -39,8 +36,8 @@ Particle p2(CHARGE, MASS, vec_3d{25., 25., 0.}, * */ vec_3d analytical_solution_particle_1(double t) { - double w_0 = T / MASS; - double w_z2 = (50. * V / 1000.) / (MASS * 500. * 500.); + double w_0 = T / CA_MASS; + double w_z2 = (50. * V / 1000.) / (CA_MASS * 500. * 500.); double w_p = (w_0 + std::sqrt(w_0 * w_0 - 2. * w_z2)) / 2.; double w_n = (w_0 - std::sqrt(w_0 * w_0 - 2. * w_z2)) / 2.; double A_p = (25. + w_n * 20.) / (w_n - w_p); @@ -184,11 +181,7 @@ void potential_resonance_wide_sweep() for (size_t i = 0; i < 3; i++) { for (size_t j = 0; j < freq_iterations; j++) { // Reset particles and give new time dependent potential. - trap.reinitialize(std::bind( - [](double f, double r, double t) { - return (25. * V / 1000.) * (1. + f * std::cos(r * t)); - }, - amplitudes[i], res[0][j], std::placeholders::_1)); + trap.reinitialize(amplitudes[i], res[0][j]); res[i + 1][j] = trap.fraction_of_particles_left(time, N, "rk4", false); } @@ -211,7 +204,7 @@ void potential_resonance_wide_sweep() * MHz. * * */ -void potential_resonance_no_interaction_narrow_sweep() +void potential_resonance_narrow_sweep() { double time = 500.; @@ -230,6 +223,12 @@ void potential_resonance_no_interaction_narrow_sweep() std::ofstream ofile; +#pragma omp parallel for + // Insert frequencies + for (size_t i = 0; i < freq_iterations; i++) { + res[0][i] = freq_start + freq_increment * i; + } + #pragma omp parallel { // Each thread creates a PenningTrap instance and reuses it throughout @@ -239,11 +238,7 @@ void potential_resonance_no_interaction_narrow_sweep() for (size_t i = 0; i < 3; i++) { for (size_t j = 0; j < freq_iterations; j++) { // Reset particles and give new time dependent potential. - trap.reinitialize(std::bind( - [](double f, double r, double t) { - return (25. * V / 1000.) * (1. + f * std::cos(r * t)); - }, - amplitudes[i], res[0][j], std::placeholders::_1)); + trap.reinitialize(amplitudes[i], res[0][j]); res[i + 1][j] = trap.fraction_of_particles_left(time, N, "rk4", false); } @@ -266,7 +261,7 @@ void potential_resonance_no_interaction_narrow_sweep() * MHz. * * */ -void potential_resonance_with_interaction_narrow_sweep() +void potential_resonance_narrow_sweep_interaction() { double time = 500.; @@ -299,11 +294,7 @@ void potential_resonance_with_interaction_narrow_sweep() for (size_t i = 0; i < 3; i++) { for (size_t j = 0; j < freq_iterations; j++) { // Reset particles and give new time dependent potential. - trap.reinitialize(std::bind( - [](double f, double r, double t) { - return (25. * V / 1000.) * (1. + f * std::cos(r * t)); - }, - amplitudes[i], res[0][j], std::placeholders::_1)); + trap.reinitialize(amplitudes[i], res[0][j]); res[i + 1][j] = trap.fraction_of_particles_left(time, N); } } @@ -354,7 +345,7 @@ int main() t1 = omp_get_wtime(); - potential_resonance_no_interaction_narrow_sweep(); + potential_resonance_narrow_sweep(); t2 = omp_get_wtime(); @@ -363,7 +354,7 @@ int main() t1 = omp_get_wtime(); - potential_resonance_with_interaction_narrow_sweep(); + potential_resonance_narrow_sweep_interaction(); t2 = omp_get_wtime(); diff --git a/src/test_suite.cpp b/src/test_suite.cpp index a2a0639..16c3f53 100644 --- a/src/test_suite.cpp +++ b/src/test_suite.cpp @@ -15,6 +15,7 @@ #include #include "PenningTrap.hpp" +#include "constants.hpp" #include "utils.hpp" /** @brief Test class for the Penning trap. @@ -81,9 +82,9 @@ public: vec_3d v{0., 0., 0.}; // Add particles to test - trap.add_particle(Particle(1., 40., vec_3d{0., 0., 0.}, v)); - trap.add_particle(Particle(1., 40., vec_3d{1., 0., 0.}, v)); - trap.add_particle(Particle(1., 40., vec_3d{0., 3., 4.}, v)); + 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)); // Test p0 and p1 vec_3d expected{-1., 0., 0.}; @@ -106,7 +107,7 @@ public: { PenningTrap trap; trap.add_particle( - Particle(1., 40., vec_3d{1., 2., 3.}, vec_3d{3., 4., 5.})); + Particle(vec_3d{1., 2., 3.}, vec_3d{3., 4., 5.})); vec_3d expected{395.58954878, -270.15871624, -57.89115348}; vec_3d result = trap.total_force_external(0); @@ -122,7 +123,7 @@ public: { PenningTrap trap; trap.add_particle( - Particle(1., 40., vec_3d{0., 0., 0.}, vec_3d{0., 0., 0.})); + Particle(vec_3d{0., 0., 0.}, vec_3d{0., 0., 0.})); vec_3d expected{0., 0., 0.}; vec_3d result = trap.total_force_particles(0); @@ -131,13 +132,13 @@ public: "with only a single particle"); trap.add_particle( - Particle(1., 40., vec_3d{1., 0., 0.}, vec_3d{0., 0., 0.})); + Particle(vec_3d{1., 0., 0.}, vec_3d{0., 0., 0.})); trap.add_particle( - Particle(1., 40., vec_3d{0., 1., 0.}, vec_3d{0., 0., 0.})); + Particle(vec_3d{0., 1., 0.}, vec_3d{0., 0., 0.})); trap.add_particle( - Particle(1., 40., vec_3d{0., 0., 1.}, vec_3d{0., 0., 0.})); + Particle(vec_3d{0., 0., 1.}, vec_3d{0., 0., 0.})); - expected = vec_3d().fill(-3473.383325); + expected = vec_3d().fill(-138935.333); result = trap.total_force_particles(0); ASSERT(close_to(expected, result), "Testing the total force of all particles on particle 0 "