/* intégration du mouvement du pendule par Runge-Kutta 4
 * version minimaliste avec qqes modifications cosmetiques
 * compiler/executer comme suit:
 *   gcc -Wall pendule_rk4_moy.c -lm -o pendule_rk4_moy
 *   ./pendule_rk4_moy > pendule.data
 * 2008 Adrian Daerr, Univ Paris Diderot
 * domaine public
 */

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

double Dtheta(double theta, double omega) { return omega; }
double Domega (double theta, double omega) { return -sin(theta); }

int main(void) {
  double t=0.0, tfinal=50.0, dt=0.05;
  double theta=0.0, omega=1.99;// conditions initiales
  double theta_int, omega_int;
  double dtheta1, domega1;
  double dtheta2, domega2;
  double dtheta3, domega3;
  double dtheta4, domega4;
  int n=0;
  while (t<tfinal) {
    t = (n++)*dt;

    // afficher l'état du système
    printf("%f %f %f\n",t,theta,omega);

    // integrer d'un pas de temps
    dtheta1 = Dtheta(theta, omega);
    domega1 = Domega(theta, omega);

    theta_int = theta + 0.5*dt*dtheta1;
    omega_int = omega + 0.5*dt*domega1;

    dtheta2 = Dtheta(theta_int, omega_int);
    domega2 = Domega(theta_int, omega_int);

    theta_int = theta + 0.5*dt*dtheta2;
    omega_int = omega + 0.5*dt*domega2;

    // ce programme contient des traces d'Omega3 !
    dtheta3 = Dtheta(theta_int, omega_int);
    domega3 = Domega(theta_int, omega_int);

    theta_int = theta + dt*dtheta3;
    omega_int = omega + dt*domega3;

    dtheta4 = Dtheta(theta_int, omega_int);
    domega4 = Domega(theta_int, omega_int);

    theta = theta + dt*(dtheta1 + 2*dtheta2 + 2*dtheta3 + dtheta4)/6;
    omega = omega + dt*(domega1 + 2*domega2 + 2*domega3 + domega4)/6;
  }
  return EXIT_SUCCESS;
}