%   Creeping flow in a wedge with specified velocity on one side and zero shear on other
%   The zero shear is at theta=0, the specified velocity is at theta=-alpha.
%   Domain is 0 < r < 1; -alpha < theta < 0
    clear all
    alpha=input(' alpha = ');
    [theta,r]=meshgrid((-1:0.05:0)*alpha, 1.e-6:0.05:1+1.e-6);
    [x,y]=pol2cart(theta,r);
    denom=alpha-sin(alpha)*cos(alpha);
    SF=r.*(alpha*cos(alpha).*sin(theta)-sin(alpha).*theta.*cos(theta))./denom;
    figure(1)
    clf reset
    contour(x,y,SF,20,'r')
    hold on
    plot([0 1],[0 0],'k',[0 x(end,1)],[0 y(end,1)],'k',x(end,:),y(end,:),'k')
    title(['Streamlines for Flow in Wedge with Zero Shear at Theta = 0; alpha = ',num2str(alpha)])
    xlabel('x')
    ylabel('y')
    axis equal
%   Calculate velocity
    figure(2)
    clf reset
    vx=zeros(size(theta));
    vy=zeros(size(theta));
    vr=(alpha*cos(alpha).*cos(theta)-sin(alpha).*(cos(theta)-theta.*sin(theta)))./denom;
    vth=-(alpha*cos(alpha).*sin(theta)-sin(alpha).*theta.*cos(theta))/denom;
    for j=1:length(theta(1,:))
        L=[cos(theta(1,j)), -sin(theta(1,j));sin(theta(1,j)), cos(theta(1,j))];
        vcart=L*[vr(1,j);vth(1,j)];
        vx(:,j)=vcart(1);
        vy(:,j)=vcart(2);
    end
    quiver(x,y,vx,vy)
    title(['Velocity for Flow in Wedge with Zero Shear at Theta = 0; alpha = ',num2str(alpha)])
    xlabel('x')
    ylabel('y')
    axis equal
%   Calculate pressure
    figure(3)
    clf reset
    P=-2*sin(alpha).*cos(theta)./(r.*denom);
    P(1,:)=0;
    continv=linspace(-100,100,200);
    contour(x,y,P,continv,'r')
    hold on
    contour(x,y,P,[0 0],'g')
    plot([0 1],[0 0],'k',[0 x(end,1)],[0 y(end,1)],'k',x(end,:),y(end,:),'k')
    title(['Pressure for Flow in Wedge with Zero Shear at Theta = 0; alpha = ',num2str(alpha)])
    xlabel('x')
    ylabel('y')
    axis equal
