#include <iostream>
#include <fstream>
#include <math.h>
#include <vector>

using namespace std;

ofstream fout;

double Uo;
double w;
double B;
double S;
double R;
double J;
double to;
double t1;
int n;
double fo;
double eps;
double dfo;

void decay_time(double *t, double *f, int len)
{
    double buf_f = f[0], buf1_t, buf2_t;
    int flg = 0;
    for (int i = 1; i < len - 1; i++)
    {
        if (fabs(f[i]) > fabs(f[i-1]) and fabs(f[i]) > fabs(f[i+1]))
        {
            if (flg == 0)
            {
                cout << f[i] << " " << t[i] << endl;
                buf_f = f[i];
                buf1_t = t[i];
                flg = 1;
            }
            else
            {
                if (fabs(buf_f) / fabs(f[i]) >= M_E)
                {
                    cout << f[i] << " " << t[i] << endl;
                    buf2_t = t[i];
                    break;
                }
            }

        }
        else if (i == 1 and fabs(f[i-1]) > fabs(f[i]))
        {
            cout << f[i-1] << " " << t[i-1] << endl;
            buf_f = f[i-1];
            buf1_t = t[i-1];
            flg = 1;
        }
    }
    cout << "decay time: " << buf2_t - buf1_t << endl;
}

double func(double f, double t,  double df)
{
    return (-Uo*cos(w*t) - B*S*sin(f)*df)*S*B*sin(f)/(R*J);
}

void RungeStep(double* df, double* f, double* t, int N)
{
    double k1, k2, k3, k4;
    double q1, q2, q3, q4;

    df[0]=dfo;
    f[0]=fo;
    t[0]=to;

    double dt = (double) (t1 - to)/(N-1);
    for (int i = 1; i < N; i++)
    {
        q1 = dt * func(f[i-1]         , t[i-1]         , df[i-1]         );
        k1 = dt * (df[i-1]         );

        q2 = dt * func(f[i-1] + k1/2.0, t[i-1] + dt/2.0, df[i-1] + q1/2.0);
        k2 = dt * (df[i-1] + q1/2.0);

        q3 = dt * func(f[i-1] + k2/2.0, t[i-1] + dt/2.0, df[i-1] + q2/2.0);
        k3 = dt * (df[i-1] + q2/2.0);

        q4 = dt * func(f[i-1] + k3    , t[i-1] + dt    , df[i-1] + q3    );
        k4 = dt * (df[i-1] + q3    );

        df[i] = df[i-1] + (q1 + 2.0*q2 + 2.0*q3 + q4)/6.0;
        f[i]  = f[i-1]  + (k1 + 2.0*k2 + 2.0*k3 + k4)/6.0;
        t[i]  = t[i-1]  + dt;
    }
}

double maxNorm(double* f1, double* f2, double* t1, double* t2, int n1)
{
    double max = -1.0;
    double diff;
    for (int i = 0; i < n1; i++)
    {
        diff = fabs( f1[i]-f2[2*i]);
        if (diff > max) max=diff;
    }
    return max;
}

void Runge(bool flag)
{
    int n1=n;
    int n2=n;
    vector<double>arr_df1, arr_f1, arr_t1;
    vector<double>arr_df2, arr_f2, arr_t2;

    arr_df1.resize(n1, 0.0);
    arr_f1.resize(n1, 0.0);
    arr_t1.resize(n1, 0.0);
    RungeStep(&arr_df1[0], &arr_f1[0], &arr_t1[0], n1);

    double diff=2.0*eps;
    while (eps<diff)
    {
        n2 = 2 * n1 - 1;
        arr_df2.resize(n2, 0.0);
        arr_f2.resize(n2, 0.0);
        arr_t2.resize(n2, 0.0);

        RungeStep(&arr_df2[0], &arr_f2[0], &arr_t2[0], n2);
        diff = maxNorm(&arr_f1[0], &arr_f2[0],  &arr_t1[0], &arr_t2[0], n1);
        diff/= 15.0;
        cout<< n1 <<" "<< n2 <<" "<< diff <<endl;

        n1 = n2;
        arr_df1 = arr_df2;
        arr_f1 = arr_f2;
        arr_t1 = arr_t2;
    }
    cout <<"The Runge's step suitable for the given epsilon="<<eps<<" is "<< (double) (t1 - to)/(n2-1) <<endl;

    decay_time(&arr_t2[0], &arr_f2[0], n2);

    if(flag)
    {
        ofstream fout;
        fout.open("Runge_Kutt.txt");
        for (int i = 0; i < n2; i++)
            fout << arr_t2[i] << "\t" << arr_f2[i] << "\t" << arr_df2[i] << endl;
        fout.close();
    }

    arr_df1.clear();
    arr_df2.clear();
    arr_f1.clear();
    arr_f2.clear();
    arr_t1.clear();
    arr_t2.clear();

}

