function [Ac,Bc,Cc,Dc] = dis2con(Ad,Bd,Cd,Dd,dt,ntrp);
% [Ac,Bc,Cc,Dc] = dis2con(Ad,Bd,Cd,Dd,dt,ntrp);
% discrete time to continuos time transformation 
% Ac must be invertible ... 
% ntrp =  'foh' ... first order hold, default
% ntrp =  'zoh' ... zero  order hold, optional

 [n, r, m] = abcddim (Ad, Bd, Cd, Dd);

 if nargin < 6, ntrp = 'foh'; end

% first solve for zero-order hold
 NumSqrt = 1e2;    % number of squqre roots in the logm calculation
 M = logm( [ Ad Bd ; zeros(r,n) eye(r) ] ) / dt; 
 Ac = M(1:n,1:n);
 Bc = M(1:n,n+1:n+r);
 Cc = Cd; 
 Dc = Dd;

if strcmp(ntrp,'foh')

  % todo: the following three lines should be computed using expm(?) ...
  invAc = inv(Ac);
  A1  = invAc * ( Ad - eye(n) );
  A2  = invAc*invAc * ( Ad - eye(n) - Ac*dt ); 
  Bc  = ( A1 + (Ad-eye(n)) * A2 / dt ) \ Bd;
  Bdp = A2*Bc;
  Dc  = Dd - Cc*Bdp/dt;

end

% ================================================================= dis2con
% HP Gavin,  2023-03-06, 2024-01-23, 2024-03-25
