Jint_Ptot conservation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "import xarray as xr\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%pylab inline\n",
    "pylab.rcParams['figure.figsize'] = (18.0, 10.0); # Large figures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# gx3v7 output; nstep output\n",
    "# Tracing reason for error message\n",
    "# Message from (lon, lat) ( 302.489,  42.997), which is global (i,j) (95, 80)\n",
    "filename = '/Users/mlevy/C1850ECO.T62_g37.passes_Pmassbal_nstep_out.pop.h.0001-01-05-24192.nc'\n",
    "ds = xr.open_dataset(filename, decode_times=False, decode_coords=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "302.488515327 42.9973933488\n"
     ]
    }
   ],
   "source": [
    "# Variables to keep in dataset\n",
    "grid_vars = ['TLAT', 'TLONG', 'z_t', 'z_t_150m', 'z_w_bot', 'dz']\n",
    "POP_vars = ['Jint_Ctot', 'Jint_Ptot', 'POP_PROD']\n",
    "vars_to_keep = grid_vars + POP_vars\n",
    "\n",
    "# Fortran error message is indexed starting from 1\n",
    "# Python indices start from 0\n",
    "lat_in = 79\n",
    "lon_in = 94\n",
    "\n",
    "raw_data = ds.isel(time=0).drop([v for v in ds.variables if v not in vars_to_keep])\n",
    "single_col = ds.isel(time=0, nlat=lat_in, nlon=lon_in).drop([v for v in ds.variables if v not in vars_to_keep])\n",
    "print single_col['TLONG'].values, single_col['TLAT'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Volumes/OriginalHD/Users/mlevy/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/ipykernel/__main__.py:1: RuntimeWarning: divide by zero encountered in log10\n",
      "  if __name__ == '__main__':\n",
      "/Volumes/OriginalHD/Users/mlevy/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/ipykernel/__main__.py:1: RuntimeWarning: invalid value encountered in log10\n",
      "  if __name__ == '__main__':\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVMAAAJRCAYAAAD4a1TiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+UFOd5J/rv0z0F9CBZDTaOxRiEpCQoxhjGGlvYOJtI\ncURuLJGxkIS1Inuc3F0lN3uSI8XBgUjHhr1OhHfi2Lvrk7vRvSe5dyOiYBtrLC2ycWS0sYODZNAM\nQthiV7IkpJYcEcEQwTTQM/PcP7qrqa6u3/12V3XP93MOB6a6q/qdHvqZt973eZ9XVBVERNSaXNoN\nICLqBQymREQGMJgSERnAYEpEZACDKRGRAQymREQGMJgSERnAYEpEZACDKRGRAX1pN8CUd7zjHbps\n2bK0m0FEXejQoUP/rKqLWrlGzwTTZcuW4eDBg2k3g4i6kIi83Oo1eJtPRGQAgykRkQEMpkREBjCY\nEhEZwGBKRGQAgykRkQEMpkREBjCYEhEZwGBKRGQAgykRkQEMpkREBjCYEhEZwGBKRGQAgykRkQEM\npkREBjCYEhEZwGBKRGQAgykRkQEMpkREBjCYEhEZwGBKRGQAgykRkQEMpkREBjCYEhEZwGBKRGQA\ngykRkQEMpkREBjCYEhEZwGBKRGQAgykRkQEMpkREBjCYEhEZwGBKRGRAX9oNIOomo2MljOw9htcm\nylhcLGDzuuUYHhxIu1mUAQymRBGNjpWw9etHUK5MAwBKE2Vs/foRAGBAJd7mE0U1svdYPZDaypVp\njOw9llKLKEsYTIkiem2iHOs4zS4MpkQRLS4WYh2n2YXBlCiizeuWo2DlG44VrDw2r1ueUosoSzgB\nRRSRPcnE2XzyklowFZHbAGwD8HMAPqiqB12PLwXwQwDbVPVPO99CombDgwMMnuQpzdv8ZwHcAuC7\nPo//GYBvdq45RETJpdYzVdUfAYCIND0mIsMAXgRwtsPNIiJKJHMTUCJyCYA/BLA97bYQEUXV1p6p\niDwO4F0eD92rqt/wOW0bgC+q6hmvXqvr+ncBuAsAli5d2kJLiYha09ZgqqofTXDadQBuFZH/CKAI\nYEZEzqnqlz2u/wCABwBgaGhIW2osEVELMpcapao/b/9bRLYBOOMVSImIsiS1MVMR+biIvArgQwD2\niMjetNpCRNSqNGfzHwbwcMhztnWmNURErcncbD4RUTdiMCUiMoDBlIjIAAZTIiIDGEyJiAxgMCUi\nMoDBlIjIAAZTIiIDGEyJiAxgMCUiMoDBlIjIAAZTIiIDGEyJiAxgMCUiMoDBlIjIAAZTIiIDGEyJ\niAxgMCUiMoDBlIjIAAZTIiIDGEyJiAxgMCUiMoDBlIjIAAZTIiIDGEyJiAxgMCUiMoDBlIjIAAZT\nIiIDGEyJiAxgMCUiMoDBlIjIAAZTIiIDGEyJiAxgMCUiMoDBlIjIAAZTIiIDGEyJiAxgMCUiMoDB\nlIjIAAZTIiIDGEyJiAxgMCUiMoDBlIjIAAZTIiIDGEyJiAxgMCUiMoDBlIjIAAZTIiIDGEyJiAxg\nMCUiMoDBlIjIAAZTIiIDGEyJiAxgMCUiMoDBlIjIAAZTIiIDGEyJiAxgMCUiMoDBlIjIAAZTIiID\nGEyJiAxgMCUiMoDBlIjIAAZTIiIDGEyJiAxgMCUiMoDBlIjIAAZTIiIDGEyJiAxgMCUiMoDBlIjI\nAAZTIiIDGEyJiAxgMCUiMoDBlIjIAAZTIiIDGEyJiAxgMCUiMoDBlIjIAAZTIiIDGEyJiAxgMCUi\nMiC1YCoit4nIURGZEZEh12PvE5F/rD1+RETmpdVOIqIo+lJ87WcB3ALgL5wHRaQPwIMAfl1VD4vI\n2wFUUmgfEVFkqQVTVf0RAIiI+6EbATyjqodrz3uzw00jIooti2OmPwtARWSviDwtIp/2e6KI3CUi\nB0Xk4IkTJzrYRCKiRm3tmYrI4wDe5fHQvar6DZ/T+gB8BMAHAEwC+I6IHFLV77ifqKoPAHgAAIaG\nhtRMq4mI4mtrMFXVjyY47VUA31XVfwYAEXkMwPsBNAVTIqKsyOJt/l4AK0WkvzYZ9QsAfphym4iI\nAqWZGvVxEXkVwIcA7BGRvQCgqqcA/BmAHwAYB/C0qu5Jq51ERFGkOZv/MICHfR57ENX0KCKirpDF\n23wioq7DYEpEZACDKRGRAQymREQGMJgSERnAYEpEZACDKRGRAQymREQGMJgSERnAYEpEZACDKRGR\nAQymREQGMJgSERnAYEpEZACDKRGRAQymREQGMJgSERnAYEpEZACDKRGRAQymREQGMJgSERnAYEpE\nZACDKRGRAQymREQGMJgSERnAYEpEZACDKRGRAQymREQGMJgSERnAYEpEZACDKRGRAQymREQGMJgS\nERnAYEpEZACDKRGRAQymREQGMJgSERnAYEpEZACDKRGRAQymREQGMJgSERnAYEpEZACDKRGRAQym\nREQGMJgSERnAYEpEZACDKRGRAQymREQGMJgSERnAYEpEZACDKRGRAQymREQGMJgSERnAYEpEZACD\nKRGRAQymREQGMJgSERnAYEpEZACDKRGRAQymREQGMJgSERnAYEpEZACDKRGRAQymREQGMJgSERnA\nYEpEZACDKRGRAQymREQGMJgSERnAYEpEZACDKRGRAQymREQGMJgSERnAYEpEZACDKRGRAQymREQG\nMJgSERnAYEpEZACDKRGRAQymREQGMJgSERmQWjAVkdtE5KiIzIjIkOP4PBF5SESOiMiPRGRrWm0k\nIooqzZ7pswBuAfBd1/FPABBVXQngWgC/JSLLOts0IqJ4+tJ6YVX9EQCIiPuhnwCYLyJ9AAoALgD4\nl862jogonsyNmarqtwCcBvA6gOMA/lRVT6bbKiKiYG3tmYrI4wDe5fHQvar6DZ9zNgHoB7AYwAIA\n3xORx1X1xx7PvQvAXQCwdOlSY+0mIoqrrcFUVT+a4LS1AB5W1QqAN0RkP4AhAE3BVFUfAPAAAAwN\nDWkrbSUiakXmbvMBPAfgBgAQkfkA1tSOERFlVpqpUR8XkVcBfAjAHhHZW3voLwDMEZFnAfwAwF+p\n6jNptZOIKIo0Z/MfBvCwx/FzAO7sfIuIiJJLLZgS0ew2OlbCyN5jeG2ijMXFAjavW47hwYG0m5UY\ngykRddzoWAlbv34E5co0AKA0UcbWrx8BgK4NqFmcgCKiHjey91g9kNrKlWmM7D2WUotax2BKRB33\n2kQ51vFuwGBKRB23uFiIdbwbMJgSUcdtXrccBSvfcKxg5bF53fKUWtQ6TkARUcfZk0yczSciatHw\n4EBXB0833uYTERnAYEpEZACDKRGRAZHGTEVkEMC/BvCvACwDMAPgZQDfA/A3qjrWrgYSEXWD0GAq\nIo8BOAngUQB/DuA1AALgcgAfBPAHIlJU1Y+1s6FERFkWpWf6SVV9w+P4i7U/u0RkkdlmEc0+WS78\nkeW2ZUXomKkdSEXk8+7H7GOqesJ804hmD7vwR2miDMXFwh+jY6W0m9a2to2OlbB2xz5cuWUP1u7Y\nl4nvtRWiGm23DxF5WlXf7zr2jKq+ry0ti2loaEgPHjyYdjOIElm7Yx9KHuvSB4oF7N9yQwotuqiV\ntvn1aN1Vo4DqCqgN1w7giedOdLwHLCKHVHWolWtEGTP9PwD8DoCrRMRZ8f5SAPtbeXEiqspy4Y+k\nbQsqs+dXNWrngeOwu3fdVpYvSmrU3wC4GcAjtb/tP9eq6qY2to1o1uh04Y84t9h+bciJBJ4XVGbP\nLxC775O7qSxflDHT06r6kqreAaCIi8F0SbsbRzRbdLLwR9wxUK+2AcC0auB5QT3aOL8kstA7jyJy\n0r6I/B6AnQDeWfvzoIj8brsaRjSbDA8O4P5bVmKgWICgOh55/y0r23J7G7cws922vEjTY0HnBfW2\nvQJ089WDr5M1cQqd/FsA16nqWaA+k/+PAP5LOxpGNNt0qvBHkjHQ4cEB3LNrPNZ5m9ct95xkck4q\nOSenrr9mEXYfKjUF+rPnpzA6Vsr8uGmcYCoAnN/lNPx/mRBRRi0uFjxn54N6gKNjJeREMO2R/eN3\nXliZPa9fHkNXLMT2R4/i1GSlfmyiXKlPRAVdL21xgulfAXhSROztmYcB/KX5JhFROwX1GL3YY6xe\ngTRsXNcdUO0hAWcAdKdPXZiaabpOuTKNe74yjj4RVGaq7cjabH/kPFMAEJH3A/hI7cvvZWlNPvNM\niaKLs6LJL880L4Iv3L4qMJB55ZNaecH8OX04Xa6g2G/hzLmpeoBMwkQubkfyTB0v9teq+usAnvY4\nRkRdJM74rFcgBaqz+XYCvl9g9prsqkwrJsrV23jn7XxSWZntj3Obv8L5hYjkAVxrtjlElDV5n7HS\nnAQn5g8PDnQk0GVltj/KCqitAP4IQEFE3kI1r1YAXADwQHubR0StiHM77/dcr0AKADMKbH/0qG+a\n1fDggO9klylZ2oQvStL+/ap6KYARVb1UVd9W+/vtqrq1A20kogTiJOcHPXcgoOfnd5tu90j9Ev5N\naGcubhKhwVRE3ikiXwLwPhH5ExF5WwfaRUQtCuo1ugUl8ifp+S0uFuo9Xfd1TRAA+7fckJlACkQb\nM/1vAA4B+M8Abqr9/ck2tomIQoTdvo+OlQJ7je7z/W7FX5soY3hwANseOVqfNHKz8oLK9MWhAEG1\nZ3vPrvGmtfamOIN1VnJOoywnvVxV71XVvar6uwAyUXKPaLaKcvseVBxEAdyza7zh/LClnNvWr/B9\nzvw5ffWhAMHFYiXtCqQFK4/rr1mUufqvkdbmi8gCEVkoIgsB5F1fE1EbuCs73Td6BGt37MPdu8ZD\nb9/DZtHdgc6eVXZ7/XQZ940ewfDggG9wnChXsHndcgwUC20LoDZ7nPSJ507Eqi/QCVFu8y9D9Tbf\n+V7buaYK4CrTjSLqBkGFj1u9/fRKOXrwwPHAc5wB9LKC5Xtb7kcBFKwcypWLK5BmFHjwwHG8eOKM\nb4oUgKbEfNOsnGDktosLBOLWCeiE0GCqqsuiXEhEVqjq0ZZbRJQRQUHRL7/y4MsnseupVxqWPG7+\n6mEA8ZY8Jpm4seuLDg8OwKPAUyTnPZZyAsD+F04GntfOQFosWNi2fkXD+5ekvkC7xUnaD/PXAN4f\n+iyiLhCWjO5bKf7J43B33iozim2PHI0VTJP0sOz6okDylUUtrOo0wu79DoT06OPWF+gEk8GUFaSo\nZwSlCgWt7PErdRH3ljtpsnu5Mo27fW6Bsy4vghfu/9VIzw2rSJWGyMWhI0j5dxqROWE1P9t9O9nO\nZPesmlaNvEupPQRTmigjJ4JSrSJV5mfziWabsD2Zrr9mkefjls8nakG/Fev1vSrvF/wu3kOipDg5\nU8MA1CfF0k6PMnmbf8HgtYiMijvDHjYm98RzJzzPu2SehTPnpxqS2K284LM3r/B8flgbnaXlrtyy\nJ/QavcA5nOIlaHIu7Nx2ilOC7zuq+kt+x1R1jenGEZkQNpnkJWxMzm8YYGKygi9uXN20HcfI3mO4\nZ9d4/Wv33vAAQttY7Lc8J5ZE/Mdqu1XQBFzY5Fxa6VFRqkbNA9AP4B0isgAXJ5reBiA7C2OJfIRN\nJvkJqvkZlJrjPC8sX9QOmvOsXGgb/QJm3EDqXv6ZRX5DGkHbp9jSSo+KMgjzW6gm7V9T+9v+8w0A\nX25f04jMSLKBXJgoWzOPjpXwqa8cDs3BLFemQ6svAfEzAvxs/MCS2GO4nTZZmWka+wzaPsWWZnpU\nlKT9/wTgP4nI76oqdyKlrtOOBO+wYYAoH/wociK4csseXFYwF/zCVlJlhbNXbv9iCno/w3JT2y3u\nHlAfBrAMjiCsqv/NfLPi4x5Q5MdrH6KClW9rLUy/fZP8FAsWzk/NtHUlUTcSIHSfKAHw4o6PtfY6\nnd4DCsDVAMZxcctnRbVEH1FmJU3wbmWNfdwhhJtWXY6hKxbWXy9sXHC2UISv5uqabUschgC8R+N0\nZYkyIs4GckCyDACnuCuYnnjuBD43fLGnPFvSoFqV9hJSpzhZwM8CeFe7GkKUJUEZAFH4TVD5cfdk\nixmfIMqCvEimti2J0zN9B4AfishTAM7bB1V1vfFWEaXMr1cZpbfptV1HXgQbrh3AE8+d8LyGojrO\naveyJgxsgdzL2j3mnUScYLqtXY0gyhq/2p35kNp2XpNdQHXJ4+5DJWy4dgC7D5U8J5rq5fqEhS68\nFAsWTpcrmShq4iVyMFXVv29nQ4iyxG/yJ2xSKGyp40NPvhJ4Db8Z69lOgKaaplkTZXfSf6j9/ZaI\n/Ivjz1si8i/tbyJR5/ltbxy07TEQPovPGfpkFMDmrx1OtSpUmChJ+x+p/X1p+5tDlA1Rig/fN3qk\n3tPMi+CO65YkrkNK4SrTmloRkyiirM0P2jTvvKqeNdgeokCd2t43LDf1vtEjDSuJplXx4IHjWHv1\nQrzx1rnMr33vVmnu8RQmypjpIfhvXtgn1QH5Laq602TDiNxayf1MEoSDclN3+izJ3P/CSVi55o9K\nTtLfEqQXZCVB30uU2/wrgx4XkUUA/h4Agym1VdLqT0mDcFAADoqLXpNIDKRm+BXlzoIot/lXqeqP\n/R5X1RMi8n+ZbRZRs6TVn6IEYXfgvP6aRdj1g1fqt+uliTI2fy3+LqMEzMkLLhga9th9qDoB5a4H\nm4WfSZQVUH8iIo+KyF0i8n4RuVxEFtf+fZeI7AHw4XY3lChsKxE/YUHYuQ2G4mLNUfe4Z2Vasf3R\n6m7m8+d4r2birpLNbv/AEmPXKlemsfPA8YafVZpblTiFBlNV/QSAuwEsAvDHAL4D4PHavxcB+F1V\nvaOdjSQCotUQ9RIWhOPsUW8X3fjjj69E3jU2ms8J7lyz1HPMdDYzXfLP3ceNs8y3nSIl7avqC6gG\nT6LUDA8O4ODLJxvSkTZcG17AJCzNKckMsd9sPwDs+sErsa9HjQRA/5w8zl6I9ksuC7P8cUrw3eJx\n+DSAI6r6hrkmEXkbHSth96FSPfHdXqI5dMXC0O1HAP80pzi5oUVHkWav2f61O/YxLSqB+a7AqQAu\nTM00bbEi8J78y8Isf5y1+f87gA8BeKL29S+imjZ1pYj8B1X9a8NtI2qQdDYfCE5z8uq5Wrnq2nzn\nLLyVE2xbH7zLaBZ6SN2o2D8HZy80vneVGUWxYGH+3L6GiUF3bYOslOGLE0z7APycqv4TAIjIT6Fa\nGPo6AN8FwGBKbdWOvZyA4Fv2uLmpfjuIkr9iwfL9GZ4uVzD+2RsbjjmLaGdpNj9OMF1iB9KaN2rH\nTooI//dQ27Wyl9Od//c/Yv8LJ+tfr716IXb+uw/Vv/bruQZ9SL2Wk3LpfXxBGwW6f7adWgGXROQ9\noETkzwEsBfDV2qFbAbwCYDOA/66q17elhRFxD6jel3QvJ3cgtbkDatQ2jOw9xvX3HWDlBCO3rfLd\nNhu4OIba6mZ6Hd0DCsC/B3ALgI/Uvv7/AOyubWOSaiCl2SHpXk5egTTouB+/WqXUHpfM62v42XqN\nmdtdwbjbyrRDnHqmWivHdwHV7+Ep7gdFnRZ3LyeTouaj+s04Uzzu3QbCxsajTka2S+Q9oETkdgBP\noXp7fzuAJ0Xk1nY1jChrok50eQXSsAr91Mw9XhplbDzNbIo4t/n3AviAnVNaK3DyOICvtaNhRKbM\n7cvh/NSM5/Eg7smOywpW4GRJEBaFjsfKS1O6k1cKm1ua+aZxdifNuZLz34x5PlEqPr/hfXCv8BRU\nJ6+u3LIHa3fsa1rbPTpWwuavHW5YA/7W+SkuFTUkKHCIABs/sKTpdn14cAD337KyvtuB+yeRdr5p\nnJ7pt0RkL4CHal9vBPCY+SYR+UtalxS4OHFV7Ldw5txUvZfpNXmx/dGjTSuZpmcU8+bk8c7+OXht\noox5VrXHy/J68WxasxRAtSas11unCt+Vbc4x86ylSUVOjQIAEdkAYG3ty++p6sNtaVUCTI3qfUlT\no9zW7tjnmdo0UCxg/5YbAADLtuzxPf+lHR+LdL3ZqhgyHDJQLGDywlTo4oZiwWpK2G8XE6lRsW7T\nVXW3qv5+7U9mAinNDkHLSeMwvZKKS0gbzZ/bF7jxYGmiHGmV2ES5konSelFF2Z3UvSspdyelVPj1\n/qL0CkfHSli7Yx+u3LIHOZ+ZdefkhbOgiZPX8SwU2ciS1ybKnuUSk8hCab2ootQzvVRV3+bx51JV\nfVsnGkkE+KcXhaUduYs/e82suycvtq1f4fnhmChXmiasTAWOXrG4WKhPFvVbrc1Rd1Ovn7Px1DX8\n0ovC0o78ku3zIhBUx/C8xl3zee8g7a7u7pxltq+3ac3SWRtgSxNl3Dd6BMODAyhXmlPSbMWCVX/P\n/JIkuqnXH2c2nyhVAz6FToLG5wD/3s2MKl50TSbZRvYeC6xLWq5MY/ujRwNnk4euWIhPfeXwrMwx\nffDAcTz54zcDV4LdtOpyfG54JQD/ycUslNaLij1T6hpJty25zGf8sxBwCxplHPbUZKUhD/XuXeNY\ntmUPrt76WL1ndsd15vY/6jb/642zgY/vPlQK7N3HzdJIG3um1DWSFjrxG1KdrMxgdKzkeX5eJHGP\nclq1vu/RE8+dSHSN2cC9lj7NugsmMJhSV0nygXMXzHDyK4xh4tbc9EZyvaibJpjC8Dafel7QJIbf\nhzlsHDYqLj4N1k0TTGEYTKnnbV633Deo+X2YvcZnkwTG2Tf1FF23TTCFYTClnjc8OIA71yyNVRjD\na0LkTo90J/Y8k1nQb+H+W6oz+fZiCq+CM90ktTFTERkBcDOqxaZfAPAbqjpRe2wrqruhTgP4PVXd\nm1Y7qbv4Fb/43PDK2BuxeY3POq9xWcGCSHVMtmDlMOmTU1mwcoH5lrNR/5xq6HGmQ2WhWn4rYhU6\nMfrCIjcC2KeqUyLyeQBQ1T8UkfegWpnqgwAWo1oz9WdVNbDEOQudzB5+AdNUIZSobWjaHjovgFa3\nKK4fy0nD11Ql8N8g0VlwpmPt6XShE5NU9duqOlX78gCAd9f+/WsA/lZVz6vqiwCeRzWwEjUtDXWu\nRjJVCCUKr9eqTGtTFgADqbfFxULbtu5OS1bGTH8TwDdr/x5AdddT26u1Y01E5C4ROSgiB0+cYD7f\nbBAUMDv54fRfVWX8pXqOPVbtN/nXrTP8bQ2mIvK4iDzr8efXHM+5F8AUgJ1xr6+qD6jqkKoOLVq0\nyGTTKaOCAmYnP5zd+oFPi1cdhKQr2rKqrRNQqvrRoMdF5JMAbgLwS46dTksAnGvw3l07RuQ7zmaP\nnZpc3x1Uyf36axYxKT+GO65bUl+Hb0u6oi2rUrvNF5FfAfBpAOtVddLx0CMAPiEic0XkSgA/g+qu\nqESBvRmT67uDxmYBLhONaza8X2nO5j8PYC6qG/MBwAFV/e3aY/eiOo46BeBuVf2m91Uu4mz+7JFk\n75+wc9yPnzx73jOdqViwMH9u36zfpmSgNoEUNXoI0FShq5PZF6HtMzCbn1owNY3BlPyEfWjvGz3i\nu7kbedu0ZikeevKVyDUMvNKdouzFFcTkhnomgikLnVDP88sA+NRXDuPuXeMptaq7/c2B44i6DMHK\ni+e4dSvZF+5fkFlI+M9KahRR2/h9OGdj0WZT4qznmj+nzzPAtZJ90cmc4qgYTKnOuelct6+TdmIa\nU7pO+2z73EpqVBYT/hlMCUD47HU344Z36fL7ZdZK9kUWE/45ZkoAgm+bujXvz+bOZ8yFVNEXAB++\neiF++Ppb9f3diwULN626HLueeqWlJaIL+i2cOT8VuL9UN8gh2q1+WE8zaXV90znFJjCYEoBs3jZ5\nSTqD6/zQes3uC6q1R+3tSl56s4zP3rzCs2rUtkeOYsLn1jWIlRN89uYVOPjyyVgz4Vmz9uqFuG1o\nKUb2HgtMERtoYxJ+FhP+mRpFAFpPU+kEk3mJ7qB8/TWLsPtQKfK1V3zmWzh7wb+QWbFWns/u2Ras\nHOZZ+frX7bBpzVJ8bnhlw/dW7Ldw5txUS71p+xeNzX5fAPhmQ3jllWYZU6PImCzeNrmZHIpw316u\n3bEv1rWtfA7VcrveTpcr9WBi/xJoZyAFLq4ycvfCtz96tKXXdofhKLPms3HSj8GUAGTjtinsFr7V\nvER3T/SJ507Uv/a7XfW7tt8Mtc0ZTLx+CbSDu61ePXlTwlaAZemXcKcwmFJdmlvtRknCDipyEvfa\nziIlpYly061s2LWDArC7R9+pcWd3WzsVxN2KBavrJy2TYGoUZUKUJOykeYlRgorfiOL113iXdgxK\nt5pnNX6sOnHLa+UEkxemGnKE2xXEg/a9EgDb1q9oy+tmHYMpZUKUW/goeYleCw9aCSpe1Y78qvrb\nTk1WcPeucaze/m2MjpWwed1y5HPt23qvWLCA2mSXM0f4soLle07Bav7oW3kJ3SBwQb8VWMNA0Z37\nN5nA23zKhKi38EFDEX5DBcV+K/EETCvjkBPlCrZ+/Qjuv2UlZtpUgn+g9v64U7XKlWmcm5qGlZeG\nnFZnhoLXGDWAwO+vf04f+uf4V80amIUTTzamRlEmmEh78kvvKhYsnD2fLD0oL4IZ1XqwCcut9GLl\ngLQ2J7Vygkvm9WFisoL+OXnfdC47dWtistLSLx8RQLX6nlemZ+qvVyxY2La+OW83K5gaRT0jajZB\n0Iy/3+38RLlS3TnUISfAZYXwoGEn1tu93CQTOmnu8lyZUfTP6cPH3nd54M4A5cpMvX5rS2lUtd9X\n7p7yRLmCzV89DKB3hwHYM6VMcieeq1Y/kH4J5MODA74907zP8lG/40GSnGPSQLGA10+XY23cJ0Do\nEtpOydIiEKeu3uqZyI+76MqpyUq9pxOUQO432+8XRJIEl2lVWG2cTAqzed3y2DugLi4WMhFIgewt\nTzaJwZQyJ25+pP0B9ZvtjzspUrDyWNDvPRM+UCzgknnJRsfy0noQjluv004dM/HaJvTyyiiOmVLL\nTG4fAcTvvTg/oH6z/Zu/ejjSBNSAY1bbfY6Vk5Z2JY3aO/RbQACErzxyWlAbHrln1zgKVg6TlXR7\np1bOu+J+r2DPlFrSjjqocXovkesHROiY2eN59WDsOmcGwK6nXgm9TtiHyrmHfL9HvmcrIc/KAV/a\nuBpf2rjmX25bAAAgAElEQVQa5yozmChXc08nY8yCDRQL1dzViERQ78nbPeBiwcL8OReHXIoFCyO3\nrerZySeAPVNqUTvqoHoVXXGye25RS7yN7D0WWj/UHZS9zpme0YDSJhfNrU2K3bNr3DMwzqjWi6Bc\nuWVPhCtG85KjSpNX4ZYo7F8oXqlqVl4ARUNv3St9zXmn4uzpj+w9hnt2jWeiXF47MJhSS9pRB9Wd\nJmXP5p8uVxJ9EIPaIkCsoipR2L9MoixECFrjH4d7XDhp+53jz0BzqprXMXcgdS+c2PzVw4Cg/ssp\nC5vftQODKbUkafGRMH5jn3avJ04Px6+NQWk6rQa51ybK+OLG1aFlDcN64VF4DXUkbb9zCarfzyDo\n/fa6U/Eaq+6VXRycOGZKidhr4O2KS07tqoOadHw2SYEUvwInUT8wORHcs2scc/tyWNBv+dYScGcg\nxJ1zFwAbrm0Oekn3vTp7Yaql8e44PeJeS5Niz5Ric9/KKeKPYyaRdHx2eHCgYauQvEhTAHJnJExe\nmPK81mWOBQRB7Jn7iXIFBSuPL25c3VCw2X2rbPeQ446hKryLsdiv9amvHPZdsJDPARdc48KVacWn\nvpJ8pVKcHnGvpUkxmGac6bQjE7yCmh1I27m6Jen47OhYCbsPlepBZVoVuw+VMHTFQgBoqkQfFAxO\nTVZi9/icAf++0SPYeeB4fWLKPX6Y5Pbc7/sfHhzAPT7bikyrYtpnZGFaNfGYptewhZWThjFTIHu7\nOJjAYJphUQompyGtzfeSjs/69Wi3PXIU56dmYo9XJhnffG2ijNGxUkMgdV7PDrZJxlCLHgsMRsdK\n+KOvP+ObZhW2LNbZptGxUsMmggv6Lc/NBoHoE1fXX7Oo52b3GUwzLKvbL7dr0ilM0n2q/Hp6SXYY\nTeqygoWRvcdCk/Htn2ucHVDPnKuOczqHEX7/K+OBy06jLCAo1X4BuBcvnJqsYPPX/IcCwiaustpJ\naBUnoDIsq9svJ61436ooxaG9ZGEppUj4z80uZj08OIDxz96IL21cHantlRltWGY6svdY6Pr9KNfN\ni1TzbT0uVplufM04ouyq0I3YM82wVnqA7RxrTXPzvST7VJku8lGw8rFv9U9NVkJvrd09tKAxTzdn\noI7yyzbKezKtGnitVnNZTV0vKxhMMyzpbW0nbqO8glrWJsvs9vgpFqzYY6Z2toLfLHmQKM93D+NE\nnZCKuxBgoFjA2fNTgUMJ9kIAv2slHdZJa5io3Xibn2FJb2vTuI1qxxp9U+3xUrDy2LZ+Be6/ZaXn\n+ni/c+xfEF+4fVWiPM4onD20KPmiXgsBgqoE2s/ftn6Fb16r1K6zed1yz5KDVj550ZK0honajcWh\ne9CVW/Z4TnQIUF8TbppfYWa/dKl292L92gNUxwLvuG4JPje8sn7svtEj9TxUp5wAM+qdP+v+HiYm\nL/huCxKHHbrcM+GliXJ9qMD+u1iwIAJMTDYutbVn890FTtzfhztVy6lYsHDTqsvx8NOlhu9r/pw8\n/vjj0beT8ZK1uxgTxaEZTHtQ3MBmQpwAbmK/p6TtCXq9VtsV9pq2gRi5pPbrA80b3fnlb8Z9H92p\nT2GsnPRcBShW2idPadxG+Y13eR3vxDBE2Pib1+u12q4oY35xizXbr++35t1d2SrJ+zg8OID5c6NP\nn7izB6iKwbQHJR1rbUWcAN6J2dwoY43u12u1XWGv6VxHf8d1SyJd0379dq95j3tOt8+8twNn83tU\nkhSiVl8PiJYu1c7ZXPdGfHP7cr63rzkRXLllT72tfu1SVIdOwsb1nO+B33UeevIVDF2xsD5eG6Vq\n/+KQWXW/58cRdxmr8zWyNv6ZFo6ZUl2nPhRJxybD2ud33Q3XDmD3oVJgClSU55kaP3Vex6vNXs8F\n2jdmClTfu81fOxxaRNt+XXvMtBPj353AMVMyppOpTUmGIaK0z2/M84nnTmDDtQOB6UL284I24CtX\npnH3rnGs3v7twPdldKyEXMCYqHNc0/1eLOi3qjP0aHxfvN6zkdtWYeTWVUaGc4YHBzBy66qGjQSL\nBQub1ixtOuacfOrV1UxJsGdKANLJAIhqdKwUWEruC7dXP9xhvcGw5Hxn5kHYzLzfjHZYT9PrtbpZ\nGml47WCiZ8oxUwKQ3SV+dnDyWz3kLBfnN+6XF4m0ykkBLItYT9Se0XYH06jbVF9WsLB2xz7PIQt7\nOKM0Ua7nudri5L0ue3sB3//xSdhv3Zy8IJ8TlB25p3Pygsq0egbEuX05FKx8ooIw3b6aKQkGUwKQ\n3SV+UYKTfVvptw2z6bX5Nq/3K8ovHysnOHvh4lJO53JfoHFs1F1jxP7avUTYawmxu30XphVwjYm6\ni0M7nZ+awfmp6LuaOr+/bl/NlASDKQFIXgcgiIkJrag949cmyp4V59vJzhV1fp85n2ImeRHMqNar\n+DuLUQON44xRawU41/FH7RF3wiXz+rpq8skUBlMCYLYSlNeKmtJEGXfvGse2R45i23rvwsJeiv1W\nU+DxsrhY6PiQxLRqU4/QK5C6Z7f9tiZpJT807eEYp4kIP69exGBKdSZyU8MmYCbKlcgVrEbHSjhz\nznsvJie7B+2X39kuA8WCb4/Q2RN1/1IKG1JJku9pastoE9IeGkoLU6PIqDhjnFGu5VWYuGDlPNOB\n/CoctYNdVcmvRzijihd3fAz7t9wQaedQ+xdCnF1FncMwSXcjNa0Xqj8lxZ4ptcw5Zhh1qqc0Ufad\nzbb5BapzlRnPdK3hwYGmzfFasfbqhdj/wknPx7T2en694aDeWZTdUqP0st35qu5rXrWoH8+fOFuf\nzc/XMgGcP6M5eQmchIrD6/uYTRhMqSVR8yq92MHCr3h1kgwDU4EUAG4bWoqX3myeFQcuFk5OMnEX\ntFuqPdRivw9BpQSBi+/dwZdPNl3z1VPn8MXbVweuVIqyOszKC0ZuDc+pdX8fsw1v86klpmaRvW79\n41a/Mr1aa2TvsdA2JFnNFWfVUJTb93JlGg89+UrgNf1e0+s8N7/9nrj6qRF7ptSSoFlkASLf9gPe\nt/7337IycoaB6Q9xaaJcDxh2MWavZPm4E3dxFki4syz83k+/XNqw2f6oObhe52d1oUdaGEypJX63\n4vYy1NXbvx1rBY3z1v+eXePQ2rW+uHF1aMAy/SEWR3umVRu2LWlF3OGLOLf9ftcMWh0WJaB6tS2r\nCz3Swtt8aknYbXAruyzbH/GwoiujYyWs3bEvVi84jFev2tQtbCvFuzevWx64b5PfNf1e847rloQO\nI/jt99SrezklxZ4ptSRsZtpUAne5Mo3tjx5tuuUHmkvTtUIQnLNpovfbygKJ4cEB3O2z/bPdi/e6\nZtBrDl2xsKEG7LnKdH39/oJ+C5+92XuRRZpbfmcRq0ZRS8LqWca9LY2jYOUxz8rFmsEvWLmGQh9O\neRG8cP+vAsh2Fa0st61bcUM9BwbTdIR9sFtJneq0/lqgXVws4PprFjWlDLmXhfrVHnAe75+Tx+SF\naSi8d0V1ilrLwO89LVg5zLPyODVZaRimcPYunbuwigCFvovfc9CupWE7knZ7tX0GUwcG03REqWfp\nLCmXVe5Sd3YO5hPPnfAtk5c0b3PTmqVNATVuxfrRsVKsBQpWXvDBZQt8FyE4X+/gyyc9q2/lc4Iv\nRKzh2m3V9hlMHRhM0+E3W+9cm379NYvwxHMnMh1MvQTdNvv1yKPMjjuHE8Kul6QNrRgoFvCT0+d8\nvwev9vTCsAOLQ1OqRsdKOHvBuxCJ/WEsTZQjbRqXRUGTTX5BLEqakddzkuRstiOfM2xJMPNN/TGY\ndhlTY1MmrjOy91ikDdi6VVC+ZNT8TL9zvV4rbs5mOypFLQ7pmTLf1B/zTLuIqU3v/K5z3+gRrN2x\nD1du2YO1O/aFXreXex5h+ZJBgTQstfaO65Y0HfNbNjp5Ycr353D9NYtCX8uWzwnWXr0w8Dn29+zV\nPvsazDf1x2DaRUythfa7zs4Dx2MF6l7teRQLVujkid8OpkB1Ft0OcvPn5Ov/zot4Tj4BF9f4FwtW\nw/FTkxXPn4NdLMUd0ufP8U7Az6FauGXTmqX1nrFINYPBXVPgc8MrsWnN0oZAPX9O3nPyydl2E7uk\ndjNOQHWRoB0z7WTzKLfrYTtvOgVNIoyOlXwTyOMQAbL039C546mfKClfSSZgok7mBD0P8B7T7aYJ\noU4zMQHFnmkXCeoJxrntj9OjDBqTM9XzyFIgBS7ueBr0Pjp7Y35a2YYk7HjQ8zghlA4G0y4StRyb\n322/vYa9NFGOPNZmnzfbRBk+GR4cwP4tN/gG1CTDIH7nuI8HPS/qNcgsBtMu4h6b8uPVA3FOOgGN\n43ph/Mbs1u7YF/EK3SlqT87kBEzUa8Xd+mQ2Tgh1GlOjukyUcmxePRCvSSf77joszce5pTDQWnX9\nbhK1J2ey4EfUa0V5Xjcv7+xGnIDqYnGW8cWZdPLiXB7azuIlndBv5TBZK3YiAnz4qoV4+vjprl4O\nSa3hCqhZLk6PqNUEb2cvrZsnMooFq2HVlirw1EunsPEDS3zX4QNoKBASVrDErZUFEu5z7aW5fvUC\nvEoUOsvrqcJz+a+zNkGxYGHb+hUNRVtKE+XA3QaIPdNZo5Vbc3cvbfA/fNvoxnWdIgD65+Rx9kLz\ne7Cg38LYZ270PO++0SOeS2L9ckadWikCEuVnZl8LaK7rauUEECRapWblBBs/uMS3aEuv9dyZGkWR\nRUnlWdBvYUG/1XTM/aHp1t+/CngGUiB4V9OHnnwl1nGnVhZaRNms0L6W13MrM5p4uW9lRgM325vN\nG+f54W1+RnSiHqQ9eeXXW3rP5Zfi+64SbWfOTWH7o0dxz67xertOx9jTKUsGYgx1OH8ecTexc2ol\n5zPqcEq7hl3Cvr9uHu5pB/ZMM8DUmvuohgcHsOHagfqywrwI3r/0Mnz/hZNNgaMyozg1Wam36+5d\n4+jrwv81dmqQe7mmTXAxn9b98/DjVbDErZWcz6jZBEG5pa0I+/6Yt9qoCz8WvafT+4/b67rtnse0\nqmcg9VOZ6a7/OM6him3rV1THEl0UF/Npo9xeA94FS9xa3TwvbJFGUG6plRNY+WQ7Glo5Cdxsj3mr\nzXibnwGdXv4XlHMalfcuStl0arKCu3eN4/e/Mo4PXbUQ8+f2ec5o27/Awt73OLP5rW6e5z43aDbf\n63Wcx5LM5tub7XE2Pxxn8zOg05XKW8057WVBu5OyUEjv4mx+jwi6Fbxv9Aiu3voYlm3Zg6u3Pob7\nRo9Evq695NNdn5RjXf7sHh2XY1JcvM3PAL9bQffGZtOq9a/DbjHdu0uWJsq4Z9c4Dr58EpvXLW+a\nzXfuZjlb2QGT+8FTErzNz7AoO396GR0r4Z5d477nfnHjagCNwcLE8tB2BmR7g77LChYq0zO++aIi\n1TG/uIsKOA44u3E5aY/zC0zu4+4c1bPnpwLPHdl7DPu33NAQOEysty/2xw9iUbhX26zdsQ9nL3i3\nVbWaG5v02kRJccy0y3nlqHrN1jqVJspNOaxR0nDCTExWAldYJeUOdmGz7ZWZ6P1jBlIyhcG0y0XN\niXRzLwpw10r120soSLHfMhKUG65ZsJqC3WU+ifdxDRQLnoHUb+KOKEhqt/kiMgLgZgAXALwA4DdU\ndUJEfhnADgBzao9tVtXerkLsY9Oapb4FNmxJc1HdNUrdzlXiZ5KqVoPyVw8ex37XstSkJsoV/PQf\nPYapWm9TUM2JNKE0UcbaHfvqY6WjYyVse+RoQ8/eXo1ma/ekVJa28qZ4UpuAEpEbAexT1SkR+TwA\nqOofisgggH9S1ddE5L0A9qpq6P+CXpyAAsJLv7U61vlSbSLLRMFnAfDhqxcaC6SdUrDy2HDtgG+F\nJKC6iupcZaatNU9bqTDVjuvMJiYmoDIxmy8iHwdwq6re6TouAN4EcLmqng+6Rq8G0zBeH5w4s+p2\nkvrkhamWJ48KVg7lBD3aLAjbbcCPyUR+U4s3Or0IpBf0UtL+bwL4psfxDQCeDguks5nXnuVxQoI9\nadVqILVy0rWBFIhWAcqLySW/ppYVc3fSdLR1zFREHgfwLo+H7lXVb9Secy+AKQA7XeeuAPB5AN4V\ne6vPuQvAXQCwdOlSv6f1POe+UEDntxUpFqzQDIKsC+qZWjnBlKpnHVeTq8n88n3jvoap61A8be2Z\nqupHVfW9Hn/sQPpJADcBuFMd4w0i8m4ADwP4N6r6QsD1H1DVIVUdWrRoUTu/la5iekY9TNqBtGDl\n8aWNq5sKWzsFVU8qWHnfCkn9Vg4Q74LYppeYmlrGyuWw6UjtNl9EfgXApwGsV9VJx/EigD0Atqjq\n/rTalxVJ0nS86pWuvXphfSggSh3ObiFSzUy4e9c4zgVMns2f09cQbO13YKBYwP23rMTnhlc2DZd8\naeNqLJg/17NafV7E+ISO15BNktcwdR2KJ83Z/OcBzEV1ggkADqjqb4vIfQC2AvhfjqffqKpvBF2v\nFyegks7Khp3nt6dRN8kJICKYjpGgX7Dysd/LpEt6qbt09QSUqv60qi5R1dW1P79dO/45VZ3vOL46\nLJD2qqRFo8POe+K5E4naM39OHnH7tJvWLG3LqigFYgXSvEii97KVSvk0u2RlNp88JJ2VDTsvyayu\nAPjjj6/E4phDBXueeb0tk2FxbqisnP/kUth7wfFHiorBNMOS9orCzkvSq7qsYDXUAIiaSpTGltAL\n+q36eGGxYAW2Ney94PgjRcWqURnmVXc0Sq8o7Dyvx62c4JJ5fb7BrzI909LqqE6xcoLP3ryiocKU\nX7aBAJF6mO7UMyIv7JlmWNJeUdh5Xo+P3LYKY5+50Xd8069+aNqcM/TFgoWR21ZFrjClAIMkGcOe\nacYl7RWFnef3eDetksmLYOwzvms6APgnsANoy8QYzV4Mph0SVrDET9LqP0nOu2/0iO9S1GLBwuly\npe1bm1h5wfw53ruHuq25akF9tZd750wA2P7o0cAxW3fVKJNYtWn2yUShExOynGfql9e5ac3SwIDa\nrjzTOG1sVdytTL60cXXT3lde1+zLAX6lAKycYAbRU6eyWv2JOqer80xnk4eefCXWcVu78ky97GxT\nEn+x38KmNUsjL2/d/uhR30CaF8GXNq7GPCvvG0iBaqV9v0DqldIV5T2NI+nPjbobg2kH+KXmhKUX\ntSvP1Eu77k9OTVaw+1AJG66N1iMLui2fVk28s4DzGl6yWP2JuguDaQf4JbiHJb6bzjM1td1HXOXK\nNPY883rLNQHyIi0HJL82mK7+1O7XoOxhMO2AO65bEuu4bfO65bBce3RYOYmUZ+r1g50oV7B6+7cx\nOlZqKqAyJ6CqkgmnJiuJa4ba7rhuSUsBycqJb3Wos+enjO31xFVTsxNn8zvAnmRKMpvftBg+Yszz\nG1KcKFew+auHAUG9GlJpogwrL23d974Vzvcr6URZsWBh2/pqMv/QFQubZvonypX6Xk+tThLZ53M2\nf3bhbH6GJd1+Imlx6AX9Fvrn9NUDwPXXLMKeZ16vBx3xqevZCnclJ/dj7hnwJN+b1/vFrT3IycRs\nPnumGWZ6AirMxGSlKQne2Xs2semek50TavfgLitYEKm2w683l+R78zqHk0RkGoNphiXdfiJo1U/Y\neV6cCejFfgtz+3KYKFdaHhawg6UzYDpfy04lcj6e5Hvz+r64tQeZxgmoDEs6kRG2bYmVk6ZtPPyu\na/dG7WpRpyYrOF2uoN/KtRRI1169sKnXed/oEdyza7z+Wvae9c6Joc3rlvsOGy/otyK/X5wkItMY\nTDPMVKGTBf0WigWroajJyK2rIl3XK69TAUzG2Ik0h2plfKA6mbRpzVLs/HcfanjO6FgJOw8cbwrQ\n7mT34cEB3LmmefNEu1pU1PeLpfXINE5AUSC/bTuiECDyTHbQxJJ7i5DRsRI2f+1ww95MVl4wcusq\nBkNKhBNQXaSTBUuiCCq84nzNXMAWyEHyIvjC7dXgZue0OrMEnnjuRP3rZW8PHgctunYd3f7o0aZN\n7irT1dVRDKaUFgbTDnDPgttjgUBwTmPS88K4czWnVfHggeN48cQZ/PD1txryL5Mm2k+rYuvXj+Dg\nyyex+1Cp4XtwvnZpohw6oXTmXDWh3g7MfktOORNPaeKYaQd0smBJFH4FVva/cNLoNiPlyjQeevKV\nllOpKjNa/56DvnfOxFOa2DPtgE4WLImi1WWdabxWlM0A3TPxrClKncSeaQe0a2O8pFotOJKGsM0A\niwWrKV/VmdLllWZFZBKDaQeYzBc1kQsZVmAla9ybAXq9J9vWr2g4xpqi1Gm8zU8g7u1j0sIXcc+L\n2i6vwitrrlqAp4+fjjy+KQD65+QTb7RXLFj1rUnm9uVwfso7bzUv0rQZIBD+nnC5KHUa80xjyuqW\nFCbaFTclqmDlMc/KJZq0svLSlN7kxZ1jGhULmVAc3LYkBVm9fTTRruHBAezfcgNe3PExfOH2VaFb\njZQr01BF5C1JnKIEUiD5+DCXi1Kn8TY/pqzePppul/t22i/0TZQrWOAofpITwN5+qWDlUI6x7NSt\nYOVx/TWLGhL+o87IDw8O4ODLJxuGMjZcm2zbbKIo2DONKatbUrSjXc6eatAe86cmKzg/NYNNa5Zi\nbt/F3mC5MhO1lnWd/fyBYgEbrh3A7kOlRDPyo2Ml7D5Uqg9VTKti96ESZ/OpbRhMY8rq7WO72xVW\niapcmcbOA8c9i6K4A6pX1Srn8+29nrwS/qMOXWR1OIZ6F2/zY8rqlhTtbpfz+n7LP/2GAhTVnqaz\nXQBw965xz+c7e5NeogxdZHU4hnoXg2kC7oLGSZleoWOqXWHXj7t1iN8MelBgDhJl6ILFn6nTGExT\n0q4iJs7rBwVqr8eBaD3bzeuWx9q+5PprFnkej3sdIPrQhde1szAcQ72LwTQlQWN6rQbTsEDt9bjX\njqV+wd1rSGHywpRvvunuQyUMXdFcWd99nbDc1n4rhz+JmDeb1eEY6l0Mpilp55heWKD2erwy0xzE\nnBM2XkHJay28Vy8z6JeE8zphG/bFqe7vvjZRu3E2PyXtTLEKC9RxArbdQw1LT7K3AYnbpjjXAMBi\nJZRZDKYpaWcqU1igjhOw8yKRU4yGBwd881Gjvubw4EBgVSumN1FWMZimpF0buo2OlXD2/FTT8bDK\nS347lsZNTzLxS2LNVQsCH2d6E2URx0xTZHpMz2/McUG/hc/evCK08pLXMb/0Jb+epomJn5feDA6W\nTG+iLGIw7SFeE0sA0D+nz3Mm3W9CyC1uilGrvySCep5Mb6Ks4m1+D2lHhkAa+8v79TzdtU2JsoTB\ntIe0I0MgjX2U/MZd7a2jibKIwbSHmM4QSGsfpTR6w0St4phpDzG96qedq7TCMOGeug2DaY8xGYRY\neYkoOgbTHtBK0ZIgcSovtasNUaQxrkvkxg31upxXbqmVF0Ab19sn2fQv6iZ9nm3ISUPhlKRtMNVG\noiDcUI+8i5ZMa1PhkiTLMKNOBPkVTnFvmteOpaCsqE9Zwdv8Lhdn/DLJWGeUMdh2tyHJ9TiuS53G\nnmmXi5ND2q5lmGm2IasbHNLsw2Da5TyLluSlOmbp0M5lmHEKp5huQ1Y3OKTZh7f5XS5O0ZJ2Tcik\n2QZW1Kes4Gw+Ec16nM0nIsoIBlMiIgMYTImIDGAwJSIygMGUiMgABlMiIgOYZ0qZwepP1M0YTCkT\n3NWf7Kr+gPcmf0RZw9t8ygRWf6Jux2BKmcDqT9TtGEwpE1j9ibodgyllAqs/UbfjBBRlAqs/Ubdj\nMKXM4PbO1M14m09EZACDKRGRAQymREQGMJgSERnAYEpEZACDKRGRAQymREQGMJgSERnAYEpEZACD\nKRGRAQymREQGMJgSERnAYEpEZACDKRGRAQymREQGpBZMRWRERJ4TkWdE5GERKboeXyoiZ0TkD9Jq\nIxFRVGn2TP8OwHtV9X0A/ieAra7H/wzANzveKiKiBFILpqr6bVWdqn15AMC77cdEZBjAiwCOptE2\nIqK4sjJm+puo9UJF5BIAfwhge6otIiKKoa17QInI4wDe5fHQvar6jdpz7gUwBWBn7bFtAL6oqmdE\nJOz6dwG4CwCWLl1qqNVk0uhYiZvk0awgqprei4t8EsBvAfglVZ2sHfsegCW1pxQBzAD4jKp+Oeha\nQ0NDevDgwTa2luIaHSth69ePoFyZrh8rWHncf8tKBlTKFBE5pKpDrVwjtd1JReRXAHwawC/YgRQA\nVPXnHc/ZBuBMWCClbBrZe6whkAJAuTKNkb3HGEyp56Q5ZvplAJcC+DsRGReR/5piW6gNXpsoxzpO\n1M1S65mq6k9HeM62DjSF2mRxsYCSR+BcXCyk0Bqi9srKbD71oM3rlqNg5RuOFaw8Nq9bnlKLiNon\ntZ4p9T57XJSz+TQbMJhSWw0PDjB40qzA23wiIgMYTImIDGAwJSIygMGUiMgABlMiIgMYTImIDGAw\nJSIygMGUiMgABlMiIgMYTImIDGAwJSIygMGUiMgABlMiIgMYTImIDGAwJSIygMGUiMgABlMiIgMY\nTImIDGAwJSIygMGUiMgABlMiIgMYTImIDGAwJSIygMGUiMgABlMiIgMYTImIDGAwJSIygMGUiMgA\nBlMiIgMYTImIDGAwJSIygMGUiMgABlMiIgMYTImIDGAwJSIygMGUiMgABlMiIgNEVdNugxEi8haA\nY2m3I4F3APjntBsRE9vcGd3YZqA7271cVS9t5QJ9plqSAcdUdSjtRsQlIge7rd1sc2d0Y5uB7my3\niBxs9Rq8zSciMoDBlIjIgF4Kpg+k3YCEurHdbHNndGObge5sd8tt7pkJKCKiNPVSz5SIKDVdH0xF\nZEREnhORZ0TkYREp1o7/sogcEpEjtb9vSLuttoA2v11EnhCRMyLy5bTb6eTX5tpjW0XkeRE5JiLr\n0mynm4jcJiJHRWRGRIYcx+eIyF/V/n8cFpFfTLGZDQLaPE9EHqq1+UcisjXNdjoFtPlOERl3/JkR\nkQIOCg8AAAZeSURBVNVpttXJr921x94nIv9Ye/yIiMwLvJiqdvUfADcC6Kv9+/MAPl/79yCAxbV/\nvxdAKe22RmjzfAAfAfDbAL6cdjsjtvk9AA4DmAvgSgAvAMin3V5Hu38OwHIA/wPAkOP4vwfwV7V/\nvxPAIQC5tNsb0uZPAvjb2r/7AbwEYFna7Q1qs+s5KwG8kHZbI77XfQCeAbCq9vXbw/5fd33PVFW/\nrapTtS8PAHh37fiYqr5WO34UQEFE5qbRRreANp9V1X8AcC61xvnwazOAX0P1A35eVV8E8DyAD6bR\nRi+q+iNV9VrM8R4A+2rPeQPABIBM5EYGtPknAOaLSB+AAoALAP6lo43zEdBmpzsA/G0n2hNVQLtv\nBPCMqh6uPe9NVZ0OulbXB1OX3wTwTY/jGwA8rarnO9yeKPzanGXONg8AeMXx2Ku1Y1l3GMB6EekT\nkSsBXAtgScptCqSq3wJwGsDrAI4D+FNVPZluq2LZCOChtBsR0c8CUBHZKyJPi8inw07oihVQIvI4\ngHd5PHSvqn6j9px7AUwB2Ok6dwWqt6U3trudrtdN3Oa0dGObgWjt9vCXqN7iHQTwMoDvAwjseZiU\npM0isgnV2/vFABYA+J6IPK6qP25fSxteP8n7bJ97HYBJVX22LY0Lfu0k7e5DdcjtAwAmAXxHRA6p\n6nf8XqcrgqmqfjTocRH5JICbAPyS1gY4asffDeBhAP9GVV9oayNdkrY5TQnbXEJjj+7dtWMdE9Zu\nn3OmANxjfy0i3wfwP022K+T1Y7cZwFoAD6tqBcAbIrIf1aGJjgTThG22fQIp9UoTtvtVAN9V1X8G\nABF5DMD7AfgG066/zReRXwHwaQDrVXXScbwIYA+ALaq6P632efFrc5YFtPkRAJ8Qkbm12+WfAfBU\nGm2MQ0T6RWR+7d+/DGBKVX+YcrPCPAfgBgCotX1N7VimiUgOwO3I2HhpiL0AVtb+n/QB+AUAwf8/\n0p5NMzAb9zyqY3bjtT//tXb8PgBnHcfHAbwz7fYGtbn22EsATgI4g+pvx/ek3d4Ibb4X1Vn8YwD+\nt7Tb6mr3x2vv43kA/wRgb+34slp7fwTgcQBXpN3WCG2eh+rwyrO1D/bmtNsa1ubaY78I4EDabUzQ\n7k2oTl4/C+A/hl2LK6CIiAzo+tt8IqIsYDAlIjKAwZSIyAAGUyIiAxhMiYgMYDAlIjKAwZQ6RkTO\ntHj+LhH56dq/XxKRd9T+/f0I594tIv0hz7lERP5CRF6olW38HyJynYgUReR3IrzGMhH5146vV4nI\nX4Z/Z9QLGEypK9SC6CWq+rz7MVX9cIRL3I3quvYg/w+qCyZ+RlWvBfAbqG5bXAQQGkxRXQhQD6Za\nrTh0tYj8VIRzqcsxmFLHSdWIiDxbK7q7sXY8JyJ/XitC/Xci8piI3Fo77RMAHvW53pna379Y601+\nrXaNnbXX+j1Ui4M8ISJP+FzjagDXAbhPVWcAQFVfVNU9AHagGhTHa+32bH/teT9fe5697v+bAG5r\n9T2j7OuKQifUc24BsBrAKlR7fj8Qke+iWshjGaq1Rt+J6lJP+zb5I6guEQ4zCGAFgNcA7AewVlX/\ns4j8PoDrtVa4wsMKAOPqXbNyC4D3qupqABCRDT7t3wLgD1T1Jse5T6FW7DtC26mLsWdKafgIgIdU\ndVpV/wnA36Na6uwjAL6qqjOq+hMAzl7kFajW8QzzlKq+WutdjqManE3za7+X19rUBsoYBlPqJhLh\nOc4C4NOIfvd1FMAqEcnHblUwAcACGLMAgyml4XsANopIXkQWAfhXqN4O7wewoTZ2+lOoVhuyvQzv\nAr9RvQXgUr8HtVrv9iCA7SIiQH12/mMe5/q13+s1Lq+1nXocgyml4WFUNys7jOo+TJ+u3dbvRrUc\n2g8BPAjgaVS36QCAf0DjHk19aOyFhnkAwLf8JqBq/i2AnwLwvIg8C+D/BfCGqr4JYH9twmkkoP3P\nAJiW6m6n9gTUB1ENvtTjWIKPMkVELlHVMyLydlR7e2tV9Se12fb/oqq/WusNjqtq5veaEpG/B7Cx\nFmyph3E2n7Lmv9d2SZgD4P+0g5CqviAib9XSnH4HQGb2jPcjIqtQ3dqYgXQWYM+UZh0ReRKAe9vv\nX1fVI2m0h3oDgykRkQGcgCIiMoDBlIjIAAZTIiIDGEyJiAxgMCUiMuD/B9CD75wGyc3eAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x15e13e790>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(np.log10(raw_data['Jint_Ctot'].values),np.log10(raw_data['Jint_Ptot'].values))\n",
    "plt.xlabel('log(Jint_Ctot)')\n",
    "plt.xlim([-22,-16])\n",
    "plt.ylim([-25,-13])\n",
    "plt.ylabel('log(Jint_Ptot)')\n",
    "plt.axes().set_aspect('equal')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}