function [Cl,Cm,cp,lambda,Cd,time] = panel1(airfoil,alpha)
% PANEL1 determina distribuzione di pressione e coefficienti aerodinamici
% di un profilo alare 2D. Metodo a pannelli.
% Tipo pannelli flat, tipo vortici linear.
% Condizione impermeabilità formulata sulla velocità, Neumann BC.
% Scritto da Pavan Andrea - 29/10/2019
%       [Cl,Cm,cp,lambda,Cd,time] = panel1(airfoil,alpha)
% INPUT:
%       airfoil: matrice Nx2 dei punti del profilo (TE->upper->LE->lower)
%       alpha: angolo di attacco profilo [°]
% OUTPUT:
%       Cl: coefficiente di portanza unitario
%       Cm: coefficiente di momento (rispetto al quarto di corda) unitario
%       cp: matrice Nx2 dei coefficienti di pressione [x,cp]
%       lambda: vettore singolarità (vorticità sui nodi)
%       Cd: coefficiente di resistenza unitario (errore commesso)
%       time: tempo impiegato per i calcoli

T0 = clock;

%angolo di attacco profilo
alphaRad = alpha*pi/180;
Rv = [cos(-alphaRad), -sin(-alphaRad); sin(-alphaRad), cos(-alphaRad)];
x = Rv*airfoil';
x = x';
vinf = [1; 0];

%punti di controllo
nP = length(airfoil)-1;
P = zeros(nP,2);
P(:,1) = mean([x(1:end-1,1)';x(2:end,1)']);
P(:,2) = mean([x(1:end-1,2)';x(2:end,2)']);

%versori normali e tangenti
li = zeros(nP,1);
t = zeros(nP,2);
n = zeros(nP,2);
for i=1:nP
    t(i,:) = x(i+1,:)-x(i,:);
    li(i) = norm(t(i,:));
    t(i,:) = t(i,:)/li(i);
    n(i,:) = [-t(i,2), t(i,1)];
end

%matrice di influenza aerodinamica
A = zeros(nP+1,nP+1);
b = zeros(nP+1,1);

%condizione di impermeabilità
vVAn = zeros(nP,nP);
vVAt = zeros(nP,nP);
vVBn = zeros(nP,nP);
vVBt = zeros(nP,nP);
for i=1:nP
    for j=1:nP
        %effetto del pannello j sul punto di controllo i
        %[cVAt,cVAn,cVBt,cVBn] = induced_panel_vortex(x(j,:),x(j+1,:),P(i,:));
        %[cVAt,cVAn,cVBt,cVBn] = induced_panel_vortex(xP1,   xP2,     x)
        l = li(j);
        r = P(i,:)-P(j,:);
        xREL = dot(r,t(j,:));
        yREL = dot(r,n(j,:));

        tmp1 = 0;
        tmp2 = -pi;
        if i ~= j
            tmp1 = log(yREL^2+(xREL+l/2)^2) - log(yREL^2+(xREL-l/2)^2);
            tmp2 = atan2(xREL+l/2,yREL) - atan2(xREL-l/2,yREL);
        end
        cVAt = -tmp2/(2*pi) + ((l+2*xREL)*tmp2 -yREL*tmp1)/(4*pi*l);
        cVAn = tmp1/(4*pi) - (0.25*(l+2*xREL)*tmp1 + yREL*tmp2 -l)/(2*pi*l);
        cVBt = -((l+2*xREL)*tmp2 -yREL*tmp1)/(4*pi*l);
        cVBn = (0.25*(l+2*xREL)*tmp1 + yREL*tmp2 -l)/(2*pi*l);
        
        %pannello vortice
        A(i,j) = A(i,j) + dot(cVAt*t(j,:)+cVAn*n(j,:), n(i,:));
        A(i,j+1) = A(i,j+1) + dot(cVBt*t(j,:)+cVBn*n(j,:), n(i,:));
        
        %tengo i coefficienti in memoria
        vVAt(i,j) = cVAt;
        vVAn(i,j) = cVAn;
        vVBt(i,j) = cVBt;
        vVBn(i,j) = cVBn;
    end
    b(i) = -dot(vinf, n(i,:));
end

%condizione di Kutta
A(nP+1,1) = 1;
A(nP+1,nP+1) = 1;
%error;

%calcolo cp
lambda = linsolve(A,b);
cp = zeros(nP+1,2);
%cp(:,1) = x(:,1);
cp(:,1) = airfoil(:,1);
cp(:,2) = 1-lambda(:).^2;

%coefficienti aerodinamici
Cl = 0;
Cd = 0;
Cm = [0, 0, 0];
x0 = Rv*[0.25; 0];
Ftot = zeros(nP,2);
for i=1:nP
    %Ftot(i,:) = 0.5*(cp(i+1,2)+cp(i,2))*li(i)*n(i,:);
    Ftot(i,:) = li(i)*n(i,:)*(3-(lambda(i+1)^2)-(lambda(i)^2)-(lambda(i)*lambda(i+1)))/3;
	Cl = Cl + dot(Ftot(i,:),[0, 1]);
	Cd = Cd + dot(Ftot(i,:),[1, 0]);
	Cm = Cm + cross([P(i,1),P(i,2),0]-[x0(1),x0(2),0], [Ftot(i,1),Ftot(i,2),0]);
end
Cm = -Cm(3);

time = etime(clock,T0);
end


