from numpy import *
from matplotlib.pyplot import *
from math import *

#1. uloha
#nasledujuci postup plati pre lubovolne Q > 1
M = 333000     #pomer M_S/M_Z
r = 149597900
Q = linspace(0,1,1000)
f = power(Q, 3) + power(Q, 2)*power(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|
plot(Q,f)
ylim([-1,3])
plot([0,1],[0,0],'k')    #nazorny graf funkcie y = 0
xlabel('Q')
ylabel('f')

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 = power(Q, 3) + power(Q, 2)*power(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 = arange(0,P)    #casovy vektor v rokoch
M = 2*pi/P*n   #stredna anomalia v radianoch
E = M*0
v = M*0
r = M*0
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*sin(E[i])
    v[i] = 2*atan(sqrt((1+e)/(1-e))*tan(E[i]/2))    #vypocet pravej anomalie
    if v[i] < 0:    #korekcia pre uhly vacsie ako pi
        v[i] = v[i]+2*pi
    r[i] = a*(1-e**2)/(1+e*cos(v[i]))    #r = vzdialenost od ohniska
plot(n,M)
plot(n,E)
plot(n,v)
xlabel('n [rok]')
ylabel('anomalia [rad]')
legend(('stredna anomalia M', 'excentricka anomalia E', 'prava anomalia v'), loc=4)
show()
figure(figsize=(20,20))
polar(v,r,'.')