{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ed177532-5b0b-4fe5-a697-b1c3b8058501",
   "metadata": {},
   "source": [
    "# Calculating transit timing variations\n",
    "\n",
    "In this example we will be using heyoka.py's [event detection](<./Event detection.ipynb>) system to find the transit times in a two-planet system. This example is directly inspired by a [notebook](https://rebound.readthedocs.io/en/latest/ipython_examples/TransitTimingVariations/) in the documentation of the [REBOUND](https://rebound.readthedocs.io) software package.\n",
    "\n",
    "We begin as usual with importing heyoka.py and a couple of extra useful modules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "fcd9fcfa-4a18-4083-815b-87e188c900b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pykep as pk\n",
    "import heyoka as hy\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d739c1a-2434-4038-93ff-d6bcaeed69f6",
   "metadata": {},
   "source": [
    "Next, we define the initial conditions for the two planets in the system. The inner planet is on a low-eccentricity orbit, while the outer planet is on a circular orbit. Both orbits have zero inclination, so that the motion is constrained to the $xy$ plane. We use adimensional units ($G=1$), the mass of the central star is $1$ and both planets have a mass of $10^{-5}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d3cb090b-76b3-46fa-917c-8b6b971b698a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initial cartesian conditions for the two planets.\n",
    "# We use pykep to convert from Keplerian elements.\n",
    "r1, v1 = pk.par2ic([1.0, 0.1, 0.0, 0.0, 0.25, 0.0])\n",
    "r2, v2 = pk.par2ic([1.757, 0.0, 0.0, 0.0, 0.0, 0.0])\n",
    "\n",
    "# Masses of the 3 bodies in the system.\n",
    "masses = np.array([1, 1e-5, 1e-5])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b3f0b64-2a0e-4eac-abfb-ec977fcc1ff8",
   "metadata": {},
   "source": [
    "We can now proceed to the definition of the system of differential equations via the ``model.nbody()`` helper:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0d29f42b-f6c8-424f-89a9-31bc703e9c75",
   "metadata": {},
   "outputs": [],
   "source": [
    "sys = hy.model.nbody(3, masses=masses)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "638e6bbd-93bf-499d-8d9d-378356ab9fa6",
   "metadata": {},
   "source": [
    "The next step is the definition of the event for the detection of the planetary transit. In this setup, we assume that the observer of the system is in the direction of the positive $x$ axis, so that a planetary transit happens when the $y$ coordinate of the planet is zero. Thus, in order to detect the transit of the inner planet, we will use the trivial event equation\n",
    "\n",
    "$$\n",
    "y_1 = 0.\n",
    "$$\n",
    "\n",
    "Because $y_1$ is zero also when the planet transits *behind* the star, we will further constrain the event by imposing a *positive* direction for the event. That is, the event will trigger only if $y_1 = 0$ and the value of $y_1$ is switching from negative to positive. Alternatively (and more robustly), we could compute the value of $x_1$ when the event triggers via dense output and ignore the event if $x_1<0$, however for this simple example specifying the event direction will suffice.\n",
    "\n",
    "Finally, we will log the trigger time of each event in a global list for further analysis.\n",
    "\n",
    "Let us see the definition of the event:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "62eb3b81-30d6-4f56-b4c0-8cf141ddc749",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Global list of transit times.\n",
    "tt_list = []\n",
    "\n",
    "\n",
    "# The callback that will be invoked\n",
    "# when a transit is detected.\n",
    "def cb(ta, time, d_sgn):\n",
    "    # Append the transit time to the\n",
    "    # global list.\n",
    "    tt_list.append(time)\n",
    "\n",
    "\n",
    "# Define the event object.\n",
    "tt_event = hy.nt_event(\n",
    "    hy.expression(\"y_1\"),  # The event equation.\n",
    "    cb,  # The callback.\n",
    "    direction=hy.event_direction.positive,  # The event direction.\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e027a4f-4313-4e06-abb5-c411f8bdcf75",
   "metadata": {},
   "source": [
    "We can now proceed to the creation of the integrator object. After creation, we will adjust the system state so that the centre of mass of the system is in the origin with zero velocity:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d779db8a-c2e0-4156-ac83-410caa6b9c56",
   "metadata": {},
   "outputs": [],
   "source": [
    "ta = hy.taylor_adaptive(\n",
    "    sys, [0.0] * 6 + list(r1) + list(v1) + list(r2) + list(v2), nt_events=[tt_event]\n",
    ")\n",
    "\n",
    "# Adjust the COM.\n",
    "com_x = np.sum(ta.state[0::6] * masses) / np.sum(masses)\n",
    "com_y = np.sum(ta.state[1::6] * masses) / np.sum(masses)\n",
    "com_z = np.sum(ta.state[2::6] * masses) / np.sum(masses)\n",
    "\n",
    "com_vx = np.sum(ta.state[3::6] * masses) / np.sum(masses)\n",
    "com_vy = np.sum(ta.state[4::6] * masses) / np.sum(masses)\n",
    "com_vz = np.sum(ta.state[5::6] * masses) / np.sum(masses)\n",
    "\n",
    "ta.state[0::6] -= com_x\n",
    "ta.state[1::6] -= com_y\n",
    "ta.state[2::6] -= com_z\n",
    "ta.state[3::6] -= com_vx\n",
    "ta.state[4::6] -= com_vy\n",
    "ta.state[5::6] -= com_vz"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc5a5c80-5ba5-4118-9c3c-02e2c24602e9",
   "metadata": {},
   "source": [
    "Let us now propagate the system state for a few time units. Whenever a transit is detected, the callback will be invoked and the global list of transit times will be updated:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b155c48e-3b29-4c44-bab3-461e9984fa51",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<taylor_outcome.time_limit: -4294967299>,\n",
       " 0.22961192655196333,\n",
       " 0.4928525973770144,\n",
       " 3344,\n",
       " None)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ta.propagate_until(1100.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b907afcc-a795-4424-b7f3-36bd10b66254",
   "metadata": {},
   "source": [
    "We can now proceed to the analysis of the transit times. First we will run a linear fit to remove the linear trend from the transit times:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9e7e4a20-6f6b-42f8-8c83-8ff78c1d314c",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = len(tt_list)\n",
    "\n",
    "A = np.vstack([np.ones(N), range(N)]).T\n",
    "c, m = np.linalg.lstsq(A, tt_list, rcond=-1)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7d1620e-2750-47ad-b023-ac29817af170",
   "metadata": {},
   "source": [
    "We are now ready to visualize the transit time variations (TTVs):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c8eb4e22-221c-4d8b-ae65-76c9613e35ca",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1oAAAHFCAYAAADrDYghAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABjpUlEQVR4nO3deXhUVZ7/8U+xJYimEJCEPYALm4KCSAAHm5HFDcSZUbA74KggrTQC6iDdQwt224j+VNxYxK0dbKAXUWxtJDaIIAkokFYkooMoSycgCAmiQEju7w+mIpVUpba7Vr1fz8PzkMqtqpNTy73fc873e3yGYRgCAAAAAJimjtMNAAAAAIBkQ6AFAAAAACYj0AIAAAAAkxFoAQAAAIDJCLQAAAAAwGQEWgAAAABgMgItAAAAADAZgRYAAAAAmIxACwAAAABMRqAFAAAAACbzVKD1/vvv67rrrlPLli3l8/n0+uuvR7zPmjVr1LNnT6Wnp6tDhw6aP3++9Q0FAAAAkNI8FWgdPXpU3bt31zPPPBPV8Tt37tTVV1+tyy+/XFu2bNEvf/lLTZw4UX/5y18sbikAAACAVOYzDMNwuhHx8Pl8WrZsma6//vqwx0ydOlXLly9XUVFR1W3jx4/XP/7xD+Xn59vQSgAAAACpqJ7TDbBSfn6+Bg8eHHTbkCFD9MILL6i8vFz169evcZ/jx4/r+PHjVT9XVlbq22+/VdOmTeXz+SxvMwAAAAB3MgxDR44cUcuWLVWnTu2LA5M60CopKVFmZmbQbZmZmTp58qQOHDigFi1a1LjPrFmzNHPmTLuaCAAAAMBjdu/erdatW9d6TFIHWpJqzEIFVkqGm52aNm2apkyZUvVzaWmp2rZtq927dysjI8O6hgIAAABwtbKyMrVp00ZnnXVWxGOTOtDKyspSSUlJ0G379+9XvXr11LRp05D3SUtLU1paWo3bMzIyCLQAAAAARJVS5Kmqg7HKyclRXl5e0G0rV65Ur169QuZnAQAAAIAZPBVofffddyosLFRhYaGkU+XbCwsLtWvXLkmnlv2NHj266vjx48fr66+/1pQpU1RUVKQXX3xRL7zwgu69914nmg8AAAAgRXhq6eBHH32kn/zkJ1U/B3KpxowZo5dfflnFxcVVQZcktW/fXm+//bYmT56sZ599Vi1bttRTTz2lf/u3f7O97QAAAABSh2f30bJLWVmZ/H6/SktLydECAAAAUlgssYGnlg4CAAAAgBcQaAEAAACAyQi0AAAAAMBkBFoAAAAAYDICLQAAAAAwGYEWAAAAAJiMQAsAAAAATEagBQAAAAAmI9ACAAAAAJMRaAEAAACAyQi0AAAAAMBkBFoAAAAAYDICLQAAAAAwGYEWAAAAAJiMQAsAAAAATEagBQAAAAAmI9ACAAAAAJMRaAEAAACAyQi0AAAAAMBkBFoAAAAAYDICLQAAAAAwGYEWAAAAAJiMQAsAAAAATEagBQAAAAAmI9ACAAAAAJMRaAEAAACAyQi0AAAAAMBkBFoAAAAAYDICLQAAAAAwGYEWAAAAAJiMQAsAAAAATEagBQAAAAAm81ygNXfuXLVv317p6enq2bOn1q5dW+vxr776qrp3764zzjhDLVq00H/+53/q4MGDNrUWAAAAQCryVKC1dOlSTZo0Sb/61a+0ZcsWXX755brqqqu0a9eukMevW7dOo0eP1m233aZPP/1Uf/rTn/Thhx/q9ttvt7nlAAAAAFKJpwKtxx9/XLfddptuv/12de7cWXPmzFGbNm00b968kMcXFBQoOztbEydOVPv27dW/f3/dcccd+uijj2xuOQAAAIBU4plA68SJE9q0aZMGDx4cdPvgwYO1fv36kPfp27ev9uzZo7fffluGYWjfvn3685//rGuuuSbs8xw/flxlZWVB/wAAAAAgFp4JtA4cOKCKigplZmYG3Z6ZmamSkpKQ9+nbt69effVV3XTTTWrQoIGysrLUuHFjPf3002GfZ9asWfL7/VX/2rRpY+rfAQAAACD5eSbQCvD5fEE/G4ZR47aAbdu2aeLEifr1r3+tTZs2acWKFdq5c6fGjx8f9vGnTZum0tLSqn+7d+82tf0AAAAAkl89pxsQrWbNmqlu3bo1Zq/2799fY5YrYNasWerXr5/uu+8+SdJFF12kRo0a6fLLL9dvf/tbtWjRosZ90tLSlJaWZv4fAAAAACBleGZGq0GDBurZs6fy8vKCbs/Ly1Pfvn1D3uf7779XnTrBf2LdunUlnZoJAwAAAAAreCbQkqQpU6bo+eef14svvqiioiJNnjxZu3btqloKOG3aNI0ePbrq+Ouuu06vvfaa5s2bpy+//FIffPCBJk6cqN69e6tly5ZO/RkAAAAAkpxnlg5K0k033aSDBw/qwQcfVHFxsbp166a3335b7dq1kyQVFxcH7al1yy236MiRI3rmmWd0zz33qHHjxho4cKBmz57t1J8AAAAAIAX4DNbQ1aqsrEx+v1+lpaXKyMhwujkAAAAAHBJLbOCppYMAAAAA4AUEWgAAAABgMgItAAAAADAZgRYAAAAAmIxACwAAAABMRqAFAAAAACYj0AIAAAAAkxFoAQAAAIDJCLQAAAAAwGQEWgAAAABgMgItAAAAADAZgRYAAAAAmIxACwAAAABMRqAFAAAAACYj0AIAAAAAkxFoAQAAAIDJCLQAAAAAwGQEWgAAAABgsnpONwAAAAAA3Kai0tDGnd9q/5Fjan5Wunq3bxLT/Qm0AAAAAOA0K7YWa+ab21Rceqzqthb+dN37kzZRPwaBFgAAAAD8nxVbi/XzRZtlVLu9pPSYpiz9R9SPQ44WAAAAAOjUcsGZb26rEWRJCnlbbZjRAgAAAOBKofKk6tbxWfZ8G3d+G7RcsLpYgi0CLQAAAACuEy5P6oHrumhotxaWPOf+I+GDrFixdBAAAACAqwTypKrPLpWUHtPPF23Wiq3Fljxv87PSTXssAi0AAAAArhFNntTMN7epojLWrKnIerdvohb+dIVbnBjLokUCLQAAAACuEU2eVHHpMW3c+a3pz123jk8PXNdFUs2gKtbMMAItG1RUGsrfcVBvFO5V/o6DlkTfAAAAQDKINk/KzHyq06/X/Q0b6NmbL1GWP3gZYZY/XY/f1D3qx6QYhsWcSOIDAAAAvCraPCmz8qnCXa9Pv6azzm6UFlTx8Oh3R6J+XGa0LORUEt/pmE0DAACAl0STJ9XCfyrwSVRt1+t3/WGLSn84oeE9WimnY9OYy8ozo2WRSEl8Pp1K4hvUJcuyvQCYTQMAAIBXnL5n1shL22rOu5/Lp+C9qwJXzQ9c1yXha2irr9cJtCwSSxJfTsempj9/IDqv/sYJzKbN+9klBFsAAABwhVATBI3PqC9JOvx9edVtWSZOGlh9ve65pYNz585V+/btlZ6erp49e2rt2rW1Hn/8+HH96le/Urt27ZSWlqaOHTvqxRdftLydTiTxBThZEhMAAADJxepUlHDL90q/L9fh78s1+crz9OTIHlo8to/WTR1o2mSB1dfrnprRWrp0qSZNmqS5c+eqX79+WrBgga666ipt27ZNbdu2DXmfG2+8Ufv27dMLL7ygc889V/v379fJkyctb6vdSXync3o2DQAAAMnB6lSUaJbvLflwt9ZNHWh6uo3V1+uemtF6/PHHddttt+n2229X586dNWfOHLVp00bz5s0LefyKFSu0Zs0avf3227ryyiuVnZ2t3r17q2/fvmGf4/jx4yorKwv6Fw87k/iqc3I2DQAAAMnBjsJuTu6ZZfX1umcCrRMnTmjTpk0aPHhw0O2DBw/W+vXrQ95n+fLl6tWrlx555BG1atVK559/vu6991798MMPYZ9n1qxZ8vv9Vf/atGkTV3uj2ezMjCS+UJycTQMAAID32ZWK4uQEgdXX654JtA4cOKCKigplZmYG3Z6ZmamSkpKQ9/nyyy+1bt06bd26VcuWLdOcOXP05z//WXfddVfY55k2bZpKS0ur/u3evTvuNg/t1kLzfhZ6szMri1E4OZsGAAAA77NrpsnpCQIrr9c9laMlST5fcPhgGEaN2wIqKyvl8/n06quvyu/3Szq1/PDf//3f9eyzz6phw4Y17pOWlqa0tLS423d6WcrmZ6VrUJcsDeqSFXRb7/ZNLCvpLv0Ynf980WZLS2ICAAAgOdk10xSYICgpPRZy9synU0GPlRMEQ7u1sOR63TOBVrNmzVS3bt0as1f79++vMcsV0KJFC7Vq1aoqyJKkzp07yzAM7dmzR+edd56pbXTTvlWB6Lx6e8wsiRlJ9aDT6gATAAAA5rBrpsktEwR16/hMLxLnmUCrQYMG6tmzp/Ly8jRixIiq2/Py8jR8+PCQ9+nXr5/+9Kc/6bvvvtOZZ54pSfr8889Vp04dtW7d2tT2uXHfKqui82i4KegEAABAdAID5SWlP6hJowY6dPSE5TNNbpggsILPMAzPbKa0dOlS5ebmav78+crJydFzzz2nhQsX6tNPP1W7du00bdo07d27V6+88ook6bvvvlPnzp3Vp08fzZw5UwcOHNDtt9+uAQMGaOHChVE9Z1lZmfx+v0pLS5WRkRHymIpKQ/1nrwq7jjXwRrSiLGX1drhhBilc0BloCZslAwAAuE+ogfJQrLqmc8u1bG2iiQ0CPDOjJUk33XSTDh48qAcffFDFxcXq1q2b3n77bbVr106SVFxcrF27dlUdf+aZZyovL0+/+MUv1KtXLzVt2lQ33nijfvvb35raLjfsW+WWGaRo9kKY+eY2DeqS5boPDgAAQKoKN1AeilUzTVYs33OSp2a0nBBN1PpG4V7dvaQw4mM9ObKHhvdoZXIL3TWDlL/joEYtLIh43OKxfZLqgwQAAOBVkVZnSVKTRvU1/dquyspw50yTXWKZ0fJMeXc3c7IspV17HESLzZIBAAC8JdLqLEn69mi5sjLSldOxacoGWbEi0DKBk/tWObmbdihO74UAAACA2DBQbg0CLRNYvat0bdz2wWCzZAAAAG9hoNwaBFomsXJX6dq47YPhZNAJAACA2DFQbg1PVR10Oyf2rXLDbtrVJeteCAAAAE6xsvS5WzYNTjZUHYwglsoiTglUHZRCfzCc2rfKC3shAAAAuJ1d2/i4ZbsgN4slNiDQisALgZbk/AeDoAoAAMB8dm/jwzVd7Qi0TOSVQEty7oPhdJAHAACQjCLtbxVIEVk3dSDBkE1iiQ3I0UoiTuymHW6UpaT0mH6+aLNjyxYBAAC8LpZtfOy+BnSKl2bcCLQQt0ibJft0arPkQV2ybP8AeOlDCAAAEIrbtvFxmtdWURFoIW5uHWXx2ocQAADgdIEB4y/2HYnq+FTY38qLq6gItBA3N46yePFDCAAAEBBqwDgcJ7bxcYKbV1HVhg2LETe3bZYc6UMonfoQVlRS/wUAALhPYMA42iBLSo39rWJZReUmBFoOqKg0lL/joN4o3Kv8HQc9e+Hvtl3EvfohBAAAqG3AOJQsf3rKrNRx4yqqaLB00GZO5w+ZWSTCbbuIe/VDCAAAEGnAOGDCT85Vv3ObpVShL7etoooWgZaNnM4fsiLIG9qtheb97JIaj5vlQPEJr34IAQAAoh0IPi/zzJQp5R4QWEVVUnos5IyfW3PVCLRs4nQSn5VB3tBuLTSoS5bj5dS9+iEEAABgwDg8t62iihY5WjZxMn/IjiIRgc2Sh/dopZyOTR15owc+hJJq5I25+UMIAADgttx3twmsosryBweabs5VY0bLJk7mD7l1vysruGkpIwAAQLS8OmtjJ7esoooWgZZNnJwOTrUiEV77EAIAAEjeGTA2s7harAKrqLyAQMsmTuYPpeKaXy99CAEAAALcPmDsdAVtLyFHyyZO5g+5Zc2vVfuHJcu+ZAAAAJI7ct9DCbehcqC42oqtxQ61zJ2Y0bKRU9PBbljza9XoB6MqAAAA1nO6grYX+QzDYPi/FmVlZfL7/SotLVVGRoYpj+nUulangpJwpeUDf3G8lWKselwAAAAEy99xUKMWFkQ8bvHYPkmdvhFLbMCMlgOcyh9yYs2vVaMfjKoAAADYJ9WKq5mBQCvF2B3kWVVaPpVK1gMAAHdwstqe01KxuFqiCLRgKatGPxhVAQAAdkr1vHAnK2h7FVUHYSmrRj/cPKpCFUQAAJIL1facraDtVcxowVJWjX64dVQl1Ue7AABINuSF/8grGyq7BYEWLGVVaXk3lKyvLlwVxMBoF1UQAQDwjkA+1gf/+w154adx+4bKbkKgBctZNfrhplEVRrsAAEgeoVaoRJJKeeFOVdD2Gs8FWnPnztWjjz6q4uJide3aVXPmzNHll18e8X4ffPCBBgwYoG7duqmwsND6hiKIVaMfbhlVoQoiAADJIdwKlUiotofqPBVoLV26VJMmTdLcuXPVr18/LViwQFdddZW2bdumtm3bhr1faWmpRo8erX/913/Vvn37bGwxTmfV6IcbRlWogggAgPfVtkIlHKrtIRxPVR18/PHHddttt+n2229X586dNWfOHLVp00bz5s2r9X533HGHbr75ZuXk5NjUUqQaN1dBBAAA0Ym0QqU6qu2hNp4JtE6cOKFNmzZp8ODBQbcPHjxY69evD3u/l156STt27NADDzwQ1fMcP35cZWVlQf/gfVaXXA9UQQz3FevTqeqDjHYBAOBesa48yfKnU+wKYXlm6eCBAwdUUVGhzMzMoNszMzNVUlIS8j5ffPGF7r//fq1du1b16kX3p86aNUszZ85MuL1wDztKrruxCiIAAIhNtCtPJvzkXPU7t5mj1fYCVRGp/BcfO/rPM4FWgM8X3AGGYdS4TZIqKip08803a+bMmTr//POjfvxp06ZpypQpVT+XlZWpTZs28TcYjrKz5LqbqiACAIDYRbtP5+RB5zsa1CTLvp1OBYt29Z/PMAxz11BZ5MSJEzrjjDP0pz/9SSNGjKi6/e6771ZhYaHWrFkTdPzhw4d19tlnq27dulW3VVZWyjAM1a1bVytXrtTAgQMjPm9ZWZn8fr9KS0uVkZFh3h8Ey1VUGuo/e1XYtdaBL8t1Uwea+qFmhAkAAO8KDNJKoVeoOL1UMNwgslvaFy2ngsVE+y+W2MAzOVoNGjRQz549lZeXF3R7Xl6e+vbtW+P4jIwMffLJJyosLKz6N378eF1wwQUqLCzUZZddZlfT4ZBYSq6bKVAFcXiPVsrp2JQgCwAADwmsUMnyBy8jdEM+VqR9O6VT+3aanYtutkCwU/06LbDiaMXWYkue1+7+89TSwSlTpig3N1e9evVSTk6OnnvuOe3atUvjx4+XdGrZ3969e/XKK6+oTp066tatW9D9mzdvrvT09Bq3IzlRch0AAMTDLft0VpcM+3ZGCnZ8OhXsDOqSZXp/291/ngq0brrpJh08eFAPPvigiouL1a1bN7399ttq166dJKm4uFi7du1yuJVwC0quAwCAeLlhn87qkmEQ2exgJ5aUDbv7z1OBliTdeeeduvPOO0P+7uWXX671vjNmzNCMGTPMbxRcKdqE1kRKrpOPBQAA7JIMg8hmBjux5nnZ3X+eC7SAaFldcj1ZKv4AAABvsGMQ2WpmBTvxVJa2u/88UwwDiIdVCa1OJXECAIDUFRhEln4cNA6IdxC5otJQ/o6DeqNwr/J3HLS8kEYg2AnXQp9ODVzXFuzEW9TCiv6rjWfKuzslUglHlo55g5mvk1Nl4wEAACTzVtU4XWJdiq+Efv6Ogxq1sCDi8ywe2ydknlcif3cs5d1ZOpgAM9+cBGzWMjOhNRkq/gAAAO8yoypiPEvvzBJYcVT9OjoryuvoRPO87KoqSaAVJzPfnE7n+hDkxSYZKv4AAABvS2QQ2ckS6wGJBDtm5HnZUVWSQCsOZr45nRxNCDw/BR1ikwwVfwAAQOpyy+qceIMdrxQFoRhGHGJ5c9Ym0d2pE01epKBDfMxI4gQAAPGxu3hDMvL66hy7i1rEixmtOJj15kxkNCHRmSg3TBl7ldVl4wEAQGisxDFHMqzOSTTPyw4EWnEw680Zb8BmxnJDt0wZe1UsH247c+DItwMAJCun0y2SidlL75y6/rCrqEW8CLTiYNabM56AzayZKK9PGbtBNB9uO0feGOUDACSbwAV8SekP+s1bRazEMYmZq3Ocvv6wo6hFvMjRioNZ60LjyfUxKz/MLVPGXl9nHfhwD+/RSjkdm9YIsuzKgSPfDgCQbFZsLVb/2as0amGBJv/xH/r26Imwx0Z7/YMfBVbnZPmDr/Wy/OlRzw5y/VE7ZrTiZMa60HhGE8yaiXJDtRanR0CsZGcOHPl2AIBkE26ZYCSsxIlNIkvvuP6IjEArAWasC401YDNrJsrpgg7Jvs7azhw48u0AAMmktgv4SNxcvMGt4l16x/VHZARaCTJjXWgsAZuZM1FOVWtJhREQO3PgyLcDACSTSBfwoTi1b5JZRSC8WMyK64/ICLRcItqAzeyZKCeqtaTCCIidOXBuybcDAMAMsV6YO7W1ilkpEF5NpeD6IzKKYXiQGcmLp6utoIMVUmEExM5NjdlAGQCQTGK9MI/3+icRZhWB8HIxCa4/ImNGy6Oc3jcgkSnuVBgBsTMHzul8OwAAzBQpTUKSmjSqr+nXdlVWhv3L7MxKgfB6KgXXH5Exo+Vhds9EBZxebvXuJYUatbBA/WevinrUJVVGQMyeeXTLcwEAYKVI2+j4JP1uxIUacbG91z8BZm21Y9bjOInrj9oxo4WYmFEtMJVGQOyceXR6lhMAALM4VbArGmalQCRLKoUbrz/cUlwkqkBrypQpMT/wf//3f6tJE2/PSCCYmVPcbv4CNZudO5a7eXd0AABi4cYLeMm8FIhkSqVw0/WHm4qL+AzDiLhNQZ06dZSTk6MGDRpE9aDr1q3T9u3b1aFDh4Qb6LSysjL5/X6VlpYqIyPD6eY4Kn/HQY1aWBDxuMVj+0T9YXPLiAMAAEA0KioN9Z+9KuJWO+umDoyYo2XG4+BHkTa6vq1ftq7skpXQ9WYssUHUSweXLVum5s2bR3XsWWedFe3DwkOsmOJ20wgIAABAJGalQKRSKoUdotno+oUPvtILH3xl2wxXVMUwXnrpJfn9/qgfdMGCBcrMzIy7UXCnZJrilk59IPN3HNQbhXuVv+OgKirj2YMeAACkGrOKQFBMwjyxbHRtV/n8qJYOpjKWDv4omaa43bR+N14suwQAwFlmnYs5pyfujcK9untJYdTHx3vdasnSwYDdu3fL5/OpdevWkqSNGzfqD3/4g7p06aJx48bF+nDwkGSZ4jajcqLTkiFQBADA68xKgTDjcdwerFndvlhXVJ1ePt+qNJaYA62bb75Z48aNU25urkpKSjRo0CB17dpVixYtUklJiX79619b0c6U49YPi9erBXp9c0ApOQJFAABgHjcOwJ5+LfvVge+1eOMulZRZ175oNroOxcry+TEvHTz77LNVUFCgCy64QE899ZSWLl2qDz74QCtXrtT48eP15ZdfWtVWRzixdNCNH5bq3BoIRmJF5UQ7BZZvhluD7KXlmwAAIHHhBmADVwFODMCGupatzor2BfpCUtTBVqzXfLHEBlEVwzhdeXm50tLSJEnvvvuuhg0bJknq1KmTioutTShz0sYvv7WlWELgDVL9jWlX0l60AlPcw3s4syt7vLy+OWAy7CIPAADMEWmljnRqpY6dBb/CXctWZ0X7whUXCcWnUxMZvdtbt+9vzIFW165dNX/+fK1du1Z5eXkaOnSoJOmf//ynmjZ13wyAWW79/YfqP3uVpYGOWz4syVyNz+uVE2MNFJP5tQQAINW5bQA2mhLrp7OifUO7tdC6qQO1eGwf3dovW9KPs2cBdtUWiDlHa/bs2RoxYoQeffRRjRkzRt27d5ckLV++XL179za9gW5idQ5MLB8Wq5a1eWHZYiIird8NLL2zcnQjEbEEisn+WgIAYCYvpkW4baVOLCXWT2d2+wIrr3I6NlXv9k0cqy0QU6BlGIbat2+vr7/+WhUVFTr77LOrfjdu3DidccYZpjfQTQIX5vf/5ROdlV5ffTqYu2TO6Q9LKhRZ8HrlxGgDxUNHT+iuPyT3awkAgFm8OjjptpU68V6jWtm+od1aaFCXLEeC6JiWDhqGofPOO0/79u0LCrIkKTs7W82bNze1cW51+Idy/fT5DaYvJXTyw+KWZYt28PLmgIFAUQo/DT79ms76zVup8VoCAJAor+THhxIYgA0XMtiRh3S6WK9R7WqfU7UFYgq06tSpo/POO08HDx60qj0RzZ07V+3bt1d6erp69uyptWvXhj32tdde06BBg3TOOecoIyNDOTk5euedd0xri9kfwEgfFklq3LC+Kg3D9Itkt63xrY0ZeUenr999cmQPLR7bR+umDnR1kBUQKVA8u1GaI68l+WAAkDqS5Tvf6wPN0QzAnr5Sx+rXLZpr2dral2xiztF65JFHdN9992nevHnq1q2bFW0Ka+nSpZo0aZLmzp2rfv36acGCBbrqqqu0bds2tW3btsbx77//vgYNGqTf/e53aty4sV566SVdd9112rBhgy6++OKE22P2vku1LWsLCMymmT2d7fSyxWiZObVv1iaDTqhtGvyNwr1RPYaZr6VXl1wAACKrnrt06OgJ/eat5PjOd0N+fKKi3ePUjnN1NNey4dqXjOLaR+v777/XyZMn1aBBAzVs2DDo999+a92Mx2WXXaZLLrlE8+bNq7qtc+fOuv766zVr1qyoHqNr16666aabot5YOVArv82kP6pOWvgcNDP3XXJi7wEv7C/lxn0i3Mju15LXBQCSVzTXJJJ3v/PfKNyru5cURjzuyZE9NLxHK+sblIDainnYfa4O9b7JykjTqN5tld2skWeKjYQSyz5aMc9ozZkzJ952JeTEiRPatGmT7r///qDbBw8erPXr10f1GJWVlTpy5IiaNAm/DvT48eM6fvx41c9lZWVRPbaZswOB2YqCHQd11x826/AP5TWOMXs2ze3V+CJN7ZvZF15n52vJ6wIAySvcxXkoXv3Od1sxiUSEW6njxLnayQIUbhJzoDVmzBgr2hHRgQMHVFFRoczMzKDbMzMzVVJSEtVjPPbYYzp69KhuvPHGsMfMmjVLM2fOjLl9Zn8A69bxqU4dX8ggK8DM6Wy3V+NLhql9u9j5WvK6AEByinU/JMmb3/luH2g2g1Pnai+naJgl5g2Ld+3aVes/q/l8wReHhmHUuC2UxYsXa8aMGVq6dGmt1RGnTZum0tLSqn+7d++uvT2yrlqK3XlTbq7G55UcMrew67XkdQGA5BTvfkiSt77zYy0m4UWcq50T84xWdnZ2rYFNRUVFQg0Kp1mzZqpbt26N2av9+/fXmOWqbunSpbrtttv0pz/9SVdeeWWtx6alpSktLa3G7b7/+2fnTI8T09lOT/WGW1+cTFP7kj2bItrxWibb6wIAOCWRi267v/MTPadGW0zCqzhXOyfmQGvLli1BP5eXl2vLli16/PHH9dBDD5nWsOoaNGignj17Ki8vTyNGjKi6PS8vT8OHDw97v8WLF+vWW2/V4sWLdc0118T9/I/f1F3/b/VuWz+ATk1nOzXVW1s1nEFdspJmat/OCn1Wv5apsOQCAFJRPBfdTnznm3VOdXqg2SoVlYYqKw01blg/bDoK52rrxFx1MJy33npLjz76qN577z0zHi6kpUuXKjc3V/Pnz1dOTo6ee+45LVy4UJ9++qnatWunadOmae/evXrllVcknQqyRo8erSeffFI33HBD1eM0bNhQfr8/quc8vbJIozPPsv0DGEhElULPptmxpM+O2ZdoquFIcrwvEpWMFfrc8B4FAJirotJQ/9mrwg6kVefEd34ynlPN5EQV61QQS9VB0wKtL774Qj169NDRo0fNeLiw5s6dq0ceeUTFxcXq1q2bnnjiCf3Lv/yLJOmWW27RV199VRXsXXHFFVqzZk2NxxgzZoxefvnlqJ4vls60ipN7FNnx3IEv83BfBIGRlnVTBypvW4ln92uK5e/02gga+2gBQPIJN5AWit3f+cl8TjVDtBUjOVfHztJAq3q5c8MwVFxcrBkzZuizzz5TYWFhzA12MzcEWpI9s0rV2TVSFOu+T070hRm8sFdZIrz6ugAAwgs3kDb9ms46u1GaY9/5yX5OTUSkIFSSGjesr2d/eon6dGjKuTpGlu6j1bhx45CV/9q0aaMlS5bE+nCIkt15U3buuRBrNRyvlgtN9qo/Xn1dAADhuTV3KdnPqYmIpmLk4R/KVcfnc/x1THYxB1qrV68O+rlOnTo655xzdO6556pevZgfDi5l554LqVINJ1X+TgBAcnHjQBrn1PAIQt0j5showIABVrQDLmPnhzRVKtelyt8JAIDV3HpOdcMyeoJQ94hrCmrHjh2aM2eOioqK5PP51LlzZ919993q2LGj2e2DQ+z8kAY2C/z5os2271Vmp1T5OwEAsJobz6nRFoayOhhzaxCaiurEeod33nlHXbp00caNG3XRRRepW7du2rBhg7p27aq8vDwr2ggHBD6k4T72Pp368jDrQxrYLDDLHxy4ZfnTk6rkaKr8nQAAWM1N59RAAbHqaRclpcf080WbtWJrcdVx/Wev0qiFBbp7SaFGLSxQ/9mrqn5vhkAQKqnGdRwDu/aKuergxRdfrCFDhujhhx8Ouv3+++/XypUrtXnzZlMb6DS3VB10ghP7I7lhyt0OqfJ3AgAQi3jOj06fU6MtNT/9mi666w/27fvF1ivWsLS8e3p6uj755BOdd955Qbd//vnnuuiii3TsWHIl1qVyoCXxIUVo0ZzUnD7xAQC8xavXHNGWmm/SqIG+PXoi5O+s2veLc7H5LC3vfs4556iwsLBGoFVYWKjmzZvH+nBwObeWdYVzojkRevVkCQBwRri9OwNL79y8vD7awmDhgizJ3GrOp3NjxchUEnOgNXbsWI0bN05ffvml+vbtK5/Pp3Xr1mn27Nm65557rGgjHMaHFAHRnAglefZkCQCwn517d1rBzOp9lFxPLjEHWtOnT9dZZ52lxx57TNOmTZMktWzZUjNmzNDEiRNNbyBgF5bD1S6aE+GM5Z9K8nn2ZAkAsJ+de3daIZoqf2c3qq9vj5ZHfCxKrieXmAMtn8+nyZMna/LkyTpy5Igk6ayzzjK9YYBkX2DDcrjIojkRlpQdr/Ux3H6yBADYz+sb7EZTav63w7vpN28VUXI9xcS1j1YAARasZFdgw3K46Jh5gnPryRIAYL9k2GA3UGq++nVL1mnXLXXq+Fy17xesF3OgtW/fPt177736+9//rv3796t60cKKigrTGofUZVdSLMvhomfmCc7NJ0sAgL2SZYPdSAXEognGkFxiDrRuueUW7dq1S9OnT1eLFi3k8yXvhSWcYWdSLMvhohfNiTAzI02ST/vK7D1ZpnLuHAB4XTRL77wy2xOpgBjVnFNLzIHWunXrtHbtWvXo0cOC5gD2JsWm4nK4eIOSaE6EM4Z1lSRbT5apnjsHIDUk+4BSMs/2hHrtknlgFj+KOdBq06ZNjeWCgJnsTIpNteVwiQYl0Z4I7TpZennfFQCIVrjv7unXdNbZjdKSJvhKxtkeBgNTm8+IMWpauXKlHnvsMS1YsEDZ2dkWNcs9Ytn9GeaIdof1xWP7JDwiVFFpqP/sVaYshzN7N3ezhQtKAi2OJShxQyn8wGsXbvbTK68LANQm3Hd3KFzAu4uZ5124RyyxQVQzWmeffXZQLtbRo0fVsWNHnXHGGapfv37Qsd9++20cTQZ+ZGdSrFuXw5nN7Ly3aDaxtnqja6/vuwIAkdT23R0Ks/nu4fVNmGGOqAKtOXPmWNwM4Ed2J8W6bTmcFZIxKPH6visAEEmk7+7quIB3j2Q87yJ2UQVaY8aMsbodQBC7k2KjWRfu5bXjyRiUJMO+KwBQm3i+k7mAd4dkPO8idlEFWmVlZTHlJx05coTNjJEwuwMbNyyHs0oyBiXJsu8KAISTyHcyF/C1szqPOBnPu4hd1DlaxcXFat68eVQP2qpVKxUWFqpDhw4JNQ7wamDjNskYlCTTvisAEEqk7+7a2H0B76Xy83ZUAkzG8y5iF1WgZRiGnn/+eZ155plRPWh5eXlCjQJgrmQNSpJ53xUAqO27OxwnLuC9VMLcrm1BkvW8i9hEVd49Ozs7qOpgNN5//321adMm7oa5BeXdkUy8dDKMhZdGUgEgVqG+u0Nxomy4l0qYO7EtSLKed1NZLLFBzPtopRoCLSQbghIA8J7q392Hjp7Qb95y9gLea/sZ2rlP5+k47yYX0/fRApA8yHsDAO8J9d09pJs1BaOiDQy8VsLcqUqAnHdTF4EW4AKMdgEAYmXFBXwsS93cWsI83DmVSoCwG4EWEIZdwQ/rtwEAbhBroQg3Bi61nVMHdcmiEiBsRY5WBORopSa7gh8vJRF7CTOEABCbePKtAveJFLjYlaMVzTlVkn6+aLOk0JUAOe8iklhigzrRPmhhYWGi7QI8IfBFXf1kExjRW7G12JTnqag0NPPNbSFPToHbZr65TRWVjIXEYsXWYvWfvUqjFhbo7iWFGrWwQP1nrzLtdQOAZBRLvlVAoIS59GOgEmB3CfNoz6mDumRp3s8uUZY/eJYty59OkAXTRR1oXXLJJerZs6fmzZun0tJSK9sEOMbO4CeekxpqZ1eQDADJJt58q8B+hk4HLrGcU4d2a6F1Uwdq8dg+enJkDy0e20frpg4kyILpos7R+uCDD/Tiiy/q/vvv1z333KMbbrhBt912m37yk59Y2T7AVnZWUHJrErFXRQqSffpxNJNlhAAQLJF8q6HdWmhQF2sqIEYr1nMqlQBhh6hntHJycrRw4UKVlJRo3rx52rNnj6688kp17NhRDz30kPbs2WNlO6vMnTtX7du3V3p6unr27Km1a9fWevyaNWvUs2dPpaenq0OHDpo/f74t7YQ32Rn8xHpSq6g0lL/joN4o3Kv8HQdTdklhuH5ghhAA4te7fRO18KfXWAIY4NOpXOVwhSICgcvwHq2U07Gp7QNabizMAcRcdbBhw4YaM2aMxowZox07duill17SggULNGPGDA0aNEhvv/22Fe2UJC1dulSTJk3S3Llz1a9fPy1YsEBXXXWVtm3bprZt29Y4fufOnbr66qs1duxYLVq0SB988IHuvPNOnXPOOfq3f/s3y9oJ77LzizpwUoum+hGVCU+prR+On6yM6jGYIQSAmgL5Vj9ftFk+hS4UYVe+VTxiOacCdkm46uB3332nV199Vb/85S91+PBhVVRUmNW2Gi677DJdcsklmjdvXtVtnTt31vXXX69Zs2bVOH7q1Klavny5ioqKqm4bP368/vGPfyg/Pz+q56TqYGqxu4JSIKdICl/9SFLSVCZMpBpgpGpSk648T0+8+0XEx1k8tg/LRQAgDC8P7EVzTnX73wD3iyU2iHsfrTVr1ujFF1/UX/7yF9WtW1c33nijbrvttngfLqITJ05o06ZNuv/++4NuHzx4sNavXx/yPvn5+Ro8eHDQbUOGDNELL7yg8vJy1a9fv8Z9jh8/ruPHj1f9XFZWZkLr4RV2j+gFkoirn9SyTtvzo//sVUmRd5TIyTua/KvFG3cpKyNd+8oYzQSAeLkh3ypekc6pBFmwW0yB1u7du/Xyyy/r5Zdf1s6dO9W3b189/fTTuvHGG9WoUSOr2ihJOnDggCoqKpSZmRl0e2ZmpkpKSkLep6SkJOTxJ0+e1IEDB9SiRc0P3KxZszRz5kzzGg7PsfuLuraTWv6Og7YV57BSrJtgVhdN/lVJ2XFNvvJ8zXn3c08uewEAt/ByoYhw51RJyt9x0HPBI7wt6kBr0KBBWr16tc455xyNHj1at956qy644AIr2xaSzxf8oTAMo8ZtkY4PdXvAtGnTNGXKlKqfy8rK1KZNm3ibC4+ye0Qv3EktGSoTmlENMNq/L7vZGbaPZrI5MpA6+Lx7Q/VzqpeXQ8Lbog60GjZsqL/85S+69tprVbduXSvbFFKzZs1Ut27dGrNX+/fvrzFrFZCVlRXy+Hr16qlp09AjNWlpaUpLSzOn0fA0N4zoJUMVJTNK5sfSDzkdm9oWJHPyBlIHn3dvSnRFBZCIqMu7v/XWW8rJyXEkyJKkBg0aqGfPnsrLywu6PS8vT3379g15n5ycnBrHr1y5Ur169QqZnwW4TaLldt3AjFm5WPvBjjLDbI4MpA4+794UaUWFdGpFRapulwLrRR1oJVic0BRTpkzR888/rxdffFFFRUWaPHmydu3apfHjx0s6texv9OjRVcePHz9eX3/9taZMmaKioiK9+OKLeuGFF3Tvvfc69ScAMQkU55BUI8jwSt6RGbNybusHTt5A6oj0eTck/XLZJ1q2JbX3OIyH1ftDsr8inBZ31UEn3HTTTTp48KAefPBBFRcXq1u3bnr77bfVrl07SVJxcbF27dpVdXz79u319ttva/LkyXr22WfVsmVLPfXUU+yhBU/xehUls/Y2cVM/mLEcEoA3RPq8S9K3R8s1eWmhJJYTRsuOpZjJkOcMb4t6H606dero97//vfx+f63HDRs2zJSGuQX7aMEtvJyEbebeJm7ohzcK9+ruJYURj3tyZA8N79HK+gYBsEy0n/cA9myKLNK+iGb1Xf6Ogxq1sCDiceyviFhYto/WmDFjav29z+ezdMNiIJW5oThHvMycjXJDPyRDkRIA0Yn1c+y1PQ7tZkYl2miZtaICiFdMgVZJSYmaN29uVVsAJDEvb4JZHSdvIHVE+ryHwvLh8Oxceh3I7/35os3srwhHRF0Mo7a9qgAgGnZUA7SD24pzALBObZ/3SBLJ/bG6UIRT7M6bCqyoyPIHz0xm+dNZ3gnLRT2j5YaqgwDgFm4qzgHAWuE+75HEu3w4mffscmLpdTKtqIC3RB1ojRkzRg0bNrSyLQAc4IbiEl7FyRtIHad/3ktKf9Bv3irSoaMnTF8+nMgGu174Pndq6bUb8nuReqKuOli3bl0VFxenXI4WVQfhBLtOlsk8auoUL1zoAEicmdVUAyoqDfWfvSrsrFkgCFk3dWCN7xUvfZ9b0XeAXWKJDWIq756KxTAItGA3u06WdpXXTSVeutABkDizP/PxliP34vc535fwKssCrX379umcc84xpZFeQaAFO9l1skxk1BShefFCB0DizJzFjmePPi9/n7MCAF5k2T5a559/fsTqg99++20sDwng/9i5t4id5XVTgZ2vHQB3MTP3J55CEV7+PidvCskupkBr5syZ8vv9VrUFSGl2niztLq+b7Lx8oQPAPeIpFOHW73Nmq4AYA62RI0emXI4WYBc7T5ZOlNd1mpUnfbde6ADwlng22HXj9zn5V8ApUQdabFgMWMvOk6VT5XWdYvVJ340XOgC8KdY9+tz2fZ5IeXog2bBhMeASdp4s4xk19So7Tvpuu9AB4G2x7NHnpu9z8lWBYHWiPbCyspJlg4CFAidL6ceTY4AVJ8vAqGmWP3iWJcufnjQjjpFO+tKpk35FZWIDSXa/dkAqq6g0lL/joN4o3Kv8HQcT/vy6VaBQxPAerZTTsWmt3x9u+T6PJV8VSAUx5WgBsFasS0bMeL5oR029yM4iFXa/dkAqIvcnPDd8n5OvCgQj0AJcxu6TZTTldb1aPcrsk36kfrD7tfPq6wLEg9yfyJwul06+KhCMQAtwIadPlqfz8giymSf9aPvBrtfOy68LECtyf7yBfFUgWNQ5WgBST2AEufryu8AI8oqtxQ61LDqBk364yy6fTgUnkU76busHt7UHsBq5P9YxM+ctUr6qIWnkpW3014//mdT5dUAAM1oAQkqGEWQzqnG5rR/c1h7ADuT+WMOKmfFw+ar+M+pLkp549wvTngtwO2a0AISULCPIiVbjcls/uK09gB2iXQb8xb7vmCmJkpUz40O7tdC6qQO1eGwfPTmyhyZfeb5Kvy/X4e/LTX8uwM2Y0QIQUjKNICdSpMJt/eC29gB2iJT7E/DM6v/VM6v/17SZkmQtOGPHzHggX7Wi0lD/2auYhUdKItACEFKyVY+Kt0iF2/rBbe0B7FDbMuBQzKhEmMwFZ+zc+sLO5wLchqWDAEIyq5CE17mtHyK1R5IaN6yvSsNg+RSSSrhlwKEkuiF5shecsXNmnFl4pDICLQAhRaoeJUUuJJEM3NYPtbUn4PAP5frp8xvUf/Yqz18QAqc7Pfdnwk861npsvPmKkZbVSfEHcG5h58w4s/BIZQRaAMJKtJBEsnBbP0Q7sp8so+/A6QLLgM/LPCuq42OdKUmFgjN2ztS7bVUAYCdytADUKpFCEpF4KdHcyn5IpD0FOw7qrj9s1uEfymscQ6I5kplVMyWpsNTNjK0v3PhcgNsQaAGIKN5CErXxYqK5Ff2QiLp1fKpTxxcyyAog0RzJKlIlQp9OzTrHOlOSKkvdwu13lWXB97CdzwW4CYEWANsFEs2rXxyZUSks1aTC6DsQilUzJVYFcG5k50y921YFAHYg0AJgKzv2b3EbK5dIpsroO9zFLct+rZgpSbWlbnbO1LttVQBgNQItAEGsvoBKtT1VrF4imUqj73AHty37tWKmhKVuAMxAoAWgih0XUKm01M2OJZKpNvoOZ7l12a8VMyXxBnBume0D4DwCLQCS7LuASpWlbnYukWT0HXZIxWW/sQZwbpvtA+Asz+yjdejQIeXm5srv98vv9ys3N1eHDx8Oe3x5ebmmTp2qCy+8UI0aNVLLli01evRo/fOf/7Sv0YBH2LlBZ6rsqWL3Xjynb+T65MgeWjy2j9ZNHWjZxV1FpaH8HQf1RuFe5e846OnNWxGdVNhfKhGBwarqfcR+dkDq8syM1s0336w9e/ZoxYoVkqRx48YpNzdXb775Zsjjv//+e23evFnTp09X9+7ddejQIU2aNEnDhg3TRx99ZGfTAdezM28qVZa6ObFE0q5Ec0btU1MqLfuNVSrO9gGIzBOBVlFRkVasWKGCggJddtllkqSFCxcqJydH27dv1wUXXFDjPn6/X3l5eUG3Pf300+rdu7d27dqltm3bhnyu48eP6/jx41U/l5WVmfiXAO5k9wVUKix1S9YlkuGWmBaXHtP4RZt1W79sXdkli7yUJJSs72kzpFqRHwDR8USglZ+fL7/fXxVkSVKfPn3k9/u1fv36kIFWKKWlpfL5fGrcuHHYY2bNmqWZM2cm2mTAU5y4gEr2PVWSsRpgbaP2AS988JVe+OArZriSkFPvaS8Ul2C2D0AonsjRKikpUfPmzWvc3rx5c5WUlET1GMeOHdP999+vm2++WRkZGWGPmzZtmkpLS6v+7d69O+52A17hVN5UYKnb8B6tlNOxqesunhIRWCIpqUa/enWJZKRR+9ORl5J8nHhPr9harP6zV2nUwgLdvaRQoxYWqP/sVa57X7l1to9cSsBZjgZaM2bMkM/nq/VfIJ/K56v5xW0YRsjbqysvL9fIkSNVWVmpuXPn1npsWlqaMjIygv4ByS4ZgwI3CCyRzPIHX1xl+dMdK4OdiFhG480uogJ3qO09/ezNF8vfsIFpF/VeKi6R6GCVFQGRV4JUIJk5unRwwoQJGjlyZK3HZGdn6+OPP9a+fftq/O6bb75RZmZmrfcvLy/XjTfeqJ07d2rVqlUETkAYqZA35QSzlki6YflUrKPx5KUkp1Dv6UNHT+g3b5lXIMVrxSUSKfJjRXEZt+53BqQan2EYrh9qLCoqUpcuXbRhwwb17t1bkrRhwwb16dNHn332WdgcrUCQ9cUXX2j16tU655xzYn7usrIy+f1+lZaWEqQhJbjhgh7B3FLlr6LSUP/Zq8Lm6ITz5MgeGt6jlWXtSlVu+ayGu6gPtCSei/r8HQc1amFBxOMWj+3jqiA+1s+qFX0X+JyGW+YbyKVbN3Ug3+1AHGKJDTxRDKNz584aOnSoxo4dqwULFkg6Vd792muvDQqyOnXqpFmzZmnEiBE6efKk/v3f/12bN2/WX//6V1VUVFTlczVp0kQNGjRw5G8B3M6uEuHRcsvFpFPcNDJd26h9bVKxCp3V3BR8WzHz5NXiErHMYFvVd1RABNzDE4GWJL366quaOHGiBg8eLEkaNmyYnnnmmaBjtm/frtLSUknSnj17tHz5cklSjx49go5bvXq1rrjiCsvbDCAxbrmYdIobl0+FW2IaihcrK3qBm4Jvqy7q3VpcIhrRDlZZ1XfRBp9/+79crVQbvALs5JlAq0mTJlq0aFGtx5y+CjI7O1seWBUJIAw3XUw6xa0j06eP2udtK9GLH3zl6ObTqTTr6bbg26qZp2TcHqE6q/ou2uDzlfyv9Ur+1yk1eAXYzTOBFoDU4baLSaeYcSFmVRASGLXP6dhUvds3cayIitOznnYHeW4Lvq2aeUqkuIRXWNV3kYLU6lJp8AqwG4EWANdx28WkUxK9ELMrCLF78+lAcBOYTavOrgtHJ4I8t+UuWTnzlOyVUK3qu1hzKVNp8AqwG4EWANdx28WkUxK5ELN76aVdRVRCBTfV2XHh6NTSVqdyl8LN3Fk982R3EG8nK/sullxKKXUGrwC7ObphMQCEYsfFpBUbhJot3o2kIy29lLy5kXC4DWxDOf3C0WxO9m+iG+PGI9LGt1ZvzB0I4of3aKWcjk2TIsgKsLLvhnZroXVTB2rx2D4andMuqvsk++AVYDdmtAC4jtWJ8E7n9cQinuVTybj0srbgpjZWXDg62b925y5FO3OXzDNPVrOy706faX4l/+uIx7uxiiPgZQRaAFzHyotJL1YzjPVCLBmXXkYKbsKx4sLR6f61K3cp1qI0btuDz0us7rtUqOIIuBGBFgBXsuJi0svVDGO5EPPyHkThxBq0WHnh6Ib+tWMGKRlnRlNVKlRxBNyIQAuAa5l9MZkqF47JOHodS9Bi5YVjRaWhykpDjRvW1+EfysM+vx39Gyr4NrPcvNMzdzBXsldxBNyIQAuAq5m5pCZVLhyTcfQ6lr2BrLpwjKbiodVBXm1BlNm5h26YuYO5yKUD7EWgBcA1rN78NZUuHJNt9DqavYFu65etK7tkWXLhGC63rzo7g7zTgygrcg+9PDNq90bSXkIuHWAfn2EY3qrva7OysjL5/X6VlpYqIyPD6eYAScuOSoAVlYb6z14V8cJx3dSBSXNR5uQFpxXP7UTFyMD7praZrMYN6+vZn16iPh3MLz8eLogKPMuzN1+s37xVFLZ9ibyvA88thZ4ZdWPxGC9VFQXgPbHEBgRaERBoAdaLdCFp5sWcFy8cvcjKi127g8f8HQc1amFBxOMWj+1j+kxBpCDPJ+nsRvX17dHQ+WJmtM9LgYud3yUAUlMssQFLBwE4yu5KgMm2pM6NrC6hb/fSJydz+6Ip4BJNkCXF3z6v5PUk+l3CckMAZiPQAuAoJyoBeuXC0Yu8XEI/HCdz+8wM3hJpnxfyehL5LvHSrB0A7yDQAuAop2YLvHDh6BaxjPQnYwn9WItCmDkzEm1w1KRRAx06esJzRSvMFO93iRc3MQfgDQRaAByVSpUAvSjWkX63lNA3M9iJpVy+2TMj0QZ506/porv+kDzl/OMRz3dJMs7AAnCPOk43AEBqC1xIhruE8enUhWqyj8a7UWCkv/oMVWCkf8XW4hr3cUPgvGJrsfrPXqVRCwt095JCjVpYoP6zV4Vsb7QCuX1Z/uB2Z/nTq2Y84umvSAJBnqQan5HTg6irL4rcvmQXz3dJLDOwABArZrQAOCoZN9dNBvGO9Du995KVy8Bqy+2zcmYk2gIuqZ57GM93iVtmYAEkJwItAI6jEqD7xJtr5WTgbMcysHC5fVbnpkUbRKV67mGs3yVumIEFkLwItAC4QqqPxrtNIiP9TgXOThbisGNmJNWDqGjF8l1ixwwsZeOB1EWgBcA1uJB0j0RH+p0InJ1cBsbMiLtE+11i9QwsZeOB1EYxDACwWUWlofwdB/VG4V7l7zioispQY+nOMqNISeBid3iPVsrp2NTyUXwngx2KunhXNIVO4mFFcRQA3sKMFgCEYcWSH6+McHuxSImThTi82F/4kdkzsJHyBSXpl8s+0Q/llcrKYDkhkKx8hmG4byjVRcrKyuT3+1VaWqqMjAynmwPAJlYEROEq4gUur9xYhtsrgWFAoI+l0MGO1X3stf6CNfJ3HNSohQVRH897BPCOWGIDAq0ICLSA1GNFQFRRaaj/7FVhizUEZlvWTR3oupFtryXzOx3seK2/YL43Cvfq7iWFUR/v5sEWAMFiiQ1YOggAp7GqRLiTFfES5bUiJU5XsPRaf8F8seYBmrX9AAB3IdACgNNYFRCxMaq9CHa8KVlmAyPlC4bi5sEWAPEh0AKA01gVEFH+G6id00s+zVRbcZRIGGwBkgfl3QHgNFYFRJT/jo8XSuEjtFheu2QshR6ubHwkDLYAyYMZLQA4jVUlwin/HbtkmuFINbG8dlblRbrB6fmCJaU/6DdvFenQ0RO2bz8AwBnMaAHAaQIBkaQas0+JBkRWbYyajJJxhiNVxPraxZIX6UWBfMERl7TW70Z0k2T+dwsAd/JMoHXo0CHl5ubK7/fL7/crNzdXhw8fjvr+d9xxh3w+n+bMmWNZGwEkBysDoqHdWmjd1IFaPLaPnhzZQ4vH9tG6qQMJsk4TzWavM9/cxjJCF4rntUulQjEMtgCpxTNLB2+++Wbt2bNHK1askCSNGzdOubm5evPNNyPe9/XXX9eGDRvUsmVLq5sJwKOqVzsb1CXLshLhVMSrnZdL4ae6eF67VCsU4/T2AwDs44lAq6ioSCtWrFBBQYEuu+wySdLChQuVk5Oj7du364ILLgh7371792rChAl65513dM0119jVZAAeQi6Qu6TSDEeyiee1syov0s0YbAFSgyeWDubn58vv91cFWZLUp08f+f1+rV+/Puz9KisrlZubq/vuu09du3aN6rmOHz+usrKyoH8Akhe5QO6TajMcySSe187KvEgAcJInAq2SkhI1b968xu3NmzdXSUlJ2PvNnj1b9erV08SJE6N+rlmzZlXlgfn9frVp0yauNgNwP3KB3IlS+M5KpKR+vK8duUsAkpGjSwdnzJihmTNn1nrMhx9+KEny+Wp+bRuGEfJ2Sdq0aZOefPJJbd68OewxoUybNk1Tpkyp+rmsrIxgC0hSqZgLVD0XzY25IW4phe+FvjJbostoE3ntyF0CkGwcDbQmTJigkSNH1npMdna2Pv74Y+3bt6/G77755htlZmaGvN/atWu1f/9+tW3btuq2iooK3XPPPZozZ46++uqrkPdLS0tTWlpa9H8EAM9KtVwgL+WiBWY4qrc3y6b2eqmvzBJYRlt9/iqwjDbamaVEXjtylwAkE59hGK5fE1NUVKQuXbpow4YN6t27tyRpw4YN6tOnjz777LOQxTAOHjyo4uLg3IohQ4YoNzdX//mf/1lrAY3TlZWVye/3q7S0VBkZGYn/MQBcI3/HQY1aWBDxuMVj+3j+4i/cRXRgrsCty7OcmFXyal8loqLSUP/Zq8LO8AYKUqybOjDq/k/FGUEAyS+W2MATVQc7d+6soUOHauzYsVqwYIGkU+Xdr7322qCAqVOnTpo1a5ZGjBihpk2bqmnT4Auj+vXrKysrK+ogC0ByS5VqZ5Fy0Xw6lYs2qEuW6y6E7Z7hcEtf2R2kWLGMltkpAKnOE4GWJL366quaOHGiBg8eLEkaNmyYnnnmmaBjtm/frtLSUieaB8CD3JILZLVUzEWLlxv6yolli6m2jBYA7OCZQKtJkyZatGhRrcdEWgUZLi8LQOpyOhfIDlxER8/pvjIrTyog2pkxSuoDgPk8E2gBgFWSvdoZF9HRc7KvzF62GMvMWLIvoyVfDIATCLQAQMmdT5LsF9FmcrKvzFy2GOvMWDIvow0XcE6/prPObpRG8AXAMp7YsBgAUk0im8ZWF7iIllRjI1mvX0Sbzcm+MmvZYrwbcSfjpsGBgLN6AFtcekx3/mGLRi0s0N1LCjVqYYH6z16lFVuLwzwSAMSOGS0AcBkriiGkQi6aWZzqK7OWLSYyM5ZMy2hrCzhDiTcPDgDCIdACABcxuxjC6Zy+iPZSnowTfWXWssVEZ8aSZRltpICzOrdvdQDAewi0AMAl7NjDyamLaCdKlifK7r4yK0+K4ienxFMZkq0OAJiJHC0AcIlYlnx5Sbg8mcAsHXkxPzIjTyowMxYuHPPpVJCb7MVPEgkk2eoAgBmY0QIAl3B6Dycr2DFL53axLplMdNliMlcQjEWkpZi1SfbZPgD2INACAJdIxiVfZpYs96J4l0wmumyR4ie1B5zhsNUBADMRaAGASyTjfldumaVzohCHlYVNouF08RM3CBdwhpJKs30A7EGgBQAOOz0IGHlpW8159/OkWfLlhlk6JwpxuGXJZLJUEExEqIDz0NET+s1bqTvbB8AeBFoA4KBQQUDjM+pLkg5/X151m1cvAp2epXNqVinVl0y6TaiAc0i31J7tA2A9Ai0AcEi4IKD0+3IZkiZfeZ6ymzXy9EWgk4UZnJxVcsuSSYTHbB8Aq1HeHQAcEE0QsOTD3br2opbK6djUk0FWgBkly+PhZLl8NyyZBAA4ixktAHBAqi0tc6Iwg5OzSk4vmQQAOI9ACwAckIpLy+xequXkrBJ7WQEAWDoIAA5gaZn1ArNK4UIZn05VH7RqVsmpJZMAAHdgRgsAHMDSMuu5YVaJvawAIHUxowUADggEAZJqzLiwtMw8bphVCiyZHN6jlecLmySiotJQ/o6DeqNwr/J3HFRFZaghBgBIHj7DMPimq0VZWZn8fr9KS0uVkZHhdHMAJBknNtON5PQNlJNlBiYZ/yYvCfc+n35NZ53dKI3XBYBnxBIbEGhFQKAFwGpuCgLcGPjB28LtFxcK7zUAbkegZSICLQCpItwFcSDko4ADYlVRaaj/7FW1bmVwOt5rANwultiAHC0AQMQNlCVp5pvbyKtBTCLtF1cd7zUAyYRACwAQ0wbKQLTi2QeO9xqAZEGgBQBIyQ2UYb1E9oHjvQbA6wi0AABsoAxLRNo0uja81wB4HYEWAKSwwN5GJaU/qEmjBmEviH06VRGODZQRi9r2iwuH9xqAZFHP6QYAAJwRqpR7KGygjEQENo3mvQYg1RBoAUAKimVvoyz2NkKChnZroUFdsoL2izt09IR+81Zw8MV7DUAyIdACgBRTWyn3gCaN6mv6tV2VleHsBspIHnXr+JTTsWnQbUO6Zblms24AMBuBFgCkmGj2Nvr2aLmyMtJrXBgDZgoVfAFAsiDQAoAU4+ZS7hWVBjMcNqPPAcAangm0Dh06pIkTJ2r58uWSpGHDhunpp59W48aNa71fUVGRpk6dqjVr1qiyslJdu3bVH//4R7Vt29aGVgOA+7i1lHuo4hwtyNkxVfWgKlSeFH0OAObwGYYRTS6046666irt2bNHzz33nCRp3Lhxys7O1ptvvhn2Pjt27FDv3r112223adSoUfL7/SoqKtKll16q5s2bR/W8ZWVl8vv9Ki0tVUZGhil/CwA4qaLSUP/Zq1RSeixknpZPp4oSrJs60LaZjXDFOQLPPu9nl3Dhn6BYq0zS5wBQUyyxgScCraKiInXp0kUFBQW67LLLJEkFBQXKycnRZ599pgsuuCDk/UaOHKn69evrf/7nf+J+bgItAMkoENhICgpunLjIDgR+4QIAJwK/ZBNLlUmJPgeAcGKJDTyxYXF+fr78fn9VkCVJffr0kd/v1/r160Pep7KyUm+99ZbOP/98DRkyRM2bN9dll12m119/vdbnOn78uMrKyoL+AUCyCextlOUPXh6Y5U+3fSYjUnEOQ1Jx6TFt3PmtbW1KJtFUmayOPgeAxHkiR6ukpCTkUr/mzZurpKQk5H3279+v7777Tg8//LB++9vfavbs2VqxYoVuuOEGrV69WgMGDAh5v1mzZmnmzJmmth8A3CjU3kZ2FkII5Av9bWtxVMc7UZwjGURTZTIc+hwA4udooDVjxoyIQc2HH34oSfL5ap74DcMIebt0akZLkoYPH67JkydLknr06KH169dr/vz5YQOtadOmacqUKVU/l5WVqU2bNpH/GADwIKfKa0ebL3Q6u4tz2MnKyn+JBEvJ3OcAYDVHA60JEyZo5MiRtR6TnZ2tjz/+WPv27avxu2+++UaZmZkh79esWTPVq1dPXbp0Cbq9c+fOWrduXdjnS0tLU1paWhStBwDEI958od7tm1jZLMfKnFtdbTGeYMmuPgeAZOZooNWsWTM1a9Ys4nE5OTkqLS3Vxo0b1bt3b0nShg0bVFpaqr59+4a8T4MGDXTppZdq+/btQbd//vnnateuXeKNBwDELNZ8oUCY88B1XSwNepwqLR8u6CwpPaafL9psSr5c7/ZN1MKfHrbKZHV29TkAJDtPFMPo3Lmzhg4dqrFjx6qgoEAFBQUaO3asrr322qCKg506ddKyZcuqfr7vvvu0dOlSLVy4UP/7v/+rZ555Rm+++abuvPNOJ/4MAEh5seYL2VGcIxDsVG9XINhZEWUOWbQqKg3l7zioZZv36JfLtoYMfgK3zXxzmyoqEysOXLeOTw9cd2p1RzRhkxMFUQAgGXmiGIYkvfrqq5o4caIGDx4s6dSGxc8880zQMdu3b1dpaWnVzyNGjND8+fM1a9YsTZw4URdccIH+8pe/qH///ra2HQBwSrT5QqNz2umqbi0sW74XWCZYUvqDfvNWUdhgx6dTwc6gLlmmtCOW3LTTK/8lmkcXqDIZatZu+jWddXajNEcKogBAMvNMoNWkSRMtWrSo1mNCbQl266236tZbb7WqWQCAGESbL3RVtxaWFelwKtiJNTctwKzKf05XmQSAVOOZQAsA4H2R8oWsLsLgVLATz15WAWZW/nOqyiQApCJP5GgBAJJDbflCVhdhcDLYiWcvK59OLe2j8h8AeBOBFgDAVoF8oSx/cPBidREGJ4OdWGfEqPwHAN7H0kEAgO2cyBdyMtiJdUYsy4bS8gAAaxFoAQAcYXe+kJPBTjR7WTVpVF/Tr+2qrAyKVABAMiDQAgCkBCeDnUBu2s8XbZZPCnr+wDP8bsSFzGABQBIhRwsAkBIiFeLw6VSwM+LiVsrp2NT0GSWnctMAAM7wGaE2n0KVsrIy+f1+lZaWKiMjw+nmAAASFGofrRY25kQFNktmLysA8J5YYgMCrQgItAAg+RDsAADiEUtsQI4WACDlsHEvAMBq5GgBAAAAgMkItAAAAADAZARaAAAAAGAyAi0AAAAAMBmBFgAAAACYjEALAAAAAExGoAUAAAAAJiPQAgAAAACTEWgBAAAAgMkItAAAAADAZARaAAAAAGAyAi0AAAAAMBmBFgAAAACYjEALAAAAAExGoAUAAAAAJiPQAgAAAACTEWgBAAAAgMkItAAAAADAZARaAAAAAGAyAi0AAAAAMBmBFgAAAACYjEALAAAAAEzmmUDr0KFDys3Nld/vl9/vV25urg4fPlzrfb777jtNmDBBrVu3VsOGDdW5c2fNmzfPngYDAAAASFmeCbRuvvlmFRYWasWKFVqxYoUKCwuVm5tb630mT56sFStWaNGiRSoqKtLkyZP1i1/8Qm+88YZNrQYAAACQijwRaBUVFWnFihV6/vnnlZOTo5ycHC1cuFB//etftX379rD3y8/P15gxY3TFFVcoOztb48aNU/fu3fXRRx/Z2HoAAAAAqaae0w2IRn5+vvx+vy677LKq2/r06SO/36/169frggsuCHm//v37a/ny5br11lvVsmVLvffee/r888/15JNPhn2u48eP6/jx41U/l5aWSpLKyspM+msAAAAAeFEgJjAMI+Kxngi0SkpK1Lx58xq3N2/eXCUlJWHv99RTT2ns2LFq3bq16tWrpzp16uj5559X//79w95n1qxZmjlzZo3b27RpE1/jAQAAACSVI0eOyO/313qMo4HWjBkzQgY1p/vwww8lST6fr8bvDMMIeXvAU089pYKCAi1fvlzt2rXT+++/rzvvvFMtWrTQlVdeGfI+06ZN05QpU6p+Pnz4sNq1a6ddu3ZF7EzEp6ysTG3atNHu3buVkZHhdHOSDv1rPfrYWvSv9ehja9G/1qOPrUcfn2IYho4cOaKWLVtGPNbRQGvChAkaOXJkrcdkZ2fr448/1r59+2r87ptvvlFmZmbI+/3www/65S9/qWXLlumaa66RJF100UUqLCzU//t//y9soJWWlqa0tLQat/v9/pR+U9khIyODPrYQ/Ws9+tha9K/16GNr0b/Wo4+tRx8r6skXRwOtZs2aqVmzZhGPy8nJUWlpqTZu3KjevXtLkjZs2KDS0lL17ds35H3Ky8tVXl6uOnWC633UrVtXlZWViTceAAAAAMLwRNXBzp07a+jQoRo7dqwKCgpUUFCgsWPH6tprrw0qhNGpUyctW7ZM0qloe8CAAbrvvvv03nvvaefOnXr55Zf1yiuvaMSIEU79KQAAAABSgCeKYUjSq6++qokTJ2rw4MGSpGHDhumZZ54JOmb79u1VVQIlacmSJZo2bZp++tOf6ttvv1W7du300EMPafz48VE/b1pamh544IGQywlhDvrYWvSv9ehja9G/1qOPrUX/Wo8+th59HDufEU1tQgAAAABA1DyxdBAAAAAAvIRACwAAAABMRqAFAAAAACYj0AIAAAAAkxFo1WLu3Llq37690tPT1bNnT61du9bpJnnWrFmzdOmll+qss85S8+bNdf3112v79u1BxxiGoRkzZqhly5Zq2LChrrjiCn366acOtdjbZs2aJZ/Pp0mTJlXdRv8mbu/evfrZz36mpk2b6owzzlCPHj20adOmqt/Tx/E7efKk/vu//1vt27dXw4YN1aFDBz344INB+x7Sv7F5//33dd1116lly5by+Xx6/fXXg34fTX8eP35cv/jFL9SsWTM1atRIw4YN0549e2z8K9yttj4uLy/X1KlTdeGFF6pRo0Zq2bKlRo8erX/+859Bj0EfhxfpPXy6O+64Qz6fT3PmzAm6nf6tXTR9XFRUpGHDhsnv9+uss85Snz59tGvXrqrf08fhEWiFsXTpUk2aNEm/+tWvtGXLFl1++eW66qqrgt5YiN6aNWt01113qaCgQHl5eTp58qQGDx6so0ePVh3zyCOP6PHHH9czzzyjDz/8UFlZWRo0aJCOHDniYMu958MPP9Rzzz2niy66KOh2+jcxhw4dUr9+/VS/fn397W9/07Zt2/TYY4+pcePGVcfQx/GbPXu25s+fr2eeeUZFRUV65JFH9Oijj+rpp5+uOob+jc3Ro0fVvXv3GluhBETTn5MmTdKyZcu0ZMkSrVu3Tt99952uvfZaVVRU2PVnuFptffz9999r8+bNmj59ujZv3qzXXntNn3/+uYYNGxZ0HH0cXqT3cMDrr7+uDRs2qGXLljV+R//WLlIf79ixQ/3791enTp303nvv6R//+IemT5+u9PT0qmPo41oYCKl3797G+PHjg27r1KmTcf/99zvUouSyf/9+Q5KxZs0awzAMo7Ky0sjKyjIefvjhqmOOHTtm+P1+Y/78+U4103OOHDlinHfeeUZeXp4xYMAA4+677zYMg/41w9SpU43+/fuH/T19nJhrrrnGuPXWW4Nuu+GGG4yf/exnhmHQv4mSZCxbtqzq52j68/Dhw0b9+vWNJUuWVB2zd+9eo06dOsaKFStsa7tXVO/jUDZu3GhIMr7++mvDMOjjWITr3z179hitWrUytm7darRr18544oknqn5H/8YmVB/fdNNNVd/DodDHtWNGK4QTJ05o06ZNVZsjBwwePFjr1693qFXJJbCxdJMmTSRJO3fuVElJSVCfp6WlacCAAfR5DO666y5dc801uvLKK4Nup38Tt3z5cvXq1Uv/8R//oebNm+viiy/WwoULq35PHyemf//++vvf/67PP/9ckvSPf/xD69at09VXXy2J/jVbNP25adMmlZeXBx3TsmVLdevWjT6PU2lpqXw+X9VMOH2cmMrKSuXm5uq+++5T165da/ye/k1MZWWl3nrrLZ1//vkaMmSImjdvrssuuyxoeSF9XDsCrRAOHDigiooKZWZmBt2emZmpkpISh1qVPAzD0JQpU9S/f39169ZNkqr6lT6P35IlS7R582bNmjWrxu/o38R9+eWXmjdvns477zy98847Gj9+vCZOnKhXXnlFEn2cqKlTp2rUqFHq1KmT6tevr4svvliTJk3SqFGjJNG/ZoumP0tKStSgQQOdffbZYY9B9I4dO6b7779fN998szIyMiTRx4maPXu26tWrp4kTJ4b8Pf2bmP379+u7777Tww8/rKFDh2rlypUaMWKEbrjhBq1Zs0YSfRxJPacb4GY+ny/oZ8MwatyG2E2YMEEff/yx1q1bV+N39Hl8du/erbvvvlsrV64MWjddHf0bv8rKSvXq1Uu/+93vJEkXX3yxPv30U82bN0+jR4+uOo4+js/SpUu1aNEi/eEPf1DXrl1VWFioSZMmqWXLlhozZkzVcfSvueLpT/o8duXl5Ro5cqQqKys1d+7ciMfTx5Ft2rRJTz75pDZv3hxzX9G/0QkUIxo+fLgmT54sSerRo4fWr1+v+fPna8CAAWHvSx+fwoxWCM2aNVPdunVrROL79++vMfqH2PziF7/Q8uXLtXr1arVu3brq9qysLEmiz+O0adMm7d+/Xz179lS9evVUr149rVmzRk899ZTq1atX1Yf0b/xatGihLl26BN3WuXPnqgI5vIcTc9999+n+++/XyJEjdeGFFyo3N1eTJ0+umqGlf80VTX9mZWXpxIkTOnToUNhjEFl5ebluvPFG7dy5U3l5eVWzWRJ9nIi1a9dq//79atu2bdV57+uvv9Y999yj7OxsSfRvopo1a6Z69epFPPfRx+ERaIXQoEED9ezZU3l5eUG35+XlqW/fvg61ytsMw9CECRP02muvadWqVWrfvn3Q79u3b6+srKygPj9x4oTWrFlDn0fhX//1X/XJJ5+osLCw6l+vXr3005/+VIWFherQoQP9m6B+/frV2JLg888/V7t27STxHk7U999/rzp1gk9JdevWrRpRpX/NFU1/9uzZU/Xr1w86pri4WFu3bqXPoxQIsr744gu9++67atq0adDv6eP45ebm6uOPPw4677Vs2VL33Xef3nnnHUn0b6IaNGigSy+9tNZzH30cgTM1ONxvyZIlRv369Y0XXnjB2LZtmzFp0iSjUaNGxldffeV00zzp5z//ueH3+4333nvPKC4urvr3/fffVx3z8MMPG36/33jttdeMTz75xBg1apTRokULo6yszMGWe9fpVQcNg/5N1MaNG4169eoZDz30kPHFF18Yr776qnHGGWcYixYtqjqGPo7fmDFjjFatWhl//etfjZ07dxqvvfaa0axZM+O//uu/qo6hf2Nz5MgRY8uWLcaWLVsMScbjjz9ubNmypariXTT9OX78eKN169bGu+++a2zevNkYOHCg0b17d+PkyZNO/VmuUlsfl5eXG8OGDTNat25tFBYWBp37jh8/XvUY9HF4kd7D1VWvOmgY9G8kkfr4tddeM+rXr28899xzxhdffGE8/fTTRt26dY21a9dWPQZ9HB6BVi2effZZo127dkaDBg2MSy65pKoUOWInKeS/l156qeqYyspK44EHHjCysrKMtLQ041/+5V+MTz75xLlGe1z1QIv+Tdybb75pdOvWzUhLSzM6depkPPfcc0G/p4/jV1ZWZtx9991G27ZtjfT0dKNDhw7Gr371q6ALUvo3NqtXrw75vTtmzBjDMKLrzx9++MGYMGGC0aRJE6Nhw4bGtddea+zatcuBv8adauvjnTt3hj33rV69uuox6OPwIr2HqwsVaNG/tYumj1944QXj3HPPNdLT043u3bsbr7/+etBj0Mfh+QzDMKydMwMAAACA1EKOFgAAAACYjEALAAAAAExGoAUAAAAAJiPQAgAAAACTEWgBAAAAgMkItAAAAADAZARaAAAAAGAyAi0AAAAAMBmBFgAgpWRnZ2vOnDlONyNqX331lXw+nwoLC51uCgAgBgRaAABL+Hy+Wv/dcsstjrTrww8/1Lhx44La+frrrzvSFgBA8qrndAMAAMmpuLi46v9Lly7Vr3/9a23fvr3qtoYNGwYdX15ervr161vernPOOcfy5/CCEydOqEGDBk43AwCSFjNaAABLZGVlVf3z+/3y+XxVPx87dkyNGzfWH//4R11xxRVKT0/XokWLdPDgQY0aNUqtW7fWGWecoQsvvFCLFy8OetwrrrhCEydO1H/913+pSZMmysrK0owZM4KOmTFjhtq2bau0tDS1bNlSEydOrPrd6UsHs7OzJUkjRoyQz+er+rm6wPK91157TT/5yU90xhlnqHv37srPzw96zh49egTdb86cOUGPecstt+j666/X7373O2VmZqpx48aaOXOmTp48qfvuu09NmjRR69at9eKLL9Zow2effaa+ffsqPT1dXbt21XvvvRf0+23btunqq6/WmWeeqczMTOXm5urAgQNB/TZhwgRNmTJFzZo106BBg0L+rQAAcxBoAQAcM3XqVE2cOFFFRUUaMmSIjh07pp49e+qvf/2rtm7dqnHjxik3N1cbNmwIut/vf/97NWrUSBs2bNAjjzyiBx98UHl5eZKkP//5z3riiSe0YMECffHFF3r99dd14YUXhnz+Dz/8UJL00ksvqbi4uOrncH71q1/p3nvvVWFhoc4//3yNGjVKJ0+ejOlvXrVqlf75z3/q/fff1+OPP64ZM2bo2muv1dlnn60NGzZo/PjxGj9+vHbv3h10v/vuu0/33HOPtmzZor59+2rYsGE6ePCgpFOzhwMGDFCPHj300UcfacWKFdq3b59uvPHGGv1Wr149ffDBB1qwYEFM7QYAxMgAAMBiL730kuH3+6t+3rlzpyHJmDNnTsT7Xn311cY999xT9fOAAQOM/v37Bx1z6aWXGlOnTjUMwzAee+wx4/zzzzdOnDgR8vHatWtnPPHEE1U/SzKWLVtWaxsC7X3++eerbvv0008NSUZRUZFhGIbxwAMPGN27dw+63xNPPGG0a9eu6ucxY8YY7dq1MyoqKqpuu+CCC4zLL7+86ueTJ08ajRo1MhYvXhz03A8//HDVMeXl5Ubr1q2N2bNnG4ZhGNOnTzcGDx4c9Ny7d+82JBnbt283DONUv/Xo0aPWvxMAYB5mtAAAjunVq1fQzxUVFXrooYd00UUXqWnTpjrzzDO1cuVK7dq1K+i4iy66KOjnFi1aaP/+/ZKk//iP/9APP/ygDh06aOzYsVq2bFnMs07hnP68LVq0kKSq541W165dVafOj6ffzMzMoBm3unXrqmnTpjUeNycnp+r/9erVU69evVRUVCRJ2rRpk1avXq0zzzyz6l+nTp0kSTt27Ki6X/X+BgBYh2IYAADHNGrUKOjnxx57TE888YTmzJmjCy+8UI0aNdKkSZN04sSJoOOqF83w+XyqrKyUJLVp00bbt29XXl6e3n33Xd1555169NFHtWbNmoSLbZx+f5/PJ0lVz1unTh0ZhhF0fHl5ea2PEXic2v6e2pzehuuuu06zZ8+ucUwgIJRq9jcAwDrMaAEAXGPt2rUaPny4fvazn6l79+7q0KGDvvjii5gfp2HDhho2bJieeuopvffee8rPz9cnn3wS8tj69euroqIi0abrnHPOUUlJSVCwZebeVwUFBVX/P3nypDZt2lQ1a3XJJZfo008/VXZ2ts4999ygfwRXAOAMAi0AgGuce+65ysvL0/r161VUVKQ77rhDJSUlMT3Gyy+/rBdeeEFbt27Vl19+qf/5n/9Rw4YN1a5du5DHZ2dn6+9//7tKSkp06NChuNt+xRVX6JtvvtEjjzyiHTt26Nlnn9Xf/va3uB+vumeffVbLli3TZ599prvuukuHDh3SrbfeKkm666679O2332rUqFHauHGjvvzyS61cuVK33nqrKUEkACB2BFoAANeYPn26LrnkEg0ZMkRXXHGFsrKydP3118f0GI0bN9bChQvVr18/XXTRRfr73/+uN998U02bNg15/GOPPaa8vDy1adNGF198cdxt79y5s+bOnatnn31W3bt318aNG3XvvffG/XjVPfzww5o9e7a6d++utWvX6o033lCzZs0kSS1bttQHH3ygiooKDRkyRN26ddPdd98tv98flA8GALCPz6i+oBwAAAAAkBCGuQAAAADAZARaAAAAAGAyAi0AAAAAMBmBFgAAAACYjEALAAAAAExGoAUAAAAAJiPQAgAAAACTEWgBAAAAgMkItAAAAADAZARaAAAAAGAyAi0AAAAAMNn/By7xRo3NRDf6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1000x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig = plt.figure(figsize=(10, 5))\n",
    "ax = plt.subplot(111)\n",
    "ax.set_xlim([0, N])\n",
    "ax.set_ylim([-0.8, 1.0])\n",
    "ax.set_xlabel(\"Transit number\")\n",
    "ax.set_ylabel(\"TTV [hours]\")\n",
    "plt.scatter(\n",
    "    range(N),\n",
    "    (np.array(tt_list) - m * np.array(range(N)) - c) * (24.0 * 365.0 / 2.0 / np.pi),\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5dc5a888-a29c-49e4-8472-5f1ffd6029be",
   "metadata": {},
   "source": [
    "The transit times deviate from a purely linear behaviour because of the perturbations induced onto the inner planet by the outer planet."
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}