{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical solution of the geodesics equation\n",
    "## Schwarzshild implementation\n",
    "\n",
    "In class, you have computed the geodesics equation in s Schwarzshild spacetime. \n",
    "After some algebraic manipulation, the equation for the orbit of a massive particle can casted into the form\n",
    "\n",
    "## $\\frac{dx^2}{d\\phi^2} - 1 + x = \\frac{3 G^2 M^2}{L^2} x^2$,\n",
    "\n",
    "where $x = \\frac{L^2}{G M r}$, $M$ is the mass parameter in the Schwarzshild metric, $L$ is the angular momentum per unit mass (a conserved quantity in s Schwarzshild spacetime) and $G$ is the gravitational constant. \n",
    "\n",
    "You have seen in class how to solve perturbatively this equation. In this tutorial we will solve the equation numerically. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before solving the orbit equation in a Schwarzshild space time, we will consider a more general problem, i.e. \n",
    "how to solve numerically a $1$st order differential equation in the form\n",
    "\n",
    "## $\\frac{dy}{dt} = f(y, y', t)$,\n",
    "\n",
    "where $f$ denotes a generic function, and $y' = \\frac{dy}{dt}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 1200x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Import some useful packages \n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Routine that solves a 1st order differential equation using the Eucler method - seen in the tutorial's sheet.\n",
    "def odeEuler(f,y0,t):\n",
    "    '''Approximate the solution of y'=f(y,t) by Euler's method.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    f : function\n",
    "        Right-hand side of the differential equation y'=f(t,y), y(t_0)=y_0\n",
    "    y0 : number\n",
    "        Initial value y(t0)=y0 wher t0 is the entry at index 0 in the array t\n",
    "    t : array\n",
    "        1D NumPy array of t values where we approximate y values. Time step\n",
    "        at each iteration is given by t[n+1] - t[n].\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    y : 1D NumPy array\n",
    "        Approximation y[n] of the solution y(t_n) computed by Euler's method.\n",
    "    '''\n",
    "    y = np.zeros(len(t))\n",
    "    y[0] = y0\n",
    "    for n in range(0,len(t)-1):\n",
    "        y[n+1] = y[n] + f(y[n],t[n])*(t[n+1] - t[n])\n",
    "    return y\n",
    "\n",
    "# Discretization of the variable t in Npoints in the range [tin, tfin]\n",
    "# Change Npoints to test the accuracy of the method\n",
    "Npoints = 10\n",
    "tin = 0.0\n",
    "tfin = 2.0\n",
    "t = np.linspace(tin, tfin,Npoints)\n",
    "\n",
    "# Set the initial condition\n",
    "y0 = 1\n",
    "\n",
    "#Define the function on the right-end side of the differential equation\n",
    "f = lambda y,t: y\n",
    "\n",
    "#Call my Eucler-solver\n",
    "y = odeEuler(f,y0,t)\n",
    "\n",
    "#Compute the analytical/exact solution \n",
    "y_true = np.exp(t)\n",
    "\n",
    "#Compare the two in a plot\n",
    "import matplotlib.gridspec as gridspec\n",
    "from matplotlib import ticker\n",
    "from matplotlib import rc\n",
    "\n",
    "plt.rcParams['mathtext.fontset'] = 'cm'\n",
    "import matplotlib.gridspec as gridspec\n",
    "from matplotlib import ticker\n",
    "from matplotlib import rc\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=(12, 8))\n",
    "plt.plot(t,  y_true, ls = '-', lw = '2', color='C0', label = 'analytical solution')\n",
    "plt.plot(t,  y, ls = '--', lw = '2', color='C1', label = 'numerical - Euler method')\n",
    "#plt.plot(t,y,'b.-',t,y_true,'r-')\n",
    "plt.legend(loc = 'upper left', fontsize = 28, frameon=True, ncol = 1)\n",
    "plt.axis([0,2,0,9])\n",
    "plt.grid(True)\n",
    "plt.title(\"Solution of $y'=y , y(0)=1$\", fontsize = 28)\n",
    "plt.tick_params(labelsize=22)\n",
    "plt.xlabel(r'$t$', fontsize = 32)\n",
    "plt.ylabel(r'$y(t)$', fontsize = 32)\n",
    "plt.show()\n",
    "plt.close(fig)\n",
    "# I could add the computation of the relative difference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Possible extensions\n",
    "You can change the parameter 'Npoints' and test how the precision of the Euler's method on the number of sample points."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# System of 2 differential equation, a Newtonian example\n",
    "In the tutorial sheet, you showed that a 2nd order ordinary differential equation can be written as a system of two 1st order ordinary differential equation. \n",
    "\n",
    "Now we consider the Newtonian geodesic equation for a systema Sun (source of gravitational field) - Mercury (test massive particle). We just set to zero the relativistic correction in the equation considered in class:\n",
    "\n",
    "## $\\frac{dx^2}{d\\phi^2}= 1 - x$.\n",
    "\n",
    "We apply the Eucler method to solve this equation. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x initial:  1.2079587324981578\n",
      "v initial:  0.0\n"
     ]
    }
   ],
   "source": [
    "# Routine that solves a 1st order differential equation using the Eucler method.\n",
    "def ode2Euler(f, v0, x0, phi):\n",
    "    '''Approximate the solution of the system\n",
    "       v'= f(x, phi)\n",
    "       x' = v\n",
    "       \n",
    "       by Euler's method.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    f : function\n",
    "        Right-hand side of the differential equation y'=f(t,y), y(t_0)=y_0\n",
    "    x0 : number\n",
    "        Initial value x(phi_in)=x0 where phi_in is the entry at index 0 in the array phi\n",
    "    v0 : number\n",
    "        Initial value v(phi_in)=v0 where phi_in is the entry at index 0 in the array phi\n",
    "    phi : array\n",
    "        1D NumPy array of phi values where we approximate x values. Angular step\n",
    "        at each iteration is given by phi[n+1] - phi[n].\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    x : 1D NumPy array\n",
    "        Approximation x[n] of the solution x(t_n) computed by Euler's method.\n",
    "    '''\n",
    "    x = np.zeros(len(phi))\n",
    "    v = np.zeros(len(phi))\n",
    "    x[0] = x0\n",
    "    v[0] = v0\n",
    "    for n in range(0,len(phi)-1):\n",
    "        v[n+1] = v[n] + f(x[n], phi[n])*(phi[n+1] - phi[n])\n",
    "        x[n+1] = x[n] + v[n]*(phi[n+1] - phi[n])\n",
    "    return x\n",
    "\n",
    "\n",
    "Npoints = 21000 #2100000 #number of points\n",
    "Norbits = 3.0  #number of orbits\n",
    "\n",
    "# Discretization of the variable t in Npoints in the range [phi_in, phi_fin]\n",
    "phi_in = 0.0\n",
    "phi_fin = 2.0*np.pi*Norbits\n",
    "phi = np.linspace(phi_in, phi_fin, Npoints)\n",
    "\n",
    "# Here I compute the numerical factors, lenghts are in unit of R0=1.0e10\n",
    "rS = 2.95e-7 #units of R0 \n",
    "L2 = 8.196e-7 #units of R0^2\n",
    "# Set the initial condition\n",
    "r0 = 4.6\n",
    "x0 = 2.0*L2/rS/r0\n",
    "print('x initial: ', x0)\n",
    "v0 = 0.0\n",
    "print('v initial: ', v0)\n",
    "\n",
    "#Define the function on the right-end side of the differential equation\n",
    "f = lambda x,ph: (1-x)\n",
    "\n",
    "#Call my Eucler-solver\n",
    "x = ode2Euler(f, v0, x0, phi)\n",
    "\n",
    "#Compute the radial coordinate from the solution x of the differential equation\n",
    "r = 2.0*(L2/rS)/x\n",
    "#r = r/np.sqrt(1.-rS/4.6)\n",
    "# Compute the x and y coordinate for a better visualisation of the orbit\n",
    "xpos = r*np.cos(phi)\n",
    "ypos = r*np.sin(phi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now we plot the Newtonian orbit "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAosAAAJwCAYAAADhvyvNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdeXxU9b3/8feHJEBYRDaDSyHIolZAqahYQcSNCi5oXapXW71aWmu1pb1aFBEoolRtrZb29mdvqy1VW+oCrazWBZeKikgAN1BZBCSyh4QQksn398ckITNzZrLNzJnl9Xw85pGcc2bOfKIkeee7mnNOAAAAgJdWfhcAAACA1EVYBAAAQFSERQAAAERFWAQAAEBUhEUAAABERVgEAABAVLl+F5CpunXr5goLC/0uAwAAoEHvvvvududcd69rhMUEKSws1LJly/wuAwAAoEFmtiHaNbqhAQAAEBVhMQozyzez283sHTPbbWb7zGydmf3DzE73uz4AAIBkoBvag5n1lrRYUl9JX0paIqlCUqGkiyUVSXrDr/oAAACShbAYxszaS3pBUh9J0yRNc85V1rveVVJXn8oDAABIKsJipLsUDIp/cc7dHX7RObdD0o6kVwUAAOADxizWY2atJX235nCGn7UAAACkAloWQ52kYBfz5865D83s65IuqDm3VdJC59ybfhYIAACQTITFUANrPq41s8clfSfs+t1m9oyka51z5UmtDAAAwAd0Q4fqUvPxDEnflvSggjOiOys4C3qzpG9K+q0v1QEAACQZYTFU7X+PXEl/dM7d5pz71Dm32zn3T0ljJTlJ3zGzo8NfbGbjzGyZmS3btm1bEssGAABIDMJiqL31Pv9D+EXn3DJJ7yr43+1Mj+uPOueGOOeGdO/uub0iAABAWiEshlpf7/N1UZ5Te75HYksBAADwHxNcQi2v93lXSV59yd1qPpYmspCKigqVlJSorKxMFRUVqq6uTuTbAVkrJydH7dq10yGHHKKOHTvKzPwuCQBSCmGxHufcZjN7S9Kpks6W9FH962bWWdLXag6XJaqOkpISbd26VZ06dVL37t3Vtm1btWrVil9iQJw55xQIBFRaWqrt27ervLxchx12GN9rAFAP3dCRptd8vNvMTqw9aWZtJf2vpE4KjltMyHqL5eXlKi4uVs+ePVVQUKD27dsrJyeHX15AApiZcnNzdeihh6pXr14qKyvT3r17G34hAGQRwmIY59y/FFwy5zBJb5nZq2b2nKRPJV2p4PI5VznnXCLef9euXerSpYvatm2biNsDiCInJ0ddunRRSUmJ36UAQEohLHpwzt0m6VJJbyi4UPdoSfsk/UrSYOfc2kS9d2lpqQ455JBE3R5ADB06dNC+ffv8LgMAUgpjFqNwzj0n6blkv28gEFBuLv9bAD/k5OQoEAj4XQYApBRaFlMQ4xMBf/C9BwCRCIsAAACIirAIAACAqBgcBwBJVFpRpeKS/dq9r1K1vd7b91Zoxee7lZfTSgWHtNW/ijbrzc92Num+Jul7I/ro6O7t1blda7XJbaUNO8rUr6CjOrTJVbcObdS5fZ5a57BmK4CmISwCQDM553QgUK0dpQf08sdfave+Sj2w6GN/apH0+yWfxuVeT373VB3WsY3a5OboyEPz1aoV4RLIZoRFIMtMmTJFU6dOlRQMO6mmttVr8uTJmjJliq+1BKqdtpbs1/rtZfr9kk/12trtvtaTLFf/4a0Gn3PWsYfpF98cpG4dWtNSCWQ4wiKArOecU2lFlTbtKtemXeU6f8I8v0tKeS999KVOnv5vz2sd2uRo6Z3nqEMbfsUAmYDvZABZpTJQrbXFe1VVnfhW1fUzxiT8PeKl38T5qgzE579JaUVAAyYvijhf2LWdXvjJCOXlMLcSSCeERQApJZ5d4/srA1pbvFfxjoWd8nNVNHlUnO/qr7XTRzfqeYUtaHVdv2Of+k1cEHF+zT3nq3UuARJIVYRFABnhQFVAm3aVq7SiKi73y8RAGA+xWkubGyT73xUaIPNyTGvuOZ+xkECKICwCSDvOOe07ENCn20pbfK+u7dvoyM75dccf7s1Pq+7jVBLtv1tTQ2RlwKn3HfNDzn1672jlMCsb8AXt/oCkHTt26Pbbb9cxxxyj/Px8FRQU6Nxzz9VzzwW3B3/88cdlZjIzrV+/3vMegUBAf/7zn3XhhRfqyCOPVJs2bdSlSxcNHTpU9957r0pKShqsY9++fXrwwQc1bNgwdevWTW3atNHhhx+uCy64QE8++WSjumg3b96sm2++WUcffbTatm2rI444QhdddJH+/W/vyQjR7N+/XzNnztQ555yjHj16qHXr1urevbtGjhypmTNnqqKiIubr33vvPX33u9/Vscceqw4dOqhNmzY68sgjNXjwYN14442aPXu25z1q/zvXnwntnFNFZUArN+3Wyk27tWrzHv37zeV6YMqdumLUcA0f0Fsn9+2h0aefqO9+62L95dHfauuWTRH3HnTUoRGP+kERibF+xpiIR1P1uXO+CifMq3uk4kx+IFMZ33CJMWTIELds2bImv+7DDz/Ucccdl4CKEM2qVat07rnnqri42PP6uHHjdNppp+n666+XJK1bt06FhYUhz1m3bp0uvvhirVq1Kur7HH744Zo7d65OPvnkqHWMGTNGn3/+edR7DB8+XHPnzlXnzp09r7/22mu64IILogbTKVOmyDnX4NI5y5cv1yWXXKKNGzdGreXYY4/V888/rz59+kRce+SRRzR+/HhVV1dHfb0U/Pd+7LHHhpyr7Xq8++67ddl3x0eMN6yurtbM++/R479/RIFAIOq9h54+XG++/mrM949WE9+D/mjReEhag4EWMbN3nXNDvK7RDY2stnv3bn3jG9+oC4rf+ta39O1vf1uHHXaYPvnkEz3yyCN69NFHVVRUFPUeW7du1emnn64vvvhCeXl5uv7663X22WersLBQFRUVeuWVV/TrX/9aX3zxhc4//3wtX75cPXv2DLnHli1bNHLkSO3YsUOSdPXVV+uaa66pq+M3v/mN3njjDb322msaM2aMXnvtNeXk5ITcY+PGjXVBsVWrVrrxxht1+eWX69BDD9WqVas0Y8YMTZkyRUOGeP4sqPPBBx9oxIgRKi0tVfv27XXTTTfptNNOU8+ePVVSUqKFCxdq5syZ+uijj/SNb3xDy5YtU6dOnepev3LlyrqgWFhYqJtvvlmDBw9W165dVVZWprVr12rJkiWaO3duxHuv215W9/mXeys8J6bcd9dtmj3rT5KkLt2668rv3KjBJw/VkP5f0a5du7R8+XI999xzap2X4/FqpLLwwNd7wrxGT06qHzT7HdZeL/zkzPgVBmQ75xyPBDxOOukk1xwffPBBs16H5vnxj3/sFNz8ws2YMSPielVVlbv44ovrniPJrVu3LuQ5F110kZPkjjjiCPfRRx95vs9nn33mCgoKnCR3zTXXRFy/4oor6u7/8MMPR1wPBALuyiuvrHvOI488EvGcyy67rO76Y489FnG9pKTEnXDCCSFfS7jq6uq65xx//PHuiy++8Px63nnnHdeuXTsnyd11110h1yZNmuQkufbt20d9vXPOlZWVuZLSMrdq025X9Pmuukdtbd8f/7OQ80Wf73Iz//z3uutDhgxx27dvj3r/jRs3Rr0WC9+DqavvnfNcr5893+QHgIZJWuaiZBq6oROEbujUV1FRoR49emj37t0aNGiQVqxY4Tn7sri4WIWFhdq/f7+k0G7o999/XwMGDJAkzZ49W5dffnnU9/vd736nm2++WXl5edq9e7fatWsnSfriiy/Us2dPVVVV6YwzztCSJUs8X79nzx4dffTR2rlzp/r27au1a9fWXdu6dauOOuooBQIBnXPOOXrhhRc87/H222/r1FNPrTsO//6fN2+eLrjggrrnRusyl6Tbb79dDzzwgI444ght3ry57vy4ceP0hz/8QYMHD9by5csjXldd7bR6y56o9z3hK8Eu9u+P/5lu+smEuvODjjpUw4YN0xtvvKE2bdro448/Vq9evaLep7n4Hkwfzem2prsa8BarG5oJLmmsdqB3fTc8/o4KJ8zTvz84OP7uybc2qnDCPN3x7Mq6c8Ul+1U4YZ5OCduB4YLfvKbCCfO0atPBX+YPvbBGhRPmqfeEeZo0Z7UkadWmPSqcME8X/Oa1kNefMv3fKpwwT8Ul++vO3fHsShVOmKcn3zo4/u3fHxSrcMI83fD4Ow1+TYmybNky7d69W5L0X//1X1GX6SgoKNCoUd5LqMyZM0eSlJeXp4svvjjm+51xxhmSpMrKStX/Q+Lll19WVVVwuZcbb7wx6us7deqkK664QpL0ySefhEy0efnll+vG79WOrfRyyimn6Pjjj496vfbr6dWrV8ygWP/r2bJlS8jYxiOOOEJSsDv77bfflnNO5fUmp8QKivWZhU5I2blzp/7zn/9Iki699NKEBEWkl+ZMmqk/SQZA4xAW0WhOCgl86W716tV1n5900kkxnxttnF9t6KusrFSbNm3qZvJ6PQYOHFj3uq1bt3rWMXTo0Jh11L9efzJN/c8bCnmnnHJK1Gu1X8+GDRtifi1mpgsvvNDz67nqqqvUunVrVVRU6PTTT9eIc7+he+5/SGs+XN3ghJdBRx1a9/lhHduGXFuxYkVdS+jw4cNj3gfZieAIJAYTXNKY1w/EP14XGRSuPrWnrj41dEJFwSFtPV///C2Rv4THn9tfO8sO6Mm3NtbdZ+BRnTxf//bEcyLO3XfpIN136aCQc+d8tcDz9cnsItq1a1fd5927d4/53GjXv/zyy2a99759++o+37lzZ93nhx12WMzX9ejRw/N1TblHQUFB1Gst/Xqqq50q2hfo/t89pim33aLdu3bqtRcX67UXF0uSDul0qIYOP1Njr/gvnT4y+G+lfkCMZdu2bXWfH3744c2qE9mj/s+SE6Yu0p7y2Iu11w+MdFUDoQiLaJRpYwdo2tgBfpeRMM3dKaK26/eoo47SggWR25hFc9RRRzWrjmhjjOufb+49pINfz6mnnqr/+7//i3mf+vcrb9tVKzftrjs3ctRonXL6cC1+fq7+s+RFLX/7TW3/slgle3Zr8fNztPj5ORo9erSefvrpRr1HOHb2QFPU34ln0pzVmrV0Q8znExyBUIRFZK36axU21KJWv1Wrvm7dutVd79+/v1q3bt3kOrp06VL3eXFxsQ455JCoz62/FmT91zXlHrG+1m7duqm4uFjFxcV1E3ei+XznPu3ad0CS1NbjevsOHXXJt67RJd+6RpKUX75Nzz//vGbOnKnPPvtM8+fP18SJE/WrX/0q5vvUr63Wli1bGvUaIFz9P3yb0uJ47dCemjZ2YMznApmKMYvIWvUnerz77rsxnxttZvvgwYMlBWdWv/zyy82qo34oe+utt2I+t/71+mMg63/+zjuhk4bCxbpe+/WsX79ea9asibheVV1dN1GlNijGUn+CSr9+/TR+/HgtW7asrit89uzZDd6jfm21LYqvvtr0xbaBcEWTR9WNcWyorXrW0o2Mb0TWIiwiaw0ZMqRuMeknnngiavdscXGxFi1a5Hlt7NixdZ//4he/aNYWZCNHjlRubrCR/09/+lPU55WUlNSFq759+4bsIjNy5Mi6Rbr//Oc/R73HsmXLQibUhKv/9dx33311n3++c59WbtqtD7Y0vGVh1/Zt6gKil86dO+trX/uaJGn79u0N3q9Wly5ddPrpp0uSnnvuOW3YELsrEWiKdU2YHFMbGvtNnN/gc4FMQFhE1mrbtq2+/e1vSwruOnL//fdHPKe6ulrf+9736tZYDHfSSSdp9OjRkoLL19x2220xA2NxcXHEWMDDDz9cl156ad09fv/730e8zjmnm266qW6Hl1tvvTXiHrVL9yxevFizZs2KuEdpaanGjRsXtTZJuuSSS+paOh9//HFN+PkvYrYibtq4QQvmPK2cVlYXEN9+ZWHI5KFwO3furGvJ7d27d8x6wk2YEFx3saKiQpdddlnIxJ5wsbZNBGJp7KzqyoCjtRFZgUW5E4RFudPDzp07dfzxx9ct/XLVVVeFbPf38MMP6z//+Y9OOeUUvf3225KCXbT11/grLi7WySefXBdOBg8erBtvvFEnnHCC2rVrp927d2v16tV68cUXtWDBAg0cODCiW3vLli0aNGiQduzYITPTtddeq6uvvlrdu3fXp59+qkceeUSvv/66JOm0007z3O5v/fr1GjRokPbu3atWrVpp3Lhxuvzyy9WpU6e67f4+/vhjDRkypO79vb7/P/jgAw097TTtrdlfeshpw3TBpVeqT/9jlJfXWrt37dSaD9/Xf5a8qLffeFWXXHJJyESVM888U2+//bZGjx6ts846S8cdd5w6d+6sPXv2qKioSDNnzqxbUPyRRx7RLbfcEvL+tV3NkydP1pQpUyLq+973vqdHH31UUnDm980336zhw4erc+fO2rVrl1asWKFnn31Wubm5zRoawPcgvDR260EmxCBdxVqUm7CYIITF9FFUVKRzzz036iSW6667TsOHD9cNN9wgKbimYPjyM5s3b9aVV16pN954o8H3GzlypF566aWI8ytXrtSYMWO0adOmqK8dNmyY5s6dGzKhpb5XXnlFF110kfbu3et5/e6775aZaerUqZIOhkXnnD7fVa7dNS2In639WD/9/nX6bM1HDX49119/fUj3+Zlnnhl1F5r6br31Vv3617+OmNncUFisrq7WbbfdpoceeihmK+6IESP0yiuvNFhHOL4H0ZDGtCTm5ZjWTh+dhGqA+GAHFyCGE044QR988IF++tOfql+/fmrTpo26deumkSNH6sknn9Rjjz2mkpKDY/VqxznWd+SRR+r111/X888/r2uuuUZ9+vRR+/btlZubq65du+qUU07RLbfcogULFkTdim/QoEH66KOP9MADD+j0009Xly5dlJeXpx49emj06NF64okn9Oqrr0YNilIwqL3//vu66aab1KtXL7Vu3VoFBQUaM2aMFi5cWBcS61u9eY9Wbd5TFxQl6eh+x+jpxa9rxm/+oHPHXKwjjvqK8vPzlZeXp4KCAg0bNky33XabXn311YhxlrNnz9YTTzyh66+/XoMHD9bhhx+uvLw8tWvXTscee6xuuOEGvfnmm3r44YebtQROq1at9Mtf/lIrVqzQD37wAx133HHq2LGj8vPzdfTRR+vss8/Www8/rL/97W9NvjfQGHRRI9vQspggtCxmlhtvvFF//OMfdeSRR8Zs+UsXzjmt2tzwtntmpoFHRobjTMb3IJqqMWs3SnRRI7XFallknUWgAeXl5Zo7d66khrfjS3VVgWp98EXjZjQf2Tk/CRUB6a/+2o2xWhJrrxEakW4Ii8h6n376qY4++mjPLtFAIKCbbrqpbomX2tnT6aa62mn1loZbEhu79R4Ab7VBkNCITEJYRNa777779Prrr+uqq67Sqaeeqh49eqi8vFxFRUX6wx/+oOXLl0sKTky58MILfa62aSoD1fqwES2JhEQgvmqDYKxZ1IRGpAvCIiDp448/9px5W2vo0KGaPXt22uxJ3NgxiYREILHW1QTBWFsLEhqR6giLyHq33367+vTpoxdeeEHr16/Xtm3bVFlZqa5du+qkk07SlVdeqauuukqtWqX+4gHVzmk1IRFIOUWTR0mSzntoidYUl3o+h9CIVMVs6ARhNjSSiZAYP3wPIhkaM4Oa0IhkYp1FIIOt3LQ7ZlA0s5h7NQNIvmljB2j9jDHqlB+9g491GpEq6IYG0pBzTu9vKVF1Az0DBEQgtdE9jXRAWATSzKZd+7Sz7EDM5xASgfSyePwISQ3Pnu6Un1sXMIFkoRsaSBOBaqeVm3bHDIp0NwPpbV0DWwnuKa+iaxpJR8sikOIas6A2O64AmaWhxb3pmkYyERZTkHMubdbzQ2K9v2WPAtXRxyW2zctR/4KOSawos7E6BFINoRGpgG7oFJOTk6OqKu+FW5E9DlRVa+Wm3VGDYu0MZ4JifAUCAeXk5PhdBhBh/YwxunZor6jXCyfM03kPLUliRcgmhMUU06FDB5WUNLw9GzLXyk279dHW6P8GBh11qAYe2SmJFWWP0tJStWvXzu8yAE+1y+3k5Xj3PK0pLmU8IxKCbugU07lzZ23atEnt27dX27Zt/S4HSdRQlzMTVxIrEAho586d6tatm9+lADGtnT5aUuyuadPBrQaBlqJlMcXk5+eroKBAGzduVHFxscrKyhQIBBhLlcGqArG7nHNaGUExQZxzqqqq0u7du7Vhwwa1b99eHTvStY/0sD7GzGmn6GESaCq2+0uQ5m73V6uiokIlJSUqKytTRUWFqqur41gdUsXmXeVR11STpKOY4ZxwOTk5ateunQ455BB17NiRyWVIS7HWZ5SYAIOGxdruj7CYIC0Ni8hs+w5U6at3L4p6/dqhvTRt7IAkVgQgE8RqTexf0KFu8W8gHGHRB4RFeHHOqe+d8xWI8m2Xl2N145EAoLlihUZaGeElVlhkzCKQJJPmrFLvO6IHxfUzxhAUAcRFrEDIWEY0FS2LCULLIurre+c8VUUZdsperwASiVZGNAYti4BPdpUdUOGE6EFx/YwxBEUACUUrI1qKsAgkyMDJizR42gue1zrl5/IXPYCkibUDTOGEeepNaEQMdEMnCN3Q2cs5p953zI96nZAIwE90S8ML3dBAkny0tSRqUMzLMX4QA/BdQ62M7DGNcLQsJggti9kn1l/rrJsIIBXRyohatCwCCVRd7Rr8gUtQBJCK1s8Yo7wc712LmPyCWrQsJggti9nh1qeW659FX3heozURQDqJFg7Z+SU7sIOLDwiLma/PHfNiLrANAOmGbunsRTc0EGeFEwiKADLP+hlj1L+gg+c1uqWzFy2LCULLYma645mVeuqdzz2v0VUDIJPQLZ1d6Ib2AWEx8wyb8ZI27S73vEZrIoBMRLd09qAbGmih3hPmERQBZB26pSERFoGYnAsui+PV/s6WfQCyweLxI6L+rCMwZgfCIhDF9tKKqLuxXDu0l4omj0pyRQDgn1iBcdKc1UmuBsnEmMUEYcxieos1kYXWRADZLFprYl6Oae300UmuBvHCmEWgCc568BWCIgBEsX7GGHXKz404XxmIvZsV0hdhEajhnFPfO+bps+1lEdfycoygCAA1iiaPYhxjFiEsApIC1U6975ivKo9RGf0LOtC1AgAeCIzZgbAISOpzZ/SJLCw+CwDRERgzH2ERWa0yUB31B9r6GWM0beyAJFcEAOmHwJjZCIvIWj/+2wr1m7jA8xrjEwGgadbPGCPzOF84YZ5OmLoo6fUgfgiLDTCze83M1Tz+x+96EB8TnlmpOSs2e14jKAJA86ybMUZ5OZGRcU95lfpN9B7ug9RHWIzBzE6WdLvkuYEH0tSgqYv0N4+lcZjxDAAtt3b6aF07tFfEeZbWSV+ExSjMrI2kxyUVS5rrbzWIl1ufek8l5VUR5zvl5zLjGQDiZNrYAYxjzCCExeh+Lumrkr4vaY/PtSAOvjppof5ZtCXiPFv3AUBiEBgzA2HRg5mdKumnkp50zv3L73rQcoOmLNK+ykDE+WuH9mLGMwAkEIEx/REWw5hZW0l/lrRT0o98LgdxUDhhnkr2R3Y9ExQBIDkIjOmNsBhpuqRjJN3inNvudzFomeMmRV8ah6AIAMlDYExfhMV6zOzrkn4saY5z7u9+14OWKZwwT+WV1RHnmfEMAP4gMKYnwmINM8uX9JikEkk/aOY9xpnZMjNbtm3btrjWh6aJtSsLAMA/BMb0Q1g86F5J/SX9xDn3RXNu4Jx71Dk3xDk3pHv37vGtDo12/N0LPc8TFAEgNRAY0wth8aBLJFVL+o6ZvVL/IekbNc+5qebc//lWJWIqnDBPZQciZz0TFAEgtRAY0wdhMVQrSSM8HgU114+uOR7iS3WIia5nAEgvBMb0QFis4ZwrdM6Z10PBpXQk6baacyf6WSsiDZrivUk9QREAUhuBMfURFpH2Bk5e5LmOIkERANIDgTG1ERaR1obNeEl7KwiKAJDuCIypi7CItHXWL1/Rpt3lEecJigCQngiMqYmw2AjOuetqxio+6HctCLprzip9tq0s4vy1Q3v5UA0AIF6iBcZ+E+cnuRLUIiwi7RyoqtZfl26MOM9ezwCQGbwCY2XA6byHlvhQDQiLSCsl+yvV/67I/Z77F3QgKAJABvEKjGuKSzVpzmofqsluhEWklUFTFkecu3ZoLy0eP8KHagAAieQVGGct3eBDJdmNsIi04TXAuVN+Li2KAJDBvAIjE16Si7CItOD1gyEvx1Q0eZQP1QAAkslr8iKBMXkIi0h5fe7w/oGwdvroJFcCAPDDtLED1Ck/N+I8gTE5CItIaQOnLFLARZ5nLUUAyC7RepIIjIlHWETKumvOKu1lGz8AQI1oP/9PmLooyZVkF8IiUlJ1tfNcS5GgCADZzev3wJ7yKpbUSSDCIlJOZaBaR98ZuVI/u7MAACSW1Ek2wiJSTr+JkYtuszsLAKA+ltRJHsIiUsqgKZHjTvJyjKAIAIjAkjrJQVhEyrj96SKVeExoYYkcAICXaWMHyDzOExjji7CIlLCz7IBmL9sUcZ4JLQCAWNbxeyLhCIvwnXNOX5v2QsR5giIAoDEYv5hYhEX4rvcdzHwGALQMgTFxCIvw1UCPCS2d8nOZ0AIAaDImvCQGYRG++cnsFZ47tETb0gkAgFiiNTSww0vLEBbhi8pAtZ5dvjniPOMUAQAtEW2HFzQfYRG+iLbwNgAALcX4xfgiLCLpRj74csQ5U/TuAwAAmorAGD+ERSTV7n0HtG77vojzrJMFAEiGSXNW+11C2iEsImmqq51O/DnrKQIAksPr98uspRt8qCS9ERaRNP3vYj1FAEBy0R3dcoRFJMWdz61UVXXoubwcY5wiACDhvBomWE6n8QiLSIon3/o84tza6aN9qAQAkG28GiZYTqfxCItIOK/mfsYpAgCSie7o5iMsIqF++OTyiHP9Czr4UAkAINsRGJuHsIiEen7lFxHnFo8f4UMlAABInfJz/S4h7RAWkTC96X4GAKSYosmjIs7RuhgbYREJMfwXL8mFneOvOQBAKqA7umkIi4i7/ZUBfTcRDScAACAASURBVL6rPOK8119zAAD4gQaMxiMsIu6OnbQw4hzdzwCAVEJ3dOMRFhFXdzyzMuIcu7QAAFIR3dGNQ1hEXD31TuTi2+zSAgBIVeZ3AWmAsIi46XMHs58BAOllHa2LDSIsIi72VwYUCJv+zOBhAEA68GrYYO/ogwiLiAuvSS3MfgYApCv2jj6IsIgWO+vBlyPO0f0MAEgnTHaJjrCIFtlfGdBn2/eFnGOwMAAgHdHQ4Y2wiBbx6n72GiwMAEA6onWRsIg461/Qwe8SAABoNq/WxUlzVvtQSeogLKLZvP7aWjx+hA+VAACQOLOWbvC7BF8RFtEstz9dFHGOsR4AgEzAZJdQhEU0mXNOs5dtCjmXl8O0FgBA5mCr2oMIi2iyr94dOall7fTRPlQCAEBieG1Vm62ti4RFNFl5ZXXIMa2KAIBMxPCqIMIimqTvnZF/VdGqCADIFtnYukhYRKPtO1ClqtBGRcZ0AAAyGkvpEBbRBF+9O3JTda8xHQAAZJLw4VbZtpQOYRGNUhWojjhHqyIAIBt4DbfKptZFwiIape/EBRHnaFUEAGSLTvm5IcfZ1LpIWESDKqoCEeeYIQYAyCZFk0dFnDvvoSU+VJJ8hEU06Ji7ItdVBAAg24QPv1pTXOpTJclFWERMP3xyecQ5WhUBANnIa/jVCVMjJ39mGsIiYnp+5Rd+lwAAQMoIb13cU17lUyXJQ1hEVHc+tyriHK2KAIBs5tW6mOljFwmLiOrJtzaGHLOtHwAA2Td2kbAIT598GfkPn239AADwbl3M5HUXCYvwdM6vQpvUaVUEAOCg8NbFTF53kbCIRqFVEQCAg7KpdZGwiAiFE+b5XQIAACkvW1oXCYsI4ZyLOMcMaAAAImXLtreERYToc+d8v0sAACBt9C/oEHKcib1zhEXUOVBVreqwhkVaFQEAiG7x+BF+l5BwhEXUOXn6C36XAABA2glfMSTTWhcJi6gTvmVR+MBdAAAQKdNXDCEsQpJ01i9fiTiXLQN3AQCIt0xqXSQsQpL02baykGNaFQEAaLxMHuNPWIT+5x9FEedoVQQAoGXOe2hJw09KA4TFeswsz8zONrNfmtlSM/vCzA6Y2WYze9rMzvS7xkR4+t1NIcfhywAAAICGhbcuriku9amS+Mr1u4AUM0JS7ZTgrZLelVQm6auSvinpm2Y2zTl3t0/1JUU2LAMAAAAah5bFUNWSnpF0hnPucOfcBc65K51zAyV9S1JA0iQzG+lrlXF0/N0L/S4BAICMcc3QniHHmTDRhbBYj3PuJefcZc651zyu/V3S4zWH1yS1sARxzqnsQCDkXCYP0AUAINHuGTvQ7xLijrDYNO/VfDzK1yriZNOucr9LAAAg41jYcbpPdCEsNk2/mo9f+FpFnAy//+WQYya2AADQcusybKILYbGRzKyHpOtqDp/xsZSEYWILAAAIR1hsBDPLlfRXSZ0kveic+5fPJbXYsF+85HcJAABkrEya6EJYbJzfSzpb0ueKMbnFzMaZ2TIzW7Zt27akFdcc4eMVmdgCAED8ZNJEF8JiA8zsYUk3KLju4tnOua3Rnuuce9Q5N8Q5N6R79+5Jq7Gp9uyr9LsEAACyzqQ5q/0uoVkIizGY2S8l3Sppm4JBca3PJcXFiT9fHHLcKZ+12QEAiLfwXrtZSzf4VEnLEBajMLP7Jf1E0g5J5zrnPvC5pLhxYcdFk0f5UgcAAEh9hEUPZjZD0m2SdikYFIt8Liluxv99hd8lAACQNfod1j70eOJ8nyppPsJiGDObJulnknYrGBTfa+AlaeW59zaHHF87tJdPlQAAkPle+MmZIceVgfD+vdTHYLV6zOwiSXfVHH4i6Raz8HXYJUkfOedmJK2wOHEu8h/otLEDfKgEAACkC1oWQ3Wp9/kQSd+J8vhG8ktrubvnpucsLAAA0ll4L17vNFtzkbBYj3PuceecNeJxpt+1NsespRtDjllbEQCAxAvvxUu3jmjCYpaoqAr4XQIAAEhDhMUs8d0/Lws59hyJCQAAEmLdfaNDjtNp+z/CYpZ4de32kON1dEEDAJA0USbMpgXCIgAAAKIiLGaB2/4RuqZ4+v5tAwBA+gqfWJouC3QTFrPAP97dFHJMFzQAAP5LlwW6CYsZLlCdHv8QAQBAaiIsZri1X+71uwQAAFDj2qE9Q44nzUn9DTMIixnu/F+/FnLMXtAAAPhn2tiBIcezlm7wqZLGIyxmuPBOaPaCBgAATUFYzGDOMV4RAIBUc0jbXL9LaBLCYgY79d4XQ477F3TwqRIAAFBr5ZRRIcepvpsLYTGDfbm3IuR48fgRPlUCAADSFWERAAAAUREWM1RFVcDvEgAAQBThu7mk8hI6hMUM9bWfvxByzJI5AACkrlReQoewmKHKDoS2LLJkDgAAaA7CYgbaX0kXNAAAqS4vx/wuoVEIixmodU7o/9ZO+em1nhMAANlg7fTRIcfnPbTEp0piIyxmoGMmLQg5Lpo8KsozAQBAqlhTXOp3CZ4IixmoMsDOLQAAID4IixmG8YoAAKSPjm1Sf6gYYTHDrNy0J+SY8YoAAKSuVVNDh4ql4nqLhMUMc8Pj74QcM14RAID0kYrrLRIWM8zeiiq/SwAAABmEsAgAAOCjVF9tkbCYQb7cu9/vEgAAQBOtnpraQ8YIixnkp7OLQo7ZDxoAgNTXPmxGdO8J83yqxBthMYO8tnZ7yDH7QQMAkH5SbbVkwiIAAACiIiwCAAD47BvHF/hdQlSExQyxe9+BkGMW4wYAIH08fNVgv0uIirCYIcLHK7IYNwAA6aNNbk7Icb+J832qJBJhMUP8+G/v+V0CAACIk8pA6kxzISxmiBT6NwUAADIIYREAACAFHNU53+8SPBEWM0BVoNrvEgAAQAs9f8swv0vwRFjMAMs37g45ZucWAADSz6HtWvtdgifCYgZ4+MU1Icfs3AIAQPorTJFt/wiLGeCNT3b4XQIAAMhQhEUAAABERVgEAABIEeccd5jfJUQgLKa58gMBv0sAAABx8sBlJ/hdQgTCYppbt70s5Lh/QQefKgEAAC3VuX3qzYgmLKa5wm7tQo4Xjx/hUyUAACDezntoid8lEBbT3cUzX/e7BAAAkCBrikv9LoGwmO7WflnW8JMAAACaibAIAACAqAiLAAAAKeQbAwr8LiEEYREAACCF/PrKwX6XECI33jc0M5PUT9Jhkg6V1FHSAUmlkjZL+sw5ty/e7wspL8f8LgEAALRQm9zUasuLS1g0s9MlfVPS1yUNkJTfwPM/lfS6pEWS/kV4jI+100f7XQIAAMgwzQ6LZtZB0nhJ35fUQ5JJcpKKJb1f87Fc0n5JrRUMkF0lfUVSoaS+kr4jaZ+Z/V3SL5xza5tbTzbatIuMDQBApgl20qaOJofFmm7mn0r6maQukpZJ+l8FWwqXO+dKGnGPNpJOlDRU0gWSrpV0nZn9WdL/OOd2NbWubPTGJ9v9LgEAACRYv4nzfe09bFJYNLMjJD0l6WuSHpX0O+fcp019U+dchaS3ah4Pm1lXSd+TdLOkMWZ2tXPupabeN9tM+ef7fpcAAAASrDLgfH3/Ro+gNLOvSSqStF5Sf+fcT5sTFL0453Y45+5VcGLM/5M018z+Ox73zmTlldV+lwAAADJco1oWzewYSX+X9H3n3DOJKqZmostkM5sr6SkzK3fOPZWo9wMAAEBsjW1ZnCXp8kQGxfqcc8slDZN0i5n1ScZ7AgAApIruHVv7XUKdxobFkc65FQmtJIxzbpuksyRtSOb7AgCSqKJUmv2d4EcAdR677hS/S6jTqLDonCtLdCFR3ne/c67Kj/cGACTBuiXSB3Okda/6XQmQUnp0aut3CXUSskS4mbUys2FmdoWZDa9ZkxEAgFAf/iv0IwBJUrcObfwuoU4itvu7TNLDCi7UXausZg3FiY1ZhxFN17+APA4gzTgnrVkU/HztouBxii1GDPiltCJ1Olbj2rJoZjdImi3pcEk7JZUouLNLB0k/kFRkZifE8z0RtHj8CL9LAICm2faRVLU/+HllubTtY3/rAVJIKu0PHbdKzKy3pJmSnpDU1znX3TnXWVInSVcquA90L0kvmNnx8XrfbMVWfwDS3prFUnUg+LkLSGsX+1sPkELyctIwLJpZuwae8iNJTzrnrnXOfVZ70jm31zn3D+fcaElnS6qS9IyZ5TerYkiSikv2+10CALTMB89JgYrg51UV0vvP+VsPAE9Nia3LzKzIzB4xs2+aWbew6+dImhjrBs65lyWdKilP0s+bVirq69W1vd8lAEDzVeyVisO2LC1ezRI6QApqdFh0zn1V0gWS3pY0StLrZrbKzGaa2eWSzDm3tRH3+VzSJZJuNLPuzaw766XSLCkAaLJ1r0q5YT/HctuwhA6QgprUIe6c+9w591fn3Djn3LEKtia+KmmEpH5m9kFNeLzEzDrHuM9KBXeFuaYlxWczxiwCSGsf/CuyFbGiVPrwn/7UA6S4SXNW+/beLRo96Zwrds7Nds79UNIKSWdKekPSGAW7rd8xsxlmdp7HGMXf1jwPzbBq0x6/SwCA5nEuuFSOXPiF4FI6Lvw8gFlL/dvQLp7rLFZJynPOPSXpKUkys6MVbH28QdKfzOwTSS/WPN6SFLX1EbFt2lXudwkA0Dz1l8wJV7U/uITOYccmtyYAUcVzXvYiBWdE13HOfeace9Q5d6Vz7qia62WS7pK0SVJ/MxtvZifGsY4WM7Orzew1M9tjZqVmtszMbjazlJnH/tALa/wuAQCap/6SOeFYQgdIOfEMP39VcNJK72hPcM4VOed+VbOMTk9JX0hqL+lBM/vczJ4zs1v9XIfRzH6r4FqRQyS9JukFSf0VXEPyaTPL8au2+vZVRvlBCwCprv6SOeFYQgdIOXHrhnbOfWpm8yX9w8zOcM41NAOjnaSAc+4eSfeYWZ6Cy+qcKekeBWdMJ5WZfVPBnWa2SjrDObe25nyBpJdravqhgtsZAkD2+vs1zd/POad17OvFq6UpnZp37+MulK78a/NeC8BTvLtVx0s6StIrNeMVY/mpgi13kiTnXKVz7nXn3D3OuaQHxRp31Hz8WW1QrKmtWNJNNYcTUqk7GgB8cc5UqcdAKa+h/Ro8BA607LqXvHZSj0HBugDEVVxDj3Num6SrJQ2QtMrMfhXepWxm7cxsqqTbJf0unu/fEmZ2lKSTJB2Q9I/w6865JZI2S+ohaWhyqwOAFNO1jzRuiTTyTik3X/Lrb2hrFXz/kROD9XTt408dQAaL+3e3c+4lSZcqODv6R5JWmlmxmb1pZm8r2MV7l6R7nXMr4v3+LTC45uP7zrloU43fCXsuAGSvVjnS12+RbnpDKji+ea2MLZHXTioYEHz/r/9QakWnD5AICfnOcs4tlHSCgjOkJam7guMRhyi4sNZPnHOptt1f7cScWAsZbQx7LgCgrpVxYnJaGWlNBJIqnusshnDOrZc02sz6SjpdUoGCQWyhcy4VV5TuUPOxLMZzarcb6JjgWgAgvbTKCbbuHXO+9I/rpB2fSJUJ2Gkqr53UtZ90+WOERCBJEhYWaznnPpH0SaLfJw6s5mOztw4ws3GSxklSz54941ETAKSX2lbGpb+TXronuESOq275fa2VlNMm2Jo49Ad0OQNJxHfbQXtrPnaI8Zzaa3u9LtYsQD7EOTeke/fucS0OANJGq1bBVsab3giOKWzpWMa8dlLBQMYmAj7hO+6g9TUfe8V4zlfCngsAiKa2lXHYeCm3bfPukds2+Ppxr9DtDPikwbBoZoVmdlkyivF479uT+Hbv1Xw83szyozzn5LDnAgBiadVK6lDQ/EkvliN17EFrIuCjxnz3bZb0fTP7QaKLqWVmeWY2S1LnZL2nc+5zScsltZZ0uUdNIxRccHyrpDeTVRcApL1lf2z+ZJfKMumdP8a3HgBN0mBYdM5VSvqmgvs+zzKzhM4ENrNjJL0lqco5d0dDz4+z+2o+/qJmFndtTYfp4ALiM5yLx2htAMgCezZJX37Usnt8+aG0Z3N86gHQZI1q169Z6uYMSV0kfWxm3zezuM6kNrNuZvYrSUWSFjjnro/n/RvDOfe0pP9VcJeWVWb2LzN7VtJaSV+VNEfSzGTXBQBpa+Vsyazh58ViFrwPAF80ehCIc65U0gWSfiXpQUmfm9m9ZnZic9/czFqb2flm9hdJnyvY/Xulc25ic+/ZUs65H0j6LwW7pEdIGqXg0j8/lPRN51zAr9rqa5vL+B0AaeDdx6Sq/S27R9X+4H2ALNa/INZiLYnVpNZB55yT9KCZPSNpsqT/kfQzM9sm6T+SVirYCrdF0nZJ+yVVSmojKV/Bhbm/IukYBfdhPllSOwWXovmVpF8450pa/mW1jHPuSUlP+l1HLI//9yn61qNL/S4DAKLbuloq2xb7ObXrJza0HmPZl1Lx+8FtBYEstHj8CN/eu1ldyc65dZKuM7O7JV0v6SpJY2sejVnU2mqe96akpyT9NUV3dUlZQ4/u6ncJABDbiiekQFX063ntpK59pXOnSi9Mjr3rS6AqeL9R9yamVgBRtWjcoXNuo6SpkqaaWU9JwyWdKKmPpCMU3BavjaR9kkoU3O7vE0nLJL3unNvVkvfPZnv3V/pdAgBEV10tFT0lVXv8rPLajaX3iJpdX6Z7tzJWV0ornpLOvYdldIAki9sklZrg+ETNAwlWzXxsAKlswxtSwCMo1rYmXv546CLbrXKkr98iHTNa+sd3pB2fRrYyBg5IG/8jFQ5LaOkAQvHnWZr6ZJvnjoMAkBqW/0U6UHbw2FpJufnB1sRxS6LvxlK768vIO4PPr7+Yd+W+4H0BJFXcWhbNbJqCk1Zek/SipGWsR5g4h7Zr7XcJAOCtcr/00b9UN4Q9WmtiNNFaGV219OG/pKoKKbdNIr8CwHeVgdSJUPFsWRyn4DIz9yg4cWWnmc01sx+Z2cA4vg8k9enu3xR6AIhp7eLgNn2NbU2MxquV0VoF7w9kuF1lB/wuoU48w+Jlkp5VcMkck3SIpAsVXBJnhZkVm9nfzOy7ZsZu8ACQqZb9KdgFXTBAuukN6es/bP6klNpWxpveCC6bc6AseH8gw5VXpsSyzpLiGBadc6855y53zhVIGijpFknPSNqhYHjsLukKSb+XtMbM1pvZH83s7HjVkM3Oe2iJ3yUAQFDrdtJ59zSvNTGa2lbG86YFu7WBDGdq4c5HcRTXLftqOefel/S+pN9KkpkNkHRmzeMMSd0k9VRwjcbrzGyFpBudc+8lop5ssKa41O8SACDoWwna06C2lfHrtyTm/kAKaZ1CO7UlpRLn3Grn3Ezn3GWSDpf0a0llkhZKOiBpsKSlZnZZMuoBAABIZam0nnLSY6tzLuCc+4mkv0qaJ+lISRMUDI9Pmdk5ya4JAAAglVz1h9TZ0tfPNs47Jf23c26nc+5+Sf0kvSTpMTM7xMe6AAAAfLW9NANnQ5tZOzP7jZmtNLNnzOz8WM+v2eqvdb3jHZLGSCqS9N/xqgsAAADNF8+WxUck3SxpgKRLJD1vZuvM7C4zOzr8yWbWS8G9o+s456pq7nFFHOsCAABAM8VzNvS5Cq612E7SaEkXS+olaaqkqWa2VtJSSVskdZF0qaRF4Tdxzm0ws9SZApTC8vNaqbwydVZ4BwAA8ef3IjrxDGW7nHPPOuf+6py7WtIRCq61WDtCs7+kb0v6mYK7vVQoOG5RZnaRmT1gZkfUPLdtHOvKWP+6ZZjfJQAAgARbN2OMr+8fz5bF+WY2wTk3Q5Kcc3sUXGfxt2Z2pIKtjcdLaiPpI0mznHM7a177v5J6SBphZj9UsPURDeh7WMeGnwQAANJKaUWV3yWEiGdYvFvBYLhA0mTn3Nu1F5xzmyX9IcZrX5V0paSTJD0t6aI41gUAAJA2KqtSa4hZPLf7q3LOfU/SdAXHKH5gZo3dk+namsdtkkY451bEq65scsLUiCGgAAAgzewpT50FuaUEbPfnnHtd0vlm1tU5t6+Rr6mS9ES8a8k2e8pTq9kaAAA0XUkK7d4iJXBR7pp1EwEAANAEF818w+8SQrBEDQAAAKIiLAIAACAqwmKaO6Rt3IedAgCAFNEp3//f84TFNPfMTV/3uwQAAJAgRZNH+V0CYTHdHdk5P+R40pzVPlUCAABa6suS/X6XEIGwmObatQ5tnp61dINPlQAAgJZ687PUW0yGsAgAAJAifvS31NuXhLAIAACAqAiLAAAAiIqwmAGOOjS/4ScBAICUVhWoDjlOhWVzJMJiRljw4+F+lwAAAFpo/Y6ykONUWDZHIixmhI5t80KO+02c71MlAACguV766Eu/S/BEWMxAlQHndwkAAKCJ7p3/kd8leCIsAgAAICrCIgAAAKIiLGaISwYf6XcJAACgmbaXVoQcXzu0l0+VRCIsZoiHrjwx5JhJLgAApI9F728NOZ42doBPlUQiLGYoJrkAAJA+Jj632u8SoiIsAgAAICrCIgAAAKIiLGaQ8Ekuk+akbpM2AAAIKi7ZH3Lcv6CDT5V4IyxmkPBJLrOWbvCpEgAA0FjPLt8ccrx4/AifKvFGWAQAAPDRLxam5s4ttQiLAAAAiIqwmGHa5eX4XQIAAMgghMUMs2zSOSHHJ0xd5FMlAACgIas27Qk5Xj9jjE+VREdYzDDtWueGHO8pr/KpEgAA0JALZ77udwkNIiwCAAAgKsIiAAAAoiIsZqCrT/1KyDHjFgEASD17yitDjlNxvKJEWMxI914yKOSYcYsAAKSeE6Yu9ruERiEsAgAAICrCYobKyzG/SwAAAFFUBar9LqHRCIsZau300SHH/SbO96kSAAAQ7uYnl4ccp+p4RYmwmDUqA87vEgAAQI1F7xf7XUKjERYBAACSyLn0asAhLGawvoe1Dzk+76ElPlUCAABqrd5cEnKc6vMMCIsZ7N8/OTPkeE1xqT+FAACAOuFb/IXPM0g1hEUAAABERVjMcOEN25PmrPalDgAAkH7jFSXCYsZbO/38kONZSzf4VAkAADjzgZdDjlN5yZxahMUMl5vD/2IAAFLFhp3lfpfQZCQJAACAJEjHLmiJsJgVrh3aK+SY3VwAAEi+m/76bshx/4IOPlXSNITFLDBt7ICQY3ZzAQAg+RaG7dqyePwInyppGsIiAABAgu2vDPhdQrMRFrNE+G4uvSfM86kSAACyz7i/LAs57pSf61MlTUdYzBLhu7nQEQ0AQPK8unZ7yHHR5FE+VdJ0hEUAAABERVisYWbHmNl4M1tgZp+Y2X4z22Nmb5rZj82std81thSzogEASL5bn3ov5DiduqAlwmJ9L0r6laQzJX0h6TlJ70o6UdJDkpaaWRffqosDZkUDAJB8/yzaEnKcTl3QEmGxvo8l3SCpu3NuuHPuKufcWZKOk/S+pMEKhsaMwl7RAAAkzpcl+/0uocUIizWcc2c75/7knCsNO79e0vdrDq9I9+7o8K5o9ooGACBxht73YshxuizEXR9hsXFqBxu0ldTVz0JaKrwrGgAAJE512IivdFmIuz7CYuP0q/l4QNJOPwuJBws7Pu+hJb7UAQBAJkvnhbjrIyw2zoSaj8875yp8rSQO1s0YE3K8prg0yjMBAEBzHTtpYcjx+rDfv+mCsNgAM7tO0pWS9km6s4HnjjOzZWa2bNu2bckoDwAApKADVdV+lxA36bXQTxRmdr+ki5rx0rOdc5tj3PdsSf9PwQ1Pvuec+zjWzZxzj0p6VJKGDBmS0uvS5OVYyNI5vSfMi2hxBAAAzTPywVf8LiFuMiIsSjpC0jHNeF1etAtmNkzSXEmtJd3qnPtrM2tLSWunj1Zhvf2hUzrZAgCQZjbvLg85TtcuaClDuqGdc9c456wZj/Ve9zOzr0uaL6m9pJ85536TzK/HL0x0AQCg5aoCmdMFLWVIWIwnMxsqaYGkjpLucs7d73NJCXPN0J4hx0x0AQCg5fpOXBByHL7GcbohLNZjZqdIWiTpEElTnHPTfS4poe4ZO9DvEgAAyCiVHq2K6b7GMWGxhpmdJGmxgkFxmnNuqs8lJUVeTuiqi73rjWMEAABNc9xdCxp+UprJlAku8fCCpE6SdkvqaWaPR3ne/zjntietqgRjogsAAPFTFfaLNJ0nttQiLB7UuebjoZK+E+N5UyRlTFj0csLURSqaPMrvMgAASCulFVV+l5AQdEPXaOkM6nQW/lfPnvLM/McOAEAiDZi8KOQ43Se21CIswtOkOav9LgEAgLSxcce+iHPpPrGlFmERkiL/+pm1dINPlQAAkH7OeODlkOPwCaTpjLAISZnz1w8AAKlg7fTRfpcQN4RF1An/K6iQZXQAAGhQpv++JCyiTib9FQQAQDJ4be2XCcvl1EdYREws0g0AQHTHT17odwkJR1hEiPC/hlikGwCA6CrCVuHOtFZFibCIRug3cb7fJQAAkHLOevAVv0tICsIiIoT/VVQZoH0RAIBwn20vCznOxFZFibCIRqJ1EQCAg743a5nfJSQNYRGewhfppnURAICDFr1fHHKcKVv7eSEswpPXIt20LgIAIN38xLsR5zJ5cwvCIqKidREAgEjzVm0NOc7kVkWJsIgYaF0EACDUj//2XsS5TG5VlAiLaAAzowEACHLOac6KLSHnMr1VUSIsohnY1QUAkI0GTF4UcS7TWxUlwiIaIfyvJtoWAQDZqOxAIOQ4G1oVJcIiGsHrr6ZCWhcBAFnEq1ctG1oVJcIiGilb/noCACBcRVUgolctm34vEhbRKLQuAgCy1TF3LYw4ly2tihJhEU3gteflpDmrfagEAIDkqKgKRJzL1D2goyEsokks7HjW0g2+1AEAQDJ4tSpmG8IimmSdx19T5z20xIdKAABIrANV1RHnsq1VUSIsohnyckLbF9cUl/pUCQAAidP/rgV+l5ASCItosrXTR0ecYxtAAEAmGTbjpYhz2diqKBEW0UzhSwawDSAAIJNs2l0ectwpP9enSvxHWESzsJQOwTiH5wAAH5xJREFUACBTef0+K5o8yodKUgNhEc3mtSApk10AAOlsZ9mBiHPZ2v1ci7CIZps2dkDEUjpMdgEApLOvTXvB7xJSDmERLeK1lA6TXQAA6ehHT70XcS7bWxUlwiLioH9Bh5BjJrsAANLR3KItIcfhvWfZirCIFls8fkTEOSa7AADSidfvLa/es2xEWERceDXTM9kFAJAONu7YF3GO7ueDCIuIG3Z2AQCkozMeeNnvElIaYRFx47WzC93RAIBU5vV7ilbFUIRFxBVrLwIA0kVloDriXPikTRAWEWdeO7vQHQ0ASDXOOfWbuCDivNekzWxHWETceTXf0x0NAEglve+IXBOY7mdvhEUkhFcz/qQ5q32oBACAUKUVVRHnwidp4iDCIhLCqxl/1tINPlQCAECoAZMXRZzzmqSJIMIiEobuaABAqmH2c9MRFpFQXt3RzI4GAPjh852Ri2/T+dwwwiISyqs7mtnRAAA/DL8/cvFttvRrGGERCUd3NADAb3Q/Nx9hEUnh1R19wtTIAcYAAMTbz54pijjXKT/Xh0rSE2ERSeHVHb2nvIrldAAACbXvQJX+/s6miPNFk0f5UE16Iiwiabya+1lOBwCQSF+9O7IXi+7npiEsIqkYvwgASBav3y/XDu3lQyXpjbCIpPMaJ8JyOgCAePIapyhJ08YOSHIl6Y+wiKTzGieypriU8YsAgLgoPxDwHKdI93PzEBbhC8YvAgAS5bi7F0acIyg2H2ERvvEaN8L4RQBASzBOMf4Ii/DNtLEDlJcTudES6y8CAJrj7F++4nmecYotQ1iEr9ZOHx1xjvUXAQBNtXvfAX26rSziPN3PLUdYhO8YvwgAaImqQLVO/PkLEecJivFBWERKYPwiAKC5+k5cEHGOoBg/hEWkhGljB3juH01gBADE4vV7InI0PFqCsIiUsXj8CM9vcCa8AAC8/PhvKzzPr6NVMa4Ii0gpXt/ge8qr2OEFABDis22lmrNic8R5up/jj7CIlOP1jc4OLwCAWlWBap31y8hGBIJiYhAWkZK8JrwwQxoAIHlPaGHh7cQhLCIlMeEFAODF6/dAp/xcFt5OIMIiUtbi8SPUKT834jyBEQCyU7Sf/0WTRyW5kuxCWERKK5o8ynOGdG8CIwBkleMnL/Q8zzjFxCMsIuV5zZB2kvpNnJ/8YgAASfc/s4tUVhGIOE9QTA7CItKC1w+EyoBjDUYAyHA/mb1CTy/fFHGeoJg8hEWkDa8fDHvKq1hSBwAy1IYdZXp2OWsp+o2wiLTi9QNi1tINBEYAyDD7KwMa8cArEedZIif5CItIO6zBCACZzTmnYydFTmi5dmgvlsjxAWERaYc1GAEgc1UFqtX7jsgJjKyl6B/CYgxmNsDMKszMmRn9nCmENRgBIPNUVzvP3Vkk1lL0E2ExCjPLlfS4pDyfS0EURZNHKS8nchVGAiMApKej7/ReEo0JLf4iLEZ3p6STJP3O70IQ3drpowmMAJABov3cJij6j7DowcwGSbpL0rOSnva5HDRg7fTRnru8EBgBIPU55wiKKY6wGMbM8iT9WVKppB/4XA4aad2MMQRGAEhDXpNZJIJiKiEsRrpL0omSfuycK/a7GDSe17aAEoERAFIVLYrpgbBYj5kNVnCs4kLn3F/8rgdNF+0HDIERAFILQTF9EBZrmFlrBbufyyWNa+Y9xpnZMjNbtm3btrjWh8ZbT5c0AKQsxiimn4wIi2Z2v5l91IzHkfVuc7ekgZJud8593pw6nHOPOueGOOeGdO/ePS5fG5qHLmkASD3OOcYopqHIVY3T0xGSjmnG6/IkycxOkvQzSa9I+n/xKwt+Wj9jjGc4LJwwjx9KAOADgmJ6yoiWRefcNc45a8Zjfc0tLlQwOBdIetnMXql9SPp1zXN61zvfN/lfJZqDLmkA8N+Bqmq6ntNYRoTF/9/enUfLVZb5Hv8+JGHQMKgoOCBEIYKESYJNey9GDA3aqAtacUIXrna4jXodUYMKQqOI9EVxOTZXbWi43Q6o2BoIoAKNLfESRkGj4oUICAEChMzDOc/9Y+8jh5OqOnXOqapdw/ezVtZOau/a9aTWObt+9b7vft8W2geYN+bPAeW+J416bMuFidW1GnVJn3KJqzhKUjutWr+J2Z+svYSfQbE3GBaBzDytXusjcHh52O2jHr+5yno1cXeddXTNlV4uXLysgmokaTCs3rCZ/U67ouY+g2LvMCxqYDRa6cUWRklqrfWbhpjzqctr7jMo9hbDogbKnWcdzY7bbXlf14WLl3HA6bUvapKkiXn/t29i71MWbfF4YFDsRYZFDZxbPnUUbz109y0eX7luM7O88UWSpmTfUxfxo5v/vMXjM6ZF3THk6m6GxXFk5tXlOMU5Vdei1jnjmDk1A2OCgVGSJmmPBQtZs3Foi8dn7zKTP3zmbyuoSK1gWNTAOuOYOTW7QxKn1pGkiWi0KstbD92dKz44r8MVqZUMixp4rictSZO3aWi44WTbZxxjx1yvMyxKNA6M3iktSbWt2bCZvT7hHIr9zrAolepd2C5cvIwjv3BNh6uRpO72sYtvYV+nxhkIhkVplHrLA/5++WpvfJGk0r6nLuI7S+6puc+g2H8Mi9IY9eZi9MYXSap/x7NzKPYvw6JUQ725GMHAKGkwjXfHs3Mo9q/IzKpr6Etz587NJUuWVF2GWqDRxdG7/CQNggceW8+Lz/xZzX22JvaHiLghM+fW2mfLojSORje+OI5RUr/b55TLDIoDzrAoNeEuxzFKGkB7LFjIuk3DNfcZFAeHYVFq0njjGJ2PUVK/WLl2U90vwrN3mWlQHDCOWWwTxyz2t3oX0RnTwvVPJfW0Q8/8Kfc/tqHmPkNi/3LMotRi9S6Ym4bq3y0oSd1ujwULDYragmFRmqR64xjBcYySesvGzcN1r1vOnyi7odvEbujBYre0pF71/JMXMlQnChgSB4fd0FKb2S0tqdeMTLJtUNR4DItSi9x11tHM3mVmzX3eLS2pm3zs+7cw6+RLa+6bMS0MinoCu6HbxG7owdaoNdGLsKSqDA0nz/947ZAIXp8Gmd3QUoeNd/OLrYySOu3uh9caFDUphkWpTRpN4n3h4mWOZZTUEcPDxdjEw86+qub+tx66u0FRDdkN3SZ2Q2u0RsHwrYfuzhnHzOlgNZIGxQ3LHuG1X/tl3f2GRI2wG1qq2F1nHW0ro6SOGbnTuV5QdMk+TYQti21iy6LqadzK+FzOOGa/DlYjqd989OJb+O6Se+ruNySqFlsWpS7SuJXxT7YySpqU9ZuG2GPBwrpB0bGJmixbFtvElkU1o1EwDOBOL+ySmjBrwUIafZobEjUeWxalLtWolTFxmh1Jjd3zyFr2aBAU7zrraIOipsyWxTaxZVETdcDpl7Ny3ea6+73gSxqxeWiYPT9xWcNjvGZoIhq1LBoW28SwqMkab8yiHwDSYNvz4wvZPFx/v9cITYbd0FIPadQ1DUWYPOD0yztYkaRucMcDq9ljQf2g6HQ4ahdbFtvElkW1wnitjE7oLfW/R9Zs5KAzrmx4jCFRU2U3dAUMi2olu6alwbN+0xB7n7Ko4TH+7qtVDIsVMCyqHQyNUv/LTGadfGnDY+xVUKsZFitgWFS7nHLJbVy4eFnd/TtuN51bPnVUByuS1Crj3bzi/KtqF8NiBQyLarfxptoxNEq9Y99TF7Fm41DDY+w5UDsZFitgWFSnjNc1bWiUutc7L7ieK3/7QMNjDInqBMNiBQyL6jRDo9Q7TvreLVx8Q+01nEcYEtVJhsUKGBZVFUOj1L1OvOgGLrvt/obHGBJVBcNiBQyLqpqhUeoe886+imUPr214jCFRVTIsVsCwqG4xXmicMS34w2f+tkPVSINjeDh53scbT4EDhkR1B8NiBQyL6jbjhUbwQ0tqhU1Dw+z1icvGPc7fN3UTw2IFDIvqVs2ERif8lSZu6f2P8Ypzrx33OEOiupFhsQKGRXW7ZkLj7F1mcsUH53WgGql33fHAao74/DXjHmdIVDczLFbAsKhe0UxodNUI6YnWbtzMC0+9vKljDYnqBYbFChgW1WvGWxFmhB98GmQrVm/g4E//dNzj/IKlXmNYrIBhUb1qvLWnRzM4ahCs3zTEnFMXsbmJj0t/J9SrDIsVMCyqHzTTRQ3eEKP+dPfDazns7KuaOtaQqF5nWKyAYVH9pNkuavBDU73tTyvW8tJ/ai4gOrG9+olhsQKGRfWrZlsbnexbvWLz0DB7f/KyprqZwS9E6k+GxQoYFtXvJjK20RYYdZvNQ8Mc/41f8as7H27qeH+G1e8MixUwLGqQNNvaCLDDttO59TQ/dNV5m4eGuW/l+qbHIYKtiBochsUKGBY1qCYUHLebzq221qiNNmwe4qWfu4rlqzY0/RwDogaRYbEChkVpYsER/JBWa/x++SqO/MJ/Tug5/uxp0BkWK2BYlB43kfGNI2751JHsuN2MNlWkfjKZ7mUwIEqjGRYrYFiU6ptoiyP4wa4nWv7Yev76sz9jeIIfYf4cSbUZFitgWJSaM5ng6LQ8g2fF6g3MP+caHl23acLPNSBK4zMsVsCwKE3ckV+4ht8vXz3h503fCpae8UqmT9uqDVWp0zKTW+9ZyRvPu451m4Yn/HynuZEmzrBYAcOiNHWTaXUcsfjk+ey647YtrEbtsmbDZn5z32Mc9/XrJn0OWw+lqTEsVsCwKLXeVMLjjK2C/1rwcp6+/TZERAur0kRs2DzEn1as5ZXn/mfTK6bUYjiUWsuwWAHDotR+UwmPIxZ94DBm7fxktpk+rQUVaURm8sCqDaxav4kjPj+xaWxqMRxK7WVYrIBhUeq8WQsW0qor2oVvfzGH7PFUtp1hiGxk/aYhblz2CA+u3sD7v31zy85rOJQ6y7BYAcOi1B1a0fo41l8/76l85fiD2X7b6czo85tqMpOHVm/krhVruHHZI3z2sqUtfw2DoVQ9w2IFDItS99rrE5eyaai9174Dd9uJd730ecza+ck896lPYsPmYXbYdnpX3LE9NJys2zTEshVr+N39q1izcYg/PrCa8395V9tf22AodSfDYgUMi1Lv6USInKiDd9+J982fzW33ruSPD65m5jbT2SqCw/d+Br++51FuvWclP/3t8glPTt1uhkKptxgWK2BYlPpLK8dD9gsDodQ/GoXF6Z0uRpJ60Z1NBKNubJmcDFfIkTSaYVGSWmSqAeuA0y9n5brNU65j9i4zueKD86Z8HkkCw6IkdQ2XqJPUjaq/LU+SJEldy7AoSZKkugyLkiRJqsuwWENEHBURP4yI+yJiY0Qsj4hrI+KkqmuTJEnqJMPiKBGxVUScBywCjgZ+D3wfuB14AfAPFZYnSZLUcd4N/URnAu8EfgW8PjP/NLIjIqYBL6qqMEmSpCoYFksR8ULgJGAF8KrMfGj0/swcAq6vojZJkqSq2A39uBOBacA3xwZFSZKkQWXL4uOOLLeXR8SuwJuA2cA6YAnwg8xcX1VxkiRJVTAsAhGxDbBX+c99gEuA7cccdndEHJuZN3S0OEmSpArZDV14ChDl388FlgIvoQiM+wH/AewGXBoRO9c7SUS8KyKWRMSSBx98sM0lS5IktV9kZtU1TFlEnA28ZhJPnZ+Z90bEs4B7y8dWAM/PzJWjzr8VcBOwP3BaZp4+3onnzp2bS5YsmURJkiRJnRURN2Tm3Fr7+qUb+lkU8yBO1Ixyu2rUY98dHRQBMnO4nH/xy8B8YNywKEmS1A/6ohs6M9+SmTGJP3eVz19F0aIIcGedlxl5fNc2/3ckSZK6Rl+ExRa5sdw+rc7+kbGKqztQiyRJUlcwLD7uB+V2fp39R5RbByJKkqSBYVh83AXAPcDciDh59I6IeB1wPDAEfLWC2iRJkiphWCxl5jrgOIpu5jMjYmlEXBwRS4DvUUyt84HMvLnKOiVJkjrJsDhKZi6mmB7nW8BMiul4dgd+BMzLzC9XWJ4kSVLH9cvUOS2TmXcCb6+6DkmSpG5gy6IkSZLq6osVXLpRRDwILKu6ji60M/BQ1UX0Md/f9vG9bR/f2/bxvW2vfnp/d8/Mp9faYVhUR0XEknrLCWnqfH/bx/e2fXxv28f3tr0G5f21G1qSJEl1GRYlSZJUl2FRnXZe1QX0Od/f9vG9bR/f2/bxvW2vgXh/HbMoSZKkumxZlCRJUl2GRUmSJNVlWFTXiIijIuKHEXFfRGyMiOURcW1EnFR1bf0iIuZExIaIyIi4rep6elVEvCAiPhgRl0XEHRGxPiJWRsR1EfGBiNi66hq7XUS8ufz9XhkRqyNiSUS8JyL8XJqkiJgREfMj4pyIWDzqWnpvRFwcES+rusZ+ExFnltfT7OfPKscsqnLlh8PXgXcCm4DrgD8DuwBzgMcyc8/qKuwPETEdWAy8CAjg9sycU21VvSki7gGeDawHlgD3UPy8/jWwLXATcERmPlxZkV0sIr4CvJvi/fsZxe/9fGB74IfAcZk5VF2FvSkijgCuLP95P3ADsAZ4IcW1FOCMzDy1gvL6TkQcQvF5tRXFNfUjmfm/qq2qPfwGp25wJkVQ/BWwZ2bOy8w3ZebLgWcCb6q0uv7xceBg4KtVF9IHfkexhvzTM/OwUT+v+wC3AwcBX6iywG4VEa+lCIr3A/tn5qsy81hgL+C3wLHAeysssZcNA98HXpqZzyzf2zdk5n7AG4Eh4JSIOLzSKvtARGwDnA8sB35UbTXtZ8uiKhURLwRuBR4F9s7Mflk2qatExP4ULWA/Br4EXIUti20REf8duJai1WzHzNxYcUldJSKWUHxpOSEz/3XMvnnA1RRB8tmZOdz5CvtXRHyD4kvOtzLz7VXX08si4nPAR4HXAK8FTsCWRaltTgSmAd80KLZHRMwALgBWU7ToqL1uKrfbAk+rspBuExHPoQiKG4Hvjd2fmdcA9wK7Aod2trqBMPKz+ZxKq+hxEfFXwIeBf8vMH1ddTydMr7oADbwjy+3lEbErRZfzbGAdRUvYDzJzfVXF9YlPAgdStOQsj4h9qi6oz+1VbjcCjll8ooPK7e2Zua7OMddTjAc9CPhlR6oaHCM/m/dVWkUPi4htKb58Pwy8v+JyOsawqMqUYz5GLl77AJdQDHAf7e6IODYzb+hocX0iIg6iGKu4aGyXn9pmQbn9SWZuqLSS7jOr3C5rcMyfxhyrFii/jL+t/Of3Kyyl130GeAHwxkHqDbMbWlV6CsUdZADnAkuBl1AExv2A/wB2Ay6NiJ0rqbCHldO3XEDRSvuuissZCBHxNuANwFqKkK4nmllu1zQ4ZnW5HfvFUZNUzoRwEbAj8LNB6TpttYh4CfAB4JLM/E7V9XSSLYuatIg4m2Jw70TNz8x7eeKXlZXA32TmyvLft0XEsRRjbPYH3gOcPpV6e0kL3luAUylC94mZeXfLiutxLXpva513PvDPQAL/IzN/N8kS+9nIl0PvrOysr1NMTXQ38JaKa+lJEbEd8C/AYwzg2G/DoqbiWRTN8RM1o9yuGvXYd0cFRQAyczgizgO+THGhG5iwyBTf24g4GPgYxZ2l/9y6svrCVH9ut1DeAf0jYGvgfZl50SRr63cjv/MzGxwzsm9Vg2PUpIj4IsUd0PdTfOG5v+KSetWZFOPp/z4zB27Mp2FRk5aZb2EK31Izc1VErKC4Y/TOOoeNPL7rZF+nF031vQVeTfH7vQtwVUSM3rdTuZ0VEVeXf39HZt4xhdfrGS14b5+g7Jq6FHgy8LHM/FKrzt2H7iq3uzc4Zrcxx2qSIuIc4H3AgxRB8Q8Vl9TLjqWYx/KEiDhhzL69y+2JEfEq4I7MfEdHq2szw6KqdiPwN9SfYmRkrOLqOvvV2D7ln1qeBMwr/96opUd1RMShwGUU4+s+mZlnV1xStxuZumXfiNiuzh3Rh4w5VpNQDrf4ELCCYojPbyouqR9sxePXzFqeV/7ZqcExPckbXFS1H5Tb+XX2H1Ful3Sglr6RmadlZtT6A4ys3nD7qMdvrrLeXhQRLwYuB3YATsvMz1RcUtcrx87eSNFdf9zY/eWk3M+h6DK9rrPV9Y+IOAv4CPAIRVC8peKSel5m7tHgmnpBedhHyscOrLLWdjAsqmoXUKyrOzciTh69IyJeBxxPsUSVS9Spa5RjQq+gCIpnZOYgjaedqs+W289FxF/WfI+IZ/D47/lZrt4yORFxBsV45UcpgqIttJoyu6FVqcxcFxHHAVcCZ5ZjQW4D9qBY6SGB99vypS5zJcU0JI8Cz42I8+scd9IgzcXWjMy8OCK+RrF6068j4qfAJorehR0o5lv9coUl9qyIeA3FJPwAdwD/c8x45RFLM/OsjhWmnmdYVOUyc3G5dvEngaMopjVZSXF36TmZeW2V9Uk1PKXc7kSxJmw9pwGGxTEy890R8QuKKbHmUSz5uRT4FvA1WxUn7amj/j63/FPLNYBhUU2LTKe7kiRJUm2OWZQkSVJdhkVJkiTVZViUJElSXYZFSZIk1WVYlCRJUl2GRUmSJNVlWJQkSVJdhkVJkiTVZViUpC4XEVtHxKyq65A0mAyLktTFImJH4DLgkDadf2ZEXB8Re7fj/JJ6n2FRkrpURGwPXA1ckpnfncJ5fhsRM2vty8zVwLuAn0TEwZN9DUn9y7AoSV0oIrYCvg38NjO/NIXzPBfYXIbCmjLzJmABcHlEPGeyryWpPxkWJak7fQg4GDhxiueZB1wz3kGZeTFwKfBvETFtiq8pqY8YFiWpy5Q3s5wOnJGZK6d4upfSRFgsLQBeBHx4iq8pqY9EZlZdgyRplIg4H3glsFtmbpziuX4HHJaZDzR5/OeBE4BZmfnYVF5bUn+wZVGSukhE7A4cD5zfgqC4K8V4xaaCYumrwFOAD07ltSX1D8OiJHWXtwHTge+04FxNjVccLTPvAJYA/1DeZCNpwHkhkKTu8mbgQeCmFpxrIuMVR/sxsCvw31pQg6QeZ1iUpBaIiD0j4sKIuCkifhcRP4mIOeW+nct9t0bEDeXfd6lxjt2B2cAvcxIDyssaTouI6yLieor5E0+LiN9ExBcj4nlNnuoX5fa1E61BUv8xLErSFEXEy4AfAv+SmQcB+wLDwC8j4kBgEfDtzNwfOBc4FvhmjVMdVm5vneDrPyMiLqLoPp4JvB54FfDHzNwHOALYHvhNRLyziVPeCCTwkonUIak/Ta+6AEnqZRHxbOBrwJGZeTdAZm6OiG8Cr6ZopfvHzFwYEc8Hzqf4ov6sGqc7sNz+YQKv/3Lg34E/A3PLMYdExJuBn5f1/Bn4+4jYGjgvIlY2WhEmM1dGxHJgTkRMy8yhZuuR1H9sWZSkqfk08OmRoDjK9qP+/rVyux5YCawF/rHGuZ5fbv/czAtHxCuAhcA9FNPj3DFq93zgZ2Oe8qly+8WIGK+x4G5gO2CvZmqR1L8Mi5I0SRHxJOBwimX5xnpRub0mM1cBZOa9wDOBnTLzkhrPGWltfKiJ157F43dMv7HGcn6HU6wr/ReZ+UdgBcXNK+N1MS8vt88YrxZJ/c2wKEmTlJlrgf3qdNO+rNz+fMxzNmTmpjqnHGmNXN/Ey38F2AE4LzOf0G1ddnc/kpkrajxvpEWxVjf4aOvK7U5N1CKpjxkWJWkKRloNR4uIpwEHlP/8+dj9DWxTbuuFyZHz7wW8ovznN2occgRbdkETETsAO5b/vG+cWkYC6/YNj5LU9wyLktR6h1NcXx8GbpnA89aW220aHlXMfxjAQ5n56xr751M7pI4EzDWMP4/j1uV2SqvISOp9hkVJar355fbqzBweuzMiZtR53shazNuNc/6RORr/X41zB0WYvLbG895Wbv+1iXWfR2oYOxZS0oAxLErSFETE1jWWxTui3F5d4/hnUn8exWXldosJu+scN63GvgMo5ldcM+Z1DwdeSXGX88fHOT/A08rtvU0cK6mPGRYlaZIi4gSK1sClEbFN+dhsYM/ykOtrPO1E4P/UOeXI1De7jfPSi4BVwN4RMXPMvi2mzImI3crXfBg4JjMfHef8o2vYovVS0mAxLErS5J1DMb7wOTx+l/GHeXzs4drRB0fEIcDfAV+qc74l5fYFjV60DHunAk8GPjpm9xNubomI/SnWh14DHJ6ZNzY6d/mcbSnulv5jjSl5JA0Yw6IkTd5yinkLX5+ZayLiHRRrO78a2Az8ZWm9iDgauAB4Q2aurHO+X1Ass3fweC+cmecCpwELIuLTEbFTORbyIOBXEbF3RHyJYuziJcCBmdnsMoIHUoTfX4x3oKT+Z1iUpMn7O4owdnpE/BrYHzg6M38OHAUcGBFLy32vA47KzNvrnSwzHwauAw4pW/cayszTgUMpJtleDPyGYqqb/wt8FbgLmJ2ZHxo7hnEch5bbKybwHEl9KjKz6hokSaWIeC9FN/WrMnPhBJ97CrA2M8+ZYg1XUqzwsovd0JJsWZSk7nIRxU0zb5rEc18GXDWVF4+InYF5wHcNipLAsChJXaW8eeUrwHER0fS6zBGxNbA3cPMUS3gXxXjFz03xPJL6hGFRkrrP54CHaG4+xBF/BSypNQl4syLiycB7KSbtXjrZ80jqL4ZFSeoy5d3S7wHeHRH7NPm0lzHFLmiKaXimU0z/I0mAYVGSulJmXkIxj+O/j0z4PY65jJmMeyIi4kXAScDxmblisueR1H+8G1qSulS5zvP5FK19b8k2XbAj4unAfwH/lJn/ux2vIal3GRYlqYuVgfEs4JHMPKsN598O+Dnwlcy8qNXnl9T7DIuSNMDKVV/2yMw/VF2LpO5kWJQkSVJd3uAiSZKkugyLkiRJqsuwKEmSpLoMi5IkSarLsChJkqS6DIuSJEmqy7AoSZKkuv4/Nit6lrXU418AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(10, 10))\n",
    "plt.plot(xpos, ypos, 'o:', lw = '2', ms = '2', color='C0', label = 'geodesic')\n",
    "plt.plot([0.0], [0.0], '*', ms = 40, color='C1')#, label = 'orbit')\n",
    "plt.legend(loc = 'upper left', fontsize = 28, frameon=True, ncol = 1)\n",
    "#plt.axis([0,2,0,9])\n",
    "#plt.grid(True)\n",
    "#plt.title(\"Newtonian orbit\", fontsize = 28)\n",
    "plt.tick_params(labelsize=22)\n",
    "plt.xlabel(r'$x(\\phi)$', fontsize = 32)\n",
    "plt.ylabel(r'$y(\\phi)$', fontsize = 32)\n",
    "#plt.gca().set_aspect('equal', adjustable='box')\n",
    "plt.show()\n",
    "#plt.savefig('/Users/francescalepori/Desktop/orbit1.pdf', format = \"PDF\", bbox_inches='tight', edgecolor=\"k\")\n",
    "plt.close(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now we add the relativistic term in the equation of the orbit\n",
    "\n",
    "We add the relativistic correction to the geodesic equation:\n",
    "\n",
    "## $\\frac{dx^2}{d\\phi^2}= 1 - x + \\frac{3 G^2 M^2}{L^2} x^2$.\n",
    "\n",
    "## Exercise:\n",
    "Modify the code above to include the relativistic correction. Do you see any difference with the Newtonian result?\n",
    "Can you explain why?\n",
    "\n",
    "Solution below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x initial:  1.2054888751304351\n",
      "v initial:  0.0\n",
      "Amplitude of GR correction:  0.02\n",
      "Done.\n"
     ]
    }
   ],
   "source": [
    "# Routine that solves a 1st order differential equation using the Eucler method.\n",
    "def ode2Euler(f, v0, x0, phi):\n",
    "    '''Approximate the solution of the system\n",
    "       v'= f(x, phi)\n",
    "       x' = v\n",
    "       \n",
    "       by Euler's method.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    f : function\n",
    "        Right-hand side of the differential equation y'=f(t,y), y(t_0)=y_0\n",
    "    x0 : number\n",
    "        Initial value x(phi_in)=x0 where phi_in is the entry at index 0 in the array phi\n",
    "    v0 : number\n",
    "        Initial value v(phi_in)=v0 where phi_in is the entry at index 0 in the array phi\n",
    "    phi : array\n",
    "        1D NumPy array of phi values where we approximate x values. Angular step\n",
    "        at each iteration is given by phi[n+1] - phi[n].\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    x : 1D NumPy array\n",
    "        Approximation x[n] of the solution x(t_n) computed by Euler's method.\n",
    "    '''\n",
    "    x = np.zeros(len(phi))\n",
    "    v = np.zeros(len(phi))\n",
    "    x[0] = x0\n",
    "    v[0] = v0\n",
    "    for n in range(0,len(phi)-1):\n",
    "        v[n+1] = v[n] + f(x[n], phi[n])*(phi[n+1] - phi[n])\n",
    "        x[n+1] = x[n] + v[n]*(phi[n+1] - phi[n])\n",
    "    return x\n",
    "\n",
    "\n",
    "# Discretization of the variable t in Npoints in the range [tin, tfin]\n",
    "Npoints = 21000 #80000000  #210000000 we get it right with this number of points, but it's very slow\n",
    "Norbits = 10.0\n",
    "phi_in = 0.0\n",
    "phi_fin = 2.0*np.pi*Norbits\n",
    "phi = np.linspace(phi_in, phi_fin, Npoints)\n",
    "\n",
    "# Here I compute the adimensional factor 3G^2 M^2 / L^2\n",
    "# Everything is in unit of R = 10^{10} m\n",
    "rS = 2.95e-7\n",
    "a = 5.79\n",
    "L2 = 0.5*rS*(1.-0.2056*0.2056)*a #8.19e-7 \n",
    "\n",
    "fact = 3.0/4.0 * rS*rS/L2\n",
    "\n",
    "# Set the initial condition\n",
    "r0 = 4.6\n",
    "x0 = 2.0*L2/rS/r0\n",
    "print('x initial: ', x0)\n",
    "v0 = 0.0 \n",
    "print('v initial: ', v0)\n",
    "#Define the function on the right-end side of the differential equation\n",
    "GRamp = 0.02#rS*rS/L2*3.0/4.0 #0.01 #rS*rS/L2*3.0/4.0\n",
    "print('Amplitude of GR correction: ', GRamp)\n",
    "\n",
    "f = lambda x,ph: (1-x + GRamp*x*x)\n",
    "\n",
    "#Call my Eucler-solver\n",
    "x = ode2Euler(f,v0, x0, phi)\n",
    "r = 2.0*(L2/rS)/x\n",
    "r = r/np.sqrt(1.-rS/r)\n",
    "\n",
    "xpos = r*np.cos(phi)\n",
    "ypos = r*np.sin(phi)\n",
    "print('Done.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Let's visualize the resulting orbit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10, 10))\n",
    "plt.plot(xpos, ypos, ls = '-', lw = '2', color='C0', label = 'geodesic')\n",
    "plt.plot([0.0], [0.0], '*', ms = 40, color='C1')#, label = 'orbit')\n",
    "plt.legend(loc = 'upper left', fontsize = 28, frameon=True, ncol = 1)\n",
    "#plt.axis([0,2,0,9])\n",
    "#plt.grid(True)\n",
    "#plt.title(\"Schwarzschild orbit\", fontsize = 28)\n",
    "plt.tick_params(labelsize=22)\n",
    "plt.xlabel(r'$x(\\phi)$', fontsize = 32)\n",
    "plt.ylabel(r'$y(\\phi)$', fontsize = 32)\n",
    "#plt.gca().set_aspect('equal', adjustable='box')\n",
    "plt.savefig('/Users/francescalepori/Desktop/orbit2.pdf', format = \"PDF\", bbox_inches='tight', edgecolor=\"k\")\n",
    "#plt.show()\n",
    "plt.close(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Can we find the perihelion orbit?\n",
    "We need:\n",
    "- find perihelion \n",
    "- find the corresponding $\\phi_{p_i}$ for each orbit $i$\n",
    "- compute the difference $\\phi_{p_{i+1}} - \\phi_{p_i}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.6\n",
      "9.424778193434236e-07 4.599999920935773\n",
      "1.8849556386868471e-06 4.599999763889636\n",
      "2.827433458918449e-06 4.599999664274163\n",
      "3.7699112773736942e-06 4.599999683079214\n",
      "4.0500000999803574e-05\n",
      "60.47386512925216\n"
     ]
    }
   ],
   "source": [
    "list_perih = list()\n",
    "phi_perih = list()\n",
    "r_old = r[0]\n",
    "print(r[0])\n",
    "orbit = 0\n",
    "for i in range(1, len(r)-1):\n",
    " if (r[i] < r_old) and (r[i] < r[i+1]):\n",
    "  list_perih.append(r[i])\n",
    "  phi_perih.append(phi[i] - orbit*2.0*np.pi)\n",
    "  orbit += 1\n",
    "  print(phi[i]- orbit*2.0*np.pi, r[i])\n",
    " r_old = r[i]\n",
    "    \n",
    "sum_angle = 0.0\n",
    "for n in range(1, int(orbit)):\n",
    " sum_angle = sum_angle + (phi_perih[n] - phi_perih[n-1])*180.0/np.pi\n",
    "\n",
    "print(sum_angle/(Norbits-1)) # degrees per orbit\n",
    "\n",
    "#In degress per century on Earth\n",
    "angle = sum_angle/(orbit) * 3600.0 * 365.0/88.0 * 100.0\n",
    "\n",
    "print(angle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Higher order methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x initial:  1.2079587324981578\n",
      "v initial:  0.0\n",
      "Amplitude of GR correction:  7.963488286969253e-08\n",
      "[1.20795873 1.20795873 1.20795873 ... 1.20795875 1.20795875 1.20795875]\n"
     ]
    }
   ],
   "source": [
    "from scipy.integrate import odeint\n",
    "\n",
    "# Discretization of the variable t in Npoints in the range [tin, tfin]\n",
    "Npoints = 80000000 #21000000\n",
    "Norbits = 5.0\n",
    "phi_in = 0.0\n",
    "phi_fin = 2.0*np.pi*Norbits\n",
    "phi = np.linspace(phi_in, phi_fin, Npoints)\n",
    "\n",
    "# Here I compute the adimensional factor 3G^2 M^2 / L^2\n",
    "# Everything is in unit of R = 10^{10} m\n",
    "rS = 2.95e-7\n",
    "L2 = 8.196e-7 \n",
    "fact = 3.0/4.0 * rS*rS/L2\n",
    "\n",
    "# Set the initial condition\n",
    "r0 = 4.6\n",
    "x0 = 2.0*L2/rS/r0\n",
    "print('x initial: ', x0)\n",
    "v0 = 0.0\n",
    "print('v initial: ', v0)\n",
    "#Define the function on the right-end side of the differential equation\n",
    "GRamp = rS*rS/L2*3.0/4.0 #0.01 #rS*rS/L2*3.0/4.0\n",
    "print('Amplitude of GR correction: ', GRamp)\n",
    "\n",
    "f = lambda x,ph: (1-x + GRamp*x*x)\n",
    "\n",
    "#Call my Eucler-solver\n",
    "\n",
    "def geod(y, phi):\n",
    "    x, v = y\n",
    "    dydphi = [v, 1-x + GRamp*x*x]\n",
    "    return dydphi\n",
    "\n",
    "y0 = [x0, v0]\n",
    "sol = odeint(geod, y0, phi)\n",
    "print(sol[:,0])\n",
    "\n",
    "x = sol[:,0]\n",
    "r = 2.0*(L2/rS)/x\n",
    "#r = r/np.sqrt(1.-rS/r)\n",
    "\n",
    "xpos = r*np.cos(phi)\n",
    "ypos = r*np.sin(phi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiUAAAIDCAYAAADbmG7AAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd3hU1dYG8HelE2ogoQuhgzSBSFeqiqAidrzqFVS42K6oKFgAlXYtqNjxU7ChF1DkKkWUjtIC0nsJnQChJJQEkqzvj0lOZlImbSZ7yvt7njzZe+bMmTdokpV99tlbVBVEREREpgWYDkBEREQEsCghIiIiD8GihIiIiDwCixIiIiLyCCxKiIiIyCMEmQ7gqyIjIzU6Otp0DCIiIrdbt27dKVWNKu55WJS4SXR0NGJjY03HICIicjsROeCK8/DyTR5EpJSIvCAia0XkrIhcFJH9IjJDRDqZzkdERORrOFKSCxGpA2ABgPoATgBYCiAFQDSAvgA2AvjTVD4iIiJfxKIkGxEpDeB3APUAvAHgDVW9Yvd8JQCVDMUjIiLyWSxKcnoFtoLka1Udmf1JVU0AkFDiqYiIiHwc55TYEZEQAI9ldCeYzEJERORvOFLiqA1sl2YOqep2EekI4JaMx44DmK+qK00GJCIi8lUsShw1z/i8W0SmAvhntudHisiPAB5U1UslmoyIiMjH8fKNo4oZn68H8BCAt2G7AycCtrtujgC4E8BHub1YRAaJSKyIxJ48ebIE4hIREfkOFiWOMv89ggB8oarDVHWvqp5V1f8BuB2AAviniNTN/mJVnayqMaoaExVV7IXtiIiI/AqLEkdJdu3Psz+pqrEA1sH279a1hDIRERH5BRYljuLs2vvzOCbz8arujUJERORfWJQ4Wm/XzmuBtMiMz+fdnIWIiMiv8O4bO6p6RERWA2gHoAeAHfbPi0gEgNYZXbfutpeWlobExEQkJSXh0qVLSE9Pd+fbEfmtwMBAlClTBhEREShVqpTpOER+jUVJTmMB/A+223//VNUNACAiYQA+AVAetnklbluv5PLlyzhw4ADCw8NRoUIF1KhRAwEBARARd70lkV9SVaSmpiIxMRGHDx9GlSpVUK5cOdOxiPwWi5JsVPUXEXkbwPMAVmeMnCQAaAugOmy3BfdXVXXH+6elpeHAgQOIjIxERESEO96CiDKICIKDg1GpUiWULl0aBw8eRGhoKEJDQ01HI/JLnFOSC1UdBuAO2HYCbg6gN4CLACYCaKWqu9313omJiQgPD2dBQlTCwsLCUL58eSQmJpqOQuS3OFKSB1WdBWBWSb9vUlISKlSoUNJvS0QAypQpg5MnT4LrDBGZwaLEw1y6dAk1atQwHYOoRKWlK85evIwzFy/j4uU0YznKhgbgbNIF1EpXBARwDhdRSWNR4mHS09MREMCrauQ7Uq6k4fDZS7iQkmo6Sr4Sk9NwKikFdV+a6/S4RzrXwXM3NkR4CH+EErkSv6M8EO+yIW+TmpaOA6cvekXh4UxBv/e+WLEfX6zIfX3FZ3o2wNPdG3CkhagIWJQQUYGpKk4kpSA+Mdkl5xMIosqGolKZEAQHesYI4fakUvjfk53wf8v3438bjxb69e/9sRvv/ZFzLvyPQzqgTe2KubyCiDKxKCGiPCWcT8GRs5eK/PoKpUJQvUIYgjyk4CioFjUrYFL/VpjUv1Wex2w8dBZPTFuPw2cK9u9z5yc5lzZaMPR6NKxStsg5iXwNixIiAgCkpyt2xSfhclrhVw+uF1UGpUP968dJy6sqYMWL3XM8rqr4ZOlevDl/Z77nuPHdZQ79oT0b4t89G7gsI5G38a+fIkRkSVfFjmNJSC3EFgblSwXjqorhCOC8pzyJCB7vWh+Pd63v8Pih0xfR7e0lSE3Pe93Fd//YhXf/2GX1n+5eH8/e2MhtWYk8DYsSIj9y7NwlvP7aa/j03f8AADYeOpPnsYEiaFS1bIlfesmcbDpq1CiMHj26RN/bna6qGI4943o7PPbX3lO4//PVeb5m0qI9mLRoj9X/cUhHtKnNhRXJd7EoIfJhqWnp2HasYCuURoSHoGZEKd79VYI61otE3IQ+Vv/Mhcto9cbveR5/5yd/We32dSvi+8fa878X+RQWJUQ+5kJKKvaePF+gY5vXKM9fah4konSIQ5FyIjEZbcctzPXYVftOo86IrPVUtr1+E9dNIa/H/4OJfEDipSuIS7iQ73GNq5ZDlXJhVt8TCxI37XXplSqXC3MoUmauO4znZ2zM9dirR/5mtTeMvAEVwkPcno/I1ViUEHmpgoyIlAkNQp3I0h5ZfFDh3dWmJu5qUxMAcPbiZVzzeu6Xeuwf3/FGL4QFB5ZIPqLiYlFC5EUKMkekarkwVLYbDSHfVCE861JPerqi+ejfcCGXfYMavzofAFA3qjQWPde1JCMSFZp3rWhE5CYJCQl44YUX0KhRI5QqVQpVqlTBDTfcgFmzbBtFT506FSICEUFcXFyu50hLS8NXX32FW2+9FTVq1EBoaCgqVqyI9u3bY9y4cUhMzH/C6cWLF/H222+jc+fOiIyMRGhoKKpVq4auN/TChA8+x9aj53J9XbXyYWhRswJa1KyAK0kJeOKJJ1C3bl2EhYWhevXquO222/DHH38U6t8kOTkZH374IXr27ImqVasiJCQEUVFR6NatGz788EOkpKQ4ff3ff/+Nxx57DI0bN0aZMmUQGhqKGjVqoFWrVnj00Ucxffr0XM+R+e+c3503O3bswNChQ3HNNdegYsWKCAsLQ926ddGjRw9MnDgRhw4dKtTX680CAgRbX++FuAl9sH98b0SEB+c4Zt/JC4gePgfRw+fgx3WHDaQkKgBV5YcbPtq0aaNFsW3btiK9jopu06ZNWqVKFQWQ68egQYN0ypQpVn///v05zrFv3z5t3rx5nucAoNWqVdM1a9Y4zXHVVVc5PUfrth10+eb9uvHQGT2YcCHHOZYtW6blypXL8/WjR4/WUaNGWf28rFu3TmvVquU0S+PGjXXPnj25vv7999/XgIAAp68HoNu3b8/x2sznRo0aleu509LSdMSIERoYGOj03F26dMnz63PGl74HU66kae0Xf3X6celyqumY5AMAxKoLfnfy8g35tbNnz6JXr16Ij48HANx333146KGHULlyZezZsweTJk3C5MmTsXFj7pMLAeD48ePo1KkTjh07huDgYAwYMAA9evRAdHQ0UlJSsGTJErz33ns4duwYbr75Zqxfvx61atVyOMfRo0fRrVs3JCQkAAB6334X+txxDypWisLBuH34fupkbFi7GuvXrMQLg+7H8uXLERjoOE/g4MGDuOWWW5CYmIiAgAA8+uijuPvuu1GhQgVs3rwZEyZMwOjRoxETE+P032Tbtm3o0qULzp8/j9KlS2PIkCHo0KEDatWqhcTERMyfPx8ffvghduzYgV69eiE2Nhbly5e3Xr9p0yYMHToU6enpiI6OxhNPPIFWrVqhUqVKuHDhAnbv3o2lS5di9uzZhfpvlenJJ5/EJ598AgCoXLkynnjiCXTu3BkRERE4c+YM1q9fj1mzZnEeDYCQoADrEs++k+fR/Z2lOY7JvLwz5vZmeKB97RLNR5SDKyobfnCkxFs988wz1l/WEyZMyPF8amqq9u3b1+Ev8OwjJbfddpsC0OrVq+uOHTtyfZ99+/ZZozEPPPBAjufvuvtu6/wvvDZBNx464/Bx/lKK3nvvvdYxkyZNynmOu+6ynp8yZUqO5xMTE7Vly5YOX0t26enp1jFNmzbVY8eO5fr1rF27VsPDwxWAvvLKKw7PvfrqqwpAS5cunefrVVUvXLigly5dyvE4nIyUzJkzx3o+JiZGT506lef5Dx48mOdzzvjD9+Co2Vucjp6kp6ebjkheBi4aKTH+y9tXP1iUeL7k5GStUKGCAtAWLVrk+YP4+PHjGhYWlmtRsmXLFuvx6dOnO32/jz76SAFocHCwXrhgu/Ry6nyy/hG7XYOCghSAtmnX0SpEdhxLdHj92bNntWLFigpA69ev7/DcsWPHrMsZPXv2zDPD6tWrnRYlv/76q/Wcs0tNqqrDhg2zijF7jz32mALQVq1aOX19XpwVJZ06dVIAGhoaqnFxcUU6f3786Xsw8dJlp8XJicRk0xHJS7iqKOHlGy8VPXyO6Qglwn6NBleLjY3F2bNnAQD/+Mc/8hzur1KlCm666aZcLzf8/PPPAIDg4GD07dvX6ftdf/31AIArV65g9u/L0LRNewDA2r+WIzU1FQDQr/+DeW5uV758edxzzz349NNPsWfPHsTFxSE6OhoAsHjxYqSl2e68GDBgQJ4Z2rZti6ZNm2Lr1q25Pp/59dSuXRvXXnttvl/PW2+9haNHj+LgwYPWJanq1asDsF0GWrNmDdq2bev0PAV1+vRp/PWXbUXTO+64A7Vr81JDcZUNC7a+x5754W/8vOGow/PXjrVNjp4y4Fp0a1S5xPOR/+HdN+S3tmzZYrXbtGnj9Ni85mHExsYCsBUaoaGh1p0juX00b97cet3ho8es9p6d2632Pb17ON1tt3379lZ78+bNubbzKyacFQmZX8+BAwecfi0igltvvdV63fHjx612//79ERISgpSUFHTq1Am33HILPv74Y2zatAnphdj8L7sNGzbYhncBXHfddUU+D+XuvftaIW5CH/zyZOcczw2YshbRw+fg/T92G0hG/oQjJV7KnSMI/uLMmazN6KKiopwem9fzJ06cKNJ7J1+6BAAIDwlCUOpF6/HKlZ3/NVq1alWrffr06Vzb+Z2jSpUqeT5X1K/n4sWsr6FRo0aYPn06HnnkESQkJGDOnDmYM8c2shcREYGePXti4MCB6NWrV6He4+TJk1a7WrVqRcpJ+WteszziJvRBSmoaGr0y3+G5zF2Mb2tZHZP6tzKUkHwZixIiFH259cxLJjVr1sS8efMcnlNV7IpPyvV1LRvVRf2aFQqdI3OkwNnjRT0HkPX1tGvXDv/3f//n9Dz26tSp49Dv27cvunfvjhkzZmD+/PlYvnw5jh8/jjNnzmDGjBmYMWMGevfujZkzZ6JUqVIFfp9MvLPG/UKDAq0/frJfLv7fxqP438ajuK5BJL55pJ2JeOSjWJSQ34qIyNoCPr8RAvu/0u1FRkZazzds2BAhISFQVWw+YlvkrEG2XeajI0ujXJjjwlYVK1a02vHx8ShXrlyeOTJvXc7+usKcw9nXGhkZifj4eMTHx6NZs2Z5HlcQZcuWxcCBAzFw4EAAwO7du/Hrr7/iww8/xL59+zB37ly8/PLLmDhxYoHOl/lvDdhuoaaSk1mc3PHxn1h/8Kz1+PLdpxA9fA5uvLoKJj/k/FZzooLgnBLyW02bNrXa69atc3ps5lyL7Fq1sg1hp6SkYPHixdh0+KxVkNirX7kMWtSskKMgAeDwy3/16tVOc9g/bz9Hxb69du1ap+dw9nzm1xMXF4ddu3Y5PU9hNWjQAEOHDkVsbKx1CWn69OkFfn2rVq2sEZJly5a5NBsVzE+Pd0LchD64tWV1h8cXbItH9PA5eP2XbYaSka9gUUJ+KyYmxlr067vvvsvzskZ8fDx+++23XJ+7/fbbrfarr4/NcY7MYsTZlvLdunVDUJDt+S+//DLP4xITE61f4vXr17fuvMk8R+Zial999VWe54iNjXWY4Ovs6xk/fnyexxVHREQEWrduDQA4depUgV9XsWJFdOrUCQAwa9YsHDhwwC35KH8f9G+Va3Hy5Z/7ET18DqbH+s8S/+RaLErIb4WFheGhhx4CYFuF9M0338xxTHp6OgYPHozk5ORcz1G5ThNc1/0GALZbeyeOGQlVRd2o3IuR+Pj4HHM1qlWrhjvuuAOA7dbeTz/9NMf7qCqGDBlirfj69NNP5zhH5i3JCxYswDfffJPjHOfPn8egQYNy/Toy9evXzxq5mTp1KiZNmuT0+P379+P77793eGzWrFkOk4izO336tDUylX0uSn6GDx8OwDYydddddzlM8M3On/a+MSWzOOlYr5LD4y/M3ITo4XOwPZ/NI4myE2eT3qjoYmJiNK8hf2e2b9+OJk2auCER5eb06dNo2rSpdUtr//79HZaZf//99/HXX3+hbdu2WLNmDQDbpY2IytURl3ABAJBw8gTuv6U7jh89AgDWhnMtW7ZEeHg4zp49iy1btmDhwoWYN28emjdvnuNy0NGjR9GiRQskJCRARPDggw/i/vvvR1RUFPbu3YtJkyZhxYoVAIAOHTrkusx8XFwcWrRogaSkJAQEBGDQoEG4++67Ub58eWuZ+Z07dyImJsZ6/9y+/7dv344OHTrg3DnbZaiuXbviwQcfxNVXX42QkBAkJCRg06ZNmD9/PhYtWoR+/fph5syZ1uu7du2KNWvWoHfv3ujevTuaNGmCiIgInDt3Dhs3bsSHH36I3bttt5ZOmjQJTz31lMP7Z16iGTVqVK6b8g0ePBiTJ08GkLXM/HXXXWctM79hwwb89NNPCAoKwuLFi53+988NvweLrvmo35CUkprj8R1v9EJYcGAuryBfISLrVLX4E4tcsQIbP7iiqzfbsGGDRkVF5bmx28MPP6xffPGF1V+0fmeOZeDXbNltrTaa30e3bt1yzbFx40atWbOm09d27txZExIS8vxaFi9erGXLls3z9SNHjizQhnzbtm3Tq6++ukBfz4ABAxxe26VLlwK97umnn851Fd3M551tyPfss8+qiHBDPg+V1wqx5LvgohVdefmG/F7Lli2xbds2PPfcc2jQoAFCQ0MRGRmJbt26Ydq0aZgyZQoSE7OGocuUdbyzpUXNCri2aX2sWLECv/76Kx544AHUq1cPpUuXRlBQECpVqoS2bdviqaeewrx58/D777/nmqNFixbYsWMH3nrrLXTq1AkVK1ZEcHAwqlatit69e+O7777DsmXLHO60ya5r167YunUrhgwZgtq1ayMkJARVqlRBnz59MH/+fLz22msF+jdp0qQJNm3ahGnTpuGuu+5C7dq1UapUKQQHB6NKlSro3Lkzhg0bhmXLluWYBzN9+nR89913GDBgAFq1aoVq1aohODgY4eHhaNy4MR555BGsXLkS77//fpFu7Q0ICMA777yDDRs24PHHH0eTJk1QtmxZlCpVCnXr1kWPHj3w/vvv44cffij0uck14ib0wfbXc65DEz18Dt7+baeBROQtePnGTXj5xnccT0zG44MHYdYP36By1er4fa1tifbmNcpzvQwfxO9B1/r74Bn0+/ivHI+vGtEDVcuHGUhE7uCqyzccKSHKQ2paOjYdPouD8WewZMFcAECL1jFoWKUsWtSswIKEqABa1YpA3IQ+uLN1TYfH249f6Dd7eFHBsSghv7d3794cEz43HT6LbccSkZaWhrEvPYczp213vTw1+BFO2CMqgnfuaZnr9hjRw+fghzUHDSQiT8QVXcnvjR8/HitWrED//v3RsnUMLoeUQ0ryJezctgU/ff81tm/eCMC2Foj9JnREVHhxE/og7tQFdH17ifXY8J82Y/hPm7FvXG8EBHAE0p+xKCECsHPnzlxvP83Uvn17TJ8+nZdsiFwgOrI04ib0wcNT1mDJzqwtHOq+NBcPtK+FMbc3d/Jq8mUsSsjvPf70UJSOrI6Vy5fg6OGDOJOQgNTUK4isVAlt2rTBvffei/79+yMggFc7iVxp6oC2SE1LR/2Xszaz/HbVQXy76iB2jbkZIUH8nvM3vPvGTXj3jXfYdPhsjseurlYOQYH8Yeiv+D1oxmdL92L8vB0Ojz3YvjbeuL14G0NSyeDdN0TFcPFyao6CJDwkCC1qVmBBQmTA4C71sH98b4fHvll1ANHD5+S5LxX5Hv70Jb+z6fBZ7Dlx3uGxq6uXQ/3KZQwlIiLAtsVA3IQ+GHnL1Q6P1xkxFzPXHTaUikoSixLyG2np6TlGR4ICAmyjI5wvQuQxBnauk2PU5PkZG7muiR/gT2LyC0fPXsLWo447ljaqUhZXVy+XxyuIyKTMUZM7WtdweDx6+BycSMx9127yfixKPBCvn7rWpsNncep8isNjLWpWQCgXQaNs+L3neSbecw02jb7R4bG24xbi2ekbDCUid2JR4mECAgKQnp5uOoZPuJKW83JNjQql0KJmBUOJyNOlp6fz1m8PVC4sOMdqsD+tP8LLOT6I330eplSpUrhw4YLpGF7v2LlL2H7M8XJNsxrlUalMqKFE5A2Sk5MRGsr/RzxV3IQ++KB/K4fHoofPwYkkXs7xFSxKPEzZsmWRlJRkOoZX23T4LE4m5bxcE8DVWCkf58+fR+nSpU3HICdubVkde8be7PBY27EL8cHC3YYSkSuxKPEw5cqVw8WLF3HmzBnTUbyOqua8XBPByzVUMMnJyTh37hzKlePkZ08XFBiQ43LOO7/v4uUcH8CixMMEBgaidu3aOHXqFI4cOYLERNtOtZyA59zFy6nYfOScw2NNq5dHpdIciqe8qSquXLmCU6dO4eDBg6hatSov33iR3NY04WJr3o3LzLtJUZeZz5SWlobExEQkJSXh0qVLnPzqxJmLl3EhJc3hsZoRpQylIW8TGBiIMmXKICIiAqVK8f8bb3T6wmW0fuN3h8dWjeiBquXDDCXyP65aZp5FiZsUtyihgsk+XNuiZnn878nOhtIQkUnZfx78587muPfaWobS+BfufUN+L/sPoE8faM2ChMiPZZ9n8uKPm3HrBysMpaGiYFFCXudKWnqOgmTjyBvRq1k1Q4mIyFPETeiDp7rXt/qbj5zjBFgvwqKEvMrBhIto8PI8h8fiJvRB+fBgQ4mIyNM8d2Mj/PqU46gpCxPvwKKEvMb8Lcdx/VuLHR7LPlxLRATYFkvc8UYvh8eih89BWjrnUXoyFiX5EJFxIqIZH8+bzuOvxs/djn99u87qlwsLYkFCRE6FBQfm+DlR76W5OHfxiqFElB8WJU6IyLUAXgDA0tqgvh+uwGfL9ln9B9vXxqbRNxlMRETeJHth0vL1BdgVz5WzPRGLkjyISCiAqQDiAcw2m8Z/RQ+fg42HsxZFm9S/Fd64vZnBRETkjbIXJje+uwyLd5wwlIbywqIkb68DuBrAvwCcy+dYcoPsE9PmPN0Zt7WsbigNEXm7uAl90KhKWas/YOpafL0yzlgeyolFSS5EpB2A5wBMU9VfTOfxR9kLkrUv90TT6uUNpSEiX/Hb0Otxb8xVVn/k7K2YMG+HwURkj0VJNiISBuArAKcB/NtwHL+UvSDZ8UYvRJXlfiRE5Br/uasFXu7dxOp/unQvnp+x0WAiysSiJKexABoBeEpVT5kO42+yFyR7x/VGWHCgoTRE5Kseu74uJvVvZfVnrjuMwd9waxDTWJTYEZGOAJ4B8LOq/td0Hn+TvSDZP743AgPEUBoi8nW3tayOqQOutfq/bY3HgClrDCYiFiUZRKQUgCkAEgE8XsRzDBKRWBGJPXnypEvz+brsBUnchD4QYUFCRO7VtVFlTB/cweov3nkSA6euNZjIv7EoyTIOQEMAz6rqsaKcQFUnq2qMqsZERUW5Np0Py60gISIqKW3rVMSPQzpa/UU7TuDp7/82mMh/sSjJ0g9AOoB/isgS+w8AmWsVD8l47P+MpfQxLEiIyBO0qR3hUJj8b+NRjP7fVoOJ/FOQ6QAeJgBAFyfP1834qFAycXxbbnNIiIhMaVM7AtMebYf7/281AGDqX3GoUi4MQ7rWM5zMf3CkJIOqRquq5PYB2y3CADAs47FrTGb1BbkVJJxDQkSmdawfickPtrH6/5m/A7M3HDGYyL+wKKESV2eEY0GyZ+zNLEiIyGPc2LQqxvVrbvX//cMGrDtw2mAi/8GihErU/Z+vgtptb7j1tZsQFMj/DYnIs9zfrhYet7tsc+cnK3Hs3CWDifwDfxtQiZm0cDf+2ptg9VeN6IHSoZzWRESe6YVejdGtUdadlB3GL0LylTSDiXwfi5ICUNWHM+aSvG06i7davvskJv6+y+r//EQnVC0fZjAREVH+pgxoi+DArMvLjV+dbzCN72NRQm4Xn5iMB7/IWiVxwh3Ncc1VvIGJiLzD7rGOdwZmn6hPrsOihNwqLV3RbtxCq39ry+q4r20tg4mIiAov+xpKLEzcg0UJuVW9l+Y69D+w2wCLiMibsDBxPxYl5DZcrZWIfM2+cY6Xcl77hau+uhKLEnILFiRE5IsCAgSrRvSw+lP+jMOGQ2cNJvItLErI5SYt3O3Q5/LxRORLqpYPw/v3ZS3sfftHfyIllbcKuwKLEnKpuFMXHG79XTmiO1drJSKf0/eaGmhVK+suwkav8FZhV2BRQi6jquj69hKrP+b2ZqhWvpS5QEREbjTr8U4OfU58LT4WJeQydUZk3WlTNjQID7SvbTANEZH7ZZ8vd9/klYaS+AYWJeQS2f9C2PzaTYaSEBGVrD1jb7baq/adxroDZwym8W4sSqjYvl9z0KHPO22IyJ8EBQZg+uAOVv/OT/6C2u88SgXGooSK5dLlNIz4abPVX/1SDydHExH5prZ1KqJlzfJW3/5yNhUcixIqliYjs2acD7q+LqqU4yZ7ROSfZj/Z2aHf4GUWJoXFooSKLPs8kpd6NzGUhIjIM9hfvr6Splix+5TBNN6HRQkVyU/rDzv0OY+EiMhm4XNdrPYDX6zm/JJCYFFChZaernh2+karH/tKT4NpiIg8S72oMmhbp6LV5/ySgmNRQoVW127n37va1ERkmVCDaYiIPI/93TgAMPibWENJvAuLEiqUez9zXBjo7btbGkpCROTZ7Pf9+m1rPM5cuGwwjXdgUUIFlnA+Bav3n7b6nEdCRJQ3EcHH/2ht9Vu98bvBNN6BRQkVWJsxf1jtLx+OMZiEiMg79G5ezaHP24SdY1FCBdJ0pOMOmN0bVzGUhIjIu2S/TfhgwkWDaTwbixLK17Fzl3DhcprV52UbIqLCmfN01sJq17+12GASz8aihPLVYfwiqz37iU5OjiQiotw0rV7eoZ998UmyYVFCTnUcv9Ch3/KqCoaSEBF5t+yjzPGJyYaSeC4WJZSnpOQrOHou65uGl22IiIrn16eyLuO0G7fQyZH+iUUJ5an56AVW+7tH2xlMQkTkG5rVcLyM03/yKkNJPBOLEsrVG79uc+h3qh9pKAkRkW+xHw6SEicAACAASURBVHVeuS8BaencGycTixLKQVXxxYr9Vp+XbYiIXOsdu9Ww673EtUsysSihHOw3jxrcpa7BJEREvunONjUd+vM2HzOUxLOwKCEHe06cd+iPuLmJoSRERL5t99ibrfaQ79YbTOI5WJSQg54Tl1rtda/0NJiEiMi3BQcGoGZEKavf6vUFTo72DyxKyDJq9haHfqUyoYaSEBH5hxUvdrfaZy5e8ftJryxKyPLVygNWm5NbiYhKxuQH21htf5/0yqKEADguefxA+1oGkxAR+Zcbm1Z16O+KTzKUxDwWJYTLqekO/TG3NzeUhIjIP20efaPVvvHdZQaTmMWihNDwlXlW+8chHQwmISLyT2XDgh3636yMM5LDNBYlfu7wmYsO/Ta1KxpKQkTk3+zn8r06e6vBJOawKPFznf+z2GpvHHWjkyOJiMjdoiuFW+1BX8caTGIGixI/tnJvgkO/fKngPI4kIqKSsGRYN6u9YFu8wSRmsCjxY/0/z9qdct+43gaTEBFRpgfb17bard/43WCSkseixE/N3nDEoR8QIIaSEBGRvTdub2a1T1+4bDBJyWNR4qf+/cMGq82F0oiIPMuoW6+22vbrSPk6FiV+aHrsIdMRiIjIiQGd6piOYASLEj/0wsxNVpujJEREnunNO1tYbX8ZLWFR4md+3XTUdAQiIiqAe669ynSEEseixM88Oe1vq81REiIiz/Z636ZWu+VrCwwmKRksSvzI2rjTpiMQEVEhPNQh2mqfu3TFXJASwqLEj9z96UqrvX881yUhIvIGg6+va7Uf8/FVXlmU+IkTickOfRGuS0JE5A1G9G5itX/38VVeWZT4ibbjFlrtnWN6GUxCRESF1bR6Oav93eoDBpO4F4sSP5Cerg790KBAQ0mIiKgo5jx9ndV+edYWg0nci0WJH6j70lyrvfj5ruaCEBGRS8SdumA6gluwKPEzdSJLm45ARERFYL9xate3l5gL4kYsSuyISLCI9BCRd0RklYgcE5HLInJERGaKSFfTGQtr6H+z9rh59ZarnRxJRESeLPvGqaqax5Hei0WJoy4A/gDwLIDaANYBmAXgNIA7ASwWkdfNxSu8WX9n7Qb8SGf/3EuBiMhX/PFsF6tdZ8RcJ0d6JxYljtIB/AjgelWtpqq3qOq9qtocwH0A0gC8KiLdjKYsoD0nkkxHICIiF6pfuYzpCG7FosSOqi5S1btUdXkuz/0XwNSM7gMlGqyIek5cZrW5WBoRkW+4q01Nq/3ftQcNJnE9FiWFk7lxTE2nR3kgLpZGROQb3r67pdV+8cfNBpO4HouSwmmQ8fmY0RQF0HPiUqv945COBpMQEZE7paalm47gMixKCkhEqgJ4OKP7o8EoBbLnxHmr3aZ2hMEkRETkahtG3mC16788z2AS12JRUgAiEgTgWwDlASxU1V/yOG6QiMSKSOzJkydLNKO9Y+cuGXtvIiJyvwrhIaYjuAWLkoL5FEAPAIfgZJKrqk5W1RhVjYmKiiqxcNl1GL/IanOCKxGRb+rWKOv3zNq40waTuA6LknyIyPsAHgFwHEAPVT1uOFKhcIIrEZFvmjKgrdW++9OVBpO4DosSJ0TkHQBPAzgJW0Gy23CkfE1fe8hqj7i5scEkREREhcOiJA8i8iZsK7smALhBVbcZjlQgL/y4yWoP7lLPYBIiInK3L/4ZY7XHz91uMIlrsCjJhYhMADAMwBnYCpKNhiMRERHl0KNJFav92bJ9BpO4BouSbETkDQAvAjgLW0Hydz4v8RgvzcpaROfHIR0MJiEiIio8FiV2ROQ2AK9kdPcAeEpEpubyMdxgzDxNW5213HCb2hUNJiEiopJiv0nfwKlrDSYpviDTATyM/W/ymIyP3CwFMMH9cYiIiJyz36Rv0Y4TBpMUH4sSO6o6FVmb7nmVF2ZmTXuZ/8x1BpMQEREVDS/f+IjpsYetduOq5QwmISKikrb4+a5W+6nvvWYqZA4sSoiIiLxcncjSVvuXjUcNJikeFiU+4JuVcVZ72qPtjOUgIiIqDhYlPuDV2Vutdsf6kQaTEBGRKV8NzFp2fkbsISdHei4WJURERD6gS8OsDfqGzdzk5EjPxaLEyx09e8lqP9ShtsEkRERExcOixMt1nLDIar/et5nBJEREZFpPu2Xnz128YjBJ0bAoISIi8hGfP9TGard8fYHBJEXDooSIiMhHiIjpCMXCosSLfbp0r9VeNqybwSRERETFx6LEi02Yt8Nq16oUbjAJERF5immPZa1XtWzXSYNJCo9FCRERkQ/pWC9rvaqHvlxjMEnhsSjxUqpqtQO8+xIiERERABYlXmuk3Squm0bfZDAJERGRa7Ao8VLfrDpgtcuEBhlMQkREnubXpzpb7aVeNK+ERQkREZGPaVajvNX+pxfNK2FRQkRERB6BRYkX+n1bvNVeOqyruSBEREQuxKLECz32dazVrl2ptMEkRETkqUbc3NhqJyV7xz44LEqIiIh80OAu9ax2r/eWG0xScCxKiIiIfNyRs5dMRygQFiVexn7RNPstqomIiLwdixIv8+4fu632Jw+0NpiEiIjItViUeJlJC7OKkuBA/ucjIqK8ffpAG6t99uJlg0kKhr/ViIiIfFSvZlWt9r2frTKYpGBYlBAREfmBnfFJpiPki0WJlyrL/W6IiMjHsCjxIvtOnrfas5/sZDAJERGR67Eo8SL9P8+6Hlg3qozBJERE5C0aVvGe3xcsSrxIfGKK6QhERORlvh7YzmrvPO7Z80pYlBAREfmwquXDrPZDX642mCR/LEqIiIj8hKePuLMoISIiIo/g8vtKRUQANABQGUAFAGUBXAZwHsARAPtU9aKr39fX2e9580qfJgaTEBERuYdLihIR6QTgTgAdATQDUCqf4/cCWAHgNwC/sEjJ3/wtx632wE51DCYhIiJyjyIXJSJSBsBQAP8CUBWAAFAA8QC2Zny+BCAZQAhshUolAFcBiAZQH8A/AVwUkf8C+I+q7gbl6t//3WC1AwLEYBIiIvI2PRpXxsIdJ0zHyFehi5KMyzPPAXgRQEUAsQA+gW3kY72qJhbgHKEArgHQHsAtAB4E8LCIfAXgeVU9U9hcvu5yarrpCERE5KVe69vUKkpOX7iMiqVDDCfKXaEmuopIdQBLAIwC8DWAhqraTlXHqOqSghQkAKCqKaq6WlXfV9UbAFQHMBJALwDbRaR7ob4KIiIiylPNiHCr/faCnQaTOFfgokREWgPYCCAOtmLkOVXd64oQqpqgquNgmyD7GYDZIjLQFecmIiKiLNNWHzQdIU8FunwjIo0A/BfAv1T1R3eFyZjwOkpEZgP4XkQuqer37no/IiIi8hwFHSn5BsDd7ixI7KnqegCdATwlIvVK4j2JiEpEynlg+j9tn4nIQUGLkm6quiH/w1xHVU8C6A7gQEm+r6e7NjrCdAQiKo79S4FtPwP7l5lOQuRxClSUqOoFdwfJ432TVTXVxHt7EvuF0/7VhQNHRF5t+y+On4nI4vIVXQFARAJgW0itOoBjAP5WVY5VFtHWo1k3NXVrVNlgEiIqFlVg12+29u7fbH3hukNEmVy+942I3AXgEIClAL6H7RbioyLygYiUc/X7+YNpa7JmSnPhNCIvdnIHkJpsa1+5BJz03FsziUxwaVEiIo8AmA6gGoDTABJhW+m1DIDHAWwUkZaufE9/8MMaz719i4gKYdcCID3N1tY0YPcCs3nIr9zRuobpCPlyWVEiInUAfAjgOwD1VTVKVSMAlAdwL2z73NQG8LuINHXV+/qDdM3/GCLyAttmAWkZW8enpgBbZ5nNQ36lXysfKkpEJDyfQ/4NYJqqPqiq+zIfVNUkVZ2hqr0B9ACQCuBHEXG6aR8RkU9JSQLitzo+Fr+FtwZTiWlft5LVjk9MNpgkb4UZKYkVkY0iMklE7hSRyGzP9wTwsrMTqOpiAO0ABAN4vXBRiYi82P5lQFCo42NBobw1mEpMcGDWr/yF2z1zc74CFyWqejVsm+etAXATgBUisllEPhSRuwGIqh4vwHkOAegH4FERiSpibiIi77Ltl5yjIinnge3/M5OH/NrSXV5elAC2gkJVv1XVQaraGLbRkWUAugBoICLbMoqUfiKS5ypfqroJtlViHyhOeCIir6BquwUY2SeIZdwirJw4RiVr1b7TpiPkqlgTXVU1XlWnq+qTADYA6ArgTwB9YLvcs1ZEJojIjbnMIfko4zgiIt9mfytwdqnJvDWYSty5S1dMR8iVKxdPSwUQnLGB3vcAICJ1YRtNeQTAlyKyB8DCjI/VALhmOhH5PvtbgbPLvDW4cuOSzUTkgVy5TslvsN2BY1HVfao6WVXvVdWaGc9fAPAKgMMAGorIUBG5xoU5ik1E7heR5SJyTkTOi0isiDyRsVItEVHh2N8KnB1vDSayuHKk5FsAa0XkE1Xdn9sBqroRwEYAE0UkGMBWAKUBvC0ijQDEAlgMYKGqbs3tHO4mIh/BttBbMmwjOldgu5X5QwA9RORuVc3jTx4i8ln/faDo+9UEhjh/Pn4LMLp80c7d5Fbg3m+L9loiD+Oyv/xVdS+AuQBmFGBNEwAIB5CmqmNUtSeAugDeAVAOwBhX5SoMEbkTtoLkOIAWqnqLqvYD0ADAdtjuGnrSRDYiMqzna0DV5kBwQX68ZZN2uXjP5yY4HKjawpaLyEe4+nLEUAA1ASzJmE/izHMAlmd2VPWKqq7IKFL6uThXQY3I+Pyiqu62yxYPYEhGdzgv4xD5oUr1gEFLgW4vAUGlAFM/BiTA9v7dXrblqcSdw8l3uPS7SlVPArgfQDMAm0VkYvYl5UUkXEReA/ACgI9d+f7FISI1AbQBcBnAjOzPq+pSAEcAVAXQvmTTEZFHCAgEOj4FDPkTqNK0aKMmxREcDlRpZnv/jk8CAfz7iHyLy/+PVtVFAO6A7W6cfwPYJCLxIrJSRNbAdmnkFQDjVHWDq9+/GFplfN6qqpfyOGZttmOJyB9ZoyYvl8yoCUdHyE+45TtJVecDaAnbHTkAEAXb8vIxsK0e9Kyqetoy83UyPh9wckzmdr11nBxDRP4gINA2WjHkT9vohbtGTYLDgSrNOTpCfsGVd984UNU4AL1FpD6ATgCqwPYLf76qnnPX+xZDmYzPF5wck7lGdNncnhSRQQAGAUCtWrVcl4yIPFfmqMmqj4FFY2y3/mp68c8rAUBgqG10pP3jLEbIL7itKMmkqnsA7HH3+7iAZHwu8nrPqjoZwGQAiImJ4brRRP4iIMA2itHoZmDGw0DCHuDKxaKfLzgcqNQAuHsKL9WQX2HpnSUp43MZJ8dkPpfk5Bgi8leZoyadhwJBYUU7R1CY7fWDlrAgIb/DoiRLXMbn2k6OuSrbsUREjgICgDJVij75VQKBslV5uYb8Ur7/14tItIjcVRJhcnnvF0rw7f7O+Nw0l80DM12b7Vgiopxivyj65ZsrF4C1X7g2D5GXKEgpfgTAv0TkcXeHySQiwSLyDUpwwz5VPQRgPYAQAHfnkqkLbAvDHQewsqRyEZGXOXcYOLGjeOc4sR04d8Q1eYi8SL5FiapeAXAngEdF5BsRyfXOE1fJ2ANnNYBUVR2R3/EuNj7j838y7hrKzFQZWQu9TVB1xdR6IvJJm6YDIvkf54yI7TxEfqZAFy0zbuG9HkBFADtF5F8i4tI7d0QkUkQmwrZh3zxVHeDK8xeEqs4E8Alsq7ZuFpFfROQnALsBXA3gZ9g25iMiyt26KUBqcvHOkZpsOw+RnynwTCpVPQ/gFgATAbwN4JCIjBORa4r65iISIiI3i8jXAA7BdtnkXlV9uajnLC5VfRzAP2C7lNMFwE2w3dL8JIA7uUMwEeXp+Bbgwknnx2SuzprfRNgLJ4B4I5ulExlTqNEOVVUAb4vIjwBGAXgewIsichLAXwA2wTaqcBTAKQDJAK4ACAVQCrYF1K4C0Ai2fWauhW234CTYip3/qGpi8b+s4lHVaQCmmc5BRF5mw3dAWmrezweHA5XqAze8Bvw+yvl6JmmptvPdNM49WcmvVS9fxFvW3axIl2BUdT+Ah0VkJIABAPoDuD3joyCLhknGcSsBfA/gWw9d5ZWIqGDS04GN3wPpV3I+l9vqrHW6ZKwCOzb3VWDTrwAbvgduGMPbg8nlOtaPNB0hV8WaF6KqBwG8BuA1EakF4DoA1wCoB6A6bMuxhwK4CCARtmXm9wCIBbBCVc8U5/2JiDzGgT+BtFwKkszRkbunOi6GlrnjcKPewIx/Agl7c46apF0GDv4FRHd2a3TyP10bRZmOkCuXTVbNKFC+y/ggIvIv678GLtttnVXQvWsc9s7JNmpy5aLtvCxKyAUupGRdWry+oWcWJRwTJCIqrivJwI5fYF29Dg637Rxc0J19M0dNhvwJVGmateOwpgPbfwFSU9wan/zDwh0nrHa5sGCDSfLmsqJERN4QkbkiMkJE2ooUdY1lyi6qbKjpCETkzO4FtuXhM++s6faybfSjsHvXZI6adHsp6w4dCbCdn6iYfv7b8xfkc2XhMAi222fHwDaB9bSIzBaRf4tIcxe+j9+5N+aq/A8iInNiv7RduinM6Eheso+aXL5gOz9RMS2yGynxVK4sSu4C8BNstwILgHIAboXtVt8NIhIvIj+IyGMiwq0vC6F/u1pWO/kKl0kh8jgh4cCNY4o2OpKXzFGTG9/IupxD5ONcOdF1OYDlACAiTQF0zfjoAiASQBSAe5Cxr4yIHAKwEMA0VV3oqhy+qEaFrP0BZ284gnuvreXkaCIqcfe5aVmjzFGTjk+55/xEHsYt8z5UdauqfqSqd6tqZQAtADwN20hKAmwjKbVgW+NkgYisE5FW7sjiaz5Zstd0BCIiIrcokcmoqrpFVT9U1bsAVAPwHoALAOYDuAygFYBVInJXSeTxZnEJRdwOnYiIyMOV+B0yqpqmqs8C+BbAHAA1AAyHrUj5XkR6lnQmIiIiMs/kbbsvARioqqdV9U0ADQAsAjBFRMoZzEVEREQGuHKdknAR+UBENonIjyJys7PjM5aYD7HrJwDoA2AjgIGuykVERERZht3UyHSEPLlypGQSgCcANAPQD8CvIrJfRF4RkbrZDxaR2rDtjWNR1dSMc9zjwlxERER+7e+DWVvNPXpdHYNJnHPZLcEAboBtrZJwAL0B9AVQG1kb9u0GsArAUQAVAdwB4LfsJ1HVA1wNNqcGlctg94nzpmMQEZEXGjZzk9UODQo0mMQ5V/7yP6OqP6nqt6p6P2y7BD8FWyECAA0BPATgRdhWf02BbV4JROQ2EXlLRKpnHBvmwlw+4Z17Wlpt+02ViIiI8rPHS/6odWVRMldEhmd2VPVcxlolnWBbk2QwbJd4JgMYCqClqh7KOPwTAM8C+FlE2sI2mkJ2WtSsYLXHzt1uMAkREZF7uLIoGQmgjojMyygsLKp6RFU/V9VnVHWIqr6vqqftDlkG24JqbQDMRMYICuVu2uqDpiMQERG5nCuXmU8FMFhEOsM2h6Q2gBhVLchqXw8C+BVAFQCzVHW/q3IRERGRd3DlRFcAgKquAHCziFQqYEGSWdB85+osRERElOW6BpGmIzjltrtcMtYdISIiIoPW290OPOk+z95mjrfeehFPr3CJiMjz3PHxX1Y7onSIkyPNY1HiRT76R2urvWTnCYNJiIiIXI9FiRcpFxZstR+estZgEiIiItdjUUJEREQegUUJERGRj1JVq31H6xoGkxQMixIiIiIf9e7vu6z2m3e2MJikYFiUeJnPHmxjtU+dTzGYhIiIPN2kRXusdlCg5//K9/yE5OCmplWt9o3vLjOYhIiIyLVYlHix0xcum45ARETkMixKiIiIfFBaetYk136tPH+SK8CihIiIyCcN/ibWar9zd0uDSQqORYkXWvhcF6u9eAdXdiUiopz+2J71+yEgQAwmKTgWJV6oXlQZqz1gKld2JSIi38CihIiIiDwCixIiIiIfMyP2kNVe/kI3g0kKh0WJl3r/vmus9p4T5w0mISIiTzNs5iarfVXFcINJCodFiZfqe03W7V09Jy41mISIiMg1WJQQERGRR2BRQkRE5EM+WbLXav85vLvBJIXHosSLTXusndVetCPeYBIiIvIU/5m/w2rXqFDKYJLCY1HixTrWi7TaA6fGOjmSiIjI87EoISIi8hHJV9KsduOqZQ0mKRoWJV4uODBr6WD7zZeIiMj/NH51vtWe+/R1BpMUDYsSL7fltZusdscJCw0mISIiT+It+93YY1Hi5UKDAq12fGKKwSRERETFw6KEiIjIB7yzYKfVnv+M9126AViU+ITZT3Sy2pMW7jaYhIiITPlg0R6r3bhqOYNJio5FiQ9oeVUFqz3x910GkxARERUdixIiIiIvN2fTMav91l0tDCYpHhYlPuKD/q2s9nS7LauJiMj3PTFtvdW+O+Yqg0mKh0WJj7i1ZXWr/YLdltVERETegkVJBhFpJCJDRWSeiOwRkWQROSciK0XkGREJMZ2RiIgou9kbjljtMbc3M5ik+FiUZFkIYCKArgCOAZgFYB2AawC8C2CViFQ0lq4APv5Ha6v96dK9To4kIiJf8e8fNljtB9rXNpik+FiUZNkJ4BEAUap6nar2V9XuAJoA2AqgFWzFicfq3bya1Z4wb4eTI4mIiDwPi5IMqtpDVb9U1fPZHo8D8K+M7j28jENERJ5i9P+2Wu2vBrY1mMQ1WJQUzN8Zn8MAVDIZJD8Lhl5vtft+9KfBJERE5G5T/4qz2l0aRpkL4iIsSgqmQcbnywBOmwySn4ZVsraq3njorMEkRETkTuk+uDM8i5KCGZ7x+VdV9apd7xKTr5iOQEREblD3pblW++9XbzCYxHVYlORDRB4GcC+AiwBeyufYQSISKyKxJ0+eLIl4udoz9mar3WL0AmM5iIioZESU9o3pjkGmA7iCiLwJ4LYivLSHqh7J60kR6QHgMwAKYLCq7szrWABQ1ckAJgNATEyMsXG1oEDWmkREvmxtXNZMgmvs9j/zdj5RlACoDqBREV4XnNcTItIZwGwAIQCeVtVvi5jNiOduaIh3Mjbn+3DRbjzZvUE+ryAiIm9x96crrfbPdjvFezuf+JNaVR9QVSnCR1xu5xORjgDmAigN4EVV/aAkvx5XeKpHVhHy9gLuHExE5CtUfW+CayafKEpcSUTaA5gHoCyAV1T1TcORXIITXomIfEOdEVkTXO2XgfAFLErsiEhbAL8BKAdgtKqONRypWHa80ctqc8IrEZHvsV8GwhewKMkgIm0ALICtIHlDVV8zHKnYwoIDTUcgIiIX+mHNQavds0llg0ncg0VJlt8BlAdwFkAtEZmax0ek4ZyF8tH9WZv03fLBcoNJiIiouIb/tNlq/98/rzWYxD185e4bV4jI+FwBwD+dHDcawCm3p3GRPi2q4YlptvaWI4lmwxARUZGdSEo2HcHtOFKSobh37Hiy5jXKW+2Z6w4bTEJEREXVduxCq71rzM1OjvReLEr8wC9Pdbbaz8/YaDAJEREVRfZ9bkKCfPPXt29+VeTU3pPnTUcgIqJCsN/n5o9nfes2YHssSvzEvnG9rXaPd5YaTEJERMVRv7Jv3QZsj0WJnwgIEIf+uUtcTI2IyBtcPXK+1X7v3msMJnE/FiV+ZOOoG612y9e4mBoRkTe4eDnNat/eqobBJO7HosSPlC/luP9g8pW0PI4kIiJP0PfDFVb78a71DCYpGSxK/Myfw7tb7cavzndyJBERmbbx8Dmr/UKvxgaTlAwWJX6mRoVSDv2UVI6WEBF5ons+W2m1+7etZTBJyWFR4oeWPN/Vajd6haMlRESeaM3+01Z7/B3NDSYpOSxK/FB0ZGmH/oWUVENJiIgoNze9u8xq39/OP0ZJABYlfmvFi92sdtNRvxlMQkRE2e2MT7La4/r5xygJwKLEb9WMCHfoHzt3yVASIiKyFz18jtV+ukcDg0lKHosSP7ZxZNa6JR3GLzKYhIiIAEDVcY+bZ29oaCiJGSxK/Fj5cMd1S1buTTCUhIiIAKDOiKw9bj7+R2uDScxgUeLn9trtidP/81UGkxAR+bfz2W466N28mqEk5rAo8XOB2fbEGTtnm6EkRET+rZndTQeLnutiMIk5LEoIcRP6WO3Pl+83mISIyD+t3ud4+bxuVBlDScxiUUIAgL7XVLfa9jO/iYjI/e6dnHX5fNeYmw0mMYtFCQEA3r+vlUOfm/UREZWMx79b59APCfLfX83++5VTDv8d1N5qc7M+IqKSMXfzcattfzndH7EoIUu7upUc+r9tPZ7HkURE5Ar2l8ufv9G/1iTJDYsScrB7bNa1zMHfrHNyJBERFceeE0kO/Se7+9fqrblhUUIOggMDUDOilNXnpFciIvfoOTFr073YV3oaTOI5WJRQDite7O7QP3fxiqEkRES+6faP/nToR5YJNZTEs7AooVz9OKSj1W75+gKDSYiIfMvl1HRsOHTW6vv75FZ7LEooV21qRzj0n5+x0VASIiLf0vCVeVb7g/6tnBzpf1iUUJ7sq/eZ6w4jLV2dHE1ERPl5Z8FOh/6tLavncaR/YlFCTr16y9VWu95Lc50cSUREzqgqPli0x+rvH9/bydH+iUUJOfVI5zoO/fFztxtKQkTk3eqMyPrD7vGu9SAiTo72TyxKKF/21fxny/YhNS3dYBoiIu+T/bLNC70aG0ri2ViUUL5EBMNuamT16788z8nRRERkLzUt3eGyzb5xvGyTFxYlVCBPdKvv0H/0q1hDSYiIvIv9H3JDezZEQAAv2+SFRQkVmP1lnD+2x+P0hcsG0xAReb77Jq906P+7J5eSd4ZFCRWYiOCj+1tb/dZv/G4wDRGRZzt+Lhmr9p22+rzbJn8sSqhQ+rSo5tDn3jhERLlrP36h1Z7y8LW826YAWJRQoWVfEvnXTUcNJSEi8kzZ/2Dr1riyoSTehUUJFcmal3tY7Sen/c3VXomIMny0eI9Dn3vbFByLxCamoQAAIABJREFUEiqSymXD0K9VDavP1V6JiICk5Ct467esNUm2vnaTwTTeh0UJFdm7917j0Of8EiLyd81HZ+2qPuymRigdGmQwjfdhUULFkn1Y8oc1Bw0lISIyK/sfZtnXd6L8sSihYls5orvVHv7TZiQlXzGYhoio5A397waHPueRFA2LEiq2auVLYXCXulbffviSiMjXbT16DrP+PmL1d47pZTCNd2NRQi4x4uYmDn3OLyEif5CerugzaYXV/+zBNggNCjSYyLuxKCGXyT5cGTOGK74SkW+ra3fnYZ3I0ripaVWDabwfixJyKfvdL0+dv4yvV8YZy0JE5E7ZR4QXP9/VTBAfwqKEXCogQLDwuS5Wf+TsrTiYcNFgIiIi1+v13jKHPie2ugaLEnK5elFl8PyNDa3+9W8t5oqvROQzvlkZhx3Hk6z+rjE3mwvjY1iUkFs82b0BytgtGsQVX4nIF+w5kYRXZ2+1+nOe7oyQIP4qdRX+S5LbbMm2vDLvyCEib5Z8JQ09J2Zdthl2UyM0rV7eYCLfw6KE3Cr7dVYWJkTkrRq/Ot9qX1WxFFdsdQMWJeR2+8f3duhfO/YPQ0mIiIom+x9Uy1/onseRVBwsSsjtRATbX89a4fBkUgqen7HRYCIiooLLXpDwThv3YVFCJaJUSCCWDutq9WeuO4wpf+43F4iIqABYkJQsFiVOiEgzEUkRERWRLabzeLvalUpj6oBrrf5rv2zDwu3xBhMREeUte0GS/VI0uR6LkjyISBCAqQCCDUfxKV0bVcart1xt9R/5KhYbD501mIiIKKfmo35z6O8aczNExFAa/8GiJG8vAWgD4GPTQXzNI53r4P52tax+34/+xN6T5w0mIiLK0mfSciSlpFr9jSNv5FokJYT/yrkQkRYAXgHwE4CZhuP4pHH9mqNjvUpWv8c7S3Hk7CWDiYiIgAe/WI2tRxOt/ooXu6F8OAfMSwqLkmxEJBjAVwDOA3jccByfNu2x9qhVMdzqd5qwCMfPJRtMRET+7F/frMPy3aes/vxnrkPNiHAnryBXY1GS0ysArgHwjKpyFqabLXuhG0qHBFr99uMXIj6RhQkRlazHv1uH+VuPW/1Zj3dE46rlDCbyTyxK7IhIK9jmksxX1a9N5/EXW1/vhcCArAlk7cYtxOEz3FmYiErGgClrMHdzVkHyw6D2aFUrwmAi/8WiJIOIhMB22eYSgEFFPMcgEYkVkdiTJ0+6NJ+v2zuut0Nh0vk/izn5lYjc7o6P/8TinVk/r38Y1B7t61Zy8gpyJ58oSkTkTRHZUYSPGnanGQmgOYAXVPVQUXKo6mRVjVHVmKioKJd8bf5k77jeKBuWtbNwj3eW8nZhInKbduP+wPqDWT9jfhzSkQWJYUH5H+IVqgNoVITXBQOAiLQB8CKAJQA+c10sKqzNo29Cl7cW40CC7fJN34/+xOcPxeCGq6sYTkZEviT7wmi/PtUZzWpwx1/TfGKkRFUfUFUpwkdcxiluha1AqwJgsYgsyfwA8F7GMXXsHufWkG60dFg3dLD7a+Wxr2Px5QouSU9ErpG9IFk6rCsLEg/hKyMlrtIk4yM34QC6ZLTLlEwc//X9oPZ4edZmfLf6IADg9V+3YcuRc5h47zWGkxGRN8tekMS+0hORZUINpaHsfGKkpLhUdXReoykAumUcttXu8Q0m8/qLsf2aOyxJ/9PfRxAz5g+DiYjIm2UvSHa80YsFiYdhUUIe7ZHOdfDVwLZW/9T5lBw/WIiInElNS891c72w4MA8XkGmsCghj9elYRSWPN/V4TEWJkRUEMfOXUL9l+c5PBY3oQ831/NQLErIK0RHlsb213s5PjZ8DtLS1VAiIvJ087ccQ4fxixwei5vQx1AaKggWJflQ1SUZ80iamc7i70qFBGL/+N4Oj9V7aS4Onebqr0Tk6PHv1uFf3653eIwFiedjUUJeRURy/GC57s3F+GbVAUOJiMjTRA+f47Bs/D/a1WJB4iVYlJBXipvQx2Ejv1d/3oK2Y3lnDpG/yz7f7NtH2mFsv+aG0lBhsSghr7X19V547oaGVv9EEu/MIfJXJxKTc3z/bxp9Izo3iDSUiIqCRQl5tad6NMDC57o4PBY9fA7Op6QaSkREJe3/lu9D23ELHR6Lm9AH5cKCDSWiomJRQl6vXlQZ7BvnOAG22ajfMCO2SPsqEpEXiR4+B2PmbHd4jPNHvBeLEvIJAQE5J8AOm7mJl3OIfFR6uub4/n6xV2MWJF6ORQn5lLgJfTCwUx2Hx6KHz0E61zMh8hnLdp1E3ZfmOjy29uWeGNK1nqFE5Cqiyh/W7hATE6OxsbGmY/itI2cvodMEx0WTPn8oBjdcXcVQIiJyhYavzMPl1HSHxzg6Yp6IrFPVmOKehyMl5JNqVCiV4wfVY1/H8nIO0f+3d+fhUlRnHse/LyCLsqgYAWXVgDsSccEFAcG4xsSJmok6iYnRxGgcE40hxijuS2I0UZOJUQf3xJjBDeMCxAVQB5QogkRQUUQBRdlFBN7541TPbfp239vdt6u7uvv3eZ5+iq7l1NuHun3fe+rUOVXq82j+mvSEZGj/bZSQ1BglJVLT5l99FIfttmnrSN/R41m59vMKRSQihbpj6nz6Z8xfM+Enw7jr1P0qFJHERbdvYqLbN8mS7XbO8YN78qvj96xQRCKSj2ytm2odSR7dvhEpQLbbOX996T36jh6PEnOR5HnpnU8aJSTfO6ifEpIap5aSmKilJLnun76A8x94dZN1vzhyF047eIcKRSQi6bK1jrx+6eF0SJtaQpJFLSUiRTph716NZhu+4rHX1WoiUmGz3l/eKCHp3rk9868+SglJnVBLSUzUUlIdxjw8i7FT52+y7j9H9ufHaXPqiEj8srWOTBl9CNtv2aEC0UihStVSoqQkJkpKqsf6DRv5YkbPfoA3Lj+Ctm3UmCgSpwmzF/O9Oxt/V6rvSHVRUpJwSkqqz+WPzubWyW9vsm5gzy48fNZBFYpIpHa5O/1+/lij9VNHH8J2ah2pOkpKEk5JSXXK9UX52NlD2XW7zhWISKT2nHnPy4yf+cEm6wZ068iTPx6W4whJOiUlCaekpLo98sr7/Oi+GY3Wv33VkZhZBSISqX7zlqxk1G+ebbR+zmWH034zdWStZkpKEk5JSW3I1vnugB27cu9pQyoQjUh1ytUCqUfxa4eSkoRTUlI7lq1Zx6BLn2q0/rZv783IXTTBn0hTDvn107z10epG69WRtbYoKUk4JSW157+nvM0lj8xutP6lC0fRtWO7CkQkkly5fl7+edGhbLl52wpEJHFSUpJwSkpqV7ap0wHeuvJIWrVSfxOpb6++t4xjbprSaP2vjhvI8Xv3qkBEUg5KShJOSUlt27jR2eGCxvfIQc3SUp8+XPkZ+1wxodH6PbbvwiM/0mP1tU5JScIpKakPi5avZchVExutb79ZK+ZcdkQFIhIprzXr1rPrRU9k3aYEvX4oKUk4JSX1Zeq8jzjx1hcbre+5VQcm/+yQCkQkEq+1n29g518+nnWbRkOuP0pKEk5JSX2658V3+MW41xqt33qLtrz8y0MrEJFIaTXVMjLjl4ey1RbqxFqPlJQknJKS+vbbCXO5fsIbWbepSVuqUa4+I6CJ80RJSeIpKRGA3zz1Br+bODfrtnlXHEGb1mrilmR7ZcEyvnpz46dpAJ47fwS9tt68zBFJEikpSTglJZJu7JS3GZNlzAbQX5mSTLc8+yZXPjYn67b/vWAk23ZuX+aIJMmUlCSckhLJ5h9zlvCdsdOybrvsa7vzH0P6lDkikU3tMeYJVq5dn3Xb65ceToe2mqNGGlNSknBKSqQpCz5ew9Br/5Fzuyb+k3JasmIt+17Z+NH2FF2P0hwlJQmnpETysWGjs2OOQdgAHjrzQPbstWUZI5J6csG4mdz74rtZtx2xe3f+cPLgMkck1UpJScIpKZFC/eCul3h81qKs21oZvHml/lqVlss14F/Kn08fwpAdupYxIqkFSkoSTkmJFOuNxSv58vXP5tx+8Vd25TsH9itjRFLt3J19r5zIhys/y7nPm1ceSWvN3SRFUlKScEpKpBRG/Ppp3s4y7XvKoz86iN2371LGiKSaXPX31/njM2/l3K5J8qRUlJQknJISKaXmOsaCxoyQ4O4X3uHCBxuPKpxOHVel1JSUJJySEonLHVPnc/HDs5rc54lzDman7p3KFJFU2k2T5vLrJ7OPIJwy+Wcj6LmVklaJh5KShFNSIuXQVOfYlGu+vgff2Kd3mSKScnB3jrlpCjMXLm9yvzu/uy8HD/hCmaKSeqakJOGUlEi5nXTrC0yZt7TJfdq1acXMMYdpBtcqNGfRCg6/4blm9/v9SXtx5B49yhCRSAMlJQmnpEQq6eKHXuOO599pdr/vD9uBnx+xSxkikkKt/XwD+14xgRU5RldN97cz9mdwn63LEJVIdkpKEk5JiSTF828u5Zt/eiGvfc8csSM/PWznmCOSbNasW8/B1z7NR6tyP7abbuaYL9Op/WYxRyWSHyUlCaekRJJo40bnyN89x5xFK/Pav2O7Njx7/gi23qJtzJHVn6nzPuLEW1/Me/9Lv7ob39q/b3wBibSAkpKEU1Ii1WDpqs8YfPmEgo7ZpUdn/vL9IXTWX+l5e/W9ZRxz05SCjjlst27818mD9eiuVAUlJQmnpESq0Yq1n3PQ1ZPy6seQ6ZxR/Tn7kP60quNRQZetWccZd7/M82813eE4m2/v34cxx+ymJESqkpKShFNSIrXi/mkLOP9vr7aojAuP2oVTDuhLm9bV/9TP0lWfcfn41xk3Y2GLypl47jB2/ELHEkUlUllKShJOSYnUskdffZ+z7p1R0jL336ErJ+7Xm0N37Ub7zVqXtOx8LFmxlgdefo97XniXhcs+LVm5ndq1YfzZQ+ndVQOXSe1SUpJwSkqk3rg7t0+Zz2WPzq50KBWxRdvW3PW9/dir91aVDkWk7EqVlLQpRTAiImbGqQf149SDGs9g7O5MmbeUix5+jbc+zD3BYNIdPbAHo4/YWcO1i8RESYmIxM7MOKj/Nkw6d3iz+y5ZsZbHZy3iqdmLeW7uR7HG1aldG0bt2o2j9ujB0AHb0K5N+W8biUgDJSUikijbdm7Pt/bvqzE5ROpQ9XeFFxERkZqgpCQLMzvMzMaZ2Qdmts7MFpvZc2Z2XqVjExERqVVKStKYWSszuwV4HDgKeAP4GzAL2An4QQXDExERqWnqU7KpK4HTgBeBE9z93dQGM2sN7FWpwERERGqdkpKIme0KnAcsBY529026/bv7BmBaJWITERGpB7p90+AMoDVwW2ZCIiIiIvFTS0mDL0fLJ8ysO/BNYADwKTAd+B93X1up4ERERGqdkhLAzNoB/aO3uwAPAp0ydltgZse6+0tlDU5ERKRO6PZNsBWQmi/8BmAOcAAhMdkDeBjoBTxmZttUJEIREZEaVxMtJWZ2LXBMEYeOdPeFbJqcLQcOdffl0fvXzOxYYAYwEDgTuCRHHKcDpwP07t27iHBERETqV00kJcB2hHFECrVZtFyZtu7+tIQEAHffGI1fchMwkhxJibvfAtwCYZbgIuIRERGpWzVx+8bdT3Z3K+I1Pzp+JeFRYIC3c5wmtb57zB9HRESkLtVEUlIiL0fLrjm2p/qSrCpDLCIiInVHSUmD/4mWI3NsHxUtp5chFhERkbqjpKTBHcB7wN5m9vP0DWZ2HHASsAH4fQViExERqXlKSiLu/ilwPOH2zJVmNsfMHjCz6cBfCY8Mn+Pu/6xknCIiIrVKSUkad3+B8Njv7UBHwmPGfYCHgGHuflMFwxMREalp5q4nV+NgZh8C70RvtwE0n048VLfxUd3GR3UbH9VtfJqq2z7u/oWWnkBJSRmY2XR337vScdQi1W18VLfxUd3GR3Ubn3LUrW7fiIiISCIoKREREZFEUFJSHrdUOoAaprqNj+o2Pqrb+Khu4xN73apPiYiIiCSCWkpEREQkEZSUiIiISCIoKSkRM5tvZp7H66ICyhyeR3lD4vxcSRBnPZjZfmY2zsyWmNlaM5trZteaWZdSf44kMrOdzOzHZvZ3M5sX1cFyM3vezM4xs7ZFlFlX162ZnWhmz0X1tsrMppvZmWZW1PermR1uZk+a2cdmtsbMXjOzX5hZu1LHnlRmtpmZjTSz68zsBTP7wMzWmdnCaKTt4UWUObaZa3JODB8lkeKoCzNrFV3306Ofg+XRz8U3CymnTaEnlpweoGEm4UxbA1+J/v2PIspeDDyeY9uHRZRXrUpaD9EPy11Aa2AKsBAYAvwUONbMDnT3JUXGWi0mAtsDawmTTU4DugH7E+riW2Y2yt0/LqLsmr9uzexm4IeE+psIfE6Y1PMmYKSZHe/uGwoo73zgGsI8W08DnwDDgMuBo81spLuvKemHSKZhwFPRvxcBLwGrgV2BrwNfN7PL3D3vP/LSTAHmZVn/QTGBVrmS1IWZtSZMansMsAJ4EmhH+Fm418z2d/ez8yrM3fWK+QWcDzjwrwKPGx4d93SlP0OF66/k9QD0BNYQvvy/mra+DfDn6HzjKv3Zy1C3E4HvAh0z1vcFXovq4Y5K/38l8UX45ejRF3j/tPXdgNnRtv8soLy9gY2EX777pa3vCDwTlXd9pT93mer2EMIfekOzbPsGsD6qjxEFlDk2OuaUSn++Sr9KXRfAuVF5s4Buaev7E5JKT/+ebeql2zfl8d1oeXtFo5B05wAdCL9wH0qtdPf1wOmEbP9rZrZrheIrC3cf6e63u/uqjPXzgR9Eb08o5jZOHUjNJv4zd5+bWunui4EzorejC7iNM5ow8ec17v5iWnmrgO8QEpYfmtmWLY484dx9krsf5+7PZdn2F8IvVYCTyxqYNBK1kpwfvT0juv4BiH4ufha9/UU+5SkpiZmZHQjsRMjs76xwONLga9HynswN7r4CeCRjv3o0I1q2B7pWMpCkMbOewGBgHWEW8U24+zOE24HdCbfBmiuvLXBE9DbbNfkW8DzQFjiy6MBrR+ra7FnRKATCrd5tgffc/dks2/9KuK25j5lt31xh6lMSv1QryWPuXuw9y25mdjHh3v9qYCbwkLsvLUWAVaQk9WBmnYEdo7fTcuw2DTgJ+FKRsdaC/tFyHVBMn5Javm5T18Usd/80xz7TCJ/9S8DUZsrbCdgc+Njd32yivAOj8u4tLNyak7o2i/lOHWFmAwm3xRYDk4Gn3H1jqYKrIqWoi9TPQtbvUndfY2azgEHRa2FThSkpiZGZbQGcEL29rQVF7QyMyVh3o5mNdvcbW1ButSlVPfSNlsuiVpFs3o2W/Qoot9aMjpaPuvtnRRxfy9dt6rp4p4l9CrmGUvu828Q+uiYBM+sOnBK9/VsRRXwry7rZZvbv7j6z6MCqUynqIt+fhUHkce3q9k28TiBkoIuAx4o4fjlwPTCU0AzcCdgLuJXQpP47MzutNKEmWqnroWO0XN3EPqk+Fp0KC7U2mNkphA6Fa4ALCjy8Hq7bUl9DuibzYGZtgLuBLsBEd3+kmUPS/RM4G9iNUN/bAUcDrxCe6pmQz+2FGlHKuijptauWEsDMriU8ylSoke7eVFPUqdHyzqgDZUHcfQYN905TZgCnmdmrwO+Aq83sziL/ko1dKeo2hnqwaFnVcyzEdd2a2Ujgj4T6+b67/6uQwmvhus1Dqa+hmrgmy+C/CI+ZLqDATq7ufkPGqtXAeDN7ivB00xBC5+WzShBnopW4Lkp67SopCbYj3NMt1Ga5NpjZAML9X4jnqZubgYsIY6PsB2TrYJQEJa/bDMXUw8po2bGJfVLbVjaxT6XFcd0eBDxE6FB5trvfXWRsuVTLdducUl9DtXJNxsbMfkv4Q28RIbFeVIpy3X2dmV1FuO7ruhNxkXVR0mtXt28Adz/Z3a2I1/wmik11cJ1c6F+aeca8EUg9hpjYJseY6ja9/GLqIVX2llGn12x6ZeybOKWuWzM7gHCbcQvCY64l7/dRLddtHuZHyz5N7FPINZTap3eJyqspZnYd4XbDh4SEZG4zhxQqNYJpNV+TpVJoXcyPliX5WVBSEoPoue1UB6KWdHBtTuoxzVVN7lX7CqqHqHNr6gmHfXLstm+0zLwNUZOiYd//Trjne6G7Xxvj6Wrhuk1dF7uZWYcc++yTsW9T5gCfAlub2Y459qmrazIluk35E2ApcKi7z47hNLVwTZZKoXXxcrTM+l1qZpsDu0dvm712lZTE40igB6GpqtEYBqVgZnsCAwj38abHcY5q0IJ6SA2YdlKWMjvTMC3AuBYFWAXMbF/gCaAzMMbdr4jxXDVx3br7AsKXcVvg+MztZjaMMIbGIsL4Is2Vt46QFEL2a3IHwngQ64DxRQdeZczsasK0D58QEpJXYjpV6inJXEME1JNC6+J5YAnQ08wOzrL9eMIt42nN9MEM8hn2Va+Ch9wdR/jSvSWPffcl/JU0J8u2s4GuWdbvT2gCd+C+Sn/eMtRnUfUAHBvV7cQs23rRMMz8MWnr2wD3UT/DzA8GlkWf99ICjqv76xY4joZh5r+Ytn5bwnDbjYaZJ3QcnEPo/J5Z3j40DDO/b9r6joR5cOpmmPnoc18WfeZPgMF5HnNVVL9XZawfRHi6pHXG+jaEVpgN0bkOq/TnLkO9FlUXhME/5wBnZSnzPBqGmd82bX3/6Ocj72Hm1dG1xMxsW+Co6G0+t242J3dnxUuB68zsdUJv8zWE/+SBhB7PU4Dvtyjg6lBsPXQh1G37zA3uvsDMTiVMyPegmU0G3if0Ou9DmKSqHur2KUI9LQN6m9nYHPud5+4fpb2v++vW3R8wsz8QhpSfaWYTaJiQrzPwIGFivnTbEOqtUSdNd59mZqMJE/JNNbNJhP+XYYRE50XyHKq72pnZMcCF0dt5wI/MLNuuc9z96rT3PQj12yNjv76EPxY/NrM3gPcItyr3IHQY30joR/VEqT5DgvWluLroTajbbBPPXg8cTGhhnmtmEwmtI6MI3783etp0Hk2qdNZWay/SJibKc//h0f6eZdtPgYcJ/R+WE77wFhOa2k8hI9Ot1Vex9RBtc2B+E2XvR/jl8SHwGeEL8FqgS6U/d5nq1vN89c04Ttdtw+c9kZBorSC0crwEnAm0yrLvGJqZrBA4nJAsfkLoZzKLkIy0q/RnLWOdpn52m3s9nXHc2Gj92Iz1/YAbCCPrLiTM6vwpoeXudvJsiamFV7F1QUNr3Zgc21sRWgJTMzqvIIwQe2Ih8VlUmIiIiEhFqaOriIiIJIKSEhEREUkEJSUiIiKSCEpKREREJBGUlIiIiEgiKCkRERGRRFBSIiIiIomgpEREREQSQUmJiIiIJIKSEhGpembW1sz6VToOEWkZJSUiUtXMrAvwd8Isu3GU39HMppnZznGULyINlJSISNUys06EicIedPf7W1DO62bWMds2d18FnA48amaDiz2HiDRPSYmIVCUzawX8GXjd3W9sQTm9gfVR8pGVu88ARgNPmFnPYs8lIk1TUiIi1eonwGDgjBaWMwx4prmd3P0B4DHgXjNr3cJzikgWSkpEpOpEnVovAS5z9+UtLO5g8khKIqOBvYBzW3hOEcnC3L3SMYiIFMTMxgJHAL3cfV0Ly/oXMNTdl+S5/2+AbwP93H1FS84tIptSS4mIVBUz6wOcBIwtQULSndCfJK+EJPJ7YCvgxy05t4g0pqRERKrNKUAb4C8lKCuv/iTp3H0eMB34QdTZVkRKRD9QIlJtTgQ+BGaUoKxC+pOkewToDhxYghhEJKKkRETKwsy+aGZ3mdkMM/uXmT1qZrtH27aJtr1qZi9F/+6WpYw+wABgqhfRIS6KYYyZPW9m0wjjj4wxs9lm9lsz2yHPoiZHy68XGoOI5KakRERiZ2bDgXHAf7v7l4DdgI3AVDMbBDwO/NndBwI3AMcCt2Upami0fLXA829rZncTbrt0BE4AjgbedPddgFFAJ2C2mZ2WR5EvAw4cUEgcItK0NpUOQERqm5ltD/wB+LK7LwBw9/VmdhvwFUKrw6XuPt7MdgTGEv5g2i5LcYOi5dwCzn8IcB/wPrB31CcEMzsRmBTF8z7wXTNrC9xiZsubGiHW3Zeb2WJgdzNr7e4b8o1HRHJTS4mIxO1y4PJUQpKmU9q//xAt1wLLgTXApVnK2jFavp/Pic3scGA88B7hsd95aZtHAhMzDrk4Wv7WzJr7o20B0AHon08sItI8JSUiEhsz2xwYQRgOPtNe0fIZd18J4O4LgR7Alu7+YJZjUq0nH+Vx7n40PKHz71mGkR9BmDfn/7n7m8BSQifW5m7NLI6W2zYXi4jkR0mJiMTG3dcAe+S4vTE8Wk7KOOYzd/88R5Gp1pW1eZz+ZqAzcIu7b3K7J7pN9Im7L81yXKqFJNvto3SfRsst84hFRPKgpEREYpVqBUlnZl2BPaO3kzK3N6FdtMyVtKTK7w8cHr29Ncsuo2h86wYz6wx0id5+0EwsqcSoU5N7iUjelJSISCWMIHz/fAy8UsBxa6Jluyb3CuOHGPCRu8/Msn0k2ZOhVCKzmubHQWkbLVs0qqyINFBSIiKVMDJaPu3uGzM3mtlmOY5LzTXToZnyU2OcvJWlbCMkLc9lOe6UaHlnHvPapGLI7KsiIkVSUiIisTKztlmGYx8VLZ/Osn8Pco9D8k60bDSwWo79WmfZtidhfJLVGecdQZjkbwFwQTPlA3SNlgvz2FdE8qCkRERiY2bfJrRuzDGzdtG6AcAXo12mZTnsDOCeHEWmHunt1cypHwdWAjubWceMbY0eBTazXtE5Pwa+5u7Lmik/PYZGrTEiUhwlJSISp+sI/T960vBUy7k09A1Zk76zme0D/BtwY47ypkfLnZo6aZRUXARsAZyfsXmTTq5mNpAw/81qYIS7v9xU2dEx7QlP57yZ5VFjESmSkhIRidNiwrgfJ7j7ajP7HmHumq8A64H/H9LdzI4C7gC+4e7Lc5Q3mTBmJITtAAAB4klEQVS8++DmTuzuNwBjgNFmdrmZbRn1VfkS8KKZ7WxmNxL6ljwIDHL3fIevH0RIsiY3t6OI5E9JiYjE6d8Iv/QvMbOZwEDgKHefBBwGDDKzOdG244DD3H1WrsLc/WPgeWCfqLWiSe5+CTCEMBjaC8BswiO8/wv8HpgPDHD3n2T2MWnGkGj5ZAHHiEgzrIiJNkVEKsbMziLc3jna3ccXeOwvgTXufl0LY3iKMOJrN92+ESkdtZSISLW5m9B59ptFHDsc+EdLTm5m2wDDgPuVkIiUlpISEakqUSfWm4HjzSzveWeiGYB3Bv7ZwhBOJ/QnuaaF5YhIBiUlIlKNriFMypfPeCIp+wHTsw3Wli8z2wI4izC42pxiyxGR7JSUiEjViZ7OORP4oZntkudhw2nhrRvC48VtCI81i0iJKSkRkark7g8SxkG5LzUwWzP2JsskfPkys72A84CTcswuLCItpKdvRKRqRfPYjCW0XpzsMX2hmdkXgCnAr9z9T3GcQ0SUlIhIlYsSk6uBT9z96hjK70CYUfhmd7+71OWLSAMlJSIiTYhGge3r7nMrHYtIrVNSIiIiIomgjq4iIiKSCEpKREREJBGUlIiIiEgiKCkRERGRRFBSIiIiIomgpEREREQS4f8ALi5sgY7+EsIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(12, 8))\n",
    "plt.plot(xpos, ypos, ls = '-', lw = '2', color='C0', label = 'geodesic')\n",
    "plt.plot([0.0], [0.0], '*', ms = 40, color='C1')#, label = 'orbit')\n",
    "plt.legend(loc = 'upper left', fontsize = 28, frameon=True, ncol = 1)\n",
    "#plt.axis([0,2,0,9])\n",
    "#plt.grid(True)\n",
    "#plt.title(\"Solution of $y'=y , y(0)=1$\", fontsize = 28)\n",
    "plt.tick_params(labelsize=22)\n",
    "plt.xlabel(r'$x(\\phi)$', fontsize = 32)\n",
    "plt.ylabel(r'$y(\\phi)$', fontsize = 32)\n",
    "plt.gca().set_aspect('equal', adjustable='box')\n",
    "plt.show()\n",
    "plt.close(fig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.6\n",
      "1.2566370770628055e-06 4.600000315654222\n",
      "2.1205750684316627e-06 4.600000162337097\n",
      "2.591813970553858e-06 4.600000057831452\n",
      "3.455751961922715e-06 4.600000073906107\n",
      "3.15000003917177e-05\n",
      "47.03522785763301\n"
     ]
    }
   ],
   "source": [
    "list_perih = list()\n",
    "phi_perih = list()\n",
    "r_old = r[0]\n",
    "print(r[0])\n",
    "orbit = 0\n",
    "for i in range(1, len(r)-1):\n",
    " if (r[i] < r_old) and (r[i] < r[i+1]):\n",
    "  list_perih.append(r[i])\n",
    "  phi_perih.append(phi[i] - orbit*2.0*np.pi)\n",
    "  orbit += 1\n",
    "  print(phi[i]- orbit*2.0*np.pi, r[i])\n",
    " r_old = r[i]\n",
    "    \n",
    "sum_angle = 0.0\n",
    "for n in range(1, int(orbit)):\n",
    " sum_angle = sum_angle + (phi_perih[n] - phi_perih[n-1])*180.0/np.pi\n",
    "\n",
    "print(sum_angle/(Norbits-1)) # degrees per orbit\n",
    "\n",
    "#In degress per century on Earth\n",
    "angle = sum_angle/(orbit) * 3600.0 * 365.0/88.0 * 100.0\n",
    "\n",
    "print(angle)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_per = min(r)\n",
    "r_af = max(r)\n",
    "\n",
    "e = (r_af - r_per)/(r_af + r_per)\n",
    "print(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}