/////plik z klasą równania

#include "SecOrderNumDiffEq.h"

SecOrderNumDiffEq::SecOrderNumDiffEq(double (*p)(double), double (*r)(double), double (*s)(double),
                                     double a, double b,double h) : p(p), r(r), s(s), h(h), a(a), b(b) {
    n=(b-a)/h + 1;
}

void SecOrderNumDiffEq::setupConditions(double alpha, double beta, double gamma, double phi, double psi, double theta) {
    this->alpha = alpha;
    this->beta = beta;
    this->gamma = gamma;
    this->phi = phi;
    this->psi = psi;
    this->theta = theta;
}

void SecOrderNumDiffEq::prepareData() {
    vector<double> u = vector<double>(n-1);
    vector<double> d = vector<double>(n);
    vector<double> l = vector<double>(n);
    right = vector<double>(n);

    u[0] = alpha/h;
    d[0] = beta - alpha/h;
    d[n-1] = phi/h + psi;
    l[n-1] = -phi/h;
    right[0] = -gamma;
    right[n-1] = -theta;
    double xi = a+h;
    for(int i=1; i<n-1; i++)
    {
        l[i] = 1/h/h + r(xi-h)/12;
        u[i] = 1/h/h + r(xi+h)/12;
        d[i] = 10*r(xi)/12 - 2/h/h;
        right[i] = -(s(xi-h) + 10*s(xi) + s(xi+h))/12;
        xi+=h;
    }
    l.erase(l.begin());
    T = ThomasMatrix(u,d,l);
}

vector<double> SecOrderNumDiffEq::solve() {
    prepareData();
    ThomasEquation eq = ThomasEquation(T);
    vector<double> y = eq.solveXfor(right);

    return y;
}

unsigned long SecOrderNumDiffEq::getN() const {
    return n;
}






//plik main
#include <iostream>
#include <cmath>
#include <fstream>
#include "ThomasEquation.h"
#include "SecOrderDiffEq.h"
#include "SecOrderNumDiffEq.h"

using namespace std;
void display(vector<double> x);
double u(double x)
{
    return (2*x*cos(2*x)+2*sin(2*x)-log(2)*sin(2*x)-2*log(cos(x))*sin(2*x))/4;
}

void command(FILE *pipe,string cmd)
{
    cmd+="\n";
    fprintf(pipe,cmd.c_str());
    fflush(pipe);
}

