#include <iostream>
#include <string>
#include <vector>
#include <fstream>
#include <numeric>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_spline.h>


extern "C" const int gsl_success() { return GSL_SUCCESS; } // It's zero, but just for clarity sake.

// Our units are {kiloparsec, solar mass, gigayear}
constexpr double G = 4.498317481097514e-06;

class Interp {
public:
    Interp(std::vector<double>& x, std::vector<double>& y)
    {
        acc    = gsl_interp_accel_alloc();
        spline = gsl_spline_alloc(gsl_interp_cspline, x.size());
        gsl_spline_init(spline, x.data(), y.data(), x.size()); 
    }
    inline double operator()(double x) const
    {
        return gsl_spline_eval(spline, x, acc);
    }

private:
    gsl_interp_accel *acc;
    gsl_spline *spline;
};

class Plummer {
public:
    Plummer(double M, double b)
        : M(M), b(b) {}
    void calc_acceleration(const double *pos, double *acc)
    {
        double r2 = (pos[0]*pos[0] + pos[1]*pos[1] + pos[2]*pos[2] + b*b);
        double r  = sqrt(r2);
        double r3_inv = 1/(r*r2);
        acc[0] = G*M*pos[0]*r3_inv;
        acc[1] = G*M*pos[1]*r3_inv;
        acc[2] = G*M*pos[2]*r3_inv;
    }

private:
    double M, b;
};

class Galaxy {
public:
    Galaxy(std::string file_name)
    {
        std::vector<double> t, M_bulge_data, b_bulge_data;

        t.push_back(0);
        t.push_back(1);
        t.push_back(2);

        M_bulge_data.push_back(9e9);
        M_bulge_data.push_back(1e10);
        M_bulge_data.push_back(2e10);

        interp_M_bulge = new Interp(t, M_bulge_data);
        std::cout << "exiting..." << std::endl;
        exit(0);
    }
    int func(double t, const double y[], double f[], void *params)
    {
        double M_bulge = (*interp_M_bulge)(t);
        Plummer plummer(M_bulge, 0);
        f[0] = y[3]; // vx -> x'
        f[1] = y[4]; // vy -> y'
        f[2] = y[5]; // vz -> z'
        plummer.calc_acceleration(y, f+3); // a -> v'
        return GSL_SUCCESS;
    }

private:
    Interp *interp_M_bulge;
    Interp *interp_b_bulge;
} galaxy("file.dat");
// Not very nice to have it as a global variable but GSL will have problem otherwise.

int jac(double t, const double y[], double *dfdy, double dfdt[], void *params)
{
    return GSL_SUCCESS;
}


int func(double t, const double y[], double f[], void *params)
{
    return galaxy.func(t, y, f, params);
}

extern "C"
int integrate(const double y0[], const double t_max, double y[])
{
    double t = 0;
    constexpr double h = 1./4096.;
    constexpr double epsabs = 1e-7;
    constexpr double epsrel = 0;
    const gsl_odeiv2_step_type *T = gsl_odeiv2_step_rk8pd;
    gsl_odeiv2_step *s    = gsl_odeiv2_step_alloc(T, 6);
    gsl_odeiv2_evolve *e  = gsl_odeiv2_evolve_alloc(6);
    gsl_odeiv2_control *c = gsl_odeiv2_control_y_new(epsabs, 0);
    gsl_odeiv2_system sys = {func, jac, 6, nullptr};
    gsl_odeiv2_driver *d  = gsl_odeiv2_driver_alloc_y_new(&sys, T, h, epsabs, epsrel);

    std::copy(y0, y0+6, y);
    int status = gsl_odeiv2_driver_apply(d, &t, t_max, y);
    return status;
}

int main()
{
    std::cout << "bye" << std::endl;
}