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

# # 3-sphere: the round metric
# 
# This worksheet demonstrates a few capabilities of
# [SageManifolds](http://sagemanifolds.obspm.fr) (version 1.0.2, as included in SageMath 8.0)
# on the example of the 3-dimensional sphere, $\mathbb{S}^3$.
# 
# Click [here](https://raw.githubusercontent.com/sagemanifolds/SageManifolds/master/Worksheets/v1.0/SM_sphere_S3_metric.ipynb) to download the worksheet file (ipynb format). To run it, you must start SageMath with the Jupyter notebook, via the command `sage -n jupyter`

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

# 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 increase the computational speed, we ask for demanding computations to be parallelly performed on 8 cores:

# In[3]:


Parallelism().set(nproc=8)


# ## $\mathbb{S}^3$ as a 3-dimensional differentiable manifold
# 
# We start by declaring $\mathbb{S}^3$ as a differentiable manifold of dimension 3 over $\mathbb{R}$:

# In[4]:


S3 = Manifold(3, 'S^3', latex_name=r'\mathbb{S}^3', start_index=1)


# The first argument, `3`, is the dimension of the manifold, while the second argument is the symbol used to label the manifold, with the LaTeX output specified by the argument `latex_name`. The argument `start_index` sets the index range to be used on the manifold for labelling components w.r.t. a basis or a frame: `start_index=1` corresponds to $\{1,2,3\}$; the default value is `start_index=0`, yielding to $\{0,1,2\}$.

# In[5]:


print(S3)


# In[6]:


S3


# ### Coordinate charts on $\mathbb{S}^3$
# 
# The 3-sphere cannot be covered by a single chart. At least two charts are necessary, for instance the charts associated with the stereographic projections from two distinct points, $N$ and $S$ say,
# which we may call the *North pole* and the *South pole* respectively. Let us introduce the open subsets covered by these two charts: 
# $$ U := \mathbb{S}^3\setminus\{N\} $$  
# $$ V := \mathbb{S}^3\setminus\{S\} $$

# In[7]:


U = S3.open_subset('U') ; print(U)


# In[8]:


V = S3.open_subset('V') ; print(V)


# We declare that $\mathbb{S}^3 = U \cup V$:

# In[9]:


S3.declare_union(U, V)


# Then we introduce the stereographic chart on $U$, denoting by $(x,y,z)$ the coordinates resulting from the stereographic projection from the North pole onto the equatorial plane:

# In[10]:


stereoN.<x,y,z> = U.chart()
stereoN


# In[11]:


stereoN.coord_range()


# Similarly, we introduce on $V$ the coordinates $(x',y',z')$ corresponding to the stereographic projection from the South pole onto the equatorial plane:

# In[12]:


stereoS.<xp,yp,zp> = V.chart("xp:x' yp:y' zp:z'")
stereoS


# In[13]:


stereoS.coord_range()


# We have to specify the **transition map** between the charts `stereoN` = $(U,(x,y,z))$ and `stereoS` = $(V,(x',y',z'))$; it is given by the standard inversion formulas:

# In[14]:


r2 = x^2+y^2+z^2
stereoN_to_S = stereoN.transition_map(stereoS, 
                                      (x/r2, y/r2, z/r2), 
                                      intersection_name='W',
                                      restrictions1= x^2+y^2+z^2!=0, 
                                      restrictions2= xp^2+yp^2+zp^2!=0)
stereoN_to_S.display()


# In the above declaration, `'W'` is the name given to the open subset where the two charts overlap: $W := U\cap V$, the condition $x^2+y^2+z^2\not=0$  defines $W$ as a subset of $U$, and the condition $x'^2+y'^2+z'^2\not=0$ defines $W$ as a subset of $V$.
# 
# The inverse coordinate transformation is computed by means of the method `inverse()`:

# In[15]:


stereoS_to_N = stereoN_to_S.inverse()
stereoS_to_N.display()


