import numpy as np
import matplotlib.pyplot as plt

# Odvodenie riesenia ukazuje ukol2_1.jpg.
M = 1/333000     # pomer M_Z/M_S
r = 149597900
# Pre Q = 1 nie je funkcia definovana, preto upravime hornu hranicu.
Q = np.linspace(0, 0.99999, 1000)
f = M*Q**2/(1 - Q)**2 + (1 + M)*Q**3 - M*Q**2 - 1
plt.plot(Q, f)
plt.ylim([-1, 3])
plt.plot([0, 1], [0, 0], 'k')    #nazorny graf funkcie y = 0
plt.xlabel('Q')
plt.ylabel('f')
plt.show()
# Vidime, ze funkcia je na intervale [0, 1] monotonne rastuca a pre urcite Q
# nadobuda nulovu hodnotu.

# Funkcia f je polynomicka piateho radu a analyticky takmer neriesitelna.
# Preto pouzijeme metodou delenia intervalov, kde vzdy porovname funkcne
# hodnoty na hraniciach intervalu a stred intervalu vyberieme ako novu hranicu
# tak aby vzdy platilo: f(a) < 0, f(b) > 0.
a = 0.0      # spodna hranica
b = 1.0      # horna hranica
Q = a/2 + b/2
# Podmienka ukoncujuca cyklus pri dosiahnuti potrebnej presnosti.
while (b - a)/Q > 10e-10:
    f = M*Q**2/(1 - Q)**2 + (1 + M)*Q**3 - M*Q**2 - 1
    if f > 0:
        b = Q    # meni sa horna hranica
    else:
        a = Q    # meni sa spodna hranica
    Q = a/2 + b/2;
d = (1 - Q)*r     # vzdialenost bodu L1 od Zeme v km
print('Vzdialenost bodu L1 od Zeme je {} km.'.format(np.round(d)))
