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

# # 3-sphere: charts, quaternions and Hopf fibration
# 
# This notebook demonstrates some differential geometry capabilities of SageMath on the example of the 3-dimensional sphere, $\mathbb{S}^3$. The corresponding tools have been developed within the  [SageManifolds](https://sagemanifolds.obspm.fr) project.
# 
# Click [here](https://raw.githubusercontent.com/sagemanifolds/SageManifolds/master/Notebooks/SM_sphere_S3_Hopf.ipynb) to download the notebook 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 7.5 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')


# and we initialize a time counter for benchmarking:

# In[3]:


import time
comput_time0 = time.perf_counter()


# To increase the computational speed, we ask for demanding computations to be parallelly performed on 8 threads:

# In[4]:


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[5]:


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[6]:


print(S3)


# In[7]:


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[8]:


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


# In[9]:


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


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

# In[10]:


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[11]:


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


# In[12]:


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[13]:


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


# In[14]:


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[15]:


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[16]:


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[17]:


S3.atlas()


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

# In[18]:


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[19]:


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[20]:


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


# We have of course

# In[21]:


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[22]:


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[23]:


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()


# From this choice of stereographic projection, the "North" pole is actually the point of coordinates $(-1,0,0,0)$ in $\mathbb{R}^4$:

# In[24]:


X4(Phi(N))


# while the "South" pole is the point of coordinates $(1,0,0,0)$:

# In[25]:


X4(Phi(S))


# We may use the embedding $\Phi$ to plot the stereographic coordinate grid in terms of the $\mathbb{R}^4$'s Cartesian coordinates: 

# In[26]:


graph_stereoN = stereoN.plot(chart=X4, mapping=Phi, 
                             ambient_coords=(X,Y,Z),
                             number_values=9,
                             color={x: 'red', y: 'green', z: 'gold'},
                             label_axes=False)
show(graph_stereoN, axes_labels=['X', 'Y', 'Z'])


# In[27]:


graph_stereoN = stereoN.plot(chart=X4, mapping=Phi, 
                             ambient_coords=(X,Y,T),
                             number_values=13, plot_points=150,
                             color={x: 'red', y: 'green', z: 'gold'},
                             label_axes=False)
pointN = N.plot(chart=X4, mapping=Phi, ambient_coords=(X,Y,T), 
                color='maroon', label_offset=0.05)
pointS = S.plot(chart=X4, mapping=Phi, ambient_coords=(X,Y,T), 
                color='maroon', label_offset=0.05)
show(graph_stereoN + pointN + pointS, axes_labels=['X', 'Y', 'T'])


# ## 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[28]:


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


# We then declare the chart $(A,(\chi,\theta,\phi))$ by specifying the intervals spanned by the various coordinates:

# In[29]:


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


# In[30]:


spher.coord_range()


# The specification of the hyperspherical coordinates is completed by providing the transition map to the stereographic chart $(A,(x,y,z))$:

# In[31]:


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()


# We also provide the inverse transition map by means of the method `set_inverse`, which (by default) performs a check that the provided formulas are indeed correct:

# In[32]:


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)


# The check is passed, modulo some lack of trigonometric simplifications in the first three lines.

# In[33]:


spher_to_stereoN.inverse().display()


# The transition map $(A,(\chi,\theta,\phi))\rightarrow (A,(x',y',z'))$ is obtained by combining the transition maps $(A,(\chi,\theta,\phi))\rightarrow (A,(x,y,z))$ and $(A,(x,y,z))\rightarrow (A,(x',y',z'))$:

# In[34]:


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


# Similarly, the transition map $(A,(x',y',z'))\rightarrow (A,(\chi,\theta,\phi))$ is obtained by combining the transition maps $(A,(x',y',z'))\rightarrow (A,(x,y,z))$ and $(A,(x,y,z))\rightarrow (A,(\chi,\theta,\phi))$:

# In[35]:


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


# At this stage, the user atlas of $\mathbb{S}^3$ is

# In[36]:


S3.atlas()


# Let us get the coordinate expression of the restriction of the embedding $\Phi$ to $A$:

# In[37]:


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


# In[38]:


Phi.display(spher, X4)


# In[39]:


Phi.display()


# ### Plots of the hyperspherical coordinate grid
# 
# First let us plot the chart $(A,(\chi,\theta,\phi))$ in terms of the stereographic chart $(U,(x,y,z))$ (notice that the "point at infinity" corresponds to $\chi\rightarrow \pi$, hence the $0.9$ truncation in the $\chi$ range)::

# In[40]:


graph = spher.plot(stereoN, 
                   number_values=7,
                   ranges={ch: (0, 0.9*pi)},
                   color={ch: 'green', th: 'blue', ph: 'red'},
                   label_axes=False)
show(graph, axes_labels=['x', 'y', 'z'])


# In terms of the stereographic coordinates $(V, (x',y',z'))$ (notice that the "point at infinity" corresponds then to $\chi\rightarrow 0$):

# In[41]:


graph = spher.plot(stereoS, 
                   number_values=7,
                   ranges={ch: (0.1, pi)},
                   color={ch: 'green', th: 'blue', ph: 'red'},
                   label_axes=False)
show(graph, axes_labels=["x'", "y'", "z'"])


# Of course we may use the embeddding $\Phi$ to get views of the hyperspherical coordinates in terms of the Cartesian coordinates $(T,X,Y,Z)$ of $\mathbb{R}^4$:

# In[42]:


graph = spher.plot(X4, mapping=Phi, 
                   ambient_coords=(X,Y,T),
                   number_values=7,
                   color={ch: 'green', th: 'blue', ph: 'red'},
                   label_axes=False)
show(graph, axes_labels=['X', 'Y', 'T'])


# In[43]:


graph = spher.plot(X4, mapping=Phi, 
                   ambient_coords=(X,Y,Z),
                   number_values=7,
                   color={ch: 'green', th: 'blue', ph: 'red'},
                   label_axes=False)
show(graph, axes_labels=['X', 'Y', 'Z'])


# ## Projection of $\mathbb{R}^4$ to $\mathbb{S}^3$