int main() {
    double h=1e-7,a=0, b=M_PI_4;

    ofstream bledy;
    bledy.open("bledy.csv");
    vector<vector<double>> rozwClassic = vector<vector<double>>(14,vector<double>(100,0));//wektor na rozwiązania numeryczne metodą klasyczną
    vector<vector<double>> rozwNumerow = vector<vector<double>>(14,vector<double>(100,0));//wektor na rozwiązania numeryczne metodą numerowa

    for(int j=0; j<14; j+=2) {
        //rozwiązanie schematem klasycznym
        /*SecOrderDiffEq classicEQ([](double x) -> double { return 1.0; }, [](double x) -> double { return 0.0; },
                               [](double x) -> double { return 4.0; }, [](double x) -> double { return tan(x); }, a, b, h);
        classicEQ.setupConditions(0.0, 1.0, 0.0, 0.0, 1.0, -1.0 / 2.0);
        vector<double> yClassic = classicEQ.solve();
*/
        //rozwiązanie schematem numerowa
        SecOrderNumDiffEq numerowEQ([](double x) -> double { return 1.0; }, [](double x) -> double { return 4.0; },
                                   [](double x) -> double { return tan(x); }, a, b, h);
        numerowEQ.setupConditions(0.0, 1.0, 0.0, 0.0, 1.0, -1.0 / 2.0);
        vector<double> yNumerow = numerowEQ.solve();

        int n = numerowEQ.getN(),k=0,step;
        if(n > 100)
            step=n/100+1;
        else
            step=1;
        double xi = a;
        //double maxC = fabs(u(xi) - yClassic[0]);
        double maxN = fabs(u(xi) - yNumerow[0]);
        for (int i = 1; i < n; i++) {
            xi += h;
            //maxC = (maxC > (fabs(u(xi) - yClassic[i]))) ? maxC : fabs(u(xi) - yClassic[i]);
            maxN = (maxN > (fabs(u(xi) - yNumerow[i]))) ? maxN : fabs(u(xi) - yNumerow[i]);
            if(i%step == 0){
                //rozwClassic[j][k]=xi;
                //rozwClassic[j+1][k]=yClassic[i];
                rozwNumerow[j][k]=xi;
                rozwNumerow[j+1][k]=yNumerow[i];
                k++;
            }
        }
        bledy<<log10(h)<<" "/*<<log10(maxC)<<" "*/<<log10(maxN)<<endl;
        h*=10;
    }
    bledy.close();

    //wpisanie tablicy z rozwiązaniami numerycznymi do pliku
    //classic
    /*ofstream rozwC;
    rozwC.open("rozw.csv");
    for(int i=0; i<100; i++)
        rozwC<<rozwClassic[0][i]<<" "<<rozwClassic[1][i]<<" "<<rozwClassic[2][i]<<" "<<rozwClassic[3][i]<<" "<<rozwClassic[4][i]<<" "
            <<rozwClassic[5][i]<<" "<<rozwClassic[6][i]<<" "<<rozwClassic[7][i]<<" "<<rozwClassic[8][i]<<" "<<rozwClassic[9][i]<<" "
            <<rozwClassic[10][i]<<" "<<rozwClassic[11][i]<<" "<<rozwClassic[12][i]<<" "<<rozwClassic[13][i]<<endl;
    rozwC.close();*/
    //numerow
    ofstream rozwN;
    rozwN.open("rozw2.csv");
    for(int i=0; i<100; i++)
        rozwN<<rozwNumerow[0][i]<<" "<<rozwNumerow[1][i]<<" "<<rozwNumerow[2][i]<<" "<<rozwNumerow[3][i]<<" "<<rozwNumerow[4][i]<<" "
             <<rozwNumerow[5][i]<<" "<<rozwNumerow[6][i]<<" "<<rozwNumerow[7][i]<<" "<<rozwNumerow[8][i]<<" "<<rozwNumerow[9][i]<<" "
             <<rozwNumerow[10][i]<<" "<<rozwNumerow[11][i]<<" "<<rozwNumerow[12][i]<<" "<<rozwNumerow[13][i]<<endl;
    rozwN.close();


    //rysowanie wykresu rozwiązania
    //classic
    FILE *plot_wykres = popen("gnuplot 2>/dev/null","w");
    command(plot_wykres,"set yrange [0:pi/4]");
    command(plot_wykres,"plot (2*x*cos(2*x)+2*sin(2*x)-log(2)*sin(2*x)-2*log(cos(x))*sin(2*x))/4 t 'analityczne',\
                         'rozw.csv' u 1:2 t '10^-7',\
                         'rozw.csv' u 3:4 t '10^-6',\
                         'rozw.csv' u 5:6 t '10^-5',\
                         'rozw.csv' u 7:8 t '10^-4',\
                         'rozw.csv' u 9:10 t '10^-3',\
                         'rozw.csv' u 11:12 t '10^-2',\
                         'rozw.csv' u 13:14 t '10^-1'");
    //numerow
    FILE *plot_wykres2 = popen("gnuplot 2>/dev/null","w");
    command(plot_wykres2,"set yrange [0:pi/4]");
    command(plot_wykres2,"plot (2*x*cos(2*x)+2*sin(2*x)-log(2)*sin(2*x)-2*log(cos(x))*sin(2*x))/4 t 'analityczne',\
                         'rozw2.csv' u 1:2 t '10^-7',\
                         'rozw2.csv' u 3:4 t '10^-6',\
                         'rozw2.csv' u 5:6 t '10^-5',\
                         'rozw2.csv' u 7:8 t '10^-4',\
                         'rozw2.csv' u 9:10 t '10^-3',\
                         'rozw2.csv' u 11:12 t '10^-2',\
                         'rozw2.csv' u 13:14 t '10^-1'");


    //rysowanie wykresu błędu
    FILE *plot_bledy = popen("gnuplot 2>/dev/null", "w");
    command(plot_bledy,"plot 'bledy.csv' u 1:2 w lp t 'klasyczny','bledy.csv' u 1:3 w lp t 'numerowa'");


    //zablokowanie terminala żeby zobaczyć wykres
    cout<<"nacisnij enter"<<endl;
    cin.get();
    pclose(plot_wykres);
    pclose(plot_wykres2);
    pclose(plot_bledy);
    return 0;
}