function eta=mketaarx(nn,par,lam,Tsamp)
%MKETAARX constructs an ETA-model structure for (multivariate) ARX-models
%	Auxiliary routine to ARX2TH
%	TH = mketaarx([NA NB NK])
%
%	TH: The resulting structure in the THETA-format (see HELP THETA)
%	NA: an ny|ny matrix whose i-j entry is the order of the
%	    polynomial (in the delay operator) that relates the
%	    j:th output to the i:th output
%	NB: an ny|nu matrix whose i-j entry is the order of the 
%	    polynomial that related the j:th input to the i:th
%	    output
%	NK: an ny|nu matrix whose i-j entry is the delay from
%	    the j:th input to the i:th output
%	(ny: the number of outputs, nu: the number of inputs)
%	
%	Each row of [NA NB NK] is thus consistent with the way
%	single output ARX structures is defined (See Help ARX)
%	
%	With TH = mketaarx([NA NB NK], PAR, LAM, T) the nominal values
%	of the free parameters are set to PAR (default zeros). LAM will
%       be the innovations covariance matrix, and T the sampling interval.

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

if nargin<4, Tsamp=[];end
if nargin<3, lam=[];end
if nargin<2, par=[];end

[ny,nc]=size(nn);
if isempty(Tsamp),Tsamp=1;end
if isempty(lam),lam=eye(ny);end
nu=(nc-ny)/2;
na=nn(:,1:ny);
if nu>0,nb=nn(:,ny+1:ny+nu);else nb=zeros(ny,1);end
if nu>0,nk=nn(:,ny+nu+1:nc);else nk=ones(ny,1);end
mna=max(max(na)');
nb1=nb+(nk-1);
mnb=max(max(nb1)');

mnk=min(min(nk)');
nx=mna*ny+mnb*nu;
C=zeros(ny,nx);
for ky=1:ny
for kr=1:ny
for kc=1:na(ky,kr)
C(ky,(kc-1)*ny+kr)=NaN;
end
end
for kr=1:nu
for kc=max(1,nk(ky,kr)):nb1(ky,kr)
C(ky,(kc-1)*nu+mna*ny+kr)=NaN;
end
end
end

if mna>0
A0=C;A1=[eye(ny*(mna-1)),zeros(ny*(mna-1),ny+nu*mnb)];
else A0=[];A1=[];
end
if mnb>0
A3=[zeros(nu,nx);[zeros(nu*(mnb-1),ny*mna),eye(nu*(mnb-1)),zeros(nu*(mnb-1),nu)]];
else A3=[];
end
A=[A0;A1;A3];
if nu>0
B=zeros(nx,nu);D=zeros(ny,nu);
bc=(nk==0).*(nb>0);
for kbc=1:ny
for kbr=1:nu
if bc(kbc,kbr), B(kbc,kbr)=NaN;end
end
end
D=B(1:ny,1:nu);
if mnb>0,B(ny*mna+1:ny*mna+nu,:)=eye(nu);end
else B=[];end
if nx==0, B=zeros(ny,nu);A=zeros(ny,ny);C=A;end
if nu>0,if mna>0,D=zeros(ny,nu);end,else D=[];end
K=zeros(nx,ny);

if mna>0;C=zeros(ny,nx);K(1:ny,1:ny)=eye(ny);end
x0=zeros(nx,1);if nx==0,x0=zeros(ny,1);end
ms=modstruc(A,B,C,D,K,x0);
if isempty(par), par=zeros(1,sum(sum(na)')+sum(sum(nb)'));end
[rms,cms]=size(ms);ms(1,cms)=-ms(1,cms);
eta=ms2th(ms,'d',par,lam,Tsamp);
eta(2,7)=39;eta(2,8)=3;

