/* Planetenbewegungen in drei Dimensionen
 * basierend auf dem 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 verwendet
 * 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 "Planeten.dat"
 */

#include "N-BodyProblem.hpp" // Header-Datei des n-Körperproblem in drei Dimensionen
#include <iostream>          // Standard Input- und Output Bibliothek in C, z.B. fprintf(...)
#include <vector>            // Vector-Container der Standardbibliothek
#include <string>            // Standard-Header-Datei für Zeichenketten
#include <map>               // Standard-Header-Datei des assoziativer Containers 'map' für Key-Value Zuordnungen
using namespace std;         // Benutze den Namensraum std

// Struktur zum georneten Speichern der Himmelskörper-Parameter
struct Body {
    double gm;          // gm=G*m: Gravitationskonstante * Masse des Körpers
    vector<double> r;   // Anfangsort in drei Dimensionen
    vector<double> v;   // Anfangsogeschwindigkeit in drei Dimensionen
};

// Hauptfunktion
int main(){
    // Definition der gm-Parameter und Anfangswerte der Himmelskörper unseres Sonnensystems in einer C++-map
    // Ephemeriden-Datei 'de421.bsp': 24. Dezember 2000 um 20:30 Uhr UTC)
    // Länge in Astronomischen Einheitn (AU = 149.597.870,7 km), Zeit in Tagen [d]
    map<string, Body> bodies_data = {
        {"sun", {0.00029591220828559, {-0.00469206964217851, -0.00455284299129719, -0.00180351606161285}, {0.00000808182754769, -0.00000364991309399, -0.00000177739810148}}},
        {"mercury", {0.00000000004912480, {0.00869081881722255, -0.41212824812482196, -0.22090270075603882}, {0.02248815298088397, 0.00281999190824476, -0.00082493351551472}}},
        {"venus", {0.00000000072434525, {0.58492431879299430, 0.39106068887197776, 0.13886288700490920}, {-0.01176860984756512, 0.01468495634135895, 0.00735174095931183}}},
        {"earth", {0.00000000088876924, {-0.06209037112534756, 0.89629552530222610, 0.38875814156074100}, {-0.01745436390436983, -0.00098293944800350, -0.00042570411094977}}},
        {"moon", {0.00000000001093189, {-0.06238174828301801, 0.89381342029917210, 0.38778197364226596}, {-0.01689334413161172, -0.00103583077929171, -0.00050101150118707}}},
        {"mars", {0.00000000009549535, {-1.65650069587999330, 0.01013256724249670, 0.04958858122190123}, {0.00024579442831315, -0.01163534295071834, -0.00534326234956417}}},
        {"jupiter barycenter", {0.00000028253459095, {1.84307596674717920, 4.32528771663242900, 1.80908857439581870}, {-0.00711591051273230, 0.00280462876522796, 0.00137544423420692}}},
        {"saturn barycenter", {0.00000008459706073, {4.71735462953881500, 7.27073386854747500, 2.80001882890842420}, {-0.00506407973083511, 0.00258153583418907, 0.00128420716954018}}},
        {"uranus barycenter", {0.00000001292024826, {15.35097038912193100, -11.59975177639601200, -5.29752281748167600}, {0.00248407812132207, 0.00261666373062273, 0.00111089476005838}}},
        {"neptune barycenter", {0.00000001524357348, {17.71921742890633700, -22.37141413515812300, -9.59788796335357700}, {0.00251753419769610, 0.00174985139276685, 0.00065354736058049}}},
        {"pluto barycenter", {0.00000000000217844, {-8.78866167066838600, -28.36872537395957700, -6.20501655731820000}, {0.00306982805604962, -0.00102661679740455, -0.00124530173577011}}}
    };

    double G = 1.0;                             // Gravitationskonstante G=1, da Massenwerte in G*m
    double a = 0;                               // Untergrenze des Zeit-Intervalls
    double b = 365.2307715;                     // Obergrenze des Intervalls in Tagen (1 Jahr aprox. 365.2307715 Tage)
    int N = 10000;                              // Anzahl der Punkte der Stützstellen 1000000

    FILE *ausgabe;                              // Deklaration eines Files fuer die Ausgabedatei
    ausgabe = fopen("Planeten.dat", "w+");      // Ergebnisse werden in die Datei 'Planeten.dat' geschrieben

    // String-Vektor der zu simulierenden Himmelskörper
//    const vector<string> set_bodies {"sun","mercury","venus","earth","moon","mars"};
    const vector<string> set_bodies {"sun","mercury", "venus","earth","moon","mars","jupiter barycenter",
              "saturn barycenter", "uranus barycenter", "neptune barycenter", "pluto barycenter"};
    const size_t n = set_bodies.size();        // Anzahl der zu simulierenden Himmelskörper

    vector<double> m(n);              // Vektor der Massen der Körper (gm-Werte, gm=G*m)
    vector<double> u_init(6*n);       // Vektor der Anfangswerte für die Orte und Geschwindigkeiten der n Körper in drei Dimensionen
    for (size_t i = 0; i < n; ++i) {                                       // for-Schleife über alle Körper
        const auto& data = bodies_data.at(set_bodies[i]);                  // Daten des i-ten Körpers als Referenz auf die C++-map
        m[i] = data.gm;                                                    // gm-Wert in den m-Vektor schreiben
        copy(data.r.begin(), data.r.end(), u_init.begin() + (i * 3));      // Orte an den Anfang von u_init schreiben (Position i * 3)
        copy(data.v.begin(), data.v.end(), u_init.begin() + (n + i) * 3);  // Geschwindigkeiten in die zweite Hälfte von u_init schreiben (Position: (n + i) * 3)
    }                                                                      // Ende for-Schleife über alle Körper

    nBody L = nBody(G,m,a,b,N,u_init); // Konstruktor 1 bilded eine Instanz der Klasse nBody
    L.solve();                         // Methode des Runge-Kutta-Verfahrens ausführen

    fprintf(ausgabe, "# 0: Index i \n# 1: t-Wert \n# 2:x_1 \n# 3:vx_1 \n# 4:y_1 \n# 5:vy_1 \n# 6:z_1 ... \n"); // Beschreibung der ausgegebenen Groessen
    for(size_t i=0; i  < L.t.size(); i += 1){                         // for-Schleife über die Zeitgitterpunkte (i += 10 -> nur jeden 10-ten Wert ausgeben)
        fprintf(ausgabe, "%3ld %19.15f ",i,L.t[i]);                   // Ausgabe Index und Zeit
        for(unsigned j=0; j  < 2*3*L.n; ++j){                         // for-Schleife zur Ausgabe der Orte und Geschw. der n-Körper
            fprintf(ausgabe, "%19.15f %19.15f ",L.r(i)[j],L.v(i)[j]); // Ausgabe xyz-Werte und Geschwindigkeiten der Körper
        }                                                             // Ende der for-Schleife; n-Körper
        fprintf(ausgabe, "\n");                                       // Neue Zeile
    }                                                                 // Ende der for-Schleife; Zeitgitterpunkte
} // Ende der Hauptfunktion