#include <iostream>
#include<cmath>
#include<iomanip>
#include<windows.h>
#include<stdio.h>
#include<fstream>


using namespace std;

double funPME(double yk1, double yk, double dt, double tk)
{
    return (yk1 - yk) / dt + (10 * tk * tk + 20) / (tk * tk + 1) * (yk1 - 1);
}

double funPMT(double yk1, double yk, double dt, double tk)
{
    return ((yk1 - yk) / dt) - (-1 * (((10 * tk * tk + 20) / (tk * tk + 1)) * (yk - 1)) - (10 * (tk + dt) * (tk + dt) + 20) / ((tk + dt) * (tk + dt) + 1) * (yk1 - 1)) / 2;
}

double fun(double tk, double yk)
{
    return -1 * ((10 * tk * tk  + 20) / (tk * tk + 1) * (yk - 1));
}

double CorrecrValue(double tk)
{
    return 1 - exp( -10 * (tk + atan(tk)));
}

double bisekcja(double a, double b, double tolerance, double yk, double dt, double tk, int ops = 2)
{
    double tmp, znak[2] = {0, 0}, f[2], ftmp, x, fx;
    int i = 0;

    //printf("\nMetoda Bisekcji.\n");

    f[0] = funPME(a, yk, dt, tk);
    f[1] = funPME(b, yk, dt, tk);

    x = (a + b) / 2; //przyblizenie pierwiastka
    tmp = (fabs(b - a)) / 2; //estymator bledu
    fx = funPME(x, yk, dt, tk); //wartosc funkcji w punkcie podzialu

    if (f[0] < 0) znak[0] = -1;
    else if (f[0] > 0) znak[0] = 1;

    if (f[1] < 0) znak[1] = -1;
    else if (f[1] > 0) znak[1] = 1;

    if (znak[0] != znak[1])
    {

        switch(ops)
        {
            case 0:
            {
                cout << "Kryterium estymatora bledu\n";

//                while (fabs(tmp) > est)
                {
                    tmp = ( b - a) / 2; //estymator bledu
                    ftmp = funPME(tmp, yk, dt, tk);
                    x = (a + b) / 2; //kolejne przyblizenia
                    fx = funPME(x, yk, dt, tk);
                    if (fx > 0)
                    {
                        if (f[0] < 0) b = x;
                        else if (f[1] < 0) a = x;
                    }

                    else if (fx < 0)
                    {
                        if (f[0] > 0) b = x;
                        else if (f[1] > 0) a = x;
                    }
                    cout<<"x"<<i++<<" = "<<x;
                    cout<<"\t estymator: "<<fabs(tmp)<<endl;


                }
            }break;

            case 1:
            {
                cout<<"Kryterium ograniczenia iteracji\n";

                //while (i++ < n)
                {

                    x = (a + b) / 2;
                    fx = funPME(x, yk, dt, tk);

                     if (fx > 0)
                    {
                        if (f[0] < 0) b = x;
                        else if (f[1] < 0) a = x;
                    }

                    else if (fx < 0)
                    {
                        if (f[0] > 0) b = x;
                        else if (f[1] > 0) a = x;
                    }
                    cout<<"x"<<i<<" = "<<x<<endl;

                }
            }break;

            case 2:
            {
                //cout << "Kryterium wiarygodnosci przyblizenia pierwiastka\n";

                while (fabs(funPME(x, yk, dt, tk)) > tolerance)
                {

                    x = (a + b) / 2;
                    fx = funPME(x, yk, dt, tk);

                    if (fx > 0)
                    {
                        if (f[0] < 0) b = x;
                        else if (f[1] < 0) a = x;
                    }

                    else if (fx < 0)
                    {
                        if (f[0] > 0) b = x;
                        else if (f[1] > 0) a = x;
                    }
                    //cout<<"x"<<i++<<" = "<<x;
                    //cout<<"\tResiduum rownania: "<<fabs(funPME(x, yk, dt, tk))<<endl;

                }
            }break;

            default:
            {
                printf("Bledny numer, program zostanie zakonczony\n");
                exit(0);
            }break;
        }
    }
    return x;
}

