/******************************************************************************
** [laplace.c]
**	 ラプラス方程式を差分法により解きます。
**
**	 使い方： laplace.exe > result.txt
******************************************************************************/
#include <stdio.h>
#include <math.h>

/*
** 計算領域の分割数 NX, NY を決めます。理論上は分割数を細かくすれば
** するほど計算誤差は減るはずですが、コンピュータが扱う浮動小数点数
** は有限値なので、あまり細かくしすぎると、桁落ちによって逆に
** 計算精度が落ちると思います。いろいろな DX, DY で結果を比較して、
** 計算精度が充分に良いところで、かつ計算速度も妥当な分割数を選ぶと
** 良いと思います。ここでは仮に 41 x 41 とします。
** NX と NY は同じ数にした方が、微分方程式の特性に合うらしいです。
*/

#define NX 41
#define NY 41


/*
** for (y のループ) { for (x のループ) { (x 方向の計算); } }
** とする場合は、Ｃ言語の配列のメモリ配置を考慮して、T(x, y)
** ではなく、T(y, x) とします。
*/
double T[NY][NX];


/*-----------------------------------------------------------------------------
** [SetBoundaryCondition]
** 境界条件を設定します。
** ２次元配列の扱い方はいろいろあります。書籍や web page 等を
** あたってみてください。
**---------------------------------------------------------------------------*/
void SetBoundaryCondition(double T[NY][NX])
{
	int x;
	int y;
	double left, top, right, bottom;

	fprintf(stderr, "left top right bottom\n");
	scanf_s("%lf%lf%lf%lf", &left, &top, &right, &bottom);

	/* 領域の底の境界値 */
	for (x = 0; x < NX; x++)
	{
		T[0][x] = bottom;
	}

	/* 領域の上の境界値 */
	for (x = 0; x < NX; x++)
	{
		T[NY-1][x] = top;
	}

	/* 左の壁の境界値 */
	for (y = 0; y < NY; y++)
	{
		T[y][0] = left;
	}

	/* 右の壁の境界値 */
	for (y = 0; y < NY; y++)
	{
		T[y][NX-1] = right;
	}
}


int main(void)
{
	const double dx = 1.0 / (NX-1);
	const double dy = 1.0 / (NY-1);

	/* 式(16)から、C1 と C2 の値を先に計算しておきます。 */
	const double C1 = 0.5 * dx*dx / (dx*dx + dy*dy);
	const double C2 = 0.5 * dy*dy / (dx*dx + dy*dy);

	int i;
	int j;

	double errorMax;
	double previousValue;

	/*
	** 境界条件を設定します。この境界条件は
	** ディリクレ条件なので、はじめに１度だけ設定
	** すればＯＫです。
	*/
	SetBoundaryCondition(T);

	/*
	** 前回の計算値と今回の計算値との最大誤差が
	** 0.00001 を下回るまで繰り返します。
	*/
	do
	{
		errorMax = 0.0;

		/* 計算領域は境界の１つ内側です。 */
		for (i = 1; i < NY-1; i++)
		{
			for (j = 1; j < NX-1; j++)
			{
				previousValue = T[i][j];

				/* ガウス・ザイデル法 */
				T[i][j] = C1*(T[i+1][j] + T[i-1][j]) + C2*(T[i][j+1] + T[i][j-1]);

				if (errorMax < fabs(T[i][j] - previousValue))
				{
					errorMax = fabs(T[i][j] - previousValue);
				}
			}
		}
	} while (errorMax > 0.00001);

	/* 全領域の計算結果を出力します。 */
	for (i = 0; i < NY; i++)
	{
		for (j = 0; j < NX; j++)
		{
			printf("%d\t%d\t%15.15f\n", j, i, T[i][j]);
		}
		printf("\n");
	}

	return 0;
}

// gnuplot>cd 'C:\hoge'
// gnuplot>set pm3d
// gnuplot>splot 'result.txt'