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

beta_k = 15;                  % Угол полураствора конуса, [deg]
M_inf = 3.397;                % Число Маха набегающего потока, [б/р]
g_0 = 9.8066;                 % Ускорение свободного падения, [м/с^2]
p_inf = 0.098 * g_0 * 10^4;   % Давление набегающего потока (статическое давление), [Па]
p_0 = 6.45 * g_0 * 10^4;      % Давление полного торможения (давление в форкамере), [Па]
T_0 = 288.15;                 % Давление полного торможения, [К]
alpha = 0;                    % Угол атаки, [deg]

k = 1.4;                      % Отношение удельных теплоемкомстей воздуха [б/р]
R = 287.05287;                % Газовая постоянная, Дж/(кг*К)
x_cm = 2/3;                   % Положение ц.м. конуса от носка к длине конуса, [б/р]

beta_k_rad = deg2rad(beta_k); % Переведем угол полураствора конуса в радианы

%% Параметры потока

T_inf = T_0 * (1 + (k - 1) / 2 * M_inf^2)^(-1);     % Температура набегающего потока, [К]
V_inf = M_inf * sqrt(k * R * T_inf);                % Скорость набегающего потока, [м/с]
q_inf = ((p_inf / (R * T_inf)) * V_inf^2) / 2;      % Скоростной напор, [Па]

%% Донное давление

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

%% Результаты эксперимента
alpha_list = [3 6 9 12];
gamma_list = 0:(pi/6):pi;
Pressue_3deg  = [0.1891 0.1917 0.2104 0.2269 0.2499 0.2746 0.2771];
Pressue_6deg  = [0.1551 0.1612 0.1867 0.2261 0.2805 0.3231 0.3376];
Pressue_9deg  = [0.1266 0.1367 0.1675 0.2246 0.2970 0.3625 0.4019];
Pressue_12deg = [0.1074 0.1158 0.1502 0.2256 0.3222 0.4170 0.4628];

% Переводим в Па
Pressue_3deg_Pa  = Pressue_3deg .* g_0 .* 10^4;
Pressue_6deg_Pa  = Pressue_6deg .* g_0 .* 10^4;
Pressue_9deg_Pa  = Pressue_9deg .* g_0 .* 10^4;
Pressue_12deg_Pa = Pressue_12deg .* g_0 .* 10^4;

% Рассчитываем коэффициенты давления
Pressue_Coef_3deg  = (Pressue_3deg_Pa - p_inf) ./ q_inf;
Pressue_Coef_6deg  = (Pressue_6deg_Pa - p_inf) ./ q_inf;
Pressue_Coef_9deg  = (Pressue_9deg_Pa - p_inf) ./ q_inf;
Pressue_Coef_12deg = (Pressue_12deg_Pa - p_inf) ./ q_inf;


% Обьединим все в одну таблицу
Pressue_Coef = cat(1, Pressue_Coef_3deg, Pressue_Coef_6deg, Pressue_Coef_9deg, Pressue_Coef_12deg);

%disp('Коэффициенты давления | 3 | 6 | 9 | 12 |'); disp(Pressue_Coef');

%% Обработка результатов эксперимента методом трапеций

Integral_1 = zeros();
Integral_2 = zeros();
Cx = zeros();
Cy = zeros();
mz = zeros();

Cxa = zeros();
Cya = zeros();
K   = zeros();

for i = 1:4
    Integral_1(i) = pi / 6 * ((Pressue_Coef(i, 1) + Pressue_Coef(i, 7)) / 2 + Pressue_Coef(i, 2) + ...
        Pressue_Coef(i, 3) + Pressue_Coef(i, 4) + Pressue_Coef(i, 5) + Pressue_Coef(i, 6));

    Integral_2(i) = pi / 6 * ((Pressue_Coef(i, 1) * cos(gamma_list(1)) + Pressue_Coef(i, 7) * cos(gamma_list(7))) / 2 + ...
        Pressue_Coef(i, 2) * cos(gamma_list(2)) + Pressue_Coef(i, 3) * cos(gamma_list(3)) + Pressue_Coef(i, 4) * cos(gamma_list(4)) + ...
        Pressue_Coef(i, 5) * cos(gamma_list(5)) + Pressue_Coef(i, 6) * cos(gamma_list(6)));

    Cx(i) = Integral_1(i) / pi + abs(C_p_d);
    Cy(i) = (Integral_2(i) / pi) * (-1 / (tan(beta_k_rad)));
    mz(i) = (Integral_2(i) / pi) * (4 / (3 * sin(2 * beta_k_rad)));

    Cxa(i) =  Cx(i) * cosd(alpha_list(i)) + Cy(i) * sind(alpha_list(i));
    Cya(i) = -Cx(i) * sind(alpha_list(i)) + Cy(i) * cosd(alpha_list(i));
    K(i)   =  Cya(i) / Cx(i);
end

%% Расчет теоретических АДХ
% Далее код для определения АДХ конуса методом последовательных приближений
% и методом местных конусов
%% Найдем угол СУ для плоского случая с помощью файла функции findThetaShockWave
theta_c_pl = findThetaShockWave(k, M_inf, beta_k);

% Цикл для определения расчетного угла полураствора конуса и 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);

%% Расчет параметров воздуха за СУ

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

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

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

C_x = zeros(1, 5);
C_y = zeros(1, 5);
m_z = zeros(1, 5);

