// MO9.cpp : Ten plik zawiera funkcję „main”. W nim rozpoczyna się i kończy wykonywanie programu.
//

#include <cmath>
#include "pch.h"
#include <iostream>
#include "Calc.h"
#include <fstream>
#define M_PI 3.1415926535897932384626433

double U(double x) { return 1.0 / 4.0 * (2 * x * cos(2 * x) + 2 * sin(2 * x) - log(2)*sin(2 * x) - 2 * log(cos(x))*sin(2 * x)); }
double p(double x) { return 1.; };
double q(double x) { return 0.; };
double r(double x) { return 4.; };
double s(double x) { return tan(x); };
double error(double x, double y) { return fabs(x - y); }

int main()
{

	double alfa = 0., fi = 0., beta = 1., psi = 1., gamma = 0., theta = 0.;
	int N = 10;

	std::vector<double> er_nor;
	std::vector<double> er_num;
	std::vector<double> h_v;
	std::vector<double> l;
	l.resize(N);
	auto d = l; auto u = l; auto b = l;

	double _a = 0;
	double _b = M_PI / 4;
	double h = 1e-14;
	double max_er;
	std::cout << h << std::endl;

	std::shared_ptr<Calc> calc = std::make_shared<Calc>(N);

	while (h < 0.15) {
		//metoda konwencjonalna
		d[0] = beta - alfa / h;
		u[0] = alfa / h;
		b[0] = -gamma;
		l[0] = p(_a) / (h*h) - q(_a) / (2.*h);

		double xi = _a;
		for (int k = 1; k < N-1 && xi < _b; k++) {
			xi += h;
			l[k] = p(xi) / (h*h) - q(xi) / (2.*h);
			d[k] = (-2.*p(xi)) / (h*h) + r(xi);
			u[k] = p(xi) / (h*h) + q(xi) / (2.*h);
			b[k] = -s(xi);

		}
		l[N - 1] = (-fi / h);
		d[N - 1] = (-fi / h + psi);
		b[N - 1] = (-theta);

		auto output = calc->solve_t(l, d, u, b);

		xi = _a;
		max_er = 0;

		for (int i = 0; i < N; i++) {
			double temp = error(output[i], U(xi));
			if (temp > max_er)
				max_er = temp;
			xi += h;

		}

		er_nor.push_back(max_er);

		max_er = 0;

		//metoda Numerowa
		xi = _a;
		d[0] = beta - alfa / h;
		u[0] = alfa / h;
		b[0] = -gamma;
		l[0] = p(_a) / (h*h) + 1. / 12.;
		for (int k = 1; k < N - 1; ++k) {
			xi += h;
			l[k] = p(xi) / (h*h) + 1. / 12.;
			d[k] = (-2.*p(xi)) / (h*h) + r(xi) * 10. / 12.;
			u[k] = p(xi) / (h*h) + 1. / 12.;
			b[k] = -(s(xi - h) + 10.*s(xi) + s(xi + h)) / 12.;
		}
		l[N - 1] = (-fi / h);
		d[N - 1] = (-fi / h + psi);
		b[N - 1] = (-theta);

		auto output1 = calc->solve_t(l, d, u, b);

		for (int i = 0; i < N; i++) {
			double temp = error(output1[i], U(xi));
			if (temp > max_er)
				max_er = temp;
			xi += h;

		}

		h_v.push_back(h);
		er_num.push_back(max_er);
		h += 5e-5;
	}

	/*
	d[0] = 30.0; u[0] = 2.0 / 3.0;
	l[0] = 3.0 / 4.0; d[1] = 20.0; u[1] = 5.0 / 6.0;
	l[1] = 7.0 / 8.0; d[2] = 10.0; u[2] = 9.0 / 10.0;
	l[2] = 11.0 / 12.0; d[3] = 10.0; u[3] = 13.0 / 14.0;
	l[3] = 15.0 / 16.0; d[4] = 20.0; u[4] = 17.0 / 18.0;
	l[4] = 19.0 / 20.0; d[5] = 30.0;
	b[0] = 94.0 / 3.0; b[1] = 173.0 / 4.0; b[2] = 581.0 / 20.0;
	b[3] = -815.0 / 28.0; b[4] = -6301.0 / 144.0; b[5] = -319.0 / 10.0;
	*/
	std::fstream file;
	file.open("data.csv", std::ios::out);
	if (file.good()) std::cout << "Plik otwary\n";
	file << "h;er_norm;er_numerov\n";
	for (int i = 0; i < h_v.size(); i++) {
		//printf("%5.5lf   %5.5lf\n", h_v[i], er[i]);
		file << std::fixed << std::setprecision(10) << (h_v[i]) << ";" << (er_nor[i]) << ";" << er_num[i] << std::endl;
	}

}



//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
#pragma once
#include <vector>
#include <iostream>
#include <iomanip>
class Calc
{
	int dim;

public:
	Calc(int _dim) :dim(_dim) {};
	~Calc();
	void solve_thomas(std::vector<std::vector<double>>, std::vector<double>);

	std::vector<double> solve_t(std::vector<double>, std::vector<double>, std::vector<double>, std::vector<double>);
};

/////////////////////////////////////////////////////////////////////////


std::vector<double> Calc::solve_t(std::vector<double> l, std::vector<double> d, std::vector<double> u, std::vector<double> b)
{
	std::vector<double> ni;
	ni.resize(dim);
	ni[0] = d[0];
	for (int i = 1; i < dim; i++) {
		ni[i] = d[i] - l[i-1] * (1.0 / ni[i-1]) * u[i-1];
		d[i] = ni[i];
	}
	for (int i = 1; i < dim; i++) {
		b[i] = b[i] - l[i-1] * (1.0 / ni[i - 1])*b[i-1];
	}

	std::vector<double> solVx; solVx.resize(dim);
	solVx[dim - 1] = 1.0 / ni[dim - 1] * b[dim - 1];
	for (int i = dim - 2; i >= 0; i--) {
		solVx[i] = (1.0 / ni[i])*(b[i] - u[i] * solVx[i + 1]);

	}
	return solVx;
}
////////////////