/* 
 
  euler.c: Integracao do oscilador harmonico simples usando o alg. de Euler
								       
  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 (segundos) */
#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)
{                                                /* osc. harm. x nao eh relevante */
   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 de Euler */
void euler(double x, double y[], double passo)    
{
   int i;
   for (i=0; i<N; i++){
    y[i] = y[i] + passo*f(x,y,i);    
   }
}
/*--------------------------------------------------------------------*/


int main()
{   
   double x, y[N];                          /* x = tempo, y[1] = pos, y[2]=vel */

   FILE *output;                            /* armazena os dados em euler.dat */
   output = fopen("euler.dat","w");
   
   y[0] = 1.0;                                /*  posicao inicial  */
   y[1] = 0.0;                                /*  velocidade inicial  */
   x = MIN;  /* instante inicial */

   /* Imprime valores iniciais do tempo, posicao, velocidade, e energia mecanica (k = m = 1)*/
   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)
   {
     euler(x, y, dt); /* atualiza valores de posicao e velocidade */

      /* imprime tempo, posicao e energia mecanica */
      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 euler.dat\n");
   fclose(output);

   return 0;

}