#include <iostream>
#include <cmath>
#include <fstream>

using namespace std;

const int number_of_points_x = 10000;
const int number_of_points_t = 1000;
const int t1 = 2;
const int t2 = 5;
const int length = 1;
const double dt1 = (double) t1 / number_of_points_t;
const double dt2 = t2 / number_of_points_t;
const double dx = (double) length / number_of_points_x;
const int point_of_start_internal_condition = number_of_points_x / 3;
const int border_of_internal_condition = number_of_points_x / 10;
double upper_diagonal[number_of_points_x];
double lower_diagonal[number_of_points_x];
double main_diagonal[number_of_points_x];
double external_function[number_of_points_x];
double intial_conditions[number_of_points_x];
double u_x_t[number_of_points_x];

void create_approximation() {
    for (int i = 0; i < number_of_points_x; i++) {
        ::upper_diagonal[i] = 1;
        ::lower_diagonal[i] = 1;
        ::main_diagonal[i] = -2;
        ::external_function[i] = 0;
        ::intial_conditions[i] = 0;
    }
    cout << "end1" << endl;
    ::upper_diagonal[number_of_points_x - 1] = 0;
    ::lower_diagonal[0] = 0;
    ::main_diagonal[0] = -1;
    ::main_diagonal[number_of_points_x - 1] = -2 + (1 - 1 / (2 * number_of_points_x)) / (1 + 1 / (2 * number_of_points_x));
    ::external_function[0] = -1 / number_of_points_x;
    cout << "end2" << endl;
    for (int i = point_of_start_internal_condition - border_of_internal_condition; i < point_of_start_internal_condition + border_of_internal_condition + 1; i++) {
        intial_conditions[i] = 1;
    }
    cout << "end3" << endl;
}

void return_u_x_t_to_intital_condition() {
    for (int i = 0; i < number_of_points_x; i++) {
        u_x_t[i] = intial_conditions[i];
    }
}
void multiplication_3dig_matrix_on_vector(int dim, double* lower_diagonal, double* main_diagonal, double* upper_diagonal, double* vector, double* result) {
    result[0] = main_diagonal[0] * vector[0] + upper_diagonal[0] * vector[1];
    for (int i = 1; i < dim - 1; i++) {
        result[i] = lower_diagonal[i] * vector[i - 1] + main_diagonal[i] * vector[i] + upper_diagonal[i] * vector[i + 1];
    }
    result[dim - 1] = lower_diagonal[dim - 1] * vector[dim - 2] + main_diagonal[dim - 1] * vector[dim - 1];
}

void calculate_for_2T() {
    return_u_x_t_to_intital_condition();
    double result[number_of_points_x];
    for (double i = 0; i < t1 + dt1; i += dt1) {
        multiplication_3dig_matrix_on_vector(number_of_points_x, lower_diagonal, main_diagonal, upper_diagonal, u_x_t, result);
        for (int j = 0; j < number_of_points_x; j++) {
            u_x_t[j] = u_x_t[j] + dt1 * (result[j] * ((double)number_of_points_x *(double)number_of_points_x) / (length * length) + external_function[j]);
            cout<<"i = " <<i << " j = " << j << "u = "<< u_x_t[j] << endl;
        }
    }
    fstream f;
    f.open("C://Users//User//Documents//EULER//euler_2T.txt", fstream::trunc | fstream::out);
    for (int i = 0; i < number_of_points_x; i++) {
        f << u_x_t[i] << endl;
    }
    f.close();
    return_u_x_t_to_intital_condition();
    cout << "end" << endl;
}

void calculate_for_5T() {
    return_u_x_t_to_intital_condition();
    double result[number_of_points_x];
    for (double i = 0; i < t2 + dt2; i += dt2) {
        multiplication_3dig_matrix_on_vector(number_of_points_x, lower_diagonal, main_diagonal, upper_diagonal, u_x_t, result);
        for (int i = 0; i < number_of_points_x; i++) {
            u_x_t[i] = u_x_t[i] + dt2 * (result[i] + external_function[i]);
        }
    }
    fstream f;
    f.open("C://Users//User//Documents//EULER//euler_5T.txt", fstream::trunc | fstream::out);
    for (int i = 0; i < number_of_points_x; i++) {
        f << u_x_t[i] << endl;
    }
    f.close();
    return_u_x_t_to_intital_condition();
}

void calculate_for_033L() {
    return_u_x_t_to_intital_condition();
    double result[number_of_points_x];
    fstream f;
    f.open("C://Users//User//Documents//EULER//euler_033L.txt", fstream::trunc | fstream::out);
    for (double i = 0; i < t2 + dt2; i += dt2) {
        multiplication_3dig_matrix_on_vector(number_of_points_x, lower_diagonal, main_diagonal, upper_diagonal, u_x_t, result);
        for (int i = 0; i < number_of_points_x; i++) {
            u_x_t[i] = u_x_t[i] + dt2 * (result[i] + external_function[i]);
        }
        f << u_x_t[point_of_start_internal_condition] << endl;
    }
    f.close();
    return_u_x_t_to_intital_condition();
}

void output_intial_conditions() {
    fstream f;
    f.open("C://Users//User//Documents//EULER//euler_0.txt", fstream::trunc | fstream::out);
    for (int i = 0; i < number_of_points_x; i++) {
        f << intial_conditions[i] << endl;
    }
    f.close();
    cout << "end5" << endl;
}
void output_x() {
    fstream f;
    f.open("C://Users//User//Documents//EULER//euler_x.txt", fstream::trunc | fstream::out);
    for (double i = 0; i < length + dx; i += dx) {
        f << i << endl;
    }
    f.close();
    cout << "end4" << endl;
}

void output_t1() {
    fstream f;
    f.open("C://Users//User//Documents//EULER//euler_t1.txt", fstream::trunc | fstream::out);
    for (double i = 0; i < t1 + dt1; i += dt1) {
        f << i << endl;
    }
    f.close();
}

void output_t2() {
    fstream f;
    f.open("C://Users//User//Documents//EULER// euler_t2.txt", fstream::trunc | fstream::out);
    for (double i = 0; i < t2 + dt2; i += dt2) {
        f << i << endl;
    }
}
int main() {
    create_approximation();
    calculate_for_2T();
    //calculate_for_5T();
    //calculate_for_033L();
    output_t1();
    //output_t2();
    output_x();
    output_intial_conditions();
}