//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  
  const double beta = 2.1;   //exposant loi gravitationnelle
  
  double dT = 0.002;   //pas d'integration

  //animation avec DISLIN
  metafl("XWIN");
  page(3000, 3000);
  disini();     //initialisation de dislin
  name("axe-X","X");
  name("axe-Y","Y");
  titlin("Deviation de la loi de Newton 1/r2",1);
  
  //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;

    //animation avec DISLIN
    //tableau de position du satellite pour animation
    double xsat[1]={xpos[i]};
    double ysat[1]={ypos[i]};

    //poisition de la Terre pour animation
    double X_M[1]={0.}, Y_M[1]={0.};
		
    //pour que l'animation ne soit pas trop lente, on ne dessine que toutes les 200 iterations
    if (i%200==0) {
      erase();        //on efface le tableau
      color("white");    //on remet un axe decrit en dessous en blanc
      graf(-2.,2.,-2.,1.,-2.,2.,-2.,1.);
      title();
      
      //commandes pour dessiner une terre en vert
      incmrk(-1);
      color("green");
      marker(21);
      curve(X_M,Y_M,1);

      //commande pour dessiner la trajectoire de la planete
      incmrk(0);
      thkcrv(2);
      color("red");
      curve(xpos, ypos, i);

      //commande pour afficher la planete
      incmrk(-1);
      color("yellow");
      marker(21);
      curve(xsat, ysat, 1);

      sendbf();
      endgrf();
    }
  }

  disfin();     //fin de DISLIN

  return 0;
}