Penning Trap Simulation Simulate particle behavior inside a Penning Trap

Loading...
Searching...
No Matches
PenningTrap.cpp
Go to the documentation of this file.
1
13#include "PenningTrap.hpp"
14
15 PenningTrap::PenningTrap(double B_0, std::function<double(double)> V_0,
16 double d, double t)
17{
18 this->B_0 = B_0;
19 this->V_0 = V_0;
20 this->d = d;
21 this->t = t;
22}
23
24 PenningTrap::PenningTrap(unsigned int i, double B_0,
25 std::function<double(double)> V_0, double d, double t)
26 : PenningTrap::PenningTrap(B_0, V_0, d)
27{
28 vec_3d r, v;
29 for (size_t j = 0; j < i; j++) {
30 r = vec_3d().randn() * .1 * this->d;
31 v = vec_3d().randn() * .1 * this->d;
32 this->add_particle(Particle(1., 40., r, v));
33 }
34}
35
36 PenningTrap::PenningTrap(std::vector<Particle> particles, double B_0,
37 std::function<double(double)> V_0, double d, double t)
38 : PenningTrap::PenningTrap(B_0, V_0, d)
39{
40 this->particles = particles;
41}
42
43 vec_3d PenningTrap::v_func(unsigned int i, unsigned int j, double dt)
44{
45 switch (i) {
46 case 0:
47 return .5 * dt * this->k_v[0][j];
48 case 1:
49 return .5 * dt * this->k_v[1][j];
50 case 2:
51 return dt * this->k_v[2][j];
52 case 3:
53 return (dt / 6.) * (this->k_v[0][j] + 2. * this->k_v[1][j] +
54 2. * this->k_v[2][j] + this->k_v[3][j]);
55 default:
56 std::cout << "Not valid!" << std::endl;
57 abort();
58 }
59}
60
61 vec_3d PenningTrap::r_func(unsigned int i, unsigned int j, double dt)
62{
63 switch (i) {
64 case 0:
65 return .5 * dt * this->k_r[0][j];
66 case 1:
67 return .5 * dt * this->k_r[1][j];
68 case 2:
69 return dt * this->k_r[2][j];
70 case 3:
71 return (dt / 6.) * (this->k_r[0][j] + 2. * this->k_r[1][j] +
72 2. * this->k_r[2][j] + this->k_r[3][j]);
73 default:
74 std::cout << "Not valid!" << std::endl;
75 abort();
76 }
77}
78
79 void PenningTrap::add_particle(Particle particle)
80{
81 this->particles.push_back(particle);
82}
83
84 vec_3d PenningTrap::external_E_field(vec_3d r)
85{
86 r(2) *= -2.;
87 double f = this->V_0(this->t) / (this->d * this->d);
88
89 return f * r;
90}
91
92 vec_3d PenningTrap::external_B_field(vec_3d r)
93{
94 return vec_3d{0., 0., this->B_0};
95}
96
97 vec_3d PenningTrap::force_on_particle(unsigned int i, unsigned int j)
98{
99 // Calculate the difference between the particles' position
100 vec_3d res = this->particles[i].r_vec - this->particles[j].r_vec;
101
102 // Get the distance between the particles
103 double norm = arma::norm(res, 2);
104
105 return vec_3d(res * this->particles[j].q / (norm * norm * norm));
106}
107
108 vec_3d PenningTrap::total_force_external(unsigned int i)
109{
110 Particle p = this->particles[i];
111
112 if (arma::norm(p.r_vec) > this->d) {
113 return vec_3d{0., 0., 0.};
114 }
115
116 vec_3d force =
117 p.q * (this->external_E_field(p.r_vec) +
118 arma::cross(p.v_vec, this->external_B_field(p.r_vec)));
119
120 return force;
121}
122
123 vec_3d PenningTrap::total_force_particles(unsigned int i)
124{
125 Particle p = this->particles[i];
126
127 vec_3d res;
128
129 for (size_t j = 0; j < this->particles.size(); j++) {
130 if (i == j) {
131 continue;
132 }
133
134 res += this->force_on_particle(i, j);
135 }
136
137 return vec_3d(res * K_E * (p.q / p.m));
138}
139
140 vec_3d PenningTrap::total_force(unsigned int i)
141{
142 if (arma::norm(this->particles[i].r_vec) > this->d) {
143 return vec_3d{0., 0., 0.};
144 }
145 return this->total_force_external(i) - this->total_force_particles(i);
146}
147
148 void PenningTrap::evolve_RK4(double dt, bool particle_interaction)
149{
150
151 std::vector<Particle> original_particles = this->particles;
152 std::vector<Particle> tmp_particles = this->particles;
153
154 vec_3d (PenningTrap::*force)(unsigned int);
155 if (particle_interaction) {
156 force = &PenningTrap::total_force;
157 }
158 else {
159 force = &PenningTrap::total_force_external;
160 }
161
162 size_t size = this->particles.size();
163
164 this->k_v = sim_arr(4, sim_cols(size));
165 this->k_r = sim_arr(4, sim_cols(size));
166
167 for (size_t i = 0; i < 4; i++) {
168#pragma omp parallel for
169 for (size_t j = 0; j < this->particles.size(); j++) {
170 this->k_v[i][j] = (this->*force)(j) / this->particles[j].m;
171 this->k_r[i][j] = this->particles[j].v_vec;
172
173 Particle *p = &tmp_particles[j];
174
175 p->v_vec = original_particles[j].v_vec + this->v_func(i, j, dt);
176 p->r_vec = original_particles[j].r_vec + this->r_func(i, j, dt);
177 }
178 this->particles.swap(tmp_particles);
179 }
180 this->t += dt;
181}
182
183 void PenningTrap::evolve_forward_euler(double dt, bool particle_interaction)
184{
185 size_t size = this->particles.size();
186 vec_3d force_res[size];
187 Particle *p;
188
189 vec_3d (PenningTrap::*force)(unsigned int);
190 if (particle_interaction) {
191 force = &PenningTrap::total_force;
192 }
193 else {
194 force = &PenningTrap::total_force_external;
195 }
196
197#pragma omp parallel for
198 for (size_t i = 0; i < size; i++) {
199 force_res[i] = (this->*force)(i);
200 }
201
202#pragma omp parallel for private(p)
203 for (size_t i = 0; i < size; i++) {
204 p = &this->particles[i];
205 p->r_vec += dt * p->v_vec;
206 p->v_vec += dt * force_res[i] / p->m;
207 }
208
209 this->t += dt;
210}
211
212 simulation_t PenningTrap::simulate(double time, unsigned int steps,
213 std::string method,
214 bool particle_interaction)
215{
216 double dt = time / (double)steps;
217
218 unsigned int size = this->particles.size();
219
220 simulation_t res{sim_arr(size, sim_cols(steps)),
221 sim_arr(size, sim_cols(steps))};
222
223 void (PenningTrap::*func)(double, bool);
224 if (method == "rk4") {
225 func = &PenningTrap::evolve_RK4;
226 }
227 else if (method == "euler") {
228 func = &PenningTrap::evolve_forward_euler;
229 }
230 else {
231 std::cout << "Not a valid method!" << std::endl;
232 abort();
233 }
234
235 for (size_t j = 0; j < steps; j++) {
236 for (size_t i = 0; i < size; i++) {
237 res.r_vecs[i][j] = this->particles[i].r_vec;
238 res.v_vecs[i][j] = this->particles[i].v_vec;
239 }
240 (this->*func)(dt, particle_interaction);
241 }
242
243 return res;
244}
245
246 void PenningTrap::write_simulation_to_dir(std::string path, double time,
247 unsigned int steps,
248 std::string method,
249 bool particle_interaction)
250{
251 if (path.back() != '/') {
252 path += '/';
253 }
254 if (mkpath(path, 0777) != 0) {
255 std::cout << "Hello" << std::endl;
256 return;
257 }
258
259 simulation_t res =
260 this->simulate(time, steps, method, particle_interaction);
261
262 std::ofstream ofile;
263
264#pragma omp parallel for private(ofile)
265 for (size_t i = 0; i < this->particles.size(); i++) {
266 ofile.open(path + "particle_" + std::to_string(i) + "_r.txt");
267 for (vec_3d &vec : res.r_vecs[i]) {
268 ofile << vec(0) << "," << vec(1) << "," << vec(2) << "\n";
269 }
270 ofile.close();
271 ofile.open(path + "particle_" + std::to_string(i) + "_v.txt");
272 for (vec_3d &vec : res.v_vecs[i]) {
273 ofile << scientific_format(vec(0), 10, 8) << ","
274 << scientific_format(vec(1), 8, 10) << ","
275 << scientific_format(vec(2), 8, 10) << "\n";
276 }
277 ofile.close();
278 }
279}
280
281 double PenningTrap::fraction_of_particles_left(double time, unsigned int steps,
282 std::string method,
283 bool particle_interaction)
284{
285 simulation_t res =
286 this->simulate(time, steps, method, particle_interaction);
287
288 int particles_left = 0;
289
290 for (Particle p : this->particles) {
291 if (arma::norm(p.r_vec) < this->d) {
292 particles_left++;
293 }
294 }
295
296 return (double)particles_left / (double)this->particles.size();
297}
298
299 vec_3d PenningTrap::get_r(int i)
300{
301 return this->particles[i].r_vec;
302}
303
304 double PenningTrap::get_t()
305{
306 return this->t;
307}
PenningTrap.hpp
A class for simulating a Penning trap.
Particle
A class that holds attributes of a particle.
Definition: Particle.hpp:21
Particle::q
double q
Charge.
Definition: Particle.hpp:23
Particle::v_vec
vec_3d v_vec
velocity
Definition: Particle.hpp:26
Particle::m
double m
Mass.
Definition: Particle.hpp:24
Particle::r_vec
vec_3d r_vec
position
Definition: Particle.hpp:25
PenningTrap
A class that simulates a Penning trap.
Definition: PenningTrap.hpp:31
PenningTrap::particles
std::vector< Particle > particles
The particles in the Penning trap.
Definition: PenningTrap.hpp:37
PenningTrap::fraction_of_particles_left
double fraction_of_particles_left(double time, unsigned int steps, std::string method="rk4", bool particle_interaction=true)
Simulate and calculate what fraction of particles are still left inside the Penning trap after the si...
Definition: PenningTrap.cpp:281
PenningTrap::B_0
double B_0
Magnetic field strength.
Definition: PenningTrap.hpp:33
PenningTrap::external_E_field
vec_3d external_E_field(vec_3d r)
Calculate E at point r.
Definition: PenningTrap.cpp:84
PenningTrap::total_force
vec_3d total_force(unsigned int i)
calculate the total force on a particle p_i.
Definition: PenningTrap.cpp:140
PenningTrap::k_r
sim_arr k_r
Definition: PenningTrap.hpp:40
PenningTrap::evolve_RK4
void evolve_RK4(double dt, bool particle_interaction=true)
Go forward one timestep using the RK4 method.
Definition: PenningTrap.cpp:148
PenningTrap::r_func
vec_3d r_func(unsigned int i, unsigned int j, double dt)
Helper for evolve_RK4 when calculating values.
Definition: PenningTrap.cpp:61
PenningTrap::total_force_particles
vec_3d total_force_particles(unsigned int i)
Calculate the total force on a particle p_i from other particles.
Definition: PenningTrap.cpp:123
PenningTrap::d
double d
Characteristic dimension.
Definition: PenningTrap.hpp:35
PenningTrap::add_particle
void add_particle(Particle particle)
Add a particle to the system.
Definition: PenningTrap.cpp:79
PenningTrap::t
double t
Current time.
Definition: PenningTrap.hpp:36
PenningTrap::force_on_particle
vec_3d force_on_particle(unsigned int i, unsigned int j)
Calculate the force between 2 particles.
Definition: PenningTrap.cpp:97
PenningTrap::V_0
std::function< double(double)> V_0
Applied potential.
Definition: PenningTrap.hpp:34
PenningTrap::evolve_forward_euler
void evolve_forward_euler(double dt, bool particle_interaction=true)
Go forward one timestep using the forward Euler method.
Definition: PenningTrap.cpp:183
PenningTrap::write_simulation_to_dir
void write_simulation_to_dir(std::string path, double time, unsigned int steps, std::string method="rk4", bool particle_interaction=true)
Simulate and write the displacement of all particles to files.
Definition: PenningTrap.cpp:246
PenningTrap::external_B_field
vec_3d external_B_field(vec_3d r)
Calculate B at point r.
Definition: PenningTrap.cpp:92
PenningTrap::simulate
simulation_t simulate(double time, unsigned int steps, std::string method="rk4", bool particle_interaction=true)
Simulate the particle system inside the Penning trap over a certain amount of time.
Definition: PenningTrap.cpp:212
PenningTrap::PenningTrap
PenningTrap(double B_0=T, std::function< double(double)> V_0=[](double t) { return 25. *V/1000.;}, double d=500., double t=0.)
Constructor for the PenningTrap class.
Definition: PenningTrap.cpp:15
PenningTrap::total_force_external
vec_3d total_force_external(unsigned int i)
Calculate the total external force on a particle.
Definition: PenningTrap.cpp:108
PenningTrap::k_v
sim_arr k_v
Definition: PenningTrap.hpp:38
PenningTrap::v_func
vec_3d v_func(unsigned int i, unsigned int j, double dt)
Helper for evolve_RK4 when calculating values.
Definition: PenningTrap.cpp:43
K_E
#define K_E
Coulomb constant. unit: .
Definition: constants.hpp:15
simulation
Typedef for PenningTrap::simulation return value.
Definition: typedefs.hpp:40
sim_cols
std::vector< arma::vec::fixed< 3 > > sim_cols
Typedef for the column of the result vector from simulating particles.
Definition: typedefs.hpp:24
vec_3d
arma::vec::fixed< 3 > vec_3d
Typedef for a fixed 3d arma vector.
Definition: typedefs.hpp:36
sim_arr
std::vector< sim_cols > sim_arr
Typedef for the result of the simulate method.
Definition: typedefs.hpp:32
mkpath
bool mkpath(std::string path, int mode=0777)
Make path given.
Definition: utils.cpp:74
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:15