#include <stdlib.h>

#include <stdio.h>

#include <math.h>



FILE *output;



    int     numx = 100;

    int     numy = 100;

    int     numz = 100;

    int     numt = 20000;



    int     i = 0, j = 0, k = 0, m = 0;

    double  x = 0.0, y = 0.0, z = 0.0, t = 0.0;



    double  space_step = 20 / (numx);

    double  time_step = 0.002;



    double  C_coef = time_step / pow(space_step,2);





void boundary(){

    double  C[numx][numy][numz];



      for(k = 0; k < numz; k++){

        for(i = 0; i < numx; i++){

            C[i][0][k] = 0.1;  C[i][numy-1][k] = 0.1;

            C[0][i][k] = 0.1;  C[numx-1][i][k] = 0.1;

            for(j = 0; j < numy; j++){

                C[i][j][numz-1] = 0.1;}

        }

      }

}





int main (){

    double  C[numx][numy][numz];

    double  New_C[numx][numy][numz];

    double  C_diffuse;



    output = fopen("C_3d.txt", "w");



    for(i = 0; i < numx; i++)

    for(j = 0; j < numy; j++)

    for(k = 0; k < numz; k++){



        C[i][j][k] = 0;

        C[numx/2][numy/2][0] = 100;

    }





    for(m = 0; m < numt; m++){



        for(i = 1; i < numx-1; i++)

        for(j = 1; j < numy-1; j++)

        for(k = 1; k < numz-1; k++){



            C_diffuse = C_coef * (C[i+1][j][k] + C[i-1][j][k] +\

                                    C[i][j+1][k] + C[i][j-1][k] +\

                                    C[i][j][k+1] + C[i][j][k-1] - 6 * C[i][j][k]);



            New_C[i][j][k] = C[i][j][k] + C_diffuse;



        }



        for(i = 1; i < numx-1; i++)

        for(j = 1; j < numy-1; j++)

        for(k = 1; k < numz-1; k++){



            C[i][j][k] = New_C[i][j][k];



        }



        boundary();



        C[numx/2][numy/2][0] = 100;



    }



    for(i = 0; i < numx; i++)

    for(j = 0; j < numy; j++)

    for(k = 0; k < numz; k++){



        x = i * space_step;   y = j * space_step;   z = k * space_step;   t = m * time_step;



        fprintf(output, "%f\t%f\t%f\t%f\t%f\n", t, x, y, z, C[i][j][k]);

    }

    fflush(output);

    fclose(output);

    

    return 0;

}