{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 22. Riemannian and pseudo-Riemannian manifolds\n",
"\n",
"This notebook is part of the [Introduction to manifolds in SageMath](https://sagemanifolds.obspm.fr/intro_to_manifolds.html) by Andrzej Chrzeszczyk (Jan Kochanowski University of Kielce, Poland)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'SageMath version 9.6, Release Date: 2022-05-15'"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"version()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let $M$ be a smooth manifold and $g$ a covariant tensor field of rank two: $g\\in T^{(0,2)}M.$ We call $g$ **symmetric** if\n",
"$$g_p(X_p,Y_p)=g_p(Y_p,X_p)\\quad \\text{for } X_p,Y_p\\in T_pM,\\ p\\in M.$$\n",
"\n",
"We call $g$ **positive definite** if \n",
"$$g_p(X_p,X_p)\\geq 0\\quad \\text{for all}\\ X_p\\in T_pM,\\ p\\in M$$ and \n",
"$$g_p(X_p,X_p)=0\\quad \\text{implies } X_p=0_p\\in T_pM.\n",
"$$\n",
"\n",
"The tensor field $g$ is **non-singular** if for all $p\\in M$\n",
"\n",
"$$g_p(X_p,Y_p)=0\\quad \\text{for all } Y_p\\in T_pM\\quad \\text{implies } X_p=0_p\\in T_pM.\n",
"$$\n",
"\n",
"Note that the positive-definiteness implies the non-singularity:
\n",
"if $g_p(X_p,Y_p)=0\\quad \\text{for all } Y_p\\in T_pM,\\ $ then \n",
"$\\ g_p(X_p,X_p)=0\\ $ and consequently $X_p=0_p$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A **pseudo-Riemannian manifold** is a smooth manifold $M$ with a non-singular, symmetric, smooth tensor field $g\\in T^{(0,2)}M$, called **metric** of $M$.\n",
"\n",
"If the metric is positive definite we use the name **Riemannian manifold**."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
\n",
"\n",
"**Example 22.1**\n",
"\n",
"Metric method in SageMath:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2-dimensional Riemannian manifold M"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
" # Riemannian manifold M\n",
"M = Manifold(2, 'M', structure='Riemannian'); M"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Riemannian metric g on the 2-dimensional Riemannian manifold M"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g = M.metric(); g # metric on M"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
\n",
"\n",
"In a Riemannian manifold, $M$, with a metric $g$, the bilinear form $g_p(\\cdot,\\cdot)$ is an **inner product on** $T_pM$. The **norm or length** of a tangent vector $X_p ∈ T_p M$ is defined by \n",
"\n",
"\\begin{equation}\n",
"\\|X_p\\| = \\sqrt{g_p (X_p , X_p )},\n",
"\\tag{22.1}\n",
"\\end{equation}\n",
"\n",
"and the **length of a curve** $\\gamma : [a, b] → M$ is defined by\n",
"\n",
"\\begin{equation}\n",
"L_\\gamma=\\int_a^b\\|\\gamma'_t\\|dt.\n",
"\\tag{22.2}\n",
"\\end{equation}\n",
"\n",
"If $M$ is a Riemannian manifold and $(x^1,\\ldots,x^n)$ are local coordinates, then the metric is given by\n",
"\n",
"\\begin{equation}\n",
"g = g_{i j} dx^i ⊗ dx^j,\n",
"\\tag{22.3}\n",
"\\end{equation}\n",
"\n",
"where $g_{ij}=g(\\frac{\\partial}{\\partial x^i},\\frac{\\partial}{\\partial x^j})$ (this is a special case of (13.7) ).\n",
"\n",
"\n",
"The **standard metric on $R^n$** with Cartesian coordinates is defined by\n",
"\n",
"\\begin{equation}\n",
"g = \\delta_{ij}dx^i\\otimes dx^j=dx^1 ⊗ dx^1 + · · · + dx^n ⊗ dx^n ,\n",
"\\tag{22.4}\n",
"\\end{equation}\n",
"\n",
"
\n",
"**Example 22.2**\n",
"\n",
"Let us start from the standard metric in 2-dimensional Euclidean space"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"g = dx⊗dx + dy⊗dy"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"E. = EuclideanSpace() # Euclidean space E^2\n",
"g = E.metric() # standard metric on E^2\n",
"g.disp()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and use (22.2) to compute the length of the curve c defined in Cartesian coordinates by\n",
"\n",
"$$x(t)=\\sin(t),\\quad y(t)=\\sin(2t)/2,\\quad t\\in (0,2\\pi),$$\n",
"\n",
"(in simple cases the **exact** integral can be computed in SageMath, but **not in the present case**)."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"6.097223470104982"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"E. = EuclideanSpace() # Euclidean space E^2\n",
"t = var('t') # symbolic variable for c \n",
"c = E.curve({E.cartesian_coordinates():[sin(t), sin(2*t)/2]},\n",
" (t, 0, 2*pi), name='c') # define curve in Cart. coord.\n",
"v=c.tangent_vector_field() # vector field of tangent vect. to c\n",
"w=v.norm().expr() # norm of v\n",
"numerical_integral(w,0,2*pi)[0] # numerical version of (22.2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If the curve is simple, for example \n",
"\n",
"$$x(t)=\\sin(2t),\\quad y(t)=\\cos(2t),\\quad t\\in (0,2\\pi),$$\n",
"\n",
"then we don't need numerical tools:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\\(\\displaystyle 4 \\, \\pi\\)"
],
"text/latex": [
"$\\displaystyle 4 \\, \\pi$"
],
"text/plain": [
"4*pi"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"%display latex\n",
"E. = EuclideanSpace() # Euclidean space E^2\n",
"t = var('t') # symbolic variable for c \n",
"c = E.curve({E.cartesian_coordinates():[sin(2*t), cos(2*t)]},\n",
" (t, 0, 2*pi), name='c') # define curve in Cart. coord.\n",
"v=c.tangent_vector_field() # vector field of tangent vect. to c\n",
"w=v.norm().expr() # norm of v\n",
"w.integral(t,0,2*pi) # exact version of (22.2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
\n",
"\n",
"**Example 22.3**\n",
"\n",
"The standard metric in 4-dimensional Euclidean space is predefined."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\\(\\displaystyle g = \\mathrm{d} {x_{1}}\\otimes \\mathrm{d} {x_{1}}+\\mathrm{d} {x_{2}}\\otimes \\mathrm{d} {x_{2}}+\\mathrm{d} {x_{3}}\\otimes \\mathrm{d} {x_{3}}+\\mathrm{d} {x_{4}}\\otimes \\mathrm{d} {x_{4}}\\)"
],
"text/latex": [
"$\\displaystyle g = \\mathrm{d} {x_{1}}\\otimes \\mathrm{d} {x_{1}}+\\mathrm{d} {x_{2}}\\otimes \\mathrm{d} {x_{2}}+\\mathrm{d} {x_{3}}\\otimes \\mathrm{d} {x_{3}}+\\mathrm{d} {x_{4}}\\otimes \\mathrm{d} {x_{4}}$"
],
"text/plain": [
"g = dx1⊗dx1 + dx2⊗dx2 + dx3⊗dx3 + dx4⊗dx4"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"%display latex\n",
"E=EuclideanSpace(4) # Euclidean space E^4\n",
"E.metric().disp() # standard metric on E^4"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we need upper indices and more general manifolds, then we can use the commands:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\\(\\displaystyle g = \\mathrm{d} {x^{0}}\\otimes \\mathrm{d} {x^{0}}+\\mathrm{d} {x^{1}}\\otimes \\mathrm{d} {x^{1}}+\\mathrm{d} {x^{2}}\\otimes \\mathrm{d} {x^{2}}\\)"
],
"text/latex": [
"$\\displaystyle g = \\mathrm{d} {x^{0}}\\otimes \\mathrm{d} {x^{0}}+\\mathrm{d} {x^{1}}\\otimes \\mathrm{d} {x^{1}}+\\mathrm{d} {x^{2}}\\otimes \\mathrm{d} {x^{2}}$"
],
"text/plain": [
"g = dx0⊗dx0 + dx1⊗dx1 + dx2⊗dx2"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"N = 3 # dimension of manifold\n",
" # variables with superscripts:\n",
"M = Manifold(N, 'M') # manifold M\n",
"X = M.chart(' '.join(['x'+str(i)+':x^{'+str(i)+'}' for i in range(N)])) # chart on M\n",
"g = M.metric('g') # metric g on M \n",
"g[0,0], g[1,1], g[2,2] = 1, 1, 1 # components of g\n",
"g.disp() # show g"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
\n",
"\n",
"### Determinant of of $[g_{ij}]$\n",
"\n",
"
\n",
"\n",
"If in the implication \n",
"$$g_p(X_p,Y_p)=0\\ \\text{for all } Y_p\\quad \\Rightarrow X_p=0,\n",
"$$\n",
"\n",
"we put $X_p=a^i\\frac{\\partial}{\\partial x^i}\\big|_p, \\quad Y_p=\\frac{\\partial}{\\partial x^j}\\big|_p$, then we obtain\n",
"\n",
"$$g_{i j} ( p)\\, a^i = 0,\\quad j=1,\\ldots,n\\quad \\Rightarrow a^i=0, \\ i=1,\\ldots,n,$$\n",
"\n",
"which means that the homogeneous system of linear equations for the unknowns $a^i$ admits only the zero solutions i.e., \n",
"\n",
"$$\\det [g_{i,j}](p)\\not=0,$$\n",
"for all $p$ in the the coordinate domain."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
\n",
"\n",
"### Pullback of a metric\n",
"\n",
"
\n",
"\n",
"If $M$ is a Riemannian manifold with a positive definite metric tensor $g$ and $ψ :N → M $ is a smooth map such that for all $p ∈\n",
"N ,\\ dψ_p$ has **maximal rank**, that is,\n",
"\n",
" $$dψ_p X_p = 0_{ψ( p)}\\quad \\text{implies}\\ \\ X_p = 0_p,$$\n",
"\n",
"then **the pullback $ψ^∗ g$ is\n",
"a positive definite metric tensor in $N$.**\n",
"\n",
"In fact since $(ψ^∗ g)_p (X_p , Y_p ) ≡ g_{ψ( p)} (dψ_p X_p , dψ_p Y_p )$, and $g$ is symmetric then $\\psi^*g$ is symmetric.\n",
"Furthermore, if $(ψ^*g)_p (X_p , Y_p ) = 0$ for all $Y_p ∈ T_p N$ , from the definition of $ψ^*g$, we have $g_{ψ( p)} (dψ_p X_p , dψ_p Y_p ) = 0$ for all $Y_p ∈ T_p N$, In particular taking $X_p = Y_p$ and using the fact that $g$ is positive definite we see that $dψ_p X_p = 0_{\\psi(p)}$, and therefore $X_p = 0_p$ .\n",
"\n",
"
\n",
"\n",
"### Immersions and embeddings\n",
"\n",
"
\n",
"\n",
"A smooth mapping with maximal rank is called an **immersion** (in other words, $ψ : N → M$ is an immersion if for all $p ∈ N$, the rank of the linear mapping $dψ_p$\n",
"is equal to the dimension of $N$ ).
\n",
"The smooth map $\\ \\psi: N\\to M\\ $ is an **embedding** if it is one-to one immersion and the image $\\psi(N)$ with the subspace topology is homeomorphic to $N$ under $\\psi$.\n",
"\n",
"
\n",
"\n",
"**Example 22.4**\n",
"The metric on the standard **sphere** $S^2$:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\\(\\displaystyle {\\psi}^*g = \\mathrm{d} {\\theta}\\otimes \\mathrm{d} {\\theta} + \\sin\\left({\\theta}\\right)^{2} \\mathrm{d} {\\phi}\\otimes \\mathrm{d} {\\phi}\\)"
],
"text/latex": [
"$\\displaystyle {\\psi}^*g = \\mathrm{d} {\\theta}\\otimes \\mathrm{d} {\\theta} + \\sin\\left({\\theta}\\right)^{2} \\mathrm{d} {\\phi}\\otimes \\mathrm{d} {\\phi}$"
],
"text/plain": [
"psi^*(g) = dth⊗dth + sin(th)^2 dph⊗dph"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"%display latex\n",
"M = Manifold(3, 'R^3') # manifold M=R^3\n",
"c_xyz. = M.chart() # Cartesian coordinates\n",
"N = Manifold(2, 'N') # manifold N=S^2\n",
" # chart on N:\n",
"c_sph.=N.chart(r'th:(0,pi):\\theta ph:(0,2*pi):\\phi') \n",
"psi = N.diff_map(M, (sin(theta)*cos(phi), \n",
" sin(theta)*sin(phi),cos(theta)),name='psi',\n",
" latex_name=r'\\psi') # embedding S^2 -> R^3\n",
"g=M.metric('g') # standard metric on R^3\n",
"g[0,0],g[1,1],g[2,2]=1,1,1 # components of g\n",
"plb=psi.pullback(g) # pullback of g\n",
"plb.display() # show metric on S^2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It is easy to check that $\\psi$ is of maximal rank:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\\(\\displaystyle 2\\)"
],
"text/latex": [
"$\\displaystyle 2$"
],
"text/plain": [
"2"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"psi.jacobian_matrix().rank()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Example 22.5**\n",
"\n",
"The metric on the **paraboloid**:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\\(\\displaystyle {\\psi}^*g = \\left( 4 \\, u^{2} + 1 \\right) \\mathrm{d} u\\otimes \\mathrm{d} u + 4 \\, u v \\mathrm{d} u\\otimes \\mathrm{d} v + 4 \\, u v \\mathrm{d} v\\otimes \\mathrm{d} u + \\left( 4 \\, v^{2} + 1 \\right) \\mathrm{d} v\\otimes \\mathrm{d} v\\)"
],
"text/latex": [
"$\\displaystyle {\\psi}^*g = \\left( 4 \\, u^{2} + 1 \\right) \\mathrm{d} u\\otimes \\mathrm{d} u + 4 \\, u v \\mathrm{d} u\\otimes \\mathrm{d} v + 4 \\, u v \\mathrm{d} v\\otimes \\mathrm{d} u + \\left( 4 \\, v^{2} + 1 \\right) \\mathrm{d} v\\otimes \\mathrm{d} v$"
],
"text/plain": [
"psi^*(g) = (4*u^2 + 1) du⊗du + 4*u*v du⊗dv + 4*u*v dv⊗du + (4*v^2 + 1) dv⊗dv"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M = Manifold(3, 'R3') # manifold M=R^3\n",
"c_xyz. = M.chart() # Cartesian coordinates\n",
"N=Manifold(2,name='R2') # manifold N; paraboloid\n",
"c_uv.=N.chart() # coordinates on N\n",
"g = M.metric('g'); # metric in R^3\n",
"g[:]=[[1,0,0],[0,1,0],[0,0,1]]; # components of g\n",
"psi = N.diff_map(M, (u, v,u^2+v^2),\n",
" name='psi',latex_name=r'\\psi') # embedding N->R^3\n",
"plb=psi.pullback(g) # pullback of g\n",
"plb.display() # show metric on paraboloid"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\\(\\displaystyle 2\\)"
],
"text/latex": [
"$\\displaystyle 2$"
],
"text/plain": [
"2"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# check if psi is of maximal rank\n",
"J = psi.jacobian_matrix()\n",
"J.rank()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Example 22.6**\n",
"\n",
"The metric on the **hyperboloid**:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\\(\\displaystyle {\\psi}^*g = \\left( 4 \\, u^{2} + 1 \\right) \\mathrm{d} u\\otimes \\mathrm{d} u -4 \\, u v \\mathrm{d} u\\otimes \\mathrm{d} v -4 \\, u v \\mathrm{d} v\\otimes \\mathrm{d} u + \\left( 4 \\, v^{2} + 1 \\right) \\mathrm{d} v\\otimes \\mathrm{d} v\\)"
],
"text/latex": [
"$\\displaystyle {\\psi}^*g = \\left( 4 \\, u^{2} + 1 \\right) \\mathrm{d} u\\otimes \\mathrm{d} u -4 \\, u v \\mathrm{d} u\\otimes \\mathrm{d} v -4 \\, u v \\mathrm{d} v\\otimes \\mathrm{d} u + \\left( 4 \\, v^{2} + 1 \\right) \\mathrm{d} v\\otimes \\mathrm{d} v$"
],
"text/plain": [
"psi^*(g) = (4*u^2 + 1) du⊗du - 4*u*v du⊗dv - 4*u*v dv⊗du + (4*v^2 + 1) dv⊗dv"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M = Manifold(3, 'R3') # manifold M=R^3\n",
"c_xyz. = M.chart() # Cartesian coordinates\n",
"N=Manifold(2,name='R2') # manifold N; hyperboloid\n",
"c_uv.=N.chart() # coordinates on N\n",
"g = M.metric('g'); # metric in R^3\n",
"g[:]=[[1,0,0],[0,1,0],[0,0,1]]; # components of g\n",
"psi = N.diff_map(M, (u, v,u^2-v^2),\n",
" name='psi',latex_name=r'\\psi') # embedding N->R^3\n",
"plb=psi.pullback(g) # pullback of g\n",
"plb.display() # show metric on hyperboloid"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\\(\\displaystyle 2\\)"
],
"text/latex": [
"$\\displaystyle 2$"
],
"text/plain": [
"2"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# check if psi is of maximal rank\n",
"J=psi.jacobian_matrix()\n",
"J.rank()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
\n",
"\n",
"### Levi-Civita connection\n",
"\n",
"
\n",
"\n",
"If $M$ is a Riemannian manifold, then there exists a unique connection $∇$ with vanishing torsion and such\n",
"that $∇_X g = 0$ for all $X ∈ \\mathfrak{X}(M)$. \n",
"Connections with this properties are called **Riemannian or Levi-Civita connections.**\n",
"\n",
"From the definition (21.13) of torsion \n",
"($T (X, Y) = ∇_X Y − ∇_Y X − [X, Y]$) it follows that the torsion vanishes iff\n",
"\\begin{equation}\n",
"[X, Y] = ∇_X Y − ∇_Y X.\n",
"\\tag{22.5}\n",
"\\end{equation}\n",
"\n",
"Recall that connections satisfying (22.5) are called **torsion free** or **symmetric**.\n",
"\n",
"\n",
"From the rule of covariant differentiation of covariant tensor fields (21.14) it follows that\n",
"\n",
"$$(∇_X g)(Y,Z)=X (g (Y,Z))-g(∇_X Y,Z)-g(Y,∇_X Z),$$\n",
"\n",
"so $∇_X g$ vanishes iff\n",
"\n",
"\\begin{equation}\n",
"X (g(Y, Z)) = g(∇_X Y, Z) + g(Y, ∇_X Z).\n",
"\\tag{22.6}\n",
"\\end{equation}\n",
"\n",
"A connection $\\ \\nabla\\ $ in a Riemannian manifold $M$ is said to be \n",
"**compatible with the metric** if (22.6) is true for all \n",
"$X,Y,Z\\in \\mathfrak{X}(M)$.\n",
"\n",
"
\n",
"\n",
"One can prove that:\n",
"\n",
"**On any Riemannian manifold $M$ there exists a unique connection for which (22.5),(22.6) hold for $X, Y, Z ∈ \\mathfrak{X}(M).$**\n",
"\n",
"
\n",
"\n",
"Let us assume that (22.5),(22.6) hold true. We have \n",
"\n",
"$$\n",
"X( g(Y, Z)) +Y( g(Z, X)) − Z (g(X, Y))\\\\\n",
"= g(∇_X Y, Z) + g(Y, ∇_X Z)\n",
"+ g(∇_Y Z, X) + g(Z, ∇_Y X)\n",
"− g(∇_Z X, Y) − g(X, ∇_Z Y)\\\\\n",
"= g(∇_X Y + ∇_Y X, Z) + g(Y, ∇_X Z − ∇_Z X) + g(X, ∇_Y Z − ∇_Z Y)\\\\\n",
"= g(∇_X Y + ∇_X Y+(\\nabla_Y X-\\nabla_X Y), Z) + g(Y, ∇_X Z − ∇_Z X) + g(X, ∇_Y Z − ∇_Z Y)\\\\\n",
"= g(∇_X Y + ∇_X Y + [Y, X], Z) + g(Y, [X, Z]) + g(X, [Y, Z])\\\\\n",
"= 2g(∇_X Y, Z) + g(Z, [Y, X]) + g(Y, [X, Z]) + g(X, [Y, Z]).\n",
"$$\n",
"\n",
"The obtained equality implies the following **Koszul formula**\n",
"\n",
"\\begin{equation}\n",
"2g(∇_X Y, Z) = X (g(Y, Z)) + Y (g(Z, X)) − Z (g(X, Y))\\\\\n",
"− g (Z, [Y, X]) − g (Y, [X, Z]) − g( X, [Y, Z]).\n",
"\\tag{22.7}\n",
"\\end{equation}\n",
"\n",
"Since $g$ is non-singular if such a $g$ exists, it defines the connection $∇_X Y$ in an unique manner."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let us assume that $\\nabla_X Y$ is defined by (22.7). We have\n",
"\n",
"$$2g(\\nabla_XY,Z)-2g(\\nabla_Y X,Z)=\\\\\n",
"X(g(Y,Z))+Y(g(Z,X))-Z(g(X,Y))\n",
"-g(Z,[Y,X])-g(Y,[X,Z])-g(X,[Y,Z])\\\\\n",
"-Y(g(X,Z))-X(g(Z,Y))+Z(g(Y,X))\n",
"+g(Z,[X,Y])+g(X,[Y,Z])+g(Y,[X,Z])\\\\\n",
"=2g([X,Y],Z),\n",
"$$\n",
"thus (22.5) is fulfilled.\n",
"\n",
"To check (22.6) let us note that\n",
"\n",
"$$2g(\\nabla_XY,Z)+2g(Y,\\nabla_XZ)=2g(\\nabla_XY,Z)+2g(\\nabla_XZ,Y)\\\\\n",
"=X(g(Y,Z))+Y(g(Z,X))-Z(g(X,Y))-g(Z,[Y,X])-g(Y,[X,Z])-g(X,[Y,Z])\\\\\n",
"+X(g(Z,Y))+Z(g(Y,X))-Y(g(X,Z))-g(Y,[Z,X])-g(Z,[X,Y])-g(X,[Z,Y])\\\\\n",
"=2X(g(Y,Z)).\n",
"$$\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
\n",
"\n",
"### Levi-Civita connection in components\n",
"\n",
"
\n",
"\n",
"\n",
"Let $g_{ij}=g\\big(\\frac{\\partial}{\\partial x^i},\\frac{\\partial}{\\partial x^j}\\big)$. Recall from the notebook 12 that $[\\frac{\\partial}{\\partial x^i},\\frac{\\partial}{\\partial x^j}]=0.$\n",
"\n",
"If we put in (22.7) .\n",
"$$X=\\frac{\\partial}{\\partial x^i},\\quad\n",
"Y=\\frac{\\partial}{\\partial x^j},\\quad\n",
"Z=\\frac{\\partial}{\\partial x^k},\n",
"$$\n",
"then \n",
"\n",
"$$2g\\Big(\\nabla_{\\frac{\\partial}{\\partial x^i}}\n",
"\\frac{\\partial}{\\partial x^j},\\frac{\\partial}{\\partial x^k}\\Big)\\\\\n",
"=\\frac{\\partial}{\\partial x^i}\n",
"\\big(g(\\frac{\\partial}{\\partial x^j},\\frac{\\partial}{\\partial x^k}\\big)\\big)\n",
"+\n",
"\\frac{\\partial}{\\partial x^j}\n",
"\\big(g(\\frac{\\partial}{\\partial x^k},\\frac{\\partial}{\\partial x^i}\\big)\\big)\n",
"-\n",
"\\frac{\\partial}{\\partial x^k}\n",
"\\big(g(\\frac{\\partial}{\\partial x^i},\\frac{\\partial}{\\partial x^j}\\big)\\big)\\\\\n",
"=\\frac{\\partial g_{jk}}{\\partial x^i}+\n",
"\\frac{\\partial g_{ki}}{\\partial x^j}-\n",
"\\frac{\\partial g_{ij}}{\\partial x^k}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using the formula \n",
"$\\nabla_{\\frac{\\partial}{\\partial x^i}}\n",
"\\frac{\\partial}{\\partial x^j}=\\Gamma^m_{ji}\\frac{\\partial}{\\partial x^m}\n",
"$ we obtain\n",
"\n",
"$$2\\Gamma^m_{ji}g_{mk}=\n",
"\\frac{\\partial g_{jk}}{\\partial x^i}+\n",
"\\frac{\\partial g_{ki}}{\\partial x^j}-\n",
"\\frac{\\partial g_{ij}}{\\partial x^k},\n",
"$$\n",
"and finally\n",
"\\begin{equation}\n",
"\\Gamma^m_{ji}=\\frac{1}{2}g^{km}\\Big(\n",
"\\frac{\\partial g_{jk}}{\\partial x^i}+\n",
"\\frac{\\partial g_{ki}}{\\partial x^j}-\n",
"\\frac{\\partial g_{ij}}{\\partial x^k}\n",
"\\Big),\n",
"\\tag{22.8}\n",
"\\end{equation}\n",
"\n",
"where $[g^{km}]$ is the matrix inverse to $[g_{mk}]$.\n",
"\n",
"Note that from (22.8) it follows that $\\Gamma_{ij}^k=\\Gamma_{ji}^m$.\n",
"\n",
"Note that a symmetric $n\\times n$ matrix has $n(n+1)/2$ independent elements and $\\Gamma^m_{ij}$ defines $n$ such matrices, so the Riemannian manifolds have $n^2(n+1)/2$ independent Christoffel symbols."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
\n",
"\n",
"### Geodesics in Riemannian manifolds\n",
"\n",
"
\n",
"\n",
"\n",
"In the case of Riemannian connections SageMath offers `integrated_geodesic`, a numerical method of finding geodesics \n",
"
\n",
"\n",
"**Example 22.7**\n",
"\n",
"Use `integrated_geodesic` method to find the geodesics on Poincaré half-plane passing through the point with coordinates $(0,1)$.\n",
"\n",
"\n",
"Let us start from the geodesic with tangent vector at $(0,1)$ parallel to $x$-axis."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAADbCAYAAACcGb3IAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAqLElEQVR4nO3deXgT1frA8e9pS0stu6xSWS5rBa+CCAIKIqAiLlxQQRRFRanWBREVBAT03osoKihFRVkUEfiJoICAIrIUxQUuq2WTHSkgi2Utpe35/XGaJi0tadI0M0nez/PM08l0knk7TftmzpxzXqW1RgghhLiYMKsDEEIIYX+SLIQQQrglyUIIIYRbkiyEEEK4JclCCCGEW5IshBBCuCXJQgghhFuSLIQQQrglyUIIIYRbkiyEEEK4JclCCCGEW5IshBBCuCXJQgghhFuSLIQQQrglyUIIIYRbkiyEEEK4FWF1AEL4m1KqAvAGcCmQAfTUWp93+f4ooK7WuptFIQphO0oq5YlQo5QaB/wbkyw2AXdqredlf08BR4HVWuubrYtSCHuRZigRUpRScUCK1vog0Cp782GXXa4EygM/+Ds2IexMkoUINZWAqdnrDwJ/AL+6fL9N9ldJFkK4kHsWIqRorVcAKKVqANcDQ3Tuttg2wAlgjQXhCWFbcmUhQlXX7K+z8mxvA6zQWmf6OR4hbE2ShQhV12LuXWx1bFBKNQCqIE1QQlxAkoUIVZcCe/Js65D9damfYxHC9iRZiFC1GqihlAoHUEr9E3gNOAastzIwIexIbnCLUPVfIBZYoJT6AzgFlAAWahl8JMQFZFCeCDnZA+9Kaq3PumzrAswB7tBaz7cqNiHsSpKFCDlKqW+BlkA1rfXp7OTxI3BCa32rtdEJYU9yz0KEomsxA/HOZt+zeAfzt9Dd0qiEsDG5shAhRynVEegIXAJUxiSOsa6TCQohcpNkIYQQwi1phhJCCOGWJAshhBBuSbIQQgjhliQLIYQQbkmyEEII4ZYkCyGEEG5JshBCCOGWJAshhBBuSbIQQgjhliQLIYQQbkmyELahlGqjlJqnlDqglNLZ04a7e05bpdQapVSaUmqnUireD6EKEXIkWQg7icFUqXuqMDsrpWoDC4AkoAmmoNG7SqlunhxUGWWypyoXQuTDn5XyZMZCcVGuk1oqpZgzZ86ci+3/4osvMnfuXDZv3vws8CxAfHw869evn+XJcVNTUylbtiypqanehC1EcbPFhxgpqyoC1qpVq7i1fXvYvx8OHoSDB3k8IoLpv/5K5qBBhJ84AcePw+nTZjlzBs6cIevsWTh/HjIzUX/9hU5LMy/YqBFERkJEhPl6ySVmiYkxX8uWhQoVoHx5s1SsCFWrmqVaNbOPEEHKn1OUy5WFKLTsKwu6dO4Mu3fDtm2wYwfs2WMe79nD0XXruDQzs8jHOgGUBVKBMkV5odKl4fLLoWZN5/KPf0D9+lCvnkk6QnjOFlcWkiyEPWRkmISwcSNs3MjX//kPN1WvTulDh8z3fCUqCl2yJDo8HMLDUSdPciI9nXJZWaReeilllDLHS083VyK+VL26SRyNG8OVV5qlcWMoVcq3xxHBRpKFCFEZGZCcDKtXO5f1680/6MIKC+NwRARnypenVsuWOc1B6w4c4I2PP+bT+fOJqFTJNBeVKmU+1UdHQ3j4BS914sSJnHsWZcq4XFtoDWlpJmmcOgV//22atY4fh2PH4K+/4NAhSEkxzWAHDsC+feY5nqhXD5o1cy5Nm0oCEa4kWYgQcfIk/PwzrFwJP/5o1k+fdvu0jMhIIho2NJ/GGzSAunWhVi3TvBMby0tDhjBv3jySk5NznvPEE0+wbt06Vq1aVejwCkwW3tAaDh92Npft2GGumBzLkSPuXyMsDK66Cq6/Hlq3NktsbNHiEoFMkoUIUmlp8NNPsGQJ/PAD/PYbuLu3UL8+5xs35kjVqpytW5cO/fvz9OjRtGvfngoVKlCjRg0GDRrEn3/+yaeffgrArl27aNy4MX379uWxxx5j1apVxMfHM336dLp1c997NjExkcTERDIzM9m2bZtvkoU7R47kNLWxcaO5otqwAc6du/jz6taFm26C9u2hXTuoVKl44xR2IslCBJE//oAFC8yyfPnFm2KqV4dWraB5c2ezS5kyLFu2jHbt2l2w+0MPPcSUKVPo3bs3u3fvZtmyZTnfW758Oc899xy///47l112GS+99BLx8Z6Ny/PplYU3zp+H3383zXG//WauvDZuNFcpBWnSBG67DTp1ghYtTA8uEawkWYgAlpkJq1bBnDkwbx5s317wvnFx0LataVa5/nqoUQNsNP7N8mSRn7//djbdrVhh1s+fz3/f8uVN0ujSxXyV+x3BxhZ/LJIsROFlZJhmpS++gLlzTdt8fmJjoWNH02Ry001mDIINWdIM5a3Tp839niVL4Pvv4X//y3+/qChz7rt2hX/9C8qV82uYolhIshABICvLfLqdMQNmzTI9gPIKDzdXDJ06maaRxo1tdeXgji2vLNw5eBAWLTLNft99B/mNPo+MhFtvhR494I475IojcNnij0mShcjf9u3wyScwdSrs3Xvh96Oj4ZZbzKfX2283I5sDVEAmC1fp6eY+0Zw58NVXpitvXjExcPfd8NBDpkkwTKaFCyCSLITNnD4NM2fCxImmN1NeJUuaT6jdu5uriACf3iKgmqEKKysLfvnF/B7/7//yTxw1a5qk8eij5v6RsDtJFsImNmyADz+Ezz6DEydyfy8szFxB3H8/3HmnmdIiyAT8lUVBMjMhKQk+/9wkjrxNVWFhJun37WuaD/MZsChsQZKFsFBGhmmyGDvW3JPIq1Ej6N3bJAmb3qD2laBNFq7OnjWdEqZMMfc4srJyf79GDUhIgD59ArpJMUjZIllIw2WoOX4cXn/dTHB3zz25E0V0NDz8sLOf/4ABliSK8ePHU7t2bUqWLMk111xDUlLSRfdPTEwkLi6O6OhoGjRokDNoT7iIjjbNhwsXmtHlw4fnHhW+dy+89JLZ1rcvbN5sWajCprTW/lqElfbs0fq557SOidHaDPdyLnFxWr/3ntbHj1sdpZ4xY4YuUaKE/uijj3RycrJ+9tlndUxMjN6zZ0+++48fP16XLl1az5gxQ+/YsUNPnz5dlypVSs+dO9ftscaNG6fj4uJ0/fr1NaBTU1N9/ePY2/nzWs+bp3WnThe+J0DrO+7QOilJ66wsqyMNdf78P13gIski2G3erHWvXlpHROT+R6CU1rffrvV339nqn0Hz5s11fHx8rm0NGzbUAwcOzHf/li1b6gEDBuTa9uyzz+rWrVsX+pipqamhmSxcbd2q9dNPa12q1IVJo2VLk1Rs9D4JMZYnCq21NEMFrU2bTP/6K64w3V8d03yXLAlPPAFbt5qR1x072mZMRHp6OmvWrOHmm2/Otf3mm2/mp/x6ZwHnzp2jZMmSubZFR0fz66+/cr6AEc/nzp3jxIkTuZaQV78+vPsu/PknjB6du4lq1SrTC+6aa0z33Lz3O0RIkGQRbDZvNvcirrzSdJ/U2f0KypeHoUNNe/X48WZabJs5cuQImZmZVKlSJdf2KlWqcPDgwXyfc8stt/Dxxx+zZs0atNasXr2aSZMmcf78eY4UMMPryJEjKVu2bM5y+eWX+/xnCVhlysDzz5vZcj/5xHR0cFi71owMb9LE3CzX0mcllEiyCBa7dpneS40bm5HWDpUqwahRJkm8+ipUrmxZiIWl8lzpaK0v2OYwdOhQOnXqxHXXXUeJEiW466676N27NwDhBXQFHTRoEKmpqTnLvn37fBp/UIiMhAcfNN2q58wxCcJhwwa46y5o2dJMPyJCgiSLQHf0KPTrZ+o9fPKJs4mgShV4+22TRF58MSDGR1SsWJHw8PALriIOHz58wdWGQ3R0NJMmTeLMmTPs3r2bvXv3UqtWLUqXLk3FihXzfU5UVBRlypTJtYgChIWZCQrXrIH58+Haa53f++UX6NDBNGWuW2dVhMJPJFkEqrQ0eOMNqFPHjJVwtM+XLw8jR5pmhOeeC6i6z5GRkVxzzTUsXrw41/bFixfTqlWriz63RIkSxMbGEh4ezowZM7j99tsJkyktfEcp6NzZJIg5c8wVrMP335tp5h96yFQKFMHJj3fThS9kZWk9e7bWtWrl7rESHa314MG26P5aFI6usxMnTtTJycm6X79+OiYmRu/evVtrrfXAgQN1r169cvbfunWrnjp1qt62bZv+5ZdfdPfu3XWFChX0rl273B4r5LvOFkVGhtaffZb/+3D4cK3PnLE6wmBieU8oraXrbGDZtEnr9u31BV1gH3lE6/37rY7OZxITE3XNmjV1ZGSkbtq0qV6+fHnO9x566CHdtm3bnMfJycn66quv1tHR0bpMmTL6rrvu0lu2bPHoeNJ1tgjS0rR+6y2ty5XL/b6sUUPrL76Q7ra+YXmi0FrLdB8B4dQpGDEC3nknd3nS9u3NfYl//tO62IJASEz3UdyOHYPXXoNx45zdtMHc07Bp77sAYou+7dKoa2dam/bhuDjT992RKGrVgtmzYfFiSRRFkJiYyBVXXMG1rjdthXcqVDAfZjZsANdxMt9/b+5vDBt28VK7wvbkysKuDhyAJ5+Er792bouKgpdfhhdeMHP9CJ+QKwsf09qMw3jmmdy1UOrVg48/hjZtrIstMMmVhciH1vDRR2bktWuiuPVW+P13eOUVSRTC3pQy4zCSk83khBERZvv27abw0hNPXDgVvrA9SRZ2snev6bP++OPO2gOVK5uR2AsWmG6ywmekGaqYxcSYGY7XrTMD+Bw++MCMDM/TRVrYmzRD2YHWMHmyGVx38qRze+/e8NZbUl+gmEkzlB9kZkJiomlGPX3auf2JJ8x4IakPfjHSDCWAv/4yl+yPPupMFLGxsGiRSSAhmCg8rWcxbdo0rrrqKi655BKqVavGww8/zNGjR/0UrSiU8HBzD2PTJtOLz+H99+Gqq0wNFWFvfuynK/JauFDrKlVy909/+GGt//7b6sgs42k9i6SkJB0WFqbHjh2rd+7cqZOSknSjRo10ly5dCn1MGWfhZ5mZWicman3JJc73fXi41iNGmBobIi/Lx1hoLYPyrHHunClE5JokKlXSuhAFe4Kdp/Us3nzzTf2Pf/wj17Z3331Xx8bGuj2WjOC22PbtplaG699B69Za791rdWR2Y3mi0FrqWfjfzp3QurXpk+7QqZMpY3rHHdbFZQPe1LNo1aoV+/fvZ8GCBWitOXToELNmzaJz585uj5eQkEBycjK//fabT+IXHqpbF1asMCVeHfN4/fgjXH21mbRQ2IokC39yTPW8erV5HBlpJgH85hszS2yI86aeRatWrZg2bRrdu3cnMjKSqlWrUq5cOd57770CjyPFj2wkIsIM2EtKgpo1zbZjx8wHpwEDnBNkCstJsvCHjAwzTXjXrs7+5XXqmApkzzxjm0p1duFJPYvk5GSeeeYZXnnlFdasWcOiRYvYtWsX8fHxBb6+FD+yoVatTHGlLl2c2956y9wML+CDgvAzP7Z5haZDh7S+8cbc7bLdu2st7eMXOHfunA4PD9ezZ8/Otf2ZZ57Rbdq0yfc5DzzwgL777rtzbUtKStKAPnDgQL7PSUtL06mpqTnLvn375J6FXWRlaT12rNYlSjj/XqpW1TopyerIrGT5/Qot9yyK2dq10KwZLFtmHkdEmGan6dNN+UqRizf1LM6cOXNB3QpHhTyt8x/aI8WPbEwpc7W9YgVUr262HTwIN90EEyZYG1uo82NmCi0zZpi5/R2fjqpV03rlSqujsj1P61lMnjxZR0RE6PHjx+sdO3bolStX6mbNmunmzZu7PZb0hrK5Q4e0btcu91X5E09onZ5udWT+ZvlVhdbSddb3srK0HjYs9xu8RQut//zT6sgChif1LLQ2XWWvuOIKHR0dratVq6bvv/9+vd+D+h4yzsLGzp/Xul+/3H9PN92k9bFjVkfmT5YnCq2lnoVvpaXBI4+YZiaH3r3NKNWSJS0LS1ycTPcRAKZMgb59IT3dPG7QwHSvrVvX0rD8xBY9YOSeha8cPWp6bjgShVKmBsWkSZIohCiq3r3Nvb9KlczjrVvhuutMj0LhF5IsfGHXLtP1zzFwLCYGvvoKnn9eusXamMw6G2BatoRffjHT94P5gHbTTeZvTRQ7aYYqqrVrzQjsQ4fM46pVzSC7pk2tjUsUmjRDBZjUVOjWDZYsMY+VMuVcn3zS2riKjy0+ccqVRVEsXw433uhMFA0bmtkzJVEIUXzKljX1XR54wDzWGhIS4NVXzbooFpIsvDV3Ltxyi3NEdqtWZl4bx5QFwvakGSqARUbCp5+aSnwOw4aZmjBZWZaFFcykGcobM2aYTzWZmebxbbfBF1/AJZdYG5fwijRDBbi33jLzSDk8/LApTZw9ODMISDNUQJo8GXr2dCaKnj3NDTZJFD7jSfGj3r17o5S6YGnUqJEfIxaWev5583fpGMk/eTLcf79MQuhjkiw8MWGCGUfhuBp7/HGYOhVKlLA2riAyc+ZM+vXrx+DBg1m7di033HADnTp1Yu/evfnuP3bsWFJSUnKWffv2UaFCBe655x63x5JmqCDSu7epVR8RYR7PnAk9ekjC8CU/jgAMbB98kHsUab9+ZrS28ClPix/lNWfOHK2UypkepDBkBHcQmT9f66go599p167BMD2I5aO3tZaJBAvno4/AdcrrF16At9+WMRQ+5k3xo7wmTpxIhw4dqCkdDUJT586m84ljIOzs2eYKIyPD2riCgCQLdz791Ewz4PDiizBqlCSKYuBN8SNXKSkpLFy4kD59+lx0Pyl+FORuvvnChPHgg877jMIrkiwu5osvTM8Kxz2K55+H11+XRFHMPCl+5GrKlCmUK1eOLq4FdPIhxY9CQMeOJmFERZnH06fDY49Jt9oikGRRkEWLTI8Kx5vrqafgzTclURSjihUrEh4efsFVxOHDhy+42shLa82kSZPo1asXkZGRF9130KBBpKam5iz79u0rcuzChjp2hFmznDe9J082H/hk4J5XJFnk56efTAlUR0+KRx4xRYskURQrb4ofOSxfvpw//viDRx991O1xpPhRCLn9djMuytGtdswY+Pe/LQ0pYPnxbnpg2LRJ63LlnL0punXTOiPD6qhChqfFjxweeOAB3aJFC4+OJcWPQsjEibl7M374odURecLynlBaS/Gj3Pbt0zo21vmG6tBB67Q0q6MKOZ4WP/r77791dHS0njBhglfHk66zIWL0aOffdliY1nPmWB1RYVmeKLSW4kdOqalwww2wcaN5fM01sHQplC5tbVyi2Ml0HyHkhRdMnRkwvaWWLjV1MezNFu3fcs8CzL2Je+5xJoo6dcw045IogpqM4A5Bo0aZjitgKlveeSfs3GltTAFCriy0NuMoPvrIPK5QwUwzXq+etXEJv5ErixCTng633mquKsCUaP35ZyhXztKwLkKuLGxh7FhnooiMhK+/lkQhRDCLjIQvv4S4OPN461a4914Z5e1GaCeLhQtNv2uHSZPg+uuti0cI4R/ly5um5ooVzePFi+G556yNyeZCN1ls2wb33eccdDd4sLMtU4QEuWcR4mrXNlOBOGaNHjcOJk60NiYbC817FidOmB4Qmzebx//6lxnpGRa6uTOUyT2LEDdpEjgGc0ZGmnLJ9uohJfcsLKG1me/JkSgaNYJPPpFEYSOeFD8CMzHg4MGDqVmzJlFRUdSpU4dJkyb5KVoR8B55xEznA+bmd7ducOiQtTHZUOj9hxw92lx6gun98NVX0kXWRjwtfgRw7733smTJEiZOnMjWrVuZPn06DRs2dHssaYYSOd5+G9q0MesHDsi05vkIrWaopCRo1845VfH8+Wb+e2EbLVq0oGnTprz//vs52+Li4ujSpQsjR468YP9FixbRo0cPdu7cSYUKFbw6pjRDCcBcTTRtapIFwKBB8N//WhuTIc1QfvXXX+bTgiNRDBkiicJmvCl+NHfuXJo1a8Ybb7xB9erVqV+/PgMGDODs2bMFHkfqWYh8ValiyhI4ZqkdOdL0mBRAqCSLrCxT/MTxiaFdOxg+3NKQxIW8KX60c+dOVq5cyaZNm5gzZw5jxoxh1qxZJCQkFHgcqWchCtSqlalZ49CrF/z5p3Xx2EhoJIt33jH1KcB8evj8cwgPtzYmUSBPih9lZWWhlGLatGk0b96c2267jbfffpspU6YUeHUh9SzERfXvD3fcYdaPHjUJQ6rshUCy+N//TNujw9SpULWqdfGIAnlT/KhatWpUr16dsmXL5myLi4tDa83+/fvzfY7UsxAXpZQplFS9unm8dKmZUyrEBXeyOHMGevZ0FjF68UVTPUvYkjfFj1q3bs2BAwc4depUzrZt27YRFhZGbGxsscYrgtill8K0ac4u9cOGwerV1sZkNT/Oh+5/CQnO+eubNtX63DlLwhCF52nxo5MnT+rY2Fh99913699//10vX75c16tXT/fp08ftsaT4kXDr5Zed/0Pq19f69GkrorC8loXWwVz86LvvnL/k6GitN2/2ewjCO54WP9q8ebPu0KGDjo6O1rGxsbp///76zJkzhT6eFD8SBUpP17pZM+f/kmeesSIKyxOF1sFa/Cg1FRo3Bkeb9XvvOUdoCpGHjLMQF7VlCzRpYupfAPzwg+lR6T8yzqLYPP+8M1F06ABPPmltPMKWZAS3KJSGDXN3p330UXC5RxYqgu/KYvFicAzqKl0aNm2CGjX8cmgRmOTKQriVlWWuJlasMI+ffhrefddfR5crC587fRoef9z5+M03JVEIIYouLMxMXx4dbR6PGwerVlkbk58FV7J45RXYvdus33hj7sQhRB7SDCU8UrcuvPaaWdca+vQxs9SGiOBphlq7Fpo1M5eLUVGwcaOURxWFIs1QotAyMkytizVrzOP//Adefrm4jyrNUD6TmQnx8c6qd0OHSqIIYJ7Us1i2bBlKqQuWLVu2+DFiETIiIuCjj5yD9V57DXbutDYmPwmOZDFxIvz6q1mPi4MXXrA2HuE1b+pZAGzdupWUlJScpZ58WBDFpUkTePZZs56W5lwPcoHfDHX0KNSvD8eOmcdLl5r7FSIgeVrPYtmyZbRr147jx49Trlw5j46VmJhIYmIimZmZbNu2TZqhROGdPGm61Dpmsp471zn5oO9JM5RPDB3qTBQ9e0qiCGDe1LNwaNKkCdWqVaN9+/YsXbq0UMdLSEggOTmZ3377zeuYRYgqXdpU13N47jk4d866ePwgsJPFhg3w4YdmPSbGdJUVAcubehbVqlVjwoQJfPnll8yePZsGDRrQvn17Vjj6w+dDih8Jn7j3XueH0x07TCmEIBZhdQBe09rMO++4qT1kCFx2mbUxCZ/wpJ5FgwYNaNCgQc7jli1bsm/fPkaPHk0bR03lPEaOHMmIESN8F7AITUrB2LHmHkZWlukZ9fDDpmZOEArcK4tvvoElS8x67drQr5+l4Yii86aeRX6uu+46tm/fXuD3pfiR8Jl//tM5nuvUKdMsHqQCM1lkZJjaFA6jRkHJktbFI3zCm3oW+Vm7di3VqlUr8PtS/Ej41IgR5h4GmJ6ZycnWxlNMArMZavJk2LzZrLdsCXffbW08wmf69+9Pr169aNasGS1btmTChAns3buX+Ph4wFwV/Pnnn3z66acAjBkzhlq1atGoUSPS09P57LPP+PLLL/nyyy/dHsu1N5QQXqtc2VTjfPll0xw1cKDpHRVs/Dgfum+cOaP1ZZc555dfudJnLy3swZN6FqNGjdJ16tTRJUuW1OXLl9fXX3+9/uabbzw6ntSzEEV2+rTW1asX1/8ly2tZ6ICsZzF6tHPQ3Z13wtdf++RlReiS6T6ET3z8MTz2mFlv0waWLTM3wYvOFuMsAitZnDxpbmYfPWp+CRs2mCJHQhSBJAvhExkZ5v/R1q3m8XffQceOvnhlWySLwLrB/d57JlEA3HefJApRJDLrrPCpiAgYPtz5eOhQ0ygVJALnyuLECXNVceyYmcRr82YzzYcQRSRXFsJnsrLgqqtM0TWARYvglluK+qpyZeGRxETntB4PPCCJQghhP2FhMGyY8/Hw4UFzdREYVxanT0OtWnDkiFxVCJ+RiQRFscjKMoP1fv/dPP7+e2jfviivKFcWhTZxokkUAN27S6IQPiETCYpiERYGgwc7H7/+unWx+JD9k8X586a7rMOgQdbFIvzCk+JHrn788UciIiK4+uqrizdAIdy5916oU8esf/89rF5tbTw+YP9kMXMmOObu6dwZrrzS2nhEsfK2+FFqaioPPvgg7T243JfeUKLYhIfnnpIoCGbEtvc9C62haVNYt848Xr7cDHYRQcvT4kcOPXr0oF69eoSHh/PVV1+xzvGeKQTpDSWKRVoa1KwJhw+bpqk//jA9Oj0n9yzcWrbMmSiuvRZuuMHKaEQx87b40eTJk9mxYwfDXHuhCGG1kiXh6afNelaWGScWwOydLMaMca737++rofPCprwpfrR9+3YGDhzItGnTiIgo3LyYUvxI+E18vHNG7I8/NrNQBCj7JoudO2HePLNevTp062ZtPMJvClv8KDMzk549ezJixAjqe9BDbuTIkZQtWzZnufzyy4scsxD5qlgRevUy6ydPwiefWBtPEdg3Wbz/vnMwy5NPQokS1sYjip2nxY9OnjzJ6tWreeqpp4iIiCAiIoJXX32V9evXExERwQ8//JDvcaT4kfArR1MUmMHFATpIz57J4uxZmDTJrEdGOmdyFEHN0+JHZcqUYePGjaxbty5niY+Pp0GDBqxbt44WLVrkexwpfiT86soroW1bs75lCxTwIcbu7Fn8aNYs59Qe99wDlSpZG4/wG0+KH4WFhdE4z2SSlStXpmTJkhdsz48UPxJ+8+STpjcnwIcfFnVEtyXsmSwmTHCuZ/+TEKGhe/fuHD16lFdffZWUlBQaN27MggULqFmzJgApKSlux1wUVkJCAgkJCTldZ4UoNl26mIp6hw/DV1+Zr5UrWx2VR+w3zmLLFoiLM+txcWZ+FekFJYqRjLMQfjFwIIwaZdbffBMGDCjsM23xD9B+9ywmT3au9+kjiUIUGxnBLfzq0Ued6xMnBtyNbntdWWRkQI0akJJiCokcOCD3K0SxkysL4Tdt2oBjrrOff4YCOmHkYYtPzPa6svj+e5MoAG6/XRKFECK49O7tXA+wMRf2ShZTpzrXH3rIujiEEKI43HMPREeb9ZkzIT3d2ng8YJ9kcfq06SUAUL483HabpeGI4Cf3LITflS5tekaBGR7w7beWhuMJ+ySLuXPhzBmzfu+9ZjCeEMVIih8JS9x/v3N9+nTr4vCQfZLFzJnO9R49rItDWM6T4kcrV66kdevWXHrppURHR9OwYUPeeecdP0YrhIc6djStJwBff+38kGxz9kgWJ07AokVmvWpVmYo8hHla/CgmJoannnqKFStWsHnzZoYMGcKQIUOY4DqwswDSDCUsERkJXbua9TNnYOFCa+MpJHt0nZ0+HXr2NOsJCTBunJ9CEnbjbfEjV127diUmJoaprh0mLkK6zgq/++47uOUWs96jh7vmKOk6m2P2bOf63XdbF4ewlLfFj1ytXbuWn376ibaOidvyIfUshOXatXM2Rc2fD+fOWRtPIVifLNLSnJdhFSvC9ddbG4+wjDfFjxxiY2OJioqiWbNmJCQk0KdPnwL3lXoWwnIlSsAdd5j1U6cCYiZa65PF0qWm2yyYgXiFrHYmgldhix+5SkpKYvXq1XzwwQeMGTOG6Re5rJd6FsIWHF1owfQGtTnr/zO7nqQ777QuDmE5T4sfuapduzYAV155JYcOHWL48OHcd999+e4bFRVFVFSUb4IWwlsdO0JUlGmCmj8fxo+39Vx41l5ZaA3ffGPWIyPNyRMhy9PiRwXRWnMuANqARYgrVcrcuwDYvx/Wr7c2HjesvbJITgZHE8CNN5qTJ0KaJ8WPwHR/rVGjBg0bNgTMuIvRo0fztGspywJI8SNhuc6dncMGFi6Eq6+2NJyLsTZZOE4SwK23WheHsA1Pix9lZWUxaNAgdu3aRUREBHXq1OH111+nb9++bo8lxY+E5Tp1cq5/+y0MGmRdLG5YO86ic2dYsMCsJyc7ix4J4UcyzkJYql49+OMP00Pq+HGIicm7hy1uZFh7ZTFrFqxYAT/+CNnNCEL4izRDCVsYOhTCwsw92wsThW3488pCCFtSSpUBUoGyWmsZoSdEPiRZiJCnzCCO0sBJLX8QQuRLkoUQQgi3rB/BLYQQwvasH8EthBDC55RS5YFhmP/zdYH/Az4H3sT0sCoP/EdrnVyY15NkIYQQQUYpFQmMB57XWh9QStUEdgF3Af2AesA3wHHgqcK8pjRDCSFE8IkHJmutD2Q/TsNcTezWWu8CwoHtQKHrusqVhRBCBJ/jWuvvXB43y/66CEBrvRDwqESfXFkIIUSQ0VrnLRPZDsgEVnr7mtJ1VgghgpxS6n/Aea11C29fQ64shBAiiGX3iroKWJZne8HlJPMhyUIIIYKIUqqSUupXpdSw7E23Yv7X/+q6D1D4IjFIshBCiGDTFrgWM5NNNNAdOACUwmyMAd4Fhnvyon65Z+Ey944QQgjPFXreMqVUaeAdIB2TIEYCZYD/AnuASOANrfUGTwLwV7JwzOophBDCc5bPiOyvcRYngVhgf/bXk0V4rV+B5hY+3w6vURp7nEtfvIYdYvDF+bTDz+GL17DDufRFHHY4l0V9DV+dS5/wS7LQWmullOOHPVmUDKmUyrLy+XZ4DdOqB1h8Ln3xGjaJwbHq9fm0w8/hi9eww7n0URyWn8uivkaec2n5GIdAvMGdaPHz7fQaRWWHn8MOMfiCXX4OO/w+fMEOP4ddXsMW/DYoT6qR+Y6cS9+S8+k7ci59x27n0p9XFueAEdlfRdHIufQtOZ++I+fSd2x1LmW6DyGEEG4F4j0LIYQQfibJQgghhFuSLIQQQrglyUIIIYRbxZoslFKDlVI/KaXOKKX+LuRzlFJquFLqgFLqrFJqmVKqUXHGGQiUUuWVUlOVUqnZy1SlVDk3z5milNJ5lp/9FLJtKKWeVErtUkqlKaXWKKVucLN/2+z90pRSO5VS8f6KNRB4cj6VUjfm8x7USqmG/ozZjpRSbZRS87L/12mlVJdCPMey92ZxX1lEAl8A73vwnBeB/pgi4tcCB4HF2ZNjhbLPgasx0w3fmr2etxpWfhYB1VyW24onPHtSSnUHxgD/AZoAScBCpVSNAvavDSzI3q8JZvK1d5VS3fwSsM15ej5dNCD3+3B7MYYZKGKA9Zj/dW5Z/t7UWhf7AvQG/i7EfgpIAV5y2RYF/A309UesdlyAOEADLVy2XZe9rcFFnjcF+Mrq+C0+d78A7+fZthkYWcD+o4DNebZ9AKyy+mexw+LF+bwx+31azurY7bxkn6Mubvax9L1pt3sWtYGqQE6hca31OWA5HhbqCDItgVSt9S+ODVrrnzGjO92dlxuVUoeVUtuUUh8ppSoXZ6B2opSKBK7B5f2U7TsKPm8t89n/W6CZUqqEbyMMLF6eT4e1SqkUpdQSpVS7Ygkw+Fn63rRbsqia/fVQnu2HXL4XiqoCh/PZfpiLn5eFwP3ATcDzmGa9H5RSUT6P0J4qAuF49n6qWsD+EdmvF8q8OZ8pwONAN6ArsBVYopRqU1xBBjFL35sezzqrlBoODHOz27Va69VeRWTkHVau8tkW8Ap7LrO/5vfzX/S8aK1nujzcpJRajSl+0hmYXfhIA56n76f89s9ve6gq9PnUWm/FJAiHVUqpy4EBwIriCS+oWfbe9GaK8nHADDf77PbidcHczAaTQVNctlfmwowaDAp7Lv8JVMnne5Xw4LxorVOUUnuAeoV9ToA7AmRy4afei72fDhawfwZw1KfRBR5vzmd+fgYe8FVQIcTS96bHyUJrfQTzpikOuzAnpCOwFnLaSdsCLxXTMS1T2HOplFoFlFVKNdda/5q9rQVQFvipsMdTSl0KXE7uRBy0tNbpSqk1mPfTHJdvdQS+LuBpq4A78my7GVittT7v+ygDh5fnMz9NCJH3oI9Z+94s5jv8NTBdPF/BVHq6Onsp5bLPFuBfLo9fwvR++hfQGNNl9ABQ2uoeC1YumPsP6zG9oK4DNgDz8uyTcy4xtXdHY26K1cL0SvkJU3krZM4lplh9OvAIplfZO8ApoGb290cCn7rsXxs4Dbydvf8j2c/vZvXPYofFi/PZD+iCuZptlP19DXS1+mexesn+G3X8T9TAc9nrNQo4l5a+N4v7ZEzJPgl5lxtd9tFAb5fHChiO+eSRhukJ1djqX6zVC1AB+Aw4kb18Rp7uiK7nEojG9JQ4nP2G2pP9+7jc6p/FgnP3JKY57xywBmjj8r0pwLI8+7cF/pe9/y4g3uqfwU6LJ+cTM27qD+AscAwzRuA2q38GOyw4uxXnXabkdy6zt1n23pQpyoUQQrhlt66zQgghbEiShRBCCLckWQghhHBLkoUQQgi3JFkIIYRwS5KFEEIItyRZCCGEcEuShRBCCLckWQghhHBLkoUQQgi3JFkIIYRwS5KFEEIIt/4f9ZI3nFjvi7MAAAAASUVORK5CYII=\n",
"text/plain": [
"Graphics object consisting of 2 graphics primitives"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"M = Manifold(2,'M',structure='Riemannian') # Riemannian manifold M\n",
"X. = M.chart('x y:(0,+oo)') # Poincare half-plane\n",
"g = M.metric() # metric on M\n",
"g[0,0], g[1,1] = 1/y^2, 1/y^2 # nonzero components of the metric\n",
"p = M((0,1), name='p') # point p\n",
"t = var('t') # symbolic variable for geodesic \n",
"v = M.tangent_space(p)((1,0), name='v') # tang. vect. at p, comp.: (1,0)\n",
"v1 = M.tangent_space(p)((-1,0), name='v1') # tang. vect. at p, comp.: (-1,0)\n",
"c = M.integrated_geodesic(g, (t, 0, 2), v, name='c') # geodesic with init.vect. v\n",
"c1 = M.integrated_geodesic(g, (t, 0, 2), v1, name='c1') # geodesic with init.vect. v1\n",
"sol = c.solve() # find the geodesic facing right\n",
"interp = c.interpolate() # interpolate the result\n",
"p0=c.plot_integrated(thickness=2) # plot the geodesic c\n",
"sol1=c1.solve() # find the geodesic facing left\n",
"interp1 = c1.interpolate() # interpolate the result\n",
"p1=c1.plot_integrated(thickness=2) # plot the geodesic c1\n",
"(p0+p1).show(aspect_ratio=1,figsize=[4,4]) # combine plots"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let us show how the other geodesics through the same point may look.\n",
"\n",
"First we repeat the six commands defining the Poincaré half-plane."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"M = Manifold(2, 'M', structure='Riemannian')\n",
"X. = M.chart('x y:(0,+oo)')\n",
"g = M.metric()\n",
"g[0,0], g[1,1] = 1/y^2, 1/y^2\n",
"p = M((0,1), name='p')\n",
"t = var('t')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we plot 6 geodesics facing right, 6 geodesics facing left,\n",
"one geodesic facing up and one geodesic pointing down."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAAFDCAYAAAAef4vuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAABcpUlEQVR4nO3dd3hU1dYG8HcnIQk10nuRoEjvKlgQsIBi7xUVRb2o94oVK1ZQVPzsDXvvYkNEVEQFkaKiSJPeS0xomZRZ3x+L8ZyTQjKZsqe8v+eZJzOTKWvOzJxZZ5e1jYiAiIiIiKomxXYARERERPGMyRQRERFRCJhMEREREYWAyRQRERFRCJhMEREREYWAyRQRERFRCJhMEREREYWAyRQRERFRCJhMEREREYWAyRQRERFRCJhMEREREYWAyRQRERFRCJhMEREREYWAyRQRERFRCJhMEREREYWAyRQRERFRCNJsB0BEic0YUw/AAwDqAygCcI6IFLr+fz+AdiJyqqUQiYhCYkTEdgxElMCMMY8DuAeaTC0AcIKIfLLnfwbAVgC/iMjR9qIkIqo6dvMRUcQYYzoAWC8iGwD023P1JtdNugCoC2BatGMjIgoXJlNEFEkNAby65/wFAJYC+Nn1/8P3/GUyRURxi2OmiChiRGQ6ABhjWgE4FMCt4h1bcDiAPABzLIRHRBQWbJkiomg4Zc/f90pcfziA6SJSHOV4iIjChskUEUVDH+jYqUWBK4wx7QE0Brv4iCjOMZkiomioD2BlieuO3PP3myjHQkQUVkymiCgafgHQyhiTCgDGmK4A7gawDcCvNgMjIgoVB6ATUTTcB6AFgM+NMUsB7ABQDcAXwmJ3RBTnWLSTiCJqT2HOTBHZ7bruJAAfAjheRD61FRsRUTgwmSKiiDLGfAmgL4CmIrJzT3L1A4A8ERlsNzoiotBxzBQRRVofaKHO3XvGTE2A7nvOtBoVEVGYsGWKiCLKGHMUgKMA1ADQCJpY/Z97sWMionjGZIqIiIgoBOzmIyIiIgoBkykiIiKiEDCZIiIiIgoBkykiIiKiEDCZIiIiIgoBkykiIiKiEDCZIiIiIgoBkykiIiKiEDCZIiIiIgoBkykiIiKiEDCZIiIiIgoBkykiihqj6hhjjO1YiIjCJS2I23JFZCIKSW5uLrKyspCbm2s7FCJKDDFxYMaWKSKKPL8f2LwZ2LJFL2/ZotcRESUAI1LpBie2TBFR1WzeDDRqhDwAWQByAdTZtAlo2NByYEQU59gyRURERBTvmEwRUcRNnDjRdghERBHDbj4iijx28xFRZLCbj4iIiCjeMZkiIiIiCgGTKSIiIqIQMJkiIiIiCgGTKSKKOM7mI6JExtl8RBR5nM1HRJHB2XxERERE8Y7JFBFVWps2bWCMKXUaOXKk7dCIiKxJsx0AEcWP2bNno7i4+N/LCxYswFFHHYXTTz/dYlRERHYxmSKiSmtYYozTuHHjkJ2djf79+1uKiIjIPnbzEVGVFBQU4LXXXsPFF18MY8oeA+rz+ZCXl4ft27dHOToiouhhMkVEVfLRRx/hn3/+wYUXXljubcaOHYusrCy0zc6OXmBERFHG0ghEVCXHHHMM0tPT8cknn5R7G5/PB5/PB7NlC2pnZ7M0AhGFW0yURuCYKSIK2sqVKzF16lR88MEHe71dRkYGMjIyAJ8vSpEREUUfu/mIKGgvvvgiGjVqhOOOO852KERE1jGZIqKg+P1+vPjiixg2bBjS0ti4TUTEZIqIgjJ16lSsWrUKF198se1QiIhiAg8riSgoRx99NIKYuEJElPDYMkVEETdx4kTbIRARRQxLIxBR5G3eDDRqxNIIRBRuMVEagS1TRERERCFgMkVEREQUAiZTRERERCFgMkVEREQUAiZTRBRxnM1HRImMs/mIKPI4m4+IIoOz+YiIiIjiHZMpIiIiohAwmSIiIiIKAZMpIiIiohAwmSIiIiIKAZMpIoo4lkYgokTG0ghEFHksjUBEkcHSCERERETxjskUERERUQiYTBERERGFgMkUERERUQiYTBFRxHE2HxElMs7mI6LI42w+IooMzuYjIiIiindMpogoKGvXrsV5552H+vXro0aNGujevTvmzJljOywiImvSbAdARPEjJycHhxxyCAYMGIAvvvgCjRo1wrJly7DPPvvYDo2IyBomU0RUaffffz9atmyJF1988d/r2rRpYy8gIqIYwG4+Iqq0SZMmoXfv3jj99NPRqFEj9OjRA88991y5t/f5fMjLy8P27dujGCURUXQxmSKiSvv777/x1FNPYb/99sOXX36Jyy+/HFdffTVeeeWVMm8/duxYZGVloW12dpQjJSKKHpZGIKJKS09PR+/evfHjjz/+e93VV1+N2bNn46effip1e5/PB5/PB7NlC2pnZ7M0AhGFG0sjEFF8adq0KTp27Oi5rkOHDli1alWZt8/IyECdOnVQu3btaIRHRGQFkykiqrRDDjkEixYt8ly3ePFitG7d2lJERET2MZkiokq75pprMHPmTNx3331YunQp3njjDTz77LMYOXKk7dCIiKzhmCkiCsqnn36K0aNHY8mSJdh3330xatQoXHrppXu/E5eTIaLIiIkxU0ymiCjiJo4bh+GjRzOZIqJwYzJFREmCLVNEFBkxkUxxzBQRERFRCJhMEREREYWAyRQRERFRCJhMEREREYWAyRQRERFRCJhMEVHETZw40XYIREQRw9IIRBR5LI1ARJHB0ghERERE8Y7JFBEREVEImEwRERERhYDJFBEREVEImEwRUcRxNh8RJTLO5iOiyONsPiKKDM7mIyIiIop3TKaIiIiIQsBkioiIiCgETKaIiIiIQsBkioiIiCgETKaIKOJYGoGIEhlLIxBR5LE0AhFFBksjEFF8GTNmDIwxnlOTJk1sh0VEZFWa7QCIKL506tQJU6dO/fdyamqqxWiIiOxjMkVEQUlLS2NrFBGRC7v5iCgoS5YsQbNmzbDvvvvirLPOwt9//13ubX0+H/Ly8rB9+/YoRkhEFF1Mpoio0g466CC88sor+PLLL/Hcc89hw4YN6NevH7Zu3Vrm7ceOHYusrCy0zc6OcqRERNHD2XxEVGU7d+5EdnY2brjhBowaNarU/30+H3w+H8yWLaidnc3ZfEQUbjExm49jpoioymrWrIkuXbpgyZIlZf4/IyMDGRkZgM8X5ciIiKKH3XxEVGU+nw8LFy5E06ZNbYdCRGQNkykiqrTrrrsO3333HZYvX45Zs2bhtNNOQ15eHoYNG2Y7NCIia9jNR0SVtmbNGpx99tnYsmULGjZsiIMPPhgzZ85E69atbYdGRGQNB6ATUeRxORkiioyYGIDObj4iijgudExEiYwtU0QUeWyZIqLIYMsUERERUbxjMkVEREQUAiZTRERERCFgMkVEREQUAiZTRBRxnM1HRImMs/mIKPI4m4+IIoOz+YiIiIjiHZMpIiIiohAwmSIiIiIKAZMpIiIiohAwmSIiIiIKAZMpIoo4lkYgokTG0ghEFHksjUBEkcHSCERERETxjskUERERUQiYTBERERGFgMkUERERUQiYTBFRxHE2HxElMs7mo9AVFABbtwJbtjinzZv1b2EhUL06kJlZ/t/y/lejBpDCfD8hcDZf4vP7gV27gN27gfz84P9mZAANGpQ+1a8PVKtm+9VR7IqJ2XxptgOgGLZlCzBvHrBihTdRKnnKy4vM82dmAu3aAfvt5z3tvz/QpAlgYuI7RJQ8RIC1a4ElS0qfli0DfL7IPG9WVtmJVuDUti3QvTtQr15knp+oAmyZIt1BrlypidO8ecD8+fp3zRrbkZWvVq2yE6399tPWDiZaUTF27FjcfPPN+O9//4tHHnmk/BuyZSp+iAAbN5adMC1Zoi1Jsap1a02qevRwTi1acH+Q2GLizWXLVLIpKgL++subOM2fD+TkBPc4xuhR4N6OFhs0ANLTvU35lWnuz8/X7oK1a4G//9auwpJ27HBiLykrCzjoIGDgQD317Amkpga/rWivZs+ejWeffRZdu3a1HQqFoqgImD0bmDZNT7NnA9u3B/cY6elAdjbQrFnF3fpldfNnZmqr1t5awLdsAbZt23scK1fq6eOPnevq19ekyp1k7b8/9wkUVmyZSnQ7dwJTpwKTJwNz5gC//67JSkXq1HF2Ph07Ao0aeZOkunWjszMqKgJWrSr7CHn5cqC4uOLHyMoCjjjCSa46deKRaoh27NiBnj174sknn8Q999yD7t27s2UqXvj9wG+/OcnT9OmVS57S0rQ7razW4JYto7c/yMnxJlgbNwJ//OEcXO3YUfHj1KgBdO0K9OoFDBmi+4Xq1SMdPUVGTOzMmUwlovXrgU8/BSZN0kSqouSpaVNvs3iPHsC++8Z+wlFYqOO5Fi/2Jll//AGsW1f+/Ro1AgYMcJKr7OzYf60xZtiwYahXrx4mTJiAI444otxkyufzwefzwWzZgtrZ2UymbBDR78i0acDXXwPffqsTRsrTsqUeQAXGJwYSptatNaGKZX6/jt0KtLwHTps27f1+NWoAxxwDnHACcNxx/FzGl5jYeTOZSgQimkBMmqTN2z//XP5t99vPmzR17w40bhy1UKNCBFi6VH84pk0DvvlGj2DL06qVk1gNGaItb1Sut956C/feey9mz56NzMzMvSZTY8aMwZ133okGADYDTKaiZcMG4IsvnNanvR1cNG7sfP4HDdIDqUSzfn3pBOvvv8u+rTFAv36aWJ1wAnDAAdGNlYLFZIpCUFgIfP+9JlCTJmmXV1maNnV2CocdBtSuHd04Y4HfDyxY4ByZf/dd+d0a1aoBQ4cCF16oiRWnZHusXr0avXv3xpQpU9CtWzcAYMtUrPD5gE8+AV56Sbv1y+sC32cf7fYeNEgTqA4dkrNlNjdXW+kmTdLttnlz2bfbf39nH9q3b+y3ziWfmPjwMpmKJ4WF+sV//33g8891Z1CWrl2dL3+vXqzVVFJRkY4fCxy1z5hRdldow4bAeedpYsVB1gCAjz76CCeffDJSXeNjiouLYYxBSkoKfD6f53//4pipyBDRz/JLLwFvvFH2RJIaNYDDD3dan7p35+DrkoqLtUU/cHD6559l365+fe0GPPVU4NhjmVjFBiZTVEkbNgDPPKOn9etL/z8tTY80TzgBOP54oE2baEcY3/LzgZkzgc8+A157Tbd3ST16aFJ1zjlJ3Q24fft2rFy50nPdRRddhAMOOAA33ngjOnfuXPYdmUyF14YN+ll96SXt4i+pRQs9EDjuOODAA3W2HVXe0qXaWjVpkvYAlNXK16oVcMUVwCWXJPU+IQYwmaK9ENEf+MceA957r3R5gH320SOjE04ABg/WGWsUuqIiYMoU/ZH6+GOt7u7GbsBS9tbN9y8mU6GrqBsvM1NbTC68UCdYsPUpPLZu1fFnkybpdi85RCAjAzj7bODKK7UngKKNyRSVIT8feOstTaLmzvX+LyUFOOkk4PLLtSWKP+SRtW2bvhcvvaS1d0oKdANeeqmOO0lSTKYi7LffgOee0268suos9eunCdQZZ/CgKtJ8Pp0h/fTT2pJd8vezb19Nqk47ja2B0cNkilxWrQKeekp3miWnLdevD4wYoUlUq1Z24kt2f/wBvPwy8OqrpbsBjQFOPx247TagvG6uJDdx3DgMHz2ayVQw5s4F7rrLW4AyoEUL4IILgGHDdIA0Rd+yZcCTTwIvvAD884/3f02aAJddpvvtZs2shJdEmEwlPRGdTfLYY7rD9Pu9/+/ZE7jqKuCss7QJn+yrqBvwtNOA228HunSxEl7MYstU5c2ZA9x5p3bpuWVmAqecoq1QAweyGy9W7NwJvP667scXLPD+Ly1N9wlXXqktiMk4azLyYmKjMpmyZdo04IYbdMfpVq2atnJceSVw8MH88sWybduAiROBBx8sXRTw1FM1qeIsQMVkqmKzZ2sS9dln3uubNQNuvFFbodiNF7tEtJr8448DH35Yekxbv37AAw8AhxxiJ77EFRM/kkymou3333XH+MUX3uubNtVuvBEjtImY4seuXTqG4oEHdGkLt5NP1qSqe3crocUMJlPlmzVLk6iS+4QWLYCbbgKGD2fLdLxZvVr3Cc8+W7pg8MknA2PHAu3b24kt8TCZSipr1uiP6ksveQctduumO8xTTuGAxXi3a5fuPO+/v/S4qhNP1Pe/Z087sdnGZKq0n37SJOrLL73Xt2wJjB4NXHyxzhSj+JWfD7zzDjBuHLBwoXN9aqoeON9xR+KtQBF9TKaSQm6utlhMmADs3u1c37IlcO+9wLnnsqhmotm9WycSjBtXui7Y8cfrD2iPHnZis4XJlGPWLE2sp0zxXt+qFXDzzTomiklUYikqAl58Ud9394FWrVo63GPUKKBmTXvxxTcmUwmtoECLbN51l7eZNysLuOUWHVjOpvvElp/vJFXutdFSUnRM3N13A3Xq2IsvmphM6Ri7m27Sz4RbmzaaRA0bxtbpRLdzJ/Dww3qAvWOHc32TJvpbcdFFrKoevJhIptgkEm4iwLvv6qrrV1/tJFLp6Xr0sWwZcP31TKSSQWamJs3Llumg1BYt9Hq/H3j0Uf2MfPhh6Vo1CWjixIm2Q7BHRGtEdejgTaT23Rd4/nlg8WKtVcZEKvHVrKklVJYuBf7zH2dG5oYN2u3XtavO4kyCfUKiYctUOP36q9YWmTXLe/3ZZ2uXXiKuxk6V5/MB//d/wJgx3i7fE07QadWJXEMsWVumli7VJUemTnWuq10buOcevZ6Fd5PbokU6Pu7DD73X9++vA9gPOMBOXPGFLVMJo6hId459+ngTqSOO0OnOb7zBRIp0HMwNN2gB0CFDnOsnTdJWqocf1s8Sxb+CAj2A6tzZm0idcoouonv11UykSGf0ffCBLrbet69z/Xff6bjKhx8ue11AijlsmQrVn3/qWIdffnGu69gRGD9efzBZJ4rKIqJrLl59tXdAavfuOiOwTx9roUVEMrVMff+9tlC7Z2+1aqVdvccfby8uim0i2kJ1443aohlw6KE6Czw721poMS4mfmTZMlVVxcVarLFnTyeRSknRgaRz5+oixEykqDyBJWj++kvHTgQ+K/PnAwcdpElWXp7VEClI27bp2KfDD3cSqdRU4NprtTWSiRTtjTHacvnrr8D//ufsE2bM0LFUTz5ZepUMihlsmaqKpUt1+vIPPzjXtW+va7cddJC1sCiOzZqlA1B/+825rlkznRE6dKi9uMIl0Vum3n0XGDlSX2dAnz76/iVbGQwKj+nT9Xdm+XLnuiOP1FUXEnl8ZfBiotWCLVPB8PuBJ57QQpuBRMoYnaU3bx4TKaq6gw7SFs7x44EaNfS6deu0NeOGG4DCQrvxhShhZ/P5fNqyeMYZTiJVu7ZOKPjpJyZSVHWHH64HV5df7lw3daqu+/nii5zxF2PYMlVZK1dqReJp05zr2rbVvuzDDrMWFiWgFSu0lePzz53rDj0UeOstoHlza2GFJBFbpv7+W5Mo9/qap5yiZS/i9X2i2DRlii4rtGaNc93QoTq+smlTe3HFBrZMxY0339SjAXcidcUV2rfNRIrCrU0b4NNPtYxCYMbXjBk6OP2rr2xGRgEff6zjJQOJVGam1ox67z0mUhR+Rx+t67oOG+Zc9+mnQKdO+lkk65hM7U1xsc6sOOccYPt2va5FCz1KePJJXQqAKBKM0UHo33/vjI/YsgU45hitU8Xp0nYUFgLXXQecdJIuFQUA7doBM2dqywEnnVCk7LOP9oR8/LGznl9Ojn4W776b3X6WMZkqzz//6HiVBx5wrjv/fGDBAuCoo6yFRUnmoIOc2aGA7jDvvBMYPBjYtMlubMlmzRqtHffQQ851p52mrVPdulkLi5LMCSfo79BppznX3X67djnv3GkvriTHZKosixcDBx8MfPGFXk5N1RoxL7+sa+sRRVP9+rrExNixzqLYU6dqt9/331sNLWlMmaKDyX/8US9Xq6Zjo955J3nWV6TY0aCBfvbGjXNaQ997DzjkEB3fS1HHAeglTZ4MnHWW04Rfr55Oex440G5cVDkiusBwYeHeT2lpQN26eqpRI366Z777TpcnWr9eL6emaqXtG26I7dcQrwPQi4u1JfCee5xulFat9IcsXmbvimiLRU6OnoqLNRks75SertX6Y/nzRI7PPtOhKIG6dA0aAO+/r7MBk0NMfFCZTAWIaOn+G25wCqN17qz9023b2o2NlAiwerV2t6xdq6d165zzgcu7dgX3uNWq6XiEQHIVOLVoofXDAqcGDWLjB2bjRt15uidEXHSRzuyJ0RXnJ44bh+GjR8dXMuXzadf+u+861w0dqi3U9erZiytARLt6Fy1yTuvWOUmT+xTsMkU1a+pA+ubNtd5Z4HzgcsuWej4Wvg+kRWJPOMGpnJ6Wpr0pl11mN67oiIkPIZMpQFsyRowAXn3Vue6kk4BXXtGaMRR9hYW6g5g3zznNn2+3Knjdut7k6oADtHWiWbPox1JcDNx1l3fg6dChwNtvO3WqYkm8tUzl5QEnn+wkrIEWwOuvd7pao2n1auDnn7Vivjt5CrSg21C3rnY19+jhnNq3j9mEPuHl5GivypQpznVXXOGdFZyYmEzFhHXrdKf588/OdbffDtxxh52dZrJatUp3ArNmaeK0YIG2DASjTh09Wm7USLsp0tLK7sZIS9Nkrawj+MCszWBkZ2uJjMCpXbvoHbG/9x5w7rm6sC6gi6V+8omOs4ol8ZRMbdig62rOn6+Xa9TQ7exenDqSRDRRmj5dx8R9/33VxsHUqVO6tbVuXefzX1ioLVYlu8B9Pm39XLcu+O9DZqaWkenRQz+LRx9t52AjWRUVATfd5J0k0b+/fn4bNLAXV2QxmbJu+XIdC7VihV6uUUOb8N2zJCgydu/WH4vJk4Evv/QuClueFi30SLht27K7HsJRqqKoSGdybtumnw93K8CiRd6ieeVp0sRJrI46SluwIumbb7QlNdBq16GDbtOWLcP+VE899RSeeuoprNjznenUqRNuv/12DKko0YiXZGrpUk0AAkt41K+vY1IiOT5KRNfumzpVvxMzZniXpSlPq1beltL27bVGWb162m2dmhp6bNu3l92lvmyZJpuBsXt706WLlvQYPFiLz2ZkhB4X7d0rr+g6kYGDrAMOAL7+OlETWyZTVi1eDAwa5Pw4tmql46O6d7caVkJbtEirek+erD8a+fll384YYP/9vd0H3bvHxg/vzp362Vm0SJd6+P57bdUM7LTK0qmTJuinnabnI9FqNX++tpxs2KCXmzfXhKpTp7A+zSeffILU1FS0a9cOAPDyyy9j/PjxmDdvHjrt7bniIZmaM0e3YSCRadVKW0vbtw//c4noe/bee3pavLj822Zm6uziQw/VxKR9e2C//WKjO3fjRm83/Lx5wJIl5d++Rg0tLzF4sJb7yM6OVqTJZ9Ys7XUJJLzZ2dptnXjr+jGZsubPPzWRCvzwJHbWbtfmzVpB/pVXvMtuuKWk6JH/4MHaUti9e3wVRM3PB2bPdrplfvyx/O6R9u2dxKpbt/AmVsuXa6tKYBDqPvtoleRDDgnfc5ShXr16GD9+PIYPH17qfz6fDz6fD2bLFtTOzo7dZOqrr3QpmB079HLnzpr0h7OauYiuvxhIoP7+u+zbZWVp4hRo3ezVK75ac7Zv19poX3+t2/CXX8ovKNm3L3DBBcCZZ2oXJIXXihW6Tw20tLZure9LYiWxMZFMQUQqe0oM8+aJNGggol9vka5dRTZutB1VYsnPF3n/fZETThBJS3O2tfvUvLnI8OEi77wjsm2b7YjDq7BQZM4ckQcfFOnXr+zXD4hkZ4vcfrvImjXhe+6NG0V693aeIzNT5OOPw/f4LkVFRfLmm29Kenq6/PHHH2Xe5o477hAA0mBPPLl6UCa5gMimTRGJK2hvvilSrZqzzQ47TCQnJ3yPv3y5yE03ibRuXfbnwBiR/v1FHnlEZP58kaKi8D13LNi8WeSNN0SGDRNp0qTsbZCeLnLaaSKTJokUFNiOOLGsXi2y//7Otm7WTOSvv2xHFU7B5DEROyVXMvXzzyL77ON8qHr1Etm61XZUieOXX0SuuEKkbt2yd5i9eomMHy+yYIGI32872uhZs0bk0UdFDj9cfzhLbpfUVP0h+eab8GyXvDyRo45yHj8lReS110J/3D1+++03qVmzpqSmpkpWVpZ89tln5d42Pz9fcnNzJW/ZsthMpp5+2vtenHSSyK5doT9ucbHIl1/qAUVKSun3PCVFZNAgkaeeElm/PvTnixd+vyaMY8eKdOlS9n6iYUOR//5X5LffbEebONavF+nUydnGjRuL/P677ajCxXoiJUmVTM2YIVK7tvNh6ttX5J9/bEcV//x+kc8/FzniiLJ3jM2aidxwQyJ9cUOzfr3Ik0+KDBxY9o9sp076/7y80J7H5xM55xxvwvbhh2F5CT6fT5YsWSKzZ8+Wm266SRo0aFBuy9S/Nm2KvWTq1Ve9ye2IEaG3CuXkaAuTuyUgcEpLEznmGJHnnrP/2mOB3689BddcI9KoUdn7j2OOEfn66+Q6+IqUzZtFund3tm39+tqCHv+sJ1KSNMnUtGkiNWs6H6LDDw/9xyrZFRSIvPJK2UeX1auLnHuuyJQpiddlEU6rV4vcdpseJZbchrVri1x1lcjKlVV//OJikcsvdx4zPV3fkzAbNGiQjBgxYu83irVk6sMPNcEMbJvrrw/tB3vpUpHLLhOpUaP0e9m8ucjdd4ts2BC28BNOQYHIp5+KnHGGSEZG6W3Yq5fI229rFzpV3bZtIn36ONs1K0tk5kzbUYXKeiIlSZFMTZum40YCH56jjhLZudN2VPErL0/k4YdFWrYsvcPbf3+RZ55hohosn0/HlBxySOltmp4ucuWVIuvWVe2xi4tFzj/febwaNUR++CGs4Q8cOFCGDRu29xvFUjL11Ve6XQPb5D//qXoitXKlyKWXehOzwGnAAJH33mMCEKycHJHHHhPZd9/S27RtW5EnnuA+PBS5ud59Ta1a8Z5QWU+kJOGTqV9/FalTx/nQDB0qsnu37ajiU36+DqguazzUwQfrkX5xse0o49+8eSKXXKKte+5tnJkpcu21VUtACgt1LJD7aHTu3CqFN3r0aJk+fbosX75cfvvtN7n55pslJSVFplTU4hUrydQPP3hbj84/v2qf27VrRUaO9CZlgR+m//xHxwVSaAoLRd56S6Rnz9L7nIYNRR5/nIPVq2r7dk32A9uzQQORxYttR1VV1hMpSehkauVKHa8T+LAcd5y2AFBw/H5tXi/rKHHoUJHvv+d4hkjYskVk9Ghv9zSgl2++OfgZkPn5Ikce6f0xWrgw6LAuvvhiad26taSnp0vDhg1l0KBBFSdSIrGRTM2bp4lkYBucdFLwrUYbN4qMGuVt7Qb0oG3MGI7DjAS/X2TqVJGjjy69D2rfXmcAch8UvF27vGNd27aN165o64mUJGwytW2bSMeOzofkwANFduywHVX8+fFHHajv3nkZI3LeeTzyjpbyfryzsnSgczDJwI4d3lINzZvrtP0oeH7sWLvJ1MKFmkC6u/vz8yt/f59PZNy48pNbzgqOjnnzRE4/vXRSNWBAogymjq6cHJHOnZ3t2Lu3tlrFF+uJlCRkMrV7tw4wD3w42rVjHalgLVtW9g5r0CDdmVH0BbqV3PWQAJEePYIb75CT453Rk51d9fFYwbDZMrVihUiLFs5r7tcvuIOrb7/1Hpy5u125b7Fj5syyxxief75O7KDKW73a+/049th4G+dnPZGShEumiot1Noi7K2PpUttRxY/CQpEHHijdCtKhg8hnn7EpPRasWCFy0UXe98cYndZf2daRjRu1e8R9NBqO2kp7YyuZysvz1tfp3r3yBTk3bPAO3ge0nMUVV2hyS3b5/VocuF0773tUq5aOp+IYzsr7/XdvF/jw4fG0v7eeSEnCJVPXXON8GGrU0CKdVDm//+6dMgto7Zenn463o5Tk8OOPIt26ed+vBg1EXnyxcjvB1atFWrVy7nveeZHdedpIpoqLvQPv99+/ci1JRUVa68td4DcwXIBdSbHH59Mu73r1vO/XYYfF86Dq6PvmG++EijFjbEdUWdYTKUmoZOrhh50PQWqq1iyhihUUiNx1l7f7yBiR//1Pp9BS7CosFJkwQY/E3T8iAwdWrqtj/nzvzLaHHopcrDaSqTvucF7bPvtU7od12bLS4wT32UcPKlgzLbZt26a1vkp2xz74IN+7ynrzTe/2e/552xFVhvVEShImmfryS28l4+eesx1RfJg7t3TrxgEHaKsHxY81a7zd24CWsHjvvYrv+847zn1SUkQmT45MjNFOpt57L7jX5fdrEVr3KgmArifHcVHxZdo0nZnmfh8POkikoir9pB580NluaWlhr0sXAdYTKUmIZGrVKi2LH3jzb7vNdkSxz+/XVgj3IsSpqToVn3W44tfkyd6BpIGxDxXNzrn11uBbcIIU1dl8v/7qnXU3fvzeb5+TI3LWWd7t1ratyHffRSY+irwdO3R9P/dBdkaGFhWOn7FAdvj9Ildf7Wy35s1j/YDCeiIlcZ9M+XxaMDLwph93HAcdViQnR+Tkk70/HF276iLFFP+2btVFk93v7377icyeXf59iot1QV73hINwd/FGq2Vq82aRNm2c11LRWLDvvxdp3bp0axSr+CeGGTNKr5N43nkslVORwkLvrPgjj4zlrlLriZTEfTLlzp7btGGtl4rMnVu6+fvGG1nMNNH4/SIvvOBtnUlL03GF5SUWubne6f/HHx/eA5NoJFMFBd6qznubpVhcrOvluRebzsrSituUWHbt0iWZ3Pu9jh1F/vzTdmSxbd0677qhsdvrYz2RkrhOpt5+23mT09P3fuSd7Px+bd52LyBat67IJ5/YjowiacmS0jM0L7yw/GKVS5Z4lwu6887wxRKNZOq665zYGzcufxD+jh2lW+8OOyy0RaUp9r39tneyRs2aIq+/bjuq2Pbtt951Jz//3HZEZbGeSEncJlMLF3q/FE89ZTui2OXzla5L1KdP1Cpfk2UFBSI33eR9//v2FVm/vuzbT5nitNakpor89FN44oh0MvXVV96Dq/ImUaxc6S1aaozOZo3dLgwKp0WLRLp08X4frrqK7//e3H+/s63q1dNad7HFeiIlcZlM7djhLcIX6fo48eyff3SafMkdRzDLaFBiePtt7+LJLVqUP05uzBjndm3bhmf8UCSTqS1bvOtwTphQ9u1mzPAuKVO7Nltnk9HOnaUPMI8/nuOoyuP3e8dU9ukTa78h1hMpictkavhw503t1IlfgPKsXOlNOjMztYYIJa85c7yz/TIzyx4jVFjorbV00UWhP3ekkim/X+TUU51Yjzqq7LFeL7zgraWWnc2p8snu+ee9M5p79y6/xTbZ5eR4F7u/5hrbEblZT6Qk7pKpzz933syaNau06n1SmDtXpGlTZ1s1aMDaUaQ2bChdlPLxx0vfbtkyb1d6ZWpW7UXESiO88IITY/36ZS/zsue5/z0NGsTJKqS+/tq7jEqbNhyYXp65c51xt8boTNjYYD2RkrhKpnJytN5F4EP/zDO2I4pNX3zhncXVrp0OLCYKyM/XgejuBOOee0p3l7/4ovP/unW1OGhVRaJlaskS72f9/fe9//f7dbaq+3VeeaWOIyMKWLBApGVL5zOyzz6sMVae8eO9vy07d9qOSCQGEimJq2TKvfM/6iiOkyrLRx95uzL69dO6O0Ql+f1apNWdaFx3nfd75fd7Z70NGlT1cgnhTqYKCrSqdSC2iy/2/r+oqPTSImPHVv35KLGtXSvSo4fzWaleXWTqVNtRxZ6iIm/L9n//azsikRhIpCRukqlPP3XevNq1OYW5LB984O3/P/XU8mvskB2FhbFXVPaBB7wJxyWXeGc2bd3qbRF+5JGqPU+4k6m77nJiys72Vnn3+UTOPNP5vzGxN+O3qIgzyGLN9u0iQ4Y4n5vMTF2qjLwWLdJtE/hu2W/Fs55IiQiMiKCSKn3DsMrJATp3Btat08vPPw8MH24llJj1/vvAWWcBRUV6+dxzgZdeAtLSrIaVFP75B/jrLz2tXAmsX++ctmwBdu3S0+7dzvuTmgqkp+upRg2gYUOgUSPntO++wH77Ae3aAW3aANWqRfY1PPcccNll+hMC6OfnlVeAlBS9PG0aMGiQnq9VC1i4EGjRIrjn2LwZaNQIeQCyAOQCqLNpk772YC1eDHTpAhQU6LacMQM4+GD9X2EhcPrpwMcf6+W0NH0tZ58d/PMEo6AAWL4cWLIEWLpUz2/apK970ybns1BQoKfiYie+GjWA6tX1b6NGQNOmQJMm+rdNG+CAA/RUp05kXwPpe3PmmcBHH+nljAw9P3iwzahiz8MPA9deq+fbtgV++w2oWdNWNMbWE3sEkXnZMWyYc6QweDC790p6+21vUbULLuARb6Rs2CDy8cfaPTZokEiTJt5WnUic0tJ0VuaFF+pA8VmzIjMt+a23vC2bl13m/a5dfrnzv1NOCf7xw9Uy5ffrtg/EcsMNzv+KikTOPtvbsvDpp8E/R0V279b6W48+KnL++bo4uPs7GKlTs2Y6xOGWW/R1sQs/MgoK9DMe2O7p6ZH5HMWzoiKRQw5xttFVV9mMxnqrlMR8y9SUKcAxx+j5rCxgwYLgj4gT2ccfA6ee6hzlXnSRtjKkptqNK1Hk5gJTpwKTJwNff62tDZVlDFC/vh6tBVoeqlfXXU+gdaKgANi+XVsvCgoq/9iZmcBhhwFHHgkcdRTQrZvTihSKSZOAU05xPk/XXguMH6+vJSdHW0c2bdL/ffIJMHRo5R87XC1Tb7yhLWcA0Lo18Mcfuo1FtHXtuef0fxkZwGefOS1qoSguBubOBb76Sj8PP/wQ3PuVkaGvs1Ytp0UyPV236+7dTsvljh3A1q3BxdaunX4OBg8GBg4EatcO7v5UtsJC/Zy9+65erlYN+Pxz3dakFi8GunfXzy6g34t+/WxEEhMtU7GbTBUUAF27AosW6eWJE4GLL45qCDHthx/0i52fr5cvuQR45pnw/Kgms7VrgXfeAT78EPjxRyexKEvDhkCHDppkdOgAZGdr10zTptpdU9nuOREgL08TjvXrgWXLdEe1ZIn+/fNPp4uwvDhOPhk47TRgwIDQunffeAM47zyny+/OO4Hbb9fzr7+u/wO8iUwlTBw3DsNHjw4tmSovoRMBrrtOux4APZj48EPg+OMr/9glFRRoAv3uu3rQsm1b+betVg3o1AnYf3/tng100TZpop+DWrU0caqMwkJg40b9HKxbp5+FhQudU0VxHHqoJsSnnw40bhzcayavoiL9vL/9tl6uVQuYPh3o0cNuXLFkwgRg1Cg936sXMGuWjYP5mEimYreb78EHnSbEvn3Zvee2YIFO3w1sn3PPjb2BzfFk61aRJ5/UVdKNKbuLJSND12+7/nqdgl9WPaNI2bVLu5Uee0y7ldzTuEue6tcXufRSvX1VvzPPPON9zEBF8ZJdbDfeWPnHDEc3n7ur8eSTnevvvNO53hiRN94I7nED/H4dTHvhhd7vV8nTvvvqbZ54QuTnn6NbDXrVKpF33hEZNUpn67pn77pPKSn6Xj33nK6EQFVTWChy4onOdm3cWGuwkSos9C7P8+yzNqKw3sUn+upjMJlav15n7QV2jnPmRPXpY9rq1d4q1kcfrbOXKDh+v44/uvBCZ2ZKyVP79jr194svYqWeivL7Rf76S8dQnXCCd5kY96lzZ5H/+z+RbduCf46HHvImKB9+KCIiT48aJb49CWcBIFcNGCB//fVXxY8XajL1009OolurlrOI8SuveF/zc88F/1o3bdKDt/bty96OtWrpGJqnnhJZujT4x4+k7dt1SZyRI3Xpn7Lir1lTZMQIkXnzbEcbn3bt0sQ1sD3btRPZuNF2VLHj22+dbdOgQdX2N6GxnkhJzCZT7ppSI0ZE9aljWk6Od4mYXr3Cs25aMikq0mV1evUqP4EaM0aTlXixY4e2Vpx+ukiNGqVfU40aOkA02MWtb73V+xhz58oxxxwj81xH6r/VqyetWraUHRUt6xRKMlVcLNKzpxPLww/r9d9/r4ODA9c/8EBwr2/RIt2/BKo6u0+1a4ucc44mkfFSYsTvF/n9dx2gXl5i1bevllFhS3Zwtm4V6dDB2Y69e3MpMzd3KZKrr472s1tPpCQmk6mffnLelH32Ce/K8vGsqEjk2GOdbZOdrbPLqHLy87UJul270j8wWVm6A5g3L/67k3fs0CVWSi4ZA+iMs7PPFvntt8o9lt/vnR3XvLl2b+7e7dmOxwLyXUW1ZkJJpl5/3Ymha1ftWli2TI+CA9eXnH24N7/8oi1NZXXp9u8v8tpr8ZNAlcfv1y7IK67wLgsUOHXooK16rAZfeStXemuunXVW/O8vwmXVKudALjVVk/rosZ5IScwlU8XFuiJ14MP66KNRedq4cMstznapV49LxFRWUZEui1LWOKOePXWx00Q9wvz9d10+pWRrlTE6zq4yYz927xY5+GDnvr16aZfnu+/+e91vgPw+f36Zd8/Pz5fc3FzJW7asaslUfr6ulxZ4/ilTdAyQu5Vg0KDKJQULF3oXRXa3Qo0apS1ViSgvT7so3WNbAqe2bbWlli1VlfP7797kNNjW0ER2zz3OdhkwIJqJpvVESmIumXr7befN6NRJj0BJF5kNbJeUFC5zUBl+v8hnn5X9A3LkkSLTpiXPUeWWLVox3N2SA2hdqZEjK65XtGGDSOvWzv32THjwu5dzeeGFMu96xx13CABpsOd2QSdTEyY4zxFYRur4453r2reveIzG+vUiw4frd8f9+ps0ERk3LnkGaPv9Ov7v0ENLfyd69dJFf6liH37o3R+zSrravVsnZwS2zRdfROuZrSdSElPJVGGhdwDo5MkRf8q48Pvv3sVcA+NFqHx//+3tEg2cjj1WB50nq507daB1/fre7VKvnsjTT++92GvJI/Jnn5WHTzrJudy8eZmD9ENqmcrJ0dgCrWlz53oXWq2ohbawUJOxOnW8r7dxY50ZGc1ZeLFm+nSRgQNLf0dOPdUZ3E/lu+MOZ5vVrRt7ExNseecdb4IenQNW64mUxFQy9dJLzptw2GHJ02qwNzt2eBPMc8/ldtkbn0/k3ntLz87r00dnnJDKzRW5/fbSY2l69dKEpTyuluOC1FQ5ulEj2eEulbC3hYSrMmbqppucxz7vPJEZM7yVxvd25Pvjjzqb0f36srL085Go3brB8vv1oLVbN+92qlVLZ3OyZ6B8xcU6kzawzXr0SO7kPKC4WKR7d2e7fPBBNJ7VeiIlMZNM+Xze5kH7CyfGhuHDvV/WeB8UG0nz5+vgZPePQrNmukwKE9CyrVunCXrJrr/bby+33IZ/5Mh/b+tr3Vpk5kyn+6xOHZ31VJZgk6m1a52kOD1dB427B//efHPZ99u9W2uBlezSGz6ck1nKU1ysB7MNG5Y+CFm40HZ0sSs313uwO2qU7YhiwyefONukU6doLG9mPZGSmEmmnnrK2fhHHx3Rp4ob7ubSmjVFFi+2HVFsKiwUue8+b/HClBSRa65h2YjK+u670q043bqVOSPnqhEjZI6rdWjXSSfJznPOce43ZkzZzxFsMnXNNc5jXnONt9u2f/+yW01++UXXySvZ2jZzZpU2S9LZtk1nRbpnOWZmaq0yDlAv2/z53vIc0RsnFLv8fhH3eMrXX4/0M1pPpCQmkqldu7QFIbDhk3lMS8CKFdolEdgmL75oO6LYtHatVi13/3h27br3rioqm8+nLVLubrTMTC2C6WrZAyBtANnm2uYjASkOtATVrVt2EhtMMrV5szMDsXp17zipRo20Rc3N79cffHdCnZ6ug8vZVRW8n34qXcD0mGPYsleeRx7xfj5ZskYnSQW2Sbt2kf4eWk+kJCaSqSefdDb6iSdG7GniRnGxHnkHtsmZZ7KbqixffeXtlkhJERk9muMWQjVnTulWqrPP1krbbm+95fy/Xj2RM85wLpc1XTyYZMpdBmTYMO8EjM8+8942J8e73AegBRUXLAjDxkhiO3dq9X/3dm3eXOSHH2xHFnv8fpEhQ5ztNHQo99l+v8gRRzjb5LXXIvls1hMpsZ5MFRVp8cnABueyMd4uz9at9ceCHH6/yP33e7siWrTQwckUHrt3i/znP94f0i5ddJakmzuB6t/feU8aNy41vu/5sWMrl0z984/TKlutmsiBBzrPcckl3tv+9ZfI/vt747z+ehaiDKepU/X9DGzftDRdk5C8Nm70bqfId23Fvm++cbZHt26RTDCtJ1JiPZly108aNCgiTxFXVq1y1iQEtPWFHPn53qWGAJHBgyuuk0RV8+673rIC9evrDjJgyxbvD0jv3s75xx/3PlZlW6buvdd5DHex0DZtvN2Hn39eOrZPP43IZkh669aV7k4fOZJdqCW9/77385js3aJ+v3efELnfM+uJlFhNpkoOUkv2ulIlm4qHD7cdUWzZtk1LZrh36GPGcGBspJVs/alWzXvUPWmS8z93ctOqlffHtjLJ1K5dTtdtSoq3cru7tMWzz3pn63XtGvy6gxScwkKR667zfv+OPpqTPEo67TRn+5x1lu1o7HMX4o7c5DLriZRYTaamT/fuDJO9j/nNN53t0bQpu/fc1q71juPJzNQvKUVHTo4OQHb/kD70kPP/kuv3Bc6//75zm8okUy+84Ny3RQvn/OWX6//9fpE77/TGccoppcdzUeS88IJ3oH+fPmwZdtuwQSdhBLZPsreWFhZ6yx6Vs+xUiKwnUmI1mXIXPHvllbA/fFzZvt37I/TRR7Yjih1Ll3rXZmvUSBdwpegqLBQZMcKbyNxyiyY469d7Z58GTgMGOPevKJny+7WWWsnHaNRIWyX9fpGrrvL+79pr2TJpw3ff6SL0gffhgANYNd3t5ZedbdOuHSfFPPaYsz3OPz8Sz2A9kRJrydSKFc5g1RYtOGB09GjnwzZ0qO1oYseSJd4WijZtWG/LJr9fu1bdCc0NN+j1TzzhXOeuuxOoVVVRMvXDD2Xf/7XXNGG64grv8z74oJ1tQOr3370lbbKzdcwn6ffBPSRhbysDJIMdO5xlodLTI9GSaT2REmvJ1K23Oh+0u+8O60PHnSVLnB+P9PS9rzWWTJYu9SZSnTppdx/Z5z7SBHQGXWGhdvmUbFm67DIRqcRsvrPOKn3fAQM0kbrsMuc6Y7RaN9m3fLm2vDChKm3+fGdcX82aImvW2I7ILvd4u/Hjw/3o1hMpsZJMFRToau2AFggsWYAv2bgXix092nY0sWHNGi0LEdgunTvr1GOKHU8/7U187rlH18NzJz2ADiLPydl7y9S6dTrl3n2/lBSR337TRC3wmCkpnHIea1av9iZU7dtzDFWAu7xIZLq34sfixc62aNcu3N3z1hMpEUEKom3SJGDDBj1/4olA06ZRDyFm/Pgj8NFHer5ZM+Dmm62GExO2bQOOOQZYuVIvd+oEfP010KiR3bjI67LLgGeecS7feiswbx5w+ul6WUT/7toFvPHG3h/rpZeAoiLv/S6+GPjsM2D8eL1sDPDaa8A554TtJVAYtGgBfPMNkJ2tlxctAo49Fti+3W5cseDuu4G6dfX8a68Bv/9uNx6b9tsPGDRIzy9dCkybZjeeCIh+MuXeAV9+edSfPmaIADfd5FweMwaoVctaODHB59ME+48/9PK++wJffcVEKlaNGAE88IBz+aqrgKOOAtLTvbd76aXyH0Ok9P9r1gR69ABGj3aue/pp4OyzQ42YIqFFC2DqVD0gBIDZszWpDiTIyapePecAWYQHy+7f+6efthdHhBgJHAlWrNI3LNeqVUDr1nq+bVtgyRIgJfr5XEz4/HPguOP0/P77awKRlmY3JptEgGHDgFdf1cuNGgE//AC0a2c3rnhXUABs3gzk5morUWGhtvKkp2vyXrcuUL9+aN/Dm24C7r9fz9eoAZx8MvD6697bTJ8OHH448gBkAcgFUGfTJj1K7dfPe9uLL9Yj+YICvXzffd7EKljFxcCWLcA//wA7dzqPm5amiVtWFtCwIVCtWtWfg4AFC4DDDwdycvTylVcCjz1mNybbdu/W/fuaNXr5+++BQw+1G5MthYVAy5bAxo363Vu/HmjQIByPbMLxIKGK7q+3u7n/oouSN5ESAe64w7l8773JnUgB+mMcSKSqV9cuHiZSlbd5s7YIzJ+vP2pLlgArVmgSUZHUVO1ub9tWd/ydOwO9egE9e2pyVJH77tPnevttTdimTdP77drl3Oatt8q+74svei/vs492fQcSnksv9bbg7s2OHcCcOcDcuboNFi8G/v5bhxX4/RXfv1EjoE0b7ZLo3Flbx/r00RYGqljnzvreHXmk/nA+/jjQsSNwxRW2I7OnenXtdbjkEr08Zoy24iWjatWA888HHnxQWy3ffTehPhvRa5kSAbp0cbpwli3TnXcymjwZGDJEz3fvrjt/ExPJtR1ffQUMHuz84L3/PnDKKXZjinU7d+pYssmTdczKX3+F/znS0jSZGDRIx8EceKAmXmXJz9fb/fijXm7ZEli92vl/w4bA5s3elqkVK4CuXYG8POd2zZoB69bp+aOO0qS6vBajoiLgp5+0lffrrzWRqkzSFKxOnYCBA/UzOmCA/kBS+V5+GbjwQj1frRrw3XdA375WQ7KqqAjo0EFbYQH9jiTr9vj1V/3NA7RF+ocfwvGoMfHjGb1kyr0RDzkEmDEjpIeLWyLazBv40XnvPeDUU+3GZNPq1doCsHWrXh4zxttqR46CAk0cXntN/+7eXf5tjdGxLC1aaItL3brapVWtmn4GfT5tydm2TVtuVq/W1q29adJEP6sXXKBJVskDgA0btEUrkAxVq6YtFC6eZOr5550jdgDIyNC4AD3Qmj27dKuQiO6AX30V+OCDilveGjfWxK5xY90GtWtrF6cxuj137tRuqU2bdBusW+cMgi9LrVrA0KHAeefpRIlkb1Euz7XXAg8/rOebNdPJCck89vGFF4Dhw/X8scfqQUIyEtEDqAUL9HJ4GlViIpmKXmkEd52Jp54K+eHi1nffOduhY8fkruBcWOgtbnfcccm9PcqzapXIjTc669aVPKWl6aLAo0aJvPGGyIIFVau6nJsrMnOmyDPP6NqQ7duX/XyBul9PPFF6KZeZM73LjZQ4eUojuFdBcJ+qVxf59Vfv4/7zj8iECd51AsuKacQIkeef1yr5VVlmZvdufe5XXxX573+1dpZ7HUD3qWlTkdtuY3mXshQWivTv72yrIUOS+7vt8+l6lYHtMXeu7YjsGTfO2Q533RWOR7ReFkEkWnWm/H5nSZC0NF1tPlmdeKLzQXr1VdvR2OVeZ61VK102hBx//ily3nlaj63kD3nDhiKXXqoLDUdybbq1a0UmTtTPbUZG6Tjq1tWEYutW5z4TJlQumcrMLPt2L7zgPNbGjVplvXbtspOuU0/V5Ts2bIjcNsjNFfngA5GLLvKuuxY4paeLXHyxFpolx/r1uhxQYDtNmGA7IrvcqwQkc92pVauc7dC1azge0XoiJVFLpubOdTZe5FaOjn1LljhFCZs3T+5ldH75xUkSUlJEZsywHVHsWLlS5IILnM9K4FStmsgZZ4h88YUe+Udbbq4mOoccUjqhqFNHVzPYsUMPno4/vuJkqqzTuefq/XNzdaWEGjVK32bAAD0Q2bEj+tugoEAT2JNOKp3kpqZqgstK/Y4vvnC2T0aGyMKFtiOyZ+dOJxmvVi25WzQPOsj5XIS+6of1REqilky5l495+umQHiquXX21sx2Seb2m/Hztkglsi9tvtx1RbNi9W+SOO0q32NSvr9dHsvUlWPPniwwb5lQuD5xatBB5+22NtYxuyb0mUy1barX0l18WadzY+7/0dE1U/vzT9it3rFkjctNNpRd5rllTv98+n+0IY8P//udsmwMPtHMgECtuusnZFrfcYjsae+6/39kODzwQ6qNZT6QkaslUx476VMZo028y2rFDj94D3RPubpFk414st3v35G6hC/juO++yHIEutHHjItuNF6oVKzTJKdlKM2SIyHPPBZdMvfqqtjqVbI278srYXtvsn3+0y7pkV2SnTiKzZtmOzr6dO71j3ZK5u2/NGucApHHj5N33uZeXOfjgUB/NeiIlUUmmlixxNlq/flV+mLj3/PPOdrj4YtvR2PPXX87Czmlp2sKRzPLzdeC4+0c4LU2vy8mxHV3l/fmnJlAlu/4OPrhyydRhh5Xu0jvllPgah7R5s67H5h6wnpKiLfPJ+qMZ8MMPzjapWTO5F0Q+9VRnW7z7ru1o7OnSxdkOoTWyWE+kRKKxNt/kyc7544+P+NPFrGefdc5fdpm9OGwSAa6+2inIeO21QLdudmOy6e+/tdZKYAo5oJfnzwceekgLWMaLDh20XMMHHzjLiuTlATNnVlxZPCNDK0MHiny2bg188YXWGwus+RYPGjQAnnhCSzr06qXX+f3APfcARxzhVMFORv36OcuJ7NwJXHed3Xhscu//3curJZuhQ53zU6bYiyNcgsi8qubYY53sM1lbIRYscLZB9+46wDYZffyxsx1atdLm/2Q1ZYp3Zlh6ushDD4kUFdmOLHTbtuksxDK68iocgD5ihEhenu1XELrCQh2Q7+7+bNRI5PvvbUdmT06OdxzdN9/YjsiO4mKR7Gz5d+jLypW2I7Jj+nTns3DWWaE8kvVWKYl4y1R+vlZnBvRotWvXiD5dzAoskwLoumPJWO28sBC4/nrn8kMPVW6pkkT07LNaAT+whtl++wGzZgGjRpVfYTxGTJ8+HccffzyaNWsGYww++uij0jeqW1c/86+9poVCKyMrS1u1nnlGC2vGu7Q04NZbtcBoq1Z63aZNWiXevaxWMtlnH2DsWOfydddFpmJ9rEtJcSrEi5RexzJZ9O2r33tAW6aKi+3GE6LIJlMzZjhVmgcPTs4korhYf1QA3cGedZbdeGx54QVdKw3QCvDJWPVdBLjrLm3mD+w4hg7VbqHA6gAxbufOnejWrRsef/zxim987rnAzz9X3FXXqZMuBXPyyeEJMpYcdJAuFzVokF4uKNDtMmGC3bhsuegip2t/zhxdzzEZnXeec/6VV3TfkGzS0nQdR0BXYpgzx248IYpsMhVolQKcjZZsfvgBWLtWzw8ZomuUJZv8fODuu53L48cnX2Itogv2upfKGTVKF4YNHJ3FgSFDhuCee+7BKRWtnejz6RIi8+cDZ5xR/lIiLVtqEjV7ti45VWL5mYRQv76OHR0xwrlu1ChdIDrZpKQADzzgXL7jDl27Ltm0aQMcfrie/+sv4LffrIZjjTsv+PZba2GEQ2QXlnInU0ccEdGnilnvvOOcT9ZWqeeecxLKE08EDj7Ybjw23Hab90fkoYf0BzVR+P3aVeledDgw0WBvVq/WAdoBGRm6oPLAgcBxxwG9eydG4p2WBjz9tA53GDNGr7vlFr3+hhushhZ1Rx+tvwfffgssWaLdXMOG2Y4q+s46C5g+Xc+/805yTsYZMMA5/803cf1diNxCxzt26NiJoiLggAOAhQuDDi7u+f1A8+a6AGxGhi4kmwjjQYJRUKDdPIGZTPPmxU2XVthMmOBNnJ5+OiFmdBpj8NXjj+PIlSv1BzGwwPFeeBY6rsyTtG6t3WLDh4djQdTYMH6890fjuee8Cz4ng++/d1pmDjgA+OMPbbVKJhs3anLt9+u4yUWLEuPAIRgi+hu5fr2Or8zJqXj2b2kxsdEi9+n9+Wen+bZ//4g9TUybPVsTKUBXmE+2RAoA3nzTSaSGDk2+ROq997yJ1GOPJUQihenT8TmAI6+8UpODkonUfvsB556Lonvvxe7774eUN9kgKwt48kkdmHz22aUTppUrtTusXTvgpJO01EK8u/564N57ncuXXw58+aW9eGw47DBvN9ekSXbjsaFxY90OgLbQLVpkNx4bjHHyg50747q7M3LJ1I8/OucPOSRiTxPTPv7YOX/iifbisEXEW0PpppvsxWLD/PnABRc4l++4A7jySmvhhMW8eTrOoX9/DHFfX62aJssTJ2qX7uLFwGuvIe3yy1H92WdhAjWk0tOd2wNAbq7WZrrySp3ltmwZsGqVzuo75hhndqOIfp/69tXn+fPPaL3iyBg9GrjmGj1fXKzjypLtx9S9P3joIXtx2OT+XXD/XiQTd37gzhviTRB1FILjri8V+kKG8albN6eWSCytqxYt337rfAYOOii56mvl5Ijsu6/z+s8/P75ff26uVvcusfjyjkaNRMaP1+rfJRUXe/cD3btLbpMmWmeqcWOR9u2d/515ZtnbZ/16kfvu04XB3fWoUlNFrr3WzmLH4VJcLHLyyc5r6tgxvl9PsPx+Z6kxQGTePNsRRd/Spc7rP/xw29HY8csv3v1A8KzXmBKJ1HIyfr8uzgqINGgQ3z8iVbV2rfMB6d3bdjR2nHmmsw3eeMN2NNHj93uXjOjTRxcxjlfffadFVl3JTH6LFnIeIBPGj5d58+bJyrIKD44f79ynQQORlSslt2lTTaaaNhVZtMhZrxLY+yLoPp+u9deihTepatdOZObMyL32SNu+3bvod7ItNfXUU85rv/RS29HYsd9++vrT0nSdx2RTUKDr1QJ6ABo864mURCyZWr7c+YIMGRLUXRPGyy872yAZVwffvFkXqQW06nF+vu2Iouell5z3vm5dXQw4HhUXi9x7r2etuaLMTLkekPQ9lcwDp2HDhnnv++uvzvtvjMiXX4qIeJMpEZH333e2VY0augDq3uzapYsKZ2Q490tLE3nkkfg9aPvrL12vLvB6PvjAdkTRs327s0B0rVqxvah3pFx1lfPef/ih7Wjs6NfP2QZbtwZ7b+uJlEikKqC7i2/17h2Rp4h506Y55486yl4ctrzxhlMzaNgwnc2YDNatA/77X+fy88/rjLR4k58PnHOOTt8PVKnu3x+pCxfiARH4SuxIXnrpJee+RUVa4Tnw/l93HZ5YsgQdO3bExk2bvM9zyinOmm27dumsvb1Vxa5eHbj9dh2oetBBzvP97386Iy4e61S1bw+4i6Befjmwdau9eKKpVi2deADoDPD337cbjw3uWkvuckLJxJ0nxGvxziAyr8q79VZm2oFukczM5GqVCejd2/kM/PGH7Wii5/TTndd93nm2o6ma3FyR/v2d12GMtgZVdt3ABx907tupk+fzX6plSkTHCbVt69zn2Wcr9zwFBSI33ujcDxAZPDg+13z0+0VOOMF5HcOH244oen76yXndAwfajib6cnKc1t8uXWxHY4e7Nf+BB4K9t/VWKZFIdfO5dwp//x3UXRPCypXO6z/iCNvRRN+SJc7r79nTdjTR8/XXzutu0EBkyxbbEQUvL0+kb1/nddSsKTJpUuXvv26ddtcEkrAS45nKTKZERKZOdZ6zfn1dLLmy3nrL2+03cKB2B8abtWu9Y8hmzbIdUXT4/Tr2DdCkYv162xFFX8+ezncmmM9+opg3L5SDUOuJlEikuvkCtSJq147PLo5Quad3HnqovThsee8953yyVH33+52p7oBWO69f3148VeHzaS2nn37Sy/XqabfD8cdX/jFuv127awBdPiXQFVeRQYOcz8rWrd6q6BU580xdKDVQx23aNL0u3pYpadZM124MGDVKf14SnTH6fgH6PfrwQ7vx2BCoNyXifP+SSYcOuhoAELe1psKfTO3cCaxYoec7dky+qraALqsR0K+fvThsce8Mk2VB47fecnYCvXrF3/IYIpr8BMb61a2r5/v0qfxjLFqkC1oDQJ063vUYK+OBB4DMTD3/+OO61ExlHX64JlS1aunlTz7xJrfx4j//0TFUgK7r+dlnduOJltNOc84nYzLVt69z3v37kSwyMrTQL6D7kcBC8HEk/JmOu/Bchw5hf/i48PPPzvkDD7QXhw0bNzqvv2vXxFkCZG+Ki70tCvffH38HEY8+qqvXAzrI+/PPg18r7J57nMHjN94Y/KLeLVs6g/cLCoBx44K7/8EHa+HDQEHQxx8HXnwxuMewrVo17wLIY8YkR+tUt25OL8a33wLbt1sNJ+rcLbju349kEsgXfD5g+XK7sVRBZJOpAw4I+8PHvOJirXwNAPvuG39dPaGaMsU5f9xx9uKIpo8+cj73hx+uXVbx5JdfdImTgFdeCX4x6hUrdOkgQD/zV19dtViuv95pXZo4ESg5+68iAwdq9fSAkSPjr1r6yScDPXro+TlzgKlT7cYTDcZoZXtAZ2S6Z0Mng9atnd+KuXOTI4EuyZ0vLF5sL44qCn8ytWSJc37//cP+8DFv6VKd4g0APXvajcUGdzI1eLC9OKLJvWTOzTfbi6Mq8vOB8893SgrccIO3y6WyHn3UaZq/+monIdrjiSeeKLs0Qkn16ztrF/p8um5fsC66yCm3sHu3LukTT+OnjNHlZgKSZakV9/7iq6/sxWGDMc7vxaZNzpquySTQzQd484g4Ef5katky53x2dtgfPua5B88F200S70ScOik1agTfuhGP5s93Jhx07gwcfbTVcII2dqwuNAvoWK9gBn4H7NrldKdlZOi4nxJGjhyJP//8E40bNar48a6+2ukmffbZqiVCDz/sdBvMmQM88kjwj2HTKacAbdro+S+/9O5XE1X//s5ajMnWMgXosIiAOB2EHZJ27Zzzcfh5D38y5e7rTIbxMiUtWOCc79zZXhw2rFihi9wCunhlYFHbRPb88875//xHjzDjxfLlOr4L0Jk0L73kjDcKxvvvA//8o+fPOgto0CC0uFq1chaAXb9ex28Fq3p1TfAC78edd8bX0X5qqtO6BmiXZ6KrXVsTegBYuBDYssVuPNHm/r344w97cdiy777OeY6ZArBypf6tX79UU39SCBzlA8k3AN89pde9EniiKix0xglVr64Vw+PJbbdpVxqgM9+qmvy/+qpz/pJLQo+r5OMEBsYH66CDdIYioOUa3JME4sGFFzotNa+9tvfK8InCXUpm5kx7cdjQsaNz3v07kiyaNnUOwAN5RBwJbzJVXKzLaQA6MycZBfp6U1OTr5vTPaU3Gbr4vv4a2LZNz594IpCVZTeeYCxapEv+AHrgc8stVXuczZt1OwB6ZBmuJProo4FAl+Bnnzm1q4J1113OQd3zzwdXbsG2xo2BY47R86tXJ0dy4d5vJNustkBJDMA7kStZpKQALVro+Xj6nu4R3mRq82ZnfENgoyQTER2ADuh4h6p0mcSzefOc88mwJqO7Hs4ZZ9iLoyrGj3dmDF1/fdUTwU8+cVpMTj89fN2caWlOjbL8fGDy5Ko9TqNGzszCwkJgwoTwxBct7s9VMtRfcu835s61F4cNWVlOF3kcjhkKi+bN9e8//zgTueJEeJOpQKsUoE12yWbbNqc+irv/NxmIOIMmW7RI/JIQIsCnn+r5zEynBSEebNsGvP66nq9TB7jiiqo/lns800knlXuzSs/mcwuMmyr5PMH63/+cYqAvvKCFhePF0KHOYPzA5y2RtWnjJPa//mo1FCsC44zXrXO64JNJs2bO+fXr7cVRBeFNpjZudM43aRLWh44Lq1Y555NtGZ21a4HcXD3fpYvdWKLB59NxPT17al2pGjVsR1R5b76prT2AjsupU6dqj1Nc7K2YvpcCtUHN5gvo31/HogFaa6mqtXcaNgTOPlvP5+YCH3xQtcexoX597To98EAdkxdPJR6qwhhn7N6aNUBent14oi3wuyGirz/ZuPMGdz4RB8LfzRcQzE4zUbg//Mk2ZizZKt9nZuoMsTlz4uvHGXDGSgGhDRhfsADIydHzAwY4g6XDJTPTGYO1erWzTFVVDB/unHe//ngwdaqOR7ztNmf9skTm3n8k29gh9+9GYGZ0MnGvmhBswV7LwptMuaeyhjo9Oh65myXdzZXJIDBWDPAWX0sG8VQCYsMGpy5Wp06htSK6Z28GFmoNN/fsrlAWgO3XT0suADpgPp5aPOLp8xUO7v2He7+SDNzDY9zDZpKFO2/YutVeHFUQ3mQqcJQKaLN/snE3SzZubC8OG9ytBslYXyxeuAdy72WMU6XMmeOcj9QalO41y9zPFyxjnDFYhYXJsURLvHKPNw2lNTIexXE3V1jUq+ecD8yUjhBjTF1jzCPGmMeNMZONMRcbYzKNMY/tue51Y0zHih9JhTeZChTuA4B99gnrQ8cFdzdnsIu8xjv3VNZk6+KMJ/n5zg57yJDQHuv3353z7urN4dS9e9nPVxWB5UpatIivQejJJtCCCMTlFPmQuH83kq1oKeBthHHnE2FmjEkH8CSAB0TkSgCXAXgewNsAHgYwCcAZAEov51CO8CZT7qbzeKq5Ey7ulrlEn81WkrtJOhnLYsSLyy/X9+rXX72tPsEScQoLtm4duQK9jRs7O9hQCxkOGKCVpVet0vUIKTYFpscDcTejK2Tu340It8zEpNq1nfOBmfGRcTmAF0Uk8MOVD8AAWCEiywGkAlgC4M3KPmB4RzO6C+slY/XzwGw2IPmSyZNP1ub5rVu9XwiKPcaE3pKUk+N83t1rapXjiSeewBNPPIGPgx1UaoyOofn5Z53gUVBQ9TFE1at7q0xTbGrcGDj2WB0/FErCH4/cPToRbJmJWe68IbLJVI6ITHFdDhQ4mwwAIvIFgC+CecDwJlO7dzvn42mqeLi4W+aSLaG46irbEVA0ubtfKlEGZOTIkRg5ciTymjULvrWhVStNpkR0hlOy1XBLNtWqadX7ZBS9lpnYVLOmc96dT4SZiLxa4qoBAIoBzKjqY4a3m8/94gNF8pJJoGKrMU59HKJE5B4cG+kCvck+KJeSh7sRIs4qgIeFO2+IYDJVhoEA5ohIlTPY8CZTBQXO+WSbzgs4hRAzM8O3rAZRLHKPD3TPwKlIYJZrMLNdk30cCSUP90F44PckmWRkOOfd+UQEGWPqAugG4NsS1wdVhC+8yZS7Om8yFJcrKfDmuz8QRInI3aUdTAX1777z/q2MZO/6oOSRmuociEcpmYgp7rwhQtX+jTENjTE/G2Pu2HPVYGgu9LP7NgD6BfO4lcp4jDEm1z24ujyFhYE7xFdRvHC57TY9Yk9JSc7XT8mjZ0/gsce0xECXLuV+3n0+H3yuNca270mG8oL5fhx4IDBunA5OPeAAfrcosT38sCYVjRol32fdXbKkoKBSrz8rK6sOgO0ilV5vqj+APgA+N8ZUB3AmgHUAagGAMaYmgEcB3BhM6KYyz2+MqQOgEtkUERERUVRliUilMk9jTG0AEwAUQBOosQDqALgPwEoA6dD6U78FE0BlkymTm5vrr/CG27Zp65TIv4NG8/Ly0LJlS6xevRp1qrigap8+fTB79uwq3Tcc9w/1MbgNYmMbhOMxbG+DUGOI9v1LtkytX78eBx54IP788080d9cTimAMsXj/WPg+2P4+xcI2CMdjcBuEcRvUrq2TTIzRcdeVWEklKysrC8G1TEVEpbr5Kh3kXj4MderUqfKHJTU1NaQfn1DvH67H4Dawuw3C8Ri2t0E4YrB9fwCoXbt2XL+GcGwDgPsEgPsEgNsAcG2DIGs0VrZFKtLCOwA9QkaOHGn1/uF6DJvPz20QnsewvQ3CEYPt+4eD7dfAbRC+x7D9/LGwHW0/f7xvg1hQqW6+ParUhJaXl4esrCzk5uaG5UguHnEbcBsA3AYAsGbNmn+b9Fsk8bJD/CxwGwDcBkBYtkFM1CEKJpmq2hMYkwFgNICxIuKr6PaJiNuA2wDgNgAAY0wDAJsBNBSRJFzJVfGzwG0AcBsAibMNIp5MEREFGGMMgNqIgQGjREThwmSKiIiIKARxMQCdiIiIKFYxmSIiIiIKAZMpIiIiohAwmSIiIiIKQVSTKWPM/saYj40xW4wxecaYH4wxA6IZQywwxhxnjJlljNm9Z1t8YDsmG4wxGcaY+cYYMcZ0tx1PNBlj2hhjJhpjlu/5HCwzxtxpjEm3HRtFhjFmtDFmtjFmuzFmkzHmI2NMe9tx2bRnm4gx5hHbsUSbMaa5MeY1Y8xWY8yuPfvCXrbjihZjTJox5h7XPvBvY8ztxpi4bOSp1HIyYfQZgMUABgLYDeB/AD41xmSLyIYox2KFMeZUAM8BuBnANGjBsS5Wg7LnAehq3d1sB2LBAdCDmcsALAXQGfq5qAngOotxUeT0B/AEgNnQfe+9AKYYYzqKyE6rkVlgjOkDYASAoBaUTQTGmLoAfgDwDYAhADYByAbwj8Wwou1GAJcDGAbgDwC9AbwIIBfA/1mMq0qiVhrBVazvcBH5fs91tQHkAThSRL6OSiAWGWPSAKwAcIeITLQcjlXGmCEAHgZwKvSL1ENE5lsNyjJjzPUArhCRtrZjocgzxjSE/oj2F5HptuOJJmNMLQBzAfwHwK0A5ovI/6wGFUXGmHEADhGRw2zHYosx5lMAG0VkuOu69wHsEpHz7UVWNdFsTtsKYCGAC4wxNfckFpcB2AhgThTjsKkngOYA/MaYecaY9caYL4wxnWwHFk3GmMbQVpjzAeyyHE4syQKwzXYQFDWBFV2T8T1/AsBnIjLVdiCWnADgF2PMu3u6fOcZYy61HVSUzQAwyBizPwAYY7oBOBTA51ajqqKodfOJiBhjjgLwMYDtAPzQRGqwiPwTrTgsC7Q4jAEwCtpKdS2A74wx+4tIwu9U91TAfgnA0yLyizGmjd2IYoMxJhvAVdDPAyW4Pd+DhwHMEJEFtuOJJmPMWdADyz62Y7GoLYAroJ+B+wAcCOBRY4xPRF6xGln03A89oPjLGFMMIBXALSLypt2wqibkliljzJg9Awj3duq9Z+fxJLRZ+zDoh+dj6JippqHGYVNltwGc7X2viLwvInMAXARdRPp0ay8gDILYBlcBqANgrOWQIyKI7eC+TzMAkwG8KyLP24mcouxxAF0BnG07kGgyxrSEjoc5T0TybcdjUQqAuSJys4jME5FnoK31V1iOK5rOBHAegHOgyfUwANcZY4ZZjaqKQh4ztWcsVIMKbrYCwCEApgCoKyJ5rvsvATBRRMaFFIhFQWyDvtBB54eJyAzX/WcBmCoit0QsyAgLYhu8BeB4aAIZkAqgGMDrIhKXX6SAym6HwA/JnkTqGwCzAFwoIv4Ih0iWGWMeA3ASdPzocsvhRJUx5iQAH0K/7wGp0P2BH0CGiBSXcdeEYoxZCeArEbnEdd0VAG4Vkeb2IoseY8xqAONE5AnXdbdCE+0D7EVWNSF38+1Z+b3C1d+NMTX2nC35Y+FHnNe7CmIbzAHgA9Ae2l8MY0w1AG0ArIxgiBEXxDa4GjrgNKAZgC+hRymzIhNd9FR2OwA6NRqaSM0BcBETqcS2p3X+MQAnAzgi2RKpPb5G6dnLLwL4C8D9yZBI7fED9HfAbX/E+e9AkGqgdD5QjDjNB6JZGuEnADkAXjbG3AUtjXApgH2hJRMSnojkGWOeBnDnnqx8JYDr9/z7XXuRRY+IrHJfNsbs2HN2mYissRCSFXtapL4FsApaCqGh/tYCyVImJAk9Ae3SOBHAdmNMkz3X54rIbnthRY+IbAfgGSNmjNkJYGuSjR2bAOBHY8zNAN6BDnsZseeULD4BcIsxZhX2zOiGjiV+wWpUVRTNAehbjDGDobVVpgGoBt2AJ4rIr9GKIwZcD6AIwKsAqkNbYwaKSI7VqCjajgbQbs+pZBJpoh8ORUFgPMy3Ja6/CDopg5KEiMw2xpwMHTt6O4DlAP4nIq/bjSyqrgJwN3QsdSNozcFnANxlM6iqilqdKSIiIqJEFJd9k0RERESxgskUERERUQiYTBERERGFgMkUERERUQiYTBERERGFgMkUERERUQiYTBERERGFgMkUERERUQiYTBERERGFgMkUERERUQiYTBERERGF4P8BU1ljF75qjRkAAAAASUVORK5CYII=\n",
"text/plain": [
"Graphics object consisting of 14 graphics primitives"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"Gr=[] # list of plots\n",
" \n",
"for y in range(-1,5): # 6 geodesics pointing facing *right*\n",
" v = M.tangent_space(p)((1,y),name='v') # initial tangent vector\n",
" c = M.integrated_geodesic(g,(t,0,2),v,name='c') # define geodesic\n",
" sol = c.solve() # find the geodesic points\n",
" interp = c.interpolate() # interpolate the result\n",
" gr=c.plot_integrated(thickness=2) # plot the geodesic\n",
" Gr=Gr+[gr] # add the plot to the list of plots\n",
" \n",
"for y in range(-1,5): # 6 geodesics pointing facing *left*\n",
" v = M.tangent_space(p)((-1,y), name='v') \n",
" c = M.integrated_geodesic(g,(t,0,2),v,name='c') # comments as above\n",
" sol = c.solve()\n",
" interp = c.interpolate()\n",
" gr=c.plot_integrated(thickness=2) \n",
" Gr=Gr+[gr]\n",
" # vertical geodesic facing up \n",
"v = M.tangent_space(p)((0,1), name='v') \n",
"c = M.integrated_geodesic(g,(t,0,2),v,name='c') # comments as above\n",
"sol = c.solve()\n",
"interp = c.interpolate()\n",
"gr=c.plot_integrated(thickness=3) \n",
"Gr=Gr+[gr]\n",
" # vertical geodesic pointing down\n",
"v = M.tangent_space(p)((0,-1), name='v') \n",
"c = M.integrated_geodesic(g,(t,0,2),v,name='c') # comments as above\n",
"sol = c.solve()\n",
"interp = c.interpolate()\n",
"gr=c.plot_integrated(thickness=2) \n",
"Gr=Gr+[gr]\n",
"\n",
"sum(Gr).show(aspect_ratio=1) # combine all plots"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As we can see (one can suspect that) geodesics in Poincare half-plane are semicircles or vertical half-lines."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
\n",
"**Example 22.8**\n",
"\n",
"Show (numerically) that one of the geodesics on the sphere $S^2$ is the \"equator\"."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"N = Manifold(2, 'N') # manifold N=S^2\n",
"c_sph.=N.chart() # spherical coordinates\n",
"g=N.metric('g') # metric on S^2\n",
"g[0,0],g[1,1]=1,sin(theta)^2 # components of g\n",
"p = N((pi/2,pi), name='p') # initial point on S^2\n",
"t = var('t') # symbolic variable for geodesic\n",
"v = N.tangent_space(p)((0,1), name='v') # tang.vector at p\n",
" # define geodesic\n",
"c = N.integrated_geodesic(g, (t, 0, 2*pi), v, name='c')\n",
"sol = c.solve() # find the geodesic points\n",
"interp = c.interpolate() # interpolate the result"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The geodesic is computed, but to show a 3-dimensional picture\n",
"we need the embedding $ \\psi: S^2\\to R^3$."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Graphics3d Object"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"M = Manifold(3, 'R^3') # manifold M=R^3\n",
"c_xyz. = M.chart() # Cartesian coordinates\n",
"psi = N.diff_map(M, (sin(theta)*cos(phi), \n",
" sin(theta)*sin(phi),cos(theta)),\n",
" name='psi',latex_name=r'\\psi') # embedding s^2 -> R^3\n",
" # plot the geodesic:\n",
"p1=c.plot_integrated(c_xyz,mapping=psi,thickness=3,color='red',\n",
" plot_points=200, aspect_ratio=1,label_axes=False)\n",
"p2=sphere(color='lightgrey',opacity=0.6) # plot sphere\n",
"(p1+p2).show(frame=False) # combine plots"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## What's next?\n",
"\n",
"Take a look at the notebook [Curvature](https://nbviewer.org/github/sagemanifolds/IntroToManifolds/blob/main/23Manifold_Curvature.ipynb)."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "SageMath 9.6",
"language": "sage",
"name": "sagemath"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.10"
}
},
"nbformat": 4,
"nbformat_minor": 4
}