function [A,B,dA,dB]=th2arx(eta)
%TH2ARX converts a THETA-format model to an ARX-model
%
%	[A,B]=th2arx(TH)
%
%	TH: The model structure defined in the THETA-format (See help theta)
%	A, B : Matrices defining the ARX-structure:
%
%	       y(t) + A1 y(t-1) + .. + An y(t-n) = 
%	       = B0 u(t) + ..+ B1 u(t-1) + .. Bm u(t-m)
%
%	       A = [I A1 A2 .. An],  B=[B0 B1 .. Bm]
%
%
%	With [A,B,dA,dB] = th2arx(TH), also the standard deviations
%	of A and B, i.e. dA and dB are computed.
%
%	See also ARX2TH, and ARX

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

[Ncap,nd]=getncap(eta);
if ~isthss(eta), [A,B]=th2poly(eta);
if nargout>2,[par,P]=th2par(eta);eta1=eta;eta1(3,1:nd)=par+sqrt(diag(P))';
[A1,B1]=th2poly(eta1);dA=abs(A-A1);dB=abs(B-B1);end,return,end

tnr=eta(2,8); if tnr~=3,error('This is not an ARX-model!'),end
etapar=eta(3,1:nd);
arg=getargth(eta);
[rarg,carg]=size(arg);
as1=arg(:,1:rarg);
nr=find(sum(as1')==0);end
[as,bs,cs,ds]=ssmodx9(etapar,-1,arg);%FIXA
[ny,nz]=size(cs);[nx,nz]=size(as);[nz,nu]=size(bs);
if isempty(nr),na=nx/ny;else na=(nr(1)-1)/ny;end
if nu>0,nb=(nx-na*ny)/nu;else nb=0;end
A=[eye(ny) -cs(:,1:na*ny)];
if nu>0,B=[ds cs(:,na*ny+1:nx)];else B=[];end
if na==0 & ~any(isnan(arg(:,nx+nu+1:nx+nu+ny))), B=ds;end

if nargout>2
p=diag(eta(4:3+nd,1:nd));p=sqrt(p')+eps;
etap1=etapar+p;
[das,dbs,dcs,dds]=ssmodx9(etap1,-1,arg);
dA=abs([eye(ny) -dcs(:,1:na*ny)]-A);
if nu>0,dB=abs([dbs(1:ny,:) dcs(:,na*ny+1:nx)]-B);else dB=[];end
end

