Update documentation

include/Particle.hpp

@@ -30,6 +30,11 @@ public:
      *
      *  @details Initialize the particle with a charge, mass, position and
      *  velocity.
+     *
+     *  @param q The charge of the particle
+     *  @param m The mass of the particle
+     *  @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);
 

include/PenningTrap.hpp

@@ -19,7 +19,6 @@
 #include "constants.hpp"
 #include "typedefs.hpp"
 
-
 #pragma omp declare reduction(+ : vec_3d : omp_out += omp_in)                  \
     initializer(omp_priv = omp_orig)
 
@@ -31,35 +30,100 @@
 class PenningTrap {
 private:
     double B_0;                        ///< Magnetic field strength
-    double V_0;                      ///< Applied potential
+    std::function<double(double)> V_0; ///< Applied potential
     double d;                          ///< Characteristic dimension
+    double t;                          ///< Current time
     std::vector<Particle> particles;   ///< The particles in the Penning trap
-    sim_arr k_v;
+    sim_arr k_v; ///< A 2D vector containing all \f$k_{i,j}\f$ where \f$j\f$ is
-    sim_arr k_r;
+                 ///< the index of a particle
+    sim_arr k_r; ///< A 2D vector containing all \f$k_{i,j}\f$ where \f$j\f$ is
+                 ///< the index of a particle
 
-    vec_3d v_func(int, int, double);
+    /** @brief Helper for evolve_RK4 when calculating \f$k_{v,i,j}\f$ values
-    vec_3d r_func(int, int, double);
+     *
+     *  @details Something
+     *
+     *  @param i Index i for \f$k_{v,i,j}\f$
+     *  @param j Index j for \f$k_{v,i,j}\f$
+     *  @param dt the step length (delta time)
+     *
+     *  @return vec_3d
+     * */
+    vec_3d v_func(unsigned int i, unsigned int j, double dt);
+
+    /** @brief Helper for evolve_RK4 when calculating \f$k_{r,i,j}\f$ values
+     *
+     *  @details Something
+     *
+     *  @param i Index i for \f$k_{r,i,j}\f$
+     *  @param j Index j for \f$k_{r,i,j}\f$
+     *  @param dt the step length (delta time)
+     *
+     *  @return vec_3d
+     * */
+    vec_3d r_func(unsigned int i, unsigned int j, double dt);
 
 public:
-    /** @brief Set B_0, V_0 and d.
+    /** @brief Constructor for the PenningTrap class
+     *
+     *  @param B_0 The magnetic field strength
+     *  @param V_0 The time dependent applied potential
+     *  @param d The characteristic dimension
+     *  @param t The starting time
      * */
-    PenningTrap(double B_0 = T, double V_0 = 25. * V / 1000., double d = 500.);
+    PenningTrap(
+        double B_0 = T,
+        std::function<double(double)> V_0 =
+            [](double t) { return 25. * V / 1000.; },
+        double d = 500., double t = 0.);
 
-    PenningTrap(int i, double B_0 = T, double V_0 = 25. * V / 1000.,
+    /** @brief Constructor for the PenningTrap class
-                double d = 500.);
+     *
+     *  @param i The number of particles to generate
+     *  @param B_0 The magnetic field strength
+     *  @param V_0 The time dependent applied potential
+     *  @param d The characteristic dimension
+     *  @param t The starting time
+     * */
+    PenningTrap(
+        unsigned int i, double B_0 = T,
+        std::function<double(double)> V_0 =
+            [](double t) { return 25. * V / 1000.; },
+        double d = 500., double t = 0.);
 
-    PenningTrap(std::vector<Particle> particles, double B_0 = T,
+    /** @brief Constructor for the PenningTrap class
-                double V_0 = 25. * V / 1000., double d = 500.);
+     *
+     *  @param particles The starting particles
+     *  @param B_0 The magnetic field strength
+     *  @param V_0 The time dependent applied potential
+     *  @param d The characteristic dimension
+     *  @param t The starting time
+     * */
+    PenningTrap(
+        std::vector<Particle> particles, double B_0 = T,
+        std::function<double(double)> V_0 =
+            [](double t) { return 25. * V / 1000.; },
+        double d = 500., double t = 0.);
 
     /** @brief Add a particle to the system
+     *
+     *  @param particle The particle to add to the Penning trap
      * */
     void add_particle(Particle particle);
 
     /** @brief Calculate E at point r
+     *
+     *  @param r The position where we want to calculate the E field
+     *
+     *  @return vec_3d
      * */
     vec_3d external_E_field(vec_3d r);
 
     /** @brief Calculate B at point r
+     *
+     *  @param r The position where we want to calculate the B field
+     *
+     *  @return vec_3d
      * */
     vec_3d external_B_field(vec_3d r);
 
@@ -67,37 +131,90 @@ public:
      *
      *  @details Calculate the force exhibited on particle p_i from
      *  particle p_j.
+     *
+     *  @param i The index of particle p_i
+     *  @param j The index of particle p_j
+     *
+     *  @return vec_3d
      * */
-    vec_3d force_on_particle(int i, int j);
+    vec_3d force_on_particle(unsigned int i, unsigned int j);
 
     /** @brief Calculate the total external force on a particle.
      *
      *  @details Calculate the total amount of force that E and B exhibits
      *  on particle p_i.
+     *
+     *  @param i The index of particle p_i
+     *
+     *  @return vec_3d
      * */
-    vec_3d total_force_external(int i);
+    vec_3d total_force_external(unsigned int i);
 
-    /** @brief Calculate the total force on a particle from other particles.
+    /** @brief Calculate the total force on a particle p_i from other particles.
+     *
+     *  @param i The index of particle p_i
+     *
+     *  @return vec_3d
      * */
     vec_3d total_force_particles(unsigned int i);
 
-    /** @brief calculate the total force on a particle.
+    /** @brief calculate the total force on a particle p_i.
+     *
+     *  @param i The index of particle p_i
+     *
+     *  @return vec_3d
      * */
-    vec_3d total_force(int i);
+    vec_3d total_force(unsigned int i);
 
     /** @brief Go forward one timestep using the RK4 method
+     *
+     *  @param dt The step length
+     *  @param particle_interaction Turn particle interactions on/off
      * */
     void evolve_RK4(double dt, bool particle_interaction = true);
 
     /** @brief Go forward one timestep using the forward Euler method
+     *
+     *  @param dt The step length
+     *  @param particle_interaction Turn particle interactions on/off
      * */
     void evolve_forward_euler(double dt, bool particle_interaction = true);
 
+    /** @brief Simulate the particle system inside the Penning trap over
+     *  a certain amount of time.
+     *
+     *  @param time The time to simulate in microseconds
+     *  @param steps The amount of steps for the whole simulation
+     *  @param method The method to use when moving forward a timestep
+     *  @param particle_interaction Turn particle interactions on/off
+     * */
     sim_arr simulate(double time, unsigned int steps,
                      std::string method = "rk4",
                      bool particle_interaction = true);
 
-    void write_simulation_to_dir(std::string path, double time, int steps,
+    /** @brief Simulate and write the displacement of all particles to files.
+     *
+     *  @param path The directory to save the data
+     *  @param time The time to simulate in microseconds
+     *  @param steps The amount of steps for the whole simulation
+     *  @param method The method to use when moving forward a timestep
+     *  @param particle_interaction Turn particle interactions on/off
+     * */
+    void write_simulation_to_dir(std::string path, double time,
+                                 unsigned int steps, std::string method = "rk4",
+                                 bool particle_interaction = true);
+
+    /** @brief Simulate and calculate what fraction of particles are still
+     *  left inside the Penning trap after the simulation.
+     *
+     *  @param time The time to simulate in microseconds
+     *  @param steps The amount of steps for the whole simulation
+     *  @param method The method to use when moving forward a timestep
+     *  @param particle_interaction Turn particle interactions on/off
+     *
+     *  @return double
+     * */
+    double fraction_of_particles_left(double time, unsigned int steps,
                                       std::string method = "rk4",
                                       bool particle_interaction = true);
 };

include/typedefs.hpp

@@ -1,12 +1,38 @@
+/** @file typedefs.hpp
+ *
+ *  @author Cory Alexander Balaton (coryab)
+ *  @author Janita Ovidie Sandtrøen Willumsen (janitaws)
+ *
+ *  @version 1.0
+ *
+ *  @brief Useful typedefs for cleaner code.
+ *
+ *  @details These typedefs make the code more readable and easy to follow 
+ *  along.
+ *
+ *  @bug No known bugs
+ * */
 #ifndef __TYPEDEFS__
 #define __TYPEDEFS__
 
 #include <vector>
 #include <armadillo>
 
+/** @brief Typedef for the column of the result vector from simulating 
+ *  particles.
+ * */
 typedef std::vector<arma::vec::fixed<3>> sim_cols;
+
+/** @brief Typedef for the row of the result vector from simulating particles.
+ * */
 typedef std::vector<arma::vec::fixed<3>> sim_rows;
+
+/** @brief Typedef for the result of the simulate method.
+ * */
 typedef std::vector<sim_cols> sim_arr;
+
+/** @brief Typedef for a fixed 3d arma vector.
+ * */
 typedef arma::vec::fixed<3> vec_3d;
 
 #endif

include/utils.hpp

@@ -45,40 +45,44 @@
 #define ASSERT(expr, msg) m_assert(expr, #expr, __METHOD_NAME__, __FILE__, \
                                    __LINE__, msg)
 
+/** @def __METHOD_NAME__
+ *  @brief Get the name of the current method/function.
+ * */
 #define __METHOD_NAME__ methodName(__PRETTY_FUNCTION__)
 
 
-/** Code stolen from https://github.com/anderkve/FYS3150
- *  Header: https://github.com/anderkve/FYS3150/blob/master/code_examples/compilation_linking/example_1/include/utils.hpp
- *  Source: https://github.com/anderkve/FYS3150/blob/master/code_examples/compilation_linking/example_1/src/utils.cpp
- * */
-
 /** @brief Turns a double into a string written in scientific format.
+ *
+ *  @details The code is stolen from https://github.com/anderkve/FYS3150.
  *
  *  @param d The number to stringify
  *  @param width The reserved width of the string
  *  @param prec The precision of the stringified number
  *
- *  @return String
+ *  @return std::string
  * */
 std::string scientific_format(double d, int width=20, int prec=10);
 
+
 /** @brief Turns a vector of doubles into a string written in scientific format.
+ *  
+ *  @details The code is stolen from https://github.com/anderkve/FYS3150.
  *
  *  @param v The vector to stringify
  *  @param width The reserved width of the string
  *  @param prec The precision of the stringified number
  *
- *  @return String
+ *  @return std::string
  * */
 std::string scientific_format(const std::vector<double>& v, 
                               int width=20, 
                               int prec=10);
 
+
 /** @brief Test an expression, confirm that test is ok, or abort execution.
  *
- *  This function takes in an expression and prints an OK message if it's
+ *  @details This function takes in an expression and prints an OK message if 
- *  true, or it prints a fail message and aborts execution if it fails.
+ *  it's true, or it prints a fail message and aborts execution if it fails.
  *
  *  @param expr The expression to be evaluated
  *  @param expr_str The stringified version of the expression
@@ -97,27 +101,48 @@ void m_assert(bool expr,
 
 /** @brief Test if two armadillo vectors are close to each other.
  *
- *  This function takes in 2 vectors and checks if they are approximately
+ *  @details This function takes in 2 vectors and checks if they are 
- *  equal to each other given a tolerance.
+ *  approximately equal to each other given a tolerance.
  *
  *  @param a Vector a
  *  @param b Vector b
  *  @param tol The tolerance
  *
- *  @return Boolean
+ *  @return bool
  * */
 bool arma_vector_close_to(arma::vec &a, arma::vec &b, double tol=1e-8);
 
 
-static inline std::string methodName(const std::string& prettyFunction)
+/** @brief Takes in the __PRETTY_FUNCTION__ string and removes the return type.
+ *
+ *  @details This function should only be used for the __METHOD_NAME__ macro,
+ *  since it takes the output from __PRETTY_FUNCTION__ and strips the return
+ *  type.
+ *
+ *  @param pretty_function The string from __PRETTY_FUNCTION__
+ *
+ *  @return std::string
+ * */
+static inline std::string methodName(const std::string& pretty_function)
 {
-    size_t colons = prettyFunction.find("::");
+    size_t colons = pretty_function.find("::");
-    size_t begin = prettyFunction.substr(0,colons).rfind(" ") + 1;
+    size_t begin = pretty_function.substr(0,colons).rfind(" ") + 1;
-    size_t end = prettyFunction.rfind("(") - begin;
+    size_t end = pretty_function.rfind("(") - begin;
 
-    return prettyFunction.substr(begin,end) + "()";
+    return pretty_function.substr(begin,end) + "()";
 }
 
+
+/** @brief Make path given.
+ *
+ *  @details This tries to be the equivalent to "mkdir -p" and creates a new
+ *  directory whenever it needs to.
+ *
+ *  @param path The path to be created
+ *  @param mode The mode/permissions for all the new directories
+ *
+ *  @return bool
+ * */
 bool mkpath(std::string path, int mode = 0777);
 
 #endif

src/Makefile

@@ -8,7 +8,7 @@ CLASSOBJS=$(CLASSSRCS:.cpp=.o)
 
 INCLUDE=../include
 
-CFLAGS=-Wall -larmadillo -llapack -std=c++11 -O3
+CFLAGS=-Wall -larmadillo -lblas -llapack -std=c++11 -O3
 OPENMP=-fopenmp
 
 # Add a debug flag when compiling (For the DEBUG macro in utils.hpp)

src/Particle.cpp

@@ -8,7 +8,6 @@
  *  @brief The implementation of the Particle class.
  *
  *  @bug No known bugs
- *
  * */
 
 #include "Particle.hpp"

src/PenningTrap.cpp

@@ -8,33 +8,28 @@
  *  @brief The implementation of the PenningTrap class.
  *
  *  @bug No known bugs
- *
- *  @todo Implement evolve_RK4
- *  @todo Implement evolve_forward_euler
  * */
 
 #include "PenningTrap.hpp"
 #include "constants.hpp"
 #include "typedefs.hpp"
 #include "utils.hpp"
-#include <cstdlib>
-#include <functional>
-#include <string>
-#include <vector>
 
-PenningTrap::PenningTrap(double B_0, double V_0, double d)
+PenningTrap::PenningTrap(double B_0, std::function<double(double)> V_0,
+                         double d, double t)
 {
     this->B_0 = B_0;
     this->V_0 = V_0;
     this->d = d;
-
+    this->t = t;
 }
 
-PenningTrap::PenningTrap(int i, double B_0, double V_0, double d)
+PenningTrap::PenningTrap(unsigned int i, double B_0,
+                         std::function<double(double)> V_0, double d, double t)
     : PenningTrap::PenningTrap(B_0, V_0, d)
 {
     vec_3d r, v;
-    for (int j = 0; j < i; j++) {
+    for (size_t j = 0; j < i; j++) {
         r = vec_3d().randn() * .1 * this->d;
         v = vec_3d().randn() * .1 * this->d;
         this->add_particle(Particle(1., 40., r, v));
@@ -42,13 +37,13 @@ PenningTrap::PenningTrap(int i, double B_0, double V_0, double d)
 }
 
 PenningTrap::PenningTrap(std::vector<Particle> particles, double B_0,
-                         double V_0, double d)
+                         std::function<double(double)> V_0, double d, double t)
     : PenningTrap::PenningTrap(B_0, V_0, d)
 {
     this->particles = particles;
 }
 
-vec_3d PenningTrap::v_func(int i, int j, double dt)
+vec_3d PenningTrap::v_func(unsigned int i, unsigned int j, double dt)
 {
     switch (i) {
     case 0:
@@ -66,7 +61,7 @@ vec_3d PenningTrap::v_func(int i, int j, double dt)
     }
 }
 
-vec_3d PenningTrap::r_func(int i, int j, double dt)
+vec_3d PenningTrap::r_func(unsigned int i, unsigned int j, double dt)
 {
     switch (i) {
     case 0:
@@ -92,7 +87,7 @@ void PenningTrap::add_particle(Particle particle)
 vec_3d PenningTrap::external_E_field(vec_3d r)
 {
     r(2) *= -2.;
-    double f = this->V_0 / (this->d * this->d);
+    double f = this->V_0(this->t) / (this->d * this->d);
 
     return f * r;
 }
@@ -102,21 +97,21 @@ vec_3d PenningTrap::external_B_field(vec_3d r)
     return vec_3d{0., 0., this->B_0};
 }
 
-vec_3d PenningTrap::force_on_particle(int i, int j)
+vec_3d PenningTrap::force_on_particle(unsigned int i, unsigned int j)
 {
+    Particle p_j = this->particles[j];
     // Calculate the difference between the particles' position
-    vec_3d res = this->particles[i].r_vec - this->particles[j].r_vec;
+    vec_3d res = this->particles[i].r_vec - p_j.r_vec;
 
     // Get the distance between the particles
-    double norm = arma::norm(res);
+    double norm = arma::norm(res, 2);
 
     // Multiply res with p_j's charge divided by the norm cubed
-    res *= this->particles[j].q / (norm * norm * norm);
 
-    return res;
+    return vec_3d(res * p_j.q / (norm * norm * norm));
 }
 
-vec_3d PenningTrap::total_force_external(int i)
+vec_3d PenningTrap::total_force_external(unsigned int i)
 {
     Particle p = this->particles[i];
 
@@ -124,13 +119,9 @@ vec_3d PenningTrap::total_force_external(int i)
         return vec_3d{0., 0., 0.};
     }
 
-    vec_3d B = this->external_B_field(p.r_vec);
+    vec_3d force =
-
+        p.q * (this->external_E_field(p.r_vec) +
-    vec_3d v_cross_B{p.v_vec(1) * B(2) - p.v_vec(2) * B(1),
+               arma::cross(p.v_vec, this->external_B_field(p.r_vec)));
-                     p.v_vec(2) * B(0) - p.v_vec(0) * B(2),
-                     p.v_vec(0) * B(1) - p.v_vec(1) * B(0)};
-
-    vec_3d force = p.q * (this->external_E_field(p.r_vec) + v_cross_B);
 
     return force;
 }
@@ -140,9 +131,8 @@ vec_3d PenningTrap::total_force_particles(unsigned int i)
     Particle p = this->particles[i];
 
     vec_3d res;
-    size_t size = this->particles.size();
 
-    for (size_t j = 0; j < size; j++) {
+    for (size_t j = 0; j < this->particles.size(); j++) {
         if (i == j) {
             continue;
         }
@@ -150,12 +140,10 @@ vec_3d PenningTrap::total_force_particles(unsigned int i)
         res += this->force_on_particle(i, j);
     }
 
-    res *= K_E * (p.q / p.m);
+    return vec_3d(res * K_E * (p.q / p.m));
-
-    return res;
 }
 
-vec_3d PenningTrap::total_force(int i)
+vec_3d PenningTrap::total_force(unsigned int i)
 {
     return this->total_force_external(i) - this->total_force_particles(i);
 }
@@ -163,43 +151,36 @@ vec_3d PenningTrap::total_force(int i)
 void PenningTrap::evolve_RK4(double dt, bool particle_interaction)
 {
 
-    DEBUG("Inside RK4");
     std::vector<Particle> original_particles = this->particles;
     std::vector<Particle> tmp_particles = this->particles;
 
-    DEBUG("Choose force func");
+    vec_3d (PenningTrap::*force)(unsigned int);
-    std::function<vec_3d(int)> force;
     if (particle_interaction) {
-        force = [this](int i) { return this->total_force(i); };
+        force = &PenningTrap::total_force;
     }
     else {
-        force = [this](int i) { return this->total_force_external(i); };
+        force = &PenningTrap::total_force_external;
     }
 
-    DEBUG("Create k_v and k_r");
+    size_t size = this->particles.size();
-    int size = this->particles.size();
 
     this->k_v = sim_arr(4, sim_cols(size));
     this->k_r = sim_arr(4, sim_cols(size));
 
-    DEBUG("Loop");
+    for (size_t i = 0; i < 4; i++) {
-    for (int i = 0; i < 4; i++) {
-        DEBUG("Inside outer loop");
 #pragma omp parallel for
-        for (int j = 0; j < size; j++) {
+        for (size_t j = 0; j < this->particles.size(); j++) {
-            DEBUG("Set k_v and k_r");
+            this->k_v[i][j] = (this->*force)(j) / this->particles[j].m;
-            this->k_v[i][j] = this->total_force(j) / this->particles[j].m;
             this->k_r[i][j] = this->particles[j].v_vec;
 
             Particle *p = &tmp_particles[j];
 
-            DEBUG("Update v and r");
             p->v_vec = original_particles[j].v_vec + this->v_func(i, j, dt);
             p->r_vec = original_particles[j].r_vec + this->r_func(i, j, dt);
         }
-        DEBUG("After inner loop");
         this->particles = tmp_particles;
     }
+    this->t += dt;
 }
 
 void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction)
@@ -208,58 +189,49 @@ void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction)
 
     Particle *p;
 
-    std::function<vec_3d(int)> force;
+    vec_3d (PenningTrap::*force)(unsigned int);
     if (particle_interaction) {
-        force = [this](int i) { return this->total_force(i); };
+        force = &PenningTrap::total_force;
     }
     else {
-        force = [this](int i) { return this->total_force_external(i); };
+        force = &PenningTrap::total_force_external;
     }
 
 #pragma omp parallel for private(p)
     for (size_t i = 0; i < this->particles.size(); i++) {
         p = &new_state[i];
-        p->v_vec += dt * force(i) / p->m;
+        p->v_vec += dt * (this->*force)(i) / p->m;
         p->r_vec += dt * this->particles[i].v_vec;
     }
 
     this->particles = new_state;
+    this->t += dt;
 }
 
 sim_arr PenningTrap::simulate(double time, unsigned int steps,
                               std::string method, bool particle_interaction)
 {
-    DEBUG("Inside simulate");
     double dt = time / (double)steps;
 
     sim_arr res(this->particles.size(), sim_cols(steps));
 
-    DEBUG("Choose function to use");
+    void (PenningTrap::*func)(double, bool);
-    std::function<void(double, bool)> func;
     if (method == "rk4") {
-        func = [this](double dt, bool particle_interaction) {
+        func = &PenningTrap::evolve_RK4;
-            DEBUG("Inside RK4 lambda");
-            this->evolve_RK4(dt, particle_interaction);
-        };
     }
     else if (method == "euler") {
-        func = [this](double dt, bool particle_interaction) {
+        func = &PenningTrap::evolve_forward_euler;
-            this->evolve_forward_euler(dt, particle_interaction);
-        };
     }
     else {
         std::cout << "Not a valid method!" << std::endl;
         abort();
     }
 
-    DEBUG("Loop");
     for (size_t j = 0; j < steps; j++) {
-        DEBUG("Inside outer loop");
         for (size_t i = 0; i < this->particles.size(); i++) {
             res[i][j] = this->particles[i].r_vec;
         }
-        DEBUG("After inner loop");
+        (this->*func)(dt, particle_interaction);
-        func(dt, particle_interaction);
     }
 
     return res;
@@ -269,7 +241,6 @@ void PenningTrap::write_simulation_to_dir(std::string path, double time,
                                           int steps, std::string method,
                                           bool particle_interaction)
 {
-    DEBUG("Inside write to dir");
     if (path.back() != '/') {
         path += '/';
     }
@@ -278,12 +249,10 @@ void PenningTrap::write_simulation_to_dir(std::string path, double time,
         return;
     }
 
-    DEBUG("Create sim_arr");
     sim_arr res = this->simulate(time, steps, method, particle_interaction);
 
     std::ofstream ofile;
 
-    DEBUG("Loop");
 #pragma omp parallel for private(ofile)
     for (size_t i = 0; i < this->particles.size(); i++) {
         ofile.open(path + "particle_" + std::to_string(i) + ".txt");
@@ -293,3 +262,19 @@ void PenningTrap::write_simulation_to_dir(std::string path, double time,
         ofile.close();
     }
 }
+
+double PenningTrap::fraction_of_particles_left(double time, unsigned int steps, std::string method, bool particle_interaction)
+{
+    sim_arr res = this->simulate(time, steps, method, particle_interaction);
+
+    int particles_left = 0;
+
+    for (Particle p : this->particles) {
+        if (arma::norm(p.r_vec) < this->d) {
+            particles_left++;
+        }
+    }
+
+    return (double) particles_left / (double) this->particles.size();
+}
+

src/main.cpp

@@ -10,6 +10,7 @@
  *  @bug No known bugs
  * */
 
+#include <cmath>
 #include <fstream>
 #include <omp.h>
 #include <string>
@@ -59,45 +60,81 @@ void simulate_single_particle_with_different_steps()
     for (int i = 0; i < 4; i++) {
         int steps = 4000 * (i + 1);
         PenningTrap trap(std::vector<Particle>{p1});
-        trap.write_simulation_to_dir("output/N_steps/RK4/" + std::to_string(steps) +
+        trap.write_simulation_to_dir("output/N_steps/RK4/" +
-                                         "_steps",
+                                         std::to_string(steps) + "_steps",
                                      time, steps, "rk4", false);
     }
 
     for (int i = 0; i < 4; i++) {
         int steps = 4000 * (i + 1);
         PenningTrap trap(std::vector<Particle>{p1});
-        trap.write_simulation_to_dir("output/N_steps/euler/" + std::to_string(steps) +
+        trap.write_simulation_to_dir("output/N_steps/euler/" +
-                                         "_steps",
+                                         std::to_string(steps) + "_steps",
                                      time, steps, "euler", false);
     }
 }
 
 void simulate_100_particles()
 {
-    PenningTrap trap(100);
+    PenningTrap trap((unsigned)100, T,
+                     [](double t) {
+                         return 25. * V / 1000. *
+                                (1. + .4 * std::cos(1.5 * t));
+                     },
+                     500., 0);
 
-    double time = 50.; // microseconds
+    double time = 500.; // microseconds
 
-    trap.write_simulation_to_dir("output/simulate_100_particles", time, N);
+    trap.write_simulation_to_dir("output/simulate_100_particles", time, N*5);
+}
+
+void simulate_100_particles_with_time_potential()
+{
+    double amplitudes[]{.1, .4, .7};
+
+    double freq_start = .2;
+    double freq_end = 2.5;
+    double freq_increment = .02;
+    size_t freq_iterations = (size_t) ((freq_end - freq_start) / freq_increment);
+
+    std::string path = "output/time_dependent_potential/";
+    mkpath(path);
+
+    std::ofstream ofile;
+
+    for (double f : amplitudes) {
+        ofile.open(path + "f_" + std::to_string(f) + ".txt");
+        #pragma omp parallel for ordered schedule(static, 1)
+        for (size_t i=0; i < freq_iterations; i++) {
+            double freq = freq_start + i*freq_increment;
+            PenningTrap trap((unsigned)100, T,
+                             [f, freq](double t) {
+                                 return (25. * V / 1000.) *
+                                        (1. + f * std::cos(freq * t));
+                             },
+                             500., 0.);
+            double res = trap.fraction_of_particles_left(500., 40000, "rk4", true);
+            #pragma omp ordered
+            ofile << freq << "," << res << "\n";
+        }
+        ofile.close();
+    }
 }
 
 int main()
 {
-    DEBUG("Start");
-    double start = omp_get_wtime();
 
-    DEBUG("Single particle");
     simulate_single_particle();
 
-    DEBUG("Two particles");
     simulate_two_particles();
 
-    DEBUG("Single particle with different steps");
     simulate_single_particle_with_different_steps();
 
-    DEBUG("100 particles");
+    double start = omp_get_wtime();
-    simulate_100_particles();
+
+    //simulate_100_particles();
+
+    simulate_100_particles_with_time_potential();
 
     double end = omp_get_wtime();