/* Integrateur pendule, version modularisée */
/* les variables du modele sont stockées dans les tableaux selon
 * la correspondance:
 *   theta (angle) <-> S[0]
 *   omega (vitesse angulaire) <-> S[1]
 * compiler/executer comme suit:
 *   gcc -Wall pendule_rk4.c -lm -o pendule_rk4
 *   ./pendule_rk4 theta_init omega_init > pendule.data
 * 2008 Adrian Daerr, Univ Paris Diderot
 * domaine public
*/

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

// nombre de variables d'etat (=nombre d'equations d'evolution)
#define NEQ 2

/* conditions initiales de notre problème,
 si possible à partir des paramètres donnés en ligne de commande */
void init(int argc, char * argv[], double *S) {
  double val;
  int conv_ok;
  S[0]=0;// conditions initiales par défaut, utilisés si les paramètres sont
  S[1]=1.99;// absents ou illisibles
  if (argc<3) {
    fprintf(stderr,"usage: %s theta_init omega_init",argv[0]);
    return;
  }
  conv_ok = sscanf(argv[1],"%lf",&val);// convertir l'argument en un double
  if (conv_ok<1)
    fprintf(stderr,"erreur à la conversion du premier argument %s\n",argv[1]);
  else
    S[0] = val;
  conv_ok = sscanf(argv[2],"%lf",&val);
  if (conv_ok<1)
    fprintf(stderr,"erreur à la conversion du deuxième argument %s\n",argv[2]);
  else
    S[1] = val;
}

/* Cette fonction contient toute la physique: les équations
   d'évolution sont utilisées pour calculer le changement d'état dS à
   partir de l'état S */
void derive(double *S, double* dS) {
  dS[0] = S[1];
  dS[1] = -sin(S[0]);
}

/* Cette fonction intègre les équations d'évolution sur un pas de
   temps h, selon un schéma de Runge-Kutta d'ordre 4, en partant de
   l'état contenu dans le tableau 'etat'. Le nouvel état du système
   est finalement stocké dans le même tableau. */
void rk4(double *etat, double h) {
  double etatTest[NEQ], k1[NEQ], k2[NEQ], k3[NEQ], k4[NEQ];
  int i;

  // k1 = f(etat)
  derive(etat,k1);

  // k2 = f(etat+0.5*dt*k1)
  for (i=0; i<NEQ; i++)
    etatTest[i] = etat[i] + 0.5*h*k1[i];
  derive(etatTest,k2);

  // k3 = f(etat+0.5*dt*k2)
  for (i=0; i<NEQ; i++)
    etatTest[i] = etat[i] + 0.5*h*k2[i];
  derive(etatTest,k3);

  // k4 = f(etat+dt*k3)
  for (i=0; i<NEQ; i++)
    etatTest[i] = etat[i] + h*k3[i];
  derive(etatTest,k4);

  // nouvel etat = ancien etat + dt/6 * (k1 + 2*k2 + 2*k3 + k4)
  for (i=0; i<NEQ; i++)
    etat[i] = etat[i] + h/6*(k1[i] + 2*k2[i] + 2*k3[i] + k4[i]);
}

/* Point d'entree du programme. Cette fonction contient
   essentiellement la boucle d'intégration et l'affichage de l'état à
   chaque pas de temps */
int main(int argc, char * argv[]) {
  // tableau contenant l'etat du systeme (voir premieres lignes du
  // fichier pour la description des elements)
  double S[NEQ];
  // intervalle d'intégration [t,tmax] et pas d'intégration
  double t=0, tmax=100, dt=0.05;
  // compteur de pas
  int n=0;
  // initialisation
  init(argc,argv,S);
  // boucle d'intégration
  while (t < tmax) {
    t = (n++)*dt;
    // afficher l'état actuel
    printf("%lf %lf %lf\n",t,S[0],S[1]);
    // intégrer d'un pas de temps dt
    rk4(S,dt);
  }
  // affichage de l'état final
  t = n*dt;
  printf("%lf %lf %lf\n",t,S[0],S[1]);
  return EXIT_SUCCESS;
}