diff --git a/jupyter_examples/Lammers24DynamicalModels.ipynb b/jupyter_examples/Lammers24DynamicalModels.ipynb
new file mode 100644
index 0000000..1fac831
--- /dev/null
+++ b/jupyter_examples/Lammers24DynamicalModels.ipynb
@@ -0,0 +1,425 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "9eb21d05",
+   "metadata": {},
+   "source": [
+    "# Simplified dynamical models"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e0baabcd",
+   "metadata": {},
+   "source": [
+    "This notebook illustrates how to set up and integrate dynamical models that include certain resonant interactions between the planets using celmech, as discussed at length in [Lammers et al. 2024](https://ui.adsabs.harvard.edu/abs/2024arXiv240317928L/abstract). First, we specify the initial conditions using a REBOUND simulation, which we use to initialize a PoincareHamiltonian object. Then, we generate the corresponding equations of motion and numerically integrate them, monitoring for close encouters between the planets."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "2882ff03",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "%matplotlib inline\n",
+    "from tqdm import tqdm\n",
+    "import rebound as rb\n",
+    "from celmech import Poincare, PoincareHamiltonian"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0ed0b647",
+   "metadata": {},
+   "source": [
+    "# Simulation setups"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "b2991344",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# get REBOUND simulation of idealized system (coplanar, equal spacings, masses, and eccentricites)\n",
+    "def initialize_sim(P_ratio, P1=0.0316, m=3e-6, e_norm=0.0, Npl=5):\n",
+    "    # calculate eccentricity\n",
+    "    e_cross = (P_ratio**(2/3) - 1)/(1 + P_ratio**(2/3))\n",
+    "    e = e_norm*e_cross\n",
+    "    \n",
+    "    # setup sim object\n",
+    "    sim = rb.Simulation()\n",
+    "    sim.units = ('yr', 'AU', 'Msun')\n",
+    "    sim.add(m=1.00)\n",
+    "    \n",
+    "    # add planets\n",
+    "    for i in range(Npl):\n",
+    "        sim.add(m=m, P=P1*(P_ratio**i), e=e, pomega=np.random.uniform(0.0, 2*np.pi), l=np.random.uniform(0.0, 2*np.pi))\n",
+    "    sim.move_to_com()\n",
+    "    \n",
+    "    return sim"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "93e3ef4e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# get PoincareHamiltonian that includes a specific set of MMRs\n",
+    "def construct_Hamiltonian(sim, res_indices):\n",
+    "    # initialize Hamiltonian\n",
+    "    pvars = Poincare.from_Simulation(sim)\n",
+    "    pham = PoincareHamiltonian(pvars)\n",
+    "    \n",
+    "    # add resonances\n",
+    "    for (j, k, ind_pl1, ind_pl2) in res_indices:\n",
+    "        pham.add_MMR_terms(j, k, max_order=k, indexIn=ind_pl1, indexOut=ind_pl2, inclinations=False)\n",
+    "    \n",
+    "    # set integration settings\n",
+    "    pham._update()\n",
+    "    pham.set_integrator('dop853', nsteps=int(1e9), atol=1e-100, rtol=1e-5)\n",
+    "    \n",
+    "    return pham"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "ac83d33b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# numerically integrate PoincareHamiltonian model to predict the system's instability time\n",
+    "def predict_t_inst(pham, t_max=1e4):\n",
+    "    # determine minimum stopping distance\n",
+    "    ps = pham.state.particles\n",
+    "    stopping_dist = ps[1].a*(ps[1].m**(1/3))\n",
+    "    Npl = len(ps) - 1\n",
+    "    P1 = ps[1].P\n",
+    "    \n",
+    "    # integrate\n",
+    "    times = np.logspace(0.0, np.log10(t_max), 10000)*P1\n",
+    "    for t in times:\n",
+    "        pham.integrate(t)\n",
+    "        \n",
+    "        for i in range(1, Npl):\n",
+    "            if (1 - ps[i+1].e)*ps[i+1].a - (1 + ps[i].e)*ps[i].a < stopping_dist:\n",
+    "                return pham.t/P1\n",
+    "        \n",
+    "    return t_max"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1888e901",
+   "metadata": {},
+   "source": [
+    "# Model (1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "2ccda7a6",
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "100%|████████████████████████████████████████████████████████████████████████████████████| 25/25 [02:37<00:00,  6.31s/it]\n"
+     ]
+    }
+   ],
+   "source": [
+    "P_ratios = []\n",
+    "log_t_insts = []\n",
+    "\n",
+    "# predict instability times for 25 systems\n",
+    "for _ in tqdm(range(25)):\n",
+    "    P_ratio = np.random.uniform(1.03, 1.11)\n",
+    "    sim = initialize_sim(P_ratio)\n",
+    "    \n",
+    "    j = int(np.round(1/(1 - 1/P_ratio))) # get index of nearest first-order MMR\n",
+    "    res_inds = [[j, 1, i, i+1] for i in range(1, 5)] # all resonances to be included\n",
+    "    \n",
+    "    pham = construct_Hamiltonian(sim, res_inds)\n",
+    "    t_inst = predict_t_inst(pham)\n",
+    "    \n",
+    "    P_ratios.append(P_ratio)\n",
+    "    log_t_insts.append(np.log10(t_inst))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "id": "b10d6d29",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle \\left( -0.4, \\  4.5\\right)$"
+      ],
+      "text/plain": [
+       "(-0.4, 4.5)"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAHhCAYAAABqR2T2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABGYUlEQVR4nO3deXwN9/7H8fcJ2ZAgYgtBU1Ttsau6llJVW1tdaRXVSzet7i23Sau91eXebq72tlVFLbcoSje0olVKENTW0lhii92JSkKS+f2hOb+Ec+IkmbNk8no+Hnk8Tma+M/M5Z3DevvOd79gMwzAEAABgYQG+LgAAAMDTCDwAAMDyCDwAAMDyCDwAAMDyCDwAAMDyCDwAAMDyCDwAAMDyCDwAAMDyCDwAAMDyCDwAAMDyCDwASpz4+HjZbDbZbDaPHaNr166y2Wzq2rWrKfsbNmyYbDabHnjgAVP2dzHDMNSsWTPZbDZNmTLFI8cASjICD1BKJCQkOEKCzWZTWFiYzp49e9nt0tPTVbFixXzbJiQkeL5gC1m3bp2mTp2qoKAgPffcc07bnDp1SkuXLtUrr7yim266SVFRUY7P253QZbPZNHbsWEnS888/rzNnzpj5FoASj8ADlFJnzpzRggULLttu4cKFstvtni/IwsaOHSvDMDRs2DDVqVPHaZvY2Fhdf/31GjdunBYuXKhDhw4V+ji33367GjVqpMOHD+u9994rbtmApRB4gFIoJCREkjR9+vTLts1tk7sNCicxMVFLliyRJD3xxBMu2xmG4XhdvXp19e3bt9DHCggI0GOPPSZJ+te//qX09PRC7wOwKgIPUAr1799fkrR06VIdPnzYZbsjR444vqwHDBjgldqs5u2335YktWvXTg0aNHDZ7uGHH9bcuXO1d+9eHT58WIsWLSrS8W677TYFBgbq+PHjmjFjRpH2AVgRgQcoha6//nrVqFFD2dnZmjVrlst2s2bNUlZWlqpXr66ePXt6sUJrOH36tObNmydJuvvuuwts++STT2rgwIEuL3m5KyIiQjfccIMkafLkycXaF2AlBB6gFCpTpozuuusuSQVf1po2bZokadCgQSpTpoxb+z537pwmTZqkbt26qWrVqgoKClKNGjV044036rPPPlNOTs5l97F//3499NBDiomJUUhIiKKiotS/f38tW7bMrRpynT17Vm+//ba6deum6tWrKygoSNWqVdP111+vKVOmKDs7u1D7K6yFCxcqMzNTkjRw4ECPHiuv3GP98ssv2rt3r9eOC/gzAg9QSt1zzz2SpKSkJG3duvWS9du2bdOGDRvytb2cvXv3qmXLlnrooYeUkJCgY8eO6fz580pNTdU333yje+65R126dNGJEydc7mPFihVq0qSJJk2apN27dyszM1OHDh3SokWL1LNnT7300ktu1ZKYmKiGDRtqzJgxSkhI0JEjR3T+/HkdPXpUS5cu1fDhw3XNNdcoNTXVrf0VxfLlyyVJ0dHRioqK8thxLtahQwfH6xUrVnjtuIA/I/AApVRsbKyaNm0qyXkvT+6yJk2aKDY29rL7O3PmjLp3767t27dLkm666SZ9+eWXWrdunebMmaMuXbpIklauXKm+ffs67V3Zs2eP+vXrJ7vdroCAAI0aNUrLli1TYmKiJk+erAYNGiguLk5fffVVgbX8+uuv6tatmw4cOKBq1aopLi5Oy5YtU1JSkr777js99NBDKlu2rNauXasBAwbo/Pnzl31/RfHTTz9Jktq2beuR/bvSsGFDVapUKV8NQKlnACgVli9fbkgyJBlTpkwxDMMwXnvtNUOSUbt2bSM7O9vRNicnx4iOjjYkGRMmTDAMwzCmTJni2H758uWX7P/JJ590rB83btwl63NycozBgwc72kyaNOmSNrfccotj/cyZMy9Zb7fbjRYtWjjaOPsnLCcnx2jevLkhyWjRooVx9OhRp5/HN998YwQEBBiSjI8//viS9V26dDEkGV26dHG6/eWkpqY6ahw/fnyR9pG7fVFq6NatmyHJaNy4cZGODVgNPTxAKTZ48GAFBARo//79+S59JCQkKCUlRQEBARo8ePBl95OZmamPP/5YktS4cWPFx8df0sZms2nSpEmqUqWKJGnixIn51h86dEgLFy6UJPXt29cxxiivsLAwffjhhwXW8tVXX2nz5s2SLoxBioyMdNruhhtu0K233ipJHpmZeP/+/Y7X1apVM33/l5N7zLx1AKUZgQcoxWrVqqVu3bpJyn9ZK/d1165dVbt27cvuZ/369Tp16pQkaejQoS4HOIeHh+v222+XdGGMUN7J9ZYvX+64zDVs2DCXx2rXrp2aNGnicn1uaLrqqqvUvHnzAuv+29/+JunCeB+zBzAfPXrU8bpy5cqm7tsdERERkiS73a5z5855/fiAvyHwAKXckCFDJElz585Venq60tPTHbdSuztYecuWLY7X7du3L7Bt3vV5t/v1118dry835qVdu3Yu161bt06S9Ntvv+V7HIazn4cffljShTvLChpIXRR59+eLwJP3mMePH/f68QF/Q+ABSrlbbrlF5cqVU1pamhYuXKgFCxbIbrcrNDTU7Vup8365V69evcC2NWrUcLrdyZMnHa8vdwmooGMcOXKkwG1dcee5YoWRd2ZqX8x4nPeYoaGhXj8+4G/K+roAAL5VoUIF3XzzzZoxY4amT5/ueMTBTTfdpLCwsELv73JPMDfyPELB1fKi7kOS49JUp06d9MEHHxS4n7zMvm28atWqjtdm9x65I/eYgYGBjju2gNKMwANAQ4YM0YwZMxyPkZDcv5wl/f94EUk6fPiwGjZs6LJt3nlv8m6X93Vqaqqio6Nd7qOgXpwqVaooNTVVR48eddx27wt5A0/e3itvyT2mq0HbQGnDJS0Auu6661SzZk1lZWU5HiVx/fXXu7193mCxZs2aAtuuXbvW6XbNmjVzvE5MTCxwHwWtz50z6Pfff/fpLMNXXnmlypUr56jF23KPmfdzBUozAg8AlSlTRvfcc4+Cg4MVHBysu+++2+1HSUhS69atHZdNpk6d6vKOp7S0NH3++eeSLty+XrNmTce6bt26OY45depUl8dat25dvsHOF8t9MKokvf76626/B7OVLVvWMePx5QKc2ex2u3777TdJUufOnb16bMBfEXgASJJee+01ZWRkKCMjQ2+++Wahtg0ODtaIESMkSVu3btWLL754SRvDMPTwww/r2LFjkuS4QypXzZo1HU9k//LLLx3BKK8zZ87o73//e4G1DBw4UFdffbUk6f3337/sAzS3bNlS5CeTX05u2Ni8ebPjmVresG7dOsc4JwIPcAFjeACY4oUXXtAXX3yh5ORkjR8/Xlu2bNHw4cMVFRWl3bt3a+LEiUpISJAkdezY0Wlw+de//qWlS5cqLS1NgwYN0ooVK3TrrbcqPDxcmzdv1oQJE/T777+rTZs2jtvPL1amTBn973//0zXXXKMzZ85oxIgRmjNnjgYNGqSrrrpKgYGBOnLkiJKSkrR48WKtWrVKTzzxhPr162f6Z9KnTx+9+OKLOnfunH766Sf16NHDZduNGzdq48aNTtcdPnxYn376ab5lt956qypUqOC0/ffffy/pwq3pHTt2LFLtgOX4cppnAN7j7NEShXG5R0sYhmHs3r3baNSoUb5HP1z806lTJ+P48eMF1hkWFuZy+7i4OCMuLs7loyVybdq0yWjQoEGBteT+vPjii5dsX9xHS+TK/TyGDRtWYLu878mdn927d7vc1xVXXGFIMkaOHFms2gEr4ZIWANPUq1dPmzZt0sSJE9WlSxdVqVJFgYGBql69um644QZNnz5dP/74Y747si7WtWtXbd26VQ888IDq1q2roKAgVa9eXX369NG3337r9LEVzjRv3lzbtm3T1KlTddNNNyk6OlohISEKCgpSzZo11bVrV40bN07r16/XCy+8YNIncKkHH3xQkjRv3jxlZGR47Di5Vq9erd27d+c7NgDJZhgFTGgBACiWM2fOqF69ejp+/LimT5+uu+++26PHGzFihCZPnqyePXvmm2YAKO3o4QEAD6pQoYKefvppSdIrr7yinJwcjx1r3759mjZtmiQ5HTgOlGYEHgDwsNGjR6tu3brasWOH07vPzPLqq6/q/PnzGjhwIIOVgYtwlxYAeFhISIg+++wzLVu2TFlZWR45hmEYqlu3ruLi4jR8+HCPHAMoyRjDAwAALI9LWgAAwPK4pPWXnJwcHTx4UGFhYZd9UjMAAPAPhmEoLS1NUVFRCghw3Y9D4PnLwYMHC3w6MwAA8F8pKSmqXbu2y/UEnr+EhYVJuvCBhYeH+7gaAADgDrvdrujoaMf3uCsEnr/kXsYKDw8n8AAAUMJcbjgKg5YBAIDlEXgAAIDlEXgAAIDlEXgAAIDlEXgAAIDlEXgAAIDlEXgAAIDlEXgAAIDlEXgAAIDlEXgAAIDlEXgAAIDlEXgAAIDlEXgAAIDlEXgAAIDlEXgAAIDlEXgAAIDlEXgAAIDlEXgAAIDlEXgAAIDlEXgAAIDlEXgAAIDlEXgAAIDlEXgAAIDlEXgAAIDlEXgAAIDlEXgAAIDlEXgAAIDlWTrwjBs3TjabTTabTS+//LKvywEAAD5i2cCzfft2vfHGG74uAwAA+AFLBh7DMDRy5EgFBgaqe/fuvi4HAAD4mCUDz+TJk/XTTz/phRdeUHR0tK/LAQAAPma5wHP06FE988wzaty4scaMGePrcgAAgB8o6+sCzDZmzBidOHFCX3zxhQIDA31dDlAq2O12paenKzQ0VOHh4b4up8Qz+/MsaH9FOZa/nu+8dUnyyxrNUthzUJj2rtr663l3l6UCz/fff68ZM2bo7rvvVpcuXXxdDlAqJCcnKyUlxfF7dHS0YmJifFhRyWb251nQ/opyLH893xfXlZe/1GiWwp6DwrR31dZfz3thWOaSVkZGhkaNGqWKFSvqzTffLPJ+7Ha7Wz+ZmZkmVg+UTHa7/ZIvmZSUFNntdh9VVLKZ/XkWtL+iHMtfz7ezuvLyhxrNUthzUJj2rtoeOnTIL897YVkm8Lz88svatWuXXnnlFVWvXr3I+4mOjlbFihUv+/Pqq6+aWD1QMqWnpxdqOQpm9udZ0P6Kcix/Pd/uHN/XNZqlsOegMMtdtU1LSyvUvv2VJS5p5c6506pVKz3wwAPF2ldKSopb1yaDg4OLdRzACnLHSri7HAUz+/Msyv6Kss7X59ud4/u6RrMU9hwUZrmrtmFhYTp06JDb7f2VJXp4HnzwQWVlZen9999XQEDx3lJ4eLhbPwQe4MLfl4unfoiOji6RAxr9gdmfZ0H7K8qx/PV8O6srL3+o0SyFPQeFae+qbc2aNf3yvBeWzTAMw9dFFFelSpV05swZRUZGXrLu9OnTysjIUIUKFVS+fHlFR0crMTHxknZ2u10VK1bU6dOnS9xJBHytpN+94W+4S6touEvLnPYl7S4td7+/LRN4Tp8+7VbbunXras+ePZcsJ/AAnuOv/1ACKPnc/f62xCWtU6dOyTAMpz/33nuvJGn8+PEyDMNp2AHgOcnJyUpKStKOHTuUlJSk5ORkX5cEoBSyROAB4J/89TZmAKUPgQeAx/jrbcwASh8CDwCP8dfbmAGUPpYYtGwGBi0DnmGFKekB+C93v78tMfEgAP8VExOjyMhI7tIC4FMEHgAelzvRHQD4CmN4AACA5RF4AACA5RF4AACA5RF4AACA5RF4AACA5RF4AACA5RF4AACA5RF4AACA5RF4AACA5RF4AACA5RF4AACA5fEsLQCy2+083BOApRF4gFIuOTlZKSkpjt+jo6MVExPjw4oAwHxc0gJKMbvdni/sSFJKSorsdruPKgIAzyDwAKVYenp6oZbDuux2u1JTUwm7sCwuaQGlWGhoaKGWw5q4rInSgB4eoBQLDw9XdHR0vmXR0dEMXC5FuKyJ0oIeHqCUi4mJUWRkJHdplVIFXdbkzwKshMADQOHh4Xy5lVJc1kRpwSUtACjFuKyJ0oIeHgAo5bisidKAwAMA4LImLI9LWgAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPLKmr3DY8eO6cCBAzp69KiOHz+u0NBQVa1aVVWrVlVMTIwCAshYAADAu4odeNLS0rRw4UIlJCTop59+0q5du1y2LV++vDp06KDOnTurT58+atWqVXEPDwAAcFk2wzCMomy4fv16vfPOO5o3b54yMjIkSe7symazOV5fddVVeuihhzR06FCVL1++KGWYxm63q2LFijp9+rTCw8N9WgsAAHCPu9/fhQ4869ev17hx47RkyRJJ/x9yatasqbZt26p169aqVq2aIiIiVLlyZaWnp+vEiRM6efKkfv/9dyUmJmrz5s06f/78hQJsNkVEROipp57So48+quDg4KK+52Ih8AAAUPJ4JPAMGzZM06dPV05OjiSpVatWGjx4sAYOHKg6deq4Xdy5c+f0448/aubMmZo/f75Onz4tm82mOnXqaPr06br22mvd3pdZCDwAAJQ87n5/F2oE8dSpU1W2bFndf//92rFjh9atW6cxY8YUKuxIUlBQkHr06KFPPvlEhw8f1rRp03TVVVdp7969+uGHHwq1L8Db7Ha7UlNTZbfbfV0KAMBNhRq0/MADD+i5555T7dq1TSsgODhYd999twYPHqw5c+YoOzvbtH0DZktOTlZKSorj9+joaMXExPiwIgCAO4o8aNlquKSFy7Hb7UpKSrpkeWxsLH9mAMBHPHJJy58tWLBAI0eOVOvWrVWzZk0FBQWpUqVKuuaaa/TOO+/o3Llzvi4RJVx6enqhlgMA/IfpEw/6yptvvqmff/5ZwcHBioqKUosWLXTo0CGtXr1aq1ev1vTp07Vs2TJVqlTJ16WihAoNDS3UcgCA/7BMD8+IESO0fPlypaWlKTk5WYmJidq/f79Wr16t2rVra/369Ro7dqyvy0QJFh4erujo6HzLoqOjuZwFACWAX43hWbhwoTZt2qQXXnjB1P3OmTNHt99+u6KionTgwAGnbRjDA3fZ7Xalp6crNDSUPysA4GMlcgzPggUL9OKLL5q+30aNGkmSzp49a/q+UfqEh4erevXqhB0AKEH8KvB4yurVqyWJZ3cBAFBKeXTQ8rRp0wrVvqAHjxZWdna2Dh06pC+//FLPPvusypcvr1dffdW0/QMAgJLDo4Fn6NCh+R4WejmGYRSqvTNvv/22xowZk2/ZTTfdpPHjx6tp06aX3d7d2XODg4N99twvAABQOB4NPEFBQYqKitLIkSPdaj9nzhynE7sVRq1atdSpUyedP39ee/fuVWpqqpYvX65Zs2bppZdeUpkyZQrc/uK7cFyJi4tTfHx8sWoFAADe4dG7tNq2bat9+/YpNTXVrfbDhg3TtGnTTH28xJo1azRy5Eht2rRJo0aN0vvvv++0Xe4o75SUFLcGo9LDAwCA7/nFXVqtW7fWsWPH8j17yNvat2+vr7/+WsHBwfrwww+1d+/eAtuHh4e79UPYAQCg5PBo4OncubPq1KmjnTt3utX+2muv1ZAhQ0yvIyoqSi1btlROTo42bdpk+v4BAIB/86uJBz2pTZs2Wr9+vebNm6dbbrnlkvVMPAgAQMnj7vd3kQct//HHH5o5c6Z27Nih7Oxs1axZU3/729904403+t3lnj179jh6dlq0aOHjagAAgLcVKfD897//1ejRo5WVlZVv+bvvvqvo6Gh98MEHuuGGG0wp0B3r16/Xl19+qXvvvVcxMTH51n377bcaM2aMsrKydOONN+rKK6/0Wl0AAMA/FPqSVmJiojp27KicnByXbcqWLau5c+eqf//+xS7QHQkJCerWrZskqUaNGqpdu7bOnTunffv26dSpU5Iu3DH29ddfKzIy0uk+uKQFAEDJ47G7tN577z3l5OTIZrOpb9++mj9/vtasWaMvvvhCw4cPV2BgoLKysjR8+HAdO3asWG/CXS1atNA777yj/v37q3z58tqxY4d27Nih0NBQ9e7dW1OmTNGqVatchh0AAGBthe7hiYmJ0d69e9W7d28tXrz4kvU//vijevXqpXPnzmn8+PF6/vnnTSvWk+jhAQCg5PFYD8+hQ4ckyeXsyX/729/0xBNPyDAMzZs3r7C7BwAAMF2hA09mZqYkqV69ei7b3HXXXZKkX3/9VefOnStaZQAAACYp8sSDBT2Tqn79+pIuPLH86NGjRT0EAACAKTwy03LeeXjS0tI8cQgAAAC3FTnw2Gw2t9oVdPs6AACANxR5puVrr71WLVq0UMuWLR0/jRs3VtmyRd4lAACARxQpnRiGoZMnT2rFihVasWKFY3lgYKAaN26sli1bOpadP3++2EUCAAAUR6Hn4fnggw+0ceNGbdy4UVu2bNHZs2fz7/CiS10BAQG68sor1bx5czVv3lwtWrRQ8+bNVbdu3eJXbyLm4QEAoORx9/u7WE9LNwxDv/32myMAbdy4UZs2bVJqauqlB7ooCIWFhalZs2Zq0aKFJk6cWNQSTEPgAQCg5PFK4HElNTVVSUlJ+YLQrl27nA5gttlsys7ONruEQiPwAABQ8vg08DiTnp6uTZs25QtBW7ZsUXp6OoEHAAAUibvf3167pSo0NFQdOnRQhw4dHMsMw9Dvv//urRIAAEApVeh5eNavX2/awW02m6666irT9gcAAOBMoQNP27ZtVbt2bY0cOVKLFy9WRkaGJ+oC4Aa73a7U1FTZ7XZflwIAfq3QY3gCAi5kpNy7rkJCQtS9e3f169dPffv2VVRUlPlVegFjeFDSJCcnKyUlxfF7dHS0YmJifFgRAHifu9/fhe7h2b9/vz744APdeOONCgkJUXp6ur766is98MADio6OVuvWrRUfH2/qpS8A+dnt9nxhR5JSUlLo6QEAF4p1l1Z6erqWLVumxYsX66uvvtLBgwcv7PSv3p8aNWqoT58+6tevn3r06KHQ0FBzqvYAenhQkqSmpmrHjh2XLG/UqJGqV6/ug4oAwDd8clv6+vXrtWjRIi1evFgbNmy4cIAScumLwIOSxG63Kykp6ZLlsbGx/PkFUKr4fB6egwcPavHixVq0aJF++OEHpaenXzjgXwGoZcuW6tevn/r166fWrVt7ooRCIfCgpGEMDwD4QeDJKyMjQ8uWLdOiRYtcXvrq27evHnzwQbVo0cLT5ThF4EFJZLfblZ6ertDQUP7cAiiV/CrwXGz9+vWO3p+kpCQZhiGbzaa4uDi98MIL3i5HEoEHAICSyO9mWs6rdevWat26teLi4hyXvhYvXqxy5cr5ohwAAGBxXuvh2bdvnypVqlRg+rLb7Tp16pTq1KnjjZIuOTY9PAAAlCwem4enqK644gq98847BbZ5//33dcUVV3ipIgAAUFp4LfAYhqHLdSb5YDgRAAAoBbwWeNyxf/9+hYWF+boMAABgMR4dtPzSSy/l+z0hIcFpu+zsbO3fv1+zZ89W+/btPVkSAAAohTw6aDn3QaPShTl3LneoqKgozZ8/X23btvVUSS4xaBkAgJLHL25LX758uaQLY3O6d++uoUOH6t57772kXZkyZRQREaFGjRrlC0kAAABm8Gjg6dKli+N1XFycunbtmm8ZAACAN/hkpuW8DMPQrl27FBoaqtq1a/usDi5pAQBQ8vjdPDxffvmlhg8frpMnTzqW7dmzR82aNVOjRo1Ut25dDR48WDk5Od4qCQAAlBJeCzzvv/++EhMTVblyZceyxx57TNu2bVO3bt3UvHlzzZ49W1OmTPFWSQAAoJTwWuDZunWr2rVr5/j99OnT+vrrr3XHHXdo2bJlWrt2ra6++mpNnjzZWyUBAIBSwmuB5+jRo6pZs6bj95UrVyorK0t33XWXJCkwMFA9e/bUrl27vFUSAAAoJbwWeMLDw3X8+HHH7wkJCQoICFDnzp0dywIDA/Xnn396qyQAAFBKeC3wNGrUSIsWLdKJEyd0+vRpzZ49W61atco3pmfv3r2qXr26t0oCAAClhNcCz+jRo3Xw4EHVqlVL0dHROnjwoEaNGuVYn52drZUrV6pFixbeKgkAAJQSHp14MK+BAwfqP//5j2NQ8u23367hw4c71n///fc6e/asbrjhBm+VBAAASgmfTzzoL5h4EACAksfvJh4EAADwFQIPAACwPK+N4ZGkc+fOacGCBUpMTNSpU6eUnZ19SRubzcbkgwAAwFReCzx79+5Vz5499ccff6igYUMEHgAAYDavBZ4xY8Zo165duueeezR8+HDVrl1bZct6tYMJAACUUl5LHD/88IOuu+46TZ061VuHBAAAkOTFQcs5OTmKjY311uEAAAAcvBZ4OnbsqO3bt3vrcAAAAA5eCzwTJkzQ8uXLNXfuXG8dEgAAQJIXx/AsWrRI3bp10x133KEuXbooNjZWFStWvKSdzWbTP/7xD2+VBQAASgGvPVoiIMC9ziSbzeZ0fh5P49ESAACUPO5+f3uth2f58uXeOhQAAEA+Xgs8Xbp08dahAAAA8uFZWgAAwPI81sMzfPhw2Ww2/fOf/1T16tU1fPhwt7bj0RIAAMBsHhu0HBAQIJvNpu3bt6thw4YMWgYAAKbz+aDl3bt3S5Jq1aqV73cAAABvMyXwDBo0SC1atNADDzzgSFd169bN1+bi3wEAALzFlEHLs2fP1vPPP6/9+/ebsbtCMwxDK1eu1FNPPaUOHTqoUqVKCgoKUlRUlAYOHMgt8QAAlHJeuy39YmfOnNGuXbvUtGlTlS1bvDJ++OEH9ejRQ9KFsUP169dX+fLltXPnTn3xxRf64osvNG7cOI0fP96M0uFH7Ha70tPTFRoaKkmO14zDAgDk5bPAs3fvXrVq1UrBwcFKT08v1r4Mw1D9+vX1+OOP684771TlypUlSefOnVN8fLxeffVVvfzyy2rfvr369u1rRvnwA8nJyUpJSXG6Ljo6WjExMV6uCADgr4p0Sev111/X8uXLZbfbi13AuXPnir2Pdu3aafv27XrggQccYUeSgoKC9M9//lO9e/eWJH300UfFPhb8g91udxl2JCklJcWUP58AAGsoUuB59tln1aNHD0VERKhRo0aO5cuXL1dycrJb+8jIyJB04Tb04goPDy/wsljPnj0lSb///nuxjwX/4E6vYHF7DgEA1lGkS1ply5ZVVlaWDMPIFyJGjx4tSQoLC1Pz5s3VsmVLtWzZUi1atFCzZs0UFBTkaPvLL79IkipVqlSM8t2TG65yx3mg5HPnXHK+AQC5ihR40tLStGHDBq1du1aJiYmaOXOmpAtjaaQLlxtWrlypn3/+2bFNmTJldNVVV6lJkybKysrS4sWLZbPZ1KZNGxPehmuGYWjOnDmSpE6dOnn0WPCe8PBwRUdHFziGh4HLAIBcpsy0nDur8pw5c5Senq6NGzdq48aN2rRpk44dO5b/gH9dwjIMQwEBAfrmm28cl5w84cMPP9TIkSMVFBSkbdu26corr3TaLnemxpSUFLe+KIODgxUcHGx2uSgk7tICgNLN3ZmWTQk8VatW1YkTJ/Trr7+qcePG+dYdPHgwXwD67bffZLfbFRMTozFjxqhPnz7FPbxLGzZsUKdOnZSRkaHXX39dTz31lMu2uR+Yu+Li4hQfH29ClQAAoKi8Gngk6dChQ6pcubJCQkLM2F2x7d69W506ddKhQ4c0aNAgffbZZwUOkKaHp2TL29ND7w4AlB5ef5ZWzZo1zdpVsR0+fFg9e/bUoUOH1KdPH3366adu3w0WHh7OF2YJc/F8PMzBAwC4mCmPlvAnJ06cUM+ePfXHH3+oS5cumjNnjgIDA31dFjzE2Xw8zMFjPrvdrtTUVD5XACWWz2Za9oQzZ87oxhtv1JYtW9S2bVstWrSIW5MtztVcO+np6fTUmYQeNABWUKgentzbuz1l//79WrVqVZG2zczM1IABA7RmzRo1adJE3377rcLCwkyuEP7GVaAl6JqDHjQAVlGowHPHHXeoWbNmpgefffv26YEHHlD9+vW1bNmyQm+fnZ2tO++8Uz/88IOuvPJKLV26VBEREabWCP+UOx9PXszBY56CetAAoCQp1CWtBg0aaOvWrbrzzjv11FNPadCgQRo8eLCaNGlS6AP/+eefmj9/vmbMmKHvv/9e2dnZjiedF9bnn3+uBQsWSLowJ9Btt93mtF3NmjU93ksF74uJiVFkZCR3aXkAPWgArKJQt6VnZWXp3Xff1YQJE3Ts2DHHnU8NGjRQhw4d1LZtW8XGxqpatWqqXLmyKleurPT0dJ04cUInT57U77//rsTERK1du1Zr165VRkaGY3bm3r1767XXXlPTpk0L/SY+/fRTDRs27LLt6tatqz179jhd5+5tbUBpwxgeAP7Mo/Pw/Pnnn5o0aZL+85//aN++fRd2VIiHgOYeskyZMhowYICeeuoptW/fvrBlmIrAA7jGPEcA/JVXJh7MycnR0qVL9fnnn2v58uUue0/yCg0NVbt27dSnTx8NGjRIUVFRRT28qQg8AACUPF6faVmSDhw4oFWrVmn//v06evSoTpw4oZCQEFWtWlVVq1ZVs2bN1KZNG7+cF4fAAwBAyeP1mZYlqVatWi4HDAMAAPiK5WZaBgAAuBiBBwAAWB6BBwAAWB6BBwAAWB6BBwAAWB6BBwAAWB6BBwAAWB6BBwAAWJ6pEw+68scff2jKlClasWKFdu7cqdOnT0uSKlasqAYNGqhr16669957i/SkdAAAgMsx9dESzkyYMEFxcXE6f/68JCkyMtIx9bPdbtexY8ckSYGBgXrxxRf17LPPerIcl3i0BAAAJY+7398evaQ1a9YsPf/882rYsKFmz56tEydO6MiRI9q1a5d27dqlI0eO6MSJE5o1a5YaNGigsWPHavbs2Z4sCQAAlEIe7eFp3769jh8/rk2bNql8+fIFtk1LS1PLli0VGRmpNWvWeKokl+jhAQCg5PGLHp6tW7fqlltuuWzYkaSwsDDdcsst2rp1qydLAgAApZBHA09QUJBjgLI77Ha7goKCPFgRAAAojTwaeDp06KDZs2dr8+bNl227adMmzZo1Sx07dvRkSQAAoBTy6BiexMREde7cWQEBARo8eLB69uypBg0aqGLFipKk06dPa+fOnVqyZIlmzpypnJwcrVy5Um3atPFUSS4xhgcAgJLH3e9vj9+WnpCQoBEjRig5OVk2m81pG8MwFBMTo48//lhdu3b1ZDkuEXgAACh53P3+9vjEg127dtVvv/2mH374QQkJCU4nHuzSpYuuu+46lSlTxtPlAACAUsjjPTwlBT08AACUPH5xWzoAAIA/8KvAs3DhQr300ku+LgMAAFiMXwWeBQsW6MUXX/R1GQAAwGL8KvAAAAB4gkfv0po2bVqh2u/atctDlQAAgNLMo4Fn6NChLufeccYwjEK1BwAAcIdHA09QUJCioqI0cuRIt9rPmTNHSUlJniwJAACUQh4NPM2aNdO+ffv0zDPPuNV+x44dBB4AAGA6jw5abt26tY4dO6aUlBRPHgYAAKBAHu3h6dy5s7777jvt3LlT0dHRl21/7bXXerIcAABQSvFoib/waAkAAEoeHi0BAADwFwIPAACwPAIPAACwPFMGLe/bt6/Q29hsNoWEhKhixYoKCgoyowwAAACnTAk8V1xxRbG2r127tjp06KChQ4eqd+/eZpQEAADgYMolLcMwivWTkpKiuXPnqm/fvurVq5dOnjxpRlkAAACSTOrhmTJliiTpgw8+0Jo1axQSEqJevXqpTZs2qlq1qiTp6NGjWrdunb777jtlZmaqXbt2+vvf/y673a4tW7Zo4cKFOnbsmJYtW6abbrpJK1asMKM0AAAA8+bhGTVqlD766CP1799f//3vf1WtWjWn7Y4cOaK///3vWrRokYYNG6aPP/5YkpSRkaFRo0Zp2rRpstlsmjlzpu644w4zSnML8/AAAFDyeHUengULFujDDz9Ux44d9cUXX7gMO5JUrVo1zZ8/Xx06dNCUKVP0+eefS5JCQkL0ySefqFWrVpKk2bNnm1EaAACAOYFn0qRJstlsevTRR2Wz2S7b3maz6bHHHpNhGPrwww//v5iAAN1///0yDEPr1q0zozQAAABzAs/mzZslSfXr13d7m9y2v/76a77lzZs3lyQdP37cjNIAAADMCTx2u13ShYHJ7sptm5aWlm95aGioJCkwMNCM0gAAAMwJPLlPQp8xY4bb20yfPl2SVKdOnXzLjxw5IkmOu7sAAACKy5TA079/fxmGoc8++0xvv/32Zdu/9dZbmjFjhmw2m/r3759v3Zo1ayRJdevWNaM0AAAAc25LP378uK6++mrHuJtWrVppyJAhat26teOOrSNHjmjdunWaPn26NmzYIMMwVK1aNW3btk0RERGOfTVp0kQ7duzQK6+8omeffba4pbmN29LdZ7fblZ6ertDQUD4rAIBPufv9bcrEg1WqVNHSpUt1ww03KDU1VRs2bNCGDRtctjcMQzVq1NC3336bL+wkJyerXbt2ateunQYOHGhGaTBZcnKyUlJSHL9HR0crJibGhxUBAHB5pk08KEmnTp1SXFycpk2bptOnTzttU7FiRQ0ZMkTx8fGqXLmyWYcuNqv18HiiF8ZutyspKemS5bGxsZb4zAAAJY+739+mBp5cmZmZWr9+vbZs2eJ4LlblypXVpEkTtWnTRsHBwWYfstisFHg81QuTmpqqHTt2XLK8UaNGql69erH3DwBAYXn1ktbFgoODdc011+iaa67xxO5RALvdni/sSFJKSooiIyOLHeRypwxwdzkAAP7ClLu04D/S09MLtbwwwsPDHVMQ5IqOji7xPWIAAOvzSA+PdOHyx5YtW3TixAlJUkREhJo2bcqlDw/zdC9MTEyMIiMjLxkfxJ1bAAB/ZmrgyX021sSJE7Vt2zanbRo3bqxHHnlE999/v1vP3ULh5PbCXDyGx8wQEh4enm9/3LkFAPB3pg1aPnnypPr166fVq1dLuhB+nB7wr5BzzTXXaNGiRapUqZIZhy82Kw1alrzX48KdWwAAX/LqoGXDMDRgwACtWrVK0oV5eW6//Xa1b99eNWrUkGEYSk1N1dq1a/X555/r2LFjWrVqlQYMGKAVK1aYUQIucnEvjKcUNGaIwAMA8Bem9PDMmDFD99xzj2w2mwYNGqRJkyYpLCzMadszZ87ooYce0vTp02Wz2fTZZ5/prrvuKm4J2r17t5YtW6a1a9dq7dq12rp1q7KzszV+/HiNGzfusttbrYfHW1z18DRs2FABAQGM6QEAeJRXe3hmzpwpSerSpYvjoaCuVKhQQVOnTtW+ffu0YsUK0wLPO++8o3feeafY+0HhOBszFBYWpt9//93xO2N6AAC+Zspt6Rs2bJDNZtPDDz/s9jaPPPKIJDntHSiKyMhI9e3bVy+99JK++eYbHk3hRTExMYqNjVWjRo3UsGFDpaWl5VufkpIiu93uo+oAADCphyf31vMrrrjC7W1y2+ZuW1wXX7aaPXu2KfuFe3LHDKWmpjpdz5geAIAvmdLDU7FiRUnSwYMH3d4mty1fgtbCbMwAAH9kSuBp2rSpJGnKlClub/PJJ5/k2xbWwGzMAAB/ZErgufXWW2UYhubPn6/4+HiXc/BIF25hj4+P1/z582Wz2XTbbbeZUQL8SN4xPbGxsQxYBgD4nCljeO6//35NnDhRO3bs0Pjx4zVv3jwNHTpU7du3V/Xq1WWz2XT48GGtWbNGU6dO1datWyVdeMr2/fffb0YJpnF3cG1wcLBfPvXdX3hrHiAAANxhSuAJDAzUN998o+7du2v37t3atm2bnn76aZftDcNQTEyMvvnmG5Ut67HHeRXJxZdjXImLi1N8fLxniwEAAKYwLW3UrVtXmzdvVnx8vCZPnqxTp045bVepUiWNGDFCL7zwgipUqGDW4U2TkpLiVs8EvTsAAJQcpnavlC9fXm+88YZeeeUVrV+/3unT0lu3bq2goCAzD2sqLsUAAGA9HrmeFBQUpI4dO6pjx46e2D0AAEChmHKXFgAAgD8rVA/Pvn37PFJEnTp1PLJfAAAAqZCBpzCPjnCXzWZTVlaW6fsFAADIZTMKmiXwIgEB5l8Bs9lsys7OLvZ+fv75Zw0YMMDx+5kzZ5SZmaly5crle6xBUlKS01vP3X28PAAA8B/ufn8XqoenMI+O8Lbz58/r+PHjlyw/e/aszp496/jdjHAFAABKlkL18FgZPTwAAJQ87n5/c5cWAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwPAIPAACwvLK+LsBK7Ha70tPTFRoaqvDwcF+XAwAA/kLgMUlycrJSUlIcv0dHRysmJsaHFQEAgFxc0jKB3W7PF3YkKSUlRXa73UcVAQCAvAg8JkhPTy/UcgAA4F0EHhOEhoYWajkAAPAuAo8JwsPDFR0dnW9ZdHQ0A5cBAPATDFo2SUxMjCIjI7lLCwAAP0TgMVF4eDhBBwAAP8QlLQAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHkEHgAAYHmWCzxff/21evTooYiICJUvX16tWrXSe++9p5ycHF+XBgAAfMRSgWfChAnq06ePvv/+e1WuXFn169fXpk2bNHr0aN18882EHgAASinLBJ7Vq1fr+eefV0BAgGbOnKk//vhDmzZt0oYNG1S9enV9+eWX+ve//+3rMgEAgA9YJvC8/PLLMgxDI0aM0F133eVY3qJFC0fQmTBhgs6fP++rEgEAgI9YIvDY7XYtW7ZMknTfffddsv62225TeHi4jh8/ruXLl3u7PAAA4GOWCDxJSUk6d+6cQkJC1KpVq0vWBwYGqm3btpKkNWvWeLs8AADgY2V9XYAZdu7cKUmqU6eOypZ1/pZiYmL0/fffO9q6Yrfb3TpmcHCwgoODC1coAADwCUsEnpMnT0qSKleu7LJN7rrctq5ER0e7dcy4uDjFx8e7VyAAAPApSwSejIwMSVJQUJDLNrm9Menp6QXuKyUlReHh4Zc9Jr07AACUHJYIPCEhIZKkc+fOuWyTmZkpSQoNDS1wX+Hh4W4FHgAAUHJYYtCyO5er3LnsBQAArMkSgadBgwaSpH379ikrK8tpm+Tk5HxtAQBA6WGJwBMbG6vAwEBlZGRow4YNl6w/f/68EhMTJUnt27f3dnkAAMDHLBF4wsPD1aNHD0nS5MmTL1k/Z84c2e12ValSRV27dvVydQAAwNcsEXgkaezYsbLZbPr44481a9Ysx/JNmzbp8ccflyQ9/fTTBd7JBQAArMlmGIbh6yLM8sorr2jcuHGSLkw0WKFCBW3ZskU5OTnq06ePFi5cqDJlyjjd1m63q2LFijp9+jR3aQEAUEK4+/1tmR4e6UIvz6JFi9S9e3cdP35cu3btUrNmzfT2228XGHYAAIC1WaqHpzjo4QEAoOQplT08AAAAzhB4AACA5RF4AACA5RF4AACA5RF4AACA5RF4AACA5RF4AACA5RF4AACA5RF4AACA5RF4UKDMzEzFx8crMzPT16XARJxX6+LcWhPntfh4tMRfeLSEc3wu1sR5tS7OrTVxXl3j0RIAAAB/IfAAAADLI/AAAADLI/AAAADLI/AAAADLI/AAAADLI/AAAADLI/AAAADLK+vrAvxF7vyLdrvdx5X4l9zPg8/FWjiv1sW5tSbOq2u5n8nl5lFmpuW/7N+/X9HR0b4uAwAAFEFKSopq167tcj2B5y85OTk6ePCgwsLCZLPZfF0OAABwg2EYSktLU1RUlAICXI/UIfAAAADLY9AyAACwPAIPAACwPAIPAACwPAKPhezevVsfffSR7r//frVo0UJly5aVzWbTyy+/XKz9rl69WgMGDFDVqlUVGhqqxo0ba/z48crIyHB7H8uWLZPNZpPNZlOPHj2KVU9p44/ndenSpRo4cKCioqIUHBysGjVqqGvXrnrjjTeKVVNp4k/n9dy5c3rnnXfUoUMHVaxYUYGBgapZs6Zuvvlm/fDDD8Wqp7Qx+7wePnxY06ZN08MPP6x27dopODhYNptNI0aMcGv77du3a/DgwapZs6ZCQkJ05ZVX6sknn9SpU6eKVE+JZsAyHn30UUPSJT/jx48v8j4/++wzo0yZMoYko1atWkZsbKwRGBhoSDLatm1r/Pnnn5fdR3p6ulG/fn1HPdddd12R6ymN/Om85uTkGKNGjXLUULt2baNt27ZGvXr1jLJlyxpVqlQpck2ljb+c1z///NPo2LGj4/j16tUzWrVqZVSqVMmx7LXXXivOWy1VzD6vb731ltP93XfffZfd9ocffjBCQ0MNSUbVqlWNVq1aGeXKlTMkGTExMcbhw4eLVFNJRQ+PhURGRqpv37566aWX9M0332jgwIHF2t+ePXt03333KTs7W6+//rpSUlK0YcMG7dy5U1dddZUSExP19NNPX3Y/L7/8snbt2qX+/fsXq57Syp/O69ixY/XBBx+oadOmWrt2rVJSUrR27Vrt3r1bx48f15QpU4pVW2niL+f13//+t1avXq2qVavql19+0e7du7V+/XodOXJE8fHxkqTnn39eu3btKlZ9pYXZ5zU8PFw9e/bU2LFjtXDhQj3yyCNubZeWlqY77rhD6enpGj16tA4cOKD169dr37596tSpk5KTk3XfffcVq7YSx9eJC55z7733Fut/Fg8++KAhybj++usvWffzzz8bkozAwMAC/5ewbds2IygoyOjdu7cxZcoUenhM4Kvz+uuvvxplypQxqlataqSmphbp2HDNV+e1Q4cOhiTj3Xffdbrfli1bGpKMSZMmFamu0q645/VicXFxbvXwvP7664Yk4+qrrzaysrLyrdu7d69RtmxZQ5Kxfv16U+oqCejhgVOGYWj+/PmS5PR/Addcc40aNWqk8+fPa+HChS73MXLkSAUEBGjixIkerRfuKc55nThxorKzs/Xoo4+qWrVqXqkX7inOeU1PT5ckxcTEON33lVdeKUnKysoys2R42BdffCFJGjp0qMqUKZNvXZ06dRxjKefOnev12nyFwAOn9u3bp0OHDkmSOnXq5LRN7vI1a9Y4XT958mT99NNPeu6551z+YwrvKs55XbRokSSpb9++2rBhgx566CH17NlTAwYM0D//+U8dOXLEg5WjIMU5r82bN5ckrVq16pJtMjMztX79eklS27ZtTasXnpWVleU4b0X999uKCDxwaufOnZKk4OBgRUVFOW2TG2Jy2+Z19OhRPfPMM6pfv76eeeYZzxWKQinqeT18+LAOHjwom82m5cuXq127dpo0aZKWLVumL7/8UmPHjlWDBg20bNkyz78JXKI4f1+fffZZVahQQW+88Yb+/e9/68CBA0pPT9fGjRs1cOBA7dmzR3fffbc6dOjg2TcB0+zZs0fnz5+X5LrnrqB/v62KwAOnTp48KUmqVKmSy2eLVa5cOV/bvMaMGaMTJ05o4sSJCg4O9lyhKJSintfc3gObzaYnnnhC7dq104YNG5SZmamtW7eqZ8+estvtGjhwoFJSUjz8LnCx4vx9bdy4sX7++Wf17NlTTz75pGrXrq1y5copNjZWv/zyi9577z1NnTrVs28Apsp7jnPP+8UK+vfbqgg8cCp3zo6goCCXbXKDTO4YgFzff/+9ZsyYoVtvvVW9evXyXJEotKKe1z///FPShYfsVqhQQV999ZViY2MVFBSkxo0ba+HChYqKipLdbtfbb7/tuTcAp4rz91W6cEksNTVVhmEoKipKLVu2VIUKFRx33m3evNkzhcMj8s655OrPREF/HqyKwAOnQkJCJF2YkMyVzMxMSVJoaKhjWUZGhkaNGqUKFSrorbfe8myRKLSintfc7SRpyJAhl/yvMTQ0VKNGjZIkffvtt6bVC/cU9bxK0owZM9S/f38dOHBACQkJOnDggJKSknT8+HGNGzdOGzZs0N/+9jft3r3bc28Apsr799XVnwlXfx6sjMADp3K/0E6dOiXDMJy2ye0Kzfvl99prr2nXrl2Ki4tT7dq1PV8oCqWo5zXv60aNGjnd7uqrr5Z0YfwAvKuo5/X8+fN64oknZBiG3n77bXXp0sWxLigoSOPHj9f111+vtLQ0TZgwwYPvAGbKe45dXbJy9ufB6gg8cKpBgwaSLvwv4ODBg07bJCcn52srSUlJSZKk119/XTVq1Mj38+ijj0qSfvrpJ8cyxnt4V1HPa7169Rxd4K7GZOUuz87ONq1euKeo53Xnzp1KTU2VJF133XVOt8u9fXndunWm1QvPqlevngIDAyX9/3m/mLM/D1ZH4IFTderUUY0aNSRJP//8s9M2ucvbt29/ybqjR48qNTU134/dbpd0oYs1dxlfjt5V1PNapkwZx23Jl/sHtFatWqbVC/cU9bympaVddt+5PUaFeXYefKts2bJq1aqVpKL9+21VBB44ZbPZdPPNN0u6MJ/OxVatWqUdO3YoMDAw3yMjFixYIMMwnP7kPnbguuuucyyrV6+eV94PLijqeZWk22+/XZI0a9Ysxy2veeXeydO9e3ezy8ZlFPW8XnnllY67ur7//nun+86daqBhw4Zmlw0PuuWWWyRJn3766SX/sdy3b5/jvBb30RclCYGnlHv77bdVr1493XnnnZese+qppxQUFKQlS5bojTfecPxPb+/evRo+fLgkacSIEY7/WcJ/eOK8jhgxQtHR0dqzZ48effRRx2DI7OxsjR07VklJSQoKCtKYMWM8/O5KL7PPa2RkpONOyscee0w//vijY925c+f0j3/8Q0uXLpUk3XPPPR57X6VdQee1qEaNGqXIyEht375djz/+uOM/KcePH9egQYOUlZWl3r17q3Xr1qYd0+959UEW8KiVK1caVapUcfwEBwcbkoxy5crlW75v3z7HNrnPZenSpYvTfU6dOtUICAhw+vTl1q1bG2fOnHG7Pp6lVTT+dF4TExON8PBwQ5JRuXJlo23btkbVqlUNSUaZMmWMTz/91BMfgSX5y3nds2ePUadOHcdTuGvVqmW0bNnSCAsLcyy7//77PfUxWI7Z53Xfvn35tst9+nlwcHC+5StXrrxk22XLlhkhISGOp6W3bt3a8bT0evXqGYcOHfLkR+F36OGxkPPnz+v48eOOn9zbDs+ePZtveWHGzQwZMkQ//fST+vbtq/T0dG3btk0xMTGKj4/XypUrVb58eU+9HfzFn85rmzZttHnzZo0YMULly5fXxo0bJV3oPl+1apXuvffeYr/f0sJfzmvdunW1adMmxcXFKTY2VqdPn9aWLVsUEhKi3r17a968efrwww9Ne99WZ/Z5zc7Ozrdd7rw5mZmZ+ZY7u8x83XXXad26dbrzzjtls9n066+/qnr16nr88ce1YcOGUtc7bzMMF/cwAgAAWAQ9PAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAAAwPIIPAB86tNPP5XNZpPNZtOePXv8Zl++PAYA8xF4ALiUkJDg+HKPj4/3dTkAUGQEHgB+jR4VAGYg8ADwqaFDh8owDBmGoXr16vnNvrzpzJkzOnDggPbv368DBw4oNTVVp06dKtST0gEUrKyvCwCA0iYnJ0f/+9//9L///U8JCQk6ffq003YBAQGqWrWqGjZsqDZt2ujWW2/VNddc4+VqAWughwcAvGjlypVq27atpk6dqhYtWmj06NGqXLmyY32TJk2UmpqqzMxMHThwQNOnT1dUVJTeeustderUSUOGDFFOTo4P3wFQMhF4ABRZfHy8Y3yNJGVkZOiNN95Qq1atFBYWprCwMLVr104TJ05UVlaW0324GqOTO2B62LBhjmVXXHGFo23uT0JCwmX3ldeWLVv08ssvq1evXqpdu7aCg4NVoUIFNWjQQPfee69++eWXYn8urvz3v/9VfHy85s2bp2+//VYvvviiXnrpJfXo0cPRZuTIkapWrZqCgoJUo0YN9ezZU7Nnz3Z8DtOnT9cnn3zisRoBq+KSFgBTpKamqlevXtq0aVO+5YmJiUpMTNSSJUu0YMECBQT47v9ZCQkJ6tat2yXLz507p127dmnXrl2aNm2ann32Wb366qumH79FixYaMWKEypQpk2/5li1bHK/79evndNvhw4drypQpkqQlS5ZoxIgRptcHWBmBB4ApbrnlFm3fvl2jR49Wv379FBERod9++03jx4/X9u3btWjRIn300UcaOXKkW/tr27atfv31Vy1cuFDjxo2TJH333XeKiorK1+6KK65wu8asrCyVL19effr0Uffu3dWoUSOFh4fryJEj2rp1q959913t3btXEyZMUMOGDfP1LpmhQ4cOlyw7fPiwtm/fLkm6+uqrXQ62rlatmuN1enq6qXUBpQGBB4Apcntxunbt6ljWqlUr9erVS40bN1ZqaqomTZrkduApX768mjZtqnXr1jmWNWzYsFh3X7Vs2VL79+9XpUqVLlnXq1cvPfzww+rbt6+WLl2qF198UUOGDLmkN8Zs3333neP1jTfe6LLdkSNHHK8LE/IAXMAYHgCmeOSRR/KFnVwRERGOnpLNmze7vCPJGyIjI52GnVxBQUF64403JEl79+7Vxo0bPV7TV1995XhdUOBZvny543Xnzp09WhNgRfTwADDF4MGDXa5r3bq14/Xu3bvVsmVLL1R0eZmZmUpNTdWZM2ccdz4ZhuFYv2nTpny1m+38+fOOHp6wsDCXQSYtLU3vv/++pAuXvW666SaP1QRYFYEHgCkaNWrkcl1ERITjdVpamjfKcenPP//Uu+++q9mzZ2vr1q0FTu537Ngxj9ayYsUK2e12SVLPnj0VGBh4SZucnBwNHz5chw4dUmhoqKZNm+a0HYCCEXgAmKJcuXIu1+W9M8uXswfv2bNH3bt31+7du91q7+nBwYsWLXK8dnY5KzU1VUOHDtW3336riIgIzZ07V23atPFoTYBVMYYHQKlxzz33aPfu3bLZbBo+fLiWLFmilJQUZWRkOB5JkTeQ5b285QmLFy92vO7du7fjmGvXrtXjjz+u+vXr68cff9TIkSO1detWp7fUA3APPTwASoUdO3Zo5cqVkqTnnntOr7zyitN2J0+e9Eo927ZtU3JysqQLd4+FhISob9++OnTokCpUqKCWLVvq448/Vt++fVW+fHmv1ARYGYEHgF/LncW5uLZu3ep4feedd7psl/c2eE+6+HJWREREvh4fAObikhYAvxYSEuJ4nZmZWeT95H20xdmzZ122++CDD4p8jMK43PgdAOYi8ADwazVr1nS8/uOPP4q8nwYNGjheT5061Wmb999/XwsWLCjyMdx17NgxxzO7qlSp4nQGZgDm4pIWAL8WGxurkJAQZWRk6B//+IfKli2revXqOe78qlWrlkJDQ93aT9OmTbVlyxa9//77OnXqlAYPHqyaNWsqJSVFn332mebOnatOnTrp559/9uh7+vzzzx2Dox988EGPz+YMgMADwM+FhYVp9OjRev3117Vhwwb16tUr3/rly5c7neH5YjabTdOnT1f37t118uRJzZo1S7NmzcrXplmzZpozZ84lz+sy08mTJ/Wvf/1L0oUZk5966imPHQvA/+OSFgC/N2HCBH300Ufq3LmzIiIiitwj0rJlS23cuFGjRo1S3bp1FRgYqIiICLVr105vvvmm1q5dm+8SmtlWrlypzp07q0aNGvrPf/6jJUuWKCwszGPHA/D/bIanJ5oAAADwMXp4AACA5RF4AACA5RF4AACA5RF4AACA5RF4AACA5RF4AACA5RF4AACA5RF4AACA5RF4AACA5RF4AACA5RF4AACA5RF4AACA5RF4AACA5RF4AACA5RF4AACA5f0fY0vRr1xnsD4AAAAASUVORK5CYII=",
+      "text/plain": [
+       "<Figure size 640x480 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "# plot\n",
+    "plt.scatter(P_ratios, log_t_insts, c='silver', s=10)\n",
+    "plt.tick_params(labelsize=16,size=8,direction='in')\n",
+    "plt.xlabel(r'Initial $\\mathcal{P}$', fontsize=20)\n",
+    "plt.ylabel(r'$\\log_{10}(t_\\mathrm{inst}/P_1)$', fontsize=20)\n",
+    "plt.title(r'Model (1)', fontsize=20)\n",
+    "plt.xlim([1.028, 1.112])\n",
+    "plt.ylim([-0.4, 4.5])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1f9aff58",
+   "metadata": {},
+   "source": [
+    "Model (1) predicts that most widely-spaced systems will survive for the full integration time. Increasing the number of systems and/or the maximum integration time requires computational resources (or a great deal of patience)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "dcd6a2b4",
+   "metadata": {},
+   "source": [
+    "# Model (2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "id": "a66cbf0f",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "100%|████████████████████████████████████████████████████████████████████████████████████| 25/25 [04:01<00:00,  9.66s/it]\n"
+     ]
+    }
+   ],
+   "source": [
+    "P_ratios = []\n",
+    "log_t_insts = []\n",
+    "\n",
+    "# predict instability times for 25 systems\n",
+    "for _ in tqdm(range(25)):\n",
+    "    P_ratio = np.random.uniform(1.03, 1.11)\n",
+    "    sim = initialize_sim(P_ratio)\n",
+    "    \n",
+    "    # get indices of bordering first-order MMRs\n",
+    "    j1 = int(np.floor(1/(1 - 1/P_ratio)))\n",
+    "    j2 = int(np.ceil(1/(1 - 1/P_ratio)))\n",
+    "    \n",
+    "    res_inds1 = [[j1, 1, i, i+1] for i in range(1, 5)] # all j1 resonances\n",
+    "    res_inds2 = [[j2, 1, i, i+1] for i in range(1, 5)] # all j2 resonances\n",
+    "    \n",
+    "    pham = construct_Hamiltonian(sim, res_inds1 + res_inds2)\n",
+    "    t_inst = predict_t_inst(pham)\n",
+    "    \n",
+    "    P_ratios.append(P_ratio)\n",
+    "    log_t_insts.append(np.log10(t_inst))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "id": "dceb3a56",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle \\left( -0.4, \\  4.5\\right)$"
+      ],
+      "text/plain": [
+       "(-0.4, 4.5)"
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAHhCAYAAABqR2T2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABHUUlEQVR4nO3dd3hUVf7H8c8E0oAECKEFBjAUkQ6hiixFEJCmYAWlCYINxa7wExRdUXdXUBZsLE3Kikq1UgIKKAQSQJqCoYQWIJQJkgSS3N8fMbOJmYRJMpOZ3Lxfz5PHyZ0z935nLjIfzj33HIthGIYAAABMzMfTBQAAALgbgQcAAJgegQcAAJgegQcAAJgegQcAAJgegQcAAJgegQcAAJgegQcAAJgegQcAAJgegQcAAJgegQdAsTN58mRZLBZZLBa3HaNLly6yWCzq0qWLS/Y3YsQIWSwWPfLIIy7Z318ZhqGmTZvKYrFozpw5bjkGUJwReIASYsOGDfaQYLFYFBQUpCtXrlz3dUlJSSpfvny2127YsMH9BZvI9u3bNW/ePPn5+emll15y2ObXX3/Vu+++qzvuuEM33HCDAgMDVaZMGd1www2699579dVXXymvpQ8tFosmTJggSXr55Zd1+fJlt7wXoLgi8AAl1OXLl7V8+fLrtluxYoVsNpv7CzKxCRMmyDAMjRgxQrVq1crx/LBhw9SwYUM9/fTTWrFihY4cOaLk5GQlJSXpyJEj+uyzz9S3b1/17t1bFy9ezPU499xzjxo2bKjTp0/r/fffd+M7AoofAg9QAgUEBEiSFixYcN22mW0yX4P8iYqK0vfffy9JeuaZZxy2OXHihCQpJCREDz/8sBYtWqQtW7Zo27Zt+vDDD3XjjTdKkr777jv169dP6enpDvfj4+Ojp556SpL0z3/+U0lJSS5+N0DxReABSqD+/ftLktasWaPTp0/n2u7MmTP2L+sBAwYUSW1mM23aNElS27ZtVb9+fYdtatasqQ8//FAnT57Uhx9+qPvvv18dOnRQmzZt9PDDDys6Olq33HKLJGnTpk1auHBhrse7++675evrq4SEhDzbASUNgQcogW677TZVq1ZNaWlpWrx4ca7tFi9erNTUVFWtWlU9evQowgrN4dKlS/riiy8kSQ888ECu7ebOnauHH35Y/v7+Dp8vU6aMZs2aZf/9888/z3VfISEh6tWrlyRp9uzZBSkbMCUCD1AClSpVSvfff7+kvC9rzZ8/X5I0ePBglSpVyql9X716VTNnzlTXrl1VuXJl+fn5qVq1arr99tv16aef5no5Jqvjx4/rscceU3h4uAICAhQWFqb+/ftr7dq1TtWQ6cqVK5o2bZq6du2qqlWrys/PT1WqVNFtt92mOXPmKC0tLV/7y68VK1YoJSVFkjRo0KBC7atJkyYKDQ2VJP3+++95ts081s8//6yjR48W6riAWRB4gBLqwQcflCTFxMRo7969OZ7ft2+foqOjs7W9nqNHj6pFixZ67LHHtGHDBp07d07Xrl1TfHy8vvnmGz344IPq3Lmzzp8/n+s+Nm7cqMaNG2vmzJk6fPiwUlJSdOrUKa1atUo9evTQa6+95lQtUVFRatCggcaPH68NGzbozJkzunbtms6ePas1a9Zo5MiRuvnmmxUfH+/U/goiMjJSkmS1WhUWFlbo/V29elVSxlidvLRv397+eOPGjYU+LmAGBB6ghGrZsqWaNGkiyXEvT+a2xo0bq2XLltfd3+XLl9WtWzft379fknTHHXdo5cqV2r59u5YuXarOnTtLyhiD0rdvX4e9K0eOHFG/fv1ks9nk4+OjsWPHau3atYqKitLs2bNVv359TZo0SV999VWetfzyyy/q2rWrTpw4oSpVqmjSpElau3atYmJi9N133+mxxx5T6dKltW3bNg0YMEDXrl277vsriB9//FGS1KZNm0LvKyYmxn63XMOGDfNs26BBA1WoUCFbDUCJZwAoESIjIw1JhiRjzpw5hmEYxltvvWVIMmrWrGmkpaXZ26anpxtWq9WQZEydOtUwDMOYM2eO/fWRkZE59v/ss8/an584cWKO59PT040hQ4bY28ycOTNHm4EDB9qfX7RoUY7nbTab0bx5c3sbR3+FpaenG82aNTMkGc2bNzfOnj3r8PP45ptvDB8fH0OS8cknn+R4vnPnzoYko3Pnzg5ffz3x8fH2GqdMmVKgfWR111132fe3dOnS67bv2rWrIclo1KhRoY8NmAE9PEAJNmTIEPn4+Oj48ePZLn1s2LBBcXFx8vHx0ZAhQ667n5SUFH3yySeSpEaNGmny5Mk52lgsFs2cOVOVKlWSJM2YMSPb86dOndKKFSskSX379rWPMcoqKChIH330UZ61fPXVV9q9e7ekjDFImeNe/qpXr1666667JMktMxMfP37c/rhKlSqF2tcXX3xhH6gcERHh1HigzGNmrQMoyQg8QAlWo0YNde3aVVL2y1qZj7t06aKaNWtedz87duywT4g3fPjwXAc4BwcH65577pGUMUbo1KlT9uciIyPtl7lGjBiR67Hatm2rxo0b5/p8Zmi68cYb1axZszzr/tvf/iYpY7yPqwcwnz171v64YsWKBd7PgQMH7J9HYGCg5s+f79SSGiEhIZIkm81mH/sDlGQEHqCEGzp0qKSMW52TkpKUlJRkv5Xa2cHKe/bssT9u165dnm2zPp/1db/88ov98fXGvLRt2zbX57Zv3y4pY6mGrMthOPp5/PHHJWUMBs5rIHVBZN1fQQPPyZMn1bt3byUmJspisWj27Nlq1KiRU6/NesyEhIQCHR8wEwIPUMINHDhQZcqUUWJiolasWKHly5fLZrMpMDDQ6Vups365V61aNc+21apVc/i6Cxcu2B9f7xJQXsc4c+ZMnq/NjTPriuVH1pmpCzLj8fnz53XbbbfpyJEjkqTp06c7vMyXm6zHDAwMzPfxAbMp7ekCAHhWuXLldOedd2rhwoVasGCBfYHKO+64Q0FBQfne3/Uutxi5LICZdXtB9yHJfmmqY8eO+uCDD/LcT1auuG08q8qVK9sf57f3KDExUb169bJPFzBlyhQ98cQT+dpH5jF9fX3td2wBJRmBB4CGDh2qhQsX2peRkJy/nCX9b7yIJJ0+fVoNGjTItW3WeW+yvi7r4/j4eFmt1lz3kVcvTqVKlRQfH6+zZ8/ab7v3hKyBJ2vv1fUkJSWpX79+ioqKkiQ999xzmjhxYr6Pn3nM3AZtAyUNl7QA6NZbb1X16tWVmppqX0ritttuc/r1WYPF1q1b82y7bds2h69r2rSp/XHml31u8no+c86g3377zaOzDNetW1dlypSx1+KMa9euadCgQfY75saOHau33367QMfPPGbWzxUoyQg8AFSqVCk9+OCD8vf3l7+/vx544AGnl5KQMm6VzrxsMm/evFzveEpMTNRnn30mKeP29erVq9uf69q1q/2Y8+bNy/VY27dvzzbY+a8yF0aVVOCw4AqlS5e2z3h8vQAnZVyKGzx4sL755htJGT1sM2fOLNCxbTabfv31V0lSp06dCrQPwGwIPAAkSW+99ZaSk5OVnJysf/zjH/l6rb+/v0aNGiVJ2rt3r1599dUcbQzD0OOPP65z585Jkv0OqUzVq1e3r8i+cuVKezDK6vLly3r44YfzrGXQoEG66aabJEmzZs267gKae/bs0apVq/JsU1CZYWP37t32NbUcMQxDo0ePts+1M2jQIM2ZM8ep288d2b59u32cE4EHyEDgAeASr7zyisLDwyVlDLIdOHCgVq9erejoaH3xxRfq1q2bfTHSDh06OAwu//znP+0DpQcPHqzHHntMkZGR2rFjh+bMmaOIiAjFxMSodevWudZRqlQp/fe//1W5cuVkGIZGjRqlXr16af78+dq6dauio6P17bff6s0331THjh3VtGlTt6031adPH0kZt73ntcTDs88+a5/8sEmTJnr55Ze1f/9+7dmzJ9efvKxbt05Sxq3pHTp0cNG7AYo3Bi0DcImgoCCtW7dOvXv31oEDB7Rs2TItW7YsR7uOHTtq5cqVDi+Z1alTRytXrlT//v2VmJiomTNn5risM2nSJEn/m2/HkaZNm2rz5s266667dPDgQX333Xf67rvvcm0fHBzs7NvMlzZt2qhhw4Y6cOCAFi1apO7duztslznvkZTR4xQREXHdfed1p9rixYslSffcc4/8/PzyWTVgTvTwAHCZOnXqaNeuXZoxY4Y6d+6sSpUqydfXV1WrVlWvXr20YMEC/fDDD9nuyPqrLl26aO/evXrkkUdUu3Zt+fn5qWrVqurTp4++/fZbh8tWONKsWTPt27dP8+bN0x133CGr1aqAgAD5+fmpevXq6tKliyZOnKgdO3bolVdecdEnkNOjjz4qKSPUJCcnu+04mX766ScdPnw427EBSBYjr38mAAAK5fLly6pTp44SEhK0YMECPfDAA2493qhRozR79mz16NEj2zQDQElHDw8AuFG5cuX0/PPPS5LeeOMNpaenu+1Yx44ds4+TcjRwHCjJCDwA4Gbjxo1T7dq1deDAAYd3n7nKm2++aZ/Lh8HKQHYMWgYANwsICNCnn36qtWvXKjU11S3HMAxDtWvX1qRJkzRy5Ei3HAMozhjDAwAATI9LWgAAwPS4pPWn9PR0nTx5UkFBQQWe3RQAABQtwzCUmJiosLAw+fjk3o9D4PnTyZMn81ydGQAAeK+4uDjVrFkz1+cJPH/KnM4+Li7ObbOuAgAA17LZbLJarfbv8dwQeP6UeRkrODiYwAMAQDFzveEoDFoGAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmZ+rAM3HiRFksFlksFr3++uueLgcAAHiIaQPP/v379c4773i6DAAA4AVMGXgMw9CYMWPk6+urbt26ebocAADgYaYMPLNnz9aPP/6oV155RVar1dPlAAAADzNd4Dl79qxeeOEFNWrUSOPHj/d0OQAAwAuU9nQBrjZ+/HidP39eX375pXx9fT1dDgCgAGw2m5KSkhQYGKjg4GBPl+OQoxqzbpPkkvdwveM42rezn587PmdvPXemCjzr1q3TwoUL9cADD6hz586eLgcAUACxsbGKi4uz/261WhUeHu7BinJyVKOkbNuyKuh7cOY4f923s5+fOz5nbz53pgk8ycnJGjt2rMqXL69//OMfBd6PzWZzqp2/v7/8/f0LfBwAQE42my1HaIiLi1NoaKjX9BbkVmNeCvIenD1O1n07+/m543P29nNnmjE8r7/+ug4dOqQ33nhDVatWLfB+rFarypcvf92fN99804XVAwCkjEtA+dnuCQWtJb+vy0/7zLbOfn7u+Jy9/dyZoocnc86dVq1a6ZFHHinUvuLi4pxKovTuAIDrZY59cXa7JxS0lvy+Lj/tM9s6+/m543P29nNnih6eRx99VKmpqZo1a5Z8fAr3loKDg536IfAAgOsFBwfnmE7EarV6xSWRTLnVmNc0KAV5D84eJ+u+nf383PE5e/u5sxiGYXi6iMKqUKGCLl++rNDQ0BzPXbp0ScnJySpXrpzKli0rq9WqqKioHO1sNpvKly+vS5cuec3JAYCSylvv9MmKu7SKbp/XO54z39+mCTyXLl1yqm3t2rV15MiRHNsJPAAAFD/Ofn+b4pLWxYsXZRiGw59hw4ZJkqZMmSLDMByGHQDA/9hsNsXHxzt91ypQHJhi0DIAwDW8eR4VoDBM0cMDACi83OZRoacHZkDgAQBI8v55VIDCMP0lrblz52ru3LmeLgMAvJ63z6MCFAY9PAAASd4/jwpQGKbv4QEAOC88PFyhoaFePwcOkF8EHgBANpkzygNmwiUtAABgegQeAABgegQeAABgegQeAABgegQeAABgegQeAABgegQeAABgegQeAABgegQeAABgegQeAABgegQeAABgegQeAABgegQeAABgegQeAABgegQeAABgegQeAABgeqU9XQAAeJLNZlNSUpICAwMVHBzs6XIAuAmBB0CJFRsbq7i4OPvvVqtV4eHhHqwIgLtwSQtAiWSz2bKFHUmKi4uTzWbzUEUA3InAA6BESkpKytd2AMUbgQdAiRQYGJiv7QCKNwIPgBIpODhYVqs12zar1crAZcCkGLQMoMQKDw9XaGgod2kBJQCBB0CJFhwcTNABSgAuaQEAANMj8AAAANMj8AAAANMj8AAAANMj8AAAANMj8AAAANMj8AAAANMj8AAAANMj8AAAANMj8AAAANMj8AAAANMj8AAAANMj8AAAANMj8AAAANMj8AAAANMj8AAAANMr7ekCAMBTbDabkpKSFBgYqODgYE+XA8CNCDwASqTY2FjFxcXZf7darQoPD/dgRQDciUtaAEocm82WLexIUlxcnGw2m4cqAuBuBB4AJU5SUlK+tgMo/gg8AEqcwMDAfG0HUPwReACUOMHBwbJardm2Wa1WBi4DJubyQcvnzp3TiRMndPbsWSUkJCgwMFCVK1dW5cqVFR4eLh8fMhYAzwsPD1doaCh3aQElRKEDT2JiolasWKENGzboxx9/1KFDh3JtW7ZsWbVv316dOnVSnz591KpVq8IeHgAKLDg4mKADlBAWwzCMgrxwx44dmj59ur744gslJydLkpzZlcVisT++8cYb9dhjj2n48OEqW7ZsQcpwGZvNpvLly+vSpUv8BQgAQDHh7Pd3vgPPjh07NHHiRH3//feS/hdyqlevrjZt2igiIkJVqlRRSEiIKlasqKSkJJ0/f14XLlzQb7/9pqioKO3evVvXrl3LKMBiUUhIiJ577jk9+eST8vf3L+h7LhQCDwAAxY9bAs+IESO0YMECpaenS5JatWqlIUOGaNCgQapVq5bTxV29elU//PCDFi1apGXLlunSpUuyWCyqVauWFixYoFtuucXpfbkKgQcAgOLH2e/vfI0gnjdvnkqXLq3Ro0frwIED2r59u8aPH5+vsCNJfn5+6t69u/7zn//o9OnTmj9/vm688UYdPXpU69evz9e+AAAAridfg5YfeeQRvfTSS6pZs6bLCvD399cDDzygIUOGaOnSpUpLS3PZvgEAAKRCDFo2Gy5pAQBQ/LjlkpY3W758ucaMGaOIiAhVr15dfn5+qlChgm6++WZNnz5dV69e9XSJAADAQ0zTw3PLLbdo8+bN8vf3V1hYmCpVqqRTp07pxIkTkqSIiAitXbtWFSpUcPh6engAACh+SlwPz6hRoxQZGanExETFxsYqKipKx48f108//aSaNWtqx44dmjBhgqfLBAAAHuBVPTwrVqzQrl279Morr7h0v0uXLtU999yjsLAwe4/PX9HDAwBA8VMse3iWL1+uV1991eX7bdiwoSTpypUrLt83AADwfl4VeNzlp59+kiTW7gIAoIRy+WrpWc2fPz9f7fNaeDS/0tLSdOrUKa1cuVIvvviiypYtqzfffNNl+wcAAMWHWwPP8OHDsy0Wej2GYeSrvSPTpk3T+PHjs2274447NGXKFDVp0uS6r7fZbE4dx9/f32PrfgEAgPxxa+Dx8/NTWFiYxowZ41T7pUuXKiYmplDHrFGjhjp27Khr167p6NGjio+PV2RkpBYvXqzXXntNpUqVyvP1VqvVqeNMmjRJkydPLlStgDvZbDYlJSUpMDCQgfji8wBKOrfepdWmTRsdO3ZM8fHxTrUfMWKE5s+f79LlJbZu3aoxY8Zo165dGjt2rGbNmuWwXeYo77i4OKf+MqSHB94sNjZWcXFx9t+tVqvCw8M9WJFn8XkA5uUVd2lFRETo3Llz2f6iKWrt2rXT119/LX9/f3300Uc6evRonu2Dg4Od+iHswFvZbLYc/8/FxcU5fbnWbPg8AEhuDjydOnVSrVq1dPDgQafa33LLLRo6dKjL6wgLC1OLFi2Unp6uXbt2uXz/gDdJSkrK13az4/MAILl5DM+QIUM0ZMgQp9s/9NBDeuihh9xSS2pqarb/AmYVGBiYr+1mx+cBQCpE4Pn999+1aNEiHThwQGlpaapevbr+9re/6fbbb/e6yz1Hjhyx9+w0b97cw9UA7hUcHCyr1ZpjzEpJHajL5wFAKmDg+fDDDzVu3LgcvSXvvfeerFarPvjgA/Xq1cslBTpjx44dWrlypYYNG5ZjIOK3336r8ePHKzU1Vbfffrvq1q1bZHUBnhIeHq7Q0FDuSvoTnweAfN+lFRUVpQ4dOig9PT3XNqVLl9bnn3+u/v37F7pAZ2zYsEFdu3aVJFWrVk01a9bU1atXdezYMV28eFFSxh1jX3/9tUJDQx3ug7W0AAAoftx2l9b777+v9PR0WSwW9e3bV8uWLdPWrVv15ZdfauTIkfL19VVqaqpGjhypc+fOFepNOKt58+aaPn26+vfvr7Jly+rAgQM6cOCAAgMD1bt3b82ZM0dbtmzJNewAAABzy3cPT3h4uI4eParevXtr9erVOZ7/4Ycf1LNnT129elVTpkzRyy+/7LJi3YkeHgAAih+39fCcOnVKknKdPflvf/ubnnnmGRmGoS+++CK/uwcAAHC5fAeelJQUSVKdOnVybXP//fdLkn755RddvXq1YJUBAAC4SIEnHsxrTap69epJylix/OzZswU9BAAAgEu4ZablrPPwJCYmuuMQAAAATitw4LFYLE61y+v2dQAAgKJQ4JmWb7nlFjVv3lwtWrSw/zRq1EilS7t1tQoAAIB8K1A6MQxDFy5c0MaNG7Vx40b7dl9fXzVq1EgtWrSwb7t27VqhiwQAACiMfM/D88EHH2jnzp3auXOn9uzZoytXrmTf4V8udfn4+Khu3bpq1qyZmjVrpubNm6tZs2aqXbt24at3IebhAQCg+HH2+zvfgScrwzD066+/2gPQzp07tWvXLsXHx+c80F+CUFBQkJo2barmzZtrxowZBS3BZQg8AByx2WyswQV4sSIJPLmJj49XTExMtiB06NAhhwOYLRaL0tLSXF1CvhF4APxVbGxsjlXW/7pAMQDP8mjgcSQpKUm7du3KFoL27NmjpKQkAg8Ar2Oz2RQTE5Nje8uWLfk7AvAizn5/F9ktVYGBgWrfvr3at29v32YYhn777beiKgEAnJaUlJTrdgIPUPzkex6eHTt2uOzgFotFN954o8v2BwCuEhgYmK/tALxbvgNPmzZtVLNmTY0ZM0arV69WcnKyO+oCAI8KDg6W1WrNts1qtdK7AxRT+R7D4+OTkZEy77oKCAhQt27d1K9fP/Xt21dhYWGur7IIMIYHgCPcpQV4N7cNWj558qRWr16tVatWaf369fbr3JkBqEWLFurXr5/69euniIiIQryFokXgAQCg+CmSu7SSkpK0du1arV69Wl999ZVOnjyZsdM/w0+1atXUp08f9evXT927d/fqa98EHgAAih+P3Ja+Y8cOrVq1SqtXr1Z0dHTGAYrJpS8CDwAAxY/H5+Epbpe+CDwAABQ/Hg88WSUnJ2vt2rVatWpVrpe++vbtq0cffVTNmzd3dzkOEXgAACh+vGriwYCAAPXt21d9+/aVlHHpK7P3JyYmRqdOndInn3yiGjVqeCzwAAAA8yqymZazioiIUEREhCZNmmS/9LV69WqVKVPGE+UAAACTK7K1tI4dO6YKFSrk2d1ks9l08eJF1apVqyhKynFsLmkBAFC8OPv9ne+Zlgvqhhtu0PTp0/NsM2vWLN1www1FVBEAV7DZbIqPj5fNZvN0KQCQqyK7pGUYhq7XmVREnU0AXCQ2NlZxcXH2361Wq8LDwz1YEQA4VmQ9PM44fvy4goKCPF0GACfYbLZsYUeS4uLi6OkB4JXc2sPz2muvZft9w4YNDtulpaXp+PHjWrJkidq1a+fOkgC4SObcWo62Mw4OgLdx66DlzIVGpYw5d653qLCwMC1btkxt2rRxV0m5YtAykD82m00xMTE5trds2ZL/hwAUGa+YhycyMlJSxticbt26afjw4Ro2bFiOdqVKlVJISIgaNmyYLSQB8F7BwcGyWq05xvAQdgB4I7cGns6dO9sfT5o0SV26dMm2DUDxFh4ertDQUCUlJSkwMJCwA8BrFdk8PLkxDEOHDh1SYGCgatas6bE6uKQFAEDx43Xz8KxcuVIjR47UhQsX7NuOHDmipk2bqmHDhqpdu7aGDBmi9PT0oioJwF8wpw4AsyqyeXhmzZql48ePq2LFivZtTz31lPbt26du3bopISFBS5YsUbdu3fTQQw8VVVkA/sScOgDMrMh6ePbu3au2bdvaf7906ZK+/vpr3XvvvVq7dq22bdumm266SbNnzy6qkgD8iTl1AJhdkQWes2fPqnr16vbfN23apNTUVN1///2SJF9fX/Xo0UOHDh0qqpIAlyuul4TymlMHAMygyC5pBQcHKyEhwf77hg0b5OPjo06dOtm3+fr66o8//iiqkgCXKs6XhAIDA/O1HQCKmyLr4WnYsKFWrVql8+fP69KlS1qyZIlatWqVbUzP0aNHVbVq1aIqCXCZ4n5JKHNOnayYUweAmRRZD8+4ceN09913q0aNGvaenFdffdX+fFpamjZt2pRtnA9QXJhhmQV3zaljs9mYpweAxxVZ4Bk0aJD+/e9/2wcl33PPPRo5cqT9+XXr1unKlSvq1atXUZUEuIxZLgkFBwe7NJQU58t8AMzF4xMPegsmHkRh8eWeHWttASgKXrGWFlCSsMxCdma4zAfAPAg8gAu5+pJQcWaWy3wAzKFIA8/Vq1e1fPlyRUVF6eLFi0pLS8vRxmKxMPkgYAKspg7AmxTZGJ6jR4+qR48e+v3335XXIS0Wi8Mg5G6M4QHcg7u0ALiT143hGT9+vA4dOqQHH3xQI0eOVM2aNVW6NFfUALPjMh8Ab1BkiWP9+vW69dZbNW/evKI6JGBK9JgAQP4VWeBJT09Xy5Yti+pwgClx6zsAFEyRLS3RoUMH7d+/v6gOB5hOcV++AgA8qcgCz9SpUxUZGanPP/+8qA4JmAormgNAwRXZJa1Vq1apa9euuvfee9W5c2e1bNlS5cuXz9HOYrHo//7v/4qqLKDYYF4bACi4Irst3cfHuc4kbktHcVRUA4kZwwMA2XndbemRkZFFdSigSBVlCGH5CgAomCILPJ07dy6qQwFFJreBxKGhoW4LI8xrAwD5V2SDlgEz8vRAYpvNpvj4eO7UAoDrcFsPz8iRI2WxWPT3v/9dVatW1ciRI516HWtpoTjx5EBixvMAgPPcNmjZx8dHFotF+/fvV4MGDRi0DNPyRPCw2WyKiYnJsb1ly5b8+QVQonh80PLhw4clSTVq1Mj2O2A2nhhInNelNAIPAOTkksAzePBgNW/eXI888oj9L9vatWtna/PX3wEzKeqBxMzJAwD545JBy0uWLNHLL7+s48ePu2J3+WYYhjZt2qTnnntO7du3V4UKFeTn56ewsDANGjSIW+JhOsHBwbJardm2Wa1WencAIBdFdlv6X12+fFmHDh1SkyZNVLp04cpYv369unfvLilj7FC9evVUtmxZHTx4UF9++aW+/PJLTZw4UVOmTHFF6YBX8MY5eVjJHYC38tht6UePHlWrVq0UFBRU6H0ZhqF69epp5syZOnfunH799VdFR0crISFBL730kiTp9ddf1+rVqwt9LMCbBAcHq2rVql4RLmJjYxUTE6MDBw4oJiZGsbGxni4JAOwKFHjefvttRUZGumTuj6tXrxZ6H23bttX+/fv1yCOPqGLFivbtfn5++vvf/67evXtLkj7++ONCHwtATqzkDsDbFSjwvPjii+revbtCQkLUsGFD+/bIyEin/1WXnJwsKeM29MIKDg7O87JYjx49JEm//fZboY8FICdPT8AIANdToMEzpUuXVmpqqgzDyBYixo0bJ0kKCgpSs2bN1KJFC7Vo0ULNmzdX06ZN5efnZ2/7888/S5IqVKhQiPKdkxmuuIMFcA/uGgPg7QoUeBITExUdHa1t27YpKipKixYtkpQxlkbK6N7etGmTNm/ebH9NqVKldOONN6px48ZKTU3V6tWrZbFY1Lp1axe8jdwZhqGlS5dKkjp27OjWYwElVeZdY3+dgNEbxhYBgOSimZYzZ1VeunSpkpKStHPnTu3cuVO7du3SuXPnsh/wz0tYhmHIx8dH33zzjf2Skzt89NFHGjNmjPz8/LRv3z7VrVvXYbvMmRrj4uKc+kva399f/v7+ri4XKNa4SwtAUXN2pmWXBJ7KlSvr/Pnz+uWXX9SoUaNsz508eTJbAPr1119ls9kUHh6u8ePHq0+fPoU9fK6io6PVsWNHJScn6+2339Zzzz2Xa9vMD8xZkyZN0uTJk11QJQAAKKgiDTySdOrUKVWsWFEBAQGu2F2hHT58WB07dtSpU6c0ePBgffrpp3kOkKaHBwCA4qfI19KqXr26q3ZVaKdPn1aPHj106tQp9enTR3PnznX6brCiXiIAAAC4n8cmHnSX8+fPq0ePHvr999/VuXNnLV26VL6+vp4uCwAAeJCpAs/ly5d1++23a8+ePWrTpo1WrVrFbbEAACB/gSfz9m53OX78uLZs2VKg16akpGjAgAHaunWrGjdurG+//dYly1YAAIDiL1+B595771XTpk1dHnyOHTumRx55RPXq1dPatWvz/fq0tDTdd999Wr9+verWras1a9YoJCTEpTUCAIDiK1+DluvXr6+9e/fqvvvu03PPPafBgwdryJAhaty4cb4P/Mcff2jZsmVauHCh1q1bp7S0NPtK5/n12Wefafny5ZIy5gS6++67HbarXr2623upAACA98nXbempqal67733NHXqVJ07d85+51P9+vXVvn17tWnTRi1btlSVKlVUsWJFVaxYUUlJSTp//rwuXLig3377TVFRUdq2bZu2bdum5ORk++zMvXv31ltvvaUmTZrk+03MnTtXI0aMuG672rVr68iRIw6fc/a2NgAA4D3cOg/PH3/8oZkzZ+rf//63jh07lrGjfCwCmnnIUqVKacCAAXruuefUrl27/JbhUgQeAACKnyKZeDA9PV1r1qzRZ599psjIyFx7T7IKDAxU27Zt1adPHw0ePFhhYWEFPbxLEXgAACh+inymZUk6ceKEtmzZouPHj+vs2bM6f/68AgICVLlyZVWuXFlNmzZV69atvXJeHAIPAADFT5HPtCxJNWrUyHXAMAAAgKeYauJBAAAARwg8AADA9Ag8AADA9Ag8AADA9Ag8AADA9Ag8AADA9Ag8AADA9Ag8AADA9Fw68WBufv/9d82ZM0cbN27UwYMHdenSJUlS+fLlVb9+fXXp0kXDhg0r0ErpAAAA1+PSpSUcmTp1qiZNmqRr165JkkJDQ+1TP9tsNp07d06S5Ovrq1dffVUvvviiO8vJFUtLAABQ/Dj7/e3WS1qLFy/Wyy+/rAYNGmjJkiU6f/68zpw5o0OHDunQoUM6c+aMzp8/r8WLF6t+/fqaMGGClixZ4s6SAABACeTWHp527dopISFBu3btUtmyZfNsm5iYqBYtWig0NFRbt251V0m5oocHAIDixyt6ePbu3auBAwdeN+xIUlBQkAYOHKi9e/e6syQAAFACuTXw+Pn52QcoO8Nms8nPz8+NFQEAgJLIrYGnffv2WrJkiXbv3n3dtrt27dLixYvVoUMHd5YEAABKILeO4YmKilKnTp3k4+OjIUOGqEePHqpfv77Kly8vSbp06ZIOHjyo77//XosWLVJ6ero2bdqk1q1bu6ukXDGGBwCA4sfZ72+335a+YcMGjRo1SrGxsbJYLA7bGIah8PBwffLJJ+rSpYs7y8kVgQcAgOLH2e9vt0882KVLF/36669av369NmzY4HDiwc6dO+vWW29VqVKl3F0OAAAogdzew1Nc0MMDAEDx4xW3pQMAAHgDrwo8K1as0GuvvebpMgAAgMl4VeBZvny5Xn31VU+XAQAATMarAg8AAIA7uPUurfnz5+er/aFDh9xUCQAAKMncGniGDx+e69w7jhiGka/2AAAAznBr4PHz81NYWJjGjBnjVPulS5cqJibGnSUBAIASyK2Bp2nTpjp27JheeOEFp9ofOHCAwAMAAFzOrYOWIyIidO7cOcXFxbnzMAAAAHlyaw9Pp06d9N133+ngwYOyWq3XbX/LLbe4sxwAAFBCsbTEn1haAgCA4oelJQAAAP5E4AEAAKbn1jE8QH7ZbDYlJSUpMDCQS4sAAJdxSeA5duxYvl9jsVgUEBCg8uXLy8/PzxVloJiLjY3Ndkef1WpVeHi4BysCAJiFSwLPDTfcUKjX16xZU+3bt9fw4cPVu3dvV5SEYsZms+WYviAuLk6hoaH09AAACs0lY3gMwyjUT1xcnD7//HP17dtXPXv21IULF1xRFoqRpKSkfG0HACA/XNLDM2fOHEnSBx98oK1btyogIEA9e/ZU69atVblyZUnS2bNntX37dn333XdKSUlR27Zt9fDDD8tms2nPnj1asWKFzp07p7Vr1+qOO+7Qxo0bXVEaionAwMB8bQcAID9cNg/P2LFj9fHHH6t///768MMPVaVKFYftzpw5o4cfflirVq3SiBEj9Mknn0iSkpOTNXbsWM2fP18Wi0WLFi3Svffe64rSnMI8PJ7HGB4AQH45+/3tksCzfPlyDRw4UDfffLN+/PHH6654bhiGbrnlFv38889avHix7rnnHklSenq62rZtq5iYGPXv31/Lli0rbGlOI/B4B+7SAgDkR5FOPDhz5kxZLBY9+eST1w07UsYdWk899ZQMw9BHH330v2J8fDR69GgZhqHt27e7ojQUM8HBwapatSphBwDgUi4JPLt375Yk1atXz+nXZLb95Zdfsm1v1qyZJCkhIcEVpQEAALgm8NhsNkkZA5Odldk2MTEx2/bMQaq+vr6uKA0AAMA1gSdzJfSFCxc6/ZoFCxZIkmrVqpVt+5kzZyTJfncXAABAYbkk8PTv31+GYejTTz/VtGnTrtv+3Xff1cKFC2WxWNS/f/9sz23dulWSVLt2bVeUBgAA4Jq7tBISEnTTTTfZx920atVKQ4cOVUREhP329DNnzmj79u1asGCBoqOjZRiGqlSpon379ikkJMS+r8aNG+vAgQN644039OKLLxa2NKdxlxYAAMWPs9/fLpl4sFKlSlqzZo169eql+Ph4RUdHKzo6Otf2hmGoWrVq+vbbb7OFndjYWLVt21Zt27bVoEGDXFEaAACA6yYelKSLFy9q0qRJmj9/vi5duuSwTfny5TV06FBNnjxZFStWdNWhC40eHgAAip8inXjwr1JSUrRjxw7t2bPHvi5WxYoV1bhxY7Vu3Vr+/v6uPmShEXgAACh+PBp4iiMCDwAAxU+RzrQMAADgzVwyaNmR+Ph47dmzR+fPn5ckhYSEqEmTJqpataq7DgkAAOCQSwNP5tpYM2bM0L59+xy2adSokZ544gmNHj3aqXW3AAAACstlY3guXLigfv366aeffpKUEX4cHvDPkHPzzTdr1apVqlChgisOX2iM4QEAoPgp0nl4DMPQgAEDtGXLFkkZ8/Lcc889ateunapVqybDMBQfH69t27bps88+07lz57RlyxYNGDBAGzdudEUJAAAAuXJJD8/ChQv14IMPymKxaPDgwZo5c6aCgoIctr18+bIee+wxLViwQBaLRZ9++qnuv//+wpagw4cPa+3atdq2bZu2bdumvXv3Ki0tTVOmTNHEiROv+3p6eAAAKH6KtIdn0aJFkqTOnTvbFwXNTbly5TRv3jwdO3ZMGzdudFngmT59uqZPn17o/QAAAPNxyW3p0dHRslgsevzxx51+zRNPPCFJiomJcUUJCg0NVd++ffXaa6/pm2++YWkKD7DZbIqPj5fNZvN0KQAAZOOSHp7MW89vuOEGp1+T2TbztYX118tWS5Ysccl+4ZzY2FjFxcXZf7darQoPD/dgRQAA/I9LenjKly8vSTp58qTTr8lsy3iZ4s9ms2ULO5IUFxdHTw8AwGu4JPA0adJEkjRnzhynX/Of//wn22tRfCUlJeVrOwAARc0lgeeuu+6SYRhatmyZJk+enOscPFLGLeyTJ0/WsmXLZLFYdPfdd7uiBHhQYGBgvrYDAFDUXDKGZ/To0ZoxY4YOHDigKVOm6IsvvtDw4cPVrl07Va1aVRaLRadPn9bWrVs1b9487d27V5LUsGFDjR492hUluIyzl2H8/f29ctV3TwgODpbVas0xhofLlQAAb+GSwOPr66tvvvlG3bp10+HDh7Vv3z49//zzubY3DEPh4eH65ptvVLq025bzKhCr1epUu0mTJmny5MnuLaYYCQ8PV2hoqJKSkhQYGEjYAQB4FZeljdq1a2v37t2aPHmyZs+erYsXLzpsV6FCBY0aNUqvvPKKypUr56rDu0xcXJxTX9b07uQUHBxM0AEAeCWXdq+ULVtW77zzjt544w3t2LHD4WrpERER8vPzc+VhXYovbQAAzMct15P8/PzUoUMHdejQwR27BwAAyBeX3KUFAADgzfLVw3Ps2DG3FFGrVi237BcAAEDKZ+DJz9IRzrJYLEpNTXX5fgEAADJZjLxmCfwLHx/XXwGzWCxKS0sr9H42b96sAQMG2H+/fPmyUlJSVKZMmWwT4MXExDi89dzZ5eUBAID3cPb7O189PPlZOqKoXbt2TQkJCTm2X7lyRVeuXLH/7opwBQAAipd89fCYGT08AAAUP85+f3OXFgAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD3TBZ6vv/5a3bt3V0hIiMqWLatWrVrp/fffV3p6uqdLAwAAHmKqwDN16lT16dNH69atU8WKFVWvXj3t2rVL48aN05133knoAQCghDJN4Pnpp5/08ssvy8fHR4sWLdLvv/+uXbt2KTo6WlWrVtXKlSv1r3/9y9NlAgAADzBN4Hn99ddlGIZGjRql+++/3769efPm9qAzdepUXbt2zVMlAgAADzFF4LHZbFq7dq0k6aGHHsrx/N13363g4GAlJCQoMjKyqMsDAAAeZorAExMTo6tXryogIECtWrXK8byvr6/atGkjSdq6dWtRlwcAADystKcLcIWDBw9KkmrVqqXSpR2/pfDwcK1bt87eNjc2m82pY/r7+8vf3z9/hQIAAI8wReC5cOGCJKlixYq5tsl8LrNtbqxWq1PHnDRpkiZPnuxcgQAAwKNMEXiSk5MlSX5+frm2yeyNSUpKynNfcXFxCg4Ovu4x6d0BAKD4MEXgCQgIkCRdvXo11zYpKSmSpMDAwDz3FRwc7FTgAQAAxYcpBi07c7nKmcteAADAnEwReOrXry9JOnbsmFJTUx22iY2NzdYWAACUHKYIPC1btpSvr6+Sk5MVHR2d4/lr164pKipKktSuXbuiLg8AAHiYKQJPcHCwunfvLkmaPXt2jueXLl0qm82mSpUqqUuXLkVcHQAA8DRTBB5JmjBhgiwWiz755BMtXrzYvn3Xrl16+umnJUnPP/98nndyAQAAc7IYhmF4ughXeeONNzRx4kRJGRMNlitXTnv27FF6err69OmjFStWqFSpUg5fa7PZVL58eV26dIm7tAAAKCac/f42TQ+PlNHLs2rVKnXr1k0JCQk6dOiQmjZtqmnTpuUZdgAAgLmZqoenMOjhAQCg+CmRPTwAAACOEHgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXiQp5SUFE2ePFkpKSmeLgUuxHk1L86tOXFeC4+lJf7E0hKO8bmYE+fVvDi35sR5zR1LSwAAAPyJwAMAAEyPwAMAAEyPwAMAAEyPwAMAAEyPwAMAAEyPwAMAAEyPwAMAAEyvtKcL8BaZ8y/abDYPV+JdMj8PPhdz4byaF+fWnDivucv8TK43jzIzLf/p+PHjslqtni4DAAAUQFxcnGrWrJnr8wSeP6Wnp+vkyZMKCgqSxWLxdDkAAMAJhmEoMTFRYWFh8vHJfaQOgQcAAJgeg5YBAIDpEXgAAIDpEXgAAIDpEXhM5PDhw/r44481evRoNW/eXKVLl5bFYtHrr79eqP3+9NNPGjBggCpXrqzAwEA1atRIU6ZMUXJystP7WLt2rSwWiywWi7p3716oekoabzyva9as0aBBgxQWFiZ/f39Vq1ZNXbp00TvvvFOomkoSbzqvV69e1fTp09W+fXuVL19evr6+ql69uu68806tX7++UPWUNK4+r6dPn9b8+fP1+OOPq23btvL395fFYtGoUaOcev3+/fs1ZMgQVa9eXQEBAapbt66effZZXbx4sUD1FGsGTOPJJ580JOX4mTJlSoH3+emnnxqlSpUyJBk1atQwWrZsafj6+hqSjDZt2hh//PHHdfeRlJRk1KtXz17PrbfeWuB6SiJvOq/p6enG2LFj7TXUrFnTaNOmjVGnTh2jdOnSRqVKlQpcU0njLef1jz/+MDp06GA/fp06dYxWrVoZFSpUsG976623CvNWSxRXn9d3333X4f4eeuih6752/fr1RmBgoCHJqFy5stGqVSujTJkyhiQjPDzcOH36dIFqKq7o4TGR0NBQ9e3bV6+99pq++eYbDRo0qFD7O3LkiB566CGlpaXp7bffVlxcnKKjo3Xw4EHdeOONioqK0vPPP3/d/bz++us6dOiQ+vfvX6h6SipvOq8TJkzQBx98oCZNmmjbtm2Ki4vTtm3bdPjwYSUkJGjOnDmFqq0k8Zbz+q9//Us//fSTKleurJ9//lmHDx/Wjh07dObMGU2ePFmS9PLLL+vQoUOFqq+kcPV5DQ4OVo8ePTRhwgStWLFCTzzxhFOvS0xM1L333qukpCSNGzdOJ06c0I4dO3Ts2DF17NhRsbGxeuihhwpVW7Hj6cQF9xk2bFih/mXx6KOPGpKM2267LcdzmzdvNiQZvr6+ef4rYd++fYafn5/Ru3dvY86cOfTwuICnzusvv/xilCpVyqhcubIRHx9foGMjd546r+3btzckGe+9957D/bZo0cKQZMycObNAdZV0hT2vfzVp0iSnenjefvttQ5Jx0003GampqdmeO3r0qFG6dGlDkrFjxw6X1FUc0MMDhwzD0LJlyyTJ4b8Cbr75ZjVs2FDXrl3TihUrct3HmDFj5OPjoxkzZri1XjinMOd1xowZSktL05NPPqkqVaoUSb1wTmHOa1JSkiQpPDzc4b7r1q0rSUpNTXVlyXCzL7/8UpI0fPhwlSpVKttztWrVso+l/Pzzz4u8Nk8h8MChY8eO6dSpU5Kkjh07OmyTuX3r1q0On589e7Z+/PFHvfTSS7n+ZYqiVZjzumrVKklS3759FR0drccee0w9evTQgAED9Pe//11nzpxxY+XIS2HOa7NmzSRJW7ZsyfGalJQU7dixQ5LUpk0bl9UL90pNTbWft4L+/W1GBB44dPDgQUmSv7+/wsLCHLbJDDGZbbM6e/asXnjhBdWrV08vvPCC+wpFvhT0vJ4+fVonT56UxWJRZGSk2rZtq5kzZ2rt2rVauXKlJkyYoPr162vt2rXufxPIoTD/v7744osqV66c3nnnHf3rX//SiRMnlJSUpJ07d2rQoEE6cuSIHnjgAbVv3969bwIuc+TIEV27dk1S7j13ef39bVYEHjh04cIFSVKFChVyXVusYsWK2dpmNX78eJ0/f14zZsyQv7+/+wpFvhT0vGb2HlgsFj3zzDNq27atoqOjlZKSor1796pHjx6y2WwaNGiQ4uLi3Pwu8FeF+f+1UaNG2rx5s3r06KFnn31WNWvWVJkyZdSyZUv9/PPPev/99zVv3jz3vgG4VNZznHne/yqvv7/NisADhzLn7PDz88u1TWaQyRwDkGndunVauHCh7rrrLvXs2dN9RSLfCnpe//jjD0kZi+yWK1dOX331lVq2bCk/Pz81atRIK1asUFhYmGw2m6ZNm+a+NwCHCvP/q5RxSSw+Pl6GYSgsLEwtWrRQuXLl7Hfe7d692z2Fwy2yzrmU25+JvP48mBWBBw4FBARIypiQLDcpKSmSpMDAQPu25ORkjR07VuXKldO7777r3iKRbwU9r5mvk6ShQ4fm+FdjYGCgxo4dK0n69ttvXVYvnFPQ8ypJCxcuVP/+/XXixAlt2LBBJ06cUExMjBISEjRx4kRFR0frb3/7mw4fPuy+NwCXyvr/a25/JnL782BmBB44lPmFdvHiRRmG4bBNZldo1i+/t956S4cOHdKkSZNUs2ZN9xeKfCnoec36uGHDhg5fd9NNN0nKGD+AolXQ83rt2jU988wzMgxD06ZNU+fOne3P+fn5acqUKbrtttuUmJioqVOnuvEdwJWynuPcLlk5+vNgdgQeOFS/fn1JGf8KOHnypMM2sbGx2dpKUkxMjCTp7bffVrVq1bL9PPnkk5KkH3/80b6N8R5Fq6DntU6dOvYu8NzGZGVuT0tLc1m9cE5Bz+vBgwcVHx8vSbr11lsdvi7z9uXt27e7rF64V506deTr6yvpf+f9rxz9eTA7Ag8cqlWrlqpVqyZJ2rx5s8M2mdvbtWuX47mzZ88qPj4+24/NZpOU0cWauY0vx6JV0PNaqlQp+23J1/sLtEaNGi6rF84p6HlNTEy87r4ze4zys3YePKt06dJq1aqVpIL9/W1WBB44ZLFYdOedd0rKmE/nr7Zs2aIDBw7I19c325IRy5cvl2EYDn8ylx249dZb7dvq1KlTJO8HGQp6XiXpnnvukSQtXrzYfstrVpl38nTr1s3VZeM6Cnpe69ata7+ra926dQ73nTnVQIMGDVxdNtxo4MCBkqS5c+fm+IflsWPH7Oe1sEtfFCcEnhJu2rRpqlOnju67774czz333HPy8/PT999/r3feecf+L72jR49q5MiRkqRRo0bZ/2UJ7+GO8zpq1ChZrVYdOXJETz75pH0wZFpamiZMmKCYmBj5+flp/Pjxbn53JZerz2toaKj9TsqnnnpKP/zwg/25q1ev6v/+7/+0Zs0aSdKDDz7otvdV0uV1Xgtq7NixCg0N1f79+/X000/b/5GSkJCgwYMHKzU1Vb1791ZERITLjun1inQhC7jVpk2bjEqVKtl//P39DUlGmTJlsm0/duyY/TWZ67J07tzZ4T7nzZtn+Pj4OFx9OSIiwrh8+bLT9bGWVsF403mNiooygoODDUlGxYoVjTZt2hiVK1c2JBmlSpUy5s6d646PwJS85bweOXLEqFWrln0V7ho1ahgtWrQwgoKC7NtGjx7tro/BdFx9Xo8dO5btdZmrn/v7+2fbvmnTphyvXbt2rREQEGBfLT0iIsK+WnqdOnWMU6dOufOj8Dr08JjItWvXlJCQYP/JvO3wypUr2bbnZ9zM0KFD9eOPP6pv375KSkrSvn37FB4ersmTJ2vTpk0qW7asu94O/uRN57V169bavXu3Ro0apbJly2rnzp2SMrrPt2zZomHDhhX6/ZYU3nJea9eurV27dmnSpElq2bKlLl26pD179iggIEC9e/fWF198oY8++shl79vsXH1e09LSsr0uc96clJSUbNsdXWa+9dZbtX37dt13332yWCz65ZdfVLVqVT399NOKjo4ucb3zFsPI5R5GAAAAk6CHBwAAmB6BBwAAmB6BBwAAmB6BBwAAmB6BBwAAmB6BBwAAmB6BBwAAmB6BBwAAmB6BBwAAmB6BB4BHzZ07VxaLRRaLRUeOHPGafXnyGABcj8ADIFcbNmywf7lPnjzZ0+UAQIEReAB4NXpUALgCgQeARw0fPlyGYcgwDNWpU8dr9lWULl++rBMnTuj48eM6ceKE4uPjdfHixXytlA4gb6U9XQAAlDTp6en673//q//+97/asGGDLl265LCdj4+PKleurAYNGqh169a66667dPPNNxdxtYA50MMDAEVo06ZNatOmjebNm6fmzZtr3Lhxqlixov35xo0bKz4+XikpKTpx4oQWLFigsLAwvfvuu+rYsaOGDh2q9PR0D74DoHgi8AAosMmTJ9vH10hScnKy3nnnHbVq1UpBQUEKCgpS27ZtNWPGDKWmpjrcR25jdDIHTI8YMcK+7YYbbrC3zfzZsGHDdfeV1Z49e/T666+rZ8+eqlmzpvz9/VWuXDnVr19fw4YN088//1zozyU3H374oSZPnqwvvvhC3377rV599VW99tpr6t69u73NmDFjVKVKFfn5+alatWrq0aOHlixZYv8cFixYoP/85z9uqxEwKy5pAXCJ+Ph49ezZU7t27cq2PSoqSlFRUfr++++1fPly+fh47t9ZGzZsUNeuXXNsv3r1qg4dOqRDhw5p/vz5evHFF/Xmm2+6/PjNmzfXqFGjVKpUqWzb9+zZY3/cr18/h68dOXKk5syZI0n6/vvvNWrUKJfXB5gZgQeASwwcOFD79+/XuHHj1K9fP4WEhOjXX3/VlClTtH//fq1atUoff/yxxowZ49T+2rRpo19++UUrVqzQxIkTJUnfffedwsLCsrW74YYbnK4xNTVVZcuWVZ8+fdStWzc1bNhQwcHBOnPmjPbu3av33ntPR48e1dSpU9WgQYNsvUuu0L59+xzbTp8+rf3790uSbrrpplwHW1epUsX+OCkpyaV1ASUBgQeAS2T24nTp0sW+rVWrVurZs6caNWqk+Ph4zZw50+nAU7ZsWTVp0kTbt2+3b2vQoEGh7r5q0aKFjh8/rgoVKuR4rmfPnnr88cfVt29frVmzRq+++qqGDh2aozfG1b777jv749tvvz3XdmfOnLE/zk/IA5CBMTwAXOKJJ57IFnYyhYSE2HtKdu/enesdSUUhNDTUYdjJ5Ofnp3feeUeSdPToUe3cudPtNX311Vf2x3kFnsjISPvjTp06ubUmwIzo4QHgEkOGDMn1uYiICPvjw4cPq0WLFkVQ0fWlpKQoPj5ely9ftt/5ZBiG/fldu3Zlq93Vrl27Zu/hCQoKyjXIJCYmatasWZIyLnvdcccdbqsJMCsCDwCXaNiwYa7PhYSE2B8nJiYWRTm5+uOPP/Tee+9pyZIl2rt3b56T+507d86ttWzcuFE2m02S1KNHD/n6+uZok56erpEjR+rUqVMKDAzU/PnzHbYDkDcCDwCXKFOmTK7PZb0zy5OzBx85ckTdunXT4cOHnWrv7sHBq1atsj92dDkrPj5ew4cP17fffquQkBB9/vnnat26tVtrAsyKMTwASowHH3xQhw8flsVi0ciRI/X9998rLi5OycnJ9iUpsgayrJe33GH16tX2x71797Yfc9u2bXr66adVr149/fDDDxozZoz27t3r8JZ6AM6hhwdAiXDgwAFt2rRJkvTSSy/pjTfecNjuwoULRVLPvn37FBsbKynj7rGAgAD17dtXp06dUrly5dSiRQt98skn6tu3r8qWLVskNQFmRuAB4NUyZ3EurL1799of33fffbm2y3obvDv99XJWSEhIth4fAK7FJS0AXi0gIMD+OCUlpcD7ybq0xZUrV3Jt98EHHxT4GPlxvfE7AFyLwAPAq1WvXt3++Pfffy/wfurXr29/PG/ePIdtZs2apeXLlxf4GM46d+6cfc2uSpUqOZyBGYBrcUkLgFdr2bKlAgIClJycrP/7v/9T6dKlVadOHfudXzVq1FBgYKBT+2nSpIn27NmjWbNm6eLFixoyZIiqV6+uuLg4ffrpp/r888/VsWNHbd682a3v6bPPPrMPjn700UfdPpszAAIPAC8XFBSkcePG6e2331Z0dLR69uyZ7fnIyEiHMzz/lcVi0YIFC9StWzdduHBBixcv1uLFi7O1adq0qZYuXZpjvS5XunDhgv75z39Kypgx+bnnnnPbsQD8D5e0AHi9qVOn6uOPP1anTp0UEhJS4B6RFi1aaOfOnRo7dqxq164tX19fhYSEqG3btvrHP/6hbdu2ZbuE5mqbNm1Sp06dVK1aNf373//W999/r6CgILcdD8D/WAx3TzQBAADgYfTwAAAA0yPwAAAA0yPwAAAA0yPwAAAA0yPwAAAA0yPwAAAA0yPwAAAA0yPwAAAA0yPwAAAA0yPwAAAA0yPwAAAA0yPwAAAA0yPwAAAA0yPwAAAA0yPwAAAA0/t/w5lUeaBBK9oAAAAASUVORK5CYII=",
+      "text/plain": [
+       "<Figure size 640x480 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "# plot\n",
+    "plt.scatter(P_ratios, log_t_insts, c='silver', s=10)\n",
+    "plt.tick_params(labelsize=16,size=8,direction='in')\n",
+    "plt.xlabel(r'Initial $\\mathcal{P}$', fontsize=20)\n",
+    "plt.ylabel(r'$\\log_{10}(t_\\mathrm{inst}/P_1)$', fontsize=20)\n",
+    "plt.title(r'Model (2)', fontsize=20)\n",
+    "plt.xlim([1.028, 1.112])\n",
+    "plt.ylim([-0.4, 4.5])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "de35da87",
+   "metadata": {},
+   "source": [
+    "Model (2) already appears to do a better job at predicting instability times. To determine exactly how accurate the predictions of the model are, one can compare with true instability times from full $N$-body simulations (see Figs. 4-9 in the paper)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f11e3faa",
+   "metadata": {},
+   "source": [
+    "# Model (2, 1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "id": "792ed7eb",
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "100%|████████████████████████████████████████████████████████████████████████████████████| 25/25 [07:04<00:00, 16.99s/it]\n"
+     ]
+    }
+   ],
+   "source": [
+    "P_ratios = []\n",
+    "log_t_insts = []\n",
+    "\n",
+    "# predict instability times for 25 systems\n",
+    "for _ in tqdm(range(25)):\n",
+    "    P_ratio = np.random.uniform(1.03, 1.11)\n",
+    "    sim = initialize_sim(P_ratio)\n",
+    "    \n",
+    "    # get indices of bordering first-order MMRs\n",
+    "    j1 = int(np.floor(1/(1 - 1/P_ratio)))\n",
+    "    j2 = int(np.ceil(1/(1 - 1/P_ratio)))\n",
+    "    j3 = j1 + j2 #second-order MMR that lies between the first-order MMRs\n",
+    "    \n",
+    "    res_inds1 = [[j1, 1, i, i+1] for i in range(1, 5)] # all j1 resonances\n",
+    "    res_inds2 = [[j2, 1, i, i+1] for i in range(1, 5)] # all j2 resonances\n",
+    "    res_inds3 = [[j3, 2, i, i+1] for i in range(1, 5)] # all j3 resonances\n",
+    "    \n",
+    "    pham = construct_Hamiltonian(sim, res_inds1 + res_inds2 + res_inds3)\n",
+    "    t_inst = predict_t_inst(pham)\n",
+    "    \n",
+    "    P_ratios.append(P_ratio)\n",
+    "    log_t_insts.append(np.log10(t_inst))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "id": "157c720e",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle \\left( -0.4, \\  4.5\\right)$"
+      ],
+      "text/plain": [
+       "(-0.4, 4.5)"
+      ]
+     },
+     "execution_count": 10,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAHhCAYAAABqR2T2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABJPklEQVR4nO3deVxVdf7H8fdFVhXcUARFDbXMDcQ9c1zKzNwqm6m03EbT9nGaZiotMGvGan5TNo62jqm5TLZo2lhpqWWWIqDmlhouqIiKy8W8oMD5/UHcAbngBe4Ch9fz8eDxuJz1c+9R79vv+Z7v12IYhiEAAAAT8/F2AQAAAO5G4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AHgMfHx8bJYLLJYLG47R9++fWWxWNS3b1+XHG/cuHGyWCx68MEHXXK8yu7ll1926ecHVBYEHqCSW79+vT0kWCwWBQcH6+LFi1fdz2azqU6dOkX2Xb9+vfsLNpGtW7dq/vz58vf319NPP+1wm59++kmvvvqqbr/9dl1zzTUKCgpSzZo1dc011+juu+/WZ599Jk9NWZiTk6Pk5GS9+eabmjBhgjp27ChfX1/79T906NBVj/Hwww8rNDRUGzZs0EcffeT+ogEP8fV2AQDK5sKFC1q+fLlGjhxZ6nYrVqyQ1Wr1UFXmNHXqVBmGoXHjxqlZs2bF1o8ZM0YLFixwuO+hQ4d06NAhffDBBxo4cKCWLl2qunXrurXeF198UfHx8RU6Rq1atTRlyhRNnTpVzz77rO644w75+PB/Y1R9/CkGqpDAwEBJ0sKFC6+6bcE2BfugbBISEvTll19Kkp544gmH2xw7dkySVL9+fT3wwANavHixNm3apC1btujNN9/UddddJ0n64osvNHToUOXl5bm15sItSYGBgerRo4datmxZ5uM8/PDDCgwM1J49e2jlgWkQeIAqZNiwYZKkNWvW6MSJEyVud/LkSfuX9fDhwz1Sm9m89tprkqRu3bqpdevWDrdp2rSp3nzzTR0/flxvvvmm7r33XvXs2VNdu3bVAw88oKSkJN14442SpI0bN2rRokVurblnz5564403lJiYqMzMTH3//ff285dFnTp1dNttt0mSZs2a5eoyAa8g8ABVyC233KLGjRsrNzdXS5YsKXG7JUuWKCcnR2FhYRowYIAHKzSH8+fP21s27rvvvhK3e++99/TAAw8oICDA4fqaNWtq7ty59t8//PBD1xZ6hYEDB2rSpEmKjY2Vr2/FeiyMGjVKkvTdd9/pp59+ckV5gFcReIAqpEaNGrr33nsllX5bq6BfyciRI1WjRg2njn3p0iXNmTNH/fr1U8OGDeXv76/GjRvrtttu0/vvv+/U7ZijR4/q4YcfVlRUlAIDAxUREaFhw4Zp7dq1TtVQ4OLFi3rttdfUr18/hYWFyd/fX40aNdItt9yiefPmKTc3t0zHK6sVK1YoOztbkjRixIgKHat9+/YKDQ2VJP38888Vrs1TBg8ebL8d+p///MfL1QAVR+ABqpj7779fkpScnKxdu3YVW797924lJSUV2fZqDh8+rJiYGD388MNav369Tp8+rcuXLys9PV2rV6/W/fffrz59+ujMmTMlHmPDhg1q166d5syZo4MHDyo7O1tpaWlauXKlBgwYoOeff96pWhISEnTttddqypQpWr9+vU6ePKnLly/r1KlTWrNmjcaPH68bbrhB6enpTh2vPNatWydJioyMVERERIWPd+nSJUmqUp1/AwICFBMTI0k83QdTqDp/+wBIkjp16qT27dtLctzKU7CsXbt26tSp01WPd+HCBfXv31979uyRJN1+++369NNPtXXrVi1btkx9+vSRlN8HZciQIQ5bVw4dOqShQ4fKarXKx8dHkydP1tq1a5WQkKB3331XrVu3VlxcnD777LNSa/nxxx/Vr18/HTt2TI0aNVJcXJzWrl2r5ORkffHFF3r44Yfl6+urLVu2aPjw4bp8+fJV3195fPvtt5Kkrl27VvhYycnJ9qfl2rRpU+HjeVK3bt0kST/88IPbPmvAYwwAldq6desMSYYkY968eYZhGMZLL71kSDKaNm1q5Obm2rfNy8szIiMjDUnGzJkzDcMwjHnz5tn3X7duXbHj/+lPf7KvnzZtWrH1eXl5xqhRo+zbzJkzp9g2d955p3394sWLi623Wq1GdHS0fRtH//Tk5eUZHTt2NCQZ0dHRxqlTpxx+HqtXrzZ8fHwMScY777xTbH2fPn0MSUafPn0c7n816enp9hpnzJhRrmMUdtddd9mPt2zZsgofr6zGjBljP//BgwfLtO/8+fPt+27ZssU9BQIeQgsPUAWNGjVKPj4+Onr0qDZs2GBfvn79eqWmpsrHx8fe6bQ02dnZeueddyRJbdu2dTiGi8Vi0Zw5c9SgQQNJ0uzZs4usT0tL04oVKyRJQ4YMsfcxKiw4OFhvvfVWqbV89tln2rFjh6T8PkgF/V6udOutt+quu+6SJM2bN6/UY5bH0aNH7a8bNWpUoWN99NFH9o7KnTt3rnB/IE8r/P4Lfy5AVUTgAaqgJk2aqF+/fpKK3tYqeN23b181bdr0qsdJTEzUuXPnJEljx44tsYNzSEiIfve730nK7yOUlpZmX7du3Tr7ba5x48aVeK5u3bqpXbt2Ja4vCE3XXXedOnbsWGrdv/nNbyTl9/dxdQfmU6dO2V/Xq1ev3MfZu3ev/fMICgrSggUL3DqlhjvUr1/f/rrw5wJURQQeoIoaPXq0pPxHnW02m2w2m/1Ramc7K+/cudP+unv37qVuW3h94f1+/PFH++ur9Xkp6BPiyNatWyXlT9VQeDoMRz+PPPKIpPzOwKV1pC6Pwscrb+A5fvy4Bg0apMzMTFksFr377rtq27atq0r0mMLvPyMjw4uVABVH4AGqqDvvvFM1a9ZUZmamVqxYoeXLl8tqtSooKMjpWyeFv9zDwsJK3bZx48YO9zt79qz99dVuAZV2jpMnT5a6b0mcmVesLAqPTG2z2cq8/5kzZ3TLLbfY562aNWuWw9t8VUHh9x8UFOTFSoCKYy4toIqqXbu27rjjDi1atEgLFy60Tytw++23Kzg4uMzHu9rtFqOECTALLy/vMSTZb0316tVLb7zxRqnHKcwVj40X1rBhQ/vrsrYeZWZm6tZbb7UPFzBjxgw9+uijLq3Pkwq//8KfC1AVEXiAKmz06NFatGiRfRoJyfnbWVLRPhonTpzQtddeW+K2hce9Kbxf4dfp6emKjIws8RilteI0aNBA6enpOnXqlP2xe28o/MVeuPXqamw2m4YOHaqEhARJ0pNPPqlp06a5vD5PKvz+CTyo6rilBVRhN910k8LDw5WTk2OfSuKWW25xev/CwWLz5s2lbrtlyxaH+3Xo0MH+uuDLviSlrS8YM2jfvn06fPhwqcdxp5YtW6pmzZr2Wpxx+fJljRgxwv7E3OTJk/Xyyy+7rUZPKfz+C19noCoi8ABVWI0aNXT//fcrICBAAQEBuu+++5yeSkLKf1S6bt26kqT58+eX+MRTZmamPvjgA0n5j6+Hh4fb1/Xr189+zvnz55d4rq1btxbp7HylgolRJXk1LPj6+qpHjx6Srh7gpPxbcSNHjtTq1asl5bewzZkzx601ekrB+2/ZsmWRaw5URQQeoIp76aWXlJWVpaysLP39738v074BAQGaMGGCJGnXrl2aPn16sW0Mw9Ajjzyi06dPS5L9CakC4eHh9hnZP/30U3swKuzChQt64IEHSq1lxIgRuv766yVJc+fO1bvvvlvq9jt37tTKlStL3aa8evfuLUnasWOHfU4tRwzD0MSJE+1j7YwYMULz5s0r1+PnBU+gtWjRolw1u0NBq17B5wFUZfThAaq55557Th9//LFSUlI0Y8YM7dy5U+PHj1dERIQOHjyo2bNn2+dS6tmzp8Pg8n//939as2aNMjMzNXLkSG3YsEF33XWXQkJCtGPHDs2cOVP79u1Tly5d7I+fX6lGjRr6z3/+oxtuuEEXLlzQhAkTtGzZMo0cOVLXXXed/Pz8dPLkSSUnJ2vVqlXatGmTnnjiCQ0dOtTln8ngwYM1ffp0Xbp0Sd9++61uvvlmh9v96U9/sg9+2L59ez3zzDP2KTpK4s7+SRcuXCg2I/uBAwfsrz/88MMiAzrGxMTY58u60v79+5Wamiop//MAqjxvDvMM4OocTS1RFlebWsIwDOPgwYNGmzZtikz9cOVPr169jIyMjFLrDA4OLnH/uLg4Iy4ursSpJQps377daN26dam1FPxMnz692P4VnVqiQMHnMW7cuBK3ad68uVN1Fv5x5OLFi/b1sbGx5a754MGDZaolLi6uxGPFx8cbkow6deoYNput3DUBlQW3tACoRYsW2r59u2bPnq0+ffqoQYMG8vPzU1hYmG699VYtXLhQ33zzTZEnsq7Ut29f7dq1Sw8++KCaN28uf39/hYWFafDgwfr8888dTlvhSMeOHbV7927Nnz9ft99+uyIjIxUYGCh/f3+Fh4erb9++mjZtmhITE/Xcc8+56BMo7qGHHpKUPz1EVlaW284jSd9//7399ZQpU9x6LmctWbJEkjR+/PgiYxMBVZXFMEoZGAMAqqkLFy6oRYsWysjI0MKFC3Xfffe57Vzx8fGaPn26WrdurT179pSp47k7bNy4Ub1795afn5/27dtXqfoVAeVFCw8AOFC7dm39+c9/liS9+OKLysvLc9u5Ch5nf+aZZ7wedqT8AROl/LnRCDswC1p4AKAEWVlZatOmjQ4fPqwlS5bonnvucfk5Ll26pLp166px48bat2+ffH29+yzJli1b1L17d9WuXVv79u3jcXSYBk9pAUAJAgMD9f7772vt2rXKyclxyzn8/f1dPh9YRZw+fVpxcXGKjY0l7MBUaOEBAACmRx8eAABgetzS+lVeXp6OHz+u4ODgco2SCgAAPM8wDGVmZioiIkI+PiW34xB4fnX8+PFSZ3kGAACVV2pqqpo2bVriegLPr4KDgyXlf2AhISFergYAADjDarUqMjLS/j1eEgLPrwpuY4WEhBB4AACoYq7WHYVOywAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPRMHXimTZsmi8Uii8WiF154wdvlAAAALzFt4NmzZ49eeeUVb5cBAAAqAVMGHsMwNGnSJPn5+al///7eLgcAAHiZKQPPu+++q2+//VbPPfecIiMjvV0OAADwMtMFnlOnTukvf/mL2rZtqylTpni7HAAAUAn4ersAV5syZYrOnDmjjz/+WH5+ft4uBwBQxVitVtlsNgUFBSkkJMRtx8/Ly5OPj4/D85RUg7trc0ZaWpoyMzMVHBys8PDwYusrQ42OmCrwfPXVV1q0aJHuu+8+9enTx9vlAACqmJSUFKWmptp/j4yMVFRUlNuO7+g8JdXg7tqckZSUpMzMTEn5wSctLU2xsbH29ZWhxpKYJvBkZWVp8uTJqlOnjv7+97+X+zhWq9Wp7QICAhQQEFDu8wAAKher1VosjKSmpio0NNQlLRWOjn/leQpeX7kuKCjIrbU5o6Blp7DMzEylpaUpPDzc7Z9fRZmmD88LL7ygAwcO6MUXX1RYWFi5jxMZGak6depc9edvf/ubC6sHAHibzWYr03JXHb/w+pK2uTJoOHtMVyqphoLl7v78KsoULTwFY+7ExsbqwQcfrNCxUlNTnUqitO4AgLmcOXPG4fKgoCCXHP9qxyltfXBwsNLS0sp8TFcqqYbg4OBSa/FkjaUxReB56KGHlJOTo7lz58rHp2KNViEhIZWi6Q0A4DlWq1UnT54strxRo0Yu+04ICQlRZGRkiX14Cs5z5TaRkZEKDw+XzWYrttyT31fh4eHFbmsV7rjs6P15usbSWAzDMLxdREXVrVtXFy5csN//LOz8+fPKyspS7dq1VatWLUVGRiohIaHYdlarVXXq1NH58+crzcUBAHhGenq69u7dW2x5mzZtKtRNwhGe0nItZ7+/TdHCI0m5ublKT08vcf2FCxd04cIFBQYGerAqAEBV4MnbMc7cSShpm8pwFyI8PNxh0ClQGWp0xBSdls+dOyfDMBz+jBkzRpI0Y8YMGYahQ4cOebdYAEClU3A7prDKdDsGFWeaFh4AACoiKipKoaGhXr9lBPcg8AAA8KvKejsGFWeKW1oAAAClMcVTWq7AU1oAAFQ9zn5/08IDAABMj8ADAABMj8ADAABMj8ADAABMj8fSAQAuVxmmQAAKI/AAAFwqJSWl2ASSUVFRXqwI4JYWAMCFrFZrsdnAU1NTZbVavVQRkI/AAwBwGZvNVqblgKcQeAAALuPJWceBsiDwAABchlnHUVnRaRkAqil3PUnFrOOojAg8AFANuftJKmYdR2XDLS0AqGZ4kgrVEYEHAKoZnqRCdUTgAYBqhiepUB0ReACgmuFJKlRHdFoGgGqIJ6lQ3RB4AKCa4kkqVCfc0gIAAKZH4AEAAKZH4AEAAKZHHx4AgFe5a4oLoDACDwDAa9w9xQVQgFtaAACvYIoLeBKBBwDgFUxxAU8i8AAAvIIpLuBJBB4AgFcwxQU8iU7LAACvYYoLeAqBBwDgVUxxAU/glhYAADA9Ag8AADA9Ag8AADA9Ag8AADA9Ag8AADA9Ag8AADA9Ag8AADA9Ag8AADA9Ag8AADA9Ag8AADA9Ag8AADA9Ag8AADA9Ag8AADA9Ag8AADA9Ag8AADA9Ag8AADA9Ag8AADA9Ag8AADA9Ag8AADA9Ag8AADA9X1cf8PTp0zp27JhOnTqljIwMBQUFqWHDhmrYsKGioqLk40PGAuB+VqtVNptNQUFBCgkJ8XY5ALyswoEnMzNTK1as0Pr16/Xtt9/qwIEDJW5bq1Yt9ejRQ71799bgwYMVGxtb0dMDQDEpKSlKTU21/x4ZGamoqCgvVgTA2yyGYRjl2TExMVGzZs3SRx99pKysLEmSM4eyWCz219ddd50efvhhjR07VrVq1SpPGS5jtVpVp04dnT9/nv8NAlWY1WpVcnJyseWdOnXi7zZgQs5+f5c58CQmJmratGn68ssvJf0v5ISHh6tr167q3LmzGjVqpPr166tevXqy2Ww6c+aMzp49q3379ikhIUE7duzQ5cuX8wuwWFS/fn09+eSTevzxxxUQEFDe91whBB7AHNLT07V3795iy9u0aaOwsDAvVATAnZz9/i7TLa1x48Zp4cKFysvLkyTFxsZq1KhRGjFihJo1a+b0cS5duqRvvvlGixcv1ieffKKMjAw9/fTTmjt3rhYuXKgbb7yxLGUBgF1QUFCZlgOoHsrUg3j+/Pny9fXVxIkTtXfvXm3dulVTpkwpU9iRJH9/f918883697//rRMnTmjBggW67rrrdPjwYX399ddlOhYAFBYSEqLIyMgiyyIjI2m5Baq5Mt3Sevjhh/X000+radOmLi/EMAwtW7ZMubm5uvfee11+/KvhlhZgLjylBVQPbuvDY1YEHgAAqh5nv79NMyjO8uXLNWnSJHXu3Fnh4eHy9/dX3bp1dcMNN2jWrFm6dOmSt0sEAABeYpoWnhtvvFHfffedAgICFBERoQYNGigtLU3Hjh2TJHXu3Flr165V3bp1He5PCw+AiuAWGuAd1a6FZ8KECVq3bp0yMzOVkpKihIQEHT16VN9//72aNm2qxMRETZ061dtlAjChlJQUJScna+/evUpOTlZKSoq3SwJwhUoVeFasWKHnn3++XPuOHTtWffv2lZ+fX5HlPXr00D/+8Q9J+be9AMCVrFZrkVGdJSk1NVVWq9VLFQFwpFIFnuXLl2v69OkuP26bNm0kSRcvXnT5sQFUbzabrUzL3cFqtSo9PZ2QBZTC5ZOHVkbff/+9JDF3FwCX8/ZAh8wbBjjHrYFnwYIFZdq+tIlHyyo3N1dpaWn69NNP9dRTT6lWrVr629/+5rLjA4D0v4EOrwwdnui4XNLttNDQUDpOA1dwa+AZO3ZskclCr8YwjDJt78hrr72mKVOmFFl2++23a8aMGWrfvv1V93e2STggIMBr834BqFyioqIUGhrq8ae0SrudRuABinJr4PH391dERIQmTZrk1PbLli1zOMtxWTRp0kS9evXS5cuXdfjwYaWnp2vdunVasmSJnn/+edWoUaPU/a8ckr4kcXFxio+Pr1CtAMwjJCTE4yHD27fTgKrErePwdO3aVUeOHFF6erpT248bN04LFixQbm6uy2rYvHmzJk2apO3bt2vy5MmaO3euw+0KnuNPTU116h8tWngAVAb04UF155bZ0suqc+fOSkpKUmpqqtMtJ67WvXt3/fe//1VUVJTeeustPfXUU2revHmJ23vjf2kAUF7eup0GVDVufSy9d+/eatasmfbv3+/U9jfeeKNGjx7t8joiIiIUExOjvLw8bd++3eXHBwBvCgkJUVhYGGEHKIVpppa4mi5duigxMVEfffSR7rzzzmLrmVoCAICqx+23tH7++WctXrxYe/fuVW5ursLDw/Wb3/xGt912W6Xr23Lo0CF7y050dLSXqwEAAJ5WrsDz5ptv6rHHHlNOTk6R5a+//roiIyP1xhtv6NZbb3VJgc5ITEzUp59+qjFjxhTrrPf5559rypQpysnJ0W233aaWLVt6rC4AAFA5lPmWVkJCgnr27Km8vLwSt/H19dWHH36oYcOGVbhAZ6xfv179+vWTJDVu3FhNmzbVpUuXdOTIEZ07d05S/hNj//3vfxUaGurwGNzSAgCg6nHbbOn//Oc/lZeXJ4vFoiFDhuiTTz7R5s2b9fHHH2v8+PHy8/NTTk6Oxo8fr9OnT1foTTgrOjpas2bN0rBhw1SrVi3t3btXe/fuVVBQkAYNGqR58+Zp06ZNJYYdAABgbmVu4YmKitLhw4c1aNAgrVq1qtj6b775RgMHDtSlS5c0Y8YMPfPMMy4r1p1o4QEAoOpxWwtPWlqaJJU4evJvfvMbPfHEEzIMQx999FFZDw8AAOByZQ482dnZkqQWLVqUuM29994rSfrxxx916dKl8lUGAADgIuUeeLC0OalatWolKX/G8lOnTpX3FAAAAC7hlpGWC4/Dk5mZ6Y5TAAAAOK3cgcdisTi1XWmPrwMAAHhCuUdavvHGGxUdHa2YmBj7T9u2beXr69b5SAEAAMqsXOnEMAydPXtWGzZs0IYNG+zL/fz81LZtW8XExNiXXb58ucJFAgAAVESZx+F54403tG3bNm3btk07d+7UxYsXix7wiltdPj4+atmypTp27KiOHTsqOjpaHTt2VPPmzStevQsxDg8AAFWPs9/fFZot3TAM/fTTT/YAtG3bNm3fvl3p6enFT3RFEAoODlaHDh0UHR2t2bNnl7cElyHwAABQ9Xgk8JQkPT1dycnJRYLQgQMHHHZgtlgsys3NdXUJZUbgAQCg6vFq4HHEZrNp+/btRULQzp07ZbPZCDwAAKBcnP3+9tgjVUFBQerRo4d69OhhX2YYhvbt2+epEgAAQDVV5nF4EhMTXXZyi8Wi6667zmXHAwAAcKTMgadr165q2rSpJk2apFWrVikrK8sddQEAALhMmfvw+PjkZ6SCp64CAwPVv39/DR06VEOGDFFERITrq/QA+vAA7mW1WmWz2RQUFMTfMQAu47ZOy8ePH9eqVau0cuVKff3117LZbPkH+jUAxcTEaOjQoRo6dKg6d+5cgbfgWQQewH1SUlKUmppq/z0yMlJRUVFerAiAWXjkKS2bzaa1a9dq1apV+uyzz3T8+PH8g/4afho3bqzBgwdr6NChuvnmmxUUFFTeU7kdgQdwD6vVquTk5GLLO3XqxN81ABXmlcfSExMTtXLlSq1atUpJSUn5J6git74IPIB7pKena+/evcWWt2nTRmFhYV6oCICZeH0cnqp264vAA7gHLTwA3MnrgaewrKwsrV27VitXrizx1teQIUP00EMPKTo62t3lOETgAdyHPjwA3KVSBZ4rJSYm2lt/kpOTZRiGLBaL4uLi9Nxzz3m6HEkEHsDdeEoLgDtUupGWC+vcubM6d+6suLg4+62vVatWqWbNmt4oB4AHhISEEHQAeI3HWniOHDmiunXrlvoPntVq1blz59SsWTNPlFTs3LTwAABQtTj7/V3mkZbL65prrtGsWbNK3Wbu3Lm65pprPFQRAACoLjwWeAzD0NUak7zQnQgAAFQDHgs8zjh69KiCg4O9XQYAADAZt3Zafv7554v8vn79eofb5ebm6ujRo1q6dKm6d+/uzpIAAEA15NZOywUTjUr5Y+5c7VQRERH65JNP1LVrV3eVVCI6LQMAUPVUisfS161bJym/b07//v01duxYjRkzpth2NWrUUP369dWmTZsiIQkAAMAV3Bp4+vTpY38dFxenvn37FlkGAADgCV4ZabkwwzB04MABBQUFqWnTpl6rg1taAABUPZVuHJ5PP/1U48eP19mzZ+3LDh06pA4dOqhNmzZq3ry5Ro0apby8PE+VBAAAqgmPBZ65c+cqISFB9erVsy/7wx/+oN27d6tfv37q2LGjli5dqnnz5nmqJAAAUE14LPDs2rVL3bp1s/9+/vx5/fe//9Xdd9+ttWvXasuWLbr++uv17rvveqokAABQTXgs8Jw6dUrh4eH23zdu3KicnBzde++9kiQ/Pz8NGDBABw4c8FRJAACgmvBY4AkJCVFGRob99/Xr18vHx0e9e/e2L/Pz89Mvv/ziqZIAAEA14bHA06ZNG61cuVJnzpzR+fPntXTpUsXGxhbp03P48GGFhYV5qiQAAFBNeCzwPPbYYzp+/LiaNGmiyMhIHT9+XJMnT7avz83N1caNGxUdHe2pkgAAQDXh1oEHCxsxYoT+9a9/2Tsl/+53v9P48ePt67/66itdvHhRt956q6dKAgAA1YTXBx6sLBh4EACAqqfSDTwIAADgLQQeAABgeh7rwyNJly5d0vLly5WQkKBz584pNze32DYWi4XBBwEAgEt5LPAcPnxYAwYM0M8//6zSug0ReAAAgKt5LPBMmTJFBw4c0P3336/x48eradOm8vX1aAMTAACopjyWOL7++mvddNNNmj9/vqdOCQAAIMmDnZbz8vLUqVMnT50OAADAzmOBp2fPntqzZ4+nTgcAAGDnscAzc+ZMrVu3Th9++KGnTgkAACDJg314Vq5cqX79+unuu+9Wnz591KlTJ9WpU6fYdhaLRc8++6ynygIAANWAx6aW8PFxrjHJYrE4HJ/H3ZhaAgCAqsfZ72+PtfCsW7fOU6cCAAAowmOBp0+fPp46FQAAQBHMpQUAAEzPbS0848ePl8Vi0V//+leFhYVp/PjxTu3H1BIAAMDV3NZp2cfHRxaLRXv27NG1115Lp2UAAOByXu+0fPDgQUlSkyZNivwOAADgaS4JPCNHjlR0dLQefPBBe7pq3rx5kW2u/B0AAMBTXNJpeenSpXrmmWd09OhRVxyuzAzD0MaNG/Xkk0+qR48eqlu3rvz9/RUREaERI0bwSDwAANWc157SunDhgrZt26acnJwKH+vrr79W79699fe//10JCQkKCwtT+/btlZmZqY8//lj9+/dn9GYAAKoxrwWew4cPKzY2VsHBwRU+lmEYatWqlebMmaPTp0/rp59+UlJSkjIyMvT0009Lkl544QWtWrWqwucCAABVT7kCz8svv6x169bJarVWuIBLly5V+BjdunXTnj179OCDD6pevXr25f7+/vrrX/+qQYMGSZLefvvtCp8LAABUPeXqtPzUU0/JYrHIYrGoVatW9uXr1q1TYGCgoqKirnqMrKwsSfmPoVfU1R4jHzBggFavXq19+/ZV+FwAAKDqKVfg8fX1VU5OjgzDKBIiHnvsMUlScHCwOnbsqJiYGMXExCg6OlodOnSQv7+/fdsffvhBklS3bt0KlO+cgnAVFBTk9nMBAIDKp1yBJzMzU0lJSdqyZYsSEhK0ePFiSfl9aaT8QYA2btyo7777zr5PjRo1dN1116ldu3bKycnRqlWrZLFY1KVLFxe8jZIZhqFly5ZJknr16uXWcwEAgMrJJSMtF4yqvGzZMtlsNm3btk3btm3T9u3bdfr06aIn/PUWlmEY8vHx0erVqzVgwICKllCit956S5MmTZK/v792796tli1bOtyuYKTG1NRUp0ZaDggIUEBAgKvLBQAAZeDRkZYbNGigM2fOqE2bNmrbtq1GjRplX3f8+PEiAeinn36S1WpVVFSUpkyZ4tawk5SUpMcff1xS/lNaJYWdwiIjI506dlxcnOLj4ytSHgAA8BCXzaWVlpamevXqKTAw0BWHq7CDBw+qV69eSktL08iRI/X++++X2kGaFh4AAKoej8+lFR4e7qpDVdiJEyc0YMAApaWlafDgwXrvvfecfhosJCSEyUMBADAZrw086C5nzpzRgAED9PPPP6tPnz5atmyZ/Pz8vF0WAADwIlMFngsXLui2227Tzp071bVrV61cuZJH0QEAQNkCT8Hj3e5y9OhRbdq0qVz7Zmdna/jw4dq8ebPatWunzz//3CXTVgAAgKqvTIHn7rvvVocOHVwefI4cOaIHH3xQrVq10tq1a8u8f25uru655x59/fXXatmypdasWaP69eu7tEYAAFB1lanTcuvWrbVr1y7dc889evLJJzVy5EiNGjVK7dq1K/OJf/nlF33yySdatGiRvvrqK+Xm5srHx6fIVBXO+uCDD7R8+XJJ+WMC/fa3v3W4XXh4uNtbqQAAQOVTpsfSc3Jy9Prrr2vmzJk6ffq0/cmn1q1bq0ePHuratas6deqkRo0aqV69eqpXr55sNpvOnDmjs2fPat++fUpISNCWLVu0ZcsWZWVl2UdnHjRokF566SW1b9++zG/ivffe07hx4666XfPmzXXo0CGH65x9rA0AAFQezn5/l2scnl9++UVz5szRv/71Lx05ciT/QGWYBLTglDVq1NDw4cP15JNPqnv37mUtw6UIPAAAVD1uDTwF8vLytGbNGn3wwQdat25dia0nhQUFBalbt24aPHiwRo4cqYiIiPKe3qUIPAAAVD0eCTxXOnbsmDZt2qSjR4/q1KlTOnPmjAIDA9WwYUM1bNhQHTp0UJcuXSrluDgEHgAAqh6Pj7QsSU2aNCmxwzAAAIC3mGrgQQAAAEcIPAAAwPQIPAAAwPQIPAAAwPRc2mkZqO6sVqtsNpuCgoJ42g8AKhECD+AiKSkpSk1Ntf8eGRmpqKgoL1YEACjALS3ABaxWa5GwI0mpqamyWq1eqggAUBiBB3ABm81WpuUAAM/yyC2tn3/+WfPmzdOGDRu0f/9+nT9/XpJUp04dtW7dWn379tWYMWPKNVM6UBkEBQWVaTkAwLNcOrWEIzNnzlRcXJwuX74sSQoNDbV35rRarTp9+rQkyc/PT9OnT9dTTz3lznJKxNQSqCj68ACA53llaokrLVmyRM8884zatWunZ599Vrfccovq1q1bZJtz587piy++0IwZMzR16lS1aNFC99xzjzvLAtwiKipKoaGhPKUFAJWQW1t4unfvroyMDG3fvl21atUqddvMzEzFxMQoNDRUmzdvdldJJaKFBwCAqsfZ72+3dlretWuX7rzzzquGHUkKDg7WnXfeqV27drmzJAAAUA25NfD4+/vbOyg7w2q1yt/f340VAQCA6sitgadHjx5aunSpduzYcdVtt2/friVLlqhnz57uLAkAAFRDbu3Dk5CQoN69e8vHx0ejRo3SgAED1Lp1a9WpU0eSdP78ee3fv19ffvmlFi9erLy8PG3cuFFdunRxV0klog8PAABVj7Pf325/LH39+vWaMGGCUlJSZLFYHG5jGIaioqL0zjvvqG/fvu4sp0QEHgAAqp5K8Vi6JPXt21c//fSTvv76a61fv97hwIN9+vTRTTfdpBo1ari7HAAAUA25vYWnqqCFBwCAqqdSPJYOAABQGVSqwLNixQo9//zz3i4DAACYTKUKPMuXL9f06dO9XQYAADCZShV4AAAA3MGtT2ktWLCgTNsfOHDATZUAAIDqzK2BZ+zYsSWOveOIYRhl2h4AAMAZbg08/v7+ioiI0KRJk5zaftmyZUpOTnZnSYCsVqtsNpuCgoIYggAAqgm3Bp4OHTroyJEj+stf/uLU9nv37iXwwK1SUlKUmppq/z0yMlJRUVFerAgA4Alu7bTcuXNnnT59usgXDOAtVqu12J/F1NRUWa1WL1UEAPAUt7bw9O7dW1988YX279+vyMjIq25/4403urMcVHM2m63E5dzaAgBzY2qJXzG1hPlZrVaHt0w7derENQeAKoqpJYArhISEFGtpjIyMJOwAQDXg9tnSgcokKipKoaGhPKUFANUMgQfVTkhICEEHAKoZlwSeI0eOlHkfi8WiwMBA1alTR/7+/q4oAwAAwCGXBJ5rrrmmQvs3bdpUPXr00NixYzVo0CBXlAQAAGDnkk7LhmFU6Cc1NVUffvihhgwZooEDB+rs2bOuKAsAAECSi1p45s2bJ0l64403tHnzZgUGBmrgwIHq0qWLGjZsKEk6deqUtm7dqi+++ELZ2dnq1q2bHnjgAVmtVu3cuVMrVqzQ6dOntXbtWt1+++3asGGDK0oDAABw3Tg8kydP1ttvv61hw4bpzTffVKNGjRxud/LkST3wwANauXKlxo0bp3feeUeSlJWVpcmTJ2vBggWyWCxavHix7r77bleU5hTG4QEAoOrx6Dg8y5cv11tvvaWePXvq448/LjHsSFKjRo30ySefqEePHpo3b54++OADSVJgYKD+/e9/KzY2VpK0dOlSV5QGAADgmsAzZ84cWSwWPf7447JYLFfd3mKx6A9/+IMMw9Bbb731v2J8fDRx4kQZhqGtW7e6ojQAAADXBJ4dO3ZIklq1auX0PgXb/vjjj0WWd+zYUZKUkZHhitIAAABc02m5YLbpU6dOOb1PwbaZmZlFlgcFBUmS/Pz8XFEanGC1Whl5GABgai5p4SmYn2jRokVO77Nw4UJJUrNmzYosP3nypCTZn+6Ce6WkpCg5OVl79+5VcnKyUlJSvF0SAAAu55LAM2zYMBmGoffff1+vvfbaVbd/9dVXtWjRIlksFg0bNqzIus2bN0uSmjdv7orSTM9qtSo9Pd3eylbWfVNTU4ssS01NLdexAACozFxyS+upp57S/PnzlZGRoSeeeEKLFi3S6NGj1blzZ/sTWydPntTWrVu1cOFCJSUlScpvxXnqqaeKHGvp0qWyWCwaMGCAK0oztZSUlCKBJTIyUlFRUU7vb7PZSlzOrS0AgJm4JPA0aNBAa9as0a233qr09HQlJSXZQ40jhmGocePG+vzzz1W/fn378pSUFHXr1k3dunXTiBEjXFGaaZXUOhMaGup0WCnoL+XscgAAqiqX3NKSpOjoaO3Zs0ePPvqoQkJCSpxGIiQkRI8++qh27dplfyKrQFRUlObNm6d58+apdevWrirNlEprnXFWSEiIvf9VgcjISFp3AACm47KRlgvLzs5WYmKidu7caZ8Xq169emrXrp26dOmigIAAV5+ywqraSMtWq1XJycnFlnfq1KnM9fOUFgCgqnL2+9stgacqqmqBR6p4Hx4AAKo6Z7+/XdKHB94RFRWl0NBQWmcAALgKtwWe9PR07dy5U2fOnJEk1a9fX+3bt1dYWJi7TlkthYSEEHQAALgKlwaegrmxZs+erd27dzvcpm3btnr00Uc1ceJEp+bdgufQlwcAYFYu68Nz9uxZDR06VN9//72k/PDj8IS/hpwbbrhBK1euVN26dV1x+gqrin14XIn+QACAqsijfXgMw9Dw4cO1adMmSfnj8vzud79T9+7d1bhxYxmGofT0dG3ZskUffPCBTp8+rU2bNmn48OHasGGDK0pABbhiTB8AACozl4zDs3jxYm3cuFEWi0WjRo1SSkqK/vWvf2n06NG65ZZbNHDgQI0ePVqzZ89WSkqK7r//fhmGoY0bN2rJkiWuKEEHDx7U22+/rYkTJyo6Olq+vr6yWCx64YUXXHJ8M3PFmD4AAFRmLmnhWbx4sSSpT58+9klBS1K7dm3Nnz9fR44c0YYNG/T+++/r3nvvrXANs2bN0qxZsyp8nOqIEZcBAGbnkhaepKQkWSwWPfLII07v8+ijj0qSw8HzyiM0NFRDhgzR888/r9WrVzM1RRkw4jIAwOxc0sJT8Oj5Nddc4/Q+BdsW7FtR06ZNK/L70qVLXXLc6oIxfQAAZuaSFp46depIko4fP+70PgXb8sVaeYSEhCgsLIxrAgAwHZcEnvbt20uS5s2b5/Q+//73v4vsC/eyWq1KT0+X1Wr1dikAAHicSwLPXXfdJcMw9Mknnyg+Pr7EMXik/EfY4+Pj9cknn8hisei3v/2tK0pAKVJSUpScnKy9e/cqOTlZKSkp3i4JAACPckkfnokTJ2r27Nnau3evZsyYoY8++khjx45V9+7dFRYWJovFohMnTmjz5s2aP3++du3aJUlq06aNJk6c6IoSXMbZFpCAgIBKOev7lRhjBwAAFwUePz8/rV69Wv3799fBgwe1e/du/fnPfy5xe8MwFBUVpdWrV8vXt3LNX3rl00oliYuLU3x8vHuLcYHSxtgh8AAAqguXpY3mzZtrx44dio+P17vvvqtz58453K5u3bqaMGGCnnvuOdWuXdtVp3eZ1NRUp4JAVWjdkVw/xg7zbQEAqiKXNq/UqlVLr7zyil588UUlJiY6nC29c+fO8vf3d+VpXcpss48XjLFz5TxZ5XmPzLcFAKiq3HI/yd/fXz179lTPnj3dcXiUkSvG2KEvEACgKqtcHWjgNhVtuaIvEACgKitT4Dly5IhbimjWrJlbjgvXYb4tAEBVVqbAU5apI5xlsViUk5Pj8uPCtVzZFwgAAE+zGKWNEngFHx+XjFNYtACLRbm5uRU+znfffafhw4fbf79w4YKys7NVs2bNIq0QycnJDh89t1qtqlOnjs6fP8+XeCl4SgsAUJk4+/1dphaeskwd4WmXL19WRkZGseUXL17UxYsX7b+7IlxVZ2Z7ig0AUD2UqYXHzGjhAQCg6nH2+9v196gAAAAqGQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPWZLryKY0gEAgPIj8FQBKSkpxSbtjIqK8mJFAABULdzSquSsVmuRsCNJqampslqtXqoIAICqh8BTydlstjItBwAAxRF4KrmgoKAyLQcAAMUReCq5kJAQRUZGFlkWGRlJx2UAAMqATstVQFRUlEJDQ3lKCwCAciLwVBEhISEEHQAAyolbWgAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPR8vV2AmVitVtlsNgUFBSkkJMTb5QAAgF8ReFwkJSVFqamp9t8jIyMVFRXlxYoAAEABbmm5gNVqLRJ2JCk1NVVWq9VLFQEAgMIIPC5gs9nKtBwAAHgWgccFgoKCyrQcAAB4FoHHBUJCQhQZGVlkWWRkJB2XAQCoJOi07CJRUVEKDQ3lKS0AACohAo8LhYSEEHQAAKiEuKUFAABMj8ADAABMj8ADAABMj8ADAABMj8ADAABMj8ADAABMj8ADAABMz3SB57///a9uvvlm1a9fX7Vq1VJsbKz++c9/Ki8vz9ulAQAALzFV4Jk5c6YGDx6sr776SvXq1VOrVq20fft2PfbYY7rjjjsIPQAAVFOmCTzff/+9nnnmGfn4+Gjx4sX6+eeftX37diUlJSksLEyffvqp/vGPf3i7TAAA4AWmCTwvvPCCDMPQhAkTdO+999qXR0dH24POzJkzdfnyZW+VCAAAvMQUgcdqtWrt2rWSpN///vfF1v/2t79VSEiIMjIytG7dOk+XBwAAvMwUgSc5OVmXLl1SYGCgYmNji6338/NT165dJUmbN2/2dHkAAMDLTDFb+v79+yVJzZo1k6+v47cUFRWlr776yr5tSaxWq1PnDAgIUEBAQNkKBQAAXmGKwHP27FlJUr169UrcpmBdwbYliYyMdOqccXFxio+Pd65AAADgVaYIPFlZWZIkf3//ErcpaI2x2WylHis1NVUhISFXPSetOwAAVB2mCDyBgYGSpEuXLpW4TXZ2tiQpKCio1GOFhIQ4FXgAAEDVYYpOy87crnLmthcAADAnUwSe1q1bS5KOHDminJwch9ukpKQU2RYAAFQfpgg8nTp1kp+fn7KyspSUlFRs/eXLl5WQkCBJ6t69u6fLAwAAXmaKwBMSEqKbb75ZkvTuu+8WW79s2TJZrVY1aNBAffv29XB1AADA20wReCRp6tSpslgseuedd7RkyRL78u3bt+uPf/yjJOnPf/5zqU9yAQAAc7IYhmF4uwhXefHFFzVt2jRJ+QMN1q5dWzt37lReXp4GDx6sFStWqEaNGg73tVqtqlOnjs6fP89TWgAAVBHOfn+bpoVHym/lWblypfr376+MjAwdOHBAHTp00GuvvVZq2AEAAOZmqhaeiqCFBwCAqqdatvAAAAA4QuABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+ABAACmR+BBqbKzsxUfH6/s7GxvlwIX4rqaF9fWnLiuFcfUEr9iagnH+FzMietqXlxbc+K6loypJQAAAH5F4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKbn6+0CKouC8RetVquXK6lcCj4PPhdz4bqaF9fWnLiuJSv4TK42jjIjLf/q6NGjioyM9HYZAACgHFJTU9W0adMS1xN4fpWXl6fjx48rODhYFovF2+UAAAAnGIahzMxMRUREyMen5J46BB4AAGB6dFoGAACmR+ABAACmR+ABAACmR+AxkYMHD+rtt9/WxIkTFR0dLV9fX1ksFr3wwgsVOu7333+v4cOHq2HDhgoKClLbtm01Y8YMZWVlOX2MtWvXymKxyGKx6Oabb65QPdVNZbyua9as0YgRIxQREaGAgAA1btxYffv21SuvvFKhmqqTynRdL126pFmzZqlHjx6qU6eO/Pz8FB4erjvuuENff/11heqpblx9XU+cOKEFCxbokUceUbdu3RQQECCLxaIJEyY4tf+ePXs0atQohYeHKzAwUC1bttSf/vQnnTt3rlz1VGkGTOPxxx83JBX7mTFjRrmP+f777xs1atQwJBlNmjQxOnXqZPj5+RmSjK5duxq//PLLVY9hs9mMVq1a2eu56aabyl1PdVSZrmteXp4xefJkew1NmzY1unbtarRo0cLw9fU1GjRoUO6aqpvKcl1/+eUXo2fPnvbzt2jRwoiNjTXq1q1rX/bSSy9V5K1WK66+rq+++qrD4/3+97+/6r5ff/21ERQUZEgyGjZsaMTGxho1a9Y0JBlRUVHGiRMnylVTVUULj4mEhoZqyJAhev7557V69WqNGDGiQsc7dOiQfv/73ys3N1cvv/yyUlNTlZSUpP379+u6665TQkKC/vznP1/1OC+88IIOHDigYcOGVaie6qoyXdepU6fqjTfeUPv27bVlyxalpqZqy5YtOnjwoDIyMjRv3rwK1VadVJbr+o9//EPff/+9GjZsqB9++EEHDx5UYmKiTp48qfj4eEnSM888owMHDlSovurC1dc1JCREAwYM0NSpU7VixQo9+uijTu2XmZmpu+++WzabTY899piOHTumxMREHTlyRL169VJKSop+//vfV6i2KsfbiQvuM2bMmAr9z+Khhx4yJBm33HJLsXXfffedIcnw8/Mr9X8Ju3fvNvz9/Y1BgwYZ8+bNo4XHBbx1XX/88UejRo0aRsOGDY309PRynRsl89Z17dGjhyHJeP311x0eNyYmxpBkzJkzp1x1VXcVva5XiouLc6qF5+WXXzYkGddff72Rk5NTZN3hw4cNX19fQ5KRmJjokrqqAlp44JBhGPrkk08kyeH/Am644Qa1adNGly9f1ooVK0o8xqRJk+Tj46PZs2e7tV44pyLXdfbs2crNzdXjjz+uRo0aeaReOKci19Vms0mSoqKiHB67ZcuWkqScnBxXlgw3+/jjjyVJY8eOVY0aNYqsa9asmb0v5Ycffujx2ryFwAOHjhw5orS0NElSr169HG5TsHzz5s0O17/77rv69ttv9fTTT5f4jyk8qyLXdeXKlZKkIUOGKCkpSQ8//LAGDBig4cOH669//atOnjzpxspRmopc144dO0qSNm3aVGyf7OxsJSYmSpK6du3qsnrhXjk5OfbrVt5/v82IwAOH9u/fL0kKCAhQRESEw20KQkzBtoWdOnVKf/nLX9SqVSv95S9/cV+hKJPyXtcTJ07o+PHjslgsWrdunbp166Y5c+Zo7dq1+vTTTzV16lS1bt1aa9eudf+bQDEV+fv61FNPqXbt2nrllVf0j3/8Q8eOHZPNZtO2bds0YsQIHTp0SPfdd5969Ojh3jcBlzl06JAuX74sqeSWu9L+/TYrAg8cOnv2rCSpbt26Jc4tVq9evSLbFjZlyhSdOXNGs2fPVkBAgPsKRZmU97oWtB5YLBY98cQT6tatm5KSkpSdna1du3ZpwIABslqtGjFihFJTU938LnClivx9bdu2rb777jsNGDBAf/rTn9S0aVPVrFlTnTp10g8//KB//vOfmj9/vnvfAFyq8DUuuO5XKu3fb7Mi8MChgjE7/P39S9ymIMgU9AEo8NVXX2nRokW66667NHDgQPcViTIr73X95ZdfJOVPslu7dm199tln6tSpk/z9/dW2bVutWLFCERERslqteu2119z3BuBQRf6+Svm3xNLT02UYhiIiIhQTE6PatWvbn7zbsWOHewqHWxQec6mkPxOl/XkwKwIPHAoMDJSUPyBZSbKzsyVJQUFB9mVZWVmaPHmyateurVdffdW9RaLMyntdC/aTpNGjRxf7X2NQUJAmT54sSfr8889dVi+cU97rKkmLFi3SsGHDdOzYMa1fv17Hjh1TcnKyMjIyNG3aNCUlJek3v/mNDh486L43AJcq/Pe1pD8TJf15MDMCDxwq+EI7d+6cDMNwuE1BU2jhL7+XXnpJBw4cUFxcnJo2ber+QlEm5b2uhV+3adPG4X7XX3+9pPz+A/Cs8l7Xy5cv64knnpBhGHrttdfUp08f+zp/f3/NmDFDt9xyizIzMzVz5kw3vgO4UuFrXNItK0d/HsyOwAOHWrduLSn/fwHHjx93uE1KSkqRbSUpOTlZkvTyyy+rcePGRX4ef/xxSdK3335rX0Z/D88q73Vt0aKFvQm8pD5ZBctzc3NdVi+cU97run//fqWnp0uSbrrpJof7FTy+vHXrVpfVC/dq0aKF/Pz8JP3vul/J0Z8HsyPwwKFmzZqpcePGkqTvvvvO4TYFy7t3715s3alTp5Senl7kx2q1SspvYi1YxpejZ5X3utaoUcP+WPLV/gFt0qSJy+qFc8p7XTMzM6967IIWo7LMnQfv8vX1VWxsrKTy/fttVgQeOGSxWHTHHXdIyh9P50qbNm3S3r175efnV2TKiOXLl8swDIc/BdMO3HTTTfZlLVq08Mj7Qb7yXldJ+t3vfidJWrJkif2R18IKnuTp37+/q8vGVZT3urZs2dL+VNdXX33l8NgFQw1ce+21ri4bbnTnnXdKkt57771i/7E8cuSI/bpWdOqLqoTAU8299tpratGihe65555i65588kn5+/vryy+/1CuvvGL/n97hw4c1fvx4SdKECRPs/7NE5eGO6zphwgRFRkbq0KFDevzxx+2dIXNzczV16lQlJyfL399fU6ZMcfO7q75cfV1DQ0PtT1L+4Q9/0DfffGNfd+nSJT377LNas2aNJOn+++932/uq7kq7ruU1efJkhYaGas+ePfrjH/9o/09KRkaGRo4cqZycHA0aNEidO3d22TkrPY9OZAG32rhxo9GgQQP7T0BAgCHJqFmzZpHlR44cse9TMC9Lnz59HB5z/vz5ho+Pj8PZlzt37mxcuHDB6fqYS6t8KtN1TUhIMEJCQgxJRr169YyuXbsaDRs2NCQZNWrUMN577z13fASmVFmu66FDh4xmzZrZZ+Fu0qSJERMTYwQHB9uXTZw40V0fg+m4+roeOXKkyH4Fs58HBAQUWb5x48Zi+65du9YIDAy0z5beuXNn+2zpLVq0MNLS0tz5UVQ6tPCYyOXLl5WRkWH/KXjs8OLFi0WWl6XfzOjRo/Xtt99qyJAhstls2r17t6KiohQfH6+NGzeqVq1a7no7+FVluq5dunTRjh07NGHCBNWqVUvbtm2TlN98vmnTJo0ZM6bC77e6qCzXtXnz5tq+fbvi4uLUqVMnnT9/Xjt37lRgYKAGDRqkjz76SG+99ZbL3rfZufq65ubmFtmvYNyc7OzsIssd3Wa+6aabtHXrVt1zzz2yWCz68ccfFRYWpj/+8Y9KSkqqdq3zFsMo4RlGAAAAk6CFBwAAmB6BBwAAmB6BBwAAmB6BBwAAmB6BBwAAmB6BBwAAmB6BBwAAmB6BBwAAmB6BBwAAmB6BB4BXvffee7JYLLJYLDp06FClOZY3zwHA9Qg8AEq0fv16+5d7fHy8t8sBgHIj8ACo1GhRAeAKBB4AXjV27FgZhiHDMNSiRYtKcyxPunDhgo4dO6ajR4/q2LFjSk9P17lz58o0UzqA0vl6uwAAqG7y8vL0n//8R//5z3+0fv16nT9/3uF2Pj4+atiwoa699lp16dJFd911l2644QYPVwuYAy08AOBBGzduVNeuXTV//nxFR0frscceU7169ezr27Vrp/T0dGVnZ+vYsWNauHChIiIi9Oqrr6pXr14aPXq08vLyvPgOgKqJwAOg3OLj4+39ayQpKytLr7zyimJjYxUcHKzg4GB169ZNs2fPVk5OjsNjlNRHp6DD9Lhx4+zLrrnmGvu2BT/r16+/6rEK27lzp1544QUNHDhQTZs2VUBAgGrXrq3WrVtrzJgx+uGHHyr8uZTkzTffVHx8vD766CN9/vnnmj59up5//nndfPPN9m0mTZqkRo0ayd/fX40bN9aAAQO0dOlS++ewcOFC/fvf/3ZbjYBZcUsLgEukp6dr4MCB2r59e5HlCQkJSkhI0Jdffqnly5fLx8d7/89av369+vXrV2z5pUuXdODAAR04cEALFizQU089pb/97W8uP390dLQmTJigGjVqFFm+c+dO++uhQ4c63Hf8+PGaN2+eJOnLL7/UhAkTXF4fYGYEHgAuceedd2rPnj167LHHNHToUNWvX18//fSTZsyYoT179mjlypV6++23NWnSJKeO17VrV/34449asWKFpk2bJkn64osvFBERUWS7a665xukac3JyVKtWLQ0ePFj9+/dXmzZtFBISopMnT2rXrl16/fXXdfjwYc2cOVPXXnttkdYlV+jRo0exZSdOnNCePXskSddff32Jna0bNWpkf22z2VxaF1AdEHgAuERBK07fvn3ty2JjYzVw4EC1bdtW6enpmjNnjtOBp1atWmrfvr22bt1qX3bttddW6OmrmJgYHT16VHXr1i22buDAgXrkkUc0ZMgQrVmzRtOnT9fo0aOLtca42hdffGF/fdttt5W43cmTJ+2vyxLyAOSjDw8Al3j00UeLhJ0C9evXt7eU7Nixo8QnkjwhNDTUYdgp4O/vr1deeUWSdPjwYW3bts3tNX322Wf216UFnnXr1tlf9+7d2601AWZECw8Alxg1alSJ6zp37mx/ffDgQcXExHigoqvLzs5Wenq6Lly4YH/yyTAM+/rt27cXqd3VLl++bG/hCQ4OLjHIZGZmau7cuZLyb3vdfvvtbqsJMCsCDwCXaNOmTYnr6tevb3+dmZnpiXJK9Msvv+j111/X0qVLtWvXrlIH9zt9+rRba9mwYYOsVqskacCAAfLz8yu2TV5ensaPH6+0tDQFBQVpwYIFDrcDUDoCDwCXqFmzZonrCj+Z5c3Rgw8dOqT+/fvr4MGDTm3v7s7BK1eutL92dDsrPT1dY8eO1eeff6769evrww8/VJcuXdxaE2BW9OEBUG3cf//9OnjwoCwWi8aPH68vv/xSqampysrKsk9JUTiQFb695Q6rVq2yvx40aJD9nFu2bNEf//hHtWrVSt98840mTZqkXbt2OXykHoBzaOEBUC3s3btXGzdulCQ9/fTTevHFFx1ud/bsWY/Us3v3bqWkpEjKf3osMDBQQ4YMUVpammrXrq2YmBi98847GjJkiGrVquWRmgAzI/AAqNQKRnGuqF27dtlf33PPPSVuV/gxeHe68nZW/fr1i7T4AHAtbmkBqNQCAwPtr7Ozs8t9nMJTW1y8eLHE7d54441yn6MsrtZ/B4BrEXgAVGrh4eH21z///HO5j9O6dWv76/nz5zvcZu7cuVq+fHm5z+Gs06dP2+fsatCggcMRmAG4Fre0AFRqnTp1UmBgoLKysvTss8/K19dXLVq0sD/51aRJEwUFBTl1nPbt22vnzp2aO3euzp07p1GjRik8PFypqal6//339eGHH6pXr1767rvv3PqePvjgA3vn6IceesjtozkDIPAAqOSCg4P12GOP6eWXX1ZSUpIGDhxYZP26descjvB8JYvFooULF6p///46e/aslixZoiVLlhTZpkOHDlq2bFmx+bpc6ezZs/q///s/SfkjJj/55JNuOxeA/+GWFoBKb+bMmXr77bfVu3dv1a9fv9wtIjExMdq2bZsmT56s5s2by8/PT/Xr11e3bt3097//XVu2bClyC83VNm7cqN69e6tx48b617/+pS+//FLBwcFuOx+A/7EY7h5oAgAAwMto4QEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKb3/w2dDFEGWOhDAAAAAElFTkSuQmCC",
+      "text/plain": [
+       "<Figure size 640x480 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "# plot\n",
+    "plt.scatter(P_ratios, log_t_insts, c='silver', s=10)\n",
+    "plt.tick_params(labelsize=16,size=8,direction='in')\n",
+    "plt.xlabel(r'Initial $\\mathcal{P}$', fontsize=20)\n",
+    "plt.ylabel(r'$\\log_{10}(t_\\mathrm{inst}/P_1)$', fontsize=20)\n",
+    "plt.title(r'Model (2, 1)', fontsize=20)\n",
+    "plt.xlim([1.028, 1.112])\n",
+    "plt.ylim([-0.4, 4.5])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a07d21cc",
+   "metadata": {},
+   "source": [
+    "Extending beyond Model (2, 1) to include higher-order MMRs can be done by simply adding the indices of the desired resonances to `res_inds`. Secular terms can be added to the Hamiltonian models using `pham.add_secular_terms`."
+   ]
+  }
+ ],
+ "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.11.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}