//etude de la trajectoire d'une balle dans l'air avec frottement
//avec la methode d'Euler
#include <iostream>
#include <cmath>
#include <fstream>
#include "dislin.h"

using namespace std;

const double DEG2RAD = M_PI/180.;

int main() {
  const int STEPS = 10000;
  double xpos[STEPS]={0.}, ypos[STEPS]={0.};   //position
  double xvel[STEPS]={0.}, yvel[STEPS]={0.};   //vitesse
  double xpos0[STEPS]={0.}, ypos0[STEPS]={0.};   //position sans resistance
  double xvel0[STEPS]={0.}, yvel0[STEPS]={0.};   //vitesse sans frottement
  double xposTheo[STEPS]={0.}, yposTheo[STEPS]={0.};   //position theorique sans resistance
       
  //declaration des conditions initiales: position, vitesse, angle
  xpos[0] = xpos0[0] = xposTheo[0] = 0.;   //position initiale pour tous les calculs
  ypos[0] = ypos0[0] = yposTheo[0] = 0.;
  double velini, theta;
  cout << "Entrez la vitesse initial (en m/s) : ";
  cin >> velini;   cout << endl;
  cout << "Entrez l'angle initial (en degres) : ";
  cin >> theta;   cout << endl; 
  xvel[0] = xvel0[0] = velini*cos(theta*DEG2RAD);    //vitesse initiale pour tous les calculs
  yvel[0] = yvel0[0] = velini*sin(theta*DEG2RAD);

  //on considere une balle avec un diametre de 7.4 cm et une masse m = 0.145kg (baseball)
  double drag = 0.1;
  cout << "Entrez le coefficient de frottement C, F = C v**2 [C < 1 kg/m] : "; 
  cin >> drag;   cout << endl;
  double area = 4.3e-3;        //section de la balle
  const double G = 9.81;       //acceleration gravitationnelle
  double m = 0.145;            //mass de la balle
  double rho = 1.2;            //densite de l'air (kg/m3)
  //resistance dans l'air
  double airRes = -0.5*drag*rho*area/m;   

  //calcul de la trajectoire
  double tau;
  cout << "Entrez un pas de temps pour l'integration delta t en s : ";
  cin >> tau;   cout << endl;
  int count = 0, count0 = 0;
  double vel, xacc, yacc, t;   //vbariables auxiliares
  
  //on boucle sur le nombre de pas maximal tant que la trajectoire
  //ne passe pas en dessous de 0 le long de l'axe y
  for (int istep=0; istep<STEPS-1; istep++) {
    //trajectorie avec resistance
    //mise a jour de la position
    xpos[istep+1] = xpos[istep] + tau*xvel[istep];
    ypos[istep+1] = ypos[istep] + tau*yvel[istep];
    //acceleration
    vel = sqrt(pow(xvel[istep],2)+pow(yvel[istep],2));   //modul de la vitesse
    xacc = airRes * vel * xvel[istep];   //nul si pas de resistance
    yacc = airRes * vel * yvel[istep];   //nul si pas de resistance
    yacc -= G;                          //ajoute a yacc de la composante gravitationnelle
    //mise a jour de la vitesse
    xvel[istep+1] = xvel[istep] + tau*xacc;
    yvel[istep+1] = yvel[istep] + tau*yacc;

    //trajectoire sans resistance
    //mise a jour de la position
    xpos0[istep+1] = xpos0[istep] + tau*xvel0[istep];
    ypos0[istep+1] = ypos0[istep] + tau*yvel0[istep];
    //acceleration
    xacc = 0.;
    yacc = -G;                          //composante gravitationnelle
    //mise a jour de la vitesse
    xvel0[istep+1] = xvel0[istep] + tau*xacc;
    yvel0[istep+1] = yvel0[istep] + tau*yacc;
    
    //trajectoire theorique sans resistance
    t = istep*tau;       //temps au pas actuel
    xposTheo[istep] = xpos[0] + xvel[0]*t;
    yposTheo[istep] = ypos[0] + yvel[0]*t - 0.5*G*t*t; 

    if (ypos[istep+1]>0.)
      count++;
    else {   //la balle a attainte le sol (avec resistance); elle ne bouge plus
      xpos[istep+1] = xpos[istep];
      ypos[istep+1] = 0.;
      xvel[istep+1] = 0.;
      yvel[istep+1] = 0.;
    }
    count0++;
        
    //on sort de la boucle si la balle atteint le sol
    if (ypos0[istep+1]<0.)
      break;
  }

  //affichage de la jetee maximale et le temps de vol
  cout << "Jetee maximale avec resistance : "  << xpos0[count+1]  << " m," << endl;
  cout << "et temps de vol : " << (count+1)*tau << " s." << endl << endl;
  cout << "Jetee maximale sans resistance : "  << xpos0[count0]  << " m," << endl;
  cout << "et temps de vol : " << count0*tau << " s." << endl << endl;

  //comparaison avec les valeurs theoriques 
  double range = 2*xvel[0]*yvel[0]/G;
  double tof = 2*yvel[0]/G;
  cout << "Jetee maximale theorique sans frottement : "  << range  << " m," << endl;
  cout << "et temps de vol : " << tof << " s." << endl << endl;

  //copie les positions dans un fichier text
  ofstream trajectory("trajectoire.dat");
  for (int i=0; i<count0; i++)
    trajectory  << tau*i << "\t" << xpos[i] << "\t" << ypos[i] << "\t" 
                << xpos0[i] << "\t" << ypos0[i] << "\t" 
                << xposTheo[i] << "\t" << yposTheo[i] << endl;

  //dessin avec dislin
  metafl("XWIN");   //XWIN ou PDF
  disini();     //initialisation de dislin
  name("jetee (m)", "X");
  name("hauteur (m)", "Y");
  titlin("Trajectoire", 1);
  graf(0.,2000.,0.,500.,0.,500.,0.,100.);
  title();

  thkcrv(10);
  color("RED");
  curve(xpos, ypos, count+1);
  color("GREEN");
  curve(xpos0, ypos0, count0);
  color("YELLOW");
  incmrk(-10);   //symboles au lieu de la curbe
  marker(21);   //symbole pour le plot
  curve(xposTheo, yposTheo, count0-1);
  color("FORE");
  
  disfin();   //fin de DISLIN
 
  return 0;
}