//******************************************************************************
// Programme de visualisation du système différentiel chaotique de Lorenz
//
// Dominique Lefebvre - TangenteX.com  Novembre 2009
//******************************************************************************

#include <cstdlib>
#include <iostream>
#include <math.h>
#include <dislin_d.h>

using namespace std;

//******************************************************************************
// Définition des constantes
//******************************************************************************
#define N        3          // nombre d'équations du système différentiel
#define PAS      0.01       // pas d'intégration temporel
#define STEP     10000      // nb de pas de calcul

//******************************************************************************
// Définition des variables globales
//******************************************************************************
double sigma, b, r, dt;

//******************************************************************************
// Définition des subroutines de service
//******************************************************************************
void InitGraphDislin(void)
{
    metafl("XWIN");
    disini();
    pagera();
    pagfll(255);
    color("black");
    winfnt("arial");
    
    return;     
}

//******************************************************************************
// Définition du système différentiel de Lorenz
// dx/dt = sigma*(y - x)
// dy/dt = r*x - y -x*z
// dz/dt = x*y - b*z
//******************************************************************************
void Derivee(double Y[], double dY[])
{
   dY[0] = sigma*(Y[1] - Y[0]);           // dx/dt
   dY[1] = r*Y[0] - Y[1] - Y[0]*Y[2];     // dy/dt
   dY[2] = Y[0]*Y[1] - b*Y[2];            // dz/dt
   return;
}

//******************************************************************************
// Routine RK4
//******************************************************************************
void RK4(double yin[], double yout[], double h)
{
   // La constante N, qui indique le nombre d'équations du système différentiel 
   // est définie globalement pour tout le programme.
   int i;
   double dydx[N], dydxt[N], yt[N], k1[N], k2[N], k3[N], k4[N];
  
   Derivee(yin, dydx); 
   for (i = 0; i < N; i++) 
    {
      k1[i] = h*dydx[i];
      yt[i] = yin[i] + 0.5*k1[i];
    }

   Derivee(yt, dydxt); 
   for (i = 0; i < N; i++)
    {
      k2[i] = h*dydxt[i];
      yt[i] = yin[i] + 0.5*k2[i];
    }   

   Derivee(yt, dydxt); 
   for (i = 0; i < N; i++)
    {
      k3[i] = h*dydxt[i];
      yt[i] = yin[i] + k3[i];
    }

   Derivee(yt, dydxt); 
   for (i = 0; i < N; i++)
    {
      k4[i] = h*dydxt[i];
      yout[i] = yin[i] + k1[i]/6. + k2[i]/3. + k3[i]/3. + k4[i]/6.;
    }
  return;
}


//******************************************************************************
// Programme principal
//******************************************************************************
int main(int argc, char *argv[])
{
    double x0, y0, z0, x , y, z, dt;
    double Yin[N], Yout[N];
    double X[STEP], Y[STEP], Z[STEP];
    int i,nbstep;
    char buffer[80];
    
    // Déclaration des paramètres du système différentiel
    // r est choisi par l'utilisateur
    sigma = 10.0;
    b = 8./3.;

    // Saisie des paramètres
    printf("Trace du papillon de Lorenz\n\n");
    printf("Parametre r: ");scanf("%lf",&r);
    printf("x0: "); scanf("%lf",&x0); 
    printf("y0: ");scanf("%lf",&y0); 
    printf("z0: "); scanf("%lf",&z0);

    // Initialisation des variables
    dt = PAS;
    nbstep = STEP;
    X[0] = x0; Y[0] = y0; Z[0] = z0;
   
    // Initialisation de la page graphique
    InitGraphDislin();
    wintit("Papillon de Lorenz - Espace de phase");
    name("Axe des x","X");
    name("Axe des y","Y");
    name("Axe des z","Z");
    labl3d("HORIZONTAL");
    height(30);           // hauteur des caractères 
    hname(30);            // hauteur des labels d'axe
    htitle(40);           // hauteur du titre
    titlin("Papillon de Lorenz",1);
    graf3d(-20.,20.,-20.0,10.,-50.,50.,-50.,25.,0.0,50.,0.,10.);
    box3d();
    title();

    // Affichage des paramètres dans le header de la page graphique
    height(40);
    sprintf(buffer," -  Parametres r: %.3lf  b: %.3lf  sigma: %.3lf",r,b,sigma);
    paghdr("Dominique Lefebvre - ",buffer,2,0); 
    

    // boucle de calcul et d'affichage de la trajectoire par RK4
    for (i=0;i<nbstep;i++)
       {
          Yin[0] = X[i];
          Yin[1] = Y[i];
          Yin[2] = Z[i];
          RK4(Yin, Yout, dt);
          X[i+1] = Yout[0];
          Y[i+1] = Yout[1];
          Z[i+1] = Yout[2];
        } 
    
    // tracé de la courbe 
    color("blue");     
    curv3d(X, Y, Z, nbstep);
    
    // Fin du programme
    disfin();
    return EXIT_SUCCESS;
}