/*
plot "D:\spline_graphic.txt" with lines title "Сплайн", "D:\original_graphic.txt" with lines title "График функции"
plot "D:\difference_graphic.txt" with lines title "Разница"
*/
#include <iostream>
#include <math.h>
#include <cmath>
#include <Windows.h>
#include <stdlib.h>
#include <fstream>
using namespace std;
const double PI = 3.141592653589793238463;

//Для хранения сетки
typedef struct Node
{
    double x, y;
}
Node;

//Задание функции
typedef double(*functiontype)(double x);
double Myfunc(double x)
{
    //return x*x*x;

    if (x > 120)
        return 100 + (cos(x / 4) / 3);
    return 100;


}
functiontype Func = &Myfunc;

//Возведение числа в натуральную степень
double degree(double x, int n)
{
    double b = x;
    for (int i = 1; i < n; i++)
        x *= b;
    return x;
}


//Возвращает значение сплайна в точке на итервале[xi,x_i+1]
double Spline(double x, double* Array_steps, double* x_massiv, Node* Array, int CountSegments)
{
    double h2, g1, g2, x1, x2, y1, y2;
    int i = 0, k = 0;

    while ((i < CountSegments) & (k == 0)) {
        if ((x <= Array[i + 1].x) & (x >= Array[i].x)) {
            k = 1;
            i--;
        }
        i++;
    }


    h2 = Array_steps[i + 1];
    g1 = x_massiv[i];
    g2 = x_massiv[i + 1];
    x1 = Array[i].x;
    x2 = Array[i + 1].x;
    y1 = Array[i].y;
    y2 = Array[i + 1].y;

    return (
        (y1 * (x2 - x) + y2 * (x - x1)) / h2 +
        ((g1 * (degree((x2 - x), 3) - h2 * h2 * (x2 - x))) + (g2 * (degree((x - x1), 3) - h2 * h2 * (x - x1)))) / (6 * h2)
        );
}

//Приблизительное значение 2 производной cплайна
double approximate_derivative_2(double x, double* Selected_steps, double* x_massiv, Node* SelectedNode, int CountSegments)
{
    double delta = 1e-5;
    return ((Spline(x, Selected_steps, x_massiv, SelectedNode, CountSegments) - 2 * Spline(x + delta, Selected_steps, x_massiv, SelectedNode, CountSegments)
        + Spline(x + 2 * delta, Selected_steps, x_massiv, SelectedNode, CountSegments)) / (delta * delta));
}

//Приблизительное значение 1 производной
double approximate_derivative_1(double x, double* Selected_steps, double* x_massiv, Node* SelectedNode, int CountSegments)
{
    double delta = 1e-5;
    return (Spline(x, Selected_steps, x_massiv, SelectedNode, CountSegments) - Spline(x - delta, Selected_steps, x_massiv, SelectedNode, CountSegments)) / delta;
}

//Равномерная сетка
void ValueUniformTable(functiontype* f, Node* Array, double Initial, double End, int CountSegments, double* Array_steps_uni)
{
    double h;
    ofstream fout4("D:\setka_uni.txt");
    Array[0].x = Initial;
    Array[0].y = (*f)(Initial);
    fout4 << Array[0].x << " " << Array[0].y << endl;
    Array_steps_uni[0] = 0;
    h = abs(End - Initial) / CountSegments;

    for (int i = 1; i <= CountSegments; i++)
    {
        Array[i].x = Array[i - 1].x + h;
        Array[i].y = (*f)(Array[i].x);
        fout4 << Array[i].x << " " << Array[i].y << endl;
        Array_steps_uni[i] = h;
    }
    fout4.close();
}

//Чебышёвская сетка
void ValueChebyshevTable(functiontype* f, Node* Array, double Initial, double End, int CountSegments, double* Array_steps_cheb)
{
    ofstream fout5("D:\setka_cheb.txt");
    Array_steps_cheb[0] = 0;
    Array[0].x = Initial;
    Array[0].y = (*f)(Array[0].x);
    fout5 << Array[0].x << " " << Array[0].y << endl;
    double k1, k2, denom;
    k1 = (End + Initial) / 2;
    k2 = (End - Initial) / 2;
    denom = 2 * (CountSegments + 1);
    for (int i = 1; i < CountSegments; i++)
    {
        Array[CountSegments - i].x = k1 + k2 * cos(((2 * i + 1) * PI) / (denom));
        Array[CountSegments - i].y = (*f)(Array[CountSegments - i].x);
        fout5 << Array[CountSegments - i].x << " " << Array[CountSegments - i].y << endl;
    }

    Array[CountSegments].x = End;
    Array[CountSegments].y = (*f)(Array[CountSegments].x);
    Array_steps_cheb[CountSegments] = Array[CountSegments].x - Array[CountSegments - 1].x;
    fout5 << Array[CountSegments].x << " " << Array[CountSegments].y << endl;
    fout5.close();

    for (int i = 1; i < CountSegments; i++)
        Array_steps_cheb[i] = Array[i].x - Array[i - 1].x;
}

