% Plot velocity field for low Reynolds flow around a sphere
% Based on analysis in Landau and Lifshitz, Section 20

 ell = 1;                             % radius of the sphere
 u = [1;0];                           % velocity field as |x| -> oo
 m = 20;                              % grid density for quiver plot

 x1 = linspace(-4,4,m);               % position, horizontal grid
 x2 = linspace(-2.5,2.5,m);           % position, vertical grid
 v1 = zeros(length(x2),length(x1));   % velocity, horizontal grid
 v2 = zeros(length(x2),length(x1));   % velocity, vertical grid

 I = eye(2);
 for k=1:length(x1), for j=1:length(x2)
    x = [x1(k);x2(j)];
    if norm(x)<ell, vv = [0;0]; 
    else
       vv = -.75*(ell/norm(x))*(I+(x*x')/(norm(x)^2))*u ...
            -.25*((ell/norm(x))^3)*(I-3*x*x'/(norm(x)^2))*u + u;
    end
    v1(j,k) = vv(1);
    v2(j,k) = vv(2);
 end, end


 figure(1), clf
 quiver(x1,x2,v1,v2,'k-');              % quiver plot
 z = ell*exp(linspace(0,2i*pi,500));  
 hold on
 fill(real(z),imag(z),.8*[1 1 1])       % plot sphere
 axis equal
 axis([min(x1) max(x1) min(x2)-.1 max(x2)+.1])