#include <fstream>
#include <iostream>
#include <sstream>
#include <vector>
#include "../grapite/etics_interface.h"
#include "../grapite/particle_manager.h"

// #include <format> // C++23
//#include <print>   // C++23
#include <cstdio>   // Legacy printf

int main(int argc, char* argv[])
{
    double _double;
    int _int;

    if (argc != 2) {
        std::cout << "Specify the input file\n";
        exit(1);
    }
    std::fstream file{argv[1], std::ios::in};
    if (!file) {
        throw std::runtime_error("Could not open file");
    }

    int n_total;
    std::string line{};
    std::getline(file, line);
    std::stringstream s{line};
    s >> n_total >> _int >> _int >> _int;
    if (!s) ("Error in header");
    printf("n_total = %d\n", n_total);
    std::vector<long long> names;

    grapite::Passiv_particles particles_j{n_total};
    grapite::Active_particles particles_i{n_total}; // This is just needed to caluclate the forces, not the coefficients
    int counter{};
    printf("Loading particles...\n");
    while (std::getline(file, line)) {
        if (counter >= n_total) throw std::runtime_error("Particle exceeded maximum");
        //std::stringstream s{line};
        s = std::stringstream{line};
        int id;
        double m, x, y, z;
        long long name;
        std::string _string;
        s >> id >> m >> x >> y >> z >> _double >> _double >> _double >> _double >> _int >> name >> _int;
        if (!s) throw std::runtime_error("Error in input");
        if (id == 0) continue; // skip the black hole
        if (x == 0) continue; // bad particle
        grapite::Particle_ext p{};
        p.m = m;
        p.t = 0;
        p.pos = vec3(x, y, z);

        particles_j.add_particle(counter, &p);
        particles_i.add_particle(counter, &p);
        names.emplace_back(name);
        counter++;
        if ((n_total > 100) && (counter % (n_total/100) == 0)) {
            printf("\r%.0f%%", counter/static_cast<double>(n_total)*100);
            fflush(stdout);
        }
    }
    n_total = counter;
    printf("\rFinished\n");

    // Initialize ETICS
    grapite::Etics_interface etics{n_total};
    printf("ETICS was compiled with n_max=%d and l_max=%d\n", etics.GetNMAX(), etics.GetLMAX());
    constexpr double length_scale = 10.0;
    etics::scf::SetLengthScale(length_scale);
    printf("ETICS scaling is %.3e\n", length_scale);

    // Calculate the expanstion from particles_j
    particles_j.nj_total = n_total;
    particles_j.nj_core  = 0;
    particles_j.flush_memory();
    etics.calculate_expansion(&particles_j, grapite::Selector::halo, 0);

    // Write coefficients to a file
    std::ofstream coeff_file("out.coeff", std::ofstream::binary);
    coeff_file.write((char*)etics.A_h, 1024*sizeof(Complex));
    coeff_file.close();

    // Calculate potentials and forces
    particles_i.ni_core  = 0;
    particles_i.ni_total = n_total;
    particles_i.flush_memory();
    std::vector<double> pot(n_total);
    std::vector<vec3> acc(n_total);
    etics.calculate_gravity(&particles_i, grapite::Selector::halo, 0, pot.data(), acc.data());

    // Write results to a file
    printf("Writing results...\n");
    std::ofstream output_file("out.dat");
    for (int i = 0; i < counter; i++) {
        constexpr int buf_size = 4096;
        char line[buf_size];
        auto [ax, ay, az] = acc[i];
        snprintf(line, buf_size, "%015lld%16.8e%16.8e%16.8e%16.8e\n", names[i], pot[i], ax, ay, az);
        output_file.write(line, strlen(line));
        if ((counter > 100) && (i % (counter/100) == 0)) {
            printf("\r%.0f%%", i/static_cast<double>(counter)*100);
            fflush(stdout);
        }
    }
    printf("\rFinished\n");

}