--- /dev/null
+## Copyright (C) 2011 Lukas F. Reichlin
+##
+## This file is part of LTI Syncope.
+##
+## LTI Syncope is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## LTI Syncope is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{scaledsys}, @var{info}] =} prescale (@var{sys})
+## Prescale state-space model.
+## Frequency response commands perform automatic scaling unless model property
+## @var{scaled} is set to @var{true}.
+##
+## @strong{Inputs}
+## @table @var
+## @item sys
+## LTI model.
+## @end table
+##
+## @strong{Outputs}
+## @table @var
+## @item scaledsys
+## Scaled state-space model.
+## @item info
+## Structure containing additional information.
+## @item info.SL
+## Left scaling factors. @code{Tl = diag (info.SL)}.
+## @item info.SR
+## Right scaling factors. @code{Tr = diag (info.SR)}.
+## @end table
+##
+## @strong{Equations}
+## @example
+## @group
+## Es = Tl * E * Tr
+## As = Tl * A * Tr
+## Bs = Tl * B
+## Cs = C * Tr
+## Ds = D
+## @end group
+## @end example
+##
+## For proper state-space models, @var{Tl} and @var{Tr} are inverse of each other.
+##
+## @strong{Algorithm}@*
+## Uses SLICOT TB01ID and TG01AD by courtesy of
+## @uref{http://www.slicot.org, NICONET e.V.}.
+## @end deftypefn
+
+## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
+## Created: June 2011
+## Version: 0.2
+
+function [retsys, varargout] = prescale (sys)
+
+ if (nargin != 1)
+ print_usage ();
+ endif
+
+ if (! isa (sys, "ss"))
+ warning ("prescale: system not in state-space form");
+ sys = ss (sys);
+ endif
+
+ [retsys, lscale, rscale] = __prescale__ (sys);
+
+ if (nargout > 1)
+ varargout{1} = struct ("SL", lscale, "SR", rscale);
+ endif
+
+endfunction
+
+
+## Scaling of state-space models, direct access to TB01ID
+%!shared Ao, Bo, Co, SCALEo, MAXREDo, Ae, Be, Ce, SCALEe, MAXREDe
+%! A = [ 0.0 1.0000e+000 0.0 0.0 0.0
+%! -1.5800e+006 -1.2570e+003 0.0 0.0 0.0
+%! 3.5410e+014 0.0 -1.4340e+003 0.0 -5.3300e+011
+%! 0.0 0.0 0.0 0.0 1.0000e+000
+%! 0.0 0.0 0.0 -1.8630e+004 -1.4820e+000 ];
+%!
+%! B = [ 0.0 0.0
+%! 1.1030e+002 0.0
+%! 0.0 0.0
+%! 0.0 0.0
+%! 0.0 8.3330e-003 ];
+%!
+%! C = [ 1.0000e+000 0.0 0.0 0.0 0.0
+%! 0.0 0.0 1.0000e+000 0.0 0.0
+%! 0.0 0.0 0.0 1.0000e+000 0.0
+%! 6.6640e-001 0.0 -6.2000e-013 0.0 0.0
+%! 0.0 0.0 -1.0000e-003 1.8960e+006 1.5080e+002 ];
+%!
+%! MAXRED = 0.0;
+%!
+%! [Ao, Bo, Co, MAXREDo, SCALEo] = sltb01id (A, B, C, MAXRED);
+%!
+%! Ae = [ 0.0000000D+00 0.1000000D+05 0.0000000D+00 0.0000000D+00 0.0000000D+00
+%! -0.1580000D+03 -0.1257000D+04 0.0000000D+00 0.0000000D+00 0.0000000D+00
+%! 0.3541000D+05 0.0000000D+00 -0.1434000D+04 0.0000000D+00 -0.5330000D+03
+%! 0.0000000D+00 0.0000000D+00 0.0000000D+00 0.0000000D+00 0.1000000D+03
+%! 0.0000000D+00 0.0000000D+00 0.0000000D+00 -0.1863000D+03 -0.1482000D+01 ];
+%!
+%! Be = [ 0.0000000D+00 0.0000000D+00
+%! 0.1103000D+04 0.0000000D+00
+%! 0.0000000D+00 0.0000000D+00
+%! 0.0000000D+00 0.0000000D+00
+%! 0.0000000D+00 0.8333000D+02 ];
+%!
+%! Ce = [ 0.1000000D-04 0.0000000D+00 0.0000000D+00 0.0000000D+00 0.0000000D+00
+%! 0.0000000D+00 0.0000000D+00 0.1000000D+06 0.0000000D+00 0.0000000D+00
+%! 0.0000000D+00 0.0000000D+00 0.0000000D+00 0.1000000D-05 0.0000000D+00
+%! 0.6664000D-05 0.0000000D+00 -0.6200000D-07 0.0000000D+00 0.0000000D+00
+%! 0.0000000D+00 0.0000000D+00 -0.1000000D+03 0.1896000D+01 0.1508000D-01 ];
+%!
+%! SCALEe = [0.1000000D-04 0.1000000D+00 0.1000000D+06 0.1000000D-05 0.1000000D-03 ];
+%!
+%! MAXREDe = 0.3488E+10;
+%!
+%!assert (Ao, Ae, 1e-4);
+%!assert (Bo, Be, 1e-4);
+%!assert (Co, Ce, 1e-4);
+%!assert (MAXREDo, MAXREDe, 1e6);
+%!assert (SCALEo, SCALEe.', 1e-4);
+
+
+## Scaling of descriptor state-space models, direct access to TG01AD
+%!shared Ao, Eo, Bo, Co, LSCALEo, RSCALEo, Ae, Ee, Be, Ce, LSCALEe, RSCALEe
+%! A = [ -1 0 0 0.003
+%! 0 0 0.1000 0.02
+%! 100 10 0 0.4
+%! 0 0 0 0.0];
+%!
+%! E = [ 1 0.2 0 0.0
+%! 0 1 0 0.01
+%! 300 90 6 0.3
+%! 0 0 20 0.0];
+%!
+%! B = [ 10 0
+%! 0 0
+%! 0 1000
+%! 10000 10000];
+%!
+%! C = [ -0.1 0.0 0.001 0.0
+%! 0.0 0.01 -0.001 0.0001];
+%!
+%! TRESH = 0.0;
+%!
+%! [Ao, Eo, Bo, Co, LSCALEo, RSCALEo] = sltg01ad (A, E, B, C, TRESH);
+%!
+%! Ae = [ -1.0000 0.0000 0.0000 0.3000
+%! 0.0000 0.0000 1.0000 2.0000
+%! 1.0000 0.1000 0.0000 0.4000
+%! 0.0000 0.0000 0.0000 0.0000];
+%!
+%! Ee = [ 1.0000 0.2000 0.0000 0.0000
+%! 0.0000 1.0000 0.0000 1.0000
+%! 3.0000 0.9000 0.6000 0.3000
+%! 0.0000 0.0000 0.2000 0.0000 ];
+%!
+%! Be = [100.0000 0.0000
+%! 0.0000 0.0000
+%! 0.0000 100.0000
+%! 100.0000 100.0000 ];
+%!
+%! Ce = [ -0.0100 0.0000 0.0010 0.0000
+%! 0.0000 0.0010 -0.0010 0.0010];
+%!
+%! LSCALEe = [ 10.0000 10.0000 0.1000 0.0100 ];
+%!
+%! RSCALEe = [ 0.1000 0.1000 1.0000 10.0000 ];
+%!
+%!assert (Ao, Ae, 1e-4);
+%!assert (Eo, Ee, 1e-4);
+%!assert (Bo, Be, 1e-4);
+%!assert (Co, Ce, 1e-4);
+%!assert (LSCALEo, LSCALEe.', 1e-4);
+%!assert (RSCALEo, RSCALEe.', 1e-4);
+
+
+## Scaling of state-space models, user function
+%!shared Ao, Bo, Co, INFOo, Ae, Be, Ce, SCALEe
+%! A = [ 0.0 1.0000e+000 0.0 0.0 0.0
+%! -1.5800e+006 -1.2570e+003 0.0 0.0 0.0
+%! 3.5410e+014 0.0 -1.4340e+003 0.0 -5.3300e+011
+%! 0.0 0.0 0.0 0.0 1.0000e+000
+%! 0.0 0.0 0.0 -1.8630e+004 -1.4820e+000 ];
+%!
+%! B = [ 0.0 0.0
+%! 1.1030e+002 0.0
+%! 0.0 0.0
+%! 0.0 0.0
+%! 0.0 8.3330e-003 ];
+%!
+%! C = [ 1.0000e+000 0.0 0.0 0.0 0.0
+%! 0.0 0.0 1.0000e+000 0.0 0.0
+%! 0.0 0.0 0.0 1.0000e+000 0.0
+%! 6.6640e-001 0.0 -6.2000e-013 0.0 0.0
+%! 0.0 0.0 -1.0000e-003 1.8960e+006 1.5080e+002 ];
+%!
+%! SYS = ss (A, B, C);
+%!
+%! [SYSo, INFOo] = prescale (SYS);
+%!
+%! [Ao, Bo, Co] = ssdata (SYSo);
+%!
+%! Ae = [ 0.0000000D+00 0.1000000D+05 0.0000000D+00 0.0000000D+00 0.0000000D+00
+%! -0.1580000D+03 -0.1257000D+04 0.0000000D+00 0.0000000D+00 0.0000000D+00
+%! 0.3541000D+05 0.0000000D+00 -0.1434000D+04 0.0000000D+00 -0.5330000D+03
+%! 0.0000000D+00 0.0000000D+00 0.0000000D+00 0.0000000D+00 0.1000000D+03
+%! 0.0000000D+00 0.0000000D+00 0.0000000D+00 -0.1863000D+03 -0.1482000D+01 ];
+%!
+%! Be = [ 0.0000000D+00 0.0000000D+00
+%! 0.1103000D+04 0.0000000D+00
+%! 0.0000000D+00 0.0000000D+00
+%! 0.0000000D+00 0.0000000D+00
+%! 0.0000000D+00 0.8333000D+02 ];
+%!
+%! Ce = [ 0.1000000D-04 0.0000000D+00 0.0000000D+00 0.0000000D+00 0.0000000D+00
+%! 0.0000000D+00 0.0000000D+00 0.1000000D+06 0.0000000D+00 0.0000000D+00
+%! 0.0000000D+00 0.0000000D+00 0.0000000D+00 0.1000000D-05 0.0000000D+00
+%! 0.6664000D-05 0.0000000D+00 -0.6200000D-07 0.0000000D+00 0.0000000D+00
+%! 0.0000000D+00 0.0000000D+00 -0.1000000D+03 0.1896000D+01 0.1508000D-01 ];
+%!
+%! SCALEe = [0.1000000D-04 0.1000000D+00 0.1000000D+06 0.1000000D-05 0.1000000D-03 ];
+%!
+%!assert (Ao, Ae, 1e-4);
+%!assert (Bo, Be, 1e-4);
+%!assert (Co, Ce, 1e-4);
+%!assert (INFOo.SL.^-1, SCALEe.', 1e-4);
+%!assert (INFOo.SR, SCALEe.', 1e-4);
+
+
+## Scaling of descriptor state-space models, user function
+%!shared Ao, Eo, Bo, Co, INFOo, Ae, Ee, Be, Ce, LSCALEe, RSCALEe
+%! A = [ -1 0 0 0.003
+%! 0 0 0.1000 0.02
+%! 100 10 0 0.4
+%! 0 0 0 0.0];
+%!
+%! E = [ 1 0.2 0 0.0
+%! 0 1 0 0.01
+%! 300 90 6 0.3
+%! 0 0 20 0.0];
+%!
+%! B = [ 10 0
+%! 0 0
+%! 0 1000
+%! 10000 10000];
+%!
+%! C = [ -0.1 0.0 0.001 0.0
+%! 0.0 0.01 -0.001 0.0001];
+%!
+%! SYS = dss (A, B, C, [], E);
+%!
+%! [SYSo, INFOo] = prescale (SYS);
+%!
+%! [Ao, Bo, Co, ~, Eo] = dssdata (SYSo);
+%!
+%! Ae = [ -1.0000 0.0000 0.0000 0.3000
+%! 0.0000 0.0000 1.0000 2.0000
+%! 1.0000 0.1000 0.0000 0.4000
+%! 0.0000 0.0000 0.0000 0.0000];
+%!
+%! Ee = [ 1.0000 0.2000 0.0000 0.0000
+%! 0.0000 1.0000 0.0000 1.0000
+%! 3.0000 0.9000 0.6000 0.3000
+%! 0.0000 0.0000 0.2000 0.0000 ];
+%!
+%! Be = [100.0000 0.0000
+%! 0.0000 0.0000
+%! 0.0000 100.0000
+%! 100.0000 100.0000 ];
+%!
+%! Ce = [ -0.0100 0.0000 0.0010 0.0000
+%! 0.0000 0.0010 -0.0010 0.0010];
+%!
+%! LSCALEe = [ 10.0000 10.0000 0.1000 0.0100 ];
+%!
+%! RSCALEe = [ 0.1000 0.1000 1.0000 10.0000 ];
+%!
+%!assert (Ao, Ae, 1e-4);
+%!assert (Eo, Ee, 1e-4);
+%!assert (Bo, Be, 1e-4);
+%!assert (Co, Ce, 1e-4);
+%!assert (INFOo.SL, LSCALEe.', 1e-4);
+%!assert (INFOo.SR, RSCALEe.', 1e-4);