diff --git a/include/PenningTrap.hpp b/include/PenningTrap.hpp
index 411b3cc..425afeb 100644
--- a/include/PenningTrap.hpp
+++ b/include/PenningTrap.hpp
@@ -31,7 +31,8 @@
 class PenningTrap {
 private:
     double B_0;                        ///< Magnetic field strength
-    std::function<double(double)> V_0; ///< Applied potential
+    double V_0; ///< Applied potential
+    std::function<double(double)> perturbation; ///< Time-dependent perturbation
     double d;                          ///< Characteristic dimension
     double t;                          ///< Current time
     std::vector<Particle> particles;   ///< The particles in the Penning trap
@@ -50,7 +51,7 @@ private:
      *
      *  @return vec_3d
      * */
-    vec_3d v_func(unsigned int i, unsigned int j, double dt);
+    vec_3d v_func(uint i, uint j, double dt);
 
     /** @brief Helper for evolve_RK4 when calculating \f$k_{r,i,j}\f$ values
      *
@@ -62,7 +63,7 @@ private:
      *
      *  @return vec_3d
      * */
-    vec_3d r_func(unsigned int i, unsigned int j, double dt);
+    vec_3d r_func(uint i, uint j, double dt);
 
 public:
     /** @brief Constructor for the PenningTrap class
@@ -74,8 +75,7 @@ public:
      * */
     PenningTrap(
         double B_0 = T,
-        std::function<double(double)> V_0
-        = [](double t) { return (25. * V) / 1000.; },
+        double V_0 = (25. * V) / 1000.,
         double d = 500., double t = 0.);
 
     /** @brief Constructor for the PenningTrap class
@@ -87,9 +87,8 @@ public:
      *  @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.; },
+        uint i, double B_0 = T,
+        double V_0 = (25. * V) / 1000.,
         double d = 500., double t = 0.);
 
     /** @brief Constructor for the PenningTrap class
@@ -102,10 +101,11 @@ public:
      * */
     PenningTrap(
         std::vector<Particle> particles, double B_0 = T,
-        std::function<double(double)> V_0
-        = [](double t) { return 25. * V / 1000.; },
+        double V_0 = (25. * V) / 1000.,
         double d = 500., double t = 0.);
 
+    void set_pertubation(double f, double omega_V);
+
     /** @brief Give all particles new positions and velocities, and change t 
      *  and V_0.
      *
@@ -113,8 +113,7 @@ public:
      *  @param t The starting time
      * */
     void reinitialize(
-        std::function<double(double)> V_0
-        = [](double t) { return 25. * V / 1000.; },
+        double f, double omega_V,
         double t = 0.);
 
     /** @brief Add a particle to the system
@@ -149,7 +148,7 @@ public:
      *
      *  @return vec_3d
      * */
-    vec_3d force_on_particle(unsigned int i, unsigned int j);
+    vec_3d force_on_particle(uint i, uint j);
 
     /** @brief Calculate the total external force on a particle.
      *
@@ -160,7 +159,7 @@ public:
      *
      *  @return vec_3d
      * */
-    vec_3d total_force_external(unsigned int i);
+    vec_3d total_force_external(uint i);
 
     /** @brief Calculate the total force on a particle p_i from other
      * particles.
@@ -169,7 +168,7 @@ public:
      *
      *  @return vec_3d
      * */
-    vec_3d total_force_particles(unsigned int i);
+    vec_3d total_force_particles(uint i);
 
     /** @brief calculate the total force on a particle p_i.
      *
@@ -177,7 +176,7 @@ public:
      *
      *  @return vec_3d
      * */
-    vec_3d total_force(unsigned int i);
+    vec_3d total_force(uint i);
 
     /** @brief Go forward one timestep using the RK4 method
      *
@@ -203,7 +202,7 @@ public:
      *
      *  @return simulation_t
      * */
-    simulation_t simulate(double time, unsigned int steps,
+    simulation_t simulate(double time, uint steps,
                           std::string method = "rk4",
                           bool particle_interaction = true);
 
@@ -216,7 +215,7 @@ public:
      *  @param particle_interaction Turn particle interactions on/off
      * */
     void write_simulation_to_dir(std::string path, double time,
-                                 unsigned int steps,
+                                 uint steps,
                                  std::string method = "rk4",
                                  bool particle_interaction = true);
 
@@ -230,12 +229,9 @@ public:
      *
      *  @return double
      * */
-    double fraction_of_particles_left(double time, unsigned int steps,
+    double fraction_of_particles_left(double time, uint steps,
                                       std::string method = "rk4",
                                       bool particle_interaction = true);
-
-    vec_3d get_r(int i);
-    double get_t();
 };
 
 #endif
diff --git a/include/constants.hpp b/include/constants.hpp
index 06942ee..bc8694a 100644
--- a/include/constants.hpp
+++ b/include/constants.hpp
@@ -18,4 +18,8 @@
 
 #define V 9.64852558*1e7 ///< 1 Volt. unit: \f$ \frac{u (\mu m)^2}{(\mu s)^2 e} \f$
 
+#define CA_MASS 40.078 ///< Mass of a single calcium ion. unit: amu
+
+#define CA_CHARGE 1. ///< Charge of a singly charged calcium ion. unit: e
+
 #endif
diff --git a/src/PenningTrap.cpp b/src/PenningTrap.cpp
index af4b19e..806901b 100644
--- a/src/PenningTrap.cpp
+++ b/src/PenningTrap.cpp
@@ -14,45 +14,52 @@
 #include "typedefs.hpp"
 #include <algorithm>
 #include <functional>
+#include <sys/types.h>
 #include <vector>
 
-PenningTrap::PenningTrap(double B_0, std::function<double(double)> V_0,
-                         double d, double t)
+PenningTrap::PenningTrap(double B_0, double V_0, double d, double t)
 {
     this->B_0 = B_0;
     this->V_0 = V_0;
     this->d = d;
     this->t = t;
+    this->perturbation = [](double t) { return 1.; };
 }
 
-PenningTrap::PenningTrap(unsigned int i, double B_0,
-                         std::function<double(double)> V_0, double d, double t)
+PenningTrap::PenningTrap(uint i, double B_0, double V_0, double d, double t)
     : PenningTrap::PenningTrap(B_0, V_0, d)
 {
     for (size_t j = 0; j < i; j++) {
-        this->particles.push_back(Particle(1., 40., vec_3d(vec_3d().randn() * .1 * this->d),
-                                    vec_3d(vec_3d().randn() * .1 * this->d)));
+        this->particles.push_back(
+            Particle(1., 40., vec_3d(vec_3d().randn() * .1 * this->d),
+                     vec_3d(vec_3d().randn() * .1 * this->d)));
     }
 }
 
 PenningTrap::PenningTrap(std::vector<Particle> particles, double B_0,
-                         std::function<double(double)> V_0, double d, double t)
+                         double V_0, double d, double t)
     : PenningTrap::PenningTrap(B_0, V_0, d)
 {
     this->particles = particles;
 }
 
-void PenningTrap::reinitialize(std::function<double(double)> V_0, double t)
+void PenningTrap::set_pertubation(double f, double omega_V)
+{
+    this->perturbation
+        = [f, omega_V](double t) { return f * std::cos(omega_V * t); };
+}
+
+void PenningTrap::reinitialize(double f, double omega_V, double t)
 {
-    this->V_0 = V_0;
     this->t = 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;
     }
 }
 