# Note that the situation is of course perfectly symmetric regarding the coordinates $(x,y,z)$ and $(x',y',z')$.
# 
# At this stage, the user's atlas has four charts:

# In[16]:


S3.atlas()


# For future reference, we store $W=U\cap V$ into a Python variable:

# In[17]:


W = U.intersection(V)
print(W)


# ### The North and South poles
# 
# $N$ is the point of $V$ of stereographic coordinates $(x',y',z')=(0,0,0)$:

# In[18]:


N = V((0,0,0), chart=stereoS, name='N')
print(N)


# while $S$ is the point of $U$ of stereographic coordinates $(x,y,z)=(0,0,0)$:

# In[19]:


S = U((0,0,0), chart=stereoN, name='S')
print(S)


# We have of course

# In[20]:


all([N not in U, N in V, S in U, S not in V])


# ## Embedding of $\mathbb{S}^3$ into $\mathbb{R}^4$
# 
# Let us first declare $\mathbb{R}^4$ as a 4-dimensional manifold covered by a single chart (the so-called **Cartesian coordinates**):

# In[21]:


R4 = Manifold(4, 'R^4', r'\mathbb{R}^4')
X4.<T,X,Y,Z> = R4.chart()
X4


# The embedding of $\mathbb{S}^3$ into $\mathbb{R}^4$ is then defined by the standard formulas relating the stereographic coordinates to the ambient Cartesian ones when considering a **stereographic projection** from the point $(-1,0,0,0)$ to the equatorial plane $T=0$:

# In[22]:


rp2 = xp^2 + yp^2 + zp^2
Phi = S3.diff_map(R4, {(stereoN, X4): 
                       [(1-r2)/(r2+1), 2*x/(r2+1), 
                        2*y/(r2+1), 2*z/(r2+1)],
                       (stereoS, X4):
                       [(rp2-1)/(rp2+1), 2*xp/(rp2+1), 
                        2*yp/(rp2+1), 2*zp/(rp2+1)]},
                  name='Phi', latex_name=r'\Phi')
Phi.display()


# ## Hyperspherical coordinates
# 
# The hyperspherical coordinates $(\chi, \theta, \phi)$ generalize the standard spherical coordinates $(\theta, \phi)$ on $\mathbb{S}^2$. They are defined on the open domain $A\subset W \subset \mathbb{S}^3$ that is the complement of the "origin meridian"; since the latter is defined by $y=0$ and $x\geq 0$, we declare:

# In[23]:


A = W.open_subset('A', coord_def={stereoN.restrict(W): (y!=0, x<0), 
                                  stereoS.restrict(W): (yp!=0, xp<0)})
print(A)


# In[24]:


spher.<ch,th,ph> = A.chart(r'ch:(0,pi):\chi th:(0,pi):\theta ph:(0,2*pi):\phi')
spher


# In[25]:


den = 1 + cos(ch)
spher_to_stereoN = spher.transition_map(stereoN.restrict(A), 
                                        (sin(ch)*sin(th)*cos(ph)/den,
                                         sin(ch)*sin(th)*sin(ph)/den,
                                         sin(ch)*cos(th)/den))
spher_to_stereoN.display()


# In[26]:


spher_to_stereoN.set_inverse(2*atan(sqrt(x^2+y^2+z^2)),
                             atan2(sqrt(x^2+y^2), z),
                             atan2(-y, -x)+pi,
                             verbose=True)


# In[27]:


spher_to_stereoN.inverse().display()


# In[28]:


spher_to_stereoS = stereoN_to_S.restrict(A) * spher_to_stereoN
spher_to_stereoS.display()


# In[29]:


stereoS_to_spher = spher_to_stereoN.inverse() * stereoS_to_N.restrict(A)
stereoS_to_spher.display()


# In[30]:


Phi.display(stereoN.restrict(A), X4)


# In[31]:


Phi.display(spher, X4)


# In[32]:


Phi.display()


