#!/usr/bin/env python
# coding: utf-8

# # Schwarzschild horizon in Eddington-Finkelstein coordinates
# 
# This Jupyter/SageMath notebook is relative to the lectures
# [Geometry and physics of black holes](https://relativite.obspm.fr/blackholes/).
# 
# The involved computations are based on tools developed through the [SageManifolds](https://sagemanifolds.obspm.fr) project.

# *NB:* a version of SageMath at least equal to 9.2 is required to run this notebook: 

# In[1]:


version()


# First we set up the notebook to display mathematical objects using LaTeX formatting:

# In[2]:


get_ipython().run_line_magic('display', 'latex')


# ## Spacetime
# 
# We declare the spacetime manifold $M$:

# In[3]:


M = Manifold(4, 'M', structure='Lorentzian')
print(M)


# and the **Eddington-Finkelstein coordinates** $(t,r,\theta,\phi)$ as a chart on $M$:

# In[4]:


X.<t,r,th,ph> = M.chart(r't r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi')
X


# The mass parameter and the metric tensor:

# In[5]:


var('m', domain='real')
assume(m>=0)


# In[6]:


g = M.metric()
g[0,0] = -(1-2*m/r)
g[0,1] = 2*m/r
g[1,1] = 1+2*m/r
g[2,2] = r^2
g[3,3] = (r*sin(th))^2
g.display()


# In[7]:


g[:]


# In[8]:


g.inverse()[:]


# In[9]:


g.christoffel_symbols_display()


# Let us check that we are dealing with a solution of Einstein's equation in vacuum: 

# In[10]:


g.ricci().display()


# ## The scalar field $u$ defining the horizon

# In[11]:


u = M.scalar_field(coord_expression={X: (1-r/(2*m))*exp((r-t)/(4*m))}, name='u')
u.display()                                     


# In[12]:


plm2 = implicit_plot(u.expr().subs(m=1)+2, (r, 0, 15), (t, -10, 15), color='green') + \
       text('$u=-2$', (11,12), color='green')
plm1 = implicit_plot(u.expr().subs(m=1)+1, (r, 0, 15), (t, -10, 15), color='green') + \
       text('$u=-1$', (6.3,12), color='green')
pl0 = implicit_plot(u.expr().subs(m=1), (r, 0, 15), (t, -10, 15), color='green') + \
       text('$u=0$', (3.1, 12), color='green') 
pl1 = implicit_plot(u.expr().subs(m=1)-1, (r, 0, 15), (t, -10, 15), color='green')  + \
       text('$u=1$', (0.8, 0.5), color='green') 
pl2 = implicit_plot(u.expr().subs(m=1)-2, (r, 0, 15), (t, -10, 15), color='green')  + \
       text('$u=2$', (-0.2, -2.5), color='green', background_color='white') 
graph = plm2+plm1+pl0+pl1+pl2
show(graph, aspect_ratio=True, axes_labels=[r'$r/m$', r'$t/m$'], figsize=8)


# In[13]:


graph.save('def_Schwarz_Hu.pdf', aspect_ratio=True, 
           axes_labels=[r'$r/m$', r'$t/m$'], figsize=8)


# In[14]:


du = diff(u)
print(du)
du.display()


# In[15]:


grad_u = du.up(g)
print(grad_u)
grad_u.display()


# Let us check that each hypersurface $u={\rm const}$ is a null hypersurface:

# In[16]:


g(grad_u, grad_u).expr()


# ## The null normal $\ell$

# In[17]:


rho = - log(-grad_u[[0]])
print(rho)
rho.display()


# In[18]:


l = - exp(rho) * grad_u
l.set_name('l', latex_name=r'\ell')
print(l)
l.display()


# In[19]:


graph_l = l.plot(ambient_coords=(r,t), ranges={r:(0.01,8), t:(0,8)}, 
                 fixed_coords={th:pi/2, ph:pi}, parameters={m:1}, color='green', 
                 scale=0.8, aspect_ratio=1)
show(graph_l, figsize=8)


# Let us check that $\ell$ is a null vector everywhere:

# In[20]:


g(l,l).expr()


# In[21]:


l_form = l.down(g)
l_form.set_name('lf', latex_name=r'\underline{\ell}')
print(l_form)
l_form.display()


# In[22]:


nab = g.connection()
print(nab)
nab


# In[23]:


nab_l_form = nab(l_form)
nab_l_form.display()


# In[24]:


nab_l_form.symmetrize().display()


# Check of the identity $\ell^\mu \nabla_\alpha \ell_\mu=0$:

# In[25]:


v = l.contract(nab_l_form, 0)
v.display()


# ## The null normal as a pregeodesic vector field

# In[26]:


nab_l = nab(l)
nab_l[:]


# In[27]:


div_l = nab_l.trace()
print(div_l)
div_l.display()


# In[28]:


div_l.expr().factor()


# In[29]:


div_l.expr().subs(r=2*m)


# In[30]:


acc_l = l.contract(0,nab_l,1)
print(acc_l)
acc_l.display()


# The non-affinity parameter $\kappa$:

# In[31]:


kappa = l(rho)
kappa.display()


# Check of the pregeodesic equation $\nabla_{\ell} \ell = \kappa \ell$:

# In[32]:


acc_l == kappa * l


# In[33]:


kappa.expr().factor()


# Value of $\kappa$ on the horizon:

# In[34]:


kappaH = kappa.expr().subs(r=2*m)
kappaH


# ## The complementary null vector field $k$

# In[35]:


k = M.vector_field(name='k')
k[0] = 1/2 + m/r
k[1] = -1/2 - m/r
k.display()


# In[36]:


g(k,k).expr()


# In[37]:


g(k,l).expr()


# In[38]:


graph_k = k.plot(ambient_coords=(r,t), ranges={r:(1,8), t:(0,8)}, number_values={r:8, t:9}, 
                 fixed_coords={th:pi/2, ph:pi}, parameters={m:1}, color='red', 
                 scale=0.8, aspect_ratio=1)
graph_lk = graph_l+graph_k
show(graph_lk, figsize=8)
graph_lk.save('def_plot_lk.pdf', figsize=8)


# In[39]:


k_form = k.down(g)
k_form.set_name('kf', latex_name=r'\underline{k}')
k_form.display()


# ## The 2-metric $q$
# 
# We define $q = g + \underline{\ell}\otimes \underline{k} + \underline{k}\otimes \underline{\ell}$:

# In[40]:


q = g + l_form*k_form + k_form*l_form
q.set_name('q')
q.display()


# In[41]:


q_up = q.up(g)
print(q_up)
q_up.display()


# ## Expansion along the null normal

# We compute $\theta_{(\ell)}$ as $\theta_{(\ell)} = q^{\mu\nu}\nabla_\mu \ell_\nu$:

# In[42]:


theta = q_up.contract(0,1,nab(l_form),0,1)
theta.expr()


# Check of the formula $\theta_{(\ell)} = \nabla\cdot\ell - \kappa$:

# In[43]:


theta == div_l - kappa


# Check of the forumla $\theta_{(\ell)} = \frac{1}{2} \mathcal{L}_{\ell} \ln \det q$:

# In[44]:


detq = M.scalar_field({X: r^4*sin(th)^2})
theta == 1/2*ln(detq).lie_der(l)


# In[45]:


theta == 1/2*l(ln(detq))


# ## Deformation rate tensor of the cross-sections
# 
# We compute $\Theta$ as $\Theta = \frac{1}{2} \mathcal{L}_{\ell} q$:

# In[46]:


Theta = 1/2 * q.lie_der(l)
Theta.set_name('Theta', latex_name=r'\Theta')
print(Theta)
Theta.display()


# ## Expansion of the cross-sections along the null normal $k$:
# 
# We compute $\theta_{(k)}$ as $\theta_{(k)} = q^{\mu\nu}\nabla_\mu k_\nu$:

# In[47]:


theta_k = q_up.contract(0,1,nab(k_form),0,1)
theta_k.expr()


# Value of $\theta_{(k)}$ at the horizon:

# In[48]:


theta_k.expr().subs(r=2*m)