# We will need some projection operator from (a subset of) $\mathbb{R}^4$ to $\mathbb{S}^3$.
# Let $\mathbb{R}^4_N$ be $\mathbb{R}^4$ minus the hyperplane $T=-1$:

# In[44]:


R4N = R4.open_subset('R4N', latex_name=r'\mathbb{R}^4_N', 
                     coord_def={X4: T!=-1})
X4N = X4.restrict(R4N)


# In[45]:


ProjN = R4N.diff_map(U, {(X4N, stereoN): 
                         [X/(1+T), Y/(1+T), Z/(1+T)]},
                     name='P_N', latex_name=r'\Pi_N')
ProjN.display()


# Let us check that once applied to an embedded point of $U\subset \mathbb{S}^3$, this projection reduces to the identity:

# In[46]:


var('a b c', domain='real')
p = S3((a,b,c), chart=stereoN)
ProjN(Phi(p)) == p


# In[47]:


q = R4((sqrt(3)/2, sin(a)*cos(b)/2, sin(a)*sin(b)/2, cos(a)/2))
X4(q)


# In[48]:


all([q in R4N, Phi(ProjN(q)) == q])


# ## Quaternions

# We consider the (division) algebra of quaternions $\mathbb{H}$ as $\mathbb{R}^4$ endowed with the following (non-commutative) product:

# In[49]:


def qprod(p,q):
    if p in R4 and q in R4:
        T1, X1, Y1, Z1 = X4(p)
        T2, X2, Y2, Z2 = X4(q)
        return R4(((T1*T2-X1*X2-Y1*Y2-Z1*Z2).simplify_full(),
                   (T1*X2+X1*T2+Y1*Z2-Z1*Y2).simplify_full(),
                   (T1*Y2-X1*Z2+Y1*T2+Z1*X2).simplify_full(),
                   (T1*Z2+X1*Y2-Y1*X2+Z1*T2).simplify_full()))
    if p in S3 and q in S3:
        a = qprod(Phi(p),Phi(q))
        if X4(a) == (-1,0,0,0):
            return N
        return ProjN(R4N(a))
    raise ValueError("Cannot evaluate qprod of {} and {}".format(p,q))


# Note that we have extended the definition of the quaternionic product to $\mathbb{S}^3$ via the embedding $\Phi$. 
# 
# ### Distinguished quaternions on $\mathbb{S}^3$
# 
# Let us introduce two special points on $\mathbb{S}^3$: $\mathbf{1}$ and $-\mathbf{1}$.

# In[50]:


One = S3((0,0,0), chart=stereoN, name='1', latex_name=r'\mathbf{1}')
X4(Phi(One))


# As we can see from the Cartesian coordinates of $\Phi(\mathbf{1})$, the point $\mathbf{1}$ is actually nothing but the "South" pole used to define the stereographic chart $(V,(x',y',z'))$:

# In[51]:


One == S


# In[52]:


minusOne = S3((0,0,0), chart=stereoS, name='-1', latex_name=r'-\mathbf{1}')
X4(Phi(minusOne))


# The point $\mathbf{-1}$ is thus nothing but the "North" pole used to define the stereographic chart $(U,(x,y,z))$:

# In[53]:


minusOne == N


# Next we introduce points $\mathbf{i}$, $\mathbf{j}$ and $\mathbf{k}$ on $\mathbb{S}^3$:

# In[54]:


I = S3((1,0,0), chart=stereoN, name='i', latex_name=r'\mathbf{i}')
X4(Phi(I))


# In[55]:


stereoS(I)


# In[56]:


J = S3((0,1,0), chart=stereoN, name='j', latex_name=r'\mathbf{j}')
X4(Phi(J))


# In[57]:


stereoS(J)


# Since $\mathbf{j}$ lies in $A$, contrary to $\mathbf{i}$, we may ask for its hyperspherical coordinates:

# In[58]:


spher(J)


# In[59]:


K = S3((0,0,1), chart=stereoN, name='k', latex_name=r'\mathbf{k}')
X4(Phi(K))


# In[60]:


stereoS(K)


# Hamilton's fundamental relations
# $$ \mathbf{i} \mathbf{j} \mathbf{k} = \mathbf{-1} $$
# $$ \mathbf{i} \mathbf{j} = \mathbf{k},\quad \mathbf{j} \mathbf{k} = \mathbf{i}, \quad \mathbf{k} \mathbf{i} = \mathbf{j}$$
# are satisfied:

# In[61]:


qprod(I, qprod(J,K)) == minusOne


# In[62]:


all([qprod(I,J) == K, qprod(J,K) == I,
     qprod(K,I) == J])


# These relations imply $\mathbf{i}^2 = \mathbf{-1}$, $\mathbf{j}^2 = \mathbf{-1}$ and $\mathbf{k}^2 = \mathbf{-1}$:

# In[63]:


all([qprod(One,One) == One, qprod(I,I) == minusOne,
     qprod(J,J) == minusOne, qprod(K,K) == minusOne])


# Let us introduce $\mathbf{-i}$, $\mathbf{-j}$ and $\mathbf{-k}$, as points of $\mathbb{S}^3$:

# In[64]:


minusI = qprod(minusOne, I)
X4(Phi(minusI))


# In[65]:


minusJ = qprod(minusOne, J)
X4(Phi(minusJ))


# In[66]:


minusK = qprod(minusOne, K)
X4(Phi(minusK))


# ### Quaternionic conjugation

# In the comments below (but not in the SageMath code), we shall identify $\mathbf{1}\in \mathbb{S}^3$ with $\Phi(\mathbf{1})\in \mathbb{R}^4$, $\mathbf{i}\in \mathbb{S}^3$ with $\Phi(\mathbf{i})\in \mathbb{R}^4$, etc. In particular, we consider $(\mathbf{1}, \mathbf{i}, \mathbf{j},\mathbf{k})$ as a basis of the quaternion algebra $\mathbb{H}$. 
# 
# The *conjugate* of a quaternion $q = T + X\mathbf{i} + Y\mathbf{j} + Z\mathbf{k}$ is $\bar{q} = T - X\mathbf{i} - Y\mathbf{j} - Z\mathbf{k}$; hence we define:
# 

# In[67]:


