function [x0,w0,y0,ypast] = resolveIC(D,x0,u0)
% Resolve initial condition for STEP or IMPULSE.
%
%   u(t) is assumed constant and equal to u0 (input offset) prior to start.
%
%   x0=[] is resolved to a steady initial condition (x0,u0,w0) satisfying
%        x0 = [q0 ; 0]
%      K q0 = B1 u0 + B2 w0    or (M+C+K) q0 = B1 u0 + B2 w0
%        w0 = C2 x0 + D21 u0 + D22 w0 .
%   The second equation expresses q'=q''=0 or q[k]=q[k+1]=q[k+2]. The 
%   The third equation says that the internal signals z(t) are steady with
%   value w0, which is required to condense the z(t<t0) history into a
%   single value w0. The value y0=[F1 G1]*x0+D11*u0+D12*w0 is the output 
%   at time t0, and this value is steady for t<t0 so ypast=y0.
%
%   When x0=x(t0) is specified, the model is required to have no state delays 
%   and be of the form
%      M q'' + C q' + K q =  sum Bj u(t-tj)
%      y = F q + G q' + sum Dj u(t-tj)
%   (see checkCompatibility). It is then enough to consider the dependence
%   of z(t) on u(t), which can be initialized to the steady value w0
%   satisfying
%      w0 = D21 u0 + D22 w0 .
%   The output at time t0 is y0=[F1 G1]*x0+D11*u0+D12*w0. This value is not
%   steady for t<t0 except for entries y0(i) that do not depend on x due to
%   F1(i,:)=G1(i,:)=0. ypast coincides with y0 for such entries and is set
%   to zero elsewhere (unknown or unsteady past value).

%   Copyright 1986-2022 The MathWorks, Inc.
if isempty(x0)
   % Compute steady-state IC for input offset u0
   opspec = ltioptions.findop;
   opspec.u = u0;
   [op,SINGULAR] = findop(D,opspec);
   if SINGULAR
      error(message('Control:analysis:step2'))
   end
   x0 = [op.q;op.dq];
   w0 = op.w;
   y0 = op.y;
   ypast = y0;
else
   % Specified x0, NO state delays
   [~,~,~,~,~,f1,~,g1,~,d11,d12,d21,d22] = getBlockData(D);
   fg1 = [f1 g1];
   nfd = size(d22,1);
   rhs = full(d21*u0);
   if any(rhs)
      w0 = matlab.internal.math.nowarn.mldivide(eye(nfd,'like',d22)-d22,rhs);
   else
      w0 = zeros(nfd,1);
   end
   y0 = fg1*x0+d11*u0+d12*w0;
   ypast = y0;
   ypast(any(fg1,2),:) = 0;
end
