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

# # Kerr-Schild coordinates on Kerr spacetime
# 
# 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](http://sagemanifolds.obspm.fr) project.
# 
# *NB:* a version of SageMath at least equal to 8.8 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')


# To speed up computations, we ask for running them in parallel on 8 threads:

# In[3]:


Parallelism().set(nproc=8)


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

# In[4]:


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


# ## 3+1 Kerr coordinates $(t,r,\theta,\phi)$

# We restrict the 3+1 Kerr patch to $r>0$, in order to introduce latter the Kerr-Schild coordinates:

# In[5]:


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


# In[6]:


X.coord_range()


# The Kerr parameters $m$ and $a$:

# In[7]:


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


# ## Kerr metric
# 
# We define the metric $g$ by its components w.r.t. the 3+1 Kerr coordinates:

# In[8]:


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


# In[9]:


g.display_comp()


# The inverse metric is pretty simple:

# In[10]:


g.inverse()[:]


# as well as the determinant w.r.t. to the 3+1 Kerr coordinates:

# In[11]:


g.determinant().display()


# In[12]:


g.determinant() == - (rho2*sin(th))^2


# ## Ingoing principal null geodesics

# In[13]:


k = M.vector_field(1, -1, 0, 0, name='k')
k.display()


# Let us check that $k$ is a null vector:

# In[14]:


g(k,k).display()


# Computation of $\nabla_k k$:

# In[15]:


nabla = g.connection()
acc = nabla(k).contract(k)
acc.display()


# In[16]:


k_form = k.down(g)
k_form.display()


# ## Kerr-Schild form of the Kerr metric

# Let us introduce the metric $f$ such that
# $$ g = f + 2 H \underline{k} \otimes \underline{k}$$
# where $H = m r / \rho^2$:

# In[17]:


H = M.scalar_field({X: m*r/rho2}, name='H')
H.display()


# In[18]:


f = M.lorentzian_metric('f')
f.set(g - 2*H*(k_form*k_form))
f.display()


# In[19]:


f[:]


# $f$ is a flat metric:

# In[20]:


f.riemann().display()


# which proves that $g$ is a Kerr-Schild metric.

# Let us check that $k$ is a null vector for $f$ as well:

# In[21]:


f(k,k).expr()


# ## Kerr-Schild coordinates $(t, x, y, z)$

# In[22]:


KS.<t,x,y,z> = M.chart()
KS


# In[23]:


X_to_KS = X.transition_map(KS, [t, 
                                (r*cos(ph) - a*sin(ph))*sin(th),
                                (r*sin(ph) + a*cos(ph))*sin(th),
                                r*cos(th)])
X_to_KS.display()


# In[24]:


R = sqrt((x^2 + y^2 + z^2 - a^2 
          + sqrt((x^2 + y^2 + z^2 - a^2)^2 + 4*a^2*z^2))/2)
R


# In[25]:


#X_to_KS.set_inverse(t, R, acos(z/R), 
#                    atan2(R*y - a*x, R*x + a*y))


# Check of the identity
# $$\frac{x^2 + y^2}{r^2 + a^2} + \frac{z^2}{r^2} = 1$$

# In[26]:


((x^2 + y^2)/(R^2 + a^2) + z^2/R^2).simplify_full()


# ### Minkowskian expression of $f$ in terms of Kerr-Schild coordinates:

# In[27]:


f.display(KS.frame())


# Equivalently, we may check the following identity:

# In[28]:


dt, dx, dy, dz = KS.coframe()[:]
f == - dt*dt + dx*dx + dy*dy + dz*dz


# In[29]:


dx.display()


# In[30]:


(dx*dx + dy*dy + dz*dz).display()


# ### Expression of $k$ and $g$ in the Kerr-Schild frame:

# In[31]:


k.display(KS.frame())


# In[32]:


k_form.display(KS.frame())


# In[33]:


g.display(KS.frame())


# Expression of the Killing vector $\partial/\partial\phi$ in terms of the Kerr-Schild frame:

