#include <math.h>
#include <stdio.h>
// Parametri della simulazione. Unità di misura SI
const double
h0= 2000, // Quota iniziale
g= 9.8, // Accelerazione di gravità
rho= 1, // Densità dell'aria
m= 80, // Massa
S= .7, // Area
Cr= 1, // Coeff. di resistenza
k= .5 * rho * S * Cr,
dt= 1; // Passo d'integrazione
/*
Modificare il numero in base al tipo d'integratore che si vuole usare:
1: Eulero
2: Eulero semi-implicito
3: Runge-Kutta 2° ordine
4: Runge-Kutta 4° ordine
*/
#define TIPO_INTEGRATORE 1
// ============================================================
// === STATO ===
struct STS
{
double p, v;
};
STS operator*(const double mul, const STS &s)
{
STS r; r.p= mul * s.p; r.v= mul * s.v;
return r;
}
STS operator+(const STS &s1, const STS &s2)
{
STS r; r.p= s1.p + s2.p; r.v= s1.v + s2.v;
return r;
}
// ============================================================
int c_DER= 0; // Contatore chiamate a DER()
// Calcolo delle derivate
STS DER(const STS &st)
{
// Forza aerodinamica
double Fa= k * st.v * st.v; // Fa = k * v^2
// Forza dovuta all'accelerazione di gravità
double Fg= m * g;
STS der; c_DER++;
der.p= st.v; // Derivata della posizione = velocità
der.v= (Fa - Fg) / m; // Derivata della velocità = accelerazione
return der; // Ritorna p' e v'
}
// ============================================================
int main(void)
{
STS st; // STATO INIZIALE
st.p= h0; // Quota
st.v= 0; // Velocità
const double k1= sqrt(g*m/k), // Questa è la velocità limite
k2= sqrt(g*k/m),
tc= acosh(exp(k2/k1 * h0)) / k2; // Tempo di caduta
double emPOS= 0, emVEL= 0; // Errori massimi
int c_stp= 0; // Numero di passi d'integrazione
for(double t= 0; ; t += dt) {
// Valori ESATTI della velocità e della quota
double v= -k1 * tanh(t * k2);
double p= h0 - m/k*log(cosh(t * k2));
double ePOS= st.p - p; emPOS= fmax(emPOS, fabs(ePOS));
double eVEL= st.v - v; emVEL= fmax(emVEL, fabs(eVEL));
STS der = DER(st); // Calcolo delle derivate di p e di v
printf("%3d) t= %7.3f", ++c_stp, t);
printf(" h= %8.3f (%8.1e)", st.p, ePOS);
printf(" v= %10.6f (%8.1e)", st.v, eVEL);
printf(" a= %f\n", der.v);
if(st.p < 0) break; // Interrompe quando il corpo arriva al suolo
// Integrazione numerica per calcolare p e v a t + dt
#if TIPO_INTEGRATORE == 1 // Eulero
#define INTEGRATORE "Eulero"
st = st + dt * der; // Equivale a: p = p + dt * p', v = v + dt * v'
#elif TIPO_INTEGRATORE == 2 // Eulero semi-implicito
#define INTEGRATORE "Eulero semi-implicito"
st.v = st.v + dt * der.v;
st.p = st.p + dt * st.v;
#elif TIPO_INTEGRATORE == 3 // RK2
#define INTEGRATORE "RK-2"
STS f2= DER(st + .5 * dt * der);
st = st + dt * f2;
#elif TIPO_INTEGRATORE == 4 // RK4
#define INTEGRATORE "RK-4"
STS f1= der;
STS f2= DER(st + .5 * dt * f1);
STS f3= DER(st + .5 * dt * f2);
STS f4= DER(st + dt * f3 );
st= st + (dt/6) * (f1 + 2 * (f2 + f3) + f4);
#else
#error Il valore per TIPO_INTEGRATORE può andare da 1 a 4.
#endif
}
printf("\nIntegratore usato: %s\n", INTEGRATORE);
printf("\nErrori massimi\n");
printf(" POS: %g m, VEL = %g m/s\n", emPOS, emVEL);
printf("\nTempo di caduta: %.3f s, ", tc);
printf("velocita' limite: %.3f m/s\n", k1);
printf("\nIterazioni: %d, chiamate a DER: %d\n", c_stp, c_DER);
return 0;
}