# ## Global left-invariant vector frame on $\mathbb{S}^3$

# The vector frames associated with the two stereographic charts are

# In[33]:


frameN = stereoN.frame()
frameN


# In[34]:


frameS = stereoS.frame()
frameS


# None of these two frames cover entirely the 3-sphere, since $U$ and $V$ are strict subsets of $\mathbb{S}^3$. Now, as it is well known, $\mathbb{S}^3$ admits global vector frames, i.e. $\mathbb{S}^3$ is a **parallelizable manifold**. Among all the spheres, it shares this remarkable property with $\mathbb{S}^1$ and $\mathbb{S}^7$. We shall use a global vector frame $(\mathbb{S}^3, (\epsilon_1, \epsilon_2, \epsilon_3))$ associated with the Lie group structure of $\mathbb{S}^3$, namely the left-invariant vector frame constructed in the worksheet 
# [3-sphere: vector fields and left-invariant parallelization](http://nbviewer.jupyter.org/github/sagemanifolds/SageManifolds/blob/master/Worksheets/v1.0/SM_sphere_S3_vectors.ipynb). We first declare this vector frame on all $\mathbb{S}^3$:

# In[35]:


E = S3.vector_frame('E', latex_symbol=r'\varepsilon')
E


# On $U$, we relate this frame to the stereographic coordinate frame, by means of the formulas obtained in the vector field worksheet mentionned above:

# In[36]:


E_U = E.restrict(U); E_U


# In[37]:


E_U[1][frameN,:,stereoN] = \
    [(x^2-y^2-z^2+1)/2, x*y+z, x*z-y]
E_U[2][frameN,:,stereoN] = \
    [x*y-z, (1-x^2+y^2-z^2)/2, x+y*z]  
E_U[3][frameN,:,stereoN] = \
    [x*z+y, y*z-x, (1-x^2-y^2+z^2)/2]


# Similarly, on $V$, we relate the global frame with the stereographic coordinate frame from the South pole:

# In[38]:


E_V = E.restrict(V); E_V


# In[39]:


E_V[1][frameS,:, stereoS] = \
    [(yp^2+zp^2-xp^2-1)/2, zp-xp*yp, -yp-xp*zp]
E_V[2][frameS,:, stereoS] = \
    [-zp-xp*yp, (xp^2-yp^2+zp^2-1)/2, xp-yp*zp]
E_V[3][frameS,:, stereoS] = \
    [yp-xp*zp, -xp-yp*zp, (xp^2+yp^2-zp^2-1)/2]


# Let us display the links between the various frames:

# In[40]:


for i in S3.irange():
    show(E[i].display(frameN))
print(" ")
for i in S3.irange():
    show(E[i].display(frameS))


# To complete the links, we introduce the change-of-frame operators: first the operator $P$ such that, on $U$, $\epsilon_i = P(\partial/\partial{x^i})$, with $x^i=(x,y,z)$:

# In[41]:


P = U.automorphism_field()
for i in S3.irange():
    for j in S3.irange():
        P[j,i] = E_U[i][j]


# In[42]:


all([E_U[i] == P(frameN[i]) for i in S3.irange()])


# We add $P$ to the known changes of frame on $U$:

# In[43]:


U.set_change_of_frame(frameN, E_U, P)


# Similarly, on $V$, we introduce the operator $P$ such that $\epsilon_i = P(\partial/\partial {x'}^i)$, with ${x'}^i=(x',y',z')$:

# In[44]:


P = V.automorphism_field()
for i in S3.irange():
    for j in S3.irange():
        P[j,i] = E_V[i][j]


# In[45]:


all([E_V[i] == P(frameS[i]) for i in S3.irange()])


# In[46]:


V.set_change_of_frame(frameS, E_V, P)


# ## Hopf coordinates

