Quando colpirà la freccia il bersaglio?

/*
* Descrizione : Distanza di caduta ed tempo di volo
* di un freccia lanciata contro un bersaglio in aria.
* Attrito viscoso.
*
* Input : modulo velocita' iniziale
* angolo di lancio in gradi
* coefficiente di attrito
* distanza e dimensioni del bersaglio
*
* Output : Scrive sul file traiettoria.dat t, x(t) e y(t) della traiettoria
* della freccia;
* Stampa sullo standard output:
* (a) la quota massima y_max nonche' t_max e x_max ad essa corrispondenti
* (b) il tempo di impatto sul bersaglio calcolato con il metodo della bisezione
* (c) l'intervallo dei valori della velocita' iniziale per cui il bersaglio
* viene colpito.
*/
#include
#include
#include

/* Macro */

#define ACC_G 9.81 /* acc. gravita' in m/s^2 */
#define GR_RD 0.0174532925199433 /* fattore di conversione gradi -> rad */
#define DELTA_T 0.001 /* risoluzione temporale traiettoria */
#define T_2 100.0 /* tempo iniziale t_2 per la bisezione */
#define SIZE 10000 /* dimensione massima vettori per la traiettoria */
#define GAMMA 0.5 /* coefficiente di attrito / massa */
#define L 26.0 /* dimensione del bersaglio */

/* Prototipi delle funzioni */

double Xt(double time, double V[2]);
double Yt(double time, double V[2]);
double FindImpact(double t_1, double t_2, double V[2]);

/* Dichiarazione stream */

FILE *fp_out;

int main(void)
{
register int i;
int N_dati;
double h = 2.0; /* altezza centro del bersaglio */
double r = 0.5; /* raggio del bersaglio */
double v_0=24.8, alpha = 45.0; /* velocita iniziale */
double t_1, t_2, t_star;
double t_max, x_max, y_max;
double x_coord[SIZE];
double y_coord[SIZE];
double V[2];
/* char line[41];*/

/* Modulo velocita', angolo di lancio e attrito */
/*printf("\n");
printf("Modulo velocita' (m/s): ");
fgets(line, sizeof(line), stdin);
sscanf(line, "%lf", &v_0);*/

/* printf("Inserisci l'angolo (gradi) : ");
fgets(line, sizeof(line), stdin);
sscanf(line, "%lf", &alpha);*/

/* Componenti velocita' */
alpha *= GR_RD; /* Gradi -> Radianti */
V[0] = v_0 * cos(alpha);
V[1] = v_0 * sin(alpha);

if (V[1] <= 0) {
printf("\n Lancio verso il basso!\n\n");
return 1;
}


/***************************************************************************
* Stampa su file dei dati t. x(t), y(t) della triettoria della freccia
* con risoluzione temporale DELTA_T. Riempimento dei vettori x_coord e y_coord
***************************************************************************/

if(!(fp_out=fopen("traiettoria.dat", "w"))){
fprintf(stderr,"il file da aprire non esiste nella destinazione specificata\n");
exit(1);
}

i=0;
while(1) {
if(Xt(DELTA_T*i, V) < L) {
fprintf(fp_out,"%f %f %f\n",DELTA_T*i, Xt(DELTA_T*i, V),Yt(DELTA_T*i, V));
/* riempimento dei vettori con i dati sulla triettoria */
x_coord[i]=Xt(DELTA_T*i, V);
y_coord[i]=Yt(DELTA_T*i, V);
}
else break;
++i;
}

N_dati = i; /* numero di punti della traiettoria */

fclose(fp_out);

/***************************************************************************
* Determinazione dell' intervallo iniziale [t_1, t_2] per la bisezione:
* t_1 < t^* < t_2
***************************************************************************/

/* Sappiamo che al tempo t=0 x(0) - L e' negativo. Dobbiamo dunque identificare un tempo
in cui x(t_2)- L sia positivo. Si fissa un t_2 iniziale pari ad un valore scelto
arbitrariamente T_2. Se verifica x(t_2) - L e' ancora minore di zero si raddoppia
t_2 finche' non sia verificata la condizione x(t_2) - L > 0. */

t_2 = T_2;
while ( (Xt(t_2, V) - L) < 0.0 ) t_2 *= 2.0;

/* Una volta identificato t_2 si cerca di fissare un tempo t_1 per cui
x(t_1) - L < 0. A tal fine si fissa progressivamente t_1 alla meta' di t_2.
Finche' x(t_1) -L e' ancora positivo si ridefinisce t_2 = t_1 e si dimezza t_1.
Alla fine della procedura si hanno due valori t_1 e t_2 tali che
* t_1 < t^* < t_2 */

t_1 = 0.5 * t_2;
while ( (Xt(t_1, V) - L) > 0.0 ) {
t_2 = t_1;
t_1 = 0.5 * t_2;
}

/***************************************************************************
* Determinazione del tempo di impatto: t^*: x(t^*) - L = 0
* con t_1 < t^* < t_2. I dati sono forniti in maniera tale che
* la freccia raggiunga il bersaglio. In generale si dovrebbe verificare
* che il tempo di impatto corrisponda a delle coordinate x(t^*) e y(t^*)
* all'interno del bersaglio. L'algoritmo e' contenuto nella funzione
* FindImpact.
***************************************************************************/

t_star = FindImpact(t_1, t_2, V);

printf("\n************************************************************************\n\n");
if(Yt(t_star, V) <= h +r && Yt(t_star, V) >= h-r )
printf("La freccia ha colpito il bersaglio alla quota di %.4f metri\n\n", Yt(t_star, V));
printf("Tempo di volo: %.8f secondi \n\n", t_star);
printf("************************************************************************\n");


/***************************************************************************
* Determinazione dell'intervallo di velocita' per cui la freccia raggiunge il
* bersaglio.
***************************************************************************/

return 0;
}
/* Fine main */

/* Funzioni */

/* Funzione coordinata x */
double Xt(double t, double V[2])
{
double coef;
coef = V[0] / GAMMA;

return coef * (1.0 - exp(-GAMMA * t));
}

/* Funzione coordinata y */
double Yt(double t, double V[2])
{
double coef;
double y;

coef = 1.0 / GAMMA;
y = coef * (1.0 - exp(-GAMMA *t)) * (V[1]+ ACC_G / GAMMA) - ACC_G / GAMMA * t;
return y;
}



/* Funzione che calcola il tempo di impatto
* FindImpact()
*
* tempo di impatto con il metodo della bisezione
*
* MAX_BIS: Numero massimo bisezioni
* EPSILON: Precisione richiesta
*/
#define MAX_BIS 40
#define EPSILON (1.e-8)
double FindImpact(double t_1, double t_2, double V[2])
{
int i;
double t_root, t_mid;
double diff_mid, dt;

dt = t_2 - t_1; /* intervallo iniziale */
/* t_root = t_2; */ /* assunzione iniziale */

for (i = 0; i < MAX_BIS; ++i) {
dt *= 0.5;
t_mid = t_2 - dt;
diff_mid = Xt(t_mid, V)-L;
if (diff_mid >= 0.0)
t_2 = t_mid;
else
t_1 = t_mid;


if (dt < EPSILON || diff_mid == 0.00) return (t_1+t_2)*0.5;
}

printf("\n Error!!! Precisione richiesta eccessiva \n");
exit (1);
}
#undef MAX_BIS
#undef EPSILON