#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>
#include <sstream>

// определяем будущие функции, чтобы написать весь основной код ниже main()
double Razgon(double H, double V1, double V2, double mass, double S);
double Podyom(double H1, double H2, double V, double mass, double S);
double RazgPod(double H1, double H2, double V1, double V2, double mass, double S);

double ro(double H);
double a(double H);
double Cy(double alpha);
double Cx(double alpha);
double P(double M);

// указываем количество двигателей
int engineCount = 2;

int main()
{
    int n;

    std::cout << "Введите количество элементарных разбиений: ";
    std::cin >> n;
    int m = n + 1;

    // начальные условия из варианта
    double Vn = 310 * 1000.0 / 3600.0;
    double Vk = 700 * 1000.0 / 3600.0;
    double Hn = 300.0;
    double Hk = 6000.0;
    double mass = 47000.0;
    double S = 127.0;
    // считаем шаг разбиения
    double Hel = (Hk - Hn) / n;
    double Vel = (Vk - Vn) / n;

    std::vector<std::vector<double>> Tr(m, std::vector<double>(n));
    std::vector<std::vector<double>> Tp(n, std::vector<double>(m));
    std::vector<std::vector<double>> Tpr(n, std::vector<double>(n));

    // заполнение матрицы разгона
    for (int i = 0; i < m; ++i)
    {
        double H_start = Hn + i * Hel;
        double H_end = Hn + (i + 1) * Hel;
        for (int j = 0; j < n; ++j)
        {
            double V_start = Vn + j * Vel;
            double V_end = Vn + (j + 1) * Vel;
            double traz = Razgon(H_start, V_start, V_end, mass, S);
            Tr[i][j] = traz;
        }
    }

    // заполнение матрицы подъема
    for (int i = 0; i < n; ++i)
    {
        double V_start = Vn + i * Vel;
        double V_end = Vn + (i + 1) * Vel;
        for (int j = 0; j < m; ++j)
        {
            double H_start = Hn + j * Hel;
            double H_end = Hn + (j + 1) * Hel;
            double tpod = Podyom(H_start, H_end, V_end, mass, S);
            Tp[i][j] = tpod;
        }
    }

    // заполнение матрицы разгона-подъема
    for (int i = 0; i < n; ++i)
    {
        double H_start = Hn + i * Hel;
        double H_end = Hn + (i + 1) * Hel;
        for (int j = 0; j < n; ++j)
        {
            double V_start = Vn + j * Vel;
            double V_end = Vn + (j + 1) * Vel;
            double tpodraz = RazgPod(H_start, H_end, V_start, V_end, mass, S);
            Tpr[i][j] = tpodraz;
        }
    }

    // вывод значений
    std::cout << std::fixed << std::setprecision(2);

    std::cout << "\nМатрица разгона [время в сек]:\n";
    std::cout << std::setw(20) << "H\\V";
    for (int j = 0; j < n; ++j)
    {
        double V_start = (Vn + j * Vel) * 3.6;
        double V_end = (Vn + (j + 1) * Vel) * 3.6;
        std::ostringstream oss;
        oss << V_start << "-" << V_end;
        std::cout << std::setw(20) << oss.str();
    }
    std::cout << "\n";

    for (int i = 0; i < m; ++i)
    {
        double H_start = Hn + i * Hel;
        double H_end = Hn + (i + 1) * Hel;
        std::ostringstream oss;
        oss << H_start << "-" << H_end;
        std::cout << std::setw(20) << oss.str();
        for (int j = 0; j < n; ++j)
        {
            std::cout << std::setw(20) << Tr[i][j];
        }
        std::cout << "\n";
    }

    std::cout << "\nМатрица подъема [время в сек]:\n";
    std::cout << std::setw(20) << "V\\H";
    for (int j = 0; j < m; ++j)
    {
        double H_start = Hn + j * Hel;
        double H_end = Hn + (j + 1) * Hel;
        std::ostringstream oss;
        oss << H_start << "-" << H_end;
        std::cout << std::setw(20) << oss.str();
    }
    std::cout << "\n";

    for (int i = 0; i < n; ++i)
    {
        double V_start = (Vn + i * Vel) * 3.6;
        double V_end = (Vn + (i + 1) * Vel) * 3.6;
        std::ostringstream oss;
        oss << V_start << "-" << V_end;
        std::cout << std::setw(20) << oss.str();
        for (int j = 0; j < m; ++j)
        {
            std::cout << std::setw(20) << Tp[i][j];
        }
        std::cout << "\n";
    }

    std::cout << "\nМатрица разгона-подъема [время в сек]:\n";
    std::cout << std::setw(20) << "H\\V";
    for (int j = 0; j < n; ++j)
    {
        double V_start = (Vn + j * Vel) * 3.6;
        double V_end = (Vn + (j + 1) * Vel) * 3.6;
        std::ostringstream oss;
        oss << V_start << "-" << V_end;
        std::cout << std::setw(20) << oss.str();
    }
    std::cout << "\n";

    for (int i = 0; i < n; ++i)
    {
        double H_start = Hn + i * Hel;
        double H_end = Hn + (i + 1) * Hel;
        std::ostringstream oss;
        oss << H_start << "-" << H_end;
        std::cout << std::setw(20) << oss.str();
        for (int j = 0; j < n; ++j)
        {
            std::cout << std::setw(20) << Tpr[i][j];
        }
        std::cout << "\n";
    }

    return 0;
}

