//etude de rebonds d'une balle sur deux plans inclines
#include <iostream>
#include <cmath>
#include <iomanip>

using namespace std;

//declarations des fonctions
double angleRebond(double theta, double alpha, double x);
double solutionGauche(double v, double theta, double alpha, double x0, double y0);
double solutionDroite(double v, double theta, double alpha, double x0, double y0);
double intersect(double v, double theta, double alpha, double x0, double y0);
double angleIncidence(double v, double theta, double alpha, double x0, double y0, double x);

//on utilise des unites MKS
const double GN = 9.81;

int main() {
  //saisi des parametre initials
  double h;
  cout << "Entrez la hauteur initiale h en metres : ";
  cin >> h;
  double x0;
  cout << "\nEntrez la position horizontale x en metres : ";
  cin >> x0;
  double alpha;
  cout << "\nEntrez l'angle d'inclination des plans en degres : ";
  cin >> alpha;
  alpha *= M_PI/180.;   //conversion de l'angle en radians
  int nrebond;
  cout << "\nEntrez le nombre de rebonds : ";
  cin >> nrebond;
  cout << endl;

  while(h <= fabs(tan(alpha)*x0)) {
    cout << "La hauteur initiale est au dessous du plan." << endl;
    cout << "Essayez des nouvelles donnees (meme inclination)" << endl;
    cout << "Entrez la nouvelle hauteur initiale h : ";
    cin >> h;
    cout << "\net la position horizontale x : " << endl;
    cin >> x0;
  }
        
  //premier rebond
  double y0, v0, t0;
  y0 = fabs(tan(alpha)*x0);
  v0 = sqrt(2.*GN*(h-y0));
  t0 = sqrt(2.*(h-y0)/GN);

  //affichage
  cout << endl << endl;
  cout << "rebond   position x0   position y0    vitesse v0         temps";
  cout << endl << endl;
  cout << setw(6) << 1
       << setw(14) << setprecision(7) << x0
       << setw(14) << setprecision(7) << y0
       << setw(14) << setprecision(7) << v0
       << setw(14) << setprecision(7) << t0 << endl;

  double thetaIn, thetaOut;
  thetaIn = -M_PI/2.;   //angle d'incidence initial en radians
  thetaOut = angleRebond(thetaIn,alpha,x0);

/*
Ensuite, il faut determiner l'intersection d'une parabole :
y = y0 + tan(thetaOut)*(x-x0) - g*(x-x0)^2/(2*v0^2*cos(thetaOut)^2)
avec la fonction qui decrit les plans y = tan(alpha)*|x|.

On doit donc resoudre deux equations quadratiques :
(1) x < 0 :
g*(x-x0)^2/(2*v0^2*cos(thetaOut)^2)-tan(thetaOut)*(x-x0)-y0-tan(alpha)*x=0
  pour y = -tan(alpha)*x   --> x < 0
ou
(2) x > 0 :
g*(x-x0)^2/(2*v0^2*cos(thetaOut)^2)-tan(thetaOut)*(x-x0)-y0+tan(alpha)*x=0
  pour y = tan(alpha)*x   --> x > 0

Les solutions de ces equations sont de la forme :
  x = (-B_i +- sqrt( B_i^2- 4*A*C ))/(2*A)     i=1,2
avec 
  A =  g/(2*v0^2*cos(thetaOut)^2), 
  B1= -x0*g/(v0^2*cos(thetaOut)^2) - tan(thetaOut) -tan(alpha)
  B2= -x0*g/(v0^2*cos(thetaOut)^2) - tan(thetaOut) +tan(alpha)
  C =  g*x0^2/(2*v0^2*cos(thetaOut)^2) + tan(thetaOut)*x0 -y0
*/
  
  double x, t;
  for(int compt=1; compt<nrebond; compt++) {
    if (fabs(thetaOut-M_PI/2.)>0.00001) {
      x = intersect(v0,thetaOut,alpha,x0,y0);
      thetaIn = angleIncidence(v0,thetaOut,alpha,x0,y0,x);
    }
    else {
      x = x0;
      thetaIn = thetaOut;
    }
    
    if (fabs(x0-x)>0.00001)
      t = fabs((x0-x) / (v0*cos(thetaOut)));
    else
      t = 2.*v0*fabs(sin(thetaOut)) / GN;

    thetaOut = angleRebond(thetaIn,alpha,x);
    x0 = x;
    y0 = tan(alpha)*fabs(x);
    v0 = sqrt(2.*GN*(h-y0));

    //affichage
    cout << setw(6) << compt+1  
         << setw(14) << setprecision(7) << x0
         << setw(14) << setprecision(7) << y0
         << setw(14) << setprecision(7) << v0
         << setw(14) << setprecision(7) << t << endl;
  }

  return 0;
}

