Create an airfoil in MATLAB

%%%%% part 1
t = .18;
tau = .349;
e = .0721;
a = .574;
k = 1.889;
c = 2;
theta = [0:1:360];
r1 = sqrt(((a*cos(theta)-a).^2)+(a.^2)*sin(theta).^2)
r2 = sqrt(((a*cos(theta)-e*a).^2)+(a.^2)*sin(theta).^2)
theta1 = atand((a*sind(theta))./(a*cosd(theta)-a)) + pi
for theta2 = atand((a*sind(theta))./(a*cosd(theta)-e*a)) + n*pi
    if theta2 <= 90
        n = 0
    elseif 90<theta2<=270
        n = 1
    elseif 270<theta2<360
        n = 2
    end
end
theta2 = atand((a*sind(theta))./(a*cosd(theta)-e*a)) + n*pi
x = ((r1.^k)/r2.^k-1).*(cos(k*theta1)*cos(k-1).*theta2+sin(k*theta1).*sin(k-1).*theta2)
z = ((r1.^k)/r2.^k-1).*(sin(k*theta1)*cos(k-1).*theta2-cos(k*theta1).*sin(k-1).*theta2)
plot(x,z)

(It should look like this but does not any help would be appreciated)

[KSSV]

clc; clear all ;
clear;
%%%%% part 1
n = 0 ;
t = .18;
tau = .349;
e = .0721;
a = .574;
k = 1.889;
c = 2;
theta = 0:1:360;
r1 = sqrt(((a*cos(theta)-a).^2)+(a.^2)*sin(theta).^2) ;
r2 = sqrt(((a*cos(theta)-e*a).^2)+(a.^2)*sin(theta).^2) ;
theta1 = atand((a*sind(theta))./(a*cosd(theta)-a)) + pi ;
for theta2 = atand((a*sind(theta))./(a*cosd(theta)-e*a)) + n*pi 
    if theta2 <= 90
        n = 0 ;
    elseif theta2>90 && theta2<=270
        n = 1 ;
    elseif theta2>270 && theta2<360
        n = 2 ;
    end
end
theta2 = atand((a*sind(theta))./(a*cosd(theta)-e*a)) + n*pi ;
x = ((r1.^k)/r2.^k-1).*(cos(k*theta1)*cos(k-1).*theta2+sin(k*theta1).*sin(k-1).*theta2) ;
z = ((r1.^k)/r2.^k-1).*(sin(k*theta1)*cos(k-1).*theta2-cos(k*theta1).*sin(k-1).*theta2) ;
plot(x,z) 

Are you syre the equations are correct? You have alos messedup with the units. Check your units of angls while using sin and cos.

[Sam Chak]

Hi @James Hall-Prior

Do you have the 4-digit NACA code?

[James Hall-Prior]

I don't and I'm not sure if I'm allowed to use it

Answer 1

[Davide Masiello]

From the wikipedia page

https://en.wikipedia.org/wiki/NACA_airfoil

t  = 0.25;      % maximum thickness as a fraction of the chord
x0 = 0:0.001:1;
y0 = 5*t*(0.2969*sqrt(x0)-0.1260*x0-0.3516*x0.^2+0.2843*x0.^3-0.1015*x0.^4);
x  = [x0,flip(x0)];
y  = [y0,flip(-y0)];
plot(x,y)
axis equal