# Despite there does not exist any coordinate chart associated to the global vector frame $(\mathbb{S}^3, (\epsilon_1, \epsilon_2, \epsilon_3))$, there exist a chart, called the
# **Hopf chart**, which has some link to it in the sense that some coordinate lines are integral curves of $\epsilon_3$. The Hopf coordinates have been introduced in the worksheet [3-sphere: charts, quaternions and Hopf fibration](http://nbviewer.jupyter.org/github/sagemanifolds/SageManifolds/blob/master/Worksheets/v1.0/SM_sphere_S3_Hopf.ipynb): 

# In[47]:


B = U.open_subset('B', coord_def={stereoN.restrict(U): 
                                   [x^2+y^2!=0, x^2+y^2+z^2!=1, 
                                    (1-x^2-y^2-z^2)*x-2*y*z!=0]})
print(B)


# In[48]:


Hcoord.<eta,alp,bet> = B.chart(r"eta:(0,pi/2):\eta alpha:(0,2*pi):\alpha beta:(0,2*pi):\beta")
Hcoord


# In[49]:


Hcoord_to_stereoN = Hcoord.transition_map(
                        stereoN.restrict(U),
                        (sin(eta)*cos(alp+bet)/(1+cos(eta)*sin(alp)),
                         sin(eta)*sin(alp+bet)/(1+cos(eta)*sin(alp)),
                         cos(eta)*cos(alp)/(1+cos(eta)*sin(alp))))
Hcoord_to_stereoN.display()


# In[50]:


Hcoord_to_stereoN.set_inverse(asin(2*sqrt(x^2+y^2)/(1+x^2+y^2+z^2)),
                              atan2(x^2+y^2+z^2-1, -2*z) + pi,
                              atan2(-y,-x) - atan2(x^2+y^2+z^2-1, -2*z),
                              verbose=True)


# In[51]:


Hcoord_to_stereoN.inverse().display()


# As explained in the worksheet [3-sphere: vector fields and left-invariant parallelization](http://nbviewer.jupyter.org/github/sagemanifolds/SageManifolds/blob/master/Worksheets/v1.0/SM_sphere_S3_vectors.ipynb), due to some lack of simplification, it's better to set the components of $\varepsilon_1$ and $\varepsilon_2$ by hand:

# In[52]:


E_U[1].add_comp(Hcoord.frame())[1, Hcoord] = sin(2*alp+bet)
E_U[1].add_comp(Hcoord.frame())[2, Hcoord] = -cos(2*alp+bet) * tan(eta)
E_U[1].add_comp(Hcoord.frame())[3, Hcoord] = cos(2*alp+bet) / (cos(eta)*sin(eta))
E_U[1].display(Hcoord.frame(), Hcoord)


# In[53]:


E_U[2].add_comp(Hcoord.frame())[1, Hcoord] = -cos(2*alp+bet)
E_U[2].add_comp(Hcoord.frame())[2, Hcoord] = -sin(2*alp+bet) * tan(eta)
E_U[2].add_comp(Hcoord.frame())[3, Hcoord] = sin(2*alp+bet) / (cos(eta)*sin(eta))
E_U[2].display(Hcoord.frame(), Hcoord)


# On the contrary, the components of $\varepsilon_3$ are particularly simple:

# In[54]:


E_U[3].display(Hcoord.frame(), Hcoord)


# ## Round metric on $\mathbb{S}^3$
# 
# The standard metric on $\mathbb{S}^3$ is that induced by the Euclidean metric of $\mathbb{R}^4$. Let us start by defining the latter:

# In[55]:


h = R4.metric('h')
h[0,0], h[1,1], h[2,2], h[3, 3] = 1, 1, 1, 1
h.display()


# The metric $g$ on $\mathbb{S}^3$ is the pullback of $h$ by the embedding $\Phi$:

# In[56]:


g = S3.metric('g')
g.set( Phi.pullback(h) )
print(g)


# Let us display $g$ in terms of the coordinate frame associated with the stereographic chart from the North pole $(U,(x,y,z))$:

# In[57]:


g.display(frameN)


# The components can be factored:

