#include <iostream>
#include <string>
#include <cmath>
#include <fstream>
#include <cstdlib>

using namespace std;
double E = 2.718281828459;
double Ka = 0.00159833;
double KET = 0.0959 / 60;

double rungeKutta2();
double rungeKutta4();
double eulerSimples();

void metodoBissecao();
void metodoCorda();
void picardPeano();

double rungeKutta4Sistema();

double funcao1(int minutos, double concAtual);
double funcaoMi(int minutos, double mi);
double funcaoMp(int minutos, double mi, double mp);

double volumePlasma = 3.5; // Em mL
double ket = 0.00159833333 * volumePlasma; // constante de eliminacao por minuto em mL/min

int main() {

    cout
            << "Calculo dos valores associados à administracao de de ceftazidima pela via intravenosa"
            << endl << endl;
    cout << "Que metodo numerico deseja utilizar para o calculo?" << endl
            << endl;
    cout << "1: Runge-Kutta segunda ordem (Concentracao, modelo linear)"
            << endl;
    cout << "2: Runge-Kutta quarta ordem (Concentracao, modelo linear)" << endl;
    cout << "3: Euler Simples (Concentracao, modelo linear)" << endl;

    cout << "4: Metodo Bissecao  (Cálculo de Ka) " << endl;
    cout << "5: Metodo Corda (Cálculo de Ka) " << endl;
    cout << "6: Metodo Picard-Peano (Cálculo de Ka) " << endl;

    cout
            << "7: Runge-Kutta quarta ordem para sistemas de equacoes(Concentracao, modelo bi-compartimental)"
            << endl;

    int selecao;
    cin >> selecao;
    if (cin.fail()) {
        cout << "Por favor introduza um inteiro.";
        system("pause");
        selecao = 0;
    }
    switch (selecao) {
    case 1:
        rungeKutta2();
        break;
    case 2:
        rungeKutta4();
        break;
    case 3:
        eulerSimples();
        break;
    case 4:
        metodoBissecao();
        break;
    case 5:
        metodoCorda();
        break;
    case 6:
        picardPeano();
        break;
    case 7:
        rungeKutta4Sistema();
        break;
    default:
        main();
    }

    system("pause");
    main();
    return 0;
}

double rungeKutta2() {
    int minutos = 0; //começa em t = 0
    int nAdministracoes = 0; //começa com a primeira administraçao
    double concAtual = 0; //Inicializar concentraçao
    double h = 1; //Passo: de 1 em 1 minuto

    ofstream file;
    file.open("Runge-Kutta2.txt");

    cout << "Metodo de Runge-Kutta 2 ordem" << endl << endl;
    file << "[";

    while (nAdministracoes < 6) {

        if ((minutos % (60 * 12)) == 0) {
            nAdministracoes++;
        }
        minutos++;
        concAtual += funcao1(minutos,
                concAtual + h * funcao1(minutos, concAtual) / 2);
        file << "[" << minutos << ", " << concAtual << "], ";
    }
    file << "]";
    cout << "Concentracao ao fim de 3 dias: " << concAtual << endl;
    return concAtual;

}
double rungeKutta4() { ///Runge-Kutta 4ªordem
    int minutos = 0; //começa em t = 0
    int nAdministracoes = 0; //começa com a primeira administraçao
    double concAtual = 0; //Inicializar concentraçao
    double h = 1; //Passo: de 1 em 1 minuto
    double xn = 0;
    double d1, d2, d3, d4;

    ofstream file;
    file.open("Runge-Kutta4.txt");
    cout << "Metodo de Runge-Kutta 4 ordem" << endl << endl;
    file << "[";

    while (nAdministracoes < 6) {

        if ((minutos % (60 * 12)) == 0) {
            nAdministracoes++;
        }
        minutos++;
        d1 = h * funcao1(minutos, concAtual);
        d2 = h * funcao1(minutos + (h / 2), concAtual + (d1 / 2));
        d3 = h * funcao1(minutos + (h / 2), concAtual + (d2 / 2));
        d4 = h * funcao1(minutos + h, concAtual + (d3));
        xn += h;
        concAtual += (d1 / 6) + (d2 / 3) + (d3 / 3) + (d4 / 6);
        file << "[" << minutos << ", " << concAtual << "], ";
    }
    file << "]";
    cout << "Concentracao ao fim de 3 dias: " << concAtual << endl;
    return concAtual;
}

