# Metoda puleni intervalu pro reseni Keplerovy rovnice
# Vymazeme hodnoty promennych z predchazejicich vypoctu a vypis na terminalu
clear;
clc;
# zadani konstant
a=187.8*1.4959787e11;
rp=0.9141*1.4959787e11;
kappa=6.672e-11;
M=1.9891e30;
t=618*24*3600;
# Definice funkce
function y=f(E)
   y=E-0.99513*sin(E)-0.00413;
endfunction
# Pozadovana presnost
presnost=1e-10;
# Zaciname merit dobu vypoctu	
t0=clock();
# Krajni 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);       # hodnota funkce ve stredu intervalu
# Formatovany vypis hodnot
printf("k=%f xk=%f f(xk)=%.20f\n", [k,xk,fk]);
# Cyklus, ktery se opakuje do dosazeni pozadovane presnosti
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);
# Formatovany vypis hodnot
printf("k=%f xk=%f f(xk)=%.20f\n", [k,xk,fk]);
endwhile
# Ukoncime mereni trvani vypoctu a vypiseme cas potrebny na prubeh cyklu 
casvypoctu=etime(clock(),t0)
# Vypocteme polohu a rychlost v pozadovanem okamziku
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