-vec_3d PenningTrap::v_func(unsigned int i, unsigned int j, double dt)
+vec_3d PenningTrap::v_func(uint i, uint j, double dt)
 {
     switch (i) {
     case 0:
@@ -71,7 +78,7 @@ vec_3d PenningTrap::v_func(unsigned int i, unsigned int j, double dt)
     }
 }
 
-vec_3d PenningTrap::r_func(unsigned int i, unsigned int j, double dt)
+vec_3d PenningTrap::r_func(uint i, uint j, double dt)
 {
     switch (i) {
     case 0:
@@ -99,7 +106,8 @@ vec_3d PenningTrap::external_E_field(vec_3d r)
 {
     r(2) *= -2.;
 
-    return vec_3d((this->V_0(this->t) / (this->d * this->d)) * r);
+    return vec_3d(
+        (this->V_0 * this->perturbation(this->t) / (this->d * this->d)) * r);
 }
 
 vec_3d PenningTrap::external_B_field(vec_3d r)
@@ -107,7 +115,7 @@ vec_3d PenningTrap::external_B_field(vec_3d r)
     return vec_3d{0., 0., this->B_0};
 }
 
-vec_3d PenningTrap::force_on_particle(unsigned int i, unsigned int j)
+vec_3d 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;
@@ -118,7 +126,7 @@ vec_3d PenningTrap::force_on_particle(unsigned int i, unsigned int j)
     return vec_3d((this->particles[j].q / (norm * norm * norm)) * res);
 }
 
