FYS3150/Project-3
3
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Refactor

Browse Source 
- Change PenningTrap so that perturbation is its own field
- Add Calcium mass and charge to constants
This commit is contained in:
Cory Balaton 2023-10-23 09:45:08 +02:00
parent 2b69094ff4
commit 8fc0bf8c13
No known key found for this signature in database
GPG Key ID: 3E5FCEBFD80F432B
3 changed files with 56 additions and 50 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

40
include/PenningTrap.hpp
View File

@ -31,7 +31,8 @@
class PenningTrap { class PenningTrap {
private: private:
double B_0; ///< Magnetic field strength 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 d; ///< Characteristic dimension
double t; ///< Current time double t; ///< Current time
std::vector<Particle> particles; ///< The particles in the Penning trap std::vector<Particle> particles; ///< The particles in the Penning trap
@ -50,7 +51,7 @@ private:
* *
* @return vec_3d * @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 /** @brief Helper for evolve_RK4 when calculating \f$k_{r,i,j}\f$ values
* *
@ -62,7 +63,7 @@ private:
* *
* @return vec_3d * @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: public:
/** @brief Constructor for the PenningTrap class /** @brief Constructor for the PenningTrap class
@ -74,8 +75,7 @@ public:
* */ * */
PenningTrap( PenningTrap(
double B_0 = T, double B_0 = T,
std::function<double(double)> V_0 double V_0 = (25. * V) / 1000.,
= [](double t) { return (25. * V) / 1000.; },
double d = 500., double t = 0.); double d = 500., double t = 0.);
/** @brief Constructor for the PenningTrap class /** @brief Constructor for the PenningTrap class
@ -87,9 +87,8 @@ public:
* @param t The starting time * @param t The starting time
* */ * */
PenningTrap( PenningTrap(
unsigned int i, double B_0 = T, uint i, double B_0 = T,
std::function<double(double)> V_0 double V_0 = (25. * V) / 1000.,
= [](double t) { return 25. * V / 1000.; },
double d = 500., double t = 0.); double d = 500., double t = 0.);
/** @brief Constructor for the PenningTrap class /** @brief Constructor for the PenningTrap class
@ -102,10 +101,11 @@ public:
* */ * */
PenningTrap( PenningTrap(
std::vector<Particle> particles, double B_0 = T, std::vector<Particle> particles, double B_0 = T,
std::function<double(double)> V_0 double V_0 = (25. * V) / 1000.,
= [](double t) { return 25. * V / 1000.; },
double d = 500., double t = 0.); 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 /** @brief Give all particles new positions and velocities, and change t
* and V_0. * and V_0.
* *
@ -113,8 +113,7 @@ public:
* @param t The starting time * @param t The starting time
* */ * */
void reinitialize( void reinitialize(
std::function<double(double)> V_0 double f, double omega_V,
= [](double t) { return 25. * V / 1000.; },
double t = 0.); double t = 0.);
/** @brief Add a particle to the system /** @brief Add a particle to the system
@ -149,7 +148,7 @@ public:
* *
* @return vec_3d * @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. /** @brief Calculate the total external force on a particle.
* *
@ -160,7 +159,7 @@ public:
* *
* @return vec_3d * @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 /** @brief Calculate the total force on a particle p_i from other
* particles. * particles.
@ -169,7 +168,7 @@ public:
* *
* @return vec_3d * @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. /** @brief calculate the total force on a particle p_i.
* *
@ -177,7 +176,7 @@ public:
* *
* @return vec_3d * @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 /** @brief Go forward one timestep using the RK4 method
* *
@ -203,7 +202,7 @@ public:
* *
* @return simulation_t * @return simulation_t
* */ * */
simulation_t simulate(double time, unsigned int steps, simulation_t simulate(double time, uint steps,
std::string method = "rk4", std::string method = "rk4",
bool particle_interaction = true); bool particle_interaction = true);
@ -216,7 +215,7 @@ public:
* @param particle_interaction Turn particle interactions on/off * @param particle_interaction Turn particle interactions on/off
* */ * */
void write_simulation_to_dir(std::string path, double time, void write_simulation_to_dir(std::string path, double time,
unsigned int steps, uint steps,
std::string method = "rk4", std::string method = "rk4",
bool particle_interaction = true); bool particle_interaction = true);
@ -230,12 +229,9 @@ public:
* *
* @return double * @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", std::string method = "rk4",
bool particle_interaction = true); bool particle_interaction = true);
vec_3d get_r(int i);
double get_t();
}; };
#endif #endif