double eulerSimples() {
    int minutos = 0; //começa em t = 0
    int nAdministracoes = 1; //começa com a primeira administraçao
    double concAtual = 0; //Inicializar concentraçao
    double h = 1; //Passo: de 1 em 1 minuto

    ofstream file;
    file.open("Euler Simples.txt");

    cout << "Metodo de Euler Simples" << endl << endl;
    file << "[";

    while (nAdministracoes < 6) {

        if ((minutos % (60 * 12)) == 0) {
            nAdministracoes++;
        }
        minutos++;

        concAtual += h * funcao1(minutos, concAtual);
        file << "[" << minutos << ", " << concAtual << "], ";
    }
    file << "]";
    cout << "Concentracao ao fim de 3 dias: " << concAtual << endl;
    return concAtual;

}
void metodoBissecao() {
    // Dados
    double Ket = (0.0959 / 60);
    double Ke = Ket;
    double tMax = 5;

    ofstream file;
    file.open("Bissecao.txt");
    // Valores Inicias
    double a, b;

    // Valor Medio
    double xn;

    // Valor da Funcao
    double fxn, fa, fb;

    cout << "Valores Iniciais" << endl;
    cout << "A: ";
    cin >> a;
    cout << "B: ";
    cin >> b;

    cout << endl << "A | B | Xn | F(Xn) | F(a) | F(b)" << endl;

    file << "[";

    double dif;
    do {
        xn = (a + b) / 2;

        fa = a * pow(E, -a * tMax) - Ke * pow(E, -Ke * tMax);
        fb = b * pow(E, -b * tMax) - Ke * pow(E, -Ke * tMax);
        fxn = xn * pow(E, -xn * tMax) - Ke * pow(E, -Ke * tMax);

        cout << a << " | " << b << " | " << xn << " | " << fxn << " | " << fa
                << " | " << fb << endl;

        file << "[" << xn << ", " << fxn << "], ";

        dif = abs(a - xn);

        if (abs(fxn) < abs(fa))
            a = xn;
        else
            b = xn;

        if (fxn < 0)
            fxn = 0 - fxn;

    } while (dif >= 0.0001);
    file << "]";

    cout << endl << " XN: " << xn << endl;
    cout << endl << " FXN: " << fxn << endl;

}
void metodoCorda() {
    // Dados
    double Ket = (0.0959 / 60);
    double Ke = Ket;
    double tMax = 5;

    ofstream file;
    file.open("Corda.txt");
    // Valores Inicias
    double a, b;

    // Valor Medio
    double xn;

    // Valor da Funcao
    double fxn, fa, fb;

    cout << "Valores Iniciais" << endl;
    cout << "A: ";
    cin >> a;
    cout << "B: ";
    cin >> b;

    cout << endl << "A | B | Xn | F(Xn) | F(a) | F(b)" << endl;

    file << "[";

    double dif;
    do {

        fa = a * pow(E, -a * tMax) - Ke * pow(E, -Ke * tMax);
        fb = b * pow(E, -b * tMax) - Ke * pow(E, -Ke * tMax);

        xn = (a * fb - b * fa) / (fb - fa);

        fxn = xn * pow(E, -xn * tMax) - Ke * pow(E, -Ke * tMax);

        cout << a << " | " << b << " | " << xn << " | " << fxn << " | " << fa
                << " | " << fb << endl;

        file << "[" << xn << ", " << fxn << "], ";

        dif = abs(a - xn);

        if (abs(fxn) < abs(fa))
            a = xn;
        else
            b = xn;

        if (fxn < 0)
            fxn = 0 - fxn;

    } while (dif >= 0.0001);
    file << "]";

    cout << endl << " XN: " << xn << endl;
    cout << endl << " FXN: " << fxn << endl;

}
void picardPeano() {
    const double KET = 0.0959 / 60;        //min^-1
    //const long double VAP = 3.5;                //L
    const double TMAX = 5;                 //min
    //const long double KE = KET * VAP;           //L^-1*min^-1
    //const long double nep = 2.7182818284590452353602874;
    double Ka;
    double Ka2 = 0.2;              //Guess
    double g;

    do {
        Ka = Ka2;
        //g = (KET*exp(-TMAX*KET)) / exp(-Ka*TMAX);
        g = -(log(KET) - KET * TMAX - log(Ka)) / TMAX;
        cout << "Ka= " << Ka << " ; g(Ka)= " << g << endl;
        Ka2 = g;

    } while (abs(Ka2 - Ka) > 0.0001);
}