# In[58]:


for i in S3.irange():
    g[frameN, i,i].factor()
g.display(frameN)


# In[59]:


g.display(frameS)


# In[60]:


for i in S3.irange():
    g[frameS, i,i].factor()
g.display(frameS)


# ### Expression of the round metric in the global frame
# 
# The expression of $g$ in terms of the global frame 
# $(\varepsilon_1, \varepsilon_2, \varepsilon_3)$ is deduced from that w.r.t. the stereographic coordinates on the two domains $U$ and $V$, where the change-of-frame formulas are known:

# In[61]:


g.display(E_U)


# In[62]:


g.display(E_V)


# We may then set the components globally:

# In[63]:


for i in S3.irange():
    g.add_comp(E)[i, i] = 1
g.display(E)


# The above shows that the global frame $(\varepsilon_1, \varepsilon_2, \varepsilon_3)$ is **orthonormal** (w.r.t. to $g$), which is not surprising since this frame is induced by 
# three vector fields of $\mathbb{R}^4$ which are clearly orthonormal for the Euclidean metric $h$, namely the three vector fields $E_{\mathbf{i}}$, $E_{\mathbf{j}}$ and $E_{\mathbf{k}}$
# in the worksheet [3-sphere: vector fields and left-invariant parallelization](http://nbviewer.jupyter.org/github/sagemanifolds/SageManifolds/blob/master/Worksheets/v1.0/SM_sphere_S3_vectors.ipynb).

# ### Expression of the round metric in terms of the hyperspherical coordinates

# In[64]:


g.display(spher.frame(), spher)


# ### Expression of the round metric in terms of the Hopf coordinates

# In[65]:


g.display(Hcoord.frame(), Hcoord)


# We note that the components of $g$ depend only on $\eta$. This implies that the tori $\eta=\mathrm{const}$ are *flat* tori.

# ### The Christoffel symbols

# This is necessary for what follows (should be fixed in a future version):

# In[66]:


for chart in S3.top_charts():
    g.restrict(chart.domain()).inverse()


# In[67]:


g.christoffel_symbols_display(stereoN)


# In[68]:


g.christoffel_symbols_display(spher)


# In[69]:


g.christoffel_symbols_display(Hcoord)


# ### The Riemann tensor

# In[70]:


Riem = g.riemann()
print(Riem)


# In[71]:


Riem.display_comp(frameN, stereoN, only_nonredundant=True)


# In[72]:


Riem.display_comp(spher.frame(), spher, only_nonredundant=True)


# In[73]:


Riem.display_comp(Hcoord.frame(), Hcoord, only_nonredundant=True)


# In[74]:


Riem.display_comp(E_U, stereoN, only_nonredundant=True)


# In[75]:


Riem.display(E_U)


# ### The Ricci tensor

# In[76]:


Ric = g.ricci()
print(Ric)


# In[77]:


Ric.display(frameN)


# In[78]:


for i in S3.irange():
    Ric[frameN, i,i].factor()
Ric.display(frameN)


# We have $Ric(g) = 2 g$:

# In[79]:


Ric.restrict(U) == 2 * g.restrict(U)


# ### The Ricci scalar

# In[80]:


R_A = g.restrict(A).ricci_scalar()
print(R_A)


# In[81]:


R_A.display()


# We note that the Ricci scalar is constant, as for any maximally symmetric space.

# ## The volume 3-form (Levi-Civita tensor)
# 
# We have

# In[82]:


eps = g.volume_form()
print(eps)
eps.display()


# In[83]:


print(eps.parent())


# ### Expression of the volume 3-form in terms of the stereographic coordinates

# In[84]:


eps.display(frameN, stereoN)


# The component $\epsilon_{xyz}$ can be factored:

# In[85]:


eps_xyz = eps[frameN, 1,2,3, stereoN]
eps_xyz


# In[86]:


eps_xyz.factor()


# In[87]:


eps.display(frameN, stereoN)