From dfb287e0b2e00ea36ff79956097a8ad51ad634f0 Mon Sep 17 00:00:00 2001
From: Cory <cory@balaton.dev>
Date: Mon, 2 Oct 2023 21:50:53 +0200
Subject: [PATCH] Simulate 100 particles using forward euler

---
 src/main.cpp | 74 +++++++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 73 insertions(+), 1 deletion(-)

diff --git a/src/main.cpp b/src/main.cpp
index 4c46902..b479649 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -10,7 +10,79 @@
  *  @bug No known bugs
  * */
 
-int main() 
+#include <filesystem>
+#include <fstream>
+#include <string>
+#include <omp.h>
+#include <sys/stat.h>
+
+#include "PenningTrap.hpp"
+
+#define PARTICLES 100
+#define N 10000
+#define CHARGE 1.
+#define MASS 40. // unit: amu
+
+void euler_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
+        trap.add_particle(Particle(CHARGE, MASS, r, v));
+    }
+
+    double time = 50.; // microseconds
+    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++) {
+            res[i][j] = trap.get_particle(i);
+        }
+        trap.evolve_forward_euler(dt);
+    }
+
+    std::cout << counter << std::endl;
+
+    arma::vec::fixed<3>* cur_row;
+    arma::vec::fixed<3> cur_elem;
+
+    mkdir("output", 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++) {
+        cur_row = res[i];
+        ofile.open("output/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.close();
+    }
+}
+
+int main()
+{
+    double start = omp_get_wtime();
+
+    euler_100_particles();
+
+    double end = omp_get_wtime();
+
+    std::cout << "Time: " << end - start << " seconds" << std::endl;
+
     return 0;
 }