#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
double force(double x, double k, double m, double alpha);
void passo(double *v, double *x, double m, double k, double alpha, double t);
int main() {
    double xmax, k, m, alpha;
    k = 9.81;
    m = 1.0;
    alpha = -1/6;
    double Thar = 2*M_PI*sqrt(m/k);
    printf("Il periodo valido per piccole oscillazioni e': %.5lf\n",Thar);
    double dt = Thar/pow(2,20);
    printf("Il passo temporale vale: %e\n",dt);
    double v, x, t, Ts[3], Test, err[30][2];
    int n, zeri;
    for (n = 0; n < 30; n++) {
        xmax = M_PI*(n+1)/180;
        err[n][0] = xmax;
        v = 0;
        t = 0;
        zeri = 0;
        x = xmax;
        do {
            passo(&v,&x,m,k,alpha,dt);
            t += dt;
            if (xmax*x < 0) {
                Ts[zeri++] = t;
            }
            xmax = x;
        } while (zeri < 3);
        Test = Ts[2] - Ts[0];
        err[n][1] = 100.0*fabs(Thar - Test)/(Thar);
    }
    FILE *f;
    f = fopen("differenza.dat","w");
    for (n = 0; n < 30; n++) {
        fprintf(f,"%lf %lf\n",err[n][0], err[n][1]);
    }
    fclose(f);
}

double force(double x, double k, double m, double alpha) {
    return (-k/m)*(x + alpha*x*x*x);
}

void passo(double *v, double *x, double m, double k, double alpha, double t) {
    double tempx = *x;
    double tempv = *v;
    tempv += force(tempx,k,m,alpha)*t;
    tempx += tempv*t;
    *x = tempx;
    *v = tempv;
}