//gcc zweiteabgabe.c -o zweiteabgabe -lm -lh2_cip -L.

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

double ExterneKraft(double t){
	double alpha=0.5;
	double w=1.5;

	return alpha*cos(w*t);
	}

double KeineKraft(double t){
	return 0;
}

void dgl(int neq, double w0, double gamma, double *y, double *f, double t,  double (*fext)(double))
{ //für unsere DGL mit neq = 2
    //sicherstellen wie in Vorlesung:(Brauch man das wirklich?
    if (neq!= 2)
    {
        printf ("neq passt nicht.\n");
        abort(); //Bedeutung?->Abbruch
    }
    //Speicherplatz muss zuvor bereitgestellt werden
    f[0] = y[1];
    f[1] = -w0*w0*y[0]-2*gamma*y[1]-(fext(t));
}

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

    dgl(neq,w0,gamma,y,f,t,fext); // Iformationen ü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 (int neq, double h, double w0, double gamma, double *y, double *f, void (*dgl) (int, double,
            double, double*, double*,double,double(double)), double **k, double t, double (*fext)(double)  ) {
	int i;
	double yTemp[neq];
	double tTemp;
	dgl(neq,w0,gamma,y,k[0],t,fext);
	for(i=0;i<neq;i++){
		yTemp[i]=y[i]+h/2*k[0][i];
		tTemp=t+h/2;
		}
	dgl(neq,w0,gamma,yTemp,k[1],tTemp,fext);
	for(i=0;i<neq;i++){
		yTemp[i]=y[i]+h/2*k[1][i];
		tTemp=t+h/2;
			}
	dgl(neq,w0,gamma,yTemp,k[2],tTemp,fext);
	for(i=0;i<neq;i++){
		yTemp[i]=y[i]+h*k[2][i];
		tTemp=t+h;
			}
	dgl(neq,w0,gamma,yTemp,k[3],tTemp,fext);
	for(i=0;i<neq;i++){
		y[i] += (h/2 * (k[0][i]+k[3][i])+h*(k[1][i]+k[2][i]))/3;
			}
}

double Leistung(double *y, double t,double (*fext)(double))
{
    double P = fext(t) * y[1];
    return P;
}

int Leistung_Mittelwert(int n, int neq, double h, double w0, double gamma, double *y, double *f, void (*dgl) (int, double,
            double, double*, double*, double, double(double)), double **k, double t, double (*fext)(double), double (*P)(double*, double, double(double)))
{   
    double x1,x2;
    int j = 0; //Definition Schrittezähler
    double Pges; // Variable zum speichern der addierten Leistungswerte
    double P_Mittelwert; //Varibel zum speichern des Ergebnisses
    x1 = y[0]; //entspricht am Anfang x0
    RK4(neq,h,w0,gamma,y,f,dgl,k,t,fext); // t = t0
    x2 = y[0];
    while (x1 * x2 > 0) //solange kein Nustellenübergang (bzw.Vorzeichenwechsel) zwischen x1 und x2 vorliegt
    {   
        RK4(neq,h,w0,gamma,y,f,dgl,k,t,fext);
        x1 = x2; //Wert übergeben
        x2 = y[0];
        t += h;
    }
    RK4(neq,h,w0,gamma,y,f,dgl,k,t,fext);
    x1 = x2;
    x2 = y[1];
    t += h;
    for (int k = 0; k < n; k++) //wird so oft durchlaufen wie die gewünschte Periodenanzahl n
    {
        Pges = 0;
        for (int i = 0; i < 2; i++)// Zwei Durchläufe, da eine Periode durch jede zweite NST begrenzt wird
        {
            while (x1 * x2 > 0)
            {
                Pges += P(y,t,fext); //Leistung am Punkt x1(erster Punkt nach der NST, wenn nicht x1 selbst die NST ist)
                t += h;
                j ++;
                x1 = x2; //Bestimmung der neuen x-Werte
                RK4(neq,h,w0,gamma,y,f,dgl,k,t,fext);
                x2 = y[0];
            }
            Pges += P(y,t,fext);
            j ++;
            RK4(neq,h,w0,gamma,y,f,dgl,k,t,fext);
            x1 = x2;
            x2 = y[1];
            t += h;
        }
        P_Mittelwert = Pges/ j;
        printf("Mittelwert der Leistung (über Periode %d): %20.5le\n", k, P_Mittelwert);
    }
    return 0;
}



int main()
{
    double h=0.01;//0.001
    double t0=0;
    double x0=1.5;
    double v0=2.75;
    double tend=10;//20
    int neq = 2;

    double *y, *f;
    double t;
    double w0=1;
    double gamma=0.1;
    double *k[4];

    int n; // Dafiniere Variable für die Periodenanzahl

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

    y = (double*)malloc(sizeof(double)*neq);
    f = (double*)malloc(sizeof(double)*neq);
    //k= (double*)malloc(sizeof(double)*neq);
    for(int i=0;i<4;i++){
    	k[i]=(double*)malloc(sizeof(double)*neq);
    }

    y[0] = x0;
    y[1] = v0;
    
    //H4:
    //FILE *fp = fopen("H6_werte", "w");
    //printf("t\tLösung x(t)\t dx(t)/dt\n");
    /*for (t = t0; t <= tend; t+=h)
    {   if (t==0){
    	printf(" %20.5le\t%20.5le\t%20.5le\n", t, x0, v0); //brauch man das überhaupt denke nicht-> wird doppelt ausgegeben!
    }

        printf(" %20.5le\t%20.5le\t%20.5le\n", t, y[0], y[1]);

        //fprintf(fp, " %20.5le\t%20.5le\t%20.5le\n", t, y[0], y[1]);

        RK4(neq,h,w0,gamma,y,f,dgl,k,t,KeineKraft);
    }
    */
    
    //H.6:
    /*for (t = t0; t <= tend; t+=h)
    {
        //printf(" %20.5le\t%20.5le\t%20.5le\n", t, y[0], y[1]);

        fprintf(fp, " %20.5le\t%20.5le\n", t, y[0]);
        
        //fprintf(fp, "%20.5le\t%20.5le\n", t, Leistung(y, t, fext));
        
        RK4(neq,h,w0,gamma,y,f,dgl,k,t,fext);
        //fprintf(fp,"%20.5le\t %20.5le\n",t, fext(t));
    }*/
    
    printf("Bitte geben sie die Periadenanzahl an:\n");
    scanf("%d", &n);
    
    Leistung_Mittelwert(n,neq,h,w0,gamma,y,f,dgl,k,t,fext,Leistung);
    
    //fclose(fp); //Schließe Datei wieder
    
    free(y); free(f);
    for(int i=0;i<4;i++){
    	free(k[i]);
    }

}