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

int M, N; //размерность сетки
int An;
double **A;
double *x;
double *b;
int *rep;
double c = 1;
double den = 1;
double alph = 1;
double dx = 0.3;
double dy = 0.3;


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

//Заполняет матрицу коэффициентов А и вектор b(x)
void mysolver(void) { //расчет значений
	int i, j;
	int k = 0;
	int id;
	for (j = 0; j < N; j++) {
		for (i = 0; i < M; i++) {
			id = j * M + i;
			if (i == 20 && j == 15)
			printf("%d\n",id );
			if (j >= An && i > M - j + An) { //за  скошенной границей
				continue;
			} else if (j >= An && i == M - j + An) { // скошенная граница
				A[k][rep[id]] = - 1;
				// x[k] = 0;
				A[k][rep[id - 1 - M]] =  1 + sqrt((dx * dx) + (dy * dy));
				b[id] = 10;
			}
			else if (j == 0 ) {//нижняя граница
				A[k][rep[id]] = 1;
				b[k] = 0;
			} else if (j == N - 1) {//Верхняя
				A[k][rep[id]]  = 1;
				b[k] = 50;
			} else if (i == 0) {//левая
				A[k][rep[id]]  = 1;
				b[k] = 100;
			} else if (i == M - 1) {//Правая
				A[k][rep[id]]  = 1;
				b[k] = 50;
			} else {
				wrtk(id, k);
			}
			k++;
		}
	}
}


void gauss(int n) {
	int i, j, k, h;
	double a;
	int p[n];
	double gl;
	int pos;
	double *tmp;
	double tm;
	//Прямой ход метода Гаусса
	for ( k = 0; k < n; k++ )  {
		gl = abs(A[k][k]);
		pos = k;
		//выбор главного элемента
		for (h = k + 1; h < n; h++)
			if (abs(A[h][k]) > gl) {
				gl = abs(A[h][k]);
				pos = h;
			}
		p[k] = pos;
		tmp = A[pos];
		A[pos] = A[k];
		A[k] = tmp;

		tm = b[pos];
		b[pos] = b[k];
		b[k] = tm;
		for ( j = k + 1; j < n; j++) {

			a = - A[j][k] / A[k][k];
			b[j] += b[k] * a;
			for ( i = k ; i < n; i++ )
				A[j][i] += A[k][i] * a;
		}
	}
	//Обратный ход
	b[n - 1] /= A[n - 1][n - 1];
	for ( i = n - 2; i >= 0; i-- ) {
		for ( j = i + 1; j < n; j++ )
			b[i] -= b[j] * A[i] [j];
		b[i] /= A[i][i];
	}
	//запись ответа и перестановка
	for (i = 0; i < n; i++)
		x[p[i]] = b[i];
	x = b;
}

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

	M = 8/dx;
	N = 6/dy;
	An = 3/dx;

	rep = (int*) malloc(sizeof(double) * N * M);

	nur = 0;
	for (j = 0; j < N; j++)
		for (i = 0; i < M; i++)
			if (j >= An && i > M - j + An)  // скошенная граница
				rep[j * M + i] = -1;
			else {
				rep[j * M + i] = nur;
				nur++;
			}

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

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

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

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

	mysolver();
	gauss(nur);
	printf("N %d  M %d  A %d\n",N, M, An );
	FILE *fd = fopen("out.txt", "w");
	//Печать решения
	for (j = 0; j < N; j++) {
		for (i = 0; i < M; i++) {
			if (j >= An && i > M - j + An)  // скошенная граница
				continue;
		fprintf(fd, "%d\t%d\t%f\n", i, j, x[rep[j * M + i]]);
		}
		fprintf(fd, "\n");
	}
	fclose(fd);

	FILE *fg = fopen("gra.gnu", "w");
	fprintf(fg, "#!/usr/bin/gnuplot -persistent\n");
	fprintf(fg, "set size square\n");
	fprintf(fg, "set output 'plot.eps'\n");
	fprintf(fg, "set cbrange [0:100]\n");
	fprintf(fg, "set xrange[0:%d]\nset yrange[0:%d]\n", M - 1, N - 1);
	fprintf(fg, "set palette rgbformulae 22,13,-31\n");
	fprintf(fg, "set pm3d map\n");
	fprintf(fg, "set pm3d flush begin ftriangles scansforwar interpolate 5,5\n");
	fprintf(fg, "splot 'out.txt' using 1:2:3 with pm3d title 'var'");
	fclose(fg);

	// int id = 620;
	// printf("%f \n", x[rep[id + M]] );
	// printf("%f \n", x[rep[id - M]] );
	// printf("%f \n", x[rep[id + 1]] );
	// printf("%f \n", x[rep[id - 1]] );
	// printf("4 %f \n", x[rep[id]] );
	// printf("%f\n", (x[rep[id + M]] + x[rep[id - M]] + x[rep[id + 1]] + x[rep[id + 1]] - 4 * x[rep[id]]) / dx / dx);

	return 0;
}