#define _USE_MATH_DEFINES

#include <iostream>
#include <cmath>
#include <fstream>
#include "calerf.h"

const double D = 1.0;
const double tMax = 2.0;
const double a = 6.0 * sqrt(D * tMax);
const double LAMBDA = 1.0;
const double h = 0.1;
const double dt = (LAMBDA * (h*h) / D);
const int N = (int)((tMax / dt) + 2);
const int M = (int)((a / h) + 1);

double *newWektor(int n);
double **newMacierz(int n, int m);
void deleteWektor(double* x);
void deleteMacierz(double** x, int n);

void Wypelnij(double* x, double delta, int k);

void ZapiszPlik(const char* fileName, double** macierz, int r, int c);
void ZapiszPLik2(char *filename, double **macierz, int r, int c, double *hPoziom);

void RozwiazanieAnalityczne(double** analit, int N, int M, double* dtPoziom, double* hPoziom, double dt, double h);
void CrankNicolsonGS(double** U, int N, int M, double* dtPoziom, double* hPoziom, double dt, double h, double LAMBDA);
void CrankNicolsonLU(double** U, int N, int M, double* dtPoziom, double* hPoziom, double dt, double h, double LAMBDA);
void GS(double* L, double* D, double* U, double* b, double* x, int n);
void LU(double* L, double* D, double* U, double* b, double* x, int n);

void BladGS(double **err, double **Analityczne, double **U, int N, int M, double *dt_kroki, double *h_kroki);
void BladLU(double **err, double **Analityczne, double **U, int N, int M, double *dt_kroki, double *h_kroki);

int main()
{
	double* dtPoziom = newWektor(N);
	double* hPoziom = newWektor(M);

	std::cout << "Zapisanie do pliku." << std::endl;

	Wypelnij(dtPoziom, dt, N);
	Wypelnij(hPoziom, h, M);

	std::ofstream file("dt.txt");
	file.setf(std::ios::fixed, std::ios::floatfield);
	file.precision(4);
	for (int i = 0; i < N; i++)
	{
		file << dtPoziom[i] << std::endl;
	}

	file.close();

	std::ofstream file2("h.txt");
	file2.setf(std::ios::fixed, std::ios::floatfield);
	file2.precision(4);
	for (int i = 0; i < M; i++)
	{
		file2 << hPoziom[i] << std::endl;
	}

	file2.close();

	double** A = newMacierz(N, M);
	RozwiazanieAnalityczne(A, N, M, dtPoziom, hPoziom, dt, h);

	double** Err = newMacierz(N, M);

	double** U = newMacierz(N, M);
	CrankNicolsonGS(U, N, M, dtPoziom, hPoziom, dt, h, LAMBDA);
	ZapiszPLik2("Crank+GS.txt", U, N, M, hPoziom);
	BladGS(Err, A, U, N, M, dtPoziom, hPoziom);
	deleteMacierz(U, N);

	U = newMacierz(N, M);
	CrankNicolsonLU(U, N, M, dtPoziom, hPoziom, dt, h, LAMBDA);
	ZapiszPLik2("Crank+LU.txt", U, N, M, hPoziom);
	BladLU(Err, A, U, N, M, dtPoziom, hPoziom);

	system("pause");
	return 0;
}

double *newWektor(int n)
{
	return new double[n];
}

double **newMacierz(int n, int m)
{
	double** macierz = new double*[n];

	for (int i = 0; i < n; i++)
	{
		macierz[i] = new double[m];
	}

	return macierz;
}

void deleteWektor(double* x)
{
	delete[] x;
}

void deleteMacierz(double** x, int n)
{
	for (int i = 0; i < n; i++)
	{
		delete[] x[i];
	}
	delete[] x;
}

void Wypelnij(double* x, double delta, int k)
{
	x[0] = 0.0;
	for (int i = 1; i<k; i++)
		x[i] = x[i - 1] + delta;
}

void RozwiazanieAnalityczne(double** analit, int N, int M, double* dtPoziom, double* hPoziom, double dt, double h)
{
	for (int i = 0; i < N; i++)
	{
		analit[i][0] = 0.0;
	}
	for (int i = 0; i < N; i++)
	{
		analit[i][M - 1] = 1.0;
	}

	for (int i = 0; i < M; i++)
	{
		analit[0][i] = 1.0;
	}

	for (int i = 1; i < N; i++)
	{
		for (int j = 1; j < M - 1; j++)
		{
			analit[i][j] = erf(hPoziom[j] / (2 * sqrt(D * dtPoziom[i])));
		}
	}

	ZapiszPLik2("RozwiazanieAnalityczne.txt", analit, N, M, hPoziom);
}

