\documentclass[a4paper,12pt]{article}
\usepackage[left=2cm,right=2cm,
    top=2cm,bottom=2cm,bindingoffset=0cm]{geometry}
\usepackage{amsmath, amssymb}
\usepackage{graphicx}
\usepackage{hyperref}
\usepackage[utf8]{inputenc}
\usepackage[russian]{babel}
\usepackage{cite}
\usepackage{listings}
\lstset{extendedchars=\true, inputencoding=utf8}
\usepackage{color}
\usepackage{float}
\usepackage{caption}

\definecolor{mygreen}{rgb}{0,0.6,0}
\definecolor{mygray}{rgb}{0.5,0.5,0.5}
\definecolor{mymauve}{rgb}{0.58,0,0.82}

\lstset{
    language=C++,
    aboveskip=3mm,
    belowskip=3mm,
    showstringspaces=false,
    columns=flexible,
    basicstyle={\small\ttfamily},
    numbers=left,
    numberstyle=\tiny\color{mygray},
    keywordstyle=\color{blue},
    commentstyle=\color{mygreen},
    stringstyle=\color{mymauve},
    breaklines=true,
    breakatwhitespace=true,
    tabsize=4
}

\begin{document}

%----------------- ТИТУЛЬНЫЙ ЛИСТ -----------------%
\begin{titlepage}
    \begin{center}
        \vspace*{2cm}
        \Large{\textbf{Объединённый отчёт: \\
        1) Полиномиальная аппроксимация \\
        2) Параметры самолёта \\
        3) Расчёт матриц времени \\
        4) Оптимальный маршрут}}\\[1.5cm]
        
        \normalsize
        \textbf{Выполнил:} \\
        \textit{Давыдов Даниил Андреевич} \\
        \vspace{0.4cm}
        \textbf{Проверил:} \\
        \textit{Малыхин Дмитрий Андреевич}\\[2cm]
        
        \textbf{Дата:} \today
    \end{center}
\end{titlepage}

\tableofcontents
\newpage


%============================================================%
%                         ВВЕДЕНИЕ                          %
%============================================================%
\section{Введение}

В данном отчёте собраны четыре тематических блока, каждый из которых представляет собой законченный небольшой проект (или часть проекта) на языке C++. Они оформлены в виде отдельных подпрограмм, сопровождаемых описательной частью, теоретическими сведениями, листингами кода, а также (где это необходимо) иллюстрациями с графиками.

Объединяемые блоки:

\begin{enumerate}
    \item \textbf{Полиномиальная аппроксимация (метод наименьших квадратов)}.
    \item \textbf{Получение параметров самолёта (расчёт плотности, скорости звука, аэродинамических коэффициентов и т.д.)}.
    \item \textbf{Формирование матриц времени разгона, подъёма и разгон-подъёма}.
    \item \textbf{Поиск оптимального маршрута (динамическое программирование по матрицам времени)}.
\end{enumerate}

В последующих разделах представлен полный код решений, сопровождаемый комментариями по работе алгоритмов и ключевым результатам.

\newpage


%============================================================%
%           1. ПОЛИНОМИАЛЬНАЯ АППРОКСИМАЦИЯ                  %
%============================================================%
\section{Полиномиальная аппроксимация}

\subsection{Постановка задачи}

Задача — реализовать метод наименьших квадратов для полиномиальной аппроксимации экспериментальных данных. Пользователь может задать требуемую степень полинома, а программа должна вычислить искомые коэффициенты, решая систему нормальных уравнений методом Гаусса.

\subsection{Основные формулы}

Метод наименьших квадратов (МНК) предполагает минимизацию функционала:

\[
S = \sum_{i=1}^{N} [y_i - P(x_i)]^2.
\]

Для полинома степени $n$:

\[
P(x) = a_0 + a_1 x + a_2 x^2 + \dots + a_n x^n,
\]

минимизация суммы квадратов приводит к системе нормальных уравнений, которая решается методом Гаусса.

\subsection{Листинг кода}

