Penning Trap Simulation Simulate particle behavior inside a Penning Trap

Loading...
Searching...
No Matches
main.cpp
Go to the documentation of this file.
1
13#include <cmath>
14#include <complex>
15#include <fstream>
16#include <omp.h>
17#include <string>
18#include <vector>
19
20#include "PenningTrap.hpp"
21#include "utils.hpp"
22
23#define PARTICLES 100
24#define N 40000
25#define CHARGE 1. // unit: e
26#define MASS 40. // unit: amu
27
28// Particles used for testing
29 Particle p1(CHARGE, MASS, vec_3d{20., 0., 20.}, vec_3d{0., 25., 0.});
30 Particle p2(CHARGE, MASS, vec_3d{25., 25., 0.}, vec_3d{0., 40., 5.});
31
38 vec_3d analytical_solution_particle_1(double t)
39{
40 double w_0 = T / MASS;
41 double w_z2 = (50. * V / 1000.) / (MASS * 500. * 500.);
42 double w_p = (w_0 + std::sqrt(w_0 * w_0 - 2. * w_z2)) / 2.;
43 double w_n = (w_0 - std::sqrt(w_0 * w_0 - 2. * w_z2)) / 2.;
44 double A_p = (25. + w_n * 20.) / (w_n - w_p);
45 double A_n = -(25. + w_p * 20.) / (w_n - w_p);
46 std::complex<double> f =
47 A_p * std::exp(std::complex<double>(0., -w_p * t)) +
48 A_n * std::exp(std::complex<double>(0., -w_n * t));
49 vec_3d res{std::real(f), std::imag(f), 20. * std::cos(std::sqrt(w_z2) * t)};
50 return res;
51}
52
55 void simulate_single_particle()
56{
57 // Initialize trap with particle 1
58 PenningTrap trap(std::vector<Particle>{p1});
59
60 double time = 50.; // microseconds
61
62 // Simulate and write results to file
63 trap.write_simulation_to_dir("output/simulate_single_particle", time, N,
64 "rk4", false);
65}
66
70 void simulate_two_particles()
71{
72 // Initialize traps with particles
73 PenningTrap trap_no_interaction(std::vector<Particle>{p1, p2});
74 PenningTrap trap_with_interaction(std::vector<Particle>{p1, p2});
75
76 double time = 50.; // microseconds
77
78 // Simulate and write results to files
79 trap_no_interaction.write_simulation_to_dir(
80 "output/simulate_2_particles/no_interaction", time, N, "rk4", false);
81 trap_with_interaction.write_simulation_to_dir(
82 "output/simulate_2_particles/with_interaction", time, N);
83}
84
88 void simulate_single_particle_with_different_steps()
89{
90 double time = 50.; // microseconds
91
92 std::ofstream ofile;
93
94 // Calculate relative error for RK4
95 std::string path = "output/relative_error/RK4/";
96 mkpath(path);
97 for (int i = 0; i < 4; i++) {
98 int steps = 4000 * std::pow(2, i);
99 double dt = time / (double)steps;
100 ofile.open(path + std::to_string(steps) + "_steps.txt");
101 PenningTrap trap(std::vector<Particle>{p1});
102 simulation_t res = trap.simulate(time, steps, "rk4", false);
103 for (int i = 0; i < steps; i++) {
104 ofile << arma::norm(res.r_vecs[0][i] -
105 analytical_solution_particle_1(dt*i))
106 << "\n";
107 }
108 ofile.close();
109 }
110
111 // Calculate relative error for forward Euler
112 path = "output/relative_error/euler/";
113 mkpath(path);
114 for (int i = 0; i < 4; i++) {
115 int steps = 4000 * std::pow(2, i);
116 double dt = time / (double)steps;
117 ofile.open(path + std::to_string(steps) + "_steps.txt");
118 PenningTrap trap(std::vector<Particle>{p1});
119 simulation_t res = trap.simulate(time, steps, "euler", false);
120 for (int i = 0; i < steps; i++) {
121 ofile << arma::norm(res.r_vecs[0][i] -
122 analytical_solution_particle_1(dt*i))
123 << "\n";
124 }
125 ofile.close();
126 }
127}
128
131 void simulate_100_particles()
132{
133 PenningTrap trap((unsigned)100);
134
135 double time = 50.; // microseconds
136
137 trap.write_simulation_to_dir("output/simulate_100_particles", time, N);
138}
139
146 void simulate_100_particles_with_time_potential()
147{
148 double time = 500.;
149
150 double amplitudes[]{.1, .4, .7};
151
152 double freq_start = .2;
153 double freq_end = 2.5;
154 double freq_increment = .02;
155 size_t freq_iterations = (size_t)((freq_end - freq_start) / freq_increment);
156
157 double res[4][freq_iterations];
158
159 std::string path = "output/time_dependent_potential/";
160 mkpath(path);
161
162 std::ofstream ofile;
163
164 double freq = freq_start;
165 for (size_t i = 0; i < freq_iterations; i++) {
166 res[0][i] = freq;
167 freq += freq_increment;
168 }
169
170#pragma omp parallel for collapse(2) num_threads(4)
171 for (size_t i = 0; i < 3; i++) {
172 for (size_t j = 0; j < freq_iterations; j++) {
173 PenningTrap trap(
174 (unsigned)100, T,
175 std::bind(
176 [](double f, double r, double t) {
177 return (25. * V / 1000.) * (1. + f * std::cos(r * t));
178 },
179 amplitudes[i], res[0][j], std::placeholders::_1),
180 500., 0.);
181 res[i + 1][j] =
182 trap.fraction_of_particles_left(time, N, "rk4", false);
183 }
184 }
185
186 ofile.open(path + "res.txt");
187 for (size_t i = 0; i < freq_iterations; i++) {
188 ofile << res[0][i] << "," << res[1][i] << "," << res[2][i] << ","
189 << res[3][i] << "\n";
190 }
191 ofile.close();
192}
193
194 int main()
195{
196 double t0 = omp_get_wtime();
197
198 // simulate_single_particle();
199
200 // simulate_two_particles();
201
202 simulate_single_particle_with_different_steps();
203
204 double t1 = omp_get_wtime();
205
206 simulate_100_particles();
207
208 //simulate_100_particles_with_time_potential();
209
210 double end = omp_get_wtime();
211
212 std::cout << "Time: " << (end - t1) << " seconds" << std::endl;
213
214 return 0;
215}
PenningTrap.hpp
A class for simulating a Penning trap.
Particle
A class that holds attributes of a particle.
Definition: Particle.hpp:21
PenningTrap
A class that simulates a Penning trap.
Definition: PenningTrap.hpp:31
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::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
T
#define T
1 Tesla. unit:
Definition: constants.hpp:17
V
#define V
1 Volt. unit:
Definition: constants.hpp:19
simulate_100_particles
void simulate_100_particles()
Simulate 100 particles over 50 .
Definition: main.cpp:131
simulate_100_particles_with_time_potential
void simulate_100_particles_with_time_potential()
Simulate 100 particles over 500 using a time dependent potential.
Definition: main.cpp:146
analytical_solution_particle_1
vec_3d analytical_solution_particle_1(double t)
The analytical solution for particle p1.
Definition: main.cpp:38
simulate_two_particles
void simulate_two_particles()
Simulate 2 particles over the period of 50 with and without particle interactions.
Definition: main.cpp:70
p1
Particle p1(CHARGE, MASS, vec_3d{20., 0., 20.}, vec_3d{0., 25., 0.})
Particle 1.
simulate_single_particle
void simulate_single_particle()
Simulate a single particle over the period of 50 .
Definition: main.cpp:55
p2
Particle p2(CHARGE, MASS, vec_3d{25., 25., 0.}, vec_3d{0., 40., 5.})
Particle 2.
simulate_single_particle_with_different_steps
void simulate_single_particle_with_different_steps()
Simulate a single particle over 50 using different amount of steps and different methods.
Definition: main.cpp:88
simulation
Typedef for PenningTrap::simulation return value.
Definition: typedefs.hpp:40
vec_3d
arma::vec::fixed< 3 > vec_3d
Typedef for a fixed 3d arma vector.
Definition: typedefs.hpp:36
utils.hpp
Function prototypes and macros that are useful.
mkpath
bool mkpath(std::string path, int mode=0777)
Make path given.
Definition: utils.cpp:74