//definition de la fonction qui calcule l'angle de rebond
double angleRebond(double theta, double alpha, double x) {
  if (x<0.)     //rebond sur le plan a gauche
    return -2.*alpha-theta;
  else if (x>0.)     //rebond sur le plan a droite
    return 2.*alpha-theta;
  else     //rebond entre les deux plans (x=0): approxime a un plan horizontal
    return M_PI-theta;
}

//definition de la fonction qui calcule l'intersection entre
//la trajectoire de la balle et le plan incline
double intersect(double v, double theta, double alpha, double x0, double y0) {
  double x1 = solutionGauche(v,theta,alpha,x0,y0);
  double x2 = solutionDroite(v,theta,alpha,x0,y0);
  if (x1>0 && x2>=0)
    return x2;
  else if (x2<0 && x1<0)
    return x1;
  else {
    cout << "Il y a quelque probleme, deux solutions !? " << endl;
    cout << "x1 = " << setw(10) << setprecision(7) << x1
         << ", x2 = " << setw(10) << setprecision(7) << x2 << endl;
    return -99999.99;
  }
}

//definition de la fonction qui calcule l'intersection a gauche (x < 0)
double solutionGauche(double v, double theta, double alpha, double x0, double y0) {
  double den = v*v*cos(theta)*cos(theta);
  double a = GN/(2.*den);
  double b1= -tan(theta) - GN*x0/den - tan(alpha);
  double c = GN*x0*x0/(2.*den) + tan(theta)*x0 -y0;
  double x1, x2;
  double disc = b1*b1-4.*a*c;

  if (disc>0) {
    x1 = (-b1 + sqrt(disc))/(2.*a);
    x2 = (-b1 - sqrt(disc))/(2.*a);
  } 
  else if(disc==0)
    x1 = x2 = -b1/(2.*a);
  else {
    //cout << "Pas de solution possible pour x < 0" << endl;
    return 99999.9;
  }
 
  //choisir une solution x < 0, qui ne soit pas la position initiale
  if (x1<0 && 0.01<fabs(x1-x0))
    return x1;
  else if (x2<0 && 0.01<fabs(x2-x0))
    return x2;
  else { 
    //cout << "Pas de solution possible pour x < 0" << endl;
    return 99999.9;
  }
}

//definition de la fonction qui calcule l'intersection a droite (x > 0)
double solutionDroite(double v, double theta, double alpha, double x0, double y0) {
  double den = v*v*cos(theta)*cos(theta);
  double a = GN/(2.*den);
  double b2= -tan(theta) - GN*x0/den + tan(alpha);
  double c = GN*x0*x0/(2.*den) + tan(theta)*x0 -y0;
  double x1, x2;
  double disc = b2*b2-4.*a*c;

  if (disc>0) {
    x1 = (-b2 + sqrt(disc))/(2.*a);
    x2 = (-b2 - sqrt(disc))/(2.*a);
  }
  else if (disc==0)
    x1 = x2 = -b2/(2.*a);
  else {
    //cout << "Pas de solution possible pour x > 0" << endl;
    return -99999.9;
  }
  
  //choisir une solution x > 0, qui ne soit pas la position initiale
  if (x1>=0 && fabs(x1-x0)>0.001)
    return x1;
  else if (x2>=0 && fabs(x2-x0)>0.001)
    return x2;
  else { 
    //cout << "Pas de solution possible pour x>0" << endl;
    return -99999.9;
  }
}

//definition de la fonction qui calcule l'angle d'incidence avant le rebond
double angleIncidence(double v0, double thetaOut, double alpha, double x0, double y0, double x) {
  double sign;
  double den = v0*v0*cos(thetaOut)*cos(thetaOut);
  double derive = -GN*(x-x0)/den + tan(thetaOut);
  
  if (fabs(x-x0)<0.0001)   //pour eviter une division par zero
    sign = 1.0;
  else
    sign = (x-x0)/fabs(x-x0);
    
  double reponse = atan2(derive,sign);   //l'angle est entre -Pi et Pi

  return reponse;
}