src/PenningTrap.cpp

@@ -103,6 +103,65 @@ arma::vec PenningTrap::total_force(int i)
 
 void PenningTrap::evolve_RK4(double dt)
 {
+    std::vector<Particle> 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; i<size; i++) {
+        k_v[i] = this->total_force(i)/this->particles.at(i).m;
+        k_r[i] = this->particles.at(i).v_vec;
+    }
+
+    for (int i=0; i<size; i++) {
+        Particle *p = &this->particles.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; i<size; i++) {
+        k_v[1*size + i] = this->total_force(i)/this->particles.at(i).m;
+        k_r[1*size + i] = this->particles.at(i).v_vec;
+    }
+
+    for (int i=0; i<size; i++) {
+        Particle *p = &this->particles.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; i<size; i++) {
+        k_v[2*size + i] = this->total_force(i)/this->particles.at(i).m;
+        k_r[2*size + i] = this->particles.at(i).v_vec;
+    }
+
+    for (int i=0; i<size; i++) {
+        Particle *p = &this->particles.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; i<size; i++) {
+        k_v[3*size + i] = this->total_force(i)/this->particles.at(i).m;
+        k_r[3*size + i] = this->particles.at(i).v_vec;
+    }
+
+    for (int i=0; i<size; i++) {
+        Particle *p = &this->particles.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)