Ниже приведён листинг программы на C++.  
Обратите внимание, что здесь степень полинома передаётся через аргумент командной строки (\texttt{argv[1]}), а данные считываются из стандартного ввода:

\begin{lstlisting}
#include <iostream>
#include <vector>
#include <cmath>
#include <limits>
#include <iomanip>
#include <string>
#include <sstream>
#include <cstdlib>

std::vector<double> polynomialFit(const std::vector<double>& x, const std::vector<double>& y, int degree) {
    int N = (int)x.size();
    int n = degree;

    std::vector<double> coeffs(n + 1, 0.0);

    std::vector<std::vector<double>> A(n + 1, std::vector<double>(n + 2, 0.0));

    // Заполнение матрицы нормальных уравнений
    for (int i = 0; i <= n; ++i) {
        for (int j = 0; j <= n; ++j) {
            double sum = 0.0;
            for (int k = 0; k < N; ++k) {
                sum += std::pow(x[k], i + j);
            }
            A[i][j] = sum;
        }
        double sum = 0.0;
        for (int k = 0; k < N; ++k) {
            sum += y[k] * std::pow(x[k], i);
        }
        A[i][n + 1] = sum;
    }

    // Прямой ход метода Гаусса
    for (int i = 0; i <= n; ++i) {
        double maxElement = std::fabs(A[i][i]);
        int maxRow = i;
        for (int k = i + 1; k <= n; ++k) {
            if (std::fabs(A[k][i]) > maxElement) {
                maxElement = std::fabs(A[k][i]);
                maxRow = k;
            }
        }

        for (int k = i; k <= n + 1; ++k) {
            std::swap(A[maxRow][k], A[i][k]);
        }

        for (int k = i + 1; k <= n; ++k) {
            double factor = -A[k][i] / A[i][i];
            for (int j = i; j <= n + 1; ++j) {
                if (i == j) {
                    A[k][j] = 0;
                } else {
                    A[k][j] += factor * A[i][j];
                }
            }
        }
    }

    // Обратный ход
    for (int i = n; i >= 0; --i) {
        coeffs[i] = A[i][n + 1] / A[i][i];
        for (int k = i - 1; k >= 0; --k) {
            A[k][n + 1] -= A[k][i] * coeffs[i];
        }
    }

    return coeffs;
}

int main(int argc, char* argv[]) {
    if (argc < 2) {
        std::cerr << "Ошибка: не указана степень полинома. Использование: " 
                  << argv[0] << " degree" << std::endl;
        return 1;
    }

    int degree = std::atoi(argv[1]);
    if (degree < 0) {
        std::cerr << "Ошибка: степень полинома должна быть неотрицательным числом." 
                  << std::endl;
        return 1;
    }

    std::vector<double> x;
    std::vector<double> y;

    std::string line;
    while (std::getline(std::cin, line)) {
        std::stringstream ss(line);
        double xi, yi;
        if (ss >> xi >> yi) {
            x.push_back(xi);
            y.push_back(yi);
        }
    }

    if (x.size() != y.size() || x.empty()) {
        std::cerr << "Ошибка: некорректные входные данные." << std::endl;
        return 1;
    }

    std::vector<double> coeffs = polynomialFit(x, y, degree);

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "DEGREE " << degree << std::endl;
    std::cout << "COEFFICIENTS ";
    for (size_t i = 0; i < coeffs.size(); ++i) {
        std::cout << coeffs[i];
        if (i != coeffs.size() - 1)
            std::cout << ", ";
    }
    std::cout << std::endl;

    return 0;
}
\end{lstlisting}

\subsection{Тестовый пример и визуализация}

Для тестовых данных:

\begin{verbatim}
-3 9
-1 3
0 0
1 3
3 9
\end{verbatim}

(пример параболической зависимости), программа, запущенная с аргументом \texttt{2} (т.е. аппроксимация квадратичным полиномом), даёт результат:

\begin{verbatim}
a[0] = 1.3714285714285723294381114
a[1] = 0.0000000000000000000000000
a[2] = 0.8571428571428570952761561
\end{verbatim}

