src/PenningTrap.cpp

@@ -52,8 +52,8 @@ vec_3d PenningTrap::v_func(unsigned int i, unsigned int j, double dt)
         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]));
+                      * (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();
@@ -71,8 +71,8 @@ vec_3d PenningTrap::r_func(unsigned int i, unsigned int j, double dt)
         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]));
+                      * (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();
@@ -116,8 +116,8 @@ vec_3d PenningTrap::total_force_external(unsigned int i)
     }
 
     return vec_3d(p.q
-           * (this->external_E_field(p.r_vec)
-              + arma::cross(p.v_vec, this->external_B_field(p.r_vec))));
+                  * (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(unsigned int i)
@@ -134,11 +134,10 @@ vec_3d PenningTrap::total_force_particles(unsigned int i)
 
 vec_3d PenningTrap::total_force(unsigned int i)
 {
-    if (arma::norm(this->particles[i].r_vec) > this->d) {
-        return vec_3d{0., 0., 0.};
-    }
-    return vec_3d(this->total_force_external(i)
-                  - this->total_force_particles(i));
+    return arma::norm(this->particles[i].r_vec) > this->d
+               ? vec_3d{0., 0., 0.}
+               : vec_3d(this->total_force_external(i)
+                        - this->total_force_particles(i));
 }
 
 void PenningTrap::evolve_RK4(double dt, bool particle_interaction)

src/main.cpp

@@ -97,6 +97,7 @@ void simulate_single_particle_with_different_steps()
     // Calculate relative error for RK4
     std::string path = "output/relative_error/RK4/";
     mkpath(path);
+#pragma omp parallel for
     for (int i = 0; i < 4; i++) {
         int steps = 4000 * std::pow(2, i);
         double dt = time / (double)steps;
@@ -114,6 +115,7 @@ void simulate_single_particle_with_different_steps()
     // Calculate relative error for forward Euler
     path = "output/relative_error/euler/";
     mkpath(path);
+#pragma omp parallel for
     for (int i = 0; i < 4; i++) {
         int steps = 4000 * std::pow(2, i);
         double dt = time / (double)steps;
@@ -137,8 +139,9 @@ void simulate_100_particles()
 
     double time = 50.; // microseconds
 
-     trap.write_simulation_to_dir("output/simulate_100_particles", time, N,
-     "rk4", false);
+    //trap.write_simulation_to_dir("output/simulate_100_particles", time, N,
+                                 //"rk4", false);
+    trap.simulate(time, N);
 }
 
 /** @brief Simulate 100 particles over 500 \f$ \mu s \f$ using a time
@@ -148,7 +151,7 @@ void simulate_100_particles()
  * MHz.
  *
  * */
-void simulate_100_particles_with_time_potential()
+void simulate_100_particles_with_time_potential_wide_sweep()
 {
     double time = 500.;
 
@@ -167,24 +170,85 @@ void simulate_100_particles_with_time_potential()
 
     std::ofstream ofile;
 
-    #pragma omp parallel for
+#pragma omp parallel for
     for (size_t i = 0; i < freq_iterations; i++) {
-        res[0][i] = freq_start+ freq_increment*i;
+        res[0][i] = freq_start + freq_increment * i;
     }
 
-    // Using num_threads() is usually bad practice, but not having it
-    // causes a SIGKILL.
-    #pragma omp parallel for collapse(2) num_threads(4)
+// Using num_threads() is usually bad practice, but not having it
+// causes a SIGKILL.
+#pragma omp parallel for collapse(2) num_threads(4)
     for (size_t i = 0; i < 3; i++) {
         for (size_t j = 0; j < freq_iterations; j++) {
-            res[i+1][j] = PenningTrap(
-                (unsigned)100, T,
-                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),
-                500., 0.).fraction_of_particles_left(time, N, "rk4", false);
+            res[i + 1][j]
+                = PenningTrap((unsigned)100, T,
+                              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),
+                              500., 0.)
+                      .fraction_of_particles_left(time, N, "rk4", false);
+        }
+    }
+
+    ofile.open(path + "res.txt");
+    for (size_t i = 0; i < freq_iterations; i++) {
+        ofile << res[0][i] << "," << res[1][i] << "," << res[2][i] << ","
+              << res[3][i] << "\n";
+    }
+    ofile.close();
+}
+
+/** @brief Simulate 100 particles over 500 \f$ \mu s \f$ using a time
+ *  dependent potential.
+ *
+ *  @details The simulation sweeps over different frequencies in [1., 1.7]
+ * MHz.
+ *
+ * */
+void simulate_100_particles_with_time_potential_narrow_sweep()
+{
+    double time = 500.;
+
+    double amplitudes[]{.1, .4, .7};
+
+    double freq_start = 1.;
+    double freq_end = 1.7;
+    double freq_increment = .001;
+    size_t freq_iterations
+        = (size_t)((freq_end - freq_start) / freq_increment);
+
+    double res[4][freq_iterations];
+
+    std::string path = "output/time_dependent_potential/";
+    mkpath(path);
+
+    std::ofstream ofile;
+
+#pragma omp parallel for
+    for (size_t i = 0; i < freq_iterations; i++) {
+        res[0][i] = freq_start + freq_increment * i;
+    }
+
+// Using num_threads() is usually bad practice, but not having it
+// causes a SIGKILL.
+#pragma omp parallel for collapse(2) num_threads(4)
+    for (size_t i = 0; i < 3; i++) {
+        for (size_t j = 0; j < freq_iterations; j++) {
+            res[i + 1][j]
+                = PenningTrap((unsigned)100, T,
+                              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),
+                              500., 0.)
+                      .fraction_of_particles_left(time, N, "rk4", false);
         }
     }
 
@@ -200,15 +264,16 @@ int main()
 {
     double t0 = omp_get_wtime();
 
-    simulate_single_particle();
+    //simulate_single_particle();
 
-    simulate_two_particles();
+    //simulate_two_particles();
 
-    simulate_single_particle_with_different_steps();
+    //simulate_single_particle_with_different_steps();
 
     simulate_100_particles();
 
-    simulate_100_particles_with_time_potential();
+    //simulate_100_particles_with_time_potential_wide_sweep();
+    //simulate_100_particles_with_time_potential_narrow_sweep();
 
     double end = omp_get_wtime();