// 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, Optimierte Version a1
 * Plotten mittels Python Jupyter Notebook AngetriebenesPendel.ipynb, oder PythonPlot_GetriebPendel_attr_a.py
 */

#include <iostream>
#include <cmath>
#include <vector>
#include <fstream>
#include <iomanip>
#include <chrono>
#include <omp.h>

// Klasse zum Berechnen der Attraktorwerte des periodisch angetriebenen Pendels
class Attraktor {
private:
    unsigned N;
    double alpha_x;
    double alpha_v;
    double A;
    double Omega;
    unsigned N_attr;
    double beta = 0.5;

     // Definition der Bewegungsgleichung des periodisch angetriebenen Pendels als private Member-Funktionen der Klasse
    void f(double t, double x, double v, double& dx, double& dv) const {
        dx = v;
        dv = A * std::cos(Omega * t) - beta * v - std::sin(x);
    }

    // Funktion um bei Ueberschlaegen des Pendels den Winkel wieder in den Bereich [-pi,pi] zu verschieben
    double wrap_angle(double angle) const {
        double wrapped = std::fmod(angle + M_PI, 2.0 * M_PI);
        if (wrapped < 0) wrapped += 2.0 * M_PI;
        return wrapped - M_PI;
    }

public:
    // Konstruktor initialisiert die Instanzvariablen
    Attraktor(unsigned N_, std::pair<double, double> alpha_, double A_, double Omega_, unsigned N_attr_)
        : N(N_), alpha_x(alpha_.first), alpha_v(alpha_.second), A(A_), Omega(Omega_), N_attr(N_attr_) {}

    // Öffentliche Member-Funktion der Klasse (Berechnung der Attraktorwerte des Pendels)
    std::vector<double> calc_attr_werte() const {
        double x = alpha_x;
        double v = alpha_v;

        double t_end = N_attr * 2.0 * M_PI / Omega;
        unsigned i_schritt = N / N_attr;
        double h = t_end / N;
        double h_half = h * 0.5;

        std::vector<double> attr_werte;
        attr_werte.reserve(1 + 2 * (N_attr - 70));
        attr_werte.push_back(A);

        double k1_x, k1_v, k2_x, k2_v, k3_x, k3_v, k4_x, k4_v;
        double x_tmp, v_tmp;

        // for-Schleife ueber die einzelnen Punkte des t-Intervalls
        for (unsigned i = 0; i < N; ++i) {
            double t = i * h;
            f(t, x, v, k1_x, k1_v);

            x_tmp = x + h_half * k1_x;
            v_tmp = v + h_half * k1_v;
            f(t + h_half, x_tmp, v_tmp, k2_x, k2_v);

            x_tmp = x + h_half * k2_x;
            v_tmp = v + h_half * k2_v;
            f(t + h_half, x_tmp, v_tmp, k3_x, k3_v);

            x_tmp = x + h * k3_x;
            v_tmp = v + h * k3_v;
            f(t + h, x_tmp, v_tmp, k4_x, k4_v);

            x += (h / 6.0) * (k1_x + 2.0 * k2_x + 2.0 * k3_x + k4_x);
            v += (h / 6.0) * (k1_v + 2.0 * k2_v + 2.0 * k3_v + k4_v);

            // Herausfiltern der Attraktorwerte
            unsigned aktueller_schritt = i + 1;
            if ((aktueller_schritt / i_schritt >= 70) && (aktueller_schritt % i_schritt == 0)) {
                attr_werte.push_back(v);
                attr_werte.push_back(wrap_angle(x));
            }
        }
        return attr_werte;
    }
};

int main() {
    // Starten der Zeitmessung
    auto start = std::chrono::high_resolution_clock::now();

    unsigned Anz_Pendel = 2000;
    double Omega = 2.0 / 3.0;
    unsigned Anz_Punkte = 150000;
    unsigned Anz_Punkte_attr = 100;
    double Anf_A = 0.95;
    double End_A = 1.6;
    double dA = (End_A - Anf_A) / (Anz_Pendel - 1);

    std::vector<std::vector<double>> attr_wert(Anz_Pendel);

    // Parallelisierte Schleife über alle Pendel mit dynamischer Verteilung der Aufgaben auf die CPU-Kerne
    #pragma omp parallel for schedule(dynamic)
    for (unsigned n = 0; n < Anz_Pendel; ++n) {
        Attraktor Loes{Anz_Punkte, {0.0, 0.0}, Anf_A + n * dA, Omega, Anz_Punkte_attr};
        attr_wert[n] = Loes.calc_attr_werte();
    }

    // Ausgabe in die Datei
    std::ofstream ausgabe("GetriebenesPendel_attr_dia_a.dat");
    ausgabe << std::fixed << std::setprecision(15);

    for (unsigned n = 0; n < Anz_Pendel; ++n) {
        ausgabe << attr_wert[n][0] << " ";
        for (size_t i = 1; i < attr_wert[n].size(); i += 2) {
            ausgabe << attr_wert[n][i] << " ";
        }
        for (size_t i = 2; i < attr_wert[n].size(); i += 2) {
            ausgabe << attr_wert[n][i] << " ";
        }
        ausgabe << "\n";
    }

    auto end = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> diff = end - start;
    std::cout << "Das Programm benoetigte: " << diff.count() << " Sekunden." << std::endl;

    return 0;
}