/* Header-Datei
 * 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
 * Berechnung der Lösung eines Systems von sechs Differentialgleichungen erster Ordnung
 * 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)
 * 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) 'geodesic'
 * Zeitentwicklung fuer unterschiedliche Eigenzeit-Werte tau in [a,b]
 * Unterklasse 'geodesic'
 * Konstruktor mit fuenf Argumenten: 'geodesic'(Masse des schwarzen Lochs 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}
 */

#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 (Eigenzeit)
    double b;                  // Obergrenze des Intervalls (Eigenzeit)
    int N;                     // Anzahl der Punkte
    vector<double> alpha;      // Anfangswerte

    vector<vector<double>> ym; // Lösungsmatrix
    vector<double> tau;        // Eigenzeit-Vektor

    // Virtuelle Funktion dgls, wird an anderer Stelle definiert
    virtual vector<double> dgls(double tau, 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_) {
        tau.push_back(a_);     // Zum Eigenzeit-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 dtau zwischen den aequidistanten Punkten des tau-Intervalls (h=dtau)
        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 tau-Intervalls
            for (int j = 0; j < n; ++j) k1[j] = h * dgls(tau[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(tau[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(tau[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(tau[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
            tau.push_back(a + (i + 1) * h); // Zum Eigenzeit-Vektor die neuen Wert eintragen
        }                                   // Ende for-Schleife ueber die einzelnen Punkte des tau-Intervalls
    }                                       // Ende des Methode 'solve()'
};  // Ende der Klasse 'dsolve'

// Sub-Klasse 'geodesic'
class geodesic : public dsolve {
public:
    double M; // Masse M des Schwarzen Lochs

    // Überschreibung der virtuellen Funktion dgls, Definition der Geodätengleichung
    vector<double> dgls(double tau, vector<double> u_vec){
        vector<double> du_dt(6);   // Die sechs DGLs erster Ordnung werden als Standard-Vektor definiert
        du_dt[0] = u_vec[3];       // Sechs DGLs erster Ordnung
        du_dt[1] = u_vec[4];       // ...
        du_dt[2] = u_vec[5];
        du_dt[3] = 2*M*u_vec[3]*u_vec[4]/(u_vec[1]*(2*M - u_vec[1]));
        du_dt[4] = -M*pow(u_vec[4], 2)/(u_vec[1]*(2*M - u_vec[1])) + M*pow(u_vec[3], 2)*(2*M - u_vec[1])/pow(u_vec[1], 3) + pow(u_vec[5], 2)*(-2*M + u_vec[1]);
        du_dt[5] = -2*u_vec[4]*u_vec[5]/u_vec[1];
        return du_dt;              // Rueckgabewert der DGLs als Standard-Vektor
    }                              // Ende: Definition der Bewegungsgleichung-Funktion

    // 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}
    geodesic(double M_ = 1.0, double a_ = 0, double b_ = 800.0, int N_ = 10000, vector<double> alpha_ = {0.0,10.0,0.0,1.208356321620407,0.0,0.041}) : dsolve(a_, b_, N_, alpha_), M(M_) {}
};  // Ende der Unterklasse 'geodesic'