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

%	 Author(s): A. Grace, P. Gahinet
%   Copyright 1986-2023 The MathWorks, Inc.
Ts0 = D1.Ts;
if Ts0==0
    % Static gain with Ts=0
    D2 = D1;  D2.Ts = Ts; return
end
CONSISTENT = (nargin>2 && strcmp(options.Consistency,'on'));

switch options.Method(1)
   case 't'
      % Tustin resampling
      w = options.PrewarpFrequency;
      if w==0
         r = Ts/Ts0;
      else
         r = tan(w*Ts/2)/tan(w*Ts0/2);
      end
      
      % Compute delays at new sampling time and round to
      % nearest sampling period
      Delay = D1.Delay;
      Delay.Input = Delay.Input * Ts0;
      Delay.Output = Delay.Output * Ts0;
      Delay.Internal = Delay.Internal * Ts0;
      [Delay,fid,fod,ffd] = ltipack.ssdata.discretizeDelay(Delay,Ts);
      
      % Round delays to nearest multiple of Ts
      if (any(fid) || any(fod) || any(ffd)) && ~CONSISTENT
         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
      tau = (r-1)/2;
      a1 = D1.a;  b1 = D1.b;  e1 = D1.e;  dx01 = D1.dx0;  x0 = D1.x0;
      nx = size(a1,1);
      nuw = size(D1.d,2);
      if isempty(e1)
         % Explicit SS
         a2 = ((1+r)*a1+(1-r)*eye(nx))/2;
         e2 = ((1+r)*eye(nx)+(1-r)*a1)/2;
         % Prevent scaling-induced "near singularity" (g330910)
         dx01 = dx01-x0;
         [aux,c2] = ltipack.util.safeMdivide(e2,[a2 r*b1 dx01],D1.c);
         if ~(allfinite(aux) && allfinite(c2))
            error(message('Control:transformation:d2d08'))
         end
         a2 = aux(:,1:nx);      % e2\a2=a2/e2
         b2 = aux(:,nx+1:nx+nuw);  % r*e2\b
         d2 = D1.d + tau * c2 * b1;  % d + (r-1)/2 * c/e2*b
         e2 = [];
         dx02 = r*aux(:,nx+nuw+1)+x0;
         y02 = D1.y0 + tau * (c2 * dx01);
      else
         % Descriptor SS
         a2 = ((1+r)*a1+(1-r)*e1)/2;
         e2 = ((1+r)*e1+(1-r)*a1)/2;  % e2
         [~,aux] = ltipack.util.safeMdivide(e2,[],D1.c); % c1/e2
         if ~allfinite(aux)
            error(message('Control:transformation:d2d08'))
         end
         b2 = r*b1;
         c2 = aux*e1;
         d2 = D1.d + tau * (aux*b1);
         dx01 = dx01-e1*x0;
         dx02 = r*dx01+e2*x0;
         y02 = D1.y0 + tau * (aux * dx01);
      end
      D2 = ltipack.ssdata(a2,b2,c2,d2,e2,Ts);
      D2.dx0 = dx02;
      D2.x0 = x0;
      D2.u0 = D1.u0;
      D2.y0 = y02;     
      D2.Delay = Delay;
      
      % Eliminate zero internal delays (after rounding)
      % Note: May result in improper model
      D2 = elimZeroDelay(D2);
      
   case 'z'
      % ZOH resampling
      r = Ts/Ts0; % sample time ratio
      FractionalResampling = (abs(round(r)-r)>sqrt(eps)*r);
      if ~FractionalResampling
         r = round(r);
      end
      
      % Check causality
      [isProper,D1] = isproper(D1,true);
      if ~isProper
         error(message('Control:transformation:NotSupportedImproperIZF'))
      end
      
      % Delay handling
      id = D1.Delay.Input/r;    rid = round(id);
      od = D1.Delay.Output/r;   rod = round(od);
      fd = D1.Delay.Internal/r; rfd = round(fd);
      tolint = 1e4*eps;
      FractionalDelay = (any(abs(id-rid)>tolint) || ...
         any(abs(od-rod)>tolint) || any(abs(fd-rfd)>tolint));
      
      % Extract data (note: E=[] from ISPROPER call)
      a1 = D1.a;  b1 = D1.b;  c1 = D1.c;  d1 = D1.d;  dx01 = D1.dx0-D1.x0;
      
      % Look for real negative poles
      if FractionalResampling
         p = eig(a1);
         RealNegPole = any(imag(p)==0 & real(p)<=0);
      end
      
      % Resample
      if FractionalDelay || (FractionalResampling && RealNegPole)
         % Fractional delays or negative real poles with fractional resampling:
         % let D2C/C2D handle it
         try
            pC2D = c2dOptions(Method='zoh',Consistency=options.Consistency);
            pD2C = d2cOptions(Method='zoh',Consistency=options.Consistency);
            D2 = c2d(d2c(D1,pD2C),Ts,pC2D);
         catch E
            error(E.identifier,strrep(E.message,'d2c','d2d'))
         end
      else
         % Proceed directly
         nx = size(a1,1);
         nuw = size(d1,2);
         % Resample
         M = [a1 b1 dx01;zeros(nuw+1,nx) eye(nuw+1)];
         if FractionalResampling
            [s,~,M] = mscale(M,'noperm','safebal');
            M = expm(r*logm(M));
            M = lrscale(M,s,1./s);
         else
            M = M^r;
         end
         if ~allfinite(M)
            % Check for overflow
            error(message('Control:transformation:d2d11'))
         elseif isreal(a1) && isreal(b1)
            M = real(M);
         end
         a2 = M(1:nx,1:nx);  
         b2 = M(1:nx,nx+1:nx+nuw);
         dx02 = M(1:nx,nx+nuw+1);

         % Create result
         D2 = ltipack.ssdata(a2,b2,c1,d1,[],Ts);
         D2.dx0 = dx02+D1.x0;
         D2.x0 = D1.x0;
         D2.u0 = D1.u0;
         D2.y0 = D1.y0;
         D2.Delay.Input = rid;
         D2.Delay.Output = rod;
         D2.Delay.Internal = rfd;
         D2.StateName = D1.StateName;
         D2.StatePath = D1.StatePath;
         D2.StateUnit = D1.StateUnit;
      end
end
