Spheroid/Spectral Components Analysis [WGM2012 Free Air].ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " * postgresql://localhost:5432/grav_29_1\n",
      "Done.\n",
      " * postgresql://localhost:5432/grav_29_1\n"
     ]
    }
   ],
   "source": [
    "%load_ext sql\n",
    "%sql postgresql://localhost:5432/grav_29_1\n",
    "%sql create extension if not exists postgis\n",
    "%sql"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plotting modules\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def raster_filter_range(raster0, g1, g2):\n",
    "    raster = raster0.copy()\n",
    "    raster.values = raster.values.astype(np.float32)\n",
    "    raster.values = gaussian_filter(raster.values,g1,mode='constant', cval=0) \\\n",
    "            - gaussian_filter(raster.values,g2,mode='constant', cval=0)\n",
    "    return raster\n",
    "\n",
    "def spectrum(raster, gammas, dgamma = 1.0):\n",
    "    rasters = []\n",
    "    for g in gammas:\n",
    "        print (g,\". \", end = '')\n",
    "        _raster = raster_filter_range(raster, g-dgamma/2, g+dgamma/2)\n",
    "        rasters.append(_raster)\n",
    "    return rasters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  25.,   75.,  125.,  175.,  225.,  275.,  325.,  375.,  425.,\n",
       "        475.,  525.,  575.,  625.,  675.,  725.,  775.,  825.,  875.,\n",
       "        925.,  975., 1025., 1075., 1125., 1175., 1225., 1275., 1325.,\n",
       "       1375., 1425.])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Gaussian filter sigma, km\n",
    "sigmaskm = np.linspace(25,1425,29)\n",
    "sigmaskm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Raster from file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><svg style=\"position: absolute; width: 0; height: 0; overflow: hidden\">\n",
       "<defs>\n",
       "<symbol id=\"icon-database\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M16 0c-8.837 0-16 2.239-16 5v4c0 2.761 7.163 5 16 5s16-2.239 16-5v-4c0-2.761-7.163-5-16-5z\"></path>\n",
       "<path d=\"M16 17c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "<path d=\"M16 26c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "</symbol>\n",
       "<symbol id=\"icon-file-text2\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M28.681 7.159c-0.694-0.947-1.662-2.053-2.724-3.116s-2.169-2.030-3.116-2.724c-1.612-1.182-2.393-1.319-2.841-1.319h-15.5c-1.378 0-2.5 1.121-2.5 2.5v27c0 1.378 1.122 2.5 2.5 2.5h23c1.378 0 2.5-1.122 2.5-2.5v-19.5c0-0.448-0.137-1.23-1.319-2.841zM24.543 5.457c0.959 0.959 1.712 1.825 2.268 2.543h-4.811v-4.811c0.718 0.556 1.584 1.309 2.543 2.268zM28 29.5c0 0.271-0.229 0.5-0.5 0.5h-23c-0.271 0-0.5-0.229-0.5-0.5v-27c0-0.271 0.229-0.5 0.5-0.5 0 0 15.499-0 15.5 0v7c0 0.552 0.448 1 1 1h7v19.5z\"></path>\n",
       "<path d=\"M23 26h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 22h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 18h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "</symbol>\n",
       "</defs>\n",
       "</svg>\n",
       "<style>/* CSS stylesheet for displaying xarray objects in jupyterlab.\n",
       " *\n",
       " */\n",
       "\n",
       ":root {\n",
       "  --xr-font-color0: var(--jp-content-font-color0, rgba(0, 0, 0, 1));\n",
       "  --xr-font-color2: var(--jp-content-font-color2, rgba(0, 0, 0, 0.54));\n",
       "  --xr-font-color3: var(--jp-content-font-color3, rgba(0, 0, 0, 0.38));\n",
       "  --xr-border-color: var(--jp-border-color2, #e0e0e0);\n",
       "  --xr-disabled-color: var(--jp-layout-color3, #bdbdbd);\n",
       "  --xr-background-color: var(--jp-layout-color0, white);\n",
       "  --xr-background-color-row-even: var(--jp-layout-color1, white);\n",
       "  --xr-background-color-row-odd: var(--jp-layout-color2, #eeeeee);\n",
       "}\n",
       "\n",
       "html[theme=dark],\n",
       "body.vscode-dark {\n",
       "  --xr-font-color0: rgba(255, 255, 255, 1);\n",
       "  --xr-font-color2: rgba(255, 255, 255, 0.54);\n",
       "  --xr-font-color3: rgba(255, 255, 255, 0.38);\n",
       "  --xr-border-color: #1F1F1F;\n",
       "  --xr-disabled-color: #515151;\n",
       "  --xr-background-color: #111111;\n",
       "  --xr-background-color-row-even: #111111;\n",
       "  --xr-background-color-row-odd: #313131;\n",
       "}\n",
       "\n",
       ".xr-wrap {\n",
       "  display: block;\n",
       "  min-width: 300px;\n",
       "  max-width: 700px;\n",
       "}\n",
       "\n",
       ".xr-text-repr-fallback {\n",
       "  /* fallback to plain text repr when CSS is not injected (untrusted notebook) */\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-header {\n",
       "  padding-top: 6px;\n",
       "  padding-bottom: 6px;\n",
       "  margin-bottom: 4px;\n",
       "  border-bottom: solid 1px var(--xr-border-color);\n",
       "}\n",
       "\n",
       ".xr-header > div,\n",
       ".xr-header > ul {\n",
       "  display: inline;\n",
       "  margin-top: 0;\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-obj-type,\n",
       ".xr-array-name {\n",
       "  margin-left: 2px;\n",
       "  margin-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-obj-type {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-sections {\n",
       "  padding-left: 0 !important;\n",
       "  display: grid;\n",
       "  grid-template-columns: 150px auto auto 1fr 20px 20px;\n",
       "}\n",
       "\n",
       ".xr-section-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-section-item input {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-item input + label {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label {\n",
       "  cursor: pointer;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label:hover {\n",
       "  color: var(--xr-font-color0);\n",
       "}\n",
       "\n",
       ".xr-section-summary {\n",
       "  grid-column: 1;\n",
       "  color: var(--xr-font-color2);\n",
       "  font-weight: 500;\n",
       "}\n",
       "\n",
       ".xr-section-summary > span {\n",
       "  display: inline-block;\n",
       "  padding-left: 0.5em;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in + label:before {\n",
       "  display: inline-block;\n",
       "  content: '►';\n",
       "  font-size: 11px;\n",
       "  width: 15px;\n",
       "  text-align: center;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label:before {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label:before {\n",
       "  content: '▼';\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label > span {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-summary,\n",
       ".xr-section-inline-details {\n",
       "  padding-top: 4px;\n",
       "  padding-bottom: 4px;\n",
       "}\n",
       "\n",
       ".xr-section-inline-details {\n",
       "  grid-column: 2 / -1;\n",
       "}\n",
       "\n",
       ".xr-section-details {\n",
       "  display: none;\n",
       "  grid-column: 1 / -1;\n",
       "  margin-bottom: 5px;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked ~ .xr-section-details {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-array-wrap {\n",
       "  grid-column: 1 / -1;\n",
       "  display: grid;\n",
       "  grid-template-columns: 20px auto;\n",
       "}\n",
       "\n",
       ".xr-array-wrap > label {\n",
       "  grid-column: 1;\n",
       "  vertical-align: top;\n",
       "}\n",
       "\n",
       ".xr-preview {\n",
       "  color: var(--xr-font-color3);\n",
       "}\n",
       "\n",
       ".xr-array-preview,\n",
       ".xr-array-data {\n",
       "  padding: 0 5px !important;\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-array-data,\n",
       ".xr-array-in:checked ~ .xr-array-preview {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-array-in:checked ~ .xr-array-data,\n",
       ".xr-array-preview {\n",
       "  display: inline-block;\n",
       "}\n",
       "\n",
       ".xr-dim-list {\n",
       "  display: inline-block !important;\n",
       "  list-style: none;\n",
       "  padding: 0 !important;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list li {\n",
       "  display: inline-block;\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list:before {\n",
       "  content: '(';\n",
       "}\n",
       "\n",
       ".xr-dim-list:after {\n",
       "  content: ')';\n",
       "}\n",
       "\n",
       ".xr-dim-list li:not(:last-child):after {\n",
       "  content: ',';\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-has-index {\n",
       "  font-weight: bold;\n",
       "}\n",
       "\n",
       ".xr-var-list,\n",
       ".xr-var-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-var-item > div,\n",
       ".xr-var-item label,\n",
       ".xr-var-item > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-even);\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-var-item > .xr-var-name:hover span {\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-var-list > li:nth-child(odd) > div,\n",
       ".xr-var-list > li:nth-child(odd) > label,\n",
       ".xr-var-list > li:nth-child(odd) > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-odd);\n",
       "}\n",
       "\n",
       ".xr-var-name {\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-var-dims {\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-var-dtype {\n",
       "  grid-column: 3;\n",
       "  text-align: right;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-var-preview {\n",
       "  grid-column: 4;\n",
       "}\n",
       "\n",
       ".xr-var-name,\n",
       ".xr-var-dims,\n",
       ".xr-var-dtype,\n",
       ".xr-preview,\n",
       ".xr-attrs dt {\n",
       "  white-space: nowrap;\n",
       "  overflow: hidden;\n",
       "  text-overflow: ellipsis;\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-var-name:hover,\n",
       ".xr-var-dims:hover,\n",
       ".xr-var-dtype:hover,\n",
       ".xr-attrs dt:hover {\n",
       "  overflow: visible;\n",
       "  width: auto;\n",
       "  z-index: 1;\n",
       "}\n",
       "\n",
       ".xr-var-attrs,\n",
       ".xr-var-data {\n",
       "  display: none;\n",
       "  background-color: var(--xr-background-color) !important;\n",
       "  padding-bottom: 5px !important;\n",
       "}\n",
       "\n",
       ".xr-var-attrs-in:checked ~ .xr-var-attrs,\n",
       ".xr-var-data-in:checked ~ .xr-var-data {\n",
       "  display: block;\n",
       "}\n",
       "\n",
       ".xr-var-data > table {\n",
       "  float: right;\n",
       "}\n",
       "\n",
       ".xr-var-name span,\n",
       ".xr-var-data,\n",
       ".xr-attrs {\n",
       "  padding-left: 25px !important;\n",
       "}\n",
       "\n",
       ".xr-attrs,\n",
       ".xr-var-attrs,\n",
       ".xr-var-data {\n",
       "  grid-column: 1 / -1;\n",
       "}\n",
       "\n",
       "dl.xr-attrs {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  display: grid;\n",
       "  grid-template-columns: 125px auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt,\n",
       ".xr-attrs dd {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  float: left;\n",
       "  padding-right: 10px;\n",
       "  width: auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt {\n",
       "  font-weight: normal;\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-attrs dt:hover span {\n",
       "  display: inline-block;\n",
       "  background: var(--xr-background-color);\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-attrs dd {\n",
       "  grid-column: 2;\n",
       "  white-space: pre-wrap;\n",
       "  word-break: break-all;\n",
       "}\n",
       "\n",
       ".xr-icon-database,\n",
       ".xr-icon-file-text2 {\n",
       "  display: inline-block;\n",
       "  vertical-align: middle;\n",
       "  width: 1em;\n",
       "  height: 1.5em !important;\n",
       "  stroke-width: 0;\n",
       "  stroke: currentColor;\n",
       "  fill: currentColor;\n",
       "}\n",
       "</style><pre class='xr-text-repr-fallback'>&lt;xarray.DataArray &#x27;z&#x27; (y: 5401, x: 10801)&gt;\n",
       "[58336201 values with dtype=float32]\n",
       "Coordinates:\n",
       "  * x        (x) float32 -180.0 -180.0 -179.9 -179.9 ... 179.9 179.9 180.0 180.0\n",
       "  * y        (y) float32 -90.0 -89.97 -89.93 -89.9 ... 89.9 89.93 89.97 90.0\n",
       "Attributes:\n",
       "    long_name:     z\n",
       "    actual_range:  [-388.63299561  990.96801758]</pre><div class='xr-wrap' hidden><div class='xr-header'><div class='xr-obj-type'>xarray.DataArray</div><div class='xr-array-name'>'z'</div><ul class='xr-dim-list'><li><span class='xr-has-index'>y</span>: 5401</li><li><span class='xr-has-index'>x</span>: 10801</li></ul></div><ul class='xr-sections'><li class='xr-section-item'><div class='xr-array-wrap'><input id='section-2fcdfd96-63b6-4c1e-9f82-47990658fefa' class='xr-array-in' type='checkbox' checked><label for='section-2fcdfd96-63b6-4c1e-9f82-47990658fefa' title='Show/hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-array-preview xr-preview'><span>...</span></div><div class='xr-array-data'><pre>[58336201 values with dtype=float32]</pre></div></div></li><li class='xr-section-item'><input id='section-05023f5b-198f-4a4b-8c86-c8c8d8c3d76a' class='xr-section-summary-in' type='checkbox'  checked><label for='section-05023f5b-198f-4a4b-8c86-c8c8d8c3d76a' class='xr-section-summary' >Coordinates: <span>(2)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>x</span></div><div class='xr-var-dims'>(x)</div><div class='xr-var-dtype'>float32</div><div class='xr-var-preview xr-preview'>-180.0 -180.0 ... 180.0 180.0</div><input id='attrs-f1f07a8e-4c5b-4403-96ce-e9843c444960' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-f1f07a8e-4c5b-4403-96ce-e9843c444960' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-e0b20733-f570-4ebf-9eda-c2f58620870d' class='xr-var-data-in' type='checkbox'><label for='data-e0b20733-f570-4ebf-9eda-c2f58620870d' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>long_name :</span></dt><dd>x</dd><dt><span>actual_range :</span></dt><dd>[-180.  180.]</dd></dl></div><div class='xr-var-data'><pre>array([-180.     , -179.96666, -179.93333, ...,  179.93333,  179.96666,\n",
       "        180.     ], dtype=float32)</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>y</span></div><div class='xr-var-dims'>(y)</div><div class='xr-var-dtype'>float32</div><div class='xr-var-preview xr-preview'>-90.0 -89.97 -89.93 ... 89.97 90.0</div><input id='attrs-c4e03d5d-7015-4d9c-ad18-6bbd299a1460' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-c4e03d5d-7015-4d9c-ad18-6bbd299a1460' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-95d97ed9-0fd8-41ec-98c8-4b986206aa6c' class='xr-var-data-in' type='checkbox'><label for='data-95d97ed9-0fd8-41ec-98c8-4b986206aa6c' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>long_name :</span></dt><dd>y</dd><dt><span>actual_range :</span></dt><dd>[-90.  90.]</dd></dl></div><div class='xr-var-data'><pre>array([-90.      , -89.96667 , -89.933334, ...,  89.933334,  89.96667 ,\n",
       "        90.      ], dtype=float32)</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-1befdf26-5e58-4bbc-a103-c6381055f061' class='xr-section-summary-in' type='checkbox'  checked><label for='section-1befdf26-5e58-4bbc-a103-c6381055f061' class='xr-section-summary' >Attributes: <span>(2)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><dl class='xr-attrs'><dt><span>long_name :</span></dt><dd>z</dd><dt><span>actual_range :</span></dt><dd>[-388.63299561  990.96801758]</dd></dl></div></li></ul></div></div>"
      ],
      "text/plain": [
       "<xarray.DataArray 'z' (y: 5401, x: 10801)>\n",
       "[58336201 values with dtype=float32]\n",
       "Coordinates:\n",
       "  * x        (x) float32 -180.0 -180.0 -179.9 -179.9 ... 179.9 179.9 180.0 180.0\n",
       "  * y        (y) float32 -90.0 -89.97 -89.93 -89.9 ... 89.9 89.93 89.97 90.0\n",
       "Attributes:\n",
       "    long_name:     z\n",
       "    actual_range:  [-388.63299561  990.96801758]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "da = xr.open_dataarray('WGM2012_Freeair_ponc_2min.grd').squeeze(drop=True)\n",
    "da"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reduce dataset size\n",
    "da = da.coarsen({'y':10, 'x':10}, boundary='trim').mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Raster to Database"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2.04 s, sys: 214 ms, total: 2.26 s\n",
      "Wall time: 2.49 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "da.to_dataframe(name='z').to_csv('data.csv', header=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "fname = !pwd\n",
    "fname = fname[0] + '/data.csv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " * postgresql://localhost:5432/grav_29_1\n",
      "Done.\n",
      "Done.\n",
      "Done.\n",
      "583200 rows affected.\n",
      "583200 rows affected.\n",
      "Done.\n",
      "Done.\n",
      "CPU times: user 15.5 ms, sys: 5.74 ms, total: 21.3 ms\n",
      "Wall time: 10.5 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "%%sql\n",
    "\n",
    "create extension if not exists postgis;\n",
    "drop table if exists data;\n",
    "create table data (lat float, lon float, z float, the_geom geography(Point, 4326));\n",
    "COPY data (lat, lon, z) FROM :fname WITH DELIMITER ',';\n",
    "UPDATE data SET the_geom = ST_SetSRID(ST_MakePoint(lon, lat), 4326);\n",
    "CREATE INDEX data_the_geom_idx ON data USING GIST (the_geom);\n",
    "ANALYZE data;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Grid to Database"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "32\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div><svg style=\"position: absolute; width: 0; height: 0; overflow: hidden\">\n",
       "<defs>\n",
       "<symbol id=\"icon-database\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M16 0c-8.837 0-16 2.239-16 5v4c0 2.761 7.163 5 16 5s16-2.239 16-5v-4c0-2.761-7.163-5-16-5z\"></path>\n",
       "<path d=\"M16 17c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "<path d=\"M16 26c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "</symbol>\n",
       "<symbol id=\"icon-file-text2\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M28.681 7.159c-0.694-0.947-1.662-2.053-2.724-3.116s-2.169-2.030-3.116-2.724c-1.612-1.182-2.393-1.319-2.841-1.319h-15.5c-1.378 0-2.5 1.121-2.5 2.5v27c0 1.378 1.122 2.5 2.5 2.5h23c1.378 0 2.5-1.122 2.5-2.5v-19.5c0-0.448-0.137-1.23-1.319-2.841zM24.543 5.457c0.959 0.959 1.712 1.825 2.268 2.543h-4.811v-4.811c0.718 0.556 1.584 1.309 2.543 2.268zM28 29.5c0 0.271-0.229 0.5-0.5 0.5h-23c-0.271 0-0.5-0.229-0.5-0.5v-27c0-0.271 0.229-0.5 0.5-0.5 0 0 15.499-0 15.5 0v7c0 0.552 0.448 1 1 1h7v19.5z\"></path>\n",
       "<path d=\"M23 26h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 22h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 18h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "</symbol>\n",
       "</defs>\n",
       "</svg>\n",
       "<style>/* CSS stylesheet for displaying xarray objects in jupyterlab.\n",
       " *\n",
       " */\n",
       "\n",
       ":root {\n",
       "  --xr-font-color0: var(--jp-content-font-color0, rgba(0, 0, 0, 1));\n",
       "  --xr-font-color2: var(--jp-content-font-color2, rgba(0, 0, 0, 0.54));\n",
       "  --xr-font-color3: var(--jp-content-font-color3, rgba(0, 0, 0, 0.38));\n",
       "  --xr-border-color: var(--jp-border-color2, #e0e0e0);\n",
       "  --xr-disabled-color: var(--jp-layout-color3, #bdbdbd);\n",
       "  --xr-background-color: var(--jp-layout-color0, white);\n",
       "  --xr-background-color-row-even: var(--jp-layout-color1, white);\n",
       "  --xr-background-color-row-odd: var(--jp-layout-color2, #eeeeee);\n",
       "}\n",
       "\n",
       "html[theme=dark],\n",
       "body.vscode-dark {\n",
       "  --xr-font-color0: rgba(255, 255, 255, 1);\n",
       "  --xr-font-color2: rgba(255, 255, 255, 0.54);\n",
       "  --xr-font-color3: rgba(255, 255, 255, 0.38);\n",
       "  --xr-border-color: #1F1F1F;\n",
       "  --xr-disabled-color: #515151;\n",
       "  --xr-background-color: #111111;\n",
       "  --xr-background-color-row-even: #111111;\n",
       "  --xr-background-color-row-odd: #313131;\n",
       "}\n",
       "\n",
       ".xr-wrap {\n",
       "  display: block;\n",
       "  min-width: 300px;\n",
       "  max-width: 700px;\n",
       "}\n",
       "\n",
       ".xr-text-repr-fallback {\n",
       "  /* fallback to plain text repr when CSS is not injected (untrusted notebook) */\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-header {\n",
       "  padding-top: 6px;\n",
       "  padding-bottom: 6px;\n",
       "  margin-bottom: 4px;\n",
       "  border-bottom: solid 1px var(--xr-border-color);\n",
       "}\n",
       "\n",
       ".xr-header > div,\n",
       ".xr-header > ul {\n",
       "  display: inline;\n",
       "  margin-top: 0;\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-obj-type,\n",
       ".xr-array-name {\n",
       "  margin-left: 2px;\n",
       "  margin-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-obj-type {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-sections {\n",
       "  padding-left: 0 !important;\n",
       "  display: grid;\n",
       "  grid-template-columns: 150px auto auto 1fr 20px 20px;\n",
       "}\n",
       "\n",
       ".xr-section-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-section-item input {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-item input + label {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label {\n",
       "  cursor: pointer;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label:hover {\n",
       "  color: var(--xr-font-color0);\n",
       "}\n",
       "\n",
       ".xr-section-summary {\n",
       "  grid-column: 1;\n",
       "  color: var(--xr-font-color2);\n",
       "  font-weight: 500;\n",
       "}\n",
       "\n",
       ".xr-section-summary > span {\n",
       "  display: inline-block;\n",
       "  padding-left: 0.5em;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in + label:before {\n",
       "  display: inline-block;\n",
       "  content: '►';\n",
       "  font-size: 11px;\n",
       "  width: 15px;\n",
       "  text-align: center;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label:before {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label:before {\n",
       "  content: '▼';\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label > span {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-summary,\n",
       ".xr-section-inline-details {\n",
       "  padding-top: 4px;\n",
       "  padding-bottom: 4px;\n",
       "}\n",
       "\n",
       ".xr-section-inline-details {\n",
       "  grid-column: 2 / -1;\n",
       "}\n",
       "\n",
       ".xr-section-details {\n",
       "  display: none;\n",
       "  grid-column: 1 / -1;\n",
       "  margin-bottom: 5px;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked ~ .xr-section-details {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-array-wrap {\n",
       "  grid-column: 1 / -1;\n",
       "  display: grid;\n",
       "  grid-template-columns: 20px auto;\n",
       "}\n",
       "\n",
       ".xr-array-wrap > label {\n",
       "  grid-column: 1;\n",
       "  vertical-align: top;\n",
       "}\n",
       "\n",
       ".xr-preview {\n",
       "  color: var(--xr-font-color3);\n",
       "}\n",
       "\n",
       ".xr-array-preview,\n",
       ".xr-array-data {\n",
       "  padding: 0 5px !important;\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-array-data,\n",
       ".xr-array-in:checked ~ .xr-array-preview {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-array-in:checked ~ .xr-array-data,\n",
       ".xr-array-preview {\n",
       "  display: inline-block;\n",
       "}\n",
       "\n",
       ".xr-dim-list {\n",
       "  display: inline-block !important;\n",
       "  list-style: none;\n",
       "  padding: 0 !important;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list li {\n",
       "  display: inline-block;\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list:before {\n",
       "  content: '(';\n",
       "}\n",
       "\n",
       ".xr-dim-list:after {\n",
       "  content: ')';\n",
       "}\n",
       "\n",
       ".xr-dim-list li:not(:last-child):after {\n",
       "  content: ',';\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-has-index {\n",
       "  font-weight: bold;\n",
       "}\n",
       "\n",
       ".xr-var-list,\n",
       ".xr-var-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-var-item > div,\n",
       ".xr-var-item label,\n",
       ".xr-var-item > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-even);\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-var-item > .xr-var-name:hover span {\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-var-list > li:nth-child(odd) > div,\n",
       ".xr-var-list > li:nth-child(odd) > label,\n",
       ".xr-var-list > li:nth-child(odd) > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-odd);\n",
       "}\n",
       "\n",
       ".xr-var-name {\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-var-dims {\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-var-dtype {\n",
       "  grid-column: 3;\n",
       "  text-align: right;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-var-preview {\n",
       "  grid-column: 4;\n",
       "}\n",
       "\n",
       ".xr-var-name,\n",
       ".xr-var-dims,\n",
       ".xr-var-dtype,\n",
       ".xr-preview,\n",
       ".xr-attrs dt {\n",
       "  white-space: nowrap;\n",
       "  overflow: hidden;\n",
       "  text-overflow: ellipsis;\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-var-name:hover,\n",
       ".xr-var-dims:hover,\n",
       ".xr-var-dtype:hover,\n",
       ".xr-attrs dt:hover {\n",
       "  overflow: visible;\n",
       "  width: auto;\n",
       "  z-index: 1;\n",
       "}\n",
       "\n",
       ".xr-var-attrs,\n",
       ".xr-var-data {\n",
       "  display: none;\n",
       "  background-color: var(--xr-background-color) !important;\n",
       "  padding-bottom: 5px !important;\n",
       "}\n",
       "\n",
       ".xr-var-attrs-in:checked ~ .xr-var-attrs,\n",
       ".xr-var-data-in:checked ~ .xr-var-data {\n",
       "  display: block;\n",
       "}\n",
       "\n",
       ".xr-var-data > table {\n",
       "  float: right;\n",
       "}\n",
       "\n",
       ".xr-var-name span,\n",
       ".xr-var-data,\n",
       ".xr-attrs {\n",
       "  padding-left: 25px !important;\n",
       "}\n",
       "\n",
       ".xr-attrs,\n",
       ".xr-var-attrs,\n",
       ".xr-var-data {\n",
       "  grid-column: 1 / -1;\n",
       "}\n",
       "\n",
       "dl.xr-attrs {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  display: grid;\n",
       "  grid-template-columns: 125px auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt,\n",
       ".xr-attrs dd {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  float: left;\n",
       "  padding-right: 10px;\n",
       "  width: auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt {\n",
       "  font-weight: normal;\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-attrs dt:hover span {\n",
       "  display: inline-block;\n",
       "  background: var(--xr-background-color);\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-attrs dd {\n",
       "  grid-column: 2;\n",
       "  white-space: pre-wrap;\n",
       "  word-break: break-all;\n",
       "}\n",
       "\n",
       ".xr-icon-database,\n",
       ".xr-icon-file-text2 {\n",
       "  display: inline-block;\n",
       "  vertical-align: middle;\n",
       "  width: 1em;\n",
       "  height: 1.5em !important;\n",
       "  stroke-width: 0;\n",
       "  stroke: currentColor;\n",
       "  fill: currentColor;\n",
       "}\n",
       "</style><pre class='xr-text-repr-fallback'>&lt;xarray.DataArray (y: 4, x: 8)&gt;\n",
       "array([[-31.74633  , -33.00165  , -33.436398 , -33.367577 , -33.757626 ,\n",
       "        -31.9953   , -31.16913  , -31.360115 ],\n",
       "       [ 26.987411 ,  -2.5499501,  -2.23793  ,   4.5407796,   5.9301405,\n",
       "         18.586128 , -20.639273 ,   9.747361 ],\n",
       "       [ -9.16997  , -11.019231 , -82.82903  ,   5.4330297, -12.44292  ,\n",
       "        -45.52271  ,  35.572807 ,  -6.545051 ],\n",
       "       [-48.02255  ,   6.5627704, -75.38486  ,  44.41182  , -51.596848 ,\n",
       "         -9.07185  , -30.33515  ,  88.0656   ]], dtype=float32)\n",
       "Coordinates:\n",
       "  * x        (x) float32 -179.9 -129.9 -79.85 -29.85 20.15 70.15 120.2 170.1\n",
       "  * y        (y) float32 -89.85 -39.85 10.15 60.15</pre><div class='xr-wrap' hidden><div class='xr-header'><div class='xr-obj-type'>xarray.DataArray</div><div class='xr-array-name'></div><ul class='xr-dim-list'><li><span class='xr-has-index'>y</span>: 4</li><li><span class='xr-has-index'>x</span>: 8</li></ul></div><ul class='xr-sections'><li class='xr-section-item'><div class='xr-array-wrap'><input id='section-5fcee954-55d4-47da-9920-5c25e6e35f72' class='xr-array-in' type='checkbox' checked><label for='section-5fcee954-55d4-47da-9920-5c25e6e35f72' title='Show/hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-array-preview xr-preview'><span>-31.75 -33.0 -33.44 -33.37 -33.76 ... 44.41 -51.6 -9.072 -30.34 88.07</span></div><div class='xr-array-data'><pre>array([[-31.74633  , -33.00165  , -33.436398 , -33.367577 , -33.757626 ,\n",
       "        -31.9953   , -31.16913  , -31.360115 ],\n",
       "       [ 26.987411 ,  -2.5499501,  -2.23793  ,   4.5407796,   5.9301405,\n",
       "         18.586128 , -20.639273 ,   9.747361 ],\n",
       "       [ -9.16997  , -11.019231 , -82.82903  ,   5.4330297, -12.44292  ,\n",
       "        -45.52271  ,  35.572807 ,  -6.545051 ],\n",
       "       [-48.02255  ,   6.5627704, -75.38486  ,  44.41182  , -51.596848 ,\n",
       "         -9.07185  , -30.33515  ,  88.0656   ]], dtype=float32)</pre></div></div></li><li class='xr-section-item'><input id='section-16d0ff6a-3a6e-4e1f-a7cf-968312a08967' class='xr-section-summary-in' type='checkbox'  checked><label for='section-16d0ff6a-3a6e-4e1f-a7cf-968312a08967' class='xr-section-summary' >Coordinates: <span>(2)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>x</span></div><div class='xr-var-dims'>(x)</div><div class='xr-var-dtype'>float32</div><div class='xr-var-preview xr-preview'>-179.9 -129.9 ... 120.2 170.1</div><input id='attrs-897879c4-c00b-4bb9-b26d-a8e24c86ff17' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-897879c4-c00b-4bb9-b26d-a8e24c86ff17' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-2a23054b-49a4-4c3e-bc0f-fb4285cc62b0' class='xr-var-data-in' type='checkbox'><label for='data-2a23054b-49a4-4c3e-bc0f-fb4285cc62b0' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([-179.85    , -129.85    ,  -79.850006,  -29.85    ,   20.15    ,\n",
       "         70.149994,  120.15    ,  170.15    ], dtype=float32)</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>y</span></div><div class='xr-var-dims'>(y)</div><div class='xr-var-dtype'>float32</div><div class='xr-var-preview xr-preview'>-89.85 -39.85 10.15 60.15</div><input id='attrs-7a4ce803-9bbb-4864-9fb2-1c998e1e8917' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-7a4ce803-9bbb-4864-9fb2-1c998e1e8917' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-1d4bc59c-e593-4a28-a970-b40c27a6b105' class='xr-var-data-in' type='checkbox'><label for='data-1d4bc59c-e593-4a28-a970-b40c27a6b105' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([-89.850006, -39.85    ,  10.15    ,  60.15    ], dtype=float32)</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-3ae706c6-07e0-43b2-8519-aff42cf41353' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-3ae706c6-07e0-43b2-8519-aff42cf41353' class='xr-section-summary'  title='Expand/collapse section'>Attributes: <span>(0)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><dl class='xr-attrs'></dl></div></li></ul></div></div>"
      ],
      "text/plain": [
       "<xarray.DataArray (y: 4, x: 8)>\n",
       "array([[-31.74633  , -33.00165  , -33.436398 , -33.367577 , -33.757626 ,\n",
       "        -31.9953   , -31.16913  , -31.360115 ],\n",
       "       [ 26.987411 ,  -2.5499501,  -2.23793  ,   4.5407796,   5.9301405,\n",
       "         18.586128 , -20.639273 ,   9.747361 ],\n",
       "       [ -9.16997  , -11.019231 , -82.82903  ,   5.4330297, -12.44292  ,\n",
       "        -45.52271  ,  35.572807 ,  -6.545051 ],\n",
       "       [-48.02255  ,   6.5627704, -75.38486  ,  44.41182  , -51.596848 ,\n",
       "         -9.07185  , -30.33515  ,  88.0656   ]], dtype=float32)\n",
       "Coordinates:\n",
       "  * x        (x) float32 -179.9 -129.9 -79.85 -29.85 20.15 70.15 120.2 170.1\n",
       "  * y        (y) float32 -89.85 -39.85 10.15 60.15"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid = da[::150,::150]\n",
    "print (grid.shape[0]*grid.shape[1])\n",
    "grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid.to_dataframe(name='z0').to_csv('grid.csv', header=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "fname = !pwd\n",
    "fname = fname[0] + '/grid.csv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " * postgresql://localhost:5432/grav_29_1\n",
      "Done.\n",
      "Done.\n",
      "32 rows affected.\n",
      "32 rows affected.\n",
      "Done.\n",
      "Done.\n",
      "CPU times: user 25.5 ms, sys: 6.97 ms, total: 32.5 ms\n",
      "Wall time: 64.1 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "%%sql\n",
    "\n",
    "drop table if exists grid;\n",
    "create table grid (lat float, lon float, z0 float, the_geom geography(Point, 4326));\n",
    "COPY grid (lat, lon, z0) FROM :fname WITH DELIMITER ',';\n",
    "UPDATE grid SET the_geom = ST_SetSRID(ST_MakePoint(lon, lat), 4326);\n",
    "CREATE INDEX grid_the_geom_idx ON grid USING GIST (the_geom);\n",
    "ANALYZE grid;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gaussian Filtering in Database"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " * postgresql://localhost:5432/grav_29_1\n",
      "Done.\n",
      "Done.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%sql\n",
    "\n",
    "DROP FUNCTION IF EXISTS gaussian_transform(sigma float, geom geography);\n",
    "CREATE OR REPLACE FUNCTION gaussian_transform(sigma float, geom geography)\n",
    "    RETURNS TABLE (z float, count bigint)\n",
    "    AS '\n",
    "with weights as (\n",
    "    select\n",
    "        z,\n",
    "        exp(-(pow(ST_Distance($2,the_geom,false),2))/(2*pow($1,2))) as weight\n",
    "    from data\n",
    "    where ST_Buffer($2,4.*$1) && the_geom and ST_DWithin($2, the_geom, 4.*$1)\n",
    ")\n",
    "select\n",
    "    sum(z*weight)/sum(weight) as z,\n",
    "    count(1) as count\n",
    "from weights\n",
    "'\n",
    "LANGUAGE SQL STABLE;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " * postgresql://localhost:5432/grav_29_1\n",
      "Done.\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "Done.\n",
      "sigmakm 25.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 75.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 125.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 175.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 225.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 275.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 325.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 375.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 425.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 475.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 525.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 575.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 625.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 675.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 725.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 775.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 825.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 875.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 925.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 975.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 1025.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 1075.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 1125.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 1175.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 1225.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 1275.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 1325.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 1375.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "sigmakm 1425.0\n",
      " * postgresql://localhost:5432/grav_29_1\n",
      "32 rows affected.\n",
      "CPU times: user 123 ms, sys: 22.1 ms, total: 145 ms\n",
      "Wall time: 3min 2s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "%sql drop table if exists gaussian_transform;\n",
    "%sql create table gaussian_transform (sigmakm int, lat float, lon float, z0 float, z float, count bigint);\n",
    "\n",
    "for sigmakm in sigmaskm:\n",
    "    print (\"sigmakm\", sigmakm)\n",
    "    %sql insert into gaussian_transform \\\n",
    "         select :sigmakm, lat, lon, z0, t.* from grid, gaussian_transform(:sigmakm*1000, the_geom) as t;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Show Gaussian Transform map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " * postgresql://localhost:5432/grav_29_1\n",
      "928 rows affected.\n",
      "Returning data to local variable gaussian\n"
     ]
    }
   ],
   "source": [
    "%sql gaussian << select * from gaussian_transform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "gaussian = gaussian.DataFrame()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><svg style=\"position: absolute; width: 0; height: 0; overflow: hidden\">\n",
       "<defs>\n",
       "<symbol id=\"icon-database\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M16 0c-8.837 0-16 2.239-16 5v4c0 2.761 7.163 5 16 5s16-2.239 16-5v-4c0-2.761-7.163-5-16-5z\"></path>\n",
       "<path d=\"M16 17c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "<path d=\"M16 26c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "</symbol>\n",
       "<symbol id=\"icon-file-text2\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M28.681 7.159c-0.694-0.947-1.662-2.053-2.724-3.116s-2.169-2.030-3.116-2.724c-1.612-1.182-2.393-1.319-2.841-1.319h-15.5c-1.378 0-2.5 1.121-2.5 2.5v27c0 1.378 1.122 2.5 2.5 2.5h23c1.378 0 2.5-1.122 2.5-2.5v-19.5c0-0.448-0.137-1.23-1.319-2.841zM24.543 5.457c0.959 0.959 1.712 1.825 2.268 2.543h-4.811v-4.811c0.718 0.556 1.584 1.309 2.543 2.268zM28 29.5c0 0.271-0.229 0.5-0.5 0.5h-23c-0.271 0-0.5-0.229-0.5-0.5v-27c0-0.271 0.229-0.5 0.5-0.5 0 0 15.499-0 15.5 0v7c0 0.552 0.448 1 1 1h7v19.5z\"></path>\n",
       "<path d=\"M23 26h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 22h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 18h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "</symbol>\n",
       "</defs>\n",
       "</svg>\n",
       "<style>/* CSS stylesheet for displaying xarray objects in jupyterlab.\n",
       " *\n",
       " */\n",
       "\n",
       ":root {\n",
       "  --xr-font-color0: var(--jp-content-font-color0, rgba(0, 0, 0, 1));\n",
       "  --xr-font-color2: var(--jp-content-font-color2, rgba(0, 0, 0, 0.54));\n",
       "  --xr-font-color3: var(--jp-content-font-color3, rgba(0, 0, 0, 0.38));\n",
       "  --xr-border-color: var(--jp-border-color2, #e0e0e0);\n",
       "  --xr-disabled-color: var(--jp-layout-color3, #bdbdbd);\n",
       "  --xr-background-color: var(--jp-layout-color0, white);\n",
       "  --xr-background-color-row-even: var(--jp-layout-color1, white);\n",
       "  --xr-background-color-row-odd: var(--jp-layout-color2, #eeeeee);\n",
       "}\n",
       "\n",
       "html[theme=dark],\n",
       "body.vscode-dark {\n",
       "  --xr-font-color0: rgba(255, 255, 255, 1);\n",
       "  --xr-font-color2: rgba(255, 255, 255, 0.54);\n",
       "  --xr-font-color3: rgba(255, 255, 255, 0.38);\n",
       "  --xr-border-color: #1F1F1F;\n",
       "  --xr-disabled-color: #515151;\n",
       "  --xr-background-color: #111111;\n",
       "  --xr-background-color-row-even: #111111;\n",
       "  --xr-background-color-row-odd: #313131;\n",
       "}\n",
       "\n",
       ".xr-wrap {\n",
       "  display: block;\n",
       "  min-width: 300px;\n",
       "  max-width: 700px;\n",
       "}\n",
       "\n",
       ".xr-text-repr-fallback {\n",
       "  /* fallback to plain text repr when CSS is not injected (untrusted notebook) */\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-header {\n",
       "  padding-top: 6px;\n",
       "  padding-bottom: 6px;\n",
       "  margin-bottom: 4px;\n",
       "  border-bottom: solid 1px var(--xr-border-color);\n",
       "}\n",
       "\n",
       ".xr-header > div,\n",
       ".xr-header > ul {\n",
       "  display: inline;\n",
       "  margin-top: 0;\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-obj-type,\n",
       ".xr-array-name {\n",
       "  margin-left: 2px;\n",
       "  margin-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-obj-type {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-sections {\n",
       "  padding-left: 0 !important;\n",
       "  display: grid;\n",
       "  grid-template-columns: 150px auto auto 1fr 20px 20px;\n",
       "}\n",
       "\n",
       ".xr-section-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-section-item input {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-item input + label {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label {\n",
       "  cursor: pointer;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label:hover {\n",
       "  color: var(--xr-font-color0);\n",
       "}\n",
       "\n",
       ".xr-section-summary {\n",
       "  grid-column: 1;\n",
       "  color: var(--xr-font-color2);\n",
       "  font-weight: 500;\n",
       "}\n",
       "\n",
       ".xr-section-summary > span {\n",
       "  display: inline-block;\n",
       "  padding-left: 0.5em;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in + label:before {\n",
       "  display: inline-block;\n",
       "  content: '►';\n",
       "  font-size: 11px;\n",
       "  width: 15px;\n",
       "  text-align: center;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label:before {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label:before {\n",
       "  content: '▼';\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label > span {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-summary,\n",
       ".xr-section-inline-details {\n",
       "  padding-top: 4px;\n",
       "  padding-bottom: 4px;\n",
       "}\n",
       "\n",
       ".xr-section-inline-details {\n",
       "  grid-column: 2 / -1;\n",
       "}\n",
       "\n",
       ".xr-section-details {\n",
       "  display: none;\n",
       "  grid-column: 1 / -1;\n",
       "  margin-bottom: 5px;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked ~ .xr-section-details {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-array-wrap {\n",
       "  grid-column: 1 / -1;\n",
       "  display: grid;\n",
       "  grid-template-columns: 20px auto;\n",
       "}\n",
       "\n",
       ".xr-array-wrap > label {\n",
       "  grid-column: 1;\n",
       "  vertical-align: top;\n",
       "}\n",
       "\n",
       ".xr-preview {\n",
       "  color: var(--xr-font-color3);\n",
       "}\n",
       "\n",
       ".xr-array-preview,\n",
       ".xr-array-data {\n",
       "  padding: 0 5px !important;\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-array-data,\n",
       ".xr-array-in:checked ~ .xr-array-preview {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-array-in:checked ~ .xr-array-data,\n",
       ".xr-array-preview {\n",
       "  display: inline-block;\n",
       "}\n",
       "\n",
       ".xr-dim-list {\n",
       "  display: inline-block !important;\n",
       "  list-style: none;\n",
       "  padding: 0 !important;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list li {\n",
       "  display: inline-block;\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list:before {\n",
       "  content: '(';\n",
       "}\n",
       "\n",
       ".xr-dim-list:after {\n",
       "  content: ')';\n",
       "}\n",
       "\n",
       ".xr-dim-list li:not(:last-child):after {\n",
       "  content: ',';\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-has-index {\n",
       "  font-weight: bold;\n",
       "}\n",
       "\n",
       ".xr-var-list,\n",
       ".xr-var-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-var-item > div,\n",
       ".xr-var-item label,\n",
       ".xr-var-item > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-even);\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-var-item > .xr-var-name:hover span {\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-var-list > li:nth-child(odd) > div,\n",
       ".xr-var-list > li:nth-child(odd) > label,\n",
       ".xr-var-list > li:nth-child(odd) > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-odd);\n",
       "}\n",
       "\n",
       ".xr-var-name {\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-var-dims {\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-var-dtype {\n",
       "  grid-column: 3;\n",
       "  text-align: right;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-var-preview {\n",
       "  grid-column: 4;\n",
       "}\n",
       "\n",
       ".xr-var-name,\n",
       ".xr-var-dims,\n",
       ".xr-var-dtype,\n",
       ".xr-preview,\n",
       ".xr-attrs dt {\n",
       "  white-space: nowrap;\n",
       "  overflow: hidden;\n",
       "  text-overflow: ellipsis;\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-var-name:hover,\n",
       ".xr-var-dims:hover,\n",
       ".xr-var-dtype:hover,\n",
       ".xr-attrs dt:hover {\n",
       "  overflow: visible;\n",
       "  width: auto;\n",
       "  z-index: 1;\n",
       "}\n",
       "\n",
       ".xr-var-attrs,\n",
       ".xr-var-data {\n",
       "  display: none;\n",
       "  background-color: var(--xr-background-color) !important;\n",
       "  padding-bottom: 5px !important;\n",
       "}\n",
       "\n",
       ".xr-var-attrs-in:checked ~ .xr-var-attrs,\n",
       ".xr-var-data-in:checked ~ .xr-var-data {\n",
       "  display: block;\n",
       "}\n",
       "\n",
       ".xr-var-data > table {\n",
       "  float: right;\n",
       "}\n",
       "\n",
       ".xr-var-name span,\n",
       ".xr-var-data,\n",
       ".xr-attrs {\n",
       "  padding-left: 25px !important;\n",
       "}\n",
       "\n",
       ".xr-attrs,\n",
       ".xr-var-attrs,\n",
       ".xr-var-data {\n",
       "  grid-column: 1 / -1;\n",
       "}\n",
       "\n",
       "dl.xr-attrs {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  display: grid;\n",
       "  grid-template-columns: 125px auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt,\n",
       ".xr-attrs dd {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  float: left;\n",
       "  padding-right: 10px;\n",
       "  width: auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt {\n",
       "  font-weight: normal;\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-attrs dt:hover span {\n",
       "  display: inline-block;\n",
       "  background: var(--xr-background-color);\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-attrs dd {\n",
       "  grid-column: 2;\n",
       "  white-space: pre-wrap;\n",
       "  word-break: break-all;\n",
       "}\n",
       "\n",
       ".xr-icon-database,\n",
       ".xr-icon-file-text2 {\n",
       "  display: inline-block;\n",
       "  vertical-align: middle;\n",
       "  width: 1em;\n",
       "  height: 1.5em !important;\n",
       "  stroke-width: 0;\n",
       "  stroke: currentColor;\n",
       "  fill: currentColor;\n",
       "}\n",
       "</style><pre class='xr-text-repr-fallback'>&lt;xarray.Dataset&gt;\n",
       "Dimensions:  (lat: 4, lon: 8, sigmakm: 29)\n",
       "Coordinates:\n",
       "  * sigmakm  (sigmakm) int64 25 75 125 175 225 275 ... 1225 1275 1325 1375 1425\n",
       "  * lat      (lat) float64 -89.85 -39.85 10.15 60.15\n",
       "  * lon      (lon) float64 -179.9 -129.9 -79.85 -29.85 20.15 70.15 120.2 170.1\n",
       "Data variables:\n",
       "    z0       (sigmakm, lat, lon) float64 -31.75 -33.0 -33.44 ... -30.34 88.07\n",
       "    z        (sigmakm, lat, lon) float64 -32.06 -32.56 -32.95 ... -2.269 2.448\n",
       "    count    (sigmakm, lat, lon) int64 2863 2863 2863 ... 151522 151661 151402</pre><div class='xr-wrap' hidden><div class='xr-header'><div class='xr-obj-type'>xarray.Dataset</div></div><ul class='xr-sections'><li class='xr-section-item'><input id='section-0bfc8345-8479-4df3-8a21-41c0028bc61f' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-0bfc8345-8479-4df3-8a21-41c0028bc61f' class='xr-section-summary'  title='Expand/collapse section'>Dimensions:</label><div class='xr-section-inline-details'><ul class='xr-dim-list'><li><span class='xr-has-index'>lat</span>: 4</li><li><span class='xr-has-index'>lon</span>: 8</li><li><span class='xr-has-index'>sigmakm</span>: 29</li></ul></div><div class='xr-section-details'></div></li><li class='xr-section-item'><input id='section-b0bc0f56-61cb-4ac1-a958-88307f2bb42a' class='xr-section-summary-in' type='checkbox'  checked><label for='section-b0bc0f56-61cb-4ac1-a958-88307f2bb42a' class='xr-section-summary' >Coordinates: <span>(3)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>sigmakm</span></div><div class='xr-var-dims'>(sigmakm)</div><div class='xr-var-dtype'>int64</div><div class='xr-var-preview xr-preview'>25 75 125 175 ... 1325 1375 1425</div><input id='attrs-191af4b4-bfe4-4943-897c-dacf57b9fb75' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-191af4b4-bfe4-4943-897c-dacf57b9fb75' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-c057afb7-60b4-4195-a24a-3e66a2570e3b' class='xr-var-data-in' type='checkbox'><label for='data-c057afb7-60b4-4195-a24a-3e66a2570e3b' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([  25,   75,  125,  175,  225,  275,  325,  375,  425,  475,  525,  575,\n",
       "        625,  675,  725,  775,  825,  875,  925,  975, 1025, 1075, 1125, 1175,\n",
       "       1225, 1275, 1325, 1375, 1425])</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>lat</span></div><div class='xr-var-dims'>(lat)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>-89.85 -39.85 10.15 60.15</div><input id='attrs-73ee8585-6e6f-4591-babc-9f068a0070c1' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-73ee8585-6e6f-4591-babc-9f068a0070c1' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-3a913352-7d64-4a85-ba91-115cc6f78049' class='xr-var-data-in' type='checkbox'><label for='data-3a913352-7d64-4a85-ba91-115cc6f78049' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([-89.850006, -39.849998,  10.15    ,  60.150002])</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>lon</span></div><div class='xr-var-dims'>(lon)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>-179.9 -129.9 ... 120.2 170.1</div><input id='attrs-0e908cbf-e601-4ea6-957c-5868d749fc58' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-0e908cbf-e601-4ea6-957c-5868d749fc58' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-f29426a2-c6d1-4ac1-98d1-4a32807c3f91' class='xr-var-data-in' type='checkbox'><label for='data-f29426a2-c6d1-4ac1-98d1-4a32807c3f91' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([-179.850006, -129.850006,  -79.850006,  -29.85    ,   20.15    ,\n",
       "         70.149994,  120.150002,  170.149994])</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-d4f53dd3-b111-45fc-a16e-1ceaf644f612' class='xr-section-summary-in' type='checkbox'  checked><label for='section-d4f53dd3-b111-45fc-a16e-1ceaf644f612' class='xr-section-summary' >Data variables: <span>(3)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span>z0</span></div><div class='xr-var-dims'>(sigmakm, lat, lon)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>-31.75 -33.0 ... -30.34 88.07</div><input id='attrs-4ff176f1-976a-4410-b4b9-a11e22171f21' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-4ff176f1-976a-4410-b4b9-a11e22171f21' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-4ba1e0a3-e47d-4ec4-89a5-7d600de3b50e' class='xr-var-data-in' type='checkbox'><label for='data-4ba1e0a3-e47d-4ec4-89a5-7d600de3b50e' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([[[-31.74633  , -33.00165  , -33.436398 , -33.367577 ,\n",
       "         -33.757626 , -31.9953   , -31.16913  , -31.360115 ],\n",
       "        [ 26.987411 ,  -2.5499501,  -2.23793  ,   4.5407796,\n",
       "           5.9301405,  18.586128 , -20.639273 ,   9.747361 ],\n",
       "        [ -9.16997  , -11.019231 , -82.82903  ,   5.4330297,\n",
       "         -12.44292  , -45.52271  ,  35.572807 ,  -6.545051 ],\n",
       "        [-48.02255  ,   6.5627704, -75.38486  ,  44.41182  ,\n",
       "         -51.596848 ,  -9.07185  , -30.33515  ,  88.0656   ]],\n",
       "\n",
       "       [[-31.74633  , -33.00165  , -33.436398 , -33.367577 ,\n",
       "         -33.757626 , -31.9953   , -31.16913  , -31.360115 ],\n",
       "        [ 26.987411 ,  -2.5499501,  -2.23793  ,   4.5407796,\n",
       "           5.9301405,  18.586128 , -20.639273 ,   9.747361 ],\n",
       "        [ -9.16997  , -11.019231 , -82.82903  ,   5.4330297,\n",
       "         -12.44292  , -45.52271  ,  35.572807 ,  -6.545051 ],\n",
       "        [-48.02255  ,   6.5627704, -75.38486  ,  44.41182  ,\n",
       "         -51.596848 ,  -9.07185  , -30.33515  ,  88.0656   ]],\n",
       "\n",
       "       [[-31.74633  , -33.00165  , -33.436398 , -33.367577 ,\n",
       "         -33.757626 , -31.9953   , -31.16913  , -31.360115 ],\n",
       "...\n",
       "        [-48.02255  ,   6.5627704, -75.38486  ,  44.41182  ,\n",
       "         -51.596848 ,  -9.07185  , -30.33515  ,  88.0656   ]],\n",
       "\n",
       "       [[-31.74633  , -33.00165  , -33.436398 , -33.367577 ,\n",
       "         -33.757626 , -31.9953   , -31.16913  , -31.360115 ],\n",
       "        [ 26.987411 ,  -2.5499501,  -2.23793  ,   4.5407796,\n",
       "           5.9301405,  18.586128 , -20.639273 ,   9.747361 ],\n",
       "        [ -9.16997  , -11.019231 , -82.82903  ,   5.4330297,\n",
       "         -12.44292  , -45.52271  ,  35.572807 ,  -6.545051 ],\n",
       "        [-48.02255  ,   6.5627704, -75.38486  ,  44.41182  ,\n",
       "         -51.596848 ,  -9.07185  , -30.33515  ,  88.0656   ]],\n",
       "\n",
       "       [[-31.74633  , -33.00165  , -33.436398 , -33.367577 ,\n",
       "         -33.757626 , -31.9953   , -31.16913  , -31.360115 ],\n",
       "        [ 26.987411 ,  -2.5499501,  -2.23793  ,   4.5407796,\n",
       "           5.9301405,  18.586128 , -20.639273 ,   9.747361 ],\n",
       "        [ -9.16997  , -11.019231 , -82.82903  ,   5.4330297,\n",
       "         -12.44292  , -45.52271  ,  35.572807 ,  -6.545051 ],\n",
       "        [-48.02255  ,   6.5627704, -75.38486  ,  44.41182  ,\n",
       "         -51.596848 ,  -9.07185  , -30.33515  ,  88.0656   ]]])</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>z</span></div><div class='xr-var-dims'>(sigmakm, lat, lon)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>-32.06 -32.56 ... -2.269 2.448</div><input id='attrs-17dd60d0-009b-44ab-989b-3d08f56f3073' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-17dd60d0-009b-44ab-989b-3d08f56f3073' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-c9ab56a4-6c76-4228-91f8-2793cdba477a' class='xr-var-data-in' type='checkbox'><label for='data-c9ab56a4-6c76-4228-91f8-2793cdba477a' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([[[-3.20607800e+01, -3.25610955e+01, -3.29491505e+01,\n",
       "         -3.31758585e+01, -3.30798394e+01, -3.24690730e+01,\n",
       "         -3.19418572e+01, -3.19866600e+01],\n",
       "        [ 2.96603334e+01, -1.44536877e+00, -1.00078404e+00,\n",
       "          3.25220169e+00,  2.91973510e+00,  1.48966582e+01,\n",
       "         -2.02865320e+01,  1.17390327e+01],\n",
       "        [-9.45877030e+00, -1.00285340e+01, -6.81756173e+01,\n",
       "         -1.97847072e+00, -1.24919366e+01, -4.39421400e+01,\n",
       "          3.71994103e+01, -8.35049082e-01],\n",
       "        [-3.00424710e+01,  1.78676267e+00, -6.62742169e+01,\n",
       "          4.87010105e+01, -5.02308778e+01, -1.18618227e+01,\n",
       "         -2.50039960e+01,  5.55986253e+01]],\n",
       "\n",
       "       [[-3.06622320e+01, -3.09007055e+01, -3.15965667e+01,\n",
       "         -3.22981795e+01, -3.24525595e+01, -3.19220243e+01,\n",
       "         -3.11633854e+01, -3.06902793e+01],\n",
       "        [ 2.52455421e+01, -1.39391814e+00, -7.31965586e-01,\n",
       "          4.65412272e-01,  1.82215457e+00,  1.38442847e+01,\n",
       "         -1.91708387e+01,  1.22357630e+01],\n",
       "        [-8.59705507e+00, -9.30841626e+00, -1.55672036e+01,\n",
       "...\n",
       "         -8.99126256e+00,  1.68733820e+00],\n",
       "        [-1.41575411e+00, -5.32024292e+00,  3.03201959e+00,\n",
       "         -1.24870506e+00, -4.66855250e-01, -1.72205421e+01,\n",
       "          1.29057632e+01, -7.42091923e-01],\n",
       "        [ 2.72660422e+00, -2.77179042e-01, -5.84946092e+00,\n",
       "          1.46227547e+01,  7.08638363e+00, -4.37036394e+00,\n",
       "         -2.43973864e+00,  2.54052491e+00]],\n",
       "\n",
       "       [[-9.96980174e+00, -9.96545742e+00, -9.92117318e+00,\n",
       "         -9.86859303e+00, -9.84491582e+00, -9.86716233e+00,\n",
       "         -9.91975865e+00, -9.96494412e+00],\n",
       "        [ 1.52687098e+00, -1.66300823e+00,  3.44717490e+00,\n",
       "         -3.13807429e-01,  6.46087919e+00,  1.05170127e+01,\n",
       "         -8.57591065e+00,  1.38048446e+00],\n",
       "        [-1.25183114e+00, -5.12180469e+00,  2.56821480e+00,\n",
       "         -1.28335572e+00, -3.24359566e-01, -1.66151671e+01,\n",
       "          1.23939185e+01, -5.31526375e-01],\n",
       "        [ 2.61384637e+00, -3.44876389e-01, -5.09536311e+00,\n",
       "          1.39701196e+01,  7.14678415e+00, -4.05531825e+00,\n",
       "         -2.26854826e+00,  2.44847429e+00]]])</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>count</span></div><div class='xr-var-dims'>(sigmakm, lat, lon)</div><div class='xr-var-dtype'>int64</div><div class='xr-var-preview xr-preview'>2863 2863 2863 ... 151661 151402</div><input id='attrs-476a640d-a320-4f71-98e6-da81514cf89e' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-476a640d-a320-4f71-98e6-da81514cf89e' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-2dd2548e-4a4c-4e41-970c-cdb6f9cfbd34' class='xr-var-data-in' type='checkbox'><label for='data-2dd2548e-4a4c-4e41-970c-cdb6f9cfbd34' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([[[  2863,   2863,   2863,   2863,   2863,   2863,   2863,   2863],\n",
       "        [    31,     31,     31,     31,     31,     31,     31,     31],\n",
       "        [    21,     21,     21,     21,     21,     21,     21,     21],\n",
       "        [    45,     45,     45,     45,     45,     45,     45,     45]],\n",
       "\n",
       "       [[  8811,   8811,   8811,   8811,   8811,   8811,   8811,   8811],\n",
       "        [   271,    271,    271,    271,    271,    271,    271,    271],\n",
       "        [   213,    213,    213,    213,    213,    213,    213,    213],\n",
       "        [   415,    415,    415,    415,    415,    415,    415,    415]],\n",
       "\n",
       "       [[ 14559,  14559,  14559,  14559,  14559,  14559,  14559,  14559],\n",
       "        [   743,    743,    743,    741,    743,    743,    743,    743],\n",
       "        [   587,    587,    587,    586,    587,    587,    586,    587],\n",
       "        [  1149,   1149,   1148,   1149,   1149,   1149,   1149,   1149]],\n",
       "\n",
       "       [[ 20287,  20287,  20287,  20287,  20287,  20287,  20287,  20287],\n",
       "        [  1457,   1455,   1456,   1456,   1457,   1457,   1457,   1456],\n",
       "        [  1139,   1138,   1139,   1136,   1138,   1138,   1135,   1139],\n",
       "        [  2239,   2239,   2237,   2239,   2239,   2239,   2238,   2237]],\n",
       "\n",
       "...\n",
       "\n",
       "       [[148333, 148333, 148333, 148333, 148333, 148333, 148333, 148333],\n",
       "        [ 86003,  87946,  86042,  87946,  85984,  85986,  87946,  86050],\n",
       "        [ 60978,  60312,  60546,  60131,  60249,  60252,  60262,  60556],\n",
       "        [132169, 132411, 132238, 132327, 132293, 132296, 132304, 132240]],\n",
       "\n",
       "       [[154049, 154049, 154049, 154049, 154049, 154049, 154049, 154049],\n",
       "        [ 94772,  96960,  94750,  96960,  96960,  96960,  96960,  94761],\n",
       "        [ 65900,  65167,  65418,  64930,  65066,  65075,  65057,  65435],\n",
       "        [138656, 138956, 138720, 138867, 138783, 138789, 138825, 138728]],\n",
       "\n",
       "       [[159995, 159995, 159995, 159995, 159995, 159995, 159995, 159995],\n",
       "        [104937, 107568, 104930, 107568, 107568, 107568, 107568, 104942],\n",
       "        [ 71028,  70145,  70474,  69915,  70073,  70087,  70045,  70507],\n",
       "        [145025, 145424, 145103, 145306, 145189, 145177, 145256, 145102]],\n",
       "\n",
       "       [[165777, 165777, 165777, 165777, 165777, 165777, 165777, 165777],\n",
       "        [121344, 121344, 121344, 121344, 121344, 121344, 121344, 121344],\n",
       "        [ 76377,  75341,  75745,  75128,  75306,  75318,  75253,  75799],\n",
       "        [151309, 151827, 151401, 151697, 151525, 151522, 151661, 151402]]])</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-efc72820-fd07-4072-96a7-20562373b480' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-efc72820-fd07-4072-96a7-20562373b480' class='xr-section-summary'  title='Expand/collapse section'>Attributes: <span>(0)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><dl class='xr-attrs'></dl></div></li></ul></div></div>"
      ],
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:  (lat: 4, lon: 8, sigmakm: 29)\n",
       "Coordinates:\n",
       "  * sigmakm  (sigmakm) int64 25 75 125 175 225 275 ... 1225 1275 1325 1375 1425\n",
       "  * lat      (lat) float64 -89.85 -39.85 10.15 60.15\n",
       "  * lon      (lon) float64 -179.9 -129.9 -79.85 -29.85 20.15 70.15 120.2 170.1\n",
       "Data variables:\n",
       "    z0       (sigmakm, lat, lon) float64 -31.75 -33.0 -33.44 ... -30.34 88.07\n",
       "    z        (sigmakm, lat, lon) float64 -32.06 -32.56 -32.95 ... -2.269 2.448\n",
       "    count    (sigmakm, lat, lon) int64 2863 2863 2863 ... 151522 151661 151402"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gaussian = gaussian.set_index(['sigmakm','lat','lon']).to_xarray()\n",
    "gaussian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.42 ms, sys: 260 µs, total: 1.68 ms\n",
      "Wall time: 1.56 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# discrete power spectrum\n",
    "power_spectrum = gaussian.std(['lat','lon'])['z'].rename('spectrum')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABFYAAAGqCAYAAAAhqMCjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsnXe4HVXVh99FQgmh9yIQBekgCEg1BKnSREEUEIxiQYpKExALoqgo+EkTRcQAUkSUXiMmIAhCiEhv0kLvJUCAJOv7Y+3hzJ3MmVPuuffcm/ze55lnzpnZe8+aPbPLrL322ubuCCGEEEIIIYQQQojWma3bAgghhBBCCCGEEEIMVqRYEUIIIYQQQgghhGgTKVaEEEIIIYQQQggh2kSKFSGEEEIIIYQQQog2kWJFCCGEEEIIIYQQok2kWBFCCCGEEEIIIYRoEylWhBCzLGb2rJm5mf2227LMaijvhRBiRszswlQ33t1tWcSMmNku6fm4ma1ecv7kdO7FbsgnOoeZHZJ71vN0Wx4x8JFipUnMbDYzeyUVrnubCH9XrjAe3SDsSrmwf2oQdriZfcnM/mRm95nZi2b2XpLtwdQgH2RmyzZI57HcNd3MvtbonlK8IwvxbqkTbn4z293M/mBm/zGzV5OcL5nZjWb2PTNbrJlrpvTmMrMDzeyWlMZbZvaQmZ1iZis3EX8eMxuZKskLzOzR3D081oIc6yTZrzazSWb2jpm9aWb/M7NzzWybZtNqcJ0xhXxutI3oxHUHK4Xy9pduy9Ms6b3cO70796d3OyvPD5jZX8zs22b2gW7LOqtjZufXKXvTU/12t5mdZmbrd1vWmREzW8TMvmNm1yWl3DtmNiX9viXl/ZfMbKluyypEN0ntSiv9h8u7LbPoHtZTUZTfpprZy6m/fL2Z/crMdjazObotsxADFnfX1uQGXAZ42hapCLcQMD0XdnyDdL+aC/u1inDfBJ7Pha3apgOXAKvVSeuxQvh/NpkHDxTi3VIS5pPAlCZkfBXYuYlrfhC4tyKdt4HRDdIYVxH/sSbv/fom8/5yYP5evmtjmrxWto3odvno1gasU8iLKcCCTcZ9NsX5bT/LbMDBwIstlOe/Ait3O787mAd18x5YOXfvn++2rEmm81soj6cAs3Vb5pllA3YGXm4y7x/utrx9nBfHpfuc3G1ZtPXZM74wPeO724w/T4v9h8u7fc+9zK91c/eyfT9cb5fc9VYvOX9yOvdinfh3p/MXdjvvSu6nme154PvA7N2WvR/y5pDcfc9Tcv7ydG5Ct2XVNjC2oYhWuAHYPv0eCfytTriPEx9O04AhwPpmNoe7v1sn/Mjc7+uLJ81sTuAsYNd0aBpwKTAWeJDocM4HLAVsmmRcGtgReAbYp+KepgBzARub2Qfd/dF6Ac1sA2DFQrwyFgbmJCqbfwBXAxOBV4AlgM8AXwbmB/5sZtu5+zV1rjkPcAWwSjo0BjgHmAxsBHw3Xe90M3va3a+tJ37u98vAhBS/FdO+pdP+OaLjcwPweLrPdYFvAx8GtgMuNbPN3H16C+nXY2vg6QZhnurAdQYrX0z7N4HhxLv3OaDhFBN3X6IP5SrFzIYB5wI7pUNTgYuBvwMPEeVkPuJ9G0WU5yWJcjOJeM8GPd3I+w4yCngp/TZgMWAb4ADi/dsXeAE4qguyzVSY2VbAn4m2dCpRdi4nBgemEXm/NrAV0fYKIWrcAOzXIMzr/SHIrIK77w/s32052uSnwHm5//MACwJrAp8AtgQWBY4GdjSz7d39uX6Xsp9w9+MIZbYQTSHFSmvklR6NFCsAVwKbERXTx4AbG4R/zt0fKDl/CjWlykRgN3d/sE5a55nZ7MAewE/qhMlzJ9ExHQHsSVSW9dgr7W8hPvSWqxPuPeD3wE/d/bGS81eZ2VWEcmIIcLKZrejuXhL2YGpKle+5+zG5c7eY2ZWEkmQ4cJKZrebuU0vSORf4HXCbuz8MMR2K1hQr9wPfI0YZite41czOBK4lFDYjgd2ByqldTfJgnXyc5Unv+m7p7/nEu7IRoWwZqL47TqOmVLkV2N3d/1cn7LnJ7HZPmivPon94wN2fLRy7zswuI5TJQ4BDzexX7q6Plt7xayI/3wW2cPd/loS5GviZmS0JfKo/hRNigPOGu8tXi2iWZ+q8L1cBx5rZKsCZwHrEgOIlZjbK3af0p5BCDFTkY6U1JhKWEtDTyqRIdm4coYSoG97MlqGmoLih5PyngL3T37uAkRVKFQDc/T13HwOsQSh3KoNT+/jfs16g9HH3+fT37AbX/7O7f61KGeDufwMuSn9XANYquebswLfS3weBn5Wkcz9wbPq7InU61e5+mruflylV2sHdt3f38+sobnD3N+lpHfTZdq8lmmY7YJH0+2xq7+YGZrZieZTuYWafA76Q/t4OjKpQqgDg7u+6+x+I8jy2j0UUvcDdbyCsjwDmJiwIRZukMpwp1s+to1R5H3d/xt0HqkJVCCEGNe5+HzEYnH3brA8c2D2JhBhYSLHSAumD+ub09yNmNl8xTJq6snb6+09qVir1FDH5jvcMihXCQiLji+njvVl5X3H3S5sIelbar2BmG9YJswNhDvguYRnQCcblfq9Qcn5UuibAmIppNWfkfn+6A3K1jbvfRW2KQNk99QvJmaOb2dXp/0oWnuofSs523Uqc/prZ8slB2Z3JIecUC0fH55jZJk1eewkz+4mZ3Zacsb5jZk+Z2d/MbPvGKbRENg3ocaL8XEC8o/lzVbLWXZnGzLbJOXHbwMKB9d5mNs7MnktOS5suC2ZmwJHprwN7ufvbzcZ39xfd/YqSdFfOyfn5dGwXM7sy5ftUyzmZtmAjM/tpckj3nJm9a2avmdk96T1Zqc49rJi71jFlYQrhP5EL/7XCuRny3sJJtQP35YKeV+JU7/AU/p70/44mZJnbzF5P4S9oFL4X5B16l1r1mdmaFs5Ws/I42cJ58cn1FIJmNip3/1vUCXN5LkyphZOZHW41x4Tz1wmzpJkdY2YTLJwXvmNmT1o4R9+23o1bT6eZh6RjW5nZXy2cfb9nZk/Wi19C3gl7pQKyERaO3t3MTk7/NzSzPye5pqT7O8tKVvmok15beVRIY14zO9jCIe8zqRy+YWb/NbNTzewTubD7p7JxcDo0vKRcuJktkoszwXLOSc1sFTP7jfVsB1bIp19Mo0TmGZ5x4fxx6dzk9H9Bi/bg3nTNZ83sUjNbtxBvKTM71sIp/1tm9kLKy9WaycsqzGxtM/uhmf3dzJ7O5fODFg72P9ogfvGe5jazwywc87+e0ppg4Wh89ibk+YjFAgRPpXfvifTufaS399oJzGzd3DPePh3b3cyuSe/pNDP7ey68WSwOcKyZ/dPMnk9l/TWLvsSvzexDLVx/RwuH7o+kd+Ydi37I3yza4HlSuEVSmbgtF/2ykjKxfyH9lc3sCDO7ymrlP1uE4E9mVjVw2oz8pasCWaqfgeyd3rlE1gkp7NG5Yx9u4pp/TmFfN7PhvZG/Ee7+DjFANC0dOtBiinM92RayWPThZot6+N30Hl1mZruamVXEnZzu67j0f02LxR0eT+/Fc+m92KBKZjOb3cy+bmZjrWef5yEz+0eSb42SeKWrAmXPmBjcA1in5Fm+mMJ+OXdsyyo5U/hjreYYf0Sj8GKA0W0nL4NtIz6MMkdG25ac3yqde4MwX/5E+v86MKQk/Gm59NYonFs9d258h+/jsZTuLen/zen/qXXCX5LO/60sfpsyHJS7v8+UnP9R7vxGDdJ6MIV7rI08aDpOk+m+ltK9sxdpjMnd+4g24t+S4l5NOH58M5detq1ciHMooZSoclp2IhVOOYnGdnKDNP4CDOtAPi+Sk/enueMXpWNPVMmawlY5UN0mJ/MnKHeAfH4L8uYd7F3dwfct7+h1d2LaW1HOW3Lh92nwfJyYzrd3nevdmsI8ClgD2f6Qws7gULgs7wm/TY1kc+DwFP7A3LF1GsjyxVzYrdvI57zz2iUqwn0rF+5bJed/QHRIq/L+gJJ4c1FzCv6TkvNDqNU9DtxYR75rqHC2l/KprL7o8d4Dc5XEzTvNPAQ4viTuky3k+cdz8c7pZTnJHEWfTPicmFrn3t4llJ6N3qW28iiXxlaEH55G7/oiKfz+TZaNRXLXmJCOXU5MJ36rJPwKJelXOefv8YxLzr/vXJcYXHikIp+3S3HWI3yXlYV7E9iwF899+ybybDox1bheGvl7+gA156Nl21gqnHoCo6nfzk4h6vBOOq9t2TEtPduqTxNWeEVZ/54Lf0id+5nh3hpcd3FigKRRWvun8Is0Efb98CX3VrWdWCFnW85rqTk6rdompLDLUWsnftYg3xai1jb8vo3nnb+f/VuId1Eu3pYV5e+VBvd8DTBfnfhZf/I4ouy8UyeNqcAeddJYmFpdWLXNUFao47w294yrthdT2OHEN6AD5zXI0yGET0UHxrZT/rV1d5OPldbJW5WMZMapNpm/lJvdfZrFSPFUYF7CkmVCnfAvE411nlG531e1K3CTnAVsAOxqZt/ynKNdM1uUWOknC9cpNs39vq/k/Kq5342WuL6fcBy7rJkN9xYsezqJma1NOB6F8nvqb5Ynntlkwn/OjcT7+FHiIwwAM/sukFkg3A2cSjhSfTml8XVCuXBAivf94oXMbA9qU3EeIxqee4kO8zKE0mWXtE2hYupZk+wGZKOD+elpZxM+TJYhfBxd18vrAPwf4bztEmJ+8eOEb6KFWkhjVO53X5Xnwwg5byB8zDxIOInOW08NJayqLibeh4eID66liU7nAcR9nWZm97v7TYVr/In4EBpB+LMpngfC+oRQ6gFc6e6vNCH/O8SUpw8SDrohFH5XF8Jl/k3OIqYIzkk4xL69Iu0vp/0k+nZK1Zq53z2cTpvZQYTCGOKj+lgi/2Yj3tVDied1opm95u7v17fuPiW1J5umsEU+Sq3uAVjPzOZ297dy1x8KbJz+ji8mYGZ7EUpdCMVZVoafJ6xH9iScKH+OWI3tS2UZkNiTyIv/EArZe4BhSc5muYuor4YCnzezi9z9whbil7EhoVx8Bvg5Mdo9J6FIPTDJeIaZPe7uZc7ke51HadTySqITPY0oUxcT7+YchLJ0a3pObT2XeGYHEx8YbxO+24qUlbMViXrrNeCHxDs3jSjvb5SE7wRDiJXMFiHanr8T9f5mhHJxOHCmma1HrLhowHeIOsmJD7LDiCl1Y8xsVXefVrxIEwwl7vFyQjn+QPq/OPF+fpOo+35sZg+4+18a3NMlRH36a8Kx/svp/xHElOYtiEGjY4uRzWxTQtk8G/H8fkXUbe8RdekR6XzbU5b7gKOIfBoLnE5Yji1IT2uyocT7fzHxbj1MPOuliX7lfkS9Nia1KROLF7Gwnvsn0Y8D+DeRF3cTiqgPAJsQ5SrjFaK9WI2aJfV+zGj9na+HhxJ5fxXxTt5HlItFiD7nAUSf5wAze8jdT6rIm1bZl6ijLyLembHEu5LnbQB3f9zMxhL1wF5m9r2K938Pog6DyLP+4lpq/uI+TqFdNbNPEu/EEKK+PYnw6/gM4adxV6K+3IqwTN3ek3ahhE2IQYtHiP7YRKIcbUmUm2HA78xsnLsXF3s4jlg9EkJpeQHwJNHfWJwot9sR9U6zHEP0sU4hvgXvoeYuIWMqgLu/aWbnEv3oncxswYr+0LZE3kD/PkvRKbqt2RlsG1F5vU0UwH+VnM+W5P1+7ti/07GDCmEXpabZvLgkrbw1yxYdvo/HUrqZxcqC1DTBnymE/WY6/iJpJKYYv43rr01txLDUsoOaFU3DZSWJyi3Lq5VazIPHOpivf83JMYMVTgvpjMmlsxVhvVRvG14S/5Zc/CeAD1Rc6yO5Z/Fj6lghEB1JJzqByxXOLUF0Vh34I3VG7OhpYbBxL/M6G4GYUDg+J7WlWc9skEazFisO/LCX8p6VS2uTDr5zKxfk/EO9Z5jCf4DqkfSFqC1vPsOICaFQei+d/01FOvlRsDKLtI4st0ysYODpmZfeF9FRnp7CHd1mPje0WEnXySwZpgKLF/I9G1V8Elimzn2/lMK8BixQOP/DdO7dYrknlDJOfFRk7/+WhTAb5u5h+8K5pamNDp4ODK1zj4fm0li/cK64zOsV9HJJziRLPs27CIXITsDSLaSTX9r8YWDRkjAb5p7RfcVy1ME8yqwzJhN+0+rJvDgwR+FY08st03OU9jFgqYqwnbZYye7vIyVhPpcL8zyxqt2yJeG+kwtXOiLeRB4sTp3R8HR+GPFB78C9dcLk72kKJW1XypfMOueJkvNGrV59p04aI+hpudMJi5Xrqe4/rFASv2jVUdd6I4VftvieljyD/6W0LqoT5ozc9Y6tSGsohfqXFpZbJvq6Ve/3UGoWOs+W3Rf9tNwyMTCRXWe7inD/SWHuafN9addiZf1cvPML5+ZNZdtTftZrm/fMpfGpkvN5C+jrgblLwnwmF+bIwjmj1ib/rsH9LFxyrCPLLROKnSyd/SrCZVZALwNztvM8tXV3k4+VFvGYW3hr+rtufl6hxbLI2QhS3sle9rs4bzP/v8y/Sn6e8wv1ZDKzOcxs9Yqt4ZxfD+3p5elv0ZIgWw3ofHd/r1FajUgj2WcQWmyIJZPLmDftJ9c5nycfppWVfjqGhX+Lz6S/E6g55+0t1xAfE/W29RrEP9Tdq/waHEo8i9uBH3iq3Us4jHgPh1J7JzL2J/L9aeAb9d4Td/8/YsQCahYELWMx7z4bgeix8lIqo9mo4875ebG94B5C6dQbmi3PczUoz1WWhi8C36x4hrj7k17hwd/dX6ZmVfEJM5u3cP55aiNTu1bUL5mT3leJD+y+4vS0X5D6Ppa+THSwnFD8dQwLFkuWDP8kRtgBzvKey1B+ldqo4oHuPqmYlocz7syv1nzM6CdofNrPTs3yJGNU2v+DWpszqk6YafRsoyBGaocTVhP7eh0n3cRHZmZBWFWG3wW+0oE241vEyHLG6kRddBHwpIV/ijFWx+9MHQ5w9xnKoLvfTFjXQCi5iml2Io++QignAQ7zcHhcirs/5znr0V5ysM84ktvX/Nzd/1ty/C/U6sBFiXx4oiTcadT8OLS1lHbKw7orc3n4uTo8/V3F6viXyvFLn9GKD3efDPwm/V3GZvQpMoqaI+bf1EnjMWL0vZOMpLr/cHH9qEAogQ+tCuDuT1S9p6ke/Gn6+0mLxRDex8w+QK3f+U9qz6Msrak+44psTePhe/DFivNTqfkxWpxQtnaLSwnlBNSpa5OVdLb4wxllYfqQl3K/i9a7exNl+w3CP2Rpn8Pdz6bmb7GqPZma0nmreMJjMYzMB9emhdPzUmuTZ7BALKTzUtX53uDutwOZL7h6z3Ixaj5bzkl9WTHIkGKlPbLCOTs9K92PEfPg3yOsVDIyB7abFJw05RUrZQU+/0FTpVz4ENUN59IVcfNk0ym2NbOFAMxsVWofsJWrAbXAaeQaAne/vE64TGnVTMcyXwHVdaLVV5jZmtQ+8N4C9qz6uO1H3qL+suCY2WyEyTXE6EnVB3leqVjsbGQm65dVfbQnsg+63nRYsg/OaYTFQpHsXR1ObTpKbzjP6ztPbpZmy/PKVJfnuo4lCcu3lqbBmdlwMxthZqtlyhtqZW42wtS6SKbMWpiw7CmmuSC16YMX9nEH4R/ESDHUVlDLyzIbNUXgP9z90Q5c85nMGR1hCfMcMdUiM+H9J/ERnidzWvcy1UrXs4lR8XycjFty596fDmRmQwhTaYhO6rhimMSotL/D3V8rnMvK8KUNPpKc5srwP9z9mYrzTZHe562IqX/ZNJE8yxD1wVgzG29mjdq7p5lxalme/MdJMf87kUdZffsG/Wfq/QYxhaW/KXXunerSbOrze4SlZ1m4VwmLS4h+Tq8xs2FmtqyZrZqr7/LvVCMHsn+qOJef6l2UN6+kq1LunkeaDjJAaLn+tnDK/MFCm5LVW3MSbVyebeB91wS/7s++UxqYXMbCsXMm6zDivYTG70OfkZTSZ6a/O6Rp+UWyj/T36FwfvVny/Zh5C+eyuvLakramSKZcrmpPbvSKlUaplb1iuZtMbbrjHkWlXj/z+7T/qJU7qt6T2vR2TQMapEix0h5FPysZ2YjKBO+52kemWFmYmjfwfPjXqWky8+TnPveHFcYVhAY6v7Ry9gH7gLv/uzRWC1isJJKNTNzCjB8febI8bKYinDP3u187Jclr95XER/x0Qqt+fwcv8UF3t4ptfEXc+xqMGH+YmPsM8LMSr+Y9Nmra9CWyBJLVVvZef72JNPYrptEK6SMys4a4tmAVkHET4f8AmlgdqAnKRl1bpT/Kc1NymtniZvYLM3uIkOtR4kMnU97klXFlipyLqXWq9ig5vyu1clv1IdJrUic8U2p+wsyWKwTZmpiGA33bWXmPqNO+DmxWouDKVpuZUGHpkI18Zx+daxTOvUNt1aG80iTzrzIlnR+fjq9raYWIKv8qyaor+9jZr4ky/PUUtqoMd6LMAPGMPZa6/zhh7bET8BPCmi+fz5sC/zSzhSuSm9BAgXw/0SZDzl9OB/Mo8zFzWxNK6E5xb9U710dMp3olp1fT/kmvXh0tC1f8cGsaM5vPzH5gZncT9dbjhBViVt/lrUeqFNfTCb9V9Xg597sob1aW32ZGf3rvk/Ki7vk2uKJB/6HRKljNtilLW6wo+AhRfh6hZ5uSbweKeZz3u1S5pHonSMq1Q8zsdqL+eIKwMMsPYGQfuFXvQ3+QtW2zU7Amt7CS3z39vSxZk/Yn+Xe8aBWWrfpVtvJRsa78YQq7SOrfldHIZ2FW9ooWttOp+YbcFnjEYpWqT1nF6md9xDnEYCeUDAJR88c10d0brnQoBiZSrLTHzdS02XnFSva7R8OQzA7vz4exWKo501je5OVOqfJmaWWa6iz9+4uNJTUtd9OkD/BshGnPNNKbfTT1WhNuZodSm/ZzJ7Gq0gxmfTmyD9FmPkLzYZqZOtQRzGxJYlpENkr6de+9c8VO0shh6GINztdj7tzvRYhpFr1JoxW2pGYZUPrRnj6csnOjSj62W6UZx6uNaLY831FSnv/c5DUaymmxLOF9hHn3CjR+djNYgKVym1ld7FicLkSt3niC8mmOneaPhKmwMaOz0GxE7xU6Nz1vFPGhlDlOXA6Y1903dPfTivV5UmpkedRM5zczdS9zjpxZo6yTy/dRaX9LUr7cSdxvfsrQeoTyF2Z0XFv3fWxAVRnuRJmZAY9lxy9x9++7+zaE7PtSay8+SG06VRnN5H8WJq+g6XUepY+GBdLfXlvztECfPIsGvN3A8iCzAKzqA+TD1fvgqiRN7bmHmN64Go37vVUWr83eE8wob1aWX2rC+rFssKBbNNOmjCLalAOJ8teIYh5nZWuql0zR6yTJou0/wC8JhU6jRTz63QI6j7s/SM2ivdi27UTtvervaUDQU+n0vlIxDba1M3hk1M/v3tQTh1KzbF6amF56MfC8md1lsbT1kiXxOkqy3Mmmqe+et54xs/WpDVB241mKDiHFShukkchs9YkNLNZHH0J4dYdyjXvRz8om1PK/3odHfqSglZUUekOm2d2AWDlhacJMtleKFTP7BvCL9PdBYCtvvEpI5hdkuJktUBkyzMGB9x1D9jlJ2z2W2oorB7r76RVRukGjVRTyjdDR1D4WG2071Enj9BbSaPedzlugnFMxCpKtXGT0fgWidlajKNIf5blSztThuYDwR/IuUSZHEqPqc+UUOXnLunqKl0xxNYycb5OkxMqmpZzXH2bdac595sdltFlMuUyWCzum4+d20ELgAXe/O233evgYaNZcvrf5MT7th1LL51FpPw7eH6XLOuObFcKU+VfJl+Hf0XwZrvLv1Iky0xB3f9vdT6XnR8dnq6K0ealO5JHRnhK6t/TLsxhopHrgPMJibTqxiscW1Bx4Z/Xd/PlofSzWQJgi3AqN2pT5CMX/vIQ1zjGEMndxeubx+vlofSRrM5wOZH50ziGmrC5HtGOz5eTNrKi6KWtG1q9c3czyK4FlgwaNpjf2Ffl+zAO53/m68q80X1euQWMFSsukNmJ3YjD7p8Tg+LvEs12d6Cs+bGZV7UanyJ7lwtRWVILas5xCvJdikKLlltvnBkL5MIzoOL1LNCxO+fKjNxLOC7PpP438q0DPUcVtqSkm+gx3v9XM7idMno9Ph6/3csdyTWFmXyRW7YFYmWDzOtM3iuSXWF4V+FdF2MxE+4lWfUy0g8XSgNdS+wD9vrv/uq+v2wfknbi95+7tmCDn07A202iKlO87NQw4I3sRUwe6yfjc722pOcnsT7ampoTc293rTdNpZhnp66gtm7gHNaXs7tQ6o306DajA74l53csRH09jiSlj2ahQ1+Ysu/tUM3uDaCMWbyJKNn3k5ZJz/yY6/cOAzczsWmoKlvG5cOOJsjIq/c/2Zf5Vejhz7Msy3Fe4+1/N7GXi3V3azIbVmWLSTP5nlnx5K7Ne51F6D14lrFb6fIS0RfJWFFWDbsMrzg001iFWIIRYLeTndcI1U9/1lmwgaREzm62B1Uoz7+hAYUdq5WV3d6/nDLcqj7OyNdTMFu0rqxULJ7mZT7AT3P3bdcLNRZctVQpcSPQXFiSmkNxqZstQ89szpo7Ve1+zVe73+8p6d59sZu8QU/SHDpT2xN3vJC2ekAaZNiF8d+1FWBb+ycxua+DLpbcy3Jj7xtobuCDJkrlf+JuHbykxSJHFSvsU/axkCpO761hiZJXOUma2Qi78W/R0ePY+qTLKLGM2tfD+3R9k1ilzpf1Z9QI2wsx2JS3/SmjVN/fqFWry3Jj7XfT0nb/G0oSvkGKcPiHNtb+aWoftWHfv9kd7uzxIbYRgk6qA9Ug+IR7uTRotsCu19/L7RKNYtWUKvQ+bWTe9+0OU5ayDsZWFY+j+JlMEOmG5Uo91K85FAtGRy6YObm5mmTIgmwb03150qNoZ1b2amrVaNvqTWTHc4e7/aVOWTpHlxToV88hJPlEyvwd3Fc+X+FnJ/Ku8nTsOtSlD61o4Ey71r5LSfJVQekPfl+G+JD+9pt47tG5m0VRGmjoyX/r7fv53MI+yNn09y60q2AJ9ZfGQ9wG1YEW4ouPRgUze8q5qOmXD+q4DZO/SMGrlewbSO9HI78ksA4axAAAgAElEQVRAIsvjKVQ7SK7K49tzv4urZzZDs2VioLwPLZXhZGmZDVJ8Pr0jo6l9w3V0pbtmSN8xmSPu55lxQDlrbzesau+6RbJiGevuXyamsEEMwrQ6cNdOfZw5sd0iKch2odbmaBrQIEeKlfa5kdoIz0jq+FfJ8FiJ4qn0d2tqFffNXu1c9Jjc7zGZI8I+5k9EJ/0dwnFcWz5DzGwHwqRtCLG04ubu/kh1rB6MpzbKMzr5fCkjbwLeKR8KpaQG7TLCWgngJHevuzTgQMdjZYur0t8tLZYxboesQ7WSmX2yMmTvyKYBvUws5Xl+1UaUn+mFuF0hTYnJlpw04Kw2P6x6Q2alaNTxj5H8gXylyfSyzt4QosO3FrXOa2+sVfJTduasGypHUvRknZKdLJbfzfxYDQQP+9kS1QtTWzGhjC9QUx6OrRMmU5qsTW2q083ec6WauwiLi6HAt6nvXyUjK8OrmVlxNZyuUKUAKQm7ALVpmc9UTPtaip4jrUXyS2EW878TeXRZ2s9LuQPDRmT31enVLfJtc9UUrzJn1QOVvFV2Vd9pn74WhJ5Lhle1RZ9nYFlLNGJobl9aVydfEkX/IHmuJnxkAXyrlXKfaLa9GCjvQ36FpGbJPsbnI6Y6jk7/r3f3h0tj9BHJae6fqH1D/l9JfZvVlYtRc7A7ULku97tVh7btPMuziFkOsxHPMWtzHiVWORSDGClW2iSZUmc+EzamNoJV5dE8s6Y4iFqnqNKxo7tfRE0bvSZwg5mtWBUnNUqNfJJUXfMJd5/b3edy9wXd/Y3GsWaQYQvCSdNQQjmypbe4Uk5SOJ2Q/q4IzKDASKOLh6W/D9KHS0qmzsFfqZnU/4FwgjXY+SmhfJgNuNDMlq0X0IJPl1hb/Ira6hx/tFh+ui5mtpmZbVwVpiTOCtRG3S/xJla58FjuNZtC9rnUIega7n4eNSdq6wDjzGz5qjhJoTh/VZgWeCj3e3SdMMfRc2SvLu4+kZq3/j2ordY0nfJlsJvleWoKscr8KXBGijcXNcu7gTJn+ffUloX/PzNbqhgg1e2ZMv116jshH5/2Q6itrDY+HyAp8rL2JaunyvyrZBxPza/AmRbLjtbFzDbvByuwj5nZ5Y3qCjObnfB7kpXvRgr2k6xkRYjkQPCb6e+DzKhY6UQe/YEYaAD4uZnVHaE3s8VsxuVBM6uc2avq6jaYSM1q5dspT4vybE97yqBu0bC+s3Cqv3k/yDKO2iIG+1k4ES/Ksiw15ftgIcvjodTq//dJ/dFTqFguO1kxZ4r4jwP1pmxhZkNz1pEZeUu1qvaimfdhL/peeZjJ23Tb5u53EdNAAY6llp/9OmhgZisTbUjmM+cWoGwa/KnUplKe1KitMLP103dDRzGzDzShBM8r2h+tG6qc7FkulwalGuKxoEnWRu1HzSL/j/3hk070LfKx0jtuIEYM58sda6RY+Rw9G5h6/lXy7Et8WH2GMPu+18wuIXx8PEgoLuYg5uWvQ2izV0lxpxGa0X4jdU4vITq504BDgGkNOqHPe/lScccTebYKcEz6uD6HWPlnQ2L1h3nSdQ6o97Gd4hXNtzOv5fOY2ejCuRtLRgHOJRydQTi/OpEYuax7UwNlbmkV7j7RzA4jvOSvDNxlZqcTI2zPEs9xWcJKZ+f0e0tyPnDc/Wkz+zIxNWRxYg7wmcQy1JOIumYpYiT0M+k6X6XcH1E99sr9bsWK6kLi2S9AWApUTYHpD75CjFZvT3RO7jezi4ny/BBhJTYHkV/rENOfMmXqe9RWJGuHy4n6YkHgl0mpcxnRAVqRWCJ2U+K5NKv4+hOhDFiXmsXAeHd/qn6Uatz9bTP7D3H/X7dYJvUuavf+Ytk8ZHd/3MzGElaBWef7oiYcZfc57v6kmX2XqNOWBSaa2bGE4m82QmH7HWpK8QMq5lrn/axkSrdxJeHGEY6FszBl/lUy+SaZ2d5E/bokcFuuDD9JlOGlqZXhlYhR6Jsb3nz7GLHE+3Zm9iDxrt6a5HmT8NuwHjHilzmkfJ5wxF2P/xCWTBPN7OfAbUQdtzVwMKGUmw7sU+zkdiKPkv+BvQhny8OB68zsT8QqFZOIsr8i0dn/DPGu5P275H2NnWxmvyBWkclkfaSB/45SUpn7A2HdtE6S63hiZa/Fiffoy+leWlKKd5F/EUs+Lw8cZGaLEW3Us8AIwnJkR1qr79rC3d3M9iXa1Tmp5e81RL22EXAE0Z+8m8EzHegiok6bBzglDbpcSbQzqxD91w1pnMcHEpbfHwK+Y2abEo4+7yb6sEun+LsRfZWTs4ju/rKZPUCUt2+a2cNE/yTrD77g7q+5+4Nm9m+i3d3DYlr3GUTZXZqwFtqNUDKuSs1ysNP8i1gAYBUzO4oo+9lqllMqpsufnmTP2rbXaNOivIIlC3314UQ9uwahgNySmg+1CcBOZdaB7v6ame1GvAvzE4PC5wKXElMqZ6P2zbITUScfQU/Lrk6wBHBteicuItqPScS7sSTRD8ssdJ8nBk9b4V/E4Ma8hALpDGqW9tPSbIUyTie+bTJ/StOBMS1eWwxE3F1bmxvR6fHc9kiD8B8phJ9CeE1v5lpGNDwvFtKot00jOqGr1knvsRTuljbvvW584KgmZcxvR1Vc64NEI1kv7tvA6Abyjm5RnhnSa+OevBfv1phcOiPaiH9Lint1C3G+RnysNLqvqcCGddLYkfhIbyZ/PteCbEaMJDhJkdhC3GxFCAeuKJx7Nh3/bUm8bXKybtDus6wj02zER/TLTebVNKJTsFJJWivnwn2+iWvvQFhO1LvWtYTCuKk0CWex0wtpfKkJOermfTpfrF/z2+EV6e5cCLtFh57Z+bk0l+hFOj9Mz7Pevb1HKFUapXNdLs6bZWWC6Azn0z6uiXQ/3eR7OR3YuRB3ntz5QzqQ56sRisZm69w7gJXrpJW1nScTHeF6z+DdRu9vb/Iol8b2RF3WKI1FSuJe0Ux44sPHgcubzO95CUVTvbT/TTi/rvuMCYs3ByY3uNaFKdzdDcK1dA8l8TcgrL/q3dNtRB3WiXtaPZfOLnXC7E2U8TJZ3iGsJZrKmwo58uWw5XwjlORZ/O2bCP+5intyYqBto0ZpEoMJt1Skk237l8Tds5nwRHv5fEXY+9P7MJk6dSbhEyMLv3rJ+ZPTuRfr3Oei1Nq/4jahwXN9Ixe2tO1s43nvUkeWetvzwA+A2ZtIexRh1dFMuvuWxK/7HJrJ88K7XLU9A3ysJN1DcmHmKTk/F7Gce1mapc8/xTNC6ZuFbbqvrm1gb5oK1DtuIApERiPHqXcRGuaM27zJ5T89+D9ilGVvwsz+QaJzN5XofD4M/I1Ys32Eu+/g7veWpzh48ND4rkOMJv6b6IxOISqlU4G13X1M1wScSXD30wgl1g+Id/kF4t16i3i3LiGmFCzr7qWj1O5+aUrjYOLD71miw/U28Djhz+Vw4MPuXuU8rsimxLsPcJn39CXR6L6epObUc+sSM+J+x92nu/sviA7cVwlHeg8R73ZWnh8iOtgHE3n+aXd/oE6SrVz7MmJE/Tzi+UwlOkrjiLplG2rTHZpJ73F61n1TaH3UpyzdvxGj9pdTe4+a4VJq9exj9Jw/3XXc/UeE5eHviXL1VtoeAH4DrObuJzWR1Ljc76J/lYy76WntML4J+S4iyvAhxHzvsjJ8GFGGe/2cG8hyDzFHf2tiesB1hAXF24Ri5DXiHs8mlB3reBNTTlP+bkqUr6cIZcrThCXKuu7+xwbxe51H7n55SuO7xGj+S+me3iCmGf8G+LiH2XiRnYEjCaefmdKg13hM+92UsAS9K93PG+k6BxGWf123/moFd7+FUBT/gRipfo/I65uJ9mxjeq7+1Nfy/IGe9e+7hMXEOcSAxUCYttgSqS3fkChPzxF5/BxhfbCnu3+KJiyn3f3plM6uRF/2SULZNIUYWPkrYWU0w/QXdz+bGDS4Kl271Ho51Q9rER/ijyS5XiHe8SOIOuTx5u68PTxWPfoYsfz3gzTZ3nosFJCf5tjX04CmEX2RJwhr/F8TSpgPuPvRXu0fEgB3H09YIe1HPJuniTyfQjzfvxN9ztXd/Td9cA8TiXrrx+laDxLtxlSibbyeGORayd1vbTXx9A03klC+3kNtSnyjeE7PKcoDwQ+c6AAWz1YIIYQY/KRpf9lc+h+6e9W0EDGLYGYvEo6DT3H3/bstjxBCtEqaIrsWcJe7V/qyEwMbM7uImAb1ErBUKwOGYuAiixUhhBAzE5mHfc1ZFkIIMVNgZh8llCogC4dBjZktTvgPAzhbSpWZBylWhBBCzBSY2fyEA14IfzpPdFMeIYQQokNkK2C+Rf0V48Tg4CBgdpKvnC7LIjqIVgUSQggxaEkjP/MSqzocRaxgAINv2VIhhBACADMbTqxcMw/hd2bXdOpUr79inBiAmNmchOPvuQifYQelU3/rhO88MXCQYkUIIcRg5gRiVYo8pyWnlUIIIcRg5JPAXwrHHqV6OXkxMFmDWIEszyvE4gRiJkJTgYQQQswMvAPcR6zUsl+XZRFCCCE6gRMr6PyRWCns9S7LI3rHc8QKVxv29QpUov/RqkBCCCGEEEIIIYQQbSKLFSGEEEIIIYQQQog2kWJFCCGEEEIIIYQQok2kWBFCCCGEEEIIIYRoEylWhBBCCCGEEEIIIdpEihUhhBBCCCGEEEKINpFiRQghhBBCCCGEEKJNpFgRQgghhBBCCCGEaBMpVoQQQgghhBBCCCHaRIoVIYQQQgghhBBCiDaRYkUIIYQQQgghhBCiTaRYEUIIIYQQQgghhGgTKVaEEEIIIYQQQggh2kSKFSGEEEIIIYQQQog2kWJFCCGEEEIIIYQQok2kWBFCCCGEEEIIIYRoEylWhBBCCCGEEEIIIdpEihUhhBBCCCGEEEKINpllFStm9piZeWGbYmaPmtlZZrZWt2UcKJjZMDM72MxuMrNXzOw9M3vezO42s7PNbB8zG95tOUVjzGw5M3vLzM4vHB+VysBjXRKtEjM7wcymmdlHui2L6D9UTzeP6umZh0FcTx+V5BvTbVmKWHCHmT1hZsO6LY8Qg4mStni6mb1hZpPM7O9m9jMzW7PbcrZDdk/dlkPMHMyyipUc1wBnpu1aYC5gT+A2M/t8NwUbCJjZksBE4DhgHeBO4ELgJmAo8AXgVOCD3ZKxHXKNxIhuy9LPHA/MDny/24K0yE+Bt4Ffd1sQ0RVUT1egenqmY7DW0wMWd3fge8AywHe6LI4Qg5WsLT4LuAp4GFgXOBz4r5ldamZLdFG+jmFm41P7M6rbsojBw9BuCzAA+Lm7j8/+pJGM3wN7AL8zs2vd/eVuCTcAOBlYGRgHfM7dX8ifNLNlgS8Ck7sgm2gBM9sY2Bk4290f6rY8reDuz5nZ74CDzGx7d7+82zKJfkX1dDWqp2cSBnM9PdBx98vNbCLwHTM71d2f77ZMQgwyerTFAGY2G7AD8Ku0v97MNnL3l7ogXzus0m0BxMyDLFYKuPvbwDeAN4H5gK27K1H3SB8vO6a/+xQ76wDu/oS7/9jdH+tX4UQ7fDvtT++qFO1zRtp/q6tSiK6jerqG6umZjsFeTw90zgDmBr7abUGEmBlw9+nufglhufIwsCJhdTcocPf73f3+bsshZg6kWCnB3d8AHkx/l8ufM7PhZnakmf3XzN5M2x1m9l0zm7sQ9qPJjOzfxWuY2a/SuffMbN7CuW3TuUtL4i2T/E08YGZvm9nraU79aDOzkvDvm7KZ2Ugzu8LMXkzzI3dqkBULUrNqamlkx8zGpOuONrO1zOzidN23zex2M/tSg/hbJ5PC58zsXTN7xszOM7M1KuIsbGZHm9l/Ur68aWYPJVk2SmFGp7mU2XN9tDBvdEQ+XIq7sJmdaOHX4V0zu7gYpo482Xz48fWOm9lcZvZjM3s45c0jZvY9MxuSwi5jZn8ws6csfEvcZWZfaOYZFK65FLAT8Cjwz1bjV6S7moWvi0lm9k56xlea2Scr4sxjMR/3kRRnkpmdbGYL5d+bYjx3vwe4HdjczFbs1D2IwYnq6fdRPa16ulG6/VZPtynfdmZ2VZLr3XStM82s7kiymX3QzP5k4UfobTO7x8wOMbMhVj2F7DzgPeDrFiPtQogO4O6vUFMMf8EKU4JSHf2TVD9OTnX/RDM70MxmL6ZXaJ9WMLNzU1vzjpndb2aHlZXhVF8fntKenMI/Y2Y3p+vPVQjfw8dKVvcDm6ZD4wrtzygz2zP9vrpefpjZGinMU2amGSKzCHrQ9Zkv7d/JDpjZIsA/gDWAV4i5hgCbAccAu5rZJ3Im6XcALwPrmNkC7v5qLv3N034oMAq4rOTc3/MCmdlmwEXA/IRW+GpgHmAD4I/AJ4C96tzPZ4F9gHuBscAiROeiihcJvxbDCCuBHzUIX8b6xNz+p9J1FyMqqzPMbG13/2YxgpmdAHwTmArcBjwJrAB8HtjJzHZ29ysLcdYGrgCWJPJ8PDCF6JjvloL9i8i3M4FdgOHAX+lpHl80lV8kyTA/0dGdAHTKvHEOIk9WS/I+BIwEfgwsbWbHET4S3krXXhrYBDjbzNzdz2nhWtsR79p1aa55rzGzHYELgDmBe5KMHyCsBz5pZj9x9+8X4sxL3OtHgdeJObrTiGe7dUqnir8TPiR2JPxJiFkb1dOqp0H1dF26VE+3It/PCP8M04EbiXdwTaKM7Gpmu7j7FYU4qwPXAwsBTxDlfQHgJ8DHqq7n7i9bTAdaP93fhE7dixCCK4m6fSGizT0PQslAtIVLEW3FeGJwf31iCtF2Zratu79bkuZawAlEWzeOaJ8+DvycqMsOyAImRcsVRDv7GlFPvAYsDqwEHElMnX224h6eJdqfbVK8awrhnwVuJvqgW5nZ8u7+v5J09kv709x9asX1xMyEu8+SG/AY4MCoknNrEZ0IBzbLHb8gHbsBWCB3fEGiY+XAeYW0LkzHd8odW5ToRNyZzp1QiHNHOr5a7ljWEZ1KzJW33LllgP+kOKMLaY1Pxx34Whv5dGIu/j3AL4FdgeUbxBuTi3cCMCR3bn2is+bAtoV4+6TjdwMrF87tRHxkvAIsmDs+DzApxTsVGFaItyiwSZ3nP6KO/KNz8l8DzFsRZkydNEal8+PrHHeikzt/7txHgHfT+3cv4aw1n3f7pXgPt/gcz0nxvtxA1seaTG8JorFy4KCStN5M57YunPt1On47sEju+HxEg5nly+g61/1UOn9lq++ytsG3oXq62XxSPa16uix8v9bTwFFVeV0i37Yp/GRgZOHcoencq8BiueNGrRz9AZg9d24l4OmcfPXem/9L57/TalnTpm1W3Khoi0vCjk1hf5L+DwMeSccOB4bmwi6UC39UIZ0xubJ8FDBb7tzIVP9OA5YpHM/qruGF9AzYGJi7cNxJ/q0Lx8dX3TOhXHfguJJz8wFvpHZiyW4/P239t3VdgK7deEklQXS8d8xVAP/JCjIxojY9FeKVStJbtU4hzzqgJ+eOfS4d+1rqBNyTO7dIus4zhfSPTXGOrXM/62aVSeF4VjFc22Y+zUF02t/LVXDZNolYrWXBknhZhfgkMGfJ+R+l82Nzx4ZQ6xStWkeek9P5A3LHvp2O/Yvch0yTz39EnfOj0/l3mwgzps75UVR32KcBq5TEuzidfxSYo3BuKDES68CyLTzHe1KcDRvI+liT6X0/hb+xzvmflzzfuYkOtFP4gMqVoelUK1Y+nM4/04yc2gb3hurpZvNJ9bTq6bLw/VpP07pi5boU/pg6529O54/MHcs+nF6mXJH2jdy7X++d+ColClZt2rSVb7SmWDkvhT01/c/K5J/rhF8q1eEv0HMwYkyKd2tZm0FYxziwV+7YZ9OxX7dwb057ipWliTb3JWCuwrkDUtwLuv3stPXvpvmlublzREN9CbEk5URi9HJ6CvdxQtt5i7s/UEzE3e8F/k2Yto3Mnbou7TfPHcubkF8HrGqxXCaE+Zrl4mVsm/Z/qXMftxOdobWK8wcTf6sTrxJ3f9fDDHw5oqL4M2GmDWGCdwRwR525zAAXuvs7JcfPTvtNcnMP1yJGfO9J+VnG9Wm/Ye7YNml/hnvUaB1kovedw8fH3f2+kuNZ/o7zglmkhznho+nvUi1ca7G075R5/KZpf2ad85mj2U0yPwTEFJ7hwCR3v7EYIT3z/za4bjZ9Y1GzGX1ViJkW1dMVqJ5WPV2HbtXTDUnv08bp75g6wf6Y9qNyx7J7utzDz1KRc5u4fNaOLN5EWCFEa2Tfllm7XNkuuvvTxBTLRYjBsyJX1mkzMoez+Tp2IqEM39vM9jWzPivj7v4U0WYvREyTzPONtD+lr64vBiZSrNTWZD8TOI0YcdkcWNfdH8+FWzrtH6U+jxTC4rFc4iRgZTPLjm8OPOruj1Cbn79F7hzM2GH/UNrfVnCilH1sTCdMrWcDFi6R7fGSY03j7k+7+8nu/nl3/zCwLDGa+U76Xa/yqJdfTySZ58rJm93jamX3mO7zghRm0Vxay6V9X3j17lW+NeDJOscnN3m+7MOsHvOn/estxKmiUXl4jBmfbxanKk8b5Xcm/xBg3qqAYqZC9XQTqJ7uE1RPz0gn8nthwu/L9Ir0ZiirNJDP3V8jpj9VkeXvAo3FFEK0yCJpnykwszbjLxVtxqopzKLMyBN1rpOV4/frWA9fJwcSVpynAM+a2f/M7Gwz2yWnQO4UJ6b9vtmB5GdtFWLw4frSWGKmRc5rS9Zkb0A7I21/B75ErGZyA1HJZEspZh3zLYjRwXod9qwy+DPh7K+KspHHt1sRuBHuPgk4ysxeIxxPbWVmwzyWQW2X7B6fouAQsoR857zTo595enM/jRSX03t5vhVeJRqs+ah22tUq7eR9VZxG95w5K51KzF8Vswaqp9tA9XRTqJ5uPU4n77nRtdqJ02w78kob1xVC1CFZEq+d/t6V9lmbcQXhgLaKMmu9luobdz/JzP5C+PvaJG1fSNsdZrapu3dEee3uN5nZf4D1zGwdd7+dmtPa33TiGmJwIcVK8zyV9h+qCJOde6pw/Dqiw74FkC0p9ncAd3/SzB4gOvPLAcsDD7l7UUM7iVhx4ccey84OFK5N+6GE74NiB3dEnXjLEh3aKdQq0klp/4y7j25BhicI7fBKxKoC/UVm/j1PnfPL1TneDZ4nOuxlo+Tt8BSwMvHOFz8uIZ579nyzUYun074qX0Y0uG4m/4t9MJ1ADH5UT5ejelr1dH/W083wEqFcnDOl91BJmLKyWimfmc1HvONVZPnb0vLkQoiGbEeUv/cIHyUQbcZKhM+VK+rE6yju/izw27RhZh8hBkXWIhzofreDlzuJmFa5n5l9j1hk4Q1qU2nFLISmAjXPP4kRkg3MbMXiSTNbhVhFYTqxGkWe/Pz9zVM6/yicXxrYvxA+z1Vp/9l2hG+HJn1YLJv271Cuid7FzOYoOb5H2t/ktWXIbiU6W2ub2QotiJotp/rlFvxuZJ3t3igXs87eynXOb1vneDeYmParVoZqnsy8sd6ysV9K+xtzz3cCsSTpsma2YTGCma1MrLZRRSb/xMpQYlZF9XQ5qqdVT5fRV/V0Q9L1bmog3+i0H587lpXb7c2sTFm2W8mxImpHhOgwZrYgseIWwFnuniku+71dLOLu/yVWvoPm669m25/ziDbx86RVj4j7l1X1LIgUK02S5vH/lciz35lZNhcaM1sA+F06d0Eyv87HfZbw9r8U8GngTnd/IRckM6eu6rD/kphP+F0z2y/nSPB9zGw1M/tMO/dXh/nN7HYz293M5i653hrEsowAFxUd+CU+APw8rS2fxVsPOCj9zSo63P09YvmyIcDFZvaxkmvOYWY7ps5dxunEKNZGwElFp5BmtqiZbVJIKutsr1Iic7PcRmilVzOzHp05M9sX2KUXaXeacWk/Q0e5TX5P3PsmZvbN/AkzG0k40AQ4Pjvu7m9Rc5Z4opktnIszLzEftlGdlMk/rjKUmCVRPa16ugTV0/1fTzfLr9L+22a2cf6EmR1E5MNr1KbkQSiL7iIcRv4qX8bM7MPAD5q4rtoRITqEmc1mZjsSde0KxBTQQ3NBTiOsVr5oZkfVaac+aGZf6IAsnzCzbYttb/KtkinRm/UR1VT74+5TiLp2GLU6VdOAZlE0Fag1vkGMeo0CHjGz8en4ZoTp23+pza0rch2wGuFkqdghH0fNgdx0eo6SAjFX3sx2Ai4klrI80szuIUxZFwDWAJYh5va3tbJEHT4KnANMSfMIJxHvzQfpOY/y23Xi/5Zw6rSDmU0gzJw3TWn8xt0vywd29xOSqf2BwL/N7E7gf4TmeOl0zeHAJ0nz9939DTP7FDF/cz/g82Z2E2HevFyKcx49zc8vIp7jOWZ2LTG3HeAwd29qRQZ3f8vMjiY+ps4xs/2IefFrpPz5BfCdZtLqB64kTDM3NzOrmEazpJndUpHORHff192fNbM9ifftBDP7CnA38VH6caLj/RN3v7oQ/0hivuu6wP/MbBzhwX1TogN9KbGUbtnHH8Q0DU/hhChD9bTq6bysqqf7v57eroF8p7v76e5+hZkdCxwG3GBm/ySUb2sAqxPvxhfc/bksort7uqfxxLLJW5nZzUT52gy4HFiPsNKaQT4zW4h41yYhixUhWuVwMxudfs9FtBUfpeYI+mLg6+7+vv8id59sZtsRZfOHwAGpzXiaWIRgFUIh82/gT72Ub03CauY1M5sIPEMsIb8+sZrds8CxTaZ1EWE190sz25La1MFflqw6+BtCmTQEGF+xYp6Y2fEBsOZzNzZaWJO9EG8e4HvAnYS57FtER/27wPCKeDum6znwyZLzt6Zztze4/hLAMcAdxEjUlHQv4wkTtOUL4ce3c58prhGV0ZHAWGIO9GSis/IMMW9/H2COkrhj0nVHE5XuZcQc7reJzsxXKFmXPhd/JNHJfoIwX38VuA84H9i9LK+JpSp/RnQc30qyPkgs2bhBIZu3UdkAACAASURBVOxs6Tnel/IwezYj0vnR6f+YJvLpy+l5TCGc4V0JfIz4IHCiks2HLz2eO39UOn9UnfNtPdOUnw5sWnIuk6nRVryX1Yl5pE+l9+Ilwuxz2wo55iUatkfTs30SOJXwJP/3dJ2tSuKtns6N7UQdoG3gb6iebuZeVU+rnh4Q9XQuTxptRxXibQ9cneR6N13rLGDVCvmWJ5SJL6Rneh9RtuZM8k4D5iqJt3+S4chWy5o2bbPqRq0tzrbpRNs2KdUHPwVWb5DG/MARwC2EgjarV24GjgbWLIQfk641uk56M9TBqV44ihj4eCLVDS8S7dkPgEVL0nFCZ1t2jW8Q7cZbuXsfVSfs3en8Lt1+Xtq6t1l6GYToKGY2Bvgi8CV3H9NdaQSAmW1ANGBnu3u9Oe1dI03beISwKljce07DwMx+RYyQ7+Dul3dBRCFmKlRPDzwGez3dbczs44QflrvdfY2S87cTFm0f9JoPCCGEaBsL57h3EFY4y3nNZ5WYxZCPFSFmEdz9FmKKwm5pLnpXMLN18r4c0rGFiRHrhYArSpQqiwNfJ0ZipVQRQsyUDOZ6uh9lm8fCEXXx+CqEPweI0e7i+e0Jy6xfSKkihOggR6f9iVKqzNrIx4oQsxaHEMvh/ZjwYN4NLgGGmtndhBn3ksSc9/kIs9D9S+J8l5jPW89HhBBCzCwM1nq6v1gCuNfMHqI29W0EsA7h4+AfwIn5CGZmRH4+SfjbEUKItkkOez9F+IRaj5gudXI3ZRLdR4oVIWYhPFZNmcEjez9zPNEYrUaMfE4l5vFfDhxfNgrq7t8CvtWfQgohRDcYrPV0P/I8saLQJwj/QvMDbxLOL88HfuuxetX7eMx7XxshhOgMHyV8d71B+Ij6tru/2V2RRLeRjxUhhBBCCCGEEEKINpGPFSGEEEIIIYQQQog2kWJFCCGEEEIIIYQQok2kWBFCCCGEEEIIIYRoEylWhBBCCCGEEEIIIdpEihUhhBBCCCGEEEKINpFiRQghhBBCCCGEEKJNpFgRQgghhBBCCCGEaBMpVoQQQgghhBBCCCHaRIoVIYQQQgghhBBCiDaRYkUIIYQQQgghhBCiTaRYEUIIIYQQQgghhGgTKVaEEEIIIYQQQggh2kSKFSGEEEIIIYQQQog2kWJFCCGEEEIIIYQQok2kWBFCCCGEEEIIIYRoEylWhBBCCCGEEEIIIdpEihUhhBBCCCGEEEKINpFiRQghhBBCCCGEEKJNpFgRQgghhBBCCCGEaBMpVoQQQgghhBBCCCHaRIoVIYQQQgghhBBCiDaRYkUIIYQQQgghhBCiTaRYEUIIIYQQQgghhGiTod246CKLLOIjRozoxqWFEKIj3H777S+6+6LdlmOwoHpfCDHYUb3fOqr7hRCDnWbr/q4oVkaMGMGECRO6cWkhhOgIZvZ4t2UYTKjeF0IMdlTvt47qfiHEYKfZul9TgYQQQgghhBBCCCHaRIoVIYQQQgghhBBCiDaRYkUIIYQQQgghhBCiTaRYEUIIIYQQQgghhGgTKVaEEEIIIYQQQggh2kSKFSGEEEIIIYQQQog2kWJFCCGEEEKImRAzm8vMbjWz/5rZPWb2o8L5E81scu7/aDN7wczuSNtXcue+aGYPpe2L/XkfQggx0BnabQGEEEIIIYQQfcI7wCfcfbKZzQ7caGZXufstZrYusGBJnD+7+/75A2a2EPBDYF3AgdvN7FJ3f6Wvb0AIIQYDslgRQgghhBBiJsSDzCJl9rS5mQ0Bfgl8p8mktgbGuvvLSZkyFtim4wILIcQgZdAoVv73P7jhhm5LIYQQor948UUYP77bUgghxODGzIaY2R3A84Ry5N/A/sCl7v5MSZSdzexOM7vQzJZJx5YGJuXCPJmOlV3va2Y2wcwmvPDCCx28EyGEaI9LLoF33+3bawwaxcqRR8Jee3VbCiGEEP3FscfCZpvB6NGhZBFCCNE67j7N3dcCPgB8zMxGAp8FTioJfhkwwt3XJKxSzmzjeqe5+7ruvu6iiy7aG9GFEKJXuMMPfwg77QSnnNK31xo0ipWPfQwefxyee67bkgghhOgPjj4avvtdOOccWGkl+OMfo4EUQgjROu7+KjAO2AxYAXjYzB4D5jazh1OYl9z9nRTldGCd9PspYJlcch9Ix4QQYkAybRrss0/0J7/8ZTjggL693qBSrADcdlt35RBCCNE/DBsGxxwDd9wBq6wSjeKoUXDffd2WTAghBgdmtqiZLZB+DwO2BG539yXcfYS7jwDecvcVUpglc9F3BLIa9xpgKzNb0MwWBLZKx4QQYsAxZQp89rNw2mlwxBFw+ukwtI+X7Rk0ipW114YhQ+DWW7stiRBCiP5ktdXCx9bvfw933QUf+Qj84AfRaAohhKhkSWCcmd0J3Eb4WLm8Ivw307LM/wW+CYwGcPeXgR+nNG4Djk7HhBBiQPHqq7D11nDRRfDrX8NPfwpmfX/dQbPc8vDhsPrqUqwIIcSsyGyzwVe+AjvuCAcfDD/+MZx3Hpx6KmyxRbelE0KIgYm73wms3SDMPLnfRwBH1Al3BnBGRwUUQogO8swzsM02Yd187rmw2279d+1BY7ECMR3o1ls1x14IIWZVFlsMzj4bxo6N/1tuCV/4Ajz/fHflEkIIIYQQ3ePBB2GjjWI14Suu6F+lCgxCxcorr0RmCSGEmHXZYouYFvT978MFF8DKK8f82enTuy2ZEEIIIYToT267DTbeGCZPhvHjY+Ctvxl0ihXQdCAhhBAw11zh6f2//4U11oCvfhVGjoR77um2ZEIIIYQQoj8YOxY22wzmmQduugnWXbc7cgwqxcqqq8Lcc0uxIoQQosYqq8ToxB//GHNq11orlml+++1uSyaEEEIIIfqK886D7baD5ZeHf/0LVlyxe7IMKsXK0KGwzjpSrAghhOiJGYweDfffD3vsAT/7WTg8v0aLgQohhBBCzFS4w/HHw+67h1+VG26AJZdsHK8vGVSKFYjpQBMnwnvvdVsSIYQQA41FF4UxY+Af/whl/DbbhPOyZ5/ttmRCCCGEEKK3TJ0K++4LhxwCu+wCV18N88/fbakGqWLlnXfCaaEQQghRxmabwZ13wo9+BH/7Wzi3/e1v5dxWCCGEEGKw8vrrsMMO0ac7/HD485/D595AYFAqVkDTgYQQQlQz55zwgx+EguWjH4VvfAM22USKeSGEEEKIwcakSdGPGzsWfv/7mPY92wDSZgwgUZpjueXC1FuKFSGEEM2w0kpw3XVw1lnw0EOw9tpw2GHw5pvdlkwIIYQQQjTi9tth/fXh8cfhqqvgK1/ptkQzMugUK2ZhtSLFihBCiGYxgz33DOe2o0fDL34Rzm2vvLLbkgkhhBBCiHpceimMHAlzzBEr/2y5ZbclKmfQKVYgFCv33gtvvNFtSYQQQgwmFl4YTj8drr8ehg2LJfp23RWefrrbkgkhhBBCiDwnnAA77QSrrQa33BL7gcqgVay4h0mQEEII0SojR8Idd8BPfhIjIausAqecAtOmdVsyIYQQQohZm6lT4YAD4NvfDsXK+PGwxBLdlqqaQalYWW+92Gs6kBBCiHaZYw448ki4++6Yt7v//rDRRqFwEUIIIYQQ/c/kyaFMOfnkWFL5wgth7rm7LVVjBqViZeGFYfnlpVgRQgjRe1ZYAa65Bs45Bx57DNZdNxryyZO7LZkQQgghxKzDs8+GVfHVV8Opp8IvfzmwVv6pYpCIOSNyYCuEEKJTmMHuu4dz2733huOPh1VXhcsu67ZkQgghhBAzPw89FJbDDzwQ/a999um2RK0xqBUrkybBM890WxIhhBAzCwsuCL/7Hdx4I8w3H+y4I+y8Mzz5ZLclE0IIIYSYObnttlCqvPEGjBsHn/xktyVqnUGtWIF4CEIIIUQn2XhjmDgRfvYzuOqqcG574olybiuEEEII0UmuvhpGjYJ5543llLPv/MHGoFWsrL02DBmi6UBCCCH6hjnmgMMPD+e2G28M3/pWOLnVinRCCCGEEL3nzDNhhx1gpZVCqfLhD3dbovYZtIqVYcNgzTWlWBFCCNG3fOhDYbVy/vnw1FMxknLggWGuKoQQQgghWsMdfv5zGD06rFWuv37gL6fciEGrWIHo3N52WzwYIYQQoq8wg899Du67D77+dTjhhHBue/HF3ZZMCCGEEGLwMG1aWAEfcUQsHHDFFTENaLAzqBUr660Hr74KDz/cbUmEEELMCiywAPzmN2GuutBC8OlPw047hTN1IYQQQghRnylTYLfd4KST4OCD4eyzY+r1zMCgVqxkjm00HUgIIUR/ssEGMGEC/OIXMHZsOLf91a9g6tRuSyaEEEIIMfB47bVY7ecvf4HjjotttkGtjejJoL6VVVeF4cOlWBFCCNH/zD47HHoo3HNPzA8++OCwpNRqdUIIIYQQNZ55BkaOhJtugnPOiT7TzMagVqwMGQLrrCPFihBCiO4xYgRcdhlceCE891ysHHTA/7N353Fazf0fx18fFWVJliluIcttK2nTQiValaLFGi32JftStpRE1myhKGRLshSVhBbcaaPdFlmSm0TIksrn98f39DP3qJm56rrmXMv7+XicRzPnnJl5z718r2u+5/v9fC6En3+OO5mIiIhIvD77DBo2hE8/DfVUTjkl7kSpkdETKxC2A73/Pvz5Z9xJREQkV5lBx46huO0FF8CgQWF70KhRKrAuIiIiuWnBgjCpsnIlvPkmNG8ed6LUyYqJldWrYf78uJOIiEiu2377UJDt3XehYkU4/nho2xY+/zzuZCIiIiIl5913w/YfM5g69e/6qNkqKyZWQNuBREQkfdStG2qt3HUXTJ4MVavC7bfDmjVxJxMRERFJrYkToWlT2GmnUFelatW4E6VesSZWzKyimZ1tZk+b2Vwz+8bMlkUfP2VmZ5lZxVSH3ZA99ghPBTWxIiKSPOk87meK0qXh0kth0SJo1gyuugrq1AlPcERERESy0ahR0KYN/Pvf8PbboRZdLih0YsXMqpvZs8CXwENAe2Br4Avgq+jjjsBg4EszG2FmB6c2csGM4cmgJlZERDZfJoz7mWaPPWD0aHjxRVixAg47DM4/P+w3FhEREckWDz8MJ54Y/j6fPBkqVYo7UcnZ6MSKmQ0D3gNqAwOABkB5d/+3u9d393ru/m+gfHTtNuBQ4D0zeyT10f9Wt24oGKgODCIimy6Txv1MdNxx4bXq4oth8OBQ3PbZZ1XcVkRSx8zKmtmMaLXhQjPrW+D6vWa2Kt/nW5nZs2a22Mymm1mVfNeujs5/ZGYtS+63EJFMcOutcPbZ0LIlvPYaVKgQd6KSVdiKlUOA9u6+r7v3cffp7v6P3eHu/md0rbe77wN0AmqlKvCG1K0b3pjOnl2SP1VEJOtkzLifqbbbDgYODKssd9sNTjoJWreGJUviTiYiWWo1cJS7HwLUAFqZWX0AM6sD7FDg/jOAH919X2AgcGt070HASUBVoBXwgJmVKplfQUTSmTv07Am9eoX3NS+9BFtvHXeqkrfRiRV3r+3uLyf6Dd19tLuX6BvsQw8N/2o7kIjIpsukcT/T1a4N06fDPfeE/cdVq8KAASpuKyLJ5cH6FSllosOjSZHbgasKfMmxwOPRx6OApmZm0fkR7r7a3ZcAi4Es7/EhIkVZtw7OOQduuw3OOw+efBK23DLuVPHI+K5AADvuCPvuq4kVERHJHKVKwUUXhe1BRx8NV18NtWrBf/4TdzIRySZmVsrM5gDfARPdfTrQAxjj7t8UuH03Qj0t3H0t8BOwU/7zkaXRuQ39vLPNbJaZzVq+fHlyfxkRSRtr1sCpp4a6KtddB4MGhfc2uarYEytm9rqZnWhmaTkHpQK2IiLJle7jfraoXBmefx7GjAm1wg4/PDz9+fHHuJOJSDZw93XuXgOoDNQ1s8bA8cB9Kfp5Q9y9jrvXycvLS8WPEJGYrV4ditSOGBFqq/TrF5rK5LJEVqzUBJ4GlpnZ3enWBaJuXVi6FJYtizuJiEjWSOtxP9u0bQsLF8Lll8PQoXDAAfD00ypuKyLJ4e4rgUnAkcC+wGIz+xzY2swWR7d9DewOYGalge2BFfnPRypH50Qkx/zxB3ToELod3nMPXFVwQ2GOSmRiZVegM/A+cCEwJ6oWfpaZbZuSdAmoG+3ynDkz3hwiIlkkrcf9bLTttnDHHTBrFlSpAp07h+r6ixcX+aUiIv9gZnlmViH6uBzQHJjt7ru4exV3rwL8FhWrBRgDdI0+7gS86e4enT8p6hq0F/BvQGvFRXLMr7+GB0Hjx4cOhxddFHei9FHsiZWoC8QId28O7A3cBFQCBgPfmNlQMzs8RTmLVKMGlC6t7UAiIsmS7uN+NqtRI9Rauf/+UOS2WjXo3x/+/DPuZCKSYXYFJpnZPGAmocbKK4XcPxTYKVrBchnQC8DdFwIjgUXAq8AF7r4upclFJK388kuoCffmm/Doo6G1svxtk4rXuvsX7n4DsBeh5dokoBsw1cwWmdklJf00s1w5qF5dEysiIqmQjuN+titVCi64IBS3bdcuFIarUQPeeivuZCKSKdx9nrvXdPfq7l7N3W/cwD3b5vv4D3c/3t33dfe67v5Zvmv93X0fd9/f3ceX1O8gIvFbuRJatAgPfZ5+Grp2Lfprcs3mdgWqAbQDGgEGfAr8BdwFfGJmh23m909I3bphK9Bff5XkTxURySlpNe7ngn/9C0aOhLFj4bffoHFjOPNMWLEi7mQiIiKS7VasgKZNYfZseO65ULRW/inhiRUzq2BmF5jZe8As4ExgAtDM3fdz92pAM+A3YFBS0xahbl346Sf45JOS/KkiItktncf9XNK6dShue9VV8NhjobjtE0+ouK2IiIikxrffwpFHhvcfL70E7dvHnSh9JdJuuamZPQUsI7Rn2xq4CtjN3U9y9zfX3xt9PAComuS8hVpfwFbbgURENl8mjPu5ZpttQlvD996DffeFLl2gWTP4+OO4k4mIiEg2+fpraNIEPv00rJpt3TruROktkRUrE4EOwIvAke5+gLvf6e4bW4y8GHhncwMm4oADQkcFTayIiCRF2o/7uap6dXjnHXjwwbA09+CDoW9fWL067mQiIiKS6b78Eo44IkyuvPpq2AokhUtkYuUywlPKzu4+paib3X2Sux+56dESV6oUHHoovPGGlkaLiCRB2o/7uWyLLeDcc+HDD6FDB+jTJ0y4TJwYdzIRERHJVOsnVb7/PrynaNQo7kSZIZGJlfLAvzZ20cyqmlnvzY+0ebp2DR0UXn017iQiIhkvI8b9XLfLLvDMMzBhQije3qIFnHACLF0adzIRERHJJF99FWqq/PgjvP461KsXd6LMkcjEyg1A9UKuV4vuidXJJ0PlymEPuoiIbJaMGPclaNEC5s+HG2+El18O22PvuAPWrIk7mYiIiKS7r78OkyrrV6rUqRN3osySyMSKFXG9LLB2M7IkxZZbwuWXw5Qp8O67cacREcloGTHuy9/KloXrr4dFi8KboyuvhJo1YerUuJOJiIhIulq2LLxv+O47eO21UF5DElPoxIqZlTezPcxsj+jUTus/L3DUADoDX6U8cTGceSbssINWrYiIJCpTx335X3vtFVatjB4Nq1aFvdJduoS2iSIiIiLrffMNHHVU+PfVV7X9Z1MVtWLlUmBJdDhwd77P8x+zgWbAQylLmoBtt4UePUKv7Q8+iDuNiEhGychxXzasXbuweuXaa2HECNh/f7j/fli3Lu5kIiIiErdvvw2TKkuXwvjxcNhhcSfKXKWLuD45+teA3oSWm/MK3OPAKuBdd/9PUtNthgsvDHvLb78dhg2LO42ISMaYHP2bceO+bNjWW8NNN4UVKz16hNfHYcPggQegfv2404mIiEgcvvsuTKp8+WWYVGnYMO5Ema3QiZWoveYUADPbE3jI3aeXRLDNlZcHZ5wBgweHQn6VK8edSEQk/WXyuC+F22+/0Dlo1Ci45BJo0CBsnR0wAHbaKe50IiIiUlKWL4emTWHJEhg3Dho3jjtR5it28Vp3755pb64vvzy0nhw4MO4kIiKZJxPHfSmcGRx/PHz4YXiNfPTRMOHyyCPh9VJERESy2/ffQ7NmsHgxvPIKNGkSd6LssNGJlQLFC9lI8cJ/HCUTu3iqVIGTToIhQ+CHH+JOIyKS3rJh3Jfi2W67sF12zhyoWhXOOgsOPxzefz/uZCIiIpIqP/wAzZvDxx+HIvdHHRV3ouxR2IqVz4HPzGzLfJ9vqIBhwSOtXHVV6IjwwANxJxERSXufkwXjvhRftWowZQoMHw6ffQZ16oQaLCtXxp1MREREkmnlyjCp8sEHoWtgs2ZxJ8ouhdVYuZFQoHBtgc8zSvXq0Lo13HtvWPZcrlzciURE0lZWjPuSGDM47TRo2xauvz48iHjuuVD8/dRTw3URERHJXKtWQZs2MH9+6JzbokXcibLPRidW3L1PYZ9nkp494Ygjwl7y88+PO42ISHrKpnFfElehAtx3H3TvHl4ru3QJtVcGDQorW0RERCTz/PEHHHssTJ8OI0eGRQeSfMUuXpvJGjUKLSVvvx3Wri36fhERkVxVqxb85z+hPtmCBVCzJlxxBfzyS9zJREREJBF//gmdOsGkSfDYY9ChQ9yJslexJ1bMbF8za1XgXD0ze9nM3jGzs5MfLznMoFcv+PzzsLxZRESKlsnjvmyeLbYIBW0/+gi6dYM774QDDwyvoa7NYSIiImlv3bqwpXfsWHjwwfCxpE4iK1ZuBXqu/8TMdgbGAy2BasCDZnZccuMlT9u24U3hrbfqTaGISDFl9Lgvm2/nneHhh2HaNKhYEU44AVq2DN0EREREJD399ReceWZ4IHLnnXDOOXEnyn6JTKzUAV7P9/nJQHmgFpAHTAcuTl605Npii9AhaO5cmDAh7jQiIhkho8d9SZ769WHmzFCDZcYMOPhguO46+O23uJOJiIhIfu5w0UVh60/fvnDZZXEnyg2JTKzkAcvyfd4KeMfdF7j7n8AI4KBkhku2U06BypVhwIC4k4iIZISMH/cleUqVgh49wvagE0+E/v2halUYMybuZCIiIgJhUqVXr1B4/sorQ7c/KRmJTKz8ClQAMLNSQENgar7rvxOeZKatLbeESy+FKVNCVWQRESlUxo/7knyVKsHw4TB5MmyzTeg00LYtLFkSdzIREZHc1r8/3HYbnHdeKIFhFnei3JHIxMpCoIuZ7QScBWwLTMx3fU9geRKzpcRZZ8EOO4T/oYmISKGyYtyX1DjiCHj//dBxb9IkOOgguOkmWL067mQiIiK5Z+DAsEKlSxe4/35NqpS0RCZWbgcOBr4DBgHvA2/lu94CeC950VJju+3gggvgpZfgww/jTiMiktayYtyX1ClTJrRi/vDDsGrl+utD/ZXXXos7mYiISO4YMiTUUunUCYYODfVFpWQV+z9ydx8LHAXcDfQFWriH/jrR08ylwGMpyJh0F10EZcvCzTfHnUREJH1l07gvqVW5Mowc+Xdx+JYt4aST4Pvv480lIiKS7Z55Bs49F1q3hqeegtKl406UmxL6j93dp/K/++vXn18BdEhWqFTLy4MLLwz7z844IyxnFhGRf8qWcV9KRosWMH9+eH3t1y9sERo8GI5TU24REZGkGzcubP1p3BhGjQo1RSUeObtI6IYbYK+94Oyz4Y8/4k4jIiKSHbbaKmwJmj0bdtsN2reH006DH3+MO5mIiEj2eOst6NgRDjkkdOgrVy7uRLktoYkVM2tgZk+Z2Qwz+9TMPitwfJqqoMm29dbhKdrHH4fqySIi8k/ZNO5LyTr44NCB74YbYMQIqFYNxo+PO5WIiEjme/99OOYYqFIlvLaWV4/G2BV7YsXMugBvAx2BssCXwBcFji9TkDFlmjcPT9EGDIAFC+JOIyKSXrJx3JeSVaYM9OkTJlh22CHs/z7zTPj557iTieQGMysbTYzPNbOFZtY3Oj80OjfPzEaZ2bbR+W5mttzM5kTHmfm+V1cz+yQ6usb1O4nkuo8/DrXMKlQIxeLz8uJOJJBYjZVrgY+AZu6+LEV5Stydd4a9aWedBW+/DaVKxZ1IRCRtZOW4LyWvVq2wNahPn1B/ZeJEGDYMmjaNO5lI1lsNHOXuq8ysDPC2mY0HLnX3nwHM7C6gBzAg+ppn3b1H/m9iZjsCNwB1AAdmm9kYd9cmP5ES9NVXYXEAhNfS3XePN4/8LZGtQHsCD2bbm+u8vNDz+9134aGH4k4jIpJWsnLcl3hstRXccgu8807YB96sGZx/PqxaFXcykezlwfr/l5WJDs83qWJAOcJkSWFaAhPd/YdoMmUi0CpFsUVkA5YvD0XiV64MXfj22y/uRJJfIhMrS4GtUhUkTqeeGmb+rr4ali6NO42ISNrI2nFf4lO/ftgbftll4YHGIYfA1H/0nRKRZDGzUmY2B/iOMDkyPTr/KPBf4ADgvnxf0jHfFqH1z8N3A77Kd8/S6NyGft7ZZjbLzGYtX7482b+OSE76+Wc4+mj44gsYOxZq1ow7kRSUyMTKQ0BnM8u6zTJm4c3d2rVwwQXgRc3Zi4jkhqwd9yVe5cqFrbhTpoTX4CZN4NJL4bff4k4mkn3cfZ271wAqA3XNrFp0vjvwL+AD4MTo9peBKu5enbAq5fFN+HlD3L2Ou9fJU/EHkc32++/Qrh3MnRtaKjdsGHci2ZBEJlZmA38AM8zsdDM70swaFzxSlDPl9t4b+vYNrapeeCHuNCIiaSGrx32JX6NG4Y3i+efD3XeHJ3DTpsWdSiQ7uftKYBL5tvC4+zpgBKFIOe6+wt1XR5cfAWpHH38N5K/mUDk6JyIptGYNnHhiWNk5fHgoAi/pKZHitW/k+/gR/rkX06JzGftk89JL4ZlnoEePUFCvQoW4E4mIxCrrx32J3zbbwP33Q/v2cPrp4UncFVeEhx1ly8adTiSzmVkesMbdV5pZOaA5cJuZ7evui6MaK+2AfNA01wAAIABJREFUD6P7d3X3b6Ivb0dYzQIwAbjZzHaIPm8BXF1iv4hIDvrrr/C6+PLL8MADcPLJcSeSwiQysdI9ZSnSROnS8PDDULcu9OwJgwfHnUhEJFZZP+5L+mjaFObPh8svD52DXnklPJ2rXbvorxWRjdoVeDza0rkFMBIYC7xlZuUJE+RzgfOi+y8ys3bAWuAHoBuAu/9gZv2AmdF9N7r7DyX2W4jkGPfw0P/JJ6F/fzjvvKK/RuJV7IkVd094j2Umql0bLrkE7rorFLVt1CjuRCIi8ciVcV/SR/ny4QFHhw5w5plQrx5ccw1cdx1suWXc6UQyj7vPAzZU5vLwjdx/NRtZieLuw4BhyUsnIhtz221w771hcuVqrQ3LCInUWMkZN94IVarA2WfD6tVF3i4iIiJJdPTRsGABdO4M/fqFCZZ58+JOJSIiknqPPw69esEpp8Add4Qi75L+EppYMbPdzWyYmS01sz/N7KjofF50/tDUxCxZ22wTugR9+CHcfHPcaURE4pMr476knx12CG8uX3oJvvkG6tQJy6HXro07mYiISGqMGwdnnAHNm8Ojj8IWWgaRMYr9X5WZ7QXMIlQNX0i+YoXuvhyoA5yZ7IBxadkyPCm75RZYuDDuNCIiJS/Xxn1JT8ceG1avdOgQtgQ1aACLFsWdSkREJLmmT4fjj4dDDoHnn9cW2EyTyBxYf+AvoBrQmVDsKr9xQFZ11R44MOz3Pvlk+OmnuNOIiJS4nBv3JT3tvDOMGAEjR8KSJaEt8y23hDaUIiIime6jj6BNG9h117BqZbvt4k4kiUpkYqUZ8IC7f8U/W24CfEHoaZ818vLCG7kPPoBOnfQGTkRyTs6N+5Lejj8+rCJt1y4Uta1fH+bOjTuViIjIplu2LOyWKFUKJkyASpXiTiSbIpGJlfLAN4Vc35LE2jdnhGbNQoeC118PxWx9Q39aiIhkp5wc9yW9VaoEzz0XjqVLQ+2VG26AP/+MO5mIiEhifvoJWreGFSvCSpV99ok7kWyqRCZWvgKqFnK9PrB48+Kkp27doE8feOwx6Ns35jAiIiUnZ8d9SX+dOoVaKyedFLr51a4NM2fGnUpERKR4Vq+G444Lr2UvvBBexyRzJTKx8gJwuplVy3fOAcysI3A8MDKJ2dJK795hgqVv31ChWUQkB+T0uC/pb6ed4Ikn4JVX4Mcfw9agnj3h99/jTiYiIrJx69bBaafB5Mnh4X3z5nEnks2VaPHapcB04EnCm+teZjaN8MZ6LnBn0hOmCTMYMiT8j/7ss+G11+JOJCKScjk97kvmaNMm1F454wy47TaoUQPeeSfuVCIiIv/kDpdcEra03nknnHJK3IkkGYo9seLuPwMNgEcILTYNaA7sDzwAHOnuf6QiZLooUwZGjYKDDgpLkFUwT0SymcZ9ySTbbx8egEycGOqtNGoEF18Mq1bFnUxERORvAwbA/ffDFVfAZZfFnUaSJZEVK7j7z+5+sbvnAZWAXYCd3P3C6A141itfHsaODf+2aRMK54mIZCuN+5JpmjWD+fOhRw+49144+GB44424U4mIiMCwYaGr3amnwq23xp1GkimhiZX83H25u3/nnnt9cipXhvHj4Zdf4OijQzVnEZFsl8vjvmSWbbcNkypTp4bVps2awTnn6PVa0o8Fp0WHmVkTM7vXzM6IO5uIJNfLL4eSEi1bhgmWLTb5L3FJRxttk2lmjTflG7r71E2PkzkOPhiefz5MrHTsGNpjbbll3KlERDadxn3JNo0ahW27N9wQ9rGPGweDB4fWliJp4jbgPMLDzsOBk4CPge5m9i937xdnOBFJjv/8B044AWrVCqUlypSJO5Ek20YnVoDJRN0fElRq06JknmbN4JFHQregs84KFZ3N4k4lIrLJJqNxX7JMuXKhoO3xx0P37mEbb5cuMHAg7Lhj3OlEOBbYG9gBWAA0cPdZZlYVeB7QxIpIhlu4EI45BvbYI5SU2HbbuBNJKhQ2sdK9wOcGXAjsBzwFLIrOVwVOJsyu35/sgOmua1f44ovwNGznneGOOzS5IiIZS+O+ZK1DD4XZs+Hmm8MxYQIMGhRWnYrE6Gd3/w74zswWufssAHdfaGa/xZxNRDbTV19Bq1ZQtmx43cnLizuRpMpGJ1bc/fH8n5vZxUAesL+7LytwrR8wDdguFSHT3fXXw/LlcNdd8P33YRWLlneJSKbRuC/ZbqutoG9f6NABTj89dPjr0CF0Z9h117jTSY7Kv+Lv6gLXVIFBJIOtWBHqqfzyS6j5VaVK3IkklRIZsHsAgwu+uQZw96XAYMKTzZxjFork3XgjDB8O7dqpvaOIZAWN+5KVDjkEpk8PLS/HjoWDDgrbeVWWWWIwzMy2BXD3cetPmtn+wOzYUonIZvn117D957PPYMwYqF497kSSaolMrOwOFLYk8dfonpxkFlauDBkCr70GRx0VVrGIiGQwjfuStUqXhp49Yd68UJC+e/fwZHHJkriTSS5x9/vc/R+P49z9I3dXZyCRDLRmDZx4IsyYAc88A403qTWAZJpEJlaWAKeaWdmCF6JzXYDPk5QrY511FrzwAsyfD4cfDp9/HnciEZFNpnFfst5++8HkyfDAAzBtGlSrFlahrlsXdzIREck07uHvwbFjw+tK+/ZxJ5KSksjEyl1ATWCmmZ1rZkdGx3nALKB6dE/OO/ZYeP31UG+lQYPQ6lFEJANp3JecsMUWcN55oXPDEUfAxReHVs2LFhX9tSLJZmaHmdlTZjbDzD41s88KHJ/GnVFENuzqq+Hxx0M9r3POiTuNlKTCugL9D3d/2My2AW4CHuDvlpwG/A5c6e4PJz9iZjr8cHjrrVAFunFjGD0amjSJO5WISPFp3Jdcs74V5tNPh8mVmjXDNt+ePVWUXkqGmZ0FPAT8CXwEfBlvIhEproED4dZbw0T99dfHnUZKWrEnVgDc/W4zexRoAewdnf4MmOjuK5MdLtNVrQr/+U+YXGnZEp56KnQgEBHJFBr3JdeYQefO0Lx5mFy5/np47jkYOhTq1Ik7neSAa4A5QEt3/z7uMCJSPE8/DZddBh07wn33hdcSyS0Jt3Fz95/c/Tl3vzU6ntOb643bffewcuXQQ+GEE2DQoLgTiYgkRuO+5KKKFUPRwdGjw9beevXgqqvgt8LKOYtsvkrAUE2qiGSO116Dbt3CVtInn4RSpYr8EslCCU+sSOJ23BEmToS2baFHD7jySli7Nu5UIiIiUpR27UKtlTPOgNtvD62ap02LO5VksQ+AHeIOISLFM2sWdOgABx0UJuLL/qPcv+QKTayUkHLl4Pnn4YIL4I474OijYcWKuFOJiIhIUbbfHoYMgTfeCG00GzaEa6+FP/+MO5lkof7A+Wb2r7iDiEjhPvkEWreGvDwYPz68VkjuSqjGimye0qXh/vuhVq1Q1KhOHXjxRahRI+5kIiIiUpSjjoJ588I++ptvDoVun3wytGgWSQZ3f8HMtgYWmdloQkv7gs2/3d37lXg4Efl///1vqKHpDhMmwK67xp1I4qYVKzE4/fRQd2XtWjjssFDUVkRERNJf+fLwyCNhyfc330Dt2mEl6rqCf/qKbAIz2w+4ESgPnAZcD/TZwCEiMfnpp9Cc5LvvYNw42G+/uBNJOtDESkzq1oXZs0NR21NPhUsvVd0VERGRTNGuHSxYAG3ahNppRx4JS5bEnUqywANAReBioBaw1waOvTf61SKSUqtXw3HHwcKFoczDoYfGnUjShSZWYlSxIrz+emjnePfdobXj8uVxpxIREZHiyMsLb6wffxzmzoXq1cNqFve4k0kGawDc4e73ufscd/9iQ0fcIUVy0bp14YH45Mnw2GNhK5DIeppYiVmZMmFSZfhwePfdsKR49uy4U4mIiEhxmEGXLjB/fnhyedZZYTXLf/8bdzLJUD8BSXvMZmZlzWyGmc01s4Vm1jc6PzQ6N8/MRpnZttH5rczsWTNbbGbTzaxKvu91dXT+IzPTn5SSU9zDw/BRo+DOO6Fz57gTSbrZ6MSKmf1lZusSPLSZZROddhq88054g3b44eHpl4hISdK4L7Lp9tgjrEK9++7wb7VqYTWLSIJGAh2S+P1WA0e5+yFADaCVmdUHLnX3Q9y9OvAl0CO6/wzgR3ffFxgI3ApgZgcBJwFVgVbAA2ZWKok5RdLazTfDoEFwxRWhgLlIQYV1BRoOaDFrCapVK/RCP+kk6NYtrGAZOFD90EWkxGjcF9kMW2wRnmi2aBEemHTqBN27h46AW28ddzrJEIOBx83sJeBeYAn/7AqEu39ZnG/m7g6sij4tEx3u7j8DmJkB5fh77D+Wv4vjjgLuj+45Fhjh7quBJWa2GKgLTEv0FxTJNI88AtddF8b1W2+NO42kq41OrLh7txLMIZG8vNCy69pr4bbbYNo0GDlS1aZFJPU07oskx4EHhtfvvn3DU84ZM+C558J5kSIsJExy1AHaFnJfsVeLRCtLZgP7AoPcfXp0/lGgNbAIuDy6fTfgKwB3X2tmPwE7Reffzfdtl0bnNvTzzgbOBthjjz2KG1MkLb38MpxzTugCNHRomEAX2ZDCVqxITEqXDrOhRxwR9m3XqgWDB2svn4iISKYoUwZuugkaNw6v34ceCg89FAofihTiRpK8ctDd1wE1zKwC8KKZVXP3Be7ePZp0uQ84EXg0ST9vCDAEoE6dOloFKRlrxgw48cRQA/O558K4LrIxmlhJY61bw5w5cPLJf1egvuceLScWERHJFC1ahNfyU04Jy8inTIF774Vy5eJOJunI3fuk8HuvNLNJhBopC6Jz68xsBHAVYWLla2B3YKmZlQa2B1bkO79e5eicSFb67DM45hjYZZewamXbbeNOJOkuocVMZna4mb1iZsvNbK2KGKZe5cowaRJcc01YflavHnzwQdypRCRXaNwX2Xy77QZvvBFeyx95JLyWf/RR3KkkXZjZU2bW0cy2ScH3zotWqmBm5YDmwEdmtm90zoB2wIfRl4wBukYfdwLejOq0jAFOiroG7QX8G5iR7Lwi6WDFCjj66NBeefx4qFQp7kSSCYo9sWJmjYFJQD1gevS1k4CZgBFmvp9IQcacV7o09O8Pr74K334Ldeqoa5CIpJ7GfZHkWf9aPn48LFsWlpY//XTcqSRNHAI8Byw3s5fN7HQz2zlJ33tXYJKZzSOM3ROBsYQCufOB+dE9N0b3DwV2iorTXgb0AnD3hYSORYuAV4ELoi1GIlnl99+hXTv44gsYMwb23z/uRJIpEtkKdC3wDaGYlgPfATe7+5tm1oJQOfz85EeU9dYvJ+7cOXQNmjQptP3aJunPN0REAI37IknXqtXf23w7dw5bg+6+W1uDcpm7V4tWkHQAjgMeBgab2TTgBWC0uy/ZxO89D6i5gUuHb+T+P4DjN3KtP9B/U3KIZIK//gr1Ldc3Dzl8g/8vEdmwRLYC1QUecfflwF/5v97dXyM8teyX3HhS0L/+Ba+/Dr17w/DhoRje/PlxpxKRLKVxXyQF1m/z7dULhgyBBg3g00/jTiVxcvfF7n6bux9GqF9yEfAbcCuw2MzmmFlvMzsk1qAiWezKK2HUKLjjDujUKe40kmkSmVjZir+LVK2O/t0u3/U5QO1khJLClSoVWjhOnAg//BAmV+67D1x110UkuTTui6RI6dJwyy0wdix89RXUrRuK1Iu4+zfu/qC7twIqAl2AT4ArgffM7DMzu9PMDog1qEgWufdeuOsuuOgiuPTSuNNIJkpkYuUbwgw67v4rsBKolu96ZUBFDEtQ06Ywbx40axYGgWOOge++izuViGQRjfsiKda6dWjpWakSNG8eVrCIrOfuP7n7U+5+PLAz0B6YDJwGnBBnNpFs8eKLcMkl0L59mFwxizuRZKJEaqzM5H/3Y74GXGpmXxAmaHoQihtKCapYMbQAGzQIrrgCqlcPhW1btow7mYhkAY37IiVgn33Cnv5TToFzzglbfAcODKtaRNZz99WE7jxjzGwLYKeYI4lkvPVjb7168NRTYWeAyKZIZMXKUOD7qFUbwDXA78BjwDDCMvGrkppOisUMevSAmTNh551DYbzLL4fVq4v+WhGRQmjcFykh228fOlBcfjncf39o9fnjj3GnkpJmZnsUcewedQzyqP6ViGyiTz6Btm1D3asxY1REXDZPsZ+FuPtEQou29Z9/Zmb7AU2BdcDb7v5T8iNKcR18cJhcufLKsIztzTfhmWfgAO3AFZFNoHFfpGSVKhWKJlatGlau1KsX3uzrdTynfE7owlaU38zsDeAGd5+b2kgi2Wf58jCBbQbjx0NeXtyJJNMVa8WKmZUzsy5mVi//eXf/1d3HuPtYvblOD+XKhSddY8bA0qVQq1bYr63CtiKSCI37IvHp3j10DVq5EurXhwkT4k4kJehG4D3C5MqrwD3RMSE6Nxu4i1BnpQ3wtpltqJ2yiGzE779Du3bw9dehpMK++8adSLJBcbcCrQYeATRwZ4i2bUNh24YNw1Ovjh1hxYq4U4lIBtG4LxKjww8Pq1D33DMUuL333rgTSQlZBFQBDnH3Nu5+WXS0JozHewPT3b0tf3dl6x1PVJHMs24dnHoqTJ8OTz8dJq9FkqFYEyvu/hfwJVA+tXEkmXbdFV59NSwrfuWVsFXotdfiTiUimUDjvkj89twT3nknPFm9+GK46iqtQM0B1wCD3H1hwQvuPh94ALgu+nwe8DDQqEQTimSwK6+EF14IBcLbt487jWSTRIrXPg6cZmZbpSqMJN8WW4RCeDNmwA47hG5BF18clsCJiBRB475IzLbdFkaNggsugNtvD9uE1qyJO5Wk0P5AYUVpv4vuWe8DYLuUJhLJEvfdFyZULr44HCLJlEgjv/8AHYA5ZvYA8AnwW8Gb3H1qkrJJEtWoAbNmwdVXwz33wOuvw5NPQk0t8heRjdO4L5IGSpUKfxBUrAg33ADffw8jR8LWW8edTFLgW+A4YFDBC2ZmQPvonvXygB9KJppI5ho9OkymHHcc3Hln3GkkGyUysTIx38f38M+K5RadU/fvNFWuHNx9N7RpA926hW4D/frBFVeoZ7uIbJDGfZE0YQa9e4fJlfPPh+bNQ9HFHXeMO5kk2VCgr5mNJYy7H0fn9wcuBpoAffLd3waYU4L5RDLOjBlw8slw6KHw1FP6u0dSI5GJldMpXvs3SXPNm4fCtueeC716wbhxMHx42MstIpKPxn2RNHPuuaEt6CmnQOPGoZZa5cpxp5Ik6g/8CzgHaFXgmgFDgJsAzKwsMJzQRUhENmDJktDUY5ddwmS0VvpJqhR7YsXdH0thDilhO+0UlhE/8QT06AHVq8OgQdC5c3gqJiKicV8kPXXsGCZUjj02dA+aMAEOOCDuVJIMUeHw88zsPuAYYK/o0ufAy+6+KN+9fwCDSzykSIb44Qc4+uhQl2r8+LDiTyRVil281syGmVm9Qq7XNbNhyYklJcEMunSBuXNDx6DTTgvL5H7QTl0RQeO+SDo78kiYMgX++AMaNgytmSV7uPsid7/N3c+LjlvXT6qYmYrVihRh9erQ9WfJklBfZf/9i/4akc2RSFegbsA+hVzfC+i6WWkkFnvtFd6c3XwzPP98mGSZMCHuVCKSBrqhcV8kbdWsGdoxly8PzZrBu+/GnUg2l5ndU8T17QC9SxMpxF9/hQ5qU6fC449DIzUklxKQyMRKUbYB1AAwQ5UqFToGTZ8OFSpAq1Zhi9Cvv8adTETSmMZ9kZjtu294OJKXBy1banIlC1xoZldu6IKZbQ2MB2qUbCSRzHLddfDMM3DLLXDSSXGnkVxRaI0VM9sDqJLv1AFm1ngDt+4InAcsTl40iUOtWqEt87XXhj7vr70W6rDU2+hmABHJJhr3RTLP7rvD5MnQpAm0aBFeu+vXjzuVbKLrgQFmtszdn1p/0szKAeOAOoSWyyKyAUOGhAmVs8+Gnj3jTiO5pKjitd2BGwhdIRy4NjoKMuCv6H7JcOXKwV13hQraXbuGwnjXXAPXXw9lysSdTkRSTOO+SAaqXDlMrhx5ZJhcmTABGjSIO5Ukyt37m1llYKiZfevur0fdf14GGgCd3H18vClF0tO4caEdfatWoSmHGnJISSpqYuUlQhVyA4YRWrxNK3CPA6uAme7+VbIDSnyOPBLmz4eLLoJ+/cJg9cQTcOCBcScTkRTSuC+SodZPrjRpErYFvfoqHHZY3KlkE1wA7AI8b2ZHA72BxsDJ7v5yrMlE0tSsWXD88aHT6ciRULrYvW9FkqPQ/8m5+1xgLoCZ7Qm84O7zSyKYpIfttw9Fn9q1g3POCVuFBgyACy+ELZJZoUdE0oLGfZHMtttu/zu5MmGCJlcyjbv/ZWYnA68DUwmrA0919+fjTSaSnj77DNq0CbWmxo6F7dQ3S2JQ7D+N3b2v3lznro4dYcECaNoULrkEmjeHL7+MO5WIpJLGfZHMtH5yZdddw+TKtIJrziRtmFnjDR1AXeBO4FfgUeC/Ba6LCLBiBRx9NKxZA+PHh3FPJA7Fnlgxs75mtqCQ6/PM7LrkxJJ0tMsu8PLLoSjU9OmhLfNjj4F73MlEJBU07otkrt12g0mTwmt369Zha6+kpcnApI0co4DtgLPynVt/v0jO+/33sKr+iy9gzBiVK5B4JbL7rD0wsZDrE4FOwE2blUjSmhmcdVZYudKtW+gR/+KLYbKlUqW404lIkmncF8lgu+0GEyeGIvQtW8Lbb8Pee8edSgpQAXCRTbBuHXTuHFbkjRwJDRvGnUhyXSITK3sBHxZy/SPgzM2LI5li773Dk7B77gkdg6pVg4ceCluGRCRraNwXyXBVqoT2y40bh25Bb78dVrFIenD3x+POIJJp3OHSS8PD3YEDoVOnuBOJJLAVKFKhkGs7AKU2I4tkmFKl4LLL4L33YM89w6DWuTP8+GPcyUQkiTTui2S4qlVDZ7///jesXFm5Mu5EIiKb7q674L77wuTKJZfEnUYkSGRiZSFw7IYumJkB7Sj8yaZkqYMOCsvw+vYNS/GqVQstHkUk42ncF8kS9erBCy/ABx9A27bw229xJxIRSdyIEXDFFaG18h13xJ1G5G+JTKwMBeqb2WNmlrf+ZPTxMKB+dI/koDJloHdvePddqFAhVOc+91xYtSruZCKyGTTui2SRFi3gySfhnXfghBNCFw0RkUwxZQp07RrqqQwfDlskuvdCJIUSabf8MPA00IXQ8m2pmS0F/gt0BUa6+4OpiSmZonZtmD0brrwyFLStXj0MgiKSeTTui2SfE06ABx+EsWPhjDPU2U9EMsPChXDccaHO4+jRULZs3IlE/ldC83zufipwEvAK8FN0jAFOcPeTkx9PMlHZsnDbbTB1aphJPvLIsAdSy45FMo/GfZHsc845cOON8MQT4V/JXmZW1sxmmNlcM1toZn2j80+Z2UdmtsDMhplZmeh8EzP7yczmREfvfN+rVfQ1i82sV1y/k+SeZcvCaviyZWH8eNhxx7gTifxTIl2BAHD3kcDIFGSRLNOwIcydC716wd13h8J5jz0GDRrEnUxEEqFxXyT7XHcdfPop9OkD++4bis9LVloNHOXuq6LJk7fNbDzwFHBqdM/ThA5v61cgvuXux+T/JmZWChgENAeWAjPNbIy7LyqJX0Jy1y+/QJs2oTnG1Kmh05lIOtqknWlmtpWZ7WZmWyY7kGSXbbYJVbvfeANWrw6TLT17wh9/xJ1MRBKhcV8ku5iFLbtNmsDpp8Nbb8WdSFLBg/UV78pEh7v7uOiaAzOAykV8q7rAYnf/zN3/BEawkeLmIsmyZk0oUjt/Pjz3HNSsGXcikY1LaGLFzGqZ2ZvAL8CXQMPofEUze8PMmqUgo2SBo46CefPgzDPDNqHatWHWrLhTiUhRNO6LZK8tt4Tnnw9PgNu3h8WL404kRTGz3mZ2erSCpLhfU8rM5gDfARPdfXq+a2WA04D8/RwbRFuHxptZ1ejcbsBX+e5ZGp3b0M8728xmmdms5cuXFzemyP9wh/POgwkTYPBgaNUq7kQihSv2xIqZ1QDeAvYBhue/5u7fAeUIxQxFNqh8+TAwvvoq/PQT1K8P118Pf/4ZdzIR2RCN+yLZb8cdw1ZdCMvtf/gh3jxSpD7Aw8AiMzuuOF/g7uvcvQZhVUpdM6uW7/IDwFR3X79m6T1gT3c/BLgPeCnRgO4+xN3ruHudvLy8or9AZANuugmGDg1/K5xxRtxpRIqWyIqVG4FlQFWgF2AFrr9BWCYoUqiWLWHBAjj11DBo1q0barGISNrRuC+SA/bZB156CT7/HDp00AOPNNcduBhYSGh7X2zuvhKYBLQCMLMbgDzgsnz3/Lx+65C7jwPKmNnOwNfA7vm+XeXonEjSPf449O4NXbpA375xpxEpnkQmVhoBD0eD7Yaa830J/CspqSTrVagQCtmOHg3ffgt16kC/fmEvpYikDY37IjmiYcPwujxlSlh+rzbM6cndH3f3+929A2FSpFBmlmdmFaKPyxGKz35oZmcCLYGT3f2vfPfvYmYWfVyX8LfCCmAm8G8z2yuqtXUSoUOcSFK9/nooHdC0KTz8cKgHJZIJEplYKUtos7kx5Tczi+Sgdu3C6pUTTggz0/XrhwJVIpIWNO6L5JCTTw7dgoYNC1t3Jb25+7pi3LYrMMnM5hEmRya6+yvAQ0AlYFqBtsqdgAVmNhe4FzgpqnG7FugBTAA+AEa6+8Ik/0qS4+bNC6vmDjgg1H/aUuXyJYMk0m75U6B2IdePAtRyTRK2007w1FPQqROce24obNu7d+geVKZM3OlEcprGfZEc06cPvPceXHQRHHwwHH543Ilyk5ltBTQB9o5OfQZMcfeE+iq6+zzgH71U3H2DfwO4+/3A/Ru5Ng4Yl8jPFymupUuhdetQk3HcONh++7gTiSQmkRUrTwOnFegA4QBmdjlhv+YTScwmOaZ9e1i4EDp2DIWq6tcPq1lEJDYa90VyTKlS4WEQtIzGAAAgAElEQVTHnnuGBx7LlsWdKPeYWRdC/ZJxwKDoGAd8bWbdYowmkhI//xyKZ//8M4wdC7vvXvTXiKSbRCZW7gDeJSwBnEp4cz3QzL4GbgMmEiqLi2yynXeGZ54Jy/+++gpq1YL+/WHt2riTieQkjfsiOahCBXjxRfjllzC5snp13Ilyh5mdCDwGrAKuBY6Ljuuic0Oje0Sywpo14aHqokXh/f8hh8SdSGTTFHtixd3/JBS8ugL4HfgD2A/4HrgKOCZ/8SuRzdGhQxhgO3QI+70bNAirWUSk5GjcF8ld1aqFYrbTpsHFF8edJqdcA3wIVHf3Ae4+JjpuAaoDn0T3iGQ8dzj77FCw9uGHoXnzuBOJbLpEVqzg7mvdfWDUm34bd9/a3Q9x9zujolYiSbPzzjBiBDz3HHzxRVi9csstWr0iUpI07ovkrk6doFevUMj24YfjTpMzDgAedfefC15w95+ARwkT3CIZr1+/MIHbpw906xZzGJHNlNDEikgcOnUKq1WOPRauuUadg0RERErKTTdBixZw4YXw/vtxp8kJa4q47sC3JRFEJJWefRZuuAG6dg1NK0QyXUITK2ZW1syuMrNpZvZtdEyLzpVLVUiRvDwYOTKsXvnyy9A5qF+/sC9TRFJH475IbitVCp58MnTwO+GEUFxSUuo9oLuZbVvwgpmVB7oTVq2IZKyZM8MKlYYNw4o4s7gTiWy+Yk+smFkeMBMYABxIqFb+dfTxAGBmdI9IynTqFGqvdOwYZrfr1oU5c+JOJZKdNO6LCISHGyNGwGefhXoI7nEnymrnAJWB+WZ2pZm1jY6rgLmEArZvmVnj/EesiUUS8PXXYRX6LrvACy/AVlvFnUgkORJZsXI7cBBwGVDR3Wu5ey2gInA54Y327cmPKPK/1ncOevFF+OYbOPTQMMny559xJxPJOhr3RQSARo3CtqBnn4UhQ+JOk73c/QOgEbAMuBV4KToGAHsCtQkd2SZFx+ToX5G099tv0K5d6Dj28sth0lYkW5RO4N62wFB3vzv/yahrxEAzqwq0T2Y4kcIcdxw0bgyXXBK2Bb30Ejz6aNgmJCJJoXFfRP5fz54wZUroElSvHtSoEXei7OTuc4HDzWxfYF8gD9VFlAz311+hnsr778OYMaHzmEg2SWRiZUvCvs+NmQWcuHlxRBKz444wfHjY933OOeGN3lVXhWJYWloostk07ovI/9tii/CaW7NmeN2dPRu22y7uVNnL3RcDi+POIZIMffvCqFFw++1wzDFxpxFJvkRmv2cCtQq5XhuYsXlxRDbNMceEzkFduoSWzDVrwvTpcacSyXga90Xkf1SsGLbjfvopnHtu3GlEJBOMGAE33gjdu8Pll8edRiQ1EplYuRzoZGYXmtn/r3Qxs9JmdjHQIbpHJBYVKsCwYTB+PKxaBYcdFgbv336LO5lIxtK4LyL/0LhxWBn69NPhEBHZmBkzwoRKw4bw4IPqACTZy7yYpd3N7E1gd2Bv4Gfgs+jS3kB54FNgaYEvc3dvWvB71alTx2fNmrWpmUWK9PPPYS/4Qw/BPvvA0KFwxBFxp5JsYmaz3b1O3DlSSeO+iGzM2rXhdXXhQpg7F/bcM+5EqZcL436yaezPbUuXhiYTZcuGCRYVq5VMVNyxP5EVK3sTarJ8CawEdoyOldG5MsBeBY69E4stkhzly4dZ8UmTQlvIJk3g/PNDFXIRKTaN+yKyQaVLwxNPwLp1YRvuunVxJxKRdPLrr6Gt8qpV6gAkuaHYxWvdvUoKc4ikRJMmMG8eXH893H03jB0b2kS2bBl3MpH0p3FfRAqz995w331hmf8dd4SVoiIi6gAkuUit2yTrbbMN3HUXvPNO+LhVq/Am8Mcf404mIiKS2bp2hU6dwgOM9wrrISYiOaNPH3j+eXUAktyyyRMrZlbBzE4xsyvNrE0yQ4mkQoMG4U3fNdeE5csHHQQvvRR3KpHMoXFfRAoyg8GDwzL/U05RwXiRXPf009CvX3iIedllcacRKTmFTqyYWXszG21mFQucrwUsBJ4AbgXGmNlEMyuTuqgim69sWejfH2bOhEqVoH17OOEE+PbbuJOJpAeN+yKSqB13hOHD4aOP4Npr404jInGZNg1OPz10DnvoIXUAktxS1IqVE4Dd3f27AucfBXYFngEuAt4AjgLOT3pCkRSoWTNMrvTvD6NHw4EHhjeFxWySJZLNNO6LSMKaNoULLoB77oG33oo7jYiUtC++gOOOg912C9uAttwy7kQiJauoiZXawKT8J6KnlgcDY9z9VHe/H2gJvEd4Qy6SEcqUCduC5syBAw4I+8Rbt4Yvv4w7mUisNO6LyCYZMACqVAlPrLUlSCR3/PILtG0Lq1fDK6/AzjvHnUik5BU1sVIJWFzgXCPACcvBAXB3B54HDkxqOpEScOCB4enavfeGf6tWhUGDQkVzkRykcV9ENsm228KwYbB4sbYEieSKdetCfaVFi2DkyPC+WiQXFTWxsqGdcYdG/75d4Px/gW02O5FIDEqVggsvhAULQpHbHj3giCPCfnGRHKNxX0Q2WZMm2hIkkkt69gyrVO65B1q0iDuNSHyKmlj5AqhZ4Fwj4Ct3L1juc3vgh2QFE4lDlSowYQI8+miYZDnkkLC0ec2auJOJlBiN+yKyWbQlSCQ3DB0Kd94ZJlMvuCDuNCLxKmpiZQLQ2cyOMbOtzewSYHdgzAburQWoOoVkPDPo1g0++ADatIGrr4Z69UKrZpEcoHFfRDZL/i1B110XdxoRSYXJk+Hcc8MqlbvvjjuNSPyKmli5HVgFjAZ+Ae4CfgLuyH+TmZUF2gJTU5BRJBa77BKqmo8aBcuWQd260KsX/P573MlEUkrjvohstiZNwh9d99wDs2bFnUZEkumTT6BjR/j3v+HZZ6F06bgTicSv0ImVaNn3ocADwGvRv7XcveATynrAO8BzqQgpEqeOHcPqla5d4dZbw/agKVPiTiWSGhr3RSRZBgyASpXgrLNg7dq404hIMvz4IxxzTFjh/fLLUKFC3IlE0kNRK1Zw9y/d/UJ3P9rde7j7kg3cM8Xd27r7jNTEFInXDjuEfaSvvx6qnzdpAuecAz/9FHcykeTTuC8iybD99nDffTBnjrYKiGSDNWvg+ONhyRJ48UXYZ5+4E4mkjyInVkTkb02bwvz5cMUV8MgjcNBBMHp03KlERETSU4cO0K4d9O4d/hgTkczkDhddBG+8AUOGQKNGcScSSS8bnVgxszPMLOGJFzMrZWZnbl4skfS19dZw++0wfTrsvDMcd1yYvf/vf+NOJrJ5NO6LSLKZwf33Q6lScP754Y8zEck8gwbBQw/BVVeFJg8i8r8KewN9J/ChmfUws52L+kZmVsnMLgU+IhQ/FMlqdeqEgnz9+4c9pgcdFNo0602jZDCN+yKSdLvvDjffDK++Cs88E3ea3GJmZc1shpnNNbOFZtY3Ov+UmX1kZgvMbJiZlYnOm5nda2aLzWyemdXK9726mtkn0dE1rt9JSt7EiXDJJWH12S23xJ1GJD2Zb+SvQDPLA/oD3YG/gFnADOBT4AfAgB2BfwP1gRrRlw4Ferv78o390Dp16vgslYiXLPLRR3DmmfD222G70ODB2nea7cxstrvXiTtHMmncF5FUWbcODjsMvvgivGZuv33ciRKXieO+mRmwjbuviiZP3gYuJozl46PbngamuvuDZtYauBBoTShSfo+71zOzHQmvCXUAB2YDtd39x8J+vsb+zPfxx1CvXpggfecd2G67uBOJlKzijv0bbY4VvUE+O5rZPhfoRBiIN2QhcBPwsLt/swl5RTLa/vuHTkEPPxyWSB58MPTpA5ddphZ0kjk07otIqpQqBQ88AIceCjfcoGK2JcXDE9RV0adlosPdfdz6e8xsBlA5+vRYYHj0de+aWQUz2xVoAkx09x+ir5kItAK0BimL/fgjtG0b3suOGaNJFZHCFKcr0Nfufr27HwjsAhwFnAicABwJVHT3g939Rr25lly2xRahU9CiRdCyJfTsCXXrwnvvxZ1MJDEa90UkFWrXhnPPDZ2C5s6NO03uiOpgzQG+I0yOTM93rQxwGvBqdGo34Kt8X740Orex8xv6eWeb2Swzm7V8+UYXMkqaW7sWTjwxFJ1+4QWoUiXuRCLpLaEihe7+nbtPdvfn3H1U1G7z+1SF+7/27jtOqvL64/jnAPaGBUvEGg32FmzRKGDDEtFEo8bee4moPyuKGo0l2IVgJ7bYQSxYAkQRCzZEEcWOFVFANIrA+f1xnpFhmF12l9m5Mzvf9+s1r925986dc+/O3rn33Oc5j0g1Wn75GILu/vvh888juXLqqfDDD1lHJtJ4Ou6LSCldeCEssQQceyzMmJF1NLXB3ae7+wZEq5RNzGydvNnXE92Aninh+/V1947u3rFdu3alWq2UWffuUVuld2+NACTSEBpuWaSZ/PGPMHo0HHIIXH55dA966qmsoxIREcnOEkvApZdGrYZ//SvraGqLu08EBhNdeDCzc4F2wMl5i30KrJD3vH2aVtd0aYH69oWrr4a//hUOPTTraESqgxIrIs2obdv4choyJPqXb7ddDFE3YULWkYmIiGTjwANh882jNefEiVlH07KZWTsza5t+XwDYjhj97TBgB2Afd89vOzQAOCCNDrQZMCl1+RwEbG9mi5vZ4sD2aZq0MEOHRouyrl0jCSoiDaPEikgZbL01jBwJZ54Jd9wBa64Jd96poZlFRKT2tGoVhWwnTIAePbKOpsVbDhhsZiOBl4gaKwOBPsAywHAze83Mcn+JR4H3gbHADcAxAKlo7QVpHS8B5+cK2UrL8f778Kc/wWqrwd13awAGkcbQv4tImcw/P/ztb1EI7PDDYd99oV+/6Lu6yipZRyciIlI+G2wARxwRCZajjoK11so6opbJ3UcCGxaZXvQaII0GdGwd824Gbi5pgFIxJk+GXXeN2kcDBlTnkOgiWVKLFZEyW289eO656Ls6bBissw784x9RfV1ERKRWnH8+LLwwnHyyWnCKZGn69Ljh9/bbcO+9sPrqWUckUn2UWBHJQOvWcPzxMTTzNtvAKafApptqaGYREakd7drBuefCoEHw6KNZRyNSu848EwYOjJt+22yTdTQi1UmJFZEMrbAC9O8fdwc++ww23jiSLN9/n3VkIiIize/YY6FDh2i1MnVq1tGI1J5+/aJI7dFHwzHHZB2NSPVqcI0VM/vPHBZx4H/Ax8ATQP/UT1NE6mEGe+wB224Lp58e3YLuvx/69IEddsg6OqllOu6LSHObd17o1Qt23hmuuy6GdxWR8njuuaj717kzXHVV1tGIVLfGFK9dFViAGO8eIDdAXtv0czzRAmYn4EhgmJnt6O669y7SAG3bRjJl332joF/XrvCXv8AVV8DSS2cdndQoHfdFpNnttFN85/XsCfvvD0stlXVEIi3fRx/B7rvDiivCfffBPPNkHZFIdWtMV6BOwA/AZcAy7r6Euy9BDNV2OfA90BFYCugFbAloED2RRvr97+G116Lf+X33wRprwE03qbCfZKITOu6LSBn06gXffQcXXJB1JCIt35QpMQLQTz/Bww/DEktkHZFI9WtMYuUKYJi7/5+7j89NdPfx7n4aMBy4wt2/cfdTgUeAP5U2XJHaMN98cN558PrrsO66cNhh0UxzzJisI5Mao+O+iJTFmmvGd93118PYsVlHI9JyzZgRLcNGjYK7744beCIy9xqTWOkCPFPP/GfSMjlPAe2bEpSIhDXWgMGD4cYbI8my3nrRVPqnn7KOTGqEjvsiUjY9e8aNhTPOyDoSkZbrnHPgoYeilVjXrllHI9JyNHZUoPpymmsAlvd8BlHUUETmQqtWcOih8Pbb8Kc/RUuWDTaAZ+q73BUpHR33RaQsll0WTj01usE+/3zW0Yi0PHfcARddFAVrTzgh62hEWpbGJFaeAo42s70LZ5jZPsBRwJN5kzcCPpyr6ETkF8ssA3feCY89Bj/+CFttFV+M336bdWTSgum4LyJl1b17JFhOOUW1xURK6YUX4kbd1lvDtdfGqJQiUjqNSaycTIwAcYeZjTOzIekxDrgd+BroDmBm8wMrAf1KHbBIrevaNfrFnnoq3HJLdBe66y6dgEqz0HFfRMpq4YWjS9CwYdFdQUTm3iefwG67wfLLR4uweefNOiKRlqfBiRV3/whYH/gHMBnYND2+S9PWT8vg7j+6exd3v6L0IYvIQgvBpZfCiBGw0koxLPOOO8L772cdmbQkOu6LSBYOOSSK2Z55JkyblnU0ItXt+++hW7f4OWCAhjMXaS6NqrGSRn44zd3XcvcF0mPNNG1CcwUpIsVtsAEMHw7XXAPPPQdrrw1//zv8/HPWkUlLoeO+iJRbmzZw4YVRW+xf/8o6GpHqNWMGHHggvPZajAC09tpZRyTScjW2eK2IVJjWreG442D0aNh55xhNYaONItEiIiJSjXbfHTbeOAq2ayQ8kabp2RPuvx8uuwx22inraERatkYlVsxsITPraWYjzWxKeow0s/PMbKHmClJE5izXb3bAAJg0CbbYAo46SsVtZe7ouC8iWTCL0Us+/hj69Mk6GpHq8+9/w/nnw8EHw8knZx2NSMvX4MSKmS0BvAicAywDvJoeywA9gBfTMiKSoT/8Ad56K75Eb7gh+qnffbeK20rj6bgvIlnadlvo0gX+9jf47rusoxGpHsOHw0EHwZZbQu/eGgFIpBwa02LlfGAN4DjgV+7+e3f/PfAr4FigA3BeySMUkUZbeGH4xz+iuO0KK8A++6i4rTSJjvsikqmLLoLx4+HKK7OORKQ6jB0Lu+4K7dvDgw/CfPNlHZFIbWhMYmVX4EZ3v97dp+cmuvt0d+8N3AzsVuoARaTpNtwQnn8err46hq5ce+04SZ06NevIpErouC8imdp00xgm9vLL1bVVZE6+/jpqqbjDY49pBCCRcmpMYiXXDLwur6RlRKSCtG4Nxx8foyvsvDOcdVYkXJ55JuvIpArouC8imevZEyZPhl69so5EpHL9+GMkIT/+OOrtrbZa1hGJ1JbGJFa+BDasZ/6GaRkRqUC54rYDB8L338NWW8Ghh8IEDZgrddNxX0Qyt956sOee0R1I31kis8sNqzxsWAxR/rvfZR2RSO1pTGLlYeBQMzvSzH55nZm1MrMjgEOAAaUOUERKa+ed4c034bTToF8/WGMNuO02FbeVonTcF5GKcO65cVPg8suzjkSk8pxxBtxzTwyrvOeeWUcjUpsak1jpAbwPXA98ZmZDzWwo8BnQO807t/QhikipLbQQXHIJvPIKrL56VI7v3Dm6C4nk0XFfRCrC2mvD3nvDNddEMVsRCX36wKWXwtFHQ/fuWUcjUrsanFhx9wlAR+DvwARg4/T4GrgY2DgtIyJVYt114dlnoW9feP31aG59zjnwv/9lHZlUAh33RaSS9OgR30+XXZZ1JCKV4dFH4dhjozXy1VdrWGWRLDWmxQruPtndz3L3td19wfRYx93PdvfJzRWkiDSfVq3g8MNhzBjYay+48MJIuAwalHVkUgl03BeRSrHGGvCXv8C118JXX2UdjUi2XnkF/vxn2GADuPtuaNMm64hEalujEiv1SX3w3yrV+kSkvJZeOgqePfVUjCTUtWskWj77LOvIpFLpuC8i5Xb22TH6iUYIklr28cewyy6w5JIxKMHCC2cdkYiULLECLAV0KOH6RCQD22wDI0fC+edD//5xh/Dqq2HatKwjkwqk476IlFWHDpH0v+46jRAktWnSpOj68/330RVoueWyjkhEoLSJFRFpIeabL2qtjBoVQ/adeCJssgm8+GLWkYmISK076yyYMiWGXxapJVOnwp/+FIMNPPBAFHUWkcqgxIqI1Gm11eCxx2IIvy++gM02g2OOgYkTs45MRERq1TrrxMXl1Vfr+0hqx/TpsN9+8PTTcOON0cJYRCqHEisiUi8z2HPPuDtywgnwz39G96A77gD3rKMTEZFadPbZMHlyDL8s0tK5w5FHwr33wuWXw4EHZh2RiBRSYkVEGmTRRaPZ9UsvwYorxl2TbbeNhIuIiEg5bbAB/OEPcMUV8N13WUdTucxsfjN70cxeN7M3zaxnmn6cmY01MzezpfKW72Rmk8zstfTokTevq5mNSa87PYvtqUXu0L073HRTdNPu3j3riESkmHoH5jKzkxuxri3mMhYRqQIbbQTDh0PfvnDmmbDeenDqqdHnfcEFs45O5paO+yJSLc46K7qo9u2ri816/AR0cfcpZjYP8KyZPQYMAwYCQ4q85hl33yV/gpm1Bq4DtgPGAS+Z2QB318hwzez88yOBeMIJ0LNn1tGISF3mNOL55Y1cnzoGiNSA1q3h6KOjj/tpp8FFF0XXoGuuiTuIUtV03BeRqrDpptC5cwy9fNxxUXhdZuXuDkxJT+dJD3f3VwHMrKGr2gQY6+7vp9fdDXQDlFhpRldeCeedBwcdFMmVhv+5RKTc5pRY6VyWKESkKi29NNx6KxxySBS13XXXeFx1Fay8ctbRSRPpuC8iVeOMM2D77aFfPzj88KyjqUyptcnLwGrAde7+whxesrmZvQ58Bpzi7m8CywOf5C0zDti0jvc7AjgCYMUVV5zL6GvXzTfDX/8aN7FuuAFaqYCDSEWrN7Hi7kPLFYiIVK+ttoJXX515Z2WttWb2A5533qyjk8bQcV9Eqsm220LHjnDppZHkb90664gqj7tPBzYws7bAg2a2jruPqmPxV4CVUtehnYCHgNUb+X59gb4AHTt2VKvGJrj33kgU7rBDtAhuM6db4SKSOeU+RaQk5pknaq2MHg077hj1V9ZfHwYPzjoyERFpqcyi1crYsXDffVlHU9ncfSIwGOhazzKT3X1K+v1RYJ5U3PZTYIW8RdunaVJijz0G++4Lv/sdPPCAuriJVAslVkSkpFZcEe6/Hx55BKZOhS5d4gTh88+zjkxERFqi3XaDNdaAv/89RlCRmcysXWqpgpktQBSfrXM8PzNb1lLhFTPbhLhWmAC8BKxuZquY2bzA3sCA5o6/1vz3v/DHP8I668DAgRoUQKSaKLEiIs1ip51g1Cjo0SMSLR06RFehadOyjkxERFqSVq2ixeRrr8FTT2UdTcVZDhhsZiOJ5MiT7j7QzE4ws3FEy5ORZnZjWn4PYFSqsXI1sLeHacBxwCBgNHBPqr0iJTJiBOyyS9SoGzQIFlss64hEpDHMM0jtd+zY0UeMGFH29xWRbIwdC8cfD48/DuuuC9dfD1tumXVUc8fMXnb3jlnHUS103BeR5vTTT7DKKrD22vDkk83zHjruN56O/Q0zahR06gQLLwzPPgvt22cdkYjkNPTYrxYrItLsVlsNHn00+gpPmgS//z0ceCB8+WXWkYmISEsw33xw4onRYuXVV7OORqThnnkmzovmnTc+v0qqiFQnJVZEpCzMYPfd4a23otDgXXdF96Brr4Xp07OOTkREqt2RR8Iii8Bll2UdiUjD3HcfbLcdLLMMPPdc3IgSkeqkxIqIlNVCC8FFF8Ebb8DGG0cXoY03huefzzoyERGpZm3bwhFHwD33wIcfZh2NSP2uuQb+/GfYaCMYNixqq4hI9VJiRUQy0aEDPPEE/Pvf0SVo883h0ENh/PisIxMRkWp10knRQvKKK7KORKS4GTPgtNPghBOgWzd4+mlYcsmsoxKRuaXEiohkxizu1rz9dozo0K8f/OY3cN116h4kIiKN17497LMP3HwzTJyYdTQis5o6FfbfP7qrHXNMdAVaYIGsoxKRUlBiRUQyt8gicOmlMHIk/Pa3cNxx0LFj9DcWERFpjL/+FaZMgRtvnPOyIuUyeTLstBPceWd0ib72WmjdOuuoRKRUlFgRkYqx5poxTOY998DXX8MWW8BBB2n0IBERabgNN4TOneHqq+Hnn7OORgQ++wy22gqGDoVbb40i/mZZRyUipaTEiohUFDPYc08YPRpOPz3u7HToECfI06ZlHZ2IiFSDk0+GTz6B++/POhKpdaNHRx25996DRx6BAw/MOiIRaQ5KrIhIRVp4Ybj44hg9aNNN4cQTo5vQM89kHZmIiFS6nXaKpPw//gHuWUcjtWro0Gh9+9NP8fv222cdkYg0FyVWRKSidegAjz8edx0nToymtPvtF81qRUREimnVKmqtjBgRQ9mKlNNHH8W5SqdO0K4dDB8ewyqLSMulxIqIVDwz+OMfoznt2WdHFf0OHaKq/tSpWUcnIiKVaP/9YfHF4aqrso5EsjRmDAwYUJ5RoiZNim7MHTrEDaGzzoKXXoJVVmn+9xaRbCmxIiJVY8EF4YIL4M03ozDhaafBeuvBoEFZRyYiIpVmwQXh8MPhwQfh44+zjkay8PTT0VKkWzdYcknYbLO4QTN0aHTPKZWff45RflZbDS65BPbaC955By68EBZdtHTvIyKVS4kVEak6v/513H169FGYMQO6doXdd4cPPsg6MhERqSTHHhs/r7su2zik/B57DHbeOc4ZHn88EiqtWsHf/x5ddJZYImrx9OoFI0c2rRaPO/TvD+usA8cfD+uuCy+/DLfdBiusUPJNEpEK1ibrAEREmmrHHaFLF7jiirgrtNZa0Yrl//4v7lSKiEhtW3HFSLz37Qs9esBCC2UdkZRD//4xwuC668ITT0RrlR12gJ49o7vOkCHw1FPw5JPQvXu8Zpll4Pe/h5VWgl/9Kh7LLx8/l1tu9vOKESPitf/9L6yxBjz8cCRyNIyySG1SYkVEqtp880V/5v32g1NPhfPPjztFvXrFybROcEREatuJJ0Ztrn/9C446KutopLndcw/su2+MJPj449C27azzF1ssugZ16xbPP/kkugw9+SQ8/zwMHAg//jj7etu2nZlwad06uiG3awe9e8Nhh0EbXVWJ1DTzDMag69ixo48YMaLs7ysiLd/QodEc9403YNtto2jhWmuV/n3M7GV371j6NbdMOu6LSFbcoX2Uh+AAACAASURBVGPHuFgeNarpCXcd9xuv3Mf+fv3g4INjiONHHoFFFmn8Otyj0O1nn9X9mDAB9tgjbuyohopIy9bQY79yqyLSomy9NbzyCvTpA+ecE8Vtjz8ezj139rtWIiLS8pnBCSfAQQdFF5DOnbOOSJrDDTfAkUdGF+H+/Zve7cssRpNafHFYe+3SxigiLZeK14pIi9OmDRx3HLz7bjTPveoq+M1v4KabotitiIjUlr32ijob11yTdSRSn2nT4tFY114LRxwRxewffli1dESk/JRYEZEWa6mlouXKiBGRWDnsMNh00+hDLSIitWP++WPo5f79NfRyJXv66RitZ8cd4eKL4bnnYOrU+l9z2WXRMrVbtxhae4EFyhOriEg+JVZEpMXbaCN45hm4/Xb49FPYfPNoEv7FF1lHJiIi5ZIrXNunT7ZxSN2WXTaK0X/8MZx5ZtRKads2uvf07AmDB8MPP8xc/oILYjTAvfaCe++NgvYiIllQYkVEaoJZjBIwZkwUm7vzzmjFcvnlc74bJiIi1W+llWDXXaMWR7FRXyR7668P118Pb74JX30FDzwQXXwmTozESpcukWjZcssoHtujBxxwANxxB8wzT9bRi0gtU2JFRGrKIotE8+I334xCtz16wJdfZh2ViIiUw3HHwTffwLBhWUcic9KuHey+O1x5ZRSl/+abGAr5pJOiDsvAgXD00XDLLTH8sYhIljQqkIjUpNVXjwJ3778PK6yQdTQiIlIOXbrABx/AiitmHYk0Vtu2sPPO8YAoRt9Kt4hFpELocCQiNW3VVbOOQEREysVMSZWWQkkVEakkOiSJiIiIiIiIiDSREisiIiIiIiIiIk2kxIqIiIiIiIiISBMpsSIiIiIiIiIi0kRKrIiIiIiIiIiINJESKyIiIiIiLZCZzW9mL5rZ62b2ppn1TNOPM7OxZuZmtlTe8mZmV6d5I81so7x5B5rZu+lxYBbbIyJSqdpkHYCIiIiIiDSLn4Au7j7FzOYBnjWzx4BhwEBgSMHyOwKrp8emQG9gUzNbAjgX6Ag48LKZDXD3b8uzGSIilU0tVkREREREWiAPU9LTedLD3f1Vd/+wyEu6Af3S654H2prZcsAOwJPu/k1KpjwJdC3DJoiIVAUlVkREREREWigza21mrwFfEcmRF+pZfHngk7zn49K0uqYXe78jzGyEmY0YP3783AUvIlIllFgREREREWmh3H26u28AtAc2MbN1mvn9+rp7R3fv2K5du+Z8KxGRipFJjZWXX375azP7aA6LLQV8XY54mkjxNV0lxwaKb25VcnyljG2lEq2nJjTwuF8qiwGTKmR9TXltQ1/TkOXmtMyc5lfy//OclPpzUO73q9TPXUOXnZvPXqV+7qr6uO/uE81sMNGFZ1Qdi30KrJD3vH2a9inQqWD6kDm9Z8Gxv1L/rqVWK9sJtbOt2s6WpzHb2qBjfyaJFXefY/razEa4e8dyxNMUiq/pKjk2UHxzq5Ljq+TYWrqGHPdLxcz6uvsRlbC+pry2oa9pyHJzWqYB86v2f6bUn4Nyv1+lfu4auuzcfPaq+XNXacysHfBzSqosAGwHXFLPSwYAx5nZ3UTx2knu/rmZDQIuMrPF03LbA2fM6f3zj/218netle2E2tlWbWfL0xzbqq5AIiLS0jxcQetrymsb+pqGLDenZUq9rypJubetVj53DV22lj97lWQ5YLCZjQReImqsDDSzE8xsHNHyZKSZ3ZiWfxR4HxgL3AAcA+Du3wAXpHW8BJyfpomICGDunnUMRVV6xkzxNV0lxwaKb25VcnyVHJtIJdL/jGRBn7uWqVb+rrWynVA726rtbHlqrcVK36wDmAPF13SVHBsovrlVyfFVcmwilUj/M5IFfe5aplr5u9bKdkLtbKu2s+Up+bZWbIsVEREREREREZFKV8ktVkREREREREREKpoSKyIiIiIiIiIiTVSRiRUz62pmY8xsrJmdnsH7r2Bmg83sLTN708xOTNOXMLMnzezd9HPxNN3M7OoU70gz26hMcbY2s1fNbGB6voqZvZDi+LeZzZumz5eej03zVy5DbG3N7D4ze9vMRpvZ5pWy/8zsr+nvOsrM7jKz+bPcd2Z2s5l9ZWaj8qY1el+Z2YFp+XfN7MBmju+y9LcdaWYPmlnbvHlnpPjGmNkOedOb5f+6WHx587qbmZvZUul52fefiIhIrcv63L6UquU6oVQq+XqjVCr5uqXUKu06qFQq4Xqq4hIrZtYauA7YEVgL2MfM1ipzGNOA7u6+FrAZcGyK4XTgaXdfHXg6PSfFunp6HAH0LlOcJwKj855fAlzh7qsB3wKHpumHAt+m6Vek5ZrbVcDj7r4GsH6KM/P9Z2bLAycAHd19HaA1sDfZ7rtbga4F0xq1r8xsCeBcYFNgE+Dc3MGjmeJ7EljH3dcD3gHOSHGsRezPtdNrrk9fyM35f10sPsxsBWB74OO8yVnsPxERkZpVIef2pVQt1wmlUsnXG6VSkdctpVah10GlcisZX09VXGKF2Iix7v6+u08F7ga6lTMAd//c3V9Jv39H/HMtn+K4LS12G7Bb+r0b0M/D80BbM1uuOWM0s/bAzsCN6bkBXYD76ogvF/d9wDZp+eaKbTFgK+AmAHef6u4TqZz91wZYwMzaAAsCn5PhvnP3/wLfFExu7L7aAXjS3b9x92+JxMdsyYZSxefuT7j7tPT0eaB9Xnx3u/tP7v4BMJb4n262/+s69h/EF8BpQH6F7rLvP5GWwsx2M7Mb0t2r7bOOR2qHma1pZn3SHeWjs45HGi3zc/tSqobrhFKp5OuNUqmC65ZSq6jroFKphOupSkysLA98kvd8XJqWidTkaUPgBWAZd/88zfoCWCb9nkXMVxIXjTPS8yWBiXkXu/kx/BJfmj8pLd9cVgHGA7ekpoM3mtlCVMD+c/dPgcuJVgyfE/viZSpn3+U0dl9l+X9zCPBY+r0i4jOzbsCn7v56wayKiE+kUhRrOpumz9Zs390fcvfDgaOAvbKIV1qORn72Rrv7UcCfgS2yiFfmSov9jq3g64RSqeTrjVKp2OuWUqui66BSKev1VCUmViqGmS0M3A+c5O6T8+d5jFOdyVjVZrYL8JW7v5zF+zdAG2AjoLe7bwh8z8ymV0B2+y815+pGHER/BSxEhbdMyPKzNidmdhbRJPaOrGPJMbMFgTOBHlnHIlIFbqXgGNiAZvtnp/kic+NWGvHZM7NdgUeAR8sbpkhxlXqdUCpVcL1RKhV73VJq1XgdVCrl+BtWYmLlU2CFvOft07SyMrN5iIPlHe7+QJr8Za6pV/r5VZpe7pi3AHY1sw+J5pRdiL6BbVOzrsIYfokvzV8MmNCM8Y0Dxrn7C+n5fcQBqxL237bAB+4+3t1/Bh4g9mel7Lucxu6rsv/fmNlBwC7AvulgVSnx/Zr4wng9/Y+0B14xs2UrJD6RilFH09mizfZTsbdLgMdyzeBFmqoxn720/AB33xHYt7yRSgm0uO/YCr9OKJVKv94olUq+bim1arkOKpWyXk9VYmLlJWB1i+rE8xIFdQaUM4DUd+wmYLS798qbNQDIVQc+EOifN/2AdNK5GTApr9lRybn7Ge7e3t1XJvbPf9x9X2AwsEcd8eXi3iMt32wZO3f/AvjEzDqkSdsAb1EZ++9jYDMzWzD9nXOxVcS+y9PYfTUI2N7MFk/Z6O3TtGZhZl2JpqG7uvsPBXHvbVFFfBWiKNSLlPH/2t3fcPel3X3l9D8yDtgofS4rYv+JVLi6msIeT5yU7WFmR2URmLR4RT97ZtbJYgSHf6IWK9Uo83P7Uqr064RSqfTrjVKp8OuWUquW66BSKe/1lLtX3APYiRhp5D3grAzef0uiqdBI4LX02InoU/Y08C7wFLBEWt6IpqvvAW8QlZbLFWsnYGD6fVXiInYscC8wX5o+f3o+Ns1ftQxxbQCMSPvwIWDxStl/QE/gbWAU8C9gviz3HXAX0c/xZ+Ik8tCm7Cui1snY9Di4meMbS5z85v4/+uQtf1aKbwywY970Zvm/LhZfwfwPgaWy2n966FHpD2BlYFTe8z2AG/Oe7w9cm3WcerS8hz57tfNornOAjLalaq4TSrjNnajA640Sbl/FXrc0w7ZW1HVQCbcr8+spSysQERGRGpSKLw70GHoRM9scOM/dd0jPzwBw94uzilFaJn32RESkpajErkAiIiKSnRbVbF+qij57IiJSlZRYERERqVFmdhcwHOhgZuPM7FCPIRWPI/oVjwbucfc3s4xTWh599kREpCVRVyARERERERERkSZSixURERERERERkSZSYkVEREREREREpImUWJHMmdntZjYt6zgay8yeNbOxzbTucWb2VHOsW0RERESkWpnZh2Y2pJnW7WZ2a3OsW1o2JVaqkJl1Tf/0FxaZt1ma95OZLVhk/uNmNsPMlipPtNXNzA4xsxOyjkNEREREpDmYWad0/ZB7TDezb81slJndlq49rMwxnWRmB5XzPUXmRpusA5AmeRaYBnQqMq9zmjcv8Dvgl1YPZtYG2BIY5e5fN3+YLcIhwLLA1VkHIiIiIiLSjO4CHgUMWAToAOwGHAA8ZWZ7uvvEMsVyEvAhcGuZ3k9kriixUoXcfYqZvQRsbGYLuvsPebM7AU8CG6bf87uTbAwsBAwpT6QiIiIiIlIlXnH32/MnmNnJwKXAyUTiZccsAhOpdOoKVL0GE61StshNSC1StgCGpkfngtd0yntt7jVrmllvM3vLzL4zsx/MbISZHZL/QjM7PjUN3KkwEDNrbWafm9mIgumbmFl/M5uQuiaNMbMzzKx1QzbQzJY3sz5m9omZTTWzT9PzpQqWuzDFtpqZXZKW+8nMXjWzrkXWu5CZXWlmX6TtHW5mnQtrvZjZuLQ/f13QPHLLgvW1N7N/pyaT35vZY2a2WkO2sTHM7Ndm9m6qv7JuwbZ3MLOr0zZ9b2ZPmtnqaZk90r74n5l9YGaHljo2EalcZra0mU0ys8OzjqU55B0H22cdS2M0pr6YmXVM3Xi3mPPSc1xXs9UHE5GWx92nu3t3osV81yLnwYul8++x6fx7vJndZWarFix3UDpWb2tm55nZR2n5kWa2d8GyDqwEbF1wDr5ywXJrmNkj6RpmkpndZ2bLlnofmNlG6Rz7LTNbMU27NcW0ZPr96xTHQ7kYzOwIMxttZj+a2dtm1q3UsUnlUGKleuWSI53ypuVapOQSKxub2UJ58zsBnublbEMkDwYApwDnADOAm8zs1Lzl7gJ+JpoCFtqe6C5zW26Cme1KHIBXBS4DTgBeAP4G3F5kHbMws1WAEcAf0/LHAncA+wLDzGzRIi+7HdgsvV+PFFN/M1shb70GPACcCAxL2/wc0B9Yr2B9xwPvAF8C++c9xuQtswjwX+An4EygN7FPHzKzkv1/mVnHFOdUYHN3f6NgkduBtYn924v4mw6y6Jt6NbHNpwKTgRvNbLNSxSZS62xm3/RTso6lDhcC44Fb4JcTUTez2Y7FZrZs3gnsbAnilNx2M/tt84dd/czsj2bWY27X4+4jgIFAr/Q9JiJSbjelnzvnJpjZYsT56THAI8S587VAF+AFM1upyHouAfYGrifO1+cF7rJZ66nsD3wNvM2s5+Dj85ZZnmiF/zFxjnsncd3Qr+mbODsz24G4dnoP2NLdPy5Y5HFgsbQtNwC7AA+m66hTieuj04ntvC9d40hL5O56VOEDWIC4mH8ub9oZwHdEF681iSTK9mleG2AK8FrBehYqsu7WRFLkG6B13vQHgf8BixYsfxdxwb9Uer4gceAbnP/6NO/UFNeWedNuB6YVLPcI8AXwq4LpmwLTgbPzpl2Y1vkQYHnTN0/TL8ibtmua1rtgvbnphXE8C4yt42/wbHrNyQXTz0jTt5mLv+844Kn0e9f0t3sGWLxgubq2/eQ0fRKwfN70ZdPn5l9Zf4b10KOlPJiZtD4l61iKxNaeqLt1YsH0z4FPiiy/T9qWn4HDiswfA3wLtMp62/Jiyh0H22cdS5HYZvt+a8i8OpbvkrZzh7mMqc7vNT300KM2Hw35HgM2SsvcnzftqnRtsH7BsisRN/NuzZt2UHr9R8BiedMXS9O+ARbIm/4hMKSOWD5M6/pzwfTr0vQOc7EvPBc3kcyZms6zFyhY7ta07HUF03ul6R+Td81E3MB14OKs/956NM9DLVaqlLv/j2gB0jGvVUonItEyzd1HA18xs0VLrjXL4IL1fJ/73czmN7MlgcWBQennb/IWvw2YH/hz3msWBboBj/jMgrg7AEsBNwOLm9lSuQdREAuilUtRKYYdiYPY1ILXvwd8UMfrr/J05ErbNpw42K+et8wf0s9eBfthAPBuXTHVYxqRmc/3n/RzdeaSmR0IPAw8AWzn7t/Wsegs204kYQAedPdPcxPd/QtgbCliE5GqcCRxIndXwfTBQPsirVI6AaOJFoOd8meY2XLEd8Iz7j6jOYKVeg0hku5HZRyHiNSmyennovBLK/B9iZbbnxacr38PPE/x8/Xe7j4p9yT93oe47ujUiHg+c/d7CqaV8hz8dOLa52bgT+naq5grC57nzsH7uXtun+HuI4l9qHPwFkqJleo2GJgH2NJmra+S819m1lnplH4OyV+BmS1iZr3M7BMiCfE10drk/LTI4nmLP5Lm53cH2pNoPZPf7G7N9LNfWlf+Y1Sat0w929WBqEZ+ZJHXjwd+Xcfr3y8y7RtgybznqxDJkPeKLDumyLQ5GefuUwumTUg/lyxcuJE2IbLhjwN7uPuP9SxbuO25BMwHRZb9tgSxiUgTpJPO62xm7ahP0vPZ/ifNbGUzu9/MJqdHfzNbxcw+NLMhDXzLPYER7v5VwfRckr1YLa5cd9JOReblvzZXY6qXmb1uZhNTP/I3zexUy6unZWZ/SF2IjikWpJm9lPqvt8mb1sHM7kjTp1rUiLrUzBZsyIabWdu0/Hs2s9//nYXNsM3ssBTbVmZ2mpm9bzPrgu1fZL2tzexci/oAP6Zt39MKar2Y2bPERUdrm7VGwH4F61vczPqm+H60qIGyceH7pmTWIGDnhu6DhjKzdmb2YvobdirYL51sZj2E/1nUJdskLdPZzIZZ1Pb6zMzOLGVcIlJRct3wc8mCdsT55PYUP1/fjuLn66OLTHsr/Vy1yLy6FDvvL9U5+B+Bi4Eb3f0od5/eiDh0Dl6jNCpQdRtM9OfrRBzkcvVVcoYCV5jZwmmZGUSyJd+/ia4mfYjmwROIrjZ/IOqi/JJ8c/efzexu4FgzW8XdPyCSLBOIpEtOrv/3yUBhLZCcT+uYnv/626i7HssPRabVddBrzv7o9R1o5/Z930k/tyW+nAY1IY4s9omIFJHXF3014g7YK8QIbkcDXcxsE3f/Li27JHHXaxni+Dwa+D1x3F9o9rUXfb9liER1seHi8+t03ZCWz7VIOZfoRvh/Zra6u7+bt2z+awE2IIbifJBIWM8L7ESMILEyUR8L4DGiFeUBRL/6/DjXBDoCvdx9Wpq2CTGq3TdE7arPgfWJ4Tc3N7POuWXr2PbFiX29PLGv3wJ+RdQB2NbMfuvunxS87FKiVWZvoivUMUA/M3vH3V/IW64PcBjwNHA5sDTwT2Y/uT6f2JebAQfmTR+WHyrRIvEz4DziQuVk4BEzW9XdpxSsczhwKPA7Zh31r8ksCkwOIrrxbpXuqua7LMV5JbF/TgGeMLODic9OX+K7em/gb2b2vrvfXYrYRKSi5GoR5m5E5s4lnyLqppRbc56Dv0h8h+1hZn096lwVVU/SRefgNUaJleo2HPiRuOM4mWhx8lLe/KHE37gT0ZrltfyuJDazy83N7j7LXUQrMppOchtwHLC/mfUjTvSvK2i1kTsJn+LuTTnxG0s0XZ+nia+vz4dEcdlfM3vXnw5Flvci08plIvAn4oS3v5nt4e4DM4xHRObOaUQT4GPd/Zfkgpm9RnQpPI0oIA7wf0R9lP3c/Y40rbeZXUrUqmqItdLP2VrouftYi5HPOuVNzv0+hGjGPZ34fslPrHwL5F94Pw38uqAr4pVmdhdwpJn1dPev3H2amd0JnGRmv3H3d/KWz7WCvA1+aV5+M9HtZZP85EJqqXMPcRFfXyH0C4k+/pu4e66lJGZ2G5HwP5dIjuRrA2zs7j+nZR8g9t1xRNdbzGy99LpHgT/kukSZ2f3Ay/krc/cnzOwAYFMvGL40TytguLufkBfjGKII417MLBaZk/tbrk0JEisWRYgfJW6QbO6zF2XM2Txvv4wB7if+Dpu4+6tp+s3AJ0QyTYkVkZYnN6pk7mbqeOJcddFGnq+vSQwakS/3fZWfoM7yHHwckRD/D/CUmXV19+czjEeqgLoCVTF3/4lIrvyWqEA9vCDBMYo4WTqVuMM5pGAVuUzqLJlTM1seOIQiUsb2LWZW5zbyRgNKcidpZ6S7hrMwswXMbJF6tutL4g7ensWaQ1toV9fr5+Dh9PPkgnXuSvE+j1OYtTtUWaV+p9sTJ+z3m9luWcUiInNtd+JEtG/B9H+m6bvnTfsD0UqjsDbK5Y14v9xx8ps65g8GfmVmuVpanYB33P2L1HLmlTQtvzXL0Pz6Ku7+v1xSxczmNbMlLPrXDyIKoeePHpT7rvilO6nF6Gn7Eon/XMJmAyJxcAcwv83ab38ocUOhvjpdrYC/EN95XxS8/jviTmSx11+XSx6kbfuY2WtS5ep0XVmwH14jkkxNcUXB8/pqBOSauS/dxPf6hcVIF0OIbSw20kXO9fn7hZn1A4blkirwyznJS6h+gEiLkro/Xg5sCTzq7sPgl+6JdwCbmNkedby22LHq6NSCM7fMYkTtqInM2vJ+CrBEabai8VKNwq2JFoVPWAmGu5eWTS1Wqt9g4o7i74g7cL9wdzezZ4hm2rll8+dPNLOngQPN7Cfi4n0l4uD2HtE0u5jbiCZ/pwBvFTaPc/cp6S7dA8AYM7uFOHFbHFiD6Le4C9H1qC5HpvnPppYxrxKf11XS9txE3JFsrAHEye9RKTnzH6L1ymHEXcw1C5Z/HuhqZtek36cTo/V8TSOlu8PLunuj/u/cfXI6AX4EuMfM/uLu9zX2/UUkc6sQ9U5m6cKSWnO8Q4y4kL/si4VFYt39KzOb2MD3y93tq6vZ8WAiQd6Z6HrYidm7k+6bfu+Ufg7JX4GZzUOMhLY/cSwtfK9fEtPu/pqZjQT2M7NzUkKmM7ACsyYXcsfhC6n7OF9fna5lgbZEi8zxdSxTWBsL6u6vn/9eufosxWpyjSG6bTbGDGI0jML3hOL98HP7d27v5C5PDN88Eti2nqKMoPoBIrVko7w6UIsQrbl3I64PniCS1vnOIlrF32Nm9xDnylPT8jsR1xYHFbzma2Io5lvS84OBFYmR6PK7+j8PHGpmFxDdYWcAD+cPvNFQZubAR+6+cmNe5+5fpLpTTwGPm9ku7j50Di+TGqXESvXLT5YU+0cfShwQpzPzLlO+fYgkyW7Ege9dYqx1o+7Eyu1EQadFqWOseHd/NLU2OZ044V6KOOEaS/TXfrOebcLdPzKzjdLrdyXucP5ANDN+CGhSYiElm3YDLiKaWe8MvEaMbHQS0Z8y3z/StL2I5s1GdH+qLyk0m9S0fSHqry1TX9xTzGxHIjF0l5nN4+6Fd7JFRPLlkgp13fH7pc6KmQ0gWqScnzd/KHCKmXWgeH0ViKE2jyZa1lyQ3vNnYiS6i5i9ZWw/otXN1kSS5gCioPgdecvkkgeXAk/WEfuEOqbnv34QdbfwKTaqURb94b2eEZaKvW/ub1lXwqihcsXkdyC6Vd1Sz7KqHyBSO/ZJjxlEi5FxxHfBXe7+eOHC7j4pteToTowa2o04po8jzpVvLPIe/0ecSx9LJK7fAfZ19zsLljuLOOYdSyTLjUhuNyqxktdKvqnn4F+ZWWciufKome3q7k1toSgtmBIrVc7dn6Wekxh3v5LZhwHLnz+eOrr9UPxgiLt/RjTxnlNsbzDzbmd9y+0H7Fdk+njiQN19Dq8/Gzi7jnnti0ybQhTmPSF/upn1JsacL1z24DrWvWUd08cy+99kQ+JL4YTZX9HguH8gCtnmTyu67XXEUG/cItLs3gc6mFmb/FYrFiPh/IZZWwZ8CKxmZq3yL7xTs+q2DXy/XAK7aNcMd//QzD4kkiad0uT8BP2zxMl15zR/ArMXJN8P+I+7z3IX08zWqCOmO4hk/gFm9hLRgvHxglGLcjVdpjWxztaXRJefRZqpThfEXdzCrjPlqNOVGx57VL1LzdlU4obKvcBNKVlf2EVNRGqEuw+hiUnRdH56QXo0xDR3P5eClvZF1vsVUWuw2LyV65g+hNm3Izf6XdFrhSLrmG0/uPsE4lw+f9pBzN4ap9592dgWM1JdVGNFao6ZLVBkWjeim1Jdd0fn1g5EvYL6ii2KSMv2EFH3pLBo6uFp+oN50x4GliPuHOY7paFvlpLTbxKj0tRlMNF15ijgfXcfl/f6icDrKYZcfZXCRMF0Zq/TtQjRArBYTF8QLUn2IBLvCzN7na4RRLPvY8xs5cJ1mNk8xep35b3HNKL46+/qqktVR7//hsjV6Top1XLJrW8DojB6oSnEcMuLFpnXFJsRLYKGz+2KUk22PYjPXR8zO3YOLxERqUY7AAPdvbDFpUhJqcWK1KKeZrYucUExmahrcDDRNPqy5nhDd7+Y6D4lIi3XNmY2f5HpX7t7H6Jry57Adamr46vEHbBDifocl+a95hKiL/stFkMPv000nf4d0T+9oS0h7gXOMbPl3P3zIvMHE8e/rSjeHWQoM5MkxU5K7yf6wN9F1KxalmgFWV8dqtuIvveXEV1EH86f6e4zzGx/otn1G2m0mbeI7pSrEXcwu1N/ovp0Yl/dn9fv/2eia+dO6XlhgmuO3P31FM8hRDHDh4im7McSf8/fMuvf5nkiadXHzB5LMQx398K6KnOUEjldgUcK6hBgZs8SdQ5WyE+ONWB7fjazzt/LoAAAAi9JREFUvYh9eW1quVJnK1cRkWrj7koaS1kosSK16L/A5sTQposRI2bcA5yTujmJiDRF1/QoNAbok9cXvSdRO+pgottKH+DcNBIPAO7+tZltSdR5OoS4WM8VK38JqK/YaL4biObPf0nrKtSQOl25xMqQIvNPBCYRCaPdie4xvYnaVYPqiGkAMfpDW2K//FS4gLu/bGYbAmcSXVaOJhLhHxLdVIvFkv/6iWa2OdHCZ8+0jly///9SR1fXBjoyrecQYp++DRxBjJjxW2b929wOrE/U6dqLaCm8P7MXrG2IzkTR2UOLzFuEaB0zqbErTcWT9yX2zxUpudIsNxlERERaKpu9Va+IiIhUIjNbkmgN8k93P6qBr+lDDC/coWDYXCmh1CJlC2CxIl2mSrH+h4lhljfLX38aRvpL4Dx3b2iNAxERESkh1VgRERGpQMXqQRFdXKBx9aB6EEPgFi3ELY1TR52uDYnk1dPNlFTpSIxi173I+rcDPqd4iyQREREpA7VYERERqUBmNpjoMvIKcSNkG2AX4DlgK3eva7hbaUapyOs+wGNEba41ia5AAJu7+8isYhMREZFsKLEiIiJSgcysO3AAUXB1AaKuxwNAz/x6LFJeqXZLT6J2yuLE0M7PEl1xXs0yNhEREcmGEisiIiIiIiIiIk2kGisiIiIiIiIiIk2kxIqIiIiIiIiISBMpsSIiIiIiIiIi0kRKrIiIiIiIiIiINJESKyIiIiIiIiIiTfT/4IlJhhakhtcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1152x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(16,6))\n",
    "\n",
    "ax = fig.add_subplot(1, 3, 1)\n",
    "ax.plot(power_spectrum.sigmakm.values, power_spectrum.values, c='blue')\n",
    "ax.set_yscale('log')\n",
    "ax.set_title(f'Power Spectrum (Log)\\n', fontsize=22)\n",
    "ax.set_ylabel('Log (Spectral density)',fontsize=18)\n",
    "ax.set_xlabel('Wavelength, km', fontsize=18)\n",
    "ax.set_yticks([])\n",
    "ax.set_yticks([], minor=True)\n",
    "\n",
    "ax = fig.add_subplot(1, 3, 2)\n",
    "ax.plot(power_spectrum.sigmakm.values, power_spectrum.values, c='blue')\n",
    "ax.set_xscale('log')\n",
    "ax.set_yscale('log')\n",
    "ax.set_title(f'Power Spectrum (LogLog)\\n', fontsize=22)\n",
    "ax.set_ylabel('Log (Spectral density)',fontsize=18)\n",
    "ax.set_xlabel('Log (Wavelength), km', fontsize=18)\n",
    "ax.set_yticks([])\n",
    "ax.set_yticks([], minor=True)\n",
    "\n",
    "ax = fig.add_subplot(1, 3, 3)\n",
    "# calculate fractality index\n",
    "slope = np.diff(np.log10(power_spectrum))/np.diff(np.log10(power_spectrum.sigmakm))\n",
    "depths = (power_spectrum.sigmakm.values[1:]+power_spectrum.sigmakm.values[:-1])/2/np.sqrt(2)\n",
    "fractal = (3 - (slope/2))\n",
    "ax.plot(depths, 1000*fractal, c='blue')\n",
    "ax.set_title(f'Density\\n', fontsize=22)\n",
    "ax.set_ylabel('ρ, kg/m³',fontsize=18)\n",
    "ax.set_xlabel('Depth, km', fontsize=18)\n",
    "\n",
    "plt.suptitle('WGM2012 Free Air Gravity Power Spectrum and Fractality Density', fontsize=28)\n",
    "fig.tight_layout(rect=[0.03, 0.0, 1, 0.9])\n",
    "plt.savefig('Spectral Components Analysis [WGM2012 Free Air].jpg', dpi=150, quality=95)\n",
    "plt.show()"
   ]
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}