include/Particle.hpp

@@ -15,15 +15,16 @@
 #include <armadillo>
 
 #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.
      * */

src/Particle.cpp

@@ -12,13 +12,11 @@
 
 #include "Particle.hpp"
 
-Particle::Particle(double q, double m,
+Particle::Particle(vec_3d r_vec, vec_3d v_vec, double q, double m)
-                   vec_3d r_vec,
-                   vec_3d v_vec)
 {
     // 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;
 }

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

src/main.cpp

@@ -18,18 +18,15 @@
 #include <vector>
 
 #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.},
+Particle p1(vec_3d{20., 0., 20.}, vec_3d{0., 25., 0.}); ///< Particle 1
-            vec_3d{0., 25., 0.}); ///< Particle 1
+Particle p2(vec_3d{25., 25., 0.}, vec_3d{0., 40., 5.}); ///< Particle 2
-Particle p2(CHARGE, MASS, 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_0 = T / CA_MASS;
-    double w_z2 = (50. * V / 1000.) / (MASS * 500. * 500.);
+    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(
+                trap.reinitialize(amplitudes[i], res[0][j]);
-                    [](double f, double r, double t) {
-                        return (25. * V / 1000.) * (1. + f * std::cos(r * t));
-                    },
-                    amplitudes[i], res[0][j], std::placeholders::_1));
                 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(
+                trap.reinitialize(amplitudes[i], res[0][j]);
-                    [](double f, double r, double t) {
-                        return (25. * V / 1000.) * (1. + f * std::cos(r * t));
-                    },
-                    amplitudes[i], res[0][j], std::placeholders::_1));
                 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(
+                trap.reinitialize(amplitudes[i], res[0][j]);
-                    [](double f, double r, double t) {
-                        return (25. * V / 1000.) * (1. + f * std::cos(r * t));
-                    },
-                    amplitudes[i], res[0][j], std::placeholders::_1));
                 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();
 

src/test_suite.cpp

@@ -15,6 +15,7 @@
 #include <string>
 
 #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(vec_3d{0., 0., 0.}, v));
-        trap.add_particle(Particle(1., 40., vec_3d{1., 0., 0.}, v));
+        trap.add_particle(Particle(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., 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 "