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

# # Sphere $\mathbb{S}^2$
# 
# This notebook demonstrates some differential geometry capabilities of SageMath on the example of the 2-dimensional sphere. The corresponding tools have been developed within
# the [SageManifolds](https://sagemanifolds.obspm.fr) project.

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

# In[1]:


version()


# First we set up the notebook to display math formulas using LaTeX formatting:

# In[2]:


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


# ## $\mathbb{S}^2$ from the manifold catalog

# The 2-sphere, with predefined charts and embedding in the Euclidean 3-space, can be obtained directly from SageMath's manifold catalog:

# In[3]:


S2 = manifolds.Sphere(2)
S2


# In[4]:


print(S2)


# In[5]:


S2.spherical_coordinates()


# In[6]:


S2.metric().display()


# ## $\mathbb{S}^2$ defined from scratch as a 2-dimensional smooth manifold

# For the purpose of introducing generic smooth manifolds in SageMath, we shall not use the above predefined object. Instead we shall construct $\mathbb{S}^2$ from scratch, by invoking the generic function `Manifold`:

# In[7]:


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


# The first argument, `2`, is the dimension of the manifold, while the second argument is the symbol used to label the manifold.
# 
# 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\}$; the default value is `start_index=0` and yields $\{0,1\}$.

# The function `Manifold` has actually many options, which are displayed via the command `Manifold?`:

# In[8]:


# Manifold?


# By default `Manifold` constructs a smooth manifold over $\mathbb{R}$:

# In[9]:


print(S2)


# In[10]:


S2


# $\mathbb{S}^2$ is in the category of smooth manifolds over $\mathbb{R}$:

# In[11]:


S2.category()


# In[12]:


print(S2.category())


# At the moment, the real field $\mathbb{R}$ is modeled by 53-bit floating-point approximations, but this plays no role in the manifold implementation:

# In[13]:


print(S2.base_field())


# In[14]:


S2.base_field() is RR


# ### Coordinate charts on $\mathbb{S}^2$
# 
# The function `Manifold` generates a manifold with no-predefined coordinate chart, so that the manifold (user) **atlas** is empty:

# In[15]:


S2.atlas()


# Let us introduce some charts. At least two charts are necessary to cover the sphere. Let us choose the charts associated with the **stereographic projections** to the equatorial plane from the North pole and the South pole respectively. We first introduce the open subsets covered by these two charts: 
# $$ U := \mathbb{S}^2\setminus\{N\}, $$  
# $$ V := \mathbb{S}^2\setminus\{S\}, $$
# where $N$ is a point of $\mathbb{S}^2$, which we shall call the **North pole**, and $S$ is the point of $U$ of stereographic coordinates $(0,0)$, which we call the **South pole**:

# To find the method to create an open subset, we type `U = S2.<TAB>` to get the list of possible methods by autocompletion:

# In[16]:


#U = S2.


# In[17]:


U = S2.open_subset('U')
print(U)


# In[18]:


V = S2.open_subset('V')
print(V)


# As an open subset of a smooth manifold, $U$ is itself a smooth manifold:

# In[19]:


print(U.category())


# <p>We declare that $\mathbb{S}^2 = U \cup V$:</p>

# In[20]:


S2.declare_union(U, V)


# The **stereographic chart** on $U$ is constructed from the stereographic projection from the North pole onto the equatorial plane: in the [Wikipedia figure](https://en.wikipedia.org/wiki/Stereographic_projection) below, the stereographic coordinates $(x,y)$ of the point $P\in U$ are the Cartesian coordinates of the point $P'$ in the equatorial plane.
# 
# ![stereographic projection](https://upload.wikimedia.org/wikipedia/commons/thumb/e/e3/Stereoprojzero.svg/241px-Stereoprojzero.svg.png)
# 
# We call this chart `stereoN` and construct it via the method `chart`:

# In[21]:


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


# The syntax `.<x,y>` in the left-hand side implies that the Python names `x` and `y` are added to the global namespace, to access to the two coordinates of the chart as symbolic variables. This allows one to refer subsequently to the coordinates by these names. Besides, in the present case, the function `chart()` has no argument, which implies that the coordinate symbols will be `x` and `y` (i.e. exactly the characters appearing in the `<...>` operator) and that each coordinate range is $(-\infty,+\infty)$. As we will see below, for other cases, an argument must be passed to `chart()`  to specify each coordinate symbol and range, as well as some specific LaTeX symbol.

# *Note:* the notation `.<x,y>` is not standard Python syntax, but a "SageMath enhanced" syntax. 
# Actually the SageMath kernel preparses the cell entries before sending them to the Python interpreter. The outcome of the preparser is shown by the function `preparse`. In the present case:

# In[22]:


print(preparse("stereoN.<x,y> = U.chart()"))


# Another example of preparsing:

# In[23]:


preparse("2^3")


# The chart created by the above command:

# In[24]:


stereoN


# In[25]:


print(stereoN)


# In[26]:


stereoN.coord_range()


# The coordinates can be accessed individually, either by means of their indices in the chart ( following the convention `start_index=1` set in the manifold's definition) or by their names as Python variables:

# In[27]:


stereoN[1]


# In[28]:


y is stereoN[2]


# The coordinates are SageMath symbolic expressions:

# In[29]:


type(y)


# In[30]:


y.parent()


# #### Stereographic coordinates from the South Pole
# 
# We introduce on $V$ the coordinates $(x',y')$ corresponding to the stereographic projection from the South pole:

# In[31]:


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


# In this case, the string argument passed to `chart` stipulates that the text-only names of the coordinates are xp and yp (same as the Python variables names defined within the `<...>` operator in the left-hand side), while their LaTeX names are $x'$ and $y'$.

# In[32]:


stereoS


# At this stage, the user's atlas on the manifold is made of two charts:

# In[33]:


S2.atlas()


# To complete the construction of the manifold structure, we have 
# to specify the <strong>transition map</strong> between the charts `stereoN` = $(U,(x,y))$ and `stereoS` = $(V,(x',y'))$; it is given by standard inversion formulas:

# In[34]:


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


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

# In[35]:


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


# <p>In the present case, the situation is of course perfectly symmetric regarding the coordinates $(x,y)$ and $(x',y')$.</p>
# <p>At this stage, the user's atlas has four charts:</p>

# In[36]:


S2.atlas()


# <p>Let us store $W = U\cap V$ into a Python variable for future use:</p>

# In[37]:


W = U.intersection(V)


# <p>Similarly we store the charts $(W,(x,y))$ (the restriction of  $(U,(x,y))$ to $W$) and $(W,(x',y'))$ (the restriction of $(V,(x',y'))$ to $W$) into Python variables:</p>

# In[38]:


stereoN_W = stereoN.restrict(W)
stereoN_W


# In[39]:


stereoN_W is S2.atlas()[2]


# In[40]:


stereoS_W = stereoS.restrict(W)
stereoS_W


# Coordinate charts are endoved with a method `plot`. For instance, 
# we may plot the chart $(W, (x',y'))$ in terms of itself, as a grid:

# In[41]:


stereoS_W.plot()


# More interestingly, let us plot the stereographic chart $(x',y')$ in terms of the stereographic chart $(x,y)$ on the domain $W$ where both systems overlap. We split the plot in four parts to avoid the singularity at $(x',y')=(0,0)$ and
# ask for the coordinate lines along which $x'$ (resp. $y'$) varies to be colored in purple (resp. cyan):

# In[42]:


graph = (stereoS_W.plot(stereoN, ranges={xp:[-6,-0.02], yp:[-6,-0.02]},
                          color={xp: 'purple', yp: 'cyan'}) 
         + stereoS_W.plot(stereoN, ranges={xp:[-6,-0.02], yp:[0.02,6]},
                          color={xp: 'purple', yp: 'cyan'})
         + stereoS_W.plot(stereoN, ranges={xp:[0.02,6], yp:[-6,-0.02]},
                          color={xp: 'purple', yp: 'cyan'})
         + stereoS_W.plot(stereoN, ranges={xp:[0.02,6], yp:[0.02,6]},
                          color={xp: 'purple', yp: 'cyan'}))
graph.show(xmin=-1.5, xmax=1.5, ymin=-1.5, ymax=1.5)


# ### Spherical coordinates
# 
# The standard **spherical coordinates** $(\theta,\phi)$ are defined on the open domain $A\subset W \subset \mathbb{S}^2$ that is the complement of the "origin meridian"; since the latter is the half-circle defined by $y=0$ and $x\geq 0$, we declare:

# In[43]:


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


# <p>The restriction of the stereographic chart from the North pole to $A$ is</p>

# In[44]:


stereoN_A = stereoN_W.restrict(A)
stereoN_A


# <p>We then declare the chart $(A,(\theta,\phi))$ by specifying the intervals $(0,\pi)$ and $(0,2\pi)$ spanned by respectively $\theta$ and $\phi$:</p>

# In[45]:


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


# In[46]:


spher.coord_range()


# <p>The specification of the spherical coordinates is completed by providing the transition map with the stereographic chart $(A,(x,y))$:</p>

# In[47]:


spher_to_stereoN = spher.transition_map(stereoN_A, 
                                        (sin(th)*cos(ph)/(1-cos(th)),
                                         sin(th)*sin(ph)/(1-cos(th))))
spher_to_stereoN.display()


# We also provide the inverse transition map:

# In[48]:


spher_to_stereoN.set_inverse(2*atan(1/sqrt(x^2+y^2)), atan2(-y,-x)+pi)


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

# In[49]:


spher_to_stereoN.inverse().display()


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

# In[50]:


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


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

# In[51]:


stereoS_to_N_A = stereoN_to_S.inverse().restrict(A)
stereoS_to_spher = spher_to_stereoN.inverse() * stereoS_to_N_A 
stereoS_to_spher.display()


# <p>The user atlas of $\mathbb{S}^2$ is now</p>

# In[52]:


S2.atlas()


# <p>Let us draw the grid of spherical coordinates $(\theta,\phi)$ in terms of stereographic coordinates from the North pole $(x,y)$:</p>

# In[53]:


spher.plot(stereoN, number_values=15, ranges={th: (pi/8,pi)})


