%   Script Taipeh.m
%
%   Demonstration der Wirkung von Nullstellen
%   am Beispiel des Schwingungsdämpfers für das
%   Hochhaus Taipeh-101
%
% J. Dastych
% J. Lunze
% 22.1.2010
% für 8. Auflage: 16.2.2010
% für 9. Auflage: 10.2.2012
% für 10. Auflage: 27.1.2014
% für 11. Auflage: 20.4.2016
% für 12. Auflage: 23.03.2020
%
%   
clear
close all
%
%    Parameter des Versuchsaufbaus
c1=100;   % N/m  Petersen, S.824:  m1/m2=0.08...0.1
c2=10;    % N/m
m1=1;     % kg
m2=0.11;   % kg
%   Dämpfung von Haus und Pendel: d1=d2=0
%   (dann heißt der Schwingungsdämpfer Tilger)
%
%                Zustandsraummodell
%    x1 = Position der Masse m1
%    x2 = Geschwindigkeit der Masse m1
%    x3 = Position der Masse m2
%    x4 = Geschwindigkeit der Masse m2
%    u  = Anregung durch Wind
%    y  = Position der Masse m1
n=4;
A=[0  1  0  0;
   -(c1+c2)/m1 0 c2/m1 0;
   0  0  0  1;
   c2/m2 0 -c2/m2 0];
b=[0; 1/m1; 0; 0];
c = [1  0  0  0];
d=0;
Hochhaus=ss(A, b, c, d);
%
%   ohne Tilger
nH=2;
AH=[0  1;  -c1/m1  0];
bH=[0; 1/m1];
cH=[1 0];
dH=0;
HausOhneTilger=ss(AH, bH, cH, dH);
% Nullstellen des Systems:
Nullstellen=tzero(Hochhaus)
% Eigenwerte des Systems;
Eigenwerte=eig(A)
% Anregung durch Sinusschwingung an den Nullstellen
% Frequenz= omega
AG=[0 abs(Nullstellen(1,1)); -abs(Nullstellen(1,1)) 0];
bG=[0; 0];
cG=[1 1];
dG=0;
Anregung=ss(AG, bG, cG, dG);
xG0=[1; -1];
Tend = 20;
T = [0:0.01:Tend];   % Zeitskala


%%    Hochhaus bei Erregung "in der Nullstelle"
%   
figure(1);
subplot(6,1,[1,2]);
U=initial(Anregung, xG0, T);
plot(T, U, 'b', [-5 Tend], [0 0], 'k:');
hold on;
plot([-5 0 0],[0 0 U(1,1)],'b');
latextitle('Hochhaus mit unged\"ampftem Schwingungstilger');
latexylabel('$$u$$ in N');
rotateY;
xleer;
axis([-2, Tend, -2.3, 2.3]);
yticks([-2 0 2]);
hold off;
% Bewegung des Hauses
subplot(6,1,[3,4]);
[Y1,T,X1]=lsim(Hochhaus, U, T);
plot(T, Y1, 'b');
hold on
plot([-5 Tend], [0 0], 'k:', [-5 0 0],[0 0 Y1(1,1)],'b');
latexylabel('$$y$$  in m');
rotateY;
axis([-2 Tend -0.1 0.1]);
xleer;
yticks([-0.1 0 0.1]);
hold off;
%  Bewegung des Pendels
subplot(6,1,[5,6]);
plot(T, X1(:,3), 'b');
hold on
plot([-5 Tend], [0 0], 'k:', [-5 0 0],[0 0 X1(1,3)],'b');
latexxlabel('$$t$$ in s');
latexylabel('$$x_3$$  in m')
rotateY;
axis([-2 Tend -0.5 0.5]);
xticks([0 5 10 15 20]);
yticks([-0.5 0.5]);
hold off;


%%  Vergleich mit Hochhaus ohne Schwingungstilger
%   
figure(2);
subplot(6,1,[1,2]);
[Y2,T,X2]=lsim(HausOhneTilger, U, T);
plot(T, Y2, 'b', T, Y1, 'b');
hold on
plot([-5 Tend], [0 0], 'k:', [-5 0 0],[0 0 Y2(1,1)]);
latexylabel('$$y$$  in m');
rotateY;
axis([-2 Tend -0.5 0.5]);
xticks([0 5 10 15 20]);
yticks([-0.5 0.5]);
latexxlabel('$$t$$ in s');
hold off;
 

