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

using namespace std;

int main() {
  const int steps = 20000; 
  const double GM = 1.;     //G_N * M_Soleil
  const double m = 1.;      //masse planete
  double xpos[steps], ypos[steps], xvel[steps], yvel[steps], time[steps]; 
  double ePot[steps], eCin[steps], eTot[steps], l[steps];
  //conditions initiales
  xpos[0] = 1.0; ypos[0] = 0.0;
  xvel[0] = 0.0; yvel[0] = 1.3;
  time[0] = 0.;
    
  double dT = 0.01;   //pas d'integration

  //methode de Runge
  double dist, vel, xacc, yacc;     //varialbes auxiliares
  //redefinition de la vitesse initiale
  dist = sqrt(xpos[0]*xpos[0] + ypos[0]*ypos[0]);
  vel = sqrt(xvel[0]*xvel[0] + yvel[0]*yvel[0]);
  xacc = -GM * xpos[0] / (dist*dist*dist);  //force ou acceleration, GM = 1
  yacc = -GM * ypos[0] / (dist*dist*dist);
  xvel[0] = xvel[0] - xacc*dT/2.;     //vitesse 
  yvel[0] = yvel[0] - yacc*dT/2.;
  for (int i=1; i<steps; i++) {
    dist = sqrt(xpos[i-1]*xpos[i-1] + ypos[i-1]*ypos[i-1]);
    vel = sqrt(xvel[i-1]*xvel[i-1] + yvel[i-1]*yvel[i-1]);
    xacc = -GM * xpos[i-1] / (dist*dist*dist);
    yacc = -GM * ypos[i-1] / (dist*dist*dist);
    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;
    time[i] = i*dT;
    ePot[i] = -GM * m / dist;
    eCin[i] = 0.5 * m * vel*vel;
    eTot[i] = ePot[i] + eCin[i];
    //L = m * (r x v) !!! r et v ne sont pas toujours orthogonaux
    l[i] = m * (xpos[i] * yvel[i] - ypos[i] * xvel[i]);
  }

  //partie graphique avec DISLIN
  metafl("XWIN");   //XWIN ou PDF
  page(3000, 3000);   //pour creer une image proportionee
  disini();     //initialisation de DISLIN

  //dessin de l'orbite
  name("axe-X [UA]", "X");
  name("axe-Y [UA]", "Y");
  titlin("Orbite elliptique",1); 
  graf(-7.,3.,-7.,1.,-5.,5.,-5.,1.);
  title();
  thkcrv(5);
  color("RED");
  curve(xpos, ypos, steps);
  //dessin du soleil
  incmrk(-1);
  color("GREEN");
  marker(21);
  double xx, yy;
  xx=0.;
  yy=0.;
  curve(&xx,&yy,1);   //la fonction s'attend un adresse memoire
  endgrf();
  erase();
  system("PAUSE");    //avant de dessiner le deuxieme graph
    
  //dessin de l'energie
  name("temps", "X");
  name("energie", "Y");
  titlin("Orbite elliptique",1);
  titlin("jaune = moment cinetique",2);
  titlin("vert = energie total",3);
  titlin("rouge / blue = energie potentiele / cinetique",4);
  color("FORE");
  graf(0., dT*steps, 0., dT*steps/5., -2., 2., -2., 0.5);
  title();
  incmrk(0);
  thkcrv(2);
  color("RED");
  curve(time, ePot, steps);
  color("BLUE");
  curve(time, eCin, steps);
  color("GREEN");
  curve(time, eTot, steps);
  color("YELLOW");
  curve(time, l, steps);
  
  disfin();    //fin de DISLIN

  return 0;
}