void CrankNicolsonGS(double** U, int N, int M, double* dtPoziom, double* hPoziom, double dt, double h, double LAMBDA)
{
	double* L = newWektor(M);
	double* Up = newWektor(M);
	double* D = newWektor(M);
	double* b = newWektor(M);
	double* u = newWektor(M);

	for (int i = 0; i < N; i++)
	{
		U[i][0] = 0.0; //Warunek brzegowy U(0, t) = 0.0.
	}
	for (int i = 0; i < N; i++)
	{
		U[i][M - 1] = 1.0; //Warunek brzegowy U(INF, t) = 1.0.
	}

	for (int i = 0; i < M; i++)
	{
		U[0][i] = 1.0; //Warunek poczatkowy U(X, 0) = 1.0.
	}

	for (int k = 1; k < N; k++)
	{
		L[0] = LAMBDA / 2.0;
		D[0] = 1.0;
		Up[0] = 0.0;
		b[0] = 0.0;

		for (int i = 1; i < M - 1; i++)
		{
			L[i] = LAMBDA / 2.0;
			D[i] = -(1 + LAMBDA);
			Up[i] = LAMBDA / 2.0;
			b[i] = -(LAMBDA / 2.0 * U[k - 1][i - 1] + (1.0 - LAMBDA) * U[k - 1][i] + (LAMBDA / 2.0) * U[k - 1][i + 1]);
		}

		L[M - 1] = 0.0;
		D[M - 1] = 1.0;
		Up[M - 1] = 0.0;
		b[M - 1] = 1.0;

		GS(L, D, Up, b, u, M);

		for (int i = 1; i < M - 1; i++)
		{
			U[k][i] = u[i];
		}
	}
}

void CrankNicolsonLU(double** U, int N, int M, double* dtPoziom, double* hPoziom, double dt, double h, double LAMBDA)
{
	double* L = newWektor(M);
	double* Up = newWektor(M);
	double* D = newWektor(M);
	double* b = newWektor(M);
	double* u = newWektor(M);

	for (int i = 0; i < N; i++)
	{
		U[i][0] = 0.0;	//Warunek brzegowy U(0, t) = 0.0.
	}
	for (int i = 0; i < N; i++)
	{
		U[i][M - 1] = 1.0; //Warunek brzegowy U(INF, t) = 1.0.
	}

	for (int i = 0; i < M; i++)
	{
		U[0][i] = 1.0; //Warunek początkowy U(X, 0) = 1.0.
	}

	for (int k = 1; k < N; k++)
	{
		L[0] = LAMBDA / 2.0;
		D[0] = 1.0;
		Up[0] = 0.0;
		b[0] = 0.0;

		for (int i = 1; i < M - 1; i++)
		{
			L[i] = LAMBDA / 2.0;
			D[i] = -(1 + LAMBDA);
			Up[i] = LAMBDA / 2.0;
			b[i] = -(LAMBDA / 2.0 * U[k - 1][i - 1] + (1.0 - LAMBDA) * U[k - 1][i] + (LAMBDA / 2.0) * U[k - 1][i + 1]);
		}

		L[M - 1] = 0.0;
		D[M - 1] = 1.0;
		Up[M - 2] = 0.0;
		b[M - 1] = U[k-1][M-1];

		LU(L, D, Up, b, u, M);

		for (int i = 1; i < M - 1; i++)
		{
			U[k][i] = u[i];
		}
	}
}

void GS(double* L, double* D, double* U, double* b, double* x, int n)
{
	double** A = newMacierz(n, n);

	for (int i = 0; i < n; i++)
	{
		for (int j = 0; j < n; j++)
		{
			A[i][j] = 0.0;
		}
	}
	for (int i = 0; i < n - 1; i++)
	{
		A[i][i] = D[i];
		A[i + 1][i] = U[i];
		A[i][i + 1] = L[i];
	}
	A[n - 1][n - 1] = D[n - 1];

	double x1;
	for (int k = 0; k<M - 1; k++)
	{
		for (int i = k + 1; i<M; i++)
		{
			x1 = A[i][k] / A[k][k];
			A[i][k] = x1;
			for (int j = k + 1; j<M; j++)
			{
				A[i][j] = A[i][j] - (x1*A[k][j]);
			}
		}
	}
	double suma;
	double *z = new double[M];

	for (int i = 0; i<M; i++)
	{
		suma = 0.0;
		for (int j = 0; j <= i - 1; j++)
		{
			suma += A[i][j] * z[j];
		}

		z[i] = b[i] - suma;
	}

	for (int i = M - 1; i >= 0; i--)
	{
		suma = 0.0;
		for (int j = i + 1; j<M; j++)
		{
			suma += A[i][j] * x[j];
		}

		x[i] = (z[i] - suma) / A[i][i];
	}
}

