//******************************************************************************
//* Programme de simulation du frottement d'un archer sur une corde de violon
//* Dominique Lefebvre  juin 2010
//* www.tangenteX.com
//******************************************************************************
#include <cstdlib>
#include <iostream>
#include <stdio.h>
#include <math.h>
#include <dislin_d.h>


using namespace std;

// Déclaration des paramètres
#define PI 3.141592

//******************************************************************************
// calcul du min et max de la table
//******************************************************************************
void MinMax(const double *Tableau, int n, double *Min, double *Max)
{
   int i;
   double mini, maxi;
   
   maxi = 0;
   mini = 1;
   for (i=0;i<n;i++)
     {
       if (Tableau[i] < mini) mini = Tableau[i];
       if (Tableau[i] > maxi) maxi = Tableau[i];              
     }  
   *Min = mini;
   *Max = maxi;  
}


//******************************************************************************
// Programme principal
//******************************************************************************
int main(int argc, char *argv[])
{
  // déclaration des variables locales
  double g,m,k,f,f0,v;
  double Omega0,Phi,A,t1,t2,t3,T;
  double x0,t0,dt,t, xt2;
  int i, npas;
  double *x, MinX, MaxX;


  // constantes et paramètres
  g = 9.81;

  // Saisie des paramètres de la simulation
  cout << "Masse d'excitation (m): "; cin >> m;
  cout << "Constante de raideur de la corde (k): "; cin >> k;
  cout << "Coefficient de frottement cinetique f: "; cin >> f;
  cout << "Coefficient de frottement statique f0: "; cin >> f0;
  cout << "Vitesse de l'archet v:"; cin >> v;

  // calcul des constantes du modèle
  Omega0 = sqrt(k/m);
  printf("Omega0: %f\n", Omega0);
  Phi = atan((Omega0*v)/(g*(f - f0)));
  printf("Phi: %f\n", Phi);
  A = -v/(Omega0*sin(Phi));
  printf("A: %f\n", A);

  // calcul des temps caractéristiques
  t1 = f0*m*g/k*v;
  printf("t1: %f\n",t1);
  t2 = t1 + (PI - 2*Phi)/Omega0;
  printf("t2: %f\n",t2);
  t3 = 2*g*(f0 - f)/(v*Omega0*Omega0) + t2;
  printf("t3: %f\n",t3);

  // calcul de la période 
  T = t3 - t1;

  // Conditions initiales
  t0 = 0;
  x0 = 0;
  dt = 0.00001;
  npas = (int)(t3/dt)+1;
  x = (double *)malloc(npas*sizeof(double));

  // calcul de la trajectoire  de la masse M
  // si t0 <= t <= t1  x(t) = v*t + x0
  // si t1 <= t <= t2  x(t) = A*cos(omega0(t - t1) + phi) + gf/omega0² + x0
  // si t2 <= t <= t3  x(t) = v(t-t2)+x(t2)
  t = t0;
  i = 0;
  while (t <= t1)
    {
	  x[i] = v*t + x0;
	  i++;
	  t += dt;
    }

   while ((t > t1) && (t <= t2))
    {
	  x[i] = A*cos(Omega0*(t - t1) + Phi) + g*f/(Omega0*Omega0) + x0;
	  i++;
	  t += dt;
    }
   xt2 = A*cos(Omega0*(t2 - t1) + Phi) + g*f/(Omega0*Omega0) + x0;

   while ((t > t2) && (t <= t3))  
    {
	  x[i] = v*(t-t2) + xt2;
	  i++;
	  t += dt;
    }

    MinMax(x,npas,&MinX,&MaxX);

    // Initialisation graphique
    metafl("XWIN");
    disini();
    pagera();           // dessine un contour au bord de la page
    pagfll(255);
    color("black");
    hwfont();           // type de font
    axspos(450,1800);
    axslen(2200,1200),
    
    name("Temps", "x");
    name("Deplacement", "y");
    
    titlin("Frottement de l'archet sur une corde",1);
    
    graf(0.0f,0.21f,0.0f,0.05f,MinX*1000,MaxX*1000,MinX*1000,0.01f);
    title();
    
    // tracé de la courbe
    t = t0;
    rlstrt(t,x0);
    t += dt;   
    for(i=0;i<npas;i++)
    {
      rlconn(t,x[i]*1000);
      t += dt;
    }
    
    // fin du programme
    disfin();
	free(x);
}