function th=thss2th(eta,iy)
%THSS2TH Converts an internal TH(SS)-format to the standard multi-input-
%	singe-output THETA-format. Auxiliary routine to th2ff and th2zp
%
%	TH = thss2th(TH_SS,IY)
%
%	TH_SS: The model defined in the TH(SS) format (See help th(ss))
%	IY:  The output number in TH_SS to be considered as the
%	     output in TH (The THETA-format handles only single-
%	     output models). Default IY=1
%	TH:  The corresponding model in the THETA-format including
%	     the transformed covariance matrix (see help theta)
%	
%	The tranformation of the covariance matrix is carried out using
%	the Gauss' approximation formula with numerical differentiation.
%	The stepsize in the differentiation is decided by the
%	m-file nuderst.

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

if nargin<2, iy=1;end
[nr,nc]=size(eta);
T=eta(1,2);nu=eta(1,3);ny=eta(1,4);

sspm=getmfth(eta);
nd=eta(1,5);
[etapar,p,lambda]=th2par(eta);if length(etapar)==0,etapar=0;end
dt=nuderst(etapar);
arg=getargth(eta);
if any(eta(2,8)==[2 3 4]) & T<0,flag=1;else flag=0;end
if eta(2,8)==3 & T<0,flaga=1;else flaga=0;end
if any(eta(2,8)==[2 3]),Tmod=-1;else Tmod=T;end
[a,b,c,d,k,x0]=feval(sspm,etapar,Tmod,arg);
if flag, [a,bk,c,d]=contth(a,[b k],c,d,abs(T),flaga);
b=bk(:,1:nu);k=bk(:,nu+1:nu+ny);end
[ny,nx]=size(c);
c=c(iy,:);
if nu>0,d=d(iy,:);else d=[];end
Apol=poly(a);

[nx,nu]=size(b);
for kk=1:nu
  Bpol(kk,:)=poly(a-b(:,kk)*c)+(d(kk)-1)*Apol;
  bs=(abs(Bpol(kk,:))>eps);
  Bpol(kk,:)=bs.*Bpol(kk,:);
end
Cpol=poly(a-k(:,iy)*c);
if T>0,ia=max(find(abs(Cpol)>eps));Cpol=Cpol(1:ia);end
if T>0,ia=max(find(abs(Apol)>eps));Apol=Apol(1:ia);end

Ncap=getncap(eta);lambb=lambda(iy,iy);
th=poly2th(Apol,Bpol,Cpol,1,ones(nu,1),lambb,T);
if ~isempty(Ncap),
   ndtemp=length(th2par(th));th(2,1)=lambb*(1+ndtemp/Ncap)/(1-ndtemp/Ncap);
else th(2,1)=lambb;
end
if nu>0
    na=th(1,4);nb=th(1,5:4+nu);nc=th(1,5+nu);
    ndd=th(1,6+nu);nf=th(1,7+nu:6+2*nu); nk=th(1,7+2*nu:6+3*nu);
else na=th(1,4);nb=0;nc=th(1,5);ndd=th(1,6);nf=0;nk=0;
end

if ~isempty(p), if norm(p)>0
  for kl=1:nd
      th1=etapar;th1(kl)=th1(kl)+dt(kl)/2;
      th2=etapar;th2(kl)=th2(kl)-dt(kl)/2;
      [A1,B1,C1,D1,K1,X1]=feval(sspm,th1,Tmod,arg);
      [A2,B2,C2,D2,K2,X2]=feval(sspm,th2,Tmod,arg);
      if flag [A1,BK1,C1,D1]=contth(A1,[B1 K1],C1,D1,abs(T),flaga);
	 [A2,BK2,C2,D2]=contth(A2,[B2 K2],C2,D2,abs(T),flaga);
	 if nu>0,B1=BK1(:,1:nu);B2=BK2(:,1:nu);end
	 K1=BK1(:,nu+1:nu+ny);K2=BK2(:,nu+1:nu+ny);
      end
      Apol1=poly(A1);Apol2=poly(A2);C1=C1(iy,:);C2=C2(iy,:);
      if nu>0,D1=D1(iy,:);D2=D2(iy,:);end
      for kk=1:nu
	B1pol(kk,:)=poly(A1-B1(:,kk)*C1)+(D1(kk)-1)*Apol1;
	B2pol(kk,:)=poly(A2-B2(:,kk)*C2)+(D2(kk)-1)*Apol2;
      end
      C1pol=poly(A1-K1(:,iy)*C1);C2pol=poly(A2-K2(:,iy)*C2);
      if na>0, diff(kl,1:na)=(Apol1(2:na+1)-Apol2(2:na+1))/dt(kl);end
      if nc>0, diff(kl,na+sum(nb)+1:na+nc+sum(nb))=... 
      (C1pol(2:nc+1)-C2pol(2:nc+1))/dt(kl);end

      sb=na;sf=na+sum(nb)+nc+ndd;
      for k=1:nu
        if nb(k)>0,diff(kl,sb+1:sb+nb(k))=...
           (B1pol(k,nk(k)+1:nk(k)+nb(k))-B2pol(k,nk(k)+1:nk(k)+nb(k)))/dt(kl);
        end
        sb=sb+nb(k);
   
      end
   end

   PTH=diff'*p*diff;
   [d,nd]=size(diff');
   if nd>0,th(4:d+3,1:d)=PTH;end
end,end