%%    Haus mit Dämpfung; Vergleich des Amplitudengangs
%   
figure(3)
d1=0; % bisher: 0.015    % Petersen, S. 824:  d1=0,  d2=0.05
                             %   "optimale Einstellungen": d2=0.168
d2=0.02; % bisher: 0.015
Ad=[0  1  0  0;
   -(c1+c2)/m1 -(d1+d2)/m1 c2/m1 d2/m1;
   0  0  0  1;
   c2/m2 d2/m2 -c2/m2 -d2/m2];
bd=[0; 1/m1; 0; 0];
cd = [1  0  0  0];
dd=0;
Hochhausged=ss(Ad, bd, cd, dd);
%
%  Haus ohne Tilger
AHd=[0  1;  -c1/m1  d1/m1];
bHd=[0; 1/m1];
cHd=[1 0];
dHd=0;
HausOhneTilgerged=ss(AHd, bHd, cHd, dHd);
%
Omega=logspace(0, 2, 2000);
subplot(2,1,1)
GBetrag=bode(Hochhausged, Omega);
semilogx(Omega(:), 20*log(GBetrag(:)), 'b');
hold on
GBetragH=bode(HausOhneTilgerged, Omega);
semilogx(Omega(:), 20*log(GBetragH(:)),'b--');
latexxlabel('$$\omega$$ in  $$\frac{\rm rad}{\rm s}$$');
latexylabel('$$|G|$$ in dB');
rotateY;
grid('on');
axis([5, 20, -150 40]);
xticks([5 10 20]);
yticks([-150 -100 -50 0]);
latextitle('Vergleich der Amplitudeng\"ange');
hold off
 

%%  Gedämpftes Hochhaus bei Erregung "in der Nullstelle"
%   
figure(4);
% Anregung durch Sinusschwingung an den Nullstellen
% Frequenz= omega
Nullstellenged=tzero(Hochhausged)
AGd=[real(Nullstellenged(1,1))      abs(imag(Nullstellenged(1,1)));
    -abs(imag(Nullstellenged(1,1))) real(Nullstellenged(1,1))];
bGd=[0; 0];
cGd=[1 1];
dGd=0;
Anregungged=ss(AGd, bGd, cGd, dGd);
xG0=[1; -1];
subplot(6,1,[1,2]);
Ud=initial(Anregungged, xG0, T);
plot(T, Ud, 'b', [-5 Tend], [0 0], 'k:');
hold on;
plot([-5 0 0],[0 0 Ud(1,1)],'b');
xleer;
axis([-2, Tend, -2.3, 2.3]);
yticks([-2 0 2]);
latextitle('Hochhaus mit ged\"ampftem Schwingungstilger');
latexylabel('$$u$$ in N');
rotateY;
hold off;
% Bewegung des Hauses
subplot(6,1,[3,4]);
[Y1d,T,X1d]=lsim(Hochhausged, Ud, T);
plot(T, Y1d, 'b');
hold on
plot([-5 Tend], [0 0], 'k:', [-5 0 0],[0 0 Y1d(1,1)],'b');
latexylabel('$$y$$ in m');
rotateY;
axis([-2 Tend -0.1 0.1]);
xleer;
yticks([-0.1 0 0.1]);
hold off;
%  Bewegung des Pendels
subplot(6,1,[5,6]);
plot(T, X1d(:,3), 'b');
hold on
plot([-5 Tend], [0 0], 'k:', [-5 0 0],[0 0 X1d(1,3)],'b');
latexxlabel('$$t$$ in s');
latexylabel('$$x_3$$ in m');
rotateY;
axis([-2 Tend -0.5 0.5]);
xticks([0 5 10 15 20]);
yticks([-0.5 0 0.5]);
hold off;


%%   Vergleich mit Hochhaus ohne Schwingungstilger
%   
figure(5);
subplot(6,1,[1,2]);
[Y2d,T,X2d]=lsim(HausOhneTilgerged, Ud, T);
plot(T, Y2d, 'b', T, Y1d, 'b');
hold on
plot([-5 Tend], [0 0], 'k:', [-5 0 0],[0 0 Y2d(1,1)],'b');
latexylabel('$$y$$ in m');
rotateY;
axis([-2 Tend -0.5 0.5]);
xticks([0 5 10 15 20]);
yticks([-0.5 0 0.5]);
latexxlabel('$$t$$ in s');
hold off;



%% Herstellung der Bilddateien *.eps
%
Epsbildc15('taipeh');