def qconj(p):
    if p in R4:
        T, X, Y, Z = X4(p)
        return R4((T, -X, -Y, -Z))
    if p in S3:
        a = qconj(Phi(p))
        if X4(a) == (-1,0,0,0):
            return N
        return ProjN(a)
    raise ValueError("Cannot evaluate qconf of {}".format(p)) 


# In particular, we have $\bar{\mathbf{1}} = \mathbf{1}$, $\bar{\mathbf{i}} = -\mathbf{i}$,  $\bar{\mathbf{j}} = -\mathbf{j}$ and  $\bar{\mathbf{k}} = -\mathbf{k}$:

# In[68]:


all([qconj(One) == One, 
     qconj(I) == minusI,
     qconj(J) == minusJ, 
     qconj(K) == minusK])


# The conjugate of an element of $\mathbb{S}^3$

# In[69]:


assume(a != 0)  # to ensure that qconj(p) is not N
p = S3((a,b,c), chart=stereoN)
stereoN(qconj(p))


# In[70]:


p = S3((a,b,c), chart=stereoS)
stereoS(qconj(p))


# In[71]:


forget(a!=0)


# ### Norm of a quaternion
# 
# The quaternionic norm $\| q\| = \sqrt{q\bar{q}}$ coincide with the Euclidean norm in $\mathbb{R}^4$, so that $\mathbb{S}^3$ can be viewed as the set of unit quaternions; hence we define:

# In[72]:


def qnorm(p):
    if p in R4:
        T, X, Y, Z = X4(p)
        return (sqrt(T^2 + X^2 + Y^2 + Z^2)).simplify_full()
    if p in S3:
        return 1
    raise ValueError("Cannot evaluate qnorm of {}".format(p)) 


# In[73]:


var('d', domain='real')
q = R4((a,b,c,d))
qnorm(q)


# Let us check that $\| q\|^2 = q\bar{q}$:

# In[74]:


R4((qnorm(q)^2,0,0,0)) == qprod(q, qconj(q))


# As elements of $\mathbb{S}^3$, $\mathbf{1}$,  $\mathbf{i}$,  $\mathbf{j}$ and  $\mathbf{k}$ have all unit norm:

# In[75]:


(qnorm(One), qnorm(I), qnorm(J), qnorm(K)) == (1, 1, 1, 1)


# ## Hopf map

# We shall define the Hopf map by considering first the map 
# $$ \begin{array}{cccc}
#    C: &  \mathbb{R}^4 & \to & \mathbb{R}^4\\
#    &  p & \mapsto & p \mathbf{k} \bar{p}
#    \end{array} $$
# The coordinate expression of $C$ is obtained as follows:

# In[76]:


p = R4((T,X,Y,Z)) # a generic point of R^4
coord_Cp = X4(  qprod(p, qprod(Phi(K), qconj(p))) )
coord_Cp


# Therefore we define $C$ as

# In[77]:


C = R4.diff_map(R4, coord_Cp, name='C')
C.display()


# The restriction of $C$ to $\Phi(\mathbb{S}^3)\subset \mathbb{R}^4$ can be viewed as the map 
# $C\circ \Phi: \mathbb{S}^3 \to \mathbb{R}^4$ :

# In[78]:


CS = C * Phi
CS.display()


# On the above coordinate expressions, we note that the codomain of $C\circ \Phi$ lies in the hyperplane $T=0$, i.e. in the set $\operatorname{Im}\mathbb{H}$ of pure imaginary quaternions.
# Moreover, if we consider a generic point $p\in U\subset\mathbb{S}^3$:

# In[79]:


p = S3((a,b,c), chart=stereoN)


# we have $\| C\circ\Phi(p) \| = 1$:

# In[80]:


qnorm(CS(p))


# For the only point of $\mathbb{S}^3$ not lying in $U$, i.e. $N = -\mathbf{1}$, we have as well

# In[81]:


qnorm(CS(N))


# Hence the codomain of $C \circ\Phi$ lies in $\Phi(\mathbb{S}^3)$. From the previous result, we conclude that it actually lies in $\Phi(\mathbb{S}^3)\cap \operatorname{Im}\mathbb{H}$, which is a 2-sphere: the 2-sphere of unit imaginary quaternions.
# 
# In particular, we have:

# In[82]:


all([CS(K) == Phi(K), CS(One) == Phi(K), CS(minusOne) == Phi(K)])


# In[83]:


all([CS(I) == Phi(minusK), CS(J) == Phi(minusK),
     CS(minusI) == Phi(minusK), CS(minusJ) == Phi(minusK)])


# On $\Phi(\mathbb{S}^3)\cap \operatorname{Im}\mathbb{H}$, the inverse embedding $\Phi^{-1}$ coincides with 
# the projector $\Pi_N$ introduced above since $(\Phi(\mathbb{S}^3)\cap \operatorname{Im}\mathbb{H})\subset \mathbb{R}^4_N$. Hence the map $H = \Phi^{-1}\circ C \circ \Phi: \mathbb{S}^3 \to \mathbb{S}^3$ can be obtained as $\Pi_N\circ C \circ \Phi$:

# In[84]:


H = ProjN * CS.restrict(S3, subcodomain=R4N)


# Note that we have used the method `restrict` with the argument `subcodomain=R4N`
# to declare that the codomain of $C\circ\Phi$ actually lies in $\mathbb{R}^4_N$, so that the composition with $\Pi_N$ is well defined. 
# We have

# In[85]:


H.display()


# Actually since neither $N$ (which has $T=-1$) nor $S$ (which has $T=1$) lie in the codomain of $C\circ\Phi$, we may safely declare that the codomain of $H$ is $W = U\cap V$:

# In[86]:


H = H.restrict(S3, subcodomain=W)
H.display()


# We have $H(\mathbf{k})=H(\mathbf{1})=H(-\mathbf{k})=H(-\mathbf{1})=\mathbf{k}$:

# In[87]:


all([H(K) == K, H(One) == K, H(minusK) == K, 
     H(minusOne) == K])


# and $H(\mathbf{i})=H(\mathbf{j})=-\mathbf{k}$:

# In[88]:


all([H(I) == minusK, H(J) == minusK])


# Let us consider the expression of $H$ in stereographic coordinates:

# In[89]:


