/* 
  
  caos.c: Integracao do oscilador harmonico caotico 
              usando Runge Kutta de  quarta ordem
								       
     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 <math.h>

#define N 2                                   /* Numero de equacoes */
#define dist 0.001                            /* passo de integracao */
#define MIN 0.0                               /* x minimo */
#define MAX 160.0                              /* x maximo */
#define omega2 1.0                            /* freq natural do sist */
#define q 0.5                                 /* const. de amortecimento */
#define omegaE 0.666666667                    /* freq da forca excitadora */ 
#define FE 1.2                                /* Amplitude da forca excitadora */
#define DomRest 0                             /* flag logica:  angulo entre 0 e 2pi */
#define ang0 0.2                            /* angulo inicial */
#define vel0 0.0                            /* vel angular inicial */

#define PI 3.141592653589793

/* funcao que diz o que estamos integrando - oscilador nao linear forcado */
double  f(double x, double y[], int i)
{
   if (i == 0) 
     return(y[1]);                   /* lado direito da primeira eq */
   else if (i == 1) 
     return(-omega2*sin(y[0]) -q*y[1] +
                      FE*sin(omegaE*x));           /* lado esquerdo da segunda eq */
   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 do meio */
          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];                          /* */
   
   FILE *output;                              /* saida em rk4.dat */
   output = fopen("caos.dat","w");
   
   x = 0;                                    /* inicializa tempo */   

   y[0] = ang0;                                /* posicao inicial  */
   y[1] = vel0;                                /* velocidade inicial  */
   

   fprintf(output, "%8.6f %10.8f %10.8f %10.8f\n", x, y[0], y[1], 
           0.5*(y[1]*y[1]+y[0]*y[0])); 
    

   for(x = MIN; x <= MAX ; x += dist)
   {
      runge4(x, y, dist);
      if(DomRest) {
         if (y[0]>PI)  y[0] = y[0] - 2.0*PI;
         if (y[0]<-PI) y[0] = y[0] + 2.0*PI;
      } 
      fprintf(output, "%8.6f %10.8f %10.8f %10.8f\n", x, y[0], y[1], 
           0.5*(y[1]*y[1]+y[0]*y[0]));       
   }
   printf("dados escritos em caos.dat\n");
   fclose(output);
   return 0;
}