# Metoda puleni intervalu pro reseni Keplerovy rovnice
clear;
clc;
# vymazeme nastaveni vykreslovani
__gnuplot_raw__ "reset;";
# nastavuje opet vykresleni plnou carou
__gnuplot_set__ data style lines;
printf("Numericke reseni Keplerovy rovnice metodou puleni intervalu\n");
# zadani konstant
a=187.8*1.4959787e11
rp=0.9141*1.4959787e11
kappa=6.672e-11
M=1.9891e30
t=618*24*3600
printf("Koeficienty rovnice\n");
[-((a-rp)/a) -sqrt(kappa*M/(a^3))*t]
# Definice funkce
function y=f(E)
   y=E-0.9951326*sin(E)-0.0041307;
#   y=E-((a-rp)/a)*sin(E)-sqrt(kappa*M/(a^3))*t;
endfunction
# Zad‚án‚í požadovan‚é p‚øesnosti
presnost=1e-6;
# Zaciname merit dobu vypoctu	
t0=clock();
# Krajn‚í body a stred intervalu
x1=0; 		# levy okraj
x2=pi; 		# pravy okraj
f1=f(x1); 	# leva funkcni hodnota
f2=f(x2); 	# prava funkcni hodnota
k=1; 		# krok vypoctu 
xk=0.5*(x1+x2);	# stred intervalu
fk=f(xk);
   # Form‚átovan‚ý v‚ýstup hodnot
printf("k=%f xk=%f f(xk)=%.20f\n", [k,xk,fk]);
   # Cyklus
while abs(x1-xk)>2*presnost 
if (f1*fk<0)
  x2=xk;
else
  x1=xk;
endif
k=k+1;
xk=0.5*(x1+x2);
fk=f(xk);
   # Form‚átovan‚ý v‚ýstup hodnot
printf("k=%f xk=%f f(xk)=%.20f\n", [k,xk,fk]);
endwhile
casvypoctu=etime(clock(),t0)
printf("Vykreslime zavislost\n");
printf("Stisknete libovolnou klavesu...\n");
pause
# vykresleni - rozsah x
xp=linspace(0,pi,100);
plot(xp,f(xp));
printf("Vykreslime detailnejsi pohled\n");
printf("Stisknete libovolnou klavesu...\n");
pause
# vykresleni - rozsah x s osou y=0
__gnuplot_set__ xzeroaxis lt -1
__gnuplot_set__ mxtics
__gnuplot_set__ mytics
__gnuplot_set__ tics out
__gnuplot_set__ grid
__gnuplot_set__ xlabel "E"
__gnuplot_set__ ylabel "f(E)"
xp=linspace(0,1,100);
plot(xp,f(xp));
printf("Pouziti vestavene funkce fsolve\n");
printf("Stisknete libovolnou klavesu...\n");
pause
format long
fsolve("f",0)
printf("Dopocitame vzdalenost (v AU) a rychlost (v km/s) v dobe objevu\n");
printf("Stisknete libovolnou klavesu...\n");
pause
format short
r=sqrt((a*cos(xk)-(a-rp))^2+(sqrt(a^2-(a-rp)^2)*sin(xk))^2)/1.5e11
v=sqrt(kappa*M*(2/(r*1.5e11)-1/a))/1e3