Hx, Hy, Hz = H.expr(stereoN, stereoN)
(Hx.factor(), Hy.factor(), Hz.factor())


# We have

# In[90]:


(Hx^2 + Hy^2 + Hz^2).simplify_full()


# which shows that the codomain of $H$ lies in the 2-sphere of equation $x^2+y^2+z^2=1$ in stereographic coordinates. This is not surprising since (i) the equation of the 2-sphere of
# unit imaginary quaternions, $\Phi(\mathbb{S}^3)\cap \operatorname{Im}\mathbb{H}$, is $T=0$ and $X^2+Y^2+Z^2=1$ and (ii) for 
# $T=0$, we have $x=X$, $y=Y$ and $z=Z$. 
# Let us construct this 2-sphere as a manifold by itself, which we call the **base 2-sphere**. This will be the codomain of the Hopf map. 

# ## The base 2-sphere (unit imaginary quaternions)

# In[91]:


S2 = Manifold(2, 'S^2', latex_name=r'\mathbb{S}^2')
print(S2)


# As for $\mathbb{S}^3$, we introduce on $\mathbb{S}^2$ two complementary stereographic coordinate systems:

# In[92]:


U2 = S2.open_subset('U_2')
V2 = S2.open_subset('V_2')
S2.declare_union(U2, V2)


# In[93]:


stereoN2.<x2,y2> = U2.chart("x2:x_2 y2:y_2")
stereoN2


# In[94]:


stereoS2.<xp2,yp2> = V2.chart(r"xp2:{x'}_2 yp2:{y'}_2")
stereoS2


# In[95]:


stereoN_to_S2 = stereoN2.transition_map(stereoS2, 
                                        (x2/(x2^2+y2^2), y2/(x2^2+y2^2)), 
                                        intersection_name='W_2',
                                        restrictions1= x2^2+y2^2!=0, 
                                        restrictions2= xp2^2+xp2^2!=0)
stereoN_to_S2.display()


# In[96]:


stereoS_to_N2 = stereoN_to_S2.inverse()
stereoS_to_N2.display()


# In[97]:


W2 = U2.intersection(V2)


# We embed the base 2-sphere $\mathbb{S}^2$ in $\mathbb{S}^3$ by considering that the North pole defining the above stereographic coordinates is $\mathbf{k}$, i.e. the point 
# $(x,y,z)=(0,0,1)$ in $\mathbb{S}^3$:

# In[98]:


Phi2 = S2.diff_map(S3, {(stereoN2, stereoN): 
                        [2*x2/(1+x2^2+y2^2), 
                         2*y2/(1+x2^2+y2^2),
                         (x2^2+y2^2-1)/(1+x2^2+y2^2)],
                        (stereoS2, stereoN): 
                        [2*xp2/(1+xp2^2+yp2^2), 
                         2*yp2/(1+xp2^2+yp2^2),
                         (1-xp2^2-yp2^2)/(1+xp2^2+yp2^2)]},
                   name='Phi2', latex_name=r'\Phi_2')
Phi2.display()


# The unit imaginary quaternions $\mathbf{i}$, $\mathbf{j}$, $\mathbf{k}$ and $-\mathbf{k}$ as elements of the base 2-sphere:

# In[99]:


I2 = S2((1,0), chart=stereoN2)
J2 = S2((0,1), chart=stereoN2)
K2 = S2((0,0), chart=stereoS2)
minusK2 = S2((0,0), chart=stereoN2)


# In[100]:


all([Phi2(I2) == I, Phi2(J2) == J,
     Phi2(K2) == K, Phi2(minusK2) == minusK])


# ### Spherical coordinates on $\mathbb{S}^2$
# 
# We introduce spherical coordinates $(\theta_2,\phi_2)$ on the base 2-sphere in the standard way (cf. the [2-sphere notebook](http://nbviewer.jupyter.org/github/sagemanifolds/SageManifolds/blob/master/Notebooks/SM_sphere_S2.ipynb)):

# In[101]:


A2 = W2.open_subset('A_2', 
                    coord_def={stereoN2.restrict(W2): (y2!=0, x2<0), 
                               stereoS2.restrict(W2): (yp2!=0, xp2<0)})
spher2.<th2,ph2> = A2.chart(r'th2:(0,pi):\theta_2 ph2:(0,2*pi):\phi_2')
spher2_to_stereoN2 = spher2.transition_map(stereoN2.restrict(A2), 
                                        (sin(th2)*cos(ph2)/(1-cos(th2)),
                                         sin(th2)*sin(ph2)/(1-cos(th2))))
spher2_to_stereoN2.set_inverse(2*atan(1/sqrt(x2^2+y2^2)), 
                               atan2(-y2,-x2)+pi)
spher2_to_stereoS2 = stereoN_to_S2.restrict(A2) * spher2_to_stereoN2
stereoS2_to_spher2 = spher2_to_stereoN2.inverse() * \
                     stereoN_to_S2.inverse().restrict(A2)


# In[102]:


A2.atlas()


# ### Projectors $\mathbb{S}^3 \to \mathbb{S}^2$
# 
# Let $W_{z\not=1}$ denote the subset of $W\subset\mathbb{S}^3$ defined by $z\not=1$:

# In[103]:


Wz1 = W.open_subset('Wz1', latex_name=r'W_{z\not=1}', 
                    coord_def={stereoN.restrict(W): z!=1})


# In[104]:


all([I in Wz1, J in Wz1, minusK in Wz1,
     K not in Wz1, One not in Wz1, minusOne not in Wz1])


# In[105]:


Proj2N = Wz1.diff_map(U2, {(stereoN.restrict(Wz1), stereoN2): 
                           [x/(1-z), y/(1-z)]},
                      name='P_2^N', latex_name=r'\Pi_2^N')
Proj2N.display()


# In[106]:


p = U2((a,b), chart=stereoN2)
Proj2N(Phi2(p)) == p


# In[107]:


assume(cos(a)!=1, cos(a)!=0)
p = U((sin(a)*cos(b), sin(a)*sin(b), cos(a)), chart=stereoN)
Phi2(Proj2N(p)) == p


# In[108]:


forget(cos(a)!=1, cos(a)!=0)


