% function out = daug(mat1,mat2,mat3,...,matN)
%
%   Diagonal augmentation of SYSTEM/VARYING/CONSTANT,
%   matrices, currently limited to 9 input matrices.
%
%            |  mat1  0     0    .   0   |
%            |   0   mat2   0    .   0   |
%      out = |   0    0    mat3  .   0   |
%            |   .    .     .    .   .   |
%            |   0    0     0    .  matN |
%
%   See also: ABV, MADD, MMULT, SBS, SEL, and VDIAG.

function out = daug(mat1,mat2,mat3,mat4,mat5,mat6,mat7,mat8,mat9)
 if nargin < 2
   disp('usage: out = daug(mat1,mat2,mat3,...,matN)')
   return
 elseif nargin == 2
   top = mat1;
   bot = mat2;
   if isempty(top)
     out = bot;
   elseif isempty(bot)
     out = top;
   else
     [toptype,toprows,topcols,topnum] = minfo(top);
     [bottype,botrows,botcols,botnum] = minfo(bot);
     if (toptype == 'syst' & bottype == 'vary') | ...
          (toptype == 'vary' & bottype == 'syst')
       error('DAUG of SYSTEM and VARYING not allowed')
       return
     end
     zur = zeros(toprows,botcols);
     zll = zeros(botrows,topcols);
     if toptype == 'vary'
       if bottype == 'vary'
%	 both varying
         code = indvcmp(top,bot);
         if code == 1
           [vdt,rpt,indvt] = vunpck(top);
           [vdb,rpb,indvb] = vunpck(bot);
           tmp = [reshape(vdt,toprows,topnum*topcols);...
                 zeros(botrows,topnum*topcols)];
           tmpp = reshape(tmp,topnum*(toprows+botrows),topcols);
           tmp = [zeros(toprows,botnum*botcols);...
                 reshape(vdb,botrows,botnum*botcols)];
           tmppp = reshape(tmp,botnum*(toprows+botrows),botcols);
	   out = zeros((toprows+botrows)*topnum+1,topcols+botcols+1);
	   out((toprows+botrows)*topnum+1,topcols+botcols+1) = inf;
	   out((toprows+botrows)*topnum+1,topcols+botcols) = topnum;
	   out(1:topnum,topcols+botcols+1) = indvt;
	   out(1:(toprows+botrows)*topnum,1:topcols+botcols) = [tmpp tmppp];
         else
           error('varying data is inconsistent')
           return
         end
       else
%	 varying constant
         [vdt,rpt,indvt] = vunpck(top);
         tmp = [reshape(vdt,toprows,topnum*topcols);...
               zeros(botrows,topnum*topcols)];
         tmpp = reshape(tmp,topnum*(toprows+botrows),topcols);
         tmp = kron(ones(topnum,1),[zur;bot]);
	 out = zeros((toprows+botrows)*topnum+1,topcols+botcols+1);
	 out((toprows+botrows)*topnum+1,topcols+botcols+1) = inf;
	 out((toprows+botrows)*topnum+1,topcols+botcols) = topnum;
	 out(1:topnum,topcols+botcols+1) = indvt;
	 out(1:(toprows+botrows)*topnum,1:topcols+botcols) = [tmpp tmp];
       end
     elseif toptype == 'cons'
       if bottype == 'vary'
%	 constant varying
         [vdb,rpb,indvb] = vunpck(bot);
         tmp = [zeros(toprows,botnum*botcols);...
               reshape(vdb,botrows,botnum*botcols)];
         tmpp = reshape(tmp,botnum*(toprows+botrows),botcols);
         tmp = kron(ones(botnum,1),[top;zll]);
	 out = zeros((toprows+botrows)*botnum+1,topcols+botcols+1);
	 out((toprows+botrows)*botnum+1,topcols+botcols+1) = inf;
	 out((toprows+botrows)*botnum+1,topcols+botcols) = botnum;
	 out(1:botnum,topcols+botcols+1) = indvb;
	 out(1:(toprows+botrows)*botnum,1:topcols+botcols) = [tmp tmpp];
       elseif bottype == 'cons'
         out = [top zur; zll bot];
       else
         [a,b,c,d] = unpck(bot);
         bb = [zeros(botnum,topcols) b];
         cc = [zeros(toprows,botnum) ; c];
         dd = [top zur ; zll d];
         out = pck(a,bb,cc,dd);
       end
     else
       if bottype == 'cons'
         [a,b,c,d] = unpck(top);
         bb = [b zeros(topnum,botcols)];
         cc = [c ; zeros(botrows,topnum)];
         dd = [d zur ; zll bot];
         out = pck(a,bb,cc,dd);
       else
         [at,bt,ct,dt] = unpck(top);
         [nrat,dum] = size(at);
         [ab,bb,cb,db] = unpck(bot);
         [nrab,dum] = size(ab);
         a = [at zeros(nrat,nrab) ; zeros(nrab,nrat) ab];
         b = [bt zeros(nrat,botcols) ; zeros(nrab,topcols) bb];
         c = [ct zeros(toprows,nrab) ; zeros(botrows,nrat) cb];
         d = [dt zur ; zll db];
         out = pck(a,b,c,d);
       end
     end
   end
 else
   exp = ['out=daug(daug('];
   for i=1:nargin-2
     exp = [exp 'mat' int2str(i) ','];
   end
   exp = [exp 'mat' int2str(nargin-1) '),mat' int2str(nargin) ');'];
   eval(exp);
 end
%
% Copyright MUSYN INC 1991,  All Rights Reserved
