/* Siehe VARTC 2026/27 (siehe http://itp.uni-frankfurt.de/~hanauske/VARTC/index.html , Vorlesung 6)
 * basierend auf VPROG 2026 (siehe https://itp.uni-frankfurt.de/~hanauske/VPROG/index.html , Vorlesung 8,9,10)
 * C++ Programm zum Loesen der Geodaetengleichung
 * in vorgegebener Schwarzschild-Raumzeit (nicht-rotierendes, ungeladenes Schwarzes Loch)
 * Die Loesung benutzt die Runge-Kutta Ordnung vier Methode (RK4)
 * Koordinaten (Koordinatenzeit: t=y_1 , Radius: r=y_2, Winkel: theta=pi/2=const,  Winkel: phi=y_3)
 * Die zeitliche Entwicklung erfolgt in Eigenzeit tau von y_1(tau), y_2(tau), y_3(tau)
 * Vewendet die Header-Datei GeodesicSchwarzschild.hpp
 * Konstruktor mit fuenf Argumenten: 'geodesic'(Masse M, Anfangswert der Eigenzeit a, Endwert der Eigenzeit b, Anzahl der Punkte N, Anfangswerte u_init)
 * u_init = alpha = {t,r,phi,d_t,d_r,d_phi}
 * Ausgabe zum Plotten mittels Python Jupyter Notebook GeodSchwarzschildII_OGRePy.ipynb
 * bzw. mittels des Python-Skrips GeodSchwarzschild.py
 */

#include "GeodesicSchwarzschild.hpp" // Header-Datei zum Loesen der Geodaetengleichung
#include <iostream>                  // Standard Input- und Output Bibliothek
#include <vector>                    // Vector-Container der Standardbibliothek
using namespace std;                 // Benutze den Namensraum std

// Hauptfunktion
int main(){
    double a = 0;                          // Untergrenze des Eigenzeit-Intervalls
    double b = 800.0;                      // Obergrenze des Eigenzeit-Intervalls
    int N = 100000;                        // Anzahl der Punkte

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

    geodesic L = geodesic(1,a,b,N,{0.0,10.0,0.0,1.208356321620407,0.0,0.041}); // Konstruktor erzeugt eine Instanz der Klasse geodesic

    L.solve();                             // Methode des Runge-Kutta-Verfahrens ausführen

    fprintf(ausgabe, "# 0: Index i \n# 1: tau \n# 2: t \n# 3:r \n# 4:phi \n# 5:d_t \n# 6:d_r  \n# 7:d_phi \n"); // Beschreibung der ausgegebenen Groessen
    for(size_t i = 0; i  < L.tau.size(); i += 10){       // for-Schleife über die Zeitgitterpunkte (nur jeden 10-ten Wert ausgeben)
        fprintf(ausgabe, "%3zu %19.15f ", i, L.tau[i]);  // Ausgabe Index und Zeit
        for (const double value : L.ym[i]) {             // Bereichsbasierte for-Schleife zur
            fprintf(ausgabe, "%19.15f ", value);         // Ausgabe von Ort und Geschw. des Probekoerpers
        }                                                // Ende der for-Schleife; n-Körper
        fprintf(ausgabe, "\n");                          // Neue Zeile
    }                                                    // Ende der for-Schleife; Zeitgitterpunkte
} // Ende der Hauptfunktion