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

/****************** 図7.6の数値解を求めるプログラム *********************

	・SOR法によりポアソン方程式を解きます．

*************************************************************************/

///////////////////////////////////
#define Pi	4.*atan(1.)
#define omg	1.0
#define LM	300
///////////////////////////////////
#define x0	(double)0
#define xM	(double)1
#define M	40
#define dx	(xM-x0)/(double)M
#define y0	(double)0
#define yM	(double)1
#define dy	(yM-y0)/(double)M
///////////////////////////////////
#define u_x0	(double)0
#define u_xM	(double)0
#define u_y0	(double)0
#define u_yM	(double)0
///////////////////////////////////

void initial_c(double u[][M+1]){
	int i, j;
	for(j=0; j<=M; j++){
		for(i=0; i<=M; i++){
			u[i][j] = 0.;
		}
	}
}

double fij(double x, double y){
	double z;
	z = sin(Pi*x)*sin(2.*Pi*y);
	return z;
}

double sij(double x, double y){
	double z;
	z = -5.*Pi*Pi*sin(Pi*x)*sin(2.*Pi*y);
	return z;
}

void cal_f(double x[], double y[], double f[][M+1]){
	int i, j;
	double xx, yy;
	for(j=0; j<=M; j++){
		for(i=0; i<=M; i++){
			xx = x[i];
			yy = y[j];
			f[i][j] = fij(xx, yy);
		}
	}
}

void cal_xy(double x[], double y[]){
	int i, j;
	for(i=0, j=0; i<=M || j<=M; i++, j++){
		x[i] = x0+(double)i*dx;
		y[j] = y0+(double)j*dy;
	}
}

void boundary_c(double u[][M+1]){
	int i, j;
	for(i=0, j=0; i<=M || j<=M; i++, j++){
		u[0][j] = u_x0;
		u[M][j] = u_xM;
		u[i][0] = u_y0;
		u[i][M] = u_yM;
	}
}

void solver_sor(double x[], double y[], double u[][M+1]){
	int i, j;
	double uli, suli;
	double xx, yy, ar, ddx, ddy, s;
	double sumx, fx, sumy, fy;
	ar = dy/(dx);
	ddx = dx*dx;
	ddy = dy*dy;
	uli = 0.;
	for(j=1; j<=M-1; j++){
		for(i=1; i<=M-1; i++){
			xx = x[i];
			yy = y[j];
			s = sij(xx, yy);
			uli = u[i][j];
			sumx = u[i-1][j]+u[i+1][j];
			fx = -2.;
			sumy = u[i][j-1]+u[i][j+1];
			fy = -2.;
			u[i][j] = omg*(ddy*s-ar*ar*sumx-sumy)
					/(ar*ar*fx+fy)
					+(1-omg)*u[i][j];
			uli -= u[i][j];
			suli += sqrt((uli/u[i][j])*(uli/u[i][j]));
		}
	}
}

int main(void){

	int l=0, j;
	double x[M+1];
	double y[M+1];
	double u[M+1][M+1];
	double f[M+1][M+1];

	initial_c(u);
	boundary_c(u);
	cal_xy(x, y);
	cal_f(x, y, f);

	while(l <= LM){

		solver_sor(x, y, u);

		if(l%10 == 0)
			for(j=0; j<=M; j++)
				printf("%5.3lf, %12.7lf\n", y[j], u[M/2][j]);

		l++;

	}

	return 0;

}