{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 1: A simple example for OBS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example we generate plane wave seismograms and P receiver functions for an ocean-bottom seismic stations. This example follows that of Figure 2 in [Audet (2016)](#references). \n",
    "\n",
    "Start by importing the necessary modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from obspy.core import Stream\n",
    "from obspy.signal.rotate import rotate_ne_rt\n",
    "from telewavesim import utils as ut\n",
    "from telewavesim import wiggle as wg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select the model file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "modfile = '../models/model_Audet2016.txt'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select the type of incident wave - options are `'P'`, `'SV'`, `'SH'`, or `'Si'`, which is an isotropic S-wave source:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "wvtype = 'P'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we use global variables to define the desired time series. Be careful to use a total length of time large enough to avoid wrap around effects. Sometimes if you see signals arriving at aberrant times, try with either a greater number of samples or higher `dt`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "npts = 3000 # Number of samples\n",
    "dt = 0.01   # Sample distance in seconds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we choose to model seismograms for an ocean-bottom station, so we specify the deployment depth. \n",
    "\n",
    "<p class=\"alert alert-warning\">\n",
    "  <strong>Warning!</strong> This one is required for OBS simulations. Otherwise an Exception will be raised and the code will stop.\n",
    "</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "dp = 2000. # Deployment depth below sea level in meters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The parameters (P-wave velocity and density) of the water layer have default values of 1500 m/s and 1027 kg/m^3, respectively, but these values can be changed here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = 1.5      # P-wave velocity in salt water (km/s)\n",
    "rhof = 1027. # Density of salt water (kg/m^3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now specify the parameters of the incident wavefield in terms of a horizontal slowness and back-azimuth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "slow = 0.06 # Horizontal slowness (or ray parameter) in s/km \n",
    "baz = 0.    # Back-azimuth direction in degrees (has no influence if model is isotropic)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read the model parameters from the model file and return a Model object. Up to here, the steps could have been performed in no particular order, except the name of the file that needs to be defined before the call to `read_model()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2800.0, 2800.0, 3200.0]\n"
     ]
    }
   ],
   "source": [
    "model = ut.read_model(modfile)\n",
    "print(list(model.rho))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following step is useful for plotting the seismograms with respect to the predicted propagation time through the model (from lowermost interface up to the seafloor). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Predicted propagation time from model:  1.1 sec\n"
     ]
    }
   ],
   "source": [
    "t1 = ut.calc_ttime(model, slow, wvtype=wvtype)\n",
    "print('Predicted propagation time from model: {0:4.1f} sec'.format(t1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we generate the plane wave seismograms by calling the `plane` module from the `utils` module and storing the output in a 3-component stream oriented in North, East and Vertical directions. This is also where we specify that we are interested in OBS synthetic data, i.e., including effects from an overlying water layer. The `obs` boolean variable defaults to `False`, so use the following call with argument `obs=True`.\n",
    "\n",
    "<p class=\"alert alert-info\">This one line below is the main time-consuming part of the code!</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "trxyz = ut.run_plane(model, slow, npts, dt, baz=baz, wvtype=wvtype,\n",
    "                     obs=True, dp=dp, c=c, rhof=rhof)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could also calculate the transfer functions between components - here since the incident wave is a P-wave, we can get the radial and transverse transfer functions as an `obspy` `Stream` object. The option `pvh=False` tells the code to use the regular Z-R-T geometry. The option `pvh=True` rotates the components to the P-SV-SH wave modes using the seismic velocity of the topmost solid layer. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "tfs = ut.tf_from_xyz(trxyz, pvh=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we can extract individual components as `obspy` `Traces` and rotate them to the Radial and Transverse directions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "ntr = trxyz[0] # North component\n",
    "etr = trxyz[1] # East component\n",
    "ztr = trxyz[2] # Vertical component\n",
    "\n",
    "# Copy to radial and transverse\n",
    "rtr = ntr.copy() # Radial component\n",
    "ttr = etr.copy() # Transverse component\n",
    "\n",
    "# Rotate to radial and transverse\n",
    "rtr.data, ttr.data = rotate_ne_rt(ntr.data, etr.data, baz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For plotting purposes, we create a new `obspy` `Stream` object that includes one transfer function and two Green's functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "strf = Stream(traces=[tfs[0],ztr,rtr])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now plot the result as wiggles, using the format displayed in the paper by Audet (GJI, 2016). First you will have to define frequency corners (defaults are `f1 = 0.1` and `f2 = 1.0` for `'P'` wavetype)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaIAAAEYCAYAAAAeWvJ8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeXwV1fn48c/JvhKysgUIKHs2EgJEZJPViiCICrUqVkTcv/oT9ypVq61SF9SqdUOtFaqtiKLVilDApRIwKDtEAoQ1JGTf7vL8/pjkGtYkkDAJPO/X676SuXfumWfm3jvPzJkz5xgRQSmllLKLl90BKKWUOrtpIlJKKWUrTURKKaVspYlIKaWUrTQRKaWUspWPXQuOioqSuLg4uxavlFKqCaxevfqgiEQ35D22JaK4uDgyMjLsWrxSSqkmYIzZ0dD3aNWcUkopW2kiUkopZStNREoppWxl2zUipVTjczgc5OTkUFFRYXco6gwXEBBAbGwsvr6+p1yWJiKlziA5OTmEhoYSFxeHMcbucNQZSkTIy8sjJyeHLl26nHJ5WjWn1BmkoqKCyMhITUKqSRljiIyMbLQz73olImPMWGPMZmPMNmPMvcd4fZgxptAYk1n9eKhRolNKNZgmIXU6NOb3rM6qOWOMN/AiMArIAVYZYxaJyIYjZl0hIuMaLTKllFJnhfqcEfUHtonIzyJSBcwHJjRtWEqplsrb25vk5GTi4+O5+OKLKSgoaND7Z8+ezZw5cwB46KGH+PLLL084/7Rp0/jggw+Oen7evHns2bOnQcs+Ge+//z69evVi+PDhtpZxpIKCAv7yl794pvfs2cPkyZMbrfzGVJ9E1AHYVWs6p/q5I6UbY9YaYz4zxvRplOiUUi1OYGAgmZmZrFu3joiICF588cWTLuuRRx5h5MiRJ/XeEyUil8t10jEd6fXXX+cvf/kLS5curdf8x1p2Q8uojyMTUfv27Y+ZsJuD+iSiY1UEHjms6xqgs4gkAc8DC49ZkDEzjDEZxpiM3NzchkWqlGpx0tPT2b17NwAlJSWMGDGClJQUEhIS+Oijjzzz/eEPf6BHjx6MHDmSzZs3e56vfbbzyCOPkJaWRnx8PDNmzOBEo0t/8MEHZGRkcOWVV5KcnEx5eTlxcXE88sgjnH/++bz//vu8+uqrpKWlkZSUxKWXXkpZWZlnmbfddhvnnXceXbt29Sx/7969DBkyxHO2t2LFCh555BFWrlzJzJkzmTVrFi6Xi1mzZpGWlkZiYiKvvPIKAMuWLWP48OH8+te/JiEh4bBYjyxj3rx53HLLLZ7Xx40bx7JlywAICQnhgQceICkpiYEDB7J//34A9u/fz8SJE0lKSiIpKYlvvvmGe++9l6ysLJKTk5k1axbZ2dnEx8cDVqOWa6+9loSEBPr27etJgPPmzWPSpEmMHTuWbt26cffddzfwEz9JInLCB5AOfF5r+j7gvjrekw1EnWie1NRUUUo1rg0bNtgdggQHB4uIiNPplMmTJ8tnn30mIiIOh0MKCwtFRCQ3N1fOOecccbvdkpGRIfHx8VJaWiqFhYVyzjnnyFNPPSUiItdcc428//77IiKSl5fnWcZvfvMbWbRo0VHz1DZ06FBZtWqVZ7pz587ypz/9yTN98OBBz/8PPPCAzJ0711Pe5MmTxeVyyfr16+Wcc84REZE5c+bIY4895lm3oqKio5bzyiuvyKOPPioiIhUVFZKamio///yzLF26VIKCguTnn38+5jarXcabb74pN998s+e1iy66SJYuXSoiIoBnvWfNmuVZ1uWXXy7PPPOMJ7aCggLZvn279OnTx1NO7ek5c+bItGnTRERk48aN0rFjRykvL5c333xTunTpIgUFBVJeXi6dOnWSnTt3HjNmkWN/34AMqSOvHPmoz31Eq4BuxpguwG5gCvDr2jMYY9oC+0VEjDH9sc608hojUSqlTt6UKVPYt29fo5XXtm1b5s+ff8J5ysvLSU5OJjs7m9TUVEaNGgVYB733338/y5cvx8vLi927d7N//35WrFjBxIkTCQoKAmD8+PHHLHfp0qU8+eSTlJWVkZ+fT58+fbj44osbFP8VV1zh+X/dunU8+OCDFBQUUFJSwpgxYzyvXXLJJXh5edG7d2/PWUdaWhq//e1vcTgcXHLJJSQnJx9V/hdffMGPP/7oOYsqLCxk69at+Pn50b9//1O+58bPz49x46w2YampqfznP/8B4KuvvuLtt98GrGt0YWFhHDp06LjlrFy5kltvvRWAnj170rlzZ7Zs2QLAiBEjCAsLA6B3797s2LGDjh07nlLcdakzEYmI0xhzC/A54A28ISLrjTEzq19/GZgM3GiMcQLlwJTqzKiUslFdSaMp1FwjKiwsZNy4cbz44ovcdtttvPvuu+Tm5rJ69Wp8fX2Ji4vz3IdSV1PgiooKbrrpJjIyMujYsSOzZ88+qXtYgoODPf9PmzaNhQsXkpSUxLx58zzVXwD+/v6e/2t2ZUOGDGH58uUsXryYq666ilmzZnH11VcfVr6I8Pzzzx+W1MCqmqu97BPx8fHB7XZ7pmuvp6+vr2dbeXt743Q661XmkU60e6697qeyjIao131EIvKpiHQXkXNE5A/Vz71cnYQQkRdEpI+IJInIQBH5pimDVko1f2FhYcydO5c5c+bgcDgoLCwkJiYGX19fli5dyo4d1mgBQ4YM4cMPP6S8vJzi4mI+/vjjo8qq2RlHRUVRUlJSr4vuoaGhFBcXH/f14uJi2rVrh8Ph4N13362zvB07dhATE8P111/Pddddx5o1a46aZ8yYMbz00ks4HA4AtmzZQmlpaZ1l1xYXF0dmZiZut5tdu3bx/fff1/meESNG8NJLLwFWY4iioqITrv+QIUM867xlyxZ27txJjx49GhRnY7Kti5+q7duhXz/w84PRo+H666HDsRrjKaVaqr59+5KUlMT8+fO58sorufjii+nXrx/Jycn07NkTgJSUFK644gqSk5Pp3LkzgwcPPqqc1q1bc/3115OQkEBcXBxpaWl1LnvatGnMnDmTwMBAvv3226Nef/TRRxkwYACdO3cmISHhhEkLrLOap556Cl9fX0JCQjxVYbVNnz6d7OxsUlJSEBGio6NZuPCYbbeOa9CgQXTp0oWEhATi4+NJSUmp8z3PPfccM2bM4PXXX8fb25uXXnqJ9PR0Bg0aRHx8PBdeeCE333yzZ/6bbrqJmTNnkpCQgI+PD/PmzTvsTOh0M3bVoPUzRjKGD4eqKvjmGwgIgHvugVmzoLquWCnVMBs3bqRXr152h6HOEsf6vhljVotIv4aUY1tfc1sAvvoKVq6ErVvh4oth9mxITIRadbVKKaXObLYlosNOgs85BxYsgKVLwRgYPhxmzoSiIrvCU0opdZo0r963hw2DtWvhrrvg1VehTx/49FO7o1JKKdWEmlciAuv60FNPwbffQlgYXHQRXHWVnh0ppdQZqvklohr9+8OaNfDww/Dee5CeDtu22R2VUkqpRtZ8ExFYTbtnz4b//Af274eBA+G77+yOSimlVCNq3omoxvDhVgIKD7f+b2C7fKXU6dNchoFoKtOnT2fDhiOHYzt1mzZtIjk5mb59+5KVldVo5T777LOeDl0BfvWrXzX4M2lqLSMRAZx7rnW/UWIiTJoEp9C1PADFxfDJJzB3Lrz2Gmza1DhxKnWWay7DQNTWmN3UvPbaa/Tu3fuUyzlyOIiFCxcyYcIEfvjhB84555xTLr/GkYno008/pXXr1o1WfmNoOYkIIDraauI9bhzccot1A2ytPpnq5aef4MYboX17696l22+3enXo1QtGjIAff2ya2JU6C9k1DATAsGHDuP/++xk6dCjPPfccubm5XHrppaSlpZGWlsbXX3/tiatmSITExET++c9/AlYHpunp6aSkpHDZZZdRUlLiKTcjI4OXXnrpsGES5s2b5+lI9G9/+xv9+/cnOTmZG264wZN0QkJCeOihhxgwYMBhvT18+umnPPvss7z22msMHz78sCEbAObMmcPs2bM9y7/nnnvo378/3bt3Z8WKFYCV2O666y7Pejz//PPMnTuXPXv2MHz4cM+ge3FxcRw8eBCAp59+mvj4eOLj43n22WcByM7OplevXlx//fX06dOH0aNHU15eXscnfYoa2l13Yz2sRZ8kh0PkxhtFQGT0aJGsrBPPX14u8u67Iuefb70nIEBk2jSRr74SOXBAZOtWkaeeEomOFvHzE3nllZOPTSkbHdYt/+23iwwd2riP22+vM4bmNAzEjTfe6JmeOnWqrFixQkREduzYIT179hQRkbvvvltur7Ve+fn5kpubK4MHD5aSkhIREfnjH/8ov//97z3lrlq1Sg4cOOAZIkJEZOzYsbJixQrZsGGDjBs3TqqqqkRE5MYbb5S33npLRKxhHBYsWHDM7fbwww971vvIIRyeeuopefjhhz3Lv/POO0VEZPHixTJixAgREfnLX/4ikyZNEofDcdj26ty5s+Tm5nrKqpmu2e4lJSVSXFwsvXv3ljVr1sj27dvF29tbfvjhBxERueyyy+Sdd945ZsyncxiI5sfHx6qai4+3zop69oQrroDLLoPkZAgNhYMHrbObL76Af/wDCgqsG2fnzIFp0yAy8pfyoqOte5euuQauvhpuuAEOHIAHH7RtFZVqqZrTMBC1h3348ssvD7u2U1RURHFxMV9++eVhvZSHh4fzySefsGHDBgYNGgRAVVUV6enph5UdHR1N165d+e677+jWrRubN29m0KBBvPjii6xevdrTH155eTkxMTGAdf3s0ksvrdd2PJFJkyYB1lAQ2dnZnvWbOXMmPj7Wbj0iIuKEZaxcuZKJEyd6egWfNGkSK1asYPz48XTp0sUzzEXtZTSVlpmIwOqB4aabYPx4ePJJmDcP/va3o+cLDraq4K67Di64ALxOUBsZHQ2LFlnz/u531n1M1afaSrU41VUtp1tzGgai9tALbrebb7/9lsDAwMPmEZGjli8ijBo1ivfee++E5V9xxRX84x//oGfPnkycOBFjDCLCNddcwxNPPHHU/AEBAXh7e9cZ94mGgoBfhmqoPUzDsdbjRKQBQ0E0ddVcy7pGdCyxsVaDg/37YcUKePllePppeOstq3FDfr51H9LIkSdOQjV8feHNN2HCBPi//4PFi5t+HZQ6A9k9DMSRRo8ezQsvvOCZzszMPObzhw4dYuDAgXz99ddsq753sayszDNwXG2TJk1i4cKFvPfee56zrxEjRvDBBx9w4MABAPLz8z3rWl9t2rThwIED5OXlUVlZySeffFKv9Xv55Zc9iSk/Px84/nAYQ4YMYeHChZSVlVFaWsqHH354zJ7PT4cWkYiysrJ48sknTzxTYCCcf75VrXbHHVYVW3q6dS9SQ3l7w7vvQlKS1avDzp0nF7hSTcnthu+/hz/+ESZPtoZV2b3bqpJetw42b4asLNi1C3JzrZaiVVVwGnvcP3IYiIyMDPr168e77757zGEgLr300jqHgbjkkkvqNQzEkebOnUtGRgaJiYn07t2bl19+GYAHH3yQQ4cOER8fT1JSEkuXLiU6Opp58+YxdepUEhMTGThwIJuO0bI2PDzcM4pp//79AWtU08cee4zRo0eTmJjIqFGj2Lt3b4Ni9fX19TRqGDdunGdbncj06dPp1KkTiYmJJCUl8fe//x2AGTNmcOGFF3oaK9RISUlh2rRp9O/fnwEDBjB9+nT69u3boDgbi23DQBhjpL7LXrFiBS+88AILFixo4qiOsHUrpKRYCWnZMuvalFJ227XLqop+803Yvt16rnt36NqVjfffT6/YWCtJORzW48jk4+trVVkHB1vXU4ODraruM5EIuFzW+nl5nbnraZPGGgaiRexZS0tL2bdvX5OU7XA4cLvdxx4Uqls3eOUVuPJK+P3v4dFHmyQGpeokYh0MPf88fPSRlWhGjLC+l2PGQPXFcDZuhC5djn5vVRVUVkJ5OZSVQUmJ1YAHrAOsVq2gdWvrUZ8q7ObK4bDWq6gISkut6ZokbIw17llQkJWAW7U6uRqTlsjptD77ykprm9QkZ29vaxsEBFi1SjYl6haRiEpKSposEb333nvccsstFB2vU9Vf/9pqeff449ZIsjbVoR6mrMyqLszPt744rVpBp07Wj6upOZ3WEXlOjrX8kBCr1aKOrts0SkutauLnn7eq2yIj4e67YcaMoxPO8RgD/v7Wo1WrX553OKzqusJC65Gfb+2YIiIgKso6U2oJRKzkc/CgtR5g7VyDg60drI+PlbidTqiosObJy7PmCw21tml4uLXuZ4qqKisZFxdbBx2VlYe/bszRVbReXtbvOSzM2h6nMUnXKxEZY8YCzwHewGsi8scjXjfVr/8KKAOmicjRA7qfpKZMRIWFhce8kPfhhx8SHBzM6NGjrZ3AihXwm99Yw1SczruSRWDDBqu/vW++sXolz8k59rwdO8KAAVaHsenpkJpqHeWcioMH4fPP4euvrU5o1661fsxHSk62bjK+6qqz5yizKe3bZzW6efVVayfbty+88QZMmVLnZ1rv1lO+vlbSiYiwvmfFxdbnnZdnXVMKDYV27ay/zbFKy+Gw4s3NtXa8fn7Qtq21Pic6uhexzg4KCqx1zc62Duyio6FNm5b7/S0vh0OHrEdNKzcfH+vzi4qyzgRrErO39y/VllVV1vylpVby2rXLegQHW9vkOEm6MS/r1JmIjDHewIvAKCAHWGWMWSQitTtbuhDoVv0YALxU/bdRlJSUUFRURFlZmedeg8aSV3NkdITJkydz2223WYkoNBT+/ncYNAhuvtk6Qm1qmzbB66/Dhx9aF5wBOne2zsji463/o6J+ORrMzraSxP/+BzUtinx9rR1Yejqcd571t2PHEy/X5YKMDPjsM+uxapW1jFatrOtlN90EvXtb5UREcOeMGcTl5XGbywXTp1tN6efMsXq/aI47r2Zu908/UfXQQ3T597+tHcSll8Jtt1nfvXpsz4CAAPLy8oiMjGxQU17PmXWrVtaZw8GDVkvULVusHVKHDs0nIZWVWff55eVZ383QUKtGICysfvEZY+2Ug4KsRFtaapW3f7/1NzzcSmiNvK9pEm63lXhyc60zH7DOamJjrc/yRAnZGCsp+fhY61pzb2VFhVVmTZLetct6LSbGSmRYSSgvL4+A6ulTVWdjBWNMOjBbRMZUT99XHcgTteZ5BVgmIu9VT28GhonIcZuKeHl5ibuO7nnKysrYvHkzKSkptG7dmltuuYWbb76ZmJgYvI5Tj11zp25RURGhoaEYY46a1+VyUVZWRk5ODsnJyVRVVfHYY49x9dVXEx4eTkFBAX369GHChAm8/fbbv7zxscfgd7/D/fbbeF111QljP5GCggKCg4Px9fU97PnCggLyFywg7r33MP/9r5VIRoywmpJfdNFxk4iIUFVVhYhYX4wDB6yEVHMG9f33vxwhxcZaCalXL4iLQ8LCcFVU4L1rF+Z//7O6UCooQIyhIjER3/Hj8bn4Ytx9+7LvwAHat29/2HYcNGgQZWVlfPfttwQtXWrdGLx5sxXv3LnQtetJb6caBw4cYM6cOTz++OOem/Xq5HQiW7fCunWYTZusH2rNNZGQEOuo+dxzrYv8ffue+pnjqRKBf/yDsunT8SspwfvqqzEPPkh5bCz+/v7H/b6vWrWKZcuWcdddd+F0Onn44YeZNm3aUf2YnQynw4G7uBi/sjLrACUgwNrZN9LOp6HcFRV4FRVZ32VjfmlscYwzmGPtW463DT2cTuuMoKTE+jwCA62duU3re0IiUFaGFBRgnM5fznyCgxu3irGiwtoepaXWdGCgtZzAQAICAoiNjT1qP3YyjRXqk4gmA2NFZHr19FXAABG5pdY8nwB/FJGV1dNLgHtEJOME5crgwYMxxuB2uwkMDCQsLIy8vDx2795NTk7OYR31HUtwcDDBwcH4+/vj5+eHt7c33t7eBAUFUVFRQXFxMd7e3rjdbtxuN5WVlZSVlXn6jKqvtm3bEhYWRoCvL69lZdGrqoqrEhLI8fVFRCgoKCAuLo7NmzfTu3dvgoODcTgclJeXU1VVddSR6ZYtWwgMDKRPnz4UFhZijKFzcTH/b/t2Eg8dYp+/Pwvbt+ffbdtyyNcXl8uFl5fXUeXUnvbz86O0tJTMzEx69epFVVUV/v7+1jZxu4krLqZPQQFpLhd9iopoW1l5VNv9XT4+rA4LY014OBnh4QR27Mj+/fvx9vamrKyM3NxcYmJiuPDCC/nuu+/44osvDnt/u3btiAoL45qiImbu24cP8I+uXVnUsycVWAnT39+fgIAA9u3b59lZ1Pyt2SYbNmygrKyMgIAASkpK8PX1JSkpia+//ppQPz8inU5C3W6CnE6C3W4CXS6CnE66RkURVVJC5M8/072sjIBaO6MSb2+KfXwQINDtJtTh8FQHOI1hS0gIa1u14n8REeTExRHdoYNnm+7evZvS0lK8vLw81RE1Ny6C1dQ2JibmsKPDyMhIioqK2L59Oy6Xy5MYoqKiMMZQWlpKSEgIDoeDsMJCZv70EwMPHmRTSAj/njSJd9atIyQkhLKyMlwuF1FRUfj7+x92LdPlchEZGUlcXByZmZk4nU7S09P59ttv8fHxOez7UY/f+WHrBtaNjNHR0eTt3s24vXv5zc6dRDocbOnalX8mJbG8vLzO32jN97b235rtWFFRgTEGHx8fz+9ERAgJCaFt27aUl5ezZ88e4vPzuWbXLtIKCijw9WVBhw4sateO0loJqKY6smYdgoKC8PPzw+1243K5MMZQXFyMMcYzX+2/tYU4HEzYs4dLd+8mwuFgS2goC9u25as2bag4zk5eRIiNjeWCCy5g5cqV5OTk0KNHDxwOB1u2bPEcINcsy8vLi+DgYM/+T0Tw8fHB6XRSWlrq2VYulwsfHx+ioqLo3bs3GzZsoO2GDczcvp2eJSVkBQfz7yFD+KC0FDcctS71deR2qP19McYQUVXFuD17GL93L1FVVewLDeXL7t1ZHBPDvvJyz/wiwvLlyxuciOrTJ9xlWNeFaqavAp4/Yp7FwPm1ppcAqccoawaQAWR06tRJqqqqpLKyUiorKyUvL0+2bNkie/bskbKyMnG73eJwOKS4uFji4uLk8ssvl0cffVT27dsneXl5UlpaKm63+5j9H9XF5XJJRUWFZGVlSdeuXWXYsGHy4osvyrZt22Tz5s2yaNEiuemmm+See+6RnJwccblcv7x5+3aRVq1EzjtPXJWVUlJSIps2bZLY2Fjp3r27bNmyRXbt2iW5ublSXl5+VIxut1tcLpe89tprMnfuXJGKCpGHHxbx9RWJihJ5/nnruSPeUx8Oh0Pef/992b59u+zateuo1ysrK2tPiGRlyaGlS6UPyPOPPVav8rOysmTcuHHyz3/+U8rKyqSiokKKioqkvLz8sPWr/PlncUyeLALiiIuTgnfflYO5ufLzzz/L1q1bxeFwiMvlOmzbFhcXy+jRo2X9+vUiIlJUWCilGRkizz0nMm6cSESE1VfgCR6lxlh9Ct5+u8i8eSIZGSKlpUevTFWV7FiyRDY/9ZTIvfeKDB5sfQYgrqAgKRo+XLLvvVd+XrrU09/Y8VRWVsquXbskKyvL8/jss8/ks88+k4qKCnE4HJ7PcNeuXbJv3z4pKiqSfbt3S+njj4s7JETcQUHinjPH6kfxOMsoKCg46e98oygtFeef/iTlISHW53rRRSIrVoicICa32y1Op1McDodUVlZKRUWFlJWVeX7jLpdLKisrPevldrulsLBQNqxfLwfee0/cQ4dan22bNiJz5ogc8VnUlNEUSg4elOtBDnXsaMXQqpXIzTeL/PjjUfO6XC7ZsWOHPPHEE7Jy5UqpqKiQzz//XD7++GNP32+1OZ1OKSwslIKCAikuLpbi4mLP/26327P/c7vdUlVVJXv27JF/P/GElAwebMXSsaPIW2+JOJ1Nsu7HVVkplW+9JeWpqSIg7sBAcV1xhcjHH4tU963HSfQ1V59ElA58Xmv6PuC+I+Z5BZhaa3oz0O5E5aampjbp9mpSf/+7temqO0EUEUlLS5NRo0Y1rJzvvhPp08cq68orRWp1THg6de3aVZYvX940hf/nPyI9eljrGB8v8sILIjt2HH/+nTutDmqnTRPp0OGXJHPuuSLTp4s89pjIq6+KvP++yL//LfLNNyI//SSSnS333XCDjLrggpOPtbhY5KOPRGbOFOnc+Zdl9+olcued1roccZBw0lavFqn+McvYsdYBTktRVGR9Dq1bW/F36SLyu99Zn8OpJkqHw/oM0tOtstu3F3nmmWMfTJwGBw8eFJfTKbJypchvfiPi72/Fdd55Im+/LVJW1vRBfPedyGWXWcuNiLAScvWBn61Wr7Z+KzUHiJGRIjNnNlki8gF+BroAfsBaoM8R81wEfAYYYCDwfV3ltuhEJGJ9Kb29rR2hiIwcOVLOO++8+r23tNTasXl5icTGiixe3ISBNgOVldbZSVLSLzv3Tp1ERowQmTTJegwaJNKu3S+vR0RYP76//rXeO2mn0ykVjZUo3G6RjRtFnn5aZNQoq1d2EAkKss7OXnlFZPfuhpe7ebPIlCkixlhH+fPnn/rO2y5FRdZR+ahR1ncZRLp2Ffm//xNZutRzhFwnl0vk++9F7r//l4OPzp1FXnqpeexwazt4UOTPfxbp1u2X7+kdd1ifa2PKy7MOuAYOtJYTFmZtn0OHGnc5jaGy0jojmjpVJDDwpBJRvXpWMMb8CngWq/n2GyLyB2PMzOqqvZerm2+/AIzFar59rZzg+hBAv379JCPjhLM0b4WFVpPlqipYuZLc6vr+2hfzj2nZMqt1WVYWzJwJf/rT4fd2nMlErIYMixfDDz/Atm2/XBiOibFaPvXrZ7XwS05uXvd1lJZaDTk++ww+/dRqTQSQlmZ1qjt+vDVo47Hq6MvK4MsvrQEYFy+2Ln7ffrt1P1AzG6DspO3dCx9/bN1su2SJdd9KUJC1fQYOtG4Oj421Gha4XFaDmqwsWL3aGn15/37rPpaxY63xwS66yGqs01yJWN+Hl1+2WrY6nVbDohtugAsvtBrENITTaf0m/vtf67uyZIn1XPfuVsfL11xzeu4TPEmczqAAACAASURBVFUlJZjQ0MZvrNBUWnwiAmuQvaFDreaeS5ZAXNzx592+3Rqy4v33reEoXnsNhg07XZGqxiQC69dbPbUvWmS1UAQrkdY0kffzs5rAbtxotVysrLTuUfntb60k1KaNvevQlEpKrPve/vtfq+XmDz9YO9UjGWPtaPv1s3beY8cePjxLS7Fvn3WrxV//at2P5OdnfQ8SE63feliYlZSdTuvep/Jyq3n8gQNW34CbNlndiTkcVnndu1stZadMsVp0Nocm8w3QJK3mmsoZkYjA2gmNHWsdvb/wAlx++S9dpIhYP8IXXrDuPfLxsY6CZ81qGfcoqPrZt8860/n4YytB7dplHfWHhlpnAoMGWd3wXHBB8z7KbyoOh3UT9p491v9eXlYi7tCh4WcOzZnLZSXfzz6D5cutG9FP1EI3LMy6X6lHD+t2ir59YcgQ696mFkwTkV22brUG5vvhB+tLlJpqJaHMTOuIJyjIOrV+4AHtCkeps0XNzaY1/d75+lpnS/7+1pnfsfq3PAOcsZ2eNnvdulk9EHzwAfzzn1ZiAutax6hRVhf94eH2xqiUOr28vKyE0xKrG08zTUSNxdvbOiuqNTSxUkqpurXg/t6VUkqdCTQRKaWUspUmIqWUUrbSRKSUUspWmoiUUkrZShORUkopW2kiUkopZSvbelYwxhRjDReh6hYFHLQ7iBZCt1X96baqP91W9ddDRBrUQ6udN7Rubmg3EGcrY0yGbqv60W1Vf7qt6k+3Vf0ZYxrcd5tWzSmllLKVJiKllFK2sjMR/dXGZbc0uq3qT7dV/em2qj/dVvXX4G1lW2MFpZRSCrRqTimllM00ESmllLKVLYnIGDPWGLPZGLPNGHOvHTG0BMaYjsaYpcaYjcaY9caY2+2OqTkzxngbY34wxnxidyzNnTGmtTHmA2PMpurvV7rdMTVHxpg7qn9764wx7xljAuyOqTkxxrxhjDlgjFlX67kIY8x/jDFbq//WOSroaU9Exhhv4EXgQqA3MNUY0/t0x9FCOIH/JyK9gIHAzbqtTuh2YKPdQbQQzwH/FpGeQBK63Y5ijOkA3Ab0E5F4wBuYYm9Uzc48YOwRz90LLBGRbsCS6ukTsuOMqD+wTUR+FpEqYD4wwYY4mj0R2Ssia6r/L8baWXSwN6rmyRgTC1wEvGZ3LM2dMaYVMAR4HUBEqkSkwN6omi0fINAY4wMEAXtsjqdZEZHlQP4RT08A3qr+/y3gkrrKsSMRdQB21ZrOQXeudTLGxAF9gf/ZG0mz9SxwN+C2O5AWoCuQC7xZXZX5mjEm2O6gmhsR2Q3MAXYCe4FCEfnC3qhahDYishesg2kgpq432JGIzDGe0zbkJ2CMCQH+CfyfiBTZHU9zY4wZBxwQkdV2x9JC+AApwEsi0hcopR7VJ2eb6msbE4AuQHsg2BjzG3ujOjPZkYhygI61pmPR093jMsb4YiWhd0XkX3bH00wNAsYbY7KxqnovMMb8zd6QmrUcIEdEas6uP8BKTOpwI4HtIpIrIg7gX8B5NsfUEuw3xrQDqP57oK432JGIVgHdjDFdjDF+WBf/FtkQR7NnjDFY9fgbReRpu+NprkTkPhGJFZE4rO/TVyKiR67HISL7gF3GmB7VT40ANtgYUnO1ExhojAmq/i2OQBt11Mci4Jrq/68BPqrrDae9920RcRpjbgE+x2qF8oaIrD/dcbQQg4CrgJ+MMZnVz90vIp/aGJM6M9wKvFt9MPgzcK3N8TQ7IvI/Y8wHwBqsFqw/oF39HMYY8x4wDIgyxuQADwN/BP5hjLkOK5lfVmc52sWPUkopO2nPCkoppWyliUgppZStNBEppZSylSYipZRSttJEpJRSylaaiJRSStlKE5FSSilbaSJSSillK01ESimlbKWJSCmllK00ESmllLKVJiKllFK2qlciMsaMNcZsNsZsM8YcNYCWMWaYMabQGJNZ/Xio8UNVSil1JqpzGAhjjDfwIjAKa0CtVcaYRSJy5PglK0RkXBPEqJRS6gxWnzOi/sA2EflZRKqwRsCc0LRhKaWUOlvUZ2C8DsCuWtM5wIBjzJdujFmLNez3Xcca7M4YMwOYAeDl5ZXat2/fhkeslFKq2Vq9evVBEYluyHvqk4jMMZ47cjS9NUBnESkxxvwKWAh0O+pNIn+leoTD0NBQycjIaEisSimlmjljzI6Gvqc+VXM5QMda07FYZz0eIlIkIiXV/38K+BpjohoajFJKqbNPfRLRKqCbMaZL9fj2U4BFtWcwxrQ1xpjq//tXl5vX2MEqpZQ689RZNSciTmPMLcDngDfwhoisN8bMrH79ZWAycKMxxgmUA1NE5MjqO6WUUuooxq58ERoaKsXFxbYsWynVMA6Hg5ycHCoqKuwORTUTAQEBxMbG4uvre9jzxpjVItKvIWXVp7GCUuosl5OTQ2hoKHFxcVTXwquzmIiQl5dHTk4OXbp0OeXytIsfpVSdKioqiIyM1CSkADDGEBkZ2WhnyJqIlFL1oklI1daY3wdNREoppWyliUgp1ewNGzaMzz///LDnnn32WW666aYGlfP4448fNn3eeeedVDyzZ89mzpw5J/XeliQzM5NPP/20yZejiUgp1exNnTqV+fPnH/bc/PnzmTp1ar3eLyK43e6jEtE333zTaDGeiTQRKaVUtcmTJ/PJJ59QWVkJQHZ2Nnv27OH8888H4KmnniItLY3ExEQefvhhzzy9evXipptuIiUlheuuu47y8nKSk5O58sorAQgJCfEs48knnyQhIYGkpCTuvdca7ebVV18lLS2NpKQkLr30UsrKyk4Y5/79+5k4cSJJSUkkJSV5Et3TTz9NfHw88fHxPPvss574evbsyfTp04mPj+fKK6/kyy+/ZNCgQXTr1o3vv/8esM6+rrrqKi644AK6devGq6++CljJddasWcTHx5OQkMCCBQsAWLZsGcOGDWPy5Mn07NmTK6+8kprbdFavXs3QoUNJTU1lzJgx7N27F7DOOO+55x769+9P9+7dWbFiBVVVVTz00EMsWLCA5ORkT/lNQkRseYSEhIhSqmXYsGGD3SHIr371K1m4cKGIiDzxxBNy1113iYjI559/Ltdff7243W5xuVxy0UUXyX//+1/Zvn27GGPk22+/9ZQRHBx8WJk1059++qmkp6dLaWmpiIjk5eWJiMjBgwc98z7wwAMyd+5cERF5+OGH5amnnjoqxssvv1yeeeYZERFxOp1SUFAgGRkZEh8fLyUlJVJcXCy9e/eWNWvWyPbt28Xb21t+/PFHcblckpKSItdee6243W5ZuHChTJgwwbOsxMREKSsrk9zcXImNjZXdu3fLBx98ICNHjhSn0yn79u2Tjh07yp49e2Tp0qXSqlUr2bVrl7hcLhk4cKCsWLFCqqqqJD09XQ4cOCAiIvPnz5drr71WRESGDh0qd955p4iILF68WEaMGCEiIm+++abcfPPNx/1MjvW9ADKkgflA7yNSSjXYlClT2LdvX6OV17Zt26Oq3o5UUz03YcIE5s+fzxtvvAHAF198wRdffEFNb/4lJSVs3bqVTp060blzZwYOHFjn8r/88kuuvfZagoKCAIiIiABg3bp1PPjggxQUFFBSUsKYMWNOWM5XX33F22+/DYC3tzdhYWGsXLmSiRMnEhwcDMCkSZNYsWIF48ePp0uXLiQkJADQp08fRowYgTGGhIQEsrOzPeVOmDCBwMBAAgMDGT58ON9//z0rV65k6tSpeHt706ZNG4YOHcqqVato1aoV/fv3JzY2FoDk5GSys7Np3bo169atY9SoUQC4XC7atWvnWcakSZMASE1NPWzZp4MmIqVUg9WVNJrCJZdcwp133smaNWsoLy8nJSUFsGp17rvvPm644YbD5s/Ozvbs/OsiIsdsjjxt2jQWLlxIUlIS8+bNY9myZQ2OW07Qe42/v7/nfy8vL8+0l5cXTqfT89qRsRlj6l2ut7c3TqcTEaFPnz58++23J3xPzfynk14jUkq1CCEhIQwbNozf/va3hzVSGDNmDG+88QYlJSUA7N69mwMHDhyzDF9fXxwOx1HPjx49mjfeeMNzDSg/Px+A4uJi2rVrh8Ph4N13360zxhEjRvDSSy8B1hlHUVERQ4YMYeHChZSVlVFaWsqHH37I4MGDG7TuH330ERUVFeTl5bFs2TLS0tIYMmQICxYswOVykZuby/Lly+nfv/9xy+jRowe5ubmeRORwOFi//qhh4w4TGhrK6eiKTRORUqrFmDp1KmvXrmXKlCme50aPHs2vf/1r0tPTSUhIYPLkycfdec6YMYPExERPY4UaY8eOZfz48fTr14/k5GRP0+xHH32UAQMGMGrUKHr27FlnfM899xxLly4lISGB1NRU1q9fT0pKCtOmTaN///4MGDCA6dOn09BBQfv3789FF13EwIED+d3vfkf79u2ZOHEiiYmJJCUlccEFF/Dkk0/Stm3b45bh5+fHBx98wD333ENSUhLJycl1thocPnw4GzZsaPLGCtrpqVKqThs3bqRXr152h3FWmj17NiEhIdx11112h3KUY30vTqbTUz0jUkopZSttrKCUUs3Y7Nmz7Q6hyekZkVKqXuyqxlfNU2N+HzQRKaXqFBAQQF5eniYjBfwyHlFAQECjlGdr1dzx2u6r0ys7O5u4uDi7w1DNWGxsLDk5OeTm5todimomakZobQy2tporLCzEy0tPyuxmjGHfvn20adPG7lCUUi1ci2o1V9edwer0crlcdoeglDpL2ZaIysvLPXcgK/vV9GqslFKnm22JyOl0cuutt9q1eHWE8vJyu0NQSp2l9AKNAqCiosLuEJRSZylNRArQMyKllH00ESlAE5FSyj6aiBSgiUgpZR/bE9HmzZvtDkGh14iUUvaxPRG9+eabdoeg0DMipZR9bE9Ep3tsdHW4mpuKNREpdbSoqCitLTgNbE9ECxYsoKioyO4wzlo1N7LqZ6DU0fLy8li7dq3dYZzxbE9EAGFhYWRlZdkdxlmpsrKS6Oho8vLy7A5FqWbF7XYDkJmZaXMkZ756JSJjzFhjzGZjzDZjzL3HeN0YY+ZWv/6jMSaloYGce+65TJgwgZycHO2D7jTKy8vj3HPP1UR0Fvvxxx8pLS21O4xG01j7j4MHD9KxY0c9SD4N6kxExhhv4EXgQqA3MNUY0/uI2S4EulU/ZgAn1YncokWL6NixI15eXhhjGD9+PLfeeis333wzbrfb8xCRBn3ZnE4nBQUFJxNSnQoLC2nXrh1lZWVHxdQS6pY3btzI4MGD2bNnzwnn27lzZ7M5QGgJ27Up5ebmNur3OSkpib/97W+NVp6d7r77bqZOndooZX377bdceuml2rL3NKjPeET9gW0i8jOAMWY+MAHYUGueCcDbYu2pvjPGtDbGtBORvacS3Mcff+z5/9xzz+XOO++s8z3nnXce33zzzVHPn3POOSxevJi1a9eyfft2Kisr8fb2JiEhgZiYGJYsWUJVVRXGGEJDQxkzZgzLly9n165dGGNwu90YY/Dy8sLpdOJ2uykvL+e9994jPz+f4OBgrr32WgYNGoTL5eKnn37ihRde4PXXX6d6u+Ht7Y2vry979uxh3bp19OvXj9DQUHx9fcnLy+Pzzz/nz3/+M6WlpRQVFZGVlcWbb77J73//e/z9/Tl06BCFhYUEBQXx6quvMmrUKNasWUN0dDRBQUGEhoayYsUK2rZtizGG9u3bA/DRRx8RFRXFiBEjiImJwd/f39P7+fz587njjjt48sknMcbw7LPPkpyczMqVKwkMDMTpdLJu3TreeecdHn/8cdq2bUunTp0wxpCfn09qaiohISFkZmZSWFiIy+WitLSUtWvX0rlzZ+Li4igqKsLlcnHo0CE6d+5Meno6ERERbNu2jcrKSvLz8+nZsye7d++mdevW7NixA6fTyerVq6msrKR///5ccMEF/Pjjj+Tk5HDLLbfw2GOP0a1bN8/BCUBOTg4bN24kMTGR7du3k5aWRlhYGKGhoTz33HP4+fkxfPhwli1bRkxMDDExMZSUlBAVFUVISAhBQUH4+PiwZMkSVq9ejb+/PxdeeCFdu3YlKCiI/Px8/Pz8CAkJweFw4OXlRVpaGuHh4VRVVbF+/Xo2btxIeHg45eXl+Pn54Xa78fLyIjIykoKCAg4dOsShQ4cIDw9HRHjyySeZNGkSQ4YM8WwvLy8v9uzZw9q1axkwYAAZGRkUFxfjcrkoLCxk7ty5hIeHc91111FVVUVQUBBBQUG0bduW1NRUtm3bRklJCU6nE6fTyc6dO8nLy6NTp04AtGvXjvz8fPLz88nKymLYsGG8+uqr9OrVC4fDwf79+4mKiuLAgQMcOnSI2NhYoqOjWb9+PcHBwXz88cfExMRw3nnnsXr1an744QduvPFGvL29cblchISE4Ha72b59O3FxcQQGBiIi+Pj44HK58PHxISgoiIMHD7JlyxZ++uknBg8eTEhICMHBwVRWVrJlyxbatWtHZGQkLpeLtLQ0CgsLycvLIzIyEofDQWlpKaWlpZSVlXn+X7p0KW63m0cffZTo6GgcDgcxMTEAeHt7ExAQQEREBG63m4qKCqqqqmjdujVFRUVs3bqV//3vf5x33nkUFRXx97//nUWLFnHvvfdy0003ERgYSKdOnYiJicHb2xs/Pz+ioqLo2bMnGzduJCAggNDQUHbv3k1hYSHvvPMO3t7eXHjhhfj7+x/28PX1xe12H/b7MMbQoUMH2rRpw08//YQxht27d1NRUUFAQIBnu5aWluLn50dQUBBt2rQhIiLCM19UVBSHDh1i7969eHl54eXl5YnV7XZTUFCAw+GgqqoKb29vKisrCQoKwsvLCz8/PyIjI6msrMTtdnsOrsvKyjzf59atW1NaWkphYSEi4vlsQ0JCaNWqFb6+vie1r69zPCJjzGRgrIhMr56+ChggIrfUmucT4I8isrJ6eglwj4hkHFHWDKwzJoDUhgSakpLCmjVrGvKWY5o1axbbtm0jLy+PiooK1q1bR1lZGR07dsTPz4+AgACcTidVVVVs374dgLZt2xISEoKfn5/nA6s5a2vXrh133nknq1evJioqio0bN3q+KPHx8fz8889UVFR4zuJ8fX0555xzWLt2LZ06deLAgQNUVFSQm5tL3759yczMJDw8nPDwcEJDQ+nRowfz5s2jXbt2BAUF8cwzzxy1TsHBwZxzzjl0796dwMBAXC4XRUVFni9f586dGTNmDKWlpWzduhWHw4Gvry8BAQHEx8czatQoPvzwQx544AE2btzI0KFDSUpKYsuWLYSHh+Nyubjxxht56qmn6N69O6+88gqTJ08mNzeXzMxM4uLi6NWrF2FhYfj6+uLr68v555/PN998Q35+vucs1u12s2rVKjIzMxk5ciT79u1j1KhR+Pj4sGPHDjp37sy+ffsIDw8nJCSEtm3b4u3tzV//+lc2b97MDTfcQIcOHRgxYgRff/01BQUFeHt7ew4OAPr06cPOnTvx9/f37GyLi4sZO3Yshw4dIj8/n8TERLZt24avry+VlZX4+fmRl5eH0+nkwIEDjBkzhtjYWJxOJ5s3byYrK4uSkhLat2+P2+2mpKQEsFoafvbZZ9x///3ccccdXHDBBfTt25eysjL8/PwA6zpDTQKJiIggPDyc1q1bk5+fj5eXF3369GHjxo1s2rSJt99+mxtuuIHbb7+dmJgYAgICqKio4De/+Q2RkZF4eXnRqlUrEhMT6dSpE19//bVnZ7N//36++uorDhw4QM+ePenRowc+Pj74+PjQqlUrz8B2VVVV7Ny5k9jYWCIiIoiIiGDMmDHs2LGDefPm4eXlRfv27Tl48CAdOnQgJCSE/Px8cnJySExMxOVy0bt3b/bv38+uXbvo2rUrgYGBrFmzxjPQZUVFBeXl5fTs2ZOtW7fidDo9v6uazyo3N5cuXbrQsWNHunbtyvr16yktLaWkpAR/f3+6du3Knj17PDu8zMxMT7wHDhwgMDCQ4OBggoODCQoK8vx//vnnExkZyZIlSygsLCQgIIDc3FyMMTgcDiorK9m7dy9hYWGeA62qqirCwsIIDw8nPT2dNWvWEBgYyODBg2nfvj2VlZWsXbsWYwzbt28nPz8fl8uFw+Fg3759rF271rNtiouLPckhPT0dYwybNm2isrLysEdVVRVeXl5s376d4OBgAgICEBG+++47ioqKSEtLIyAggHbt2nmSc1FREb6+vgQHB1NRUUFpaSn79+8nOzub+Ph4AgMD2bVrF23atKF9+/aICC6XC7fbTWVlJSJCdHS05zda87mUlJQgIlRVVZGbm4u/vz8+Pj4EBwdjjCE4OBg/Pz8cDgcFBQUEBwfTunVrwKqdMMZQXFxMcXExDoeDxx57rMHjEXl2kMd7AJcBr9Wavgp4/oh5FgPn15peAqTWUa7U59GjRw8BZNCgQZKXlyf79++XvLw8KSsrk7KyMqmsrBSn0ykul0vcbrfs3r1bHnzwQXE4HOJ0OqW4uFg+/PBDadOmjdx5551ypPfff1+szXA4Hx8fueuuu4563m5ut1sGDRokvXv3lp9++kncbnejll9YWNigMocMGSIDBw5s0DJ+//vfS2hoqOzdu7de8xcXF8tbb73VoGWcLosXLxYvLy/57LPPTrmsqqoq6dmzp9x2223y8ccfi4hIRUVFvd9fUlIid9xxh5SUlJxyLMoehw4dapTvkp2ADKkjrxz5qE9jhRygY63pWODICwr1madO7du351//+helpaWeADdt2sSuXbv429/+RkREBDExMURERBAYGEhgYCB+fn6eo+Ka6qhHH30UHx8fvL29CQkJYfz48Rw6dIjg4OCjlpmamsr06dOPet7pdBIWFtbQVWhyNdWD0dHRxMfHN/pQ661atWpQmQEBATgcjgYto2PHjp4jx/oICQnh6quvbtAyTpekpCTcbjepqQ06wT8mX19funfvzvfff8+IESMA8Pf3r/f7g4ODefrpp4/5PVctQ+vWrRk7dqzdYZx29blGtAroZozpAuwGpgC/PmKeRcAt1dePBgCF0sDrQ2vWrKFv377HfO1Ux0X38vKidevWhISEHPValy5dePXVV496/r777mPUqFGntNym4na78fGpz0fX9CIjIxvceCAwMBCg0ZOoHWquw0VHRzdKeREREaxZs4aAgIBGKU+plqDOvZmIOI0xtwCfA97AGyKy3hgzs/r1l4FPgV8B24Ay4NqGBPGf//znuEmosURHRx8zER3P448/3oTRnLrmshN//fXXG3xGNHHiRHbu3NlEEZ1eNfXjjSUiIoL8/Pxm8/kqdTrU67BaRD7FSja1n3u51v8C3HyyQYwcOfJk31pv0dHRZ0yVRfv27encubPdYQB4qkgbwt/fn44dO9Y9YwvRkAOcukRERFBWVtZo5SnVEthevyOn6d6Umia6Z4IFCxboEfMZqqapsVJnk2bRxc/p0NCqueZMk9CZq0ePHnaHoNRpd9YkoqFDh9K1a1e7w1DqhLp37253CEqddrZWzb3zzjunbVlXXHHFaVuWUierbdu2uFwuu8NQ6rSy9YyouTRBVqo58fI6ayoqlAJsTkT6g1NKKaWJSCmllK00ESmllLKVJiKllFK2sjUT6P0wSiml6hyPqMkWbEwxoEMf1k8UcNDuIFoI3Vb1p9uq/nRb1V8PEQltyBvsbD+9WRo6eNJZyhiToduqfnRb1Z9uq/rTbVV/xpiMuuc6nF6kUUopZStNREoppWxlZyL6q43Lbml0W9Wfbqv6021Vf7qt6q/B28q2xgpKKaUUaNWcUkopm2kiUkopZStNREoppWyliUgppZStNBEppZSylSYipZRSttJEpJRSylaaiJRSStlKE5FSSilbaSJSSillK01ESimlbKWJSCmllK00ESmllLKVJiKllFK20kSklFLKVpqIlFJK2UoTkVJKKVtpIlJKKWUrTURKKaVspYlIKaWUrTQRKaWUspUmIqWUUrbSRKSUUspW9UpExpixxpjNxphtxph7j/H6MGNMoTEms/rxUOOHqpRS6kzkU9cMxhhv4EVgFJADrDLGLBKRDUfMukJExjVBjEoppc5g9Tkj6g9sE5GfRaQKmA9MaNqwlFJKnS3qPCMCOgC7ak3nAAOOMV+6MWYtsAe4S0TWHzmDMWYGMAMgODg4tWfPng2PWCmlVLO1evXqgyIS3ZD31CcRmWM8J0dMrwE6i0iJMeZXwEKg21FvEvkr8FeAfv36SUZGRkNiVUop1cwZY3Y09D31qZrLATrWmo7FOuvxEJEiESmp/v9TwNcYE9XQYJRSSp196pOIVgHdjDFdjDF+wBRgUe0ZjDFtjTGm+v/+1eXmNXawSimlzjx1Vs2JiNMYcwvwOeANvCEi640xM6tffxmYDNxojHEC5cAUETmy+k4ppZQ6irErX+g1IqXODA6Hg5ycHCoqKuwORZ1GAQEBxMbG4uvre9jzxpjVItKvIWXVp7GCUkodV05ODqGhocTFxVFdQ6/OcCJCXl4eOTk5dOnS5ZTL0y5+lFKnpKKigsjISE1CZxFjDJGRkY12FmxbIiorK7Nr0UqpRqZJ6OzTmJ+5bYkoL08b1SmllNIzIqXUGcDb25vk5GTi4+O5+OKLKSgoaND7Z8+ezZw5cwB46KGH+PLLL084/7Rp0/jggw9OOt6WYtmyZXzzzTdNvhy9RqSUavECAwPJzMxk3bp1RERE8OKLL550WY888ggjR45sxOhaLk1ESil1EtLT09m9ezcAJSUljBgxgpSUFBISEvjoo4888/3hD3+gR48ejBw5ks2bN3uer32288gjj5CWlkZ8fDwzZsygrttdtm3bxsiRI0lKSiIlJYWsrCxEhFmzZhEfH09CQgILFiwArJ380KFDufzyy+nevTv33nsv7777Lv379ychIYGsrCxPPDNnzmTw4MF0796dTz75BLAaiVx77bUk/zf/0wAAG0VJREFUJCTQt29fli5dCsC8efOYNGkSY8eOpVu3btx9992e+L744gvS09NJSUnhsssuo6SkBIC4uDgefvhhz3batGkT2dnZvPzyyzzzzDMkJyezYsWKU/pcTkSbbyulGtWUKVPYt29fo5XXtm1b5s+fX695XS4XS5Ys4brrrgOse10+/PBDWrVqxcGDBxk4cCDjx49nzZo1zJ8/nx9++AGn00lKSgqpqalHlXfLLbfw0EPW8GpXXXUVn3zyCRdffPFxl3/llVdy7733MnHiRCoqKnC73fzrX/8iMzOTtWvXcvDgQdLS0hgyZAgAa9euZePGjURERNC1a1emT5/O999/z3PPPcfzzz/Ps88+C0B2djb//e9/ycrKYvjw4Wzbts1z1vfTTz+xadMmRo8ezZYtWwDIzMzkhx9+wN/fnx49enDrrbcSGBjIY489xpdffklwcDB/+tOfePrppz3rFxUVxZo1a/jLX/7CnDlzeO2115g5cyYhISHcdddd9dr+J0sTkVKqUdU3aTSm8vJykpOTyc7OJjU1lVGjRgHW/S73338/y5cvx8vLi927d7N//35WrFjBxIkTCQoKAuD/t3fv4VHWZ8LHv3cmk0ySyQkCCYE0WEQODSAnF6ogwuqLuxQRtchrPb7brqWoe/By267vpVerddfat6/2wuqurxtFUKwLrIeqWA9QQFtAUFFAkEMgJIZwSCaTmWQyc79/zGQ2hBwmQHiScH+u67kyz8xvfr/7+c3Mcz+nPL+5c+e2We/777/Po48+Sn19PceOHeNb3/pWu4nI5/NRXl7OtddeC0STIMD69etZuHAhLpeL/Px8Lr/8cjZt2kRWVhaTJ09m0KBBAAwbNoyrrroKgDFjxsT3cAC++93vkpSUxPDhw/nmN7/Jzp07Wb9+PXfddRcAI0eOpLi4OJ6IZs2aRXZ2NgCjR4/mwIEDnDhxgi+++IJLL70UgMbGRqZOnRpvY/78+QBMnDiRlStXdqn/z5QlImNMr9d8jqimpoY5c+awZMkS7r77bpYtW8aRI0fYsmULbreboUOHxv/3pbPLj4PBIIsWLWLz5s0UFRXx4IMPdvh/M+0dtuvocF5qamr8cVJSUnw+KSmJpqam+GutYxWRhOt1uVw0NTWhqlx55ZW8+OKLHb6nufy5ZOeIjDF9RnZ2Nk888QSPPfYYoVCImpoaBg4ciNvt5v333+fAgegIBdOnT2fVqlUEAgF8Ph+vvfbaKXU1J528vDzq6uo6vUouKyuLIUOGsHr1agAaGhqor69n+vTprFixgnA4zJEjR1i3bh2XXHJJl5brd7/7HZFIhK+++oq9e/cyYsQIpk+fzrJlywD48ssvKSsrY8SIEe3WMWXKFDZs2MCePXuA6JXLzXtQ7cnMzMTn83Up1tNhicgY06eMHz+ecePG8dJLL3HTTTexefNmJk2axLJly2gejHPChAksWLCAiy++mOuuu45p06adUk9OTg7f//73GTNmDPPmzWPy5Mmdtr106VKeeOIJxo4dy7e//W0qKyu59tprGTt2LOPGjWPmzJk8+uijFBQUdGmZRowYweWXX87VV1/NU089hcfjYdGiRYTDYcaMGcOCBQsoLS09aU+otQEDBlBaWsrChQsZO3YsU6ZMYefOnR22+53vfIdVq1Z1+8UKjt30NDMzU89FpjXGdK8dO3YwatQop8Pos2677TbmzJnD9ddf73Qop2jrsz+dm57aHpExxhhH2cUKxhjTg5WWljodQrezPSJjzBmzcTDPP2fzM7dEZIw5Ix6Ph6NHj1oyOo80j0fU/L9SZ8oOzRljzsiQIUM4dOgQR44ccToUcw41j9B6NlgiMsacEbfbfVZG6TTnr15xaE5VueKKK5wOwxhjTDfoFYkoEAjwwQcfOB2GMcaYbtArElF3DqJ34sQJli9f3m31G2OM6VivSER+v7/b6t64cSM33XRTt9VvjDGmY+d9IrJLTo0xxlm9IhF156E5Y4wxzuoVicj2iIwxpu867xNRY2Njt9VtjDGmc44mokT3RrozEbV32O/5559n9+7d3dZuX/HCCy/w+OOPOx2GMaYXcywRNTQ08PzzzydU1olEdOutt7Jq1apua7evWLFixXlxd2BjTPdxLBGFQiFuu+22hMrW1NR0WxyBQKDd1yKRSLe121eEw2FcLpfTYRhjejHHzxE1NDR0WqY5EYXD4bPefkdX5HVHe31NJBIhKcnxr5ExphdLaA0iIrNFZJeI7BGRH7fxuojIE7HXPxWRCYkGsHDhwk7PFdXU1JCdnU1lZWWi1SbM7/cjIqc8n5SURDAYPOvt9TXnYo9IVXn44Ye7tY2+oKmpiVtuuYW6ujqnQzGmSzpNRCLiApYAVwOjgYUiMrpVsauB4bHpB8BvEw1g1apVJCUl8cgjj7Bnzx4aGxtPSkyhUIjNmzdTU1PDkCFD+OCDD/D5fKeU60x7Zb/66iuGDx9+yiG6AQMGdDnxtV4B9PZLw6urq+PL4PP5uP32208pU19f3+3LuXbtWu6//358Pl+3ttMTrVmzJuFD0x9//DFLly7lo48+6lIbBw8eTOjIxOlqaGjo1vpbCgQC59X3JBKJsGHDBqfDOGPS2UpERKYCD6rq/4jN/wRAVR9pUeZp4ANVfTE2vwuYoaoVHdTb6dpr3rx5vPHGG4RCoUSWpUNz5szhkksuYc+ePZSXl1NRUcEXX3wBwC9/+UuWL1/OyJEjCYfDHDlyhKysLLZu3cqMGTPweDyoanyF2/xYVQmHw6xdu5YDBw6c1F56ejr19fXMnz+flJQUUlNTSU6OjrpRXV3N9OnT+eSTT6itraVfv354PB4KCwv56quvcLlcJCUlkZKSwqeffsqUKVPIzc3l+PHjuFwuCgsL+fWvf01FRQUlJSWkpqZSX1/PoUOHePrppzlx4gS1tbVUVFQQCAQYMGAAAEePHo0ncFUlEonQ2NjIxo0bGTlyJAUFBQwaNIhIJMLBgwdZu3YtxcXF5OTksHfvXioqKrjqqqvwer1EIhHKysooKCggNzeXffv24fF4cLvdZGdnk5GRQSgUIhwO09jYGF+eESNGUFhYiNvt5sMPPyQ5OZmvv/6ad999l+LiYi644AJcLhdut5uSkhL8fj9PP/00+fn5lJWVcfnllzN16lQ+/vhjQqEQkUiEhoYG/H4/999/Py6Xiw0bNjB16lTefvttwuEwJSUlvPXWW4TDYXJycjh8+DDp6elMmjSJ4uJi1q1bRzAYREQQEQYOHEhycjJ+v5/a2lqCwWC8v9LS0lBVmpqaiEQibNy4kcWLF7Np0yZcLhdNTU243W4uvPBCLrroItavX09DQwNut5v09HTS09PJyclh//79uFwuPB7PSRtVKSkpZGRksGfPHmpqali/fj25ublcdNFF8UOgVVVV1NXVMXz4cJKSkuJx79ixg3nz5rF8+XLGjRtHRkYGKSkpZGVlkZ6ezrFjx2hqaqJfv3643e54P6xZs4aMjAyuuOIKAoEAVVVVfP3111x44YWEQiFcLhdpaWlkZ2cTCAS48cYbKS8vZ/PmzaSlpZGVlUVZWVn8e9v8+c2cOZPXX3+d119/nWAwyLRp0xg6dCglJSXs2rWLiooKgsEgaWlpJCUlEYlE8Hg88en48eOEQiEyMzPZt28fXq+XcDh8yqHgYDBIZWUlJ06cIBKJ4Pf7mTRpEh6PB5/PR0NDAy6Xi6ysLDIyMkhOTo4nx/r6eurq6rjgggsoLCzE5XJx+PBhQqEQx48fp6qqiqqqKjIyMsjMzEREKCwsJDU1FY/HQ2pqKqpKIBCgoaGB1NRUAHJycrjuuusQEVavXk11dTVut5vGxkZCodBJ65CWU25uLuPHj6esrIxAIEAwGCQzMxOfz0cwGKS2tha/3x//TVVUVKCqeL1eBg8ejN/vj3/HmpqaCIfDRCIRtm3bRiQSYd68eYwZM4YtW7ZQXV3NgQMHyMvLIy0tjdTUVNxuN0A8xubvl8fjob6+nmPHjpGWlkZ6enq8fpfLRSQSiZ/KWLJkyRZVndSV9XMiieh6YLaq/k1s/mbgL1R1cYsyrwP/oqrrY/PvAv+kqptb1fUDontMABPbaq+oqIi5c+cyefJkhg0bRigUYubMmaeUGzZsGMOGDaOwsJC8vDwGDBiA2+2msrKS9957j/nz5yMi8T2qV199lWuuuYa77rqLSCRCIBBARBg0aBATJ06MD+wVDodJT09n1KhRHD16NP5Fbv6xt56qqqqorKxk/vz5AIwePZo77riD/Px8ioqK4luCIoLL5cLlcjF06FDuu+8+Xn75ZZYvX052djZut5sDBw4wYMAAMjMzcblcpKSkkJ+fz8qVK+PL6vP5uPHGGzl69CjPPvssWVlZiAipqal8+OGHPPzww8yePZv+/fuzaNGieDL7wx/+wJw5cygoKCA1NZWkpCSSkpJwu92MGDGCPXv24PP52L17NwMHDmTo0KEUFxezd+9eAoFAPOHt27ePUCiEiFBUVITX60VV2bdvHxA9jFpfX8/evXspKirC7XaTmppKJBKhqamJX/ziF7zxxhsUFRVRWlqKqpKTk0NJSQn79+/H5/OhqoRCIV5++WUGDx7MDTfcQHFxMV9++SVlZWXs2rWLfv36kZOTQzgcxufzsWHDBp588kny8/NZtGgRDzzwACtWrMDj8fDMM8/wwx/+kMGDBxMMBikoKKC8vJxf/epXbN26leeeey6eYADKysoIhUJ4vV4GDBhAWloaED1cGwgE4itbESEpKYkrr7ySFStWkJWVRWpqKqFQiKVLl7Jp0yZ+9rOf4fF4aGpqor6+nvr6eiorKxk2bBiqSlVVFbm5ufGRLv1+P4cPH2b8+PFkZWXxjW98g1AoFO8XgNTUVLKzsykvLz9pJda/f39yc3Opqanh8OHDBINBwuEwtbW1nDhxgvz8fJKTkzlx4gSNjY0MHjyY2tpaJkyYQF1dHTt37iQ5OZn8/Hzy8vIoLy/H7XYTDofx+/3s37+fnJwc7r33XoqLi7nnnnuoqanh+PHjjBo1inA4HJ/q6uq44447KC0tZcaMGXi9XrZv387u3bvZunUrkydPprCwkLS0NAKBQPxcY0NDA8FgkEAgQE5ODm63m9raWoqLiwkGg/GVXksej4f+/fvHf0eNjY3s3buXxsZGvF4vKSkpqCq1tbXU1dURiUTi30uPx4PX66WsrIzy8nLq6+u58MILcblc9OvXL/57bN4rbWhooLq6Op7Imjdgmlfkzb/3qqoqfv7znyMi/PSnP6WgoIBQKERKSgput7vd9cnhw4fZsGEDEyZMiCfk2tpavF4vGRkZeL1evF5vPBmnpaXRr18/jh07xpEjR/B6vfj9fmpqakhOTo5v/BYXF5OUlMSaNWvYunUrl112GQMHDqSoqIijR4/Gl6d5oz85OTm+cdCcaDMyMsjJySEYDBIMBklOTo5vfDX/JgDGjh3b5UTUbmZusQdwA/BMi/mbgd+0KvMGcFmL+XeBiZ3Uq83T448/rg0NDdoev9+vQ4YM0TfffFMjkUi75TpSUlKiv/nNb07rvYmYPn26Tp8+PeHygUBAy8rKTqutWbNm6cyZM9t8ze/3azgc1lAodFp1d7eKigrdvn27fv7552e97hkzZuill16qqqrbtm3rtHwkEtG6urozbrej725f9N5772lVVZXTYfR4b775pr7zzjtOh3HOAZu1k7zSekpkhNZDQFGL+SHA4dMo06aKigoKCgo6LJOens7BgwcTqa5dkydPpri4+Izq6IjP56OwsDDh8h6Ph6Kios4LtqG0tLTdY+7p6ekAPfZKtoKCgk4/79OVlZXF/v37ARg3blyn5UWEjIyMM243JSXljOvoTWyQysTMnj3b6RB6jUQS0SZguIhcAJQDNwL/s1WZV4HFIvIS8BdAjXZwfqjZnXfe2W0rpdYeeugh+vfv3231BwKB+CGc7na2xonva44dOxZPxMaY3qPTRKSqTSKyGHgbcAHPqurnInJn7PWngN8DfwXsAeqBUy+vasOsWbNON+4u68reyuk4duwYubm53dqG6diCBQsYPbr1BZ3GmJ4ukT0iVPX3RJNNy+eeavFYgR91tfHmKzT6gurqavLy8pwO47y2ePHizgsZY3ocR08k9KVENHfuXKZNm+Z0GMYY0+t0evl2dykpKdFt27bFLy80xhjT+4lIly/fdmyPyOPxWBIyxhjj/E1PjTHGnN8sERljjHGUJSJjjDGOskRkjDHGUZaIjDHGOMoSkTHGGEdZIjLGGOMox/6hVUR8wC5HGu998oBqp4PoJayvEmd9lTjrq8SNUNXMrrzByf8o3dXV/749X4nIZuurxFhfJc76KnHWV4kTkc2dlzqZHZozxhjjKEtExhhjHOVkIvo3B9vubayvEmd9lTjrq8RZXyWuy33l2MUKxhhjDNihOWOMMQ6zRGSMMcZRjiQiEZktIrtEZI+I/NiJGHoDESkSkfdFZIeIfC4i9zgdU08mIi4R2SoirzsdS08nIjki8oqI7Ix9v6Y6HVNPJCJ/H/vtbReRF0XE43RMPYmIPCsiVSKyvcVz/UTkHRHZHfub21k95zwRiYgLWAJcDYwGForI6HMdRy/RBPyjqo4CpgA/sr7q0D3ADqeD6CUeB95S1ZHAOKzfTiEig4G7gUmqWgK4gBudjarHKQVmt3rux8C7qjoceDc23yEn9oguAfao6l5VbQReAq5xII4eT1UrVPXj2GMf0ZXFYGej6plEZAjw18AzTsfS04lIFjAd+H8AqtqoqiecjarHSgbSRCQZSAcOOxxPj6Kq64BjrZ6+Bngu9vg5YF5n9TiRiAYDB1vMH8JWrp0SkaHAeOBPzkbSY/1f4D4g4nQgvcA3gSPAf8QOZT4jIhlOB9XTqGo58BhQBlQANaq6xtmoeoV8Va2A6MY0MLCzNziRiKSN5+wa8g6IiBf4T+DvVLXW6Xh6GhGZA1Sp6hanY+klkoEJwG9VdTzgJ4HDJ+eb2LmNa4ALgEIgQ0S+52xUfZMTiegQUNRifgi2u9suEXETTULLVHWl0/H0UJcCc0VkP9FDvTNF5AVnQ+rRDgGHVLV57/oVoonJnOwvgX2qekRVQ8BK4NsOx9QbfC0igwBif6s6e4MTiWgTMFxELhCRFKIn/151II4eT0SE6HH8Har6f5yOp6dS1Z+o6hBVHUr0+/SeqtqWaztUtRI4KCIjYk/NAr5wMKSeqgyYIiLpsd/iLOyijkS8Ctwae3wr8F+dveGc331bVZtEZDHwNtGrUJ5V1c/PdRy9xKXAzcBnIrIt9txPVfX3DsZk+oa7gGWxjcG9wO0Ox9PjqOqfROQV4GOiV7BuxW71cxIReRGYAeSJyCHgAeBfgJdF5H8RTeY3dFqP3eLHGGOMk+zOCsYYYxxlicgYY4yjLBEZY4xxlCUiY4wxjrJEZIwxxlGWiIwxxjjKEpExxhhHWSIyiEh/EdkWmypFpLzFfIqIbOymdu+OjYWzrIMydW09PoM2z9qyiMjfxvprm4jsFZHbzlbdTjmT/mnv8xGRNBFZGxsC5ozFvpPrYnfENn2A/UOrOYmIPAjUqepj56CtncDVqrqvgzJ1qupt/bgnEJElwGeq+pSITADeUdX+TseViNgta0RVI23Nn2adbX4+IvIjIFlVHz/tgE+t8wGiw8m0uxFjeg/bIzKdEpE6ERkaG83zmdholctE5C9FZENsJMZLWpT/noj8Oban8HRbW8Ii8hTR4QhejY2C+Q+xereLyN8lENMp5UXkPhG5O/b41yLyXuzxrOaboDZvtceWZ4eI/HtsBM41IpIWe+1/x5b1HYmOynlvO2GM4b/vPXaI6C2rWsbYaTyx+dUisiUWxw9iz/2riCxqUeZBEfnHLvRvW3U2L/OTRG9bM63VfFGL/umo/VPq7sRNxO43JiIZIvKGiHwS++wWdLRMInKLiHwaK7+0RZ2rY/WavkBVbbIpPgEPAve2eq4OGEr0fltjiG7AbAGeJTqsxzXA6ljZUcBrgDs2/yRwSztt7QfygInAZ0AG4AU+B8Y3t90yjtjfNssTHcX2d7EyfwT+DLiJ3v/qb1vV0bw8F8fmXwa+B0wCtgFpQCawu3V/tIjnOJAf64OHgBdavd5pPLHX+sX+pgHbgf6x5VnboswXwDcS7d926hxKdLymKS36ID7fqn/abL+9ult/Vi3elwJUtpi/Dvj3FvPZ7S0T8C1gF5DXst3YYxdwxOnfi01nZ7JjrKYr9qnqZwAi8jnR4YBVRD4julKD6B2KJwKbokd7SKPz28BfBqxSVX+s7pXANKI3mexK+d8CE0UkE2ggupU/Kfba3e0sT/PNZLfEliEP+C9VDcTqfq2tAESkiGgSfBsIEU0yP2pVbEuC8dwtItfGHhcBw1X1IxEZKCKFwADguKqWSfSGwYn07yl1ApXAAVX9qEW51vMAqOrWttrvoO6jbfUT0f5sOfrrZ8BjIvKvwOuq+kcRubmdZcoGXlHV6lhM8ZFAVTUsIo0ikqnR0YtNL2aJyHRFQ4vHkRbzEf77uyTAc6r6ky7U29ZgiV0ur6ohiY5JdDuwEfgUuAIYRtu372+5PGGiK8BEYxlLNBHPPimw6PmQ78dm/4roXl+78YjIDKLj3kxV1XoR+QDwxN7/CnA9UEB0nCVIoH87qdPfqnjr+ZZOab+TutsSaPm6qn4pIhOJ9s0jIrKG6J7lKcsUO6zZ0UnsVCDYweuml7BzROZsexe4XkQGAohIPxEp7uQ964B5Eh33JQO4luihrNMpvw64N/b3j8CdwDZVTfSqnPXAd0TEI9GRcf+6nXJjgE9aP6mqS1T14th0OIF4sonubdSLyEiih/OavUR0fKXriSYFSKx/O6qzK9pqv0t1q+pxwCUinli8hUC9qr5AdBjuCR0s07vAd0Wkf/PzzfXGnmsesM70cpaIzFmlql8A9wNrRORT4B1gUCfv+RgoJXp460/AM6ra3mG5zsr/Mdbeh6r6NdEt5o6SWuu6NxEd2OsToiNybgZq2ig6hugeTmc6i+ctIDnWVz8H4ofJNDpOVyZQrqoVsecS6d926+yKtto/zbrXED2cCtF++7NEx9f6Z+Ch9pYp1v7DwFoR+QRoOTjkFYCNy9VH2OXbxrQiIl5VrRORdKJ7Mj+IJT9zGkRkPPAPqnrzWaxzJfATVd11tuo0zrFzRMac6t9EZDTRcxvPWRI6M7ELH94XEZeqhs+0PomOKrvaklDfYXtExhhjHGXniIwxxjjKEpExxhhHWSIyxhjjKEtExhhjHGWJyBhjjKMsERljjHGUJSJjjDGO+v9u7e9f4VRUjQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Set frequency corners in Hz\n",
    "f1 = 0.1\n",
    "f2 = 1.0\n",
    "\n",
    "# Plot as wiggles\n",
    "wg.pw_wiggles_Audet2016(strf, t1=t1, tmax=10., f1=f1, f2=f2,\n",
    "                        ftitle='audet2016', scale=2.e-7, save=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voilà! \n",
    "\n",
    "You can play with these parameters and change the model file (while respecting its format) and see how the seismograms change. You can also compare the run time difference between the Python and Fortran implementations, but for this simple case with only one plane wave simulation, the difference is minor.\n",
    "\n",
    "## References\n",
    "* Audet, P. (2016). Receiver functions using OBS data: promises and limitations from numerical modelling and examples from the Cascadia Initiative. Geophysical Journal International, 205, 1740-1755. https://doi.org/10.1093/gji/ggw111"
   ]
  }
 ],
 "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
}