# Let $W_{z\not=-1}$ denote the subset of $W\subset\mathbb{S}^3$ defined by $z\not=-1$:

# In[109]:


Wzm1 = W.open_subset('Wzm1', latex_name=r'W_{z\not=-1}', 
                     coord_def={stereoN.restrict(W): z!=-1})


# In[110]:


all([I in Wzm1, J in Wzm1, K in Wzm1,
     minusK not in Wzm1, One not in Wzm1, 
     minusOne not in Wzm1])


# In[111]:


Proj2S = Wzm1.diff_map(V2, {(stereoN.restrict(Wzm1), stereoS2): 
                            [x/(1+z), y/(1+z)]},
                       name='P_2^S', latex_name=r'\Pi_2^S')
Proj2S.display()


# In[112]:


p = V2((a,b), chart=stereoS2)
Proj2S(Phi2(p)) == p


# In[113]:


assume(cos(a)!=-1, cos(a)!=0)
p = U((sin(a)*cos(b), sin(a)*sin(b), cos(a)), chart=stereoN)
Phi2(Proj2S(p)) == p


# In[114]:


forget(cos(a)!=-1, cos(a)!=0)


# ### Hopf map
# 
# We are now in position to define the Hopf map as a map $\mathbb{S}^3 \to \mathbb{S}^2$.
# To give its coordinate expressions, we have to consider that $\mathbb{S}^3$ is covered by
# two charts, `stereoN` = $(U,(x,y,z))$ and `stereoS` = $(V,(x',y',z'))$, and $\mathbb{S}^2$ is covered by two charts: 
# - `stereoN2` = $(U_2,(x_2,y_2))$, the domain of which contains 
#   all points of $\mathbb{S}^2$ but $\mathbf{k}$
# - `stereoS2` = $(V_2,(x'_2,y'_2))$, the domain of which contains 
#   all points of $\mathbb{S}^2$ but $-\mathbf{k}$.
# 
# First we search for all the points $p\in U$ such that $H(p)\in U_2$, i.e. such that
# $H(p)\not=\mathbf{k}$, or equivalently, $z(H(p))\not= 1$, or again $H(p)\in W_{z\not=1}$.
# On the chart $(U,(x,y,z))$, the expression of $z(H(p))$ is 

# In[115]:


Hx, Hy, Hz = H.expr(stereoN, stereoN)
Hz


# The condition $z(H(p))\not=1$ is

# In[116]:


Hz.numerator() - Hz.denominator() != 0


# which is equivalent to $x^2+y^2 \not= 0$. We define thus the subdomain
# $D_1 = U \setminus (H^{-1}(\mathbf{k})\cap U)$ as

# In[117]:


D1 = U.open_subset('D_1', coord_def=({stereoN: x^2+y^2!=0}))
stereoN_D1 = stereoN.restrict(D1)


# By construction, the restriction of $H$ to $D_1$ has $W_{z\not=1}$ as codomain and we declare the Hopf map on $D_1$ by considering the image points as points of $\mathbb{S}^2$ via $\Pi_2^N$:

# In[118]:


hD1 = Proj2N * H.restrict(D1, subcodomain=Wz1)
hD1.display()


# We have $\mathbf{i}\in D_1$ and $\mathbf{j}\in D_1$; we can check that 
# $h(\mathbf{i}) = h(\mathbf{j}) = -\mathbf{k}$:

# In[119]:


all([hD1(I) == minusK2, hD1(J) == minusK2])


# Let us now consider the points $p\in U$ such that $H(p)\in V_2$, i.e. such that
# $H(p)\not=-\mathbf{k}$, or equivalently, $z(H(p))\not= -1$, or again $H(p)\in W_{z\not=-1}$.
# The condition $z(H(p))\not= -1$ is equivalent to $s\not=0$ with

# In[120]:


s = ((Hz.numerator() + Hz.denominator())/2).simplify_full()
s


# Since

# In[121]:


(s-(x^2+y^2+z^2-1)^2).simplify_full()


# the condition $s\not =0$ is equivalent to 

# In[122]:


(x^2+y^2+z^2-1)^2 + 4*z^2 != 0


# i.e. to ($x^2+y^2\not= 1$ or $z\not= 0$). Hence we introduce the subset 
# $D_2 = U \setminus (H^{-1}(-\mathbf{k})\cap U)$ by

# In[123]:


D2 = U.open_subset('D_2', coord_def=({stereoN: (x^2+y^2!=1, z!=0)}))
stereoN_D2 = stereoN.restrict(D2)


# By construction, the restriction of $H$ to $D_2$ has $W_{z\not=-1}$ as codomain and we declare the Hopf map on $D_2$ by considering the image points as points of $\mathbb{S}^2$ via $\Pi_2^S$:

# In[124]:


hD2 = Proj2S * H.restrict(D2, subcodomain=Wzm1)
hD2.display()


# We have $\mathbf{k}\in D_2$, $-\mathbf{k}\in D_2$ and $\mathbf{1}\in D_2$; we can check that 
# $h(\mathbf{k}) = h(-\mathbf{k}) = h(\mathbf{1}) = \mathbf{k}$:

# In[125]:


all([hD2(K) == K2, hD2(minusK) == K2, hD2(One) == K2])


# Since $H^{-1}(\mathbf{k})\cap H^{-1}(-\mathbf{k})=\emptyset$, we have $U=D_1\cup D_2$:

# In[126]:


U.declare_union(D1, D2)


# Similarly let us consider the points $p\in V$ such that $H(p)\in U_2$, i.e. such that
# $H(p)\not=\mathbf{k}$, or equivalently, $z(H(p))\not= 1$, or again $H(p)\in W_{z\not=1}$.
# On the chart $(V,(x',y',z'))$, the expression of $z(H(p))$ is 

# In[127]:


Hx, Hy, Hz = H.expr(stereoS, stereoN)
Hz


# The condition $z(H(p))\not=1$ is

# In[128]:


Hz.numerator() - Hz.denominator() != 0


# which is equivalent to ${x'}^2+{y'}^2 \not= 0$. We define thus the subset
# $D_3 = V \setminus (H^{-1}(\mathbf{k})\cap V)$ as

# In[129]:


