function fr = evalfr(D,s)
%EVALFR  Evaluates frequency response at a single (complex) frequency.

%   Author(s):  P. Gahinet  5-13-96
%   Copyright 1986-2007 The MathWorks, Inc.
[ny,nu] = size(D.num);
fr = zeros(ny,nu);
Td = getIODelay(D,'total');  % total I/O delays
Ts = D.Ts;
if isnan(s)
   fr(:) = NaN;
elseif isinf(s)
   % Compute feedthrough
   if Ts==0
      % limit for s->+Inf (|s|->Inf with angle(s)->0)
      for ct=1:ny*nu
         if Td(ct)~=0  
            % Delay -> no feedthrough
            fr(ct) = 0;
         elseif D.den{ct}(1)==0
            % Improper
            fr(ct) = Inf;
         else
            fr(ct) = D.num{ct}(1)/D.den{ct}(1);
         end
      end
   else
      % limit for |z|->Inf
      for ct=1:ny*nu
         iLD = find(D.den{ct}~=0,1)-1;
         if iLD>Td(ct)
            % Improper
            fr(ct) = Inf;
         elseif iLD==Td(ct)
            fr(ct) = D.num{ct}(1)/D.den{ct}(1+iLD);
         else
            fr(ct) = 0;
         end
      end
   end
else
   % Response at finite point
   for ct=1:ny*nu
      Denom = polyval(D.den{ct},s);
      if Denom==0
         fr(ct) = Inf;
      else
         fr(ct) = polyval(D.num{ct},s) / Denom;
      end
   end
   if any(Td(:))
      if Ts
         fr = s.^(-Td) .* fr;
      else
         fr = exp(-s*Td) .* fr;
      end
   end
end
