%1. uloha
%nasledujuci postup plati pre lubovolne Q > 1
M = 333000;     %pomer M_S/M_Z
r = 149597900
Q = 0:0.001:1;
f = Q.^3+Q.^2./(1-Q).^2/M-1;    %pracovna funkcia (Q^3+Q^2/(1-Q)^2/M-1 = 0) ...
    ...odvodena na cviceni, kde Q = a/r; a = |SL1|; r = |SZ|
plot(Q,f)    %vykresli priebeh funkcie v zavislosti na parametri Q
hold on
plot([0,1],[0,0],'k')    %nazorny graf funkcie y = 0
xlabel('Q')
ylabel('f')

a = 0;      %spodna hranica
b = 1;      %horna hranica
Q = a/2+b/2;
while (b-a)/Q > 10^-6   %podmienka ukoncujuca cyklus pri dosiahnuti potrebnej presnosti
    f = Q^3+Q^2/(1-Q)^2/M-1;
    if f > 0
        b = Q;    %meni sa horna hranica
    else
        a = Q;    %meni sa spodna hranica
    end
    Q = a/2+b/2;
end
d = (1-Q)*r     %vzdialenost bodu L1 od Zeme v km


%2. uloha
a = 17.8;    %velka polos
e = 0.967;    %excentricita
P = a^(3/2);    %perioda
n = 0:P;    %casovy vektor v rokoch
M = 2*pi*n/P;    %stredna anomalia v radianoch
for i = 1:size(n,2)    %riesenie Keplerovej rovnice pre excentricku anomaliu E
    E(i) = M(i);
    for c = 1:10
        E(i) = M(i)+e*sin(E(i));
    end
end
v = 2*atan(sqrt((1+e)/(1-e))*tan(E/2));    %vypocet pravej anomalie
for i = 1:size(n,2)    %korekcia pre uhly vacsie ako pi; cyklus je ekvivalentny ...
    ... prikazu v = 2*pi.*(v<0)+v
    if v(i) < 0
        v(i) = v(i)+2*pi;
    end
end
[n;M;E;v]';    %pripadne vypisanie tabulky jednotlivych excentricit
figure
plot(n,M,'r',n,E,'b',n,v,'g')
xlabel('n [rok]')
ylabel('anomalia [rad]')
legend('stredna anomalia M', 'excentricka anomali E', 'prava anomalia v',...
    'Location', 'SouthEast')
r = a*(1-e^2)./(1+e*cos(v));    %r = vzdialenost od ohniska
figure
polar(v,r,'.')
