% test_mor
%
% ---------------------------------------------------------------------
%  T. Stykel, TU Berlin, 2009-01-12


clear all;
k   = 200;
n   = 2*k-1;
eta = 1/2;
tol = 1e-6;

% Systemmatrizen
[A, B, C] = stab(k, eta);

% Balanciertes Abschneiden
tic;
[Ar1, Br1, Cr1, hsv, errb] = btmor(A, B, C, tol);
toc
r = size(Ar1,1)

% Arnoldi-Verfahren, s0=0 
tic;
[Ar2, Br2, Cr2] = arnmor(A, B, C, 0, r);
toc

% Arnoldi-Verfahren, s0=1000*i
tic;
[Ar3, Br3, Cr3] = arnmor(A, B, C, 1000*i, r);
toc

% Frequenzg"ange G(iw), Gr1(iw), Gr2(iw)
wmin  = -3;
wmax  = 8;
numw  = 100;
wh    = (wmax-wmin)/numw;
w     = logspace(wmin,wmax,numw);
i     = sqrt(-1);
for j = 1:numw
   iwj = i*w(j);
   G   = C*((iwj*speye(n)-A)\B);
   Gr1 = Cr1*((iwj*eye(r)-Ar1)\Br1);
   Gr2 = Cr2*((iwj*eye(r)-Ar2)\Br2);
   Gr3 = Cr3*((iwj*eye(r)-Ar3)\Br3);
   
   nrmG(j) = norm(full(G));
   nrmGr1(j) = norm(full(Gr1));
   nrmGr2(j) = norm(full(Gr2));
   nrmGr3(j) = norm(full(Gr3));
   
   errGr1(j) = norm(full(G-Gr1));
   errGr2(j) = norm(full(G-Gr2));
   errGr3(j) = norm(full(G-Gr3));
end

figure(1); clf;
semilogy(1:length(hsv), hsv, '-*');
xlabel('j');
xlim([1 length(hsv)]);
title('Hankel-Singul"arwerte');

figure(2); clf;
semilogx(w,nrmG(:),'k', w,nrmGr1(:),'b-.', w, nrmGr2(:),'r--',w, nrmGr3(:),'m:','LineWidth',2);
%loglog(w, nrmG(:), 'k', w, nrmGr1(:), 'b-.', w, nrmGr2(:),'r:','LineWidth', 1.5);
xlabel('Frequenz w');
ylabel('Magnitude');
title('Frequenzg"ange');
legend('Originalsystem', 'BalTrunc', 'Arnoldi, s_0=0','Arnoldi, s_0=1000i');

figure(3); clf;
loglog(w, errb*ones(numw,1), 'k', w, errGr1(:), 'b-.', w, errGr2(:),'r--',w, errGr3(:),'m:','LineWidth', 2);
xlabel('Frequenz w');
ylabel('Magnitude');
title('Fehler und Fehlerschranke');
legend('Fehlerschranke', 'BalTrunc', 'Arnoldi, s_0=0', 'Arnoldi, s_0=1000i');



