|
#include <stdlib.h>
|
|
#include <stdio.h>
|
|
#include <time.h>
|
|
#include <math.h>
|
|
|
|
#define KAPPA 1.3
|
|
#define MASSA 0.9
|
|
|
|
double energy(double x, double v) {
|
|
return 0.5 * MASSA * v * v + 0.5 * KAPPA * x * x;
|
|
}
|
|
|
|
void updSemiEuler(double *x, double *v, double dt) {
|
|
*v += dt * (-KAPPA * (*x)) / MASSA;
|
|
*x += dt * (*v);
|
|
}
|
|
|
|
void analytic(double t, double x0, double v0, double* x, double* v) {
|
|
double w;
|
|
w = sqrt(KAPPA / MASSA);
|
|
*x = x0 * cos(w * t) + (v0 / w) * sin(w * t);
|
|
*v = -w * x0 * sin(w * t) + v0 * cos(w * t);
|
|
}
|
|
|
|
void input(double *dt) {
|
|
do {
|
|
printf("Immetti il valore di dt, che deve essere compreso tra 0.001 e 0.01\n");
|
|
scanf("%lf", dt);
|
|
if(*dt < 0.001 || *dt > 0.01) {
|
|
printf("Valore fuori dall'intervallo [0.001,0.01], riprova!\n");
|
|
}
|
|
} while(*dt < 0.001 || *dt > 0.01);
|
|
}
|
|
|
|
int main() {
|
|
srand48(time(NULL));
|
|
int i, nsteps;
|
|
double VMIN = 0.5, VMAX = 1.5, x0 = 1.0;
|
|
double dt, xA, vA, eneA, eneSE;
|
|
double v0 = VMIN + drand48() * (VMAX - VMIN);
|
|
double xSE = x0;
|
|
double vSE = v0;
|
|
input(&dt);
|
|
nsteps = 2.0 * M_PI * sqrt(MASSA / KAPPA) / dt;
|
|
|
|
FILE *out_traj = fopen("traiettoria.dat", "w");
|
|
if(out_traj == NULL) {
|
|
printf("ERRORE NELL'APERTURA FILE\n");
|
|
exit(0);
|
|
}
|
|
|
|
FILE *out_E = fopen("energia.dat", "w");
|
|
if(out_E == NULL) {
|
|
printf("ERRORE NELL'APERTURA FILE\n");
|
|
exit(0);
|
|
}
|
|
|
|
for(i = 0; i < nsteps; i++) {
|
|
updSemiEuler(&xSE, &vSE, dt);
|
|
if(i % 10 == 0) {
|
|
analytic(dt * i, x0, v0, &xA, &vA);
|
|
eneSE = energy(xSE, vSE);
|
|
eneA = energy(xA, vA);
|
|
fprintf(out_traj, "%.5f %.6f %.6f\n", dt * i, xA, xSE);
|
|
fprintf(out_E, "%.5f %.6f %.6f\n", dt * i, eneA, eneSE);
|
|
}
|
|
}
|
|
|
|
fclose(out_traj);
|
|
fclose(out_E);
|
|
|
|
return 0;
|
|
}
|