function th=poly2th(a,b,c,d,f,LAM,T)
%POLY2TH constructs a "theta-matrix" from given polynomials
%
%	TH = poly2th(A,B,C,D,F,LAM,T)
%
%	TH: returned as a matrix of the standard format, describing the model
%	A(q) y(t) = [B(q)/F(q)] u(t-nk) + [C(q)/D(q)] e(t)
%	The exact format is given in HELP THETA.
%
%	A,B,C,D and F are entered as the polynomials. A,C,D and F start with
%	1 , while B contains leading zeros to indicate the delay(s) for discr.
%	time models. For multi-input systems B and F are matrices with the 
%	number of rows equal to the number of inputs. For a time series, B and
%	F are entered as []. POLY2TH is thus the inverse of TH2POLY.
%
%	LAM is the variance of the noise term e, and T is the sampling interval
%	A negative value of T indicates a time-continuous model. Then the poly-
%	nomials are entered in descending powers of s. Example: A=1,B=[1 2;0 3]
%	C=1;D=1;F=[1 0;0 1]; T=-1 corresponds to the time-continuous system
%	Y = (s+2)/s U1 + 3 U2.
% 
%	Trailing C,D,F, LAM, and T can be omitted, in which case they are
%	taken as 1's (if B=[], then F=[]).

%	L. Ljung 10-1-86
%       Revised 4-21-91
%	Copyright (c) 1986-90 by the MathWorks, Inc.
%	All Rights Reserved.

[nu,nb1]=size(b); nu=nu(1);

if nargin<7, T=1;end
if nargin<6, LAM=1;end
if nargin<5, f=ones(nu,1);end
if nargin<4, d=1;end
if nargin<3, c=1;end
if a(1)~=1, error('The A-polynomial should be monic (start with "1")'),end 
if T>0,for ku=1:nu,if f(ku,1)~=1,error('Each of the F-polynomials should be monic (start with "1") for discrete time models!'),end,end,end 
na=length(a)-1;nc=length(c)-1;nd=length(d)-1;
if nu>0,ib=b~=0; if nb1>1,nk=sum(cumsum(ib')==0);else nk=(cumsum(ib')==0);end
	if T>0,	ib=(b(:,nb1:-1:1)~=0); 
	if nb1>1,nb=-sum(cumsum(ib')==0)-nk+nb1;else nb=-(cumsum(ib')==0)-nk+nb1;end
	nb=max(nb,zeros(1,nu));
	else nb=nb1-nk;end,end
if nu==0, nb=0;nk=0;end
[nuf,nf1]=size(f);
if nuf~=nu, error('f and b must have the same number of rows!'),return,end
if nf1==1 | nf1==0
   nf=zeros(1,nu);
else
   if T>0
      ih=(f(:,nf1:-1:1)~=0); nf=-sum(cumsum(ih')==0)+nf1-1;
   else 
      ih=f~=0;
      for ku=1:nu 
        nf(ku)=nf1-min(find(ih(ku,:)~=0));
        if f(ku,nf1-nf(ku))~=1
         error('For continuous time systems the first non zero element of the rows of f must be a "1", corresponding to a monic F-polynomial')
        end
      end
   end
end
n=na+sum(nb)+nc+nd+sum(nf);
th=zeros(3,max([7 6+3*nu n]));
if nu>0
    th(1,1:6+3*nu)=[LAM T nu na nb nc nd nf nk];
else
    th(1,1:6)=[LAM T nu na nc nd];
end

if na>0, th(3,1:na)=a(2:na+1);end
 if nc>0, th(3,na+sum(nb)+1:na+nc+sum(nb))=c(2:nc+1);end
if nd>0, th(3,na+nc+1+sum(nb):na+nc+nd+sum(nb))=d(2:nd+1);end

sb=na;sf=na+sum(nb)+nc+nd;
for k=1:nu
   if nb(k)>0,th(3,sb+1:sb+nb(k))=b(k,nk(k)+1:nk(k)+nb(k));end
   if nf(k)>0,
     if T>0,th(3,sf+1:sf+nf(k))=f(k,2:nf(k)+1);
     else   th(3,sf+1:sf+nf(k))=f(k,nf1-nf(k)+1:nf1);
     end
   end
   sb=sb+nb(k);
   sf=sf+nf(k);
end
ti=fix(clock); ti(1)=ti(1)/100;
th(2,2:6)=ti(1:5); th(2,7)=6;
if T<0 & length(c)>length(conv(d,a)),disp('This model has differentiated white measurement noise. Expect problems!'),end
if T<0 & length(c)<length(conv(d,a))
disp('This model has no white  component in the measurement noise. Expect problems with sampling and state space representations!'),end
