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 project.

Click here to download the notebook file (ipynb format). To run it, you must start SageMath with the Jupyter interface, via the command sage -n jupyter

NB: a version of SageMath at least equal to 9.2 is required to run this notebook:

version()

'SageMath version 9.3.beta2, Release Date: 2020-11-24'

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

%display latex

and we initialize a time counter for benchmarking:

import time
comput_time0 = time.perf_counter()

$\mathbb{S}^2$ as a 2-dimensional differentiable manifold

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

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\}$.

print(S2)

2-dimensional differentiable manifold S^2

S2

$$\newcommand{\Bold}[1]{\mathbf{#1}}\mathbb{S}^2$$

The manifold is a Parent object:

isinstance(S2, Parent)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\mathrm{True}$$

in the category of smooth manifolds over $\mathbb{R}$:

S2.category()

$$\newcommand{\Bold}[1]{\mathbf{#1}}\mathbf{Smooth}_{\Bold{R}}$$

Coordinate charts on $\mathbb{S}^2$

The sphere cannot be covered by a single chart. At least two charts are necessary, for instance the charts associated with the stereographic projections from the North pole and the South pole respectively. Let us 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:

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

Open subset U of the 2-dimensional differentiable manifold S^2

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

Open subset V of the 2-dimensional differentiable manifold S^2

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

S2.declare_union(U, V)

Then we declare the stereographic chart on $U$, denoting by $(x,y)$ the coordinates resulting from the stereographic projection from the North pole:

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

The expression .<x,y> in the left-hand side means that the Python variables x and y are set to the two coordinates of the chart. This allows one to refer subsequently to the coordinates by their names. 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.

stereoN

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left(U,(x, y)\right)$$

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:

stereoN[1]

$$\newcommand{\Bold}[1]{\mathbf{#1}}x$$

y is stereoN[2]

$$\newcommand{\Bold}[1]{\mathbf{#1}}\mathrm{True}$$

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

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'$.

stereoS

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left(V,({x'}, {y'})\right)$$

At this stage, the user's atlas on the manifold has two charts:

S2.atlas()

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left(U,(x, y)\right), \left(V,({x'}, {y'})\right)\right]$$

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

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

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left\{\begin{array}{lcl} {x'} & = & \frac{x}{x^{2} + y^{2}} \\ {y'} & = & \frac{y}{x^{2} + y^{2}} \end{array}\right.$$

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

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

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left\{\begin{array}{lcl} x & = & \frac{{x'}}{{x'}^{2} + {y'}^{2}} \\ y & = & \frac{{y'}}{{x'}^{2} + {y'}^{2}} \end{array}\right.$$

In the present case, the situation is of course perfectly symmetric regarding the coordinates $(x,y)$ and $(x',y')$.

At this stage, the user's atlas has four charts:

S2.atlas()

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left(U,(x, y)\right), \left(V,({x'}, {y'})\right), \left(W,(x, y)\right), \left(W,({x'}, {y'})\right)\right]$$

Let us store $W = U\cap V$ into a Python variable for future use:

W = U.intersection(V)

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:

stereoN_W = stereoN.restrict(W)
stereoN_W

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left(W,(x, y)\right)$$

stereoS_W = stereoS.restrict(W)
stereoS_W

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left(W,({x'}, {y'})\right)$$

We may plot the chart $(W, (x',y'))$ in terms of itself, as a grid:

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

graphSN1 = stereoS_W.plot(stereoN, ranges={xp:[-6,-0.02], yp:[-6,-0.02]})
graphSN2 = stereoS_W.plot(stereoN, ranges={xp:[-6,-0.02], yp:[0.02,6]})
graphSN3 = stereoS_W.plot(stereoN, ranges={xp:[0.02,6], yp:[-6,-0.02]})
graphSN4 = stereoS_W.plot(stereoN, ranges={xp:[0.02,6], yp:[0.02,6]})
show(graphSN1+graphSN2+graphSN3+graphSN4,
     xmin=-1.5, xmax=1.5, ymin=-1.5, ymax=1.5)

Spherical coordinates

The standard spherical (or polar) 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:

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

Open subset A of the 2-dimensional differentiable manifold S^2

The restriction of the stereographic chart from the North pole to $A$ is

stereoN_A = stereoN_W.restrict(A)
stereoN_A

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left(A,(x, y)\right)$$

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

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

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left(A,({\theta}, {\phi})\right)$$

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

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

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left\{\begin{array}{lcl} x & = & -\frac{\cos\left({\phi}\right) \sin\left({\theta}\right)}{\cos\left({\theta}\right) - 1} \\ y & = & -\frac{\sin\left({\phi}\right) \sin\left({\theta}\right)}{\cos\left({\theta}\right) - 1} \end{array}\right.$$

We also provide the inverse transition map:

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

Check of the inverse coordinate transformation:
  th == 2*arctan(sqrt(-cos(th) + 1)/sqrt(cos(th) + 1))  **failed**
  ph == pi + arctan2(sin(ph)*sin(th)/(cos(th) - 1), cos(ph)*sin(th)/(cos(th) - 1))  **failed**
  x == x  *passed*
  y == y  *passed*
NB: a failed report can reflect a mere lack of simplification.

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

spher_to_stereoN.inverse().display()

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left\{\begin{array}{lcl} {\theta} & = & 2 \, \arctan\left(\frac{1}{\sqrt{x^{2} + y^{2}}}\right) \\ {\phi} & = & \pi + \arctan\left(-y, -x\right) \end{array}\right.$$

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'))$:

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

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left\{\begin{array}{lcl} {x'} & = & -\frac{\cos\left({\phi}\right) \cos\left({\theta}\right) - \cos\left({\phi}\right)}{\sin\left({\theta}\right)} \\ {y'} & = & -\frac{\cos\left({\theta}\right) \sin\left({\phi}\right) - \sin\left({\phi}\right)}{\sin\left({\theta}\right)} \end{array}\right.$$

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

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

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left\{\begin{array}{lcl} {\theta} & = & 2 \, \arctan\left(\sqrt{{x'}^{2} + {y'}^{2}}\right) \\ {\phi} & = & \pi - \arctan\left(\frac{{y'}}{{x'}^{2} + {y'}^{2}}, -\frac{{x'}}{{x'}^{2} + {y'}^{2}}\right) \end{array}\right.$$

The user atlas of $\mathbb{S}^2$ is now

S2.atlas()

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left(U,(x, y)\right), \left(V,({x'}, {y'})\right), \left(W,(x, y)\right), \left(W,({x'}, {y'})\right), \left(A,(x, y)\right), \left(A,({x'}, {y'})\right), \left(A,({\theta}, {\phi})\right)\right]$$

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

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

Points on $\mathbb{S}^2$

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

N = V.point((0,0), chart=stereoS, name='N') ; print(N)
S = U.point((0,0), chart=stereoN, name='S') ; print(S)

Point N on the 2-dimensional differentiable manifold S^2
Point S on the 2-dimensional differentiable manifold S^2

Since points are Sage Element's, the corresponding Parent being the manifold subsets, an equivalent writing of the above declarations is

N = V((0,0), chart=stereoS, name='N') ; print(N)
S = U((0,0), chart=stereoN, name='S') ; print(S)

Point N on the 2-dimensional differentiable manifold S^2
Point S on the 2-dimensional differentiable manifold S^2

Moreover, since stereoS in the default chart on $V$ and stereoN is the default one on $U$, their mentions can be omitted, so that the above can be shortened to

N = V((0,0), name='N') ; print(N)
S = U((0,0), name='S') ; print(S)

Point N on the 2-dimensional differentiable manifold S^2
Point S on the 2-dimensional differentiable manifold S^2

N.parent()

$$\newcommand{\Bold}[1]{\mathbf{#1}}V$$

S.parent()

$$\newcommand{\Bold}[1]{\mathbf{#1}}U$$

We have of course

N in V

$$\newcommand{\Bold}[1]{\mathbf{#1}}\mathrm{True}$$

N in S2

$$\newcommand{\Bold}[1]{\mathbf{#1}}\mathrm{True}$$

N in U

$$\newcommand{\Bold}[1]{\mathbf{#1}}\mathrm{False}$$

N in A

$$\newcommand{\Bold}[1]{\mathbf{#1}}\mathrm{False}$$

Let us introduce some point at the equator:

E = S2((0,1), chart=stereoN, name='E')

The point $E$ is in the open subset $A$:

E in A

$$\newcommand{\Bold}[1]{\mathbf{#1}}\mathrm{True}$$

We may then ask for its spherical coordinates $(\theta,\phi)$:

E.coord(spher)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\frac{1}{2} \, \pi, \frac{1}{2} \, \pi\right)$$

which is not possible for the point $N$:

try:
    N.coord(spher)
except ValueError as exc:
    print('Error: ' + str(exc))

Error: the point does not belong to the domain of Chart (A, (th, ph))

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

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

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\mathbb{R}^3,(X, Y, Z)\right)$$

The embedding of $\mathbb{S}^2$ into $\mathbb{R}^3$ is then defined by 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$:

Phi = S2.diff_map(R3, {(stereoN, cart): 
                       [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, cart): 
                       [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')

Phi.display()

$$\newcommand{\Bold}[1]{\mathbf{#1}}\begin{array}{llcl} \Phi:& \mathbb{S}^2 & \longrightarrow & \mathbb{R}^3 \\ \mbox{on}\ U : & \left(x, y\right) & \longmapsto & \left(X, Y, Z\right) = \left(\frac{2 \, x}{x^{2} + y^{2} + 1}, \frac{2 \, y}{x^{2} + y^{2} + 1}, \frac{x^{2} + y^{2} - 1}{x^{2} + y^{2} + 1}\right) \\ \mbox{on}\ V : & \left({x'}, {y'}\right) & \longmapsto & \left(X, Y, Z\right) = \left(\frac{2 \, {x'}}{{x'}^{2} + {y'}^{2} + 1}, \frac{2 \, {y'}}{{x'}^{2} + {y'}^{2} + 1}, -\frac{{x'}^{2} + {y'}^{2} - 1}{{x'}^{2} + {y'}^{2} + 1}\right) \end{array}$$

Phi.parent()

$$\newcommand{\Bold}[1]{\mathbf{#1}}\mathrm{Hom}\left(\mathbb{S}^2,\mathbb{R}^3\right)$$

print(Phi.parent())

Set of Morphisms from 2-dimensional differentiable manifold S^2 to Euclidean space R^3 in Category of smooth manifolds over Real Field with 53 bits of precision

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

$$\newcommand{\Bold}[1]{\mathbf{#1}}\mathrm{True}$$

$\Phi$ maps points of $\mathbb{S}^2$ to points of $\mathbb{R}^3$:

N1 = Phi(N) ; print(N1) ; N1 ; N1.coord()

Point Phi(N) on the Euclidean space R^3

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left(0, 0, 1\right)$$

S1 = Phi(S) ; print(S1) ; S1 ; S1.coord()

Point Phi(S) on the Euclidean space R^3

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left(0, 0, -1\right)$$

E1 = Phi(E) ; print(E1) ; E1 ; E1.coord()

Point Phi(E) on the Euclidean space R^3

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left(0, 1, 0\right)$$

$\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. The latter is then computed by means of the transition map $(A,(x,y))\rightarrow (A,(\theta,\phi))$:

Phi.expr(stereoN_A, cart)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\frac{2 \, x}{x^{2} + y^{2} + 1}, \frac{2 \, y}{x^{2} + y^{2} + 1}, \frac{x^{2} + y^{2} - 1}{x^{2} + y^{2} + 1}\right)$$

Phi.expr(spher, cart)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\cos\left({\phi}\right) \sin\left({\theta}\right), \sin\left({\phi}\right) \sin\left({\theta}\right), \cos\left({\theta}\right)\right)$$

Phi.display(spher, cart)

$$\newcommand{\Bold}[1]{\mathbf{#1}}\begin{array}{llcl} \Phi:& \mathbb{S}^2 & \longrightarrow & \mathbb{R}^3 \\ \mbox{on}\ A : & \left({\theta}, {\phi}\right) & \longmapsto & \left(X, Y, Z\right) = \left(\cos\left({\phi}\right) \sin\left({\theta}\right), \sin\left({\phi}\right) \sin\left({\theta}\right), \cos\left({\theta}\right)\right) \end{array}$$

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$:

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

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:

graph_stereoN = stereoN.plot(chart=cart, mapping=Phi, number_values=25, 
                             label_axes=False)
graph_stereoN

and then have a view with the stereographic coordinates from the South pole superposed (in green):

graph_stereoS = stereoS.plot(chart=cart, mapping=Phi, number_values=25, 
                             color='green', label_axes=False)
graph_stereoN + graph_stereoS