function [op,SINGULAR] = findop(D,opspec)
% Compute operating condition.

%   Copyright 2023 The MathWorks, Inc.
K = D.k;
nq = size(K,1);
[ny,nu] = iosize(D);
nfd = numel(D.Delay.Internal);
opspec = validate_q(opspec,nq,ny,nu,nfd);
M = D.m;
if isempty(M)
   M = eye(nq,'like',K);
end
C = D.c;
F = D.f;
G = D.g;
if D.Ts~=0
   % dq[k] = q[k+1]-q[k]  (for consistency with dx[k])
   % d2q[k] = dq[k+1]-dq[k] = q[k+2]+q[k]-2*q[k+1] (so that d2q[k]=0 <=> dq[k+1]=dq[k])
   % With these definitions, 
   %   M q[k+2] + C q[k+1] + K q[k] = M d2q[k] + (C+2*M) dq[k] + (M+C+K) q[k]
   K = K+C+M;
   C = C+2*M;
   F = F+G;
end
[op,SINGULAR] = ltipack.numerics.findopSO(M,C,K,D.b,F,G,D.d,...
   opspec.q,opspec.u,opspec.w,opspec.dq,opspec.d2q,opspec.y);
if D.Ts~=0
   % Return dq=q[k+1] and d2q=q[k+2] 
   op.d2q = op.d2q+2*op.dq+op.q;
   op.dq = op.dq+op.q;
end
