/* 
  
  pincare.c: Integracao do oscilador harmonico caotico 
              usando Runge Kutta de  quarta ordem
             construcao das secoes de poincare 
								       
     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 60000.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 0.5                                /* Amplitude da forca excitadora */
#define DomRest 1                             /* flag logica:  angulo entre 0 e 2pi */
#define ang0 0.2                            /* angulo inicial */
#define vel0 0.0                            /* vel angular inicial */
#define nmin 20                             /* num min de periodos antes do inicio */
                                            /* do map de poincare  */

#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 direito 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];                          /* */
   int n; 
   
   FILE *output;                              /* saida  */
   output = fopen("poincare.dat","w");
   
   x = 0;                                    /* inicializa tempo */   
   n = 0;

   y[0] = ang0;                                /* posicao inicial  */
   y[1] = vel0;                                /* velocidade inicial  */
   /* conte numero de periodos */
   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;
      } 
      /* em fase: 2*pi = n * w_d * t  */
      if (fabs(x-n*2.0*PI/omegaE)< dist/2.0) {
        if (n>nmin) fprintf(output, "%8.6f %10.8f %10.8f \n", x, y[0], y[1]);
        n = n+1;
      }      
   }
   printf("dados escritos em poincare.dat\n");
   fclose(output);
   return 0;
}