Ниже приведён пример графика аппроксимации данных полиномом второй степени:

\begin{figure}[H]
    \centering
    \includegraphics[width=0.6\textwidth]{../approximation/approximation.png}
    \caption{Пример аппроксимации параболы полиномом 2 степени.}
\end{figure}

\newpage


%============================================================%
%                2. ПАРАМЕТРЫ САМОЛЁТА                      %
%============================================================%
\section{Параметры самолёта}

\subsection{Описание задачи}

В данном разделе приводится код, который вычисляет ординаты аппроксимирующих функций: плотность воздуха \(\rho(H)\), скорость звука \(a(H)\), коэффициент подъёмной силы \(Cy(\alpha)\), лобовое сопротивление \(Cx(\alpha)\) и силу тяги \(P(M)\). Программа также формирует сетку по высоте, скорости, углу атаки, числу Маха и выводит результаты в табличном виде.

\subsection{Листинг программы}

\begin{lstlisting}
#include <iostream>
#include <cmath>
#include <vector>
#include <iomanip>

using namespace std;

// плотность воздуха 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 19.0500000000000397903932026 
           - 12.7708333333333623471617102 * alpha 
           + 3.2104166666666742457891814 * pow(alpha, 2) 
           - 0.3541666666666675733488034 * pow(alpha, 3) 
           + 0.0145833333333333717285463 * pow(alpha, 4);
}

// сила тяги двигателей P(M)
double P(double M)
{
    return 47589.0511051998255425132811069 
           - 26183.6352466731186723336577415 * M 
           - 79976.5836117977305548265576363 * pow(M, 2) 
           + 62418.6378453120341873727738857 * pow(M, 3);
}

int get_k()
{
    int k;
    cout << "Введите количество разбиений (от 1 до 99): ";
    cin >> k;

    if (k <= 0 || k >= 100)
    {
        cout << "Ошибка: количество разбиений должно быть от 1 до 99.\n";
        return get_k();
    }
    return k;
}

int main()
{
    int k = get_k();

    vector<double> M(k + 1);
    vector<double> H(k + 1);
    vector<double> V(k + 1);
    vector<double> alpha(k + 1);

    // Задаём диапазоны
    V[0] = 310 / 3.6;   // м/с
    V[k] = 700 / 3.6;   // м/с

    H[0] = 300;
    H[k] = 6000;

    cout << "Введите начальное число Маха: ";
    cin >> M[0];
    cout << "Введите конечное число Маха: ";
    cin >> M[k];

    // Равномерное разбиение
    for (int i = 0; i < k; ++i)
    {
        H[i + 1] = H[i] + (H[k] - H[0]) / k;
        V[i + 1] = V[i] + (V[k] - V[0]) / k;
        M[i + 1] = M[i] + (M[k] - M[0]) / k;
    }

    cout << "Введите начальный угол атаки (градусы): ";
    cin >> alpha[0];
    for (int i = 0; i < k; ++i)
    {
        alpha[i + 1] = alpha[i] + 1.0;
    }

    // Выводим таблицу
    cout << fixed;
    cout << "\nРезультаты:\n";
    cout << left << setw(7) << "№"
         << right << setw(10) << "H (м)"
         << right << setw(15) << "V (м/с)"
         << right << setw(12) << "alpha (°)"
         << right << setw(12) << "M"
         << right << setw(12) << "a (м/с)"
         << right << setw(12) << "ro (кг/м^3)"
         << right << setw(12) << "Cy"
         << right << setw(10) << "Cx"
         << right << setw(15) << "P (кгс)"
         << endl;

    cout << string(107, '-') << endl;

    for (int i = 0; i <= k; ++i)
    {
        cout << left << setw(5) << (i + 1);
        cout << right << setw(10) << setprecision(0) << H[i];
        cout << right << setw(10) << setprecision(2) << V[i];
        cout << right << setw(12) << setprecision(2) << alpha[i];
        cout << right << setw(10) << setprecision(4) << M[i];
        cout << right << setw(10) << setprecision(2) << a(H[i]);
        cout << right << setw(15) << setprecision(6) << ro(H[i]);
        cout << right << setw(10) << setprecision(4) << Cy(alpha[i]);
        cout << right << setw(10) << setprecision(4) << Cx(alpha[i]);
        cout << right << setw(15) << setprecision(2) << P(M[i]);
        cout << endl;
    }

    return 0;
}
\end{lstlisting}

