//chute d'un objet dans l'air
#include <iostream>
#include <cstdlib>
#include <cmath>
#include <fstream>
//#include "dislin.h"

using namespace std;

int main() {
  //parametres
  const double G = 9.81;
  const double Rho0 = 1.2;
  const double H0 = 10000.;

  const int STEPS = 50000;
  double ypos[STEPS];   //position
  double yvel[STEPS];   //vitesse
  double temps[STEPS];
  
  //conditions intiales
  double h = 10000.;
  cout << "Entrez l'hauteur initiale (m) : ";
  cin >> h;
  ypos[0] = h;        //position initiale                
  yvel[0] = 0.;       //vitiesse initiale

  //resistance de l'air
  double drag0 = 0.1; 
  cout << "Entrez le coefficient de resistance dans l'air : ";      
  cin >> drag0;

  //calcul de la trajectoire (chute)
  double dT;
  cout << "Entrez un pas de temps pour l'integration, dT (sec) : ";
  cin >> dT;
  
  double drag, yacc;
  //methode de Runge
  drag = drag0 * Rho0 * exp(-ypos[0]/H0);
  yacc = drag * pow(yvel[0],2) - G;
  yvel[0] = yvel[0] - yacc*dT/2.;
  //on boucle tant que l'objet ne touche pas le sol
  int count = 1;
  bool term = false;
  for (int istep=1; istep<STEPS; istep++) {
    //acceleration
    drag = drag0 * Rho0 * exp(-ypos[istep-1]/H0);
    yacc = drag * pow(yvel[istep-1],2);   //frottement, nul si pas de resistance
    yacc = yacc - G;          //addition de la composante gravitationnelle
    //mise a jour de la vitesse
    yvel[istep] = yvel[istep-1] + yacc*dT;
    //mise a jour de la position
    ypos[istep] = ypos[istep-1] + yvel[istep]*dT;
    temps[istep] = istep*dT;

    //velocite termianle : < 0.01% da variation
    if ((!term) && abs((yvel[istep]-yvel[istep-1])/yvel[istep])<0.0001) {
      cout << "Veleocite terminale : " << yvel[istep]
        << " m/s attainte apres " << count*dT << " s a "
        << ypos[istep] << " m." << endl;
      term = true;  
    }    
    
    //on sort de la boucle si l'objet atteint le sol
    if (ypos[istep]<0.) break;

    count++;
  }
  cout << "Velocite finale au sol : " << yvel[count] << " m/s apres " 
    << count*dT << " s de chute." << endl;
  if (!term)
    cout << "On n'a pas attainte la velocite terminale :-(" << endl;

   //initialisation de DISLIN
   metafl("XWIN");     //XWIN ou PDF
   disini();
   //dessin de la trajectoire
   name("temps (s)", "X");
   name("hauteur (m)", "Y");     
   graf(0.,500.,0.,100.,0.,10000.,0.,1000.);
   thkcrv(10);
   color("RED");
   curve(temps, ypos, count);
   color("FORE");
   endgrf();
   //dessin de la vitesse
   newpag();
   name("temps (s)", "X");
   name("vitesse (m/s)", "Y");     
   graf(0.,500.,0.,100.,-100.,0.,-100.,20.);
   thkcrv(10);
   color("BLUE");
   curve(temps, yvel, count);
   //fin de DISLIN
   disfin();

  system("PAUSE");
  return 0;
}