From 58c9212e5c8a48310d5f7419c59873f201a51bae Mon Sep 17 00:00:00 2001 From: Cory Date: Sat, 14 Oct 2023 03:12:23 +0200 Subject: [PATCH] Update documentation --- include/Particle.hpp | 5 ++ include/PenningTrap.hpp | 167 ++++++++++++++++++++++++++++++++++------ include/typedefs.hpp | 26 +++++++ include/utils.hpp | 59 ++++++++++---- src/Makefile | 2 +- src/Particle.cpp | 1 - src/PenningTrap.cpp | 125 +++++++++++++----------------- src/main.cpp | 65 ++++++++++++---- 8 files changed, 322 insertions(+), 128 deletions(-) diff --git a/include/Particle.hpp b/include/Particle.hpp index cf3232c..4b5e716 100644 --- a/include/Particle.hpp +++ b/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); diff --git a/include/PenningTrap.hpp b/include/PenningTrap.hpp index eb7bff7..f6cf2c0 100644 --- a/include/PenningTrap.hpp +++ b/include/PenningTrap.hpp @@ -19,8 +19,7 @@ #include "constants.hpp" #include "typedefs.hpp" - -#pragma omp declare reduction(+ : vec_3d : omp_out += omp_in) \ +#pragma omp declare reduction(+ : vec_3d : omp_out += omp_in) \ initializer(omp_priv = omp_orig) /** @brief A class that simulates a Penning trap. @@ -30,36 +29,101 @@ * */ class PenningTrap { private: - double B_0; ///< Magnetic field strength - double V_0; ///< Applied potential - double d; ///< Characteristic dimension - std::vector particles; ///< The particles in the Penning trap - sim_arr k_v; - sim_arr k_r; + double B_0; ///< Magnetic field strength + std::function V_0; ///< Applied potential + double d; ///< Characteristic dimension + double t; ///< Current time + std::vector 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 - vec_3d v_func(int, int, double); - vec_3d r_func(int, int, double); + /** @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 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 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., - double d = 500.); + /** @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 V_0 = + [](double t) { return 25. * V / 1000.; }, + double d = 500., double t = 0.); - PenningTrap(std::vector particles, double B_0 = T, - double V_0 = 25. * V / 1000., double d = 500.); + /** @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 particles, double B_0 = T, + std::function 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,39 +131,92 @@ 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); + std::string method = "rk4", + bool particle_interaction = true); - void write_simulation_to_dir(std::string path, double time, int steps, - std::string method = "rk4", + /** @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 diff --git a/include/typedefs.hpp b/include/typedefs.hpp index 42c00bb..f0ceca2 100644 --- a/include/typedefs.hpp +++ b/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 #include +/** @brief Typedef for the column of the result vector from simulating + * particles. + * */ typedef std::vector> sim_cols; + +/** @brief Typedef for the row of the result vector from simulating particles. + * */ typedef std::vector> sim_rows; + +/** @brief Typedef for the result of the simulate method. + * */ typedef std::vector sim_arr; + +/** @brief Typedef for a fixed 3d arma vector. + * */ typedef arma::vec::fixed<3> vec_3d; #endif diff --git a/include/utils.hpp b/include/utils.hpp index 2210880..c9c9485 100644 --- a/include/utils.hpp +++ b/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& 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 - * true, or it prints a fail message and aborts execution if it fails. + * @details This function takes in an expression and prints an OK message if + * 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 - * equal to each other given a tolerance. + * @details This function takes in 2 vectors and checks if they are + * 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 begin = prettyFunction.substr(0,colons).rfind(" ") + 1; - size_t end = prettyFunction.rfind("(") - begin; + size_t colons = pretty_function.find("::"); + size_t begin = pretty_function.substr(0,colons).rfind(" ") + 1; + 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 diff --git a/src/Makefile b/src/Makefile index 69f309e..f678e43 100644 --- a/src/Makefile +++ b/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) diff --git a/src/Particle.cpp b/src/Particle.cpp index 0ff18e8..2caad68 100644 --- a/src/Particle.cpp +++ b/src/Particle.cpp @@ -8,7 +8,6 @@ * @brief The implementation of the Particle class. * * @bug No known bugs - * * */ #include "Particle.hpp" diff --git a/src/PenningTrap.cpp b/src/PenningTrap.cpp index 0c50836..d0478ec 100644 --- a/src/PenningTrap.cpp +++ b/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 -#include -#include -#include -PenningTrap::PenningTrap(double B_0, double V_0, double d) +PenningTrap::PenningTrap(double B_0, std::function 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 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 particles, double B_0, - double V_0, double d) + std::function 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 v_cross_B{p.v_vec(1) * B(2) - p.v_vec(2) * B(1), - 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); + vec_3d force = + p.q * (this->external_E_field(p.r_vec) + + arma::cross(p.v_vec, this->external_B_field(p.r_vec))); 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 res; + return vec_3d(res * K_E * (p.q / p.m)); } -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 original_particles = this->particles; std::vector tmp_particles = this->particles; - DEBUG("Choose force func"); - std::function 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; } - DEBUG("Create k_v and k_r"); - int size = this->particles.size(); + size_t 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 (int i = 0; i < 4; i++) { - DEBUG("Inside outer loop"); + for (size_t i = 0; i < 4; i++) { #pragma omp parallel for - for (int j = 0; j < size; j++) { - DEBUG("Set k_v and k_r"); - this->k_v[i][j] = this->total_force(j) / this->particles[j].m; + for (size_t j = 0; j < this->particles.size(); j++) { + this->k_v[i][j] = (this->*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 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"); - std::function func; + void (PenningTrap::*func)(double, bool); if (method == "rk4") { - func = [this](double dt, bool particle_interaction) { - DEBUG("Inside RK4 lambda"); - this->evolve_RK4(dt, particle_interaction); - }; + func = &PenningTrap::evolve_RK4; } else if (method == "euler") { - func = [this](double dt, bool particle_interaction) { - this->evolve_forward_euler(dt, particle_interaction); - }; + func = &PenningTrap::evolve_forward_euler; } 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"); - func(dt, particle_interaction); + (this->*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(); +} + diff --git a/src/main.cpp b/src/main.cpp index 35ec896..b07ed93 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -10,6 +10,7 @@ * @bug No known bugs * */ +#include #include #include #include @@ -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{p1}); - trap.write_simulation_to_dir("output/N_steps/RK4/" + std::to_string(steps) + - "_steps", + trap.write_simulation_to_dir("output/N_steps/RK4/" + + 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{p1}); - trap.write_simulation_to_dir("output/N_steps/euler/" + std::to_string(steps) + - "_steps", + trap.write_simulation_to_dir("output/N_steps/euler/" + + 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"); - simulate_100_particles(); + double start = omp_get_wtime(); + + //simulate_100_particles(); + + simulate_100_particles_with_time_potential(); double end = omp_get_wtime();