
clear ; close;
disp(' ')
disp('  Rotational Flow  ')
x = linspace(-2,2, 17);
y = x;
[X,Y]= meshgrid(x,y);

R = (X.^2 + Y.^2).^.25 +eps;
U = -Y./R;
V = X./R;
quiver(X,Y,U,V)
axis image

pause
disp(' ')
disp(' Radial vector field ')
clear U V
U = X./R;
V = Y./R;
quiver(X,Y,U,V)
axis image

pause
disp(' ')
disp(' Velocity field for flow around a cylinder  ')
a = .75;
clear U V R;
denom = (X.^2 + Y.^2).^2 + eps; 
U = 1+ a^2*(Y.^2 - X.^2)./denom;
U = U.*(1-(denom < a^4));
V = -2*a^2*X.*Y./denom;
V = V.*(1-(denom < a^4));
quiver(X,Y,U,V)
t = linspace(0, 2*pi,101);
hold on
fill(a*cos(t), a*sin(t), 'c')
axis image
hold off
pause

disp(' ')
disp(' Electric field for positive line charges at (x,y) = (-1,0) and at ')
disp(' (x,y) = (1,0) ')

clear x y X Y U V R

x = linspace(-2,2,11);
y = x;
[X,Y] = meshgrid(x,y);
R1 = sqrt((X+1).^2 + Y.^2);
R2 = sqrt((X-1).^2 + Y.^2);
U1 = ((X+1)./R1.^2).*(R1 > .3);
U2 = ((X-1)./R2.^2).*(R2 > .3);
U = U1+U2;
V1 = (Y./R1.^2).* (R1> .3);
V2 = (Y./R2.^2).*(R2 > .3);
V = V1 + V2;
quiver(X,Y,U,V)
axis image
pause
close