% function [delta,lowerbnd,upperbnd] = wcperf(sysg,blk,pertsize)
function [delta,lowerbnd,upperbnd] = wcperf(sysg,blk,pertsize,startalpha)

    if nargin == 0
        disp('usage: [delta,lowerbound,upperbound] = wcperf(sysg,blk,pertsize)')
        return
    elseif nargin == 3
        startalpha = [.3 1 1/.3];
    end

 perctol = 0.001;
 itermax = 10;
    

 [nblk,dum] = size(blk);
 [mtype,mrows,mcols,mnum] = minfo(sysg);
 repblks = find(blk(:,2)==0);
 ablk = abs(blk)
 if length(repblks)>0
    ablk(repblks,2) = ablk(repblks,1);
 end
 rdimblk = sum(ablk(:,1));
 cdimblk = sum(ablk(:,2));
 if rdimblk<mcols & cdimblk<mrows
    blk = [blk;[mcols-rdimblk mrows-cdimblk]];
    nblk = nblk + 1;
 elseif rdimblk~=mcols | cdimblk~=mrows
    error('dimension error in Uncertainty Structure')
    return
 end
 
 
 cperf = blk(nblk,1);
 rperf = blk(nblk,2);
 if blk(nblk,2) == 0
   rperf = blk(nblk,1);
 end

% FIRST POINT
 nominal = pkvnorm(sel(sysg,mrows-rperf+1:mrows,mcols-cperf+1:mcols));
% LAST POINT
 if any(blk(1:nblk-1,1)>0)
    [bnds,dvec,sens,pvec] = mu(sel(sysg,1:mrows-rperf,1:mcols-cperf),blk(1:nblk-1,:),'csw');
    xdestab = 1/pkvnorm(sel(bnds,1,2)); 
 else
    xdestab = inf;
    disp('Warning: Uncertainty Structure is All Real - There may be problems')
 end

    if pertsize>=xdestab
        disp('Perturbation Size is destabilizing')
        delta = dypert(pvec,blk(1:nblk-1,:),bnds,[1:nblk-1]);
        return
    end

    xl = 0;
    yl = nominal;
    xu = 0;
    yu = nominal;

    xmax = xdestab;;

    [out,delta] = onealpha(sysg,blk,startalpha,pertsize,perctol);
    xl = [xl;out(:,1)]
    yl = [yl;out(:,2)] 
    xu = [xu;out(:,3)];
    yu = [yu;out(:,4)]; 
 
 if isempty(delta)
    dist = inf;
    bestdelta = []; 
 else
    ndelta = hinfnorm(delta);
    ndelta = ndelta(1);
    bestdelta = delta;
    dist = abs(ndelta-pertsize);
 end
 cnt = 0;

 while cnt<itermax & isempty(delta)
    p = polyfit(xl,yl,min([4 length(xl)-2]));
    newalpha = 1/pertsize/polyval(p,pertsize);
    [out,delta] = onealpha(sysg,blk,newalpha,pertsize,perctol);
    if ~isempty(delta)
        ndelta = hinfnorm(delta);
        ndelta = ndelta(1);
        if abs(ndelta-pertsize)<dist
            bestdelta = delta;
            dist = abs(ndelta-pertsize);
        end
    end
    xl = [xl;out(:,1)];
    yl = [yl;out(:,2)]; 
    xu = [xu;out(:,3)];
    yu = [yu;out(:,4)]; 
    cnt = cnt+1;
end
delta = bestdelta;
[xls,index] = sort(xl);
lowerbnd = vpck(yl(index),xls);
[xus,index] = sort(xu);
upperbnd = vpck(yu(index),xus);
