{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Isometric embedding of the extremal Kerr throat\n", "\n", "This Jupyter/SageMath notebook is relative to the lectures\n", "[Geometry and physics of black holes](https://relativite.obspm.fr/blackholes/)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'SageMath version 9.4.beta3, Release Date: 2021-06-21'" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "version()" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "%display latex" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left( r, {\\mu} \\right) \\ {\\mapsto} \\ \\sqrt{-{\\mu}^{2} + 1} \\sqrt{r^{2} - \\frac{2 \\, {\\left({\\mu}^{2} - 1\\right)} r}{{\\mu}^{2} + r^{2}} + 1}\$" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left( r, {\\mu} \\right) \\ {\\mapsto} \\ \\sqrt{-{\\mu}^{2} + 1} \\sqrt{r^{2} - \\frac{2 \\, {\\left({\\mu}^{2} - 1\\right)} r}{{\\mu}^{2} + r^{2}} + 1}$$" ], "text/plain": [ "(r, mu) |--> sqrt(-mu^2 + 1)*sqrt(r^2 - 2*(mu^2 - 1)*r/(mu^2 + r^2) + 1)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r = var('r')\n", "mu = var('mu', latex_name=r'\\mu')\n", "rho(r, mu) = sqrt(1 - mu^2)*sqrt(r^2 + 1 + 2*r*(1 - mu^2)/(r^2 + mu^2))\n", "rho" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{\\sqrt{-{\\mu}^{2} + 1} {\\mu}^{4} r + 2 \\, \\sqrt{-{\\mu}^{2} + 1} {\\mu}^{2} r^{3} + \\sqrt{-{\\mu}^{2} + 1} r^{5} + {\\left({\\mu}^{2} - 1\\right)} \\sqrt{-{\\mu}^{2} + 1} r^{2} - {\\left({\\mu}^{4} - {\\mu}^{2}\\right)} \\sqrt{-{\\mu}^{2} + 1}}{{\\left({\\mu}^{4} + 2 \\, {\\mu}^{2} r^{2} + r^{4}\\right)} \\sqrt{\\frac{r^{4} + {\\left({\\mu}^{2} + 1\\right)} r^{2} + {\\mu}^{2} - 2 \\, {\\left({\\mu}^{2} - 1\\right)} r}{{\\mu}^{2} + r^{2}}}}\$" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{\\sqrt{-{\\mu}^{2} + 1} {\\mu}^{4} r + 2 \\, \\sqrt{-{\\mu}^{2} + 1} {\\mu}^{2} r^{3} + \\sqrt{-{\\mu}^{2} + 1} r^{5} + {\\left({\\mu}^{2} - 1\\right)} \\sqrt{-{\\mu}^{2} + 1} r^{2} - {\\left({\\mu}^{4} - {\\mu}^{2}\\right)} \\sqrt{-{\\mu}^{2} + 1}}{{\\left({\\mu}^{4} + 2 \\, {\\mu}^{2} r^{2} + r^{4}\\right)} \\sqrt{\\frac{r^{4} + {\\left({\\mu}^{2} + 1\\right)} r^{2} + {\\mu}^{2} - 2 \\, {\\left({\\mu}^{2} - 1\\right)} r}{{\\mu}^{2} + r^{2}}}}$$" ], "text/plain": [ "(sqrt(-mu^2 + 1)*mu^4*r + 2*sqrt(-mu^2 + 1)*mu^2*r^3 + sqrt(-mu^2 + 1)*r^5 + (mu^2 - 1)*sqrt(-mu^2 + 1)*r^2 - (mu^4 - mu^2)*sqrt(-mu^2 + 1))/((mu^4 + 2*mu^2*r^2 + r^4)*sqrt((r^4 + (mu^2 + 1)*r^2 + mu^2 - 2*(mu^2 - 1)*r)/(mu^2 + r^2)))" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "diff(rho(r, mu), r).simplify_full()" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\sqrt{\\frac{{\\mu}^{2} r^{12} - 2 \\, {\\left({\\mu}^{2} - 1\\right)} r^{11} + 2 \\, {\\left(2 \\, {\\mu}^{4} + {\\mu}^{2}\\right)} r^{10} + 2 \\, {\\mu}^{10} - 2 \\, {\\left(3 \\, {\\mu}^{4} - {\\mu}^{2} - 2\\right)} r^{9} + {\\left(6 \\, {\\mu}^{6} + 4 \\, {\\mu}^{4} + 9 \\, {\\mu}^{2} - 4\\right)} r^{8} - 3 \\, {\\mu}^{8} - 2 \\, {\\left(5 \\, {\\mu}^{6} - {\\mu}^{4} - 3 \\, {\\mu}^{2} - 1\\right)} r^{7} + {\\left(4 \\, {\\mu}^{8} + 9 \\, {\\mu}^{6} + 9 \\, {\\mu}^{4} - {\\mu}^{2} - 1\\right)} r^{6} + 3 \\, {\\mu}^{6} - 2 \\, {\\left(5 \\, {\\mu}^{8} - 6 \\, {\\mu}^{4} + 2 \\, {\\mu}^{2} - 1\\right)} r^{5} + {\\left({\\mu}^{10} + 10 \\, {\\mu}^{8} + 5 \\, {\\mu}^{6} - 5 \\, {\\mu}^{4} + 5 \\, {\\mu}^{2} - 1\\right)} r^{4} - {\\mu}^{4} - 2 \\, {\\left(2 \\, {\\mu}^{10} + {\\mu}^{6} - 5 \\, {\\mu}^{4} + 2 \\, {\\mu}^{2}\\right)} r^{3} + {\\left(7 \\, {\\mu}^{10} - 9 \\, {\\mu}^{8} + 13 \\, {\\mu}^{6} - 7 \\, {\\mu}^{4} + 2 \\, {\\mu}^{2}\\right)} r^{2} - 2 \\, {\\left(3 \\, {\\mu}^{10} - 6 \\, {\\mu}^{8} + 4 \\, {\\mu}^{6} - {\\mu}^{4}\\right)} r}{r^{12} - 10 \\, {\\mu}^{2} r^{9} + 2 \\, {\\left(2 \\, {\\mu}^{2} + 1\\right)} r^{10} - 2 \\, r^{11} + 6 \\, {\\mu}^{8} r^{2} + 3 \\, {\\left(2 \\, {\\mu}^{4} + 4 \\, {\\mu}^{2} - 1\\right)} r^{8} + {\\mu}^{8} - 2 \\, {\\left(9 \\, {\\mu}^{4} + 2 \\, {\\mu}^{2} - 1\\right)} r^{7} + 4 \\, {\\left({\\mu}^{6} + 6 \\, {\\mu}^{4} - 2 \\, {\\mu}^{2}\\right)} r^{6} - 2 \\, {\\left(7 \\, {\\mu}^{6} + 6 \\, {\\mu}^{4} - 3 \\, {\\mu}^{2}\\right)} r^{5} + {\\left({\\mu}^{8} + 20 \\, {\\mu}^{6} - 6 \\, {\\mu}^{4}\\right)} r^{4} - 2 \\, {\\left(2 \\, {\\mu}^{8} + 6 \\, {\\mu}^{6} - 3 \\, {\\mu}^{4}\\right)} r^{3} - 2 \\, {\\left(2 \\, {\\mu}^{8} - {\\mu}^{6}\\right)} r}}\$" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\sqrt{\\frac{{\\mu}^{2} r^{12} - 2 \\, {\\left({\\mu}^{2} - 1\\right)} r^{11} + 2 \\, {\\left(2 \\, {\\mu}^{4} + {\\mu}^{2}\\right)} r^{10} + 2 \\, {\\mu}^{10} - 2 \\, {\\left(3 \\, {\\mu}^{4} - {\\mu}^{2} - 2\\right)} r^{9} + {\\left(6 \\, {\\mu}^{6} + 4 \\, {\\mu}^{4} + 9 \\, {\\mu}^{2} - 4\\right)} r^{8} - 3 \\, {\\mu}^{8} - 2 \\, {\\left(5 \\, {\\mu}^{6} - {\\mu}^{4} - 3 \\, {\\mu}^{2} - 1\\right)} r^{7} + {\\left(4 \\, {\\mu}^{8} + 9 \\, {\\mu}^{6} + 9 \\, {\\mu}^{4} - {\\mu}^{2} - 1\\right)} r^{6} + 3 \\, {\\mu}^{6} - 2 \\, {\\left(5 \\, {\\mu}^{8} - 6 \\, {\\mu}^{4} + 2 \\, {\\mu}^{2} - 1\\right)} r^{5} + {\\left({\\mu}^{10} + 10 \\, {\\mu}^{8} + 5 \\, {\\mu}^{6} - 5 \\, {\\mu}^{4} + 5 \\, {\\mu}^{2} - 1\\right)} r^{4} - {\\mu}^{4} - 2 \\, {\\left(2 \\, {\\mu}^{10} + {\\mu}^{6} - 5 \\, {\\mu}^{4} + 2 \\, {\\mu}^{2}\\right)} r^{3} + {\\left(7 \\, {\\mu}^{10} - 9 \\, {\\mu}^{8} + 13 \\, {\\mu}^{6} - 7 \\, {\\mu}^{4} + 2 \\, {\\mu}^{2}\\right)} r^{2} - 2 \\, {\\left(3 \\, {\\mu}^{10} - 6 \\, {\\mu}^{8} + 4 \\, {\\mu}^{6} - {\\mu}^{4}\\right)} r}{r^{12} - 10 \\, {\\mu}^{2} r^{9} + 2 \\, {\\left(2 \\, {\\mu}^{2} + 1\\right)} r^{10} - 2 \\, r^{11} + 6 \\, {\\mu}^{8} r^{2} + 3 \\, {\\left(2 \\, {\\mu}^{4} + 4 \\, {\\mu}^{2} - 1\\right)} r^{8} + {\\mu}^{8} - 2 \\, {\\left(9 \\, {\\mu}^{4} + 2 \\, {\\mu}^{2} - 1\\right)} r^{7} + 4 \\, {\\left({\\mu}^{6} + 6 \\, {\\mu}^{4} - 2 \\, {\\mu}^{2}\\right)} r^{6} - 2 \\, {\\left(7 \\, {\\mu}^{6} + 6 \\, {\\mu}^{4} - 3 \\, {\\mu}^{2}\\right)} r^{5} + {\\left({\\mu}^{8} + 20 \\, {\\mu}^{6} - 6 \\, {\\mu}^{4}\\right)} r^{4} - 2 \\, {\\left(2 \\, {\\mu}^{8} + 6 \\, {\\mu}^{6} - 3 \\, {\\mu}^{4}\\right)} r^{3} - 2 \\, {\\left(2 \\, {\\mu}^{8} - {\\mu}^{6}\\right)} r}}$$" ], "text/plain": [ "sqrt((mu^2*r^12 - 2*(mu^2 - 1)*r^11 + 2*(2*mu^4 + mu^2)*r^10 + 2*mu^10 - 2*(3*mu^4 - mu^2 - 2)*r^9 + (6*mu^6 + 4*mu^4 + 9*mu^2 - 4)*r^8 - 3*mu^8 - 2*(5*mu^6 - mu^4 - 3*mu^2 - 1)*r^7 + (4*mu^8 + 9*mu^6 + 9*mu^4 - mu^2 - 1)*r^6 + 3*mu^6 - 2*(5*mu^8 - 6*mu^4 + 2*mu^2 - 1)*r^5 + (mu^10 + 10*mu^8 + 5*mu^6 - 5*mu^4 + 5*mu^2 - 1)*r^4 - mu^4 - 2*(2*mu^10 + mu^6 - 5*mu^4 + 2*mu^2)*r^3 + (7*mu^10 - 9*mu^8 + 13*mu^6 - 7*mu^4 + 2*mu^2)*r^2 - 2*(3*mu^10 - 6*mu^8 + 4*mu^6 - mu^4)*r)/(r^12 - 10*mu^2*r^9 + 2*(2*mu^2 + 1)*r^10 - 2*r^11 + 6*mu^8*r^2 + 3*(2*mu^4 + 4*mu^2 - 1)*r^8 + mu^8 - 2*(9*mu^4 + 2*mu^2 - 1)*r^7 + 4*(mu^6 + 6*mu^4 - 2*mu^2)*r^6 - 2*(7*mu^6 + 6*mu^4 - 3*mu^2)*r^5 + (mu^8 + 20*mu^6 - 6*mu^4)*r^4 - 2*(2*mu^8 + 6*mu^6 - 3*mu^4)*r^3 - 2*(2*mu^8 - mu^6)*r))" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "zp(r, mu) = sqrt( (r^2 + mu^2)/(r - 1)^2 - (diff(rho(r), r))^2 ).simplify_full()\n", "zp(r, mu)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{2 \\, r^{7} + 4 \\, r^{5} - 4 \\, r^{4} + 2 \\, r^{3} - r^{2} + 2 \\, r - 1}{{\\left(r^{3} + r + 2\\right)} {\\left(r - 1\\right)}^{2} r^{3}}\$" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{2 \\, r^{7} + 4 \\, r^{5} - 4 \\, r^{4} + 2 \\, r^{3} - r^{2} + 2 \\, r - 1}{{\\left(r^{3} + r + 2\\right)} {\\left(r - 1\\right)}^{2} r^{3}}$$" ], "text/plain": [ "(2*r^7 + 4*r^5 - 4*r^4 + 2*r^3 - r^2 + 2*r - 1)/((r^3 + r + 2)*(r - 1)^2*r^3)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "zp2 = (2*r^7 + 4*r^5 - 4*r^4 + 2*r^3 - r^2 + 2*r - 1)/(r^3*(r^3 + r + 2)*(r - 1)^2)\n", "zp2" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathrm{True}\$" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\mathrm{True}$$" ], "text/plain": [ "True" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bool(zp(r, 0)^2 == zp2)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\newcommand{\\Bold}[1]{\\mathbf{#1}}{\\mu}^{2} r^{12} - 2 \\, {\\left({\\mu}^{2} - 1\\right)} r^{11} + 2 \\, {\\left(2 \\, {\\mu}^{4} + {\\mu}^{2}\\right)} r^{10} + 2 \\, {\\mu}^{10} - 2 \\, {\\left(3 \\, {\\mu}^{4} - {\\mu}^{2} - 2\\right)} r^{9} + {\\left(6 \\, {\\mu}^{6} + 4 \\, {\\mu}^{4} + 9 \\, {\\mu}^{2} - 4\\right)} r^{8} - 3 \\, {\\mu}^{8} - 2 \\, {\\left(5 \\, {\\mu}^{6} - {\\mu}^{4} - 3 \\, {\\mu}^{2} - 1\\right)} r^{7} + {\\left(4 \\, {\\mu}^{8} + 9 \\, {\\mu}^{6} + 9 \\, {\\mu}^{4} - {\\mu}^{2} - 1\\right)} r^{6} + 3 \\, {\\mu}^{6} - 2 \\, {\\left(5 \\, {\\mu}^{8} - 6 \\, {\\mu}^{4} + 2 \\, {\\mu}^{2} - 1\\right)} r^{5} + {\\left({\\mu}^{10} + 10 \\, {\\mu}^{8} + 5 \\, {\\mu}^{6} - 5 \\, {\\mu}^{4} + 5 \\, {\\mu}^{2} - 1\\right)} r^{4} - {\\mu}^{4} - 2 \\, {\\left(2 \\, {\\mu}^{10} + {\\mu}^{6} - 5 \\, {\\mu}^{4} + 2 \\, {\\mu}^{2}\\right)} r^{3} + {\\left(7 \\, {\\mu}^{10} - 9 \\, {\\mu}^{8} + 13 \\, {\\mu}^{6} - 7 \\, {\\mu}^{4} + 2 \\, {\\mu}^{2}\\right)} r^{2} - 2 \\, {\\left(3 \\, {\\mu}^{10} - 6 \\, {\\mu}^{8} + 4 \\, {\\mu}^{6} - {\\mu}^{4}\\right)} r\$" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}{\\mu}^{2} r^{12} - 2 \\, {\\left({\\mu}^{2} - 1\\right)} r^{11} + 2 \\, {\\left(2 \\, {\\mu}^{4} + {\\mu}^{2}\\right)} r^{10} + 2 \\, {\\mu}^{10} - 2 \\, {\\left(3 \\, {\\mu}^{4} - {\\mu}^{2} - 2\\right)} r^{9} + {\\left(6 \\, {\\mu}^{6} + 4 \\, {\\mu}^{4} + 9 \\, {\\mu}^{2} - 4\\right)} r^{8} - 3 \\, {\\mu}^{8} - 2 \\, {\\left(5 \\, {\\mu}^{6} - {\\mu}^{4} - 3 \\, {\\mu}^{2} - 1\\right)} r^{7} + {\\left(4 \\, {\\mu}^{8} + 9 \\, {\\mu}^{6} + 9 \\, {\\mu}^{4} - {\\mu}^{2} - 1\\right)} r^{6} + 3 \\, {\\mu}^{6} - 2 \\, {\\left(5 \\, {\\mu}^{8} - 6 \\, {\\mu}^{4} + 2 \\, {\\mu}^{2} - 1\\right)} r^{5} + {\\left({\\mu}^{10} + 10 \\, {\\mu}^{8} + 5 \\, {\\mu}^{6} - 5 \\, {\\mu}^{4} + 5 \\, {\\mu}^{2} - 1\\right)} r^{4} - {\\mu}^{4} - 2 \\, {\\left(2 \\, {\\mu}^{10} + {\\mu}^{6} - 5 \\, {\\mu}^{4} + 2 \\, {\\mu}^{2}\\right)} r^{3} + {\\left(7 \\, {\\mu}^{10} - 9 \\, {\\mu}^{8} + 13 \\, {\\mu}^{6} - 7 \\, {\\mu}^{4} + 2 \\, {\\mu}^{2}\\right)} r^{2} - 2 \\, {\\left(3 \\, {\\mu}^{10} - 6 \\, {\\mu}^{8} + 4 \\, {\\mu}^{6} - {\\mu}^{4}\\right)} r$$" ], "text/plain": [ "mu^2*r^12 - 2*(mu^2 - 1)*r^11 + 2*(2*mu^4 + mu^2)*r^10 + 2*mu^10 - 2*(3*mu^4 - mu^2 - 2)*r^9 + (6*mu^6 + 4*mu^4 + 9*mu^2 - 4)*r^8 - 3*mu^8 - 2*(5*mu^6 - mu^4 - 3*mu^2 - 1)*r^7 + (4*mu^8 + 9*mu^6 + 9*mu^4 - mu^2 - 1)*r^6 + 3*mu^6 - 2*(5*mu^8 - 6*mu^4 + 2*mu^2 - 1)*r^5 + (mu^10 + 10*mu^8 + 5*mu^6 - 5*mu^4 + 5*mu^2 - 1)*r^4 - mu^4 - 2*(2*mu^10 + mu^6 - 5*mu^4 + 2*mu^2)*r^3 + (7*mu^10 - 9*mu^8 + 13*mu^6 - 7*mu^4 + 2*mu^2)*r^2 - 2*(3*mu^10 - 6*mu^8 + 4*mu^6 - mu^4)*r" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "s = (zp(r, mu)^2).numerator()\n", "s.collect(r)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\$\\newcommand{\\Bold}[1]{\\mathbf{#1}}{\\left({\\mu}^{2} r^{2} + r^{4} - 2 \\, {\\mu}^{2} r + {\\mu}^{2} + r^{2} + 2 \\, r\\right)} {\\left({\\mu}^{2} + r^{2}\\right)}^{3} {\\left(r - 1\\right)}^{2}\$" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}{\\left({\\mu}^{2} r^{2} + r^{4} - 2 \\, {\\mu}^{2} r + {\\mu}^{2} + r^{2} + 2 \\, r\\right)} {\\left({\\mu}^{2} + r^{2}\\right)}^{3} {\\left(r - 1\\right)}^{2}$$" ], "text/plain": [ "(mu^2*r^2 + r^4 - 2*mu^2*r + mu^2 + r^2 + 2*r)*(mu^2 + r^2)^3*(r - 1)^2" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "s = (zp(r, mu)^2).denominator()\n", "s.factor()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Specific choice of $\\mu := \\cos\\theta$" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "#mu0 = sqrt(3)/2 # theta = pi/6\n", "mu0 = 0 # theta = pi/2" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "zpf = fast_callable(zp(r, mu0), vars=[r])\n", "\n", "def zz(r0, r1, verbose=False):\n", " numint = numerical_integral(zpf, r0, r1, algorithm='qags')\n", " error = numint[1]\n", " if verbose:\n", " print(\"Error = {}\".format(error))\n", " if error > 1e-3: \n", " print(\"Warning: error = {}\".format(error))\n", " return numint[0] " ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Error = 1.6403995233359683e-14\n" ] }, { "data": { "text/html": [ "\$\\newcommand{\\Bold}[1]{\\mathbf{#1}}1.4775405364069212\$" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}1.4775405364069212$$" ], "text/plain": [ "1.4775405364069212" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "zz(2, 3, verbose=True)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Error = 5.314212959789082e-11\n" ] }, { "data": { "text/html": [ "\$\\newcommand{\\Bold}[1]{\\mathbf{#1}}2.6279854005699983\$" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}2.6279854005699983$$" ], "text/plain": [ "2.6279854005699983" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "zz(1.5, 3, verbose=True)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Error = 2.628708897735641e-09\n" ] }, { "data": { "text/html": [ "\$\\newcommand{\\Bold}[1]{\\mathbf{#1}}9.33512707880594\$" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}9.33512707880594$$" ], "text/plain": [ "9.33512707880594" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "zz(1.001, 3, verbose=True)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Error = 8.104079286898108e-07\n" ] }, { "data": { "text/html": [ "\$\\newcommand{\\Bold}[1]{\\mathbf{#1}}13.941287346979056\$" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}13.941287346979056$$" ], "text/plain": [ "13.941287346979056" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "zz(1.00001, 3, verbose=True)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAGSCAYAAACrGxZhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAzo0lEQVR4nO3deZhU1Z3/8fdhVSMgKAIKiguioBIjKhgioglucVySjP50HMYsYwxRsxqXLJqZCdnGbBqdGKMxm9EYtxgVEwVccAFRQBAXcGcRRVmUtc/vj1NtN00vVd1Vfau63q/nuU/BXep+6+ZGP55z7rkhxogkSZLaV6esC5AkSapGhjBJkqQMGMIkSZIyYAiTJEnKgCFMkiQpA4YwSZKkDBjCJEmSMmAIkyRJyoAhTJIkKQOGMEmSpAwYwiRJkjJQNiEshLBjCGFAgcd8PYTw8Tae91NtOV6SJKk1yiKEhRAGA1cAqws89DLg7BDC0a0876HAhNYcK0mS1BZdivVFIYSfAF/K/fU94DVgJRBz63YEBgEvxRgH1zuuB/B74NQY46pCzhlj3BRCOAN4IIQwK8a4tMCyTwf+UOAxkiRJbRZijC3vlc8XhfASsAb4GnBvjHFDvW07AQ8DvYHDYoxP1dt2NfBYjPHqNpz7W8CQGOO/F3BMF+A5YHiM8d3WnluSJKk1itISFkL4ILARGBNjfKvBtl7A3cAA4OgGAewAYBxwdhtLuBxYGELYM8b4fJ7HHA08aACTJElZKNaYsH8BvtxIAOsO3AYMB86IMd7f4LjzgWtjjBvbcvIY4wrgPuAzBRxmV6QkScpMsULYPsDf6q8IIXQC/giMJQW0Gxts7wmcCPy9SDVMA07KZ8cQwrbAKOAfRTq3JElSQYoVwk6PMdY0WHc5cDLw/Rjjzxs55qPAJmB2ww0hhBEhhN+FEB4MIXw8hNA7hPDTEMKVIYRbct2fDU0DhoYQ+uVR78nAHQ1b4EIIY0II14cQpocQjgshdAohTAwh/CKEcFWunlG5fU/L1XNlbv0ReZxXkiQJKFIIaxjAQgjfJo3zui7GeGEThx0KPB5j3NTItvOAM0ljyX4DXA38CLgROI7Gux1fzH2OyKPkLboicy13/0masuIh4FrgZ8C8GOM5McbPA08Df8z9vrUxxrNjjGeTukJvyn2HJElSi4oeGkIInwMuJXUzfq6ZXfchTWPR8Pg9gNdzrVQ7AX2A78UYXwO2J80ldkcj3/c2qWVtcAv19QN2jTE+2mDTQcCsmB4X3QnoC9zZYBzbSmA3YGmM8a/11i/N1dm3uXNLkiTVKmoICyGcAFwJPAb8awsD7ncmBaeG+gE35/48hvQE4xMAMca/xBj7xBgnNzwoF57eAXq1UOapwJ8bWd8duCX3548A98QY726wz/7AIuBXDdbvQ5ob7c0Wzi1JkgQUMYSFEMYANwAvAMfFGNc02H5tg0O2JYWmzcQYH44xzgoh7ADsC0wpoIx1pDDVnNNJDww0PO+0GOOLIYQ9gYHAvfW35+YVOxSYErecXG088EBbn/KUJEnVoyghLIQwHLgdWAEcFWNc3mD7GOrGbNXaBHRr5mvHAYHCQlgfmnn1UQhhL9IEtQua+Y7aAfYNz3sQKThutj6EsB8wBPhLAXVKkqQq1+YQFkIYBNwDdAaOiTG+2MhuX6euq6/W26TQ1JRxpJatR/Ks4wOkVrAtxpnVk8/cYONytc1qsP7w3OeUButPI9X5l1wdhcxVJkmSqlSbZswPIfQhBbAdSAHsqUb2+RQwLMbYcCqKF2k5hD0SY1ybZzm1U1M018p1CnVhqimHA9MamXJjHLAwxvhyg/UnA7fHGFfkWsUG5leuJEmqZq1uCQshbE2aoHUojcyGH0IYEkL4HqnlqWErGKT5wXZr4rv7A3tTWFfkaNLTi3Ob+M5RwMsxxiVNfUEIYR+gf8PzhhC6khsP1shhOwBTc9NTfIM0P5okSVKz2tIS9t+k4LMemBRCmJRb343UwvWB3N8jcFMjx08DvhtC2C7G+HaDbf1I0z7cuMVRTTsM+HsjLVi18umK3B54Hfhrg/XbkV5O/vtGjjkX+CypBe2qGKNPSEqSpBaFLR/0a6cTp5aj14HPxhj/1tL+eXzX88AXY4xbvAYp92Tj88B+McZVbTmXJElSMWQ2w3uuxerXpHm72uqTpJaqu5rYPh541AAmSZLKRdav2fkZcHQIodGxYQW4ADivkfm7auXTFSlJktRuMuuOfL+AED4PHBxj/HQrj78Q2D7G+LUmtn8AmAfsGWPc0PpKJUmSiifrljBijFcBnUMIEwo9NoRwFOnpzG80s9tJpAH7BjBJklQ2Mg9hOZ8BRuWmiCjEYuDMGOOmZvYZC1zf6sokSZJKIPPuSEmSpGpULi1hkiRJVcUQJkmSlAFDmCRJUgYMYZIkSRkwhEmSJGXAECZJkpQBQ5gkSVIGDGGSJEkZMIRJkiRlwBAmSZKUAUOYJElSBio+hIWkZwghZF2LJElSvrpkXUAR9ADeAXoBKws81reXS5KUvapsSKn4ljBJkqRKZAiTJEnKgCFMkiQpA4YwSZKkDBjCJElSUa1YAQ8/DNdck/6sxnWEpyMlSVI7ixGWLoX589Myb17d55IlaZ8QYOhQGDMm21rLVYixMmdpCCFMBCaSWvOGAr1ijE5RIUlSEdXUwCuvbBm05s+va+Xq2hWGDIFhw2CffdIybBjstRdsvXVep6nKKSoqNoTVCiH0JDdPmCFMkqTW2bgRFi7cMmjNnw/vvpv22XrrupBVG7T22Qf22CMFsTYwhFUiQ5gkSflbuxaefXbLlq3nnoP169M+vXrVBaz6n7vsAp1KM5rcEFaJDGGSJG1p7Vp45hl4+um6Zd681NpVU5P26ddvy6C1zz7Qv38az9WOqjKEOTBfkqQKtn49LFiwedh6+ml4/vm6sLXzzjB8OBx//Objtvr0ybb2amcIkySpAmzYkLoMG4at555L47kgtWDtuy8cc0wKXcOHp9C13XaZlq4m2B0pSVIZ2bgRXnhhy7C1YEEKYgB9+6awVRu0apcKbtmqyu5IQ5gkSRmoqYFFi2Du3M3D1jPPwLp1aZ8+fRoPW337Zlt7CRjCKpEhTJJUzmKE116DOXPSUhu65s+H995L+2y33ZZBa/jwNHC+nQfIZ6U6fmUDmYawEMKFwMnA3sB7wMPAN2KMCwr4DkOYJKksrFyZQlZt4Kpdaic13XbbFK4atm7ttFPVhK2mVOWvzzqE3Q3cADxOekjgf4D9gGExxjV5fochTJLUrjZsSGO0Goatl15K2zt3Tq/r2W+/zZdddy3ZPFuVzhCWtRBCX2AZMDbGOC3PYwxhkqSSiDG9sqdh2HrmmbpB8jvvnALW/vvXha2994bu3bOtvcJUZQgrtykqeuU+32pqhxBCd6D+rd2jpBVJkqrC229vGbbmzoV33knbe/RIAevQQ+Gss+oCV+/emZatClY2LWEhhADcBvSOMX6kmf0uAb7TyCZbwiRJLaqpSVNAPPXU5svLL6ftXbqklqyGXYm77FL147ZKqSqvbDmFsCuA44AxMcZXm9mvsZawVzGESZIaWL0aZs/ePGzNmQNrcqOO+/WDESPqlv33T2O5unXLtu4qZAjLrIgQfgGcCBwWY1xU4LGOCZOkKhdjGhTfsHXrhRfS9trWrfqBa8SIFMJUFqoyhGU6JizXBfkL4CTg8EIDmCSp+rz3XhqrVT9szZ5dN3arT58UsI4/vi5sDRvmQHmVn6wH5l8BnAacAKwKIfTPrX8nxvhedmVJkrIWIyxeDE8+uXngevbZNK6rUyfYa6/UhXj00XWBa+edHbulypD1PGFNnfzMGON1eX6H3ZGSVOFqamDhQpg1Ky1PPJE+ly1L23v23LIrcfhw2GabbOtW0VRlbC6LMWFtYQiTpMqyYUN6ZU/9wPXkk7BqVdq+885wwAGbL7vuautWB1eV/+tm3R0pSerA3n03jdeqH7jmzq17QfWQISlkHXdcXeDqgC+nlhplS5gkqShWrNiyO3HBgtTV2KVL6j780IfqwtaIEWkCVIkqbQkzhEmSCvb663VBq/az9r2J22yTAlb9wDV8uE8nqlmGsEpkCJOk0lqyBGbOhBkz6j4XL07b+vTZPGwdcEDqYuzcOduaVXEMYZXIECZJxbNs2eZha+ZMeO21tG377WHkSDjwwLrPQYMcMK+iqMq7yBAmSVXqjTe2bOF6NffSuD59Ng9bI0f67kSVVFXeWRUbwkIIE4GJQCdgKIYwSWrSm29uGbhqX1i93XZbBq7Bgw1caldVebdVbAirZUuYJG3u7bdTyKpdZs6EF19M23r1SkGrfujafXcDlzJXlXeg84RJUgVbty7Nw/XYY2l59NE0LQSk6R8OPBA+8Ym6wLXHHul1P5KyZ0uYJFWImhp4/vkUtGpD15NPwvr10LUrfPCDcPDBdcteexm4VDGqsiXMECZJZWrJkrqw9dhj8PjjqasRUsCqDVuHHJLm5XIeLlUwQ1glMoRJ6ghWr05jt+qHrtqB8zvumILWIYek0DVyJPTunW29UpFVZQhzTJgktbNNm9L7Ex99tK5rcd681N24zTYpZJ1ySl1Ll3NxSR2TLWGSVGLLl8Mjj8D06enzscdSy1fnzrDffpuP49pnn/SeRanKVOV/ZhjCJKmINm5MrVzTp9eFrueeS9t23BFGj07LqFGpxesDH8i2XqlMGMIqkSFMUpbeeGPLVq41a1Jr1ogRm4eu3XazW1FqQlX+P8MQJkl52rgR5szZvJXr+efTtn796gLX6NFpTq5ttsm2XqmCVGUIc+SBJDVh1aoUtB56CB58MP25tpXrgx+EY46pC1277morl6TC2BImSTmvvZbCVm3oeuqp9MRinz7w4Q+n5dBD01iurbfOulqpQ6nK/4Sp2BDmC7wltUVNDTz9dF3geuihuvcr7rEHjBmTQteYMTB0qDPPSyVmCKtEtoRJysd776VB87Wha/r0NPt8587woQ/VBa4Pfxj698+6WqnqGMIqkSFMUmOWL09hq3Z54gnYsAF69kxjuGpD18EHO02EVAaqMoQ5MF9Sh7BkCUydCtOmpc+nn07rBw1KYevf/z0Fr333Ta1fkpQ1W8IkVaSXX64LXNOmwbPPpvVDhsDYsXDYYWnZddds65SUl6psCTOESSp7McILL2weumoH0Q8fvnnoGjAg01IltU5VhjC7IyWVnRhh/vzNQ9frr6d5uD74QTjxxBS4PvIR2GGHrKuVpNYxhEnKXIzp/Yr33ZeWKVPS64A6d05zcp1+emrt+vCHYbvtsq5WkorD7khJmXjppbrQdf/9aaLUzp3T04rjxsHhh6enGLfdNutKJbUDuyMlqVSWLElhqzZ4LVyYuhcPOAD+3/+DI45ITzH26JF1pZLUPgxhkkrizTfTeK7a0DV/flo/fDgcd1wKXYcdll4JJElZCSF8KsZ4UybntjtSUjGsWZNC1z//mVq8nnwyjfXac88UuI44InUx9uuXdaWSylAm3ZEhhEOBi2KMH8/i/LaESWqVTZtg5ky49960PPxwmpF+4EA48kj40pfS2K5Bg7KuVJKadDrwh6xObkuYpLwtWlQXuv75T1ixIo3hGjcOxo+Hj30sTZYaqnKIraQ2aPd/aoQQugDPAcNjjO+29/mhglvCQggTgYlAp6xrkTqqt99OXYuTJ6fg9cILdU8wnntuCl0HHwxdu2ZdqSQV7GjgwawCGNgSVtk/XiqyDRvgkUfqWrseewxqatK4rtqWrnHjoFevrCuV1MFk0RL2J+C3Mca72/vc79dgCJOq26uvwl13peUf/4BVq9ITi0cemULXxz4GgwdnXaWkDq5dQ1gIYVtgDjAkxrixPc9dX8V2R0pqnfXr4cEHU+i6+26YOxc6dYJRo+D88+Goo+BDH0rdjpJUrkIII4CvAbsB3wceAr4DdAf6A5fGGJ9s4vCTgTsaBrAQwhjgP4EhwH8DdwFnA3sDXYF9ga/FGB8JIZwGfCR36H7At2OM9xXyGwxhUhV46aUUuO66Kw2oX70a+veHo4+Gb30LPvpR5+uSVHHOA84ELgB+A0zLrdsLuAd4FTiniWNPB75df0UIoRMpgE0AfgRcC/wZ+GuM8YrcPv8H/DGEcB0wN8Z4dm79d4GbQgh9Y4w1+f4AQ5jUAa1bBw88UNfNOH9+atk69FC46CI45hjYf//UAiZJlSaEsAfweoxxYwhhJ6AP8L0Y42shhNHAauCOJo7tB+waY3y0waaDgFkxxpj7zr7AnTHG++vts5LU8rY0xvjXeuuX5mrom/tzfr/DMWFSx7BsGdx5J9xxR3qacc0a2GmnFLiOOSaN8fLl15LKVEFjwnKTrL4XY5wVQpgNvBVjPDzPY88D+sQYv9Ng/WHAyzHGF0MIrwBPxxiPbrDPPaSuyj1ivQAVQrgc+DTQs5AxZraESRUqRpgzJ4WuO+5ITzJCGtt18cXp1UD77eecXZI6nhjjwwAhhB1I47S+W8DhpwNnNPKd03LfuScwEPhp/e25ecUOBW6KW7ZgjQceKHSQf1mEsBDCF4CvAwOAp4EvxRgfyLYqqfysXQtTpqTQ9be/wcsvw7bbpsH0Z58Nxx4LfftmXaUktZtxpFa0KfnsHELYi9QLuKCZ3Y7IfTb8zoOAbRuuDyHsR2od+1E+NdSXeQgLIZxCSptfID3ZcBZwVwhhWIzx5Sxrk8rB0qV13Yz33pu6GQcPhhNOgOOPTy/B7t496yolKRPjgHXAI3nun89risYBbwOzGqw/PPc5pcH603I1/AUghPCZGOM1+RST+ZiwEMKjwBO1Txjk1s0Hbo0xXpjH8Y4JU4cSI8ybB7feCrffnroZQ4DRo1Po+vjHYfhwuxkldSit+idaLi8sLWA82DPA4THGJc3ssxh4LMZ4QoP1k0ljwfZosH4B8FSM8V9zrWInxxgvzaeeTFvCQgjdgANJ83vUN5nU7ypVhZqaFLZuuSUtzz1X1834hS/YzShJDYUQ+pPm7/pznvuPIg28by6A7UOaY2xKg/VdSbmksXPtAEzNTXHxDdI0GXnJujtyB6AzWz7OuZR0EbYQQuhOmoitVo/WnDiEEN55553WHCoVxYYNaRqJv/0tdTcuWQLbb58C1//8D4wdC1ttVbf/ykLbeSWpQvTq1asnsKqRAe/N6UfKCzfmuX8+XZHbA68Df22wfjtgDfD7Ro45F/gsqbvyqhjjm3nWk213ZG4ejteAQ2OM0+utvxg4I8a4dyPHXEKaEbehgroj63VjSpKk7LVmWFFeck82Pg/sF2NcVYpztEbWLWHLgU1s2eq1I01PdjYJuKze33uQZsUt1KqWWsJWrlzJoEGDeOWVV+jZs2crTgEHHXQQjz/+eLsfm+W523rdKvV3N3X8O++kCVNvvz3NVr92LeyzT934rv33h1WrvNdaw3ut8OP951r13WttPb49rluvXr16AaUMR+OBR8spgEHGISzGuD6EMBP4GHBLvU0fA25r4ph1pKcQAAitHJ1cSJNnz549W/0Pq86dO2dybNbnhtZft0r+3bXHr1yZnma88cb0uqD169P8XZdeCiedBEOGbH5c7W3svdY61XyvtZb3WutU4r3W1uPb47qVqgWsnny6Ittd1i1hkFq1fhdCmAFMJ723aRfgqkyrKpKJEydmcmzW526LSv3dq1bBwQf/lBNPTMFr3br0ROMPfgCf/CQMHNjqr86L91r7n7tSr3lbVervrtRrlvXxWV63YgghfAAYA/xHxqVsIfMpKuD9yVrPJ03WOhf4cu3MtXkcW7IpKlauXEmvXr1455132vxfAdWkmq7bqlVpYP2NN6Yux3XrUovXv/5rCl6DBuX3PdV0zYrJ61Y4r1nreN1ap4DrVrJJd0II/wZ8uP5UWOWiHFrCiDH+Evhl1nU01L17d77zne/Q3ZkwC9LRr9vq1elpxhtvhL//PY3xOuSQ9ETjJz8Ju+5a+Hd29GtWKl63wnnNWsfr1jplct3GAr/JsoCmlEVLWFs4Wavaw4YNcM898Ic/wG23wXvvwUEH1bV4DR6cdYWSVNGqcvppQ5jUhJoaePjhFLxuugnefBP23RdOPx1OOQV22y3rCiWpw6jKEFYW3ZFSOZk7NwWvP/0JXnopjev67GfhtNPSdBKSJBWDIUwCXnstBa8//AFmz4bevVNX42mnwZgx0KlT1hVKkjoaQ5iq1tq1aXzXddfB5MnQrRv8y7/Af/0XHH10+rskSaVS1f99P23aNI4//nh22mknQgjceuutLR4zdepUDjzwQLbaait23313rrqqQ0xnlrdCr9mUKVMIIWyxPPPMM+1TcAMxphdlf+ELMGAAnHpqeifjVVeldzf++c8piBUzgE2aNImDDjqIHj16sOOOO3LiiSeyYMGCFo+r9nutNdet3O639nbllVey//77vz8x5ujRo7nrrruaPaba7zMo/LpV+33WmEmTJhFC4Etf+lKz+3m/ba6qQ9iaNWsYMWIEl19+eV77L1q0iGOPPZaPfOQjzJo1i4suuohzzz2Xm2++ucSVlo9Cr1mtBQsWsHjx4veXIQ2njC+xxYvhRz9KA+sPOSS9Qujss+GZZ+Chh+Bzn4NevUpz7qlTpzJx4kQeeeQR7r33XjZu3Mj48eNZs2ZNk8d4r7XuutXK+n7LysCBA/n+97/PjBkzmDFjBkcccQQnnHACTz/9dKP7e58lhV63WtV6nzX0+OOP86tf/Yr9Wxg06/3WiBhjRS7ARGAe8AzpKceerfie9wHxlltuic05//zz4957773ZurPOOiuOGjWq2eM6qnyu2f333x+BuGLFinapqb61a2O86aYYjz02xk6dYuzePcZTT43x7rtj3Lix3ct537JlyyIQp06d2uQ+3mtbyue6ZXm/lavevXvHX//6141u8z5rWnPXzfuszqpVq+KQIUPivffeG8eOHRvPO++8Jvdt4X7LPFdksVRsS1iM8YoY4zDg4PY65/Tp0xk/fvxm64466ihmzJjBhg0b2quMinTAAQcwYMAAjjzySO6///6Snmv+fPjyl2GnneBTn0pTS1xxRepu/NOf4KijoHPnkpbQrNoXx/fp06fJfbzXtpTPdavVnvdbudq0aRM33HADa9asYfTo0Y3u4322pXyuWy3vs/RKo+OOO46PfvSjLe7b3P0WQuhaqhrLmQPzC7BkyRL69eu32bp+/fqxceNGli9fzoABAzKqrHwNGDCAX/3qVxx44IGsW7eO3/3udxx55JFMmTKFww47rGjnWbsWbr4Z/u//4IEHYIcd4Mwz4TOfgX32Kdpp2izGyFe+8hXGjBnDvvvu2+R+3muby/e6tdf9Vs7mzJnD6NGjWbt2Ldtuuy233HILw4YNa3Rf77M6hVw377Pkhhtu4IknnuDxxx/Pa//m7jdgB2Bx0Yssc4awAoWw+XxyMXWNbrFeydChQxk6dOj7fx89ejSvvPIKP/7xj4vyD6v58+FXv4Lrr4e33oJx41Jr10knQTm+XeSLX/wis2fP5sEHH2xxX++1Ovlet1Lfb5Vg6NChPPnkk7z99tvcfPPNTJgwgalTpzYZKLzPkkKum/cZvPLKK5x33nlMnjyZrbbaKu/jmrrfqNLJ0yu2OzIL/fv3Z8mSJZutW7ZsGV26dGH77bfPqKrKM2rUKJ577rlWH792bZrPa+xYGDYMfv97+PSnYcECuO++9MRjOQawc845h9tvv53777+fgQMHNruv91qdQq5bY9p6v1Wabt26seeeezJy5EgmTZrEiBEj+NnPftbovt5ndQq5bo2ptvts5syZLFu2jAMPPJAuXbrQpUsXpk6dys9//nO6dOnCpk2btjimufsNeLN9Ki8vtoQVYPTo0dxxxx2brZs8eTIjR46ka9eq7M5ulVmzZrWqm+OVV9JUEldfDW+8Uf6tXrVijJxzzjnccsstTJkyhd3yeN+R91rrrltjWnu/dRQxRtatW9foNu+zpjV33RpTbffZkUceyZw5czZbd+aZZ7L33nvzjW98g86NDLxt7n6bPn16dQ5CzPrJgLYuQE9a+XTkqlWr4qxZs+KsWbMiEC+77LI4a9as+NJLL8UYY7zgggviGWecEWstXLgwbrPNNvHLX/5ynDdvXrzmmmti165d41/+8pdYLQq9Zj/5yU/iLbfcEp999tk4d+7ceMEFF0Qg3nzzzXmdr6Ymxvvui/Hkk9MTjj16xHjuuTE+80xJfl5JnH322bFXr15xypQpcfHixe8v77777vv7eK9tqTXXra33W6W78MIL47Rp0+KiRYvi7Nmz40UXXRQ7deoUJ0+eHGP0PmtKodet2u+zpjR8OrLA+y3zPJHFknkBbf4BbQhhtY8ZN1wmTJgQY4xxwoQJcezYsbG+KVOmxAMOOCB269YtDh48OF555ZWxmhR6zX7wgx/EPfbYI2611Vaxd+/eccyYMfHOO+9s8TyrVsV45ZUxDh+e7tJhw2L85S9jXLmyRD+shBq7XkC89tpr39/He21Lrblurb3fOopPf/rTcdddd43dunWLffv2jUceeeT7QSJG77OmFHrdqv0+a0rDEFbg/ZZ5nshiCTFW9li4EEJP4B2gV4xxZYGHV/aP74BefRV+8Ys02H7lSjjhBPjiF1PXY5WNE5akalKV/4R3TJjKwhNPwGWXpdcGbbMN/Od/pvC1665ZVyZJUmkYwpSZmhq4884UvqZMSYHrRz9Kc3v16JF1dZIklZYhTO1u7Vq47jr4yU/g2Wdh1Ci46SY48UTo4h0pSaoS/itP7WblyjTFxGWXpSkmTjophbEW3gwiSVKHVLEhLIQwkfQSbyecLXPLl8PPf54G3K9ZAxMmwPnnw5AhWVcmSVJ2fDpSJfPaa/C//5ve5whw1lnwla9AKyY9lyR1bD4dKRXDq6/C974H11yTnnT86lfh3HPTS7UlSVJiCFPRLF4Mkyallq9tt4VLL4UvfAF69sy6MkmSyo8hTG22dCn84Adw5ZWw1Vbw7W/DOecYviRJao4hTK22fDn88Idw+eXQtStccAF86UvQq1fWlUmSVP4MYSrY6tVpjq8f/jD9/atfTQPue/fOti5JkiqJIUx527AhDba/5BJYsSK9VujCCx1wL0lSazjHlloUI9x8M+y7bxpoP348LFiQpp8wgEmS1DqGMDXriSfgsMPgk5+E3XaDWbPg+uth8OCsK5MkqbIZwtSoZcvgc5+DkSNT1+PkyXD33TBiRNaVSZLUMTgmTJtZvz497XjppdC5c3rd0Oc/74u1JUkqNv/VqvfdfXeaYuK559Irhr77Xcd8SZJUKhXbHRlCmBhCmAc8lnUtle6119KYr2OOgQED0rivX/7SACZJUin5Au8qtmkTXHEFfPOb6R2PP/0pnHIKhKp8jaokKUNV+W+eim0JU9vMmgWHHJK6H//t3+CZZ+DUUw1gkiS1F0NYlVm3Dr71LTjooDT56vTpqetxu+2yrkySpOriwPwqMmMGnHlmavX69rfTux67dcu6KkmSqpMtYVVg3Tq4+GIYNSqFrpkzUwgzgEmSlB1bwjq4uXPhtNNS69ell8L550PXrllXJUmSbAnroGJMk66OHAk1Nakr8uKLDWCSJJULQ1gHtGwZfPzjcM456dVDjz8O+++fdVWSJKm+zEJYCGFwCOGaEMKiEMJ7IYQXQgiXhhAcqdQG99+fAteMGXDnnfCLX8DWW2ddlSRJaijLlrC9c+c/CxgOfBn4PPC9DGuqWDU18P3vw0c/CvvuC7Nnw7HHZl2VJElqSlnNmB9C+Dpwdoxx9wKOqfoZ899+GyZMgNtvT+O+al++LUlShajKqcLL7enIXsBbze0QQugOdK+3qkdJKypzTz0Fn/gEvPkm3HFHGgsmSZLKX9kMzA8h7AGcA1zVwq4Xklq+apdXS1xa2br1Vjj0UOjZM839ZQCTJKlyFD2EhRAuCSHEFpaRDY7ZCbgbuCnG+OsWTjGJ1GJWuwws9m8odzHCD38IJ5+cxn09+CDsnncHriRJKgdFHxMWQtgB2KGF3V6MMa7N7b8TcD/wKPAfMcaaAs9XVWPC1q+Hz38err0WvvnNNP6rU9m0Z0qS1CqOCSuGGONyYHk++4YQdiYFsJnAmYUGsGrz9ttw0knw8MNw/fVwxhlZVyRJklors4H5uRawKcDLwNeAviGkIBxjXJJVXeVq8WI4+mh49VX45z9hzJisK5IkSW2R5dOR44E9c0vDwfVV2SzZlIUL4WMfSy/ifuABGDYs64okSVJbldU8Ya3R0ceEzZkD48enJyAnT4Zdd826IkmSiq4qG18c0l3GZs+GceNgwIDUAmYAkySp4yi3yVqVM3s2HHFECl7/+Af07p11RZIkqZhsCStDBjBJkjo+Q1iZee65NAjfACZJUsdmCCsjr7+eBuH36QP33GMAkySpIzOElYkVK+Coo2DTpvQU5A4tvXNAkiRVtIodmB9CmAhMpAMEyQ0b4BOfSC1hDz0EgwZlXZEkSSo15wnLWIxw1llw3XVpDNhhh2VdkSRJ7a4q5wmr2JawjuJnP4Orr4bf/MYAJklSNan4rrxKdvfd8NWvwte/DmeemXU1kiSpPdkdmZGXX4YDDoBRo+D226Fz56wqkSQpc1XZHWkIy8CGDTB2LLz2GsyalaakkCSpilVlCHNMWAYuuggefzy9D9IAJklSdTKEtbO77oIf/xguuyx1RUqSpOpkd2Q7WrEC9t0X9t8f/v53CFXZ+CpJ0haq8t+IPh3Zjs47D9asSVNSGMAkSapudke2k9tug9/9Dn77Wxg4MOtqJElS1uyObAcrV8LQoXDQQSmM2QomSdJmqvLfjHZHtoNLLklB7IorDGCSJCmp2BAWQpgYQpgHPJZ1Lc2ZMwd+/nP49rd9MbckSapjd2QJxQiHHw5Ll8Ls2dCtWynPJklSxarKfiIH5pfQX/8K06bBPfcYwCRJ0uZsCSuRjRvTnGCDB6cXdUuSpCbZEqbiuf56WLAA/vjHrCuRJEnlyJawEli7FoYMgUMPhT//uRRnkCSpQ6nKlrCKfTqynF1zDbz+OvzXf2VdiSRJKle2hBXZxo2pFWzUKPjTn4r97ZIkdUhV2RLmmLAi+8tf4MUX05ORkiRJTbElrIhihA99CPr2hcmTi/nNkiR1aLaEqW3uuw+efBLuvTfrSiRJUrmzJayIPvUpmDcP5s71HZGSJBWgKv+t6dORRbJkCdx6K5x1lgFMkiS1rGJDWLm9wPu666BLFzjjjKwrkSRJlcDuyCKoqUnTUowZA7/9bTG+UZKkqlKVfUgOzC+Chx6ChQtTa5gkSVI+KrY7spzccAMMGgQf/nDWlUiSpEphCGujjRvhppvglFOgk1dTkiTlydjQRvffD2+8AaeemnUlkiSpkhjC2ujPf4Y990wz5UuSJOXLENYGNTXwt7/BySc7N5gkSSqMIawNnngCli6FY4/NuhJJklRpyiKEhRC6hxCeDCHEEMIHs64nX3feCb16waGHZl2JJEmqNGURwoAfAq9nXUSh7rwTjjoKunbNuhJJklRpMg9hIYRjgPHA17KupRBvvAGPP25XpCRJap1MZ8wPIfQDrgZOBN7N85juQPd6q3oUv7KWTZuWPo88MouzS5KkSpdZS1gIIQDXAVfFGGcUcOiFpHdF1i6vFr+6lk2dCrvvDgMHZnF2SZJU6YoewkIIl+QG2De3jATOAXoCkwo8xSSgV70lkxg0bRqMHZvFmSVJUkcQYozF/cIQdgB2aGG3F4EbgOOB+gV0BjYBf4gxTsjzfD1JLWK9YowrCyy3VT9+xQrYfnu49lqYkFeVkiSpGVU522bRx4TFGJcDy1vaL4RwLvDNeqt2Au4BTgEeLXZdxfTAAxAjHHZY1pVIkqRKldnA/Bjjy/X/HkJYnfvjCzHGTMZ55evRR2HAABg8OOtKJElSpcp8iopKNHMmHHigryqSJEmtVzYhLMb4YowxxBifzLqW5sQIM2akECZJktRaZRPCKsXLL8ObbxrCJElS2xjCCjRzZvo0hEmSpLYwhBVo5kzo3x922inrSiRJUiUzhBVozhwYMSLrKiRJUqUzhBVo/nwYPjzrKiRJUqUzhBVg7VpYuBD22SfrSiRJUqWr2BAWQpgYQpgHPNZe53z2WaipgWHD2uuMkiSpo6rYEBZjvCLGOAw4uL3OOW9e+rQlTJIktVXFhrAszJ8P/fpB795ZVyJJkiqdIawACxfCnntmXYUkSeoIDGEFWLQIdtst6yokSVJHYAgrwKJFMHhw1lVIkqSOwBCWp7Vr4fXXbQmTJEnFYQjL08svp09DmCRJKgZDWJ4WLUqfdkdKkqRiMITladEi6NwZBg3KuhJJktQRGMLy9OqrMGAAdOmSdSWSJKkjMITlacmSFMIkSZKKwRCWp8WLDWGSJKl4KjaEtfcLvBcvhv792+NMkiSpGlRsCGvvF3jbHSlJkoqpYkNYe9q0CZYutSVMkiQVjyEsD2+8ATU1toRJkqTiMYTlYcmS9GkIkyRJxWIIy8OyZemzb99s65AkSR2HISwPK1akz+23z7YOSZLUcRjC8vDWW+mVRT16ZF2JJEnqKAxheXjrLejTB0LIuhJJktRRGMLyUBvCJEmSisUQlocVK6B376yrkCRJHYkhLA+2hEmSpGIzhOXBECZJkoqtYkNYe77A2xAmSZKKrWJDWHu+wNsxYZIkqdgqNoS1p1WroGfPrKuQJEkdiSGsBTHC6tWw7bZZVyJJkjoSQ1gL1q5NQewDH8i6EkmS1JEYwlqwenX6tCVMkiQVkyGsBYYwSZJUCoawFqxZkz7tjpQkScVkCGuBLWGSJKkUMg9hIYTjQgiPhhDeCyEsDyH8Neua6jOESZKkUuiS5clDCJ8ArgYuAu4DArBfljU1ZHekJEkqhcxCWAihC/Az4OsxxmvqbVqQUUmNsiVMkiSVQpbdkR8CdgZqQgizQgiLQwh3hRCGN3dQCKF7CKFn7QL0KGWRq1dDp06w1ValPIskSao2WYaw3XOflwD/DXwcWAFMDSE097rsC4F36i2vlrBG1qxJXZEhlPIskiSp2hQ9hIUQLgkhxBaWkfXO/T8xxptjjDOBM4EIfKqZU0wCetVbBhb7N9RXG8IkSZKKqRRjwi4Hbmhhnxep60acV7syxrguhLAQ2KWpA2OM64B1tX8PJW6iWrvWrkhJklR8RQ9hMcblwPKW9gshzCSFqaHAg7l1XYHBwEvFrqu11q2D7t2zrkKSJHU0mT0dGWNcGUK4Crg0hPAKKXh9Pbf5pqzqamjdOujWLesqJElSR5PpPGGk0LUR+B2wNfAocESMcUWmVdWzfr0tYZIkqfgyDWExxg3A13JLWbI7UpIklULmry0qd4YwSZJUCoawFhjCJElSKRjCWmAIkyRJpWAIa4EhTJIklYIhrAXr1ztFhSRJKj5DWAtsCZMkSaVQsSEshDAxhDAPeKyU5zGESZKkUqjYEBZjvCLGOAw4uJTnMYRJkqRSqNgQ1l4MYZIkqRQMYS0whEmSpFIwhLVgwwbokvUbNiVJUodjCGvBpk2GMEmSVHyGsBZs2gSdO2ddhSRJ6mgMYS0whEmSpFIwhLXAECZJkkrBENYCQ5gkSSoFQ1gLNm40hEmSpOIzhDUjRqipMYRJkqTiM4Q1o6YmfTpFhSRJKraKDWHt8QLvTZvSpy1hkiSp2Co2hLXHC7wNYZIkqVQqNoS1B0OYJEkqFUNYMzZuTJ+GMEmSVGyGsGbYEiZJkkrFENYMQ5gkSSoVQ1gzakOYU1RIkqRiM4Q1w5YwSZJUKoawZhjCJElSqRjCmmEIkyRJpWIIa4ZTVEiSpFIxhDXDljBJklQqhrBmGMIkSVKpVGwI8wXekiSpklVsCGvPF3g7T5gkSSq2ig1h7cGWMEmSVCqGsGb4dKQkSSoVQ1gzbAmTJEmlYghrRk1N+uzkVZIkSUVmvGhGjOkzhGzrkCRJHY8hLA+GMEmSVGyGsGbUtoRJkiQVW6YhLISwVwjhthDC8hDCyhDCQyGEcVnWVJ/dkZIkqVSybgm7E+gCHAEcCDwJ/C2E0D/LohoyhEmSpGLLLISFEHYA9gS+H2OcHWN8DrgA2AYYnlVd9dkdKUmSSiXLF/K8CcwH/j2E8ASwDjgLWArMbOqgEEJ3oHu9VT1KVaDdkZIkqVQyC2ExxhhC+BhwG7AKqCEFsKNjjG83c+iFwHdKX2EdQ5gkSSq2ondHhhAuCSHEFpaRIYQA/BJYBnyE9CLu20hjwgY0c4pJQK96y8Bi/4ZadkdKkqRSKUVL2OXADS3s8yJpMP7Hgd4xxpW59V/ItY5NAL7f2IExxnWkrksAQgmbqeyOlCRJpVL0EBZjXA4sb2m/EMI2uT/WNNhUQ/ZPbW7GECZJkooty7AzHVgB/DaEMCI3Z9iPgN1IU1dkzu5ISZJUKpmFsFyL2dHAtsB9wAxgDHBCjPGprOqqz+5ISZJUKllOUUGMcQZwVJY15MMQJkmSiq2sxl6VG7sjJUlSqRjCmmF3pCRJKhVDWB4MYZIkqdgMYc2wO1KSJJWKIawZdkdKkqRSMYTlwRAmSZKKrWJDWAhhYghhHvBYqc5hd6QkSSqVig1hMcYrYozDSC/+LtE50qctYZIkqdgqNoS1J0OYJEkqNkNYM+yOlCRJpWIIa4bdkZIkqVQMYXkwhEmSpGIzhDXD7khJklQqhrBm2B0pSZJKxRCWB0OYJEkqNkNYM+yOlCRJpWIIa4bdkZIkqVQMYXkwhEmSpGIzhDXD7khJklQqFRvC2vMF3raESZKkYqvYENYeL/CuZQiTJEnFVrEhrD3YHSlJkkrFENYMuyMlSVKpGMLyYAiTJEnFZghrht2RkiSpVAxhzbA7UpIklYohLA+GMEmSVGyGsGbYHSlJkkrFENYMuyMlSVKpGMLyYAiTJEnFZghrht2RkiSpVAxhzbA7UpIklUrFhrD2eIF33blKfQZJklRtKjaEtccLvO2OlCRJpVKxIaw92B0pSZJKxRCWB0OYJEkqNkNYM+yOlCRJpWIIa4bdkZIkqVQMYXkwhEmSpGIzhDXD7khJklQqhrA82BImSZKKrWQhLIRwcQjh4RDCuyGEt5vYZ5cQwh0hhDUhhOUhhJ+HELqVqqZC2RImSZJKpUsJv7sbcBMwHfhMw40hhM7AncAbwBhge+C3QADOKWFdeXNgviRJKpWShbAY43cAQgj/0cQu44FhwKAY4+u5fb8KXBdCuDjGuLJUtRXKECZJkootyzFho4G5tQEs5x6gO3BgUweFELqHEHrWLkCPUhVod6QkSSqVLENYf2Bp/RUxxhXA+ty2plwIvFNvebVUBRrCJElSqRQUwkIIl4QQYgvLyAK+srGYE5pYX2sS0KveMrCA8xXMrkhJklQKhY4Juxy4oYV9Xszzu5YAh9RfEULoDXSlQQtZfTHGdcC6esfkebrC2RImSZJKpaAQFmNcDiwv0rmnAxeHEAbEGBfn1o0nBayZRTpHm8RoS5gkSSqNkj0dGULYBegD7AJ0DiF8MLfp+RjjamAyMA/4XQjh67l9fwxc7ZORkiSpowuxRH1uIYTrgAmNbBoXY5yS22cX4JfAEcB7wB+Br+W6HPM9T0/SAP1erQhvzf74tWthzRrYfvsCv1WSJBWiKps8ShbC2kspQ5gkSWoXVRnCfHekJElSBgxhkiRJGTCESZIkZcAQJkmSlIGSTVFRIapyIKAkScpeR3g6MpBe4r0qVvqPkSRJVaPiQ5gkSVIlckyYJElSBgxhkiRJGTCESZIkZcAQJkmSlAFDmCRJUgYMYZIkSRkwhEmSJGXg/wN9hOFauLhfNwAAAABJRU5ErkJggg==\n", "text/plain": [ "Graphics object consisting of 1 graphics primitive" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r_0 = 2\n", "rmin = 1.0001\n", "plot(lambda r: zz(r_0, r), (r, rmin, 4), \n", " axes_labels = [r'$r/m$', r'$Z(r)/m$'])" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "rhof = fast_callable(rho(r, mu0), vars=[r])" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAGSCAYAAACrGxZhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA1eUlEQVR4nO3deZyd4/n48c9lS1skqmqvb1VRVWsljTWJtsQWIpYIYo8lal9b1f6+35b221JVWymKUlWK2qktQRaJWKNUUWIJjUiCZLK4f3/cJ/2OMTOZM3NmnnPmfN6v13mdmec8zznXPJ6Ma677fq47UkpIkiSpay1WdACSJEn1yCRMkiSpACZhkiRJBTAJkyRJKoBJmCRJUgFMwiRJkgpgEiZJklQAkzBJkqQCmIRJkiQVwCRMkiSpACZhkiRJBSg0CYuIFSNilTKPOTkidu7g5+7ZkeMlSZI6qrAkLCK+DFwIfFDmoecCR0bEwHZ+7hbAAe05VpIkqVKW6MjBEbE78CPgy0BPYA7wPPBxaZflgKWAe4GzUkovl45bFvgDMDSlNKucz0wpLYiI/YHRETEppTS1zLD3Ba4t8xhJkqSKipRSx98k4lTgZ8ApKaVfNHltHeAWYGVgy5TS8xFxGTA+pXRZBz7zh8DaKaXhZRyzBPAPYP2U0kft/WxJkqSOqtRw5Nal57uavpBSehG4FPg8cHJEbAIMAK7s4GdeAOwSEV8t45iBwCMmYJIkqWgdTsIiYnFyEvZWSunZFnb7XOn588ApwJUppfkd+dyU0nTgAeCQMg5zKFKSJFWFSlTCvkmeD/a3Vvb5bun5QWA34M4KfC7AKGBwW3aMiGWAvrQepyRJUpeoRBI2oPTcbHITEdsA/YFHgLeBBcDTzey3UURcExGPRMTOEfH5iDgvIi6OiJsjYuNm3n4UsG5ErNSGOHcHbmtagYuIrSLi6ogYExE7RcRiETEyIn4TEZeU4ulb2ndYKZ6LS9u3bcPnSpIkfUqnJWER0aN0F+ONwG/J1bA+wOMppQXNvM+xwEHA3cAVwGXAL4AbgJ1oftjx1dLzRm2I81NDkRGxGDCC3LLiUfI8tV8Dk1NK30spHQE8B1wXEWcCc1JKR6aUjiQPhf659B6SJEll6WiLiiWALYEPgXMjYuFLi5feeyKwRUrppdL+6wFvNPM+awFvppTmR8SqwPLklhZvRMTm5F5itzUTwvvkytqXFxHnSsB/pZTGNXmpNzAppZRKn/tF4I6U0oON9pkJrAlMTSn9pdH2qaU4v1j6WpIkqc06lISRk5hlgD+mlIa1Yf/VgFea2b4ScFPp663IdzA+AZBSupFcTfuUUvI0A+i1iM8dCvypme09gJtLX28N3JNSurvJPhuWYr60yfb1gNnAtEV8tiRJ0qd0dCht4VDk6Dbuvwwwo+nGlNJjKaVJEbEC8A3goTJiaCAnU63ZF7iumc8dlVJ6tdTmYnXgvsavlyp9WwAPpU83VNsOGN3RuzwlSVJ96uokbAG5g35r7xeUl4QtTytLH5WaxUZK6YVW3mPhBPumn7uw0veJ7RGxAbA2LVToJEmSFqXdSVhELEWuEr1HnrzeFu+Tk6aWDCBXtsa2MYalyVWwT80za6QtvcEGlGKb1GR7/9LzQ022DyPHeWMpjnJ6lUmSJHVoTti3yE1Y/9bMUF1LXmXRSdjYlNKcNr7fwtYUrVW59ub/kqmW9AdGpZQ+brJ9APBySum1Jtt3B/6aUppeqoqt3rZwJUmSso4MR+5Qeh5TxjFPk+80/JSIWBn4GuUNRW5Ovnux2U79pf5er6WU3m7pDUp3bK7c9HMjYklK88GaOWwF4OFSe4pTyUsoSZIktVlZSVhELBkR4yLin8Bppc2nRcRTEbFXG95iFLBBRCzXzGsrkVs93FBGSNsAdzZTwVqoLUORXwDeBP7SZPty5NYbf2jmmGOAPch3XF6ZUvIOSUmSVJZo+0hiBT4sV47eBA5NKd1egfd6CTg6pfSpZZBKdza+BGyQUprVkc+SJEmqtC7t9l6qWP2O3Lero/YgV6ruauH17YBxJmCSJKkaFbHkzq+BgRHR7NywMpwGHNvKTQFtGYqUJEkqRJcOR/7nQyOOAPqklA5u5/GnA19IKZ3UwutLA5OBr6aU5rU/UkmSpM5RyOLTKaVLgMUj4oByj42I7YF1yXcltmQwecK+CZgkSapKhSRhJYcAfUstIsrxFnBQSmlBK/v0A65ud2SSJEmdrJDhSEmSpHpXZCVMkiSpbpmESZIkFcAkTJIkqQAmYZIkSQUwCZMkSSqASZgkSVIBTMIkSZIKYBImSZJUAJMwSZKkApiESZIkFcAkTJIkqQA1m4RF1jMiouhYJEmSyrVE0QF0wLLADKAXMLPMY121XJKk6lC3xZSarYRJkiTVMpMwSZKkAtRcEhYRIyNiMjC+6FgkSZLaK1KqzelREdGT0pywlJJzwiRJqiJvvw3Tp8N66y1yV+eESZIkVcKECdC7Nxx2GNRoradLmIRJkqSKue462HprWHVVuOEGsJFUy0zCJElShy1YAKeeCvvuC3vtBQ8/nBMxtayW+4RJkqQqMGMGDBsGd98N55wDxx9vBawtTMIkSVK7vfgiDBoEU6fCnXfC9tsXHVHtcDhSkiS1y913Q58+sNhiMH68CVi5TMIkSVJZUoJf/AJ22ilPwh87FtZeu+ioao9JmCRJarPZs2H4cDjllDwR/5ZboGfPoqOqTc4JkyRJbfLGGzB4MDz7LFx/Pey9d9ER1TaTMEmStEhjx+YEbMkl4ZFHYNNNi46o9jkcKUmSWnXlldCvH6y1Fjz+uAlYpZiESZKkZs2bB8ccAwcfDAcdBA88ACutVHRU3UfNDUdGxEhgJCaQkiR1mn//O3e+Hz0aLr4Yjjii6Ii6n0g1urJmRPQEZgC9Ukozyzy8Nn9oSZK6wFNPwW67wYcfwk035TYUnahue+tbTZIkSf/x5z/DFlvA8svDhAmdnoDVNZMwSZLExx/DD36QhyB33TUPQ66xRtFRdW81NydMkiRV1owZsN9+cMcd8L//Cyed5ALcXcEkTJKkOrZwAe63384LcA8cWHRE9cPhSEmS6tRdd/3fAtyPP24C1tVMwiRJqjMpwc9/nhfg3mYbF+AuikmYJEl15KOPYNgwOO20PBHfBbiLU/EkLCJOj4gUEectYr9+ETExIuZExMsRYRs4SZI60b/+BVttBX/9a25F8T//k4ciVYyKnvqI6A2MAJ5exH5rAncCo4FNgLOA8yNiSCXjkSRJ2cMPw2abwfTpMGYM7LFH0RGpYklYRCwDXAscBkxfxO5HAK+llI5LKT2fUvodcAVwUqXikSRJef7XRRfBd74DG26YJ+BvuGHRUQkqWwm7ELgjpfS3Nuy7OXBvk233AJtFxJLNHRARPSKi58IHsGzHwpUkqXtraIARI2DkyPy45x5YYYWio9JCFekTFhFDgU2B3m08ZGVgapNtU0vxrAC81cwxpwM/am+MkiTVk7ffhiFD8tJDV14JBx5YdERqqsNJWER8Cfg1sF1KaU4ZhzZdRDta2L7Q2cC5jb5fFphSxudJklQXJkzIC3B//HGeC9a3b9ERqTmVGI78JrAiMDEi5kfEfKAfcEzp+8WbOeZtcjWssRWB+cC05j4kpdSQUpq58AHMqkDskiR1K9dck++AXH31nIyZgFWvSiRh9wMbABs3ekwgT9LfOKW0oJljxgDfbbJtO2BCSmleBWKSJKmuzJ8PJ54Iw4fnPmAPPQSrrlp0VGpNh4cjU0qzgGcbb4uID4FpKaVnS9+fDayWUhpe2uUS4OiIOBe4jDxR/xBgn47GI0lSvXnvPRg6FB54AM4/H44+2gW4a0FXLeC9CrDGwm9SSq9ExI7Ar4CRwJvAMSmlm7ooHkmSuoVnn4Vdd4UZM+Dee2HbbYuOSG0VKbU0D766ldpUzAB6leaIlaM2f2hJkhq5+WbYf39Ya628/NCaaxYdUbvUbc3OxQokSaoxH38MP/4x7L477LADPPpozSZgda2rhiMlSVIFzJqVJ9/feiv85Cfw/e87/6tWmYRJklQj/vnPPP/rtddyErbLLkVHpI5wOFKSpBpw333QuzfMnQvjxpmAdQcmYZIkVbGU4NxzYeBA+Na3YPx4WG+9oqNSJZiESZJUpWbPzvO/TjwRTjoJbr8dlluu6KhUKc4JkySpCk2ZAoMH5z5g110H+9jOvNupuSQsIkaSG7xaxZMkdUujR8Mee0CPHrn9xKabFh2ROoPNWiVJqhIpwcUXw7HHwpZbwg03wIorFh1Vp6vbBhtWkyRJqgINDXDYYTByJBx1VL4bsg4SsLpWc8ORkiR1N2+8AUOGwJNPwu9/DwccUHRE6gomYZIkFejRR3MCtuSSeS5Y795FR6Su4nCkJEkF+e1vYcAAWGcdmDDBBKzemIRJktTFGhpgxAg44gg4/HC4/35YaaWio1JXczhSkqQu9Oabuf3ExIlwxRVw0EFFR6SimIRJktRFxozJ878WWyzP/+rTp+iIVCSHIyVJ6gKXXQb9+sFaa+X5XyZgMgmTJKkTzZ2b536NGAGHHprnf628ctFRqRo4HClJUid5++08/+vxx3Ml7NBDi45I1cQkTJKkTjBuHOy+e16K6OGHoW/foiNStam54ciIGBkRk4HxRcciSVJzLr8cttkGvvzlfBekCZia4wLekiRVyNy5cPzxcNFFuf/X+efDUksVHVXVq9sFvB2OlCSpAqZOzfO/xo3LnfBHjCg6IlU7kzBJkjpo/Pg8/2vBAnjoIdhii6IjUi2ouTlhkiRVkyuvzPO/vvSlPP/LBExtZRImSVI7zJsH3/seHHww7L9/roCtumrRUamWOBwpSVKZ3nkH9twTHnsMLr44T8KPup1ervYyCZMkqQwTJ8LgwflOyAcfhK22Kjoi1SqHIyVJaqOrr4Ytt4RVVsnJmAmYOsIkTJKkRZg3D447Dg44APbdN3fAX221oqNSrXM4UpKkVrz7Luy1FzzyCFxwARx1lPO/VBkmYZIkteCJJ/L8rzlz4P77cysKqVIcjpQkqRnXXpvnf624IkyYYAKmyjMJkySpkXnz4IQTYL/9YO+9YfTo3IhVqrSaG46MiJHASEwgJUkV9s47ef7Xo4/mxbePPtr5X+o8kVIqOoZ2iYiewAygV0ppZpmH1+YPLUnqNOPHw5AhuRL25z/D1lsXHVHdqNs012qSJKnuXX55TrpWXz33/zIBU1cwCZMk1a2Ghrzk0KGH5jUgH3rI/l/qOjU3J0ySpEp44w3YY4/chuLyy3MSJnUlkzBJUt0ZPTovwL3kkrkJa+/eRUekeuRwpCSpbqSU73rcdltYb708/8sETEUxCZMk1YWPPoLhw+HYY+GYY+C++3IjVqkoDkdKkrq9V16B3XeHF16A666DffYpOiLJJEyS1M3de29OupZbDsaOhQ03LDoiKXM4UpLULaUEP/sZ7LAD9OkDjz9uAqbqYhImSep2Zs3Kdz+efjp8//tw++2w/PJFRyV9ksORkqRu5YUXYPBgmDIFbr4Zdtut6Iik5tVcJSwiRkbEZGB80bFIkqrLX/+ahx5TymtBmoCpmtVcEpZSujCl9HWgT9GxSJKqw8cfw5lnwq67wne+kxOwr32t6Kik1jkcKUmqadOnw377wV13wVlnwWmnQUTRUUmLZhImSapZzzyT53+9915OwrbfvuiIpLarueFISZIArr8e+vaFZZbJyw+ZgKnWmIRJkmrK/Plw0km5AevgwfDYY7DmmkVHJZXP4UhJUs14913Ye28YNQrOOy+vAen8L9UqkzBJUk2YMCGv/9jQAPffD/36FR2R1DEOR0qSqt4VV8BWW8Gqq+b5XyZg6g5MwiRJVWvuXDjySDjkEBg+HB5+GFZfveiopMrocBIWEUdGxNMRMbP0GBMRO7Syf/+ISM08bKsnSfqPN96A/v1zFezSS/OjR4+io5IqpxJzwqYApwEvlb4/ALg1IjZJKT3XynHrAjMbff9uBWKRJHUDDz2UJ+AvuWSehP+tbxUdkVR5Ha6EpZRuSyndmVJ6sfT4AfAB0HcRh76TUnq70WNBR2ORJNW2lOCcc/LSQ+uvD088YQKm7quic8IiYvGIGAosDYxZxO6TIuKtiLg/Iga04b17RETPhQ9g2UrELEmqDrNmwV575R5gJ50E994LK65YdFRS56lIi4qI2ICcdH2GXAUbnFKa3MLubwEjgIlAD2B/4P6I6J9SGtXKx5wO/KgS8UqSqsvzz+f2E2+8ATfdlL+WurtIKXX8TSKWAtYAlgOGAIcC/VpJxJoefxuQUkqDWtmnBzlpW2hZ8ny0Ximlmc0f1aKO/9CSpIq48UY46CBYYw34y19g3XWLjkhdrG7b7VZkODKlNDel9FJKaUJK6XTgKeDYMt5iLLD2Ij6jIaU0c+EDmNWBkCVJBZs/H04+GfbcE3baCcaNMwFTfemsjvnBJ6tWi7IJeZhSklQHpk6FoUNh9Gj41a/g2GNdfkj1p8NJWEScBdwFvE4eIhwK9AcGll4/G1gtpTS89P1xwKvAc8BSwH7kIcwhHY1FklT9xoyBPfaABQvgwQdh662LjkgqRiUqYSsB1wCrADOAp4GBKaX7Sq+vQp4vttBSwC+B1YDZ5GRsp5TSnRWIRZJUpVKCCy+EE06APn3ghhvyMkRSvarIxPwilNpUzMCJ+ZJU9T76CA4/HP7whzz0+Itf5EasEnU8Mb+z5oRJkgTASy/BkCH5+brrYJ99io5Iqg4u4C1J6jS33QabbQazZ+e7H03ApP9jEiZJqrgFC+CMM2DQoLwI9+OPwze+UXRUUnVxOFKSVFHTpsGwYfC3v8HZZ8Mpp8Bi/skvfYpJmCSpYiZMyO0nPvwwr/347W8XHZFUvfzbRJJUEb/7HWy5ZV50e+JEEzBpUWouCYuIkRExGRhfdCySJJgzBw49FA47LK8BOXp0XgdSUuvsEyZJardXX83Dj889BxdfDAceWHREqkH2CZMkqRz33JMn4PfqBY89BptsUnREUm2pueFISVKxPv4YfvIT2GEH+Na38mR8EzCpfFbCJElt9v77sP/+cMcdcOaZ+WH7Cal9TMIkSW3y1FN5+aFp0+D222HHHYuOSKpt/v0iSVqk3/8e+vaFZZfN7SdMwKSOMwmTJLVo9uz/az2x7755Av5XvlJ0VFL34HCkJKlZL7+c2088/zxccUVOxCRVjkmYJOlTbrsNhg+HL3wBxoyBjTcuOiKp+3E4UpL0H/Pnw+mnw6BB0K9fbj9hAiZ1DithkiQApk6FffaBUaPg5z+Hk0+GqNte5lLnMwmTJPHII7DXXrkR6/335yqYpM7lcKQk1bGU4NxzoX9/+OpXYdIkEzCpq9RcEhYRIyNiMjC+6FgkqZbNnAl77gknnggnnAAPPACrrFJ0VFL9iJRS0TG0S0T0BGYAvVJKM8s8vDZ/aEmqkGeeyd3vp07NjVgHDy46ItWxup15WHOVMElSx1xzTV54+7OfzXc/moBJxTAJk6Q6MWcOHHFE7v+19965/9faaxcdlVS/vDtSkurAq6/m7vfPPguXXQaHHGL7CaloJmGS1M3deSfstx8st1xe+3HTTYuOSBI4HClJ3daCBfDDH8JOO8FWW8HEiSZgUjWxEiZJ3dC778KwYbntxFlnwamnwmL+2S1VFZMwSepmHnssd7+fNw/uuw+23bboiCQ1x7+LJKmbSAnOOSd3vP/yl+GJJ0zApGpmEiZJ3cD06bnf10kn5e73Dz4Iq61WdFSSWuNwpCTVuAkT8vDj9Olw660waFDREUlqCythklSjUoKLLoItt4QVVsiLb5uASbWj5pIwF/CWJJg1C/bZB0aOhMMPh9Gj8zwwSbXDBbwlqcY8/TTsuSe89RZcfnn+Wqphdbt2Q81VwiSpnl15ZV58+zOfyXPBTMCk2mUSJkk14KOP4KCD4OCD8xJEY8fCOusUHZWkjvDuSEmqcn//e654vfwyXHUVDB9edESSKsFKmCRVseuug802y+tAjh9vAiZ1JyZhklSF5syBI4+EffeF3XbLCdj66xcdlaRKcjhSkqrMP/+Zhx8nT4bf/hYOOwyibu8fk7ovkzBJqiI335wn4K+wAowZA5tsUnREkjqLw5GSVAXmzs1rPu6+O3z72zBxogmY1N1ZCZOkgr32Guy9d+77dd55cMwxDj9K9cAkTJIKdOedsP/+sPTSeemhvn2LjkhSV3E4UpIKMH8+nH467LRTTrwmTTIBk+qNlTBJ6mKvv54X3x47Fs4+G045BRbzT2Kp7tRcEhYRI4GRWMWTVIPuuCM3XP3c5+Dhh2HLLYuOSFJRai6RSSldmFL6OtCn6Fgkqa3mzYOTT4add4YttoAnnzQBk+pdzVXCJKnW/OtfMHRovvvxl7/MrSi8+1GSSZgkdaJbb4UDD4Revbz7UdIn1dxwpCTVgrlz4bjj8rqP/ft796OkT7MSJkkV9vLLufnqU0/Br38N3/uew4+SPs0kTJIq6Kab4OCD89qPjz0Gm21WdESSqpXDkZJUAXPmwNFHwx57wHbbwRNPmIBJap2VMEnqoJdegr32gueegwsvhCOPdPhR0qJ1uBIWEUdGxNMRMbP0GBMROyzimH4RMTEi5kTEyxFxREfjkKQiXH89bLopzJqVO+AfdZQJmKS2qcRw5BTgNGCz0uMB4NaIWL+5nSNiTeBOYDSwCXAWcH5EDKlALJLUJWbPhiOOyMsP7bxzHn7cZJOio5JUSyKlVPk3jXgPODmldHkzr/0cGJRSWq/RtkuAjVJKm5fxGT2BGUCvlNLMMkOs/A8tqW688EIefnzxRTj/fDj0UKtfUgfU7b+eik7Mj4jFI2IosDQwpoXdNgfubbLtHmCziFiylffuERE9Fz6AZSsStCSV4Q9/gG9+ExoaYNw4OOwwEzBJ7VORJCwiNoiID4AG4BJgcEppcgu7rwxMbbJtKvkmgRVa+ZjTyZWvhY8pHQpaksrw0UdwyCGw//6w++55CaINNyw6Kkm1rFKVsBeAjYG+wMXAVRHx9Vb2bzocGC1sb+xsoFejx+rtilSSyvTMM9C7N/zxj3DllXD11bDMMkVHJanWVSQJSynNTSm9lFKakFI6HXgKOLaF3d8mV8MaWxGYD0xr5TMaUkozFz6AWZWIXZJakhJcfDH06QOLL56rXwceWHRUkrqLzmrWGkCPFl4bA3y3ybbtgAkppXmdFI8klWX69Nx49aijcgf8cePg663V9yWpTB1u1hoRZwF3Aa+TJ8sPBfoDA0uvnw2sllIaXjrkEuDoiDgXuIw8Uf8QYJ+OxiJJlfDIIzBsWO799Ze/wODBRUckqTuqRCVsJeAa8ryw+4FvAQNTSveVXl8FWGPhzimlV4AdyYnak8APgWNSSjdVIBZJarcFC+AnP4F+/WCNNfIC3CZgkjpLp/QJ6wr2CZNUSW+8AfvtBw8/DGecAWeeCUu4sJvUFeq2yYu/YiTVvdtvzxPue/SA+++HAQOKjkhSPeisifmSVPUaGuC442CXXWDzzfPwowmYpK5iJUxSXXrxRRg6FJ57Ds47D445xs73krqWlTBJdefqq2HTTeGDD2DMGDj2WBMwSV3PJExS3Zg1Ky87dMABMGQITJyYkzFJKoLDkZLqwsSJefjx7bfhmmvynZCSVCQrYZK6tZTgV7/KE+979oQnnjABk1Qdai4Ji4iRETEZGF90LJKq27vv5jsfTzgBjj4aHnsM1l676KgkKbNZq6Ru6d5789yv+fPhqqtgxx2LjkhSC+r2tpiaq4RJUmsaGuDEE2H77WGDDeDpp03AJFUnJ+ZL6jb+/nfYZ5/c++ucc3Ij1sX8U1NSlfLXk6SalxJcdlluNzF7Nowbl+eBmYBJqmb+ipJU06ZNyz2/RozIPcAmToRNNik6KklaNIcjJdWsBx/Midfs2fCXv8DgwUVHJEltZyVMUs2ZNw9OPx2+/W1YZ5288LYJmKRaYyVMUk156SUYNgwmTYKzz4aTToLFFy86Kkkqn0mYpJqQUu73dfTRsMoqufFq795FRyVJ7edwpKSq9/77ed3Hgw6CPffMSw+ZgEmqdVbCJFW10aPzWo8zZsD118PeexcdkSRVhpUwSVVp3jw480zo3x/WWCNPvjcBk9Sd1FwlLCJGAiMxgZS6rRdfzNWvJ56AH/843wm5RM39tpKk1rmAt6SqsbDz/fHHw6qrwrXXQp8+RUclqZO5gLckFemdd2DXXeHww3MVbNIkEzBJ3ZsFfkmFu/12OOSQXAn7619hl12KjkiSOp+VMEmF+fBDOPLInHT17g3PPGMCJql+WAmTVIjHH8/Djq+/DhdfnIcho25nhkiqR1bCJHWp+fPhJz+BLbaAnj3z3K8jjjABk1R/rIRJ6jIvvwz77w9jx8L3v5/7gC25ZNFRSVIxTMIkdbqF6z5+73vwxS/CqFGw5ZZFRyVJxXI4UlKnmjYtr/e4cN3HJ580AZMksBImqRPdc09Ovhoa4MYbYciQoiOSpOphJUxSxX3wQW49MXAgbLBBbj1hAiZJn2QlTFJFPfooHHAAvPUWXHSRdz5KUkushEmqiIYGOO002GYbWHFFeOqpXA0zAZOk5tVcJSwiRgIjMYGUqsZTT+XWE3//O/z0p3DyybD44kVHJUnVreYSmZTShSmlrwMu7SsVbP58OPvsvORQRO6Cf9ppJmCS1BY1l4RJqg7/+EceejzjDDjxRBg/HjbaqOioJKl2mIRJKktKecL9xhvDO+/A6NG5GtajR9GRSVJtMQmT1GZTpuS2EyNH5jsgn3wyrwEpSSpfzU3Ml9T1UoLrroOjj4bPfQ7uvhu2377oqCSptlkJk9Sqd9+FvfaC/faDHXeEZ581AZOkSrASJqlFN92Ue30tWAB/+lNOxiRJlWElTNKn/PvfMHQo7LEHbLUVTJ5sAiZJlWYlTNIn3HxzXmpo3jy49lrYZx+73ktSZ7ASJgmAadNg331h991h881z9WvYMBMwSeosVsIkceutcPjhMHcuXHNNTsZMviSpc1kJk+rYe+/lNR932w369IHnnst3QZqASVLnq7lKmAt4S5Vx220wYgTMmQNXXZWTMZMvSeo6kVIqOoZ2iYiewAygV0ppZpmH1+YPLVXA9Olw3HFw9dWw007w29/CaqsVHZWkOla3f/7VXCVMUvvdcUeufn34IVx5ZV56yOqXJBXDIT2pDkybBsOHw847w0Yb5a73Bx5oAiZJRbISJnVjKcGNN+Y1H+fOtfolSdXESpjUTb31Vu75tddeuev9889b/ZKkamIlTOpmUoLf/x5OOAF69MiVsCFDio5KktSUlTCpG3n1Vdh+ezj4YBg0KHe9NwGTpOpkEiZ1AwsWwPnnwze+AX//O9x5Z+79tfzyRUcmSWqJSZhU455/HrbZBo49Ns/5eu452GGHoqOSJC1Kh5OwiDg9Ih6PiFkR8U5E3BIR6y7imP4RkZp5fK2j8Uj1Yt48+OlPYeON4d13YdQouOACWHbZoiOTJLVFJSph/YALgb7Ad8mT/e+NiKXbcOy6wCqNHv+oQDxSt/fEE3mtxx/9CI4/Hp56CrbeuuioJEnl6PDdkSmlgY2/j4iDgHeAbwKjFnH4Oyml9zsag1QvPvwQ/t//g3PPhfXXh3Hj4JvfLDoqSVJ7dMacsF6l5/fasO+kiHgrIu6PiAGt7RgRPSKi58IH4KCL6srdd+eJ9+efnxOxxx83AZOkWlbRJCwiAjgXeCSl9Gwru74FjACGALsDLwD3R8Q2rRxzOnnB7oWPKRUJWqpyU6fCsGF5sv1XvgLPPAM/+AEstVTRkUmSOiJSSpV7s4gLgZ2ArVJKZSVJEXEbkFJKg1p4vQfQo9GmZcmJWK+U0swyQ63cDy11kpTgiivg5JNhscXyEOT++9vxXlK3U7e/1SpWCYuI3wCDgAHlJmAlY4G1W3oxpdSQUpq58AHMameoUtV74QUYMAAOPRR22SX3/ho+3ARMkrqTSrSoiIi4gDysuG1K6ZV2vtUm5GFKqW41NMB//zdsuCFMmQL33Zebrq6wQtGRSZIqrRJrR14IDAN2BWZFxMql7TNSSrMBIuJsYLWU0vDS98cBrwLPAUsB+5Hnh7nAiurW6NEwYgS89BKccgqccQZ89rNFRyVJ6iyVSMKOLD0/1GT7QcDvS1+vAqzR6LWlgF8CqwGzycnYTimlOysQj1RT3n8fTj0VLr0U+vbNPcA22KDoqCRJna2iE/O7UqlNxQycmK8alRJcey2ceCLMmQNnnw1HHJEn4UtSHanb2a7+upcKMHlynni///75+fnn4aijTMAkqZ74K1/qQh9+CKedBhttBG++CffeC9dfD6uuWnRkkqSuVok5YZIWISW49VY49lh455285uPJJ0OPHos+VpLUPVkJkzrZyy/nXl+DB+dlh557Lt/5aAImSfXNJEzqJA0N8JOf5IW2n34abr4Zbr89Lz0kSZLDkVIn+NvfYOTIXAU78UT44Q9h6aWLjkqSVE1qrhIWESMjYjIwvuhYpKbeeAOGDoXvfhdWWQWeegp+9jMTMEnSp9knTKqAhoa8wPZPf5oTrl/+Evbbz7UeJakN6vY3pcORUgeklOd5HX88/OtfcMwxcOaZ0KtX0ZFJkqpdzQ1HStXihRdgxx1h0KA82f7pp+Gcc0zAJEltYxImlWnmzLzA9gYb5ETsllvgnntgvfWKjkySVEscjpTa6OOP4Q9/yIttz5iRhx1POgk+85miI5Mk1SIrYVIbTJgAW24JBxwA/frlCtgZZ5iASZLazyRMasXUqXDoodCnT1738aGH8lqPX/pS0ZFJkmqdw5FSM+bMgfPOg7POgiWWgN/8Bg4/PH8tSVIl+L8UqZGU4E9/gtNOy41XR47Mc7+WX77oyCRJ3Y1JmFQydmzu9zV2LOy6K9x7L6yzTtFRSZK6K+eEqe79618wbBhsvnkehrz//tx2wgRMktSZTMJUt2bOhO9/H9ZdFx58EK64It8Fue22RUcmSaoHDkeq7sybB7/7Hfz4xzBrVm68esopsMwyRUcmSaonNVcJi4iRETEZGF90LKotKcGNN8L66+cJ9wMH5n5f//3fJmCSpK4XKaWiY2iXiOgJzAB6pZRmlnl4bf7QardRo3K1a9w42GEH+NnPYMMNi45KkgRE0QEUpeYqYVI5nn0Wdtkld7lfsAAeeADuvNMETJJUPJMwdUtTpsDBB8NGG8HkybnL/bhxMGBA0ZFJkpQ5MV/dyrvv5qHGiy6CZZeFX/8aRoyApZYqOjJJkj7JJEzdwowZcM458KtfQQSceiqccAL07Fl0ZJIkNc8kTDXtww/zuo7/+78wezZ873t5Av4KKxQdmSRJrTMJU01qaIBLL4Wf/hTeew8OOwx+8ANYddWiI5MkqW2cmK+aMndubrS6zjpw3HG53cSLL8KFF5qASZJqi5Uw1YS5c+HKK+Hss/Naj3vuCffcA1/7WtGRSZLUPiZhqmoNDXD55fmOxylTYK+94Pbb4RvfKDoySZI6xuFIVaU5c+CCC2CttfJk+623zo1Xr7/eBEyS1D2YhKmqzJ6de3t95Stw7LG5uerkyXDttfD1rxcdnSRJlVNzw5ERMRIYiQlktzJrVr7b8Re/gH//G/bbL9/tuPbaRUcmSVLncAFvFeqdd+D88/PdjR98APvvn5OvtdYqOjJJUhep2wW8a64Spu7h5Zfhl7/MdzwuvnheWuj44+FLXyo6MkmSuoZJmLrUpEnw85/Dn/8MX/hCrnoddRQsv3zRkUmS1LVMwtTpUoIHHsjJ1333wZpr5qWGDjoIPvvZoqOTJKkYJmHqNA0N8Kc/5bsdn3gCNt4Y/vhH2GMPWMIrT5JU5/xfoSrunXfgkkvgootg6lQYODB3t//udyHqdvqlJEmfZBKminnqqVz1uvbaXOk64AA45hiXFpIkqTkmYeqQuXPh5pvh4ovh4Ydh9dXhf/4HDj3UyfaSJLXGJEzt8uqrubnq5Zfn4cdttsnzvwYPhiWXLDo6SZKqn0mY2mzBArjrrlz1uusuWHbZPOR4xBEuKSRJUrlMwrRIb72VK16XXgqvvw7f/CZcdhkMHQpLL110dJIk1SaTMDVr7ly4/fbc0f6uu2CppWCffeDII2GzzYqOTpKk2mcSpk948smceF17LUybBr1758aq++wDyy1XdHSSJHUfNZeERcRIYCSwWNGxdBdvvAE33ABXXZXbTKy0Uu5mf+CBsP76RUcnSVL3FCmlomNol4joCcwAeqWUZpZ5eG3+0BU0bRrceGPuYD9qVL6jceedc/I1cKAd7SVJXaZu23j7v9o6MmsW3HJLTrzuuy+v6fjtb8MVV8BuuzncKElSVzIJ6+amTcsT7G+5Be6+G+bMgS23hPPOgz33hBVXLDpCSZLqk8OR3dBrr8Gtt+ZO9qNG5f5em2+eG6nuvTessUbREUqS9B8OR6p2zZsHY8fmStddd8GkSXmO17bbwoUXwqBBsMoqRUcpSZIasxJWo15/He65Jyddf/sbzJwJK6wA228PO+4IO+0EvXoVHaUkSYtUt5Uwk7Aa8e67eWjxoYfgwQfhuedgscWgb998N+MOO8Cmm+ZtkiTVEJOwWtOdk7CU4JVX8hDjY4/Bww/Ds8/m1776VejXL1e8vvMd+Pzni41VkqQOMgmrNd0pCZs2LXeqHzcuJ15jx+bKF8Daa+ekq3///Lz66kVGKklSxZmEtfsNIk4Hdge+BswGHgNOTSm9sIjj+gHnAusDbwL/m1K6pIzPrbkkbN48eOEFePrp3Jn+6afz48038+s9e8K3vpWHGPv2zV9/4QtFRCpJUpcxCWv3G0TcDVwPPE6+2/KnwAbA11NKH7ZwzJrAs8BlwG+BLYGLgH1SSje18XOrLglLKVe13ngDpkyBV1+Ff/wDXnopP7/ySk7EILeJ2HDD/Nhoo/y8zjrO6ZIk1R2TsIq9YcQXgXeAfimlUS3s83NgUEppvUbbLgE2Silt3sIxPYAejTYtC0yhE5KwF1/MydTcudDQkB8Lv549G6ZPh/fey4/p02Hq1Jx0TZmS91loySXhK1/J87jWXjs/NtggP+xOL0kSUMdJWGf0CVvYGOG9VvbZHLi3ybZ7gEMiYsmU0rxmjjkd+FEF4lukkSNz24eW9OqVJ8Qvv3x+rLoq9OmT52utvjqstlp+Xmkl12CUJEnNq2glLCICuBX4fEpp61b2exH4fUrprEbbtgAeBVZNKb3VzDFdVgn75z/z8j49esBSS33yuUcPWHzxMj9NkiS1xEpYhVwAbAhs1YZ9myZC0cL2vDGlBuA/g3053+sca63VaW8tSZIEVDAJi4jfAIOAbVJKUxax+9vAyk22rQjMB6ZVKiZJkqRq1eF78SK7gNymYtuU0ittOGwM8N0m27YDJrQwH0ySJKmiImLPIj+/Eg0RLgT2A4YBsyJi5dLjswt3iIizI+LqRsdcAvxXRJwbEetFxMHAIcAvKxCPJElSq0pz0Q8oMoZKDEceWXp+qMn2g4Dfl75eBVhj4QsppVciYkfgV8BIcrPWY9raI6wC6nYSoCRJAmBf4NoiA6jlZYuCfIfkrFSrP4QkSepyEbEE8A9g/ZTSR0XFUbNdrEqJV7mtKSRJkgYCjxSZgEFl5oRJkiTVksKHIqGGhyMlSZLKFRHLAM8Aa6eU5hcZi5UwSZJUUyJio4i4JiIeiYidI+LzEXFeRFwcETdHxMatHL47cFvTBCwitoqIqyNiTETsFBGLRcTIiPhNRFxS+qy+pX2HlT7r4tL2bdvzc5iESZKkWnMsuQvD3cAVwGXAL4AbgJ3Iba9a8qmhyIhYDBhBblnxKHAl8GtgckrpeymlI4DngOsi4kxgTkrpyJTSkcADwJ9L71EWkzBJklQzImIt4M1SJWtVYHngrJTSG8AXgA+A21o4diXgv1JK45q81BuYVLrpb1Xgi8AdKaUHG+0zE1gTmJpS+kuj7VNLMXyx3J+lZu+OlCRJdWklYGFf0a3Idzk+AZBSuhG4sZVjhwJ/amZ7D+Dm0tdbA/eklO5uss+GwCvApU22rwfMph3LLloJkyRJNSOl9FhKaVJErAB8g083i2/NvsB1zbznqJTSqxHxVWB14L7Gr5f6im0BPNRMb9LtgNHtmeRvEiZJkmrRAPIKOA+1ZeeIWIfcFeKFVnZbOMG+6Xv2BpZpuj0iNgDWpvXqW4tMwiRJUi0aADQAY9u4f1t6gw0A3gcmNdnev/T8UJPtw0ox3AgQEa3dEPApzgmTJEm1aAAwNqU0p437783/JVMt6Q+MSil93MxnvZxSeq3J9t2Bv6aUppeqYqu3MRagDpOwRmtOSpKk6lDWOtARsTLwNZqfZN/c/n2B11JKb7eyz3rAynx6yHFJ8nyw5j5rBeDhUnuKU8mtM9qs7pIwcgI2o+ggJEnSf/SivPWgVyK3hrihjfu3ZSjyC8CbwF+abF8O+BD4QzPHHAMcSq6gXZJSKusOybpbtqiNlbBlgSnksuKsdn7UeKBPAccW+dkdPW9F/twdPb5ez1tRn130v9GOHl+v561W/43X83nrqt9tZVXCylG6s/ElYIOUUnv/+3WKuquElf4jt5pt5zwNyBdFOZl54/f4uIhji/zsjp63In/ujh5fr+etqM8u+t9oR4+v1/NWq//G6/m8Ffm7rYK2A8ZVWwIG3h3ZmS4s6NiiP7sjivy5O3p8vZ63ov+bFfnZnreuPbYaji/qs2v5Wq0GbRmKLETdDUe2RUT0JM8b61Vw9l5TPG/t43krn+esfTxv7eN5a59qOG8RsTQwGfhqSmleETG0xkpY8xqA/1d6Vtt53trH81Y+z1n7eN7ax/PWPtVw3gYDd1ZjAgZWwiRJUjcVEZcBV6SUxhQdS3NMwiRJkgrgcKQkSVIBTMIkSZIKYBImSZJUAJMwSZKkAtRdEhYR20TEbRHxZkSkiNitDcf0i4iJETEnIl6OiCO6INSqUu55i4j+pf2aPr7WRSFXhYg4PSIej4hZEfFORNwSEeu24bi6vebac8683iAijoyIpyNiZukxJiJ2WMQxdXudLVTuefNaa17p322KiPMWsV/dX3ON1V0SBiwNPAUc3ZadI2JN4E5gNLAJcBZwfkQM6bQIq1NZ562RdYFVGj3+UeG4ql0/csfpvsB3yUuF3VtqINgsr7nyz1kj9Xy9TQFOAzYrPR4Abo2I9Zvb2evsP8o6b43U87X2CRHRGxgBPL2I/bzmmqjrFhURkYDBKaVbWtnn58CglNJ6jbZdAmyUUtq886OsPm08b/2BB4HPp5Te75LAakBEfBF4B+iXUhrVwj5ec4208Zz1x+vtUyLiPeDklNLlzbzmddaCRZy3/nit/UdELAM8ARwFnAE8mVI6roV9veaaqMdKWLk2B+5tsu0eYLOIWLKAeGrNpIh4KyLuj4gBRQdTBXqVnt9rZR+vuU9qyzlbyOsNiIjFI2IouYLdUpNKr7Mm2njeFvJayy4E7kgp/a0N+3rNNbFE0QHUgJWBqU22TSWfuxWAt7o8otrwFrk8PRHoAewP3B8R/VuqZnR3ERHAucAjKaVnW9nVa66kjHPm9QZExAbk5OEzwAfkivXkFnb3Oisp87x5rZWUEtZNgd5tPMRrrgmTsLZpOmYbLWxXSUrpBeCFRpvGRMSXgJOAuvpF1cgFwIbAVm3Y12sua9M583r7jxeAjYHlgCHAVRHRr5WEwussa/N581rLSj/zr4HtUkpzyjjUa64RhyMX7W1y9t7YisB8YFrXh1PTxgJrFx1EESLiN8AgYEBKacoidveao+xz1py6u95SSnNTSi+llCaklE4n30xzbAu7e52VlHnemlN31xrwTfL1MjEi5kfEfPJNNceUvl+8mWO85pqwErZoY4BdmmzbDphQrauyV7FNqLNyc2k47TfAYKB/SumVNhxW19dcO89Zc+ruemtGkIfMmlPX19kitHbemlOP19r9wAZNtl0J/B34eUppQTPHeM01UXdJWOlOjq822rRmRGwMvJdSei0izgZWSykNL71+CXB0RJwLXEaeWHgIsE8Xhl24cs9bRBwHvAo8BywF7Ecu89fbrcgXAsOAXYFZEbHwr8AZKaXZAF5zn1L2OfN6g4g4C7gLeB1YFhgK9AcGll73OmtGuefNay1LKc0CPjFPMyI+BKYtnL/pNbdodZeEkfvAPNjo+3NLz1cBB5L7vayx8MWU0isRsSPwK2Ak8CZwTErppi6JtnqUdd7Iv5x+CawGzCb/wtoppXRnp0daXY4sPT/UZPtBwO9LX3vNfVLZ5wyvN4CVgGvI52YGuWfTwJTSfaXXvc6aV9Z5w2utHF5zi1DXfcIkSZKK4sR8SZKkApiESZIkFcAkTJIkqQAmYZIkSQUwCZMkSSqASZgkSVIBTMIkSZIKYBImSZJUAJMwSZKkApiESZIkFcAkTJIkqQD/H14i8AANjD/+AAAAAElFTkSuQmCC\n", "text/plain": [ "Graphics object consisting of 1 graphics primitive" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(rhof, (1, 4),\n", " axes_labels = [r'$r/m$', r'$P(r)/m$'])" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "rmin = 1.00001\n", "r_0 = 2\n", "\n", "def xe(r, ph):\n", " return rhof(r)*cos(ph)\n", "\n", "def ye(r, ph):\n", " return rhof(r)*sin(ph)\n", "\n", "def ze(r, ph):\n", " return zz(r_0, r)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "Graphics3d Object" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rmax = 10\n", "graph = parametric_plot3d([xe, ye, ze], (rmin, rmax), (0, 2*pi), \n", " color='lightsteelblue')\n", "graph" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Constant $r$ curve:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "def const_r(r0, color='blue', thickness=5):\n", " return parametric_plot3d([lambda ph: xe(r0, ph), lambda ph: ye(r0, ph), \n", " lambda ph: ze(r0, ph)], (0, 2*pi), \n", " color=color, thickness=thickness)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Adding the ergosurface, i.e. the curve $r=m(1 + \\sin\\theta)$:" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "graph += const_r(1 + sqrt(1 - mu0^2), thickness=7)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Adding the curves $r/m = 1 + 10^{-5},\\ldots, 1 + 10^{-1}$:" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "for i in range(1, 6):\n", " graph += const_r(1 + 10^(-i), color='black')" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "Graphics3d Object" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "show(graph, axes_labels_style=dict(fontsize=12, fontfamily='Liberation Sans'))" ] } ], "metadata": { "kernelspec": { "display_name": "SageMath 9.4.beta3", "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.5" } }, "nbformat": 4, "nbformat_minor": 4 }