\subsection{Примеры графиков}

Ниже приведены примеры графиков, полученных на основе рассмотренных функций (плотность, скорость звука и т.д.):

\begin{figure}[H]
    \centering
    \includegraphics[width=0.5\textwidth]{../TY-134-graphs/images/1.png}
    \caption{График зависимости плотности воздуха $\rho(H)$ от высоты $H$.}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[width=0.5\textwidth]{../TY-134-graphs/images/2.png}
    \caption{График зависимости скорости звука $a(H)$ от высоты $H$.}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[width=0.5\textwidth]{../TY-134-graphs/images/3.png}
    \caption{График зависимости $C_y(\alpha)$ от угла атаки $\alpha$.}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[width=0.5\textwidth]{../TY-134-graphs/images/4.png}
    \caption{График зависимости $C_x(\alpha)$ от угла атаки $\alpha$.}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[width=0.5\textwidth]{../TY-134-graphs/images/5.png}
    \caption{График зависимости силы тяги $P(M)$ от числа Маха $M$.}
\end{figure}

\newpage


%============================================================%
%                3. МАТРИЦЫ ВРЕМЕНИ                         %
%============================================================%
\section{Формирование матриц времени}

\subsection{Постановка задачи}

Данный раздел посвящён расчёту временных матриц, описывающих три ключевых режима движения летательного аппарата:
\begin{itemize}
    \item разгон;
    \item подъём;
    \item разгон-подъём.
\end{itemize}

Формируются соответствующие матрицы (размера \((n+1)\times n\), \(n\times (n+1)\) или \(n\times n\)), где каждый элемент — это время выполнения соответствующего элементарного шага.

\section{Исходные данные}
Для расчётов были заданы следующие исходные параметры:
\begin{itemize}
    \item Начальная скорость: $V_n = 310 \, \text{км/ч}$.
    \item Конечная скорость: $V_k = 700 \, \text{км/ч}$.
    \item Начальная высота: $H_n = 300 \, \text{м}$.
    \item Конечная высота: $H_k = 6000 \, \text{м}$.
    \item Масса летательного аппарата: $m = 47\,000 \, \text{кг}$.
    \item Площадь крыла: $S = 127 \, \text{м}^2$.
    \item Угол атаки: $\alpha = 4^\circ$.
    \item Ускорение свободного падения: $g = 9{,}81 \, \text{м/с}^2$.
\end{itemize}
Данные параметры соответствуют характеристикам типичного коммерческого самолёта среднего класса, что позволяет провести анализ в условиях, приближенных к реальным эксплуатационным ситуациям.

\section{Методика расчётов}
Для вычисления временных затрат на разгон, подъём и разгон-подъём используются следующие формулы:

Для проекций сил на оси скоростной системы координат получим следующие два уравнения:

\begin{equation}
    m \frac{dV}{dt} = P \cos(\alpha + \varphi_p) - X - m g \sin\theta,
\end{equation}

\begin{equation}
    m V \frac{d\theta}{dt} = P \sin(\alpha + \varphi_p) + Y - m g \cos\theta.
\end{equation}

Где:
\begin{itemize}
    \item $m$ — масса тела;
    \item $V$ — скорость;
    \item $\theta$ — угол наклона траектории;
    \item $P$ — сила, приложенная под углом $(\alpha + \varphi_p)$;
    \item $X$ — сила бокового сопротивления;
    \item $Y$ — подъемная сила;
    \item $g$ — ускорение свободного падения.
\end{itemize}