# ## Points on $\mathbb{S}^2$
# 
# To create a point on $\mathbb{S}^2$, we use SageMath's ***parent / element*** syntax, i.e. the call operator `S2(...)` acting on the parent `S2`, with the point's coordinates in some chart as argument. 
# 
# For instance, we declare the **North pole** (resp. the **South pole**) as the point of coordinates $(0,0)$ in the chart $(V,(x',y'))$ (resp. in the chart $(U,(x,y))$):

# In[54]:


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


# In[55]:


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


# In[56]:


N.parent()


# In[57]:


S.parent()


# <p>We have of course</p>

# In[58]:


N in S2


# In[59]:


N in U


# In[60]:


N in V


# In[61]:


N in A


# Let us introduce some point $p$ of stereographic coordinates $(x,y) = (1,2)$:

# In[62]:


p = S2((1,2), chart=stereoN, name='p')


# $p$ lies in the open subset $A$:

# In[63]:


p in A


# #### Charts acting on points

# By definition, a chart maps points to pairs of real numbers (the point's coordinates): 

# In[64]:


stereoN(p)  # by definition of p


# In[65]:


stereoS(p)


# In[66]:


spher(p)


# In[67]:


stereoS(N)


# In[68]:


#stereoN(N)   ## returns an error


# ## Maps between manifolds: the embedding of $\mathbb{S}^2$ into $\mathbb{R}^3$
# 
# Let us first declare $\mathbb{R}^3$ as the 3-dimensional Euclidean space, denoting the Cartesian coordinates by
# $(X,Y,Z)$:

# In[69]:


R3.<X,Y,Z> = EuclideanSpace(name='R^3', latex_name=r'\mathbb{R}^3', metric_name='h')
cartesian = R3.cartesian_coordinates()
cartesian


# As an Euclidean space, `R3` is considered by Sage as a smooth manifold:

# In[70]:


print(R3.category())


# The embedding $\Phi: \mathbb{S}^2 \longmapsto \mathbb{R}^3$ is then defined via the method `diff_map` by providing the standard formulas relating the stereographic coordinates to the ambient Cartesian ones when considering the **stereographic projection** from the point $(0,0,1)$ (North pole) or $(0, 0, -1)$ (South pole) to the equatorial plane $Z=0$:

# In[71]:


Phi = S2.diff_map(R3, {(stereoN, cartesian): 
                       [2*x/(1+x^2+y^2), 2*y/(1+x^2+y^2),
                        (x^2+y^2-1)/(1+x^2+y^2)],
                       (stereoS, cartesian): 
                       [2*xp/(1+xp^2+yp^2), 2*yp/(1+xp^2+yp^2),
                        (1-xp^2-yp^2)/(1+xp^2+yp^2)]},
                  name='Phi', latex_name=r'\Phi')


# In[72]:


Phi.display()


# In[73]:


Phi.parent()


# In[74]:


print(Phi.parent())


# In[75]:


Phi.parent() is Hom(S2, R3)


# <p>$\Phi$ maps points of $\mathbb{S}^2$ to points of $\mathbb{R}^3$:</p>

# In[76]:


N1 = Phi(N)
print(N1)
N1


# In[77]:


cartesian(N1)


# In[78]:


S1 = Phi(S)
print(S1)
S1


# In[79]:


cartesian(S1)


# In[80]:


p1 = Phi(p)
print(p1)
p1


# In[81]:


cartesian(p1)


# $\Phi$ has been defined in terms of the stereographic charts $(U,(x,y))$ and $(V,(x',y'))$, but we may ask its expression in terms of spherical coordinates. This triggers a computation involving the transition map $(A,(x,y))\rightarrow (A,(\theta,\phi))$:

# In[82]:


Phi.display(stereoN_A, cartesian)


# In[83]:


Phi.display(spher, cartesian)


# Let us use $\Phi$ to draw the grid of spherical coordinates $(\theta,\phi)$ in terms of the Cartesian coordinates $(X,Y,Z)$ of $\mathbb{R}^3$:

# In[84]:


graph_spher = spher.plot(chart=cartesian, mapping=Phi, number_values=11, 
                         color='blue', label_axes=False)
graph_spher


# <p>We may also use the embedding $\Phi$ to display the stereographic coordinate grid in terms of the Cartesian coordinates in $\mathbb{R}^3$. First for the stereographic coordinates from the North pole:</p>

# In[85]:


graph = stereoN.plot(chart=cartesian, mapping=Phi, number_values=25, 
                     label_axes=False)
graph


# <p>and then have a view with the stereographic coordinates from the South pole superposed (in green):</p>

# In[86]:


graph += stereoS.plot(chart=cartesian, mapping=Phi, number_values=25, 
                      color='green', label_axes=False)
graph


# We may add the points $N$, $S$ and $p$ to the graphic, thanks to the method `plot` of points:

# In[87]:


graph += N.plot(chart=cartesian, mapping=Phi, color='red', 
                label_offset=0.05)
graph += S.plot(chart=cartesian, mapping=Phi, color='green', 
                label_offset=0.05)
graph += p.plot(chart=cartesian, mapping=Phi, color='blue', 
                label_offset=0.05)
graph


# ## Tangent spaces
# 
# The **tangent space** to the manifold $\mathbb{S}^2$ at the point $p$ is

# In[88]:


Tp = S2.tangent_space(p)
print(Tp)
Tp


