/* 
  
     lorenz.c
     Estudo do sistema de lorenz
     dx/dt = sigma(y-x)
     dy/dt = -xz + rx -y
     dz/dt = xy -bz
        ->  trajetorias :  z X t para vários valores de r
        ->  projecao em 2D: z X x
        ->  secao de Poincare  z X x (quando a trajetoria
                                     cruza o plano y = 0)
     MODO = 0  trajetorias e projecao  (imprime todos os pontos) 
     MODO = 1  secao de Poincare ( imprime so quando y_anterior*y_posterior<0)
     

     a) "Projects in Computational Physics"
        Landau,Paez,Bordeianu   
	John Wiley and Sons, New York, 2007  
     b) "Computational Physics",
         N. Giordano, Prentice Hall, 1997                     
*/
     
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define N 3                                   /* Numero de equacoes */
#define dist 0.001                            /* passo de integracao */
#define MIN 0.0                               /*  tempo  minimo */
#define MAX 100.0                              /* tempo  maximo */
 
#define cteR 25.0                               /* r = 5, 10 ou 25 */
#define cteB 2.66666667                        /* b = 8/3  */
#define cteS 10.0                              /* sigma = 10  */

#define MODO 1                                 /* 0: z X t 
                                                  1: borboleta z X x
                                                  2: secao Poincare */

#define PI 3.141592653589793

/* o integrando - sistema de Lorenz */
double  f(double x, double y[], int i)
{
   if (i == 0) 
     return(cteS*(y[1]-y[0]));                    
   else if (i == 1) 
     return(-y[0]*y[2]+cteR*y[0]-y[1]); 
   else if (i == 2) 
     return(y[0]*y[1]-cteB*y[2]);        
   else{ 
     printf("Numero de eq. incorreto!\n");
     return -999999.;
   }

}
/* funcao do algoritmo de Runge-Kutta */
void runge4(double x, double y[], double passo)
{
   double h=passo/2.0,                        /* ponto medio */
          t1[N], t2[N], t3[N],                /*  variaveis temporarias */
          k1[N], k2[N], k3[N],k4[N];          /*  variaveis temporarias */
   int i;
 
   for (i=0; i<N; i++){
       k1[i]=passo*f(x, y, i);
       t1[i] = y[i]+0.5*k1[i];
       }
   for (i=0; i<N; i++){
       k2[i]=passo*f(x+h, t1, i);
       t2[i] = y[i]+0.5*k2[i];
       }
   for (i=0; i<N; i++) {
       k3[i]=passo*f(x+h, t2, i);
       t3[i] = y[i]+ k3[i];
       }
   for (i=0; i<N; i++) k4[i] = passo*f(x + passo, t3, i);
      
   for (i=0; i<N; i++) y[i] += (k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;
}
/*--------------------------------------------------------------------*/


int main()
{
 
   double x, y[N];                          /* */
   double y_ant;
   
   FILE *output;                          /* arquivo de dados de saida  */
   FILE *fpout;                           /* arquivo com macro p/ o gnuplot */

   output = fopen("lorenz.dat","w");
   
   x = 0;                                 /* inicializa o tempo */   

   /* condicoes iniciais */
   y[0] = 1.0;              
   y[1] = 0.0;
   y[2] = 0.0;

   for(x = MIN; x <= MAX ; x += dist)
   {  
      y_ant = y[1];
      runge4(x, y, dist);
      if (MODO==2){   /* secao de Poincare */
	if (y_ant*y[1]<0) {  /* fura o plano y=0 */
            fprintf(output, "%8.6f %10.8f %10.8f %10.8f\n", x, y[0], y[1], y[2]); 
            y_ant = y[1];
         }
      } 
      else {  /* Modo 0 e 1, guarda todos os pontos */
            fprintf(output, "%8.6f %10.8f %10.8f %10.8f\n", x, y[0], y[1], y[2]); 
      } 
       
   }
   fclose(output);


   /* plots */

  fpout=fopen("mplot.plt", "w+");   
  fprintf(fpout, "set term post eps enhanced color\n");
  fprintf(fpout, "set out 'Grafico.eps'\n");

  if (MODO==0) {
     fprintf(fpout, "set title ' Z(t) contra t '\n");
     fprintf(fpout, "set xlabel 'Tempo'\n");
     fprintf(fpout, "set ylabel 'Z(t) '\n");
     fprintf(fpout, "set key left top\n");
     fprintf(fpout, "plot 'lorenz.dat'  using 1:4 w lines lw 3 title 'r = 25' \n"); 
  }
  if (MODO ==1) {
     fprintf(fpout, "set title ' Z(t) contra X(t) '\n");
     fprintf(fpout, "set xlabel 'X'\n");
     fprintf(fpout, "set ylabel 'Z '\n");
     fprintf(fpout, "set key left top\n");
     fprintf(fpout, "plot 'lorenz.dat'  using 2:4 w lines lw 3 title 'r = 25' \n"); 
  }
  if (MODO ==2) {
     fprintf(fpout, "set title ' Secao de Poincare '\n");
     fprintf(fpout, "set xlabel 'X'\n");
     fprintf(fpout, "set ylabel 'Z '\n");
     fprintf(fpout, "set key left top\n");
     fprintf(fpout, "plot 'lorenz.dat'  using 2:4 pt 7  title 'r = 25' \n"); 
  }

  fclose(fpout);  
  system("gnuplot  'mplot.plt'");
  system("evince Grafico.eps");
  system("rm mplot.plt");
  system("rm lorenz.dat");


   return 0;
}