function thd=thc2thd(thc,T)
%THC2THD  Converts a continuous time model to discrete time	
%
%	THD = thc2thd(THC,T)
%	
%	THC: The continuous time model, specified in THETA format 
%		(See help theta)
%	T: The sampling interval
%	THD: The discrete time model, in THETA format
%
%	Note that the covariance matrix is not translated to discrete time
%	for input-output type models
%
%	See also THD2THC.

%	L. Ljung 10-2-90
%	Copyright (c) 1990 by the MathWorks, Inc.
%	All Rights Reserved.

Told=gett(thc);
if nargin<2,error('The sampling time T must be specified!'),end
if Told>0,error('This model is already denoted as discrete-time!'),end
if T<0,error('The sampling interval must be positive!'),end
lamscale=abs(Told)/T;
[nrth,ncth]=size(thc);[p,covp,lam]=th2par(thc);nnd=length(p);nu=thc(1,3);
if isthss(thc),rarg=thc(1,6);ny=thc(1,4);
thd=thc; thd(1,2)=T;
thd(4+nnd+rarg:nrth,1:ny)=lamscale*lam;
thd(1:2,1)=lamscale^ny*thc(1:2,1);
return,end 
if nrth>3+nnd,delays=thc(nrth,1:nu);else delays=zeros(1,nu);end

nnk=[round(delays/T) 0];
[a,b,c,d,k]=th2ss(thc);nu=thc(1,3);
[m,n] = size(a);
[m,nb] = size(b);
[m,nk] = size(k);

s = expm([[a b k]*T; zeros(nb+nk,n+nb+nk)]);
a = s(1:n,1:n);
bk = s(1:n,n+1:n+nb+nk);
dk = [d 1];

den = poly(a);

[mc,nc] = size(c);
for iu=1:nu+1
b1 = bk(:,iu);
d1 = dk(:,iu);

	num1 = poly(a-b1*c) + (d1 - 1) * den;
	num(iu,1:length(num1)+nnk(iu))=[zeros(1,nnk(iu)) num1];
end

thd=poly2th(den,num(1:nu,:),num(nu+1,:),1,ones(nu,1),lamscale*lam,T);
thd(2,1)=thc(2,1)*lamscale;
end