C_xa = zeros(1, 5);
C_ya = zeros(1, 5);
K_k = zeros(1, 5);

dalpha = 3;                                             % Шаг угла атаки, [deg]
dgamma = deg2rad(0.1);                                  % Шаг гамма, [rad]
i = 1;                                                  % Счетчик цикла

while alpha <= 12
    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
    C_x(1, i) = Int_1 + abs(C_p_d);
    C_y(1, i) = -1 / tan(beta_k_rad) * Int_2;

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

    C_xa(1, i) = C_x(1, i) * cos(deg2rad(alpha)) + C_y(1, i) * sin(deg2rad(alpha));
    C_ya(1, i) = - C_x(1, i) * sin(deg2rad(alpha)) + C_y(1, i) * cos(deg2rad(alpha));

    K_k(1, i) = C_ya(1, i) / C_xa(1, i);

    alpha = alpha + dalpha;

    i = i + 1;
end

al = 0:dalpha:12;
table = cat(1, al, C_x, C_y, m_z, C_xa, C_ya, K_k);            % Обьединяем вектора в таблицу
table = table';                                         % Транспонируем таблицу для вывода

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


%% Вывод результатов обработки эксперимента

% Обьединим полученные данные в одну таблицу
table_data = cat(1, alpha_list, Cx, Cy, mz, Cxa, Cya, K);
table_data = table_data';

disp('Вывод результатов обработки эксперимента | Cx | Cy | mz | Cxa | Cya | K |'); disp(table_data);

%% Графики

p_1_1 = polyfit(al,C_x,2);
p_1_2 = polyfit(alpha_list,Cx,2);

x_all = linspace(0,12);
y_1_1 = polyval(p_1_1,x_all);
y_1_2 = polyval(p_1_2,x_all);
figure
plot(al, C_x, "or", alpha_list, Cx, "sb"); grid on;
hold on
plot(x_all, y_1_1, "--r", x_all, y_1_2, "b");
hold off
ylim([0 0.3]);
title('Cx');
xlabel('alpha, [град]');
ylabel('Cx, [б/р]');
legend({'Теория','Эксперимент'},'Location','southeast');

p_1_1 = polyfit(al,C_y,2);
p_1_2 = polyfit(alpha_list,Cy,2);

x_all = linspace(0,12);
y_1_1 = polyval(p_1_1,x_all);
y_1_2 = polyval(p_1_2,x_all);
figure
plot(al, C_y, "or", alpha_list, Cy, "sb"); grid on;
hold on
plot(x_all, y_1_1, "--r", x_all, y_1_2, "b");
hold off
title('Cy');
xlabel('alpha, [град]');
ylabel('Cy, [б/р]');
legend({'Теория','Эксперимент'},'Location','northwest');

p_1_1 = polyfit(al,m_z,2);
p_1_2 = polyfit(alpha_list,mz,2);

x_all = linspace(0,12);
y_1_1 = polyval(p_1_1,x_all);
y_1_2 = polyval(p_1_2,x_all);
figure
plot(al, m_z, "or", alpha_list, mz, "sb"); grid on;
hold on
plot(x_all, y_1_1, "--r", x_all, y_1_2, "b");
hold off
title('mz');
xlabel('alpha, [град]');
ylabel('mz, [б/р]');
legend({'Теория','Эксперимент'},'Location','northeast');

p_1_1 = polyfit(al, C_xa, 2);
p_1_2 = polyfit(alpha_list, Cxa, 2);

x_all = linspace(0,12);
y_1_1 = polyval(p_1_1, x_all);
y_1_2 = polyval(p_1_2, x_all);
figure
plot(al, C_xa, "or", alpha_list, Cxa, "sb"); grid on;
hold on
plot(x_all, y_1_1, "--r", x_all, y_1_2, "b");
hold off
title('Cxa');
xlabel('alpha, [град]');
ylabel('Cxa, [б/р]');
legend({'Теория','Эксперимент'},'Location','northwest');

p_1_1 = polyfit(al, C_ya, 2);
p_1_2 = polyfit(alpha_list, Cya, 2);

x_all = linspace(0,12);
y_1_1 = polyval(p_1_1, x_all);
y_1_2 = polyval(p_1_2, x_all);
figure
plot(al, C_ya, "or", alpha_list, Cya, "sb"); grid on;
hold on
plot(x_all, y_1_1, "--r", x_all, y_1_2, "b");
hold off
title('Cya');
xlabel('alpha, [град]');
ylabel('Cya, [б/р]');
legend({'Теория','Эксперимент'},'Location','northwest');

p_1_1 = polyfit(al, K_k, 2);
p_1_2 = polyfit(alpha_list, K, 2);

x_all = linspace(0,12);
y_1_1 = polyval(p_1_1, x_all);
y_1_2 = polyval(p_1_2, x_all);
figure
plot(al, K_k, "or", alpha_list, K, "sb"); grid on;
hold on
plot(x_all, y_1_1, "--r", x_all, y_1_2, "b");
hold off
title('K');
xlabel('alpha, [град]');
ylabel('K, [б/р]');
legend({'Теория','Эксперимент'},'Location','northwest');

%% Функции

% Касательная составляющая скорости на СУ
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

% Цикл для определения расчетного угла полураствора конуса и 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 [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

% Фиктивный угол полураствора
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