function thc=thd2thc(th,delays,nop)
%THD2THC  converts a model to continuous-time form.
%
%	THC = thd2thc(TH)
%
%	TH: A discrete time model defined in the format described by HELP THETA
%	THC: The continuous-time counterpart, stored in THETA format
%
%	If the model TH has extra delays from some inputs (nk>1), these are
%	first removed. They should be appended to the continuous-time model
%	as a pure time-delay (deadtime) (this is done by TH2FF).
%	To have the delays approximated by THC, enter THC = thd2thc(TH,'del');
%
%	The covariance matrix P of TH is translated by the use of numerical
%	derivatives. The step sizes used for the differences are given by the
% 	m-file nuderst. To inhibit the translation of P, which saves time,
%	enter THC = thd2thc(TH,c,'nop'), where c is 'del' or 'nodel'. 
%
%	CAUTION: The transformation could be ill-conditioned. Negative
%	real poles, poles in the origin and in 1 may cause problems. See the
%	manual.
%	See also THC2THD

%	L. Ljung 10-1-86,4-20-87,1-13-90
%	Copyright (c) 1986-90 by the MathWorks, Inc.
%	All Rights Reserved.


T=th(1,2);
if T<=0, error('This model is already denoted as continuous time!'),end
if isthss(th), thc=th;thc(1,2)=-T;
if thc(2,8)==3, if ~exist('minreal')
disp('WARNING: Transformation of multioutput arx-models will be supported only if you have the CONTROL SYSTEMS TOOLBOX'),
end,end
return,end
nu=th(1,3);nk=th(1,7+2*nu:6+3*nu);
nn=th(1,3:6+2*nu);na=nn(2);nb=nn(3:2+nu);nc=nn(3+nu);nddd=nn(4+nu);nf=nn(5+nu:4+2*nu);
% *** Set up default values ***
if nargin<3,nop=[];end
if nargin<2,delays=[];end

ndd=na+sum(nb)+nc+nddd+sum(nf);

if isempty(delays),delays='nodel';end

if delays(1)=='d',delays=0;else delays=1;end
[par,P]=th2par(th);
if norm(P)==0 | isempty(P) | ~isempty(nop), nkder=0;else nkder=[0:ndd];end
dt=nuderst(par);
for kder=nkder
th1=th;if kder>0,th1(3,kder)=th(3,kder)+dt(kder);end	

[a,b,c,d,f]=th2poly(th1);
b(nu+1,1:length(c))=c;
f(nu+1,1:length(d))=d;
nk(nu+1)=0;nb(nu+1)=nc+1;nf(nu+1)=nddd;
ss=1;
thb=[];thcc=[];thd=[];thf=[];ncn=0;ndn=0;
for k=[nu+1,1:nu]
den=conv(a,f(k,1:nf(k)+1)); 
if delays~=0,nkmod=max(nk(k),1);else nkmod=1;end
num=b(k,nkmod:nb(k)+nk(k));
nd=length(den);, nn=length(num);
n=max(nd,nn);
de=zeros(1,n);, de(1:nd)=den;
if kder==0
if abs(de(n))<eps, disp('WARNING: Pole in the origin. Result may be unreliable!'),end
rp=roots(de);
if any(abs(rp-1))<eps disp('WARNING: Pure integration. Result may be unreliable!'),end
if any((rp==real(rp)).*(real(rp)<0)),disp('WARNING: Pole on the negative real axis. Result may be unreliable!'),end
end
nume=zeros(1,n);, nume(1:nn)=num;
A=[-de(2:n);eye(n-2,n-1)]; % Transform to state-space
B=eye(n-1,1);
C=nume(2:n)-nume(1)*de(2:n);
D=nume(1); 
s=real(logm([[A B];zeros(1,n-1) eye(1)]))/T;  % Sample
Ac=s(1:n-1,1:n-1);
Bcc=s(1:n-1,n);
ff=poly(Ac);                % Transform to i/o form
bb=poly(Ac-Bcc*C)+(D-1)*ff;
if nk(k)>0,bb=bb(2:length(bb));nkn(ss)=1;else nkn(ss)=0;end

if k<nu+1,
	nbn(ss)=length(bb);nfn(ss)=length(ff)-1;ss=ss+1; 
	thb=[thb bb]; thf=[thf ff(2:length(ff))];
end
if k==nu+1, 
	ncn=length(bb)-1;ndn=length(ff)-1;
	thcc=bb(2:length(bb)); thd=ff(2:length(ff));
end

end
thh=[thb thcc thd thf];
if kder==0,
[ntr,ntc]=size(th);ntcc=max([ntc length(thh)]);
thc=zeros(3+length(thh),ntcc);
thc(1:ntr,1:ntc)=th;
thc(1,2)=-T;
if nu>0,thc(1,3:6+3*nu)=[nu 0 nbn ncn ndn nfn nkn];
else thc(1,3:6)=[nu 0 ncn ndn];end
thc(3,1:length(thh))=thh;thh0=thh;
end
if kder>0,
dif(:,kder)=(thh-thh0)'/dt(kder);
end
end

if ~isempty(dif),Pc=dif*P*dif';else Pc=zeros(length(thh),length(thh));end
thc(4:length(thh)+3,1:length(thh))=Pc;
if delays, thc(length(thh)+4,1:length(nk))=(max(nk,1)-1)*T;end

