/*
	Author: Alvin R. Lebeck
	Date: 5/19/92
	Description:
		Meshgenerations with thompsons solver
		modified serial version of W.Gentzsch Fh Regensburg FRG
*/

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

	/* SIZE has to be one greater than what the fortran N would be*/
#ifndef SIZE
#define SIZE 1026
#endif
#define N (SIZE-1)

#define ALFA 0.1
#define RELFA 0.98
#define EPS 0.5e-7
#define H 1.0/(N-1)
#define REL 2.0/RELFA

/* some macros for an optimization for CPROF */

/*
#define xa(i,j) (xy[2*i][j]) 
#define ya(i,j) (xy[(2*i)+1][j]) 
*/
#define xa(i,j) (x[i][j]) 
#define ya(i,j) (y[i][j]) 

int LMAX = 100;

/*
double aa[SIZE][SIZE],xy[2*SIZE][SIZE],rx[SIZE][SIZE],ry[SIZE][SIZE];
*/
double aa[SIZE][SIZE],x[SIZE][SIZE],y[SIZE][SIZE],rx[SIZE][SIZE],ry[SIZE][SIZE];
double d[SIZE][SIZE],dd[SIZE][SIZE];

main(argc,argv)
int argc;
char **argv;
{
	int i,j,k,i1p,j1p,i2m,j2m,ll,ixcm,jxcm,iycm,jycm,irxm,jrxm,irym,jrym;	
	int ip,im,jp,jm,l,m;
	double dxcm,dycm,rxm,rym,r,dmax;
	double xx,yx,xy,yy,a,b,c,qi,qj,pxx,qxx,pyy,qyy,pxy,qxy;

	/* initialize x and y with FACTOR ALFA */


	if (argc == 2)
		LMAX = atoi(argv[1]);

	printf("LMAX = %d\n",LMAX);

	for (i = 1; i <= N; ++i)
	{
		xa(i,1) = (float) (i - 1) / (float) (N-1);
	}

	for (i = 1; i <= N; ++i)
	{
		xa(i,N) = xa(i,1);
		xa(1,i) = 0.0;
		xa(N,i) = 1.0;
		ya(i,1) = 0.0;
		ya(1,i) = xa(i,1);
		ya(N,i) = ALFA*xa(i,1);
	}

	for (i = 1; i <= N; ++i)
	{
		ya(i,N) = (1.0 - xa(i,1)) * ya(1,N) + xa(i,1) * ya(N,N);
	}

	for (i = 2; i <= N-1; ++i)
		for (j = 2; j <= N-1; ++j)
		{
			xa(i,j) = 0.9 * xa(i,1);
			ya(i,j) = 0.9*((1.0-xa(i,1))*ya(1,j)+xa(i,1)*ya(N,j));
		}
	
	printf("*** 2-D ITERATION BEHAVIOR ***\n");
	printf("    I   J   X-COR      I   J   Y-COR      I   J  X-RES      I   J   Y-RES\n");

	i1p = 2;
	j1p = 2;
	i2m = N-1;
	j2m = N-1;

	ll = 0;
	dmax = 1.0;

	while ((ll < LMAX) && (dmax > EPS))
	{
		ixcm = jxcm = iycm = jycm = irxm = jrxm = irym = jrym = 0;
		dxcm = dycm = rxm = rym = 0.0;

		/* if first time through don't add corrections */
		/* else update the entire first row */

		if (ll > 0)
		{
			l = 0;
			for (j = j1p; j < j2m; ++j)
			{
				++l;
				xa(1,j) = xa(1,j) + rx[1][l];
				ya(1,j) = ya(1,j) + ry[1][l];
			}
		}

		/* i loop is outer in C */
		for (i = i1p; i <= i2m; ++i)
		{
			ip = i+1;
			im = i-1;

			/* if not first time through add */
			/* corrections to the new column */
			if (ll > 0)
			{
				if (i == i1p)
				{
					l = 0;
					for (j = j1p; j <= j2m; ++j)
					{
						++l;
						xa(i,j) = xa(i,j) + rx[i][l];
						ya(i,j) = ya(i,j) + ry[i][l];
					}
				} 

				l = 0;
				for (j = j1p; j <= j2m; ++j)
				{
					++l;
					xa(ip,j) = xa(ip,j) + rx[ip][l];
					ya(ip,j) = ya(ip,j) + ry[ip][l];
				}
			}

			/* j is inner loop in C */
			m = 0;
			for (j = j1p; j <= j2m; ++j)
			{
				m++;
				jp = j+1;
				jm = j-1;
				xx = xa(ip,j)-xa(im,j);
				yx = ya(ip,j)-ya(im,j);
				xy = xa(i,jp)-xa(i,jm);
				yy = ya(i,jp)-ya(i,jm);
				a  = 0.25 * (xy*xy+yy*yy);
				b  = 0.25*(xx*xx+yx*yx);
				c  = 0.125 *(xx*xy+yx*yy);
				qi = 0.0;
				qj = 0.0;

				/* as commented out in spec version
				qi = a*0.5;
				qj = b*0.5;
				*/

				aa[i][m] = -b;
				dd[i][m] = b+b+a*REL;
				pxx = xa(ip,j)-2.0*xa(i,j)+xa(im,j);
				qxx = ya(ip,j)-2.0*ya(i,j)+ya(im,j);
				pyy = xa(i,jp)-2.0*xa(i,j)+xa(i,jm);
				qyy = ya(i,jp)-2.0*ya(i,j)+ya(i,jm);
				pxy = xa(ip,jp)-xa(ip,jm)-xa(im,jp)+xa(im,jm);
				qxy = ya(ip,jp)-ya(ip,jm)-ya(im,jp)+ya(im,jm);
				
/*
     CALCULATE RESIDUALS (EQUIVALENT TO RIGHT HAND SIDES OF EQUS.)
*/

				rx[i][m] = a*pxx+b*pyy-c*pxy+xx*qi+xy*qj;
				ry[i][m] = a*qxx+b*qyy-c*qxy+yx*qi+yy*qj;
				
			}


			for (j = 1; j <= m; ++j)
			{
			    if(fabs(rx[i][j]) > fabs(rxm))
			    {
				rxm  = rx[i][j];
				irxm = i;
				jrxm = j;
			    }


			    if(fabs(ry[i][j]) > fabs(rym)) 
			    {
				rym  = ry[i][j];
				irym = i;
				jrym = j;
			    }
			}

			d[i][1] = 1.0/dd[i][1];

			for (j = 2; j <= m; ++j)
			{
				r = aa[i][j]*d[i][j-1];
				d[i][j] = 1.0/(dd[i][j]-aa[i][j-1]*r);
				rx[i][j] = rx[i][j] - rx[i][j-1]*r;
				ry[i][j] = ry[i][j] - ry[i][j-1]*r;
			}

			rx[i][m] = rx[i][m]*d[i][m];
			ry[i][m] = ry[i][m]*d[i][m];

		} /* end i loop */

		for (i = i1p; i <= i2m; ++i)
		{
			k = 0;
			for (j = 2; j <= m ; ++j)
			{
				k = m-j+1;
				rx[i][k] = (rx[i][k]-aa[i][k]*rx[i][k+1])*d[i][k];
				ry[i][k] = (ry[i][k]-aa[i][k]*ry[i][k+1])*d[i][k];
			}
		}

		++ll;
		dmax = (fabs(rxm) > fabs(rym)) ? fabs(rxm) : fabs(rym);
		printf("%d %d %d %0.3e %d %d %0.3e %d %d %0.3e %d %d %0.3e\n",
			ll,ixcm,jxcm,dxcm,iycm,jycm,dycm,irxm,jrxm,rxm,
			irym,jrym,rym);
	}
}