//pendule elastique
#include <iostream>
#include <cmath>
#include <ctime>
#include "dislin.h"     //bilbiotheque DISLIN

using namespace std;

int main() {
  const int STEPS = 10000;
  double dT = 0.005;    
  double angle[STEPS], longueur[STEPS], temps[STEPS];
  double vitA[STEPS], vitL[STEPS], accA[STEPS], accL[STEPS]; 
  double energie[STEPS];

  //parametres du pendule  
  double L0 = 1.;   //longueur a repos du ressort en metres
  double k = 2.;    //constante elastique du ressort en N/m
  double g = 9.81;  //acceleration gravitatinelle
  double m = 0.1;   //masse

  cout <<"**********************************************************\n";
  cout << "Nous avons un pendule avec un ressort au lieu d'une corde.\n"
       << "La longueur a repos du ressort est " << L0 << " m et\n" 
       << "la constante elastique est " << k << " N/m.\n"
       << "La masse est " << m << " kg.\n"
       << "On assume que l'oscillation demarre depuis l'etat de repos\n"
       << "a un certain angle et elongation initiales.\n\n";
  cout << "Quelle est l'angle initial du pendule ? ";
  cin >> angle[0];
  cout << "Et l'elongation initiale du ressort en metres ? ";
  cin >> longueur[0];

  int opt;
  do {
    cout << "Que type de graphique souhaitez-vous ?" << endl; 
    cout << "1 pour dessiner un graphique" << endl;
    cout << "2 pour une animation du mouvement du pendule" << endl;
    cin >> opt;   cout << endl;
  } while (opt!=1 && opt!=2);

  temps[0] = 0.;
  angle[0] *= M_PI/180.;
   
  //ces vecteurs vont etre utilises avec DISLIN
  double x[STEPS], y[STEPS];
  x[0] = sin(angle[0])*longueur[0];
  y[0] = -cos(angle[0])*longueur[0]; 
  
  vitA[0]=0.;
  vitL[0]=0.;
  accA[0] = -1.*sin(angle[0])*g/longueur[0]; 
  accL[0] = cos(angle[0])*g - k*(longueur[0]-L0)/m;
 
  energie[0] = m*g*y[0] + k*(longueur[0]-L0)*(longueur[0]-L0)/2.;
  cout << "L'energie initiale du systeme est : " << energie[0] << endl;

  //initialisation de DISLIN
  metafl("XWIN");   //XWIN ou PDF
  disini();
  
  //methode de Runge
  for (int i=0; i<STEPS-1; i++) {   
    vitA[i+1] = vitA[i] + accA[i]*dT;
    vitL[i+1] = vitL[i] + accL[i]*dT; 
    angle[i+1]= angle[i] + vitA[i+1]*dT;
    longueur[i+1]=longueur[i] + vitL[i+1]*dT;
    temps[i+1]=(i+1)*dT;
   
    accA[i+1] = -1.*sin(angle[i+1])*g/longueur[i+1]
              - 2.*(vitA[i+1]*vitL[i+1])/longueur[i+1];
    accL[i+1] = cos(angle[i+1])*g - k*(longueur[i+1]-L0)/m
              + longueur[i+1]*vitA[i+1]*vitA[i+1]; 
    
    energie[i+1] = m*vitL[i+1]*vitL[i+1]/2.
                 + m*longueur[i+1]*longueur[i+1]*vitA[i+1]*vitA[i+1]/2. 
                 + m*g*y[i+1] + k*(longueur[i+1]-L0)*(longueur[i+1]-L0)/2.;
  
    if (abs(longueur[i+1]<=1.0E-6)) {
      longueur[i+1]=1.0E-6;
      cout << "Attention au conditions initiales et aux parametres du pendule !" << endl;
    }

    //tableaux de 1 element pour DISLIN
    double xplot[1], yplot[1];
    xplot[0] = sin(angle[i+1])*longueur[i+1];
    yplot[0] = -cos(angle[i+1])*longueur[i+1];
    x[i+1] = xplot[0];
    y[i+1] = yplot[0];
 
    double lignex[2]={0., x[i+1]};
    double ligney[2]={0., y[i+1]};   
    //dessin du pendule
    if (opt==2) {
      erase();   //efface la fenetre du graphique
      color("WHITE");
      name("axe-X", "X");
      name("axe-Y", "Y"); 
      titlin("Animation du pendule elastique", 1);
      graf(-4.*L0, 4.*L0, -4.*L0, 1.*L0, -5.*L0, 1.*L0, -5.*L0, 1.*L0);
      title();

      incmrk(-1);
      color("RED");
      marker(21);       
      curve(xplot,yplot,1);
      incmrk(0);
      thkcrv(1);
      color("BLUE");
      curve(x, y, i);
      thkcrv(5);
      color("YELLOW");
      curve(lignex, ligney, 2);
      sendbf();
      endgrf();
    }
    	
  }   //fin de la boucle

  //dessin de la trajectoire  
  if (opt==1) {
    name("temps (s)", "X");
    name("energie, angle, elongation", "Y"); 
    titlin("pendule elastique", 1);
    graf(0., STEPS*dT, 0., 1000*dT, -200.*L0, 350.*L0, -200.*L0 , 100.*L0); 
    title();
    
    for (int i=0; i < STEPS; i++) {
      angle[i] *= 180/M_PI;
      longueur[i] *= 100.;
      energie[i] = energie[i]*10. -100.;
    }
    
    thkcrv(10);
    color("RED");
    curve(temps, longueur, STEPS);
    color("BLUE");
    curve(temps, angle, STEPS);
    color("GREEN");
    curve(temps, energie, STEPS);
  }

  disfin();     //fin de DISLIN
 
  return 0;
}