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

using namespace std;

/*
 Sistem se transformise u njemu ekvivalentan tako sto se Xi' zamenjuje sa Vi, a Xi'' sa Vi'
 i onda izraze Vi' preko Xi i Vi, pri cemu je neophodno dodati jos i jednacine Xi' = Vi.
 Dakle, sistem glasi

 X1' = V1,
 X2' = V2,
 X3' = V3,
 X4' = V4,
 V1' = (68.67 - Fg - 6000(X1 - X3) - 5000(X1 - X2) - 250(V1 - V3) - 600(V1 - V2))/7,
 V2' = (58.86 + 5000(X1 - X3) - 9000(X2 - X3) - 25000(X3 - X4) - 250(V1 - V3) - 1700(V3 - V4))/6,
 V3' = (137.34 - 6000(X1 - X3) - 9000(X2 - X3) - 25000(X3 - X4) - 250(V1 - V3) + 1700(V3 - V4))/14;
 V4' = (441.45 + 25000(X3 - X4) + 1700(V3 - V4))/7.
 */

typedef long double real;
typedef real vec[4][2];

/*
 Ako je promenljiva A tipa vec, onda je

 A[0][0] = X1,
 A[1][0] = X2,
 A[2][0] = X3,
 A[3][0] = X4,
 A[0][1] = V1,
 A[1][1] = V2,
 A[2][1] = V3,
 A[3][1] = V4
 */

void f(real t, vec arg, vec &der) {
/*
 Parametar t je visak, jer ne figurise u sistemu. Ovde je stavljen samo da bi se skolski napisao kod,
 tj. onako kako izgleda u opstem slucaju.

 arg je vektor koji se sastoji od X1, X2, X4, X4, V1, V2, V3 i V4,
 dok je der vektor koji se sastoji od X1', X2', X3', X4', V1', V2', V3' i V4'.

 Ova funkcija izracunava der na osnovu t i arg
 */
    const real Ac = 2;
    const real a = 1e6;
    const real b = 1.56;
    const real c = 2e4;
    const real d = 0.73;
    const real e = 1;
    real Fg = Ac*(a*exp(b*log(arg[0][0])) + c*exp(d*log(arg[0][0]))*exp(e*log(arg[0][1])));

    for (int i = 0; i < 4; ++i) {
        der[i][0] = arg[i][1];
    }

    der[0][1] = (68.67 - Fg
                 - 6000*(arg[0][0] - arg[2][0])
                 - 5000*(arg[0][0] - arg[1][0])
                 - 250*(arg[0][1] - arg[2][1])
                 - 600*(arg[0][1] - arg[1][1]))/7;
    der[1][1] = (58.86
                 + 5000*(arg[0][0] - arg[2][0])
                 - 9000*(arg[1][0] - arg[2][0])
                 + 600*(arg[0][1] - arg[2][1]))/6;
    der[2][1] = (137.34
                 - 6000*(arg[0][0] - arg[2][0])
                 - 9000*(arg[1][0] - arg[2][0])
                 - 25000*(arg[2][0] - arg[3][0])
                 - 250*(arg[0][1] - arg[2][1])
                 - 1700*(arg[2][1] - arg[3][1]))/14;
    der[3][1] = (441.45
                 + 25000*(arg[2][0] - arg[3][0])
                 + 1700*(arg[2][1] - arg[3][1]))/45;
}

void combine(vec x, vec k, real h, vec &result) {
/*
 Pomocna funkcija koja racuna result kao x + h*k
 */
    for (int i = 0; i < 4; ++i) {
        for (int j = 0; j < 2; ++j) {
            result[i][j] = x[i][j] + h*k[i][j];
        }
    }
}

void runge_kutta(real t0, real t1, vec x, real h) {
/*
 Funkcija koja opisuje Runge-Kuta algoritam
 */
    ofstream out("rezultati.txt");

    for (real t = t0; t <= t1; t += h) {
        vec k1, k2, k3, k4, tmp;

        for (int i = 0; i < 4; ++i) {
            out << "X" << i + 1 << "(" << t << ") = " << x[i][0] << ",\t";
        }

        out << ";\n";

        for (int i = 0; i < 4; ++i) {
            out << "X\'" << i + 1 << "(" << t << ") = " << x[i][1] << ",\t";
        }

        out << ".\n\n";

        f(t, x, k1);
        combine(x, k1, h/2, tmp);
        f(t + h/2, tmp, k2);
        combine(x, k2, h/2, tmp);
        f(t + h/2, tmp, k3);
        combine(x, k3, h, tmp);
        f(t + h, tmp, k4);

        for (int i = 0; i < 4; ++i) {
            for (int j = 0; j < 2; ++j) {
                x[i][j] += (k1[i][j] + 2*k2[i][j] + 2*k3[i][j] + k4[i][j])*h/6;
            }
        }
    }
}

int main() {
/*
 Glavna funkcija
 */
    real t0, t1, h;
    vec x;

    cout << "Unesi pocetnu vrednost za t " << flush;
    cin >> t0;
    cout << "Unesi krajnju vrednost za t " << flush;
    cin >> t1;
    cout << "Unesi korak " << flush;
    cin >> h;
    cout << "Unesi pocetni uslov\n" << endl;

    for (int i = 0; i < 4; ++i) {
        cout << "X" << i + 1 << "(" << t0 << ") = " << flush;
        cin >> x[i][0];
    }

    for (int i = 0; i < 4; ++i) {
        cout << "X'" << i + 1 << "(" << t0 << ") = " << flush;
        cin >> x[i][1];
    }

    runge_kutta(t0, t1, x, h);

    return EXIT_SUCCESS;
}