%% Thin airfoil theory
%
% For an airfoil whose camberline gradient is defined by a cosine series truncated after the 4th coefficient:
% \frac{d y_{c}}{dx} = A_0 + \sum\limits_{n=1}^\infty A_n \cos{n\theta}
% the camber line is given by the following expression (assume constant of integration is null):
% \frac{y_{c}}{c} = A_{0} \frac{1-\cos \theta}{2} + A_{1} \frac{\sin^2 \theta}{4} + A_{2} \Big( \frac{\cos \theta}{2} - \frac{\cos^3 \theta}{3} \Big) + A_{3} \Big( \frac{3 \cos^2 \theta}{4} - \frac{\cos^4 \theta}{2} \Big)
% With boundary conditions: $y_c(0) = 0 \textrm{ and } y_c(\pi) = 0:$
% A_3=-\frac{2}{3}A_2\\
% A_2=3A_0
%

clear all
close all

% play with the values of A0 and A1:
A0 = 0.15; 
A1 = -0.15; 
A2 = 3*A0;
A3 = -2/3*A2;


theta=linspace(0,pi,101);

% different contributions to the camberline:
yc0 = A0*(1-cos(theta))/2;
yc1 = A1*sin(theta).^2;
yc2 = A2*(cos(theta)/2 - cos(theta).^3/3);
yc3 = A3*(3*cos(theta).^2/4 - cos(theta).^4/2);
yc = yc0 + yc1 + yc2 + yc3;

% plots of the camber line
figure, hold on
set(gcf,'DefaultLineLineWidth',2)
tiledlayout(5,1)
ax1 = nexttile; plot(theta,yc0), ylabel('y_{c,0}')
ax2 = nexttile; plot(theta,yc1), ylabel('y_{c,1}')
ax3 = nexttile; plot(theta,yc2), ylabel('y_{c,2}')
ax4 = nexttile; plot(theta,yc3), ylabel('y_{c,3}')
ax5 = nexttile; plot(theta,yc), xlabel('\theta'),ylabel('y_{c}')
linkaxes([ax1 ax2 ax3 ax4 ax5],'xy')
xticks([0 pi/4 pi/2 3*pi/4 pi])
xticklabels({'0','\pi/4','\pi/2','3\pi/4','\pi'})
ax1.XLim = [0 pi];
ax1.YLim = [-0.4 0.4];

alpha = 9.5; % angle of attack in degree

% different contributions to the vortex sheet strength
k0 = (deg2rad(alpha) - A0)*(cos(theta)+1)./sin(theta);
k1 = A1*sin(theta);
k2 = A2*sin(2*theta);
k3 = A3*sin(3*theta);
kdividedby2U = k0 + k1 + k2 + k3;

% plots of the sheet strength
figure, hold on
set(gcf,'DefaultLineLineWidth',2)
tiledlayout(5,1)
bax1 = nexttile; plot(theta,k0), ylabel('k_{0}')
bax2 = nexttile; plot(theta,k1), ylabel('k_{1}')
bax3 = nexttile; plot(theta,k2), ylabel('k_{2}')
bax4 = nexttile; plot(theta,k3), ylabel('k_{3}')
bax5 = nexttile; plot(theta,kdividedby2U), 
xlabel('\theta'),ylabel('k/(2U_{\infty})')
linkaxes([bax1 bax2 bax3 bax4 bax5],'xy')
xticks([0 pi/4 pi/2 3*pi/4 pi])
xticklabels({'0','\pi/4','\pi/2','3\pi/4','\pi'})
bax1.XLim = [0 pi];
bax1.YLim = [-2 5];

% Let focus solely on the camberline and the vortex sheet strength:

clear all
close all

alpha = 10; % angle of attack in degree

% play with the values of A0 and A1:
A0 = 0.1; 
A1 = 0.2; 
A2 = 3*A0;
A3 = -2/3*A2;

theta = linspace(0,pi,101);
xoverc = (1-cos(theta))/2;  

% different contributions to the camberline:
yc0 = A0*(1-cos(theta))/2;
yc1 = A1*sin(theta).^2;
yc2 = A2*(cos(theta)/2 - cos(theta).^3/3);
yc3 = A3*(3*cos(theta).^2/4 - cos(theta).^4/2);
yc = yc0 + yc1 + yc2 + yc3;

% different contributions to the vortex sheet strength
k0 = (deg2rad(alpha) - A0)*(cos(theta)+1)./sin(theta);
k1 = A1*sin(theta);
k2 = A2*sin(2*theta);
k3 = A3*sin(3*theta);
kdividedby2U = k0 + k1 + k2 + k3;

% 
figure, hold on
set(gcf,'DefaultLineLineWidth',2)
tiledlayout(2,1)
ax1 = nexttile; 
plot(xoverc,yc), 
xlabel('x/c'),ylabel('y_{c}')
axis equal
title(sprintf('this camberline at \\alpha = %.1f° yields',alpha))

ax2 = nexttile; 
plot(xoverc,kdividedby2U), 
xlabel('x/c'),ylabel('k/(2U_{\infty})')
linkaxes([ax1 ax2],'x')
ax1.XLim = [0 1];

cn = 2*pi * deg2rad(alpha) + pi*(A1 - 2*A0);
cm1over4 = -pi/4 * (A1 - A2);
title(sprintf('a lift coefficient: C_n = %.2f and a moment coefficient C_{m,1/4} = %.3f',cn,cm1over4))
% the camberline:
yc = @(A0,A1,theta) A0.*(1-cos(theta))/2 + A1.*sin(theta).^2 + 3*A0.*(cos(theta)/2 - cos(theta).^3/3) - 2*A0.*(3*cos(theta).^2/4 - cos(theta).^4/2);
% the vortex sheet strength:
kdividedby2U = @(A0,A1,theta,alpha) (deg2rad(alpha) - A0).*(cos(theta)+1)./sin(theta) + A1.*sin(theta) + 3*A0.*sin(2*theta) - 2*A0.*sin(3*theta);

theta = linspace(0,pi,101);
xoverc = (1-cos(theta))/2; 

figure, hold all
for A1=-0.2:0.1:0.4
plot(xoverc,kdividedby2U(A0,A1,theta,alpha))
end
xlabel('x/c'), ylabel('k/(2U_{\infty})') 
cn =@(A0,A1,A2,alpha) 2*pi * deg2rad(alpha) + pi*(A1 - 2*A0);
cm1over4 =@(A0,A1,A2,alpha) -pi/4 * (A1 - A2);
cmLE = @(A0,A1,A2,alpha) -1/4*( (2*pi * deg2rad(alpha) + pi*(A1 - 2*A0)) + pi*(A1-A2) );
xcp = @(A0,A1,A2,alpha) 1/4 + pi/4/( 2*pi * deg2rad(alpha) + pi*(A1 - 2*A0) ).*(A1-A2);

% keep in mind that for this example: A2 = 3*A0; A3 = -2/3*A2 = -2*A0;

% for 
A0 = 0.1; 
% and 
A1 = -0.2:0.1:0.4;
% and 
alpha = 10;

figure, plot(A1, cn(A0,A1,3*A0,alpha)), xlabel('A_1'), ylabel('c_n')



%
%
%
%