//etude d'une orbite avec la methode de Runge
//avec frottement
//avec animation :-)
#include <iostream>
#include <cstdlib>
#include <cmath>
#include "dislin.h"

using namespace std;

int main() {
  //declaration de toutes les constantes que nous allons utiliser dans les caculs
  const double RT = 6378.;
  const double GM = 398645.;   //G * masse_terre
  double alt = 1000.;
  const double R0 = RT + alt;
  const double V0 = sqrt(GM/R0);
  const double CX = 0.0015;   //coefficient de frottement

  //declaration du nombre de points et du delat t
  const double dT = 1.;     //pas d'integration
  const int steps = 50000;
  
  //declaration des tableaux pour enregistrer les points calcules pour le satellite
  double xpos[steps], ypos[steps], xvel[steps], yvel[steps]; 
  //declaration des conditions initiales pour le satellite
  xpos[0] = R0; ypos[0] = 0.;
  xvel[0] = 0.; yvel[0] = V0;

  //animation avec DISLIN
  metafl("XWIN");
  page(3000, 3000);
  disini();     //initialisation de dislin
  name("axe-X","X");
  name("axe-Y","Y");
  
  //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 / (normR*normR*normR) * (xpos[0]-CX*ypos[0]) ;   //force ou acceleration
  yacc = - GM / (normR*normR*normR) * (ypos[0]+CX*xpos[0]);    
  xvel[0] = xvel[0] - xacc*dT/2.;   //vitesse 
  yvel[0] = yvel[0] - yacc*dT/2.;

  //boucle de calcul des vitesses et positions avec la methode de Runge
  int nsteps = 0;
  int norbits = 0;
  for (int i=1; i<steps; i++) {
    normR = sqrt(xpos[i-1]*xpos[i-1] + ypos[i-1]*ypos[i-1]);
    xacc = - GM / (normR*normR*normR) * (xpos[i-1]-CX*ypos[i-1]) ;   //acceleration
    yacc = - GM / (normR*normR*normR) * (ypos[i-1]+CX*xpos[i-1]);
    
    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;

    //tableaux des positions 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 100 iterations
    if (i%100==0) {
      erase();    //on efface le tableau
      color("white");    //on remet un axe decrit en dessous en blanc
      graf(-7500., 7500., -7500., 2500., -7500., 7500., -7500., 2500.);
      
      //commandes pour dessiner la terre en vert
      incmrk(-1);
      color("green");
      marker(21);
      curve(X_M,Y_M,1);
      
      //commande pour dessiner la trajectoire du satellite avec frottements en rouge
      incmrk(0);
      thkcrv(2);
      color("red");
      curve(xpos, ypos, i);

      //commande pour afficher le satellite avec frottements en jaune
      incmrk(-1);
      color("yellow");
      marker(21);
      curve(xsat, ysat, 1);

      sendbf();
      endgrf();
    }
    
    nsteps++;
    if (ypos[i-1]*ypos[i] < 0.) norbits = norbits + 1;
    if (normR < RT) {   //le satellite est tombe sur la Terre
      cout << "Le satellite est tombe sur la Terre !" << endl;
      break;
    }
  }
  cout << "temps : " << int(nsteps*dT) << " seconds" << endl;
  cout << "orbites : " << norbits/2;
  if (norbits%2 == 1) cout << " et demi";
  cout << endl;

  //fin de DISLIN
  disfin();

  return 0;
}