//gcc hausaufgabe2.c -o hausaufgabe2 -lm -lh2_cip -L.
//./hausaufgabe2
//Pauline Maas, Ariane Ufer

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "h2.h"

//Globale Variablen definieren, damit nicht so viele Werte übergeben werden müssen
//Hier: Variablen haben Werte aus Aufgabenstellung

double g = 0.1; //gamma
double a = 0.5; //alpha
double w = 1.5; //omega
double x0 = 1.5; //Startwerte
double v0 = 2.75;
int neq = 2;  //Anzahl Gleichungen

//Definiere die externe Kraft für H5
double ExterneKraft(double t) {

	return a * cos(w * t); //Rückgabe der Kraft
}

//Damit man eine RK4 Funktion für alle Aufgaben verwenden kann hier für H4 eine Funktion die als Kraft Null zurückgibt
double KeineKraft(double t) {
	return 0;
}

//rechte Seite der DGL
void dgl(double w0, double *y, double *f, double t, double (*fext)(double))
{
    f[0] = y[1]; //Definition der Hilfsfunktionen
	f[1] = -w0 * w0 * y[0] - 2 * g * y[1] - (fext(t));
}

void euler( double h, double w0, double *y, double *f,
		void (*dgl)(double, double*, double*, double, double(double)), double t,
		double (*fext)(double)) {
	int i;

	dgl(w0, y, f, t, fext); // Informationen über rechte Seite holen

	//Euler-Verfahren anwenden und x_n+1 berechnen
    for (i = 0; i < neq; i++) {
		y[i] += h * f[i];
	}
}

void RK4(double h, double w0, double *y, double *f,
		void (*dgl)(double, double*, double*, double, double(double)),
		double t, double (*fext)(double)) {
	int i;
	double yTemp[neq]; //Temporäre Variable für y und t, die beim RK4 zur Berechnung von k verändert werden
	double tTemp;
    double *k[4]; //Pointer Array um die 4 zweidimensionalen k's abzuspeichern
    for(int i=0;i<4;i++){ //Speicherplatz für die k-Werte zur Verfügung stellen (2 dimensionaler Array)
    	k[i]=(double*)malloc(sizeof(double)*neq);
    }
	//k's berechnen: Hier k[0]=k1, k[1]=k2 usw
	dgl(w0, y, k[0], t, fext); // dgl berechnet k1
	//for Schleife um mit den k1 die nach RK4 veränderten y und t zu berechnen
	for (i = 0; i < neq; i++) {
		yTemp[i] = y[i] + h / 2 * k[0][i];
		tTemp = t + h / 2;
	}
	dgl(w0, yTemp, k[1], tTemp, fext); //dgl berechnet mit verändertem y und t k2
	//for Schleife berechnet neue y und t für k3 nach RK4
	for (i = 0; i < neq; i++) {
		yTemp[i] = y[i] + h / 2 * k[1][i];
		tTemp = t + h / 2;
	}
	dgl(w0, yTemp, k[2], tTemp, fext); //k3 mit Hilfe von Dgl und neuen y und t berechnen
	//y und t für k4 berechnen
	for (i = 0; i < neq; i++) {
		yTemp[i] = y[i] + h * k[2][i];
		tTemp = t + h;
	}
	dgl(w0, yTemp, k[3], tTemp, fext); //k4 berechnen
	//mit allen k's y nach der RK4 vorschrift berechnen
	for (i = 0; i < neq; i++) {
		y[i] += (h / 2 * (k[0][i] + k[3][i]) + h * (k[1][i] + k[2][i])) / 3;
	}
    for(int i=0;i<4;i++){ // Speicherplatz wieder freigeben
    	free(k[i]);
    }
}

//berechnet analytische lösung der Differentialgleichung für H4
double exactH4(double w0, double t) {
	return (exp(-g * t)
			* (x0 * cos(sqrt(w0 * w0 - g * g) * t)
					+ v0 * sin(sqrt(w0 * w0 - g * g) * t)));
}

//berechnet analytische Lösung der DGL für H5
double exactH5(double w0, double t) {
	//berechnen der Koeffizienten siehe Pdf
	double A = a / (sqrt(pow((w * w - w0 * w0), 2) + 4 * g * g * w * w));
	double D = atan(2 * g * w / (w * w - w0 * w0)); //D ist äquivalent zu delta
	double C1 = x0 - A * cos(D);
	double C2 = (v0 - w * A * cos(D) + g * C1) / sqrt(w0 * w0 - g * g);

	return (exp(-g * t)
			* (C1 * cos(sqrt(w0 * w0 - g * g) * t)
					+ C2 * sin(sqrt(w0 * w0 - g * g) * t)) + A * cos(w * t + D));

}

//berechnet Leistung P(t) für H.6
double Leistung(double *y, double t,double (*fext)(double))
{
    double P = fext(t) * y[1];
    return P;
}

