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

// для отладки
void printTimeMatrices(const std::vector<std::vector<double>> &Tr,
                       const std::vector<std::vector<double>> &Tp,
                       const std::vector<std::vector<double>> &Tpr,
                       double Hel, double Vel, double Hn, double Vn, int m, int n)
{
    // вывод значений
    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";
    }
}
void printDebugInfo(const std::vector<std::vector<double>> &S, 
                   const std::vector<std::vector<int>> &Path,
                   double Hel, double Vel, double Hn, double Vn, int m, int n)
{
    std::vector<std::string> H_labels;
    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;
        H_labels.push_back(oss.str());
    }

    std::vector<std::string> V_labels;
    for (int j = 0; j < m; ++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;
        V_labels.push_back(oss.str());
    }

    // Вывод матрицы S
    std::cout << "\nМатрица минимальных времён (S) [время в сек]:\n";
    std::cout << std::setw(20) << "H\\V";
    for (const auto& V_label : V_labels)
    {
        std::cout << std::setw(20) << V_label;
    }
    std::cout << "\n";
    for (size_t i = 0; i < H_labels.size(); ++i)
    {
        std::cout << std::setw(20) << H_labels[i];
        for (size_t j = 0; j < V_labels.size(); ++j)
        {
            std::cout << std::setw(20) << S[i][j];
        }
        std::cout << "\n";
    }
    
    // Вывод матрицы Path
    std::cout << "\nМатрица маршрутов (Path):\n";
    for (int i = 0; i < m; ++i)
    {
        for (int j = 0; j < m; ++j)
        {
            std::cout << std::setw(3) << Path[i][j] << " ";
        }
        std::cout << "\n";
    }
}

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

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_area = 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_area);
            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_area);
            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_area);
            Tpr[i][j] = tpodraz;
        }
    }

    // Открываем файл для записи
    std::ofstream outfile("output.txt");
    if (!outfile.is_open()) {
        std::cerr << "Не удалось открыть файл для записи.\n";
        return -1;
    }

    // Устанавливаем формат вывода
    outfile << std::fixed << std::setprecision(2);

    // Создание и инициализация матрицы S (минимальные времена)
    std::vector<std::vector<double>> S(m, std::vector<double>(m));

    // Создание матрицы Path для восстановления пути
    // 0 - Разгон, 1 - Подъем, 2 - Разгон-Подъем
    std::vector<std::vector<int>> Path(m, std::vector<int>(m, -1));

    for (int i = m - 2; i >= 0; --i)
    {
        S[i][m - 1] = S[i + 1][m - 1] + Tp[i][m - 1];
        Path[i][m - 1] = 1; // Подъем
    }

    for (int j = m - 2; j >= 0; --j)
    {
        S[m - 1][j] = S[m - 1][j + 1] + Tr[m - 1][j];
        Path[m - 1][j] = 0; // Разгон
    }

    printDebugInfo(S, Path, Hel, Vel, Hn, Vn, m, n);
    printTimeMatrices(Tr, Tp, Tpr, Hel, Vel, Hn, Vn, m, n);

    // Заполнение оставшейся части матрицы S
    for (int i = m - 2; i >= 0; --i)
    {
        for (int j = m - 2; j >= 0; --j)
        {
            // std::cout << "Обработка ячейки (" << i << ", " << j << ")\n";
            double option0 = S[i][j + 1] + Tr[i][j];
            // std::cout << "Разгон: " << S[j][i + 1] << "; "<< Tr[i][j] << std::endl;

            double option1 = S[i + 1][j] + Tp[i][j];
            // std::cout << "Подъем: " << S[i + 1][j] << "; "<< Tp[i][j] << std::endl;

            double option2 = S[i][j] + Tpr[j][i];
            // std::cout << "Разгон-подъем: " << S[i][j] << "; "<< Tpr[j][i] << std::endl;

            // Находим минимальное из трех вариантов
            double min_time = option0;
            int path_choice = 0; // Разгон

            if (option1 > 0 && option1 < min_time)
            {
                min_time = option1;
                path_choice = 1; // Подъем
            }

            if (option2 > 0 && option2 < min_time)
            {
                min_time = option2;
                path_choice = 2; // Разгон-Подъем
            }

            S[i][j] = min_time;
            Path[i][j] = path_choice;
        }
    }

    printTimeMatrices(Tr, Tp, Tpr, Hel, Vel, Hn, Vn, m, n);
    printDebugInfo(S, Path, Hel, Vel, Hn, Vn, m, n);

    // Восстановление оптимального пути и запись в файл
    outfile << "\nОптимальный маршрут по времени:\n";
    int current_i = 0; // Начинаем с начального узла (низкая высота, низкая скорость)
    int current_j = 0;
    double total_time = S[current_i][current_j];
    outfile << "Общее время: " << total_time << " сек\n";
    outfile << "Последовательность действий (начиная с H=" << Hn << ", V=" << Vn * 3.6 << " км/ч):\n";

    // Создадим секцию с координатами для удобства чтения в Python
    outfile << "\nКоординаты оптимального пути (H_start, V_start, H_end, V_end, Action):\n";
    outfile << "H_start,V_start,H_end,V_end,Action\n";

    while (Path[current_i][current_j] != -1)
    {
        double H_current = Hn + current_i * Hel;
        double V_current = Vn + current_j * Vel;
        double H_next, V_next;
        std::string action_str;

        switch (Path[current_i][current_j])
        {
        case 0:
            {
                H_next = H_current;
                V_next = V_current + Vel;
                action_str = "Разгон";
                outfile << H_current << "," << V_current * 3.6 << "," << H_next << "," << V_next * 3.6 << "," << action_str << "\n";
                current_j += 1;
            }
            break;
        case 1:
            {
                H_next = H_current + Hel;
                V_next = V_current;
                action_str = "Подъем";
                outfile << H_current << "," << V_current * 3.6 << "," << H_next << "," << V_next * 3.6 << "," << action_str << "\n";
                current_i += 1;
            }
            break;
        case 2:
            {
                H_next = H_current + Hel;
                V_next = V_current + Vel;
                action_str = "Разгон-Подъем";
                outfile << H_current << "," << V_current * 3.6 << "," << H_next << "," << V_next * 3.6 << "," << action_str << "\n";
                current_i += 1;
                current_j += 1;
            }
            break;
        }
        outfile << H_current << "," << V_current * 3.6 << "," << H_next << "," << V_next * 3.6 << "," << action_str << "\n";
        // std::cout << "\n" << current_i <<  "; " << current_j << "; " << action << "\n";
    }

    outfile.close();

    std::cout << "Данные успешно записаны в файл 'output.txt'.\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;
}