#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <stdlib.h>

using namespace std;

const double a = 0.05; //częstość narodzin ofiar lub współczynnik przyrostu ofiar
const double b = 0.002; //częstość umierania ofiar na skutek drapieżnictwa
const double c = 0.06; //częstość narodzin drapieżników lub współczynnik przyrostu drapieżników
const double d = 0.004; //częstość umierania drapieżników lub współczynnik ubywania drapieżników

const double h = 0.1; //UŻYWAM h ZAMIAST step

const int l_krokow = 4000;

double dxdt (double x, double y)
{
    return (a - b * y) * x;
}

double dydt (double x, double y)
{
    return (c * x - d) * y;
}

int main (int argc, char** argv)
{
    // x - ofiary, y - drapieżniki

    double y0h = 1;     // h zamiast step

    for(int i=0; i<100; i++)        //ten pierwszy for idzie po ilosciach drapieznikow (czyli    analizuje sie caly czas zmiany w czasie ale przy roznych startowych wartosciach drapieznikow)
    {
        double x0 = 100;    // ilosc ofiar
        double y0 = 10+i*y0h;        //zastanawiam sie skad taki wzor… (po co to mnozyc jeszcze razy 1…)

        double x=x0;
        double y=y0;

        ofstream af;    //  strumien do zapisu wsp X
        stringstream afN;   //strumien do utworzenia nazwy pliku
        afN << "./txt/xt" << setfill('0') << setw(4) << i << ".txt"; //tworzy plik w folderze z programem/txt/xt’iteracja’.txt
        string afName = afN.str();  // .str() zwraca string z przekopiowana zawartoscia

        ofstream bf;
        stringstream bfN;
        bfN << "./txt/yt" << setfill('0') << setw(4) << i << ".txt"; //tworzy plik w folderze z programem/txt/yt’iteracja’.txt
        string bfName = bfN.str();

        af.open(afName.c_str());
        bf.open(bfName.c_str());

        double kx1, ky1, kx2, ky2, kx3, ky3, kx4, ky4;

        for (double i1 = 0; i1 < l_krokow * h; i1 += h) // Funkcja zapisująca współrzędne do plików af i bf.
        {
            af << i1 << ' ' << x << endl;
            bf << i1 << ' ' << y << endl;

            x0 = x;
            y0 = y;
/*
            x += h * dxdt(x0,y0);
            y += h * dydt(x0,y0);
*/

            kx1 = h * dxdt(x0,y0);
            ky1 = h * dydt(x0,y0);

            kx2 = h * dxdt(x0 + (h/2.0), y0 + (kx1/2.0));
            ky2 = h * dydt(x0 + (h/2.0), y0 + (ky1/2.0));

            kx3 = h * dxdt(x0 + (h/2.0), y0 + (kx2/2.0));
            ky3 = h * dydt(x0 + (h/2.0), y0 + (ky2/2.0));

            kx4 = h * dxdt(x0 + h, y0 + kx3);
            ky4 = h * dxdt(x0 + h, y0 + ky3);

            x += (kx1 + 2.0*kx2 + 2.0*kx3 + kx4)/6.0;
            y += (ky1 + 2.0*ky2 + 2.0*ky3 + ky4)/6.0;


        }

        af.close();
        bf.close();


    }

    if(argc > 1)
    {
        system("rm ./png/* ./txt/*");
    }

    return 0;

}