//chute d'un objet dans l'air
#include <iostream>
#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 double M = 90.;
  const double k1 = 2.;
  const double k2 = 5.;
  
  const int STEPS = 75000;
  double ypos[STEPS];   //position
  double yvel[STEPS];   //vitesse
  double temps[STEPS];
  
  //conditions intiales
  double h = 3500.;
  cout << "Entrez l'hauteur initiale (m) : ";
  cin >> h;
  ypos[0] = h;        //position initiale                
  yvel[0] = 0.;       //vitiesse initiale

  //calcul de la trajectoire (chute)
  double dT;
  cout << "Entrez un pas de temps pour l'integration, dT (sec) : ";
  cin >> dT;
  
  double drag1, drag2, yacc;
  //methode de Runge
  drag1 = k1/M * 1./50. * Rho0 * exp(-ypos[0]/H0);
  drag2 = k2/M * 1./50. * Rho0 * exp(-ypos[0]/H0);
  yacc = drag1 * yvel[0] + drag2 * yvel[0]*yvel[0]- 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
    if (temps[istep-1] < 60.) {     //parachute ferme
      drag1 = k1/M * 1./50. * Rho0 * exp(-ypos[istep-1]/H0);
      drag2 = k2/M * 1./50. * Rho0 * exp(-ypos[istep-1]/H0);
      yacc = drag1 * yvel[istep-1] + drag2 * yvel[istep-1]*yvel[istep-1];
    } else {     //parachute ouvert
      drag1 = k1/M * Rho0 * exp(-ypos[istep-1]/H0);
      drag2 = k2/M * Rho0 * exp(-ypos[istep-1]/H0);
      yacc = drag1 * yvel[istep-1] + drag2 * yvel[istep-1]*yvel[istep-1];
    } 
    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.,-200.,0.,-200.,50.);
  thkcrv(10);
  color("BLUE");
  curve(temps, yvel, count);
  //fin de DISLIN
  disfin();

  return 0;
}