#include "stdafx.h"
#include <iostream>
#include <cstdlib>
#include <math.h>

using namespace std;

double nu = 1 * 1e-6;
double y = 73 * 1e-3;
double ro = 1e3;

double Patm = 101325;

double fr = 26500;
double w = 2 * 3.1415 * fr;

double P0 = Patm;

double r0 = 0.1;
double h = r0 / 8.86;

double dt = 1e-6;

double c = 1481;
double k = 1.4;

double Pdr = 1.5 * P0;

double P(double t)
{
    //return -Pdr * cos(w * t);

    return sin(w*t);
}

double p_gas(double r, double t)
{
    return (P0 + 2 * y / r0) * pow((pow(r0, 3) - pow(h, 3)) / (pow(r, 3) - pow(h, 3)), k);
}

double f(double t, double r, double v)
{
    return v;
}

double g(double t, double r, double v)
{
    return (1 / r) * ((1 / ro) * (p_gas(r, t) - P(t) - P0 + (r / c) * (p_gas(r, t + dt) - p_gas(r, t)) / dt - 4 * nu * v / r - 2 * y / r) - 3 * v * v / 2);
}

int main()
{
    freopen("out.txt", "w", stdout);

    double cc = 5;

    double tmin = 0, tmax = 1 * cc;

    double N = 1e3 * cc;
    double h = (tmax - tmin) / N;

    double t0 = tmin;
    double v0 = 0;

    double k[4], l[4];

    //cout << "t\tv\tr\tp" << endl;

    for (int i = 0; i <= N; i++)
    {
        cerr << i << endl;

        cout << P(t0) << endl;
        //cout << t0 << "\t" << v0 << "\t" << r0 << "\t" << P(t0) << endl;

        k[0] = h * f(t0, r0, v0);
        l[0] = h * g(t0, r0, v0);

        k[1] = h * f(t0 + h / 2, r0 + k[0] / 2, v0 + l[0] / 2);
        l[1] = h * g(t0 + h / 2, r0 + k[0] / 2, v0 + l[0] / 2);

        k[2] = h * f(t0 + h / 2, r0 + k[1] / 2, v0 + l[1] / 2);
        l[2] = h * g(t0 + h / 2, r0 + k[1] / 2, v0 + l[1] / 2);

        k[3] = h * f(t0 + h / 2, r0 + k[2] / 2, v0 + l[2] / 2);
        l[3] = h * g(t0 + h / 2, r0 + k[2] / 2, v0 + l[2] / 2);

        t0 += h;
        r0 += (k[0] + 2 * k[1] + 2 * k[2] + k[3]) / 6;
        v0 += (l[0] + 2 * l[1] + 2 * l[2] + l[3]) / 6;
    }

    return 0;
}