// функции из 2 домашнего задания

// плотность воздуха ro(H)
double ro(double H)
{
    return 1.3661924862989573034545085 - 0.0001418065280049439870457 * H + 0.0000000067724452175195049 * pow(H, 2) - 0.0000000000000901906541734 * pow(H, 3);
}

// скорость звука a(H)
double a(double H)
{
    return 342.8873180441046321920151030 - 0.0035032436716794758557481 * H - 0.0000000518023700664785279 * pow(H, 2);
}

// коэффициент подъёмной силы Cy(alpha)
double Cy(double alpha)
{
    return -0.2234285714285730306549738 + 0.1278571428571436963128605 * alpha - 0.0035714285714287135643785 * pow(alpha, 2) + 0.0000000000000000077098821 * pow(alpha, 3);
}

// коэффициент лобового сопротивления Cx(alpha)
double Cx(double alpha)
{
    return 0.2550000000000246513920388 - 0.1408333333333484094218591 * alpha + 0.0275000000000030324354139 * pow(alpha, 2) - 0.0016666666666668671342028 * pow(alpha, 3);
}

// сила тяги двигателей P(M)
double P(double M)
{
    return 44294.8064442556933499872684479 + 47830.2664371878781821578741074 * M - 126839.0779774459660984575748444 * pow(M, 2) + 95497.6562483121379045769572258 * pow(M, 3);
}

// коэффициенты Cy (нужно отдельно аппроксимировать функцию Cy(alpha) как 2 порядка и взять коэффициенты)
double cy0 = -0.1019999999999999795718963;
double cy1 = 0.0849999999999999922284388;

// функции из 3 домашнего задания
double Razgon(double H, double V1, double V2, double mass, double S)
{
    double a_H = a(H);
    double ro_H = ro(H);
    double g = 9.81;

    double V = (V1 + V2) / 2;
    double M = V / a_H;

    double P_M = P(M) * engineCount;

    double numerator = mass * g - P_M / 57.3 - cy0 * ro_H * V * V * S / 2.0;
    double denominator = P_M / 57.3 + cy1 * ro_H * V * V * S / 2.0;
    double alpha_deg = numerator / denominator;
    double alpha = alpha_deg * M_PI / 180.0;

    double Cx_alpha = Cx(alpha_deg);

    double traz = (mass * (V2 - V1) / (P_M * cos(alpha) - Cx_alpha * ro_H * pow(V, 2) / 2 * S * sin(alpha)));
    return traz;
}

double Podyom(double H1, double H2, double V, double mass, double S)
{
    double H = (H1 + H2) / 2;
    double a_H = a(H);
    double g = 9.81;
    double ro_H = ro(H);

    double M = V / a_H;
    double P_M = P(M) * engineCount;

    double numerator = mass * g - P_M / 57.3 - cy0 * ro_H * V * V * S / 2.0;
    double denominator = P_M / 57.3 + cy1 * ro_H * V * V * S / 2.0;
    double alpha_deg = numerator / denominator;

    double Cx_alpha = Cx(alpha_deg);

    double teta = ((P_M - Cx_alpha * ro_H * pow(V, 2) / 2 * S) * 57.3) / (mass * g);
    double tpod = ((57.3 * (H2 - H1)) / (V * teta));
    return tpod;
}

double RazgPod(double H1, double H2, double V1, double V2, double mass, double S)
{
    double H = (H1 + H2) / 2;
    double V = (V1 + V2) / 2;
    double a_H = a(H);
    double g = 9.81;
    double ro_H = ro(H);

    double M = V / a_H;
    double P_M = P(M) * engineCount;

    double numerator = mass * g - P_M / 57.3 - cy0 * ro_H * V * V * S / 2.0;
    double denominator = P_M / 57.3 + cy1 * ro_H * V * V * S / 2.0;
    double alpha_deg = numerator / denominator;
    double alpha = alpha_deg * M_PI / 180.0;

    double Cx_alpha = Cx(alpha_deg);

    double k = (V2 - V1) / (H2 - H1);

    double teta = (P_M * cos(alpha) - Cx_alpha * ro_H * pow(V, 2) / 2 * S) / (mass * (k * V + g));
    double tpodraz = ((1 / (k * teta)) * log(V2 / V1));
    return tpodraz;
}