/** @file main.cpp
 *
 *  @author Cory Alexander Balaton (coryab)
 *  @author Janita Ovidie Sandtrøen Willumsen (janitaws)
 *
 *  @version 0.1
 *
 *  @brief The main program for this project
 *
 *  @bug No known bugs
 * */

#include <fstream>
#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;
}