/* Header-Datei
 * Das n-Körperproblem in drei Dimensionen
 * n gravitativ wirkende Körper bewegen sich durch ihre gegenseitige Newtonsche Anziehungskraft
 * Berechnung der Lösung eines Systems von 2*3*n Differentialgleichungen erster Ordnung
 * mittels Runge-Kutta Ordnung vier Methode
 * Verfahren zur Loesung der DGL ist in einer abstrakten Basisklasse ausgelagert
 * Berechnung der Lösung innerhalb einer Methode (Klassenfunktion) Runge-Kutta-Verfahren
 * Überschreibung der virtuellen Funktion der Bewegungsgleichungen in der abgeleiteten Klasse (Subklasse) 'nBody'
 * Zeitentwicklung der fuer unterschiedliche t-Werte in [a,b]
 * Unterklasse 'nBody'
 * Konstruktor 1: 'nBody'(Grav.Konstante G, Massenvektor der Körper, Anfangszeit a, Endzeit b, Anzahl der Punkte N, Anfangswerte u_init)
 * Konstruktor 2: Elliptische Lagrangesche homografischen Lösungen (Exzentrizität exz in [0,1])
 *                'nBody'(Anfangszeit a, Endzeit b, Anzahl der Punkte N, Exzentrizität exz)
 * Ausgabe zum Plotten in Datei
 */

#include <cmath>               // Bibliothek für mathematisches (e-Funktion, Betrag, ...)
#include <vector>              // Vector-Container der Standardbibliothek
using namespace std;           // Benutze den Namensraum std

// Abstrakte Basisklasse 'dsolve'
class dsolve {
public:
    // Daten-Member der Klasse
    double a;                        // Untergrenze des Zeit-Intervalls
    double b;                        // Obergrenze des Intervalls
    int N;                           // Anzahl der Punkte
    vector<double> alpha;            // Anfangswerte

    vector<vector<double>> ym;       // Lösungsmatrix
    vector<double> t;                // Zeit-Vektor

    // Virtuelle Funktion dgls, wird an anderer Stelle definiert
    virtual vector<double> dgls(double t, vector<double> u_vec) = 0;

    // Konstruktor mit vier Argumenten (Initialisierung der Daten-Member)
    dsolve(double a_, double b_, int N_, vector<double> alpha_) :a(a_), b(b_), N(N_), alpha(alpha_) {
        t.push_back(a_);                // Zum Zeit-Vektor die Anfangszeit eintragen
        ym.push_back(alpha_);           // Zur Loesungs-Matrix die Anfangswerte eintragen
    }

    // Methode der Runge-Kutta-Methode
    void solve() {
        double h = (b - a)/N;                   // Abstand dt zwischen den aequidistanten Punkten des t-Intervalls (h=dt)
        int n = alpha.size();                   // Anzahl der Anfangswerte bzw. DGL-Gleichungen (hier n=2)
        vector<double> k1(n),k2(n),k3(n),k4(n); // Deklaration der vier Runge-Kutta Parameter fuer jede Loesung
        vector<double> y(n);                    // Temporaerer Hilfsvektor

        for (int i = 0; i < N; ++i) {                                         // for-Schleife ueber die einzelnen Punkte des t-Intervalls
            for (int j = 0; j < n; ++j) k1[j] = h * dgls(t[i], ym[i])[j];     // for-Schleife ueber die n Runge-Kutta k1-Parameter
            for (int j = 0; j < n; ++j) y[j] = ym[i][j] + k1[j] / 2;          // Temporaerer Hilfsvektor
            for (int j = 0; j < n; ++j) k2[j] = h * dgls(t[i] + h / 2, y)[j]; // for-Schleife ueber die n Runge-Kutta k2-Parameter
            for (int j = 0; j < n; ++j) y[j] = ym[i][j] + k2[j] / 2;          // Temporaerer Hilfsvektor
            for (int j = 0; j < n; ++j) k3[j] = h * dgls(t[i] + h / 2, y)[j]; // for-Schleife ueber die n Runge-Kutta k3-Parameter
            for (int j = 0; j < n; ++j) y[j] = ym[i][j] + k3[j];              // Temporaerer Hilfsvektor
            for (int j = 0; j < n; ++j) k4[j] = h * dgls(t[i] + h, y)[j];     // for-Schleife ueber die n Runge-Kutta k4-Parameter
            for (int j = 0; j < n; ++j) y[j] = ym[i][j] + (k1[j] + 2 * k2[j] + 2 * k3[j] + k4[j]) / 6; // Temporaerer Hilfsvektor

            ym.push_back(y);              // Zum Loesungs-Vektor den neuen Wert eintragen
            t.push_back(a + (i + 1) * h); // Zum Zeit-Vektor die neue Zeit eintragen
        }                                 // Ende for-Schleife ueber die einzelnen Punkte des t-Intervalls
    }                                     // Ende des Methode 'solve()'
};  // Ende der Klasse 'dsolve'

// Sub-Klasse 'nBody'
class nBody : public dsolve {
public:
    double G;           // Gravitationskonstante G
    vector<double> m;   // Vektor der Massen der Körper
    unsigned n;         // Anzahl der Körper

