diff --git a/.clang-format b/.clang-format new file mode 100644 index 0000000..d4d13ad --- /dev/null +++ b/.clang-format @@ -0,0 +1,10 @@ +UseTab: Never +IndentWidth: 4 +TabWidth: 4 +AccessModifierOffset: -4 +IndentAccessModifiers: false +AllowShortFunctionsOnASingleLine: false +BreakBeforeBraces: Custom +BraceWrapping: + AfterFunction: true + BeforeElse: true diff --git a/Doxyfile b/Doxyfile index 591e8ba..eeba27f 100644 --- a/Doxyfile +++ b/Doxyfile @@ -1070,7 +1070,7 @@ EXAMPLE_RECURSIVE = NO # that contain images that are to be included in the documentation (see the # \image command). -IMAGE_PATH = +IMAGE_PATH = ./images # The INPUT_FILTER tag can be used to specify a program that doxygen should # invoke to filter for each input file. Doxygen will invoke the filter program diff --git a/README.md b/README.md index 9d655f4..254c153 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,8 @@ # Penning trap simulation +![](./images/100_particles.gif) + + ## Credits The Doxygen theme used here is diff --git a/docs/100_particles.gif b/docs/100_particles.gif new file mode 100644 index 0000000..3dde053 Binary files /dev/null and b/docs/100_particles.gif differ diff --git a/docs/Particle_8hpp_source.html b/docs/Particle_8hpp_source.html index 76031eb..dbb28de 100644 --- a/docs/Particle_8hpp_source.html +++ b/docs/Particle_8hpp_source.html @@ -128,7 +128,7 @@ $(document).ready(function(){initNavTree('Particle_8hpp_source.html',''); initRe
Particle::q
double q
Charge.
Definition: Particle.hpp:21
Particle::r_vec
arma::vec::fixed< 3 > r_vec
position
Definition: Particle.hpp:23
Particle::m
double m
Mass.
Definition: Particle.hpp:22
-
PenningTrap
A class that simulates a Penning trap.
Definition: PenningTrap.hpp:25
+
PenningTrap
A class that simulates a Penning trap.
Definition: PenningTrap.hpp:29
diff --git a/docs/PenningTrap_8cpp.html b/docs/PenningTrap_8cpp.html index ceb79d5..62355e1 100644 --- a/docs/PenningTrap_8cpp.html +++ b/docs/PenningTrap_8cpp.html @@ -104,12 +104,9 @@ $(document).ready(function(){initNavTree('PenningTrap_8cpp.html',''); initResiza

The implementation of the PenningTrap class. More...

-
#include "utils.hpp"
-#include "PenningTrap.hpp"
+
#include "PenningTrap.hpp"
#include "constants.hpp"
-#include <algorithm>
-#include <stdexcept>
-#include <omp.h>
+#include "utils.hpp"

Go to the source code of this file.

Detailed Description

diff --git a/docs/PenningTrap_8cpp_source.html b/docs/PenningTrap_8cpp_source.html index 0be7c11..993de8e 100644 --- a/docs/PenningTrap_8cpp_source.html +++ b/docs/PenningTrap_8cpp_source.html @@ -102,156 +102,202 @@ $(document).ready(function(){initNavTree('PenningTrap_8cpp_source.html',''); ini
Go to the documentation of this file.
1
-
16#include "utils.hpp"
-
17#include "PenningTrap.hpp"
-
18#include "constants.hpp"
-
19#include <algorithm>
-
20#include <stdexcept>
-
21#include <omp.h>
-
22
-
23#pragma omp declare reduction( + : arma::vec : omp_out += omp_in ) \
-
24 initializer( omp_priv = omp_orig )
-
25
-
26PenningTrap::PenningTrap(double B_0, double V_0, double d)
-
27{
-
28 this->B_0 = B_0;
-
29 this->V_0 = V_0;
-
30 this->d = d;
-
31}
-
32
-
33void PenningTrap::add_particle(Particle particle)
-
34{
-
35 this->particles.push_back(particle);
-
36}
-
37
-
38arma::vec PenningTrap::external_E_field(arma::vec r)
-
39{
-
40 arma::vec::fixed<3> res;
-
41 double f = this->V_0/(this->d*this->d);
-
42 res(0) = r(0);
-
43 res(1) = r(1);
-
44 res(2) = -2.*r(2);
-
45
-
46 return f*res;
-
47}
-
48
-
49arma::vec PenningTrap::external_B_field(arma::vec r)
-
50{
-
51 arma::vec::fixed<3> res;
-
52 res(0) = 0.;
-
53 res(1) = 0.;
-
54 res(2) = this->B_0;
+
16#include "PenningTrap.hpp"
+
17#include "constants.hpp"
+
18#include "utils.hpp"
+
19
+
20PenningTrap::PenningTrap(double B_0, double V_0, double d)
+
21{
+
22 this->B_0 = B_0;
+
23 this->V_0 = V_0;
+
24 this->d = d;
+
25}
+
26
+
27void PenningTrap::add_particle(Particle particle)
+
28{
+
29 this->particles.push_back(particle);
+
30}
+
31
+
32arma::vec PenningTrap::external_E_field(arma::vec r)
+
33{
+
34 arma::vec::fixed<3> res;
+
35 double f = this->V_0 / (this->d * this->d);
+
36 res(0) = r(0);
+
37 res(1) = r(1);
+
38 res(2) = -2. * r(2);
+
39
+
40 return f * res;
+
41}
+
42
+
43arma::vec PenningTrap::external_B_field(arma::vec r)
+
44{
+
45 arma::vec::fixed<3> res{0., 0., this->B_0};
+
46
+
47 return res;
+
48}
+
49
+
50arma::vec PenningTrap::force_on_particle(int i, int j)
+
51{
+
52 // Calculate the difference between the particles' position
+
53 arma::vec::fixed<3> res =
+
54 this->particles.at(i).r_vec - this->particles.at(j).r_vec;
55
-
56 return res;
-
57}
+
56 // Get the distance between the particles
+
57 double norm = arma::norm(res);
58
-
59arma::vec PenningTrap::force_on_particle(int i, int j)
-
60{
-
61 // Calculate the difference between the particles' position
-
62 arma::vec::fixed<3> res = this->particles.at(i).r_vec
-
63 - this->particles.at(j).r_vec;
+
59 // Multiply res with p_j's charge divided by the norm cubed
+
60 res *= this->particles.at(j).q / (norm * norm * norm);
+
61
+
62 return res;
+
63}
64
-
65 // Get the distance between the particles
-
66 double norm = arma::norm(res);
-
67
-
68 // Multiply res with p_j's charge divided by the norm cubed
-
69 res *= this->particles.at(j).q/(norm*norm*norm);
+
65arma::vec PenningTrap::total_force_external(int i)
+
66{
+
67 Particle p = this->particles.at(i);
+
68
+
69 arma::vec::fixed<3> B = this->external_B_field(p.r_vec);
70
-
71 return res;
-
72}
-
73
-
74arma::vec PenningTrap::total_force_external(int i)
-
75{
-
76 Particle p = this->particles.at(i);
-
77
-
78 arma::vec::fixed<3> v_cross_B;
+
71 arma::vec::fixed<3> v_cross_B{p.v_vec(1) * B(2) - p.v_vec(2) * B(1),
+
72 p.v_vec(2) * B(0) - p.v_vec(0) * B(2),
+
73 p.v_vec(0) * B(1) - p.v_vec(1) * B(0)};
+
74
+
75 arma::vec force = p.q * (this->external_E_field(p.r_vec) + v_cross_B);
+
76
+
77 return force;
+
78}
79
-
80 arma::vec::fixed<3> B = this->external_B_field(p.r_vec);
-
81
-
82 v_cross_B(0) = p.v_vec(1)*B(2) - p.v_vec(2)*B(1);
-
83 v_cross_B(1) = p.v_vec(2)*B(0) - p.v_vec(0)*B(2);
-
84 v_cross_B(2) = p.v_vec(0)*B(1) - p.v_vec(1)*B(0);
+
80arma::vec PenningTrap::total_force_particles(int i)
+
81{
+
82 Particle p = this->particles.at(i);
+
83
+
84 arma::vec res(3);
85
-
86 arma::vec force = p.q
-
87 *(this->external_E_field(p.r_vec) + v_cross_B);
-
88
-
89 return force;
-
90}
-
91
-
92arma::vec PenningTrap::total_force_particles(int i)
-
93{
-
94 Particle p = this->particles.at(i);
+
86 for (int j = 0; j < this->particles.size(); j++) {
+
87 if (i == j) {
+
88 continue;
+
89 }
+
90
+
91 res += this->force_on_particle(i, j);
+
92 }
+
93
+
94 res *= K_E * (p.q / p.m);
95
-
96 arma::vec res(3);
-
97
-
98 for (int j=0; j < this->particles.size(); j++) {
-
99 if (i == j) {
-
100 continue;
-
101 }
-
102
-
103 res += this->force_on_particle(i, j);
-
104 }
-
105
-
106 res *= K_E*(p.q/p.m);
-
107
-
108 return res;
-
109}
+
96 return res;
+
97}
+
98
+
99arma::vec PenningTrap::total_force(int i)
+
100{
+
101 return this->total_force_external(i) - this->total_force_particles(i);
+
102}
+
103
+
104void PenningTrap::evolve_RK4(double dt)
+
105{
+
106 std::vector<Particle> tmp_particles = this->particles;
+
107
+
108 arma::vec::fixed<3> *k_v = new arma::vec::fixed<3>[this->particles.size()*4];
+
109 arma::vec::fixed<3> *k_r = new arma::vec::fixed<3>[this->particles.size()*4];
110
-
111arma::vec PenningTrap::total_force(int i)
-
112{
-
113 return this->total_force_external(i) - this->total_force_particles(i);
-
114}
-
115
-
116void PenningTrap::evolve_RK4(double dt)
-
117{
-
118
-
119}
+
111 int size = this->particles.size();
+
112
+
113 for (int i=0; i<size; i++) {
+
114 k_v[i] = this->total_force(i)/this->particles.at(i).m;
+
115 k_r[i] = this->particles.at(i).v_vec;
+
116 }
+
117
+
118 for (int i=0; i<size; i++) {
+
119 Particle *p = &this->particles.at(i);
120
-
121void PenningTrap::evolve_forward_euler(double dt)
-
122{
-
123 std::vector<Particle> new_state = this->particles;
+
121 p->v_vec = tmp_particles.at(i).v_vec + (dt/2)*k_v[i];
+
122 p->r_vec = tmp_particles.at(i).r_vec + (dt/2)*k_r[i];
+
123 }
124
-
125 Particle *p;
-
126
-
127 #pragma omp parallel for private(p)
-
128 for (int i=0; i < this->particles.size(); i++) {
-
129 p = &new_state.at(i);
-
130 p->v_vec += dt*this->total_force(i)/new_state.at(i).m;
-
131 p->r_vec += dt*this->particles.at(i).v_vec;
-
132 }
+
125
+
126 for (int i=0; i<size; i++) {
+
127 k_v[1*size + i] = this->total_force(i)/this->particles.at(i).m;
+
128 k_r[1*size + i] = this->particles.at(i).v_vec;
+
129 }
+
130
+
131 for (int i=0; i<size; i++) {
+
132 Particle *p = &this->particles.at(i);
133
-
134 this->particles = new_state;
-
135}
-
136
-
137arma::vec PenningTrap::get_particle(int i)
-
138{
-
139 return this->particles.at(i).r_vec;
-
140}
-
141
-
142double PenningTrap::get_d()
-
143{
-
144 return this->d;
-
145}
+
134 p->v_vec = tmp_particles.at(i).v_vec + (dt/2)*k_v[1*size + i];
+
135 p->r_vec = tmp_particles.at(i).r_vec + (dt/2)*k_r[1*size + i];
+
136 }
+
137
+
138 for (int i=0; i<size; i++) {
+
139 k_v[2*size + i] = this->total_force(i)/this->particles.at(i).m;
+
140 k_r[2*size + i] = this->particles.at(i).v_vec;
+
141 }
+
142
+
143 for (int i=0; i<size; i++) {
+
144 Particle *p = &this->particles.at(i);
+
145
+
146 p->v_vec = tmp_particles.at(i).v_vec + dt*k_v[2*size + i];
+
147 p->r_vec = tmp_particles.at(i).r_vec + dt*k_r[2*size + i];
+
148 }
+
149
+
150
+
151 for (int i=0; i<size; i++) {
+
152 k_v[3*size + i] = this->total_force(i)/this->particles.at(i).m;
+
153 k_r[3*size + i] = this->particles.at(i).v_vec;
+
154 }
+
155
+
156 for (int i=0; i<size; i++) {
+
157 Particle *p = &this->particles.at(i);
+
158
+
159 p->v_vec = tmp_particles.at(i).v_vec + dt*(k_v[i] + k_v[size + i] + k_v[2*size + i] + k_v[3*size + i])/6;
+
160 p->r_vec = tmp_particles.at(i).r_vec + dt*(k_r[i] + k_r[size + i] + k_r[2*size + i] + k_r[3*size + i])/6;
+
161 }
+
162
+
163 delete [] k_v;
+
164 delete [] k_r;
+
165}
+
166
+
167void PenningTrap::evolve_forward_euler(double dt)
+
168{
+
169 std::vector<Particle> new_state = this->particles;
+
170
+
171 Particle *p;
+
172
+
173#pragma omp parallel for private(p)
+
174 for (int i = 0; i < this->particles.size(); i++) {
+
175 p = &new_state.at(i);
+
176 p->v_vec += dt * this->total_force(i) / new_state.at(i).m;
+
177 p->r_vec += dt * this->particles.at(i).v_vec;
+
178 }
+
179
+
180 this->particles = new_state;
+
181}
+
182
+
183arma::vec PenningTrap::get_particle(int i)
+
184{
+
185 return this->particles.at(i).r_vec;
+
186}
+
187
+
188double PenningTrap::get_d()
+
189{
+
190 return this->d;
+
191}
PenningTrap.hpp
A class for simulating a Penning trap.
Particle
A class that holds attributes of a particle.
Definition: Particle.hpp:19
Particle::v_vec
arma::vec::fixed< 3 > v_vec
velocity
Definition: Particle.hpp:24
Particle::q
double q
Charge.
Definition: Particle.hpp:21
Particle::r_vec
arma::vec::fixed< 3 > r_vec
position
Definition: Particle.hpp:23
Particle::m
double m
Mass.
Definition: Particle.hpp:22
-
PenningTrap::particles
std::vector< Particle > particles
The particles in the Penning trap.
Definition: PenningTrap.hpp:30
-
PenningTrap::total_force_external
arma::vec total_force_external(int i)
Calculate the total external force on a particle.
Definition: PenningTrap.cpp:74
-
PenningTrap::B_0
double B_0
Magnetic field strength.
Definition: PenningTrap.hpp:27
-
PenningTrap::total_force_particles
arma::vec total_force_particles(int i)
Calculate the total force on a particle from other particles.
Definition: PenningTrap.cpp:92
-
PenningTrap::external_B_field
arma::vec external_B_field(arma::vec r)
Calculate B at point r.
Definition: PenningTrap.cpp:49
-
PenningTrap::force_on_particle
arma::vec force_on_particle(int i, int j)
Calculate the force between 2 particles.
Definition: PenningTrap.cpp:59
-
PenningTrap::evolve_forward_euler
void evolve_forward_euler(double dt)
Go forward one timestep using the forward Euler method.
Definition: PenningTrap.cpp:121
-
PenningTrap::d
double d
Characteristic dimension.
Definition: PenningTrap.hpp:29
-
PenningTrap::add_particle
void add_particle(Particle particle)
Add a particle to the system.
Definition: PenningTrap.cpp:33
-
PenningTrap::V_0
double V_0
Applied potential.
Definition: PenningTrap.hpp:28
-
PenningTrap::PenningTrap
PenningTrap(double B_0=T, double V_0=25.*V/1000., double d=500.)
Set B_0, V_0 and d.
Definition: PenningTrap.cpp:26
-
PenningTrap::total_force
arma::vec total_force(int i)
calculate the total force on a particle.
Definition: PenningTrap.cpp:111
-
PenningTrap::external_E_field
arma::vec external_E_field(arma::vec r)
Calculate E at point r.
Definition: PenningTrap.cpp:38
-
PenningTrap::evolve_RK4
void evolve_RK4(double dt)
Go forward one timestep using the RK4 method.
Definition: PenningTrap.cpp:116
+
PenningTrap::particles
std::vector< Particle > particles
The particles in the Penning trap.
Definition: PenningTrap.hpp:34
+
PenningTrap::total_force_external
arma::vec total_force_external(int i)
Calculate the total external force on a particle.
Definition: PenningTrap.cpp:65
+
PenningTrap::B_0
double B_0
Magnetic field strength.
Definition: PenningTrap.hpp:31
+
PenningTrap::total_force_particles
arma::vec total_force_particles(int i)
Calculate the total force on a particle from other particles.
Definition: PenningTrap.cpp:80
+
PenningTrap::external_B_field
arma::vec external_B_field(arma::vec r)
Calculate B at point r.
Definition: PenningTrap.cpp:43
+
PenningTrap::force_on_particle
arma::vec force_on_particle(int i, int j)
Calculate the force between 2 particles.
Definition: PenningTrap.cpp:50
+
PenningTrap::evolve_forward_euler
void evolve_forward_euler(double dt)
Go forward one timestep using the forward Euler method.
Definition: PenningTrap.cpp:167
+
PenningTrap::d
double d
Characteristic dimension.
Definition: PenningTrap.hpp:33
+
PenningTrap::add_particle
void add_particle(Particle particle)
Add a particle to the system.
Definition: PenningTrap.cpp:27
+
PenningTrap::V_0
double V_0
Applied potential.
Definition: PenningTrap.hpp:32
+
PenningTrap::PenningTrap
PenningTrap(double B_0=T, double V_0=25.*V/1000., double d=500.)
Set B_0, V_0 and d.
Definition: PenningTrap.cpp:20
+
PenningTrap::total_force
arma::vec total_force(int i)
calculate the total force on a particle.
Definition: PenningTrap.cpp:99
+
PenningTrap::external_E_field
arma::vec external_E_field(arma::vec r)
Calculate E at point r.
Definition: PenningTrap.cpp:32
+
PenningTrap::evolve_RK4
void evolve_RK4(double dt)
Go forward one timestep using the RK4 method.
Definition: PenningTrap.cpp:104
constants.hpp
Library of constants.
K_E
#define K_E
Coulomb constant. unit: .
Definition: constants.hpp:15
utils.hpp
Function prototypes and macros that are useful.
diff --git a/docs/PenningTrap_8hpp.html b/docs/PenningTrap_8hpp.html index 8e1232f..306bc89 100644 --- a/docs/PenningTrap_8hpp.html +++ b/docs/PenningTrap_8hpp.html @@ -107,8 +107,9 @@ $(document).ready(function(){initNavTree('PenningTrap_8hpp.html',''); initResiza

A class for simulating a Penning trap. More...

#include <armadillo>
-#include "constants.hpp"
+#include <omp.h>
#include "Particle.hpp"
+#include "constants.hpp"

Go to the source code of this file.

diff --git a/docs/PenningTrap_8hpp_source.html b/docs/PenningTrap_8hpp_source.html index b8a0d43..02b2aed 100644 --- a/docs/PenningTrap_8hpp_source.html +++ b/docs/PenningTrap_8hpp_source.html @@ -106,60 +106,64 @@ $(document).ready(function(){initNavTree('PenningTrap_8hpp_source.html',''); ini
13#define __PENNING_TRAP__
14
15#include <armadillo>
-
16
-
17#include "constants.hpp"
+
16#include <omp.h>
+
17
18#include "Particle.hpp"
-
19
-
25class PenningTrap {
-
26private:
-
27 double B_0;
-
28 double V_0;
-
29 double d;
-
30 std::vector<Particle> particles;
-
31
-
32public:
-
35 PenningTrap(double B_0 = T, double V_0 = 25.*V/1000., double d = 500.);
-
36
-
39 void add_particle(Particle particle);
+
19#include "constants.hpp"
+
20
+
21#pragma omp declare reduction( + : arma::vec : omp_out += omp_in ) \
+
22 initializer( omp_priv = omp_orig )
+
23
+
29class PenningTrap {
+
30private:
+
31 double B_0;
+
32 double V_0;
+
33 double d;
+
34 std::vector<Particle> particles;
+
35
+
36public:
+
39 PenningTrap(double B_0 = T, double V_0 = 25.*V/1000., double d = 500.);
40
-
43 arma::vec external_E_field(arma::vec r);
+
43 void add_particle(Particle particle);
44
-
47 arma::vec external_B_field(arma::vec r);
+
47 arma::vec external_E_field(arma::vec r);
48
-
54 arma::vec force_on_particle(int i, int j);
-
55
-
61 arma::vec total_force_external(int i);
-
62
-
65 arma::vec total_force_particles(int i);
+
51 arma::vec external_B_field(arma::vec r);
+
52
+
58 arma::vec force_on_particle(int i, int j);
+
59
+
65 arma::vec total_force_external(int i);
66
-
69 arma::vec total_force(int i);
+
69 arma::vec total_force_particles(int i);
70
-
73 void evolve_RK4(double dt);
+
73 arma::vec total_force(int i);
74
-
77 void evolve_forward_euler(double dt);
+
77 void evolve_RK4(double dt);
78
-
79 arma::vec get_particle(int i);
-
80
-
81 double get_d();
-
82};
-
83
-
84#endif
+
81 void evolve_forward_euler(double dt);
+
82
+
83 arma::vec get_particle(int i);
+
84
+
85 double get_d();
+
86};
+
87
+
88#endif
Particle.hpp
A class that holds the properties of a particle.
Particle
A class that holds attributes of a particle.
Definition: Particle.hpp:19
-
PenningTrap
A class that simulates a Penning trap.
Definition: PenningTrap.hpp:25
-
PenningTrap::particles
std::vector< Particle > particles
The particles in the Penning trap.
Definition: PenningTrap.hpp:30
-
PenningTrap::total_force_external
arma::vec total_force_external(int i)
Calculate the total external force on a particle.
Definition: PenningTrap.cpp:74
-
PenningTrap::B_0
double B_0
Magnetic field strength.
Definition: PenningTrap.hpp:27
-
PenningTrap::total_force_particles
arma::vec total_force_particles(int i)
Calculate the total force on a particle from other particles.
Definition: PenningTrap.cpp:92
-
PenningTrap::external_B_field
arma::vec external_B_field(arma::vec r)
Calculate B at point r.
Definition: PenningTrap.cpp:49
-
PenningTrap::force_on_particle
arma::vec force_on_particle(int i, int j)
Calculate the force between 2 particles.
Definition: PenningTrap.cpp:59
-
PenningTrap::evolve_forward_euler
void evolve_forward_euler(double dt)
Go forward one timestep using the forward Euler method.
Definition: PenningTrap.cpp:121
-
PenningTrap::d
double d
Characteristic dimension.
Definition: PenningTrap.hpp:29
-
PenningTrap::add_particle
void add_particle(Particle particle)
Add a particle to the system.
Definition: PenningTrap.cpp:33
-
PenningTrap::V_0
double V_0
Applied potential.
Definition: PenningTrap.hpp:28
-
PenningTrap::total_force
arma::vec total_force(int i)
calculate the total force on a particle.
Definition: PenningTrap.cpp:111
-
PenningTrap::external_E_field
arma::vec external_E_field(arma::vec r)
Calculate E at point r.
Definition: PenningTrap.cpp:38
-
PenningTrap::evolve_RK4
void evolve_RK4(double dt)
Go forward one timestep using the RK4 method.
Definition: PenningTrap.cpp:116
+
PenningTrap
A class that simulates a Penning trap.
Definition: PenningTrap.hpp:29
+
PenningTrap::particles
std::vector< Particle > particles
The particles in the Penning trap.
Definition: PenningTrap.hpp:34
+
PenningTrap::total_force_external
arma::vec total_force_external(int i)
Calculate the total external force on a particle.
Definition: PenningTrap.cpp:65
+
PenningTrap::B_0
double B_0
Magnetic field strength.
Definition: PenningTrap.hpp:31
+
PenningTrap::total_force_particles
arma::vec total_force_particles(int i)
Calculate the total force on a particle from other particles.
Definition: PenningTrap.cpp:80
+
PenningTrap::external_B_field
arma::vec external_B_field(arma::vec r)
Calculate B at point r.
Definition: PenningTrap.cpp:43
+
PenningTrap::force_on_particle
arma::vec force_on_particle(int i, int j)
Calculate the force between 2 particles.
Definition: PenningTrap.cpp:50
+
PenningTrap::evolve_forward_euler
void evolve_forward_euler(double dt)
Go forward one timestep using the forward Euler method.
Definition: PenningTrap.cpp:167
+
PenningTrap::d
double d
Characteristic dimension.
Definition: PenningTrap.hpp:33
+
PenningTrap::add_particle
void add_particle(Particle particle)
Add a particle to the system.
Definition: PenningTrap.cpp:27
+
PenningTrap::V_0
double V_0
Applied potential.
Definition: PenningTrap.hpp:32
+
PenningTrap::total_force
arma::vec total_force(int i)
calculate the total force on a particle.
Definition: PenningTrap.cpp:99
+
PenningTrap::external_E_field
arma::vec external_E_field(arma::vec r)
Calculate E at point r.
Definition: PenningTrap.cpp:32
+
PenningTrap::evolve_RK4
void evolve_RK4(double dt)
Go forward one timestep using the RK4 method.
Definition: PenningTrap.cpp:104
constants.hpp
Library of constants.
T
#define T
1 Tesla. unit:
Definition: constants.hpp:17
V
#define V
1 Volt. unit:
Definition: constants.hpp:19
diff --git a/docs/animate__100__particles_8py_source.html b/docs/animate__100__particles_8py_source.html new file mode 100644 index 0000000..80feb43 --- /dev/null +++ b/docs/animate__100__particles_8py_source.html @@ -0,0 +1,178 @@ + + + + + + + +Penning Trap Simulation: src/animate_100_particles.py Source File + + + + + + + + + + + + + + + + +
+
+
+ + + + + +
+
Penning Trap Simulation +
+
Simulate particle behavior inside a Penning Trap
+
+
+ + + + + + + +
+
+ +
+
+
+ +
+ +
+
+ + +
+
+
+
+
+
Loading...
+
Searching...
+
No Matches
+
+
+
+
+ +
+
animate_100_particles.py
+
+
+
1import matplotlib.pyplot as plt
+
2import numpy as np
+
3from mpl_toolkits.mplot3d import Axes3D
+
4from matplotlib import animation
+
5
+
6def get_data(files):
+
7 res = []
+
8 for file in files:
+
9 arr = [[], [], []]
+
10 with open(file, encoding="utf8") as f:
+
11 lines = f.readlines()
+
12
+
13 for line in lines:
+
14 xi,yi,zi = map(float, line.strip().split(","))
+
15 arr[0].append(xi)
+
16 arr[1].append(yi)
+
17 arr[2].append(zi)
+
18 res.append(arr)
+
19
+
20 return np.array(res)
+
21
+
22def update(num, lines, arr):
+
23 for line, a in zip(lines, arr):
+
24 line.set_data(a[:2, num])
+
25 line.set_3d_properties(a[2, num])
+
26
+
27
+
28
+
29def animate():
+
30 plt.style.use("dark_background")
+
31 fig = plt.figure()
+
32 ax = fig.add_subplot(projection="3d")
+
33
+
34
+
35 arr = get_data([f"output/p{i}_RK4.txt" for i in range(100)])
+
36
+
37 arr = arr[:,:,::10]
+
38
+
39 N = len(arr[0][0])
+
40
+
41 lines = [ax.plot(*a[:,1], "o")[0] for a in arr]
+
42
+
43 ax.set_xlim3d([-500.0, 500.0])
+
44 ax.set_xlabel('X')
+
45
+
46 ax.set_ylim3d([-500.0, 500.0])
+
47 ax.set_ylabel('Y')
+
48
+
49 ax.set_zlim3d([-500.0, 500.0])
+
50 ax.set_zlabel('Z')
+
51
+
52 ani = animation.FuncAnimation(fig, update, N, fargs=(lines, arr),
+
53 interval=1,
+
54 blit=False)
+
55
+
56
+
57 ani.save("100_particles.gif", writer=animation.FFMpegFileWriter(fps=30))
+
58 plt.show()
+
59
+
60
+
61if __name__ == "__main__":
+
62 animate()
+
63
+
64
+
+
+ + + + diff --git a/docs/annotated.html b/docs/annotated.html index 7f9374f..fdee9b2 100644 --- a/docs/annotated.html +++ b/docs/annotated.html @@ -105,6 +105,7 @@ $(document).ready(function(){initNavTree('annotated.html',''); initResizable(); +
 CParticleA class that holds attributes of a particle
 CPenningTrapA class that simulates a Penning trap
 CPenningTrapTest
diff --git a/docs/annotated_dup.js b/docs/annotated_dup.js index 1f9034f..574cc6b 100644 --- a/docs/annotated_dup.js +++ b/docs/annotated_dup.js @@ -1,5 +1,6 @@ var annotated_dup = [ [ "Particle", "classParticle.html", "classParticle" ], - [ "PenningTrap", "classPenningTrap.html", "classPenningTrap" ] + [ "PenningTrap", "classPenningTrap.html", "classPenningTrap" ], + [ "PenningTrapTest", "classPenningTrapTest.html", null ] ]; \ No newline at end of file diff --git a/docs/classPenningTrap.html b/docs/classPenningTrap.html index 7809dc2..af62e25 100644 --- a/docs/classPenningTrap.html +++ b/docs/classPenningTrap.html @@ -167,7 +167,7 @@ Private Attributes

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.

-

Definition at line 25 of file PenningTrap.hpp.

+

Definition at line 29 of file PenningTrap.hpp.

Constructor & Destructor Documentation

◆ PenningTrap()

@@ -203,7 +203,7 @@ Private Attributes

Set B_0, V_0 and d.

-

Definition at line 26 of file PenningTrap.cpp.

+

Definition at line 20 of file PenningTrap.cpp.

@@ -226,7 +226,7 @@ Private Attributes

Add a particle to the system.

-

Definition at line 33 of file PenningTrap.cpp.

+

Definition at line 27 of file PenningTrap.cpp.

@@ -248,7 +248,7 @@ Private Attributes

Go forward one timestep using the forward Euler method.

-

Definition at line 121 of file PenningTrap.cpp.

+

Definition at line 167 of file PenningTrap.cpp.

@@ -270,7 +270,7 @@ Private Attributes

Go forward one timestep using the RK4 method.

-

Definition at line 116 of file PenningTrap.cpp.

+

Definition at line 104 of file PenningTrap.cpp.

@@ -292,7 +292,7 @@ Private Attributes

Calculate B at point r.

-

Definition at line 49 of file PenningTrap.cpp.

+

Definition at line 43 of file PenningTrap.cpp.

@@ -314,7 +314,7 @@ Private Attributes

Calculate E at point r.

-

Definition at line 38 of file PenningTrap.cpp.

+

Definition at line 32 of file PenningTrap.cpp.

@@ -347,7 +347,7 @@ Private Attributes

Calculate the force between 2 particles.

Calculate the force exhibited on particle p_i from particle p_j.

-

Definition at line 59 of file PenningTrap.cpp.

+

Definition at line 50 of file PenningTrap.cpp.

@@ -366,7 +366,7 @@ Private Attributes
-

Definition at line 142 of file PenningTrap.cpp.

+

Definition at line 188 of file PenningTrap.cpp.

@@ -386,7 +386,7 @@ Private Attributes
-

Definition at line 137 of file PenningTrap.cpp.

+

Definition at line 183 of file PenningTrap.cpp.

@@ -408,7 +408,7 @@ Private Attributes

calculate the total force on a particle.

-

Definition at line 111 of file PenningTrap.cpp.

+

Definition at line 99 of file PenningTrap.cpp.

@@ -431,7 +431,7 @@ Private Attributes

Calculate the total external force on a particle.

Calculate the total amount of force that E and B exhibits on particle p_i.

-

Definition at line 74 of file PenningTrap.cpp.

+

Definition at line 65 of file PenningTrap.cpp.

@@ -453,7 +453,7 @@ Private Attributes

Calculate the total force on a particle from other particles.

-

Definition at line 92 of file PenningTrap.cpp.

+

Definition at line 80 of file PenningTrap.cpp.

@@ -480,7 +480,7 @@ Private Attributes

Magnetic field strength.

-

Definition at line 27 of file PenningTrap.hpp.

+

Definition at line 31 of file PenningTrap.hpp.

@@ -506,7 +506,7 @@ Private Attributes

Characteristic dimension.

-

Definition at line 29 of file PenningTrap.hpp.

+

Definition at line 33 of file PenningTrap.hpp.

@@ -532,7 +532,7 @@ Private Attributes

The particles in the Penning trap.

-

Definition at line 30 of file PenningTrap.hpp.

+

Definition at line 34 of file PenningTrap.hpp.

@@ -558,7 +558,7 @@ Private Attributes

Applied potential.

-

Definition at line 28 of file PenningTrap.hpp.

+

Definition at line 32 of file PenningTrap.hpp.

diff --git a/docs/classPenningTrapTest-members.html b/docs/classPenningTrapTest-members.html new file mode 100644 index 0000000..1b1f2e1 --- /dev/null +++ b/docs/classPenningTrapTest-members.html @@ -0,0 +1,121 @@ + + + + + + + +Penning Trap Simulation: Member List + + + + + + + + + + + + + + + + +
+
+ + + + + + +
+
Penning Trap Simulation +
+
Simulate particle behavior inside a Penning Trap
+
+
+ + + + + + + +
+
+ +
+
+
+ +
+ +
+
+ + +
+
+
+
+
+
Loading...
+
Searching...
+
No Matches
+
+
+
+
+ +
+
PenningTrapTest Member List
+
+
+ +

This is the complete list of members for PenningTrapTest, including all inherited members.

+ + + + + + +
test_external_B_field() (defined in PenningTrapTest)PenningTrapTestinlinestatic
test_external_E_field() (defined in PenningTrapTest)PenningTrapTestinlinestatic
test_force_on_particle() (defined in PenningTrapTest)PenningTrapTestinlinestatic
test_total_force_external() (defined in PenningTrapTest)PenningTrapTestinlinestatic
test_total_force_particles() (defined in PenningTrapTest)PenningTrapTestinlinestatic
+
+ + + + diff --git a/docs/classPenningTrapTest.html b/docs/classPenningTrapTest.html new file mode 100644 index 0000000..8eda947 --- /dev/null +++ b/docs/classPenningTrapTest.html @@ -0,0 +1,273 @@ + + + + + + + +Penning Trap Simulation: PenningTrapTest Class Reference + + + + + + + + + + + + + + + + +
+
+ + + + + + +
+
Penning Trap Simulation +
+
Simulate particle behavior inside a Penning Trap
+
+
+ + + + + + + +
+
+ +
+
+
+ +
+ +
+
+ + +
+
+
+
+
+
Loading...
+
Searching...
+
No Matches
+
+
+
+
+ +
+
+Static Public Member Functions | +List of all members
+
PenningTrapTest Class Reference
+
+
+ + + + + + + + + + + + +

+Static Public Member Functions

static void test_external_E_field ()
 
static void test_external_B_field ()
 
static void test_force_on_particle ()
 
static void test_total_force_external ()
 
static void test_total_force_particles ()
 
+

Detailed Description

+
+

Definition at line 20 of file test_suite.cpp.

+

Member Function Documentation

+ +

◆ test_external_B_field()

+ +
+
+ + + + + +
+ + + + + + + +
static void PenningTrapTest::test_external_B_field ()
+ 		+inlinestatic
+
+ +

Definition at line 59 of file test_suite.cpp.

+ +
+
+ +

◆ test_external_E_field()

+ +
+
+ + + + + +
+ + + + + + + +
static void PenningTrapTest::test_external_E_field ()
+ 		+inlinestatic
+
+ +

Definition at line 22 of file test_suite.cpp.

+ +
+
+ +

◆ test_force_on_particle()

+ +
+
+ + + + + +
+ + + + + + + +
static void PenningTrapTest::test_force_on_particle ()
+ 		+inlinestatic
+
+ +

Definition at line 70 of file test_suite.cpp.

+ +
+
+ +

◆ test_total_force_external()

+ +
+
+ + + + + +
+ + + + + + + +
static void PenningTrapTest::test_total_force_external ()
+ 		+inlinestatic
+
+ +

Definition at line 95 of file test_suite.cpp.

+ +
+
+ +

◆ test_total_force_particles()

+ +
+
+ + + + + +
+ + + + + + + +
static void PenningTrapTest::test_total_force_particles ()
+ 		+inlinestatic
+
+ +

Definition at line 108 of file test_suite.cpp.

+ +
+
+
The documentation for this class was generated from the following file: +
+
+ + + + diff --git a/docs/classes.html b/docs/classes.html index 77925d2..7b8f02e 100644 --- a/docs/classes.html +++ b/docs/classes.html @@ -105,7 +105,7 @@ $(document).ready(function(){initNavTree('classes.html',''); initResizable(); })
P
-
Particle
PenningTrap
+
Particle
PenningTrap
PenningTrapTest
diff --git a/docs/dir_68267d1309a1af8e8297ef4c3efbcdba.html b/docs/dir_68267d1309a1af8e8297ef4c3efbcdba.html index 2f63e67..81efd6b 100644 --- a/docs/dir_68267d1309a1af8e8297ef4c3efbcdba.html +++ b/docs/dir_68267d1309a1af8e8297ef4c3efbcdba.html @@ -104,6 +104,8 @@ $(document).ready(function(){initNavTree('dir_68267d1309a1af8e8297ef4c3efbcdba.h + + @@ -113,8 +115,6 @@ Files - - diff --git a/docs/dir_68267d1309a1af8e8297ef4c3efbcdba.js b/docs/dir_68267d1309a1af8e8297ef4c3efbcdba.js index f72184d..13a8a58 100644 --- a/docs/dir_68267d1309a1af8e8297ef4c3efbcdba.js +++ b/docs/dir_68267d1309a1af8e8297ef4c3efbcdba.js @@ -1,9 +1,9 @@ var dir_68267d1309a1af8e8297ef4c3efbcdba = [ + [ "animate_100_particles.py", "animate__100__particles_8py_source.html", null ], [ "main.cpp", "main_8cpp.html", null ], [ "Particle.cpp", "Particle_8cpp.html", null ], [ "PenningTrap.cpp", "PenningTrap_8cpp.html", null ], - [ "test.py", "test_8py_source.html", null ], - [ "test_suite.cpp", "test__suite_8cpp.html", null ], + [ "test_suite.cpp", "test__suite_8cpp.html", "test__suite_8cpp" ], [ "utils.cpp", "utils_8cpp.html", "utils_8cpp" ] ]; \ No newline at end of file diff --git a/docs/files.html b/docs/files.html index f3ab17f..b07cc8d 100644 --- a/docs/files.html +++ b/docs/files.html @@ -109,10 +109,10 @@ $(document).ready(function(){initNavTree('files.html',''); initResizable(); }); - - - - + + + +

Files

file  animate_100_particles.py [code]
 
file  main.cpp [code]
 The main program for this project.
 
file  PenningTrap.cpp [code]
 The implementation of the PenningTrap class.
 
file  test.py [code]
 
file  test_suite.cpp [code]
 The test suite for the project.
 
 PenningTrap.hppA class for simulating a Penning trap
 utils.hppFunction prototypes and macros that are useful
 ▼ src
 main.cppThe main program for this project
 Particle.cppThe implementation of the Particle class
 PenningTrap.cppThe implementation of the PenningTrap class
 test.py
 animate_100_particles.py
 main.cppThe main program for this project
 Particle.cppThe implementation of the Particle class
 PenningTrap.cppThe implementation of the PenningTrap class
 test_suite.cppThe test suite for the project
 utils.cppImplementation of the utils
diff --git a/docs/globals.html b/docs/globals.html index 78504e0..9107eef 100644 --- a/docs/globals.html +++ b/docs/globals.html @@ -99,8 +99,11 @@ $(document).ready(function(){initNavTree('globals.html',''); initResizable(); })
Here is a list of all documented file members with links to the documentation:
@@ -188,7 +186,7 @@ Janita Ovidie Sandtrøen Willumsen (janitaws)
-

Definition at line 22 of file main.cpp.

+

Definition at line 20 of file main.cpp.

@@ -204,30 +202,11 @@ Janita Ovidie Sandtrøen Willumsen (janitaws)
-

Definition at line 21 of file main.cpp.

+

Definition at line 19 of file main.cpp.

Function Documentation

- -

◆ euler_100_particles()

- -
-
- - - - - - - -
void euler_100_particles ()
-
- -

Definition at line 26 of file main.cpp.

- -
-

◆ main()

@@ -245,6 +224,25 @@ Janita Ovidie Sandtrøen Willumsen (janitaws)

Definition at line 77 of file main.cpp.

+ + + +

◆ simulate_100_particles()

+ +
+
+ + + + + + + +
void simulate_100_particles ()
+
+ +

Definition at line 24 of file main.cpp.

+
diff --git a/docs/main_8cpp_source.html b/docs/main_8cpp_source.html index 3cc7c06..5542b5a 100644 --- a/docs/main_8cpp_source.html +++ b/docs/main_8cpp_source.html @@ -102,65 +102,65 @@ $(document).ready(function(){initNavTree('main_8cpp_source.html',''); initResiza
Go to the documentation of this file.
1
-
13#include <filesystem>
-
14#include <fstream>
-
15#include <string>
-
16#include <omp.h>
-
17#include <sys/stat.h>
+
13#include <fstream>
+
14#include <omp.h>
+
15#include <sys/stat.h>
+
16
+
17#include "PenningTrap.hpp"
18
-
19#include "PenningTrap.hpp"
-
20
-
21#define PARTICLES 100
-
22#define N 10000
-
23#define CHARGE 1.
-
24#define MASS 40. // unit: amu
-
25
-
26void euler_100_particles()
-
27{
-
28 PenningTrap trap;
-
29
-
30 // Add particles inside trap
-
31 for (int i=0; i < PARTICLES; i++) {
-
32 arma::vec r = arma::vec(3).randn() * 0.1 * trap.get_d(); // random initial position
-
33 arma::vec v = arma::vec(3).randn() * 0.1 * trap.get_d(); // random initial velocity
+
19#define PARTICLES 100
+
20#define N 10000
+
21#define CHARGE 1.
+
22#define MASS 40. // unit: amu
+
23
+
24void simulate_100_particles()
+
25{
+
26 PenningTrap trap;
+
27
+
28 // Add particles inside trap
+
29 for (int i = 0; i < PARTICLES; i++) {
+
30 arma::vec r = arma::vec(3).randn() * 0.1 *
+
31 trap.get_d(); // random initial position
+
32 arma::vec v = arma::vec(3).randn() * 0.1 *
+
33 trap.get_d(); // random initial velocity
34 trap.add_particle(Particle(CHARGE, MASS, r, v));
35 }
36
37 double time = 50.; // microseconds
-
38 double dt = time / (double) N;
+
38 double dt = time / (double)N;
39
40 auto res = new arma::vec::fixed<3>[PARTICLES][N];
41
42 int counter = 0;
43
44 // Get the path of all particles
-
45 for (int j=0; j < N; j++) {
-
46 #pragma omp parallel for
-
47 for (int i=0; i < PARTICLES; i++) {
+
45 for (int j = 0; j < N; j++) {
+
46#pragma omp parallel for
+
47 for (int i = 0; i < PARTICLES; i++) {
48 res[i][j] = trap.get_particle(i);
49 }
-
50 trap.evolve_forward_euler(dt);
+
50 trap.evolve_RK4(dt);
51 }
52
53 std::cout << counter << std::endl;
54
-
55 arma::vec::fixed<3>* cur_row;
+
55 arma::vec::fixed<3> *cur_row;
56 arma::vec::fixed<3> cur_elem;
57
58 mkdir("output", 0777);
-
59
-
60 std::ofstream ofile;
-
61
-
62 // Write particle paths to file
-
63 #pragma omp parallel for private(cur_row, cur_elem, ofile)
-
64 for (int i=0; i < PARTICLES; i++) {
-
65 cur_row = res[i];
-
66 ofile.open("output/p" + std::to_string(i) + ".txt");
-
67 for (int j=0; j < N; j++) {
-
68 cur_elem = cur_row[j];
-
69 ofile << cur_elem(0) << ","
-
70 << cur_elem(1) << ","
-
71 << cur_elem(2) << "\n";
+
59 mkdir("output/simulate_100_particles", 0777);
+
60
+
61 std::ofstream ofile;
+
62
+
63// Write particle paths to file
+
64#pragma omp parallel for private(cur_row, cur_elem, ofile)
+
65 for (int i = 0; i < PARTICLES; i++) {
+
66 cur_row = res[i];
+
67 ofile.open("output/simulate_100_particles/p" + std::to_string(i) + ".txt");
+
68 for (int j = 0; j < N; j++) {
+
69 cur_elem = cur_row[j];
+
70 ofile << cur_elem(0) << "," << cur_elem(1) << "," << cur_elem(2)
+
71 << "\n";
72 }
73 ofile.close();
74 }
@@ -170,7 +170,7 @@ $(document).ready(function(){initNavTree('main_8cpp_source.html',''); initResiza
78{
79 double start = omp_get_wtime();
80
-
81 euler_100_particles();
+
81 simulate_100_particles();
82
83 double end = omp_get_wtime();
84
@@ -180,9 +180,9 @@ $(document).ready(function(){initNavTree('main_8cpp_source.html',''); initResiza
88}
PenningTrap.hpp
A class for simulating a Penning trap.
Particle
A class that holds attributes of a particle.
Definition: Particle.hpp:19
-
PenningTrap
A class that simulates a Penning trap.
Definition: PenningTrap.hpp:25
-
PenningTrap::evolve_forward_euler
void evolve_forward_euler(double dt)
Go forward one timestep using the forward Euler method.
Definition: PenningTrap.cpp:121
-
PenningTrap::add_particle
void add_particle(Particle particle)
Add a particle to the system.
Definition: PenningTrap.cpp:33
+
PenningTrap
A class that simulates a Penning trap.
Definition: PenningTrap.hpp:29
+
PenningTrap::add_particle
void add_particle(Particle particle)
Add a particle to the system.
Definition: PenningTrap.cpp:27
+
PenningTrap::evolve_RK4
void evolve_RK4(double dt)
Go forward one timestep using the RK4 method.
Definition: PenningTrap.cpp:104
diff --git a/docs/navtreeindex0.js b/docs/navtreeindex0.js index 75ec38f..27aef05 100644 --- a/docs/navtreeindex0.js +++ b/docs/navtreeindex0.js @@ -1,13 +1,14 @@ var NAVTREEINDEX0 = { -"Particle_8cpp.html":[4,0,1,1], -"Particle_8cpp_source.html":[4,0,1,1], +"Particle_8cpp.html":[4,0,1,2], +"Particle_8cpp_source.html":[4,0,1,2], "Particle_8hpp.html":[4,0,0,1], "Particle_8hpp_source.html":[4,0,0,1], -"PenningTrap_8cpp.html":[4,0,1,2], -"PenningTrap_8cpp_source.html":[4,0,1,2], +"PenningTrap_8cpp.html":[4,0,1,3], +"PenningTrap_8cpp_source.html":[4,0,1,3], "PenningTrap_8hpp.html":[4,0,0,2], "PenningTrap_8hpp_source.html":[4,0,0,2], +"animate__100__particles_8py_source.html":[4,0,1,0], "annotated.html":[3,0], "bug.html":[1], "classParticle.html":[3,0,0], @@ -32,6 +33,7 @@ var NAVTREEINDEX0 = "classPenningTrap.html#a83cc7a04dae009ed8e75453c38fdf8be":[3,0,1,7], "classPenningTrap.html#a978a2e2026a4b5308d5184712f81704b":[3,0,1,5], "classPenningTrap.html#ab891b839c40eeb38d8a2a636c70090a9":[3,0,1,3], +"classPenningTrapTest.html":[3,0,2], "classes.html":[3,1], "constants_8hpp.html":[4,0,0,0], "constants_8hpp.html#a0acb682b8260ab1c60b918599864e2e5":[4,0,0,0,1], @@ -50,20 +52,24 @@ var NAVTREEINDEX0 = "globals_func.html":[4,1,1], "index.html":[], "index.html#autotoc_md1":[0], -"main_8cpp.html":[4,0,1,0], -"main_8cpp_source.html":[4,0,1,0], +"main_8cpp.html":[4,0,1,1], +"main_8cpp_source.html":[4,0,1,1], "pages.html":[], -"test_8py_source.html":[4,0,1,3], "test__suite_8cpp.html":[4,0,1,4], "test__suite_8cpp_source.html":[4,0,1,4], "todo.html":[2], "utils_8cpp.html":[4,0,1,5], -"utils_8cpp.html#a58565270b643b24e3132f38c653e0199":[4,0,1,5,0], -"utils_8cpp.html#acd2a9c7a7d5a7fe9163be8c4cc110746":[4,0,1,5,1], +"utils_8cpp.html#a58565270b643b24e3132f38c653e0199":[4,0,1,5,2], +"utils_8cpp.html#a5d2e1e032fd19614f2fbb58149a7b02a":[4,0,1,5,0], +"utils_8cpp.html#acd2a9c7a7d5a7fe9163be8c4cc110746":[4,0,1,5,3], +"utils_8cpp.html#aff5e07c3c1d321709b0cc38e999f427b":[4,0,1,5,1], "utils_8cpp_source.html":[4,0,1,5], "utils_8hpp.html":[4,0,0,3], -"utils_8hpp.html#ad54b96a1074f9df4dc892a41d115b72d":[4,0,0,3,1], -"utils_8hpp.html#adfb618b2fdff47ef30a4a2b62c04f384":[4,0,0,3,2], -"utils_8hpp.html#aecc1f7a8a2493b9e021e5bff76a00a5b":[4,0,0,3,0], +"utils_8hpp.html#a0c95c4791692b06f8811905a76dbd772":[4,0,0,3,2], +"utils_8hpp.html#a2cc3a2cdb635bac3c8b02e89d4d6af38":[4,0,0,3,3], +"utils_8hpp.html#a73d4f21ad937dbc50a0c0538c78fd4f9":[4,0,0,3,0], +"utils_8hpp.html#ad54b96a1074f9df4dc892a41d115b72d":[4,0,0,3,4], +"utils_8hpp.html#adfb618b2fdff47ef30a4a2b62c04f384":[4,0,0,3,5], +"utils_8hpp.html#aecc1f7a8a2493b9e021e5bff76a00a5b":[4,0,0,3,1], "utils_8hpp_source.html":[4,0,0,3] }; diff --git a/docs/search/all_0.js b/docs/search/all_0.js index 8003164..4fc815a 100644 --- a/docs/search/all_0.js +++ b/docs/search/all_0.js @@ -1,4 +1,6 @@ var searchData= [ - ['add_5fparticle_0',['add_particle',['../classPenningTrap.html#a6e9776ff5b149f01080800716455d7c8',1,'PenningTrap']]] + ['add_5fparticle_0',['add_particle',['../classPenningTrap.html#a6e9776ff5b149f01080800716455d7c8',1,'PenningTrap']]], + ['arma_5fvector_5fclose_5fto_1',['arma_vector_close_to',['../utils_8hpp.html#a0c95c4791692b06f8811905a76dbd772',1,'arma_vector_close_to(arma::vec &a, arma::vec &b, double tol=1e-8): utils.cpp'],['../utils_8cpp.html#a5d2e1e032fd19614f2fbb58149a7b02a',1,'arma_vector_close_to(arma::vec &a, arma::vec &b, double tol): utils.cpp']]], + ['assert_2',['ASSERT',['../utils_8hpp.html#a73d4f21ad937dbc50a0c0538c78fd4f9',1,'utils.hpp']]] ]; diff --git a/docs/search/all_7.js b/docs/search/all_7.js index a26e115..959abdd 100644 --- a/docs/search/all_7.js +++ b/docs/search/all_7.js @@ -1,5 +1,6 @@ var searchData= [ ['m_0',['m',['../classParticle.html#aedcc7e1bc53b0e2b1a4a07c9a1b47563',1,'Particle']]], - ['main_2ecpp_1',['main.cpp',['../main_8cpp.html',1,'']]] + ['m_5fassert_1',['m_assert',['../utils_8hpp.html#a2cc3a2cdb635bac3c8b02e89d4d6af38',1,'m_assert(bool expr, std::string expr_str, std::string func, std::string file, int line, std::string msg): utils.cpp'],['../utils_8cpp.html#aff5e07c3c1d321709b0cc38e999f427b',1,'m_assert(bool expr, std::string expr_str, std::string f, std::string file, int line, std::string msg): utils.cpp']]], + ['main_2ecpp_2',['main.cpp',['../main_8cpp.html',1,'']]] ]; diff --git a/docs/search/all_8.js b/docs/search/all_8.js index c435adb..3ec818e 100644 --- a/docs/search/all_8.js +++ b/docs/search/all_8.js @@ -7,5 +7,6 @@ var searchData= ['penning_20trap_20simulation_4',['Penning trap simulation',['../index.html',1,'']]], ['penningtrap_5',['PenningTrap',['../classPenningTrap.html',1,'PenningTrap'],['../classParticle.html#aa797d319549dc2a0beb06cdbfd430232',1,'Particle::PenningTrap()'],['../classPenningTrap.html#a81add5063bd5f29d77f0b5c299c6c560',1,'PenningTrap::PenningTrap()']]], ['penningtrap_2ecpp_6',['PenningTrap.cpp',['../PenningTrap_8cpp.html',1,'']]], - ['penningtrap_2ehpp_7',['PenningTrap.hpp',['../PenningTrap_8hpp.html',1,'']]] + ['penningtrap_2ehpp_7',['PenningTrap.hpp',['../PenningTrap_8hpp.html',1,'']]], + ['penningtraptest_8',['PenningTrapTest',['../classPenningTrapTest.html',1,'']]] ]; diff --git a/docs/search/classes_0.js b/docs/search/classes_0.js index c907a03..463bbc3 100644 --- a/docs/search/classes_0.js +++ b/docs/search/classes_0.js @@ -1,5 +1,6 @@ var searchData= [ ['particle_0',['Particle',['../classParticle.html',1,'']]], - ['penningtrap_1',['PenningTrap',['../classPenningTrap.html',1,'']]] + ['penningtrap_1',['PenningTrap',['../classPenningTrap.html',1,'']]], + ['penningtraptest_2',['PenningTrapTest',['../classPenningTrapTest.html',1,'']]] ]; diff --git a/docs/search/defines_0.js b/docs/search/defines_0.js index c6466cb..97ce9d9 100644 --- a/docs/search/defines_0.js +++ b/docs/search/defines_0.js @@ -1,4 +1,4 @@ var searchData= [ - ['debug_0',['DEBUG',['../utils_8hpp.html#aecc1f7a8a2493b9e021e5bff76a00a5b',1,'utils.hpp']]] + ['assert_0',['ASSERT',['../utils_8hpp.html#a73d4f21ad937dbc50a0c0538c78fd4f9',1,'utils.hpp']]] ]; diff --git a/docs/search/defines_1.js b/docs/search/defines_1.js index cb8673f..c6466cb 100644 --- a/docs/search/defines_1.js +++ b/docs/search/defines_1.js @@ -1,4 +1,4 @@ var searchData= [ - ['k_5fe_0',['K_E',['../constants_8hpp.html#a4e451456ad7e9276ed0afa42826e7ccb',1,'constants.hpp']]] + ['debug_0',['DEBUG',['../utils_8hpp.html#aecc1f7a8a2493b9e021e5bff76a00a5b',1,'utils.hpp']]] ]; diff --git a/docs/search/defines_2.js b/docs/search/defines_2.js index 6eef01f..cb8673f 100644 --- a/docs/search/defines_2.js +++ b/docs/search/defines_2.js @@ -1,4 +1,4 @@ var searchData= [ - ['t_0',['T',['../constants_8hpp.html#a0acb682b8260ab1c60b918599864e2e5',1,'constants.hpp']]] + ['k_5fe_0',['K_E',['../constants_8hpp.html#a4e451456ad7e9276ed0afa42826e7ccb',1,'constants.hpp']]] ]; diff --git a/docs/search/defines_3.js b/docs/search/defines_3.js index 06f433f..6eef01f 100644 --- a/docs/search/defines_3.js +++ b/docs/search/defines_3.js @@ -1,4 +1,4 @@ var searchData= [ - ['v_0',['V',['../constants_8hpp.html#af40a326b23c68a27cebe60f16634a2cb',1,'constants.hpp']]] + ['t_0',['T',['../constants_8hpp.html#a0acb682b8260ab1c60b918599864e2e5',1,'constants.hpp']]] ]; diff --git a/docs/search/defines_4.js b/docs/search/defines_4.js new file mode 100644 index 0000000..06f433f --- /dev/null +++ b/docs/search/defines_4.js @@ -0,0 +1,4 @@ +var searchData= +[ + ['v_0',['V',['../constants_8hpp.html#af40a326b23c68a27cebe60f16634a2cb',1,'constants.hpp']]] +]; diff --git a/docs/search/functions_0.js b/docs/search/functions_0.js index 8003164..72dba03 100644 --- a/docs/search/functions_0.js +++ b/docs/search/functions_0.js @@ -1,4 +1,5 @@ var searchData= [ - ['add_5fparticle_0',['add_particle',['../classPenningTrap.html#a6e9776ff5b149f01080800716455d7c8',1,'PenningTrap']]] + ['add_5fparticle_0',['add_particle',['../classPenningTrap.html#a6e9776ff5b149f01080800716455d7c8',1,'PenningTrap']]], + ['arma_5fvector_5fclose_5fto_1',['arma_vector_close_to',['../utils_8hpp.html#a0c95c4791692b06f8811905a76dbd772',1,'arma_vector_close_to(arma::vec &a, arma::vec &b, double tol=1e-8): utils.cpp'],['../utils_8cpp.html#a5d2e1e032fd19614f2fbb58149a7b02a',1,'arma_vector_close_to(arma::vec &a, arma::vec &b, double tol): utils.cpp']]] ]; diff --git a/docs/search/functions_3.js b/docs/search/functions_3.js index 2dc79c5..30ce46b 100644 --- a/docs/search/functions_3.js +++ b/docs/search/functions_3.js @@ -1,5 +1,4 @@ var searchData= [ - ['particle_0',['Particle',['../classParticle.html#a00e108823877a25513ccae7cac011b4c',1,'Particle']]], - ['penningtrap_1',['PenningTrap',['../classPenningTrap.html#a81add5063bd5f29d77f0b5c299c6c560',1,'PenningTrap']]] + ['m_5fassert_0',['m_assert',['../utils_8hpp.html#a2cc3a2cdb635bac3c8b02e89d4d6af38',1,'m_assert(bool expr, std::string expr_str, std::string func, std::string file, int line, std::string msg): utils.cpp'],['../utils_8cpp.html#aff5e07c3c1d321709b0cc38e999f427b',1,'m_assert(bool expr, std::string expr_str, std::string f, std::string file, int line, std::string msg): utils.cpp']]] ]; diff --git a/docs/search/functions_4.js b/docs/search/functions_4.js index bf61385..2dc79c5 100644 --- a/docs/search/functions_4.js +++ b/docs/search/functions_4.js @@ -1,4 +1,5 @@ var searchData= [ - ['scientific_5fformat_0',['scientific_format',['../utils_8hpp.html#adfb618b2fdff47ef30a4a2b62c04f384',1,'scientific_format(double d, int width=20, int prec=10): utils.cpp'],['../utils_8hpp.html#ad54b96a1074f9df4dc892a41d115b72d',1,'scientific_format(const std::vector< double > &v, int width=20, int prec=10): utils.cpp'],['../utils_8cpp.html#acd2a9c7a7d5a7fe9163be8c4cc110746',1,'scientific_format(double d, int width, int prec): utils.cpp'],['../utils_8cpp.html#a58565270b643b24e3132f38c653e0199',1,'scientific_format(const std::vector< double > &v, int width, int prec): utils.cpp']]] + ['particle_0',['Particle',['../classParticle.html#a00e108823877a25513ccae7cac011b4c',1,'Particle']]], + ['penningtrap_1',['PenningTrap',['../classPenningTrap.html#a81add5063bd5f29d77f0b5c299c6c560',1,'PenningTrap']]] ]; diff --git a/docs/search/functions_5.js b/docs/search/functions_5.js index 83f26c2..bf61385 100644 --- a/docs/search/functions_5.js +++ b/docs/search/functions_5.js @@ -1,6 +1,4 @@ var searchData= [ - ['total_5fforce_0',['total_force',['../classPenningTrap.html#a83cc7a04dae009ed8e75453c38fdf8be',1,'PenningTrap']]], - ['total_5fforce_5fexternal_1',['total_force_external',['../classPenningTrap.html#a045adb85e97a8e0c38fb36d2fd74eeee',1,'PenningTrap']]], - ['total_5fforce_5fparticles_2',['total_force_particles',['../classPenningTrap.html#a135db7d991d9b88b00fde4edee71f78d',1,'PenningTrap']]] + ['scientific_5fformat_0',['scientific_format',['../utils_8hpp.html#adfb618b2fdff47ef30a4a2b62c04f384',1,'scientific_format(double d, int width=20, int prec=10): utils.cpp'],['../utils_8hpp.html#ad54b96a1074f9df4dc892a41d115b72d',1,'scientific_format(const std::vector< double > &v, int width=20, int prec=10): utils.cpp'],['../utils_8cpp.html#acd2a9c7a7d5a7fe9163be8c4cc110746',1,'scientific_format(double d, int width, int prec): utils.cpp'],['../utils_8cpp.html#a58565270b643b24e3132f38c653e0199',1,'scientific_format(const std::vector< double > &v, int width, int prec): utils.cpp']]] ]; diff --git a/docs/search/functions_6.js b/docs/search/functions_6.js new file mode 100644 index 0000000..83f26c2 --- /dev/null +++ b/docs/search/functions_6.js @@ -0,0 +1,6 @@ +var searchData= +[ + ['total_5fforce_0',['total_force',['../classPenningTrap.html#a83cc7a04dae009ed8e75453c38fdf8be',1,'PenningTrap']]], + ['total_5fforce_5fexternal_1',['total_force_external',['../classPenningTrap.html#a045adb85e97a8e0c38fb36d2fd74eeee',1,'PenningTrap']]], + ['total_5fforce_5fparticles_2',['total_force_particles',['../classPenningTrap.html#a135db7d991d9b88b00fde4edee71f78d',1,'PenningTrap']]] +]; diff --git a/docs/search/searchdata.js b/docs/search/searchdata.js index ca77c16..906104d 100644 --- a/docs/search/searchdata.js +++ b/docs/search/searchdata.js @@ -3,10 +3,10 @@ var indexSectionsWithContent = 0: "abcdefkmpqrstuv", 1: "p", 2: "cmptu", - 3: "aefpst", + 3: "aefmpst", 4: "bdmpqrv", 5: "p", - 6: "dktv", + 6: "adktv", 7: "bpt" }; diff --git a/docs/test__suite_8cpp.html b/docs/test__suite_8cpp.html index 2bab26f..787939b 100644 --- a/docs/test__suite_8cpp.html +++ b/docs/test__suite_8cpp.html @@ -99,6 +99,7 @@ $(document).ready(function(){initNavTree('test__suite_8cpp.html',''); initResiza
+Classes | Functions
test_suite.cpp File Reference
@@ -106,9 +107,19 @@ $(document).ready(function(){initNavTree('test__suite_8cpp.html',''); initResiza

The test suite for the project. More...

- +
#include "PenningTrap.hpp"
+#include "utils.hpp"
+#include <iomanip>
+#include <sstream>
+#include <string>
+

Go to the source code of this file.

+ + + +

+Classes

class  PenningTrapTest
 
@@ -139,7 +150,7 @@ Janita Ovidie Sandtrøen Willumsen (janitaws)

Functions

int main ()
-

Definition at line 12 of file test_suite.cpp.

+

Definition at line 135 of file test_suite.cpp.

diff --git a/docs/test__suite_8cpp.js b/docs/test__suite_8cpp.js new file mode 100644 index 0000000..5bace7d --- /dev/null +++ b/docs/test__suite_8cpp.js @@ -0,0 +1,4 @@ +var test__suite_8cpp = +[ + [ "PenningTrapTest", "classPenningTrapTest.html", null ] +]; \ No newline at end of file diff --git a/docs/test__suite_8cpp_source.html b/docs/test__suite_8cpp_source.html index 056e5c3..3fab92c 100644 --- a/docs/test__suite_8cpp_source.html +++ b/docs/test__suite_8cpp_source.html @@ -102,10 +102,151 @@ $(document).ready(function(){initNavTree('test__suite_8cpp_source.html',''); ini
Go to the documentation of this file.
1
-
12int main()
-
13{
-
14 return 0;
-
15}
+
13#include "PenningTrap.hpp"
+
14#include "utils.hpp"
+
15
+
16#include <iomanip>
+
17#include <sstream>
+
18#include <string>
+
19
+
20class PenningTrapTest {
+
21public:
+
22 static void test_external_E_field()
+
23 {
+
24 PenningTrap trap;
+
25
+
26 // Vector containing inputs and expected results
+
27 std::vector<std::pair<arma::vec, arma::vec>> tests;
+
28
+
29 tests.push_back(
+
30 std::make_pair(arma::vec{0., 0., 0.}, arma::vec{0., 0., 0.}));
+
31
+
32 tests.push_back(std::make_pair(arma::vec{10., 0., 0.},
+
33 arma::vec{96.4852558, 0., 0.}));
+
34
+
35 tests.push_back(std::make_pair(arma::vec{10., 0., 0.},
+
36 arma::vec{96.4852558, 0., 0.}));
+
37
+
38 tests.push_back(std::make_pair(arma::vec{0., 10., 0.},
+
39 arma::vec{0., 96.4852558, 0.}));
+
40
+
41 tests.push_back(std::make_pair(arma::vec{0., 0., 10.},
+
42 arma::vec{0., 0., -192.9705116}));
+
43
+
44 arma::vec result;
+
45 arma::vec v;
+
46 std::stringstream msg;
+
47 for (int i = 0; i < tests.size(); i++) {
+
48 v = tests.at(i).first;
+
49 result = trap.external_E_field(v);
+
50
+
51 msg.str("");
+
52 msg << "Testing the external E field at (" << std::setprecision(2)
+
53 << v(0) << "," << v(1) << "," << v(2) << ").";
+
54
+
55 ASSERT(arma_vector_close_to(result, tests.at(i).second), msg.str());
+
56 }
+
57 }
+
58
+
59 static void test_external_B_field()
+
60 {
+
61 // No point in testing at different points since it's not dependent
+
62 // on position.
+
63 PenningTrap trap;
+
64 arma::vec expected{0., 0., T};
+
65 arma::vec result = trap.external_B_field(arma::vec{0., 0., 0.});
+
66 ASSERT(arma_vector_close_to(expected, result),
+
67 "Testing the external B field at (0,0,0)");
+
68 }
+
69
+
70 static void test_force_on_particle()
+
71 {
+
72 PenningTrap trap;
+
73 arma::vec v{0., 0., 0.};
+
74
+
75 // Add particles to test
+
76 trap.add_particle(Particle(1., 40., arma::vec{0., 0., 0.}, v));
+
77 trap.add_particle(Particle(1., 40., arma::vec{1., 0., 0.}, v));
+
78 trap.add_particle(Particle(1., 40., arma::vec{0., 3., 4.}, v));
+
79
+
80 // Test p0 and p1
+
81 arma::vec expected{-1., 0., 0.};
+
82 arma::vec result = trap.force_on_particle(0, 1);
+
83 ASSERT(arma_vector_close_to(expected, result),
+
84 "Testing the force on a particle at (0,0,0) from a "
+
85 "particle at (1,0,0).");
+
86
+
87 // Test p0 and p2
+
88 expected = arma::vec{0, -.024, -.032};
+
89 result = trap.force_on_particle(0, 2);
+
90 ASSERT(arma_vector_close_to(expected, result),
+
91 "Testing the force on a particle at (0,0,0) from a "
+
92 "particle at (0,3,4).");
+
93 }
+
94
+
95 static void test_total_force_external()
+
96 {
+
97 PenningTrap trap;
+
98 trap.add_particle(
+
99 Particle(1., 40., arma::vec{1., 2., 3.}, arma::vec{3., 4., 5.}));
+
100
+
101 arma::vec expected{395.58954878, -270.15871624, -57.89115348};
+
102 arma::vec result = trap.total_force_external(0);
+
103 ASSERT(arma_vector_close_to(expected, result),
+
104 "Testing the total external force on a particle at "
+
105 "(1,2,3) with velocity (3,4,5)");
+
106 }
+
107
+
108 static void test_total_force_particles()
+
109 {
+
110 PenningTrap trap;
+
111 trap.add_particle(
+
112 Particle(1., 40., arma::vec{0., 0., 0.}, arma::vec{0., 0., 0.}));
+
113
+
114 arma::vec expected{0., 0., 0.};
+
115 arma::vec result = trap.total_force_particles(0);
+
116 ASSERT(arma_vector_close_to(expected, result),
+
117 "Testing the total force of all particles on particle 0 "
+
118 "with only a single particle");
+
119
+
120 trap.add_particle(
+
121 Particle(1., 40., arma::vec{1., 0., 0.}, arma::vec{0., 0., 0.}));
+
122 trap.add_particle(
+
123 Particle(1., 40., arma::vec{0., 1., 0.}, arma::vec{0., 0., 0.}));
+
124 trap.add_particle(
+
125 Particle(1., 40., arma::vec{0., 0., 1.}, arma::vec{0., 0., 0.}));
+
126
+
127 expected = arma::vec(3, arma::fill::value(-3473.383325));
+
128 result = trap.total_force_particles(0);
+
129 ASSERT(arma_vector_close_to(expected, result),
+
130 "Testing the total force of all particles on particle 0 "
+
131 "with 3 other particles.");
+
132 }
+
133};
+
134
+
135int main()
+
136{
+
137 PenningTrapTest::test_external_E_field();
+
138 PenningTrapTest::test_external_B_field();
+
139 PenningTrapTest::test_force_on_particle();
+
140 PenningTrapTest::test_total_force_external();
+
141 PenningTrapTest::test_total_force_particles();
+
142 return 0;
+
143}
+
PenningTrap.hpp
A class for simulating a Penning trap.
+
Particle
A class that holds attributes of a particle.
Definition: Particle.hpp:19
+ +
PenningTrap
A class that simulates a Penning trap.
Definition: PenningTrap.hpp:29
+
PenningTrap::total_force_external
arma::vec total_force_external(int i)
Calculate the total external force on a particle.
Definition: PenningTrap.cpp:65
+
PenningTrap::total_force_particles
arma::vec total_force_particles(int i)
Calculate the total force on a particle from other particles.
Definition: PenningTrap.cpp:80
+
PenningTrap::external_B_field
arma::vec external_B_field(arma::vec r)
Calculate B at point r.
Definition: PenningTrap.cpp:43
+
PenningTrap::force_on_particle
arma::vec force_on_particle(int i, int j)
Calculate the force between 2 particles.
Definition: PenningTrap.cpp:50
+
PenningTrap::add_particle
void add_particle(Particle particle)
Add a particle to the system.
Definition: PenningTrap.cpp:27
+
PenningTrap::external_E_field
arma::vec external_E_field(arma::vec r)
Calculate E at point r.
Definition: PenningTrap.cpp:32
+
T
#define T
1 Tesla. unit:
Definition: constants.hpp:17
+
utils.hpp
Function prototypes and macros that are useful.
+
arma_vector_close_to
bool arma_vector_close_to(arma::vec &a, arma::vec &b, double tol=1e-8)
Test if two armadillo vectors are close to each other.
Definition: utils.cpp:59
+
ASSERT
#define ASSERT(expr, msg)
A prettier assertion function.
Definition: utils.hpp:45
diff --git a/docs/utils_8cpp.html b/docs/utils_8cpp.html index 74ad93f..c2ffd2e 100644 --- a/docs/utils_8cpp.html +++ b/docs/utils_8cpp.html @@ -118,6 +118,12 @@ Functions std::string scientific_format (const std::vector< double > &v, int width, int prec)  Turns a vector of doubles into a string written in scientific format.
  +void m_assert (bool expr, std::string expr_str, std::string f, std::string file, int line, std::string msg) + Test an expression, confirm that test is ok, or abort execution.
+  +bool arma_vector_close_to (arma::vec &a, arma::vec &b, double tol) + Test if two armadillo vectors are close to each other.
+ 

Detailed Description

Implementation of the utils.

@@ -129,6 +135,122 @@ Janita Ovidie Sandtrøen Willumsen (janitaws)

Definition in file utils.cpp.

Function Documentation

+ +

◆ arma_vector_close_to()

+ +
+
+ + + + + + + + + + + + + + + + + + + + + + + + +
bool arma_vector_close_to (arma::vec & a,
arma::vec & b,
double tol = 1e-8 
)
+
+ +

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.

+
Parameters
+ + + + +
aVector a
bVector b
tolThe tolerance
+
+
+
Returns
Boolean
+ +

Definition at line 59 of file utils.cpp.

+ +
+
+ +

◆ m_assert()

+ +
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
void m_assert (bool expr,
std::string expr_str,
std::string func,
std::string file,
int line,
std::string msg 
)
+
+ +

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.

+
Parameters
+ + + + + + + +
exprThe expression to be evaluated
expr_strThe stringified version of the expression
funcThe function name of the caller
fileThe file of the caller
lineThe line number where this function is called from
msgThe message to be displayed
+
+
+ +

Definition at line 40 of file utils.cpp.

+ +
+

◆ scientific_format() [1/2]

diff --git a/docs/utils_8cpp.js b/docs/utils_8cpp.js index f2eda7a..36c50fa 100644 --- a/docs/utils_8cpp.js +++ b/docs/utils_8cpp.js @@ -1,5 +1,7 @@ var utils_8cpp = [ + [ "arma_vector_close_to", "utils_8cpp.html#a5d2e1e032fd19614f2fbb58149a7b02a", null ], + [ "m_assert", "utils_8cpp.html#aff5e07c3c1d321709b0cc38e999f427b", null ], [ "scientific_format", "utils_8cpp.html#a58565270b643b24e3132f38c653e0199", null ], [ "scientific_format", "utils_8cpp.html#acd2a9c7a7d5a7fe9163be8c4cc110746", null ] ]; \ No newline at end of file diff --git a/docs/utils_8cpp_source.html b/docs/utils_8cpp_source.html index 11f53d7..069cd7c 100644 --- a/docs/utils_8cpp_source.html +++ b/docs/utils_8cpp_source.html @@ -111,15 +111,60 @@ $(document).ready(function(){initNavTree('utils_8cpp_source.html',''); initResiz
18 return ss.str();
19}
20
-
21std::string scientific_format(const std::vector<double>& v, int width, int prec)
+
21std::string scientific_format(const std::vector<double> &v, int width, int prec)
22{
23 std::stringstream ss;
-
24 for(double elem : v) {
+
24 for (double elem : v) {
25 ss << scientific_format(elem, width, prec);
26 }
27 return ss.str();
28}
+
29
+
30static void print_message(std::string msg)
+
31{
+
32 if (msg.size() > 0) {
+
33 std::cout << "message: " << msg << "\n\n";
+
34 }
+
35 else {
+
36 std::cout << "\n";
+
37 }
+
38}
+
39
+
40void m_assert(bool expr, std::string expr_str, std::string f, std::string file,
+
41 int line, std::string msg)
+
42{
+
43 std::string new_assert(f.size() + (expr ? 4 : 6), '-');
+
44 std::cout << "\x1B[36m" << new_assert << "\033[0m\n";
+
45 std::cout << f << ": ";
+
46 if (expr) {
+
47 std::cout << "\x1B[32mOK\033[0m\n";
+
48 print_message(msg);
+
49 }
+
50 else {
+
51 std::cout << "\x1B[31mFAIL\033[0m\n";
+
52 print_message(msg);
+
53 std::cout << file << " " << line << ": Assertion \"" << expr_str
+
54 << "\" Failed\n\n";
+
55 abort();
+
56 }
+
57}
+
58
+
59bool arma_vector_close_to(arma::vec &a, arma::vec &b, double tol)
+
60{
+
61 if (a.n_elem != b.n_elem) {
+
62 return false;
+
63 }
+
64
+
65 for (int i = 0; i < a.n_elem; i++) {
+
66 if (std::abs(a(i) - b(i)) >= tol) {
+
67 return false;
+
68 }
+
69 }
+
70 return true;
+
71}
+
arma_vector_close_to
bool arma_vector_close_to(arma::vec &a, arma::vec &b, double tol)
Test if two armadillo vectors are close to each other.
Definition: utils.cpp:59
scientific_format
std::string scientific_format(double d, int width, int prec)
Turns a double into a string written in scientific format.
Definition: utils.cpp:14
+
m_assert
void m_assert(bool expr, std::string expr_str, std::string f, std::string file, int line, std::string msg)
Test an expression, confirm that test is ok, or abort execution.
Definition: utils.cpp:40
utils.hpp
Function prototypes and macros that are useful.
diff --git a/docs/utils_8hpp.html b/docs/utils_8hpp.html index 0789c66..dea4601 100644 --- a/docs/utils_8hpp.html +++ b/docs/utils_8hpp.html @@ -107,10 +107,11 @@ $(document).ready(function(){initNavTree('utils_8hpp.html',''); initResizable();

Function prototypes and macros that are useful. More...

-
#include <string>
-#include <vector>
+
#include <armadillo>
#include <iomanip>
#include <sstream>
+#include <string>
+#include <vector>

Go to the source code of this file.

@@ -119,6 +120,11 @@ Macros + + + + +
#define DEBUG(msg)
 Writes a debug message.
 
#define ASSERT(expr, msg)
 A prettier assertion function.
 
#define __METHOD_NAME__   methodName(__PRETTY_FUNCTION__)
 
@@ -128,6 +134,12 @@ Functions + + + + + +

Functions

std::string scientific_format (const std::vector< double > &v, int width=20, int prec=10)
 Turns a vector of doubles into a string written in scientific format.
 
void m_assert (bool expr, std::string expr_str, std::string func, std::string file, int line, std::string msg)
 Test an expression, confirm that test is ok, or abort execution.
 
bool arma_vector_close_to (arma::vec &a, arma::vec &b, double tol=1e-8)
 Test if two armadillo vectors are close to each other.
 

Detailed Description

Function prototypes and macros that are useful.

@@ -140,6 +152,58 @@ Janita Ovidie Sandtrøen Willumsen (janitaws)

Definition in file utils.hpp.

Macro Definition Documentation

+ +

◆ __METHOD_NAME__

+ +
+
+ + + + +
#define __METHOD_NAME__   methodName(__PRETTY_FUNCTION__)
+
+ +

Definition at line 48 of file utils.hpp.

+ +
+
+ +

◆ ASSERT

+ +
+
+ + + + + + + + + + + + + + + + + + +
#define ASSERT( expr,
 msg 
)
+
+Value:
m_assert(expr, #expr, __METHOD_NAME__, __FILE__, \
+
__LINE__, msg)
+
m_assert
void m_assert(bool expr, std::string expr_str, std::string func, std::string file, int line, std::string msg)
Test an expression, confirm that test is ok, or abort execution.
Definition: utils.cpp:40
+
+

A prettier assertion function.

+

This macro calls the m_assert function which is a more informative assertion function than the regular assert function from cassert.

+ +

Definition at line 45 of file utils.hpp.

+ +
+

◆ DEBUG

@@ -157,13 +221,129 @@ Janita Ovidie Sandtrøen Willumsen (janitaws)

Writes a debug message.

-

This function writes a debug message that includes the filename, line number, and a custom message. The function is wrapped in an ifdef that checks if DBG is defined, so one can choose to display the debug messages by adding the -DDBG flag when compiling.

+

This macro writes a debug message that includes the filename, line number, and a custom message. The function is wrapped in an ifdef that checks if DBG is defined, so one can choose to display the debug messages by adding the -DDBG flag when compiling.

-

Definition at line 35 of file utils.hpp.

+

Definition at line 36 of file utils.hpp.

Function Documentation

+ +

◆ arma_vector_close_to()

+ +
+
+ + + + + + + + + + + + + + + + + + + + + + + + +
bool arma_vector_close_to (arma::vec & a,
arma::vec & b,
double tol = 1e-8 
)
+
+ +

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.

+
Parameters
+ + + + +
aVector a
bVector b
tolThe tolerance
+
+
+
Returns
Boolean
+ +

Definition at line 59 of file utils.cpp.

+ +
+
+ +

◆ m_assert()

+ +
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
void m_assert (bool expr,
std::string expr_str,
std::string func,
std::string file,
int line,
std::string msg 
)
+
+ +

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.

+
Parameters
+ + + + + + + +
exprThe expression to be evaluated
expr_strThe stringified version of the expression
funcThe function name of the caller
fileThe file of the caller
lineThe line number where this function is called from
msgThe message to be displayed
+
+
+ +

Definition at line 40 of file utils.cpp.

+ +
+

◆ scientific_format() [1/2]

diff --git a/docs/utils_8hpp.js b/docs/utils_8hpp.js index 7b91aee..34c7b36 100644 --- a/docs/utils_8hpp.js +++ b/docs/utils_8hpp.js @@ -1,6 +1,9 @@ var utils_8hpp = [ + [ "ASSERT", "utils_8hpp.html#a73d4f21ad937dbc50a0c0538c78fd4f9", null ], [ "DEBUG", "utils_8hpp.html#aecc1f7a8a2493b9e021e5bff76a00a5b", null ], + [ "arma_vector_close_to", "utils_8hpp.html#a0c95c4791692b06f8811905a76dbd772", null ], + [ "m_assert", "utils_8hpp.html#a2cc3a2cdb635bac3c8b02e89d4d6af38", null ], [ "scientific_format", "utils_8hpp.html#ad54b96a1074f9df4dc892a41d115b72d", null ], [ "scientific_format", "utils_8hpp.html#adfb618b2fdff47ef30a4a2b62c04f384", null ] ]; \ No newline at end of file diff --git a/docs/utils_8hpp_source.html b/docs/utils_8hpp_source.html index 5959341..33ddc52 100644 --- a/docs/utils_8hpp_source.html +++ b/docs/utils_8hpp_source.html @@ -105,25 +105,54 @@ $(document).ready(function(){initNavTree('utils_8hpp_source.html',''); initResiz
15#ifndef __UTILS__
16#define __UTILS__
17
-
18#include <string>
-
19#include <vector>
-
20#include <iomanip>
-
21#include <sstream>
-
22
-
31#ifdef DBG
-
32 #define DEBUG(msg) std::cout << __FILE__ << " " << __LINE__ << ": " \
-
33 << msg << std::endl
-
34#else
-
35 #define DEBUG(msg)
-
36#endif
-
37
-
51std::string scientific_format(double d, int width=20, int prec=10);
-
52
-
61std::string scientific_format(const std::vector<double>& v,
-
62 int width=20,
-
63 int prec=10);
-
64
-
65#endif
+
18#include <armadillo>
+
19#include <iomanip>
+
20#include <sstream>
+
21#include <string>
+
22#include <vector>
+
23
+
32#ifdef DBG
+
33 #define DEBUG(msg) std::cout << __FILE__ << " " << __LINE__ << ": " \
+
34 << msg << std::endl
+
35#else
+
36 #define DEBUG(msg)
+
37#endif
+
38
+
45#define ASSERT(expr, msg) m_assert(expr, #expr, __METHOD_NAME__, __FILE__, \
+
46 __LINE__, msg)
+
47
+
48#define __METHOD_NAME__ methodName(__PRETTY_FUNCTION__)
+
49
+
50
+
64std::string scientific_format(double d, int width=20, int prec=10);
+
65
+
74std::string scientific_format(const std::vector<double>& v,
+
75 int width=20,
+
76 int prec=10);
+
77
+
90void m_assert(bool expr,
+
91 std::string expr_str,
+
92 std::string func,
+
93 std::string file,
+
94 int line,
+
95 std::string msg);
+
96
+
97
+
109bool arma_vector_close_to(arma::vec &a, arma::vec &b, double tol=1e-8);
+
110
+
111
+
112static inline std::string methodName(const std::string& prettyFunction)
+
113{
+
114 size_t colons = prettyFunction.find("::");
+
115 size_t begin = prettyFunction.substr(0,colons).rfind(" ") + 1;
+
116 size_t end = prettyFunction.rfind("(") - begin;
+
117
+
118 return prettyFunction.substr(begin,end) + "()";
+
119}
+
120
+
121#endif
+
arma_vector_close_to
bool arma_vector_close_to(arma::vec &a, arma::vec &b, double tol=1e-8)
Test if two armadillo vectors are close to each other.
Definition: utils.cpp:59
+
m_assert
void m_assert(bool expr, std::string expr_str, std::string func, std::string file, int line, std::string msg)
Test an expression, confirm that test is ok, or abort execution.
Definition: utils.cpp:40
scientific_format
std::string scientific_format(double d, int width=20, int prec=10)
Turns a double into a string written in scientific format.
Definition: utils.cpp:14
diff --git a/images/100_particles.gif b/images/100_particles.gif new file mode 100644 index 0000000..53eb275 Binary files /dev/null and b/images/100_particles.gif differ diff --git a/include/PenningTrap.hpp b/include/PenningTrap.hpp index 1b1752e..3601634 100644 --- a/include/PenningTrap.hpp +++ b/include/PenningTrap.hpp @@ -13,9 +13,13 @@ #define __PENNING_TRAP__ #include +#include -#include "constants.hpp" #include "Particle.hpp" +#include "constants.hpp" + +#pragma omp declare reduction( + : arma::vec : omp_out += omp_in ) \ + initializer( omp_priv = omp_orig ) /** @brief A class that simulates a Penning trap. * diff --git a/include/utils.hpp b/include/utils.hpp index c8d3e4a..7a03149 100644 --- a/include/utils.hpp +++ b/include/utils.hpp @@ -15,15 +15,16 @@ #ifndef __UTILS__ #define __UTILS__ -#include -#include +#include #include #include +#include +#include /** @def DEBUG(msg) * @brief Writes a debug message * - * This function writes a debug message that includes the filename, + * This macro writes a debug message that includes the filename, * line number, and a custom message. The function is wrapped in an ifdef * that checks if DBG is defined, so one can choose to display the debug * messages by adding the -DDBG flag when compiling. @@ -35,6 +36,18 @@ #define DEBUG(msg) #endif +/** @def ASSERT(expr) + * @brief A prettier assertion function. + * + * This macro calls the m_assert function which is a more informative + * assertion function than the regular assert function from cassert. + * */ +#define ASSERT(expr, msg) m_assert(expr, #expr, __METHOD_NAME__, __FILE__, \ + __LINE__, msg) + +#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 @@ -62,4 +75,47 @@ 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. + * + * @param expr The expression to be evaluated + * @param expr_str The stringified version of the expression + * @param func The function name of the caller + * @param file The file of the caller + * @param line The line number where this function is called from + * @param msg The message to be displayed + * */ +void m_assert(bool expr, + std::string expr_str, + std::string func, + std::string file, + int line, + std::string msg); + + +/** @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. + * + * @param a Vector a + * @param b Vector b + * @param tol The tolerance + * + * @return Boolean + * */ +bool arma_vector_close_to(arma::vec &a, arma::vec &b, double tol=1e-8); + + +static inline std::string methodName(const std::string& prettyFunction) +{ + size_t colons = prettyFunction.find("::"); + size_t begin = prettyFunction.substr(0,colons).rfind(" ") + 1; + size_t end = prettyFunction.rfind("(") - begin; + + return prettyFunction.substr(begin,end) + "()"; +} + #endif diff --git a/man_pages/man3/Particle.3 b/man_pages/man3/Particle.3 index db9056a..f37639d 100644 --- a/man_pages/man3/Particle.3 +++ b/man_pages/man3/Particle.3 @@ -1,4 +1,4 @@ -.TH "Particle" 3 "Mon Oct 2 2023" "Penning Trap Simulation" \" -*- nroff -*- +.TH "Particle" 3 "Sun Oct 8 2023" "Penning Trap Simulation" \" -*- nroff -*- .ad l .nh .SH NAME diff --git a/man_pages/man3/Particle.cpp.3 b/man_pages/man3/Particle.cpp.3 index 2799e38..7a93b45 100644 --- a/man_pages/man3/Particle.cpp.3 +++ b/man_pages/man3/Particle.cpp.3 @@ -1,4 +1,4 @@ -.TH "src/Particle.cpp" 3 "Mon Oct 2 2023" "Penning Trap Simulation" \" -*- nroff -*- +.TH "src/Particle.cpp" 3 "Sun Oct 8 2023" "Penning Trap Simulation" \" -*- nroff -*- .ad l .nh .SH NAME diff --git a/man_pages/man3/Particle.hpp.3 b/man_pages/man3/Particle.hpp.3 index 2a5a85e..abb26a7 100644 --- a/man_pages/man3/Particle.hpp.3 +++ b/man_pages/man3/Particle.hpp.3 @@ -1,4 +1,4 @@ -.TH "include/Particle.hpp" 3 "Mon Oct 2 2023" "Penning Trap Simulation" \" -*- nroff -*- +.TH "include/Particle.hpp" 3 "Sun Oct 8 2023" "Penning Trap Simulation" \" -*- nroff -*- .ad l .nh .SH NAME diff --git a/man_pages/man3/PenningTrap.3 b/man_pages/man3/PenningTrap.3 index 42c26cf..2087bb5 100644 --- a/man_pages/man3/PenningTrap.3 +++ b/man_pages/man3/PenningTrap.3 @@ -1,4 +1,4 @@ -.TH "PenningTrap" 3 "Mon Oct 2 2023" "Penning Trap Simulation" \" -*- nroff -*- +.TH "PenningTrap" 3 "Sun Oct 8 2023" "Penning Trap Simulation" \" -*- nroff -*- .ad l .nh .SH NAME @@ -85,7 +85,7 @@ 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\&. .PP -Definition at line \fB25\fP of file \fBPenningTrap\&.hpp\fP\&. +Definition at line \fB29\fP of file \fBPenningTrap\&.hpp\fP\&. .SH "Constructor & Destructor Documentation" .PP .SS "PenningTrap::PenningTrap (double B_0 = \fC\fBT\fP\fP, double V_0 = \fC25\&.*\fBV\fP/1000\&.\fP, double d = \fC500\&.\fP)" @@ -93,7 +93,7 @@ Definition at line \fB25\fP of file \fBPenningTrap\&.hpp\fP\&. .PP Set B_0, V_0 and d\&. .PP -Definition at line \fB26\fP of file \fBPenningTrap\&.cpp\fP\&. +Definition at line \fB20\fP of file \fBPenningTrap\&.cpp\fP\&. .SH "Member Function Documentation" .PP .SS "void PenningTrap::add_particle (\fBParticle\fP particle)" @@ -101,63 +101,63 @@ Definition at line \fB26\fP of file \fBPenningTrap\&.cpp\fP\&. .PP Add a particle to the system\&. .PP -Definition at line \fB33\fP of file \fBPenningTrap\&.cpp\fP\&. +Definition at line \fB27\fP of file \fBPenningTrap\&.cpp\fP\&. .SS "void PenningTrap::evolve_forward_euler (double dt)" .PP Go forward one timestep using the forward Euler method\&. .PP -Definition at line \fB121\fP of file \fBPenningTrap\&.cpp\fP\&. +Definition at line \fB167\fP of file \fBPenningTrap\&.cpp\fP\&. .SS "void PenningTrap::evolve_RK4 (double dt)" .PP Go forward one timestep using the RK4 method\&. .PP -Definition at line \fB116\fP of file \fBPenningTrap\&.cpp\fP\&. +Definition at line \fB104\fP of file \fBPenningTrap\&.cpp\fP\&. .SS "arma::vec PenningTrap::external_B_field (arma::vec r)" .PP Calculate B at point r\&. .PP -Definition at line \fB49\fP of file \fBPenningTrap\&.cpp\fP\&. +Definition at line \fB43\fP of file \fBPenningTrap\&.cpp\fP\&. .SS "arma::vec PenningTrap::external_E_field (arma::vec r)" .PP Calculate E at point r\&. .PP -Definition at line \fB38\fP of file \fBPenningTrap\&.cpp\fP\&. +Definition at line \fB32\fP of file \fBPenningTrap\&.cpp\fP\&. .SS "arma::vec PenningTrap::force_on_particle (int i, int j)" .PP Calculate the force between 2 particles\&. Calculate the force exhibited on particle p_i from particle p_j\&. .PP -Definition at line \fB59\fP of file \fBPenningTrap\&.cpp\fP\&. +Definition at line \fB50\fP of file \fBPenningTrap\&.cpp\fP\&. .SS "double PenningTrap::get_d ()" .PP -Definition at line \fB142\fP of file \fBPenningTrap\&.cpp\fP\&. +Definition at line \fB188\fP of file \fBPenningTrap\&.cpp\fP\&. .SS "arma::vec PenningTrap::get_particle (int i)" .PP -Definition at line \fB137\fP of file \fBPenningTrap\&.cpp\fP\&. +Definition at line \fB183\fP of file \fBPenningTrap\&.cpp\fP\&. .SS "arma::vec PenningTrap::total_force (int i)" .PP calculate the total force on a particle\&. .PP -Definition at line \fB111\fP of file \fBPenningTrap\&.cpp\fP\&. +Definition at line \fB99\fP of file \fBPenningTrap\&.cpp\fP\&. .SS "arma::vec PenningTrap::total_force_external (int i)" .PP Calculate the total external force on a particle\&. Calculate the total amount of force that E and B exhibits on particle p_i\&. .PP -Definition at line \fB74\fP of file \fBPenningTrap\&.cpp\fP\&. +Definition at line \fB65\fP of file \fBPenningTrap\&.cpp\fP\&. .SS "arma::vec PenningTrap::total_force_particles (int i)" .PP Calculate the total force on a particle from other particles\&. .PP -Definition at line \fB92\fP of file \fBPenningTrap\&.cpp\fP\&. +Definition at line \fB80\fP of file \fBPenningTrap\&.cpp\fP\&. .SH "Member Data Documentation" .PP .SS "double PenningTrap::B_0\fC [private]\fP" @@ -165,25 +165,25 @@ Definition at line \fB92\fP of file \fBPenningTrap\&.cpp\fP\&. .PP Magnetic field strength\&. .PP -Definition at line \fB27\fP of file \fBPenningTrap\&.hpp\fP\&. +Definition at line \fB31\fP of file \fBPenningTrap\&.hpp\fP\&. .SS "double PenningTrap::d\fC [private]\fP" .PP Characteristic dimension\&. .PP -Definition at line \fB29\fP of file \fBPenningTrap\&.hpp\fP\&. +Definition at line \fB33\fP of file \fBPenningTrap\&.hpp\fP\&. .SS "std::vector<\fBParticle\fP> PenningTrap::particles\fC [private]\fP" .PP The particles in the Penning trap\&. .PP -Definition at line \fB30\fP of file \fBPenningTrap\&.hpp\fP\&. +Definition at line \fB34\fP of file \fBPenningTrap\&.hpp\fP\&. .SS "double PenningTrap::V_0\fC [private]\fP" .PP Applied potential\&. .PP -Definition at line \fB28\fP of file \fBPenningTrap\&.hpp\fP\&. +Definition at line \fB32\fP of file \fBPenningTrap\&.hpp\fP\&. .SH "Author" .PP diff --git a/man_pages/man3/PenningTrap.cpp.3 b/man_pages/man3/PenningTrap.cpp.3 index fecf53b..f850299 100644 --- a/man_pages/man3/PenningTrap.cpp.3 +++ b/man_pages/man3/PenningTrap.cpp.3 @@ -1,4 +1,4 @@ -.TH "src/PenningTrap.cpp" 3 "Mon Oct 2 2023" "Penning Trap Simulation" \" -*- nroff -*- +.TH "src/PenningTrap.cpp" 3 "Sun Oct 8 2023" "Penning Trap Simulation" \" -*- nroff -*- .ad l .nh .SH NAME @@ -7,17 +7,11 @@ src/PenningTrap.cpp \- The implementation of the \fBPenningTrap\fP class\&. .SH SYNOPSIS .br .PP -\fC#include 'utils\&.hpp'\fP -.br \fC#include 'PenningTrap\&.hpp'\fP .br \fC#include 'constants\&.hpp'\fP .br -\fC#include \fP -.br -\fC#include \fP -.br -\fC#include \fP +\fC#include 'utils\&.hpp'\fP .br .SH "Detailed Description" diff --git a/man_pages/man3/PenningTrap.hpp.3 b/man_pages/man3/PenningTrap.hpp.3 index 4d434f2..c16a94b 100644 --- a/man_pages/man3/PenningTrap.hpp.3 +++ b/man_pages/man3/PenningTrap.hpp.3 @@ -1,4 +1,4 @@ -.TH "include/PenningTrap.hpp" 3 "Mon Oct 2 2023" "Penning Trap Simulation" \" -*- nroff -*- +.TH "include/PenningTrap.hpp" 3 "Sun Oct 8 2023" "Penning Trap Simulation" \" -*- nroff -*- .ad l .nh .SH NAME @@ -9,10 +9,12 @@ include/PenningTrap.hpp \- A class for simulating a Penning trap\&. .PP \fC#include \fP .br -\fC#include 'constants\&.hpp'\fP +\fC#include \fP .br \fC#include 'Particle\&.hpp'\fP .br +\fC#include 'constants\&.hpp'\fP +.br .SS "Classes" diff --git a/man_pages/man3/PenningTrapTest.3 b/man_pages/man3/PenningTrapTest.3 new file mode 100644 index 0000000..07c0b75 --- /dev/null +++ b/man_pages/man3/PenningTrapTest.3 @@ -0,0 +1,56 @@ +.TH "PenningTrapTest" 3 "Sun Oct 8 2023" "Penning Trap Simulation" \" -*- nroff -*- +.ad l +.nh +.SH NAME +PenningTrapTest +.SH SYNOPSIS +.br +.PP +.SS "Static Public Member Functions" + +.in +1c +.ti -1c +.RI "static void \fBtest_external_E_field\fP ()" +.br +.ti -1c +.RI "static void \fBtest_external_B_field\fP ()" +.br +.ti -1c +.RI "static void \fBtest_force_on_particle\fP ()" +.br +.ti -1c +.RI "static void \fBtest_total_force_external\fP ()" +.br +.ti -1c +.RI "static void \fBtest_total_force_particles\fP ()" +.br +.in -1c +.SH "Detailed Description" +.PP +Definition at line \fB20\fP of file \fBtest_suite\&.cpp\fP\&. +.SH "Member Function Documentation" +.PP +.SS "static void PenningTrapTest::test_external_B_field ()\fC [inline]\fP, \fC [static]\fP" + +.PP +Definition at line \fB59\fP of file \fBtest_suite\&.cpp\fP\&. +.SS "static void PenningTrapTest::test_external_E_field ()\fC [inline]\fP, \fC [static]\fP" + +.PP +Definition at line \fB22\fP of file \fBtest_suite\&.cpp\fP\&. +.SS "static void PenningTrapTest::test_force_on_particle ()\fC [inline]\fP, \fC [static]\fP" + +.PP +Definition at line \fB70\fP of file \fBtest_suite\&.cpp\fP\&. +.SS "static void PenningTrapTest::test_total_force_external ()\fC [inline]\fP, \fC [static]\fP" + +.PP +Definition at line \fB95\fP of file \fBtest_suite\&.cpp\fP\&. +.SS "static void PenningTrapTest::test_total_force_particles ()\fC [inline]\fP, \fC [static]\fP" + +.PP +Definition at line \fB108\fP of file \fBtest_suite\&.cpp\fP\&. + +.SH "Author" +.PP +Generated automatically by Doxygen for Penning Trap Simulation from the source code\&. diff --git a/man_pages/man3/bug.3 b/man_pages/man3/bug.3 index 0362f46..f93d6e3 100644 --- a/man_pages/man3/bug.3 +++ b/man_pages/man3/bug.3 @@ -1,4 +1,4 @@ -.TH "bug" 3 "Mon Oct 2 2023" "Penning Trap Simulation" \" -*- nroff -*- +.TH "bug" 3 "Sun Oct 8 2023" "Penning Trap Simulation" \" -*- nroff -*- .ad l .nh .SH NAME diff --git a/man_pages/man3/constants.hpp.3 b/man_pages/man3/constants.hpp.3 index 2002a1b..ed6ff65 100644 --- a/man_pages/man3/constants.hpp.3 +++ b/man_pages/man3/constants.hpp.3 @@ -1,4 +1,4 @@ -.TH "include/constants.hpp" 3 "Mon Oct 2 2023" "Penning Trap Simulation" \" -*- nroff -*- +.TH "include/constants.hpp" 3 "Sun Oct 8 2023" "Penning Trap Simulation" \" -*- nroff -*- .ad l .nh .SH NAME diff --git a/man_pages/man3/main.cpp.3 b/man_pages/man3/main.cpp.3 index 4cab7a6..e02bbe4 100644 --- a/man_pages/man3/main.cpp.3 +++ b/man_pages/man3/main.cpp.3 @@ -1,4 +1,4 @@ -.TH "src/main.cpp" 3 "Mon Oct 2 2023" "Penning Trap Simulation" \" -*- nroff -*- +.TH "src/main.cpp" 3 "Sun Oct 8 2023" "Penning Trap Simulation" \" -*- nroff -*- .ad l .nh .SH NAME @@ -7,12 +7,8 @@ src/main.cpp \- The main program for this project\&. .SH SYNOPSIS .br .PP -\fC#include \fP -.br \fC#include \fP .br -\fC#include \fP -.br \fC#include \fP .br \fC#include \fP @@ -40,7 +36,7 @@ src/main.cpp \- The main program for this project\&. .in +1c .ti -1c -.RI "void \fBeuler_100_particles\fP ()" +.RI "void \fBsimulate_100_particles\fP ()" .br .ti -1c .RI "int \fBmain\fP ()" @@ -78,29 +74,29 @@ Definition in file \fBmain\&.cpp\fP\&. .SS "#define CHARGE 1\&." .PP -Definition at line \fB23\fP of file \fBmain\&.cpp\fP\&. +Definition at line \fB21\fP of file \fBmain\&.cpp\fP\&. .SS "#define MASS 40\&." .PP -Definition at line \fB24\fP of file \fBmain\&.cpp\fP\&. +Definition at line \fB22\fP of file \fBmain\&.cpp\fP\&. .SS "#define N 10000" .PP -Definition at line \fB22\fP of file \fBmain\&.cpp\fP\&. +Definition at line \fB20\fP of file \fBmain\&.cpp\fP\&. .SS "#define PARTICLES 100" .PP -Definition at line \fB21\fP of file \fBmain\&.cpp\fP\&. +Definition at line \fB19\fP of file \fBmain\&.cpp\fP\&. .SH "Function Documentation" .PP -.SS "void euler_100_particles ()" - -.PP -Definition at line \fB26\fP of file \fBmain\&.cpp\fP\&. .SS "int main ()" .PP Definition at line \fB77\fP of file \fBmain\&.cpp\fP\&. +.SS "void simulate_100_particles ()" + +.PP +Definition at line \fB24\fP of file \fBmain\&.cpp\fP\&. .SH "Author" .PP Generated automatically by Doxygen for Penning Trap Simulation from the source code\&. diff --git a/man_pages/man3/test_suite.cpp.3 b/man_pages/man3/test_suite.cpp.3 index 35cfa4b..fa522ef 100644 --- a/man_pages/man3/test_suite.cpp.3 +++ b/man_pages/man3/test_suite.cpp.3 @@ -1,4 +1,4 @@ -.TH "src/test_suite.cpp" 3 "Mon Oct 2 2023" "Penning Trap Simulation" \" -*- nroff -*- +.TH "src/test_suite.cpp" 3 "Sun Oct 8 2023" "Penning Trap Simulation" \" -*- nroff -*- .ad l .nh .SH NAME @@ -7,6 +7,24 @@ src/test_suite.cpp \- The test suite for the project\&. .SH SYNOPSIS .br .PP +\fC#include 'PenningTrap\&.hpp'\fP +.br +\fC#include 'utils\&.hpp'\fP +.br +\fC#include \fP +.br +\fC#include \fP +.br +\fC#include \fP +.br + +.SS "Classes" + +.in +1c +.ti -1c +.RI "class \fBPenningTrapTest\fP" +.br +.in -1c .SS "Functions" .in +1c @@ -46,7 +64,7 @@ Definition in file \fBtest_suite\&.cpp\fP\&. .SS "int main ()" .PP -Definition at line \fB12\fP of file \fBtest_suite\&.cpp\fP\&. +Definition at line \fB135\fP of file \fBtest_suite\&.cpp\fP\&. .SH "Author" .PP Generated automatically by Doxygen for Penning Trap Simulation from the source code\&. diff --git a/man_pages/man3/todo.3 b/man_pages/man3/todo.3 index c342666..a66e064 100644 --- a/man_pages/man3/todo.3 +++ b/man_pages/man3/todo.3 @@ -1,4 +1,4 @@ -.TH "todo" 3 "Mon Oct 2 2023" "Penning Trap Simulation" \" -*- nroff -*- +.TH "todo" 3 "Sun Oct 8 2023" "Penning Trap Simulation" \" -*- nroff -*- .ad l .nh .SH NAME diff --git a/man_pages/man3/utils.cpp.3 b/man_pages/man3/utils.cpp.3 index 910bd0a..3d60a86 100644 --- a/man_pages/man3/utils.cpp.3 +++ b/man_pages/man3/utils.cpp.3 @@ -1,4 +1,4 @@ -.TH "src/utils.cpp" 3 "Mon Oct 2 2023" "Penning Trap Simulation" \" -*- nroff -*- +.TH "src/utils.cpp" 3 "Sun Oct 8 2023" "Penning Trap Simulation" \" -*- nroff -*- .ad l .nh .SH NAME @@ -21,6 +21,14 @@ src/utils.cpp \- Implementation of the utils\&. .RI "std::string \fBscientific_format\fP (const std::vector< double > &v, int width, int prec)" .br .RI "Turns a vector of doubles into a string written in scientific format\&. " +.ti -1c +.RI "void \fBm_assert\fP (bool expr, std::string expr_str, std::string f, std::string file, int line, std::string msg)" +.br +.RI "Test an expression, confirm that test is ok, or abort execution\&. " +.ti -1c +.RI "bool \fBarma_vector_close_to\fP (arma::vec &a, arma::vec &b, double tol)" +.br +.RI "Test if two armadillo vectors are close to each other\&. " .in -1c .SH "Detailed Description" .PP @@ -51,6 +59,51 @@ No known bugs Definition in file \fButils\&.cpp\fP\&. .SH "Function Documentation" .PP +.SS "bool arma_vector_close_to (arma::vec & a, arma::vec & b, double tol = \fC1e\-8\fP)" + +.PP +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\&. +.PP +\fBParameters\fP +.RS 4 +\fIa\fP Vector a +.br +\fIb\fP Vector b +.br +\fItol\fP The tolerance +.RE +.PP +\fBReturns\fP +.RS 4 +Boolean +.RE +.PP + +.PP +Definition at line \fB59\fP of file \fButils\&.cpp\fP\&. +.SS "void m_assert (bool expr, std::string expr_str, std::string func, std::string file, int line, std::string msg)" + +.PP +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\&. +.PP +\fBParameters\fP +.RS 4 +\fIexpr\fP The expression to be evaluated +.br +\fIexpr_str\fP The stringified version of the expression +.br +\fIfunc\fP The function name of the caller +.br +\fIfile\fP The file of the caller +.br +\fIline\fP The line number where this function is called from +.br +\fImsg\fP The message to be displayed +.RE +.PP + +.PP +Definition at line \fB40\fP of file \fButils\&.cpp\fP\&. .SS "std::string scientific_format (const std::vector< double > & v, int width = \fC20\fP, int prec = \fC10\fP)" .PP diff --git a/man_pages/man3/utils.hpp.3 b/man_pages/man3/utils.hpp.3 index 76556bd..a67beb7 100644 --- a/man_pages/man3/utils.hpp.3 +++ b/man_pages/man3/utils.hpp.3 @@ -1,4 +1,4 @@ -.TH "include/utils.hpp" 3 "Mon Oct 2 2023" "Penning Trap Simulation" \" -*- nroff -*- +.TH "include/utils.hpp" 3 "Sun Oct 8 2023" "Penning Trap Simulation" \" -*- nroff -*- .ad l .nh .SH NAME @@ -7,14 +7,16 @@ include/utils.hpp \- Function prototypes and macros that are useful\&. .SH SYNOPSIS .br .PP -\fC#include \fP -.br -\fC#include \fP +\fC#include \fP .br \fC#include \fP .br \fC#include \fP .br +\fC#include \fP +.br +\fC#include \fP +.br .SS "Macros" @@ -23,6 +25,13 @@ include/utils.hpp \- Function prototypes and macros that are useful\&. .RI "#define \fBDEBUG\fP(msg)" .br .RI "Writes a debug message\&. " +.ti -1c +.RI "#define \fBASSERT\fP(expr, msg)" +.br +.RI "A prettier assertion function\&. " +.ti -1c +.RI "#define \fB__METHOD_NAME__\fP methodName(__PRETTY_FUNCTION__)" +.br .in -1c .SS "Functions" @@ -35,6 +44,14 @@ include/utils.hpp \- Function prototypes and macros that are useful\&. .RI "std::string \fBscientific_format\fP (const std::vector< double > &v, int width=20, int prec=10)" .br .RI "Turns a vector of doubles into a string written in scientific format\&. " +.ti -1c +.RI "void \fBm_assert\fP (bool expr, std::string expr_str, std::string func, std::string file, int line, std::string msg)" +.br +.RI "Test an expression, confirm that test is ok, or abort execution\&. " +.ti -1c +.RI "bool \fBarma_vector_close_to\fP (arma::vec &a, arma::vec &b, double tol=1e\-8)" +.br +.RI "Test if two armadillo vectors are close to each other\&. " .in -1c .SH "Detailed Description" .PP @@ -67,14 +84,74 @@ No known bugs Definition in file \fButils\&.hpp\fP\&. .SH "Macro Definition Documentation" .PP +.SS "#define __METHOD_NAME__ methodName(__PRETTY_FUNCTION__)" + +.PP +Definition at line \fB48\fP of file \fButils\&.hpp\fP\&. +.SS "#define ASSERT(expr, msg)" +\fBValue:\fP.PP +.nf + m_assert(expr, #expr, __METHOD_NAME__, __FILE__, \\ + __LINE__, msg) +.fi + +.PP +A prettier assertion function\&. This macro calls the m_assert function which is a more informative assertion function than the regular assert function from cassert\&. +.PP +Definition at line \fB45\fP of file \fButils\&.hpp\fP\&. .SS "#define DEBUG(msg)" .PP -Writes a debug message\&. This function writes a debug message that includes the filename, line number, and a custom message\&. The function is wrapped in an ifdef that checks if DBG is defined, so one can choose to display the debug messages by adding the -DDBG flag when compiling\&. +Writes a debug message\&. This macro writes a debug message that includes the filename, line number, and a custom message\&. The function is wrapped in an ifdef that checks if DBG is defined, so one can choose to display the debug messages by adding the -DDBG flag when compiling\&. .PP -Definition at line \fB35\fP of file \fButils\&.hpp\fP\&. +Definition at line \fB36\fP of file \fButils\&.hpp\fP\&. .SH "Function Documentation" .PP +.SS "bool arma_vector_close_to (arma::vec & a, arma::vec & b, double tol = \fC1e\-8\fP)" + +.PP +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\&. +.PP +\fBParameters\fP +.RS 4 +\fIa\fP Vector a +.br +\fIb\fP Vector b +.br +\fItol\fP The tolerance +.RE +.PP +\fBReturns\fP +.RS 4 +Boolean +.RE +.PP + +.PP +Definition at line \fB59\fP of file \fButils\&.cpp\fP\&. +.SS "void m_assert (bool expr, std::string expr_str, std::string func, std::string file, int line, std::string msg)" + +.PP +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\&. +.PP +\fBParameters\fP +.RS 4 +\fIexpr\fP The expression to be evaluated +.br +\fIexpr_str\fP The stringified version of the expression +.br +\fIfunc\fP The function name of the caller +.br +\fIfile\fP The file of the caller +.br +\fIline\fP The line number where this function is called from +.br +\fImsg\fP The message to be displayed +.RE +.PP + +.PP +Definition at line \fB40\fP of file \fButils\&.cpp\fP\&. .SS "std::string scientific_format (const std::vector< double > & v, int width = \fC20\fP, int prec = \fC10\fP)" .PP diff --git a/src/Makefile b/src/Makefile index f0cd401..52c16b7 100644 --- a/src/Makefile +++ b/src/Makefile @@ -27,7 +27,7 @@ all: test_suite main main: main.o $(LIBOBJS) $(CLASSOBJS) $(CC) $^ -o $@ $(CFLAGS) $(DBGFLAG) -I$(INCLUDE) $(OPENMP) -test_suite: test_suite.o $(LIBOBJS) +test_suite: test_suite.o $(LIBOBJS) $(CLASSOBJS) $(CC) $^ -o $@ $(CFLAGS) $(DBGFLAG) -I$(INCLUDE) $(OPENMP) # Rule for object files diff --git a/src/PenningTrap.cpp b/src/PenningTrap.cpp index 59be75c..5d9a16b 100644 --- a/src/PenningTrap.cpp +++ b/src/PenningTrap.cpp @@ -13,15 +13,9 @@ * @todo Implement evolve_forward_euler * */ -#include "utils.hpp" #include "PenningTrap.hpp" #include "constants.hpp" -#include -#include -#include - -#pragma omp declare reduction( + : arma::vec : omp_out += omp_in ) \ - initializer( omp_priv = omp_orig ) +#include "utils.hpp" PenningTrap::PenningTrap(double B_0, double V_0, double d) { @@ -38,20 +32,17 @@ void PenningTrap::add_particle(Particle particle) arma::vec PenningTrap::external_E_field(arma::vec r) { arma::vec::fixed<3> res; - double f = this->V_0/(this->d*this->d); + double f = this->V_0 / (this->d * this->d); res(0) = r(0); res(1) = r(1); - res(2) = -2.*r(2); + res(2) = -2. * r(2); - return f*res; + return f * res; } arma::vec PenningTrap::external_B_field(arma::vec r) { - arma::vec::fixed<3> res; - res(0) = 0.; - res(1) = 0.; - res(2) = this->B_0; + arma::vec::fixed<3> res{0., 0., this->B_0}; return res; } @@ -59,14 +50,14 @@ arma::vec PenningTrap::external_B_field(arma::vec r) arma::vec PenningTrap::force_on_particle(int i, int j) { // Calculate the difference between the particles' position - arma::vec::fixed<3> res = this->particles.at(i).r_vec - - this->particles.at(j).r_vec; + arma::vec::fixed<3> res = + this->particles.at(i).r_vec - this->particles.at(j).r_vec; // Get the distance between the particles double norm = arma::norm(res); // Multiply res with p_j's charge divided by the norm cubed - res *= this->particles.at(j).q/(norm*norm*norm); + res *= this->particles.at(j).q / (norm * norm * norm); return res; } @@ -75,16 +66,13 @@ arma::vec PenningTrap::total_force_external(int i) { Particle p = this->particles.at(i); - arma::vec::fixed<3> v_cross_B; - arma::vec::fixed<3> B = this->external_B_field(p.r_vec); - v_cross_B(0) = p.v_vec(1)*B(2) - p.v_vec(2)*B(1); - v_cross_B(1) = p.v_vec(2)*B(0) - p.v_vec(0)*B(2); - v_cross_B(2) = p.v_vec(0)*B(1) - p.v_vec(1)*B(0); + arma::vec::fixed<3> 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)}; - arma::vec force = p.q - *(this->external_E_field(p.r_vec) + v_cross_B); + arma::vec force = p.q * (this->external_E_field(p.r_vec) + v_cross_B); return force; } @@ -95,7 +83,7 @@ arma::vec PenningTrap::total_force_particles(int i) arma::vec res(3); - for (int j=0; j < this->particles.size(); j++) { + for (int j = 0; j < this->particles.size(); j++) { if (i == j) { continue; } @@ -103,7 +91,7 @@ arma::vec PenningTrap::total_force_particles(int i) res += this->force_on_particle(i, j); } - res *= K_E*(p.q/p.m); + res *= K_E * (p.q / p.m); return res; } @@ -115,7 +103,65 @@ arma::vec PenningTrap::total_force(int i) void PenningTrap::evolve_RK4(double dt) { + std::vector tmp_particles = this->particles; + + arma::vec::fixed<3> *k_v = new arma::vec::fixed<3>[this->particles.size()*4]; + arma::vec::fixed<3> *k_r = new arma::vec::fixed<3>[this->particles.size()*4]; + int size = this->particles.size(); + + for (int i=0; itotal_force(i)/this->particles.at(i).m; + k_r[i] = this->particles.at(i).v_vec; + } + + for (int i=0; iparticles.at(i); + + p->v_vec = tmp_particles.at(i).v_vec + (dt/2)*k_v[i]; + p->r_vec = tmp_particles.at(i).r_vec + (dt/2)*k_r[i]; + } + + + for (int i=0; itotal_force(i)/this->particles.at(i).m; + k_r[1*size + i] = this->particles.at(i).v_vec; + } + + for (int i=0; iparticles.at(i); + + p->v_vec = tmp_particles.at(i).v_vec + (dt/2)*k_v[1*size + i]; + p->r_vec = tmp_particles.at(i).r_vec + (dt/2)*k_r[1*size + i]; + } + + for (int i=0; itotal_force(i)/this->particles.at(i).m; + k_r[2*size + i] = this->particles.at(i).v_vec; + } + + for (int i=0; iparticles.at(i); + + p->v_vec = tmp_particles.at(i).v_vec + dt*k_v[2*size + i]; + p->r_vec = tmp_particles.at(i).r_vec + dt*k_r[2*size + i]; + } + + + for (int i=0; itotal_force(i)/this->particles.at(i).m; + k_r[3*size + i] = this->particles.at(i).v_vec; + } + + for (int i=0; iparticles.at(i); + + p->v_vec = tmp_particles.at(i).v_vec + dt*(k_v[i] + k_v[size + i] + k_v[2*size + i] + k_v[3*size + i])/6; + p->r_vec = tmp_particles.at(i).r_vec + dt*(k_r[i] + k_r[size + i] + k_r[2*size + i] + k_r[3*size + i])/6; + } + + delete [] k_v; + delete [] k_r; } void PenningTrap::evolve_forward_euler(double dt) @@ -124,11 +170,11 @@ void PenningTrap::evolve_forward_euler(double dt) Particle *p; - #pragma omp parallel for private(p) - for (int i=0; i < this->particles.size(); i++) { +#pragma omp parallel for private(p) + for (int i = 0; i < this->particles.size(); i++) { p = &new_state.at(i); - p->v_vec += dt*this->total_force(i)/new_state.at(i).m; - p->r_vec += dt*this->particles.at(i).v_vec; + p->v_vec += dt * this->total_force(i) / new_state.at(i).m; + p->r_vec += dt * this->particles.at(i).v_vec; } this->particles = new_state; diff --git a/src/animate_100_particles.py b/src/animate_100_particles.py new file mode 100644 index 0000000..19f10ca --- /dev/null +++ b/src/animate_100_particles.py @@ -0,0 +1,71 @@ +import matplotlib.pyplot as plt +import numpy as np +from matplotlib import animation +from mpl_toolkits.mplot3d import Axes3D + + +def get_data(files): + res = [] + for file in files: + arr = [[], [], []] + with open(file, encoding="utf8") as f: + lines = f.readlines() + + for line in lines: + xi, yi, zi = map(float, line.strip().split(",")) + arr[0].append(xi) + arr[1].append(yi) + arr[2].append(zi) + res.append(arr) + + return np.array(res) + + +def update(num, lines, arr): + for line, a in zip(lines, arr): + line.set_data(a[:2, num - 1 : num]) + line.set_3d_properties(a[2, num]) + + +def animate(): + plt.style.use("dark_background") + fig = plt.figure() + ax = fig.add_subplot(projection="3d") + + arr = get_data([f"output/simulate_100_particles/p{i}.txt" for i in range(100)]) + + arr = arr[:, :, ::10] + + N = len(arr[0][0]) + + lines = [ax.plot(*a[:, 1], "o")[0] for a in arr] + + ax.set_title("100 particles inside a Penning trap") + plt.figtext( + 0.5, + 0.01, + "100 randomly generated particles " + "evolving over a time of 50 microseconds.", + fontsize=12, + horizontalalignment="center", + ) + + ax.set_xlim3d([-500.0, 500.0]) + ax.set_xlabel("X (micrometers)") + + ax.set_ylim3d([-500.0, 500.0]) + ax.set_ylabel("Y (micrometers)") + + ax.set_zlim3d([-500.0, 500.0]) + ax.set_zlabel("Z (micrometers)") + + ani = animation.FuncAnimation( + fig, update, N, fargs=(lines, arr), interval=1, blit=False + ) + + ani.save("../images/100_particles.gif", writer=animation.FFMpegFileWriter(fps=50)) + # plt.show() + + +if __name__ == "__main__": + animate() diff --git a/src/main.cpp b/src/main.cpp index b479649..1c55279 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -10,9 +10,7 @@ * @bug No known bugs * */ -#include #include -#include #include #include @@ -23,52 +21,54 @@ #define CHARGE 1. #define MASS 40. // unit: amu -void euler_100_particles() +void simulate_100_particles() { PenningTrap trap; // Add particles inside trap - for (int i=0; i < PARTICLES; i++) { - arma::vec r = arma::vec(3).randn() * 0.1 * trap.get_d(); // random initial position - arma::vec v = arma::vec(3).randn() * 0.1 * trap.get_d(); // random initial velocity + for (int i = 0; i < PARTICLES; i++) { + arma::vec r = arma::vec(3).randn() * 0.1 * + trap.get_d(); // random initial position + arma::vec v = arma::vec(3).randn() * 0.1 * + trap.get_d(); // random initial velocity trap.add_particle(Particle(CHARGE, MASS, r, v)); } double time = 50.; // microseconds - double dt = time / (double) N; + double dt = time / (double)N; auto res = new arma::vec::fixed<3>[PARTICLES][N]; int counter = 0; // Get the path of all particles - for (int j=0; j < N; j++) { - #pragma omp parallel for - for (int i=0; i < PARTICLES; i++) { + for (int j = 0; j < N; j++) { +#pragma omp parallel for + for (int i = 0; i < PARTICLES; i++) { res[i][j] = trap.get_particle(i); } - trap.evolve_forward_euler(dt); + trap.evolve_RK4(dt); } std::cout << counter << std::endl; - arma::vec::fixed<3>* cur_row; + arma::vec::fixed<3> *cur_row; arma::vec::fixed<3> cur_elem; mkdir("output", 0777); + mkdir("output/simulate_100_particles", 0777); std::ofstream ofile; - // Write particle paths to file - #pragma omp parallel for private(cur_row, cur_elem, ofile) - for (int i=0; i < PARTICLES; i++) { +// Write particle paths to file +#pragma omp parallel for private(cur_row, cur_elem, ofile) + for (int i = 0; i < PARTICLES; i++) { cur_row = res[i]; - ofile.open("output/p" + std::to_string(i) + ".txt"); - for (int j=0; j < N; j++) { + ofile.open("output/simulate_100_particles/p" + std::to_string(i) + ".txt"); + for (int j = 0; j < N; j++) { cur_elem = cur_row[j]; - ofile << cur_elem(0) << "," - << cur_elem(1) << "," - << cur_elem(2) << "\n"; + ofile << cur_elem(0) << "," << cur_elem(1) << "," << cur_elem(2) + << "\n"; } ofile.close(); } @@ -78,7 +78,7 @@ int main() { double start = omp_get_wtime(); - euler_100_particles(); + simulate_100_particles(); double end = omp_get_wtime(); diff --git a/src/test.py b/src/test.py deleted file mode 100644 index a215cad..0000000 --- a/src/test.py +++ /dev/null @@ -1,63 +0,0 @@ -import matplotlib.pyplot as plt -import numpy as np -from mpl_toolkits.mplot3d import Axes3D -from matplotlib import animation - -def get_data(files): - res = [] - for file in files: - arr = [[], [], []] - with open(file) as f: - lines = f.readlines() - - for line in lines: - xi,yi,zi = map(lambda x: float(x), line.strip().split(",")) - arr[0].append(xi) - arr[1].append(yi) - arr[2].append(zi) - res.append(arr) - - return np.array(res) - -def update(num, lines, arr): - for line, a in zip(lines, arr): - line.set_data(a[:2, num]) - line.set_3d_properties(a[2, num]) - - - -def animate(): - fig = plt.figure() - ax = fig.add_subplot(projection="3d") - - - arr = get_data([f"output/p{i}.txt" for i in range(100)]) - - arr = arr[:,:,::10] - - N = len(arr[0][0]) - - lines = [ax.plot(*a[:,1], "o")[0] for a in arr] - - ax.set_xlim3d([-500.0, 500.0]) - ax.set_xlabel('X') - - ax.set_ylim3d([-500.0, 500.0]) - ax.set_ylabel('Y') - - ax.set_zlim3d([-500.0, 500.0]) - ax.set_zlabel('Z') - - ani = animation.FuncAnimation(fig, update, N, fargs=(lines, arr), - interval=1, - blit=False) - - - # ani.save("100_particles.gif", writer=animation.FFMpegFileWriter(fps=30)) - plt.show() - - -if __name__ == "__main__": - animate() - - diff --git a/src/test_suite.cpp b/src/test_suite.cpp index 41f8b9c..81ea5c1 100644 --- a/src/test_suite.cpp +++ b/src/test_suite.cpp @@ -9,7 +9,135 @@ * * @bug No known bugs * */ -int main() + +#include "PenningTrap.hpp" +#include "utils.hpp" + +#include +#include +#include + +class PenningTrapTest { +public: + static void test_external_E_field() + { + PenningTrap trap; + + // Vector containing inputs and expected results + std::vector> tests; + + tests.push_back( + std::make_pair(arma::vec{0., 0., 0.}, arma::vec{0., 0., 0.})); + + tests.push_back(std::make_pair(arma::vec{10., 0., 0.}, + arma::vec{96.4852558, 0., 0.})); + + tests.push_back(std::make_pair(arma::vec{10., 0., 0.}, + arma::vec{96.4852558, 0., 0.})); + + tests.push_back(std::make_pair(arma::vec{0., 10., 0.}, + arma::vec{0., 96.4852558, 0.})); + + tests.push_back(std::make_pair(arma::vec{0., 0., 10.}, + arma::vec{0., 0., -192.9705116})); + + arma::vec result; + arma::vec v; + std::stringstream msg; + for (int i = 0; i < tests.size(); i++) { + v = tests.at(i).first; + result = trap.external_E_field(v); + + msg.str(""); + msg << "Testing the external E field at (" << std::setprecision(2) + << v(0) << "," << v(1) << "," << v(2) << ")."; + + ASSERT(arma_vector_close_to(result, tests.at(i).second), msg.str()); + } + } + + static void test_external_B_field() + { + // No point in testing at different points since it's not dependent + // on position. + PenningTrap trap; + arma::vec expected{0., 0., T}; + arma::vec result = trap.external_B_field(arma::vec{0., 0., 0.}); + ASSERT(arma_vector_close_to(expected, result), + "Testing the external B field at (0,0,0)"); + } + + static void test_force_on_particle() + { + PenningTrap trap; + arma::vec v{0., 0., 0.}; + + // Add particles to test + trap.add_particle(Particle(1., 40., arma::vec{0., 0., 0.}, v)); + trap.add_particle(Particle(1., 40., arma::vec{1., 0., 0.}, v)); + trap.add_particle(Particle(1., 40., arma::vec{0., 3., 4.}, v)); + + // Test p0 and p1 + arma::vec expected{-1., 0., 0.}; + arma::vec result = trap.force_on_particle(0, 1); + ASSERT(arma_vector_close_to(expected, result), + "Testing the force on a particle at (0,0,0) from a " + "particle at (1,0,0)."); + + // Test p0 and p2 + expected = arma::vec{0, -.024, -.032}; + result = trap.force_on_particle(0, 2); + ASSERT(arma_vector_close_to(expected, result), + "Testing the force on a particle at (0,0,0) from a " + "particle at (0,3,4)."); + } + + static void test_total_force_external() + { + PenningTrap trap; + trap.add_particle( + Particle(1., 40., arma::vec{1., 2., 3.}, arma::vec{3., 4., 5.})); + + arma::vec expected{395.58954878, -270.15871624, -57.89115348}; + arma::vec result = trap.total_force_external(0); + ASSERT(arma_vector_close_to(expected, result), + "Testing the total external force on a particle at " + "(1,2,3) with velocity (3,4,5)"); + } + + static void test_total_force_particles() + { + PenningTrap trap; + trap.add_particle( + Particle(1., 40., arma::vec{0., 0., 0.}, arma::vec{0., 0., 0.})); + + arma::vec expected{0., 0., 0.}; + arma::vec result = trap.total_force_particles(0); + ASSERT(arma_vector_close_to(expected, result), + "Testing the total force of all particles on particle 0 " + "with only a single particle"); + + trap.add_particle( + Particle(1., 40., arma::vec{1., 0., 0.}, arma::vec{0., 0., 0.})); + trap.add_particle( + Particle(1., 40., arma::vec{0., 1., 0.}, arma::vec{0., 0., 0.})); + trap.add_particle( + Particle(1., 40., arma::vec{0., 0., 1.}, arma::vec{0., 0., 0.})); + + expected = arma::vec(3, arma::fill::value(-3473.383325)); + result = trap.total_force_particles(0); + ASSERT(arma_vector_close_to(expected, result), + "Testing the total force of all particles on particle 0 " + "with 3 other particles."); + } +}; + +int main() { + PenningTrapTest::test_external_E_field(); + PenningTrapTest::test_external_B_field(); + PenningTrapTest::test_force_on_particle(); + PenningTrapTest::test_total_force_external(); + PenningTrapTest::test_total_force_particles(); return 0; } diff --git a/src/utils.cpp b/src/utils.cpp index 0e93a31..b55fa34 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -18,11 +18,54 @@ std::string scientific_format(double d, int width, int prec) return ss.str(); } -std::string scientific_format(const std::vector& v, int width, int prec) +std::string scientific_format(const std::vector &v, int width, int prec) { std::stringstream ss; - for(double elem : v) { + for (double elem : v) { ss << scientific_format(elem, width, prec); } return ss.str(); } + +static void print_message(std::string msg) +{ + if (msg.size() > 0) { + std::cout << "message: " << msg << "\n\n"; + } + else { + std::cout << "\n"; + } +} + +void m_assert(bool expr, std::string expr_str, std::string f, std::string file, + int line, std::string msg) +{ + std::string new_assert(f.size() + (expr ? 4 : 6), '-'); + std::cout << "\x1B[36m" << new_assert << "\033[0m\n"; + std::cout << f << ": "; + if (expr) { + std::cout << "\x1B[32mOK\033[0m\n"; + print_message(msg); + } + else { + std::cout << "\x1B[31mFAIL\033[0m\n"; + print_message(msg); + std::cout << file << " " << line << ": Assertion \"" << expr_str + << "\" Failed\n\n"; + abort(); + } +} + +bool arma_vector_close_to(arma::vec &a, arma::vec &b, double tol) +{ + if (a.n_elem != b.n_elem) { + return false; + } + + for (int i = 0; i < a.n_elem; i++) { + if (std::abs(a(i) - b(i)) >= tol) { + return false; + } + } + return true; +}