//berechnet den Leistungsmittelwert über eine Anzahl von Perioden für H.6
int Leistung_Mittelwert(int n, double h, double w0, double *y, double *f, void (*dgl) (double,double*,
        double*,double, double(double)), double t, double (*fext)(double),double (*P)(double*, double, double(double)))
{   
    double x1,x2;
    int j; //Definition Schrittezähler innerhalb einer Periode
    double Pges; // Variable zum speichern der addierten Leistungswerte
    double P_Mittelwert[n]; //Array zum speichern des Ergebnisses
    double P_M_add = 0; //Hilfsvariablen zur Berechnung des Mittelwerts der einzelnen Mittelwerte
    double Ges_Mittelwert;
    double s_i; // Hilfsvariablen zur Berechnung der Standardabweichung
    double s_add = 0;
    double Ges_s;
    
    x1 = y[0]; //entspricht am Anfang x0
    RK4(h,w0,y,f,dgl,t,fext); // mit t = t0 berechnet nächsten x-Wert
    x2 = y[0]; 
    printf("Omega_0 = %20.5le\n",w0);
    while (x1 * x2 > 0) //solange kein Nustellenübergang (bzw.Vorzeichenwechsel) zwischen x1 und x2 vorliegt
    {   
        RK4(h,w0,y,f,dgl,t,fext);
        x1 = x2; //Wert übergeben
        x2 = y[0]; // neuen Wert übernehmen
        t += h;
    }
    RK4(h,w0,y,f,dgl,t,fext); // einen weiteren Schritt berechnen, so dass x1 den ersten Wert nach der NST annimmt (wenn nicht sogar x1 = 0 wird)
    x1 = x2;
    x2 = y[1];
    t += h;
    for (int l = 0; l < n; l++) //wird so oft durchlaufen wie die gewünschte Periodenanzahl n
    {
        Pges = 0;
        j = 0;
        for (int i = 0; i < 2; i++)// Zwei Durchläufe, da eine Periode durch jede zweite NST begrenzt wird
        {
            while (x1 * x2 > 0) // überprüft auf Vorzeichenwechsel (NST)
            {
                Pges += P(y,t,fext); //Leistung am Punkt x1 wird berechnet und zu Pges hinzugefügt
                t += h;
                j ++; // zählt die Anzahl der berechneten Leistungswerte
                x1 = x2; //Bestimmung der neuen x-Werte
                RK4(h,w0,y,f,dgl,t,fext);
                x2 = y[0];
            }
            Pges += P(y,t,fext); //den letzten Wert vor der NST berechnen
            j ++;
            RK4(h,w0,y,f,dgl,t,fext); //und einen Schritt weiter gehen, um aus dem Nullstellenübergang herauszukommen 
            x1 = x2;
            x2 = y[1];
            t += h;
        }
        P_Mittelwert[l] = Pges/ j;
        printf("Mittelwert der Leistung (Periode %d): %20.5le\n", l+1 , P_Mittelwert[l]);
        P_M_add += P_Mittelwert[l];
    }
    Ges_Mittelwert = P_M_add/n;
    for (int i = 0; i < n; i++) // Berechnung Standardfehler
    {   
         s_i = pow (P_Mittelwert[i]-Ges_Mittelwert,2);
         s_add += s_i;
    }
    Ges_s = sqrt(s_add/n*(n-1));
    printf("Mittelwert über alle %d Perioden:%20.5le\n",n,Ges_Mittelwert);
    printf("Standardfehler:%20.5le\n",Ges_s);
    return 0;
}


int main() {
	double h = 0.001; //Schrittweite h
	double t0 = 0; //Startzeit

	double tend = 1; //Endzeit einstellbar

	double *y, *f; //Pointer auf denen die Verktoren y und f gespeichert werden
	double t; //Variable ,die momentane Zeit beschreibt
	double w0 = 1; //omega null
	int n; // Dafiniere Variable für die Periodenanzahl

	/*printf("Bitte geben Sie h, x0, v0 und tend, sowie omega0 und gamma ein:\n");
	 scanf("%le %le %le %le %le %le", &h, &t0, &x0, &v0, &tend, &w0, &g);*/

	y = (double*) malloc(sizeof(double) * neq); //Speicherplatz für y allocieren(y hat dim neq (hier 2))
	f = (double*) malloc(sizeof(double) * neq); //Speicherplatz für f allocieren (dim 2 siehe y)
	
    //startwerte zuweisen
	y[0] = x0;
	y[1] = v0;
    
	//H4:

     printf("t\tLösungen von H.4: x(t)\tdx(t)/dt\n");
	 //Für jeden Zeitpunkt RK4 ausführen
	 for (t = t0; t <= tend; t += h) {
	 printf(" %20.5le\t%20.5le\t%20.5le\n", t, y[0], y[1]);//Ergebnisse als Tabelle ausgeben (y[0] ist x, y[1] ist xPunkt)
     //Oder:
	 //printf(" %20.5le\t%20.5le\n", t, abs(y[0]-exactH4(w0,t))); //Vergleich von numerischer und analytischer Lösung für x

	 RK4(h, w0, y, f, dgl, t, KeineKraft); //RK4 ohne anregende Kraft
	 }
	 
	//H5:

	printf("t\tLösung von H.5: x(t)\tdx(t)/dt\n");
	for (t = t0; t <= tend; t += h) {
		printf(" %20.5le\t%20.5le\t%20.5le\n", t, y[0], y[1]); //Ergebnisse in Tabelle ausgeben
		//Oder:
		//printf(" %20.5le\t%20.5le\n", t, y[0]- exactH5(w0,t)); //Vgl numerische und analytische Lösung



		RK4(h, w0, y, f, dgl, t, ExterneKraft); //Rk4 mit Anregung, die in ExterneKraft gespeichert ist
	}
	
	//H.6:
	
	printf("Bitte geben sie die Periodenanzahl an:\n");
    scanf("%d", &n);
    
    for (double w0_scan = 1; w0_scan < 3.2; w0_scan += 0.2) //stellt den scan für w0 dar
    {
        Leistung_Mittelwert(n, h, w0_scan, y, f, dgl, t0, fext, Leistung);
    }
    
    
    //Speicherplatz frei machen
	free(y);
	free(f);

}