    // Überschreibung der virtuellen Funktion dgls
    // Definition der Bewegungsgleichung des n-Körper Problems in drei Dimensionen
    vector<double> dgls(double t, vector<double> u_vec){
        unsigned dim = 3 * n;            // Die n-Körper bilden 3*n Dimensionen
        vector<double> du_dt(2 * dim);   // Das System von 2*3*n DGLs erster Ordnung werden als Standard-Vektor definiert

        for(unsigned i=0; i < dim; ++i){          // 1. Block von 3*n-DGLs erster Ordnung
            du_dt[i] = u_vec[dim + i];            // drdt = v
        }                                         // Ende 1. Block
        for(unsigned i=0; i < n; ++i){                    // 2. Block von 3*n-DGLs erster Ordnung
            double res_x = 0.0, res_y = 0.0, res_z = 0.0; // x,y,z-Werte der DGL des i-ten Körpers
            double xi = u_vec[3 * i];                     // x-Posiion von Körper i
            double yi = u_vec[3 * i + 1];                 // y-Posiion von Körper i
            double zi = u_vec[3 * i + 2];                 // z-Posiion von Körper i
            for(unsigned j=0; j < n; ++j){                   // Schleife j über alle Körper
                if (i != j) {                                // j-ter Körper wirkt auf i-ten Körper
                    double xj = u_vec[3 * j];                // x-Posiion von Körper j
                    double yj = u_vec[3 * j + 1];            // y-Posiion von Körper j
                    double zj = u_vec[3 * j + 2];            // z-Posiion von Körper j
                    // alle j-Körper wirken gravitativ auf den i-ten Körper (in drei Dimensionen, xyz)
                    double r2 = (xj-xi)*(xj-xi) + (yj-yi)*(yj-yi) + (zj-zi)*(zj-zi);
                    double f = G * m[j] / (r2 * sqrt(r2));
                    res_x += f * (xj-xi);
                    res_y += f * (yj-yi);
                    res_z += f * (zj-zi);
                }                                            // Ende if
            }                                                // Ende Schleife j über alle Körper
            du_dt[dim + 3 * i]     = res_x;                  // DGL des i-ten Körpers in x-Richtung
            du_dt[dim + 3 * i + 1] = res_y;                  // DGL des i-ten Körpers in y-Richtung
            du_dt[dim + 3 * i + 2] = res_z;                  // DGL des i-ten Körpers in z-Richtung
        }                                                 // Ende 2. Block
        return du_dt;                        // Rueckgabewert der DGLs als Standard-Vektor
    }                                        // Ende: Definition der Bewegungsgleichung-Funktion

    // 1. Konstruktor mit sechs Argumenten: 'nBody'(Grav.Konstante G, Massenvektor der Körper, Anfangszeit a, Endzeit b, Anzahl der Punkte N, Anfangswerte u_init)
    nBody(double G_ = 1.0, vector<double> m_ = {1,1,1}, double a_ = 0, double b_ = 45.0, int N_ = 10000, vector<double> alpha_ = {1.0,0.0,0.0,-0.5,0.8,0.0,-0.5,-0.8,0.0,0.0,0.5,0.0,0.2,-0.1,0.0,-0.2,-0.4,0.0}) : dsolve(a_, b_, N_, alpha_), G(G_), m(m_) {n = static_cast<unsigned>(m_.size());}

    // 2. Konstruktor mit vier Argumenten: 'nBody'(Anfangszeit a, Endzeit b, Anzahl der Punkte N, Exzentrizität exz)
    nBody(double a_ = 0, double b_ = 45.0, int N_ = 10000, double exz_ = 1.0) : dsolve(a_, b_, N_, alpha_Lagr(exz_)), G(1.0), m({1,1,1}), n(3) {}

    // Methode der Ausgabe der Orte der Körper am Stützpunkt i
    vector<double> r(size_t i) {
        vector<double> w(3 * n, 0.0);
        if (i < t.size()) {
            for(unsigned j=0; j < 3*n; ++j){
                w[j] = ym[i][j];
            }
        }
        return w;
    }
    // Methode der Ausgabe der Geschwindigkeiten der Körper am Stützpunkt i
    vector<double> v(size_t i) {
        vector<double> w(3 * n, 0.0);
        if (i < t.size()) {
            for(unsigned j=0; j < 3*n; ++j){
                w[j] = ym[i][j + 3*n];
            }
        }
        return w;
    }

private:
    // Elliptische Lagrangesche homografischen Lösungen (Exzentrizität exz)
    static vector<double> alpha_Lagr(double exz) {
        double c1 = 1.0 / sqrt(3.0);
        double c2 = sqrt(1.0 - exz);
        double c3 = sqrt(3.0*(1.0 - exz)) / 2.0;
        return {c1,0.0,0.0,-c1/2.0,0.5,0.0,-c1/2.0,-0.5,0.0,0.0,c2,0.0,-c3,-c2/2.0,0.0,c3,-c2/2.0,0.0};
    }

}; // Ende der Unterklasse 'nBody'