/* Steady 2-dimensional groundwater flow

        uniform, anisotopic hydraulic conductivity
        changing 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
              Number of effective grids is NUMBER_OF_HORIZONTAL_GRID-2
                                           NUMVER_OF_VERTICAL_GRID-1
                              (i=0 is water table not imaginary grid)

             0 1 2             NUMBER_OF_HORIZONTAL_GRID-1
             *-o-o-..........o-*

 */

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

#define         NUMBER_OF_HORIZONTAL_GRID       102
#define         NUMBER_OF_VERTICAL_GRID         11
#define         PI      3.141592653
#define         CYCLE   5.25                    /* sin関数の周期 */
#define         AMP     0.1                     /* sin関数の振幅 */

double          head[NUMBER_OF_VERTICAL_GRID][NUMBER_OF_HORIZONTAL_GRID];

void    main()
{
        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;
        double  Kx,Ky;
        double  a,b;
        double  dx,dy;
        double  rad;

/* parameters */

        Kx=100.0;
        Ky=1.0;
        dx=1.0;
        dy=1.0;

/* initial water table */

        water_table[0]=0.0;     /* imaginary grid */

        for(i=1;i<=NUMBER_OF_HORIZONTAL_GRID-2;i++){

                water_table[i]= (double)(i-1)*0.0101;

/*
                rad=(double)(2.0*(PI)*(CYCLE))/(double)(NUMBER_OF_HORIZONTAL_GRID-2)*(double)(i-1);
                water_table[i]= ((double)(i-1)*0.0092) + sin(rad)*(double)(AMP);
*/
        }
        water_table[NUMBER_OF_HORIZONTAL_GRID-1]=1.0;    /* imaginary grid */


/* judge for convergence */

        judge=0.0001;

/* over relaxation factor */

        omega=1.8;

/* uniform aquifer */

        a=Kx/(dx*dx);
        b=Ky/(dy*dy);

/* iteration */

        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=(b*head[i-1][j]+b*head[i+1][j]+a*head[i][j-1]+a*head[i][j+1])/(2*a+2*b);

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

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

                if(judge>sum_of_correction)     break;

        }/* wend */

        fprintf(stderr,"   ");
        for(j=1;j<NUMBER_OF_HORIZONTAL_GRID-1;j++){
                fprintf(stderr,"%6d",j);
        }
        fprintf(stderr,"\n");

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

                fprintf(stderr,"%2d:",i);
                for(j=1;j<NUMBER_OF_HORIZONTAL_GRID-1;j++){
                        printf("%6.3lf",head[i][j]);
                }
                printf("\n");
        }


}/* end of main */