//******************************************************************************
//* Programme de  d'utilisation de la méthode d'Euler pour résoudre
//* une EDO de premier ordre.
//* Calcul de la chute d'un grelon
//*
//* Dominique Lefebvre - TangenteX.com octobre 2008
//******************************************************************************


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

//******************************************************************************
//* Calcul de la valeur de la dérivée en x
//* l'équation différentielle est y' = k*x²
//******************************************************************************
double Derivee(double x)
{
 double dx,g,k;
 
    //* definition desconstantes
    g = 9.81;
    k = 1.56e-2;
 
    dx = g-k*(x*x);
    return dx;
}
//******************************************************************************
//* Programme principal
//******************************************************************************
int main(int argc, char *argv[])
{
 
 //* Declaration des variables
 int i, N;
 double h;
 double *x, *y;
 FILE  *fp;

 //* Déclaration des constantes de calcul
  h = 0.1;						// pas pour le calcul
  N = 1000; 
 
 //* Allocation de la mémoire pour les tableaux de calcul
 x = (double *)malloc((unsigned)(N+1)*sizeof(double));
 y = (double *)malloc((unsigned)(N+1)*sizeof(double));
 
 
 //* Determination des conditions initiales
 x[0] = 0.0;                      // abscisse initiale                      
 y[0] = 0.0;                      // ordonnée initiale						
 
 //* Calcul
 for (i=0; i<N; i++)
 {
     y[i+1] = y[i] + Derivee(y[i])*h;         // on applique le schéma d'Euler
     x[i+1] = x[i] + h;
 }

 // Sauvegarde des résultats dans un fichier texte pour tracer par GnuPlot
 // Ce fichier se nomme grelon.dat. Il est cree dans le repertoire courant 
 // (celui dans lequel est lance le programme) 
 // S'il n'existe pas, il est cree. S'il existe, son contenu est ecrase.
 fp=fopen("grelon.dat","w+");		
   for(i=0;i<N; i++)
     fprintf(fp,"%f %f\n",x[i], y[i]);
 fclose(fp);

 //* Liberation de la memoire et sortie
 free(x);
 free(y);

 system("PAUSE");	
 return 0;
}