//pendule physique avec la methode de Runge
#include <iostream>
#include <cmath>
#include "dislin.h"

using namespace std;

int main() {
  const int STEPS = 50000;
  double dt = 0.002;
  double angle[STEPS], omega[STEPS], time[STEPS]; 
  angle[0] = 5.*M_PI/180.;   //angle initial en radians
  omega[0] = 0.;   //vitesse angulaire initiale
  time[0] = 0.;

  //constantes du pendule physique  
  double g = 9.81;
  double L = 40.;
  double gamma = 0.04;
  double force = 0.2;
  double freq = 2.;

  double accG, accAm, accEn, acc;
  //vitesse initiale pour la methode de Runge
  accG = -g/L * sin(angle[0]);
  accAm = -gamma * omega[0];
  accEn = force * sin(freq*time[0]);
  acc = accG + accAm + accEn;
  omega[0] = omega[0] - dt/2. * acc;
  for (int i=1; i<STEPS; i++) {
    accG = - g/L * sin(angle[i-1]);
    accAm = - gamma * omega[i-1];
    accEn = force * sin(freq*time[i-1]);
    acc = accG + accAm + accEn;
    omega[i] = omega[i-1] + acc*dt;
    angle[i] = angle[i-1] + omega[i]*dt;
    time[i] = i*dt;
  }
  for (int i=0; i<STEPS; i++)     //conversion de l'angle en degrees
    angle[i] *= 180./M_PI;
    
  //initialisation de dislin
  metafl("XWIN");   // XWIN ou PDF
  disini();
  name("temps (s)", "X");
  name("angle (degrees)", "Y");
  graf(0.,100.,0.,10.,-20.,20.,-20.,5.);
  //desin de la trajectoire
  thkcrv(10);
  color("RED");
  curve(time, angle, STEPS);
  //fin de DISLIN
  disfin();
 
  return 0;
}