PREM_normal_modes_example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import subprocess # We are using python 3.5 features! \n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import earth_model.earth_model as earth_model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some Mineos wrapping stuff\n",
    "\n",
    "Nothing much to see here - probably belongs in a module"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "MINEOS_BIN_FILE = '../cig/bin/minos_bran' # set for your setup!\n",
    "\n",
    "def run_mineos(earth_model, model_name, mode_type, \n",
    "               lmin, lmax, wmin, wmax, nmin, nmax,\n",
    "               eps=1E-10, wgrav=10, model_spacing=25):\n",
    "    \"\"\"\n",
    "    Run Mineos\n",
    "    \n",
    "    Given a tabulated Earth model, set up a Mineos model file\n",
    "    and a Mineos parameter file, and run the minos_bran code\n",
    "    to compute mode frequencies, eigenvectors etc.\n",
    "    \n",
    "    Arguments:\n",
    "    earth_model: an Earth model (e.g. a Prem instance)\n",
    "    model_name: a string describing the name of the model\n",
    "    mode_type: 'radial', 'toroidal', 'spheroidal' or 'inner_core_toroidal'\n",
    "    lmin: minimum angular order\n",
    "    lmax: maximum angular order\n",
    "    wmin: minimum frequency (mHz)\n",
    "    wmax: maximum frequency (mHz)\n",
    "    nmin: minimum dispersion branch number\n",
    "    nmax: maximim dispersion branch number\n",
    "    eps: controls the accuracy of the Runge-Kutta integration scheme\n",
    "    wgrav: the frequency in millihertz (mHz) above which gravitational\n",
    "           terms are neglected\n",
    "    model_spacing: maximum spacing (in km) used to represent the model\n",
    "    \"\"\"\n",
    "    if mode_type == 'radial':\n",
    "        code_letter = 'S' # this does l=0 spheroidal modes\n",
    "        jcom = 1\n",
    "    elif mode_type == 'toroidal':\n",
    "        code_letter = 'T'\n",
    "        jcom = 2\n",
    "    elif mode_type == 'spheroidal':\n",
    "        code_letter = 'S'\n",
    "        jcom = 3\n",
    "    elif mode_type == 'inner_core_toroidal':\n",
    "        code_letter = 'T' # should I call this something else?\n",
    "        jcom = 4\n",
    "    else:\n",
    "        raise ValueError(\"invalide mode type\")\n",
    "    \n",
    "    name = 'mineos_run'\n",
    "    param_file_name = name + '_param.txt'\n",
    "    model_file_name = name + '_model.txt'\n",
    "    text_output_file_name = name + '_text_out' + code_letter + '.txt'\n",
    "    unformatted_output = name + 'eig_out' + code_letter + '.dat'\n",
    "    \n",
    "    # Tabulate and write out model\n",
    "    table = prem.tabulate_model_outwards(model_spacing)\n",
    "    write_mineos_model_file(table, model_file_name, model_name)\n",
    "    \n",
    "    # Run Mineos, piping in commands\n",
    "    mineos = subprocess.Popen([MINEOS_BIN_FILE],\n",
    "                              stdin=subprocess.PIPE,\n",
    "                              universal_newlines=True,\n",
    "                              bufsize=0)\n",
    "    # Write param file\n",
    "    mineos.stdin.write(model_file_name + \"\\n\")\n",
    "    mineos.stdin.write(text_output_file_name + \"\\n\")\n",
    "    mineos.stdin.write(unformatted_output + \"\\n\")\n",
    "    mineos.stdin.write(\"{} {}\\n\".format(eps, wgrav))\n",
    "    mineos.stdin.write(\"{}\\n\".format(jcom))\n",
    "    mineos.stdin.write(\"{} {} {} {} {} {}\\n\".format(lmin, lmax, \n",
    "                                                    wmin, wmax, \n",
    "                                                    nmin, nmax))\n",
    "    mineos.stdin.close()\n",
    "    mineos.wait()\n",
    "    \n",
    "    modes = read_mineos_mode_file(text_output_file_name)\n",
    "    return modes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def write_mineos_model_file(table, mineos_model_filename, mineos_model_name):\n",
    "    \"\"\"\n",
    "    Write out a formatter Mineos model\n",
    "    \"\"\"\n",
    "    f = open(mineos_model_filename, 'w')\n",
    "    f.write(mineos_model_name + \"\\n\")\n",
    "    f.write(\"   0    1.00000  1\\n\") # isotropic tabulated format 1s ref\n",
    "    n = table.shape[0]\n",
    "    nic = np.argwhere(table.vs==0)[0][0] # This is the last inner core index. \n",
    "    noc = np.argwhere(table.vs==0)[-1][0] + 1 # This is the last outer core index (assuming no ocean).\n",
    "    f.write(\"{} {} {}\\n\".format(n, nic, noc))\n",
    "    for r, rho, vp, vs, qk, qm in zip(table.radius, table.density, table.vp, table.vs, \n",
    "                              table.qkappa, table.qshear):\n",
    "        f.write(\"{:#8.0f} {:#9.2f} {:#9.2f} {:#9.2f} {:#9.1f} {:#9.1f} {:#9.2f} {:#9.2f} {:#9.5f}\\n\".format(\n",
    "                                                                   r*1000.0, rho, \n",
    "                                                                   vp*1000.0, vs*1000.0,\n",
    "                                                                   qk, qm, vp*1000.0, vs*1000.0, 1.0))\n",
    "    f.close()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_mineos_mode_file(filename):\n",
    "    \"\"\"\n",
    "    Read output data from a Mineos run\n",
    "    \n",
    "    Returns results in dictionary with mode names as keys\n",
    "    \"\"\"\n",
    "    read_modes = {}\n",
    "    in_modes = False\n",
    "    f = open(filename, 'r')\n",
    "    for line in f:\n",
    "        if not in_modes:\n",
    "            words = line.split()\n",
    "            if not words:\n",
    "                # blank line\n",
    "                continue\n",
    "            elif words[0] == 'mode' and words[1] == 'phs':\n",
    "                # One blank line then into the modes\n",
    "                in_modes = True\n",
    "        else:\n",
    "            words = line.split()\n",
    "            if not words:\n",
    "                continue\n",
    "            this_mode = {}\n",
    "            n = words[0]\n",
    "            letter = words[1].upper()\n",
    "            l = words[2]\n",
    "            this_mode['n'] = int(n)\n",
    "            this_mode['type'] = letter\n",
    "            this_mode['l'] = int(l)\n",
    "            this_mode['phs_vel'] = float(words[3])\n",
    "            this_mode['w'] = float(words[4])              \n",
    "            this_mode['t'] = float(words[5])\n",
    "            this_mode['grp_vel'] = float(words[6])\n",
    "            this_mode['q'] = float(words[7])\n",
    "            this_mode['raylquo'] = float(words[8])\n",
    "            read_modes[n+letter+l] = this_mode\n",
    "                \n",
    "    return read_modes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup and run a model with Mineos\n",
    "\n",
    "And write in Mineos format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prem = earth_model.Prem()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_modes = run_mineos(prem, 'a model name', mode_type='spheroidal', lmin=1, lmax=250, \n",
    "               wmin=0.0, wmax=166.0, nmin=0, nmax=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mode 0S2 has period 3230.644 and is 3233.115 secs in PREM Table V\n"
     ]
    }
   ],
   "source": [
    "print(\"Mode 0S2 has period\", s_modes['0S2']['t'], \"and is 3233.115 secs in PREM Table V\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_modes = run_mineos(prem, 'a model name', mode_type='toroidal', lmin=1, lmax=250, \n",
    "               wmin=0.0, wmax=166.0, nmin=0, nmax=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mode 0T2 has period 2631.841 and is 2636.38 secs in PREM Table V\n"
     ]
    }
   ],
   "source": [
    "print(\"Mode 0T2 has period\", t_modes['0T2']['t'], \"and is 2636.38 secs in PREM Table V\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl4AAAJRCAYAAACUd7lTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOy9f5xdVX3v/VlzkhMCQr13aJWCSSAv0EBGIBnAbUqy41CsBH9AqJTqE/BXKmqv2JcIsXLhis74gPXh+nArxF7QXH0ohVysir3FTrNDYA6UgVJiQcVUKAEtMC0afs2enPk+f5yZyQRm5qw1s9f5rrXP5/16zevMjz17r732j/VZ31/LiAgIIYQQQoh/OrQbQAghhBDSLlB4EUIIIYS0CAovQgghhJAWQeFFCCGEENIiKLwIIYQQQloEhRchhBBCSIuYp90AGw499FBZsmSJdjMIIYQQQppy//33PysivznV36IQXkuWLMHg4KB2MwghhBBCmmKMeXy6v9HVSAghhBDSIii8CCGEEEJaBIUXIYQQQkiLiCLGaypGRkawe/duvPzyy9pNIS3ggAMOwBFHHIH58+drN4UQQgiZNdEKr927d+Pggw/GkiVLYIzRbg7xiIhgaGgIu3fvxpFHHqndHEIIIWTWROtqfPnll9HZ2UnR1QYYY9DZ2UnrJiGEkOiJVngBoOhqI3itCSGElIGohVdMXHDBBbj11lu97PuMM87Ac88996rfX3HFFfjyl7884//abDNXvvGNb+ATn/iE12MQQgghMRBtjFc7ISIQEXR0TK2Tf/CDH7S4RYQQQgiZDbR4zZIXXngB69atw/HHH4/ly5fj5ptvBtCosn/JJZfg5JNPxsknn4yf/exnE/9z55134q1vfSuOOuqo/axfV199NU466SS8+c1vxuWXXw4AeOyxx7Bs2TJ87GMfw4oVK/DEE0/gpptuQldXF5YvX45LLrlk4v+XLFmCZ599FgDwxS9+EW984xtx2mmn4Sc/+cnENl//+tdx0kkn4fjjj8f69evx4osvznh+F1xwAS688EKsXbsWRx11FLZv344PfvCDWLZsGS644IKJ7aZr04033ohjjjkGa9aswd133z3x+2eeeQbr16/HSSedhJNOOmnib9u3b8cJJ5yAE044ASeeeCL27NljfS0IIYSQaBi3phT9BeAGAE8D+NErfv/HAH4C4J8BXGWzr5UrV8orefjhh1/1u6YMDIj09jY+58itt94qH/7whyd+fu6550REZPHixfKFL3xBRES++c1vyrp160RE5Pzzz5dzzjlH6vW6/PM//7MsXbpURET+9m//Vj7ykY/I6Oio1Ot1WbdunWzfvl1+/vOfizFGarWaiIg8+eST8oY3vEGefvppGRkZkbVr18ptt902ccxnnnlGBgcHZfny5fLCCy/Ir371K1m6dKlcffXVIiLy7LPPTrT1T//0T+WrX/2qiIhcfvnlE9tM5vzzz5dzzz1XRkdH5Tvf+Y4cfPDB8tBDD0m9XpcVK1bIP/7jP07bpqeeemri98PDw/LWt75VPv7xj4uIyHnnnSc7duwQEZHHH39c3vSmN4mIyJlnnil33XWXiIjs2bNHRkZGXtWmWV1zQgghpMUAGJRpNI1PV+M3AFwLYMv4L4wxawG8G8CbRWTYGPNbHo+/P7Ua0NMD5DlQrQL9/UCSzHp3XV1d+PSnP41LLrkEZ555Jk499dSJv5133nkTn5/61Kcmfv+e97wHHR0dOPbYY/Fv//ZvAIA77rgDd9xxB0488UQAwPPPP49HH30UixYtwuLFi/GWt7wFAHDfffchTVP85m821tx83/vehzvvvBPvec97Jva/Y8cOnHXWWTjwwAMBAO9617sm/vajH/0In/vc5/Dcc8/h+eefx9vf/vam5/jOd74Txhh0dXXhda97Hbq6ugAAxx13HB577DE8/vjjU7YJwH6/P/fcc/HTn/4UAPB3f/d3ePjhhyeO8etf/xp79uzBqlWr8Cd/8id43/veh7PPPhtHHHFE0/YRQgghseHN1SgidwL491f8+kIAXxKR4bFtnvZ1/FeRZQ3RVa83PrNsTrs75phjcP/996OrqwubNm3C5z//+Ym/Tc7Am/z9ggULJr5vCOLG56ZNm/Dggw/iwQcfxM9+9jN86EMfAgAcdNBBr9q+GdNl/11wwQW49tprsXPnTlx++eVWpRnG29vR0bFf2zs6OrB3794Z2zRdO0ZHR1Gr1SbO98knn8TBBx+MSy+9FH/xF3+Bl156CW95y1vw4x//uGn7CCGEkNhodYzXMQBONcbca4zZbow5qWVHTtOGpatSaXym6Zx299RTT+HAAw/E+9//fnz605/GAw88MPG38Xivm2++GUkTq9rb3/523HDDDXj++ecBAE8++SSefvrVevSUU07B9u3b8eyzz6Jer+Omm27CmjVr9ttm9erVuO222/DSSy9hz549+N73vjfxtz179uCwww7DyMgIvv3tb8/6vG3adMoppyDLMgwNDWFkZAS33HLLxP+cfvrpuPbaayd+fvDBBwEAu3btQldXFy655BJ0d3dTeBFCCCklrc5qnAfgPwF4C4CTAPyVMeYomcJ0YozZCGAjACxatGjuR06Shnsxyxqiaw5uRgDYuXMnLr74YnR0dGD+/Pn42te+NvG34eFhnHLKKRgdHcVNN900435OP/10PPLIIxMC7TWveQ2+9a1voVKp7LfdYYcdhr6+PqxduxYigjPOOAPvfve799tmxYoVOPfcc3HCCSdg8eLF+7k/r7zySpxyyilYvHgxurq6Cglen6lNV1xxBZIkwWGHHYYVK1agXq8DAL761a/i4x//ON785jdj7969WL16Na677jpcc8012LZtGyqVCo499li84x3vmHP7CCGEkNAwti6sWe3cmCUAvi8iy8d+/j9ouBqzsZ93AXiLiDwz0366u7tlcHBwv9898sgjWLZsmYdWz40lS5ZgcHAQhx56qHZTSkeo15wQQgiZjDHmfhHpnupvrXY1fgfA2wDAGHMMgCqAZ1vcBkIIIYQQFby5Go0xNwFIARxqjNkN4HI0SkzcYIz5EYAcwPlTuRlj5rHHHtNuAiGEEEICxZvwEpHzpvnT+30dkxBCCCEkZFi5nhBCCCGkRVB4EUIIIYS0CAovQgghhLQFtRrQ19f41ILCa5Y899xz+PM///NC9/nWt751yt9fcMEF+y2qPdtt5soVV1yBL3/5y16PQQghhPhgfOXAyy5rfGqJLwqvWTIb4TVeRHQ6BgYG5tIkQgghhExDwSsHzhoKr1ly6aWXYteuXTjhhBNw8cUXQ0Rw8cUXY/ny5ejq6ppYNijLMqxduxZ/+Id/OLHI9Fe+8hUsX74cy5cvxzXXXDOxz9e85jUAGusyfuITn8Cxxx6LdevW7beE0Oc//3mcdNJJWL58OTZu3Nh0Dcc0TfGpT30Kq1evxrJly3Dffffh7LPPxtFHH43Pfe5zE9tN16YvfvGLeOMb34jTTjsNP/nJTyZ+v2vXLvze7/0eVq5ciVNPPXViiZ9bbrkFy5cvx/HHH4/Vq1fPtnsJIYSQQil45cDZIyLBf61cuVJeycMPP/yq3zVjYECkt7fxOVd+/vOfy3HHHTfx86233iqnnXaa7N27V375y1/KG97wBnnqqadk27ZtcuCBB8q//Mu/iIjI4OCgLF++XJ5//nnZs2ePHHvssfLAAw+IiMhBBx0kIiJbt26d2NeTTz4pv/EbvyG33HKLiIgMDQ1NHPP973+/fPe73xURkfPPP39im8msWbNGPvOZz4iIyDXXXCOHHXaYPPXUU/Lyyy/L4YcfLs8+++y0bRr//QsvvCC/+tWvZOnSpXL11VeLiMjb3vY2+elPfyoiIvfcc4+sXbtWRESWL18uu3fvFhGR//iP/5h7R09iNtecEEIIGadIHTATAAZlGk3T6rUa1Rj37eZ5Q+n29895ucb9uOuuu3DeeeehUqngda97HdasWYP77rsPhxxyCE4++WQceeSRE9udddZZOOiggwAAZ599Nnbs2IETTzxxYl933nnnxL5++7d/G29729sm/rZt2zZcddVVePHFF/Hv//7vOO644/DOd75zxra9613vAgB0dXXhuOOOw2GHHQYAOOqoo/DEE09M26bR0VGcddZZOPDAA/fbz/PPP4+BgQH8/u///sQxhoeHAQCrVq3CBRdcgPe+9704++yzZ9+hhBBCSMEkSbFj/2xoG+E1lW+3yM6XGVx+44Km2XaTMca86ncvv/wyPvaxj2FwcBBveMMbcMUVV+Dll19uuq8FCxYAADo6Oia+H/957969M7ZpqnaMjo7ita99LR588MFX/e26667Dvffei9tvvx0nnHACHnzwQXR2djZtIyGEEDIbarXGmJ6m+qLKhraJ8Srat3vwwQdjz549Ez+vXr0aN998M+r1Op555hnceeedOPnkk1/1f6tXr8Z3vvMdvPjii3jhhRdw22234dRTT33VNn/5l3+Jer2OX/ziF9i2bRsATIisQw89FM8//3xhWYzTtWn16tW47bbb8NJLL2HPnj343ve+BwA45JBDcOSRR+KWW24B0BCT//RP/wSgEft1yimn4POf/zwOPfRQPPHEE4W0kRBCCHkloWQqutA2Fq8kabgXi1LFnZ2dWLVqFZYvX453vOMduOqqq1Cr1XD88cfDGIOrrroKr3/96yeCzsdZsWIFLrjggglR9uEPf3g/NyMAnHXWWfj7v/97dHV14ZhjjsGaNWsAAK997WvxkY98BF1dXViyZAlOOumkuZ2ERZvOPfdcnHDCCVi8ePF+AvHb3/42LrzwQnzhC1/AyMgI/uAP/gDHH388Lr74Yjz66KMQEfT09OD4448vpI2EEELIK/HtzfKBsXV9adLd3S2Dg4P7/e6RRx7BsmXLlFpENOA1J4QQMhnf8duzxRhzv4h0T/W3trF4EUIIIaRcFO3NagUUXoQQQgiJlhAyFV1om+B6QgghhMRBCGsq+iJqi5eITFnugJSPGGIRCSGEzJ1Q47aKIlqL1wEHHIChoSEOyG2AiGBoaAgHHHCAdlMIIYR4JpQ1FX0RrcXriCOOwO7du/HMM89oN4W0gAMOOABHHHGEdjMIIYR4Zrzu5rjFS21NRU9EK7zmz58/sQwPIYQQQspBjJmKLkQrvAghhBBSTmLLVHQh2hgvQgghhMRFmbMVbaHFixBCCCHeKXu2oi20eBFCCCHEO2XPVrSFwosQQggh3hnPVqxU9LIVQ3B10tVICCGEEO9oZyuG4uqk8CKEEEJIS9DMVpzK1anRFroaCSGEEDJrQnDf2RCCqxOgxYsQQgghsyQU950N2q7OcSi8CCGEEDIrQnHf2RJCYVa6GgkhhBAyK0Jw38Xi6hyHFi9CCCGEzApt911Mrs5xKLwIIYQQMmuYqegGXY2EEEIIeRUxuPBCcHW6QosXIYQQQvYjFheetqtzNlB4EUIIIWQ/tF14tZq9mAohU9EFCi9CCCGE7Me4C2/c4tVKF14s1rbZQuFFCCGEkP3QdOFpW9t8Q+FFCCGEkFeh5cLTtLa1AgovQgghpE1wiZ3SIsaAeRcovAghhJA2IKbYqdgC5l1gHS9CCCGkDZgqdqrVxFAbzDe0eBFCCCFtgHbsVEwWN59QeBFCCCFtgHbsVNmzFW2h8CKEEEIixzZoXjN2StviFgoUXoQQQkjExOLC07a4AWFkdVJ4EUIIIRETkwtP0+IWikBlViMhhBASMeMuvEpFL2g+hkzFELI6AVq8CCGEkKjRdOGFYkWyIZQYMwovQgghJHK0XHixuTm1Y8wACi9CCCEkSEIIBG9GKFYkW0KoiE/hRQghhARGLC68EKxIMQjUyVB4EUIIIYERmwuv3TMVXWBWIyGEEBIY2pmKQBzZiqFkKrpAixchhBASGNouvFgsSbHFmAEUXoQQQkjLcIlH0nThxeLq1Baos8Gb8DLG3ADgTABPi8jyV/zt0wCuBvCbIvKsrzYQQgghoRCLFQnQtyR5E6gBROL7tHh9A8C1ALZM/qUx5g0AfhfAv3o8NiGEEBIUsViRgJIWZQ1E+XoLrheROwH8+xR/+n8AfAaA+Do2IYQQEhohBMy7kCTApk2t1ybeAuYDicRvaYyXMeZdAJ4UkX8yxrTy0IQQQogqIcQjBeBpa4o3N6e2/3QMI+LP8GSMWQLg+yKy3BhzIIBtAE4XkV8ZYx4D0D1djJcxZiOAjQCwaNGilY8//ri3dhJCCCFzIQZBE4inzQpv/dmiC2WMuV9Euqf6WystXksBHAlg3Np1BIAHjDEni8gvX7mxiGwGsBkAuru76ZYkhBASJLEImthizLwEzAewZlDLhJeI7ATwW+M/N7N4EUIIITEQi6AJwdNWuMEpFtU7CZ/lJG4CkAI41BizG8DlIvI/fR2PEEII0SAEQWODdoyZF40Ui+qdhDfhJSLnNfn7El/HJoQQQlqFtqBxoXRFWWNRvZNg5XpCCCFkCmIKHYohuN+LRopJ9Y5B4UUIIYS8gphCh2JpqzeNpK16HfFWQJUQQgiJlUBqbVoRU1uti7LWakBfX+OzZNDiRQghhLyCmEKHtNsaU6ZiCC5ZCi9CCCHkFcQUOlS6dRU9ZSqG4pKl8CKEENJW2Fo9Ygod0mprTJmKoVSeoPAihBDSNoRi9bAhBLdYM2LKVNR2yY5D4UUIIaRtCMXq0YxYBGJMmYqhuI8pvAghhLQNoVg9mhGCQPTikvVgxoup3hpA4UUIIaQEuIiEEKwezdAWiF4sbh526rzLAPy3FF6EEEKixnXwDcHq0QxtgejF4uZhp067DMR/ywKqhBBCoiamAqIuWBcb9cC4xa1SKdDi5mGnTrsM5EahxYsQQkjUaLvlXAjA02WFF4ubh5067TKQG8WIiMqBXeju7pbBwUHtZhBCCAmUGARNIJ6u4omh88dpUVuNMfeLSPdUf6PFixBCSLCUqdhpTJmKTjv0oCa96aMAbhQKL0IIIUFSNguRtqcrluV9YsxUdIHCixBCSJCEYCEqklJmKnpQkzFmKrpA4UUIISRItC1EPtD0dMWyvI9TOyNU5xRehBBCWkrZip0CcXi7YlneJ8ZMRReY1UgIIaRlROgZakoZzwmA+vI+ujudG8xqJIQQEgQReoaaEsI5xZCt6E2gBpCp6AIr1xNCCGkZXiqiK6N9TuOC5rLLGp+1WgE79VDl3WmXtRrQ11fQyYQFLV6EEELmTBnjtmzRPqdYshWtd1la320DCi9CCCFzooyLVANu7jtmKxa4yxB8tx6h8CKEEDInyjhOxmR0cdJIymrSapcRZiq6QOFFCCFkTpRxnIxNTFoJmliW99H23XqGwosQQsicKOM4GYKYLFzQxLS8jy/fbQClJyi8CCGETIv2OKmFtpj0Ypxq9+V9ajXU0k3IRlYhnb8JSdanctNSeBFCCJkS7XFSG00x6cXV2ebL+9S2PIqe/AfIUUU1z9G/5VYkFF6EEEJCIbY4J1sC8DY1xZurs42X98mwBjmqqGMecggyrIHG5afwIoQQMiUhxDkVTSxWPGfjlKKatNZyyv7bdMNiVG+sI8/rqFY7kG5Y3NLjj0PhRQghbUY7FzuNyYpnLWhiyVYEVP23SQL0b6uo388UXoQQ0kaUtdipLdpWPC9iRjtbUdt363D8EO5nCi9CCGkjYrL4+EDTiufNzamZrajtuw0kU9EFCi9CCGkjtC0+vohheR9volczW5GZis5QeBFCSEmwER9ljNvSNrrY4ix6FdWk9X3CTEVnKLwIIaQEuIiPEOJciiQW96nzmorKAfNW9wkzFZ2h8CKEkBIQi/jwQUzuU2vRqx0w7wIzFZ2g8CKEkICxtVDEJD6KpozuU/XlfQDVbMXa5p3Itg4hXd+JZGPXjNvGZsGl8CKEkEBxdR+WTnwgnrUiC9co2sv7KAbO1TbvRM8fLUWOZajekaMfO5uKr5ig8CKEkEBxtVBoi4+iiSVo3lvNK83lfRR919nWIeRYti9gfusQko0tOXRLoPAihJBAaWf3IRBP3Fo0Na/goOUUb750fSeqd+TIIahiBOn6zpYduxVQeBFCiALtWvrBhViEZwg1r2JwdVofemMX+mEf4xUbRkS029CU7u5uGRwc1G4GIYQUQgCGj2jQXo3GFqt2eiwREc3yPm2CMeZ+Eeme6m+0eBFCSIuJxYXmixiqzLuiWfMqJlenS7ZiWaHwIoSQFhOLC80HAYz91ngxDnlQkiG4Om0oe7aiLRRehBBSEC6lD9o1disWa5+zQFR04UWzvE/JsxVtofAihJACcB2oY3GhFU0s1j4ngRiAGS+K5X1CyFYMIMaNwosQQgogFkuONrFY+5wEoqeLH4ur0/rQ2tmKAQhkgMKLEEKaYjMAxmLJ8UXZAuadBKKHix+Tq9NpeZ+NXXruxUBmRxRehBAyA7YDYCyWHB8EYkgoHGuB6OHix+LqjCpgPpDZEYUXIYTMgMsAGIMlxweBGBKssTYOKZrxQnB12hBVwHwgsyMKL0IImYFAJslBE1MfWRuHlM142q5OW4IImHchgNkRhRchpC1h6YfiiKmPrI1DAZjxNF2dtqgHzANOlskAkhq5ZBAhpP0oa0ySD0IYqIpE2+JVtv5Ux+E6tfK5V1kyyBhzA4AzATwtIsvHfnc1gHcCyAHsAvABEXnOVxsIIWQqAjBmREEZBaq1cciDFSmmTEUgkuV9HB7mUJ57n67GbwC4FsCWSb/7IYBNIrLXGPN/A9gE4BKPbSCEtBG241RMMUmahDJQFU2CGhJkAFIAratyG0umIhBRtqLDwxzKc+9NeInIncaYJa/43R2TfrwHwDm+jk8IaS9cxqmYYpI0CWWgssEpU1FJ0MSSqQhElK3o8DCH8txrBtd/EMDNiscnhJQI13EqgOSm4AlloGqGk5ZSFDSxZCoCAWQrugTMI0GGpJn9EkAYz72K8DLG/CmAvQC+PcM2GwFsBIBFixa1qGWEkBBh5fhiKVuVeSctpXyjxJCpCChnKwYaMF8ULRdexpjz0Qi675EZUipFZDOAzUAjq7FFzSOEBAYrxxdLjANVM5y0lKcbJZZ1FZ1Et9byPhEGzLvQUuFljPk9NILp14jIi608NiEkTlg5vlhiHKia4aylCr5RnMSs5pqKNaBnbR15blCtCvq3VcK89hEGzLvgs5zETWikjBxqjNkN4HI0shgXAPihMQYA7hGRj/pqAyEkXJiBqENZ+9M6U9ED1mJW2dyYbXkc+fDhqKOCfHgE2ZbdSJLFLTu+NREGzLvgM6vxvCl+/T99HY8QEg/MQPSDjZiNrT+tBLqyoLEWs8rmxhTbUcU5+wLmsR3AhpYdH4DbkhGWfRObpZtLBhFCWg4zEIvHVczG0J/W56QsaKzFrHZg/4aj0X/DGchGViGdfzeSDX0tPb5r0HwskwNXKLwIIYXCDEQdyhi7ZX1OAdxQVmJW29yYJEiyPiRZBqR9rT++5QUtYwLIZCi8CCGFwQxEPQLQHoVjfU4hZCp6cKG5YN1WTXOn5QUt4yRiMhRehJDCYAaiHmUUs0kC9F8zqZZUMkMtKe1MRc3lfWLJVrS8Scs4iZgMhRchpCnMQNSjbMVOnajVkFzUgyTPgR1VoKt1gsbJ6qJsookmWxGwukm9TiICCB6j8CKEzAgzEPUoa6yL9dinKGicJhHay/toZyu6LO+j6ZEN5IGi8CKEzAgzEPUoY6yL09inKGicJhHay/toZivGtLxPIA8UhRchbQwzEMOmjH3vNPZpCxqXSYTm8j6a2YoOF1Rd9wTyQFF4EdKmMANRlzIWO7UhTYHqvDryUaA6D0jTysz/oClotJf3cbEOaZmaY1reJ5AHisKLkDaFGYh6lLHYqS0JauiXTciwCqncjQR9QAuX+LHue+3lfTIgHxbURw3yYUGWmTDvAwcxE4TuCeCBovAipGQwAzF81F0ummQZkvpdSGQ7UK+0PgMws+x75YuUdu5EdXQpcsxHdXQEaecuADOU0/CB5cuksUpmYrVKpreg+YjMwhRehJQIZiDGQRlFr/XYp50BaHt47eV9hr6P/o7bkY2eirRjB5KhdWip8LJ8magHzKs3wB0KL0IiwHZQYwaiLi6p8mUSvU5jXwAB81aH175IaYpkwZVI8nvGhN/VrT2+5ctE3Xqr3gB3KLwICRyXQa2MlpRYcJ14l0n0OscjKZ+89eGVMxW1hZ/Ny0T9naPeAHcovAhRxOYl7BoEXyZLSkxEOPEujBDikbysq1gw0WQqjh/b4mWi/s5Rb4A7FF6EKGH7Enad0JXJkhITEU68C0M7HimWdRWjyVQcwzZoXj1g3qEBIcThU3gRooSthSTCCV3paNeaW4DlQKUcj+RkbVQ0TYZgGbTOVNSMWfd08FDi8Cm8CCkYH+UcaMXSo51rblmfu7LqjGVdRW3LoMvNrOo693TwUMIBKLwIKRCWcygfobysNXByjSmqTqdnSfPBiyRTEVB2nXs6eCjhABRehFhSdCA8UD4LSRkJ5WWtQRCuMUuSsYgk2JTx1MpW1J5tOdzMqk31dHDt7h/HiIjOkR3o7u6WwcFB7WaQNsbWkhVKDAFpjkuQbQgBuSr09aH2uUmusS+sAzZtamkTrPpe+cGL6rnXvpk9HF/7lKbCGHO/iHRP9TdavEhbU3Rh0lBmVGRm2rnmlhPKrjHr66TsDw4iW9HD8j6F40GhRiV6x6DwIm2Lr8KkbTtIR0Q7x20B8RTxtL5O2ssQabtkY1nex8ODF+OzTOFFSgkLk5KZaOe4rZiKeFpfJ+UHVD1b0fJlpi5SPDx4MT7LFF4kGmxn6SxM2r6061qJLgThFrPEOVtRI2AeUHfJRrO8j4cHL8ZnmcKLRIHLLJ3xWO0J47bsUHeLOeKUrVggMS38bXt8b830UGXeU+H6IKDwIqoUHdwOMB6rXVF3o0SCulsMDoOq9vI+LveT8stEbXkfBsw7Q+FFvGDzYvUZ3E5LVrmwuZ/U3SgBEMPyPk6DqqKaTlOgOq+OfBSozgPStNKS4+5HDMv7MGDeGQovUji2LwGfwe20ZJUH2/up3QV3LMv7OA2qmsv7oIZ+2YQMq5DK3UjQh5YWYIhleThvbH8AACAASURBVB8GzDtD4UWsKdotyOB2YoOrQG/Xe8QpaD6GTEVAVyRmGZL6XUhkO1CvtN7s4nDjqwqV2ALmA6i2SuFF1NyC7W6hIHaUffZbFLEEzScJ0H/NTmRbh5Cu70SSNGmjVrai9o0Xy/I+4w2wPKhL5nHh5xFI8BiFV0kpuvSCL7dgO1so2h2WfiiWEILmrajVkFzUgyTPgR1VoCvQ5X20b7wkQe2ae60FqjehUuD5q+ueQILHKLwCoOg143yUXqBbkBQJSz/YE0stKet2Kg9+Tof3ceO5BMxf1NV4RnYA/V0tfgY8qCR13aNtxRyDwssRF0tSkRYnl219lF7QnvyRcqH+Ao6EWGpJObVTe3mfVDFbMZaAeU8NUNc9gQxkFF4oXiT5sDi5bOur9EI7Wx2IPSz9UBzOY5/SQ+psRdJc3kczWzGWgHlPDfB26SOrttr2wsuHSPJhcXLZ1vXmDuA+JCWBpR+KJYhaUhY4j9FaAfOAbrZibAHzlg1Q1T3qgWPutL3w8iGSfFqcaJ0iIeOahNHO96jNYKVeS8oS50zFgonG1RlbwLxFA9R1j7pP1p22F14+RJJPi1O7D1ZEB9t3tbp7JBKsBytN64wLypmKsbg6GTDvgQhfOm0vvHyJJAokUhZc3tXq7pFIsB6sYhlUlEffEFydNrMTdZHCgHkWUA0FiiRCpsf1Xc3nqTnWg1UASjaGYqPark7b2Ym6SIkpYN61ETYHVveLNqDwIqSNYQaiDk6DlaKSjabYqLKr03Z2ot1Nrg1QrzJfdEepmxwbUHgR0qYwA7F4nOKWUUOCDECKEAPmgbFxSnn9R6s+1R5QHbMVvTStYJWkahzydfBAZpEUXoS0KcxALBansSIQl0cztNd/tO4m7QHVMVuxcMoWNO/r4IHMIim8CCkZzEDUwWms0LbQWKK9/qN1NykPqOrZimULmvd48IadOVG1M1N4EVIimIGoh9NYoax6Y1n/0amblDIVgQB0dNmC5j0dPBRDM4UXISWCGYh+sCp26howH8O6isrqXDVb0aGj1K3HDtdJfXUdHxH7sQjkMSi8CIkAug/1cLUiWr/IlVSvU8A8oKvONbMVHUbpIKzHFtdJ3eLjowExCeQxKLwICRy6D3VxmiUHUJyxGdoB84DbuopqJoo0Ra3yO8hGVyGt3I2kySgdQ0kFdYuPjwbEJpBB4UVI8NB9qIv1LFndnGCHdsB8LOsq1pCgx/Qjh0HVCPpRaW0wtof7Sd3i46MBIQhkRyi8CFGEBUz1cAk1sZolq5sTLFEOmHfqJl8mCtvlffZWUBcg31sO65C6xcdDA9QF8iyg8CJECRYw1cPVmGA1Sw5AIVuJSeUbSn1dxZIu7xNNlfmCG6AukGcBhRchBWP7DnKZ0IZgHi8TXoxTIdSScslWVLqhVDMVgVIu71PKKvOWqAvkWUDhRUiBuLyDYnxhlAVvfa8oaJyzFbXQXlfRISYoFuuQqpfb18Et+0ldIM8CCi9CLLF5D7hasWJ7YZSFMva9drZiFJmKUI4J8mQdKl2V+VoNtXQTspFVSOdvQpL1NRVf1t0YQOYxhRchFti+L13fQXQhFou3UJMAXtbN0MxWjCVTEVCOCfIkOlUnEj4C5rc8ip78B8hRRTXP0b/lViRFWdECyDz2JryMMTcAOBPA0yKyfOx3/xnAzQCWAHgMwHtF5D98tYGQZhQdj1VGS0oseHunBvKybopitmIsmYpAXNahUlaZtyDDGuSooo55yCHIsKYYq2Qgmcc+LV7fAHAtgC2TfncpgH4R+ZIx5tKxny/x2AZCpsVXPBatWDp4e6dqu8ZsB19F1R9LpuL4oWOwDqnrfcUGpBsWo3pjHXleR7XagXTD4oJ2nAYRWOtNeInIncaYJa/49bsBpGPffxNABgov4gHGY5UL1XpnmkU8a0DP2jry3KBaFfRvqzQXXwo3ZyyZiuOoTo5iCJhXbkCSAP3bKsW/bwN5kbc6xut1IvILABCRXxhjfqvFxydtAOOxyoV6vTPFl3W25XHkw4ejjgry4RFkW3YjSQqa/VtiZXGLKFPRGwXHAaobZzw1wIv3Ut0n60awwfXGmI0ANgLAokWLlFtDQoDxWO2Jq2XSy3VUelmn2I4qzkEOQRUjSLEdwIaWHd/a26TtjtWuXu7BLaf+fvIRNO9qwbXeaQQxmJNotfD6N2PMYWPWrsMAPD3dhiKyGcBmAOju7pZWNZCECeOx2hf1mb8iyYaj0X/DGWNp9Xcj2dDX0uNb66l2zlScaICfbEW1gHkPDfBiwVX3ybrTauH1XQDnA/jS2Odft/j4JDB8VXmnJSt8Cl8rMTJsl/dJsj4kWQakM9cy8oG1nvJ5kSw6Sl2cOzRAtTKJdpV5HxZc9Ys/C0TEyxeAmwD8AsAIgN0APgSgE0A/gEfHPv+zzb5WrlwpJC4GBkR6exufM22zcKFIpdL4LGpbEj7tfj0HBkQWLtgrFVOXhQv2Bn3+Ns+y14Nb3iiq7bRsgPp939vbODjQ+Oztbe3xBwZkoLpGes1nZaC6prgOUL/4rwbAoEyjaXxmNZ43zZ96fB2ThIHtpIpWrPYlQu9AoYQQNG9LghoSZGgkpLf4ImUZasMrGkVhh3c0LH8zlIkI3S2nft8rB8w7WXAjC5h3IdjgehIeRbsFmVXYvgThHVD0+YQQNG916squqVrnmegZ/WSjgvlojv7OXVEv76N+3/sKmHfpJpsXeYQB8y5QeBGrl7CP4HZascqJzf2kfu2VX+yaQfNOp65sosmGupB3jC383VFBNtTVOuHl4dzV7/vxRhQZMJ95uEXUTYN+ofAqKbYzWm23IK1Y5cJlUFe99p5e7F5cLgXjdOrKJpo0BaoLzNjhTSmW9/F23ytZcNMUqM6rIx8FqvOANK0UtFNtk7g/KLwio2jrFN2CpEiimah6eLHHUmXe6dR9mWgsRYKqhcjh4OqeMcUGJKihXzYhwyqkcjcS9GHOsYBBmAb9QeHlCZfJh6Z1im5BYoPtPRrNRNXDDR1LwLzzqRctEGs11NJNY27WTUiymS1+qkHzlgdXn3BoWnCzDEn9LiSyHahXmh9b3TSoD4UXihdJLpMPbesU3YKkGa7uw2gEesE3tHbAvAuamYq1LY+iJ/9BI2A+z9G/5VYkrbxRyhg078uCa9NNLsdWNw2GQdsLLx8iyWXyEYJ1ioKKzITrZLpd7yftKvOA/bqKmoNfhjXIUUUd85BDkGFNa6VfGYPmfVhwM8tucjm2umkwDNpeePkQSS4TAFqniCY2A7X6bD4AYqgyH8u6iumGxajeWEee11GtdiDd0GJ3rKcq86ouUQ8NcI4FtDk2XyYAKLy8iCRXgUTrFNHAdqBWn80r4xQ0r/iAWusp5cEvSYD+bZXgrUPqXjHtkic+nvt2f5mM0fbCy5dIcnn/UkwRDVwMH+18j8YSNG+tp3wNfg7WGfX7yaIB6l4xzYD5MaxjAdVNg3HR9sILoEgi5aJ0GYjKxBI0nyRA/zU7kW0dQrq+E0nSNfPGipmKXii4jpX686EZMO+ysbppMD4ovAgpEaXNQFQkhKB5K2o1JBf1IMlzYEcV6GrdAFjGTEX150MzYN5lY3XTYHxQeBFSIto+A9HB6hFDlXnAoZ2KA2BsmYqqpaRiCZi33VjdNBgfFF6ERADdhxY4WD1iqTLvZMhRvPixZSqqecY8HrzwNVJtN1Y3DbqhtLLSflB4ERI4dB9a4mD1iCVg3smQo3jxY8lUBJQ9Yx4D5lXXSI3EdB5KOBqFFyGB0/buQ1scrB6xBMw7G7E8XPza5kkB+xunD9j3dt8V7BdUtQp7OrgXPReKSimQUMLRKLwIUYQFTAskSVC75l6rrL5YAuadMhU9UNu8Ez1/tBQ5lqF6R45+7JxRfBXfgJIFzXs6uNM7wlbIhqJSCiSYd6mIBP+1cuVKIaRsDAyILFwoUqk0PgcGZt62t3fmbdodl/6c+IfQO9X5pIql9/RtUsGIACIV5NJ7+raWHl96exvnDjQ+e3un3VT9cio3wOrwri8dxXvPBZeub9VlAjAo02iaGS1expgEwPsBnArgMAAvAfgRgNsBfEtEfuVbGBISGz4mlG3rPnQgy4B8WFAfNciHBVlmggyYH8fqPlG2OqTrO1G9I9/nkl3f2bJjNxqQWpko1L1inhpQeAKk60sngoBR164P4V06rfAyxvwNgKcA/DWALwJ4GsABAI4BsBbAXxtjviIi321FQwmJAZeXQDBm75KQdu5EdXQpcsxHdXQEaecuAK11zdlifZ9oL++zsQv9sIvx8tMAu8Ff3SvmoQFetFyaAvPmAaOjjc9m91MIKqUJ6td+Fsxk8fq/ROTZV/zueQAPjH39mTHmUG8tIyRCSjihDAKrVPmh76O/43Zko6ci7diBZGgdQhVe1vdJADdJsrELycaWH3ZSA5oP/uqTGA8N8CYoRPb/jBznrg+gnsS0wmtcdBljLgPwDRF5YvxvxpiNIrJ5CmFGSGnxEQgfwYRSHRfrULLgSiT5PWOdf3XL22pLmgLVeXXko0B1HpCmlek3VsxU9EbBg5+6PvXQAC9iMssaSk6k8RmwecglmdW669V90mNMF/w1/oWGi/FhAGsn/e6BZv9X5BeD64k2DITXwyG+Op7OHxiQgeoa6TWflYHqmpa2d+D6h2QhXpAKRmQhXpCB6x9q2bEbDYgnaFv7fir88JH0vbdmOr1M5gZmG1w/xpMA3g3gFmPMrSJyNQDjSQcS0lIYCB8+2tYhL2QZkvpdSGQ7UK+0dnmfrUPIsWzf8j5bh1rrSnT0oal5hmIJmHfZsbpp0A5vblZ1n3QDqzpeIvKvxpg1AL5mjLkFwEK/zSLEPwyEj4MENfSPXopMViEdvRsJvgS0diVAa6wHVc3lfSLJVASUPUOxBMyP73jt2n073rYt6pmht8cjEOFpI7wGAUBEXgbwAWPMxwGs9NoqQuaIzQDIQHhdrEXKli1IRu5EgjuBkcbPIV4Ap0FVc3mfSDIVAeWMtZgC5rdsAYaHG98PDwf7jNji/HgoLjw+G5oKLxH5yCt+/h8A/oe3FhEyR2wHQAbC6xFKjGuROA+qijeUl0xFD4OfqqU5loD5yCh4BagoXyYz1fHaCWDafFMRebOXFhEyDUXHY9GKpYeTSNmwAbjhBmBkBJg/v/FzgDjFonlCLVvR0+Dn7RktfPS3w9v5RPKMeLlNIizkNZPF68yxT4NGpfoz/DeHkKnxFY9FK5YOzgHzWRa8Qk5QQ79sQoZVSOVuJOhDK2PRVNdV9BgwX/gzqhw07+WdE8kz4kUjRWhGnKmO1+Pj3xtjhif/TEhR+MoqpCVLD6vkKleREoNCVsxUBJSzFWMJmAfiC5qPKHapGV40UoQvfKusRkJ84DOrMIJ3UCmxvqbKIsUVq/FPeeatmq0YS8A8EE/QvEumojLeip2WSHROZqYYrxWTflxojDkRk+p3icgDPhtG4oZZhe2J9TWNyD1gLSaVb1Jv2YoF+9DUL30sQfORZCp6WaRa3Szql5ksXn826ftfAvjKpJ8FwNu8tIgEi+37l1mF5cP22ltf04iUdJYB+bCgPmqQDwuyzMxcJqLgc3Ga+BedrehhAFQPmB9vRAxB8xHgxdqnbhb1y0wxXmtb2RASNi7vX2YVlguXa+90TSNR0mnnTlRHlyLHfFRHR5B27kKrFt+u1YCetXXkuUG1KujfVgk2HqqdA+YBD+cUSaaiF2ufulnUL01jvIwxFQDrACyZvL2IfGW6/yFxUbRbkFmF5cJ18pmghgQZgBStzOrzRTL0ffR33I5s9FSkHTuQDK1Dq4RXtuVx5MOHo44K8uERZFt2I0kWt+TYAKwfZnXPUEwB87ZEkqnoPIEu0dJGs8UmuP57AF4GsBPAqN/mkKLQdAuW/JlpO5wmn+qjlT0uy/skC65Ekt8z1gFXt6iFQIrtqOKcfQHz2A6ghZYPy4dZ3TMUS8A8EFXQuGqx0xLPym2E1xEslhoONg9CCG7BEj8zpaLwyafLaKW2+nFEy/tsOBr9N5yBbGQV0vl3I9nQ17Jj72tE84dZ3TMUS8B8ZBMTFjv1g43w+htjzOkicof31pQIlzGlaOsU3YLEBi+TT9sbSnkAcgqYB/Ru/iRBkvUhyTIg7SuuDQWLXq/atExV5rOskaE4Otr4DFh4eNFI6go9DGyE1z0AbjPGdKCxRK0BICJyiNeWtZCiRZLLmOLDOkW3ILHBy4vV9oZSnvlqBsyPY/3eKVr0eVzep/BLqCzQCz+nzs6G6AIan50trLXmCIud+sNGeP0ZGhGyO0Vk2rUbY8WHSHIZU3xYp+gWJDZ4m3za3FDKM1/NgHlAOVvRUfQqeoS9CXS1cxoaAjo6GqKro6Pxc4thsVN9bITXowB+VEbRBfgRSS5jii/rFO/t9sXLi7VotGe+igHzgHK2osNLRz0kyYNAV13eJ02BBQvUJhwsdhoGNsLrFwAyY8zfABge/2VZykn4EEmuAonWKVIUXl6svvBwcCf3naLwU81WdDh39VhoD9fJyznZPnjK952Xc1e/SeLDRnj9fOyrOvZVKnyJJJcxhWKKFIW3d6Cqv8kOZ/ed4oPnLVuxbMv7AIVfJy/n5PLgKd53LHYaBiYGD2J3d7cMDg5qN4MQNYrOfHU+eASuhL4LH8dl1x2OOuahghFc+dGnsOlrLSw26krRYjaA6u26O1U8fAkXtNbfadwYY+4Xke6p/jbTItmbAfy/IrJzir8dBOBcAMMi8u3CWkoIeRXeluyxJRJXgnqxUTiOP0VbPhyuE5f38XALjxsxlIwZXqpuKJXyKDszuRr/HMBlxpguAD8C8AyAAwAcDeAQADcAoOgixDOuuqfwd2AkrgTtYqPq6ypyeR+dgHmgsU293hBd9XrLJyftbOmOkZkWyX4QwHuNMa8B0A3gMAAvAXhERH7SovYRUmpsE6FUdY92BiIsx7/EU7FRS9TXVbS8TuoGzFiW93ERHsoPKYPm46JpcL2IPA8g898UQtqLSBKh9jVC6aXrZEnSDFwOwNVpc/5lFPJBBMxrZsmmDJqPCZusRkKIByJJhHKn4EBbdUuSJeqZipZ40wiKQVZezslVeGhmyfo4/yBmfOWEWY2EFIxqBqI2Hk6qduEW9Fx3TmN5H4yg/6O3Ivlaiy1JtkSSqVg4HtupmjBXtuxLbzslr2RWWY2T/nm5iPyo+GYRUj7UMxBdKfol7CEuRDto3gnFTEVAcUz1FA+kHjSvCAPmy4uNq/E6Y0wVwDcA/H8i8pzfJhESL+oZiC74eAk7LkdjNfYpB80DQG3zTmRbh5Cu70SysYWLaceyvI+neCDVoHllkcKA+fJiE1z/O8aYowF8EMCgMeYfANwoIj/03jpCAiKKDEQXfLyELc14MVWZr23eiZ4/Woocy1C9I0c/drZOfDmYRVXHVE/mW9WgeWWRwoD58mIVXC8ijxpjPgdgEMBXAZxojDEAPisi/9tnAwkJgagyEG3x9RK2EEmxBMwDQLZ1CDmWoY55yCHItg4h2VjAjmNZ3kexiKZq0LynDvW2iL3NjqN6QZUXmxivNwP4AIB1AH4I4J0i8oAx5rcB1ABQeJHSU8oMRMWXcBClFyxJ13eieke+r63rO+e+Uw9uLC+Xs4xV5m07ykOHelvE3jW4NIoXVHmxsXhdC+DraFi3Xhr/pYg8NWYFc8YY8ykAHwYgAHYC+ICIvDybfREyF2wHgNJa6JVewjEFzCcbu9CPgmO8HJS8qkjx4G4LImDetqMK7lBv3kvGbkWFjfA6A8BLIlIHAGNMB4ADRORFEflfrgc0xhwO4L8AOFZEXjLG/BWAP0AjeJ+QlhFdBmIkxFBl3pVkY1cx7sVxYlnep4xV5hXxNoEr7cywnNgIr78DcBqA58d+PhDAHQDeOsfjLjTGjIzt76k57IuQV2Ez+EeVgRgJsVSZBxQzFQFrJa9uyPAw41CvMg+olZPwNoHjzDAqbITXAWPLBgFoLCFkjDlwtgcUkSeNMV8G8K9orP14h4jcMdv9EfJKbCe/pZ0kKtYoiiVo3lumYsHuLvWAect2uqAaMA+ox605daeqn5n4wkZ4vWCMWSEiDwCAMWYlGoJpVhhj/hOAdwM4EsBzAG4xxrxfRL71iu02AtgIAIsWLZrt4UiJsH0H2U5+SzlJVHa5xBI07yVTkQHz4QfMA/HErUXiPiXu2Aivi9AQR+PuwMMAnDuHY54G4Oci8gwAGGP+Nxpuy/2El4hsBrAZaCwZNIfjkRLg8g5ymfyWbpLoyzdlOarGEjTvJVPRU9+3dcC8C7YdFUvcmrqfuaQEsGqBTQHV+4wxbwLwRgAGwI9FZGQOx/xXAG8Zc1e+BKAHjfpghEyLazmH0lmybPHhm6rVUFtz6ZiYuhTJ9i/N2PkxBM17yVR07Hu1938swmOcojsqlri10sZCKFKroZZuGnuXbUKS6byjrBbJNsa8FcASTBJqIrJl1gc15r+hYTXbC+AfAXxYRIan256LZJeXtl5Q2hcFD1S1s65Cz3c+gRxVVJGj/z3XIrntM3Pebymx7Hv1+7noe8RniYhIHnwuaB0+tQu3oOe6c/a9yz56K5Kv+QmFmOsi2f8LwFIADwKoj/1aAMxaeInI5QAun+3/k3LAcg6eKLr20FPHIEd1XzzUU8cg1O5XzVQErPvem4VIqcq8t+dT2d3mJW6NAfNqZFiz/7sMa1TeZTYxXt1o1NxinBWxhuUcykP6oaWo/sOkeKgPLdVu0pR4XVOxYMuDFy9STNl6trhmK4ZuxYvIghcVltc+3bAY1RvryPM6qtUOpBt0Mq5thNePALwewC88t4UETtFuQYYw6GO1vJuPeCgPeF1TMYZsxTIGzdt2lIeGMmA+EhzitpIE6N9WUfec2AivQwE8bIz5BwATcVgi8i5vrSLB4fJea+tyDhHhUuy08MrtHvCSqQjEs7xPTEHzRXeUh4YyYD4OalseRU/+g0bcVp6jf8utSGbKvA7Ac2IjvK7w3QiiS9FuwbYu5+ALD0G2sRQ7tcWbZS6W5X1iydbz0VEeMkq9TAw52yycUOK2nBCRpl8AFgM4bez7AwEcbPN/RX2tXLlSiBsDAyK9vY3PZtstXChSqTQ+p9vedjvX4xMLXDvfdrcf/aYsxAtSQS4L8YIMfPSbhey3lFjc0L29jUsEND57e1t3bJ8UfnhfHWXZUE+Pk/p1Kh0u13PBXqmYvbJwwd5guh/AoEynqab7w8QGwEcA3Adg19jPRwPob/Z/RX5ReLnh8mJxeQfyvaKE40BlfZ0GBmSgukZ6zWdloLom6As7cP1D0nv6Nhm4/qECd1rsDe1lQPclujWfZW/Kxw4vuk/5nEqH47spxLFpJuFl42r8OICTAdw7ZiF71BjzW8Xa3YgtdAu2IQ4XycmLE0mxUy/ZigyYL5VL1IUgFukmMxJj3JYLNsJrWERyYwwAwBgzD406XmQGXEJyNLMFGXIQAQ4XKcuAfFhQHzXIhwVZZma+phG8sbxkK3oaKBkwnxWfWVBwfKOXdx6D5gslyrgtB2yE13ZjzGcBLDTG/C6AjwH4nt9mtZaiRZLLjDKEbMEIxl5ieZHSzp2oji5Fjvmojo4g7dwFYAbrUASVsb1kK8ayvE87B8w77jeaRbrbHYsLFUq9LW9M54Mc/wLQgUac1y0Abh373jT7vyK/fMZ4ubjmbbd1iSFwjbFiGEF7Yh3D0NsrAx2rpBeXykDHqtLcUJoxXhF1kzXRBMxb7pcB85HgELsVe9djLjFeIjIK4OtjX6XDxexuu63LjJJuQdIMJ2NCmiJZcCWS/J6xG+rq6XccUVyKlzpilmaPIFxzBe8yBpeoy369XCP1YLjy4RK7VWZPjM1ajT/HFDFdInKUlxa1GB8iyUUg0S3Y5liMlE6DissNpRyX4mVdRQ9iJhbXnKpO8DUrtNwvA+bjoOyxW9ZMZwob/wLQOenrcAAXAfh8s/8r8st3OQkXk2bs5k8SEJb+Ea91apRu6IHrHxqrITbSqCFWhAvRo08wBtdcWeuI2VJ4M8voY/aEbShAqDW3fIA5uhqHXvGra4wxdwH4r8VKQD1crEi0OJHCyDLUhlcgGz0V6fCORlmHKW6uBDX0yyZkWIVU7kaCPqCoeaLSDR1TpiIQh2suFsvcbJqg4j5lbIcVLuVekiSMtRK1sXE1rpj0YweAbgAHe2sRIW1CrfNM9Ix+shHvMJqjv3PX1HIqy5DU70Ii24F6pRQuj7bOVAS8DOqx1BEDYN353nSfmporH66TKHapXTmJP5v0/V4AjwF4r5fWENJGZENdyDvGam51VJANdU0tvEpYI8jLuooOykN9QB9vb8EjUAyWOZfOZ9B8+HhbnL7E2Lga17aiIYSUBduxN02B6gIz9v43049pJXV5lC5T0dOAXjbLnEvnM2heF5sEGG+L05cYG1fjn8z0dxH5SnHNISRuXMbeJAH6r5n0wkpmeGFFZJ+PIVsxlgE9CONM0feeQ+ezyrweTrFbPiZRJcbG1dgN4CQA3x37+Z0A7gTwhK9GERIrWeawZE+thuSiHiR5DuyoAl3xuzzael3Fsi7vUzRJgto199pNOMCgeS28JMAQAHbC61AAK0RkDwAYY64AcIuIfNhnwwiJEacle0ro8ogpWzGGAb2MmYq1GtBzUVfj8DuA/q5yxNeVDcZu+cNGeC0CkE/6OQewxEtrCAkYm/d6MvR99Hfc3igR0bEDydA6TCu8fLk8FK0Z2tmK3k5dKQsuqkxFzcMH4ZONA9tQAMZueWS6Al/jXwD+FMA/AbgCwOUAHgTw2Wb/V+SX7wKqhDTDupaia9HFoqs+BlD0UWtdRa/r9XnYsVpdUl8dpbn2pbfqseXCS+FiMiWYYwHVLxpj/gbAqWO/+oCI/KMXFUhIoFjP0rXXgPJkzXAJmNfKVvRmyClb0LwPM1qthlq6S+vw2wAAIABJREFUCdnIKqTzNyHJ+qbdbyzxdWWEcVthYONqBIADAfxaRG40xvymMeZIEfm5z4YR0gqcSj/MqyMfBarzgDStTL+xZvyIhwHIS8A8EEemoqcdqwfNF3yPuix+7OHwDJi3hHFbYWBTTuJyNDIb3wjgRgDzAXwLwCq/TSPEL06lH1yW7VHOGCt6APIyS44lU9HTjssWNO918WNWmbeCNbfiwcbidRaAEwE8AAAi8pQxhksGkehxsjpkmd2yPSEE+RY8AHmZJceSqehpx2ULmk83LEb1xjryvI5qtQPphsXF7DiE5ykCWHMrLmyEVy4iYowRADDGHOS5TYTMGZtJspPVwXbjEpaI8DJLDmFdRQ87Va1ooBjnlCSeFj8u4fPkA8ZuxYWN8PorY8z1AF5rjPkIgA8C+LrfZhEye2wnyU5WB9uNYwvytVQKhc+SHTrfi9HDw07VjTOefK2qnr7YniclGLsVFzZZjV82xvwugF+jEef1X0Xkh95bRsgscZkkOw0WNhvHFOTrkInmBcvO92L08LDTIIwzRQfN14CetXXkuUG1KujfVimFmIwF1twqJzMKL2NMBcDfishpACi2iCpOGYi2k2QfPqxIgnxdM9HsdxxBtqKHnXozzigma2RbHkc+fDjqqCAfHkG2ZTeSpID4LVaZb4prNjFjt+JhRuElInVjzIvGmN8QkV+1qlGEvBLnxadtJsnqviFdvGSixZKt6GGnXtqpfI+m2I4qztnnwsJ2ABvmttM2f+5sYdxWebGJ8XoZwE5jzA8BvDD+SxH5L95aRdoG24mvqxvHapIchG9IDy+ZaA596sXo4cmSohbnpHyPJhuORv8NZ4y5o+9GsqFv7jtt8+fO2nLPuK3SYiO8bh/7IqRQXCa+QbibNOtzuRzbNmDeRyaaZZ/GEjDvcbd2ePJfWt9OSYIk60OSZUBaUAxgGwfMO1nuGbdVWqYVXsaYRSLyryLyzVY2iLQPrkHwqu4mzdHX5diOAfOFW2gs+zSWgHmPu7XDw43vfCsXfZO0ccB8lgH5sKA+apAPC7LMzGxBZdxWKZnJ4vUdACsAwBizVUTWt6ZJpAwUXkcLnmJsbXeqOfq6uO98Bcy7YNGnQVgwdXer5r/0diszYL4paedOVEeXIsd8VEdHkHbuAkBLVrsxk/Ayk74/yndDSHnwUkdrNo0ocsea7hGHY3tbuqXg/lS3YGrvVtGC6rTuqC0MmLd6RJKh76O/43Zko6ci7diBZGgdKLzaj5mEl0zzPWlTfATCe5n4+hgENN0jDsf2EjDvaVB1uvbK6/WVKWjead1RWxgwb/eIpCmSBVciye8Zm0Rd3fK2En1mEl7HG2N+jYbla+HY9xj7WUTkEO+tIy3BZkxTD4R3wdcgoOgeqe18DbIsQdr5mqa6o/CAecf+LDwHwWPQvFqYkWbQfJbZrTvqgvpDr4t17FYbx7eRfUwrvESkAPszCR3bMU09EN6Fkg0CzoUUi9aHDv3pRSN5ENLqnjHNoHkfz4f6Q6+LU+xWm8a3kX3YlJMgEVK0WzCIQHhbSjYIqBdSdOhPL8ZGD0KhjAHm1ufk+nwou3k1sT51xm4RByi8SogPt2B0WqZEg0AQhRQt+9OLsdHDzeelndpV5lOHc7J9PtRNg3o4nTpjt4gDFF6RYTMD8+UWLJGW2YdqoI8dMRVSdNJIitYhLxMJ5QDzMp6TJtZxW0CEM1OiCYVXANiOP7YzsKjcgi4ULZICmM3XNtsJKm+FFD0IT6v7KYCAec1YOBeczgk1JMgApJhzpiJQunhJF5xrbkXzIiXaUHg54iKSihRTgP3ks5STLx8DtfJs3jVovvgG2Pdp4fqMAfNWOJ1T2cqoeIQ1t4gmFF4oXiT5EFOA2+SzdJMvHyJJeTavHjRv2adeBE1MAfMuaAXMu27cxlXmWXOLaNP2wsuHSPIppko4+bSjhCnw6QnP7R80f8JzLT2+bZ96ETSxBMyPoxQL6HROthurmwZ1sY7dausXLvFJ2wsvHyLJp5gq2eTTHl8vQQ8dap2C/tpH0G+uRiarkZo7kbz2TADvKbQtM2LZp94ETQwB84CqUHE6J9uNgzAN6sGaW0SbthdePkQSxZQnIuioWg3oWVtHnhtUq4L+bZWZU9APuBJJfu/YzffllrYVgFWfOgsaD9Yh1VJSnoSKl3Oy2bikAfOsuUViwYiEvwxjd3e3DA4Oett/BBUFSCT0Xfg4LrvucNQxDxWM4MqPPoVNX5thvcSy3XwerEPqnrGynlOJ7jv1JARCXoEx5n4R6Z7qb21v8QKiMKSQSEixHVWcsy9uC9sBbJj+H3zcfA6DagzZiuqeMQ8+TKdzcrlIbVpl3qk/GbtFlKHwIsQSqxT0DUej/4YzkI2sQjr/biQb+lrZROcSETFkK3rzjClm9lmfk8tFamNLTpoC1Xl15KNAdR6Qpk2WGi6Z8CRxQeFFdPHl8ih4v9axW0mCJOtDkmVA2tf6l7vD1D+WbEUvBgplkWJ9Ti4XSd006AerCQ9q6JdNyLAKqdyNBH1AEQVkCfEAhRfRw9fg52G/2ZbHkQ8fjjoqyIdHkG3ZjSSZJnZLczbtYB5ysiQp130qfJfKAfOA5Tm5XKQSBs1bP8pZhqR+FxLZDtQrpRGdpJxQeBE9fM3QHfZrO1A6x25p4WAest40gOV9CseDSPHSTS7mvhLGLlk/yiUUnaS8UHgRPXy9LC336zJQqsduueBgHrLa1INAVg9H0g6Y90XJYpesY7dKKDpJeaHwInr4ella7tc1E0o1dsuBwi1JXN7HCi/zCHWF6gfr5EuX2K2SiU5SXlSElzHmtQD+AsByAALggyJS02gLUcbXy9Jiv94yoXz40Cz3qe7usiSITMWCce4mm7YGoVCLxekeZewWKSFaFq//DuD/iMg5xpgqgAOV2kFKilomlA/l47BP53Faqe5TbJmKhXeTbVtLGLvkdI+W8PwJabnwMsYcAmA1gAsAQERyAHmr20HKi2omlA8LhcM+nTMVFYPmY8pULLybbNtawtglp3u0hOdPdAlh0QYNi9dRAJ4BcKMx5ngA9wP4pIi8oNAW4gvFu1s1E0p5n07jVNmC5j1ZR7zouTQF5s0DRkcbn80WiY1AcLgI7v5rdiLbOoR0fSeSpMk6iZGcPwmfUEImNYTXPAArAPyxiNxrjPnvAC4FcNnkjYwxGwFsBIBFixa1vJFkDihbUqzHXx+z6QD2aT1OlS1o3pN1xJu3a3yd3AjWy22G61qJyUU9SPIc2FEFusqRMEDCJ5SQSQ3htRvAbhG5d+znW9EQXvshIpsBbAYai2S3rnlkzviypNhUjofjjNrHbNrDPmtIkCFBijlHoe2jjEHzHvrei57LssbzIdL4jDxo3OmRD2X0I21HKCGDLRdeIvJLY8wTxpg3ishPAPQAeLjV7SAe8WFJcakcX7IZtbMBUbHKfExB86rF+EMZAQrC6XRKdu5EHyc3dwAhg1pZjX8M4NtjGY3/AuADSu0gPnC4u71Ujo9pRm3RAU6nE0AQQwxB8966KbYRwAKrDGFXK3Mk507Cx/VZDiFkUEV4iciDALo1jk1ahMXd7a1yfCwzassOcDodj5l9auNkLLFoMY4ATbA+JVcrcwTnTuIgpnn2OKxcT+wpePR1emASh8rxscyoLTvA6XRiWYPQhVhi0bIMGB5uZCoOD8cxAjTB+hmNcfQjQVN4MlVAUHgROxxHX5uHxusDE8OMOk1Rq/wOstFVSCt3I2lSJsLqdDyIlCDG1Bhi0To7G6ILaHx2dhawU12sn9EYRz8SLE7ekEjm2ZOh8CJ2OIy+tg+N0wOjbnYpnhoS9Jh+5DCoGkE/KtNnLCpGgqtnKnraZeHafGgI6OhoiK6OjsbPkWMduxXj6EeCxXWyF8M8ezIUXsQOB+uMy0Nj/cCEYHbx4WrdW0FdgHzvDKekLDpjyVRU1+ZpCixYEIXVx/pWdondim30I8FSdgMqhRexwsU64+Wh0X4SPYzq1qcUgOiMIVNRvZsisfq4LlKtfe+RcmGdJRv+ozRrKLyI1YNgbZ3BLFyIMaTfuwxAludkfUoeRadatqKHc/KqzTdvBrZuBdavBzZunH67CKw+TlpKe8JDSoVr7JaXRymExRpFJPivlStXCvHDwIDIwoUilUrjc2Bgbtt5OXgIqHbU2H57ewvtI/Xu93ROBe9S5PrrRRo15htf119f4M5bj/N199KppB3p7W3cd0Djs7e3xQ1o4UsPwKBMo2lo8SoptqLedvbrxeAUkxvDtgOyDLXhFchGT0U6vKNR/qKIc/Iw/fPW/R6W91FcMahh6XrlzzNZvRQpvNjp+D+E+lySIIim9EMoY850iiykL1q83HAR9apWj4EBkQULRIxpfJZgRj1w/UOyEC9IBSOyEC/IwPUPNfkHPWtCLFZMdctcJBYv635S71BSJlxvJy+vPNud0uJFZoN1PJZDVqHT7LdoRPb/jJxsqAt5h6A+apB3VJANdc1cIqJs2YplDJoft27ZxHgpwmKnRAP10g8RFv2i8AoAWzOt7f3lZM7VXFA6yxpPq0jjswQDQJoC1QVmrO9N8Mv7FP4SjClo3qWjNm4MVnCNw2KnRAP120ld+blD4eWJosUU4CkeS3P2q/7EulHbPMkyuHFqy6BT33N5H61dBtBR9jgl/rLYKSkQ1dIP0QSOzYLpfJAhffmO8XLxOdts6+JGdsnyiCUmZ2K/tj73CDKmnGO3rHdc7Pl7yxqKJQPRlt5ekY6ORkd1dCikV9nh9HgydosUiHr8r3rg2NwAY7ymx2Xia7utixHJRax7mVn42GkQxVqKJds6hBzLUMc85BBkW4eQTOd5KtvyPmWsMh/JuopOBmnGbpECUb2dInQfutD2wsvl+tpu61NMebm/it5pCQeAdH0nqnfkyCGoYgTp+mkGagbMa+3SjUjWVXQS0jG6XEjLUffg2TSg5Pdy2wsvl+tru20QYkqTmB4a2yrzG7vQj+YxXvqKggHz1qOK8rqKhcfPMHaLNEE9AdC2ASW/l41EkMbf3d0tg4OD3vbv4hkKYbWBKIiho2o11NJNyEZWIZ1/N5Ksb+5t9WTxUu9ODw0ofJeufa/YqU5NVb/4pCz09QGXXdaYF1YqwJVXAps2tVMDWocx5n4R6Z7qb21v8QLcLASls061MbUtj6In/wFyVFHNc/RvuRXJXC+uh5maV+9lmarMRxQXYt1U9WA4UiaCsDTH4g3xCIUXKZ5IBovsl29Cjuq+gPlfvmn6YqcuFDyge/Neli1oPqKXunVTA3BdkzhQK/2g7r+MDwovUjyRDBbp63+MKpbvC5h//Y8BnDz9Pyi5fLzpibIFzScJcM01+yrMB3jPjWM9/kQkJokeqonkEVmaQ4HCi9gTkTnZava34Wj033DGvhivDX0z71DJlONtkhhT0LwNtRpw0UWNg+/YAXR1qcRuFZpUQwsBscCrVbzNMxB9QOFF7IjInGzd1CRBkvUhyTIgbRJYr2zF81ZGJIYq80Dxi5R6wFmbqwXDkVhQnesyA9EbFF7EjojMyVkG5MNjC1UPC7LMTN8U23Z6mtWpJ6x5uE6F79LLIqXF4/SIRBIHSfRQn+u63NCcHDhB4UXsiMicnHbuRHV0KXLMR3V0BGnnLgDT1N2yJaZsxRhKP7hgOwAoz7ydHpFI4iCJHt7muhGFjJQVCi9iRwDmZGvPzND30d9xO7LRU5F27EAytA5zFl5AHNmKZctUBNyXgvDQuMIzxjiokSaoL/8VwDu/rFB4EXsUzclOg3+aIllwJZL8nrE31tUtbastXl6sZctUBNQHAKexCjUkyACkwEzFSTiotS0uoX3qy3/RhegFCi8SBU7vC5c3lqIPzcuLNaZMRfUANzus7z1X0yAHtbbD6y3CDMRooPAqKz4GNcWBMk2B6rw68lGgOg9I08rcd6ruQ/Mw9saSqejS98rXicVOSVF4LfvADMRooPAqIz7WrFNegzBBDf2yCRlWIZW7kaAP07pybNsa00BZeIEot90WLhBd+t7TdSrc5UNrQlujanBiBmJUUHiVEZeHUFGkOGm5LENSvwuJbAfqlWIG6lhKRHgUvVEs7+PhOnlx+dCa0LZ4MzgxA7GUUHiVEZeHUFGkOGk5HwN1LCUiPFl8VA1+Ln3v4To5nzuLnZIZ8GJwYgZiaaHwKiMuD6EnkVK42d3XQF3wQOlFzHiazaoHzbv0fcHXyencA4gFJGETRIYyRX80GBHRbkNTuru7ZXBwULsZ5aVg35hT3PTmnci2DiFd34lkYwG1tpSJqSiql91GJFKsz72vD7jsssYAWKkAV14JbNrUolYSTVyeD+ttbTeM6Fkir8YYc7+IdE/1N1q8iJ7Vp1ZDclEPkjwHdlSBrvhfLN5iODwEzDvu1g7X+MKCxaSXc2f8TFviJQ6Q7kMCCi/iQOFxnj6DjJTrcxUew2GJ+iTZ9uKHcO6qlSxJ6Hh5PdF9SEDhRQCrAcjLRM2XJUG59IU1Ht7sqgHzgP3F1z53FjslTXB+PbGAKbGEwqusFBxH4DxRs1k6xZclwcOg7kXLlbXKvI1I0T53dYVKNCl87U0WMCUOUHiVEReVkGWoDa9oLCg9vAPJNAOQtywwH5YE7dIXii4sL+91H6rTQ5as0y5peWhbvLyeXF4QtJ62PRReZcThJVDrPBM9o59Ejiqqozn6O3dNaZ9yGtS0rQke1If1OB2AC0s1YN4Fy4Z6GShpeWhbvEyiKOSJAxReZSRNUav8DrLRVUgrdyOZ4SWQDXUh7xDURw3yjgqyoa7pHINxZYEVrD6sx2lPIkV1PWnl6+lloARoeSgZhWskZiAST1B4lZAaEvSYfuQwqBpBPyrTiqk0BaoLzNi7xRQzpvp6CcVQfiCA5W0Kx/V6FnydvFkbSWnwopGYgUg8QeEVGTYFR7MMyPdWUBcg39s83MDLRK3ol1CtBqxdu+/Num1bmAHzHjpU23MLwP56euhUbWsjCR8vyT8hWO5JKaHwioja5p3o+aOlyLEM1Tty9GPnlOIrTYHqvDryUaA6D0jTyoz7tdZImv6uLVuA4eHG98PDjZ9bXXpBab0+r+//oq+pS3xhBNZGEgdekn/oPiSeoPAKANsBKNs6hBzLUMc85BBkW4eQbHz1dglq6JdNyLAKqdyNBH2Ydlbn0siSuXFiWa/P2/vfxzlZdmos1kaiT+EZrcxAJMpQeClTqwE9a+vIc4NqVdC/rTLtc56u70T1jhw5BFWMIF3fOfWGWYakfhcS2Q7UK8W4XLTdOBs2ADfcAIyMAPPnN36eIyFkaioZ0Rr4OCfLTo3F2kh0cYrdsnEfArSMEnUovDxhbcXa8jjy4cNRRwX58AiyLbuRJIun3DbZ2IV+WCwq7ePFov2ySpJGhxZszdDM1FQ3Ivq6phadGou1kehiLdCZgUgigsILxa9A7/IOSLEdVZyzz4qF7QCmt+YkXc8jGaoBXen0jfTxYgnhZeVgzSg8HK2MQfOK1zQEayPRQ33dV1pGiSYiEvzXypUrxRcDAyILF4pUKo3PgYG5b9vbK1LpGBWg8dnbO3MDBqprpNd8Vgaqa4ppQJvj1E0DA40LptCX3i6nj3Ny2Gfhh+d9XypcL+fA9Q9J7+nbZOD6h4rbKSGeATAo02iatrd4uUyUbLdNO3eiOroUOeajOjqCtHMXgGncgkmCJOtrLNWT9jXxSzo0to3x4p7wgBeDk49zctin0+EVl1Yieji9xmo1JBf1IMlzYEcV6GIGIomfthdeLrEmttsmQ99Hf8ftjfUPO3YgGVqHaYUXYG/29hWTo1kmwsOxvbknPFC4x8PHOTns05vopWuoNDiVu2EGIikhbS+8XCZK1tumKZIFVyLJ7xkb+a9ufWNt0bT6eDq2y3XyIWRLt7yPwz5jEr2kWKwNmC7lbrSTegjxQNsLL8BtomS1bUxmb80B0PHYTgU3bVLLPVwn9QQ8T4kVtWvu3ZdNm0xvvdUWvUQHp/s+cyh3E9O7lBBLKLx84cPsrVjw0gtpCsybB4yONj5nOLZz7JBLanmB18mbjlVc/LlWA3ou6mp05w6gv6sAryAH1FKRZUA+LKiPGuTDgiwzxYluuhBJyaDwigkfo7r2ACiy/+c0OJ26ohXPi45VNqM5dyeLnZYKm8vpmlBE0U3aGQqvmFAseOmFLGuM5iKNzxlGdKdTV7TieRlTfApJi1GVxU7bF+tlDX0lFBFSQtSElzGmAmAQwJMicqZWO6KibDNFhxHd6dSV+6nwMcVnNqvFqMpip+2LtQvRV0IRISVE0+L1SQCPADhEsQ3uqKasoVwzRYegbcAyYH7SvovuJ7VL70tIZhlqwysaVorhHY1actPs27o7GTRfKqxdiGWbFBLiERXhZYw5AsA6AF8E8CcabZgVdKPYY6FSnIK2lfve2+EV46FqnWeiZ/STyFFFdTRHf+euZnK2ORyAo8D6tnNxIZZpUkiIR7QsXtcA+AyAg5WOPzticqNoF0W1UCmxBMx7O7x20PxQF/KOMTdSRwXZUNf0wksxq5IUS60G9KytI88NqlVB/7ZK62sSEtLGtFx4GWPOBPC0iNxvjEln2G4jgI0AsGjRoha1rgmxuFG0LXOWKiWWgHlvh1cWk2kKVBeYsXMy05+T9v1ECiXb8jjy4cNRRwX58AiyLbuRJIun3pgWTEIKR8PitQrAu4wxZwA4AMAhxphvicj7J28kIpsBbAaA7u7umWsNtApfL6GirVPaljlLlRJbwHzhh1cWk9bnpH0/EWusSj9gO6o4BzkEVYwgxXYAG6bfKS2YhBSKkSb1k7wevGHx+nSzrMbu7m4ZHBxsTaNajfKixr6obd65L2h+4wxB89rJCtp4OH/rvrfeof79RJpj7UKs1VBLNyEbWYV0/t1Isj5eT0IKxhhzv4h0T/U31vHSpoRFUa2D5gMY0L3oPs0q85t3ouePliLHMlTvyNGPnXMXX3Q3RYG1CzFJkGR9jSzWlKKLkFajKrxEJAOQabZBnbIVRYWDllR2YXnRfdoB81uHkGMZ6piHHIJs6xCSjU3ayyrzpcDJhcjrSYgaHdoNaHvGrQlXXlkaF06aAtV5dVRMHdV59em15LjorFRUYpym0n1h7tSedH0nqshRwUhj8F3fOf3G4yLxsssan7Va6xpKrKlt3om+t2eobd4543bJhqPRX/3/27v/WMnOu77jn++d3VkDSUlzcdpiO2sSWSWGxTa5pTlau55kI2Rs0zhsUONCFxXEgkRFVhBFuaFVU1ZoW7lCW6Eg1SUJ2RYIKOsY17Go3cHHubbGTtbGYR0cJzHExUoauxvSxED37N777R9zrnx3MzP3eXbnmeecmfdLiu69c+fHM3O8OZ/7PN/zfW7WYXu/+t2bVRy4akYjBBCDpcYmaMlfn6G1Q4UG6vuqSu1Vzx9RoSMa2fA08xJWksnG3AXzB/eor8AaL4rmGy9q6ZglRKAVCF5N0IIC86gTQFmqOLumwh+Szi5NPqFnDJ1Jcl+iMBlVNnZwz+TlxU1taY+ywKKXjlvyRxywyAheucXWBGUKaVEngOVlaWNj+P3GxvDnhoo6T2Wqh4ouG4sZJ0Xz2YTMIPf2L6t7f/Vy3dakpWMArUDwyi1muSdj4XbUCeDUKWlpaRi6lpaGP8/Y1PNpxs8+eKPiCxknMyRZhM4gRy0dA2gFglduMcs9iWpyQkJK1Amg15N27cq2hJUkI2WshwreqDjzOBEuZgY5eOkYQCsQvHKLWe5JUJMTE1KCTwCZl7CSZI9E9VBBoTdmo2LqtvIKnGplCRFYXASvVFI00UwQaKKXsVqwUXKS7JHgsw8OvTEbFVO3lc85HeFXJ3aEZwkRWFwErxRS1gMFBprQjBS8jNWALvOhorJHxjAZPDNXFNLRo9Lx49L+/Y0OvYtscOwL2lfdp0pddatK/WMfUzHhOLCECCwmglcKLerIHryM1bLaoaDskTlMBs/MDQbSoUPDO66tSXvG7cGEZAICeqkbVan7ct2WbhzVvQ7AgqNzfQpt6sje66nY9YRWO3eo2PXE+LFmfk/S8Nx35MgUG6xn7jJfFFL/6Ekd3leqf/Tk+CyVeZwLr15CPPIrL2nQWx37H2DvwG51d9lwx4ZdS+odGLFPIoCFx4xXCrF1NhHLXSF33dyyp9qQujukXq9z8WNtwsbb056cylgwv3nH4tA+FVUlrXWlPWPeFAXzWYUuIRaF1H+wQ3kdgIkIXlKapqShdTYRiSL0rsFb9sSONaMkK505C+al8DdFwXwagf/uY5YQW/BPCUBmBK/cReNlqcHpHxzWWJ1eG+6zNq6GpAwMH2WpYv3h4ZY9653ppJS21ENtytRlPiog9nrSjh3DRrM7dkx+U5zRpyviv+fegd3qfnhdVbWubpclRAAXh+CVuxB++Vbt23jXcBljo1J/+dmxf00Hh48US1OZP6foKxVzdfjvRX707ud+xWxE/MHDEiKAaSJ4parzCdiHTZLKU3tULdV9tJY6Kk/tmbiMERQ+UtSYNaDOKHjSJ2NIjProy3I4Rvfh14ZfKdoKgTOdMX/wSEw4ApgegleKOp/AfdikOs/ssjrP2LZ5ptBAhUpJPU2lbiu4cCxNnVGSPb8zF80Hn6QbEGbnSsRMZ8wfPAAwTQQvKerP2ZCZrKh92HIvocXMDk35z/7otxOTfHIWzWcc50KLWD6M/YMHAKaF4KXwZcHQmazYfdiyLqFlnHWJejuxKS1X0Xzmcc6tgDAbs3xI5gWQy8IHr5hlwfKDz547k/XBZ0feN9k+bClCUsYzUNTbyVzcHzzWspROnx5eqXj6NHVb0xAYZmOXD8m8AHJY+OAVsyzY++7Pq6stM1nf/fmxz5tkH7ZUISnTGSjq7WSuhwoe6/LyMHRJw6/Lk2c7ESAwdLN8CKANFj54xSwLFu+5Qf17b1J59nr1djys4j2W3mWWAAAV60lEQVT/YYYj3RzEfP2ZHvx2GlDcHzTWU6ekpaVh6FpaGv6Mi9PradC5XuXGXvU6j6gYk6hYPgTQBuYt6B+0srLiJ06cSPb8oTVewzsHnqljzuhJLu3LK2bbnNZsQxQy1twNedsktPXDQNr35nVVlanbdfUf7PCRAmg0M3vc3VdG/W7hZ7ykyGXBkGmPmJPvHJ6og99S5veepLifaZcwMa0fSqk629G6S9VZyuYAtNtS7gHMpVFn9GncN8ZgIB05Mvw6zfsGCH5Lqd57oM2ysU7nAov7xykKaXWVdDBJxOcZdZwAoOGY8UohphA8RdF45hm3rFsbRYgu7g/dV3HRBe6EEFK3JTGJCGC+ELxSiDlTpDirxKyhJWjTUBRS/+iWurliTN1cwjNqkj2y2Vdxe4FBfqBC+6yvSqauufrq0PoBwEIgeKUSc6aY9lkl0YxbTMF8cWifiqqS1rrSngmzaAnOqEm6zJcl+yqGCAzy1G0BWFQEr3mUYMYtKsxkbnaapMv8ou+rGBpQA5cQF/3jBLC4CF7SXLZzmPaMW1SWynxWjeoyH7NP5aIWGkUE1NAlxEX+OAEsNoJXqpYGcxbmej2pu2Nd1YbU3SH1ep3xd858Vg1++diAuKiFRhEBNWYJcVE/TgCLjeCVYllsDntzFRqo76sqtVc9f0SFjkiTyqET1W5Ntcs80y7BVyCGBlSWEAFgMoJXijNF5hqnWEGBpixVrD+swh+S1jszf09zmGXzi2gKOzj62PZXqYosCwDbIXilOFO06M/+4ECT+T1FZ1m299le4Ic6GEj7Du0ZfkxrUn/PFFt0AMCCIXhJ0z9TpPyzf8q1Y2UpVadd6xum6rSrLG3002aeyojKfaGBqmUzk1MX+KEu+scEANNE8EolxZ/9KbrML59Ud+P1qrRT3Y0z6i0/K2lCw9NMZ9yo3BeaFFo0MxklontsyBLivH5MAJADwatNUnSZP3Wv+kufULlxg3pLaypO3aKxwSuz4NwXmhTmsSAppvVD4BLiPH5MAJALwatNYrrM37lly56DE4JUr6di12EV1aP1c94x9WFPkqTrRmwD2XlKErGtHyLamM3TxwQAuRC82iS0y/ydJ7Xv516vSm9Q9/5KfZ0cH74yTmdEr5xOvZ9Ei0R0jqf1AwA0F8GrbUK6zB8/pUpv0Lp2qJKrPH5KxcGLe84UolZOY1PaPDWwjXnvtH4AgEYjeMVqwQm9d+3X1b2/UiVXV2fUu/bruYc0UtSMS0xKm7c2ERHvndYPANBsBK8YDTihh+S+4lVPq293qPR/op59UsWrbpV02wxHGSZqxiUmpc1b/4OI9z5vbx0A5g3BK0bms1pMs9PiksMqqsfqE/V/nNkYt451qrXtMSmtTcVLQUk6/L236a0DwCIieMVoQvf2FjQ7jZoYTFEw35bipZjWDypUqlBPE3fIbM1bB4BFRfCKURTS0aPS8ePS/v1TO6uFtn5oS7PT4InBlEu3bSheitmyJ+JjasNbB4BFRfCKMRhIhw4Nz4Bra9KebSqXQ54yovVDW5qdBk8MzmtB0pRbP8zrxwQAi4jgFSPBGTCq9UPmZqehgpe7YpduW3BFaYrWD9RtAcD8IHjFiOkcHzrpsX/53NYP+5fH37lFBTxBy10x76cBV5QGSdD6oUWHHQCwDYJXjNDO8TGTHgf3qP/s3Srv+pp6P/ZqFQe3afuQsYAn2fY+IU/WlvW2RK0fqNsCgPlA8JKiEkXI1WWxHdmL3/jnKqpK+o2udFszO7In3d4nRBPW22j9AAC4SASvwUCD3qrKM3vV27mqojxy0TNZ89iRPen2PiFyr7fR+gEAMAULH7wGx76gfdV9qtRVt6rUP/YxFWPOgKHhoyik/tGTQfvltaUje7IwGSPnehutHwAAU7DwwavUjarUffmqQt04dpYiOHwMBioO7RsuH651pT2Tr2xrQ0f2ZNv7tOFKRYnWDwCAqVj44NU7sFvdD6+rqtbV7S6pd2D32PsGh4/Ys2/mjuzZtvdpwpWKMW8+4D1RtwUAmGThg1dRSP0HO+G71tQVPJpUwZPy7Dvltalk2SdknLmnhyLffEjtFnVbAIBJZh68zOwKScck/X1JG5LudPf/NOtxbBWcZUJP1Im2Fkohumh+nq5UjO25Fd4XtcmHHACQUY4Zr7OSftndnzCzV0p63MwecPc/yzCWOKEn6gRbC6USU7c2d1cqJuq5BQDAODMPXu7+FUlfqb//ppk9LekySfmC15T31mvTWTpZ3VrMAFJ8NoE9t0K27JHyT84BAOZD1hovM7tS0nWSHss2iNg1pJCU0rKzdFD2adO+ioHHNHTLHin/5BwAYD5kC15m9gpJxyUdcvdvjPj9QUkHJem1r31tuoGkuAKxAWfpqeeeNu2rGHhMU118CgDAOFmCl5nt1DB0/Y673zXqPu5+p6Q7JWllZcWTDSbV7FTmPRWzXako5V9qDTymLZuYBADMgRxXNZqkD0p62t1/fdav/y0aMDs1bdG5py1XK0b03Aqp3ZrDQw8AaDhzTzeZNPIFza6XtCbppIbtJCTpfe5+37jHrKys+IkTJ2YxvDymHHyiZrxSTY9NO8zF7JXYgL6sAIDFZWaPu/vKqN/luKrxYUk269dtrAQpIWomJ2Z6LCZMTXupNWKcuVc6AQAYZ+E712eXKCUE557QZcHc00gRy5fUbgEAmorglVtkSsh2tWKqaaQp123Vd6V2CwDQSDOv8boQyWu8YtJMiv5Ugc+ZddIpxYtTtwUAmEONqvFqnJgzeqqzf+C6YNZ9FVNMI1G3BQBYMASvmDN65rN/sn0VI5b7pvp+qdsCACwYglfMGT3z2T/Jvoo51/Co2wIALBiCV8wZPdHZf+pdGmICYsai+Zi9EiW27AEAtB/BS4o7o0/57J9kwikmIKaYxQt8U9RtAQAWzVLuASy6UeFjosFAOnJk+HUaNkPa4cPTW2YMfFObma/ToW4LALAYmPHKLGrCKXR6LHYaLUHR/KBzvcqNvep1HlEx5k1RtwUAWDQEr8ySbO+TeQ1voEL7rK9Kpq65+upo3KtTtwUAWCQEr0SSbGsYOj2W+erLspSqsx2tu1SdpXYLAIBNBC9p6s1GBwNp35vXVVWmbtfVf7Az2+19Uq3hBX5O9NwCAGA0gleCywrLY8+pOn2Z1tVRdfqMymPPqyh2Tx7DtEPStNfwBgMNeqsqz+xVb+eqivLI2OendgsAgNEIXgnqoXp6SF29Q5VcXZ1RTw9JOjD6zk3YsijA4NgXtK+6T5W66laV+sc+pmLCa1O7BQDAt6KdRGRPg5BuDsWBq9Tv3qzD9n71uzerOHDV+DvH9JOI7j0xPaVuVKWu1rVDlXaq1I0ze20AAOYFM14R62LBE05FoaI8oqIspd74JTlJrdmyqHdgt7ofXldVravbXVLvwISlUwAAMBLBS8P2B6UK9aSxbQ+kRF0aGrBlUUiNWVFI/Qc71G0BAHARzN1zj2FbKysrfuLEiSTPnaTEKufG07HOKZp/ZGLRPAAA2J6ZPe7uK6N+t/AzXjGzWMETTrFTYymuagwUWzQPAAAu3MIHr9iyqaCr9WKeNPPs2LlF865SN05cbgUAABdu4YNXdNlUyOxUzJNm3t6HonkAAGZn4YOXFNFzKmZ2KvRJE12pGLp6SdE8AACzQ/CKkWJ2KsGVirGrlzQ7BQBgNgheMVL10Zpy8ilLqTrtWt8wVaddZWkEKwAAGoDgFaMlmxD2lk+qu/F6Vdqp7sYZ9ZaflbQn97AAAFh4BK9YLViXK07dq/7SJ1Ru3KDe0pqKU7eI4AUAQH4ELymuj1bofRP15gp62l5Pxa7DKqpH6yXRO6b2+gAA4MIRvFK0rk/Umytmr8g2LIkCALBolnIPILtRVype7H1jnjPRUFUU0uoqoQsAgAYheG1eqdjpbH+lYuh9Y54zdqg71tWxdXV3rE/tokoAADAbLDXGLMuF3jfRUl+hgfq+qlJ71fNHVOiIxAY/AAC0BsFLirtSMfS+Ec8ZXIdflirWH1bhD0nrnZlvLwQAAC4OwSuzqDr8VA1cAQDATBC8MovahYirFQEAaDWCV6wp9+eKnsRqQQNXAAAwGsErRoL+XEUh9Y+eVHn8lHr7l1UUdJgHAGBeEbxiRK0LBhoMVBzap6KqpLWutGc6zVYBAEDz0McrRkR/rsFAOnJk+HWiRM1WAQBA8zDjFSOwuJ0rFQEAwCgEr1gBxe1cqQgAAEYheCXAlYoAAGAUglcCXKkIAABGIXilwJWKAABgBK5qlCIuQQzElYoAAGAEZrwim6IGNa7nSkUAADACwSviEsTgjMaVigAAYASCV8TsVFlK1WnX+oapOu0qS5vcJoLABQAAtiB4RcxO9ZZPqrvxelXaqe7GGfWWn5XEFYsAACAMwUsKnp0qTt2r/tInVG7coN7SmopTt4jgBQAAQhG8YvR6KnYdVlE9Wi9L3pF7RAAAoEUIXjEomgcAABeB4BWLonkAAHCBaKAKAAAwIwQvAACAGckSvMzsJjN7xsy+aGbvzTEGAACAWZt58DKzjqQPSPoRSVdLut3Mrp71OAAAAGYtx4zXD0n6orv/ubtXkj4q6W0ZxgEAADBTOYLXZZL+csvPz9e3AQAAzLUcwctG3Obfciezg2Z2wsxOvPjiizMYFgAAQFo5gtfzkq7Y8vPlkr58/p3c/U53X3H3lUsvvXRmgwMAAEglR/D6tKSrzOx7zKwr6Z2S7skwDgAAgJmaeed6dz9rZv9K0v+Q1JH0IXf/7KzHAQAAMGtZtgxy9/sk3ZfjtQEAAHKhcz0AAMCMELwAAABmhOAFAAAwIwQvAACAGSF4AQAAzAjBCwAAYEYIXgAAADNC8AIAAJgRghcAAMCMELwAAABmxNw99xi2ZWYvSnruAh/+XZL+zxSHg+njGLUDx6kdOE7NxzFqh4s5Trvd/dJRv2hF8LoYZnbC3VdyjwPjcYzagePUDhyn5uMYtUOq48RSIwAAwIwQvAAAAGZkEYLXnbkHgG1xjNqB49QOHKfm4xi1Q5LjNPc1XgAAAE2xCDNeAAAAjTC3wcvMbjKzZ8zsi2b23tzjwZCZfcjMXjCzp7bc9moze8DMvlB//bs5x7jozOwKM3vQzJ42s8+a2bvq2zlODWJml5jZp8zsM/Vx+nf17RynBjKzjpn9iZndW//McWoYM/uSmZ00syfN7ER929SP01wGLzPrSPqApB+RdLWk283s6ryjQu23Jd103m3vldR396sk9eufkc9ZSb/s7m+Q9CZJv1D/++E4NctpSW9x92skXSvpJjN7kzhOTfUuSU9v+Znj1Exvdvdrt7SRmPpxmsvgJemHJH3R3f/c3StJH5X0tsxjgiR3/6Skr51389skfaT+/iOSbpvpoHAOd/+Kuz9Rf/9NDU8Wl4nj1Cg+9FL94876fy6OU+OY2eWSbpH0W1tu5ji1w9SP07wGr8sk/eWWn5+vb0Mz/T13/4o0POlLek3m8aBmZldKuk7SY+I4NU69fPWkpBckPeDuHKdmOirpPZI2ttzGcWoel3S/mT1uZgfr26Z+nHZc7BM0lI24jcs3gQhm9gpJxyUdcvdvmI36Z4Wc3H1d0rVm9ipJHzez7889JpzLzG6V9IK7P25mvdzjwUR73f3LZvYaSQ+Y2edSvMi8zng9L+mKLT9fLunLmcaC7X3VzP6BJNVfX8g8noVnZjs1DF2/4+531TdznBrK3b8uqdSwfpLj1Cx7Jf1TM/uShmUvbzGz/yaOU+O4+5frry9I+riGZUtTP07zGrw+LekqM/seM+tKeqekezKPCePdI+mn6u9/StIfZhzLwrPh1NYHJT3t7r++5VccpwYxs0vrmS6Z2bdJequkz4nj1Cjuvurul7v7lRqei/7Y3X9SHKdGMbPvMLNXbn4v6YclPaUEx2luG6ia2c0arqt3JH3I3X8t85Agycx+T1JPw13fvyrp30q6W9IfSHqtpP8l6cfd/fwCfMyImV0vaU3SSb1ck/I+Deu8OE4NYWY/oGGxb0fDP6L/wN1/1cyWxXFqpHqp8d3ufivHqVnM7HUaznJJwzKs33X3X0txnOY2eAEAADTNvC41AgAANA7BCwAAYEYIXgAAADNC8AIAAJgRghcAAMCMELwAAABmhOAFAAAwIwQvAEmZ2dvNzM3sexM9/0spnjfgdd9vZu++yOf4NjN7qN7s+q1m9l/NrGtmnzSzed1LF1hoBC8Aqd0u6WENt0vJyoai/3/vQh8X8PiflnRXvdn1NZI+4+6VpL6kf3ahrweguQheAJIxs1douEnwz6gOXmZ2pZk9bWb/xcw+a2b313sNbj7m35jZ58zsATP7PTN7d/2Yp7bc591m9v4Rr3e3mT1eP+/B817vNyU9IemKEY/7JTN7qv7foXGPM7NfMbNnzOx/SvqHWx7/k2b2KTN70sz+cz2Dte3rSvoJvbz32zWSnqy/v7v+HYA5Q/ACkNJtkv7I3T8v6Wtm9oP17VdJ+oC7f5+kr0vaL0lmtlJ/f52kH5O0Evl6P+3ub6wf94v1PmvSMCQdc/fr3P25rQ8wszdK+peS/rGkN0n6WTO77vzHabi/6Du3jO0f1Y9/g4azU3vd/VpJ63o5NE163a6k17n7l+qbrpH0mfr7pzafH8B8IXgBSOl2SR+tv/9o/bMk/YW7b87uPC7pyvr76yX9obv/rbt/U9J/j3y9XzSzz0h6VMMZpqvq259z90fHPOZ6SR93979295ck3SXphhGPu6G+39+4+zck3VPfvk/SGyV92syerH9+XcDrfpeGoVNmtlPS33H3FyWpXnqszOyVMW8eQPNRvAkgiXq26S2Svt/MXFJHkkv6TUmnt9x1XdLmUqONebqzOvcPxUtGvF5P0lslFe7+N2ZWbrnfX08a6oTfnf84H/P4j7j76nnjuXKb1/3bLeO7WtLT5/1+l6T/N+HxAFqIGS8AqbxDw2W23e5+pbtfIekvJF0+4TEPS/pRM7ukrg+7pb79q5JeY2bLZrZL0q0jHvudkv6qDl3fq+GyYYhPSrrNzL7dzL5D0tslrY2539vrKxFfKelH69v7kt5hZq+RJDN7tZnt3u5F3f2vJHXM7BKdW9+1GVpfdPczge8BQEsw4wUgldsl/fvzbjsu6X3jHuDunzazezSsdXpO0glJ/9fdz5jZr0p6TMPw9rkRD/8jST9vZn8q6RkNlxu35e5PmNlvS/pUfdNvufuf1DNW59/v9zUMSM+pDmfu/mdm9q8l3V9fuXhG0i9I+t8BL3+/hkud12x5fUl6s6T7QsYPoF3MfdTMOQDkYWavcPeXzOzbNZxlOujuT+QeVwp1Ef8vufu/OO/2uyStuvszeUYGIBVmvAA0zZ1mdrWG9U8fmdfQJUn1zNqDZtapC+o3r3a8m9AFzCdmvAAAAGaE4noAAIAZIXgBAADMCMELAABgRgheAAAAM0LwAgAAmBGCFwAAwIwQvAAAAGaE4AUAADAj/x9LFbPOwaLjGwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sls = []\n",
    "sfs = []\n",
    "for mode_name in s_modes:\n",
    "    l = int(s_modes[mode_name]['l'])\n",
    "    if l < 50:\n",
    "        sls.append(l)\n",
    "        sfs.append((1.0/s_modes[mode_name]['t'])*1E3)\n",
    "    \n",
    "tls = []\n",
    "tfs = []\n",
    "for mode_name in t_modes:\n",
    "    l = int(t_modes[mode_name]['l'])\n",
    "    if l < 50:\n",
    "        tls.append(l)\n",
    "        tfs.append((1.0/t_modes[mode_name]['t'])*1E3)    \n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10,10))\n",
    "\n",
    "ax.plot(sls, sfs, 'r.', label='spheroidal modes')\n",
    "ax.plot(tls, tfs, 'b.', label='toroidal modes')\n",
    "\n",
    "ax.set_xlabel('Angular order ($l$)')\n",
    "ax.set_ylabel('Frequency (mHz)')\n",
    "ax.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}