-vec_3d PenningTrap::total_force_external(unsigned int i)
+vec_3d PenningTrap::total_force_external(uint i)
 {
     Particle *p = &this->particles[i];
 
@@ -132,7 +140,7 @@ vec_3d PenningTrap::total_force_external(unsigned int i)
            + arma::cross(p->v_vec, this->external_B_field(p->r_vec))));
 }
 
-vec_3d PenningTrap::total_force_particles(unsigned int i)
+vec_3d PenningTrap::total_force_particles(uint i)
 {
     vec_3d res;
 
@@ -141,10 +149,10 @@ vec_3d PenningTrap::total_force_particles(unsigned int i)
             res += this->force_on_particle(i, j);
     }
 
-    return vec_3d(res * K_E * (this->particles[i].q));
+    return vec_3d(res * (K_E * this->particles[i].q));
 }
 
-vec_3d PenningTrap::total_force(unsigned int i)
+vec_3d PenningTrap::total_force(uint i)
 {
     if (arma::norm(this->particles[i].r_vec) > this->d) {
         return vec_3d{0., 0., 0.};
@@ -159,7 +167,7 @@ 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)(unsigned int)
+    vec_3d (PenningTrap::*force)(uint)
         = particle_interaction ? &PenningTrap::total_force
                                : &PenningTrap::total_force_external;
 
@@ -197,7 +205,7 @@ void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction)
     vec_3d force_res[size];
     Particle *p;
 
-    vec_3d (PenningTrap::*force)(unsigned int)
+    vec_3d (PenningTrap::*force)(uint)
         = particle_interaction ? &PenningTrap::total_force
                                : &PenningTrap::total_force_external;
 
@@ -220,13 +228,12 @@ void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction)
     this->t += dt;
 }
 
-simulation_t PenningTrap::simulate(double time, unsigned int steps,
-                                   std::string method,
+simulation_t PenningTrap::simulate(double time, uint steps, std::string method,
                                    bool particle_interaction)
 {
     double dt = time / (double)steps;
 
-    unsigned int size = this->particles.size();
+    uint size = this->particles.size();
 
     simulation_t res{sim_arr(size, sim_cols(steps)),
                      sim_arr(size, sim_cols(steps))};
@@ -255,8 +262,7 @@ simulation_t PenningTrap::simulate(double time, unsigned int steps,
 }
 
 void PenningTrap::write_simulation_to_dir(std::string path, double time,
-                                          unsigned int steps,
-                                          std::string method,
+                                          uint steps, std::string method,
                                           bool particle_interaction)
 {
     if (path.back() != '/') {
@@ -293,7 +299,7 @@ void PenningTrap::write_simulation_to_dir(std::string path, double time,
     }
 }
 
-double PenningTrap::fraction_of_particles_left(double time, unsigned int steps,
+double PenningTrap::fraction_of_particles_left(double time, uint steps,
                                                std::string method,
                                                bool particle_interaction)
 {
@@ -317,9 +323,9 @@ double PenningTrap::fraction_of_particles_left(double time, unsigned int steps,
 
     int particles_left = 0;
 
-    // A reduction is perfect here
-    #pragma omp parallel for reduction(+:particles_left)
-    for (size_t i=0; i < this->particles.size(); i++) {
+// A reduction is perfect here
+#pragma omp parallel for reduction(+ : particles_left)
+    for (size_t i = 0; i < this->particles.size(); i++) {
         if (arma::norm(this->particles[i].r_vec) < this->d) {
             particles_left++;
         }