import sympy as sp
phi = (1+sp.sqrt(5))/2
# golden ratio identities
print("phi =", sp.N(phi), " phi^2 =", sp.simplify(phi**2), "= phi+1 ?", sp.simplify(phi**2-phi-1))
print("r_isco = 3+sqrt5 =", sp.simplify(3+sp.sqrt(5)), " ; 2 phi^2 =", sp.simplify(2*phi**2), " equal?", sp.simplify(3+sp.sqrt(5)-2*phi**2))
# psi_isco = (3-sqrt5)/2
psi_isco = (3-sp.sqrt(5))/2
print("psi_isco = (3-sqrt5)/2 =", sp.N(psi_isco))
print("  2 - phi =", sp.simplify(2-phi), "=", sp.N(2-phi))
print("  phi^{-2} =", sp.simplify(phi**-2), "=", sp.N(phi**-2))
print("  (3-sqrt5)/2 == 2-phi+ (1/2)? check 2-phi=",sp.N(2-phi))
# 3-sqrt5 = 2(2-phi)?  2-phi = 2 - (1+sqrt5)/2 = (3-sqrt5)/2. So psi_isco = 2-phi = phi^{-2}? 
print("  is psi_isco == 2-phi?", sp.simplify(psi_isco-(2-phi)))
print("  is psi_isco == phi^-2?", sp.simplify(psi_isco-phi**-2))
# redshift z(r) = e^{u} - 1 for static emitter; here surface redshift uses sqrt(-gtt)=e^{-psi/2}
# gravitational redshift between emitter at r and observer at infinity:
# 1+z = sqrt(-gtt(inf))/sqrt(-gtt(r)) = 1 / e^{-psi/2} = e^{psi/2} = e^{u}   (psi=2u)
r,u = sp.symbols('r u', positive=True)
print("\n1+z = e^{psi/2} = e^{u}, u=GM/c^2 r. z = e^u - 1.  (so z=e^u-1)")
# z at photon sphere r=2 -> u=1/2 ; at ISCO r=3+sqrt5 -> u=1/(3+sqrt5)=(3-sqrt5)/4
zph = sp.exp(sp.Rational(1,2))-1
print("z(photon sphere, u=1/2) = e^{1/2}-1 =", sp.N(zph))
u_isco = sp.simplify(1/(3+sp.sqrt(5)))
print("u_isco =", sp.simplify(u_isco), "=", sp.N(u_isco))
z_isco = sp.exp(u_isco)-1
print("z(ISCO) = e^{(3-sqrt5)/4}-1 =", sp.N(z_isco), " exponent (3-sqrt5)/4=", sp.N(u_isco))
