#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <math.h>
#include <string.h>
using namespace std;
typedef struct {
	double N; //Parameter anzahl der Solitonen
	double d; //Schrittweiten von Zeit und Ort
	double h;
	double xmin, xmax, tmin, tmax; //Intervalle
	double x; //momentanes x und t
	double t;
	int jmax; //Maximale Anzahl der Schritte
	int nmax;

} soli_param;
char dateiname[256]="plotu1";
double *u0; //Pointer auf denen u(n=0,j) gespeichert wird
double *uj1; //Pointer auf denne u(n=1,j) gespeichert wird
//u0 mit Glg (3) berechnet
void ustart(soli_param *sp) {
	double N = sp->N; //Variablen aus Struct holen
	double jmax = sp->jmax;
	double h = sp->h;
	double x = sp->xmin;
    
	for (int j = 0; j <= jmax; j++) {

		u0[j] = -N * (N + 1) / (pow(cosh(x), 2));
		x=x+h;

	}

}
//analytische Lösungen berechnen
double u_analytic(soli_param *sp){
    double x=sp->x; //Variablen aus Struct holen
    double t=sp->t;
    double N= sp->N;
    //für verschiedene N
    if (N==1){
        return (-2/pow(cosh(x-4*t),2));
    }
    else if (N==2){
        return (-12*(3+4*cosh(2*x-8*t)+cosh(4*x-64*t))/(pow((3*cosh(x-28*t)+cosh(3*x-36*t)),2)));
    }
    else{
        printf("Keine analytische Lösung für dieses N");
        return 0;
    }
}
//uj1 nach Glg (8) berechnen
void ujeins(soli_param *sp) {
	double jmax = sp->jmax; //Variablen aus Struct holen
	double d = sp->d;
	double h = sp->h;
    double x=sp->xmin;
    //für alle j
	for (int j = 0; j < jmax; j++) {
        //Randbed: u0(j-2) fällt weg wenn j=1 ist 
		if (j == 1) {
			uj1[j] = (d * 6 * u0[j] * (u0[j + 1] - u0[j - 1]) / (2 * h))
					- (d * (u0[j + 2] + 2 * u0[j - 1]-2*u0[j+1]) / (2 * pow(h, 3)))
					+ u0[j];

		}
		//Randbed: u0(j-2) und u0(j-1) fällt weg wenn j=0 ist 
		else if (j == 0) {
			uj1[j] = (d * 6 * u0[j] * (u0[j + 1]) / (2 * h))
					- (d * (u0[j + 2]-2*u0[j+1]) / (2 * pow(h, 3))) + u0[j];

		}
		//Randbed: u0(j+2) fällt weg wenn j=jmax-2 ist 
		else if (j == jmax - 2) {
			uj1[j] = (d * 6 * u0[j] * (u0[j + 1] - u0[j - 1]) / (2 * h))
					- (d * (-2*u0[j+1]+2 * u0[j - 1] - u0[j - 2]) / (2 * pow(h, 3)))
					+ u0[j];

		}
		//Randbed: u0(j+2) und u0(j+1) fällt weg wenn j=jmax -1 ist 
		else if (j == jmax-1) {
			uj1[j] = (d * 6 * u0[j] * (-u0[j - 1]) / (2 * h))
					- (d * (2 * u0[j - 1] - u0[j - 2]) / (2 * pow(h, 3)))
					+ u0[j];
		}
		else {
			uj1[j] = ((d * 6 * u0[j] * (u0[j + 1] - u0[j - 1])) / (2 * h))
					- (d * (u0[j + 2] -2*u0[j+1]+ 2 * u0[j - 1] - u0[j - 2])
							/ (2 * pow(h, 3))) + u0[j];
		}
		

	}

}
void nextu(soli_param *sp) {
	int jmax = sp->jmax;
	int nmax = sp->nmax;
	double d = sp->d;
	double h = sp->h;
	double * unm = u0; //unm=u(n-1) wird auf u0 gesetzt 
	double * un = uj1; //analog dazu auf uj1
	double unp[jmax]; //unp=u(n+1) ist also im ersten schritt in der Schleife u(n=2)
    double t=sp->tmin;
    double x=sp->xmin;
    int n;
    //Berechnen für alle n bis tmax erreicht ist
	for ( n = 2; n < nmax; n++) {
        
		for (int j = 0; j < jmax; j++) {
            //Randbedingungen werden genauso benutzt wie in ujeins
			if (j == 1) {
				
                unp[j]=(2*(d/h) *(un[j+1]+un[j]+un[j-1])*(un[j+1]-un[j-1]))-(d/((h*h*h))*(un[j+2]-2*un[j+1]+2*un[j-1]) ) +unm[j];  

			} else if (j == 0) {
                unp[j]=(2*(d/h) *(un[j+1]+un[j])*(un[j+1]))-(d/(h*h*h)*(un[j+2]-2*un[j+1]) ) +unm[j];  

			} else if (j == jmax - 2) {
                unp[j]=(2*(d/h) *(un[j+1]+un[j]+un[j-1])*(un[j+1]-un[j-1]))-(d/((h*h*h))*(-2*un[j+1]+2*un[j-1]-un[j-2]) ) +unm[j];  

			} else if (j == jmax - 1) {
                unp[j]=(2*(d/h) *(un[j]+un[j-1])*(-un[j-1]))-(d/((h*h*h))*(2*un[j-1]-un[j-2]) ) +unm[j];  

			} else {
				unp[j]=(2*(d/h) *(un[j+1]+un[j]+un[j-1])*(un[j+1]-un[j-1]))-(d/((h*h*h))*(un[j+2]-2*un[j+1]+2*un[j-1]-un[j-2]) ) +unm[j];                     
			}
		}
		//Variablen fuer den nächsten Schritt tauschen (komponentenweise, da Pointer)
		for (int j = 0; j < jmax; j++) {
    		  unm[j]=un[j];
    		  un[j]=unp[j];
            }
		

	}
    //Speichern der Werte
    soli_param diffp= *sp; //Neue Variablen in struct um neues x und t u_analytic mitzgeben
    diffp.t=sp->tmax;
    diffp.x=sp->xmin;
    double diff;
	FILE *fp= fopen("numplot0_1_2.csv", "w"); //Datei öffnen
	for(int j=0;j<jmax;j++){
		diff= fabs(un[j]-u_analytic(&diffp)); //Differenz berechnen
		fprintf(fp, "%f\t%f\t%f\t%f\n",diffp.x, un[j],u_analytic(&diffp),diff); //In datei ausgeben
        printf("%f\t%f\t%f\t%f\n",diffp.x, un[j],u_analytic(&diffp),diff);
        diffp.x=diffp.x+h; //neues x berechnen
	}
	fclose(fp);
}
//analytische Lösung fuer H10.1 speichern
void plotu_analytic(soli_param * sp){
	FILE *fp =fopen(dateiname, "w"); //Dateiöffnen
	double f; //Variable zum speichern des Wertes
	double x=sp->xmin;
    double xmax=sp->xmax;
    double xmin=x;
    soli_param ua= *sp; //neuer struct zum ändern der Parameter
    ua.x=x;
    ua.t=sp->t;
	while(ua.x<=xmax){
		f=u_analytic(&ua); //rechnen
		fprintf(fp, "%f\t%f\n", ua.x,f); //in Datei schreiben
		ua.x=ua.x+0.1;

	} 
}	

int main(void) {
    soli_param sp; //definieren der Variablen im struct
    sp.N = 1;
    sp.h=0.1;
    sp.d=pow(sp.h,3)/10;
    sp.tmax=1;
    sp.tmin=0;
    sp.xmin=-30;
    sp.xmax=30;
    sp.jmax=(sp.xmax-sp.xmin)/sp.h;
    sp.nmax=(sp.tmax-sp.tmin)/sp.d;
    //H10.1:
    /*sp.t=-1;
    sp.xmin=-20;
    sp.xmax=20;
	int i=0;
	while(sp.t<=1){
		sprintf(&dateiname[0],"%s%d","plotu1",i);
		plotu_analytic(&sp);
		sp.t+=0.1;
		i++;
	}*/
	//H10.2
	sp.N=2;
    sp.xmin=-30;
    sp.xmax=30;
    u0 = (double*) malloc(sp.jmax * sizeof(double)); //speicher allocieren
    ustart(&sp); //ustart aufrufen um u0 zu berechnen
    uj1 = (double*) malloc(sp.jmax * sizeof(double));
    ujeins(&sp); //ujeins aufrufen um uj1 zu berechnen
    nextu(&sp); //berechnet u für t=1
    free(u0); //speicher befreien
    free(uj1);
}