function ID_err = SysID_Err(Ac,Bc,C,D , A_id,B_id,C_id,D_id, dt, FigNo );
% ID_err = SysID_Err(Ac,Bc,C,D , Aid,Bid,Cid,Did, dt, FigNo )
% Compute the ratio of H2 norms of the difference of two LTI systems:
%  (Ac,Bc,C,D)       - a  "reference"  model system in continuous-time
%  (Aid,Bid,Cid,Did) - an "identified" model system in discrete-time
%  dt is the sample time for the identified discrete-time system
% if FigNo > 0, a Bode plot of the difference of "ref" - "id" is plotted

 if nargin < 10, FigNo = 0; end

 [n,r,m] = abcddim ( Ac  ,Bc   ,C   ,D    );
 n_id    = abcddim ( A_id,B_id,C_id,D_id );

% continuous-time to discrete-time transformation of the reference system
 AB = expm([ Ac Bc ; zeros(r,n+r) ]*dt);  A = AB(1:n,1:n);  B = AB(1:n,n+1:n+r);

% put the "id" system in parallel with the "reference" system
 Ae_id = [ A_id          , zeros(n_id,n) ;
           zeros(n,n_id) , A              ];
 Be_id = [ B_id ;
           B    ];
 Ce_id = [ C_id , -C ];

 De_id = zeros(m,r);

% H2 norm of the difference of outputs "y_ref" - "y_id"
 model_err  = sqrt(real(trace( Ce_id*dliap(Ae_id,Be_id*Be_id')*Ce_id' )/dt));
% Hw norm of the "ref" system
 model_norm = sqrt(real(trace( C*dliap(A,B*B')*C' )/dt));

 ID_err = model_err / model_norm;

 if FigNo
   w = 2*pi*logspace(-1,log10(1/(5*dt)),200);
   mag_err = mimoBode(Ae_id,Be_id,Ce_id,De_id,w,dt,0);
   figure(FigNo)
     clf
     hold on
     for jj=1:r
       loglog( w/2/pi , mag_err(:,:,jj) );
     end
     title('model output error (id - reference)')
     xlabel('frequency, Hz')
     axis('tight')
     hold off
 end

% ----------------------------------------------------------------------------
% H.P. Gavin, 2017-11-07