void EulerStep(double* df, double* f, double* t, int N)
{
    df[0]=dfo;
    f[0]=fo;
    t[0]=to;
    double dt = (double) (t1 - to)/(N-1);
    for (int i = 1; i < N; i++)
    {
        df[i]=df[i-1]  + dt*func(f[i-1], t[i-1], df[i-1]);
        f[i]=f[i-1] + dt*df[i];
        t[i]=t[i-1] + dt;
    }
}

void Euler(bool flag)
{
    int n1=n;
    int n2=n;
    vector<double>arr_df1, arr_f1, arr_t1;
    vector<double>arr_df2, arr_f2, arr_t2;

    arr_df1.resize(n1, 0.0);
    arr_f1.resize(n1, 0.0);
    arr_t1.resize(n1, 0.0);
    EulerStep(&arr_df1[0], &arr_f1[0], &arr_t1[0], n1);

    double diff=2.0*eps;
    while (eps<diff)
    {
        n1 = n2;
        n2 = 2 * n1 - 1;
        arr_df1 = arr_df2;
        arr_f1 = arr_f2;
        arr_t1 = arr_t2;
        arr_df2.resize(n2, 0.0);
        arr_f2.resize(n2, 0.0);
        arr_t2.resize(n2, 0.0);
        EulerStep(&arr_df2[0], &arr_f2[0], &arr_t2[0], n2);
        diff = maxNorm(&arr_f1[0], &arr_f2[0],  &arr_t1[0], &arr_t2[0], n1);
        diff/= 3.0;
        cout<< n1 <<" "<< n2 <<" "<< diff <<endl;
    }
    cout <<"The Euler's step suitable for the given epsilon="<<eps<<" is "<< (double) (t1 - to)/(n2-1) <<endl;

    decay_time(&arr_t2[0], &arr_f2[0], n2);

    if(flag)
    {
        ofstream fout;
        fout.open("Euler.txt");
        for (int i = 0; i < n2; i++)
            fout << arr_t2[i] << "\t" << arr_f2[i] << "\t" << arr_df2[i] << endl;
        fout.close();
    }

    arr_df1.clear();
    arr_df2.clear();
    arr_f1.clear();
    arr_f2.clear();
    arr_t1.clear();
    arr_t2.clear();

}

void InitDef()
{
    Uo=10;
    w=0;
    B=5;
    S=0.01;
    R=100;
    J=0.001;
    to=0;
    t1=1000;
    n=21;
    fo=M_PI/2 - 0.1;
    dfo=0;
    eps=0.001;
    //   koef = 0.0;
}

/* void InitUser()
{
    cout <<"Enter amplitude of voltage: ";
    cin >> Uo;
    cout<<"Enter frequency:   ";
    cin>>w;
    cout<<"Enter B:  ";
    cin>>B;
    cout<<"Enter S:  ";
    cin>>S;
    cout<<"Enter resistance:  ";
    cin>>R;
    cout<<"Enter moment of inertia:  ";
    cin>>J;
    cout<<"Enter t0:  ";
    cin>>to;
    cout<<"Enter t1:  ";
    cin>>t1;
    cout<<"Enter initial angle:  ";
    cin>>fo;
    cout<<"Enter epsilon:  ";
    cin>>eps;
    cout<<"Enter angular velocity:  ";
    cin>>dfo;
    cout <<"Thank you!"<< endl;
}
 */

int main()
{
    bool flag = true; //true - print to files, false - do not print
    InitDef();
  //  InitUser();
    Runge(flag);
    Euler(flag);
    return 0;
}