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;

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<Particle> 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,7 +70,7 @@ 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.)
+        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:
@@ -79,7 +79,7 @@ vec_3d PenningTrap::v_func(uint i, uint j, double dt)
     }
 }
 
-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,7 +89,7 @@ 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.)
+        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:
@@ -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
+    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,8 +166,8 @@ void PenningTrap::evolve_RK4(double dt, bool particle_interaction)
     std::vector<Particle> original_particles = this->particles;
     std::vector<Particle> tmp_particles = this->particles;
 
-    vec_3d (PenningTrap::*force)(uint) =
-        particle_interaction ? &PenningTrap::total_force
+    vec3 (PenningTrap::*force)(uint) = particle_interaction
+                                           ? &PenningTrap::total_force
                                            : &PenningTrap::total_force_external;
 
     size_t size = this->particles.size();
@@ -203,11 +201,10 @@ 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
+    vec3 (PenningTrap::*force)(uint) = particle_interaction
+                                           ? &PenningTrap::total_force
                                            : &PenningTrap::total_force_external;
 
     // Calculating the force for each particle is independent and therefore
@@ -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';

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<double> f =
         A_p * std::exp(std::complex<double>(0., -w_p * t))
         + A_n * std::exp(std::complex<double>(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];
 

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]
 

src/scripts/plot_particles_left.py

@@ -1,7 +1,40 @@
 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:
+    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 = []
@@ -14,17 +47,22 @@ def main():
                 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])
+
+        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.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")
+        fig.savefig(output, bbox_inches="tight")
 
 if __name__ == "__main__":
     main()

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<std::pair<vec_3d, vec_3d>> tests;
+        std::vector<std::pair<vec3, vec3>> 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;
 }

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<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) {