# $T_p \mathbb{S}^2$ is a vector space over $\mathbb{R}$ (represented here by Sage's symbolic ring SR):

# In[89]:


print(Tp.category())


# <p>Its dimension equals the manifold's dimension:</p>

# In[90]:


dim(Tp)


# In[91]:


dim(Tp) == dim(S2)


# The tangent space at $p$ is the **fiber over** $p$ of the **tangent bundle** $T\mathbb{S}^2$:

# In[92]:


Tp is S2.tangent_bundle().fiber(p)


# The vector space $T_p \mathbb{S}^2$ is endowed with bases inherited from the coordinate frames defined around $p$:

# In[93]:


Tp.bases()


# On the contrary, since $(V,(x',y'))$ is the only chart defined so far around the point $N$, 
# we have a unique predefined basis in $T_N \mathbb{S}^2$:

# In[94]:


T_N = S2.tangent_space(N)
T_N.bases()


# To shorten some writings, there is the concept of default basis:

# In[95]:


Tp.default_basis()


# An element of $T_p\mathbb{S}^2$ is constructed via SageMath's *parent/element* syntax, i.e. via the call method of the parent:

# In[96]:


v = Tp((-2, 3), name='v')
print(v)


# Equivalently, one can use the method `tangent_vector` of manifolds:

# In[97]:


v == S2.tangent_vector(p, -2, 3, name='v')


# One has of course:

# In[98]:


v in Tp


# In[99]:


v.parent()


# The vector $v$ expanded in the default basis of $T_p \mathbb{S}^2$:

# In[100]:


v.display()


# Expansion in other bases:

# In[101]:


v.display(Tp.bases()[1])


# In[102]:


v.display(Tp.bases()[2])


# Tangent vectors are endowed with a method `plot`:

# In[103]:


graph += v.plot(chart=cartesian, mapping=Phi, scale=0.2, width=0.5)
graph


# ### Differential of a smooth map
# 
# The differential of the map $\Phi$ at the point $p$ is

# In[104]:


dPhi_p = Phi.differential(p)
print(dPhi_p)
dPhi_p


# In[105]:


dPhi_p.domain()


# In[106]:


dPhi_p.codomain()


# In[107]:


dPhi_p.parent()


# The image by $\mathrm{d}\Phi_p$ of the vector $v\in T_p\mathbb{S}^2$ introduced above is

# In[108]:


dPhi_p(v)


# In[109]:


print(dPhi_p(v))


# In[110]:


dPhi_p(v) in R3.tangent_space(Phi(p))


# In[111]:


dPhi_p(v).display()


# ## Algebra of scalar fields
# 
# The set $C^\infty(\mathbb{S}^2)$ of all smooth functions $\mathbb{S}^2\rightarrow \mathbb{R}$ has naturally the structure of a commutative algebra over $\mathbb{R}$. $C^\infty(\mathbb{S}^2)$ is therefore returned by the method `scalar_field_algebra()` of manifolds:

# In[112]:


CS = S2.scalar_field_algebra()
CS


# <p>Since the algebra internal product is the pointwise multiplication, it is clearly commutative, so that $C^\infty(\mathbb{S}^2)$ belongs to Sage's category of commutative algebras:</p>

# In[113]:


print(CS.category())


# <p>The base ring of the algebra $C^\infty(\mathbb{S}^2)$ is the field $\mathbb{R}$, which is represented here by Sage's Symbolic Ring (SR):</p>

# In[114]:


CS.base_ring()


# <p>Elements of $C^\infty(\mathbb{S}^2)$ are of course (smooth) scalar fields:</p>

# In[115]:


print(CS.an_element())


# <p>This example element is the constant scalar field that takes the value 2:</p>

# In[116]:


CS.an_element().display()


# <p>A specific element is the zero one:</p>

# In[117]:


f = CS.zero()
print(f)


# <p>Scalar fields map points of $\mathbb{S}^2$ to real numbers:</p>

# In[118]:


f(N), f(S), f(p)


# In[119]:


f.display()


# <p>Another specific element is the algebra unit element, i.e. the constant scalar field 1:</p>

# In[120]:


f = CS.one()
print(f)


# In[121]:


f(N), f(S), f(p)


# A generic scalar field is defined by its coordinate expression in some chart(s); for instance:

# In[122]:


f = S2.scalar_field({stereoN: 1/(1+x^2+y^2)}, name='f')
f.display()


# We see that Sage has used the transition map between the two stereographic charts on $W$ to express $f$ in terms of the coordinates $(x',y')$ on $W$. Let us this expression to extend $f$ to the whole of $V$: 

# In[123]:


f.add_expr_by_continuation(stereoS, W)


# Then $f$ is well defined in all $\mathbb{S}^2 = U \cup V$:

# In[124]:


f.display()


# In[125]:


f(N)


# In[126]:


f.parent()


# <p>Scalar fields map the manifold's points to real numbers:</p>

# In[127]:


f(N)


# In[128]:


f(S)


# In[129]:


f(p)


# <p>We may define the restrictions of $f$ to the open subsets $U$ and $V$:</p>

# In[130]:


fU = f.restrict(U)
fU.display()


# In[131]:


fV = f.restrict(V)
fV.display()


# In[132]:


fU(p), fU(S)


# In[133]:


