% micek odrazejici se od stolu
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 mechanickeho oscilatoru\n");
% definice poèáteèních podmínek (v zakladních jednotkách SI)
t0=0;
% pocatecni vychylka
y0=input("Zadejte pocatecni vychylku y0: ");
% pocatecni rychlost
v0=input("Zadejte pocatecni rychlost v0: ");
% soucinitel tlumeni
b=input("Zadejte soucinitel tlumeni b: ");
k=1;            % tuhost 
m=1;            % hmotnost
dt=0.01;	% integracni krok
mez=8*pi*sqrt(m/k);       % mez casu - 4 kmity
% -----------------------------------------------
% vlastní numerická integrace - klasická Eulerova metoda
% -----------------------------------------------
cas_spusteni = clock();
t(1)=t0;
y(1)=y0;
v(1)=v0;
i=1;
 while (t(i)<=mez) 
   i=i+1;
   a=-k/m*y(i-1)-b/m*v(i-1);
   v(i)=v(i-1)+a*dt;
   y(i)=y(i-1)+v(i-1)*dt;
   t(i)=t(i-1)+dt;
 endwhile  
% -----------------------------------------------
celkovy_cas = etime(clock(), cas_spusteni);
disp(['Program bezel ',num2str(celkovy_cas),' sekund.']);
% -----------------------------------------------
% vykreslení zavislosti
% nastaveni obrazku
__gnuplot_set__ grid
__gnuplot_set__ xzeroaxis lt -1
__gnuplot_set__ nokey
__gnuplot_set__ tics out
__gnuplot_set__ mxtics
__gnuplot_set__ mytics
__gnuplot_set__ xlabel "t(s)"
__gnuplot_set__ ylabel "y(cm)"
__gnuplot_set__ title "Kmity oscilatoru"
plot(t,y);
printf("Vykreslime zavislost rychlosti na case\n");
printf("Stisknete libovolnou klavesu...\n");
pause
__gnuplot_set__ xlabel "t(s)"
__gnuplot_set__ ylabel "v(cm/s)"
plot(t,v);
