//Programme pour calculer la trajectoire d'une balle
//de golf dans l'air avec la methode d'Euler
//----------------------------------------------------
#include <iostream>
#include <cstdlib>
#include <cmath>
#include <fstream>
#include "dislin.h"

using namespace std;

const double G = 9.81;    //acceleration gravitationnelle m/s2
const double RHO = 1.2;   //densite de l'air en kg/m3

int main() {
  //sesie des conditions initiales
  double height, velocity, angle;
  cout << "Entrez la hauteur initiale (m) : ";
  cin >> height;
  cout << "Entrez la vitesse initiale (m/s) : ";
  cin >> velocity;
  cout << "Entrez l'angle initial (degrees) : ";
  cin >> angle;

  //vecteurs de position, vitesse, acceleration
  //on etudiera 4 trajectoires differentes
  double rIni1[2], rNew1[2];
  double rIni2[2], rNew2[2];
  double rIni3[2], rNew3[2];
  double rIni4[2], rNew4[2];
  double vIni1[2], vNew1[2];
  double vIni2[2], vNew2[2];
  double vIni3[2], vNew3[2];
  double vIni4[2], vNew4[2];
  double a1[2], a2[2], a3[2], a4[2];
  
  //initialisation : conditions initiales
  rIni1[0] = rIni2[0] = rIni3[0] = rIni4[0] = 0;
  rIni1[1] = rIni2[1] = rIni3[1] = rIni4[1] = height;         
  vIni1[0] = vIni2[0] = vIni3[0] = vIni4[0] = velocity*cos(angle*M_PI/180.);    
  vIni1[1] = vIni2[1] = vIni3[1] = vIni4[1] = velocity*sin(angle*M_PI/180.);
  rNew1[0] = rNew2[0] = rNew3[0] = rNew4[0] = rIni1[0];                      
  rNew1[1] = rNew2[1] = rNew3[1] = rNew4[1] = rIni1[1];
  vNew1[0] = vNew2[0] = vNew3[0] = vNew4[0] = vIni1[0];
  vNew1[1] = vNew2[1] = vNew3[1] = vNew4[1] = vIni1[1];

  //declaration des parametres
  double cDragNormal = 7.0;   //coef. frottement pour une balle normale (~7.0/v)
  double cDragSmooth = 0.5;   //coef. frottement pour une balle lisse
  double cSpin = 0.25;        //coefficient de spin en 1/s
  double xSection = 1.45e-3;  //section pur un balle de diametre 0.043 m
  double m = 0.046;           //masse en kg
  
  int count;
  double tau;
  cout << "Entrez un pas d'integration, tau (sec) : ";
  cin >> tau;
  int maxStep;
  cout << "Entrez le nombre d'iterations : ";
  cin >> maxStep;

  //allocation dynamique des vecteurs pour stocker les trajectoires
  double *xNormal = new double[maxStep];     //trajectoire avec spin, balle normale        
  double *yNormal = new double[maxStep];
  double *xNoSpin = new double[maxStep];     //trajectoire sans spin, balle normale
  double *yNoSpin = new double[maxStep];
  double *xSmooth = new double[maxStep];     //trajectoire avec spin, balle lisse
  double *ySmooth = new double[maxStep];
  double *xExtraS = new double[maxStep];     //trajectoire avec extra spin (+50%), balle normale        
  double *yExtraS = new double[maxStep];

  for (int iStep=0; iStep < maxStep; iStep++) {
    //sauve la position a chaque debut d'iteration
    xNormal[iStep] = rNew1[0];
    yNormal[iStep] = rNew1[1];
    xNoSpin[iStep] = rNew2[0];
    yNoSpin[iStep] = rNew2[1];
    xSmooth[iStep] = rNew3[0];
    ySmooth[iStep] = rNew3[1];
    xExtraS[iStep] = rNew4[0];
    yExtraS[iStep] = rNew4[1];

    //temps au pas actuel
    double t = iStep*tau; 

    //calcul de l'acceleration (4 fois)
    double normV1 = sqrt(vNew1[0]*vNew1[0] + vNew1[1]*vNew1[1]);
    a1[0] = -(1/m)*cDragNormal/normV1*RHO*xSection*pow(normV1,2)*(vNew1[0]/normV1);   //drag en x
    a1[0] -= cSpin*vNew1[1];                                                            //spin en x
    a1[1] = -(1/m)*cDragNormal/normV1*RHO*xSection*pow(normV1,2)*(vNew1[1]/normV1);   //drag en y
    a1[1] += cSpin*vNew1[0];                                                            //spin en y
    a1[1] -= G;                                                                    //acc. gravitationnelle

    double normV2 = sqrt(vNew2[0]*vNew2[0] + vNew2[1]*vNew2[1]);
    a2[0] = -(1/m)*cDragNormal/normV2*RHO*xSection*pow(normV2,2)*(vNew2[0]/normV2);   //drag en x
    a2[1] = -(1/m)*cDragNormal/normV2*RHO*xSection*pow(normV2,2)*(vNew2[1]/normV2);   //drag en y
    a2[1] -= G;                                                                    //acc. gravitationnelle

    double normV3 = sqrt(vNew3[0]*vNew3[0] + vNew3[1]*vNew3[1]);
    a3[0] = -(1/m)*cDragSmooth*RHO*xSection*pow(normV3,2)*(vNew3[0]/normV3);   //drag en x
    a3[0] -= cSpin*vNew3[1];                                                   //spin en x
    a3[1] = -(1/m)*cDragSmooth*RHO*xSection*pow(normV3,2)*(vNew3[1]/normV3); //drag en y
    a3[1] += cSpin*vNew3[0];                                                   //spin en y
    a3[1] -= G;                                                                //acc. gravitationnelle

    double normV4 = sqrt(vNew4[0]*vNew4[0] + vNew4[1]*vNew4[1]);
    a4[0] = -(1/m)*cDragNormal/normV4*RHO*xSection*pow(normV4,2)*(vNew4[0]/normV4);   //drag en x
    a4[0] -= cSpin*vNew4[1]*1.5;                                                      //extra spin en x
    a4[1] = -(1/m)*cDragNormal/normV4*RHO*xSection*pow(normV4,2)*(vNew4[1]/normV4);   //drag en y
    a4[1] += cSpin*vNew4[0]*1.5;                                                      //extra spin en y
    a4[1] -= G;                                                                    //acc. gravitationnelle
        
    //mise a jour des positions et vitesses
    rNew1[0] += tau*vNew1[0];
    rNew1[1] += tau*vNew1[1];
    vNew1[0] += tau*a1[0];
    vNew1[1] += tau*a1[1];

    rNew2[0] += tau*vNew2[0];
    rNew2[1] += tau*vNew2[1];
    vNew2[0] += tau*a2[0];
    vNew2[1] += tau*a2[1];

    rNew3[0] += tau*vNew3[0];
    rNew3[1] += tau*vNew3[1];
    vNew3[0] += tau*a3[0];
    vNew3[1] += tau*a3[1];

    rNew4[0] += tau*vNew4[0];
    rNew4[1] += tau*vNew4[1];
    vNew4[0] += tau*a4[0];
    vNew4[1] += tau*a4[1];

    count = iStep;

    if (rNew1[1]<0.) rNew1[1] = 0.;
    if (rNew2[1]<0.) rNew2[1] = 0.;
    if (rNew3[1]<0.) rNew3[1] = 0.;
    if (rNew4[1]<-2.) {
      //xNormal[iStep+1] = rNew1[0];      
      //yNormal[iStep+1] = rNew1[1];
      xExtraS[iStep+1] = rNew4[0];      
      yExtraS[iStep+1] = rNew4[1];
      break;                            
    }
  }

  //affiche la distance maximale et le temps de vol
  cout << "Distance parcourue (avec frottements, spin, balle normale: "  << rNew1[0]  << " metres." << endl;
  cout << "Temps de vol: " << count*tau << " secondes." << endl;

  //comparaison avec les valeurs theoriques 
  double range = 2*(vIni1[0]*vIni1[0] + vIni1[1]*vIni1[1])*cos(angle*M_PI/180.)*sin(angle*M_PI/180.)/G;
  double tof = 2*sqrt(vIni1[0]*vIni1[0] + vIni1[1]*vIni1[1])*sin(angle*M_PI/180.)/G;
  cout << "Distance maximale theorique (sans frottements, sans spin) : "  << range  << " metres." << endl;
  cout << "Temps de vol theorique : " << tof << " secondes." << endl;

  //copie les positions dans un fichier text
  ofstream trajectory("golf.dat");
  for (int i=0; i<(count+1); i++) {
    trajectory << xNormal[i] << " " << yNormal[i] << " "
               << xNoSpin[i] << " " << yNoSpin[i] << " "
               << xSmooth[i] << " " << ySmooth[i] << " "
               << xExtraS[i] << " " << yExtraS[i] << endl;
  }
    
  //DISLIN
  metafl("XWIN");   //XWIN ou PDF
  disini();
  name("distance (m)","X");
  name("hauteur (m)","Y");
  titlin("Trajectoire d'une balle de golf",1);
  color("white");
  thkcrv(10);
  lintyp(0);
  graf(0., xNormal[count]*1.1, 0., 50., 0., yExtraS[int(count/2)]*1.3, 0., 20.);
  title();
  color("magenta");
  curve(xNormal,yNormal,count);
  color("cyan");
  curve(xExtraS,yExtraS,count);
  color("green");
  curve(xNoSpin,yNoSpin,count);
  color("yellow");
  curve(xSmooth,ySmooth,count);
  disfin();
    
  //libere la memoire allouee pendant l'execution du programme
  delete [] xNormal, yNormal;
  delete [] xNoSpin, yNoSpin;
  delete [] xSmooth, ySmooth;
  delete [] xExtraS, yExtraS;
  
  system("PAUSE");  
  return 0;
}