/* 
 
  rk4.c: Integracao do oscilador harmonico simples usando Runge Kutta de
         quarta ordem
								       
  adaptado de  "Projects in Computational Physics" by Landau and Paez   
	        John Wiley and Sons, New York     
                  
                                                                     
*/

       
#include <stdio.h>

#define N 2                                   /* Numero de equacoes */
#define dt 0.1                                /* passo de integracao */
#define MIN 0.0                               /* tempo (x) minimo */
#define MAX 60.0                              /* tempo (x) maximo */
#define omega2 1                              /* T = 2*pi = 6.28 s */
 
/* funcao que diz o que estamos integrando - oscilador harmonico */
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*y[0]);           /* 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];                          /* */
   
   FILE *output;                              /* save data in rk4.dat */
   output = fopen("rk4.dat","w");
   
   y[0] = 1.0;                                /* initial position  */
   y[1] = 0.0;                                /* initial velocity  */
   x = MIN; /* instante 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 += dt)
   {
      runge4(x, y, dt);
      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]));   /* position vs. time */
   }
   printf("data stored in rk4.dat\n");
   fclose(output);
   return 0;

}
/*-----------------------end of main program--------------------------*/