
disp(' ')
disp(' Rotational flow about the z axis ')

x = linspace(-2,2,6);
y = x;
z = linspace(-2,2,6);
[X,Y,Z] = meshgrid(x,y,z);
R = sqrt(X.^2 + Y.^2)+eps;
U = -Y;
V = X;
W = 0*X;
quiver3(X,Y,Z,U,V,W)
axis image
view(30, 55)
rotate3d on

pause

disp(' ')
disp(' Gravitational vector field  ')

clear x y z X Y Z U V W R
x = linspace(-1,1,6);
[X,Y,Z] = meshgrid(x,x,x);

R = sqrt(X.^2 +Y.^2 +Z.^2)+eps;
U = -(X./R.^3).*(R > .45);
V = -(Y./R.^3).*(R > .45);
W = -(Z./R.^3).*(R >.45);
quiver3(X,Y,Z,U,V,W);
axis image
rotate3d on

pause

disp(' ')
disp(' Electric field for two positive point charges at (x,y,z) = (0,0,0) ')
disp(' and at (x,y,z) = (2,0,0). Charge at (0,0,0) is twice that at (2,0,0) ')
clear x y z X Y Z U V  W R
x = linspace(0,2,6);
y = linspace(-1,1,6);
z = y;
[X,Y,Z] = meshgrid(x,y,z);
R0 = sqrt(X.^2 + Y.^2 + Z.^2)+eps;
R1 = sqrt((X-2).^2 + Y.^2 + Z.^2)+eps;
U0 = (X./R0.^3).*(R0 > .45);
V0 = (Y./R0.^3).*(R0 > .45);
W0 = (Z./R0.^3).*(R0 > .45);

U1 = ((X-2)./R1.^3).*(R1 > .45);
V1 = (Y./R1.^3).*(R1 > .45);
W1 = (Z./R1.^3).*(R1 > .45);

U = U0 + .5*U1;
V = V0 + .5*V1;
W = W0 + .5*W1;

quiver3(X,Y,Z,U,V,W)
axis image
rotate3d on
pause

close


