{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example use of PyDSTool\n",
    "Install: pip3 install PyDSTool"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import PyDSTool as dst\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# name\n",
    "DSargs = dst.args(name='AMOC')\n",
    "# parameters\n",
    "DSargs.pars = { 'mus': 6.25,\n",
    "                'F': 1.1  }\n",
    "# rhs of the differential equation\n",
    "DSargs.varspecs = {'y': 'F - y*(1 + mus*(1-y)**2)'}\n",
    "# initial conditions\n",
    "DSargs.ics      = {'y': 1.5}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEWCAYAAAB1xKBvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGL5JREFUeJzt3XuwHnWd5/H3h5MbBgzBBIyBACqrqKOgp9ApnBFXxei6ouNlw3jBKa1YrozOWDW7zKVkBmu2HJ0ZLyuKmZVVZxW8OxlXRcbL4GXRnCAiF4GICMegCYRbuCf57h9Pp3w4POekc+k85yTvV9VT5+lf/7r72+nkfNL966efVBWSJO3IAcMuQJI0MxgYkqRWDAxJUisGhiSpFQNDktSKgSFJasXAkCS1YmBILSX5TpLbkszta/t4kkry0gl939+0v6Gv7Ygkn0pya5K7k/woyUsmLJckb0tyRdNnPMnnkvxO5zso7YCBIbWQ5Gjg94ACXjph9rXA6X19ZwGvAn7e13Yo8D3gAeDJwCLgfcCnk7yyb10fAN4OvA04FPgPwJeB/7Qn90faFbOGXYA0Q7weuAT4Ib1w+FzfvH8FXptkYVXdBiwHLgcO7uvzp8Bm4I1Vta1pOz/JMuAfknwBeDzwVuB3q+pHfct+qosdknaWZxhSO6+n94v7U8ALkxzeN+8+YDWwoq/vJycs/wLgC31hsd1ngWX0ziSeB4xPCAtp2jAwpB1I8mzgKOCzVbWW3qWmP5zQ7ZPA65MsAJ5D7zJSv0XAzQNWf3Pf/EdN0keaFgwMacdOB75RVbc005+mb8wCoKq+BywG/gr4SlXdO2EdtwBLBqx7Sd/8WyfpI00LjmFIU0hyIPBqYCTJr5vmucAhSZ42ofv/Ad4JPHfAqv4NeEWSv5lwWerVwE30Bs63AeckGa2qsT25H9Ke4BmGNLWXAVuBJwHHN6/jgO/SG6vo90F6YxUXD1jP+4BHAh9L8ugk85KcBvwl8GfVcx3wYXqD4ScnmdP0W5HkzE72TtoJBoY0tdOB/11VN1bVr7e/gA8Br6HvLL2qNlXVN2vAl8xU1a3As4F5wFX0Lj+9A3hdVX2mr+vbmnWfA9xOb7zk5fTuxJKGKn6BkiSpDc8wJEmtdBYYSY5M8u0kVye5MsnbB/RJkg8mWZfk8iRP75t3epLrmtfpE5eVJO1dnV2SSrIEWFJVlyY5GFgLvKyqrurr82Lgj4EXA88EPlBVz2weozAGjNJ7FMNa4BnNp2glSUPQ2RlGVd1cVZc27+8CrgaWTuh2KvDJ5g6RS+jdqrgEeCFwUTOIeBtwEb3HLUiShmSvfA6jeXDbCfSew9NvKb170Lcbb9omax+07pXASoD58+c/44lPfOIeqVmS9gdr1669paoWt+nbeWAkOQj4AvAnVXXnxNkDFqkp2h/eWLUKWAUwOjpaY2N+3kmS2kryy7Z9O71LKslsemHxqar64oAu48CRfdNHAOunaJckDUmXd0kF+BhwdVX94yTdVtN7YFuSPAu4o6puBi4ETkmyMMlC4JSmTZI0JF1ekjoJeB3w0ySXNW1/Qe9RzlTVucBX6d0htQ64B/ijZt6mJO8C1jTLnV1VmzqsVZK0A50FRvP0zkFjEf19it4Xxgyadx5wXgelSZJ2gZ/0liS1YmBIkloxMCRJrRgYkqRWDAxJUisGhiSpFQNDktSKgSFJasXAkCS1YmBIkloxMCRJrRgYkqRWDAxJUisGhiSpFQNDktSKgSFJasXAkCS1YmBIklrp7Ctak5wHvATYUFVPGTD/z4DX9NVxHLC4+T7vG4C7gK3Alqoa7apOSVI7XZ5hfBxYPtnMqnpvVR1fVccDfw78e1Vt6uvy3Ga+YSFJ00BngVFVFwObdtix5zTg/K5qkSTtvqGPYSR5BL0zkS/0NRfwjSRrk6wcTmWSpH6djWHshP8MfH/C5aiTqmp9ksOAi5L8rDljeZgmUFYCLFu2rPtqJWk/NfQzDGAFEy5HVdX65ucG4EvAiZMtXFWrqmq0qkYXL17caaGStD8bamAkWQA8B/iXvrb5SQ7e/h44BbhiOBVKkrbr8rba84GTgUVJxoGzgNkAVXVu0+3lwDeq6u6+RQ8HvpRke32frqqvd1WnJKmdzgKjqk5r0efj9G6/7W+7HnhaN1VJknbVdBjDkCTNAAaGJKkVA0OS1IqBIUlqxcCQJLViYEiSWjEwJEmtGBiSpFYMDElSKwaGJKkVA0OS1IqBIUlqxcCQJLViYEiSWjEwJEmtGBiSpFYMDElSKwaGJKkVA0OS1EpngZHkvCQbklwxyfyTk9yR5LLm9c6+ecuTXJNkXZIzu6pRktRel2cYHweW76DPd6vq+OZ1NkCSEeAc4EXAk4DTkjypwzp5cOs27ntwa5ebkKQZr7PAqKqLgU27sOiJwLqqur6qHgAuAE7do8VN8OR3XsgHvnldl5uQpBlv2GMYv5vkJ0m+luTJTdtS4Ka+PuNN20BJViYZSzK2cePGXSrioHmz2Hzfll1aVpL2F8MMjEuBo6rqacD/BL7ctGdA35psJVW1qqpGq2p08eLFu1TIQXNnsfl+A0OSpjK0wKiqO6tqc/P+q8DsJIvonVEc2df1CGB9l7XMnzuLuzzDkKQpDS0wkjw6SZr3Jza13AqsAY5NckySOcAKYHWXtRw8dxab73+wy01I0ow3q6sVJzkfOBlYlGQcOAuYDVBV5wKvBN6SZAtwL7CiqgrYkuQM4EJgBDivqq7sqk7ojWFsvOv+LjchSTNeZ4FRVaftYP6HgA9NMu+rwFe7qGuQg+bO4he33L23NidJM9Kw75KaFhzDkKQdMzCAg+c5hiFJO2Jg0Lskdd+D29iydduwS5GkacvAoBcYAHff7+NBJGkyBga/DYy7vCwlSZMyMOjdVgv4aW9JmoKBwW/PMHyelCRNzsDAMwxJasPAoO8Mw8CQpEkZGHhJSpLaMDDwkpQktWFgAPPnGBiStCMGBjByQJg/Z4Q77zUwJGkyBkbj8EfO4zd33jfsMiRp2jIwGksOmcf6O+4ddhmSNG0ZGI0lCw5k/e0GhiRNxsBoPGbBPDbcdT8P+sRaSRrIwGgsOeRAqnAcQ5Im0VlgJDkvyYYkV0wy/zVJLm9eP0jytL55NyT5aZLLkox1VWO/xxxyIAA332FgSNIgXZ5hfBxYPsX8XwDPqaqnAu8CVk2Y/9yqOr6qRjuq7yEes2AegOMYkjSJWV2tuKouTnL0FPN/0Dd5CXBEV7W0scQzDEma0nQZw3gj8LW+6QK+kWRtkpVTLZhkZZKxJGMbN27c5QIOmjuLBQfO5sZN9+zyOiRpX9bZGUZbSZ5LLzCe3dd8UlWtT3IYcFGSn1XVxYOWr6pVNJezRkdHa3dqecKjD+bqm+/cnVVI0j5rqGcYSZ4K/C/g1Kq6dXt7Va1vfm4AvgScuDfqefJjHsnPbr6Lrdt2K3ckaZ80tMBIsgz4IvC6qrq2r31+koO3vwdOAQbeabWnPeUxC7j3wa384pbNe2NzkjSjdHZJKsn5wMnAoiTjwFnAbICqOhd4J/Ao4MNJALY0d0QdDnypaZsFfLqqvt5Vnf2evPSRAFzxqzt5/GEH741NStKM0eVdUqftYP6bgDcNaL8eeNrDl+je4xcfxNxZB3D5+B287ISlwyhBkqat6XKX1LQwa+QARo9eyHev2/W7rSRpX2VgTPDcJxzGdRs2c5O310rSQxgYEzzvuMMB+NbPNgy5EkmaXgyMCY5ZNJ9jDzuIL/74V8MuRZKmFQNjgNc+6yh+ctPtXHbT7cMuRZKmDQNjgD94+lLmzxnh3O/8fNilSNK0YWAMcPC82az8/cfx9St/zQ9+fsuwy5GkacHAmMSbn/NYjlh4IP/t85dz+z0PDLscSRo6A2MS82aP8MHTTuA3d97Hm/95Lfc8sGXYJUnSUBkYU3j6soX8/auexpobNrFi1SXceKufzZC0/zIwduDU45ey6nWj/OKWu3n++/6dv/2/V7Hpbi9RSdr/pGrfeZT36OhojY118xXg62+/l/dddC2fv3SckYTfO3YRzzvucEaPXsixhx3MyAHpZLuS1KUka9t+FbaBsZPWbbiLz60d5ys/uZlfNd//PWfkAI5YeCBHHvoIliyYx4IDZ7PgEbNZcOBsDpo7izkjBzBnVu81u3k/+4ADSOi9yMPfA8mE9/y2jwaLfzTaD40cEB7TfM30zjIw9oKq4pe33sOlN97Gtb/ZzI2b7uaXt97Dhrvu5457H+SBLdv2Sh2StOiguYz91fN3admdCYyhf0XrTJWEoxfN5+hF8wfOv+/Brdx+z4Nsvn8LD27dxgNbtvHA1m08uGUb92/dxpatRVXR+3K/oqr3Rea9n/3T9dC2fSff9zj/aKa2L/3nUA81b/bIXtmOgdGRebNHePSCvXMQJWlv8C4pSVIrBoYkqZVOAyPJeUk2JLlikvlJ8sEk65JcnuTpffNOT3Jd8zq9yzolSTvW9RnGx4HlU8x/EXBs81oJfAQgyaHAWcAzgROBs5Is7LRSSdKUdhgYSc7Y1V/WVXUxsGmKLqcCn6yeS4BDkiwBXghcVFWbquo24CKmDh5JUsfanGE8GliT5LNJlid79KNRS4Gb+qbHm7bJ2h8mycokY0nGNm7cuAdLkyT122FgVNVf0btk9DHgDcB1Sf5Hksftge0PCp+aon1QfauqarSqRhcvXrwHSpIkDdJqDKN6n/j5dfPaAiwEPp/kPbu5/XHgyL7pI4D1U7RLkoakzRjG25KsBd4DfB/4nap6C/AM4BW7uf3VwOubu6WeBdxRVTcDFwKnJFnYjJ+c0rRJkoakzSe9FwF/UFW/7G+sqm1JXjLVgknOB04GFiUZp3fn0+xm+XOBrwIvBtYB9wB/1MzblORdwJpmVWdX1VSD55KkjvnwQUnaj+3Mwwf9pLckqRUDQ5LUioEhSWrFwJAktWJgSJJaMTAkSa0YGJKkVgwMSVIrBoYkqRUDQ5LUioEhSWrFwJAktWJgSJJaMTAkSa0YGJKkVgwMSVIrBoYkqRUDQ5LUSqeBkWR5kmuSrEty5oD570tyWfO6NsntffO29s1b3WWdkqQdm9XVipOMAOcALwDGgTVJVlfVVdv7VNWf9vX/Y+CEvlXcW1XHd1WfJGnndHmGcSKwrqqur6oHgAuAU6fofxpwfof1SJJ2Q5eBsRS4qW96vGl7mCRHAccA3+prnpdkLMklSV422UaSrGz6jW3cuHFP1C1JGqDLwMiAtpqk7wrg81W1ta9tWVWNAn8IvD/J4wYtWFWrqmq0qkYXL168exVLkibVZWCMA0f2TR8BrJ+k7womXI6qqvXNz+uB7/DQ8Q1J0l7WZWCsAY5NckySOfRC4WF3OyV5ArAQ+H99bQuTzG3eLwJOAq6auKwkae/p7C6pqtqS5AzgQmAEOK+qrkxyNjBWVdvD4zTggqrqv1x1HPDRJNvohdq7+++ukiTtfXno7+mZbXR0tMbGxoZdhiTNGEnWNuPFO+QnvSVJrRgYkqRWDAxJUisGhiSpFQNDktSKgSFJasXAkCS1YmBIkloxMCRJrRgYkqRWDAxJUisGhiSpFQNDktSKgSFJasXAkCS1YmBIkloxMCRJrRgYkqRWOg2MJMuTXJNkXZIzB8x/Q5KNSS5rXm/qm3d6kuua1+ld1ilJ2rFZXa04yQhwDvACYBxYk2R1VV01oetnquqMCcseCpwFjAIFrG2Wva2reiVJU+vyDONEYF1VXV9VDwAXAKe2XPaFwEVVtakJiYuA5R3VKUlqocvAWArc1Dc93rRN9Ioklyf5fJIjd3JZkqxMMpZkbOPGjXuibknSAF0GRga01YTpfwWOrqqnAv8GfGInlu01Vq2qqtGqGl28ePEuFytJmlqXgTEOHNk3fQSwvr9DVd1aVfc3k/8EPKPtspKkvavLwFgDHJvkmCRzgBXA6v4OSZb0Tb4UuLp5fyFwSpKFSRYCpzRtkqQh6ewuqarakuQMer/oR4DzqurKJGcDY1W1GnhbkpcCW4BNwBuaZTcleRe90AE4u6o2dVWrJGnHUjVwaGBGGh0drbGxsWGXIUkzRpK1VTXapq+f9JYktWJgSJJaMTAkSa0YGJKkVgwMSVIrBoYkqRUDQ5LUioEhSWrFwJAktWJgSJJaMTAkSa0YGJKkVgwMSVIrBoYkqRUDQ5LUioEhSWrFwJAktWJgSJJa6TQwkixPck2SdUnOHDD/HUmuSnJ5km8mOapv3tYklzWv1V3WKUnasVldrTjJCHAO8AJgHFiTZHVVXdXX7cfAaFXdk+QtwHuA/9LMu7eqju+qPknSzunyDONEYF1VXV9VDwAXAKf2d6iqb1fVPc3kJcARHdYjSdoNXQbGUuCmvunxpm0ybwS+1jc9L8lYkkuSvKyLAiVJ7XV2SQrIgLYa2DF5LTAKPKeveVlVrU/yWOBbSX5aVT8fsOxKYCXAsmXLdr9qSdJAXZ5hjANH9k0fAayf2CnJ84G/BF5aVfdvb6+q9c3P64HvACcM2khVraqq0aoaXbx48Z6rXpL0EF0Gxhrg2CTHJJkDrAAecrdTkhOAj9ILiw197QuTzG3eLwJOAvoHyyVJe1lnl6SqakuSM4ALgRHgvKq6MsnZwFhVrQbeCxwEfC4JwI1V9VLgOOCjSbbRC7V3T7i7SpK0l6Vq4LDCjDQ6OlpjY2PDLkOSZowka6tqtE1fP+ktSWrFwJAktWJgSJJaMTAkSa0YGJKkVgwMSVIrBoYkqRUDQ5LUioEhSWrFwJAktWJgSJJaMTAkSa0YGJKkVgwMSVIrBoYkqRUDQ5LUioEhSWrFwJAktWJgSJJa6TQwkixPck2SdUnOHDB/bpLPNPN/mOTovnl/3rRfk+SFXdYpSdqxzgIjyQhwDvAi4EnAaUmeNKHbG4HbqurxwPuAv2uWfRKwAngysBz4cLM+SdKQdHmGcSKwrqqur6oHgAuAUyf0ORX4RPP+88DzkqRpv6Cq7q+qXwDrmvVJkoZkVofrXgrc1Dc9Djxzsj5VtSXJHcCjmvZLJiy7dNBGkqwEVjaTm5Ncs4v1LgJu2cVlp5t9ZV/2lf0A92U62lf2A3ZvX45q27HLwMiAtmrZp82yvcaqVcCqnSvt4ZKMVdXo7q5nOthX9mVf2Q9wX6ajfWU/YO/tS5eXpMaBI/umjwDWT9YnySxgAbCp5bKSpL2oy8BYAxyb5Jgkc+gNYq+e0Gc1cHrz/pXAt6qqmvYVzV1UxwDHAj/qsFZJ0g50dkmqGZM4A7gQGAHOq6ork5wNjFXVauBjwD8nWUfvzGJFs+yVST4LXAVsAd5aVVu7qrWx25e1ppF9ZV/2lf0A92U62lf2A/bSvqT3H3pJkqbmJ70lSa0YGJKkVvb7wNjR40tmkiQ3JPlpksuSjA27np2R5LwkG5Jc0dd2aJKLklzX/Fw4zBrbmmRf/jrJr5pjc1mSFw+zxjaSHJnk20muTnJlkrc37TPuuEyxLzPxuMxL8qMkP2n25W+a9mOaRyxd1zxyac4e3/b+PIbRPG7kWuAF9G7lXQOcVlVXDbWwXZTkBmC0qmbch5GS/D6wGfhkVT2laXsPsKmq3t2E+cKq+u/DrLONSfblr4HNVfX3w6xtZyRZAiypqkuTHAysBV4GvIEZdlym2JdXM/OOS4D5VbU5yWzge8DbgXcAX6yqC5KcC/ykqj6yJ7e9v59htHl8ifaCqrqY3p1y/fofHfMJev/Ap71J9mXGqaqbq+rS5v1dwNX0nrgw447LFPsy41TP5mZydvMq4D/Se8QSdHRc9vfAGPT4khn5l6hRwDeSrG0emTLTHV5VN0PvHzxw2JDr2V1nJLm8uWQ17S/j9GueJH0C8ENm+HGZsC8wA49LkpEklwEbgIuAnwO3V9WWpksnv8v298Bo/QiSGeKkqno6vScEv7W5NKLp4SPA44DjgZuBfxhuOe0lOQj4AvAnVXXnsOvZHQP2ZUYel6raWlXH03sKxonAcYO67ent7u+BsU89gqSq1jc/NwBfYuY/4fc3zbXn7degNwy5nl1WVb9p/pFvA/6JGXJsmmvkXwA+VVVfbJpn5HEZtC8z9bhsV1W3A98BngUc0jxiCTr6Xba/B0abx5fMCEnmN4N5JJkPnAJcMfVS017/o2NOB/5liLXslu2/YBsvZwYcm2Zw9WPA1VX1j32zZtxxmWxfZuhxWZzkkOb9gcDz6Y3JfJveI5ago+OyX98lBdDcRvd+fvv4kr8dckm7JMlj6Z1VQO+RL5+eSfuS5HzgZHqPaf4NcBbwZeCzwDLgRuBVVTXtB5Mn2ZeT6V32KOAG4M3bxwGmqyTPBr4L/BTY1jT/Bb1r/zPquEyxL6cx847LU+kNao/Q+0//Z6vq7OZ3wAXAocCPgddW1f17dNv7e2BIktrZ3y9JSZJaMjAkSa0YGJKkVgwMSVIrBoYkqRUDQ+pQkkOS/Ndh1yHtCQaG1K1DAAND+wQDQ+rWu4HHNd+18N5hFyPtDj+4J3WoeTLqV7Z/L4Y0k3mGIUlqxcCQJLViYEjdugs4eNhFSHuCgSF1qKpuBb6f5AoHvTXTOegtSWrFMwxJUisGhiSpFQNDktSKgSFJasXAkCS1YmBIkloxMCRJrfx/nzf1o4bOtIoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "DSargs.tdomain = [0,30]                         # set the range of integration.\n",
    "ode  = dst.Generator.Vode_ODEsystem(DSargs)     # an instance of the 'Generator' class.\n",
    "traj = ode.compute('polarization')              # integrate ODE\n",
    "pts  = traj.sample(dt=0.1)                      # Data for plotting\n",
    "\n",
    "# PyPlot commands\n",
    "plt.plot(pts['t'], pts['y'])\n",
    "plt.xlabel('t')                                 \n",
    "plt.ylabel('y')                                \n",
    "plt.ylim([0,2])                                 \n",
    "plt.title(ode.name)                             \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VPW9//HXZ2ay73uAhD1hFUQCqAgCbmhd61Kxtupta21rrbW3v9r23va2vba9XW7vra1XrXWpu+JSXFoFRVFZJCggWwJhS0gg+77PfH5/zDgETICETCZhPs/HA5M58z3nfA4jeed8v99zjqgqxhhjDIAj2AUYY4wZPCwUjDHG+FkoGGOM8bNQMMYY42ehYIwxxs9CwRhjjJ+FgjFBJiI3i8j7x3h/nogUDGRNJnRZKJhBTUTeEZEaEYk4avmjIqIicvlRy//Ht/zmLsuyRORJEakSkSYR+VBELj1qPRGRO0Rki69NiYg8LyKnBfQAu+Grf/ynr1X1PVWdcKLtRSTXV3uliNSJyGYRuUtEnIGu3Qx9Fgpm0BKR0cA8QIHLu2lSCNzUpb0LuBYo6rIsGXgfaAemAKnAH4CnROSaLtv6X+A7wB1AMpALvAx8rr+OZyCIyDhgHVAMnKaqCXj/TvKAuGDWZoYGCwUzmH0ZWAs8Spcf/l28AswVkSTf68XAZuBglzbfBRqBr6jqQVVtUdWngXuA3/vOEHKAbwFLVPVtVW1T1WZVfVJVf91dYb4zmP8UkdUi0igir4hIiu+MpF5E1vtCDREZ7ftt3nXU+l/tZrurfN9u8m33CyKyQERKTvDv7GfAalW9S1XLAFS1QFVvUNVaEYkUkSd8Z021vjozTnDbJgRYKJjB7MvAk74/F3Xzw6sVWAZc36X9345qcwHwgqp6jlr+HDAS7xnBeUCJqn7Yy/quB74EjADGAWuAR/CeaWwHftrL7aGq833fTlfVWFV9tpebOB9Yeoz3bwISgGwgBbgNaOltnebUZaFgBiUROQcYBTynqhvwdgnd0E3TvwFfFpEE4Fy8XT5dpQJl3axX1uX9lB7aHM8jqlqkqnXAP4AiVV2hqp3A88CMPmzzZB3vWDp8bcarqltVN6hq/cCUZoYCCwUzWN0EvKmqlb7XT9FNF5Kqvg+kAf8GvKqqR//WWwkM62b7w7q8X9VDm+M51OX7lm5ex/ZhmyfreMfyOPAG8IyIlIrIb0QkbGBKM0OBhYIZdEQkCrgOOFdEDorIQbxjA9NFZHo3qzwBfI/Pdh0BrACuFpGj/1+/Du9gbCHwFpAlInn9dQxHafJ9je6yLDNA+1oBXN3Tm6raoao/U9XJwNnApXi73YwBLBTM4HQl4AYmA6f7/kwC3qP7H2B/xDt2sKqb9/4AxAN/FZFM30DrEuDHwPfVaydwH/C0b1A33NfuehG5+2QPRlUrgAPAjSLiFJF/wTsG0ZNDwNg+7u6nwNki8lsRyQQQkfG+weVEEVkoIqf5pqfW4+1OcvdxX+YUZKFgBqOb8PbX7/fNGDqoqgeBPwFf7DqLB0BVq1X1Le3m4SCqWgWcA0QC2/B2r9wFfOmoQdw7fNv/M1CLdwzjKrwznPrD14Dv+/Y/BVh9jLb/ATzmmx10XW92oqpFwFnAaGCriNQBLwD5QAPeM5SleANhO/Au3jMtYwAQe8iOMcaYT9mZgjHGGD8LBWOMMX4WCsYYY/wsFIwxxvi5jt9kcElNTdXRo0cHuwxjjBlSNmzYUKmqacdrN+RCYfTo0eTn5we7DGOMGVJEZN+JtLPuI2OMMX4WCsYYY/wsFIwxxvhZKBhjjPGzUDDGGOMXsFAQkYdFpFxEtvTw/hd9DxTf7HukYXe3RDbGGDOAAnmm8CjeZ+b2ZA9wrqpOA34BPBjAWowxxpyAgIWCqq4Cqo/x/mpVrfG9XAtkBaoWgO2NLfxnUSkNnXbreGOM6clgGVP4Ct5n3HZLRG4VkXwRya+oqOjTDva3tvOn/eUUNrX2tUZjjDnlBT0URGQh3lD4QU9tVPVBVc1T1by0tONepd2tCTGRABQ0WygYY0xPgnqbCxGZBjwEXOx7QlbAZEeGE+kQO1MwxphjCNqZgoiMBF7E+1jEwkDvzynC+OhICiwUjDGmRwE7UxCRp4EFQKqIlOB9oHgYgKreD/wESAHuExGATlXNC1Q94O1CWlvbGMhdGGPMkBawUFDVJcd5/6vAVwO1/+5MionkhUM1VLV3khI+5G4Qa4wxARf0geaBNDMhBoAN9U1BrsQYYwankAqF6XHRuAQ21DcHuxRjjBmUQioUop0OJsdGsc7GFYwxplshEwqqSlNTEfOT4sivb7Irm40xphshEwplB19g7boLmRtTT6fCO9UNwS7JGGMGnZAJhZTkcwEHw5v+SaLLyZtVdcEuyRhjBp2QCYWIiDSSkuZQVb6MS1ITeK2ijkbrQjLGmCOETCgAZGZcSUvLPi6LLaXZ7eHl8tpgl2SMMYNKSIVCRsZlhIWlkFj5AJNjInmguBy3arDLMsaYQSOkQsHpjCA7+yaqa1Zxa3oTO5vbeOlQzfFXNMaYEBFSoQCQnXUzERGZjD70c06LjeTnRaXUdHQGuyxjjBkUQi4UXK4YcnP+naamLdwZvYLqjk7+taAYj3UjGWNM6IUCQHr6YrJGfInIQ3/gW8mlvFZRx892laIWDMaYEBeSoQCQk/Mj0lIvYFblt7kmdh8PlFTw7e37abJpqsaYEBYyoVBRUcFjjz3Gnj17UFUcjnCmTr2XEcOv58qGu7ghbAUvHKpmwfoCXjxUY7OSjDEhKWQeKlBXV+cPhszMTKZOncr48eOZkPsLUpLnE7Xrl4zXt3m07dt8c1s79+zax1WZaSxOTeS0uCgiHCGTn8aYECZDrR89Ly9P8/Pz+7RuR0cHGzduZNOmTZSUlAAQERFBRkYGycnxJCbuwMMa1nkSWcl5bJHpeHAShptxYU2MivAwNiqcYVGRJIdFkBoeQXJ4NDFhUUS5wol0hhHldBHhEBzep8kZY8ygICIbTuTpliEVCl3V1tayb98+iouLKS8vp7q6msZG7y21w8KbSUw8iDOxiQNxGeyLzKLUOYKDDKOKVFSOf9bgUDeC4sCDoF3+eHB0+f7T5d078WDpeRu9cfLbsCg0JnDOqfyEv1z3nT6te6KhEDLdR0dLTEwkMTGR6dOn+5e11jRRvbWM+qIKGg7V0lzaSK56cOOhU9rxhBfSGrGB+qh2GsOEhnAHzWFO2p1Cp1NodzjpFCftDhceETwC8OlXUJEjIsIj3q990s8/ffsnUiwSjAmkzIaWgO8jZEPhU6pK645qmtaW0VpYAwrxMWGkZGXTnFnPnj0bKdnzCc0d9XgihPDYEbQ0xaGeRBLCkxgd7iHm4D5iKvcQ1VZNzIgkKjM6WB++h62xNVTHC8nDxjAmexrjknMYkzCG9NZG0ra8TNL2f+BsbwBnBGTNgmHTIDUHUsZDbCbEpkFkIlhXlDFmgIR0KLTtq6f21d10FDfgiAsnbkE2UVNSqGgq4Z+P3E/5niJiU1LJOXc+DbXZFBc4cONk0tnJpJWsxvXSH3B0tBI7fz4xN1/JspT9PLjnSRraG5iVOYurx17GvKx5pEalendYXwZv/BC2vgThsTD5Sjjtahh5NoRFBvcvwxhjCNFQULeH+jf30bCqBGdcOEnX5BA9Ix1Fee/px8h/9SVik1O46LbvEBY9mfeeLcLd4WHm4mwmDW+g8sf/SkdxMQlXXUXqN26jNMHDN977AdsKtnFu1rncPuN2JiZPPHKnnyyFV+8CdxvM/z6c+U2ITg7OX4AxxvQg5ELB09ZJ1ZM7aCusIWZ2JgmfG4sjwklbczPL/vuX7P9kI9MvuJh5N9zMprcOsf7ZAjLHxnPeTZOR/JWUffVHuNLSGPm3x4iZPZv1B9dz5+t34hAHvz/391w4+sIjd6gKb/8C3vs9ZJ8JV94HKeOCc/DGGHMcIRUKnjY3lQ9tof1AA0lX5xAzKxOA9pZmXvzVTzlYVMhFt32HKQvO54MXdrFpRTETzx7Gghsm0PiP1yj9fz8gauYZZN17L66kJNaVreObK77JiLgR3HfefWTFZR25Q1V440ew9j4448twye/BFR6EIzfGmBMTMqGgbg9VT26nvaSBlBsnETXF28/vcbv5++/uoWxXAZfe+QNy58xlwz/3smlFMdMWZnHOdTk0vvsupXf/kOjZs8l+4H4ckZFsrdrKHW/fwcj4kTxy0SMkRiZ+dqerfucNhDm3weJf24CxMWbQC5nLdJs2HKKtsIakq3L8gQCw6smH2b9lExd+/Q5y58xl35Yq1r68m5xZGZxzbQ4d+/dT+v3/R+SECWTf92cckZHUtdXx3ZXfJSEigQcueKD7QCh8A1b+J0z7ggWCMWbICJkzhZi8TFzJkUSOT/Iv27V+LRte+zszFl/G1AXn01jTyopHtpEyIpZFX5qIdnZQcsd3EIeDEX/8I46YGFSVH7//YypbKnn84sdJj07/7M7qSuDFr0HmNLjsfy0QjDFDRsDOFETkYREpF5EtPbwvIvJHEdklIptF5IxA1QIgDjkiEFoaG1j+lz+RNnos537pK6gq7zxVQGeHm8W3TsUV7qTqgQdpKyhg2K9+RXjWCACWFS3j3ZJ3+V7e95iSOuWzO1L1zTLqgOseg7CoQB6WMcb0q0B2Hz0KLD7G+xcDOb4/twL/F8BaPuPdvz1Ea2MDi79xJ06Xi6KPKtj3SRVzLh9LYkY0rYWFVD74IPGXXUbcooUA1LTW8Lv83zEjfQZLJi7pfsOfLIWdb8B5P4HksQN4RMYYc/ICFgqqugqoPkaTK4C/qddaIFFEhgWqnq7KdhWw9d23yLv0KtJHj6W9tZP3niskbWQc0xZmoaoc+tWvcEZHk/GjH/rX+8OGP9DY0chPzvwJju7uf9TeBG/+GEbMhNm3DsShGGNMvwrmQPMIoLjL6xLfss8QkVtFJF9E8isqKk5qp6rKu4//leiEROZcdR0Am98uprmunfnX5+JwOmh6/32a16wl9VvfxJXk7XLaVbOLvxf9nS9O/CLjk8Z3v/E1f4bGQ96BZYfzpOo0xphgCGYodDf62u192VT1QVXNU9W8tLS0k9rprvy1HNixjbnX3Uh4VDStjR18/OZ+xkxPJXNsAup2U/7b3xGWnU3S9df71/vjx38k2hXNV0/7avcbbqyAD/4XJl0G2bNPqkZjjAmWYIZCCZDd5XUWUBrIHaoqa5c+Q2LmMKYuvACADW/so6PNzZwrvP3/DW+8QVthIenfvRMJ915otqliEyuLV3LL1Fu6n34K8P4foKMFzvtpIA/BGGMCKpihsAz4sm8W0plAnaqWBXKHezduoHxvEbOvvBaH00lrUwdbVh0gZ3YGKcNjUVUqH/wL4WPHErf48Bj5Q5sfIjEikRsn3dj9hpurYcOjcNo13rucGmPMEBWw6xRE5GlgAZAqIiXAT4EwAFW9H3gduATYBTQDtwSqFt8+WfvSc8SlpDF5nnc20ZZ3D9DZ5uaMC0cB0PTee7Tt2MGwX/4S8T1+c3ftbt4peYdvTP8G0WHR3W/8w79ARxPM7dvDL4wxZrAIWCioag9zNv3vK/CtQO3/aCXbt1BasI1Ft3wdpyuMzg43m1cWM3JKMikjYgGoevAvuIYNI+HSz/nXe3Tro0Q6I7l+4vXdb7i9CdbdD7mLIaOb6xaMMWYICZnbXIRHRpF75jlMXeS9i2nB2oO0NHQww3eW0LJlK835+aTccrN/LKG8uZxXd7/KFeOvIDmyh9tcb3wKWqph7p0DchzGGBNIIXObi4yx47nsu3cD3q6kTW8VkzYyjhG53oHj2mefQaKiSLjqKv86zxY8S6enk5sm39T9RlW9XUfDz4BRZwX8GIwxJtBC5kyhq9KdtdQcbOa0BVmICO6GBupefY2ESz+HMy4OgA5PBy/ufJF5WfPIjs/ufkN734fKApjVwzRVY4wZYkIyFLauOkBEtIvxed6b2dW98gra0kLidV/wt3mn+B0qWyq5Lve6nje0/iHvM5Snfj7QJRtjzIAIuVBorm+n6OMKJpyZSVi4E1Wl9plniZwyhajTpvrbPV/wPJkxmZwz4pzuN1RfBjtehRk32k3vjDGnjJALhR1ryvC4lSnzvHfUaNm4kbbCQhKvP3yWsL9+P2vK1nB1ztU4e7pdxUd/A08nzPrKQJRtjDEDIqRCQT3K1vcOMDwnkeRhMQDUPvMsjpgYEi65xN9uaeFSnOLk8zk9dAu5O2HDIzD+fLsTqjHmlBJSoVC8o5r6ylamzB8OgLu2lvp//pOEKy7HEeMNiXZ3Oy/vepmF2Qu7f4AOQOE/oaEM8uwswRhzagmpUNi6qpTI2DDGne4bYP7739G2NhK/cLjraMW+FdS01XBt7rU9b+ijxyBuOORcGOiSjTFmQIVMKDTVtrFncyWTzh6GM8yBqlLz7HNEnX46kRMm+Ns9X/g8WbFZnDn8zO43VFsMO5d7B5idIXOZhzEmRIRMKBzYWQPAlHnerqPm9etp3737iLOE3bW7yT+Uz7UTru3+IToAHz/h/XrGlwJarzHGBEPI/KqbOyuT7InJRMV5b2FR+8yzOOLjib/48N1Ql+5cisvh4opxV3S/EY8bPn4cxp8HiSMHomxjjBlQIXOmAPgDobOqivrly0m48gockZEAtLnbWFa0jPNGnkdKVEr3G9i1AuoPwBk93PbCGGOGuJAKhU/VvfQSdHSQ1KXraPm+5dS11XFN7jU9r7jhUYhJhwkXB75IY4wJgpALBfV4qHnueaLz8ogYN86/fGnhUkbGjWR2Zg+P0qwvg8I3YMYXwRk2QNUaY8zACrlQaPpgNR3795O45PDzEXbX7WbDoQ1cnXt1zwPMG58AdcMZXx6gSo0xZuCFXCjUPPMMzpQU4i+4wL/shcIXjjPA7PHe1mLMfLuC2RhzSgupUOgoK6Nx5UoSr77a/yCdTweYF2Uv6nmAuegtqN1vA8zGmFNeSIVC7fPPgyqJ1x2+Hfabe9+ktq322APM6x6A2AyYdPkAVGmMMcETMqGgHR3UPr+UmPnzCM/y3iFVVXly+5OMSRjDmcN6uIK5qgh2LYe8fwFX+ABWbIwxAy9kQqHh7ZV0VlSQdP3hAeZNFZvYWrWVL078IiLS/YofPgiOMJh5ywBVaowxwRMyoRA1fRppd91F7Pz5/mVPbH+CuPA4Lht3WfcrtTXAx0/ClKsgLmOAKjXGmOAJmVAIy8wk9davIU7vQ3MONh1kxb4VXJ1zNdFh0d2vtPFpaG+AOV8fwEqNMSZ4QiYUjvbU9qdQlCUTl3TfwN0Ba+6FrFmQlTewxRljTJCEZCjUtdXxbMGzXDTqIobHDu++0ZYXvNNQ531vYIszxpggCslQeHL7kzR3NvO1aV/rvoHHA+/9N6RPgZyLBrY4Y4wJopALhcb2Rp7Y/gTnjTyPnKSc7hsVvAaVBTDvLnCE3F+RMSaEBfQnnogsFpECEdklInd38/5IEVkpIh+LyGYRuSSQ9YD3LKGhveEYZwluWPkr7+0sJl8Z6HKMMWZQCVgoiIgT+DNwMTAZWCIik49q9m/Ac6o6A7geuC9Q9QBUtVTxyNZHWJi9kCkpU7pvtPk5KN8Ki/7dHrdpjAk5gTxTmA3sUtXdqtoOPAMcfcc5BeJ93ycApQGsh/s33U9rZyvfnfnd7ht0tMLKe2D4DDtLMMaEpED+KjwCKO7yugSYc1Sb/wDeFJFvAzHA+d1tSERuBW4FGDmyb4/B3Fu3l6WFS7km9xrGJIzpvtHa+6CuGK74k40lGGNCUiB/8nV33wg96vUS4FFVzQIuAR4X+ewDDVT1QVXNU9W8tLS0PhVT3FBMRkwGt02/rfsGNXvh3d/AxEth7II+7cMYY4a6QJ4plADZXV5n8dnuoa8AiwFUdY2IRAKpQHl/FzMvax6vDX8Np8P52TdV4fXvgzjg4v/q710bY8yQEcgzhfVAjoiMEZFwvAPJy45qsx84D0BEJgGRQEWgCuo2EMA7uLzzTVj4I0jICtTujTFm0AtYKKhqJ3A78AawHe8so60i8nMR+fTBBN8DviYim4CngZtV9egupsCq3gOvfQ9GngVzeuhaMsaYEBHQOZeq+jrw+lHLftLl+23A3EDWcEydbfDi17zdRp9/0KagGmNCXuj+FFSFV78LJevh2scgsW+zmowx5lQSuvMu3/8DbHwSzr0bptg1CcYYA6EaCmvvh7d+BlOvgXN/EOxqjDFm0AitUFCFD/4I//wBTLoMrrrfLlIzxpguQmdMwd0J//g+5D8Mk6+Azz8EzrBgV2WMMYNK6PyavPEJbyDMvROueRRc4cGuyBhjBp3QOVOY8SXvhWnju729kjHGGELpTMHhtEAwxpjjCJ1QMMYYc1wWCsYYY/wsFIwxxvhZKBhjjPGzUDDGGONnoWCMMcYvZELB41E8noF9VIMxxgw1IRMK6/ZUc85/vc0vXt3GR/trGOhn+RhjzFAQMlc0R4Y5mDw8nsfX7OOv7+9hVEo0N501mmvzsoiLtHsgGWMMgAy135jz8vI0Pz+/z+vXt3awfOshnvpwPxv21RAX4eK2BeP4yjljiAzr4RnOxhgzxInIBlXNO267UAuFrjYV13Lv2ztZsb2c4QmR3HPVaSycmN4v2zbGmMHkREMhZMYUujM9O5GHbprFM7eeSVxkGLc8up4fvriZlnZ3sEszxpigOG4oiMjtIpI0EMUEy5ljU1j27bl8/dyxPLO+mOsfXEN5fWuwyzLGmAF3ImcKmcB6EXlORBaLiAS6qGCIcDn54cWTeODGmewsb+SKP3/A3sqmYJdljDED6rihoKr/BuQAfwVuBnaKyC9FZFyAawuKC6dk8vxtZ9HW6WHJX9ayr8qCwRgTOk5oTEG9o9EHfX86gSRgqYj8JoC1Bc2U4Qk8+dU5tHa4+dJfP6S6qT3YJRljzIA4kTGFO0RkA/Ab4APgNFX9BjATuDrA9QXNpGHxPHzzLA7Wt/KNJzbQ3ukJdknGGBNwJ3KmkAp8XlUvUtXnVbUDQFU9wKUBra4fVbR3cE9RKZ29uNXFjJFJ/ObqaazbU83v3iwIYHXGGDM4nMiYwk9UdV8P723v/5ICY3VtI/fuL+f+4vJerXfljBHcMGckf3lvN2uKqgJUnTHGDA4BvU7BN1upQER2icjdPbS5TkS2ichWEXkqULVcnpbIJakJ/HbvQXY29W666b99bhKjU2L41+c30dTWGaAKjTEm+AIWCiLiBP4MXAxMBpaIyOSj2uQAPwTmquoU4M4A1sOvc7OIdjj4XkFxr26IFx3u4nfXTuNAbQt/WrkrUCUaY0zQBfJMYTawS1V3q2o78AxwxVFtvgb8WVVrAFS1d307vZQeEca/jRvOh3VNvFJR16t1Z45K5pqZWTz03m52VzQGqEJjjAmuQIbCCKC4y+sS37KucoFcEflARNaKyOLuNiQit4pIvojkV1RUnFRR1w9LZlJMJP9ZVEqru3czin6weCKRYU5++fqOk6rBGGMGq0CGQndXPh/dZ+PCe2HcAmAJ8JCIJH5mJdUHVTVPVfPS0tJOqiinCD8bP4L9re38rbSyV+umxUXw9fljWbH9EJtLak+qDmOMGYwCGQolQHaX11lAaTdt/q6qHaq6ByjAGxIBNT85jrMTY7lvfwVtnt6dLdx09mgSo8P4w/LCAFVnjDHBE8hQWA/kiMgYEQkHrgeWHdXmZWAhgIik4u1O2h3AmvzuHJXBwfYOnjtY3av14iLDuHX+WFYWVLCp2M4WjDGnloCFgqp2ArcDbwDbgedUdauI/FxELvc1ewOoEpFtwErg+6o6IBcDzEuKZUZcNPfuK+/VBW0AXz5rNHERLh7+YE+AqjPGmOAI6HUKqvq6quaq6jhVvce37Cequsz3varqXao6WVVPU9VnAllPS8sB//ciwu2j0tnf2s6KqvpebSc2wsV1s7J5bXMZB+vsFtvGmFNHyDxkp6zsRVavmU9T0+HeqYtSEhgeEcYjB3o34Axw89mj8ajy+Nq9/VilMcYEV8iEQnLyPESclJY951/mcghfGp7CuzUNFDX37jf+7ORozp+UwdMfFtvN8owxp4yQCYWIiDRSUxZRVvYiHs/hW2F/cVgKLoG/Hej9UMaS2SOpbmrn7R2H+rNUY4wJmpAJBYDhw79AR0cVlZUr/cvSI8K4KDWBpYdq6OjlgPO8nFTS4yJ4Pr+kv0s1xpigCKlQSEmZT0REJqVlzx6x/NqMZKo6Onm3pqFX23M5HVw9M4t3Civsmc7GmFNCSIWCiJPMzCuprn6f9vbD1ycsSokjyeVkaS+vWQC4dmYWbo/y4scHjt/YGGMGuZAKBYCM9EtRdVNe8U//snCHg8vTE/lnZR0Nne5ebW9sWiynZyfyyqajL9Y2xpihJ+RCITZ2ItHR4zh06NUjll+bmUyrR3m9l3dPBfjcacPYWlrPvqqm/irTGGOCIuRCQUTIyLiU2toPaWs7PGtoZnw0IyLCeL2y97euuPi0TABe+6Ss3+o0xphgCLlQAG8XEiiHyl/3LxMRLklL4J3qBpp62YWUlRTN9OxEXttsoWCMGdpCJhT279/PE088QVtbGzExY4mJyaWiYvkRbS5JS6TNo7xV3btZSACXWheSMeYUEDKh4PF42LVrF0VFRQCkpZ5PXV0+HR2Hu4tmJ8SQEubi9Yq+dyH9Y8vB/inYGGOCIGRCITs7m8jISAoKCgBITTsfVTeVVe/42zhFWJwaz/Kq+l4/lS0rKZpJw+J5e3tAnyhqjDEBFTKh4HQ6ycnJYefOnXg8HuLjTiM8PJ3KihVHtLskLZEmt4f3enkhG8D5k9LJ31dNbXP78RsbY8wgFDKhAJCbm0tzczMlJSWIOEhNXURV9So8njZ/m7mJsUQ5HH0aV1g0MR2PwjsFJ/ccaWOMCZaQCoXx48fjcDgoLPQ+SjMt9Xzc7iZqatb620Q6HcxPjuWtqnpUe3cvpOlZiaTGhvPWDutCMsYMTSEVClFRUYwcOdIfCklJZ+N0RlNR+fYR7c5Ljqe4tZ2dzW3dbaZHDoewcEI67xaU09HLMQljjBkMQioUwNuFVF5eTk1NDU4FJgJmAAAXmUlEQVRnBElJZ1FdteqIs4JFKfEAvNXLJ7IBnDcpg/rWTvL31vRbzcYYM1BCLhQmTJgA4D9bSEmeT0vrflpa9vrbZEWGMzEmsk+hMC8nlXCnw56xYIwZkkIuFFJSUkhOTmbnzp2+1/MBqKpadUS781LiWVfX1Osb5MVEuJgzNpmVNthsjBmCQiYUPG2dtGyrQj1Kbm4ue/fupb29naiokURFjaaq+qhQSI6nQ7VPU1MXTkhnV3kjxdXN/VW+McYMiJAJhZatVVT9bRsdpY3k5OTQ2dnJnj17AO/ZQk3NWtzuwwPLsxJiiHM6WNGHLqSFE9MBWFlgs5CMMUNLyIRCZG4SCLQW1DBq1CjCwsIOdyElz8fjaaW2br2/fZhDODc5jpXVDb2emjomNYbRKdG8bVNTjTFDTMiEgjM2nLARsbQWVONyuRg3bhw7d+5EVUlKmoPDEU71UeMKi1LiKWvrYHtT7x+1uXBiOmuKqmhp792YhDHGBFPIhAJA5IRk2osbcDd1kJOTQ11dHeXl5Tid0SQmzP7MuMKi5L5PTV04IZ22Tg9rd1f1S+3GGDMQQiwUkkChbWcNOTk5AF1mIZ1LU9NOWlsPP1YzMyKMKbF9m5o6e0wyUWFO60IyxgwpIRUK4VlxOKJdtBbUEB8fT2ZmZjdTU989Yp3zkuNZX99EfS+npkaGOZk7PpWVBeW9HpMwxphgCWgoiMhiESkQkV0icvcx2l0jIioieQGtxyFE5ibRWliDepScnBz2799PS0sL0dHjiIwY/plQWJQSj1thVR9ukLdwYholNS3sKm/sr0MwxpiAClgoiIgT+DNwMTAZWCIik7tpFwfcAawLVC1dRU5IxtPUQccB79RUVaWoqAgRISXlXKprVuPxHL71dV58DPEuB29X921cAWxqqjFm6AjkmcJsYJeq7lbVduAZ4Ipu2v0C+A3Q+yk+fRDhn5paTVZWFlFRUYe7kFIX4nY3UVt7eGqqyyGcmxTP21W9n5o6PDGKiZlxNq5gjBkyAhkKI4DiLq9LfMv8RGQGkK2qrx5rQyJyq4jki0h+RUXfbx+hqjhjwgjPiqO1oAaHw8H48eP9D95JTjoLkfBuupDiONjewbY+TE1dMCGd/L011Ld29LluY4wZKIEMBelmmf9XbRFxAH8Avne8Danqg6qap6p5aWlpfSqmaMM6HrjtyzTV1hA5IYn2kgbcje3k5OTQ3NxMaWkpTmc0SUlzjnhEJxyemvp2H2YhLZqYTqdHeX9nZZ/qNsaYgRTIUCgBsru8zgJKu7yOA6YC74jIXuBMYFmgBpvjUtJoqq1h98friZyQ7JuaWsv48eMREX8XUmrKApqbi2hp2e9fNyMijKmxUX2amnrGyETiI12stC4kY8wQEMhQWA/kiMgYEQkHrgeWffqmqtapaqqqjlbV0cBa4HJVzQ9EMWmjxhCXksbuDR8SNiIWR0wYLQXVREdHk5WVdfhW2ikLAD5ztnBeindqal1HZ6/263I6mJ+bxjuFFXg8NjXVGDO4BSwUVLUTuB14A9gOPKeqW0Xk5yJyeaD22xMRYezM2ezd/DHuzg4ic5No6zI1taysjIaGBqKjR3vvmvqZLqQ479TUmt5PL104IZ2Khja2lvb+TMMYYwZSQK9TUNXXVTVXVcep6j2+ZT9R1WXdtF0QqLOET42fOZvOtjb2b91E5IQkPM2dtJc0kJubC8CuXbsASE1d6Ltraot/3ZnxMSS4nH2amnruhDREbGqqMWbwC6krmrOmTCMsMordGz4kIufwXVMzMjKIi4s7ogvJ42mjpmatf12XQ5ifFMfbVfW9npqaGhvBtKxEm5pqjBn0QioUXGFhjJ42g6INH+KIdhGeHUfrjmpEhJycHIqKiujs7CQpcRZOZywVFW8esf75KfEcau9kY0NLD3vo2cIJaWwqqaWqse34jY0xJkhCKhQAxuXNobG6ivI9RURNSaHjQCOd1a3k5ubS3t5OcXExDkcEaannUV7xJh7P4esLLkyNxyXwSnltr/e7aGI6qvBuoT2m0xgzeIVcKIw9YxYOp5OCNe8RdZr3moeWTyoYM2YMTqfT34WUnn4xnZ21R3QhJYW5mJ8UxysVtb3uQpo6PIHU2HB7drMxZlALuVCIiotn1LQZ7Fi9CmdiOGHZcTRvriQiIoJRo0ZRWFiIqpKcPB+nM5by8tePWP+y9ESKW9t73YXkcAiLJqazckc5rR324B1jzOAUMqFQX9nCuld243F7mDj3XBoqKygt3EH0aaneLqTKFiZPnkxVVRVlZWU4nRGkpi6ionL5EV1IF6cmECbSpy6kK04fQWNbJ8u3HerPQzPGmH4TMqFQWdJI/mt72b+1mvF5c3CFhbNj9btETUsFoHlzBVOmTMHpdLJp0yYAMtIvoaOjhurq9/zbSfR1If29vAZ3L7uQzhybQmZ8JC99fKD/DswYY/pRyITCqNNSiIoPZ9sHpYRHRTN25mwKVr+HxDoJH5NA04ZDREZEMmHCBD755BPcbjcpKQsIC0uhtPS5I7b1hWHJHGjr4N1ePmPB6RCumDGcdwsrbBaSMWZQCplQcDodTDwzk72fVNFU18bUhRfQ0lDPznWriZ2TibuqlbbdtUyfPp3m5mZ27tyJwxHGsGGfp7JqJW1thweIF6fGkxzm5Mmy3j9/+fMzsnB7lJc3lh6/sTHGDLCQCQWAyXOHox5lx5oyRk+bQUJGJpuW/4OoKak4ol00fXiQ8ePHExcXx4cffgjA8GHXodpJ2cEX/dsJdzi4LjOZNyrrqGjv3S2xJ2TGMWNkIo+v2Wv3QjLGDDohFQqJGdGMyE1ky7sH8Hhg2nmLKdm+haqDxUSfkUHL1ipo6mT27Nns3r2b8vJyYmLGkpg4h5KSx48YcL5xeAqdCo8d6P3Zwi1zx7C3qtlue2GMGXRCKhQAZlw4isaaNgo/PMTUhRfgiohg/d+XEnvWMFClYdUBZs6cicvlYu1a7zUKo0Z+jba2Mg4desW/nfHRkVyYEs/DBypocvduiunFUzPJjI/k4Q/29OuxGWPMyQq5UBg5JZmUEbF8/OY+omLjmX7BJWz/4F0aOmqInp5O07oyIjSM008/nY0bN1JTU0NKygJiYyawb/+DqHr827pjVAbVHW6eLO3d2UKY08Etc0fzwa4qPtxT3d+HaIwxfRZyoSAizLx4FDUHm9mxtoxZl30ep9PFmqVPE7cgC+300PD2fubPn4/D4WDlypWICKNHf5Ompp0cPPiSf1t5CTHMTYzlf/Yd6vVzFr581mjS4yL47Rs7en11tDHGBErIhQLA+JnpZI6NZ83LuwmLiGPmpVey/b2VHKreQ8ysTBrXlBHV4mTOnDls3ryZAwcOkJ5+CfHxM9hV9Fs6Ow8/U+E/xg+npsPNf+/r3QVpUeFOvn1eDuv31vD6Jwf7+xCNMaZPQioU2ktKAO/ZwjnX5dJS387qF3cx58rriEtNY8VD9xGzcDgS7qTmpV2cc/Y5xMXF8fLLL+N2e8jN/Xfa2yvYVfQb/zZPi4vmhmHJ/LWkgo/rm3tVz5JZ2UwdEc9Pl22htrm9X4/VGGP6ImRCoW7ZMoouvoSmdd6pphmj45lxwUi2vldK8Y4GLvjqt6gq2c+qFx4l8fKxtO+tp311OZdddhkVFRUsX76chPjpjMz+CgcOPHnEbbX/fdxwMsLD+Oa2vTR2nvigs8vp4DdXT6e2uYPvL91sU1SNMUEXMqEQu3Ah4dnZHLjjDtr37gVgzhVjSR8Vx4pHthGdlMusy69m0/J/UFixnugz0ml4ez9ZzYnMmTOHdevWkZ+fz7hx/0pc3BS2bvtX6us3A95bX/x58ij2t7bzlS17afd4jlHJkSYPj+fHn5vE8m2H+N2bBYE4dGOMOWEhEwrOuDiy/+8+EGHfl75Ma0EhTpeDS74xjcjYMF65dyPj8i5jXN4c3n74fvZGFxI+Op7q5wqZlz6D8ePH8+qrr5Kfv5Fp0x4kLCyJjZv+hdq6DQCcmRjL7ydk825NA1/ZspemXpwx3Hz2aJbMHsl97xTxX/+0gWdjTPCETCgAhI8axagnHgcR9i5ZQs3zzxOdEM4Vd55OZHQYy+79hJHTbmDszNmsfPwBNna8gysrmtrnd7I4dha543N4/fXX+cfra5g8+S+4XPF89NGN7N//MKpurh+Wwq9zs3irqp7PfbSTTQ0nNsYgItxz5VRumDOS/3uniFseXU95fWuA/zaMMeazZKj9VpqXl6f5+fkntY2OQ4covftumtesJXLyZFK+/nVcs+ey8qld7P2kiqRh0cTGb2LXh68RGRPPoslfJLY8Folz8cmIctbt20hkZCRnnTWNmNjnqK1dRUxMLqNH3UZa2mJW1bZx5479VLR38vmMJG7LTmNKbBQicsy6VJUn1u7jP1/bjkOEm+eO5obZI8lOjj6p4zXGGBHZoKp5x20XiqEAoB4PdcuWUfmnP9NRUoIzIYHos8+mctRcth5MobbWgziqcbKaxqpCUqOymJm5mERSqHDUsyF2LyXtFYjAxNwmUjPWAqU4HNGkJM/HGT+Hxxun8myVkxaPMi4qggXJcZwRH83k2CiyI8OJdTm7rW1vZRO/X17Iq5tLUYUZIxM5a2wKM0YmMTYthuykaMJdIXWSZ4w5SRYKJ0g7O2lavZr6116jae06Og8dQhFqkiZQnjmTyuQptDo6cbdvx9Oxm0Snk+yYiYyIzsEdEcFO50EOOKqokgYSk8pITdtPctIBIiK9XUcNxLLOM5cPOZudjhzaJcK/7xh3E7GeZqI8rUR52oj0tOH0eHCqGxceOlqhuiKSuppIGpvDUf30TEMJc3lwudzer04PIorD4fsj4BAPCIji/crhz9l7wuJ9Lf7/mFBgH/XQNlGruf/OO/q07omGgqtPWz9FHGw6yOrS1WxwbGDXvL3sOa2J2ConWZVKZk0haXUFxNVBQnsc4ZoNjmQ6XQkUt9Wwp+YDxCnERcQx0hXJeFcKnXUZtJVMo9rpwRPViCO+HImqZXrkQWZGPI66OjkUkcIhVwrVrmSqJIUmiaYlLIoWoqmTeNy4cOPEjYvOOCfuNCeduKATtNGNNLuh2U1bu4e2DoV2D7gVcSt4AE+XrwBdM1+7LPjMcnOqs0AY+tzDAt9DEJKh8P6B93ls62OsLfPe8C45MpkJYYlc1ephmLuWlHg3GhPDwWHDaA/PpC12GI644cTHpxIdE09EdCLR0fFERUcR7oogPCycsLBwwsIiCHeF43KF4XB6Pzzx/VP8dDzh03+Y4nAc+dr3vsO3nOOMPxhjTCCEVCjUtNbw09U/ZWXxSjKiM7j99NtZFD+e8SvuQQ68QxFZPN5xFdvGXMSCOTO5ITeN2IiQ+isyxoS4kPmJV1RbxLfe+hYVzRXcNfMubpx8I2ElH8FT19Hihh+2f5PdmRfzq2umM2V4QrDLNcaYoAhoKIjIYuB/ASfwkKr++qj37wK+CnQCFcC/qOq+QNRS1VKFW908svgRpqVNg4pCeOpa6h2JXFL3XWZMO52l1063WT3GmJAWsNlHIuIECoELgBJgPbBEVbd1abMQWKeqzSLyDWCBqn7hWNs9mdlHbe42IpwR4O6AB86ls/4g5zf8lFHjJvHwzbNwOqwf3xhzajrR2UeB/LV4NrBLVXerajvwDHBF1waqulJVP73sdy2QFcB6vIEAsPpeKN/KbyNupz5yOP/zhdMtEIwxhsCGwgiguMvrEt+ynnwF+Ed3b4jIrSKSLyL5FRUVJ1dVaz188D+UD1vIA4cmcPfiiSTFhJ/cNo0x5hQRyFDo7lfvbvuqRORGIA/4bXfvq+qDqpqnqnlpaWknV1X+X6G1jt+1XUVWUhSfP+NYOWWMMaElkKFQAmR3eZ0FlB7dSETOB34MXK6qbQGsB1Tho8dpHnYmz5WmcsvcMbicNrBsjDGfCuRPxPVAjoiMEZFw4HpgWdcGIjIDeABvIJQHsBavAx9BdRGros/HIXDF6cMDvktjjBlKAhYKqtoJ3A68AWwHnlPVrSLycxG53Nfst0As8LyIbBSRZT1srn9sfRF1hnNv2STOHpdKamzE8dcxxpgQEtDrFFT1deD1o5b9pMv35wdy/59RtJLWYbPZukv4xbyMAd21McYMBaHTod5wCMq3UhDjnaZ71rjUIBdkjDGDT+iEwu53AFjeNom0uAjGpcUEtx5jjBmEQicUJlwM1z/NsoOpzB6dfNynoBljTCgKnVCIjKdu1AUU17YxZUR8sKsxxphBKXRCAdhRVg/ApGEWCsYY052QCoVtvlCYbKFgjDHdCqlQ2FHWQHJMOOlxdn2CMcZ0J6RCYW9VE2NSY2yQ2RhjehBSoVBc3cyo5Ohgl2GMMYNWyIRCW6ebsvpWsi0UjDGmRyETCiU1LajCqBQLBWOM6UnIhML+au8D3kbamYIxxvQoZEIhLsLFRVMyGJ1qt7cwxpieBPQuqYNJ3uhk8kYnB7sMY4wZ1ELmTMEYY8zxWSgYY4zxs1AwxhjjZ6FgjDHGz0LBGGOMn4WCMcYYPwsFY4wxfhYKxhhj/ERVg11Dr4hIBbCvj6unApX9WE4w2bEMTqfKsZwqxwF2LJ8apappx2s05ELhZIhIvqrmBbuO/mDHMjidKsdyqhwH2LH0lnUfGWOM8bNQMMYY4xdqofBgsAvoR3Ysg9OpciynynGAHUuvhNSYgjHGmGMLtTMFY4wxx2ChYIwxxi9kQkFEFotIgYjsEpG7g13PyRCRvSLyiYhsFJH8YNfTGyLysIiUi8iWLsuSRWS5iOz0fU0KZo0noofj+A8ROeD7XDaKyCXBrPFEiUi2iKwUke0islVEvuNbPqQ+l2Mcx5D7XEQkUkQ+FJFNvmP5mW/5GBFZ5/tMnhWR8H7fdyiMKYiIEygELgBKgPXAElXdFtTC+khE9gJ5qjrkLsgRkflAI/A3VZ3qW/YboFpVf+0L7CRV/UEw6zyeHo7jP4BGVf1dMGvrLREZBgxT1Y9EJA7YAFwJ3MwQ+lyOcRzXMcQ+FxERIEZVG0UkDHgf+A5wF/Ciqj4jIvcDm1T1//pz36FypjAb2KWqu1W1HXgGuCLINYUkVV0FVB+1+ArgMd/3j+H9hzyo9XAcQ5KqlqnqR77vG4DtwAiG2OdyjOMYctSr0fcyzPdHgUXAUt/ygHwmoRIKI4DiLq9LGKL/s/go8KaIbBCRW4NdTD/IUNUy8P7DBtKDXM/JuF1ENvu6lwZ1d0t3RGQ0MANYxxD+XI46DhiCn4uIOEVkI1AOLAeKgFpV7fQ1CcjPsVAJBelm2VDuN5urqmcAFwPf8nVlmOD7P2AccDpQBvw+uOX0jojEAi8Ad6pqfbDr6atujmNIfi6q6lbV04EsvL0dk7pr1t/7DZVQKAGyu7zOAkqDVMtJU9VS39dy4CW8/8MMZYd8/cGf9guXB7mePlHVQ75/yB7gLwyhz8XXb/0C8KSqvuhbPOQ+l+6OYyh/LgCqWgu8A5wJJIqIy/dWQH6OhUoorAdyfCP34cD1wLIg19QnIhLjG0RDRGKAC4Etx15r0FsG3OT7/ibg70Gspc8+/QHqcxVD5HPxDWr+Fdiuqv/d5a0h9bn0dBxD8XMRkTQRSfR9HwWcj3eMZCVwja9ZQD6TkJh9BOCbhvY/gBN4WFXvCXJJfSIiY/GeHQC4gKeG0rGIyNPAAry3AD4E/BR4GXgOGAnsB65V1UE9iNvDcSzA20WhwF7g65/2yQ9mInIO8B7wCeDxLf4R3v74IfO5HOM4ljDEPhcRmYZ3INmJ95f351T1575//88AycDHwI2q2tav+w6VUDDGGHN8odJ9ZIwx5gRYKBhjjPGzUDDGGONnoWCMMcbPQsEYY4yfhYIx/UBEEkXkm8Guw5iTZaFgTP9IBCwUzJBnoWBM//g1MM53v/7fBrsYY/rKLl4zph/47sr56qfPVjBmqLIzBWOMMX4WCsYYY/wsFIzpHw1AXLCLMOZkWSgY0w9UtQr4QES22ECzGcpsoNkYY4yfnSkYY4zxs1AwxhjjZ6FgjDHGz0LBGGOMn4WCMcYYPwsFY4wxfhYKxhhj/P4/saajZn61trUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.clf()                                                                     \n",
    "for i, y0 in enumerate(np.linspace(0.1,1.2,20)):\n",
    "    ode.set( ics = { 'y': y0 } )                # Initial condition\n",
    "    tmp = ode.compute('pol%3i' % i).sample()    \n",
    "    plt.plot(tmp['t'], tmp['y'])\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('y')\n",
    "plt.title(ode.name + ' multi ICs')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LP Point found \n",
      "LP Point found \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/dijkstra1_henk/anaconda3/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:107: MatplotlibDeprecationWarning: Using pyplot.axes(ax) with ax an Axes argument is deprecated. Please use pyplot.sca(ax) instead.\n",
      "  warnings.warn(message, mplDeprecation, stacklevel=1)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Bifurcation diagram AMOC')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABakAAAE3CAYAAACgmc3qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3WmYVdWd/v3vj6kYRNHgXExiAVLKWIqIKKLGOUg6Jqgdh0TtxL9Da2LUJzHRNppotGOi0ShRYxKcEiWC0dbGIQ6AUqAgICLFoCiIoEZFBoH1vKhDdVEWo3Vq1/D9XNe56uy119nnPofCF7eLtSOlhCRJkiRJkiRJWWiSdQBJkiRJkiRJUuNlSS1JkiRJkiRJyowltSRJkiRJkiQpM5bUkiRJkiRJkqTMWFJLkiRJkiRJkjJjSS1JkiRJkiRJyowltSRJkiRJkiRpq0TEsxHxRkRMjYgXI6J7bnxUbnx6RNwVEc03dy1LakmSJEmSJEnStjg1pdQbuAf4VW5sFNAD2A9oBZy1uYtYUkuSJEmSJEmSvozngL0BUkqPpRzgZaBwcy+2pJYkSZIkSZIkfRknAK9VHsht8/Ft4H829+JmeQolSZIkSZIkSWrYRkXECmA+cH6Vc7cCz6WUnt/cRSypJUmSJEmSJEnb4tSUUmnVwYj4GbAz8B9bchFLakmSJEmSJElSjYiIs4CjgMNTSuu25DXuSS1JkiRJkiRJqim/B3YFJkTEqxHx0829IMpvsihJkiRJkiRJUu1zJbUkSZIkSZIkKTOW1JIkSZIkSZKkzFhSS5IkSZIkSZIyY0ktSZIkSZIkScqMJbUkSZIkSZIkKTOW1JIkSZIkSZKkzFhSS5IkSZIkSZIyY0ktSZIkSZIkScqMJbUkSZIkSZIkKTOW1JIkSZIkSZKkzFhSS5IkSZIkSZIyY0ktSZIkSZIkScqMJbUkSZIkSZIkKTOW1JIkSZIkSZKkzFhSS5IkSZIkSZIyY0ktSZIkSZIkScqMJbUkSZIkSZIkKTOW1JIkSZIkSZKkzFhSS5IkSZIkSZIyY0ktSZIkSZIkScqMJbUkSZIkSZIkKTOW1JIkSZIkSZKkzFhSS5IkSZIkSZIyY0ktSZIkSZIkScqMJbW2WkT8PiKuqHT8/Yh4LyI+jYiv1KVseXyfIRGxsNLxjIgYku/3lSRJkiRJkhoaS2p9QUTMj4gVudL5w4j4R0R0WH8+pfS9lNLVubnNgf8GvppS2i6ltKwWc54RES9UHqucrTallIpTSs/W9vvWlIjoEhHrIuLWas6l3P+EaFZprFlELImIVGXu8RHxckQsj4hlETEqIgqrzNk9Iu6MiEUR8UlEzIqIqyKiTf4+oSRJkiRJkuoqS2ptzAkppe2A3YH3gJs3Mm9XoCUwY1vepHLxqS+qxe/nNOBDYEREFFRz/iPgmErHx+bmV4iIbwD3Ar8B2gPFwCrghYjYMTdnJ2AC0AoYmFJqCxwJtAO61uQHkiRJkiRJUv1gSa1NSimtBP4G9Fw/FhF/jIifR0Q34I3c8EcR8XREdM6tvK286vbZiDgr9/yMiHgxIn4dER8AV+bGz46I13Mra2dGRL/c+GURUVZpfHhufB/g98DA3Irvjypnq/TeZ0fEnIj4ICLGRMQelc6liPheRLyZWzH+u4iI6r6HiGiVu/aHETET2L/K+fkRcUTu+QERMSEiPsqtFr4lIlpUmvvViHgjIv4VEbdGxD839f1ERNfcd7ssIpbmVie3q/Lel0TEtNwK5jsjYteIeDz3vY1bXxJvwmnAT4DPgROqOf/n3JzK8/9UKUMANwI/TymNSimtSCktBs4CPgUuyk29GPgE+PeU0nyAlNLbKaULU0rTNpNRkiRJkiRJDZAltTYpIloD3wImVj2XUppN+WpZgHYppaFbeNkBwFxgF+CaiDiJ8rL6NGB74GvA+m1DyoDBwA7AVcBfImL3lNLrwPeACbltRtpRRUQMBX4BfJPyFeELgPurTDue8sK5d27eURvJ/DPKV/p2zc05fROfby3lpWx7YCBwOHBuLlN7ykv/y4GvUF7yH1Tl9Rt8P0DkPscewD5AB3LlfiX/RvmK5G6Ul8yPA/9fLkMT4IKNhY2IwUAh5d/Ng2xYRq/3d+CQiGiXK8gHA49UOt8d6Aj8tfKLUkrrgIdy2QCOAB7OjUuSJEmSJEmW1Nqov+dWJ39MecH4qxq89rsppZtTSmtSSisoX217fUppUio3J6W0ACCl9NeU0rsppXUppQeAN4EDtvB9TgXuSilNSSmtorwYHhgRnSvN+WVK6aOU0lvAM0CfjVzrm8A1KaUPUkpvA7/d2JumlCanlCbmPt984Hbg0NzpY4EZKaWHU0prctdZvKnvJ/d9/G9KaVVK6X3K9wA/tMprbk4pvZdSegd4HngppfRK7nOPBvpu/GvidODxlNKHlG/XcUxE7FJlzkpgLOX/w2IEMCY3tl773M9F1Vx/UaXzX9nIHEmSJEmSJDVSltTamBNzq5MLgPOAf0bEbjV07berHHegfMX0F0TEaRHxam7rjI+Affm/wnNz9qB89TQAKaVPKV+hvWelOZUL4s+A7TZxrcq5F2xkHhHRLSIejYjFEfExcG2lzBtcJ6WUgIVVLrHB9xMRu0TE/RHxTu56f+GL38F7lZ6vqOa42s8VEa2Ak4BRuTwTgLeAU6qZ/ifKV1lvsNVHztLcz92red3ulc4v28gcSZIkSZIkNVKW1NqklNLalNLDlG9hcfAWvGR57mfrSmNVy+1U5fhtqrlpXkR0AkZSXpJ/JVeaT6d8+4vqrlPVu0CnStdrQ/lK3nc287rqLKK8TF+v4ybm3gbMAopSSttTvu3G+syLKN9aY32mqHycU/Vz/SI31it3vX+vdL0vazjlW6zcmivVF1Ne4le35cfzlBfMuwIvVDn3BuVl+0mVByOiCeVbkTyVGxoHDM+NS5IkSZIkSZbU2rQoNwzYEXh9c/Nz21G8A/x7RDSNiO9QTQFdxR+AH0ZE/9z77Z0rqNtQXs6+n8tyJuUrqdd7DyisfFPCKu4FzoyIPhFRQPmK5pfW37BvKz0IXB4RO0ZEIXD+Jua2pXyblE8jogfw/Urn/gHsFxEnRvnNJf8fXyzxq7vep5TfnHJP4JJtyL8xpwN3AftRvtVJH2AQ0Cci9qs8Mbfq+wTga7nnVc/9EPhJRJySu9HkbpT/2W4P/Do39b9zx/fk/oyJiD0j4r8jolcNfi5JkiRJkiTVE5bU2pixEfEp5WXrNcDpKaUZW/jasykvUpdRfmPF8ZuanFL6a+497gU+ofwmfTullGYCNwITKC+k9wNerPTSp4EZwOKIWEoVKaWngCsov3HfIsrL8hFb+BmquoryLT7mAU8Cf97E3B9Svl3GJ5SvBH+gUqallK82vp7y76cnUAqs2sx79wP+RXnJ/fA2foYN5Arvw4GbUkqLKz0mA/9DNTeHTCnN2NjvQW7P8G9TftPIpcBMoBUwKKW0LDfnA8pvFPk58FJEfEL5Kut/AXNq4nNJkiRJkiSpfokqCyIl1aLcthcLgVNTSs9knUeSJEmSJEmqba6klmpZRBwVEe1yW5Cs3696YsaxJEmSJEmSpExYUku1byBQRvmWGCcAJ6aUVmQbSZIkSZIkScqG231IkiRJkiRJkjLjSmpJkiRJkiRJUmaaZR1ga7Vv3z517tw56xiSVK9Mnjx5aUpp56xzSJIkSZIkVVXvSurOnTtTWlqadQxJqlciYkHWGSRJkiRJkqrjdh+SJEmSJEmSpMxYUkuSJEmSJEmSMmNJLUmSJEmSJEnKjCW1JEmSJEmSJCkzltSSJEmSJEmSpMxYUkuSJEmSJEmSMmNJLUmSJEmSJEnKjCW1JEmSJEmSJCkzltSSJEmSJEmSpMzkraSOiLsiYklETN/MvP0jYm1EfCNfWSRJkiRJkiRJdVM+V1L/ETh6UxMioilwHfBEHnNIkiRJkiRJkuqovJXUKaXngA82M+184CFgSb5ySJIkSZIkSZLqrsz2pI6IPYHhwO+zyiBJkiRJkiRJylaWN068Cbg0pbR2cxMj4pyIKI2I0vfff78WokmSJEmSJEmSakOzDN+7BLg/IgDaA8dGxJqU0t+rTkwp3QHcAVBSUpJqNaUkSZIkSZIkKW8yK6lTSl3WP4+IPwKPVldQS5IkSZIkSZIarryV1BFxHzAEaB8RC4GfAc0BUkruQy1JkiRJkiRJyl9JnVI6eSvmnpGvHJIkSZIkSZKkuivLGydKkiRJkiRJkho5S2pJkiRJkiRJUmYsqSVJkiRJkiRJmbGkliRJkiRJkiRlxpJakiRJkiRJkpQZS2pJkiRJkiRJUmYsqSVJkiRJkiRJmbGkliRJkiRJkiRlxpJakiRJkiRJkpSZZlkHkCTVvJQSS5YsoaysjLlz52YdR5IkSZIkaaMsqSWpnlq3bh3vvPMOb775Jm+++SZz5syhrKysopguKCiga9eudO3aNeuokiRJkiRJG2VJLUl1WEqJxYsXVxTR6x+zZ8+mrKyMdu3aUVRURFFREXvvvTcHHHBARTG9ww47VFznvvvuy/BTSJIkSZIkbZwltSRlLKXE0qVLmT179hfK6Dlz5tCqVauKIrqoqIhvfetbFaV027Zts44vSZIkSZL0pVhSS1It+fzzzykrK+ONN95g1qxZFT9nzZpFSolu3bpRVFREt27dGD58eEUR3a5du6yjS5IkSZIk5Y0ltSTVsA8++OALJfQbb7zB/PnzKSwspEePHnTv3p2BAwdy5pln0r17d3beeWciIuvokiRJkiRJtc6SWpK2QUqJd955hxkzZjBjxgxef/31ikJ69erVdO/enR49etCjRw9OO+00evTowd57701BQUHW0SVJkiRJkuoUS2pJ2oTKZfTMmTMrSumZM2fSqlUriouL6dmzJ/369eOUU06hR48e7Lbbbq6KliRJkiRJ2kKW1JJEeRn97rvvVpTQ64vomTNn0rJlS3r27ElxcTElJSWcdtppFBcX85WvfCXr2JIkSZIkSfWeJbWkRueTTz7htddeY9q0aUybNo2pU6cyY8YMCgoKKC4utoyWJEmSJEmqRZbUkhqsdevWMW/evIoiev3PxYsX07NnT3r37k2vXr345je/yb777kv79u2zjixJkiRJktTo5K2kjoi7gOOBJSmlfas5fypwae7wU+D7KaWp+cojqWFbvzq6chn92muvseOOO1aU0SNGjODaa6+lqKiIpk2bZh1ZkiRJkiRJ5Hcl9R+BW4A/beT8PODQlNKHEXEMcAcwII95JDUQH3zwAVOmTKG0tJTJkyczZcoUFi9eTHFxMb169aJ3796MGDGCXr16seOOO2YdV5IkSZIkSZuQt5I6pfRcRHTexPnxlQ4nAoX5yiKp/qpaSE+ePJmlS5fSt29f+vfvz4knnsjVV1/t6mhJkiRJkqR6qq7sSf1d4PGsQ0jK1pYU0j//+c8pKiqiSZMmWceVJEmSJElSDci8pI6IwygvqQ/exJxzgHMAOnbsWEvJJOXTypUrmTx5MhMmTOCll16ykJYkSZIkSWqkIqWUv4uXb/fxaHU3Tsyd7wWMBo5JKc3ekmuWlJSk0tLSGssoqXYsXLiQCRMmMGHCBMaPH89rr73GPvvsw0EHHcSAAQMoKSmxkM6jiJicUirJOockSZIkSVJVma2kjoiOwMPAt7e0oJZUP6xevZpXX321opCeMGECK1asYODAgRx00EFcd911lJSU0KZNm6yjSpIkSZIkKWN5K6kj4j5gCNA+IhYCPwOaA6SUfg/8FPgKcGtEAKxxlZ9UP7333nsbrJJ+5ZVX6Nq1KwMHDuS4447j5z//OXvvvTe5v+uSJEmSJElShbxu95EPbvchZW/JkiU8++yzPP300zz99NO8//77HHjggRUrpQ844AC23377rGOqErf7kCRJkiRJdVXmN06UVPd99NFHPPfccxWl9FtvvcUhhxzC0KFDOffcc9l3333dS1qSJEmSJEnbxJJa0hcsX76cF198saKUfv311xk4cCBDhw7lD3/4A/369aNZM//zIUmSJEmSpC/PlkkSq1at4qWXXqoopadMmUK/fv0YOnQoN9xwAwMGDKCgoCDrmJIkSZIkSWqALKmlRuqdd95h7NixjBkzhueff5599tmHoUOH8pOf/IRBgwbRpk2brCNKkiRJkiSpEbCklhqJlBLTp0/nkUce4ZFHHqGsrIxjjz2WM888k3vvvZd27dplHVGSJEmSJEmNkCW11ICtWbOGF154oaKYXrduHcOGDeO6665j8ODBNG/ePOuIkiRJkiRJauQsqaUG5tNPP+WJJ57gkUce4bHHHqNTp04MGzaM0aNH06tXLyIi64iSJEmSJElSBUtqqQFYtGgRY8eO5ZFHHuH555/nwAMPZNiwYVxzzTV06NAh63iSJEmSJEnSRllSS/XUypUrefjhhxk5ciSvvvoqRx99NN/+9re599572WGHHbKOJ0mSJEmSJG0RS2qpnpkxYwYjR47kL3/5C/369ePcc8/la1/7GgUFBVlHkyRJkiRJkraaJbVUDyxfvpwHH3yQkSNHsmDBAs4880wmTZpEly5dso4mSZIkSZIkfSmW1FIdNmXKFEaOHMkDDzzAoEGDuOyyyzj22GNp1sy/upIkSZIkSWoYbLqkOubjjz/m3nvvZeTIkSxbtozvfve7TJs2jcLCwqyjSZIkSZIkSTXOklqqA1JKTJw4kZEjRzJ69GgOP/xwrr32Wo444giaNm2adTxJkiRJkiQpbyyppYyNHz+eCy+8kH/961+cddZZzJo1i1133TXrWJIkSZIkSVKtsKSWMrJo0SJ+9KMf8cwzz3D99ddz8sknExFZx5IkSZIkSZJqVZOsA0iNzerVq/nVr37Ffvvtx5577smsWbM45ZRTLKglSZIkSZLUKLmSWqpFTzzxBBdeeCFdu3ZlwoQJFBUVZR1JkiRJkiRJypQltVQL5s6dy0UXXcSMGTO46aabOP7447OOJEmSJEmSJNUJbvch5dFnn33GFVdcwQEHHMCBBx7IjBkzLKglSZIkSZKkSvJWUkfEXRGxJCKmb+R8RMRvI2JOREyLiH75yiLVtpQSDz74ID169GDOnDm8+uqrXH755RQUFGQdTZIkSZIkSapT8rndxx+BW4A/beT8MUBR7jEAuC33U6rXpk+fzgUXXMDSpUv585//zKGHHpp1JEmSJEmSJKnOyttK6pTSc8AHm5gyDPhTKjcRaBcRu+crj1Qbrr/+eoYOHcq//du/MWXKFAtqSZIkSZIkaTOyvHHinsDblY4X5sYWZRNH+nJee+01brjhBl599VX22GOPrONIkiRJkiRJ9UKWN06MasZStRMjzomI0ogoff/99/McS9p6KSXOP/98rrzySgtqSZIkSZIkaStkWVIvBDpUOi4E3q1uYkrpjpRSSUqpZOedd66VcNLWeOCBB/joo4/4j//4j6yjSJIkSZIkSfVKliX1GOC0KHcg8K+Uklt9qN759NNPueSSS7jlllto2rRp1nEkSZIkSZKkeiVve1JHxH3AEKB9RCwEfgY0B0gp/R54DDgWmAN8BpyZryxSPl1zzTUMGTKEgw8+OOsokiRJkiRJUr2Tt5XUKaWTU0q7p5Sap5QKU0p3ppR+nyuoSeX+X0qpa0ppv5RSab6ySPkye/ZsRo4cyfXXX8+QIUMoLd3w1/jZZ59lhx12oG/fvuyzzz5cddVVAPzv//4v/fv3Z7/99qN///48/fTTWcSXJEmSJEmSMpe3ldRSQ5dS4sILL+Syyy5j99133+i8wYMH8+ijj7J8+XL69OnD8ccfT/v27Rk7dix77LEH06dP56ijjuKdd96pxfSSJEmSJElS3ZDlntRSvTZ27Fjmz5/PBRdcsEXz27RpQ//+/SkrK6Nv377sscceABQXF7Ny5UpWrVqVz7iSJEmSJElSnWRJLW2DFStW8J//+Z/89re/pUWLFlv0mmXLljFx4kSKi4s3GH/ooYfo27cvBQUF+YgqSZIkSZIk1Wlu9yFtgxtuuIF+/fpx5JFHbnbu888/T9++fWnSpAmXXXbZBiX1jBkzuPTSS3nyySfzGVeSJEmSJEmqsyyppa00f/58brrpJqZMmbJF89fvSV3VwoULGT58OH/605/o2rVrTceUJEmSJEmS6gW3+5C2xKhR0LkzNGlC6549GXnYYXTq1GmbL/fRRx9x3HHH8Ytf/IJBgwbVXE5JkiRJkiSpnrGkljZn1Cg45xxYsABSYpcVKxj++OPl41Ucd9xxFBYWUlhYyEknnbTRS95yyy3MmTOHq6++mj59+tCnTx+WLFmSz08hSZIkSZIk1UmRUso6w1YpKSlJpaWlWcdQY9K5c3lBXVWnTjB/fm2nkbZJRExOKZVknUOSJEmSJKkqV1JLm/PWW1s3LkmSJEmSJGmLWVJLm9Ox49aNS5IkSZIkSdpiltTS5lxzDbRuveFY69bl45IkSZIkSZK+FEtqaXNOPRXuuKN8D+oI6NSJdMcd5eOSJEmSJEmSvhRLamlLnHpq+U0S163j1h/9iKvefDPrRJIkSZIkSVKDYEktbaVDDz2U22+/nc8//zzrKJIkSZIkSVK9Z0ktbaXi4mK6devG3//+96yjSJIkSZIkSfWeJbW0DS699FLWrl2bdQxJkiRJkiSp3muWdQCpPjr22GNJKbFq1SoKCgqyjiNJkiRJkiTVW66klrbRLbfcwiWXXJJ1DEmSJEmSJKles6SWttGJJ57IqFGjWL58edZRJEmSJEmSpHrLklraRh06dGDw4MHce++9WUeRJEmSJEmS6q28ltQRcXREvBERcyLismrOd4yIZyLilYiYFhHH5jOPVNMuvfRSdt9996xjSJIkSZIkSfVW3krqiGgK/A44BugJnBwRPatM+wnwYEqpLzACuDVfeaR8GDhwIMceeyzPPvts1lEkSZIkSZKkeimfK6kPAOaklOamlFYD9wPDqsxJwPa55zsA7+Yxj5QXixcv5tvf/jZ333131lEkSZIkSZKkeqdZHq+9J/B2peOFwIAqc64EnoyI84E2wBF5zCPlxR577MG4ceMYOnQoBQUFnHLKKVlHkiRJkiRJkuqNfK6kjmrGUpXjk4E/ppQKgWOBP0fEFzJFxDkRURoRpe+//34eokpfTvfu3XnyySeZNm1a1lEkSZIkSZKkeiWfJfVCoEOl40K+uJ3Hd4EHAVJKE4CWQPuqF0op3ZFSKkkpley88855iit9OcXFxfzyl79kypQpjB07Nus4kiRJkiRJUr2Qz5J6ElAUEV0iogXlN0YcU2XOW8DhABGxD+UltUulVa+tW7eO7373uzz++ONZR5EkSZIkSZLqvLyV1CmlNcB5wBPA68CDKaUZEfFfEfG13LQfAGdHxFTgPuCMlFLVLUGkeqWkpIQxY8Zw+umnM27cuKzjSJIkSZIkSXVa1LdOuKSkJJWWlmYdQ9qs559/nrZt29KnT5+so0hExOSUUknWOSRJkiRJkqra7ErqiDgvInasjTBSQzJ48GD69OnDD37wA8aPH591HEmSJEmSJKlO2pLtPnYDJkXEgxFxdEREvkNJDcmRRx7JiSeeyKRJk7KOIkmSJEmSJNU5my2pU0o/AYqAO4EzgDcj4tqI6JrnbFKDcPTRR3PnnXdy/PHHM2/evKzjSJIkSZIkSXXKFt04MXczw8W5xxpgR+BvEXF9HrNJDcYJJ5zAY489RqdOnfj000+zjiNJkiRJkiTVGVuyJ/UFETEZuB54EdgvpfR9oD/wb3nOJzUY/fv3Z+3atZSUlDBy5Mis40iSJEmSJEl1QrMtmNMe+HpKaUHlwZTSuog4Pj+xpIapefPmjBkzhhNOOIHp06dz44030qzZlvw1lCRJkiRJkhqmLdmT+qdVC+pK516v+UhSw9atWzdeeukllixZwrvvvpt1HEmSJEmSJClTW7QntaSa1a5dO+677z46dOjARRddxKxZs7KOJEmSJEmSJGXCklrKUESw3377ccghh/DEE09kHUeSJEmSJEmqdZbUUsa+853v8PDDD3PGGWdQWlqadRxJkiRJkiSpVllSS3XAwQcfzNSpU+nfvz9Tpkxh9erVWUeSJEmSJEmSaoUltVRH7LLLLkQEv/nNbzj88MNZsmRJ1pEkSZIkSZKkvLOkluqYu+++myFDhjBgwACmT5+edRxJkiRJkiQpr5plHUDShpo0acLVV19Nr169aNu2LWvXrqVp06ZZx5IkSZIkSZLywpJaqqNOOukkAM466yw6duzIT37yE5o08R8/SJIkSZIkqWGx8ZLquKuvvprHH3+cb37zmyxfvjzrOJIkSZIkSVKNsqSW6rjdd9+dZ599lrZt2zJ69Ois40iSJEmSJEk1yu0+pHqgoKCAu+66C4C///3v7LTTThxyyCEZp5IkSZIkSZK+PFdSS/VERBARtGnThpNOOonbb78960iSJEmSJEnSl+ZKaqmeOfLII3nhhRcYNmwYy5cv5+KLL846kiRJkiRJkrTN8lpSR8TRwG+ApsAfUkq/rGbON4ErgQRMTSmdks9MUkNQVFTExIkT+eyzz1i2bBkpJdq3b591LEmSJEmSJGmr5W27j4hoCvwOOAboCZwcET2rzCkCLgcGpZSKgf/MVx6podl+++3ZbbfdeOSRRzjggAN47bXXso4kSZIkSZIkbbV8rqQ+AJiTUpoLEBH3A8OAmZXmnA38LqX0IUBKaUke80gN0ne+8x1atmzJ0KFDuf322/n617+edSRJkiRJkiRpi+Xzxol7Am9XOl6YG6usG9AtIl6MiIm57UEkbaVTTjmFxx9/nBUrVmQdRZIkSZIkSdoq+VxJHdWMpWrevwgYAhQCz0fEvimljza4UMQ5wDkAHTt2rPmkUgNQUlJCSUkJ999/P4sWLeKiiy7KOpIkSZIkSZK0WflcSb0Q6FDpuBB4t5o5j6SUPk8pzQPeoLy03kBK6Y6UUklKqWTnnXfOW2CpITjooIP43e9+xw033JB1FEmSJEmSJGmz8llSTwKKIqJLRLQARgBjqsz5O3AYQES0p3z7j7l5zCQ1eB07duTZZ5/l9ttv55577sk6jiRJkiRJkrTy3SzSAAAZa0lEQVRJedvuI6W0JiLOA54AmgJ3pZRmRMR/AaUppTG5c1+NiJnAWuCSlNKyfGWSGovCwkL++c9/0qpVK1asWEGrVq2yjiRJkiRJkiRVK1Kquk103VZSUpJKS0uzjiHVGyNGjKBnz5789Kc/zTqKMhQRk1NKJVnnkCRJkiRJqiqfN06UVAfcdNNNHH744axdu5Yrr7ySiOruaSpJkiRJkiRlI597UkuqA3bbbTeeeeYZHnnkEaZOnZp1HEmSJEmSJGkDltRSI7DLLrswadIk+vTpw8svv0x92+ZHkiRJkiRJDZcltdRING/enM8//5zvfe97XHrppRbVkiRJkiRJqhMsqaVGpHnz5owbN45x48bxgx/8wKJakiRJkiRJmbOklhqZnXbaiaeeeopPPvmEFStWZB1HkiRJkiRJjZwltdQI7bjjjowcOZLPPvuMdevWZR1HkiRJkiRJjZgltdSIDR8+nEcffTTrGJIkSZIkSWrELKmlRuzcc8/lxhtvzDqGJEmSJEmSGjFLaqkR+8Y3vsG8efMoLS3NOookSZIkSZIaKUtqqRFr3rw599xzD7vtttsWv2bIkCFfKLWfffZZdthhB/r27cs+++zDVVddBcCyZcs47LDD2G677TjvvPNqNLskSZIkSZIaBktqqZE77LDDWLlyJYsXL/5S1xk8eDCvvPIKpaWl/OUvf2Hy5Mm0bNmSq6++mhtuuKGG0kqSJEmSJKmhsaSWxK233lpje1O3adOG/v37U1ZWRps2bTj44INp2bJljVxbkiRJkiRJDY8ltSQuvPBC7rrrLj7++OMvfa1ly5YxceJEiouLayCZJEmSJEmSGjpLakl06tSJI488kjvvvHObr/H888/Tt29fvvrVr3LZZZdZUkuSJEmSJGmLNMs6gKS64Ve/+hWtW7fe5tcPHjyYRx99tAYTSZIkSZIkqTFwJbUkADp06MA777zDiy++mHUUSZIkSZIkNSKW1JIqrLr7bvYaOpTUpAl07gyjRlU777jjjqOwsJDCwkJOOumkTV6zc+fOXHzxxfzxj3+ksLCQmTNn5iG5JEmSJEmS6qtIKWWdYauUlJSk0tLSrGNIDc+oUaRzziE+++z/xlq3hjvugFNPzS6XakRETE4plWSdQ5IkSZIkqSpXUksq9+Mfb1hQA3z2Gfz4x9nkkSRJkiRJUqOQ15I6Io6OiDciYk5EXLaJed+IiBQRrvKTsvLWW1s3LkmSJEmSJNWAvJXUEdEU+B1wDNATODkielYzry1wAfBSvrJI2gIdO27duCRJkiRJklQD8rmS+gBgTkppbkppNXA/MKyaeVcD1wMr85hF0uZccw2rmzXbcKx1a7jmmmzySJIkSZIkqVHIZ0m9J/B2peOFubEKEdEX6JBSejSPOSRtgbcPOYTzCgpYs+eeEAGdOnnTREmSJEmSJOVds81P2WZRzViqOBnRBPg1cMZmLxRxDnAOQEe3HpDy4sorr2TXCy6g2bXXZh1FkiRJkiRJjUg+S+qFQIdKx4XAu5WO2wL7As9GBMBuwJiI+FpKqbTyhVJKdwB3AJSUlCQk1aiZM2cyduxYZs+enXUUSZIkSZIkNTL53O5jElAUEV0iogUwAhiz/mRK6V8ppfYppc4ppc7AROALBbWk/Pvxj3/Mj370I9q1a5d1FEmSJEmSJDUyeVtJnVJaExHnAU8ATYG7UkozIuK/gNKU0phNX0FSbZgwYQKTJ0/mvvvuyzqKJEmSJEmSGqF8bvdBSukx4LEqYz/dyNwh+cwi6YtSSlx22WVceeWVtGzZMus4kiRJkiRJaoTyud2HpDru8ccf5/333+e0007LOookSZIkSZIaKUtqqZFat24dl19+Oddeey3NmuX1H1VIkiRJkiRJG2VJLTVS9957L23atGHYsGFZR5EkSZIkSVIj5vJJqRFasWIFV1xxBffccw8RkXUcSZIkSZIkNWKupJYameeee46+ffsydOhQDjnkkKzjSJIkSZIkqZFzJbXUSHz44YdceumlPPbYY9x8880MHz4860iSJEmSJEmSK6mlhi6lxAMPPEBxcTHNmzdnxowZFtSSJEmSJEmqM1xJLTVgCxYs4Nxzz2XBggU89NBDDBw4MOtIkiRJkiRJ0gZcSS01QGvWrOHXv/41/fv3Z9CgQUyZMsWCWpIkSZIkSXWSK6mlBuaVV17h7LPPZvvtt2fChAkUFRVlHUmSJEmSJEnaKFdSSw3E8uXLueSSSzj66KM577zzeOqppyyoJUmSJEmSVOdZUksNwP/8z/+w7777snjxYqZPn84ZZ5xBRGQdS5IkSZIkSdost/uQ6qmUEhMmTODXv/41kydP5vbbb+erX/1q1rEkSZIkSZKkrWJJLdUzn3zyCaNGjeK2225jxYoVfO973+Oee+6hdevWWUeTJEmSJEmStpoltVRPTJ8+ndtuu4377ruPww47jBtvvJGhQ4fSpIm79kiSJEmSJKn+sqSW6rBVq1bx0EMPcdtttzF37lzOPvtsXnvtNfbcc8+so0mSJEmSJEk1wpJaqoPmzZvH7bffzt13302vXr246KKLOOGEE2jevHnW0SRJkiRJkqQaZUkt1RFr167l8ccf57bbbuOll17itNNO4/nnn6dbt25ZR5MkSZIkSZLyxpJaytjs2bP561//ysiRI9l11135/ve/z9/+9jdatWqVdTRJkiRJkiQp7yyppVqWUmLKlCmMHj2a0aNH8+GHH3LiiSfy0EMP0b9//6zjSZIkSZIkSbUqryV1RBwN/AZoCvwhpfTLKucvBs4C1gDvA99JKS3IZyYpC2vXruWFF16oKKZbtGjB8OHD+cMf/sCAAQNo0qRJ1hElSZIkSZKkTOStpI6IpsDvgCOBhcCkiBiTUppZadorQElK6bOI+D5wPfCtfGWSatPKlSsZN24co0ePZuzYsRQWFjJ8+HD+8Y9/UFxcTERkHVGSJEmSJEnKXD5XUh8AzEkpzQWIiPuBYUBFSZ1SeqbS/InAv+cxj5R3H3/8Mf/4xz8YPXo0Tz75JL169eLrX/86V1xxBZ07d846niRJkiRJklTn5LOk3hN4u9LxQmDAJuZ/F3g8j3mkGpdSoqysjHHjxjFmzBheeOEFBg8ezPDhw7nlllvYZZddso4oSZIkSZIk1Wn5LKmr28sgVTsx4t+BEuDQjZw/BzgHoGPHjjWVT9om7733Hk899RRPPfUU48aNY82aNRx++OGcfvrp3H///Wy//fZZR5QkSZIkSZLqjXyW1AuBDpWOC4F3q06KiCOAHwOHppRWVXehlNIdwB0AJSUl1RbdUr588skn/POf/2TcuHE89dRTLFy4kEMPPZQjjjiCH/7wh/To0cP9pSVJkiRJkqRtlM+SehJQFBFdgHeAEcAplSdERF/gduDolNKSPGaRttjq1auZOHFiRSk9depUBgwYwOGHH86dd95Jv379aNYsn391JEmSJEmSpMYjb01bSmlNRJwHPAE0Be5KKc2IiP8CSlNKY4BfAdsBf82tRH0rpfS1fGWSqrNmzRqmTp3KM888w1NPPcWLL75I9+7dOeKII7jqqqsYNGgQrVq1yjqmJEmSJEmS1CBFSvVr94ySkpJUWlqadQzVYx999BETJ07kxRdfZPz48UyaNIkOHTowZMgQjjjiCIYMGcKOO+6YdUypRkXE5JRSSdY5JEmSJEmSqnLPAjVoKSXmzJnD+PHjGT9+PC+++CILFixg//3356CDDuKHP/whBx54oKW0JEmSJEmSlBFLajUoK1euZPLkyRWrpMePH09BQQGDBg3ioIMO4pxzzqFXr140b94866iSJEmSJEmSsKRWPZZSoqysjNLSUiZNmsSECROYOnUq++yzDwcddBAnn3wyN998Mx06dMg6qiRJkiRJkqSNsKRWvZBSYsGCBZSWllY8Jk+eTNu2bSkpKaGkpIRrrrmG/fffn+222y7ruJIkSZIkSZK2kCW16pyUEgsXLvxCId2iRYuKQvriiy+mf//+7LrrrlnHlSRJkiRJkvQlWFIrU+sL6VdffXWDUjqlxP77709JSQnnnXce/fv3Z4899sg6riRJkiRJkqQaZkmtWrNq1SpmzJjB1KlTKx7Tpk2jWbNm9OnTh5KSEs466yx+//vfU1hYSERkHVmSJEmSJElSnllSKy8WL168QRk9depUysrK6Nq1K71796Z3794cc8wx9O7dm9122y3ruJIkSZIkSZIyYkmtL2XVqlXMnj2badOmbVBIr169uqKMPuKII/jBD35Az549admyZdaRJUmSJEmSJNUhltTaIp9//jmzZ89mxowZGzzmz59P586d2W+//ejduzfnn38+vXv3drsOSZIkSZIkSVvEklobWLNmDXPmzPlCGV1WVkaHDh0oLi6muLiYb3zjG/zsZz+jW7duFBQUZB1bkiRJkiRJUj1lSd1IrV69mjlz5jBr1ixef/31ijL6zTffZPfdd68oo0844QQuv/xyunfvTqtWrbKOLUmSJEmSJKmBsaRu4D744ANmzZr1hcdbb71Fp06d6NGjB927d+eoo47i4osvZp999qFNmzZZx5YkSZIkSZLUSFhSNwBr165lwYIF1ZbRq1atokePHhWPM888kx49etC1a1datGiRdXRJkiRJkiRJjVyDKqmHDBnCokWLaNmyJdtttx133XUX3bt355ZbbuGmm26irKyM999/n/bt22cddaullFiyZAmzZ8/mzTffZPbs2RXPy8rK2HnnnSuK6L59+3LyySfTo0cPdtttN29gKEmSJEmSJKnOalAlNcCoUaMoKSnhjjvu4JJLLmHMmDEMGjSI448/niFDhmQdb7M+/PDDihK66s8WLVrQrVs3ioqK6NatGyNGjKBbt27svffebLfddllHlyRJkiRJkqSt1uBK6vUOOeQQbrrpJgD69u2bcZr/k1Ligw8+YM6cOZSVlTFnzpyKx5tvvsmqVasqSuiioiKOO+44ioqKKCoqYqeddso6viRJkiRJkiTVqAZbUo8dO5b99tsvk/dOKbFo0SLKysq+UESXlZWRUmLvvfeueAwdOpSzzz6boqIidt11V7fnkCRJkiRJktRoNLiS+tRTT6VVq1Z07tyZm2++OW/vs3LlSubNm8fcuXOZO3cuZWVlFT/nzZtH27Zt2WuvvSqK6BNOOKHi+U477WQRLUmSJEmSJEk0wJJ6/Z7UX9a6detYtGhRRRFd9efSpUvp2LEjXbt2Za+99mKvvfZiyJAh7LXXXnTp0oW2bdvWwKeRJEmSJEmSpIYtryV1RBwN/AZoCvwhpfTLKucLgD8B/YFlwLdSSvPzmWm9lBJLly5l/vz5zJs37ws/FyxYQLt27ejSpUtF8TxkyBDOPPNMunTpQmFhIU2bNq2NqJIkSZIkSZLUYOWtpI6IpsDvgCOBhcCkiBiTUppZadp3gQ9TSntHxAjgOuBbm7ru5MmT6dy5M9dccw2nnnrqRuetL6EXLFjArbfeysMPP8zHH39MYWEhrVq14vPPP6dFixZ06dKFzp0706VLF/bdd1+OP/74irHWrVt/6e9BkiRJkiRJkrRxkVLKz4UjBgJXppSOyh1fDpBS+kWlOU/k5kyIiGbAYmDntIlQEZEAWrduzXXXXUf//v2ZP38+CxYsqHisPy4oKKBTp04VpXPlR6dOndhhhx3y8tklqa6JiMkppS+/F5IkSZIkSVINy+d2H3sCb1c6XggM2NiclNKaiPgX8BVg6eYu/tlnn3HhhRdSUlJCp06d6NSpE8XFxRx33HEVx9tvv30NfRRJkiRJkiRJUj7ks6SOasaqrpDekjlExDnAOVXH161bx8svvzz55Zdf3raE+dWeLSjb6wBz1qz6kLM+ZARz1rTuWQeQJEmSJEmqTj5L6oVAh0rHhcC7G5mzMLfdxw7AB1UvlFK6A7gDICJK68M/WTdnzTJnzakPGcGcNS0iSrPOIEmSJEmSVJ0mebz2JKAoIrpERAtgBDCmypwxwOm5598Ant7UftSSJEmSJEmSpIYlbyupc3tMnwc8ATQF7kopzYiI/wJKU0pjgDuBP0fEHMpXUI/IVx5JkiRJkiRJUt2Tz+0+SCk9BjxWZeynlZ6vBE7aysveUQPRaoM5a5Y5a059yAjmrGn1JackSZIkSWpkwt01JEmSJEmSJElZyeee1JIkSZIkSZIkbVKdKqkj4uiIeCMi5kTEZdWcL4iIB3LnX4qIzpXOXZ4bfyMijso458URMTMipkXEUxHRqdK5tRHxau5R9UaStZnxjIh4v1KWsyqdOz0i3sw9Tq/62lrO+etKGWdHxEeVztXKd5l7r7siYklETN/I+YiI3+Y+x7SI6FfpXK18n1uQ8dRctmkRMT4ielc6Nz8iXst9l6X5yriFOYdExL8q/dn+tNK5Tf6+1HLOSyplnJ77fdwpd642v88OEfFMRLweETMi4sJq5mT++ylJkiRJkrQxdWa7j4hoCswGjgQWApOAk1NKMyvNORfolVL6XkSMAIanlL4VET2B+4ADgD2AcUC3lNLajHIeBryUUvosIr4PDEkpfSt37tOU0nY1nWsbMp4BlKSUzqvy2p2AUqAESMBkoH9K6cMsclaZfz7QN6X0ndxx3r/LSu99CPAp8KeU0r7VnD8WOB84FhgA/CalNKCWv8/NZTwIeD2l9GFEHANcmVIakDs3n/Lfh6U1nWsbcg4BfphSOr7K+Fb9vuQ7Z5W5JwAXpZSG5o7nU3vf5+7A7imlKRHRlvLfsROr/H3P/PdTkiRJkiRpY+rSSuoDgDkppbkppdXA/cCwKnOGAffknv+N/7+9uwu146riAP5fmFixgRoMftA2EqGICEJRRKmgqFQRTFvsQwTrBwVRFPFZQaG+CIKPYvHjRbRFa6sRWm0giqBGrUWJRZFqRUMKhQRaYyUQXT7MRMZrbnJS7505Mb/fy53ZZ/Y966y7uQ/r7FmTvLmqahy/u7tPd/djSR4df98icXb3D7r76fH0SJJrtimWZxzjebw1yaHuPjkWqg4leduaxPmuDF9GzK67f5Tk5HkuuSlDMbO7+0iS543Fw9nyeaEYu/snk+LjEuvybBwXyuVm/pd1fdEuMs4l1+bj3f3wePzXJL9NcvWGyxZfnwAAAACbWaci9dVJ/jI5P5b/LrT8+5ruPpPkySTPX3HunHFO3Z7kgcn5c6rqoao6UlU3b0eAWT3Gd463/t9TVdde5NytsPJ71dAyZV+Sw5PhOXK5qs0+y5z5vBgb12UnebCqfllVH1gopqnXVdWvq+qBqnrFOLaWuayq52Yo7H5rMrxIPmtogXR9kp9teOlSW58AAADAZWTH0gFM1DnGNvYi2eyaVeZulZXfq6reneE2+jdMhvd29/GqemmSw1V1tLv/sECM301yV3efrqoPZtih/qYV526Vi3mvA0nu2dDCZY5crmod1uZKxnY0tyd5/WT4hjGXL0hyqKp+N+4kXsLDSV7S3afGNhXfTnJd1jCXo3ck+XF3T3ddz57PqtqVoVD+se5+auPL55iylusTAAAAuPys007qY0munZxfk+T4ZtdU1Y4kV2W4HX+VuXPGmap6S5JPJNnf3afPjnf38fHnH5P8MMOux9lj7O4Tk7i+mORVq86dM86JA9nQTmGmXK5qs88yZz4vqKpemeRLSW7q7hNnxye5fCLJfdm+djkX1N1Pdfep8fj+JDurak/WLJcT51ubs+SzqnZmKFB/rbvvPccll8T6BAAAAC5P61Sk/kWS66pqX1U9O0Ph5+CGaw4mee94fGuSwz08+fFgkgNVdUVV7cuw6/LnS8VZVdcnuTNDgfqJyfjuqrpiPN6T5IYk2/HQt1VifPHkdH+GPrZJ8v0kN46x7k5y4zi2HVb5m6eqXpZkd5KfTsbmyuWqDiZ5Tw1em+TJ7n488+bzvKpqb5J7k9zW3b+fjF85PnAvVXXlGONvlohxjOFFY6/5VNVrMvyfOpEV18ucquqqDHdKfGcyNms+x1x9OcNDMT+3yWVrvz4BAACAy9fatPvo7jNV9ZEMBZJnJflKdz9SVXckeai7D2YoxHy1qh7NsIP6wDj3kar6RoYi5ZkkH97QFmLuOD+bZFeSb461tj939/4kL09yZ1X9M0Ph7TPdveWF1RVj/GhV7c+Qr5NJ3jfOPVlVn85QEEySOza0MZg7zmR4KN3d4xcSZ82Sy7Oq6q4kb0yyp6qOJflUkp3j5/hCkvuTvD3DQzufTvL+8bXZ8rlCjJ/M0MP98+O6PNPdr07ywiT3jWM7kny9u7+3HTGuGOetST5UVWeS/D3JgfFvf871smCcSXJLkge7+2+TqbPmM8MXNLclOVpVvxrHPp5k7yTWxdcnAAAAwGbqP+t+AAAAAAAwn3Vq9wEAAAAAwGVGkRoAAAAAgMUoUgMAAAAAsBhFagAAAAAAFqNIDQAAAADAYnYsHQD8P6iqfyQ5Ohm6ubv/tFA4AAAAAHDJqO5eOga45FXVqe7etXQcAAAAAHCp0e4DAAAAAIDF2EkNW2BDu4/HuvuWJeMBAAAAgEuFIjVsAe0+AAAAAOCZ0e4DAAAAAIDFKFIDAAAAALAYRWoAAAAAABajJzUAAAAAAIuxkxoAAAAAgMUoUgMAAAAAsBhFagAAAAAAFqNIDQAAAADAYhSpAQAAAABYjCI1AAAAAACLUaQGAAAAAGAxitQAAAAAACzmX2+7J9ntIWLlAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Prepare the system to start close to a steady state\n",
    "ode.set(pars = {'F': 0} )       # Lower bound of the control parameter \n",
    "ode.set(ics =  {'y': 0} )       # Close to one of the steady states present for this parameter\n",
    "\n",
    "PC = dst.ContClass(ode)         # Set up continuation class\n",
    "\n",
    "PCargs = dst.args(name='EQ1', type='EP-C')     # 'EP-C' stands for Equilibrium Point Curve. The branch will be labeled 'EQ1'.\n",
    "PCargs.freepars     = ['F']                    # control parameter(s) (it should be among those specified in DSargs.pars)\n",
    "PCargs.MaxNumPoints = 100                      # The following 3 parameters are set after trial-and-error\n",
    "PCargs.MaxStepSize  = 0.1\n",
    "PCargs.MinStepSize  = 1e-5\n",
    "PCargs.StepSize     = 2e-2\n",
    "PCargs.LocBifPoints = 'LP'                     # detect limit points / saddle-node bifurcations\n",
    "PCargs.SaveEigen    = True                     # to tell unstable from stable branches\n",
    "PC.newCurve(PCargs)\n",
    "PC['EQ1'].forward()\n",
    "PC.display(['F','y'], stability=True, figure=3)  # stable and unstable branches as solid and dashed curves, resp.\n",
    "plt.xlim(0,2)\n",
    "plt.ylim(0,1.5)\n",
    "plt.title(\"Bifurcation diagram AMOC\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}