double rungeKutta4Sistema() { ///Runge-Kutta 4ªordem
    int minutos = 0; //começa em t = 0
    int nAdministracoes = 0; //começa com a primeira administraçao
    double mi = 0; //Inicializar concentraçao
    double mp = 0;
    double h = 1; //Passo: de 1 em 1 minuto
    double d1, d2, d3, d4;
    double l1, l2, l3, l4;

    ofstream file1;
    ofstream file2;
    ofstream file3;
    file1.open("Runge-Kutta4 Mi.txt");
    file2.open("Runge-Kutta4 Mp.txt");
    file3.open("Runge-kutta4 conc.txt");
    cout << "Metodo de Runge-Kutta 4 ordem para sistemas de equações lineares"
            << endl << endl;
    file1 << "[";
    file2 << "[";
    file3 << "[";

    while (nAdministracoes < 6) {
        if ((minutos % (60 * 12)) == 0) {
            nAdministracoes++;
        }
        minutos++;
        d1 = h * funcaoMi(minutos, mi);
        d2 = h * funcaoMi(minutos + (h / 2), mi + (d1 / 2));
        d3 = h * funcaoMi(minutos + (h / 2), mi + (d2 / 2));
        d4 = h * funcaoMi(minutos + h, mi + d3);

        l1 = h * funcaoMp(minutos, mi, mp);
        l2 = h * funcaoMp(minutos + (h / 2), mi + (d1 / 2), mp + (l1 / 2));
        l3 = h * funcaoMp(minutos + (h / 2), mi + (d2 / 2), mp + (l2 / 2));
        l4 = h * funcaoMp(minutos + h, mi + d3, mp + l3);

        mi += (d1 / 6) + (d2 / 3) + (d3 / 3) + (d4 / 6);
        mp += (l1 / 6) + (l2 / 3) + (l3 / 3) + (l4 / 6);

        file1 << "[" << minutos << ", " << mi << "], ";
        file2 << "[" << minutos << ", " << mp << "], ";

        cout << "t= " << minutos << " ;mi= " << mi << " ;mp= " << mp << endl;
        file2 << "[" << minutos << ", " << mp << "], ";
        file1 << "[" << minutos << ", " << mi << "], ";
        file3 << "[" << minutos << ", " << (mp / volumePlasma) << "], ";
    }
    file1 << "]";
    file2 << "]";
    file3 << "]";
    return 0;
}

double funcao1(int minutos, double concAtual) {
    int dt = 0;
    if ((minutos % (60 * 12)) == 0)
        dt = 250;
    if (minutos == 0)
        dt = 250;
    return (dt - (ket * concAtual)) / volumePlasma;
}

double funcaoMi(int minutos, double mi) {
    int dt = 0;
    if ((minutos % (60 * 12)) == 0)
        dt = 250;

    if (minutos == 0)
        dt = 250;

    return (dt - (Ka * mi));
}

double funcaoMp(int minutos, double mi, double mp) {
    return (Ka * mi - KET * mp);
}