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

using namespace std;

int main() {
  const int STEPS = 5000;
  double dt = 0.01;
    
//initialisation de DISLIN
  metafl("XWIN");     //XWIN ou PDF
  disini();
  name("temps (s)", "X");
  name("angle (degrees)", "Y");     
  graf(0.,50.,0.,10.,-180.,180.,-180.,30.);
  
//premiere iteration
  double angle[STEPS], omega[STEPS], time[STEPS]; 
  double omega0, acc;
  //conditions initiales
  angle[0] = 90.;   //angle initial
  omega0 = 0.;
  omega[0] = omega0 + angle[0]*dt/2.;     //vitesse initiale
  time[0] = 0.;
  for (int i=0; i<STEPS-1; i++) {
    acc = - sin(angle[i]/180.*M_PI)*180./M_PI;
    omega[i+1] = omega[i] + acc*dt;
    angle[i+1] = angle[i] + omega[i+1]*dt;
    time[i+1] = time[i] + dt;
  }
//dessin de la trajectoire
  thkcrv(10);
  color("RED");
  curve(time, angle, STEPS);

//deuxieme iteration
  double angle2[STEPS], omega2[STEPS], time2[STEPS];
  //conditions initiales
  angle2[0] = 135.;   //angle initial
  omega0 = 0.;
  omega2[0] = omega0 + angle2[0]*dt/2.;     //vitesse initiale
  time2[0] = 0.;  
  for (int i=0; i<STEPS-1; i++) {
    acc = - sin(angle2[i]/180.*M_PI)*180./M_PI;
    omega2[i+1] = omega2[i] + acc*dt;
    angle2[i+1] = angle2[i] + omega2[i+1]*dt;
    time2[i+1] =  time2[i] + dt;
  }
//dessin de la trajectoire  
  color("GREEN");
  curve(time2, angle2, STEPS);
	
//troisieme iteration
  double angle3[STEPS], omega3[STEPS], time3[STEPS];
  //conditions initiales
  angle3[0] = 175.;   //angle initial
  omega0 = 0.;
  omega3[0] = omega0 + angle3[0]*dt/2.;     //vitesse initiale
  time3[0] = 0.;  
  for (int i=0; i<STEPS-1; i++) {
    acc = - sin(angle3[i]/180.*M_PI)*180./M_PI;
    omega3[i+1] = omega3[i] + acc*dt;
    angle3[i+1] = angle3[i] + omega3[i+1]*dt;
    time3[i+1] = time3[i] + dt;
  }
//dessin de la trajectoire       	 
  color("BLUE");
  curve(time3, angle3, STEPS);
	
//quatrieme iteration
  double angle4[STEPS], omega4[STEPS], time4[STEPS];
  //conditions initiales
  angle4[0] = 30.;   //angle initial
  omega0 = 0.;
  omega4[0] = omega0 + angle4[0]*dt/2.;     //vitesse initiale
  time4[0] = 0.;  
  for(int i=0; i<STEPS-1; i++) {
    acc = - sin(angle4[i]/180.*M_PI)*180./M_PI;
    omega4[i+1] = omega4[i] + acc*dt;
    angle4[i+1] = angle4[i] + omega4[i+1]*dt;
    time4[i+1] = time4[i] + dt;
  }
//dessin de la trajectoire
  color("WHITE");
  curve(time4, angle4, STEPS);
  color("FORE");
  
  disfin();     //fin de DISLIN
 
  return 0;
}