fU.parent()


# In[134]:


fV.parent()


# In[135]:


CU = U.scalar_field_algebra()
fU.parent() is CU


# <p>A scalar field on $\mathbb{S}^2$ can be coerced to a scalar field on $U$, the coercion being simply the restriction:</p>

# In[136]:


CU.has_coerce_map_from(CS)


# In[137]:


fU == CU(f)


# The arithmetic of scalar fields (operations in the algebra $C^\infty(\mathbb{S}^2)$):

# In[138]:


g = f*f - 2*f
g.set_name('g')
g.display()


# ## Vector fields
# 
# The set $\mathfrak{X}(\mathbb{S}^2)$ of all smooth vector fields on $\mathbb{S}^2$ is a module over the algebra $C^\infty(\mathbb{S}^2)$. It is obtained by the method `vector_field_module()`:

# In[139]:


XS = S2.vector_field_module()
XS


# In[140]:


print(XS)


# In[141]:


XS.base_ring()


# In[142]:


XS.category()


# <p>$\mathfrak{X}(\mathbb{S}^2)$ is not a free module:</p>

# In[143]:


isinstance(XS, FiniteRankFreeModule)


# <p>because $\mathbb{S}^2$ is not a parallelizable manifold:</p>

# In[144]:


S2.is_manifestly_parallelizable()


# <p>On the contrary, the set $\mathfrak{X}(U)$ of smooth vector fields on $U$ is a free module:</p>

# In[145]:


XU = U.vector_field_module()
isinstance(XU, FiniteRankFreeModule)


# <p>because $U$ is parallelizable:</p>

# In[146]:


U.is_manifestly_parallelizable()


# <p>Due to the introduction of the stereographic coordinates $(x,y)$ on $U$, a basis has already been defined on the free module $\mathfrak{X}(U)$, namely the coordinate basis $(\partial/\partial x, \partial/\partial y)$:</p>

# In[147]:


XU.print_bases()


# In[148]:


eU = XU.default_basis()
eU


# <p>Similarly</p>

# In[149]:


XV = V.vector_field_module()
eV = XV.default_basis()
eV


# From the point of view of the open set $U$, `eU` is also the default vector frame:

# In[150]:


eU is U.default_frame()


# Similarly:

# In[151]:


eV is V.default_frame()


# `eU` is also the default vector frame on $\mathbb{S}^2$ (although not defined on the whole $\mathbb{S}^2$), for it is the first vector frame defined on an open subset of $\mathbb{S}^2$:

# In[152]:


eU is S2.default_frame()


# In[153]:


S2.frames()


# Let us introduce a vector field on $\mathbb{S}^2$ by providing its components in the
# frame `eU`:

# In[154]:


v = S2.vector_field(1, -2, frame=eU, name='v')
v.display(eU)


# In[155]:


v.parent()


# On $W$, we can express $v$ in terms of the $(x',y')$ coordinates:

# In[156]:


v.restrict(W).display(stereoS.restrict(W).frame(), stereoS.restrict(W))


# We extend the definition of $v$ to $V$ thanks to the above expression:

# In[157]:


v.add_comp_by_continuation(eV, W, chart=stereoS)
v.display(eV)


# At this stage, the vector field $v$ is defined on the whole manifold $\mathbb{S}^2$: it has expressions in each of the two frames `eU` and `eV`, which cover $\mathbb{S}^2$.

# According to the hairy ball theorem, $v$ has to vanish somewhere. This occurs at the North pole:

# In[158]:


vN = v.at(N)
print(vN)


# In[159]:


vN.display()


# <p>$v|_N$ is the zero vector of the tangent vector space $T_N\mathbb{S}^2$:</p>

# In[160]:


vN.parent()


# In[161]:


vN.parent() is S2.tangent_space(N)


# In[162]:


vN == S2.tangent_space(N).zero()


# <p>On the contrary, $v$ is non-zero at the South pole:</p>

# In[163]:


vS = v.at(S)
print(v)


# In[164]:


vS.display()


# <p>Let us plot the vector field $v$ is terms of the stereographic chart $(U,(x,y))$, with the South pole $S$ superposed:</p>

# In[165]:


v.plot(chart=stereoN, chart_domain=stereoN, max_range=4, 
       number_values=5, scale=0.5, aspect_ratio=1) \
+ S.plot(stereoN, size=30, label_offset=0.2)


# <p>The vector field appears homogeneous because its components w.r.t. the frame $\left(\frac{\partial}{\partial x}, \frac{\partial}{\partial y}\right)$ are constant:</p>

# In[166]:


v.display(stereoN.frame())


# <p>On the contrary, once drawn in terms of the stereographic chart $(V, (x',y'))$, $v$ does no longer appears homogeneous:</p>

# In[167]:


v.plot(chart=stereoS, chart_domain=stereoS, max_range=4, scale=0.02, 
       aspect_ratio=1) \
+ N.plot(chart=stereoS, size=30, label_offset=0.2)


# <p>Finally, a 3D view of the vector field $v$ is obtained via the embedding $\Phi$:</p>

# In[168]:


graph_v = v.plot(chart=cartesian, mapping=Phi, chart_domain=spher, 
                 number_values=11, scale=0.2)
graph_spher + graph_v


