% 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 s harmonickym buzenim\n");
% definice poèáteèních podmínek (v zakladních jednotkách SI)
t0=0;
% pocatecni vychylka
y0=0;
% pocatecni rychlost
v0=0;
% soucinitel tlumeni
b=0.06;
% amplituda budici sily
Fm=0.3;
% uhlova frekvence budici sily
om=0.7;
k=1;            % tuhost 
m=1;            % hmotnost
dt=0.02;	% integracni krok
mez=60*pi/om;   % mez casu - 30 budicich kmitu
% -----------------------------------------------
% 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)+Fm/m*sin(om*t(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), Fm(N)"
__gnuplot_set__ title "Kmity mechanickeho oscilatoru s harmonickym buzenim"
plot(t,y,t,Fm*sin(om*t));
% -----------------------------------------------
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);
% %%%%%%%%%%%%%%%%%%%%%%%% End of file %%%%%%%%%%%%%%%%%%%%%%%%