Кинематические уравнения, определяющие проекции скорости на оси нормальной земной системы координат для поступательного движения, имеют вид:

\begin{equation}
    \frac{dH}{dt} = V \sin\theta,
\end{equation}

\begin{equation}
    \frac{dL}{dt} = V \cos\theta.
\end{equation}

Здесь $H$ и $L$ — проекции перемещения на вертикальную и горизонтальную оси соответственно.

\subsection{Время разгона $t_{\text{раз}}$}
\[
t_{\text{раз}} = \frac{m \cdot (V_2 - V_1)}{P(M) \cos\alpha - C_x(\alpha) \cdot \dfrac{\rho(H) \cdot V^2}{2} \cdot S}
\]
где:
\begin{itemize}
    \item $V_1$, $V_2$ -- начальная и конечная скорости (м/с).
    \item $V = \dfrac{V_1 + V_2}{2}$ -- средняя скорость (м/с).
    \item $m$ -- масса аппарата (кг).
    \item $P(M)$ -- сила тяги двигателей, зависящая от числа Маха $M$ (Н).
    \item $M = \dfrac{V}{a(H)}$ -- число Маха.
    \item $a(H)$ -- скорость звука на высоте $H$ (м/с).
    \item $\alpha$ -- угол атаки (рад).
    \item $C_x(\alpha)$ -- коэффициент лобового сопротивления.
    \item $\rho(H)$ -- плотность воздуха на высоте $H$ (кг/м$^3$).
    \item $S$ -- площадь крыла (м$^2$).
\end{itemize}

\subsection{Время подъёма $t_{\text{под}}$}
\[
t_{\text{под}} = \frac{57{,}3 \cdot (H_2 - H_1)}{V \cdot \theta}
\]
где:
\begin{itemize}
    \item $H_1$, $H_2$ -- начальная и конечная высоты (м).
    \item $V$ -- скорость (м/с).
    \item $\theta$ -- угол набора высоты (град), определяется как:
    \[
    \theta = \dfrac{(P(M) - C_x(\alpha) \cdot \dfrac{\rho(H) \cdot V^2}{2} \cdot S) \cdot 57{,}3}{m \cdot g}
    \]
    \item $g$ -- ускорение свободного падения (м/с$^2$).
\end{itemize}

\subsection{Время разгон-подъёма $t_{\text{раз-под}}$}
\[
t_{\text{раз-под}} = \frac{1}{k \cdot \theta} \ln{\left( \frac{V_2}{V_1} \right)}
\]
где:
\begin{itemize}
    \item $k = \dfrac{V_2 - V_1}{H_2 - H_1}$ -- коэффициент изменения скорости относительно высоты.
    \item $\theta$ -- угол набора высоты (град), аналогичный предыдущей формуле.
\end{itemize}

\subsection{Дополнительные зависимости}
Для точности расчётов используются следующие эмпирические зависимости:

\textbf{Плотность воздуха $\rho(H)$}:
\[
\rho(H) = 1{,}36619 - 0{,}000141807 \cdot H + 6{,}77245 \times 10^{-9} \cdot H^2 - 9{,}01907 \times 10^{-14} \cdot H^3
\]

\textbf{Скорость звука $a(H)$}:
\[
a(H) = 342{,}887 - 0{,}00350324 \cdot H - 5{,}18024 \times 10^{-8} \cdot H^2
\]

\textbf{Коэффициент подъёмной силы $C_y(\alpha)$}:
\[
C_y(\alpha) = -0{,}22343 + 0{,}12786 \cdot \alpha - 0{,}00357 \cdot \alpha^2
\]

\textbf{Коэффициент лобового сопротивления $C_x(\alpha)$}:
\[
C_x(\alpha) = 19{,}05 - 12{,}7708 \cdot \alpha + 3{,}2104 \cdot \alpha^2 - 0{,}35417 \cdot \alpha^3 + 0{,}01458 \cdot \alpha^4
\]

