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

Coryab/create tests #4

Merged
coryab merged 10 commits from coryab/create-tests into develop 2023-10-07 20:01:08 +00:00
Conversation 0 Commits 10 Files Changed 8 +17 -24
1 changed files with 17 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files
Showing only changes of commit 8a2100a334 - Show all commits

41
src/PenningTrap.cpp
View File

@ -32,20 +32,17 @@ void PenningTrap::add_particle(Particle particle)
arma::vec PenningTrap::external_E_field(arma::vec r) arma::vec PenningTrap::external_E_field(arma::vec r)
{ {
arma::vec::fixed<3> res; 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(0) = r(0);
res(1) = r(1); 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 PenningTrap::external_B_field(arma::vec r)
{ {
arma::vec::fixed<3> res; arma::vec::fixed<3> res{0., 0., this->B_0};
res(0) = 0.;
res(1) = 0.;
res(2) = this->B_0;
return res; return res;
} }
@ -53,14 +50,14 @@ arma::vec PenningTrap::external_B_field(arma::vec r)
arma::vec PenningTrap::force_on_particle(int i, int j) arma::vec PenningTrap::force_on_particle(int i, int j)
{ {
// Calculate the difference between the particles' position // Calculate the difference between the particles' position
arma::vec::fixed<3> res = this->particles.at(i).r_vec arma::vec::fixed<3> res =
- this->particles.at(j).r_vec; this->particles.at(i).r_vec - this->particles.at(j).r_vec;
// Get the distance between the particles // Get the distance between the particles
double norm = arma::norm(res); double norm = arma::norm(res);
// Multiply res with p_j's charge divided by the norm cubed // 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; return res;
} }
@ -69,16 +66,13 @@ arma::vec PenningTrap::total_force_external(int i)
{ {
Particle p = this->particles.at(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); 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); arma::vec::fixed<3> v_cross_B{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); 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); p.v_vec(0) * B(1) - p.v_vec(1) * B(0)};
arma::vec force = p.q arma::vec force = p.q * (this->external_E_field(p.r_vec) + v_cross_B);
*(this->external_E_field(p.r_vec) + v_cross_B);
return force; return force;
} }
@ -89,7 +83,7 @@ arma::vec PenningTrap::total_force_particles(int i)
arma::vec res(3); 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) { if (i == j) {
continue; continue;
} }
@ -97,7 +91,7 @@ arma::vec PenningTrap::total_force_particles(int i)
res += this->force_on_particle(i, j); res += this->force_on_particle(i, j);
} }
res *= K_E*(p.q/p.m); res *= K_E * (p.q / p.m);
return res; return res;
} }
@ -109,7 +103,6 @@ arma::vec PenningTrap::total_force(int i)
void PenningTrap::evolve_RK4(double dt) void PenningTrap::evolve_RK4(double dt)
{ {
} }
void PenningTrap::evolve_forward_euler(double dt) void PenningTrap::evolve_forward_euler(double dt)
@ -118,11 +111,11 @@ void PenningTrap::evolve_forward_euler(double dt)
Particle *p; Particle *p;
#pragma omp parallel for private(p) #pragma omp parallel for private(p)
for (int i=0; i < this->particles.size(); i++) { for (int i = 0; i < this->particles.size(); i++) {
p = &new_state.at(i); p = &new_state.at(i);
p->v_vec += dt*this->total_force(i)/new_state.at(i).m; p->v_vec += dt * this->total_force(i) / new_state.at(i).m;
p->r_vec += dt*this->particles.at(i).v_vec; p->r_vec += dt * this->particles.at(i).v_vec;
} }
this->particles = new_state; this->particles = new_state;