%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% Find the pressure for the airfoil found by airfoil.m %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set the density if it has not yet been set
if exist('rho') ~= 1
  rho=1
end

% set the pressure at infinity if it has not yet been set
if exist('pinf') ~= 1
  p_inf=0
end

% find the z-derivative of F
dF_dz=U*(exp(-i*alpha)-r0^2*exp(i*alpha)./z.^2)+i*Gamma/(2*pi)./z;

% find the z-derivative of zeta
dzeta_dz=1.-1./(z-eps).^2;

% find the zeta derivative of F
dF_dzeta=dF_dz./dzeta_dz;

% the magnitude of the velocity is the magnitude of dF_dzeta
vel=abs(dF_dzeta);

% find the pressure
p=p_inf+0.5*rho*U^2-0.5*rho*vel.^2;

% plot it in blue
hold on
contour(xi,eta,p,floor(min(min(p))):.25:ceil(max(max(p))),'b')
axis equal
hold off
