function tab = run_example_OdOP_IETCTA_ex3()

points = 41;

val_omega = linspace(0,4,points);


tab_max_beta = draw_figure(val_omega);
calculate_complexities_omega(val_omega(21),tab_max_beta(:,21));

tab = [];

%________________________________________________________________________________________________________
function tab_max_beta = draw_figure(val_omega)

tab_max_beta = [];
for method=1:2
    for kappa=1:4
        t=clock;
        param.kappa=kappa;
        ini_upper = 1;
        y = [];
        for i=1:size(val_omega,2)
            max_beta = bissection(method,param,val_omega(i),ini_upper,0.05);
            fprintf('omega=%.3f -> %.3f\n ',val_omega(i),max_beta);
            y = [y max_beta];
        end
        tab_max_beta = [tab_max_beta;y];
        save data_ex3.mat val_omega tab_max_beta;
        t=etime(clock,t);
        fprintf('finished method=%d -> kap = %d \t %.2f\n ',method,kappa,t);
    end
    yalmip('clear');
    pack;
end

f=figure(1);
plot(val_omega,tab_max_beta(1:4,:));
axis([0 4 0 1]);
xlabel  ome;
ylabel  bet;
pos = get(f,'Position');
pos(3) = 840;
pos(4) = 630;
set(f,'Position',pos);

f=figure(2);
plot(val_omega,tab_max_beta(5:8,:));
axis([0 4 0 1]);
xlabel  ome;
ylabel  bet;
pos = get(f,'Position');
pos(3) = 840;
pos(4) = 630;
set(f,'Position',pos);



%________________________________________________________________________________________________________
function calculate_complexities_omega(omega,val_beta)
%Complexities for omega=2

cont=1;
times = 10;
fprintf('Method    K \t L \t Time\n');
for method=1:2
    for kappa=1:4
        param.kappa = kappa;
        [Ai,dot_alphas] = system_matrices(val_beta(cont),omega);
        cont = cont + 1;
        t=0;
        for j=1:times
            Info = run_method(method,Ai,dot_alphas,param);
            t = t + Info.cpusec;
        end
        t = t/times;
        fprintf('C%d(k=%d)   %d \t %d \t %.2f - %d\n',method,kappa,Info.K,Info.L,t,Info.feas);
    end
end

%Method    K 	 L 	 Time
%C1(k=1)   17 	 18 	 0.06 - 1
%C1(k=2)   62 	 30 	 0.06 - 1
%C1(k=3)   135 	 42 	 0.10 - 1
%C1(k=4)   236 	 54 	 0.15 - 1
%C2(k=1)   33 	 30 	 0.07 - 1
%C2(k=2)   126 	 60 	 0.12 - 1
%C2(k=3)   279 	 84 	 0.27 - 1
%C2(k=4)   492 	 108 	 0.66 - 1




%________________________________________________________________________________________________________
function max_beta = bissection(method,param,omega,upper,lower)

% bissection precision
maxgap = 0.001;

ub = upper;
lb = lower;

oncefeas = 0;
gap = ub - lb;
max_beta = 0;
while(gap > maxgap)

    beta = (ub + lb)/2;

    [Ai,H] = system_matrices(beta,omega);
    Info = run_method(method,Ai,H,param);

    if Info.feas == 1
        lb = beta;
        oncefeas = 1;
        max_beta = beta;
        %fprintf('%.3f ',beta);
    else
        ub = beta;
    end
    gap = ub-lb;
end

%________________________________________________________________________________________________________
function Info = run_method(method,Ai,H,param)


switch method
    case 1
        Info = stabLPV_OdOP_IETCTA_C1b_yal(Ai,param.kappa,H);
    case 2
        Info = stabLPV_OdOP_IETCTA_C3b_yal(Ai,param.kappa,H);

end

%________________________________________________________________________________________________________
function [Ai,H] = system_matrices(beta,omega)

A1 = [0 1;
    -(1 + 3*beta^2) -0.05];
A2 = [0 1;
    -(1 - beta^2) -0.05];
A3 = A2;
Ai = [A1 A2 A3];

q=omega/4;

H = [q -3*q q;
    3*q -q -q;
    -4*q 4*q 0];
