function [Dc,gic] = d2c(Dd,options)
%D2C  Conversion of discrete state-space models to continuous time.

%   Author(s): P. Gahinet
%   Copyright 1986-2012 The MathWorks, Inc.
method = options.Method(1);

switch method
   case 'm'
      % Use zpkdata algorithm (Note: zpk conversion never fails for
      % discrete models)
      Dc = ss(d2c(zpk(Dd),options));
      gic = [];  % not supported
      
   case {'z','f'}
      % Check causality and get rid of E matrix
      [isProper,Dd] = isproper(elimZeroDelay(Dd),true);
      if ~isProper
         error(message('Control:transformation:NotSupportedImproperIZF'))
      end
      % Zero- or first-order hold approximation
      [Dc,gic] = utInvDiscretize(Dd,options);
      
   case 't'
      % Tustin approximations
      % Handle prewarp      
      w = options.PrewarpFrequency;
      if w == 0
          r = 1/Dd.Ts;
      else
          % Handle prewarping
          r = w/2/tan(w*Dd.Ts/2);
      end

      dx0d = Dd.dx0;  x0d = Dd.x0;
      nx = size(Dd.a,1);
      nuw = size(Dd.d,2);
      if isempty(Dd.e)
         % Explicit SS
         % Prevent scaling-induced "near singularity" (g330910)
         m = eye(nx)+Dd.a;
         [X,Y] = ltipack.util.safeMdivide(m, (2*r)*[Dd.a-eye(nx),Dd.b,dx0d-x0d], Dd.c);
         if ~(allfinite(X) && allfinite(Y))
            error(message('Control:transformation:d2c06'))
         end
         ac = X(:,1:nx);         % 2/T*(Ad+I)\(Ad-I)
         bc = X(:,nx+1:nx+nuw);  % 2/T*(Ad+I)\B
         cc = 2*Y;            % 2*C/(Ad+I)
         dc = Dd.d - Y*Dd.b;  % D-C/(Ad+I)*B
         ec = [];
         dx0 = X(:,nx+nuw+1); % 2/T*(Ad+I)\(dx0d-x0d)
         y0 = Dd.y0-(Dd.c*dx0)/(2*r);
         gic = [m , Dd.b]/2;
      else
         % Descriptor SS
         ad = Dd.a;  ed = Dd.e;  bd = Dd.b;
         ec = (ad+ed)/2;
         ac = r*(ad-ed);
         bc = r*bd;
         [~,Y] = ltipack.util.safeMdivide(ec,[],Dd.c/2); % Cd/(Ad+Ed)
         if ~allfinite(Y)
            error(message('Control:transformation:d2c06'))
         end
         cc = 2*(Y*ed);
         dc = Dd.d - Y*bd;
         dx0d = dx0d-ed*x0d;
         dx0 = r*dx0d;
         y0 = Dd.y0-Y*dx0d;
         gic = matlab.internal.math.nowarn.mldivide(ed,[ec , bd/2]);
      end
      
      % Zero out columns of GIC where there are input delays
      iud = find(Dd.Delay.Input);
      gic(:,nx+iud) = 0;

      % Update delays and sample time
      % Discard state names because of bilinear transformation
      Dc = ltipack.ssdata(ac,bc,cc,dc,ec,0);
      Dc.dx0 = dx0;
      Dc.x0 = x0d;
      Dc.u0 = Dd.u0;
      Dc.y0 = y0;
      Ts = Dd.Ts;
      Dc.Delay.Input = Ts * Dd.Delay.Input;
      Dc.Delay.Output = Ts * Dd.Delay.Output;
      Dc.Delay.Internal = Ts * Dd.Delay.Internal;
end


