# Version vom 17.5.23
# Spuren der Geodäten in Schwarzschild in der (r, \phi)-Ebene


import matplotlib.pyplot as plt
#https://matplotlib.org/

import numpy as np
#https://numpy.org/

from scipy.integrate import odeint
#https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html


### Differentialgleichung zeitartige für Geodätische in Schwarzschild
### u=1/r als Funktion von phi (also die Spur von der Geodäten (ohne die Auflösung in Zeitrichtung und theta) )
### u''=3/2*u^2-u +1/(2*a) 
### zeitartig: a>0 (um so größer a um so kleiner der Drehimpuls (bei festem Schwarzschildradius)  
### lichtartig: a=0
 
def schwarzgeotime(z, phi, a):
    u, y = z
    dzdt = [y, 3./2*u**2-u+a]
    return dzdt


#Anfangswerte ( Start an einem Punkt der Bahn mit r'(\phi)=0)
z0 = [0.35, 0]
#z1 = [0.35, 0]


phi = np.linspace(0, 54,10000)

sol = odeint(schwarzgeotime, z0, phi, args=(0.125,))
#sol2 = odeint(schwarzgeotime, z1, phi, args=(0.125,))


print(phi[1], sol[1])

#Plotten der Lösung r=r(\phi) - in Polarkoordinaten
fig, ax = plt.subplots(subplot_kw={'projection': 'polar'})


u1, y1= zip(*sol)
ax.plot(phi, 1./np.array(u1), color='blue')

#u2, y2= zip(*sol2)
#ax.plot(phi, 1./np.array(u2), color='blue')


plt.show()