function fresp = evalfr(D,s)
%EVALFR  Evaluates frequency response at a single (complex) frequency.
%                                   -1
%       FRESP =  D + C * (X * E - A)  * B   .

%   Copyright 1986-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;
   e = D.e;
   if isempty(e)
      e = eye(size(D.a),'like',d);
   end
   
   % Sizes
   [rs,cs] = size(d);
   tau = D.Delay.Internal;
   nfd = length(tau);
   ny = rs-nfd;
   nu = cs-nfd;
   
   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
      K = [D.a-s*e , D.b(:,nu+1:nu+nfd) ; ...
         D.c(ny+1:ny+nfd,:) , d(ny+1:ny+nfd,nu+1:nu+nfd)-diag(etau)];
      b = [D.b(:,1:nu) ; d(ny+1:ny+nfd,1:nu)];
      c = [D.c(1:ny,:) , d(1:ny,nu+1:nu+nfd)];
      d = d(1:ny,1:nu);
   else
      K = D.a-s*e;
      b = D.b;
      c = D.c;
   end
   
   % Evaluate response
   if rs>=cs
      aux = ltipack.util.safeMdivide(K,b,[]);
      FINITE = allfinite(aux);
      if FINITE
         fresp = full(d - c * aux);
      else
         fresp = inf(ny,nu);
      end
   else
      [~,aux] = ltipack.util.safeMdivide(K,[],c);
      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
   % |s|=NaN or Inf
   if issparse(D) || isnan(s)
      fresp = nan(iosize(D));
   else
      % Compute feedthrough term
      if Ts==0 || isempty(D.e)
         % Continuous or explicit discrete: Get rid of delays using exp(-s*tau)->0
         % and 1/z^n->0
         D = localRemoveDelays(D);
      else
         % Discrete descriptor. Map all delays to poles at z=0 to assess properness
         D = elimDelay(D);
      end
      % Compute feedthrough of delay-free model
      [isProper,D] = isproper(D);
      if isProper
         fresp = D.d;
      else
         fresp = inf(size(D.d));
      end
   end
end


function D = localRemoveDelays(D)
% Replaces nonzero delays by zero
Delay = D.Delay;
% Replace nonzero input and output delays by zero
iu = find(Delay.Input>0);
iy = find(Delay.Output>0);
D.b(:,iu) = 0;  D.d(:,iu) = 0;
D.c(iy,:) = 0;  D.d(iy,:) = 0;
% Replace nonzero internal delays by zero and close LFT
nfd = numel(Delay.Internal);
if nfd>0
   [nr,nc] = size(D.d);
   D.Delay.Internal = zeros(0,1);
   D.Delay.Input = zeros(nc,1);
   D.Delay.Output = zeros(nr,1);
   D = lft(D,createGain(D,diag(Delay.Internal==0)),...
      nc-nfd+1:nc,nr-nfd+1:nr,1:nfd,1:nfd);
end
