% after velocity computation

Ix=0.5*real(pt*conj(ux));
Iy=0.5*real(pt*conj(uy));

for k1=1:3:100
   l1=1+(k1-1)/3;
   for k2=1:3:100
      l2=1+(k2-1)/3;
      amp=sqrt(Ix(k1,k2)^2+Iy(k1,k2)^2);
      ux3(l1,l2)=sqrt(sqrt(amp))*Ix(k1,k2)/amp;
      uy3(l1,l2)=sqrt(sqrt(amp))*Iy(k1,k2)/amp;
   end
   x2(l1)=x(k1);
   y2(l1)=x(k1);
end

figure(1)
H=quiver(x2,y2,ux3,uy3);
axis equal
axis([0 2 0 2])
xlabel('x');
ylabel('y');
title('Intensity')
   
% to store the images for web posting (bitmaps)
filnme='spher';
eval(['print -dbitmap ' filnme ' -f']);
