function[cond] = condNumSolvent(coeffs,S,C,R,nodes)
% This function approximates the condition number for a solvent S

% Input:  coeffs  -  cell array containing coefficients of matrix polynomial
%              S  -  nxn matrix solvent S
%              C  -  center of circle
%              R  -  radius of circle
%          nodes  -  number of nodes of the trapezoid rule to approximate the contour integral
% Output:    cond -  condition number of the solvent S

len = length(coeffs)-1; BB = [];

% C is the center of the circle and R the radius. Last parameter is the number of nodes
[M,B] = condNumIntegrals(coeffs,S,C,R,nodes);

for i=len+1:-1:1
    BB = [BB norm(coeffs{i},'fro')*B(:,:,i)];
end

if(det(M)<eps)
    cond = inf;
else
    cond = norm( M\BB )/norm(S,'fro');
end
end

function[M,B] = condNumIntegrals(poly,S,C,R,nodes)
% This function computes the INTEGRALS FOR THE CONDITION NUMBER,
% approximating the contour integral with the trapezoid sum and using a
% circle as contour.

n = length(S); len = length(poly)-1; M = zeros(n*n,n*n);
for jj = 1:len+1, B(:,:,jj) = zeros(n*n,n*n); end

for j=0:nodes-1
    expon = exp(1i*j*2*pi/nodes); phi = (C+R*expon);
    % Compute T(phi)
    T = zeros(n,n);
    for m = 0:len, T = T + phi^m*poly{m+1}; end
    M = M + expon*kron( inv(phi*eye(n,n)-S).', T/(phi*eye(n,n)-S) );
    
    for ii = 1:len+1
        B(:,:,ii) = B(:,:,ii) + expon*kron(phi^(ii-1)*(inv(phi*eye(n,n)-S)).', eye(n,n));
    end
end
M = (R/nodes)*M; B = (R/nodes)*B;
end