function D2 = d2d(D1,Ts2,options)
% Resample discrete state-space model to Ts2.

%   Copyright 2020 The MathWorks, Inc.
Ts1 = D1.Ts;
if Ts1==0
   % Static gain with Ts=0
   D2 = D1;  D2.Ts = Ts2; return
end

if strcmp(options.Method,'tustin')
   % Tustin resampling
   w = options.PrewarpFrequency;
   if w==0
      r = Ts2/Ts1;
   else
      r = tan(w*Ts2/2)/tan(w*Ts1/2);
   end
   m1 = D1.m;  c1 = D1.c;  k1 = D1.k;  b1 = D1.b;  f1 = D1.f;  g1 = D1.g;
   if isempty(m1)
      m1 = speye(size(k1));
   elseif nnz(m1)==0
      % Use first-order Tustin when M=0
      D2 = mechss(d2d(sparss(D1),Ts2,options));
      return
   end

   % Compute delays at new sampling time and round to
   % nearest sampling period
   Delay = D1.Delay;
   Delay.Input = Delay.Input * Ts1;
   Delay.Output = Delay.Output * Ts1;
   Delay.Internal = Delay.Internal * Ts1;
   [Delay,fid,fod,ffd] = ltipack.ssdata.discretizeDelay(Delay,Ts2);
   
   % Round delays to nearest multiple of Ts2
   if any(fid) || any(fod) || any(ffd)
      warning(message('Control:transformation:RoundedDelay'))
   end
   Delay.Input = Delay.Input + round(fid);
   Delay.Output = Delay.Output + round(fod);
   Delay.Internal = Delay.Internal + round(ffd);
   
   % Compute matrices for resampled model
   rp = (r+1)/2;   rm = (r-1)/2;
   m2 = rp^2*m1+rp*rm*c1+rm^2*k1;
   c2 = ((r^2-1)/2)*(m1+k1)+((r^2+1)/2)*c1;
   k2 = rm^2*m1+rp*rm*c1+rp^2*k1;
   b2 = r^2*b1;
   [~,aux] = ltipack.util.safeMdivide(m2,[],(r-1)*f1+(r+1)*g1);
   if ~allfinite(aux)
      error(message('Control:transformation:d2d08'))
   end
   alpha = (1+r^(-2))/2;  beta = (1-r^(-2))/2;  gamma = rm/2/r^2;
   f2 = alpha*f1+beta*g1-gamma*aux*(m2+k2);
   g2 = beta*f1+alpha*g1-gamma*aux*c2;
   d2 = D1.d+(rm/2)*aux*b1;
   D2 = ltipack.mechdata(m2,c2,k2,b2,f2,g2,d2,Ts2);
   D2.Delay = Delay;
   D2.StateInfo = D1.StateInfo;

   % Eliminate zero internal delays (after rounding)
   % Note: May result in improper model
   D2 = elimZeroDelay(D2);
else
   error(message('Control:transformation:d2d04'))   
end
