%
% Phasenportraet einer gedaempften Schwingung
% mit Vektorfeld und einer Trajektorie
%
% PT2-Glied
printf("Parameter fuer das PT2-Glied:\n");
d=0.15
%T = 1./ (2. * pi)
T = 1.
K = 1.
% Impulsantwort auf einen Dirac-Stoß
% Nach Wiki PT2-Glied
% K/( sqrt(1-d^2) ) * exp( -1. * (d / T) * t ) * sin( sqrt(1-d^2)/T * t + arctan( sqrt(1-d^2)/d ) )
% Nach Wiki linear gedämpfte Schwingung
% d = delta / omega_0
% T = 1. / omega_0
printf("Umrechnung fuer die linear gedaempfte Schwingung:\n");
delta = d / T
omega_0 = 1. / T
omega_d = sqrt(omega_0^2 - delta^2)
printf("Startwerte:\n");
x_0 = 1;
v_0 = 0;
printf("Berechnung:\n");
perioden = 6.
werte_pro_periode = 40.
periodendauer = 2.*pi/omega_0
t = 0:periodendauer/werte_pro_periode:periodendauer*perioden;
% Formel: x = e^(-1. * delta * t) * ( (v_0 + delta * x_0) / omega_d * sin(omega_d * t) + x_0 * cos(omega_d * t) );
printf("Berechne Zeitverhalten\n");
part_1 = e.^(-1. * delta .* t);
part_2 = (v_0 + delta * x_0) / omega_d .* sin(omega_d .* t);
part_3 = x_0 .* cos(omega_d .* t);
x = part_1 .* (part_2 .+ part_3);
printf("Plotter\n");
printf("Zeitdiagramm\n");
% TODO Normiert auf Schwingungen
plot(t/periodendauer, x);
title("Zeitdiagramm")
xlabel("Zeit in Perioden");
ylabel("Amplitude");
print("ausgabe/zeitdiagramm_zum_phasenportraet.svg", "-dsvg");
printf("Weiter mit beliebiger Taste\n");
pause
printf("Berechne Trajektorie\n");
printf("Trajektorie im Phasenportraet\n");
% part_1 = e.^(-1. * delta .* t);
% part_2 = (v_0 + delta * x_0) / omega_d .* sin(omega_d .* t);
% part_3 = x_0 .* cos(omega_d .* t);
part_1_ableit = -1. * delta * e.^(-1. * delta .* t);
part_2_ableit = (v_0 + delta * x_0) / omega_d .* omega_d .* cos(omega_d .* t);
part_3_ableit = x_0 .* omega_d .* -1. .* sin(omega_d .* t);
% x_ableit = v = Geschwindigkeit
x_ableit = part_1_ableit .* (part_2 .+ part_3) + part_1 .* (part_2_ableit .+ part_3_ableit);
v = x_ableit;
plot(x, v);
%[exit phasenportraet]=system(["echo '" "Phasenporträt mit Trajektorie" "'| iconv --from-code=UTF-8 --to-code=ISO-8859-15"], 1);
phasenportraet = "Phasenportraet mit Trajektorie";
title(phasenportraet)
xlabel("Amplitude (Auslenkung)");
ylabel("Ableitung (Geschwindigkeit)");
text(x_0 - 0.1, v_0, "Startpunkt ->", "HorizontalAlignment","right", "color", "blue");
text(0.2, -0.75, "<- Zeit", "color", "blue");
print("ausgabe/phasenportraet_mit_trajektorie.svg", "-dsvg");
printf("Weiter mit beliebiger Taste\n");
pause
printf("Phasenportraet als Vektorfeld\n");
[x_mesh, v_mesh] = meshgrid (-1.6:.2:+1.6);
% Vektor: v(x, v)
% v_x = v (bekannt)
% v_y = a (zu ermitteln)
a_mesh = -1. .* (2 .* delta .* v_mesh .+ omega_0 .* x_mesh);
% Randbedingungen: x_v und v_
quiver(x_mesh, v_mesh, v_mesh, a_mesh);
hold on;
plot(x, v, "r");
hold off;
% [exit phasenportraet]=system(["echo '" "Phasenporträt mit Vektorfeld" "'| iconv --from-code=UTF-8 --to-code=ISO-8859-15"], 1);
phasenportraet = "Phasenportraet mit Vektorfeld";
title(phasenportraet)
xlabel("Amplitude");
ylabel("Ableitung");
axis ([-1.5, +1.5, -1.5, +1.5 ]);
print("ausgabe/phasenportraet_mit_vektorfeld.svg", "-dsvg");
printf("Weiter mit beliebiger Taste\n");
pause
% EOF