#include <armadillo>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <ios>
#include <string>

arma::vec* general_algorithm(
    arma::vec* sub_diag,
    arma::vec* main_diag,
    arma::vec* sup_diag,
    arma::vec* g_vec
)
{

}


arma::vec* special_algorithm(
    double sub_sig,
    double main_sig,
    double sup_sig,
    arma::vec* g_vec
)
{

}

void error(
    std::string filename,
    arma::vec* x_vec,
    arma::vec* v_vec,
    arma::vec* a_vec
)
{
    std::ofstream ofile;
    ofile.open(filename);

    if (!ofile.is_open()) {
        exit(1);
    }

    for (int i=0; i < a_vec->n_elem; i++) {
        double sub = (*a_vec)(i) - (*v_vec)(i);
        ofile << std::setprecision(8) << std::scientific << (*x_vec)(i)
            << std::setprecision(8) << std::scientific << std::log10(std::abs(sub))
            << std::setprecision(8) << std::scientific << std::log10(std::abs(sub/(*a_vec)(i)))
            << std::endl;
    }
    
    ofile.close();

}

double analytic_solution(double x) {
    return 100*std::exp(-10*x);
}

void build_array(
    int n_steps, 
    arma::vec* sub_diag, 
    arma::vec* main_diag, 
    arma::vec* sup_diag,
    arma::vec* g_vec
) 
{
    sub_diag->resize(n_steps-2);
    main_diag->resize(n_steps-1);
    sup_diag->resize(n_steps-2);

    sub_diag->fill(-1);
    main_diag->fill(2);
    sup_diag->fill(-1);

    g_vec->resize(n_steps-1);

    double step_size = 1./ (double) n_steps;
    for (int i=0; i < n_steps-1; i++) {
        (*g_vec)(i) = analytic_solution((i+1)*step_size);
    }

}

int main() 
{

    arma::vec sub_diag, main_diag, sup_diag, g_vec;
    int n_steps;
    // Timing
    for (int i=1; i <= 6; i++) {
        n_steps = std::pow(10, i);
      

        // construct arrays

    }

}