//Составляет СЛАУ с учётом нач. данных и заполняет этим массив matrix_coeffs - массив вида[[a1,b1,c1,d1], [a2,b2,c2,d2],..., [an,bn,cn,dn]]
void get_matrix_coeffs(functiontype* f, double** matrix_coeffs, double* Array_steps, Node* Array, double Initial, double End, int CountSegments)
{
    int i;
    double h1, h2, A, B, y1, y2;
    int matrix_size = CountSegments + 1;

    //A = approximate_derivative_2(&Func, Array[0].x);
    //B = approximate_derivative_1(&Func, Array[CountSegments].x);
    A = 12;
    B = -15;

    //первое уравнение
    matrix_coeffs[0][0] = 0;
    matrix_coeffs[0][1] = 2;
    matrix_coeffs[0][2] = 1;
    matrix_coeffs[0][3] = 6 * (y2 - y1 - A * h2) / (h2 * h2);

    for (i = 1; i < matrix_size - 1; i++)
    {
        h1 = Array_steps[i];
        h2 = Array_steps[i + 1];
        matrix_coeffs[i][0] = h1;
        matrix_coeffs[i][1] = 2 * (h1 + h2);
        matrix_coeffs[i][2] = h2;
        matrix_coeffs[i][3] = 6 * ((Array[i + 1].y - Array[i].y) / h2 - (Array[i].y - Array[i - 1].y) / h1);
    }

    //n+1 уравнение
    h2 = Array_steps[matrix_size - 1];
    y1 = Array[matrix_size - 2].y;
    y2 = Array[matrix_size - 1].y;
    matrix_coeffs[matrix_size - 1][0] = 0;
    matrix_coeffs[matrix_size - 1][1] = 1;
    matrix_coeffs[matrix_size - 1][2] = 0;
    matrix_coeffs[matrix_size - 1][3] = B;
}

/*Принимает на вход коэффициенты СЛАУ в виде массива matrix_coeffs - массив вида[[a1,b1,c1,d1], [a2,b2,c2,d2],..., [an,bn,cn,dn]].
   Заполняет Двумерный массив прогоночных коэффициентов coeffs_massiv (0 - Кси, 1 - Эта). Находит решение в виде массива x_massiv*/
void tridiagonal_matrix_algorithm(double** matrix_coeffs, double* Array_steps, Node* Array, double* x_massiv, int CountSegments)
{
    int i;
    double denominator;
    double K, E, K_prev, E_prev;
    int matrix_size = CountSegments + 1;

    double** coeffs_massiv;
    coeffs_massiv = new double* [matrix_size];
    for (i = 0; i < matrix_size; i++)
        coeffs_massiv[i] = new double[2];

    //Прямой ход
    coeffs_massiv[0][0] = -(matrix_coeffs[0][2]) / matrix_coeffs[0][1];
    coeffs_massiv[0][1] = (matrix_coeffs[0][3]) / matrix_coeffs[0][1];
    K_prev = coeffs_massiv[0][0];
    E_prev = coeffs_massiv[0][1];

    for (i = 1; i <= matrix_size - 1; i++)
    {
        denominator = matrix_coeffs[i][0] * K_prev + matrix_coeffs[i][1];
        K = -(matrix_coeffs[i][2]) / denominator;
        E = (matrix_coeffs[i][3] - matrix_coeffs[i][0] * E_prev) / denominator;
        K_prev = K;
        E_prev = E;
        coeffs_massiv[i][0] = K;
        coeffs_massiv[i][1] = E;
    }

    //Обратный ход
    x_massiv[matrix_size - 1] = coeffs_massiv[matrix_size - 1][1];
    for (i = matrix_size - 2; i >= 0; i--)
    {
        x_massiv[i] = coeffs_massiv[i][0] * x_massiv[i + 1] + coeffs_massiv[i][1];
    }
}

