Rayleigh_drag.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rayleigh drag\n",
    "Display closeup of Rayleigh drag (if any).\n",
    "\n",
    "Based on `/home/156/aek156/notebooks/github/aekiss/notebooks/bathymetry-closeup-KaraStr-Rayleighdamping.ipynb`\n",
    "https://github.com/aekiss/notebooks/blob/master/bathymetry-closeup-KaraStr-Rayleighdamping.ipynb\n",
    "\n",
    "and\n",
    "https://github.com/COSIMA/ACCESS-OM2-1-025-010deg-report/blob/master/figures/Rayleigh_drag/Rayleigh_drag.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/g/data3/hh5/public/apps/miniconda3/envs/analysis3-19.01/lib/python3.6/site-packages/dask/config.py:168: YAMLLoadWarning: calling yaml.load() without Loader=... is deprecated, as the default Loader is unsafe. Please read https://msg.pyyaml.org/load for full details.\n",
      "  data = yaml.load(f.read()) or {}\n",
      "/g/data3/hh5/public/apps/miniconda3/envs/analysis3-19.01/lib/python3.6/site-packages/distributed/config.py:20: YAMLLoadWarning: calling yaml.load() without Loader=... is deprecated, as the default Loader is unsafe. Please read https://msg.pyyaml.org/load for full details.\n",
      "  defaults = yaml.load(f)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "import re\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "figdir = ''\n",
    "NOSYNC = 'NOSYNC/' # where to save figs we don't want to sync\n",
    "if not os.path.exists(NOSYNC):\n",
    "    os.makedirs(NOSYNC)\n",
    "def savefigure(fname):\n",
    "    plt.savefig(os.path.join(figdir, fname+'.png'),dpi=200, bbox_inches=\"tight\")  # comment out to disable saving\n",
    "    #plt.savefig(os.path.join(figdir, fname+'.pdf'),dpi=200, bbox_inches=\"tight\")  # comment out to disable saving\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = '/g/data/ik11/inputs/access-om2/input_rc/mom_1deg/'\n",
    "path = '/scratch/v45/aek156/access-om2/archive/1deg_jra55_ryf_v2.0.0rc7/output000/ocean/'\n",
    "ftfile = path + 'field_table'\n",
    "# gridfile = xr.open_dataset(path + 'ocean_hgrid.nc')\n",
    "# gridfile = xr.open_dataset(path + 'ocean_grid.nc')\n",
    "# geolon_c = geolon_c.x\n",
    "# geolat_c = geolat_c.y\n",
    "# kmu = gridfile.kmu\n",
    "# hu = gridfile.hu  # hu has correct u-cell mask (kmu uses t-cell mask)\n",
    "\n",
    "geolon_c = xr.open_dataset(path + 'ocean-2d-geolon_c.nc').geolon_c\n",
    "geolat_c = xr.open_dataset(path + 'ocean-2d-geolat_c.nc').geolat_c\n",
    "kmu = xr.open_dataset(path + 'ocean-2d-kmu.nc').kmu\n",
    "hu = xr.open_dataset(path + 'ocean-2d-hu.nc').hu  # hu has correct u-cell mask (kmu uses t-cell mask)\n",
    "\n",
    "damp_table = ''\n",
    "with open(ftfile,\"r\") as f:\n",
    "    for line in f:\n",
    "        if line.startswith('\"rayleigh_damp_table\"'):\n",
    "            break\n",
    "    for line in f:\n",
    "        if line.startswith('\"rayleigh\"'):\n",
    "            damp_table += line\n",
    "        else:\n",
    "            break\n",
    "\n",
    "damp_table_ii = [ int(i)-1 for i in re.findall('itable=(\\d+)', damp_table) ]  # subtract 1 to index from 0\n",
    "damp_table_jj = [ int(j)-1 for j in re.findall('jtable=(\\d+)', damp_table) ]  # subtract 1 to index from 0\n",
    "damp_time = [ int(j) for j in re.findall('rayleigh_damp_table=(\\d+)', damp_table) ]\n",
    "ktable_1 = [ int(j) for j in re.findall('ktable_1=(\\d+)', damp_table) ]\n",
    "ktable_2 = [ int(j) for j in re.findall('ktable_2=(\\d+)', damp_table) ]\n",
    "\n",
    "# make damp time map; use kmu to get land mask\n",
    "dtm = hu.where(np.isnan(hu), 0)  # zero out ocean points\n",
    "for i, j, dt in zip(damp_table_ii, damp_table_jj, damp_time):\n",
    "    if not np.isnan(dtm[j][i]):\n",
    "        dtm[j][i] = dt  # put damp_time value in at every ocean point that has damping\n",
    "damp_time_map = dtm\n",
    "\n",
    "# make map of number of non-damped cells in each damped column\n",
    "undamped = hu.where(np.isnan(hu), 0)  # zero out ocean points\n",
    "for i, j, kt1, kt2 in zip(damp_table_ii, damp_table_jj, ktable_1, ktable_2):\n",
    "    if not np.isnan(undamped[j][i]):\n",
    "        undamped[j][i] = max(0, kmu[j][i] - kt2) + kt1 - 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAicAAAFNCAYAAADb1898AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XvclHWd//HXW1Q8hwq6CCLaUnlIMJHsoOEhJXPTNBV/lpTukmlla7XqtqW5y/7sZ0crD6yauJVKnk1RkVSq9QQeQXJFRSUIBEQpFQU+vz+u743Dzczccx+ue66ZeT8fj3nMzPc6zOe67jl87u/1PSgiMDMzMyuKDeodgJmZmVkpJydmZmZWKE5OzMzMrFCcnJiZmVmhODkxMzOzQnFyYmZmZoXi5MSsRpJC0t/XsN5+kp6ucZ+jJc3vRkznSvplV7fvxuvWfIw99Hp/lbRLL73WIZJuqmG9H0o6pTdiMms1Tk6sqUiaJ+mN9GP2F0lXStqiN2OIiN9HxHt78zV7W57HKOleSf/Y7vW2iIjn8ni9Mv4TOL+G9S4AviVp45zjMWs5Tk6sGf1DRGwBjAD2As6uczzWICTtA7wrIh7oaN2IWAj8CfhU7oGZtRgnJ9a0IuIvwJ1kSQoAkj4p6VFJr0l6SdK5Jctuk/SV0n1IekLSke33LamvpO9LelHSIkmXSNo0LVvnUo2kD6TXXCHpN5KulfQf7fb3dUmLJS2U9IVKxyRpZ0n3pX1NBfq3W/6bVGP0qqTpknYvWXalpIskTUk1S3+U9HeSfizpFUl/krRXyfrzJJ0t6am0/BeSNqlwjPMkfSOdr1fTMW5Ssvxf0rEtkPSPlS6RSZoA7Af8LMX4s1S+dv0uHMcOkq6X9LKk5yV9tdL5BT4B3FeyrST9KP1tXk3Ht0fJ+vcCn6yyPzPrAicn1rQkDSb7sZlbUvw34ESgH9mPypdKko9JwGdLth8ODAJuL7P77wHvIUt8/j6t950yMWwM3AhcCWwDXA18ut1qfwe8K+3jZODnkraucFi/BmaSJSX/Doxrt3wKMAzYDngE+FW75ccC/5a2Xwncn9brD1wH/LDd+icAhwLvTsf7bxXiatv3GGBnYE/g8wCSxgBnAAeTnauPVdpBRHwL+D3w5XQp58tVXqvD45C0AXAr8DjZ+T0I+JqkQyvs9/1AaVuaQ4D9yY69H3AcsLRk+RxgeKXjMbOucXJizegmSSuAl4DFwDltCyLi3oh4MiLWRMQTZMlC24/lzcAwScPS888B10bEW6U7lyTgn4B/johlEbGCrJ3C2DKx7AtsCFwYEW9HxA3AQ+3WeRs4Ly2/HfgrsF57DklDgH2Ab0fEyoiYTvbDu1ZEXBERKyJiJXAuMFzSu0pWuTEiZkbEm2RJ05sRcVVErAauJbsMVupnEfFSRCwDJgDHlznGNhdGxIK07q28U2N1LPCLiJgdEa8D362yj1rVehz7AAMi4ryIeCu1W/kvyv+tIEtAVpQ8fxvYEngfoIiYky7ntFmRtjGzHuTkxJrRkRGxJTCa7Edl7aUPSR+UdE+q4n8VOKVtefpBnwx8Nv3HfTzw32X2PwDYDJgpabmk5cAdqby9HYA/x7ozbL7Ubp2lEbGq5PnrQLlGvDsAr0TE30rKXig5tj6Szpf0rKTXgHlpUemln0Ulj98o87z965bG+kKKoZK/VDiGHdrtp/3xd0Wtx7ETsEPb3yn9rf4V2L7Cfl8hS0YAiIjfAT8Dfg4skjRR0lYl628JLO/WkZjZepycWNOKiPvILqd8v6T418AtwI4R8S7gEkAlyyeRXco4CHg9Iu4vs+slZD+Au0dEv3R7V2qE295CYFCqbWmzYxcPaSGwtaTNS8qGlDz+P8ARZJdP3gUMTeWlr91ZpbEOARZ0YR8LgcEV9llOT06V/hLwfMnfqV9EbBkRh1VY/wmySzjvBBNxYUTsDeyeln2zZPGuZJeMzKwHOTmxZvdj4OOS2i4xbAksi4g3JY0i+0FfKyUja4AfUL7WhIhYQ3Zp4EeStgOQNKhCO4b7gdXAlyVtKOkIYFRXDiQiXgBmAN+VtLGkjwL/ULLKlmTtL5aS1ez8Z1dep53TJA2WtA1ZjcO1XdjHZOALknaVtBll2ua0swjoqTFNHgJek3SmpE1T7dIeynrllHM7JW1iJO2Tats2Imuv9CbZ37PNx8ja+ZhZD3JyYk0tIl4GrgK+nYpOBc5LbVK+Q/bD2d5VZA0jqw1udiZZQ9sH0iWUuynTTiS1VzmKrKHrcrIGt78lSyK64v8AHwSWkbWluapd3C8AfwaeAjrsDluDXwN3Ac+l239UX319ETEFuBC4h+yctdVGVToHPwE+k3reXNjpiNd97dVkCdwI4HmyWq/LyGqWyq3/CPCqpA+moq3IEtFXyM7tUlJNnKSBwG5AhwO2mVnnaN1L4WYm6URgfER8NKf9PwhcEhG/yGP/PUXSPOAfI+LuHt7vrsAsoG+7tjaFIOkQ4NSIWK8Lebv1fgA8GxEX9U5kZq1jw3oHYFYk6bLDqUCP/eBI+hhZ99QlZO1Z9iRrQNsyJH0auA3YnKwb9q1FTEwAIuIustqijtb7ei+EY9aSfFnHLEltRl4ma/Pw6x7c9XvJGk2+Cnwd+Ey77qit4Itk5/ZZsjYbX6pvOGZWZL6sY2ZmZoXimhMzMzMrFCcnZmZmVigN3SBW/fsHQ4eufb53/UKpama9AyijGc9VvY4pz79vUf9O1XR0PhrxmIqqiN8t3ZXn+6PXztfMmUsiotyI0bkYozGxhCWd3m4mM++MiDE5hNRtDd3mRCNHBjNmrH1e1CPpzvCceWnGc1WvY8rz71vUv1M1HZ2PRjymoirid0t35fn+6LXzJc2MiJG99XIjNTJmMKPjFdsRvRtnZzR0zYmZmZlBdCXzKvB/Ck5OzMzMGlyzJSe5NYiVdIWkxZJmlZRdIOlPkp6QdKOkfql8qKQ3JD2WbpfkFZdllOPNzHpHdHBrRP5u6ZpQ529FlmdvnSuB9g1tpgJ7RMSewP8CZ5csezYiRqTbKTnGZWZm1jQCJyc1i4jpZJOTlZbdVTJk9QOsO426mZmZdVYXEpOWTU5qcBLrTjW+s6RHJd0nab96BWVmZtZomi05qUuDWEnfAlYBv0pFC4EhEbFU0t7ATZJ2j4jXymw7HhgPwJAhvRSxmZlZcRU92eisXq85kTQOOBw4IdIgKxGxMiKWpsczySYHe0+57SNiYkSMjIiRDOi1MW7MzMwKyzUn3SBpDHAm8LGIeL2kfACwLCJWS9oFGAY815uxmZmZNaK2BrHNJLfkRNLVwGigv6T5wDlkvXP6AlMlATyQeubsD5wnaRXZdOqnRMSysju2wmuyz0iHOuqy2Z3z0Wrn0npWtfem31tNpAFqQjort+QkIo4vU3x5hXWvB67PKxYzM7Nm5uTEzMzMCsXJiZmZmRWKkxMzMzMrjGZsEFvPQdjMzMysu3IaIVbSPElPpjnvZqSybSRNlfRMut+6ZP2zJc2V9LSkQ0vK9077mSvpQqUeMdU4OTEzM2twOY5zckCa825ken4WMC0ihgHT0nMk7QaMBXYnm1fvIkl90jYXkw2eOizd2s+7tx4nJ2ZmZlarI4BJ6fEk4MiS8mvSoKrPA3OBUZIGAltFxP1p4NWrSrapqKmSE0+1bV1V7b3j90/P6uhc++/Q2qLKzSrLqeYkgLskzUxTxwBsHxELAdL9dql8EPBSybbzU9mg9Lh9eVVuEGtmZtbAutEgtn9bW5JkYkRMLHn+kYhYIGk7ssFT/1RlX+UiiCrlVTk5MTMza3BdTE6WlLQlWX+fEQvS/WJJNwKjgEWSBkbEwnTJZnFafT6wY8nmg4EFqXxwmfKqmuqyjpmZWcvJobeOpM0lbdn2GDgEmAXcAoxLq40Dbk6PbwHGSuoraWeyhq8PpUs/KyTtm3rpnFiyTUWuOTEzM2twOYxzsj1wY+r1uyHw64i4Q9LDwGRJJwMvAscARMRsSZOBp4BVwGkRsTrt60vAlcCmwJR0q8rJiZmZWYPr6eQkIp4DhpcpXwocVGGbCcCEMuUzgD068/pOTszMzBpYM44Q6+TEzBpGR9+/7m5au47OVSP+1jVizD3FyYmZmZkVR+dGfG0ITk7MzMwanJMTMzMzKxQnJ2ZmZlYYbhBrZmZmhePkxMzMzIrDDWLNzMysaJycmJkVVLXvZ4+BYs2s2ZITT/xnZmZmheKaEzMzswbm3jpmZmZWOE5OzMzMrDjcW8fMzMyKxsmJmZmZFUqzJSe59daRdIWkxZJmlZRdIOlPkp6QdKOkfiXLzpY0V9LTkg7NKy4zM7Nm0tYgtrO3IsuzK/GVwJh2ZVOBPSJiT+B/gbMBJO0GjAV2T9tcJKlPTwajDm5m1tz8HWDNzMlJjSJiOrCsXdldEbEqPX0AGJweHwFcExErI+J5YC4wKq/YzMzMmkYXEpOWTU5qcBIwJT0eBLxUsmx+KjMzM7MONFtyUpcGsZK+BawCftVWVGa1sqNNSxoPjAdgyJA8wjMzM2soRU82OqvXa04kjQMOB06IiLYEZD6wY8lqg4EF5baPiIkRMTIiRjJgQL7BmpmZFZwbxHaTpDHAmcCnIuL1kkW3AGMl9ZW0MzAMeKg3YzMzM2tUzZac5HZZR9LVwGigv6T5wDlkvXP6AlMlATwQEadExGxJk4GnyC73nBYRq/OKzczMrGk0QLLRWbklJxFxfJniy6usPwGYkFc8HclzqvVq2zfZ+6lp+e/UO7r7WfPfqed053uru39H6zwnJ2ZmZlYoTk7MzMysMNoaxDaTeo5zYmZmZrYe15yYmZk1uGarOXFyYmZm1sjcW8fMzMyKxsmJmZmZFYqTE7OCarLPppm1k+dnvDtjs9T7u6cZe+s4OTEzM2twTk7MzMysONwg1szMzIrGyYmZmZkVipMTMzMzKww3iDUzM7PCcXLSgvKcHryjbZvs/WZmLcDfW73MDWLNzMysaJycmJmZWaE0W3KyQb0DMDMzs65raxDb2VstJPWR9Kik36bn20iaKumZdL91ybpnS5or6WlJh5aU7y3pybTsQkkdvrqTEzMzswaXV3ICnA7MKXl+FjAtIoYB09JzJO0GjAV2B8YAF0nqk7a5GBgPDEu3MR29qJMTMzMzW4+kwcAngctKio8AJqXHk4AjS8qviYiVEfE8MBcYJWkgsFVE3B8RAVxVsk1FbnNiZmbWyLreW6e/pBklzydGxMSS5z8G/gXYsqRs+4hYCBARCyVtl8oHAQ+UrDc/lb2dHrcvr8rJiZmZWYPrYnKyJCJGllsg6XBgcUTMlDS6hn2ViyCqlFfl5KTgqv0F8xx/pckafluDqOf7rtprd+ezZNYbcuit8xHgU5IOAzYBtpL0S2CRpIGp1mQgsDitPx/YsWT7wcCCVD64THlVbnNiZmbWwPLorRMRZ0fE4IgYStbQ9XcR8VngFmBcWm0ccHN6fAswVlJfSTuTNXx9KF0CWiFp39RL58SSbSpyzYmZmVmD68VxTs4HJks6GXgROAYgImZLmgw8BawCTouI1WmbLwFXApsCU9KtKicnZmZmjSzn4esj4l7g3vR4KXBQhfUmABPKlM8A9ujMazo5MTMza3DNNkKskxMzM7MG5+TEzMzMCqOtQWwzya23jqQrJC2WNKuk7BhJsyWtkTSypHyopDckPZZul+QVVzOJDm557dvMrNmoG7ciyHH4+rrIs+bkSuBnZEPVtpkFHAVcWmb9ZyNiRI7xmJmZNZ8GSDY6K7fkJCKmSxrarmwOQA0TEpqZmVmNmi05KdIgbDunaZnvk7RfvYMxMzNrFL6sk4+FwJCIWCppb+AmSbtHxGvtV5Q0nmzqZRgypHejNDMzKxg3iM1JmmJ5aXo8E3gWeE+FdSdGxMiIGMmAAb0ZppmZWSG55iQHkgYAyyJitaRdyMbkf67OYZmZmRVfAyQbnZVbciLpamA00F/SfOAcYBnwU2AAcJukxyLiUGB/4DxJq4DVwCkRsSyv2MzMzJqJk5MaRcTxFRbdWGbd64Hr84rFzHpWd8a7Kep3qMfwMSuOQlzWMTMzs65zzYmZmZkVRjP21nFyYmZm1sjcINbMzMyKpmWTE0mbA29GxOoc4zEzM7NOapnkRNIGwFjgBGAfYCXQV9LLwO3AxIh4pleiNDMzs4paJjkB7gHuBs4GZkXEGgBJ2wAHAOdLujEifpl/mMXWnfdEI3Zf7CjmJvuM1FUzdtm1xufvgGJptQaxB0fE2+0L0+Bo1wPXS9oot8jMzMysY63UILYtMZG0PTCILDlbEBGL2q9jZmZm9dMyyYmkvYCLgXcBf07FgyUtB06NiEd6IT4zMzPrQMskJ8AvgC9GxIOlhZL2TcuG5xmYmZmZ1aaVkpPN2ycmABHxQOpWbGZmZnXWag1ip0i6DbgKeCmV7QicCNyRd2BmZmZWgxZrEPtVSZ8AjiBrECtgPvDziLi9l+IzMzOzDrRMcgIQEVOAKb0US0tqxjFSqsXVZJ8fM7P1vvPq8T3XUsmJmZmZFVsztjnZoN4BmJmZmZVyzYmZmVmDa7aakw6TE0k7A18BhpauHxGfyi8sMzMzq0kr9dYpcRNwOXArsCbfcMzMzKyzWjE5eTMiLsw9EjMzM+uSVkxOfiLpHOAuYGVboefWsW5bsgTmzIGhQ2HHHesdjZlZQ2rG3jq1JCfvBz4HHMg7l3UiPbc66ui9WMRxUEKwhjWcyZlcxmXsxm48zdMcwiFczuVsGpvWO8Qe153vjKJ+3xQ1rlbTjOMkWde0YnLyaWCXiHgr72CsNfycn/NH/shc5rIt2/I6r/N5Ps83+SY/42f1Ds/MrLE0YYPYWsY5eRzol3cg1jou5VIu4AK2ZVsANmMzfsJP+G/+m7fecg5sZtZZoc7fiqyWmpPtgT9Jeph125y4K7F1yRKWsBM7rVO2PduzilW8+eabbLzxxnWKzMysMRU92eisWpKTc3KPwlrKaEZzNVfzTb65tuxWbuW9vJctt9yyjpGZmTWelmwQGxH3Sdoe2CcVPRQRi/MNy5rZuZzLaEbzMi/zcT7ODGbwI37E1VyN1GSfMDOzXtBsyUmHbU4kHQs8BBwDHAs8KOkzeQdmzet9vI8HeZA1rOF7fI95zGMa0ziIg+odmplZ4+lCe5OiJzO1XNb5FrBPW22JpAHA3cB11TaSdAVwOLA4IvZIZccA5wK7AqMiYkbJ+mcDJwOrga9GxJ2dPhpbR17vve50QdTajXcCvg/ANGBiD+y7qNq6T1/HddzADWzABhzDMRzJkSjHb4iCf/cUjrvlWlcV4bPW018lkjYBpgN9yXKF6yLiHEnbANeSTWkzDzg2Il5J25T9HZe0N3AlsClwO3B6RFT92NTSW2eDdpdxlta43ZXAmHZls4CjyA54LUm7AWOB3dM2F0nqU8NrmBVeEIxjHBdwAWMYw4EcyDmcw2mcVu/QzKxJ5FBzshI4MCKGAyOAMZL2Bc4CpkXEMLL/Lc+CDn/HLwbGA8PSrX1usJ5aak7ukHQncHV6fhxZ5lNVREyXNLRd2Zx0EO1XPwK4JiJWAs9LmguMAu6vIT6zQvsDf+ABHuAJnmBTskHmxjKW9/E+xj82nhEjRtQ5QjNrZHk0iE01G39NTzdKtyD7vR6dyicB9wJnUuF3XNI8YKuIuB9A0lXAkcCUaq/fYQ1IRHwTuBTYExgOTIyIM2s+wtoMAl4qeT4/lZk1vLu5m2M5dm1iArAFW3AUR3H33XfXMTIzaxZ5tDmR1EfSY8BiYGpEPAhsHxELAdL9dmn1Sr/jg9Lj9uVVVa05SVUyd0bEwcANHR9Kl5U7TWWvR0kaT1Y9BEOG5BiSWc/oRz9mMWu98oUsZPjWw+sQkZkZAP0lzSh5PjEi2poAEhGrgRGS+gE3Stqjyr4q/Y7X/PteqmrNSQrsdUnv6mhH3TQfKJ35bTCwoEJMEyNiZESMZMCAnMMy677jOZ6buIk/8se1ZXdzN/dyL0cffXQdIzOzptD13jpL2n5P021iud1HxHKyyzdjgEWSBgKk+7Y2qZV+x+enx+3Lq6qlzcmbwJOSpgJ/Kwn2qzVsW6tbgF9L+iGwA1mDmYd6cP9mdfN3/B2/4lcczdHszM6sYhULWMB1XEe/fp4Zwsy6L4feOgOAtyNiuaRNgYOB75H9Xo8Dzk/3N6dNyv6OR8RqSStSY9oHgROBn3b0+rUkJ7elW6dIupqs0Ux/SfPJRppdloIaANwm6bGIODQiZkuaDDwFrAJOS7U2Zk1hDGN4gRf4I39kAzbgI3yEjdio3mGZWZPIYVSCgcCk1LxjA2ByRPxW0v3AZEknAy+SjYFGB7/jX+KdrsRT6KAxLFRJTiRNi4iDgN260gA2Io6vsOjGCutPACZ09nWs9xWhT39X1C3utVdX+wIHll9kDa3ae6uof+OOPg/dibvato36/VFkOfXWeQLYq0z5Uig/Ymal3/E0plm19irrqVZzMlDSx4BPSbqGdu+piHikMy9kZmZm+Sj6iK+dVS05+Q7Z4CqDgR+2Wxa0/xfQzMzMel8DDEffWRWTk4i4DrhO0rcj4t97MSYzMzPrhJZJTto4MTEzMyu2lktOzMzMrLjyaBBbb05OzMzMGlzLJCdpWuSKImJZz4djra7JPl9Avl1Jm/F8NZsOu+wK1rCGu7iL27mdzdiMEziB9/N+VNR+yFYsTdggttrw9TOBGen+ZeB/gWfS45n5h2Zm1vzWsIYTOZFv8A12ZEc2ZEM+zse5lEvrHZo1kDwm/qunar11dgaQdAlwS0Tcnp5/gmwYWzMz66YpTOFJnmQGM9iETQA4iZPYm71h6Wdg223rHKE1gqInG51VdeK/ZJ+2xAQgIqYAH8svJDOz1nErt3ISJ61NTAB2YRdGMxqmTq1fYNYw2hrEtkTNSYklkv4N+CXZOfgssDTXqMzMWsQmbMIKVqxXvoIVsOmmdYjIGlHRk43OqqXm5HiyifpuTLcBqczMzLrpBE7gEi5hPvPXlv2O3/EkT8Ihh9QxMmsYXag1KXoyU8sgbMuA0yVtERF/7YWYzMxaxj7swxmcwft5P4dyKMtZzkxm8ht+wwGuObEW1WHNiaQPS3qKbBpkJA2XdFHukZmZtYgzOIPZzOZgDuYLfIF5zMvanJjVqOVqToAfAYcCtwBExOOS9s81KrMuaMbxRAr+/WE9ZBaz+Gf+melMZ0M25DiO4wf8ANi63qH1uI4+p37Pd03Rk43OqqXNCRHxUrui1TnEYmbWchazmIM5mKM5mtd4jRd4gY3ZmCM5EsKjsFnHmrG3Ti3JyUuSPgyEpI0lfQOYk3NcZmYt4Qqu4B/4B07hFPrSl/705yIuYiEL4cEH6x2eNYhWTE5OAU4DBgHzgRHpuZmZddMzPMMoRq1TtgEbsA/7wDPP1CkqayhN2Funw+QkIpZExAkRsX1EbBcRn4UynfLNzKzT9mRP7uGedcre5m1+z+9hzz3rFJU1mpZLTiTdK2loyfN9gIdzjMnMrGWMYxz/w//wHb7Dn/kzs5nNcRzHB/gADB9e7/CsQbRccgL8X+AOSadKmgBMBL6Qb1hmZq2hH/2YznRe4AWGM5zDOIxd2ZVruKbeoVmDaMYGsbUMwnanpFOAqcASYK+I+EvukZk1CPensO4awhDO4Azex/vYnM05hmPWmWunHqr9dvk9XzxFTzY6q5bLOt8GfgrsD5wL3CvpkznHZWbWEoLgm3yTwziMl3mZx3iM3dmdG7ih3qFZo2jCBrG1DMLWHxgVEW8A90u6A7gMuC3XyMzMWsB93MeN3MgsZrF1GnTtK3yFgzgIXjsYttqqzhFaIyh6stFZtfTWOT0i3pC0eXr+QkR8PP/QzMya33Vcxxf54trEBGAv9uJDfAjuuquOkVkjabaak1ou63woza0zJz333DpmZmYF0YwNYmvprfNjsrl1lgJExONk7U/MzKybPsNnuJRLeYVX1pY9yqPcz/1wyCF1jMwaSbMlJ7W0OSEiXpLWORLPrWNm1gM+xsf4NJ9mD/bgOI5jOcu5iZu4jMs42u1NrBYNkGx0lufWMTOrIyEu4AKmMIUBDGAEI3iKpziKo+odmjWQVqw5OQX4Ce/MrXMXNcytI+kK4HBgcUTskcq2Aa4FhgLzgGMj4pU0Au0c4Om0+QMRcUonjsOsw6nWizg2Q8G/H6yXvMiLfJ/vczu3sxmbsYhFfJtvQ53HOqkkz89atW39eWkdXZpbJyKW1rDvK4Ex7crOAqZFxDBgWnre5tmIGJFuTkzMrCUsZzn7sz9DGcrjPM4UpjCHOYxlbL1DswbSMjUnkn5KlSQ2Ir5abccRMb10Tp7kCGB0ejwJuBc4s+Mwzcya0yQm8WE+zHmcB8AgBnEt1/Ju3g2PP+75daxDbb11mkm1mpMZwEyyesUPAM+k2wi63iB2+4hYCJDutytZtrOkRyXdJ2m/Lu7fzKyhPMETHMAB65RtxEbsx37wxBN1isoaTbPVnFRMTiJiUkRMAoYBB0TETyPip8BBZAlKT1oIDImIvYAzgF9LKttMXdJ4STMkzeDll3s4DDOz3jWMYTzEQ+uUrWEND/MwDBtWp6isoTTh8PW19NbZAdiy5PkWqawrFkkaCJDuFwNExMq2diwRMRN4FnhPuR1ExMSIGBkRIxkwoIthmJkVw0mcxG/5LRdzMStZyRKWcCqnMpCB8MEP1js8axCtmJycDzwq6UpJVwKPAP/Zxde7BRiXHo8DbgaQNEBSn/R4F7Lamue6+BpmZg1jO7ZjKlO5gRvYiq3YiZ14i7e4mZtBBf8FscJotuSkw67EEfELSVOAthT+rIj4S0fbSbqarPFrf0nzgXPIEp3Jkk4GXgSOSavvD5wnaRVZe5ZTImJZZw/GrJqCfxatVQXswR5MZSorV66kT58+bLjhhm2Lqirqe7paXEXs0t/omrFBbE0jxAJ9gJfT+u+R9J6ImF5tg4g4vsKig8qsez1wfY2xmJk1pb59+9Y7BGtQLZecSPoecBwwG1iTigOompyYmZlZL2iAyzSdVUvNyZHAeyNiZd7BmJmZWee1YnLyHLAR4OTEzMysgJotOamlt87rwGOSLpV0Ydst78DMzMysY21tNOXcAAARa0lEQVQNYnuyt46kHSXdI2mOpNmSTk/l20iaKumZdL91yTZnS5or6WlJh5aU7y3pybTsQqnjbmi11Jzckm5mZmZWQDnUnKwCvh4Rj0jaEpgpaSrwebI58s6XdBbZHHlnStoNGAvsTjYW2t2p88xq4GJgPPAAcDvZvHtTqr14LV2JJ3X50MzMzCxfOTSITVPMtE03s0LSHGAQlefIOwK4JrVPfV7SXGCUpHnAVhFxP4Ckq8jasnYtOZE0OSKOlfQkZbqmR8SetR+mWfNqsku9Zrny5yUfebY5SZP47gU8SLs58iS1zZE3iKxmpM38VPZ2ety+vKpqNSenp/vDa4jdzMzM6qSLyUl/STNKnk+MiImlK0jagmwcsq9FxGtVmouUWxBVyquqlpz8BSAiXqgYiaSI8IB/ZmZmddKNEWKXRMTISgslbUSWmPwqIm5IxYskDUy1JmvnyCOrEdmxZPPBwIJUPrhMeVXVeuvcI+krkoa0C3ZjSQdKmsQ78+SYmZlZk0g9ai4H5kTED0sWlZ0jL5WPldRX0s5kc+Q9lC4BrZC0b9rniSXbVFSt5mQMcBJwdXqh5cAmZEPZ3wX8KCIeq/E4zczMLA/5jBD7EeBzwJOS2n7r/5UKc+RFxGxJk4GnyHr6nJZ66gB8CbgS2JSsIWzVxrBQJTmJiDeBi4CLUtVOf+CNiFje2SM0MzOz/OTQW+cPVG6/vN4ceWmbCcCEMuUzgD068/o1TfwXEW+TuhSZmZlZsTTbCLG1zkpsZmZmBeXkxMxaXl5d9Jrs+9WsV3Sjt05hVeytI+kP6X6FpNfa3V6V9LykU3svVDMzM1tPF+bVKXoyU61B7EfT/ZbllkvaFvgfskazZmZmVidFTzY6q8uXdSJiqaTRPRiLmZmZdYGTkxJt4+ubmZlZ/Tg5MTMzs8JoxgaxTk7MzMwaWQM0cO0sJydmVhgddVFusu9fsx7j5MTMzMwKxcmJmZmZFYbbnJiZmVnhODkxMzOz4mjCBrEVh683MzMzqwfXnJiZmTU415zUSNIVkhZLmlVSto2kqZKeSfdblyw7W9JcSU9LOjSvuMzMzJpNs038l+dlnSuBMe3KzgKmRcQwYFp6jqTdgLHA7mmbiyT1yTE2M2tAkePNrFG19dZxclKDiJgOLGtXfAQwKT2eBBxZUn5NRKyMiOeBucCovGIzMzNrJk5Oumf7tskC0/12qXwQ8FLJevNTmZmZmVXThcSk6MlJURrEljtNZWtaJY0HxgMwZEiOIZmZmTWGoicbndXbNSeLJA0ESPeLU/l8YMeS9QYDC8rtICImRsTIiBjJgAG5BmtmZtYImq3mpLeTk1uAcenxOODmkvKxkvpK2hkYBjzUy7GZmZk1nGZsEJvbZR1JVwOjgf6S5gPnAOcDkyWdDLwIHAMQEbMlTQaeAlYBp0XE6rxiMzMzayZFTzY6K7fkJCKOr7DooArrTwAm5BWPmZlZU2qAmpDOKkqDWDMzM+siJydmZmZWKE5OzMzMrDDaGsQ2EycnZmZmDc7JiZmZmRWHG8SamZlZ0TRbctLbg7CZmZmZVeWaEzNbT9mJrcyssJqt5sTJiZmZWQNzbx0zMzMrHCcnZmZmVhzurWNmZmZF4+TEzMzMCsXJiZmZmRWGG8SamTWoJvvuNluHkxMzMzMrDjeINTMzs6JptuTEw9ebmZk1uFDnbx2RdIWkxZJmlZRtI2mqpGfS/dYly86WNFfS05IOLSnfW9KTadmFkjp8dScnZmZmDaytQWxPJyfAlcCYdmVnAdMiYhgwLT1H0m7AWGD3tM1FkvqkbS4GxgPD0q39Ptfj5MTMzKzB5ZGcRMR0YFm74iOASenxJODIkvJrImJlRDwPzAVGSRoIbBUR90dEAFeVbFOR25yYmZk1sq43iO0vaUbJ84kRMbGDbbaPiIUAEbFQ0napfBDwQMl681PZ2+lx+/KqnJyYmZk1uC4mJ0siYmQPhVAugqhSXpWTE7MC6/AT3A1N1ri/JXX0/vDfuHX0Ym+dRZIGplqTgcDiVD4f2LFkvcHAglQ+uEx5VW5zYmZmZrW6BRiXHo8Dbi4pHyupr6SdyRq+PpQuAa2QtG/qpXNiyTYVuebEzMysgeU1fL2kq4HRZG1T5gPnAOcDkyWdDLwIHAMQEbMlTQaeAlYBp0XE6rSrL5H1/NkUmJJuVTk5MTMza3B5JCcRcXyFRQdVWH8CMKFM+Qxgj868tpMTMzOzRubh683MzKxonJyYmZlZoTRbclKX3jqSTpc0S9JsSV9LZedK+rOkx9LtsHrEZtabooObWXf4vdUachy+vm56veZE0h7APwGjgLeAOyTdlhb/KCK+39sxmZmZNbKiJxudVY/LOrsCD0TE6wCS7gM+XYc4zMzMGl8D1IR0Vj0u68wC9pe0raTNgMN4Z1S5L0t6Ik3TvHXlXZiZmVmbZrus0+vJSUTMAb4HTAXuAB4nG7DlYuDdwAhgIfCDcttLGi9phqQZvPxy7wRtZmZWYE5OekBEXB4RH4iI/cmmY34mIhZFxOqIWAP8F1mblHLbToyIkRExkgEDejNsMzOzwnGD2B4iabuIWCxpCHAU8KG2iYTSKp8mu/xjZmZmHSh6stFZ9Rrn5HpJ2wJvk42//4qk/5Y0giwJnAd8sU6xmZmZNY4GqAnprLokJxGxX5myz9UjFrNW5bEuzJqHkxMzMzMrFCcnZmZmVhhtDWKbSV1665iZmZlV4poTMzOzBtdsNSdOTszMzBqZe+uYmZlZ0Tg5MTMzs0JxcmJmneLxRIqho79DXt/t/vtb3pqxt46TEzMzswbn5MTMzMyKww1izczMrGicnJiZmVmhODkxMzOzwnCDWDMzMyscJydmZk2oWpffJvvet2bjBrFmZmZWNE5OzMzMrFCcnJiZmVlhuEGsmZmZFYvbnJiZmVnRNFtyskG9AzAzMzMr5ZoTMzOzBtdsNSdOTsxq4GnvW1tHf/9G/F1oxmNqZU5OzMzMrDDcW8fMzMyKxb11zMzMrGicnJiZmVmhODkxMzOzQnFyYmZmZoXRjA1i6zIIm6TTJc2SNFvS11LZNpKmSnom3W9dj9jMzDorqtzMcpcaxHb21uFupTGSnpY0V9JZ+R/IO3o9OZG0B/BPwChgOHC4pGHAWcC0iBgGTEvPzczMrAM9nZxI6gP8HPgEsBtwvKTd8j+STD1qTnYFHoiI1yNiFXAf8GngCGBSWmcScGQdYjMzM2s4OdScjALmRsRzEfEWcA3Z73SvqEdyMgvYX9K2kjYDDgN2BLaPiIUA6X67OsRmZmbWcHJITgYBL5U8n5/KekWvN4iNiDmSvgdMBf4KPA6sqnV7SeOB8enpSqRZPR9l0+oPLKl3EA1inXPVZG3N8uD3Vu18rjqnEc/XTr36ajNn3onUvwtbbiJpRsnziRExMT0u97XXa82o6tJbJyIuBy4HkPSfZBnZIkkDI2KhpIHA4grbTgQmpm1nRMTIXgq74fl81c7nqnN8vmrnc9U5Pl8di4gxOex2PtlVjTaDgQU5vE5Z9eqts126HwIcBVwN3AKMS6uMA26uR2xmZmbGw8AwSTtL2hgYS/Y73SvqNc7J9ZK2Bd4GTouIVySdD0yWdDLwInBMnWIzMzNraRGxStKXgTuBPsAVETG7t16/Xpd19itTthQ4qJO7mtjxKlbC56t2Pled4/NVO5+rzvH5qpOIuB24vR6vrQgPE2RmZmbFUZc2J2ZmZmaVNExyImlHSfdImpOGvT89lXvY+3YkbSLpIUmPp3P13VTuc1WBpD6SHpX02/Tc56oCSfMkPSnpsbZuiD5flUnqJ+k6SX9K318f8vlan6T3pvdU2+01SV/zuWpNDZOckI2F8vWI2BXYFzgtDaXrYe/XtxI4MCKGAyOAMZL2xeeqmtOBOSXPfa6qOyAiRpR08fT5quwnwB0R8T6yKTvm4PO1noh4Or2nRgB7A68DN+Jz1ZIaJjmJiIUR8Uh6vILsAz4ID3u/nsj8NT3dKN0Cn6uyJA0GPglcVlLsc9U5Pl9lSNoK2J80rlNEvBURy/H56shBwLMR8QI+Vy2pYZKTUpKGAnsBD+Jh78tKlykeIxvMbmpE+FxV9mPgX4A1JWU+V5UFcJekmWnEZvD5qmQX4GXgF+my4WWSNsfnqyNjyca/Ap+rltRwyYmkLYDrga9FxGv1jqeoImJ1qh4dDIxKs0FbO5IOBxZHxMx6x9JAPhIRHyCbrfQ0SfvXO6AC2xD4AHBxROwF/A1flqgqDfj1KeA39Y7F6qehkhNJG5ElJr+KiBtS8aI03D3Vhr1vVakK+V5gDD5X5XwE+JSkeWSzbh4o6Zf4XFUUEQvS/WKyNgGj8PmqZD4wP9VcAlxHlqz4fFX2CeCRiFiUnvtctaCGSU4kiey67ZyI+GHJIg97346kAZL6pcebAgcDf8Lnaj0RcXZEDI6IoWRVyb+LiM/ic1WWpM0lbdn2GDiEbKZxn68yIuIvwEuS3puKDgKewuermuN555IO+Fy1pIYZhE3SR4HfA0/yTtuAfyVrdzIZGEIa9j4iltUlyIKQtCdZw7E+ZAno5Ig4L00Z4HNVgaTRwDci4nCfq/Ik7UJWWwLZJYtfR8QEn6/KJI0ga2y9MfAc8AXS5xKfr3VI2gx4CdglIl5NZX5vtaCGSU7MzMysNTTMZR0zMzNrDU5OzMzMrFCcnJiZmVmhODkxMzOzQnFyYmZmZoXi5MTMzMwKxcmJWU4k/U+F8k0l3ZfmP9pB0nWd3O/nJf2sZ6Ks+BrzJPUvU36MpDmS7sn59SXpQklzJT0h6QOpfGNJ0yVtmOfrm1l9OTkxy0lEfLjCopOAG9L8Rwsi4jO9GVc3nQycGhEHlBbmkCx8AhiWbuOBiyGb1ReYBhzXw69nZgXi5MQsJ5L+WmHRCaQhuCUNlTQrPf68pBsk3SHpGUn/r2RfX5D0v5LuI5sPqK18gKTrJT2cbh9J5RdK+k56fGiqbdigyvrbSrorzZx7KaAyx/Md4KPAJZIuSPH+RtKtZLMUK5XPkvSkpOPSdqNTTdHkdAznSzpB0kNpvXeXOUdHAFdF5gGgX9v8KsBN6RyaWZNy1ahZL0ozru4SEfMqrDIC2AtYCTwt6afAKuC7wN7Aq8A9wKNp/Z8AP4qIP0gaAtwJ7Eo28+3Dkn4PXAgcFhFrJFVa/xzgD2mag0+S1VasIy07kGyI/xmSPg98CNgzIpZJOjrFPxzon15/etp8eHqdZWRDuF8WEaMknQ58Bfhau5cbRDaMeZv5qWwh2Vw++1Q4f2bWBJycmPWu/sDyKsunlcwp8hSwU9rm3oh4OZVfC7wnrX8wsFs2LyYAW0naMiJWSPonYDrwzxHxbLX1gf2BowAi4jZJr9R4PFNL5jn5KHB1RKwmm0n2PrIk4jXg4YhYmOJ/FrgrbfMkcADrW6/mBogU32pJb7UdZ41xmlkDcXJi1rveADapsnxlyePVvPMZrTQJ1gbAhyLijTLL3g8sBXboaP2UrHRloq2/le6mynqlx7Wm5Pkayn8PzQd2LHk+GFhQ8rwv8GbtYZpZI3GbE7NeFBGvAH0kVUtQ2nsQGJ3ahWwEHFOy7C7gy21P0gy4SNoJ+DrZJaJPSPpgtfXJalhOSGWfALbuzHGV7OO41AtpAFltzENd2A/ALcCJqR3LvsCrJTUv2wIvR8TbXdy3mRWckxOz3ncX2SWQmqQf5XOB+4G7gUdKFn8VGJm62z4FnKKsGuRysrYhC8h62FyWEqL11k/7+S6wv6RHgEPIpqbvrBuBJ4DHgd8B/xIRf+nCfgBuJ2ubMhf4L+DUkmUHpOVm1qQU0ZWaXDPrKkl7AWdExOfqHUsjknQDcHZEPF3vWMwsH645MetlEfEocI+kPvWOpdGk3k43OTExa26uOTEzM7NCcc2JmZmZFYqTEzMzMysUJydmZmZWKE5OzMzMrFCcnJiZmVmh/H+rvVQLAtDHiAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# damping rate plots\n",
    "margin = 15  # grid cells\n",
    "if len(damp_time) > 0:\n",
    "    imin = min(damp_table_ii) - margin\n",
    "    imax = max(damp_table_ii) + margin\n",
    "    jmin = min(damp_table_jj) - margin\n",
    "    jmax = max(damp_table_jj) + margin\n",
    "    plt.figure(figsize=(10,5))\n",
    "    plt.imshow(damp_time_map, origin='lower', cmap=plt.get_cmap('cool'))\n",
    "    plt.colorbar()\n",
    "    plt.scatter(damp_table_ii,damp_table_jj, facecolors='none', edgecolors='k')\n",
    "    plt.xlim(imin,imax)\n",
    "    plt.ylim(jmin,jmax)\n",
    "    plt.title('Rayleigh damping time (s)')\n",
    "    plt.xlabel('i (indexed from 0)')\n",
    "    plt.ylabel('j (indexed from 0)')\n",
    "#     savefigure('Rayleigh_damp_time')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhQAAAFNCAYAAABYE/ZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xm4ZGV57/3vj9lGBIQWmQcDKqKAtkSjIopRNEZURCGoOJwgrzgdo0dJXgU9IYe8JnHWSARBjSgyiQcUGiMQYxAbQZlFFOyWlmYUtQVpuN8/ajVUN3vXrt1r165de38/11VXVT1rrVp3Pbuq+u5nPUOqCkmSpDbWGnYAkiRp9JlQSJKk1kwoJElSayYUkiSpNRMKSZLUmgmFJElqzYRCmqGSvCHJ94Ydx0SSHJ3ky1P8mqu89ySV5E+m8hySppYJhSRJas2EQpIktWZCIU2h1Zvmk5yY5O+bx/skWZLkb5IsS7I0yRu79t0syVlJ7k5yCfC41V7740kWN9svTfKcrm1HJ/l6ki8n+W2SK5LskuTI5lyLk7ywa/8LkvyfJJck+U2SbyR5dNf2ZyT5fpK7kvw4yT5d23ZMcmFznoXA5hPUyf5JLm/iviHJfk35xkmOb+rhV0n+PsnafdTxS5Jc3Zz/V0neM9ExkgbPhEKaXo8FNga2Bt4MfDrJps22TwP3AFsCb2pu3X4I7AE8GvgK8PUkG3Rt/0vgS8CmwGXAuXS+41sDHwY+t9rrvb45x1bACuATAEm2Bs4G/r4513uA05LMb477CnApnUTifwOHjvdmk+wFfBF4L7AJsDdwY7P5pOa8fwLsCbwQ+B/jvVaX44G3VNVGwG7Af/RxjKQBM6GQptd9wIer6r6qOgf4HfD45n/mBwAfrKrfV9WVdP7BfVBVfbmqbq+qFVX1z8D6wOO7dvnPqjq3qlYAXwfmA8dW1X3AV4EdkmzStf+XqurKqvo98AHg1U0crwXOqapzquqBqloILAJekmQ74OnAB6rq3qq6CPhmj/f7ZuCEqlrYvNavquraJFsALwbe1bzfZcBHgYP6rMNdkzyqqu6sqh/1cYykATOhkKbX7c0/+CstBx5J5x//dYDFXdtu6j6wuVRyTXOJ4i46LR3dlxtu6Xr8B+C2qrq/6znNuVZa/VzrNq+3PXBgc7njruZcz6bTcrIVcGeThIwZ52q2BW4Yo3z75nxLu87xOeAxPV5rpQOAlwA3NZdentnHMZIGbJ1hByDNMsuBeV3PHwss6eO4W+k0/28LXNuUbbdyY9Nf4n3AvsBVVfVAkjuBtIh1267H29H5n/9tdBKNL1XVX69+QJLtgU2TbNiVVGwHjLds8WJW6wvSVX4vsPlqCdaEquqHwP5J1gXeBpyy2nuRNAS2UEhT63Lgr5Ks3XQ+fG4/BzUtCacDRyeZl2RXVu2bsBGdhONWYJ0kHwQe1TLW1ybZNck8On0sTm3i+DLwl0le1LyPDZoOpdtU1U10Ln98KMl6SZ5Np+/GeI4H3phk3yRrJdk6yROqailwHvDPSR7VbHtckp711ZzzkCQbN5dy7gbu73WMpOlhQiFNrXfS+Qf2LuAQ4MxJHPs2Opckfg2cCHyha9u5wLeAn9K5xHAPq16yWBNfas7za2AD4B0AVbUY2B/4WzoJzGI6nSpX/l78FfCnwB3AUXQ6XY6pqi4B3kinf8RvgAvpXO6ATqfQ9YCrgTuBU+lcVpnI64Abk9wNHE6nz4ekIUvVeC2VkmarJBcAX66qzw87Fkmzgy0UkiSpNRMKSZLUmpc8JElSa7ZQSJKk1kwoJElSayM9sdW8bF6bsMODz7d62vBi6eXmS4cdwcPNxroa1nsa5N93pv6depmoPkbxPc1UM/G3pa1Bfj6mq76WcultVTV/4j2nRpLakslX3FIuparaTI63ahyj3Idiqyyow1j04POjZ+hbOXrK/lxTZzbW1bDe0yD/vjP179TLRPUxiu9pppqJvy1tDfLzMV319SFyaVUtmJ6zdRKKo8edrHZ8R5MpTSi85CFJklozoZAkSa0NrA9FkhOAlwLLqmq3puwjdKYl/iOdFQjfWFV3JdkBuAa4rjn84qo6fFCxyWZ6aTaY6Ls2ipdEesXsb8vMNsgWihOB/VYrWwjsVlVPobMmwZFd226oqj2am8mEJEkjZGAJRVVdRGfxoO6y87qWKr4Y2GZQ55ckSdNnmH0o3kRn9cSVdkxyWZILkzxnWEFJkqTJG8o8FEn+DlgB/HtTtBTYrqpuT/I04MwkT6qqu8c49jDgMICN2W66QpYkST1MewtFkkPpdNY8pJpJMKrq3qq6vXl8KZ0Om7uMdXxVHVdVC6pqwTymbd4QSZLUw7QmFEn2A94HvKyqlneVz0+ydvN4J2Bn4OfTGZskSVpzgxw2ejKwD7B5kiXAUXRGdawPLEwCDw0P3Rv4cJIVwP3A4VV1x5gvrBlvFIeqtTHIoXtzrS41tXp9Nv1saaoNLKGoqoPHKD5+nH1PA04bVCySJGmwnClTkiS1ZkIhSZJaM6GQJEmtmVBIkqTWTCgkSVJrJhSSJKm1oUy9PSgue6s1NdGYfD8/U6fVvBz+HWY9/8ajyxYKSZLUmgmFJElqzYRCkiS1ZkIhSZJaM6GQJEmtmVBIkqTWZtWwUUmzm8N7p85EdTWKy5uPYsyziS0UkiSpNRMKSZLUmgmFJElqzYRCkiS1ZkIhSZJaM6GQJEmtmVBIkqTWnIdC0qzRax4C56iQBssWCkmS1JoJhSRJas2EQpIktWZCIUmSWjOhkCRJrZlQSJKk1gaWUCQ5IcmyJFd2lX0kybVJfpLkjCSbdG07MsnPklyX5EWDikuSJE29Qc5DcSLwKeCLXWULgSOrakWSfwSOBN6XZFfgIOBJwFbA+Ul2qar7pyqYXuPTwTHq0mznb4A0WANroaiqi4A7Vis7r6pWNE8vBrZpHu8PfLWq7q2qXwA/A/YaVGySJGlqDbMPxZuAbzWPtwYWd21b0pRJkqQRMJSEIsnfASuAf19ZNMZuYzZAJjksyaIki5Zz66BClCRJkzDtCUWSQ4GXAodU1cqkYQmwbddu2wA3j3V8VR1XVQuqasE85g82WEmS1JdpTSiS7Ae8D3hZVS3v2nQWcFCS9ZPsCOwMXDKdsUmSpDU3sFEeSU4G9gE2T7IEOIrOqI71gYVJAC6uqsOr6qokpwBX07kUcsRUjvCQJEmDNbCEoqoOHqP4+B77HwMcM6h4JjLIZY97HT/RUDbNDP6dpkfr75p/pynT5nfLIbhzkzNlSpKk1kwoJElSayYUkiSpNRMKSZLUmgmFJElqzYRCkiS1ZkIhSZJaG+Ty5dK0cg4CaXYb5He8zdwZ/vZ02EIhSZJaM6GQJEmtmVBIkqTWTCgkSVJrJhSSJKk1EwpJktSaw0b7MMileic61uFIkkaNv1tzky0UkiSpNRMKSZLUmgmFJElqzYRCkiS1ZkIhSZJaM6GQJEmtmVBIkqTWnIdihus1T8VA58dwHLmGYJifu17nbvNdkuYKWygkSVJrJhSSJKk1EwpJktSaCYUkSWrNhEKSJLVmQiFJklob2LDRJCcALwWWVdVuTdmBwNHAE4G9qmpRU74DcA1wXXP4xVV1+KBimy0GOZStzXBVSRo1/q61N8gWihOB/VYruxJ4JXDRGPvfUFV7NDeTCUmSRsjAWiiq6qKm5aG77BqAxFRQkqTZZCb1odgxyWVJLkzynGEHI0mS+jdTpt5eCmxXVbcneRpwZpInVdXdq++Y5DDgMICN2W6aw5QkSWOZES0UVXVvVd3ePL4UuAHYZZx9j6uqBVW1YB7zpzNMSZI0jhmRUCSZn2Tt5vFOwM7Az4cblSRJ6tcgh42eDOwDbJ5kCXAUcAfwSWA+cHaSy6vqRcDewIeTrADuBw6vqjsGFZskSZpagxzlcfA4m84YY9/TgNMGFYukqdVmDpSZOt7fJcqldmbEJQ9JkjTaTCgkSVJrJhSSJKk1EwpJktRa3wlFkg1XDu2UJEnqNm5CkWStJH+V5Owky4BrgaVJrkrykSQ7T1+YkiRpJus1bPS7wPnAkcCVVfUAQJJHA88Djk1yRlV9efBhzmxthsGN4lC1iWKeqcMCR9FsHJ6p0edvgMbSK6F4QVXdt3phM+HUacBpSdYdWGSSJGlkjJtQrEwmkmwBbA0UcHNV3bL6PpIkaW4bN6FIsifwWWBj4FdN8TZJ7gLeWlU/mob4JEnSCOh1yeMLwFuq6gfdhUme0WzbfZCBSZKk0dFr2OiGqycTAFV1MbDh4EKSJEmjplcLxbeSnA18EVjclG0LvB749qADkyRJo6NXp8x3JHkxsD+dTpkBlgCfrqpzpik+SZI0AnouX15V3wK+NU2xzEmzcQ6LXnE5Pl3SbLP6b96H5ujvnGt5SJKk1kwoJElSayYUkiSptZ59KACS7Ai8Hdihe/+qetngwpIkSaNkwoQCOBM4Hvgm8MBgw5EkSaOon4Tinqr6xMAjkSRJI6ufhOLjSY4CzgPuXVnoWh5qazm3cSvXsAk7sDHbDjscSVIL/SQUTwZeBzyfhy55VPNcQzTRnA4zdZ6K4gEW8j4u4/PMZ1du4zoexwt5GcdzTD1i2OFNuVZzjczQ8ewzNa65ZjbOY6PR1U9C8Qpgp6r646CD0dxwCZ9mMf/F2/kZ89iM+1jOmbyBhbyXY/jUsMOTJK2BfoaN/hjYZNCBaO64lM/x53yEeWwGwLrMYz8+zk/4En/8o3mrJI2iflootgCuTfJDVu1D4bBRrZHl3MYmbL9K2SPZggdYwT333MN66603pMgkSWuqn4TiqIFHoTllB/bhCk7mWbz3wbLr+Cab8Xg22mijIUYmSVpTEyYUVXVhki2ApzdFl1TVssGGpdlsH47mRPZhObeyE3/OzSziYj7KAZxMYm8/SRpFE/ahSPJq4BLgQODVwA+SvGrQgWn22pwn8D/4AcUD/Bf/yF3cyOv5Djux77BDkyStoX4uefwd8PSVrRJJ5gPnA6f2OijJCcBLgWVVtVtTdiBwNPBEYK+qWtS1/5HAm4H7gXdU1bmTfjdaxaCG9k3FcLNN2J4X8k/tX2hEFA9wNadyDacT1mJXDuQJvJwP1eBaZBzaOTkOwdSa8rvW0c8oj7VWu8Rxe5/HnQjst1rZlcArgYu6C5PsChwEPKk55jNJ1u7jHNKMVxRncijf5yP8CfuxI8/nAo7iHI4YdmiSNGX6aaH4dpJzgZOb568BzpnooKq6KMkOq5VdA4x1nXx/4KtVdS/wiyQ/A/YC/ruP+KQZ7Zd8jyVczOH8hHXpTNy1GwfxKZ7A5Zcfxh577DHkCCWpvQlbGqrqvcDngKcAuwPHVdX7pjiOrYHFXc+XNGXSyPs557Mrr34wmQBYj0fyRF7J+eefP8TIJGnq9GyhaC47nFtVLwBOH2AcY12BGvOqZJLDgMMANma7AYYkTY0N2IRlXPmw8t+xlE033X0IEUnS1OvZQlFV9wPLk2w84DiWwCqrQ20D3DxOTMdV1YKqWjCP+QMOS2rvyRzMdZzJL/mvB8t+zvncyAUccMABQ4xMkqZOX8uXA1ckWQj8fmVhVb1jCuM4C/hKkn8BtgJ2pjNUVRp5j+SxvJJ/5xQOYFN25AFW8Ftu5kBOZZNNnNVe0uzQT0JxdnOblCQnA/sAmydZQmfGzTuATwLzgbOTXF5VL6qqq5KcAlwNrACOaFpHpFnhT9iPd3ETi/kvwlpsy7NYm3WHHZYkTZlxE4ok36mqfYFd16QTZlUdPM6mM8bZ/xjgmMmeR9NvVMdcDzvudVifHXn+cIPQQPT6bM3UOSom+j60ibvXscP+HmpwerVQbJnkucDLknyV1TpOVtWPBhqZJEkaGb0Sig8C76fTQfJfVttW4H+1JElSx7gJRVWdCpya5ANV9b+nMSZJkjRi+pnYymRCkiT11M+aHJIkST2ZUEiSpNZ6DRt9dK8Dq+qOqQ9Hc91sHFI2yGGDs7G+Zpt+/kbFA9zAeVzPOazLPJ7MIWzBkwcfnNRDkv2AjwNrA5+vqmN77d+rheJSYFFzfyvwU+D65vGlUxKtJM1xxQOcwes5j/fwKLZlLdbhS/w5i/jcsEPTHNas5fVp4MXArsDBSXbtdUyvUR47Ni/6r8BZVXVO8/zFwAumKmhJmsuu51ss4woOYxHrsAEAe/ImjuNp7MqrmMdmQ45Qc9RewM+q6ucAzXxU+9OZ0XpM/fShePrKZAKgqr4FPLdloJIk4Kd8kz1404PJBMCm7MQO7MPPWTjEyDTHbQ0s7nq+pCkbVz9redyW5P8FvkxnQqvXArevaYSSpIeswwb8kd8+rPxefss6PGIIEWkuSbKo6+lxVXXcyk1j7N6zR1g/LRQH01nM64zmNr8pkyS19GQOYRH/yt0sebDsF/wHy7iCx/HCIUamuaCqFnTdjuvatATYtuv5NsDNvV5rwhaKZjTHO5M8sqp+t0YRS5LGtDVP55m8m8/yZB7Hi7iHu1jKpRzI11nXFgoNzw+BnZPsCPwKOAj4q14HTJhQJPkz4PPAI4HtkuwOvKWq3to+XknSM3k3u3EQ13MO67ERr+Y01mPDYYelOayqViR5G3AunWGjJ1TVVb2O6acPxUeBFwFnNSf5cZK92wYrTbXZON+D80zMDcu4knP5n9zERazFOjyJ1/BC/plHsOmwQ5tyE31P/czPHM2AjHMm3LHR10yZVbV4taL7JxOUJGlsv2cZX+QFPJEDeD938y5uYm3W42u8nOrdB06aUfpJKBY3lz0qyXpJ3gNcM+C4JGlOuIwT2IW/ZAGHsw7rM4/N+Qs+w29Zyq/4wbDDk/rWT0JxOHAEnfGnS4A9mueSpJZu53q2Zq9VysJabM3TuZ3rhxSVNHn9LF9+W1UdUlVbVNVjquq1MMagaUnSpG3BU7iR765Sdj/3cRP/yRY8ZUhRSZM3YUKR5IIkO3Q9fzqd4SSSpJb24FAW832+ywe5m1+xjKs4ldewJU/lsew+7PCkvvUzyuP/AN9O8gk6lz1eArxxoFFJ0hyxAZvwRi7iu3yAf2V31mNDnsxreS4fGHZo0qT0M7HVuUkOBxYCtwF7VtWvBx6ZNCIGOVxVc8PGbMczeDeb8QTWY0N25cBV1vYYhl7DN/3Mayz9XPL4APBJYG/gaOCCJH8x4LgkaU4oivN4L1/hJSznVn7N5XyGJ3ENpw87NGlS+rnksTmwV1X9AfjvJN+mM3Pm2QONTJLmgJu4kGs5g/+HKx+cyGov3s4X2ZedeAHr86ghRyj1p59RHu+sqj8k2bB5flNV/fngQ5Ok2e9qTuVpvGWVWTG3ZE+25ZncwHlDjEyanH4ueTwzydU0k1kl2T3JZwYemSRJGhn9TGz1MTpredwOnbU86PSnkCS1tCuv4lI+xx+488GypVzGYv7b5cs1UvrpQ0FVLU5W6fLrWh6SNAW257k8gVfwWXbjSbyGe7iLazmTl/F5+09opPSTUKyylgfwDlzLQ5KmRAgv5CPszuv4KWezCTuwL//AI3nssEOTJqWfhOJw4OM8tJbHefSxlkeSE4CXAsuqarem7NHA14AdgBuBV1fVnc1MnNcA1zWHX1xVh0/ifUgTLns8E8fOu1SzAH7DL/k+/8T1nMO6zON33MJz+cDQ56IYzyC/a72O9fsys63RWh5VdXsfr30isN9qZe8HvlNVOwPfaZ6vdENV7dHcTCYkzQn3cBdfYG82YQcO58ccwre4jWs4lYOGHZo0KeO2UCT5JDBurlhV7+j1wlV1UfcaII39gX2axycBFwDvmzhMSZqdLucktuXPeB4fBuBRbM2r+Bqf4HH8mh+7nodGRq8WikXApcAGwFOB65vbHqx5p8wtqmopQHP/mK5tOya5LMmFSZ6zhq8vSSPlFn7CDjxvlbK1WZfteQ638JMhRSVN3rgJRVWdVFUnATsDz6uqT1bVJ4F96SQVU2kpsF1V7Qm8G/hKkjG7Nyc5LMmiJIuWc+sUhyFJ02szduZXXLJKWfEAv+KHbMbOQ4pKmrx+5qHYCtio6/kjm7I1cUuSLQGa+2UAVXXvyn4ZVXUpcAOwy1gvUFXHVdWCqlowj/lrGIYkzQx78iau5//yQz7LCu5lObdxNm9lI7Zka/502OFJfesnoTgWuCzJiUlOBH4E/MManu8s4NDm8aHANwCSzE+ydvN4JzqtIj9fw3NI0sjYkMfwOhZyLadzLI/iY2zP/fyRg/gGwWENGh39LF/+hSTfggdT5ff3s3x5kpPpdMDcPMkS4Cg6yckpSd4M/BI4sNl9b+DDSVbQ6Z9xeFXdMdk3I/XikDPNRJ1hkrsBC7n33ntZe+21WWedvuYcnLGfaZc+n5v6+9TC2sCtzf67JNmlqi7qdUBVHTzOpn3H2Pc04LQ+Y5GkWWn99dcfdgjSGpswoUjyj8BrgKuAB5riAnomFJIkae7op4Xi5cDjq+reQQcjSZJGUz+dMn8OrDvoQCRJ0ujqp4ViOXB5ku8AD7ZSTDRTpiRJmjv6SSjOam6SJElj6mfY6EnTEYgkSRpdvRYHO6WqXp3kCsZYJKyqnjLQyKQRMVPnApBmIr8vs1evFop3NvcvnY5AJEnS6OqVUPwaoKpuGm+HJKkq5z2TJGmO6zVs9LtJ3p5ku+7CJOsleX6Sk3hoXQ5JkjSH9Wqh2A94E3Bykh2Bu4AN6EzDfR7w0aq6fPAhSpKkmW7chKKq7gE+A3wmybrA5sAfququ6QpOkiSNhr4WB6uq+4ClA45FkiSNqH6m3pYkSeqp3+XLJelBRw9obJdzFEija9wWiiTfa+5/m+Tu1W6/SfKLJG+dvlAlSdJM1atT5rOb+43G2p5kM+D7dDpuSpKkOWyN+1BU1e3APlMXiiRJGlWtOmVWlSM/JEmSozwkSVJ7JhSSJKk1h41KmjEmGo7qsFJp5rKFQpIktWZCIUmSWjOhkCRJrZlQSJKk1kwoJElSayYUkiSptYElFElOSLIsyZVdZY9OsjDJ9c39pl3bjkzysyTXJXnRoOKSJElTb5DzUJwIfAr4YlfZ+4HvVNWxSd7fPH9fkl2Bg4AnAVsB5yfZparuH2B8kkbMoJZNB+e4kNoaWAtFVV0E3LFa8f7ASc3jk4CXd5V/tarurapfAD8D9hpUbJIkaWpNdx+KLVYuKNbcP6Yp3xpY3LXfkqZMkiSNgJnSKXOsxsYxGzeTHJZkUZJFy7l1wGFJkqR+THdCcUuSLQGa+2VN+RJg2679tgFuHusFquq4qlpQVQvmMX+gwUqSpP5Md0JxFnBo8/hQ4Btd5QclWT/JjsDOwCXTHJskSVpDAxvlkeRkYB9g8yRLgKOAY4FTkrwZ+CVwIEBVXZXkFOBqYAVwhCM8JEkaHQNLKKrq4HE27TvO/scAxwwqHkmSNDgzpVOmJEkaYSYUkiSpNRMKSZLUmgmFJElqzYRCkiS1ZkIhSZJaM6GQJEmtDXL5ckkjapDLhEuanWyhkCRJrZlQSJKk1kwoJElSayYUkiSpNRMKSZLUmgmFJElqzWGjkuaEozPsCKTZzRYKSZLUmgmFJElqzYRCkiS1ZkIhSZJaM6GQJEmtmVBIkqTWTCgkSVJrzkMhzWCDXEbceRlG30SfD//Gmk62UEiSpNZMKCRJUmsmFJIkqTUTCkmS1JoJhSRJas2EQpIktTaUYaNJ3gn8NRDg36rqY0mObspubXb726o6ZxjxSdNlkMNCpV6fL4eUaqpNe0KRZDc6icNewB+Bbyc5u9n80ar6p+mOSZIktTOMFoonAhdX1XKAJBcCrxhCHJIkaYoMow/FlcDeSTZLMg94CbBts+1tSX6S5IQkmw4hNkmStAamPaGoqmuAfwQWAt8GfgysAD4LPA7YA1gK/PNYxyc5LMmiJIuWP9jdQpIkDdNQRnlU1fFV9dSq2hu4A7i+qm6pqvur6gHg3+j0sRjr2OOqakFVLZjH/OkMW5IkjWMoCUWSxzT32wGvBE5OsmXXLq+gc2lEkiSNgGGtNnpaks2A+4AjqurOJF9KsgdQwI3AW4YUmyRJmqShJBRV9Zwxyl43jFikuco5MCRNJWfKlCRJrZlQSJKk1kwoJElSayYUkiSpNRMKSZLUmgmFJElqzYRCkiS1NqyJraQ5w/keZoaJ/g5HZzjnlWYLWygkSVJrJhSSJKk1EwpJktSaCYUkSWrNhEKSJLVmQiFJklpz2Kgk0Xt456CGlEqziS0UkiSpNRMKSZLUmgmFJElqzYRCkiS1ZkIhSZJaM6GQJEmtmVBIkqTWnIdC6oNLUM9tw1r6fJBm43vScNlCIUmSWjOhkCRJrZlQSJKk1kwoJElSayYUkiSpNRMKSZLU2lCGjSZ5J/DXQIB/q6qPJXk08DVgB+BG4NVVdecw4pOkyXBYsTSEFooku9FJJvYCdgdemmRn4P3Ad6pqZ+A7zXNJkjQChnHJ44nAxVW1vKpWABcCrwD2B05q9jkJePkQYpMkSWtgGAnFlcDeSTZLMg94CbAtsEVVLQVo7h8zhNgkSdIaSNX0X/xL8mbgCOB3wNXAH4A3VtUmXfvcWVWbjnHsYcBhzdPd6CQo6s/mwG3DDmJEWFeTY331z7qanFGsr+2rav50nSxJAZeuwaFPq6opm2R9KAnFKgEk/wAsAd4J7FNVS5NsCVxQVY+f4NhFVbVgOuKcDayv/llXk2N99c+6mhzra3QMZdhoksc099sBrwROBs4CDm12ORT4xjBikyRJkzes1UZPS7IZcB9wRFXdmeRY4JTmcsgvgQOHFJskSZqkoSQUVfWcMcpuB/ad5EsdNzURzRnWV/+sq8mxvvpnXU2O9TUiht6HQpIkjT6n3pYkSa2NTEKRZNsk301yTZKrmum7SfLoJAuTXN/cP2yo6VyTZIMklyT5cVNXH2rKratxJFk7yWVJ/m/z3LoaR5Ibk1yR5PIki5oy62scSTZJcmqSa5vfr2daXw+X5PHNZ2rl7e4k77KuRsfIJBTACuBvquqJwDOAI5LsilN2j+Ve4PlVtTuwB7BfkmdgXfXyTuCarufWVW/Pq6o9uobzWV/j+zjw7ap6Ap3lBq7B+nqYqrqu+UztATwNWA6cgXU1MkYmoaiqpVX1o+bxb+l8KbfGKbvlw4YdAAAFz0lEQVQfpjp+1zxdt7kV1tWYkmwD/AXw+a5i62pyrK8xJHkUsDdwPEBV/bGq7sL6msi+wA1VdRPW1cgYmYSiW5IdgD2BH+CU3WNqmvAvB5YBC6vKuhrfx4D/BTzQVWZdja+A85Jc2sxcC9bXeHYCbgW+0FxS+3ySDbG+JnIQnfmJwLoaGSOXUCR5JHAa8K6qunvY8cxUVXV/03S4DbBXs8qrVpPkpcCyqlqTaWvnqmdV1VOBF9O59Lj3sAOawdYBngp8tqr2BH6PTfY9JVkPeBnw9WHHoskZqYQiybp0kol/r6rTm+Jbmqm6ae6XDSu+mahpXr0A2A/raizPAl6W5Ebgq8Dzk3wZ62pcVXVzc7+MzjXuvbC+xrMEWNK0EAKcSifBsL7G92LgR1V1S/PcuhoRI5NQJAmd65DXVNW/dG1yyu7VJJmfZJPm8SOAFwDXYl09TFUdWVXbVNUOdJpZ/6OqXot1NaYkGybZaOVj4IV0FuizvsZQVb8GFidZuS7RvnQWRLS+xncwD13uAOtqZIzMxFZJng38J3AFD13r/ls6/ShOAbajmbK7qu4YSpAzRJKn0Om8tDadpPGUqvpwM925dTWOJPsA76mql1pXY0uyE51WCeg053+lqo6xvsaXZA86HX7XA34OvJHme4n1tYok84DFwE5V9ZumzM/WiBiZhEKSJM1cI3PJQ5IkzVwmFJIkqTUTCkmS1JoJhSRJas2EQpIktWZCIUmSWjOhkAYkyffHKX9Ekgub9Va2SnLqJF/3DUk+NTVRjnuOG5NsPkb5gc0S3N8d8PmT5BNJfpbkJ0me2pSvl+SiJOsM8vySJs+EQhqQqvqzcTa9CTi9WW/l5qp61XTG1dKbgbdW1fO6CwfwD/yLgZ2b22HAZ6GzWiedJaxfM8Xnk9SSCYU0IEl+N86mQ2imD06yQ5Irm8dvSHJ6km8nuT7J/9f1Wm9M8tMkF9JZf2Rl+fwkpyX5YXN7VlP+iSQfbB6/qPlf/Vo99t8syXnNipifAzLG+/kg8GzgX5N8pIn360m+SWf10TTlVya5IslrmuP2aVpkTmnew7FJDklySbPf48aoo/2BL1bHxcAmK9dzAM5s6lDSDGKzoTSNmpUUd6qqG8fZZQ9gT+Be4LoknwRWAB8Cngb8BvgucFmz/8eBj1bV95JsB5wLPJHOipY/TPKfwCeAl1TVA0nG2/8o4HvNFO1/QadVYBXNtufTmZ58UZI3AM8EnlJVdyQ5oIl/d2Dz5vwXNYfv3pznDjrTT3++qvZK8k7g7cC7Vjvd1nSmYF5pSVO2lM7aIU8fp/4kDYkJhTS9Ngfu6rH9O11rGFwNbN8cc0FV3dqUfw3Ypdn/BcCunbXzAHhUko2q6rdJ/hq4CPifVXVDr/2BvYFXAlTV2Unu7PP9LOxaV+HZwMlVdT+dFSIvpPMP/93AD6tqaRP/DcB5zTFXAM/j4R7WQgJUE9/9Sf648n32GaekATOhkKbXH4ANemy/t+vx/Tz0HR1v0Z21gGdW1R/G2PZk4HZgq4n2bxKMNVnY5/fdL9Njv+739UDX8wcY+3doCbBt1/NtgJu7nq8P3NN/mJIGzT4U0jSqqjuBtZP0SipW9wNgn6afw7rAgV3bzgPetvJJs7IlSbYH/obO5ZMXJ/nTXvvTack4pCl7MbDpZN5X12u8phm9Mp9Oq8cla/A60Fmy+vVNv4xnAL/pauHYDLi1qu5bw9eWNAAmFNL0O4/O5YG+NP+QHg38N3A+8KOuze8AFjRDK68GDk+nueF4On0dbqYzMuPzTRLzsP2b1/kQsHeSHwEvpLNM9GSdAfwE+DHwH8D/qqpfr8HrAJxDp6/Fz4B/A97ate15zXZJM4jLl0vTLMmewLur6nXDjmUUJTkdOLKqrht2LJIeYguFNM2q6jLgu0nWHnYso6YZJXOmyYQ089hCIUmSWrOFQpIktWZCIUmSWjOhkCRJrZlQSJKk1kwoJElSa/8/2k1hiT7DDKkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# missing cell plots\n",
    "margin = 15  # grid cells\n",
    "maxundamped = int(undamped.max())\n",
    "\n",
    "if len(damp_time) > 0:\n",
    "    imin = min(damp_table_ii) - margin\n",
    "    imax = max(damp_table_ii) + margin\n",
    "    jmin = min(damp_table_jj) - margin\n",
    "    jmax = max(damp_table_jj) + margin\n",
    "    plt.figure(figsize=(10,5))\n",
    "    plt.imshow(undamped, origin='lower', cmap=plt.get_cmap('rainbow',maxundamped+1))\n",
    "    plt.colorbar(ticks=range(maxundamped+1))\n",
    "    plt.clim(-0.5, maxundamped+0.5)\n",
    "    plt.scatter(damp_table_ii,damp_table_jj, facecolors='none', edgecolors='k')\n",
    "    plt.xlim(imin,imax)\n",
    "    plt.ylim(jmin,jmax)\n",
    "    plt.title('undamped cells')\n",
    "    plt.xlabel('i (indexed from 0)')\n",
    "    plt.ylabel('j (indexed from 0)')\n",
    "#     savefigure('Rayleigh_undamped_cells')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:analysis3-19.01]",
   "language": "python",
   "name": "conda-env-analysis3-19.01-py"
  },
  "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}