%% Thin airfoil theory for 4-digit NACA airfoils

% 4-digit NACA airfoils


m=2/100;
p=4/10;

x1=0:0.01:p;
x2=p+0.01:0.01:1;
x=[x1,x2];
theta = acos(1-2*x);

yc1 = m/p^2*(2*p*x1-x1.^2); % yc and x normalised by c
yc2 = m/(1-p)^2*( (1-2*p) + 2*p*x2 -x2.^2); % yc and x normalised by c
yc=[yc1,yc2];

figure, plot(x, yc), axis equal, axis([0 1 0 0.2])
xlabel('x/c'), ylabel('y/c'), grid on
title(sprintf('NACA%i%iXX airfoil',m*100,p*10))

thetap= acos(1-2*p);

A0 = m/(pi*p^2)*((2*p-1)*thetap + sin(thetap)) + m/(pi*(1-p)^2)*((2*p-1)*(pi-thetap)-sin(thetap));
A1 = 2*m/(pi*p^2) * ( (2*p-1)*sin(thetap) + 0.25*sin(2*thetap) + thetap/2) - ...
    ( 2*m/pi/(1-p)^2 * ((2*p-1)*sin(thetap) + 0.25*sin(2*thetap) - 0.5*(pi-thetap)) );
A2 = 2*m/(pi*p^2) * ( (2*p-1)*0.5*sin(2*thetap) + sin(thetap) - 2/3*sin(thetap)^3) - ...
    ( 2*m/pi/(1-p)^2 * ((2*p-1)*0.5*sin(2*thetap) + sin(thetap) - 2/3*sin(thetap)^3) );

figure, hold all
plot(x, gradient(yc,x),'k-'), 
plot(x, A0 + A1*cos(theta) +  A2*cos(2*theta),'r-')
xlabel('x/c'), ylabel('dy_c/dx'), grid on
title(sprintf('NACA%i%iXX airfoil',m*100,p*10))
legend('original','truncated sine series')

% Now let's look at the truncated vortex sheet strength for this 4-digit NACA airfoil:

% the vortex sheet strength:
kdividedby2U = @(alpha) (deg2rad(alpha) - A0).*(cos(theta)+1)./sin(theta) + A1.*sin(theta) + A2.*sin(2*theta);
 
figure, hold all
for alpha=0:1:10
plot(x,kdividedby2U(alpha))
end
xlabel('x/c'), ylabel('k/(2U_{\infty})') 
title(sprintf('NACA%i%iXX airfoil',m*100,p*10))

% Now let's look at the lift and moment polar for this 4-digit NACA airfoil:

cn =@(alpha) 2*pi * deg2rad(alpha) + pi*(A1 - 2*A0);
cm1over4 =@(alpha) -pi/4 * (A1 - A2);
cmLE = @(alpha) -1/4.*( (2*pi * deg2rad(alpha) + pi*(A1 - 2*A0)) + pi*(A1-A2) );
xcp = @(alpha) 1/4 + pi/4./( 2*pi * deg2rad(alpha) + pi*(A1 - 2*A0) ).*(A1-A2);

alpha=0:1:10;
figure, plot(alpha, cn(alpha))
figure, plot(alpha, cmLE(alpha))
figure, plot(alpha, xcp(alpha))