# <p>Similarly, let us draw the first vector field of the stereographic frame from the North pole, namely $\frac{\partial}{\partial x}$:</p>

# In[169]:


stereoN.frame()


# In[170]:


ex = stereoN.frame()[1]
ex


# In[171]:


graph_ex = ex.plot(chart=cartesian, mapping=Phi, chart_domain=spher,
                   number_values=11, scale=0.4, width=1, 
                   label_axes=False)
graph_spher + graph_ex


# <p>For the second vector field of the stereographic frame from the North pole, namely $\frac{\partial}{\partial y}$, we get</p>

# In[172]:


ey = stereoN.frame()[2]
ey


# In[173]:


graph_ey = ey.plot(chart=cartesian, mapping=Phi, chart_domain=spher,
                   number_values=11, scale=0.4, width=1, color='red', 
                   label_axes=False)
graph_spher + graph_ey


# We may combine the two graphs, to get a 3D view of the vector frame associated with the stereographic coordinates from the North pole:

# In[174]:


graph_frame = graph_spher + graph_ex + graph_ey \
              + N.plot(cartesian, mapping=Phi, label_offset=0.05, size=5) \
              + S.plot(cartesian, mapping=Phi, label_offset=0.05, size=5)
graph_frame + sphere(color='lightgrey', opacity=0.4)


# The same scene rendered with Tachyon:

# In[175]:


show(graph_frame + sphere(opacity=0.5), viewer='tachyon', figsize=10)


# ### Vector fields acting on scalar fields
# 
# $v$ and $f$ are both fields defined on the whole sphere (respectively a vector field and a scalar field). By the very definition of a vector field, $v$ acts on $f$:

# In[176]:


vf = v(f)
print(vf)
vf.display()


# Values of $v(f)$ at the North pole, at the South pole and at point $p$:

# In[177]:


vf(N)


# In[178]:


vf(S)


# In[179]:


vf(p)


# ## 1-forms
# 
# A 1-form on $\mathbb{S}^2$ is a field of linear forms on the tangent spaces. For instance it can be the differential of a scalar field:

# In[180]:


f.display()


# In[181]:


df = diff(f)
print(df)


# In[182]:


df.display()  # display w.r.t. the default frame


# In[183]:


df.display(eV)


# In[184]:


df.display(spher.frame())


# In[185]:


df.display(spher.frame(), spher)  # asking for the components to be shown in the spherical chart


# In[186]:


print(df.parent())


# In[187]:


df.parent()


# <p><span id="cell_outer_146">The 1-form acting on a vector field:</span></p>

# In[188]:


print(df(v))
df(v).display()


# <p>Let us check the identity $\mathrm{d}f(v) = v(f)$:</p>

# In[189]:


df(v) == v(f)


# <p>Similarly, we have $\mathcal{L}_v f = v(f)$:</p>

# In[190]:


f.lie_derivative(v) == v(f)


# ## Curves in $\mathbb{S}^2$
# 
# In order to define curves in $\mathbb{S}^2$, we first introduce the field of real numbers $\mathbb{R}$ as a 1-dimensional smooth manifold with a canonical coordinate chart:

# In[191]:


R.<t> = manifolds.RealLine()
print(R)


# In[192]:


print(R.category())


# In[193]:


dim(R)


# In[194]:


R.atlas()


# <p>Let us define a <strong>loxodrome of the sphere</strong> in terms of its parametric equation with respect to the chart <span style="font-family: courier new,courier;">spher</span> = $(A,(\theta,\phi))$</p>

# In[195]:


c = S2.curve({spher: [2*atan(exp(-t/10)), t]}, (t, -oo, +oo), name='c')


# <p>Curves in $\mathbb{S}^2$ are considered as morphisms from the manifold $\mathbb{R}$ to the manifold $\mathbb{S}^2$:</p>

# In[196]:


c.parent()


# In[197]:


c.display()


# <p>The curve $c$ can be plotted in terms of stereographic coordinates $(x,y)$:</p>

# In[198]:


c.plot(chart=stereoN, aspect_ratio=1)


# <p>We recover the well-known fact that the graph of a loxodrome in terms of stereographic coordinates is a <strong>logarithmic spiral</strong>.</p>
# <p>Thanks to the embedding $\Phi$, we may also plot $c$ in terms of the Cartesian coordinates of $\mathbb{R}^3$:</p>

# In[199]:


graph_c = c.plot(mapping=Phi, max_range=40, plot_points=200, 
                 thickness=2, label_axes=False)
graph_spher + graph_c


# <p>The <strong>tangent vector field</strong> (or <strong>velocity vector</strong>) to the curve $c$ is</p>

# In[200]:


vc = c.tangent_vector_field()
vc


# <p>$c'$ is a vector field <em>along</em> $\mathbb{R}$ taking its values in tangent spaces to $\mathbb{S}^2$:</p>

# In[201]:


print(vc)


# <p>The set of vector fields along $\mathbb{R}$ taking their values on $\mathbb{S}^2$ via the differential mapping $c: \mathbb{R} \rightarrow \mathbb{S}^2$ is denoted by $\mathfrak{X}(\mathbb{R},c)$; it is a module over the algebra $C^\infty(\mathbb{R})$:</p>

# In[202]:


