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("Balisticke krivky pro ruzne elevacni uhly");
% definice poèáteèních podmínek (v zakladních jednotkách SI)
t0=0;
x0=0;
y0=0;         
printf("Pocatecni rychlost:\n");
v0=200         % poèáteèní rychlost
printf("Elevacni uhlu pro ktere bude vykreslena zavislost:\n");
eluhel=[10 15 20 30 37.9 45 60 75]     % elevaèní úhly ve stupních
alpha=eluhel*pi/180;
% -----------------------------------------------
% krok integrace
dt = 0.05;
% parametry prostredi
R=0.04;        	% polomìr støely
m=2.5;  	% hmotnost støely
ro=1.3;        	% hustota vzduchu 
C=0.48;        	% koeficient odporu
g=9.81; 	% tíhové zrychlení
S=pi*R^2;
k=C*ro*S/2;
cas_spusteni = clock();
% cyklus pro hodnoty elevacnich uhlu
clearplot;
hold on % pro vice obrazku do jednoho grafu
__gnuplot_set__ zeroaxis lt -1
__gnuplot_set__ yrange [-5:]
__gnuplot_set__ nokey
__gnuplot_set__ grid
__gnuplot_set__ tics out
__gnuplot_set__ mxtics
__gnuplot_set__ mytics
__gnuplot_set__ xlabel "x(m)"
__gnuplot_set__ ylabel "y(m)"
__gnuplot_set__ title "Balisticke krivky pro ruzne elevacni uhly"
% -----------------------------------------------
for j=1:columns(alpha)
clear t
clear x
clear y
vx=v0*cos(alpha(j));
vy=v0*sin(alpha(j));
v=sqrt(vx^2+vy^2);
% -----------------------------------------------
% vlastní numerická integrace - klasická Eulerova metoda
% -----------------------------------------------
t(1)=t0;
x(1)=x0;
y(1)=y0;
i=1;
 while (y(i)>=0)  % testujeme dopad na rovinu y=0)
   i=i+1;  
   x(i)=x(i-1)+vx*dt;
   y(i)=y(i-1)+vy*dt;
   ax=-k*v*vx/m;
   ay=-g-k*v*vy/m;
   vx=vx+ax*dt;
   vy=vy+ay*dt;
   v=sqrt(vx^2+vy^2);
   t(i)=t(i-1)+dt;
 endwhile  
% -----------------------------------------------
% vykreslení zavislosti
plot(x,y); 
endfor
hold off
celkovy_cas = etime(clock(), cas_spusteni);
disp(['Program bezel ',num2str(celkovy_cas),' sekund.']);
