//******************************************************************************
// Programme de visualisation du système différentiel chaotique de Rossler
//
// 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     60000      // nb de pas de calcul

//******************************************************************************
// Définition des variables globales
//******************************************************************************
double a, b, c, 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 Rossler
// dx/dt = -y - z
// dy/dt = x + ay
// dz/dt = b + z(x - c)
//******************************************************************************
void Derivee(double Y[], double dY[])
{
   dY[0] = -Y[1] - Y[2];           // dx/dt
   dY[1] = Y[0] + a*Y[1];          // dy/dt
   dY[2] = b + Y[2]*(Y[0] - c);    // 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
    a = 0.2;
    b = 0.2;
    c = 5.7;

    // Saisie des paramètres
    printf("Trace de l'attracteur de Rossler\n\n");
    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("Attracteur de Rossler - 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("Attracteur de Rossler",1);
    graf3d(-15.,15.,-15.0,5.,-15.,15.,-15.,5.,0.0,30.,0.,10.);
    box3d();
    title();

    // Affichage des paramètres dans le header de la page graphique
    height(40);
    sprintf(buffer," -  Parametres a: %.3lf  b: %.3lf  c: %.3lf",a,b,c);
    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;
}