D3 = V.open_subset('D_3', coord_def=({stereoS: xp^2+yp^2!=0}))
stereoS_D3 = stereoS.restrict(D3)


# By construction, the restriction of $H$ to $D_2$ has $W_{z\not=1}$ as codomain and we declare the Hopf map on $D_3$ by considering the image points as points of $\mathbb{S}^2$ via $\Pi_2^N$:

# In[130]:


hD3 = Proj2N *  H.restrict(D3, subcodomain=Wz1)
hD3.display()


# We have $\mathbf{i}\in D_3$ and $\mathbf{j}\in D_3$; we can check that 
# $h(\mathbf{i}) = h(\mathbf{j}) = -\mathbf{k}$:

# In[131]:


all([hD3(I) == minusK2, hD3(J) == minusK2])


# Finally, let us consider the points $p\in V$ such that $H(p)\in V_2$, i.e. such that
# $H(p)\not=-\mathbf{k}$, or equivalently, $z(H(p))\not= -1$, or again $H(p)\in W_{z\not=-1}$.
# The condition $z(H(p))\not= -1$ is equivalent to $s\not=0$ with

# In[132]:


s = ((Hz.numerator() + Hz.denominator())/2).simplify_full()
s


# Since 

# In[133]:


(s-(xp^2+yp^2+zp^2-1)^2).simplify_full()


# the condition $s\not=0$ is equivalent to

# In[134]:


(xp^2+yp^2+zp^2-1)^2 + 4*zp^2 == 0


# i.e. to (${x'}^2+{y'}^2\not= 1$ or $z'\not= 0$). Hence we introduce the subset 
# $D_4 = V \setminus (H^{-1}(-\mathbf{k})\cap V)$ by

# In[135]:


D4 = V.open_subset('D_4', coord_def=({stereoS: (xp^2+yp^2!=1, zp!=0)}))
stereoS_D4 = stereoS.restrict(D4)


# By construction, the restriction of $H$ to $D_4$ has $W_{z\not=-1}$ as codomain and we declare the Hopf map on $D_4$ by considering the image points as points of $\mathbb{S}^2$ via $\Pi_2^S$:

# In[136]:


hD4 = Proj2S * H.restrict(D4, subcodomain=Wzm1)
hD4.display()


# We have $-\mathbf{1}\in D_4$ and we can check that 
# $h(-\mathbf{1}) = \mathbf{k}$:

# In[137]:


hD4(minusOne) == K2


# Since $H^{-1}(\mathbf{k})\cap H^{-1}(-\mathbf{k})=\emptyset$, we have $V=D_3\cup D_4$:

# In[138]:


V.declare_union(D3, D4)


# #### Declaration of the Hopf map
# 
# We construct the Hopf map $h:\mathbb{S}^3\to \mathbb{S}^2$ from the coordinate expressions of its restriction to $D_1$, $D_2$, $D_3$ and $D_4$, as obtained above:

# In[139]:


h = S3.diff_map(S2, name='h')
h.add_expression(stereoN_D1, stereoN2, hD1.expr(stereoN_D1, stereoN2))
h.add_expression(stereoN_D2, stereoS2, hD2.expr(stereoN_D2, stereoS2))
h.add_expression(stereoS_D3, stereoN2, hD3.expr(stereoS_D3, stereoN2))
h.add_expression(stereoS_D4, stereoS2, hD4.expr(stereoS_D4, stereoS2))


# In[140]:


h.display(stereoN_D1, stereoN2)


# In[141]:


h.display(stereoN_D2, stereoS2)


# In[142]:


h.display(stereoS_D3, stereoN2)


# In[143]:


h.display(stereoS_D4, stereoS2)


# Let us check that $h(\mathbf{1})=h(-\mathbf{1})=h(\mathbf{k})=h(-\mathbf{k})=\mathbf{k}$
# and $h(\mathbf{i})=h(-\mathbf{i})=h(\mathbf{j})=h(-\mathbf{j})=-\mathbf{k}$:

# In[144]:


all([h(One)==K2, h(minusOne)==K2, 
     h(K)==K2, h(minusK)==K2,
     h(I)==minusK2, h(minusI)==minusK2,
     h(J)==minusK2, h(minusJ)==minusK2])


# ### Expression of the Hopf map in spherical coordinates

# In[145]:


D1A = A.intersection(D1)
spherD1A = spher.restrict(D1A)
spherD1A


# In[146]:


stereoND1A = stereoN_D1.restrict(D1A)
stereoND1A.add_restrictions((y!=0, x<0))


# In[147]:


D1A.atlas()


# In[148]:


spher_to_stereoND1A = spher_to_stereoN.restrict(D1A)
stereoN_to_spherD1A = spher_to_stereoN.inverse().restrict(D1A)
spher_to_stereoSD1A = spher_to_stereoS.restrict(D1A)
stereoS_to_spherD1A = stereoS_to_spher.restrict(D1A)


# In[149]:


h.expr(stereoND1A, stereoN2)  # necessary
h.display(spherD1A, stereoN2)


# In[150]:


hA = h.restrict(D1A, subcodomain=A2)
hA.display(spherD1A, spher2)


# ## Hopf coordinates
# 
# The Hopf coordinates are coordinates $(\eta,\alpha,\beta)$ on $\mathbb{S}^3$ which are related to the Cartesian coordinates on $\mathbb{R}^4$ (via the embedding $\Phi$) by
# $$
# \begin{equation} \tag{1}
#     \left\{ \begin{array}{lcl}
#         T & = &\cos\eta \sin\alpha \\
#         X & = &\sin\eta \cos(\alpha+\beta) \\        
#         Y & = &\sin\eta \sin(\alpha+\beta) \\
#         Z & = &\cos\eta \cos\alpha 
#       \end{array} \right . 
# \end{equation}
# $$
# and whose range is $\eta\in(0,\pi/2)$, $\alpha\in (0, 2\pi)$, $\beta\in (0, 2\pi)$. They are
# defined in $D_1$ minus the points for which $X^2+Y^2+T^2=1$ (limit $\alpha\rightarrow 0$
# or $2\pi$) or $TX-YZ=0$ (limit $\beta\rightarrow 0$ or $2\pi$). In terms of the stereograĥic coordinates, this corresponds to $x^2+y^2+z^2=1$ or to $x(1-x^2-y^2-z^2)-2yz=0$. Hence we declare the domain $B$ of Hopf coordinates as 