\textbf{Сила тяги двигателей $P(M)$}:
\[
P(M) = 44\,294{,}81 + 47\,830{,}27 \cdot M - 126\,839{,}08 \cdot M^2 + 95\,497{,}66 \cdot M^3
\]

\subsection{Листинг программы}

\begin{lstlisting}
#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>
#include <sstream>

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;
        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;
}

// ------------------- Реализация функций ------------------- //

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

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

// Коэффициенты Cy (упрощённый пример)
double cy0 = -0.1019999999999999795718963;
double cy1 = 0.0849999999999999922284388;

// Коэффициент подъёмной силы (пример для расчётов)
double Cy(double alpha) {
    // Можно использовать более точную аппроксимацию
    return cy0 + cy1 * 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);
}

// Время разгона
double Razgon(double H, double V1, double V2, double mass, double S_area) {
    double a_H = a(H);
    double ro_H = ro(H);
    double g = 9.81;

    double V = (V1 + V2) / 2.0;
    double M = V / a_H;
    double P_M = P(M) * engineCount;

    // Пример вычисления "эффективного" alpha через баланс сил:
    double numerator = mass * g - P_M / 57.3 - cy0 * ro_H * V * V * S_area / 2.0;
    double denominator = P_M / 57.3 + cy1 * ro_H * V * V * S_area / 2.0;
    if (std::fabs(denominator) < 1e-9) return 1e9;

    double alpha_deg = numerator / denominator;
    double alpha = alpha_deg * M_PI / 180.0;

    double Cx_alpha = Cx(alpha_deg);
    double cos_alpha = cos(alpha);
    double sin_alpha = sin(alpha);

    double denominator_time = P_M * cos_alpha 
                              - Cx_alpha * ro_H * pow(V, 2) / 2.0 * S_area * sin_alpha;
    if (std::fabs(denominator_time) < 1e-9) return 1e9;

    double traz = mass * (V2 - V1) / denominator_time;
    return traz;
}

// Время подъёма
double Podyom(double H1, double H2, double V, double mass, double S_area) {
    double H = (H1 + H2) / 2.0;
    double a_H = a(H);
    double ro_H = ro(H);
    double g = 9.81;

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

    // Оценка alpha
    double numerator = mass * g - P_M / 57.3 - cy0 * ro_H * V * V * S_area / 2.0;
    double denominator = P_M / 57.3 + cy1 * ro_H * V * V * S_area / 2.0;
    double alpha_deg = (std::fabs(denominator) < 1e-9) ? 0.0 : numerator / denominator;

    double Cx_alpha = Cx(alpha_deg);
    double teta = ((P_M - Cx_alpha * ro_H * pow(V, 2) / 2 * S_area) * 57.3) 
                  / (mass * g);
    if (std::fabs(teta) < 1e-9) return 1e9;

    double dp = H2 - H1;
    double tpod = (57.3 * dp) / (V * teta);
    return tpod;
}

// Время разгон-подъема
double RazgPod(double H1, double H2, double V1, double V2, double mass, double S_area) {
    double H = (H1 + H2) / 2.0;
    double V = (V1 + V2) / 2.0;
    double a_H = a(H);
    double ro_H = ro(H);
    double g = 9.81;

    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_area / 2.0;
    double denominator = P_M / 57.3 + cy1 * ro_H * V * V * S_area / 2.0;
    if (std::fabs(denominator) < 1e-9) return 1e9;

    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_area) 
                  / (mass * (k * V + g));
    if (std::fabs(teta) < 1e-9) return 1e9;

    double tpodraz = (1.0 / (k * teta)) * log(V2 / V1);
    return tpodraz;
}
\end{lstlisting}

\subsection{Пример результата}

На рисунке ниже приведён пример консольного вывода с частично заполненными таблицами для заданных пользователем \(\textit{n}\).

\begin{figure}[H]
    \centering
    \includegraphics[width=1\textwidth]{../time/Screenshot.png}
    \caption{Пример фрагмента вывода программы (матрицы разгона, подъёма и разгон-подъёма).}
\end{figure}

\newpage