//Запись в 3 файла: значений сплайна, интерполируемой функции и их разницы в countdots точках
void tables_in_file(functiontype* f, Node* Array, int Count_dots, int CountSegments, double Initial, double End, double* Array_steps, double* x_massiv)
{
    int i = 0;
    double step = (End - Initial) / (Count_dots - 1), y_value1, y_value2, x_value;

    ofstream fout1("D:/spline_graphic.txt");
    ofstream fout2("D:/original_graphic.txt");
    ofstream fout3("D:/difference_graphic.txt");
    x_value = Initial;
    for (i = 0; i < CountSegments; i++)
    {
        while (x_value < Array[i + 1].x)
        {
            fout1 << x_value << " ";
            fout2 << x_value << " ";
            fout3 << x_value << " ";
            y_value1 = Spline(x_value, Array_steps, x_massiv, Array, CountSegments);
            y_value2 = (*f)(x_value);
            fout1 << y_value1 << endl;
            fout2 << y_value2 << endl;
            fout3 << abs(y_value2 - y_value1) << endl;
            x_value += step;
        }
    }
    fout1 << End << " ";
    fout2 << End << " ";
    fout1 << Spline(End, Array_steps, x_massiv, Array, CountSegments) << endl;
    fout2 << (*f)(End) << endl;
    fout1.close();
    fout2.close();
    //cout << endl << "Запись в файл осуществлена!" << endl;
}

int main()
{
    setlocale(LC_ALL, "RUS");
    functiontype Func = &Myfunc;
    int CountSegments = 5, i, Countdots = 250;

    /*---------------------------------Входные данные-------------------------------------*/
    double Initial = 100, End = 130;


    /*---------------------------------Сетки-------------------------------------*/
    double* Array_steps_cheb = new double[CountSegments + 1];
    double* Array_steps_uni = new double[CountSegments + 1];
    double* Selected_steps = new double[CountSegments + 1];

    //Равномерная сетка
    Node* ArrayUniformNodes = new Node[CountSegments + 1];
    ValueUniformTable(&Func, ArrayUniformNodes, Initial, End, CountSegments, Array_steps_uni);
    //Чебышёв
    Node* ArrayChebyshevNodes = new Node[CountSegments + 1];
    ValueChebyshevTable(&Func, ArrayChebyshevNodes, Initial, End, CountSegments, Array_steps_cheb);

    Node* SelectedNode = new Node[CountSegments + 1];

    //Array_steps_cheb  ArrayChebyshevNodes  ArrayUniformNodes  Array_steps_uni

    SelectedNode = ArrayChebyshevNodes;
    Selected_steps = Array_steps_cheb;

    /*---------------------------------СЛАУ - matrix_coeffs -------------------------------------*/
    double** matrix_coeffs;
    matrix_coeffs = new double* [CountSegments + 1];
    for (i = 0; i < CountSegments + 1; i++)
        matrix_coeffs[i] = new double[3];

    get_matrix_coeffs(&Func, matrix_coeffs, Selected_steps, SelectedNode, Initial, End, CountSegments);

    /*---------------------------------Прогонка - запись в x_massiv -------------------------------------*/
    double* x_massiv = new double[CountSegments + 1];
    tridiagonal_matrix_algorithm(matrix_coeffs, Selected_steps, SelectedNode, x_massiv, CountSegments);

    /*---------------------------------Графики -------------------------------------*/
    tables_in_file(&Func, SelectedNode, Countdots, CountSegments, Initial, End, Selected_steps, x_massiv);

    /*---------------------------------Вывод в консоль значений -------------------------------------*/

    cout << "Проверка А: " << approximate_derivative_2(Initial, Selected_steps, x_massiv, SelectedNode, CountSegments) << endl;
    cout << "Проверка B: " << approximate_derivative_1(End, Selected_steps, x_massiv, SelectedNode, CountSegments) << endl;
    //Шаги
    cout << endl;
    for (i = 1; i <= CountSegments; i++) {
        cout << "h_" << i << ": " << Selected_steps[i] << endl;
    }

    cout << endl << "Cетка:  " << endl;
    for (int i = 0; i < CountSegments + 1; i++)

        cout << "(" << SelectedNode[i].x << ":" << SelectedNode[i].y << ")" << endl;
    cout << endl;


    //Гаммы
    for (i = 0; i < CountSegments + 1; i++)
        cout << "g_" << i << ": " << x_massiv[i] << endl;
    cout << endl;

    //Значение сплайна в узлах
    for (i = 0; i < CountSegments; i++)
        cout << "Spline in x_" << i << " " << Spline(SelectedNode[i].x, Selected_steps, x_massiv, SelectedNode, CountSegments) << endl;
    cout << "Spline in x_" << CountSegments << " " << Spline(SelectedNode[CountSegments].x, Selected_steps, x_massiv, SelectedNode, CountSegments) << endl;


    cout << endl;
    system("pause");
    return 0;
}