double bisekcja2(double a, double b, double tolerance, double yk, double dt, double tk, int ops = 2)
{
    double tmp, znak[2] = {0, 0}, f[2], ftmp, x, fx;
    int i = 0;
    ops = 2;

    //printf("\nMetoda Bisekcji.\n");

    f[0] = funPMT(a, yk, dt, tk);
    f[1] = funPMT(b, yk, dt, tk);

    x = (a + b) / 2; //przyblizenie pierwiastka
    tmp = (fabs(b - a)) / 2; //estymator bledu
    fx = funPMT(x, yk, dt, tk); //wartosc funkcji w punkcie podzialu

    if (f[0] < 0) znak[0] = -1;
    else if (f[0] > 0) znak[0] = 1;

    if (f[1] < 0) znak[1] = -1;
    else if (f[1] > 0) znak[1] = 1;


    if (znak[0] != znak[1])
    {

        switch(ops)
        {
            case 0:
            {
                cout << "Kryterium estymatora bledu\n";

//                while (fabs(tmp) > est)
                {
                    tmp = ( b - a) / 2; //estymator bledu
                    ftmp = funPME(tmp, yk, dt, tk);
                    x = (a + b) / 2; //kolejne przyblizenia
                    fx = funPME(x, yk, dt, tk);
                    if (fx > 0)
                    {
                        if (f[0] < 0) b = x;
                        else if (f[1] < 0) a = x;
                    }

                    else if (fx < 0)
                    {
                        if (f[0] > 0) b = x;
                        else if (f[1] > 0) a = x;
                    }
                    cout<<"x"<<i++<<" = "<<x;
                    cout<<"\t estymator: "<<fabs(tmp)<<endl;


                }
            }break;

            case 1:
            {
                cout<<"Kryterium ograniczenia iteracji\n";

                //while (i++ < n)
                {

                    x = (a + b) / 2;
                    fx = funPME(x, yk, dt, tk);

                     if (fx > 0)
                    {
                        if (f[0] < 0) b = x;
                        else if (f[1] < 0) a = x;
                    }

                    else if (fx < 0)
                    {
                        if (f[0] > 0) b = x;
                        else if (f[1] > 0) a = x;
                    }
                    cout<<"x"<<i<<" = "<<x<<endl;

                }
            }break;

            case 2:
            {
                //cout << "Kryterium wiarygodnosci przyblizenia pierwiastka\n";
                while (fabs(funPMT(x, yk, dt, tk)) > tolerance)
                {

                    x = (a + b) / 2;
                    fx = funPMT(x, yk, dt, tk);

                    if (fx > 0)
                    {
                        if (f[0] < 0) b = x;
                        else if (f[1] < 0) a = x;
                    }

                    else if (fx < 0)
                    {
                        if (f[0] > 0) b = x;
                        else if (f[1] > 0) a = x;
                    }
                    //cout<<"x"<<i++<<" = "<<x;
                    //cout<<"\tResiduum rownania: "<<fabs(funPME(x, yk, dt, tk))<<endl;

                }
            }break;

            default:
            {
                printf("Bledny numer, program zostanie zakonczony\n");
                exit(0);
            }break;
        }
    }
    return x;
}


void BME(double yk, double dt)
{
    fstream file, file2;

    file.open("BME_Stabilna.txt");
    file2.open("BME_blad.txt");

    file << "t\tyk\ty(t)\n";
    file2 << "log10|t|\t\tlog10|blad|\n";

    double tk = 0.0 + dt, yk1 = 0;
    int i = 0;

    while (i++ < 100)
    {
        yk1 = yk + fun(tk, yk) * dt;
        yk = yk1;
        tk += dt;
        cout << "przyblizenie nr "<< i << " = " << yk1 << endl;
        cout << "Dokladna wartosc w wezle = " << CorrecrValue(tk) << endl;
        cout << "Wezel przyblizenia tk = " << tk << endl;
        cout << "Blad przyblizenia = " << fabs(CorrecrValue(tk) - yk) << endl;
        cout << endl;

        file << tk << "\t" << CorrecrValue(tk) << endl;//"\t" << CorrecrValue(tk) << endl;
        file2 << log10(tk) << "\t\t" << log10(fabs(CorrecrValue(tk) - yk)) << endl;
    }

    file.close();
    file2.close();
}

void PME(double yk, double dt)
{
    double tk = 0.0 + dt, yk1 = 1;
    int i = 0;
    fstream file, file2;

    file.open("PME.txt");
    file2.open("PME_blad.txt");
    file2 << "log10|t|\t\tlog10|blad|\n";

    file << "t\ty(t)\t\tyt\n";


    for (i = 0; i < 100; i++)
    {
        yk1 = bisekcja(-189.0, 231.0, 1e-12, yk, dt, tk);
        cout << "Przyblizenie nr " << i + 1 << " Wynosi: " << yk1 << endl;
        cout << "Dokladna wartosc w wezle  = " << CorrecrValue(tk) << endl;
        cout << "Blad przyblizenia = " << fabs(CorrecrValue(tk) - yk1) << endl;
        cout << "Residuum rownania nieliniowego wynosi = " << fabs(funPME(yk1, yk, dt, tk)) << endl;
        cout << "Wezel: " << tk << endl << endl;

        file << setw(6) << tk << "\t" << CorrecrValue(tk) << endl;
        file2 << log10(tk) << "\t\t" << log10(fabs(CorrecrValue(tk) - yk)) << endl;

        yk = yk1;
        tk += dt;


    }

    file.close();
    file2.close();
}

void PMT(double yk, double dt)
{
    double tk = 0.0 + dt, yk1 = 1;
    int i = 0;
    fstream file, file2;

    file.open("PMT.txt");
    file2.open("PMT_blad.txt");

    file << "t\ty(t)\t\tyt\n";
    file2 << "log10|t|\t\tlog10|blad|\n";

    for (i = 0; i < 100; i++)
    {
        yk1 = bisekcja2(-137.0, 230.0, 1e-12, yk, dt, tk);
        cout << "Przyblizenie nr " << i + 1 << " Wynosi: " << yk1 << endl;
        cout << "Dokladna wartosc w wezle  = " << CorrecrValue(tk) << endl;
        cout << "Blad przyblizenia = " << fabs(CorrecrValue(tk) - yk1) << endl;
        cout << "Residuum rownania nieliniowego wynosi = " << fabs(funPMT(yk1, yk, dt, tk)) << endl;
        cout << "Wezel: " << CorrecrValue(tk) << endl << endl;

        file << setw(6) << tk << "\t" << yk1 << endl;
        file2 << log10(tk) << "\t\t" << log10(fabs(CorrecrValue(tk) - yk)) << endl;

        yk = yk1;
        tk += dt;


    }
    file.close();
    file2.close();
}

int main()
{
    double y = 0.0, dt = 0.0001;

    cout << "Bezposrednia Metoda Eulera: \n";
    BME(y, dt);

    cout << "Posrednia Metoda Eulera: \n";
    PME(y, dt);

    cout << "Posrednia Metoda Trapezow: \n";
    PMT(y, dt);


    return 0;
}