function varargout = hasDelayDynamics(D,PZFLAG)
% Checks if dynamics (poles, zeros, or both) depend on internal delays.
%
%   [DP,DY] = hasDelayDynamics(D,'pole') checks if the poles depend on 
%   nonzero internal delays (DP=true). When DP=false, the model equations  
%   are of the form
%      dx = dx0 + A (x-x0) + sum_j Bj (u(t-tj)-u0)
%       y = y0 + C (x-x0) + sum_j Cj (x(t-tj)-x0) + Dj (u(t-tj)-u0)
%   and can be discretized for simulation. DY=false further indicates that
%   there are no state delays in the output equation, i.e.
%       y = y0 + C (x-x0) + sum_j Dj (u(t-tj)-u0)
%   For static models, DP=DY=false when the output equation is of the form
%       y = y0 + sum_j Dj (u(t-tj)-u0) .
%   
%   DZ = hasDelayDynamics(D,'zero') checks if the transmission zeros depend 
%   on nonzero internal delays (DZ=true). When DZ=false, the internal delays 
%   are equivalent to input and output delays.
%
%   DPZ = hasDelayDynamics(D) returns true when either the poles or zeros
%   depend on nonzero internal delays.

%   Copyright 1986-2023 The MathWorks, Inc.
if nargin<2
   varargout{1} = localDP(D) || localDZ(D);
elseif strcmp(PZFLAG,'pole')
   [varargout{1:nargout}] = localDP(D);
else
   varargout{1} = localDZ(D);
end

%------------------------------------
function [DP,DY] = localDP(D)
% Checks if poles depend on nonzero internal delays.
% Note: Offsets can be ignored since the structure of (D22,C2,B2) is what
% makes the finite expansion of w(t)-w0 possible. 
xDelay = D.Delay.Internal;
i2 = find(xDelay>0);
if isempty(i2)
   DP = false;  DY = false;
else
   [ny,nu] = iosize(D);
   d22 = D.d(ny+i2,nu+i2);
   if ltipack.util.isNilpotent(d22~=0)
      % Independence requires D22 nilpotent
      nx = size(D.a,1);
      i1 = find(xDelay==0);
      b2 = [D.b(:,nu+i2);D.d(ny+i1,nu+i2)];
      c2 = [D.c(ny+i2,:),D.d(ny+i2,nu+i1)];
      DP = nx>0 && ~isempty(smreal(d22,c2,b2,[]));
      if nargout>1
         d12 = D.d([1:ny,ny+i1'],nu+i2);
         DY = nx>0 && ~isempty(smreal(d22,c2,d12,[]));
      end
   else
      DP = true;  DY = true;
   end
end

%------------------------------------
function DZ = localDZ(D)
% Checks if zeros depend on internal delays
% REVISIT: Only sufficient condition for DPZ=0, misses 1/(1+exp(-tau*s))
nfd = numel(D.Delay.Internal);
if nfd>0
   iod = getIODelay(D);
   delta = diff(diff(iod,1,1),1,2);
   DZ = ~allfinite(iod) || norm(delta,'fro')>1e-6*norm(iod,'fro');
else
   DZ = false;
end
