% numericke reseni balistickeho problemu
%
clear;
clc;
# vymazeme nastaveni vykreslovani
__gnuplot_raw__ "reset;";
# nastavuje opet vykresleni plnou carou
__gnuplot_set__ data style lines;
__gnuplot_set__ xrange [:]
__gnuplot_set__ yrange [:]
% -----------------------------------------------
printf("Balisticka krivka pomoci preddefinovanych procedur\n");
% -----------------------------------------------
% definice pocatecnich podminek
t0=0;
x0=0;
y0=0;         
%v=200;       % pocatecni rychlost
%eluhel=32;   % elevacni uhel ve stupnich
v=input('zadej pocatecni rychlost v m/s: ');
alpha=input('elevacni uhel ve stupnich: ')*pi/180;
vx=v*cos(alpha);
vy=v*sin(alpha);
% -----------------------------------------------
% parametry prostredi a strely
R=0.04;        	        % polomer strely
global m=2.5; 		% hmotnost strely
ro=1.3;        	        % hustota vzduchu 
% C=0.48;        	% koeficient odporu
C=input('zadej koeficient odporu C (<0,1>): ');
global g=9.81;		% tihove zrychleni
S=pi*R^2;
global k=C*ro*S/2;
% -----------------------------------------------
% vlastni numericka integrace - preddefinovana procedura programu octave
% -----------------------------------------------
% resi jako ctyrrozmerny vektor xp o slozkach
% xp(1)=x, xp(2)=vx, xp(3)=y, xp(4)=vy
% piseme diff. rovnice
% xdot(1)=1. der. x
% xdot(2)=2. der. x
% xdot(3)=1. der. y
% xdot(4)=2. der. y
function [xdot]=fce(xp,t)    
 global k m g % prebira globalne definovane konstanty
  xdot=zeros(4,1);
  xdot(1)=xp(2);
  xdot(2)=-k/m*sqrt(xp(2).^2+xp(4).^2).*xp(2);  
  xdot(3)=xp(4);
  xdot(4)=-g-k/m*sqrt(xp(2).^2+xp(4).^2).*xp(4);
endfunction  
% -----------------------------------------------
% meze integrace
dm=0;
hm=2*vy/g;
% -----------------------------------------------
cas_spusteni = clock();
% cas jako nezavisla promenna
t=linspace(dm,hm,200);
yp=lsode("fce",[x0,vx,y0,vy],t);
x=yp(:,1);
y=yp(:,3).*(yp(:,3)>-1);
% -----------------------------------------------
celkovy_cas = etime(clock(), cas_spusteni);
disp(['Program bezel ',num2str(celkovy_cas),' sekund.']);
% -----------------------------------------------
% vykresleni zavislosti
% vykresleni osy x
__gnuplot_set__ xzeroaxis lt -1
% nastaveni rozsahu na ose y na cisla vetsi nez -5
__gnuplot_set__ yrange [-5:]
__gnuplot_set__ nokey
__gnuplot_set__ tics out
__gnuplot_set__ grid
__gnuplot_set__ mxtics
__gnuplot_set__ mytics
% popisky obrazku
__gnuplot_set__ xlabel "x(m)"
__gnuplot_set__ ylabel "y(m)"
__gnuplot_set__ title "Balisticka krivka"
plot(x,y);
