{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Example 2: Lower crustal anisotropy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this example we generate P receiver functions for a model that includes a lower-crustal anisotropic layer. This example follows that of Figure 2 in [Porter et al. (2011)](#references), which uses the Raysum software developed by [Frederiksen and Bostock (2000)](#references])\n", "\n", "Start by importing the necessary modules" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from obspy.core import Stream\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_Porter2011.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\n", "\n", "

\n", " Danger! Using 'SH' will not work properly for modeling receiver functions as the code will think you want plane-wave displacements (see below). Do not use 'SH' if you want S-wave receiver functions.\n", "

" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "wvtype = 'P'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we use variables to define the desired time series. \n", "\n", "
\n", " Warning! Be careful to use a total length of time large enough to avoid wrap around effects. Sometimes if you see signals arriving at aberrant (early) times, try with either (or both) a greater number of samples or higher sample distance.\n", "
" ] }, { "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": [ "Now specify the parameters of the incident wavefield in terms of a horizontal slowness and back-azimuth. In this example the slowness won't change, so we can pass it as a global variable now. The back-azimuth will range from 0 to 360 degrees with a 10-degree increment, so we define a `np.ndarray` for this variable and do not yet pass it as a global variable." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "slow = 0.06 # Horizontal slowness (or ray parameter) in s/km \n", "baz = np.arange(0., 360., 10.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Read the model parameters 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": 6, "metadata": {}, "outputs": [], "source": [ "model = ut.read_model(modfile)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we need to loop through back-azimuth values, we will initialize empty `Stream` objects to store the traces from the output of the main routine. " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "trR = Stream(); trT = Stream()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now the main loop over back-azimuths where all calculations are done - self explanatory\n", "\n", "Remember that the `obs` boolean variable defaults to `False`, so if you want to change to `True`, either explicitely set them prior to this step or use the following call to `ut.run_plane()` with argument `obs=True`. Here we are not simulating OBS seismograms, so we don't need to specify anything.\n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# Loop over range of data \n", "for bb in baz:\n", " # Calculate the plane waves seismograms\n", " trxyz = ut.run_plane(model, slow, npts, dt, bb, wvtype=wvtype, obs=False)\n", "\n", " # Then the transfer functions in Z-R-T coordinate system\n", " tfs = ut.tf_from_xyz(trxyz, pvh=False)\n", "\n", " # Append to streams\n", " trR.append(tfs[0]); trT.append(tfs[1])\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The result of the previous loop is a set of transfer functions. To get receiver functions, we simply filter the `Stream` objects using some frequency corners" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "36 Trace(s) in Stream:\n", "\n", "... | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:29.990000Z | 100.0 Hz, 3000 samples\n", "...\n", "(34 other traces)\n", "...\n", "... | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:29.990000Z | 100.0 Hz, 3000 samples\n", "\n", "[Use \"print(Stream.__str__(extended=True))\" to print all Traces]" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Set frequency corners in Hz\n", "f1 = 0.01\n", "f2 = 1.0\n", "\n", "# Filter to get wave-like traces\n", "trR.filter('bandpass',freqmin=f1, freqmax=f2, corners=2, zerophase=True)\n", "trT.filter('bandpass',freqmin=f1, freqmax=f2, corners=2, zerophase=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now plot the result as wiggles, using the format displayed in the paper by Porter et al. (2011). We also need to define 'stacked traces' that represent the average of all recceiver functions to be displayed." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Stacking ALL traces in streams\n", "\n", "Plotting Wiggles by baz\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU8AAAEcCAYAAABH1zarAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydeXycVdn3v2cmk0nStOmS7qW0tGXpvqFs+rwgKIJSkUdFBFnFBRUUReHlUVTABRdEUdkfUBZRERCECigFZWmbtGmapUnTpmn2fZlmMuv1/nHue+bOZCZJk5nMTN75fT7nM/fc9332a86cc37nui4lImSQQQYZZHB0sCW7ABlkkEEG6YjM4JlBBhlkMAZkBs8MMsgggzEgM3hmkEEGGYwBmcEzgwwyyGAMyAyeGWSQQQZjQGbwnARQStUqpc42rm9RSj04ynj/q5S6PbGlyyCDyYnM4JkkGAOeWynlUko1GwNZ/njTFZE7ReSaeJQxg+TBkAszBC2y4lJKfSbZ5csgM3gmGx8VkXxgPbABuDnJ5ckgRSAi+WYA6jBkxQiPR76vlMqa+FIeHdKhjEeDzOCZAhCRZmArehBFKXW+UmqXUqpXKXVYKXWb9X2l1GVKqUNKqQ6l1P+NeHabUuoPlu9/Mma2PUqpN5RSqyagShkkGEqp25VSf1RKPamU6gMuVUqdqpR6RynVrZRqUkrdo5RyGO9nKaVEKfV5pdR+pVSXUuoeS3rHG/LRo5RqV0o9Ydx/UCn1o4i8X1RKfdW4XqSU+qtSqk0pdVApdd0IZTxFKVVsyHaLUuouy/unW8q/Wyn1/sS24viQGTxTAEqpRcCHgf3GrSPAZ4HpwPnAF5VSHzPeXQn8FrgMWADMAhYNk/xLwApgDlAMDJm1ZJC2uBB4AigA/gj4geuBQuB04Fzg8xFxzgM2oVc6l5p75cAdwIvADLQ83WvcfwK4WCmlAJRSs4CzgD8qpezAC8AOYCFwDvBNpdQHhinjr4C7RGQasBz4s5HuMcDzwHeBmcC3gWeM/FISmcEzuXjW+Ec+DLSiBQcReV1ESkUkKCJ7gCeB/zLi/Dfwgoi8ISIe4H+AYKwMRORhEekz3r0NWKeUKkhclTKYQPxbRP5myIlbRHaIyLsi4heRA8D9hOXGxA9FpEdEaoHXMVY7gA9YAswXkQER+Y9x/3XAAZxqfP8k8KaItACnANOMfXaviOwHHgIujlVGI58VSqlZhly+a7z3WeB5EdlqvPsyUIL+A0hJZAbP5OJjIjIV+D/AiegZA0qp9yql/mUshXqAL5jP0LPNw2YCInIE6IiWuFLKrpT6kVKqRinVC9QajwqjvZ9B2uGw9YtS6kRjSd1s9Pf3GdrXzZbrfsAkKW9ED5I7lVKlSqnLAUQkiJ4xftp47xLCq5djgcXGMrtbKdUN3ATMi1VG4EpgJbBPKbVdKXWeJa1PR6R1ClreUxKZwTMFICLbgP8FfmrcegK9hDlGRAqA3wHKeNYEHGPGVUrloZfu0XAJsAU4G71sWmJGi1/pM0giIk2i3QfsBZYby+LvMMq+FpEmEblGROYD1wH3K6WWGo+fBD5pfN8I/NW4fxioFpHpljBVRD4aq4wisk9ELkZvI/0M+ItSKsdI65GItKaIyF2kKDKDZ+rgbuAcpdR6YCrQKSIDSqn3oAdBE38GPqKUOkMplY2eXcTqx6mABz0zzQPuTFjpM0gFTAV6gCNKqZMYut8ZE0qpTyqlFhpfu9GDXgBARHYY6d4P/F1Eeo333ga8SqkblVI5xkpnjVJq0zD5XKaUKjRmtD1GPkHg98CFSqlzjHRylFJnKqUyM88MhoeItAGPofcwvwR839gP/Q7wtOW9MvTM4An0LLQLqI+R7GPAIaABKAfeSVT5M0gJ3AhcDvShZ6F/PIq47wV2KKWOAM8A14lIneX5k+gVzBPmDRHxowmo96C3hNqNfKcNk895QIUh2z8FPmXsl9aiyaX/AdrQx7NuJIXHKJUxhpxBBhlkcPRI2VE9gwwyyCCVkRk8M8gggwzGgMzgmUEGGWQwBmQGzwwyyCCDMSAzeGaQQQYZjAFpbeWksLBQlixZkuxijA9VVXDkCKxaBdnZyS5NWqKoqKhdRGbHO91JIV8ZjAvDyVZaD55Llixh586dyS7G2PHKK/DBD+rr974XHnggueVJUyilDiUi3bSXrwzGjeFkK7NsTyZutxhxf+wxqI911j2DDDJINWQGz2Rh2zZ4443wd68XfvSj2O9nkEEGKYXM4JkMiMC3vjX0/gMPQHX1xJcngwwyOGpkBs9k4Kc/hXffHXrf64UrrwSfL/FlENF7rnfeqT+DMU2CZpBBBlGQ1oSRzweNjXocgPBnLEQ+VyocTNhskJWlA0AgAH6/DsHg4DTMeGYaIvp9j0cT6JHjUVYWrGn7J7Y//Qk2boxeSLcb7riD2s9+h46uof9tubmQk6PTslkem+XOyQGnM/zMLFMwqJ873d2wfTv8/vdQXq5f+stfYNkyuOIKBk7/AFWHnEPGb7tdp5udHW6frCxwOPSzYDBcb5dLt1ck8vJgxgyYNk2nY7fr8g0M6Hhut/7/MNvarJcZ7PZwO/t8Oj+vN3ozZpBBopHWg6fDAQtS1GBVMAjd3dDSAu3t+ge/ciXY1p0FZ28fMf58D7S06YFoxgyYNw8KC+Nwmsk5XTP8JssfgRx0OTs69ACVkwP5+XrgVKO0AmoOiH19enBzOKCgQKcVCaX0oJqXN7bqTMQkPYMMoiGtB89Uhs0GM2fqMBY4nfr0UjKQlQVz5449vlJ6hpybG78yxYLDkfg8MsggGjJ7nhlkkEEGY0Bm8MwggwwyGAMyg2cGGWSQwRiQGTwzyCCDDMaAzOCZQQYZZDAGpLUPI6VUG9rBWQb/f+PYRFhVyshXBgwjW2k9eGaQQQYZJAuZZXsGGWSQwRiQGTwzyCCDDMaAzOCZQQYZZDAGZAbPDDLIIIMxIDN4ZpBBBhmMAWltGGTmzJly3HHHJbsY44bX6yV7Ejh/S1Y9EuUAbjLIV0a2xodJ6wBu6dKlk8JBl8fjwel0JrsY40ay6pEoB3CTQb4ysjU+TFoHcN5JYgm3rq4u2UWICyZLPUxMBvmaLH2SivVI68EzKyutJ84hzB2P8cwUwmSph4nJIF+TpU9SsR5pPXgGAoFkFyEu6O7uTnYR4oLJUg8Tk0G+JkufpGI90nrwtNnSuvga27aR88ILyS5FXJATzc9GGmMyyNdk6ZNUrEf6r0vSGfX1sGULLF2qfVZceWWyS5RBBhmMEmn91xpMd3e5X/kK9PQwMGMG3HgjtLUlu0TjwsDAQLKLEFekvXwxefokFeuRsMFTKZWjlNqulCpRSpUppb5n3L9NKdWglNpthPMscW5WSu1XSu1TSn1opDzsdnuiip94PPMMPPssANMPHICuLrj++iQXanyYPn16sosQV6S1fBmYLH2SivVI5MzTA5wlIuuA9cC5SqlTjGe/EJH1Rvg7gFJqJXAxsAo4F/iNUmpY6fVHcw6eDqipgc99LvS1xfTh/uST8NBD8cvH54MnnoDPfx7uuw/6++OXdhS0tLQkNP2JRtrKlwWTpU9SsR4J2/MUbSjUZXx1GGE446FbgKdExAMcVErtB94DvB0rgt2ejcsV6+nYYLNBjlOw+TzaAXl2Nt6AnXgd+ZtSWYS67FLweGDKFAAWv/VW6Jobb4RAAPelnyMQHKWjdAuUgimdh+G55+C3v4VDxhnfxx+HH/0IvvAF5FMXc6Tw2PhUyIJZsxbHvT+SiRE1WiJt4SqFiL4dDOpg3MZu17JFIDD4QWQYLyLKtPiYY8aVjDU5pUD5ffpPuK8Penv1tcMB06fD7Nm4VR6xDilkZUF2NtiUhBM26xxZ94gCLF60aEicKM0fPY3RtKvlXRHw+xnxN5/QPU+llF0ptRtoBV4RkXeNR19WSu1RSj2slJph3FsIHLZErzfuRaZ5rVJqp1JqZ3NzA11d7bS3N9Ha2kBfXxf19TV4vW4OHChHqSAVFcVkZUF5eRFKQVlZMV5vkMrKcjo63FRX19DS0kVjYwOdnU3097dTd/gQLr+fykOH8ItQUVFCXl44jb17ixCBqqpSvF4PdXXV9Pf30txcR0dHKy0trdTU1FFV1ctbb1Xz2mseXnutlMOHoQigooKibdugp4eSb36TyosuovLll3E1N1O7ezftH/84h+qa2batgVdf7eLNN2uorHRTWlqO2x1k795igkEoLS3iyBHYtauY1tYgZWXleL1uarxeuj7zGRr+8Q+aqqtpr62ldu9eXKWlVG7ZQmDRQsrKSigthVdeKeLpp2Hr1iLefht27CilpcXD/v3V9Pb20thYR3t7K+3trTQ01NHT08vBg9V4PB6qqkpxOKCioginE4qL/43NBpWVJfj9fg4cqKSvz8WhQ7W0tLTT0NBEfX0DHR1d1NbW0NfnpqqqHJ9P18nn03XyeqG0tBiPJ0hVVTn9/W4OH67B5eqira2BlpYmDh5sZ9u2Wn75y/iO1lb5amhooL29naamJhoaGujq6qKmpga32015eTlBEYp37QKlKCouBnRfgC53IBAud1NTA01NTbR3dVHb0IDL46Fy/378IpSUluo0iooAQp+lpaV4PB6qq3Vf1NXV0draSmtrK3V1dfT29lJdrfuitLRUxy0uDpdHKf7z1lv4/X4qKytxuVzU1tYOX6dgkOLiYpSC4uIibDbYvVvXqaKiHLffT017O11TptBQUEDTwoW0L1xIrQiuYJDa2krcbj+7d5fQ3KzT8Pm0jOTkQFlZKR6vl+qaGnpdLuoOH6a1rW1onbxeSvfuBZuNol27qNq/P1SnkhItX/v2VXLkiJavjo4odRoYoLyiIlQna9sWFxcTDAYpLy/H7XZTc+AAXd3dNDQ00NzcRE9PO62ttcPLykRYkldKTQf+CnwFaAPa0bPQHwDzReQqpdS9wNsi8gcjzkPA30XkL7HS3bx5s6S1+tzPfgbf+Ia+PvVU+M9/4jP7+P8MSqkiEdkc73TTXr4yGDeGk60JYdtFpBt4HThXRFpEJCAiQeAB9NIc9EzTusZYBDQOl25/gvfwEopDh+C73wWg6IYb4O234f77k1yo8cH8V58sSGv5MjBZ+iQV65FItn22MeNEKZULnA1UKqXmW167ENhrXD8PXKyUciqllgIrgO3D5ZGXlxf/gk8EAgG4/HI4cgSATXffre9/85uwf39883K74d13oacnvulGwaZNmxKex0QibeXLgsnSJ6lYj0TOPOcD/1JK7QF2oPc8XwB+opQqNe6fCXwNQETKgKeBcuBl4DoRGVY/Li1nBl6vPgy/bVvoVtENN+iLvj447zyorR1/Po2NcOutsGABnHIKzJ8PX/oSVFaOP+0YSMXZwXiQlvIVgYT2ic8Hra1MBEuYirKV1t4z163bLFu37oxJVFrvWVnQQCAcRDQL6nCA0wk5OZCXB/agD0SQLAden2JgQMuKlSw187DmHwzq944c0ZO9jg4w1XJXrIDNzlJU1b5wAi+8AKWlmmU32d2CAmqXn82OnTrfadNg1iyYMUOT8g6Hzs/v13l5PLpcDgfMng2z3PWwYwdRaU+lYM0a2mceT2WlHqfb2jQbPGcOzJun88rPN5hR2+D2M2Gz6WC36wDhtjHfjWzvyPayEq3R2F3zWWQ+Xq9u1/374Z134K67ErfnuWPHzuG3oSPo6Ej5Al3+rCzIsgsq4NcdZwqeWTmb7aj2u2MyzdF+zzHY6UHxItKODEqB3SYod79u/IYGaG7WlVRKC81JJ9GRNZfKSqir04ueqVNh7lz93z1rlv7uyJLoJw5iFcJaWCNEPrLWJWpbWNKPGs9M0GYjENAy5nLBnDmxZSutB89Vq1ZJWVlZsosxdrz9NpxxBiWf+xzrpk6Fu+5KdonGhZKSEtatWzfh+SaKMEp7+SJ5fRJvpKJspbVueyoaCxg1/H59eD0YZNWjj+rvl14KaSzoq1atSnYR4oq0li8Dk6VPUrEeaa3b7vF4kl2EsePOO/VyHdi/ZYsePK++Wq/D0xT74012JRlpLV8GJkufpGI90nrwTFvfLC+8AN//fujrojfe0BdFRfDlL8cvn5oabXxk7Vr47Gc1655ALFq0KKHpTzTSVr4smCx9kor1SOtl+8CAn6am0Wu3mWSGqTpnfjf37LOyNGmUmwvZDkEhBLHh8egNZJ8vvNdvIjLvQEC/19+vCaPOzjBxMHcubFZFqIcfhnPP1TddLtpXrCC/qUl/b2iABx/k4FlXU1YertD06Trk52tiyMzL69XBJIwKC2GWpxH1/HPw8sv6weLFuiC33w4bN8KWLXhOWEv1wSzq6gYTYCbpNG2aJs/s9jChY92kj+Q6ABob25kzJ38Q+RMI6DaLJNnMOLG08iLfNfMCXd+uLt1UiVQ/9/n8I2v4WQoshAkjk1wU0eV2ODQBp3yGIAWDgwXPaOgRVQ4ZSq4NIkoimJT2tjby8/OHxIlM3/yMJPvM34hSkJ0VRPX2wOHDUF4OJSWaGZo2TcvVf/0XjfnHU1IylKucOhWWLNGylZtjEGdmxmYHWxHBPLa3tpKfmxt63yxfrN9iqC0iKhqTNLMQRn6//v2OZH85rQmj9evXy+7du5NdjGERDIZPckydGvFj6OrSzHdBAYU/+Qmcf/6guP39mkm3DpjxRn+/Pq8fDGrBLiwcKsejRXt7O4WFhfEt4DDweGD3bjjllMQQRukgXyNhovsE9Im7vXu1eBcUaHO18+ePT3kuGfWASUwYpcPAb7PpP+aouOkmaGjAt2ABfOELUFGhR0oDeXk6JBJ5eXDSSfFJyzfB+7VOJ7z3vYlLPx3kayRMdJ+AniScemp800xGPUZCWu95pjW2bQuZnwtmZ2ur8rfemuRCjQ+TwXjwZMNk6ZNUrEdaD55p62OmrU0TOMbMJq+1Vd+/5x546aX45lVZCQ88oM+UJhiTQZ3RirSVLwsmS5+kYj3SWjrS0ljt/v1w9tl6o91A5/HH6wsR+MQntC3O8cDt1mmce65ek197LZx2GqxerS05lZYmhGnp7OyMe5rJRFrKVwTi3ider97QvP9++PSn9UmOM8/UR++qq+OblwWpKFtpvecJjkFsO8TelI6mcmZlK02DtSbj7nTq/UqTOfV69XhjMqgmrOyeyQBa2faOjrBNjuXLYfOK2ajXXgsncNVVLKiu1gPavHmhRA8cgJ07GaKeORzbnpVlqGc6fKjTT4fTT4/ddP39tHunsWuXJk37+3UeixbpYkyfrvdDs7IGs+3WeptksdnmTucCTIPfZhwr82zC1ESM1mfR2GCzb8z3fD7dpocO6bInCg6Hg2DQUk4sAmOFQUdLUAiKCqnNmqcMTLY9JwdsfkuHhfQ2s0LHFiQoQxpEqXCe5vMh7YSgbBHHIpRiwfz5oXhCdDbfZJuVTYXes6qZmidSnI4gNpdLC8jHP65DRGKHDmnvMv/8pyblp02DlSv13vTatXDssTC9QHQ7WI+6WNVTo9D9CwoLdYOabLuljCbMJGzK0hbWikZpv0HMvAiibPh8muRtbx9ePtKabV+9erXs3bt35BdTFY8+CldcQfmll7Kyuxv+9rdkl2hcKC8vZ+XKlROeb6LUM9Nevkhen8QbqShbab1sT2v1ufb2kCHkE594Qh+cf+aZJBdqfDjxxBOTXYS4Iq3ly8Bk6ZNUrEdaD55pbTLsi18MrQt2m1pF112nTXylKdL9TGQk0lq+DEyWPknFeqT14JmKDNyocMst8Oc/h75uvOcefdHcDB/7WHwMFweDmrm/4AK9mXnaafDrX+s8EoSNphfQSYK0lS8LEtYnZWVw222alLziCi1rCdwCTEXZSmvCyOXqPyr1zFgeAa0wN/cdDr0BbW5K+/1hAiDmpjv6uderCW9TPdNU81q+HDZRhJoxA37yE31zzx6KTjqJTSYzA/DKKxzY+N+DCKPCwrA9T1Pl2kpmmYTRzJlQKG2oV1/R+ovve58OoAv1hz9olY9ly5AlSymrL+CNNzSxVVAACxdqbZCZMzVhZLZDNPuJkYTRvn1FLF++aVCbmKSWlbiOpp5p5Qqita2VT/D5tBZLQ4PWK0gU+vv7Q4RRVNmKKGykeqbPF1ZtNAkje9BigFWpIYRRJBk0nL1N6zsh0sOq1qgURbt3s2nTphBZFMsWrTUfk6exyrzdDrnZAWwNhzWTefCgFsYPfEBH2rsXvF5qV3+EPzxpH/L/n5+vD86vXg1z5wh2rzusomqtv1nBCB3qotJSNq1bN0idM4htkH1dU8szJF8xGmo4e56CwudX9PbCSN6O05owSnsHXRUVWid4YEDb8jSdwWVwVMg4gEst9PfD1q16PHU6YfNmffDD6Ux2yY4eSSGMlFI5SqntSqkSpVSZUup7xv2ZSqlXlFLVxucMS5yblVL7lVL7lFIfGimPtN6TEtEqmQMDFH/1q9oZnOljPU1hunedLEhr+TKQjD7Jy4MLL4T/+R+tgXzWWeMfOFNRthK55+kBzhKRdcB64Fyl1CnAt4HXRGQF8JrxHaXUSuBiYBVwLvAbpZR9uAzSek/qvvvAMEW3/te/1n/Xn/98kgs1Pqxfvz7ZRYgr0lq+DEyWPknFeiRs8BQN0zOUwwgCbAEeNe4/CnzMuN4CPCUiHhE5COwn7JY4KgYGBuJe7gnBjh3w9a+HvlZecom+2LoV7rgjvnnt3q21QV55RW8+JhCVCXQulwykrXxZkLA+8Xjg8ce1mvEtt8THaeEwSEXZSihhZMwci4DlwL0i8q5Saq6INAGISJNSao7x+kLgHUv0euPeMOlnj4kwCsePTlSYmkYmURKpbRErHauGkds9VAMyKwuO792Juv56TdoA+P0sLSvTahgATz4J8+Zx+INX09s7NB+TxDHJGJOYMAmj6dNhtqMb9eYb8NRTg1Vwpk6FM86AE07QDPz8+fQs30R961CjvzabYdc0wglcZHnMDXoRyM9fSkvL4PsmgWa1gxrZV8P1WzRSw9TgiubfLp5wOLJDvM5o5ctKGJlki8kLOZ3gsFkYNFN1ygyG5kwkomlgRZVBIogWm42lixeHiJBI+bX2k1WGrZphpmzZ7ZCf48deWwP/+hf88Y/hY3VFRZptv/pqWrZcS3tvbCPSublwzCLB4T2i28D8sZmkkQmzDkahl86fr39UFgZR7Fn4A2qQc0GTNLLZjPY4WsJI2fB6Ndk7koZRQo8qiUhARNYDi4D3KKVWD/N6VD5zyEtKXauU2qmU2tne3oTD0Y5IE35/Aw5HF319NUyd6qa9vZzCwiD19cXMmQOHDxcxdy40NBQze3aQjo5ypk5109tbQ1ZWFz5fA4FAEzZbO319tdjtLhoaKsnL83PoUAmzZlnTKGLePGhrK2XGDA+9vdXk5PTi8dRht7eSk9NKTk4dhYW99PdX4/d7aG0tZcECKLYpeOstih57DIqKKLnoIg4vXUrlPffgevddal94gfYtW4Am+voa6OnRdZo5041IOYsXB2lpKWbBAmhtLeK44+DIkWJOPDGI36/rdKCzg673vY+Gu+6i6dVXad+2jdoXX8T16qtUXnMN/iuvpGTZMjj1VKoPlWK3w4EDRbS3Q0dHKYWFHuz2ambO7MXvr8Nma8Vma8Xnq2PKlF5crmpmzvTQ1lbKnDlQV1fEnDlw4MC7zJ0LTU0lTJ/up62tktxcFz5fLVOmtON0NuFwNDBtWhc+Xw2zZ7txucpZuDBIW1sxxxyj63TMMdDWVszChUF6e8uZOdON211DTk4XwWADDkcTs2a1I1JLV1d83d5a5au5uYmurnZaWppobGygq6uLmpoa3G435eXlBIPB0F5cUVERSkFJSTFZWUEOHSonO9tNZ2cNIl309DTQ3t5Ee1cXtS0tuJSisr4ev8NBSUUF2O0UFRejFBQX67T27i3F6/VQXV1Nb28vdXV1tLW10tbWSn19HS5XLzU11fh8HvbuLQ2x6zgcFO3dCw4HO3ftwh8IsG9fJQMDLhoaauntbaejo4nW1gZ6e7s4eLCGgQFdJ5Egu3cX43BARUURBQVQW1vM1KlB9tVU4V68mJpzzqHrmWdo+Mc/wvL117/iuuIK2noO8K9/+XnssRJWr4YHHyzivPPgnXeKcLvB5SolKF6qm5rotdmo6+qita+P1q4u6urr6XW5qK6pwRMIULpvH+TmUlRRQWNnJ0Xl5ZCTQ0llJX6l2Fe1D4/HRXNzLS5XO52d4TodOFCDe2CA8spKgkDx7t26fYx+2rWrGJEgFRXluN1uag4coKunh8bGBjo79TiQnV07vKxMFNuulPoucAT4HPB/jFnnfOB1ETlBKXUzgIj80Hh/K3CbiMQ0B7RhwwbZtWvXBJQ+QbjlFvjhD+lasYIZ2dmwa1firB5PALq6upgxY8bIL8YZiWLb016+SF6fiOijZMEgHHPMyLP2kZCKspVItn22Umq6cZ0LnA1UAs8DlxuvXQ6YJoSeBy5WSjmVUkuBFcD24fJIRRt/o0ZpKfz0pwD0z5mjDx0b39MVk4GdtiKt5ctAsvpEKW1kZvHi8Q+ckJqylcg9z/nAo8a+pw14WkReUEq9DTytlLoaqAM+ASAiZUqpp4FywA9cJyIJ3tVKEvr7tZthwzq2zSRyvvc9+OAHYdOmJBZu7JgM9i8nGxLaJ+++C3/6k3ZOdNVVCXV7kIqylbDBU0T2ABui3O8APhAjzh3AqOlmFY+/tIlGVxd88pOwZ0/oluPIEX3h8cBHPgIvvqgPz48HR47Aq6/C889r0mjePDjlFFi1SqsSBQJa7WPjxridXnak8ZZDNKSlfEUg7n3i98Nrr8EvfqFPh5j48Y/h7rvhoovim5+BVJSttFbP9HqDo7bnCdFV2iI1wsxnkeYFrZ43IxGp2maywb29eqy0qmdunN6FuvNObTwW4Je/xGWzUfjrX4f9F9ls1NToE00ejx7rIu15mky2aTvS1HSbOhVm57ux1+wH0zfScDh8mHLvcvbu1WlNnapVM6dNC6uCWj1oRqu7+aylxYXPVzjIRKPZHlYPmtY40dowVh5WG6H9/Vr1tb5++OqNB4FAEJ9vKCMdDZEmKK2qjSapHvKg6RJJLQUAACAASURBVPeFj2JYXYMqFZVtj5ZXJEJsu6mvaAiEq7eXwlmzCIoKPYpk261qjWYdTLkyDwZkZUHBFD/2A9VaQH7wAx0icHhfPzf/II99+wbfz86Gj34UPvxhOOlEIXugV2dgbRhTaKyUv1EXV2cnhebM1qLzGrQ7QmqwUetkMu6WzovJtqNVbD0e/ZttbBy+H9JaPXPjxo2SipoHo8ZLL8F55+GaP5/8886DBx9MdonGBZfLRb7Fgd1EIVGEUdrLF8npE58PHntMn5Tr6dE7UV/6kv4vHytSUbZSbyPhKOBN8KHvhKK/X0sUUP/+98PDD8Obbya5UONDfSKngUlAWsuXgWT0icMBV1+t9TK2b4fbbx/fwAmpKVtpPXg609HSgImbbgppZSx/7jm9bLjmGr1XmaZYvnx5sosQV6S1fBlIaJ80Nem9zwnwgJCKspXWg2faqs/dfz/ce2/oa9nlxsmtqqpBLHxcEAjojcFIG2GdnfogXhxRVlYW1/SSjbSVLwsS0idHjmg14hUrtJrxBRdo03QWp4bxRirKVlrvea5bt1m2bt0ZdSM/8nss1UIrYWR17haZZqTTuFh5WQmjvj49RvX06HgrVsB63w7Uo/8bTqi5WbsiXrMmnMgpp3DgjM+G7HlaCSOTxIm0NWoSE3l5MHN6kOzGWm13cceO8Gy2oEAzQt3d2sMVwLJlBC65jP29c2hp0cXKydFqdKZ6pqmqOhys7WLyH3a7vh/NcZ61TWMhkgS0EkYejybkWlvhvPMSs+e5adNmeeednYMIlVDBIwTMKh9WdV6rrcmQFqJNBj+0CmMUUmMkWQaLamYEYWR2RBBbSE6sZTL7KUQ4KRWSYa9Xt7NJGE2f4iNr727tmTWaN8vp02m89CaefX161PZUCjZsgA3rgjiPdOoM7Pawx0XTeKxZyQibnqHON4UpKwtycwnaHVHVf4cjjaJ2HiDKNogw2rRpGNkSkbQNJ510kqQ1GhpEZsyQnTfcIPLYY8kuzbixc+fOpOQL7JSMfEVFMvrkhRdEjj02/HeilMjll4vs2zf2NFNRttJ65pn2xmr/+7/hL3/R17Nmwb59+jODo0LGGHLqYWAA/vpX7X54yxZtiyYdMWnZ9lRU2Ro1nn02NHAW3XCD9oPxta8luVDjQ1FRUbKLEFektXwZSGif+HxaoSPyUCd66+fTn9a8aDwGzlSUrbQePNPWWO2BA5pZN7Dp7rv1xe9/r325TxTc7rgmtylN1UpjIW3ly4KE9clbb8HatVojbuVK+Pa3h9pgjCNSUbbSWsOov99Nuk0OcvfuQF12qWaUcnMBKL3sMtb8/vf6ha9+Ffx+3JdcPSyZMhKcWQHsbc163dTWpm8Gg9qr1c6d2iZjQwOceCJccw2y5WO4Z4zvMN6+faWccMKakV9ME/T3u/F6B3M5sfgG89PkNqx8kAmr5ovdDkqCQ1mOKPYmR8Ig529WL4VKUVpVxeo1awmKChF31ngWU6Ihwsi05RlZ/mzlI2vXDq3M8fTT+gVDhrnnHiguxvfIH/DNmMNwcNgCOLrbNBtlt+s0cnLC7KS1YY26lJaVseb44zXJZKo92WyaITVIIyuXFFm/Ye2xSphQCorC642uSTik3Uez52n4GVoAuIFaEUkJczPr12+SN9/U0/lYDRPN2LF53+qwL5LYi2R6R8MMm+9EY9tBq2cef7ylrMEgnH02nupqnP/8p6bjDTQ1aRuzHo9WlZw1Sxs6zs8Ps+2RhCRo2cvLG5u6usulx1mPR6dvas05HDrdSDXJyHbxej04HM5B6pnWtjXbyay/tb2Ha1frp1lnr1e3cU8PbN6cKLZ9k4x3uWj9IScDHo9n0HlVa3+NdILCZzj69Pl0/+fnj1yP9na48kp44YXB97OytO2Qa67R5hSsdo/HUo9EQkSf5GhqgpNOGgPbDhQAtwClwD7g38BO4DDwJ+DMWHEnKqxevTrO3NoE4+67RUCqLrxQ5H3vEwkGk12icaGqqiop+ZIgtj3t5UuS0yfBoMjvfidSWKiH6rPOEikrG1+aqShbw/33/NkYKN8nIieIyBkisllEjgF+BGwxzMolDVlZabzrcPiwdi8IzC0u1qqZDz+c5EKND3Pnzk12EeKKtJYvAwnvk9bWIb6xlNK+DBsadHjttbCXmbEiFWUr5uApIueIyO9FpDvKsyIRuUFEHkps8YZHINFObBKFQACuuEKv64Hu447T92+8UZNJE4U46253dw8RlbRG2sqXBQnrk74+rcA+bx4sXQovvzzklezs8eu0m0hF2RqRbVdKbYwSlimlkv63nIoGUkeE16s3f/75z9CtnK4ufdHTA+edp2eliYLPp49InX223qhfuBBuu21k+1ujQE5OzvjLl0JIS/mKQNz7JBjU8rNmjV4piWjZOf98bQEkQefGU1G2RjMA/gbYCOxBO2lbbVzPUkp9QUT+kcDyDQu/P+zA72hhZT5hsGabiWheCyOvI1UHTcLTJIy6usLqmcuXw3pVgTr/fC1soE3PuFzwzW+GmcaKCmr9i9hZpBgY0BqVVsLI4RhqY9Nqb9S0xWja0bTZtGrn/Kku7O++pSt67bU6mHj7bTj5ZDrzF9PUFHZUaNWcC7GyMdRWe3rCp5+iqVZGqgNGU2OM1rbWa7OOAcMJZZxPWw1CMBgmhGOx7ZHEmZVxt5JkkTDVNW1KBrNKFrZ9NGTTIPLR2qBmJoEAImH7olY1RqvTymhs+8CADiZhVFjgw7Fru07/Jz+JWp72PY18+1cLh3h+BU18XnIJnP5eP86uZt15DodmOE2WMysLQQ2Vg2BQvz8woH9cZsfn5mqd5SlTCNgcIbm32lEdVMdorHsE2z4woGW5pSV2uxvxht80B54CVlm+rwQeAY4Ddo8UP5FhzZo18doXTg527BCx2+XQmWeK3HJLskszbhw6dCgp+ZIgwijt5UuS0yeVlSKrVln/TkSuvFKkrW3saaaibI1mXXKiiIRMmohIObBBRIbdnFNKHaOU+pdSqkIpVaaUut64f5tSqkEptdsI51ni3KyU2q+U2qeU+tBIBbMf7XmHVII5+wsEmH7gANx1l3YCl8aYPj26QYh0RVrLl4Fk9MkJJ2h7NN/7HnzsY3o79OGHobBw7GmmomyNZvDcp5T6rVLqv4zwG6BKKeUEhrOd5gduFJGTgFOA65RSJuf2CxFZb4S/AxjPLgZWAecCvzGcx8XOIIEaDQnHj3+sXQ0DLRs36rXGVVclVEtjCHp7Q6RVPNAy4jonvZDW8mUgoX3y4otw2mmaOOroGPQoNxe+8x2t3/6hEadBIyMVZWs0g+cVwH7gBuBrwAHjng84M1YkEWkSkWLjug+oABYOk88W4CkR8YjIQSPP9wxXsOzs7FEUPwXx0kvw3e+Gvi5+7TV9sX271jBKNDo6tG+jwkK9mXrllXGxxbh48eI4FC51kLbyZUFC+qSvT3tB+MhH9F75ww/rk+8J9HGfirI1ImEkIm5jtvmCiERaAHCNJhOl1BK0J813gdOBLyulPos+dH+jiHShB9Z3LNHqGX6wpb9/IKR5eLQk30iE0XCaNIPrFv6MRRiZto3nz4f1/p2o++6Dc84xK0HVaaexZvdu/b22Fh58kENnX01ZeXhne9o0TRiZ9jxNjR+r465AQO+jd3YSss1pxbRpcNq0vdju+rEeQM8+Wz9oaYHrr4dPfAL/hz9KbcdU2tp0etnZg7XnYpFGwSDU1VUxf/6aEGdhqsqZPr6cTp2G0xnWMLGSW1ZNpFjtHEnMJHJy6HYPMDAQrnO0ekeWyRqsiHxmte+pGMyOjCRrQ+6ZjRGpJme3U7VvH6tWrx1CGFnJFNNeq0IQVEiDyySMvF7df7On+3DseAseeUQztR/+8ODC3HknXbfcxdtNS2K2aUEBnLzOS3Z7oxbW7OwQ4UNODpLlGOSozizrvn1VrD3xBE2umismv18L07RpMG0aMiUfj88WIrnEMPlpaslZ+zFWAwcCus59fWGt5lgYzVGlC4DdwMvG9/VKqedHimeJnw/8BbhBRHqB3wLLgPVAE/Az89Uo0YfUUCl1rVJqp1JqZ29vN0q14/c34fM14HB00ddXw9SpbtrbyyksDFJfX8ycOXD4cBFz5kB9fTGFhUFaW8vJzXXT0VGDSBdudwN+fxN2eztudy15eS46OiqZOdNPU1MJc+dCQ0MR8+aFP9vaSpk+3UNvbzU5Ob14PHXY7a04na3k5tZx/PG9HHtsNaee6mHevFI2bIBim4Jnn6XoBz+Ap56i5D3v4aS776byuutwPf00tb/5De0f+xjZzmZOOqmBVau6WLGihk2b3Nhs5Rx7bJCWlmLmztV1KiyEQ4eKmTkzSEdHOXPnuikoqGHDhi4KChpwuZpoaGgnGKzlhBNcVGVn4X/kEUp++EP4+991Of7+d4puvRUuuYSK+lqOOcZDYWE1a9f2Mm9eHXPntlJQ0IpSdUyb1suRI9XMnOmhvb2U+fOhqamIRYvA6fSydCkcOVLCsmV+bLZKjj3WRU5OLQUF7djtup+CwS5aWmqw2900NpYza1aQ5uZiFi2CtrYiFi+G9vZiFi0K4nKVM3u2G6+3hvz8Luz2BrKzm8jPb8dmq2XBglH9f48aVvnq6enG5Wqnvb2J5uYGuru7qKmpwe12U15ejkiQXbuKsdlg9+4isrKgtLSYrKwgNTXlBINu6utr6O/voqOjge7uJvr722lrqyUQcFFTU0kg4Kdkzx5QiiLD2VxxcRE2G5SVleLzeaipqcbl6uXw4Tra2lppa2vl8OE6+vp6qa6uxuP1UlpeDg4HRXv3Qm4uRRUV4HQSFCEY9HPgQCXBoIvOzlp8vnaOHGmiu7uB/v4uamtrGBhwU15RgUiQkpJicnPhwIEi5s+Hjo5i5s0LUl1bjfs976Hmf/6Hrscfp+GBB2h66CHaH3uM2t/8Btcjj9CSO8CePX7+9rcSzj8fXnmliPPPh23biqithSlTShGHUO3z0btwIXUOB602G61uN3VNTfT19XLwYDUiHqqqSsnJgb17i1i7dg1FpaUwcyYlPT34TzqJytxcXIsXU+v30+710tzaQkdHAyJd9PbqcaC5uZy8vCBlZcVkZcGuXVrdtnjXLoIilFdU4B4YoObAAbq6umhubqC3t4msrHamTq0dXlhiMUlmAIrQqpq7LPf2jBTPeM8BbAW+HuP5EmCvcX0zcLPl2Vbg1OHST3tjtV/6kghoY8jHHity5EiySzQupKLB2vGEtJcvSV6fvPhiWD3z7LNF6uvHl14qytZo9jz9ItIz8muDoZRSwENAhYj83HJ/vuW1C4G9xvXzwMVKKadSaimwAtg+XB5pbTLs7bfhd78DDJN0hw7pHfaJRH293mc92j2PGEhFs2HjQVrLl4GE9snAAPzhD1qGInDeeVBTAyUl8I9/aF2M8SAVZWs0g+depdQlgF0ptUIp9SvgrVHEOx24DDgr4ljST5RSpUqpPWjC6WsAoo9DPQ2Uo7cIrhORYfXj0tZYbWendvRmbPAV3XCDvv+LX8DWrRNThrvugmXL4L3vhVNOgcrKcSeZigZrx4O0lS8LEtYnO3fCunVw2WVahq6/fojjwmnTtMnPeFiUSknZijUlNQOQB9wB7EATPHcAOSPFm4iwadOmOE/SJwA1NSLr18fiGESmTBF55pnE5d/TI3LRRUPzzckRufVWkdbWxOWdIJCgZXtayleiUVkp8uUvi2RlDZWhM84QaW4eXTq9vVrWurtF/P6jK0Nvr0htrS5LTY1IR4dIIHD0dRkFhpOttPZhdNJJq+Tf/y4L9Z4J83q4f7xItt3KDlvtVlpZv0jEsjVpZdut3jOXL4f1x3ahfBaDHF/5CiWLF7Puiitg9uzQ7XrPbIp22ejv1+qZM2Zott2ixRZSdTTVHX0+ze4fOKANff/tb4NPIF1wAfzoO/3Myh7mbKfDweEjMzl4cDAZmp8fnXG3tl9NTQnz569jYEDHdbv1fYcjnMbUqeG0FII/oPB4wqyuyQYP14fWdlcKFi1KjD3PlStXya5dZUNV+0whITrjbr1vIpIIN+2lhth2i+BFk7dIWRuSqSl8ZjA6p6SigjVr14dUdU0ZH+TN07S9asQzDQJHsu2FBT4cfVG8ZlrQ7ZvCZV/Mj7aS5/TT4YYb4LSNA2QdsewEKqUFIicHcWSHymn1CltdXcKGtatRvT3RXXPbbDBtGl7lxOXS5Q4EdP1ycsIqxoPa29KPZhqmZ1aXSx9IWbkytmzFHDyVUn8jCtttQkQuiPVsohAPY7VJxR//CBdfjD8nh6wPfxieeSauyQeD2qDyvn36GN54zYKNBL/fnxQzbolyAJf28kVy+sTths99Dh5/XH/Py4Nf/lKfpR/rEj4VZWu4Pc+foo8RHURbkH/ACC7CJE9S4fF4kl2EsaO7W/8NA/u3bNGqGM89F9csbDY4+WS9vRp14Ny7F+69V5v/jgP2798fl3RSBWktXwaS0Se5udod12OP6a3QPXu0Bfnx7H2mpGzFWs+bAXhjNPeSETZs2BCvrY2JxyWXhBZcffPn6+v58yduz/GPfxTJztb5LlggEoejIH19fXEo2NGDBO15prV8GUhonzz8sLYAsmWLSHV14vKR1JSt0bDts5VSx5lfjGNEs4d5f8KQtrrHt90GTzwR+tq+xnCa1tSkNyYTbfj1iSe0bTDTGHJjI7z//VpzZBxoj9MMNlWQtvJlQUL6pLwcPvhBbYuhrEyvmNav19PNBCEVZWs0mwhfA15XSplWlJYA18Z+fSJho6trMKkTi9ix3jeJjkjCyOr8zXSwZu7DRyMxIm0DRhJGvb2awOnuDhNGG4JFqNxc+NGPdKSKCvJ9PvjhD8OJvfQSDWd8KkQY5edrwqigQBMv1s1vh0MHpfRY2Nys9zl37YreFpde0Muaxka4446hD1tb4bXXaFr5AUpLNdHldOp8zbxzcsJ5Wz0v6g3+fA4eDNe7q0tvvjuduvyzZ2t1+pkzweHr1wXOykJy8+gf0HX1eAbbnIxmF9SqWuhwjCwlY4VStjHZ8xycxmAex2pqMxZhZLVtGs2maeheNMLIWgi7nfwp+SG1TKuMRxJGIZuZNltIPdMk8kx7nrOm+XCU7IR//xs+8AEdrGhqom93DT9/bllUr7bZ2VodfvMqN/amep2406lZxClTQl4wI+3R6v7OJ+gPYuvr0XqT7e2a1bHbw8I1axberDz6+vQjv1+XOzdXB/N3o5ChhnuVQlAhO7FHjoQdN8ZErCmpNQBOYJ0RnKOJMxFh7dq1cZ+mTyj27xfJy5PGk08W+fnP45p0dbXIF74QXpnn5Yn85jeJ9THX2NiYuMSHAQlatqe9fEly+qSkROS44waP7BdeKDIek5ypKFsxl+1KqTMsA6xHREqM4DGeT1NKrR5hbM5gOHzxi9DfTzA7WzuDi4NlIxPLl8Nvf6uTfO45OHhQZzdo5vTOO5oWLSmJS57B0Ti7zmBCkYw+WbtWn6H/zGf09RNPaM8d4zGMlIqyNdyy/SKl1E/Q2j5FQBuQAyxHawYdC9yY8BIOg7T2MfO//6tdcAB5ra16nXDttdpcXRydfM+dq7dRh+CNN7ShxYEBeOoprdl02mnjymsyqDNakdbyZSDhfbJ3rxay2YNpkBkztOZmvJCKsjWc98yvAeejLR99AvgB8HW0zvl9IvJ+EdkxIaWMgbTd0N+9G77yldDXzuOP1xdbt8IPfpD4/Kur4cIL9cAJeoPo3HPh9dfHlWxn5/AHqNMNaStfFiSsT7Zvh1NP1Y7gFi2Cb32LqBudcUIqylZaaxitW7dRtm0rHqS5MZrqDEcYWZ+LDN5kj5aGeQ3h932+sHaNFXY7rOgrRn31K1r1yMjYnZ9PrtVz2nXX0fLxL9LZPfS/zekMkzambczsbO1IbMCjaGzUG/2xMNUxwKKvf1IPoJFwOOCmm2g9+xI6uobmbbOFN96zs8N2ErOydF07OtwMDOQO65StoADmFgawdxkb/llZUFBAIL+A3j5Fv8EjmW1n1faKJIucTpMMSMwh+Q0bNsq77xYPsuUZy56n+RlJdJn3rRpsECaM7HaLPU+LdlEkwRntcwhhZA1GQ/V7vDgcuVHzN9sxZKfUiGdq93g8YQLPbocZU/04KvbAk0/C3/8+1Pjq8uW4f3wPtWppzDbNzoalc/uxtTTpxHNyBtnzDNqyQmSR1xtWJhJxU1iYS1Z/ryY229o0M6mUVr2bMwdmzcKXMzVk7tPjGexfzvzdKIkwFhqFMHK7tXguWTIGDaN0wOrVq2Xv3pQ4rz8EIvqPuLMzzLYvWaKNJYTg9cIpp1C+di0rr7tOn2g3EAjofcr+fi1X06frgWckJQufT6tm/vKX+ty9FdOnw6uvwkgGagIB2L9fq6c5ndrY/KxZmvUfbkehvLyclStX0tenTz81NWn5zsnRhqAXL9Z1iDcSpWGUyvI1Wph9MlEIBuHXv4abbhr8J75kCdx6K3zqU1qOjhYTXQ8Tw8pWLCYpHULaG2647TYRkIDNJrJ6tYjXG9fkt24VWbpU/8UuWyayZ0/EC2++qY0tvvpqXPILJMg4w0ggYxgkJpLVJ3v2iJx8ssiMGSLf/a5If//40ktF2Ur6ADiekNbGaisrRZxOEZCir35Vd8Udd8Q9G49H5J13RNzuiAfl5SJTp+p8c3NF3n573HkVFRWNO42xIFGDZ1rLl4Fk9Um8kYqyNaplu1LqNPTh+NCiUUQei8OseFzYvHmz7Ny5M9nFOHoMDMAZZ+jT7FZkZ8Obb8J7hvV7N3643XqLwOrqePZs+M9/YMWKxOadACRq2Z628jVRqK+HF16AefPg/PMTq7GQJIzVMIgZ+fdoIyFnACcbIe6COhakpbHa3l646KJBA2fIGLLXG/ZImEh8+ctDfcS3tWkVzX//e8zJprsFokikpXxFIO59IqLPB19xBRx3nD48fOGF+mDxww8PZUnjhKj1cLsTyvCPhBFnnkqpCmCljGaKOsHYsGGzvP76ziHe9kaClW2P9MpoPjfV6kybgpHPrMGEWNh2Uz3TShgtWwabZx5AtVvc8t13nxbGX/9aM0MAdjtN8zawc5ed/n7NFFrVM027miHm3aapVJ/KZv9+eP55ePbZ6J4or77cz7WbhvlB2Ww0z9/Av97MorlZpz9njj7KN3Om3uw3GffcXCNvrxfsdtyBbFpb4fBh7YKhpkZr0eXmasLgxBPhhBNg0WwPtobDWn/TZoPCQmTOXI4Ecujt1RNzq5dHq+1QU63Q9OqZlyvYs2wJMkm3Wd55Z+cgpt0qLybM60jZM9+1emY13xnEdMeLbbcWwmDbg6IGMe3mY5NtH/Q7EAl5zzRZZ1NF0m6HaVMCZB2q0f0WA+65S/jFE3OHGAhzOuHjH9dzg2UL3Po34PPpjszLCwlVAHtIPdPKtjudWvYd3iNaqDo6wg8djpDerz8nn95ezZR7PLqNTbY9pJ4ZybYbDWCe2vH5dFy3GxYsGAdhBPwJmD/Se8kIab8n9eqr4T3Pq6+Oa9Lbt4usWDH413XOOXoPNISiIpH3vz9uhFEq7kuNJ6S9fEny+uTll0XWrtVO4G64QaShYXzppaJsDTdo/g3tlO1fQBfam+XzZogVzxL/GCNuBVAGXG/cnwm8AlQbnzMscW4G9gP7gA+NlEdas6Fud2h0C9hsIkqJbNsW1yy6u0U+9Sndyx//eIRzztZWbQIPRPLzRfbuHXd+qciIjiektXwZSFafxBupKFvDDX7/NVyIFc8Sfz6w0bieClQBK4GfAN827n8b+LFxvRIoQRshWQrUAPbh8li1atUENF+CcOONoSlh2aWX6uvly0USYLewpyfKzQsvHDwtPf54ka6uceVTVlY2rvhjRaIGz7SWLwPJ6pN4IxVlazj1zG0isg04z7y23osVzxK/SUSKjes+Ywa6ENgCPGq89ijwMeN6C/CUaCMkB40Z6LC0c3Z29kjFSE089hj8POSNmaUvv6wv9u/Xp4jjbMF80MF80ErHkSfoq6q0Enxv75jzWbo0tmZJOiJt5cuChPVJf792G/ONb8Dtt2sbiAlEKsrWaOx5ngN8K+Leh6Pciwml1BJgA/AuMFdEmkAPsEqpOcZrC4F3LNHqjXsx4fX66Osj6ob48OVhiMpdpHk/GOyEykw3GtlkQiIIo56eKISRfRfqnXe0ERCA1lYaFy1imdfiFO6RR2i54HPsKLZz5MhQwsgkbEzSxubV3q4Cjhyqauz85z+63NFwzMIg55eVwec/H/2FBx6g6b+/wivbskNE5rRpsGCB3pOfNi1chlxnEOXu1wxPVhb1re3k5i2npkZrf0aWYf587a128UwXtgP7oaVFN+K8eXDMMXjyZtDVpW2k+HyD9/NNcsNu15v+Tqcmr7Ilca4yfD4fgUB0ctCKkWTOlAtTPq11irTlab5vlTcTg4ii4T7Nl202GhoaWbJk2ZDfxhA5tmRolte0A2qY+STXGcTe0gjFxbBtm+4o0OzMfffBqlV4L/o0f99eSHPz0HZYuFCfkJtbMKAduZneAU19X4cDsdlD+ZoBoLm5kRUrlpHlc+s/+O5unb+IFkjD6GzAmUd/v/79+Xy6jU0HcKbzQiXBwbrYRmOIhMk9k7AaoWNjLru/CJQC/cAeSzgIPB4rXpR08tFWmT5ufO+OeN5lfN4LXGq5/xBwUZT0rkW7QN45f/58aWtrk8bGRqmvr5fOzk7Zv3+/9Pf3S1lZmQQCgdBG807DzURRUZEEAgEpKyuT/v5+2b9/v3R2dkp9fb00NjZKW1ubHDx4UPr6+qSiokJ8Pp/s3r17UBrm5549e2RgYECqqqqkp6dHDh06JC0tLdLS0iKHDh2Snp4eqaqqkoGBAdljqPcMSqO5WXbfcIO0rVolFU89JX19fXLw4MFx12nHjv1yzTWdcvrp9XLyyY2yiiHWywAAIABJREFUcmWbfPCDB2XRoj555RVLnQ4flp333COyZ8+QOlVWVkltbY9UVY2+Tm+99ZaIiOzevVtcLp+8/nqFvPhin2zdelD27UtcPxHHZfukki8Reffdd8Xn80lFRUXc5GukOvX3++TPf94t55wj8oMf7JSbbhL5xz92SjA49jp1dnaGyrN79+4Jq9NwsjWcEBWgD8Y/iTY/Z4aZRyGIDjTR9HXLvX0Y7D16X3SfcX0zcLPlva3AqcOln/bGaj/9aRGQ+tNPF5k1S6StLW5JB4Mi3/++Od/Q4aGHLC94vSIbN+oHCxaItLSMO8/6+vpxpzEWxHPwlMkkX5K8Pok3UlG2htvz7BGRWvTyXCwhXyk1ollTpZQyZo8VIvJzy6PngcuN68uB5yz3L1ZKOQ0/SSuAKN6fJwleeklbpwFsXq8+t/bVr8YteaW0feWnntIKS488ol3OhPDDH+rlF2grHp/97Oj2PIbBZLB/OdkwWfokFesxmj3PF9GDpkIbQ16Knj2uGiHe6cBlQKlSardx7xbgR8DTSqmrgTq0rVBEpEwp9TRQDviB60RkWHUFFUejwROKw4fhyitDXx3m3tGTT8KZZ2rr7nHCpz6lwyBUVg71YbR1K9x8c9i30hjgmGTqeWkrXxZMlj5JxXqMOHiKyBrrd6XURiAG2zAo3r/RA240fCDaTRG5A4jimSw6AoFgyJ5vusBZsRv1mUv0hrfTCYBryRIKa2r0CzfcAEeO4P38Vwgq+8jpZQtqwA1+P8GcPLzBkf8Pc667Tk9NjfxDuPtuCATwfud2gg5n9MgGlIJsmx/l6tO781lZ9PW5yM8vHDH/7IFebPsqtG60UppJWLIE34w5BIKjH7CysiDLc2TU7x8tAoEgweDwZJGJaJN2q4aRldMZoqFmIWtGSitmhpHfjcT7+lzMnFkY6/FQEsqi5WS1k6sUZNkF2xHD3mBVlVbxPXxYMzLHHw/r1xNYtRafI7bV95DceAZ0wqbKmOFN0NRwsjrLA+judjFjRiF28WuCcmAgfColO1uXISeHoN0RIhutsJKOIZLO2rAR2l1W+74x6yJjWKoppYpFZONRR4wzNmzYKNu3F4fYdmtlzQ6PBgvBFmIbrc1gfjdZt0gmPpKpN+OanW4aUzXZ9q6uMNu+Zk2EwH7oQ7gOHSL/b3/TAmigs1Ovql2usD3P6dMHs+05OWGDzGZyBw/Co4/CXXcR1SjxTTfBj39sfAkE4E9/go9+NKwailZzf+kl7f8oL0+PbYsWadshZhny8oa2b1+fC683n8pKbSx/1y44dEi/u2qV3j44+WTN3Efrm2CQkGqd1zvIpu8g473mb8XpNPsgMYZBNm7cKMXm1kaawuVykT8WA5pjhIhWy33nHW2+YWBAq8CvXatVcxcuPDr7IeZvr79/9PUQCXsLNcfmsaz69dg+PvXMr1vCN4AngK0jxZuIkPaHmO+9VwSk4lOfEjn99Li6tiwtFVm8eDBhdMopEeqZX/yifnD22SJx0OCoqKgYdxpjAZlD8jGRrD6JN1JRtkYzHk+1BCd6D3TL0Y/j8YczctmZTmhshFtuAWD5c89pc3D33x+35Fev1gaSTjhBf1+6VJ9pDp373rpVu9cEbV4+mh/3o8Ty5cvHnUYqIa3ly8Bk6ZNUrMeIg6eIfM8S7hCRx0UkJXYaB9Jtw9NEMKip754eAMouNw4ffPObei8pTjjmGHj3XXjwQf05f77xwO2GL31p8Mvf+x7885/jyq8s0sxdmiNt5cuCydInqViP0djz3KyU+qtSqlgptccME1G4kZCbm5vsIhw9/H5tA3Hr1tCtdffdpy/6+uC88/SmUZxQUABXXx3hGfa22+DAgcEvBgLaZth//jPmvNatW3d0EUSisyMpgrSUrwgcdZ+kKFKxHqM5qvQ48E20tlFKeZ4/cqR/kHqmubkcSf5EYrSEkWlXMIom1xD1TCtD6fUOVc8MBjVhdHLePmxnnqmPJAG89hpFS5aw6dhjw97dKipoyV3C9qKh6pn5+SFikZwcyM0x2Hafj2DuFCqqs3j88aFjo4nNm+EbJ58cOmM6BM3NNBzw8PxWZ4gwWrRIz2LnzNGEUX6+Jo2cNp9meAy2fefheuYvPJmKCk12lZaGTS7Ong2nnALvfS8cV9iLraIszLYvWABLl+KbNY+OTkVf31DCyAxWwig/H7IGXMMLyTgwWYwhbxrJ418aICXrEWsz1AzAv0d6J1kh7U2GFReLZGXpcfdb34pbsh6PTi7SWm5ennadFMLrr2tW6Ykn4pa3FQMDIgcOiDQ1xZULGwIyJukySBCGk63REEbfVUo9qJT6tFLq42ZI6Ig+SqT1zCAY1MZB/H7thuNnP4M98dkNyc7WZ93vvXfw/YceChNItLXBJz+pzyNdeWVY22gciHSV4HRqomrevJHPSaYi0lq+DEwW1yipWI/RDJ5XAuuBc4GPGuEjiSzUaJGXF/swbsrj5z8Hw7nYprvv1nsEV10VXufGAV/6EjzxhHZ/cf/9cPHFlofXXw+trfra49FqSAaBNVak3LJqnEhr+TIwWfokJesRa0pqBqB0pHeSFVauXBnfOfpE4dVXRRyO0Hp6z1VXhdfWV1+d2DWuiMiLLw5d05t+Oob4KB499gxxDD8xIEHL9rSVLwuS1SfxRirK1mgIo3eUUitFpDyho/gY4HBoh2EmWWRVg4sFc/kYqSUEQ7WSIgkjM56VNDJhjkCR9jxNsgi0A7XNjhJsv/oVnHWWvtnfz/FVVfChD+nv9fXwq1/RfNF17Cmzh+KaJgutGka5uZDjFG3P0+8nkJ3Lvpos6upi13/qVDj9r38N5xeJO++k8apb2VM52BBwdra252mSVlOmgNNuqGe63WC3s3zJUhoatC3PaKd8lNKk2XGz+warZxqEkbdgNh2dapCGkdnOVtLIShg5fIlbWuc4ncOrqkWDKYCRcWJ4EYyU12jqkrGyiLyOhFJw/PHHD/6LjMg/WhohlVHzR2UY9AzasuhzKRobtVZmpPpiQQGsXR1kSl1FbBWjRYvwZU9hYEDHN0lAs29tNhikn2lkcvySJUhQ8AcUAwNa5Ex7m1aNs2yHoPyGMU6ripERhKFtbm2SQT/kETyBjtZ75jK0HU8PWl9dRGTtsBEnAGvWrJHS0tJkF2NYBAL6BFIwqBnzQb+Fvj5Yt47q9etZ8fnPDxnQTA9+eXmWw+0jQAT+9S/tXbiiYujzP/wBPvMZy8vbtsHppw/RmTt0SP9ATPXM2bNHVnGrrq5mxYoV+P2wb19YPTM3Vx/aP/lk3QbxRqLUM9NBvkaC2Sfxhtut5cuq2r5kSeL2thNVj5EwrGzFmpKagcG2PENhpHgTEdavXx/PGfrE4ytfEQHpOfZYzXq7XHFLuq9P5EMfGrwqv+aaiJe+8AX94IIL4qKe2RPVWVLiQYKW7WkvX5K8Pok3UlG2Ys4llFKm55u+GCHpCIwwrU5pbN8eosO7jztOs9633hq35PPz4bnnwiTRRz8awb4//TT87nf6+vnn9cH5caK7u3vcaaQS0lq+DEyWPknFegy3EHvC+CxCuyUosoSdCS7XqJCKBlJHhe5uvXY29nNyurr0/V/+El58MW7ZOJ2abd+/H5591rL0b2kZqp55++16QB0HcnJy/l97bx7X1nXm/7+PFsAgY8BgDMYYbMDYGIORk0yatkmnzTJppm2WmSZd0jRJky5pm+Q732k7nc5kftM2XZO0nS5JM2kybSbd26TNnnzTpEnsxAhDwDbGYAzGwhjMvgi0nN8fRxddCQGykZAuL71fr/vSwtVZdB4d7j3PeT7Pkj6faBjWvnTEdUxGRpSGQxRIRNtaSEn+Cv9jqZRys/9ROzYvXxNXGN3dcPHFakYLRUq193KJk5gezUkTNA/cfLNSrg+t+yMfUVeji3ndokFfn9rX2tys9pwmWTm8+ipcd51aKN+wQe2Vu+sutTVvvsyEsSDGdw6LetuFEC8C35VSPqV77wEp5S0xbVkEeL2+ZcmeqR+DJWfP3GxDPPFE4EO3347L41EinLoA9D6nN/Lsme5p5W33Z8987DGVdiOcnX76kz6+soh6U+9xD8+8aOXECVV+QUFk2TMnx6fomVGh+QcOQEuLuvBYtUo5Wnft8mfPXDuByecL9NfjgeFhpldlJVT2TJ9mWPotGaExueGea+i3cYQqC/s7I3V64UGedn15YVJoan8OfdRjMoFrSgllB3ncw/wAgpz7miGHpM/0pqbT22fC4VB6r88+G9hVsWmTupm67ho3uVu2qH3M99wzt1EjI7gy1jISPnkmZpMkXPpM18gIMjuHaa+FkRH12xofV0212ZRdrlnjt8nJCfUD1JSstR9LSgrSbAn79QrhF0k+k/SZ8y2GagdwFHgZ+Hfdew2LfW45DsMv6P/61wGH0Qc+ENWi33xTyqKiYIfReeepkMlZ+vqk/NrXpDxxIip1JuKi/lIOw9uXjO2YjI1J+corUr71Vuy3JieibUWyqDOMSpuRL4T4kxBiTQSfQQjxkBDilBCiRffeXUKIE0KIRv9xue5vXxJCtAshDgsh5tmEGIxnOW8Bos3IiIryAfrq6tSi5B//GLXizzkH9uxRV3oA27crv9CsROXICFx0EXz5y/D2t6vb6CXSF4UyEglD25efWI6JzQbveEdIdoQYkYi2FcnkKaSUHinlp4HfAa8C6yL43MOokM5Q7pVS1vqPpwCEENuBa1FJ5S4DfiTE4gl8UiLd/JiIfPaz0NsLQPGLL6r3PvWpQMhkFCgqUnuVn3xS6Xmu04/aTTcFNoJ2dsL73qful5dAcfGiSVUNhaHty09cx6SjA15/ffFkQBGQiLYVyeT5E+2JlPJh4AbgucU+JKV8BRiMsB3vB34ppZyWUnYC7cC5i33IsGK1X/0q/Pznsy/brrlGPTl5Uk1imvc9CqSlKYnQoPQvP/wh/O53wSe++Sa8971qkfYsaYuikHMiYFj70hGXMRkcVLs5KipUAMbWrWoRfvrs16cT0bYiyZ55P4AQYh0q9XA/cNcS6rxNCHE9arvT/5FSDgEbgL26c3r8781BCHELcAtAYWEhXV0DuN1uPB4fKSnpjI8Pkp1dyMmTnRQXV9LR0Uh5eR1HjjgoL7dz5EgD5eW1dHe3UlhYyuCgk8zMHCYnJzGZTFgsVqamxsnKyqWvr4eCgjI6Ow9QWlpDe7uDsjI7HR2qrGPHmtm4sYK+vm6ys/MZHx/GYknD64XxcRceTxZDQ31MTBRjNrexaVM1DD3PbrMZx8MPY3c6aUpLo+qpp2h94AGKRkcZSE/H9swz9O98J0c6wOVKJy1tkDVrCoFOiooq6e1tZNs21afdu+281eSgdvt2DrQfBTazb58TpzOHtLRJfD4THo+VVavGGR3Npbqqm/fNzHDgoYeoOXkSR2Ehdqdz9rH5+efJOvcK9tY76evLx2YbJjMzjawsSE93kZ2dxdRUHxs3FtPT00ZN1TaltbhrF1NTM3R3w6FDTfT1VTE52U5fXxGZmQPYbDY2bnRTUOCjKN/CcNcxCn0+OqemqCwtpXFggIpt59DY6KCw0E53dwMFBbX097eydm0pY2NOMjJy8HgmsVhMpKVZ8XjGKVi/dgmmuLh9DfT14fZ68fl8pGdkMDg4SGFhIZ2dnVRWVtLY2Ejdrl3qO6iro2H/fmp37qS1rY3SkhKcvb3kZGUxOTGBSQisZjPjk5Pk5ubSc/IkW8rKOXDgADU1NTgcajwd9fXY6+pobmmhoryc7uPHyc/PZ3h4mDS/QPPUlIs1a7Lo6+ujqKiYjo42tm2rprHRQW2tncZGpX/p83rxuFy0d3RQtH49A4OD2NLTcft8+KQk3WZjcGiIgoJAn5qaGqmrrsaxfz/2igoaDh2idutWDvY48bEFh8PJiRNz7cvlyuXdFx6j7sAeDpxzDjUbNwbsy2TC/vDDNFdVsXHHObS1dZOeno/bPYzNlkZqKvh8LnKy19DX20txfj5tHR1Ub9qEo60N+44d1O/bR+W2c2hsbCIlpYrTp9uxWouwWgfIzbWxapWbFIuHDOljcGCAwvR0OoeHqdy8mcauLjVOb73l/34aqK6upa2tldLSUk6edJKTna3GSUqsUjK+mFDOfIuh2oFSUToCTKBCNL1Ay2Kf83+2RH8ukA+YUVe8XwMe8r//Q+AjuvP+G7h6sfK3bdsWgyXiZaSjQ8qMDFl/++1S3ndf1Ip1uaT81rekTE0NdhgVF0vpdEatmjnU19cHvfZ6pRwZCXFSxQBi5DAyvH3JuWMSDTweKX//e+WABCnT0lSw2vHjUa9qllj0IxIWsq1IJsAmYC2w3//6XcADi31Ohpk85/sb8CXgS7q/PQucv1j5hherveSSwMyWkSFlV1dUi6+vl7KkRBW/aZOUbW0hJ+zdK+Utt0jZ1BTVepebWE2ehrevZeDEiahGFSccC9lWJGuebinlacAkhDBJKV9C6XueMUKIAt3LKwHNE/8EcK0QIlUIUQqUA28uVp6hxWofeQSeU0vHjttvV86aW2+NahV2u9qD/vjj0NQEQboKr76qUoE88IDytu/Zs+T6ElGwdikY2r78xHpMCgvVvt+wHDsGf/1rVIIuEtK25ptVtQN4AbABPwAeA74HvB7B5x4DegE3ag3zJuDnqFxIb6EmzALd+V8GOoDDwN8tVr408pVBY6OUNlvwPbV2/Od/xr7+I0ekXLs2uN7MTClffjn2dccAkleeiYXPJ+V3viNlSoqyrbo6KffsiXerzoqFbCsSSboMYAq1TvlhYA3wqFRXo3Fl27YquWfPgaDgjUj+yekDLLSwRX0COe3vUgbeDy1XL4mo7XHT2uF2M6tXqMdshi0jDYjP3hYIj/R6abrySmq0qCMh4LOf5dTVn2JoZO6NgdWqPOhWq4rK0KIzTELimhb09i4cGGGzuNhwxz8qzbhQUlLgC1+g/+IPMTg8t24hdLqJKYE2WCzqe9q/v4l162oWdKquXg35uV7MQwMqRMRigTVr8NrWMDommJwMSDGG+561hHAWSyDCatWq2EjSVW3fLg80Nc0fUgYLh/roI4z0YXCgyrNakSbz7FtzNCWlTktU96ivKtyhYTZDS0sTtTuqgkVvtfo1IU2TKag6gVRGPD2tDo8HLBbcGVl0nbAsGPW4Nm2C3G9/AZ5/PvgPJhNccQVcfz0TpTvoPSnwelUTMjLUkZYGFrNfj1OL8PFnVmjq6mLnrjpGJ8wMDMyNnktJgbVrYfUqD2JslNksghaLCtFLT4e0NKTFOuf3rE/siNerPjc1BePjiE2b5rWtRSfPoJOFuEJK+eeIPxBj7Ha7TMjLedTgjI2pXUdaeGZJico+OcvMDJx/Pp72diwvvQR1dbN/8nhUBszJyeDwzMW2Hs7MqOzB3/ueul3Xs2YNvPCCyqC5EB6PEjQ+fVpNUGvXKiHnzMyFNT09Hg9ms4XRURWW2durkmumpakQz+JiY+l5JrJ9RYrH48FiiUTzPDKGh+E3v4G771bbg/V88pPwzW8qO1kM7f+J2RzZBnt9P3w+9bsaHVWfzcpSth2LjfpL0vPUHyRIWKZ2VFVVRefaPF78x39ICfLQBz8o5c6dUrrdUS3+6aeVowikLC2NvV/o0KFDsa1gHojRbbvh7UvGbkympqS8+24pN2+W8qKLpPzrX2NSzSyJaFtnqrmVUDkQDR0BcvgwfP3rABS98opSGAonpLAELrtMVfP663DwoMqEEMRDDynBz1dfjUp9RUVFUSknUTC0ffmJ1ZikpcEXv6iCiF56Sfkc5zAxoZxGUSARbetMJ89bAYQQqYuduBwYNvbY5VISNP7FwYHqavX+v/1bVFIA60lNhfPPV8YexB13qBDNX/1Ked3vv3/JdQ0MDCy5jETCsPalI25j8tJLUFamck9fey0ssR2JaFuLTp5CiIe051LKN4UQNuCpBT6ybBhSrHZ0FK6+WiXH8mPr6VFPpqdVLGV9DLWmpVQJju67L/Cex6MWrG68cUlGbguKATU+hrSvEOIyJt/9rtKsPXlSvf7Vr6CqSv2Dnpo6qyIT0bYiWUk+IYT4sZTyU0KIbOBJ4KcxbldEeL2SwcGAXOIZ+L7meNs1GT/937X35/O2a4+hso1ut7KR0VEV5jsyEtDzrMs8hfjXfw2k3PjpT3GPj8P3v688Q/6CThz3Ud9gYmpKeSL1ep6pqcrTrR0Wi2rv9LSyV4cD/vAH5mTRTE2FH9w9QfWHP6zLAhfCiRM4Z3JpblbtT0lR9WpZMzVPvz6bpcmk+t3f72Z0VH1Oc5S5XKrenBylB5qXp55b3ZOz3lBfWjqTLtOsp12Tnwz3XWuPmsc9JG9dVJE+n2qQ1sn5PBLh3N/6BkOwt1tLG2m1zupLaszxuIeWM4+3Xb/bRPu72QwzM+5gYVrtRH1WSbM5aKeJQKp+T0+rAXS7lbc9cy1vNlh46CG1DBRKdjb84O5xtrztbWp/Zzja25koqeJEr2nWNlavDmjEWi0S4fWnrXW7VbuFwD0+ji97LcMjgv5+9T9+fFz1MTs7YFcZKW6V0VUThbVYtC0ZkJo6x9se+jsWPr+3fWJC/XAXYr7FUP0BfBMlELKPCMIml+vYuXNnlJeHl5kXXpASZM8FF0h5661RLdrtlvLRR6XcsUPObrWLdYRbT09PbCuYB2LkMDK8fcnYjInHI+V//ZcKitOm6/PPl7KzM+pVzZKItjXvlacQ4irdyzeBr/gfpRDiKinl7xeelmOPoW+rXC4lQQeknzqlIn2uvx7e9raoFG+xwIc+pLIhjIzEbiuHnnTtynmFYGj78hOLMTGb4TOfUUkF779fLW1+9KPK5mJFItrWQtbx97rjCmA/YNW9jjuGXtD/8pfVZkpgsKJC/QO/6Sa1sTOKaPvg5kycQ0Pwuc+pX8Azz0SlrsElyNklIoa2Lz+xHJPiYvja1+DjHw8zcQ4MwDXXKIfRImlfIiEhbWu+S1IjHLt27YrmFfry8bOfSd3ylJzMzQ28vvJKKWdmYlu/wxHYAKodN94o5eDgkoqdnJyMTvvOEGJ0225Y+9IRlzF5662AIo12XHHFkmSXEtG2IkkA9wjweSnlsP91Nioh3I2xndYXx+Wa4fTpuWFpMmQxOBzhwjO16Dl9KOBiYZ+LOYyGhgIOo7IyqPXUI/buhU98Qn2ov5/ObdvYrvdy/8//0HPxx3HsV04Um00thusdRlpknXYIoeodGlKRSUePhm/r1W/rJeeRR+CSS+ae8M1vwsc/zvH0rXR2qj6kpKg6NWdRSkqw/0T//R071sm6ddtxudRnp6YCSb60MlavDpQlkHi8YtYvoTmLwiXjC+1H6HcfC2a08ER9hxdKAKc/QpFyrsPG7/HTnz4nPHP2TQjnONIfen8UqGZ3dHSyY3tlIJma1oYQh5EUgZtQIf2OMpcrMDBWKzNr8thbb6G1df7v7OJzhin9xcPK2x7Kt78NV1+Nq+5t9Jy0MDUVSOQ364y0SOW00drqT8DX2dlJxfads07YwUHVNKtV3Vnl5KjHVDGjPEkuF7Pxn/qYYn/SvdAhmv2KfV415uPjczPMzh3Thf/74peiW+y9eByGF244eVLK3FzpNZmkfOyxqBbtdEr5k59Ieeml6iLgqqukfOONqFYxB6/XG9sK5oGkMMi8RHtMfD4pn3tOyq1bg6duIaT8+tdjlwguEW0rkhVxk/9q0z9Dixwi2+IUcwwvGXb77TAwQONtt6n1xyiu6xQUKIW7Z55RMci/+x2cu2hik6XR2NgY2wqWGcPbF9EfEyHURWVDg9oubDLBxo0qDfGXvhS7O4FEtK1IVJWuR4kV/9b/1j8AX5NS/nz+Ty0Pu3fvlvWx3FAeS/78Z+Ws0XPDDSrXS5IzIlbCIIa2r2ViaEjdcofdb/vMM7Bvn3KEFhYue9uiwUK2teiVp5Tyf4BrgD7gFHBVIkycYOArg64uFc3jx3H77erJww/DL36xPG347W9Vcq63vx0efTQqGQ6NrkAUimHtS0esxyQ7O8zE6XarbXh/93cq5LiqCh57bEn1JKJtRSxJp0sAB4CUsnuB05cFu3233LPHWFcG1qZ6xIc/NL9gQkoK3HMP7hs+gYyBDov1cAvin/8vaOmONXbuhDvvxPfev8ezanXU640lqakxuvKsq5P1e/YEHEbh9n1qv59Q743+fSHm/s3vsNFHGAX5hUIdT6H3w/7Xoc6iUEwmMJv8UTua00hzGPmjnCJ1GHlz8vCKhVfszMOnMX/og/DKK3P/eOmlcOedeM9/O17T3EtVIZR/R/h0jjV/W6XFisc719ETVLcZzB5dVJSUqkBNfNbvXZ3PnyeEv+8ul9KT7O9HVFfPa1uReNvfB3wXKERdeW4CDqFyrMeVycnJoFDZhTzioXaoF0SGQHhmOHtdyNuuP08rZyFv+/bdu0FLo+rzwaWX0lBbS90nP6niN/30O9W6ksulbovWrg32todzAHu9apvo4KC6uD10SB2Dg2q73fveBxdeuAOefnr+zvjgZA+zgrMpKSqyTfO0WyzBuxH0P94jRxooKambdZJq4a7ab1Qv3myxqM95PIEoPL0gdbjvVkM/WcRyK+akFj+4EOE84WeAmO+jEW4lCBLyJTic2GyGxsYG6urqwGRdMJY1qCZhUgMeoiYjPVD/pook/tWv5paxaxf88Y9rKX7hhQXbLD3Q36d+I1arsm/NroUgsIVEx/6GBmpr65icVL+nsTE19qmpSj80M9PfvdTURcdsoR0cCFNAPDk/f8FyIvG2n3UCuFgfhveG/uAHUoLytl94Ybxbs2QS0SO6lMPw9iVjNyZPPy1lfn7gUvq666ScmIhJVVLKxLStmCWAE0I8JIQ4JYRo0b2XI4R4XghxxP+o9+J/SQjRLoQ4LIS4NIJ24XK5IjktMTlxQkUZAa0f+hC8/LJa8zQwrQuGJkiBAAAgAElEQVRtADQghrYvP7Eak8sug5YW+MY3VMaCRx8N6NoEEaVbg0S0rUgmz2G/DN0rwKNCiO8BkXwjDwOXhbz3ReBFKWU58KL/NUKI7cC1qKWAy4AfCSHMLIJhxWp9PuVZHx0FoFQLj7zjjqiJx0ZElO95S0tLo1pevDGsfemI5Zjk5sIXvqCWg+bcCjscKn1reroKhB8bW1JdiWhbkUye7wcmgTuAZ1AZLv9+wU8AUspXgNCNi+8HHvE/fwT4gO79X0opp6WUnUA7sOiuRLc/OZShcLtVdJFuXch53nnqyfCw0vM8cSJ29Xs86lLhssvUYmZxMXz1q9DXt+SinU5nFBqYOBjSvkKIy5j85CdK4KahQdn7j34E1dUq+dFZ7upIRNta1GEkpZzwP/UJIZ4ETvvXAs6GfCllr7/cXr8HH2ADsFd3Xo//vTkIIW4BbgFYv76A9vYBPB43Pp+P1NR0xscHyckp5OTJToqLK+noaKSsrI4jRxyUl9tpb29gy5Zajh9vZcOGUk6fdpKZmcPU1CRgwmSyMjU1TmZmLqdP91BYWEZX1wFKS2tob3ewZYud9nYHZWV2urqa2bChgv7+brKy8pmYGMZiScPjgclJF15vFoODfUxOFmM2t1FaWo134GV2v/vdOC6/HPv0NE2nTrHh9ddp/eUvKfL5GDCbsTU1caLfRHuHj+npdDIyBsnKKgQ62bixEqdT9amjQ7Wno6OBkpJaenpayc0tZXDQSUqK6pOUJsxmK1brOJvXp9PbsI+yiQkO3HwzNddfjyM1Ffv0NI7XXsO+ezf7+kZITa3g9OluVq3Kx+cbJi0tDasVvF4XGRlZjI72kZtbTG9vGxs3VnP0qINNm+x0dZ3CYtnC8eNNbNhQxalT7eTkFDE+PkBamg2fzw34SE9X45SdXYjT2UlhYSWdnY1s2lTHsWMOSkrsdHU1sGlTLU5nK3l5pYyMOMnIyMHlmkQIE0JYmZgYx2TKPUtTDI/evgoKChg4eRK314tPStIzMhgcHKSwsJDOzk4qKytpbGykbtcuHA0N2HftomH/fmp37KC1rY3SjRtx9vWRs2YNk1NTmMxmrFYr4xMT5Obm0nPyJJs3l3Hw4AGqq2vYv9+B3W6nocGBva6O5pYWKsrL6e7uJj8/n+FhNRaglhSysrPp6+tj48Zi2tra2Latmv37HVRV2WlpcbBzp53+/tMUF2+i40gbRfn5DAwMYFu1CrfHg08I0m02BkdGKCjcMNunpsb91O3YgWP/fuylpTQcOULtli0cOHWagcEt7N3r5MSJHDIyJvF6TbjdVlavHmdkJJfrr25jV2YmBx59lJqZmYB9+R+bn3uO0vMu4FD7SVJT83G7h7HZ0khNBZ/PRU5OFv39fRRv3EhbWxvV1dXUOxooKdnMG2+o3/ChQ03k5lZx+nQ769YVMTMzQHa2DaQb38w06SkpDA4PU5iXR6fTSeXWrTQeOEBdXZ0aJ7udhoYGampqaW1tpaSklN5eJ9nZObimJjB5vVinpxk/dWphY5lvMRT4G+AvwO+BXUALcBLlcb9svs+FlFECtOheD4f8fcj/+EPgI7r3/5sIdEMNr7e4b5+UZrPS8/zyl+PdmiWTiJqLSzkMb18yNmMyMCDlDTcEnEUgZVGRlG++GfWqZklE21rotv2/gK8DjwH/D7hZSrkeeCdwd4T/yEPpE0IUAPgftam9B9ioO68ISLzr9Gji9cItt4DXi2lmBr71rfDy3AZiJehfrjRiMSZr16pAuMcfV3fnN94I+/fDOedEvapZEtG2FmqRRUr5nJTyN8BJKeVeACnlUtxeTwAf8z//GPC47v1rhRCpQohSoBwlvLwgItbqvrHkW99SFgdYtZQBN98cnAvEYFhjmRMjDhjavvzEckze9z547TX47/9WzqM5/OUv8OMfB3IZLYFEtK2FJk/9ym5o1qZF1zyFEI8Be4CtQogeIcRNwDeAi4UQR4CL/a+RUh4Afg0cRDmlPiOlXHQW8UUhpDAuPPusClvzM66lVd2zR3ncl4Njx5TTKgqOIo3x8fGolZUIGNa+dMRlTHp7lePzXe+CT38aystVUrglOOAS0bbmDc8UQniBCVTwwSqUxx3/6zQpZdz/FezcWSefe64hVOYwIkIjjKRcPMJFW+EJV5ZGaITRyEigzLw8v57n//cfs2mHmZxkPC0Nm/625NprOXrRjbS3B95atUrFEWdkqAid+SKjLBa1OyQzw4t5oE9lgRsYCDReyxDX2ak+YLXCFVcg3/d++q2FnDqlovGEUH/Sks2Fy4Gm1akV7XKNk5Zmm03Qpp0fGrm1UFTXfN9rOM3UyUm1QeHCC2MTnllXWysb9u0LCKYuZFxaJ/RhhfrOaQbn/2KkMM1JRKYlXxMmMfdLmU9HVKfrqUV1aeUKAdPT46xZY8OELxCeqX1eJwgbFJ6JVF+wy6Xs1J9IbWb1WhqbzQwPz/817CrsI+8b/wfCOVuKi+H66/Hs/hv6R1JwuQL52UK1YrU+aX0ZGxsnNdWG2x3QfTX5g4FmE8f5dKK8Uga+b2385hm2oK8TiZiZVobV24uoq5vftuZbDDXCUVVVtfQV4Rjj8Ug5PCzl0FAYrcORESk3bZKHPvhBJZIYwtCQkvyMZeRGNDl06FBc6iVGDiMj2NdiRHtMpqdVYsFQoXiTSel5xioQKBFtK/FWYc+A1MXijhMAs1nFpIfNI/TFL0JXF2WPP66cRyEqPllZKrw2AXNfhaWsrCzeTYgqRrCvxYj2mKSkqMSCLS1w553qzmT7drW8+aUvhddOiQaJaFuGnjwNHT73+usq9SBw4GMfU2uQd90V1yYtlQMHDsS7CVHF0PblJ1ZjkpGhljHHxtRE+o53xKSaWRLRtgw9ea5atSreTTg7BgfhIx+ZXQyt8U+ifPe7c6XiDERNTU28mxBVDGtfOmI9JlbrPMuJe/fCtdcqvdjvfIcg+bOzIBFtKyHSaZwthhSrPXYMrroq4LBBiSHb77tPTaZXXaX0vi4LlQU4Q6amVHjcnj3Q0aEyrxUXq8fhYeVImpqCd74Trrxycem1CHA4VITMSsGQ9hXCso/J6Cj80z/Bgw8GPDGvvQb33qt2klx7LWi7S86AZe+HlOGdXjoiFkNORGpqdstnn62fV1JxPgel9rdw3vZQ76f+8/N527Xztb9rnuBQPc/ycqgtGlDePI3bboM33oCnnlLueH9hR12F1Ncrh+eaNWpjst7brtWn18EUQnke12b7SBk+FfGeUe8qG+39a+jrU2WmpQU8oFqmzsXWsrS+653KZnOgjZqTN5LdC6Hfq/ao9dnrT3A4Oqrs+/LLYySGbLfL+r17A4Yyn7io/rne067fvqHLVukT5iDxYr0tzgohL0aoxx0RVK1WrjYOJuZp0+wJugHWtjPMzKgv2uMBiwW3LZvGZjP/+78qEUHozpS//Vv4/tfHyTaNzN9uIfBm5zI0kcLMjKpak+C0WoN/j/psoKE7CqQMeOqtZl/wLoIQwV69qPjsEOrO1dvhrLd9ZAQ8HkRR0cr0tm/bti16brV48OijUoKsv/12la/d4NTX18elXmLkbTe8fcnYjUljo5QXXaSmndWrpbzvvth52qVMTNsy9JWnoRN0DQ1BZWXwrcETT8xNCpdkUZIJ4OLH0aPqhmm1sTK3RMySEsAlMlNLXISOK5/+9OzE2awlg/vkJ+H06Tg2amk0NzfHuwlRxdD25SfWY7J58wIT59iYijaKAoloW4aePNNCcqwYhrvugl/+cvZlxW/9WZ2dTuW8iVYo2siIqudf/kUlnnnzzeiVHYaKioqYlR0PDGtfOuIyJk1NKvA9O1ulHK6uVgHwS/hnlIi2ZWhv+9TUDL29wQnJzgR9ArVwC/j6v2kL2KFODr1DQztHH545OBjsMNrlrUdYrUqAGODgQbrz8ynPyQkU9tRTdNj/kfp6tV6fmamEF0IdRl5vIIGalhAxMxPybFOYGhvUXtKZGfWhsTF4/nl1WCzqwxkZ8O53c9C0g5YWVVZmpqpnzRq1OV8LmdOHR87Xd6ezm4KCcn0UIlIGwun06/nzlRVKqMMIVFlTU+oivafnzMb8TJiZnlbfU6QGpvdq6DNVap4bqxVpTZlNeKf1S5+c02Ri4dBMdH/XOz38DiPNaaR9rLOzm4qKckxCF3+sdxgFVexHK2RmRh1eL5jNuNPXcLDVxHPPETZEMzsbPv33x0l/+mk47zx1aJw8CT/9KVx0Eb6t2xidsuJ2BycHDA3N1IeaHjvWTUlJeZBj0mr1h2V63IHvWfvh+vukd2rNyU6qfXcyxNekhWcuJsA832KoEY7a2trorQzHg8OHpVy1So5s2iTlvffGuzVLZmRkJC71EiOHkeHtS8ZmTCYnpfzRj5SGJ0gphNL37OuLelWzJKJtGfq23Wtg+TZArXFOTTG8eTP867/C8ePxbtGSGF5IMcKAGN6+iM2YrFoFn/qU2j78yitqy/LPfgbr1i3+2bMlEW3L0JNnIgqkRsyDD8JLLwGQNjQEExPKIg3MSlgj1GNo+/ITyzFJSVFhmZs2xayKWRLRtoxvHUZk/364/fa57z/5pBJJjiYHD6qUxi+/HPVsmUmSzMvBg/DDH8IPfgCNjfFuTUwwtMPI7fbhdAYHf8y3ph9uO6veB6D/u9kciKrRL1zrpQLDlaN3GE1PMweTCTYP1iNu+wwUFKg3vV5cW7YEvOAPPgi5uTgvu5GJibllWK2BeGKvV9Wlv7tcswbyrMOIv76iPO319YEGZ2ercMytW1XIUkEBo+V2+obCp9jVFvL1UR+h/da+eynh1CnX7Dmak8nnUz4HvQ5uON/LQr4YfT0w//cbbXyac2UBZ9Ecp5dmKJpXSPu8xQKpqXiwzDrQ9BFAeqdRUMELGXSIw0izP72e58SEa9aZGPRHrWy9Z1RfttYHzcDMZtypNk70mhbUNC6mm9R77lYaDVo93/8+VFWp0OMLLsBXXMKEy4zHE/itWSwBm9H3Q/vNjY66sNmCg4csFjD5PMEeMp0DLLSr2uOiDiP3jPLyLrLNytCb5Hft2iX3+1NZJBpe71xv++bNav6aZXoazjmH0clJMn/5S9gd2Is7NQWHDs0NzwwXgq4ZWWiUXSRICc3N0NWlwjHz89WRna0mzjNhdHSUzMzMM/vQWeJ2q22yHR2xE0NOZPuKlGiPicsFr74KjzwCv/61+scIynT/8z/h0kvPfNdLJCynbelZaJO8oa88PQl8G2o2qwkoaLIM5atfheZm+q68kswbb1RCHhY1JKtWQV1dZHWF7jQ5E4SAnTvVsVT6+vqWzcCtVtiwQR2xIpHtK1KiPSZpafCe96jj3nuVjufGjXDuubGZNDWW07YiJS5rnkKIY0KIZiFEoxCi3v9ejhDieSHEEf/jQtMOAClnemmUSBw8OLu+Wfzii+ry75574tyopVFcXBzvJkQVQ9uXn1iOSW4uXHON2s4Z61x5iWhb8XQYvUtKWau7JP4i8KKUshx40f96QQwrVjs1pfQ8/fc8bddco97/t3+L/uK61wttbUvWU4yEtra2mNexnBjWvnSslDFJxH4kkrf9/cAj/uePAB9Y7AOGFKsdHVUL57q1tOqHHlJPpqfhve+Ft95aej2dnfCVr6j72q1b1WXChz8Mf/xjVDNm6qmuro5JufHCkPYVwrKPic8HTz8NN9ygFkIvuAA+/3n13hJCgxPRtuLiMBJCdAJDqBTG90spHxBCDEsps3TnDEkp59y6CyFuAW4BWLcu3/7KKy34fG58Ph9paemMjw+Sk1PIyZOdFBdX0tHRSFlZHUeOOCgrs9Pe3sDmzbUcP97KunWlDA87sdlymJ6exGQyYbVamZkZJy8vl76+HoqLyzh8+ABbttRw6JCDLVvstLerx66uZjZsqGBgoJusrHwmJoaxWNLweJSX0+3OYni4j/HxYiyWNkpKqmHsL+xOseIYG8O+ejVNb7zBzNgYq3fvpmj1agbcbmxpafRkbObIER/T0+nYbINkZRUihOrTiRONlJbW0dHhYMMGO93dDeTm1jI83EppaSnTp4+QMzzEpNeLCbCaTIx7veRarfRMT1O2ahUHZmbYad/N/6s/wuHDdtLTHYyO2tm4sZnc3ApMpm5yc/OZnh4mJUXtsZuZcZGensXISB95ecX09raxeXM1R444qKiwU1//PLt2XczRo02UllbhdLaTk1PE0NAAFosNr9eNlD5SU9OZnFR9GhjopLCwkq6uRkpK6jh2zMGmTXa6uhooLq6lt1eN08iIk9Wrc3C5JpHSxPS0lf7+cY4ezeW221ZHzWGkt6/8/Hx7S0sLbreyr/T0dAYHByksKOBo5zEqKytpampk1646Ghoc1NXZ2d/goHbHDlpbWyldvx5nfz85a9Yw6fFgSklFmFMYHh4nMzOX/v4eSkrK6Og4wM6dNTQ2KsHfhgb12NzcTEVFBd3d3eSvW8fw8DBp/gndNTVFVlYWfX19FG/aRFtbG1VV1TgcDmpr7TQ2qse//OUFLrzwIo4ebadowwYGBgawpafj9nhUnzIyGBwaCvRp61aa3mqirqYGR0MD9qoqGpqbqa2q4mBnF9aUzbS0OPnrX3MYG5vE4zHhclkpKhrnggty2bbhIOXuGQ5MTFBjs83auWNsDHtWFs1AWXUNR4/1kJOTz/j4MOnpaZhMMD3tmu3Thg3FtLe3sXVrNW+95ZhdFqipsXPwYBPbt1fRebSdosICBvr7sa1eHeiTzcbg4CDr1xfS2dnJ1q2VvPWWGqf9DQ7su+00NDRQW1ND6+HDlJSU4nQ6ycrKYWpqEouQWGdmGO/ro/Rtb5vXtuI1eRZKKZ1CiHXA88BngScimTz1GF4y7NlnA4rxn/gEPPBAfNtjUJKSdMuPlGqFae9etTvjssuMk6jwTEg4SToppdP/eAr4A3Au0CeEKADwPy6sgY/B0yRMTSlZOlQaDh58UKUrMDAOhyPeTYgqhrYvP7EaEyFg1y4VFHfVVbGfOBPRtpZ98hRCZAghVmvPgUuAFuAJ4GP+0z4GPL5YWelG/lf3L/+ilGRB5S+SEm66aU76YSOxkvIXgcHty89KGZNE7Ec8rjzzgVeFEE3Am8CTUspngG8AFwshjgAX+18viGHFah9+GO67b/Zl0623qieHD8P110cvjHL/frVwv3mz0lf8058izmt0NjQ1NcWs7HhgWPvSEbcxGRhQ9tfRERWbS0TbMnSE0c6ddvnss445ubkWkj8MPbTztcguqzWQkMpkCsgaatF2+iRx+rq0Ry08c3JSRRadPh2s52mX9YifPhDQVDx1Cs/ICJby8kChb387nRfewL59AT3P0ARwMFdy0WJRTvW1M72Ip55UieVCxzc/H979bti5k9Mil4YGFV2k9SEvD9avh6wsdStmsQQnt9OjhRUGwiY9CBEcdxEuhDQ0InCxRH2hyfrc7kACUKcTPvGJ2Kx52u126dDWPMNsZAzXVny+gMiqFp5pNkNKCjIllRm3YGYmOGpTC00MmwBuoaRzIQYcatcAXq8Hq9USKDec4YfWo4+P1IRYTSa8llRO9QtaW9WcGMr27XD+hm7E439U+5a1ejIywG5XkRglJfhWr2HGLbRig/uua6ZeetTjUf3QJ8oTUhe/GaLlOV84driB00JbZ78Sr0fp3w4MICoqEsthFC127NghW1pa4t2Ms6e3F3bsoPXii6m84gq19zMOjIyoyT4nZ2kZiFtbW6msrIxewyLg1CnIz4/N5Gl4+yJ2Y3LypIouOnBA/WO/+GI1ecZqs3w8bAtWcHim4SNAPv95GByk6JVXlJDC5ZerGWyZWbNGHUul6CzycS+VWGpIGt6+iN2YrF+vUrAvF/GwrcVIpE3yZ4yhY4//9Cf4zW8AGKiuVmtEd94Z50YtjYGBgXg3IaoY2r78rJQxScR+GHryNKxYbVeX8qz7sWmJeB55BH7xi+jVI6XSCP3nf4ZHH425jpvNZotp+cuNYe1Lx0oZk0Tsh6Fv26WUMdP3NZskQvrwCXNQEqkll9voQFx3rdKqM5sBcK9ePfucW26BqSm8H78ZydktIFlOOeHxx+FHP1K6dhpf+ALceCN84AN4q3YizdEd/ulp94rSW5aa5yICLU+NWSeG5iWTMij523yO5yDN0sW0PMNVHqJLqTEz4w50IZyjKVw9oR4bfwE+S8qivwXzqV7Ea6+qDJpDQ2o9qLRUhQiXlEBeHr7UVXMSs2n916rXa3FCwLaCvid9RkZ9kj4hIncYhfneTNKrnACLpP4w9OTp9ar+zWcDoQ5EvT2EZsvUsvilpirZLWESgBkTwR53TRhZX77eaamJ/2re9sHBwBiUlUGt3Q5HjgROvvhifGazmuR0HnfncaVj7HLN9bbrxZA1b7vHo97Py4M1hYVq9/ICaT18bjjaoao9dkyVlZenQuHz8pTNp6YGC9SGOoFnBWn9Rj8x4WN0NPg3qjlsI/G2h/629fXotUpnZtTvsrMT9u2bt4vRYYEJLOyfhCmwbSP0T/7P6LOdamLAQWVF4nUJc044zWYpfYH3FlMMDy0o5MpbSGXXx48rUbCmJrXjITNTySdeeCFs3lwA//AP6pgHfanzCWNr1xIBm/Bpao26gua/M1iwiyF/nNsGs0pGP29Cev/njOxtN7xY7fe/D5//PEPl5WSvX69SZcRa2yuGDA0Nkb2ggGlsiFV4puHti/iNSbRJRNsy9KKOoRf0e3pUxkxgsKIC/vpXlYLQwAwODsa7CVHF0PblZ6WMSSL2w9CTpzXMrZEh8Png4x9XG3GBwjfeUO/feae6h44mMzPq/ipcQqQoU1hYGPM6lhPD2peOlTImidgPQ0+eM1oCFSPh8cCtt8ILL8y+1akpK42MqL2eJ04svZ6BAZXmo6gIamvVxrw771QLhTGiM4ZlxwND2lcIcRsTKVkwU9wZkoi2ZWiHkdmcxsmTC6+Dhwv58nqDnRja+n5KisodlJ4OZp8bpERarMy4BS5XIOIu1KGxkMPo9OmAw6i8HHZbDyAuvFCtrgO88AKVhw6pW3ZtRbypiWMzheyrF0EOo6wssNkCvojQ8EzNYbR2qgdef115N0NTe+zZA243AzkVtLaqubS/X1WthWauXavWykMzZ4b2W5/5ESAvr5JTpwIONS0JY6iTLfQ70zNfhk59PdPTyhHX3q4k0WLFQrnC53MVCHSdDjEwn1AZI7XVgFBn2HzfyaKN8H8onMNt69bKuRsGIvXm6w8h8Akzk5PKpk+cUM6i06fVb6a0FLZtg/Upg0qjobtb/WC0rIIFBbOG5TNZFrSHcCHUZWWVeL0h54aGmob5LvTM6a6/X3O87fh/xIuINxvaYbR9+3Z58ODBeDfj7Kmvh7/5Gxo+8xnqMjLg61+Pd4uWRENDA3WRZq2LEio+OjYOI8PbF/EZk1gQr34s5DAy9ORpaLFarxfOOSeQjsNqVc+rquLbLgOSFENOEitWrLfd0GK13/727MTpuP12dYtz003Rk6OLA4koWLsUDG1fflbKmCRiPww9eRpWrPb551VyNj92TdvzjTfgjjuiV09nJ3zuc8phdOONSiYshiSiYO1SMKx96VgpY5KI/TC0w2h8fDLIYQTzRxnpo4s0h5Hm/PHLLZKWppxFq9IkJs9MwGHkNTM9HYjkidRhNDwcvOaclQV20YD43vcCDqPJSRouuYS6V19Vrw8dgp/+lGPvuZn2jkBn0tLO0GH0xz/Cn/+sGpybqxbw77wT3vMe5LXXcVwUh9WpTUsLjmQymYLX5fWRcHqHkZRw+HADmzfXBS32+3zhNVDDR8MEP2rnhjqmpqeVo6K/n5iy0JWn9IUP25x1ZOhD2PwilD4p5kS2BR1If2RbaGWL6Hsu0Kb9+xvYVbsrfLmLETLoEqVFOjWldtmNjio7t1qVbebmgs08pXZ6jI4q20tNVUZrsynPUkoKPjk3wgfp77uUs2HJ+m7v39/Arl11wefqTwoNJwxnXPOF2eq/N61sj2c2Nfh8JNc8Y8z0dMDbnpcXElE2Ogo7d+Lr6cH05JNw6aVBnx0aCoRnZmREv20zMyoL8diYsvG8POVlP9sgJ5/Pt6xiGkNDalPBFVck1zznY7nHJFbEqx8rds3T5XLFuwmLkpqqdmrk54cJxf3iF6Gri9brrlN7P0OudLKz1Q6PWEycoK62N25UIrZbtqhJeinRoa2trdFrXARkZ6s097HCCPa1GMs9JrEiEfuRcJOnEOIyIcRhIUS7EOKLC51raLHaPXvg/vsBKH3mGSVT9+//HudGLY3S0tJ4NyGqGNq+/KyUMUnEfiTU5CmEMAM/BP4O2A5cJ4TYPt/57ihGMCwrQ0Pw0Y/OLn45zztPvX/PPcqZFE0aGlSyuVdeiW65YXA6nTGvYzkxrH3pWCljkoj9SKjJE5W/vV1KeVRKOQP8Enj/fCdb5mhUGYDubrjkkqAMWjltbeqJzwdXXw1PPbW0Olwu5TD6279VibfuuEM5qM45Bx57TK34x4CcOKQQiSWGtK8QVsqYJGI/Es06NgDHda97gPP0JwghbgFuAcjPL6CtbQCv142UPlJT0xkfH2Tt2kJOnuykuLiSjo5GysvraGtzsGWLnfb2BjZurKWnp5WsrFJGR51kZOTg802SlmbCZrMyMz1GXuZqenp7KSsvp/nQYbZsqaGlRZVx5Ih67O5uZsOGCvr7u8nOzmdiYhizOQ2vFyYmXMzMZDE01MfYWDEpKW1s2lSNmDnO7t//HkdnJ/bSUpq6u8lOT6fP7aYoO5uB8XFsaWk4m3o50u5jejqdjIxBsrIKAdWnEyca2by5jo4OB4WFdrq7G8jNrWV4uJWyslIm+4+QU1nJ5L33YjKZsJrNjE9Pk2uz0TM0xJbRMfbWtzE+XsPYmINTp+zk5TlIS7OzalUzBQUVTE11k5OTz+TkMCkpaUgJMzafDGkAAAZ+SURBVDMu0tOzGBnpY926Ynp72ygrq6atzUFlpZ39+1vYseMddHQ0UVxchdPZTnZ2EaOjA6Sk2ILGaXJS9WlgoJPCwkqOHWukpKSOzk4HJSV2uroaKC6upbe3lfz8UkZGnKxZk4PLNYnPZ2J83EpPzzj19blRNUC9fRUUFDAwMIDb7cbn85Gens7g4CCFhYV0dnZSWVlJY2MjdXV1OBwO7HY7DQ0N1NbW0traSmlpKU6nk5ycHCYnJ9VYWK2Mj4+Tm5tLT08PZWVlHDhwgJqamtkytMfm5mYqKiro7u4mPz+f4eHh2ZBRl8tFVlYWfX19FBcX09bWRnV19ZwyDh06xLnnnkt7eztFRUUMDAxgs9kM16f169dz9OhR7HY7TU1NVFVVLUufFrSVRPK2CyH+AbhUSnmz//VHgXOllJ8Nd35NTY1MxHzOZ0pvby8FBQXxbsaSiVc/YhVhtBLsK2lbS8NI3vYeYKPudREw72KHMLBwsJ6VIH0GK6cfGivBvlbKmCRiPxLtytMCtAHvBk4A+4APSSkPzHP+GHB4+VoYM3KBxEsPeObEqx+bpJR50S50hdhX0raWxry2lVBrnlJKjxDiNuBZwAw8NN/E6edwLG7XlhshRH2yHwmJ4e1rpYxJIvYjoSZPACnlU8AS3c1JkiRJElsSbc0zSZIkSQyB0SfPB+LdgCiR7EdishL6sxL6AAnYj4RyGCVJkiSJUTD6lWeSJEmSxAXDT55CiLuEECeEEI3+4/J4t+lMOBMhlERGCHFMCNHsHwNj67j5SdpWYpCotmX423YhxF3AuJTyO/Fuy5niF0JpAy5GBQjsA66TUhou65gQ4hiwW0q5EvYUAknbShQS1bYMf+VpcM5ICCVJkjMgaVsxZqVMnrcJId4SQjwkhMiOd2POgHBCKBvi1JalIoHnhBAOv7jGSiFpW/EnIW3LEJOnEOIFIURLmOP9wI+BLUAt0At8N66NPTPCBU8bdR3lAillHUqL9TNCiHfGu0GRkLQtQ5CQtpVwEUbhkFK+J5LzhBA/Bf4c4+ZEkzMSQklkpJRO/+MpIcQfULeNsVdgXiJJ20p8EtW2DHHluRBCCL1O1ZVAS7zachbsA8qFEKVCiBTgWuCJOLfpjBFCZAghVmvPgUsw1jiEJWlb8SeRbcsQV56L8C0hRC3qluQYcGt8mxM5ZyGEkqjkA3/wS7hZgP+VUj4T3yZFhaRtxZ+EtS3Db1VKkiRJknhg+Nv2JEmSJIkHyckzSZIkSc6C5OSZJEmSJGdBcvJMkiRJkrMgOXkmSZIkyVmQnDyTJEmS5CxYCfs8ExIhxFrgRf/L9YAX6Pe/npRSvi0Gde4CPqPlvY9CebcBE1LKn0WjvCTRIWlbiUFyn+cysFzSZkKI3wBflVI2Ram8dOA1KeWuaJSXJPokbSt+JG/b44AQYtz/eJEQ4mUhxK+FEG1CiG8IIT4shHjTL/66xX9enhDid0KIff7jgjBlrgZ2asYthLhQJ+K7Xxfi9n/9ZbwlhPgP3eev97/XJIT4OYCUchI4JoQ4N/bfSpJokLSt5SN52x5/aoBtwCBwFHhQSnmuEOLzwGeB24HvAfdKKV8VQhSjQu62hZSzm+CY339C3Wa9JoSwAS4hxCVAOUpYQQBP+BVqTgNfRqnXDAghcnTl1APvAN6Maq+TLAdJ24ohyckz/uyTUvYCCCE6gOf87zcD7/I/fw+w3R/fC5AphFgtpRzTlVNAYN0L4DXgHiHEo8DvpZQ9fgO/BNjvP8eGMvga4LeaUreUclBXzimgcundTBIHkrYVQ5KTZ/yZ1j336V77CIyPCThfSjm1QDlTQJr2Qkr5DSHEk8DlwF4hxHtQVwR3Synv139QCPE55td6TPOXncR4JG0rhiTXPI3Bc8Bt2gu/0k8oh4Ay3TlbpJTNUspvom6PKlG3ZDf6b7UQQmwQQqxDeW7/0e/FJeTWqoIEkQBLEhOStnWWJCdPY/A5YLd/0f0g8MnQE6SUrcAabfEeuF0oRfQm1H/3p6WUzwH/C+wRQjQDvwVW+6XKvga87D//Hl3RFwAvxKxnSeJN0rbOkuRWpRWEEOIOYExK+WCUytsF3Cml/Gg0yktiXJK2NZfklefK4scEr3MtlVzgK1EsL4lxSdpWCMkrzyRJkiQ5C5JXnkmSJElyFiQnzyRJkiQ5C5KTZ5IkSZKcBcnJM0mSJEnOguTkmSRJkiRnwf8PAgjriR2K2I4AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Stack over all traces\n", "trR_stack, trT_stack = ut.stack_all(trR, trT, pws=True)\n", "\n", "# Plot as wiggles\n", "wg.rf_wiggles_baz(trR, trT, trR_stack, trT_stack, 'test', btyp='baz',\n", " scale=1.e3, tmin=-5., tmax=8., save=False, ftitle='porter2011',\n", " wvtype='P')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And VoilĂ ! \n", "\n", "Now try the same example but setting `wvtype = 'SV'` to get S receiver functions for the same model. Such an example would correspond to a core-refracted shear wave (such as SKS) with no incident transverse component. In this case the receiver functions are unstable (spectral division by zero-valued SH component) and are not computed for the transverse component. Setting `wvtype = 'Si'` will produce a transverse component receiver function, which is more realistic for incident S waves propagating through the mantle only. Be careful with the slowness values!!!\n", "\n", "Did you notice the boolean `pvh=False` in the code above? It sets whether or not the seismograms are rotated to the P-SV-SH wave modes. Setting it to `True` will essentially make the zero-lag signal on the radial. component disappear, since we know the exact value of the seismic velocities at the surface. This may not always be true so use with caution when comparing with real data!\n", "\n", "## References\n", "\n", "* Frederiksen, A.W., & Bostock, M.G. (2000). Modelling teleseismic waves in dipping anisotropic structures. Geophysical Journal International, 141, 401-412. https://doi.org/10.1046/j.1365-246x.2000.00090.x\n", "* Porter, R., Zandt, G., & McQuarrie, N. (2011). Pervasive lower-crustal seismic anisotropy in Southern California: Evidence for underplated schists and active tectonics. Lithosphere, 3(3), 201-220. https://doi.org/10.1130/L126.1" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.3" } }, "nbformat": 4, "nbformat_minor": 2 }