import numpy as np
import matplotlib.pyplot as plt
import math

#1. uloha
#nasledujuci postup plati pre lubovolne Q > 1
M = 333000     #pomer M_S/M_Z
r = 149597900
Q = np.linspace(0, 1, 1000)
f = Q**3 + Q**2/(1-Q)**2/M - 1
#f je pracovna funkcia (Q^3+Q^2/(1-Q)^2/M-1 = 0) odvodena na cviceni, kde Q = a/r; a = |SL1|; r = |SZ|
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()

a = 0.0      #spodna hranica
b = 1.0      #horna hranica
Q = a/2 + b/2
while (b - a)/Q > 10e-10:   #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
    Q = a/2 + b/2;
d = (1 - Q)*r     #vzdialenost bodu L1 od Zeme v km
print('vzdialenost bodu L1 od Zeme je', d, 'km')



#2. uloha
a = 17.8    #velka polos
e = 0.967    #excentricita
P = a**(1.5)    #perioda
n = np.arange(0, P)    #casovy vektor v rokoch
M = 2*np.pi/P*n   #stredna anomalia v radianoch
E = M*1
v = M*1
r = M*1
for i in range(len(n)):    #riesenie Keplerovej rovnice pre excentricku anomaliu E
    E[i] = M[i]
    for c in range(10):
        E[i] = M[i] + e*np.sin(E[i])
    v[i] = 2*math.atan(np.sqrt((1 + e)/(1 - e))*np.tan(E[i]/2))    #vypocet pravej anomalie
    if v[i] < 0:    #korekcia pre uhly vacsie ako pi
        v[i] = v[i] + 2*np.pi
    r[i] = a*(1 - e**2)/(1 + e*np.cos(v[i]))    #r = vzdialenost od ohniska
plt.figure()
plt.plot(n,M)
plt.plot(n,E)
plt.plot(n,v)
plt.xlabel('n [rok]')
plt.ylabel('anomalia [rad]')
plt.legend(('stredna anomalia M', 'excentricka anomalia E', 'prava anomalia v'), loc=4)
plt.show()
plt.figure()
plt.polar(v,r,'.')
plt.show()