%% Perfil simétrico
clear all
close all

%% definição do problema
U0 = 1;

alpha = 5*pi/180;

C = 1;
N = 50;

x = 0:1/N:1;


% y = x*0; gammaex = @(x) 4*U0*alpha*(1-x)./sqrt(1-(1-2*x).^2); % Placa Plana
% t = 0.01; R = (.5^2 + t^2)/2/t; y = -R + t + sqrt(R^2 - (x-0.5).^2); gammaex = 0;% arco Ciscunferencia
 t = 0.05; y = -4*t*x.*(x-1); gammaex = @(x) 2*U0*(2*alpha*(1-x)./sqrt(1-(1-2*x).^2) +4*t*sqrt(1-(1-2*x).^2))


xy1 = [x(1:end-1)' y(1:end-1)'];
xy2 = [x(2:end)' y(2:end)'];

figure
plot(x,y,'o-')
axisd equal

xf= x; yf = y;


%% Solução do sistema
[gamma,xy0,xyc,c] = LumpedVortex(xy1,xy2,alpha,U0);

figure
plot(xy0(:,1),gamma,'o',xy0(:,1),gammaex(xy0(:,1)).*c,'-k')

%% Plot do campo de velocidades
% x = -1: 0.02:2;
% y = -0.5:0.02:0.5;

beta  = 1.1;
n1 = 30;
eta   = 0:1/(n1-1):1;
H  = 1;
x1 = H*(beta+1 -(beta-1)*((beta+1)/(beta-1)).^(1-eta))...
    ./(((beta+1)/(beta-1)).^(1-eta)+1);

n2 = floor(1/min(diff(x1)));
x2 = 0:1/(n2-1):1;

x = unique([-x1(end:-1:1) x2 x1+1]);

y = unique([ -x1(end:-1:1) x1 ]);



[x,y] = meshgrid(x,y);

u = x*0;
v = u;
for i=1:numel(x)
    V = vort2D(xy0,[x(i) y(i)],gamma);
    
    u(i) = V(1)+U0*cos(alpha);
    v(i) = V(2)+U0*sin(alpha);
end

figure
quiver(x,y,u,v); hold on
plot(xy0(:,1),xy0(:,2),'o');
hold off
axis equal

figure
hold on
streamline(x,y,u,v,-ones(size(-1:0.1:1)),-1:0.1:1)
plot(xy0(:,1),xy0(:,2),'o',xyc(:,1),xyc(:,2),'*');
hold off

%% Plot do campo de pressão

p = 1-(u.^2 + v.^2)/U0.^2;
p(p<-3) = -3;

figure
hold on
contourf(x,y,p,100)
contour(x,y,p,100)
plot(xf,yf,'k','LineWidth',2)
hold off