function fresp = evalfr(D,s)
%EVALFR  Evaluates frequency response at a single (complex) frequency.
%                                                  -1
%       FRESP =  D + (F + s G) * (M s^2 + C s + K )  * B   .

%   Copyright 2020 The MathWorks, Inc.
Ts = D.Ts;
% Compute H(s)
if isfinite(s)
   % Jointly close the integrator and delay loops (the A-sI or D22-exp(-tau*s)
   % blocks could be singular, see g337144)
   d = D.d;
   k = D.k;
   m = D.m;
   f = D.f;
   if isempty(m)
      m = speye(size(k));
   end
   if nnz(k)>0 || nnz(f)>0
      % Solve for [q;w]
      MCK = m*s^2+D.c*s+k;
      FG = f + s*D.g;
   else      
      % Solve for [s*q;w] to avoid 0/0 cancellation when s=0
      MCK = m*s+D.c;
      FG = D.g;
   end
   
   % Sizes
   [rs,cs] = size(d);
   tau = D.Delay.Internal;
   nfd = length(tau);
   ny = rs-nfd;
   nu = cs-nfd;
   
   % Form linear system
   if nfd>0
      % Note: K may be invertible when its (1,1) or (2,2) blocks are not
      if Ts==0
         etau = exp(tau*s);
      else
         etau = s.^tau;
      end
      MCK = [MCK , D.b(:,nu+1:nu+nfd) ; ...
         FG(ny+1:ny+nfd,:) , diag(etau)-d(ny+1:ny+nfd,nu+1:nu+nfd)];
      b = [D.b(:,1:nu) ; -d(ny+1:ny+nfd,1:nu)];
      FG = [FG(1:ny,:) , -d(1:ny,nu+1:nu+nfd)];
      d = d(1:ny,1:nu);
   else
      b = D.b;
   end
   
   % Evaluate response by solving for [q;w] or [s*q;w]
   if rs>=cs
      aux = ltipack.util.safeMdivide(MCK,b,[]);
      FINITE = allfinite(aux);
      if FINITE
         fresp = full(d + FG * aux);
      else
         fresp = inf(ny,nu);
      end
   else
      [~,aux] = ltipack.util.safeMdivide(MCK,[],FG);
      FINITE = allfinite(aux);
      if FINITE
         fresp = full(d + aux * b);
      else
         fresp = inf(ny,nu);
      end
   end
   
   if FINITE
      % Add delay contribution
      id = D.Delay.Input;
      od = D.Delay.Output;
      if any(id) || any(od)
         if Ts==0
            fresp = fresp .* exp(-s*(id'+od));
         else
            fresp = fresp .* s.^(-(id'+od));
         end
      end
   end
else
   fresp = nan(iosize(D));
end
