%% Исходные данные

N = 4;              % Номер варианта
c = 80;             % Высота профиля крыла, [мм]
b = 1000;           % Длина профиля крыла, [мм]
T_st = 405;         % Температура охлаждаемой стенки, [К]
Re_krit = 5 * 10^6; % Критическое число Рейнольдса
k = 1.4;            % Отношение удельных теплоемкостей воздуха

% Высота полета (геометрическая), [м]
h = (10 + 0.5 * N) * 10^3;

% Число Маха
M_inf = 4 + 0.15 * N;

% Угол атаки
alpha_deg = 2 + 0.1 * N;    % [deg]
alpha = deg2rad(alpha_deg); % [rad]

%% Определим параметры стандартной атмосферы по ГОСТ 4401-81

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

%% Угол полураствора клина

beta_k = atan(c / b); % [deg]
beta_k_deg = rad2deg(beta_k);     % [rad]

%% Области 1, 2

% Определим угол СУ
theta_1_deg = findThetaShockWave(k, M_inf, (beta_k_deg - alpha_deg)); % [deg]
theta_2_deg = findThetaShockWave(k, M_inf, (beta_k_deg + alpha_deg)); % [deg]
theta_1 = deg2rad(theta_1_deg); % [rad]
theta_2 = deg2rad(theta_2_deg); % [rad]

% Давление за СУ
p_1 = Pressure_Behind_SW(p_inf, k, M_inf, theta_1_deg);
p_2 = Pressure_Behind_SW(p_inf, k, M_inf, theta_2_deg);

% Плотность за СУ
rho_1 = Density_Behind_SW(rho_inf, k, M_inf, theta_1);
rho_2 = Density_Behind_SW(rho_inf, k, M_inf, theta_2);

% Температура за СУ
T_1 = Temperature_Behind_SW(T_inf, p_inf, rho_inf, p_1, rho_1);
T_2 = Temperature_Behind_SW(T_inf, p_inf, rho_inf, p_2, rho_2);

% Число Маха за СУ
M_1 = Mach_Behind_SW(k, M_inf, theta_1);
M_2 = Mach_Behind_SW(k, M_inf, theta_2);

% Угол Маха
mu_1 = asin(1 / M_1);
mu_2 = asin(1 / M_2);
mu_1_deg = rad2deg(mu_1);
mu_2_deg = rad2deg(mu_2);

% Давление полного торможения
p_0_1 = Pressure_Full_Brake(p_1, k, M_1);
p_0_2 = Pressure_Full_Brake(p_2, k, M_2);

% Плотность полного торможения
rho_0_1 = Density_Full_Brake(rho_1, k, M_1);
rho_0_2 = Density_Full_Brake(rho_2, k, M_2);

% Температура полного торможения
T_0_1 = Temperature_Full_Brake(T_1, k, M_1);
T_0_2 = Temperature_Full_Brake(T_2, k, M_2);

% Шапка таблицы
heading = ["area", "beta, [deg]", "theta, [deg]", "p, [Pa]", "rho, [kg/m^3]",...
    "T, [K]", "M, [-]", "p_0, [Pa]", "rho_0, [kg/m^3]", "T_0, [K]"];


% Создадим 2 строки, в которые запишем, полученные результаты
area_1 = [1, round((beta_k_deg - alpha_deg), 2), round(theta_1_deg, 2), ...
    round(p_1, 2), round(rho_1, 4), round(T_1, 2), round(M_1, 2), ...
    round(p_0_1, 2), round(rho_0_1, 4), round(T_0_1, 2)];

area_2 = [2, round((beta_k_deg + alpha_deg), 2), round(theta_2_deg, 2), ...
    round(p_2, 2), round(rho_2, 4), round(T_2, 2), round(M_2, 2), ...
    round(p_0_2, 2), round(rho_0_2, 4), round(T_0_2, 2)];

% Таблица
table = cat(1, heading, area_1, area_2);

%% Области 3, 4. Течение Прандтля-Майера (прямая задача)

% Угол поворота потока
beta_turn = 2 * beta_k;

% Фиктивные углы
omega_1 = Fictious_Angle_Omega(k, M_1);
omega_2 = Fictious_Angle_Omega(k, M_2);
omega_3 = omega_1 + beta_turn;
omega_4 = omega_2 + beta_turn;

% Числа Маха в областях 3, 4 через omega_3 и omega_4
M_3 = Mach_Through_Omega(k, omega_3, M_1);
M_4 = Mach_Through_Omega(k, omega_4, M_2);

% Угол Маха
mu_3 = asin(1 / M_3);
mu_4 = asin(1 / M_4);
mu_3_deg = rad2deg(mu_3);
mu_4_deg = rad2deg(mu_4);

% Параметры торможения в области 3 и 4 такие же, как в области 1 и 2
p_0_3 = p_0_1;
p_0_4 = p_0_2;
rho_0_3 = rho_0_1;
rho_0_4 = rho_0_2;
T_0_3 = T_0_1;
T_0_4 = T_0_2;

