#include <iostream>
#include <cmath>
#include <vector>

using namespace std;

int nc = 25;
int width = 8, height = 6, r = 3;
double h = 0.25;

double *gauss(double **matr, double *y, int n)
{
	double d, s;
	double *x = new double [n], *b = new double [n];
	double **a = new double *[n];
 	for (int i = 0; i < n; i++)
 		a[i] = new double [n];
 	for (int i = 0; i < n; i++)
 		for (int j = 0; j < n; j++)
 			a[i][j] = matr[i][j];
 	for (int i = 0; i < n; i++)
 		b[i] = y[i];

 	for (int k = 1; k < n; k++)
 		for (int j = k + 1; j < n; j++)
 		{
 			d = a[j][k] / a[k][k];
 			for (int i = k; i < n; i++)
 				a[j][i] = a[j][i] - d * a[k][i];
 			b[j] = b[j] - d * b[k];
 		}

 	for (int k = n - 1; k >= 1; k--)
 	{
 		d = 0;
 		for (int j = k + 1; j < n; j++)
 		{
 			s = a[k][j] * x[j];
 			d = d + s;
 		}
 	x[k] = (b[k] - d) / a[k][k];
 	}
 	delete b;
 		for (int i = 0; i < n; i++)
 			delete [] a[i];
 	delete [] a;
 	return x;
}

int main() {
 	int ht_h = height / h + 1, wh_h = width / h + 1, r_h = r / h + 1;
 	int n = wh_h * ht_h;
 	int i_0 = r_h, j_0 = wh_h - 1 - r_h;
 	double *y = new double [n], *x = new double [n];
 	double **a = new double *[n];
 	for (int i = 0; i < n; i++)
		a[i] = new double [n];
 	FILE *gp = popen("gnuplot -persist", "w");
 	if (gp == NULL)
 	{
 		printf("Ошибка открытия трубы gnuplot\n");
 		return 3;
 	}
 	for (int i = 0; i < ht_h; i++)
 		for (int j = 0; j < wh_h; j++)
 		{
 			int y_coor = i_0 - i;
 			int x_coor = j - j_0;

 		if (j == 0)
 		{
 			a[i * wh_h + j][i * wh_h + j] = 1;
 			y[i * wh_h + j] = 0;
 		} else if (i == ht_h - 1)
 		{
 			a[i * wh_h + j][i * wh_h + j] = -1 / h - 1;
 			a[i * wh_h + j][(i - 1) * wh_h + j] = 1 / h;
 			y[i * wh_h + j] = 0;
 		} else if ( i < i_0 && j > j_0 && ((pow(y_coor + 1, 2) + pow(x_coor, 2) >= pow(r_h, 2)
|| pow(y_coor, 2) + pow(x_coor + 1, 2) >= pow(r_h, 2)) && pow(y_coor, 2) + pow(x_coor, 2) < pow(r_h,2)))
 		{
 			double l, m;
 			if ((l = (sqrt(r_h * r_h - pow(x_coor, 2)) - y_coor)) > 1)
				l = 1;
 			if ((m = (sqrt(r_h * r_h - pow(y_coor, 2)) - x_coor)) > 1)
				m = 1;
 			a[i * wh_h + j][i * wh_h + j] = -1 /(m * h * h) - 1 / (1 * h * h) - 1;
 			a[i * wh_h + j][i * wh_h + j + 1] = 1 / (m * (m + 1) * h * h);
 			a[i * wh_h + j][i * wh_h + j - 1] = 1 / ((m + 1) * h * h);
 			a[i * wh_h + j][(i + 1) * wh_h + j] = 1 / ((l + 1) * h * h);
 			a[i * wh_h + j][(i - 1) * wh_h + j] = 1 / (l *(l + 1) * h * h);
 		}
 		else if (i < i_0 && j > j_0 && ((pow(y_coor - 1, 2) + pow(x_coor, 2) < pow(r_h, 2) ||
 pow(y_coor, 2) + pow(x_coor - 1, 2) < pow(r_h, 2)) &&  pow(y_coor, 2) + pow(x_coor, 2) >= pow(r_h, 2)))
 		{
 			a[i * wh_h + j][i * wh_h + j] = 1;
 			y[i * wh_h + j] = 200;
 		}
 		else if (i < i_0 && j > j_0 && (pow(y_coor - 1, 2) + pow(x_coor, 2) >= pow(r_h, 2) &&
 pow(y_coor, 2) + pow(x_coor - 1, 2) >= pow(r_h, 2) && pow(y_coor, 2) + pow(x_coor, 2) > pow(r, 2)))
 		{
 			a[i * wh_h + j][i * wh_h + j] = 1;
 			y[i * wh_h + j] = 0;
 		}
 		else if (j == wh_h - 1)
 		{
 			a[i * wh_h + j][i * wh_h + j] = 1;
 			y[i * wh_h + j] = 200;
		}
		else if (i == 0)
 		{
 			a[i * wh_h + j][i * wh_h + j] = 1;
 			y[i * wh_h + j] = 200;
 		}
 		else
 		{
 			a[i * wh_h + j][i * wh_h + j] = -2 / (h * h) - 2 / (h * h) - 1;
 			a[i * wh_h + j][i * wh_h + j + 1] = 1 / (h * h);
 			a[i * wh_h + j][i * wh_h + j - 1] = 1 / (h * h);
 			a[i * wh_h + j][(i + 1) * wh_h + j] = 1 / (h * h);
 			a[i * wh_h + j][(i - 1) * wh_h + j] = 1 / (h * h);
 		}
	}

	for (int k = 1; k <= nc; k++)
 	{
 		for (int i = 0; i < ht_h; i++)
 			for (int j = 0; j < wh_h; j++)
 			{
 				int y_coor = i_0 - i;
 				int x_coor = j - j_0;
 				if ((j > 0 && i > 0 && j < wh_h - 1 && i < ht_h - 1 && (i >= i_0 || j <= j_0)) || (i < i_0 && j > j_0 && pow(y_coor, 2) + pow(x_coor, 2) < pow(r_h, 2)) )
 					y[i * wh_h + j] = -x[i * wh_h + j];
 			}
 		x = gauss(a, y, n);
 	}

 	for (int i = 0; i < ht_h; i++)
 	{
 		for (int j = 0; j < wh_h; j++)
 			cout << x[i * wh_h + j] << " ";
 	cout << endl;
 	}
 	fprintf(gp, "set pm3d map\n");
 	fprintf(gp, "set autoscale fix\n");
 	fprintf(gp, "set cbtics 25\n");
 	fprintf(gp, "splot '-' matrix\n");
 	for (int i = ht_h - 1; i >= 0; i--)
 	{
 		for (int j = 0; j < wh_h; j++)
 		{
			if (x[i * wh_h + j] > 400)
				x[i * wh_h + j] = 200;
			if (x[i * wh_h + j] > 300)
				x[i * wh_h + j] = 197;
			if (x[i * wh_h + j] > 200)
				x[i * wh_h + j] = 195;
 			fprintf(gp, "%lf ", x[i * wh_h + j]);
 		}
 		fprintf(gp, "\n");
 	}
 	fprintf(gp, "e\n");
 	fclose(gp);
 	return 0;
}