First pass at documenting the cartopy project_segment algorithm

index.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How does cartopy project line segments?\n",
    "\n",
    "Cartopy goes through a number of steps to transform line segments from one coordinate system to a projection.\n",
    "\n",
    "**Note**: *This notebook is intended to illustrate the projection process for those wanting to maintain the cartopy codebase. It is not inteded to be definitive, nor to document public Cartopy APIs. If you are looking for usage documentation, please visit the cartopy docs.*\n",
    "\n",
    "Take the following example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAAC1CAYAAAD86CzsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAB8lJREFUeJzt3VuIVeUfx+HvaEfDDhpllnmIlMoITCvsRFhkJN0YZUbdSEEURVDRRXSAuugAUUkFRRFkEBRFJfWnpAwMCicKI1JJ6SREU5KVWon7f7EYxpyDZs5vnJnngcXgvGv2fkX2xzXvXmvttlarFQBqjBjoCQAMJ6ILUEh0AQqJLkAh0QUotF9fg3Pnzm11dHRUzQVgSGhvb/9fq9Wa29NYn9Ht6OjIypUr+2dWAENUW1vbkb2NWV4AKCS6AIVEF6CQ6AIUEl2AQqILUEh0AQqJLkAh0QUoJLoAhUQXoJDoAhQSXYBCogtQSHQBCokuQCHRBSjU5ydH/PhjsnhxctBBvW+jRiUTJyaHHFI1ZYB9w5YtSUdH960vba1Wq/fBtpmtZPc+rmf8+OTEE5OpU5uvndsJJzRxBtiX/f13zwHta9u8ubdHa2tvtVozexzZW9Ht/TGSCRO6IrxjlCdPTg444D89PMBua7WS1auTZcuS5cuTb7/tCuivv+7NZxrA6PZlxIhk0qR/HhlPnZrMnp0cemi/PS0wjHz/fRPZzm3Dhopn3cPojhs3szV//sps3Zpet40bk2++SbZv33vTPfDAZN685Oqrk0susTwB7L6NG5P33++K7OrV/fdc+++fHHlk9+2pp/YwujNnzmztzkew//VXsn59snZt17ZmTfP1u++aQ/o9ddhhyeWXJwsXJuefn4wcueePBQw9mzcnK1Z0Rba9fc+aM2JEMnZszxHtbRs9ullC3VlbWz9Hty9btyZff/3PIHduP/zw7x5r/PhkwYImwDNm9PyXBYa+zz5Lli5N3nsv+eij5sBvd40alZx7bjJnTnLmmcnRRzcBPfzwvXdQN6DR7csffzRB7jwqXru2Wdxet27XPzttWhPfq65q1oKBoa+9PbnrruSdd3b/Z0aObOI6Z06znXVWs4TZn/bZ6Pak1Uo+/jhZsiR5+eXkp592/TOzZjXrv1demYwb1/9zBGp9+WVy993Jq6/u3v6nntoV2fPOq39jflBFd0fbtjVrNEuWJK+9lvz+e9/7jxjRhPfRR5tfGYDBbf365N57kxdf7PvN+kmTkgsvbCJ7wQUD//oftNHd0ebNyZtvJi+9lLz9dnMic2/GjEkee6w5+rXuC4PPhg3J/fcnzz7b+2v9nHOSa69tQjtlSu38dmVIRHdHv/ySvPJKE+Dly3vf79JLk6efTo47rm5uwJ7r6EgefLC5/cDWrT3vM2NG8sADycUX77sHVX1Fd1De8GbMmOT665MPPmiuKHnooeTkk7vvt3RpcsopyTPP/LfT1oD+tWlTs4wwZUryyCM9B/ekk5qDrZUrk7lz993g7sqgjO6OJkxIbr89+fzz5n/Ind+V3LSpCfRFFzXrQ8C+Y/v25PHHm9jed1/y22/d95k8OXnhhWTVqmT+/MEb206DPrqd9tsvueOOJr5nn919fNmyZPr05Ikn9u7Vc8Ce2bgxueyy5JZbkp9/7j5+zDHJk08mX33VrN0OlQujhkx0O02blnz4YRPXnW83uXlzcvPNzSkk/XlpINC3Tz9NTj+9WQLc2dixycMPN+fw33DD0Lsp1pCLbtKcOnbTTc2vI3PmdB9fsSI57bRmOWLbtvr5wXDVajVnJMye3X25b/ToZl133brkttuSgw8ekCn2uyEZ3U6TJyfvvtv8I+98cvSffyZ33tlcnbJq1cDMD4aTLVuSRYuS665rXn87OuOM5IsvknvuGfp3GBzS0U2aRfdFi5orWubN6z7e3v7vLikE9szrryfPP9/9+zfe2CwJHn98/ZwGwpCPbqdjj03eeKO5um3s2K7vz5qV3HrrwM0LhosFC5Jrrun686hRzetx8eL+vxfCvmTYRDdpjnoXLmyOeq+4orkX5nPPNWc+AP2rra25WGn69OYN708+aV6Pw82wzM1RRzU301mzpvmkCqDGqFHJW28lRxwx9NduezMso9tJcKHexIkDPYOBNayWFwAGmugCFBJdgEKiC1BIdAEKiS5AIdEFKCS6AIVEF6CQ6AIUEl2AQqILUEh0AQqJLkAh0QUoJLoAhUQXoJDoAhQSXYBCogtQSHQBCokuQCHRBSgkugCFRBegkOgCFBJdgEKiC1BIdAEKiS5AIdEFKCS6AIVEF6CQ6AIUEl2AQqILUEh0AQqJLkAh0QUoJLoAhUQXoJDoAhQSXYBCogtQSHQBCokuQCHRBSgkugCFRBegkOgCFBJdgEKiC1BIdAEKiS5AIdEFKCS6AIVEF6CQ6AIUEl2AQqILUEh0AQqJLkAh0QUoJLoAhUQXoJDoAhQSXYBCogtQSHQBCokuQCHRBSgkugCFRBegkOgCFBJdgEKiC1BIdAEKiS5AIdEFKCS6AIVEF6CQ6AIUEl2AQqILUEh0AQqJLkAh0QUoJLoAhUQXoJDoAhQSXYBCogtQSHQBCokuQCHRBSgkugCFRBegkOgCFBJdgEKiC1BIdAEKiS5AIdEFKCS6AIVEF6CQ6AIUEl2AQqILUEh0AQqJLkAh0QUoJLoAhUQXoJDoAhQSXYBCogtQSHQBCokuQCHRBSgkugCFRBegkOgCFBJdgEKiC1BIdAEKiS5AobZWq9X7YFvbysK5AAwVHa1Wa25PA31GF4C9y/ICQCHRBSgkugCFRBegkOgCFPo/BzYRGBuAEPAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import cartopy.crs as ccrs\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "plt.figure()\n",
    "ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))\n",
    "ax.set_global()\n",
    "\n",
    "# Plot geodetic vertices, and let cartopy figure out how to\n",
    "# interpolate it & cut it just right.\n",
    "ax.plot(\n",
    "    [60, -60], [50, 50],\n",
    "    transform=ccrs.Geodetic(), color='blue', lw=4);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How did cartopy know that the line should be interpolated into two lines, each made up of multiple straight-line segments (giving the appearance of a curved line)? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Underneath the hood, cartopy has a (private) Interpolator object which knows how interpolate between two points (either as a straight line, or as a great circle, depending on the coordinate systems in use).\n",
    "\n",
    "We can hook into this interpolator (for debug purposes only: See cartopy.geodesic for a similar public API)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cartopy.trace import _Testing as TraceTesting\n",
    "import cartopy.crs as ccrs\n",
    "\n",
    "# Keep a hold of references to the coord systems.\n",
    "# If they get GCed, then you *will* segfault.\n",
    "source, dest = ccrs.Geodetic(), ccrs.PlateCarree(central_longitude=180)\n",
    "\n",
    "interp = TraceTesting.interpolator(source, dest)\n",
    "\n",
    "# Define a helper to pretty print the points that come out.\n",
    "def print_pt(pt):\n",
    "    print(f'{pt[\"x\"]:.2f}, {pt[\"y\"]:.2f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This interpolator knows how to calculate a point given a line start and end, and the fraction along the line desired (referred to a \"t\"):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t=0.0: point= -60.00, 50.00\n",
      "t=0.25: point= -38.18, 61.94\n",
      "t=0.5: point= -0.00, 67.27\n",
      "t=0.75: point= 38.18, 61.94\n",
      "t=1.0: point= 60.00, 50.00\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "for i in np.linspace(0, 1, 5):\n",
    "    print(f't={i}: point=', end=' ')\n",
    "    print_pt(TraceTesting.interp_t_pt(\n",
    "        interp, {'x': 120, 'y': 50}, {'x': -120, 'y': 50},\n",
    "        i))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The interpolator also knows how to project from the source to the destination coordinate system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-60.00, 50.00\n"
     ]
    }
   ],
   "source": [
    "print_pt(TraceTesting.interp_prj_pt(interp, {'x': 120, 'y': 50}))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The task of identifying whether a segment must be interpolated and/or cut is iterative.\n",
    "\n",
    "We start by setting up the full segment, and iterating a ``t_min`` and ``t_max`` until we find a values that produce straight (enough) lines *in both* projected and non-projected space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "i: 0 to 0.9; straight enough?: no\n",
      "i: 0 to 0.8; straight enough?: no\n",
      "i: 0 to 0.7; straight enough?: no\n",
      "i: 0 to 0.6; straight enough?: no\n",
      "i: 0 to 0.5; straight enough?: no\n",
      "i: 0 to 0.4; straight enough?: yes\n"
     ]
    }
   ],
   "source": [
    "from cartopy.trace import _Testing as TraceTesting\n",
    "\n",
    "line_start = {'x': 60, 'y': 50}\n",
    "line_end = {'x': -60, 'y': 50}\n",
    "\n",
    "# We choose a threshold that isn't as high as the one defined\n",
    "# in cartopy so that we can see the process in action.\n",
    "threshold = dest.threshold * 10\n",
    "\n",
    "i, good = 1.0, False\n",
    "while i > 0 and not good:\n",
    "    i -= 0.1\n",
    "    \n",
    "    good = TraceTesting.straight_and_within(\n",
    "        line_start, line_end, 0.0, i,\n",
    "        interp, threshold, dest.domain\n",
    "    )\n",
    "    print(f'i: 0 to {i:.1f}; straight enough?: {\"yes\" if good else \"no\"}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To visualise what is going on here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgQAAAEJCAYAAADrdD/VAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3XlclPX+Pv5rxAVRVLJciiIscRmZGcQF3GBM0UL9YCGKlhLuaZj9VNRCLS1LPWrmOfTT3FqRMI+e1DITcsMjLrhAKqckJTcwwAXZX98/bhlZZmBGkM3r+XjMg2Xuue83dDxzcd/v+32pRARERET0aKtT1QMgIiKiqsdAQERERAwERERExEBAREREYCAgIiIiMBAQERERgLqlPTlw4EBJSUmprLEQERHRQ3Ts2LGfRGSgsedKDQQpKSk4evTowxkVERERVSqVSvW4qed4yYCIiIgYCIiIiIiBgIiIiFDGHAIiouogJycHSUlJyMzMrOqhENUI1tbWsLe3R7169cx+DQMBEVV7SUlJsLW1xbPPPguVSlXVwyGq1kQEN27cQFJSEhwdHc1+HS8ZEFG1l5mZiebNmzMMEJlBpVKhefPmFp9RYyAgohqBYYDIfA/y74WBgIhqpbg4oFMn5WNtl5iYiG+++cai13z44YcPdKyVK1ciIyPjgV5L1RsDARHVOnfuAC+9BMTHA97eyte1WWUFgry8PAaCWoyBgIhqncBA4Pp1QAS4dg0YO7Z8+7tz5w68vb2h1WrRqVMnbN68GQBw7NgxeHh4wNXVFQMGDMCVK1cAADExMdBoNHB3d8fMmTPRqVMnAMDGjRvh4+ODwYMHw9HREatXr8by5cvh4uICNzc3/P333wCA33//HQMHDoSrqyt69+6Ns2fPAgACAgIQFBSEHj16oE2bNoiIiAAAzJ49G/v374dOp8OKFSuKjP3KlSvo06cPdDodOnXqhP3792P27Nm4e/cudDodRo0aBQDw8fGBq6sr1Go11qxZY3h948aNMW/ePHTv3h0ffPABLl++DL1eD71eX75fKlU/ImLy4erqKkREVS0+Pt7w+bRpIh4eph9OTiJ16ogocUB51KmjfN/Ua6ZNK/34ERERMm7cOMPXaWlpkp2dLe7u7nL9+nUREQkLC5PXX39dRETUarUcPHhQRESCg4NFrVaLiMiGDRvkueeek5s3b8r169elSZMmEhoaKiIib731lqxYsUJERPr27Svnz58XEZHDhw+LXq8XEZExY8aIr6+v5OXlSVxcnDz33HMiIhIZGSne3t5Gx75s2TJZtGiRiIjk5ubKzZs3RUSkUaNGRba7ceOGiIhkZGSIWq2WlJQUEREBIJs3bzZs5+DgIMnJyaX/wqhaKPzvpgCAo2LiPZ+3HRJRrXLhApCfX/R7+fnK91u3frB9Ojs7Y8aMGQgODsagQYPQu3dvnDlzBmfOnEH//v0BKKfTW7dujbS0NNy6dQs9evQAAIwcORI//PCDYV96vR62trawtbVF06ZNMXjwYMMxTp06hdu3b+PQoUMYNmyY4TVZWVmGz318fFCnTh107NgR165dK3PsXbt2RWBgIHJycuDj4wOdTmd0u1WrVmHr1q0AgEuXLiEhIQHNmzeHlZUVXnnlFQt/Y1QTMRAQUY2ycmXpz69fDwQFFZ03YGMDrF4NvP76gx3TyckJx44dw86dOzFnzhx4eXlh6NChUKvViI6OLrJtampqqftq0KCB4fM6deoYvq5Tpw5yc3ORn5+PZs2aITY2tszXK3/wla5Pnz7Yt28fduzYgddeew0zZ87E6NGji2wTFRWFPXv2IDo6GjY2NvD09DTcsmZtbQ0rK6syj0M1H+cQEFGtEhioTCS0tla+trYGBg9+8DAAAJcvX4aNjQ1effVVzJgxA8ePH0e7du2QnJxsCAQ5OTmIi4uDnZ0dbG1tcfjwYQBAWFiYRcdq0qQJHB0d8d133wFQ3vRPnjxZ6mtsbW1x69Yto8/9+eefaNGiBcaPH4+xY8fi+PHjAIB69eohJycHAJCeng47OzvY2Njg7NmzhrFbeiyq2RgIiKjWWb8eaNECUKmAli2BdevKt7/Tp0+jW7du0Ol0+OCDD/Duu++ifv36iIiIQHBwMLRaLXQ6HQ4dOgQAWLduHSZMmAB3d3eICJo2bWrR8b7++musW7cOWq0WarUa27ZtK3V7jUaDunXrQqvVlphUGBUVBZ1OBxcXF2zZsgXTpk0DAEyYMAEajQajRo3CwIEDkZubC41Gg5CQELi5uZk81oQJE/Diiy9yUmEtpCrtlFOXLl3k6NGjlTgcIqKSfvvtN3To0MGi18TFAcOHA5s3A2r1QxqYCbdv30bjxo0BAB999BGuXLmCTz75pHIHQY88Y/9uVCrVMRHpYmx7ziEgolpJrQbOnKmaY+/YsQOLFy9Gbm4uHBwcsHHjxqoZCJEFGAiIiCrY8OHDMXz48KoeBpFFOIeAiIiIGAiIiIiIgYCIiIjAQEBERERgICAiMktgYCBatGhhKCoyRkQQFBSE559/HhqNxrAIUHlt3LgRly9fNnv7qKgow5oIljh69CiCgoJK3SYxMdHk78DScVaWBQsWYNmyZVU9jGqPgYCIyAwBAQH48ccfS91m165dSEhIQEJCAtasWYPJkydXyLErMhDk5uaafF2XLl2watUqi8dXoLoGAjIPAwERkRn69OmDxx57rNRttm3bhtGjR0OlUsHNzQ1paWmGSuQCiYmJ6NChA8aPHw+1Wg0vLy/cvXsXABAbGws3NzdoNBoMHToUqampiIiIwNGjRzFq1CjodDrDtgVWrVqFjh07QqPRYMSIEUhMTMRnn32GFStWQKfTYf/+/QgICMDbb78NvV6P4OBgHDlyBD169ICLiwt69OiBc+fOAVCCxKBBgwAAycnJ6N+/Pzp37oyJEyfCwcEBKSkpAJQip+LjL2ucllY6i4ihOtrZ2dlQOV14jAAwdepUwzoPO3fuRPv27dGrVy8EBQUV2S4+Ph6enp5o06ZNuUJPrWaqBlFYf0zVzccfi+zda/y5vXuV56lWKlLjWlb/ccHD0VFEo1Eejo6lb1tW//E9Fy5cMFQZG+Pt7S379+83fN23b1+JiYkpsQ8rKys5ceKEiIgMGzZMvvzySxERcXZ2lqioKBERCQkJkWn3xuXh4VFiPwVat24tmZmZIiKSmpoqIiLz58+XpUuXGrYZM2aMeHt7S25uroiIpKenS05OjoiI/Pzzz/Lyyy+LSNEa5SlTpsiHH34oIiK7du0SAJKcnFzq+Esbp6WVzhEREdKvXz/Jzc2Vq1evytNPPy2XL18uUfU8ZcoU2bBhg9y9e1fs7e3ljz/+EBGRESNGGLabP3++uLu7S2ZmpiQnJ8tjjz0m2dnZRsdZm7D+mKpcaioQEQFcvQq0agX4+gJ2dhWw465dAT8/IDwcKLyOemTk/e8TFbC1BU6dUj7XaCrlkGJkKXiVSlXie46OjoYaYldXVyQmJiI9PR1paWnw8PAAAIwZM6ZIBbIpBX0EPj4+8PHxMbndsGHDDK2F6enpGDNmDBISEqBSqQwlR4UdOHDAUIc8cOBA2BX6R2xs/KV5kErnAwcOwN/fH1ZWVmjZsiU8PDwQExODJk2aGD3G2bNn0aZNGzg6OgIA/P39sWbNGsPz3t7eaNCgARo0aIAWLVrg2rVrsLe3L3XcjxoGAqowIsC8ecCyZYCVFZCRodTOBgUBM2YA77+vlM08ML1eedMvHAoKhwGWrTwayuo/LhAZCfTte/81lfC/D3t7e1y6dMnwdVJSEp588skS2xWuMLaysipxet0SO3bswL59+7B9+3YsXLgQcXFxRrdr1KiR4fOQkBDo9Xps3boViYmJ8PT0LLG9sXDzoON/kEpnU8evW7cu8vPzDV8X1DSXNl5jYy5tLsWjinMIqMLMmwcsXw5kZipd9CLKx8xM5fvz5lXAQQqHgnnzKjYMLFmivIkUFhmpfJ9qnpgYYO9e5RETUymHHDJkCL744guICA4fPoymTZuidevWZr22adOmsLOzw/79+wEAX375peFsganK4fz8fFy6dAl6vR5LlixBWloabt++XWZFcXp6Op566ikAMNmz0KtXL4TfO+u2e/dupKamlvkzmDrug1Q69+nTB5s3b0ZeXh6Sk5Oxb98+dOvWDQ4ODoiPj0dWVhbS09Pxyy+/AADat2+PP/74w3C2omDOAZmPZwioQqSmKmcG7oX1EjIylOcnTgRSUoDz54GEBOD335XX1KsH1K17/1H465LP6dHTdTK6L1yIY94hOJWoR90vTW1r/n4btesKWz8/qIydfaCaZ9as+59XQGD09/dHVFQUUlJSYG9vj/feew9jx47FZ599BgCYNGkSXnrpJezcuRPPP/88bGxssGHDBouOsWnTJkyaNAkZGRlo06aN4fUBAQGYNGkSGjZsiOjoaDRs2BCAMrnv1VdfRXp6OkQE06dPR7NmzTB48GD4+vpi27Zt+PTTT0scZ9asWRgzZgyWL1+OvgVnUYqZP38+/P39sXnzZnh4eKB169awtbXF7du3TY7f1DgBpdJ58uTJWLRoEXJycjBixAhotVqT+xo6dCiio6Oh1WqhUqmwZMkStGrVCgDg5+cHjUaDtm3bwsXFBQDQsGFD/Otf/8LAgQPx+OOPo1u3bmX8tqk41h9ThVi7Fpg+XTkjYIpKpbzxtmsHtG2rPJ5/XrmskJt7/5GTY/zzgq+fvRCJ0T/44Vf1ZPSJC8W/PMJx5gl9ma8r67nMTKBPXiS+yPLDfztPhkdcKC79IxzPjNGj0NlGqgIPUn9M5ZOVlQUrKyvUrVsX0dHRmDx5sslT/tVFQe20iGDKlClo27Ytpk+fXtXDqjKsP6ZKkZUF/PGH8ld+QgIQFlZ6GCjw7rvlvHRQ8Ff7znB46/VApB7BFXjZ4OpVPW78f5Mx4JuF+L5TCEI+0eOPaUp40WqVh06nfHziCTN2uGSJMhmy+CTImJiif8ESVTMXL16En58f8vPzUb9+faxdu7aqh1SmtWvXYtOmTcjOzoaLiwsmTpxY1UOqURgIyKScHCAx8f6bfsFp/oQE4MoV4Jln7v+l37at0j1v6pIBoJwJMPNyqnHGJhAam2hYDq1+iwR2hwIhIXg5NBQvh+uR6a5HXBxw8qTy+M9/lI82NkUDglYLODkpEyoNit8ZwcsQVEO0bdsWJ06cqOphWGT69OmP9BmB8mIgIOTlAQcPKndoFbzhJyQAly4BTz55/w3fyQnw9lY+Ojgo1+ALpKYC9+5QKvU4ZtxFZVpMjPE3/YJQEBNTvkBQPHDo9YCfH6zDw+Gq18PV9f6mIsDFi0owiI0FvvtOOftx5QqgVt8PCFqtHi7rw9HYzw+YPBkIDeUdEURULTEQPKJEgGPHgG++UU73t2oFdO+uvPH366d8dHSE2dfO7eyUWwuXL1cmEBZnYwO8/TbQrFk5Bl3aKfaCN/DyKB44SgkaKpUSihwcgCFD7n//1i0lWBWcTfjqK+DMGT0+rDcZQQsX4tc+IUhN10N7AXj22XLehklEVIEYCB4xCQlKCPjmG+Uv9pEjlbuy2rcv/77ff1/5WHwdgrw8JQwUPF9tGQscFgYNW1ugZ0/lUSBvTyTELxRxr4Sgy85QLPpIj6lJety6BXh5AZMmKbfLMxwQUVViIHgEXL2qnAX45hvlNPfw4cCXXyqXtyvyTUilAhYuVN78C69UOGxYOc8M1GSRkbDy9wO2hEN9bxLkYj8/LA4PR3InPSIilLszsrKUWzIDAoAylssnInoouDBRLZWeDmzcCPTvD3TooFznXrQISEoCPvkE6Nbt4f1FamcHjB8PhIQoHx/ZMACUehniiSeUaQUnTwLr1wMnTgDPPaeEgv/+V7msQ9XHjz/+iHbt2uH555/HRx99ZHSbP//8Ey+88AI0Gg08PT2RlJRUIcdm/XH5sP7YPAwEtUhmpjKxz9dXuQNg2zZgwgTg8mUlHHh5KesAUCWaNcv4JMhClydUKuUSw5dfKpd01Gpg1CjA1RVYswa4fRslV1HkCoqVKi8vD1OmTMGuXbsQHx+Pb7/9FvHx8SW2mzFjBkaPHo1Tp05h3rx5mDNnToUcn/XHVBkYCGq4vDxlDsC4ccodAatWAQMHKrcLbt2qnK4vtFgYVXOPPw7MnKnc4vnRR8CuXUq4++RQV+S+4qcEgYK7Ibp2rerhVl+FA1QFhKcjR47g+eefR5s2bVC/fn2MGDEC27ZtK7FdfHw8XnjhBQCAXq83ug3rj1l/XG2ZqkEU1h9XW/n5IkePirz9tsiTT4q4uIgsXSpy6VJVj4wehkuXRObNE/FtvldEuZIgWT+aqIGupSyuP9ZoROrVE3nmGeWjRlOu+uPvvvtOxo4da/j6iy++kClTppTYzt/fX1auXCkiIlu2bBEAkpKSUmQb1h+z/riysP64Fvv7b2D1amVyYE6OcofAnj3KHAGqveztgffeA3J7AfBSvjdiBNB2gjIRsU2bKh1e9WRnp6yCdfGicoqlnP3bYmat8bJlywx/sfbp0wdPPfUU6hq5Tsf6Y9YfV0cMBDXE998DU6cCL72kzAfo3p23qT1SIiNRd6Sfcn0IQPgrflh7IRzdu+vRpYty66K39yMyR8Sc+uOCyyohIcpiUPPnl2udCnNrjZ988kl8//33AJQ3wS1btqBp06YltmP9sen9CeuPqwznEFRz164p/782Z44yOf3zzwE3N4aBR07huxX0etTdEo7JXWJw8aJypujjj5UzBQsXAsnJVT3YaqDg9/X++/cXlyqHrl27IiEhARcuXEB2djbCwsIwpPCKVPekpKQY3qwWL16MwMBAs4/B+uOiWH9c+RgIqikR4OuvAY1GWTEwNhbo1auqR0VVpvjdCvfuVGjYEHjtNeDQIaVj4dIlwNkZODO6YifV1TiFf1/F7up4EHXr1sXq1asxYMAAdOjQAX5+flCr1QCAefPmYfv27QCUCW/t2rWDk5MTrl27hnfeecei42zatAkzZ86ERqNBbGws5t1rAiuoFS4+Wa+g/tjZ2RkuLi5F6o+3bt1qmFRY3KxZszBnzhz07NkTeXl5Rscyf/587N69G507d8auXbsM9celMTVOQKk/XrduHbRaLdRqtdEJl4UNHToUGo0GWq0Wffv2NdQfP/3004b641GjRhmtP+7Vqxdatmxp9OwMmcb642ror7+UU8CJicr96ZxMTpaIjgb+6RuJ0Bt+qD9tMhqsr/n9Caw/rnysP675LK0/5hmCakREuSSg0wFduihdAwwDZCl3d+Czc3r8qp6MBksW4s+XJtfoMEBV4+LFi+jatSu0Wi2CgoJqTP2xTqeDWq1Geno6648t9ChMQaoRLlxQFhFKTQV++UW5VED0oBrHRGLQxVCcHx6C5l+F4v/P0iNgk97ssioi1h8/eniGoIrl5wOffqqcCejfHzh8mGGAKsC9SXVOYe/DKiIc1mdi0K0bcPp0VQ+MiKorniGoJHFxSqnQ5s3K0rSAshrd2LHKpYKDB4F27ap2jFSLFJpE12yoHqN99MjfqLQqzp6tFCrV4Z8DRFQI/y+hEty5o6wfEB+v3Cuenq5M+u7RQ7mlcN8+hgF6uFQq4PXXgSNHgH//G3jhBWXNHiKiAgwElSAwELh+XTkTcPWqchvh7t3KWd033+RfalR5HB2BqCil76JLF+Crr9iqSEQKvhU9ZOvXAzt2KE2EgNJ7f+eOspjMvRU2iSqVlRUQHKyE0o8+Ar7WLkH6vyMfzfUKLGBO/fHFixeh1+vh4uICjUaDnTt3VsixV65ciYyMDLO3//e//220jbEs27dvN/mzFSheLlSYpeOsLAEBAYbSJDKNgeAhmzNHCQCFZWcr3yeqSjodcPQocLNdVzQd2leZYMD7XI0yt/540aJF8PPzw4kTJxAWFoY33nijQo5fkYGgtCV7hwwZgtmzZ1s8vgLVNRCQeRgIHrLFi4FCS4gDAGxslL/MiKqatTVQQe9Z1Uvh+mOg3Gc/zK0/VqlUuHnzJgBleWBjfQdRUVHw9PSEr68v2rdvj1GjRhnW4f/ll1/g4uICZ2dnBAYGIisrC6tWrcLly5eh1+uhN7KexOzZsw31xzNmzMChQ4ewfft2zJw5EzqdDr///js8PT0xd+5ceHh44JNPPsF//vMfdO/eHS4uLujXr5+hUGjjxo2YOnUqAKWu2M3NDV27dsW8efPQuHFjwzFv375dYvxljfPYsWPw8PCAq6srBgwYgCtXrgAAPD09ERwcjG7dusHJycmwsmJmZiZef/11wyqMkff+exYeIwAMGjQIUVFRAIB169bByckJnp6eGD9+fJHt9u3bV6JimYriXQYPWWAg8NNPwPbtymUDa2tg8GBlghdRtRATYyhNQkxM9V/E6K23lLW8S5OaCrz77v0q0N9+Uz43dQpfpyu1NOmvv/7C008/bfja3t4e//3vf0tst2DBAnh5eeHTTz/FnTt3sGfPHqP7O3HiBOLi4vDkk0+iZ8+eOHjwILp06YKAgAD88ssvcHJywujRoxEaGoq33noLy5cvR2RkJB5//PEi+/n777+xdetWnD17FiqVCmlpaWjWrBmGDBmCQYMGwdfX17BtWloafv3113u/nlQcPnwYKpUKn3/+OZYsWYJ//OMfRfY9bdo0TJs2Df7+/vjss8/KHH9QUJDJcebk5ODNN9/Etm3b8MQTT2Dz5s145513sH79egDKWYsjR45g586deO+997Bnzx7885//BACcPn0aZ8+ehZeXF86fP2/yv9Hly5excOFCHD9+HLa2tujbty+0Wq3h+StXruDAgQM4e/YshgwZUuR3QwqeIagE69cDLVooM71btgTWravqEREVUrDufwWs+V9t2NkpAeDUKeXRoUO5KpDNrT/+9ttvERAQgKSkJOzcuROvvfZakWa+At26dYO9vT3q1KkDnU6HxMREnDt3Do6OjnBycgKg1B/v27ev1HE1adIE1tbWGDduHL7//nvY2NiY3Hb48OGGz5OSkjBgwAA4Oztj6dKlRhsSo6OjDXXFI0eOLHP8pTl37hzOnDmD/v37Q6fTYdGiRUhKSjI8//LLLwMoWqV84MABvPbaawCU4iIHB4dSA8GRI0fg4eGBxx57DPXq1StRHW2sYpmK4hmCStCokfKHScE6BMUvIRCRBcypPwaUywR9+95/TSXUH69btw4//vgjAMDd3R2ZmZlISUlBixYtimxnrIq3rPpeY+rWrYsjR47gl19+QVhYGFavXo29BWd7iilcf/zmm2/i7bffxpAhQxAVFYUFCxZYdFxLq4RFBGq1GtHR0aXur/C+Hmb98YP8rh8FPENQSdRq4MyZ+4sSEdFDFBmpLPKxd6/y8PMrOqfAQubWHz/zzDOGOt7ffvsNmZmZeOKJJ8w6Rvv27ZGYmIj//e9/AMyrP759+zbS09Px0ksvYeXKlYbyIUvqjzdt2mR0Gzc3N2zZsgUAEBYWZtbPYOq47dq1Q3JysiEQ5OTkGD0rUVifPn3w9ddfAwDOnz+Pixcvol27dnj22WcRGxtrqH4+cuQIAOWsxa+//orU1FTk5uYaxk7mYyAgotrn3tLNhksh4eHK9x6QufXH//jHP7B27VpotVr4+/tj48aNRi8tGGNtbY0NGzZg2LBhcHZ2Rp06dTBp0iQAwIQJE/Diiy+WmKx369YtDBo0CBqNBh4eHlixYgUAYMSIEVi6dClcXFzw+++/lzjWggULMGzYMPTu3bvE9f4CK1euxPLly9GtWzdcuXLFrCphU+OsX78+IiIiEBwcDK1WC51Oh0OHDpW6rzfeeAN5eXlwdnbG8OHDsXHjRjRo0AA9e/aEo6MjnJ2dMWPGDHTu3BkA8NRTT2Hu3Lno3r07+vXrh44dO7L+2EKsPyaiao/1x5UvIyMDDRs2hEqlQlhYGL799lujd1ZUJwX1x7m5uRg6dCgCAwMxdOjQqh5WlbG0/phzCIiIqIRjx45h6tSpEBE0a9bMcEdAdbZgwQLs2bMHmZmZ8PLygo+PT1UPqUZhICAiohJ69+6NkydPVvUwLLJs2bKqHkKNxjkERERExEBAREREDAREREQEBgIiIiICAwERkVnMqT+ePn06dDoddDodnJyc0KxZswo59saNG3H58mWzt4+KiirzPn9jjh49iqCgoFK3SUxMRKdOnYw+Z+k4K8uCBQs44dAMDAREVLsUbzoEyt12aG798YoVKxAbG4vY2Fi8+eabhjX6y6siA0Fpywx36dIFq1atsnh8BaprICDzMBAQUe3StWvRpYoLljHu2vWBd2lu/XFh3377Lfz9/Ut8PzExER06dMD48eOhVqvh5eWFu3fvAgBiY2Ph5uYGjUaDoUOHIjU1FRERETh69ChGjRoFnU5n2LbAqlWrDPXHI0aMQGJiIj777DOsWLECOp0O+/fvR0BAAN5++23o9XoEBwfjyJEj6NGjB1xcXNCjRw+cO3cOgBIkBg0aBABITk5G//790blzZ0ycOBEODg5ISUkBoASk4uMva5y///47Bg4cCFdXV/Tu3Rtnz54FAAQEBCAoKKhENbGIYObMmejUqROcnZ2xefPmEmMEgKlTp2Ljxo0AgJ07d6J9+/bo1asXgoKCimwXHx8PT09PtGnTplyhp1YTEZMPV1dXISKqavHx8fe/mDZNxMOj9IdGI1KvnsgzzygfNZrSt582rdTjf/fddzJ27FjD11988YVMmTLF5PaJiYnSqlUryc3NLfHchQsXxMrKSk6cOCEiIsOGDZMvv/xSREScnZ0lKipKRERCQkJk2r1xeXh4SExMjNFjtW7dWjIzM0VEJDU1VURE5s+fL0uXLjVsM2bMGPH29jaMJz09XXJyckRE5Oeff5aXX35ZREQiIyPF29tbRESmTJkiH374oYiI7Nq1SwBIcnJyqeMvbZx9+/aV8+fPi4jI4cOHRa/XG8bm6+sreXl5EhcXJ88995yIiEREREi/fv0kNzdXrl69Kk8//bRcvny5yBgLxrlhwwa5e/eu2Nvbyx9//CEiIiNGjDBsN3/+fHF3d5fMzExJTk6Wxx57TLKzs42OszYp8u/mHgBHxcR7PhcmIqLax856srgYAAAStklEQVQOaN0auHgReOaZclUfA+bXHxcICwuDr68vrKysjD7v6OgInU4H4H7lb3p6OtLS0gyFRmPGjClR4WuMRqPBqFGj4OPjU+rKfMOGDTOMJz09HWPGjEFCQgJUKhVycnJKbH/gwAFs3boVADBw4EDYFfodGht/aW7fvo1Dhw4V+XmysrIMnxurJj5w4AD8/f1hZWWFli1bwsPDAzExMWjSpInRY5w9exZt2rSBo6MjAMDf3x9r1qwxPO/t7Y0GDRqgQYMGaNGiBa5duwZ7e/tSx/2oYSAgoprFnPrjgssEISFAaCgwf36l1B8XCAsLwz//+U+TzxevDy5+et0SO3bswL59+7B9+3YsXLjQZItg4frjkJAQ6PV6bN26FYmJifD09CyxvbEQ9KDjz8/PR7NmzQxtjKXtr+C4po5fEfXH5lQ2P4o4h4CIapeCMBAeDrz/vvKxkuqPAeDcuXNITU2Fu7u7Rcdo2rQp7OzssH//fgDm1R8XVADr9XosWbIEaWlpuH37tkX1xwXX34vr1asXwsPDAQC7d+9GampqmT+DqeM2adIEjo6O+O677wAob95lLYvcp08fbN68GXl5eUhOTsa+ffvQrVs3ODg4ID4+HllZWUhPTzfUTbdv3x5//PGH4WxFwZwDMh8DARHVLoWrj4FKrT8GlMmEI0aMMLv2uLBNmzZh5syZ0Gg0iI2Nxbx58wAoE+8mTZpUYrJeXl4eXn31VTg7O8PFxQXTp09Hs2bNMHjwYGzdutUwqbC4WbNmYc6cOejZsyfy8vKMjmX+/PnYvXs3OnfujF27dqF169awtbUtdfymxgkAX3/9NdatWwetVgu1Wl3mpMyhQ4dCo9FAq9Wib9++WLJkCVq1aoWnn34afn5+hkslLi4uAICGDRviX//6FwYOHIhevXqhZcuWrD+2EOuPiajaY/1x5cvKyoKVlRXq1q2L6OhoTJ482eQp/+qioP5YRDBlyhS0bdsW06dPr+phVRnWHxMRUbldvHgRfn5+yM/PR/369bF27dqqHlKZ1q5di02bNiE7OxsuLi6YOHFiVQ+pRuEZAiKq9niGgMhylp4h4BwCIiIiYiAgIiIiBgIiIiICAwERERGBgYCIyCzm1B8DQHh4ODp27Ai1Wo2RI0dWyLFZf1w+rD82DwMBEdUuVVh/nJCQgMWLF+PgwYOIi4vDSnOWWTYD64+pMjAQEFHtUoX1x2vXrsWUKVMMRUAtWrQosQ3rj1l/XG2ZqkEU1h8TUTVRU+qP/+///k9mzpwpPXr0kO7du8uuXbtKbMP6Y9YfVxbWHxMRVVH9cW5uLhISEhAVFYWkpCT07t0bZ86cQbNmzYpsx/pj1h9XRwwERFSzVOP6Y3t7e7i5uaFevXpwdHREu3btkJCQgK7FLlew/tj0/oT1x1WGcwiIqHapwvpjHx8fRN47TkpKCs6fP482bdqYdQzWHxfF+uPKx0BARLVLFdYfDxgwAM2bN0fHjh2h1+uxdOlSNG/e3OzjsP74PtYfVz6WGxFRtcdyo8rH+uOaj/XHRERUbqw/fvTwDAERVXs8Q0BkOdYfExERkcUYCIiIiIiBgIiIiBgIiIiICAwERFSbGGs6LFDOxkNz6o83btyIJ554AjqdDjqdDp9//vkDH6/4fll//OBYf2weBgIiqj2KNx0WKGfjobn1xwAwfPhwxMbGIjY2FuPGjXug4xXH+mOqDAwERFR7FKxKaKz+uPDqhRYyt/7YHKw/Zv1xtWWqBlFYf0xE1URNqT/esGGDtGrVSpydneWVV16RixcvltiG9cesP64sltYf8wwBEdU+heuPW7eutPrjwYMHIzExEadOnUK/fv0wZswYo/szt/543759ZY6tYE3/r776CnXrml58tnj98bBhw9CpUydMnz7daEPigQMHMGLECAAVW3+s0+kwceJEXLlyxfC8JfXHphirPy6soP748ccfN9QfU1FcupiIapZqXH9cuMho/PjxCA4ONro/1h+b3p+w/rjK8AwBEdUuVVh/XPiv3u3bt1u03DLrj4ti/XHlYyAgotrD2ARCYxMNLWRu/fGqVaugVquh1WqxatUqk2+2prD++D7WH1c+lhsRUbVndrnRkiXKrYXGLg9ERgIxMcCsWRU/wFqI9cc1H+uPiejRVdqbvV5frnkEjxrWHz96eIaAiKo91h8TWY71x0RERGQxBgIiIiJiICAiIiJOKiSiWig1FYiIAK5eBVq1Anx9y71YIVGtxzMERFRriCiLEz75JDB9urJA4fTpytchIcrztVFiYiK++eYbi17z4YcfPtCxxo0bZ7LpsUBAQIChpKiwBxlnZWncuLHZ2+7duxeDBw+Gs7Mz3N3dsXLlSpPrOcTGxmLnzp0WjcWcqu0CERERUKlUqIgbABgIiKjWmDcPWL4cyMwE7txRAsCdO8rXy5crz9dGFRkIRKTI0sDFff755+jYsaNFxypQnQOBuUJDQ7FkyRIsXrwYp0+fxp49e5CRkYERI0YYXT7Z0kBgSdX2rVu3sGrVKnTv3v2Bf57CGAiIqFZITQWWLQMyMow/n5GhPJ+WZvm+79y5A29vb2i1WnTq1MmwLO6xY8fg4eEBV1dXDBgwwLB0cUxMDDQaDdzd3Q0VvoCyTLCPjw8GDx4MR0dHrF69GsuXL4eLiwvc3Nzw999/A7C8Knj27NnYv38/dDodVqxYUWTsV65cQZ8+faDT6dCpUyfs378fs2fPxt27d6HT6TBq1ChDJfMbb7yBzp0749KlS5g8eTK6dOkCtVqN+fPnG/bn6elp+Gt03bp1cHJygqenJ8aPH4+pU6cattu3b59F4wSApUuXomvXrtBoNIZjWloXXXyMKSkpePbZZ+/9byDDsMrh8OHD0b179yJ/Wb/zzjvQarVwc3MzWn6UkJCA8PBw/PDDD4b/po0aNcLcuXPRvn37EmdFsrOzMW/ePGzevBk6nc6s5ZQtqdoOCQnBrFmzYG1tXeZ+zWKqBlFYf0xE1UThGldT7cdOTiJ16ogo5wWMP+rUUbazsP1YIiIiZNy4cYav09LSJDs7W9zd3eX69esiIhIWFiavv/66iIio1Wo5ePCgiIgEBweLWq0WEaUe+bnnnpObN2/K9evXpUmTJhIaGioiIm+99ZasWLFCRCyvCi5eCVzYsmXLZNGiRSIikpubKzdv3hQRkUaNGhm2uXDhgqhUKomOjjZ878aNG4bXeHh4yMmTJ0XkfsXxX3/9JQ4ODnLjxg3Jzs6WXr16GSqhH2ScP/30k4wfP17y8/MlLy9PvL295ddffy13XXRycrI4ODiIiMjSpUtlwoQJIiJy+vRpsbKyMmwHQLZv3y4iIjNnzpSFCxeWGOOcOXNk9+7dkpeXJ2+88YZ07txZ5s+fL0FBQfL333/LkCFDSrxmw4YNRaqy9+7dK1qttsTD3d1dRMyv2j5+/LihttpU7bSl9cecVEhEtUJ2NlDKmW4AyvPZ2Zbv29nZGTNmzEBwcDAGDRqE3r1748yZMzhz5gz69+8PQDnV27p1a6SlpeHWrVvo0aMHAGDkyJH44YcfDPvS6/WwtbWFra0tmjZtisGDBxuOcerUqSJVwQWysrIMnxurCi5N165dERgYiJycHPj4+Bhqi4tzcHCAm5ub4evw8HCsWbMGubm5uHLlCuLj46HRaAzPHzlyBB4eHnjssccAKPXK58+ff+Bx7t69G7t37zZ0E9y+fRsJCQl45plnzK6LLvw7M+bAgQOYNm0aAKBTp05Ffp769etj0KBBhmP8/PPPJV5/8uRJzJkzB//5z39Qr149HDt2DMuXL0diYiLs7OxKLZQqoNfrS10CWsyo2s7Pz8f06dMt7sooCwMBEdUoptqP165VJhDeuWP6tY0aAXPnAuPHW3ZMJycnHDt2DDt37sScOXPg5eWFoUOHQq1WIzo6usi2ZbUCFq7hrVOnjuHrOnXqIDc394GqgkvTp08f7Nu3Dzt27MBrr72GmTNnYvTo0SW2K1yPfOHCBSxbtgwxMTGws7NDQECAoWbY3GNbOk4RwZw5c0osN5yYmGhx3XLhiuTC4y5tHPXq1TO88ZqqRxYRWFlZ4ezZsxg4cCAA4MUXX8SpU6eQlZVVZJymREZGGu1XsLGxwaFDh8yq2r516xbOnDljqK2+evUqhgwZgu3bt6NLF6OLEJqFcwiIqFbw9QVMTPQ2yMsDyvgj0qjLly/DxsYGr776KmbMmIHjx4+jXbt2SE5ONgSCnJwcxMXFwc7ODra2tjh8+DAAICwszKJjPUhVcGl1x3/++SdatGiB8ePHY+zYsTh+/DgA5Q0wJyfH6Gtu3ryJRo0aoWnTprh27Rp27dpVYptu3brh119/RWpqKnJzc7Fly5Yyf7bSxjlgwACsX78et2/fBgD89ddfuH79usl9lVYX/eyzz+LYsWMAUOS6fuFK5/j4eJw+fbrMMRfm7OyM6OhotGvXDrt37wYA/PTTTxARfPzxx/D19S3zZy44Q1D8cejQIQDmVW03bdoUKSkpSExMRGJiItzc3ModBgAGAiKqJezsgBkzABsb48/b2CjPN2tm+b5Pnz6Nbt26QafT4YMPPsC7776L+vXrIyIiAsHBwdBqtdDpdIb/U1+3bh0mTJgAd3d3iIjFNbyWVgVrNBrUrVsXWq22xGS9qKgo6HQ6uLi4YMuWLYZT5hMmTDBUCBen1Wrh4uICtVqNwMBA9OzZs8Q2Tz31FObOnYvu3bujX79+6NixY5k/Z2nj9PLywsiRI+Hu7g5nZ2f4+vqWeQreVF30jBkzEBoaih49eiAlJcWw/RtvvIHk5GRoNBp8/PHH0Gg0Fv23GTNmDObOnQtvb2/cvXsXrq6uSEtLQ1xcHBo3bozAwMASr9Hr9YiPjzd7UqG5VdsPhanJBcJJhURUTRibHGVMfr7Iu++KWFuLNGokolIpH62tle/n5z/kgd5z69Ytw+eLFy+WoKCgyjlwJSv4OXNycmTQoEHy/fffV/GISpebmyt3794VEZH//e9/4uDgIFlZWRbtY+nSpfLyyy/Ln3/+KSIiGRkZ8tVXX8nFixcrfLzlxUmFRPTIUqmAhQuBt98uulLhsGEPdmbgQe3YsQOLFy9Gbm4uHBwcKnzyV3WxYMEC7NmzB5mZmfDy8oKPj09VD6lUGRkZ0Ov1yMnJgYggNDQU9evXt2gfM2bMwM6dOzF+/Hhcv34dTZs2hb+/P5566qmHNOrKw/pjIqr2WH9MZDnWHxMREZHFGAiIqEYo7WwmERX1IP9eGAiIqNqztrbGjRs3GAqIzCAiuHHjhsVLGnNSIRFVe/b29khKSkJycnJVD4WoRrC2toa9vb1Fr2EgIKJqr169enB0dKzqYRDVarxkQERERAwERERExEBAREREYCAgIiIiMBAQERERGAiIiIgIDAREREQEBgIiIiICAwERERGBgYCIiIjAQEBERERgICAiIiIwEBAREREYCIiIiAgMBERERAQGAiIiIgIDAREREYGBgIiIiMBAQERERGAgICIiIjAQEBERERgIiIiICAwEREREBAYCIiIiAgMBERERgYGAiIiIwEBAREREYCAgIiIiMBAQERERGAiIiIgIDAREREQEBgIiIiICAwERERGBgYCIiIjAQEBERERgICAiIiIwEBAREREYCIiIiAgMBERERAQGAiIiIgIDAREREYGBgIiIiMBAQERERGAgICIiIjAQEBERERgIiIiICAwEREREBAYCIiIiAgMBERERgYGAiIiIwEBAREREYCAgIiIiMBAQERERGAiIiIgIDAREREQEBgIiIiICAwERERGBgYCIiIjAQEBERERgICAiIiIwEBAREREYCIiIiAgMBERERAQGAiIiIgIDAREREYGBgIiIiMBAQERERGAgICIiIjAQEBERERgIiIiICAwEREREBAYCIiIiAgMBERERgYGAiIiIwEBAREREYCAgIiIiMBAQERERGAiIiIgIDAREREQEBgIiIiICAwERERGBgYCIiIjAQEBERERgICAiIiIwEBAREREYCIiIiAgMBERERAQGAiIiIgIDAREREYGBgIiIiMBAQERERGAgICIiIjAQEBERERgIiIiICAwEREREBAYCIiIiAgMBERERgYGAiIiIwEBAREREYCAgIiIiMBAQERERGAiIiIgIDAREREQEBgIiIiICAwERERGBgYCIiIjAQEBERERgICAiIiIwEBAREREYCIiIiAgMBERERAQGAiIiIgIDAREREYGBgIiIiMBAQERERGAgICIiIjAQEBERERgIiIiICAwEREREBAYCIiIiAgMBERERAVCJiOknVaqjlTgWIiIierhSRGSgsSdKDQRERET0aOAlAyIiImIgICIiIgYCIiIiAgMBERERgYGAiIiIAPw/AlqVGnGkv/gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 648x648 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import cartopy.crs as ccrs\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "def points_to_verts(points):\n",
    "    return np.array([[point['x'], point['y']] for point in points])\n",
    "\n",
    "\n",
    "def visualise_iteration(line_start, line_end, source, dest):\n",
    "    interp = TraceTesting.interpolator(source, dest)\n",
    "\n",
    "    plt.figure(figsize=(9, 9))\n",
    "    ax = plt.axes(projection=dest)\n",
    "    \n",
    "    verts = points_to_verts([line_start, line_end])\n",
    "\n",
    "    # Plot these geodetic vertices, and figure out how to\n",
    "    # interpolate it & cut it just right.\n",
    "    ax.plot(\n",
    "        verts[:, 0], verts[:, 1],\n",
    "        transform=source, color='blue', lw=1)\n",
    "\n",
    "    p0 = TraceTesting.interp_prj_pt(interp, line_start)\n",
    "    ax.plot(p0['x'], p0['y'],\n",
    "            marker='d',\n",
    "            color='blue',\n",
    "            transform=dest, label='segment start')\n",
    "\n",
    "    i, good = 1.0, False\n",
    "    while i > 0 and not good:\n",
    "        good = TraceTesting.straight_and_within(\n",
    "            line_start, line_end, 0.0, i,\n",
    "            interp, threshold, dest.domain\n",
    "        )\n",
    "        p1 = TraceTesting.interp_t_pt(interp, line_start, line_end, i)\n",
    "        ax.plot(p1['x'], p1['y'],\n",
    "            marker='o' if good else 'x', markersize=(1-i)*10 + 3,\n",
    "            color='blue' if good else 'red',\n",
    "            transform=dest, label=f'segment straight enough @ t={i:.1f}' if good else f'{i:.1f} not straight enough')\n",
    "        i -= 0.1\n",
    "\n",
    "    plt.legend()\n",
    "    ax.set_global()\n",
    "\n",
    "ax = visualise_iteration({'x': 60, 'y': 50}, {'x': 180, 'y': 50}, source, dest)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fact that the \"straightness\" metric is based on both source and projected coordinate systems is a key piece of being able to identify wrap-arounds and cuts.\n",
    "The following short line segment demonstrates this clearly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgQAAAEJCAYAAADrdD/VAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAE+1JREFUeJzt3X9sVfX9x/HXbWm9a6ld5SvK7ChVg8ptby8WpZcpbWXDSq0pGYypm0CxLLKlbAkEywYuYQ5/EGGLEUMGuhkXYKBRBAZDRdksEepAaIc0pl2oYn9pQWihve37+8ddLy39eYv0F89HcnN77z3nns8tbc6Tc8/tx2FmAgAAV7aQ/h4AAADofwQBAAAgCAAAAEEAAABEEAAAABEEAABA0rCuHszIyLDq6uq+GgsAALiMCgsLd5lZRkePdRkE1dXVOnjw4OUZFQAA6FMOh+P/OnuMtwwAAABBAAAACAIAAKBuziEAgIGgsbFR5eXlOnfuXH8PBRgUnE6nYmNjFRYW1uN1CAIAA155ebmioqI0ZswYORyO/h4OMKCZmWpqalReXq74+Pger8dbBgAGvHPnzmnEiBHEANADDodDI0aMCPqIGkEAYFAgBoCe683vC0EAYEgqKpISEvzXQ11ZWZn++te/BrXO73//+15ta82aNaqrq+vVuhjYCAIAQ87Zs9K0aVJxsZSZ6b89lPVVEDQ1NREEQxhBAGDIycmRKislM6miQpo379Ke7+zZs8rMzFRSUpISEhK0adMmSVJhYaFSU1OVnJyse++9VydPnpQkHThwQG63W16vV4sXL1ZCQoIk6eWXX1Z2draysrIUHx+v559/Xs8995zGjx+vlJQUffnll5KkTz/9VBkZGUpOTtbdd9+tY8eOSZLmzJmjvLw8TZo0STfeeKO2bNkiSXr88ce1b98+eTwerV69us3YT548qcmTJ8vj8SghIUH79u3T448/rvr6enk8Hj388MOSpOzsbCUnJ8vlcmndunWB9YcPH67ly5dr4sSJevLJJ/X5558rPT1d6enpl/ZNxcBjZp1ekpOTDQD6W3FxceDrhQvNUlM7v4wdaxYSYubPAf8lJMR/f2frLFzY9fa3bNlijz76aOB2bW2tNTQ0mNfrtcrKSjMz27hxo82dO9fMzFwul/3rX/8yM7MlS5aYy+UyM7OXXnrJbrrpJjt9+rRVVlba1VdfbWvXrjUzs1/+8pe2evVqMzO755577Pjx42Zmtn//fktPTzczs9mzZ9uMGTOsqanJioqK7KabbjIzs3fffdcyMzM7HPuqVavsd7/7nZmZ+Xw+O336tJmZRUZGtlmupqbGzMzq6urM5XJZdXW1mZlJsk2bNgWWi4uLs6qqqq6/YRgQWv/etJB00DrZ5/OxQwBDSmmp1Nzc9r7mZv/9o0b17jkTExO1aNEiLVmyRPfff7/uvvtuHT16VEePHtUPfvADSf7D6aNGjVJtba2+/vprTZo0SZL00EMP6a233go8V3p6uqKiohQVFaXo6GhlZWUFtvHxxx/rzJkz+uCDDzRz5szAOufPnw98nZ2drZCQEI0bN04VFRXdjv2OO+5QTk6OGhsblZ2dLY/H0+Fyf/zjH/X6669Lkk6cOKGSkhKNGDFCoaGh+uEPfxjkdwyDEUEAYFBZs6brxzdskPLy2p43EBEhPf+8NHdu77Y5duxYFRYWaseOHcrPz9fUqVM1ffp0uVwuFRQUtFn2q6++6vK5rrrqqsDXISEhgdshISHy+Xxqbm7Wt7/9bR06dKjb9f3/4eva5MmT9f7772v79u366U9/qsWLF+uRRx5ps8zevXu1Z88eFRQUKCIiQmlpaYGPrDmdToWGhna7HQx+nEMAYEjJyfGfSOh0+m87nVJWVu9jQJI+//xzRURE6Cc/+YkWLVqkjz76SLfccouqqqoCQdDY2KiioiLFxMQoKipK+/fvlyRt3LgxqG1dffXVio+P19/+9jdJ/p3+4cOHu1wnKipKX3/9dYeP/fe//9XIkSOVm5urefPm6aOPPpIkhYWFqbGxUZJ06tQpxcTEKCIiQseOHQuMPdhtYXAjCAAMORs2SCNHSg6HdN110vr1l/Z8R44c0Z133imPx6Mnn3xSv/nNbxQeHq4tW7ZoyZIlSkpKksfj0QcffCBJWr9+vebPny+v1yszU3R0dFDbe/XVV7V+/XolJSXJ5XLpjTfe6HJ5t9utYcOGKSkpqd1JhXv37pXH49H48eO1detWLVy4UJI0f/58ud1uPfzww8rIyJDP55Pb7dayZcuUkpLS6bbmz5+v++67j5MKhyBHV4ecJkyYYAcPHuzD4QBAe//5z3902223BbVOUZE0a5a0aZPkcl2mgXXizJkzGj58uCTpqaee0smTJ/WHP/yhbweBK15HvzcOh6PQzCZ0tDznEAAYklwu6ejR/tn29u3btXLlSvl8PsXFxenll1/un4EAQSAIAOAbNmvWLM2aNau/hwEEhXMIAAAAQQAAAAgCAAAgggAAAIggAIBBry+nP3700UdVXFzc5TJz5swJTLzUWm/G2VdaPibaE++8846ysrKUmJgor9erNWvWqKmpqcNlDx06pB07dgQ1lpycHI0cOTIwKVZHzEx5eXm6+eab5Xa7A39w6lIQBAAwyH2TQWBmar54MohW/vSnP2ncuHFBbavFQA6Cnlq7dq2eeeYZrVy5UkeOHNGePXtUV1enH//4xx3+KeneBMGcOXP097//vctldu7cqZKSEpWUlGjdunV67LHHgtpGRwgCAENSQYG0cqX/+lIN9emPy8rKdNttt2nBggW6/fbbdeLECT322GOaMGGCXC6XnnjiicDzpaWlqeUP1q1fv15jx45VWlqacnNz9Ytf/CKw3Pvvvx/UOCXp2Wef1R133CG32x3YZsvYcnNz5XK5NHXqVNXX10vy72xTUlLkdrs1ffr0wDwSrcdYXV2tMWPGSJLq6ur0ox/9SG63W7NmzdLEiRPV+o/v/frXv1ZSUpJSUlI6nDiqpKREmzdv1ltvvRX4N42MjNTSpUt16623tjsq0tDQoOXLl2vTpk3yeDyBn5vuTJ48Wddcc02Xy7zxxht65JFH5HA4lJKSotra2sDPX691Ng2iMf0xgAEimOmPU1PNPJ4LUyCHhPhvd7X8lT79cWlpqTkcDisoKAjc1zIdss/ns9TUVDt8+LCZmaWmptqBAwfss88+s7i4OKupqbGGhga766677Oc//3mvx7lr1y7Lzc215uZma2pqsszMTHvvvfestLTUQkND7d///reZmc2cOdNeeeUVMzNLTEy0vXv3mpnZsmXLbOH//iFbxmhmVlVVZXFxcWZm9uyzz9r8+fPNzOzIkSMWGhoaWE6Svfnmm2ZmtnjxYluxYkW7Mebn59vu3butqanJFixYYLfffrs98cQTlpeXZ19++aU98MAD7dZ56aWXAt8XM7N33nnHkpKS2l28Xm+b9UpLSwM/Nx3JzMy0ffv2BW7fc889gdfSgumPAVzxTp26MAVyc7P/dpDTCbRxJUx/HBcX12YOg82bN2vdunXy+Xw6efKkiouL5Xa7A49/+OGHSk1NDfxPdubMmTp+/Hivx7l7927t3r1b48ePl+T/888lJSUaPXq04uPjA+NOTk5WWVmZTp06pdraWqWmpkqSZs+e3eZ71pF//vOfgbkcEhIS2rye8PBw3X///YFt/OMf/2i3/uHDh5Wfn69t27YpLCxMhYWFeu6551RWVqaYmJgeTfqUnp7e6UyWwbAO3p5wOByX9JwEAYBBpbvpjyX/2wRTpkgNDVJ4uPTqq5LX2/ttDvXpjyX/oe8WpaWlWrVqlQ4cOKCYmBjNmTMnMB1yT7cd7DjNTPn5+frZz37W5v6ysrI2zxUaGhp4y6Azw4YNC5wH0XrcXY0jLCwssEMNDQ2Vz+frcIyhoaE6duyYMjIyJEn33XefPv74Y50/f77NODvz7rvv6le/+lW7+yMiIgKTY/VEbGysTpw4EbhdXl6u73znOz1evyOcQwBgyPF6pbffllas8F9fSgxIQ3/644udPn1akZGRio6OVkVFhXbu3NlumTvvvFPvvfeevvrqK/l8Pm3durXb19bVOO+9915t2LBBZ86ckSR99tlnqqys7PS5oqOjFRMTo3379kmSXnnllcDRgjFjxqiwsFCS2ryvf9ddd2nz5s2SpOLiYh05cqTbMbeWmJiogoIC3XLLLdq9e7ckadeuXTIzPf3005oxY0a3r7nlCMHFl2BiQJIeeOAB/eUvf5GZaf/+/YqOjtaoUaOCeo6LEQQAhiSvV8rPv/QYkIb+9McXS0pK0vjx4+VyuZSTk6Pvfe977Za54YYbtHTpUk2cOFHf//73NW7cuG5fZ1fjnDp1qh566CF5vV4lJiZqxowZ3R6C//Of/6zFixfL7Xbr0KFDWr58uSRp0aJFWrt2rSZNmqTq6urA8gsWLFBVVZXcbreefvppud3uoP5tZs+eraVLlyozM1P19fVKTk5WbW2tioqKNHz4cOXk5LRbJz09XcXFxUGdVPjggw/K6/Xqk08+UWxsrNb/b/7uF198US+++KIkadq0abrxxht18803Kzc3Vy+88EKPX0dnmP4YwIDXm+mP+9OVMv1xy+v0+XyaPn26cnJyNH369P4eVqeamprU2Ngop9OpTz/9VFOmTNHx48cVHh7e4+dYtWqVCgoKtHr1ao0ePVr19fV67bXXNHnyZH33u9+9jKMPHtMfA0A/u1KmP/7tb3+rPXv26Ny5c5o6daqys7P7e0hdqqurU3p6uhobG2VmWrt2bVAxIPmPPuzYsUO5ubmqrKxUdHS0HnzwQd1www2XadR9hyMEAAa8wXaEABgIgj1CwDkEAACAIAAwOPTko2sA/Hrz+0IQABjwnE6nampqiAKgB8xMNTU1cjqdQa3HSYUABrzY2FiVl5erqqqqv4cCDApOp1OxsbFBrUMQABjwwsLCFB8f39/DAIY03jIAAAAEAQAAIAgAAIAIAgAAIIIAAACIIAAAACIIAACACAIAACCCAAAAiCAAAAAiCAAAgAgCAH2sqEhKSPBfAxg4mNwIuIKcPStVVEhffOG/rqiQrrpKuvZaaeRI//W110qRkZLDcXm2P22adOKElJnpj4LIyG9+OwCCRxAAA0xBgbR3r5SWJnm93S9/7tyFnfsXX1zY2be+bvna55Ouv1667jr/9ciR0vnzUlXVhUtlpWR2IQ5ah0LrS+v7o6J6FhA5OReev6JCmjdP2rjxUr9jAL4JBAEwgBQUSFOmSA0NUliY9MIL/h1vZzv4L76Q6uv9y7Te0V9/vXTrrf6oaH1fT3fcZ8+2jYSWUKiqkj75pP39jY3dB8SHH0rbtvkDRvJfb9smbdjgDwUA/cthZp0+OGHCBDt48GAfDge4sq1cKS1bJjU1+W+PGiUlJbXd2bfewV93nRQTc3kO7wejvr59QLSOiKoqaedO/xGKi7UED4DLz+FwFJrZhI4e4wgBMICkpUnh4f4jBOHh0tatPXvboL9961vS6NH+S2c2bJDy8vxHH1pEREhPPXX5xwege3zKABhAvF7p7belFSv814MhBnoqJ8d/IqHT6b/tdEpZWdLcuf07LgB+vGUAoM+cPSuNG+f/lMHo0XzKAOhrXb1lwBECAH0mMlLascMfBdu3EwPAQMI5BAD6lMslHT3a36MAcDGOEAAAAIIAAAAQBAAAQAQBAAAQQQAAAEQQAAAAEQQAAEAEAQAAEEEAAABEEAAAABEEAABABAEAABBBAAAARBAAAAARBAAAQAQBAAAQQQAAAEQQAAAAEQQAAEAEAQAAEEEAAABEEAAAABEEAABABAEAABBBAAAARBAAAAARBAAAQAQBAAAQQQAAAEQQAAAAEQQAAEAEAQAAEEEAAABEEAAAABEEAABABAEAABBBAAAARBAAAAARBAAAQAQBAAAQQQAAAEQQAAAAEQQAAEAEAQAAEEEAAABEEAAAABEEAABABAEAABBBAAAARBAAAAARBAAAQAQBAAAQQQAAAEQQAAAAEQQAAEAEAQAAEEEAAABEEAAAABEEAABABAEAABBBAAAARBAAAAARBAAAQAQBAAAQQQAAAEQQAAAAEQQAAEAEAQAAEEEAAABEEAAAABEEAABABAEAABBBAAAARBAAAAARBAAAQAQBAAAQQQAAAEQQAAAAEQQAAEAEAQAAEEEAAABEEAAAABEEAABABAEAABBBAAAARBAAAAARBAAAQAQBAAAQQQAAAEQQAAAAEQQAAEAEAQAAEEEAAABEEAAAABEEAABABAEAABBBAAAARBAAAAARBAAAQAQBAAAQQQAAAEQQAAAAEQQAAEAEAQAAEEEAAABEEAAAABEEAABABAEAABBBAAAARBAAAAARBAAAQAQBAAAQQQAAAEQQAAAAEQQAAEAEAQAAEEEAAABEEAAAABEEAABABAEAABBBAAAARBAAAAARBAAAQAQBAAAQQQAAAEQQAAAAEQQAAEAEAQAAEEEAAABEEAAAABEEAABABAEAABBBAAAARBAAAAARBAAAQAQBAAAQQQAAAEQQAAAAEQQAAEAEAQAAEEEAAABEEAAAABEEAABABAEAABBBAAAARBAAAAARBAAAQAQBAAAQQQAAAEQQAAAAEQQAAEAEAQAAEEEAAABEEAAAABEEAABABAEAABBBAAAARBAAAAARBAAAQAQBAAAQQQAAAEQQAAAAEQQAAEAEAQAAEEEAAABEEAAAABEEAABABAEAABBBAAAARBAAAAARBAAAQAQBAAAQQQAAAEQQAAAAEQQAAEAEAQAAEEEAAABEEAAAABEEAABABAEAABBBAAAARBAAAAARBAAAQAQBAAAQQQAAAEQQAAAAEQQAAEAEAQAAEEEAAABEEAAAAEkOM+v8QYfjYB+OBQAAXF7VZpbR0QNdBgEAALgy8JYBAAAgCAAAAEEAAABEEAAAABEEAABA0v8DIIffuYlTOCcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 648x648 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgQAAAEJCAYAAADrdD/VAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3X98zvXi//HHzI+F+ZFDqbFGIbNr18yPTdguCTV8psNYCk3kR1F9zdBhnH7o4CA5h1sS6tZpRD45oSPZQuYYWrEl+8iSOGzaxmJs8/r+cbWrzX5jPzjP++22264fr+v9el3rdnU9va/39X46GWMQERGR/241qnoBIiIiUvUUCERERESBQERERBQIREREBAUCERERQYFAREREgJol3dmvXz+TmppaWWsRERGRCnTgwIF/GWP6FXVfiYEgNTWV/fv3V8yqREREpFI5OTn9obj79JGBiIiIKBCIiIiIAoGIiIigQCAiIiIoEIiIiAgKBCIiIoICgYiIiKBAICIiIigQiIiICAoEIiIiggKBiIiIoEAgIiIiKBCIiIgICgQiIiKCAoGIiIigQCAiIiIoEIiIiAgKBCIiIoICgYiIiKBAICIiIigQiIiICAoEIiIiggKBiIiIoEAgIiIiKBCIiIgICgQiIiKCAoGIiIigQCAiIiIoEIiIiAgKBCIiIoICgYiIiKBAICIiIigQiIiICAoEIiIiggKBiIiIoEAgIiIiKBCIiIgICgQiIiKCAoGIiIigQCAiIiIoEIiIiAgKBCIiIoICgYiIiKBAICIiIigQiIiICAoEIiIiggKBiIiIUEogOHAAmjaFyEjIzKysJYmIiEhlK3UPQWoqzJsHfn4KBSIiItXCvHkQHV3wtuhoeOyxom+fN6/UTZbpI4OsLDh2DObPL/NSyy0hATp0sP8WERGREnTuDCEhv7/5R0fbr/fuXfTtnTuXukknY0zxdzp1MrDfcb1JE0hJASenG3oahfz6K7RvDz/9BC1b2kNBvXo3dw4REZFb0rx59jd0m81+PToa4uJ+DwXjx8OyZbBunX3MbyEgd+x4WL6MH+au4wd3G2fPwogRTgeMMZ2KmqZmedZ07hzUqWM/rqBZM7jrrqJ/511u2hRq1y59u2FhcPYsGANnzsDo0RAVVZ6ViYiI3Abyv/lf88Zv1q7j4kVwGRHCd3PWkZRuo2WX8fi+8gpbO81k1TIbZ2bD2bM2wjLHE/76KyysN5M1f7M53p9LUq49BE2b2v8Vn5Jif+M+e/b33/kv5/1OSQFX14Ih4drfcXGweDFcuvT7vHXrwltv2YOCiIjI7Sr3jXlkPNCZtDS4+u84fmzWmW5vhhDrPZ4uB5bxunUdX1y10erHaKLO9gJgePMd/OBu4+Ea0UQcDGF/p/H4fb2MvS/Z9xB4JEfTMjwEp/HjcVqeb88B4ORU/B6CMgcCFxeYOhXmzCn7E716FdLSig4LeZe3bIHs7MKPbdbMfr+IiMjtKDcXHnWJZluO/Y3+//ns4FRbGyN/mEW/fa9wsP9Mfh77Z5o1g5bHomk+3D6OHTvsv0NCCn1MwPTpMHdu4dt/u37DgcDFBVq3hr17oX79m/sHefddmDTJfhxBnrp1YelSePrpmzuXiIhIdZK7PRrnR4p4o89/XEDebfkvP/44DBv2+3EFYH/znz8fwsML3x4XB1On3lggaNp0P+PH27d/s8NAnqFDYdMm+7cZXFzgf/5HxxCIiMh/gbxjBsD+pg2FjyHIfxsUeIMvr+sOBJ06dTL79+8v9v6bRd8yEBERqXglBYJqcerievXsxxK0bw+bNysMiIiIVLZyfe2wInl6wuHDVb0KERGR/07VYg+BiIiIVC0FAhEREVEgEBEREQUCERERQYFAREREUCAQERERFAhEREQEBQIRERFBgUBERERQIBAREREUCERERAQFAhEREUGBQERERFAgEBERERQIREREBAUCERERQYFAREREUCAQERERFAhEREQEBQIRERFBgUBERERQIBAREREUCERERAQFAhEREUGBQERERFAgEBERERQIREREBAUCERERQYFAREREUCAQERERFAhEREQEBQIRERFBgUBERERQIBAREREUCERERAQFAhEREUGBQERERFAgEBERERQIREREBAUCERERQYFAREREUCAQERERFAhEREQEBQIRERFBgUBERERQIBAREREUCERERAQFAhEREUGBQERERFAgEBERERQIREREBAUCERERQYFAREREUCAQERERFAhEREQEBQIRERFBgUBERERQIBAREREUCERERAQFAhEREUGBQERERFAgEBERERQIREREBAUCERERQYFAREREUCAQERERFAhEREQEBQIRERFBgUBERERQIBAREREUCERERAQFAhEREUGBQERERFAgEBERERQIREREBAUCERERQYFAREREgJrlfUB2djYnT54kKyurItYjcltycXHBzc2NWrVqVfVSRESKVO5AcPLkSVxdXbnvvvtwcnKqiDWJ3FaMMZw7d46TJ0/i4eFR1csRESlSuT8yyMrKokmTJgoDImXk5OREkyZNtFdNRKq16zqGQGFApHz0mhGR6q5SDipMSIAOHey/b3fJycn84x//KNdjXn/99euaa/HixVy8ePG6HisiIpJfhQeCX3+Fxx6DxEQICrJfv51VViDIzc1VIBARkZumwgNBWBicPQvGwJkzMHr0jW3v119/JSgoCG9vbzp06MDatWsBOHDgAAEBAfj6+tK3b19Onz4NQFxcHBaLBX9/f8LDw+nQoQMAq1evJjg4mAEDBuDh4cHSpUtZuHAhPj4++Pn58csvvwBw7Ngx+vXrh6+vLz169ODIkSMAjBo1ikmTJtGtWzdatWrF+vXrAZg2bRq7du3CarWyaNGiAms/ffo0PXv2xGq10qFDB3bt2sW0adO4dOkSVquV4cOHAxAcHIyvry+enp68/fbbjsfXr1+fWbNm0bVrV1577TVOnTqFzWbDZrPd2B9VRETEGFPsj6+vr7lWYmKi4/LkycYEBBT/06aNMTVqGGOPA/afGjXstxf3mMmTC01ZwPr1680zzzzjuJ6enm6uXLli/P39zdmzZ40xxkRFRZmnn37aGGOMp6en+eqrr4wxxkRERBhPT09jjDGrVq0yrVu3NufPnzdnz541DRo0MMuWLTPGGPPCCy+YRYsWGWOM6dWrlzl69Kgxxpi9e/cam81mjDFm5MiRZvDgwSY3N9ckJCSY1q1bG2OMiY6ONkFBQUWufcGCBebVV181xhiTk5Njzp8/b4wxpl69egXGnTt3zhhjzMWLF42np6dJTU01xhgDmLVr1zrGubu7m5SUlJL/YFJt5H/tiIhUBWC/KeY9v9xfOyyP48fh6tWCt129ar+9efPr26aXlxdTpkwhIiKC/v3706NHDw4fPszhw4d55JFHAPvu9ObNm5Oens6FCxfo1q0bAE888QSffvqpY1s2mw1XV1dcXV1p2LAhAwYMcMzx7bffkpmZyZ49exgyZIjjMZcvX3ZcDg4OpkaNGrRv354zZ86UuvbOnTsTFhZGdnY2wcHBWK3WIsctWbKEjRs3AvDTTz+RlJREkyZNcHZ25o9//GM5/2IiIiKlu6FAsHhxyfe/+y5MmlTwuIG6dWHpUnj66eubs02bNhw4cIAtW7Ywffp0+vTpw6BBg/D09CQ2NrbA2LS0tBK3VadOHcflGjVqOK7XqFGDnJwcrl69SqNGjYiPjy/18fbgVbKePXuyc+dONm/ezFNPPUV4eDgjRowoMCYmJobt27cTGxtL3bp1CQwMdHxdzcXFBWdn51LnERERKa8KPYYgLMx+IKGLi/26iwsMGHD9YQDg1KlT1K1blyeffJIpU6Zw8OBB2rZtS0pKiiMQZGdnk5CQQOPGjXF1dWXv3r0AREVFlWuuBg0a4OHhwUcffQTY3/S/+eabEh/j6urKhQsXirzvxx9/pFmzZowZM4bRo0dz8OBBAGrVqkV2djYAGRkZNG7cmLp163LkyBHH2ss7l4iISHlU+EGF774LzZqBkxPcdResXHlj2zt06BBdunTBarXy2muv8ac//YnatWuzfv16IiIi8Pb2xmq1smfPHgBWrlzJ2LFj8ff3xxhDw4YNyzXfBx98wMqVK/H29sbT05NPPvmkxPEWi4WaNWvi7e1d6KDCmJgYrFYrPj4+bNiwgcmTJwMwduxYLBYLw4cPp1+/fuTk5GCxWJg5cyZ+fn7FzjV27FgeffRRHVQoIiI3zKmkXd2dOnUy+/fvL3Dbd999x4MPPliuSRISYOhQWLsWPD2va53XLTMzk/r16wPwxhtvcPr0ad58883KXYQI1/faERG5mZycnA4YYzoVdV+FHlSYx9MTDh+ujJkK27x5M3PnziUnJwd3d3dWr15dNQsRERGpxiolEFSloUOHMnTo0KpehoiISLVWKacuFhERkepNgUBEREQUCERERESBQERERLhFA0FYWBjNmjVzFBUVxRjDpEmTuP/++7FYLI6TAN2o1atXc+rUqTKPj4mJcZwToTz279/PpEmTShyTnJxc7N+gvOusLLNnz2bBggVVvQwREbnGLRkIRo0axWeffVbimK1bt5KUlERSUhJvv/0248ePvylz38xAkJOTU+zjOnXqxJIlS8q9vjzVNRCIiEj1dEsGgp49e3LnnXeWOOaTTz5hxIgRODk54efnR3p6uqMSOU9ycjIPPvggY8aMwdPTkz59+nDp0iUA4uPj8fPzw2KxMGjQINLS0li/fj379+9n+PDhWK1Wx9g8S5YsoX379lgsFoYNG0ZycjLLly9n0aJFWK1Wdu3axahRo3jppZew2WxERESwb98+unXrho+PD926deP7778H7EGif//+AKSkpPDII4/QsWNHnn32Wdzd3UlNTQXsRU7Xrr+0dZa30tkY46iO9vLyclRO518jwHPPPec4z8OWLVto164d3bt3Z9KkSQXGJSYmEhgYSKtWrW4o9IiIyM1zY+cheOEFKKb4p4ATJ8DV1X75wgVo2bL4sVZr6a1JZfDzzz/TokULx3U3Nzd+/vlnml9Ts5iUlMSHH37IihUrCAkJYcOGDTz55JOMGDGCt956i4CAAGbNmsWcOXNYvHgxS5cuZcGCBXTqVPhET2+88QbHjx+nTp06pKen06hRI8aNG0f9+vWZMmUKYD+V8tGjR9m+fTvOzs6cP3+enTt3UrNmTbZv386MGTPYsGFDge3OmTOHXr16MX36dD777DPefvvtUtdf0jrHjh3L8uXLeeCBB/j3v//NhAkT2LFjBwCnT59m9+7dHDlyhIEDBzJ48GA+/vhj4uPj+eabb0hNTaVz58707Nmz2L99VlYWzz77LDt37sTDw4PQ0NAC9x85coTo6GguXLhA27ZtGT9+PLVq1Sp2eyIiUvEq58RErq7w7bf2yxZLpUxZ1CmZnZycCt3m4eHhqCH29fUlOTmZjIwM0tPTCQgIAGDkyJEFKpCLk9dHEBwcTHBwcLHjhgwZ4mgtzMjIYOTIkSQlJeHk5OQoOcpv9+7djjrkfv360bhx4xLXX5LrqXTevXs3oaGhODs7c9dddxEQEEBcXBwNGjQoco4jR47QqlUrPDw8AAgNDS0QYoKCgqhTpw516tShWbNmnDlzBjc3txLXLSIiFevGAkFZ/yUfHQ29ev3+mEoo43Fzc+Onn35yXD958iT33HNPoXH5K4ydnZ0L7V4vj82bN7Nz5042bdrEK6+8QkJCQpHj6tWr57g8c+ZMbDYbGzduJDk5mcDAwELjS+qbKO/6r6fSubj5a9asydWrVx3X82qaS6uCvnbNJR1LISIilaNyjiGIi4MdO+w/cXGVMuXAgQN57733MMawd+9eGjZsWOjjguI0bNiQxo0bs2vXLgDef/99x96C4iqHr169yk8//YTNZmPevHmkp6eTmZlZakVxRkYG9957L0CxPQvdu3dn3bp1AGzbto20tLRSn0Nx815PpXPPnj1Zu3Ytubm5pKSksHPnTrp06YK7uzuJiYlcvnyZjIwMvvjiCwDatWvHDz/84NhbkXfMgYiIVF+V85HB1Km/X74JewdCQ0OJiYkhNTUVNzc35syZw+jRo1m+fDkA48aN47HHHmPLli3cf//91K1bl1WrVpVrjjVr1jBu3DguXrxIq1atHI8fNWoU48aN44477iA2NpY77rgDsB/c9+STT5KRkYExhhdffJFGjRoxYMAABg8ezCeffMJbb71VaJ6pU6cycuRIFi5cSK+8vSjXiIyMJDQ0lLVr1xIQEEDz5s1xdXUlMzOz2PUXt06wVzqPHz+eV199lezsbIYNG4a3t3ex2xo0aBCxsbF4e3vj5OTEvHnzuPvuuwEICQnBYrHwwAMP4OPjA8Add9zB3//+d/r168cf/vAHunTpUspfW0REqlql1B/Ljbl8+TLOzs7UrFmT2NhYxo8fX+wu/+oir3baGMPEiRN54IEHePHFF6t6WVVKrx0RqWpVXn8sN+bEiROEhIRw9epVateuzYoVK6p6SaVasWIFa9as4cqVK/j4+PDss89W9ZJERKQE2kMgUkn02hGRqlbSHoJb8sREIiIicnMpEIiIiIgCgYiIiCgQiIiICLdoIPjss89o27Yt999/P2+88UaRY3788UcefvhhLBYLgYGBnDx58qbMrfrjG6P6YxGR6umWCwS5ublMnDiRrVu3kpiYyIcffkhiYmKhcVOmTGHEiBF8++23zJo1i+nTp9+U+VV/LCIit6OKDwTz5tm7DMD+e968G9rcvn37uP/++2nVqhW1a9dm2LBhfPLJJ4XGJSYm8vDDDwNgs9mKHKP6Y9Ufi4jIb4wxxf74+vqaayUmJv5+ZfJkYwICSv6xWIypVcuYli3tvy2WksdPnlxozvw++ugjM3r0aMf19957z0ycOLHQuNDQULN48WJjjDEbNmwwgElNTS0w5vjx48bZ2dl8/fXXxhhjhgwZYt5//31jjDFeXl4mJibGGGPMzJkzzeTf1hUQEGDi4uKKXFvz5s1NVlaWMcaYtLQ0Y4wxkZGRZv78+Y4xI0eONEFBQSYnJ8cYY0xGRobJzs42xhjz+eefm8cff9wYY0x0dLQJCgoyxhgzceJE8/rrrxtjjNm6dasBTEpKSonrL2mdvXr1MkePHjXGGLN3715js9kcaxs8eLDJzc01CQkJpnXr1sYYY9avX2969+5tcnJyzH/+8x/TokULc+rUqQJrzFvnqlWrzKVLl4ybm5v54YcfjDHGDBs2zDEuMjLS+Pv7m6ysLJOSkmLuvPNOc+XKlSLXebsp8NoREakCwH5TzHt+xZ+psHFjaN4cTpyAli3t12+AKWOt8YIFCxz/Yu3Zsyf33nsvNWsWfrqqP1b9sYiIVEb9cXQ0hITAzJmwbBlERt5QwVFZa43vuecePv74Y8D+JrhhwwYaNmxYaJzqj4vfnlH9sYjIf42KP4YgLg7WrYM//9n++wbrjzt37kxSUhLHjx/nypUrREVFMXDgwELjUlNTHW9Wc+fOJSwsrMxzqP64INUfi4jc/io+EEyd+vseAZutYBXydahZsyZLly6lb9++PPjgg4SEhODp6QnArFmz2LRpE2A/4K1t27a0adOGM2fO8PLLL5drnjVr1hAeHo7FYiE+Pp5Zs2YBv9cKX3uwXl79sZeXFz4+PgXqjzdu3Og4qPBaU6dOZfr06Tz00EPk5uYWuZbIyEi2bdtGx44d2bp1q6P+uCTFrRPs9ccrV67E29sbT0/PIg+4zG/QoEFYLBa8vb3p1auXo/64RYsWjvrj4cOHF1l/3L17d+66664i986IiEj1oXKjW4Dqj28Peu2ISFVT/fEtTvXHIiJS0bSHQKSS6LUjIlVN9cciIiJSIgUCERERUSAQERERBQIRERHhFg0EZak/PnHiBDabDR8fHywWC1u2bLkpcy9evJiLFy+Wefz//u//FtnGWJpNmzYV+9zyXFsulF9511lZRo0a5ShNEhGR6uOWCwRlrT9+9dVXCQkJ4euvvyYqKooJEybclPlvZiAo6ZS9AwcOZNq0aeVeX57qGghERKR6qtz6Y7jhCuSy1h87OTlx/vx5wH564KL6DmJiYggMDGTw4MG0a9eO4cOHO87D/8UXX+Dj44OXlxdhYWFcvnyZJUuWcOrUKWw2G7Yi+himTZvmqD+eMmUKe/bsYdOmTYSHh2O1Wjl27BiBgYHMmDGDgIAA3nzzTf75z3/StWtXfHx86N27t6NQaPXq1Tz33HOAva7Yz8+Pzp07M2vWLOrXr++YMzMzs9D6S1vngQMHCAgIwNfXl759+3L69GkAAgMDiYiIoEuXLrRp08ZxZsWsrCyefvppx1kYo3/775l/jQD9+/cnJiYGgJUrV9KmTRsCAwMZM2ZMgXE7d+4sVLEsIiJV68ZOTPTCC1DaGfPS0uBPf4K8719/9539cnG78K3WEkuTfv75Z1q0aOG47ubmxr///e9C42bPnk2fPn146623+PXXX9m+fXuR2/v6669JSEjgnnvu4aGHHuKrr76iU6dOjBo1ii+++II2bdowYsQIli1bxgsvvMDChQuJjo7mD3/4Q4Ht/PLLL2zcuJEjR47g5OREeno6jRo1YuDAgfTv35/Bgwc7xqanp/Pll1/+9udJY+/evTg5OfHOO+8wb948/vrXvxbY9uTJk5k8eTKhoaEsX7681PVPmjSp2HVmZ2fz/PPP88knn9C0aVPWrl3Lyy+/zLvvvgvY91rs27ePLVu2MGfOHLZv387f/vY3AA4dOsSRI0fo06cPR48eLfa/0alTp3jllVc4ePAgrq6u9OrVC29vb8f9p0+fZvfu3Rw5coSBAwcW+NuIiEjVqPg9BI0b2wPAt9/afx588IYqkMtaf/zhhx8yatQoTp48yZYtW3jqqacKNPPl6dKlC25ubtSoUQOr1UpycjLff/89Hh4etGnTBrDXH+/cubPEdTVo0AAXFxeeeeYZPv74Y+rWrVvs2KFDhzounzx5kr59++Ll5cX8+fOLbEiMjY111BU/8cQTpa6/JN9//z2HDx/mkUcewWq18uqrr3Ly5EnH/Y8//jhQsEp59+7dPPXUU4C9uMjd3b3EQLBv3z4CAgK48847qVWrVqHq6KIqlkVEpGpVfP0x2D8m6NXr98dUQv3xypUr+eyzzwDw9/cnKyuL1NRUmjVrVmBcUVW8pdX3FqVmzZrs27ePL774gqioKJYuXcqOHTuKHJu//vj555/npZdeYuDAgcTExDB79uxyzVveKmFjDJ6ensTGxpa4vfzbqsj64+v5W4uIyM1X8XsIoqMhJAR27LD/hIQUPKagnMpaf9yyZUtHHe93331HVlYWTZs2LdMc7dq1Izk5mf/7v/8DylZ/nJmZSUZGBo899hiLFy92lA+Vp/54zZo1RY7x8/Njw4YNAERFRZXpORQ3b9u2bUlJSXEEguzs7CL3SuTXs2dPPvjgAwCOHj3KiRMnaNu2Lffddx/x8fGO6ud9+/YB9r0WX375JWlpaeTk5DjWLiIi1VfFB4K4OFi3zr5XwGazX46Lu+7NlbX++K9//SsrVqzA29ub0NBQVq9eXeRHC0VxcXFh1apVDBkyBC8vL2rUqMG4ceMAGDt2LI8++mihg/UuXLhA//79sVgsBAQEsGjRIgCGDRvG/Pnz8fHx4dixY4Xmmj17NkOGDKFHjx6FPu/Ps3jxYhYuXEiXLl04ffp0maqEi1tn7dq1Wb9+PREREXh7e2O1WtmzZ0+J25owYQK5ubl4eXkxdOhQVq9eTZ06dXjooYfw8PDAy8uLKVOm0LFjRwDuvfdeZsyYQdeuXenduzft27dX/bGISDWncqNbwMWLF7njjjtwcnIiKiqKDz/8sMhvVlQnefXHOTk5DBo0iLCwMAYNGlTVy6pSeu2ISFVT/fEt7sCBAzz33HMYY2jUqJHjGwHV2ezZs9m+fTtZWVn06dOH4ODgql6SiIiUQIHgFtCjRw+++eabql5GuSxYsKCqlyAiIuVwy52pUERERG4+BQIRERFRIBAREREFAhEREeEWDQRlqT9+8cUXsVqtWK1W2rRpQ6NGjW7K3KtXr+bUqVNlHh8TE1Pq9/yLsn//fiZNmlTimOTkZDp06FDkfeVdZ2WZPXu2DjgUEamGKjYQXNt0CDfcdljW+uNFixYRHx9PfHw8zz//vOMc/TfqZgaCkk4z3KlTJ5YsWVLu9eWproFARESqp4oNBJ07FzxVcd5pjDt3vu5NlrX+OL8PP/yQ0NDQQrcnJyfz4IMPMmbMGDw9PenTpw+XLl0CID4+Hj8/PywWC4MGDSItLY3169ezf/9+hg8fjtVqdYzNs2TJEkf98bBhw0hOTmb58uUsWrQIq9XKrl27GDVqFC+99BI2m42IiAj27dtHt27d8PHxoVu3bnz//feAPUj0798fgJSUFB555BE6duzIs88+i7u7O6mpqYA9IF27/tLWeezYMfr164evry89evTgyJEjAIwaNYpJkyYVqiY2xhAeHk6HDh3w8vJi7dq1hdYI8Nxzz7F69WoAtmzZQrt27ejevTuTJk0qMC4xMZHAwEBatWp1Q6FHRERuImNMsT++vr7mWomJib9fmTzZmICAkn8sFmNq1TKmZUv7b4ul5PGTJxeaM7+PPvrIjB492nH9vffeMxMnTix2fHJysrn77rtNTk5OofuOHz9unJ2dzddff22MMWbIkCHm/fffN8YY4+XlZWJiYowxxsycOdNM/m1dAQEBJi4ursi5mjdvbrKysowxxqSlpRljjImMjDTz5893jBk5cqQJCgpyrCcjI8NkZ2cbY4z5/PPPzeOPP26MMSY6OtoEBQUZY4yZOHGief31140xxmzdutUAJiUlpcT1l7TOXr16maNHjxpjjNm7d6+x2WyOtQ0ePNjk5uaahIQE07p1a2OMMevXrze9e/c2OTk55j//+Y9p0aKFOXXqVIE15q1z1apV5tKlS8bNzc388MMPxhhjhg0b5hgXGRlp/P39TVZWlklJSTF33nmnuXLlSpHrvN0UeO2IiFQBYL8p5j2/4k9M1LgxNG8OJ05Ay5Y3VH0MZa8/zhMVFcXgwYNxdna3yMagAAAN2UlEQVQu8n4PDw+sVivwe+VvRkYG6enpjkKjkSNHFqrwLYrFYmH48OEEBweXeGa+IUOGONaTkZHByJEjSUpKwsnJiezs7ELjd+/ezcaNGwHo168fjfP9DYtaf0kyMzPZs2dPgedz+fJlx+Wiqol3795NaGgozs7O3HXXXQQEBBAXF0eDBg2KnOPIkSO0atUKDw8PAEJDQ3n77bcd9wcFBVGnTh3q1KlDs2bNOHPmDG5ubiWuW0REKlbF1x/nfUwwcyYsWwaRkZVSf5wnKiqKv/3tb8Xef2198LW718tj8+bN7Ny5k02bNvHKK68U2yKYv/545syZ2Gw2Nm7cSHJyMoGBgYXGFxWCrnf9V69epVGjRo42xpK2lzdvcfPfjPrjslQ2i4hIxavYYwjywsC6dfDnP9t/V1L9McD3339PWloa/v7+5ZqjYcOGNG7cmF27dgFlqz/OqwC22WzMmzeP9PR0MjMzy1V/nPf5+7W6d+/OunXrANi2bRtpaWmlPofi5m3QoAEeHh589NFHgP3Nu7TTIvfs2ZO1a9eSm5tLSkoKO3fupEuXLri7u5OYmMjly5fJyMhw1E23a9eOH374wbG3Iu+YAxERqb4qNhDkrz6GSq0/BvvBhMOGDStz7XF+a9asITw8HIvFQnx8PLNmzQLsB96NGzeu0MF6ubm5PPnkk3h5eeHj48OLL75Io0aNGDBgABs3bnQcVHitqVOnMn36dB566CFyc3OLXEtkZCTbtm2jY8eObN26lebNm+Pq6lri+otbJ8AHH3zAypUr8fb2xtPTs9SDMgcNGoTFYsHb25tevXoxb9487r77blq0aEFISIjjoxIfHx8A7rjjDv7+97/Tr18/unfvzl133aX6YxGRak71x7eAy5cv4+zsTM2aNYmNjWX8+PHF7vKvLvLqj40xTJw4kQceeIAXX3yxqpdVpfTaEZGqpvrjW9yJEycICQnh6tWr1K5dmxUrVlT1kkq1YsUK1qxZw5UrV/Dx8eHZZ5+t6iWJiEgJtIdApJLotSMiVa2kPQS35KmLRURE5OZSIBAREREFAhEREVEgEBEREW7RQFCW+mOAdevW0b59ezw9PXniiSduytyqP74xqj8WEamebtv646SkJObOnctXX31FQkICi8tymuUyUP2xiIjcjm7b+uMVK1YwceJERxFQs2bNCo1R/bHqj0VE5DfF1SCaW7z++H/+539MeHi46datm+natavZunVroTGqP1b9cWVS/bGIVDX+G+uPc3JySEpKIiYmhpMnT9KjRw8OHz5Mo0aNCoxT/bHqj0VE5DauP3Zzc8PPz49atWrh4eFB27ZtSUpKovM1H1eo/rj47RnVH4uI/Ne4beuPg4ODif5tntTUVI4ePUqrVq3KNIfqjwtS/bGIyO3vtq0/7tu3L02aNKF9+/bYbDbmz59PkyZNyjyP6o9/p/pjEZHbn8qNbgGqP7496LUjIlVN9ce3ONUfi4hIRdMeApFKoteOiFQ11R+LiIhIiRQIRERERIFAREREFAhERESECg4EmZn2ExM2bQo1ath/R0bab79dJScn849//KNcj3n99deva65nnnmmyKbH/EaNGuUoKcrvetZZWerXr1/msTt27GDAgAF4eXnh7+/P4sWLiz2fQ3x8PFu2bCnXWspStb169WqaNm2K1WrFarXyzjvvlGsOEZHqoMICQWYm+PnZm45TU8EY++958+y3366h4GYGAmNMgVMDX+udd96hffv25ZorT3UOBGW1bNky5s2bx9y5czl06BDbt2/n4sWLDBs2rMjTJ5c3EJS1ahtg6NChxMfHEx8fzzPPPHPdz0lEpKpUWCCYPx+OHYPfTm/vkJVlv33+/Ovb7q+//kpQUBDe3t506NDBcVrcAwcOEBAQgK+vL3379uX06dMAxMXFYbFY8Pf3d1T4gv1fdcHBwQwYMAAPDw+WLl3KwoUL8fHxwc/Pj19++QUof1XwtGnT2LVrF1arlUWLFhVY++nTp+nZsydWq5UOHTqwa9cupk2bxqVLl7BarQwfPtxRyTxhwgQ6duzITz/9xPjx4+nUqROenp5ERkY6thcYGEje10JXrlxJmzZtCAwMZMyYMTz33HOOcTt37izXOu3//ebTuXNnLBaLY87y1kVfu8bU1FTuu+8+AC5evOg4y+HQoUPp2rUr+b/i+vLLL+Pt7Y2fn5+jZCm/pKQk1q1bx6effur4b1qvXj1mzJhBu3btCu0VuXLlCrNmzWLt2rVYrdYynU65rFXbIiK3heJqEE0Z6o9Laj+uVcsY+36Bon9q1bqu9mOzfv1688wzzziup6enmytXrhh/f39z9uxZY4wxUVFR5umnnzbGGOPp6Wm++uorY4wxERERxtPT0xhjzKpVq0zr1q3N+fPnzdmzZ02DBg3MsmXLjDHGvPDCC2bRokXGmPJXBV9bCZzfggULzKuvvmqMMSYnJ8ecP3/eGGNMvXr1HGOOHz9unJycTGxsrOO2c+fOOR4TEBBgvvnmG2PM7xXHP//8s3F3dzfnzp0zV65cMd27d3dUQl/POv/1r3+ZMWPGmKtXr5rc3FwTFBRkvvzyyxuui05JSTHu7u7GGGPmz59vxo4da4wx5tChQ8bZ2dkxDjCbNm0yxhgTHh5uXnnllUJrnD59utm2bZvJzc01EyZMMB07djSRkZFm0qRJ5pdffjEDBw4s9JhVq1YVqMresWOH8fb2LvTj7+9vjCl71faqVavM3Xffbby8vMwf//hHc+LEiSL/rqo/FpGqRlXUHxfR4luu+4vj5eXFlClTiIiIoH///o5a48OHD/PII48A9l29zZs3Jz09nQsXLtCtWzcAnnjiCT799FPHtmw2G66urri6utKwYUMGDBjgmOPbb7+9rqrgknTu3JmwsDCys7MJDg521BZfy93dHT8/P8f1devW8fbbb5OTk8Pp06dJTEzEYrE47t+3bx8BAQHceeedgL1e+ejRo9e9zm3btrFt2zZHN0FmZiZJSUm0bNnyptVF7969m8mTJwPQoUOHAs+ndu3a9O/f3zHH559/Xujx33zzDdOnT+ef//wntWrV4sCBAyxcuJDk5GQaN25cYqFUHpvNVuIpoE0Zq7YHDBhAaGgoderUYfny5YwcOZIdO3aUOr+ISHVyQ4GgpPbjpk3txwyUdH9MTPnnbNOmDQcOHGDLli1Mnz6dPn36MGjQIDw9PYmNjS0wtrRWwPw1vDVq1HBcr1GjBjk5OddVFVySnj17snPnTjZv3sxTTz1FeHg4I0aMKDQufz3y8ePHWbBgAXFxcTRu3JhRo0Y5aobLOnd512mMYfr06YVON5ycnFzuuuX8Fcn5113SOmrVquV44y2uHtkYg7OzM0eOHKFfv34APProo3z77bdcvny5wDqLEx0dXWS/Qt26ddmzZ0+Zq7bzl2aNGTOGiIiIUucWEaluKuwYggkTwMWl6PtcXGD8+Ovb7qlTp6hbty5PPvkkU6ZM4eDBg7Rt25aUlBRHIMjOziYhIYHGjRvj6urK3r17AYiKiirXXNdTFVxS3fGPP/5Is2bNGDNmDKNHj+bgwYOA/Q0wu5hdJufPn6devXo0bNiQM2fOsHXr1kJjunTpwpdffklaWho5OTls2LCh1OdW0jr79u3Lu+++S+ZvR37+/PPPnD17tthtlVQXfd9993HgwAGAAp/r5690TkxM5NChQ6WuOT8vLy9iY2Np27Yt27ZtA+Bf//oXxhj+8pe/MHjw4FKfc94egmt/9uzZA5S9ajvveBWATZs26fTEInJLqrBAEB4OrVsXDgUuLvbbw8Ovb7uHDh2iS5cuWK1WXnvtNf70pz9Ru3Zt1q9fT0REBN7e3litVsf/1FeuXMnYsWPx9/fHGFPuGt7yVgVbLBZq1qyJt7d3oYP1YmJisFqt+Pj4sGHDBscu87FjxzoqhK/l7e2Nj48Pnp6ehIWF8dBDDxUac++99zJjxgy6du1K7969ad++fanPs6R19unThyeeeAJ/f3+8vLwYPHhwqbvgi6uLnjJlCsuWLaNbt26k5ttlNGHCBFJSUrBYLPzlL3/BYrGU67/NyJEjmTFjBkFBQVy6dAlfX1/S09NJSEigfv36hIWFFXqMzWYjMTGxzAcVlrVqe8mSJXh6euLt7c2SJUtYvXp1mZ+HiEh1UaHlRpmZ9m8TLFsG585Bkyb2PQPh4VCOr5rfkLwaXoA33niD06dP8+abb1bO5JUo73nm5OQwaNAgwsLCGDRoUFUvq1i5ublkZ2fj4uLCsWPHePjhhzl69Ci1a9cu8zYWLFhAbGwsixYtomXLlly6dImPP/6Ynj170qJFiwpc/fVRuZGIVLUqqz+uXx/mzLH/VJXNmzczd+5ccnJycHd3v23/9TZ79my2b99OVlYWffr0ITg4uKqXVKKLFy9is9nIzs7GGMOyZcvKFQbAvvdhy5YtjBkzhrNnz9KwYUNCQ0O59957K2jVIiK3L9Ufi1QSvXZEpKqp/lhERERKdF2BoCxfXROR3+k1IyLVXbkDgYuLC+fOndP/4ETKyBjDuXPncCnue7giItVAuQ8qdHNz4+TJk6SkpFTEekRuSy4uLri5uVX1MkREilXuQFCrVi08PDwqYi0iIiJSRXRQoYiIiCgQiIiIiAKBiIiIUMqJiZycnPYXe6eIiIjcalKNMf2KuqPEQCAiIiL/HfSRgYiIiCgQiIiIiAKBiIiIoEAgIiIiKBCIiIgI8P8Br3y82svj/LUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 648x648 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "visualise_iteration({'x': 20, 'y': 50}, {'x': -20, 'y': 50}, source, ccrs.PlateCarree())\n",
    "visualise_iteration({'x': 20, 'y': 50}, {'x': -20, 'y': 50}, source, ccrs.PlateCarree(180))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is is now obvious that one input line segment may produce multiple projected line segments. Some of these will be connected to one another, while others will be entirely disjoint.\n",
    "\n",
    "To identify the salient projected segments we perform successive binary searches. The cartopy implementation searches for the first t after t_min that would produce an appropriately straight segment. In pseudo code, one of these binary searches does the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bisecting from 0.15 to (0.15 .. 1.0)\n",
      "Bisecting from 0.15 to (0.15 .. 0.575)\n",
      "Bisecting from 0.15 to (0.3625 .. 0.575)\n",
      "Bisecting from 0.15 to (0.46875 .. 0.575)\n",
      "Bisecting from 0.15 to (0.46875 .. 0.521875)\n",
      "Bisecting from 0.15 to (0.4953125 .. 0.521875)\n",
      "Bisecting from 0.15 to (0.4953125 .. 0.50859375)\n",
      "Bisecting from 0.15 to (0.4953125 .. 0.501953125)\n",
      "Bisecting from 0.15 to (0.4986328125 .. 0.501953125)\n",
      "Bisecting from 0.15 to (0.4986328125 .. 0.50029296875)\n",
      "\n",
      "First cut t from 0.15 is in range (0.499462890625 .. 0.50029296875)\n"
     ]
    }
   ],
   "source": [
    "def bisect(line_start, line_end, t_start, interpolator, threshold, domain, DEBUG=True):\n",
    "    \"\"\"\n",
    "    The objective of this function is to find a range of t values which\n",
    "    contain a value of t that maximises the length of the straight line segment.\n",
    "    \"\"\"\n",
    "    # Track the min and max of t_current\n",
    "    t_min = t_start\n",
    "    t_max = 1.0\n",
    "    t_current = (t_min + t_max) * 0.5\n",
    "    \n",
    "    # Iterate until we have a particular confidence in t_current.\n",
    "    while t_max - t_min > 1e-3:\n",
    "        if DEBUG: print(f'Bisecting from {t_start} to ({t_min} .. {t_max})')\n",
    "        \n",
    "        # Determine if the current value of t_current produces a good\n",
    "        # line based on the straightness criteria.\n",
    "        good = TraceTesting.straight_and_within(\n",
    "            line_start, line_end, t_start, t_current,\n",
    "            interpolator, threshold, domain\n",
    "        )\n",
    "        # The line was good, so we could potentially grow it a\n",
    "        # little to get a longer straight line.\n",
    "        if good:\n",
    "            t_min = t_current\n",
    "        else:\n",
    "            t_max= t_current\n",
    "        t_current = (t_min + t_max) * 0.5\n",
    "\n",
    "    return t_min, t_max\n",
    "\n",
    "t_start = 0.15\n",
    "t_min, t_max = bisect(\n",
    "    {'x': 20, 'y': 50}, {'x': -20, 'y': 50},\n",
    "    t_start, interp, threshold, dest.domain)\n",
    "print()\n",
    "print(f'First cut t from {t_start} is in range ({t_min} .. {t_max})')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's implement the full segment splitting in pseudo code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.0, 0.5], [0.5009765625, 0.9990253448486328], []]\n"
     ]
    }
   ],
   "source": [
    "import shapely.geometry as sgeom\n",
    "segments = [[]]\n",
    "\n",
    "t_current = 0.0\n",
    "\n",
    "def point_inside(point, region):\n",
    "    p = TraceTesting.interp_t_pt(interp, line_start, line_end, t_min)\n",
    "    t_min_inside = dest.domain.contains(sgeom.Point([p['x'], p['y']]))\n",
    "\n",
    "\n",
    "line_start, line_end = {'x': 20, 'y': 50}, {'x': -20, 'y': 50}\n",
    "while t_current < 1:\n",
    "    t_min, t_max = bisect(\n",
    "        line_start, line_end,\n",
    "        t_current, interp, threshold, dest.domain, DEBUG=False)\n",
    "\n",
    "    p_min = TraceTesting.interp_t_pt(\n",
    "        interp, line_start, line_end, t_min)\n",
    "    t_min_inside = dest.domain.covers(\n",
    "        sgeom.Point([p_min['x'], p_min['y']]))\n",
    "    \n",
    "    p_current = TraceTesting.interp_t_pt(\n",
    "        interp, line_start, line_end, t_current)\n",
    "    t_current_inside = dest.domain.covers(\n",
    "        sgeom.Point([p_current['x'], p_current['y']]))\n",
    "\n",
    "    if t_min_inside:\n",
    "        if not segments[-1]:  # add_point if empty.\n",
    "            segments[-1].append(t_current)\n",
    "        if t_min == t_current:\n",
    "            # Our confidence range hasn't moved on. So we need to jump on to the next segment.\n",
    "            t_current = t_max\n",
    "            if t_current_inside:\n",
    "                segments.append([])\n",
    "        else:\n",
    "            # We have refined a new segment, so let's add it and move on.\n",
    "            segments[-1].append(t_min)\n",
    "            t_current = t_min\n",
    "    else:  # NB: in the actual implementation we also handle NANs.\n",
    "        if t_min != t_current:\n",
    "            t_current = t_min\n",
    "        else:\n",
    "            t_current = t_max\n",
    "            if t_current_inside:\n",
    "                lines.new_line()\n",
    "\n",
    "print(segments)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}