%% Начальные данные

M_inf_0 = 2;                                % Число Маха, [б/р]
H_0 = 3.5;                                  % Добавка к геопотенциальной высоте, [км]
beta_k_0 = 7.5;                             % Добавка к углу полураствора конуса, [deg]
x_cm = 0.65;                                % положение ц.м., [м]
d_m = 1.4;                                  % Диаметр миделя, [м]
N = 4;                                      % Номер варианта
P = 1;                                      % Параметр
k = 1.4;                                    % Отношение удельных теплоемкомстей воздуха [б/р]
R = 287.05287;                              % Газовая постоянная, Дж/(кг*К)

%% Геопотенциальная высота и число Маха

H = (H_0 + 0.2 * N) * 10^3;                 % Геопотенциальная высота, [м]
M_inf = M_inf_0 + (P - 1);                  % Число Маха, [б/р]

%% Геометрия конуса

beta_k = beta_k_0 + 2.5 * (N - 2 * P);      % Угол полураствора конуса, [deg]
l_k = d_m / (2 * tand(beta_k));             % Длина острого конуса, [м]
% x_cm = x_cm_s * l_k;                      % Координата ц.м., [м]
l_k_1 = 2/3 * l_k;                          % Длина конуса с затуплением, [м]

beta_k_rad = deg2rad(beta_k);               % Значение угла полураствора конуса в радианах, [rad]

%% Параметры атмосферы из файла-функции atmosphere_heopotencial

[g_inf, T_inf, p_inf, rho_inf, a_inf] = atmosphere_heopotencial(H);

T_0 = T_inf * (1 + (k - 1) / 2 * M_inf^2);  % Температура торможения, [К]
V_inf = M_inf * a_inf;                      % Скорость набегающего потока, [м/с]
q_inf = rho_inf * V_inf^2 / 2;              % Скоростной напор, [Па]
a_krit_sqr = (2 * k * R * T_0) / (k + 1);   % Критическая скорость звука, [м/с]

%% Найдем угол СУ для плоского случая с помощью файла функции findThetaShockWave

theta_c_pl = findThetaShockWave(k, M_inf, beta_k);
% theta_c_pl = deg2rad(theta_c_pl);

%% Цикл для определения расчетного угла полураствора конуса и V_r, V_th при alpha = 0

[V_r, V_th, beta_n, theta_c] = Find_Beta_V_r(theta_c_pl, beta_k_rad, V_inf, M_inf, T_0, R, k);

% Вывод данных цикла
disp('Tangent_Component_Of_Velocity'); disp(V_r);
disp('Normal_Component_Of_Velocity'); disp(V_th);
disp('Угол конуса расчетный'); disp(rad2deg(beta_n));
disp('Расчетный угол СУ для острого конуса'); disp(rad2deg(theta_c));

%% Расчет параметров воздуха за СУ (Оформить в функцию)

% Скорость на поверхности конуса
V_k = V_r;

% Температура на поверхности конуса
T_k = Temperature_Surf(T_0, V_k, k, R);

% Число Маха
M_k = V_k / (sqrt(R * T_k * k));

% Расчет давления
p_0_sh = p_inf * ((k + 1) / (2 * k * M_inf^2 * sin(theta_c)^2 - (k - 1)))^(1 / (k - 1)) * ...
    (((k + 1) * M_inf^2 * sin(theta_c)^2 * ((k - 1) * M_inf^2 + 2)) / (2 * (2 + (k - 1) * M_inf^2 * sin(theta_c)^2)))^(k / (k - 1));

p_k = p_0_sh * (1 + (k - 1) / 2 * M_k^2)^(- k / (k - 1));

% Плотность воздуха на поверхности конуса
rho_k = p_k / (R * T_k);

% Коэффициент давления
C_p_k = (p_k - p_inf) / q_inf;

% Вывод
disp('Расчетный коэффициент давления'); disp(C_p_k);
disp('Расчетное число Маха'); disp(M_k);

