CST_airfoil(wl,wu,dz,N);

function [coord] = CST_airfoil(wl,wu,dz,N);

% Description : Create a set of airfoil coordinates using CST parametrization method 
% Input  : wl = CST weight of lower surface
%          wu = CST weight of upper surface
%          dz = trailing edge thickness
% Output : coord = set of x-y coordinates of airfoil generated by CST


% Create x coordinate
x=ones(N+1,1);y=zeros(N+1,1);zeta=zeros(N+1,1);
for i=1:N+1
    zeta(i)=2*pi/N*(i-1);
    x(i)=0.5*(cos(zeta(i))+1);
end

% N1 and N2 parameters (N1 = 0.5 and N2 = 1 for airfoil shape)
N1 = 0.5;
N2 = 1;

zerind = find(x(:,1) == 0); % Used to separate upper and lower surfaces

xl= x(1:zerind-1); % Lower surface x-coordinates
xu = x(zerind:end); % Upper surface x-coordinates

[yl] = ClassShape(wl,xl,N1,N2,-dz); % Call ClassShape function to determine lower surface y-coordinates
[yu] = ClassShape(wu,xu,N1,N2,dz);  % Call ClassShape function to determine upper surface y-coordinates

y = [yl;yu]; % Combine upper and lower y coordinates

coord = [x y]; % Combine x and y into single output

%% Function to calculate class and shape function
function [y] = ClassShape(w,x,N1,N2,dz);


% Class function; taking input of N1 and N2
for i = 1:size(x,1)
    C(i,1) = x(i)^N1*((1-x(i))^N2);
end

% Shape function; using Bernstein Polynomials
n = size(w,2)-1; % Order of Bernstein polynomials

for i = 1:n+1
     K(i) = factorial(n)/(factorial(i-1)*(factorial((n)-(i-1))));
end

for i = 1:size(x,1)
    S(i,1) = 0;
    for j = 1:n+1
        S(i,1) = S(i,1) + w(j)*K(j)*x(i)^(j-1)*((1-x(i))^(n-(j-1)));
    end
end

% Calculate y output
for i = 1:size(x,1)
   y(i,1) = C(i,1)*S(i,1) + x(i)*dz;
end