function [A,B,C,D,E,dx,x,u,y,Delay] = interpolateR(t,p,S,INFO)
% Fast interpolation of A,B,C,D,E,dx,x,u,y data over rectangular grid 
% (only for IMETHOD=linear/nearest and EMETHOD=linear/nearest/clip).

%   Copyright 2022 The MathWorks, Inc.

% NOTE: Avoid renaming or changing this function (breaks load/save).
gvecs = INFO.GridVectors;

% Construct x value
tidx = INFO.tidx;
p = p(:);
if isempty(tidx)
   x = p;
else
   x = [p(1:tidx-1); t; p(tidx:end)];
end
if numel(x)~=INFO.pdim
   error(message('Control:ltiobject:ssInterpolant16',INFO.pdim-numel(tidx)))
end

% Interpolate 
% Note: S is empty iff nothing varies
if ~isempty(S)
   EMETHOD = strrep(INFO.Extrapolation,'clip','boundary');
   C = matlab.internal.math.gridStructsInterp(...
      gvecs,{S},num2cell(x),INFO.Interpolation,EMETHOD);
   Si = C{1};
end

% Matrix data
% A
if INFO.A.nv==0
   A = INFO.A.Nominal;
elseif INFO.A.nf==0
   A = Si.A;
else
   A = INFO.A.Nominal;
   A(INFO.A.iv) = Si.A;
end
% B
if INFO.B.nv==0
   B = INFO.B.Nominal;
elseif INFO.B.nf==0
   B = Si.B;
else
   B = INFO.B.Nominal;
   B(INFO.B.iv) = Si.B;
end
% C
if INFO.C.nv==0
   C = INFO.C.Nominal;
elseif INFO.C.nf==0
   C = Si.C;
else
   C = INFO.C.Nominal;
   C(INFO.C.iv) = Si.C;
end
% D
if INFO.D.nv==0
   D = INFO.D.Nominal;
elseif INFO.D.nf==0
   D = Si.D;
else
   D = INFO.D.Nominal;
   D(INFO.D.iv) = Si.D;
end
% E
if INFO.E.nv==0
   E = INFO.E.Nominal;
elseif INFO.E.nf==0
   E = Si.E;
else
   E = INFO.E.Nominal;
   E(INFO.E.iv) = Si.E;
end

% Offset data
% dx
if INFO.dx.nv==0
   dx = INFO.dx.Nominal;
elseif INFO.dx.nf==0
   dx = Si.dx;
else
   dx = INFO.dx.Nominal;
   dx(INFO.dx.iv,:) = Si.dx;
end
% x
if INFO.x.nv==0
   x = INFO.x.Nominal;
elseif INFO.x.nf==0
   x = Si.x;
else
   x = INFO.x.Nominal;
   x(INFO.x.iv,:) = Si.x;
end
% u
if INFO.u.nv==0
   u = INFO.u.Nominal;
elseif INFO.u.nf==0
   u = Si.u;
else
   u = INFO.u.Nominal;
   u(INFO.u.iv,:) = Si.u;
end
% y
if INFO.y.nv==0
   y = INFO.y.Nominal;
elseif INFO.y.nf==0
   y = Si.y;
else
   y = INFO.y.Nominal;
   y(INFO.y.iv,:) = Si.y;
end

% Delays
if INFO.DELAY
   CT = INFO.CT;  % true->continuous time
   % Input
   if INFO.uDelay.nv==0
      uDelay = INFO.uDelay.Nominal;
   elseif INFO.uDelay.nf==0
      if CT
         uDelay = max(0,Si.uDelay);
      else
         uDelay = max(0,round(Si.uDelay));
      end
   else
      uDelay = INFO.uDelay.Nominal;
      if CT
         uDelay(INFO.uDelay.iv,:) = max(0,Si.uDelay);
      else
         uDelay(INFO.uDelay.iv,:) = max(0,round(Si.uDelay));
      end
   end
   % Output
   if INFO.yDelay.nv==0
      yDelay = INFO.yDelay.Nominal;
   elseif INFO.yDelay.nf==0
      if CT
         yDelay = max(0,Si.yDelay);
      else
         yDelay = max(0,round(Si.yDelay));
      end
   else
      yDelay = INFO.yDelay.Nominal;
      if CT
         yDelay(INFO.yDelay.iv,:) = max(0,Si.yDelay);
      else
         yDelay(INFO.yDelay.iv,:) = max(0,round(Si.yDelay));
      end
   end
   % Internal
   if INFO.xDelay.nv==0
      xDelay = INFO.xDelay.Nominal;
   elseif INFO.xDelay.nf==0
      if CT
         xDelay = max(0,Si.xDelay);
      else
         xDelay = max(0,round(Si.xDelay));
      end
   else
      xDelay = INFO.xDelay.Nominal;
      if CT
         xDelay(INFO.xDelay.iv,:) = max(0,Si.xDelay);
      else
         xDelay(INFO.xDelay.iv,:) = max(0,round(Si.xDelay));
      end
   end
   Delay = struct('Input',uDelay,'Output',yDelay,'Internal',xDelay);
else
   Delay = [];
end

