#include <stdio.h>
#include <stdlib.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);
float h(float x);
float gx0(float t);
float gxfim(float t);

int main(void)
{
    FILE *outf;

    float t, t0, tfim, x, x0, xfim,
      deltat, deltax, alfa, c, v;

    float **solucao;

    int i, j, nx, nt;


    outf = fopen("questao2.dat", "w");

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

    t0   = 0.0;
    tfim = 5.0;

    x0   = 0.0;
    xfim = 5.0;

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


    solucao = (float **) malloc(nt*sizeof(float*));
    for (i = 0; i < nt; i++)
    {
       solucao[i] = (float*) malloc (nx*sizeof(float));
    }

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

    /* Parametros fisicos e variavel auxiliar */

    v = 1;

    c = v * deltat/deltax;

    /* Condicao inicial */

    x = x0;

    for (j = 0; j < nx; j++)
    {
        solucao[0][j] = f(x);
        x += deltax;
    }

    /* Contorno */
    t = t0;
    for (i = 0; i < nt-1; i++)
    {
        t += deltat;
        solucao[i][0]      = gx0(t);
        solucao[i][nx - 1] = gxfim(t);
        /* Meio */
        x = x0;
        for (j = 1; j < nx - 1; j++)
        {
            x += deltax;

            if (i == 0)
                solucao[1][j] = (1 - c*c)*f(x) + (c*c/2)*(f(x+deltax)+f(x-deltax)) + deltat*h(x);
            else
                solucao[i+1][j] = 2*(1 - c*c)*solucao[i][j] + (c*c)*(solucao[i][j+1]+solucao[i][j-1]) - solucao[i-1][j];

        }
    }

    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;
    }

    for (i = 0; i < nt; i++)
    {
       free(solucao[i]);
    }
    free(solucao);

}

float f(float x)
/* Funcao que da' a condicao inicial */
{
    if (x <= 0 || x >= 5.0)
        return 0;
    return exp(-((x-2.5)*(x-2.5)/0.04));
}

float h(float x)
/* dfi/dt em t=0 */
{
    return(0.0);
}

float gx0(float t)
{
   return(0.0);
}

float gxfim(float t)
{
   return(0.0);
}