gistfile1.txt

{
 "metadata": {
  "name": "effect_of_dispersion.ipynb"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 5,
     "metadata": {},
     "source": [
      "location of the project"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!pwd"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "/home/thomas/RESEARCH/ViscoElasticity/test_cases/elasticInversionOfViscoData\r\n"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "heading",
     "level": 5,
     "metadata": {},
     "source": [
      "Importing necessary packages"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "import pandas as pd"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 16
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Results of the elastic inversion of viscoelastic data"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Viscoelastic data have been inverted with an elastic code. Leading to Vp and Vs results that are not too bab but not as good as the results obtained for an elastic inverison of elastic data performed with the same velocity model and acquisition settings and the same tools (codewise).\n",
      "\n",
      "Here the aim of the study is to check if the elastic results are equivalent to the velocities obtained after correcting the velocity model for dispersion. In viscoelastic theory, the velocity is frequency dependent, i.e. dispersive. The dispersion curves depend on the velocity model defined at a single frequency and attenuation (Liu et al. 1976, Kjartansson 1979).\n",
      "\n",
      "In the definition of dispersion which frequency should be chose ? Should it be the dominant frequency of the source or the dominant frequency of the attenuation mechanism ? If using Kjartansson's constant Q velocity dispersion equation, it makes sense that w0 is the source frequency, as Q is constant with frequency. If using one SLS to model attenuation, Q is then frequency dependent and the it's minimum is set a the resonance frequency of the attenuation mechanism (the one SLS). Then should w0 be fdom or the relaxation mechanism ? I guess Kjartansson's dispersion model is then irrelevant ? Which one should be adopted ? Let's look into Mavko et al.\n",
      "\n"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "True model"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# geometry, acquisition parameters\n",
      "nx, nz = 1501, 301\n",
      "fdom = 7.5\n",
      "w0 = 2*np.pi*fdom\n",
      "# load Vp\n",
      "true_VP_frame = pd.read_table('tb.VP', header=None, sep='\\s+')\n",
      "true_VP = true_VP_frame.values[:nz,2]\n",
      "dz = true_VP_frame.values[1,1] - true_VP_frame.values[0,1]\n",
      "# load Qp\n",
      "true_QP_frame = pd.read_table('tb.QP', header=None, sep='\\s+')\n",
      "true_QP = true_QP_frame.values[:nz,2]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# plot true Vp and Qp models\n",
      "depth = np.arange(0, nz*dz, dz)\n",
      "plt.plot(depth, true_VP, label='P-waves velocity')\n",
      "plt.legend(loc='lower right')\n",
      "plt.figure()\n",
      "plt.plot(depth, true_QP, label='P-waves attenuation')\n",
      "plt.ylim((0,100))\n",
      "plt.legend(loc='best')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 18,
       "text": [
        "<matplotlib.legend.Legend at 0xcdce910>"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAD9CAYAAABTJWtQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9UVHX+P/DnoIJrzslfAX7CARFjGPzBoMOwmTByymPs\nSYz+AA3Wo6yfDXHdUtOMtrA95XEzET1H5FPRcVOk1DimaZhsww83Bzj+aAMkSYy2RXH8kYOoiLy/\nf7jerySC2nB/4PNxzpzkcu+d95tLryfv9/0xOiGEABEREQAPpRtARETqwVAgIiIJQ4GIiCQMBSIi\nkjAUiIhIwlAgIiLJXYXC9evXYTab8cwzzwAAXC4X4uLiYDAYMGPGDDQ3N0vrrlu3DqNHj4bJZEJZ\nWZm0vKamBuHh4QgMDER6erqbu0FERO5wV6GQlZUFk8kEnU4HAMjOzobBYMDx48fh5+eHjRs3AgCa\nmpqwYcMGFBUVITs7GwsXLpT2sXjxYixbtgwVFRUoLi5GZWVlD3SHiIh+jW5D4d///jf27NmDP/zh\nD7h5n1t5eTlSUlLg5eWFuXPnwuFwAAAcDgemTZsGg8GA6OhoCCGkUURtbS0SEhIwdOhQxMfHS9sQ\nEZF6dBsKL730Et555x14ePz/VSsqKmA0GgEARqMR5eXlAG6EQkhIiLRecHAwHA4H6urq4O3tLS03\nmUw4ePCg2zpBRETu0berb+7evRve3t4wm82w2+3S8nt5MsbNKadbdbV9Z+sTEVHX3PXEoi5HCv/8\n5z/x2WefYeTIkZg5cyb+8Y9/IDk5GRaLBTU1NQBunEC2WCwAAKvViurqamn7Y8eOwWKxICgoCKdP\nn5aWV1dXIzIy8o7vK4Tola833nhD8Tawf+wf+9f7Xu7UZSi8/fbb+PHHH1FfX4/8/HzExMTgo48+\ngtVqRW5uLi5fvozc3FypwEdERKCwsBANDQ2w2+3w8PCAXq8HcGOaKT8/H06nEwUFBbBarW7tCBER\n/Xr3dJ/Czamd1NRUNDQ0IDg4GD/99BNeeOEFAICPjw9SU1MRExOD+fPnIysrS9p29erV+Nvf/gaL\nxYLJkydj4sSJbuwGERG5g064e+zxK+l0OrcPh9TCbrfDZrMp3Ywew/5pG/unXe6smwwFIiKNc2fd\n5GMuiIhIwlAgIiJJl/cpEBGpwfLlwJo1SrfiwcBQICLVO3kSeO89ICFB6ZaoU//+7tsXQ4GIVO/S\nJWDQIMDLS+mW9H48p0BEqnfpEjBggNKteDAwFIhI9VpagIceUroVDwaGAhGp3qVLDAW5MBSISPU4\nfSQfhgIRqR5HCvJhKBCR6vGcgnwYCkSkakJw+khODAUiUrWrV4G+fW+8qOcxFIhI1Th1JC+GAhGp\nGqeO5MVQICJV45VH8mIoEJGqcfpIXgwFIlI1Th/Ji6FARKrG6SN5MRSISNU4fSQvhgIRqRqnj+TF\nUCAiVeP0kbwYCkSkai0tHCnIiaFARKrGkYK8GApEpGoMBXkxFIhI1Th9JC+GAhGpGkcK8mIoEJGq\nMRTkxSeUE6nMf/4DNDUp3Qr1OHWK00dyYigQqcz06YDLBfzmN0q3RB08PIDRo5VuxYODoUCkMhcu\nAHv3shCSMnhOgUhleLUNKYmhQKQyPLFKSmIoEKmIEAwFUhZDgUhFrl0DdDqgXz+lW0IPqi5D4cqV\nK7BarQgLC0NkZCQyMzMBABkZGfDz84PZbIbZbMbevXulbdatW4fRo0fDZDKhrKxMWl5TU4Pw8HAE\nBgYiPT29h7pDpG0cJZDSdEII0dUKLS0tGDBgAK5evYoJEyagoKAAeXl50Ov1WLRoUYd1m5qaEBUV\nhX379qG+vh4vvfQSDh06BACIjY3F7Nmz8eSTTyIuLg5r167FxIkTb2+QTodumkTUa/3734DVCvz0\nk9ItIS1xZ93sdvpowH8vg2hubkZbWxu8vLwAoNMGOBwOTJs2DQaDAdHR0RBCoLm5GQBQW1uLhIQE\nDB06FPHx8XA4HG7pAFFvwiuPSGnd3qfQ3t4Os9mMqqoqrF27FgaDAQCwfv16bNu2Dc8++yzmz58P\nvV6P8vJyhISESNsGBwfD4XDA398f3t7e0nKTyYQtW7YgLS2t0/fMyMiQ/m2z2WCz2e6ze0Tawukj\nuht2ux12u71H9t1tKHh4eODo0aM4efIkYmNjMWnSJKSmpuL111/HxYsX8fLLLyMnJwdLlizpdPSg\n0+luW9bdMOfWUCB6kPDziOlu/PKP5RUrVrht33d99VFAQABiY2PhcDjg7e0NnU6Hhx9+GGlpaSgo\nKAAAWK1WVFdXS9scO3YMFosFQUFBOH36tLS8uroakZGRbusEUW/BzyMmpXUZCk6nExcuXAAAnD17\nFvv27UNcXBwaGxsBAG1tbcjLy0NsbCwAICIiAoWFhWhoaIDdboeHhwf0ej0AwGg0Ij8/H06nEwUF\nBbBarT3ZLyJN4vQRKa3L6aPGxkbMnj0b169fh6+vL5YsWYLhw4fj97//PY4cOQJPT09ERUUhNTUV\nAODj44PU1FTExMTA09MTOTk50r5Wr16NpKQkLF++HImJiZ1eeUT0oOP0ESmt20tS5cZLUulB9n//\nB1RUAO+9p3RLSEtkvSSViOTD6SNSGkOBSEU4fURKYygQqQivPiKlMRSIVITTR6Q0hgKRinD6iJTG\nUCBSEU4fkdIYCkQqwukjUhpDgUhFOH1ESmMoEKkIRwqkNIYCkYrw8xRIaQwFIhXhSIGU1u3nKRD1\npG3bgHPnlG6FejidDAVSFh+IR4ppawO8vIA//EHplqjHgAHA3/4G9OundEtIS9xZNxkKpJiffwZG\njAAuXlS6JUTaxqekUq/AG7WI1IehQIrhNflE6sNQIMXwShsi9WEokGIYCkTqw1AgxfBGLSL1YSiQ\nYjhSIFIfhgIphqFApD4MBVIMp4+I1IehQIrhSIFIfRgKpBjep0CkPgwFUgzvaCZSH4YCKYbTR0Tq\nw1AgxXD6iEh9GAqkGE4fEakPQ4EUw+kjIvVhKJBiOH1EpD4MBVIMp4+I1IehQIrh9BGR+jAUSDGc\nPiJSH4YCKYbTR0Tqw1AgxXD6iEh9ugyFK1euwGq1IiwsDJGRkcjMzAQAuFwuxMXFwWAwYMaMGWhu\nbpa2WbduHUaPHg2TyYSysjJpeU1NDcLDwxEYGIj09PQe6g5pCaePiNSny1Do378/vvrqKxw5cgTF\nxcX44IMPcPz4cWRnZ8NgMOD48ePw8/PDxo0bAQBNTU3YsGEDioqKkJ2djYULF0r7Wrx4MZYtW4aK\nigoUFxejsrKyZ3tGqtbWduPl6al0S4joVt1OHw3476Rvc3Mz2tra4OXlhfLycqSkpMDLywtz586F\nw+EAADgcDkybNg0GgwHR0dEQQkijiNraWiQkJGDo0KGIj4+XtqEH081Rgk6ndEuI6FZ9u1uhvb0d\nZrMZVVVVWLt2LQwGAyoqKmA0GgEARqMR5eXlAG6EQkhIiLRtcHAwHA4H/P394e3tLS03mUzYsmUL\n0tLS3N0f1Tp/Hli2DLh2TemWqMOVK5w6IlKjbkPBw8MDR48excmTJxEbG4tJkyZBCHHXb6Dr5E/B\n7rbPyMiQ/m2z2WCz2e76/dSqvh4oLARu6doDb+5cpVtApE12ux12u71H9t1tKNwUEBCA2NhYOBwO\nWCwW1NTUwGw2o6amBhaLBQBgtVqxf/9+aZtjx47BYrFAr9fj9OnT0vLq6mpERkbe8b0yemHlbG8H\nHnkEmDNH6ZYQkdb98o/lFStWuG3fXZ5TcDqduHDhAgDg7Nmz2LdvH+Li4mC1WpGbm4vLly8jNzdX\nKvAREREoLCxEQ0MD7HY7PDw8oNfrAdyYZsrPz4fT6URBQQGsVqvbOqEF7e2ABy8AJiKV63Kk0NjY\niNmzZ+P69evw9fXFkiVLMHz4cKSmpiIpKQnBwcEIDw/HqlWrAAA+Pj5ITU1FTEwMPD09kZOTI+1r\n9erVSEpKwvLly5GYmIiJEyf2bM9UhqFARFqgE/dygkAGOp3uns5ZaMU//wksWXLjv0RE7uTOusm/\nXWVy/TpHCkSkfixTMmlvB/r0UboVRERdYyjIhOcUiEgLWKZkwlAgIi1gmZIJQ4GItIBlSiYMBSLS\nApYpmTAUiEgLWKZkwktSiUgLWKZkwpECEWkBy5RMeJ8CEWkBQ0EmHCkQkRawTMmEoUBEWsAyJROG\nAhFpAcuUTBgKRKQFLFMyYSgQkRawTMmE9ykQkRawTMmEl6QSkRYwFGTC6SMi0gKWKZkwFIhIC1im\nZMJQICItYJmSCUOBiLSAZUomDAUi0gKWKZnwklQi0gKWKZlwpEBEWsAyJRPep0BEWsBQkAlHCkSk\nBSxTMmEoEJEWsEzJhKFARFrAMiUThgIRaQHLlEx4SSoRaQHLlEw4UiAiLWCZkgkvSSUiLWAoyIQj\nBSLSApYpmTAUiEgLWKZkwlAgIi1gmZIJQ4GItKDbMvXjjz9iypQpCA0Nhc1mQ15eHgAgIyMDfn5+\nMJvNMJvN2Lt3r7TNunXrMHr0aJhMJpSVlUnLa2pqEB4ejsDAQKSnp/dAd9SLoUBEWtC3uxX69euH\nzMxMhIWFwel0IiIiAs888wx0Oh0WLVqERYsWdVi/qakJGzZsQFFREerr67Fw4UIcOnQIALB48WIs\nW7YMTz75JOLi4lBZWYmJEyf2TM9UhvcpEJEWdBsKvr6+8PX1BQAMGzYMoaGhqKioAAAIIW5b3+Fw\nYNq0aTAYDDAYDBBCoLm5GQMHDkRtbS0SEhIAAPHx8XA4HA9MKPCSVCLSgnv627Wurg5VVVWwWq0A\ngPXr1yMyMhKrVq2Cy+UCAJSXlyMkJETaJjg4GA6HA3V1dfD29paWm0wmHDx40B190AROHxGRFnQ7\nUrjJ5XIhISEBmZmZeOihh5CamorXX38dFy9exMsvv4ycnBwsWbKk09GDTqe7bVln692UkZEh/dtm\ns8Fms91tM1WLoUBE7mK322G323tk3zrRVXX+r2vXruF3v/sdYmNj8eKLL972/aNHj2L+/Pk4cOAA\ndu3ahf379yMrKwsAEBYWhtLSUuj1egQGBuLEiRMAgHfffRf9+/dHWlpaxwbpdF0GhlYtWAAYjTf+\nS0TkTu6sm93+7SqEQEpKCsaMGdMhEBobGwEAbW1tyMvLQ2xsLAAgIiIChYWFaGhogN1uh4eHB/R6\nPQDAaDQiPz8fTqcTBQUF0jTUg4AjBSLSgm6njw4cOIDNmzdj3LhxMJvNAIC3334bW7duxZEjR+Dp\n6YmoqCikpqYCAHx8fJCamoqYmBh4enoiJydH2tfq1auRlJSE5cuXIzEx8YE5yQwwFIhIG+5q+khO\nvXX6aN48wGIB/vd/lW4JEfU2sk4fkXtwpEBEWsAyJRPep0BEWsBQkAlHCkSkBSxTMmEoEJEWsEzJ\nhKFARFrAMiUThgIRaQHLlEz4lFQi0gKWKZlwpEBEWsAyJRNekkpEWsBQkAlHCkSkBSxTMmEoEJEW\nsEzJhKFARFrAMiUThgIRaQHLlEwYCkSkBSxTMuF9CkSkBSxTMuElqUSkBQwFmXD6iIi0gGVKJgwF\nItIClimZMBSISAtYpmTCUCAiLWCZkglDgYi0gGVKJrwklYi0gGVKJhwpEJEWsEzJhPcpEJEWMBRk\nwpECEWkBy5RMGApEpAUsUzJhKBCRFrBMyYShQERawDIlE16SSkRawDIlE44UiEgLWKZkwktSiUgL\nGAoy4UiBiLSAZUomDAUi0gKWKZkwFIhIC1imZMJQICItYJmSCUOBiLSg2zL1448/YsqUKQgNDYXN\nZkNeXh4AwOVyIS4uDgaDATNmzEBzc7O0zbp16zB69GiYTCaUlZVJy2tqahAeHo7AwECkp6f3QHfU\ni/cpEJEWdFum+vXrh8zMTFRVVWH79u147bXX4HK5kJ2dDYPBgOPHj8PPzw8bN24EADQ1NWHDhg0o\nKipCdnY2Fi5cKO1r8eLFWLZsGSoqKlBcXIzKysqe65nK8JJUItKCbkPB19cXYWFhAIBhw4YhNDQU\nFRUVKC8vR0pKCry8vDB37lw4HA4AgMPhwLRp02AwGBAdHQ0hhDSKqK2tRUJCAoYOHYr4+HhpmwcB\np4+ISAv63svKdXV1qKqqQkREBObMmQOj0QgAMBqNKC8vB3AjFEJCQqRtgoOD4XA44O/vD29vb2m5\nyWTCli1bkJaWdtv7ZGRkSP+22Wyw2Wz30kxVYigQkbvY7XbY7fYe2fddh4LL5UJCQgIyMzMxcOBA\nCCHu+k10Ot1ty7ra/tZQ6C0YCkTkLr/8Y3nFihVu2/ddlalr167hueeeQ3JyMuLi4gAAFosFNTU1\nAG6cQLZYLAAAq9WK6upqadtjx47BYrEgKCgIp0+flpZXV1cjMjLSbR1RO4YCEWlBt2VKCIGUlBSM\nGTMGL774orTcarUiNzcXly9fRm5urlTgIyIiUFhYiIaGBtjtdnh4eECv1wO4Mc2Un58Pp9OJgoIC\nWK3WHuqW+jAUiEgLdKKbeaCysjJERUVh3Lhx0jTQypUrMWnSJCQlJeHw4cMIDw/H5s2bMXDgQABA\nVlYW1q9fD09PT+Tk5GDy5MkAbowOkpKScP78eSQmJmLlypW3N0inu6epKa3w8gJ+/hno31/plhBR\nb+POutltKMitt4ZCv37ApUuAp6fSLSGi3saddZMTGjLhfQpEpAUMBZnwnAIRaQHLlAxujuo6uTKX\niEhVGAoy4CiBiLSCpUoGDAUi0gqWKhnwCalEpBUsVTLgSIGItIKlSga8HJWItIKhIAOOFIhIK1iq\nZMBQICKtuKfPU6D7w1CgWw0ZMgTnz59XuhmkQYMHD8a5c+d69D0YCjJgKNCtzp8/3yuf70U9r7PP\npnE3lioZMBSISCtYqmTA+xSISCtYqmTAS1KJSCsYCjLg9BERaQVLlQwYCkTa5eHhgRMnTvyqfYwZ\nMwYlJSVualHPYqmSAUOBtCQgIAADBgyAXq+HxWLBX/7yF7S2tirdLE379ttvERUVBQDIyMhAcnKy\nwi26M5YqGTAUSEt0Oh12794Nl8uF999/H5s2bcKnn36qdLNIJixVMmAokFaNHz8eTz/9NHbv3n3b\n9+rr6zF48GDp63nz5sHHx0f6Ojk5GVlZWQCADz/8ECaTCYMGDcKMGTOwa9cuab2QkBB8/vnn0tdt\nbW145JFHcOTIEQDA999/j6VLl8Lf3x/z5s1DdXW1tO7OnTsxZcoUDBo0CIGBgcjLy7utnf/5z38w\nYMCADjcMHj58GI888giuX78OACgrK8Pzzz+PkSNHYsWKFXA6nZ3+PFpaWpCdnY2xY8di6tSpHfoB\nAF999RWSkpIwZMgQhIaGSn0ICAhAUVERvvjiC6xcuRIff/wx9Ho9zGYztm/fjokTJ3bYz5o1azBj\nxoxO29DjhMqosEm/Wk2NEI89pnQrSC3U/jseEBAg9u/fL4QQ4tChQ2LEiBFix44dna5rMBjEoUOH\nhBBCPPbYY2LUqFGipqZG+t6RI0eEEEJ8/vnn4sSJE6K1tVVs3bpV/OY3vxGXLl0SQgjx5ptviuef\nf17a5+7du4XJZBJCCNHW1ia8vb3Fhx9+KC5evCg2bdok/Pz8hBBCtLa2Cn9/f3Hw4EEhhBCnTp0S\nVVVVnbYzJiZGvPfee9LXS5YsEampqUIIIY4ePSr8/PzEl19+Kc6dOyf+9Kc/iVmzZknr6nQ68f33\n3wshhHj99dfFlClTRG1trSgqKhIBAQHiq6++kn5W3t7eIi8vT7S2toq6ujrxww8/SD/ToqIiIYQQ\nGRkZIjk5Wdr/1atXxZAhQ6SfmxBChIWFiU8//fS2ftzpd8edv1Oq++1U+/8w96OqSgijUelWkFrc\nze/4jQ9x/XWv++Xv7y8GDhwoBg8eLCZPnixWrlwp2traOl03OTlZrFmzRjQ2Norg4GCxbNkysXHj\nRnHixAkxaNCgO77HE088IbZv3y6EEOL48eNCr9eLy5cvCyGEmDVrlvjrX/8qhBBi37594qmnnuqw\nbVhYmCgvLxfXrl0Tw4cPF5s3b5YC5k7ef/99ERMTI4QQor29XYwYMUKUlpYKIYR49dVXxVtvvSWt\n63Q6xbBhw6Q+3xoK48ePF4WFhdK66enpYuHChUIIIZYuXSr+/Oc/d/r+t4bCG2+8IZKSkjp8/4UX\nXhDp6elCCCG+/fZbMXjwYNHa2nrbfuQIBU5qyID3KdC9ckcs3C+dToedO3fi3LlzKCkpwSuvvII+\nffrg7bffhl6vh16vx/z58wEA0dHRsNvtKC0tRVRUFKKjo1FcXIySkhJMnjxZ2mdZWRlmzpwJf39/\nDBo0COXl5fjmm28AAEFBQQgJCcFnn32GlpYW7Nq1C7NmzQIA7N+/H6WlpRg8eLD0qqurQ0lJCfr2\n7YsdO3Zg+/bt8PPzQ0pKCurr6zvtU3x8PL7++mucOnUKJSUl8PDwwBNPPCG9x8qVK6X9BwUFoaWl\nBYcOHeqwD5fLhW+++QYTJkyQlk2YMAGlpaUAALvdjkmTJt3Xz3z27NnS1NdHH32EhIQE9OvX7772\n9WsxFGTAcwrUG7z66qtwuVxwuVzYsGEDgBuhUFpaCrvdDpvNhieeeAIHDhxAcXExbDYbAEAIgT/+\n8Y+Ijo7GoUOHcOHCBURERHR4/tPMmTOxdetW7Ny5EyaTCYGBgQCAmJgY2Gw2nD9/Xnq5XC4sXrwY\nAPDb3/4WBQUFOHnyJPr164elS5d22vbBgwdj6tSp+Pjjj5GXl4eZM2dK34uJicFrr73W4T0uXboE\ni8XSYR96vR7jxo1DZWWltKyyslK6qmjKlCkoKyvr9ufYt2/f2559FRkZCU9PT5SUlGDr1q2KXp3E\nUiUDhgL1VkFBQejfvz82b96M6Oho6PV6eHt7Y8eOHYiOjgYAtLa24syZM/Dx8UH//v3x4YcfwuFw\ndNhPYmIiCgsLsXHjRjz//PPS8ieffBL/+te/8Pe//x3nz5/HlStXYLfb8dNPP6GpqQk7d+7EpUuX\n0KdPH/Tv3x96vf6ObZ01axY2bdqEHTt2SCMR4MYJ8ZycHOzbtw+tra34+eefsW3btk73ERcXh3fe\neQffffcd7HY7tm7dKp0QTkxMxCeffIJPPvkEra2tqKurQ0NDw237mDBhAqqrq3H16tUOy5OTk7Fg\nwQJ4enri8ccf7+Yn33NYqmTAUKDezGazYdiwYXj00UelrwEgPDwcAODl5YWsrCy8+eabCAoKQmVl\nJRITEzvsw9fXF48//ji+/vprJCQkSMv79OkDu92O2tpaTJgwAQaDAe+++y6EEGhvb0dmZiYeffRR\nGI1GnDt3DitWrLhjO6dPn466ujoMHz4cY8eOlZabTCZs2rQJn3zyCfz8/DB27FgUFhZK37/1yaRL\nly7FjBkzEB8fj7feegtr1qyRwi8sLAxbtmzBp59+Cm9vb8THx3f6iPTo6Gg89thjGDlyZIerjpKT\nk1FVVYWkpKRuf+Y9SSd+OY5RmE6n63WPFa6sBF544cZ/iXrj7zj9ei0tLfD19cXhw4cxatSoTte5\n0++OO3+n+PerDDhSIKLurF27Fk899dQdA0Eu/JAdGfDR2UTUlYCAAAwZMgT5+flKN4WhIAdekkpE\nXTl58qTSTZDw71cZcPqIiLSCpUoGDAUi0gqWKhkwFIhIK3hOQQYMBbrV4MGDO1z7TnS3bn0qbU9h\nKMiAoUC3OnfunNJNILqjbkvV3Llz4ePj0+EOwIyMDPj5+cFsNsNsNmPv3r3S99atW4fRo0fDZDJ1\neA5ITU0NwsPDERgYiPT0dDd3Q91uXpJqt9uVbkqPYv+0jf0j4C5CYc6cOfjiiy86LNPpdFi0aBEO\nHz6Mw4cP4+mnnwYANDU1YcOGDSgqKkJ2djYWLlwobbN48WIsW7YMFRUVKC4u7vBQqd7u5iWpvf2X\nkv3TNvaPgLsIhcmTJ3c6j9XZLdUOhwPTpk2DwWBAdHQ0hBBobm4GANTW1iIhIQFDhw5FfHz8bQ/E\n6s04fUREWnHf5xTWr1+Pbdu24dlnn8X8+fOh1+tRXl6OkJAQaZ3g4GA4HA74+/vD29tbWm4ymbBl\nyxakpaV1uu9nnrnfVqnTqVPA//yP0q0gIroLd/NJPPX19WLMmDHS16dPnxbt7e3iwoULYt68eeKd\nd94RQtz4FKKNGzdK6yUkJIiioiJx/PhxERkZKS3fs2fPbZ88dBMAvvjiiy++7vHlLvc1Urj5V//D\nDz+MtLQ0zJ8/H0uWLIHVasX+/ful9Y4dOwaLxQK9Xo/Tp09Ly6urqxEZGdnpvgWfHklEpJj7mulu\nbGwEALS1tSEvLw+xsbEAgIiICBQWFqKhoQF2ux0eHh7Sh14YjUbk5+fD6XSioKAAVqvVTV0gIiJ3\n6XakMHPmTBQXF8PpdGLEiBFYsWIF7HY7jhw5Ak9PT0RFRSE1NRUA4OPjg9TUVMTExMDT0xM5OTnS\nflavXo2kpCQsX74ciYmJHT5cgoiIVMJtE1G/UnFxsTAajSIoKEisW7dO6ebcN39/fzF27FgRFhYm\nLBaLEEKIixcviunTp4sRI0aIuLg44XK5pPWzsrJEUFCQCAkJEaWlpUo1u1Nz5swR3t7eHc4n3U9f\nqqurhdlsFiNHjhSvvvqqrH3oSmf9e+ONN8Sjjz4qwsLCRFhYmNizZ4/0Pa31r6GhQdhsNmEymUR0\ndLTYsmWLEKL3HMM79a83HMPLly+LiIgIMX78eGG1WsWaNWuEEPIcO9WEQlhYmCguLhYnT54UwcHB\n4syZM0o36b4EBASIs2fPdli2atUqsWDBAnHlyhWRlpYmnZg/ffq0CA4OFj/88IOw2+3CbDYr0eQ7\nKikpEYcOHepQNO+nL08//bTIz88XTqdTTJo0SVRUVMjel8501r+MjAzx7rvv3rauFvvX2NgoDh8+\nLIQQ4syZM2LkyJHi4sWLveYY3ql/veUYXrp0SQghxJUrV0RoaKj47rvvZDl2qrh6/ueffwYAREVF\nwd/fH1OnTtX0fQziFyfLy8vLkZKSAi8vL8ydO1fqW2f3dbhcLiWa3KnO7lG5l76o/R6V3n4Pjq+v\nL8LCwgCcR7EtAAADAElEQVQAw4YNQ2hoKCoqKnrNMbxT/4DecQwHDBgAAGhubkZbWxu8vLxkOXaq\nCIWKigoYjUbpa5PJhIMHDyrYovun0+kQExODGTNm4LPPPgPQsX9GoxHl5eUAbhzIX97XcfN7anUv\nfXE4HKirq7vtHhW1H9v169cjMjISq1atkkL6TvfgaKV/dXV1qKqqQkRERK88hjf7d/MClt5wDNvb\n2zF+/Hj4+PhgwYIFMBgMshw7VYRCb3LgwAEcPXoUK1euxKJFi3Dq1Kl7usxW7U/P/LV9uZftlZCa\nmor6+noUFhbi+++/ly6W6KzdWumfy+VCQkICMjMzMXDgwF53DG/t30MPPdRrjqGHhweOHj2Kuro6\nbNiwAYcPH5bl2KkiFCwWC44dOyZ9XVVVdcf7GNRu+PDhAICQkBBMnz4du3btgsViQU1NDYAbDwa0\nWCwAAKvViurqamnbm/d1qNm99iUoKOiu71FRA29vb+h0OukenIKCAgDa7d+1a9fw3HPPITk5GXFx\ncQB61zHsrH+97RgGBAQgNjYWDodDlmOnilB4+OGHAQAlJSU4efIkvvzyS03ex9DS0iINVc+cOYPC\nwkJMmzYNVqsVubm5uHz5MnJzc6WD0tV9HWp1P33R0j0qvekeHCEEUlJSMGbMGLz44ovS8t5yDO/U\nv95wDJ1OJy5cuAAAOHv2LPbt24e4uDh5jp07zpK7g91uF0ajUYwaNUpkZWUp3Zz7cuLECTF+/Hgx\nfvx4ERMTIz744AMhRNeXka1du1aMGjVKhISEiJKSEqWa3qnExEQxfPhw4enpKfz8/ERubu599aWq\nqkqYzWYREBAgXnnlFSW60qmb/evXr5/w8/MTH3zwgUhOThZjx44VEyZMEC+99FKHK8m01r/S0lKh\n0+nE+PHjpcsz9+7d22uOYWf927NnT684ht98840wm81i3LhxYurUqWLTpk1CiPurJffaN50QKplA\nIyIixali+oiIiNSBoUBERBKGAhERSRgKREQkYSgQEZGEoUBERJL/B15wTQ0kLNM7AAAAAElFTkSu\nQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x440eb10>"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD9CAYAAAC7iRw+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHA9JREFUeJzt3XlQFGfeB/DvoAaTiIIa0RUHVlEGUA4tDhVlJK6rVkVY\nrawaZU1gLSVaHnu4KYMra1KyxkQ0yQp5s47xiForauEZIvs6gq7O4BW3cIzgERRRIXgMAovI8/7B\na68kaGBopofp76dqCrp7uvv38Ex9Z3i6p1sjhBAgIiLVcFG6ACIisi8GPxGRyjD4iYhUhsFPRKQy\nDH4iIpVh8BMRqcxzgz8hIQGenp4YMmSINM9qtSI2NhZarRZxcXGorKyUln388ccYOHAgAgICcOzY\nsbarmoiIbPbc4H/rrbfw1VdfNZqXnp4OrVaLwsJCeHl5ISMjAwBw584drF+/Hv/85z+Rnp6OBQsW\ntF3VRERks+cG/6hRo+Dh4dFontlsRmJiIlxdXZGQkACTyQQAMJlMGD9+PLRaLaKjoyGEgNVqbbvK\niYjIJi0e48/Pz4dOpwMA6HQ6mM1mAA3B7+/vLz3Pz89PWkZERI6jY0tXaMkVHjQaTbPmERHRT5Pr\nCjst/sQfFhYGi8UCALBYLAgLCwMARERE4MKFC9LzLl68KC37ISGEUz7+938FfHyWK15HWz6WL2f7\n2vPDmdvnzG0TQt5LqrU4+CMiImAwGFBdXQ2DwYDIyEgAQHh4OLKzs1FcXAyj0QgXFxe4ubnJWiwR\nEbXec4N/+vTpGDFiBC5duoR+/fph48aNSEpKQnFxMfz8/FBSUoK5c+cCADw9PZGUlISYmBi8/fbb\nWLdunV0aQERELfPcMf7t27c3OT8rK6vJ+QsXLsTChQtbX1U75u6uV7qENqXX65UuoU2xfe2XM7dN\nbhoh9+DRT+1Qo5F9vMpRHDkCrFjR8JOISE5yZicv2UBEpDIMfiIilWHwExGpDIOfiEhlGPxERCrD\n4CciUhkGPxGRyjD4iYhUhsFPRKQyDH4iIpVh8BMRqQyDn4hIZRj8REQqw+AnIlIZBj8Rkcow+ImI\nVIbBT0SkMgx+IiKVYfATEakMg5+ISGUY/EREKsPgJyJSGQY/EZHKMPiJiFSGwU9EpDIMfiIilWHw\nExGpDIOfiEhlGPxERCrD4CciUhkGPxGRyjD4iYhUhsFPRKQyDH4iIpVh8BMRqYzNwf/5559jxIgR\nGDZsGBYtWgQAsFqtiI2NhVarRVxcHCorK2UrlIiI5GFT8FdUVGDlypU4fPgw8vPzcenSJWRnZyM9\nPR1arRaFhYXw8vJCRkaG3PUSEVErdbRlpRdffBFCCNy/fx8AUFVVBXd3d5jNZiQnJ8PV1RUJCQlI\nTU2VtVgicl7FxUBhodJVqIPNwZ+eng4fHx+4urpiwYIFiIiIQH5+PnQ6HQBAp9PBbDY3uX5KSor0\nu16vh16vt6UMInIiixYBly8DPXsqXYljuHvXiHv3jG2ybZuCv6ysDElJSbhw4QI8PDzw+uuvY//+\n/RBCNGv9p4OfiAgA7t8H1qwBXn1V6Uochf7/Hw00mr/ItmWbxvjNZjMiIyPh6+uLHj164PXXX0de\nXh7CwsJgsVgAABaLBWFhYbIVSkTOzWoF3NyUrkIdbAr+UaNG4dSpU6ioqMB//vMfHDp0COPGjUNE\nRAQMBgOqq6thMBgQGRkpd71E5KSsVqBLF6WrUAebgr9r165ITk7Gr371K0RFRSE4OBhjxoxBUlIS\niouL4efnh5KSEsydO1fueonISfETv/1oRHMH5uXaoUbT7GMB7c2RI8CKFQ0/iahlunUDvvsOcHdX\nuhLHJGd28pu7RKQ4IYDKSg712AuDn4gUV10NvPAC0NGm8wyppRj8RKQ4ju/bF4OfiBTH4LcvBj8R\nKY7Bb18MfiJSHIPfvhj8RKQ4ntFjXwx+IlIcP/HbF4OfiBTH4LcvBj8RKY7Bb18MfiJSHC/QZl8M\nfiJSXGUlP/HbE4OfiBTHoR77YvATkeIY/PbF4CcixXGM374Y/ESkOH7ity9eBJVIAXo9cPas0lU4\njocPgbVrla5CPRj8RAooKABOnwZeeUXpShxDhw4c6rEnBj+RAqxWoG9f4MUXla6E1Ihj/ER29ugR\nUFcHdO6sdCWkVgx+Ijt7cgaLRqN0JaRWDH4iO+O3VElpDH4iO+Opi6Q0Bj+RnTH4SWkMfiI747dU\nSWkMfiI74xg/KY3BT2RnHOohpTH4ieyMwU9KY/AT2RmDn5TG4CeyM47xk9IY/ER2xrN6SGkMfiI7\n41APKY3BT2RnDH5SGoOfyM4Y/KQ0Bj+RnfHgLinN5uB/+PAhZs2ahUGDBiEgIAAmkwlWqxWxsbHQ\narWIi4tDZWWlnLUSOQUe3CWl2Rz8y5cvh1arxfnz53H+/HnodDqkp6dDq9WisLAQXl5eyMjIkLNW\nIqfAoR5Sms3Bn5OTg6VLl6Jz587o2LEjunXrBrPZjMTERLi6uiIhIQEmk0nOWomcAoOflGbTPXdv\n3LiBmpoaJCUlwWKxYPLkyViwYAHy8/Oh0+kAADqdDmazWdZiqf2pqwMWL24Y16YG9+4x+ElZNgV/\nTU0NLl26hNWrV2Ps2LGYM2cO/vGPf0AI0az1U1JSpN/1ej30er0tZVA7UF4ObN4MrF2rdCWO47XX\neJN1+mlGoxFGo7FNtq0RzU3rH/D394fFYgEAHDp0CJs3b0ZtbS2Sk5MRGhqK06dPIzU1FZmZmY13\nqNE0+w2ivTlyBFixouEnNSgsBMaPBy5fVroSovZNzuy0eYx/4MCBMJlMqK+vx4EDBzB27FhERETA\nYDCguroaBoMBkZGRshRJ7RfHs4kcj01DPQDw4Ycf4je/+Q1qamowduxYTJs2DfX19Zg5cyb8/Pww\ndOhQrFq1Ss5aqR3iOetEjsfm4B80aBBOnjz5o/lZWVmtKoicCz/xEzkefnOX2hSDn8jxMPipTTH4\niRwPg5/aFC9PQOR4GPzUpnhwl8jxMPipTXGoh8jxMPipTTH4iRwPg5/aFIOfyPEw+KlNVVby4C6R\no2HwU5viJ34ix8PgpzbF4CdyPAx+alMMfiLHw+CnNsXz+IkcD4Of2hQ/8RM5HgY/tRkhGj7xv/yy\n0pUQ0dMY/NRmqqoAV1ego80X/yaitsDgpzbDYR4ix8TPYjIrLwf27VO6Csdw6xaDn8gRMfhlNGBA\nw+N//kfpShzHG28oXQER/ZBGyHXb9ubuUMY7xRMRqYWc2ckxfiIilWHwExGpDIOfiEhlGPxERCrD\n4CciUhkGPxGRyjD4iYhUhsFPRKQyDH4iIpVh8BMRqQyDn4hIZRj8REQqw+AnIlIZBj8Rkcow+ImI\nVIbBT0SkMgx+IiKVaVXwP378GKGhoXjttdcAAFarFbGxsdBqtYiLi0NlZaUsRRIRkXxaFfzr1q1D\nQEAANBoNACA9PR1arRaFhYXw8vJCRkaGLEUSEZF8bA7+Gzdu4ODBg/jtb38r3QfSbDYjMTERrq6u\nSEhIgMlkkq1QIiKSh83Bv3jxYqxevRouLv/dRH5+PnQ6HQBAp9PBbDa3vkIiIpJVR1tW2r9/P3r1\n6oXQ0FAYjUZpfnPvAJ+SkiL9rtfrodfrbSmDiMhpGY3GRvkqJ41oblo/ZenSpdiyZQs6duyImpoa\nPHjwAJMnT0ZVVRWSk5MRGhqK06dPIzU1FZmZmY13qNE0+w2CiIgayJmdNg31rFy5EtevX8fVq1ex\nY8cOxMTEYMuWLYiIiIDBYEB1dTUMBgMiIyNlKZKIiOQjy3n8T87qSUpKQnFxMfz8/FBSUoK5c+fK\nsXkiIpKRTUM9rdohh3qIiFpM8aEeIiJqvxj8REQqw+AnIlIZBj8Rkcow+ImIVIbBT0SkMgx+IiKV\nYfATEakMg5+ISGUY/EREKsPgJyJSGQY/EZHKMPiJiFSGwU9EpDIMfiIilWHwExGpDIOfiEhlGPxE\nRCrD4CciUhkGPxGRyjD4iYhUhsFPRKQyDH4iIpVh8BMRqQyDn4hIZRj8REQqw+AnIlIZBj8Rkcow\n+ImIVIbBT0SkMgx+IiKVYfATEakMg5+ISGUY/EREKsPgJyJSGQY/EZHK2BT8169fx5gxYxAYGAi9\nXo9t27YBAKxWK2JjY6HVahEXF4fKykpZiyUiotazKfg7deqEtLQ0FBQUIDMzE8nJybBarUhPT4dW\nq0VhYSG8vLyQkZEhd71ERNRKNgV/7969ERISAgDo2bMnAgMDkZ+fD7PZjMTERLi6uiIhIQEmk0nW\nYomIqPU6tnYDRUVFKCgoQHh4ON566y3odDoAgE6ng9lsbnKdlJQU6Xe9Xg+9Xt/aMoiInIrRaITR\naGyTbWuEEMLWla1WK/R6Pf785z9LY/uXLl1C586dUVVVBX9/f3z33XeNd6jRoBW7JCJSJTmz0+az\neh49eoQpU6YgPj4esbGxAICwsDBYLBYAgMViQVhYmCxFEhGRfGwKfiEEEhMTMXjwYCxatEiaHxER\nAYPBgOrqahgMBkRGRspWKBERycOmoZ5jx45h9OjRCAoKgkajAQCkpqZi5MiRmDlzJs6ePYuhQ4di\n69at6NKlS+MdcqiHiKjF5MzOVo3x27RDBj8RUYs5xBg/ERG1Twx+IiKVYfATEakMg5+ISGUY/ERE\nKsPgJyJSGQY/EZHKMPiJiFSGwU9EpDIMfiIilWHwExGpDIOfiEhlGPxERCrT6lsvyqV79+64e/eu\n0mUQNeLh4YGKigqlyyCSlcNclpmXayZHxNclOQpelpmIiGzG4CciUhkGPxGRyjD4iYhUhsFPTiMp\nKQnvv/++0mUQOTwGfzP5+PjgpZdegpubG8LCwrBs2TLU1tYqXVabcnFxwZUrV6Rpo9GIfv36KVjR\nf33xxRcYNWpUo3np6elITk5WqCKi9oPB30wajQb79++H1WrF3//+d2zatAm7d+9Wuqw2x1MZiZwP\ng98GwcHBmDBhAvbv3/+jZVevXoWHh4c0PXv2bHh6ekrT8fHxWLduHQBg48aNCAgIgLu7O+Li4rBv\n3z7pef7+/jhw4IA0XVdXh1deeQXnzp0DAFy+fBlLliyBt7c3Zs+ejQsXLkjPzcrKwpgxY+Du7o7+\n/ftj27ZtTbbDbDZj+PDh8PDwwPDhw/Hpp5+irq4OADB69GiprV27dsXmzZsxceJE3Lx5E25ubuja\ntStu3boFADh06BAmTZoEPz8/pKWlobKyEgBw7do1uLi4YPfu3fD390dQUBC2bt0q7f/NN9/EsmXL\npOkf/kfx17/+Fb6+vujRowdmzJiBvLw8AIDFYkFSUhJOnDgBNzc3dO/e/ZnbmzJlCgYOHIgPPvig\n0RcEXVxc8OWXXyIkJAS+vr5IS0tr8m9E5JSEnT1rlwqU0iI+Pj4iJydHCCHEmTNnRL9+/cSuXbua\nfK5WqxVnzpwRQggxaNAgMWDAAGGxWKRl586dE0IIceDAAXHlyhVRW1srtm/fLl588UXx8OFDIYQQ\nK1asEDNmzJC2uX//fhEQECCEEKKurk706tVLbNy4UTx48EBs2rRJeHl5CSGEqK2tFd7e3uLkyZNC\nCCFu3bolCgoKmqzz9OnTwmQyibq6OnH8+HHh7e0tDh8+LC3XaDTi8uXL0rTRaJT280RWVpYICgoS\nJ06cEDdv3hS//vWvxdKlS4UQQly9elVoNBoxdepUUVxcLLKzs4Wrq6uorq4WQgjx5ptvimXLlknb\nOnLkSKPt79y5U5SWloqqqiqxZs2aRsu++OILERUV1aiWp7d35coV4e7uLnbu3Clu3Lghpk6dKmbN\nmtWobTExMcJisYjTp08LNzc3UVRU9KO/kaO/Lkk95HwttqvgB+R52MLb21t06dJFeHh4iFGjRonU\n1FRRV1fX5HPj4+PFmjVrRGlpqfDz8xN/+tOfREZGhhRGzxIVFSUyMzOFEEIUFhYKNzc3KSTfeOMN\n8d577wkhhPj666/FL37xi0brhoSECLPZLB49eiT69Okjtm7dKr2JNNe7774r5s+fL03/MPh/GMxP\n6vryyy+l6bNnz0pvUE+C//Tp09JyPz8/8dVXXwkhGoI6OTn5udt/or6+XvTr10+cOnVKCCHExo0b\nnxv8a9asEdOnT5eWFRYWih49eojHjx9LbXv6jfuXv/ylSE9P/9F+GfzkKOR8LbaroR65ot8WGo0G\nWVlZqKioQG5uLt555x106NABK1euhJubG9zc3PD2228DAKKjo2E0GpGXl4fRo0cjOjoaR48eRW5u\nbqMDkseOHcP06dPh7e0Nd3d3mM1mnD9/HgDg6+sLf39/7N27F1VVVdi3bx/eeOMNAEBOTg7y8vLg\n4eEhPYqKipCbm4uOHTti165dyMzMhJeXFxITE3H16tUm21RSUoK5c+ciKCgIXbt2RVpamrT/5srJ\nyUFSUpJUx5gxY3Dt2jXcuXNHek5ISIj0e58+fVBSUtKsbe/duxeTJ0/Gz372M3Tv3h2lpaXNru9f\n//oXhg0bJk37+vqirq4OBQUFra6LqL1rV8HviJYuXQqr1Qqr1Yr169cDaAj+vLw8GI1G6PV6REVF\n4fjx4zh69Cj0ej2AhoOmc+bMQXR0NM6cOYN79+4hPDy80cHU6dOnY/v27cjKykJAQAD69+8PAIiJ\niYFer8fdu3elh9Vqxe9//3sAwPDhw7Fnzx5cu3YNnTp1wpIlS5qs/f3338ejR49w8OBB3L9/H4sX\nL0Z9fb203MXFpVE9HTp0+NHB3piYGHz++eeNann48CF69er1k3+7vn374vbt29L02bNnpd8fPnyI\n2bNnY9asWbh48SIqKirQt29faf9N1fK0kSNH4tSpU9J0YWEhOnTogMDAwGeuo9FofrJmImfA4G8D\nvr6+6Ny5M7Zu3Yro6Gi4ubmhV69e2LVrF6KjowEAtbW1KCsrg6enJzp37oyNGzfCZDI12s60adOQ\nnZ2NjIwMzJgxQ5o/duxY/Pvf/8bmzZtx9+5d1NTUwGg0oqSkBHfu3EFWVhYePnyIDh06oHPnznBz\nc2uyzps3b6J79+7o0aMHjEYjNm/e3Gj5sGHDGoVncHAwysvLUVpaKs2Lj4/HBx98gGPHjuHx48co\nKyvD3r17m/V3evXVV3H48GEUFhbi1KlT2LRpk7TMarWisrISffr0QX19PVJTU3Hz5s1GtRUWFkoH\nkoGGN9MnbwaTJk1CdnY2du/ejZKSEixfvhyvvfYaXFyafsk/vS6Rs2PwtxG9Xo+ePXuib9++0jQA\nDB06FADg6uqKdevWYcWKFfD19cWpU6cwbdq0Rtvo3bs3RowYgRMnTmDq1KnS/A4dOsBoNOLbb7/F\nsGHDoNVq8dFHH0EIgfr6eqSlpaFv377Q6XSoqKjAX/7ylyZrTElJwblz5+Dl5YXVq1dj/vz5jT71\n/uEPf8CHH34IDw8PZGZmomvXrliyZAlGjx6N7t2749atW5gwYQJWrFiBTz/9FK+88gqGDx8Os9ks\nbeN5n6KjoqIwc+ZMvPrqq1i4cCHmzZsnPb93795ITU1FfHw8goODUVtbi6ioKGndgIAAxMXFITAw\nUPrvQqPRSOv3798fO3fuxJYtWxAdHY2goCCsWbPmmXU9vS6Rs+NlmYmeg69LchS8LDMREdmMwU9E\npDIMfiIilXGYe+56eHjw4Bo5nKcvv0HkLBzm4C4RET2bQx/czc3Nhb+/PwYOHIhPPvlE7s07PKPR\nqHQJbYrta9+cuX3O3Da5yR78CxcuxGeffYacnBz87W9/Q3l5udy7cGjO/uJj+9o3Z26fM7dNbrIG\n//379wE0XNLX29sb48aN+9G3UYmISFmyBn9+fj50Op00HRAQgJMnT8q5CyIiaiVZD+7m5ORgw4YN\n2L59OwAgIyMDJSUleO+99/67Q565Q0RkE7niWtbTOcPCwvDHP/5Rmi4oKMD48eMbPYdn9BARKUvW\noZ5u3boBaDiz59q1azh8+DAiIiLk3AUREbWS7F/gWrt2LebMmYNHjx5hwYIF6Nmzp9y7ICKiVpD9\ndM7o6GhYLBYUFRVhwYIFjZY5wzn+Pj4+CAoKQmhoKMLDwwE0XDs+NjYWWq0WcXFxja4R//HHH2Pg\nwIEICAjAsWPHlCr7mRISEuDp6YkhQ4ZI82xpj8ViwdChQ9G/f3+8++67dm3D8zTVvpSUFHh5eSE0\nNBShoaE4dOiQtKy9te/69esYM2YMAgMDodfrsW3bNgDO0YfPapuz9F9NTQ0iIiIQEhKCyMhIpKWl\nAbBT38l2E8dmCAkJEUePHhXXrl0Tfn5+oqyszJ67l4WPj4/4/vvvG81btWqVmD9/vqipqRHz5s0T\nq1evFkIIcfv2beHn5ye+++47YTQaRWhoqBIlP1dubq44c+aMGDx4sDTPlvZMmDBB7NixQ5SXl4uR\nI0eK/Px8u7elKU21LyUlRXz00Uc/em57bF9paak4e/asEEKIsrIy8fOf/1w8ePDAKfrwWW1zpv57\ncl/smpoaERgYKC5dumSXvrPbRdqc6Rx/8YMD1GazGYmJiXB1dUVCQoLULpPJhPHjx0Or1SI6OhpC\nCFitViVKfqZRo0b96Ho0LWnPk08j3377LaZOnYoePXpg8uTJDtO3TbUPaPokg/bYvt69e0v3Du7Z\nsycCAwORn5/vFH34rLYBztN/L730EgCgsrISdXV1cHV1tUvf2S34neUcf41Gg5iYGMTFxUm3GHy6\nbTqdTroDlclkgr+/v7Sun59fo7tTOaqWtMdkMqGoqKjRPXbbQ99+8skniIyMxKpVq6Q3Y7PZ3K7b\nV1RUhIKCAoSHhztdHz5p25OTRZyl/+rr6xEcHAxPT0/Mnz8fWq3WLn3HyzK30PHjx/HNN98gNTUV\nv/vd73Dr1q0WnaLaHr7H0Nr2tGR9JSQlJeHq1avIzs7G5cuX8dlnnwFouu720j6r1YqpU6ciLS0N\nXbp0cao+fLptL7/8slP1n4uLC7755hsUFRVh/fr1OHv2rF36zm7BHxYWhosXL0rTBQUFiIyMtNfu\nZdOnTx8AgL+/PyZNmoR9+/YhLCwMFosFQMNBlrCwMABAREQELly4IK178eJFaZkja2l7fH19cfv2\nbWn+hQsXHLpve/XqBY1Gg27dumHevHnYs2cPgPbbvkePHmHKlCmIj49HbGwsAOfpw6ba5mz9BzSc\nNDJx4kSYTCa79J3dgt8ZzvGvqqqS/q0sKytDdnY2xo8fj4iICBgMBlRXV8NgMEh/9PDwcGRnZ6O4\nuBhGoxEuLi5wc3NTsgnNYkt7dDodduzYgfLycuzZs8eh+7a0tBQAUFdXh23btmHixIkA2mf7hBBI\nTEzE4MGDsWjRImm+M/Ths9rmLP1XXl6Oe/fuAQC+//57fP3114iNjbVP38lxZLq5jEaj0Ol0YsCA\nAWLdunX23LUsrly5IoKDg0VwcLCIiYkRGzZsEEII8eDBAzFp0iTRr18/ERsbK6xWq7TO2rVrxYAB\nA4S/v7/Izc1VqvRnmjZtmujTp4944YUXhJeXlzAYDDa1p6CgQISGhgofHx/xzjvvKNGUJj1pX6dO\nnYSXl5fYsGGDiI+PF0OGDBHDhg0TixcvbnSWVntrX15entBoNCI4OFiEhISIkJAQcejQIafow6ba\ndvDgQafpv/Pnz4vQ0FARFBQkxo0bJzZt2iSEsC1PWto+u9+IhYiIlMWDu0REKsPgJyJSGQY/EZHK\nMPiJiFSGwU9EpDIMfiIilfk/Pinn2/9h7OUAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0xb9b0650>"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Comparison of 'elasticaly obtained' results with true velocity and 'dispersed' velocity"
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "attenuation model"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "I modelled the viscoelastic data using a single SLS mechanism, which means that Q is frequency dependent and that for a given Q and a given w0 (intrinsic or set at w0 = 2*Pi*fdom ?), te whole Q 'spectrum', i.e Q(w), is defined. It is also equivalently defined by TauSigma and TauEpsilon the caracteristic relaxation times for the attenuation mechanism, here the assumption is made that stress relaxation is the same (TauSigma) for P- and S-wave, but since the study focuses on P-waves, who cares ? "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# getting tauSigma and tauEpsilon\n",
      "tauSigma = (np.sqrt(true_QP**2+1) - 1)/(w0*true_QP)\n",
      "tauEpsilon = 1/(w0**2*tauSigma)\n",
      "# computing dispersed velocity in the sense of Gillian\n",
      "dispersed_VP = true_VP/np.sqrt((1+w0**2*tauSigma*tauEpsilon)/(1+w0**2*tauSigma**2))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# load inversion results\n",
      "inverted_VP_frame = pd.read_table('winv.mp50', header=None, sep='\\s+')\n",
      "inverted_VP = inverted_VP_frame.values[:,2].reshape((301,1501), order='F')[:,nx/2]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# plot inversion results versus dispersion vel. and true vel.\n",
      "plt.plot(depth, true_VP, label='true P-wave vel.')\n",
      "plt.plot(depth, inverted_VP, label='inverted P-wave vel.')\n",
      "plt.plot(depth, dispersed_VP, label='vel. corrected for dispersion')\n",
      "plt.legend(loc='best')\n",
      "plt.title('Comparison of inversion results versus true model and the model corrected for dispersion')\n",
      "plt.xlabel('depth (m)')\n",
      "plt.ylabel('velocity (m/s)')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 21,
       "text": [
        "<matplotlib.text.Text at 0xb9aa350>"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAEXCAYAAACOIh8aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XlYVdX6wPHvOQwOiIIDgjI5gqgIKqA5IWYOt9Sw6+y1\nNCusLLXhZv4U773NlkPlkDmVU5qZZbc0NZzF2RxQQVEUQUVQRhnOWb8/vO44cAAt9DC8n+fxkbPH\nd+2z99rvXnvtfXRKKYUQQgghRBmmt3QAQgghhBAlkYRFCCGEEGWeJCxCCCGEKPMkYRFCCCFEmScJ\nixBCCCHKPElYhBBCCFHmVeiEZcWKFfTq1cvSYWgMBgMvvfQSDRo0YPDgwYXGv/fee4wdO9YCkRUt\nLCyM//znP5YO409bunQpXbp0sXQY4i+6n+/x6aef5v/+7/8ecER36PV6zp8//8CmtxRPT0+2bt1a\n4nQXLlxAr9djNBofQlRlQ0n74rZt2+jatSv29vYcO3as1NcfHh7OyJEjAYiLi8Pe3p6y9HYSe3t7\nLly48ECWfU8Jy88//0y3bt2oV68ezs7O9O3bl927dz+QgErT8OHD2bRpk6XD0Ozdu5ft27cTHR3N\nN998U2j8W2+9xcKFCy0QWdHmzZvHlClTLB1GqSkvJ4y/KiIiAjc3N0uHYRE6nQ6dTmfpMAgODmbR\nokWWDuNPKSvb8EF40N/Lu+++y+jRo0lLS6NNmzalvvz834u7uztpaWll6rtKS0vD09PzgSy7xIRl\n4cKFDB8+nA4dOrBv3z6io6MZPXq02RNuWWIwGCwdQiG7du3C19cXOzs7S4eiycvLs3QIwMP9vix9\nNVJW9s2yEseDYOnvGChTJ5HyzFwd9Vf23Qf5vSil2L17N506dfpT85eHlipL1hvFJixpaWm8+eab\nvP7663zwwQc0adIEe3t7nnrqKebMmQPc2ZmWL19Ohw4d6NixIytWrNB2sIiICFxdXZk/fz6NGjXC\nx8eHbdu2sWPHDgICAmjRogUrV67U1hceHs6QIUMYO3Yszs7OPP/881y6dEkb//7779O0aVPq1KnD\n8OHD2blzpzZu6dKldO7cmalTp+Lh4UF4eHihprv3338fPz8/atWqha+vLydPngQgMzOTefPm0bp1\nax577DF+/PHHQsv917/+RcOGDenduzd79+4tcpvFx8czZcoUPD09eeaZZzhy5AgA06ZNY9q0aaxd\nuxZ7e3uWLFlSaN78TX13m1q/++47WrRoga+vL8uXLwfgypUrVK9enZSUFG3eI0eOUK9ePW1n2rVr\nF8OHD6dRo0ZMnz6dpKQkbVq9Xs9XX32Fv78/3t7eZrfNqVOngMLN6xEREQwcOJBmzZrx4YcfmsSg\n1+tZsWIFfn5+NG3alJkzZxa5ncLDwxk6dChhYWG4uLiwbNkysrKyWLRoEYGBgXTu3Jm1a9dqJ56E\nhASGDx9OgwYNqFevHkOGDDHZTvkP9KKuoLp27QpAmzZtsLe3Z+3atWRlZfHss8/i6elJnTp16Nq1\nq9mTXVhYGK+//rrJsP79+2tlTElJYdasWbRs2ZI+ffqwefPmYst66tQpQkNDcXJywtnZmUmTJmnb\nt2DLiKenJ9u2bQMocr78MjIy6NOnD1euXMHe3p6aNWuSkJBQKI6lS5ea/X7zr7+4chX09NNPM3Hi\nRAYNGkTdunV5/PHHyczMZOrUqbi7uxMaGsrZs2e16Ys6VgDS09N59913cXV1pUePHiQmJpqsKyEh\ngX//+980bdqUwYMHExkZWWRc+Z07d46QkBDq1q2Lr68vH3zwAenp6Sbbev78+XTs2BF3d3fCw8PJ\nzc3Vxv/3v/8lMDAQb29v1q5dW+R63n77bXbu3MlLL72Evb0948eP18bt3bu3yGOkuOM2v7v7/bff\nfou3tzeNGjVi1apVREVF0blzZxo1asTs2bNN5vnhhx/o2bMnrVu3Zv78+WRmZprE9Nhjj+Hp6cmn\nn35aaH0///wz/fr1w8vLi5kzZ5pss+IkJSXxySef0Lp1a+rWrcvLL79cYjx3y7Z27VpatWrFo48+\nyrJly+jcuTPTpk3Dw8OD6dOnk5eXx5o1awgJCcHPz49FixaRk5OjLf/QoUOEhYXh5ORE06ZN2bRp\nU5HfS3H7U0n74l3Z2dnY29uTnZ2Nv78/zZo1A4rfz59++mkmTJjAoEGDqFOnDhEREYWWe/36dd54\n4w2cnZ0ZOHCgybYvWP9t2LCB7t274+DgQOPGjbXz693z2JQpU7QuCVFRUdpyiqt7786bf9snJiaa\nrY/BtBW7NM+tAKhi7NmzR+l0OnXixIkip1m8eLHy9fVVBw4cUIcOHVJ+fn5qyZIlSimlfvvtN2Vj\nY6NefvlllZSUpP79738rZ2dnFRoaqmJiYtS2bduUnZ2dysnJUUopNW3aNGVjY6NmzJihrl27pl55\n5RXVoUMHbV1r165VCQkJKjMzU33yySfK1dVVG7dkyRJlY2Oj3nrrLXXz5k2VlZWllixZojp37qyU\nUurEiRPKy8tLxcfHK6WUOn36tEpISFBKKTV16lTVvXt3debMGbV161bl6empfvvtN225tra2avr0\n6So5OVlNmzZNW6Y5Xbt2VS+99JK6du2aWrRokapZs6bKyspSSikVHh6uRo4cWeS84eHhasSIEUop\npWJjY5VOp1ODBw9WcXFxatOmTapKlSraskJCQtTChQu1eV977TUVFhamlFLq2LFjytXVVf36668q\nOTlZvfzyy2rYsGHatDqdTnXp0kUdPXpU3b59u9ht8/TTT6v/+7//U0opdf78eeXg4KDWrl2rLl++\nrAYPHqxGjRplstyQkBAVFRWlDh06pOzt7VVMTIzZst79rufMmaOysrJUVlaWmjBhghoyZIiKjY1V\nR48eVa1atVKbN2/Wyvf666+rzMxMlZ2drXbv3m2ynQwGg7bs4OBgtWjRIu37y/996XQ6de7cOe3z\nZ599poYNG6Zu3bql8vLy1K5du8zGu2PHDuXm5qZ9Tk5OVtWqVdO205NPPqnGjx+vEhMT1Y4dO1SD\nBg1UdHR0kWV96qmn1Jw5c1ROTo7KyMhQ+/btU0rdOWby79dKKeXp6am2bt2qlFJFzldQREREoeWY\niyP/92tu/cWVq6BRo0apWrVqqQ0bNqgrV66ooKAg5ePjo9577z2VnJysnn/+eTV69Ght+uKOlUmT\nJqlevXqp8+fPqw0bNqj69eurLl26aPP6+/urd999V6WkpKiNGzcqR0dHlZ6erpS6s89OmTLFbIwx\nMTFqy5YtKicnRx07dky1bdvW5Djy9PRUbdq0Ufv371dnz55Vnp6easuWLUoppY4fP67q1Kmjfvjh\nB3Xu3DnVt2/fQvtTfvn3w7uKO0ZKOm7zu7vfDxs2TMXHx6slS5YoOzs79be//U0dPXpUHTt2TNWs\nWVPFxcUppZTatm2bcnd3V7/++qs6e/as6tGjh5o2bZpSSqmUlBRVvXp1tWzZMhUfH6+GDx+ubGxs\ntH1uw4YNytfXV+3du1dduXJFDRo0SE2ePNkkjvzHX379+vVTI0eOVNHR0So7O1s7voqL5+4yBwwY\noM6dO6fV5QXr99mzZ6uQkBB14sQJFRMTo4KDg9UXX3yhlFLq2rVrqkaNGmrBggUqIyNDxcfHq9On\nTxf5vRS3P5W0LxZUcJ8obj8fNWqUsrOzUytXrlS5ubnq9u3bhZY3cOBANWLECJWQkKCWLl2q7Ozs\ntPNI/u2fk5OjPDw8tDohMTFRnTx5Uin1x/lx0qRJ6tq1a+r9999XDRo00NZRXN1rbtsXVR8XLH9p\nnluVUqrYhOWbb75RNWvWLHYB/fr1UwsWLNA+L1y4UPXr108pdafys7KyUklJSUoppS5fvqx0Op36\n4YcftOmbNWumIiIilFJ3KlR3d3dtXHp6uqpataq6du1aofUajUbl5uamDh48qBW+WrVqKjs7W5sm\n/8nq2LFjysPDQ0VERBQ6uNq0aaM2bdqkfX777bfV+PHjtWU4Ojpq81y5ckXZ2NhoO3N+169fV9Wq\nVTMZ16lTJ/Xdd99p5bubkJiTf/zdHfHQoUPaeC8vL/Xzzz8rpZT68ssvVUhIiMm22Llzp1JKqcmT\nJ6t33nlHmy8pKUnVrVtX5eXlKaXu7FBfffWVNr64bZP/hPbJJ5+ooUOHauOio6NVnTp1tHl0Op1a\nt26dNr5Xr15q3rx5RZa1cePG2mej0ag8PT21ClYppWbOnKnGjRunlFJq4sSJasSIEerChQsmy/mr\nCcucOXPUY489pk6dOmU2zvzxubu7qx07diillPriiy9Ujx49lFJKpaamKhcXF5WZmalN/8orr6gP\nP/zQbFmVUio0NFS99tprKjEx0WR4SQlLUfMVZG455uIoeHLPP5+5cr366qtauQoaNWqUduwrpdQ7\n77yj6tWrp33evXu38vDwUEqVfKz4+PhoZVZKqZEjR2rf49mzZ5WXl5fJugcMGKDWrFljtkzFWbhw\noXr88ce1z56enurjjz/WPj///PPqzTffVEop9cEHH5hccGzdurXEhOXLL780GVbcMVLUcZubm1to\n2Xf3+8OHDyullMrNzVXVq1dXc+bM0abp2bOndvE4fvx49dZbb2njfv31V+Xr66uUulPP5z8Bnzt3\nTul0Om37Dxs2TK1YsUIbf+TIEeXj42MSh7mE5ebNm6p69epa/Z9fcfHcXebdY00p8/V7p06dTE6U\n69evV3379lVKKTV37lzVv3//QutVqvD3UtL+VNy+aE7+faKk/XzUqFFaPW5Obm6uqlWrlsk+1qVL\nF7MJS25urnJxcVHLly9XGRkZJstZsmSJyQWvUko1aNBAHTp0qMS619y2L6o+Llj+0jq33lXsLSE3\nNzfS0tK0Wyfm7Nmzh3bt2mmf27VrZ3KrxsXFhTp16gBQv359AJOOSPXr1yc+Pl777Ovrq/1tZ2dH\nkyZN2L9/P3CnCTE0NJQGDRpQu3ZtEhIS+P3337Xp27Rpg62trdk4fX19eeedd/jnP/9Jw4YNmTp1\nKpmZmaSlpfH7778XW4aWLVui1+u18uTl5XH16tVC69i3bx+NGzc26aPSvn17du3aZTame+Hn56f9\n7eLiwpUrVwAIDQ1l7969JCYmsmPHDvR6PZ07dwZgy5YtvPfeezg6OuLo6EjTpk3JzMzk8OHD2rKC\ngoJK3DYFFfyumzZtSl5ensn+UTDe/N9tQfljOH36NHFxcfj6+mpxT5s2TevcPXnyZFxdXenYsSOP\nPPII33//fckb7x6MGTOG4OBgHn/8cVq3bl1kZzydTseQIUNYtWoVACtXrmT48OHAnWb869ev06BB\nAy32xYsXm3zv+csKMHPmTDIzM2nVqhW9e/dm+/bt9xTvn52vqDiKY65cixYtKnJ/1ul0Jse2k5MT\nLVu2NPl8d38o7lhJS0sjKirKZF/y9/fX/t6yZQuxsbFaTI6OjmzdutXkmC1Keno6r7zyCgEBAdSq\nVYsJEyaY1CFgug87OztrMe/fv7/ImIpirr9EUcd0Ucdt/lsIBd3d3tbW1tSuXbtQ3Xp32ebq6ePH\nj5OamkpkZKTJfI0bN6ZWrVra5y1bthAWFqbF1b17dy5cuMC1a9eKLfvu3bvx8PDQ6v/8ioonLS1N\nG1ZwX81fv2dkZLBnzx7+9re/aXE9/fTT7NmzB7hza7O4fiT5v5fi9qeS9sWSlHRO0Ol0xR6TUVFR\nGI1GGjdurA1r27at2dvW1tbWrFu3jm+//RZXV1fGjBlDbGysNr5Zs2ZUrVrVpBx79+7lzJkzXLx4\nsci6FwqfW++lPi7Nc+tdxSYsrVq1wtHRkR9++KHIaTp16sTBgwe1zwcPHtT6CvwZ+R8DS09P59y5\ncwQFBZGRkcHYsWMZNWoUp0+fJjk5mYYNG5p8cdbW1sUue/jw4ezdu5d9+/axefNmlixZgr29Pb6+\nvqVShg4dOnD+/HkyMjK0YQcOHLjnxzHvpzOYo6Mjjz32GN988w0rV65k6NCh2riQkBCmTJlCSkqK\n9i8jI4OAgABtmoLbyty2Kajgdx0dHY2VlZXJSeley6TT6bCystI+e3l54erqyqlTp7SYb926xdGj\nRwGoU6cO7733HleuXGHq1KkMHz6clJQU6tWrh42NjXZfOS8vj+PHjxcZj16vN9lnqlevzltvvcW5\nc+dYvHgxEydO1PrvFDR06FC+/fZbLl68yP79+xk4cCAAHTt2pF69ely9elWLPTU1lQ0bNpgtK9zp\n3f/555+TmJjIoEGDGDp0KEajkYYNG5KcnKz1RUpKSuLy5cslzleQlZVVoUrNXBwNGzY0qSDynxxL\nKpc55ipSc4o7Vuzt7fH29jaJ5fDhw9q+FBISQpMmTUz279TUVK1fXXE+//xzzpw5w5o1a7h58yYz\nZ868546OgYGBhWIqjpWV1X11oryX4/bPMldP+/r6UrNmTYKCgrTjDO7087l165ZJXAsXLiwUl5OT\nU7HrfOSRR7h48SI3bty4p3hat26Nvb29NqxgHZX/s52dHUFBQWzatEmL6ebNm1qfuu7duxeZWBf8\nXorbn0raF0tyL+eEgsdkft7e3uj1es6dO6cNO3ToUJHr79ixI+vXr+fChQvY2NjwxhtvaOOio6PJ\nysrSPh85coSOHTvSvHlz3Nzciqx7ofB3UVR9nF9pnlvvKjZhsbe354MPPuDjjz9m8uTJnDt3jtTU\nVNavX88rr7wC3Ol4OH/+fA4dOsSRI0eYP38+AwYM+NMBJSYmMnPmTK5fv87UqVPx9/enbt26pKWl\nkZ6ejouLC0ajUdtY9+rgwYNERkaSm5tLtWrVsLa21g6O/v3789FHH3H27FkiIiJYtWrVnypD3bp1\nCQgIYPLkyVy7do2lS5dy8uTJe34XzL1U9vmnGTZsGMuWLWPdunUMGzZMGz5y5EgWLFjA5s2bycnJ\n4datW8V2ECxu26g7tw0B6NevH5s2beK7774jPj6eadOm8cQTT2gZsrlYiypTweF6vZ7Bgwfz5ptv\nalcV586dY8eOHQCsXbuWy5cvYzQasbOzw87ODisrK+zs7OjQoQNffPEFycnJvPfeeyZXaQW1a9fO\n5AD66aefiImJ0ZZra2trchWSn5+fH3Xr1uXZZ5+ld+/e1KxZEwAHBwc6d+7M5MmTuXjxIgaDgRMn\nTmjrMbcNli9fzvXr11FKYWdnR40aNYA7V0F169ZlyZIlXL9+nWnTpplUTkXNV1CbNm1ISkoiISGh\nyG0O0KNHD3799Veio6M5ePAgy5Yt08aVVK6C7jVZgZKPlb59+zJjxgxiY2PZuHEjW7du1Zbv5eVF\njRo1mDFjBomJieTm5nLgwAFOnz5dYhxXrlzB0dERJycnDhw4wGeffXbPMffp04dffvmFn376ifPn\nzzNr1qxip2/Xrh1HjhwpNp78x8j9Hrf34u6y+/fvz6pVq9i2bRsxMTF89NFHWh332GOPcfjwYVas\nWMGVK1eYPn26yQlq5MiRfPjhh+zatQuDwcD169eLvYi9y8HBgZ49ezJx4kRiYmK4ffu21gJiLp4n\nn3zyvso2cuRIpk6dyuHDhzEajcTHx2udwgcOHEhERASLFi0iIyOD+Ph4zpw5AxT+Xkran8zti/eq\npP28pGPGxsaGRx99VOvounz5cpNEIr9r166xYcMGMjIysLKyomrVqiYJoNFoZNq0aVy/fp2PPvoI\nuNNaU1Lda05R9XFBpXVuvavEx5qfffZZli9fzu7duwkKCqJ58+YsXrxYu6IfMWIEEyZMYNy4cYSF\nhfHqq69qTeVQ+Aq7uMxUp9MxcOBATp06RatWrUhPT2f16tXAnabZ9957j5EjR9KmTRtycnK0WyB3\n5zW3rrvDUlNTee6556hduzbdu3cnMDCQESNGAPDGG28wYMAAQkNDeeedd/jkk0/o1q1bscstyooV\nK6hevToBAQFERESwdetWqlWrVuSyioq3qPXkH9avXz9iYmJwcXGhdevW2nAfHx+WLVvGmjVrcHV1\npXXr1ibvoym43OK2Tf6YGjduzNq1a/n666/p1q0bvr6+fPLJJ0Uut7jymhsXHh5O9+7dCQsLo3bt\n2vz973/XWk4OHjxIhw4dcHR0JDw8nHnz5mkJw/vvv8/evXtp3bo1RqPRpCm44Hpee+01ZsyYgaOj\nI2vXriU6OpqePXtSq1Ytxo4dy3/+8x+T5teChg0bxrZt20wSRID58+fj4eHBU089Rb169XjuuedI\nTU0tsqybNm2iVatW1K9fn+XLl/PFF19oid+8efNYvHgxgYGB+Pr64urqek/z5VezZk3eeOMNunbt\nqt0+NRdH586dGTFiBD169OCVV17hxRdfNJmmuHIVZG7/Le7YKe5YmTZtGp06daJz587MmjWL8ePH\nm8z7/fffk5ubS48ePXBxceGtt97SnhApbr+bMGECWVlZeHh4MGnSJMaNG3fPx2Tr1q1ZsmQJ06dP\np2/fvowaNarYeUeMGEFMTAz16tXj1VdfLXH5JR235uYtyd1pgoODmTlzJu+++y4DBgygf//+2lNv\nDg4ObNq0iSVLltCxY0cCAwNN9rk+ffrwr3/9i88++4x69erRsWNH7TZ9SXEsWrSIVq1a8fjjj+Pm\n5saaNWtKjMfcMs19p2PHjmX06NFMnTqV2rVr07NnT+0ptHr16rF161b27duHh4cH3bt3Jy4uDjD/\nvRS3P5nbF+9lm9/1V84JAHPnzsXJyQk/Pz/Wr19PWFiY2fUZjUZmzpxJw4YN8fb2Jjk5menTp2vT\nBQUFYWNjQ5s2bThw4EChJxmLqnvNxVhcfZx/2tI8twLo1P1cFj1g06dPJyYmhq+//trSoQghhBAV\nwtKlS1m0aNE99fMqy8rUq/nLUO70lxgMBvz9/XniiSeAO52P+vfvj7u7OwMGDDB5jn7OnDk0a9YM\nHx8fk3uuUVFRtG3blsaNG/P2228/9DIIIYQQZUmZSlgqyuugZ8+ejY+Pj1aWefPm4e7uTnR0tPYi\nPbhzz3Hu3Lls3bqVefPmmTQ1Tpo0iTfffJMDBw6wffv2IvsNCCGEEMWpKOfWMpWwTJs2ja+++srS\nYfwlly9f5r///S/PPvus1mK0f/9+xowZQ5UqVRg9erT2FsXIyEh69+6Nu7s73bp1Qymltb6cOXOG\nwYMHU6dOHUJDQ+/5TZ5CCCFEfqNGjSq2E215UaYSlopgwoQJfPTRRyYdIQ8cOKC9At/b21vrsBYZ\nGUmLFi206by8vIiMjCQmJsbkkUEfHx/27dv3kEoghBBClD3Fv7hE3JeNGzfi5OSEv7+/yW9C3E/f\nHHPNdkXNXxGa+IQQwhIqSp/JykRaWErRnj17+OGHH2jUqBFDhw5l27ZtjBw5koCAAO2HpqKiorQX\nQQUFBZm8pOz06dMEBATQtGlTk5d5nTp1ig4dOphd5933OFTEf9OmTbN4DFI2KZ+Ur+L9E+WTJCyl\n6N133+XSpUvExsayevVqQkJC+PrrrwkKCmLx4sVkZWWxePFiLfkIDAxk06ZNxMXFERERgV6v1170\n4+3tzerVq0lKSmL9+vX39Up1IYQQoqKRhOUBunvLJiwsjLi4OLy8vIiPj+eFF14A7vzWR1hYGCEh\nIYwbN87k5+BnzJjBhx9+SEBAAF26dKF9+/YWKYMQQghRFpSpF8eJ+6PT6Sp082ZERATBwcGWDuOB\nqMhlAylfeVfRy1fR686KShKWckwOOiGEuH9Sd5ZP8pSQEEJUILduwZdfwv9+8LvSSjAe54zhZ0uH\nIUqRJCxCCFGB7N0Ln38OTz1l6UgejDxuc1V/iFRdHLbUxMUQSHXqkUsWN/SnuK7/nav6w5y2Xo1P\n3nD0ytbSIYtSIgmLEEJUIHmxl1iW+wZdLlSsJpaM3EzO3DhD3K047GzsqGFrR44hh+SsZKz11uQY\ncqhha49D1Vo4VHWgoX0ANWyvmF2WvMGqfJKERQghKhDbM7/jlnkGnvqnpUP5y5RSRCVFsfX8Vo4k\nHqFHqx70ajKButXratMYlZGkzBvUruaItf4eT2lr1z6giMWDJAmLEEJUIIab6STXbornoEGWDuWe\nGYwGLty8QOzNWOJT47mSdoVzKefYFruNajbVGPPUGD7x/xGHqg6F5tUDToUXWbzBg0sjbPGQScIi\nhBAViPFWGsbq9pYOo1hZuVnsjNvJlvNb2Ba7jVPXT+Fk50Qjx0Y0tG9Iw5oNaefSjpcCX6JN/Tby\nMyQCkIRFCCEqFJWahqpRw9JhmDAYDRxOOMyW81vYEruFyMuR+Dn70bNxT2b1noW/sz92tnaWDlOU\ncZKwCCFEBaLS08HBsi0sd/uebL+wna2xW/ntwm+41HDh0caPMqHDBLp6dKVmlZoWjVGUP5KwCCFE\nBaLPSEPnUeehrtOojBy/epztF7ez/eJ2dlzcgb2tPd08u9Hfqz9z+syhgX2DhxqTqHjkTbflmLyt\nUQhR0C+Nw8jpbMUvg41kG7Lxd/ank1snWtdvfe9P0dyDq+lX2XxuM5vObeLX879Sq0otgj2D6ebR\nja4eXXGr5VZq6yptUneWT5KwlGNy0Akh8svOy+bHll7s9s/D7ZWJ2NnYcTDhILvjdnM59TLtG7TH\n08GTKtZVAMjIySAjN4PsvGycazjjXsudlvVa4ufsRyPHRuh1f/w+7q3btziUcIhfz//KpphNxN6M\nJaRRCL2a9KJXk154OHhYqtj3TerO8kkSlnJMDjohRH6jN4xm4Ks/4zF2Fq0mmz66m5yVTOTlSK6k\nXSHHkINCYWdjh52tHbZWtiSmJ3Lh5gVOXj/J0cSjpGSl0KR2E6paV+VaxjWupl+ljXMbQhqF0LtJ\nb4Jcg0q1xeZhkrqzfCqfe5sQQggTa06uYfel3YxN8sKuXuE+LLWr1aZPsz73vLzkrGRiU2LJNmRT\nu1ptmtVuhpXeqjRDFuK+SMIihBDlXFZuFq9tfo1VA1dRbdIEqtT7608J1a5Wm9rVapdCdEKUDn3J\nkwghhCjLPt3/KQENA+jk3omqeelUdyrbL44T4s+QFhYhhCjHsvOymblvJr+O/BWlwM6YRnWnsvXi\nOCFKg7QEJzPoAAAgAElEQVSwCCFEObbqxCp86/vSyqkV2dlgTxq2daSFRVQ8krAIIUQ5pZTik72f\nMKnjJADS0xQ1SAd7SVhExSMJixBClFNbY7diVEZ6Nu4JQHrSbfJ0NmAtd/tFxSMJixBClFOf7P2E\niR0nar9mnHE1jUy99F8RFZMkLEIIUQ6dun6KwwmHGdZ6mDbs9rU0sqzldpComCRhEUKIcmjWvlmM\nCxhHVeuq2rDbSelk20jCIiomudEphCi3cnPh4kVLR/HwJWbGs+bEt2x64jQxMfmGR6dRu4okLKJi\nkoRFCFFuzZ0L4eFQp/Cb6Cu064HvoM97luGLnEyGd81Mw6+29GERFZMkLEKIciv9Siqf9D/CM89Y\nOpKH53xKLK/+soKvnvwahxe3m47cfgCOSwuLqJgkYRFClFut9i6k+/HZENvI0qGUmlxjHrfzbmOt\nt6aqdVV0+cYZlJHkxCMssG+IQ9TH5hcwdOhDiVOIh00SFiFEuWWVdpOzXZ6l/Q9TLR3KX5JjyGHD\n6Q0sOrKIPZf24OngSXLWVfKMefRs0pNeTXrhWNWRaRHT8K0/iEX9FoFOV/KChahA5CmhUnT79m2C\ngoLw8/OjQ4cOzJw5E4Dw8HBcXV3x9/fH39+fn3/+WZtnzpw5NGvWDB8fH3bt2qUNj4qKom3btjRu\n3Ji33377oZdFiPJAn5mGvlbZ7rNhMBq4mn6VaxnXSM5KJis3C6UUGTkZ7I/fz5tb3sRtphufH/ic\nEb4juPraVX4P+53LEy+zZ8weOrt15ruo75i5bybPtn2WRf0Wae9dEaIy0SmllKWDqEgyMzOpXr06\n2dnZtGvXjvXr17Ny5Urs7e2ZOHGiybTXrl2ja9eubN68mdjYWCZMmMDhw4cB6Nu3L6NGjeLRRx+l\nf//+zJo1i/bt25vMr9PpkK9PVGab3cZQf0AHGrw/gM3nNtO0dlOCXIMsGlOOIYeICxFsOLOBXXG7\niL4RjZ2tHTp05BpzuZ13m+y8bKpYV6F5neb0atKLZ9s+S/M6zS0ad2UidWf5JLeESln16tUBSE9P\nJy8vjypVqgCYPTgiIyPp3bs37u7uuLu7o5QiPT2dGjVqcObMGQYPHgxAaGgokZGRhRIWISo769vp\nXLaNo+fcljzi9gjHrh6jo2tHFj6xEDtbuwe23lxDLhm5GWTkZJCRm0FyVjKRlyPZGbeTrbFb8a7r\nzQCvASzqtwjvut7UsDVtBTIqIzp00lIixH2QhKWUGY1G/P39OXnyJLNmzcLd3R2ATz/9lLVr1/Lk\nk08ybtw47O3t2b9/Py1atNDm9fLyIjIyEg8PD5yc/nhc0cfHhxUrVvDiiy8+9PIIUZZZ511j+Y3N\nrH5pHSGNQsjMzWTcT+MI+SqEn4b9RN3qdf/S8g1GA/su72Nj9EaOJBzh5PWTXE2/ilEZsbO1w87G\nDjtbO2pVqUW7Bu3o59WPz/p+hnMN52KXq9fJ3Xgh7pckLKVMr9dz7NgxLly4QN++fenUqRNhYWFM\nnTqV1NRUXn/9dRYsWMBrr71mttXF3BVXcU2X4eHh2t/BwcEEBweXRjGEKPOUUlhVPUbnpn8npFEI\nANVtqrOk/xImb5tM58WdWT94PS3qtTCZ53TSaX678Bv7Lu/DoAzUt6uPW0033Gq54VbTjarWVfn9\n6u9sOreJzec242LvQj+vfrwU+BIt67WkgX0DbK1spXWkHImIiCAiIsLSYYi/SPqwPECvvfYaTZs2\n5YUXXtCGHTt2jHHjxrF7925+/PFHtmzZwuzZswHw8/Nj586d2Nvb07hxY86fPw/Axx9/TNWqVQu1\nsMh9WFGZ/XT2Jxp2Gkj9r7bj0qdwv5WFhxby1ta36N20N40dG3P2xlkiLkRQzaYa3T2708mtE7ZW\ntiSmJ3Ip9dKdf7cukW3IxquOF481eYxeTXrh4eBhgdKJB0nqzvJJWlhKUVJSEtbW1jg4OHDjxg02\nb97MpEmTSEhIwMXFhby8PFauXEnfvn0BCAwM5PXXXycuLo7z58+j1+uxt7/z0idvb29Wr17No48+\nyvr165k1a5YliyZEmZJjyGHCpgn8nFqbmi61zU4ztt1YnvB6gg2nN5CYnkjvpr15/9H38XTwfLjB\nCiFKhSQspSghIYFRo0ZhMBhwdnbmtddew8XFhX/84x8cPXoUW1tbunbtSlhYGAD169cnLCyMkJAQ\nbG1tWbBggbasGTNmMGLECN566y2GDBkiHW6FyGf2vtk0rd2MGjkHqeZU9JtdnWs483z75x9iZEKI\nB0VuCZVj0qwpKqMLNy/Q/ov2bB60F+9Gbaiedg1qlO13sYiyRerO8km6qgshyg2jMvL8xueZ2HEi\nTjSiCtlg9+AeXxZClB2SsAghyo0Pd39Iek46rz/yOpnXM8jS28kr6oWoJKQPixCiXPjxzI/MiZxD\n5LOR2FjZkHk1jUwre+RmkBCVgyQsQogyb1fcLkb/MJqNQzfiVssNgNvX07htLemKEJWF3BISQpRp\ney/tJfSbUFaGrjT5naDspDSybYt+QkgIUbFIwiKEKLN2XtxJ/9X9+erJr+jZpKfJuJzkdHKqSMIi\nRGUhCYsQokxa/vtyBq4ZyIrQFfRu2rvQ+LyUNPKqScIiRGUhfViEEGXK7bzbTNk2hXVR6/ht1G+0\ndGppdjrjrTQM1aQPixCVhbSwCCHKhMT0RGbvm43P5z7E3oxl/7P7i0xWAAyp6RjtpIVFiMpCWliE\nKCc++gi+/trSUZQ+g3UqCW1eIbXBBuwT+1L7/NdE3+hEj/8UP1/ouTRcgyRhEaKykIRFiHJizx54\n5hno0cPSkZSe1JwURu/sRo/aHXml5Xlq2jrc87z1P0+jTj25JSREZSG/JVSOye9hVC5zm89khFpO\nzZqWjqR0KOB8ynlsrGxwq+nGfb+v9vJlmDwZXnnlAUQnKjKpO8snaWERopxokbSTlGFDqPlMiKVD\nKRXfnVrHt6e+ZemApeisq/y5hfj4lG5QQogySxIWIcqJKjlpKN820K6dpUP5y45fPc7zW79g16Rd\nVKnrbelwhBDlgCQsQpQTVXPTOJxzgNfXLqS9S3te6fAKVa2rFjvP9YzrrDi+grhbcXjX9eZvzf5G\nw5oN/9T6lVJczbjK9YzruNVyw6Hqvfc3yS8pM4kB3wxgdu/ZeEuyIoS4R9KHpRyT+7B/Xq4hl8T0\nRO13acqDU7bNeP7FJEb/4xN+PPsjOy7uoI1zG25k3qCqdVVG+I7gGb9nsLO1w2A0sPDwQqb+NpW/\nNf8bLeq24Pi14/x09ie863rzlM9TPOL2CF51vHCs5gjcSUiSs5KJT4snPjWe+LR4zqecJyY5hujk\naGKSY6hiVQUnOycupV6ieZ3mDG45mBG+I2hg36DY2POMedy8fZNdcbuYuGkiI9uMZHrw9Iex2YQo\nROrO8kkSlnJMDro/b07kHJYdW8ah5w5ZOpR7ohTE1ajOr8te59mn7pzoL926xMnrJ3GycyIlK4W5\nB+ey8+JOHnF7hN+v/k7Dmg35vO/n+Nb31ZaTY8hhW+w21kWt40jCEc7cOEOOIQdrvTV5xjzsbOxo\nYN9A+9fYsTFNazelWe1mNK3dVEtucg257L60m+W/L2dd1DqCGgYxqs0oejbpiV6nJ+p6FDvjdrIz\nbid7L+0lNTuVmlVq0sa5Da8EvcIA7wEW2Y5CgNSd5ZUkLOWYHHR/Xrsv2nEs8Ri3/nkLO1s7S4dT\nohupGejq10AfHYODa5Mip4tNieVQwiEaOzbG39kfna74Z2+UUmQbsjEYDVjprUq8xWROZm4m66PW\ns+zYMg5eOYhRGWlauyldPLrQxb0Lndw64WTnVGIsQjwsUneWT5KwlGNy0P05J66doPfy3rjYu/Dx\nYx/T1aOrpUMq0dJ96xj+yN+xycqCKn/yiRohBCB1Z3klr+YXlU7EhQgeb/44ndw6sffSXkuHc0/+\ne+o7FHpJVoQQlZYkLKLSSUxPpKF9Qzq4dmBf/D5Lh1MipRS/X9hKhpW8hl4IUXlJwiIqnasZV6lf\noz7tXNpxJOGIpcMp0fmU89jdVtzWV5BX3AohxJ8gCYuodK6mX6W+XX08HDy4knaFXEOupUMq1o6L\nO/C1aku2rbSwCCEqL0lYRKWTmJ5I/Rr1sbWyxbmGM5dSL1k6pGLtjNtJc0NLSViEEJWaJCyi0rma\ncaeFBaCRYyMu3LxQaJp9l/exP37/Q47MvJ1xO3G73ZTcqpKwCCEqL0lYRKWilLpzS6jG/xIWh0bE\npsSaTHM44TA9v+7JjD0zLBGiicT0RG5k3sAhrTZ5VWtYOhwhhLAYSVhEpZKanYqNlQ3VbaoD/0tY\nbpomLB/u/pBn/J4hMj7SEiGa2HlxJ53dO2O8lYmxurSwCCEqL0lYRKWS/3YQgKeDZ6GE5XDCYZ5r\n9xxp2Wkkpic+7BBN7IzbSRf3LqjUNIw1JGERQlRekrCISiX/7SAo3Ifl1u1bJKQn0KJuCwIbBhJ5\n+cG3suyK28WwdcP41/Z/kWPIMRm34+IOunh0QZeehq6G3BISQlRekrCUotu3bxMUFISfnx8dOnRg\n5syZAKSlpdG/f3/c3d0ZMGAA6enp2jxz5syhWbNm+Pj4sGvXLm14VFQUbdu2pXHjxrz99tsPvSwV\n1dWMqzjXcNY+F+zDcjTxKL71fbHSWxHkGsT+KyV3vM3KzeLSrT/3pNGv534l9JvQO2/dvbyXwd8O\n1l4ZfvHmRS6lXqKtS1vISEdXS1pYhBCVlyQspahq1ar89ttvHD16lO3bt7No0SKio6OZN28e7u7u\nREdH4+rqyvz58wG4du0ac+fOZevWrcybN4/x48dry5o0aRJvvvkmBw4cYPv27Rw8eNBSxapQEtMT\nTW4JNazZkPScdJKzkgE4lHDoToIABDQI4ED8gWKXt/X8Vpp/1hzf+b4MXTeU7Lzse45lz6U9DP9u\nOOsGrePFwBf5fvD3XLx5kXkH5wEw/9B8RrUZha2VLVaZaVhJwiKEqMSsLR1ARVO9+p3OnOnp6eTl\n5VGlShX279/PlClTqFKlCqNHj+a9994DIDIykt69e+Pu7o67uztKKdLT06lRowZnzpxh8ODBAISG\nhhIZGUn79u0tVq6HbfwXqzl25RQKIwojoLS///hXeJiRXPK4TS5Z5JFl8nc2aWRwjW5MI1zLQ/Q4\n4serHxyhMT34jkM04lHCIyEVf3ZzlGnhCh2Ff2k4lm18yxAGsho3HuH74/+gzfFB/J1vscKm2PJd\nIIK1DOJJvmbrki5sBaAKnfiGN37qzMafFDtYzDPsJHwvtL6ehpWDJCxCiMpLEpZSZjQa8ff35+TJ\nk8yaNQt3d3cOHDiAt7c3AN7e3uzff+c2Q2RkJC1atNDm9fLyIjIyEg8PD5ycnLThPj4+rFixghdf\nfLHQ+sLDw7W/g4ODCQ4OfjAFe8gWnP4XXjY9qFutPjr0gA4d+nz/Cn6+M40VNlhTDWuqYkM1rKmm\n/V8Fe6pQk2rUMVmXC21J4DCNCOEC2+lGOAD2NAAUaVyhJg216Y+ylBOs5honGMgqGtMDgFBW8A2h\nrGckoaxAj1WhchkxsI+Z7OZDnuIbGtHdZHwdmvF31rKTd+lGOHVoDoCPWxpuraQPixB/RkREBBER\nEZYOQ/xFkrCUMr1ez7Fjx7hw4QJ9+/alU6dO9/Uz5jpd4Sv54ubPn7BUJK2vZ7DAvyXtmjUseeIi\nGYD0//0rWtdYHQev/MCw1vbE/ZbNnH5R6HSnATD+5kq/5gsIaBgAQGR8JHEHF7DMfwxtG4ymmnUW\n8KO2rP8zPMO/tv+L2na9eDngJaz0fyQtx67+zpeHv6S9bQ2+7vAu9e3STeY1dTc5/d/4Hy+Ai7Sw\nCPFnFLyYmz59uuWCEX+aJCwPiKenJ3379iUyMpKAgACioqLw9/cnKiqKgIA7J7+goCC2bNmizXP6\n9GkCAgKwt7fn6tWr2vBTp07RoUOHh14GS1r2cyIuUUvAxankif+i9jlpEH+E3N03eeV2VXRXF2rj\nhl5PwzpyDdQ+gkEZuXV+CzMaBFD76kngZKFl2QLhxoYcTDjIrh/+gWtNtzsvq8u4SlZeFm/V88Gl\nhgu6cxvuL0g3N/Dx+WsFFUKIckwSllKUlJSEtbU1Dg4O3Lhxg82bNzNp0iRSU1NZvHgxH374IYsX\nL9aSj8DAQF5//XXi4uI4f/48er0ee/s7V9He3t6sXr2aRx99lPXr1zNr1ixLFu2hszIqYqb8h8AB\nPR/4uqob8wj9sC717G4zrdt7+PmO0MYdO76KtafW8t3g71hx7CtWnTDy8/Cfi12eNRCkFD+e/ZF1\n57dgpbeii3sXnmj+BDZWxfdtEUIIYZ5O3c/9ikokMTGRM2fOoNfrad68OfXr1y9xnuPHjzNq1CgM\nBgPOzs4MHz6cf/zjH6SlpTFixAiOHDlC27ZtWb58OTX+906N2bNn8+mnn2Jra8uCBQvo0qULcKdV\nZcSIEaSkpDBkyBCto25+Op3uvm43lSfRjjbcWr6Z9n/rXvLEpeB00mnWnFzD+KDxOFR10IYnpifi\n87kPJ8edZMA3A5jSZQpPeD3xUGISQjwYFbnurMgkYcknISGB+fPns3LlSqpVq0aTJk1QSnH+/Hky\nMzMZPnw4L7zwAi4uLpYOFajYB905B2vSvonAr1dnS4fC+J/Hs/HsRpzsnNg9erdJvxQhRPlTkevO\nikxuCeXzwgsvMGLECE6cOEGVKlVMxuXk5PD999/zwgsvsGHDffY/EPfNSoG1bdm4ffLPzv/kRtYN\nPuvzmSQrQghhIdLCUo5V5KuESzWtSP/vflp0bmfpUIQQFUxFrjsrMnnTrRlr1qwhNTUVgLlz5zJ2\n7FhiYmIsHFXloi9DLSxCCCEsTxIWM/79739Ts2ZNjh8/zldffUVISAivvvqqpcOqVKyUwspa7lgK\nIYS4QxIWM2xs7lzZL126lHHjxjF06FCuXLli4agqFysj2FSxtXQYQgghyghJWMxo06YNI0eOZOPG\njQwaNIjbt29jMBgsHValopcWFiGEEPlIp9si/Pbbb7Ro0QJnZ2cSEhI4fvw4jz32mKXDMlGRO47d\nrKoj98xl6nn8lVfzCyFEYRW57qzI5BI2n3bt2tG5c2f69OlDcHAwVatWBcDFxaXMvHulstArsLKS\n3VMIIcQd0sKST25uLrt27eKXX34hIiKC2rVr07t3b/r06UPz5s0tHV4hFfkqIcNWh7qSTI26jpYO\nRQhRwVTkurMik4SlGPHx8fzyyy9s2rSJmJgYOnTowNy5cy0dlqYiH3S3rXXoU9KxtbezdChCiAqm\nItedFZkkLPfIaDSyd+9eOnXqZOlQNBX5oMux1mGdfht91SolTyyEEPehItedFZl0EjDjxIkTLFiw\ngL1795KdnQ3c2cF///13C0dWOeQZjFgZQW8ju6cQQog7pIXFjE6dOvHcc8/RsWNHbG3/eBeIp6en\n5YIyo6JeJdzOyaNqFRswGkGns3Q4QogKpqLWnRWdXMIWYejQoSbJinh4crJzsQX0kqwIIYT4H2lh\nMWPPnj3MmjWL3r17U6tWLeBORh4aGmrhyExV1KuE5Bup2DvVwsZQ8comhLC8ilp3VnTSwmLGqlWr\nOHbsGDY2NiatLGUtYamocnJyMUrjihBCiHwkYTHjl19+4eTJk3JLyEJys3MwSMIihBAiH/ktITO6\nd+/O3r17LR1GpWXIlRYWIYQQpqSFxYxdu3bx5Zdf0rBhQxwcHAB5rPlhys3JwyAdboUQQuQjCYsZ\n//3vfy0dQqVmyM3FIG1/Qggh8pGEJZ+bN2/i4OBQ7PtWbt26pT05JB6M3Nw8uSUkhBDChCQs+Qwa\nNAgnJyeGDRtGixYtcHd3RynFxYsXOX36NCtXruT69ets3rzZ0qFWaIbcXLklJIQQwoS8h6WAI0eO\n8MUXX3D06FHOnTsHQJMmTfDz8+O5557D39/fwhH+oaK+S+BgxAEaPt4Bl3SDpUMRQlRAFbXurOgk\nYSnHKupBt3/LXho82QXXtDxLhyKEqIAqat1Z0UnXRlHm5ObILSEhhBCmJGERZY5BOt0KIYQoQBIW\nUeYY83Ix6CVjEUII8QdJWMyYOHEiJ0+etHQYlZYhJxej3BISQgiRjyQsZrRo0YLnnnuOwMBA5s+f\nz61bt+5pvkuXLtG9e3datmxJcHAwK1euBCA8PBxXV1f8/f3x9/fn559/1uaZM2cOzZo1w8fHh127\ndmnDo6KiaNu2LY0bN+btt98u3QKWcUaDQRIWIYQQJuQpoWKcPn2apUuXsnLlSjp37syLL75Ip06d\nipw+MTGRxMRE/Pz8SEpKIjAwkGPHjvHJJ59gb2/PxIkTTaa/du0aXbt2ZfPmzcTGxjJhwgQOHz4M\nQN++fRk1ahSPPvoo/fv3Z9asWbRv395k/ora0/3npd/i/tpIWiZlWToUIUQFVFHrzopOWliKYDAY\nOH36NFFRUdSrV482bdrw73//m3HjxhU5j7OzM35+fgDUrVuXli1bcuDAAQCzB0dkZCS9e/fG3d2d\nbt26oZQiPT0dgDNnzjB48GDq1KlDaGgokZGRD6CUZZNROt0KIYQoQN50a8aECRP48ccfCQkJ4e23\n3yYwMBCAN998Ex8fn3taRkxMDCdPniQoKIidO3fy6aefsnbtWp588knGjRuHvb09+/fvp0WLFto8\nXl5eREZG4uHhgZOTkzbcx8eHFStW8OKLLxZaT3h4uPZ3cHAwwcHBf67QZYgxL09uCQkhSk1ERAQR\nERGWDkP8RZKwmOHr68t//vMf7OzsCo3bs2dPifOnpaUxePBgZs6ciZ2dHWFhYUydOpXU1FRef/11\nFixYwGuvvWa21UVn5kRdXNNl/oSlopBX8wshSlPBi7np06dbLhjxp8ktITO+/vrrQslKjx49AHBw\ncCh23tzcXAYOHMjIkSPp378/AE5OTuh0OmrVqsWLL77I+vXrAQgKCuLUqVPavKdPnyYgIICmTZty\n9epVbfipU6fo0KFDqZStPLjT6VZ2TSGEEH+Qs0I+WVlZ3Lhxg6SkJJKTk7V/p0+fJi0trcT5lVKM\nGTOGVq1a8eqrr2rDExISAMjLy2PlypX07dsXgMDAQDZt2kRcXBwRERHo9Xrs7e0B8Pb2ZvXq1SQl\nJbF+/XqCgoIeQInLJmOePNYshBDClNwSymfBggXMnj2bK1eu0K5dO224h4eHSQJSlN27d7N8+XJ8\nfX21H0l89913WbVqFUePHsXW1pauXbsSFhYGQP369QkLCyMkJARbW1sWLFigLWvGjBmMGDGCt956\niyFDhhR6QqgiU3nyWLMQQghT8lizGZ9++ikvv/yypcMoUUV9NG/1f2bhPi+cR+JvWjoUIUQFVFHr\nzopOWljy2bZtGyEhITRo0IDvvvuu0PjQ0FALRFX5qLw8DHq5WymEEOIPkrDks337dkJCQvjxxx/N\nPq0jCcvDYczLQ8ktISGEEPlIwpLP3Ufdli5datlAKjljXp481iyEEMKEtLubMXnyZFJSUrTPKSkp\nTJkyxYIRVTKGPJQ81iyEECIfOSuY8dNPP+Ho6Kh9dnR0ZOPGjRaMqHKRHz8UQghRkCQsZjg5OXHl\nyhXtc3x8vEkCIx4sZZBOt0IIIUxJHxYzxowZQ9++fRkyZAhKKVavXs3kyZMtHValoaTTrRBCiALk\nMtaMIUOG8P3332NlZYW1tTXff/89gwcPtnRYlYfRKK/mF0IIYUJaWIrg6elJ165dAWjUqJGFo6lc\n7rSwSMIihBDiD5KwmBEREcHYsWNp3rw5ANHR0SxcuJBu3bpZOLLKQSkDRr3cEhJCCPEHSVjM+Oij\nj9i4cSNeXl4AnD17lldffVUSlodE5eVh1FlZOgwhhBBliLS7m5GSkoKzs7P2uX79+ty8Kb9r89AY\npNOtEEIIU9LCYsaoUaPo06cPTz31FEop1q9fz9NPP23psCoPoxGjPNYshBAiH0lYzHj++efp2LEj\nGzduRKfTMW/ePFq3bm3psCoPedOtEEKIAiRhKYKvry++vr6WDqNyMhikhUUIIYQJSVjyqVGjhtlf\naQbQ6XSkpqY+5IgqJ2U0yHtYhBBCmJCEJZ/09HRLhyAAjAaQFhYhhBD5yFmhCAkJCaxYsQKA69ev\nExsba+GIKhGjEaNeHmsWQgjxB0lYzPjiiy8YOnQo06dPByAnJ4cRI0ZYOKpKRDrdCiGEKEDOCmZ8\n/fXXbN68GTs7OwAaNmxIWlqahaOqPHRGo9wSEkIIYULOCmbUqlULfb4TZlxcHK6urhaMqJIxylNC\nQgghTMlZwYxRo0YxfPhwbt68yfTp03n88cd59tlnLR1WpaEzGlDyan4hhBD5yFNCZvz9738nICCA\ndevWYTQa+emnn3Bzc7N0WJWH0YiSFhYhhBD5SMJixscff8yQIUOYNGmSpUOplHRGA8pKWliEEEL8\nQS5jzUhLS+Oxxx6jc+fOfPbZZ1y9etXSIVUqOqMB5CkhIYQQ+chZwYzw8HBOnjzJ559/TkJCAl27\ndqVHjx6WDqvyMBpR8h4WIYQQ+UjCUgwnJyecnZ2pU6cO169ft3Q4lYbOaABJWIQQQuQjCYsZc+fO\nJTg4mB49epCUlMSXX37J77//bumwKg2d0SCdboUQQpiQs4IZly5dYtasWZw6dYrp06fj4+Nzz/N1\n796dli1bEhwczMqVK4E7fWL69++Pu7s7AwYMMPnNojlz5tCsWTN8fHzYtWuXNjwqKoq2bdvSuHFj\n3n777dItYBmnUwqk060QQoh8JGEx47333sPPz+++57OxsWHmzJmcPHmSb7/9lilTppCWlsa8efNw\nd3cnOjoaV1dX5s+fD8C1a9eYO3cuW7duZd68eYwfP15b1qRJk3jzzTc5cOAA27dv5+DBg6VWvrJO\nZzCAvIdFCCFEPpKwlCJnZ2ct0albty4tW7bkwIED7N+/nzFjxlClShVGjx5NZGQkAJGRkfTu3Rt3\nd74k34MAACAASURBVHe6deuGUkprfTlz5gyDBw+mTp06hIaGavNUBjpllMeahRBCmJD3sDwgMTEx\nnDx5ksDAQJ555hm8vb0B8Pb2Zv/+/cCdhKVFixbaPF5eXkRGRuLh4YGTk5M23MfHhxUrVvDiiy8W\nWk94eLj2d3BwMMHBwQ+mQA+RThnkt4SEEKUmIiKCiIgIS4ch/iJJWB6AtLQ0Bg8ezMyZM6lRowZK\nqXueV6fTFRpW3Pz5E5aKQm80gl52TSFE6Sh4MTd9+nTLBSP+NLmMLWW5ubkMHDiQkSNH0r9/fwAC\nAgKIiooC7nSmDQgIACAoKIhTp05p854+fZqAgACaNm1q8rK6U6dO0aFDh4dYCsvSGY3S6VYIIYQJ\nSVhKkVKKMWPG0KpVK1599VVteFBQEIsXLyYrK4vFixdryUdgYCCbNm0iLi6OiIgI9Ho99vb2wJ1b\nR6tXryYpKYn169cTFBRkkTJZgk4Z0cl7WIQQQuQjCUsp2r17N8uXL2fbtm34+/vj7+/PL7/8QlhY\nGHFxcXh5eREfH88LL7wAQP369QkLCyMkJIRx48Yxe/ZsbVkzZszgww8/JCAggC5dutC+fXtLFeuh\n0ylpYRFCCGFKp+6ng4UoU3Q63X31jykvvm3VlLQ23XhmxSJLhyKEqIAqat1Z0UkLiyhz9EYjOmlh\nEUIIkY8kLKLMufMeFnlKSAghxB/krCDKHL1S6CtxC0vt2rVJSUmxdBhClHuOjo4kJydbOgxRSiRh\nEWWOvpI/1pySkiL314UoBebeayXKL7klJMocnTKik1tCQggh8pGERZQ5enmsWQghRAGSsIgyR68U\nemtpYRFCCPEHSVhEmaNXRvTyplshhBD5SMIiyhy9UaGTFhZRiQUHB7Nokbw4UYj8JGERZY5eScJS\nVnl6erJt27aHtr78v7FVt25dBgwYwPfff//Q1m8pOp1OnnARogBJWESZo5cfPyyzSnqleV5eXqmv\ns2HDhqSlpREfH0+vXr0YPXo0t2/fLvX1CCHKNklYRJljpYzobWwsHYYoYOTIkcTFxfHEE09gb2/P\njBkzuHDhAnq9nrVr19KqVSt69uzJ9u3bcXNzM5nX09OTrf/f3p3HVVWt/wP/nENiQISIohjCETQG\nZTjiYVIEcQRL1FIGozCyK+rVnNOrgtbNn6Y5pmiBQ4HTTY0URUEPoCJg4JAICg50cQCKEFFien5/\n8HVfD0MioZxzfN6vl6+Xe1p7PXu55WGtvfdKSBCWjxw5glGjRsHCwgJr1qzBgwcPnnr+9u3b4+OP\nP0ZlZSVOnz7dYHtoaCimT58OAKiqqoKOjg7mzZsHAHj06BFeffVV/PHHHwCAcePGwcjICN27d8es\nWbNw/fp1AEBqaiqMjIwUkrIDBw7Azs7umetuZWWFw4cPC8vV1dXo3Lkzzp8/DwDIy8vDvHnzYGpq\nikmTJiErK+up14CxlxknLEzpiAg8l5AS+u6772BiYoJDhw6hrKwMc+bMEbZFR0cjJiYGR48ebbQH\n5skhjpiYGHz66adYuHAh5HI5zp49i+XLlz/1/I8ePUJ4eDi0tbXh5ubWYLuHhwfkcjkAID09HUZG\nRkhKSgIApKSkwMrKCh06dAAAjBw5Erm5uUhLS0NxcTGWLFkCAHBycoKOjo5CchUdHY0JEyY8c90D\nAgKwa9cuYTkuLg6Ghoawt7dHTU0NXF1dYW1tjV9++QVubm4YPnz4U68BYy8zTliY0tGgWmjwMyxN\nEola509rmjVrFszMzNC+ffun7rtnzx7Mnz8fzs7OMDIywoIFC/7yuZTbt29DX18f5ubmOH36NL77\n7jtoamo22M/Z2RnXrl3D77//juTkZAQHB6OgoADl5eVITEyEu7u7sG9QUBB0dHRgZGSEJUuWIDY2\nFrW1tQAAf39/IdEoKyvDkSNH4O/v/8x1DwgIQExMjDB8FR0dLZRz4sQJ2NnZISgoCLq6unj//ffR\nqVMnpKenP/X6Mfay4oSFKR0REX/p9i8Qtc6f1uTk5NTsfePj4xESEgJ9fX3o6+tj0KBBuHnzJgoL\nCxvdv1u3bigpKcHt27cRHR0t9ET07t0burq60NXVxenTp6GlpYV+/fohMTERSUlJcHd3h6urK06f\nPi0sP7Zq1SoMGTIEHTt2hEwmwx9//IFbt24BqEtY9u/fj8rKSuzfvx8ODg7CENez1N3c3BxWVlaI\niYnBw4cP8dNPPyEgIEAoJzk5WShHX18fubm5Qo8QY6wh/qnAlI4GEfewKCkNDY1Gh3xeeaK93njj\nDfz++++oqamBhoYGiouL8d///lfY7unpiTFjxmD8+PF/qy6XL19usM7d3R0JCQnIzMyETCaDu7s7\njh49irS0NAwcOBBA3XMqX331FeLi4vDmm2/i119/xZtvvinEZW1tDVNTUxw5cgTR0dFCktGSuj/u\nrampqYG1tTXMzMyEci5evIgjR478rWvA2MuEe1iY0tGoJYjbccKijBwcHPDzzz//5T69evVCp06d\nsG3bNhQVFSE0NFThFd3AwECsXLkSp06dQk1NDYqKihATE9Mq9XN3d8fOnTvRu3dvtGvXDh4eHvj2\n229hZmYGAwMDAEBBQQF0dHRgaGiIO3fuCM+vPCkgIABr165FcnIyxo0b1+K6+/n5IS4uDuHh4cJz\nMAAwZMgQXLp0CTt37kRJSQkqKiogl8tRUFAg7MMTYDKmiBMWpnRERBDzkJBSmjx5Mg4dOoSOHTvi\nq6++AtD4jLibN29GZGQkHB0dYWtrC2NjY2Gbl5cXli1bho0bN6Jz585wcXFBWlpak+d8lu+RuLi4\noKKiQuhNsbKygpaWlrAMAKNHj4anpyfs7e3x9ttvw9fXt8E5/P39kZSUhMGDB6Njx44trnvXrl3h\n6uqKlJQU+Pr6Cus1NDQgl8uRk5MDBwcHmJiYYPXq1QpJyuM6JScnQ1dXt9nXgDF1JSJO41XW076J\noaoudtbC3f+3A8OC/96QgapS13Zl7EVr6l7ie0w1cQ8LUzo8+SFjjLH6OGFhSkeDCBr8DAtjjLEn\ncMLClA73sDDGGKuPExamdMSk+JosY4wxxgkLUzoaRNDQ5LmEGGOM/Q8nLEzp1D3DwgkLY4yx/+GE\nhSkdHhJijDFWHycsTOlo1BLEmpywMMYY+x9OWJjSEQNox0NCSqlPnz4qO0GfWCzG9evX27oaSkMu\nlwuTOjKmCjhhaUUffvghunTpAhsbG2FdWFgYjI2NIZVKIZVKFSY7W79+PXr16gVra2ucOnVKWH/l\nyhX07dsXZmZm+Ne//vVCY1AGGrX8WrOy+uWXXxQ+c/+ibN++HW5ubs+tfA8PD2hpaUFXVxc2NjaY\nM2cOfv/99+d2PsbYs+OEpRVNnDgRR48eVVgnEokwa9YsZGZmIjMzE15eXgCAwsJCbNq0CQkJCdi8\neTOmT58uHDN79mzMnz8f6enpSExMxLlz515oHG1Ng4B27TXbuhpMSVRXVz/3c4hEInz99dcoKyvD\n/v37IZfLsXXr1ud+XsZY83HC0orc3Nygr6/fYH1jc1akpqZixIgRMDExgbu7O4gIDx48AADk5OTA\n19cXBgYGGDt2LFJTU5973ZWJmIgfulVSEokEJ06cAFDXe+jv749p06aha9euGD9+PK5cuQIAWLFi\nhcIsxwAwY8YMzJgxAwDw8OFDREREwNHREQMGDMC+ffuE+2T79u0YMGAAQkNDYWpqCj8/P4SEhCAl\nJQW6urrCZITV1dXYu3evMJFhREQEKisrhfPFxsbC0dERlpaW2LdvX7Nj7NWrF9555x0cOnSowbaK\nigpoaWkJvS///ve/0a5dO+HeXbx4MWbOnAkAOHz4MKRSKfT09DB06FDs3LlTKMfLywtff/21Qtl2\ndnY4ePAgAODOnTv47LPP0LNnT/j6+jb5f8CePXsgk8kU1q1ZswY+Pj7NukaMqRJOWF6ADRs2wNnZ\nGStWrEBZWRkAIC0tDVZWVsI+FhYWSE1NRW5uLgwNDYX11tbWOHv2bJNlF+hqqN0fnSpApwPPTquM\n6s9qvH//ftjZ2eHKlSvQ09PDF198AQDw8/NDbGys8IO8pqYG+/btw4QJEwAAixYtQnx8PPbu3Yuv\nv/4ay5YtQ3x8vFBuWloaqqqqcPHiRXz//fcIDw+Hi4sLysrKhGRh06ZN2LJlCzZs2IAffvgB33//\nPXbs2AGgbujq/fffx+LFixEbG4vt27c/NbbHCVNOTg7+85//YOzYsQ32efXVV+Ho6Ai5XA4ASExM\nhEQiEYZ0ExMT4eHhAQB47bXX8P333+P333/HnDlzMG3aNOTm5gIAAgICsGvXLqHcrKws5OfnY+TI\nkQCAkSNH4pVXXsG5c+fw/vvvw8vLS7iWTxo1ahRycnKEcgEgOjpauM5/dY0YUzX8a+xzFhISgiVL\nluD+/fuYO3cutmzZgjlz5jQ5g2h9T5tR9N/jPhL+3k/qAJnU4e9Xuo1V6uuhh1Gntq6G0hItbfjv\npCUo9O/PVmthYYFJkyYBAIKDgzFmzBgAgKmpKfr27YsDBw4gMDAQJ06cgLa2NhwdHUFEOHDgAJKS\nkoSHPoODg3Hw4EEMHToUQN1r7WFhYdDUrBsabOw+2Lt3L1auXInevXsDqOvB+eabbzBp0iTExsbC\n29sbb7/9NoC6YdYnnx9rcC2IMH36dMyZMwc9e/bEyJEj8dFHHzW6r7u7OxITE+Hj44NLly5h4cKF\nQqJy7tw54Rkfd3d34Zjhw4fDx8cHP/74I2bPno3Ro0cjJCQEv/76K7p3746oqCi88847aNeuHa5d\nu4aHDx9iwYIFAOqSF3d3dxw5cqRBr5WWlhZ8fHywa9cuLF68GNeuXUNOTg5GjRr11Gv0MpHL5UKS\nyVQXJyzP2ePeEj09PUydOhVTpkzBnDlz4OTkpPAbZXZ2NmQyGXR1dXHv3j1hfVZWFpydnZssf1Pk\nludXeaaUWiPRaC12dnbC37t27Yp79+6htrYWYrFY6EUIDAxU+K0/Ozsb+fn5sLW1FY6tra1Fjx49\nFMp9nKw0pry8HGfOnBF6JIC6pONx0p+WlgZXV1dhm1Qq/cs4RCIRNmzYgA8//FBhfVRUFCZPngwA\nGDhwIA4fPgx3d3fMmjULGRkZsLGxwZAhQxAcHIzU1FT07NlTGBa+fPkyVq1ahTNnzuDu3buorKyE\nWFzXqa2rq4uRI0di165dmDdvHnbv3o1vv/0WABAfH48bN24oDC/X1NSge/fuDRIWoK63Zvbs2Vi8\neDGio6MxZswYvPrqq0+9Ri8TDw8PoecLAJYuXdp2lWEtxkNCz9mdO3cA1I0lR0dHw9vbGwDg6OiI\nuLg45OfnQy6XQywWQ1e3bhjE0tISu3fvRnFxMQ4cOAAnJ6c2qz9jLfXuu+9CLpejoKAABw8eREBA\nAIC6XhljY2NkZWWhpKQEJSUlKC0txfnz54Vj6z/DpKGhodDLoqOjAycnJ8TFxQll/PHHHygpKQFQ\nd39lZmYK+2dkZLQohgkTJqCsrAxlZWU4fPgwAMDFxQU5OTk4cOAAPDw8YGVlhfz8fMTGxir8UJwz\nZw6MjY2RmJiI0tJSvPPOOwox+Pv7Y9euXUhJSUFFRQUGDRoEAPD09IS5ubkQV0lJCe7fv4/169c3\nWschQ4agqKgIFy5cwO7du4Xr/LRrxJiq4YSlFfn7+8PV1RU5OTno3r07IiMjMX/+fNja2sLZ2RlV\nVVUICQkBAHTp0gUhISHw9PTElClTsG7dOqGcVatWYeXKlZDJZHBzc0O/fv3aKiTGWqxz587w8PBA\nUFAQzMzMYGFhAaDueyi+vr6YP38+rly5gtraWuTl5f3l910cHBxw7do1hec4AgMDsWTJEmRkZKC2\nthYFBQU4duwYAMDb2xtHjx7F4cOHcf36daxdu/ap9X3a8Otj2tracHBwwNdffy0M+7i6uiI8PFxh\nGOj27dvo1KkT9PT0EBMTg5iYGIVyvL29cevWLYSGhsLPz09Yb2Fhgddeew2rVq3C3bt3UVVVhfT0\ndGRnZzdan3bt2mHcuHGYM2cOSkpKhGG1p10jxlQOMZXFzaeelLldJRIJJSQkEBFRWFgYBQYGCttu\n3LhBYrGYampqhHXfffcdiUQiWrVqlUI55eXlFBkZSe7u7qSnp0dSqZT27NlDRETbt28nNze3Buf+\n+OOPycTEhDp37kxERJWVlbRnzx4aOXIk6enpkZWVFW3YsEHY/6effiKZTEYWFha0b98+EovFlJeX\n12hcHh4eFBER0ezrsGDBAtLW1qbKykoiItq4cSOJxWIqLCwU9jl+/Di5urpSp06dyM/Pj+bOnatw\nvYiIgoODSSwW07lz5xTWFxQU0BdffEHW1tZkYGBAgwcPpgsXLjRZn+TkZBKJRDRt2jSF9X91jU6e\nPEndu3cX9vXy8qLly5c3+xqogqbuJWW+x1jTRETN/LWCKR2RSNTs3wqZ6uB2Zax1NHUv8T2mmnhI\niDHGGGNKjxMWxhhjjCk9TlgYY4wxpvQ4YWGMMcaY0uOEhTHGGGNKjxMWxhhjjCk9TlgYY4wxpvQ4\nYWGMMcaY0uOEhTH23G3fvh1ubm5tXY02JZFIkJCQ0Oi2mpoaTJs2Dd26dYOvr+9zOb9YLMb169cB\n1M0i//nnnz+X87TE8uXLX7oZpNmz49maGWPsKeRyOQIDA/Hrr7+2uAyRSNTkTMkpKSlITEzEtWvX\noKOj0+JzNNfmzZuf+zmexYIFC9q6CkwFcA8LY0wt1dbWKixXV1e3UU2e7tSpU7C1tW1RsqLMcT2m\nCnVkyo8TFsZYs6xYsQLjxo1TWDdjxgzMmDEDAPDw4UNERETA0dERAwYMwL59+1o8X8vJkyfx3nvv\noWPHjujduzcyMzMBAAUFBVi0aBEkEgkmTpworAeAoKAgzJw5E+PHj4eBgQFOnjwJiUSCTZs2wdXV\nFR06dBBmhp43bx5MTU0xadIkZGVlCWU8ePAAW7duhbOzMzp27IixY8fi4cOH8PLywu3bt6Grq4vX\nX38dd+/eBQAcOXIEo0aNgoWFBdasWaMwm3RKSgqGDRsGiUSCDRs2NBlraGgoQkNDsW/fPujq6mLb\ntm0AgJiYGAwdOhQ2NjYIDw/Hw4cPAQA3b96EWCzGvn370KdPH4XZmZ8UGxsLR0dHWFpaYt++fQrb\ngoKCsHjxYgB17fbRRx9BIpHAwMAAAwcOFPaTSCTYuHEjHB0dYW5ujvDwcFRVVQnbL168iMmTJ8PE\nxASzZ89Gfn6+wrFPXvuamhpERETAxcUFenp6sLS0xIkTJwAAYWFhCAwMFI7NyMhAUFAQJBIJFi9e\njNu3byuUGx4eDhcXF5iYmCAsLEyhTkyNteHEi+xv4uZTT8rarrdu3SJtbW0qKysjIqLq6moyMjKi\n1NRUIiKaOXMm+fn50Y0bN+j8+fPUp08fOnbsGBERbdu2jQYMGNCs82RkZJChoSFFR0dTZWUl5ebm\n0q1bt4iIaODAgTRt2jQqLCykiIgIev311+nRo0dERPTBBx+Qjo4ORUdHU1VVFVVUVJBEIiFra2tK\nSkqiiooKqq6uJkNDQ9q2bRvdv3+fduzYQcbGxsK5p0+fTiNGjKDMzEyqrq6mpKQkIiKSy+UK+xER\n/fjjj2Rra0spKSl0+/ZtGj9+PC1cuJCIiEpKSkhbW5t27NhBBQUFNGHCBGrXrp0w03V99We+PnHi\nBJmYmNDx48fp6tWrNHjwYAoNDSWiulmxRSIRjR49mvLy8qiioqJBeZcuXSIDAwOKiYmhvLw88vb2\nJpFIJMxWHRQURIsXLyaiupmmAwICqLS0lKqrq+nUqVNCOaampmRhYUHJycl0/vx5kkqlFB4eTkRE\nxcXFpK+vTwcPHqTS0lL64osvyNXVVTi2/rUvKioiY2Njunr1KhHV/Xt6XJ+wsDB67733iKhuJu/X\nXnuNvvnmGyosLKTp06eTu7u7Qrl2dnaUlpZGV69eJYlEQvHx8Y1e16buJWW9x9hf41ZTYXzTqaen\ntivQOn9aYMCAAbRz504iIjp27BiZm5sTEVFtbS1JJBLKz88X9l2zZg1NmTKFiJ4tYZk3bx7NmDGj\nwfqioiLS0tKiBw8eCOv69+9P+/fvJ6K6hMXT01PhGIlEQsuWLROWjx07RkOHDlXYx97entLS0qim\npoa6detGmZmZDc598uTJBglLQEAARUVFCcuZmZlkbW1NRER79uwhNzc3YVteXh6JRKImE5bQ0FDh\nBzZRXeK0YMECYfn48eNka2tLRP9LWB4nU41ZsWKFQgKUkJDQZMKyfv16GjZsGGVlZTUoRyKRCPsR\nEW3ZsoXeeustIiLaunUrTZo0Sdj2OBm8d++ecOyT1764uJgMDAzo0KFDVFlZ2WT8+/fvJxcXF2Fb\neXk5aWtrU3FxsVDu6tWrhe3/+Mc/aP78+Y1eB05Y1AsPCTGmalorZWmBgIAA7Nq1CwAQHR2NCRMm\nAACys7ORn58PW1tb6OvrQ19fH6GhoTh9+vQzn0Mul6N///4N1p89exZmZmYKz3n069cPp06dAlD3\nUKuTk1OD455cFx8fj+TkZKGO+vr6yM3NRVJSErKzs1FaWgp7e/tm1TM+Ph4hISFCOYMGDcLNmzdx\n7949pKamws7OTtjXzMwMenp6zb4GZ86cgYODg7Ds4OCAS5cuoaysrNG46ktLS1OIQyqVNtiH/u/f\nQHBwMDw8PPDWW2/BxsYGERERCvvVLyclJQVAXfxRUVFC/J06dUJ5eTmSkpIaraOBgQG+++47rFmz\nBkZGRvjkk09QVFTUoF6nT59G3759hWVtbW306tULZ86cabRORkZGKCgoaPJaMPXBCQtjrNneffdd\nyOVyFBQU4ODBgwgICAAAWFhYwNjYGFlZWSgpKUFJSQlKS0tx/vz5Zz7HoEGDhCTkSc7Ozrh+/TrK\ny8uFdenp6QqvS2toaDQ47pVX/vcypKenJzw8PIQ6lpSUoKysDLNnz4alpSX09PQUnot5slyql+R5\nenrim2++USirvLwcXbp0gZOTk0LseXl5KC0tbTLm+m8P9e/fH+fOnROWz507BxsbG+jq6jYaV32O\njo4KcWRkZDS5r7a2NhYsWIC8vDxERkZi1qxZCs/11C/H1dVViP/9999XiP/Bgwd49913m6yjl5cX\n4uPjkZWVhRs3bmDlypUN6jNgwAD8/PPPwnJ5eTmuXbsmnLe++u3C1BcnLIyxZuvcuTM8PDwQFBQE\nMzMzWFhYAKj7xoevry/mz5+PK1euCA+3PvnbdnP5+flh79692Lt3LyorK5Gbm4v8/Hx06tQJMpkM\nCxcuRGFhIbZv347Lly9j+PDhAJr3g2vIkCG4dOkSdu7ciZKSElRUVAgJmFgsxrhx47Bo0SKcP38e\nVVVVQv3t7OxQXFyMO3fuCGUFBgZi5cqVOHXqFGpqalBUVISYmBgAwLBhw5CRkYGoqCjcvn0bS5cu\n/csEo37dfXx8sGvXLpw4cQK5ubn48ssvMWbMmGZfQ29vbxw9ehSHDx/G9evXsXbt2ibPd+jQIeTm\n5qK2thY6OjrQ1NTEq6++Kuz3ww8/4PTp07h48SK2bt2Kt956CwAwfvx47N+/HwcPHkR5eTnKy8tx\n+PBhhQePn3T16lWcOHECf/75JzQ1NdG+fXuFBOyxoUOH4vLly4iMjERhYSEWLVoEmUwGAwODZsfP\n1BMnLIyxZxIQEICEhAShd+WxsLAwDBo0CCEhIejYsSPGjRsnvE1T/xskffr0EYaW6rO3t0dUVBT2\n798PQ0NDjB07FiUlJQCAqKgoaGtrQyaTQS6XIyEhAVpaWo2eozEaGhqQy+XIycmBg4MDTExMsHr1\nauEV6M8//xxvv/02goOD0aVLF6xfvx4A8Prrr2PevHkYOHAgOnbsiLt378LLywvLli3Dxo0b0blz\nZ7i4uCAtLQ0A0KFDB8TFxWHbtm1wcXGBo6MjjI2Nm6xX/bp7eHhgzZo1+OKLLzB69Gj4+Phg7ty5\nCvv/lT59+mDbtm1YunQpvL298cEHHygc8+T5cnNzMXToUOjp6WHSpEn4/PPPYWZmJuw3depUzJo1\nC6NHj0ZwcDCCgoIAAPr6+oiLi8PJkyfx5ptvolevXti5c2eTdfvzzz+xYMECdO7cGf369UOHDh0w\nc+bMBvXR0dHBiRMnkJiYCJlMBi0tLURFRTX72jH1JSLuT1NZIpGIu0PVELcrUxY9evRAREQEPD09\n27oqLdLUvcT3mGriHhbGGGOMKT1OWBhjjDGm9HhISIVxt6Z64nZlrHXwkJB64R4WxhhjjCk9TlgY\nY4wxpvQ4YWGMMcaY0mv6S0aMsTahr6/P35VgrBXo6+u3dRVYK+KHblUYPzjGGGPPjv/vVE08JNSK\nPvzwQ3Tp0gU2NjbCurKyMvj4+MDExASjR49W+Gz1+vXr0atXL1hbWyvMnXLlyhX07dsXZmZm+Ne/\n/vVCY1Amcrm8ravw3KhzbADHp+rUPT6mmjhhaUUTJ07E0aNHFdZt3rwZJiYmuHbtGoyNjREeHg4A\nKCwsxKZNm5CQkIDNmzdj+vTpwjGzZ8/G/PnzkZ6ejsTERIVJ0F4m6vyfpjrHBnB8qk7d42OqiROW\nVuTm5tZgzDQtLQ3BwcFo3749PvzwQ6SmpgIAUlNTMWLECJiYmMDd3R1EJPS+5OTkwNfXFwYGBhg7\ndqxwDGOMMfay4oTlOUtPT4elpSUAwNLSUpgcLTU1FVZWVsJ+FhYWSE1NRW5uLgwNDYX11tbWOHv2\n7IutNGOMMaZk+C2h5+xZHuxq7M2Qpx2v7m+TLF26tK2r8Nyoc2wAx6fq1D0+pno4YXnOZDIZrly5\nAqlUiitXrkAmkwEAnJycEB8fL+yXnZ0NmUwGXV1d3Lt3T1iflZUFZ2fnRsvmp9wZY4y9LHhI6Dlz\ncnJCZGQkHj16hMjISCH5cHR0RFxcHPLz8yGXyyEWi6Grqwugbuho9+7dKC4uxoEDB+Dk5NSWrIdz\nwQAACIFJREFUITDGGGNtjhOWVuTv7w9XV1dcvXoV3bt3x7Zt2xASEoL8/HxYWFigoKAAkydPBgB0\n6dIFISEh8PT0xJQpU7Bu3TqhnFWrVmHlypWQyWRwc3NDv3792iokxhhjTDkQUzmJiYlkaWlJPXv2\npPXr17d1dVrM1NSUbGxsyN7enmQyGRER3b9/n0aNGkXdu3cnHx8fKisrE/Zft24d9ezZk6ysrCg5\nObmtqt2kiRMnkqGhIfXp00dY15J4srKySCqVUo8ePWjhwoUvNIamNBZbaGgovfHGG2Rvb0/29vYU\nGxsrbFOl2IiI8vPzycPDg6ytrcnd3Z2ioqKISH3ar6n41KUNHz16RI6OjmRnZ0dOTk701VdfEZH6\ntB+rwwmLCrK3t6fExES6efMmWVhYUFFRUVtXqUUkEgn99ttvCutWrFhB06ZNo4qKCpo6dSp9+eWX\nRER07949srCwoFu3bpFcLiepVNoWVf5LSUlJlJGRofBDvSXxeHl50e7du6m4uJj69+9P6enpLzyW\n+hqLLSwsjFavXt1gX1WLjYjozp07lJmZSURERUVF1KNHD7p//77atF9T8alTG5aXlxMRUUVFBfXu\n3ZuuXr2qNu3H6vCQkIopLS0FAAwcOBCmpqYYNmyYSn+nheo9OPws360pKytriyo3SZ2/w9NYbEDj\nD36rWmwA0LVrV9jb2wMAOnXqhN69eyM9PV1t2q+p+AD1aUNtbW0AwIMHD1BdXY327durTfuxOpyw\nqJgnv+sCqPZ3WkQiETw9PTF69GjExMQAeLbv1jzepszU/Ts8GzZsgLOzM1asWCEkkGlpaSodW25u\nLi5fvgxHR0e1bL/H8T1+mF9d2rC2thZ2dnbo0qULpk2bBhMTE7Vsv5cZJyyszZw+fRoXLlzA8uXL\nMWvWLNy9e/dvf7dG2fzdeJ7l+BctJCQEN27cQFxcHPLy8rBlyxYAjddZVWIrKyuDr68v1qxZg9de\ne03t2u/J+HR0dNSqDcViMS5cuIDc3Fxs2rQJmZmZatd+LztOWFSMTCZDdna2sHz58uUmv9Oi7IyM\njAAAVlZWGDVqFH766SfhuzUAGny3JisrSzj28XdrlN2zxtOzZ89mf4enrRkaGkIkEkFPTw9Tp07F\ngQMHAKhubFVVVXjnnXcQGBgIHx8fAOrVfo3Fp25tCAASiQTe3t5ITU1Vq/ZjnLCoHD09PQBAUlIS\nbt68iePHj6vkd1oePnwodD8XFRUhLi4OI0aMaNF3a5SZOn+H586dOwCA6upqREdHw9vbG4BqxkZE\nCA4ORp8+ffDJJ58I69Wl/ZqKT13asLi4GH/88QcA4LfffsOxY8fg4+OjNu3H/s+LfMKXtQ65XE6W\nlpZkbm5O69ata+vqtMj169fJzs6O7OzsyNPTkyIiIojor19DXLt2LZmbm5OVlRUlJSW1VdWb5Ofn\nR0ZGRqSpqUnGxsYUGRnZonguX75MUqmUJBIJffrpp20RSgOPY2vXrh0ZGxtTREQEBQYGko2NDTk4\nONDMmTMV3vhSpdiIiJKTk0kkEpGdnZ3wiu+RI0fUpv0aiy82NlZt2vDixYsklUrJ1taWhg0bRjt2\n7CCilv1/oozxsToiIh6kY4wxxphy4yEhxhhjjCk9TlgYY4wxpvQ4YWGMMcaY0uOEhTHGGGNKjxMW\nxhjCwsKwevXqFh174cIFHDlypEVlDRs2THgdtTn27t2LL7/88pnryBhTfZywMMb+1leDMzMzERsb\n+8xlZWZmwtDQEB06dGj2uUaPHo09e/bwF0gZewlxwsLYS2rXrl3o27cvBgwYgPz8fGF9QUEB5s6d\nCxcXF3zwwQe4ceMGACAoKAizZs2Cg4MDPDw8kJqaiqqqKixZsgR79uyBVCrF3r17AQB5eXkYNGgQ\nbG1tsXv37kbPv3XrVgQEBAAAbt68CWtra3z88cfo1asX/vnPf+LMmTPo378/+vfvj19++QUAoKmp\nCalUiuPHjz/PS8MYU0KcsDD2EiouLkZoaChiY2MRHR2NuLg4oWdkyZIl8PPzQ0pKCnx9fbFy5Urh\nuHPnzkEul2PlypWYNGkS2rVrh88++wx+fn7IzMzE+PHjQUQ4deoUDhw4gMOHD2PRokWN1uH8+fMK\nE9BlZ2djwoQJuHz5Mo4dO4a1a9fi5MmT+Pjjj4U5boC6qRwyMjKe05VhjCmrV9q6AoyxF+/xVAhd\nu3YFAAwZMgRA3SfaY2NjG00IRCIRxowZA11dXTg6OoKIUFBQACJSGKIRiUTw8fFBhw4d0KFDB2ho\naKCwsFBhFlygrhfGxMREWH7jjTfg7u4OAOjXrx8GDx4MTU1NuLi4YMOGDcJ+5ubmOHjwYOtdDMaY\nSuCEhbGXkEgkavQ5kNraWojFYpw9exbt27dvsL3+MSKRqNFnVp58LkVTUxMVFRWN1uPJ8uof83je\nLE1NTfz5558KdVSFmboZY62Lh4QYewkNHz4cx44dw7179/Drr78iISEBQF1y4O3tjc2bN6OmpgZE\nhIsXLwKoSy5+/PFHPHjwAOnp6RCLxejWrRtMTU1RVFT0zHUwNzfHrVu3nvm469evw8LC4pmPY4yp\nNk5YGHsJGRgYYOnSpfDy8oK/vz+GDx8ubFu6dCnu3r2Lfv36oU+fPoiJiQFQ15vi4OAAd3d3zJ07\nF9988w0AwMXFBWVlZQoP3TanB8Te3h7Z2dnCcv1jnlx+8u/Z2dmQSqUtiJoxpsp48kPGWLNMnDgR\nb7/9NsaOHdsq5WVmZmLVqlWIiopq9jF//vknXF1dce7cOR4WYuwlwz0sjLE2IZVKUVRUhNLS0mYf\n8+OPP8Lf35+TFcZeQtzDwhhjjDGlxz0sjDHGGFN6nLAwxhhjTOlxwsIYY4wxpccJC2OMMcaUHics\njDHGGFN6nLAwxhhjTOn9f0kaW/yXx/K5AAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x4a9ac10>"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Should check the dispersion curve obtained for a constanst Q with Kartjansson's equation and the value of Q obtained for Gillian's equation\n",
      "\n",
      "Let's take fdom=7.5 Hz and Q = 60 and Vp = 3300 (2nd layer values)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Kartajansson's dispersion curve\n",
      "f = np.linspace(0,100,1000) #Hz\n",
      "w = 2*np.pi*f #rad/s\n",
      "fdom = 50 #Hz\n",
      "w0 = 2*np.pi*fdom  #rad/s\n",
      "Vp0 = 3300.0 #m/s\n",
      "Qp = 60.0  \n",
      "\n",
      "vp_de_omega = Vp0*(w/w0)**(np.arctan(1/Qp)/np.pi)\n",
      "#vp_de_omega =Vp0*(w/w0)**(1/Qp)\n",
      "plt.plot(f, vp_de_omega)\n",
      "plt.ylim((3200, 3400))\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 22,
       "text": [
        "(3200, 3400)"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD9CAYAAABazssqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHa1JREFUeJzt3X9Q1XW+x/EnJyJgYxUxwFaQFOJHtXDwx8E1lcxM2VXn\nag62/uCmNXvVbj8226al1nS2nF3bAW3Wlpp1a1qL29TlrtOuSpbHbMuDrjd1QRFakOwaBv5IFDTg\ne//4xgEUQeXHOZzv6zHzne85x8M57+9nnPfrfH/7GYZhICIilmTzdAEiIuI5CgEREQtTCIiIWJhC\nQETEwhQCIiIWphAQEbGwTkOgoaEBh8NBSkoKaWlp5OTktPv33/3ud9hsNk6cOOF+bd26dcTFxZGU\nlMTHH3/sfv3gwYOkpqYyfPhwsrOze3gxRETkmhhdOHv2rGEYhtHQ0GDcdtttRllZmWEYhlFVVWXc\ne++9RkxMjFFbW2sYhmFUV1cb8fHxxpEjRwyn02nY7Xb350ybNs3Iz883ampqjHHjxhm7d+/u6qtF\nRKSXdbk5KDg4GIC6ujoaGxu54YYbAPj5z3/Ob3/723bvdblcTJ06lejoaCZOnIhhGNTV1QFQWlpK\nZmYmYWFhzJo1C5fL1dN5JiIiV6nLEGhubiY5OZmIiAgefvhhoqKi+Mtf/sLQoUP54Q9/2O69RUVF\nJCYmup/Hx8fjcrkoLy8nPDzc/XpSUhK7du3qwcUQEZFr4d/VG2w2G/v27aOyspKMjAxGjhzJ6tWr\nKSwsdL/H+O7KE0YHV6Dw8/O75LWO3ne594qISNcu11e7csVHB8XExJCRkcE//vEPKioqSE5O5pZb\nbuHo0aOMHDmS6upqHA4HJSUl7r85dOgQo0ePJjY2lurqavfrJSUlpKWlXXZBNBmsWLHC4zV4y6Sx\n0FhoLDqfuqPTEKipqeHUqVMA1NbWUlhYyOzZs6murqaiooKKigqGDh3K3r17iYiIYMyYMWzdupWq\nqiqcTic2m42QkBAAEhISyM/Pp6amhoKCAhwOR7cKFxGR7ut0c9CxY8fIysqiqamJyMhIli9fzpAh\nQ9q9p+0mnIiICJYsWcKkSZMICAggLy/P/W8vvvgi8+fP5+mnn2bu3LmMGjWqhxdFRESulp/R3XWJ\nHuTn59ftVRtf4XQ6SU9P93QZXkFj0Upj0Upj0ao7vVMhICLSz3Wnd+qyESIiFqYQEBGxMIWAiIiF\nKQRERCxMISAiYmEKARERC1MIiIhYmEJARMTCFAIiIhamEBARsTCFgIiIhSkEREQsTCEgImJhCgER\nEQtTCIiIWJhCQETEwhQCIiIWphAQEbEwhYCIiIUpBERELEwhICJiYQoBERELUwiIiFiYQkBExMIU\nAiIiFqYQEBGxMIWAiIiFKQRERCxMISAiYmEKARERC1MIiIhYmEJARMTCFAIiIhamEBARsTCFgIiI\nhSkEREQsTCEgImJhCgEREQtTCIiIWJhCQETEwhQCIiIWphAQEbEwhYCIiIUpBERELEwhICJiYQoB\nEREL6zQEGhoacDgcpKSkkJaWRk5ODgDPPvssycnJpKSksGDBAmprawGorKwkKCgIu92O3W5n6dKl\n7s86ePAgqampDB8+nOzs7F5cJBERuVJ+hmEYnb3h3LlzBAcHc/78eUaOHElBQQGRkZGEhIQAsGrV\nKhobG1m1ahWVlZVMnz6dAwcOXPI5GRkZZGVlMXnyZGbOnElubi6jRo1qX4yfH12UIyIiF+lO7+xy\nc1BwcDAAdXV1NDY2EhgY6A6AxsZGzp49S2BgYJdfVFpaSmZmJmFhYcyaNQuXy3VNBYuISM/x7+oN\nzc3N2O12iouLyc3NJSoqCoDs7Gzy8vKIj49n+/bt7vdXVFSQkpKCw+Fg6dKlJCcnU15eTnh4uPs9\nSUlJbNy4kWXLll3yfc8995z7cXp6Ounp6d1YPBER3+N0OnE6nT3yWV1uDmpRWVlJRkYGGzduxG63\nA+amopbt+zk5OVy4cIGzZ88SGhrK5s2beeqpp9i/fz9lZWUsXLiQTz/9FIDNmzfz5ptv8sYbb7Qv\nRpuDRESuWq9uDmoRExNDRkZGu804wcHBLFq0yN3cAwICCA0NBWDatGn4+/tTXl5OXFwc1dXV7r8r\nKSkhLS3tmgoWEZGe02kI1NTUcOrUKQBqa2spLCxk5syZlJWVAeY+gbfeeotZs2a539/U1ATA3r17\nqa+vJzY2FoCEhATy8/OpqamhoKAAh8PRawslIiJXptN9AseOHSMrK4umpiYiIyNZvnw5Q4YM4b77\n7qO0tJSgoCDS09N56KGHAPjoo4/41a9+hb+/P7GxseTl5bk/68UXX2T+/Pk8/fTTzJ0795Ijg0RE\npO9d8T6BvqB9AiIiV69P9gmIiIjvUQiIiFiYQkBExMIUAiIiFqYQEBGxMIWAiIiFKQRERCxMISAi\nYmEKARERC1MIiIhYmEJARMTCFAIiIhamEBARsTCFgIiIhSkEREQsTCEgImJhCgEREQtTCIiIWJhC\nQETEwhQCIiIWphAQEbEwhYCIiIUpBERELEwhICJiYQoBERELUwiIiFiYQkBExMIUAiIiFqYQEBGx\nMIWAiIiFKQRERCxMISAiYmEKARERC1MIiIhYmEJARMTCFAIiIhbm7+kCREQEmpuhrg5OnYLTp1vn\nbR9fbt4dfoZhGD2zCN3n5+eHF5UjInLFDAMaGuDkSXM6dar18cXPWx63beTffAPBwTBwIAwY0Dpv\n+7ij1wYMgNtvv/beqRAQEflOczOcOXP5pt1ZQz95Emw2szmHhprT5R63PG/b3L//fbjuumuruzu9\nUyEgIj6nsdFsyrW1cOKEOW/7uGXetrGfPNn6a7yrJt5RQw8NhcBAzyyvQkBEfFJzs/lLu6MG3tm8\nrs5szGFhMGiQOW/7uGXetpGHhpq/yv374Z5ShYCIeL2Whv711+ZUU3Pp/OJf7KdPw403dt3ML/63\nAQPMTTNWoRAQkT53/vyljbyj5t7y+MQJs6EPHgw33dTx/OImHxraP3+Z9zWFgIh0W1OT2bCPH4fq\n6vZTSzNv29wbGszGfbmmfvFrgwfD9dd7eil9k0JARDp04ULHTb26+tLXT540f3lHREB4uDlvmcLD\nL23qAwaAn5+nl1BAISBiKc3N5q/xY8fg//7PnB87Bl99dWmjr6szm3bbhn5xc295PHiwNr30VwoB\nER/Qsjnm4ube8rhlXl1tHlM+ZAjcfLM5HzIEIiMvbfKhodbaQWpVvRYCDQ0NTJw4kfPnzxMYGEhm\nZiaPP/44zz77LJs2bcLPz4877riD3NxcwsLCAFi3bh0vvfQS119/Pa+88gp33nknAAcPHmTevHmc\nOnWK+++/n+eff75HF0TEWxmGefz50aPwxRfm/MsvL230X39tHtbYtrF39DgyEm64wdNLJd6kV9cE\nzp07R3BwMOfPn2fkyJEUFBQQGRlJSEgIAKtWraKxsZFVq1Zx/PhxJkyYQGFhIRUVFTz++OPs3bsX\ngIyMDLKyspg8eTIzZ84kNzeXUaNG9diCiHiCYZjb0o8ebT+1NPuWCSAqCoYONec339w6tTT5iAgI\nCPDs8kj/1J3e2eUWwODgYADq6upobGwkMDDQHQCNjY2cPXuWAQMGAOByuZg6dSrR0dFER0djGAZ1\ndXXceOONlJaWkpmZCcCsWbNwuVyXhICIt6mrgyNHWqeLm/3Ro+YRLy3NfehQc7rzzvavff/7nl4S\nkY51GQLNzc3Y7XaKi4vJzc0lKioKgOzsbPLy8oiPj8fpdAJQVFREYmKi+2/j4+NxuVwMGzaM8PBw\n9+tJSUls3LiRZcuWXfJ9zz33nPtxeno66enp17hoIp0zDPPkpZYGX1l56eNz5yA6GoYNM6eoKEhP\nb23uP/gBfPebSKTPOJ1Od9/tri5DwGazsW/fPiorK8nIyGDcuHHY7Xaef/55srOzyc7O5he/+AU5\nOTkdro74dXAMWWerLW1DQKS7Tp+Gzz+Hf/2rtbG3bfaGYTb3mJjWRj92bOtrN92kwyDF+1z8A3nl\nypXX/FlXfEBYTEwMGRkZuFwu7HY7YG4qWrRoEQ899BAADoeDbdu2uf/m0KFDjB49mpCQEKqrq92v\nl5SUkJaWds1Fi7RobjZ3qrY0+s8/b/+4oQGGD4cRI8ymHhsLd9/d2uQHDlSTF2vrNARqamrw9/dn\n4MCB1NbWUlhYyBNPPEFZWRlxcXE0Njby1ltvMWvWLADGjBnDk08+SVVVFf/617+w2Wzu/QcJCQnk\n5+czefJkCgoKyM3N7f2lE5/Q2AgVFXD4cGuTb2n0FRXmSUsjRrROGRmtj8PD1eRFOtNpCBw7doys\nrCyampqIjIxk+fLlDBkyhPvuu4/S0lKCgoJIT093rwlERESwZMkSJk2aREBAAHl5ee7PevHFF5k/\nfz5PP/00c+fO1U5haccwzEMkS0vNZl9a2vq4osI8LPLWW81f8iNGwF13mfPhw+F73/N09SL9l04W\nkz51/rzZ2A8dat/wDx82f7HHx5vNPj6+9XFsLAQFebpyEe+lM4bF61y4YDb24uL2U2WluS0+MbF9\no4+PN68aqU03IldPISAe09hoNvt//tNs8iUl5ryiwtz5mpQEt93WOt16q852FelpCgHpE6dPw/79\nsG8ffPaZOR08aJ71evvt7Zt9fLyavUhfUQhIjzIM87IHLY2+pelXV5vNPiUFkpPN+R13mDcKERHP\nUQhIt1RXw+7drdOePeaVJ1NTW5t9Soq5g/a66zxdrYhcTCEgV+z0abPJt236dXUwahSMHm1Oo0aZ\nl0PQTlqR/kEhIB0yDPOkqr//vXU6cgTs9taGP3q0eby9Gr5I/6UQEMA8LHPv3taG/8kn5p2ixo1r\nnZKTdfcoEV+jELCob781N+d8+CFs3w4ul7ndvm3Tj47Wr3wRX6cQsIimJvjf/zUb/ocfmr/2Y2PN\nSyjcdReMH29eR0dErEUh4MO+/BI2b4YtW+CDD8wdtnfdBZMmwcSJMGiQpysUEU9TCPiQCxfMbfmb\nN5vTl1/ClCkwbZo5j4z0dIUi4m0UAv3cqVPw179CQQFs2wZxcWbTnzYNxozRsfki0jmFQD/01Vfw\nl7/Af/83fPqpuWnn3/4Nfvxj84bjIiJXSiHQT3z1Fbz9NvzXf5kXWps2zWz8U6fqPrUicu0UAl7s\nm2/MzTxvvglFRTBjBmRmmrc41AXWRKQnKAS8TFMTbN0Kr79uzidOhHnz4Cc/geBgT1cnIr5GIeAl\nqqpgwwZzGjIEFi2COXN0GKeI9K7u9E5dQKCbmprgvfcgL888Y/enPzWf//CHnq5MRKRrCoFrdOaM\n+Yt/3ToYPBiWLYN33tHmHhHpXxQCV+nIEbPxv/aauXP3z3+GsWM9XZWIyLWxebqA/uLzz2HxYvNG\nK35+5tU6335bASAi/ZvWBLpw+DA8/7x5Ru+yZVBWph29IuI7FAKXcewY/OpX8D//A488AuXlMHCg\np6sSEelZ2hx0kbo6WLHCvKH6wIHmmsCzzyoARMQ3KQS+YxiwcSPceqv5q/8f/4A1ayA01NOViYj0\nHm0Owvy1v3Qp1NaaF3RLS/N0RSIifcPSawLffgurVsGPfmRevXP3bgWAiFiLZdcESkthwQIICzNv\n2RgV5emKRET6nuXWBAwD1q+HO++EBx6Av/1NASAi1mWpNYEzZ8yLulVUwMcfQ3y8pysSEfEsy6wJ\nlJaCw2Ee6qkAEBExWSIEtm6F8ePh5z+HV1+FwEBPVyQi4h18fnPQ66/DU0+Zd/caN87T1YiIeBef\nDoHVq+GVV2D7dkhM9HQ1IiLex2dDYOVK8yqff/873Hyzp6sREfFOPhkCLQHw4YcQEeHpakREvJfP\nhUBuLuTng9OpABAR6YpPhcC775oXffvkEwWAiMiV8JkQ+PRT+I//MA8HHTbM09WIiPQPPnGeQHU1\nzJkDf/qTeftHERG5Mn6GYRieLqKFn58fV1tOUxNMmWLe6/fXv+6lwkREvNi19M4W/X5N4PnnzfnK\nlZ6tQ0SkP+rXawL79sE995iXgv7BD3qxMBERL2bJNYFvvzUvBf2b3ygARESuVb8Ngd/9DsLD4d//\n3dOViIj0X/1yc9CxY3DHHVBUBMOH90FhIiJerNc2BzU0NOBwOEhJSSEtLY2cnBwAnnzySRITE0lN\nTeWxxx6jvr4egMrKSoKCgrDb7djtdpYuXer+rIMHD5Kamsrw4cPJzs6+pmJbPPOMeXMYBYCISPd0\nuSZw7tw5goODOX/+PCNHjqSgoIDKykruvvtuAH72s5+RlpbG4sWLqaysZPr06Rw4cOCSz8nIyCAr\nK4vJkyczc+ZMcnNzGTVqVPtiriDNDhyAyZPNm8QMHHi1iysi4nt6dcdwcHAwAHV1dTQ2NhIYGMg9\n99yDzWbDZrNx7733smPHji6/qLS0lMzMTMLCwpg1axYul+uaCn7hBVi+XAEgItITurxsRHNzM3a7\nneLiYnJzc4m66K7sr776Kg8++KD7eUVFBSkpKTgcDpYuXUpycjLl5eWEh4e735OUlMTGjRtZtmzZ\nJd/33HPPuR+np6eTnp7ufl5WBtu2mfcIEBGxKqfTidPp7JHP6jIEbDYb+/bto7KykoyMDMaNG4fd\nbgdg1apVhISEMGfOHABuvvlmvvjiC0JDQ9m8eTMLFixg//79l6ymdLba0jYELvbb38KyZRASciWL\nJiLimy7+gbyyG2fLXvEhojExMWRkZFBUVATAa6+9xtatW/nzn//sfk9AQAChoaEATJs2DX9/f8rL\ny4mLi6O6utr9vpKSEtLS0q6q0JoaeOcd+M//vKo/ExGRTnQaAjU1NZw6dQqA2tpaCgsLmTFjBlu2\nbGHNmjVs2rSJwDZ3ba+pqaGpqQmAvXv3Ul9fT2xsLAAJCQnk5+dTU1NDQUEBDofjqgp9/XWYMQPC\nwq7qz0REpBOdHh104MABsrKyaGpqIjIyknnz5rFw4ULi4uK4cOECgwYNAmDs2LGsX7+ed999lxUr\nVuDv709sbCyPPPIIEyZMAMxf//Pnz+fkyZPMnTuX1atXX1rMZfZwGwYkJJhXCf3Rj3pq0UVEfEN3\njg7qFyeL7dxp3ivgn/8EPz8PFCYi4sV8/tpBb78NP/2pAkBEpKd5/ZpAczMMHWreM/jWWz1Tl4iI\nN/PpNYFPPoGbblIAiIj0Bq8PgXffhdmzPV2FiIhv8voQ2LoVfvxjT1chIuKbvDoEjh6F48fhuxOU\nRUSkh3l1CLz/Ptx9N9i8ukoRkf7Lq9vr9u1mCIiISO/w6hDYtQvGjvV0FSIivstrzxOorYVbboGT\nJ+G66zxcmIiIF/PJ8wSKimD0aAWAiEhv8toQ2LULrvJCoyIicpW8NgQ++wxSUz1dhYiIb/PaECgp\ngdtu83QVIiK+zSt3DNfXw6BB8M03cP31nq5KRMS7+dyO4cOHYcQIBYCISG/zyhAoLoakJE9XISLi\n+7wyBLQ/QESkb3hlCJSVQVycp6sQEfF9XhkCR45ATIynqxAR8X1eGwLDhnm6ChER3+d1h4jW1xsM\nHAjnzukS0iIiV8KnDhGtqjJvLK8AEBHpfV7XaquqIDra01WIiFiD14VAdTUMGeLpKkRErMHrQuD4\ncQgP93QVIiLW4HUhUF2tEBAR6SteFwJaExAR6TsKARERC/PKEIiI8HQVIiLW4JUhoDUBEZG+4XUh\nUFMDYWGerkJExBq87rIRNptBYyP4+Xm6GhGR/sGnLhvxve8pAERE+orXhUBIiKcrEBGxDoWAiIiF\nKQRERCxMISAiYmEKARERC1MIiIhYmEJARMTCFAIiIhamEBARsTCFgIiIhXldCNx4o6crEBGxDq8L\ngcBAT1cgImIdnYZAQ0MDDoeDlJQU0tLSyMnJAeDJJ58kMTGR1NRUHnvsMerr691/s27dOuLi4khK\nSuLjjz92v37w4EFSU1MZPnw42dnZl/3OgIDuLpJvcDqdni7Ba2gsWmksWmksekanIRAYGMj27dv5\n7LPP2LFjB3/84x8pKytjypQpFBcXs2fPHs6ePcubb74JwPHjx1m/fj0ffPABL7/8Mo888oj7s554\n4gmeeuopdu/ezY4dO9izZ0+H33nDDT24dP2Y/oO30li00li00lj0jC43BwUHBwNQV1dHY2MjgYGB\n3HPPPdhsNmw2G/feey87duwAwOVyMXXqVKKjo5k4cSKGYVBXVwdAaWkpmZmZhIWFMWvWLFwuV4ff\npzUBEZG+02UINDc3k5ycTEREBA8//DBRUVHt/v3VV19l+vTpABQVFZGYmOj+t/j4eFwuF+Xl5YS3\nuWdkUlISu3bt6vD7tCYgItJ3/Lt6g81mY9++fVRWVpKRkcG4ceOw2+0ArFq1ipCQEObMmQPQ4Z1t\n/Dq4Q0xnd8C5+27dUabFypUrPV2C19BYtNJYtNJYdF+XIdAiJiaGjIwMioqKsNvtvPbaa2zdupUP\nPvjA/R6Hw8G2bdvczw8dOsTo0aMJCQmhurra/XpJSQlpaWmXfIcX3elSRMQSOt0cVFNTw6lTpwCo\nra2lsLCQGTNmsGXLFtasWcOmTZsIbHNM55gxY9i6dStVVVU4nU5sNhsh3539lZCQQH5+PjU1NRQU\nFOBwOHpxsURE5Ep0eqP5AwcOkJWVRVNTE5GRkcybN4+FCxcSFxfHhQsXGDRoEABjx45l/fr1AKxd\nu5aXXnqJgIAA8vLyGD9+PGD++p8/fz4nT55k7ty5rF69ug8WT0REOmV4iR07dhgJCQlGbGyssW7d\nOk+X06eqqqqM9PR0IykpyZg4caKxceNGwzAM45tvvjFmzJhhREVFGTNnzjTOnDnj4Ur7RmNjo5GS\nkmL85Cc/MQzDuuNQV1dnLFy40IiLizMSExONXbt2WXYsXnnlFWPs2LFGamqq8eijjxqGYZ3/Fw88\n8IARHh5u3H777e7XOlv2tWvXGrGxsUZiYqKxc+fOLj/fa84YfvTRR8nLy2Pbtm38/ve/p6amxtMl\n9Znrr7+enJwciouLeeedd3jmmWc4c+YML7/8MtHR0ZSVlTF06FD+8Ic/eLrUPrF27VqSkpLcBxVY\ndRxWrFhBdHQ0+/fvZ//+/SQkJFhyLE6cOMELL7zA+++/z+7duzl8+DBbt261zFg88MADbNmypd1r\nl1v2zs7VuhyvCIHTp08DMGHCBIYNG8aUKVMuex6BL4qMjCQlJQWAwYMHc9ttt7F7926KiopYvHgx\nN9xwA4sWLbLEmBw9epS//e1vPPjgg+4DBaw4DgDbtm3jl7/8JYGBgfj7+zNgwABLjkVQUBCGYXD6\n9Gnq6+s5d+4cAwcOtMxYjB8/ntDQ0HavXW7ZOzpX68yZM51+vleEwO7du0lISHA/7+w8Al9XXl5O\ncXExY8aMaTcuCQkJFBUVebi63vf444+zZs0abLbW/5pWHIejR4/S0NDAkiVLcDgc/OY3v6G+vt6S\nYxEUFMTLL79MTEwMkZGRjBs3DofDYcmxaHG5ZXe5XJecq9XVuHhFCIjpzJkzZGZmkpOTw4033mi5\nQ2bfe+89wsPDsdvt7ZbdauMA5nW7Dh8+zOzZs3E6nRQXF/P2229bciy+/vprlixZQklJCZWVlXz6\n6ae89957lhyLFlez7B2dq9WWV4TA6NGjOXTokPt5cXFxh+cR+LJvv/2W2bNns2DBAmbOnAmY43Lw\n4EHAvADf6NGjPVlir/vkk0/YtGkTt9xyC/fffz8ffvghCxYssNw4AMTGxhIfH8/06dMJCgri/vvv\nZ8uWLZYci6KiItLS0oiNjSUsLIw5c+awc+dOS45Fi8stu8PhoKSkxP2+lnO1OuMVITBgwAAAPvro\nIyorK3n//fctdR6BYRgsXryY22+/nccee8z9usPhYMOGDdTX17NhwwafD8YXXniBL774goqKCvLz\n85k0aRJvvPGG5cahRVxcHC6Xi+bmZv76178yefJkS47F+PHj2bNnDydOnOD8+fNs3ryZKVOmWHIs\nWlxu2Ts7V+uyevBIpm5xOp1GQkKCMWLECGPt2rWeLqdP7dy50/Dz8zOSk5ONlJQUIyUlxdi8ebNl\nDoHriNPpNKZPn24YhnUOBbxYaWmp4XA4jOTkZOOJJ54w6urqLDsWf/rTn4wJEyYYo0aNMp555hmj\nqanJMmMxd+5cY8iQIUZAQIAxdOhQY8OGDZ0ue25urjFixAgjMTHR+Oijj7r8/E5PFhMREd/mFZuD\nRETEMxQCIiIWphAQEbEwhYCIiIUpBERELEwhICJiYf8PrYB/pkufS/4AAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0xcdf6390>"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Another dispersion model\n",
      "tauSigma = (np.sqrt(Qp**2+1) - 1)/(w0*Qp)\n",
      "tauEpsilon = 1/(w0**2*tauSigma)\n",
      "Vp_dispersive = Vp0/np.sqrt((1+w0**2*tauSigma*tauEpsilon)/(1+w0**2*tauSigma**2))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Vp_dispersive"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 24,
       "text": [
        "3272.842769431034"
       ]
      }
     ],
     "prompt_number": 24
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Vp0*(7.5*2*np.pi/w0)**(np.arctan(1/Qp)/np.pi)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 25,
       "text": [
        "3266.9566586050573"
       ]
      }
     ],
     "prompt_number": 25
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Using dispersion relation in time domain code only makes sense if fdom is different from the attenuation mechanism central frequency, if they are the same, then there is no reason to introduce dispersion as we are looking for smthing."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}