{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Critical curve of a Kerr black hole onto the EHT image of M87*\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.3.beta8, Release Date: 2021-03-07'" ] }, "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": [], "source": [ "import matplotlib.image as mpimg\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Functions $\\ell_{\\rm c}(r_0)$ and $q_{\\rm c}(r_0)$ for critical null geodesics\n", "\n", "We use $m=1$ and denote $r_0$ simply by $r$." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "a, r = var('a r') " ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left( a, r \\right) \\ {\\mapsto} \\ -\\frac{a^{2} {\\left(r + 1\\right)} + {\\left(r - 3\\right)} r^{2}}{a {\\left(r - 1\\right)}}$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left( a, r \\right) \\ {\\mapsto} \\ -\\frac{a^{2} {\\left(r + 1\\right)} + {\\left(r - 3\\right)} r^{2}}{a {\\left(r - 1\\right)}}\n", "\\end{math}" ], "text/plain": [ "(a, r) |--> -(a^2*(r + 1) + (r - 3)*r^2)/(a*(r - 1))" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lsph(a, r) = (r^2*(3 - r) - a^2*(r + 1))/(a*(r -1))\n", "lsph" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left( a, r \\right) \\ {\\mapsto} \\ -\\frac{{\\left({\\left(r - 3\\right)}^{2} r - 4 \\, a^{2}\\right)} r^{3}}{a^{2} {\\left(r - 1\\right)}^{2}}$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left( a, r \\right) \\ {\\mapsto} \\ -\\frac{{\\left({\\left(r - 3\\right)}^{2} r - 4 \\, a^{2}\\right)} r^{3}}{a^{2} {\\left(r - 1\\right)}^{2}}\n", "\\end{math}" ], "text/plain": [ "(a, r) |--> -((r - 3)^2*r - 4*a^2)*r^3/(a^2*(r - 1)^2)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "qsph(a, r) = r^3 / (a^2*(r - 1)^2) * (4*a^2 - r*(r - 3)^2)\n", "qsph" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The radii $r_+$ and $r_-$ of the two horizons:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "rp(a) = 1 + sqrt(1 - a^2)\n", "rm(a) = 1 - sqrt(1 - a^2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Critical radii $r_{\\rm ph}^{**}$, $r_{\\rm ph}^*$, $r_{\\rm ph}^+$, $r_{\\rm ph}^-$, $r_{\\rm ph}^{\\rm ms}$ and $r_{\\rm ph}^{\\rm pol}$ " ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}a \\ {\\mapsto}\\ \\cos\\left(\\frac{2}{3} \\, \\pi + \\frac{2}{3} \\, \\arcsin\\left(a\\right)\\right) + \\frac{1}{2}$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}a \\ {\\mapsto}\\ \\cos\\left(\\frac{2}{3} \\, \\pi + \\frac{2}{3} \\, \\arcsin\\left(a\\right)\\right) + \\frac{1}{2}\n", "\\end{math}" ], "text/plain": [ "a |--> cos(2/3*pi + 2/3*arcsin(a)) + 1/2" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rph_ss(a) = 1/2 + cos(2/3*asin(a) + 2*pi/3)\n", "rph_ss" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}a \\ {\\mapsto}\\ 4 \\, \\cos\\left(\\frac{4}{3} \\, \\pi + \\frac{1}{3} \\, \\arccos\\left(-a\\right)\\right)^{2}$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}a \\ {\\mapsto}\\ 4 \\, \\cos\\left(\\frac{4}{3} \\, \\pi + \\frac{1}{3} \\, \\arccos\\left(-a\\right)\\right)^{2}\n", "\\end{math}" ], "text/plain": [ "a |--> 4*cos(4/3*pi + 1/3*arccos(-a))^2" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rph_s(a) = 4*cos(acos(-a)/3 + 4*pi/3)^2\n", "rph_s" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}a \\ {\\mapsto}\\ 4 \\, \\cos\\left(\\frac{1}{3} \\, \\arccos\\left(-a\\right)\\right)^{2}$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}a \\ {\\mapsto}\\ 4 \\, \\cos\\left(\\frac{1}{3} \\, \\arccos\\left(-a\\right)\\right)^{2}\n", "\\end{math}" ], "text/plain": [ "a |--> 4*cos(1/3*arccos(-a))^2" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rph_p(a) = 4*cos(acos(-a)/3)^2\n", "rph_p" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}a \\ {\\mapsto}\\ 4 \\, \\cos\\left(\\frac{1}{3} \\, \\arccos\\left(a\\right)\\right)^{2}$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}a \\ {\\mapsto}\\ 4 \\, \\cos\\left(\\frac{1}{3} \\, \\arccos\\left(a\\right)\\right)^{2}\n", "\\end{math}" ], "text/plain": [ "a |--> 4*cos(1/3*arccos(a))^2" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rph_m(a) = 4*cos(acos(a)/3)^2\n", "rph_m" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We add the radius of the marginally stable orbit:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}a \\ {\\mapsto}\\ -{\\left(-a^{2} + 1\\right)}^{\\frac{1}{3}} + 1$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}a \\ {\\mapsto}\\ -{\\left(-a^{2} + 1\\right)}^{\\frac{1}{3}} + 1\n", "\\end{math}" ], "text/plain": [ "a |--> -(-a^2 + 1)^(1/3) + 1" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rph_ms(a) = 1 - (1 - a^2)^(1/3)\n", "rph_ms" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "as well as the radius of outer and inner polar orbits:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}a \\ {\\mapsto}\\ 2 \\, \\sqrt{-\\frac{1}{3} \\, a^{2} + 1} \\cos\\left(\\frac{1}{3} \\, \\arccos\\left(-\\frac{a^{2} - 1}{{\\left(-\\frac{1}{3} \\, a^{2} + 1\\right)}^{\\frac{3}{2}}}\\right)\\right) + 1$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}a \\ {\\mapsto}\\ 2 \\, \\sqrt{-\\frac{1}{3} \\, a^{2} + 1} \\cos\\left(\\frac{1}{3} \\, \\arccos\\left(-\\frac{a^{2} - 1}{{\\left(-\\frac{1}{3} \\, a^{2} + 1\\right)}^{\\frac{3}{2}}}\\right)\\right) + 1\n", "\\end{math}" ], "text/plain": [ "a |--> 2*sqrt(-1/3*a^2 + 1)*cos(1/3*arccos(-(a^2 - 1)/(-1/3*a^2 + 1)^(3/2))) + 1" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rph_pol(a) = 1 + 2*sqrt(1 - a^2/3)*cos(1/3*arccos((1 - a^2)/(1 - a^2/3)^(3/2)))\n", "rph_pol" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}a \\ {\\mapsto}\\ 2 \\, \\sqrt{-\\frac{1}{3} \\, a^{2} + 1} \\cos\\left(\\frac{2}{3} \\, \\pi + \\frac{1}{3} \\, \\arccos\\left(-\\frac{a^{2}}{{\\left(-\\frac{1}{3} \\, a^{2} + 1\\right)}^{\\frac{3}{2}}} + \\frac{1}{{\\left(-\\frac{1}{3} \\, a^{2} + 1\\right)}^{\\frac{3}{2}}}\\right)\\right) + 1$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}a \\ {\\mapsto}\\ 2 \\, \\sqrt{-\\frac{1}{3} \\, a^{2} + 1} \\cos\\left(\\frac{2}{3} \\, \\pi + \\frac{1}{3} \\, \\arccos\\left(-\\frac{a^{2}}{{\\left(-\\frac{1}{3} \\, a^{2} + 1\\right)}^{\\frac{3}{2}}} + \\frac{1}{{\\left(-\\frac{1}{3} \\, a^{2} + 1\\right)}^{\\frac{3}{2}}}\\right)\\right) + 1\n", "\\end{math}" ], "text/plain": [ "a |--> 2*sqrt(-1/3*a^2 + 1)*cos(2/3*pi + 1/3*arccos(-a^2/(-1/3*a^2 + 1)^(3/2) + 1/(-1/3*a^2 + 1)^(3/2))) + 1" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rph_pol_in(a) = 1 + 2*sqrt(1 - a^2/3)*cos(1/3*arccos((1 - a^2)/(1 - a^2/3)^(3/2)) + 2*pi/3)\n", "rph_pol_in" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}r_{\\rm ph}^{**} = -0.477673658836338$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}r_{\\rm ph}^{**} = -0.477673658836338\n", "\\end{math}" ], "text/plain": [ "r_{\\rm ph}^{**} = -0.477673658836338" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}r_{\\rm ph}^{\\rm ms} = 0.539741795874205$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}r_{\\rm ph}^{\\rm ms} = 0.539741795874205\n", "\\end{math}" ], "text/plain": [ "r_{\\rm ph}^{\\rm ms} = 0.539741795874205" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}r_{\\rm ph}^{*} = 0.658372153864346$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}r_{\\rm ph}^{*} = 0.658372153864346\n", "\\end{math}" ], "text/plain": [ "r_{\\rm ph}^{*} = 0.658372153864346" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}r_- = 0.687750100080080$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}r_- = 0.687750100080080\n", "\\end{math}" ], "text/plain": [ "r_- = 0.687750100080080" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}r_+ = 1.31224989991992$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}r_+ = 1.31224989991992\n", "\\end{math}" ], "text/plain": [ "r_+ = 1.31224989991992" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}r_{\\rm ph}^+ = 1.38628052846298$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}r_{\\rm ph}^+ = 1.38628052846298\n", "\\end{math}" ], "text/plain": [ "r_{\\rm ph}^+ = 1.38628052846298" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}r_{\\rm ph}^- = 3.95534731767268$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}r_{\\rm ph}^- = 3.95534731767268\n", "\\end{math}" ], "text/plain": [ "r_{\\rm ph}^- = 3.95534731767268" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}r_{\\rm ph}^{\\rm pol} = 2.49269429554008$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}r_{\\rm ph}^{\\rm pol} = 2.49269429554008\n", "\\end{math}" ], "text/plain": [ "r_{\\rm ph}^{\\rm pol} = 2.49269429554008" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}r_{\\rm ph}^{\\rm pol,in} = -0.399338575773941$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}r_{\\rm ph}^{\\rm pol,in} = -0.399338575773941\n", "\\end{math}" ], "text/plain": [ "r_{\\rm ph}^{\\rm pol,in} = -0.399338575773941" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "a0 = 0.95\n", "# a0 = 1\n", "show(LatexExpr(r'r_{\\rm ph}^{**} = '), n(rph_ss(a0)))\n", "show(LatexExpr(r'r_{\\rm ph}^{\\rm ms} = '), n(rph_ms(a0)))\n", "show(LatexExpr(r'r_{\\rm ph}^{*} = '), n(rph_s(a0)))\n", "show(LatexExpr(r'r_- = '), n(rm(a0)))\n", "show(LatexExpr(r'r_+ = '), n(rp(a0)))\n", "show(LatexExpr(r'r_{\\rm ph}^+ = '), n(rph_p(a0)))\n", "show(LatexExpr(r'r_{\\rm ph}^- = '), n(rph_m(a0)))\n", "show(LatexExpr(r'r_{\\rm ph}^{\\rm pol} = '), n(rph_pol(a0)))\n", "show(LatexExpr(r'r_{\\rm ph}^{\\rm pol,in} = '), n(rph_pol_in(a0)))" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "def alpha(a, th_obs, r0):\n", " if a == 1:\n", " ell = - r0^2 + 2*r0 + 1\n", " else:\n", " ell = lsph(a, r0)\n", " return - ell / sin(th_obs)\n", "\n", "def Theta(a, th_obs, r0):\n", " if a == 1:\n", " ell = - r0^2 + 2*r0 + 1\n", " q = r0^3 * (4 - r0)\n", " else:\n", " ell = lsph(a, r0)\n", " q = qsph(a, r0)\n", " return q + cos(th_obs)^2 * (a^2 - ell^2/sin(th_obs)^2)\n", "\n", "def beta(a, th_obs, r0, eps_theta=1):\n", " return eps_theta * sqrt(Theta(a, th_obs, r0,))" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "def r0_bounds(a, th_obs, outer=True):\n", " r\"\"\"\n", " Return `(r0_min, r0_max)`\n", " \"\"\"\n", " if outer:\n", " r1 = n(rph_p(a))\n", " r2 = n(rph_m(a))\n", " r3 = rph_pol(a)\n", " else:\n", " r1 = n(rph_ss(a))\n", " r2 = 0\n", " r3 = rph_pol_in(a)\n", " #\n", " # Computation of rmin:\n", " try:\n", " if a == 1:\n", " th_crit = n(asin(sqrt(3) - 1))\n", " if n(th_obs) < th_crit or n(th_obs) > n(pi) - th_crit:\n", " rmin = find_root(lambda r: Theta(a, th_obs, r), r1, r3)\n", " else:\n", " rmin = 1\n", " else:\n", " rmin = find_root(lambda r: Theta(a, th_obs, r), r1, r3)\n", " except TypeError: # special case th_obs = pi/2\n", " rmin = r1 \n", " #\n", " # Computation of rmax:\n", " try:\n", " rmax = find_root(lambda r: Theta(a, th_obs, r), r3, r2)\n", " except TypeError: # special case th_obs = pi/2\n", " rmax = r2 \n", " #\n", " return (rmin, rmax)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "def shadow_plot(a, th_obs, orientation=0, outer=True, color=None, number_colors=5, \n", " range_col=None, thickness=2, linestyle='-', plot_points=200,\n", " legend='automatic', legend_loc=(1.02, 0.36), fill=True, \n", " fillcolor='grey', draw_NHEKline=True, \n", " draw_spin=False, spin_arrow_length=7, spin_arrow_options=None,\n", " frame=True, axes=True, axes_labels='automatic', gridlines=True):\n", " if a==0:\n", " # Case a = 0:\n", " rs = 3*sqrt(3)\n", " if color is None:\n", " color = 'black'\n", " if legend == 'automatic':\n", " legend = None\n", " g = parametric_plot((rs*cos(x), rs*sin(x)), (x, 0, 2*pi), \n", " color=color, thickness=thickness,\n", " linestyle=linestyle, legend_label=legend, \n", " fill=fill, fillcolor=fillcolor, frame=frame, \n", " axes=axes, gridlines=gridlines)\n", " else:\n", " # Case a != 0\n", " rmin, rmax = r0_bounds(a, th_obs, outer=outer)\n", " if rmin > 0:\n", " rmin = 1.00000001*rmin\n", " rmax = 0.99999999*rmax\n", " else:\n", " rmin = 0.9999999*rmin\n", " rmax = 1.0000001*rmax\n", " print(\"rmin : \", rmin, \" rmax : \", rmax)\n", " co = cos(orientation)\n", " so = sin(orientation)\n", " fa = lambda r: co*alpha(a, th_obs, r) - so*beta(a, th_obs, r)\n", " fb = lambda r: so*alpha(a, th_obs, r) + co*beta(a, th_obs, r)\n", " fam = lambda r: co*alpha(a, th_obs, r) - so*beta(a, th_obs, r, eps_theta=-1)\n", " fbm = lambda r: so*alpha(a, th_obs, r) + co*beta(a, th_obs, r, eps_theta=-1)\n", " if range_col is None:\n", " range_col = r0_bounds(a, pi/2, outer=outer)\n", " rmin_col, rmax_col = range_col \n", " print(\"rmin_col : \", rmin_col, \" rmax_col : \", rmax_col)\n", " dr = (rmax_col - rmin_col) / number_colors\n", " rm = rmin_col + int((rmin - rmin_col)/dr)*dr\n", " r1s = rmin\n", " r_ranges = []\n", " while rm + dr < rmax:\n", " col = hue((rm - rmin_col)/(rmax_col - rmin_col + 0.1))\n", " r2s = rm + dr\n", " r_ranges.append((r1s, r2s, col))\n", " rm += dr\n", " r1s = r2s\n", " if color is None:\n", " col = hue((rm - rmin_col)/(rmax_col - rmin_col + 0.1))\n", " else:\n", " col = color\n", " r_ranges.append((r1s, rmax, col))\n", " g = Graphics()\n", " legend_label = None # a priori\n", " if a == 1 and draw_NHEKline:\n", " th_crit = asin(sqrt(3) - 1)\n", " if th_obs > th_crit and th_obs < pi - th_crit:\n", " # NHEK line\n", " alpha0 = -2/sin(th_obs)\n", " beta0 = sqrt(3 - cos(th_obs)**2 *(6 + cos(th_obs)**2))/sin(th_obs)\n", " alpha1 = co*alpha0 - so*beta0\n", " beta1 = so*alpha0 + co*beta0\n", " alpha2 = co*alpha0 + so*beta0\n", " beta2 = so*alpha0 - co*beta0\n", " if legend == 'automatic':\n", " legend_label = r\"$r_0 = m$\"\n", " if color is None:\n", " colNHEK = 'maroon'\n", " else:\n", " colNHEK = color\n", " g += line([(alpha1, beta1), (alpha2, beta2)], color=colNHEK, \n", " thickness=thickness, linestyle=linestyle,\n", " legend_label=legend_label)\n", " if fill:\n", " g += polygon2d([(fa(rmax), fb(rmax)), (alpha1, beta1), (alpha2, beta2)], \n", " color=fillcolor, alpha=0.5)\n", " for rg in r_ranges:\n", " r1s, r2s = rg[0], rg[1]\n", " col = rg[2]\n", " if legend:\n", " if legend == 'automatic':\n", " if draw_NHEKline and abs(r1s - 1) < 1e-5:\n", " legend_label = r\"${:.2f}\\, m < r_0 \\leq {:.2f}\\, m$\".format(\n", " float(r1s), float(r2s))\n", " else:\n", " legend_label = r\"${:.2f}\\, m \\leq r_0 \\leq {:.2f}\\, m$\".format(\n", " float(r1s), float(r2s))\n", " else:\n", " legend_label = legend\n", " g += parametric_plot((fa, fb), (r1s, r2s), plot_points=plot_points, color=col, \n", " thickness=thickness, linestyle=linestyle,\n", " legend_label=legend_label, \n", " frame=frame, axes=axes, gridlines=gridlines)\n", " g += parametric_plot((fam, fbm), (r1s, r2s), plot_points=plot_points, color=col, \n", " thickness=thickness, linestyle=linestyle)\n", " if fill:\n", " g += parametric_plot((fa, fb), (rmin, rmax), fill=True, fillcolor=fillcolor, \n", " thickness=0)\n", " g += parametric_plot((fam, fbm), (rmin, rmax), fill=True, fillcolor=fillcolor, \n", " thickness=0)\n", " if draw_spin:\n", " if not spin_arrow_options:\n", " spin_arrow_options = {}\n", " if 'color' not in spin_arrow_options:\n", " spin_arrow_options['color'] = color\n", " g += arrow2d((0,0), (-so*spin_arrow_length, co*spin_arrow_length), \n", " **spin_arrow_options)\n", " # end of case a != 0\n", " g.set_aspect_ratio(1)\n", " if axes_labels:\n", " if axes_labels == 'automatic':\n", " g.axes_labels([r\"$(r_{\\mathscr{O}}/m)\\; \\alpha$\", \n", " r\"$(r_{\\mathscr{O}}/m)\\; \\beta$\"])\n", " else:\n", " g.axes_labels(axes_labels)\n", " if legend:\n", " g.set_legend_options(handlelength=2, loc=legend_loc)\n", " return g" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "rmin : 1.00000001000000 rmax : 3.99999996000000\n", "rmin_col : 1 rmax_col : 4.0\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAGFCAYAAAAB9K+8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB+qklEQVR4nO2dd3yURfrAv7ObbCq99yqKKOqB2Nud/Tw8y9mxd9TT09OfZ7nzDg/UsyBgw4LSQYggPYJIhxAMCYaaACGEhISQhPRkd35/bBKyIWWTfdvuzvfz2U+S2XlnnifzvvvszDzzPEJKiUKhUCgU3mAzWwCFQqFQ+A/KaCgUCoXCa5TRUCgUCoXXKKOhUCgUCq9RRkOhUCgUXqOMhkKhUCi8RhkNhUKhUHiNMhoKhUKh8JoQswWwOkIIAXQHTpgti0KhUOhIKyBDNnHiWxmNpukOpJsthEKhUBhAT+BwYxWU0WiaEwCHDh2idevWZsviEykpKQwYMMBsMXSnoKCAXr16BcSYeUOwjGuw6AnG61r9zODFiooyGl7SunVrv/8A2r17N+edd57ZYhhGIIyZNwTLuAaLnmBtXdVGeBDRtWtXs0VQ6ECwjGuw6AnW1lXNNIKIkBA13KZx/DgcOwbdukFUlLssMxO2bIHWreHKK0/WXbsW8vNh+HCo/vAoLYWjR91127b1aDpYxjVY9ARr66pmGkHE/v37zRYhsKishJQU9wd6NatXw5/+BC++WFMkpUSOGAGnnUbp2rWcOHGC/Px8Cn/6CW6+mYqXXyYnJ6fmVfHCC/CnP1G4apW7XmEhZWvWQJ8+yIsu8pThtdfo+eKLsH79ybKyMsjIAJdLX/0NJpjuXyvral1zphFCiB7AO8ANQASwB3hEShlvqmAmcPnll5stgn9SWAi//grl5fCHP9QUy0GDEPv3UxgbS/HQoZSVlRGSnEy3RYsoHDqU5C1bqKiooLKykvNCQ4mMiGBXcjJ54eEAtM7LY+DgwRR27syeHTtq2j29UyeizjiDlGPHyP/1VwDaJSRwdmgohUDCmjWEhITgcDgYvGgRfRMTOZKUhLNPH8LCwojcvp2oa65BnnYaYs+ek3ps2AAOB5x5JkRGGvKv05Jgun+trGtAGw0hRDtgPfAzbqNxFBgA5JkolmnMmDGDV1991WwxrE1uLmzaBL/7Xc3SkOvHH7Hdcw8V553Hoe+/p6SkhOLiYk5r04bWoaHs37SJY1XLCWGdOlHw4ouUdO9OcXFxTbPbJkwAu92jq4KzzmLbJ5+cIsLuV145pez4iBGsWbEC4XQiXS7Ky8spLy8n5a67+C0ykr5t21K2bx8AHTdtYojNRkFUFPvi44mIiCAyMpIezz5L6LZtuGbPxnbHHe6Gs7MhKcmtb51lL6sRTPevlXUVgZy5TwgxDrhESnmZD220BvLz8/P93hPH6XRir/PBFYgUFBTQpk0bmhwzKd0fmp07nywaMQIRF8exjz/m6NVXU1hYiNy1i6EvvUT+2Wez8/XXa+qGFBZSGRkJNnNXeesbV1FZib2oiMo2bWrKznrtNVrv3Mn2Dz5AnHUW0dHRdIqNpcNzzyEvuACxadPJBjIzoUsXEMIoNZokWO5fMF7X6mcGaCOlLGisbqDvaYwEtgoh5gohjgohfhVCPNbYBUKIMCFE6+oX7lOSDVJcDHl5GkqsI2PHjjVbBOuwaxf07o085xyOZmWxd+9e4uPjOdy7N8W9epGVnU1WVhZFRUUU9+rFptmzPQwGQGV0tOkGA9zfSusiQ0I8DAbAjrffZsP8+RT16UNhYSGZmZkczcigpFs3Mrp3Jz4+nr1793I0Kwt5zjnuTftdu4xSo0mC6f61sq7m3/H60h94CtgLXAd8BnwshLi/kWteBfJrvdIBUlNTeffddykrK2PMmDEAXHfdSqKjJe3aSU47rZjZs1cxf/58kpOT+fjjjykoKKipO2bMGHJzc5k0aRKJiYksXLiQxYsXs23bNj7//HOys7M96hYXF/P++++ze/du5syZw8qVK9mwYQNTpkzh0KFDHnWdTidjx47l4MGDTJ06lXXr1rF69WpmzZrF3r17ee+99yittVk7ZswYMjMzmTx5MnFxcSxbtoyYmBh27NjBhAkTyM/P92g/Ly+PiRMnkpSUxA8//MCSJUuIj4/niy++ICsry6NuSUkJ77//Pnv37mX27NmsWrWK9evX891335GWlsbbb7+Ny+VizJgxuFwu3n77bdLS0vjuu+9Yv349q1atYvbs2ezdu5f333+fkpISj/azsrL44osviI+PZ8mSJfzwww8kJSUxceJE8vLyauoC5OfnM2HCBPbOm8eeW24h8aWX2LhxIx9++CEbjxyh/MgR5LFj/O/ZZ9m/fz+TJ09m9ciRjHvwQVa0b09SUhIrVqwgMzOT6dOn43Q6mTp1KgBTp07l6NGjLFu2jB07dhAfH8/PP/9MWloac+bMobS01KNuXl4eCxcuZM+ePWzcuJH169eTkpJCTEwMBQUFHnWLioqYN28eBw4cYM2aNWzZsoVdu3axePFicnJyPOqWl5cTEhLC4cOHWblyJQkJCSQmJhIbG0tmZiYzZsw4Re7snByWLl1KcnIyizt04J3HH2fVyJF8/fXX7Nu3j3efew55/DgVOTlsPnqU9957j7i4OH7929/Yf9NN7Jk6lU8//ZRjx455jM2JEyf46KOPSE5OZt68eaxYsYLNmzfz9ddfk5GR4VG3oqKCd955h9TUVKZPn86aNWtYs2YN06dPJzU1lXfeeYeKigqPax5++GG+/vprNm/ezIoVK5g3bx7Jycl89NFHnDhxwqPusWPH+PTTT0lISODHH39k0aJFJCQk8Nlnn5GTk+NRt7CwkA8//JBdu3Yxd+5cYmNj2bRpE9988w3p6ekedSsrKxk3bhwHDhxg2rRprF27ltWrVzNz5kxSUlJO+YwYM2YMGRkZfPXVV8TFxbF8+XKvPiNGjhxp6GfE+PHjG/lI9CTQl6fKga1SyotrlX0MnC+lvKiBa8KAsFpFrYD0+pY6IiOhpOTk323awOTJ8Je/aKeDlowZM4bX63xbDkQKsrJo07Ur+fn5tGrVivJx4wj7xz8ouOgifh07lup7vtXu3RT16YOramPaX5k6dSqjRo3StE1beTkRaWkUDRxYU3be//0fbTZvpuTttwl/9VWEEG4Psrw86NhR0/7rI1juXzBeV7U8dZIjQHKdsp1A74YukFKWSSkLql80cqy+stLz7/x8uOMOuP12yMpqudB68cADD5gtgv48+yyyf3/APTvcsmULv/buTdbVV5N2443U/pJ04vTT/d5gAFx33XWat+lyODwMBsDBm28m409/IrFXLzZs2MDOnTvJW7AA2aUL3Hmn5jLUJSju3yqsrGugG431wOl1ygYBB7VovEOH+svnzYMhQ2DGDPdeq1VYuXKl2SJoi5Q1a+5SSvf5h6wsRJXX0uHDhykpKaG0Rw92vvYaOZdeaqa0urFt2zZD+sm96CL2/O1vlPTqRUVFBVlZWeQuXYpwuTjudHL06FGcTqe78rJlbldlDQm4+7cRrKxrQLvcAh8CG4QQ/wDmACOAx6tePtOrl9vJpD6OHYN774XZs+Gzz9x7imYzaNAgs0XQjsJCuOAC5M6dHN6wgQybjeLiYiJvuonSESPg7383W0LD6Nmzp2l9p913H0f/8AeQktLkZOx2O13Lyxl4440QFYVIS4N27TTpK6Du3yawsq4BPdOQUsYBtwB3AzuAN4DnpZTTtWjfG7f2hQvdZ6m+/db8WUdJ7Q0Yf0NKOOyO2OxyuThaXExhRAQuh4PjK1fWnIko7t2bE4MHmymp4ZSVlZnaf2m3bpR27w64XUVPJCdT0qMHBaedRkpuLkVFRe6Kq1a5Q6G0EL++f5uJlXUN9JkGUspFwCI92q7j0dggeXnw4IMwZw58/jmY9cXw2LFj5nTsKwcOwC23IA8dYv/atRzJzaWiooKIF16gvG1bnNWxnIKUgoJG9y0Np+Dss9ny3XeEnDhB5aFDHDp0iLZhYQy97TZEURFiyxY499xmt+u3928LsLKuAT3T0BtvjUY1S5a49zq+/NKcWcc555xjfKctpeofJKUkLyqK8sxMXCdOkLdqFRUVFQCU9OgR9AYDsGaOCSGorOVtWJqaSmHXrpR17Mi+iIiT36RTU92xsrzAr+5fH7Gyrspo+EBLoi4UFMBjj8F119WsthjGokW6TLi05dgxeOUV5PXXczQri23btpHw228k/fOfbJw7l4IhQ8yW0HJsqn2S26KU9ujBtk8/Zeunn5J+5AibN28mKTGRyptvdnu7bdzYZBt+cf9qhJV1DfjlKT1p7kyjNrGxcM01sGOHcYeKn3rqKWM68gFXRQVi/HhEWRnp33/PiSojceKMM0yWzLqMHDnSbBG8pvYp9cJdu3BmZiKKijjWrh2dpHSf/ZCy3vAl/nD/aoWVdVUzDR/wxWgA7NwJBzVx/vWODz74wLjOvCU7G+bPx+VycfjwYTbv38++Rx8l6b//pSDINrRbyty5c80WoUWUde7MphkzSHj/fZIzM9m8eTOHDx9GPvYYPPccHDniUd+S969OWFnXgD4RrgWNBSz87jvw5QzOaae5jxlYIHyROaSlIc88E8rL2TZrFifat9ek2aKiIm666SYWLVpElNrz8CvCMzO54N57ES4Xzs2bsY8YYbZIQYE6EW4Qvsw0Ro6En3821mDUjstkNlJKjoaHUzB4MIX9+iFzc80WyW+pjisVCJR26cL2997jwKhRbK6s5PDhw7hcLli/nv+99prZ4hmGlZ7VuqiZRhM0NtNYvRquuqp57Z15Jnz0kXs/w2iys7Pp1KmT8R1Xs2IFTJhA/pdfsi89nRMnTrjdMqOjNQ3BHWwzjby8PNpaPBeGL0Q5nQy74w5wubCtXg0W9izSCqOfVTXT0AAhxGghRDKwpaE6zXlO7Xb4+GPYvt0cgwEwf/58czoGKC1FPvggLFpE9ltvceKEO6RXZatWlsrZ4I+sWbPGbBF0RaanU9ymDTk2G7lduxIMX3RNfVabQHlPNYCUchIwqXqmUV+d5ixPOZ3QqhWYmS/+/PPPN7ZDlwtsNlwuF2mZmZQ+8ghRu3aRef31xsoR4JwR4J5lxX36sHXyZDK2bKH7zp20y8xk4IABRL35Jtx3nzvrYIBh+LPaDNRMwweammnUTcP8f//nPqdhFkfqeKPoytq1cM45FCxaRFxcHAcOHCDzD38gZfRoj0NfCt+x8ulhzbDba6KMHj9+nIPvvgsffoi8/HL/yYLWDAx9VpuJMho+0Lo1OBynlvfvDz/8AN9/71melQX//a8hotVLTQRSI/r69lvYsQPX669bOo5OIOByucwWwRBq65l39tlk/f73pN11FxnFxQG3ZGXks9pclNHwAbvdnT+jmuhoGDcOkpPh5pvhhhvgj3/0vObDD2HfPmPlrKZv3776diAlUkqysrLYetttHLr9dnb8+9/69qmga9euZotgCLX1LO/UiZ1vvMH+e+9lz549/PrrrxQlJro3DPfsMVFKbdD9WfUBZTR85Kuv3LGkxo9336uvvAJhtfL+ffCB5z5GeTm89JLxcgKsXbtWn4aLi+H553E+/TS//fYbO3fupCQiQi1FGURiYqLZIhhCvXpWOVEUFBRQ/Nxz8NNPyGefNVgy7dHtWdUAtRHuIw4HPPJIw+8PGuQ+3Fr7gOeCBfDTT3D11frLV5u7775bn4bj4mD8eOxA8YUXQp8++vSjqJc//OEPZotgCE3pue/JJxGVlRx69FH6Hj9OO43yeJiBbs+qBqiZhgG88QbUdbl+/vlT08XqzSeffKJ5m06nk91du3LggQdIHDeOYmUwDGfBggVmi2AITelZ1rUrO8aMIb9TJ7Zv386ePXtwfvgh/Oc/xj9sPqLHs6oV6nBfEzR2uK85TJ4Mj9fJFzhxIowe7Zt8ppCXB//8J4Uvv8xv6emW2+gOtsN9ivpx5ORwwX33YS8rg5kz4a67zBbJsqjDfRbk4YdPzTvzxhvuSOBGoVVoAnnrrfDxxxQ99JDlDEYwEkhhRBqjuXqWd+jAnhdf5OhVV5EyfLhfeZlZOYyIMhoGYbe7N8trc/w4/Otfxsnw3HPP+dyG0+nkwOOPU9yrF+m33KKBVApfufXWW80WwRCaracQZF1zDclvvsmh9HS2bdtGcUEB/PvfkF/veV3LoMWzqhfKaBjI5ZfDX/7iWfbpp/Dbb8b0P2XKlJZdWFEBu3dTUlLCtm3bONi1K3HffBN0ubityrJly8wWwRB81bOwsJCjzz8P//wn8rLL3GEaLEqLn1UDUEbDYN59F8LDT/7tdMILLxiT/vXqlrhrHTsGV1+N67LLSFqyhKKiIgCk3a6xdIqWMmzYMLNFMAQt9MwdNoySbt3IuP9+rGsyWvisGoQyGgbTty/8/e+eZbGx8OOP+ve9a9euZl8jw8Ioz8nBVVSEIy1NB6kUvpIWJOOihZ4FQ4YQ98037B0+nPj4ePeXoPR091kjC9GSZ9UolNEwgVdegR49PMv+9jcoK9O33+Z6ErlcLvYcPsy2N95g2yefkHfeeTpJpvCF8NpT1wBGKz1dVadvi4uLSdiwgcprr4WLLoL9+zVpXwus7PWnjEYDeBMavaVERbmXqWqTknLqRrnWtPc2M97s2TjnzGHHjh0cOXKE0q5d1fkLC+OLK7g/oYeejkOHcGVmUpmRYaklV6+fVRNQRqMBpJSTpJRnArrkm7z7bveXm9qMGQOZmXr05sarcBOxsXDXXYhRoyiLj9dPGIVmpKSkmC2CIeihZ1H//mz9/HMS3n6bpLw8Ki1yCNDKoWGU0TAJIU6dWZw4AXpmtLzhhhuarFNy0UXkXnYZh0eOpKh3b/2EUWjGBRdcYLYIhqCXnuWdOlE4aBC5ubls27aN0jVr4OWXTfWu8uZZNQtlNEzk/PPhwQc9y775BrZu1ae/r7/+uv43qly3CgsL+TUxkaQ33yRl9Gj34RKF5Vm6dKnZIhiCEXqWHjuGuO02eO8999TfJBp8Vi2ACiPSBFqFEWmII0fcQQ0LC0+WXXKJO4eRIVlQnU547DFKe/Vi6zXXWGZ67gsqjIjCFzqtXk3PefMoi4mh88CBZotjCCqMiB/RrRu8/rpn2fr1MHu29n3VG5pg4UL45hvC3n6bUAt5jyi8R4UR0ZbsK6/k148/Jjk9nbS0NHeCp6rzSUZh5TAiaqbRBHrPNMDtanvmmZCaerKsZ0/YvfvUlLG+9VNGWO1kH0B+fj7H//53Cnv0IOeKK7TrzESCbaZRXl6Oo74UkgGGWXqenpRE13HjEIsWgUFu5/U9q3qiZhp+RlgYvP++Z1l6+qluub4yvnrnXUpwOsnPzycxMZED99wTMAYjGJk3b57ZIhiCKXq6XER98gkiIwNp4D7DeL39731AGQ2LcPPNUDfHzDvvgJaHfW+77Ta3wXj+eSruuYfEbdssnYtY4R2XX3652SIYgil62mxsf+89Uh97jB2PPGLY83LbbbcZ0k9LUEbDIggBH30EtlojUlrq9vzTii1btsD27chJkwidM4dWv/6qXeMK07ByyAktMUtPZ3Q0affcw7G8PJKSknBWVoLOsmzZovmZYs1QRsNCnHUWPPWUZ9ns2W5PKi3o1q0bpWecwa5x49jz/PPk/e532jSsMJUOHTqYLYIhWEHPvLw8jj79NHLoUPjhB9366datm25t+4oyGhbjrbegbmrjv/5Vg3NGUiKlJDExkazhw8m4+WYfG1RYBZstOB5jS+jpcmHfswdRUYFTR29Du4XPSFlgFBS16dDBbThq8+uv4FN4/c2bkddeS3xsLMUWi+ap8J2srCyzRTAES+hps7HzjTdI+u9/2X7FFbqdazp48KAu7WqBMhoW5Mkn3S64tfnHP1qYbMzpRI4ahfjpJ0Zu26aJfAprcfbZZ5stgiFYRU9pt3PsoosoKCggMTERZ2kpJCRo2sell16qaXtaooyGBQkNdW+K1+bo0RZGNbDbSR8/npyLL2ZC165aiKewGKtWrTJbBEOwop4ncnIovOEG5CWXwLp1mrU7a9YszdrSGmU0GkDP0OjecM01MHKkZ9n48bB3b/PaycrKIiUigh1vv82fR43STkCFZbjzzjvNFsEQrKqns7AQWVmJS8O84y+99JJmbWmNMhoNoHdodG/43//cs45qKirgxRe9vPizzyjcto3du3fXFM2YMUNbARWWIFjG1Yp6SoeDHf/5D79++CG7+vVDqwgb48aN06QdPVBhRJrAiDAijfHyy+6Am7VZtgyuu66Ri5YtgxtuoLJVK7Z88w3lFnBVNJJgCyOisA49evRgYMeOCKcTOnY0WxyvUWFEAojXX4fOnT3LXnjBPetoCHnWWRSeey5ZV13lYTCCJbBdsBEs4+oPemYnJFBx8cXwpz/5lHfcygELldGwOK1bw9ixnmU7d8KnnzZ8zYGKCuLfe499o0d7lF9//fU6SKgwm2AZV3/QM6S4GHH4MM7UVDh0qMXtPFg30Y6FUEbDD3jwQRg2zLPsn/+EnJw6FbOyyM3N5eDBg8iQEGSdiKBb9crupDCVYBlXf9CzuE8fEseNI278ePK6dGlxO7GxsRpKpS3KaPgBNtupLrh5eW7DUcNPPyH79iXfo9CTPn366CGewmSCZVz9Rc8TZ55Jadeu/Pbbb5SUlEB5ebPbGDx4sA6SaYMyGn7CpZfCXXd5ln32GSQluX+Xc+YgSktxNDIlVqfBA5NgGVd/07OiooLUGTOQp53W7DMcJ06c0Ekq31FGw4945x2IiDj5t8vljkslJWS8+SY7//EPUp54osHrC2vnlFUEDMEyrv6oZ4fZsxFpachGVgDqIy8vTx+BNEAZDT+id2945RXPsp9/htmzy0hJTSXrmmtw1bYqdejfv7/OEirMIFjG1R/13PPCC6TdeScHPvywWddZJWRKfQSV0RBCvCqEkEKIj8yWpaX8/e/Qq5f79y5k8hpj+L8XKygtbfrazZs36yucwhSCZVz9UU9XRASpTz7Jwdxcjh496vV1S5cu1VEq3wgaoyGEOB94HEg0WxZfiIw8mQZ2PH9lDG/wRsZf+f77nk1ee9NNN+ksncIMgmVc/V3P3bt3Uzp9Onz7bZN1H3vsMQMkahlBYTSEENHAdOAx4LjJ4vjMnXfCJZfAIm7iCF2ZyDNMndqHnBxHo9d9//33BkmoMJJgGVd/17NVXBzh992HfOIJ+O23Rut+VNdd0kIERRgRIcS3QK6U8gUhxGogQUr5vJfXmhpGpCG2bpWMGAGhspxywgC47rpM/u//giP1Z2OoMCIKS+JycdYbbyBOP532n32GcDT+Jc9IVBiRWggh7gJ+B7zqZf0wIUTr6hfQSlcBW4KU9O6dzQ03ZNYYDIDly7uyc2fD4vpDGAZF8wmWcfV7PW02dvz73yTdfz+Zx441WlWFETEJIUQvYDxwn5TSi61iwG1c8mu90gFSU1N59913KSsrqxnQMWPGkJGRwVdffUVcXBzLly9n/vz5JCcn8/HHH1NQUOBRNzc3l0mTJpGYmMjChQtZvHgx27Zt4/PPPyc7O9ujbnFxMe+//z67d+9mzpw5rFy5kg0bNhDz9tuUDhvGlEcf5ZFHUgkNLfEQfsyYTiQmJpGQkMCqVatIT09n1qxZlJeXU1q1Wz516lRyc3NZtGgRu3btYsuWLaxdu5b9+/czf/58CgsLax7QqVOnUlhYSExMDKmpqaxbt45Nmzaxe/dufvzxR3Jzcz3qlpWVMWfOHNLT01m1ahXbtm0jKSmJ5cuXk5WVxbRp03C5XEydOhWXy8W0adPIyspi+fLlJCUlsW3bthq558yZQ1lZmUf7ubm5/Pjjj+zevZtNmzaxbt06UlNTiYmJ8ZAb3DOO+fPns3//ftauXcuWLVvYtWsXixYtOkXu8vJyZs2aVSN3QkICSUlJrFixgszMTKZPn47T6fS45ujRoyxbtowdO3YQHx/Pzz//TFpaGnPmzKG0tNSjbl5eHgsXLmTPnj1s3LiR9evXk5KSQkxMDAUFBR51i4qKmDdvHgcOHGDNmjU1ci9evJicnJxT5K6srOTw4cOsXLmShIQEEhMTiY2NJTMzkxkzZpwid3Z2NkuXLiU5OZmtW7eyevVq0tLSmDt3LiUlJR518/PzWbhwIfv27WPjxo1s2LCBffv2sWDBAvLz8z3qFhcX8/3333PgwAF++eUX4uLiSE5OZsmSJafIXVlZycyZM8nIyCA2Npbt27ezfft2YmNjycjIYObMmVRWVnpcc+mll7JkyRKSk5OJi4vjl19+4cCBA3z//fcUFxefIveCBQvYt28fGzZsYOPGjezbt4+FCxeeIndJSQlz584lLS2N1atXs3XrVpKTk1m6dCnZ2dkedZ1OJzNmzCAzM5PY2FgSExNJSEhg5cqVHD58mJkzZ1JeXu5xTU5ODosXL3Y/a/HxrFmzhtjYWN57911OJCTU+xlx2WWXtfgzYsqUKRw6dMijrtPpZOzYsRw8eJCpU6eybt06Vq9ezaxZs9i7dy/jx49v+pOxioBenhJC/BmIAWpn2LYDEnABYVJKZ51rwqDW13f3TCPdMstTt98O8+aRc8kl7BgzhlmzevH55wM8qvzjH8lcc82pnhoLFy5kZN0kHQFIsC1PBcu4BpKe9pIShvzvf7TbsgWxfTvUcSf+7LPPePLJJw2TRy1PnWQlcDZwbq3XVtyb4ufWNRgAUsoyKWVB9Quw1NHM4v/9j4yRI0l56ikAbr01nR49PE/KfvHFAEpKTk1Mf2bdHLKKgCBYxjWQ9HQ6HNhycqC0FDZuPOX9Cy+80ASpvCOgjYaU8oSUckftF1AEHKv63e9ILSxkzwsvUNKjBwAOh+Tpp1M86uTkhDFzZq9Trs3OzjZERoWxBMu4BpSedju7Xn2V+EmTOH7jjae8nZ6eboJQ3hHQRiOgKCkhLy+PnFNC28JFFx1j+PBcj7JZs3qTmRnuURbIS5HBTLCMa6DpWdq1K4WDBrF7924qKys93rOyrkFnNKSUV3rrbmsZiouRgwZR+fDDhNQTf0cIGD16HzbbyRutosLGZ595rpN2rpvNSREQBMu4BqqepaWlHFy/Hp5/viYibq9ep64UWIWgMxp+yY8/ItLTidq0CWdYWL1V+vYt5uabD3uU/fJLZxIS2tT8/VsTB4oU/kmwjGvA6ul00v2ee2D8eHdUUmBjPfscVkEZDT9A3nEHOydPZs/zzyNDQxus9+CDB2jd2jMP7MSJp+Gs2u6/8sordZRSYRbBMq4Bq6fdzv6HH+bEkCE4b74ZgDvuuMNkoRpGGQ0/IDs7m6yBAzk+YkSj9Vq3ruShh/Z7lKWkRLNkSTfA7bKoCDyCZVwDWc+jv/898ePHsz8yEoBPG8vnbDIBfU5DC0wNI1JaihSCrYmJFBUVeXWJ0yl47LFh7N8fXVPWpk0506ZtITq6spErA4dgO6ehCCyGDRtGq4gICAkxrE91TiNQ+OwzXH37Ej1vnteX2O2S0aP3eZTl5zv47rs+/h+GQVEvwTKuQaGn08nxf/+b4507Qz2eklZAGQ0LI+fMwZ6Zia2iounKtRg2LI9LL/X0aZ8/vwcXXHC/luIpLMJtt91mtgiGEAx6CqB9TAztjh+HTz4xW5x6UUbDwuTHxLDz//6PrGuvbfa1Tz6ZQmioq+Zvp9PGf//bQUvxFBZhyZIlZotgCMGgp7Tb2fPCC8y+8koqXn7ZbHHqRRkNC5OenU3WddfhasDNtjF69CjlL3855FF28OBZbNrUXivxFBZh+PDhZotgCMGiZ8GQIZQ88ACpBw+aLUq9KKPRAEKI0UKIZGCL4Z2Xl1NSUlLv6e/mcO+9abRvX+ZR9sknA6msFD61q7AWBy364aI1waKnlJLDhw9TWFjIiROWCn0HKKPRIFLKSVLKM4HG/Vz14L774PLLaZWc7FMzkZFOHnvM0wX30KFIfvihh0/tKqxFZJWbZqAT6HqGhoYSFhZGSEgIbdq0oaioiJSUFJASmrmvqSfKaFiNggLkokVEbN2KS4PMXtdem8npp3t60E2Z0pe8vIYPCSr8i2BxKw5EPUNCQggLC8Nut1NRUUFZWRmRkZHYbDZcLheVcXFUXHQR/Oc/ZotagzIaVqN1a7LXr2f3iy9SNHCgz83ZbPDss54uuEVFIXz9dT+f21ZYgwMHDpgtgiEEip42m61mRlFtKJxOJ0IIoqOjOXHiRI2u4UeOELp5M3LSJCgpabxhg1BGw4IcBo7cdJNm7Q0ZUsDVV2d5lC1a1I19+6IbuELhT4xoIlJAoODPelYbitDQUJxOJ2VlZVRWViKEe38xJCSEiIgICqsCklbrmnPZZaQ+8ghZy5ZBRIRp8tdGGQ2LUVxcTH5+vubtPv54KiEhJzfFpRRMnDgQFRDA/1m6dKnZIhiCv+kphKgxFC6Xi7KyMioqKmoMRTUOh4OQkBCKi08mU6vRVQjS7ruP1IoKnM5TcsaZgjIaVuLvf8f54INE7t/fdN1m0qlTGaNGZXiUbd/eljVrOmrel8JYRo0aZbYIhuAveoaFheFwOJBS1hiKhoiIiEBKSWlpqUd5XV3Ly8s5fPgw5OXpIXKzUEbDKpSVIb/6ilZz5+I4flyXLpzOd+jSxfPm/OyzAZSVqdvAnwmK8BpYW0+Hw4HD4UAIQVlZGeVVeTEaIyoqivLy8nqNyim6Op2EvPgiskcP2L1bK7FbhPq0sAqhoRRNncqh228n75xzdOli1KjbefJJz9SwmZkRzJ3bU5f+FMZw9913my2CIVhNT4fDQVhYGEIIysvLKS8v9zrjXqtWrSguLm5wyekUXe12HBkZiOJimDHDV9F9QhkNq2CzkXn66aSMHg12uy5dzJ07lyuuyGbo0DyP8unT+5Cd7bt7r8Ic5s6da7YIhmAFPUNDQwkPD8dms1FeXk5ZWVmzUrMKIWjVqhUnTpxo9Lr6dE194gkSP/6Yytdfb5HsWqGMhkWQUpKdnd10RR+44oorEMLtgivEyRu2tNTO5Mn9G7lSYWWuuOIKs0UwBLP0rHuWorS0FJfL1fSFdbDZbERGRnp1yrs+XYt79yb37LPJOHKk2X1riTIaVmDTJsrefhuRmqprN9XpMgcOLOTGGz1vvNjYriQnG5wvRKEJAZsGtQ5G6mm32wkPDyckJITKysqasxQtJTQ0FIfD4XVenMZ0PXToEM7SUjApxIgyGlZgyhTC33iDXnPm6NpN586da35/5JH9REV5JmWaMGEgLfgCpTCZ2uMayOitZ+1Dd5WVlZSWllJZ6Xvisup9j7oeUo3RmK5tfv4ZOXAg/OtfPsvWEpTRsAKXXEL++eeTc9llhnXZrl0F999/wKNs167WxMZ2MUwGhcJsqg2Fw+Go99Cdr0RGRuJ0Or3ypvIWl8NByOHDyHnzTIlJpYyGBSi74w5+ffddjusc+vno0aMef99yy2F69Sr2KJs8uT/FxfpsxCv0oe64Bipa6SmEqHGRrT50V15erpmhqCY6OpqSkpIWzVYa0zX3/PP57Z//JPuXXyDU+Bhyymg0gJGh0XNzc/XuAoAhQ4Z4/B0aKnn6ac+4VMeOhTF9em9D5FFoQ91xDVR80bPaUISFhSGlrHGR1YtWrVpRWFjYLM+q2jSqqxBkX3kl6ceOtVA631BGowEMC42ekEBeZqauXVTzyy+/nFJ24YW5jBjhefPNnduLjIxwQ2RS+E594xqItETPuoairKys6Yt8oHbQQV/wRteCggIKCgqgGXslWqCMhplUVCAvu4xBF15IxOHDunf3l7/8pd7yp59OwW4/uQNeUWHjs88G6C6PQhsaGtdAwxs9pZQ1eSlqn6UwgrpBB33BG12j9u/Hfu21cOONPvfXHJTRMJOUFIiIwBUeTkm3brp3N3PmzHrL+/Qp5pZbPI3W2rWd2Latre4yKXynoXENNBrTs3YCo+pw4y05S9FS6gs66AvejGllVBSRcXHINWvg0KEm62uFaOmaW7AghGgN5Ofn59O6tfbnGDIOH2b/1q1UtGmjedvN4cSJEEaNGkF+/smT4f37F/LFF/HY7f51jxQVFXHTTTexaNGigEzco3ATEhKC3W6nsrLS1Aiw4eHhOJ3ORgMT6kWX2Fha/fGP9Lz4Yp/aKSgooI37M6iNlLKgsbpqpmEy+QUFhhmMxgK+tWpVycMPe0bXTU2NZtEi/WdACt+wciA/LZk6dSp2u73eBEZmERkZ2WDQQV/wdkyzrrmGdJutxRvuLUEZDZMpKGjUqGvKDTfc0Oj7f/zjEfr391yP/frrfpw4EaKnWAofaWpc/Z1qQzFy5Mia09lanqVoKa1ataKkpESXZbDmjGlpaSnHdYqMXR/KaJjFkSPI4cPp8e67GJUJacuWxr2H7XZ45hlPF9yCglC+/bavjlIpfKWpcfVHaicwqjYU69evN91QVFPtIaXXN/zmjGnY0aPw5JMwcqQustRFGQ2z2LIFER9P219/BYMehL59+zZZ57zz8rj8cs/AiTExPThwIFInqRS+4s24+gPVhqJuAqNqQ2EFPauDDmrhIdUYzdJVCNrNmwc//gj79jVd30eU0TCLESPImTiRAw8+aFiX3gZLe/LJFEJDT065XS7BpEkqNaxV8XZcrUpdQ9HQoTuz9awOOqiVh1RjNEfXsk6dSH38cbJnz4b++kerVkbDLLp1I+cPfyDn8ssN69Lbm71bt1LuvNPThW/r1vZs2tRBD7EUPmLEh5jW1E5g5G2mOzP1bEnQQV9orq6H7ryTQwMHgk3/j3RlNEzE6G9Offr08bruPfek0aGD56GoTz4ZQEWFNdaUFSdpzriahZTylEx3zU1gZJaekZGRVFZW6hp2pC4t0bWgoICSkhIdpPFEGQ0zqKxEzpoFO3ZgZCzyrVu3el03IsLJ44975vdIT48kJqaH1mIpfKQ542o0tRMYtcRQ1MYMPauDDhrt1tsSXR3HjlHy9tswYYIOEp1EGQ0zSElB3H035z75pKHd3tjMcANXX53F4MGeLsHffdeX3FzjI2sqGqa546o3ISEhpyQw0sIt1Wg9o6OjfQo66Ast0bXVrl20f/ttGDcOdDRyymiYQXExlcOGkT9kiCFrkNXMmzevWfVtNnjmmb0eZUVFIXz9dT8txVL4SHPHVQ/qHrrTKoFRbYzSUwhBVFSU7h5SjdESXXPPP59jF1xA+csvg8b/+9q0OIyIECIKcEgpjTtVYiBCiNHAaNyG9XStw4gcOXKE3bt3a9aenowdewYrVnSt+VsIyeefx3PaaeY9VI2hwogYg81mq0lepEc+CjMICQkhNDTUkL0BvejXr1+z90R0DSMihAgRQkwB8oAcIUScEGJ4nTp9hRBnCSGim9u+VdA7NLpRkTdr09JwE489lkp4+MnprpSCiROVC65VMDKMSO1Mdy6Xi9LSUo+zFHqit56hoaHY7XZLGAxfdM3Ozm66kg+0ZG3kRSAKeBYYB4QBG4QQ1wshbhZCpAIpwHbgmBBiiRDiEs0kDhDMMBojW3hitGPHcu6996BHWWJiW1av7qSFWAofaem4ekvtQ3e1M90ZjZ56RkREAOY8l/Xhi67FubmU//ADHDzYZN2W0BKjMVhK+Rcp5WdSyteklEOBPwNfADOBzbiNyWTgV+AaYI0Q4h2NZPZ/HniAnnfdRZukJEO7Xb16dYuvveOOdLp29fwG9vnnAygtVdtiZuPLuDZGtYssYJqhqI1eekZFRdWcPrcKvug6eMwYHLfcAtOnaydQLVryxJ/yn5VSLgF+BmZIKe+uMiZPSikvBLoCrwGPCyHG+iZugBAfT1RCAjaDH0Jf0mU6HC6eeirFoywrK5zZs3v5KpbCR7RM9+rrWQo90SOtbXR0NMXFxYbm3vAGX3TNHT6c8i5dIFyf7JstMRrJQoj6NCoDTvHzklIek1KOA84GbhJCnNeCPgOLL75g95gxFA4wNjteY8nqveGyy3I47zxPv4eZM3tz9GiYT+0qfMPXca0v051VDEVtfNWzLma61DaFL7pm3ngjm2bPxvnXv2oo0UlaYjQmAy8JIepe6wIa/O9LKdOBO4FHW9BnYHHxxWRddhkVbdsa2q2vm5VCwOjR+7DZTg5zWZmdL77QP96NomFaMq61D92ZkemuJWi12W5U0EFf8EVXGRKCS0ry8vK0E6gWzTYaUspCYCIwUQjhqPXW28BHTVybDPitR5VWuFwuUx7QTp1837geMKCIm27K8ChbubILSUnaZzVUeIe341rbUFQfujMzgVFz0eL+DQkJMSzooC9ooWtubi7oYBhbtIsppYwHPgGmCiEurSo7JKXc5cXlwZ3Rp7gY56JFtElIMLzr5ORkTdp56KEDREd7bm1NnHiakRFRFLVobFyrD92FhoZaItOdL/h6/xoddNAXfNU14tAhelx3HZx5pub5elrs+iKl3AGMAoYLIaYJIf4mhDhXNDKvEkI8ACxuaZ8BwaFDhN58M2e//rrhXV9xxRWatNO2bQUPPHDAo2zPnlYsX961/gsUulJ3XKvPUtROYGTUWQo98eX+jYiIoLKy0lIeUo3h67Na1qkT4enpyPR0OHSo6QuagU/+klLKcinlR1LK+4CfgMuBr4QQ3wghXhFCXCKE6CqEGCyEeA0YKqWcoYXg3iCEeLXq8OEJIcRRIcQPQojTjeq/Xux2XOeea/gmOMDChQs1a+vPf86gd2/PKL2TJ/enqMiuWR8K71i4cGG9ZykCwVDUpqX3b1RUFKWlpX41w/L1WXWFh7P93XfJTk6G3r01kspNi8OINNmwEINxG5HLgcuAncC7wAYppSFHLoUQy4BZQBzuZbG3cXtxnSml9CouuRCiNZCvZRiRwsJCS0cm9ZYtW9rzyitDPcruvjvtlOi4RhMsYUSEEDXJi8w+Q2FFrBBDymy6devG6ac3/T1Z1zAi3iKl3Cml/FxKea+UsjfwJNAT+FAI8a0Q4r9CiOv0DDUipbxeSjlFSvmblHI78BDQGximV59eymVKv1qHYRgxIpcLLzzmUfb99z05fDhC034UntQ9dPfVV1+ZLJExNOf+tdlsRERE+K3B0OpZzc/P16Sd2hh2nFdKuV9K+W3Vob8HgAlAO+BBo2QA2lT9zG2oghAiTAjRuvoFtDJGNP25/fbbNW/z6af3Ybef3AGvqLDx6afGL70FOtWGor6zFHqMqxXxVs/qzX+re0g1hlZj2mbOHFyjRkF6uibtgYmh0aWUR6SUs6SUE43or2qD/gNgXdUmfkO8CuTXeqUDpKam8u6771JWVsaYMWMAGDNmDBkZGXz11VfExcWxfPly5s+fT3JyMh9//DEFBQUedXNzc5n25pu4rrySrqNHs3HjRvbs2cPChQvJy8ur+XYxdepUSktLmTNnDmlpafz888/Ex8ezY8cOli1bxtGjRz3qOp1Opk+fTmZmJitWrCApKYmEhARWrVpFeno6s2bNory8nH/961811+Tm5rJo0SJ27drFli1bWLt2Lfv372f+/PkUFhZ6tF9YWEhMTAypqamsW7eOTZs2sXv3bn788Ueiog4zeHCsxz9w/fqOLF1ayapVq9i2bRtJSUksX76crKwspk2bhsvlYurUqbhcLqZNm0ZWVhbLly8nKSmJbdu21cg9Z84cysrKPGTJzc3lxx9/ZPfu3WzatIl169aRmppKTEyMh9zgXqaaP38++/fvZ+3atWzZsoVdu3axaNEicnNzPdotLy9n1qxZpKens2rVKhISEkhKSmLFihVkZmYyffp0nE6nxzVHjx5l2bJl7Nixg/j4eH7++WfS0tKYM2cOpaWlHnXz8vJYuHAhe/bsYePGjaxfv56UlBRiYmIoKCjwqFtUVMS8efNIT09n48aNbN26lcTERObPn3/K2JeXl/Of//yHw4cPs3LlShISEkhMTCQ2NpbMzExmzJhxitzZ2dksXbqU5ORktm7dyurVq0lLS2Pu3LmUlJR41M3Pz2fhwoXs27ePjRs3smHDBvbt28eCBQvIz8/3qFtcXMz333/PgQMH+OWXX4iLiyM5OZklS5aQk5PjUbeyspKZM2eSkZFBbGws27dvZ/v27cTGxpKRkcHMmTOprKz0uGb27NksWbKE5ORk4uLi+OWXXzhw4ADff/89xcXFTJ06FYfDwdSpU8nMzGTBggXs27ePDRs2sHHjRvbt28fChQtPkbukpIS5c+eSlpbG6tWr2bp1K8nJySxdupTs7OxTnrUZM2aQmZlJbGwsiYmJJCQksHLlSg4fPszMmTMpLy/3uCYnJ4fFixfXPGtr1qzhwIEDzJs3j6KiIo+6BQUFxMTEMGXKFNavX+/zZ0T3hQuxTZuGc/16xo4dy8GDB5k6dSrr1q1j9erVzJo1i7179zJ+/PhGPhI90W1Pw2oIISYBfwQurTpo2FC9MNxBGKtpBaRrtqexfTucey5l7duz0eA8CGlpafTWeFMMoLDQzqhRF5CXd/LYTt++RXz55VbsduPvL3/e0wgJCanJSeHtxq1e42o1mtIzIiKCiooKzfN4mIFWY9o9JoZOUtLuySfd7rcNYIk9DSshhJgAjASuasxgAEgpy6SUBdUv4ISmwvTpQ9m337L7pZc0bdYbUlP12aCOjnbyyCP7PcoOHIhi4cJuuvQXaNSXwKg5nj56javVaEzP6qCDgWAwQLsxzbjlFtIeeaRRg9FcAtpoCDcTgVuB30sp9zd1je60bYu84w5yL7rI8K6jo/U7jH/DDUcYONDTvn7zTT/y84P7LGdD2Gy2U1KiVlZWtshFVs9xtRIN6RkVFWXJoIO+oOWYnjih7ffegDYawCTgPuAe4ETVmZGuQghT3XvsdnPOMkRGRurWtt0Ozzyzz6PsxIlQvv22r259+hu1D905nc6alKi+nqXQc1ytRH16RkVFUVRUZMmgg76g5Zg6S0spXb8eNAr4GOhG4yncHlOrgSO1XneaJpGUhGzfTru4OITBp1MP6pSUpZpzzsnnyis9b8wFC3qwf79/7StoiRGH7vQeV6tQW8/qoINFRV4dt/I7tBzTs199lfBLL4XF2gTjCGijIaUUDbymmCmXuPRSznn5ZcKOHWu6soYMHz686Uo+8sQTKTgcJ9fjXS7BpEkDgi41rMPhqDl4p3cCIyPG1QpU61mdx9ufXWqbQssxLe7VC2ebNqCRgdXMaAghOgkhrhZCPCqEeFEI8Q8hxHNCiNuEEGdp1Y/fIwScdRZFAwYYPtNYtmyZ7n107VrGnXd6xrqJj2/Phg0ddO/bbGofuisvLzfslLYR42oFli1bVhN00CppWfVCyzFNefJJdq5dC888o0l7PrncCiH6A48AtwCnAwJ3IqZ83EmZ2nPSfTUfWAJ8J6Vc7oPMhqJHGBGA+Ph4zTeorEJJiY3777+AnJyTnsvdu5fwzTdbcDj0n3IY5XIrpcThcGCz2aioqAiojVgrEhERQXl5uV/FkLIKERERXHDBBQ2+r7vLbdWs4ltgK25j8Rnu0BztpZShUsqOUsoeUsoI3PkzBgCjcbuvfiKE2C6E0Cbkqp/icDiarqQxWocRaYiICBdPPOGZGjYjI4J583oa0r/ehIaGEh4ebpkERkaNq5lERUXxxRdfBI3B0HpMS0pKNLtHmz3TEEJcB7wPzAQmSSnzmnm9AP4CvAKsBF6RFnZ90GumsWfPHjIyMpquqCFOp9Mwzy0p4dlnz+O339rUlEVEVDJt2hbat9d32UaPmUZISEhN8iKrfXAZOa5GI4So2fAOZD3rorWufb/+ml6pqdg/+gjq2S/RbaYhhPgj8CxwmZTy7eYaDADpZo6UchjuGFDBEW2tNgsX0vuhh+j35ZeGdjt79mzD+hICnn3W0wW3pCSEyZP7GSaDr9Q9dGfVBEZGjquR1PWQClQ960NrXVvv3Il9/XrQIPlbc5enBgMjpZTHfe4ZkFKOA2YF3UZ5Xh7hGzbQeudOQ7v9/e9/b2h/p59+guuvP+JRtnx5V3bvtm4MyPoSGGlxlkJPjB5XI7Db7TgcDg+X2kDUsyG01jX99tvJ+fBD+MMffG6rWUZDSvk/KaWmi7dSyhVNBBAMPC6/nJLPPyf1sccM7TYpKcnQ/gAee2w/EREnQztIKZgwYaClXHCrDYW/JjAyY1z1xOFwYLfbT0nLGmh6NobWuuZecAE5118P/Xyf6Qf0OQ3L0rcvoQ8/zIkzzjC02y5duhjaH0D79uWMGuV5UOm339qwalVnw2WpTXUCI4fDgdPp1P0shZ6YMa56ER4ejsvlqncsAknPptBD15ISbXLf6W40hBAPCyHeEEK007svLRFCjBZCJANb9Gg/JCSkxqffKMzy8LnttnS6d/e8YT//vD8lJcZ/Z6nvLIW/zCgaIlBcfSMjIxsNOhgoenqD1roKpxPbb7/B0qU+t6X5U1vljlvTrpTyayAG+FgIcZXW/emFlHKSlPJMYIQuHaSn023LFqL27Wu6rkYcM/gEejUOh+Sppzz1zM4OZ9YsY8J5V88ohBCnJDAKBMwaVy2pDjrY2LgEgp7eorWu9uJizhk1Cm68EXw8Sa/lifCrhBA5QCZQIIRYWjXLaFO1Z3E/7pSvCoB336Xv88/TNTa26boacYbBy2G1ueSSYwwb5pkwcdasXmRlaT/bqv7gqZ3prry8PKAMRW3MHFctqA462BT+rmdz0FrXylatKO7VC9fQoZDbYOJSr9BypvFv3FFlnwWmAP2BL4EjQoiZuA1GRw3782/OPZfKwYMp7Wzc2v6aNWsM66suQsDo0SnYbCc/uMvL7Xz+uXapYUNDQ2tcZAHTD90ZhZnj6gvVeby9DTror3q2BD103fLddxSuXQs9fTtkq1nmPiHEJ1LKp+uUnY07LPktuA3Gs1LKmZp0aBB6He4D98bU5s2bNW2zMcrLy005iV6b8eNP44cfetQp+5WhQ/Nb1J7dbq/JSVF9hsKfM/e1BCuMa3OpPizZnBhS/qhnS9FL17POOouOHU/97m6ZzH1SyiQp5atSyjOqQov4lcHQm/DwcEJDQw3rzwqHox58cD+tWnkGapw4cSDNOTNXXwIjKx66MworjGtzqHZEaG7QQX/T0xf00lWLQI9aGo21Qog/a9hewCOEoE2bNtgMcvUcNWqUIf00Rps2lTz00AGPsr17W7FsWeOpYfVKYBQIWGFcvcWXPN7+pKev6KFr12XLaH/HHTBhgk/taGY0qmYRZwkhHtWqzYBn/nzOuO46zhg3zpDurBLYbuTIDPr29VzH/uqrfhQWesbaMSKBUSBglXFtiqioKEpLS1u8z+QvemqBHrqGHT1KxIYNsH27T+1o6T01AngI+FwIsU8I8aUQ4j4hRGCENtWDdu0Iycig1e7dhnR3ww03GNJPU9jtktGjPV1wjx93MG1aH4AaQ2FEAqNAwCrj2hBCCE3SslpdTy3RQ9ecSy4h/Z13fM6roeXy1PtAbNXPfcAdwHfAQSHEXiHEF0KIyzXsz/+56CLk4sUkGPQNyshN96YYPvw4F1+c41E2b15P0tMjlaFoJlYa17o010OqMaysp9booWvRgAHkXn89nHuuT+1oaTR2SCmflFK+LKW8HmgHXAy8DqQCdwPBs5PlDeHhiBtvpG23xtfztWLAAO3cW7XguecOEhJycqmistLGJ5/0N1Ei/8Rq41pNddBBrdKyWlVPPdBL15bsJdVFS6Ph8dVQSumUUm6SUo6VUl6H24icrWF/AUP79u0N6aegoFFPOkOoTmBks9no0uUEt9+e7vH+xo0d2bLFryLOmI4VxrUu1RkN6wYd9AUr6qkXeuhqKy0lPCkJNm70rR2N5AHY3Vg2PillpZQyp6H3gxaXi44zZnD2q68Skt+yswreouUD3Byq42xVp0WtvRl6330HadfOcynqk08GUlmpNru9xaxxbYjw8HCcTicVFRVNV24GVtNTT/TQNTI9nTMffBBuu82ndrQ0Gp8CtwohrtGwzcDHZsP+3Xd02LSJTjqfeO3d25hYT+CZwKj6LEV9XjNRUU4efTTVo+zgwSgWLOhulKh+j5Hj2hSRkZG65fG2kp56o4euldHRlHbpAr16+dSOlkbjeuAuYJkQYokQ4gUhxHkath+4vPgi+a+/Tu6FF+raTXx8vK7tV5+lqG0ovFlDvf76TAYNOuFRNmVKX/LzjTv46M/oPa7eUh10UK/QLVbR0wj00LW0a1fi5s4FHzfZtTQab+D2lpoNnIXbi2qrEOKYECJGCPGcEMI3E2cgeodG9+C++4h8803KdY5Ddf3112veZu0ERtV5KZp76M5mg2ee2etRVlgYytdf99VY2sBEj3FtDrXzeOuJ2XoaiV66amHQtTQae6WUf5dS3iOl7A2cBjwFLAcuAD4CNmjYn67oHhq9DqGhobpviM+fP1+TdmonMKo+dOdrXoqzzy7g97/P8ihbtKg7KSmBHzvKV7Qa15ZQ/aVBKw+pxjBTT6PRS1ctYg1qGbDwEuBxYBsQI6VMq/P+GUAHKeV6TTo0CD0DFnrgdJK3cCHHFy7k4AMP6NdPCxFC1MTJ0usMxdGjYdx//wjKyk6eDD/vvOO8//52mmOPgi1goVmEhobW5ChRWB9baSlD3nqL9m3bIhYsgFoBEU0JWFhlDB4ENgKn7OJIKXf5m8EwlNxc2tx1F/2mTKHVrl26dNGS0ATVMwopZU1eCr3o3LmMu+/2+K7Br7+2Y906FVG/McwIrxEWFlZzTxiFCiPiI1LSYdMmxLJl4INnW4iGIiHd0xb99wACkU6dEA8+SEFeHpWtWunSxc0339xkHSlljY99RUWF4d8i77zzEEuWdOPo0fCask8/HcAFF+TicAR+boyW4M24aklERIQpuUqM1tNM9NBVOhzs+vvfGTR4MDYfwq43a6YhhLhdCLFQCHGnECK86SsUzeLzzwmdMoWSHj2artsCVq5c2eB71QmM7HY7FRUVpiUwCg938cQTKR5lR45EMHeuCmHWEI2Nq9ZERkZSUlJiyr1hpJ5mo4eu0m4n88Yb4YEHwIeUDM0yGlLK74GxwGVAshDiWyHEtbVzgit8IyIigg4dOujS9tChQz3+rj50Z7ahqMtVV2Vz9tl5HmXTpvUhJyc4EvA0l7rjqheRkZGGbHg3hFF6WgE9dfU1SnSzP+yllBullM/g9o6ag3sfY68Q4iMhxPk+SaMAoJfdTt9vvkFovDSUmZmJ3W6vSWBUbSislsBICHjmmX0IcdJJo7TUzuTJKi5VfWRmZuravhCCiIgIUw0G6K+nldBDV1FZSfTu3YgdO8AHB6gWzxCqYkstllLeAwwF4oH/CCGShRD/FEKc1mKpghmXizY33UTf776jS2ysJk1Wn852OBxUVlb6RQKjQYMKueEGzwdnxYqu7Nypz36PP2Oz6TfRrw46WFJSolsf3qKnnlZDD11D8/MZ/uST7ii3Pjz7mkgmpSySUk6tim57JXAcmCaE2Fx1qE/fU2uBhM2GeOYZKkaM8GlvozqBUWhoaM3p7DZt2ljaUNTlkUdSiYz0PFE+YcJpWGAFzVLotZwZGhqKzWbTJEWoFuilpxXRRVeXi7LOnaGnb/uDmpszKeVRKeXHUsoLgPuA9sAvQohlQohRQgjlON8Uf/0rIZs24bz00mZd1lSmu106ufLqRfv2Fdx//0GPsp07W7NyZReTJLImeoxrWFgYLpdL86CDvuBv968v6KFreadOxMfEwMGDTVduBF3ne1LKvVLKf0kpBwP/BM4HdgghZgohbhJC2JtoIjix2xFC0LdvX6+qh4WF1fjNN3Y6+/LL/S8H1q23ptOzp+da+hdf9KekRN061Wg9rnoGHfQFf7x/W4peutrtvj83hi0SSik3SymfAwbijlF1F7BPCDFJCHGBUXL4Ex3atKHfTz/RdcmSU96rdpEVQlBWVubVEsKPP/6oh5i6EhoqefppTxfcnJwwZswInoinTaHluFZ7SGkVKUJL/PH+bSl66aqF0dAsjEiLOhciErgFuFxK+YRpgjSCYWFE6mPGDLj3Xipat2bj1KnYO3TAZrNRXl5uyYdaL6SEl18eytatJ2NzhYa6+O67LXTtemreARVGpPlYxUNKoR9tEhPpHxNDm6uvhtde83jPlDAiLUFKWSylnG5Vg2E6d9wBF19MxWmn0XfmzBoX2ZYaDH8NwyAEjB69D5vtpN4VFTY++0y54ILv41rtXWd1g+Gv929L0EPXiMOHabN6Naz3LZpT8PiwNRNDQ6M3REgIvPYakfHx9Pz+exw5viU+/Mtf/qKRYMbTt28xf/7zYY+yX37pTEJCW3MEshC+jGtoaCh2u90vsuL58/3bXPTQNf+sszjy1lvw1FM+taOL0RBCnCGECNOjbaMwOjR6g9x4I/z1r2R98w3lHX0L3Ldw4UKNhDKHBx88QOvWnt48EyYMxGL7tYbT0nGtHYjSH/D3+7c56KFrSa9eFN97L/zpTz61o7nREEI8BwwGooUQs4QQW4QQ64QQ/xZCdNO6v6Dgo4/ofPfdhIX5Zocv1DkzoN60alXJww/v9yhLTY1m8eLgvq1aMq4RERFUVlZ6lVnRKvj7/dsc9NI11IeYU9XoMdPoK6WMAf4BvC+lHAFcgTsB05dCiBd16DPgsdvt9O/fH1tJCaKFD3pKSkrTlSzOTTcdoV+/Qo+yr7/uR2GhpgGb/YrmjquZQQd9IRDuX2/RQ9eolBQi0tLAx6VIPYxGdVQ5m5QyDmpCjiyTUv4ROCGEGKtDvwFPlw0buPDBB+newqmr4d5fOmC3S0aP3udRlp/v4Ntv+5gkkfk0Z1zNDjroC4Fw/3qLHrqe8c47dLrsMli1yqd29DAavYUQXYFNQoia3ZzqSLhSyi+AQ0KIP+jQd2Bz9CiOo0fpumxZiwKO+bq8ZRWGDcvjssuyPcpiYnqQlhZpkkTm4s24BoJLbaDcv96gh67O8HBk69bQzbflXD2MxixgHrAYqKgKXvgC8Ggtw/EJcI0OfQc2jz4Kn3xCxty5LQo4lp6eroNQ5vDkkymEhp5cXnE6bXzyyQATJTKPpsY1JCTEMkEHfSGQ7t+m0EPXhI8/xnnsmDtgoQ/oEXtqBpALrAW2SinfklJ+KKX8QkpZexH1uNZ9Bzx2Ozz1FP0GD27Ryc7f/e53OghlDt27l/KXvxzyKNu8uQObNrVv4IrApbFxdTgcNVED/J1Aun+bQg9d7XY7ISEhPkW4Bf3OadwBFOCOM/WCEKK+XUr/CbdqMRwOB3369KF1M+PiL1++XEepjOfee9No397zw/CTTwZSURFct1ZD4xoeHk5lZaWlgg76QqDdv42hh67h4dokW9UtjEjVUtTfcXtRlQLTgXVAFnARUCml/EiXzjXE1DAiDSEl8rbbEDEx7Pj3v8m57DKzJTKNZcu68M47gz3KHn30V7788ndBHUYkIiKC0tLSoAo3o2iYjmvX0mf5clrdcw88/fQp71sijIiU0iWlfAfoCvytqq97gIeAnUYaDCHE00KI/UKIUiFEvBDCvz9lhUCceSYyNJTwI0e8viwQwzBce20WZ5zheY/PmBFcnlR1x7XapTbQDEYg3r8NobWurXbvptX69bBjh89tmRqw0AiEEHcCU4GngfXAE8CjwJlSyjQvrrfeTAOgsBAyMtheUsLx495tDzmdTk2iXFqN335rzTPP1F4DLgDaBM1Mo3pchRCEh4f7/YZ3QwTq/VsfWusatX8/A3NyaHfhhVDPyoQpMw0hxFAhxHgLutL+DfhKSvmllHKnlPJ54BDgWwAWs4mOhkGDGDhwoNfZ+GbNmqWzUOYwZEgB11xzak7l1NTANxjgHtfqoIOBajAgcO/f+tBa16J+/ZCPPFKvwWguWh6j/RfwR+BewLcgSRohhHAAw4Bxdd5aAVzcwDVhQG0n6VbgtsRWpV1JCRWLFpF5ww2N1rvooosoKioySCpjuffeJNasCaOszI57pgFz57ahf/8McwUzgMsuu4ySkhJL36NaEMj3b1300NXpdDZ4jzTn3tFseUoI8QkwB+ggpZynSaM+IoToDhwGLpFSbqhV/g/gASnl6fVc8y/cWQYVCoUi2GhyeUrLmUYZsFtK6f3OrHHUtYyinrJqxgIf1Pq7FZB+6NAha+1p1OaJJyjbs4cdjzxCcZ+GN4F/++03hgwZYqBgxlJWZuM//xlMQkII0IuXXorlyiv9/3xCQ4SHh1NeXk5SUlJAj2s1gX7/1kZLXdtv3syQ//wHhgyBDRvqrVNQUECvXr28ak9Lo/EW8IUQ4m0p5XYN2/WFHMCJ24OrNp1xu/6egpSyDLcBBKjZL2jdurV1jcZXXyHDw4mOj0c0MqV1OBwBvTEcFQXvvXeQnTslzz0HV15ZFrD6VseQioiICPhxrSZY9ARtdW3lcBB22mmEnX8+aPAZpqXL7SW4o9nGCSFWCSHeFEJcXrWvYApSynIgnlNDllyDO+puYBAZibDZGDhwYKPVcnNzDRLIPEJCJP37B/a6d90YUsEwrhA8eoK2uuZceik5P/8MX3+tSXtaGo3/Az4GvgQ64N4X+BnIqzIibwghztCwP2/5AHfcq4eFEIOFEB8CvYHPTJBFV9pFRzNo+XJ6NeB5MWjQIIMlUmhJQy61wTKuwaInaK9rdHS0z+FDqtHSaOwHxkopn5ZSnoN7CeguYAru5aG3AN9i8rYAKeVs4HngTSABuBy4UUp50GhZdCc2lu7jxtHv668JyzzVBXXdunUmCKXQArvdTmhoaL1pWYNlXINFT9BQVykRQriNhkZo6T01FPgrbm+leXX3NYQQXYAuUspETTo0CMse7qsPKeHOO8k++2ySL74YWedwUHl5OQ6HaauFhlFUVMRNN90UMIf7HA4HLperwSx7wTKuwaInaKdru/h4zvjf/wi780746KMG65lyuE9KmSilfAT3zOKUxAZSyix/Mxh+hxAwZw7tXn2V0IiIU96ePXu2CUIpfCEsLKzJtKzBMq7Boidop2urnTvdqw5Z9fr9tIhmzTSqZguFUkpNdxqFEKdLKXdr2aZW+NVMoxZHjhxh986d7j9suoUYsySBMtNQQQcVvmIvLmZwUREd+/aF889vsJ6eMw0X8LUQwrfUT7UQQowCntOqPYWbrocOMfz55+m2eHFNWTAFfPN3mhN0MFjGNVj0BO10dUZGEnHjjY0ajObSLKMhpcwGXgHmVHkjtfgrrBCijxDiK+D3KKOhOWLTJqKTkug9cyY4nQD88Y9/NFkqRVO0JC1rsIxrsOgJ2ukaEhJCZKS2aZCb/aEvpTwA3ABcCOwSQrwqhDhHeBE1TwjRWghxoxBiKhAHbJZSPiSldDZXDkUTPPUUvPIKh2bMcGf8AzY0cBpUYQ1sNluL0rIGy7gGi56gja7ttm5lQEwMYtcuDSQ6SYtOhEspC4HHhRDn4U609AZQKYSIw+09lQfkAw6gfdWrH3A2kA18hTs0eY6vCigaIDQUxo2jd2kpmVu24HK5mjz8pzCPkJCQFqdlDZZxDRY9QRtduy5bRpeVKyEyEv7zHw2kcuNTGBEp5a/APUKIVrhPWV8EDMZtHKJxh/DIAw4A3+Nehlov1c6eYYSHh9O7d2+ObN4c8FFQ/RVvPKQaI1jGNVj0BG10Pf6739EhNJSQ66/XQKKTaBJ7Skp5Aphf9VJYjN6TJtH7ww/Zf/fdMGKE2eIoaqGFh1RLZif+SLDoCdromjNyJKe/845mJ8GrCS5fzGYghBgthEgGtpgti6/YpMTmdHJOPafEFeYRERGhSVrWnj17aiSRtQkWPUEbXdu1a+d1grbmoIxGA0gpJ0kpzwT8/6v5G2/AihV897vfNV1XYQjVBkMLfv31V03asTrBoif4rmur5GTah4drJI0nWoZGV1iVNm3gmmt4dvBg9u3bZ7Y0QY0QQvO0rNdee61mbVmZYNETfNPVVlbGuX/7GzYhIDkZ+vXTUDI10wgqpkyZQve2bemwcaPZogQljQUd9IWYmBhN27MqwaIn+KZreEYGle3bIzp1gr59tROqCs0CFgYq/hpGpF7y8pDnnAOHDhH/xRcUBqgLoxXDiISGhiKlbLGHlELRHHr36kX/qCjoWjf/XP2YErBQYX3GTJyIuPBCnN27ExJE7otmExYWhtPp1M1gBEt4jWDRE3zXtWOnTl4bjOai9jSCiGeeeQYqK7GFh1P222+g4bq6on6MCDp4yy236Na2lQgWPaHlugqnE0dkJK1atdJYopOomUYQMW3aNOjYEVt0NAMGDDBbnIBHK5fapoiNjdW1fasQLHpCy3XtO2UKvxs1CjFzpsYSnUQZjSDiiiuuqPm9Y8eO9N61i661ouAqtKE66KCWHlKNcc455xjSj9kEi57Qcl07rF9PWEqKxtJ4ooxGEJFS+2b6+Wf6P/UUp02ciCM72zyhAgybzUZoaKhhBgMgIyPDsL7MJFj0hJbruuPzz5HffgsjR2os0UnUnkYQ4ZE+8sor4corKejbF2c9Wf4Uzac66GB5ebnh/QYDwaIntFzXDn37Iq65RmNpPAmeUVDQpUuXk38IAT/9RLTLhW3LFpwVFeYJFgA4HA5dPaQao127dob3aQbBoie0XFePZ1wn1PJUEBEfH+9ZUHXYrK8OB4CCifDwcCoqKnA6zUkLs2fPHlP6NZpg0ROar2vbX3/lnDfeoPXPP+sk0UmU0Qgibr755nrLuwNnvfsuXVasMFagAMAKebwvueQS0/o2kmDRE5qva9fly2m3bh0sXaqTRCdRRiOImDx5cr3lYsYMOi5dSv/JkxFqmcprjPSQaozFQeIBFyx6QvN1PTBqFJUvvggPP6yTRCdRYUQaQAgxGhiN27CeHhBhRBqirAwefJCUm2/mkE6nSI1E7zAiQggcDkdQ5XdQWJt27dr55JKswohoQECFRq9izJgx9b8RFgYzZ9L95pux2dQt0Rh2u52QkBBLGYxgCa8RLHpC83Xt1q2bTpKcivqECCJefPHFRt+PiIigZ8+e2Cyw5GJFQkNDAaiw2BLeHXfcYbYIhhAseoL3urbavZtBkybR0cAEa8poBBGffPJJk3X6rljBRXffTdtt2wyQyH+oDjpolodUYyxYsMBsEQwhWPQE73XtvmAB3b//Htv77+ss0UnUOY0gYqQXp0Rtv/2GLT+f7j/+SJ7K9Ae4XWq1zoGhJRdffLHZIhhCsOgJ3ut69Oqr6RwWhv2JJ3SW6CRqphFEbPNm9vDWW8hJkzg0bpz+AlkcKWWNS62VCZbzC8GiJ3ivq/3aa7HHxICBBlUZjSCiU6dOTVdq3x7x9NOcNniw/gJZHKu41DZF27ZtzRbBEIJFT/Be1549e+orSD0ooxFEhIWFeV23devWdOncmegg+nZXjc1mIywszPIzjGqqN+gDnWDRE5rWteMvv9B/4ULa2u0GSXQSZTSCiJTmhEwuKuL0Rx9l2NNPE3nwoH5CWYyQkBBsNpulXGqbIliivwaLntCEri4X/b75ht4ffghffmmcUFUooxFEXHnlld5XjorC1rkz0uEgSuf4/FbB4XDgcrn8Lo/3ueeea7YIhhAsekLjugopybrrLuT55xtyArwuymgEEc0+HDVhAuzezYkbb9RHIAtRHXTQ5XKZLUqzCZaMdsGiJzSuq7TbCXv2WcSWLeA+xW0oKoxIEwghWgP5gRBGxOVytejE99GjR0lOTtZBIn1obhgRf9nwboiWjqu/ESx6QuO6OhwOLrzwQk3/FyqMiKJe/vvf/7bous6dO9MlP5+O69ZpLJH5hIeH+7XBAJg+fbrZIhhCsOgJDes6YNIkTktNxSaEwRKdRM00mkDNNIC4OOQll+B0ONg8dSoVFk+G481MI5CCDgbLN/Bg0RPq17XN9u2c9/zzyNBQxP790KOHZv2pmYYGCCFGCyGSgS1my6IVY8eObdmFw4YhzjmH0uHDsRmcylQPrBh00BdmzJhhtgiGECx6Qv26FvfqxYlHH0WMHq2pwWguKoxIA0gpJwGTqmcaZsujBaNGjWrZhTYbrFyJIzwc55Yt4GfeRbUJCQlBSmm5oIO+cI3OOaGtQrDoCfXrKrp2Jerzz93Po4momUYQsXr16pZf3Lo1DoeDPn36aCaP0YSFheFyuSwZdNAXEhISzBbBEIJFT6hf1759+1piec58CRSGMWDAAJ/b6NG+PQNnzaLjmjUaSGQc4eHhlJWV+aVLbVN0797dbBEMIVj0BE9dO61axZD33qObRb7sKKMRRGixhm/77DN6fv45AydNwuYHewJSSstHqfWVQFpqa4xg0RNO6iqcTvpPnkynJUsQ335rslRu1J5GEJGdne17I089BT/8wNGbbsLlcPjens4EusEAyMvLM1sEQwgWPeGkrtJuJ+2ddzj9xx/hb38zV6gqlNEIIn6nRX6MyEhYu5aOxcXsj4vDyi7bgeJS2xSDBg0yWwRDCBY9wVPX7rfcAvffb6I0nqjlqSBi4cKFmrUVGRlJjx49wIJ7BCEh7u9C5QHgHuwNGzZsMFsEQwgWPcGtq620lK5du9KqVSuzxfFAGY0g4umnn9a0vX6pqZz/2GO0TkrStF1fCA0NDcjN7sa4+eabzRbBEIJFT4D7zziDi+66i4E//2y2KKegjEYQ8b7GeYTtc+cSlZpK3+++07TdlhIeHk5lZWXQGY05c+aYLYIhBIueAM4JEwjNzydk0yazRTkFFUakCQIpjIjmHD2KfOcdEm68kXwTksHUpvaGd3MDFioUViM6PJxhW7Yg7rwTunTRvT8VRgQQQvQVQnwlhNgvhCgRQqQIId4SQljf5UcnxowZo22DnTsj3n+ffsOGadtuM/GnLHt60OyQ935KsOgJsGDxYsRzzxliMJpLwBoN4Azc+j0BDAFeAJ4EWhbqNQB47LHHdGm3bdu2dOrUCcexY7q03xA2my1oPKQa449//KPZIhhCwOspJe03bqRHt24888wzZkvTIAFrNKSUy6SUD0kpV0gpU6WUC4H/AbeaLZtZLFiwQJ+Gy8o4/Z13uPCee4g4fFifPupgs9mw2+1B4yHVGOvXrzdbBEMIdD07r1zJ0H/8g4HPPceCmBizxWmQgDUaDdAGyG2sghAiTAjRuvoFWMvfzQeG6bWM5HAQcvQotvJy2m/erE8ftQgNDQWC64RwYwTL+YVA19NWVoaMjERceinDzj/fbHEaJGiMhhBiAPAs8FkTVV/FHdW2+pUOkJqayrvvvktZWVnN3sCYMWPIyMjgq6++Ii4ujuXLlzN//nySk5P5+OOPKSgo8Kibm5vLpEmTSExMZOHChSxevJht27bx+eefk52d7VG3uLiY999/n927dzNnzhxWrlzJhg0bmDJlCocOHfKo63Q6GTt2LAcPHmTq1KmsW7eO1atXM2vWLPbu3ct7771HaWkp48ePr7kmMzOTyZMnExcXx7Jly4iJiWHHjh1MmDCB/Px8j/bz8vKYOHEiSUlJ/PDDDyxZsoT4+Hi++OILsrKyGPP22/Dxx3z12GMcvvVW5syZQ3p6OqtWrWLbtm0kJSWxfPlysrKymDZtGi6Xi6lTp+JyuZg2bRpZWVksX76cpKQktm3bxqpVq0hPT2fOnDmUlZXVrGVPnTqVwsJCYmJi2LlzJ5s2bWLdunWkpqYSExNDYWGhx7p3UVER8+fPZ//+/axdu5YtW7awa9cuFi1aRG5urke75eXlzJo1q0buhIQEkpKSWLFiBZmZmUyfPh2n0+lxzdGjR1m2bBk7duwgPj6en3/+mbS0NObMmUNpaalH3by8PBYuXMiePXvYuHEj69evJyUlhZiYGAoKCjzqFhUVMW/ePA4cOMCaNWtq5F68eDE5OTmnyP3DDz9w+PBhVq5cSUJCAomJicTGxpKZmcmMGTNOkTs7O5ulS5eSnJzM1q1bWb16NWlpacydO5eSkhKPuvn5+SxcuJB9+/axceNGNmzYwL59+1iwYAH5+fkedYuLi/n+++85cOAAv/zyC3FxcSQnJ7NkyZJT5K6srGTmzJlkZGQQGxvL9u3b2b59O7GxsWRkZDBz5kwqKys9rjl48CBLliwhOTmZuLg4fvnlFw4cOMD3339PcXHxKXIvWLCAffv2sWHDBjZu3Mi+fftYuHDhKXKXlJQwd+5c0tLSWL16NVu3biU5OZmlS5eSnZ3tUdfpdDJjxgwyMzOJjY0lMTGRhIQEVq5cyeHDh5k5cybl5eUe1+Tk5LB48WJ27drFli1bWLNmDQcOHGDevHkUFRXV1P2sooL8dev4tE0bNm3aZOhnRPVngzf4nfeUEOJfwD+bqHa+lHJrrWu6A78Av0gpH22i/TAgrFZRKyA9ELynfvjhB/785z/r2kdWVhY7d+7UpW1vQ4IEm/fUunXruPTSS80WQ3cCWc/w8HCGDx9eczDViGe1NoHuPTURGNzEa0d15SqD8TOwEXi8qcallGVSyoLqF3BCcw1MQosot03RuXNn2glBGw0P/AVD0EFfCJbor4GoZ1hWFkP//nfOdDhqDAYY86y2FL8zGlLKHCnlriZepQBCiB7AamAb8JCUMrhOfdXhl19+0b0PsXMnZ992G2e99hoh+drkrlIGo3G2b99utgiGEIh6nvbxx7TfupXWr7ziUW7Es9pS/M5oeEvVDGM1cAh4CegkhOgqhOhqqmAmct999+nfyaBB2Lp3x9WpE2E+uuDabDZCQ0OD3qW2KYIlo10g6pnxj38g//hH+Pxzj3JDntUWErBGA7gWGAj8Hvdm9pFar6Bk4sSJ+ncSEgKLFsH27ZSedpoPzYQghFAeUl4QY2H3TC0JND3tdjun/f73iEWLoH9/j/cMeVZbiN9thBuNCiPScg4ePMj+/fubfZ3D4fAphlSwbYQr/IvQ/HzCMzLoeeutdLHIie9A3whXtBDNw4g0Qa+ePem+bh3t4uK8viYsLIzy8vKgCzroC8ESXiMg9HQ6GTxmDL/761/psnJlg9WMflabgzIaQcSzzz5raH+2yZMZ9MYbDProI4QXJ7fDwsLU/kULuOWWW8wWwRACQU9bRQWidWsIDYWhQxusZ/Sz2hyU0QgivjM6hPl998HAgeSNHIloomp4eLgyGC1kxYoVZotgCIGgp71VKyIWLULExcFZZzVYz/BntRmodK9BxFVXXWVsh9HRkJxMq7IyZHw81LN/JoQgNDRUudT6wHnnnWe2CIbgz3rai4txRUVx5plnEh4RAWee2Wh9w5/VZqBmGkHE3r17je80NJTo6Gi6det2ylt2u10FHdSA9PR0s0UwBH/VMyQ/n+GPPcZ5c+bQzktnGlOeVS9RRiOIiIiIMK3vfk4nZ7/xBq127wZOBh2srKw0TaZAISwsrOlKAYC/6tlxwwYiMjJovXQpFDTqmFSDmc9qU6jlqSCiQ4cOpvUd+u67dFi3DvuJEyR/8gkVFRUod29tCBZXcH/Vs/Tuu3ENHYrt3HOhXTuvrjHzWW0KNdMIIhISEszrfMwY5J//TNGll9I6NlYZDA3Zt2+f2SIYgt/pKSWRkZEMGTIE2733wpAhXl9q6rPaBGqm0QBCiNHAaALIsP7pT38yr/MePRA330yPhx6iU7t25A0bRmV0tHnyBBAXX3yx2SIYgj/p2XXpUjpt2kRUTEzNUmxzMPVZbYKA+UDUGinlJCnlmcAIs2XRii+//NJcAe65B4YPJ/ehh3CGh5srSwCxePFis0UwBH/RMzQ/n4ETJtBhzRrCZ89uURumP6uNoMKINIEKI6IxLhdlFRXExcXptgmuwogozEQIwXnFxbResgQ+/hhEU6eUzEeFEVHUiyVCE9hshIWF0b9/f3C5ECogoc8ERHgNL7C8nlIihODMM8+k9Q03wIQJLTYYlnhWG0DNNJogkGYapaWlhFtkWUimpFB0550c79ePlNGjNW072GYa5eXlOBwOs8XQHSvrGb1vH6eNH0/pd9/RRYNDiEY/q2qmoaiXCRMmmC1CDWL3bqLj4+m2dCmheXlmi+PXzJ8/32wRDMGyekrJ6e++S5sdO+jywQeaNGmlZ7UuymgEEUbmHG6SG2+EDz4gb9UqKtq2NVsavyZQ82bXxbJ6CkHhlClw++2gUR4MSz2rdVBGI4iIj483WwRPXniBjiNG1BtiROE9e/bsMVsEQ7CanqLKkWPgwIF0u/himDsX3Es8PmO5Z7UWymgEEV27WjPT7cCBA+l45Ahd/cSl0mq0b9/ebBEMwUp6Ru/ezYgHHuDMkhJ69uypeftWfVZBHe4LKkJCrDnc9tRUhjzyCJSVUdq1K3nDhpktkl9ht9vNFsEQrKRn3+++IyIjg4gJE+CGGzRv36rPKqiZRlDRktSrhjBwIOLuu6m47DIKBw0yWxq/48iR4Eh7bxU9hRAwbRo89xzMmqVLH5Z9VlEzjaDi8ssvN1uE+hECPv0UhxD0ycwkJSXFbIn8iqGNZIALJMzWM/LgQcr692fIkCHupbLx43Xry7LPKmqmEVTMmDHDbBEaxuGA0FB69epFt27daLN9OzaVZ8MrVjaSazqQME1PKenz3Xec/9BDDEtNNWRvxcrPqjIaQcTLL79stgheMeinnzj3hRcY/J//IJxOs8WxPHfddZfZIhiCmXpGFhQgpCQyNdWQ/qz8rCqj0QBCiNFCiGRgi9myaMXYsWPNFsErxJAh7plHt25IP4jbYzZW/laqJWbp2bZdO9pPmwaLFsFbbxnSp5WfVRVGpAkCKYyIX5GcTMXAgWxPTKSwsLBZlwZbGBGF9kSlpNDlp58ofvNNBp1+OjZbYH+/VmFEFPVi5SBop3DmmYQ6HAwdOpTIiAg6/fwzuFxmS2VJLB/ITyOM0jPkxAnO++tf6T1rFmesWWOKwbDys6qMRhDxwAMPmC1Cs3E4HAybOZMh//43gz76yGxxLMl1111ntgiGYJSeon17yt54A665Bv7yF0P6rIuVn1VlNIIIf/WysQ8fjgwJofycc8wWxZJs27bNbBEMQU89Q/PycOTkEB0dzbBhw4h6+WVYuhRMOoVu5WdVndMIIgb568G5++5DXHwxPXv3JjcxkYKCRpdcgw49wlhYEb30bJ2czJA338TZrx/h69djqw6/buIJdCs/q2qmEUSUlJSYLULL6d+fkJAQzjnnHDpGRnLahx+qkOpVlJWVmS2CIeilp7NdO0LLyogsKsKWk6NLH83Fys+qmmkEEceOHTNbBJ+x2+0M+eQTxMKFtNq7l22TJvlFOk09CZaZl5Z62srLcTkctGrVijMvuABb374wdChYxNvOys+qmmkEEecEyJ6A+L//g8GDKRkzJugNBsCAAQPMFsEQNNFTSrouW8aFd93FgNJSzjvvPCIiIuCiiyxjMMDaz6oyGkHEokWLzBZBG844A5KS6HLrrZx99tmEhITgyM6GID1ztGnTJrNFMASt9Oyydi2O48fpNX++Zc9fWPlZVYf7miCQDvcVFxcTGRlpthiaU5yaSsiFF5J/xhnsevllCoQIqsN9Vsr9rict1lNKhMuFtNvp3r07AyIisM+eDc8/DxYNQW70s6oO9ynq5QON8hdbjcjERELz8mh15AgyNNRscQxn7ty5ZotgCC3RMyI9naF//zsDp03j3HPPZdCgQdh79YKXXrKswQBrP6tqptEEgTTTCGi2boWICI60b8/27du54YYbWLRwIVGtWpktmcJEOq1dy5A330S2bo04eBBUPvp6UTMNRb1YOTSBzwwfDkOG0K1bN84991wAzn71VSIOHzZXLgNQYUQ8CcnPB6Bt27b0feEFeOstRHy8XxkMKz+raqbRBIE008jOzqZTp05mi6E7BdnZtOncmXwg57HHSLvnHrNF0pW8vDza+tEHYktpSk9HTg6n/+9/RB06RP7GjXTp3ds44TTG6GdVzTQ0IBBDo8+fP99sEYwhLMz98/HH6fzOO0RHR7v/DtDcHGvWrDFbBENoSk9XdDRt9u8n7OhRuvh59kcrP6vW3QkyGSnlJGBS9UzDbHm04PzzzzdbBGN57z2iW7dm2LBhpB86RPRtt1HUuzf7H3oIZ7UhCQDOOOMMs0UwhLp62ktK6LBuHUevuYYuXbrQr18/QmbMgN694bTTTJJSG6z8rCqjEUQcOXLEbBFMQQhBr4MHYetW2iQlkf6XvwSU0bDy6WEtqa2nrayMEfffT1hODn0vv5zIwYPdb/zhDyZJpy1WflbV8lQQ4QzQ5RmvuOwyWLEC24QJnH7ddTVLVpEHDvj9oUBXkOQZcVVW1vzeuksXuOkmGDiQSIse0PMFKz+raqYRRPTt29dsEczlmmsAaAcMGzaM7E2b6PjYYxSceSZJ48bhjIgwV74W0rVrV7NF0Benk94zZzJswQLSLrqIHsOHuzfEx4+HyEhLn7doKVZ+VgPPRCsaZO3atWaLYBmEEHQ+eBBhtxPWvj02P/aMS0xMNFsEfbHb6fLrr7TKyWHIL7+c9KBq3TogDQZY+1lVLrdNEEgut7m5ubQ3KamMkVS7D3o1ZocOQXk5zr59OXz4MIf37GHgW2+RMXIkx4cP94uAiAUFBX5/b9bgctE+Lo4uK1aw5//+jy59+tCzZ08it26lcOdOoh9+GILg1L/Rz6pyuVXUyyeffGK2CNajVy8YMAC73U7v3r25ICmJTuvWMWjCBL/Z61iwYIHZImiGkJLTP/yQLqtWcVF6OoMGDXLHYLr8cj7Kzg4KgwHWflbVTKMJAmmmESw0a6ZRl/R0+N//kGeeSc4tt5Cenk5+Xh69p0/n6FVXUdqjhz5CByHC6aTdli20j48n5ZlnaN+hA927d6f9tGnukB+jR0OQhH03GzXTUNSLlUMTWIaePeGjjxCPP06nTp0477zzuKCwkP5ffcWwJ59EVFSYLeEp+GsYEVtJCWe99RY9583jIoeDs88+mw4dOiD++lf44INTDEYw3b9W1jUwd5EU9fLcc8+ZLYJfEtG1K1x/Pfa+fRl01llkZmaSn59P7+nTKerXj9zhw5HVeaVN4NZbbzWtb28Jz8yk59y52EtK2P/663Tu3JkuXbpge+IJsNtxdOvWZBvBdP9aWdegmGkIIcKEEAlCCCmEONdsecxiypQpZovgn5x/Pixdim3SJLp168Z5553Hhd270++rrzj7tdcIrZ2G1ITl3mXLlhneZ1PYi4sJzc2t+TtMCHrOn0/X2FguGjyYgQMH0qpVK7fb7AcfQL9+TbYZTPevlXUNlpnGu0AGYN0cigZw9dVXmy2Cf1PrEFl4eDg88wwcOcK5N95ITk4OOTk5dP3Xvwg/coSD999P/tlnGyLWsGHDDOnHW3rOmUP/L74g5/bbKR43jg4dOrgNREICwoe0qsF0/1pZ14CfaQghbgCuBV4yWxaz2bVrl9kiBA69e8PHH8PcuURGRtK7d29+N3Qo3davp/3WrbRr04bQKk+fiEOH6P7DD0SlpuoiSlpami7tekPv6dM5529/IyI9nfDwcLp160bHESOwOZ10Pn6cvn37ug0GwNixMHIktHApL5juXyvrGtAzDSFEF2Ay8Geg2MtrwoCwWkUBk8UnGFKfmordjli3DmJj6XvXXfQJDaW4uJjK5ctpM348x666iqQ336yp3mHjRkq7dKG4Tx+k3d7ibo1I9erIzqZrbCyiooKDDzyAEIKoqCi6JiQQ+euv/K6oiNALL3RX7tXLHbZFY8+nYLp/raxrwM40hBACmAJ8JqXc2oxLX8Ud1bb6lQ6QmprKu+++S1lZWY1nw5gxY8jIyOCrr74iLi6O5cuXM3/+fJKTk/n4448pKCjwqJubm8ukSZNITExk4cKFLF68mG3btvH555+TnZ3tUbe4uJj333+f3bt3M2fOHFauXMmGDRuYMmUKhw4d8qjrdDoZO3YsBw8eZOrUqaxbt47Vq1cza9Ys9u7dy3vvvUdpaSkxMTE112RmZjJ58mTi4uJYtmwZMTEx7NixgwkTJpCfn+/Rfl5eHhMnTiQpKYkffviBJUuWEB8fzxdffEFWVpZH3ZKSEt5//3327t3L7NmzWbVqFevXr+e7774jLS2Nt99+G5fLxZgxY3C5XLz99tukpaXx3XffsX79elatWsXs2bPZu3cv77//PiUlJR7tZ2Vl8cUXXxAfH8+SJUv44YcfSEpKYuLEieTl5Xl4neTn5zNhwgR27NhBTEwMy5YtIy4ujsmTJ5OZmenRbmlpKe+99x579+5l1qxZrF69mnXr1jF16lQOHjzI2LFjcTqdHtccOnSIKVOmsGHDBlauXMmcuXPZbbPxfkUFxZWVvP3220RFRbFq3z7Kr7yS37p1IywsjIyMDFITEznrtdc4/5FHKD5ypMYDasu779J+3jw2fPEFBw4cYM2aNWzZsoVdu3axePFicnJyaupOnTqV8vJyNm3axOHDh1m5ciUJCQkkJiYSGxtLZmYmM2bMwOl0elyTnZ3N0qVLSU5OZuvWrfyyahVpaWnMnTuXkpISUt58k8H/+Q+/jh1Lfn4+P/74I3n79tF/8mR6zJ6N3WZjx44d9OvXj6UDBsBnn/Hp3r2cOHGCjz76iOQDB5iXkMCKFSvYvHkzX3/9NRkZGR7/u4qKCt555x1SU1OZPn06a9asYc2aNUyfPp3U1FTeeecdKioqTrnPv/76azZv3syKFSuYN28eycnJfPTRR5w4ccKj7rFjx/j0009JSEjgxx9/ZNGiRSQkJPDZZ5+Rk5PjUbewsJAPP/yQXbt2MXfuXGJjY9m0aRPffPMN6enpHnUrKysZN24cBw4cYNq0aaxdu5bVq1czc+ZMUlJSNPuMOHbsmKGfEePHj/f6A9LvzmkIIf4F/LOJaucDFwN3ApdLKZ1CiL7AfuA8KWVCI+3XN9NID4RzGl999RWPPPKI2WLojk/nNIwiJQVGjYKcHCp++43i4mKKi4uJ+utfaT13Locff5z9o0ZRWVlJSH4+5z3/PBVt2pDwwQc1eyvtN28m6sABFp84wdBHHwXAVlpKx6oQFEerYm0BdPzlF9rs2EHuiBEcrwq7HZaZybCnn8ZWXs7GpUsJCwsjPDyc3m+/Tds5cyh+6SXEv/9NeHg4orISHngAzjkHnn/+ZM4SAwmW+xeM17U55zT8cXlqIjCriToHgNeBC4Ey4RkKYqsQYrqU8oH6LpRSlgFl1X8LPwgj4S033HCD2SIoqhkwADZsACAUaNOmjfuhvfRSKCqixx//SI9LL8XpdFKxYwfhBw7gat2a0wcPpqKigsrKSjpOmEDr+fP54+OPU9qpE1JK7JmZDP7vf5EhIbjuvhubzYbNZqPr7t20/f57WvfsSadBgwgJCcExYACO48cBuOzcc6F67+Hhh2HYMCKvugqqgziGhsKMGcb/n2oRTPevlXX1O6MhpcwBcpqqJ4R4DrfhqKY7sBz37GOzPtJZm6+//prXX3+96YoK83juOferCrvdjn3AAIiNxVZaSrfa5xluuAGio9lSUsLdQ4a4y7p3h2uuQYSEcNZZZ52se8cd0KcPbf7wB9p0736yPCEBunWD2vlFrrvO/bIYwXT/WllXv1ueaineLk/Vc50KI+Jn+MXylEJhIVQYEUW9WDk0gaLlBMu4BoueYG1dg2am0VICaaZRVlZGmAkbmEYTbDONYBnXYNETjNdVzTQU9dIctzqF/xAs4xoseoK1dfW7jXCzKCho1Pj6Bddee21A6NEU1ToGg64QPOMaLHqC8bo2py+1PNUEQogeVB3wUygUigCnp5TycGMVlNFogqqT5d2BE/W8vQUY0cKmW3ptS69rhdv49aR+XfTq14z/kS+6miGvL9cGi65m3b9mXGuWrq2ADNmEUVDLU01Q9Q+s1/IKIVxNbRo1REuv9eG66l9PGNyvGf+j6l+brasZ8vpybbDoatb9a8a1Jurq1TVqI9w3JplwrS99+oIZ8pqhq1nyKl31RemqEWp5Kkiodh3GC5c6f0fpGngEi55gfV3VTCN4KAPeolZcrQBG6Rp4BIueYHFd1UxDoVAoFF6jZhoKhUKh8BplNBQKhULhNcpoKBQKhcJrlNFQKBQKhdcooxFkCCH6CiG+EkLsF0KUCCFShBBvCSEcZsumBUKIp6t0KxVCxAshLjNbJq0RQrwqhIgTQpwQQhwVQvwghDjdbLmMoEp3KYT4yGxZ9EAI0UMIMU0IcUwIUSyESBBCDDNbrtoooxF8nIF73J8AhgAvAE8C/zVTKC0QQtwJfAS8DZwHrAWWCiF6mymXDlyB+wDXhcA1uCM7rBBCRJkqlc4IIc4HHgcSzZZFD4QQ7YD1QAVwA3Am8CKQZ6JYp6BcbhUIIf4OPCWl7G+2LL4ghNgMbJNSPlWrbCfwg5TyVfMk0xchRCfgKHCFlHKN2fLogRAiGtgGPI07jXOClPJ5U4XSGCHEOOASKaWlZ8dqpqEAaAPkmi2EL1Qtrw0DVtR5awVwsfESGUqbqp9+PYZNMAlYLKX8yWxBdGQksFUIMbdq2fFXIcRjZgtVF2U0ghwhxADgWeAzs2XxkY6AHciqU54FdDVeHGOoisL8AbBOSrnDbHn0QAhxF/A7IGBni1X0B54C9gLX4X4mPxZC3G+qVHVQRiNAEEL8q2qDsLHX8DrXdAeWAXOllF+aI7nm1F1vFfWUBRITgaHA3WYLogdCiF7AeOA+KWWp2fLojA338uo/pJS/Sik/BybjNiSWQYVGDxwmArOaqHOg+pcqg/EzsBH35qK/kwM4OXVW0ZlTZx8BgRBiAu4ljcullIGaKGwY7jGMrxUy3A5cLoR4BgiTUjrNEk5jjgDJdcp2AreZIEuDKKMRIEgpc3B/cDZJVTbCn4F44CEppUtP2YxASlkuhIjH7U0UU+uta4AF5kilD1VLUhOAW4ArpZT7TRZJT1YCZ9cp+wbYBbwTQAYD3J5TdV2nBwEHTZClQZTRCDKqZhirgTTgJaBT9Tc4KWWmeZJpwgfAVCHEVk7OoHrj//s1dZkE3APcDJwQQlTPrvKllCXmiaU9UsoTgMdejRCiCDgWgHs4HwIbhBD/AObgzr73OBZbCVBGI/i4FhhY9aq7pCFOre4/SClnCyE6AG8C3XB/2NwopbTUNzUNqF7jXl2n/CFgiqGSKDRDShknhLgFGIv7Ht4PPC+lnG6uZJ6ocxoKhUKh8BrlPaVQKBQKr1FGQ6FQKBReo4yGQqFQKLxGGQ2FQqFQeI0yGgqFQqHwGmU0FAqFQuE1ymgoFAqFwmuU0VAoFAqF1yijoVAoFAqvUUZDEbQIIToLIbqZLMNfdGz770KIm/RqXxGcKKOhCEqEEH1xB/4rNFGGi4EHdOziA+ApIcT1OvahCDJUwEJF0CGEaAVMA+6qiqJqFvcCugWjk1I6hRCjgLVCiF+llAGZV0RhLGqmoQhGPgC+NTNxkRAiBLgRnXN9SClzcSfnek/PfhTBgzIaiqBCCHEecBXuRD5mcj3uvN7FBvQ1EfiTEGKgAX0pAhxlNBTBxsvAN1LKSpPl0HVpqjZSyuPAKuARI/pTBDYqn4YiaBBCtMadL/xiKeWvJsoRDSQBpxllvIQQfwWeklKe0czregFPAg6gLfCTlHK29hIq/AW1Ea4IJq4GnEBi3TeEEOfgTn/bDxiHO1/zP4EwoCvwlpQyQSM5bgV+rGswhBCX4k7teRowBliKO0vfGUAocBbwkpRykxDiHuCyqkvPBt6UUq5qpM81wEdCiC7ebogLIZ4Cngdul1ImCSEigJVCiKFSyteq6lwIVEgp471pUxEASCnVS72C4gX8D/i5gfe+xv0l6nXgKPA90AP3/kc5MEFDOZYDF9QpswHf4U65+78qGSYAV9Wq8zmQijsV6K21yv8NHANsjfTZDpDAtV7K+AxQAQyrU/5n3Ib39FoydTZ7bNXLuJfa01AEE4OBw3ULhRADgAzp/ubfHWgP/FdKeRjogPssx4/1XCeEEPcIIX4UQnwihPhICPF8YwIIIboAfaSUm+u8dT7wq5RSVsnQCVgspfy5Vp0C3DOhLCnl/FrlWVUyd2qk6zzcH/Z9G5OvSsYzcHtbfSVPnUFswW3gbhVC2IEOUsqjTbWpCBzU8pQimOgB7K+nvAswr+r3S3F7NW0DkFJ+j3vW4YEQIhKYDXQEbpNSZlSV3yaEOFNKmdyADHdVXVeXMCCm6vfLgOVSymV16gytkv+LOuWDgRLcs416kVJKIUQ+0KahOrV4EQgHPqvnveyqn32Bh4GpXrSnCCCU0VAEE9FAft1CKeUGACFER9z7Bv9urJGqb9hLcS/5XArcWLW2L4H/0vgp83uBUfXIsKaq7YFAT+CjOn2GABcDc6tmI7W5Flgrm95UL8NtnJridiBT1rOHI6WsEEKAewZ2K+6zJoogQi1PKYIJJ24voIa4Cveewuom2nkJOA+4CfgD7i/yzwMzgYlSypL6LhJCDMLtsbi7kbZ/X/Wzrgzn4zZ6HuVCiLNxb5yfMhuqh/Y0ETZFCNEBt5fUtibauhH4Rz0GTBHgKKOhCCbycH9wNsRVuL+Nb2qoQtWy1GvAv6WUacA5ssoFVUq5BejfSPvenM24qkrOui7BV1b9XF2n/J4qmb+vkq/esxhCiCjcs4xT9nTqUFH1s7F9Conb+8s0t2WFeajlKUUwcYCmjcYmKWVpI3Vuwf0h/XHV3zWxq6qWkCrqu6iKOzn54d8QVwJrpJSuemRLrTJUtbkVWCilPF416+jZQLtdqn42NstBSlkghEjAvVdzCkKIx3F7k9mq/g6VUjamsyLAUDMNRTCRiNv76BSEEF1xn4dY3UQbI4DpUsryqr+LqowFwDW4T17X1/6FQJqUMrOhhoUQg3GfCVldpzwU935GfbJ1BH4RQtiAV3CHDKmPi3B7X+1oqP9avABcIYQ4t5YMw4UQX+E2kpOBQcK9uTGj6vyGIkhQMw1FMLEG+LcQoq2UMq/Oe11wu67OaaKNnrj3Lqr5stYG9LPAYw1c583SVAcgA5hfp7wtUIQ7Mm9dngMexT1D+UxK2ZAH1eXAknpmMKcgpVwthPg98JwQwol7LygBeF5KeUIIsQh4H/e5kvEN7eEoAhMVRkQRNFR9G88AHpVSLmphG1/jjl21tk75E0BrKeUp0WSrZiL7gLOlCaHYq/TeBzwjpVxidP+KwEItTymChqpv2V/iPivRUmJwL/UANQf8RgP96zMYVVwLbDbDYFRxO+6ZylKT+lcEEGqmoQgqhBCdgJ3A+VLK+g76edPG33AffnPh3liPlVLGNlJ/OjBbSrmwJf35ihBiG+6YVY3FplIovEIZDUXQIYR4EhghpXzYgL6igGRgoBleRkKIV3GH+njJ6L4VgYlanlIEHVLKzwC7EELP/NzV3IJ7A9oMg3EdcDpuryqFQhOU0VAEK48AF1a5uerJFbi9jMzgCPCQlNJpUv+KAEQtTykUCoXCa9RMQ6FQKBReo4yGQqFQKLxGGQ2FQqFQeI0yGgqFQqHwGmU0FAqFQuE1ymgoFAqFwmuU0VAoFAqF1yijoVAoFAqvUUZDoVAoFF6jjIZCoVAovEYZDYVCoVB4zf8DhiogMlYYLNoAAAAASUVORK5CYII=\n", "text/plain": [ "Graphics object consisting of 9 graphics primitives" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "shadow_plot(1, pi/2, orientation=pi/6, fill=True, color='red', \n", " number_colors=1, thickness=1.5, linestyle=':', \n", " draw_spin=True, spin_arrow_options={'color': 'blue', 'width': 3},\n", " legend=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Half-width of the image in units of $m/r_{\\mathscr{O}}$:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "fsize = 17.3275984058939 m/r_obs\n" ] } ], "source": [ "EHT_field_of_view = 127 # Field of view for EHT image of M87* (in microarcseconds)\n", "scale_M87 = 3.66467403690525 # m/r for M87* (in microarcseconds)\n", "fsize = EHT_field_of_view / 2 / scale_M87\n", "extent = (-fsize, fsize, -fsize, fsize)\n", "print(\"fsize =\", fsize, \"m/r_obs\")" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "def empty_plot(size, frame=True, axes=False, axes_labels='automatic', \n", " gridlines=False):\n", " g = Graphics()\n", " g._extra_kwds['frame'] = frame\n", " g._extra_kwds['axes'] = axes\n", " g._extra_kwds['gridlines'] = gridlines\n", " if axes_labels:\n", " if axes_labels == 'automatic':\n", " g.axes_labels([r\"$(r_{\\mathscr{O}}/m)\\; \\alpha$\", \n", " r\"$(r_{\\mathscr{O}}/m)\\; \\beta$\"])\n", " else:\n", " g.axes_labels(axes_labels)\n", " g.set_aspect_ratio(1)\n", " g.set_axes_range(-size, size, -size, size)\n", " return g" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Default resolution of SageMath images:" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}100$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}100\n", "\\end{math}" ], "text/plain": [ "100" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dpi = Graphics.SHOW_OPTIONS['dpi']\n", "dpi" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Default resolution of Matplotlib images:" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}72.0$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}72.0\n", "\\end{math}" ], "text/plain": [ "72.0" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from matplotlib import rcParams\n", "rcParams['figure.dpi']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The inclination of the jet axis with respect to the line of sight is $\\iota = 17^\\circ$ and its orientation\n", "is $\\Theta = -72^\\circ$ from the North axis. We have $\\theta_{\\mathscr{O}} = \\iota$ or $\\theta_{\\mathscr{O}} = \\pi - \\iota$. Assuming that the emitting matted in the EHT image rotates in the same direction as the black hole implies $\\theta_{\\mathscr{O}} = \\pi - \\iota$:" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "inc_jet = 17. # inclination of jet axis w.r.t. line of sight (in degrees)\n", "Theta_jet = -72. # orientation of jet axis in the plane of sky (in degress)\n", "th_obs = (180. - inc_jet)/180.*pi\n", "orient = (180. + Theta_jet)/180.*pi\n", "image_parameters = [(0, th_obs, orient, 'magenta'), (1., th_obs, orient, 'green')]" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "gimage = 'M87_EHT'" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "M87_EHT resolution: (474, 474, 4)\n", "figsize: (4.74, 4.74)\n", "a/m = 0 theta_obs = 0.905555555555556*pi\n", "a/m = 1.00000000000000 theta_obs = 0.905555555555556*pi\n", "rmin : 1.89727388775268 rmax : 2.90008539695632\n", "rmin_col : 1 rmax_col : 4.0\n" ] }, { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "/* global mpl */\n", "window.mpl = {};\n", "\n", "mpl.get_websocket_type = function () {\n", " if (typeof WebSocket !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof MozWebSocket !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert(\n", " 'Your browser does not have WebSocket support. ' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.'\n", " );\n", " }\n", "};\n", "\n", "mpl.figure = function (figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = this.ws.binaryType !== undefined;\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById('mpl-warnings');\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent =\n", " 'This browser does not support binary websocket messages. ' +\n", " 'Performance may be slow.';\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = document.createElement('div');\n", " this.root.setAttribute('style', 'display: inline-block');\n", " this._root_extra_style(this.root);\n", "\n", " parent_element.appendChild(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message('supports_binary', { value: fig.supports_binary });\n", " fig.send_message('send_image_mode', {});\n", " if (fig.ratio !== 1) {\n", " fig.send_message('set_dpi_ratio', { dpi_ratio: fig.ratio });\n", " }\n", " fig.send_message('refresh', {});\n", " };\n", "\n", " this.imageObj.onload = function () {\n", " if (fig.image_mode === 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " };\n", "\n", " this.imageObj.onunload = function () {\n", " fig.ws.close();\n", " };\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "};\n", "\n", "mpl.figure.prototype._init_header = function () {\n", " var titlebar = document.createElement('div');\n", " titlebar.classList =\n", " 'ui-dialog-titlebar ui-widget-header ui-corner-all ui-helper-clearfix';\n", " var titletext = document.createElement('div');\n", " titletext.classList = 'ui-dialog-title';\n", " titletext.setAttribute(\n", " 'style',\n", " 'width: 100%; text-align: center; padding: 3px;'\n", " );\n", " titlebar.appendChild(titletext);\n", " this.root.appendChild(titlebar);\n", " this.header = titletext;\n", "};\n", "\n", "mpl.figure.prototype._canvas_extra_style = function (_canvas_div) {};\n", "\n", "mpl.figure.prototype._root_extra_style = function (_canvas_div) {};\n", "\n", "mpl.figure.prototype._init_canvas = function () {\n", " var fig = this;\n", "\n", " var canvas_div = (this.canvas_div = document.createElement('div'));\n", " canvas_div.setAttribute(\n", " 'style',\n", " 'border: 1px solid #ddd;' +\n", " 'box-sizing: content-box;' +\n", " 'clear: both;' +\n", " 'min-height: 1px;' +\n", " 'min-width: 1px;' +\n", " 'outline: 0;' +\n", " 'overflow: hidden;' +\n", " 'position: relative;' +\n", " 'resize: both;'\n", " );\n", "\n", " function on_keyboard_event_closure(name) {\n", " return function (event) {\n", " return fig.key_event(event, name);\n", " };\n", " }\n", "\n", " canvas_div.addEventListener(\n", " 'keydown',\n", " on_keyboard_event_closure('key_press')\n", " );\n", " canvas_div.addEventListener(\n", " 'keyup',\n", " on_keyboard_event_closure('key_release')\n", " );\n", "\n", " this._canvas_extra_style(canvas_div);\n", " this.root.appendChild(canvas_div);\n", "\n", " var canvas = (this.canvas = document.createElement('canvas'));\n", " canvas.classList.add('mpl-canvas');\n", " canvas.setAttribute('style', 'box-sizing: content-box;');\n", "\n", " this.context = canvas.getContext('2d');\n", "\n", " var backingStore =\n", " this.context.backingStorePixelRatio ||\n", " this.context.webkitBackingStorePixelRatio ||\n", " this.context.mozBackingStorePixelRatio ||\n", " this.context.msBackingStorePixelRatio ||\n", " this.context.oBackingStorePixelRatio ||\n", " this.context.backingStorePixelRatio ||\n", " 1;\n", "\n", " this.ratio = (window.devicePixelRatio || 1) / backingStore;\n", " if (this.ratio !== 1) {\n", " fig.send_message('set_dpi_ratio', { dpi_ratio: this.ratio });\n", " }\n", "\n", " var rubberband_canvas = (this.rubberband_canvas = document.createElement(\n", " 'canvas'\n", " ));\n", " rubberband_canvas.setAttribute(\n", " 'style',\n", " 'box-sizing: content-box; position: absolute; left: 0; top: 0; z-index: 1;'\n", " );\n", "\n", " var resizeObserver = new ResizeObserver(function (entries) {\n", " var nentries = entries.length;\n", " for (var i = 0; i < nentries; i++) {\n", " var entry = entries[i];\n", " var width, height;\n", " if (entry.contentBoxSize) {\n", " if (entry.contentBoxSize instanceof Array) {\n", " // Chrome 84 implements new version of spec.\n", " width = entry.contentBoxSize[0].inlineSize;\n", " height = entry.contentBoxSize[0].blockSize;\n", " } else {\n", " // Firefox implements old version of spec.\n", " width = entry.contentBoxSize.inlineSize;\n", " height = entry.contentBoxSize.blockSize;\n", " }\n", " } else {\n", " // Chrome <84 implements even older version of spec.\n", " width = entry.contentRect.width;\n", " height = entry.contentRect.height;\n", " }\n", "\n", " // Keep the size of the canvas and rubber band canvas in sync with\n", " // the canvas container.\n", " if (entry.devicePixelContentBoxSize) {\n", " // Chrome 84 implements new version of spec.\n", " canvas.setAttribute(\n", " 'width',\n", " entry.devicePixelContentBoxSize[0].inlineSize\n", " );\n", " canvas.setAttribute(\n", " 'height',\n", " entry.devicePixelContentBoxSize[0].blockSize\n", " );\n", " } else {\n", " canvas.setAttribute('width', width * fig.ratio);\n", " canvas.setAttribute('height', height * fig.ratio);\n", " }\n", " canvas.setAttribute(\n", " 'style',\n", " 'width: ' + width + 'px; height: ' + height + 'px;'\n", " );\n", "\n", " rubberband_canvas.setAttribute('width', width);\n", " rubberband_canvas.setAttribute('height', height);\n", "\n", " // And update the size in Python. We ignore the initial 0/0 size\n", " // that occurs as the element is placed into the DOM, which should\n", " // otherwise not happen due to the minimum size styling.\n", " if (width != 0 && height != 0) {\n", " fig.request_resize(width, height);\n", " }\n", " }\n", " });\n", " resizeObserver.observe(canvas_div);\n", "\n", " function on_mouse_event_closure(name) {\n", " return function (event) {\n", " return fig.mouse_event(event, name);\n", " };\n", " }\n", "\n", " rubberband_canvas.addEventListener(\n", " 'mousedown',\n", " on_mouse_event_closure('button_press')\n", " );\n", " rubberband_canvas.addEventListener(\n", " 'mouseup',\n", " on_mouse_event_closure('button_release')\n", " );\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband_canvas.addEventListener(\n", " 'mousemove',\n", " on_mouse_event_closure('motion_notify')\n", " );\n", "\n", " rubberband_canvas.addEventListener(\n", " 'mouseenter',\n", " on_mouse_event_closure('figure_enter')\n", " );\n", " rubberband_canvas.addEventListener(\n", " 'mouseleave',\n", " on_mouse_event_closure('figure_leave')\n", " );\n", "\n", " canvas_div.addEventListener('wheel', function (event) {\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " on_mouse_event_closure('scroll')(event);\n", " });\n", "\n", " canvas_div.appendChild(canvas);\n", " canvas_div.appendChild(rubberband_canvas);\n", "\n", " this.rubberband_context = rubberband_canvas.getContext('2d');\n", " this.rubberband_context.strokeStyle = '#000000';\n", "\n", " this._resize_canvas = function (width, height, forward) {\n", " if (forward) {\n", " canvas_div.style.width = width + 'px';\n", " canvas_div.style.height = height + 'px';\n", " }\n", " };\n", "\n", " // Disable right mouse context menu.\n", " this.rubberband_canvas.addEventListener('contextmenu', function (_e) {\n", " event.preventDefault();\n", " return false;\n", " });\n", "\n", " function set_focus() {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "};\n", "\n", "mpl.figure.prototype._init_toolbar = function () {\n", " var fig = this;\n", "\n", " var toolbar = document.createElement('div');\n", " toolbar.classList = 'mpl-toolbar';\n", " this.root.appendChild(toolbar);\n", "\n", " function on_click_closure(name) {\n", " return function (_event) {\n", " return fig.toolbar_button_onclick(name);\n", " };\n", " }\n", "\n", " function on_mouseover_closure(tooltip) {\n", " return function (event) {\n", " if (!event.currentTarget.disabled) {\n", " return fig.toolbar_button_onmouseover(tooltip);\n", " }\n", " };\n", " }\n", "\n", " fig.buttons = {};\n", " var buttonGroup = document.createElement('div');\n", " buttonGroup.classList = 'mpl-button-group';\n", " for (var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " /* Instead of a spacer, we start a new button group. */\n", " if (buttonGroup.hasChildNodes()) {\n", " toolbar.appendChild(buttonGroup);\n", " }\n", " buttonGroup = document.createElement('div');\n", " buttonGroup.classList = 'mpl-button-group';\n", " continue;\n", " }\n", "\n", " var button = (fig.buttons[name] = document.createElement('button'));\n", " button.classList = 'mpl-widget';\n", " button.setAttribute('role', 'button');\n", " button.setAttribute('aria-disabled', 'false');\n", " button.addEventListener('click', on_click_closure(method_name));\n", " button.addEventListener('mouseover', on_mouseover_closure(tooltip));\n", "\n", " var icon_img = document.createElement('img');\n", " icon_img.src = '_images/' + image + '.png';\n", " icon_img.srcset = '_images/' + image + '_large.png 2x';\n", " icon_img.alt = tooltip;\n", " button.appendChild(icon_img);\n", "\n", " buttonGroup.appendChild(button);\n", " }\n", "\n", " if (buttonGroup.hasChildNodes()) {\n", " toolbar.appendChild(buttonGroup);\n", " }\n", "\n", " var fmt_picker = document.createElement('select');\n", " fmt_picker.classList = 'mpl-widget';\n", " toolbar.appendChild(fmt_picker);\n", " this.format_dropdown = fmt_picker;\n", "\n", " for (var ind in mpl.extensions) {\n", " var fmt = mpl.extensions[ind];\n", " var option = document.createElement('option');\n", " option.selected = fmt === mpl.default_extension;\n", " option.innerHTML = fmt;\n", " fmt_picker.appendChild(option);\n", " }\n", "\n", " var status_bar = document.createElement('span');\n", " status_bar.classList = 'mpl-message';\n", " toolbar.appendChild(status_bar);\n", " this.message = status_bar;\n", "};\n", "\n", "mpl.figure.prototype.request_resize = function (x_pixels, y_pixels) {\n", " // Request matplotlib to resize the figure. Matplotlib will then trigger a resize in the client,\n", " // which will in turn request a refresh of the image.\n", " this.send_message('resize', { width: x_pixels, height: y_pixels });\n", "};\n", "\n", "mpl.figure.prototype.send_message = function (type, properties) {\n", " properties['type'] = type;\n", " properties['figure_id'] = this.id;\n", " this.ws.send(JSON.stringify(properties));\n", "};\n", "\n", "mpl.figure.prototype.send_draw_message = function () {\n", " if (!this.waiting) {\n", " this.waiting = true;\n", " this.ws.send(JSON.stringify({ type: 'draw', figure_id: this.id }));\n", " }\n", "};\n", "\n", "mpl.figure.prototype.handle_save = function (fig, _msg) {\n", " var format_dropdown = fig.format_dropdown;\n", " var format = format_dropdown.options[format_dropdown.selectedIndex].value;\n", " fig.ondownload(fig, format);\n", "};\n", "\n", "mpl.figure.prototype.handle_resize = function (fig, msg) {\n", " var size = msg['size'];\n", " if (size[0] !== fig.canvas.width || size[1] !== fig.canvas.height) {\n", " fig._resize_canvas(size[0], size[1], msg['forward']);\n", " fig.send_message('refresh', {});\n", " }\n", "};\n", "\n", "mpl.figure.prototype.handle_rubberband = function (fig, msg) {\n", " var x0 = msg['x0'] / fig.ratio;\n", " var y0 = (fig.canvas.height - msg['y0']) / fig.ratio;\n", " var x1 = msg['x1'] / fig.ratio;\n", " var y1 = (fig.canvas.height - msg['y1']) / fig.ratio;\n", " x0 = Math.floor(x0) + 0.5;\n", " y0 = Math.floor(y0) + 0.5;\n", " x1 = Math.floor(x1) + 0.5;\n", " y1 = Math.floor(y1) + 0.5;\n", " var min_x = Math.min(x0, x1);\n", " var min_y = Math.min(y0, y1);\n", " var width = Math.abs(x1 - x0);\n", " var height = Math.abs(y1 - y0);\n", "\n", " fig.rubberband_context.clearRect(\n", " 0,\n", " 0,\n", " fig.canvas.width / fig.ratio,\n", " fig.canvas.height / fig.ratio\n", " );\n", "\n", " fig.rubberband_context.strokeRect(min_x, min_y, width, height);\n", "};\n", "\n", "mpl.figure.prototype.handle_figure_label = function (fig, msg) {\n", " // Updates the figure title.\n", " fig.header.textContent = msg['label'];\n", "};\n", "\n", "mpl.figure.prototype.handle_cursor = function (fig, msg) {\n", " var cursor = msg['cursor'];\n", " switch (cursor) {\n", " case 0:\n", " cursor = 'pointer';\n", " break;\n", " case 1:\n", " cursor = 'default';\n", " break;\n", " case 2:\n", " cursor = 'crosshair';\n", " break;\n", " case 3:\n", " cursor = 'move';\n", " break;\n", " }\n", " fig.rubberband_canvas.style.cursor = cursor;\n", "};\n", "\n", "mpl.figure.prototype.handle_message = function (fig, msg) {\n", " fig.message.textContent = msg['message'];\n", "};\n", "\n", "mpl.figure.prototype.handle_draw = function (fig, _msg) {\n", " // Request the server to send over a new figure.\n", " fig.send_draw_message();\n", "};\n", "\n", "mpl.figure.prototype.handle_image_mode = function (fig, msg) {\n", " fig.image_mode = msg['mode'];\n", "};\n", "\n", "mpl.figure.prototype.handle_history_buttons = function (fig, msg) {\n", " for (var key in msg) {\n", " if (!(key in fig.buttons)) {\n", " continue;\n", " }\n", " fig.buttons[key].disabled = !msg[key];\n", " fig.buttons[key].setAttribute('aria-disabled', !msg[key]);\n", " }\n", "};\n", "\n", "mpl.figure.prototype.handle_navigate_mode = function (fig, msg) {\n", " if (msg['mode'] === 'PAN') {\n", " fig.buttons['Pan'].classList.add('active');\n", " fig.buttons['Zoom'].classList.remove('active');\n", " } else if (msg['mode'] === 'ZOOM') {\n", " fig.buttons['Pan'].classList.remove('active');\n", " fig.buttons['Zoom'].classList.add('active');\n", " } else {\n", " fig.buttons['Pan'].classList.remove('active');\n", " fig.buttons['Zoom'].classList.remove('active');\n", " }\n", "};\n", "\n", "mpl.figure.prototype.updated_canvas_event = function () {\n", " // Called whenever the canvas gets updated.\n", " this.send_message('ack', {});\n", "};\n", "\n", "// A function to construct a web socket function for onmessage handling.\n", "// Called in the figure constructor.\n", "mpl.figure.prototype._make_on_message_function = function (fig) {\n", " return function socket_on_message(evt) {\n", " if (evt.data instanceof Blob) {\n", " /* FIXME: We get \"Resource interpreted as Image but\n", " * transferred with MIME type text/plain:\" errors on\n", " * Chrome. But how to set the MIME type? It doesn't seem\n", " * to be part of the websocket stream */\n", " evt.data.type = 'image/png';\n", "\n", " /* Free the memory for the previous frames */\n", " if (fig.imageObj.src) {\n", " (window.URL || window.webkitURL).revokeObjectURL(\n", " fig.imageObj.src\n", " );\n", " }\n", "\n", " fig.imageObj.src = (window.URL || window.webkitURL).createObjectURL(\n", " evt.data\n", " );\n", " fig.updated_canvas_event();\n", " fig.waiting = false;\n", " return;\n", " } else if (\n", " typeof evt.data === 'string' &&\n", " evt.data.slice(0, 21) === 'data:image/png;base64'\n", " ) {\n", " fig.imageObj.src = evt.data;\n", " fig.updated_canvas_event();\n", " fig.waiting = false;\n", " return;\n", " }\n", "\n", " var msg = JSON.parse(evt.data);\n", " var msg_type = msg['type'];\n", "\n", " // Call the \"handle_{type}\" callback, which takes\n", " // the figure and JSON message as its only arguments.\n", " try {\n", " var callback = fig['handle_' + msg_type];\n", " } catch (e) {\n", " console.log(\n", " \"No handler for the '\" + msg_type + \"' message type: \",\n", " msg\n", " );\n", " return;\n", " }\n", "\n", " if (callback) {\n", " try {\n", " // console.log(\"Handling '\" + msg_type + \"' message: \", msg);\n", " callback(fig, msg);\n", " } catch (e) {\n", " console.log(\n", " \"Exception inside the 'handler_\" + msg_type + \"' callback:\",\n", " e,\n", " e.stack,\n", " msg\n", " );\n", " }\n", " }\n", " };\n", "};\n", "\n", "// from http://stackoverflow.com/questions/1114465/getting-mouse-location-in-canvas\n", "mpl.findpos = function (e) {\n", " //this section is from http://www.quirksmode.org/js/events_properties.html\n", " var targ;\n", " if (!e) {\n", " e = window.event;\n", " }\n", " if (e.target) {\n", " targ = e.target;\n", " } else if (e.srcElement) {\n", " targ = e.srcElement;\n", " }\n", " if (targ.nodeType === 3) {\n", " // defeat Safari bug\n", " targ = targ.parentNode;\n", " }\n", "\n", " // pageX,Y are the mouse positions relative to the document\n", " var boundingRect = targ.getBoundingClientRect();\n", " var x = e.pageX - (boundingRect.left + document.body.scrollLeft);\n", " var y = e.pageY - (boundingRect.top + document.body.scrollTop);\n", "\n", " return { x: x, y: y };\n", "};\n", "\n", "/*\n", " * return a copy of an object with only non-object keys\n", " * we need this to avoid circular references\n", " * http://stackoverflow.com/a/24161582/3208463\n", " */\n", "function simpleKeys(original) {\n", " return Object.keys(original).reduce(function (obj, key) {\n", " if (typeof original[key] !== 'object') {\n", " obj[key] = original[key];\n", " }\n", " return obj;\n", " }, {});\n", "}\n", "\n", "mpl.figure.prototype.mouse_event = function (event, name) {\n", " var canvas_pos = mpl.findpos(event);\n", "\n", " if (name === 'button_press') {\n", " this.canvas.focus();\n", " this.canvas_div.focus();\n", " }\n", "\n", " var x = canvas_pos.x * this.ratio;\n", " var y = canvas_pos.y * this.ratio;\n", "\n", " this.send_message(name, {\n", " x: x,\n", " y: y,\n", " button: event.button,\n", " step: event.step,\n", " guiEvent: simpleKeys(event),\n", " });\n", "\n", " /* This prevents the web browser from automatically changing to\n", " * the text insertion cursor when the button is pressed. We want\n", " * to control all of the cursor setting manually through the\n", " * 'cursor' event from matplotlib */\n", " event.preventDefault();\n", " return false;\n", "};\n", "\n", "mpl.figure.prototype._key_event_extra = function (_event, _name) {\n", " // Handle any extra behaviour associated with a key event\n", "};\n", "\n", "mpl.figure.prototype.key_event = function (event, name) {\n", " // Prevent repeat events\n", " if (name === 'key_press') {\n", " if (event.which === this._key) {\n", " return;\n", " } else {\n", " this._key = event.which;\n", " }\n", " }\n", " if (name === 'key_release') {\n", " this._key = null;\n", " }\n", "\n", " var value = '';\n", " if (event.ctrlKey && event.which !== 17) {\n", " value += 'ctrl+';\n", " }\n", " if (event.altKey && event.which !== 18) {\n", " value += 'alt+';\n", " }\n", " if (event.shiftKey && event.which !== 16) {\n", " value += 'shift+';\n", " }\n", "\n", " value += 'k';\n", " value += event.which.toString();\n", "\n", " this._key_event_extra(event, name);\n", "\n", " this.send_message(name, { key: value, guiEvent: simpleKeys(event) });\n", " return false;\n", "};\n", "\n", "mpl.figure.prototype.toolbar_button_onclick = function (name) {\n", " if (name === 'download') {\n", " this.handle_save(this, null);\n", " } else {\n", " this.send_message('toolbar_button', { name: name });\n", " }\n", "};\n", "\n", "mpl.figure.prototype.toolbar_button_onmouseover = function (tooltip) {\n", " this.message.textContent = tooltip;\n", "};\n", "mpl.toolbar_items = [[\"Home\", \"Reset original view\", \"fa fa-home icon-home\", \"home\"], [\"Back\", \"Back to previous view\", \"fa fa-arrow-left icon-arrow-left\", \"back\"], [\"Forward\", \"Forward to next view\", \"fa fa-arrow-right icon-arrow-right\", \"forward\"], [\"\", \"\", \"\", \"\"], [\"Pan\", \"Left button pans, Right button zooms\\nx/y fixes axis, CTRL fixes aspect\", \"fa fa-arrows icon-move\", \"pan\"], [\"Zoom\", \"Zoom to rectangle\\nx/y fixes axis, CTRL fixes aspect\", \"fa fa-square-o icon-check-empty\", \"zoom\"], [\"\", \"\", \"\", \"\"], [\"Download\", \"Download plot\", \"fa fa-floppy-o icon-save\", \"download\"]];\n", "\n", "mpl.extensions = [\"eps\", \"jpeg\", \"pdf\", \"png\", \"ps\", \"raw\", \"svg\", \"tif\"];\n", "\n", "mpl.default_extension = \"png\";/* global mpl */\n", "\n", "var comm_websocket_adapter = function (comm) {\n", " // Create a \"websocket\"-like object which calls the given IPython comm\n", " // object with the appropriate methods. Currently this is a non binary\n", " // socket, so there is still some room for performance tuning.\n", " var ws = {};\n", "\n", " ws.close = function () {\n", " comm.close();\n", " };\n", " ws.send = function (m) {\n", " //console.log('sending', m);\n", " comm.send(m);\n", " };\n", " // Register the callback with on_msg.\n", " comm.on_msg(function (msg) {\n", " //console.log('receiving', msg['content']['data'], msg);\n", " // Pass the mpl event to the overridden (by mpl) onmessage function.\n", " ws.onmessage(msg['content']['data']);\n", " });\n", " return ws;\n", "};\n", "\n", "mpl.mpl_figure_comm = function (comm, msg) {\n", " // This is the function which gets called when the mpl process\n", " // starts-up an IPython Comm through the \"matplotlib\" channel.\n", "\n", " var id = msg.content.data.id;\n", " // Get hold of the div created by the display call when the Comm\n", " // socket was opened in Python.\n", " var element = document.getElementById(id);\n", " var ws_proxy = comm_websocket_adapter(comm);\n", "\n", " function ondownload(figure, _format) {\n", " window.open(figure.canvas.toDataURL());\n", " }\n", "\n", " var fig = new mpl.figure(id, ws_proxy, ondownload, element);\n", "\n", " // Call onopen now - mpl needs it, as it is assuming we've passed it a real\n", " // web socket which is closed, not our websocket->open comm proxy.\n", " ws_proxy.onopen();\n", "\n", " fig.parent_element = element;\n", " fig.cell_info = mpl.find_output_cell(\"

\");\n", " if (!fig.cell_info) {\n", " console.error('Failed to find cell for figure', id, fig);\n", " return;\n", " }\n", " fig.cell_info[0].output_area.element.one(\n", " 'cleared',\n", " { fig: fig },\n", " fig._remove_fig_handler\n", " );\n", "};\n", "\n", "mpl.figure.prototype.handle_close = function (fig, msg) {\n", " var width = fig.canvas.width / fig.ratio;\n", " fig.cell_info[0].output_area.element.off(\n", " 'cleared',\n", " fig._remove_fig_handler\n", " );\n", "\n", " // Update the output cell to use the data from the current canvas.\n", " fig.push_to_output();\n", " var dataURL = fig.canvas.toDataURL();\n", " // Re-enable the keyboard manager in IPython - without this line, in FF,\n", " // the notebook keyboard shortcuts fail.\n", " IPython.keyboard_manager.enable();\n", " fig.parent_element.innerHTML =\n", " '

';\n", " fig.close_ws(fig, msg);\n", "};\n", "\n", "mpl.figure.prototype.close_ws = function (fig, msg) {\n", " fig.send_message('closing', msg);\n", " // fig.ws.close()\n", "};\n", "\n", "mpl.figure.prototype.push_to_output = function (_remove_interactive) {\n", " // Turn the data on the canvas into data in the output cell.\n", " var width = this.canvas.width / this.ratio;\n", " var dataURL = this.canvas.toDataURL();\n", " this.cell_info[1]['text/html'] =\n", " '

';\n", "};\n", "\n", "mpl.figure.prototype.updated_canvas_event = function () {\n", " // Tell IPython that the notebook contents must change.\n", " IPython.notebook.set_dirty(true);\n", " this.send_message('ack', {});\n", " var fig = this;\n", " // Wait a second, then push the new image to the DOM so\n", " // that it is saved nicely (might be nice to debounce this).\n", " setTimeout(function () {\n", " fig.push_to_output();\n", " }, 1000);\n", "};\n", "\n", "mpl.figure.prototype._init_toolbar = function () {\n", " var fig = this;\n", "\n", " var toolbar = document.createElement('div');\n", " toolbar.classList = 'btn-toolbar';\n", " this.root.appendChild(toolbar);\n", "\n", " function on_click_closure(name) {\n", " return function (_event) {\n", " return fig.toolbar_button_onclick(name);\n", " };\n", " }\n", "\n", " function on_mouseover_closure(tooltip) {\n", " return function (event) {\n", " if (!event.currentTarget.disabled) {\n", " return fig.toolbar_button_onmouseover(tooltip);\n", " }\n", " };\n", " }\n", "\n", " fig.buttons = {};\n", " var buttonGroup = document.createElement('div');\n", " buttonGroup.classList = 'btn-group';\n", " var button;\n", " for (var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " /* Instead of a spacer, we start a new button group. */\n", " if (buttonGroup.hasChildNodes()) {\n", " toolbar.appendChild(buttonGroup);\n", " }\n", " buttonGroup = document.createElement('div');\n", " buttonGroup.classList = 'btn-group';\n", " continue;\n", " }\n", "\n", " button = fig.buttons[name] = document.createElement('button');\n", " button.classList = 'btn btn-default';\n", " button.href = '#';\n", " button.title = name;\n", " button.innerHTML = '';\n", " button.addEventListener('click', on_click_closure(method_name));\n", " button.addEventListener('mouseover', on_mouseover_closure(tooltip));\n", " buttonGroup.appendChild(button);\n", " }\n", "\n", " if (buttonGroup.hasChildNodes()) {\n", " toolbar.appendChild(buttonGroup);\n", " }\n", "\n", " // Add the status bar.\n", " var status_bar = document.createElement('span');\n", " status_bar.classList = 'mpl-message pull-right';\n", " toolbar.appendChild(status_bar);\n", " this.message = status_bar;\n", "\n", " // Add the close button to the window.\n", " var buttongrp = document.createElement('div');\n", " buttongrp.classList = 'btn-group inline pull-right';\n", " button = document.createElement('button');\n", " button.classList = 'btn btn-mini btn-primary';\n", " button.href = '#';\n", " button.title = 'Stop Interaction';\n", " button.innerHTML = '';\n", " button.addEventListener('click', function (_evt) {\n", " fig.handle_close(fig, {});\n", " });\n", " button.addEventListener(\n", " 'mouseover',\n", " on_mouseover_closure('Stop Interaction')\n", " );\n", " buttongrp.appendChild(button);\n", " var titlebar = this.root.querySelector('.ui-dialog-titlebar');\n", " titlebar.insertBefore(buttongrp, titlebar.firstChild);\n", "};\n", "\n", "mpl.figure.prototype._remove_fig_handler = function (event) {\n", " var fig = event.data.fig;\n", " fig.close_ws(fig, {});\n", "};\n", "\n", "mpl.figure.prototype._root_extra_style = function (el) {\n", " el.style.boxSizing = 'content-box'; // override notebook setting of border-box.\n", "};\n", "\n", "mpl.figure.prototype._canvas_extra_style = function (el) {\n", " // this is important to make the div 'focusable\n", " el.setAttribute('tabindex', 0);\n", " // reach out to IPython and tell the keyboard manager to turn it's self\n", " // off when our div gets focus\n", "\n", " // location in version 3\n", " if (IPython.notebook.keyboard_manager) {\n", " IPython.notebook.keyboard_manager.register_events(el);\n", " } else {\n", " // location in version 2\n", " IPython.keyboard_manager.register_events(el);\n", " }\n", "};\n", "\n", "mpl.figure.prototype._key_event_extra = function (event, _name) {\n", " var manager = IPython.notebook.keyboard_manager;\n", " if (!manager) {\n", " manager = IPython.keyboard_manager;\n", " }\n", "\n", " // Check for shift+enter\n", " if (event.shiftKey && event.which === 13) {\n", " this.canvas_div.blur();\n", " // select the cell after this one\n", " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", " IPython.notebook.select(index + 1);\n", " }\n", "};\n", "\n", "mpl.figure.prototype.handle_save = function (fig, _msg) {\n", " fig.ondownload(fig, null);\n", "};\n", "\n", "mpl.find_output_cell = function (html_output) {\n", " // Return the cell and output element which can be found *uniquely* in the notebook.\n", " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", " // IPython event is triggered only after the cells have been serialised, which for\n", " // our purposes (turning an active figure into a static one), is too late.\n", " var cells = IPython.notebook.get_cells();\n", " var ncells = cells.length;\n", " for (var i = 0; i < ncells; i++) {\n", " var cell = cells[i];\n", " if (cell.cell_type === 'code') {\n", " for (var j = 0; j < cell.output_area.outputs.length; j++) {\n", " var data = cell.output_area.outputs[j];\n", " if (data.data) {\n", " // IPython >= 3 moved mimebundle to data attribute of output\n", " data = data.data;\n", " }\n", " if (data['text/html'] === html_output) {\n", " return [cell, data, j];\n", " }\n", " }\n", " }\n", " }\n", "};\n", "\n", "// Register the function which deals with the matplotlib target/channel.\n", "// The kernel may be null if the page has been refreshed.\n", "if (IPython.notebook.kernel !== null) {\n", " IPython.notebook.kernel.comm_manager.register_target(\n", " 'matplotlib',\n", " mpl.mpl_figure_comm\n", " );\n", "}\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "

" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%display plain\n", "%matplotlib notebook\n", "frame = True\n", "axes_labels = [r\"$(r_{\\mathscr{O}}/m)\\; \\alpha_{\\rm s}$\", \n", " r\"$(r_{\\mathscr{O}}/m)\\; \\beta_{\\rm s}$\"]\n", "test_size = False\n", "raw = False\n", "img = mpimg.imread(gimage + '.png')\n", "print(gimage, \" resolution:\", img.shape)\n", "# size of the image in inches:\n", "figsize = (float(img.shape[0])/dpi, float(img.shape[1])/dpi)\n", "print(\"figsize: \", figsize)\n", "gcrit = Graphics()\n", "for a0, th_obs, orient, color in image_parameters:\n", " print(\"a/m =\", a0, \" theta_obs =\", th_obs)\n", " # Critical curve as a SageMath graphics object from shadow_plot\n", " if th_obs == 0:\n", " th_obs = 0.001\n", " if not raw:\n", " g0 = shadow_plot(a0, th_obs, orientation=orient, fill=False, \n", " color=color, number_colors=1, thickness=2., \n", " linestyle=':', draw_spin=True,\n", " frame=frame, axes=False, \n", " axes_labels=axes_labels, gridlines=False, \n", " legend=False)\n", " gcrit += g0\n", " gcrit.axes_labels(g0.axes_labels())\n", " else:\n", " gcrit += empty_plot(fsize, axes_labels=axes_labels, frame=frame)\n", " if test_size:\n", " axes_labels_bck = gcrit.axes_labels()\n", " gcrit += point((-fsize, 0), size=60, color='red', zorder=100)\n", " gcrit += point((fsize, 0), size=60, color='red', zorder=100)\n", " gcrit += point((0, -fsize), size=60, color='red', zorder=100)\n", " gcrit += point((0, fsize), size=60, color='red', zorder=100)\n", " gcrit.axes_labels(axes_labels_bck)\n", "# Adding the jet axis\n", "ta = tan((90 + Theta_jet)/180.*pi)\n", "gcrit += arrow2d((7, 7*ta), (16, 16*ta), color='lightskyblue', \n", " linestyle='--', thickness=2)\n", "gcrit += text(\"jet\", (11, 11*ta+1.5), rotation=17, color='lightskyblue',\n", " fontsize=14)\n", "gcrit.axes_labels(g0.axes_labels())\n", "gcrit.set_axes_range(-fsize, fsize, -fsize, fsize)\n", "# Matplotlib figure corresponding to gcrit:\n", "options = gcrit.SHOW_OPTIONS.copy()\n", "options.update(gcrit._extra_kwds)\n", "options['figsize'] = figsize\n", "options['axes_pad'] = 0\n", "options.pop('dpi') # strip meaningless options for matplotlib\n", "options.pop('transparent') #\n", "options.pop('fig_tight') #\n", "fcrit = gcrit.matplotlib(**options)\n", "# Adding the M87* image onto it:\n", "ax = fcrit.axes[0]\n", "ax.imshow(img, extent=extent)\n", "# Save result to png and pdf\n", "ext = '_raw' if raw else '_crit'\n", "fcrit.savefig(gimage + ext + '.png', dpi=dpi, pad_inches=0, \n", " bbox_inches='tight', transparent=True)\n", "fcrit.savefig(gimage + ext + '.pdf', pad_inches=0, \n", " bbox_inches='tight')\n", "# Only for display in the current notebook:\n", "fig = plt.figure(figsize=figsize, dpi=dpi, frameon=False)\n", "imgc = mpimg.imread(gimage + '_crit.png')\n", "if axes_labels:\n", " imgcp = plt.imshow(imgc)\n", "else:\n", " imgcp = plt.imshow(imgc, extent=extent)" ] } ], "metadata": { "kernelspec": { "display_name": "SageMath 9.3.beta8", "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 }