%% Эмпирические соотношения
% Коэффициент давления (погрешность порядка 5%)
C_p_math = 2 * 10^(-3) * (0.8 + 1 / M_inf^2) * beta_k^(1.7);

% Угол СУ для острого конуса
theta_c_math = asind(sqrt(1 / M_inf^2 + (k + 1) / 2 * sin(beta_k_rad)^2));

% Вывод эмпирических соотношений
disp('Эмпирический коэффициент давления (погрешность 5%)'); disp(C_p_math);
disp('Эмпирический угол СУ для острого конуса'); disp(theta_c_math);

%% Коэффициент продольной силы
% Донное давление
p_d = p_inf * (1 / M_inf - 0.1);

% Донный коеффициент давления
C_p_d = (p_d - p_inf) / q_inf;

C_x = C_p_k + abs(C_p_d);

% Вывод
disp('Коэффициент продольной силы Cx'); disp(C_x);

%% Метод местных конусов

% Создадим строки, заполненные нулями, чтобы записывать туда полученные АДХ
Cx = zeros(1, 41);
Cy = zeros(1, 41);
mz = zeros(1, 41);
Cxa = zeros(1, 41);
Cya = zeros(1, 41);
K = zeros(1, 41);

dalpha = 0.02 * beta_k;                                 % Шаг угла атаки, [deg]
dgamma = 0.1;                                           % Шаг гамма, [rad]
alpha = 0;
i = 1;                                                  % Счетчик цикла

while alpha <= (0.8 * beta_k)
    gamma = 0;
    Int_1 = 0;                                          % Нижний предел первого интеграла
    Int_2 = 0;                                          % Нижний предел второго интеграла
    while gamma < pi
        beta_f = Beta_F(beta_k_rad, deg2rad(alpha), gamma);      % Фиктивный угол полураствора
        C_p_math_a = Pressue_Coefficient(beta_f, M_inf);
        Int_1 = Int_1 + C_p_math_a * dgamma / pi;
        Int_2 = Int_2 + C_p_math_a * cos(gamma) * dgamma / pi;
        gamma = gamma + dgamma;
    end
    Cx(1, i) = Int_1 + abs(C_p_d);
    Cy(1, i) = -1 / tan(beta_k_rad) * Int_2;

    mz_0 = 4 / (sin(2 * beta_k_rad) * 3) * Int_2;
    mz(1, i) = mz_0 + Cy(1, i) * x_cm;

    Cxa(1, i) = Cx(1, i) * cos(deg2rad(alpha)) + Cy(1, i) * sin(deg2rad(alpha));
    Cya(1, i) = - Cx(1, i) * sin(deg2rad(alpha)) + Cy(1, i) * cos(deg2rad(alpha));

    K(1, i) = Cya(1, i) / Cxa(1, i);

    alpha = alpha + dalpha;

    i = i + 1;
end

al = 0:dalpha:(alpha - dalpha);
table = cat(1, al, Cx, Cy, mz, Cxa, Cya, K);            % Обьединяем вектора в таблицу
table = table';                                         % Транспонируем таблицу для вывода

% Вывод
disp('Метод местных конусов | alpha | Cx | Cy | mz | Cxa | Cya | K |');
disp(table);

%% Графики
figure
plot(al, Cx, "o"); grid on;
title('Cx');
xlabel('alpha, [град]');
ylabel('Cx, [б/р]');

figure
plot(al, Cy, "o"); grid on;
title('Cy');
xlabel('alpha, [град]');
ylabel('Cy, [б/р]');

figure
plot(al, mz, "o"); grid on;
title('mz');
xlabel('alpha, [град]');
ylabel('mz, [б/р]');

figure
plot(al, Cxa, "o"); grid on;
title('Cxa');
xlabel('alpha, [град]');
ylabel('Cxa, [б/р]');

figure
plot(al, Cya, "o"); grid on;
title('Cya');
xlabel('alpha, [град]');
ylabel('Cya, [б/р]');

figure
plot(al, K, "o"); grid on;
title('K');
xlabel('alpha, [град]');
ylabel('K, [б/р]');