%============================================================%
%               4. ОПТИМАЛЬНЫЙ МАРШРУТ                      %
%============================================================%
\section{Поиск оптимального маршрута}

\subsection{Идея алгоритма}

Чтобы найти маршрут (последовательность действий) по высоте и скорости, который позволит минимизировать общее время, удобно воспользоваться методом динамического программирования:
\begin{itemize}
    \item Каждая точка на сетке \((i, j)\) соответствует некоторой высоте \(H_i\) и скорости \(V_j\).
    \item Из точки \((i, j)\) можно перейти в точку \((i, j+1)\) (разгон), \((i+1, j)\) (подъём) или \((i+1, j+1)\) (разгон-подъём).
    \item Стоимость перехода берётся из заранее рассчитанных матриц \(Tr\), \(Tp\) и \(Tpr\).
    \item Находим сумму времени (стоимость) от текущей точки до финишной (\((m-1, n-1)\)).
\end{itemize}
В результате получаем матрицу \(S\) минимальных времён и матрицу \(\textit{Path}\), по которой можно восстановить конкретную последовательность шагов.

\subsection{Код программы}

Ниже приведён пример реализации (кусочно) главной части программы, которая после расчёта матриц времени строит матрицу \(S\) с помощью динамического программирования, а затем восстанавливает путь:

\begin{lstlisting}
#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include <iomanip>
#include <sstream>
#include <string>

// ... [здесь предполагается подключение функций Tr, Tp, Tpr, ro, a, Cx, Cy, P, etc.] ...

int main()
{
    // Предполагается, что мы уже имеем:
    //  - Матрицы Tr, Tp, Tpr
    //  - n, m, шаги Hel, Vel, диапазоны Hn..Hk, Vn..Vk

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

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

    // Установка целевого узла (верхний правый угол) в 0
    S[m - 1][n - 1] = 0.0;
    Path[m - 1][n - 1] = -1; // Целевой узел

    // Заполнение правого столбца (j = n - 1)
    for (int i = m - 2; i >= 0; --i)
    {
        S[i][n - 1] = S[i + 1][n - 1] + Tp[n - 1][i];
        Path[i][n - 1] = 1; // Подъем
    }

    // Заполнение нижней строки (i = m - 1)
    for (int j = n - 2; j >= 0; --j)
    {
        S[m - 1][j] = S[m - 1][j + 1] + Tr[m -1][j];
        Path[m - 1][j] = 0; // Разгон
    }

    // Заполнение остальной части матрицы S
    for (int i = m - 2; i >= 0; --i)
    {
        for (int j = n - 2; j >= 0; --j)
        {
            // 0) Разгон
            double option0 = S[i][j + 1] + Tr[i][j];

            // 1) Подъем
            double option1 = S[i + 1][j] + Tp[j][i];

            // 2) Разгон-подъем
            double option2 = S[i + 1][j + 1] + Tpr[j][i];

            // Выбираем минимум
            double min_time = option0;
            int path_choice = 0; // Разгон

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

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

            // Обновляем
            S[i][j] = min_time;
            Path[i][j] = path_choice;
        }
    }

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

    // Запись итоговой матрицы S в файл
    outfile << "Матрица минимальных времён (S):\n";
    for (int i = 0; i < m; ++i)
    {
        for (int j = 0; j < n; ++j)
        {
            outfile << std::setw(8) << S[i][j] << " ";
        }
        outfile << "\n";
    }

    // Восстановление самого пути
    outfile << "\nОптимальный маршрут:\n";
    double total_time = S[0][0];
    outfile << "Общее время: " << total_time << " сек\n";

    int ci = 0;
    int cj = 0;
    outfile << "Шаги (H, V): \n";

    while (Path[ci][cj] != -1)
    {
        outfile << "(" << (Hn + ci*Hel) << ", " << (Vn + cj*Vel)*3.6 
                << " км/ч) -> ";

        int action = Path[ci][cj];
        if (action == 0) {
            // Разгон
            cj++;
        } else if (action == 1) {
            // Подъем
            ci++;
        } else {
            // Разгон-подъем
            ci++;
            cj++;
        }
    }
    // Финальная точка
    outfile << "(" << (Hn + (m-1)*Hel) << ", " 
            << (Vn + (n-1)*Vel)*3.6 << " км/ч)\n";

    outfile.close();
    std::cout << "Данные записаны в файл output.txt\n";

    return 0;
}
\end{lstlisting}

