%Kapitza vagy inverz pendulum

%alap paraméterek:
 Xinit= -3; %kezdeti értékek, ezektõl függés
 Yinit=0.5;
 tinit=0; %kezdeti idõ
 Max_Time= 500; %eddig megy a szimuláció
 %instabil egyensúlyi pont körülre áll be a rendszer
 wk=10;%1.5; % a függõleges oszcilláció frekvenciája a felfüggeszésnél
 a=20;%6;%0.5; % a függõleges oszcilláció amplitúdója a felfüggeszésnél
 b=0.1;
 
 x=-20:0.05:20;
 y=-4:0.05:4;
 [X,Y] = meshgrid(x,y); %hálót készít x,y vektorok alapján (az energiák kiszámításához az adott pontokban)
 energia_inga=@(X,Y) 0.5.*Y.^2-cos(X)+1;
 Z = energia_inga(X,Y); %energiát kiszámol, Z-ben tárol
 
 %egyenlet:
 egyenlet=@(t, xy) [xy(2); -b*xy(2)-sin(xy(1))+a*cos(wk*t)*sin(xy(1))];
 [t, xy]=ode45(egyenlet, [tinit, Max_Time], [Xinit; Yinit]); %diffegyenletet kiszámol
 
 figure;
subplot(2,2,1); %ábrák elhelyezése egy mezõn belül, számolás balról jobbra, fentrõl le
hold on; %több ábra egymásra
energiaszintek=[0, 0.1 0.25, 0.5, 1, 1.4, 2, 2.4, 3, 4]; %ezeken a szinteken számolom ki az energiát
contour(X,Y,Z, energiaszintek); %energia szintvonalak ábrázolása
plot(xy(:,1), xy(:,2), 'b'); %x-y ábrázolás
xlabel('x');
ylabel('y');
hold off; %ez után már felülír ábrát

subplot(2,2,3);
plot(xy(:,1), t, 'b'); %x-t ábrázolás
xlabel('x');
ylabel('t');

subplot(2,2,2);
plot(t, xy(:,2), 'b'); %t-y ábrázolás
xlabel('t');
ylabel('y');

subplot(2,2,4);
plot3(xy(:,1), xy(:,2), t, 'b'); %x-y-t 3D-s ábrázolás
xlabel('x');
ylabel('y');
zlabel('t');