function [bc,fc]=contin(th,ku,delays)
%CONTIN   converts a model to continuous-time form. USE RATHER THD2THC!
%
%	[BC,FC] = contin(TH)
%
%	TH: A discrete time model defined in the format described by HELP THETA
%
%	BC and FC are returned as the numerator and denominator polynomials,
%	respectively, of a continuous-time counterpart. Rows k contain
%	the model associated with the k:th input.
%
%	With [BC,FC] = contin(TH,KU) the continuous-time models associated with
%	the input numbers given by the entries of the row vector KU are
%	computed (default is all the inputs). To obtain a continuous-time
%	version of the noise model, let KU include the value zero. The noise 
%	model will then be given in the corresponding row of [BC,FC].
%
%	When the model TH contains extra delays (nk>1) these are approximated 
%	when forming the continuous time model. With a third, positive argument
%	[BC,FC]=contin(TH,KU,1), the extra delays are removed before forming 
%	the model. They should then be appended to the continuous model as a
%	dead-time.

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

T=th(1,2); nu=th(1,3);nk=th(1,7+2*nu:6+3*nu);
nb=th(1,5:4+nu);nf=th(1,7+nu:6+2*nu);nc=th(1,5+nu);nd=th(1,6+nu);
% *** Set up default values ***

if nargin<3,delays=0;end
if nargin<2, ku=1:nu;end
if length(ku)==1,if ku<0,ku=1:nu;end,end
if delays<0,delays=0;end


[a,b,c,d,f]=polyform(th);

if length(find(ku==0))>0
     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)=nd;	
end

ss=1;
for k=ku
   if k==0, k=nu+1;end
   den=conv(a,f(k,1:nf(k)+1)); 
   if delays~=0,nkmod=max(nk(k),1);else nkmod=1;end %Corr 89-07-27
   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;
    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;
    bc(ss,1:length(bb))=bb;
    fc(ss,1:length(ff))=ff;
    ss=ss+1;
end