\subsection{Вывод программы}
\begin{figure}[H]
    \centering
    \includegraphics[width=1\textwidth]{../best-path/console.png}
    \caption{Полученные матрицы времени, пути}
    \label{fig:route_plot}
\end{figure}

\subsection{Пример полученного графика}
\begin{figure}[H]
    \centering
    \includegraphics[width=1\textwidth]{../best-path/route_plot.png}
    \caption{График оптимального маршрута полёта}
    \label{fig:route_plot}
\end{figure}

\subsection{Результаты и визуализация}

После выполнения кода формируется файл \texttt{output.txt}, в котором:
\begin{itemize}
    \item Приведена итоговая матрица \(S\) (минимальные времена для каждой ячейки сетки).
    \item Записано общее время движения по оптимальному маршруту.
    \item Указана последовательность шагов (тип действия: разгон, подъём, разгон-подъём), а также координаты \((H, V)\).
\end{itemize}

Для наглядности можно использовать дополнительный Python-скрипт, который считывает координаты из \texttt{output.txt} и строит график в осях \((H, V)\), показывая, где выполнялся разгон, где подъём, а где их комбинация.

\newpage

%============================================================%
%                         ЗАКЛЮЧЕНИЕ                        %
%============================================================%
\section{Заключение}

В данной работе последовательно рассмотрены четыре взаимосвязанные подпрограммы (проекты):
\begin{enumerate}
    \item Программа для полиномиальной аппроксимации (метод наименьших квадратов).
    \item Код для вычисления основных аэродинамических параметров самолёта (плотность воздуха, скорость звука и т.д.) и вывода таблиц результатов.
    \item Формирование матриц времени разгона, подъёма и разгон-подъёма.
    \item Определение оптимального маршрута (минимизация времени) с помощью динамического программирования.
\end{enumerate}

Каждый из этих модулей может использоваться автономно, но вместе они позволяют не только аппроксимировать входные зависимости (например, аэродинамические характеристики), но и проводить более комплексный расчёт траектории самолёта, исходя из сформированных матриц времени. Оптимальный путь по времени достигается пошаговым выбором между режимами разгона, подъёма или их комбинации.

\newpage

%============================================================%
%                         ИСТОЧНИКИ                         %
%============================================================%
\begin{thebibliography}{9}
\bibitem{gauss} 
  \textit{Метод Гаусса для чайников.} 
  [Электронный ресурс]. URL: \url{http://www.mathprofi.ru/metod_gaussa_dlya_chainikov.html} (дата обращения: 10.12.2024).

\bibitem{mnk} 
  \textit{Метод наименьших квадратов.} 
  [Электронный ресурс]. URL: \url{http://www.mathprofi.ru/metod_naimenshih_kvadratov.html} (дата обращения: 10.12.2024).

\bibitem{habr_article}
  Н. Иванов. \textit{Аппроксимация данных и МНК.} --- 
  [Электронный ресурс]. URL: \url{https://habr.com/ru/articles/672540/} (дата обращения: 10.12.2024).

\bibitem{lsq_notes}
  \textit{Method of Least Squares.} 
  [Электронный ресурс]. Williams College. URL: \url{https://web.williams.edu/Mathematics/sjmiller/public_html/BrownClasses/54/handouts/MethodLeastSquares.pdf} (дата обращения: 10.12.2024).

\bibitem{dp_path} 
  Cormen T. H., Leiserson C. E., Rivest R. L., Stein C. 
  \textit{Introduction to Algorithms.} --- The MIT Press, 2009.

\end{thebibliography}

\end{document}