%   Plot iso % velocity contours in boundary layer
%   fp is f', the normalized velocity along surface
%   Data in betad.dat and solnd.dat
    clear all
    clf reset
    load betad.dat;
    load solnd.dat;
    velval=0.1:0.1:0.9;
    x=linspace(0.02,1,50);
    y=linspace(0.02,4,50);
    [X,Y]=meshgrid(x,y);
    SF=zeros(size(X));
    cont=linspace(0,6,60);
    for ibeta=1:length(betad)
        beta=betad(ibeta);
        angle=90*(2-beta);
        eta=solnd(:,3*ibeta-2);
        f=solnd(:,3*ibeta-1);
        fp=solnd(:,3*ibeta);
%       Iso % velocity contours
        for jfp=2:length(fp)    % Condition fp to be monotonic
            if fp(jfp)<=fp(jfp-1)
                fp(jfp)=fp(jfp-1)+1.e-6;
            end
        end
        etaval=interp1(fp,eta,velval);
        figure(ibeta)
        clf reset
        for ivel=1:length(velval)
            yval=sqrt(2-beta).*etaval(ivel).*x.^((1-beta)/(2-beta));
            plot(x,yval,'-k')
            hold on
            xlabel('x')
            ylabel('y')
            title(['Iso % velocity Contours, angle = ',num2str(angle), ' deg'])
            axis([0 1 0 4])
        end
 %      Iso streamfunction contours
 %      For each value of (x,y), compute eta, interpolate for f, and compute SF
        for ix=1:length(x)
            xfac1=x(ix)^(1/(2-beta));
            xfac2=x(ix)^((1-beta)/(2-beta));
            etaval=y./(sqrt(2-beta)*xfac2);
            SF(:,ix)=xfac1.*interp1(eta,f,etaval');
        end
        figure(10+ibeta)
        clf reset
        contour(X,Y,SF,cont,'-k')
        xlabel('x')
        ylabel('y')
        title(['Streamlines, angle = ',num2str(angle), ' deg'])
    end