/* Two-dimensional laplace equation -Laplace2.c

        isotopic hydraulic conductivity
        equivalent dx and dy

        vertical direction : i
        horizontal direction : j

        water table (constant head) i=0

        imaginary grid for impervious boundary condition
                i=NUMBER_OF_VERTICAL_GRID-1     : bottom
                j=0                             : left
                j=NUMVER_OF_HORIZONTAL_GRID-1   : right
        (note: there is no NUMVER_OF_HORIZONTAL_GRIDth grid)

 */

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

#define         NUMBER_OF_HORIZONTAL_GRID       13
#define         NUMBER_OF_VERTICAL_GRID         11

double          head[NUMBER_OF_VERTICAL_GRID][NUMBER_OF_HORIZONTAL_GRID];

void    main()
{
        FILE    *fp;
        int     i,j;
        int     number_of_horizontal_grid;
        int     number_of_vertical_grid;
        double  water_table[NUMBER_OF_HORIZONTAL_GRID];
        double  sum_of_correction, correction;
        double  old_head;
        double  judge;
        int     iteration;
        double  omega;          /* over relaxation factor */
        double  headij;

/* initial water table */

        water_table[0]=0.0;     /* imaginary grid */
        water_table[1]=0.0;
        water_table[2]=0.1;
        water_table[3]=0.2;
        water_table[4]=0.3;
        water_table[5]=0.4;
        water_table[6]=0.5;
        water_table[7]=0.6;
        water_table[8]=0.7;
        water_table[9]=0.8;
        water_table[10]=0.9;
        water_table[11]=1.0;
        water_table[12]=1.0;    /* imaginary grid */

/* judge for convergence */

        judge=0.0001;

/* over relaxation factor */

        omega=1.6;

/* iteration */
printf("Get into iteration\n");

        iteration=0;

        while(1){

                sum_of_correction=0.0;

                for(i=1;i<NUMBER_OF_VERTICAL_GRID-1;i++){              /* vertical scan */

                     /* left boundary */
                        head[i][0]=head[i][1];
                     /* right boundary */
                        head[i][NUMBER_OF_HORIZONTAL_GRID-1]=head[i][NUMBER_OF_HORIZONTAL_GRID-2];

                        for(j=1;j<NUMBER_OF_HORIZONTAL_GRID-1;j++){      /* horizontal scan */

                             /* upper boundary */
                                head[0][j]=water_table[j];
                             /* bottom boundary */
                                head[NUMBER_OF_VERTICAL_GRID-1][j]=head[NUMBER_OF_VERTICAL_GRID-2][j];

                                old_head=head[i][j];
                                headij=(head[i-1][j]+head[i+1][j]+head[i][j-1]+head[i][j+1])/4.0;

                                head[i][j]=old_head+omega*(headij-head[i][j]);

                                correction=fabs(old_head-head[i][j]);
                                sum_of_correction += correction;
                        }
                }
                printf("IT:%4d  Sum. of Cor.%lf\n",++iteration, sum_of_correction);

                if(judge>sum_of_correction)     break;

        }/* wend */

        printf("   ");
        for(j=0;j<NUMBER_OF_HORIZONTAL_GRID;j++){
                printf("%5d",j);
        }
        printf("\n");

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

                printf("%2d:",i);
                for(j=0;j<NUMBER_OF_HORIZONTAL_GRID;j++){
                        printf("%5.2lf",head[i][j]);
                }
                printf("\n");
        }


}/* end of main */