%
%  w05axi.m - velocity field (djm: 08 feb 06)
%
%  - based on w01flow.m

%  initialize workspace & set constants (try different values of B)
clear;  close all

L = 3.0;  ds = 2*L/45;
A = 1.0;  M = 4*pi/4/pi;

%  set x & y-vectors &  grid matrices ("help meshgrid")
r = -L:ds:+L;  z = -L:ds:+L;
[rr,zz] = meshgrid(r,z);
mag = sqrt(rr.^2 + zz.^2);

%  calculate psi
psi = (A/2)*rr.^2 + M*(1 - zz./mag);

%  calculate velocities (U = -psi_z/r, W = psi_r/r)
UU =     M*rr./(mag.^3);
WW = A + M*zz./(mag.^3);

vel = sqrt(UU.^2 + WW.^2);

%  plot vector field ("help quiver") outside unit circle ("help gt")
clf;  hold on

%  calculate pressure/density
pressure = (A^2 - (UU.^2 + WW.^2))/2;
pcolor(rr,zz,pressure);  shading interp
colormap(winter)
caxis([-0.8 0.5]);  colorbar 

mask = (sqrt(rr.^2 + zz.^2)>0.75);
Um = UU.*mask;  Wm = WW.*mask;

quiver(rr,zz, Um, Wm,0.5,'k')
quiver(rr,zz,-Um,-Wm,0.5,'k.')
axis([-L +L -L +L]);  axis square;

%  plot level curves of psi (two versions, try both)
contour(rr,zz,psi,[-4:0.5:4]+2*M,'w--','linewidth',2);

%  plot one special contour
contour(rr,zz,psi,[1 1]*2*M,'w','linewidth',2);
plot([0 0],[-L L],'w','linewidth',2);

%  label plot
title(['\bf an axisymmetric flow - with pressure'])
xlabel('\bf r-axis')
ylabel('\bf z-axis')

text(-0.13,-1,'H','fontsize',24)
text(-0.12,0,'L','fontsize',24)