my_error-space-diagram-single-parameterized.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Error-space diagram - example for 1 input pair, 1 derived var, 1 region"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load needed packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import rpy2\n",
    "\n",
    "# necessary to workaround the error: /u01/anaconda2/lib/libreadline.so.6: undefined symbol: PC\n",
    "# https://www.rittmanmead.com/blog/2016/07/using-r-jupyter-notebooks-big-data-discovery/\n",
    "import readline\n",
    "\n",
    "%load_ext rpy2.ipython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/homel/cnangini/anaconda/lib/python2.7/site-packages/rpy2/rinterface/__init__.py:186: RRuntimeWarning: Le chargement a nécessité le package : oce\n",
      "\n",
      "  warnings.warn(x, RRuntimeWarning)\n",
      "/homel/cnangini/anaconda/lib/python2.7/site-packages/rpy2/rinterface/__init__.py:186: RRuntimeWarning: Le chargement a nécessité le package : gsw\n",
      "\n",
      "  warnings.warn(x, RRuntimeWarning)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[1] ‘3.1.1’\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "#libraries previously installed using:\n",
    "library(seacarb)\n",
    "library(oce)\n",
    "library(gsw)\n",
    "\n",
    "#library(Hmisc)   #for minor ticks in plot\n",
    "\n",
    "packageVersion('seacarb')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Get routines and specify input vars"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 Read in new and modified functions "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The commands below read in routines that have been modified (derivnum and errors) or are completely new relative to the standard version of the seacarb distribution.  You can get these new and modified routines from github:\n",
    "\n",
    "git clone https://github.com/jamesorr/seacarb-git.git\n",
    "\n",
    "Then just change the directory in the 4 lines of the cell just below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "#branch buffesm1 of seacarb-git\n",
    "source(\"/homel/cnangini/PROJECTS/seacarb-git/R/errhalf.R\")\n",
    "source(\"/homel/cnangini/PROJECTS/seacarb-git/R/errmid.R\")\n",
    "source(\"/homel/cnangini/PROJECTS/seacarb-git/R/errors.R\")\n",
    "source(\"/homel/cnangini/PROJECTS/seacarb-git/R/derivnum.R\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Specify input variables and their corresponding uncertainties (errors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 892,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "#User inputs\n",
    "\n",
    "# Specify flag\n",
    "#   menu_flag = 15  #var1=ALK, var2=DIC\n",
    "   menu_flag = 8   #var1=pH, var2=ALK\n",
    "#  menu_flag = 25  #var1=pCO2, var2=DIC\n",
    "#   menu_flag = 21  #var1=pCO2, var2=pH\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 893,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "NULL\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "# max_error1  = 20\n",
    "# max_error2  = 10\n",
    "# redn = 1\n",
    "# numpts = 20\n",
    "\n",
    "# seq(0., max_error1, (max_error1/numpts)*(1.0/redn) )\n",
    "\n",
    "# seq(0., max_error2, (max_error2/numpts)*(1.0/redn) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 909,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "  temp = -0.49    #C\n",
    "  salt = 33.96    #psu\n",
    "  press = 0       #bar\n",
    "  Phos = 2.e-6    #(umol/kg)\n",
    "  Sil = 60.e-6    #(umol/kg)\n",
    " \n",
    "  numpts = 20\n",
    "  redn = 1  #0.5 #50% reduction\n",
    "  \n",
    "\n",
    "\n",
    "  if (menu_flag == 15){#var1=ALK, var2=DIC\n",
    "      \n",
    "      scalefactor1 = 1e+6 #ALK\n",
    "      scalefactor2 = 1e+6 #DIC\n",
    "      \n",
    "      max_error1  = 20\n",
    "      max_error2  = 10\n",
    "      \n",
    "      # Input variables:\n",
    "      var1 = 2295 * 1e-6  #(umol/kg)\n",
    "      var2 = 2155 * 1e-6  #(umol/kg)\n",
    "      \n",
    "      # Uncertainties in input variables      \n",
    "      var1_e <- seq(0., max_error1, (max_error1/numpts)*(1.0/redn) ) * 1e-6\n",
    "      var2_e <- seq(0., max_error2, (max_error2/numpts)*(1.0/redn) ) * 1e-6\n",
    "      \n",
    "      var1_e_soa   <- 2 #umol/kg\n",
    "      var2_e_soa   <- 2 #umol/kg\n",
    "\n",
    "      var1_e_soa2  <- c(var1_e_soa, var1_e_soa)\n",
    "      var2_e_soa2  <- c(var2_e_soa, var2_e_soa)\n",
    "\n",
    "      # for plot\n",
    "      x <- var2_e*1e+6  ;  y <- var1_e*1e+6      \n",
    "      xlim <- c(0,max_error2)  ; ylim <-c(0,max_error1)\n",
    "      levels1 <- c(1,seq(2,20,by=1))\n",
    "      xlabel <- expression(paste(sigma[italic(\"C\")[T]],\" (\",mu,\"mol kg\"^{-1},\")\",sep=\"\"))\n",
    "      ylabel <- expression(paste(sigma[italic(\"A\")[T]],\" (\",mu,\"mol kg\"^{-1},\")\",sep=\"\"))\n",
    "      \n",
    "  } else if (menu_flag == 8) {#var1=pH, var2=ALK\n",
    "      scalefactor1 = 1 #pH\n",
    "      scalefactor2 = 1e+6 #ALK\n",
    "      \n",
    "      max_error1  = 0.03\n",
    "      max_error2  = 20\n",
    "      \n",
    "      # Input variables:\n",
    "      var1 = 8.111823\n",
    "      var2 = 2295  * 1e-6  #(umol/kg)\n",
    "      \n",
    "      # Uncertainties in input variables      \n",
    "      var1_e <- seq(0, max_error1, (max_error1/numpts) * (1/redn) )\n",
    "      var2_e <- seq(0., max_error2, (max_error2/numpts)*(1.0/redn) ) * 1e-6\n",
    "      \n",
    "      var1_e_soa   <- c(0.003, 0.01)\n",
    "      var2_e_soa   <- 2 #umol/kg\n",
    "      \n",
    "      var1_e_soa2  <- var1_e_soa\n",
    "      var2_e_soa2  <- c(var2_e_soa, var2_e_soa)\n",
    "      \n",
    "      # for plot\n",
    "      x <- var1_e           ;  y <- var2_e * 1e+6\n",
    "      xlim <- c(0,max_error1) ; ylim <- c(0,max_error2) \n",
    "      levels1 <- c(4.2, seq(4,7,by=1))\n",
    "      xlabel <- expression(paste(sigma[pH],\" (total scale)\",sep=\"\"))\n",
    "      ylabel <- expression(paste(sigma[italic(\"A\")[T]],\" (\",mu,\"mol kg\"^{-1},\")\",sep=\"\"))\n",
    "  \n",
    "  } else if (menu_flag == 21) {#var1=pCO2, var2=pH\n",
    "      scalefactor1 = 1 #pCO2\n",
    "      scalefactor2 = 1 #pH\n",
    "      \n",
    "      # Input variables:\n",
    "      var1 = 330.5511 #uatm\n",
    "      var2 = 8.1\n",
    "      \n",
    "      # Uncertainties in input variables      \n",
    "      var1_e <- seq(0,maxlim,1/redn)\n",
    "      var2_e <- seq(0,0.03,0.0015/redn)\n",
    "      \n",
    "      var1_e_soa   <- 2\n",
    "      var2_e_soa   <- c(0.003, 0.01)\n",
    "      \n",
    "      var1_e_soa2  <- c(var1_e_soa, var1_e_soa)\n",
    "      var2_e_soa2  <- var2_e_soa\n",
    "      \n",
    "      # for plot\n",
    "      x <- var2_e           ;  y <- var1_e\n",
    "      xlim <- c(0,0.03)  ; ylim <- c(0,maxlim)\n",
    "      levels1 <- c(7,seq(0,maxlim,by=1))\n",
    "      xlabel <- expression(paste(sigma[pH],\" (total scale)\",sep=\"\"))\n",
    "      ylabel <- expression(paste(sigma[pCO[2]],\" (\",mu,\"atm\",\")\",sep=\"\"))    \n",
    "      \n",
    "\n",
    "   } else if (menu_flag == 25) {#var1=pCO2, var2=DIC\n",
    "      scalefactor1 = 1 #pCO2\n",
    "      scalefactor2 = 1e+6 #DIC\n",
    "      \n",
    "      # Input variables:\n",
    "      var1 = 330.5511 #uatm\n",
    "      var2 = 2155 * 1e-6  #(umol/kg)\n",
    "      \n",
    "      # Uncertainties in input variables      \n",
    "      var1_e <- seq(0,maxlim,1/redn)\n",
    "      var2_e <- seq(0., maxlim, 1/redn) * 1e-6\n",
    "      \n",
    "      var1_e_soa   <- 2\n",
    "      var2_e_soa   <- 2 #umol/kg\n",
    "      \n",
    "      var1_e_soa2  <- c(var1_e_soa, var1_e_soa)\n",
    "      var2_e_soa2  <- c(var2_e_soa, var2_e_soa)\n",
    "      \n",
    "      # for plot\n",
    "      x <- var1_e           ;  y <- var2_e * 1e+6\n",
    "      xlim <- c(0,maxlim)  ; ylim <- xlim\n",
    "      levels1 <- c(4.7,seq(1,maxlim,by=1))\n",
    "      xlabel <- expression(paste(sigma[pCO[2]],\" (\",mu,\"atm\",\")\",sep=\"\"))\n",
    "      ylabel <- expression(paste(sigma[italic(\"C\")[T]],\" (\",mu,\"mol kg\"^{-1},\")\",sep=\"\"))\n",
    "      \n",
    "  }\n",
    " \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 Default uncertainties in equilibrium constants (pK0, pK1, pK2, pKb, pKw, pKa, pKc, Bt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 910,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "# Estimates for errors in constants in errors.R routine (from Orr et al., 2017)\n",
    "#   epKstd  <- c(0.004, 0.015,  0.03, 0.01,  0.01, 0.02, 0.02, 0.01)\n",
    "epKstd  <- c(0.002, 0.01,  0.02, 0.01,  0.01, 0.01, 0.01, 0.00)\n",
    "# epKstd  <- c(0.004, 0.015,  0.03, 0.01,  0.01, 0.02, 0.02)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 Compute derived carbonate system vars with std seacarb routine `carb`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 911,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "# Compute derived carbonate system variables with standard seacarb routine 'carb':\n",
    "\n",
    "# Southern Ocean\n",
    "vars <- carb  (flag=menu_flag, var1=var1, var2=var2, S=salt, T=temp, Patm=1, P=press, Pt=Phos, Sit=Sil, \n",
    "                    k1k2='w14', kf='dg', ks=\"d\", pHscale=\"T\", b=\"u74\", gas=\"potential\", warn='n')\n",
    "\n",
    "pH <- vars$pH\n",
    "pCO2 <- vars$pCO2\n",
    "\n",
    "# Compute H+ from pH, i.e., pH = -log10[H+], where [H+] is the hydrogen ion concentration in mol/kg\n",
    "H = 10^(-1*vars$pH)\n",
    "\n",
    "# Keep only key columns to be consistent with output from 'errors.R' routine (called below)\n",
    "vars <- data.frame(H, vars[,c('pH','CO2','fCO2','pCO2','HCO3','CO3','OmegaAragonite','OmegaCalcite')] )\n",
    "vars\n",
    "\n",
    "# Duplicate rows in *vars* until same as number of members of error vector ALK_e\n",
    "numerrs <- length(var1_e)\n",
    "vars <- vars[rep(row.names(vars), numerrs), ]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4 Use 1-D error vectors to build 2-D error array (to plot contours in DIC-ALK space)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 912,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "if (menu_flag == 15) dat <- expand.grid(var2_e, var1_e)\n",
    "    \n",
    "else if (menu_flag == 8 || menu_flag == 9 || menu_flag == 21 || menu_flag == 24 || menu_flag == 25) {\n",
    "    dat <- expand.grid(var1_e,  var2_e)\n",
    "}\n",
    "\n",
    "#   dat_pCO2Ct <- expand.grid(pCO2_e,  DIC_e) #flag 25"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 913,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "      Var1    Var2\n",
       "1   0.0000 0.0e+00\n",
       "2   0.0015 0.0e+00\n",
       "3   0.0030 0.0e+00\n",
       "4   0.0045 0.0e+00\n",
       "5   0.0060 0.0e+00\n",
       "6   0.0075 0.0e+00\n",
       "7   0.0090 0.0e+00\n",
       "8   0.0105 0.0e+00\n",
       "9   0.0120 0.0e+00\n",
       "10  0.0135 0.0e+00\n",
       "11  0.0150 0.0e+00\n",
       "12  0.0165 0.0e+00\n",
       "13  0.0180 0.0e+00\n",
       "14  0.0195 0.0e+00\n",
       "15  0.0210 0.0e+00\n",
       "16  0.0225 0.0e+00\n",
       "17  0.0240 0.0e+00\n",
       "18  0.0255 0.0e+00\n",
       "19  0.0270 0.0e+00\n",
       "20  0.0285 0.0e+00\n",
       "21  0.0300 0.0e+00\n",
       "22  0.0000 1.0e-06\n",
       "23  0.0015 1.0e-06\n",
       "24  0.0030 1.0e-06\n",
       "25  0.0045 1.0e-06\n",
       "26  0.0060 1.0e-06\n",
       "27  0.0075 1.0e-06\n",
       "28  0.0090 1.0e-06\n",
       "29  0.0105 1.0e-06\n",
       "30  0.0120 1.0e-06\n",
       "31  0.0135 1.0e-06\n",
       "32  0.0150 1.0e-06\n",
       "33  0.0165 1.0e-06\n",
       "34  0.0180 1.0e-06\n",
       "35  0.0195 1.0e-06\n",
       "36  0.0210 1.0e-06\n",
       "37  0.0225 1.0e-06\n",
       "38  0.0240 1.0e-06\n",
       "39  0.0255 1.0e-06\n",
       "40  0.0270 1.0e-06\n",
       "41  0.0285 1.0e-06\n",
       "42  0.0300 1.0e-06\n",
       "43  0.0000 2.0e-06\n",
       "44  0.0015 2.0e-06\n",
       "45  0.0030 2.0e-06\n",
       "46  0.0045 2.0e-06\n",
       "47  0.0060 2.0e-06\n",
       "48  0.0075 2.0e-06\n",
       "49  0.0090 2.0e-06\n",
       "50  0.0105 2.0e-06\n",
       "51  0.0120 2.0e-06\n",
       "52  0.0135 2.0e-06\n",
       "53  0.0150 2.0e-06\n",
       "54  0.0165 2.0e-06\n",
       "55  0.0180 2.0e-06\n",
       "56  0.0195 2.0e-06\n",
       "57  0.0210 2.0e-06\n",
       "58  0.0225 2.0e-06\n",
       "59  0.0240 2.0e-06\n",
       "60  0.0255 2.0e-06\n",
       "61  0.0270 2.0e-06\n",
       "62  0.0285 2.0e-06\n",
       "63  0.0300 2.0e-06\n",
       "64  0.0000 3.0e-06\n",
       "65  0.0015 3.0e-06\n",
       "66  0.0030 3.0e-06\n",
       "67  0.0045 3.0e-06\n",
       "68  0.0060 3.0e-06\n",
       "69  0.0075 3.0e-06\n",
       "70  0.0090 3.0e-06\n",
       "71  0.0105 3.0e-06\n",
       "72  0.0120 3.0e-06\n",
       "73  0.0135 3.0e-06\n",
       "74  0.0150 3.0e-06\n",
       "75  0.0165 3.0e-06\n",
       "76  0.0180 3.0e-06\n",
       "77  0.0195 3.0e-06\n",
       "78  0.0210 3.0e-06\n",
       "79  0.0225 3.0e-06\n",
       "80  0.0240 3.0e-06\n",
       "81  0.0255 3.0e-06\n",
       "82  0.0270 3.0e-06\n",
       "83  0.0285 3.0e-06\n",
       "84  0.0300 3.0e-06\n",
       "85  0.0000 4.0e-06\n",
       "86  0.0015 4.0e-06\n",
       "87  0.0030 4.0e-06\n",
       "88  0.0045 4.0e-06\n",
       "89  0.0060 4.0e-06\n",
       "90  0.0075 4.0e-06\n",
       "91  0.0090 4.0e-06\n",
       "92  0.0105 4.0e-06\n",
       "93  0.0120 4.0e-06\n",
       "94  0.0135 4.0e-06\n",
       "95  0.0150 4.0e-06\n",
       "96  0.0165 4.0e-06\n",
       "97  0.0180 4.0e-06\n",
       "98  0.0195 4.0e-06\n",
       "99  0.0210 4.0e-06\n",
       "100 0.0225 4.0e-06\n",
       "101 0.0240 4.0e-06\n",
       "102 0.0255 4.0e-06\n",
       "103 0.0270 4.0e-06\n",
       "104 0.0285 4.0e-06\n",
       "105 0.0300 4.0e-06\n",
       "106 0.0000 5.0e-06\n",
       "107 0.0015 5.0e-06\n",
       "108 0.0030 5.0e-06\n",
       "109 0.0045 5.0e-06\n",
       "110 0.0060 5.0e-06\n",
       "111 0.0075 5.0e-06\n",
       "112 0.0090 5.0e-06\n",
       "113 0.0105 5.0e-06\n",
       "114 0.0120 5.0e-06\n",
       "115 0.0135 5.0e-06\n",
       "116 0.0150 5.0e-06\n",
       "117 0.0165 5.0e-06\n",
       "118 0.0180 5.0e-06\n",
       "119 0.0195 5.0e-06\n",
       "120 0.0210 5.0e-06\n",
       "121 0.0225 5.0e-06\n",
       "122 0.0240 5.0e-06\n",
       "123 0.0255 5.0e-06\n",
       "124 0.0270 5.0e-06\n",
       "125 0.0285 5.0e-06\n",
       "126 0.0300 5.0e-06\n",
       "127 0.0000 6.0e-06\n",
       "128 0.0015 6.0e-06\n",
       "129 0.0030 6.0e-06\n",
       "130 0.0045 6.0e-06\n",
       "131 0.0060 6.0e-06\n",
       "132 0.0075 6.0e-06\n",
       "133 0.0090 6.0e-06\n",
       "134 0.0105 6.0e-06\n",
       "135 0.0120 6.0e-06\n",
       "136 0.0135 6.0e-06\n",
       "137 0.0150 6.0e-06\n",
       "138 0.0165 6.0e-06\n",
       "139 0.0180 6.0e-06\n",
       "140 0.0195 6.0e-06\n",
       "141 0.0210 6.0e-06\n",
       "142 0.0225 6.0e-06\n",
       "143 0.0240 6.0e-06\n",
       "144 0.0255 6.0e-06\n",
       "145 0.0270 6.0e-06\n",
       "146 0.0285 6.0e-06\n",
       "147 0.0300 6.0e-06\n",
       "148 0.0000 7.0e-06\n",
       "149 0.0015 7.0e-06\n",
       "150 0.0030 7.0e-06\n",
       "151 0.0045 7.0e-06\n",
       "152 0.0060 7.0e-06\n",
       "153 0.0075 7.0e-06\n",
       "154 0.0090 7.0e-06\n",
       "155 0.0105 7.0e-06\n",
       "156 0.0120 7.0e-06\n",
       "157 0.0135 7.0e-06\n",
       "158 0.0150 7.0e-06\n",
       "159 0.0165 7.0e-06\n",
       "160 0.0180 7.0e-06\n",
       "161 0.0195 7.0e-06\n",
       "162 0.0210 7.0e-06\n",
       "163 0.0225 7.0e-06\n",
       "164 0.0240 7.0e-06\n",
       "165 0.0255 7.0e-06\n",
       "166 0.0270 7.0e-06\n",
       "167 0.0285 7.0e-06\n",
       "168 0.0300 7.0e-06\n",
       "169 0.0000 8.0e-06\n",
       "170 0.0015 8.0e-06\n",
       "171 0.0030 8.0e-06\n",
       "172 0.0045 8.0e-06\n",
       "173 0.0060 8.0e-06\n",
       "174 0.0075 8.0e-06\n",
       "175 0.0090 8.0e-06\n",
       "176 0.0105 8.0e-06\n",
       "177 0.0120 8.0e-06\n",
       "178 0.0135 8.0e-06\n",
       "179 0.0150 8.0e-06\n",
       "180 0.0165 8.0e-06\n",
       "181 0.0180 8.0e-06\n",
       "182 0.0195 8.0e-06\n",
       "183 0.0210 8.0e-06\n",
       "184 0.0225 8.0e-06\n",
       "185 0.0240 8.0e-06\n",
       "186 0.0255 8.0e-06\n",
       "187 0.0270 8.0e-06\n",
       "188 0.0285 8.0e-06\n",
       "189 0.0300 8.0e-06\n",
       "190 0.0000 9.0e-06\n",
       "191 0.0015 9.0e-06\n",
       "192 0.0030 9.0e-06\n",
       "193 0.0045 9.0e-06\n",
       "194 0.0060 9.0e-06\n",
       "195 0.0075 9.0e-06\n",
       "196 0.0090 9.0e-06\n",
       "197 0.0105 9.0e-06\n",
       "198 0.0120 9.0e-06\n",
       "199 0.0135 9.0e-06\n",
       "200 0.0150 9.0e-06\n",
       "201 0.0165 9.0e-06\n",
       "202 0.0180 9.0e-06\n",
       "203 0.0195 9.0e-06\n",
       "204 0.0210 9.0e-06\n",
       "205 0.0225 9.0e-06\n",
       "206 0.0240 9.0e-06\n",
       "207 0.0255 9.0e-06\n",
       "208 0.0270 9.0e-06\n",
       "209 0.0285 9.0e-06\n",
       "210 0.0300 9.0e-06\n",
       "211 0.0000 1.0e-05\n",
       "212 0.0015 1.0e-05\n",
       "213 0.0030 1.0e-05\n",
       "214 0.0045 1.0e-05\n",
       "215 0.0060 1.0e-05\n",
       "216 0.0075 1.0e-05\n",
       "217 0.0090 1.0e-05\n",
       "218 0.0105 1.0e-05\n",
       "219 0.0120 1.0e-05\n",
       "220 0.0135 1.0e-05\n",
       "221 0.0150 1.0e-05\n",
       "222 0.0165 1.0e-05\n",
       "223 0.0180 1.0e-05\n",
       "224 0.0195 1.0e-05\n",
       "225 0.0210 1.0e-05\n",
       "226 0.0225 1.0e-05\n",
       "227 0.0240 1.0e-05\n",
       "228 0.0255 1.0e-05\n",
       "229 0.0270 1.0e-05\n",
       "230 0.0285 1.0e-05\n",
       "231 0.0300 1.0e-05\n",
       "232 0.0000 1.1e-05\n",
       "233 0.0015 1.1e-05\n",
       "234 0.0030 1.1e-05\n",
       "235 0.0045 1.1e-05\n",
       "236 0.0060 1.1e-05\n",
       "237 0.0075 1.1e-05\n",
       "238 0.0090 1.1e-05\n",
       "239 0.0105 1.1e-05\n",
       "240 0.0120 1.1e-05\n",
       "241 0.0135 1.1e-05\n",
       "242 0.0150 1.1e-05\n",
       "243 0.0165 1.1e-05\n",
       "244 0.0180 1.1e-05\n",
       "245 0.0195 1.1e-05\n",
       "246 0.0210 1.1e-05\n",
       "247 0.0225 1.1e-05\n",
       "248 0.0240 1.1e-05\n",
       "249 0.0255 1.1e-05\n",
       "250 0.0270 1.1e-05\n",
       "251 0.0285 1.1e-05\n",
       "252 0.0300 1.1e-05\n",
       "253 0.0000 1.2e-05\n",
       "254 0.0015 1.2e-05\n",
       "255 0.0030 1.2e-05\n",
       "256 0.0045 1.2e-05\n",
       "257 0.0060 1.2e-05\n",
       "258 0.0075 1.2e-05\n",
       "259 0.0090 1.2e-05\n",
       "260 0.0105 1.2e-05\n",
       "261 0.0120 1.2e-05\n",
       "262 0.0135 1.2e-05\n",
       "263 0.0150 1.2e-05\n",
       "264 0.0165 1.2e-05\n",
       "265 0.0180 1.2e-05\n",
       "266 0.0195 1.2e-05\n",
       "267 0.0210 1.2e-05\n",
       "268 0.0225 1.2e-05\n",
       "269 0.0240 1.2e-05\n",
       "270 0.0255 1.2e-05\n",
       "271 0.0270 1.2e-05\n",
       "272 0.0285 1.2e-05\n",
       "273 0.0300 1.2e-05\n",
       "274 0.0000 1.3e-05\n",
       "275 0.0015 1.3e-05\n",
       "276 0.0030 1.3e-05\n",
       "277 0.0045 1.3e-05\n",
       "278 0.0060 1.3e-05\n",
       "279 0.0075 1.3e-05\n",
       "280 0.0090 1.3e-05\n",
       "281 0.0105 1.3e-05\n",
       "282 0.0120 1.3e-05\n",
       "283 0.0135 1.3e-05\n",
       "284 0.0150 1.3e-05\n",
       "285 0.0165 1.3e-05\n",
       "286 0.0180 1.3e-05\n",
       "287 0.0195 1.3e-05\n",
       "288 0.0210 1.3e-05\n",
       "289 0.0225 1.3e-05\n",
       "290 0.0240 1.3e-05\n",
       "291 0.0255 1.3e-05\n",
       "292 0.0270 1.3e-05\n",
       "293 0.0285 1.3e-05\n",
       "294 0.0300 1.3e-05\n",
       "295 0.0000 1.4e-05\n",
       "296 0.0015 1.4e-05\n",
       "297 0.0030 1.4e-05\n",
       "298 0.0045 1.4e-05\n",
       "299 0.0060 1.4e-05\n",
       "300 0.0075 1.4e-05\n",
       "301 0.0090 1.4e-05\n",
       "302 0.0105 1.4e-05\n",
       "303 0.0120 1.4e-05\n",
       "304 0.0135 1.4e-05\n",
       "305 0.0150 1.4e-05\n",
       "306 0.0165 1.4e-05\n",
       "307 0.0180 1.4e-05\n",
       "308 0.0195 1.4e-05\n",
       "309 0.0210 1.4e-05\n",
       "310 0.0225 1.4e-05\n",
       "311 0.0240 1.4e-05\n",
       "312 0.0255 1.4e-05\n",
       "313 0.0270 1.4e-05\n",
       "314 0.0285 1.4e-05\n",
       "315 0.0300 1.4e-05\n",
       "316 0.0000 1.5e-05\n",
       "317 0.0015 1.5e-05\n",
       "318 0.0030 1.5e-05\n",
       "319 0.0045 1.5e-05\n",
       "320 0.0060 1.5e-05\n",
       "321 0.0075 1.5e-05\n",
       "322 0.0090 1.5e-05\n",
       "323 0.0105 1.5e-05\n",
       "324 0.0120 1.5e-05\n",
       "325 0.0135 1.5e-05\n",
       "326 0.0150 1.5e-05\n",
       "327 0.0165 1.5e-05\n",
       "328 0.0180 1.5e-05\n",
       "329 0.0195 1.5e-05\n",
       "330 0.0210 1.5e-05\n",
       "331 0.0225 1.5e-05\n",
       "332 0.0240 1.5e-05\n",
       "333 0.0255 1.5e-05\n",
       "334 0.0270 1.5e-05\n",
       "335 0.0285 1.5e-05\n",
       "336 0.0300 1.5e-05\n",
       "337 0.0000 1.6e-05\n",
       "338 0.0015 1.6e-05\n",
       "339 0.0030 1.6e-05\n",
       "340 0.0045 1.6e-05\n",
       "341 0.0060 1.6e-05\n",
       "342 0.0075 1.6e-05\n",
       "343 0.0090 1.6e-05\n",
       "344 0.0105 1.6e-05\n",
       "345 0.0120 1.6e-05\n",
       "346 0.0135 1.6e-05\n",
       "347 0.0150 1.6e-05\n",
       "348 0.0165 1.6e-05\n",
       "349 0.0180 1.6e-05\n",
       "350 0.0195 1.6e-05\n",
       "351 0.0210 1.6e-05\n",
       "352 0.0225 1.6e-05\n",
       "353 0.0240 1.6e-05\n",
       "354 0.0255 1.6e-05\n",
       "355 0.0270 1.6e-05\n",
       "356 0.0285 1.6e-05\n",
       "357 0.0300 1.6e-05\n",
       "358 0.0000 1.7e-05\n",
       "359 0.0015 1.7e-05\n",
       "360 0.0030 1.7e-05\n",
       "361 0.0045 1.7e-05\n",
       "362 0.0060 1.7e-05\n",
       "363 0.0075 1.7e-05\n",
       "364 0.0090 1.7e-05\n",
       "365 0.0105 1.7e-05\n",
       "366 0.0120 1.7e-05\n",
       "367 0.0135 1.7e-05\n",
       "368 0.0150 1.7e-05\n",
       "369 0.0165 1.7e-05\n",
       "370 0.0180 1.7e-05\n",
       "371 0.0195 1.7e-05\n",
       "372 0.0210 1.7e-05\n",
       "373 0.0225 1.7e-05\n",
       "374 0.0240 1.7e-05\n",
       "375 0.0255 1.7e-05\n",
       "376 0.0270 1.7e-05\n",
       "377 0.0285 1.7e-05\n",
       "378 0.0300 1.7e-05\n",
       "379 0.0000 1.8e-05\n",
       "380 0.0015 1.8e-05\n",
       "381 0.0030 1.8e-05\n",
       "382 0.0045 1.8e-05\n",
       "383 0.0060 1.8e-05\n",
       "384 0.0075 1.8e-05\n",
       "385 0.0090 1.8e-05\n",
       "386 0.0105 1.8e-05\n",
       "387 0.0120 1.8e-05\n",
       "388 0.0135 1.8e-05\n",
       "389 0.0150 1.8e-05\n",
       "390 0.0165 1.8e-05\n",
       "391 0.0180 1.8e-05\n",
       "392 0.0195 1.8e-05\n",
       "393 0.0210 1.8e-05\n",
       "394 0.0225 1.8e-05\n",
       "395 0.0240 1.8e-05\n",
       "396 0.0255 1.8e-05\n",
       "397 0.0270 1.8e-05\n",
       "398 0.0285 1.8e-05\n",
       "399 0.0300 1.8e-05\n",
       "400 0.0000 1.9e-05\n",
       "401 0.0015 1.9e-05\n",
       "402 0.0030 1.9e-05\n",
       "403 0.0045 1.9e-05\n",
       "404 0.0060 1.9e-05\n",
       "405 0.0075 1.9e-05\n",
       "406 0.0090 1.9e-05\n",
       "407 0.0105 1.9e-05\n",
       "408 0.0120 1.9e-05\n",
       "409 0.0135 1.9e-05\n",
       "410 0.0150 1.9e-05\n",
       "411 0.0165 1.9e-05\n",
       "412 0.0180 1.9e-05\n",
       "413 0.0195 1.9e-05\n",
       "414 0.0210 1.9e-05\n",
       "415 0.0225 1.9e-05\n",
       "416 0.0240 1.9e-05\n",
       "417 0.0255 1.9e-05\n",
       "418 0.0270 1.9e-05\n",
       "419 0.0285 1.9e-05\n",
       "420 0.0300 1.9e-05\n",
       "421 0.0000 2.0e-05\n",
       "422 0.0015 2.0e-05\n",
       "423 0.0030 2.0e-05\n",
       "424 0.0045 2.0e-05\n",
       "425 0.0060 2.0e-05\n",
       "426 0.0075 2.0e-05\n",
       "427 0.0090 2.0e-05\n",
       "428 0.0105 2.0e-05\n",
       "429 0.0120 2.0e-05\n",
       "430 0.0135 2.0e-05\n",
       "431 0.0150 2.0e-05\n",
       "432 0.0165 2.0e-05\n",
       "433 0.0180 2.0e-05\n",
       "434 0.0195 2.0e-05\n",
       "435 0.0210 2.0e-05\n",
       "436 0.0225 2.0e-05\n",
       "437 0.0240 2.0e-05\n",
       "438 0.0255 2.0e-05\n",
       "439 0.0270 2.0e-05\n",
       "440 0.0285 2.0e-05\n",
       "441 0.0300 2.0e-05\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "dat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.5 Specify state-of-art (soa) errors for vars: from error-propagation manuscript (Orr et al. 2017, Table 1) \n",
    "\n",
    "These will be plotted as crosses on the error-space diagrams"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Compute derived variables (with `carb`) and propagated errors (with `errors`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Computations of derived vars and their errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 914,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "# Compute variables\n",
    "vars <- carb  (flag=menu_flag, var1=var1, var2=var2, S=salt, T=temp, Patm=1, P=press, Pt=Phos, Sit=Sil, \n",
    "                    k1k2='w14', kf='dg', ks=\"d\", pHscale=\"T\", b=\"u74\", gas=\"potential\", warn='n')\n",
    "\n",
    "H = 10^(-1*vars$pH)         # H+ concentration (mol/kg)\n",
    "vars <- data.frame(H, vars) # Add H+ as new column to vars data frame\n",
    "\n",
    "# Absolute errors: propagated uncertainties\n",
    "if (menu_flag == 15) {\n",
    "    absEt <- errors (flag=menu_flag, var1=var1, var2=var2, S=salt, T=temp, \n",
    "                      Patm=1, P=press, Pt=Phos, Sit=Sil, \n",
    "                      evar1=dat$Var2, evar2=dat$Var1, eS=0, eT=0, ePt=0, eSit=0, epK=epKstd, #HUOM!!!                \n",
    "                      k1k2='w14', kf='dg', ks=\"d\", pHscale=\"T\", b=\"u74\", gas=\"potential\", warn='no')\n",
    "    \n",
    "} else if (menu_flag == 8 || menu_flag == 21 || menu_flag == 25) { ##NB: SWAP dat$Var2 and dat$Var1 !!!!!!\n",
    "    absEt <- errors (flag=menu_flag, var1=var1, var2=var2, S=salt, T=temp, \n",
    "                      Patm=1, P=press, Pt=Phos, Sit=Sil, \n",
    "                      evar1=dat$Var1, evar2=dat$Var2, eS=0, eT=0, ePt=0, eSit=0, epK=epKstd, #HUOM!!!\n",
    "                      k1k2='w14', kf='dg', ks=\"d\", pHscale=\"T\", b=\"u74\", gas=\"potential\", warn='no')\n",
    "}\n",
    "\n",
    "\n",
    "# Keep only key columns in vars for consistency with columns in absEt\n",
    "vars <- vars[,colnames(absEt)] \n",
    "vars\n",
    "\n",
    "# Duplicate rows in *vars* until same as number of members of error vector ALK_e\n",
    "numerrs <- length(dat$Var1)\n",
    "vars <- vars[rep(row.names(vars), numerrs), ]\n",
    "\n",
    "#Relative errors (in percent)\n",
    "relEt <- 100* absEt / vars            #Total relative error (from constants and other input variables)\n",
    "# relEt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Define simpler names of variables to be plotted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 915,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "#   Relative changes:\n",
    "    rCO3  <- relEt$CO3\n",
    "    rH    <- relEt$H\n",
    "    rpCO2 <- relEt$pCO2\n",
    "    rCO2  <- relEt$CO2\n",
    "    #rpCO2   <- relEt$pCO2 ###DUPLICATE\n",
    "    rOmegaA <- relEt$OmegaAragonite\n",
    "    rHCO3   <- relEt$HCO3\n",
    "    H       <- absEt$H\n",
    "\n",
    "    rHCO3H <- sqrt(rHCO3^2 + rH^2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Compute other parts of error-space diagrams"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Constants-pair CURVE (along which propagated error from constants equals that from the input pair)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 916,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "# Critical curves (Constants-Pairs curve)\n",
    "# ===============\n",
    "# Southern Ocean\n",
    "errcirc <- errhalf(flag=menu_flag, var1=var1, var2=var2, S=salt, T=temp, \n",
    "                   Patm=1, P=press, Pt=Phos, Sit=Sil,\n",
    "                   epK=epKstd,\n",
    "                   k1k2='l', kf='dg', ks=\"d\", pHscale=\"T\", b=\"u74\", gas=\"potential\", warn=\"n\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Balanced-pair LINE (along which both members of input pair contrbute equally to propagated error)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 917,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "sigyspct <- seq(0,100,by=1) # in percent  ##ALWAYS STAYS THE SAME, error will never go > 100%\n",
    "\n",
    "\n",
    "# At-Ct pair\n",
    "# ----------\n",
    "# Southern Ocean\n",
    "errm <- errmid(flag=menu_flag, var1=var1, var2=var2, S=salt, T=temp, \n",
    "                   Patm=1, P=press, Pt=Phos, Sit=Sil,\n",
    "                   sigyspct, epK=epKstd,\n",
    "                   k1k2='l', kf='dg', ks=\"d\", pHscale=\"T\", b=\"u74\", gas=\"potential\", warn=\"n\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 918,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " [1]  0.000000  1.633234  8.156225 16.250375 24.220851 32.006991 39.549539\n",
       " [8] 46.791091 53.676535 60.153467 66.172596 71.688111 76.658036 81.044547\n",
       "[15] 84.814260 87.938486 90.393447 92.160458 93.226074 93.567929 93.582182\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "# Apply scale factors when necessary\n",
    "sig1   <- data.frame(errcirc[1]) * scalefactor1\n",
    "sig2   <- data.frame(errcirc[2]) * scalefactor2\n",
    "    \n",
    "sigm1   <- data.frame(errm[1]) * scalefactor1\n",
    "sigm2   <- data.frame(errm[2]) * scalefactor2\n",
    "\n",
    "sig2$CO3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Plot results: contours of calculated errors from input errors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1 Create file 'prettyaxes.R' to beautify subsequent plots (call it in cell below)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 919,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "cat(' box()\\n',\n",
    "    'axis(1,labels=FALSE, tcl=-0.35) ; axis(1, lwd=0, line=-0.4, labels=TRUE)\\n',\n",
    "    'axis(2,labels=FALSE,tcl=-0.35)  ; axis(2, lwd=0, line=-0.4, las=2, labels=TRUE)\\n',\n",
    "    'axis(3,labels=FALSE,tcl=-0.35) ; axis(4,labels=FALSE,tcl=-0.25)\\n',\n",
    "   file='prettyaxes.R')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2 Function to strip NaNs and add a 0 at beginning of a vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 920,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "zenon <- function(x){\n",
    "       # Get rid of NaN's\n",
    "       x <- c(0.0, x[!is.na(x)]) \n",
    "       return(x)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.3 Function to make plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 921,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "plterrcontour <- function(sigcritXa, sigcritYa, xlab, ylab, subtitle, xlim, ylim,\n",
    "                           sighpXa, sighpYa,\n",
    "                           sigmidXa, sigmidYa,\n",
    "                           sigsoaX, sigsoaY,\n",
    "                           x, y, za, levels1,\n",
    "                           methoda='flattest'){\n",
    "    \n",
    "#   Region 'a' is Southern Ocean  ;  Region 'b' is Tropics:\n",
    "    # 'contour' command must come after any 'plot', 'lines', or 'points' commands \n",
    "    plot (sigcritXa, sigcritYa, xlab=xlab, ylab=ylab, plot.title=title(main=subtitle),\n",
    "          xlim=xlim, ylim=ylim, col=\"black\", type='l', lwd=4, lty='solid',xaxt='n',yaxt='n',\n",
    "          cex.lab=1.4, ann=F, xaxs='r', yaxs='r')\n",
    "#   Plot critical point\n",
    "    if ( is.null(sighpXa) ) {\n",
    "        points(sigcritXa[11], sigcritYa[11], col=\"black\", cex=1.7, pch=19)\n",
    "    } else {\n",
    "        points(sighpXa, sighpYa, col=\"black\", cex=1.7, pch=19)\n",
    "    }    \n",
    "#   Plot mid-line (where e1 = e2 on each isoline [and between isolines])\n",
    "    lines(sigmidXa, sigmidYa, col=\"black\", lwd=1, lty='solid')\n",
    "#   symbols(x=sigcritXb[11], y=sigcritYb[11], circles=1, inches=1/3, ann=F, bg=\"steelblue2\", fg=NULL, add=TRUE)\n",
    "#   Plot point(s) at state-of-the-art (soa) for absolute error in each member of the input pair\n",
    "#   - pH has 2 points (random & overall error); other 3 input variables only have 1 point\n",
    "    points(sigsoaX, sigsoaY, col=\"black\", cex=1.7, pch=4)\n",
    "    axis(1,labels=TRUE, tcl=-0.35, cex.axis=1.3) \n",
    "    axis(2, tcl=-0.35, labels=FALSE, las=0) \n",
    "    title(xlab = xlab, cex.lab = 1.4, line = 2.9)\n",
    "    title(ylab = ylab, cex.lab = 1.4, line = 1.9)\n",
    "\n",
    "#   Add 2 colored text strings near origin to indicate the propagated error at the origin (from constants only)\n",
    "    a0 <- za[1,1]\n",
    "    if (a0 < 1) {ndiga <- 2} else {ndiga <- 3}\n",
    "    mtext( format(a0, digits=ndiga), side=1, line=2.5, at=0, srt=0,  cex=1.05, col='blue')\n",
    "\n",
    "#  'contour' command MUST come after any 'plot', 'lines', or 'points' commands \n",
    "#   For 'method', use 'simple, edge or flattest (flattest looks best but seems to fail often when lines are vertical)\n",
    "    contour(x, y , za, levels=levels1, method=methoda, col=\"blue\", lty=1, lwd=2, \n",
    "            yaxt='n', , labcex=1.0, add=TRUE)\n",
    "    axis(2, tcl=-0.35, labels=FALSE) \n",
    "\n",
    "    axis(2, tcl=-0.35, labels=FALSE, las=0) ; \n",
    "    axis(2, lwd=0, line=-0.4, las=2, labels=TRUE, cex.axis=1.4)\n",
    "    axis(3,labels=FALSE,tcl=-0.35) ; axis(4,labels=FALSE,tcl=-0.25)\n",
    "    done <- c('error contour plot')\n",
    "#return(NULL)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.4 Error-space diagram of relative error in CO3  for At-Ct input pair\n",
    "\n",
    "Combine (1) contour plot of CO3 with (2) constants-pair curve, (3) balanced-pair line, ane (4) symbol indicating current state-of-the-art uncertainty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 922,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAHgCAIAAADytinCAAAgAElEQVR4nO3deXwTdf4/8Pc0Tdqm\nd2tPjpazQLkPqYLcqOU+lIUFFikiiqDgF1DcrrCCcskqAspvUVYUFESoIiuXIIcIKCKKJC09KbSl\noS2F3lfm98dotrRNyD2T+byej33sI5lM0ndjefXd93xmwvE8TwAAID1uYhcAAABNQ0ADAEgUAhoA\nQKIQ0AAAEoWABgCQKAQ0AIBEIaABACQKAQ0AIFEIaAAAiUJAAwBIFAIaAECiENAAABKFgAYAkCgE\nNACARCGgAQAkCgENACBRCGgAAIlCQAMASBQCGgBAohDQAAAShYAGAJAoBDQAgEQhoAEAJMpd7AIA\n7O/mzZtEFBYWxnGc858OYC8cz/Ni1wBgZ0KwVlRUeHp6Ov/pAPaCEQcAgEQhoEFuDHMJLy+v5ORk\nIrpx48bEiRNDQkJCQ0Nnz55dXFws7JCRkTFq1Ch/f/+QkJBJkybduHGjyacDiAUzaJCbmJiYlJQU\nImrXrp1KpaqoqOjXr19ubu7IkSOrqqo++OADrVZ76tQpNze36dOn//DDD0OHDr1z586ePXsKCgqO\nHz/e4OlifzfANMygQYbqD5H/3//7f88+++zcuXM3b95MRI899tiRI0e+/fbbwYMHe3l5VVdXHz9+\nvEePHomJie7u7m+//TbHcZhBg0RgxAEy9+uvvxLRe++9JyTvkSNHiEir1bq5uQ0fPpyIhgwZ0rt3\nbz8/v9dffx3LNkBSMOIAmbtz5w4RRURE+Pn5GTZWVFQQUVJS0vbt27dv3/7999+vWrXqzJkzJ0+e\nFK1QgEYQ0CBbdXV1RBQTE0NECQkJK1euJKLz58/n5eV169aN5/m///3v4eHhp0+f/uWXX3r27Hn+\n/Hm9Xu/m5lb/6QAiwgwaZMjHx6esrGzq1Knr16/X6/UxMTGVlZV/+ctfqqqq9u7dGxgYePny5YiI\niD59+ly4cGHEiBElJSWnT58eOHDgiRMnGjw9LCxM7O8G2KVYvny52DUA2JlKpfrxxx9TU1Ofeuqp\nli1bjhgx4urVq9999921a9dGjBjx2WefRUVFEdGgQYPS0tJ++OGHvLy8oUOHbt26NSAgoMHTg4KC\nxP5ugF3ooAEAJAqrOAAAJAoBDQAgUQhoAACJQkADAEgUAhoAQKIQ0AAAEoWABgCQKAQ0AIBEIaAB\nACQKAQ0AIFEIaAAAiUJAAwBIFAIaAECiENAAABKFgAYAkCgENACARCGgAQAkSs4Bfe7cObFLsNhv\nv/1WXl4udhWWuX37dkpKithVWMwVfzx+/PFHvV4vdhWWuX79ek5OjthVWKauru6nn34SuwoieQf0\n0qVLxS7BYtu2bUtNTRW7CstcuHDhiy++ELsKi7nij8eaNWuKi4vFrsIyR48ePXbsmNhVWKawsPCt\nt94SuwoieQc0AIBLQ0ADAEgUAhoAQKLcxS7AMoWFhfv27eN53pyd09LSpk+f7uiS7OvChQtpaWmB\ngYFiF2KBvLy8wsLC5ORksQuxjFardcUfj+eee06lUoldiAXS09M5jjt69KjYhVigsrLyypUr//73\nv83Z2cPD469//atSqXREJS7WQR87duzEiRNm7hwXF+fIWuygtrb2+PHj9be0bdvW19dXrHqsExwc\nHBUVJXYVFuvZs6fYJVggOzt7x44d2dnZZ8+eFbsWy0RERERERIhdhWU8PDzMT4+tW7dmZ2c7qBIX\n66CJqF+/fs8884w5e5q5m4h+/vnnyMjI9evXi10ISN3BgwdPnTpFRLGxsZ988onY5cD//Pjjj457\ncRfroGVGo9F07NhR7CrABajVauFGRUWFuJWAMyGgxaTVajt16iR2FeACvLy8hBsudx4T2AIBLSaN\nRoOABnMYAhodNFMQ0GIqLi4OCAgQuwpwAQhoNiGgRVNRUeHh4SF2FeAaDDNojDiY4qSA/s9//jNw\n4EA/P7/WrVuvXLmypqbG8NCZM2f69evn7+//8MMPf//9986pRwquXr3aoUMHsasA14AOmk3OCOjV\nq1cnJCSEhoZu2rRpxowZa9euXbFihfCQVqsdMmRIVFTU22+/HRUVNXToUJc738FqWMIB5kNAs8nh\n66ArKireeOONcePGff755xzHEVHXrl2nTZu2bNkyhUKxcePGkSNH7ty5k+O4mTNnjh8/ftOmTZs2\nbXJ0VVKg1WqHDRsmdhXgGjw9Pd3c3PR6fVVVVV1dnUKhELsicAaHd9CZmZmlpaUjR44U0pmIBg0a\nVF5eLlwidu/evQkJCcJDHMclJCS44oUrraPRaGJjY8WuAlyGp6encKOyslLcSsBpHN5BR0dHa7Xa\nFi1aGLacPn3ay8srLCysurpap9N169bN8FD37t11Ol1tba27u+ud4mipgoKC4OBgsasAl6FWq4Uj\nhOXl5d7e3mKXA87g8BxUq9X1D4WdO3du1qxZL7zwgoeHx40bN4iofkg98MADPM/funXLcPL+zp07\nDQNrIrp161bbtm3nzp3r6LIdrbq62kFXV5GluXOpfXtasEDsOkRlbAxdVETx8bRtG7Hw99ilS/Ts\ns3T4MPn7i13Kn65cubJ58+awsDAi6tSp0+jRo+344s5bZldUVDR//vz+/fuPHj165cqVRCRclM4w\n+jBsqa6uNmzp2rXrI/W4ubk57rokzpSamtq+fXuxq3AZ587RvReVYpGxgNbp6Mcf6fJlMWpyutRU\nOn+eJPURWvn5+ZGRkb169erVq5fd12U5aZJw4MCB2bNne3t7JyUlGX7DhISEEFFhYWHz5s2FLYWF\nhUQUHh5ueGKXLl22bt1quPvzzz8L+7g6LOGwiL8/3b0rdhFiM7YUWshtRubSwpkDVVVi11GPUqns\n1avX4MGDHfHizuig9+3bN2bMmPj4+N9//71+/+/p6RkSEnLlyhXDFo1GExoaysLpGzjJ2yJ+fgho\nox20cOyQkdV3wqWwJRXQDuXwgK6srJwzZ87f/va3Dz/80HAY2mDChAk7duww3N25c+fEiRMdXZIU\n4DJJFkFAk/GAFjYzEtBC81ZvCCpzDh9xHDt2rKCgQKFQrFmzpv722bNnBwcHP//883FxcSEhIcOG\nDTty5EhSUtK5c+ccXZIU3Lx5s/4kB0xDQNP9RhyMBDRrHbTDAzotLY2Itm3b1mD7uHHjgoODu3Tp\ncvDgwb///e8ffvhht27dDh061LlzZ0eXJLra2lqcaGARBDQZ76CVSnJ3Z2sGjQ7abl588cUXX3zR\nxA4DBgw4ffq0o8uQlPT09LZt24pdhSvx9aWqKqqqIgYOTxhl4mxvT0900PKEq9mJAEs4LOXnR0Ss\nN9EmLmjn5cVKQEtwFYdDIaBFgCUclhLOSmA8oNFBE3sjDgS0CLCEw1JCB33njth1iMpEQHt5sTKD\nxogDHC4nJ8dwbg6YAyMOul9Ao4OWJQS0s+n1+vpnt4M5ENCEGTQRoYMGR8vMzGzVqpXYVbgYBDRh\nBk1EOEgIjoYlHFZAQBNGHEREpFQSx2HEAQ6DJRxWQEATDhISERHHkVKJgAaHwRIOK3h5kVJJJSVi\n1yEqEzNodkYcROThgREHOEx2dnZUVJTYVbgenO2NDlqgUqGDBsfgeR6rOKzj64uAxgyaiMjDAwEN\njoH22Wq+vhhxmBpxMNVBY8QBDoElHFbz82M9oDHiEGDEAY6CJRxWw4jDREB7eFBdHdXUOL0mMSCg\nwVGwhMNqGHGYPpOQmLlmP2bQ4CiZmZmtW7cWuwqXhBGH6TMJiZnPjUUHDQ7B83xtba2bG95za2DE\ngYAW4CAhOERubi4uYmc1X18qLSWeF7sO8SgUCpVKRUR6vb7q3ohiasSBDhocAks4bOHrS3o9lZWJ\nXYeojDXRTHXQmEGDQ2AJhy2Ey3EwPoY2dpyQqYBGBw0OgSUctvD1JWI+oI110BhxyBUC2nnwYd62\nEAIaxwmFGyyPOBDQ4BDV1dXu7u5iV+GqMOIgjDiICKs4wBHy8/PDw8PFrsKFYcRBGHEQEQ4SgiNg\nCYeNMOIgjDiICCMOcAQs4bCRENClpWLXISpDQLM84sAnqoD9YQmHjTDioHozaJY7aIw4wP6uXr3a\nvn17satwYWo1ubmhg8aIg5RKqqujujqx63AKBLSTVFdXC+fpgnU4jry9We+gjY04VCriOFbWNgj/\njBhpohHQzlBQUBAcHCx2FS5PuBwHy4yNODiOPDxY6aCFgGbk4tcIaGfAANoucElo0xe0YySglUoi\ndNBgR1hjZxc+PghoBDRGHGBv6KDtAiMOYzNoQkDLFALaGVJSUjp06CB2FS4PHbSxGTQReXriIKEM\nIaCdoaKiwlNYCQU2QAeNEQchoMG+iouL/f39xa5CDtBBI6AJqzjAvjCAthd00CY+2JudZXZYxQH2\nhCUc9uLrS2VlTH8sITpowogD7AsdtL34+LD+sYQIaEJAg31ptVp00Hbh40PE9vWSTIw4ENCyhIB2\nuNLSUm9vb7GrkANccRQdNCGgwY5KSkp8hVwBmwkdNAKamgpoDw+21kFjFQfYAeYbdoQRB0YchFUc\nYEdYwmFH6KDRQRM6aLAjLOGwIyGgWV7F4eHh4ebmRkTV1dV1916yHh20LCGgHQsfRWhH6KDJeBPt\n4UF6PRN9JTposJu7d+/6+fmJXYVMIKDJZEATMTHlQAcN9lFWVmb45wS2ExYrIqCFGw2OE7IT0Fhm\nB/aRkpKCI4R25OlJ7u5Mz6AJH+z9ZweNEQfYCks47M7HBx00RhzEcQhosBmWcNgdAtrYiIOdDpqI\nlEqMOMBmGo0mNjZW7CpkBQFtbMTBTgdNREolOmiwWVFRUWBgoNhVyAoCGiMOIlKp0EGDbaqqqjyE\nfzRgPwhow2enVd47zmBtxIEOGmySkpISExMjdhVy4+PD+ioOlbDKjKjm3ohCBy1LCGhHwRIOR/D2\nZr2DNgR09b0RxVQHrVKhgwbbYAmHI2DEoRSWATcKaKY6aKziAFvhKhyOgA4aIw7CDBpsp9PpQkND\nxa5CbjCDNtZBCyMORgIaM2iwSXV1teEfEtiRtzfV1LDyj7NJ6KAJHTTYKC0trV27dmJXIUPC9ZJY\nbqIxgyYENNgISzgcBAFtCOgGHbTQWCOgZQYB7RA4Qugg+FAVYyMOlYo4DgEtNwhoh8AaOwfBJaGN\njTg4jpRKVgIaBwnBJnl5eZGRkWJXIUMYcRjroImlz41FBw3Wq62t5ThO7CrkCQFtrIMmBLQcIaDt\nLyMjo23btmJXIU8IaHTQhIB2kLi4uOTk5Ppb9u7dy91r3LhxzizJEbCEw3FwkNB0B83IZJadGbS7\n077Svn37zp8/32BjWlqar6/vq6++atgig+XDGo2mZ8+eYlchTzhIiA6aWOqgnRHQGzdu3Lp16+XL\nlxs/lJqa2rFjx1deecUJZTiNVqudNm2a2FXIE0Ycxq5mRwhoOXLGiKN169Z/+9vf1q1b1/ihtLS0\n9u3bO6EGZ7p+/XqLFi3ErkKevLzIzY3pgDZ2ogohoOXIGQE9cuTIRYsWLVq0qPFDqamp165d69mz\np1qt7tSp07vvvqvX651QkuMI9WMVh4NwHKnVTAc0Omhi6XrQzptBN1ZWVpabm8vz/Lp169q3b//D\nDz8kJiZWVVUtXry4/m6XLl0y3G7wOWwSlJWV1apVK7GrkDNvb6YD2nQHzcg7I6nrQfM8X1JScvv2\nbSLy9vY2/Aa1CzEDWqFQfPfddx06dAgPDyeiPn36+Pv7v/TSS4sWLTJ0oKtWrVqxYoXhKVVVVT7C\ngXypwhIOR2M8oE130EVFTi9IDEol1dURz5MU/lItLS3917/+tWPHDiLq37//Cy+8YMcXFzOgPT09\nBw0aVH/LqFGjZs6cWf80vKVLly5dutSwQ8+ePQsLC51ZpKVwFQ5HU6upvFzsIsSDZXZEpFQSz1Nt\nLUnhmr6+vr7Lli0bPHiwI15czBNVMjIyDh06VH/oXFdXR0R+fn7iFWUrXIXD0dBBCzcajzhUKlZm\n0O7uRMTEGFrMgNbpdPHx8fv37zds+fzzzzt37izxIYZpWVlZ0dHRYlchZ+ighRuMd9BEVFsrdh2O\nJ+aIo2/fvgMHDpwxY8bixYtbtWp1/vz5LVu2JCUliViSjXie1+v1bm44gd6B1GrKzxe7CPGgg6Y/\nA5qFDlrMgOY47vDhw8uXL9+zZ096enpsbOxXX30VHx8vYkk2un79esuWLcWuQuYYH3GggyYEtIPw\nPN9gi4eHx6pVq1atWuXMMhwHSzicgPERBzpoYimg8ce4PWEJhxOggxZuoINGQINlsITDCdBBCzea\n7KBraqjRn6kyhIAGa6Snp7dp00bsKmROraaKCiZiqEkmOmiVinieiSYay+zAGnV1dQqFQuwqZE6t\nJp4nyZ/z7yimLzdKxERAs7PMDgFtN7m5ufgcQidQq4kYvuKou7u7cCGEmpqaBkfdhehm4TghRhxg\nMSzhcA7hktDMjqE5jjN2vSQhoNnpoBHQYAEs4XAOoYNmNqDJ+AXtWBtxIKDBAljC4RyMjzjI+AXt\nMOKQHwS03aSmpsrgAxWlj/ERB6GDRkCDFaqrq5VSuPqh3GHEgQ4ay+zAMjqdLiwsTOwqmICANrYU\nmrWDhFhmB+bCEUKn8fIiYjugjS2FZieg0UGDZbDGzmmEDprZE1UIHTQ6aLAUlnA4DUYc6KDRQYNl\nrl69GhMTI3YVTEBAG+ughVUcLBwkRAcNlqmsrPQQ/n2Ag3l4kELB9IgDHTQ6aLBAUVFRUFCQ2FUw\nxMsLHTQRZtDooMEcWMLhZF5e6KCJ0EEz0EFb+ZFXFRUVOp1OqVQ+8MADhh8XZmEJh5Mxfs1+dNA4\nk7AJlZWVO3funDp1avPmzdVqdXR0dLNmzTw8PDp27Pjss88eOnRIr9c7rlApwxIOJ2M8oNFBCx00\nRhx/KC8vX758ecuWLT/66KPu3bt/+umn169fLy8vLykpyczMfPvtt8PDwxMTE9u0afPBBx/U1dU5\numipSU5O7tChg9hVMITxEQc6aIw47tGnT59x48b98ssvzZo1a/CQj49PdHT0448/vnz5co1Gs379\n+m+++Wbfvn0OKFW6ysvL1cLiL3AKdNDCDWY7aI4jd3cmOmizAvq7774LDQ29726dOnX68MMPdTqd\nzVW5kjt37vj5+YldBVvUaiopEbsI8Ri7WJJCQQoFE30lESmVTAS0WSMOc9LZup1lAANo52O8gzZ2\nuVEiUqmYOFGFiNzdmfhVhGV2tsISDudjfAZtrIMmIpWKiREHoYMGM6GDdj6cqCLcaLKDZiSgme6g\nW7Zs6d4UJxfnEhDQzqdWo4MmMtJBsxBbRKwcJGw6oJOTk3v06PHll1/W3svJxbmEkpISHx8fsatg\nC+MjDmPL7IhIqWSlg2Z6xKFWq2fMmIHFCfdVWlrqLXxGHjiRlxdVVopdhHgw4iDGO2gimjdv3oAB\nA4jowIEDTqzHxSQnJ+MIofN5eZFez8pyhcaMrYMmjDhk5/4HCUePHu2EOlwUlnCIwtOTiOEPVcGI\ngxDQYA4cIRSF8LGEzAY0OmhifBUHmAkXGhUF4wGNDprQQYM5iouLAwICxK6COYwHtOkOmpGAZnoV\nB5ijoqLCU5iGgnMxHtAmOmimRhwIaCKiqKgoJ9ThilJSUnCVUVEwHtAmOmiMOGTm/icHZmVlGW4n\nJyc33kGlUgUHB/v5+XEcZ8fKpA9LOMTCeECjgyYid3cm1sJbdva2iTxSKBT9+/efP3/+xIkTba7K\nNWi12uHDh4tdBYsYD2h00IRVHE3Kyclp3br11KlTDxw4cPHixYMHD06bNq1Dhw6nTp06cOBA9+7d\nJ02axM7V+jUaTWxsrNhVsIjxgEYHTRhxNGnu3LkDBw7ctm2bYcvjjz8+c+bM7du3f/DBB48//nhw\ncPD69esnTJhg7zqlqKCgIDg4WOwqWMR4QGMVB2EVR5NOnTo1bNiwBhuHDRv29ddfC7cfffTR3377\nzT6lSVtVVRU+zlwsCGjhBtZBy55lAd28efMLFy402HjhwgVDI/n777+HhYXZpzRpS01Nbd++vdhV\nMEpY3MjCMaImmbhYklKJEYesWDbiWLp06fTp0/39/adPnx4ZGZmXl/fJJ59s2LDhP//5T1lZ2X/+\n859//OMfr776qoNqlRQs4RAROmjhBsufqIKAbsKUKVO8vb1ff/31f/7znzzPE1FsbOyuXbsmTZqU\nmpq6ffv2N95449lnn3VMqdKi0WiEq/2B8ymVpFCgg0YHLXYRjmfxh6SMGTNmzJgxVVVVeXl5YWFh\nXl5eRFRWVtauXbuffvrJARVKlFarZeRXkTR5erIb0KZn0CzEFjET0BYfJBRueHh4REdHC+m8f/9+\nBleb3bx5Mzw8XOwq2MVyQJu+WJJeT3V1Tq/J6RQKBHQj8fHxJ06cMNzNysoaM2bM2LFjAwMD7VyX\ntNXW1uITGsXF8oeqmD5RhYiJKYe7OxO/hywL6GnTpo0YMeLYsWPV1dWrVq3q1KnTiRMnNmzYwNRw\ng4jS0tLatm0rdhVMQwdNRjpoYiagWeigLWsD33//fXd391GjRrVo0SI1NXXKlCnr16+PiIhwUHGS\nhSUcovP0xCoO1jtoBHRDbm5umzZtUiqVGzZsWLdu3aJFixxUlsRpNJoHH3xQ7CqYhg6a0EEjoAUN\nLmI3Z86c7OzsFStWdO7cOTo6WtjI1IU3tVrtU089JXYVTENAEzpoBLTA2J/z8fHxhtvCsmhG5OTk\nNG/eXOwqmMbyiMPNzc3d3b22tpbn+QbHq4WAZuFcFWEVB8+TvC9ybNZBQt4Mji5UOurq6li78rUE\nsbyKg4xPOYTpNCMdNBHp9WLX4WD4yCuLZWZmtm7dWuwqWMdyB03GjxOyM+JQKIhI/lMOBLTFsIRD\nClieQZPxDpqdgBY6aAQ0NKTRaDp16iR2FaxjPKDRQSOgoWlarRYBLToEtHADHbS8IaAtlp2djU86\nFx3jAW1spR0CWmYQ0JbR6/V6vR6rOETHeECjgxYCWvaX47DmRJUmMXKiyrVr1wzn5oCIPDxIr6ea\nGvqzlWSLsYOEQmyxENCMrOKw6USV+hhZCo0lHBJh+NQrNgPa9EFC2ccWoYOuj5HwNQeWcEiEhwcR\nUVUV+fqKXYoYsMyOkQ7ayhl0VVVVVlZWBXunCmAJh0Qw/rmxWGbHSAdtcUDv2bOnW7duXl5erVq1\nUqvVXbp0+eKLLxxRmTThNEKJMHTQbMIMGh10E3bt2vWXv/xl7NixV69eLS8vT0tLGzdu3KRJk/bs\n2eOg+iSF5/m6ujo3Nyx9ER86aOEG4zNo2X+nll0Peu3atfPnz3/99deFu23atFmxYsXt27fXrFnz\n5JNPOqA8ablx40aLFi3ErgKI0EFjBq0gwoijgatXr8bFxTXY+PDDD6ekpNivJOnCEg7pQAct3GjQ\nQbMz4mCkg7YsoNu1a3f27NkGG8+ePdu+fXv7lSRdWMIhHeighRtNdtCyjy1CB92kl19+edOmTa+9\n9lpGRkZlZWVGRsayZcs2b968ZMkSc54eFxfX+JyXM2fO9OvXz9/f/+GHH/7+++8tqsfJsIRDOtBB\nCzcYX8Uh+19FlgX05MmTd+/enZSU1KZNGy8vrzZt2uzdu1c4cnjf5+7bt+/8+fMNNmq12iFDhkRF\nRb399ttRUVFDhw4156xFsaSnp+PDvCUCHbRwg9lVHIwss7PsICERPfnkk08++WRFRUV+fn5YWJiX\nl9d9n7Jx48atW7devny5yYdGjhy5c+dOjuNmzpw5fvz4TZs2bdq0ydKqnKO6urr+xwuBiNBBCzeY\nXcWBZXameHl5RUdHm5PORNS6deu//e1v69ata/zQ3r17ExIShGsPcRyXkJAg2VXVeXl5kZGRYlcB\nf0AHLdxo0EG7uZGbGzpo+bC4H6yqqsrOzq5r9MaYuFjSyJEjR44cSUSLFy+uv726ulqn03Xr1s2w\npXv37jqdrsHnYEoElnBICjpo4UaTH+zNQkAz0kFbloO7d+9OSEgoLy9v/JAV1+vQ6XREFBwcbNjy\nwAMP8Dx/69atiIgIYcv69etXrlxp2KGkpMTHx8fSL2QXWMIhKYx30MYuN0rMBLR0OujS0tJ//etf\nu3btIqLevXvPnj3bji9uWUAvWrRo6NChW7duDQsLs/1rC5le/9rKwpb6P3NTpkypP1hYunRp45bB\nObRa7eDBg0X50tCY0EEzG9DGLthPRO7u8u8rSUodtJeX15QpU4QTRHztfe0uywK6pKRk/vz5dkln\nIgoJCSGiwsLC5s2bC1sKCwuJKDw83LBPZGTklClTDHfXrVsn7ON8V69eZWS5t0sQOkgENOMBLYUO\nWqFQREREOOgSPZYdJHzwwQebXIxhHU9Pz5CQkCtXrhi2aDSa0NBQD+HPV4mprq42/F0JonN3J4WC\n3YA2PeJAQMuGZR30+++/P3jwYHd390GDBjVIK+s+UWXChAk7dux47LHHhLs7d+6cOHGiFa/jaLdu\n3RL6fZAOlYoapRMrjK3iIGY6aOnMoB3KsoC+fPnyrVu3XnzxxcYPWXdR/+effz4uLi4kJGTYsGFH\njhxJSko6d+6cFa/jaDhCKEEeHuigMeIQuw4Hs/hU78ceeyw/P59vxLov36VLl4MHD/70009Tpky5\nePHioUOHOnfubN1LNVBRQcuXU1mZXV4Ma+ykiOWARgfNSEBb1kEXFRXNmzcvNDTUui/WZI4PGDDg\n9OnT1r2gCbm59Prr5OVFL79sh1fTarUJCQl2eCGwH5YDGh00IwFtWQc9bNiw3377zUGl2FebNjRy\nJP3rX2SXj+VKSUlh5GPLXQhm0MR8By3779SyDnrBggUzZ86040FCh/q//6PBg2nHDrJ95XhFRYWn\nsPIWJAMdNDHcQeMgYROExdh2PEjoUIMGUd++9NZbNGsW2fIxVUVFRYGBgfarC+yD5YA23UGzcCYh\nRhxNaHxs0MaDhI62cCFdvUr799v0IrgMtLdOvBkAACAASURBVDSxHNDooBHQ5jp27FjjaydJxBNP\nUJs21NR19CyAJRzShBk04UQViQaP3dghoIOCgqR5dgkRKRQ0fz798AP98IP1L4IOWprQQRM6aLkH\ntGUz6LS0tEOHDtXfUl1dffToUUesk7OXp5+mFSto3TpKSrLyFbRaLTpoCfLwoIICsYsQiekZNAtX\nYcUqjiZcunRp/vz5DTa6ubktXLjQfiXZmbc3zZlDq1aRVkvWxWxZWZm3t7e96wJboYMmdNBy76At\nG3GMHj361q1bhYWFRUVFt2/fLi4uLigomDZt2rx58xxUn13MnUtKJW3das1z79y5Y/dLCIJdeHiw\nO4M2cbEkRgKa48jNTf4BbVkH7eHh0fhSc8uXL588ebI0r6EhaNaMNm4k685/TE5OxgBamnCQkBju\noIlIoUBAm6GysvL8+fM6nc7qU8Cd4JlnrHwilnBIlkqFEUcTHbRCgYCWD8sC+vfffxc+2cXgxo0b\nBw4caN68eWBgYGlpaU5OTrt27dxsOS1EYrRa7YQJE8SuAprA8ojDdAct+9gSKBSk14tdhINZFtCp\nqanvvPNO4+2hoaFdunS5evUqz/Nr1qxZsmSJncoTn0ajSUxMFLsKaALLIw500ESYQTcyfvz40tLS\nxtt5ntfpdL/++mt2dvbYsWPtVJsk3L1718/PT+wqoAksB7TpZXayjy0BRhzm4jguLCzs0Ucftcur\nSUdZWZlarRa7Cmia8PHVPE/1PnaYFVhmR2wEtFnD4osXL5r5cjzPm7+z8733HuXkWLB/cnIyjhBK\nlkpFPM/EhYEaM9FBszPiQED/YcGCBU8++eTp06f1xmfyVVVVSUlJjzzyyLZt2+xXnp3V1FBmJn38\nMf373/TRR/cPayzhkDKhiWRzyqFQKIRD8Xq9vsGVcDDikBOzRhwnT57ctWvXCy+8cPv27eHDh8fF\nxbVq1SooKKiurq6wsDA5Ofns2bNHjx7t2bPnm2++OWDAAEcXbbVWrSgqivr3JyI6eZJWr6Zx42jo\nUKP7a7XaESNGOK08sAjLAU1EKpWqsrKSiGpqahTCeXVEhA5aXswKaI7jpkyZMnny5J9//vm///3v\np59+mpqaqtPplEplaGho586dBw4c+MYbb7Ru3drR5dpoxAjatIkCAojn6aGHaOPG++yv0WgWL17s\nlNLAYowHtFKpFAK6urq6/qdJoIOWEwsOEnIc17t37969ezuuGoc6doyyskijoQ0byMvLrKfgUv1S\nxnhAGztOyNRBQqyDlgmep5gYGjqU+valX3+luLj7P6WioqLxee0gHcJ/HGYD2thxQnZGHFgHLR8c\nR82bExF17mzuU65evSrBD1oEA6GDxNnejTto2ceWgIURh3zOybY7LOGQOMZHHMbO9hZiS6ofQmdP\nCGi5OX6cvv7a3J3xQSoSx3hAG1ZuNF5mRwxcKJkQ0PKjUFDjqfLZs3T1ahM7azQaBLSUCQHN5okq\nRGS4JFmDsxMY+agRQkCb78knn7TL6zjawIHU+HR0Nzc6cYLef7/hdolfQBUY76CNBTRTHbTsV3HY\nJ6AbXxBAso4fJ53uni19+9Izz5BSSfUvA1VdXW2Y8YE0Cf99ENBNdtAsBDQLqzjsE9BfffUVdy+7\nvKwjqFT03nv0zjsNxxqhoaTR/O9uampq+/btnVwbWAQjDuFGgxk0OwGNDtpcY8eO5e9ll5d1hP79\nKSKCXnyRMjLovfdozx4qLaWaGjp3jnr0+N9uWMIhfUIHzWxAGw4SMjuDZqGDts86aNf6CBWViioq\n6PHHiYgyM2nvXqqspOeeo/ojDY1G069fP7EqBHNgxCHcwAxaxuwT0Pv27bPL6zjHmDH0zTf0xBNE\nRK1aUatWTeyj1Wpnz57t5MLAIhhxCDcwg5YxOwQ0z/NFRUUZGRmZmZmZmZkvv/yy7a/pUMHBFBBw\nn33y8vIiIyOdUg5YifERBwKahZParQzojRs3CnEs5HJJSUl4eHhkZGSzZs3sW5+DDBtm6tHa2lrX\nGtqwCSMO4UaTBwlln1yEDtqEF154YeDAgU899VTLli1btmzZrl27vLw8+1YmovT09LZt24pdBdwH\n4yMOYwcJMYOWEyv7xG+++ebOnTtJSUktWrSQX5ZhCYdLwIhDuMHsiAMdtFHx8fGPPfbY7t27x44d\nGx8fb9+aRHH7Nhmu/KzRaHr16iVqOXB/GHEIN1gOaHTQxp/p5jZlypRff/21bdu24eHhiYmJugan\n6LmOvDwKCaEvv/zjLi6T5BLQQQs3mF0HjWtx3J9SqXzuuefS0tK8vb27d+8+b948u5TlZOHhFB1N\nmzf/cff69estWrQQtSK4P4WC3NwQ0OyeSYgO2lze3t5Lly69cuWKj4+PXV7QyTiOpk+nY8coPZ30\ner3Ez1YHA6WS3YDGQUJ00A2lpqY+8MADxh4NDAxcvXq1zSWJY9YscnOjjz+mzMzMVk2euwLSo1Jh\nBo0ZtJxZFtB1dXWFhYVE9PPPP8fFxQUGBo4ZMyYjI8MxtTlV8+bUsye98w5duYIlHC6D5Q4aM2h0\n0E3TaDRDhgy5cOGCWq0+cOBAr169srOz7V6Z8w0fTnfv0pEjt3GE0FUgoAkdtKxZE9Bz585du3bt\nnTt3cnJybty40adPn3/84x92r8z5liwhlYpOngxGQLsKlgPa2EdeCZtln1zExokqlq2DFj5I+JVX\nXnlcuBYcUWRk5L///e/BgwfbvzSn8/engQPp+++7RUfjKhyuQalk4m/5JqGDRgfdUEBAABE98sgj\n97yEm5swmJaBUaP4iormaWm4EIdrcHdnt4NGQCOgGwoMDIyIiPj+++/rb3z33Xe7dOli16pE07Vr\nHhGdPSt2HWAelkccCGgEdEMcx02cOPG55547efJkaWnpb7/99vTTT69fv/65555zUH1OdvfuZY4j\nWRzyZAJGHNQooIXNjAS07L9Ni6/FkZiY+OWXXw4aNMiwZeHChdOmTbNnUeJJTf3d13dwbq5K7ELA\nLCyPOEwfJJR9chEbHbTFAR0WFnb69OnXXnvtypUrHTt2nDlz5tChQx1RmSi0Wm1kJJ+bK3YdYB6M\nOMjIiEP2yUVYxWFMdHT0xx9/bPdSpCA9PT06WpWTI3YdYB6MOAgzaFnDcoX/4Xm+tra2eXMOHbSr\nYHnEgYBGQLMlNze3efPmzZpRfj67fZlrwYiDENCyhoD+H+GDVCIjSa+nmzfFrgbMwPKIAwcJWVjF\ngYD+H41G06lTJ+FjbzHlcAkYcRA6aFlDQP+P8EEqQkDjOKFLwIiDGF4HjVUcbElNTW3Xrt3t20QI\naBfB8ogDy+zQQbOlurpaqVSGhJBKhRGHa8CIgzDikDUE9B/y8/PDw8OJiOMoLAwHCV2Duzu7HbSx\ng4TsjDgQ0AwRjhAKt8PDKT9f3HLALBhxEDpoWUNA/0FYYyfcDg1FQLsGljto0wcJZZ9cRMRxxPPE\n82LX4UgI6D8ISziE22FhCGjXgIAmtg8Skty/UwT0H65evRoTEyPcDgsjnU7mv5nlAQFNDC+zQ0Az\npLKy0sPDQ7gdGkrV1XTnjrgVwf2xHNCmDxLKO7YEwncq70YKAU1EVFBQEBwcbLgbFkZEmHK4AIWC\n3YDGiIPjiOT+nSKgie4dQBNRaCgRAtoVuLsz8bd8kzghnzDiQEDLXv0lHPRnB63TiVYPmEmhoLo6\nmf+Ra4yhg+bv/f5ZiC0BC98pApqoUQeNEYercHcnIkanHMZm0BxHHCfz2BIgoB1u79693L3GjRvn\n/DKSk5M7dOhguBscTO7u6KBdAMsBbWwGTWxcRYjYCGiRL5aUlpbm6+v76quvGra0a9fO+WWUl5er\n1WrDXTc3Cg5GB+0CENDUVECzcKFkQkA7QWpqaseOHV955RURayguLvb392+wEeequAQENKGDlvV3\nKvKIIy0trX379uLW0GAALQgNxYjDBSCgyUgHLe/YEmAdtMOlpqZeu3atZ8+earW6U6dO7777buOf\nNkdrsIRDgA7aJbAc0MYOEhIzIw4W1kGLOeIoKyvLzc3leX7dunXt27f/4YcfEhMTq6qqFi9ebNin\nqKjo8OHD9e/aPcG1Wu2kSZMabBTO9gaJYzmgMeKQyIhDr9fn5eVlZGQQUUBAQFBQkB1fXMyAVigU\n3333XYcOHYQLMffp08ff3/+ll15atGiRYRH+tm3b3njjDcNTSkpKfHx87FtGkyOOBx6gkhKqqCAv\nL/t+NbAnIaBZ6BYbw4hDIgFdXl7+2WefnTx5koh69er1zDPP2PHFxQxoT0/PQYMG1d8yatSomTNn\n5uXlRUZGClsWLVq0aNEiww49e/YsLCy0bxl3795tHPrCid9FRSR8RCFIEzvXPm4MAS2RGbSPj89L\nL700ePBgR7y4mDPojIyMQ4cO1f/xEqZpfn5+TquhpKTE19e38XbhzxR7/y4AO0NAE8MBzcIMWsyA\n1ul08fHx+/fvN2z5/PPPO3fubPchhglarbbxEUKq10GDlAkBzeYM2vRBQnnHlkAiHbRDiTni6Nu3\n78CBA2fMmLF48eJWrVqdP39+y5YtSUlJzqyhySUchA7aRWAGTeigZf2dihnQHMcdPnx4+fLle/bs\nSU9Pj42N/eqrr+Lj451Zg1arHTt2bOPtQkCjg5Y4jDiI4YBmoYMWeR20h4fHqlWrfv3119LS0vPn\nzzs5nenez4qtTxhxoIOWOAQ0MRzQLHTQrF/Nrri4OCAgoPF2tZo8Pen2bedXBBZAQFNTM2imrmaH\nDlq2ysrKvIyvcw4KQgctdSwHtOEgITpoGWM6oFNSUpo8QigIDsYMWupYXsVhesQh775SgA5a5owt\n4RCgg5Y+ljtozKDRQctckyd5G6CDlj4ssyOGAxodtMxpNJrY2FhjjwYFIaClDh00MXyiCjpomSso\nKAgW1tM1JTgYIw6pYzmgcZAQHbScVVVVeXh4mNghKIiqqqiszGkVgcVYDmiMONBBy1lKSkpMTIyJ\nHXAyofRhFQcxHNDooOXM2DmEBoGBRIRzVSQNHTQ1FdCMnKiCDlrOTC/hICLhg2Tv3HFSPWAFIaDl\n/U/UGBwkFAIaHbQ83beDFgL67l0n1QNWEDKKzQ7axEFCjpN5bAkQ0HKm0+nCwsJM7CB8bAA6aClD\nB02YQSOg5ae6utrd/T6XWkUHLX0sd9CmZ9Dyji0BOmjZSktLa9eunel90EFLHzpoYriDRkDLlumr\ncAjUanJ3p5IS51QE1kAHTUYuNyrv2BIgoGXrvkcIBX5+6KAljeUOGmcSYpmdbN13jZ3A3x8zaElD\nB00Mz6BxkFC2cnJymjVrdt/d/PwQ0JLGcgeN60FjxCFPtbW1hh9u0zDikDh00MT8mYQIaLnJyMho\n27atOXtixCFx6KDJyJmE8o4tAQJansxZwiFABy1xLHfQps8kZOGXFgJansxcwkHooCUPHTQxP4OW\n9399FgPazCUchIOEksdyB40ZNFZxyNP169dbtmxpzp5+flRZSVVVjq4IrIQOmpjvoOX9nTIX0Hq9\nnud5Tvhvez++vkSEkwmlS8goxgO6yTMJWXhPENAylJWV1apVKzN39vYmIiovd2A9YAuWAxqXG0VA\ny5D5SziISK0mInwsoXSxcJjIGJxJiICWIfOXcNCfAY0OWrJY7qAR0AJ5f6fMBbT5SzgIIw5XwMgB\nscYQ0OYdSHJtzAV0Zmam+TNodNDSx8iV2xrD5UYx4pAbnuf1er2ZF+IgBLQrYGTFQmM4SIiAlhvz\nV0ALcJBQ+jDiwDpoGWMroC1awkHooF0BRhzMnkkoQEDLh0VLOAgHCV0BAprxEYe8sRXQFi3hIHTQ\nrgABjYOEMsZWQKenp7dp08b8/VUqcndvOqDz8/PPnj3b5LOKi4tPnDhhVYEu4+uvv26cC/d9yBGY\n+nO+PhwkREDLTW1trbu7u0VPUaubDmi1Wp2YmHjs2LEG24uLi8ePH2/mtT5c1/Xr1xMSEhoH8Tvv\nvPPll1+av07GdowcEGsMIw4WAtqytHJpubm5kZGRlj7LWED7+vp+/fXXEyZMKC8vHz16tLCxoKBg\nwoQJa9aseeihh2ysVuLmzp3r6en517/+defOnYbfeStXrrx58+YHH3zgzN9PGHEgoGWMoYC2dAmH\nwNvb6DI7tVq9b9++CRMmENHo0aOFdF69erXs01mQkJBARFOnThUy+o033rh58+bGjRud/NcDAprZ\ngBbI+ztlK6AtOkIoMNZB//noHxldUlKyZcuW1atXP/zwwzZV6VIMGd25c+f8/HznpzMhoJk/SChv\nDM2gLV3CIfDyoooKUzuo1eqtW7fOmzdv9OjRTKWzICEhoaam5pNPPnn77bdFmbzjICGzHTQLIw6G\nAjo1NbV9+/aWPsvT8z6fqFJYWDh16tQ9e/Z8++23Bw4csL4+1/Tmm29GREQsWbJk+vTptbW1zi+A\nkTBqDCMOgby/U4ZGHFVVVUql0tJnqVRUWWn00cLCwvHjxwuTjbi4OGEePWrUKFvqdCFvvvlmTk7O\npk2bOI7jOG7atGk7duywdJ2MjZgKowbc3NyEdG5whRkW/vYnNr5NVjponU4XFhZmxRNVKqqubvoh\nIZ1XrVolTDa8vb337du3ceNGRvro+ulMRLNmzRo2bNi0adOc3EczHtDCjcZNNDvviby/U1YC2roj\nhGQ8oEtKSiZNmrRs2bJ+/foZNnp7e+/evXv9+vWnTp2yulSX8N5776WkpLz77rv1585PP/10XFzc\nnDlznFkJApoaBTQj7wk6aPmwbo0dEalUTc+gy8vLV69ePXTo0AbbAwICkpKSnHmmhihatWq1bds2\nw3EqgwULFjzxxBNOPpOQhTBqkrGFHCwkl4G8/+vLPEcMrFvCQcY76LCwsD59+jT5lICAgP79+1vx\ntVxIfHx843S+70OOwHJAm1jIwcJ7wsLvIVYCOiUlJSYmxoonengYnUGDFLAc0MYw9Z7I+ztlJaAr\nKys9PT2teKKJg4QgBUyFkZlYaC2JjW+TiYAuKioKCgqy7rkIaIlj9kxC09j5pSXv75SJgLZ6CQch\noCWP5Q7asISGv/ctYOQ9QQctE1Yv4SDjqzhAIhgJI4uwkFwG8v6vz0RAW72Eg4iUSqqtxR/R0sVy\nQBvroEnusSVg4fcQEwGdnJxsdQft4UFEmHJIF8sBbQxT74m8v1MmArqsrEwtfLyg5VQqIgS0hDEV\nRmZiobUkNr5N+Qf0nTt3/Pz8rH46AlriWA5oxkccAnl/p/IPaFsG0ISAljyWA9oYRt4TdNByYMsS\nDiISrnaAg4SSxcK/UmNMLLNjh7x/Fck/oG3soIWAduLFfwAA/sBEQKODljd591DAMvkH9N27d205\nSChcLwwBLVlM/TnfgImDhCAPMg/o0tJSb29vW14BHbT0IZ0aY+E9YeF3s8wD2pZTVAQIaIlj4V+p\npZh6T+T9q0jmAW3jEg5CQIOEYcQhezIPaBuXcBAC2hUgnUCuZB7QGo0mNjbWlldAQEscU3/ON4AO\nWvZkHtC3b98OCAiw5RUQ0AAgFvED+syZM/369fP393/44Ye///57O75yRUWFdR9zVR8CWvrQPjbG\nwnvCwh9PIge0VqsdMmRIVFTU22+/HRUVNXTo0OTkZHu9eEpKSocOHWx8EQS0xLHwr9QYnOpNcv9V\nJHJAb9y4ceTIkTt37kxISPj000/j4+M3bdpkrxe3fQkHIaBdgbz/iQLLRA7ovXv3JiQkCI0Ax3EJ\nCQlffPGFvV7c9iUchICWPKa6RWCNmAFdXV2t0+m6detm2NK9e3edTldbW2uX17fls2INENAgWVjF\nIXtiBrROpyOi4OBgw5YHHniA5/lbt24ZtixfvlxVz6VLl+7cuWPm65eXl4eEhNhYpFJJhDZNwhSK\nP36JMsjtz+/c7d63wM2NifdEuE6O8P8iKikpWbJkyfDhw4cPH7569Wo7vzovnuzsbCIqLy83bCkt\nLSWirKwsY0/ZvXv35s2bnVLdHyoq+I8+4quqnPk1wQKnTvGXLoldhEjmzJlDRAMGDGiwPSWFP3xY\nlIqcSq/nd+zgb98WuYxZs2alpaU56MXF/D0rtLeFhYWGLcLt8PBw0WpqxNOTZsz443NVQIIeeYTq\nDcnYsmXLlqysrOPHjzfY3r49PfqoKBU5FcfR1Klk23kOUidmQHt6eoaEhFy5csWwRaPRhIaGegif\npA0A9xMVFaUQ/Y98cBiRJ1UTJkzYsWOH4e7OnTsnTpwoYj0AANLhLu6Xf/755+Pi4kJCQoYNG3bk\nyJGkpKRz586JWxIAgESI3EF36dLl4MGDP/3005QpUy5evHjo0KHOnTuLWxIAgESI3EET0YABA06f\nPi12FQAAksPAakkAANeEgAYAkCgENACARCGgAQAkCgENACBRHO9S18E6cuTIvHnz/Pz8zNlZq9U6\nuh67q62tVSgUnEtdnEmv1/M873Lns9XU1CiFS2G5Dlesua6ujohc7sdDpVK1adPGnD3v3r174sSJ\nyMhIh9ThoGt8SMGgQYPELsFiL7744iVXu/bPkSNHVq5cKXYVFnPFH48JEyYUFhaKXYVlPvzww+3b\nt4tdhWXy8/MnTZokdhU8L+7FkgAAwAQENACARCGgAQAkSvxTvR0nLCxM7BIs5ufn5+npKXYVllGr\n1WYetpUUV/zxCAoKcrmDhD4+Pu7uLpYzSqUyMDBQ7CqIXG4VBwAAOzDiAACQKAQ0AIBEIaABACQK\nAQ0AIFEIaAAAiUJAAwBIFAIaAECiENAAABLlGgHN8/z69eujo6OjoqLWrVtn7OQaE7tZ95Bka66u\nrhauSmrwyiuvSKFmg759+yYnJ1vx4pKq2UHvs73KPn78+JAhQwICAqKjo1esWFFbW2vRi0uqZsn+\nSGs0msceeywwMDAiIuKZZ565e/eupS9u6zcgfbt27fLx8dm2bduHH37o7e29a9cuS3ez7iHJ1pyS\nkkJEGzdu/ORPFy9elELNPM/X1NRs3ryZiLRarRUvLqmaHfQ+26XskydPuru7T5ky5dChQ1u2bAkK\nClq1apVFLy6pmqX5I11aWtq8efN+/folJSXt27cvJiZm9uzZlr64LVwjoAcNGrR27Vrh9po1awYP\nHmzpbtY9JNmav/nmG6VSWVtba5dS7Vjzp59+ariWSIOwk+z7bKJmB73Pdil78uTJXbt2raurE+5+\n/PHHbdq0sejFJVWzNH+kv/rqKyLKz88X7h48eDAsLMzSF7eFCwS0Xq9XqVSXL18W7v72228qlUr4\nFA8zd7PuIcnWzPP8xo0b27VrV1tbm52dXVFRYWO19qqZ5/nCwsLLly9fvny5QdhJ9n02UTPvmPfZ\nXmVPmzbtn//8p2Hn7777ztvb2/wXl1TNvFR/pI8fP75w4cL6BUdERFj04jZygYAuKioiotu3bwt3\nCwsLiai4uNj83ax7SLI18zy/YMGC8PDwZs2aEZG7u/vTTz9dXl4ues31NQg7yb7PJmrmHfM+271s\nnufr6uqmTJkSHx9v0bOkUzMv+R/p33//fe/evb179/6///s/S79TW7jAQULhjfDx8RHu+vr6EpHw\ndpi5m3UPSbZmItLpdEFBQZ999llJScmRI0cOHz68YsUK0Wu2/cUlVTM55n22e9k3btwYP3780aNH\n33nnHfOfJamaSfI/0lOnTp04cWJGRsaCBQssenFb2TfvHaGgoIAa/aYqKCgwfzfrHpJszY3t2LGj\nVatWotdcH93bjUr2fTZRc2N2eZ/tWHZtbe1bb73l5eX1+OOPZ2dnW/Tikqq5Man9SOv1ep1ON3fu\n3JiYGIte3EYu0EEHBgYqFIrc3Fzhbl5enru7e+PLaZvYzbqHJFtzYzExMbdu3RK9ZttfXFI1N2aX\n99n8ekzvVlNTM27cuNWrV+/Zs+fgwYMtWrSw6MUlVXNjEvmRPnfu3NmzZ4mI47iQkJA33ngjJSVF\np9OZ/+I2coGAdnNzGzBgwLfffivcPXr06IABA9zcGlZuYjfrHpJszXq9vmvXrsuWLTO8zvnz52Nj\nY0Wv2fYXl1TNDnqf7VX2unXrzp8/f+HChZEjR1rx4pKqWbI/0keOHJk1a5ZerxceunnzJsdxwljD\nQe9zQ/ZtyB1k165dQUFBBw4c+PLLLwMCAuqvN5wxY0Zubu59d7PuIcnWvGLFCqVS+Y9//OPQoUNr\n16718fH56quvpFCzATW1Dlqa77OJmh30Ptul7DZt2gwbNuyTe5n/zUqtZmn+SP/yyy9ubm7Tp08/\nfPjwvn37evTo8cQTT5jz4vbiGgHN8/zatWujoqKio6Pfeuut+tsb/KMytpvVD0mzZr1ev2HDhu7d\nu3t7e3fv3v2LL76QTs1N7m/ms6RWs+PeZxvLrqys5DjORL8lwbfadM2S/ZH++uuvhw0bFhQUFBER\nMWfOHMPQ2fSz7AWfSQgAIFEuMIMGAGATAhoAQKIQ0AAAEoWABgCQKAQ0AIBEIaABACQKAQ0AIFEI\naAAAiUJAAwBIFAIaXAbP8wMGDLh48aJwd9OmTbdv3zbniebvyXFcg8+NdQRzvopOp4uNjS0vL3d0\nMSBlCGhwGbt27QoLC+vZs6dwd/78+fn5+eY80fw9pSM0NHTcuHHr1q0TuxAQEwIaXAPP8+vWrZs/\nf77YhTjP3LlzN2/eXFlZKXYhIBqXD+isLPLzo4KCph/leerZ855Hf/uNHnqIvL3poYfo11+dUyNb\neJ4vKirKzMxM/lNqaqrtL/vzzz8XFRX1799fuCtcF61jx471L5C2e/fuXr16+fr69u7de/fu3cb2\nLCwsXLhwYWxsrFqtDg0NHTNmTEpKyn0LyMjIePLJJ0NDQ319fQcPHvzLL78YHtq/f/+DDz7o6+v7\n4IMPbt++3dKvotVq4+PjAwMDw8PDJ0+ebJh+NGvWrFu3bvv27bPkfQJ5ccQl8pymtpZ/5BGeiL91\nq4lHP/+c/8tfGj7aoQOfmMjrdHxiK5+joQAABxJJREFUIt+xo9MqZcLNmzfnz58fERHR4GcsODjY\n9hdfvnz5008/bbgrhP6hQ4dSU1OFLZ9//jnHcfPmzdu5c+e8efM4jtuzZ0/jPevq6mJjY729vRcu\nXPjJJ5+sXr26e/fuPXr0EF6EjHzqlV6v79Chw8CBAzds2LBly5bhw4e3aNGirq6O5/n//ve/CoVi\nwYIFn3322bPPPktEGzduNP+rZGVl+fr6zpgxY8eOHW+99Vbz5s39/PwMj65du3bq1Km2v3vgolw7\noNes4Z96ymhAP/ccP2fOPY9WVPBEfF4ez/N8Xh5PxNvp892Bz8rKatGiRatWrQ4ePFhSUnLw4MFB\ngwZpNBp7vf6QIUO2bt1af0uDMO3cuXNiYqLhbmJiYteuXRvveevWrX79+n344YeGPX/66SeO45p8\nTYPr168TkeH67qWlpbNmzRIuDdynT59ly5YZ9hw9enS3bt3M/yqzZ89+6aWXDLvl5OQEBARMnDhR\nuHvy5MnmzZvf770B2XLhgL50iW/Xjr9zx2hACxo82rcvv3Qpn5fHv/oq/+CDTiiTFdOnTyeiEydO\nGLa8//77/v7+P//8s4lnKRQKE3fra9GixbFjx+pvqR9zwmqHtLQ0w6Opqakcx1VWVvLGY7egoODb\nb7+dPHmy4U9JY3tWVlaGhoYOGjRo3759RUVF9bdzHFf/65aUlOTn55v/VWJiYnbu3KmtZ+jQoeHh\n4cKjwi+G8vJyY28LyJurBnRFBd+1Ky+kgUUB/csvPNEf/7t40eF1skNYXFFSUmLYUlpaqlAo+vbt\na+JZ5ge0h4fHL7/8Un9L/ZjLzMwkoop6fxAJkZ2Zmck3it1r166NHz8+IiIiJCRk2LBhwq+Wxq/Z\ngEajmTRpko+PD8dxAwYMOHv2bJNf19Kv4unp2XjwqFKp6n8XWVlZxt4WkDdXPUj42ms0fDgNHGjx\nE2fNopdfptxcWrKEZs1yQGWs6tixIxFlZGQYtnh7e4eFhZlz/M0cSqWyrq7O2KNhYWFEZPiIZcNt\nYXt9PM8PGTIkPDz80KFD+fn5R48eXbp0qTkFdOzYcffu3UVFRWfOnAkJCRkxYoTQVhOR8DHPgmvX\nrh06dEiv15v5VUJCQs6fP9/gn2VVVZXwqPAtK5VKcyoEGRLjt4Id9Ov3v0ZY+N/77ze9p7EZdG4u\nZtD2JHy8cf3jeLm5uRzHjR49WrhLRG+++WZYWFj//v2Fxpa3pINu3br14cOH62+he7vd2NjY1157\nzXA3MTGxc+fOjffMyckhouLiYsOeH3/8Md2vg9bpdC1btjQ8S0jknJwcnuc7der05ptvGvacPHly\nx44dzf8qY8eOfeGFFwy7ZWdnx8bGLly4ULgrdOhVVVXG3haQN1cN6PrqR3BSEl9ZafRRnuc7dvxj\nFccrr/CdOjmvSBYsXLiQiBYvXpydnZ2VlfXoo4/6+/unp6cLjxLRsmXLamtrlyxZMm7cOMNGMzuG\nMWPGbNiwof4WpVL5z3/+8+TJk8Ld3bt3cxz3wgsvNFjF0WDPmpqa0NDQJ5544osvvti/f//zzz8v\nLDvZv39/TU2NsYDW6/UxMTEDBw784IMPtm7dOmTIkE6dOun1euHrKhSKl156adeuXcIy7S1btpj/\nVS5evKhQKCZNmvTRRx+tXLmyXbt2Pj4+ly5dEh49ePBghw4drP0PAi5PbgHdeB7dYMulS/yDD/Je\nXnzfvvyvvzqvSEZs3779oYce8vHxadmy5fTp069du2Z4iIiEsE5JSQkKChI2mt9Br1+/fvLkyfW3\nvPzyyz4+PoGBgYYtn332WY8ePby9vXv06PH5558b2/PMmTO9evXy8vLq0KHDwoULi4uLH3vssYCA\ngDt37hgLaJ7nf//99+HDhwcEBPj7+48aNcqwvI/n+U8//bR79+5qtbpjx45btmwRgtv8r3LmzJn+\n/fv7+PiEh4dPmDCh/qh92bJlzzzzjLH3BGQPn+oNTsJxXHp6euvWrdPT03v16lVcXExE7u7utbW1\nhn0a3K0vPT09Li4uJydHpVI5qWIJ6N2798qVKx9//HGxCwFxuOpBQnBFH3/8sV6v/+CDDwYMGGDp\nc9u0adO9e/dvvvnGEYVJk1arLSwsfPTRR8UuBESDgAbnKS8vDw8PP3Xq1IYNG6x4emJi4vr169n5\nm++tt9569dVX3dzwj5RdGHGAkwjn0dn4IvPmzUtISDBc0E7GdDrds88+u2fPHoVCIXYtIBoENDiJ\nXQIagCn46wmc5Ouvvxa7BAAXg6YGAECi0EEDAEgUAhoAQKIQ0AAAEoWABgCQKAQ0AIBEIaABACQK\nAQ0AIFEIaAAAiUJAAwBIFAIaAECiENAAABKFgAYAkCgENACARCGgAQAkCgENACBRCGgAAIlCQAMA\nSBQCGgBAohDQAAAS9f8BOIcM9dVRsJcAAAAASUVORK5CYII=\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "\n",
    "  options(repr.plot.width=7.0, repr.plot.height=7)\n",
    "\n",
    "\n",
    "  dim(rCO3) <- c(length(var2_e), length(var1_e))  #my one\n",
    "\n",
    "    \n",
    "#   par(mfrow=c(1,1), mar=c(4,3.7,1,1.5)+0.1, oma=c(2,0,0,0.5), mgp=c(2,1,0), font.main=1)\n",
    "\n",
    "  subtitle <- \"test\"\n",
    "\n",
    "    \n",
    "  za <- rCO3\n",
    "   \n",
    "# x <- var2_e           ;  y <- var1_e\n",
    "if (menu_flag == 15) {\n",
    "    sigcritXa <- sig2$CO3  ;  sigcritYa <- sig1$CO3  #sig2_ is x axis; sig1_ is y axis\n",
    "    plterrcontour (sigcritXa, sigcritYa, xlabel, ylabel, subtitle, xlim, ylim,\n",
    "#                    sig1hp_AtCt$CO3, sig2hp_AtCt$CO3,\n",
    "                   NULL, NULL,\n",
    "                   zenon(sigm2$CO3), zenon(sigm1$CO3),    \n",
    "                   var2_e_soa2, var1_e_soa2,\n",
    "                   x, y, za, levels1,\n",
    "                   'flattest')\n",
    "} else if (menu_flag ==8 || menu_flag == 25) {\n",
    "    sigcritXa <- sig1$CO3  ;  sigcritYa <- sig2$CO3  #sig2_ is x axis; sig1_ is y axis\n",
    "    plterrcontour (sigcritXa, sigcritYa, xlabel, ylabel, subtitle, xlim, ylim,\n",
    "#                    sig1hp_AtCt$CO3, sig2hp_AtCt$CO3,\n",
    "                   NULL, NULL,\n",
    "                   zenon(sigm1$CO3), zenon(sigm2$CO3),\n",
    "                   var1_e_soa2, var2_e_soa2,\n",
    "                   x, y, za, levels1,\n",
    "                   'flattest')\n",
    "} else if (menu_flag == 21) {\n",
    "     za <- t(za)\n",
    "#     sigcritXa <- sig1$CO3  ;  sigcritYa <- sig2$CO3  #sig2_ is x axis; sig1_ is y axis\n",
    "    sigcritXa <- sig2$CO3  ;  sigcritYa <- sig1$CO3  #sig2_ is x axis; sig1_ is y axis\n",
    "    plterrcontour (sigcritXa, sigcritYa, xlabel, ylabel, subtitle, xlim, ylim,\n",
    "#                    sig1hp_AtCt$CO3, sig2hp_AtCt$CO3,\n",
    "                   NULL, NULL,\n",
    "                   zenon(sigm2$CO3), zenon(sigm1$CO3),\n",
    "                   var2_e_soa2, var1_e_soa2,\n",
    "                   x, y, za, levels1,\n",
    "                   'flattest')\n",
    "}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Figure:*** Percent relative uncertainty in computed $[\\text{CO}_3^{2-}]$ from the $A_\\text{T}-C_\\text{T}$ input pair as a function of uncertainties in each member of that input pair. Each of those members is indicated by the x- and y-axis labels. Propagated relative error (in percent) is shown for the Southern Ocean (blue solid lines). Sublots without pH as an input member also include the Constants-Pair curve (thick black line), where the total propagated error from the constants is equal to that from the input pair. The balanced-pair line (thin black line) indicates where each member of the input pair contributes equally to total propagated error. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 804,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1] 21\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "length(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 616,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1] 21\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "length(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 414,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1] 121\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "length(za)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 490,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " [1]  0  2  4  6  8 10 12 14 16 18 20\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "var1_e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 491,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " var1_e1  var1_e2  var1_e3  var1_e4  var1_e5  var1_e6  var1_e7  var1_e8 \n",
       "       0        2        4        6        8       10       12       14 \n",
       " var1_e9 var1_e10 var1_e11 \n",
       "      16       18       20 \n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "unlist(list(\"var1_e\" = var1_e))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 494,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1] 21\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "menu_flag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 497,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "menu_flag   var1_e1   var1_e2   var1_e3   var1_e4   var1_e5   var1_e6   var1_e7 \n",
       "       21         0         2         4         6         8        10        12 \n",
       "  var1_e8   var1_e9  var1_e10  var1_e11 \n",
       "       14        16        18        20 \n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "unlist(list(\"menu_flag\"=menu_flag, \"var1_e\"=var1_e))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 498,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "$menu_flag\n",
       "[1] 21\n",
       "\n",
       "$var1_e\n",
       " [1]  0  2  4  6  8 10 12 14 16 18 20\n",
       "\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "list(\"menu_flag\"=menu_flag, \"var1_e\"=var1_e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 499,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       " [1]  0  2  4  6  8 10 12 14 16 18 20\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "list(\"menu_flag\"=menu_flag, \"var1_e\"=var1_e)[[\"var1_e\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 500,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1] 21\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "list(\"menu_flag\"=menu_flag, \"var1_e\"=var1_e)[[\"menu_flag\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}