r = 0.02:0.02:3;
theta = 0:2*pi/30: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);

%% 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;

%% Uniforme mais vortice

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);

h = contour(r.*cos(theta),r.*sin(theta),psi_uf(r,theta,U0)+psi_v(r,theta,k),[psi1 psi2]);
ind = find(h(1,:) == psi1);
len = h(2,ind);
hh = h(:,ind + (1:len-1));
index = find(hh(1,:) > -1.5 & hh(1,:) < 1.5);
plot(hh(1,index),hh(2,index),'Linewidth',2); hold on

ind = find(h(1,:) == psi2,1,'last');
len = h(2,ind);
hh = h(:,ind + (1:len-1));
index = find(hh(1,:) > -1.5 & hh(1,:) < 1.5);
plot(hh(1,index),hh(2,index),'Linewidth',2);

line([-1.5 -1.5],[-0.6 0.6],'Linewidth',2);
line([1.5 1.5],[-0.6 0.6],'Linewidth',2);

contour(r.*cos(theta),r.*sin(theta),psi_uf(r,theta,U0)+psi_v(r,theta,k),30,'k');


hold off
axis([-2 2 -1 1])
axis equal