three_arm_survey.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Three-arm Antenna Position Calibration\n",
    "\n",
    "Author: Tim Molteno. tim@elec.ac.nz.\n",
    "\n",
    "The antennas are laid out on three arms which are placed approximately 120 degrees apart.  This notebook will use measurements to infer the actual layout of the antennas.\n",
    "\n",
    "Arm 1 will point geographically north (NOT magnetic). This arm forms the y-axis and will be used as a reference. We will assume that the x positions along this arm are approximately zero.\n",
    "\n",
    "The other reference point is the center of the array. This is assumed to be the point [0,0]. Make sure that a marker is present at this reference point at the same height as the level of the antennas (otherwise the distances will be wrong)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import minimize "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = [0,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "A 0      210.0\n",
       "A 1      410.0\n",
       "A 2      384.0\n",
       "A 3      638.0\n",
       "A 4      212.0\n",
       "A 5      753.0\n",
       "A 6     1261.0\n",
       "A 7     1436.0\n",
       "A 8     1707.0\n",
       "A 9     1864.0\n",
       "A 10    2199.0\n",
       "A 11    2386.0\n",
       "A 12     937.0\n",
       "A 13    1084.0\n",
       "A 14    1621.0\n",
       "A 15    1924.0\n",
       "A 16    2091.0\n",
       "A 17    2408.0\n",
       "A 18     914.0\n",
       "A 19    1254.0\n",
       "A 20    1575.0\n",
       "A 21    1696.0\n",
       "A 22    1856.0\n",
       "A 23    2076.0\n",
       "Name: 0, dtype: float64"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# pip install pandas\n",
    "# pip install odfpy\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "#pd.read_excel('three_arm_measurements.ods', engine='odf')\n",
    "cols = [f\"A {i}\" for i in range(24)]\n",
    "data = pd.read_excel('three_arm_measurements.ods', \"Sheet1\", usecols=cols)\n",
    "radius = data.loc[0]\n",
    "radius"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finding the antennas\n",
    "\n",
    "This is done by measuring the distance from each antenna to x0 as well as each other antenna. The distance $m_0$ are the measurements to $x_0$. The array m_ij are the distances between antennas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_ant = 24\n",
    "m_0 = np.zeros(24)\n",
    "m_ij = np.zeros((24,24))\n",
    "\n",
    "radius = data.loc[0]\n",
    "\n",
    "for i in range(24):\n",
    "    m_ij[i,:] = data.loc[i+1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following are the measured distances from [x1, x0, x2] from the reference points in millimeters. Note that their order must be the same as the order of the variable called 'reference_points'. In this case, they are x1,x0,x2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[   0.,  198.,  509.,  756.,  365.,  872., 1054., 1225., 1489.,\n",
       "        1656., 1984., 2183., 1046., 1189., 1713., 2023., 2187., 2506.,\n",
       "        1027., 1372., 1691., 1814., 1969., 2192.],\n",
       "       [ 198.,    0.,  673.,  894.,  545., 1021.,  855., 1027., 1293.,\n",
       "        1459., 1790., 1981., 1169., 1313., 1823., 2127., 2293., 2613.,\n",
       "        1173., 1506., 1823., 1945., 2098., 2316.],\n",
       "       [ 509.,  673.,    0.,  258.,  524., 1001., 1467., 1634., 1893.,\n",
       "        2056., 2384., 2575.,  557.,  701., 1229., 1538., 1705., 2026.,\n",
       "        1152., 1488., 1807., 1928., 2082., 2305.],\n",
       "       [ 756.,  894.,  258.,    0.,  772., 1214., 1639., 1802., 2054.,\n",
       "        2213., 2534., 2721.,  299.,  444.,  973., 1284., 1447., 1768.,\n",
       "        1357., 1682., 1992., 2113., 2261., 2479.],\n",
       "       [ 365.,  545.,  524.,  772.,    0.,  535., 1385., 1558., 1819.,\n",
       "        1984., 2312., 2507., 1066., 1209., 1735., 2037., 2205., 2525.,\n",
       "         697., 1042., 1366., 1490., 1644., 1867.],\n",
       "       [ 872., 1021., 1001., 1214.,  535.,    0., 1775., 1938., 2187.,\n",
       "        2348., 2668., 2859., 1477., 1608., 2104., 2405., 2562., 2873.,\n",
       "         159.,  507.,  831.,  956., 1109., 1333.],\n",
       "       [1054.,  855., 1467., 1639., 1385., 1775.,    0.,  172.,  437.,\n",
       "         600.,  931., 1127., 1865., 1982., 2435., 2710., 2860., 3157.,\n",
       "        1905., 2199., 2484., 2598., 2740., 2946.],\n",
       "       [1225., 1027., 1634., 1802., 1558., 1938.,  172.,    0.,  264.,\n",
       "         432.,  758.,  955., 2016., 2133., 2574., 2844., 2991., 3285.,\n",
       "        2065., 2355., 2636., 2745., 2887., 3088.],\n",
       "       [1489., 1293., 1893., 2054., 1819., 2187.,  437.,  264.,    0.,\n",
       "         164.,  492.,  692., 2262., 2364., 2796., 3058., 3198., 3484.,\n",
       "        2314., 2593., 2867., 2971., 3112., 3307.],\n",
       "       [1656., 1459., 2056., 2213., 1984., 2348.,  600.,  432.,  164.,\n",
       "           0.,  328.,  523., 2415., 2522., 2935., 3195., 3335., 3613.,\n",
       "        2469., 2742., 3013., 3118., 3254., 3447.],\n",
       "       [1984., 1790., 2384., 2534., 2312., 2668.,  931.,  758.,  492.,\n",
       "         328.,    0.,  193., 2726., 2827., 3227., 3474., 3609., 3884.,\n",
       "        2785., 3047., 3307., 3414., 3548., 3737.],\n",
       "       [2183., 1981., 2575., 2721., 2507., 2859., 1127.,  955.,  692.,\n",
       "         523.,  193.,    0., 2913., 3011., 3404., 3647., 3778., 4048.,\n",
       "        2976., 3236., 3493., 3595., 3723., 3913.],\n",
       "       [1046., 1169.,  557.,  299., 1066., 1477., 1865., 2016., 2262.,\n",
       "        2415., 2726., 2913.,    0.,  143.,  673.,  978., 1144., 1463.,\n",
       "        1614., 1924., 2223., 2338., 2485., 2697.],\n",
       "       [1189., 1313.,  701.,  444., 1209., 1608., 1982., 2133., 2364.,\n",
       "        2522., 2827., 3011.,  143.,    0.,  529.,  837.,  999., 1321.,\n",
       "        1743., 2048., 2342., 2467., 2498., 2807.],\n",
       "       [1713., 1823., 1229.,  973., 1735., 2104., 2435., 2574., 2796.,\n",
       "        2935., 3227., 3404.,  673.,  529.,    0.,  305.,  468.,  789.,\n",
       "        2233., 2519., 2799., 2906., 3043., 3245.],\n",
       "       [2023., 2127., 1538., 1284., 2037., 2405., 2710., 2844., 3058.,\n",
       "        3195., 3474., 3647.,  978.,  837.,  305.,    0.,  162.,  483.,\n",
       "        2523., 2798., 3068., 3177., 3308., 3506.],\n",
       "       [2187., 2293., 1705., 1447., 2205., 2562., 2860., 2991., 3198.,\n",
       "        3335., 3609., 3778., 1144.,  999.,  468.,  162.,    0.,   nan,\n",
       "          nan,   nan,   nan,   nan,   nan,   nan],\n",
       "       [2506., 2613., 2026., 1768., 2525., 2873., 3157., 3285., 3484.,\n",
       "        3613., 3884., 4048., 1463., 1321.,  789.,  483.,   nan,    0.,\n",
       "          nan,   nan,   nan,   nan,   nan,   nan],\n",
       "       [1027., 1173., 1152., 1357.,  697.,  159., 1905., 2065., 2314.,\n",
       "        2469., 2785., 2976., 1614., 1743., 2233., 2523.,   nan,   nan,\n",
       "           0.,   nan,   nan,   nan,   nan,   nan],\n",
       "       [1372., 1506., 1488., 1682., 1042.,  507., 2199., 2355., 2593.,\n",
       "        2742., 3047., 3236., 1924., 2048., 2519., 2798.,   nan,   nan,\n",
       "          nan,    0.,   nan,   nan,   nan,   nan],\n",
       "       [1691., 1823., 1807., 1992., 1366.,  831., 2484., 2636., 2867.,\n",
       "        3013., 3307., 3493., 2223., 2342., 2799., 3068.,   nan,   nan,\n",
       "          nan,   nan,    0.,   nan,   nan,   nan],\n",
       "       [1814., 1945., 1928., 2113., 1490.,  956., 2598., 2745., 2971.,\n",
       "        3118., 3414., 3595., 2338., 2467., 2906., 3177.,   nan,   nan,\n",
       "          nan,   nan,   nan,    0.,   nan,   nan],\n",
       "       [1969., 2098., 2082., 2261., 1644., 1109., 2740., 2887., 3112.,\n",
       "        3254., 3548., 3723., 2485., 2498., 3043., 3308.,   nan,   nan,\n",
       "          nan,   nan,   nan,   nan,    0.,   nan],\n",
       "       [2192., 2316., 2305., 2479., 1867., 1333., 2946., 3088., 3307.,\n",
       "        3447., 3737., 3913., 2697., 2807., 3245., 3506.,   nan,   nan,\n",
       "          nan,   nan,   nan,   nan,   nan,    0.]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for i in range(n_ant):\n",
    "    for j in range(n_ant):\n",
    "        if not np.isnan(m_ij[i,j]):\n",
    "            m_ij[j,i] = m_ij[i,j]\n",
    "\n",
    "m_ij\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the Initial Guess Points\n",
    "\n",
    "Initial Guesses are from JSON queried from the telescope API. These are converted to millimeters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import requests\n",
    "import json\n",
    "\n",
    "def get_data(path):\n",
    "    server = \"https://tart.elec.ac.nz/rhodes\"\n",
    "\n",
    "    r = requests.get('{}/{}'.format(server, path))\n",
    "    return json.loads(r.text)\n",
    "\n",
    "def get_pos():\n",
    "    return np.array(get_data('api/v1/imaging/antenna_positions'))\n",
    "\n",
    "if False:\n",
    "    with open('original_rhodes_positions.json', \"r\") as f:\n",
    "        data = f.read()\n",
    "    cal_json = json.loads(data)\n",
    "    current_pos = cal_json['antenna_positions'] # get_pos()\n",
    "    current_pos = np.array(current_pos)\n",
    "\n",
    "if False:\n",
    "    current_pos = get_pos()\n",
    "    json_result = {}\n",
    "    json_result[\"antenna_positions\"] = current_pos.tolist()\n",
    "    print(json.dumps(json_result, indent=4, separators=(',', ': ')))\n",
    "\n",
    "current_pos = np.array([[ 0.        ,  0.21036822,  0.        ],\n",
    "       [ 0.        ,  0.40894464,  0.        ],\n",
    "       [ 0.33652171, -0.1942909 ,  0.        ],\n",
    "       [ 0.55195215, -0.31866972,  0.        ],\n",
    "       [-0.18186395, -0.1049992 ,  0.        ],\n",
    "       [-0.65058666, -0.37561638,  0.        ],\n",
    "       [ 0.        ,  1.26747133,  0.        ],\n",
    "       [ 0.        ,  1.44403899,  0.        ],\n",
    "       [ 0.        ,  1.71019648,  0.        ],\n",
    "       [ 0.        ,  1.86051529,  0.        ],\n",
    "       [ 0.        ,  2.20143587,  0.        ],\n",
    "       [ 0.        ,  2.38582938,  0.        ],\n",
    "       [ 0.81594596, -0.47108662,  0.        ],\n",
    "       [ 0.94307482, -0.5444845 ,  0.        ],\n",
    "       [ 1.40247705, -0.8097205 ,  0.        ],\n",
    "       [ 1.66872406, -0.96343829,  0.        ],\n",
    "       [ 1.80943056, -1.04467522,  0.        ],\n",
    "       [ 2.08400692, -1.20320196,  0.        ],\n",
    "       [-0.78741312, -0.45461318,  0.        ],\n",
    "       [-1.08399029, -0.62584208,  0.        ],\n",
    "       [-1.36391741, -0.78745808,  0.        ],\n",
    "       [-1.46598234, -0.8463853 ,  0.        ],\n",
    "       [-1.60875644, -0.92881597,  0.        ],\n",
    "       [-1.80004538, -1.03925669,  0.        ]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0 210.36822\n",
      "0.0 408.94464000000005\n",
      "336.52171 -194.2909\n",
      "551.95215 -318.66972\n",
      "-181.86395 -104.9992\n",
      "-650.5866599999999 -375.61638000000005\n",
      "0.0 1267.47133\n",
      "0.0 1444.03899\n",
      "0.0 1710.19648\n",
      "0.0 1860.5152899999998\n",
      "0.0 2201.4358700000003\n",
      "0.0 2385.82938\n",
      "815.94596 -471.08662000000004\n",
      "943.0748199999999 -544.4845\n",
      "1402.47705 -809.7204999999999\n",
      "1668.72406 -963.43829\n",
      "1809.43056 -1044.67522\n",
      "2084.0069200000003 -1203.2019599999999\n",
      "-787.4131199999999 -454.61318\n",
      "-1083.99029 -625.84208\n",
      "-1363.91741 -787.45808\n",
      "-1465.98234 -846.3853\n",
      "-1608.75644 -928.81597\n",
      "-1800.04538 -1039.25669\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuIAAALJCAYAAADxiEbKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABI4ElEQVR4nO3df5yVdZ3//8eLH+oU4kgKwYyGlY0IJgir7Mci1HBMXUUsV5dNDV2zT/mtLAw/fvdb7C6Jmql9NFtdNVor0kIkMzDFWXdN119D4Y8lLUkY8Fc6CToiDO/vH3OYzsCAM0dm3jNnHvfb7dzmnNd1nTOv82pmenrxPtcVKSUkSZIkda9+uRuQJEmS+iKDuCRJkpSBQVySJEnKwCAuSZIkZWAQlyRJkjIwiEuSJEkZGMQlqReIiF9GxBk72P69iPjHDr5WXUScvfO6e9vv99GIWLGD7ftGxPqI6N9dPUlSTxCeR1yS8oiIlcDZKaW7O/m8MwvP+0iJ37cOuDml9G+lPP+dKvV9S1K58Yi4JEmSlIFBXJJ6gIg4MyL+KyK+FRGvRsSzEfGJou11EXF2RIwCvgf8dWE5R2Nh+/cj4l8K9/eMiDsi4qXCa90REdUd7OMbEfHTiPhJRKyLiMci4uCi7aMKvTRGxBMRcULRtmMj4snC8xoi4quF+uSIWF24/+/AvsDPC/1fEBEjIyJFxIDCPiMiYlFEvBIRz0TEP2zV3y0R8YPC93kiIiYUbf9a4Xuvi4gVEXFUCf9zSFK3MIhLUs9xGLAC2Au4FLghIqJ4h5TSU8C5wAMppUEppcp2XqcfcBPwPlpCbxNwdSf6OBG4FRgC/AhYGBEDI2Ig8HPgLmAocB7ww4ioKTzvBuCzKaXdgTHA0q1fOKX0aeA54G8K/V/azvf/MbAaGAF8EvjmVoH6BGA+UAks2vLeCn18AfirQg+1wMpOvG9J6lYGcUnqOf6YUro+pdQMzAOGA8M6+yIppT+llH6WUnojpbQOmAN8rBMv8WhK6acppY3At4HdgImF2yBgbkrprZTSUuAO4LTC8zYCB0bE4JTSqymlxzrbe0TsA3wE+FpK6c2U0jLg34BPF+32XymlOwtz+ndgyxH7ZmDXQg8DU0orU0q/72wPktRdDOKS1HM8v+VOSumNwt1BnX2RiHhXRPxrRPwxIl4D7gMqO3FWklVFfWzmL0enRwCrCrUt/ghUFe6fDBwL/DEi/iMi/rqzvRe+xyuF/4Bo73tA0ZyAN4DdImJASukZ4EvAN4AXI2J+RIwooQdJ6hYGcUnqfd7udFdfAWqAw1JKg4FJhXps/ylt7LPlTkT0A6qBNYXbPoXaFvsCDQAppYdTSifSsmxlIXBLCf2vAYZExO7tfY+3k1L6UeFsMu8rfJ9LOvI8ScrBIC5Jvc8LQHVE7LKd7bvTsi68MSKGAF/v5OuPj4hphQ9PfgnYADwI/DfwOnBBYc34ZOBvgPkRsUtETI+IPQpLWl6jZanI9vp/f3sbUkqrgF8DF0fEbhHxYeAs4Idv13RE1ETEkRGxK/AmLTPYXg+SlJ1BXJJ6n6XAE8DzEfFyO9uvBCqAl2kJ0Is7+fq3A38LvErL2uxpKaWNKaW3aPmg5CcKr/1d4PSU0v8UnvdpYGVhOcy5wN9v5/UvBv7fwplXvtrO9tOAkbQcHb8N+HpK6Vcd6HtXYG6ht+dpOTL/fzrwPEnKwgv6SJJaRcQ3gA+mlLYXoiVJO4lHxCVJkqQMDOKSJElSBi5NkSRJkjLwiLgkSZKUwYDcDeSy1157pZEjR+Zuo9d5/fXXefe73527jV7HuZXGuZXO2ZXGuZXGuZXO2ZWmt83t0UcffTmltPfW9T4bxEeOHMkjjzySu41ep66ujsmTJ+duo9dxbqVxbqVzdqVxbqVxbqVzdqXpbXOLiD+2V3dpiiRJkpSBQVySJEnKwCAuSZIkZWAQlyRJkjIwiEuSJEkZGMQlSZKkDLIF8YjYJyLujYinIuKJiPhiof6NiGiIiGWF27FFz7kwIp6JiBURUVtUHx8RywvbvhMRkeM9SZIkSR2V8zzim4CvpJQei4jdgUcj4leFbVeklL5VvHNEHAicCowGRgB3R8SHUkrNwLXAOcCDwJ3AMcAvu+l9SJIkSZ2W7Yh4SmltSumxwv11wFNA1Q6eciIwP6W0IaX0LPAMcGhEDAcGp5QeSCkl4AfA1K7tXpIkSXpnoiW7Zm4iYiRwHzAGOB84E3gNeISWo+avRsTVwIMppZsLz7mBlqPeK4G5KaWPF+ofBb6WUjq+ne9zDi1Hzhk2bNj4+fPnd+0bK0Pr169n0KBBudvodZxbaZxb6ZxdaZxbaZxb6ZxdaXrb3I444ohHU0oTtq5nv8R9RAwCfgZ8KaX0WkRcC/wzkApfLwdmAO2t+047qG9bTOk64DqACRMmpN50adSeorddUrancG6lcW6lc3alcW6lcW6lc3alKZe5ZT1rSkQMpCWE/zCltAAgpfRCSqk5pbQZuB44tLD7amCfoqdXA2sK9ep26pIkSVKPlfOsKQHcADyVUvp2UX140W4nAY8X7i8CTo2IXSNiP2B/4KGU0lpgXURMLLzm6cDt3fImJEmSpBLlXJpyOPBpYHlELCvU/g9wWkSMpWV5yUrgswAppSci4hbgSVrOuPL5whlTAD4HfB+ooGXduGdMkSRJUo+WLYinlP6L9td337mD58wB5rRTf4SWD3pKkiRJvYJX1pQkSZIyMIhLkiRJGRjEJUmSpAwM4pIkSVIGBnFJkiQpA4O4JEmSlIFBXJLK1IwZMxg6dChjxvzl7K633noro0ePpl+/fjzyyCMZu5MkGcQlqUydeeaZLF68uE1tzJgxLFiwgEmTJmXqSpK0Rc4ra0qSutCkSZNYuXJlm9qoUaPyNCNJ2oZHxCVJkqQMDOKSVEYW1jdw+NylLG/4M4fPXcpdTzyfuyVJ0na4NEWSysTC+gYuXLCcpo3NsA80NDZxyeJVbHxzU+7WJEntMIhLUpm4bMmKlhBeZMOmZl5ZvyFTR5KkHTGIS1KZWNPY1ObxS4suZcNzy2lueo3q6mpmz57NkCFDOO+883jppZc47rjjGDt2LEuWLMnUsST1bQZxSSoTIyoraCgK43ufcAEAVZUV3D/ryNb6SSed1O29SZK25Yc1JalMzKytoWJg/za1ioH9mVlbk6kjSdKOeERcksrE1HFVQMtacVhHVWUFM2trWuuSpJ7FIC5JZWTquCqmjquirq6O86ZPzt2OJGkHXJoiSZIkZWAQlyRJkjIwiEuSJEkZGMQlSZKkDAzikiRJUgYGcUmSJCkDg7gkSZKUgUFckiRJysAgLkmSJGVgEJckSZIyMIhLkiRJGRjEJUmSpAwM4pIkSVIGBnFJkiQpA4O4JEmSlIFBXJIkScrAIC5JkiRlYBCXJEmSMjCIS1KZmjFjBkOHDmXMmDGttVdeeYUpU6aw//77M2XKFF599dWMHUpS32YQl6QydeaZZ7J48eI2tblz53LUUUfx9NNPc9RRRzF37txM3UmSDOKSVKYmTZrEkCFD2tRuv/12zjjjDADOOOMMFi5cmKEzSRIYxCWpT3nhhRcYPnw4AMOHD+fFF1/M3JEk9V0GcUmSJCkDg7gklZGF9Q0cPncpyxv+zOFzl3LXE8+32T5s2DDWrl0LwNq1axk6dGiONiVJGMQlqWwsrG/gwgXLaWhsAqChsYlLFq/gtTc3te5zwgknMG/ePADmzZvHiSeemKVXSRIMyN2AJGnnuGzJCpo2Nrc+fmnRpWx4bjmbm16jurqa2bNnM2vWLE455RRuuOEG9t13X2699daMHUtS32YQl6QysaZwJHyLvU+4AIAAnp17XGv9nnvu6c62JEnb4dIUSSoTIyorOlWXJOVlEJekMjGztoaKgf3b1CoG9mdmbU2mjiRJO+LSFEkqE1PHVQEta8VhHVWVFcysrWmtS5J6FoO4JJWRqeOqmDquirq6Os6bPjl3O5KkHXBpiiRJkpSBQVySJEnKwCAuSZIkZWAQlyRJkjIwiEuSJEkZGMQlSZKkDAzikiRJUgYGcUmSJCkDg7gkSZKUgUFckiRJysAgLkmSJGVgEJckSZIyMIhLkiRJGRjEJUmSpAwM4pIkSVIGBnFJkiQpA4O4JEmSlIFBXJIkScrAIC5JkiRlYBCXJEmSMjCIS5IkSRkYxCVJkqQMDOKSJElSBgZxSZIkKQODuCRJkpRBtiAeEftExL0R8VREPBERXyzUh0TEryLi6cLXPYuec2FEPBMRKyKitqg+PiKWF7Z9JyIix3uSJEmSOirnEfFNwFdSSqOAicDnI+JAYBZwT0ppf+CewmMK204FRgPHAN+NiP6F17oWOAfYv3A7pjvfiCRJktRZ2YJ4SmltSumxwv11wFNAFXAiMK+w2zxgauH+icD8lNKGlNKzwDPAoRExHBicUnogpZSAHxQ9R5IkSeqRoiW7Zm4iYiRwHzAGeC6lVFm07dWU0p4RcTXwYErp5kL9BuCXwEpgbkrp44X6R4GvpZSOb+f7nEPLkXOGDRs2fv78+V35tsrS+vXrGTRoUO42eh3nVhrnVjpnVxrnVhrnVjpnV5reNrcjjjji0ZTShK3rA3I0UywiBgE/A76UUnptB8u729uQdlDftpjSdcB1ABMmTEiTJ0/udL99XV1dHc6t85xbaZxb6ZxdaZxbaZxb6ZxdacplblnPmhIRA2kJ4T9MKS0olF8oLDeh8PXFQn01sE/R06uBNYV6dTt1SZIkqcfKedaUAG4Ankopfbto0yLgjML9M4Dbi+qnRsSuEbEfLR/KfCiltBZYFxETC695etFzJEmSpB4p59KUw4FPA8sjYlmh9n+AucAtEXEW8BzwKYCU0hMRcQvwJC1nXPl8Sqm58LzPAd8HKmhZN/7LbnoPkiRJUkmyBfGU0n/R/vpugKO285w5wJx26o/Q8kFPSZIkqVfwypqSJElSBgZxSZIkKQODuCRJkpSBQVySJEnKwCAuSZIkZWAQlyRJkjIwiEuSJEkZGMQlSZKkDAzikiRJUgYGcUmSJCkDg7gklamrrrqKMWPGMHr0aK688src7UiStmIQl6Qy9Oyzz3L99dfz0EMP8Zvf/IY77riDp59+OndbkqQiBnFJKkN//OMfmThxIu9617sYMGAAH/vYx7jttttytyVJKmIQl6QytN9++3Hffffxpz/9iTfeeIM777yTVatW5W5LklRkQO4GJEk7z8L6Bi5bsoJT96lk80EncOhHJlO1954cfPDBDBjgn3xJ6kn8qyxJZWJhfQMXLlhO08Zm2Ac27X8EFQd+nPOnHcRDt15DdXV17hYlSUUM4pJUJi5bsqIlhBc0v95I07sr+Zef3Mf6hQt44IEHMnYnSdqaQVySysSaxqY2j19a+E02N61jTb/+/PLH/8aee+6ZqTNJUnsM4pJUJkZUVtBQFMbfO/1SAKoqKzjqqCNztSVJ2g7PmiJJZWJmbQ0VA/u3qVUM7M/M2ppMHUmSdsQj4pJUJqaOqwJa1orDOqoqK5hZW9NalyT1LAZxSSojU8dVMXVcFXV1dZw3fXLudiRJO+DSFEmSJCkDg7gkSZKUgUFckiRJysAgLkmSJGVgEJckSZIyMIhLkiRJGRjEJUmSpAwM4pIkSVIGBnFJkiQpA4O4JEmSlIFBXJIkScrAIC5JkiRlYBCXpDJ1xRVXMHr0aMaMGcNpp53Gm2++mbslSVIRg7gklaGXXnqJ73znOzzyyCM8/vjjNDc3M3/+/NxtSZKKGMQlqUxt2rSJpqYmNm3axBtvvMGIESNytyRJKmIQl6QytPfee/PVr36Vfffdl+HDh7PHHntw9NFH525LklTEIC5JZWjdunXcfvvtPPvss6xZs4bXX3+dm2++OXdbkqQiBnFJKiML6xs4fO5SbvvVf/LEul25f/VbDBw4kGnTpvHrX/86d3uSpCIDcjcgSdo5FtY3cOGC5TRtbGbie/bm1Wef5Gs/eZiUEvfccw8TJkzI3aIkqYhHxCWpTFy2ZAVNG5sBGPnBD/GumsP5w/XnMf3Yj7J582bOOeeczB1Kkop5RFySysSaxqY2jys/Op3Kj04ngH+fe1yepiRJ2+URcUkqEyMqKzpVlyTlZRCXpDIxs7aGioH929QqBvZnZm1Npo4kSTvi0hRJKhNTx1UBLWvFYR1VlRXMrK1prUuSehaDuCSVkanjqpg6roq6ujrOmz45dzuSpB1waYokSZKUgUFckiRJysAgLkmSJGVgEJckSZIyMIhLkiRJGRjEJUmSpAwM4pIkSVIGBnFJkiQpA4O4JEmSlIFBXJIkScrAIC5JkiRlYBCXJEmSMjCIS5IkSRkYxCVJkqQMDOKSJElSBgZxSZIkKQODuCRJkpSBQVySJEnKwCAuSZIkZWAQlyRJkjIwiEuSJEkZGMQlSZKkDAzikiRJUgYGcUmSJCkDg7gkSZKUgUFckiRJyiBrEI+IGyPixYh4vKj2jYhoiIhlhduxRdsujIhnImJFRNQW1cdHxPLCtu9ERHT3e5EkSZI6I/cR8e8Dx7RTvyKlNLZwuxMgIg4ETgVGF57z3YjoX9j/WuAcYP/Crb3XlCRJknqMrEE8pXQf8EoHdz8RmJ9S2pBSehZ4Bjg0IoYDg1NKD6SUEvADYGqXNCxJkiTtJANyN7AdX4iI04FHgK+klF4FqoAHi/ZZXahtLNzfur6NiDiHliPnDBs2jLq6up3feZlbv369cyuBcyuNcyudsyuNcyuNcyudsytNucytJwbxa4F/BlLh6+XADKC9dd9pB/VtiyldB1wHMGHChDR58uSd0G7fUldXh3PrPOdWGudWOmdXGudWGudWOmdXmnKZW+414ttIKb2QUmpOKW0GrgcOLWxaDexTtGs1sKZQr26nLkl91nPPPcfYsWNbb4MHD+bKK6/M3ZYkqUiPOyIeEcNTSmsLD08CtpxRZRHwo4j4NjCClg9lPpRSao6IdRExEfhv4HTg/3Z335LUk+y7774sW7YMgObmZqqqqjjppJPyNiVJaiNrEI+IHwOTgb0iYjXwdWByRIylZXnJSuCzACmlJyLiFuBJYBPw+ZRSc+GlPkfLGVgqgF8WbpIk4J577uEDH/gA73vf+3K3IkkqkjWIp5ROa6d8ww72nwPMaaf+CDBmJ7YmSWVj/vz5nHZae39uJUk59bg14pKkneett95i0aJFfOpTn8rdiiRpKz1ujbgkqXQL6xu4bMkKTt1nHRfNXcqkXVdyyCGHMGzYsNytSZK24hFxSSoTC+sbuHDBchoamwBoaGzimn+bR83/8mLDktQTGcQlqUxctmQFTRubWx9v3vgmr/+hnof4YMauJEnb49IUSSoTawpHwrfoN3A39vnij3nxzUwNSZJ2yCPiklQmRlRWdKouScrLIC5JZWJmbQ0VA/u3qVUM7M/M2ppMHUmSdsSlKZJUJqaOqwJa1orDOqoqK5hZW9NalyT1LAZxSSojU8dVMXVcFXV1dZw3fXLudiRJO+DSFEmSJCkDg7gkSZKUgUFckiRJysAgLkmSJGVgEJckSZIyMIhLkiRJGRjEJUmSpAwM4pIkSVIGBnFJkiQpA4O4JEmSlIFBXJIkScrAIC5JkiRlYBCXJEmSMjCIS5IkSRkYxCVJkqQMDOKSJElSBgZxSZIkKQODuCSVqcbGRj75yU9ywAEHMGrUKB544IHcLUmSigzI3YAkqWt88Ytf5JhjjuGnP/0pb731Fm+88UbuliRJRQziklSGXn/9de677z6+//3vA7DLLruwyy675G1KktSGS1MkqQytXbuWvffem8985jOMGzeOs88+m9dffz13W5KkIgZxSSojC+sbOHzuUlas/TMPP/ooBx55MvX19bz73e9m7ty5uduTJBUxiEtSmVhY38CFC5bT0NhE5ZD30H/QXtz09EAW1jfwyU9+ksceeyx3i5KkIq4Rl6QycdmSFTRtbAZgcOWeDBi8F689/0cuW7ILU968jwMPPDBzh5KkYgZxSSoTaxqb2jwe8vFzefmOb/FS8yb2/l8Hc9NNN2XqTJLUHoO4JJWJEZUVNBSF8V2GvZ/hZ1xJVWUFC2cdmbEzSVJ7XCMuSWViZm0NFQP7t6lVDOzPzNqaTB1JknbEI+KSVCamjqsCWtaKwzqqKiuYWVvTWpck9SwGcUkqI1PHVTF1XBV1dXWcN31y7nYkSTvg0hRJkiQpA4O4JEmSlIFBXJIkScrAIC5JkiRlYBCXJEmSMjCIS5IkSRkYxCVJkqQMDOKSJElSBgZxSZIkKQODuCRJkpSBQVySJEnKwCAuSZIkZWAQlyRJkjIwiEuSJEkZGMQlSZKkDAzikiRJUgYGcUmSJCkDg7gkSZKUgUFckiRJysAgLkmSJGVgEJckSZIyMIhLkiRJGRjEJUmSpAwM4pIkSVIGBnFJkiQpA4O4JEmSlIFBXJIkScrAIC5JkiRlYBCXJEmSMjCIS5IkSRkYxCVJkqQMDOKSJElSBgZxSZIkKQODuCRJkpSBQVySJEnKwCAuSZIkZZA1iEfEjRHxYkQ8XlQbEhG/ioinC1/3LNp2YUQ8ExErIqK2qD4+IpYXtn0nIqK734skSZLUGbmPiH8fOGar2izgnpTS/sA9hcdExIHAqcDownO+GxH9C8+5FjgH2L9w2/o1JUmSpB4laxBPKd0HvLJV+URgXuH+PGBqUX1+SmlDSulZ4Bng0IgYDgxOKT2QUkrAD4qeI0mSJPVIA3I30I5hKaW1ACmltRExtFCvAh4s2m91obaxcH/r+jYi4hxajpwzbNgw6urqdm7nfcD69eudWwmcW2mcW+mcXWmcW2mcW+mcXWnKZW49MYhvT3vrvtMO6tsWU7oOuA5gwoQJafLkyTutub6irq4O59Z5zq00zq10zq40zq00zq10zq405TK33GvE2/NCYbkJha8vFuqrgX2K9qsG1hTq1e3UJUmSpB6rJwbxRcAZhftnALcX1U+NiF0jYj9aPpT5UGEZy7qImFg4W8rpRc+RJEmSeqSsS1Mi4sfAZGCviFgNfB2YC9wSEWcBzwGfAkgpPRERtwBPApuAz6eUmgsv9TlazsBSAfyycJMkSZJ6rKxBPKV02nY2HbWd/ecAc9qpPwKM2YmtSZIkSV2qJy5NkSRJksqeQVySJEnKwCAuSZIkZWAQlyRJkjIwiEuSJEkZGMQlSZKkDAzikiRJUgYGcUmSJCkDg7gkSZKUgUFckiRJysAgLkmSJGVgEJckSZIyMIhLkiRJGRjEJUmSpAwM4pIkSVIGBnFJkiQpA4O4JEmSlIFBXJIkScrAIC5JkiRlYBCXJEmSMjCIS5IkSRkYxCVJkqQMDOKSJElSBgZxSZIkKQODuCRJkpSBQVySJEnKwCAuSZIkZWAQlyRJkjIwiEuSJEkZGMQlSZKkDAzikiRJUgYGcUmSJCkDg7gkSZKUgUFckiRJysAgLkmSJGVgEJckSZIyMIhLkiRJGRjEJUmSpAwM4pIkSVIGBnFJkiQpA4O4JEmSlIFBXJIkScrAIC5JkiRlYBCXJEmSMjCIS5IkSRkM6MhOEbEnMAJoAlamlDZ3aVeSJElSmdtuEI+IPYDPA6cBuwAvAbsBwyLiQeC7KaV7u6VLSZIkqczs6Ij4T4EfAB9NKTUWb4iI8cCnI+L9KaUburA/SZIkqSxtN4inlKbsYNujwKNd0pEkSZLUB3R0jfiHgZHF+6eUFnRRT5IkSVLZe9sgHhE3Ah8GngC2fEgzAQZxSZIkqUQdOSI+MaV0YJd3IkmSJPUhHTmP+AMRYRCXJEmSdqKOHBGfR0sYfx7YAASQUkof7tLOJEmSpDLWkSB+I/BpYDl/WSMuSZIk6R3oSBB/LqW0qMs7kSRJkvqQjgTx/4mIHwE/p2VpCuDpCyVJkqR3oiNBvIKWAH50Uc3TF0qSJEnvwNsG8ZTSZ7qjEUmSJKkv6cgFffYDzmPbK2ue0HVtSZIkSeWtI0tTFgI30LJG3LOmSJIkSTtBR4L4myml73R5J5IkSVIf0pEgflVEfB24i7ZnTXmsy7qSJEmSylxHgvhBtFzQ50j+sjQlFR5LkiRJKkFHgvhJwPtTSm91dTOSJElSX9GvA/v8Bqjs4j4kSZKkPqUjR8SH0XJ1zYdpu0bc0xdKkiRJJepIEP96l3chSZIk9THbDeIREanFf7zdPl3TmiRJklS+drRG/N6IOC8i9i0uRsQuEXFkRMwDzuja9iRJkqTytKOlKccAM4AfFy5z3wjsBvSn5ZziV6SUlnV1g5IkSVI52m4QTym9CXwX+G5EDAT2AppSSo3d1JskSZJUtjryYU1SShuBtV3ciyRJktRndOQ84llExMqIWB4RyyLikUJtSET8KiKeLnzds2j/CyPimYhYERG1+TqXJEmS3l6PDeIFR6SUxqaUJhQezwLuSSntD9xTeExEHAicCoymZW37dyOif46GJUmSpI542yAeEV8oPvKc2YnAvML9ecDUovr8lNKGlNKzwDPAod3fniRJktQxHTki/l7g4Yi4JSKOiYjo6qYKEnBXRDwaEecUasNSSmsBCl+HFupVwKqi564u1CRJkqQeKTpyPZ5C+D4a+AwwAbgFuCGl9PsuayxiREppTUQMBX4FnAcsSilVFu3zakppz4i4BnggpXRzoX4DcGdK6WdbveY5wDkAw4YNGz9//vyuar9srV+/nkGDBuVuo9dxbqVxbqVzdqVxbqVxbqVzdqXpbXM74ogjHi1aat2qo2dNSRHxPPA8sAnYE/hpRPwqpXTBzm219XuuKXx9MSJuo2WpyQsRMTyltDYihgMvFnZfDexT9PRqYE07r3kdcB3AhAkT0uTJk7ui9bJWV1eHc+s851Ya51Y6Z1ca51Ya51Y6Z1eacplbR9aI/z8R8ShwKXA/cFBK6XPAeODkrmgqIt4dEbtvuU/L0fjHgUX85WqeZwC3F+4vAk6NiF0LFx/aH3ioK3qTJEmSdoaOHBHfC5iWUvpjcTGltDkiju+athgG3FZYjj4A+FFKaXFEPAzcEhFnAc8Bnyr08kRE3AI8ScsR+8+nlJq7qDdJkiTpHXvbIJ5S+v92sO2pndtO6+v+ATi4nfqfgKO285w5wJyu6EeSJEna2Xr6ecQlSZKksmQQlyRJkjIwiEtSmZoxYwZDhw5lzJgxuVuRJLXDIC5JZerMM89k8eLFuduQJG2HQVySytSkSZMYMmRI7jYkSdthEJckSZIyMIhLkiRJGRjEJamMLKxv4PC5S1ne8GcOn7uUu554PndLkqTt6MiVNSVJvcDC+gYuXLCcpo3NsA80NDZxyeJVbHxzU+7WJEnt8Ii4JJWJy5asaAnhBS8tupQ/3vRlVj/7DNXV1dxwww0Zu5Mkbc0j4pJUJtY0NrV5vPcJFwAQwLNzj8vQkSRpRzwiLkllYkRlRafqkqS8DOKSVCZm1tZQMbB/m1rFwP7MrK3J1JEkaUdcmiJJZWLquCqgZa04rKOqsoKZtTWtdUlSz2IQl6QyMnVcFVPHVVFXV8d50yfnbkeStAMuTZEkSZIyMIhLkiRJGRjEJUmSpAwM4pIkSVIGBnFJkiQpA4O4JEmSlIFBXJIkScrAIC5JkiRlYBCXJEmSMjCIS5IkSRkYxCVJkqQMDOKSJElSBgZxSZIkKQODuCRJkpSBQVySJEnKwCAuSZIkZWAQlyRJkjIwiEuSJEkZGMQlSZKkDAzikiRJUgYGcUkqU4sXL6ampoYPfvCDzJ07N3c7kqStGMQlqQw1Nzfz+c9/nl/+8pc8+eST/PjHP+bJJ5/M3ZYkqYhBXJLK0P/8z//wwQ9+kPe///3ssssunHrqqdx+++2525IkFTGIS1IZevnll9lnn31aH1dXV9PQ0JCxI0nS1gziklSGUkrb1CIiQyeSpO0xiEtSGVlY38Dhc5fyGhX87L7fsLC+5Sj46tWrGTFiRObuJEnFBuRuQJK0cyysb+DCBctp2tjMvgfuz/oXV/GVG37FxtM/xvz58/nRj36Uu0VJUhGPiEtSmbhsyQqaNjYD0L9/f4ZMOZfnfnQRn/7E4ZxyyimMHj06c4eSpGIeEZekMrGmsanN44oP/BVVH/grArjoouPyNCVJ2i6PiEtSmRhRWdGpuiQpL4O4JJWJmbU1VAzs36ZWMbA/M2trMnUkSdoRl6ZIUpmYOq4KaFkrDuuoqqxgZm1Na12S1LMYxCWpjEwdV8XUcVXU1dVx3vTJuduRJO2AS1MkSZKkDAzikiRJUgYGcUmSJCkDg7gkSZKUgUFckiRJysAgLkmSJGVgEJckSZIyMIhLkiRJGRjEJUmSpAwM4pIkSVIGBnFJkiQpA4O4JEmSlIFBXJIkScrAIC5JkiRlYBCXJEmSMjCIS5IkSRkYxCVJkqQMDOKSJElSBgZxSZIkKQODuCRJkpSBQVySJEnKwCAuSZIkZWAQlyRJkjIwiEuSJEkZGMQlSZKkDAzikiRJUgZlE8Qj4piIWBERz0TErNz9SJIkSTtSFkE8IvoD1wCfAA4ETouIA/N2JUmSJG1fWQRx4FDgmZTSH1JKbwHzgRMz9yRJkiRtV6SUcvfwjkXEJ4FjUkpnFx5/GjgspfSFrfY7BzgHYNiwYePnz5/f7b32duvXr2fQoEG52+h1nFtpnFvpnF1pnFtpnFvpnF1petvcjjjiiEdTShO2rg/I0UwXiHZq2/wXRkrpOuA6gAkTJqTJkyd3cVvlp66uDufWec6tNM6tdM6uNM6tNM6tdM6uNOUyt3JZmrIa2KfocTWwJlMvkiRJ0tsqlyD+MLB/ROwXEbsApwKLMvckSZIkbVdZLE1JKW2KiC8AS4D+wI0ppScytyVJkiRtV1kEcYCU0p3Anbn7kCRJkjqiXJamSJIkSb2KQVySJEnKwCAuSZIkZWAQlyRJkjIwiEtSN2pubmbcuHEcf/zxuVuRJGVmEJekbnTVVVcxatSo3G1IknoAg7gkdZPVq1fzi1/8grPPPjt3K5KkHsAgLknd5Etf+hKXXnop/fr5p1eSZBCXpG5xxx13MHToUMaPH5+7FUlSD1E2V9aUpJ5mYX0Dly1ZwZrGJjY++CPefOpe7rzzTt58801ee+01/v7v/56bb745d5uSpEwM4pLUBRbWN3DhguU0bWwGYMDE6ez90dO5eNpBVP75ab71rW8ZwiWpj3NpiiR1gcuWrGgN4Vs0bWzmsiUrMnUkSeppPCIuSV1gTWPTduuTJx/H5MmTu7chSVKP4xFxSeoCIyorOlWXJPU9BnFJ6gIza2uoGNi/Ta1iYH9m1tZk6kiS1NO4NEWSusDUcVUArWdNGVFZwczamta6JEkGcUnqIlPHVRm8JUnb5dIUSZIkKQODuCRJkpSBQVyS1CusWrWKI444glGjRjF69Giuuuqq3C1J0jviGnFJUq8wYMAALr/8cg455BDWrVvH+PHjmTJlCgceeGDu1iSpJB4RlyT1CsOHD+eQQw4BYPfdd2fUqFE0NDRk7kqSSmcQlyT1OitXrqS+vp7DDjssdyuSVDKXpkiSeqyF9Q3bnIv94/vvwcknn8yVV17J4MGDc7coSSUziEuSeqSF9Q1cuGA5TRubAWhobGLWrfVU3PstPj19OtOmTcvcoSS9My5NkST1SJctWdEawgFSSqxe9G1e6Pcezj///IydSdLOYRCXJPVIaxqb2jze0PAkrz9xLy8//Rhjx45l7Nix3HnnnZm6k6R3zqUpkqQeaURlBQ1FYXy36tG872t3UFVZwf2zjszYmSTtHB4RlyT1SDNra6gY2L9NrWJgf2bW1mTqSJJ2Lo+IS5J6pKnjqgC2OWvKlrok9XYGcUlSjzV1XJXBW1LZcmmKJEmSlIFBXJIkScrAIC5JkiRlYBCXJEmSMjCIS5IkSRkYxCVJkqQMDOKSJElSBgZxSZIkKQODuCRJkpSBQVySpE548803OfTQQzn44IMZPXo0X//613O3JKmX8hL3kiR1wq677srSpUsZNGgQGzdu5CMf+Qif+MQnmDhxYu7WJPUyHhGXJKkTIoJBgwYBsHHjRjZu3EhEZO5KUm9kEJckqZOam5sZO3YsQ4cOZcqUKRx22GG5W5LUC7k0RZKkt7GwvoHLlqxgTWMTIyormFlbw7Jly2hsbOSkk07i8ccfZ8yYMbnblNTLeERckqQdWFjfwIULltPQ2EQCGhqbuHDBchbWN1BZWcnkyZNZvHhx7jYl9UIGcUmSduCyJSto2tjc+rj5jT/z+ro/t9Sbmrj77rs54IADMnYoqbdyaYokSTuwprGpzePm9a/w8i+u4Pm0mb+6+d2ccsopHH/88Zm6k9SbGcQlSdqBEZUVNBSF8V2G7seIz3yHqsoK7p91ZMbOJPV2Lk2R9I6NHDmSgw46iLFjxzJhwoTc7Ug71czaGioG9m9TqxjYn5m1NZk6klQuPCIuaae499572WuvvXK3Ie10U8dVAWxz1pQtdUkqlUFckqS3MXVclcFb0k7n0hRJ71hEcPTRRzN+/Hiuu+663O1IktQreERcUkmKL3Cy12mXcNGnDud/VQ1kypQpHHDAAUyaNCl3i5Ik9WgeEZfUaVtf4OSlze/mwgXL+XXDRk466SQeeuih3C1KktTjGcQldVrxBU42v/Ummze8QdPGZub+fBl33XWXl/qWJKkDXJoiqdOKL3DS/EYjLy34l5b65s18/cvncMwxx+RqTZKkXsMgLqnTii9wMrDyvYyYcTUAVZUVXOQFTiRJ6hCXpkjqNC9wIknSO+cRcUmd5gVOJEl65wzikkriBU4kSXpnXJoiSZIkZWAQl/qQGTNmMHTo0DanF1y2bBkTJ05k7NixTJgwwXOAS5LUTQziUh9y5plnsnjx4ja1Cy64gK9//essW7aMf/qnf+KCCy7I1J0kSX2LQVzqQyZNmsSQIUPa1CKC1157DYA///nPjBgxIkdrkiT1OX5YU+rjrrzySmpra/nqV7/K5s2b+fWvf527JUk7wYwZM7jjjjsYOnQojz/+OAAzZ87k5z//Obvssgsf+MAHuOmmm6isrMzbqNSHeURcKnML6xs4fO5S9pv1Cw6fu5S7nni+zfZrr72WK664glWrVnHFFVdw1llnZepU0s7U3lK0KVOm8Pjjj/Pb3/6WD33oQ1x88cWZupMEBnGprC2sb+DCBctpaGwiAQ2NTVyyeAWvvbmpdZ958+Yxbdo0AD71qU/5YU2pTLS3FO3oo49mwICWfwyfOHEiq1evztGapAKDuFTGLluygqaNzW1qGzY18/L6Da2PR4wYwX/8x38AsHTpUvbff/9u7VFSHjfeeCOf+MQncrch9WmuEZfK2JrGpjaPX1p0KRueW05z02tUV1cze/Zsrr/+er74xS+yadMmdtttN6677rpM3Up6pxbWN7S54u0ZB72r3f3mzJnDgAEDmD59ejd3KKmYQVwqYyMqK2goCuN7n9ByasKqygrun3Vka/3RRx/t9t4k7VxblqJt+VewlqVoq9hYtBQNWpaj3XHHHdxzzz1ERI5WJRW4NEUqYzNra6gY2L9NrWJgf2bW1mTqSFJX6chStMWLF3PJJZewaNEi3vWu9o+WS+o+HhGXytjUcVUAbf6pemZtTWtdUvnoyFK0iy++mA0bNjBlyhSg5QOb3/ve93K0KwmDuFT2po6rMnhLfUBHlqJ5elKpZ+lxS1Mi4hsR0RARywq3Y4u2XRgRz0TEioioLaqPj4jlhW3fCRe9SZL6GJeiSb1PTz0ifkVK6VvFhYg4EDgVGA2MAO6OiA+llJqBa4FzgAeBO4FjgF92b8uSJOXjUjSp9+mpQbw9JwLzU0obgGcj4hng0IhYCQxOKT0AEBE/AKZiEJck9TEuRZN6l54axL8QEacDjwBfSSm9ClTRcsR7i9WF2sbC/a3r24iIc2g5cs6wYcOoq6vb+Z2XufXr1zu3Eji30ji30jm70ji37bvkkkt48MEHqays5KabbgJaLgp0//33k1LiPe95D1/72tfYa6+9Mnfau/gzV5pymVuWIB4RdwPvbWfTRbQsM/lnIBW+Xg7MANpb9512UN+2mNJ1wHUAEyZMSJMnT+5s631eXV0dzq3znFtpnFvpnF1pnNv29evXj0GDBnH66ae3zuiQQw5h8ODB1NXV8dvf/pa7777bs7B0kj9zpSmXuWUJ4imlj3dkv4i4Hrij8HA1sE/R5mpgTaFe3U5dkiTtJJMmTWLlypVtaoMHD269//rrr3uBIKmTetzSlIgYnlJaW3h4EvB44f4i4EcR8W1aPqy5P/BQSqk5ItZFxETgv4HTgf/b3X1LktQXXXTRRVx//fUMHTqUe++9N3c7Uq/S44I4cGlEjKVleclK4LMAKaUnIuIW4ElgE/D5whlTAD4HfB+ooOVDmn5QU5Kkd2BhfcM2Z2AZu+e2+82ZM4cpU6bwwAMPcPXVVzN79uzub1bqpXpcEE8pfXoH2+YAc9qpPwKM6cq+JEnqKxbWN3DhguU0bWw53tXQ2MSFC5bz5YmV233O3/3d33HccccZxKVO6HEX9JEkSXldtmRFawjfomljM/963x/a1J5++unW+4sWLeKAAw7olv6kctHjjohLkqS81jQ2bVN7adGlrH5uObFhHdXV1cyePZs777yTFStW0NTUxIEHHugZU6ROMohLkqQ2RlRW0LBVGN/7hAuoqqzg/llHttbOOussoHxOJSd1N5emSJKkNmbW1lAxsH+bWsXA/sysrcnUkVSePCIuSZLamDqu5QLVW581ZUtd0s5hEJckSduYOq7K4C11MZemSJIkSRkYxKUSzJgxg6FDhzJmzF9OX/+b3/yGv/7rv+aggw7ib/7mb3jttdcydihJkno6g7hUgjPPPJPFixe3qZ199tnMnTuX5cuXc9JJJ3HZZZdl6k6SJPUGBnGpBJMmTWLIkCFtaitWrGDSpEkATJkyhZ/97Gc5WpMkSb2EQVzaScaMGcOiRYsAuPXWW1m1alXmjiRJUk9mEJd2khtvvJFrrrmG8ePHs27dOnbZZZfcLUmSpB7M0xdKHbSwvqHNOXXPOOhdbbYfcMAB3HXXXQD87ne/4xe/+EWONiVJUi/hEXGpAxbWN3DhguU0NDaRgIbGJi5ZvILX3tzUus+LL74IwObNm/mXf/kXzj333EzdSpKk3sAj4lIHXLZkBU0bm1sfv7ToUjY8t5zNTa9RXV3N7NmzWb9+Pddccw0A06ZN4zOf+UyudiVJUi9gEJc6YE1jU5vHe59wAQABPDv3uNb6F7/4xe5sS5Ik9WIuTZE6YERlRafqkiRJb8cgLnXAzNoaKgb2b1OrGNifmbU1mTqSJEm9nUtTpA6YOq4KoM1ZU2bW1rTWJUmSOssgLnXQ1HFVBm9JkrTTuDRFkiRJysAgLkmSJGVgEJckSZIyMIhLkiRJGRjEJUmSpAwM4pIkSVIGBnFJkiQpA4O4JEmSlIFBXJIkScrAIC5JkiRlYBCXJEmSMjCIS5IkSRkYxCVJkqQMDOKSJElSBgZx9QirVq3iiCOOYNSoUYwePZqrrroKgFdeeYUpU6aw//77M2XKFF599dXMnUqSJO0cBnH1CAMGDODyyy/nqaee4sEHH+Saa67hySefZO7cuRx11FE8/fTTHHXUUcydOzd3q5IkSTuFQVw9wvDhwznkkEMA2H333Rk1ahQNDQ3cfvvtnHHGGQCcccYZLFy4MGOXkiRJO49BXD3OypUrqa+v57DDDuOFF15g+PDhQEtYf/HFFzN3J0mStHMYxJXNwvoGDp+7lP1m/YLD5y5lYX0D69ev5+STT+bKK69k8ODBuVuUJKlDZsyYwdChQxkzZsw22771rW8REbz88ssZOlNPZhBXFgvrG7hwwXIaGptIQENjE7NureejU45n+vTpTJs2DYBhw4axdu1aANauXcvQoUMzdi1JUvvOPPNMFi9evE191apV/OpXv2LffffN0JV6OoO4srhsyQqaNja3Pk4psXrRt3mh33s4//zzW+snnHAC8+bNA2DevHmceOKJ3d6rJElvZ9KkSQwZMmSb+pe//GUuvfRSIiJDV+rpBuRuQH3TmsamNo83NDzJ60/cy1t7j2Ts2LEAfPOb32TWrFmccsop3HDDDey7777ceuutGbqVJKnzFi1aRFVVFQcffHDuVtRDGcSVxYjKChqKwvhu1aN539fuoKqygvtnHdlm33vuuae725MkaYcW1jdw2ZIVrGlsYkRlBTNraxi751+2v/HGG8yZM4e77rorX5Pq8Vyaoixm1tZQMbB/m1rFwP7MrK3J1JEkSR3T3uecLlywnLueeL51n9///vc8++yzHHzwwYwcOZLVq1dzyCGH8Pzzz2//hdXnGMS1U2zvypi33noro0ePpl+/fjzyyCOt+08dV8XF0w6iqrKCAKoqK7h42kFMHVeV6R1IktQxW3/OCaBpYzP/et8fWh8fdNBBvPjii6xcuZKVK1dSXV3NY489xnvf+97ublc9mEtTtFNsuTLmIYccwrp16xg/fjxTpkxhzJgxLFiwgM9+9rPbPGfquCqDtySp19n6c04ALy26lNXPLSc2rKO6uprZs2dz1llnZehOvYlBXDvF8OHDWy+8U3xlzClTpmTuTJKknWvrzzkB7H3CBe1+zmmLlStXdkNn6m1cmqKdrvjKmJIklRs/56SdxSPiKll7nxj/+P57eGVMSVJZ27Kscuv/D3S5pTrLIK6SbPnE+JYPq2y5MmbFvd/i00VXxpQkqRz5OSftDC5NUUk6emVMSZIktc8grpJs78qYLz/9GGPHjmXs2LHceeed3HbbbVRXV/PAAw9w3HHHUVtbm6ljSZKknsWlKSpJZ66MedJJJ3V3e5IkST2eR8RVEj8xLkmS9M54RFwl8RPjkiRJ74xBXCXzE+OSJEmlc2mKJEmSlIFBXJIkScrAIN6HrFq1iiOOOIJRo0YxevRorrrqKgBmzpzJAQccwIc//GFOOukkGhsb8zYqSZLUBxjE+5ABAwZw+eWX89RTT/Hggw9yzTXX8OSTTzJlyhQef/xxfvvb3/KhD32Iiy++OHerkiRJZc8g3ocMHz6cQw45BIDdd9+dUaNG0dDQwNFHH82AAS2f2504cSKrV6/O2aYkSVKfYBDvo1auXEl9fT2HHXZYm/qNN97IJz7xiUxdSZIk9R2evrDMLaxv2OZc3x/ffw9OPvlkrrzySgYPHty675w5cxgwYADTp0/P2LEkSVLf4BHxMrawvoELFyynobGJBDQ0NjHr1no+OuV4pk+fzrRp01r3nTdvHnfccQc//OEPiYh8TUuSpC4xY8YMhg4dypgxY1pr3/jGN6iqqmLs2LGMHTuWO++8M2OHfY9BvIxdtmQFTRubWx+nlFi96Nu80O89nH/++a31xYsXc8kll7Bo0SLe9a535WhVkiR1sTPPPJPFixdvU//yl7/MsmXLWLZsGccee2yGzvoug3gZW9PY1ObxhoYnef2Je3n56cfa/JfvF77wBdatW8eUKVMYO3Ys5557bqaOJUlSV5k0aRJDhgzJ3YaKuEa8jI2orKChKIzvVj2a933tDqoqK7h/1pGtdf/rV5Kkvuvqq6/mBz/4ARMmTODyyy9nzz33zN1Sn+ER8TI2s7aGioH929QqBvZnZm1Npo4kSVJP8rnPfY7f//73LFu2jOHDh/OVr3wld0t9ikfEy9jUcVUA25w1ZUtdkiSVt63PnnbGQW0/CzZs2LDW+//wD//A8ccf390t9mkG8TI3dVyVwVuSpD5oy9nTtpy4oaGxiUsWr2Ljm5ta91m7di3Dhw8H4LbbbmtzRhV1PYO4JElSGdr67GkvLbqUDc8tZ3PTa1RXVzN79mzq6upYtmwZEcHIkSP513/914wd9z0GcUmSpDK09dnT9j7hAgACeHbucQCcddZZ3d2WivhhTUmSpDI0orKiU3V1vyxBPCI+FRFPRMTmiJiw1bYLI+KZiFgREbVF9fERsbyw7TtRuPxjROwaET8p1P87IkZ289uRJEnqcTx7Ws+X64j448A04L7iYkQcCJwKjAaOAb4bEVt+gq4FzgH2L9yOKdTPAl5NKX0QuAK4pMu7lyRJ6uGmjqvi4mkHUVVZQQBVlRVcPO0gT+LQg2RZI55SegqgcFC72InA/JTSBuDZiHgGODQiVgKDU0oPFJ73A2Aq8MvCc75ReP5PgasjIlJKqYvfhiRJUo/m2dN6tp62RrwKWFX0eHWhVlW4v3W9zXNSSpuAPwPv6fJO34FVq1ZxxBFHMGrUKEaPHs1VV10FwD/+4z/y4Q9/mLFjx3L00UezZs2azJ1KkiSpq3TZEfGIuBt4bzubLkop3b69p7VTSzuo7+g57fV0Di3LWxg2bBh1dXXbaaNr/elPf+K0007jQx/6EG+88Qaf/exn2WOPPZg4cSJHHXUUAD/72c8499xzOf/887P0uD3r16/PNrfezLmVxrmVztmVxrmVxrmVrhxmd8kll/Dggw9SWVnJTTfd1FpfsGABCxcupF+/fkycOJFzzz13p33PcpgbdGEQTyl9vISnrQb2KXpcDawp1KvbqRc/Z3VEDAD2AF7ZTk/XAdcBTJgwIU2ePLmEFne+Qw45hKqqKqZMmdJae+CBB9i0aRM9pcct6urqelxPvYFzK41zK52zK41zK41zK105zK5fv34MGjSI008/vfW93HvvvTz55JM8/fTT7Lrrrrz44osMHTp0p33Pcpgb9LylKYuAUwtnQtmPlg9lPpRSWgusi4iJhbOlnA7cXvScMwr3Pwks7U3rw1euXEl9fT2HHXYYABdddBH77LMPP/zhD/mnf/qnzN1JkiTt2KRJkxgyZEib2rXXXsusWbPYddddAXZqCC8nuU5feFJErAb+GvhFRCwBSCk9AdwCPAksBj6fUtpySajPAf8GPAP8npYPagLcALyn8MHO84FZ3fZGOmlhfQOHz13KfrN+weFzl/Lj+3/HySefzJVXXsngwYMBmDNnDqtWrWL69OlcffXVmTuWJEnqvN/97nf853/+J4cddhgf+9jHePjhh3O31CPlOmvKbcBt29k2B5jTTv0RYEw79TeBT+3sHne2hfUNXLhgeeulZlf/aR3/cPppfGrq8UybNm2b/f/u7/6O4447jtmzZ3d3q5IkSTu0sL6By5asYE1jEyMqKzjjoHe12b5p0yZeffVVHnzwQR5++GFOOeUU/vCHP7R3xrw+ractTSlbly1Z0RrCU0r86ZdX0W/Pan6398da93n66adb7y9atIgDDjig2/uUJEnakS0HFxsam0hAQ2MTlyxewWtvbmrdp7q6mmnTphERHHroofTr14+XX345X9M9VJYj4n3Rmsam1vsbGp7k9SfuZeDeI3n4irMZO38w3/zmN7nhhhtYsWIF/fr1433vex/f+973MnYsSZK0reKDi1ts2NTMK+s3tD6eOnUqS5cuZfLkyfzud7/jrbfeYq+99uruVns8g3g3GVFZQUMhjO9WPZr3fe0OoOUqV/fPOhKAY489Nlt/kiRJHVF8cBHgpUWXsuG55TQ3vUZ1dTWzZ89mxowZzJgxgzFjxrDLLrswb948l6W0wyDeTWbW1rRZIw5QMbA/M2trMnYlSZLUOcUHFwH2PuECoO3BRYCbb76523vrbVwj3k2mjqvi4mkHUVVZQdDyw3rxtIO87KwkSepVZtbWUDGwf5uaBxdL4xHxbjR1XJXBW5Ik9WpbskzxWVNm1taYcUpgEJckSVKneHBx53BpiiRJkpSBQVySJEnKwCAuSZIkZWAQlyRJkjIwiEuSJEkZGMQlSZKkDAzikiRJUgYGcUmSJCkDg7gkSZKUgUFckiRJysAgLkmSJGVgEJckSZIyMIhLkiRJGRjEJUmSpAwM4pIkSeoTZsyYwdChQxkzZkxr7W//9m8ZO3YsY8eOZeTIkYwdO7bb+hnQbd9JkiRJyujMM8/kC1/4Aqeffnpr7Sc/+Unr/a985Svsscce3daPQVySJEl9wqRJk1i5cmW721JK3HLLLSxdurTb+nFpiiRJkvq8//zP/2TYsGHsv//+3fY9PSIuSZKksrWwvoHLlqxgTWMTIyorOOOgd7W7349//GNOO+20bu3NIC5JkqSytLC+gQsXLKdpYzMADY1NXLJ4FRvf3NRmv02bNrFgwQIeffTRbu3PIC5JkqSydNmSFa0hfIsNm5p5Zf2GNrW7776bAw44gOrq6u5szzXikiRJKk9rGpvaPH5p0aU8/+9fpemlVVRXV3PDDTcAMH/+/G5flgIeEZckSVKZGlFZQUNRGN/7hAsAqKqs4P5ZR7bWv//973d3a4BHxCVJklSmZtbWUDGwf5taxcD+zKytydRRWx4RlyRJUlmaOq4KoM1ZU2bW1rTWczOIS5IkqWxNHVfVY4L31lyaIkmSJGVgEJckSZIyMIhLkiRJGRjEJUmSpAwM4pIkSVIGBnFJkiQpA4O4JEmSlIFBXJIkScrAIC5JkiRlYBCXJEmSMjCIS5IkSRkYxCVJkqQMDOKSJElSBgZxSZIkKQODuCRJkpSBQVySJEnKwCAuSZIkZWAQlyRJkjIwiEuSJEkZGMQlSZKkDAzikiRJUgYGcUmSJCkDg7gkSZKUgUFckiRJysAgLkmSJGUQKaXcPWQRES8Bf8zdRy+0F/By7iZ6IedWGudWOmdXGudWGudWOmdXmt42t/ellPbeuthng7hKExGPpJQm5O6jt3FupXFupXN2pXFupXFupXN2pSmXubk0RZIkScrAIC5JkiRlYBBXZ12Xu4FeyrmVxrmVztmVxrmVxrmVztmVpizm5hpxSZIkKQOPiEuSJEkZGMQlSZKkDAziahURn4qIJyJic0RMKKqPjIimiFhWuH2vaNv4iFgeEc9ExHciIgr1XSPiJ4X6f0fEyAxvqVtsb26FbRcWZrAiImqL6n1+bluLiG9EREPRz9mxRds6Nce+LCKOKczpmYiYlbufniYiVhZ+ZpZFxCOF2pCI+FVEPF34umfR/u3+7PUFEXFjRLwYEY8X1To9q772e7qdufn37W1ExD4RcW9EPFX4/9QvFurl/TOXUvLmjZQSwCigBqgDJhTVRwKPb+c5DwF/DQTwS+AThfr/Br5XuH8q8JPc7y/D3A4EfgPsCuwH/B7o79y2O8dvAF9tp97pOfbVG9C/MJ/3A7sU5nZg7r560g1YCey1Ve1SYFbh/izgksL97f7s9YUbMAk4pPjvfymz6mu/p9uZm3/f3n5uw4FDCvd3B35XmE9Z/8x5RFytUkpPpZRWdHT/iBgODE4pPZBafvJ/AEwtbD4RmFe4/1PgqB79X6TvwA7mdiIwP6W0IaX0LPAMcKhz67RS5thXHQo8k1L6Q0rpLWA+LfPTjhX/3s2j7e/jNj973d9eHiml+4BXtip3alZ98fd0O3PbHudWkFJam1J6rHB/HfAUUEWZ/8wZxNVR+0VEfUT8R0R8tFCrAlYX7bO6UNuybRVASmkT8GfgPd3VbA/ROoOCLfNxbtv3hYj4beGfdrf882Mpc+yrtjcr/UUC7oqIRyPinEJtWEppLbSEAWBooe48t9XZWfl7+hf+feugwrLMccB/U+Y/cwNyN6DuFRF3A+9tZ9NFKaXbt/O0tcC+KaU/RcR4YGFEjKbln3y2tuV8mDva1uuUOLftzaDPzG1rO5ojcC3wz7S8338GLgdmUNoc+ypn8vYOTymtiYihwK8i4n92sK/z7Dh/T3fMv28dFBGDgJ8BX0opvbaDfxQui9kZxPuYlNLHS3jOBmBD4f6jEfF74EO0/FdmddGu1cCawv3VwD7A6ogYAOxBx/+prscpZW78ZQZbbJlPn5nb1jo6x4i4Hrij8LCUOfZV25uVClJKawpfX4yI22hZavJCRAxPKa0t/LP2i4Xdnee2Ojsrf0+BlNILW+779237ImIgLSH8hymlBYVyWf/MuTRFbysi9o6I/oX77wf2B/5Q+CeidRExsbCO+XRgy9HhRcAZhfufBJYW1mr1JYuAU6PlTCj70TK3h5xb+wp/YLc4CdhyxoFS5thXPQzsHxH7RcQutHzgd1HmnnqMiHh3ROy+5T5wNC0/Z8W/d2fQ9vdxm5+97u26x+nUrPw9beHft7dXeJ83AE+llL5dtKm8f+Zyf1rUW8+50fLHYTUtR79fAJYU6icDT9Dy6eTHgL8pes4EWv6g/B64mr9crXU34FZaPjzxEPD+3O+vu+dW2HZRYTYrKPrUtnNrd47/DiwHfkvLH9jhpc6xL9+AY2k528DvaVk6lb2nnnKj5WwyvyncntgyH1o+h3EP8HTh65Ci57T7s9cXbsCPaVmauLHwN+6sUmbV135PtzM3/769/dw+QssSkt8Cywq3Y8v9Z85L3EuSJEkZuDRFkiRJysAgLkmSJGVgEJckSZIyMIhLkiRJGRjEJUmSpAwM4pKkVhFRERH/seXaATv5te8uurS3JPV5BnFJUrEZwIKUUnMXvPa/A/+7C15Xknolg7gk9QER8VcR8duI2K1whcknImJMO7tOp3AVuoiYXDg6fktE/C4i5kbE9Ih4KCKWR8QHCvt9PyKujYh7I+IPEfGxiLgxIp6KiO8XvfYi4LQuf7OS1EsMyN2AJKnrpZQejohFwL8AFcDNKaXHi/eJiF1ouZrryqLywcAo4BXgD8C/pZQOjYgvAucBXyrstydwJHAC8HPgcOBs4OGIGJtSWpZSerVwOer3pJT+1FXvVZJ6C4+IS1Lf8U/AFFou/3xpO9v3Ahq3qj2cUlqbUtpAy+Wi7yrUlwMji/b7eWq5VPNy4IWU0vKU0mZaLiVfvN+LwIh39jYkqTwYxCWp7xgCDAJ2B3ZrZ3tTO/UNRfc3Fz3eTNt/Vd3Qzj7t7bdb4ftIUp9nEJekvuM64B+BHwKXbL0xpfQq0D8i2gvp71hEBPBeYGVXvL4k9TYGcUnqAyLidGBTSulHwFzgryLiyHZ2vQv4SBe1MR54MKW0qYteX5J6lWhZ0idJEkTEOOD8lNKnu+C1rwIWpZTu2dmvLUm9kUfEJUmtUkr1wL1dcUEf4HFDuCT9hUfEJUmSpAw8Ii5JkiRlYBCXJEmSMjCIS5IkSRkYxCVJkqQMDOKSJElSBv8/YIL4kj6N8CMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x864 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def i_y(i):\n",
    "    return 2*i+1\n",
    "\n",
    "def i_x(i):\n",
    "    return 2*i\n",
    "\n",
    "\n",
    "initial_guess = np.zeros(2*n_ant)\n",
    "\n",
    "for i in range(n_ant):\n",
    "    x,y = current_pos[i][0:2]*1000\n",
    "    initial_guess[i_x(i)] = x\n",
    "    initial_guess[i_y(i)] = y\n",
    "    \n",
    "    print(x,y)\n",
    "initial_guess\n",
    "\n",
    "pos_i = current_pos*1000\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(12,12))\n",
    "\n",
    "ax.scatter(pos_i[:,0], pos_i[:,1])\n",
    "for i in range(0,pos_i.shape[0]):\n",
    "    ax.text(pos_i[i,0], pos_i[i,1],f\"{i}\")\n",
    "\n",
    "ax.grid(True)\n",
    "ax.set_title('Initial positions')\n",
    "ax.set_xlabel('x (mm)')\n",
    "ax.set_ylabel('y (mm)')\n",
    "plt.savefig('initial_positions.png')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Criteria for Optimality\n",
    "\n",
    "The function below is minimized when the positions (in variable x) are consistent with the measured distances m[i,j]. The initial value of this function is more than 3 million.\n",
    "\n",
    "Note that the x input is a 1D vector of with 48 entries as [p0.x, p0.y, p1.x, p1.y]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "662006.2323637417\n"
     ]
    }
   ],
   "source": [
    "def dist(a,b):\n",
    "    return np.sqrt((a[0]-b[0])**2 + (a[1]-b[1])**2)\n",
    "\n",
    "\n",
    "def p(x, i):\n",
    "    return [x[i_x(i)],x[i_y(i)]]\n",
    "\n",
    "max_err = 1200\n",
    "def bounds():\n",
    "    ret = [(0,0)]*48\n",
    "    for i in range(24):\n",
    "        _x,_y = pos_i[i,0], pos_i[i,1]\n",
    "        ret[i_x(i)] = (_x - max_err, _x + max_err)\n",
    "        ret[i_y(i)] = (_y - max_err, _y + max_err)\n",
    "    # x-component of the north arm set to zero.\n",
    "    for i in [0,1,6,7,8,9,10,11]:\n",
    "        xval = pos_i[i,0]\n",
    "        ret[i_x(i)] = ( - 2, + 2)\n",
    "    \n",
    "    return ret\n",
    "\n",
    "bnd = bounds()\n",
    "\n",
    "center = [0,0]\n",
    "\n",
    "def radius_model(x):\n",
    "    ret = []\n",
    "    \n",
    "    for i in range(n_ant):\n",
    "        p_i = p(x,i)\n",
    "        ret.append(dist(center, p_i))\n",
    "        \n",
    "    return np.array(ret)\n",
    "\n",
    "non_nan_ij = []\n",
    "non_nan_indices = []\n",
    "for i in range(n_ant):\n",
    "    for j in range(n_ant):\n",
    "        if not np.isnan(m_ij[i,j]):\n",
    "            non_nan_ij.append(m_ij[i,j])\n",
    "            non_nan_indices.append([i,j])\n",
    "non_nan_ij = np.array(non_nan_ij)\n",
    "non_nan_indices = np.array(non_nan_indices)\n",
    "\n",
    "def m_ij_model(x):\n",
    "    ret = []\n",
    "    for i in range(n_ant):\n",
    "        p_i = p(x,i)\n",
    "\n",
    "        for j in range(n_ant):\n",
    "            p_j = p(x,j)\n",
    "            if not np.isnan(m_ij[i,j]):\n",
    "                ret.append(dist(p_j, p_i))\n",
    "                \n",
    "    return np.array(ret)\n",
    "    \n",
    "def radius_residual(x):\n",
    "    predict = radius_model(x)\n",
    "    return predict - radius\n",
    "\n",
    "def m_ij_residual(x):\n",
    "    predict = m_ij_model(x)\n",
    "    return predict - non_nan_ij\n",
    "\n",
    "\n",
    "def f(x):\n",
    "    ret = 0\n",
    "    for i in range(n_ant):\n",
    "        p_i = p(x,i)\n",
    "        ret += (dist(center, p_i) - m_0[i])**2\n",
    "        for j in range(n_ant):\n",
    "            p_j = p(x,j)\n",
    "            if not np.isnan(m_ij[i,j]):\n",
    "                ret += (dist(p_j, p_i) - m_ij[i,j])**2\n",
    "    return ret\n",
    "\n",
    "def f(x):\n",
    "    # This is the final criterion for optimality. The sum of squared residuals.\n",
    "    ret = np.sum(radius_residual(x)**2)\n",
    "    ret += np.sum(m_ij_residual(x)**2)\n",
    "    return ret\n",
    "\n",
    "print(f(initial_guess))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the initial value of this optimality condition is about 600,000."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "      fun: 18555.484476675192\n",
       " hess_inv: <48x48 LbfgsInvHessProduct with dtype=float64>\n",
       "      jac: array([ 3.13426428e+01,  2.54658314e-03,  3.24180293e+00,  3.27417833e-03,\n",
       "       -3.63797592e-04, -1.45519037e-03,  2.18278555e-03,  1.09139278e-03,\n",
       "       -3.63797592e-04,  4.36557731e-03, -5.09316629e-03, -5.09316629e-03,\n",
       "        1.45519153e-02, -1.09140518e-03, -4.72937248e-03, -7.27603456e-04,\n",
       "        1.01863407e-02, -2.54661210e-03,  7.27595758e-04,  3.63801728e-04,\n",
       "        4.36557457e-03, -8.36743975e-03, -3.31710910e+00,  3.63801728e-03,\n",
       "       -3.63797592e-03, -7.27595184e-04,  5.09316629e-03, -1.81898796e-03,\n",
       "       -1.09140518e-03, -1.81898796e-03, -2.54661210e-03, -8.36734462e-03,\n",
       "        1.45520691e-03,  6.18455907e-03, -7.63983629e-03, -9.45884493e-03,\n",
       "        0.00000000e+00,  3.63797592e-04,  1.45520691e-03,  2.18278555e-03,\n",
       "        7.27603456e-04, -2.54658314e-03, -2.91041382e-03,  1.81898796e-03,\n",
       "        7.27603456e-04, -1.01863326e-02, -4.00181901e-03,  6.54843111e-03])\n",
       "  message: 'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'\n",
       "     nfev: 3234\n",
       "      nit: 62\n",
       "     njev: 66\n",
       "   status: 0\n",
       "  success: True\n",
       "        x: array([-2.00000000e+00,  2.07931748e+02, -2.00000000e+00,  4.06403029e+02,\n",
       "        3.45776649e+02, -1.64594379e+02,  5.79393553e+02, -2.72438636e+02,\n",
       "       -1.75629176e+02, -1.12436427e+02, -6.27158479e+02, -3.99822685e+02,\n",
       "       -1.49862541e+00,  1.26055081e+03,  1.35879335e+00,  1.43325171e+03,\n",
       "        1.67338546e+00,  1.69723404e+03,  5.34856329e-01,  1.86233906e+03,\n",
       "        2.19682211e-01,  2.19173586e+03,  2.00000000e+00,  2.38774332e+03,\n",
       "        8.50544014e+02, -3.98618229e+02,  9.73091742e+02, -4.68925948e+02,\n",
       "        1.46282967e+03, -6.83921297e+02,  1.73978604e+03, -8.16475497e+02,\n",
       "        1.88796177e+03, -8.86273483e+02,  2.17693451e+03, -1.02533613e+03,\n",
       "       -7.62698663e+02, -4.86369341e+02, -1.05686361e+03, -6.68820930e+02,\n",
       "       -1.33158606e+03, -8.39309842e+02, -1.43828022e+03, -9.03764619e+02,\n",
       "       -1.54836211e+03, -1.00336532e+03, -1.75585941e+03, -1.10625269e+03])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res = minimize(f, initial_guess, bounds=bnd)\n",
    "res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ant 0:\tJac: 031.34, \t000.00\n",
      "Ant 1:\tJac: 003.24, \t000.00\n",
      "Ant 2:\tJac: -00.00, \t-00.00\n",
      "Ant 3:\tJac: 000.00, \t000.00\n",
      "Ant 4:\tJac: -00.00, \t000.00\n",
      "Ant 5:\tJac: -00.01, \t-00.01\n",
      "Ant 6:\tJac: 000.01, \t-00.00\n",
      "Ant 7:\tJac: -00.00, \t-00.00\n",
      "Ant 8:\tJac: 000.01, \t-00.00\n",
      "Ant 9:\tJac: 000.00, \t000.00\n",
      "Ant 10:\tJac: 000.00, \t-00.01\n",
      "Ant 11:\tJac: -03.32, \t000.00\n",
      "Ant 12:\tJac: -00.00, \t-00.00\n",
      "Ant 13:\tJac: 000.01, \t-00.00\n",
      "Ant 14:\tJac: -00.00, \t-00.00\n",
      "Ant 15:\tJac: -00.00, \t-00.01\n",
      "Ant 16:\tJac: 000.00, \t000.01\n",
      "Ant 17:\tJac: -00.01, \t-00.01\n",
      "Ant 18:\tJac: 000.00, \t000.00\n",
      "Ant 19:\tJac: 000.00, \t000.00\n",
      "Ant 20:\tJac: 000.00, \t-00.00\n",
      "Ant 21:\tJac: -00.00, \t000.00\n",
      "Ant 22:\tJac: 000.00, \t-00.01\n",
      "Ant 23:\tJac: -00.00, \t000.01\n"
     ]
    }
   ],
   "source": [
    "err = res.jac.reshape((24,2))\n",
    "rms = err # np.sqrt(np.sum(err**2, axis=1))\n",
    "for i in range(24):\n",
    "    print(f\"Ant {i}:\\tJac: {rms[i,0]:06.2f}, \\t{rms[i,1]:06.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The optimized positions are now known. The final value of the function is 32. Far closer to zero than 3 million!\n",
    "\n",
    "We can recover the x,y coordinates by reshaping the array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-2.00000000e+00,  2.07931748e+02],\n",
       "       [-2.00000000e+00,  4.06403029e+02],\n",
       "       [ 3.45776649e+02, -1.64594379e+02],\n",
       "       [ 5.79393553e+02, -2.72438636e+02],\n",
       "       [-1.75629176e+02, -1.12436427e+02],\n",
       "       [-6.27158479e+02, -3.99822685e+02],\n",
       "       [-1.49862541e+00,  1.26055081e+03],\n",
       "       [ 1.35879335e+00,  1.43325171e+03],\n",
       "       [ 1.67338546e+00,  1.69723404e+03],\n",
       "       [ 5.34856329e-01,  1.86233906e+03],\n",
       "       [ 2.19682211e-01,  2.19173586e+03],\n",
       "       [ 2.00000000e+00,  2.38774332e+03],\n",
       "       [ 8.50544014e+02, -3.98618229e+02],\n",
       "       [ 9.73091742e+02, -4.68925948e+02],\n",
       "       [ 1.46282967e+03, -6.83921297e+02],\n",
       "       [ 1.73978604e+03, -8.16475497e+02],\n",
       "       [ 1.88796177e+03, -8.86273483e+02],\n",
       "       [ 2.17693451e+03, -1.02533613e+03],\n",
       "       [-7.62698663e+02, -4.86369341e+02],\n",
       "       [-1.05686361e+03, -6.68820930e+02],\n",
       "       [-1.33158606e+03, -8.39309842e+02],\n",
       "       [-1.43828022e+03, -9.03764619e+02],\n",
       "       [-1.54836211e+03, -1.00336532e+03],\n",
       "       [-1.75585941e+03, -1.10625269e+03]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pos = res.x.reshape((24,2))\n",
    "pos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuIAAALJCAYAAADxiEbKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABAVUlEQVR4nO39f7iddXkn+r9vwoZsCyYoGsyPGXB0sAgZkejxW9pzQtGJtQNGz5TieKqe1oHpaOmcb7+ZGmsppWPF5ky5ynG0Muqljkxp7GCMrZ6MRWJLh1Z+hCYizZQqliSgYJtIZjYQ4uf7x1qBnbAT9kr23s/+8Xpd17rW2vf6se91s/fOm2d9nuep1loAAICpdVzXDQAAwFwkiAMAQAcEcQAA6IAgDgAAHRDEAQCgA4I4AAB0QBAHmCGqam9VvXgCXueqqvrMRPQ0ka9bVV+qqrdPZE8A09nxXTcAwMGq6v4ki5LsH1X+x621k6bo+783yb9M8oIku5P8WWvtpyf4e1yV5CWttf/jQK219hMT+T0ApjtbxAGmp4taayeNuuyaim/a3yL9M0le2w/+K5LcPBXfG2CuEcQBZoiqalX1kv7tT1bVf6iqP6qqR6vqL6rqH4167O9U1QNV9f2qurOqfmyc3+ZVSTa11v4mSVprD7XWrh/1uouramNV/V1V3VdV//Iwva6sqh2H1O6vqtdW1euTvDfJT/eX2/xl//7NVfXO/u3jqup9VfXtqvpuVX26qhb07zu9P4u3V9XfVtUjVfUr450jwHQhiAPMXG9J8utJTklyX5L3j7rv9iSvSPK8JP85yWerav44XvPPk7ytqtZU1YqqmnfI/b+XZEeSxUn+eZLfrKoLB2m6tfb/JvnNJL/f39r/T8Z42Dv6lwuSvDjJSUk+dMhjfjTJmUkuTHJlVf3wIH0AdE0QB5ieNlTV7v5lw2Eec1Nr7WuttSeT3JBe8E6StNY+01r7Xmvtydbav09yYnqh9Yhaa59J8gtJViX5apLvVtV7kqSqlqUXfn+5tfZYa+3uJB9LbynLRHtrkt9urX2ztbY3ydokl1bV6H2bfr21NtJa+8skf5lkrEAPMG3ZWRNgelrdWvvjZ3nMQ6Nu/8/0thonSarql5K8M70t1y3Jc5OcOp5v3Fq7IckNVTWUZHX/9pb0dtz8u9bao6Me/u301pFPtMX91x79fY5PbyfWAw77/gFmAlvEAWaZ/nrwX05ySZJTWmsLk+xJUoO8TmttX2vts0m2Jjk7ya4kz6uqk0c97B8k2TnG0/9HkueM6mleekdheerln+Xb70ryDw/5Pk8m+c643wDANCeIA8w+J6cXWh9OcnxVXZneFvFnVVXvqKqfrKqT+ztM/kSSlyf5i9baA0n+W5IPVNX8qlqe5OfSWxZzqP+eZH7/tYaSvC+95TEHfCfJ6VV1uH+Hfi/J/1VVZ1TVSXl6TfmT43kfADOBIA4w+2xK8qX0wvC3kzyW5IFxPvf76R3R5G/TW4ryW0l+vrV2a//+tyQ5Pb0t1p9L8muttS8f+iKttT1J/nV6a8h3preFfPRRVD7bv/5eVd01Rh+fSPKfkvxJkm/138MvjPM9AMwI1dqzfToIAABMNFvEAQCgA4I4AAB0QBAHAIAOCOIAANCBOXtCn1NPPbWdfvrpXbcxY/yP//E/8kM/9ENdtzGjmNlgzGtwZjY4MxuMeQ3OzAYzV+Z15513PtJae8Gh9TkbxE8//fTccccdXbcxY2zevDkrV67suo0ZxcwGY16DM7PBmdlgzGtwZjaYuTKvqvr2WHVLUwAAoAOCOAAAdEAQBwCADszZNeIAAEy8ffv2ZceOHXnsscee9bELFizIvffeOwVdTY358+dn6dKlGRoaGtfjBXEAACbMjh07cvLJJ+f0009PVR3xsY8++mhOPvnkKepscrXW8r3vfS87duzIGWecMa7nWJoCAMCEeeyxx/L85z//WUP4bFNVef7znz+uTwIO6CyIV9Wyqrqlqu6tqnuq6hf79auqamdV3d2/vGHUc9ZW1X1Vtb2qVo2qn1dV2/r3XVdz7b88AMA0Mlej2KDvu8ulKU8m+aXW2l1VdXKSO6vqy/37rm2t/d+jH1xVZyW5NMnLkyxO8sdV9Y9ba/uTfCTJZUn+PMkXk7w+yZem6H0AAMDAOtsi3lp7sLV2V//2o0nuTbLkCE95Y5IbW2uPt9a+leS+JK+uqhcleW5r7bbWWkvy6SSrJ7d7AACmqx/5kR951se8853vzDe+8Y0kyW/+5m8O/PyTTjrp6JobZVqsEa+q05Ocm+Qv+qV3V9XWqvpEVZ3Sry1J8sCop+3o15b0bx9aBwBgDvpv/+2/PetjPvaxj+Wss85K8swgPp7nT4TOj5pSVScl+S9J/k1r7ftV9ZEkv5Gk9a//fZKfTTLWopt2hPpY3+uy9JawZNGiRdm8efMx9z9X7N2717wGZGaDMa/BmdngzGww5jU4M+sdkvDRRx8d12P379+fG2+7L79zy/156PuP57TnnphfvOD0/OTZi46phxe96EV58MEH86d/+qf5wAc+kOc///n5xje+kVe84hX52Mc+lqrKG97whvy7f/fv8vnPfz4jIyNZvnx5Xvayl+XjH//4U8/fu3dv3vKWt2T37t3Zt29ffvVXfzU/+ZM/+dT3Get9PvbYY+P+Geg0iFfVUHoh/IbW2k1J0lr7zqj7/2OSP+x/uSPJslFPX5pkV7++dIz6M7TWrk9yfZKsWLGirVy5ckLex1ywefPmmNdgzGww5jU4MxucmQ3GvAZnZsm999477kMS3njbffn1L96XkX37kyQPfv/x/PoX78v8+cNZfe6xLXA4+eST85znPCdbt27NPffck8WLF+f888/P1q1b86M/+qOZN29efuiHfii//du/neuvvz5bt259xvOHh4ezcePGPPe5z80jjzyS17zmNfnpn/7pp3bKHOt9zp8/P+eee+64euzyqCmV5ONJ7m2t/fao+otGPexNSb7ev70xyaVVdWJVnZHkpUm+1lp7MMmjVfWa/mu+Lcnnp+RNAABw1H7nlvufCuEHjOzbn3Wbtk/Y93j1q1+dpUuX5rjjjssrXvGK3H///eN+bmst733ve7N8+fK89rWvzc6dO/Od73zn2Z84Tl1uET8/yc8k2VZVd/dr703ylqp6RXrLS+5PcnmStNbuqar1Sb6R3hFX3tU/YkqS/HySTyYZTu9oKY6YAgAwzT30/cfHrO/aPTJh3+PEE0986va8efPy5JNPjvu5N9xwQx5++OHceeedGRoayumnnz7QccKfTWdBvLV2a8Ze3/3FIzzn/UneP0b9jiRnT1x3AABMttOee2IeHCOML144PKV9DA0NZd++fc84Nf2ePXvywhe+MENDQ7nlllvy7W9/e0K/77Q4agoAAHPPL15weoaH5h1UGx6alzWrzpzSPi677LIsX748b33rWw+qv/Wtb80dd9yRFStW5IYbbsjLXvayCf2+nR81BQCAueknz16U+fOHs27T9uzaPZLFC4ezZtWZx7yj5t69e5MkK1euPGjn2Q996ENP3R59ZJMPfvCD+eAHP/iM55966qm57bbbjvg9joUgDgBAZ1afu+SYg/dMZWkKAAB0QBAHAIAOCOIAANABQRwAADogiAMAQAcEcQAAZpXrrrsuP/zDP5xTTjkl11xzzVG/zkknnTSBXT2TwxcCADCrfPjDH86XvvSlnHHGGV23ckS2iAPMMLdv/Ggeuuolya6789BVL8ntGz/adUsAR2/r+uTas5OrFvaut64/ppf7V//qX+Wb3/xmLr744lx77bV597vfnSR5xzvekSuuuCI/8iM/khe/+MX5gz/4gyS9E/NceOGFeeUrX5lzzjknn//854/1HY2bIA4wg9y+8aM5+8735bQ8nFRyWh7O2Xe+TxgHZqTj7/1c8oUrkj0PJGm96y9ccUxh/Hd/93ezePHi3HLLLTnllFMOuu/BBx/Mrbfemj/8wz/Me97zniTJ/Pnz87nPfS533XVXbrnllvzSL/1SWmvH8rbGTRAHmEGW3bUuw/XEQbXheiLL7lrXUUcAR+/EP70m2TdycHHfSHLz1ZPy/VavXp3jjjsuZ511Vr7zne8kSVpree9735vly5fnta99bXbu3PnUfZPNGnGAGeSFrbcl/Jn1R6a+GYBjVI/uGvuOPTsm5fudeOKJT90+sNX7hhtuyMMPP5w777wzQ0NDOf300/PYY49Nyvc/lC3iADPId+sFh6mfOsWdABy7dvLise9YsHTKetizZ09e+MIXZmhoKLfccku+/e1vT9n3FsQBZpAHXrkmI+2Eg2oj7YQ88Mo1HXUEcPQe/7H3JEPDBxeHhpMLr5yyHt761rfmjjvuyIoVK3LDDTfkZS972ZR9b0tTAGaQV118eW5Pb614WvJQXpAHzluTV118edetAQzsyR9+UzJ/fm9N+J4dvS3hF16ZLL/kmF73/vvvT9I7Uso73vGOJMknP/nJgx6zd+/eJMmpp56a2267bczXOfCYySKIA8wwr7r48uTiy/NXmzfntH9xX07ruiGAY7H8kmMO3jOVpSkAANABQRwAgAk1Vcfhnm4Gfd+COAAAE2b+/Pn53ve+N+fCeGst3/ve9zJ//vxxP8cacQAAJszSpUuzY8eOPPzww8/62Mcee2yg4DrdzZ8/P0uXjv/Qi4I4AAATZmhoKGeccca4Hrt58+ace+65k9zR9GVpCgAAdEAQBwCADgjiAADQAUEcAAA6IIgDAEAHBHEAAOiAIA4AAB0QxAEAoAOCOAAAdEAQBwCADgjiAADQAUEcAAA6IIgDAEAHBHEAAOiAIA4AAB0QxAEAoAOCOAAAdEAQBwCADgjiADPMhi07c/41X8m2nXty/jVfyYYtO7tuCYCjcHzXDQAwfhu27Mzam7ZlZN/+ZFmyc/dI1t60LUmy+twlHXcHwCBsEQeYQdZt2p7X7f9qbj3hipxT38qtJ1yR1+3/atZt2t51awAMSBAHmEFWfP/LuWboY1l63CNJJUuPeyTXDH0sK77/5a5bA2BAgjjADLL2hM/mOfXEQbXn1BNZe8JnO+oIgKMliAPMIIvyyEB1AKYvQRxgBqkFSweqAzB9CeIAM8mFVyZDwwfXhoZ7dQBmFIcvBJhJll/Su7756t71gmW9EH6gDsCMIYgDzDTLL+ldNm9O3vL1rrsB4ChZmgIAAB0QxAEAoAOCOAAAdEAQBwCADgjiAADQAUEcAAA6IIgDAEAHBHEAAOiAIA4AAB0QxAEAoAOCOAAAdEAQBwCADgjiAADQAUEcAAA6IIgDAEAHBHEAAOiAIA4AAB0QxAEAoAOCOAAAdEAQBwCADgjiAADQAUEcAAA6IIgDAEAHOgviVbWsqm6pqnur6p6q+sV+/XlV9eWq+uv+9SmjnrO2qu6rqu1VtWpU/byq2ta/77qqqi7eEwAAjFeXW8SfTPJLrbUfTvKaJO+qqrOSvCfJza21lya5uf91+vddmuTlSV6f5MNVNa//Wh9JclmSl/Yvr5/KNwIAAIPqLIi31h5srd3Vv/1oknuTLEnyxiSf6j/sU0lW92+/McmNrbXHW2vfSnJfkldX1YuSPLe1dltrrSX59KjnAADAtFS97NpxE1WnJ/mTJGcn+dvW2sJR9/19a+2UqvpQkj9vrX2mX/94ki8luT/JNa211/brP5bkl1tr/2yM73NZelvOs2jRovNuvPHGyXxbs8revXtz0kkndd3GjGJmgzGvwZnZ4MxsMOY1ODMbzFyZ1wUXXHBna23FofXju2hmtKo6Kcl/SfJvWmvfP8Ly7rHuaEeoP7PY2vVJrk+SFStWtJUrVw7c71y1efPmmNdgzGww5jU4MxucmQ3GvAZnZoOZ6/Pq9KgpVTWUXgi/obV2U7/8nf5yk/Svv9uv70iybNTTlybZ1a8vHaMOAADTVpdHTakkH09yb2vtt0fdtTHJ2/u3357k86Pql1bViVV1Rno7ZX6ttfZgkker6jX913zbqOcAAMC01OXSlPOT/EySbVV1d7/23iTXJFlfVT+X5G+T/FSStNbuqar1Sb6R3hFX3tVa299/3s8n+WSS4fTWjX9pit4DAAAclc6CeGvt1oy9vjtJLjzMc96f5P1j1O9Ib0dPAACYEZxZEwAAOiCIAwBABwRxAADogCAOAAAdEMQBAKADgjgAAHRAEAcAgA4I4gAA0AFBHAAAOiCIAwBABwRxAADogCAOAAAdEMQBAKADgjgAAHRAEAeYYW7f+NE8dNVLkl1356GrXpLbN36065YAOAqCOMAMcvvGj+bsO9+X0/JwUslpeThn3/k+YRxgBhLEAWaQZXety3A9cVBtuJ7IsrvWddQRAEdLEAeYQV7YHj5M/ZEp7gSAYyWIA8wg360XHKZ+6hR3AsCxEsQBZpAHXrkmI+2Eg2oj7YQ88Mo1HXUEwNE6vusGABi/V118eW5Pb614WvJQXpAHzluTV118edetATAgQRxghnnVxZcnF1+ev9q8Oaf9i/tyWtcNAXBULE0BAIAOCOIAANABQRwAADogiAMAQAcEcQAA6IAgDgAAHRDEAQCgA4I4AAB0QBAHAIAOCOIAANABQRwAADogiAMAQAcEcQAA6IAgDgAAHRDEAQCgA4I4AAB0QBAHAIAOCOIAM8yGLTtz/jVfybade3L+NV/Jhi07u24JgKNwfNcNADB+G7bszNqbtmVk3/5kWbJz90jW3rQtSbL63CUddwfAIGwRB5hB1m3antft/2puPeGKnFPfyq0nXJHX7f9q1m3a3nVrAAxIEAeYQVZ8/8u5ZuhjWXrcI0klS497JNcMfSwrvv/lrlsDYECCOMAMsvaEz+Y59cRBtefUE1l7wmc76giAoyWIA8wgi/LIQHUApi9BHGAGqQVLB6oDMH0J4gAzyYVXJkPDB9eGhnt1AGYUhy8EmEmWX9K7vvnq3vWCZb0QfqAOwIwhiAPMNMsv6V02b07e8vWuuwHgKFmaAgAAHRDEAQCgA4I4AAB0QBAHAIAOCOIAANABQRwAADogiAMAQAcEcQAA6IAgDgAAHRDEAQCgA4I4AAB0QBAHAIAOCOIAANABQRwAADogiAMAQAcEcQAA6IAgDgAAHRDEAQCgA4I4AAB0QBAHAIAOCOIAANABQRwAADogiAMAQAcEcQAA6ECnQbyqPlFV362qr4+qXVVVO6vq7v7lDaPuW1tV91XV9qpaNap+XlVt6993XVXVVL8XAAAYRNdbxD+Z5PVj1K9trb2if/liklTVWUkuTfLy/nM+XFXz+o//SJLLkry0fxnrNQEAYNroNIi31v4kyd+N8+FvTHJja+3x1tq3ktyX5NVV9aIkz22t3dZaa0k+nWT1pDQMAAATpOst4ofz7qra2l+6ckq/tiTJA6Mes6NfW9K/fWgdAACmreptRO6wgarTk/xha+3s/teLkjySpCX5jSQvaq39bFX9hyS3tdY+03/cx5N8McnfJvlAa+21/fqPJfm3rbWLxvhel6W3hCWLFi0678Ybb5zstzdr7N27NyeddFLXbcwoZjYY8xqcmQ3OzAZjXoMzs8HMlXldcMEFd7bWVhxaP76LZo6ktfadA7er6j8m+cP+lzuSLBv10KVJdvXrS8eoj/Xa1ye5PklWrFjRVq5cOWF9z3abN2+OeQ3GzAZjXoMzs8GZ2WDMa3BmNpi5Pq9ptzSlv+b7gDclOXBElY1JLq2qE6vqjPR2yvxaa+3BJI9W1Wv6R0t5W5LPT2nTAAAwoE63iFfV7yVZmeTUqtqR5NeSrKyqV6S3NOX+JJcnSWvtnqpan+QbSZ5M8q7W2v7+S/18ekdgGU7ypf4FAACmrU6DeGvtLWOUP36Ex78/yfvHqN+R5OwJbA0AACbVtFuaAgAAc4EgDjDDbNiyM+df85Vs27kn51/zlWzYsrPrlgA4CtPuqCkAHN6GLTuz9qZtGdm3P1mW7Nw9krU3bUuSrD7XKRQAZhJbxAFmkHWbtud1+7+aW0+4IufUt3LrCVfkdfu/mnWbtnfdGgADEsQBZpAV3/9yrhn6WJYe90hSydLjHsk1Qx/Liu9/uevWABiQIA4wg6w94bN5Tj1xUO059UTWnvDZjjoC4GgJ4gAzyKI8MlAdgOlLEAeYQWrB0oHqAExfgjjATHLhlcnQ8MG1oeFeHYAZxeELAWaS5Zf0rm++une9YFkvhB+oAzBjCOIAM83yS3qXzZuTt3y9624AOEqWpgAAQAcEcQAA6IAgDgAAHRDEAQCgA4I4AAB0QBAHAIAOCOIAANABQRwAADogiAMAQAcEcQAA6IAgDgAAHRDEAQCgA4I4AAB0QBAHAIAOCOIAANABQRwAADogiAMAQAcEcQAA6IAgDgAAHRDEAWaYDVt25vxrvpJtO/fk/Gu+kg1bdnbdEgBHQRAHmEE2bNmZWz/34fz+//yXOae+ld//n/8yt37uw8I4wAwkiAPMIHf/0fW5uq7P0uMeSSpZetwjubquz91/dH3XrQEwIEEcYAZ55xOfyXPqiYNqz6kn8s4nPtNRRwAcLUEcYAZZfNz3BqoDMH0J4gAzyGPDpw1UB2D6EsQBZpDn/MTVeXLe/INqT86bn+f8xNUddQTA0Tq+6wYAGMDyS3p/uG/uB+8Fy3L8hVcmyy/psisAjoIgDjDTLL+kd9m8OXnL17vuBoCjZGkKAAB0QBAHAIAOCOIAANABQRwAADogiAMAQAcEcQAA6IAgDgAAHRDEAQCgA4I4AAB0QBAHAIAOCOIAANABQRwAADogiAMAQAcEcQAA6IAgDgAAHRDEAQCgA4I4AAB0QBAHAIAOCOIAANABQRwAADogiAMAQAcEcQAA6IAgDgAAHRDEAQCgA4I4AAB0QBAHAIAOCOIAANABQRwAADogiAMAQAcEcQAA6IAgDgAAHRDEAQCgA4I4AAB0QBAHAIAOCOIAANCBToN4VX2iqr5bVV8fVXteVX25qv66f33KqPvWVtV9VbW9qlaNqp9XVdv6911XVTXV7wUAAAbR9RbxTyZ5/SG19yS5ubX20iQ3979OVZ2V5NIkL+8/58NVNa//nI8kuSzJS/uXQ18TAACmlU6DeGvtT5L83SHlNyb5VP/2p5KsHlW/sbX2eGvtW0nuS/LqqnpRkue21m5rrbUknx71HAAAmJaO77qBMSxqrT2YJK21B6vqhf36kiR/PupxO/q1ff3bh9afoaouS2/LeRYtWpTNmzdPbOez2N69e81rQGY2GPManJkNzswGY16DM7PBzPV5TccgfjhjrftuR6g/s9ja9UmuT5IVK1a0lStXTlhzs93mzZtjXoMxs8GY1+DMbHBmNhjzGpyZDWauz6vrNeJj+U5/uUn619/t13ckWTbqcUuT7OrXl45RBwCAaWs6BvGNSd7ev/32JJ8fVb+0qk6sqjPS2ynza/1lLI9W1Wv6R0t526jnAADAtNTp0pSq+r0kK5OcWlU7kvxakmuSrK+qn0vyt0l+Kklaa/dU1fok30jyZJJ3tdb291/q59M7Astwki/1LwAAMG11GsRba285zF0XHubx70/y/jHqdyQ5ewJbAwCASTUdl6YAAMCsJ4gDAEAHBHEAAOiAIA4AAB0QxAEAoAOCOAAAdEAQBwCADgjiAADQAUEcAAA6IIgDAEAHBHEAAOiAIA4AAB0QxAEAoAOCOAAAdEAQBwCADgjiAADQAUEcAAA6IIgDAEAHBHEAAOiAIA4AAB0QxAEAoAOCOAAAdEAQBwCADgjiAADQAUEcAAA6IIgDAEAHBHEAAOiAIA4AAB0QxAEAoAOCOAAAdEAQBwCADgjiAADQAUEcAAA6IIgDAEAHBHEAAOiAIA4AAB0QxAEAoAOCOAAAdEAQBwCADgjiAADQAUEcAAA6IIgDAEAHBHEAAOiAIA4AAB0QxAEAoAPHj+dBVXVKksVJRpLc31r7waR2BQAAs9xhg3hVLUjyriRvSXJCkoeTzE+yqKr+PMmHW2u3TEmXAAAwyxxpi/gfJPl0kh9rre0efUdVnZfkZ6rqxa21j09ifwAAMCsdNoi31l53hPvuTHLnpHQEAABzwHjXiC9Pcvrox7fWbpqkngAAYNZ71iBeVZ9IsjzJPUkO7KTZkgjiAABwlMazRfw1rbWzJr0TAACYQ8ZzHPHbqkoQBwCACTSeLeKfSi+MP5Tk8SSVpLXWlk9qZwAAMIuNJ4h/IsnPJNmWp9eIAwAAx2A8QfxvW2sbJ70TAACYQ8YTxP+qqv5zki+ktzQlicMXAgDAsRhPEB9OL4D/01E1hy8EAIBj8KxBvLX2f05FIwAAMJeM54Q+ZyT5hTzzzJoXT15bAAAwu41nacqGJB9Pb424o6YAAMAEGE8Qf6y1dt2kdwIAAHPIeIL471TVryX5rzn4qCl3TVpXAAAwy40niJ+T3gl9fjxPL01p/a8BAICjMJ4g/qYkL26tPTHZzQAAwFxx3Dge85dJFk5yHwAAMKeMZ4v4ovTOrnl7Dl4j7vCFAABwlMYTxH9t0rsAAIA55rBBvKqq9Xz12R4zOa0BAMDsdaQ14rdU1S9U1T8YXayqE6rqx6vqU0nePrntAQDA7HSkpSmvT/KzSX6vf5r73UnmJ5mX3jHFr22t3T3ZDQIAwGx02CDeWnssyYeTfLiqhpKcmmSktbZ7inoDAIBZazw7a6a1ti/Jg5PcCwAAzBnjOY44AAAwwaZtEK+q+6tqW1XdXVV39GvPq6ovV9Vf969PGfX4tVV1X1Vtr6pV3XUOAADP7lmDeFW9e3TgnWIXtNZe0Vpb0f/6PUlubq29NMnN/a9TVWcluTTJy9PbyfTDVTWvi4YBAGA8xrNF/LQkt1fV+qp6fVXVZDd1BG9M8qn+7U8lWT2qfmNr7fHW2reS3Jfk1VPfHgAAjE+N53w8/fD9T5P8n0lWJFmf5OOttb+ZtMaqvpXk75O0JB9trV1fVbtbawtHPebvW2unVNWHkvx5a+0z/frHk3yptfYHh7zmZUkuS5JFixadd+ONN05W+7PO3r17c9JJJ3XdxoxiZoMxr8GZ2eDMbDDmNTgzG8xcmdcFF1xw56gVHk8Z71FTWlU9lOShJE8mOSXJH1TVl1tr/3ZiW33K+a21XVX1wiRfrqq/OsJjx9pK/4z/w2itXZ/k+iRZsWJFW7ly5YQ0Ohds3rw55jUYMxuMeQ3OzAZnZoMxr8GZ2WDm+rzGs0b8iqq6M8lvJfmzJOe01n4+yXlJ/vfJaqy1tqt//d0kn0tvqcl3qupF/b5elOS7/YfvSLJs1NOXJtk1Wb0BAMCxGs8a8VOTvLm1tqq19tn+McXTWvtBkn82GU1V1Q9V1ckHbqe3LObrSTYmeXv/YW9P8vn+7Y1JLq2qE/tnAX1pkq9NRm8AADARnnVpSmvtyiPcd+/EtvOURUk+198v9Pgk/7m19v9W1e1J1lfVzyX52yQ/1e/jnqpan+Qb6S2deVdrbf8k9QYAAMdsXGvEp1pr7ZtJ/skY9e8lufAwz3l/kvdPcmsAADAhpu0JfQAAYDYTxAEAoAOCOAAAdEAQBwCADgjiAADQAUEcAAA6IIgDzDC3b/xoHrrqJcmuu/PQVS/J7Rs/2nVLABwFQRxgBrl940dz9p3vy2l5OKnktDycs+98nzAOMAMJ4gAzyLK71mW4njioNlxPZNld6zrqCICjJYgDzCAvbA8fpv7IFHcCwLESxAFmkO/WCw5TP3WKOwHgWAniADPIA69ck5F2wkG1kXZCHnjlmo46AuBoHd91AwCM36suvjy3p7dWPC15KC/IA+etyasuvrzr1gAYkCAOMMO86uLLk4svz19t3pzT/sV9Oa3rhgA4KpamAABABwRxAADogCAOAAAdEMQBAKADgjgAAHRAEAcAgA4I4gAA0AFBHAAAOiCIAwBABwRxAADogCAOAAAdEMQBAKADgjgAAHRAEAcAgA4I4gAA0AFBHAAAOiCIAwBABwRxAADogCAOAAAdEMQBAKADgjgAAHRAEAcAgA4I4gAA0AFBHAAAOiCIAwBABwRxgBnm9o0fzUNXvSTZdXceuuoluX3jR7tuCYCjIIgDzCC3b/xozr7zfTktDyeVnJaHc/ad7xPGAWYgQRxgBll217oM1xMH1YbriSy7a11HHQFwtARxgBnkhe3hw9QfmeJOADhWgjjADPLdesFh6qdOcScAHCtBHGAGeeCVazLSTjioNtJOyAOvXNNRRwAcreO7bgCA8XvVxZfn9vTWiqclD+UFeeC8NXnVxZd33RoAAxLEAWaYV118eXLx5fmrzZtz2r+4L6d13RAAR8XSFAAA6IAgDgAAHRDEAQCgA4I4AAB0QBAHAIAOCOIAANABQRwAADogiAMAQAcEcQAA6IAgDgAAHRDEAQCgA4I4AAB0QBAHAIAOCOIAANABQRwAADogiAMAQAcEcQAA6IAgDgAAHRDEAQCgA4I4AAB0QBAHAIAOCOIAANABQRwAADogiAMAQAcEcQAA6IAgDgAAHZg1QbyqXl9V26vqvqp6T9f9AADAkcyKIF5V85L8hyQ/keSsJG+pqrO67QoAAA5vVgTxJK9Ocl9r7ZuttSeS3JjkjR33BAAAh1Wtta57OGZV9c+TvL619s7+1z+T5H9prb37kMddluSyJFm0aNF5N95445T3OlPt3bs3J510UtdtzChmNhjzGpyZDc7MBmNegzOzwcyVeV1wwQV3ttZWHFo/votmJkGNUXvG/2G01q5Pcn2SrFixoq1cuXKS25o9Nm/eHPMajJkNxrwGZ2aDM7PBmNfgzGwwc31es2Vpyo4ky0Z9vTTJro56AQCAZzVbgvjtSV5aVWdU1QlJLk2yseOeAADgsGbF0pTW2pNV9e4km5LMS/KJ1to9HbcFAACHNSuCeJK01r6Y5Itd9wEAAOMxW5amAADAjCKIAwBABwRxAADogCAOAAAdEMQBAKADgjgAAHRAEAcAgA4I4gAA0AFBHAAAOjBrzqwJ0JUNW3Zm3abt2bV7JIsXDmfNqjOz+twlXbcFwDRnizjAMdiwZWfW3rQt533/y/nTE67In468Ka/a8L/m9o0f7bo1AKY5W8QBjsG6Tdvzuv1fzTVDH8tz6okkyZI8kufd9avJ6ackyy/puEMApitbxAGOwa7dI/m3x69/KoQfMJzHk5uv7qgrAGYCQRzgGCxeOJzF9cjYd+7ZMbXNADCjCOIAx2DNqjPzYE4d+84FS6e2GQBmFEEc4BisPndJdp33bzOSEw++Y2g4ufDKbpoCYEYQxAGO0asuvjzDb/5QsmBZkupdX3SdHTUBOCJHTQGYCMsvEbwBGIgt4gAA0AFBHIDubV2fXHt2ctXC3vXW9V13BDDpLE0BoFtb1ydfuCLZN9L7es8Dva8Ty32AWc0WcQC6dfPVT4fwA/aNOCESMOsJ4gB063AnPnJCJGCWE8QB6NbhTnzkhEjALCeIA9Cp2//RLzghEjAnCeIAdGbDlp152+3/ML/8xM9lxw9OzQ9aZWc7Nbef8+t21ARmPUdNAaAz6zZtz8i+/dmYH83GJ370qfqSbwznzy7usDGAKWCLOACd2bV7ZKA6wGwiiAPQmcULhweqA8wmgjgAnVmz6swMD807qDY8NC9rVp3ZUUcAU8cacQA6s/rcJUl6a8V37R7J4oXDWbPqzKfqALOZIA5Ap1afu0TwBuYkS1MAAKADgjgAAHRAEAcAgA4I4gAA0AFBHAAOtXV9cu3ZyVULe9db13fdETALOWoKAIy2dX3yhSuSff2ze+55oPd1kiy/pLu+gFnHFnEAGO3mq58O4QfsG+nVASaQIA4Ao+3ZMVgd4CgJ4gAw2oKlg9UBjpIgDgCjXXhlMjR8cG1ouFcHmECCOACMtvyS5KLrkgXLklTv+qLr7KgJTDhHTQGAQy2/RPAGJp0gDsCct2HLzqzbtD27do9k8cLhrFl1Zlafu6TrtoBZThAHYE7bsGVn1t60LSP79idJdu4eydqbtiWJMA5MKmvEAZjT1m3a/lQIP2Bk3/6s27S9o46AuUIQB2BO27V7ZKA6wEQRxAGY0xYvHB6oDjBRBHEA5rQ1q87M8NC8g2rDQ/OyZtWZHXUEzBV21gRgTjuwQ6ajpgBTTRAHYM5bfe4SwRuYcpamAABAB2wRB56Vk50AwMQTxIEjcrITAJgclqYAR+RkJwAwOQRx4Iic7AQAJocgDhzRWCc1ufi4W3Pb/F9MrlqYXHt2snX91DcGc93W9b3fP7+HMGMJ4sARHXqyk4uPuzUfHPpYTsvDSVqy54HkC1cIATCVtq7v/d7teSB+D2HmEsSBI1p97pJ84M3nZMnC4VSS957w2QzXEwc/aN9IcvPVnfQHc9LNV/d+70bzewgzjqOmAM/qoJOdXPXWsR+0Z8fUNQRz3eF+3/wewoxiizgwmAVLB6sDE8/vIcwKgjgwmAuvTIYO2YFzaLhXB6aG30OYFQRxYDDLL0kuui5ZsCxJ9a4vuq5XB6aG30OYFawRBwa3/BL/4EPX/B7CjGeLOAAAdEAQBwCADliaArPQhi07s27T9uzaPZLFC4ezZtWZTx9+EACYFgRxmGU2bNmZtTdty8i+/UmSnbtHsvambUkijAPANGJpCswy6zZtfyqEHzCyb3/WbdreUUfAsdqwZWfOv+YrOeM9f5Tzr/lKNmzZObHfYOv65Nqzk6sW9q63rp/Y1wfGZIs4zDK7do8MVAemt0n/lGvr+uQLVyT7+n8j9jzQ+zpJ8sJjf33gsGwRh1lm8cLhgerA9Dbpn3LdfPXTIfyAfSO9OjCpBHGYZdasOjPDQ/Oe+vri427Nn514RW597M0+coYZaNI/5dqzY7A6MGEsTYFZ5sBH1es2bc+K738515zw8Qzn8d6doz9ydiIQmBEWLxzOzjFC94R9yrVgae9vw1h1YFJNuy3iVXVVVe2sqrv7lzeMum9tVd1XVduratWo+nlVta1/33VVVd10D9PD6nOX5M/e8+P5nRd84ekQfoCPnGFGOfRTriQZHpqXNavOnJhvcOGVydAhoX5ouFcHJtW0C+J917bWXtG/fDFJquqsJJcmeXmS1yf5cFUd+Mv0kSSXJXlp//L6DnqG6cdHzjDjrT53ST7w5nOyZOFwKsmShcP5wJvPmbjDkS6/JLnoumTBsiTVu77oOp+awRSYSUtT3pjkxtba40m+VVX3JXl1Vd2f5LmttduSpKo+nWR1ki911ShMGz5yhllh9blLJvc8AMsvEbyhA9Va67qHg1TVVUnekeT7Se5I8kuttb+vqg8l+fPW2mf6j/t4emH7/iTXtNZe26//WJJfbq39szFe+7L0tpxn0aJF5914442T/n5mi7179+akk07quo0ZZVrMbOTve0G8/eDpWh3X2+I1fEp3fY1hWsxrhjGzwZnZYMxrcGY2mLkyrwsuuODO1tqKQ+udbBGvqj9OctoYd/1KestMfiNJ61//+yQ/m2Ssdd/tCPVnFlu7Psn1SbJixYq2cuXKQVufszZv3hzzGsy0mdnW9b014Xt29LaEX3hlsvxNXXf1DNNmXjOImQ3OzAZjXoMzs8HM9Xl1EsQPbL1+NlX1H5P8Yf/LHUmWjbp7aZJd/frSMepA4iNnAJimpt3OmlX1olFfvinJ1/u3Nya5tKpOrKoz0tsp82uttQeTPFpVr+kfLeVtST4/pU0DAMCApuPOmr9VVa9Ib3nJ/UkuT5LW2j1VtT7JN5I8meRdrbUDpxr7+SSfTDKc3rpxO2oCwBg2bNmZdZu2Z9fukSxeOJw1q86c3B1BgcOadkG8tfYzR7jv/UneP0b9jiRnT2ZfADDTbdiyM2tv2paRfb3tWDt3j2TtTduSRBiHDky7pSkAwORYt2n7UyH8gJF9+7Nu0/aOOoK5TRAHgDli1+6RgerA5BLEAWCOWLxweKA6MLkEcQCYI9asOjPDQ/MOqg0PzcuaVWd21BHMbdNuZ00AYHIc2CHTUVNgehDEAWAOWX3uEsEbpglLUwAAoAOCOAAAdEAQBwCADgjiAADQAUEcAAA6IIgDAEAHHL4QjmDDlp2OtwsATApBHA5jw5adWXvTtozs258k2bl7JGtv2pYkwjgAcMwsTYHDWLdp+1Mh/ICRffuzbtP2jjoCAGYTQRwOY9fukYHqAACDEMThMBYvHB6oDgAwCEEcDmPNqjMzPDTvoNrw0LysWXVmRx0BALOJnTXhMA7skOmoKQDAZBDE4QhWn7vk4OC9dX1y7dXJnh3JgqXJhVcmyy/prkEAYMYSxGG8tq5PvnBFsq+/s+aeB3pfJ8I4ADAwa8RhvG6++ukQfsC+kV4dgLlh6/rk2rOTqxb2rreu77ojZjBbxGG89uwYrA7A7OKTUSaYLeIwXguWDlYHYHbxySgTTBCH8brwymTokGOIDw336gDMfj4ZZYIJ4jBeyy9JLrouWbAsSfWuL7rOx5EAc4VPRplg1ojDIJZfIngDzFUXXnnwGvHEJ6McE1vEAQDGwyejTDBbxAEAxssno0wgW8QBAKADgjgAAHRAEAcAgA4I4gAA0AFBHAAAOiCIAwBABwRxAADogCAOAAAdEMQBAKADzqzJlNuwZWfWbdqeXbtHsnjhcNasOjOrz13SdVsAAFNKEGdKbdiyM2tv2paRffuTJDt3j2TtTduSRBgHAOYUS1OYUus2bX8qhB8wsm9/1m3a3lFHAADdEMSZUrt2jwxUBwCYrSxNYUotXjicnWOE7sULhzvoBgCeZh8mppot4kypNavOzPDQvINqw0PzsmbVmR11BABP78O0c/dIWp7eh2nDlp1dt3Z0tq5Prj07uWph73rr+q47YgyCOFNq9blL8oE3n5MlC4dTSZYsHM4H3nyOLQ4AdGpW7cO0dX3yhSuSPQ8kab3rL1whjE9DlqYw5Vafu0TwBmBamVX7MN18dbLvkL73jfTqyy/ppifGZIs4ADDnHW5fpRm5D9OeHYPV6YwgzvRgLRsAHZpV+zAtWDpYnc5YmsJRmdA9yw+sZTvwMdqBtWyJj9AAmBIH/g2bFUdNufDKg/9dTZKh4V6daUUQZ2ATfnZMa9kAmAZmzT5MB/7tvPnq3nKUBUt7Idy/qdOOIM7AjrRn+VH9AbOWDQAm1vJLBO8ZwBpxBjbhe5ZbywYAzEGCOAOb8D3LL7yyt3ZtNGvZAIBZThBnYBO+Z/nyS5KLrksWLEtSveuLrvORGgAwq1kjzsAmZc9ya9kAgDlGEOeozJo9ywGAwWxd74gsE0QQBwBgfJz7Y0JZIw4AwPgc6dwfDEwQZ2I4RT0AzH7O/TGhBHGO3YGPqfY8kKQ9/TGVMA4As4tzf0woQZxj52MqAJgbnPtjQgniHDsfUwHA3ODcHxPKUVNmuQ1bdk7s8b7HsmBpf1nKGHUAYHZx7o8JY4v4LLZhy86svWlbdu4eSUuyc/dI1t60LRu27JzYb+RjKgCAgQnis9i6Tdszsm//QbWRffuzbtP2if1GPqYCABiYpSmz2K7dIwPVj4mPqQAABmKL+Cy2eOHwQHUAAKaOID6LrVl1ZoaH5h1UGx6alzWrzuyoIwAADrA0ZRY7cHSUST9qCgAw6abkSGhMKUF8llt97hK/pAAwwx04EtqBgzAcOBJaEv/Oz2CWpgAATHNTdiQ0ppQgDgAwzU3pkdCYMoI4AMA050hos1MnQbyqfqqq7qmqH1TVikPuW1tV91XV9qpaNap+XlVt6993XVVVv35iVf1+v/4XVXX6FL8dAIBJ5Uhos1NXW8S/nuTNSf5kdLGqzkpyaZKXJ3l9kg9X1YGfuo8kuSzJS/uX1/frP5fk71trL0lybZIPTnr3s83W9cm1ZydXLexdb13fdUcAwCirz12SD7z5nCxZOJxKsmThcD7w5nPsqDnDdXLUlNbavUnS36g92huT3NhaezzJt6rqviSvrqr7kzy3tXZb/3mfTrI6yZf6z7mq//w/SPKhqqrWWpvktzE7bF2ffOGKZF9/jdmeB3pfJ86UCQDTiCOhzT7TbY34kiQPjPp6R7+2pH/70PpBz2mtPZlkT5LnT3qns8XNVz8dwg/YN9KrAwAwaSZti3hV/XGS08a461daa58/3NPGqLUj1I/0nLF6uiy95S1ZtGhRNm/efJg25pDT3jn2f6UkGTWfvXv3mteAzGww5jU4MxucmQ3GvAZnZoOZ6/OatCDeWnvtUTxtR5Jlo75emmRXv750jPro5+yoquOTLEjyd4fp6fok1yfJihUr2sqVK4+ixVnm2nf3lqMcasGy5C1ff+rLzZs3x7wGY2aDMa/BmdngzGww5jU4MxvMlMxr6/reJ/17diQLliYXXjltlt9Ot6UpG5Nc2j8Syhnp7ZT5tdbag0kerarX9I+W8rYknx/1nLf3b//zJF+xPnwAF16ZDB1y6KOh4V4dAGAmO7Av3J4HkrSn94WbJgem6GRnzap6U5L/J8kLkvxRVd3dWlvVWrunqtYn+UaSJ5O8q7V24DRSP5/kk0mG09tJ80v9+seT/Kf+jp1/l95RV6alDVt2Zt2m7dm1eySLFw5nzaozu9/p4sD/EU7T/1MEADhqR9oXbhpkna6OmvK5JJ87zH3vT/L+Mep3JDl7jPpjSX5qonucaBu27Mzam7Y9dXranbtHsvambUkyPcL4NPhhBACm3rTcUDhR9uwYrD7FptvSlFlr3abtT4XwA0b27c+6Tds76ggAmOsObCjcuXskLU9vKNywZWfXrU2MBUsHq08xQXyK7No9MlAdAGCyzfoNhdN8XzhBfIosXjg8UB0AYLLN+g2Fyy9JLrqudzS4VO/6ouumzZLcTtaIz0VrVp150BrxJBkempc1q87ssCsAYC5bvHA4O8cI3bNqQ+E03hfOFvEpsvrcJfnAm8/JkoXDqSRLFg7nA28+Z/bsDAEAzDhrVp2Z4aF5B9VsKJw6tohPodXnLhG8AYBp40AumbVHTZnmBHEAgDnMhsLuWJoCAAAdEMQBAKADgvh0sHV9cu3ZyVULe9db13fdEQAAk8wa8a5tXZ984YpkX//QQXse6H2dTNtD7QAAcOxsEe/azVc/HcIP2DfSqwMAMGsJ4l3bs2OwOgAAs4Ig3rUFSwerAwAwKwjiXbvwymTokNPIDg336gAAzFqCeNeWX5JcdF2yYFmS6l1fdJ0dNQEAZjlHTZkOll8ieAMAzDG2iAMAQAcEcQAA6IAgDgAAHRDEAQCgA4I4AAB0QBAHAIAOCOIAANABQRwAADrghD4AAEx7G7bszLpN27Nr90gWLxzOmlVnZvW5S7pu65gI4gAATGsbtuzM2pu2ZWTf/iTJzt0jWXvTtiSZ0WHc0hQAAKa1dZu2PxXCDxjZtz/rNm3vqKOJIYgDADCt7do9MlB9phDEAQCY1hYvHB6oPlMI4gAATGtrVp2Z4aF5B9WGh+ZlzaozO+poYthZEwCAae3ADpmOmgIAAFNs9blLZnzwPpSlKQAA0AFBHAAAOiCIAwBABwRxAADogCAOAAAdEMQBAKADgjgAAHRAEAcAgA4I4gAA0AFBHAAAOiCIAwBABwRxAADogCAOAAAdEMQBAKADgjgAAHRAEAcAgA4I4gAA0AFBHAAAOiCIAwBABwRxAADogCAOAAAdEMQBAKADgjgAAHRAEAcAgA4I4gAA0IFqrXXdQyeq6uEk3+66jxnk1CSPdN3EDGNmgzGvwZnZ4MxsMOY1ODMbzFyZ1z9srb3g0OKcDeIMpqruaK2t6LqPmcTMBmNegzOzwZnZYMxrcGY2mLk+L0tTAACgA4I4AAB0QBBnvK7vuoEZyMwGY16DM7PBmdlgzGtwZjaYOT0va8QBAKADtogDAEAHBHEAAOiAIE6SpKp+qqruqaofVNWKUfXTq2qkqu7uX3531H3nVdW2qrqvqq6rqurXT6yq3+/X/6KqTu/gLU2qw82rf9/a/nvfXlWrRtXn7LwOVVVXVdXOUT9Xbxh130Dzm4uq6vX9+dxXVe/pup/ppKru7/+c3F1Vd/Rrz6uqL1fVX/evTxn1+DF/3marqvpEVX23qr4+qjbwfObS7+NhZuZv2GFU1bKquqWq7u3/O/mL/bqfs7G01lxckuSHk5yZZHOSFaPqpyf5+mGe87Uk/58kleRLSX6iX//XSX63f/vSJL/f9fubwnmdleQvk5yY5Iwkf5Nk3lyf1xjzuyrJ/2+M+sDzm2uXJPP6c3lxkhP68zqr676myyXJ/UlOPaT2W0ne07/9niQf7N8+7M/bbL0k+V+TvHL03/Wjmc9c+n08zMz8DTv8vF6U5JX92ycn+e/9ufg5G+NiizhJktbava217eN9fFW9KMlzW2u3td5vy6eTrO7f/cYkn+rf/oMkF862/4s9wrzemOTG1trjrbVvJbkvyavn+rwGcDTzm2teneS+1to3W2tPJLkxvblxeKN/xz6Vg3/3nvHzNvXtTZ3W2p8k+btDygPNZ679Ph5mZocz52fWWnuwtXZX//ajSe5NsiR+zsYkiDMeZ1TVlqr6alX9WL+2JMmOUY/Z0a8duO+BJGmtPZlkT5LnT1WzHXvqvfcdmIt5PdO7q2pr/2PfAx9RHs385prDzYieluS/VtWdVXVZv7aotfZg0gsJSV7Yr5tlz6Dz8fvY42/Ys+gvtTw3yV/Ez9mYju+6AaZOVf1xktPGuOtXWmufP8zTHkzyD1pr36uq85JsqKqXp/cx0aEOHAvzSPfNGEc5r8O991k/r0MdaX5JPpLkN9J7n7+R5N8n+dkc3fzmGrM4svNba7uq6oVJvlxVf3WEx5rlkfl9PDx/w55FVZ2U5L8k+Tette8f4YPeOT0zQXwOaa299iie83iSx/u376yqv0nyj9P7P9Olox66NMmu/u0dSZYl2VFVxydZkPF/rDdtHM288vR7P+DAXGb9vA413vlV1X9M8of9L49mfnPN4WZEktbarv71d6vqc+ktNflOVb2otfZg/+Pu7/YfbpY9g85nzv8+tta+c+C2v2HPVFVD6YXwG1prN/XLfs7GYGkKR1RVL6iqef3bL07y0iTf7H+s9GhVvaa/nvltSQ5sJd6Y5O392/88yVf667vmgo1JLq3ekVDOSG9eXzOvg/X/CB/wpiQHjkZwNPOba25P8tKqOqOqTkhvB9+NHfc0LVTVD1XVyQduJ/mn6f1sjf4de3sO/t17xs/b1HY9LQw0H7+P/oYdSf/9fTzJva213x51l5+zsXS9t6jL9Lik94dkR3pbv7+TZFO//r8nuSe9PZrvSnLRqOesSO+Pz98k+VCePlPr/CSfTW+Hi68leXHX72+q5tW/71f6M9meUXt4z+V5jTG//5RkW5Kt6f0RftHRzm8uXpK8Ib0jEfxNekulOu9pOlzSO5LMX/Yv9xyYTXr7XNyc5K/7188b9Zwxf95m6yXJ76W35HBf/2/Yzx3NfObS7+NhZuZv2OHn9aPpLSHZmuTu/uUNfs7GvjjFPQAAdMDSFAAA6IAgDgAAHRDEAQCgA4I4AAB0QBAHAIAOCOIAPKWqhqvqqwfOHzDBr/3Ho04FDjDnCeIAjPazSW5qre2fhNf+T0n+9SS8LsCMJIgDzAFV9aqq2lpV8/tnoLynqs4e46FvTf/sdVW1sr91fH1V/fequqaq3lpVX6uqbVX1j/qP+2RVfaSqbqmqb1bV/1ZVn6iqe6vqk6Nee2OSt0z6mwWYIY7vugEAJl9r7faq2pjk3yUZTvKZ1trXRz+mqk5I78yu948q/5MkP5zk75J8M8nHWmuvrqpfTPILSf5N/3GnJPnxJBcn+UKS85O8M8ntVfWK1trdrbW/75/G+vmtte9N1nsFmClsEQeYO65O8rr0Thv9W2Pcf2qS3YfUbm+tPdhaezy900z/1359W5LTRz3uC613quZtSb7TWtvWWvtBeqeaH/247yZZfGxvA2B2EMQB5o7nJTkpyclJ5o9x/8gY9cdH3f7BqK9/kIM/VX18jMeM9bj5/e8DMOcJ4gBzx/VJfjXJDUk+eOidrbW/TzKvqsYK6cesqirJaUnun4zXB5hpBHGAOaCq3pbkydbaf05yTZJXVdWPj/HQ/5rkRyepjfOS/Hlr7clJen2AGaV6S/oAIKmqc5P8f1trPzMJr/07STa21m6e6NcGmIlsEQfgKa21LUlumYwT+iT5uhAO8DRbxAEAoAO2iAMAQAcEcQAA6IAgDgAAHRDEAQCgA4I4AAB04P8PxaqRcMQCFiUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x864 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(12,12))\n",
    "\n",
    "ax.scatter(pos_i[:,0], pos_i[:,1], label='initial')\n",
    "ax.scatter(pos[:,0], pos[:,1], label='final')\n",
    "\n",
    "ax.grid(True)\n",
    "ax.set_title('Final Solution')\n",
    "ax.set_xlabel('x (mm)')\n",
    "ax.set_ylabel('y (mm)')\n",
    "ax.legend()\n",
    "plt.savefig('final_positions.png')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAALJCAYAAACp99XTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA9e0lEQVR4nO3dfZxdVX0v/s8KEXTkQSImDcRkWosQEzACFfTnjSANitIqoFUcWzFIbm1tVRSF5vYq3o6ilUJUqh0bFHVKqoJIkUYUjFqVUpCAgCKoeQIEkccQnlm/P85J3BMSzMPMnJnM+/165bXPXnudvb9z1iF8smadfUqtNQAAQMu4ThcAAAAjiYAMAAANAjIAADQIyAAA0CAgAwBAg4AMAAANAjIw7Eopny6l/H1j/22llNtKKatLKc8spfx/pZQb2/uv6WCpg660fLaUclcp5fIOXH/Aa781fUspf1dK+ddNPNfnSin/sKl1bq1SytT2+2e7J+mzupTyB8NVEzB6FPdBBgZTKWVZkklJHk3yWJLrk3w+SV+t9fEN9H9KknuTHFRrvbrddkmSC2qtC4ar7uFSSvlfSc5Jslet9f5O17OpSikHJ/lirXXKFj7/c0lW1Vr/zyCWtTnXX5JW/ZsU6IGxzQwyMBT+pNa6U5JpSU5N8r4kCzfSd1KSpya5rtE2bb39TVZKGb8lzxtG05Is21g4HgX1A2zzBGRgyNRa76m1XpDk9UneXEqZmfz21+2llOcmuaHd/e5SyqWllJ8n+YMk/9H+FfgOpZRdSikLSym3llJubj93u/a5ji2lfL+Ucnop5c4kH2g/52OllBXtpRufLqU8rd3/4FLKqlLKu0spt7fP+Za1NZdSnlZKOa2UsryUck8p5b8azz2olPKDUsrdpZSr27OqadTxi1LKfaWUX5ZSetZ/PUopxyX51yQvav9spzTqeV8p5VdJPtuu/4xSyi3tP2eUUnZYr/73Nup/TSnllaWUn5VS7iyl/N3GxqS51GETXou14/T0JP+ZZPd23atLKbuXUj5QSvlio/+XSym/ar9u3y2lzNiU90ljDD/Rfu5PSymHNo7vXkq5oP2z3VRKOb5x7IWllCtKKfe2x/qf2u3dpZRaShlfSulN8r+SfLJd+yfbfWop5Q/bj3cppXy+lPLr9tj/n1LKuEZ9/9V+T93VHt/DN+VnA0YnARkYcrXWy5OsSiukNNt/lmRtiHpGrfVltdbnJFmR1iz0jrXWh5KcndaSjT9M8oIkhyV5a+NUByb5RZKJSXqTfCTJc5PMaj9njyT/t9H/95Ls0m4/LsmZpZRd28c+lmT/JC9OMiHJe5M8XkrZI8nXk/xDu/09Sc4tpTyrHSA/nuTw9sz5i5Ms3cDrsDDJXyb5Yftne3+jnglpzS7PSzI/yUHt+p+f5IVJmksTfi+tWfe1P9dnkrypXff/SvJ/y6avrX2y12Jt3fcnOTzJLe26d6y13rKBc/1nkj3TGocfJenfxBqS347hbknen+S8UsqE9rFz0nr/7J7ktUk+1AjQC5IsqLXunOQ5Sb60/olrrfOTfC/J29u1v30D1/9EWq/DHyR5aZK/SPKWxvED0/rH3G5JPppkYSmlbMbPB4wiAjIwXG5JKwRullLKpLTC2TtrrffXWm9PcnqSNzTPXWv9RK310SQPJjk+ybtqrXfWWu9L8qH1+j+S5IO11kdqrRclWZ1kr/aM4dwk76i13lxrfazW+oN2SH9TkotqrRfVWh+vtX4zyRVJXtk+5+NJZpZSnlZrvbXWujlLRB5P8v5a60O11geS9LTru73W+uskpyT58/Xq7621PpJkUVqhbUGt9b72da9Lsu8mXnuDr8Vm1L5OrfWsdg0PJflAkueXUnbZxKffnuSMdh3/nlYYfVUp5dlJXpLkfbXWB2utS9OahV/7ejyS5A9LKbvVWlfXWi/b3LpL67cRr09ycrv+ZUlOy8DXfHmt9TO11sfS+gfb5LSWBwHbIAEZGC57JLlzC543LclTktzaXtpwd5J/SWuWcq2VjcfPStKV5MpG/8Xt9rV+0w7Ta61JsmNaQfOpSX6+kTpet/ac7fO+JMnk9gzr69OaHb61lPL1Usrem/Ez/rrW+mBjf/ckyxv7y9ttzfofaz9+oL29rXH8gfbPsyk29lpsllLKdqWUU0spPy+l3JtkWfvQbpt4ipvrwE+Nr/2Zd0+y9h86zWN7tB8fl9ZvC35aSvmfUsoRm1t7u8bt88TXfI/G/q/WPqi1rmk/3OzXCRgdBGRgyJVS/iitsPFfW/D0lUkeSrJbrfUZ7T8711qb61ubweqOtALijEb/XWqtmxJm7khrBvo5G6njC41zPqPW+vRa66lJUmv9Rq11Tloziz9Na9nDplr/dkK3pBXI15rabuuk33XLozcmeXWSP05rqUJ3u31TlyHssd6ShbU/8y1JJpRSdlrv2M1JUmu9sdZ6TFr/YPpIkq+0l7xsTv13pDUTvf5rfvMm1g5sYwRkYMiUUnZuz+gtSusWWz/e3HPUWm9NcnGS09rnG1dKeU4p5aUb6f94WuH09FLKxHYde5RSXr4J13o8yVlJ/qn9wbDtSikvan9A7otJ/qSU8vJ2+1PbH3KbUkqZVEr503YweyitZQqPPcmlfpdzkvyf9vrm3dJaZ/zF3/GcoXZbkmc+yZKJndL62X+T1gz+hzbz/BOT/G0p5SmllNclmZ7WkpaVSX6Q5MPt13zftGaN+5OklPKmUsqz2mN3d/tcG3rtb0trffETtGfjv5Skt5SyUyllWpIT0vnXHOgQARkYCv9RSrkvrVnX+Un+KQM/8LS5/iKtX4Ffn+SuJF9Ja6Z2Y96X5KYkl7V/3f+tbPq62vck+XGS/0lrSchHkoxrB7VXJ/m7JL9O62c7Ma2/R8cleXdas513pvUhr7/a9B/vCf4hrfXN17Rr+VG7rWNqrT9NK7j/or3EZPf1unw+rWUJN6c1Tpu7Fvi/0/qA3x1pfdDytbXW37SPHZPWjPQtSb6a1nrtb7aPvSLJdaWU1Wl9YO8N6y1XWWtBkte270Lx8Q0c/5sk96f1QcH/SvJvaf1jCRiDfFEIAB1VSjk2yVtrrS/pdC0AiRlkAAAYQEAGAIAGSywAAKDBDDIAADSM73QBTbvttlvt7u7udBmb7P7778/Tn76h220yFhj/sc34j23GH++B0e/KK6+8o9b6rA0dG1EBubu7O1dccUWny9hkS5YsycEHH9zpMugQ4z+2Gf+xzfjjPTD6lVKWb+yYJRYAANAgIAMAQIOADAAADQIyAAA0CMgAANAgIAMAQIOADAAADQIyAAA0CMgAANAgIAMAQIOADAAADQIyAAA0CMgAANAgIAMAQIOADAAADQIyAAA0CMgAANAgIAMAQIOADAAADQIyAAA0CMgAANAgIAMA0FFz587NxIkTM3PmzHVtr3/96zNr1qzMmjUr3d3dmTVr1rDVM37YrgQAABtw7LHH5u1vf3v+4i/+Yl3bv//7v697/O53vzu77LLLsNUjIAMA0FGzZ8/OsmXLNnis1povfelLufTSS4etHkssAAAYsb73ve9l0qRJ2XPPPYftmgIyAADDq78/6e5Oxo1rbfv7N9r1nHPOyTHHHDNspSWWWAAAMJz6+5N585I1a1r7y5e39nt7n9D10UcfzXnnnZcrr7xyWEs0gwwAwPCZP/+34XitNWuSj33sCV2/9a1vZe+9986UKVOGqbgWARkAgOGzYsUTmo5J8qKbb84NN9yQKVOmZOHChUmSRYsWDfvyisQSCwAAhtPUqa1lFQ3nJMm0acl6d7L43Oc+N1xVDWAGGQCA4dPbm3R1DWzr6trgGuROEZABABg+PT1JX19rxriU1ravr9U+QlhiAQDA8OrpGVGBeH1mkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAKBh7ty5mThxYmbOnDmg/ROf+ET22muvzJgxI5/+9Kc7VB3DQUAGAGg49thjs3jx4gFt3/72t/O1r30t11xzTa677rq8/vWv71B1DAcBGQCgYfbs2ZkwYcKAtk996lM56aSTssMOOyRJdt11106UxjARkAEAfoef/exn+d73vpcDDzwwL33pS/PTn/600yUxhARkAID+/qS7Oxk3rrU9//wBhx999NHcddddueyyy/KP//iPOeWUU1Jr7USlDAMBGQAY2/r7k3nzkuXLk1pb25NPTu65Z12XKVOm5KijjkopJS984Qszbty43HHHHR0smqEkIAMAY9v8+cmaNQPbHnwwue22dbuvec1rcumllyZpLbd45JFHsttuuw1nlQwjARkAGNtWrBiwe0ySFyW54ZFHMmXKlCxcuDBz587NL37xi8ycOTNveMMbctJJJ6WU0pFyGXrjO10AAEBHTZ3aWlbRds7aB9OmJcuWrWv/4he/uO7xkiVLhqU0OsMMMgAwtvX2Jl1dA9u6ulrtjEkCMgAwtvX0JH19rRnjUlrbvr5WO2OSJRYAAD09AjHrmEEGAIAGARkAABoEZAAAaBCQAQCgQUAGAIAGARkAABoEZAAAaBCQAQCgQUAGAIAGARkAABoEZAAAaNjkgFxKOauUcnsp5dpG2wdKKTeXUpa2/7yycezkUspNpZQbSikvH+zCAQBgKGzODPLnkrxiA+2n11pntf9clCSllOcleUOSGe3n/HMpZbutLRYAAIbaJgfkWut3k9y5id1fnWRRrfWhWusvk9yU5IVbUB8AAAyrwViD/PZSyjXtJRi7ttv2SLKy0WdVuw0A2MbNnTs3EydOzMyZM9e1feADH8gee+yRWbNmZdasWbnooos6WCE8uVJr3fTOpXQnubDWOrO9PynJHUlqkv+XZHKtdW4p5cwkP6y1frHdb2GSi2qt527gnPOSzEuSSZMm7b9o0aKt+4mG0erVq7Pjjjt2ugw6xPiPbcZ/bDP+T+7qq6/O0572tHz4wx/OZz/72STJ5z73uTztaU/L61//+g5XNzi8B0a/Qw455Mpa6wEbOjZ+a05ca71t7eNSymeSXNjeXZXk2Y2uU5LcspFz9CXpS5IDDjigHnzwwVtT0rBasmRJRlO9DC7jP7YZ/7HN+D+5gw8+OMuWLcvHP/7xda/TkiVLsuOOO24zr5v3wLZtq5ZYlFImN3aPTLL2DhcXJHlDKWWHUsrvJ9kzyeVbcy0AYHT75Cc/mX333Tdz587NXXfd1elyYKM25zZv5yT5YZK9SimrSinHJfloKeXHpZRrkhyS5F1JUmu9LsmXklyfZHGSv661Pjbo1QMAo8Lb3va2/PznP8/SpUszefLkvPvd7+50SbBRm7zEotZ6zAaaFz5J/94kvVtSFAAwyvT3J/PnJytWJFOnJu9854DDkyZNWvf4+OOPzxFHHDHMBcKm8016AMDW6e9P5s1Lli9Pam1tTz45ueeedV1uvfXWdY+/+tWvDrjDBYw0W/UhPQCAzJ+frFmzbveYJEsefDB3rFqVKVOm5JRTTsmSJUuydOnSlFLS3d2df/mXf+lcvfA7CMgAwNZZsWLA7jlrH5SSrFqVJDnuuOOGtybYCpZYAABbZ+rUzWuHEU5ABgC2Tm9v0tU1sK2rq9UOo5CADABsnZ6epK8vmTattaxi2rTWfk9PpyuDLWINMgCw9Xp6BGK2GWaQAQCgQUAGAIAGARkAABoEZAAAaBCQAQCgQUAGAIAGARkAABoEZAAAaBCQAQCgQUAGAIAGARkAABoEZAAAaBCQAQCgQUAGAIAGARkAABoEZAAAaBCQAQCgQUAGAIAGARkAABoEZAAAaBCQAQCgQUAGAIAGARkAABoEZAAAaBCQAQCgQUAGAIAGARkAABoEZAAAaBCQAQCgQUAGAIAGARkAABoEZAAAaBCQAQCgQUAGAIAGARkAABoEZAAAaBCQAQCgQUAGgK00d+7cTJw4MTNnznzCsY997GMppeSOO+7oQGXAlhCQAWArHXvssVm8ePET2leuXJlvfvObmTp1ageqAraUgAwAW2n27NmZMGHCE9rf9a535aMf/WhKKR2oCthSAjIADIELLrgge+yxR57//Od3uhRgMwnIALAl+vuT7u5k3LjW9vzz1x1as2ZNent788EPfrBT1QFbQUAGgM11553JvHnJ8uVJra3tyScn99yTJPn5z3+eX/7yl3n+85+f7u7urFq1Kvvtt19+9atfdbhwYFOM73QBADDq3HxzsmbNwLYHH0xuuy1Jss8+++T2229fd6i7uztXXHFFdtttt+GsEthCZpABYHM9/PCA3WOSvCjJDY88kilTpmThwoUdKQsYHGaQAWBzbb/9gN1z1j6YNi1ZtuwJ3ZdtoA0YucwgA8Dm2mOPpKtrYFtXV9Lb25l6gEElIAPA5powIenra80Yl9La9vUlPT2drgwYBJZYAMCW6OkRiGEbZQYZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaNjkgl1LOKqXcXkq5ttE2oZTyzVLKje3tro1jJ5dSbiql3FBKeflgFw4AAENhc2aQP5fkFeu1nZTkklrrnkkuae+nlPK8JG9IMqP9nH8upWy31dUCAMAQ2+SAXGv9bpI712t+dZKz24/PTvKaRvuiWutDtdZfJrkpyQu3rlQAABh647fy+ZNqrbcmSa311lLKxHb7Hkkua/Rb1W57glLKvCTzkmTSpElZsmTJVpY0fFavXj2q6mVwGf+xzfiPbcYf74Ft29YG5I0pG2irG+pYa+1L0pckBxxwQD344IOHqKTBt2TJkoymehlcxn9sM/5jm/HHe2DbtrV3sbitlDI5Sdrb29vtq5I8u9FvSpJbtvJaAAAw5LY2IF+Q5M3tx29O8rVG+xtKKTuUUn4/yZ5JLt/KawEAwJDb5CUWpZRzkhycZLdSyqok709yapIvlVKOS7IiyeuSpNZ6XSnlS0muT/Jokr+utT42yLUDAMCg2+SAXGs9ZiOHDt1I/94kvVtSFAAAdIpv0gMAgAYBGQAAGgRkAABoEJABAKBBQAYAgAYBGQAAGgRkAABoEJABAKBBQAYAgAYBGQAAGgRkAABoEJABAKBBQAYAgAYBGQAAGgRkAABoEJABAKBBQAYAgAYBGQAAGgRkAABoEJABAKBBQAYAgAYBGQAAGgRkAABoEJABAKBBQAYAgAYBGQAAGgRkAABoEJABAKBBQAYAgAYBGQAAGgRkAABoEJABAKBBQAYAgAYBGQAAGgRkAABoEJABAKBBQAYAgAYBGQAAGgRkAABoEJABAKBBQAYAgAYBGQAAGgRkAABoEJABAKBBQAa22Ny5czNx4sTMnDlzXdvf//3fZ999982sWbNy2GGH5ZZbbulghQCw+QRkYIsde+yxWbx48YC2E088Mddcc02WLl2aI444Ih/84Ac7VB0AbBkBGdhis2fPzoQJEwa07bzzzuse33///SmlDHdZALBVxne6AGDbM3/+/Hz+85/PLrvskm9/+9udLgcANosZZGDz9Pcn3d3JuHGt7fnnP6FLb29vVq5cmZ6ennzyk58c7goBYKsIyMCm6+9P5s1Lli9Pam1tTz45ueeeDXZ/4xvfmHPPPXeYiwSArSMgA5tu/vxkzZqBbQ8+mNx227rdG2+8cd3jCy64IHvvvfdwVQcAg0JABjbdihUDdo9J8qIkNzzySKZMmZKFCxfmpJNOysyZM7Pvvvvm4osvzoIFCwa9jA3dXu7EE0/M3nvvnX333TdHHnlk7r777kG/LgBjg4AMbLqpUwfsnpPk1iSPTJuWVatW5bjjjsu5556ba6+9Ntdcc03+4z/+I3vsscegl7Gh28vNmTNn3XWf+9zn5sMf/vCgXxeAsUFABjZdb2/S1TWwraur1T6MNnR7ucMOOyzjx7duzHPQQQdl1apVw1oTANsOARnYdD09SV9fMm1aUkpr29fXah9BzjrrrBx++OGdLgOAUcp9kIHN09Mz4gJxU29vb8aPH5+eEVwjACObGWRg5Fv/3sv9/RvsdvbZZ+fCCy9Mf3+/b/ADYIuZQQZGtrX3Xl57e7nly1v76617Xrx4cT7ykY/kO9/5TrrWXycNAJvBDDIwsm3g3svHrFmTF73nPbnhhhvW3V7u7W9/e+67777MmTMns2bNyl/+5V92qGAARjszyMDItt69l5PW7eXy+OOtP23HHXfc8NUEwDbNDDIwsq137+Xf2Q4AW0lABka2EXLvZQDGDgEZGNlGyb2XAdh2WIMMjHwj/N7LAGxbzCADAECDgAwAAA0CMgAANAjIAADQICADAECDgAwAAA0CMgAANAjIAADQICADAECDgAwAAA0CMgAANAjIAADQICADAECDgAwAAA0CMgAANAjIAADQICADAECDgAwAAA0CMgAANAjIAADQICADAECDgAwAAA0CMgAANAjIAADQICADAECDgAwAAA0CMgAANAjIAADQICADAECDgAwAAA0CMgAANAjIAADQICADAECDgAwAAA0CMgAANAjIAADQICADAECDgAwAAA0CMgAANAjIAADQICADAECDgAwAAA0CMgAANAjIAADQMH4wTlJKWZbkviSPJXm01npAKWVCkn9P0p1kWZI/q7XeNRjXAwCAoTKYM8iH1Fpn1VoPaO+flOSSWuueSS5p7wMAwIg2lEssXp3k7Pbjs5O8ZgivBQAAg2KwAnJNcnEp5cpSyrx226Ra661J0t5OHKRrAQDAkCm11q0/SSm711pvKaVMTPLNJH+T5IJa6zMafe6qte66gefOSzIvSSZNmrT/okWLtrqe4bJ69ersuOOOnS6DDjH+Y5vxH9uMP94Do98hhxxyZWNp8ACD8iG9Wust7e3tpZSvJnlhkttKKZNrrbeWUiYnuX0jz+1L0pckBxxwQD344IMHo6RhsWTJkoymehlco338H3zwwcyePTsPPfRQHn300bz2ta/NKaec0umyRo3RPv5sHeOP98C2bauXWJRSnl5K2Wnt4ySHJbk2yQVJ3tzu9uYkX9vaawGDZ4cddsill16aq6++OkuXLs3ixYtz2WWXdbosAOi4wZhBnpTkq6WUtef7t1rr4lLK/yT5UinluCQrkrxuEK4FDJJSyrpfDz7yyCN55JFH0v7vGADGtK0OyLXWXyR5/gbaf5Pk0K09PzB0Hnvssey///656aab8td//dc58MADO10SAHScb9KDsaa/P+nuTsaNy3bPeU6WnnhiVq1alcsvvzzXXnttp6sDgI4TkGEs6e9P5s1Lli9Pam1t583LM77+9Rx88MFZvHhxpysEgI4TkGEsmT8/WbMmSfLrJHcnyZo1eeDkk/Otb30re++9dweLA4CRQUCGsWTFinUPb01ySJJ9k/zRypWZM2dOjjjiiE5VBgAjxqDcBxkYJaZObS2rSCsYX7W2fdq05P/+305VBQAjihlkGEt6e5OuroFtXV2tdgAgiYAMY0tPT9LX15oxLqW17etrtQMASSyxgLGnp0cgBoAnYQYZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFAZpOtXLkyhxxySKZPn54ZM2ZkwYIFnS4JAGDQje90AYwe48ePz2mnnZb99tsv9913X/bff//MmTMnz3ve8zpdGgDAoDGDzCabPHly9ttvvyTJTjvtlOnTp+fmm2/ucFUAAINLQGaLLFu2LFdddVUOPPDATpcCADCoBGSeXH9/0t2djBvX2vb3Z/Xq1Tn66KNzxhlnZOedd+50hQAAg0pAZuP6+5N585Lly5Nak+XL88jxx+foF784PT09OeqoozpdIQDAoBOQ2bj585M1a9bt1iTHPfBApi9blhNOOKFzdQEADCEBmY1bsWLA7veTfCHJpffdl1mzZmXWrFm56KKLOlIaAMBQcZs3Nm7q1NbyiraXpDWLnGnTkqVLO1QUAMDQMoPMxvX2Jl1dA9u6ulrtAADbKAGZjevpSfr6WjPGpbS2fX2tdgCAbZQlFjy5nh6BGAAYU8wgAwBAg4AMAAANAjIAADQIyAAA0CAgAwBAg4AMAAANAjIAADQIyAAA0CAgAwBAg4AMAAANAjIAADQIyAAA0CAgAwBAg4AMAAANAjIAADQIyAAA0CAgAwBAg4AMAAANAjIAADQIyAAA0CAgAwBAg4AMAAANAjIAADQIyAAA0CAgAwBAg4AMAAANAjIAADQIyAAA0CAgAwBAg4AMAAANAjIAADQIyAAA0CAgAwBAg4AMAAANAjIAADQIyAAA0CAgAwBAg4AMAAANAjIAADQIyAAA0CAgAwBAg4AMAAANAjIAADQIyAAA0CAgAwBAg4AMAAANAjIAADQIyAAA0CAgAwBAg4AMAAANAjIAADQIyAAA0CAgAwBAg4A8yBYsWJCZM2dmxowZOeOMMzbYZ+7cuZk4cWJmzpy5ru3LX/5yZsyYkXHjxuWKK64YpmoBAFifgDyIrr322nzmM5/J5ZdfnquvvjoXXnhhbrzxxif0O/bYY7N48eIBbTNnzsx5552X2bNnD1e5AABsgIA8iH7yk5/koIMOSldXV8aPH5+XvvSl+epXv/qEfrNnz86ECRMGtE2fPj177bXXcJUKAMBGCMiDaObMmfnud7+b3/zmN1mzZk0uuuiirFy5stNlAQCwGcZ3uoBtyfTp0/O+970vc+bMyY477pjnP//5GT/eSwwAMJqYQd5a/f1Jd3cyblzS3Z3jnvrU/OhHP8p3v/vdTJgwIXvuuecT+qS/v8NFAwCwMQLy1rjzzmTevGT58qTWZPny3H788Ul/f1asWJHzzjsvx4wf/4Q+mTcvOf/8TlcPAMAG+P3/1rj55mTNmgFNRz/wQH7zlrfkKdOn58wzz8yuxx33hD7HrFmTJe95T+4oJVOmTMkpp5ySCRMm5G/+5m/y61//Oq961asya9asfOMb3xjOnwYAgAjIW+fhh5/Q9L0kefTR5OqrWw0rVjyhzzlJ8vjjrT8NRx555KCXCADA5hnyJRallFeUUm4opdxUSjlpqK83rLbffsPtU6du+PHG+gAAMGIMaUAupWyX5Mwkhyd5XpJjSinPG8prbo3Fixdnr732yh/+4R/m1FNP/d1P2GOPpKtrYFtXV9Lb+9v93t7f3QcAgBFjqGeQX5jkplrrL2qtDydZlOTVQ3zNLfLYY4/lda97Xe64447ssMMOOeecc3L99dc/+ZMmTEj6+pJp05JSWtu+vqSn57d9enp+dx8AAEaMUmsdupOX8tokr6i1vrW9/+dJDqy1vr3RZ16SeUkyadKk/RctWjRk9TyZ6667Lp/4xCdywgkn5MMf/nD++I//OEnS8yRBdvXq1dlxxx2Hq0RGGOM/thn/sc344z0w+h1yyCFX1loP2NCxof6QXtlA24BEXmvtS9KXJAcccEA9+OCDh7ikDbvjjjvykpe8JIcddlg+/vGPZ/bs2fnv//7vPFk9S5YsedLjbNuM/9hm/Mc244/3wLZtqJdYrEry7Mb+lCS3DPE1t8iGZtJL2VC+BwBgWzbUM8j/k2TPUsrvJ7k5yRuSvHGIr7l5/uqvkr6+THnssaxMkgcfTJKsWrUqu+++e0dLAwBg+A1pQK61PlpKeXuSbyTZLslZtdbrhvKam+Wv/ir51KeSJH+U5MYkK/v7U3fdNYsWLcq//du/dbQ8AACG35B/UUit9aIkFw31dbZIX9+6h+OTfDLJXyS5+a678v53vzszZszoVGUAAHTIkH9RyIj22GMDdr+Q5MG0PkX4qU99KgsXLuxEVQAAdNDY/qrp7bYbEJLPabavWtWRkgAA6KyxPYM8b97mtQMAsM0b2wH5n/85edvbWjPGSWv7tre12pPcfffdee1rX5u9994706dPzw9/+MMOFgsAwHAY20ssklYYbgfi9b3jHe/IK17xinzlK1/Jww8/nDVr1gxzcYPjscceywEHHJA99tgjF154YafLAQAY0cb2DPKTuPfee/Pd7343xx13XJJk++23zzOe8YzOFrWFFixYkOnTp3e6DACAUUFA3ohf/OIXedaznpW3vOUtecELXpC3vvWtuf/++ztd1mZbtWpVvv71r+etb31rp0sBABgVBOSm/v6kuzsZNy6PHn54fnTllXnb296Wq666Kk9/+tNz6qmndrrCzfbOd74zH/3oRzNunKEGANgUUtNa/f2tu1csX57Umim/+lWm1JoDb7opSfLa1742P/rRjzpc5CZqB/0LS8nEb3wj+//0p52uCABg1BCQ15o/P2l8CO/3kjy71tzw3vcmSS655JI873nP61Bxm6ER9L+f5ILVq9P953+eN7zmNbn00kvzpje9qdMVAgCMaALyWitWPKHpE0l6brkl++67b5YuXZq/+7u/2+BT586dm4kTJ2bmzJnr2u68887MmTMne+65Z+bMmZO77rprqCofqBH0P5xkVZJltWbRU5+al73sZfniF784PHUAAIxSAvJaU6c+oWlWkiumTcs111yT888/P7vuuusGn3rsscdm8eLFA9pOPfXUHHroobnxxhtz6KGHDt/65Q0E/STJbbcNz/UBAEY5AXmt3t6kq2tgW1dXq/13mD17diZMmDCg7Wtf+1re/OY354YbbsgXvvCFLFiwILNmzcrOO++cM844YxALX88Ggn6SHDxtmnsgAwBsAgF5rZ6epK8vmTYtKaW17etrtW+B2267LZMnT85ee+2V6667Lk972tNy5ZVXpqurK0ceeeQgF9+wFUEfAAABeaCenmTZsuTxx1vbJwvH/f3Jj3+cjBvXujXc+ef/ztNfcsklec5znpNp06YNUsEbMMhBHwBgrBGQt8TaO0U8/HBSa+vWcCefnNxzz7oukyZNyq233pokufXWWzNx4sQsWrQoxxxzzNDXtzlBHwCAAQTkLbHeLeGSJA8++NsPwvX3509/9aucvfvuSXd3zn7nO3PEEUfkggsuyOte97rhrxcAgE0mIG+J9e4UcUySFyW54ZFHMmXChCycOzcn3Xtvvplkz+XL882vfCUveOCBPPOZz8zLXvayzJw5M8ccc0wefPDBTlQPAMCTEJC3xHp3ijgnya1JHpk2Lat23jnHPfxwnpnkkiQ3Jrnk8cfzpc9+NnfeeWeuuOKKXHvttXnssceyaNGi4a8dAIAnJSBviSe7U8QG7kO8Jsn3HnwwO+ywQx544IE8+uijWbNmTXbffffhqRcAgE0mIG+JtXeK2H77J94pYgP3Ie5Kcte0aTnxxBMzderUTJ48ObvssksOO+yw4a8dAIAnJSBvqZ6eZJ99nniniA3NLie5695787W+vvzyl7/MLbfckvvvv9/XPgMAjEAC8mBbO7v8zGcOaP7WXXfl92+4Ic+6+OI85SlPyVFHHZUf/OAHHSoSAICNEZCHQk9PsuOOA5qmJrnsscey5uSTU2vNJZdckunTp3emPgAANkpAHirrfVjvwCSvTbLfypXZZ5998vjjj2fevHkdKQ0AgI0b3+kCtllTp7a+Ya/hlCSnTJuWXHttZ2oCAOB3MoM8VJ7sVnAAAIxYAvJQWfthvWnTnngrOAAARixLLIZST49ADAAwyphBBgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgY3+kCgA3r7u7OTjvtlO222y7jx4/PFVdc0emSAGBMEJBhBPv2t7+d3XbbrdNlAMCYYokFAAA0CMgwQpVScthhh2X//fdPX19fp8sBgDFDQIaRor8/6e5Oxo1Lurvz/Xe/Oz/60Y/yn//5nznzzDPz3e9+t9MVAsCYICDDSNDfn8yblyxfntSaLF+e3d/3vqS/PxMnTsyRRx6Zyy+/vNNVAsCYICDDSDB/frJmzbrd+5Pct2ZNMn9+7r///lx88cWZOXNm5+oDgDHEXSxgJFixYsDubUmOTJLly/PoC1+YN77xjXnFK17RicoAYMwRkGEkmDq1tbyi7Q+SXJ0k06Yl113XqaoAYEyyxAJGgt7epKtrYFtXV6sdABhWAjKMBD09SV9fa8a4lNa2r6/VDgAMK0ssYKTo6RGIAWAEMIMMAAANAjIAADQIyAAA0CAgAwBAg4AMAAANAjKMQnPnzs3EiRMHfP300qVLc9BBB2XWrFk54IADcvnll3ewQgAYvQRkGIWOPfbYLF68eEDbe9/73rz//e/P0qVL88EPfjDvfe97O1QdAIxuAjKMQrNnz86ECRMGtJVScu+99yZJ7rnnnuy+++6dKA0ARj1fFALbiDPOOCMvf/nL8573vCePP/54fvCDH3S6JAAYlcwgwzbiU5/6VE4//fSsXLkyp59+eo477rhOlwQAo5KADKNFf3/S3Z2MG9fann/+gMNnn312jjrqqCTJ6173Oh/SA4AtJCDDaNDfn8yblyxfntTa2p58cnLPPeu67L777vnOd76TJLn00kuz5557dqpaABjVrEGG0WD+/GTNmnW7xyRZ8uCDuWPVqkyZMiWnnHJKPvOZz+Qd73hHHn300Tz1qU9NX19f5+oFgFFMQIbRYMWKAbvnrH1QSrJq1br2K6+8cvhqAoBtlCUWMBpMnbp57QDAFhOQYTTo7U26uga2dXW12gGAQSUgw2jQ05P09SXTprWWVUyb1trv6el0ZQCwzbEGGUaLnh6BGACGgRlkAABoEJABAKBBQAYAgAYBGQAAGgRkAABoEJABAKBBQAYAgAYBGQAAGgRkAABoEJABAKBBQAYAgAYBGQAAGgRkAABoEJABAKBBQIYxaO7cuZk4cWJmzpy5ru3qq6/Oi170ouyzzz75kz/5k9x7770drBAAOkdAhjHo2GOPzeLFiwe0vfWtb82pp56aH//4xznyyCPzj//4jx2qDgA6S0CGMWj27NmZMGHCgLYbbrghs2fPTpLMmTMn5557bidKA4COE5CBJMnMmTNzwQUXJEm+/OUvZ+XKlR2uCAA6Q0CGsaK/P+nuTsaNa23PP3/A4bPOOitnnnlm9t9//9x3333ZfvvtO1ElAHTc+E4XAAyD/v5k3rxkzZrW/vLlycknJ7vttq7L3nvvnYsvvjhJ8rOf/Sxf//rXO1EpAHScGWQYC+bP/204XuvBB5Pbblu3e/vttydJHn/88fzDP/xD/vIv/3I4KwSAEUNAhrFgxYoBu8ckeVGSGx55JFOmTMnChQtzzjnn5LnPfW723nvv7L777nnLW97SkVIBoNMssYCxYOrU1rKKtnPWPpg2LVm2bF37O97xjmEtCwBGIjPIMBb09iZdXQPburpa7QDAAAIyjAU9PUlfX2vGuJTWtq+v1Q4ADGCJBYwVPT0CMQBsAjPIAADQICADAECDgAwAAA0CMgAANAjIAADQICADAECDgAwAAA0CMgAANGxVQC6lfKCUcnMpZWn7zysbx04updxUSrmhlPLyrS8VAACG3mDMIJ9ea53V/nNRkpRSnpfkDUlmJHlFkn8upWw3CNcCthErV67MIYcckunTp2fGjBlZsGBBkuTOO+/MnDlzsueee2bOnDm56667OlwpAGPNUC2xeHWSRbXWh2qtv0xyU5IXDtG1gFFo/PjxOe200/KTn/wkl112Wc4888xcf/31OfXUU3PooYfmxhtvzKGHHppTTz2106UCMMaMH4RzvL2U8hdJrkjy7lrrXUn2SHJZo8+qdtsTlFLmJZmXJJMmTcqSJUsGoaThsXr16lFVL4PL+A+Ota/hs571rFx00UVZtGhRTj/99CxZsiR77rln3vWud+Xwww/vbJEbYPzHNuOP98C2rdRan7xDKd9K8nsbODQ/rRB8R5Ka5P8lmVxrnVtKOTPJD2utX2yfY2GSi2qt5z7ZtQ444IB6xRVXbP5P0SFLlizJwQcf3Oky6BDjP3iWLVuW2bNn59prr83UqVNz9913rzu26667jshlFsZ/bDP+eA+MfqWUK2utB2zo2O9cYlFr/eNa68wN/PlarfW2WutjtdbHk3wmv11GsSrJsxunmZLklq39QYBRrr8/6e5Oxo1rbfv7s3r16hx99NE544wzsvPOO3e6QgDY6rtYTG7sHpnk2vbjC5K8oZSyQynl95PsmeTyrbkWMMr19yfz5iXLlye1JsuX55Hjj8/RL35xenp6ctRRRyVpLbW69dZbkyS33nprJk6c2MmqARiDtvZDeh8tpfy4lHJNkkOSvCtJaq3XJflSkuuTLE7y17XWx7byWsBoNn9+smbNut2a5LgHHsj0ZctywgknrGv/0z/905x99tlJkrPPPjuvfvWrh7tSAMa4rfqQXq31z5/kWG+S3q05P7ANWbFiwO73k3whyT733ZdZs2YlST70oQ/lpJNOyp/92Z9l4cKFmTp1ar785S8Pe6kAjG2DcRcLgN9t6tTW8oq2l6Q1i5xp05KlSwd0veSSS4azMgAYwFdNA8Ojtzfp6hrY1tXVageAEURABoZHT0/S19eaMS6lte3ra7UDwAhiiQUwfHp6BGIARjwzyAAA0CAgwxi3cuXKHHLIIZk+fXpmzJiRBQsWJEm+/OUvZ8aMGRk3blxG0zdcAsDWssQCxrjx48fntNNOy3777Zf77rsv+++/f+bMmZOZM2fmvPPOy//+3/+70yUCwLASkGGMmzx5ciZPbn0p5k477ZTp06fn5ptvzpw5czpcGQB0hiUWwDrLli3LVVddlQMPPLDTpQBAxwjIMBb19yfd3cm4ca1tf39Wr16do48+OmeccUZ23nnnTlcIAB0jIMNY09+fzJvX+la7WpPly/PI8cfn6Be/OD09PTnqqKM6XSEAdJSADGPN/PnJmjXrdmuS4x54INOXLcsJJ5zQuboAYIQQkGGsWbFiwO73k3whyaX33ZdZs2Zl1qxZueiii/LVr341U6ZMyQ9/+MO86lWvystf/vKOlAsAw81dLGCsmTq1tbyi7SVpzSJn2rRk6dIBXY888sjhrAwARgQzyDDW9PYmXV0D27q6Wu0AgIAMY05PT9LX15oxLqW17etrtQMAlljAmNTTIxADwEaYQQYAgAYBGQAAGgRkAABoEJABAKBBQAYAgAYBGQAAGgRkAABoEJABAKBBQAYAgAYBGQAAGgRkgFFg5cqVOeSQQzJ9+vTMmDEjCxYsSJL8/d//ffbdd9/MmjUrhx12WG655ZYOVwow+gnIAKPA+PHjc9ppp+UnP/lJLrvsspx55pm5/vrrc+KJJ+aaa67J0qVLc8QRR+SDH/xgp0sFGPUEZIBRYPLkydlvv/2SJDvttFOmT5+em2++OTvvvPO6Pvfff39KKZ0qEWCbMb7TBQCweZYtW5arrroqBx54YJJk/vz5+fznP59ddtkl3/72tztcHcDoZwYZYKTq70+6u5Nx41rb/v6sXr06Rx99dM4444x1s8e9vb1ZuXJlenp68slPfrKjJQNsCwRkgJGovz+ZNy9ZvjypNVm+PI8cf3yOfvGL09PTk6OOOuoJT3njG9+Yc889twPFAmxbBGSAkWj+/GTNmnW7NclxDzyQ6cuW5YQTTljXfuONN657fMEFF2TvvfcezioBtknWIAOMRCtWDNj9fpIvJNnnvvsya9asJMmHPvShLFy4MDfccEPGjRuXadOm5dOf/vSwlwqwrRGQAUaiqVNbyyvaXpLWLHKmTUuWLl3X/spXvnK4KwPY5lliATAS9fYmXV0D27q6Wu0ADCkBGWAk6ulJ+vpaM8altLZ9fa12AIaUJRYAI1VPj0AM0AFmkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZALbSypUrc8ghh2T69OmZMWNGFixYkCQ58cQTs/fee2fffffNkUcembvvvruzhQKbREAGgK00fvz4nHbaafnJT36Syy67LGeeeWauv/76zJkzJ9dee22uueaaPPe5z82HP/zhTpcKbAIBGQC20uTJk7PffvslSXbaaadMnz49N998cw477LCMH9/6Tq6DDjooq1at6mSZwCYSkAFgEC1btixXXXVVDjzwwAHtZ511Vg4//PAOVQVsDgEZALZEf3/S3Z2MG9fa9vdn9erVOfroo3PGGWdk5513Xte1t7c348ePT4+vDodRQUAGgM11553JvHnJ8uVJrcny5Xnk+ONz9ItfnJ6enhx11FHrup599tm58MIL09/fn1JKB4sGNtX4ThcAAKPOzTcna9as261JjnvggUxftiwnnHDCuvbFixfnIx/5SL7zne+kq6urA4UCW0JABoDN9fDDA3a/n+QLSfa5777MmjUrSfKhD30of/u3f5uHHnooc+bMSdL6oN6nP/3p4a0V2GwCMgBsru23H7D7krRmkTNtWrJ06br2V77ylcNZFTBIrEEGgM21xx7J+ksmurqS3t7O1AMMKgEZADbXhAlJX19rxriU1ravL3GXCtgmWGIBAFuip0cghm2UGWQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgQkAEAoEFABgCABgEZAAAaBGQAAGgotdZO17BOKeXXSZZ3uo7NsFuSOzpdBB1j/Mc24z+2GX+8B0a/abXWZ23owIgKyKNNKeWKWusBna6DzjD+Y5vxH9uMP94D2zZLLAAAoEFABgCABgF56/R1ugA6yviPbcZ/bDP+eA9sw6xBBgCABjPIAADQICADAECDgLwFSimvKKXcUEq5qZRyUqfrYWiVUp5dSvl2KeUnpZTrSinvaLdPKKV8s5RyY3u7a6drZeiUUrYrpVxVSrmwvW/8x5BSyjNKKV8ppfy0/XfBi7wHxo5Syrvaf/9fW0o5p5TyVOO/bROQN1MpZbskZyY5PMnzkhxTSnleZ6tiiD2a5N211ulJDkry1+0xPynJJbXWPZNc0t5n2/WOJD9p7Bv/sWVBksW11r2TPD+t94L3wBhQStkjyd8mOaDWOjPJdkneEOO/TROQN98Lk9xUa/1FrfXhJIuSvLrDNTGEaq231lp/1H58X1r/Y9wjrXE/u93t7CSv6UiBDLlSypQkr0ryr41m4z9GlFJ2TjI7ycIkqbU+XGu9O94DY8n4JE8rpYxP0pXklhj/bZqAvPn2SLKysb+q3cYYUErpTvKCJP+dZFKt9dakFaKTTOxgaQytM5K8N8njjTbjP3b8QZJfJ/lse5nNv5ZSnh7vgTGh1npzko8lWZHk1iT31FovjvHfpgnIm69soM298saAUsqOSc5N8s5a672drofhUUo5IsnttdYrO10LHTM+yX5JPlVrfUGS++PX6WNGe23xq5P8fpLdkzy9lPKmzlbFUBOQN9+qJM9u7E9J61ctbMNKKU9JKxz311rPazffVkqZ3D4+OcntnaqPIfX/JfnTUsqytJZUvayU8sUY/7FkVZJVtdb/bu9/Ja3A7D0wNvxxkl/WWn9da30kyXlJXhzjv00TkDff/yTZs5Ty+6WU7dNaqH9Bh2tiCJVSSlprD39Sa/2nxqELkry5/fjNSb423LUx9GqtJ9dap9Rau9P67/3SWuubYvzHjFrrr5KsLKXs1W46NMn18R4YK1YkOaiU0tX+/8GhaX0Wxfhvw3yT3hYopbwyrTWJ2yU5q9ba29mKGEqllJck+V6SH+e3a1D/Lq11yF9KMjWtv0BfV2u9syNFMixKKQcneU+t9YhSyjNj/MeMUsqstD6kuX2SXyR5S1qTTN4DY0Ap5ZQkr0/rrkZXJXlrkh1j/LdZAjIAADRYYgEAAA0CMgAANAjIAADQICADAECDgAwAAA0CMsAIV0p5WinlO6WU7Ybg3N9qf1MYAG0CMsDINzfJebXWx4bg3F9I8ldDcF6AUct9kAE6pJTy/5LcUWtd0N7vTXJbrfXj6/X7QZI31lqXtb+s5JQktyWZldbX3v44yTuSPC3Ja2qtPy+lfC7JA0n2TjItrS+2eHOSFyX571rrse1z75rke7XWmUP5swKMJmaQATpnYdpfVVtKGZfWV1n3Nzu0v9L+D2qtyxrNz08rEO+T5M+TPLfW+sK0vuntbxr9dk3ysiTvSvIfSU5PMiPJPu1vhkut9a4kO7S/GRCACMgAHdMOvb8ppbwgyWFJrqq1/ma9brsluXu9tv+ptd5aa30oyc+TXNxu/3GS7ka//6itXxP+OK2Z6R/XWh9Pct16/W5PsvtW/0AA24jxnS4AYIz71yTHJvm9JGdt4PgDSZ66XttDjcePN/Yfz8C/1x/aQJ8N9Xtq+zoAxAwyQKd9NckrkvxRkm+sf7C9BGK7Usr6IXlQlFJKWuF82VCcH2A0EpABOqjW+nCSbyf50pPcpeLiJC8ZohL2T3JZrfXRITo/wKjjLhYAHdT+cN6Pkryu1nrjRvq8IMkJtdY/H4LrL0hyQa31ksE+N8BoZQYZoENKKc9LclOSSzYWjpOk1npVkm8PxReFJLlWOAYYyAwyAAA0mEEGAIAGARkAABoEZAAAaBCQAQCgQUAGAICG/x8dq7WiPm2MiAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x864 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Differences plot\n",
    "diff = pos - pos_i[:,0:2]\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(12,12))\n",
    "\n",
    "ax.scatter(diff[:,0], diff[:,1], color='red')\n",
    "\n",
    "for i in range(0,diff.shape[0]):\n",
    "    ax.text(diff[i,0], diff[i,1],f\"{i}\")\n",
    "\n",
    "ax.grid(True)\n",
    "ax.set_title(\"Differences from initial position\")\n",
    "ax.set_xlabel('x (mm)')\n",
    "ax.set_xlabel('y (mm)')\n",
    "plt.savefig('differences.png')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "A 0     -2.058633\n",
       "A 1     -3.592050\n",
       "A 2     -1.047261\n",
       "A 3      2.249717\n",
       "A 4     -3.463294\n",
       "A 5     -9.235294\n",
       "A 6     -0.448300\n",
       "A 7     -2.747641\n",
       "A 8     -9.765133\n",
       "A 9     -1.660864\n",
       "A 10    -7.264131\n",
       "A 11     1.744161\n",
       "A 12     2.319761\n",
       "A 13    -3.815255\n",
       "A 14    -6.187632\n",
       "A 15    -2.155132\n",
       "A 16    -5.363327\n",
       "A 17    -1.683721\n",
       "A 18    -9.420326\n",
       "A 19    -3.287352\n",
       "A 20    -0.972542\n",
       "A 21     2.658438\n",
       "A 22   -10.961463\n",
       "A 23    -0.707905\n",
       "Name: 0, dtype: float64"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Residual plot\n",
    "rres = radius_residual(res.x)\n",
    "rres"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEWCAYAAAB/tMx4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAXOklEQVR4nO3dfbRddX3n8fenRBEIDyKaRrAGLYoKPhEVH8rcFFupUKGzZBWLFiyW6RQVXTA2Yn3q1Jauqi1dVkfUjkxFIyKtVJYWBpuhOsqYCDbE4ICCAmLwgYChiEa/88feGU5u7s09N7k3J7+c92utrHP2w9n7+zvn5nP2+e2zfydVhSSpPb8w6gIkSdvHAJekRhngktQoA1ySGmWAS1KjDHBJapQBrm1KsjbJxKjrGKUkv5XktiQbkzxjHvdzXpIPbmP5rUleOAf7OT3J53d0Oxo9A3yMTRUIk/9zV9VTqmrlDNtZkqSSLJinUkftncCrq2phVV03Xzupqj+rqlfN1/a1+zHAtcvbBd4YHgusHWbFXaBWjREDXNs0eJSe5NlJViW5N8n6JO/uV7umv93QdzM8N8kvJPnjJN9KcleS/5Fk/4Ht/m6/7AdJ3jxpP29LcmmSjyS5Fzi93/cXk2xIcmeS9yR56MD2KskfJrkpyY+S/Nckj+8fc2+SSwbXn9TGKWtNsmeSjcAewFeTfGOax1eSs5LcBNzUzzshyfV9vf87yVMH1v+jJHf0dX49ybED7f7IwHqvGHiO3jRpnx9O8qcD0xNJbh+YXp7kG/0+vpbkt6apPUn+qm/3PUn+LckRU62rXY8Brtm4ALigqvYDHg9c0s8/pr89oO9m+CJwev9vGfA4YCHwHoAkTwbeC5wKLAb2Bw6etK8TgUuBA4CLgZ8BrwcOAp4LHAv84aTHHAccBRwNvAG4sN/HY4AjgJdN064pa62qB6pqYb/O06rq8dM+M3AS8BzgyUmeCfwd8J+ARwDvBy7v3xCeCLwaeFZV7Qu8CLh18sb65+h9wCuAR/fbOWQb+5/sG8Cv0D23bwc+kmTxFOv9Ot3r9wS65/q3gR/MYj8aIQNc/9gfJW5IsoEuWKfzU+CXkxxUVRur6kvbWPdU4N1V9c2q2gi8ETil72J4KfBPVfX5qvoJ8BZg8qA8X6yqf6yqn1fV/VW1uqq+VFWbqupWulD8D5Me8xdVdW9VrQVuAK7s938P8BlguhOQ26p1WH9eVT+sqvuB3wfeX1XXVtXPquoi4AG6N5afAXvSBf1DqurWqprqyP6lwKer6pqqegB4M/DzYYupqk9U1Xf65+/jdJ8Mnj3Fqj8F9gUOB1JV66rqzlm0WyNkgOukqjpg8z+2PqoddAbdkdqNSb6c5IRtrPto4FsD098CFgCL+mW3bV5QVf/O1kd9tw1OJHlCkk8n+W7frfJndEfjg9YP3L9/iumFTG1btQ5rsN7HAudMemN8DPDoqroZeB3wNuCuJCuSPHqamgafo/uYxZFx30V1/cD+j2Dr54uq+hzdJ6O/BdYnuTDJfsPuR6NlgGtoVXVTVb0MeBTwF8ClSfZh66NngO/QBdlmvwRsogvVOxnoDkiyF10XwRa7mzT9PuBG4LC+C+c8INvfmqFrHdZgvbcB7xh8Y6yqvavqYwBV9dGqekG/z6J7Lie7ky70AUiyN1s+R/cBew9M/+LAuo8FPkDXVfOI/o35BqZ5vqrqb6rqKOApdG/Q/2W4JmvUDHANLcnLkzyyqn4ObOhn/wz4Ht3H+8cNrP4x4PVJDk2ykO6I+eNVtYmub/s3kzyvP7H4dmYO432Be4GNSQ4H/vNctWuGWrfHB4A/SPKc/iThPkmOT7Jvkicm+dUkewI/pvtk8LMptnEpcEKSF/TP0Z+w5f/X64EXJzkwyS/SHdVvtvlN9XsASV5JdwS+lSTP6ut8CN2bwo+nqUe7IANcs3EcsLb/ZsYFwClV9eO+C+QdwBf6j+xH053E+3u6b6jcQhcMrwHo+6hfA6ygO9L8EXAXXT/xdM4Ffqdf9wPAx+ewXdPWuj2qahVdP/h7gLuBm+lOkkLX/30+8H3gu3SfZs6bYhtrgbOAj9I9R3cDtw+s8vfAV+lOgF7JwPNRVV8D3gV8ke5TxJHAF6Ypdz+65/Nuuq6jH9B9710NiD/ooFHrj3o30HWP3DLicqRmeASukUjym0n27vvQ3wmsYYqv00mangGuUTmR7uThd4DD6Lpj/DgozYJdKJLUKI/AJalRO3XgnYMOOqiWLFmyxbz77ruPffbZZ2eWscuw7ePZdhjv9tv22bd99erV36+qR06ev1MDfMmSJaxatWqLeStXrmRiYmJnlrHLsO0Toy5jZMa5/bZ9YtaPS/KtqebbhSJJjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY3aqVdiqrNk+RUAfPi48bycWNLc8AhckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJatRQAZ7k9UnWJrkhyceSPCzJgUmuSnJTf/vw+S5WkvSgGQM8ycHAa4GlVXUEsAdwCrAcuLqqDgOu7qclSTvJsF0oC4C9kiwA9ga+A5wIXNQvvwg4ac6rkyRNK1U180rJ2cA7gPuBK6vq1CQbquqAgXXurqqtulGSnAmcCbBo0aKjVqxYscXyjRs3snDhwh1qRGvW3HEPAIfuv8fYtX2zcXzdB41z+2377Nu+bNmy1VW1dPL8BTM9sO/bPhE4FNgAfCLJy4fdcVVdCFwIsHTp0pqYmNhi+cqVK5k8b3d3+vIrAPjwcfuMXds3G8fXfdA4t9+2T8zZ9obpQnkhcEtVfa+qfgpcBjwPWJ9kMUB/e9ecVSVJmtEwAf5t4OgkeycJcCywDrgcOK1f5zTgU/NToiRpKjN2oVTVtUkuBb4CbAKuo+sSWQhckuQMupA/eT4LlSRtacYAB6iqtwJvnTT7AbqjcUnSCHglpiQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYCP0Jo77mHJ8itY0v/AgyTNhgEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRC0ZdwDjxp9MkzSWPwCWpUQa4JDXKAJekRhngktQoA1ySGjVUgCc5IMmlSW5Msi7Jc5McmOSqJDf1tw+f72IlSQ8a9gj8AuCzVXU48DRgHbAcuLqqDgOu7qclSTvJjAGeZD/gGOBDAFX1k6raAJwIXNSvdhFw0vyUKEmaSqpq2yskTwcuBL5Gd/S9GjgbuKOqDhhY7+6q2qobJcmZwJkAixYtOmrFihVbLN+4cSMLFy7coUa0Ys0d92wxvWgvWH9/d//Ig/cfQUWjM06v+1TGuf22ffZtX7Zs2eqqWjp5/jABvhT4EvD8qro2yQXAvcBrhgnwQUuXLq1Vq1ZtMW/lypVMTEwM246mTb4S85wjN/GuNd3FsLeef/woShqZcXrdpzLO7bftE7N+XJIpA3yYPvDbgdur6tp++lLgmcD6JIv7jS8G7pp1VZKk7TZjgFfVd4Hbkjyxn3UsXXfK5cBp/bzTgE/NS4WSpCkNO5jVa4CLkzwU+CbwSrrwvyTJGcC3gZPnp0RJ0lSGCvCquh7Yqv+F7mhckjQCXokpSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjFoy6gHGwZPkVoy5B0m7II3BJapQBLkmNMsAlqVEGuCQ1ygCXpEYNHeBJ9khyXZJP99MHJrkqyU397cPnr0xJ0mSzOQI/G1g3ML0cuLqqDgOu7qclSTvJUAGe5BDgeOCDA7NPBC7q718EnDSnlUmStilVNfNKyaXAnwP7AudW1QlJNlTVAQPr3F1VW3WjJDkTOBNg0aJFR61YsWKL5Rs3bmThwoU71Ihd3Zo77ply/qK9YP393f0jD95/J1Y0euPwum/LOLffts++7cuWLVtdVUsnz5/xSswkJwB3VdXqJBOz3XFVXQhcCLB06dKamNhyEytXrmTyvN3N6dNciXnOkZt415ruJbj11ImdWNHojcPrvi3j3H7bPjFn2xvmUvrnAy9J8mLgYcB+ST4CrE+yuKruTLIYuGvOqpIkzWjGPvCqemNVHVJVS4BTgM9V1cuBy4HT+tVOAz41b1VKkrayI98DPx/4tSQ3Ab/WT0uSdpJZjUZYVSuBlf39HwDHzn1JkqRheCWmJDXKAJekRhngktQoA3wXs2T5Ff6Cj6ShGOCS1CgDXJIaZYBLUqMMcElq1Kwu5NH88cSlpNnyCFySGmWAS1KjDHBJapQBvovygh5JMzHAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1asGoC9idzcVPom3exq3nH7/D25K0e/EIXJIaZYBLUqMMcElqlAEuSY0ywCWpUTMGeJLHJPmXJOuSrE1ydj//wCRXJbmpv334/JcrSdpsmCPwTcA5VfUk4GjgrCRPBpYDV1fVYcDV/bQkaSeZMcCr6s6q+kp//0fAOuBg4ETgon61i4CT5qlGSdIUUlXDr5wsAa4BjgC+XVUHDCy7u6q26kZJciZwJsCiRYuOWrFixRbLN27cyMKFC7en9l3emjvu2ebyRXvB+vuH29aRB+8/BxXtOnbn130Y49x+2z77ti9btmx1VS2dPH/oAE+yEPhfwDuq6rIkG4YJ8EFLly6tVatWbTFv5cqVTExMDFVDa2a6EvOcIzfxrjXDXQy7u12JuTu/7sMY5/bb9olZPy7JlAE+1LdQkjwE+CRwcVVd1s9en2Rxv3wxcNesq5IkbbdhvoUS4EPAuqp698Ciy4HT+vunAZ+a+/IkSdMZ5vP784FXAGuSXN/POw84H7gkyRnAt4GT56XCBs3FIFaSNJMZA7yqPg9kmsXHzm05kqRheSWmJDXKAJekRhngktQoA7wRS5Zf4clRSVswwCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRM/4qvXZtg7/Sc+v5x4+wEkk7m0fgktQoA1ySGmWAS1KjDHBJapQnMaex+eTgTCcGPYkoaVQ8ApekRhngktQoA1ySGmWAS1Kjxuok5rAnJrf12FHbVh070j5J7fEIXJIaZYBLUqMMcElqVDN94JP7d6fqCx6273dbfcWTtzubvu9dpZ98JvaVSzvHfF/o5xG4JDVqhwI8yXFJvp7k5iTL56ooSdLMtjvAk+wB/C3wG8CTgZclefJcFSZJ2rYdOQJ/NnBzVX2zqn4CrABOnJuyJEkzSVVt3wOTlwLHVdWr+ulXAM+pqldPWu9M4Mx+8onA1ydt6iDg+9tVRPts+/ga5/bb9tl7bFU9cvLMHfkWSqaYt9W7QVVdCFw47UaSVVW1dAfqaJZtH8+2w3i337bPXdt3pAvlduAxA9OHAN/ZsXIkScPakQD/MnBYkkOTPBQ4Bbh8bsqSJM1ku7tQqmpTklcD/wzsAfxdVa3djk1N270yBmz7+Brn9tv2ObLdJzElSaPllZiS1CgDXJIaNZIAT/KXSW5M8m9J/iHJAQPL3thfmv/1JC8aRX07wzgNQ5DkMUn+Jcm6JGuTnN3PPzDJVUlu6m8fPupa50uSPZJcl+TT/fRYtD3JAUku7f+/r0vy3DFq++v7v/cbknwsycPmuu2jOgK/Cjiiqp4K/F/gjQD9pfinAE8BjgPe21+yv1sZw2EINgHnVNWTgKOBs/r2LgeurqrDgKv76d3V2cC6gelxafsFwGer6nDgaXTPwW7f9iQHA68FllbVEXRf9DiFOW77SAK8qq6sqk395JfovkMO3aX4K6rqgaq6BbiZ7pL93c1YDUNQVXdW1Vf6+z+i+098MF2bL+pXuwg4aSQFzrMkhwDHAx8cmL3btz3JfsAxwIcAquonVbWBMWh7bwGwV5IFwN5018nMadt3hT7w3wM+098/GLhtYNnt/bzdzbi0cytJlgDPAK4FFlXVndCFPPCoEZY2n/4aeAPw84F549D2xwHfA/573330wST7MAZtr6o7gHcC3wbuBO6pqiuZ47bPW4An+Z9938/kfycOrPMmuo/XF2+eNcWmdsfvOY5LO7eQZCHwSeB1VXXvqOvZGZKcANxVVatHXcsILACeCbyvqp4B3Mdu2F0ylb5v+0TgUODRwD5JXj7X+5m3X+Spqhdua3mS04ATgGPrwS+jj8vl+ePSzv8vyUPowvviqrqsn70+yeKqujPJYuCu0VU4b54PvCTJi4GHAfsl+Qjj0fbbgdur6tp++lK6AB+Htr8QuKWqvgeQ5DLgecxx20f1LZTjgD8CXlJV/z6w6HLglCR7JjkUOAz4P6OocZ6N1TAESULXD7quqt49sOhy4LT+/mnAp3Z2bfOtqt5YVYdU1RK61/lzVfVyxqPt3wVuS/LEftaxwNcYg7bTdZ0cnWTv/u//WLpzP3Pa9pFciZnkZmBP4Af9rC9V1R/0y95E1y++ie6j9mem3krb+iOyv+bBYQjeMdqK5k+SFwD/CqzhwX7g8+j6wS8BfonuD/7kqvrhSIrcCZJMAOdW1QlJHsEYtD3J0+lO3j4U+CbwSroDx3Fo+9uB36bLsuuAVwELmcO2eym9JDVqV/gWiiRpOxjgktQoA1ySGmWAS1KjDHBJapQBrpFJ8rMk1/dX6P7T4KiUs9jG0iR/M82yW5MctJ21vS3JudMse12S392e7c6wzxP6r55JQzHANUr3V9XT+9HafgicNdsNVNWqqnrt3Jc2tX5got8DPjoPm7+C7qrNvedh29oNGeDaVXyRfkCvJI9P8tkkq5P8a5LD+/kn90frX01yTT9vYmCM7UckubIfOOn99GPOJFmS5IbNO0pybpK39fd/P8mX+21+cojw/FXgK5tH00yyMslfJbmmH+/6WUku68d7/tOB/d/YD+Z0Q5KLk7wwyRf69Z4N0A8psZJuiAlpRga4Rq4fH/1YHhxO4ELgNVV1FHAu8N5+/luAF1XV04CXTLGptwKf7wdOupzuareZXFZVz+q3uQ44Y4b1nw9MHpjqJ1V1DPDf6C6NPgs4Aji9v+IS4JfpxsZ+KnA48DvAC/r2nTewrVXArwxRtzR/g1lJQ9gryfXAErpQvKofsfB5wCe6ISSAbtgFgC8AH05yCXAZWzsG+I8AVXVFkruHqOGI/kj5ALrLnP95hvUXs+UPM8CDbzxrgLWbhwtN8k26Qcs20A1stKafv5ZuUP9Ksoau/ZvdRTd6nTQjj8A1SvdX1dOBx9KNlXEW3d/khr5vfPO/JwH04+X8MV0oXj9wdDtoqrEhNrHl3/rDBu5/GHh1VR0JvH3SsilrnmKdB/rbnw/c3zy9YNI6k9cbXGdzbffPUIMEGODaBVTVPXQ/P3UuXXjdkuRk6EYyTPK0/v7jq+raqnoL8H22HJIX4Brg1H7d3wA2/97geuBRfR/5nmzZx7wvcGc/3O2pQ5S7jq47ZL48AbhhxrUkDHDtIqrqOuCrdEOungqckeSrwFoe/Lm5v0yypj8heU2//qC3A8ck+Qrw63SjvVFVPwX+hG70w08DNw485s39/KsmzZ/OZ+i6aubLMrpvo0gzcjRCaZaS/APwhqq6aY63uwj4aFUdO5fb1e7LAJdmqf+BgkVVdc0cb/dZwE+r6vq53K52Xwa4JDXKPnBJapQBLkmNMsAlqVEGuCQ1ygCXpEb9PxTgXbZCQA4DAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# ij resdual\n",
    "ijres = m_ij_residual(res.x)\n",
    "plt.hist(ijres, bins='fd');\n",
    "plt.title(\"Histogram of residuals\")\n",
    "plt.xlabel(\"Residual (mm)\")\n",
    "plt.grid(True)\n",
    "plt.savefig(\"residual_histogram.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The API expects 3D coordinates (with a z value which is zero in this case). Therefore we add a column of zeros."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "95th percentile of residuals is 6.60 mm\n",
      " res[2,22] = -10.5\n",
      " res[3,22] = -11.2\n",
      " res[4,22] = -7.5\n",
      " res[5,22] = -7.7\n",
      " res[8,13] = 10.0\n",
      " res[12,22] = -11.0\n",
      " res[13,8] = 10.0\n",
      " res[13,18] = -7.1\n",
      " res[13,19] = -8.2\n",
      " res[13,20] = -7.7\n",
      " res[13,21] = -16.7\n",
      " res[13,22] = 79.5\n",
      " res[14,22] = -14.9\n",
      " res[15,22] = -14.5\n",
      " res[18,13] = -7.1\n",
      " res[19,13] = -8.2\n",
      " res[20,13] = -7.7\n",
      " res[21,13] = -16.7\n",
      " res[22,2] = -10.5\n",
      " res[22,3] = -11.2\n",
      " res[22,4] = -7.5\n",
      " res[22,5] = -7.7\n",
      " res[22,12] = -11.0\n",
      " res[22,13] = 79.5\n",
      " res[22,14] = -14.9\n",
      " res[22,15] = -14.5\n"
     ]
    }
   ],
   "source": [
    "absres = np.abs(ijres)\n",
    "p95 = np.percentile(absres, 95)\n",
    "\n",
    "print(f\"95th percentile of residuals is {p95:4.2f} mm\")\n",
    "ijres_matrix = np.zeros_like(m_ij)\n",
    "for r, i, j in zip(ijres, non_nan_indices[:,0], non_nan_indices[:,1]):\n",
    "    if np.abs(r) > p95:\n",
    "        print(f\" res[{i},{j}] = {r:3.1f}\")\n",
    "    ijres_matrix[i,j] = r\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Output best estimates\n",
    "\n",
    "Write some json files..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-2.000e-03,  2.080e-01,  0.000e+00],\n",
       "       [-2.000e-03,  4.060e-01,  0.000e+00],\n",
       "       [ 3.460e-01, -1.650e-01,  0.000e+00],\n",
       "       [ 5.790e-01, -2.720e-01,  0.000e+00],\n",
       "       [-1.760e-01, -1.120e-01,  0.000e+00],\n",
       "       [-6.270e-01, -4.000e-01,  0.000e+00],\n",
       "       [-1.000e-03,  1.261e+00,  0.000e+00],\n",
       "       [ 1.000e-03,  1.433e+00,  0.000e+00],\n",
       "       [ 2.000e-03,  1.697e+00,  0.000e+00],\n",
       "       [ 1.000e-03,  1.862e+00,  0.000e+00],\n",
       "       [ 0.000e+00,  2.192e+00,  0.000e+00],\n",
       "       [ 2.000e-03,  2.388e+00,  0.000e+00],\n",
       "       [ 8.510e-01, -3.990e-01,  0.000e+00],\n",
       "       [ 9.730e-01, -4.690e-01,  0.000e+00],\n",
       "       [ 1.463e+00, -6.840e-01,  0.000e+00],\n",
       "       [ 1.740e+00, -8.160e-01,  0.000e+00],\n",
       "       [ 1.888e+00, -8.860e-01,  0.000e+00],\n",
       "       [ 2.177e+00, -1.025e+00,  0.000e+00],\n",
       "       [-7.630e-01, -4.860e-01,  0.000e+00],\n",
       "       [-1.057e+00, -6.690e-01,  0.000e+00],\n",
       "       [-1.332e+00, -8.390e-01,  0.000e+00],\n",
       "       [-1.438e+00, -9.040e-01,  0.000e+00],\n",
       "       [-1.548e+00, -1.003e+00,  0.000e+00],\n",
       "       [-1.756e+00, -1.106e+00,  0.000e+00]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result = np.zeros((n_ant, 3))\n",
    "result[:,:-1] = np.round(pos/1000.0, 3)\n",
    "result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{\n",
      "    \"antenna_positions\": [\n",
      "        [\n",
      "            -0.002,\n",
      "            0.208,\n",
      "            0.0\n",
      "        ],\n",
      "        [\n",
      "            -0.002,\n",
      "            0.406,\n",
      "            0.0\n",
      "        ],\n",
      "        [\n",
      "            0.346,\n",
      "            -0.165,\n",
      "            0.0\n",
      "        ],\n",
      "        [\n",
      "            0.579,\n",
      "            -0.272,\n",
      "            0.0\n",
      "        ],\n",
      "        [\n",
      "            -0.176,\n",
      "            -0.112,\n",
      "            0.0\n",
      "        ],\n",
      "        [\n",
      "            -0.627,\n",
      "            -0.4,\n",
      "            0.0\n",
      "        ],\n",
      "        [\n",
      "            -0.001,\n",
      "            1.261,\n",
      "            0.0\n",
      "        ],\n",
      "        [\n",
      "            0.001,\n",
      "            1.433,\n",
      "            0.0\n",
      "        ],\n",
      "        [\n",
      "            0.002,\n",
      "            1.697,\n",
      "            0.0\n",
      "        ],\n",
      "        [\n",
      "            0.001,\n",
      "            1.862,\n",
      "            0.0\n",
      "        ],\n",
      "        [\n",
      "            0.0,\n",
      "            2.192,\n",
      "            0.0\n",
      "        ],\n",
      "        [\n",
      "            0.002,\n",
      "            2.388,\n",
      "            0.0\n",
      "        ],\n",
      "        [\n",
      "            0.851,\n",
      "            -0.399,\n",
      "            0.0\n",
      "        ],\n",
      "        [\n",
      "            0.973,\n",
      "            -0.469,\n",
      "            0.0\n",
      "        ],\n",
      "        [\n",
      "            1.463,\n",
      "            -0.684,\n",
      "            0.0\n",
      "        ],\n",
      "        [\n",
      "            1.74,\n",
      "            -0.816,\n",
      "            0.0\n",
      "        ],\n",
      "        [\n",
      "            1.888,\n",
      "            -0.886,\n",
      "            0.0\n",
      "        ],\n",
      "        [\n",
      "            2.177,\n",
      "            -1.025,\n",
      "            0.0\n",
      "        ],\n",
      "        [\n",
      "            -0.763,\n",
      "            -0.486,\n",
      "            0.0\n",
      "        ],\n",
      "        [\n",
      "            -1.057,\n",
      "            -0.669,\n",
      "            0.0\n",
      "        ],\n",
      "        [\n",
      "            -1.332,\n",
      "            -0.839,\n",
      "            0.0\n",
      "        ],\n",
      "        [\n",
      "            -1.438,\n",
      "            -0.904,\n",
      "            0.0\n",
      "        ],\n",
      "        [\n",
      "            -1.548,\n",
      "            -1.003,\n",
      "            0.0\n",
      "        ],\n",
      "        [\n",
      "            -1.756,\n",
      "            -1.106,\n",
      "            0.0\n",
      "        ]\n",
      "    ]\n",
      "}\n"
     ]
    }
   ],
   "source": [
    "json_result = {}\n",
    "json_result[\"antenna_positions\"] = result.tolist()\n",
    "print(json.dumps(json_result, indent=4, separators=(',', ': ')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The position file above is in a format suitable for uploading to the telescope api. Also for using as the calibrated_antenna_positions.json file in TART/software/containers/telescope_web_api/config_data.\n",
    "\n",
    "\n",
    "    tart_upload_antenna_positions --api https://tart.elec.ac.nz/rhodes --pw=xxxx --file rhodes_antenna_positions.json\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}