//etude d'une orbite avec la methode de Runge
//la loi gravitationelle n'est pas necessairement 1/r^2
#include <iostream>
#include <cmath>
#include "dislin.h"

using namespace std;

int main() {
  const int steps = 50000;
  const double GM = 1.;
  double xpos[steps], ypos[steps], xvel[steps], yvel[steps];
  //conditions initiales
  xpos[0] = 2.0; ypos[0] = 0.0;
  xvel[0] = 0.0; yvel[0] = 0.5;

  //nouvelle loi gravitationelle  
  double beta = 2.1;   //exposant loi gravitationnelle
  
  double dT = 0.002;     //pas d'integration

  //methode de Runge
  //vitesse initiale en t0 - dt/2
  double normR, xacc, yacc;
  normR = sqrt(xpos[0]*xpos[0] + ypos[0]*ypos[0]);
  xacc = - GM / pow(normR, beta) * xpos[0] / normR;   //force ou acceleration, GM = 1
  yacc = - GM / pow(normR, beta) * ypos[0] / 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 / pow(normR, beta) * xpos[i-1] / normR;   //force ou acceleration, GM = 1
    yacc = - GM / pow(normR, beta) * ypos[i-1] / 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 avec DISLIN
  metafl("XWIN");   //XWIN ou PDF
  page(3000,3000);   //pour dessiner une image proportionee
  disini();     //initialisation de dislin
  name("axe-X","X");
  name("axe-Y","Y");
  graf(-2.,2.,-2.,1.,-2.,2.,-2.,1.);
  thkcrv(5);
  color("red");
  curve(xpos, ypos, steps);
  disfin();     //fin de DISLIN

  return 0;
}