#include <stdio.h>
#include <math.h>

#define PI 3.1415926

/* Programa de demonstracao que implementa o metodo de resolucao de equacao */
/* de conveccao  por um metodo explicito */

float f(float x);
   
int main(void)
{	
FILE *outf;
	
float solucao[600][1200],
      t, t0, tfim, x, x0, xfim, xpulso, 
      deltat, deltax, alfa, c, v;

int i, j, nx, nt;

/* if ((*/ outf = fopen("conveccao.dat", "w"); /* ) = NULL) */
/*   { */
/*   printf("\nProblemas na abertura do arquivo\n"); */
/*   } */

deltat = 0.001; /* Discretizacao do tempo */
deltax = 0.01;   /* Discretizacao do espaco */

t0   = 0.0;
tfim = 0.5;

x0   = 0.0;
xfim = 10.0;

nt = (int) ((tfim - t0)/deltat);
nx = (int) ((xfim - x0)/deltax);

fprintf(stdout, " Numero de intervalos temporais %d e espaciais %d", nt, nx);

/* Parametros fisicos e variavel auxiliar */

v = 4;

c = v * deltat/deltax;

xpulso = (xfim - x0)/2.0;

/* Condicao inicial */

x = x0;

for (j = 1; j < nx - 2; j++)
   {
   solucao[0][j] = f(x - xpulso);

   x += deltax;
   }

t = t0;

for (i = 0; i < nt; i++)
   {
   t += deltat;
   
   for (j = 1; j < nx - 1; j++)
       {
       solucao[i + 1][j] = solucao[i][j]  - c * (solucao[i][j] - solucao[i][j - 1]);
       }
   }	   

x = x0;

for (j = 0; j < nx; j++)
   {
   fprintf(outf, "%f ", x);
   for (i = 0; i < nt; i++)
      {
      fprintf(outf, "%f ", solucao[i][j]);
      }
   fprintf(outf, "\n");
   x += deltax;
   }
   
}

float f(float x)
/* Funcao que da' a condicao inicial */
   {
    if (x >= 0.0) 
       if (x <= 1)  return(/*sin(PI * x)*/ 1.0);
       else         return(0.0);
    else            return(0.0);
    
   }