%   Script AUTOBEW.M
%
%   Längsbewegung eines Fahrzeugs
%   (Illustration der Bewegungsgleichung)
%
% J. Lunze
% 29.4.2002
% für 2. Auflage: 5.2.2007
% für 3. Auflage: 22.12.2011
% für 5. Auflage: 20.12.2019
%
%
clear
close all
%               PARAMETER
m=1000;         % kg  Masse
c=100;          % kg/s Widerstandsbeiwert
g=9.81;         % kgm/s/s  Erdbeschleunigung
                % Kräfte in N,   also Eingangsgröße u=Antriebskraft in N
                % Zeit in Sekunden
                % Weg in m
                % Geschwindigkeit in m/s,
                %         also Ausgangsgröße y=Geschwindigkeit in m/s
A=-c/m;
B=1/m;
C=1;
D=0;
Auto=ss(A, B, C, D);
%%    Eigenbewegung bei unterschiedlichen Anfangsgeschwindigkeiten
%
figure(1)
subplot(4,1,[1,2]);
for x0=5:5:25
    [Y, T]=initial(Auto, x0);
    plot(T, Y, 'b');
    hold on
end
latextitle('Eigenbewegung und erzwungene Bewegung des Fahrzeugs');
latexylabel('$$v$$ in $$\frac{\rm m}{\rm s}$$');
rotateY;
xleer;
axis([0 40 0 25]);
yticks([0 10 20]);
%
%  Erzwungene Bewegung bei unterschiedlicher Eingangsamplitude
subplot(4,1,[3,4]);
[Y, T]=step(Auto);
for u=1000:250:2000
    plot(T, u*Y, 'b');
    latextext(40.4, max(u*Y), num2str(u));
    hold on
end
%latextitle('Erzwungene Bewegung des Fahrzeugs');
latextext(41, 23, '$$u$$ in N:');
latexylabel('$$y$$ in $$\frac{\rm m}{\rm s}$$');
rotateY;
latexxlabel('$$t$$ in s');
axis([0 40 0 20]);
yticks([0 10 20]);
xticks([0 10 20 30 40 50]);


%% Überlagerung von erzwungener und Eigenbewegung
%
figure(2)
subplot(4,1,[1,2]);
[Y, T]=step(Auto);
Y=1250*Y;
for x0=0:5:25
    [Yfrei, T]=initial(Auto, x0, T);
    plot(T, Y+Yfrei, 'b');
    hold on
end
latextitle('Fahrzeugbewegung bei unterschiedlicher Anfangsgeschwindigkeit');
latexxlabel('$$t$$ in s');
latexylabel('$$y$$ in $$\frac{\rm m}{\rm s}$$');
rotateY;
axis([0 40 0 25]);
yticks([0 10 20]);
xticks([0 10 20 30 40 50]);


%%   Erweiterung des Modells um die Wirkung der Hangabtriebskraft d=alpha
%   zwei Eingangsgrößen: u1=u (Antriebskraft)
%                        u2=d (Neigungswinkel)
figure(3)
Berw=[1/m -g];
Derw=[0 0];
Autoerw=ss(A, Berw, C, Derw);
Ts=0.04;
Tend=40;
kend=Tend/Ts;
Td=[0:Ts:Tend];
U=zeros(Tend/Ts+1, 2);
U(:,1)=1250;
%                Steigung erhöht sich auf 6% und bleibt dort konstant
for k=1:200
    U(k+Tend/2/Ts,2)=0.06*k/200;
end
for k=200+Tend/2/Ts:kend+1
    U(k,2)=0.06;
end
Yd=lsim(Autoerw, U, Td);
%
subplot(4,1,[1,2]);
plot(Td, U(:,2),  'b')
latextitle('Bewegung des Fahrzeugs am Anstieg');
latexylabel('$$\alpha$$ in rad');
rotateY;
xleer;
axis([0 40 -0.01 0.07]);
yticks([0 0.03 0.06]);
%
subplot(4,1,[3,4]);
plot(Td, Yd, 'b')
hold on
plot(T, Y, 'b--');
latexxlabel('$$t$$ in s');
latexylabel('$$y$$ in $$\frac{\rm m}{\rm s}$$');;
rotateY;
axis([0 40 0 20]);
yticks([0 10 20 30]);
xticks([0 10 20 30 40]);


%%   Unterschiedliche Darstellung der Bewegung
%
figure(4)
Un=zeros(Tend/Ts+1, 1);
for k=1:kend/2
    Un(k)=1250;
end
x0=5;
[X, T]=initial(Auto, x0);  % Erzwungene + freie Bewegung
Yn=lsim(Auto, Un, Td, x0);
x02=Yn(kend/2);
X2=initial(Auto, x02, Td); % Neuberechnung als freie Bewegung
for k=1:kend/2-1
    Xneu(k)=0;
end
for k=1:kend/2+1
    Xneu(k+kend/2)=X2(k);
end
%
subplot(8,1,[1,2])
plot(T, X, 'b')
latextitle('Bewegung des Fahrzeugs am Anstieg');
latexylabel('$$y_{\rm frei}$$ in $$\frac{\rm m}{\rm s}$$');
rotateY;
xleer;
axis([0 40 -1 6]);
yticks([0 5]);
%
subplot(8,1,[3,4])
plot(Td, Un,  'b');
xleer;
latexylabel('$$u$$ in N');
rotateY;
axis([0 40 -100 1400]);
yticks([0 1000]);
%
subplot(8,1,[5,6])
plot(Td, Yn, 'b', [kend/2*Ts, kend/2*Ts], [x02, x02], 'bo');
%latexxlabel('$$t$$ in s');
xleer;
latexylabel('$$y$$ in $$\frac{\rm m}{\rm s}$$');
rotateY;
axis([0 40 -4 19]);
yticks([0 10]);
latextext(22, 15, '$$\tilde{x}_0$$');
%
subplot(8,1,[7,8])
plot(Td, Xneu, 'b')
latexxlabel('$$t$$ in s');
%ylabel('x','Interpreter','latex','FontName', 'Times', 'String','$$\tilde{y}$$ in m/s  ');
latexylabel('$$y$$ in $$\frac{\rm m}{\rm s}$$');
rotateY;
axis([0 40 -4 17]);
yticks([0 10]);
xticks([0 10 20 30 40 50]);

%%   EPS-Bilder
%
Epsbildc15('autobew');
