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

int mmax = 30;           //globale Variable, Maximale Anzahl der Schritte

double InfInt(int n, double a, double z,double eps,double (*func)(double, double, double));

double f1(double x, double a, double z) {   // Zu integrierende Funktion
    return (exp(-(x * x) / (a * a)) / sqrt(x * x + z * z));
}

double f2(double x, double a, double z){
    return((exp(-(x*x)/(a*a))*z)/pow((x*x)+(z*z),1.5));
}

double trapez(int n, double a, double z, double (*func)(double, double, double),
        double x1, double x2) {  //Trapeze
    double sum;
    double h = (x2 - x1) / (n);
    //Trapezsumme für Grenzen berechnen
    sum = (0.5 * func(x1, a, z)) + (0.5 * func(x2, a, z)); //in der Summe Funtion aufrufen um Werte an Grenzen x1 und x2 zu bekommen
    for (int i = 1; i < n; i++) {
        sum += func(x1 + (h * i), a, z);
    }
    return sum * h;
}

double IntRomb(double x1, double x2, int n, double a, double z,
        double (*func)(double, double, double), int eps) {
    int i, k;   //Variablen für Schleifen
    double ergebnis; //zum speichern des Ergebnis
    double **T; //Trapezsummen (später als 2-dim Array)
    double *h; //Schrittweiten Array

    h = (double*) malloc((mmax + 1) * sizeof(double)); //Speicherplatz für Array allozieren
    T = (double**) malloc(((mmax + 1)) * sizeof(double*));

    for (i = 0; i <= mmax; i++) {
        T[i] = (double*) malloc(((mmax + 1)) * sizeof(double)); //Speicherplatz für 2-dim Array allozieren

    }

    h[0] = (x2 - x1) / (n); //Erste Schrittweite berechnen

    T[0][0] = trapez(n, a, z, func, x1, x2); //erstes Trapez berechnen

    for (i = 1; i <= mmax; i++) {
        n = n * 2;           //Trapeze für k=0 berechnen, indem man n vergrößert
        T[i][0] = trapez(n, a, z, func, x1, x2);
        h[i] = h[i - 1] / 2.0;   //h berechnen für i>1 indem man durch 2 teilt

        for (k = 1; k <= i; k++) { //mit Formel aus VL T für k>1 berechnen

            T[i][k] = (-h[i - k] * h[i - k]
                    / ((h[i] * h[i]) - (h[i - k] * h[i - k])) * T[i][k - 1])
                    + (h[i] * h[i] / (h[i] * h[i] - h[i - k] * h[i - k])
                            * T[i - 1][k - 1]);
        }
        //printf("%d",i);
        if (abs(T[i][i] - T[i - 1][i - 1]) <= eps) { //Abschätzen ob das Ergebnis genau genug ist wenn ja aus Schleife raus
            ergebnis = T[i][i];

            break;
        }

    }
    if (i==mmax){
        printf("Fehler, Integral ist ungenau"); // Fehler falls Abbruchbedingung der vorherigen Schleife nicht greift
    }
    //Speicherplatz befreien
    for (i = 0; i <= mmax; i++) { //Speicherplatz in Schleife freigeben
        free(T[i]);
    }
    free(T);
    free(h);
    return ergebnis; //Rückgabewert ist Ergebnis
}

double potenzial(int n, double a, double z,
        double (*func)(double, double, double), double eps) {
    double V; //Variable Potential
    double alpha = 1/sqrt(M_PI); //normierung alpha
    double Q = 1; //Gesamtladung Q
    V = alpha * (Q / a) * InfInt(n, a, z, eps,func); //Potential berechnen
    return V;
}
    
double InfInt(int n, double a, double z, double eps,double (*func)(double, double, double)) {
    double x1 = -1; //Grenzen definieren
    double x2 = 1;
    double R=IntRomb(x1,x2,n,a,z,func,eps); //Variable zum speichern des Ergebnis, Integral für erste Grenzen ausrechnen
    double lastR; //Vergleichsvariable definieren
    do {
        lastR = R;//letztes Ergebnis speichern
        x1 = x1*2;  //Grenzen vergrößern
        x2 = x2*2;

        R = IntRomb(x1, x2, n, a, z, func, eps); //neues Ergebnis berechnen

    } while ((lastR/R<(1-eps)) || (lastR/R>(1+eps))); //Integralwerte vergleichen ,sodass deren Änderung kleiner als eps ist
    return R;
}

void efeld(double a, double z, int n, double eps,
        double *dfp,double (*func)(double, double, double)) {
    double f[2]; //Array zum speichern der Funktionswerte
    double h =0.5* z; //h relativ groß wird in Schleife verkleinert
    double lastdf;
    f[0] = potenzial( n, a, z +(h), func, eps); //f(z+h) berechnen
    f[1] = potenzial(n, a, z - h, func, eps); //f(z-h) berechnen
    *dfp = -1* (f[0] - f[1]) / (double)(2 * h);
    do {
        h=h/2; //h kleiner machen
        lastdf=*dfp; //Alte Ableitung auf neuer Variable speichern
        f[0] = potenzial( n, a, z +(h), f1, eps); //f(z+h) berechnen
        f[1] = potenzial(n, a, z - h, f1, eps); //f(z-h) berechnen


        *dfp = -1* (f[0] - f[1]) / (2 * h); //Differenzenquotient
    } while((lastdf/(*dfp)<(1-eps)) || (lastdf/(*dfp)>(1+eps)));  //Werte vergleichen ,sodass deren Änderung kleiner als eps ist
}


int main(void) {
    double z, a;
    int n=30;
    double eps = 1e-8; //genauigkeit epsilon
    	double df; //Variable zum speichern des Differenzenquotienten
    //Eingabe der Variablen:
    printf("Bitte geben Sie den Abstand der Punktladung von der x-Achse z ein\n ");
    scanf("%lf", &z);
    printf("Bitte geben Sie a ein\n");
    scanf("%lf", &a);



    printf("Unendliches Inegral %lf\n", InfInt(n, a, z, eps,f1)); //Integral berechnen mit InfInt
    printf("Das Potenzial ist %f\n", potenzial(n,a,z,f1,eps)); //Potenzial aus Integral berechnen
    //H2.1 Berechnung des Efeldes mit Hilfe numerischer Ableitung
    efeld(a,z,n,eps,&df,f1);        //E-feld berechnen
    printf("Das numerisch berechnete E-Feld an der Stelle %3.2lf ist:%f\n", z, df);
    //H2.2 Berechnung des Efeldes mit Hilfe analytischer Ableitung
    printf("Das anlytisch berechnete E-Feld an der Stelle %3.2lf ist:%f\n", z, potenzial(n,a,z,f2,eps));
    double diff=fabs(df - potenzial(n,a,z,f2,eps));
    printf("Differenz %1.10f:\n", diff);

}