% numericke reseni rychlosti lyzare na svahu s odporem prostredi
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("Rychlost lyzare na svahu s odporem prostredi\n")
% definice po‚è‚áte‚èn‚ích podminek
t(1)=0;
s(1)=0;     
v(1)=0;         % po‚è‚áte‚èn‚í rychlost
m=90;           % hmotnost lyža‚øe
uhel=15;        % ‚úhel svahu ve stupn‚ích
alpha=uhel*pi/180;
f=0.055;        % sou‚èinitel smykov‚ého t‚øen‚í
vm=18;          % mezni rychlost
g=9.81;  	% t‚íhov‚é zrychlen‚í
K=m*g*(sin(alpha)-f*cos(alpha))/vm^2;
dt = 1;         % krok integrace
% vlastni cyklus
cas_spusteni = clock();
% -----------------------------------------------
% vlastn‚í numerick‚á integrace - metoda Feynmannova
% -----------------------------------------------
   i=1;
   a=g*(sin(alpha)-f*cos(alpha))-K*v^2/m;
 while (v<=0.9999*vm)  % 
   v5=v(i)+a*dt/2;
   i=i+1;
   s(i)=s(i-1)+v5*dt;
   v(i)=v5+a*dt/2;
   t(i)=t(i-1)+dt;
   a=g*(sin(alpha)-f*cos(alpha))-K*v(i)^2/m;
 endwhile  
% -----------------------------------------------
celkovy_cas = etime(clock(), cas_spusteni);
disp(['Program bezel ',num2str(celkovy_cas),' sekund.']);
# % vykreslen‚í z‚ávislosti v na t
# nastaveni obrazku
__gnuplot_set__ nokey
__gnuplot_set__ grid
__gnuplot_set__ tics out
__gnuplot_set__ mxtics
__gnuplot_set__ mytics
__gnuplot_set__ tics out
__gnuplot_set__ xlabel "t(s)"
__gnuplot_set__ ylabel "v(m/s)"
__gnuplot_set__ title "Zavislost rychlosti lyzare na case"
plot(t,v);
