function [fc,lambda,L]=ghm(E,A,B,C,D,tol,opt)
% use generalized Hamiltonian method (GHM) and its symmetric variant to check
% the passivity of an immittance 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 GHM 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 test of immittance descriptor systems 
% based on generalized Hamiltonian methods",  IEEE Trans. Circuits and Systems II: 
% Express Briefs, vol. 57, no. 1, pp. 61-65, Jan 2010 
% [2] Z. Zhang, C.-U. Lei and N. Wong, "GHM: a generalized Hamiltonian method
% for passivity test of impedance/admittance descriptor systems”, IEEE/ACM Int. 
% Conf. on Computer-Aided Design  (ICCAD),  pp. 767-773, San Jose, CA, Nov. 2009


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
rd=rank(D+D.');
if rd<m % perform equivalent model conversion
    alpha=0.5+0.5*rand(1,1); % generate a random number in [0.5,1] for model conversion
    Im=eye(m,m);
    Da=alpha*Im-D;
    E=blkdiag(E,zeros(m,m));
    A=blkdiag(A,Da\Im);
    B=[B;Im]; C=[C, Im];
    D=alpha*Im;
end

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


if ~opt % use full size GHM
    Ahat=A-B*((D+D.')\C);     Rhat=B*((D+D.')\B.');
    Phat=C.'*((D+D.')\C);     J=[Ahat,-Rhat; Phat, -Ahat.'];
    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 GHM
    Jh=A-B*(D\C);  Kh=E*(A\E);
    lambda=eig(Jh,Kh);
    id=find(isinf(lambda));
    lambda(id)=[];
    
    N=length(lambda);
    imagLambda=imag(lambda); 
    for k=1:N %find the 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}=checkPR(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}=checkPR(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=checkPR(E,A,B,C,D,f) % check positive realness
if length(f)==1
    omega=2*pi*f;
    H=C*((E*1i*omega-A)\B)+D;
    G=H+H';
    lambda=real(eig(G));
    if min(lambda)>=0
        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=real(eig(G));
    if min(lambda)>=0
        s='passive';
    else
        s='non-passive';
    end
end

end