%% Exemplo de fontes simétricas
%
% Este exemplo visa mostras o escoamento resultante da utilização de duas
% fontes posicionandas simetricamente em relação ao eixo x, na posição x=0
% e distantes H do eixo x.

clear all
close all

% Escoamento potencial Fonte
Psi = @(q,r,th) q*th/2/pi;
Vr  = @(q,r,th) q/2/pi./r;



% mudança de variável polar para cartesiana da velocidae
u   = @(Vr,th) Vr.*cos(th);
v   = @(Vr,th) Vr.*sin(th);



%% definição da malha
%
% A malha, pontos onde a função corrente e a velocidades serão calculados,
% é definido no limite -2h + 2h em x e y.
H = 1;

X = -2:0.01:2;

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

%% calculo do escoamento induzido pela ponte 1, acima da abscissa
%
% Para promover este calculo, as coordenadas globais cartesianas são 
% escritas em função das coordenadas locais polares.
q = 1;

r1  = sqrt(x.^2 + (y-H).^2);
th1 = atan2(y-H,x);

Psi1 = Psi(q,r1,th1);
Vr1  = Vr(q,r1,th1);
u1   = u(Vr1,th1);
v1   = v(Vr1,th1);

uu1 = @(x,y) u(Vr(q,sqrt(x.^2 + (y-H).^2),atan2(y-H,x)),atan2(y-H,x));
vv1 = @(x,y) v(Vr(q,sqrt(x.^2 + (y-H).^2),atan2(y-H,x)),atan2(y-H,x));

%% calculo do escoamento induzido pela ponte 2, abaixo da abscissa
%
% Idem à fonte 1
q = 1;

r2  = sqrt(x.^2 + (y+H).^2);
th2 = atan2(y+H,x);

Psi2 = Psi(q,r2,th2);
Vr2  = Vr(q,r2,th2);
u2   = u(Vr2,th2);
v2   = v(Vr2,th2);

uu2 = @(x,y) u(Vr(q,sqrt(x.^2 + (y+H).^2),atan2(y+H,x)),atan2(y+H,x));
vv2 = @(x,y) v(Vr(q,sqrt(x.^2 + (y+H).^2),atan2(y+H,x)),atan2(y+H,x));

%% Escoamento final - Sobreposição de soluções
%
% dfd
Psit = Psi1 + Psi2;

uut = @(x,y) uu1(x,y) + uu2(x,y) +1 ;
vvt = @(x,y) vv1(x,y) + vv2(x,y);

ut   = uut(x,y);
vt   = vvt(x,y);

ind = 1:10:length(X);
figure
quiver(x(ind,ind),y(ind,ind),ut(ind,ind),vt(ind,ind))

rr = 0.05;
i = (0:1/20:1);
xstart = [x(ind,1)' rr*cos(2*pi*i)   rr*cos(2*pi*i)   ];
ystart = [y(ind,1)' rr*sin(2*pi*i)+H rr*sin(2*pi*i)-H ];

h=streamline(x,y,ut,vt,xstart,ystart);
set(h, 'Color', 'black');

axis equal


figure
a=stream2(x,y,ut,vt,xstart,ystart);

s=0.05;

hold on
for i=1:length(a)
    arco=sum(sqrt(sum(diff(a{i}).^2,2)));
    
    ind=1:floor(size(a{i},1)/floor(arco/0.1)):size(a{i},1);
    plot(a{i}(:,1),a{i}(:,2),'k');
    quiver(a{i}(ind,1),...
           a{i}(ind,2),...
          s*uut(a{i}(ind,1),a{i}(ind,2)),...
          s*vvt(a{i}(ind,1),a{i}(ind,2)),0,'r');
end
hold off
axis equal

%% Campo de Pressão
% 
% Aqui a pressão é obtida relativa a pressão atmosférica do escoamento
% longe das fontes. Para isso aplica-se a equação de bernoulli

rho  = 1;
p = 1-(ut.^2+vt.^2) ;

% Lembrando que no centro das fontes existem singularidades, os pontos com
% pressão abaixo de minn são igualados a minn.

minn = -8.;
N = 100;
p(p < minn) = minn;

figure
hold on
contourf(x,y,p,N)
contour(x,y,p,N)
h=streamline(a);
set(h, 'Color', 'black');
hold off
axis equal
axis([-2 2 -2 2])