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

 [n, r, m] = abcddim (Ac, Bc, Cc, Dc);

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

 % continuous-time to discrte-time conversion ...
 if strcmp(lower(ntrp),'zoh')      % zero-order hold on inputs
   M    = [ Ac Bc ; zeros(r,n+r) ];
 else                              % first-order hold on inputs
   M    = [ Ac Bc zeros(n,r) ; zeros(r,n+r) eye(r) ; zeros(r,n+2*r) ];
 end
 eMdt = expm(M*dt);                % matrix exponential
 Ad   = eMdt(1:n,1:n);             % discrete-time dynamics matrix
 Bd   = eMdt(1:n,n+1:n+r);         % discrete-time input matrix
 if strcmp(lower(ntrp),'zoh')
   Cd   = Cc;
   Dd   = Dc;
 else
   Bdp  = eMdt(1:n,n+r+1:n+2*r)    % discrete-time input matrix
   Bd = Bd  + (Ad - eye(n)) * Bdp / dt; 
   Cd = Cc;
   Dd = Dc  + Cc*Bdp/dt;
 end

% ======================================================================= con2dis
% HP Gavin, 2020-11-01, 2021-07-19, 2024-01-23