# In[34]:


X.frame()[3].display(KS.frame())


# ## Plots

# In[35]:


ap = 0.9  # value of a for the plots
rmax = 3


# In[36]:


rcol = 'green'       # color of the curves (th,ph) = const
thcol = 'red'        # color of the curves (r,ph) = const
phcol = 'goldenrod'  # color of the curves (r,th) = const
coordcol = {r: rcol, th: thcol, ph: phcol}


# In[37]:


opacity = 1
surf_shift = 0.03


# Numerical values of the event and Cauchy horizons:

# In[38]:


rHp = 1 + sqrt(1 - ap^2)
rCp = 1 - sqrt(1 - ap^2)
rHp, rCp


# In[39]:


X_plot = X.plot(KS, fixed_coords={t: 0, ph: 0}, ambient_coords=(x,y,z), 
                parameters={a: ap}, number_values=11,
                color=coordcol, thickness=2, max_range=rmax, 
                label_axes=False)


# In[40]:


X_to_KS(t, r, th, ph)


# In[41]:


xyz_n = [s.subs({a: ap, ph: 0}) for s in X_to_KS(t, r, th, ph)[1:]]
xyz_n


# In[42]:


xyz_H = [s.subs({r: rHp}) for s in xyz_n]
xyz_H


# In[43]:


xyz_C = [s.subs({r: rCp}) for s in xyz_n]
xyz_C


# In[44]:


xyz_n[1] += surf_shift  # small adjustment to ensure that the surface does not cover
                        # the coordinate grid lines


# In[45]:


graph = parametric_plot3d(xyz_n, (r, 0, rmax), (th, 0, pi), color='ivory',
                          opacity=opacity)
graph += X_plot
graph += parametric_plot3d(xyz_H, (th, 0, pi), color='black', thickness=6)
graph += parametric_plot3d(xyz_C, (th, 0, pi), color='blue', thickness=6)
graph += line([(0.03,0,0), (0.03,ap,0)], color='red', thickness=6)
graph


# In[46]:


xyz_n = [s.subs({a: ap, ph: pi}) for s in X_to_KS(t, r, th, ph)[1:]]
xyz_H = [s.subs({r: rHp}) for s in xyz_n]
xyz_C = [s.subs({r: rCp}) for s in xyz_n]
X_plot_pi = X.plot(KS, fixed_coords={t: 0, ph: pi}, ambient_coords=(x,y,z), 
                parameters={a: ap}, number_values=11,
                color=coordcol, thickness=2, max_range=rmax, 
                label_axes=False)


# In[47]:


xyz_n[1] += surf_shift  # small adjustment to ensure that the surface does not cover
                        # the coordinate grid lines
graph_pi = parametric_plot3d(xyz_n, (r, 0, rmax), (th, 0, pi), color='ivory',
                             opacity=opacity)
graph_pi += X_plot_pi
graph_pi += parametric_plot3d(xyz_H, (th, 0, pi), color='black', thickness=6)
graph_pi += parametric_plot3d(xyz_C, (th, 0, pi), color='blue', thickness=6)
graph_pi += line([(0.03,0,0), (0.03,-ap,0)], color='red', thickness=6)
graph_pi += line([(0,0,-1.1*rmax), (0,0,1.1*rmax)], color='green', thickness=4)
ph_slice = graph + graph_pi
show(ph_slice)


# In[48]:


X.plot(KS, fixed_coords={t: 0, r: 1}, ambient_coords=(x,y,z), 
       parameters={a: ap}, number_values=11,
       color=coordcol, label_axes=False)


# The BH event horizon:

# In[49]:


H_plot = X.plot(KS, fixed_coords={t: 0, r: rHp}, ambient_coords=(x,y,z), 
                parameters={a: ap}, number_values=11, color='grey', 
                label_axes=False)
H_plot


# In[50]:


