\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{Расчёт матриц времени \\ \small Программа для расчёта разгона, подъёма и разгон-подъёма}}\\[1.5cm]
        
        \normalsize
        \textbf{Выполнил:} Давыдов Даниил Андреевич\\
        \textbf{Проверил:} Малыхин Дмитрий Андреевич\\[2cm]
        
        \textbf{Дата:} \today
    \end{center}
\end{titlepage}

\tableofcontents
\newpage

\section{Введение}
Цель данной работы заключается в разработке программы для расчёта временных матриц, описывающих три ключевых режима движения летательного аппарата:
\begin{itemize}
    \item Разгон.
    \item Подъём.
    \item Разгон-подъём.
\end{itemize}
Эти матрицы позволяют оценить временные затраты на выполнение различных этапов движения в зависимости от параметров скорости, высоты, массы аппарата и аэродинамических характеристик.

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

\subsection*{1. Время разгона $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*{2. Время подъёма $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*{3. Время разгон-подъёма $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
\]

\section{Реализация программы}
Программа реализована на языке C++ и состоит из следующих основных компонентов:

\begin{enumerate}
    \item \textbf{Функции расчёта аэродинамических параметров:}
    \begin{itemize}
        \item \texttt{double ro(double H);} -- вычисляет плотность воздуха на высоте $H$.
        \item \texttt{double a(double H);} -- вычисляет скорость звука на высоте $H$.
        \item \texttt{double Cx(double alpha);} -- вычисляет коэффициент лобового сопротивления при угле атаки $\alpha$.
        \item \texttt{double Cy(double alpha);} -- вычисляет коэффициент подъёмной силы при угле атаки $\alpha$.
        \item \texttt{double P(double M);} -- вычисляет силу тяги двигателей при числе Маха $M$.
    \end{itemize}
    
    \item \textbf{Функции расчёта времени:}
    \begin{itemize}
        \item \texttt{double Razgon(double H, double V1, double V2, double mass, double S);} -- рассчитывает время разгона.
        \item \texttt{double Podyom(double H1, double H2, double V, double mass, double S);} -- рассчитывает время подъёма.
        \item \texttt{double RazgPod(double H1, double H2, double V1, double V2, double mass, double S);} -- рассчитывает время разгон-подъёма.
    \end{itemize}
    
    \item \textbf{Главная функция:}
    \begin{itemize}
        \item Ввод исходных данных и параметров расчёта.
        \item Расчёт матриц времени для разгона, подъёма и разгон-подъёма.
        \item Вывод результатов в виде таблиц.
    \end{itemize}
\end{enumerate}

Код программы представлен ниже:

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

// плотность воздуха 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
double cy0 = -0.1019999999999999795718963;
double cy1 = 0.0849999999999999922284388;

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;
}
\end{lstlisting}

\section{Результаты}
Результаты расчётов представлены в виде матриц времени для каждого режима движения.

\begin{figure}[H]
    \centering
    \includegraphics[width=1\textwidth]{Screenshot.png}
\end{figure}

\end{document}