% Matlab file to visualize the
% Velocity profile for the co-rotating vortex pair

clear
clc
H=5;
%[x] = meshgrid(0:H/10:H);
%x = [-H:H/20:H];
%y = 0; %[-10:0.1:10];
[x,y] = meshgrid(-H:H/40:H);


i=sqrt(-1);
z = x+i*y;
zC = x-i*y;
r0=3;
b = r0;
Gamma = 4;

f = (-i)*Gamma/pi*(z./((z).^2 - b^2));
% left vortex
f1 = Gamma/pi/i*1./(z+b);
%right vortex
f2 = Gamma/pi/i*1./(z-b);
rC = r0/3;

xmax = -(r0+rC);
umax = Gamma/pi/i*1/(xmax+b);

size(z)
size(zC)

for a = 1:size(x,2)
    for b = 1:size(y,2)
        if( (z(a,b) - r0)*(zC(a,b)-r0) < rC^2 )  %right vortex
            f2( a,b ) = -umax*zC(a,b)/rC+umax*r0/rC;
        elseif( (z(a,b) + r0)*(zC(a,b)+r0) < rC^2 )
            f1( a,b ) = -umax*zC(a,b)/rC-umax*r0/rC ; % left vortex
        end
    end
end

f=f1+f2;

imag(f);
subplot(2,2,1);
contour(x,y,real(f),40)
subplot(2,2,2);
contour(x,y,-imag(f),40)
subplot(2,2,3);
contour(x,y,sqrt(real(f).^2+imag(f).^2),40)


rC = 0.3*r0;
%Schwartzkopff
%f = Gamma/2/pi/x.^2*(1-e^(-alpha*(x/rC).^2))
%plot(x, f, 'o')

