#include <stdlib.h>
#include <stdio.h>
#include <string.h>

int M, N; //размерность сетки
double *prev; //массивы значений в предыдущий момент времени
double **A;
double *x;
double T = 500.0; //граничное условие
// double C = 1.0;
// double R = 5.0;
// double h = 1.0;
int done = 0;
double sec = 20;
double d_time = 0.5;
double c = 1;
double den = 1;
double alph = 1;
double dx = 0.5;
double dy = 0.5;

//ВНИМАТЕЛЬНО ПРОЧИТАЙ ОПИСАНИЕ МЕТОДА ГАУССА

//Берет номер узла и записывает общее оравнение для него
void wrtk(int k, int i, int j) {
	double a = alph / c / den;
	A[k][(j + 1) * M + i] =  a / dy / dy;
	A[k][(j - 1) * M + i] = a / dy / dy;
	A[k][j * M + i + 1] =  a / dx / dx;
	A[k][j * M + i - 1] = a / dx / dx;
	A[k][j * M + i] =  - (1 / d_time) - (2 * a / dx / dx) - (2 * a / dy / dy);
	x[k] =  - prev[j * M + i] / d_time;
}

//Заполняет матрицу коэффициентов А и вектор b(x)
void mysolver(void) { //расчет значений
	int i, j;
	int k = 0;
	for (j = 0; j < N; j++) {
		for (i = 0; i < M; i++) {
			if (j == 0 ) {//нижняя граница
				A[k][j * M + i] = 1;
				A[k][(j + 1) * M + i] = -1;
			} else if (j == N - 1 && i > M / 2) {//Верхняя правая граница
				A[k][j * M + i] = 1;
				A[k][(j - 1) * M + i] = -1;
			} else if (j == N - 1) {//Верхняя правая
				A[k][j * M + i] = 1;
				x[k] = T;
			} else if (i == 0) {//левая
				A[k][j * M + i] = 1;
				x[k] = T;
			} else if (i == M - 1) {//Правая
				A[k][j * M + i] = 1;
				x[k] = T;
			} else {
				wrtk(k , i , j);
			}
			k++;
		}
	}
}


// Вот тут будь осторожнее. Нормальный метод Гаусса использует три массива A, b, x.
// Я использую только А и x. Массив х играет роль b. После решения, ответы записываются в x. Сделано только для экономии памяти. Можно было бы использовать
// и А и b и х. Но до обратного хода вектор х нахрен не нужен. А когда мы записываем ответ в x[i], b[i] нам уже больше не пригодится, поэтому для
// экономии туда можно записаь результат. Почему ипользую A и х, а не А и b, хз. Разницы нет.

void gauss(int n) {

	int i, j, k;
	double a;
	//Прямой ход метода Гаусса
	for ( k = 0; k < n; k++ )
		for ( j = k + 1; j < n; j++) {
			a = - A[j][k] / A[k][k];
			x[j] += x[k] * a;
			for ( i = k ; i < n; i++ )
				A[j][i] += A[k][i] * a;
		}

	//Обратный ход
	x[n - 1] /= A[n - 1][n - 1];
	for ( i = n - 2; i >= 0; i-- ) {
		for ( j = i + 1; j < n; j++ )
			x[i] -= x[j] * A[i] [j];
		x[i] /= A[i][i];
	}
}

int main(int argc, char* argv[]) {
	int n_time; //количество циклов
	int i, j, k;
	int nur;

	M = 12/dx;
	N = 4/dy;
	n_time = sec * d_time;
	nur = M * N;
	FILE *gp=popen("gnuplot -persist","w");

	prev = (double*) malloc(sizeof(double) * nur);
	A = (double**) malloc(sizeof(double*) * nur);
	x = (double*) malloc(sizeof(double) * nur);

	for (k = 0; k < nur; k++) {
		A[k] = (double*) malloc(sizeof(double) * nur);
	}

	memset(prev, 10, sizeof(double) * nur);//Установка начальной температуры


	for (k = 0; k < n_time; k++) {

		memset(x, 0, sizeof(double) * nur ); // Инициалицируем х и А нулями

		for (i = 0; i < nur; i++)
			memset(A[i], 0, sizeof(double) * nur);
		mysolver();

		gauss(nur);

		FILE *fd = fopen("output.txt", "w");
		memcpy(prev, x, sizeof(double) * nur);
		//Печать решения
		for (j = 0; j < N; j++) {
			for (i = 0; i < M; i++) {

			fprintf(fd, "%.4f ", prev[j * M + i]);
			}
		}
		fclose(fd);
		fprintf(gp, "unset key\n");
		fprintf(gp, "set tic scale 0\n");
		fprintf(gp, "set palette rgb 33,13,10\n");//color
		fprintf(gp, "set cbrange [0:100]\n");//-10000:10000
		fprintf(gp, "unset cbtics\n");
		fprintf(gp, "set xrange [-0.5:%f]\n",24-0.5);
		fprintf(gp, "set yrange [-0.5:%f]\n",8-0.5);
		fprintf(gp, "set pm3d interpolate 0,0\n");
		fprintf(gp, "set pm3d map\n");
		//fprintf(gp,"do for[i=0:%d] {\n",n);

		//fprintf(gp,"5 4 3 1 0\n2 2 0 0 1\n0 0 0 1 0\n0 0 0 2 3\n0 1 2 4 3\n");
		fprintf(gp, "splot \"output.txt\" matrix with pm3d\n");
		fflush(gp);
	}//k



	return 0;
}