//etude d'une orbite avec frottement dans la haut atmosphere
#include <iostream>
#include <cmath>
#include "dislin.h"

using namespace std;

int main() {
  const double RT = 6378.;   //en km
  const double GM = 398645.;   //G_N * masse_terre
  double altIni = 100.;
  double altFin = 50.;
  const double R0 = RT + altIni;
  const double V0 = sqrt(GM/R0);
  const double CX = 0.0001;     //coefficient de frottement
  
  const int steps = 50000; 
  double xpos[steps], ypos[steps], xvel[steps], yvel[steps]; 
  //conditions initiales
  xpos[0] = R0; ypos[0] = 0.;
  xvel[0] = 0.; yvel[0] = V0;
  
  double dT = 2.;     //pas d'integration
  
  //methode de Runge
  double dist, xacc, yacc;     //variables auxiliares
  //redifinition de la vitesse initiale
  dist = sqrt(xpos[0]*xpos[0] + ypos[0]*ypos[0]);
  xacc = -GM / pow(dist,3) * (xpos[0] - CX*ypos[0]);   //acceleration
  yacc = -GM / pow(dist,3) * (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;
  //boucle
  for (int i=1; i<steps; i++) {
    dist = sqrt(xpos[i-1]*xpos[i-1] + ypos[i-1]*ypos[i-1]);
    xacc = -GM / pow(dist,3) * (xpos[i-1] - CX*ypos[i-1]);
    yacc = -GM / pow(dist,3) * (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 (dist < RT) {   //le satellite est tombe sur la Terre
      cout << "Le satellite est tombe sur la Terre !" << endl;
      break;
    }
  }
   
  cout << "apres " << int(nsteps*dT) << " seconds" << endl;
  cout << "et " << norbits/2 << " orbites";
  if (norbits%2 == 1) cout << " et demi";
  cout << endl;
  
  //partie graphique vec 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;
}