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

# # Slicing of Kerr spacetime by spacelike submanifolds
# 
# This notebook demonstrates a few capabilities of SageMath in computations regarding embedded submanifolds on the specific example of the 3+1 slicing of Kerr spacetime. The corresponding tools have been developed within the  [SageManifolds](https://sagemanifolds.obspm.fr) project (version 1.3, as included in SageMath 8.3).
# 
# Click [here](https://raw.githubusercontent.com/sagemanifolds/SageManifolds/master/Worksheets/v1.3/SM_submanifold_Kerr_slicing.ipynb) to download the notebook 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 8.3 is required to run this notebook:

# In[1]:


version()


# The first cell defines some notebook options. Parallelism is usefull to speedup calculation of the connection.

# In[2]:


get_ipython().run_line_magic('display', 'latex')
Parallelism().set(nproc=8)


# ## Manifolds declaration
# 
# Next cells are used to define the manifolds. In our case, $M$ will be the Kerr spacetime, of which $N$ is a single spacelike submanifold of dimension 3 (hypersurface).

# In[3]:


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


# In[4]:


N = Manifold(3, 'N', ambient=M, structure="Riemannian")
print(N)


# Note the `ambient=M`, used to declare $N$ as a submanifold of $M$.
# 
# We will use a single chart, corresponding to the so-called **Kerr coordinates** (*NB*: they differ from the standard Boyer-Lindquist coordinates):

# In[5]:


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


# ### Spacetime metric

# The metric depends on two quantities: the mass of the black-hole and its angular momentum. These are declared as symbolic variables. $\rho$ is a shortcut notation.

# In[6]:


m = var('m') # mass
a = var('a') # ratio J/m


# In[7]:


rho = sqrt(r^2+a^2*cos(th)^2)


# We set the metric components in the Kerr coordinates introduced above:

# In[8]:


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


# In[9]:


g[:]


# Let's compute the Levi-Civita connection:

# In[10]:


nab = g.connection()


# List of Christoffel symbols:

# In[11]:


nab.display(only_nonredundant=True)


# Ricci tensor:

# In[12]:


nab.ricci()[:]


# This proves that the metric given in cell [8] is indeed solution of the vacuum Einstein equation, as expected of a Kerr metric.
# 
# It's now time to foliate it along the $t$ axis.

# ## Slicing:
# 
# To compute 3+1 quantities, we must first declare an embedding from $N$ to $M$. because we want our foliation to encompass the whole manifold, the embedding must depend on a parameter $\tau$.

# In[13]:


tau = var(r'tau') # foliation parameter


# Let's also declare a chart on $N$.

# In[14]:


CN.<r0,th0,ph0> = N.chart(r'r0:(0,+oo) th0:(0,pi):\theta_0 ph0:(-pi,pi):\phi_0')


# The embedding used is a very simple one, it simply sets $t = \tau$.

# In[15]:


phi = N.diff_map(M, {(CN,C):[tau,r0,th0,ph0]})
phi.display()


# The two partial inverses are also quite simple:

# In[16]:


phi_inv = M.diff_map(N, {(C,CN):[r,th,ph]})
phi_inv.display()


# In[17]:


phi_inv_tau = M.scalar_field({C:t})
phi_inv_tau.display()


# The following line registers the embedding, as well as the inverses. This cannot be done during the creation of the manifolds because the embedding `phi` needs exisiting domain and codomain to be created.

# In[18]:


N.set_embedding(phi, inverse=phi_inv, var=tau, t_inverse={tau: phi_inv_tau})


# Computing the adapted chart is needed for computing the shift and extrinsic curvature. This create another chart, which in our case is identical to `CN` because of the form of the embedding.

# In[19]:


T = N.adapted_chart()
T


# The induced metric on $N$ (or first fundamental form) is then given by:

# In[20]:


N.induced_metric()[:]


# And its inverse:

# In[21]:


N.induced_metric().inverse()[:]


# The normal vector is given by

# In[22]:


N.normal().display()


# Let's check that it is correcly normalized:

# In[23]:


g(N.normal(), N.normal()).display()


# The lapse and shift are also easily computed: 

# In[24]:


N.lapse().display()


# In[25]:


N.shift().display()


# The extrinsic curvature (or second fundamental form) is long to compute because of failed simplifications (see the last component).

# In[26]:


N.second_fundamental_form()[:]


# ## Example of 3+1 projection
# 
# It is also possible to project any tensor on $M$ orhogonaly on the submnifold or on the normal vector. This can be useful for example to decompose an energy-impulsion tensor.
# 
# Let's construct a fourth-order tensor from a simple vector field.

# In[27]:


v = M.vector_field()
v.add_comp()[:] = [1,1,1,1]


# The fourth-order tensor $v\otimes v\otimes v\otimes v$ is then projected twice on the submanifold (on the indices 1 and 3) and twice along the normal vector (indicies 0 and 2). The result is the following second-order tensor on $M$:

# In[28]:


N.mixed_projection(v*v*v*v, [0,2]).display()


# In[ ]: