Project-3/include/PenningTrap.hpp
2023-10-17 17:01:50 +02:00

228 lines
7.4 KiB
C++

/** @file PenningTrap.hpp
*
* @author Cory Alexander Balaton (coryab)
* @author Janita Ovidie Sandtrøen Willumsen (janitaws)
*
* @version 0.1
*
* @brief A class for simulating a Penning trap.
*
* @bug No known bugs
* */
#ifndef __PENNING_TRAP__
#define __PENNING_TRAP__
#include <armadillo>
#include <omp.h>
#include "Particle.hpp"
#include "constants.hpp"
#include "typedefs.hpp"
#pragma omp declare reduction(+ : vec_3d : omp_out += omp_in) \
initializer(omp_priv = omp_orig)
typedef struct simulation {
sim_arr r_vecs;
sim_arr v_vecs;
} simulation_t;
/** @brief A class that simulates a Penning trap.
*
* This class simulates a Penning trap. It can take in a number of particles
* and simulate how they would behave inside a Penning trap.
* */
class PenningTrap {
private:
double B_0; ///< Magnetic field strength
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; ///< A 2D vector containing all \f$k_{i,j}\f$ where \f$j\f$ is
///< 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
/** @brief Helper for evolve_RK4 when calculating \f$k_{v,i,j}\f$ values
*
* @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 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,
std::function<double(double)> V_0 =
[](double t) { return 25. * V / 1000.; },
double d = 500., double t = 0.);
/** @brief Constructor for the PenningTrap class
*
* @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.);
/** @brief Constructor for the PenningTrap class
*
* @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);
/** @brief Calculate the force between 2 particles.
*
* @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(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(unsigned int i);
/** @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 p_i.
*
* @param i The index of particle p_i
*
* @return vec_3d
* */
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
* */
simulation_t simulate(double time, unsigned int steps,
std::string method = "rk4",
bool particle_interaction = true);
/** @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);
};
#endif