function [y] = naca4_generator(name,n)
% NACA4_GENERATOR genera un profilo NACA a 4 cifre
% Scritto da Pavan Andrea - 09/12/2020
%       [y] = naca4_generator(name,n)
% INPUT:
%       name: nome profilo (es: '4412')
%       n: numero pannelli (es: 100)
% OUTPUT:
%       y: matrice nx2 dei punti

m = str2double(name(1))/100;
p = str2double(name(2))/10;
t = str2double(name(3:4))/100;
c = 1;
Nd = floor(n/2);  %numero punti dorso
Nv = n-Nd;  %numero punti ventre

yt = @(x) 5*t*c*(0.2969*sqrt(x/c) -0.1260*(x/c) -0.3516*(x/c).^2 +0.2843*(x/c).^3 -0.1036*(x/c).^4);  %equazione metà spessore
yc1 = @(x) m*x.*(2*p-x/c)/(p^2);  %equazione camber per 0<=x<=pc
yc2 = @(x) m*(c-x).*(1-2*p+x/c)/((1-p)^2);  %equazione camber per pc<=x<=c
dyc1 = @(x) 2*m*(p-x/c)/(p^2);  %equazione derivata camber per 0<=x<=pc
dyc2 = @(x) 2*m*(p-x/c)/((1-p)^2);  %equazione derivata camber per pc<=x<=c
if abs(p)<eps || abs(m)<eps
    yc1 = @(x) 0;
    yc2 = @(x) 0;
    dyc1 = @(x) 0;
    dyc2 = @(x) 0;
end

id = 1:(Nd-2+1);
xd = 0.5*c + 0.5*c*cos((2*id-1)*pi/(2*(Nd-2+1)));  %nodi di Chebyshev al dorso
xd = [0; flipud(xd(:)); c];
iv = 1:(Nv-2);
xv = 0.5*c + 0.5*c*cos((2*iv-1)*pi/(2*(Nv-2)));  %nodi di Chebyshev al ventre
xv = [0; flipud(xv(:)); c];

xd1 = xd(xd<=p*c);
xd2 = xd(xd>p*c);
xv1 = xv(xv<=p*c);
xv2 = xv(xv>p*c);

ytd1 = yt(xd1);
ytd2 = yt(xd2);
ytv1 = yt(xv1);
ytv2 = yt(xv2);

ycd1 = yc1(xd1);
ycd2 = yc2(xd2);
ycv1 = yc1(xv1);
ycv2 = yc2(xv2);

thetad1 = atan(dyc1(xd1));
thetad2 = atan(dyc2(xd2));
thetav1 = atan(dyc1(xv1));
thetav2 = atan(dyc2(xv2));

xd1 = xd1 - ytd1.*sin(thetad1);
xd2 = xd2 - ytd2.*sin(thetad2);
yd1 = ycd1 + ytd1.*cos(thetad1);
yd2 = ycd2 + ytd2.*cos(thetad2);
xv1 = xv1 + ytv1.*sin(thetav1);
xv2 = xv2 + ytv2.*sin(thetav2);
yv1 = ycv1 - ytv1.*cos(thetav1);
yv2 = ycv2 - ytv2.*cos(thetav2);
xd = [xd1; xd2];
yd = [yd1; yd2];
xv = [xv1; xv2];
yv = [yv1; yv2];

y = zeros(n,2);
y(1:Nd,1) = flipud(xd(2:end));  %x dorso
y((Nd+1):end,1) = xv;  %x ventre
y(1:Nd,2) = flipud(yd(2:end));  %y dorso
y((Nd+1):end,2) = yv;  %y ventre
y(1,:) = [1, 0];
y(end,:) = [1,0];