vc.parent()


# In[203]:


vc.parent().category()


# In[204]:


vc.parent().base_ring()


# <p>A coordinate view of $c'$:</p>

# In[205]:


vc.display()


# <p>Let us plot the vector field $c'$ in terms of the stereographic chart $(U,(x,y))$:</p>

# In[206]:


show(vc.plot(chart=stereoN, number_values=30, scale=0.5, color='red') +
     c.plot(chart=stereoN), aspect_ratio=1)


# <p>A 3D view of $c'$ is obtained via the embedding $\Phi$:</p>

# In[207]:


graph_vc = vc.plot(chart=cartesian, mapping=Phi, ranges={t: (-20, 20)}, 
                   number_values=30, scale=0.5, color='red', 
                   label_axes=False)
graph_spher + graph_c + graph_vc


# ## Riemannian metric on $\mathbb{S}^2$
# 
# The standard metric on $\mathbb{S}^2$ is that induced by the Euclidean metric of $\mathbb{R}^3$. The latter
# is obtained by:

# In[208]:


h = R3.metric()
h.display()


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

# In[209]:


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


# Note that we could have defined $g$ intrinsically, i.e. by providing its components in the two frames `stereoN` and `stereoS`. Instead, we have chosen to get it as the pullback of $h$, since the pullback computation is implemented in SageMath.
# 
# The metric is a symmetric tensor field of type (0,2):

# In[210]:


print(g.parent())


# In[211]:


g.tensor_type()


# In[212]:


g.symmetries()


# The expression of the metric in terms of the default frame on $\mathbb{S}^2$ (stereoN):

# In[213]:


g.display()


# We may factorize the metric components:

# In[214]:


g.apply_map(factor, frame=eU, keep_other_components=True)
g.display()


# A matrix view of the components of $g$ in the manifold's default frame:

# In[215]:


g[:]


# In[216]:


g[1,1]


# Display in terms of the vector frame `eV` =  $(V, (\partial_{x'}, \partial_{y'}))$:

# In[217]:


g.apply_map(factor, frame=eV, keep_other_components=True)
g.display(eV)


# Expression of the metric tensor in terms of spherical coordinates:

# In[218]:


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


# The metric acts on vector field pairs, resulting in a scalar field:

# In[219]:


print(g(v,v))


# In[220]:


g(v,v).parent()


# In[221]:


g(v,v).display()


# The **Levi-Civita connection** associated with the metric $g$:

# In[222]:


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


# As a test, we verify that $\nabla_g$ acting on $g$ results in zero:

# In[223]:


nabla(g).display()


# The nonzero Christoffel symbols of $g$ (skipping those that can be deduced by symmetry on the last two indices) w.r.t. two charts:

# In[224]:


g.christoffel_symbols_display(chart=stereoN)


# In[225]:


g.christoffel_symbols_display(chart=spher)


# $\nabla_g$ acting on the vector field $v$:

# In[226]:


print(nabla(v))


# In[227]:


nabla(v).display(stereoN.frame())


# In[228]:


nabla(v)[:]


# In[229]:


nabla(v)[1,2]


# ## Curvature
# 
# The Riemann tensor associated with the metric $g$:

# In[230]:


Riem = g.riemann()
print(Riem)
Riem.display()


# The components of the Riemann tensor in the default frame on $\mathbb{S}^2$:

# In[231]:


Riem.display_comp()


# <p>The components in the frame associated with spherical coordinates:</p>

# In[232]:


Riem.display_comp(spher.frame(), chart=spher)


# In[233]:


print(Riem.parent())


# In[234]:


Riem.symmetries()


# The Riemann tensor associated with the Euclidean metric $h$ on $\mathbb{R}^3$ is identically zero:

# In[235]:


h.riemann().display()


# The Ricci tensor and the Ricci scalar:

# In[236]:


Ric = g.ricci()
Ric.display()


# In[237]:


print(g.inverse())


# In[238]:


g.inverse().display()


# In[239]:


R = g.ricci_scalar()
R.display()


# Hence we recover the fact that $(\mathbb{S}^2,g)$ is a Riemannian manifold of constant positive curvature.
# 
# In dimension 2, the Riemann curvature tensor is entirely determined by the Ricci scalar $R$ according to
# $$ R^i_{\ \, jlk} = \frac{R}{2} \left( \delta^i_{\ \, k} g_{jl} - \delta^i_{\ \, l} g_{jk} \right)$$
# Let us check this formula here, under the form $R^i_{\ \, jlk} = -R g_{j[k} \delta^i_{\ \, l]}$:

# In[240]:


delta = S2.tangent_identity_field()
Riem == - R*(g*delta).antisymmetrize(2,3)


# Similarly the relation $\mathrm{Ric} = (R/2)\; g$ must hold:

# In[241]:


Ric == (R/2)*g


# ### Manifold orientation and volume 2-form