%% Функции

% Касательная составляющая скорости на СУ
function [V_r_c] = SW_Tangent_Component_Of_Velocity(V_inf, theta_c)
    V_r_c = V_inf * cos(theta_c);
end

% Нормальная составляющая скорости на СУ
function [V_th_c] = SW_Normal_Component_Of_Velocity(V_inf, theta_c, k, M_inf)
    V_th_c = -SW_Tangent_Component_Of_Velocity(V_inf, theta_c) * tan(theta_c) *...
        ((2 + (k - 1) * M_inf^2 * sin(theta_c)^2)/((k + 1) * M_inf^2 * sin(theta_c)^2));
end

% Следующая касательная составляющая скорости за СУ
function [V_r] = Tangent_Component_Of_Velocity(V_r_previous, V_th_previous, dtheta)
    V_r = V_r_previous + V_th_previous * dtheta;
end

% Следующая нормальная составляющая скорости за СУ
function [V_th] = Normal_Component_Of_Velocity(theta_c, V_r_previous, V_th_previous, dtheta, T_0, k, R)
    a_sqr = k * R * T_0 - ((k - 1) / 2) * (V_r_previous^2 - V_th_previous^2);
    dv_dth = (-V_th_previous / tan(theta_c) + V_r_previous * (V_th_previous^2 / a_sqr - 2)) / ...
        (1 - (V_th_previous^2 / a_sqr));
    V_th = V_th_previous + dv_dth * dtheta;
end

% Максимальная скорость
function [V_max_sqr] = V_Maximum_Sqr(T_0, k, R)
    V_max_sqr = (2 * k * R * T_0) / (k - 1);
end

% Температура на поверхности конуса
function [T_k] = Temperature_Surf(T_0, V_k, k, R)
    T_k = T_0 * (1 - V_k^2 / V_Maximum_Sqr(T_0, k, R));
end

% Цикл для определения расчетного угла полураствора конуса и V_r, V_th
function [V_r, V_th, beta_n, theta_c] = Find_Beta_V_r(theta_c_pl, beta_k_rad, V_inf, ...
    M_inf, T_0, R, k)
    coef = 0.85;                                                           % Вводим коэффициент умешения
    dtheta = -0.00001;                                                     % Задаем шаг
    beta_n = deg2rad(theta_c_pl);                                          % Определяем новую переменную, чтобы просто войти в цикл
    while (beta_n - beta_k_rad) / beta_k_rad >= 0.00001
        % Угол СУ в радианах
        theta_c = coef * deg2rad(theta_c_pl);

        % Компоненты скорости за СУ
        V_r_c = SW_Tangent_Component_Of_Velocity(V_inf, theta_c);
        V_th_c = SW_Normal_Component_Of_Velocity(V_inf, theta_c, k, M_inf);

        % Принимаем начальные условия
        V_r = V_r_c;
        V_th = V_th_c;

        theta = theta_c;
        while V_th < 0
            V_r = Tangent_Component_Of_Velocity(V_r, V_th, dtheta);
            V_th = Normal_Component_Of_Velocity(theta, V_r, V_th, dtheta, T_0, k, R);
            theta = theta + dtheta;
        end
        beta_n = theta;
        disp(rad2deg(beta_n))
        coef = coef - 0.00001;
    end
end

% Фиктивный угол полураствора
function [beta_f] = Beta_F(beta_k_rad, alpha, gamma)
    beta_f = beta_k_rad - alpha * cos(gamma) - 0.5 * alpha^2 * sin(gamma)^2 / tan(beta_k_rad);
end

% Эмпирическая зависимость для определения коэффициента давления
function [C_p_math_a] = Pressue_Coefficient(beta_f, M_inf)
    C_p_math_a = 4 * sin(beta_f)^2 * (2.5 + 8 * sqrt(M_inf^2 - 1) * sin(beta_f)) / (1 + 16 * sqrt(M_inf^2 - 1) * sin(beta_f));
end