From 4a0ee772761b7c96acc6ff4b330c01eaeddeaf1a Mon Sep 17 00:00:00 2001
From: Yohai Meiron <sahmes@gmail.com>
Date: Thu, 9 Apr 2020 20:36:47 -0400
Subject: [PATCH] first commit

---
 .gitignore |   3 ++
 main.cpp   | 121 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 plot.py    |  24 +++++++++++
 3 files changed, 148 insertions(+)
 create mode 100644 .gitignore
 create mode 100644 main.cpp
 create mode 100644 plot.py

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..6f8fe4b
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,3 @@
+.*
+main
+libmain.so
diff --git a/main.cpp b/main.cpp
new file mode 100644
index 0000000..fef566d
--- /dev/null
+++ b/main.cpp
@@ -0,0 +1,121 @@
+#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;
+}
\ No newline at end of file
diff --git a/plot.py b/plot.py
new file mode 100644
index 0000000..d2df0db
--- /dev/null
+++ b/plot.py
@@ -0,0 +1,24 @@
+import numpy as np
+import ctypes, subprocess
+
+recompile   = True
+
+if recompile:
+    p = subprocess.Popen('g++ -shared -o libmain.so -fPIC main.cpp -lgsl'.split())
+    p.wait()
+    if p.returncode != 0: raise RuntimeError(p.returncode)
+
+libmain = ctypes.CDLL('./libmain.so')
+def integrate(y0, t_max):
+    y0 = (ctypes.c_double*6)(*y0) 
+    y  = (ctypes.c_double*6)()
+    status = libmain.integrate(y0, ctypes.c_double(t_max), y)
+    return np.array(y), status
+
+gsl_success = libmain.gsl_success()
+
+#t_array = linspace(plot_tmin, plot_tmax, plot_points)
+#r_array = empty_like(t_array)
+
+res = integrate([10,0,0,0,200,0], 0.1)
+print(res, gsl_success)
\ No newline at end of file