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

# # Kerr spacetime
# 
# This worksheet demonstrates a few capabilities of [SageManifolds](http://sagemanifolds.obspm.fr) (version 1.0, as included in SageMath 7.5) in computations regarding Kerr spacetime.
# 
# Click [here](https://raw.githubusercontent.com/sagemanifolds/SageManifolds/master/Worksheets/v1.0/SM_Kerr.ipynb) to download the worksheet file (ipynb format). To run it, you must start SageMath within the Jupyter notebook, via the command `sage -n jupyter`

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

# In[1]:


version()


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

# In[2]:


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


# We also define a viewer for 3D plots (use `'threejs'` or `'jmol'` for interactive 3D graphics):

# In[3]:


viewer3D = 'threejs' # must be 'threejs', jmol', 'tachyon' or None (default)


# Since some computations are quite long, we ask for running them in parallel on 8 cores:

# In[4]:


Parallelism().set(nproc=8)


# ## Spacetime manifold
# 
# We declare the Kerr spacetime as a 4-dimensional diffentiable manifold:

# In[5]:


M = Manifold(4, 'M', r'\mathcal{M}')
print(M)


# Let us use the standard **Boyer-Lindquist coordinates** on it, by first introducing the part $\mathcal{M}_0$ covered by these coordinates and then declaring a chart `BL` (for *Boyer-Lindquist*) on $\mathcal{M}_0$, via the method `chart()`, the argument of which is a string expressing the coordinates names, their ranges (the default is $(-\infty,+\infty)$) and their LaTeX symbols:

# In[6]:


M0 = M.open_subset('M0', r'\mathcal{M}_0')
BL.<t,r,th,ph> = M0.chart(r't r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi') 
print(BL) ; BL


# In[7]:


BL[0], BL[1]


# <h2>Metric tensor</h2>
# 
# <p>The 2 parameters $m$ and $a$ of the Kerr spacetime are declared as symbolic variables:</p>

# In[8]:


var('m, a', domain='real')


# <p>Let us introduce the spacetime metric:</p>

# In[9]:


g = M.lorentzian_metric('g')


# <p>The metric is set by its components in the coordinate frame associated with Boyer-Lindquist coordinates, which is the current manifold's default frame:</p>

# In[10]:


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


# <p>A matrix view of the components with respect to the manifold's default vector frame:</p>

# In[11]:


g[:]


# <p>The list of the non-vanishing components:</p>

# In[12]:


g.display_comp()


# <h2>Levi-Civita Connection</h2>
# 
# <p>The Levi-Civita connection $\nabla$ associated with $g$:</p>

# In[13]:


nabla = g.connection() ; print(nabla)


# <p>Let us verify that the covariant derivative of $g$ with respect to $\nabla$ vanishes identically:</p>

# In[14]:


nabla(g) == 0


# <p>Another view of the above property:</p>

# In[15]:


nabla(g).display()


# <p>The nonzero Christoffel symbols (skipping those that can be deduced by symmetry of the last two indices):</p>

# In[16]:


g.christoffel_symbols_display()


# <h2>Killing vectors</h2>
# <p>The default vector frame on the spacetime manifold is the coordinate basis associated with Boyer-Lindquist coordinates:</p>

# In[17]:


M.default_frame() is BL.frame()


# In[18]:


BL.frame()


# <p>Let us consider the first vector field of this frame:</p>

# In[19]:


xi = BL.frame()[0] ; xi


# In[20]:


print(xi)


# <p>The 1-form associated to it by metric duality is</p>

# In[21]:


xi_form = xi.down(g) ; xi_form.display()


# <p>Its covariant derivative is</p>

# In[22]:


nab_xi = nabla(xi_form) ; print(nab_xi) ; nab_xi.display()


# <p>Let us check that the Killing equation is satisfied:</p>

# In[23]:


nab_xi.symmetrize() == 0


# <p>Similarly, let us check that $\frac{\partial}{\partial\phi}$ is a Killing vector:</p>

# In[24]:


chi = BL.frame()[3] ; chi


# In[25]:


nabla(chi.down(g)).symmetrize() == 0


# <h2>Curvature</h2>
# 
# <p>The Ricci tensor associated with $g$:</p>

# In[26]:


Ric = g.ricci() ; print(Ric)


# <p>Let us check that the Kerr metric is a solution of the vacuum Einstein equation:</p>

# In[27]:


Ric == 0


# <p>Another view of the above property:</p>

# In[28]:


Ric.display()


# <p>The Riemann curvature tensor associated with $g$:</p>

# In[29]:


R = g.riemann() ; print(R)


# <p>Contrary to the Ricci tensor, the Riemann tensor does not vanish; for instance, the component $R^0_{\ \, 123}$ is</p>

# In[30]:


R[0,1,2,3]


# <h3>Bianchi identity</h3>
# 
# <p>Let us check the Bianchi identity $\nabla_p R^i_{\ \, j kl} + \nabla_k R^i_{\ \, jlp} + \nabla_l R^i_{\ \, jpk} = 0$:</p>

# In[31]:


DR = nabla(R) ; print(DR)  #long (takes a while)


# In[32]:


for i in M.irange():
    for j in M.irange():
        for k in M.irange():
            for l in M.irange():
                for p in M.irange():
                    print DR[i,j,k,l,p] + DR[i,j,l,p,k] + DR[i,j,p,k,l] ,


# <p>If the last sign in the Bianchi identity is changed to minus, the identity does no longer hold:</p>

# In[33]:


DR[0,1,2,3,1] + DR[0,1,3,1,2] + DR[0,1,1,2,3] # should be zero (Bianchi identity)


# In[34]:


DR[0,1,2,3,1] + DR[0,1,3,1,2] - DR[0,1,1,2,3] # note the change of the second + to -


# ### Kretschmann scalar
# 
# The tensor $R^\flat$, of components $R_{abcd} = g_{am} R^m_{\ \, bcd}$:

# In[35]:


dR = R.down(g) ; print(dR)


# The tensor $R^\sharp$, of components $R^{abcd} = g^{bp} g^{cq} g^{dr} R^a_{\ \, pqr}$:

# In[36]:


uR = R.up(g) ; print(uR)


# The Kretschmann scalar $K := R^{abcd} R_{abcd}$:

# In[37]:


Kr_scalar = uR['^{abcd}']*dR['_{abcd}']
Kr_scalar.display()


# <p>A variant of this expression can be obtained by invoking the<span style="font-family: courier new,courier;"> factor()</span> method on the coordinate function representing the scalar field in the manifold's default chart:</p>

# In[38]:


Kr = Kr_scalar.coord_function()
Kr.factor()


# <p>As a check, we can compare Kr to the formula given by R. Conn Henry, <a href="http://iopscience.iop.org/0004-637X/535/1/350/fulltext/">Astrophys. J. <strong>535</strong>, 350 (2000)</a>:</p>

# In[39]:


Kr == 48*m^2*(r^6 - 15*r^4*(a*cos(th))^2 + 15*r^2*(a*cos(th))^4 
              - (a*cos(th))^6) / (r^2+(a*cos(th))^2)^6


# <p>The Schwarzschild value of the Kretschmann scalar is recovered by setting $a=0$:</p>

# In[40]:


Kr.expr().subs(a=0)


# <p>Let us plot the Kretschmann scalar for $m=1$ and  $a=0.9$:</p>

# In[41]:


K1 = Kr.expr().subs(m=1, a=0.9)
plot3d(K1, (r,1,3), (th, 0, pi), viewer=viewer3D, axes_labels=['r', 'theta', 'Kr'])


# In[ ]: