#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>
#define xPoczatkowe 0.1
using namespace std;
class macierz {
	public:
		double* D;
		double* U;
		double* L;
		int size;
		macierz(double * przekatna, double *przekatnaU,double *przekatnaL,int size) {
			this->size = size;
			D = new double[size];
			copy(przekatna,przekatna+size,D);
			U = new double[size-1];
			copy(przekatnaU,przekatnaU+size-1,U);
			L = new double[size-1];
			copy(przekatnaL,przekatnaL+size-1,L);
		}
		~macierz() {
			delete [] D;
			delete [] L;
			delete [] U;
		}
};

double wzorAnalityczny(double x) {
	return sin(x)-cos(x)+log(1.0/tan(x/2.0))*sin(x);
}

double ctg(double x) {
	return 1.0 / tan(x);
}
void rozwiaz(macierz *m, double *b, double *x, int N) {
	for (int i = N - 1; i >= 0; --i) {
		x[i] = (b[i] - m->L[i] * x[i + 1]) / m->D[i];
	}
}
void algorytmThomasa(macierz *m, double *b, double *x, int N) {
	for (int i = 1; i < N; ++i)
		m->D[i] = m->D[i] - (m->U[i] * m->L[i - 1]) / m->D[i - 1];
	for (int i = 1; i < N; ++i)
		b[i] = b[i] - (m->U[i] * b[i - 1]) / m->D[i - 1];
	rozwiaz(m,b,x,N);
}

double max(double *w, int size) {
	double max = fabs(w[0]);
	for (int i = 0; i < size; i++) {
		if (fabs(w[i]) > max) {
			max = fabs(w[i]);
		}
	}
	return max;
}

double dyskretyzacjaNumerowa(double h, int N) {
	//inicjalizacja
	double xP = xPoczatkowe;
	double xP_kopia = xPoczatkowe;
	double *L = new double[N];
	double *D = new double[N];
	double *U = new double[N];
	double *b = new double[N];
	double *x = new double[N];
	double *BLEDY = new double[N];
	//warunki brzegowe
	U[0] = 0.0 / h;
	D[0] = 1.0 / h;
	b[0] = -0.5961798034812780886272;
	L[0] = 0;
	//uzupełnianie macierzy bez brzegow
	for (int i=1; i<N-1; i++) {
		L[i] = 1.0/(h*h)+1.0/12.0;
		D[i] = (-2.0*1.0)/(h*h)+1.0*(10.0/12.0);
		U[i] = 1.0/(h*h)+1.0/12.0;
		b[i] = -ctg(xP+i*h-h)/12.0-(10.0/12.0)*ctg(xP+i*h)-ctg(xP+i*h+h)/12.0;
	}
	//warunki brzegowe
	L[N-1] = 0.0/h;
	D[N-1] = 0.0/h+1.0;
	b[N-1] = 1.0;
	U[N-1] = 0;
	//rozwiazywanie
	macierz *m = new macierz(D,U,L,N);
	algorytmThomasa(m, b, x, N);
	//bledy
	for(int i=0; i<N; i++,xP+=h) {
		BLEDY[i]=fabs(x[i]-wzorAnalityczny(xP));
	}
	if (N == 60) {
		fstream wyjscie;
		wyjscie.open("numerow.txt", fstream::out);
		cout << endl << "            Dyskretyzacja   Numerowa " << endl;
		cout << "    i         punkt            x[n]           U(x)  \n" << endl;
		for (int i = 0; i < N; i++) {
			wyjscie << xP_kopia << " " << x[i] << " " << wzorAnalityczny(xP_kopia) << " " << endl;
			if (i < 25)
				cout << setw(5) << i << setw(15) << xP_kopia << setw(15) << x[i] << setw(15) << wzorAnalityczny(xP_kopia) << endl;
			xP_kopia += h;
		}
		wyjscie.close();
	}

	delete[] L;
	delete[] D;
	delete[] U;
	delete[] x;
	delete[] b;
	double maxError = max(BLEDY, N);
	delete[] BLEDY;
	return maxError;
}

double dyskretyzacjaKonwencjonalna(double h, int N) {
	double *L = new double[N];
	double *D = new double[N];
	double *U = new double[N];
	double *b = new double[N];
	double *x = new double[N];
	double *BLEDY = new double[N];
	double xP = xPoczatkowe;
	double xP_kopia = xPoczatkowe;
	U[0] = 0.0/h;
	D[0] = 1.0/h;
	b[0] = -0.5961798034812780886272;
	L[0] = 0;
	for (int i=1; i<N-1; i++) {
		L[i] = 1.0/(h*h)-0.0/(2.0*h);
		D[i] = (-2.0*1.0)/(h*h)+1.0;
		U[i] = 1.0/(h*h)-0.0/(2.0*h);
		b[i] = -ctg(xP+i*h);
	}
	L[N-1] = 0.0/h;
	D[N-1] = 0.0/h+1.0;
	b[N-1] = 1.0;
	U[N-1] = 0;
	macierz *m = new macierz(D,U,L,N);
	algorytmThomasa(m, b, x, N);
	for (int i=0; i<N; i++,xP += h) {
		BLEDY[i] = fabs(x[i]-wzorAnalityczny(xP));
	}
	if (N == 60) {
		fstream wyjscie;
		wyjscie.open("konwencjonal.txt", fstream::out);
		cout << endl << "    Dyskretyzacja Konwencjonalna   " << endl;
		cout << "    i         punkt            x[n]           U(x)  \n" << endl;
		for (int i = 0; i < N; i++) {
			BLEDY[i] = fabs(x[i] - wzorAnalityczny(xP_kopia));
			wyjscie << xP_kopia << " " << x[i] << " " << wzorAnalityczny(xP_kopia) << " " << endl;
			if (i < 25)
				cout << setw(5) << i << setw(15) << xP_kopia << setw(15) << x[i] << setw(15) << wzorAnalityczny(xP_kopia) << endl;
			xP_kopia += h;
		}
		wyjscie.close();
	}
	delete[] L;
	delete[] D;
	delete[] U;
	delete[] x;
	delete[] b;
	double maxError = max(BLEDY, N);
	delete[] BLEDY;
	return maxError;
}

int main() {
	fstream wyjscie;
	wyjscie.open("bledy.txt", fstream::out);

	for (int N = 10; N < 70000; N += 50) {
		double krok = (M_PI_2-xPoczatkowe)/(N-1);
		double konwencjonalna = dyskretyzacjaKonwencjonalna(krok, N);
		double numerowa = dyskretyzacjaNumerowa(krok, N);
		wyjscie << log10(krok) << " " << log10(konwencjonalna) << " " << log10(numerowa) << endl;
	}
	wyjscie.close();
	return 0;
}