% matematicke kyvadlo s libovolnou pocatecni vychylkou
clear;
clc;
# vymazeme nastaveni vykreslovani
__gnuplot_raw__ "reset;";
# nastavuje opet vykresleni plnou carou
__gnuplot_set__ data style lines;
# obnoveni rozsahu os
__gnuplot_set__ xrange [:]
__gnuplot_set__ yrange [:]
printf("Kmity sekundoveho matematickeho kyvadla s libovolnou vychylkou\n")
% definice poèáteèních podmínek (v základních jednotkách SI)
phi0=input("Zadejte ruzne pocatecni vychylky ve stupnich jako matici [a b c ...]\n")*pi/180;
clearplot
% nastaveni obrazku
__gnuplot_set__ xzeroaxis lt -1
__gnuplot_set__ grid
__gnuplot_set__ nokey
__gnuplot_set__ mxtics
__gnuplot_set__ mytics
__gnuplot_set__ tics out
__gnuplot_set__ xlabel "t(s)"
__gnuplot_set__ ylabel "phi(rad)"
__gnuplot_set__ title "Matematicke kyvadlo s libovolnou vychylkou"
% -----------------------------------------------
% vlastni cyklus
cas_spusteni = clock();
x=v=[];
hold on
for j=1:columns(phi0)
clear t;
clear ph;
t(1)=0;
ph(1)=phi0(j);
dph(1)=0;         % poèateèní derivace (úhlová rychlost)
mez=4;         % mez integrace - 4 male kmity 
h = mez/500;   % krok integrace
g=9.81;
l=g/(4*pi^2);  % delka sekundoveho kyvadla
% -----------------------------------------------
% vlastní numerická integrace - RK metoda 2. øádu
% -----------------------------------------------
i=2;
while t<mez
k1=-g/l*sin(ph(i-1));
php=ph(i-1)+dph(i-1)*2*h/3+k1*2*h^2/9;
k2=-g/l*sin(php);
ph(i)=ph(i-1)+dph(i-1)*h+(k1+k2)*h^2/4;
dph(i)=dph(i-1)+(k1+3*k2)*h/4;
t(i)=t(i-1)+h;
i=i+1;
endwhile
% ukladame vychylku a rychlost do matice pro pozdejsi pouziti...
x=[x;ph];
v=[v;dph];
% -----------------------------------------------
plot(t,ph);
endfor
hold off
celkovy_cas = etime(clock(), cas_spusteni);
disp(['Program bezel ',num2str(celkovy_cas),' sekund.']);
printf("Vykreslime zavislost uhlove frekvence na case\n")
printf("Stisknete libovolnou klavesu...\n")
pause
clearplot
__gnuplot_set__ xlabel "t(s)"
__gnuplot_set__ ylabel "omega(rad/s)"
for j=1:columns(phi0)
hold on
plot(t,v(j,:))
hold off
endfor
printf("Vykreslime zavislost uhlove frekvence na vychylce\n")
printf("Stisknete libovolnou klavesu...\n")
pause
clearplot
__gnuplot_set__ xlabel "phi(rad)"
__gnuplot_set__ ylabel "omega(rad/s)"
for j=1:columns(phi0)
hold on
plot(x(j,:),v(j,:))
hold off
endfor