# In order to introduce the volume 2-form associated with the metric $g$, we need to define an orientation on $\mathbb{S}^2$ first. We choose the orientation so that the vector frame $(\partial/\partial x', \partial/\partial y')$ of the stereographic coordinates from the South pole is right-handed. This is somewhat natural, because the triplet $(\partial/\partial x', \partial/\partial y', n)$, where $n$ is the unit outward normal to $\mathbb{S}^2$, is right-handed with respect to the standard orientation of $\mathbb{R}^3$. On the contrary the triplet
# $(\partial/\partial x, \partial/\partial y, n)$ formed from stereographic coordinates from the North pole is left-handed (see the above plot). Actually, we can check that $(\partial/\partial x, \partial/\partial y)$
# and $(\partial/\partial x', \partial/\partial y')$ lead to two opposite orientations, because the transition map
# $(x, y) \mapsto (x', y')$ has a negative Jacobian determinant:

# In[242]:


stereoN_to_S.jacobian_det()


# We define the orientation via the method `set_orientation()` with a list of right-handed vector frames, whose domains form an open cover of $\mathbb{S}^2$. We therefore provide `eV` = $(\partial/\partial x', \partial/\partial y')$ (domain: $V$) and the "reversed" frame $(\partial/\partial y, \partial/\partial x)$ on $U$:

# In[243]:


reU = U.vector_frame('f', (eU[2], eU[1]))
reU[1].display(eU), reU[2].display(eU)


# In[244]:


S2.set_orientation([eV, reU])


# The **volume 2-form** or **Levi-Civita tensor** associated with $g$ is then

# In[245]:


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


# Notice the minus sign in the the above expression, which reflects the fact that the default frame $(\partial/\partial x, \partial/\partial y)$ is left-handed. On the contrary, we have

# In[246]:


eps.display(eV)


# A nicer display is obtained by factorizing the components:

# In[247]:


eps.apply_map(factor, frame=eV, keep_other_components=True)
eps.display(stereoS.frame())


# The frame associated with spherical coordinates is right-handed and we recover the standard expression of the volume 2-form:

# In[248]:


eps.display(spher.frame(), chart=spher)


# The exterior derivative of the 2-form $\epsilon_g$:

# In[249]:


print(diff(eps))


# Of course, since $\mathbb{S}^2$ has dimension 2, all 3-forms vanish identically:

# In[250]:


diff(eps).display()


# ### Non-holonomic frames
# Up to know, all the vector frames introduced on $\mathbb{S}^2$ have been coordinate frames. Let us introduce a non-coordinate frame on the open subset $A$. To ease the manipulations, we change first the default chart and default frame on $A$ to the spherical coordinate ones:

# In[251]:


A.default_chart()


# In[252]:


A.default_frame()


# In[253]:


A.set_default_chart(spher)
A.set_default_frame(spher.frame())
A.default_chart()


# In[254]:


A.default_frame()


# We introduce the new vector frame $e =  \left(\frac{\partial}{\partial\theta}, \frac{1}{\sin\theta}\frac{\partial}{\partial\phi}\right)$:

# In[255]:


spher.frame()[:]


# In[256]:


d_dth, d_dph = spher.frame()[:]
e = A.vector_frame('e', (d_dth, 1/sin(th)*d_dph))
print(e)
e


# In[257]:


(e[1].display(), e[2].display())


# <p>The new frame is an orthonormal frame for the metric $g$:</p>

# In[258]:


g(e[1],e[1]).expr()


# In[259]:


g(e[1],e[2]).expr()


# In[260]:


g(e[2],e[2]).expr()


# In[261]:


g[e,:]


# In[262]:


g.display(e)


# In[263]:


eps.display(e)


# It is non-holonomic, since its structure coefficients are not identically zero:

# In[264]:


e.structure_coeff()[:]


# In[265]:


e[2].lie_derivative(e[1]).display(e)


# <p>while we have of course</p>

# In[266]:


spher.frame().structure_coeff()[:]


# ## Using SymPy as the symbolic backend
# 
# By default, the symbolic backend used in calculus on manifolds is SageMath's one (Pynac + Maxima), implemented via the symbolic ring `SR`. We can choose to use [SymPy](https://www.sympy.org/) instead:

# In[267]:


S2.set_calculus_method('sympy')


# In[268]:


F = 2*f
F.display()


# In[269]:


F.expr()


# In[270]:


type(F.expr())


# Back to Sage's default:

# In[271]:


S2.set_calculus_method('SR')


# In[272]:


F.expr()


# In[273]:


type(F.expr())


# ## Going further
# 
# See the notebooks [Smooth manifolds, charts and scalar fields](https://nbviewer.org/github/sagemanifolds/SageManifolds/blob/master/Worksheets/JNCF2018/jncf18_scalar.ipynb) and [Smooth manifolds, vector fields and tensor fields](https://nbviewer.org/github/sagemanifolds/SageManifolds/blob/master/Worksheets/JNCF2018/jncf18_vector.ipynb) from the lectures [Symbolic tensor calculus on manifolds](https://sagemanifolds.obspm.fr/jncf2018/). Many example notebooks are 
# provided at the [SageManifolds page](https://sagemanifolds.obspm.fr/examples.html).
# 
# See also the series of notebooks by Andrzej Chrzeszczyk: [Introduction to manifolds in SageMath](https://sagemanifolds.obspm.fr/intro_to_manifolds.html), as well as the tutorial videos by Christian Bär: [Manifolds in SageMath](https://www.youtube.com/playlist?list=PLnrOCYZpQUuJlnQbQ48zgGk-Ks1t145Yw).