%% Matlab solution for EX04 - 2
clear
close all

load('NACA0015.mat') % Load pressure and airfoil data
%%
% Airfoil parameters:
rho = 1.2; % Density [kg/m^3]
U_inf = 30; % Flow velocity [m/s]
c = 300; % Chord length [mm]

% Normalize pressure:
norm = 0.5 * rho * U_inf^2;
cp_down = p_down / norm;
cp_up = p_up / norm;

%% Calculate aerodynamic coefficients:
for i = 1:length(alpha)
    Cn(i) = (trapz(x_down,cp_down(:,i)) - trapz(x_up,cp_up(:,i))) / c; %#ok<*SAGROW>
    Cl(i) = Cn(i) * cosd(alpha(i));
    Cm_025c(i) = (trapz(x_up,cp_up(:,i).*(x_up/c-0.25)) - trapz(x_down,cp_down(:,i).*(x_down/c-0.25))) / c;
end

% Gradients of dCn_dalpha and dCm_dalpha at linear parts:
linfit1 = polyfit(alpha(15:40),Cn(15:40),1);
dCn_dalpha = linfit1(1);
linfit2 = polyfit(alpha(15:43),Cm_025c(15:43),1);
dCm_025c_dalpha = linfit2(1);

[CLmax, idx] = max(Cl);
alpha_CLmax = alpha(idx);

AC = 0.25 - dCm_025c_dalpha / dCn_dalpha; % aerodynamic centre

%% Plots:
figure
plot(alpha,Cn), grid on,
hold on, plot(alpha(10:45),polyval(linfit1,alpha(10:45)))
xlabel('\alpha'), ylabel('C_n')

%%
figure
plot(alpha,Cl), grid on
xlabel('\alpha'), ylabel('C_l')
%%
figure
plot(alpha,Cm_025c), grid on,
hold on, plot(alpha(7:47),polyval(linfit2,alpha(7:47)))
xlabel('\alpha'), ylabel('C_{m,c/4}')