Lines on the Hyperbolic Poincare Disk.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lines in the Poincaré Disk Model of Hyperbolic Geometry\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Someone I follow on social media made [this cool graphic](https://society6.com/product/a-hypergeometric-transformation_print#s6-3862893p4a1v45). I wondered how it would look if she used arced lines like in the Poincaré disk. She responded that doing computer graphics with hyperbolic geometry is kind of a pain in the butt. This workbook is to respond to that: affirmatively, because when I actually sat down and started working it out, it did become kind of a pain. It is possible, though I wouldn't want to try to do it in LaTeX (which is what she was using)<sup>1</sup>. Using complex numbers for the calculations made it a bit easier. \n",
    "\n",
    "Let's load `numpy` and `matplotlib` and do all the other initialization stuff."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "np.seterr(divide='ignore');#so it's not constantly yelling at us for dividing by zero"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's define the disk as $\\{z \\in \\mathbb{C} : \\left|z\\right| \\leq 1\\}$. Points on the unit circle $\\left|z\\right| = 1$ are called *ideal points* or *points at infinity*. The \"lines\" (geodesics) of the Poincaré disk are all circular arcs that meet the unit circle at right angles, i.e. all circles orthogonal to the unit circle, as well as all diameters of the disk. (It can be convenient to consider those lines to be circles with infinite center and radius.)\n",
    "\n",
    "Here's a basic task in geometry: draw a line between two points $a$ and $b$. On the Poincaré disk, that means we need a circle that is orthogonal to the unit circle and passes through $a$ and $b$. Most texts on hyperbolic geometry tell you how to do this, but they give a compass and straightedge construction, which isn't very useful if you're trying to plot it on a computer screen. You would rather have a formula for the (euclidean) center and the radius, and then you can tell your software, \"Hey, go plot a circle.\"\n",
    "\n",
    "You can translate the steps of the compass and straightedge construction into complex rational formulas. The formula for the center of the circle is fairly nice<sup>2</sup>:\n",
    "$$\n",
    "q = \\frac{a\\left(1+\\left|b\\right|^2\\right)-b\\left(1+\\left|a\\right|^2\\right)}{a\\overline b-\\overline a b}\n",
    "$$\n",
    "$q$ is undefined when $a\\overline b-\\overline a b = 0$. That is true iff $a$ and $b$ both lie on a diameter, in which case the line is just a straight Euclidean line. If we want to confuse people, we can refer to that line as a circle of infinite radius. The point is that if $\\left|q\\right|$ is really big, the arc is instinguishable from a line. We'll use this fact later to avoid computational issues when $a$ and $b$ are very close to sharing a diameter.\n",
    "\n",
    "The radius is just $r = \\left|a-q\\right| = \\left|b-q\\right|$. Let's create a function that gives us the center and radius from the two points $a$ and $b$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ortho_euc_circle(a, b):\n",
    "    \"\"\" Gives the center and radius of the circle that passes through \n",
    "    a and b and is also orthogonal to the unit circle.\"\"\"\n",
    "    numerator = a*(1+abs(b)**2)-b*(1+abs(a)**2)\n",
    "    denominator= a*b.conjugate()-a.conjugate()*b \n",
    "    center = numerator/denominator\n",
    "    radius = np.abs(a-center)            \n",
    "    return center, radius"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's pick a couple points and try it out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "a = 0.5\n",
    "b = 0.5j\n",
    "pts = np.array([a, b])\n",
    "\n",
    "center, radius = ortho_euc_circle(a,b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#a couple functions for things we'll have to do over and over again\n",
    "def plot_circle_arc(center=0, radius=1, lower_lim=0, upper_lim=2*np.pi, num=360, color='k'):\n",
    "    t = np.linspace(lower_lim, upper_lim, num=360)\n",
    "    unitcircle = radius*np.exp(1j*t)+center\n",
    "    plt.plot(unitcircle.real, unitcircle.imag,c=color)\n",
    "    \n",
    "def initialize_plot():\n",
    "    fig = plt.figure(figsize=(5,5))\n",
    "    #unit circle\n",
    "    plot_circle_arc()\n",
    "    plt.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUkAAAE4CAYAAADW9AHMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XnU3eO9/vH3RRIzoaYaY4hZRFBiVlTMVKuhaA2llKOl\ndYrT/pz21OnR1mlNVdqaWsfQQSuKllLz0NAghhARjRYx1hAhfH5/3HckIdl5hv3se3+/+3qttVe6\n1MpzZT1xPff3e0+KCMzMbPbmKR3AzKyduSTNzBpwSZqZNeCSNDNrwCVpZtaAS9LMrAGXpJlZAy5J\nM7MGXJJmZg24JM3MGnBJmpk14JI0M2vAJWlm1oBL0sysAZekmVkDLkkzswZckmZmDbgkzcwacEma\nmTXgkjQza8AlaWbWgEvSzKwBl6SZWQMuSTOzBlySZmYNuCTNzBpwSZqZNeCSNDNroF/pAGbtQkLA\nQGBpYKmZfl0EmD9/5su/vgu8A7ydf30LeBl4MX9eAJ4DJkcQLf2DWFMpwt8/6ywSCwJD8md1YLWZ\nPu8Bz+fP5Px5FZhKKsK38v+eBxgA9M+/zg8sAXxkps9H8z9/CpiQP48CY4AHI3i1r/+s1nsuSas1\niXmBocDWwCbAhsDKwCOkshoHjM+fJyN4pclff1FgELBK/qxDKuf1SAU8BrgDuA0YHcHUZn596z2X\npNVKfmQeAowAtgE2B54BbgHuAu4HHongnWIhAYl5SCPXDYEtgC2BNYHRwM3AH4B7I3ivVEZLXJJW\neRLzAdsCu+fPNFLJ3AzcEsHkYuG6IY86NwN2AHYlvQ+9FrgGuDaC1wrG61guSaukPBLbBjgA2Jv0\n+Px74GrSSLHyf7ElBgG7kIp/c1JhXgpcF8HbBaN1FJekVYrEYOBQ4LOkWeRLgMsieKZosD4msSTw\nKWA/0vvMK4FzI/hb0WAdwCVpbS9PvuwMHA0MAy4ELongwZK5SpFYEfgccDjpfeuPgSsieKtosJpy\nSVrbkliIVATHkNYdnoXL4H0S/UjvLo8k/fA4BzgzgheLBqsZ77ixtiOxmMRJwJPAcGBkBB+L4GIX\n5AwRTIvgdxGMIM2OrwA8LnG6xPKF49WGS9LahsSiEt8irVlcC9g2gn0juKdwtLYXwbgIDiMtfxLw\noMRZEssUjlZ5LkkrTqK/xNGkhd0rAZtGcFAEjxSOVjkRTIrgK6Q1l+8AD0ucIrFI4WiV5ZK0YiQk\nsTcwFtgD2CmCz0cwvnC0yotgci7LjUmL1h+X+GKeBLNu8MSNFSGxGmmiYXng+AiuLxyp1iSGAmeS\nDug4MoLRhSNVhkeS1lISAyROBu4GbgQ2dEH2vbyechvSD6Zr8vvKgYVjVYJL0lpGYhPS3unNgI0j\nOK30HupOEsF7EVxIOmSjH2lyZ0TZVO3Pj9vW5/J6vpOALwH/Rlrr6L94hUlsD/yctN3xqxG8XjhS\nW/JI0vqUxBrA7aSTboZFcLkLsj1EcCNpydB8wBiJzQtHaksuSeszEvuRCvISYETd91dXUQSvRnAw\ncBzwW4nj83Fzlvlx25pOYgDwA9J+630iGFM4knVBPnXoCtJ+8IObfQBxVXkkaU0lsQLpHMeVSJMz\nLsiKiOApYCtgEjBaYkjZRO3BJWlNI7ExaWnP1cDeHolUTwRTIzgG+AZwo8QupTOV5sdtawqJPYCf\nAV+I4KrSeaz38kTOr4FTIzizdJ5SXJLWaxLHACcCe0Zwb+k81jwSq5Cuj7gB+HIn3rnjkrQey7Og\npwJ7AbtEMKFwJOsDeWfO74C/kyZ0OmoDgEvSeiQX5A9J5xjuFMELhSNZH5JYAPgV6WShkZ10rqcn\nbqzb8kkyPwE+Bmzvgqy/CKaQLlx7GxiVT43vCC5J65Z8S+HPgcHAJzyD3TnyDY37kR67fy8xf+FI\nLeHHbeuy/Ih9Fum2vp0jeLNwJCsgP0lcAiwKfLLu19t6JGnd8W1gU2B3F2TniuBd0m2N7wK/yAeY\n1JZL0rpE4qvAPqQ92P8qncfKyjPcnwEGAufWeb+3S9LmSmIk6VrXHT1JY9PlGe5Pkq6z/XrhOH2m\n1sNk6z2JzYAzgB0imFQ6j7WXCF6X2B24U2J8BFeUztRsnrixOZJYCbgLODyCUaXzWPuS2IC0K2eP\nCO4snaeZ/Lhts5XXwV0NfN8FaXOTT3v6PPAriY8WjtNUHknah+SX8BcBAXzeJ4lbV0mcAmxH2mQw\nrXCcpvBI0mbnENLL+KNckNZN3wamkPb014JHkjaLfNDqjcDWETxSOo9Vj8SSwGjg2Docm+eRpL0v\nv4e8EviKC9J6Ki8T+wzwkzq8n/RI0t4ncQawRAQHlM5i1SfxbWAjYNcqv7ZxSRoAEtuR9uOuH8HL\npfNY9Un0Jy0h+0kE55XO01MuSUNiUeAB4MgIri2dx+pDYh3gFmDTCMaXztMTLklD4hxgQASHlc5i\n9SNxAvBx0slRlSscl2SHyzccjgLW9mO29YX82H0/8J8RXFk6T3e5JDtYPkD3LuCcCC4sHMdqTGJL\n4DJgnaqdIuUlQJ3tMNJx/BeXDmL1FsFtwPXAt0pn6S6PJDuUxOLAo6QrGMaUzmP1lxeZPwoMj+Dx\n0nm6yiPJznUC8HsXpLVKXmT+Ayq2ZdEjyQ4ksRzwILCBz4i0VpJYEHgM+HQEd5XO0xUuyQ4kcS7w\nrwhOKJ3FOo/EIaRj1bapwpIgl2SHkRgM3AGsGcFLpfNY58m3LT4EHBPBDaXzzI3fSXaeE0hLflyQ\nVkS+bfG/gZNKZ+kKjyQ7iMTypHeRa/hCLyspLzAfB+zf7tc9eCTZWb4MXOyCtNLylbSnUYHRpEeS\nHSKvixwPDI3g6dJ5zCTmByaQbuIcWzrPnHgk2TkOBq51QVq7yPd2nw8cWTpLIx5JdoB8sddjwMER\n3F46j9l0EiuQjulbOYLXSueZHY8kO8N2wFTS0h+ztpE3M9wE7XsavkuyM3wROLcKC3etI50NHJWf\neNqOH7drTmIp4HHS48yrpfOYfVAuxyeBT0Zwf+k8H+SRZP3tQ5qwcUFaW8pPOJcC+5fOMjsuyfob\nSTrs1Kyd/RLYL29ZbCsuyRrLO2yGANeVzmLWSAQPA5OBrUtn+SCXZL19Grgqgqmlg5h1wf8B+5YO\n8UEuyXrbC/hN6RBmXXQ1sGu7zXK7JGtKYjFgI+DPpbOYddGjwDvA+qWDzMwlWV87ALdF8GbpIGZd\nkWe5RwG7lc4yM5dkfe0CXFs6hFk3XU2blaQXk9dQfqczCdi2SrfSmeWTgV4Alo3g9dJ5wCPJuloF\nEPBE6SBm3ZFPBroPGF46y3QuyXragvQ+0o8JVkW3AluVDjGdS7KetgAfiWaVdQtttKjcJVlPLkmr\nsjuATST6lQ4CLsnakVgYWA0YUzqLWU/kw3cnAWuWzgIuyTpaF3g0X7RkVlV/A4aWDgEuyTpan3Qc\nvlmVjQE2KB0CXJJ1NIR0t7ZZlXkkaX1mfVySVn1jgXVKhwCXZB2tTroZ0azKJgFL5h04Rbkka0Ri\nALA08EzpLGa9EcG7wN+BlUtncUnWywrAsxFMKx3ErAmeBFYtHcIlWS8rA0+VDmHWJBNwSVqTrQQ8\nXTqEWZP8A1imdAiXZL0sCTxfOoRZk7xA+jtdlEuyXpYAXiodwqxJXgCWKh3CJVkvH8ElafUxGY8k\nrck8krQ6eYn0g78ol2S9LAb8q3QIsyaZAixQOoRLsl76A2+XDmHWJFOBAaVDuCTrpR/4iDSrjanA\nfKVDuCTrpT94t43VxjS8TtKazCNJq5O26Ke2CGFNE/h7avUxDXildAj/B1UvbTEbaNYkbfFk5JKs\nlylQ/vw9syaZD5ekNZlHklYniwKvlg7hkqyXKcCCpUOYNclA/E7Smuwl0tZEszpwSVrTPU+6vsGs\nDhbDJWlN9jxtsPjWrEmWAl4sHcIlWS8eSVqdrEi6DKwol2S9PAt8tHQIsyZxSVrTTQBWKR3CrElW\nwiVpTfYssKDEYqWDmDWBR5LWXBEE8ASwWuksZr0hsTBpOduk0llckvXzBLB66RBmvbQ2MC6Cd0sH\n6Vc6gDXdY8A6AJLmJ72jXBVYjrSkYinSdq/586c/aX/sO6RDTl8mLbt4AXiG9J7zqYh4raV/Cut0\n6wJjS4cAl2RtSFoa2By+MRh23Era+mBgWeBpYDyp8CaT3vG8CryVP9NIfw8GkEpzIOmGuqHAbsAg\nYJCkV4EH8uc+4LaIKP4oZLXlkrTekbQQsCMwAtiaNFK8E24cD//eD9iKNALs9UnlkkR6ib4BMAQY\nCZwp6Q3gFuBa4PqI8E2N1izrAz8uHQJAEVE6g3WRpEWBfYBPkUrwHuAPwM3AmIh4V0KkPdxrRDC5\nD7MIWBPYFtgF2AYYA1wBXBkRz/XV17Z6k5iH9Lpn3Qj+WTyPS7K9SZqHNGI8hDRq/DNwGXBdRMz2\nGCmJm4DvRnB9C3POD2xPGmXuDtwN/BS4KiKKnwlo1SGxFnBdBINKZwHPbrctSYtJ+jLwKPBd4CZg\nlYjYOyIun1NBZncCW7Qi53QR8VZEXBMRB5Ie/S8CjgYmSvqWJO8pt64aTvo73BZckm1G0pKSvk2a\nbNkMOBgYFhHnduOd382kx+AiIuLNiLg0IrYBPkGaUX9U0jmSvCPI5sYlaR8maVFJ3wHGkU7y2TQi\nRkbE7dH9dyJ3AMOk8gfwRsRDEXEksBZpVv2vkn4sadnC0ax9bQ3cVjrEdC7JwiT1l/QlUjkuDwyN\niMMjYnxPf88IXict1RnepJi9FhHPRcSJwBrAm8DY/BhevMitfUgMIu20+VvhKO9zSRYkaTPgr8An\ngRER8fmIeLpJv/1NwA5N+r2aJiJejIjjgWGkwhwraffCsax97AT8MYL3SgeZziVZgKRFJJ0D/IY0\nKbNDRDT7J+co0mLwthQREyNiJPAF4PuSfuvJHSOVZMtWZXSFS7LFJG1BepSYD1gnIv6vB+8cu+Ie\nYGmJVfvg926aiLiBtED9YWCMpH0KR7JCJAYA2wF/LJ1lZi7JFpE0b561/jVwfEQcGhF9dn9HPhjg\namCPvvoazRIRUyPiZGAv4FRJP/e7yo60PfBIBG21EcEl2QKSliI9QmxBmpi5qkVf+vfAni36Wr0W\nEXcBG5H2kd8paXDhSNZa+5J2bLUVl2QfkzQMGA3cC3wiIp5t4Zf/EzBUYrkWfs1eiYjXgQNJ+3Zv\nlzSicCRrgfyovQfwq9JZPsgl2Yck7UYaQR4XESc247CJ7ohgCvBbYP9Wft3eiuRcYG/ggrxEyupt\nB9KjdtudLOWS7COSjgTOB3aPiJI/HS8mjcwqJyJuJ72iOFrS9/KhGlZPBwCXlw4xOz7gog9IOgE4\nAtgxIp4sm4V5SAfn7hHBmJJZekrSEqTj2P4GHBURxU+rtuaR+AhpG+6qEbTdcXseSTaRklNIJ/Zs\nXbogAfKi3EtIe8ArKe9Z34G0+PwiSfMWjmTNdRAwqh0LEjySbCpJJ5He/23fTucpSqxMOk185bxl\nsZIkLUBaJD8ROCwi2mZXhvVMPv90LHBkBH8pnWd2PJJskjy5cCjpEbttChIggonAX6jou8npImIK\naUnTmsAP/Y6yFrYE5iWdcN+WXJJNIOnTwImk7YXFT1Keg7OAo/NP7srKS4R2IZ3MfkLhONZ7xwNn\n5OuQ25JLspckfQw4B9gtIiaUztPATUCQdjVUWj5weHfSrPenSuexnpFYE9gcuKB0lkZckr0gaUXS\nOsRD++CAiqbKP6l/QBrxVl6+qXEP4BxJG5XOYz1yPPDjCN4sHaQRT9z0kKQBpPd8v4uI75bO0xUS\n/UnnVn42gjtK52mGPJI8DdjYtzVWh8QypKtJ+vTCumZwSfaQpNOBwcCeVZpllTgc2DuCnUtnaZb8\nvVgD2KNK34tOJvG/wLwR/FvpLHPjkuyBfEjsmaS7Zyo1epGYD3gc+FQE95TO0wyS+pNG9ZdHxI9K\n57HGJFYkXT+8TgStPMugR1yS3SRpSdLVCCMjom2XLTQi8UXS3d07tvOsYndIWg24i7SI/5HSeWzO\nJM4DXoyoxvtxl2Q3SboceCYijiudpafyu8kHgOMj+EPpPM0i6QjSSeebtfowEesaicGkmxDXaNcd\nNh/k2e1ukLQnMBQ4uXSW3ojgHdIaw+9L6ywraZN85mXVnQe8Au3/nquDnQacXpWCBJdkl+WTsn9E\nOmBhSuk8TTAKJr0H2z5NOndyoqSRpUP1Rr4G4yjgJEkrlM5js5LYiXRVx+mls3SHH7e7KF+9MDhf\nXlV5aeQ47Gm4Zn5YD3gRYAqwckS09ZKMuZH0n8DaEbFv6SyW5EN1HwS+GsHVpfN0h0eSXZAXjR9F\nWvxaF4PgvqlwGekJCIC30z+vvP8Bhucre609HAs8QTqgpFI8kuwCST8DnouIk0pnaZb8DnIiLLxA\nOoTlIOAvtRhJAkg6mBlH1vkveUESq5CuLxkeweOl83SXR5JzIWlt0j7h0+b271ZJLsJD4PUp8JU3\n4LyAjxxeh4LMLgYGkg7DsELyoc8/B/6nigUJHknOlaRLgEci4tTSWfpCHlEOgje+AQs+GlGfk3Uk\n7Qt8Bdjco8kyJI4EPgdska85rhyXZAOSph9Wu1pf3pHdDiSWIl2P8NkIbi4cpynyCeYPA0dExM2F\n43QciUGkx+ytI6jsAn8/bjd2PPDTuhckQD5k4DDgIomBpfM0Q74L57vAv5fO0mkk+gG/AL5X5YIE\njyTnSNIiwNPAehHxTOk8rSJxFrB4BJ8tnaUZJM1P+j5uHhFPlM7TKSROBTYCds73LFWWR5Jzth9w\nUycVZHYCMFTiC6WDNENEvEWaODiydJZOkReNHwQcWPWCBI8k50jSfcCJEXF96Sytlk+MvhXYrQ4n\nBUmavgRl+YiYWjpPnUksB4wG9qvLu22PJGdD0nrAUqTteh0ngseAw4Er84ROpeVrNcZCfc7QbEcS\nC5BO6j+7LgUJLsk5GQlc1skHuEZwFfBL4LJ8alDV/RLq8Z61HeUL5s4HJgDfKRynqfy4/QH5mtLH\ngc9ExOjSeUqSmBf4HfAccFiVz56UtATwFOmR+7XCcWpH4uukM0q3bvc7a7rLI8kPW490D/B9pYOU\nlhf/jqQWx8PFS8DdwA6ls9SNxCeBo4E961aQ4JKcnZ2Ba7xDI4ngdWBX4FCJA0rn6aVRwG6lQ9SJ\nxHbAucAeEdRyJYhL8sN2Aa4tHaKd5HtIdgV+IFV6L/QoYNf8SsV6SWIj4HJg34j6Pnm5JGeSD9bd\nBOozM9csETxMuuf6Qqmaj6wRMR54C1izdJaqy8vERgGH12kme3ZckrPaBHgwIt4oHaQdRXA3sA9w\nqcQ2pfP00K3AVqVDVJnEWsANwIl5FUStuSRntSVwW+kQ7SyCW0mTOVdKlSybW4CtS4eoKol1gRuB\n/4jgwsJxWsIlOavNSDe5WQMR/Jm05vA3FXxHeQ9pT7F1k8QGpBHk1yK4qHSeVnFJzmoI6dJ0m4sI\n/kR6R3mBxH6l83TDI8AqkhYoHaRKJLYArgeOjeDS0nlaySWZSVoM+AjwZOksVRHBnaR1h9+TOKp0\nnq6IiLeBx4D1S2epCol9SNsNPxfBFaXztJpLcoZ1gYc7eStiT0TwIOkd37ES/5t36bS7scBapUNU\ngcSxpKuUd4qg4w57AZfkzFYBxpcOUUURPEl6n7s+cLXEYoUjzc2TpO+3zYFEf4kzgSNIVy/cXzpT\nKS7JGQYBE0uHqKoIXibtVpoA3CmxWuFIjTwJrFo6RLuSWJY0gz0I2Dyis/+7cEnOsDIuyV6J4J0I\nvgScRSrKfUpnmoOJpO+3fYDEcNLZm38m7cWu/dUlc9OvdIA2sjTwbOkQdRDBORL3ApfnRedfi6Cd\nDrudDCxZOkQ7yVe/HgucCBwSwajCkdqGR5IzLAG8VDpEXURwLzAMWAG4XWKNwpFm9gIuyffl08Sv\nBz4NbOaCnJVLcoYlgBdLh6iT/Ki2D3AhcIfEl/OIpbSXSMu9Op7E3qRjAW8jnQXpJXAf4EN3M0lP\nAdvlo/6tySRWBy4AgvQ494SkpUiTA09FxOTWZZGAd4EBETGtVV+3nUgsA/yQdF7BgXnNq81GO/xU\nbxf9gXdKh6irCJ4AtgV+A9wlXXspDJhIukdooqSRrcsSAUwF5mvV12wXEpI4FHiQdFL7EBdkYx5J\nZpKeA4ZExHOls9SddNkmsOhdsOo8cAxpOzBTgJVbNaKUFMCgiOiYFQ0SQ4AzgAWBL0R4C25XeCQ5\nwzxQ3TtcqmU/YNfX4GvAT4ArgVWnkR69W2lKi79eERLLSpxPGrVfCQx3QXadS3KGt4D5S4foEE8B\nA9KZreuSzhS5ZxF45Yi8kLkV3gRqfW6oxEISJ5O2Yb4CrBnB2fnuIusil+QMUwCfDNMC+ZH6EGAK\nvPUq/NcUGHE4LPYvYKzEqRKL99XXlzQP6X3kW331NUqSWFjiBNI22w2Aj0XwNS8M7xmX5AxvAguV\nDtEpIuIy0q6XHYGVI+49P4LjgA2BZYDxEj+QWKEPvvwiwBsRUasRlcSiEieSynEYsEME+0b4TILe\ncEnO4AXGLRYRkyPi3pknayJ4OoJDSSMgAQ9IXCixXhO/9EDg1Sb+fkVJrC7xQ9K++fWB7SIYGcFD\nhaPVgktyhudIWxOtDUTw9zyyXA0YB/xR4laJA6VevxZZiopvHJCYV2IniVGk0/SnAEMj2D9f2mZN\n4pKc4XnSY561kQhejuBU0qP56cD+wN8lfiixkURProddEfh7M3O2isTaEt8lHdJxKnAVsFIEJ0ZU\n88/U7lySM/hkmDaWTxj6bQQ7k3aJvEa683mcxHckhnSjMCtVkvlx+niJe0hHmM0DjIhgowh+GtEZ\nS5lKcUnO8ASweukQNncRTIjgG8Bg0s2NA4CrgcfyCHPEXB7JB9HGx+JJ9JPYVOK/JB4i7ateE/gm\nadR4gt83to533GSS1gKujojBpbNY9+WDMzYEdiId/juUVC43A3cAf50+4pJ0HXB2RFxdJu2s8pUX\nQ4Dt8mcrYBJwDelx+u4IfK1IIS7JTNIA0oLbJSPizdJ5rHckBgLbkwpnc9Kq9YeAu+CoA+Hjh8Kn\n/hTB6y3ONYD0xDIM2Jh0ve1Q4B+kg25vAm6O4PlW5rI5c0nORNJo4EsRcVfpLNZcEgsCG8NL28Af\nvgmffQC0Numg5cdIj99Pz/SZDLwMvNzVA4PzO9GBwLKkScBlgI+SrooYnD8rkN6H3geMzp/78vUX\n1oZckjORdD5wf0ScUzqL9Q1JnwD+IyK2zo+5q5HKayXSxN30X5cEFs+fd0jrKt8BppGOWZtGWse5\nAOnAiIVI21pfIxXvczP9+iTweP5MiODtVvxZrTl8fcOs7gW2BFyS9TWc9I6SvId5XP7MVh4dLkga\nIfbLn3nzB9JOremfKX53WD8uyVn9BThZksJD7LoaDvy4q/9yBEE6CKPWh2HYnHkJ0KzGkZaTDCqc\nw/qApPlIkzi3ls5i1eGSnEkePd4E7FA6i/WJLYGxEeEL36zLXJIfNgrYvXQI6xM7kW4FNOsyz25/\ngKSBpCUgH40Iv4eqiXz51zhgv4j4a+k8Vh0eSX5ARLwC3A2MKJ3Fmmoo6e/76NJBrFpckrN3KXBA\n6RDWVPsCV3rVgnWXH7dnQ9KipEfu1SKi0ucOGkial7SjZpeIeKB0HqsWjyRnIyL+BfyBdHahVd8I\n4BkXpPWES3LOzgW+lC+Nsmr7AnBe6RBWTS6AObsVmIrXTFaapFVJ6yMvL53FqsklOQf5Bf8ZwJdL\nZ7FeOQ44LyJaeiSa1YcnbhqQND/pxPI9I8JLRypG0pKktZHrRsQ/S+exavJIsoGIeAs4DfhG6SzW\nI18DrnBBWm94JDkXkhYgXfa+u0eT1SFpOeBBYEhEPFM6j1WXS7ILJB1BunDq416MXA2SzgHejIiv\nls5i1eaS7AJJ/YAxwNfb5fIomzNJ65OuXl0nIl4onceqzSXZRZJ2Js12D4kI33PcpvJBFjeR3kX6\nhHnrNU/cdFFEXEsaTZ5cOos1NJJ01cJPSgexevBIshvyZMADwDYRMbZ0HpuVpKVJ3589IuKe0nms\nHjyS7IaI+AdpJHlxvqfb2svZwIUuSGsml2T3nQf8E/h/pYPYDJIOBNYFTikcxWrGj9s9IGkZ4G+k\nU65vLhyn40laE7gN2N4n/VizeSTZAxHxHPA54FJJy5fO08nyYv8rgJNdkNYXPJLsBUknA7uRJnLe\nLp2n0+TlPpcC7wIHeqG/9QWXZC/ksyZ/DbwKHOz/SFsr/5Dak/RDymtXrU/4cbsXIuI90l046wLf\nLByno0gaCXwR2MsFaX2pX+kAVRcRb0jaHbhT0jMR8dPSmepO0k7Aj4Ad87Issz7jkmyCiHg2/4d7\nk6QpEfHL0pnqStLmwC9II0hP1Fifc0k2SUSMk/QJ4AZJb0fElaUz1Y2krUjvgA+IiNtL57HO4JJs\noogYK2kEcJ2khSPigtKZ6kLSdqR7avaPiBtK57HO4ZJssogYI2lb4I+SFo+I00tnqjpJnyZtOdzX\ni/et1bwEqI9IWhG4HrgBOC4iphWOVEmSvkK6zGu3iBhTOo91HpdkH5I0ELgMEPCZiHilcKTKyJew\nnQkMB3aJiKcLR7IO5XWSfSiX4m7Ao8C9kjYsHKkS8ij8L6RzIYe7IK0kl2Qfi4hpEXEs6cbFP0o6\nKm+ns9mQtBdwL/Ab0jvI1wpHsg7nx+0WkrQGaYZ2EnCEF0LPIGlh4HRge9ISnzsLRzIDPJJsqYgY\nB2wK3AfcL+lAjyrfvz/oIaA/sKEL0tqJR5KFSBoG/Ax4BTgmIh4qHKnl8jFz3wc+BnwxIv5UOJLZ\nh3gkWUhE3AdsAvyKtJ3xjHxHS+1JWkjSKaT7aCYA67sgrV25JAvKkzpnk04REvCIpG9JWqxwtD4h\naQFJxwDjgDWAYRFxUkS8WTia2Ry5JNtARDwfEccAGwODgPGS/lvSsmWTNYekRfOi8CeAHUi3Ge4f\nERMLRzObK5dkG4mICRFxEOkd3SLAw5IukLRZFSd4JA2W9CPgKWAz0q6ZPSNidNlkZl3niZs2JmlJ\n4GDgcODrGl0QAAACHUlEQVRN4OfAle28dCi/KtgXOAhYE/gpcE5ETCoazKyHXJIVkK+J2I50Cvqe\nwIOkCZ9rI+KJktkA8muB3UnZtiLtV78IuM53/1jVuSQrRtJ8wE7AXsAI4A3gT6QrVe8AJvb1XTu5\nFDcjFfd2wErAdcBVpOJ+tS+/vlkruSQrLL+nHEKaDBkObJH/rwdIi7PHAk8CTwOTujOqk9QfWApY\nBRhMmo1eD9gIWIC0dfBm4CZgdES80/s/kVn7cUnWSC7NFUllth5padEqpJHecsBrwMv58zrpKtZ3\ngSAV34LAQqRyHAi8SJp0eZy0bOcRYDTwlG+GtE7hkuwQkuYFFp/psxAwL+ngZZEmhqZ/JgMvRsS7\nZdKatQ+XpJlZA14naWbWgEvSzKwBl6SZWQMuSTOzBlySZmYNuCTNzBpwSZqZNeCSNDNrwCVpZtaA\nS9LMrAGXpJlZAy5JM7MGXJJmZg24JM3MGnBJmpk14JI0M2vAJWlm1oBL0sysAZekmVkDLkkzswZc\nkmZmDbgkzcwacEmamTXgkjQza8AlaWbWgEvSzKwBl6SZWQMuSTOzBlySZmYNuCTNzBpwSZqZNeCS\nNDNrwCVpZtaAS9LMrAGXpJlZAy5JM7MGXJJmZg24JM3MGnBJmpk14JI0M2vAJWlm1oBL0sysAZek\nmVkDLkkzswZckmZmDbgkzcwacEmamTXgkjQza8AlaWbWgEvSzKyB/w/VXhy77MXBrwAAAABJRU5E\nrkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x74e4dd8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "initialize_plot()\n",
    "plot_circle_arc(center, radius, color='b')\n",
    "plt.scatter(np.real(pts), np.imag(pts), c='k')\n",
    "plt.show() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Well, we have our circle... but we have *all* of the circle, including the part outside of the disk! What if we just want the segment between the two points? Well, We can determine the angle from the center of the blue circle to the two points, and just plot that. We just need to make sure our angles are for the smaller arc of the circle, not the larger one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def smaller_arc(arc_angle):\n",
    "    result = arc_angle.copy()\n",
    "    if abs(result[0]-result[1]) > np.pi:\n",
    "        if result[0]<result[1]:\n",
    "            result[0] += 2*np.pi\n",
    "        else:\n",
    "            result[1] += 2*np.pi\n",
    "    return np.sort(result)\n",
    "    \n",
    "def arc_angle(pts, center):\n",
    "    angles = np.angle(pts-center)\n",
    "    return smaller_arc(angles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUkAAAE4CAYAAADW9AHMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGYlJREFUeJzt3Xv0pmO9x/H3BzMMMYPQiXGIDuMwRNu5MOMsqUghVEZC\n7V2WilrtziUVpZyiTSWN2EmKkTIjlJoRk5HTNFIOYzBGTuPw3X9c1+x+Odx+v/k9z3Pd93N/Xms9\n6wlr5TP8fOa67vs6KCIwM7Pnt0TpAGZmdeaSNDOr4JI0M6vgkjQzq+CSNDOr4JI0M6vgkjQzq+CS\nNDOr4JI0M6vgkjQzq+CSNDOr4JI0M6vgkjQzq+CSNDOr4JI0M6vgkjQzq+CSNDOr4JI0M6vgkjQz\nq+CSNDOr4JI0M6vgkjQzq+CSNDOr4JI0M6vgkjQzq+CSNDOr4JI0M6vgkjQzq+CSNDOr4JI0M6vg\nkjQzq+CSNDOr4JI0M6vgkjQzq+CSNDOr4JI0M6vgkjQzq+CSNDOr4JI0M6vgkjQzq+CSNDOr4JI0\nM6vgkjQzq+CSNDOr4JI0M6vgkjQzq+CSNDOr4JI0M6vgkjQzq+CSNDOr4JI0M6vgkjQzq+CSNDOr\n4JI0M6vgkjQzq+CSNDOrsFTpAGYvRtIywHLAMgM+SwNPA08CC/P3E8CDEfFEoajWhxQRpTNYi0la\nHlgbWGfA5xXAqsAq+Xtp4GHgcVIRLvpeAhgJjMjfywBjSKV5PzAPuA+4A/grMCd/3xYR9/fi12fN\n55K0npAkYA1g4wGf8cDKwGzg9vyZDfwdmEsquPuABTHIH9T891k+//+uTCrZNfNnrfxZF3gEuB64\nIX/+CNwy2L+PtYdL0roil9VrgTcB2+bvJYAZwHUDPnMi4pkC2dYANgQ2yp/NSFP63wJX5u8ZEfFU\nL7NZ/bgkrWMkvQSYCOwB7Ao8BkwDpubP7DqP1CStDmwNbEMq9pcDlwIXA5d4it5OLkkbFkkrAnsD\ne5EK5nfARcDPI2J2yWzDJemVpLLfDdgO+DPwY2ByRNxTMpv1jkvShkzS0qTy2B+YAEwBfkIabT1U\nMlu35F/z9sC+wFtIzzDPAS7o11+zJS5JGzRJ6wEfBA4AZgLfB86PiPlFg/WYpFGk0eV+wJuB84GT\nI2J6yVzWHS5JqyRpSdKo8XBgE+AM4NSImFMyV11IWg14H3AocC/wHeBHXqvZP1yS9rwkjQTeA3wc\neAA4ifQs7vGiwWoq/2ayC3AEsAHwDdJvJg8XDWbD5pK0f5Onku8DjgZuAr4QEdPKpmoWSeNJv7ns\nAJwMnOg3483lvdsGgKQRkg4jLeieCLwjInZyQQ5dRPwpIvYFtiQtI7pZ0jGSlisczRaDS7LllLyV\n9CLm7cDuEbFnRFxbOFrjRcStEXEIsAVpwfotkiZJ8pkJDeLpdotJ2gQ4AViRNL2+pM6LvZtO0mbA\nV0ijyyMi4vLCkWwQXJItlA+V+BzwLuCTwJkR8XTZVO2Qt0S+BTgRuBr4aETcXTaVVfF0u0Xy1Ppt\nwCxgBWBcRJzuguydSC4ExpFOJbpB0ofy23GrIY8kW0LSKsCpwOuAD0TE1MKRDJD0WuC0/IcHNX0r\nZz/ySLIFJO1GOhbsVmC8C7I+IuIvpF07PwV+L+nQPCW3mvBIso/lU3mOB3YCDvRynnqT9HrgbNIZ\nmgdGxNzCkQyPJPuWpNcBfwBGARu5IOsvImaRlgvNAGZI2rpwJMMl2Zck7UM6x/G4iDgwIhaUzmSD\nExFPRsSxwCTgfElHefpdlqfbfUTSCOA4YE/SjpkZhSPZMEhaAzgPuAs4ICL+WThSK3kk2SckjQF+\nCbwG2NQF2XwR8TfSKenzgN/mk9Otx1ySfUDSWOAq0vrHPSLigcKRrEMiYiFp6v0D4BpJbygcqXVc\nkg0naVPSzo3TIuJDXhjef/IC9ONJx7BdkvfaW494o32DSZpIukLgkIj4aek81l0R8VNJdwIXSRod\nEWeVztQGLsmGkrQH6ZTwvSLit6XzWG9ExHRJ2wNTJC0fESeVztTvXJINJGlv0knhu0XEH0rnsd6K\niL9I2hb4VS7KL5XO1M9ckg0jaT/SLpodI+L60nmsjIiYk4vyMklLRsTnS2fqVy7JBpG0F6kgd8i7\nM6zFIuIuSROAaZL+GREnlM7Uj1ySDZH/YzgV2NkFaYtExN2SdiAV5SMRcXrpTP3GJdkAkrYAfgS8\nzYvE7dki4m95pcMVkh6OiHNLZ+on3pZYc5LWA64EDo6IX5TOY/UlaQPgcuDtEXFl6Tz9wovJa0zS\nysDFwLEuSHsxETET2B84T9K6pfP0C48ka0rSSOBSYHpEHFU6jzWHpEnAUcAWvu97+FySNZSPxjoN\nWI20WNxbDW1IJB0HbAZMjIinSudpMk+362kSsDnwbhekLaZPAAsBr58cJo8ka0bSxsAUYOuIuLl0\nHmsuSS8FpgMfyjc02mLwSLJGJI0mHbJ6pAvShisi5gH7AKdLWqd0nqbySLIm8nPI84C5EfHB0nms\nf0g6AjiY9CJnYek8TeORZH0cBLwa+EjhHNZ/vg3cA3yqdJAm8kiyBvJdJtNJe7JvKJ3H+o+klwN/\nAvaMiN+VztMkHkkWlqfZZwAnuCCtWyLibuBw4GxJy5XO0yQuyfI+AIwGvlI6iPW3iPgJcC3whdJZ\nmsTT7YIkvQK4AdgmIm4qncf6X14WdCPpNKnrSudpAo8kyzqedIGXC3KIJK0iaTNJq5TO0iR5WdAx\nwCmSliydpwlckoVI2g7YCk99hkzSvsAdwGXAHfmPbfC+BzwJHFI6SBN4ul1APrziT6TTff63dJ4m\nySPHO2D3UbAOcCLAY8DYiLivaLgGGXCs2uvz6NJegEeSZRwG/A3wNbBDtyawEGYCnyQVJQvzn7dB\nyseqnUf6h2gVPJLsMUnLA7eRTmfxkp8h+tdIklFp3f3OwE6PQXgkOUSSVgNmAZtFxOzSeerKI8ne\n+ygwxQW5eHIRvhd4DL75EKwScMaZLsihi4h7Sc8rfFJQBY8ke0jSqsBNwKYR8dfSeZosjyjXhGlj\nYJvvA+Mi8AGzQyTpJcAtwB4RMb10njpySfaQpK8DIyLiyNJZ+onECcDKERxQOksT5QMwJkTEW0tn\nqSOXZI/kRby3AOtHxF2l8/QTieVIqwU+FsEFpfM0jaRRwGxgx/xCxwbwM8neOQK4wAXZeRE8ArwH\n+I7EaqXzNE1EPAacAHy8dJY68kiyB/KBAnPwaeNdJfFFYH1gzwj8gz0EklYAbgc2j4jbS+epE48k\ne+P9wFQXZNf9N7AG6YBZG4KIWACcgs8zfQ6PJLtM0hLArcB+Psev+yTWB34DbBuB98QPgaRXklbp\nj42Ih0vnqQuPJLtvArAA+H3pIG0QwZ9JNwVOlli2dJ4miYh/kLYqepXAAC7J7jsUOCU8ZO+lM4Dr\nyRu7bUi+AxyeD4M2XJJdlc+L3B44p3SWNskvbQ4DtpXYv3SehrkCELBt4Ry14ZLsroOByX6+03sR\nPAzsDXxDYoPSeZoiz3jOIC2pMvzipmvydOXPwKSIuKp0nraSeDdpb/IbI/CRYIMw4AXOKyLi8dJ5\nSvNIsnvWB14CXFM6SJtFcA7pSLDJEiNK52mC/ALnOmC30lnqwCXZPfsCP46IZ0oHMY4BHge+VjpI\ng5wD7Fc6RB14ut0Feap9G/AOX7ZUDxJjSMuwvhbBaaXz1J2kMaSDoV/Z9mfqHkl2x0ZAkA5dsBqI\nYD6wB/BZiV1L56m7iJhPelQ0sXSW0lyS3bErcLHXRtZLBLcAewFnSWxaOk8D/BzYvXSI0lyS3bEr\n8IvSIey5IriGtJf+ZxJrl85TcxcDu+Wtta3V6l98N0haEdgQmFo6iz2/CC4kLQu6RGLV0nnqKt97\nMw/aPep2SXbejsA0ry+rtwi+A5wLTJFYqXSeGruE9DPdWi7Jznsz6ZAAq79Pk/5dXSKxQukwNTUV\n2KZ0iJJckp23FfDb0iHsxeU93kcB04GL8zUQ9u+uAraQtFTpIKW4JDsory1bCy/9aYxclIeT1rVe\n6OPV/l1E3E+653x86SyluCQ7awvgjxHxZOkgNngRPEN643038AuJ5QtHqptptPhUIJdkZ21Jmp5Y\nw0TwNHAg6UbLKXmHjiXX0uI33C7JzhoPzCgdwhZPHlEeStq++GuJlxaOVBd/Iu0iayWXZGdtQDpi\nyhoqP6P8L9LSl6kSqxeOVAc3AWvn+7lbxyXZIflKzpeSLnm3BosgIjgG+B/g6rYf2hsRC4GbScf/\ntY5LsnPWB2ZFxNOlg1hnRPBV4GjgcontS+cp7HpaOuV2SXbOOODG0iGssyL4EfBO4Fyp1ecr3gK8\nunSIElySnbMWcHvpENZ5EfyGdKHbFyU+I7Xyv5vZ0M4DQdr4L7tb1iQturU+lO/z/g/SPernt3At\n5WzSQKB1XJKdMxaYUzqEdU8E95BGlPcD10isUzhSL3kkacO2Ji7JvhfBE8AhwMmkN99tOSFnHrB0\nXsXRKi7JDsib/1cD7iqdxbovLxH6NrAP8D2Jz0v09QEQ+ZT9ucAqpbP0mkuyM8YAD3n5T7tEMBXY\nBHgjaYfOKwtH6rb7oH27kFySnbES8EDpENZ7EdwL7AxcCkyX2KVwpG6ah0vSFpNLssUieCaCL5Cm\n36dKnNSnZ1POw9NtW0wuSSOCaaRdKaOB6yQ2Lxyp0x6A9l114ZLsjBWAVl/gbkkED0ZwAPAJ4KcS\nX5QYWTpXhzwOLFM6RK+5JDtjBLCwdAirjwjOJ40qxwEzJLZa9NckrSJpM0lNm7o+AX1T+IPW18sW\nemgp4KnSIaxeIrhX4q3A3sB5Ej+D9f8AfIv0m+pISe+NiHOLBh28J6B9F6Z5JNkZSwG+ssGeI6+p\nnAy8HhaMhMtOh31GkZ5bjgLObNCIcgwtPOTCJdkZI/BI0ipEMB9Gnwx7PwJHAisv+ksLSbu1mmAb\n4B2lQ/Sap9udEYBKh7DamwNXLfmsa6xH0pztrOcD15QO0WseSXbG46Spk9kLioj7gPcCjwEP5e/3\n5j/fBK189u6RZGc8hkvSBiEizpV0OflAlAYVJKRRr0vSFotL0gYtF2OTynGRFWhm7mHxdLszXJLW\nBqOB+aVD9JpLsjPmgy+zt743hvQstVVckp0xF1i1dAizLhuDR5K2mOYCq0ryMiDrZysBD5YO0Wsu\nyQ6IiEdJO25at2XLWmV14O+lQ/SaS7Jz7iVd4WDWdyQtTZpu31M6S6+5JDvnDtKNiWb96FXAXRHx\nTOkgveaS7JzbaOHmf2uNNYA7S4cowSXZOS5J62dr0Zw95h3lkuwcl6T1s3HAjaVDlOCS7JzbgPVK\nhzDrknHArNIhSnBJds7NwFhJy5YOYtYFHkna8ETEQuAmYMPSWcw6SdIKpIXkcwpHKcIl2VkzgE1K\nhzDrsI2BG9q4/Adckp3mkrR+tAUtPJF8EZdkZ00HNisdwqzDWl2SiojSGfqGpBHA/cCaEfFA6Txm\nw5UPbbkX2CQiWrdvGzyS7KiIeBK4mmfd9GTWYGsDT7S1IMEl2Q1XAG8unMGsUyaQfqZbyyXZeVfg\nkrT+sRNwaekQJfmZZIfl55JzgddFROuOlbL+kX+W7wPWi4i5pfOU4pFkh+XnkpcCu5fOYjZMmwOz\n21yQ4JLslp8Be5YOYTZMu9DyqTZ4ut0VklYkHcL78oh4pHQes6HKS39uAd4VEX8snackjyS7ICIe\nBP4A7Fg6i9liGg8sSdog0Wouye6ZDLyrdAizxbQPMDk81fR0u1vylHsOMDYiWndXsTVXnmrfBuwd\nETNK5ynNI8kuyVPuXwHvKJ3FbIi2BJ4CrisdpA5ckt11NvCe0iHMhugQ4LueaieebneRpJGky9y3\niohbS+cxezGSxpAeE7V6AflAHkl2UT6t/Azg8NJZzAbp3cBlLsh/8UiyyySNJR3GOzYi/lk6j9kL\nyS9srgOOjogppfPUhUeSXRYRdwBT8bNJq7/tgZGkF46WuSR741vAkZL8z9vq7Gjg+LbeZfNC/B9t\nb1wBPAq8pXAOs+claSNgA+CHpbPUjUuyB/JSis8Dn8rPfczq5ijgmxHxROkgdeOS7J0LSc97dikd\nxGwgSeuSfi5PKZ2ljlySPZKf83g0aXX0GeAb3j77/FySvfUTYDSwa+kgZgCSNiS91T6xdJa6ckn2\nUEQ8DXwM+KqkpUrnMQM+B3zZa3hfmEuy934O3AO8v3QQazdJ2wAb42eRlbzjpgBJmwAXA6+JiAWl\n81j7SFp0oO6XI+Lc0nnqzCPJAvIZfVOAY0tnsdY6FJgP/Lh0kLrzSLIQSS8DZgLbR8TM0nmsPSS9\nFJgF7OCfvRfnkixI0qHAgcDW3gpmvSLpu8CjEfGh0lmawNPtsk4HngEmlQ5i7SBpR2Ai8MnSWZrC\nI8nCJI0j7e3eOCL+XjiO9TFJK5Ae8UyKiNbfpz1YLskakHQssB2wo6fd1i2STgGWiggvPxsCT7fr\n4SvAKODDpYNYf5I0kbTT66OlszSNR5I1IWlt4Pf4bbd1WF5JMQPYPyJ+XTpP03gkWRMRMZt06Ok5\nkpYtncf6Q140/kPgdBfk4vFIskby6UBn5z98j6/0tOGS9ClgB9KayKdL52kil2TN5FHkNaTf+U8q\nnceaS9IOpFHkGyLiH6XzNJVLsoYkrUMqyr0i4qrSeax5JL0auArYNyJ+UzpPk/mZZA1FxO3AQcBk\nSWsUjmMNI2k0cBHwaRfk8HkkWWOSPgIcTNq2+FDpPFZ/+UXNRcDsiDiidJ5+4JKssfwi55vA64Bd\nI2Jh4UhWY/nn5STgNcAuEfFk4Uh9wSVZc3lkcAHpWKuD/MbbXoikzwC7A9v5nNLO8TPJmsvLNt5N\nGh183ZeI2fORdCTwLtII0gXZQS7JBoiIR0hXfr4J+KKL0gaStB9pI8KOETG3dJ5+48uoGiIiHszH\nXP0GeJx0Dai1nKT9geOAiRExp3CcvuSSbJCImCdpAnCFpKeBL/gZZXtJOph0l/uEiJhVOk+/8oub\nBpL0ctIdOVOAo1yU7SNpEvApUkHeXDpPP3NJNpSklUg3Lt5EOkT1qcKRrAfy8+iPAYeR9mPfVjhS\n33NJNpik5UjLgx4B9ouIxwpHsi7Ky8G+BWxFWjfr/dg94LfbDZbfer+F9CLnijwNtz6UDz65AFgX\n2MYF2TsuyYaLiCeA/YCfA7+XtHHhSNZhkl5FugdpPrCb10H2lkuyD0TyOeAjwBRJbyudyTpD0puA\na4HzSTuuvDW1x/xMss9IegNpWnYe8Anv322m/ILmw8DHgQMi4rLCkVrLJdmHJK0MnAWsTDpP8I7C\nkWwIJI0BTgHWA97mReJlebrdhyLiftILnfOBayXtWTiSDZKkbYHrgXnAVi7I8jyS7HOStgR+AEwD\n/jMi5heOZM9D0kjgv0nnh74/Ii4um8gW8Uiyz0XE1cCGpLWUMyXtUjiSPYukTUnXCW8EjHdB1otH\nki0iaXvgDNIhGUdHxLzCkVpN0kuAz5KWcB0F/MBbTOvHI8kWyfcubwgsAGZJOizv4rAeUrI78GfS\ny7VxEfF9F2Q9eSTZUpI2IB31vzxweERcUzhSK+R/7l8DxgJHeGlP/Xkk2VIRMRN4M3A8cJ6kn0h6\nbdlU/UvSyySdBvwK+BmwvguyGVySLZZ36pxDWo93LTBN0pm+xrZzJK0i6UvAjaTHHK+NiJO8yL85\nXJJGRDwaEceRyvJu4DpJJ+cL7m0xSFpV0nHAzcAYYJOIOCoiHiwczYbIJWn/LyLmR8SxpCts7wd+\nJ2lyXqJigyBpPUnfJJ3zuSywUUQc5l1PzeWStOeIiLkR8UlgLeBq4AJJV0jaV9LShePVjqQlJO0s\n6RfAlaRp9QYRcURE3Fk4ng2T327bi5I0AtgLmERaQnQWcFpE3Fo0WGGSVgcOAA4iLdY/ETg3Ih4v\nmcs6yyVpQyJpXeD9pGK4HfgxcF5E3FUyV69IWp70G8aBwHhgMnA28Duvc+xPLklbLHl0OQF4J+kw\njRtIhXFxvz1/k/QyYA/grcA2pH3wZwEXedTY/1ySNmySlgF2At6evx8ALgF+CVzZtLt38q9nc2A7\nYCLpRdYlwIXALyPioYLxrMdcktZRkpYANgZ2BnbJ//sG0gugq4GrI+LucgmfS9JqwBuATYE3AW8E\nZpH2uF8OTPWJ4O3lkrSuyjc6bgpsOeCzkLRveWb+vpH0fPP+bj7Xy9fwrjvgM55UjssCM4DppKn0\nlb5HxhZxSVpP5WsJXgVsAKw/4LMmsDTwt/y5k7RW88EBnwXAk8DTwFP5G2AUqegWfUYDLwNWG/A9\nFhgB3DrgM5NUjH/1Sxd7IS5Jqw1JKwCrA2vkz0qk3Sor5s9oYClgyQHfAI8+67MAuAe4d8D3ncBc\nl6ENlUvSzKyCd9yYmVVwSZqZVXBJmplVcEmamVVwSZqZVXBJmplVcEmamVVwSZqZVXBJmplVcEma\nmVVwSZqZVXBJmplVcEmamVVwSZqZVXBJmplVcEmamVVwSZqZVXBJmplVcEmamVVwSZqZVXBJmplV\ncEmamVVwSZqZVXBJmplVcEmamVVwSZqZVXBJmplVcEmamVVwSZqZVXBJmplVcEmamVVwSZqZVXBJ\nmplVcEmamVVwSZqZVXBJmplVcEmamVVwSZqZVXBJmplVcEmamVVwSZqZVXBJmplVcEmamVVwSZqZ\nVXBJmplVcEmamVVwSZqZVXBJmplVcEmamVVwSZqZVXBJmplVcEmamVVwSZqZVXBJmplVcEmamVVw\nSZqZVXBJmplVcEmamVVwSZqZVXBJmplVcEmamVVwSZqZVXBJmplVcEmamVVwSZqZVXBJmplV+D+K\nV0mc2yCeUwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x75126d8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "angle = arc_angle(pts, center)\n",
    "\n",
    "initialize_plot()\n",
    "#\"line\" between a and b\n",
    "plot_circle_arc(center, radius, lower_lim=angle[0], upper_lim=angle[1], color='b')\n",
    "plt.scatter(np.real(pts), np.imag(pts), c='k')\n",
    "plt.show() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, but what if we want the line to extend out to the limit points? We could just calculate the whole circle and get rid of the part that's outside the unit circle. But we wind up doing a lot of processing just to throw it away, and unless we use a bajillion steps when we plot the circle, there tends to be a gap where the two circles meet. The straighter the line, the worse it is. \n",
    "\n",
    "Instead, we'll calculate the angles from the center of the orthogonal circle to where the two circles intersect. The system of equations is:\n",
    "$$\n",
    "\\left|z\\right| = 1\n",
    "$$\n",
    "$$\n",
    "z = r \\exp \\left(i\\theta\\right) + q,\n",
    "$$\n",
    "and we solve for $\\theta$. Without getting into all the tedious algebra, making the substitution $\\theta = 2 \\arctan(w)$ results in a quadratic in $w$. The formula for $\\theta$ is:\n",
    "$$\n",
    "\\theta = 2 \\arctan\\left(\\frac{2 r \\Im \\left(q\\right) \\pm \\sqrt{4 r^2 - \\left(1-\\left|q\\right|^2+r^2\\right)^2}}\n",
    "{1 + 2r \\Re \\left(q\\right) -r^2  - \\left|q\\right|^2}\\right)\n",
    "$$\n",
    "Note that these are the angles with respect to q, not the origin. Let's try it out..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ext_angle(center, radius):\n",
    "    absc_sq = abs(center)**2\n",
    "    denom = 1 + 2*radius*center.real- radius**2 -absc_sq\n",
    "    numer_1 = 2*radius*center.imag        \n",
    "    numer_2 = np.sqrt(4*radius**2-(1-absc_sq+radius**2)**2)\n",
    "    angle = 2*np.arctan((numer_1 +np.array([-1,1])*numer_2)/denom)\n",
    "    return smaller_arc(angle)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUkAAAE4CAYAAADW9AHMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHIRJREFUeJzt3XnYnuOd//H3B0nEltirlqCoLdbErtYQe7UYqtafZZSa\n/tRoVXVay2+mxqganWaEOhjViuVXRWnUNImdJpYQSwiqtUQQlESI7/xxXpl5OuKSJ7nv+7yu+/q8\njuM+njp6ePIJj0/O87rORRGBmZnN3UK5A5iZVZlL0syshEvSzKyES9LMrIRL0syshEvSzKyES9LM\nrIRL0syshEvSzKyES9LMrIRL0syshEvSzKyES9LMrIRL0syshEvSzKyES9LMrIRL0syshEvSzKyE\nS9LMrIRL0syshEvSzKyES9LMrIRL0syshEvSzKyES9LMrIRL0syshEvSzKyES9LMrIRL0syshEvS\nzKyES9LMrIRL0syshEvSzKyES9LMrIRL0syshEvSzKyES9LMrIRL0syshEvSzKyES9LMrIRL0sys\nhEvSzKyES9LMrIRL0syshEvSzKyES9LMrIRL0syshEvSzKyES9LMrIRL0syshEvSzKyES9LMrIRL\n0syshEvSzKyES9LMrIRL0sysxCK5A5h9GkmLAosDi/b49ANmAx8As4qv7wNvRsT7maJaF1JE5M5g\nDSZpSWBN4HM9Pp8FVgCWL772A94BZpKKcM7XhYC+QJ/i66LAQFJpvg5MA14DXgCeA54vvj4TEa93\n4vdn9eeStI6QJGA1YFP4xoHw0bpw0fLAssAU4NniMwX4EzCVVHCvAW/HPP6gFr/OksX3XZZUsqsX\nnzWKz9rAu8AjwKPF5w/A0/P661hzuCStLYqyWhfYAfhC8XUhYAJ8byYcuQqseQjwfER8lCHbasBG\nwMbFZyhpSn8XcGfxdUJEfNjJbFY9LklrGUlLAMOAfYA9gRnAOGBs8ZkSESGxBnBXBCtnCzsXklYF\ntgO2JxX7SsBvgVuA2zxFbyaXpC0QSUsDBwL7kwrmPuAm4OaImDL3vweRptEbRfBSp7L2lqSVSWW/\nF7AT8BhwDTAqIl7Jmc06xyVpvSapH6k8vgrsCowGriONtt6at+/BrcCICG5sW9AWKn7POwMHA/uS\nnmFeDdwwr79nqyeXpM0zSesAXwMOAyYC/wFcHxHTe/+9OAtYKILvtjZl+0nqTxpdHgrsCFwP/DQi\nxufMZe3hxeRWStLCkvaRdBvpZcYMYPOI2DEiLpufgizcD2zTsqAdFBEzIuK6iNif9HJqCnCDpAck\nHVmMOq1LeCRpcyWpL3A48G3gDeBi0rO4ma35/gwgLfVZLoLaL/6WtDCwB3ASMBj4EfDvEfFO1mC2\nwDyStL8iqb+kk4BnSC9kjo6ILSLiylYVJEAEbwFPAlu06nvmFBGzI+LmiBhOmooPAaZIOkvSspnj\n2QJwSRoAkvpIOoG0oHsYcEBE7B4R49r4y44hPdPrKhHxcEQcTHqcsBLwlKTvSFo8czSbDy7JhlPy\nRdKLmC8De0fEfhHxQAd++TF0YUnOERGTI+JYYGvSgvWnJR0nyWcm1IifSTaYpM2AC4GlgdNIS3g6\n9gPR47nkChHM6NSvm4ukocAPSaPLkyLijsyRbB54JNlAkpaUdCFwK2kZzyYRcWun9y0XzyUfootH\nkz1FxIPALqSXYZdJulrSSplj2adwSTZIMbX+EjAJWArYICJGRsTsjLF+Q3rR0QiR3AhsQDqV6FFJ\nJxdvx62CPN1uCEnLA/8OrAf8bUSMzRwJAInBwK+BNSNo3A+jpHWBS4q/PPKTtnJaPh5JNoCkvUjH\ngk0mTa0rUZCFx4CFSeXdOBHxJOlxw6+A+yUdX5xSZBXhkWQXK07lOR/YHTiizct55pvET4HnIjgv\nd5acJK0PXEk6/OOIiJiaOZLhkWTXkrQe8CDQH9i4qgVZuAE4IHeI3CJiEmm50ARggqTtMkcyPJLs\nSpIOAn4CnBYRl+fO82kkFgFeAraM4LnceapA0p7A5cA/A//iE9PzcUl2EUl9gPOA/Ug7ZiZkjjTP\nJEYAU5o+5e5J0mrAtaQ/QA6LiL9kjtRInm53CUkDSesePw8MqVNBFkYBB+UOUSUR8UfSKenTgLuK\nk9Otw1ySXUDSIOBu0vrHfSLijcyR5sc4YFWJtXIHqZKImAUcB1wF3Ctp88yRGsclWXOShgD3AJdE\nxMmZF4bPtwg+BH5BOp7NeigWoJ9POobttmKvvXWIn0nWmKRhpCsEjo2IX+XOs6AkNibdj7NGBLUs\n+3YrRpI3AadHxBW58zSBR5I1JWkf4OfA/t1QkAARPEJaI7hL7ixVVVwRsTNwdnHup7WZS7KGJB0I\nXArsFRF35c7TYj8DjsodosqKXTpfAL4h6fTcebqdp9s1I+lQ0i6a4RHxSO48rSaxDOnOmDUjqOML\nqI6R9FngduAXEXFO7jzdyiPJGpG0P6kgd+nGggQoivEm4OjcWaouIl4iXel7hKRv5M7TrTySrAlJ\nu5Je0gyv4RrIXpEYSlo3uZZf4Hy6YtH5OODciBiZO0+38UiyBiRtTVoe8+VuL0iACB4EXgH2zp2l\nDopF58OA70s6OHeebuOSrDhJ65CO0ToiIu7MnaeDLgJOzh2iLiJiMjAcuEjS9rnzdBNPtyusuIr0\nPuCHEXFp7jydJNGXdHL37hFMzBynNiTtRjpubfuiOG0BeSRZUZL6AtcBNzatIAEimAVcDPx97ix1\nEhGjge8Bt/i+79bwSLKCipOpLwFWJC0Wb+TLC4mBpHvAh/gItd6RdB4wFBgWER/mzlNnHklW03HA\nVsBXmlqQABFMJ93L49Fk750OzAK8fnIBeSRZMZI2BUYD20XEU7nz5CaxAvAksEEEL+fOUyeSlgPG\nAycXNzTafPBIskIkDSAdsvp1F2QSwVTS3eDfzJ2lbiJiGumMzpGSPpc7T115JFkRxXPIa4GpEfG1\n3HmqRGJlYCKwYQQv5c5TN8VBGEcBWxfnU1oveCRZHUcCawGnZM5RORH8mXTwxZm5s9TUT0iL8/3P\nbz54JFkBxbay8aQ92Y/mzlNFEsuRnk1uGcGzufPUjaSVgIeB/SLivtx56sQjycyKafZlwIUuyE8W\nwTTSLpyzcmepo4h4GTgRuFLS4rnz1IlHkplJOoH0vGgbr2crJ7EkMBkYHsHDufPUkaSrgGkR4VOD\n5pFLMqPiPMBHSVvInsidpw4k/hY4GNgpAv/w9lKxLOhx0mlSD+XOUweebud1PukCLxfkvBsJLA2T\nj5Y0VNLyuQPVSbEs6DvACEkL585TBy7JTCTtBGwLnJs7S52k8yXP/BX0GwmL3Q684OPBeu1y4APg\n2NxB6sDT7QyKwyseBs6IiP+fO0+dFCPHF+Ca/vAYcDbADGBQRLyWNVyNSBoM3AGsX4wu7RN4JJnH\nCcAfSedEWu+sDsxK27lPBtYk/TWr54tUPxExkbR54bu5s1SdR5IdJmlJ4BnS6Sxe8tNL/zOSpH8q\nyl2B4TMgPJLsJUkrApOAoRExJXeeqvJIsvO+CYx2Qc6fogiPBmbAj96CFQJ+drkLsvci4lXgx/ik\noFIeSXaQpBWAJ4AhEeHzERdAMaJcHcYMgB2uIu3r9rO1XpK0BPA0sE9EjM+dp4pckh0k6QKgT0R8\nPXeWbiJxAbB8BIflzlJHxQEYu0bEF3NnqSKXZIcUi3ifBjYs7ku2FpFYHHgEODXCL8N6S1J/YAqw\nW/FCx3rwM8nOOQm4wQXZehG8CxwOjJBYMXeeuomIGcCFwLdzZ6kijyQ7oDhQ4Hl82nhbSZwLbATs\n6y2LvSNpKdJ9QltFhE9Z6sEjyc44Bhjrgmy7HwArA/8nd5C6iYi3gRH4PNOP8UiyzSQtRDq55lCf\n49d+EusDY4HtIvAfSr0gac4J8IMi4p3cearCI8n22xV4G7g/d5AmiGAScAZwrcRiufPUSUT8mbRV\n0asEenBJtt/xwIjwkL2TRpJGRBflDlJD/wacWBwGbbgk26o4L3Jn4OrcWZqkeGlzPLCd5FFRL40B\nBHwhc47KcEm211HAKD/f6bwI/gIcCFwgsWHuPHVRzHguIy2pMvzipm2K6cpjwHERcXfuPE0l8VXS\nW+8tIng9d5466PEC57MRMTN3ntw8kmyfDYElgHtzB2myCK4Crie9yOmTO08dFC9wHgL2yp2lClyS\n7XMwcE1EfJQ7iHE6MBO4IHeQGrkaODR3iCrwdLsNiqn2M8ABvmypGiQGAvcBF0RwSe48VSdpIOlg\n6JWb/kzdI8n22BgI8LWnVRHBdGBf4CyJPXLnqbqImE56VDQsd5bcXJLtsSdwi9dGVksETwP7A1dI\nDMmdpwZuBvbOHSI3l2R77An8JncI+7gI7iXdEvhric/lzlNxtwB7FVtrG6vRv/l2kLQ06SSasbmz\n2NxFcCNwFnCbxAq581RVce/NNGj2qNsl2Xq7AeO8vqzaIhgB/AIYLbF07jwVdhvpZ7qxXJKttyPp\nkACrvn8g/bu6TWKp3GEqaiywfe4QObkkW29b4K7cIezTFXu8TwXGA7cU10DYX7sb2FrSIrmD5OKS\nbKFibdkaeOlPbRRFeRLpzM8bfbzaX4uI10n3nG+SO0suLsnW2hr4Q0R8kDuIzbsIPiK98X6ZNKJc\nMnOkqhlHg08Fckm21jak6YnVTASzgSNJI8rRxQ4dSx6gwW+4XZKttQkwIXcImz9FUR5POkX+PyWW\nyxypKh4m7SJrJJdkaw0mHTFlNVU8o/y/wK3AWIlVMkeqgieANYv7uRvHJdkixZWcy5EuebcaiyAi\nOAO4HLhHYnDuTDlFxCzgKWjm4cUuydbZEJgUEbNzB7HWiOB84DTgDomdc+fJ7BEaOuV2SbbOBsDj\nuUNYa0XwS+Ag4BcSX8mdJ6OngbVyh8jBJdk6awDP5g5hrRfBGNKFbv8o8QOpkf/dTAHWzB0ihyb+\ny26X1UmLbq0LRfA4sAWwC3BdA9dSTiENBBrHJdk6g4Dnc4ew9ongVVJJvk56odOkkZVHkrbAVscl\n2fUieB84DhgB3Cuxe+ZInTIN6Fes4mgUl2QLFJv/VwReyp3F2q9YIvQT0gudyyTOkejqAyCKU/an\nAsvnztJpLsnWGAi85eU/zRLBWGBzYEvSMqHPZo7Ubq9B83YhuSRbYxngjdwhrPOK55TDgduB8RLD\nM0dqp2m4JG0+uSQbLILZEZwDHAKMlLi4S8+mnIan2zafXJI2Zz3lRsAA4CGJrfImark3SD/rjeKS\nbI2lgEZf4G5JBG9GcBjwHdIhvmdL9M2dq0VmAovmDtFpLsnW6APMyh3CqiOC60h7nTcmjSq3nfP/\nSVpe0lBJdZu6vg9dU/jzrKuXLXTQIsCHuUNYtUTwisR+wJeBURI3wfoPABeT/lDtK+noiPhl1qDz\n7n1o3oVpHkm2xiKAr2ywjynWVF4HbABv94PfXQoH9Sc9t+wP/KxGI8qBNPCQC5dka/TBI0krEcF0\nGPBv8DfvpseVYyhOHptF2q1VB9sDB+QO0WmebrdGAModwirvebhrYdgMOAa4DbhpcXjjL3ljzbPr\ngXtzh+g0jyRbYyZp6mT2iSLiNeBo+GgGXPIWrDcDBo+Gb42TOK0G19k28tm7S7I1ZuCStHlQvKQZ\nBAyD6YMittoL2A4YCkyW+FqFlwz1xSVp88klafMsIl6LiAeLkSURPBXBgcC+wD7AkxKHSyycNejH\nLQW8lTtEp7kkW8MlaQssgvER7EG6//tY4AmJYyT65U323wYA03OH6DSXZGtMB19mb60RwTjgC6S3\nOwcAz0qcIrFE3mTptKvMGTrOJdkaU4EVcoew7lGsrxwXwXDSNHwrYIrEuRKrZoo1EI8kbT5NBVaQ\n5GVA1nIRTIjgIGBbYAngEYnrJHaQOrr0bBngzQ7+epWgdOCwLShJ7wCrRETjpiPWWcUlZIcDJ5F2\nel0CXB3R3pOoJL0KbBoRjTqB3yPJ1nmVdIWDWVtF8E5xfcT6wCmkEeYUiVESe7TjrbikfqTp9iut\n/t5V55JsnRdI69/MOqJ4bvm7CA4hXff6e+D7wB8lzpfYsoXT8VWAlyLioxZ9v9pwSbbOMzRw879V\nQ3GO5U8j2BIYBrwHXAE816LCXA14sRVZ68Yl2TouSauECCZF8D1gPdLi9BmkwnxRYqTE/sVzzd5Y\ng4ZemeySbB2XpFVKMR2fGMGZpMLcBXgcOAF4SeIOib+XGDoPV+JuUPy9jeNTgFrnGWCd3CHM5iaC\nAJ4qPhcWC9N3BnYHLgdWlbgHGFt8HopgZo9vsQHpfLfG8RKgFpHUl7TQdrmIeC93HrPekFiOtMtn\nh+Lr54EngD8AD8LQc+DpbSPeejZjzCxcki0kaTxwYkTclzuL2YKQ6A9sAgyBmdvAkwfB5EUjDmzc\nCfyebrfWBNKJqi5Jq7UIZpAO2L1X6v8osHpENK4gwS9uWm1OSZp1k61p4Inkc7gkW2s86fBUs27S\n6JL0M8kWktQHeJ00NWnrPlqzTigObXkV2Cwi/pQ7Tw4eSbZQ8czmHtKtcmbdYE3g/aYWJLgk22EM\nsGPmDGatsisNXR85h0uy9cbgkrTusTvw29whcvIzyRYrnktOBdaLiMYdK2Xdo/hZfg1YJyKm5s6T\ni0eSLVY8l/wtsHfuLGYLaCtgSpMLElyS7fJrYL/cIcwW0B40fKoNnm63haSlSYfwrhQR7+bOY9Zb\nxdKfp4FDIuIPufPk5JFkG0TEm8CDwG65s5jNp02AhUkbJBrNJdk+o4BDcocwm08HAaPCU01Pt9ul\nmHI/DwyKiMbdVWz1VUy1nwEOjIgJufPk5pFkmxRT7t8BB+TOYtZL2wAfAg/lDlIFLsn2upJ0P7JZ\nnRwLXOqpduLpdhsVp5X/Cdg2IibnzmP2aSQNJD0mavQC8p48kmyjiJgFXAacmDuL2Tz6CnC7C/J/\neCTZZpIGkQ7jHRQRf8mdx+yTFC9sHgJOi4jRufNUhUeSbRYRL5Bun/OzSau6nYG+pBeOVnBJdsa/\nAl+X5H/eVmWnAedHxEe5g1SJ/6PtjDHAe8C+mXOYzZWkjYHBwM9zZ6kal2QHFEspzgHOLJ77mFXN\nqcBFEfF+7iBV45LsnBtJz3v2yB3ErCdJa5N+LkfkzlJFLskOKZ7zeDRpVfQD4EfePjt3LsnOug4Y\nAOyZO4gZgKSNSG+1f5w7S1W5JDsoImYD3wL+WdIiufOYAWcD/+Q1vJ/MJdl5NwOvAMfkDmLNJml7\nYFP8LLKUd9xkIGkz4Bbg8xHxdu481jyS5hyo+08R8cvcearMI8kMijP6RgNn5M5ijXU8MB24JneQ\nqvNIMhNJnwEmAjtHxMTceaw5JC0HTAJ28c/ep3NJZiTpeOAIYDtvBbNOkXQp8F5EnJw7Sx14up3X\nSOAj4LjcQawZJO0GDAO+mztLXXgkmZmkDUh7uzeNiD9ljmNdTNJSpEc8x0VE4+/TnlcuyQqQdAaw\nE7Cbp93WLpJGAItEhJef9YKn29XwQ6A/8He5g1h3kjSMtNPrm7mz1I1HkhUhaU3gfvy221qsWEkx\nAfhqRPxn7jx145FkRUTEFNKhp1dLWix3HusOxaLxnwMjXZDzxyPJCilOB7qy+MvDfaWnLShJZwK7\nkNZEzs6dp45ckhVTjCLvJf3Jf3HuPFZfknYhjSI3j4g/585TVy7JCpL0OVJR7h8Rd+fOY/UjaS3g\nbuDgiPh97jx15meSFRQRzwJHAqMkrZY5jtWMpAHATcA/uCAXnEeSFSbpFOAo0rbFt3LnseorXtTc\nBEyJiJNy5+kGLskKK17kXASsB+wZEbMyR7IKK35eLgY+D+wRER9kjtQVXJIVV4wMbiAda3Wk33jb\nJ5H0A2BvYCefU9o6fiZZccWyja+QRgcX+BIxmxtJXwcOIY0gXZAt5JKsgYh4l3Tl5w7A/3NRWk+S\nDiVtRNgtIqbmztNtfBlVTUTEm8UxV78HZpKuAbWGk/RV4DxgWEQ8nzlOV3JJ1khETJO0KzBG0mzg\nXD+jbC5JR5Huct81IiblztOt/OKmhiStRLojZzRwqouyeSQdB5xJKsincufpZi7JmpK0DOnGxSdI\nh6h+mDmSdUDxPPpbwAmk/djPZI7U9VySNSZpcdLyoHeBQyNiRuZI1kbFcrB/BbYlrZv1fuwO8Nvt\nGiveeu9LepEzppiGWxcqDj65AVgb2N4F2TkuyZqLiPeBQ4GbgfslbZo5krWYpFVI9yBNB/byOsjO\nckl2gUjOBk4BRkv6Uu5M1hqSdgAeAK4n7bjy1tQO8zPJLiNpc9K07FrgdO/frafiBc3fAd8GDouI\n2zNHaiyXZBeStCxwBbAs6TzBFzJHsl6QNBAYAawDfMmLxPPydLsLRcTrpBc61wMPSNovcySbR5K+\nADwCTAO2dUHm55Fkl5O0DXAVMA74RkRMzxzJ5kJSX+D7pPNDj4mIW/Imsjk8kuxyEXEPsBFpLeVE\nSXtkjmT/i6QhpOuENwY2cUFWi0eSDSJpZ+Ay0iEZp0XEtMyRGk3SEsBZpCVcpwJXeYtp9Xgk2SDF\nvcsbAW8DkySdUOzisA5SsjfwGOnl2gYR8R8uyGrySLKhJA0mHfW/JHBiRNybOVIjFP/c/wUYBJzk\npT3V55FkQ0XERGBH4HzgWknXSVo3b6ruJekzki4Bfgf8GtjQBVkPLskGK3bqXE1aj/cAME7Sz3yN\nbetIWl7SPwKPkx5zrBsRF3uRf324JI2IeC8iziOV5cvAQ5J+Wlxwb/NB0gqSzgOeAgYCm0XEqRHx\nZuZo1ksuSftvETE9Is4gXWH7OnCfpFHFEhWbB5LWkXQR6ZzPxYCNI+IE73qqL5ekfUxETI2I7wJr\nAPcAN0gaI+lgSf0yx6scSQtJGi7pN8CdpGn14Ig4KSJezBzPFpDfbtunktQH2B84jrSE6ArgkoiY\nnDVYZpJWBQ4DjiQt1v8x8MuImJkzl7WWS9J6RdLawDGkYngWuAa4NiJeypmrUyQtSfoD4whgE2AU\ncCVwn9c5dieXpM2XYnS5K/A3pMM0HiUVxi3d9vxN0meAfYAvAtuT9sFfAdzkUWP3c0naApO0KLA7\n8OXi6xvAbcCtwJ11u3un+P1sBewEDCO9yLoNuBG4NSLeyhjPOswlaS0laSFgU2A4sEfxvx8lvQC6\nB7gnIl7Ol/DjJK0IbA4MAXYAtgAmkfa43wGM9YngzeWStLYqbnQcAmzT4zOLtG95YvH1cdLzzdfb\n+VyvuIZ37R6fTUjluBgwARhPmkrf6XtkbA6XpHVUcS3BKsBgYMMen9WBfsAfi8+LpLWab/b4vA18\nAMwGPiy+AvQnFd2czwDgM8CKPb4OAvoAk3t8JpKK8Tm/dLFP4pK0ypC0FLAqsFrxWYa0W2Xp4jMA\nWARYuMdXgPf+1+dt4BXg1R5fXwSmugytt1ySZmYlvOPGzKyES9LMrIRL0syshEvSzKyES9LMrIRL\n0syshEvSzKyES9LMrIRL0syshEvSzKyES9LMrIRL0syshEvSzKyES9LMrIRL0syshEvSzKyES9LM\nrIRL0syshEvSzKyES9LMrIRL0syshEvSzKyES9LMrIRL0syshEvSzKyES9LMrIRL0syshEvSzKyE\nS9LMrIRL0syshEvSzKyES9LMrIRL0syshEvSzKyES9LMrIRL0syshEvSzKyES9LMrIRL0syshEvS\nzKyES9LMrIRL0syshEvSzKyES9LMrIRL0syshEvSzKyES9LMrIRL0syshEvSzKyES9LMrIRL0sys\nhEvSzKyES9LMrIRL0syshEvSzKyES9LMrIRL0syshEvSzKyES9LMrIRL0syshEvSzKyES9LMrIRL\n0syshEvSzKyES9LMrIRL0sysxH8BxTZz+3M+rowAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x752f358>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "angle = ext_angle(center, radius)\n",
    "\n",
    "initialize_plot()\n",
    "#\"line\" between a and b\n",
    "plot_circle_arc(center=center, radius=radius, lower_lim=angle[0], upper_lim=angle[1], color='b')\n",
    "plt.scatter(np.real(pts), np.imag(pts), c='k')\n",
    "plt.show() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Much better! All we had to do was, well, derive a really ugly formula.\n",
    "\n",
    "In real life, we'll need handling for the case where the points line on a diameter. Here's a function that does most of the goofy handling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_hyp_line(a,b,extend=True, straighten = 1E10): \n",
    "    \"\"\"Plots hyperbolic lines.\n",
    "    Input: \n",
    "        a, b: Points to plot the line between.\n",
    "        extend: Whether to extend the lines to the limit points.\n",
    "        straighten: If abs(center) winds up being above this value, plot a straight line instead\n",
    "    \"\"\"\n",
    "    pts = np.array([a,b])\n",
    "    center, radius = ortho_euc_circle(a, b)\n",
    "    absc = abs(center)\n",
    "    if absc < straighten:#plot a circular arc\n",
    "        if extend:\n",
    "            angle = ext_angle(center, radius)\n",
    "            plot_circle_arc(center, radius, lower_lim=angle[0], upper_lim=angle[1], color='r')\n",
    "        angle = arc_angle(pts, center)\n",
    "        plot_circle_arc(center, radius, lower_lim=angle[0], upper_lim=angle[1], color='b')\n",
    "    else: #plot a straight line instead\n",
    "        if extend and (a != 0 or b != 0):\n",
    "            #normalize the points so they're on the unit circle\n",
    "            anorm = a/abs(a)\n",
    "            bnorm = b/abs(b)\n",
    "            #unless one is zero, in which case use the mirror image of the other one. \n",
    "            #if both are zero the line is degenerate, so don't worry about it\n",
    "            if a == 0:\n",
    "                anorm = -bnorm\n",
    "            elif b == 0:\n",
    "                bnorm = -anorm \n",
    "            plt.plot([a.real, b.real], [a.imag, b.imag], c='r')\n",
    "        plt.plot(pts.real, pts.imag, c='b')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's test it out on a representitive set of points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "n=5\n",
    "\n",
    "#random points on the unit disk\n",
    "pts = (np.random.uniform(0,1,size=n)*\n",
    "        np.exp(2j*np.random.uniform(0,np.pi,size=n)))\n",
    "\n",
    "#move two points to the unit circle\n",
    "pts[:2] = pts[:2]/abs(pts[:2])\n",
    "#put one point on a diameter with another one\n",
    "pts[-1] = -pts[0]*np.random.uniform(0,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUEAAAEzCAYAAACv5LH7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXm4zdX3x1/L5c7XdE2Z5yFEZYqUhFBJCtHwbZCkiBCK\n61yzkKFIaB40ShJCEgqlWRERDcrsGu5kWL8/9kc/leEO55zPGfbreT6Pp3vP2ft95b7P3nvttZao\nKhaLxRKu5HFbgMVisbiJNUGLxRLWWBO0WCxhjTVBi8US1lgTtFgsYY01QYvFEtZYE7RYLGGNNUGL\nxRLWWBO0WCxhjTVBi8US1lgTtFgsYY01QYvFEtZYE7RYLGGNNUGLxRLWWBO0WCxhjTVBi8US1lgT\ntIQ8IlJUROqLSFG3tVgCD2uClqBDDHlF5Lz/fkXkFmAHsBTY4fy3xfI3YsvrW9xCROKBCkAZoITz\nFHeegkABIL/zZywQCeRz/jyJ+RBX4ARwHEgHjgCHnT/TgCbX0iNCiGAB0wAygWuAb4AUtb8AYU9e\ntwVYQhsRiQQqAxc6Tw2gEsb84oDtwG/An8Bfzn+vAw4AKc5zCDiKMbBjwDFVPSkigjHCCMy/5Sgg\nAYh3nvpAvTw0ic1D2ilJeYDZQDFH3y/AT86zyXl+VNW/32AJbawJWryGY3h1gEuBes6f1YFfgY3A\nj8BCYAvG7HblZiXmvPeE82QCqRjzPKXnF2D8ceKIYvepLx8DLgP2YlaYFYFqznM98AhQRUS2AV9j\nVoxfAF9YYwxNrAlacoyI5AcaA02ByzGmtxVYD3wJzAK+c8s8VHWPiNydScycKI6cMsq7VXWP85KD\nwFfO8zciEoVZtV7sPJ2AWiKyAfgMWA2sUNV9fvpRLD7EnglasoyI5MWs8Fo5Tx2M4a0GVgFrVfWQ\newrPTAtZdvICNn33Cr1anmaA2UJEYjA/e2PgCozp/wx8BCwDVqpqurc0W/yHXQlazomIJGACCe2B\ntpjzuyVAMrA6GLaICpKH+P05NUAA5+dc5TzjnK1/A+BqYBhmpfgR8D7wgaruPutgloDCXpGx/AcR\nyS8i/xORD4A/gG6Y1V5tVa2jqgNUdWkwGKBHPHKSCDKJOuHNcVU1U1VXq2qyqjbBBHvmYT4oNovI\nShF5QESKe3Nei/exK0ELACISjfkF7gq0BD4GXgG6qmqKm9pyw3Hylo3gOEeJK+LLeVR1L/AS8JJz\npngN0BkYJSLrgdeAt1T1sC91WLKPXQmGOSJSW0SmAL8DDwKLgfKq2l5V5wSzAQJkElktLydOHqKA\n37JFVDVDVeer6q3ABcAM4AbgNxF5XkSucK73WAIAa4JhiIjEiMjdIrIOc2UlBaivqs1VdbaqHjjP\nEEFDBlGVhJOZe0ks5Mb8qpqmqm+r6g2YazgbgKeBn0TkIREp6IYuy/9jTTCMEJEyIjIak0bWARPc\nKK+qSar6i7vqfEMGUeWPkS/lFyrEDpPkKDe1qOouVZ0I1ALuBhoBv4jIDBGp5aa2cMaaYBggInVF\nZA7wLSb9rImqXqeqC1XVqwGDQCOdmNKRHNtfnF3HjhLXxm09YC55O0GVLpj7iDuBpSKyUESutFtl\n/2JNMIQRkaYishD4AHN5ubyq9lHVLS5L8xupxJSIJXVvBX75Yx+J7d3W829U9U9VHY5JI3wXc8F8\nrYjcmJUCEZbcY/+SQxDn4P0T4AXMtY1KqjohEC8y+5pUYovGkrqrHDuWb6d8c7f1nA1VTVfVWZjc\n6rHAEOArEWlvV4a+xZpgCCEiDURkCcb8ngOqqerMcM5kOEJ8oWjSfy/OrunfUqf0MEnO57amc6Gq\nJ1T1XUx2ShLmIvaXInKdNUPfYE0wBBCRKiIyF3jHeaqr6ouqetxlaa5zmISEaNJ3jNHBX5bgr/T9\nFO7ltqas4JwbzgcuAUYA44CPRaS+u8pCD2uCQYyIFBaRScAaTPmpKqr6jKpmuiwtYEihQEwUGVsB\nLubrZZuper/bmrKDY4bvYvK0XwXeE5E5IlLBZWkhgzXBIEREIkTkQUztuyjgQlUdF87b3rOxn8L5\nIsn8CaAsvw5dR8NKQ2VE0KWyqepx58ywKub/+3oR8TiFHSy5wJpgkCEiDTH17W4GmqtqT5usf2aG\nSXLMYRIkH8e3AozRwd9exHd/7qDcU25ryymqekRVkzElvmoCP4hIO3temHOsCQYJIlJIRGZgrlFM\nBK5S1Q0uywpoMomsWogDJz067O+7kJfwlWcFzdq7fXE6t6jqr6raEbgPeByzTS7tsqygxJpgECAi\n7YDvMRWUL1TVV21vjPOTSWT1wuz/x/noZO0zqwh7U/+g1BS3dHkTVV2KOS/8CvhaRHrY+4XZw/5l\nBTAiUkREXsOs/Lqq6gOqetBtXcFCBlEVCpCS+u+vN2XVwEW06ZYkw13JJ/Y2TsEGD9AMuBMTRa7k\npqZgwppggCIi12FWfzuBOqq60mVJQUc60WXjOfKfC+JT9KEZVfj5r81Ufd0NXb5CVX8AmgDvYbJO\nutmzwvNjTTDAEJFYEZkOPAl0VNX+qvqf1Yzl/KQRUyqeI2esiNOQtbcupWXLR+TxVv7W5Uucy9ZP\nYFaFDwDzRKSYbUB/dqwJBhAicjEmxzc/UFdVV7ssKahJI6ZYLKlnjJyP04GfXMsH7yyl5VuBnkWS\nE5xVYUPgB8yVmt+wDejPiDXBAEAM3TG9O0aq6m3BXsw0EDhKXGI06X+e7fsV+eWWCE7qJqov8qcu\nf+Fcmp8ExFalZlR7ZhUAYoDn7Irw/7Em6DIiEocpy94LuFxVX3VZUshwhPgC0aT/drbve3TYiZYs\nuWolV1z1gEwb6k9tfqR8FNHp+XkOJeLU144B5d2TFFhYE3QREakKfI65+tJQVX9yWVJIkUKB+GjS\nt53rNWN08NedeLP3W3T0DJDx1/tLm7+IJW57a57LfxLhfe499eUEoLCLsgIKa4IuISKtMO0bpwJ3\n2eCH9zlIwagoMn8+3+um6EPT2jF/1mt0fXeQjG3sD23+4maeWvcZzWUft6Sf5EQKkIa5cvWyiNzl\nsryAwJqgn3HO/x4CXsREf5+xF5+9j0c8so/EiHxO3vD5mK3dejRn+Xtv0mnFYBlTz9f6fI1HPHK7\nvPjpUlpW6MSbU3awrSymi2A5VR2AaSA/SESmikhYd50U+/vnP0QkHzANE7W7QVW3u6sodBkqIy6Y\nRN+dRzQ+W/fkusprC1Zz+TW38Hqbx/WRZb7S50s84pGvqfP5bkrUu4A/N8zVDrXP9DoRKQC8CRwH\nOqvqEb8KDRDsStBPiEg8MB8ojenxsd1dRaGNkzKX7XqKr2nX667mo1df4o4PH5IpPXyhzZcMlRGJ\nq2i6TYm4OJXY1FpsaHS21zo3EK7DXMhfJSKl/CY0gLAm6AdEpATwCaa3b7tw/cT1JxlEVS7IwRyV\nFnte77qzA3OT59Bl2h3y0hse8QRF1sUgGdv4HW76LZE9cWu4LE9bFl4zXJOOnus9qnoM6A68AawJ\nx6531gR9jJPD+Rmm+kt3W+3ZP6QTXbEAKTn+sJmuPYffxfMtPqdBuxVc9dujMvpCb+rzJh7xyD3y\n7DOz6ba6AetWr6VxoZt4Z8IYHZyly/ZO4daxwGDgIxE56+oxFLEm6ENEpCZmBThWVUfaAIj/SCe6\ndDxHcnXhfJwO/LgTb5YozL7fZtPt+3tl1pOBtiocJGMbL6f5H6u5/H938ELX7VSoUYdvN87QHo9k\ndyznjurdwPsi0sL7agMTa4I+QkTqAR8BA1R1ptt6wo00p9VmbscZrkkpc7XDZbfy6r2f0uSeebRP\n6SOT7z3/O33LUBlR/GZ5e+VMuq8uwa7NnXmj6H6KXr+bYokX83XTnI6rqh8AHYDXRORG7ykOXKwJ\n+gARuQxYiNn+znFbTzhyqtWmt8abpH2f68RbCZfy5etv0unpxvLZ7t4ytY+/V4aPyajKneTN5TPo\n8edBCpbuzsxGb2rHZkeJa72Aa7tcx4L2wzUpVytgVV0FtAZmiMhN3lEewKiqfbz4YK6/7AZau60l\nnJ+r+GhHN2Y+6Yuxh5Jc4E6ee6ECWzMvZMPR23nxjUcZWdlXP8swhkkvpvRuxeItBdl/8hoWbRrA\nuFanvj+E4cUr8nNGN2ZO9+a8mBL+fwE3uf3/05ePvSfoRZx2iB9gMkA+cFtPOFNfvtjfiLWjn9Re\nE3w1h0eSI/ZQZNAmqt+7nnrlavLDnkpsXVmCv2bHcfRDj3py/Ms1VEYU30fi/b9S9qZvqFsjltTj\n9fliSTl2DB6tj/5w+mtby+LNisiHek2V3P9U/0RE6gKLgQdU9R1vjx8IWBP0Es4/lg+Bbqr6vtt6\nwp0qsiX9Bt7rOkH7z/XHfENkZKk9FH14B+Wu/Y6LKh8jX55KbN1Xkp1bCpCyOZbUzTGkbcjL8b15\nOHkwDyePniCi8Akiih4jX+mjxNU+QnzVfSRW/YUKZX6ndFQNNqZU46dPS7JzVjxH3juTqfaQGZM+\n5JoHbuOVMiN0qNe2/6dz2r/t/6nqYl/M4SbWBL2AiFTBRIF7q+rbbuuxQBHZe7I7M2uN1kd/9Pfc\nHvFIBlF1UyjQfh+JTQ5SsOxBCibuIzE+jZiIDKLyHCOfRJN+MpbUE3EcPVaEfQcKcuCvghzcVJj9\ni+I5Mu98d/weldGXzabbp7fxyv+e0Idf9uXPJCKNMRWrb9RQq3Pp9n482B+gJLANEwRxXY99lCQ8\nUXnJ1CQ8+dzW4sOfMd/FfJlyC68t9NecmNzj3cDFbv/83nxsdDgXiEhBzDZhltprMAHDMfJVLshB\nTdZhx9zW4it+otqiCE6cqMZmv5X/UtPZ7n5ggYiU99e8vsaaYA4RkUjgHWA5MNZlOZbTyCSyWiEO\nZJ7/lcFJH5l873KaX9WCZS1P76nsD9QER8YBC0UkJLr1WRPMAU4HrxnAEeBhdfYKlsAgg6iKZ2q1\nGQoMkZFl5tJhegfmPj5GB3/phgZVnYrZAc0VkaBuYg/WBHPKYEzD666q6tdPYsv5SSe6XAKHD7ut\nwxespdEntfl+8wztMdhlKf2B/cDMYG/raU0wmzipRD2A61X1nNE7izukEVM6niP73dbhbbrJ7Ok7\nKFfqEr66wm0tzof/7UBtoI/LcnKFNcFsICIXAjOBDqq60209ljPjtNrc47YObzJQxl05j/Y92jH/\njhE6dJ/begDUtIRoDzwiIi3d1pNTrAlmEScSPA9TEGG923osZ+cocYkxpP3ltg5vMUySYxbRZkEb\nFr03Ufu94bae01HVX4HOwCsiUtltPTnBmmAWEJE8wMvAh6r6gstyLOfhCPEFz9VqM9j4gZofxnE0\nvRLbbnZby5lQ1ZXAcOAtEYl2W092sSaYNfoBicDDbguxnJ9D5I+LImO72zq8QR+Z3HMVTZu0YNlV\n/r4Ok02mAz8Bk90Wkl2sCZ4HpyxWf6CLmlLklgDnIAWjo8g4b6vNQGeIjCz1LjdOac+8USN06Aa3\n9ZwL55pYd6C5iHR1W092sCZ4DpzLoHMwKXE73NZjyRr7SMwbSeYmt3XklvXU+7gGG39+Ru9LcltL\nVlDVQ0BHYIqTTx8UWBM8C87dp5nAe6r6ntt6LFljqIwoqggRnAjqwEhPmT52E9XL1+eLZm5ryQ6q\n+i3mfPAVp8VswGNN8OzcClwIDHRbiCXrZBJZPZF9J3JTy89tHpXRtd/m5gE38N6DviqP5WOmAQeA\nx9wWkhWsCZ4BESkDPAHcpqo5attocYdMIqvktNVmIOARj3zClcuasuqzKfpQUBblUNWTmIZNPZwz\n9YDGmuC/cK7DPA9MUdWv3dZjyR7pRJfPz6GgzeT5hQovpVAgvhYbWrmtJTc4yQQPAC+ISIzbes6F\nNcH/cg+QH1MpwxJkpBNdJretNt2iv0y4+gOuvbUtCzsn67A0t/XkFqfizHfAMLe1nAtrgqchIhcA\nozEl8m2T9CAknegL4jgaEGll2WGYJEctodW8tiyc+7g+ssBtPV7kQeAuEbnUbSFnw5rgP3kSeEZV\nv3NbiCVnHCWuuDdbbfqLjdT4IIa0zIr80tltLd5EVXdh7tk+F6jRYmuCDiJyA6Yixki3tVhyzlHi\nCkWT/ofbOrJDH5l8xwqaNW/BsmsCPCskp7wC7AJ6uS3kTFgTBEQkFpgC3G+jwcHNYRISokn/1W0d\nWWWgjCs8n3azOjB36ih9LCQLczjZJA8Cj4pISbf1/BtrgoZHgM9VdbnbQiy5I4UCsVFkbHNbR1bZ\nQK2PK7F15wztEdQ1+c6Hqm4GngF81gc6p4S9CYpIBcwyvb/bWiy55wCFIiPJ3Oy2jqzQW6YO+JqL\nazbh02Zua/ETo4DGItLMbSGnE/YmiPlkmuTURbMEMR5JjjhAoTzBYIJOcYQx7Zk3zKOesMhLd4qw\nPgI84dzHDQgCRogbiEgToD4w0W0tltxzjLxVEjiswXDH7ksu/agWG7ZM156j3NbiZ94CMoDb3BZy\nirA1QadAwuNAkqoG/C+N5fxkEFW1MPsDvtzZg/LUkA3UqlyHb692W4u/cYIk/YBRTkDSdcLWBIEb\ngARMxWhLCJBBVKWCpAT0B5pHPOXm0sFzI+8OHKuDwrJPjap+BqwFHnJbC4SpCYpIBDAGGGhbZoYO\nGUSVS+BQQLfa/JQmH13Klz9O1d7hfgQzBOgrIgXcFhKWJohpDLMfWOy2EIv3SCe6VBxHD7it42zc\nL0+P3EKVcnX5Jmg7s3kLVf0JWEQAtOsMOxN0VoFJgMc5n7CECKnEFg/UVpuPyajK82g/uD3z+gZp\njUBfMBzoJSKF3RQRdiaIWQXuBZa5LcTiXVKJLRIboK0213DZsoas+3ay9nnKbS2BgqpuxbSxdXU1\nGFYm6NxNGgok21Vg6HGYhELRpAVcq83uMnP8dsqXrMf6Fm5rCUDGAveLSIJbAsLKBIHrgKPYVWBI\ncpiEgGu1OVDG1ZhH+4dv5N0HhujI/W7rCTRU9WdgOXCvWxrCzQQHABPsKjA0CbRWmyLIGi5b2pRV\n6ydqv1lu6wlUOsKsETBsgki8G/OHjQmKSCOgNPC221osvmE/hfNFkvmT2zpOcS8zJ/3JBUVr8kNQ\nl8r3GSLVEZn5JrxZBQ49C3e4ISNsTBBzS32SrRgdmiTJ8EIZRJGX47+7rQUgSYbXmEf7XtfyQY/h\nmhSU5f59goggcgUi84FPgD+AarfAfZvgHieTy7+SwmFnKCKlgO+B8k6DaEuIMUjGNnqJO1bv1JJ5\n3dYCcLV8tCOeI7ve0xsauK0lIBDJC9yIqdZUCNPN8UWclFUnaPkT8D8no8RvhMtK8F5gjjXA0CWD\nqGqFOJDhtg6AnjJ95DYqXlCdTa3d1uI6InGI9AI2Y9LkxgA1UJ3BaTn7TpvOabhQfTrkTdDpa9Ad\nmO62FovvyCCqQgFSXG+1OURGlptH+0HXsaDfOB0YvtFgkRKIjAS2A82AW1G9HNV5nD1V9XngGqfh\nmd8IeRPEFErYoqo/uC3E4jucVpuur/TXU29JXb7Z9KT2etJtLa4gUgOR2cBGzLb3MlRvQnXN+d6q\nqinAXOB2H6v8B+FggncB9npCiJNGjOutNh+UpwZuoFalS/jqGjd1+B0T7LgSkfeBFcCvQBVUH8Dc\nA8wOz+HnAElIm6CzrG4MvOu2FotvSSOmWCypu92af6iMKD6fdiNv4L1hI3VIUHW7yzEieRHpBKwD\nZgLvA+VRHY7q3hyOugY4CTTxksrzEtImiKle+46qun5WZPEtR4gvHE26a/X5vqHuh9XZtGOaPhD6\nlaJF4hHpDWzBBDJGYYIdM8llgWInkeFZ4O7cC80aIWuCznL6f8ALLkux+IEjxOd3q9Vmb5l6/xfU\nr30Za0J7GyxyASKjgF+AK4AuqDZF9T1MdNdbvAa0F5FoL455VkLWBIGamMrRn7otxOJ7UigQF036\nVn/P+4g8XnAB102+kXefSNZhfp/fL4hciMizwI9AAaARqjejutYX06nqTuAboI0vxv83oWyCnYC3\nbJ5weOC02tzi73m3UGVhOXbseVrvH+DvuX2KCXY0Q2QBpsDBdkyw40FMCSxf8zpwix/mCU0TdLbC\nnYA33dZi8T0eSY7YR2KefBzb5M95+8qkWz+lSaPLWR06l6JNsKMz8DkwA3gPqIDqiFwEO3LCXKC1\n+KGoQkiaIFAbiAa+cFuIxfccI2+5OI7qcE3yWwBsovSLXkzrWe2YP3uEDt3gr3l9hgl2PAT8DDwA\njAAuRHVWboMdOUGN4a4FfP4BE6om2B6Ya7fC4UEmkdUKs9+vhTG+oe47CRxOLc3v9/lzXq9jgh2j\nMdvdy4HOqF6B6nwvBztywnygna8nCVUTvA5zZ8kSBphWmwf9tlp5TEZduZC2bZqy6iaPeoLzg1ak\nJiLPAT9gAogNUe2I6jqXlZ3OfKCtmOILPiPkTFBESgCVgdVua7H4hwyiKiRw+Ig/5vKIR5bT/J2W\nLF00Uft94o85vYYJdlyFyAfAR8A2TLCjl5+CHdlCVX8DduDji9MhZ4LAtcASVT3mthCLf0gjpnQ8\nR/zSanMnJZ86QKG4C/jzJn/M5xVMsKMLsB54GpNBVR7Vkai6mmqYBRZgfqd9RiiaYBtgodsiLP4j\njZjicaT6PHKZJMOrvMcNPVqzuMck7Zvu6/lyjUgCIn0wwY4egAcT7JiNauDrNywBfNqnOaRM0CnM\n2Ayz1LeECanEFokh1eetNtfR8IOGrPtusvZ50ddz5QqRkoiMwWR2XAZ0QvVKVN8PgGBHdvkcqCgi\nRX01QUiZIFAH2Kuq4ZHAbgHgCPGFokn3aVn9XvLkwxuoVbEO37b15Ty5QqQWIs8DG4A4oAGqnVH9\n3GVlOcY51voEuNpXc4SaCTbHrgLDjsMkxPuy1aaTGjf2Bt4bM0KH/umreXKECXY0R2QRsBSz9a2M\nam9Ut7mszlssxYcmGBD9GLxIc0x1WksYcZCCMVFk+Cy6uZVK75dn++7p2nOor+bINqZiekdMz44Y\nYAJwYxCd9WWHVUBPXw0eMitB5zzwMuzVmLBjP4n5osjY6Iux+8nEGz7hyiZNWXW9L8bPNibY0Rez\n4usOJAE1UX02RA0QTJO0kiKS6IvBQ8YEgapAiqr6/IDcEjgkyfCEVGLIy/Ed3h7bI8kRy2jxcjvm\nvz5ck7729vjZQqQUIuMwmR0NgZtRbYbqgiAMdmQLNT1J1mEKJHudUDLBRphcQ0sYcYx81RLZd9IX\nmRtbqfhyBCcoy29+7XnxD0RqI/ICZjUUDdRD9RZUwy0vfjUmrc/rhJIJXoYpzW0JIzKIqlqQg15v\ntTlERl60iDa3tGTpHR4ddrbuaL7BBDuuRmQx5p7cZkyw4yFUf/GrlsDhc6CeLwYOpcDIpUBg39+y\neB2n1Waqt8ddR8P5TVm1dpwOnOftsc+KCXZ0wgQ7ojDBjhtQDYh+yi7zNXCJiIi3C6OExErQSbC+\nEPjObS0W/5JOdFlvt9rsI5MHbKJ6qaps9nkFEwBE8iPyMLAV6AY8BtRC9TlrgAZV3QUcBcp7e+yQ\nMEGgGvC7qvolid4SOKQRUyqeI15rcp4kwxPe5/pRbVg0dpwO9G0qngl2PI7J7KgPdED1KlQXhnqw\nI4d8BVzi7UFDxQTrYnoSWMKMVGKLxZK6y1vjbabq3Av4c99M7e67O4EiFyHyIibYEQlcimoXVNf7\nbM7Q4BtMVphXCRUTrI35B2UJM1KJLRxDmleyOAbKuCuX0vLqy1hzszfG+wcm2NESkQ+BxcBGoBKq\nfVDd7vX5QpONQA1vDxoqgZGqwBy3RVj8z2ESCkSRketWmx7xyCpavdWKJR+O1wHe61Bogh2dMcGO\nfJhgRzt71pcjNgHVvT1oqJhgNcw1AkuYcYj8sdGk5/rayJ9cMHYPRfO3YJl3VoEi+YF7gT6YJuWP\nAovtWV+u+AmoLCJ5VdVr7RSCfjssIhFARcw/NEuYcYBCUblttTlURhRfwHX92rJwYK6bNYmUQWQ8\nJthxKdAe1eY22JF7VDUV2IWXI8RBb4JAWUz5LK/fFbMENh7xyD4SIyLJzFXe8Hdc9P6F/Lhjij40\nJceDiNRB5GXgW8wO61JUu6L6ZW60Wf7DL1gT/A/lMH8xljDjOHnLRpPOcE1KyekYD8mUGz+jcb0m\nfJr9O4Em2NEKkSWYauYbMMGOvjbY4TO242UTDIUzwTLAb26LsPgfp9XmMSiQo/d3lLfybKbZC9ex\n4A2Pen7I8htFIvn/YEcEJtgxxwY7/MIOzMLHa4TCStCaYJiSQVSlAqTkuHxUHEdnKJK3EWuzViBB\npAAiAzBd2u4EBgK1UX3BGqDf2I7dDv+HskCur0hYgo8Mosrn51COsoSGyMjyC2l7TyuW9LhPnzl3\npNEEOyZgzK8u5orL1aguxst5rJbz8gdQ0psDhoIJlgACq+S5xS+km1abB3Py3m+pM78e6zdN0P4v\nn/VFJrPjVLAjArgE1VtR/Spnii1eYBdQ3JsDhoIJFgX2uC3C4n9SiSkRm4NWm31lUqd1NKx5Ed/9\nNxhigh0tTsvs+IH/D3Z4vXCrJdvsBop5c8BQCIwUwZpgWJJKbLHC7M/W9RiPeGQFN8xqy8I5Y3XQ\n//cl+WfPjmhsZkegshco7M0L06Figj5vvG0JPI4QX7AkO7PVXvVXyk47Rr685dhxFwAi8ZjyVX0x\nh+5JgL3YHKCo6nEROQgUxqwKc01Qb4edbJECQI7OhSzBzWESEqJJz9IWVUSKVpNLWi+k7X2tWfxQ\nMp5EREZjjK8J0NFpUB7yPTtCgENAfm8NFtQmCMQC6U4jFkuYkUKBLLXaFJFbgB3FuGNBS56V2xnQ\nFVORpADQENWOwdygPAw5BCR4a7Bg3w7HYarNWsKQ/RTOF0nmT+d6jYgUBXmuKZ1jmpFCL5J5Gi5/\nD2onqZ7zvZaAxasrwWA3wXjAVpMOQ4ZJcsxhHpV8HD/rSrCXPBl5M6MnbaBlzF8k8CezKYeSBqlA\n/iQ/6rV4lcNYE/wbuxIMUzKJrFqIAyfP1Amuj0xO3EPR6au4oUMR9qTXZ/jxQbyfty7KMfOSSMxZ\noCU4ScMEOBddAAAgAElEQVRE8L1CsJtgJGCvMIQhmURWL8z+zNOvjPWXCdW3UfGZldzWtBYb/urM\nGw+M55GXDsH2TlDwGKRj/s3crar2WlXwkoHpxucVgj0wkgewkbww5PRWmw/LEy3ayKIfZtPtxwyi\nit7LrLYrtFnJ8TpgJvBwflj7ockxbwmUU9XXXRVvyS3HMFW6vUKwrwStCYYp6USXLcOvJ6+QjN9/\n5I6SzVn+aR8m1zbVYK41LxIpAzwMNHBWfnb1FxqUwYupc8FughFYEwwrPJIc8Sclkr+i/r37SMx7\nFR+/05KlPYfqiDNdnH0cmI7qNn/rtPiU5pijjbHeGCzYTVDcFmDxDwNkfP7dFJu6mtu6xJJ6/BK+\n3HOMvBte0DvP3BNE5EqgMXCPX4Va/MHzwGpvDRbsJpiJF88GLIHHIBlb4RcqzFjBHS0q8/O+dswf\nXICUSWu47KezpsyZTKLJwABs24VQxKs7wGA3Qa9GiSyBwwAZ32QT1Z9azX11GvD5tjt5oeM4HTjX\nZLhBY/msYAxpZ6sjeSfm/uhbfpJr8S8xmGsyXsGaoCWg6CuTOn/HReO+5p6yzVix/kGeajBCh66H\na/7xukPkj4si47+9ZUyryxHA9bbgachiTfA0rAmGAB7xyB6KDvqC+v1/p3PB5ixfdBlrmo7UIb9B\nhzO+5yAFo6LI/PkM3xoMLLFd3kIaa4KncQSTOmcJQoZJcszvlH7iU265UxG5gpUvtWZx3/P1/jWt\nNgdG5Pt33rBIBUzD84t8KNviPtYETyMFKCgionbrEzQMkrElf6PM0x9z77Wl+f1QS5aOTmTfSI96\nsvT/8AQRJSI4wQgd+u97f+OAyaju9L5qSwCRH1NEwSsEtQmqaoaIHMOU1LI5xAHOYBlz8SaqP/0J\n3RvU5ZvfujDnzona7xVokK1xnJS54//YBIg0BRoC//OqaEsgUhjY763Bgj1tDkxB1Zw1nrX4hX4y\n8fpWsmTzDHp8eZy80fcy68rl2rycMcDsk0FU5YIc/P9WmyJ5gEnAIFS9tk2yBCxeNcGgXgk6HAAS\nAbsFCiA84pH9FO71JZc+toXbi1zFx8t7Mr3tKH3sZ7g+V2OnE12xACmnl1C7HZNPanOCQxwRicEk\nSdgzwdP4C9N283u3hVhgijwU+SMXjv2Mm7ofJS7yCla+0ZaFvR7TUV5rgZBOdOl4jqQAp3qEjAY6\n2CsxYUEisM+bMYBQMME/gQvcFhHuDJURiTsoN/0THu5QiANpTfj0yRLsGnKmen+5Je2frTYHAstR\nXefteSwBSSm8vOsLBRP8C2uCrpEkw2tspMaMFfRoWp1Nf93Aew9O1d7PQF2fzZlKbNGCHNyKSDmg\nJ1DHZ5NZAo1SQLY6DJ6PUDDBP4GybosINwbJ2Ku/p/bUT+ld4zLWbLqTF9qO1wGLoanP5z5KXKES\n/LUTU0XkSVR/9/mklkDB6yYYCtHh7UA5t0WEC/1k4j3NZMUfs+m2NJr0/T2ZXnuRtrnQGKB/OExC\nQjNWKHA5MN5f81oCgtLYleB/2AZUcFtEKOOR5Ig9FEleS6Neu+kU24wV7zVkXc9xOtArza+zy0EK\nxLZjfgtgKKr2fmh4UQGY680BQ8UEK9qsEe8zRgbn/5nKUz/h1i6RZB5vwqezS7LzkWQd5mpfl8as\niYwmPQJ42U0dFleoDJwpZzzHSCj4hojsA6rb5jneYZgkV9hM1Rkfc1WL8mzfdwlfjS3G7klZTWvz\nJR/ItXFV2HykHL9eG6UZC93WY/EfIiKYVNlyqnrAW+OGwkoQzCdDNWwPiVwxREY2/p7a01bSu86l\nfLmtC3M6TtK+c6GR29L+5mK+HrSZqlTVLdYAw4+iQKY3DRBCIzAC8ANwodsigpUBMv6WFrJs+3R6\nrgaO3cOzDZZpi8rGAAMIkfxF2fPgJPpkui3F4grVgS3eHjRUVoI/ADXdFhFMeMQjByg0aB0N++/g\ntgJX8smi7sy8fKwOCuTrJo/8QcmtOylV0W0hFle4CB9khoWSCbZxW0Qw8Lg8EvszlSeuptOdx8kr\nTfj0pVYsOW8NP9cRKQXc/yJ3vp/A4aJuy7G4Qm3gO28PGiomuAGobSPEZ2e0PFpqIzWmf0zva4uz\n63BTVo3eT+FRz+tdwdKydAQwawtVqsVzxGsVRCxBRW3gNW8PGiomeOryZCkgkLdzfidZhl3yPbWn\nf0y/BrXY8FsH5t45VXu/Ape6LS3riNTGdFSvmkbMwlhSbQAszBBTLq0WPtgOh0RgxFn9rQfqu60l\nUHhMRrVrI4u2TKLv+jRiou/gpWaf6JXljAEGHY8Do1BNOUpcYgxpf7ktyOJ3qgN7VNXru4BQWQkC\nfAHUA951W4hbeMQjR4jvvY6Gj26ie5FmrFjejdltJmj/n81CKggRaQFUAW4AOEJ8wWjSf3NXlMUF\nGgI+qRQUSia4Hujttgg3mCIPRf5EtTGfcuN9KRSIbMqqNy5jTa/H9RGv1fBzBbMFGg8MRjUT/m61\nud1VXRY3aAB87ouBQ8kE1wKvikheVT3uthh/8Lg8kriRGtNX0KdDPEfSGrH2SUGHztTuofLz3wqk\nA2+f+sJBCkZHkeHVtClLUNAQeMkXA4eMCarqXhH5HVNbLqR7zibJ8Go/Ue2Zj+l/RRW27GrLwl7T\neeCZ77V26ETGRaIwEeHbT68YvY/EvJFk/uieMIu/EZECQFXga1+MHxKBkdNYhT8K2rnEIBnb/DpZ\n8P1Uem88SMHiXXmt7afa5IJp+sAMVULHAA33ARtQXXXqC0NlRFFFiOCEK9VrLK5xOfC5qqaf95U5\nINRMcCVwhdsivE1fmXRnM1nx2yzuXRZJZko3Ztf5UK+pMVn7+K2Gn18RSQAeBR47/cuZRFZPZN+J\nQCjkYPErVwIrfDV4yGyHHVYAT4pIhKp6vbeFP/FIcsRuig79nIZ9/qJjfDNWzG/E2vpjdVA4XA/p\nC3yE6renfzGTyCqm1WYZl2RZXKIZ0N9Xg4eUCarqThHZibkvuNZtPTmhv0yI30fi5NV0vS0vx7UJ\nnz5Xij/6J+uw8OinK1IEeIgzdGRPJ7p8fg4FdnqfxauISGHMHUGfNdIKKRN0+BC4hiAzwX4yscwf\nlJrxCbe2LsuvB6/hw2FF2Pt4GG79BgOvo7r1399IJ7rM3602LeFCK+ATVfVZId9QNcFk5wl4+suE\nS7dQ5emV3F2vHut33MqrXSdo/zfMjYAwQ6QscCcmPeo/pBFTMo6j+/yqyeI2bQGf1o4MtcAImAhx\nLTHbqoClr0xq11KWbnmWe74QNO8DTLt8qbasYAwwbBkGzED1zzN9M5XYYrGk7vKzJotLOPnCrYFF\nvpwn5FaCqpouIsuA64AXXJbzDzzikT0UffArLhmyja5Fruajjx7miWuG6oht0N5tee4iUgNoh0mR\nOyNHiStUjN2f+U+UxWXqA3tVdbsvJwk5E3R4F7gZF01QRIoC5YHtSXgO/kGp0eu46f6jxEU2Y8Ub\n17Gg12M6KrjT2rzLCGA8qmf9OzlMQkI06b/6UZPFXW7Cy53lzkSomuACYJqIxKkLLRlF5BbguWhi\nM+txV/yL3EEhDqY1ZdVTdfh2yH36TKiktXkHkfqYRiZ3nOtlKRSIjSJjm39EWdzEaap0M9DB53OF\nag1SEVkKzFTVt3w0QQKm6XtZ5ykGJB6BC9bBTak0yNOfF50XnwDyZp4kjwqq0aQfjyX1eD6OZeTj\nWGosqSn5ObQnniN/5OPYZkE3FOLAZyN0aHgEAczxxZuozjzXy4rJ7hP38OwlY3Twt+d6nSX4EZGL\nMTnjlX1dKDlUV4IArwNdgNyZoEgh4GJMxLKm89QAYoAdpz27gJQMaNAQ5AM2UYp+/MFvnECOVqbh\n0JpU25aX41HpRMdnEFVQkUSgSCqxxQ+TUGwPRS8/TEK7fSTG/07pyOfl95Ml2ZlShL2/F2HvhpLs\n/LgkO9/prVNDp7KyKZVVFnj+XC/zSHLEAR7NE0nmZv8Is7hMR+Adf1SKD+WVYEGMOZXPcos+swSv\njsk/buQ8ZYBvMRVtf3CeH4Hdfyf2m0+tAcA1B+DV2nDvHxB92shpmF6pWa6I3F1m5o0io/5J8lx1\nhPj6uyh+4Q7Kld5Budhy7Egry6/by7FjZWl+fy5Jh/ukxJDPMX/f64CJqJ4zKv6YjKo+gx4/7tPE\nULzRYDkNJyr8C9BO/5U15AtCdiWoqgedKPFNwOyzvtAEMFqe9hwHPgHWAFMxSfz/PcMTEURaAo9g\nVoaTgfsLqab8IfIZ8ByQCUQCd2e3MbxTDmuN8/xNkgwvkEZMh90Uu+4nql37Djd1f0G2Ha/K5i0V\n+GVBZX6e1E8nBktq3fVAFFlYrWcQVbUQB45Bou9VWdymKZDiDwOEEF4JAohIB6C3qjb71zfKADdi\nDl3rYnKOlwJLgJ8511+KSD7MUn0AkA+YALx2qujnaXP/HR3OrgFmB48kRxwl9sZdFL99E9WbbKJ6\n4kV8t7sSWxcVZc/jE7R/YJadMp/2XwHDUH3vfC/vJU/2XUPjYev10oK+F2dxExF5FtioqhP8Ml+I\nm2AU8BvQWGE/cAsmAlkZeB8Tfl9KVkr0iMQD9wAPY5bq44FFqAZUt7Z+MrFMCgUGbKNi+6+4pExN\nfthXiw2vF2P3sIAKtIjcDAwEGpzzQ8ehu8ycvIUqN32sV9nqCSGMiMRimqXVUtWd/pgzZLfDAKqa\n0Upk4UhzZaYEsBiTTrcM1WNZGkSkONALU99uBdAR1YA9g5uo/X7DtBno3UcmJx4iv+db6tyyieo9\n18mSLXX5ZszjDHwxK8bjM0QiMP8f+mdVRzrRpeI4mrWzXUsw0wlY4y8DhBBfCQKME2n7O8zpApUa\nq+7N8htFqgL9MP9T5gBPoBq0Zd2HyMi626g4ehVNWxbiQEZ9vni5FH88MlyTDvtdjEhX4EGgSVZN\n8GZ5e2UeTh57Uztd7VtxFjcRkbXAKFV9329zhroJAojIKmCyqr6ThRc3wgQ7LgdmAE+hGjKVjIdJ\ncr6/KDH0G+r23EbFQi1Zurwqm+/xqMc/mRgieTHR9ftR/Sirb2sjizYWZ/dXL+j/bvWdOIubiEhd\nYD5QwZ/1QMPlusHTmJXHmRHJg8j1GLOcA3wMVEA1KZQMECBZhx17Ru9LWqcNi9zJC+13U6zik/Ta\n3lHeWvWYjKrqBwm3AX8Cy7PzJtNqM8222gxt7gNm+bsgcriY4NtAVRGp84+vikQhcjewAfAA04Aq\nqD6JC+l2/ma8Dnh/mbao1I3ZVx0lLnE6PTd1lLdWJcuw0j6Z0ETWk4Ch2T2TPEyCbbUZwjjFUztz\nrutsPiIsTFDN9ZVpmIrFIFIQkYGYKG8nTOCjHqqvn/FOYIgzTgd+slDbXngfzzROocAF03hgxx3y\n0ltJMjzOy1PdBWxFdWV232hbbYY89wHz9Sxl1HyKqobFAxQpAwf3wzSF/QqvKNRxW1cgPg8zoWND\n1uwpz7bMHkz3eGVciFL4VaFRTt6fQIo+xogybv/d2Mf7Dyah4A/gIjfmD4uVIJi+xE/A9jVQD7gY\n1dvw0430YGOi9ntrrTYq2paFnsW0fvQKWblzkIxtlMthuwHfo5rttgdJMrxQBlHk5fjvudRgCUxu\nAX5U1e/cmDwsosOnEJGKwOdAJVW1vSqyQJIML7CRGu8up3mz1iz+oDI/35ysw7LX70EkBvgZaIfq\nl9nVMEjGNnqJO1bv1JIhfa81HHHyhDcAfVR1iRsawmYlCKCq2zClunu6rSVYGK5JKW9px+Z381yL\nLVRp8i437h0o467J5jA9gM9zYoAAGURVK8x+nzXasbjKTcBhTNqqK4SVCTqMAfqIiLcP/UOa8Tpg\neVsWFb2YrxfOptui/8mLb3skOeK8bzRpUI9g+ofkiAyiKthWm6GHswocAoxQF7ekYWeCqvojsBJ4\nwG0twYZHh514Uf/X+S6ev+Yb6rZaSJvd/WVCzfO87T7gM3Jx3uO02jyU0/dbApbrMRWHP3BTRNiZ\noEMSMEBECrgtJBiZoP2Xtmde0XLs+OFVbv3uIZny0BlfaM4CB2D6h+SYNGIusK02Qwsx+eMjgWQ3\nV4EQpiaoqhsxnz793NYSrCTrsIy3tOMVN/Ju0hy6PNFVXvvwDNvjbsAXqH6Tm7nSiCkWS2pIZe5Y\n6Io5C5zvtpCwNEGHZOABESnmtpBgZrr2HHUHLzX8kQsbL+eqHUNlRFEARKIxpbJytQoEOEJ84WjS\n/VZVxOJbnBJ3w4FBbq8CIYxNUFV/AV7BpMtZcsEE7b++M2+UiiY9fS4ddgyVEfUw2SHforo+t+Mf\nJiG/bbUZUvTA3AvMduaQLwhbE3RIBm4WkVpuCwl2BuuYQ435rEptvl/xLPesW0vD4ZhP+1xziPxx\n0aRv9cZYFncRkUTgMWCQ21pOEdYmqKr7gVHARKfPqSUXeNSjr+stbR/miSXX836Rnky7whvjHqBQ\nZCSZW7wxlsV1koE3VfV7t4WcIqxN0GE6pn/wtW4LCQlE8vVnYvWhjHjiLTqOu1+ezlWfCI8kR+wj\nMU8+jm3ylkSLO4hIbUzBkhzfGfUFYW+CasrsPwRMdfobWHLH7cDW3jq1XxfmdJlLh77dZPY5ewqf\ni2PkLRfHUR2uSfaydBDj7LQmA8NVNaCuO4W9CQKo6ofAF8CjbmsJakzV6MdwzgKnau83ujCn9WJa\n336bvPJyTobMJLJaYfaHXXmzEOQ2oBCmWntAYU3w/+kL9BCRGm4LCWK6Ar9xWtRvsvZZegcvtVpO\n8y53y3PZXhFmEFWpIAfTvKrS4lecYMh4oLsGYL1Oa4IOarpbDQdmODmNluxgMgCGcIaI8Gh9dPkt\nvN52EW1uv0+emZ6dYU3e8GH/N4OyeJPxwBvqhetSvsD+sv+TaZgCj/e7LSQI6QzsxvRn+Q9P6MNL\nOvNG+3m0v6+bzM7yBeo0YkrHceSgt0Ra/IuIXA20wHxABiTWBE9DTYOXu4Bkp/agJSuYlbM5CzxH\nBsBk7bPgTl64ex7tH+srk7plZeg0YorHkZr1VqmWgEFE8gPPYbbBAbuatyb4L1R1EzAOeNZui7NM\nOyCNLNSEG6cDX7yZtz2v0fWZYZLc8nyvTyW2SAypf3lDpMXvPAEsVtXFbgs5F/aX/Mw8AUQDvd0W\nEvCYqw+PAmPOtQo8nRnaY3grlrzyMrd/0FcmVTrXa48QXyiadFtWP8gQkbbA1UB/t7WcD2uCZ8DZ\nFt8GPPafNp2Wf9McSADezc6bXtbb/3cxX3+7mss/HybJUWd73WES4m2rzeBCRIpjWmfeFcjb4FNY\nEzwLqroVU2rrNTF18SxnZjAwDtWT2X3jFaxskoeTuoFan53tNQcoFBNFhs0bDhKcI6QXgedUdYXL\ncrKENcFz8zLwHWZ7bPk3Ig2AKsCrOXn7Qzol8zLW1P+cBhf1kBlTzvSaAxTOF0XGxtzItPiVh4H8\nBFF1JmuC58CpddYDaCkiXd3WE4AMBiZgUg9zxGTt80sn3rz9bW7u1Vumtj79e0kyPCGVGPJyfEeu\nlVp8jpgPxUeAroF4KfpsWBM8D05rzpuBKSJyodt6Agbzd9EYeDa3Q03Ufq+3Z96cpbScO0YG5z/1\n9WPkq5bIvpMe9bheeNNybpzixG9jrsNsd1lOtrAmmAXUlId/BHhHRBLc1hMgDAKmoJrqjcFK8/tt\nRdh7YD31/r5snUFU1YIctK02AxwxOeOvA6+o6jy39WQXa4JZRFWfB1YDL4T9/UGR8pjSY9lKgTsX\nHvVoPdZfsYbL6jwg0/qDSZkrQIpXTNbiU0ZhusYNdVtITgjvX+bs8yBQggCrh+YCA4CZqHo1nW2S\n9t3aiTeT3uXGsYNkbOl0osvaVpuBjYjciqkR2MW5WhZ0SAD0OQkqnDtQ64CBqvqG23r8jvn5NwI1\nUN3liynayKJNeTl+IpZ9fxwlqvQHdL1SVff4Yi5LzhGRxsA8oLmqbnBbT06xK8FsouYX/wbgKRGp\n77YeF+gDvOYrAwSoyuZWX1Cv+kn2Xn2Ub6sAO0TkFl/NZ8k+Yo5E3gbuDGYDBLsSzDEicj3wDNDU\nuVgd+ogUBLYCl+LDCKCIFG1Ix537GJm3BGNZzfNgcpPL2RWh+4j5d7AamKWqZ7zfGUzYlWAOUdX3\nMU1jFotIUbf1+IkewEJfGqBD+XW8fbQcW7mIv+upZgLlfTyv5Tw42VPzgWXAVJfleAW7EswlIjIS\naIk5FwndPhimYfYvQGtUv/PtVFIU2FEJYtYCdYE/7ErQdZyrMG9h/l/cpjlIlQxE7Eow9wwFfgDm\niUi022J8yK3Ad742QADH6O7eCmnPQcZIc/3ibmuA7uE0SpoOxGPOAUPCAMGuBL2CmNLyczDlt27S\nXKSRBSTmXuQGoBeqH/lvWinaGGquhDcj4GoCqFdtOOEY4CTgMqBFMFSGyQ7WBL2EiEQCc4EjwK3B\nemfqjIhcC4wELslqzUAvz98HaImq7Q3tZxwDHAO0Aq5W1QMuS/I6djvsJVQ1E5NjXAR40Tk/CRUG\nAONdMUDD00ANRK50af5wZhjQFmgVigYI1gS9iqqmY0rNFwNeEZF8LkvKPeYuZAXMgbg7qGYAIwjS\ntKxgRAwjMdkgLVU1ZPu8WBP0MmoKCrTDVFue42yTg5n+wKTclMvyEq8AlRG5zGUdIY+zBZ6IyQ+/\nUn14MT4QsGeCPkLMlZI3gbxAR/VStRW/YjrufQ5UIBAOw0XuB65F9Tq3pYQqTnGQacAlQOtQ3QKf\njl0J+gg1W7ibgb3AUhEp7LKknNAXmBUQBmh4HrgY2/fFJzgf3HOACzFb4JA3QLAm6FOcqzJ3AWuA\nlSJS2mVJWUckEXM3MHCyAsyZ63Sgl9tSQg0RKQAsBiKAa1Q1bKr3WBP0Map6UlX7Y5rPfBpE3evu\nB95F9U+3hfyLWcBNBOfKOiARkZLASsyl/85OgC9ssCboJ1R1PKY69VIx9+4CF5P58iDmcDywUN0N\nvA/c47aUUEBELgHWYrbBvULqfmsWsSboR5z6gzcAs0SktxOFC0TuANaj+qPbQs7Ck0BPwr3Cdy4R\nkZuBD4E+qjpWwzRKaqPDLuDUYluA+QR+MKC2H8ZYfgR6EKh9Y82Hx7dA74DVGMA4EeDHgHuB9qr6\nlcuSXMV+krqA042rEaY/60oRKeOuon9wDZAKfOK2kLNiPrlfxKxYLdnAqQU4D2gDNAx3AwRrgq6h\nqkeAzphMjM9FpLnLkk7RB5jsYopcVnkNuBGROLeFBAtOUG49sB1opoEX9HIFa4IuoobxwG2YNLtk\nV3OORWoCFwGB3zvF/AKvwZyxWs6BkwJ3N6YQapKq9nZy3S1YEwwI1JSnugRoAix38T5hb+BpJ1c3\nGHgTuMltEYGMs/19HXPxvZmqvuaypIDDmmCAoKp/Yc7jFgNfOpE7/yFSBOgIzPDrvLnjfaAFIrFu\nCwlERKQJ8A2wG2igqj+4LCkgsSYYQKjqCVUdDbQHRovIqyJSyE/Td8dcjt7tp/lyj+o+4EtMrTuL\ng4hEi8hY4B2gt6r2UtW0870vXLEmGICo6hrgYmA/8L2ItPbphKbSzQNAMHYOexfzoWEBnDawXwJV\ngDqqOt9lSQGPvScY4IjI1cCzmLSmfj7psyHSFbgH1au9PravEakErAJKBUFE22eIORJIwuSqPwS8\nEa6Xn7OLXQkGOE7QpBawB9ggInd4NdPEjNUXmOy1Mf3LNiADqOG2ELcQkTaYHjDlMKu/160BZh27\nEgwiRORSTAGBg5izng1eGLQx8BJQlWDtICYyC/ge1cCpeOMHRKQUpgHSpUBPVf3QZUlBiV0JBhGq\n+iXQANPQabmITPFC4KQvMCVoDdCwDGjhtgh/ISIxIjIE+A7YAtSyBphzrAkGGap6XFWfwhS+jAQ2\nikjPHPUzESkHNMcUKw1mVgGNCdyCFF7BufTcCdiI6UlfT1Ufs5Hf3GG3w0GOiNQFxgEVMUnxb2X5\nPEhkPCCYeofBjchvwFWo/uy2FF/gpFWOBqIwVV8CN7c7yLAmGCKISAtgLKDAEGDJOc1QJB6TQ1oP\nU9AhuBF5G3PP8VW3pXgT58rLaEzHv6GYqG8wH10EHHY7HCKo6jLMeeEEzGH5GhFpc45I8m3AqpAw\nQMM6oKHbIryFiDQUkfcx9yDfBmqo6hxrgN7HmmAI4ZTyfwNzpeYJYDywTkRuFJGIv19ojPFBTHHS\nUOFLTP510OKc+V0pIksxedGLgMqq+oy63/I0ZLHb4RDGKZ55I6asf2GMMb6oppbhVKB2yFwwFimO\nCRgkBtvP5AS1bsZE6gthtr+v2kov/sGaYBjgbImbYBqpN14PB+LgxeomTzk0MD/jHqAmQdIsXExH\nv26YVflWzDHGgnDs8+EmdjscBjh1C1eravtIuHw7HKoHD4vIByJy3T+2ysGK+TT/EXN1KGA5bcv7\nKsb4agI3qGozVX3PGqD/sSvBMMXJNe2Maa1ZEngFeFFVN7oqLDeIzAS+RvVpt6X8G6dGZFfgbuAk\n8Azwsqrud1WYxa4EwxVVTVXV51W1AdAa82/hIxH5XER6Ob1og43tQFm3RZxCRAqIyF0ishzTGKoy\nZvtbU1WnWAMMDOxK0PI3Tmn/FkAX4HrM9vIt4B1V/d1NbVlC5DagDaq3uidBimJK/ncALgc+Al4G\nFgZUV0HL31gTtJwREYkCrsZUm24H/AosBD4A1gXk2ZXIlcBIVJv6b0oRTF+WNs5TB9PLdy7G+A77\nS4slZ1gTtJwXZ4XYCLgWaAuUAj4+7dkUEKWbRCoAK1At59tppDxwBXAVpiXCUcydvg+B5TaXN7iw\nJmjJNs4h/1WY4gvNMYUcVmGyNtYCX7liBCYVcBeqXmvDKabqdm2gPtAYuBKTv7sSWIFJTwzJfOVw\nwZqgJVc428EKGINo5DwXApswwYDvTv2pqnt9LQZTYDUhJx3znM5stU57LsVsdbcBX2AM/hNgc0Cs\nfB5+f9MAAAGqSURBVC1ewZqgxeuISAzmbOyi0/68CDgObMbUwNuCMZffnecPrwQORHYBdTDd+/6t\nKxIogbkSVA6o5DwVMT05CgA/YKo0bwC+wqxqj+RalyVgsSZo8QvOirEYxmyqAFWB8kBp5ykJHAb2\nAvucZy9wCEgDUp0nHVMp56Tzp2K241FA9F/Qtwss+Nh8vSAmDa0wxvwKYtpP7sQEerZijHgr8DPw\nqy1QEH7kdVuAJTxwto+7nGf1v7/v5DkXcZ5E5ykCxAOxzlMIiAbEefI4f2ZizDHjTfg0xRjcDuDA\nac+fwN6AjGpbXMWuBC0WS1hjM0YsFktYY03QYrGENdYELRZLWGNN0GKxhDXWBC0WS1hjTdBisYQ1\n1gQtFktYY03QYrGENdYELRZLWGNN0GKxhDXWBC0WS1hjTdBisYQ11gQtFktYY03QYrGENdYELRZL\nWGNN0GKxhDXWBC0WS1hjTdBisYQ11gQtFktYY03QYrGENdYELRZLWGNN0PJ/7dSBAAAAAIAgf+tB\nLohgTYLAmgSBNQkCaxIE1iQIrEkQWJMgsCZBYE2CwJoEgbUAedG6Wh00ohAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7b830f0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "initialize_plot()\n",
    "axes = plt.gca()\n",
    "axes.set_xlim([-1.1,1.1])\n",
    "axes.set_ylim([-1.1,1.1])\n",
    "\n",
    "for i in range(n):\n",
    "    plot_hyp_line(pts[i-1],pts[i])\n",
    "plt.scatter(np.real(pts), np.imag(pts), c='k')\n",
    "plt.show() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<sup>1</sup> When I actually tried to recreate the graphic using arced lines, it just looked like a mess, so I wound up doing all this work for nothing, but c'est la vie.\n",
    "\n",
    "<sup>2</sup> <a href=\"https://mathworld.wolfram.com/OrthogonalCircles.html\">Two circles intersect at right angles if the radii of the two circles and the distance between their centers satisfy the Pythagorean theorem</a>:\n",
    "$$\n",
    "r^2_1 + r^2_2 = d^2\n",
    "$$\n",
    "Let the first circle be the unit circle centered at the origin, so $r_1 = 1$. Let $q$ be the center of the second circle. Since the first circle is centered at the origin, the distance between their centers is just $d = |q|$. Let $r_2 = |q - z|$ where $z$ is any point on the second circle. Therefore,\n",
    "$$\n",
    "1 + |q-z|^2 = |q|^2\n",
    "$$\n",
    "and expanding and simplifying,\n",
    "$$\n",
    "1 + |z|^2 - q\\bar{z} - \\bar{q}z = 0\n",
    "$$\n",
    "Solving for $\\bar{q}$ gives\n",
    "$$\n",
    "\\frac{1 + |z|^2 - q\\bar{z}}{z} = \\bar{q}\n",
    "$$\n",
    "$q$ and $\\bar{q}$ stay the same regardless of $z$. Choose two points $a$ and $b$ on the circle, then:\n",
    "$$\n",
    "\\frac{1 + |a|^2 - \\bar{a}q}{a} = \\frac{1 + |b|^2 - \\bar{b}q}{b}\n",
    "$$\n",
    "Then solve for $q$.\n",
    "$$\n",
    "q = \\frac{a\\left(1+\\left|b\\right|^2\\right)-b\\left(1+\\left|a\\right|^2\\right)}{a\\overline b-\\overline a b}\n",
    "$$\n",
    "\n",
    "$a$ and $b$ can be outside the unit circle and $q$ still gives the correct center of the orthogonal circle, although it can't really be interpreted in terms of the Poincare disk in that case."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}