//etude d'une orbite avec la methode de Runge
//avec frottement
#include <iostream>
#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.001;     //coefficient de frottement
  
  const int steps = 60000; 
  double xpos[steps], ypos[steps], xvel[steps], yvel[steps]; 
  xpos[0] = R0; ypos[0] = 0.0;   //conditions initiales
  xvel[0] = 0.0; yvel[0] = V0;
  
  double dT = 1.;     //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 / (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.;
  
  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]) ;   //force ou 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;
    
    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;
  
  //partie graphique - DISLIN
  metafl("XWIN");   //XWIN ou PDF
  page(3000, 3000);
  disini();     //initialisation de dislin
  name("axe-X","X");
  name("axe-Y","Y");
  graf(-7500.,7500.,-7500.,2500.,-7500.,7500.,-7500.,2500.);
  thkcrv(2);
  color("yellow"); 
  curve(xpos, ypos, nsteps);  

  disfin();     //fin de DISLIN

  return 0;
}