// Das periodisch angetriebenen Pendel
/* Beispiel eines schwingenden Systems aus dem Bereich der Mechanik
 * welches bei gewissen Parameterkonstellationen deterministisch chaotische Bewegungen zeigt
 * Berechnung der Loesung eines Systems von Differentialgleichungen
 * mittels Runge-Kutta Verfahren (4.Ordnung)
 * Attraktordiagramm vieler Pendel mit unterschiedlichen Amplitudenwerten (Feigenbaumdiagramm)
 * Klasse attraktor optimiert mit OpenMP
 * Plotten mittels Python Jupyter Notebook AngetriebenesPendel.ipynb, oder PythonPlot_GetriebPendel_attr.py
 */

#include <iostream>   // Standard Input- und Output Bibliothek
#include <cmath>      // Bibliothek für mathematisches (e-Funktion, Betrag, ...)
#include <vector>     // Vector-Container der Standardbibliothek
#include <omp.h>      // OpenMP zum parallelen Rechnen
#include <ctime>      // Für time()
using namespace std;  // Benutze den Namensraum std

// Klasse zum Berechnen der Attraktorwerte des periodisch angetriebenen Pendels
class Attraktor {
    // Private Instanzvariablen (Daten-Member) der Klasse
    unsigned N = 10000;           // Anzahl der Punkte in die das t-Intervall aufgeteilt wird
    vector<double> alpha = {0,0}; // Zwei Anfangswerte (Ort und Geschwindigkeit des Pendels) bei t=a
    double A = 0.9;               // Amplitude A
    double Omega = 2.0 / 3.0;     // Frequenz Omega
    unsigned N_attr = 80;         // Anzahl der Punkte im Attraktordiagramm (minus 50 Punkte fuer die Einschwingphase)

    // Definition der Bewegungsgleichung des periodisch angetriebenen Pendels als private Member-Funktionen der Klasse
    vector<double> f(double t, vector<double> u_vec, double A, double Omega) {
        double beta {0.5};       // Stokessche Reibungskoeffizient
        vector<double> du_dt(2);
        du_dt[0] = u_vec[1];
        du_dt[1] = A * cos(Omega * t) - beta * u_vec[1] - sin(u_vec[0]);
        return du_dt;
    }

public:
    // Konstruktor initialisiert die Instanzvariablen
    Attraktor(unsigned N_, vector<double> alpha_, double A_, double Omega_, unsigned N_attr_) : N(N_),alpha(alpha_),A(A_),Omega(Omega_),N_attr(N_attr_){}

    // Öffentliche Member-Funktion der Klasse (Berechnung der Attraktorwerte des Pendels)
    vector<double> calc_attr_werte () {
        vector<double> k1(2), k2(2), k3(2), k4(2); // Deklaration der vier Runge-Kutta Parameter
        vector<double> y = alpha;                  // Aktueller Zustand des Pendels
        vector<double> y_temp(2);                  // Temporärer Hilfsvektor
        double t_end = N_attr * 2.0*M_PI/Omega;    // Obergrenze des t-Intervalls [0,t_end] als Vielfaches der Zeitabstaende t_ni = 2*M_PI/Omega
        int i_schritt = N / N_attr;                // Anzahl der Iterationen fuer einen Zeitabstand t_ni = 2*M_PI/Omega
        double h = t_end / N;                      // Abstand dt
        double t = 0;                              // Aktuelle Zeit
        vector<double> attr_werte;                 // Speichert nur die relevanten Attraktorwerte für das konkrete Pendel
        attr_werte.push_back(A);                   // Erster Eintrag ist die Amplitude A der äußeren Kraft

        // for-Schleife ueber die einzelnen Punkte des t-Intervalls
        for (unsigned i = 0; i < N; ++i) {
            vector<double> f_val = f(t, y, A, Omega);
            for (int j = 0; j < 2; ++j){ k1[j] = h * f_val[j];}
            for (int j = 0; j < 2; ++j){ y_temp[j] = y[j] + k1[j] / 2.0;}
            f_val = f(t + h / 2.0, y_temp, A, Omega);
            for (int j = 0; j < 2; ++j){ k2[j] = h * f_val[j];}
            for (int j = 0; j < 2; ++j){ y_temp[j] = y[j] + k2[j] / 2.0;}
            f_val = f(t + h / 2.0, y_temp, A, Omega);
            for (int j = 0; j < 2; ++j){ k3[j] = h * f_val[j];}
            for (int j = 0; j < 2; ++j){ y_temp[j] = y[j] + k3[j];}
            f_val = f(t + h, y_temp, A, Omega);
            for (int j = 0; j < 2; ++j){ k4[j] = h * f_val[j];}
            for (int j = 0; j < 2; ++j){ y[j] = y[j] + (k1[j] + 2.0 * k2[j] + 2.0 * k3[j] + k4[j]) / 6.0;}
            t = (i + 1) * h;

            // Herausfiltern der Attraktorwerte
            unsigned aktueller_schritt = i + 1;
            if ( (aktueller_schritt / i_schritt >= 70) && ( (i + 1) % i_schritt == 0) ) {
                attr_werte.push_back(y[1]); // Fuellen der Attraktordiagramm-Werte (Winkelgeschwindigkeiten des Pendels)
            }
        }
        return attr_werte;
    }   // Ende der Member-Funktionen (Berechnung der Attraktorwerte des Pendels)
};      // Ende der Klasse

// Hauptfunktion
int main() {
    int startTime = time(NULL);     // Starten der Zeitmessung
    unsigned Anz_Pendel = 2000;     // Anzahl der Pendel im Attraktordiagramm
    double Omega = 2.0 / 3.0;       // Frequenz Omega
    unsigned Anz_Punkte = 150000;   // Anzahl der Punkte im t-Intervall
    unsigned Anz_Punkte_attr = 100; // Anzahl der Punkte im Attraktordiagramm

    double Anf_A = 0.95;                            // Anfangswert Amplitude
    double End_A = 1.6;                             // Endwert Amplitude
    double dA = (End_A - Anf_A) / (Anz_Pendel - 1); // Schrittweite dA der Amplitude A im Attraktordiagramm

    // Deklaration einer double Vektor-Matrix zum speichern der Attraktorwerte der unterschiedlichen Pendel
    vector<vector<double>> attr_wert(Anz_Pendel);

    // Parallelisierte Schleife über alle Pendel
    #pragma omp parallel for
    for (unsigned n = 0; n < Anz_Pendel; ++n) {
        // Pendel-Konstruktor (Anzahl der Zeit-Punkte: Anz_Punkte, Anfangsort und Geschwindigkeit: {0,0}, Amplitude A und Frequenz Omega der ausseren periodischen Kraft, Anzahl der Punkte im Attraktordiagramm)
        Attraktor Loes {Anz_Punkte, {0, 0}, Anf_A + n*dA, Omega, Anz_Punkte_attr};
        attr_wert[n] = Loes.calc_attr_werte(); // Berechnung der Attraktorwerte des Pendels
    }

    // Ausgabe in die Datei
    FILE *ausgabe = fopen("GetriebenesPendel_attr_dia.dat", "w+");
    for (unsigned n = 0; n < Anz_Pendel; ++n) {
        for (size_t i = 0; i < attr_wert[n].size(); ++i) {
            fprintf(ausgabe, "%19.15f ", attr_wert[n][i]);
        }
        fprintf(ausgabe, "\n");
    }
    fclose(ausgabe);

    cout << "Das Programm benoetigte: " << time(NULL) - startTime << " Sekunden." << endl;
}