4
include/constants.hpp
View File

@ -18,4 +18,8 @@
#define V 9.64852558*1e7 ///< 1 Volt. unit: \f$ \frac{u (\mu m)^2}{(\mu s)^2 e} \f$ #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 #endif

60
src/PenningTrap.cpp
View File

@ -14,45 +14,52 @@
#include "typedefs.hpp" #include "typedefs.hpp"
#include <algorithm> #include <algorithm>
#include <functional> #include <functional>
#include <sys/types.h>
#include <vector> #include <vector>
PenningTrap::PenningTrap(double B_0, std::function<double(double)> V_0, PenningTrap::PenningTrap(double B_0, double V_0, double d, double t)
double d, double t)
{ {
this->B_0 = B_0; this->B_0 = B_0;
this->V_0 = V_0; this->V_0 = V_0;
this->d = d; this->d = d;
this->t = t; this->t = t;
this->perturbation = [](double t) { return 1.; };
} }
PenningTrap::PenningTrap(unsigned int i, double B_0, PenningTrap::PenningTrap(uint i, double B_0, double V_0, double d, double t)
std::function<double(double)> V_0, double d, double t)
: PenningTrap::PenningTrap(B_0, V_0, d) : PenningTrap::PenningTrap(B_0, V_0, d)
{ {
for (size_t j = 0; j < i; j++) { for (size_t j = 0; j < i; j++) {
this->particles.push_back(Particle(1., 40., 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))); vec_3d(vec_3d().randn() * .1 * this->d)));
} }
} }
PenningTrap::PenningTrap(std::vector<Particle> particles, double B_0, 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) : PenningTrap::PenningTrap(B_0, V_0, d)
{ {
this->particles = particles; 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->t = t;
this->set_pertubation(f, omega_V);
for (size_t i = 0; i < this->particles.size(); i++) { for (size_t i = 0; i < this->particles.size(); i++) {
this->particles[i].r_vec = vec_3d().randn() * .1 * this->d; 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) { switch (i) {
case 0: 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) { switch (i) {
case 0: case 0:
@ -99,7 +106,8 @@ vec_3d PenningTrap::external_E_field(vec_3d r)
{ {
r(2) *= -2.; 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) 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}; 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 // 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 - 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); 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]; 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)))); + 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; vec_3d res;
@ -141,10 +149,10 @@ vec_3d PenningTrap::total_force_particles(unsigned int i)
res += this->force_on_particle(i, j); 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) { if (arma::norm(this->particles[i].r_vec) > this->d) {
return vec_3d{0., 0., 0.}; 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> original_particles = this->particles;
std::vector<Particle> tmp_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 = particle_interaction ? &PenningTrap::total_force
: &PenningTrap::total_force_external; : &PenningTrap::total_force_external;
@ -197,7 +205,7 @@ void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction)
vec_3d force_res[size]; vec_3d force_res[size];
Particle *p; Particle *p;
vec_3d (PenningTrap::*force)(unsigned int) vec_3d (PenningTrap::*force)(uint)
= particle_interaction ? &PenningTrap::total_force = particle_interaction ? &PenningTrap::total_force
: &PenningTrap::total_force_external; : &PenningTrap::total_force_external;
@ -220,13 +228,12 @@ void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction)
this->t += dt; this->t += dt;
} }
simulation_t PenningTrap::simulate(double time, unsigned int steps, simulation_t PenningTrap::simulate(double time, uint steps, std::string method,
std::string method,
bool particle_interaction) bool particle_interaction)
{ {
double dt = time / (double)steps; 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)), simulation_t res{sim_arr(size, sim_cols(steps)),
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, void PenningTrap::write_simulation_to_dir(std::string path, double time,
unsigned int steps, uint steps, std::string method,
std::string method,
bool particle_interaction) bool particle_interaction)
{ {
if (path.back() != '/') { 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, std::string method,
bool particle_interaction) bool particle_interaction)
{ {
@ -317,9 +323,9 @@ double PenningTrap::fraction_of_particles_left(double time, unsigned int steps,
int particles_left = 0; int particles_left = 0;
// A reduction is perfect here // A reduction is perfect here
#pragma omp parallel for reduction(+:particles_left) #pragma omp parallel for reduction(+ : particles_left)
for (size_t i=0; i < this->particles.size(); i++) { for (size_t i = 0; i < this->particles.size(); i++) {
if (arma::norm(this->particles[i].r_vec) < this->d) { if (arma::norm(this->particles[i].r_vec) < this->d) {
particles_left++; particles_left++;
} }