/*  
  eq_volterra.c: Integracao das eqs nao-lineares de Volterra usando Runge Kutta de
  quarta ordem
                                                                   
*/

       
#include <stdio.h>

#define N 2                                   /* Numero de equacoes */
#define dt 0.01                                /* passo de integracao (anos) */
#define MIN 0.0                               /* tempo (x) minimo */
#define MAX 20.0                              /* tempo (x) maximo */


/* funcao que diz o que estamos integrando - equacao de Volterra */
double  f(double x, double y[], int i)
{
   if (i == 0) 
     return( y[0] - 0.01*y[0]*y[1]);           /* lado direito da primeira eq */
   else if (i == 1) 
     return(-y[1] + 0.01*y[0]*y[1]);           /* 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];                          /* tempo (x) e populacao de coelhos (y[0]) e raposas (y[1]) */
   
   FILE *output;                              /* escreve dados em volterra.dat */
   output = fopen("volterra.dat","w");

   /* Solicita ao usuario as populacoes iniciais de coelhos e raposas */   
   printf("Entre com as populacoes inicias de coelhos (c(t=0)) e de raposas (r(t=0)):");
   scanf("%lf %lf",&y[0],&y[1]);
   
   fprintf(output, "%8.6f %10.8f %10.8f\n", x, y[0], y[1]); /* imprime valores iniciais no arquivo de saida */

   /* laco for sobre a variavel tempo */
   for(x = MIN; x <= MAX ; x += dt)
   {
      runge4(x, y, dt);
      fprintf(output, "%8.6f %10.8f %10.8f\n", x, y[0], y[1]);   /* tempo, pop. coelhos, pop. raposas */
   }
   printf("dados escritos em volterra.dat\n");
   fclose(output);
   return 0;

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