# In[151]:


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


# The limiting surface $x(1-x^2-y^2-z^2)-2yz=0$, where $\beta\rightarrow 0$ or $2\pi$:

# In[152]:


beta_zero = implicit_plot3d(x*(1-x^2-y^2-z^2)-2*y*z==0, 
                            (x,-3,3), (y,-3,3), (z,-3,3))
show(beta_zero)


# We define the Hopf coordinates and provide the transition map to the stereographic coordinates $(x,y,z)$:

# In[153]:


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


# In[154]:


Hcoord_to_stereoN = Hcoord.transition_map(
                        stereoN.restrict(B),
                        (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[155]:


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))


# Note that the test of the inverse coordinate transformation is passed, modulo some lack of trigonometric simplifications.

# In[156]:


Hcoord_to_stereoN.inverse().display()


# ### Embedding $\Phi$ in terms of Hopf coordinates

# In[157]:


PhiB = Phi.restrict(B)
PhiB.display()


# The expression of $\Phi$ in terms of the Hopf coordinates can be simplified by means of the method `trig_reduce`:

# In[158]:


PhiB.expression(Hcoord, X4)[1].factor().trig_reduce()


# In[159]:


PhiB.expression(Hcoord, X4)[2].factor().trig_reduce()


# Hence the recover the expression (1) above.

# ### Expression of the Hopf map in terms of Hopf coordinates
# 
# The expression of $h$ in terms of stereographic coordinates on $B$ is

# In[160]:


h.display(stereoN.restrict(B), stereoN2)


# Let us ask for the expression in terms of Hoopf coordinates:

# In[161]:


h.display(Hcoord, stereoN2)


# We notice that the image point in $\mathbb{S}^2$ is independent of the value of $\alpha$.  
# 
# The expression of $h$ is even simpler in terms of spherical coordinates on $\mathbb{S}^2$:

# In[162]:


hB = h.restrict(B, subcodomain=A2)
hB.display(Hcoord, spher2)


# We are facing some lack of simplification: 
# 
# - $\operatorname{atan}(\sin\eta/\cos\eta)$ should simplify to $\eta$
# - the `atan2` function should simplify to $\operatorname{atan2}(\sin\beta,-\cos\beta) = \pi-\beta$, so that $\phi_2=\beta$
#   
# Hence the right-hand side of the above formula simplifies to $(\theta_2,\phi_2)=(2\eta,\beta)$. We enforce it by the method `add_expression`:

# In[163]:


hB.add_expression(Hcoord, spher2, (2*eta, bet))
hB.display(Hcoord, spher2)


# The expression of the Hopf map in terms of the Hopf coordinates is thus very simple, which justifies the name given to these coordinates.
# We also recover a very simple expression when asking to express $C\circ \Phi$ (from which $h$ has been constructed) in terms of the Hopf coordinates: 

# In[164]:


CS.restrict(B).display(Hcoord, X4)


# The right-hand side should simplify to
# $$(T,X,Y,Z) = (0, \sin(2\eta)\cos\beta, \sin(2\eta)\sin\beta, \cos(2\eta))$$
# We can (partially) obtain this by means of `trig_reduce()`:

# In[165]:


for cp in CS.restrict(B).expr(Hcoord, X4):
    show(cp.trig_reduce().simplify_full())


# ## The Hopf map as a fibration of $\mathbb{S}^3$
# 
# The above results show that the image by $h$ of a point of Hopf coordinates $(\eta,\alpha,\beta)$ is independent of $\alpha$. Since $h$ is surjective, this means that for any point 
# $p\in\mathbb{S}^2$, the preimage set $h^{-1}(p)$ corresponds to a fixed value of $(\eta,\beta)$, with $\alpha$ taking all values in the range $(0,2\pi)$. From Eq. (1), the projection of  $h^{-1}(p)$ in the $(T,X)$-plane is a circle of radius $\cos\eta$ centered on $(0,0)$, while its projection in the $(X,Y)$-plane is a circle of radius $\sin\eta$ centered on $(0,0)$. We conclude that $h^{-1}(p)$ is a great circle of $\mathbb{S}^3$, sometimes called a *Hopf circle*. 
# 
# It follows that $\mathbb{S}^3$ has the structure of a **fiber bundle** over $\mathbb{S}^2$ with 
# $\mathbb{S}^1$ fibers. The Hopf map $h:\mathbb{S}^3\to\mathbb{S}^2$ is then nothing but the projection map of this bundle.

# We can get a first view of the fibers by plotting the Hopf coordinates in terms of the stereographic ones, for a fixed value of $\eta$ ($\eta=\pi/4$): the lines along which $\alpha$ varies while $\beta$ is kept fixed, hence the fibers $h^{-1}(p)$, are plotted in green. They are
# indeed circles (remember that stereographic projection preserves circles):

# In[166]:


graph = Hcoord.plot(stereoN, 
                    fixed_coords={eta: pi/4},
                    color={alp: 'green', bet: 'orange'},
                    number_values=9, label_axes=False)
show(graph, axes_labels=['x','y','z'])


# Note that all the green circles are linked.
# The same plot, but in terms of the coordinates $(X,Y,T)$ of $\mathbb{R}^4$ via the embedding $\Phi$:

# In[167]:


graph = Hcoord.plot(X4, mapping=PhiB, ambient_coords=(X,Y,T), 
                    fixed_coords={eta: pi/4},
                    color={alp: 'green', bet: 'orange'},
                    number_values=9, label_axes=False)
show(graph, axes_labels=['X','Y','T'])


# We may fix $\beta$ instead of $\eta$, in plotting the Hopf coordinates in terms of the stereographic ones. The fibers are still the green circles, the red lines being lines along which $\eta$ varies at fixed $\alpha$. Again, note that all the green circles are linked.

# In[168]:


graph = Hcoord.plot(stereoN, 
                    fixed_coords={bet: pi/2}, 
                    ranges={eta: (0.25, 1.5)},
                    color={eta: 'red', alp: 'green'},
                    number_values=9, plot_points=150, 
                    label_axes=False)
