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

using namespace std;

int main() {
//solution exacte (premiere iteration)
  double angle[10000], omega[10000], time[10000];
  //conditions intiales 
  angle[0] = 5.*M_PI/180.;   //angle initial en radians
  omega[0] = 0.;   //vitesse angulaire initiale
  time[0] = 0.;

//initialisation de DISLIN
  metafl("XWIN");   //XWIN ou PDF
  disini();
  name("temps (s)", "X");
  name("angle (degrees)", "Y");
  graf(0.,50.,0.,10.,-25.,25.,-25.,5.);

//solution exacte
  int steps = 500;
  double dt = 0.1;
  //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
  thkcrv(10);
  color("RED");
  curve(time, angle, steps);

//premiere iteration
  double angle2[10000], omega2[10000], time2[10000]; 
  angle2[0] = 5.;   //angle initial 
  omega2[0] = 0.;   //vitesse angulaire initiale
  time2[0] = 0.;
  steps = 500;
  dt = 0.1;    
  for(int i=0; i<steps-1; i++) {
    angle2[i+1] = angle2[i] + omega2[i]*dt;
    omega2[i+1] = omega2[i] - angle2[i]*dt;
    time2[i+1] = time2[i] + dt;
  }
//dessin de la trajectoire	 
  color("GREEN");
  curve(time2, angle2, steps);
	
//deuxieme iteration
  double angle3[10000], omega3[10000], time3[10000]; 
  angle3[0] = 5.;   //angle initial 
  omega3[0] = 0.;   //vitesse angulaire initiale
  time3[0] = 0.;    
  steps = 1000;
  dt = 0.05;    
  for(int i=0; i<steps-1; i++) {
    angle3[i+1] = angle3[i] + omega3[i]*dt;
    omega3[i+1] = omega3[i] - angle3[i]*dt;
    time3[i+1] = time3[i] + dt;
  }
//dessin de la trajectoire	 
  color("BLUE");
  curve(time3, angle3, steps);
	
  disfin();     //fin de DISLIN
 
  return 0;
}