% numerické øešení pohybu kosmické sondy
% 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("Pohyb sondy v soustave Zeme-Mesic\n")
RZ=6.378e6;	% polomìr Zemì
t(1)=0;
xsZ(1)=0;
ysZ(1)=10*RZ;
vx=2.7e3;	% x-ová složka startovní rychlosti
vy=1.84e3;	% y-ová složka startovní rychlosti
xM=60.13*RZ;
yM=0;
mez= 1.2e6;     % mez integrace
h = mez/6000;	% krok integrace
% parametry prostøedí, konstanty
k=6.67259e-11;	% gravitaèní konstanta
MZ=5.98e24;	% hmotnost Zemì
MM=7.374e22;	% hmotnost Mìsíce
% -----------------------------------------------
% vlastní numerická integrace - RG metoda 2. øádu
% -----------------------------------------------
t0=clock();
i=2;
while (t<mez & xsZ(i-1)^2+ysZ(i-1)^2>RZ^2)
k1x=-k*xsZ(i-1)*MZ/(xsZ(i-1)^2+ysZ(i-1)^2)^1.5-k*MM*(xsZ(i-1)-xM)/((xsZ(i-1)-xM)^2+ysZ(i-1)^2)^1.5;
k1y=-k*ysZ(i-1)*MZ/(xsZ(i-1)^2+ysZ(i-1)^2)^1.5-k*MM*ysZ(i-1)/((xsZ(i-1)-xM)^2+ysZ(i-1)^2)^1.5;
xp=xsZ(i-1)+vx*2*h/3+k1x*2*h^2/9;
yp=ysZ(i-1)+vy*2*h/3+k1y*2*h^2/9;
k2x=-k*xp*MZ/(xp^2+yp^2)^1.5-k*MM*(xp-xM)/((xp-xM)^2+yp^2)^1.5;
k2y=-k*yp*MZ/(xp^2+yp^2)^1.5-k*MM*yp/((xp-xM)^2+yp^2)^1.5;
xsZ(i)=xsZ(i-1)+vx*h+(k1x+k2x)*h^2/4;
ysZ(i)=ysZ(i-1)+vy*h+(k1y+k2y)*h^2/4;
vx=vx+(k1x+3*k2x)*h/4;
vy=vy+(k1y+3*k2y)*h/4;
t(i)=t(i-1)+h;
i=i+1;
endwhile  
% -----------------------------------------------
celkovy_cas = etime(clock(), t0);
disp(['Program bezel ',num2str(celkovy_cas),' sekund.']);
% nastaveni obrazku
clearplot
__gnuplot_set__ nokey
__gnuplot_set__ grid
__gnuplot_set__ mxtics
__gnuplot_set__ mytics
__gnuplot_set__ tics out
__gnuplot_set__ xlabel "x(m)"
__gnuplot_set__ ylabel "y(m)"
__gnuplot_set__ title "Pohyb druzice v soustave Zeme-Mesic"
axis('equal');
% vykreslení zavislosti
hold on
% Vykresleni Zeme
u=linspace(0,2*pi,100);
plot(cos(u),sin(u),'b')
% Vykresleni Mesice
plot(1738000/RZ*cos(u)+xM/RZ,1738000/RZ*sin(u),'g')
% vykresleni trajektorie
plot(xsZ/RZ,ysZ/RZ,'r');
hold off