X.plot(KS, fixed_coords={t: 0, r: 0}, ambient_coords=(x,y,z), 
       parameters={a: ap}, number_values=11, color=coordcol, 
       label_axes=False)


# In[51]:


rzero = X.plot(KS, fixed_coords={t: 0, r: 0}, ambient_coords=(x,y), 
               parameters={a: ap}, number_values={th: 17, ph: 13}, color=coordcol)
rzero += circle((0,0), ap, color='brown', thickness=3) 
show(rzero, xmin=-1, xmax=1, ymin=-1, ymax=1, axes=False, frame=True,
     axes_labels=[r'$x/m$', r'$y/m$'])


# In[52]:


rzero.save('ksm_rzero_disk.pdf', xmin=-1, xmax=1, ymin=-1, ymax=1, 
           axes=False, frame=True, axes_labels=[r'$x/m$', r'$y/m$'])


# In[53]:


Xth_pi2 = X.plot(KS, fixed_coords={t: 0, th: pi/2}, ambient_coords=(x,y,z), 
                 parameters={a: ap}, number_values=11, color=coordcol,
                 max_range=rmax, thickness=1.5, label_axes=False)
Xth_pi2


# In[54]:


X.plot(KS, fixed_coords={t: 0, th: pi/3}, ambient_coords=(x,y,z), 
       parameters={a: ap}, number_values=11, color=coordcol,
       max_range=rmax, thickness=1.5, label_axes=False)


# In[55]:


graph = X.plot(KS, fixed_coords={t: 0, th: pi/6}, ambient_coords=(x,y,z), 
               parameters={a: ap}, number_values=11, color=coordcol,
               max_range=rmax, thickness=1.5, label_axes=False) \
        + Xth_pi2 \
        + circle((0,0,0), ap, color='pink', fill=True) \
        + circle((0,0,0), ap, color='brown', thickness=6, linestyle='dashed')
show(graph)


# ## Plot of the $r<0$ domain

# In[56]:


KS2.<t,xp,yp,zp> = M.chart(r"t xp:x' yp:y' zp:z'")
KS2


# We use replace $r$ by $-r$ in the transformation formulas, because in what follows, we consider that
# $r$ is positive:

# In[57]:


X_to_KS2 = X.transition_map(KS2, [t, 
                                 (-r*cos(ph) - a*sin(ph))*sin(th),
                                 (-r*sin(ph) + a*cos(ph))*sin(th),
                                 -r*cos(th)])
X_to_KS2.display()


# In[58]:


X_plot2 = X.plot(KS2, fixed_coords={t: 0, ph: 0}, ambient_coords=(xp,yp,zp), 
                 parameters={a: ap}, number_values=11,
                 color=coordcol, thickness=2, max_range=rmax, 
                 label_axes=False)


# In[59]:


surf_shift = 0
xyz_n2 = [s.subs({a: ap, ph: 0}) for s in X_to_KS2(t, r, th, ph)[1:]]
xyz_n2[1] -= surf_shift 


# In[60]:


graph2 = parametric_plot3d(xyz_n2, (r, 0, rmax), (th, 0, pi), color='pink',
                           opacity=opacity)
graph2 += X_plot2


# In[61]:


X_plot2_pi = X.plot(KS2, fixed_coords={t: 0, ph: pi}, ambient_coords=(xp,yp,zp), 
                parameters={a: ap}, number_values=11,
                color=coordcol, thickness=2, max_range=rmax, 
                label_axes=False)


# In[62]:


xyz_n2 = [s.subs({a: ap, ph: pi}) for s in X_to_KS2(t, r, th, ph)[1:]]
xyz_n2[1] += surf_shift


# In[63]:


graph2_pi = parametric_plot3d(xyz_n2, (r, 0, rmax), (th, 0, pi), color='pink',
                              opacity=opacity)
graph2_pi += X_plot2_pi


# In[64]:


ph_slice2 = graph2 + graph2_pi
ph_slice2


# In[65]:


ph_slice + ph_slice2


# In[ ]: