%%
% Nonlinear Dynamical Systems -- 2023 (P-ITEEA-0037)
% 1. gyakorlati óra
% 
% 2023. szeptember 15.-én élőben írt programkód
% 
% Az inga mozgásegyenlete 

m = 1;    % Az ingára rögzített tömegpont tömege [kg]
g = 9.8;  % Gravitácios gyorsulás [m/s^2]
l = 9.8;  % Az inga hossza [m]
b = 0;    % Közegellenállási együttható [...vezessétek le]
w = 1.2;  % Periodikus gerjesztés körfrekvenciája
a = 1;    % Periodikus gerjesztés amplitúdója

% Az inga mozgásegyenletének jobb oldala (vektormező)
f = @(t,x) [
    x(2)
    -b/m*x(2)-g/l*sin(x(1)) + a*cos(w*t)
    ];

% A vektormező kirajzolásához szükséges rácspontok
lx1 = linspace(-7*pi, 3*pi, 50);
lx2 = linspace(-3,3,40);
[x1, x2] = meshgrid(lx1, lx2);

% Az inga összenergiája
E = m*g*l*(1 - cos(x1)) + m*l^2*x2.^2 / 2;

% A gerjesztés nélküli mozgásegyenlet vektormezőjének kirajzolásához kell. 
% A vektormező két koordinátáját külön-külön kiszámolom.
f1 = x2;
f2 = -b/m*x2-g/l*sin(x1);

% A vektormező vektorait normalizálom, hogy a megjelenítéskor minden vektor azonos
% hosszúságú legyen.
r = sqrt( f1.^2 + f2.^2 );
f1 = f1 ./ r;
f2 = f2 ./ r;

% Megoldom a mozgásegyenletet egy numerikus ODE szolverrel (megoldóval)
x0 = [ 0 ; 1 ];
[t,x] = ode45(f,[0 100],x0);

% Vizualizálás következik

fig = figure(1);
hold off
quiver(x1,x2,f1,f2),
hold on
contour(x1,x2,E)

Pl = plot(x(:,1),x(:,2),'LineWidth',3);
plot(x(1,1),x(1,2),'.','Color',Pl.Color,'MarkerSize',50)

fig.Position(3:4) = [1719 436];
