function [fc,lambda,L]=sghm(E,A,B,C,D,tol,opt)
% use S-parameter generalized Hamiltonian method (SGHM) and its symmetric variant 
% to check the passivity of a S-parameter descriptor system
% the transfer function is H(s)=C(sE-A)^(-1)B+D
% tol is a threshold to remove the numerical noise in eigenvalue
% computation
% opt=1: symmetric system, then utilize half-size SGHM test
% opt=0 or if opt is not given, use full-size GHM

% output: fc store the crossover frequency (in Hz)
%         L stores the passivity information for each frequency range

%%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Please cite the following references:
% [1] Z. Zhang and N. Wong, "Passivity check of S-parameter descriptor systems
% via S-parameter generalized Hamiltonian methods," IEEE Trans. Advanced Packaging,
% vol. 33, no. 4, pp. 1034-1042, Nov. 2010
% [2] Z. Zhang and N. Wong, "An extension of the generalized Hamiltonian method
% to S-parameter descriptor systems ”, IEEE/ACM Asia and South Pacific Design 
% Automation Conference, pp.43-47, Taipei, Jan. 2010.


if nargin<7
    opt=0;
end
m=size(C,1); n=size(B,2);
if m~=n
    error('the transfer matrix should be square!\n');
end

if isempty(D)
    D=zeros(m,n);
end
Im=eye(m,m);
rd=rank(D*D.'-Im);
if rd<m % perform equivalent model conversion
    alpha=0.5*rand(1,1); % generate a random number in [0,0.5] for model conversion
    Da=alpha*Im-D;
    E=blkdiag(E,zeros(m,m));
    A=blkdiag(A,Im);
    B=[B;Im]; C=[C, Da];
    D=alpha*Im;
end

fc=[]; a=1/(2*pi);
L.frange=[];
L.passivity=[];

S=D*D.'-Im; R=D.'*D-Im;
if ~opt % use full size SGHM
    Ahat=A-B*(D.')*(S\C);     Rhat=B*(R\B.');
    Phat=C.'*(S\C);     J=[Ahat,-Rhat; Phat, C.'*D*(R\B.')-A.'];
    K=blkdiag(E,E.');
    lambda=eig(J,K);  
    id=find(isinf(lambda));
    lambda(id)=[];  
    N=length(lambda);
    realLambda=real(lambda); 
    imagEig=[];
    for k=1:N %find the purely imaginary eigenvalues
        if (abs(realLambda(k)))<tol
            imagEig=[imagEig; 1i*imag(lambda(k))];
        end
    end
    
   k=1;  % now identify all crossover frequencies (in Hz);   
   while (k<=length(imagEig))
       if abs(imagEig(k)+imagEig(k+1))<tol
           fc=[fc; abs(imagEig(k))*a];
           k=k+2;
       else
           fc=[fc; abs(imagEig(k))*a];
           k=k+1;
       end
   end
    
else % for symmmetric system, use half-size SGHM
    X=A-B*D*(S\C); Y=-B*(R\C);    
    Jh=X+Y;  Kh=E*((Y-X)\E);
    lambda=eig(Jh,Kh);
    id=find(isinf(lambda));
    lambda(id)=[];
    
    N=length(lambda);
    imagLambda=imag(lambda); 
    for k=1:N %find the purely positive-real eigenvalues
        if  ((abs(imagLambda(k)))<=0.01*tol)&&(real(lambda(k))>0)
            fc=[fc;realsqrt(real(lambda(k)))*a];
        end
    end   
end

if isempty(fc);
    L.frange=[0,inf];
    L.passivity{1}=checkBR(E,A,B,C,D,50);
    fprintf('the system is %s in the whole frequency band!!\n\n',L.passivity{1});
else
    fc=sort(fc);
    nfc=length(fc);
    L.frange=zeros(nfc+1,2);
    L.frange(1,:)=[0,fc(1)];
    L.frange(nfc+1,:)=[fc(nfc),inf];
    for k=2:nfc
        L.frange(k,:)=[fc(k-1),fc(k)];
    end
    
    for k=1:(nfc+1)
        L.passivity{k}=checkBR(E,A,B,C,D,L.frange(k,:));
        fprintf('frequency band: [%d, %d],     %s!!\n\n',L.frange(k,1), L.frange(k,2), L.passivity{k});
    end
    
end
end


function s=checkBR(E,A,B,C,D,f) %check bounded realness
if length(f)==1
    omega=2*pi*f;
    H=C*((E*1i*omega-A)\B)+D;
    G=H*H';
    lambda=svd(G);
    if max(lambda)<=1
        s='passive';
    else
        s='non-passive';
    end
else
    if isinf(f(2))
        omega=2*pi*max(1.5*f(1),f(1)+50);
    else
        omega=pi*sum(f);
    end
    H=C*((E*1i*omega-A)\B)+D;
     G=H*H';
     lambda=svd(G);
    if  max(lambda)<=1
        s='passive';
    else
        s='non-passive';
    end
end

end