%% Função que fornece os pontos de um perfil naca 4 dígitos
%
% [xy]=NACA4Digits(digit) retorna 99 pontos x e y do perfil definido pelo
% digito
%
% [xy]=NACA4Digits(digit,n) retorna 2*n-1 pontos x e y do perfil definido pelo
% digito, sendo metade no extradorso e a outra no intradorso.
%
% [xy]=NACA4Digits(digit,opt) retorna 2*n-1 pontos x e y do perfil definido pelo
% digito, sendo metade no extradorso e a outra no intradorso. Na sequencia
% executa a opção selecionada
%       opt = 'Plota' - Plota o perfil
%
% [xy]=NACA4Digits(digit,n,opt) retorna 99 pontos x e y do perfil definido pelo
% digito. Na sequencia executa a opção selecionada
%       opt = 'Plota' - Plota o perfil
function[xy,x,yt,yc]=NACA4Digits(digit,n,opt)

switch nargin
    case 1
        n = 50;
        opt = 0;
    case 2
        if isstr(n)
            opt = n;
            n = 50;
        else
            opt = 0;
        end
    case 3
    otherwise
        disp('pân')
end
   
M  = digit(1);
P  = digit(2);
TT = digit(3);

% Calcula a distribuição de espessura 
x  = 0.5-0.5*cos(linspace(0,1,n)*pi);
yt = 5*TT/100*(0.2969*sqrt(x) - 0.1260*x -0.3516*x.^2 + 0.2843*x.^3 - 0.1036*x.^4); % Zero thickness te

% Calcula a linha de arqueamento (Camber)
if P == 0 || M == 0
    yc = yt*0;
    theta = yc;
else
    p = P/10;
    m = M/100;
    
    yc1 = m/p^2*(2*p*x-x.^2);
    yc2 = m/(1-p)^2*((1-2*p)+2*p*x-x.^2);
    yc  = (x<p).*yc1 + ~(x<p).*yc2;
    
    % Calcula a derivada do camber
    dyc1dx = 2*m/p^2*(p-x);
    dyc2dx = 2*m/(1-p)^2*(p-x);
    dycdx  = (x<p).*yc1 + ~(x<p).*yc2;
    
    % Calcula o theta
    theta = atan(dycdx);
end


% Calcula os pontos do extradorso e intradorso
xu =  x - yt.*sin(theta);
yu = yc + yt.*cos(theta);
xl =  x + yt.*sin(theta);
yl = yc - yt.*cos(theta);

% Junta tudo em um vetor
xy = [xu(end:-1:2)' yu(end:-1:2)'; xl' yl' ];
xy(1,:) = [1 0];
xy(end,:) = [1 0];
% xy(1 ,:) = [1 0];
% xy(end ,:) = [1 0];

switch opt
    case 'Plota'
        plot(xu,yu,'b')
        hold on
        plot(xl,yl,'r')
        plot(x,yc,'k')
        
        grid on
        axis equal
        
        axis([0 1 -TT*1.5/100 TT*1.5/100])
        title(sprintf('Perfil NACA- %1.2f %1.2f %2.2f',digit))
end