{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Kerr-Schild coordinates on Kerr spacetime\n",
"\n",
"This Jupyter/SageMath notebook is relative to the lectures\n",
"[Geometry and physics of black holes](https://luth.obspm.fr/~luthier/gourgoulhon/bh16/).\n",
"\n",
"The involved computations are based on tools developed through the [SageManifolds](http://sagemanifolds.obspm.fr) project.\n",
"\n",
"*NB:* a version of SageMath at least equal to 8.8 is required to run this notebook: "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'SageMath version 9.1.beta8, Release Date: 2020-03-18'"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"version()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First we set up the notebook to display mathematical objects using LaTeX formatting:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"%display latex"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To speed up computations, we ask for running them in parallel on 8 threads:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"Parallelism().set(nproc=8)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Spacetime \n",
"\n",
"We declare the spacetime manifold $M$:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"4-dimensional Lorentzian manifold M\n"
]
}
],
"source": [
"M = Manifold(4, 'M', structure='Lorentzian')\n",
"print(M)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3+1 Kerr coordinates $(t,r,\\theta,\\phi)$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We restrict the 3+1 Kerr patch to $r>0$, in order to introduce latter the Kerr-Schild coordinates:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"Chart (M, (t, r, th, ph))"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X. = M.chart(r't r:(0,+oo) th:(0,pi):\\theta ph:(0,2*pi):\\phi')\n",
"X"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"t: (-oo, +oo); r: (0, +oo); th: (0, pi); ph: (0, 2*pi)"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X.coord_range()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Kerr parameters $m$ and $a$:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"m = var('m', domain='real')\n",
"assume(m>0)\n",
"a = var('a', domain='real')\n",
"assume(a>=0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Kerr metric\n",
"\n",
"We define the metric $g$ by its components w.r.t. the 3+1 Kerr coordinates:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"g = (2*m*r/(a^2*cos(th)^2 + r^2) - 1) dt*dt + 2*m*r/(a^2*cos(th)^2 + r^2) dt*dr - 2*a*m*r*sin(th)^2/(a^2*cos(th)^2 + r^2) dt*dph + 2*m*r/(a^2*cos(th)^2 + r^2) dr*dt + (2*m*r/(a^2*cos(th)^2 + r^2) + 1) dr*dr - a*(2*m*r/(a^2*cos(th)^2 + r^2) + 1)*sin(th)^2 dr*dph + (a^2*cos(th)^2 + r^2) dth*dth - 2*a*m*r*sin(th)^2/(a^2*cos(th)^2 + r^2) dph*dt - a*(2*m*r/(a^2*cos(th)^2 + r^2) + 1)*sin(th)^2 dph*dr + (2*a^2*m*r*sin(th)^2/(a^2*cos(th)^2 + r^2) + a^2 + r^2)*sin(th)^2 dph*dph"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g = M.metric()\n",
"rho2 = r^2 + (a*cos(th))^2\n",
"g[0,0] = -(1 - 2*m*r/rho2)\n",
"g[0,1] = 2*m*r/rho2\n",
"g[0,3] = -2*a*m*r*sin(th)^2/rho2\n",
"g[1,1] = 1 + 2*m*r/rho2\n",
"g[1,3] = -a*(1 + 2*m*r/rho2)*sin(th)^2\n",
"g[2,2] = rho2\n",
"g[3,3] = (r^2+a^2+2*m*r*(a*sin(th))^2/rho2)*sin(th)^2\n",
"g.display()"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"g_t,t = 2*m*r/(a^2*cos(th)^2 + r^2) - 1 \n",
"g_t,r = 2*m*r/(a^2*cos(th)^2 + r^2) \n",
"g_t,ph = -2*a*m*r*sin(th)^2/(a^2*cos(th)^2 + r^2) \n",
"g_r,t = 2*m*r/(a^2*cos(th)^2 + r^2) \n",
"g_r,r = 2*m*r/(a^2*cos(th)^2 + r^2) + 1 \n",
"g_r,ph = -a*(2*m*r/(a^2*cos(th)^2 + r^2) + 1)*sin(th)^2 \n",
"g_th,th = a^2*cos(th)^2 + r^2 \n",
"g_ph,t = -2*a*m*r*sin(th)^2/(a^2*cos(th)^2 + r^2) \n",
"g_ph,r = -a*(2*m*r/(a^2*cos(th)^2 + r^2) + 1)*sin(th)^2 \n",
"g_ph,ph = (2*a^2*m*r*sin(th)^2/(a^2*cos(th)^2 + r^2) + a^2 + r^2)*sin(th)^2 "
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g.display_comp()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The inverse metric is pretty simple:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"[-(a^2*cos(th)^2 + 2*m*r + r^2)/(a^2*cos(th)^2 + r^2) 2*m*r/(a^2*cos(th)^2 + r^2) 0 0]\n",
"[ 2*m*r/(a^2*cos(th)^2 + r^2) (a^2 - 2*m*r + r^2)/(a^2*cos(th)^2 + r^2) 0 a/(a^2*cos(th)^2 + r^2)]\n",
"[ 0 0 1/(a^2*cos(th)^2 + r^2) 0]\n",
"[ 0 a/(a^2*cos(th)^2 + r^2) 0 -1/(a^2*sin(th)^4 - (a^2 + r^2)*sin(th)^2)]"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g.inverse()[:]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"as well as the determinant w.r.t. to the 3+1 Kerr coordinates:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"M --> R\n",
"(t, r, th, ph) |--> a^4*cos(th)^6 - (a^4 - 2*a^2*r^2)*cos(th)^4 - r^4 - (2*a^2*r^2 - r^4)*cos(th)^2"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g.determinant().display()"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"True"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g.determinant() == - (rho2*sin(th))^2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Ingoing principal null geodesics"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"k = d/dt - d/dr"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"k = M.vector_field(1, -1, 0, 0, name='k')\n",
"k.display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us check that $k$ is a null vector:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"g(k,k): M --> R\n",
" (t, r, th, ph) |--> 0"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g(k,k).display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Computation of $\\nabla_k k$:"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"0"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nabla = g.connection()\n",
"acc = nabla(k).contract(k)\n",
"acc.display()"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"-dt - dr + a*sin(th)^2 dph"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"k_form = k.down(g)\n",
"k_form.display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Kerr-Schild form of the Kerr metric"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us introduce the metric $f$ such that\n",
"$$ g = f + 2 H \\underline{k} \\otimes \\underline{k}$$\n",
"where $H = m r / \\rho^2$:"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"H: M --> R\n",
" (t, r, th, ph) |--> m*r/(a^2*cos(th)^2 + r^2)"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"H = M.scalar_field({X: m*r/rho2}, name='H')\n",
"H.display()"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"f = -dt*dt + dr*dr - a*sin(th)^2 dr*dph + (a^2*cos(th)^2 + r^2) dth*dth - a*sin(th)^2 dph*dr + (a^2 + r^2)*sin(th)^2 dph*dph"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f = M.lorentzian_metric('f')\n",
"f.set(g - 2*H*(k_form*k_form))\n",
"f.display()"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"[ -1 0 0 0]\n",
"[ 0 1 0 -a*sin(th)^2]\n",
"[ 0 0 a^2*cos(th)^2 + r^2 0]\n",
"[ 0 -a*sin(th)^2 0 (a^2 + r^2)*sin(th)^2]"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f[:]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$f$ is a flat metric:"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"Riem(f) = 0"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f.riemann().display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which proves that $g$ is a Kerr-Schild metric."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us check that $k$ is a null vector for $f$ as well:"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"0"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f(k,k).expr()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Kerr-Schild coordinates $(t, x, y, z)$"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"Chart (M, (t, x, y, z))"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"KS. = M.chart()\n",
"KS"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"t = t\n",
"x = (r*cos(ph) - a*sin(ph))*sin(th)\n",
"y = (a*cos(ph) + r*sin(ph))*sin(th)\n",
"z = r*cos(th)"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X_to_KS = X.transition_map(KS, [t, \n",
" (r*cos(ph) - a*sin(ph))*sin(th),\n",
" (r*sin(ph) + a*cos(ph))*sin(th),\n",
" r*cos(th)])\n",
"X_to_KS.display()"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"sqrt(-1/2*a^2 + 1/2*x^2 + 1/2*y^2 + 1/2*z^2 + 1/2*sqrt(4*a^2*z^2 + (a^2 - x^2 - y^2 - z^2)^2))"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"R = sqrt((x^2 + y^2 + z^2 - a^2 \n",
" + sqrt((x^2 + y^2 + z^2 - a^2)^2 + 4*a^2*z^2))/2)\n",
"R"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"#X_to_KS.set_inverse(t, R, acos(z/R), \n",
"# atan2(R*y - a*x, R*x + a*y))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Check of the identity\n",
"$$\\frac{x^2 + y^2}{r^2 + a^2} + \\frac{z^2}{r^2} = 1$$"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"1"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"((x^2 + y^2)/(R^2 + a^2) + z^2/R^2).simplify_full()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Minkowskian expression of $f$ in terms of Kerr-Schild coordinates:"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"f = -dt*dt + dx*dx + dy*dy + dz*dz"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f.display(KS.frame())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Equivalently, we may check the following identity:"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"True"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dt, dx, dy, dz = KS.coframe()[:]\n",
"f == - dt*dt + dx*dx + dy*dy + dz*dz"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"dx = cos(ph)*sin(th) dr + (r*cos(ph) - a*sin(ph))*cos(th) dth - (a*cos(ph) + r*sin(ph))*sin(th) dph"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dx.display()"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"dx*dx+dy*dy+dz*dz = dr*dr - a*sin(th)^2 dr*dph + (a^2*cos(th)^2 + r^2) dth*dth - a*sin(th)^2 dph*dr + (a^2 + r^2)*sin(th)^2 dph*dph"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(dx*dx + dy*dy + dz*dz).display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Expression of $k$ and $g$ in the Kerr-Schild frame:"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"k = d/dt - cos(ph)*sin(th) d/dx - sin(ph)*sin(th) d/dy - cos(th) d/dz"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"k.display(KS.frame())"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"-dt - cos(ph)*sin(th) dx - sin(ph)*sin(th) dy - cos(th) dz"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"k_form.display(KS.frame())"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"g = -(a^2*cos(th)^2 - 2*m*r + r^2)/(a^2*cos(th)^2 + r^2) dt*dt + 2*m*r*cos(ph)*sin(th)/(a^2*cos(th)^2 + r^2) dt*dx + 2*m*r*sin(ph)*sin(th)/(a^2*cos(th)^2 + r^2) dt*dy + 2*m*r*cos(th)/(a^2*cos(th)^2 + r^2) dt*dz + 2*m*r*cos(ph)*sin(th)/(a^2*cos(th)^2 + r^2) dx*dt + ((2*m*r*cos(ph)^2 - a^2)*sin(th)^2 + a^2 + r^2)/(a^2*cos(th)^2 + r^2) dx*dx - 2*(m*r*cos(ph)*cos(th)^2*sin(ph) - m*r*cos(ph)*sin(ph))/(a^2*cos(th)^2 + r^2) dx*dy + 2*m*r*cos(ph)*cos(th)*sin(th)/(a^2*cos(th)^2 + r^2) dx*dz + 2*m*r*sin(ph)*sin(th)/(a^2*cos(th)^2 + r^2) dy*dt - 2*(m*r*cos(ph)*cos(th)^2*sin(ph) - m*r*cos(ph)*sin(ph))/(a^2*cos(th)^2 + r^2) dy*dx + ((2*m*r*sin(ph)^2 - a^2)*sin(th)^2 + a^2 + r^2)/(a^2*cos(th)^2 + r^2) dy*dy + 2*m*r*cos(th)*sin(ph)*sin(th)/(a^2*cos(th)^2 + r^2) dy*dz + 2*m*r*cos(th)/(a^2*cos(th)^2 + r^2) dz*dt + 2*m*r*cos(ph)*cos(th)*sin(th)/(a^2*cos(th)^2 + r^2) dz*dx + 2*m*r*cos(th)*sin(ph)*sin(th)/(a^2*cos(th)^2 + r^2) dz*dy + ((a^2 + 2*m*r)*cos(th)^2 + r^2)/(a^2*cos(th)^2 + r^2) dz*dz"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g.display(KS.frame())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Expression of the Killing vector $\\partial/\\partial\\phi$ in terms of the Kerr-Schild frame:"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"d/dph = -(a*cos(ph) + r*sin(ph))*sin(th) d/dx + (r*cos(ph) - a*sin(ph))*sin(th) d/dy"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X.frame()[3].display(KS.frame())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plots"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [],
"source": [
"ap = 0.9 # value of a for the plots\n",
"rmax = 3"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [],
"source": [
"rcol = 'green' # color of the curves (th,ph) = const\n",
"thcol = 'red' # color of the curves (r,ph) = const\n",
"phcol = 'goldenrod' # color of the curves (r,th) = const\n",
"coordcol = {r: rcol, th: thcol, ph: phcol}"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [],
"source": [
"opacity = 1\n",
"surf_shift = 0.03"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Numerical values of the event and Cauchy horizons:"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"(1.43588989435407, 0.564110105645933)"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"rHp = 1 + sqrt(1 - ap^2)\n",
"rCp = 1 - sqrt(1 - ap^2)\n",
"rHp, rCp"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [],
"source": [
"X_plot = X.plot(KS, fixed_coords={t: 0, ph: 0}, ambient_coords=(x,y,z), \n",
" parameters={a: ap}, number_values=11,\n",
" color=coordcol, thickness=2, max_range=rmax, \n",
" label_axes=False)"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"(t,\n",
" (r*cos(ph) - a*sin(ph))*sin(th),\n",
" (a*cos(ph) + r*sin(ph))*sin(th),\n",
" r*cos(th))"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X_to_KS(t, r, th, ph)"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"[r*sin(th), 0.900000000000000*sin(th), r*cos(th)]"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xyz_n = [s.subs({a: ap, ph: 0}) for s in X_to_KS(t, r, th, ph)[1:]]\n",
"xyz_n"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"[1.43588989435407*sin(th), 0.900000000000000*sin(th), 1.43588989435407*cos(th)]"
]
},
"execution_count": 42,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xyz_H = [s.subs({r: rHp}) for s in xyz_n]\n",
"xyz_H"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
"[0.564110105645933*sin(th),\n",
" 0.900000000000000*sin(th),\n",
" 0.564110105645933*cos(th)]"
]
},
"execution_count": 43,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xyz_C = [s.subs({r: rCp}) for s in xyz_n]\n",
"xyz_C"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [],
"source": [
"xyz_n[1] += surf_shift # small adjustment to ensure that the surface does not cover\n",
" # the coordinate grid lines"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Graphics3d Object"
]
},
"execution_count": 45,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"graph = parametric_plot3d(xyz_n, (r, 0, rmax), (th, 0, pi), color='ivory',\n",
" opacity=opacity)\n",
"graph += X_plot\n",
"graph += parametric_plot3d(xyz_H, (th, 0, pi), color='black', thickness=6)\n",
"graph += parametric_plot3d(xyz_C, (th, 0, pi), color='blue', thickness=6)\n",
"graph += line([(0.03,0,0), (0.03,ap,0)], color='red', thickness=6)\n",
"graph"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {},
"outputs": [],
"source": [
"xyz_n = [s.subs({a: ap, ph: pi}) for s in X_to_KS(t, r, th, ph)[1:]]\n",
"xyz_H = [s.subs({r: rHp}) for s in xyz_n]\n",
"xyz_C = [s.subs({r: rCp}) for s in xyz_n]\n",
"X_plot_pi = X.plot(KS, fixed_coords={t: 0, ph: pi}, ambient_coords=(x,y,z), \n",
" parameters={a: ap}, number_values=11,\n",
" color=coordcol, thickness=2, max_range=rmax, \n",
" label_axes=False)"
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"