//chute d'un objet dans l'air
#include <iostream>
#include <cmath>
#include <fstream>

using namespace std;

const double G = 9.81; //acceleration gravitationnelle m/s^2

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

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

  //calcul de la trajectoire (chute)
  double tau;
  cout << "Entrez un pas de temps pcour l'integration, tau (sec) : ";
  cin >> tau;

  double yacc;
  //methode de Runge
  yacc = drag * pow(yvel[0],2) - G;
  yvel[0] = yvel[0] - yacc*tau/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    
    yacc = drag * pow(yvel[istep-1],2);   //nul si pas de frottement
    yacc = yacc - G;          //addition de la composante gravitationnelle
    //mise a jour de la vitesse
    yvel[istep] = yvel[istep-1] + yacc*tau;
    //mise a jour de la position
    ypos[istep] = ypos[istep-1] + yvel[istep]*tau;
    temps[istep] = istep*tau;

    //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*tau << " 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*tau << " s de chute." << endl;
  if (!term)
    cout << "On n'a pas attainte la velocite terminale :-(" << endl;
  
  //sortie des positions dans un fichier text
  ofstream trajectory("chute.dat");
  for (int i=0; i<count; i++)
    trajectory << ypos[i] << "\t" << yvel[i] << "\t" << temps[i] << endl;

  return 0;
}