# Metoda tecen pro reseni Keplerovy rovnice
# Vymazeme hodnoty promennych z predchazejicich vypoctu a obsah terminalu
clear;
clc;
# Definice funkce
function y=f(E)
   y=E-0.99513*sin(E)-0.00413;
endfunction
# Pozadovana presnost
presnost=1e-6;
# Zaciname merit dobu vypoctu	
t0=clock();
# Krajní body
a=0; 		# levy okraj
b=pi; 		# pravy okraj
fa=f(a); 	# leva funkcni hodnota
fb=f(b); 	# prava funkcni hodnota
k=1; 		# krok vypoctu 
x(k)=(a*fb-b*fa)/(fb-fa); # aproximace øe¹ení
fk=f(x(k));       # hodnota funkce ve stredu intervalu
# Formátovaný výstup hodnot
printf("k=%f x(k)=%f f(xk)=%.20f\n", [k,x(k),fk]);
clear -all
t0=clock();
   # definice funkce
function y=f(E)
 global a=187.8*1.4959787e11;
 global rp=0.9141*1.4959787e-11;
 kapa=6.672e-11;
 M=1.9891e30;
 t=618*24*3600;
   y=E-((a-rp)/a)*sin(E)-sqrt(kapa*M/(a^3))*t;
endfunction
   # první derivace funkce
function y = g(E)
   global a rp
   y=1-((a-rp)/a)*cos(E);
endfunction;
   # druhá derivace funkce
function y = h(E)
global a rp
   y=((a-rp)/a)*sin(E);
endfunction;
   # Zadání požadované pøesnosti
presnost=1e-6;
   # krajní body
a=0;
b=pi;
fa=f(a);
ha=h(a);
k=1;
if (fa*ha>0)
x1=a;
else
x1=b;
endif
f1=f(x1);
g1=g(x1);
xk=x1-f1/g1;
fk=f(xk);
printf("k=%f xk=%f f(xk)=%f\n", [k,xk,fk]);
#cyklus
while (abs(xk-x1)>presnost)
k=k+1;
x1=xk;
f1=f(x1);
g1=g(x1);
xk=x1-f1/g1;
fk=f(xk);
printf("k=%f xk=%f f(xk)=%f\n", [k,xk,fk]);
endwhile
casvypoctu=etime(clock(),t0)
