Numerical Analysis 2

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

double k, r_a, z_a, phi_a, R, Q;

double EZ(double integral);
double EPolar(double integral);
double integralZ(double phi);
double integralPolar(double phi);
double trapezoid(double a, double b, double n, double (*func)(double));
double simpson(double a, double b, double n, double (*func)(double));
void results(double, double, double (*field)(double), double (*integral)(double), double (*method)(double,  double,  double,  double (*)(double)));

int main(){
    k = 9 * pow(10, 9);
    Q = pow(10, -6);
    r_a = 0.1;
    z_a = 0.4;
    R = 0.05;
    phi_a = 35;

    double (*E[2])(double) = {EZ, EPolar};
    double (*I[2])(double) = {integralZ, integralPolar};
    double (*trap)(double, double, double, double (*)(double)) = trapezoid;
    double (*simp)(double, double, double, double (*)(double)) = simpson;

    for (int problem = 0; problem < 2; problem++){
        if (problem == 0){
            printf("Ερώτημα 1:\Καρτεσιανή συνιστώσα της έντασης ηλεκτρικού πεδίου.\n");
        } else{
            printf("Ερώτημα 2:\tΠολική συνιστώσα της έντασης ηλεκτρικού πεδίου.\n");
        }
        printf("Υπολογισμός έντασης πεδίου με τη μέθοδο του τραπεζίου.\n");
        results(0, 2*M_PI, *E[problem], *I[problem], *trap);
        printf("Υπολογισμός έντασης πεδίου με τη μέθοδο του Simpson.\n");
        results(0, 2*M_PI, *E[problem], *I[problem], *simp);
    }

        printf("All ok.");
    return 0;
}


double EZ(double integral){
    return k * Q * z_a * integral / (2 * M_PI);
}

double integralZ(double phi){
    double a = sqrt(pow(r_a, 2) + pow(R, 2) + pow(z_a, 2));
    return (1 / pow((pow(a, 2) - 2 * r_a * R * cos(phi - phi_a)), 1.5));
}

double EPolar(double integral){
    return k * Q * integral / (2 * M_PI);
}

double integralPolar(double phi){
    double a = sqrt(pow(r_a, 2) + pow(R, 2) + pow(z_a, 2));
    return (sqrt(pow(r_a, 2) + pow(R, 2) - 2 * r_a * R * cos(phi - phi_a)) / pow((pow(a, 2) - 2 * r_a * R * cos(phi - phi_a)), 1.5));
}

double trapezoid(double a, double b, double n, double (*func)(double)){
    double h, s, x_i;
    h = (b - a) / n;
    s = 0.0;
    for (int i=0; i<=n; i++){
        x_i = a + i * h;
        if (i ==0 || i==n){
            s += func(x_i)/2;
        } else{
            s += func(x_i);
        }
    }
    return s*h;
}

double simpson(double a, double b, double n, double (*func)(double)){
    double h, s, x_0, x_i, x_n;
    h = (b - a) / n;
    s = 0.0;
    for (int i=0; i<=n; i++){
        x_i = a + i * h;
        if (i==0 || i==n){
            s += func(x_i);
        } else if (i%2==0){
            s += 2*func(x_i);
        } else{
            s += 4*func(x_i);
        }
    }
    return h*s/3;
}

void results(double a, double b, double (*field)(double), double (*integral)(double), double (*method)(double,  double,  double,  double (*)(double))){
    double integral_0, integral_1, error = 1.0;
    int n = 2;
    for (int i = 0; error >= 0.5*pow(10,-8); i++){
        integral_0 = method(a, b, n, integral);
        integral_1 = method(a, b, n *= 2, integral);
        error = fabs(integral_0 - integral_1);
    }
    printf("FIELD: %.8f kNC^-1\n", field(integral_0)/1000);
}