//*******************************************************************************
//*
//* Programme d'essai de la méthode d'Euler sur l'équation différentielle
//* du pendule simple.
//* 
//* Rappel : le système différentiel traité est d(theta)2/dt2 = -(g/l)*sin(theta)
//*
//* Dominique Lefebvre - TangenteX.com - Avril 2006
//*
//* NOTA : ce programme est a but didactique. Outre les limitations bien connues
//*        de la méthode d'Euler, la programmation n'est pas tres recherchee. En
//*        particulier, je ne ferais un usage tres modere des pointeurs afin de  
//*        ne pas compliquer les choses
//*******************************************************************************

//*******************************************************************************
//* Inclusion des headers standards
//*******************************************************************************  
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <malloc.h>

//*******************************************************************************
//*  Declaration des constantes
//*******************************************************************************
#define N 10					// Nombre de cycles de calcul
#define P 1000				// Nombre de pas de calcul par cycle
#define G 9.81				// Acceleration de la pesanteur
#define PI 3.141592654


//*******************************************************************************
//* Déclaration des variables globales.
//* La variable l est donc partagée par toutes les routines du programme
//*******************************************************************************

double l;        // Longueur du pendule (en m)


//*******************************************************************************
//*  Routine de description du système différentiel a integrer
//*  Ici, en l'occurence, il s'agit du système du pendule simple
//*******************************************************************************
void Derivee(double X[], double DX[])
{
  DX[0] = X[1];
  DX[1] = -G*sin(X[0])/l;
}

//*******************************************************************************
//*  Routine d'implémentation de la méthode d'Euler du premier ordre
//*******************************************************************************
void Euler(double x[], double y[], double DX[], double dt, int i)
{
   x[i+1] = x[i] + DX[0]*dt;
   y[i+1] = y[i] + DX[1]*dt;
}


//*******************************************************************************
//*   Corps du programme principal
//*******************************************************************************
void main()
{

 //* Declaration des variables
 int i;
 double theta,h, X[2], DX[2];
 double *x, *y, *t;
 FILE  *fp;

 
 //* Allocation de la mémoire pour les tableaux de calcul
 x = (double *)malloc((unsigned)(N*P+1)*sizeof(double));
 y = (double *)malloc((unsigned)(N*P+1)*sizeof(double));
 t = (double *)malloc((unsigned)(N*P+1)*sizeof(double));

 //* Initialisation des constantes
 l = G/(4*PI*PI);					// j'utilise une longueur qui m'arrange
 theta = 2*PI*sqrt(l/G);            // periode du pendule
 h = theta/P;						// pas temporel pour le calcul

 //* Determination des conditions initiales
 x[0] = 60*PI/180;					// angle initial du pendule de 60° converti en radians
 y[0] = 0.0;						// vitesse initiale du pendule nulle
 
 //* Calcul
 for (i=0; i<N*P; i++)
 {
   t[i] = i*h;									// incrementation du temps
   
   // Calcul des derivees, qui sont retournees dans DX
   X[0] = x[i];
   X[1] = y[i];
   Derivee(X,DX);
   
   // Application de la méthode d'Euler, qui calcule x(i+1) et y(i+1) en fonction
   // de x(i), y(i), dx(i) et dy(i)
   Euler(x, y, DX, h, i);
 }

 // Sauvegarde des résultats dans un fichier texte pour tracer par GnuPlot
 // Ce fichier se nomme EulerC.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("EulerC.dat","w+");		
   for(i=0;i<N*P; i++)
     fprintf(fp,"%f %f %f\n", t[i], x[i]*180./PI, y[i]);
 fclose(fp);

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