% compute extremal interval matrices based on component variations
% only supports buck converter
function [sys] = matrices_extrema(R,L,C)
    idx = 0;
    
    % component tolerances
    Ctol = 0.05;
    Cv(1) = C * (1 - Ctol);
    Cv(2) = C * (1 + Ctol);

    Ltol = 0.05;
    Lv(1) = L * (1 - Ltol);
    Lv(2) = L * (1 + Ltol);

    Rtol = 0.05;
    Rv(1) = R * (1 - Rtol);
    Rv(2) = R * (1 + Rtol);

    [a,b] = intervalTimes(Rv(1),Rv(2),Cv(1),Cv(2));

    [a01l,a01u]=intervalDiv(1,1,Lv(1),Lv(2));
    [a10l,a10u]=intervalDiv(1,1,Cv(1),Cv(2));
    [a11l,a11u]=intervalDiv( 1, 1, a, b);
    
    for li = 1 : 2
        for ci = 1 : 2
            for ri = 1 : 2
                idx = idx + 1;
                sys(idx).Ao = [0, -(1/Lv(li) ); (1/Cv(ci)), -(1/(Rv(ri) *Cv(ci) ))];
                sys(idx).Bo = [0; 0];
                
                sys(idx).Ac = [0, -(1/Lv(li) ); (1/Cv(ci)), -(1/(Rv(ri) *Cv(ci) ))];
                sys(idx).Bc = [(1/Lv(li)); 0];
            end
        end
    end

    reig = [];
    for ii = 1 : idx
        reig = [reig; eig(sys(ii).Ac), ii*ones(2,1)];
    end

    % compute entrywise min and max
    Amin = inf*ones(2,2);
    Amax = -inf*ones(2,2);
    Bmin = inf*ones(2,1);
    Bmax = -inf*ones(2,1);
    for i = 1 : idx
        for j = 1 : 4
            Amin(j) = min( Amin(j), sys(i).Ac(j));
            Amax(j) = max( Amax(j), sys(i).Ac(j));
        end

        for j = 1 : 2
            Bmin(j) = min( Bmin(j), sys(i).Bc(j));
            Bmax(j) = max( Bmax(j), sys(i).Bc(j));
        end
    end

    ridx = 0;
    largest = [];
    smallest = [];
    relation_gt = zeros(idx,idx); % matrix representing the greater than relation
    relation_lt = zeros(idx,idx); % matrix representing the lesser than relation
    relation_eq = zeros(idx,idx); % matrix representing the equality relation (should be the diagonal)
    ntype = 2; % norm type

    for ii = 1 : idx
        sys(ii).ltc = 0; % larger than count
        sys(ii).gtc = 0; % larger than count
        for jj = 1 : idx
            % matrices are equal
            if prod(size(sys(ii).Ac)) == msum(sys(ii).Ac == sys(jj).Ac) || norm(sys(ii).Ac) == norm(sys(jj).Ac)
                [num2str(ii), ' equal to ', num2str(jj)]
                relation_eq(ii,jj) = 1;
            %elseif size(sys(ii).Ac,1) == sum(eig(sys(ii).Ac) > eig(sys(jj).Ac))
            % try matrix norm for comparison
            elseif norm(sys(ii).Ac,ntype) < norm(sys(jj).Ac,ntype)
                [num2str(ii), ' larger than ', num2str(jj)]
                relation_lt(ii,jj) = 1;
                sys(ii).ltc = sys(ii).ltc + 1;
            elseif prod(size(sys(ii).Ac)) == msum(sys(ii).Ac > sys(jj).Ac)
                msum(sys(ii).Ac > sys(jj).Ac)
                %ii
                %jj
            else
                ridx = ridx + 1;
                rsys(ridx).Ac = sys(ii).Ac;
                relation_gt(ii,jj) = 1;
                sys(ii).gtc = sys(ii).gtc + 1;
                %break;
            end
        end

        % if larger than count is one less than the number of matrices, then it must be the largest matrix
        if sys(ii).ltc == (idx - 1)
            largest = [largest; ii]; % should always be size 1 at most, otherwise bug
        end

        % if larger than count is one less than the number of matrices, then it must be the largest matrix
        if sys(ii).gtc == (idx - 1)
            smallest = [smallest; ii]; % should always be size 1 at most, otherwise bug
        end
    end
end