%% Escoamento potencial usando coordenadas polares
%
% Aqui estão exemplificados os exemplos de escoamento potêncial elementares
% e suas combinações. Esses exemplos fazem parte das notas de aula referente
% ao capítulo 6 do livro.
%

r = 0.02:0.001:5;
theta = 0:2*pi/100:2*pi;

[r,theta]=meshgrid(r,theta);

%% Escoamento uniforme

U0 = 1;

psi_uf = @(r,theta,U0) U0*r.*sin(theta);
phi_uf = @(r,theta,U0) U0*r.*cos(theta);

Vr_uf  = @(r,theta,U0) U0*cos(theta);
Vt_uf  = @(r,theta,U0) -U0*sin(theta);

contour(r.*cos(theta),r.*sin(theta),psi_uf(r,theta,U0),30,'b'); hold on
contour(r.*cos(theta),r.*sin(theta),phi_uf(r,theta,U0),(-3:0.2:3),'r'); hold off
axis equal

legend('Corrente','Equipotênciais')

%% Fonte Sorvedouro
%

q     = 1;

psi_f = @(r,theta,q)  q/2/pi.*theta;
phi_f = @(r,theta,q) -q/2/pi.*log(r);

Vr_f  = @(r,theta,q) q/2/pi./r;
Vth_f = @(r,theta,q) 0;

contour(r.*cos(theta),r.*sin(theta),psi_f(r,theta,q),30,'k'); hold on
contour(r.*cos(theta),r.*sin(theta),phi_f(r,theta,q),(-1:0.02:1),'k--'); hold off
axis equal
legend('Corrente','Equipotênciais')

%% Vortice

k     = 1;

phi_v = @(r,theta,k)  k/2/pi.*theta;
psi_v = @(r,theta,k) -k/2/pi.*log(r);

Vth_v  = @(r,theta,k) k/2/pi./r;
Vr_v   = @(r,theta,k) 0;

contour(r.*cos(theta),r.*sin(theta),psi_v(r,theta,k),(-1:0.02:1),'k'); hold on
contour(r.*cos(theta),r.*sin(theta),phi_v(r,theta,k),30,'k:'); hold off
axis equal
legend('Corrente','Equipotênciais')

%% Dipolo

lambda = 1;
psi_d  = @(r,theta,lambda) -lambda/2/pi*sin(theta)./r;
phi_d  = @(r,theta,lambda)  lambda*cos(theta)./r/2/pi;

n=20;
v = 1 *(-n:n)/n/2;
contour(r.*cos(theta),r.*sin(theta),psi_d(r,theta,q),v,'k'); hold on
contour(r.*cos(theta),r.*sin(theta),phi_d(r,theta,q),v,'k:'); hold off
axis equal
legend('Corrente','Equipotênciais')

%% Cilindro

a = sqrt(lambda/U0);
Estag = psi_uf(a,0,U0)+psi_d(a,0,lambda);
contour(r.*cos(theta),r.*sin(theta),psi_d(r,theta,lambda)+psi_uf(r,theta,U0),-5:0.25:5,'k'); hold on
plot(cos(0:2*pi/100:2*pi)*a,sin(0:2*pi/100:2*pi)*a,'b')
contour(r.*cos(theta),r.*sin(theta),psi_d(r,theta,lambda)+psi_uf(r,theta,U0),Estag*[1 1],'k'); hold off
axis equal
axis([-3 3 -3 3])


%% Cilindro rotativo
contour(r.*cos(theta),r.*sin(theta),psi_d(r,theta,lambda)+psi_uf(r,theta,U0)+psi_v(r,theta,k),32,'k'); hold on
contour(r.*cos(theta),r.*sin(theta),psi_d(r,theta,lambda)+psi_uf(r,theta,U0)+psi_v(r,theta,k),[1 1]*psi_v(1,0,k),'r');
contour(r.*cos(theta),r.*sin(theta),psi_d(r,theta,lambda)+psi_uf(r,theta,U0)+psi_v(r,theta,k),-1:0.1:1,'r'); hold off

axis equal
axis([-3 3 -3 3])


%% Uniforme mais vortice
contour(r.*cos(theta),r.*sin(theta),psi_uf(r,theta,U0)+psi_v(r,theta,k),30,'k');
%contour(r.*cos(theta),r.*sin(theta),psi_uf(r,theta,U0)+psi_v(r,theta,k),[0 0],'r'); hold off
axis([-1 1 -1 1])
axis equal

x = -1.5; y = 0.6;
r1 = sqrt(x^2 + y^2); theta1 = atan2(y,x);
psi1 = psi_uf(r1,theta1,U0)+psi_v(r1,theta1,k);

x = -1.5; y = -0.6;
r1 = sqrt(x^2 + y^2); theta1 = atan2(y,x);
psi2 = psi_uf(r1,theta1,U0)+psi_v(r1,theta1,k);

contour(r.*cos(theta),r.*sin(theta),psi_uf(r,theta,U0)+psi_v(r,theta,k),30,'k');
axis equal
axis([-3 3 -3 3])

%% Fonte mais Uniforme

v =1;
contour(r.*cos(theta),r.*sin(theta),psi_f(r,theta,q)+psi_uf(r,theta,U0),'k'); hold on
contour(r.*cos(theta),r.*sin(theta),psi_f(r,theta,q)+psi_uf(r,theta,U0),v,'k'); hold off
axis equal
axis([-3 3 -3 3])

%% Sorvedouro mais vortice

v =1;
contour(r.*cos(theta),r.*sin(theta),-psi_f(r,theta,q)+psi_v(r,theta,k),'k'); hold on
contour(r.*cos(theta),r.*sin(theta),-psi_f(r,theta,q)+psi_v(r,theta,k),v,'k'); hold off
axis equal
axis([-3 3 -3 3])