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

using namespace std;

int main() {
  const int steps = 50000;
  const double GM = 1.;   //G * M
  
  double xpos[steps], ypos[steps], xvel[steps], yvel[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] = 0.5;

  double dT = 0.002;   //pas d'integration
    
  double normR, xacc, yacc;

  //methode de Runge
  //velocite initiale dans t0 - dt/2
  normR = sqrt(xpos[0]*xpos[0] + ypos[0]*ypos[0]);
  xacc = - GM * xpos[0] / (normR*normR*normR);   //force ou acceleration, GM = 1
  yacc = - GM * ypos[0] / (normR*normR*normR);
  xvel[0] = xvel[0] - xacc*dT/2.;   //vitesse 
  yvel[0] = yvel[0] - yacc*dT/2.;
  
  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);
    
    xvel[i] = xvel[i-1] + xacc*dT;   //vitesse 
    yvel[i] = yvel[i-1] + yacc*dT;
    
    xpos[i] = xpos[i-1] + xvel[i]*dT;   //position
    ypos[i] = ypos[i-1] + yvel[i]*dT;
  }
  
  //partie graphique
  metafl("XWIN");   // XWIN ou PDF
  page(3000,3000);
  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();
 
  return 0;
}