//etude du pendule mathematique avec la methode de Runge
//comparaison methode d'Euler et methode de Runge
#include <iostream>
#include <cmath>
#include "dislin.h"

using namespace std;

int main() {
//initialisation de DISLIN 
  metafl("XWIN");     //XWIN ou PDF
  disini();
  name("temps (s)", "X");
  name("angle (degrees)", "Y");
  graf(0.,100.,0.,10.,-15.,15.,-15.,3.);

//methode d'Euler    
  const int STEPS = 1000;
  double dt = 0.1;
  double angleE[STEPS], omegaE[STEPS], timeE[STEPS]; 
  //conditions initiales
  angleE[0] = 5.;
  omegaE[0] = 0.;
  timeE[0] = 0.;
  double acc;    
  for(int i=0; i<STEPS-1; i++) {
    angleE[i+1] = angleE[i] + omegaE[i] * dt;
    acc = - angleE[i];
    omegaE[i+1] = omegaE[i] + acc * dt;
    timeE[i+1] = timeE[i] + dt;
  }
//dessin de la trajectoire
  thkcrv(10);
  color("RED");
  curve(timeE, angleE, STEPS);

//solution exacte
  double angle[STEPS], omega[STEPS], time[STEPS];
  //conditions intiales 
  angle[0] = 5.*M_PI/180.;   //angle initial en radians
  omega[0] = 0.;   //vitesse angulaire initiale
  time[0] = 0.;
  //omega0 = sqrt(g/L) = 1
  //phi0 = 0
  double omega0 = 1.;
  for(int i=0; i<STEPS-1; i++) {
    time[i+1] = (i+1) * dt;
    angle[i+1] = angle[0] * cos(omega0*time[i+1])*180./M_PI;
  }
//dessin de la trajectoire
  color("GREEN");
  curve(time, angle, STEPS);
  	
//methode de Runge
  double angleR[STEPS], omegaR[STEPS], timeR[STEPS];
  //conditions initiales
  angleR[0] = 5.;
  omega0 = 0.;
  //vitesse initiale, l'acceleration est negative
  acc = -angleR[0];
  omegaR[0] = omega0 -acc*dt/2.;
  timeR[0] = 0.;        
  for (int i=0; i<STEPS-1; i++) {
    acc = - angleR[i];
    omegaR[i+1] = omegaR[i] + acc*dt;
    angleR[i+1] = angleR[i] + omegaR[i+1]*dt;
    timeR[i+1] = timeR[i] + dt;
  }
//dessin de la trajectoire     	 
  color("BLUE");
  curve(timeR, angleR, STEPS);
	
  disfin();     //fin de DISLIN
  
  return 0;
}