function [a,b,c,d,k,x0,da,db,dc,dd,dk,dx0]=eta2ss(eta)
%ETA2SS  Auxiliary routine to TH2SS
%
%	[A,B,C,D,K,X0]=eta2ss(TH)
%
%	A,B etc: The state space matrices corresponding to TH. 
%
%	with
%	[A,B,C,D,K,X0,dA,dB,dC,dD,dK,dX0]=eta2ss(TH)
%	
%	also the standard deviations dA, dB etc are obtained

%	L. Ljung 7-10-88
%	Copyright (c) 1988-90 The MathWorks, Inc.
%	All Rights Reserved

[nr,nc]=size(eta);
T=eta(1,2);

sspm=getmfth(eta);

nd=eta(1,5);
carg=eta(1,7);
etapar=eta(3,1:nd);if nd==0,etapar=0;end
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);
[nx,nu]=size(b);[ny,nx]=size(c);
if flag,
	[a,bk,c,d]=contth(a,[b k],c,d,abs(T),flaga);
	b=bk(:,1:eta(1,3));k=bk(:,eta(1,3)+1:eta(1,3)+eta(1,4));
end
input=1;if isempty(b),input=0;end
if nargout>6
  if nd==0,da=zeros(nx,nx);db=zeros(nx,nu);dc=zeros(ny,nx);dd=zeros(ny,nu);
           dk=zeros(nx,ny);dx0=zeros(nx,1);
  else   
  if (eta(2,8)==1 & T<0) | (any(eta(2,8)==[2 3]) & T>0)
      p=diag(eta(4:3+nd,1:nd));p=sqrt(p')+eps;
      etap1=etapar+p;
      [a1,b1,c1,d1,k1,x01]=feval(sspm,etap1,Tmod,arg);
      if flag, 
         [a1,bk1,c1,d1]=contth(a1,[b1 k1],c1,d1,abs(T),flaga);
	 b1=bk1(:,1:eta(1,3));k1=bk1(:,eta(1,3)+1:eta(1,3)+eta(1,4));
      end
      da=abs(a1-a);
      db=abs(b1-b);
      dc=abs(c1-c);
      dd=abs(d1-d);
      dk=abs(k1-k);
      dx0=abs(x01-x0);
  else
      [etapar,p]=th2par(eta);
      if isempty(p) | norm(p)==0,
      	  [nx,nu]=size(b);[ny,nx]=size(c);	
          da=zeros(nx,nx);db=zeros(nx,nu);dc=zeros(ny,nx);
	  dd=zeros(ny,nu);dk=zeros(nx,ny);dx0=zeros(nx,1);
      return,end
      dt=nuderst(etapar); 
      for kk=1:length(dt)
        etap1=etapar;etap1(kk)=etap1(kk)+dt(kk);
	[a1,b1,c1,d1,k1,x01]=feval(sspm,etap1,Tmod,arg);
	if flag, 
	  [a1,bk1,c1,d1]=contth(a1,[b1 k1],c1,d1,abs(T),flaga);
	  b1=bk1(:,1:eta(1,3));k1=bk1(:,eta(1,3)+1:eta(1,3)+eta(1,4));
	end
	da=a1-a;cola(:,kk)=da(:)/dt(kk);
	dc=c1-c;colc(:,kk)=dc(:)/dt(kk);
	dk=k1-k;colk(:,kk)=dk(:)/dt(kk);dx0=x01-x0;colx0(:,kk)=dx0(:)/dt(kk);
	if input
	   db=b1-b;colb(:,kk)=db(:)/dt(kk);dd=d1-d;cold(:,kk)=dd(:)/dt(kk);
        end
      end
      pa=cola*p*cola';da(:)=sqrt(diag(pa));
      pc=colc*p*colc';dc(:)=sqrt(diag(pc));
      pk=colk*p*colk';dk(:)=sqrt(diag(pk));
      px0=colx0*p*colx0';dx0(:)=sqrt(diag(px0));
      if input
        pb=colb*p*colb';db(:)=sqrt(diag(pb));
        pd=cold*p*cold';dd(:)=sqrt(diag(pd));
      else db=[];dd=[];
      end
   end
   end
end