% Давление
p_3 = Pressure_Behind_PM(p_0_3, k, M_3);
p_4 = Pressure_Behind_PM(p_0_4, k, M_4);

% Плотность
rho_3 = Density_Behind_PM(rho_0_3, k, M_3);
rho_4 = Density_Behind_PM(rho_0_4, k, M_4);

% Температура
T_3 =  Temperature_Behind_PM(T_0_3, k, M_3);
T_4 =  Temperature_Behind_PM(T_0_4, k, M_4);

% Создадим 2 строки, в которые запишем, полученные результаты
area_3 = [3, round(rad2deg(beta_turn), 2), "---", ...
    round(p_3, 2), round(rho_3, 4), round(T_3, 2), round(M_3, 2), ...
    round(p_0_3, 2), round(rho_0_3, 4), round(T_0_3, 2)];

area_4 = [4, round(rad2deg(beta_turn), 2), "---", ...
    round(p_4, 2), round(rho_4, 4), round(T_4, 2), round(M_4, 2), ...
    round(p_0_4, 2), round(rho_0_4, 4), round(T_0_4, 2)];

% Обновляем таблицу
table = cat(1, table, area_3, area_4);

%% Области 5 и 6

[theta_5, theta_6, delta] = Find_Theta_Delta_Bottom(k, M_3, M_4, p_3, p_4, beta_k_deg);


%% Функции

% Давление за СУ
function [p] = Pressure_Behind_SW(p_inf, k, M_inf, theta)
    p = p_inf * ((2 * k) / (k + 1) * M_inf^2 * sind(theta)^2 - (k - 1) / (k + 1));
end

% Плотность за СУ
function [rho] = Density_Behind_SW(rho_inf, k, M_inf, theta)
    rho = rho_inf * (0.5 * (k + 1) * M_inf^2 * sin(theta)^2) / ...
        (1 + 0.5 * (k - 1) * M_inf^2 * sin(theta)^2);
end

% Температура за СУ
function [T] = Temperature_Behind_SW(T_inf, p_inf, rho_inf, p, rho)
    T = T_inf * (p / p_inf) * (rho_inf / rho);
end

% Число Маха
function [M] = Mach_Behind_SW(k, M_inf, theta)
    M = sqrt((2 + (k - 1) * M_inf^2) / (2 * k * M_inf^2 * sin(theta)^2 - (k - 1)) + ...
        (2 * M_inf^2 * cos(theta)^2) / (2 + (k - 1) * M_inf^2 * sin(theta)^2));
end

% Давление полного торможения
function [p_0] = Pressure_Full_Brake(p, k, M)
    p_0 = p * (1 + (k - 1) / 2 * M^2)^(k / (k - 1));
end

% Плотность полного торможения
function [rho_0] = Density_Full_Brake(rho, k, M)
    rho_0 = rho * (1 + (k - 1) / 2 * M^2)^(1 / (k - 1));
end

% Температура полного торможения
function [T_0] = Temperature_Full_Brake(T, k, M)
    T_0 = T * (1 + (k - 1) / 2 * M^2);
end

% Фиктивный угол омега
function [omega] = Fictious_Angle_Omega(k, M)
    omega = sqrt((k + 1) / (k - 1)) * atan(sqrt((k - 1) / (k + 1) * (M^2 - 1))) - ...
        atan(sqrt(M^2 - 1));
end

% Числа Маха в областях 3, 4 через omega_3 и omega_4
function [Mach] = Mach_Through_Omega(k, omega_3_4, M_1_2)
    myfun = @(M)(omega_3_4 - Fictious_Angle_Omega(k, M));
    Mach = fsolve(myfun, M_1_2);
end

% Давление в 3й и 4й области
function [p] = Pressure_Behind_PM(p_0, k, M)
    p = p_0 * (1 + (k - 1) / 2 * M^2)^(-k / (k - 1));
end

% Плотность в 3й и 4й области
function [rho] = Density_Behind_PM(rho_0, k, M)
    rho = rho_0 * (1 + (k - 1) / 2 * M^2)^(-1 / (k - 1));
end

% Температура в 3й и 4й области
function [T] = Temperature_Behind_PM(T_0, k, M)
    T = T_0 * (1 + (k - 1) / 2 * M^2)^(-1);
end

% Функция вычисления углов СУ в 5й и 6й области
function [theta_5, theta_6, delta] = Find_Theta_Delta_Bottom(k, M_3, M_4, p_3, p_4, beta_k_deg)

    % Пусть x = [x(1) -> theta_5, x(2) -> theta_6, x(3) -> delta]
    System_Of_Equations = @(x) [x(1) - findThetaShockWave(k, M_3, beta_k_deg - x(3));
                                x(2) - findThetaShockWave(k, M_4, beta_k_deg + x(3));
                                Pressure_Behind_SW(p_3, k, M_3, x(1)) - Pressure_Behind_SW(p_4, k, M_4, x(2))];
    % Решение системы уравнений
    solution = fsolve(System_Of_Equations, [0, 0, 0]);

    theta_5 = solution(1);
    theta_6 = solution(2);
    delta   = solution(3);
end