function sysout = dprescl(sys,midomega,FSflag,BLTflag)
%DPRESCL  Diagonal state/uncertainty scaling which preserves
%   uncertain input/output map.
%
%   SYSOUT = DPRESCL(SYS)  yields a system whose input/output
%   properties (possibly uncertain) are the same as SYS.  The
%   numerical conditioning of SYSOUT is usually better than
%   that of SYS, improving the accuracy of additional computations
%   performed with SYSOUT.
%

% Copyright 2003-2004 The MathWorks, Inc.

%   XXXMINOR doesn't work for arrays, more clear with error messages
%   to alert the user of this.  Order of arguments needs to be thought
%   through: 

if nargin<3
   FSflag = 0;
   BLTflag = 1;
elseif nargin<4
   BLTflag = 1;
end
tublk = zeros(0,2);
szsys = size(sys);

switch class(sys)
case 'uss'
   if length(szsys)>2
      error('DPRESCL only applies to single systems.');
   end
   [M,Delta] = lftdata(sys);
   SN = pvget(sys,'StateName');
   Nx = length(SN);
   Ts = pvget(sys,'Ts');
   [Af,Bf,Cf,Df] = ssdata(M); % all 2-d doubles
   if Ts==0 & BLTflag==1
      if nargin==1 || isempty(midomega)
         eigA= eig(Af);
         idx = find(real(eigA)>0 & abs(imag(eigA))<1e-6*abs(real(eigA)));
         if isempty(idx)
            midomega = 1;
         else
            midomega = 1.2*max(real(eigA(idx)));
         end
      end
      [As,Bs,Cs,Ds] = ssdata(c2d(M,2/midomega,'tustin'));
   else
      [As,Bs,Cs,Ds] = ssdata(M);
   end
   [ublk,pflag] = pvget(sys,'OldBlkStructure');
   if pflag==1
      nblk = size(ublk,1);
      for i=1:nblk
         if ublk(i,2)==0
            if FSflag==1
               tublk = [tublk;abs(ublk(i,:))];
            else
               tublk = [tublk;ones(abs(ublk(i,1)),2)];
            end
         else
            tublk = [tublk;abs(ublk(i,:))];
         end
      end
   else
      error('XXX: Need OldBlkStructure surrogate');
   end
case 'ss'
   if length(szsys)>2
      error('DPRESCL only applies to single systems.');
   end
   SN = pvget(sys,'StateName');
   [Af,Bf,Cf,Df] = ssdata(sys);
   Ts = sys.Ts;
   Nx = size(Af,1);
   if Ts==0 & BLTflag==1
      if nargin==1 || isempty(midomega)
         midomega = 1;
      end
      [As,Bs,Cs,Ds] = ssdata(c2d(sys,2/midomega,'tustin'));
   else
      [As,Bs,Cs,Ds] = ssdata(sys);
   end
otherwise
   error([class(sys) ' is invalid for DPRESCL']);
end
            
if FSflag==1
   stateblk = [Nx 0];
else
   stateblk = ones(Nx,2);
end
blk = [stateblk;tublk;szsys([2 1])];
W = [As Bs;Cs Ds]; % states/perts/IO
[bounds,info] = mussv(W,blk,'s');
[dl,dr] = mussvunwrap(info);

W = dl*([Af Bf;Cf Df]/dr);
A = W(1:Nx,1:Nx);
B = W(1:Nx,Nx+1:end);
C = W(Nx+1:end,1:Nx);
D = W(Nx+1:end,Nx+1:end);
nC = norm(C);
nB = norm(B);
if nC>0 & nB>0
   rat = sqrt(nB/nC);
   B = B/rat;
   C = C*rat;
end

if FSflag==1
   SN = repmat({''},[Nx 1]);
else
   for i=1:Nx
      SN{i} = [SN{i} '(units modified)'];
   end
end

switch class(sys)
case 'uss'
   UmatData = lft(Delta,[D C;B A]);
   UmatData = UmatData([szsys(1)+1:end 1:szsys(1)],[szsys(2)+1:end 1:szsys(2)]);
   sysout = uss(UmatData,SN,Ts,pvget(sys,'dynamicsys'));
case 'ss'
   sysout = ss(A,B,C,D,sys);
   sysout.StateName = SN;
otherwise
   error('Invalid')
end