%
% 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