% pohyb rakety s promennou hmotnosti
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("Pohyb rakety s promennou hmotnosti\n");
% definice poèáteèních podmínek (v zakladních jednotkách SI)
t0=0;
r0=6.378e6; 	% pocatecni vzdalenost od stredu Zeme r=Rz
v0=0;	        % pocatecni rychlost
g=9.81;
q=0.12e3;       % hmotnost paliva spalena za sekundu
m0=13e3;         % pocatecni hmotnost rakety
mk=4.25e3;      % konecna hmotnost rakety
c=2.2e3;        % rychlost spalenych plynu
% -----------------------------------------------
dt=0.05;	% integracni krok
#mez=30;
% -----------------------------------------------
% vlastní numerická integrace - klasická Eulerova metoda
% -----------------------------------------------
cas_spusteni = clock();
t(1)=t0;
x(1)=r0;
v(1)=v0;
a(1)=g;
m(1)=m0;
i=1;
% -----------------------------------------------
while (x(i)-r0>=0) 
   i=i+1;
   % po spaleni paliva je zmena hmotnosti rovna 0
   if (m(i-1)<=mk)
   q=0;
   endif
   m(i)=m(i-1)-q*dt;
   a(i)=(c/m(i-1))*q-g*x(i-1)/r0;
   v(i)=v(i-1)+a(i)*dt;
   x(i)=x(i-1)+v(i)*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__ xlabel "t(s)"
__gnuplot_set__ ylabel "r(m)"
__gnuplot_set__ title "Pohyb rakety s promennou hmotnosti"
plot(t,x-r0);
printf("Vykreslime zavislost hmotnosti na case\n");
printf("Stisknete libovolnou klavesu...\n");
pause
__gnuplot_set__ xlabel "t(s)"
__gnuplot_set__ ylabel "m(t)"
plot(t,m/1e3);
printf("Vykreslime zavislost rychlosti na case\n");
printf("Stisknete libovolnou klavesu...\n");
pause
__gnuplot_set__ xlabel "t(s)"
__gnuplot_set__ ylabel "v(m/s)"
plot(t,v);
printf("Vykreslime zavislost gravitacniho zrychleni na case\n");
printf("Stisknete libovolnou klavesu...\n");
pause
__gnuplot_set__ xlabel "t(s)"
__gnuplot_set__ ylabel "a_g(m/s^2)"
plot(t,g./x.^2*r0^2);