void LU(double* L, double* D, double* U, double* b, double* x, int n)
{
	double** A = newMacierz(n, n);

	for (int i = 0; i < n; i++)
	{
		for (int j = 0; j < n; j++)
		{
			A[i][j] = 0.0;
		}
	}
	for (int i = 0; i < n - 1; i++)
	{
		A[i][i] = D[i];
		A[i + 1][i] = U[i];
		A[i][i + 1] = L[i];
	}
	A[n - 1][n - 1] = D[n - 1];

//Dekompozycja LU, metoda eliminacji Gaussa
	double x1;
	for (int k = 0; k<M - 1; k++)
	{
		for (int i = k + 1; i<M; i++)
		{
			x1 = A[i][k] / A[k][k];
			A[i][k] = x1;
			for (int j = k + 1; j<M; j++)
			{
				A[i][j] = A[i][j] - (x1*A[k][j]);
			}
		}
	}

//Rozwiązywanie układu równań
	double suma;
	double *z = new double[M];

//Ly=b => y=...
	for (int i = 0; i<M; i++)
	{
		suma = 0.0;
		for (int j = 0; j <= i - 1; j++)
		{
			suma += A[i][j] * z[j];
		}

		z[i] = b[i] - suma;
	}

//Ux=y => x=...
	for (int i = M - 1; i >= 0; i--)
	{
		suma = 0.0;
		for (int j = i + 1; j<M; j++)
		{
			suma += A[i][j] * x[j];
		}

		x[i] = (z[i] - suma) / A[i][i];
	}
}

void BladGS(double **err, double **Analityczne, double **U, int N, int M, double *dt_kroki, double *h_kroki)
{
	std::ofstream file1;
	std::ofstream file2;

	file1.open("Blad_max_co_krok_t_GS.txt");
	file2.open("Blad_max_od_kroku_h_GS.txt", std::ios::app);

	file1.setf(std::ios::fixed, std::ios::floatfield);
	file2.setf(std::ios::fixed, std::ios::floatfield);

	file1.precision(15);
	file2.precision(15);

	double max = 0.0;

	for (int i = 0; i < N; i++)
	{
		max = 0.0;
		for (int j = 0; j < M; j++)
		{
			err[i][j] = fabs(U[i][j] - Analityczne[i][j]);

			if (err[i][j] > max)
				max = err[i][j];
		}
		if (i == (N - 1))
		{
			file2 << log10(h) << "\t" << log10(max) << std::endl;
		}
		if (dt_kroki[i] != 0.0)
		{
			file1 << dt_kroki[i] << "\t";
			file1 << max << std::endl;
		}

	}

	ZapiszPlik("BladGS.txt", err, N, M);
}

//Obliczanie bledu dla CN+LU.
void BladLU(double **err, double **Analityczne, double **U, int N, int M, double *dt_kroki, double *h_kroki)
{
	std::ofstream file1;
	std::ofstream file2;

	file1.open("Blad_max_co_krok_t_LU.txt");
	file2.open("Blad_max_od_kroku_h_LU.txt", std::ios::app);

	file1.setf(std::ios::fixed, std::ios::floatfield);
	file2.setf(std::ios::fixed, std::ios::floatfield);

	file1.precision(15);
	file2.precision(15);

	double max = 0.0;

	for (int i = 0; i < N; i++)
	{
		max = 0.0;
		for (int j = 0; j < M; j++)
		{
			err[i][j] = fabs(U[i][j] - Analityczne[i][j]);

			if (err[i][j] > max)
				max = err[i][j];
		}
		if (i == (N - 1))
		{
			file2 << (log10(h)) << "\t" << (log10(max)) << std::endl;
		}


		if (dt_kroki[i] != 0.0)
		{
			file1 << dt_kroki[i] << "\t";
			file1 << max << std::endl;
		}
	}

	ZapiszPlik("BladLU.txt", err, N, M);
}

void ZapiszPlik(const char* fileName, double** macierz, int r, int c)
{
	std::ofstream file(fileName);

	file.setf(std::ios::scientific, std::ios::floatfield);
	file.precision(4);

	for (int i = 0; i < r; i++)
	{
		for (int j = 0; j < c; j++)
		{
			file << macierz[i][j] << "\t";
		}
		file << std::endl;
	}

	file.close();
}

void ZapiszPLik2(char *filename, double **macierz, int r, int c, double *hPoziom)
{
	std::ofstream file(filename);

	file.setf(std::ios::scientific, std::ios::floatfield);
	file.precision(4);
	file << std::endl;
	for (int i = 0; i < c; i++)
	{
		file << std::endl;
		file << hPoziom[i] << "\t";
		for (int j = 0; j < c; j++)
		{
			file << macierz[j][i] << "\t";
		}
	}
	file.close();
}