show(graph, axes_labels=['x','y','z'])


# If we vary the three coordinates $(\eta,\alpha,\beta)$, we get the following plots, where $\mathbb{S}^3$ appears as filled by the grid of Hoopf coordinates:

# In[169]:


graph = Hcoord.plot(X4, mapping=PhiB, ambient_coords=(X,Y,Z), 
                    color={eta: 'red', alp: 'green', bet: 'orange'},
                    number_values=7, label_axes=False)
show(graph, axes_labels=['X','Y','Z'])


# In[170]:


#graph = Hcoord.plot(X4, mapping=PhiB, ambient_coords=(X,Y,T), 
#                    color={eta: 'red', alp: 'green', bet: 'orange'},
#                    number_values=7, label_axes=False)
# show(graph, axes_labels=['X','Y','T'])


# In[171]:


#graph = Hcoord.plot(X4, mapping=PhiB, ambient_coords=(X,Z,T), 
#                    color={eta: 'red', alp: 'green', bet: 'orange'},
#                    number_values=7, label_axes=False)
#show(graph, axes_labels=['X','Z','T'])


# ## Fibers of the Hopf fibration
# 
# For a point $p\in\mathbb{S}^2$, identified by its spherical coordinates $(\theta_2,\phi_2)=(2\eta,\beta)$, the fiber $h^{-1}(p)$ can be seen as a curve in $\mathbb{S}^3$:

# In[172]:


R.<t> = RealLine()
def fiber(eta, bet):
    return S3.curve({Hcoord: (eta, t, bet)}, (t, 0, 2*pi))


# For instance, the fiber above the point $(\theta_2,\phi_2)=(\pi/3,\pi/4)$ is 

# In[173]:


F = fiber(pi/6, pi/4)
F


# In[174]:


F.display()


# In[175]:


graph_F = F.plot(chart=stereoN, color='green', plot_points=100)
graph_F


# ### Plot of the fibers
# 
# Let us plot the fibers for three values of $\eta$: $\eta=\pi/6$ (turquoise), $\eta=\pi/4$ (gold)
# and $\eta=5\pi/12$ (red). For each value of $\eta$, we note that the fibers fill a torus.

# In[176]:


graph= Graphics()
etas = {pi/6: ['turquoise', (0, 2*pi), 30], 
        pi/4: ['gold', (0,2*pi), 30],
        5*pi/12: ['red', (0, 2*pi), 30]}
for eta_v, param in etas.items():
    color = param[0]
    beta_min, beta_max = param[1]
    nb = param[2]
    db = (beta_max - beta_min)/(nb-1)
    betas = [beta_min + db*k for k in range(nb)]
    for beta_v in betas:
        F = fiber(eta_v, beta_v)
        F.coord_expr(stereoN.restrict(B))
        graph += F.plot(chart=stereoN, color=color, plot_points=150,
                        label_axes=None)
graph


# In[177]:


#show(graph, viewer='tachyon', figsize=24, aspect_ratio=1, frame=False)


# A "top" view of the fibers, obtained by projection to the $(x,y)$-plane (note that for clarity, we have reduced the number of fibers from 30 to 12):

# In[178]:


graph= Graphics()
etas = {pi/6: ['turquoise', (0, 2*pi), 12], 
        pi/4: ['gold', (0,2*pi), 12],
        5*pi/12: ['red', (0, 2*pi), 12]}
for eta_v, param in etas.items():
    color = param[0]
    beta_min, beta_max = param[1]
    nb = param[2]
    db = (beta_max - beta_min)/(nb-1)
    betas = [beta_min + db*k for k in range(nb)]
    for beta_v in betas:
        F = fiber(eta_v, beta_v)
        F.coord_expr(stereoN.restrict(B))
        graph += F.plot(chart=stereoN, ambient_coords=(x,y), 
                        color=color, plot_points=150)
show(graph, aspect_ratio=1)


# The same fibers, but viewed in $\mathbb{R}^4$ (via the embedding $\Phi$), in terms of the coordinates $(T,X,Z)$:

# In[179]:


graph = Graphics()
for eta_v, param in etas.items():
    color = param[0]
    beta_min, beta_max = param[1]
    nb = param[2]
    db = (beta_max - beta_min)/(nb-1)
    betas = [beta_min + db*k for k in range(nb)]
    for beta_v in betas:
        F = fiber(eta_v, beta_v)
        F.coord_expr(stereoN.restrict(B))
        graph += F.plot(chart=X4, ambient_coords=(X,Z,T), 
                        mapping=Phi, color=color, 
                        plot_points=200, label_axes=None)
show(graph, axes_labels=['X', 'Z', 'T'])


# or in terms of the coordinates $(X,Y,Z)$:

# In[180]:


graph = Graphics()
for eta_v, param in etas.items():
    color = param[0]
    beta_min, beta_max = param[1]
    nb = param[2]
    db = (beta_max - beta_min)/(nb-1)
    betas = [beta_min + db*k for k in range(nb)]
    for beta_v in betas:
        F = fiber(eta_v, beta_v)
        F.coord_expr(stereoN.restrict(B))
        graph += F.plot(chart=X4, ambient_coords=(X,Y,Z), 
                        mapping=Phi, color=color, 
                        plot_points=100, label_axes=None)
show(graph, axes_labels=['X', 'Y', 'Z'])


# In this notebook, we have used 18 charts on $\mathbb{S}^3$:

# In[181]:


S3.atlas()


# In[182]:


len(S3.atlas())


# corresponding actually to 4 different coordinate systems:

# In[183]:


S3.top_charts()


# In[184]:


print("Total elapsed time: {} s".format(time.perf_counter() - comput_time0))


# For a follow-up, see the [notebook devoted to vector fields](http://nbviewer.jupyter.org/github/sagemanifolds/SageManifolds/blob/master/Notebooks/SM_sphere_S3_vectors.ipynb) on $\mathbb{S}^3$.
# 
# You may also visit [Niles Johnson's page](http://nilesjohnson.net/hopf.html) for a very nice 
# animated 3D visualization of the Hopf fibration, also constructed with SageMath.