
%% ----------------------------------------------------------------
clf;
rotate3d on;
hold on;

%% ----------------------------------------------------------------
%n = 4;
n = 20;
%n = 100;

xlo = -2.0;
xhi =  2.0;
nx  =  n;
dx  = (xhi - xlo) / nx;

ylo = -2.0;
yhi =  2.0;
ny  =  n;
dy  = (yhi - ylo) / ny;

%% ----------------------------------------------------------------
for x = xlo:dx:xhi
	y = [ylo:dy:yhi];
	[u,v] = beanfunc(x,y);
	plot(u,v,'b');
end

for y = ylo:dy:yhi
	x = [xlo:dx:xhi];
	[u,v] = beanfunc(x,y);
	plot(u,v,'r');
end

x = xlo + (xhi-xlo)/2;
y = [(ylo+yhi)/2:dy:yhi];
[u,v] = beanfunc(x,y);
plot(u,v,'b', 'LineWidth', 3);

y = ylo + (yhi-ylo)/2;
x = [(xlo+xhi)/2:dx:xhi];
[u,v] = beanfunc(x,y);
plot(u,v,'r', 'LineWidth', 3);

t = [0:0.1:6.29];
for r = 1:.2:3
	x = r*cos(t);
	y = r*sin(t);
	[u,v] = beanfunc(x,y);
	plot(u,v,'g');
end

axis([xlo xhi ylo yhi]);

%% ----------------------------------------------------------------
figure;
[x,y]=meshgrid(xlo:dx:xhi, ylo:dx:yhi);
[u,v]=beanfunc(x,y);
quiver(x,y,u,v);
