//etude d'une orbite avec la methode d'Euler
#include <iostream>
#include <cmath>
#include "dislin.h"

using namespace std;

int main() {
  const int steps = 10000; 
  const double GM = 1.;   //G * M
  double xpos[steps], ypos[steps], xvel[steps], yvel[steps], time[steps];
  //conditions initiales pour orbite circulaire
  //xpos[0] = 1.0; ypos[0] = 0.0;
  //xvel[0] = 0.0; yvel[0] = 1.0;
  //conditions initiales pour orbite elliptique
  //si V0 < 1 -> ellipse
  //si V0 > 1 -> hyperbole
  xpos[0] = 1.0; ypos[0] = 0.0;
  xvel[0] = 0.0; yvel[0] = 1.;
  time[0] = 0.;
    
  double dT = 0.005;   //pas d'integration
    
  double normR, xacc, yacc;

//methode d'Euler    
  for (int i=1; i<steps; i++) {
    normR = sqrt(xpos[i-1]*xpos[i-1] + ypos[i-1]*ypos[i-1]);
    xacc = - GM * xpos[i-1] / (normR*normR*normR);  //force ou acceleration, GM = 1
    yacc = - GM * ypos[i-1] / (normR*normR*normR);
    xpos[i] = xpos[i-1] + xvel[i-1]*dT;   //position
    ypos[i] = ypos[i-1] + yvel[i-1]*dT;
    xvel[i] = xvel[i-1] + xacc*dT;     //vitesse 
    yvel[i] = yvel[i-1] + yacc*dT;
    time[i] = i*dT;
  }

  metafl("XWIN");   //XWIN ou PDF
  page(3000, 3000);   //pour creer une image proportionee
  disini();     //initialisation de dislin
  name("axe-X", "X");
  name("axe-Y", "Y"); 
  graf(-2.,2.,-2.0,1.,-2.,2.,-2.,1.);
  thkcrv(5);
  color("red");
  curve(xpos, ypos, steps);

  disfin();   //fin de DISLIN

  return 0;
}