{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ICPertFLRW\n",
    "\n",
    "Example of the initial conditions generated by [ICPertFLRW](https://github.com/robynlm/ICPertFLRW) and [ICPertFLRW_GRH](https://github.com/robynlm/ICPertFLRW_GRH)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import aurel\n",
    "from aurel.solutions import LCDM as sol\n",
    "#from aurel.solutions import EdS as sol\n",
    "from aurel.solutions import ICPertFLRW as IC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4th order finite difference schemes are defined\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "Cosmological constant set to AurelCore.Lambda = 1.04e-07"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define grid and classes\n",
    "L = 1821\n",
    "N = 64\n",
    "grid = {\n",
    "    't': 1.0,\n",
    "    'Lx' : L, 'Ly' : L, 'Lz' : L,\n",
    "    'xmin': -L / 2, 'ymin': -L / 2, 'zmin': -L / 2,\n",
    "    'Nx' : N, 'Ny' : N, 'Nz' : N,\n",
    "    'dx' : L / N, 'dy' : L / N, 'dz' : L / N,\n",
    "}\n",
    "fd = aurel.FiniteDifference(grid)\n",
    "rel = aurel.AurelCore(fd, Lambda=sol.Lambda)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define perturbation amplitude and wavelength\n",
    "amp = [0.01138486133517004756]*3\n",
    "lamb = [L]*3\n",
    "\n",
    "# Define initial conditions\n",
    "Rc = IC.Rc_func(fd.x, fd.y, fd.z, amp, lamb)\n",
    "rel.data[\"gammadown3\"] = IC.gammadown3(sol, fd, grid['t'], Rc)\n",
    "rel.data[\"Kdown3\"] = IC.Kdown3(sol, fd, grid['t'], Rc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "Calculated gammaup3: $\\gamma^{ij}$ Spatial metric with spatial indices up"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "Calculated s_Gamma_udd3: ${}^{(3)}{\\Gamma^{k}}_{ij}$ Christoffel symbols of spatial metric with mixed spatial indices"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "Calculated s_Ricci_down3: ${}^{(3)}R_{ij}$ Ricci tensor of spatial metric with spatial indices down"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "Calculated s_RicciS: ${}^{(3)}R$ Ricci scalar of spatial metric"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "Calculated Ktrace: $K = \\gamma^{ij}K_{ij}$ Trace of extrinsic curvature"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "Calculated Kup3: $K^{ij}$ Extrinsic curvature with spatial indices up"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "Calculated rho_n_fromHam: $\\rho^{\\{n\\}}$ Energy density in the $n^\\mu$ frame computed from the Hamiltonian constraint"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'y/L')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAHaCAYAAADhZFb0AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAWYNJREFUeJztvV2sXNl137nPqY97SUndFBuWYbUalkl/TPxgWGw2kIEygG2xbc/EmcASyc6HY8+M1eTAmZdxJNJ0MiPZHpgiYWASv4zItvwwmIfpJmUbk8CDgGzHGdiJxyYpGBEmNiBSTtqtxLD50VRE3lsf5wzWuV30ZdX61z37nlNV51b9fkKBrXV37dq1z65Tq9Ze/72SPM/zAAAAALDCpIseAAAAAMCiwSECAACAlQeHCAAAAFYeHCIAAABYeXCIAAAAYOXBIQIAAICVB4cIAAAAVh4cIgAAAFh5cIgAAABg5cEhAoA9wYMHD8LLL7+86GEAwJKCQwSwR7h48WLhECRJEt7//vcX/339+vWwKrzxxhvh0KFDix4GACwpOESwEM6ePVt8odsXuz0OHz5c6xe89VUlmlD1+aOIxosvvli8N3NiqnD16tVw5MiRcO3atWDlB8+dOxdu3LhR/P9V4cqVK+H06dMTDqKtn+1O4vY2e4k7d+4UY7f3ZQ/771u3bu2qL/sc2fNPnDhRrEGbF9WXfRatnT3H1qxhbW0MZgdYGay4K8CiOHXqlBUXzq9du1Zrv4cOHcoPHDgg/37mzJlKz4/BXqvqR83m6f79+/m82Wme5oW9d3U9ZrWG5snt27eLNTd+jY8cORL9vi5cuJBfuXLlKdulS5eKOfKu52h9jj+OHTu2kDUHsCjai3bIYLU5cOBA8e/Bgwdr7ff27ds7/hqv8vwY7Ne5/dquOk+XL18OZ86cCfNkp3ma53bZqVOn5rqG5olFYixSM3ovIy5cuFD87f79+6X6sciO9XH8+PGn7DZ3tqZtHb700ksTf7fo2x/8wR8UESKLaB47dqyISAKsEmyZwUqy262IRWFOlX1hjrY0Vm2eLl26FF555ZWwjNgc28Nz+MwxGW2Zlp0n5TiaczVqM445SPZ3+5s53ThDsIrgEMHKYbkSTYl8lMW+GO2X/8c+9rGVmydzAu2xrF/Sr7/++kRkaDuWSG5tyl4zi/BM62uVEvEBYsAhgpXCvuD3aqKo/YK3SIJFilZpnmyrcHyLZ5kwB2Waei7GiTHHyq7dvCOJAMsAOUTQSOyGbtEQ+9du8KasMmdg9MVg+RB28x9tA4xvL9lzRs/b/sVquRKWazJ+po399yg/Rz1/O9bX6EvHxmK/ymeZ32NjuXnzZvGeLQ/Eto92EzEpM+6y81QGU34ZJ0+eLF5rFAmxfs2xs+0dbwtnOxYdee2118KyYnMxzSEaXYcy2BrZaR1NWzf2GRu14YgDWDkWls4NsE3hcvPmzakKIvv7uNrGVDlKBTVN2TVS3JQZl1LxjKtvrL2Nx8PGXeWjZnMwUpnZa1tfpj6KJXbcZeZpJ0wZ5imV7P2ov3nqqyprqOnYPJiiSzH6DFRVfJnyzPqx6zo+f/Y3s4/m0P4/KjNYNdgyg0Yzik5YZGiUYDrCtlFUsmnVM4TU8+2X+vnz5wvV03YsUmW/rC26Uie2bXXv3r0iimLRlVF0xn7Jx2yLzHvcIywyNJ4fY9dyFIWaljtj2PVd5u0yY17bW3b9LfIznnQ9ityZfRQ9sjkfnWEEsCrgEMGewAvz2418UUm/3peYbTHUKde3LSVzHsa3i0YOgh3M2MRxb3+t8S9Us9kXrX35jju4HuYILuKgRXvN0aGhu3k07XBI22a1z8qbb7458Td1Lcxuz6l6ZATAXoEcItgTNCWfwX5Jlz0TpgrmONgXkUVwxqMoo7mIiSzMa9zjr2kRou28+uqrxb875Q0ZI2d3EdfexldmjPNkp2jatHm06NAoBy32h4jlcM37/CuARYBDBHuCJh66Z9s5VjpjtOVgW1t1YV9ghrdd9NxzzxX/7laGPstxj7P9C9hed/TaZag7OjRK0G+afH/W18AichYZ2o1jaZ87ZPqwKrBlBhCJ5b/Ytsgot8d+Pdv2Qp1Om+UI2ReY9yV29+7d4t/YL7h5jHuaM2LRobJbZbPIH9rNFuM8OHr06NRonzlxu42SWS6cRRmVE0ieEMBfgkME8O6Xb5l8JGtnUQtLCFYnAtf15a2+BM1ZKutU1D3usvPkffGa4zW+DaXOVLLXsMhJndtloz7LYOOaVw6ROS3TIkQ27tjrbdgY7H2MP3e787VTcn4TI2oAswKHCOBdykRKbCvLvlS9L6jtX2pVT3k2R8D78rYvL+s79nDGOscdG1EyJ8r6NWdskaU6rM+yDpZFVSznarePmPwjuyaj07g97HrEqiYt/8ye413v7YpCi8CphPrRFuNunDGAvQgOEawcXlJy2eiBtbMtDs++/UutqpTavoQ8x8S+zOwXe+yX1G7GXWWexlVlXn2s0ReucqLqjMCN6oU1EZsXe4wfiTCaB69Y67Q1Zs+xa6e2G0dbroY5neroCltr6vBTgGUEhwgWyuimvpuk0mlOxzTHxJwJ+8LwvoB2er59SVu0Y9xuXyqjM33G8z524ySN+tqe0Gpf6BZ58KTTO7GbcZeZp50wZ8j68L5ULWrlRZtGFdvLOl47rSF733bqeVOUih52tILN0fj1sUigF1kbHQsw7uSNSrtY5Xr7d/xh12O7E2qOmP3/cafI/r9dn92sNYC9CiozWAh2c7ab9yjR1W7wdnO2f0dbCHbzHv3dEj/Nbg6B3cCt3fjfRl8o259nX4T2K3hcNmxKJ+vDtgtMtTWKRuz0fHsNa29tbEti9KVtf7fnWr/2N+t71NfoC8j6sihN2e2Ur371q8U8jVRZ1t9upNOx4y4zTzFbZRap2H6WjfVl82vX35Nzm8y7TA6O9WnjGzmN9pzRVqO9J3OQttf1arJDZGvf3otd71FxVpsnc4bUGVye0zhab9PODhp3Tu0ajHLMDJs3c1Rt/e1W6g+wF0nsuOpFDwIAYIRFPvgyBoB5w5YZADQGixpZFA1nCADmDQ4RADQG2060bR8AgHnDlhkANAa2ywBgURAhAoBGYMnRbJcBwKLAIQKARmCKqqZViQeA1QHZ/RzJsix87WtfC+973/tCkiSLHg5Ao/joRz8afvAHfzA8fPhw0UMBkFiWyde//vXwwQ9+MKQpMYVlAodojpgz9MILLyx6GAAAUJG33norfOhDH1qYU/Z93//d4X3v2xf+2T9t5gnsexEcojlikSHjiy++HN7T7jz1t7Vub6L9erfv9rPW3XTt3fXJPoyO03dnzW/bWvdfM+2K9mKMqWNPO0O3bWj79qSV+fa2Y09q0gbkfuQuH0z+EsyH4tfhoOWas76w955eC8bQsW217br24Ybfvr/pt+87/fQ2/LabvTXXviHGuOn0/agv+hj4fTwa+LemDWF/NJyc2w3HVtgz/7ptDv1rv+Esz76/NOPtzrJVK1nJX5RdBaE9u4pXd8QfOml5u2q77l+esNby39B66k/iemvyAu13bEXb9sC17xf29fbkfezRsB9+7A//2ZP7+SL49V//9fDHf/QfwuZmv8i9o95cPeAQzZHRNpk5Q+MO0Xp78ibg2Yw14Vgoe6fbKmUz2l2/j3TNv6u1umnp9mlHfCu0/btu0hJ2r31dketMOUST9lx8gQanbdF1yx9klkzah4loK97oIPOvZz8vb98c+reDTubb21mntD0RbdPctwdhT3J/LHlw3qd473KxKGfYsSsHItqelXeIxKdHtlcOURrjEKXV7eIWEdZS5fjEOUT7HOfHs23Z/cHsF/eafa2d7+fzpt/vh5/5mb8fPvNzHw9f//pG+PSnPxlu3rzD9l0NMIMAAAB7hP/t0idDu52G/+4nvy/8/f/h5XD//jfC//5//NSih7UU4BABAADsAd55553wCz/36+EXL/yt0G63wvp6N/wvv3gy/M//6Gp4/Pjxooe358EhAgAA2AN87uKPh7/y3c+HH/mRjzyxvfK3/mr4pm96X/jHv/zfLnRsywAOEQAAQMP50z/90/DL//ifh4u/9Lefyl+y3KELv/R3wud+8f8Kf/EXf7HQMe51SKpeAKa4UUmms8BLCvVs0+xtkWwchGIndxQ+ecdXciTKLpIrEy9hUrStLanaez+qbd//WGUR9uGmvz4GSjUWae859s1NpSbrRtkfOaq0x1JN5tu/IdRkj6XKbHIdPnaUgcaGuG4bESqz3rC8amy6+mzyCblIk3aabtkjf+36SdX+e+94jaclVTvNhbgyDMW9Zije51AkDWcR9zE1VzF462Ee/KP/6e+Ev/FffyS89NLhib99//d/d/hr/8V3hZ//hb8bfvmf/POFjG8ZIEIEAADQYP7wD/8wvP5//l6RL6T43MW/HV67/C/CV77ylbmObZnAIQIAAGgwZ87+vfDf/9Sx8G3f9gHZ5ru/+/nwY3/vr4VzP/t35zq2ZQKHCAAAoKFcu3Yt/P7/ezv87D/8mzu2tbOJ/u/f/MPwe7/3e3MZ27KBQwQAANBAhsNhcfDiuX/4N8PBg+/dsf23fMv7w//4D/6r8A8+9RNFeQ+Ig6TqBWDlCuQJvfNKqlYJwZH2TJSvaHUml1YqkqeV3U2eVsnW4hTbaFSSuPP+vUTroovIpOqhk3WqSndEJ0/3ytvrSJ5WCdSzTJ7eap+WbzuMS5bd9JKqxXLrZVnp5Glj6HxpDcQXmUq2jsVLoG6LU5d7WVyydddJfB6oJGlpVye9J434ynukBj4D7MDFBw8eFQcwluVTn/7r4bVLvxV+4zd+I/zoj/7oTMe3bBAhAgAAaBh20KIduGiJ1HYAY1ne+971Yuvs7NmfKsp8QHlwiAAAABrG//pP/pvwgQ88Uxy8GIuV9Wi10nDptVdnMrZlBYcIAACgQfz5n/95uHD+nxYHLu6maKuV9Th/4ZXw85/99fDw4cOZjHEZwSECAABoED//Cz9WHLRoBy7ulh/5G0fCf/ZXPliU+4By4BABAAA0iH/52/82/OQnv69SH1be45Ovfl/47X/x/9U2rmUHldkCeGQKmhIqsySpR82gjrGPaZvnvu/cFkqrfDiZzJcK1U8q1FdSZdYuf3a+msOYOSnaD1rlVWbifSqV2cBRlA36SmUm7KK9Ksex6bRXarLHfWGPKMfxuCY12Tec67DVz+T13BBKsEeDOJXZhqMqkmoyoRDr56r95IsORIGJmarMxOe7k/jz3RHtPTVdX2z5SDVZ5Nv07k2ZKEWi8dehN/KNYU2K1qkcCK1WErKsWp0QyyNKwjfVNqplhwgRAAAArDxEiAAAABpGng9DnotwZmkWVIl2j4JDBAAA0DiywimqQi62asEHhwgAAKBhZPkwZBUjRFUdqlWDHCIAAABYeYgQLYCNQTsk+dNTnzpqqLpUZjHEqsykosqpCdYSiqJWu3ots1kr8mJqmQ3F+xwKhVSUykwo1VTNMk9NphRlsWqyjWF55dgjMVcxtcm27P71eeQslcexajKhHPNURZtiK2JT/KLvi1yOQTJpH4q2WVLP9kfqfJZbwb8+7VypzHz7WtIuVa/NyMXvcXWvUbVKM3cZ+uskjYwBJI7K7rFYy7WT25ZZxQiRUCyCDw4RAABAw8hD9aRqtsziwCECAABoGHk2DHlWUWWGQxQFOUQAAACw8qxMhOjixYvh7t274bnnngu3b98OL7/8cjh+/HilPl988cVw8+bN2sYIAACwxdCOyK/WBbL7KFbCITp9+nQ4fPhwuHDhwhObOUT37t0Lp06d2nWft27dqnGUAAAA9R3MSA5RHEvvEJnTcvny5ZCPSRTMObIIz24cIuvzxo0bux7TxrAVkjFVh6cyi93njOmjLvVVJlVpk/ZsOIxSX6Wp/+smddRnnvJsN3hqMiNzlCWekk61nfY+vTpkA6EmkzXLItRkRXtHOebZpqnMVB0yW9+TbVUNsjRKfabqk3mKMk95tjW+LEpl9tjJ4+gF/4tqM+m59l4yWdvPGDr9DJKBTLKtA0851R5TvY5oOaoxoytqMQ7zyfU2EH1n4utH1TLLIjI8UlHKLBHrSvbj2DayOanMrI5Z5q+b0uAQRbH0OUSXLl0KR44cmbCPbFevXo3u8/XXXw+vvPJKLeMDAACAxbP0DtH169fDoUOH3L8dOHAgXLt2LToX6dy5czWNDgAAwGNry6zagwhRDEu/ZXbnzp1w7Ngx928HDx6M2vqyrTJzrsyRKsPm5mbxGPHw4cPSrwUAACtMsWVGUvU8WeoI0YMHD6b+3RybndqMb5XFKNPOnz8fnn322SePF154ofRzAQBghcnfdYiqPIgQRbH0EaK62M1WmbX/6Z/+6aciROYUbRTJuE/7oqmT6CjyAkMSmpM8nalj9p3k5PbQ/7WTtvwkxVZblO5wEnRnnVTtlRNQSdKqnMlAJCH3naRlWaIjMqlaJko7ZToeqSTpiOTpom/n/W+teacPMd+PxX38UUQ5DpU8/Ugk928IRc9jJ1FaJU9vJhuuvR/+MlK8nWHol0q0NrKakqpTp0yHSp5uBVFCJllz7QPnC3jgJFobuciezkQStsIrx5Gom6fqQzyhlZVfy7WT1yG7xyGKYakdop22tspGh2K3ykasra0VDwAAAGg2S71lthN2DlEZJyd2qwwAAKASeRaSbFDtQQ5RFEsdITLM4THHR0WIjh49OvX5Jsu3CJEdxLidUTL2yG7nGsVGkAAAAGaWVG19QGmW3iE6efLkVCWZnVg9DYsMedGh0UnVds4RAADATJKqq/YBpVn6LbMTJ04Ujst4vpCdT2QoST4AAACsDksfITKHxyI8JoHfXsvM/vvKlSsT21xW88ywArDTUNtwZdgwFc6YaqnliByUUKIlSlq0hWojcZRGiSjzoex1qLVUqQulJlNqLW+MSmWm3o9S02mV2aR9KI7wV+ozpTIbOGotpSbriT56Q6UQE/04Y/ds09Rk2u4oc4b+vGp7iLM7ZTdUKY4YNdmWfaO0mqwXHsWpzPJJlVnmKM8Ku8gHyYNvT8Tv3TSZtKdCTdZKRIkONcZkciyZGJ9E3IJS8dn0ynSoEh2qpEdLqcwc86YYR90k+TAkFVVm5BDFsfQOkWGOj8nmz549+6TavW15eVth5iDZgY0Kq4tmp1uPSn5YPTTLQ2LrDAAAasOcmao5QDhEUayEQ2ScOXOmVLubN29O/bsVg91NQVgAAICyJNmwUIpVzyFa+syY2mCmAAAAYOVZmQgRAADA6m2ZEfcoCw4RAADAMsruhbAAfHCIFkAvSycUIC3nh0BLKKRaQt2TCnmGV/ssFeqRtC6VmaPKksouVSctra4yiyVGZaZUczFqMqPvtFdqsr5Qk6maZT2lBHOUcEo1ZuvV7UMoeTa9+k+RarJNqTLzr7NXt+xxVl1NVvSdfGNyfOGx27aX+yqzQe6rzDJHZTYUyqK8plpmSUQtM3WfGCblVWZ5UtO4lcrM+RpLhEY3Rk2m2m/OycdIMjuputrcoTKLg1gaAAAArDxEiAAAAJoGsvu5g0MEAADQNOrYMiOHKAocIgAAgKZBhGjukEMEAAAAKw8RogXQyybrlHlqho5QVSiVWVup0pwaZy0RSm0JdZMiRpUWqzJLRc22lifJEz+kYmuZKYbOnEuVmbg+MXalJpPqMzEWVZ/MrWUmVWb+XKmaTra+y9i2XlPYVXuxbjedX8K94Ku1NoXKTNUn8xRlSk3Wz3312SBTtcwm7bmoUJ7XVMsscWqZDRP/2reSNf81xWdTDESMQ6ll/XG3c3HfG6sLabTE2uyIz4Mopxg6zlD681KZ5dW3zIgQxYFDBAAA0DTM+a/qEJFDFAUOEQAAQMMoqt1zDtFcwSECAABYxggRDlEUJFUDAADAykOECAAAoGnkeeVzhDiHKA4cogVgyppkTOXUScrVhDLaom6XUhq1HXtLqJhUyFCptZS9DlJV48x5P7Hji66rFlPLLFJl1nOuRT+i7TS7WhOevS/eu1qHA6Uyc+qWKWVOP49Tk/Vz/wmbTv2vWDVZL/jKMU9RptRk/czvYyhUZpkzbs82TX0WhMpMfZo9dVcqapmpOmSyrlrEnkPqqN2m2dtKCeeozDqObdq66oi6id66nZfKjC2z+YNDBAAA0DQ4mHHukEMEAAAAKw8RIgAAgIZRHMxYNcIjtpjBB4cIAACgabBlNndwiBbAILeE4Z2TS7siyVUls/ZFMmLHS6AVycadVL2mSGhUya/Bb19HUnXLKRtQV3K3SqoeOnMYm1StEqW9vtV8D8X4VPL0ICKpWq0r+Zpiyj27TKqWdiEcEDf4vpPk20v6ou1mlH3glNeQpTik3U/wzvLJMeYqqVomT0cmVTsJxzJ5Oon7Qk2cz32iEpaTjmtvBd/eS7quvZN3yq8Tde+MWJ9zS6q26E5VlZj4HIEPOUQAAACw8hAhAgAAaBqF7J4conmCQwQAANAw6qh2Ty2zOHCIAAAAmkZWRw4RDlEM5BABAADAykOEaAEMszChiRk66gelZhikcWqggWNvS7WSUJMJFZdSVKVO+6FQPMQqxDwlmOrDG4cqxaH6VnZPHTbNnolyAt4cqj7U9Ym59qq9aquUOUNp9/qIU5MNRe5DX5SvGDgqqWHw1VrD4KvPho7iSynBho7ybKvtoHQfW3ZHfRatMoslrdy3p1RT89LK21HzPRTqQHU9vWuv1slQjFuqGr37cj5P2T05RPMEhwgAAGAZt8xwiKLAIQIAAFjGgxnJIYqCHCIAAABYeYgQAQAANI1Cdl8twoPsPg4cIgAAgKZBDtHcwSFqSC2zYV5+LStFkbQ7iqXcqQc2rfSNrqKUlFZxKWWXUlQpvFpmCjUnddQy02qyJM7uzKGc77z6uAu7p5qTCjv/Nb01O1rfE334TUMu/jIQLzoQMzN0apkNEqUyE0owqT6bbJ8LFZNSman6ZJ6iLBfjq+vLzatPpoSeuciqyHKhLnUUZd78TZ1vqSYT19O5FmqdqHWl1mHufDbn5mLUUcsMhygKHCIAAIBlLN1BUnUUJFUDAADAykOECAAAoGnkU/bIY/qA0uAQAQAANI1aTqpmyywGHCIAAIClLO5KiCgGHKIFYEt0fJl6yz6LFA5IBZLbViiKQvUaX8oe27fCG3tdyXBqXqLeew1zGNu3vp7l14pcVyHO7q+38uPY6iNODZQ5yqncUZ4VbZVd/Jr2+slVW6E+U7XCXHukEioWT1HmKc+2iHyfecR1EHOYObXJpvbjXfvIOYxZnwi3lhccIgAAgKZRyO6r5hDhvcWAQwQAANA0zJmpmgOEQxQFDhEAAEDTIEI0dziHCAAAAFYeIkQNZpYCgXyGycaqvUyeFmbVPlF1BmZITOK3ahszh7HznTd8HcYmbC8Clfgc1zbW3nQik61rmcP5EyNgmVvQpVCZcQ7RPMEhAgAAWMYtM2T3UeAQAQAANI2MnOp5Qw4RAAAArDxEiAAAAJayltnunn/x4sVw9+7d8Nxzz4Xbt2+Hl19+ORw/fnwmfTx48CCcP3+++PfOnTvh3r174dy5c9GvVwc4RAAAAI10iCr2sYvnnz59Ohw+fDhcuHDhic2cGXNUTp06VWsf5gSdPXu2aHfgwIHCduvWrfDiiy8WDtGVK1fCPMEhajBpnNAoitiu00hll9c+Vh22CDVZHWNRbWPmMHa+k4avwyRp/p59EjEa3TbW3nTi3k89czh/0oh1O8vP2oQzU9UhiryF3rp1K1y+fHmi9Ik5LOaklHGIYvqwyNB2Z8g4cuRIYTNH6fr16+HYsWNhXjRnRQIAAMDCuHTpUuGQjDOyXb16tdY+7L/NSRpn5ATNO0KEQwQAANDUKuBVHxFcv349HDp0yP2bRXGuXbtWax/WzrbRvHaG97dZgkMEAADQyFJmSbVH5AGvd+7cCQcPHnT/ZvYbN27U2oc5R/fv33e33YyXXnopzBNyiAAAAJYxqdrOdsyy8PDhw6fMa2trxWM7Dx48mNqVRW12alNHH6NtN2tbNom7LogQAQAANI0sqeXx9ttvh2efffaphyUzN5Xr168Xj9dee+2pZOt5QIRoAVgQMynhmcYoH4r2EaXClIopEZvOSjkVY4/te5YKNtm3sA+d0LN873kNcxV5ffT1LL9WYpVgcn2WfL1p9kSMXNnTfHI0SdLy2wZhT5RyarJ9otqK10yc8Y1GM06e+GGB2kSXzoVWii9pV+/TmRdv/qbNt7o+sh/v2seun4j1OTeVWU08//zz4ctf/vJTtvHokHFgB+ejTGSnjj5OnDhRRIg4hwgAAACK/B/LA6rWSQhpmoZnnnmm8nju3bsnk6Xr6sOcITuUcd5bZSNwiAAAAJrGaNurCpFJ1QcOHJDKLovuHD16dGZ92LlDlkR95syZsChWxiGqehR5k44XBwCAFSDSoan6/JMnT05Vktn35iz6sIMc7bt53Bky+zyjRSvhEFU9irxpx4sDAADUzYkTJwonxL7ztucDWZKzUebU6Ng+zG5tvchQmZyjOll6ldnoGPHxyTbnxhylMkw7XtxO2hxdaAAAgFqoegZRcQ5R3EseO3as+JE/rkKz7zr74T+eNG2BBnvstg/bbbHvYdu1saCDPez/28OCFlVzlmJZ+ghRmWPEd9r2sjb2sIumjhePqbfSTrYe22lFqH5aQm4i7WkWUW/Lf02tkSmvekojxjfNHqPKUq+ZRYaS1Vg81GFoaizeHMr5Vqqs2LnN0tLrR69D3z6+tos+/KZS9dMWL9oWaq2Wo0Bq5/7trZX49jR0SrcfKgWb6DtPhr7dOWhGqcmU+iwWVzkmxp2ouZLtW9XnW3wtyevpXXvxCVLrSqvSytlmgt1HnM9pXB/xz79y5UqRYmLOySjFxBwU73vSnBvvEMayfZjTY06RBS08tu/qzIOld4gseuM5RNuPEd/JITIv1dsTXdTx4gAAsOTUkVS9Sx/6TMnE5ps3b1bqYzzIsGiW3iEy71NFb8oeRa7qt+x0vPjm5mbxGDF+WigAAAA0g6XOIarrGPHdHi9ue6jbTwd94YUXdv1aAACwOuR2DGte8bHnjpFcLEvtEC36eHGT5b/zzjtPHm+99dbcxwkAAHt1yyyt+MAhimGpt8zqOEa8yvHiXgE9AACA+ZxUjUMUw1I7RLM6irzq8eKtdOvxtG1SWtJxbEY7UmXmtVdtO+kwSq3UbvntU0cRE6smSyPaR9dDi1SIDSPUHqoPpWzz5nAo1CGd3J/vvhhfzFpRbdU67IkpaTlSuI5o2xGyuZ74IugIdVc7L69uakk1mW9PHXsrWYtTk0UoxHK5OrNIe/n6ZFpNpuZEKfjWSveh5rsVrT5rlV4nLaEyU+vQW7d5PWI/aCBL7xDVcRR5044XBwCAJWcBpTtWnaXPIbJjxE1pVuUo8jLHiwMAANRJ5aTqyIMZV52ld4hse8vk8eP5QjFHkTfteHEAAFhyKidUp7s6mHGVWfots+3HiG8/9XLaUeTjB0aNjhe3vmzLbLsTNPobAAAA7F2W3iGq4yjyuo8Xt9IGnbGt3fH/PzVJWiS5dlTis2NXbfVrxiVbd5xEYZWALROzVbJ1ROkOZVeJz8ru9ZNFJnIrvLNC1Hyr5G51PYfi/Xjt+6JtKxPrTSWoOmadVK3som+ZbD6ZRNvN/eTcvkiIHoa+b08m7bmY7zz4102ROO8nc95L0bdIqI9OqnYSjlWStEyeTv05bDv2tpjvToizq+vZcUp3dJK0cvK0ss9rF8q2u6qqzEgAj2MlHKKqR5E37XhxAABYbkZ5QNUgqTqGlXGIAAAA9gwWPaxa3JWDGaPAIQIAAGgYtt1VecsM2X0UpKADAADAykOECAAAoHFUzyEiQhQHDtEC6LZC6I4pxTzl2JpScMkSC1lpe7c1EGMblFaNTWvfnqHKrCXa16EyUwyHrZmpzOpQwSkxiVKZDTzlochXUOuwLxQ7tr4n2iq1m7gM3dQfS18o3tYcNdQw77ptB0KtlYnyGsoeE3NPHCWUMcw3J7vI/dtyLiRDubj6XomOwu4osDzl2bQSJZ6azOgk+yZs3WS/27YbfPtavi7s3dLXXikg1bpS63CRKrMnxV2rgEMUBQ4RAABAw6C46/whhwgAAABWHiJEAAAAS3gOEbXM4sAhAgAAaBp5HTlEbALFgEMEAADQNLLqOUSU7ogDh2gBdNOtx3bWHJXZuBLtiV0pvtLydqVI67aVyqy8mqxo7/QjVWbiNZXKzFOOJWKuYlE3IG8smfj1FltXzSOLVcGJX4JDVRPNUYgNWkmUUs1bs6p9X/xQ7fvipjAQU9UXKqGhszcwEGqtgVArZUKtlSfOHEZ+TyVi3C1njMN8UEudNDkWR/HWkrXM/Pphqj6ZpyhbC5PKs92oybri62rNUc0pNdn4ffeJXaxDrz27UMsLDhEAAEDDsKLP1DKbLzhEAAAAy1i6g1pmUeAQAQAANIw8pCGvmBRtUSYoDynoAAAAsPIQIQIAAGgaluNfdcuLLbMocIgWQDfNJupDrbecemNCZbUu1FrrbVVvbNK+1u5HqcmU+qzb6ZevZSb6aIlxt4RCKkZlFlvLTO25e+2HWasWlVmM+ixWlabsXig9E+H1oVDTDaUqzWur1GHl+5jW3rYXxsnE7S0X9dBiULW/UkfxZPSFWmuYT35+suB/prKaapl5Y0yDP76WGHcnrJWuT6bUZPuk3VeZ7Uv967nurK11sTbXhJpsXdqde828fIwaDmaEOHCIAAAAlrCWGUnVcZBDBAAAACsPESIAAIBG1jJDZTZPcIgAAACaRg1bZiRVx4FDtAAsgXp9Iql6WDp5OqZER9GPk0DdrSl5utMRpT6c15QlOpyE8mnJ1kmSzb10h/dLLR2IcYvrIEuRDPzE1RhkkniEfeCU85iWEK1m3EuIlgnbkcnTKh/auz6qbSZKesTUZEgjEpaNlkhaHiaTn5Nh8D8nmVdCZBekXukO8VWgxq2Sqr0EalWKQyVPr4syIuuiHIdn10nSyp6Xt9cggphbtfvaRrMa4BABAAA0DCvPV9khQqUWBUnVAAAAsPIQIQIAAFjGHCIiRFHgEAEAADQMy4urrDLDIYoChwgAAKCJOURVD2Y0hwifqDQ4RAvAFGb7WiVUZkKVtV/Y93V6rt0r07HeFW2lmkyU6FAqM6e9Uo2lrThVltd+1iqzbNgqPW6v7VZ7oTJzxp46Srot++zKf8SeZaJe0SsXokqIDIVyKBOlF6Z9eUz0EZkimYprnzj2du5f47Yo6dFLfEWVpygbJP5nKg/1qMwSR2XWFso7pT7r5r76zFOUdUUfshSHUpMJtaOnHNsvvtn2KZWZuH+Mq4G38D+bsPfBIQIAAGgYdcjui18tRIhKg0MEAADQOOo4hwhvKAYcIgAAgIaRZzUUZ7Xnc7hOaZgqAAAAWHmIEAEAADQMSnfMHxyiBWCKsnFV2T5HObbPUYcp1dg0u6coU2qyrlCfqZplbaU+65ZXmbWEai4RKi5PlVWXykrW/nJC17lQkw0H/scqbfntvbEvQk2mlGCxmhpXZabaRvatzmXJ2nXUIVPtJztvi3G0hL0jVFkDpz7ZMBfqRaE8jCV1xthylGfT1HQd0X7NqUO2lkSqxqTKzDWHfc61Xxd7H/uF0nNfW9mdazGsR+23MzVUu6/4/FUDhwgAAKBh2A8L9SOlynEUoMF9BAAAgJWHCBEAAEATc4jqOKkaSoNDBAAA0DRqOZgRhygGHCIAAIAm1jJDZTZXcIgWgNUo2zemsFhvlVeZqZplMfXJ1tY2o2qWddb6lVVmqVCqyRpfSmUm1GoxqqzYG00+aJVWmaVD/322+kJ95tRLknMilEazVJ8p1BzGJILGXodMnryblFaNJco+9FMqE6fvltjK6AhVT99O2XPtk2t5ILR3eU0Zst77aYt00o6ozdYRyrGOM7ldqRpLotRkqj6Zpyjb186j1GRKfbbfq5tYU005aB44RAAAAMt4DhFbZlHgEAEAADQMq0OWVT6HCIcoBhwiAACAhmEKM1Rm8wWHCAAAoHGgMps3HMwIAAAAKw8RogWwvz0I+8cUFp6iLFZNpuxefTJZs2wtzt5y1GSFvTMsrzLrVK9lFhyl1q7I0sq1zDKhJssi1HSpeu81qczqUKXVEY7Paro1efXJlGpM9iHeTstRTnUy/9r3MqEmE1sfQydHZCDUZHmYocpMSO+89250xGR5irKuuAxrQk0ma5Yp9ZnzuYpVk3m1JJU9T/y2dUNx1/mDQwQAANAwsjzuCAsPcojiwCECAABoHMju5w05RAAAALDyECECAABoGOQQzR8cogWw3u5PJAju725WTp5W5Ti6TkJ0bPJ0WyVVi5IeXqK0TJ5Wdi95WiVbi7bRiORXL4FanRGSqqRqVbrDscvSHWpOYu0LKPXRlFtWGplA7FWYUNVjOqJMRV9kj/ctUaRk8rTTdMseGf738qG9ROtpydMd0XnHad6NTJ5eb/lvdF2sZS+BOjZ5+j3Cvt8rnzS3pGpyiPa8Q/Tmm2+GL33pS+HQoUPhwIEDxb8f/vCH634ZAACApYXSHUuQQ/Sxj30sfOpTnwrPPvts+PznPx8OHz4c2m0CUQAAALCCSdXmGL3xxhvhk5/8ZG1VmgEAAFalltkoSrTbB8Qx89DNpUuXwpUrV2b9MgAAAEtDHTlEu33+xYsXw927d8Nzzz0Xbt++HV5++eVw/Pjxmfbx4MGD8Oqrr4ZXXnkl+rXqYi57WUePHp3HywAAACwFtUR5dvH806dPF6kuFy5ceGIzZ+bevXvh1KlTtfdx4sSJcPDgweK/r169WjhEiyLKIfqt3/qt8AM/8APRL2LJ1fCX7O9shv1jyghPORarJltb3yytHJNqsnXfnkaW7kgd9ZlSmSnJTqKUVt6x/LNWmQ0md5dzVRpCqeaEyizZ7JQuRRKrGou179U9/iRplW7bEuVZPDWZLN0hOtdqMmV3SsIIxZdSmamMBCGaEyqzEKcmi7Cr0h1STSZVZko5Nnn/2C9K3+yLUZOJkkpZ4rddBm7duhUuX748keZijs2LL75YyiGK7WO0g3Tnzp3ieXsmh2i3W18j7w8AAADmk0OU7yLF5ciRIxP2kc0iOPPoY09EiMyD2w33798Pi2YRe6IAAAC7wRyaeZ9DdP36ddeZGe30XLt2bcfvvDr62BMOkb2R7/iO74jeArMQ2iKZ954oAABAFWzHqY5ziLIsCw8fPnzKvra2Vjy8oMexY8fkTs+NGzd2fM06+lgU0UnVFhmJJVGb2Uu6JwoAANAE3n777eJcwO185jOfCZ/97GcnVF7TsEDITm3q6GPP5BBZGMy2v8zjLPuwCMpHPvKRsChWfU8UAAD27pZZlYf9hH/++efDO++889Tj3Llzi357ez9CZGU4xj3NnRiV71gUTdwTXev2wno731FRFqsm6wq7pxxrR7SdpjLz1GTKnnREAaixuXjSXqhNXDc+pnDT1MJQYizZpMIlH4q+B0plJuyOSsazFcMTyjulSlPqu72qMotBLomh/95T8YxWMmnviOnuiTXRF4MZOv0MlJos1IM3wrZS2EXULDO6zmdWqczWxJpdV3XIhH3dUZTVoSZTNSbzQX9+SdVS/1e+jzRNwzPPPLNj2wM7pMKUiezU0ceeiRC99tpru3qR3T6vDmw/U6ncYvZEd9PH5uZmsXe7/QEAALAT5lpWVpnVeFr1vXv3Kh+hU0cfjXGIYqNDVZ9XlUXviZ4/f75476PHCy+8UGLUAACw6tSyZRYZVjxw4EDhtHjY91yZQ5br6KPRDtFXv/rV8Cu/8itEOCKxfdrt+7ZvvfXWoocEAADgcvLkyanH65iyeh59NNoh+rZv+7ZCYWWnVP/QD/1Q+LVf+7WwF1j0nqjJGm3vdvsDAABgJ+rYLovNQTpx4kShqh7/XrM8WkPJ6evuo/FbZlZ0zXJlPv/5z4ff//3fD9/+7d9e1Byxch57lVXYEwUAgL2HOTPVt8ziHKJjx44VAiFL9xg/YsYqVYx/19nZfPao0seIkQOlttsaeQ6RRYs+97nPFY8vfelLhSTdDi20CbB/P/zhD4cm0cQ90fVuf0Jl5inKYtVknX0bpZVjrX1+Hy2pMvOVFYlSma07KqnJclNbfXTFh7bdiijGFJUOp8mz0uqzZOC3zXuiDllX9O2oZ2Qdt0iVma5lJvrZo6TO+/RsW/ZWXC0zbymLs9WU+qovcjmGTi0zIYILdV2xNOa9C5Wiep9tp/2a6GM9QjU2ze4pytZbg8pqMqX+HaT9PXQwY/xzrly5UlRlOHv27JOqDKPveO+70RMcxfRhbWyLbXSAs/1/U25bv+Zf7Jlq93a+kEWMjC9+8YvhzJkzRb6RRY7ssMImbBHZfuY0JVnZPdGqfQAAAOwFzpw5U6rdzZs3K/exvfrDoqnpZ3UIn/jEJ8Ibb7xR7BOaomqUb2TJ2Itk1fdEAQBg71FsmVV8VD3HaNWozSEaYc7QKN/IHCQ72dq2lBaVb7TIPVEAAIC9klS96lTaMivjHH36058uHpZv9PrrrxdbabbFZPuJ3/u93xvmwbz3RAEAAKpg6T/Vq93XNpyVIMnHK5bOgTfffLNIlvrZn/3ZuTlFTcDOcTIn8dZ/+VfDezvtHROoY5OnO/tVovSkvbVPJE+LZOtkTZSdWMvLJ0p3RZJ0p+Pb07R8srVqG4tToqNgMCzfti+SLnt+Umjem5zDfNO/Eeab/m+Y7PFk5Wpj+Ljr2zcm2/cf+X30H6+79p7Th7Hp2Dc3/bYbvW6U/XFf2AedUrai74E/hxtDf332ssm11XOSoY1Np60xFF9sfaefPDKpWlWhUVVrvBGq+tsdkRDdEgnrXjmOruijm0YmVbfLJ0qr5Ol9nV7p5Gll/0+DQfjP3/yXxdlys8qT/Z7v+Z7w8cfvDy89+3ylfv6fe/8u/M63hPC7v/u7tY1tmZlphEjxsY99rHgAAADAJLWU3iieT5ioLJV/VlviNAAAANS8ZVbDA+boENl5AV/4wheqdgMAAADv0qTCrqtCLYkXX/nKV8LP/MzPhF/6pV+i3hkAAACsXg6Rqa/sDCLDksxMSWanTr700kvh4x//eB1jBAAAWCksWb6qymzr+eQQzc0hGjlD288gMuzEaosY3b17tziDaJXUZDvRXe+Ftc5wR0VZHWqywv6eyX5SUaIj8UpumH2faw6JUo6tOWqgrlCTtf1lmLcj1GezVpk59mQg1GRC3RQ6ovxJ27kWbVWiw1faqDIVqqTHKqBWhJ3O4tESc9Vy1GeeTZWuMAZClTZwpGBKkVaXHthTlCnVWFvZZUmPSXs3shRHN/Xt+4XKbM1TmdWgJlMllfriM1g/dZwjxLbZwlVmf/InfxIuX75cSOvtdOerV68W5/XYgYdWBqMJJT0AAACaSnHadOUIUW3DWQkqO0Tnzp17coKzlekwJ8jKXNjxRnYIox1eaDXPtm+p2cGGbKcBAADA0jhEIwfI6nqZE2RlLqxkx/attPEtNXOMbDvtU5/6VNWXBwAAWDpsi7RqhIcA0ZwdItsSszIWn/vc54qIkDk9O1GmDQAAwKqS15BDRC2zOTtER44cKQq5lsXk+V/84hepAQYAADA1QlS1lhkO0VwdIssRisEUZ6Pco1Wlu7YZumPqLE9RVoeaTCnK0v2iNtk+8QFaE8qpdb+2VOhO2vPuWj0qM2X3+kh8rVGSR6qvPEWZVJmJuW1vlq7NlrSECjD1XzOfm/Jl7yCVd4K22J9InY0Hpb7qi1pmfbEOPUWZUpnVlSCbRqjMlL0jFHmeXanG1ttxKjOlHPNUZnWoyQq7c6/tCXVc3eQ1KAvZMpuzQzSS2ZfFEqxHSdYAAAAAK1vcFQAAAKbn/5j0vlofEAMOEQAAQMOoox4ZOURx4BABAAAsYekOHKI4aqp3AAAAALB3IUK0ANrrvdAeU5mZbRZqMqUok2qyfUI1tu4rxPL1fZVVZnlLqcm6ceozr23q15xKsjiliFe3LB8IJdiwH6Wm81RziarN5kmEiopF4tqH2anPclGfK8/T0r9U67LXQSJq0Hm1z1pC8tUWKrOOsA88lZloW1c+SBJRx03XMquuMpM1zhzVmFKTKeXYmqgb2BUqs+6ab+849k4yH5VZHdecHKI4cIgAAAAahm2XVd4yq200qwEOEQAAQMMwZyaroQ8oDzlEAAAAsPIQIQIAAGhiLTNUZnMFh2gBtNZ6ob32dHCutW8ygbq1r1c5eVomUKvk6f37opKnpb0zmUCdd9f9tjJ5WozRSZRWydOxyGRrx56opGplF8njbkkPkVQtk60FtSRbqyTfmMRnmYAt7KL9LElEAnGaTF63VuZvZrSG/i21r5KTnfefi4RllU+itkWSiJImSU2lO7qtQem2Kkna62Nq2Q0ngVqV4lBJ1V7y9KKTqu2zgex+vuAQAQAANIyillkNfUB5yCECAACAlYcIEQAAQMOwSvfI7ucLDhEAAEDDyGqQ3Vd9/qqBQwQAANA4qqvMAknVUeAQLap0x7jKzFGOpY7yzEjWfZVDVDmOutRka/t9u6Moy7r76lGZtf0SIG4fqb/Ek2wQF2IebFZWmaVKCVeDQi72tuepz5TCLs1FCRmlBPNUabFKtVhVWsSNX7XVKjOndMegHZWU2RElV/rZ5LUXVUEKGXYdeKVIxPBCR5TdUOozTyHWbfuftY5Qk6n2MeU4OqJtjJpMllQSqjnY++AQAQAANAy2zOYPDhEAAEDDsEgmW2bzBYcIAACgkSqzan0QIYoDhwgAAGAJD2aEODiYEQAAAFYeIkQLIF3rhXRMZZauTaoikjVVm0x0vObXygrrazNTk2Xr7/HbO4oyqTLrvidKfZW3/Jpofh/+Es+FykyRtDecPoTar/eNqPeTOPa6fqkkouaWF4tPMl+ZE8RcpQO/fT6cHH3q2Ix2FqsyS0vbVf5F1cPupqGUagNHTabUWmqbI/b9eOo4NUa13pSiqi3UZ53WsLKarKsUYtLu1E9z7qdT1WTCbvfqCducNqLyGtZq1S23VQOHCAAAoGEUDlENfUB5cIgAAAAaBiqz+UMOEQAAAKw8RIgAAAAaBgczzh8cIgAAgKaRb51FBPMDh2gBpN1BSNee3ttNXJWZUIl0Re2r9W5p5VhdarJs/b1+P939pWzTVGN5W6jJknZpNVk0qsaZM5ZkMKk82xqLf30yVVctopZZ9B63UJl56rNk6CuHwkCo6YTKzFOUtQbiPYoaZ5lSpQ1bpdsrhU4mlGoxxCi4jJa4DgPn/WSiZlnlfJKpKjN/3G1HNTZNfebVMlN9aDWZ/xlsS3u/lM1odYV9rXz7Vj6fuIutA7UWYvqA8pBDBAAAACsPESIAAIAlLN3BllscOEQAAABNLN2BQzRXcIgAAAAaeTBjtRwg/KE4yCECAACAlYcI0QKwumXjKrN0fVKJkXTFr4M1X00WuuXteWeyvllh7wrFl6xDtr+8yqzz3rjaZEJllrgqM1HHLRZRzyvPJxUuuTeOYizVP1aqTpqyB1VXTahqgqcQk2oyX02XKlXacLLvXKjMcqEaa3VEDTqnbyNz1GpKlZVH10+rrtTph1ZptZpSx820lplomyZZ6ZplSlHWETXLlPqs046rZdZ21rhnM1od/zVToWDzakymc1KZ2XZX5S2zugazIuAQAQAANAxqmc0fHCIAAICGYQozVGbzhRwiAAAAWHmIEAEAADSQqgEeAkRx4BAtgKQzDElnLBHSy7lUJTq6fgJx3l0rbVfJ05lIns5ik6qdBGqVVJ2IpOqkta906Y6kptIduSjdkUQkVec1jCUT41DJ06lKwhZJyGEw2X/i2Ar6oo+emCsncVUlreYqmVXZRUmP1mByzjOR9N0eDkonZhdjdJKW6yqjMcwikqojS3qoMiJJKJ9UrUp0qIRozy7bqmRrYW+1h6Xtcv3UYE+H+fxKd1RcZ3Wt01UBhwgAAKBhoDKbPzhEAAAADQOV2fwhqRoAAABWHiJEAAAAy7hlRogoChwiAACAhsGW2fzBIVoErXxi5t0yHR1RjqLdjrLnrcl+8rZf5kPau+8Rfa+XtseqyWR7T2WmymhEkidCDeWozKL7lsoxx94VqrFBL+66OdfeSLy1otaVWIeJGmPfuQ0LhVASq/rp+2NsOcqkoSgXkrZEuRAxxlmqzDwl2DBC7bbVSfXXVGqyWlRmkaqxtOW/Zipe07v2sg+xrtQ6dNdtO98zBzNWff6qsRIO0cWLF8Pdu3fDc889F27fvh1efvnlcPz48dLPf/DgQTh//nzx7507d8K9e/fCuXPnovoAAACA5rL0DtHp06fD4cOHw4ULF57YzCEyp+bUqVM7Pt+coLNnzxbPP3DgQGG7detWePHFFwuH6MqVKzMdPwAArCYEeObLUqvMzHG5fPlyOHPmzFN2c27MUSqDRYa2O0PGkSNHCtvVq1fD9evXax83AACsNtm2bbNdPxb9JvYYSx0hunTpUuG8jDOymUOz07aXtbGHbbVt59ixY8W/FiEa/TcAAEAt1KAy222I6WLFNJPYPup4vTpYaofIojeeQ2RYxOfatWs7TvqhQ4fCjRs33OcbtvUGAACwDJyumGYS20cdr1cXS+0QWQK0it4cPHjQdXTGMadJbccZL730knzu5uZm8Rjx8OHD4t+klYdkXOjSdpQvqVCbtDtR9uAokGJVZiH1lTl521eCBcceqyZLU2F3VWai7lskee5/JDJHZRYbjlZKNc+eDDf8trHXTbZ31opYP0n6l2v46fb++kwcNVAiVD9JKuptCUVRjD1VyimhbspEnbTM+RyqvlOhBFP2OoitZRajJktrsMvrkJavTTatH3+9Ra4fuQ4nX9Pu33PbMquhj92kmeRjoSlzVixvtoyDEtNHHa9XJ0ubQ2TJ0NOwCM9ObXbajrM+pl0wyz969tlnnzxeeOGFXb8eAACsDuYjVM0hit1yu1QizaTOPup4vTpZWodo1ltx9njttdeeSrYex6T577zzzpPHW2+9NddxAgDA3iSv6RHD9evXizSRaWkmdfZRx+vVydI6RNMcFaNKdOjEiROFZ7tT/tHa2lp45plnnnoAAADMiyzLinSN7Y/tqRzjaSYHDx6slGYS00cdr7cyOUSWaBWTtGwTaB6l8ji3Y/2Waec5Qxb5mffeJgAArFjpjqq1zEIIb7/9dpGysZ3PfOYz4bOf/WztaSYPIvqYdVrL0jlE41L3WGxClUNlE3306NGo/uyARkuiHj/XCAAAoInFXZ9//vnw5S9/eWL3AvaYQ1SVkydPTg25mbSvLJYJb2ckjDtDZo+OFqXOZmXqKEU85VnRVux0xqjPhGostNeiapYFUUPMrS0W01aoybbat2ZWy0zhzXguXlPZk4j2ar6TtlB8pY/9vpVyzF0Tal2pdSiUU043Sbu8Qmiq6qddvr2qfZWIGmdKleXW/hJ9Z6IOmVJxeQox1TarQU1mpE571YdUn6n2nuIrYl637HG1zLxrH7NOCrvq21m3yZxqmdVV3DVN01LpGgdqSDOJ6WOWaS27ZWlziEbbWybrG5/Y0enSZQ9UtPbWhxcZWsRFAwAAmCf37t3b0Ymps486Xi+WpY4QmcNjic+j8hsj7L/thOnxybacpfGtOkv6soOjrC/bMtvuBI3+BgAA0LRq97FbbgdqSDOJ6aPutJaqLLVDZJjjY8eCmzMzOhbcnBhPIWYXZzzj3bbVzPGxrTGP7Y4WAABAHexGNu/1Me80k5MRfdSZ1lIHS+8QGWWToG/evFl7YjcAAMCiVGaxaSaXL18uojPbd1Bi0kxi+qjj9epkJRyixuElVSdp+SRXZY8gV0nVsaRiCaWTSbtJGpdsrMpxeO3Tmt5PFpHFKMct3mfuzMnUOZz39Yxdb4lq79nErVna49JJvQRdmUBbh93Pwd1FAnE+94RPd65mOO6ZXofYpHK1rmLWZ2QS+6JVZvNOMzkW0Ufs680aHCIAAACoJc0kto+YtrMGhwgAAKBh1CW7T+acZjIi5ry+ppzth0MEAADQMKwCfFZxz8z62I1DtKrgEAEAADSMRajMVp2lPpgRAAAAoAxEiBaBlTxQZQ8qkAvVT4wCKZeqMZbKXJFKNaFsE92oa++tldpW5AzW9m7LVMyq7zrUV6oExlCU6JglqhSHohb1WUOu5VIfzFjXYFYEvuUAAAAaSF7RpcEhigOHCAAAoGGYwmzepTtWHXKIAAAAYOUhQgQAALCE5xBVff6qgUMEAADQMOwMIXtU7KSu4awEOEQLkw/U322S+50mmSi85LYduPZc2GFGiPlW10ehrr1aK7XgJj7Uo5zKZ6jAiulbtY0dX7YARVnMONIZvv+mXMumQoRo/pBDBAAAACsPESIAAIBl3DJDeB8FDhEAAEDDYMts/uAQAQAANPIcoqrFXWsbzkpADhEAAACsPESIFub6j9k81U8mAp7KXoP6KPoHhVI9Zf3JvpWCLRH2vPzyrGFK3n1NNZZh+bYRc7Jlr67gi1ESSmLXm1KqeeZMqH6kPa2sKspF33XZy44j1q7a1qVIi5nZWtRkM74OUYoyta5i1mc+z8IdVUt3ECKKAYcIAABgD/xujgV3KA4cIgAAgIZh0Z2MCNFcwSECAABo4vm9FbOiqxaHXTVIqgYAAICVhwgRAABA46ieVE0WURw4RE3JlvNim4NhnOpn4KuYEs+uVEmDTb+P9oZrz9vrvt1RYCVSweXbM2GfZVjTU5Opsahxh8j36bVPhhtR10ddT/faF/30I9aVWofiZut0kw/8q5YPW3GKokH59pnqO08rK6qGom+lBBsKdZPbdw0KOyNJ/OszdNq3Uv/aq7HIvp15SUXfer7915TX07n2MeuksKu+nXWbD+bjZGQ15BBVff6qgUMEAACwlEnVEAM5RAAAALDyECECAABoGBbd4WDG+YJDBAAA0DA4h2j+4BAtgHyYFI/tJIPypTuiEmWLJMCe08ekbZo9l0m7Itk6aZeyTUOd0ur1k0T2HV+6w0mqFonP+fCx37mYK3cO5XzHXbcg2/fLryuZbC2SZcfW9pYtrSfJNcKeiYTgoUi4HWYiUdrpR/Uda/eSlmMSsHeDSoiOIRN9eHb93v35TsW6SlvD0tc+ev3IdegkVQ/nlFSd5CFLqp1VnddwrVcJcogAAABg5SFCBAAAsIRbZsju48AhAgAAaBhb7hDlXecJDhEAAEDDMFeGCNF8IYcIAAAAVh4iRIvAVDiDp1UNeW8yNJp0lZqsI+y+QioZ9iurlZLeN/zXTH3VRkgdlZlj2w2uyqyuvrPy5TWUmkypz5IIu5rv6OvmXHu5VsT6CX2hXuyJX59emQ5VSqHvX7dM2YVabTholy71oPvw2w+cvgei71i7pyhT5T/yUL60yDQ1WeJUUI8pWxKLHIewt1JRPkfNoXN90uEgal2lwh463udkXqU7MlRmcwaHCAAAoJFJ1RUdomLLrJ7jGlYBHCIAAICGsZVSjUM0T8ghAgAAgJWHCBEAAEDjqC6733o+cY+y4BABAAA0jFpKdyC7jwKHaAHk/VbIW2O1zLrOwu/5aovQ8VU/SXvTb992VFktX6mWStWYUOzUoO5Syi5VV8yrW5anQnkXS9YvP5ZBpJqs/598e+/RhC3t+Qo2ZU964jV7m+XtPaFIU+twqNd3adVYTXZPIRajGptqd9RNsWqyvlK85REqM2FX7dMIdZdqq/qepcosTVUtM2WfnNuWWj+iD7WuEsee9aselliOrQyiaq/FOURx4BABAAA08qRq8aujNPNx3pYFNhcBAABg5SFCBAAAsITnEFV9/qqBQwQAANAwtjKISKqeJzhEAAAAjaN6DlFV2f6qgUO0ALLNTsiSMfWGo35IW/6HIWn7datCWyjE2p3yijShJksi7VFEqsy8WmZe7bQ6x+LVMktiVWaOmkzZk1g1WV9cz16vvH2zF1WzLNsQyqnNTinbVNWYo1QzBr1OaXu/77fti1qAUiEWUcusN2yXrlmm+skia5bF4qrMRCShLe5Bqq5a1XEUdnF9kjRCNSeUaol4P6qWWeKtZafuJCwHOEQAAABLW7oDyoJDBAAA0DjqOqkayoJDBAAA0DC2UqrJIZonnEMEAAAAKw8RogWQ9dohS572RRMnqTpvC+9e2JOWnxSbpI7f69lq9JzzbPKXTaYSlrvDqOTkvLVefiAq2VqNReCOxXmPRdveN/yhyETpSXu64feRbIrE7I3HUfawMblWcpEsmm+KJF+VKO3YhyIZeij6UO0H/fL2gUiU7Qt7T/Tdc5Kq+yJ5WpXoGGSt0snW6jf9TEt3iNccOqVFjHYqkq2dsdRR/sNIRV0v733qMh8iqVq0D05y9nBOSdWW/1N9y4wcohhwiAAAABpGnmchy9kymyc4RAAAAI2sZUaEaJ6QQwQAAAArDxEiAACAxkG1+3mDQwQAANAwig2zvGJx14rPXzVWwiG6ePFiuHv3bnjuuefC7du3w8svvxyOHz9eqc8XX3wx3Lx5c1fPzTa7IQs7q8yStvh1kPoKqSTti/ZJOeXZLkgjVGZKlZUPRMmIdte1y7IjXh9CZZZEqszCYPI1EzFuZdcqs40ZqsnEXG1OrpVcdJFv+nMoy3FsdkvZjIGw96W9U1o5pkp3xKjJlKJMtZWlO4TSqu+ozzKR9lFHuQwjcfJKnFtEQUck9apSJHWVF6mqmkuEIk3aVVkQxz6/0h3Vc4hs1UB5lt4hOn36dDh8+HC4cOHCE5s5RPfu3QunTp3adZ+3bt2qcZQAAACwSJbaITKn5fLlyyHPn/aSzTmyCM9uHCLr88aNGzWOEgAAYFJ2n1eV3bNlFsVSq8wuXboUjhw5MmEf2a5evRrd5+uvvx5eeeWVWsYHAAAwrbhrlf8hu49jqR2i69evh0OHDrl/O3DgQLh27Vp0LtK5c+dqGh0AAICgiBBVfOAQRbHUW2Z37twJx44dc/928ODBqK0v2yoz58ocqbJsbm4WjxEPHz4s/VwAAFj1gxk5qXqeLK1D9ODBg6l/N8dmpzbjW2XbE7PLcP78+fBzP/dzE/bBRjcMxmsEOTVzlPJBqS1ypT4LvurJb1sTjqIsVSozoSZT9pBOyqHytFXL+0nEGL33E6sy89Rkhb2/WV1N9kjYHwsF3+PJNZRvCDXZ4zXXPtwQyjHH7tmmq8mEXSnEepPtN0VbZVcKsc3BZPu+UFmpPlR7T30WW7NMxQDU2vfuH55Sa9q4O879qhhjRD00+X6EXY0xpq20p+XbD3tVzwaCprLUW2Z1sdutMnvOO++88+Tx1ltvzWR8AACwXJgYqOqWmW27QXmWNkK009ZW2ejQbrbKRqytrRUPAACAOKh2P28a7RDZ+UF2XlBZLC/IEqVVIvV2rN8y7XazVQYAAFAFyx/aq7L7izUchhzbhwU5Xn311UIFvtuDlxvtENkkVMGiOsqhssk7evTo1OebLN8iRHYQ43ZGydgjuzlMu4kgAQAALBOnazgMOaaPEydOFMGQ0Xd2lWNxGu0QVeXkyZNTlWQ2wdMwL9PzNEcnVds5RwAAALPKIarYS9hrhyHfiuzjypUrT1Tl9rwqLLVDZJ6jTZBFg7ZHcOx8IkNJ8mfNcNNUZq0S9XjEYhYKD0UaBpWUZ0X7TLymsCddp1bW0K+1lrd81U+Q6jPR3mur1GdKTSZIBs7YlZpMvM+kJ+qK9XrVa5NFqMmM7NHkRz8TSrDhY2X38+MGG2uV1WSeaszY3FwrrRzbUH04qrGivbD3hpNrqOfUIJumyhoou6OoknXCQj14Gq6WuKe0xT1I1Wbz7MN0WIvKbJbEqMz6c1KZ5bXkEGWNOwz5+A7bWXX0sVuWWmVmDo9NnMnfxz1N8yrHt7ksRGePnYjJawIAAFjIwYxjUZa9cBjy9ZoPVI5hqR0iwxwfS8o6e/ZskaRl21328DxMm+xpidYWbbKo06jkh4XvxvOLAAAAVpE7d+48yefZ7WHIdfSxW5Z6y2zEmTNnSrW7efPm1L/b3uVuCsICAADMe8us2DTLsokqCbM4EuZBDYch132gcixLHyECAADYk6U78mHFRxbefvvt8Oyzzz71GE8jgRWKEAEAAOwtLDpUNUKUheeffz58+ctffso6iwODD9RwGHJdByrvFhyiBVDUMhtTqSRJ+VpmdeApz6aqz7I8Tn3mqrLEa7bbUWqyJEZllvhB0CRWzuq8H1d5Nu19RqjMwoa4Dpv9ymoypSgbfmPdbTt0VGPT6pN5dqkmi1WZiTpknqJMqcke98VrCuXYxqC6yiymlplScImPYDSp031LjK8lVGaqltnAsQ+9FywiIPNXmUXXOPNUZv29VcssTdPwzDPP7JnDkGfdxzRwiAAAAJpGPv/SHbcXfBhyXX3sFnKIAAAAmphDVPEx74MZT548WajEdnsYcl197BYcIgAAgMaRb8sjqvKYHydOnChOmh7P9Yk5DLmOPnYLDhEAAAA04jDk2D5GjByoKgcnk0O0AHqbazIpc9GokhZJJspRDEWC4WDSnohk4yCSqoNKnk7TcjZRpmAqESVKYpOqQ0+035xMoM57/jhyUdEj3yifPK0SqFXydP+RsD/2k7B7Tj+bom9VikOV3VB2L1FaJVU/GrTjkqq90h2OzeiLhOBBJuwRSdV1HTqceEnVIqlYle5Q77PjJGcPWiJJXHw6ByIJW6FKgMyKXl98vuumuOAVIzxzPqnaMKfFDkG2w5BHleqnHYbsHcIY04e1sS02iyqN/r8lelu/sfVGcYgAAAAaxlYOUDLXpOomHYZctg+LHNUFDhEAAEDjqCFCNOccor0ODhEAAEDjyBey5bXKkFQNAAAAKw8RIgAAgIaWd63aB5QHh2gB9Da6YXO44KkXR/WnuSgvkQ1Kq8mMZLAxaewLlVVHlOhIxVjardIqs2hkKZJh+bbqffb8ucp7kzetfFMojTaFmuyxr9YaPhYqM0f1VYeaTCnKZqkmK+yOosyzFX0LlZmnJjN6zmelJ1RjmxElOoy+04/aJVHZIKqkhxJrpSWVZ0ZHlA9qiRdd80p3iPc+VHPVSqNKfcw7S2bTKeVSN//m3zxde2z35OFf/at/XVNfyw9bZgAAAA3ij/7o374b3akS4dl6/r//9/+uxpEtNzhEAAAADeK7vuu73pXcZxUVakl44YUXah7d8oJDBAAA0DD+43/8DxWiRFvPe/Dg/gxGtrzgEAEAADSMb/7mb373K3oY6RSNokNpePbZZ2c4wuUDhwgAAKCBfOMbX3/3v2IdIhMyiDo/IEFltgA2e2uhk81v6nNXyeIrNlpCPZOKul2qnlfq1Tjr+Uq1pCvqobXT8vKZpCbfPle1zJwb0iArrRorEG8z25hUreSbvkIqE/ahqlkm1GcDp30dajKlKItVkz3q+X1L5Zij2nwcqSbT9sm1talqkwl7Xy0Jp/0wUmUWi/dJEeXGQk98rDpCltZ3PptrQqk2FC+qvvqjXAJZD6563bONwfxqp+3fvz984QtfCD/5kz/5bk5RUio6ZLXAul3/swUaIkQAAAAN5Sd+4iciXMKtNp/4xCdmOqZlBYcIAACgobRarfCbv/mb78YI8x2jQ7/zO78TEnWwFEwFhwgAAKDB/PAP/3AJGf6WzP6jH/3oHEe2XOAQAQAANBiL+Ny6dXOKDH/L/sd//EcLGN3ygEMEAADQcD7ykY9MiRJtRYe+8zu/cwEjWx5QmS2AjV4ntDNfLTML8jwtr7YQ9YVyR2ljpMIehpPqs6QjankJCU7SErKsNMK1VwWdstiCUZOmfCj6Hog57LdKK8eUmizb7JZWjW3Zy6vM6lCTKeVYXWqyR8LuKcqUauyxWLOemkzVJ+uJddIbxqrMJm2DfLZlOr0RtsVSbonPT0e07zpTLmuZqTlRH03Rj2dXbetgc44qs3H+5E++Gj784Q+/uxqSp6JDf/Znf7awcS0LRIgAAAD2AN/6rd+6LUo02j6zH45p+MAHPrDo4e15cIgAAAD2CPfv39vmDOVjBzhCFdgyAwAA2CMcOHDg3VjG1n7rr/7qrxYHOEJ1cIjmSJ6/680PJk9sHqSTtr5jM3oit6ab+va2Y28lfhJEW9hTkVyTvvueytgTlRzQUjlEeQ05RGGGOUSirbCrXKlsc9KeiSSVbNO3D4R92PMHM3Dsvb7fttf31+HmQJzs7ORYqLyLRzKfxzWHR2INbTjJOPLkaZEnp+y9mBwiYR80PIdIHLAd/BnU7b3UHX3ytJgscQ9SH87c+cDlib9mg7BniX/ifu6cxP+Nd22j+/kisLIca2tb+Xc//uM/vrBxLBs4RHPk61/fCmv+yL/+7UUPBQAAKt7PF1U81cpyfO1rXwtpmhYHN0I9JPki3dwVI8uyYhG/733va+RJog8fPgwvvPBCeOutt8Izzzyz6OE0FuapHMxTOZinvTVP9pVpztAHP/jBwiGB5YEI0RyxD8+HPvSh0HTsZsONeWeYp3IwT+VgnvbOPC0qMgSzBfcWAAAAVh4cIgAAAFh5cIjgCaZa+MxnPvNEvQA+zFM5mKdyME/lYJ5g1pBUDQAAACsPESIAAABYeXCIAAAAYOXBIQIAAICVB4cIAAAAVh4OZoRw8eLFcPfu3fDcc8+F27dvh5dffjkcP368Up8vvvhiuHnzZlgmqs7TgwcPwvnz54t/79y5E+7duxfOnTtXea734nqZxZprGqyXcnD/gcZgKjNYXU6dOpVfuHDhKduxY8fyS5cuVepz2ZZW1Xm6f/9+0Yf9O+LmzZvFPB0/fjxfpfUyizXXNFgv5eD+A02CVbPCjG6wZe1l+zxy5MhS3ZDqmKczZ8489eU2wr4MrI9r167lqzAPs1hzTYP1Ug7uP9A0yCFaYS5duhSOHDkyYR/Zrl69Gt3n66+/Hl555ZWwTNQxT9bGwvjjHDt2rPj3ypUrYRXmYRZrrmmwXsrB/QeaBg7RCnP9+vVw6NAh928HDhwI165di84FsByHZaOOebLnWw6I93zD+9syzkPda66JsF7Kwf0HmgYO0QpjiZoHDx50/2b2GzdulO7r1q1bxc1tdMNeJuqYJ7u5379/350346WXXgpNp455qHPNNRXWSzm4/0DTwCFaUUy5Mg27sezUZjxUvWzql1nMk7dtYH2cOnUqLPs8zHoumwDrpRzcf6CJILuHyhCq3v2WgT0sH4RftrATrBcf7j9QF0SIVpSdbqhlf50te6i6rnnyOHHiRPGLfy/8sq1jHmY5l02B9VIO7j/QRIgQ7WEOHz4clVxp+/KWm6ASGbdj/ZZpZ6HqCxcuhCbThHnyvtzsV+1e3/qoOg9199F0WC/lWKb7D+wdcIj2MHaqaxXsV5VyFOwX2tGjR6c+32Sx9gvt9OnTT9lHyZAju92wFvkLbtHzNM7Zs2eLpNgzZ86EvUQd81D3XDYR1ks5VuX+A3sHHKIV5uTJk1OVHHaE/jQsdO+F7+1GZDcqC+8vA1XnaTuXL18uShSMf7mZvem//uuYhzrnsqmwXsrB/QeaBjlEK4yF4e3GMb5fb4mb2w+BW3Xqmidrb314v/T3Qv5MHfOwCmuO9VKOVVgLsLdI7LjqRQ8CFntTsr367fvw9svMfmWN//qyXJwyW1DWp4Wzl2lpVZ0nO3PF2ttNfhS+H30R2N+8fpZ1vcT0sVdhvZSD+w80CRwiKF1t2koJjBKOPSyMb38bHblvR/BbHsCyhK6rzJPdzO2LTGGVub0yBsu6Xqh2vzrrZSe4/0BTwCECAACAlYccIgAAAFh5cIgAAABg5cEhAgAAgJUHhwgAAABWHhwiAAAAWHlwiAAAAGDlwSECAACAlQeHCAAAAFYeHCIAmAnTTloGAGgaOEQAUDtWRmFUpBMAYC+AQwQAtWP1o06ePPmkKKkV3LRaVEmSFI9RAU+cJgBoCu1FDwAAlm+rzIpwjqq0279XrlwpHKP3v//9RXVzVaATAGBRECECgNqjQxb9GWe7gwQA0DRwiACgVq5evRqOHz++6GEAAESBQwQAtXHr1q1w5MiRRQ8DACAaHCIAKLh48eJTic9nz5598rfDhw8XNvvX2w7babsMAKDpJHme54seBAA0B3NoTDZ///79J/k+5izdvn27cHimYQ6TtVOYU2URpJs3b9Y+bgCAKhAhAoCnMKfHnBaTyhumDrt79+6OzhC5QwCwl0F2DwATmEzeoj0WKbKIz4ULF3Z8zuuvv16qHQBAEyFCBAAT2FlBo3ygl156acf2FkWy84fseVXhsEYAWAQ4RADgcvTo0cLB2Z5crXjjjTfCK6+8UptSDQBg3uAQAYDcArPk53v37u2oHLNo0qlTp2qJDtURZQIAiAWHCAAmsKjQuXPnCpXZm2++ObVYq22XbS/VUQVzrDjJGgAWAQ4RADzFaIts5JiY4syiP6Y6szyhccxZGinSpmGO0zQsCmVKNSJEALAIOIcIAJ44QubcmONizpCdQ2SYE2QHNo7sx44dC6+99toTh8n+Nu1cIXveq6++WvQzyg8yeb5FlWw7zv5+48aNJw4TtyQAWAQ4RACwa8zBOX/+fCHTBwDYy7BlBgCVEq8p1QEAywARIgDYNTuV6gAA2CsQIQKAXWGqM8snAgBYBnCIAGBXUNkeAJYJtswAYNcJ1SbJBwBYBnCIAAAAYOVhywwAAABWHhwiAAAAWHlwiAAAAGDlwSECAACAlQeHCAAAAFYeHCIAAABYeXCIAAAAYOXBIQIAAICw6vz/yjusol51YX0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Calculate nonlinear density contrast\n",
    "delta = rel[\"rho_n_fromHam\"] / sol.rho(grid['t']) - 1\n",
    "\n",
    "# and plot\n",
    "plt.rcParams.update({\n",
    "    \"text.usetex\": True,\n",
    "    \"font.family\": \"lmodern\",\n",
    "    'font.size': 16.0,\n",
    "})\n",
    "plt.figure()\n",
    "plt.pcolor(fd.xarray/L, fd.yarray/L, delta[:, :, int(N/4)], cmap='inferno')\n",
    "plt.colorbar(extend='both')\n",
    "plt.gca().set_aspect(\"equal\")\n",
    "plt.title(r\"Initial $\\delta$ at z/L = 0.25\")\n",
    "plt.xlabel(\"x/L\")\n",
    "plt.ylabel(\"y/L\")"
   ]
  },
  {
   "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.13.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}