# Version vom 17.5.23
# 
# Ableiten der Schwarzschild-Metrik aus den Bedingungen:
# Eine statische (= alle Metrikkoeffizienten sind zeitunabhaengig) Raumzeit
# (Signatur = (-,+,+,+)) auf dem R^4, die invariant unter raeumlichen 
# Rotationen ist und die Einsteinvakuumgleichungen (Ric=0) erfuellt.

from sympy import symbols, diag, sin, Function, simplify, dsolve
#https://de.wikipedia.org/wiki/SymPy
#https://www.sympy.org/en/index.html 

from einsteinpy.symbolic import MetricTensor, RicciTensor
#https://einsteinpy.org/

syms = symbols('t r theta phi')
t, r, theta, phi = syms

#Allgemeine Form einer statischen (räumlich-)rotationsinvarianten Lorentz-Metrik auf R^4  
g = MetricTensor(diag(-Function('A')(r), Function('B')(r), r**2, r**2*sin(theta)**2).tolist(), syms)

#Ziel ist es alle A und B oben so zu finden, dass Ric=0 ist.
ricci = RicciTensor.from_metric(g)          # Ricci-Tensor als (0,2)-Tensor
uric = ricci.change_config(metric=g)      # Ricci-Tensor als (1,1)-Tensor

print('Ricci (1,1)-Tensor:', uric.tensor())

# Ist uric Null, dann auch ricci. Wir arbeiten mit uric.
# Die letzten beiden Eintraege sind gleich:
print('\n Ric^3_3 - Ric^2_2=', simplify(uric[3,3]-uric[2,2]))
 
# Damit zumindest die  Differenz der ersten beiden Einträge Null ist, berechnen wir: 
print('\n (rB(r))^2*A(r)*(Ric^0_0-Ric^1_1)=', simplify(r*Function('B')(r)**2*Function('A')(r)*(uric[0,0]-uric[1,1])))

# Dies ist -(AB)'
# also muss AB konstant sein - hier gleich klein a
a = symbols('a')

# Einsetzen in uric[2,2]
ode = simplify(r*uric[2,2].subs({Function('B')(r): a/Function('A')(r)}))
print('\n Für B(r)=a/A(r) ist r*Ric^2_2=', ode)

# Loesen (und Ausgeben) von ode=0
print('\n Dies ist gleich Null falls ', dsolve(ode, Function('A')(r)))




c = symbols('c')
print('\nProbe:') 
g2 = MetricTensor(diag(-c/r-a , a/(c/r+a), r**2, r**2*sin(theta)**2).tolist(), syms)
ricci2 = RicciTensor.from_metric(g2)
uric2 = ricci2.change_config(metric=g2)   
print('\n Ric^0_0=', simplify(uric2[0,0]))
print('\n Ric^1_1=', simplify(uric2[1,1]))

print('\n Ric^2_2=', simplify(uric2[2,2]))
print('\n Ric^3_3=', simplify(uric2[3,3]))   # R