ch2-2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 第二章 建模仿真原理：奶茶问题（Part 2）\n",
    "\n",
    "### 奶茶问题（Part 2）\n",
    "\n",
    "##### 问题描述：\n",
    "假设我在家泡了杯热茶，还带了杯从冰箱里拿出来的牛奶，开车去上班，我希望上班的时候能喝上热奶茶。那我是在家出发的时候把牛奶倒进热茶里，还是在开车上班的途中某个时间点把牛奶倒进热茶里呢？\n",
    "\n",
    "##### 一些假设：\n",
    "1.刚泡的热茶温度为90摄氏度，有300mL。\n",
    "\n",
    "2.牛奶从冰箱里拿出来的温度为5摄氏度，有50mL。\n",
    "\n",
    "3.车里的环境温度为22摄氏度。\n",
    "\n",
    "4.茶泡在保温杯里，从家到单位路程为30分钟，茶的温度降到了70摄氏度。\n",
    "\n",
    "5.牛奶瓶不保温，15分钟后就会升温到20摄氏度，接近车内环境温度。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "引用上一个笔记本的代码"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_func(state, t, system):\n",
    "    \"\"\"Update the thermal transfer model.\n",
    "    \n",
    "    state: State (temp)\n",
    "    t: time\n",
    "    system: System object\n",
    "    \n",
    "    returns: State (temp)\n",
    "    \"\"\"\n",
    "    r, T_env, dt = system.r, system.T_env, system.dt\n",
    "    \n",
    "    T = state.T\n",
    "    T += -r * (T - T_env) * dt\n",
    "    \n",
    "    return State(T=T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_simulation(system, update_func):\n",
    "    \"\"\"Runs a simulation of the system.\n",
    "    \n",
    "    Add a TimeFrame to the System: results\n",
    "    \n",
    "    system: System object\n",
    "    update_func: function that updates state\n",
    "    \"\"\"\n",
    "    init = system.init\n",
    "    t_0, t_end, dt = system.t_0, system.t_end, system.dt\n",
    "    \n",
    "    frame = TimeFrame(columns=init.index)\n",
    "    frame.row[t_0] = init\n",
    "    ts = linrange(t_0, t_end, dt)\n",
    "    \n",
    "    for t in ts:\n",
    "        frame.row[t+dt] = update_func(frame.row[t], t, system)\n",
    "    \n",
    "    return frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_system(T_init, r, volume, t_end):\n",
    "    \"\"\"Makes a System object with the given parameters.\n",
    "\n",
    "    T_init: initial temperature in degC\n",
    "    r: heat transfer rate, in 1/min\n",
    "    volume: volume of liquid in mL\n",
    "    t_end: end time of simulation\n",
    "    \n",
    "    returns: System object\n",
    "    \"\"\"\n",
    "    init = State(T=T_init)\n",
    "                   \n",
    "    return System(init=init,\n",
    "                  r=r, \n",
    "                  volume=volume,\n",
    "                  temp=T_init,\n",
    "                  t_0=0, \n",
    "                  t_end=t_end, \n",
    "                  dt=1,\n",
    "                  T_env=22)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 热茶和牛奶的最佳混合时间\n",
    "前面我们对热茶的冷却过程进行了模拟，一般情况下，我们并不知道参数r的取值，但我们能通过测量估计出它的值。给定一个初始温度，一个最终温度，以及间隔的时间，我们能通过尝试，找出r值。\n",
    "下面我们使用一种更好的方法来找出r值，我们采用二分法搜索。\n",
    "我们使用的ModSim库当中有个函数root_bisect就是用来搜索非线性方程的根。举一个简单的例子，假设我们要求解如下多项式的根：\n",
    "$$f(x)=(x-1)(x-2)(x-3)$$\n",
    "我们可以写一个函数："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func(x):\n",
    "    return (x-1) * (x-2) * (x-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "root_bisect函数的第一个参数是待求解的函数，第二个参数是根所在的区域。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>converged</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>root</th>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "converged    True\n",
       "root            1\n",
       "dtype: object"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res = root_bisect(func, [0.5, 1.5])\n",
    "res.root"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "如果我们给定另外一个求解区域，会得到另一个解。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>converged</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>root</th>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "converged    True\n",
       "root            2\n",
       "dtype: object"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res = root_bisect(func, [1.5, 2.5])\n",
    "res.root"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们想要找到最终温度为70摄氏度时对应的$r$，所以我们定义了一个误差函数，以$r$作为输入参数，输出是70摄氏度和最终温度的差值。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def error_func1(r):\n",
    "    \"\"\"Runs a simulation and returns the `error`.\n",
    "    \n",
    "    r: heat transfer rate, in 1/min\n",
    "    \n",
    "    returns: difference between final temp and 70 C\n",
    "    \"\"\"\n",
    "    system = make_system(T_init=90, r=r, volume=300, t_end=30)\n",
    "    results = run_simulation(system, update_func)\n",
    "    T_final = get_last_value(results.T)\n",
    "    return T_final - 70"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "选取$r=0.01$，我们发现茶有点热。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.2996253904030937"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "error_func1(r=0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "选取$r=0.02$，我们发现茶有点凉。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-10.907066281994297"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "error_func1(r=0.02)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "使用root_bisect求解70摄氏度对应的$r$值。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>converged</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>root</th>\n",
       "      <td>0.0115431</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "converged         True\n",
       "root         0.0115431\n",
       "dtype: object"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res = root_bisect(error_func1, [0.01, 0.02])\n",
    "r_tea = res.root"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "如果用估计的$r$值来运行仿真，可以看到最终温度为70度。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "69.99999985860761"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tea = make_system(T_init=90, r=r_tea, volume=300, t_end=30)\n",
    "results = run_simulation(tea, update_func)\n",
    "T_final = get_last_value(results.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 练习\n",
    "重复以上的过程，估计$r_{milk}$值。牛奶温度初始为5摄氏度，在15分钟之后达到20摄氏度。\n",
    "在你使用root_bisect函数之前，可以先用试凑的方法，看看能够估计得多接近。以下程序帮助你开始这个过程："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "18.499850754390966"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r_milk = 0.1\n",
    "milk = make_system(T_init=5, r=r_milk, volume=50, t_end=15)\n",
    "results = run_simulation(milk, update_func)\n",
    "T_final = get_last_value(results.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 混合液体\n",
    "下面这个函数表示将两种液体混合的过程，计算混合液体的温度。输入两种待混合液体的System对象，输出为混合后液体的System对象。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mix(system1, system2):\n",
    "    \"\"\"Simulates the mixture of two liquids.\n",
    "    \n",
    "    system1: System representing coffee\n",
    "    system2: System representing milk\n",
    "    \n",
    "    returns: System representing the mixture\n",
    "    \"\"\"\n",
    "    assert system1.t_end == system2.t_end\n",
    "    \n",
    "    V1, V2 = system1.volume, system2.volume\n",
    "    T1, T2 = system1.temp, system2.temp\n",
    "    \n",
    "    V_mix = V1 + V2\n",
    "    T_mix = (V1 * T1 + V2 * T2) / V_mix\n",
    "    \n",
    "    return make_system(T_init=T_mix,\n",
    "                       r=system1.r,\n",
    "                       volume=V_mix,\n",
    "                       t_end=30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "mix函数需要System对象有temp这个系统变量。make_system函数初始化这个变量；随后的函数保证当我们运行仿真的时候能够不断更新它。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_and_set(system):\n",
    "    \"\"\"Run a simulation and set the final temperature.\n",
    "    \n",
    "    system: System\n",
    "    \n",
    "    returns: TimeFrame\n",
    "    \"\"\"\n",
    "    results = run_simulation(system, update_func)\n",
    "    system.temp = get_last_value(results.T)\n",
    "    return results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 立即混合\n",
    "下一步我们仿真在上班前在热茶中立即倒入牛奶去上班的这种情况。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "90"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tea = make_system(T_init=90, r=r_tea, volume=300, t_end=30)\n",
    "tea.temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "milk = make_system(T_init=5, r=r_milk, volume=50, t_end=30)\n",
    "milk.temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "77.85714285714286"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mix_first = mix(coffee, milk)\n",
    "mix_first.temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "61.4285713124277"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mix_results = run_and_set(mix_first)\n",
    "mix_first.temp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 最后混合\n",
    "然后我们会看在到达上班地点后再混合奶茶的情况。我们将分别对热茶和牛奶进行仿真。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "69.99999985860761"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tea_results = run_and_set(tea)\n",
    "tea.temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "21.279350309321327"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "milk_results = run_and_set(milk)\n",
    "milk.temp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下图是看起来的冷却结果。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving figure to file figs/chap16-fig01.pdf\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAEPCAYAAACdhMnXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXwTdf4/8NfkTnrfFFqOFguVG2q5lxtcRQ4RkQX8ieCCuIIIK6gooKKIVlcOd4VFVFC5+9UVFRVvRQRBoCAtR7H3faVN0lzz+yPt0NAWArQJbV/PxyOPJDOTyXus7YvPZz7zGUEURRFERERuJPN0AURE1PIwfIiIyO0YPkRE5HYMHyIicjuGDxERuZ3C0wXcKJPJhKSkJISEhEAul3u6HCIiAmCz2ZCfn4+uXbtCo9HUWt/kwycpKQnTpk3zdBlERFSH999/H3FxcbWWN/nwCQkJAeA4wFatWnm4GiIiAoCcnBxMmzZN+ht9uSYfPtVdba1atUJERISHqyEioprqOx3CAQdEROR2DB8iInI7hg8REbkdw4eIiNyO4UNERG7H8AGQW2RAib7S02UQEbUYTX6o9Y0682cR/rn2BwCAr5cKbVv5oG2YD9q28pVe+3mrPVwlEVHz0uLDJyxAh563hCA5rRhlFWYknS9E0vlCp238vdU1QskRTJFhPvD1UnmoaiKipq3Fh0+ArwbPzx0AURRRUGJCWm4Z0nL0jkduGdJz9Sgpr0TJuUqcOFfg9Fl/H7UjkMJ8EFkVTpFsKRERXVWLD59qgiAgJECLkAAt+nQOk5bb7SIKSoxIy9UjLacMf+bokZ7reJToK1GiryOUvNWIrGolRYbVDCUVBEFw96EREd10GD5XIZMJCA3UITRQh7jY+kLpUitJaimVV+LkeedQ8tGppECKDPNGZKgjoAJ9NQwlImpRGD7Xqb5QEkUR+SVGpFeFUnUgpeXqoTeYcepCIU5dcD6npNMoHIEUWiOYwnwQGqCDTMZQIqLmh+HTwARBQGiADqEBOqfuO1EUUVRmkoIoPbe8KqDKoDdYkPxnMZL/LHbal0opR0Soo4UU2cpbCqfwYC8o5BwlT0RNF8PHTQRBQJCfFkF+WvSMCZWWi6KI0nJzVbddOTKqwikjT4+iskpcyCzFhcxSp30p5ALCg70QUd1SCvVGRJgPIkK9oVHxR0pENz/+pfIwQRDg76OGv08Iund0vu9FudEihVF6jUdesbGq5VSOgyeznT4TGqBFhNSF5y0FFIeFE9HNhOFzE/PWKtG5fSA6tw90Wm6qtCIj39FKSs9zdN9l5OmRlV+BvGIj8oqNOHomz+kzft4qRIQ6WkcRNYIpxF/L80pE5HYMnyZIo1agY4Q/Okb4Oy232uzILqhARl7VOaU8PTJy9cjIK0dpuRml5bUHO6hVcrQJ8XacWwq7FE6tg72gUtZ9EygiohvF8GlGFHJZ1Wg5H/Tvdml59QW0GXn6qkAqR0aeI5xK9HWfVxIEICxQV6O15C295kW0RHSjGD4tQM0LaHt1CnVaV24wO4IoV4/MfEcoZeTpkV1oQE7V48gfuU6f8dGpagSSt6PlFOaDsEAdR+ERkUsYPi2ct05V53kli9WOnEJHF54jkMql7jy9wYw/Lhbhj4tFTp+Ryxyj8Kq78apbSm1CveGj44AHIrqE4UN1UioudeHVVH29UnUrKbNGMOWXGKWgOnTKeX++XqpLraSq5zah3mgVxGuWiFoihg9dk5rXK10+NNxktjoGPOSWO0bjVbWasvLLUVZhxunUIpxOrd1aahXk5RRIbUIcD86FR9R8MXyowWhUCnRo7YcOrf2clouiiMJSk6OVVBVKmXnlyMwvR36JEZn5jteX89IqEVEVSK1DvBAR4oM2od4ID/aCmiPxiJo0hg81OkEQEOyvRbC/Fj1inFtLlRYbsqrCJ6MqkKqDqcJoQXJaMZLTii/bHxDir5VaSK1rtJiC/bWQ87olopsew4c8Sq2U19taKimvlILI0X1Xgcx8PXIKDdLFtMdS8p0+p1TIpEEPjoeXI5xCvOHrxW48opsFw4duSoIgIMBHgwAfDbpGBzuts9ocI/GyCiqkcHIEUzmKykzSzQAv56VVSmHUOti7xmsv6DRKdx0aEcHN4fP777/jhRdeQGpqKgIDA/H3v/8dkydPhtlsxvPPP4/9+/dDJpNh5syZmDNnjjtLoyZEIZdVDeP2AW51XmcwWZBVUFHVlVcdTnpkFVSgwmhBSloJUtJKau0z0FeN8OCqbrzgqlAK8UJ4EGd6IGoMbgsfu92OefPmYcmSJRg/fjxOnDiBadOmoVu3bti3bx9SU1Px5ZdfQq/XY/bs2QgLC8OECRPcVR41EzqNss6ph6pnD3e0kqpaS1UhlVVQgaKyShSVVdaafqj6/FLrYG+Eh3g5hRMvqiW6fm4Ln9LSUhQWFkIURYiiCEEQoFAooFQqkZiYiNWrV8PPzw9+fn6YNWsWtm/fzvChBnNp9nA1ukQFOa2rvivt5YGUmV+O3KJL55d+P+t8fkkmExAWqEPrYC+EB3uhdbCjtdQ62BuhAVrIGUxE9XJb+AQEBGD69OlYunQpnnrqKdhsNjz99NMICQlBfn4+OnbsKG3boUMHpKSkuKs0auFq3pW2VyfndVabHblFBqdAys6vQGZBOQpKjMguqEB2QUWtfSrkjmAKrxFI1SEVEqDjiDxq8dza7aZSqZCQkIDRo0fj2LFjePTRR+Hr6wsA0Gg00rZarRYmk8ldpRHVSyGXSSPnLme22JBdWIGs/ApkF1S3miqQVVDuuK4pvwKZ+RXAH5fvU0BYoJfjnFJVi8nxzGCilsNt4fPFF1/g2LFjWLJkCQAgPj4ekyZNQmJiIgCgsrJS2tZoNEKn07mrNKLrolLK0a6VL9q18q21zlRpdQRTVTdedoHjdXZBhTQ9UV0X1lYHU3UYhdd4hAbwHBM1H24Ln5ycHJjNZucvVygQGBiIkJAQXLhwAWFhYQCA1NRUp244oqZGo657tgfAOZiyq8OpqgV1pWCSyQSEBeicAik82DEir1WQDkoFR+VR0+G28Bk4cCASEhKwY8cO3HvvvTh16hR27tyJF154AeHh4diwYQM6deoEg8GAzZs34/7773dXaURudbVgyikyILvAubWUXVjhOMdU6HiNZOfPCQIQ7K9FeNClQKoZTho1L+mjm4vb/o+85ZZbsH79erzxxhtYs2YNgoODsWjRIowcORKDBw/G6tWrMXbsWNjtdkyZMgVTp051V2lENw2NWoH24b5oH167K89ssSG3yCCFUlZBOXKqgimvyID8YiPyi404ca6g1mcDfNRoVSOQWgU5uvVaBXnBR6fkzA/kdoIoiqKni7gRGRkZGDFiBA4cOICIiAhPl0PkEVabHXnFBmn0XXVrKbugAjmFBlht9no/66VRSIEkPQc5noP8NJBxAARdh6v9bWZbnKgZUMhlVcO5a4/Ks9lFFFZ12eUUXgqmnAIDsgsrUGGy4lxGKc5llNb6rFIhQ1igrkYw6aRwCgvUcfYHum4MH6JmTl7jOqYetzjPKl4980NOYXUgVSCr6jmn0ICS8krpBoGXEwQgyFeDsKBLgx5a1XjmRK50JQwfohas5swPl99KHQCMlVbkSC0mw6WQKqxAXrERBaUmFJSaak1LBABataJGIHkhPEiHsKpwCvHXQangsPGWjOFDRPXSXmFkns1mR36JsSqQDMiVgskRUgaTFalZZUjNKqv1WVnV6LxWVd134cFeaBXohbCqsOIgiOaP4UNE10Uul0mtmp6XrRNFEXqD5VKrqbACuYUGRzAVOYaNV8+ZVxedRiGdawoL1KFV4KVWU2gAzzU1BwwfImpwgiDA10sFXy8VYtoG1FpvsdqRX3wpjKpbS7mFjkEQV2o1AUCQn+ayYNIhLNARTgE+HKHXFDB8iMjtlApZ1T2Tao/Oq9lqyi0y1HrOKzaisNSEwnrONSkVMoQGaBFWM5yqu/QCdfDWqdxxiHQVDB8iuqlcrdVUfa4pt6rVlFtkcHrtuG9T1aSudfDSKBAaqENYdShJLScdwgJ0nA3CTfhfmYialJrnmnogpNZ6Y6W1KpCqgqnI0b2XWxVOFVfp0vPzVknBVLMFFRaoQ2iAlnPoNRCGDxE1K9orTFEkiiLKKsxSKEmPGl16peVmlJab67zdOgAE+mqkMAqrun6q+nWwv5Yzj7uI4UNELYYgCPDzVsPPW11nl57dLqJYb6pqKRmQV+zo0ssrNiCnyICCEiOKykwoKjPhj4tFtT4vE4BAP63USgqtOucUGugYpcdwuoThQ0RURSYTEOSnRZCfttbt1gHH+aaCUhPyihzdeDlVLaf8YiNyiwwoLDWioMTxOFXX/gUgyF+L0IBL4RQWcKn1FOSnbTEX3zJ8iIhcJJfLpC62bgiutd5itaOw1BFEedVdesWO13lFBhSWmaTZx+sKp+opi0JrtJZCA6rCKlCHEH9ts7nGieFDRNRAlIpLgyHqYrHaHRfYVnfpVQdTVcupqPTSlEWnU2t36wGO22NUh1FogBYhAY4wDKkKKW0TGa3XNKokImoGlAqZdE+lulhtVeF0WSjlFxuRW+w451Ssr0SxvhLJacV17sNHp0RIjRaT82vtTTPhK8OHiOgmoZBfueVks9mlrjtHKF0655RXbEB+iRF6gwV6QykuZNa+RQYAaFRyhFS1mEJqnH8KqQqnIF8N5G4YFMHwISJqIuRyWVVY6OocEGG3iygtr3S0nIqNyK96ziu+FFAGkxXpueVIz619mwygetCFBqEBOtw9tCPiu7RqlGNh+BARNRMymYAAXw0CfDXo1K7ubcqNFuTXCKPLn4v1ldKgiAAftWfD57vvvsMPP/yApKQkFBYWQiaTISQkBN26dcPQoUPRt2/fRimOiIgalrdWCW9t3bfJAACL1Yb8EiOKSk2IalP3Ng3hiuHzySefYN26dSgtLcWAAQMwePBg+Pv7w263o7i4GGfOnMFjjz2GwMBAzJ07F3fddVejFUpERI1PqZDXe0v2hlRv+MyaNQtarRYrV65EfHw8ZLK6T0DZ7XZ8//33+PDDD7F3715s2bKl0YolIqLmod7wmT9/Pnr06HHVHchkMgwdOhRDhw7F77//3qDFERFR81TveDpXgudyPXtefj9DIiKi2q44mDstLQ2LFy9GZmam0/KlS5di4cKFtZYTERG5ot7wuXjxIqZMmYLk5GQYjc73We/bty/Onz+Pe++9F+np6Y1eJBERNS/1hs8bb7yBnj17IjExER07dnRaN3HiROzevRsxMTFYu3ZtoxdJRETNS73hc/jwYTzyyCNQKOoek6BSqfDII4/g0KFDjVYcERE1T/WGT0VFBfz8rnyBUVhYGPR6fYMXRUREzVu94RMdHX3VodO///47IiIiGrwoIiJq3uoNn8mTJ+P1119HRkZGnevT09ORkJCA8ePHN1pxRETUPNV7kemUKVPw66+/4s4778SYMWPQrVs3+Pj4oLS0FCdPnsSXX36Jv/zlL5g5c6Y76yUiombginO7JSQk4H//+x/27NmD9evXQ6/Xw9/fH927d8crr7yC0aNHu6tOIiJqRq46q/Vdd93FCUOJiKhB1XvO59///jdMJpPLOzIYDFi/fn2DFEVERM1bveGjUqkwduxYvPbaazh58mS9Ozh9+jRWr16NO+64A2q1ulGKJCKi5uWKt1S4/fbbsWnTJtx///1Qq9Xo2LEjAgICIIoiioqKcPbsWVgsFowfPx7btm3jsGsiInLJFc/5tGnTBitWrMDixYvxyy+/4NSpUygsLIQgCOjcuTMeeugh9O/fHxqNxqUvy8vLw4oVK3Do0CGo1Wrce++9eOyxx2A2m/H8889j//79kMlkmDlzJubMmdMgB0hERDcfl26j7e3tjZEjR2LkyJE39GXz5s1Dly5d8PPPPyMvLw8zZsxAdHQ0UlJSkJqaii+//BJ6vR6zZ89GWFgYJkyYcEPfR0REN6cr3lKhIR0/fhzp6elYtmwZ1Go1IiMjsXXrVvTt2xeJiYmYO3cu/Pz8EBERgVmzZmH79u3uKo2IiNzMbeGTlJSEmJgYrF+/HoMHD8bIkSPx5ZdfQqPRID8/32nm7A4dOiAlJcVdpRERkZu51O3WEEpLS/Hbb78hPj4eBw4cwIULFzB79mwEBgYCgNN5I61We03DvImIqGlxW/ioVCp4e3vj0UcfBQB07twZ99xzDxITEwEAlZWV0rZGoxE6nc5dpRERkZu53O128uRJLFq0CHfffTdycnLw9ttv48svv3T5i6KiomA0GmE2m6VlNpsNfn5+CAkJwYULF6TlqamptW5gR0REzYdL4fPDDz9gxowZUKvVOHfuHCwWC4xGIxYuXIi9e/e69EUDBw5EYGAgXn75ZZjNZiQnJ2P37t248847MW7cOGzYsAFFRUXIyMjA5s2bMW7cuBs6MCIiunm51O32r3/9C0uWLMHUqVPx2WefAQAeeeQRBAYGYtOmTbj77ruvug+1Wo1t27bh+eefx+DBg6FSqTB79myMGTMGQ4cOxerVqzF27FjY7XZMmTIFU6dOvbEjIyKim5ZL4XP+/HkMGjSo1vJBgwbhpZdecvnLIiMjsXHjxlrL1Wo1li9fjuXLl7u8LyIiarpc6nYLCwvD6dOnay3/5ZdfEB4e3uBFERFR8+ZSy+fvf/87nn32WaSnp8Nut+P7779HZmYmPvjgAzz55JONXSMRETUzLoXPpEmTEBQUhLfeegtqtRoJCQmIiorC6tWrcfvttzd2jURE1My4FD6bNm3CuHHj8OGHHzZ2PURE1AK4dM7nP//5j9NFoERERDfCpfAZOnQotm7ditLS0sauh4iIWgCXut3+/PNP7Nu3D9u2bYO3t3etO5b++OOPjVIcERE1Ty6Fz5QpUzBlypTGroWIiFoIl8Jn8uTJjV0HERG1IC6Fz6JFi664PiEhoUGKoevz0ksvYefOnQgNDcXWrVsxZ84cXLx4ETNnzsT8+fM9XR4RUS0uhY9M5jwuwWazIT09HSkpKZgxY0ajFEau27ZtGzZv3ox+/frh448/hl6vx+HDh6FQuO2OGURE18Slv06vvPJKncv/85//IC0trUELaumOHTuG1atXIyUlBaGhoVi0aBGGDx+ON998E3v27IHJZEKvXr2wbNkyhIeHIy4uDlarFXPmzEFQUBDy8/NhtVpx2223Ye/evdBqtXj++edx5MgR+Pj4YPr06XjggQcAOP4RsWnTJuzatQvl5eXo168fli9fLt3gj4iosdzQP43Hjh2LcePG4cUXX2yoehrdyv/+giN/5Lrt++Jiw7B8dj+Xti0qKsJDDz2EhQsXYsqUKTh8+DDmzJmDO++8EydPnsS2bdsQGhqKV155BXPnzkViYiKOHTuGTp06Yfv27YiNjcXevXvx7rvv4qOPPoLNZsOkSZMQFxeH77//HtnZ2ZgzZw78/f0xYcIEvPfee/j444+xZcsWhIaG4uWXX8bChQvx7rvvNvJ/FSJq6Vy+mVxdEhMT4eXl1VC1tHjffPMNwsLCMG3aNCgUCvTv3x8ffPABvvrqK8ybNw+RkZFQq9V44oknkJWVhRMnTlxxf0lJSUhLS8PSpUuhVqvRvn17zJw5E9u3bwcA7Ny5E//4xz/Qtm1baDQa/POf/8Thw4dx8eJFNxwtEbVkLrV86rqdgsFggMFgwNKlSxu8qMbkaivEEwoLC2vNEt61a1cYjUa0bt1aWqZSqRAaGoqcnJwr7i8zMxNGoxH9+l06ZrvdDn9/fwBAVlYWnn76aTz77LPSeoVCgczMTLRv374BjoiIqG4uhc+CBQuc3guCAKVSiW7duiEqKqpRCmuJQkNDkZvr3CX4zjvvQBRFZGZmomfPngAAs9mM3NxcBAUFXXV/QUFBThcBFxUVwWQySeufffZZDB48WFqfkpLC4CGiRudSt1thYSHuuusuTJ48GZMnT8Y999yD8ePHIzQ0FGvWrGnsGluMIUOGIDc3F7t27YLNZsPBgwexdu1azJs3D//+97+Rnp6OyspKrFmzBgEBAejdu/cV99e9e3d4e3vjzTffhNlsRlFREebNm4e1a9cCACZMmIANGzYgOzsbNpsNGzduxLRp06RwIiJqLPW2fPLy8lBRUQEAeOONN9C3b1+pu6baH3/8gW3btuGJJ55o3CpbiICAAGzcuBEvvfQSVq9ejbCwMCQkJGDQoEGwWCy4//77UVpait69e2PLli1QqVRX3J9KpcLGjRvx4osvYvDgwRAEASNHjsRTTz0FwHGfJqvVimnTpqGkpAQxMTHYvHkzfH193XG4RNSCCaIoinWt2L9/PxYsWABBEAAANTcTBEF6P3HixGu6lXZDy8jIwIgRI3DgwAFERER4rA4iIrrkan+b6235jBkzBl9++SXsdjvGjBmDHTt2ICAgQFovCAJ0Ot1VzzsQERFd7ooDDiIjIwEAp06dglwur3Mbi8UCpVLZ8JUREVGz5dJot9LSUmzcuBHnzp2D3W4H4OiGM5vNOHv2LH799ddGLZKIiJoXl0a7PfPMM/jss88QERGBQ4cOoV27dpDJZPjtt98wd+7cxq6RiIiaGZdaPr/88gvWr1+P/v3749ixY5g0aRK6du2KV199FceOHWvsGomIqJlxqeVjNpulCw+jo6Nx6tQpAMCkSZNw9OjRRiuOiIiaJ5fCp3379lLIREVF4fjx4wAAk8kEo9HYeNUREVGz5FK326xZs7B06VLYbDbccccdGD9+PARBwPHjxxEXF9fYNRIRUTPjUstnwoQJ2LJlC6KjoxEVFYU333wTpaWl6NmzJ1atWtXYNdJlsrKy0KtXL+j1emRkZKBTp04oKytzek1EdDNzqeXz8MMPY/HixYiOjgYADB482GkySnKv1q1bSwM9SktLPVwNEdG1c6nlc/To0avOI0Y3LiMjA3Fxcdi6dSsGDhyI+Ph4bN26Fe+//z4GDx6Mvn374p133nGphWOz2bBw4UJMmzZNmqOPiOhm4VLL54EHHsCSJUswc+ZM6YZmNXXo0KFRimuJ9Ho9kpKS8M033+Crr77C4sWLMX78eBw4cAA//PADHn30UcTHx19xH3a7HU899RTy8/OxadMm6HQ6N1VPROQal8LnjTfeAABpxFvNyUYFQcAff/zRSOU1jgurJtW7Lvivc+DbezQAoOzoFyj47K16t416eo/0OmPzP2HOuXDFbVw1b948qFQq9O/fHzabDffffz9UKhWGDRsGm80GmezKDdaVK1fi6NGj+Oyzzxg8RHRTcil8vvjii8aug2qovnVF9Xx6Pj4+ACCFTj0TkUuys7Oh1+tx+PBhDBkypBErJSK6Pi6FT9u2bQEAqampSE1NRf/+/VFcXOx0a+emxNXWiG/v0VIr6GoiZr1yIyU5qW5ZXq8NGzZg3759ePbZZ7Fv3z54e3s3UGVERA3DpQEH5eXlePjhh/HXv/4V//jHP1BQUIBVq1ZhwoQJtW77TJ6nVCoxffp0hISE8E6zRHRTcil81qxZg/Lycnz99dfSYIOlS5dCp9PhxRdfbNQC6frIZDI8//zz2Lt3L3755RdPl0NE5MSlbrdvvvkG//73v5262SIjI7Fs2TLMnDmz0YpraSIiIpCcnCy99/X1dXoPQHpf/Vxzm8u3j42NRVJSUmOXTUR0zVxq+VRUVECr1da5zmq1XtMXlpWVYejQodi7d6/0fv78+ejTpw8GDx6MPXuufXQYERE1LS6Fz6BBg7BhwwbYbDZpWXFxMdasWYOBAwde0xcuX77c6TzRihUrIJPJ8OOPP2Ljxo1ISEjgzemIiJo5l8Jn2bJlyMjIQL9+/WAymTB79mwMGzYMxcXFeOqpp1z+ssTERJSXlyMmJgYAYDQasX//fixYsABarRaxsbGYPHkyduzYcX1HQ0RETYJL53xCQ0Oxc+dO/PTTTzh37hxsNhuio6MxePDgq17wWC09PR3r16/H9u3bMXv2bADAxYsXIQiCdK8gwHHLhq+//vraj4SIiJoMl8KnWmRkJIxGI+RyOW655RaXg8dms+Gf//wnlixZgpCQEGl5RUUFVCqV03UtWq0WJpPpWsoiIqImxqXwKSsrw2OPPYaDBw9CLpdDFEXY7XaMGjUKL7/8cr2DEaq9+eab6NChA0aPdr5gU6fTwWw2S9P0AI6uOE4JQ0TUvLnUdHn22WdRXFyMXbt24cSJEzh58iR27NiBtLQ0l67z2bdvH/bv34+4uDjExcUhJSUFK1euxNatWyGKIjIyMqRtL1y4gI4dO17/ERER0U3PpfD5/vvv8cILL6Br166QyWSQyWTo3r07nnvuOezfv/+qn//8889x9OhRHDlyBEeOHEFMTAyWL1+Ol156CSNHjkRCQgIqKipw5swZ7Nq1C+PGjbvhAyMiopuXS+ETEBBQ571jrFbrVbvcrua5556DTCbD8OHD8dBDD2Hu3LmcDJOIqJlz6ZzPokWL8Mwzz+Dxxx9Hnz59oFAocOrUKbz00kuYPn06UlNTpW1dubfPRx99JL328/PDa6+9dh2lExFRUyWIV5ufH0Dnzp0vfaDGvXxqLvPUvX0yMjIwYsQIHDhwABEREW79biIiqtvV/jbzfj5EROR213Q/HyIioobgUvicOXMGq1evxrlz52A2m2ut51xsRER0LVwKn6VLl0Kr1WLBggXS/XyIiIiul0vhc/HiRezevZsXfxIRUYNw6Tqf3r174+zZs41dCxERtRAutXyef/55TJ06Fd9//z0iIyNrTSg6d+7cRimOiIiaJ5fC580330R+fj6OHj1a6zoeQRAYPkREdE1cCp9PP/0U69atw8iRIxu7HiIiagFcOufj6+vr0rQ5RERErnApfBYvXowXXngBZ86cgcFggNlsdnoQERFdC5e63V5++WWUlJRg4sSJda5393xuRETUtLkUPq+88kpj10FERC2IS+HTv39/6XV5eTm8vLyk2a2JiIiulUvnfABgy5YtGDRoEOLj45GRkYEnn3wSq1atgtVqbcz6iIioGXIpfN5++2289957WLx4MVQqFQBgxIgR+Oyzz3gjOCIiumYuhc/OnTuxYsUKTJgwQepuG7XW8gQAAByKSURBVDlyJFavXo1PPvmkUQskIqLmx6Xwyc7ORlRUVK3lbdq0QWlpaYMXRUREzZtL4RMbG4vPPvus1vLt27cjNja2wYsiIqLmzaXRbkuWLMHs2bNx5MgRWCwWrFu3DhcuXMD58+exadOmxq6RiIiaGZdaPr169cL+/fsRGxuLIUOGoKysDPHx8fj0008RFxfX2DUSEVEzU2/LZ/369Zg1axa0Wi0AIDg4GAsXLnRbYURE1HzV2/LZsGEDDAaDO2shIqIWot7wEUXRnXUQEVELcsUBBzk5OaisrLzqTlq3bt1gBRERUfN3xfC55557rvhhURQhCAJntSYiomtyxfB577334O/v765aiIiohag3fARBQFRUFIKCgtxZDxERtQAccEBERG5Xb/hMnDgRarXanbUQEVELUW+320svveTOOoiIqAVx+WZyREREDYXhQ0REbsfwISIit2P4EBGR27l0P5+G8tNPPyEhIQEXL15EUFAQZs2ahfvuuw9lZWVYtmwZfvrpJ+h0Ojz22GOYNGmSO0sjImqSRJsFotUK0W6FaLMBdpv0WqZQQuEXIm1nykwB7HaINqvj2W6DaHd8RhPZGQrfYACAKescrKV58I4d0Gh1uy18srOz8eijj+Lll1/GiBEjkJSUhNmzZ6NNmzZITEyETCbDjz/+iIsXL2LWrFmIjIxEfHy8u8ojIqpFtNsgWi1Vf+AvPSt8AiFT6wAAlqJsmAsyINqsEG0WwGZ1bGu3AoIcfnG3S/sr+n4H7BWljm3t1qrPWAGbFV6x/eHTfRgAwJR+BgWfvwXRZqva1ubYb1WoRM75lxQUuXsTYEg5XGf92qieCJ/6DADAbjIge+uz9R5r2KR/Svs0XvgdhrNHmkf4ZGZmYuzYsRg1ahQAoHv37oiPj8fBgwexf/9+fPLJJ9BqtYiNjcXkyZOxY8cOhg8RSRxBYIZoMQNyBeQaLwCAzaBHZfY5iBYzRKsZdmul9Fq0muF7253StqW/fgJTxhnH+upAsZpht5qhDr8FoXc9Iu3zz389CIj2OmsJvXsxvGP7AwDK/ziI4m/fr3M7mcbLKXzKT3wLa2lendsqQyKl13ZrJcx5afX/t7BZL32HUgNBqYEglwMyOQSZQnot9w6UthMUKmgiY6u2kdd4lkGQKyD3ubStunVHCPLGjQe3hU9cXJzTXU9LSkpw5MgRxMXFQRAEtG/fXloXFRWFr7/+2l2lEVEDEUXR8UffbILdbITdbILCNwhyrQ8AoDL3IkzpZyBaTLCbTRAtlbBbTBDNJkAQEDpuvrSv7PdXwFyY6QgKS6WjVVHFr994BI24HwBgzruInO0v1FuT162DpPAxZaag4o+DdW4n1/pKrwW5whE8ggyCQglBroQgV0BQqCAoFJApVNK2yqBw6Dr2AeQK523lSshUGqfvCBg8GXZLZdV6xwNyBQSZAsrAcGk7TesYtJmdAEEmr9pGDkGmhCB3vBdq7Dd0wmP1/0BqkKm1aH1//f+datJF9YQuqqdL214vt57zqabX6/Hwww+jR48e6NKlC1QqFQRBkNZrtVqYTCZPlEbUIol2G+yVBsjUOse/hgFUZp+HpTAL9kqDI0gqDbBXGmE3G6DwDUbgkKkAALvZhPT/zIdYFTaXtxZCxs2HT7chAADTn0ko/HJL3UXIFU7hY60ogU1f5LSJoFQ7HnLlpY95+UEb1aMqGFQQFGrIlCoISsf76uABAL+4O+DVqW/VdsqqhxoyhVLqRgMAQaVBhyd3Sv8trsS7c394d+5/1e0AwKfHcJe2k6m1UIe1d2nbpsrt4ZOamop58+ahY8eOePXVV3H+/HmYzWbp9gwAYDQaodPprrInIqomiiJEsxE2UznsxgrYTeXSa2VAGLTtuwEAzIVZKPrqHdhMFY4wMVXAXlnhaHkAaDM7QfqjV/bb59Afr7sHQtUqWgofQamCTV94aaVcAZlKC5lKA0GlhUx5aZouVVh7+PQeXdVVpHZso9RAplJDUGqc/g60mvIUBEEGQaGWgqTmP1KlfYa0RfjU+s9l1KSJ7OzSdoIgAMLVg4eun1vD5/Dhw5g3bx7uu+8+PP744xAEAe3atYMoisjIyEBkpKPP88KFC+jYsaM7SyO66dgqSmHOT4PNUAa7UQ+bsRw2ox52QxnspgqE3bsUguC4WiJry1JUZp+rcz/e3YdJ4QObBYZzv9WxlQCZWgvRapaWqNvEwG6phEytg0ythUxV9azWQe4dcOmTggxtH33rUojUaJVcTtuuK7Tturp0/Eq/UJe2o6bJbeGTlpaGOXPmYOHChZgxY4a03MvLCyNHjkRCQgJWrVqF9PR07Nq1i3PLUbMhWi0QFJf+IOtPfgtbeQlshlLYDGWwVZTBXvXav/9E+PYZAwAwXPgd+R+vrXe/dpMBcq03AECm0Tn++Gu8INd6QabxhkzjBZnGB9q2sdJnFP6hCJu81LFOrYO86llQa6Ugq+bbaxR8e41y6RirR0kRucpt4fP++++joqICr732Gl577TVp+d/+9jc899xzWLlyJYYPHw6VSoW5c+diyJAh7iqN6JpVn/yu/le+KTMFxtQTsJUXw1ZRAmt5CWwVjodc4422j74lfbbwq3dhN5TVuV9rje4rpX8YNO26QK71gUzrU+PZGzKtDwTlpZPere5b5tL5CZlKC6+Y267rmIkaktvC58knn8STTz5Z7/qagUTkKTXPOditZpQf/wbW8iLY9MVOz3ZDGcIm/RNenfsBcJxIL/7uwzr3aZc5/5r59hwB0WaFXOcLmc4Pcp0v5F5VzzW6szSRndF6+nMu1e1K8BDdTDwy2o3IkyxFWTDnpcOqL4C1rADWskLHc2kBlIGtpD/4AgQUfL6x7p0IMtgrDdJbTWQs/AfcDbm3P+ReVQ9vfyi8/CGonQfPBA6b3mjHRtRUMHyo2RDtNtjKi2EpyYW1JB/WkjxYSvNgLc1H8OgHoQptBwAo/XUfyn77vM59CPJLLQhBoYTvbXdApvaCwicQcp9AKLwDIfcJgFzn69Ta0ETGOi7gIyKXMHyoyRBFEXajHpbiHFiKcyBTauDVyTELhqU0D+lv/gOw2+r8rKUoRwofdeuO0JX2gcI3GHLfYCh8g6DwDYbCLxiKGleEA0Dw6FmNe1BELRTDh24qoigCgHTeRX/yWxjOHYWlKAfW4mynri51RGcpfBTeAYAoQu7lD4V/KBT+oVD6OZ4VfiFQt4qWPufTfZg0hxYReQbDhzzCbqmEpTALlsJMmAszYSnMdLwvzkbr6c9BHe4Ii8qs86g4/ZP0OUGlhTIwHMqAMKhbRV1aLlei/RPvO017QkQ3L4YPNRpRFGGrKIWlMBMQ7dKFjpaSXKRvmFfv5ywluVL4eHcdDHXraCgDwqEMaAWZzrfOq9wBMHiImhCGDzUYc34aTBnJMOelwZzveFRfz6JuE4M2DzguHFb4BkNQaqDwDYIyqA1UwW2gDKp6BIZLk1ACgKZNDDRtYjxyPETUeBg+dE1EmxWWwkxU5v0JS34afHqMkGbjLTv2FcoO73PaXlDroAqOgDr80nRJgkyO9ovf47UpRC0Yw4euqPpCy8rcVJhzLsCcl+Y0tb0yKEIKH227rrAb9VCFtoMqJBKq0HaQ+wTV2U3G4CFq2Rg+BACwVxpQmZOKypwLEM0mBAyeDMAREoVfveM04aTCPwyqsPZQhbSFqsa0716d4qXRZ0REV8LwaaEsJXkwnj8GU9ZZVGadhaUgE0DVMGelBv4D73bcyEomh1//CZBrvKAK6wB1WHvIatwfhYjoejB8mjlRFGErK4Ap6yyUAeFQt+oAADCmnnCeOkaugCqkHdStOjiGMNvtQFXXWOBfpniidCJqxhg+zYxot6EyJxWm9NMwpZ9BZUYybBUlAAC/+LFS+GjaxsK761+gbn0L1K07Qh3WwWnafyKixsTwaeLslkqnOzxmvbcMlZkpTtvINN5Qt75Fml4GAFRBbRA6foFbayUiqsbwaWJspgqY0v9wPNL+QGX2eUTOfQPKgFYAAHWrKNgMZdBE3gpt21ioIzpDGRhe74WZRESewPBpAmzGcpQe+hjG1BOozD4PiPYaawWY8/6Uwido9IMI5jBmIrrJMXxuMqJohzn3IqwledKNygSFEqW/fOy4vkYmh7pNZ2jb3uqYxj+ik9PoM14/Q0RNAcPnJmApzYPx/O8wXjwB48Uk2I16CGod2sfcBkEmh0ypRuDI/welfyg0bW+FTKX1dMlERDeE4eNBxj9PoWD/f2HJT3NarvANhrZDd9grDdI8Z35xf/VEiUREjYLh4yY2ox7G879DUCil7jS5zgeW/DQIKi20HbpD16EHtB26QxHQigMEiKhZY/g0ElEUYclPg+HcbzCcOwpTRjIg2qEOj5bCRxkcifDpz0ETEQNBzmtsiKjlYPg0gvJTP6Dou+2wFudcWiiTQ9OuG7xuiYMoihAEAYIgQNuui+cKJSLyEIbPDRJFOyozUyBT6aAKbetYKMhgLc6BTOcLXcc46G7pDV2HHpCpdZ4tlojoJsHwuQ6i3QZTxhlU/PELKs78Alt5Eby7DUXouEcBALqOfRzdaZGdOfSZiKgODJ9rYMo6h/LjX6Mi+ZA0XxrgGJ2m8A+V3stUGnanERFdAcPnKkTRDkGQAQAMKb+i7Oh+AI572njF9odX5/5Qh0dzdBoR0TVg+NTBZihD+emfUH7yO3jdOgD+fccBAHy6D4Vot8H71kFQhbVn4BARXSeGTxW71QzD2d9QfvI7GM4fBew2xwqZXAofZWBrBA2f4cEqiYiaB4YPgNIjn6P4uw9gN1U4FggyaKN7wafbEOhieFtoIqKGxvCBY6YBu6kCqlZR8Ok2BF63DoLC29/TZRERNVsMHwC6W+IQ8ffXoQpp6+lSiIhaBJmnC7gZyJRqBg8RkRsxfIiIyO0YPkRE5HYMHyIicjuGDxERuR3Dh4iI3I7hQ0REbtfkr/Ox2RzT4OTk5FxlSyIicpfqv8nVf6Mv1+TDJz8/HwAwbdo0D1dCRESXy8/PR7t27WotF0RRFD1QT4MxmUxISkpCSEgI5HLeuI2I6GZgs9mQn5+Prl27QqPR1Frf5MOHiIiaHg44ICIit2P4EBGR2zF8iIjI7Rg+RETkdgwfIiJyO4YPERG5HcOHiIjcrsWHz5kzZzBlyhT07NkTd911F06cOOHpkhrE7t270aVLF/Tq1Ut6JCYmerqs63bixAn0799fem82m/HMM88gPj4e/fr1w1tvveXB6m5MXcfWtWtXp5/dgw8+6MEKr81PP/2Eu+++G71798aoUaOwfft2AEBZWRnmz5+PPn36YPDgwdizZ4+HK7029R1Xeno6Onfu7PTzevrppz1c7bX55ptvcNddd6FXr14YOXKkdGyN+nsmtmCVlZXisGHDxC1btohms1n85JNPxLi4OFGv13u6tBu2cuVK8ZVXXvF0GTfMbreLO3fuFPv06SP26dNHWv7qq6+K06ZNE0tKSsT09HRxzJgxYmJiogcrvXb1HduJEyfEAQMGeLCy65eVlSX26tVL/OKLL0SbzSYeP35cvO2228Tvv/9eXLhwobhgwQLRYDCIp0+fFvv37y8eOnTI0yW75ErH9dlnn4mTJ0/2dInXLTc3V+zatav47bffiqIoiklJSWK3bt3EpKSkRv09a9Etn19//RUWiwUPPPAAlEol7rzzTnTs2BGffvqpp0u7YadOnUJsbKyny7hha9euxYcffoiHH37YaXliYiLmzp0LPz8/REREYNasWdK/1pqK+o7t1KlT6Ny5s4equjGZmZkYO3YsRo0aBZlMhu7duyM+Ph4HDx7E/v37sWDBAmi1WsTGxmLy5MnYsWOHp0t2SX3HdfTo0Sb98wKA0NBQHDx4EEOGDIHdbkdJSQnkcjm8vLwa9fesRYfPuXPnEB0d7bQsKioKKSkpHqqoYdhsNiQnJ+Ojjz7CoEGDMGrUKGzcuBFiE5xJ6b777sPevXvRtWtXaVlZWRny8/PRsWNHaVmHDh2a3M+trmMDgNOnT6OoqAh33XUXBgwYgPnz5yM3N9dDVV6buLg4PPfcc9L7kpISHDlyBKGhoRAEAe3bt5fWNaXftfqO69Zbb8Xp06dx9uxZjBkzBoMGDcLTTz+NsrIyD1Z77by9vWE0GtGtWzc8+OCDmDZtGgIDAxv196xFh4/BYKg14Z1Wq4XRaPRQRQ2jqKgIXbt2xYQJE/D1119L/8L+4IMPPF3aNQsLC6u1zGAwAIDTz06r1cJkMrmtroZQ17EBjmPp3bs33n33XXz++efQaDR45JFH3FzdjdPr9Xj44YfRo0cPdOnSBSqVCoIgSOub4s8McD6uESNGwNfXFwMHDsTu3buRmJiI7OxsPPPMM54u85qp1WocO3YMu3fvxp49e/Dee+8BaLzfsyZ/S4UbodPpUFlZ6bTMaDRCp9N5qKKGERISgm3btknvY2NjMX36dHzxxRfN4tYTWq0WAJx+ds3h51btySefdHq/dOlS9O/fH9nZ2QgPD/dQVdcmNTUV8+bNQ8eOHfHqq6/i/PnzMJvNEEVRCqCm+DO7/LhkMhlef/11ab2Pjw8WLlyIv/3tb7BarVAoms6fWJlMBpVKhW7duuHee+9FUlISgMb7PWvRLZ/o6GikpqY6Lbtw4YJTM7MpOnv2LNauXeu0zGKxQK1We6iihuXn54eQkBBcuHBBWpaamtrkf27V3njjDZw/f156b7FYAKDJ/PwOHz6Me++9FyNHjsTatWuhVqvRrl07iKKIjIwMabum9rtW13EZjUasWbMGBQUF0nYWiwUKhaLJ3OLl119/xd133+20zGw2w9fXt1F/z1p0+PTt2xeiKOKdd96BxWLBvn37kJycjFGjRnm6tBvi6+uLLVu2YOfOnbDb7UhKSsLWrVtr/Q/WlI0bNw4bNmxAUVERMjIysHnzZowbN87TZTWI5ORkrF69GmVlZSgrK8OqVaswdOhQBAYGerq0q0pLS8OcOXMwf/58LFq0SGrleHl5YeTIkUhISEBFRQXOnDmDXbt2NZmfWX3HpdVq8fPPP+O1116DyWRCXl4eEhISMHHiRKcuxptZbGwscnNzsWXLFthsNhw9ehR79uzBPffc07i/Zw0yZq4JS05OFu+77z6xZ8+e4tixY8Wff/7Z0yU1iJ9//lmcOHGi2LNnT3HYsGHitm3bPF3SDfnll1+chiObTCZxxYoVYv/+/cW+ffuKr732mmi32z1Y4fW7/NiKiorExx9/XIyPjxf79OkjPv7442JJSYkHK3Tdiy++KMbExIg9e/Z0eqxZs0YsKSkRFy5cKMbHx4uDBg0S3333XU+X67IrHVd6ero4Z84cMS4uToyPjxdXrlwpmkwmT5d8TZKSksT77rtP7N27t3jnnXeK+/fvF0WxcX/PeDM5IiJyuxbd7UZERJ7B8CEiIrdj+BARkdsxfIiIyO0YPkRE5HYMHyIicjuGDzU7S5cuRadOnep9rFu3DocOHUKnTp1QUVHhtrpsNhvuvffeWrNquKpTp0745ptvGqSWoqIi/O9//2uQfW3evBnr169vkH1Ry8HrfKjZ0ev10uSHqampmDFjBnbt2iXNi6bT6aBUKlFaWorg4GC3XYn+zjvv4MKFC06zI1+L/Px8+Pn5QaVS3XAtTz75JCoqKmpNw3Q9zGYz7rjjDmzatAkdOnS44f1Ry8CWDzU7Pj4+CAkJQUhICPz9/QEAgYGB0jIvLy+oVCqEhIS4LXgqKyvx1ltvYfr06de9j5CQkAYJHgANensNlUqFCRMmNOm7yZL7MXyoRbq8261Tp07Yt28fxo0bh27dumHq1KnIyMjAypUr0bt3b/zlL3/B//3f/0mfLy8vl24v3Ldv36vec2ffvn3w9fVFTEwMACAjIwOdOnXC119/jdGjR6N79+6YO3cucnJysHDhQvTs2ROjRo3Cjz/+KO2jZrfbjBkzsHbtWsybNw89evTAmDFjsGvXLmnb4cOHO81sXv19KSkpWLduHRITE7F//3506tQJgGMyzFdffRUDBw5Enz59MGvWLKcJJY8cOYJ77rkH3bt3x6BBg7BmzRrYbDZp/ejRo/HJJ5+gqKjoun4e1PIwfIiqJCQk4KmnnsLOnTuRnZ2Nu+++G97e3ti9ezdGjx6N5cuXS2H17LPPIjU1Ff/973+xdetWCIKA2bNnw2q11rnvb7/9FoMHD661/I033sCrr76Kt99+G4cOHcK4cePQq1cv7NmzB7GxsVi2bFm99W7atAkDBw5EYmIiBgwYgBUrViA/P/+qx/nggw/ir3/9K4YNGyaF29q1a/Hdd9/h9ddfx86dO9GhQwfMmDEDer0eNpsNjzzyCIYOHYpPP/0Ua9aswa5du7B3715pnzExMQgMDMTBgwev+v1EAMOHSPK3v/0N/fr1Q2xsLIYOHQqdTofHH38cUVFReOCBB2AymZCRkYH09HTs27cPCQkJ6N69O2JiYvDKK68gIyMDP/zwQ537TkpKqnMq+rlz56J79+6Ii4tDXFwcoqOjcf/99yM6OhrTpk1DdnY2ysvL69xnfHw8pk2bhqioKCxatAhWqxXJyclXPU4vLy9oNBqp69FkMuGdd97BihUrEB8fj+joaCxbtgze3t746KOPoNfrUVJSgtDQUERERGDAgAHYvHkzBgwY4LTfjh07SveAIbqapnOnI6JG1q5dO+m1VqtFmzZtpHNC1ffSMZvNyMrKAgDcfvvtTp83Go1ITU3FsGHDau27sLAQAQEBV/3OmtvU/M661Dy57+3tDQD1tryuJC0tDWazGbNmzXI6B1ZZWYnU1FT4+/tjzpw5eOaZZ7Bu3ToMGTIEd955J7p37+60H39/fxQWFl7z91PLxPAhqnL5XSdlsro7Bmw2G5RKJRITE2sNWPDz86vzMzKZrM6T/K5+Z12USmWtZfUNJKh5fqa+dW+//TaCgoKc1lWH2uOPP46JEyfiwIED+O677/Dggw9iwYIFmDt3rtN+msoN78jz2O1GdI2ioqJgsVhgNBrRrl07tGvXDiEhIVizZg0uXrxY52eCg4PdejJeqVRCr9dL79PT053W1wzNtm3bQqFQoKioSDqetm3bYt26dThx4gSysrKwfPlytGnTBrNnz8bWrVvx0EMP4eOPP3baZ3FxMYKDgxv3wKjZYPgQXaOoqCgMHz4cTzzxBI4cOYLz589jyZIlOH78OKKiour8TJcuXXDmzBm31ditWzfs2LEDp0+fxvHjx/Gvf/3LKXB0Oh0yMzORmZkJLy8vTJ06FatWrcJ3332HP//8E8899xwOHDiA6OhoBAYGYv/+/Vi1ahX+/PNPnDp1Cj/99BO6devm9J3Jycm1lhHVh+FDdB1efvlldO3aFY888gjuuece6PV6vP322/D19a1z+6FDh+Lw4cNuq2/hwoVo3749pkyZgieeeAKPPvqoU5fexIkTUVBQgDvuuAP5+fl44oknMGbMGDz11FMYN24cTp8+jU2bNqFt27bQaDR46623kJKSggkTJmDmzJmIiYlxGol39uxZVFRU1BqEQFQfznBA5AZGoxHDhw/Hf//7X3Tp0sXT5TS4hIQEFBYW4sUXX/R0KdREsOVD5AZarRazZs3C+++/7+lSGpzJZMInn3yC2bNne7oUakIYPkRu8sADDyA5Odlp5oDmYNu2bZgwYUK957uI6sJuNyIicju2fIiIyO0YPkRE5HYMHyIicjuGDxERuR3Dh4iI3I7hQ0REbvf/Af9+swOQ03SRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(tea_results.T, label='tea')\n",
    "plot(milk_results.T, '--', label='milk')\n",
    "\n",
    "decorate(xlabel='Time (minutes)',\n",
    "         ylabel='Temperature (C)',\n",
    "         loc='center left')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "当我们混合热茶和牛奶，看看温度怎么样。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "63.039907065852425"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mix_last = mix(coffee, milk)\n",
    "mix_last.temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.6113357534247257"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mix_last.temp - mix_first.temp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "有没有可能在30分钟的某个时间点混合是最佳的呢？\n",
    "\n",
    "为了找到这个结果，我们编写下面的函数，以牛奶加入时间作为输入变量，返回最终温度。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_and_mix(t_add, t_total):\n",
    "    \"\"\"Simulates two liquids and them mixes them at t_add.\n",
    "    \n",
    "    t_add: time in minutes\n",
    "    t_total: total time to simulate, min\n",
    "    \n",
    "    returns: final temperature\n",
    "    \"\"\"\n",
    "    tea = make_system(T_init=90, r=r_tea, volume=300, t_end=t_add)\n",
    "    tea_results = run_and_set(tea)\n",
    "    \n",
    "    milk = make_system(T_init=5, r=r_milk, volume=50, t_end=t_add)\n",
    "    milk_results = run_and_set(milk)\n",
    "    \n",
    "    mixture = mix(tea, milk)\n",
    "    mixture.t_end = t_total - t_add\n",
    "    results = run_and_set(mixture)\n",
    "\n",
    "    return mixture.temp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们尝试一些混合时间"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "61.4285713124277"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "run_and_mix(t_add=0, t_total=30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "62.722755084017976"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "run_and_mix(t_add=15, t_total=30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "63.039907065852425"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "run_and_mix(t_add=30, t_total=30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "然后再进行混合时间的扫描。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "sweep = SweepSeries()\n",
    "for t_add in linspace(0, 30, 11):\n",
    "    sweep[t_add] = run_and_mix(t_add, 30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面是最终的仿真结果。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving figure to file figs/chap16-fig02.pdf\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAEPCAYAAACdhMnXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeVxU9f7H8RfbDLuKAm6kAikqqCDiilpgVoqlKZlL+VO8Ji6lpealxeVqalmZS6ZXc00r3MrlKqKl5m6KIoohkAgEKAHKNsCc3x9e50qIjAsMy+f5ePjIOefM93wONPP2nPM936+RoigKQgghRAUyNnQBQgghah4JHyGEEBVOwkcIIUSFk/ARQghR4SR8hBBCVDhTQxdQmeTl5REZGYm9vT0mJiaGLkcIIaqsoqIi0tLScHd3x9zcvMR6CZ97REZGMnToUEOXIYQQ1cbGjRvx9vYusVzC5x729vbAnR9W/fr1DVyNEEJUXX/++SdDhw7Vfa/+nYTPPe5eaqtfvz6NGzc2cDVCCFH1lXYLQzocCCGEqHASPkIIISqchI8QQogKJ+EjhBCiwkn4CCGEqHASPkIIIYpRFIVb+be5npWMVqstl31IV2shhKiBCooKSMtJJ/X2DVJu3yAl+8adv//3v7mFeQAEugcwsPWLT3z/FRo+qampzJgxgxMnTqBWqwkMDOTtt98mLS2NGTNmcPLkSYyMjOjZsycffPABNjY2JdpQFIUvvviC77//noKCAgYMGMDUqVMxNb1zKLt37+bzzz8nLS2NDh06MG/ePOrWrVuRhymEEAanKApZ+bdIuX2D1OySAZOek4FC6XOJWpiaU9/Gnlb2ruVSX4WGT3BwMK1bt+bo0aOkpqYyfPhwXFxc2LlzJ3Z2dhw+fJjCwkImTJjAokWLeP/990u08d133xEWFsa2bdtQqVSMGzeO5cuXM378eGJiYggJCWHlypV4eHjwySefMGnSJNatW1eRh1kj5eXlcfv2berVq2foUoSoMTRFBaRl3yw1YPIL80t9r5GREfaWdXG0qoeDdT0crerhaF0Ph//+11plhZGRUbnVXmHhExERQUJCAps2bcLMzAwnJyfWr1+PWq2md+/eAKhUKpKTk8nNzaV27dr3bWf79u288cYbuuFvJkyYwLRp0xg/fjw//vgjzz77rG4coXfeeYcOHToQHx9P06ZNK+Q4K4vVq1ezbNkyTExMGD58OLGxsXz22WeP1ebw4cPx8/NjxIgRJdYNGzaMN998E39//8fahxDifxRFITP/Fim304oFzN3/pudmPPD9VmYW/w0W++IBY12PepZ2mBobbgDlCgufyMhImjdvzpIlS9i6dStqtZohQ4YwcuRI3TaTJ09m9+7duLq6MmTIkPu2ExMTg4uLi+61s7MzqampZGRkEBMTg7u7u26dhYUFDRo04MqVKzUufDZv3sz06dN55ZVXKmR/6enpFbIfIaqrvII8ErKSuZaRyLXMJK5l3vnvrfzbpb7HxMiYelb3P3txsK6LtcqqAo/g4VRY+GRmZnLmzBl8fHwIDw8nNjaWoKAg7O3tCQgIAODjjz9mxowZ/POf/2TChAls3LixRDs5OTlYWFjoXt8dqjsvL6/Eurvrc3Nzy/HIKp/evXtz7do1Zs6cyalTp2jUqBGXLl1i2bJlLF68mPj4eHJycjh+/DgNGjRg8uTJujOWvXv3smLFChISElAUBX9/f2bNmoWZmVmp+xs3bhxJSUlMnjyZCRMmMHr0aMLDw1m0aBFJSUm4urry4Ycf0qpVKwBatGjB7Nmz+eqrr8jIyGDYsGF4eHgwf/58MjIyCAwMZNq0abptp0yZwpo1a8jPz2fAgAFMmTJFd49PiKqmSFtE8q3U/4ZLItcykkjITCIl+8Z9t7c0s6CBjUPJgLG2p65FbUwMePbyOCrsE6xSqbC2tmbChAkAuLm5MXDgQMLCwnTho1arUavVTJkyheeee46MjIwSl98sLCzIy8vTvb77d0tLyxLr7q63sqq86V8e9u7dy7PPPss///lP/P39Wbx4cbH1e/bsYdmyZXz55Zd89dVXzJgxAz8/P5KSkpg6dSqrVq3C29ubP/74g1dffZX9+/fzwgsvlLq/pUuXFtvfhQsXePfdd1m2bBkdOnRg165djBw5kn379mFra6urYdeuXcTExDBo0CC6d+/Ojh07+OOPPxg0aBAvv/wyLVq0AGD//v3s2LGDvLw8goKCqFu3Lv/4xz/K7wcoxBOgKArpuRm6gLl7JpOY9SeF2sIS25sYm9DYpj5OtRvxVK2GPFWrEU/Vbkhdizrleu/FUCosfJydncnNzUWj0aBSqYA7kw0BvPLKK0yZMoVOnToBoNFoMDU1xdLSskQ7rq6uxMXF0b59ewBiY2Oxt7fH1tZWt+6u3NxckpOTi12me5I+PrSUs8mR5dL233k2cGd693FPpC0PDw969uwJQL9+/Vi6dCnZ2dnY29uzc+dOnJycuHXrFunp6dSpU4fU1NSHaj80NJR+/frRuXNnAF566SW+/fZb/vOf/xAYGAjA6NGjsbS0pE2bNlhaWvLqq69ibW1N69atcXBw4Pr167rwmTx5sq7H4siRI1m/fr2Ej6hUcjS591wqSyQhM4lrmUlka3Luu729Vd1iAfNUrUY0sHE06D2YilZh4dO1a1fs7OyYP38+06ZNIy4ujtDQUGbMmMGvv/7K4sWLcXNzQ1EU5s+fz0svvaQLqXv169eP1atX07lzZywsLFi8eDEvvfQSAH379mXIkCGcOHECT09PFi5cSMuWLWnWrFlFHWaVcG/X87uXr7RaLWZmZmzZsoXQ0FDMzc1p1aoV+fn5KErp3THvJykpiRMnTrBr1y7dssLCQpKSknSv7z2jNTExKdat3tjYuNg+mzRpovu7o6MjN27c//KEEOWtsKiQpFspurOYu/dnbuTc/56ntcqqRMg0rtUASzOL+25fk1RY+KjVajZs2MDs2bPx9fVFpVIRFBRE79698fX1Zd68ebzwwguYmJjw/PPP88477+je6+npycyZM+nXrx+vvfYaN2/eZPDgweTl5fH888/z1ltvAXfuD8ydO5ePPvqIlJQU2rZty6JFi8rtmJ7UmUhlsWvXLn766Se2bNmCo6MjcCfsH5aDgwNvvPFGsd9hfHx8sW7YD3MZISUlRVdPUlKSTPQnKoRW0XI9M5nLN2KIvhFLfMZ1km6lUKQtKrGtmbEpjW0b4PTfgLkbNnXMa1XLS2ZPQoXetXVycmLFihUllltaWjJr1ixmzZp13/edPXtW93djY2MmTpzIxIkT77tt7969dV23xcO5desWJiYmqFQqCgoK2Lx5M9HR0RQUFJT5XjMzM27dugXAyy+/zMSJE/H396dNmzYcP36cMWPGsGLFCt2l1YexaNEivvjiCzIyMvjmm28YNmzYQ7chRFnyCzVcTY/n8o2rRN+4SvSNWHIKindWMsIIR2v7Emcz9a3tq+yNf0ORLkNCp3///pw4cQJ/f39UKhVeXl707duX33//vcz3DhgwgJkzZxIXF8fkyZN5//33ef/997l+/ToODg7MmjXrkYIHoGnTpgQEBJCfn8/QoUMZOnToI7UjxL0y87KIvhF7J2zSYoj96xpFSvFxzOpZ2uFWz4UW9VxwsWtC41oNMDdVG6ji6sVIedgL+tXY9evX8fPzIzw8XKbRriRatGjB9u3badmypaFLEVWYoigk30rh8o3YO5fR0q6SfLt4RxojjGhSuxEt6rngZn8ncOpZ2hmo4qqvrO9TOfMRQlQ7hUWFxP51jcs3ruouo/39YU21iQrXuk1xq+eKm70LT9dtJh0BKpCEjxCiyrutyebKjThd54CY9HgKiorfq6xlbotbPRfdZbSmdZxqVNfmykbCR1Rq0dHRhi5BVDKKopCWk87ltBjdWU1CZlKJ7RrZ1setnist6jnjZu+Ko1U96XlWiUj4CCEqvRs56ZxJvEBU2u9cvhHDX7mZxdabGpviYtfkzv2aei60qOeMjdraQNUKfUj4CCEqHUVR+CPjOqcSIzideJ64jIRi661VVjSv56y7jOZs1wSVSenjD4rKR8JHCFEpFGqLuJT2O6cSIziTeJ60e0YNUJuqaVu/Je3qt8KtnisNbR0xNjI2YLXicUn4CCEMJqcgl3PJUZxOjOBsciTZ9zzUWcvcFu+GbejQqA3ujm5yZlPNSPgIISpUek4Gp5MiOJV4nsjU6GLD1TSyrU+HRm3xbtgG17pN5eymGpPwEUKUK0VRSMhM0t2/ufrXH7p1RkZGuNVzwbtRWzo0aksDGwcDVioqkoSPEOKJK9IWcSkthtNJ5zmdGEFq9k3dOpWJGW3qt6JDwza0b+iBrbnNA1oS1ZWEjxDiicgtyCPizyhOJUbwW3JksblsaqltaN/QA+9GbWnj6IbKtOR0KaJmkfARQjyyv3IzOZ14ntNJEVxIiS42Q2dDG8f/Xk5rw9N2zTA2lvs34n8kfIQQelMUhetZybr7NzHp8bp1RhjRoq6zLnAa2sq8S6J0Ej5CiDLdzs/m5/jjhF89QuKtP3XLzUzMaOPoRodGbfFq6EFtc1sDVimqEgkfIcR9KYrC7zfjCLt6mKMJZ3QDddqorWnfwAPvRm1oU7+lzG8jHomEjxCimLyCPI5cO8W+mEPEZ1zXLW9bvyW9XLrTvqGHzNopHpuEjxACgGsZiey7eojD8SfJLcwDwEZlxTPOXfB38aW+tb2BKxTViYSPEDWYpqiA4wm/EXb1MNE3ruqWt6jnwnMu3eno5CnD2ohyIeEjRA30561U9sce4WDsUW5psgGwMDXHt6kPz7l056najQxcoaju9A4frVZLbGwsN2/exMTEhHr16tGkSZOHmpwpNTWVGTNmcOLECdRqNYGBgbz99tvcvHmTOXPmcOzYMRRFoXv37oSEhFCrVq0SbXh6ehZ7XVhYiEaj4dChQzg6OjJ16lT27NmDqen/Du3HH3/EyclJ7zqFqI6KtEWcSbpA2NVDRPx5Sbe8ae3GPOfanW5PdcDczNyAFYqapMzwOXXqFOvWrePo0aNkZ2frlhsZGWFra4uvry9DhgzBy8urzJ0FBwfTunVrjh49SmpqKsOHD8fFxYVdu3ZhY2NDeHg4hYWFTJ06lVmzZrFw4cISbZw9e1b398LCQl5//XV8fHxwdHQEICoqiqVLl9K9e3e9fgBCVHfpORnsjz1CeOwR3SRsZiZmdHFqz3Ou3XG1ayozfIoKV2r4xMfH8+GHH5KcnIy/vz9LlizBxcWF2rVroygK6enpXL58mVOnTjF58mScnJyYNWsWzZo1u297ERERJCQksGnTJszMzHBycmL9+vWoVCr27NnD+PHjsbS0BGDQoEHMnj27zOJXrlxJYWEhEydOBCAvL4/Y2Fhatmz5KD8LIaoNraLlQsplwmIOczrpPFpFC0ADGwd6uXSnZ9NOWKutDFylqMlKDZ8pU6Ywbtw4evbsed/1DRo0oEGDBjzzzDNMmTKF/fv3M2XKFEJDQ++7fWRkJM2bN2fJkiVs3boVtVrNkCFDGDlyJMuWLSu27f79+8sMkJSUFL7++ms2bdqkG7bj0qVLmJiY8MEHH3Du3DkaNGjAxIkTeeaZZx7YlhDVRVb+bX6OO0rY1SOk3E4DwMTImE6NvXjO1ZfWDi3kLEdUCqWGz/fff6/3/6RGRkb06tULf3//UrfJzMzkzJkz+Pj4EB4eTmxsLEFBQdjb2xMQEKDbbtWqVYSHh7Np06YH7nPNmjX4+voWC6ns7Gy8vb0ZN24cbm5uHDhwgLfffpvvvvsONzc3vY5FiKpGURSib8QSdvUQxxN+o+C/46vVtayDv3M3nnXuSh2LkvdPhTCkUsPnbvBkZGRga2tbbFDAixcv0rBhQ+rUqXPf99yPSqXC2tqaCRMmAODm5sbAgQMJCwsjICCAgoICZs+ezcGDB1m7di0uLi6ltlVUVMT27dv59NNPiy3v1q0b3bp1073u3bs3W7duJTw8XMJHVDs5Bbkcjj9J2NXDXMtMBO6Mr+bZwJ1eLr54NXCXwTxFpfXADgerV69m8eLFrF69ulgvs/nz53Pu3DmmTp3KsGHD9NqRs7Mzubm5aDQaVKo7w6kXFd2ZwfD27duMHTuW7OxsQkNDdZ0HSnO300Hnzp2LLT9w4ACZmZn0799ft6ygoAC1Wob/ENXHtYxE/vP7zxy+dor8wnwAbNXWPOvcFX/nbjhY1zNwhUKUrdTw2blzJ5999hlvvfUWLVq0KLZu2bJlbNq0ifnz52Nvb0/v3r3L3FHXrl2xs7Nj/vz5TJs2jbi4OEJDQ5kxYwaTJ09Gq9WyceNGLCwsymzr3LlztG3btsS/6rRaLXPmzMHFxYXWrVuze/duzp49y7/+9a8y2xSiskvLvsl3F37i8B8nUVAAaGX/NL1cffFp1A4zeRhUVCGlhs+aNWuYMmUKb7zxRol11tbWjB49msLCQlavXq1X+KjVajZs2MDs2bPx9fVFpVIRFBREkyZN+OWXX1CpVHTp0kW3vY2NDYcOHSIpKYk+ffqwcuVKvL29AUhMTMTBoeR0u/7+/kyaNInJkydz48YNnJ2dWb58OQ0bNtTrhyFEZXQr/zbbov7Df2J+oVBbiImxCX7NuvL80z1pXKuBocsT4pEYKYqi3G+Fl5cX27dv56mnnir1zVevXiUwMJAzZ86UW4EV6fr16/j5+REeHk7jxo0NXY6o4TSFGnb/fpDtl/aSU5ALQLenOjDYo59cWhOVXlnfp6We+VhYWJCTk1PaauDOPRu5nyLEk1WkLeKX+ON8H7mT9NwMANo4tmRo2/40qyMjdYjqodTw8fb2ZseOHQ/sJbZt2zZatWpVLoUJUdMoisKZpPN8e34H17OSAWhW24mhbfvTpr48OC2ql1LDZ8yYMQwZMgQbGxtGjBihG30A7jxP880337BhwwZWr15dIYUKUZ1duRHLhoitXP7vyNL2VnV5zaMfXZ7yxthIukuL6qfU8GnVqhWLFi1i+vTpfPXVVzRr1gwbGxuysrKIi4vDzs6Ozz//nA4dOlRkvUJUK0lZf/Lt+R2cTDwH3Jk/55XWL9LLxVd6r4lq7YHP+fTo0YP9+/dz4MABLl++TFZWFnXq1KFNmzZ069ZN7vcI8Yj+ys3kh8idHIg7ilbRojIxo28LP/q1eA5LVdmPGwhR1ZUaPkVFRZiYmGBpaUnfvn3p27dvmY3dfY8Q4v5yCnL58fI+dkUfIL9Ig7GRMf7O3Rjo3gc7i9qGLk+IClNq+AwcOJDRo0fzwgsvlDnGW2FhIbt27eKbb75h+/btT7xIIaq6gqICwq4eZkvUHm7l3wbAp1E7XmvzEo1s6xu4OiEqXqnhs2zZMmbPns3cuXN59tln6d69O66urtSpU6fElAp79+6lTZs2LF26tCJrF6LS0ypajl47zeYLP5KafRO4M0X1sLb9aVGv9PELhajuSg2fBg0asGzZMi5evMiGDRv46KOPuHnzZrGzoLp169K9e3e+/vprPDw8KqRgIaqK839eYuP5bcT9lQBAI9v6DG3zMu0btpFpDUSNV+ZMpq1bt+bjjz8GICkpiRs3bmBsbIy9vX2ZA4AKURPF/ZXAxohtnE+5M1V1HYtavOoeQI+mnTAxlnuiQoAe4XOvhg0byjhpQpQi9fYNNl/4kSPXTgFgaWbByy1788LTz6A2VRm4OiEql4cKHyFESVn5t9l6cTd7rx6iSFuEqbEpz7v2oH+r57FRWxu6PCEqJQkfIR5RXmE+u68cYMflfeQW5GGEEd2bdCTQIwAHq7qGLk+ISk3CR4iHpCgKh+JP8O357fyVlwlAu/qtGNKmP03ryGjoQujjocInOzubP/74A1dXVzQaDdbWcklB1Cy3NdmsOPUtx6//BoBznacY1rY/7o4yTbsQD0Ov8NFoNMydO5fvv/8egL179/LJJ5+Qn5/Pp59+io2NTbkWKURlEJV6hcUn1nAz5y/MTdX8n2cgPZp1koE/hXgEen1qFi1axNmzZ9m4caNuPLegoCASExN13bCFqK4KtUVsOr+DmQe/4GbOXzxt15QFvUN4xrmLBI8Qj0ivM589e/awYMECPD09dcvatGnD7NmzCQ4OLrfihDC0P2+n8eWx1cSkx2OEEQNaPc/A1n0xled1hHgseoVPeno69vb2JZZbW1uTl5f3xIsSwtDudipY9dtm8grzqWtZhwkd/49WDk8bujQhqgW9rhl4e3uzbt26Yss0Gg3Lli3Dy8urXAoTwlByNLl8eXw1S0+uJa8wn05OXnzSO0SCR4gnSK8zn5CQEIKCgvj111/RaDRMnz6d+Ph4TExMZCZTUa1cTrvK4uOrSctJR22qZqRnID2bdZax2IR4wvQKn6ZNm7Jnzx527tzJ77//TlFREQEBAfTr1w8LC/0nvkpNTWXGjBmcOHECtVpNYGAgb7/9Njdv3mTOnDkcO3YMRVHo3r07ISEh1KpVq0QbGo0GLy8vzMz+N8ujp6enLgR3797N559/TlpaGh06dGDevHnUrSsP/IkHK9IWsSVqD1uidqMoCi51mjCx80ga2DgYujQhqiW9wqd///7MmzePAQMGPNbOgoODad26NUePHiU1NZXhw4fj4uLCrl27sLGxITw8nMLCQqZOncqsWbNYuHBhiTaio6OpVasWv/76a4l1MTExhISEsHLlSjw8PPjkk0+YNGlSiUuGQtwrNfsmi4+tJvpmLEYY8ZLbc7zqHoCpiTyDLUR50evTlZaWVuxM41FERESQkJDApk2bMDMzw8nJifXr16NSqdizZw/jx4/H0tISgEGDBjF79uz7tnPx4kXc3O7/QN+PP/7Is88+i7e3NwDvvPMOHTp0ID4+nqZNmz5W/aJ6OvLHKVae+ZbcgjzqWNRiQscR8sCoEBVAr/AJDAzkzTff5NVXX8XJyQlzc/Ni63v06FFmG5GRkTRv3pwlS5awdetW1Go1Q4YMYeTIkSxbtqzYtvv376dly5b3bScqKor09HQCAgK4efMm3t7ehISE4OjoSExMDO7u7rptLSwsaNCgAVeuXJHwEcXkFuSx6rfNHIo/AUCHRm15s8MwGQhUiAqiV/h89dVXAHz66acl1hkZGXHp0qUy28jMzOTMmTP4+PgQHh5ObGwsQUFB2NvbExAQoNtu1apVhIeHs2nTpvu2Y2FhgZeXF+PGjcPU1JR//etfjBs3jtDQUHJyckrcgzI3Nyc3N1efwxQ1xO834/jy2GpSsm+gMjFjhOcg/Jy7SacCISqQXuFz+fLlx96RSqXC2tqaCRMmAODm5sbAgQMJCwsjICCAgoICZs+ezcGDB1m7di0uLvefYnj69OnFXr/33nt07tyZ5ORkLCwsSjx3lJeXh5WV1WPXL6o+rVbL9st7+T5yJ1pFS9PajZnYeSSNbRsYujQhahy9wqesMwd9erw5OzuTm5uLRqNBpbozsVZRUREAt2/fZuzYsWRnZxMaGvrAGVIXLVpE3759deFUUFAAgFqtxtXVlbi4uGJ1Jycnlxpkoua4kZPO4uNruJT2OwB9m/vxWpuXMDN5vHuZQohHo1f4eHp6PvCShD6X3bp27YqdnR3z589n2rRpxMXFERoayowZM5g8eTJarZaNGzeWGWTR0dFERkbqesLNmTOHnj17YmdnR9++fRkyZAgnTpzA09OThQsX0rJlS5o1a6bPYYpq6njCb3x9agPZBbnUMrdlfMc3aFu/laHLEqJG0yt8/t5VuaioiGvXrrF69WqmTJmi147UajUbNmxg9uzZ+Pr6olKpCAoKokmTJvzyyy+oVCq6dOmi297GxoZDhw6RlJREnz59WLlyJd7e3syZM4d//etf9OrVi6KiInr06KHrGdeiRQvmzp3LRx99REpKCm3btmXRokX6/ixENZNXkMc3Z3/gYNxRALwauDPWZzi1zG0NXJkQwkhRFOVR33zs2DE++eQTtm7d+iRrMpjr16/j5+dHeHg4jRvLpGBV2dX0P/jy2GqSb6diZmzK8Hav0Nu1h3QqEKKClPV9+lhP0Tk4OBATE/M4TQjxRGkVLT9d3s/mCzsoUrQ41WrIW51G8lTtRoYuTQhxD73C55dffimx7NatW6xevbrUBz6FqGjpORksObGGyNRoAJ5/uifD2vRHZaoycGVCiL/TK3zGjBlTYpmZmRkeHh58+OGHT7woIR7WyevnWH5qA7c12diqrQn2eQOvhu5lv1EIYRAV9pyPEOUhv1DD2nOh7L96GIB29VsR7PM6tS1KDkorhKg89JrPx8/Pj4yMjBLLU1JS6Ny58xMvSgh9/JFxnff2fcz+q4cxNTbljXYDea/7OAkeIaqAUs98wsPDOXPmDACJiYl8+eWXJcZ0u3btWvlWJ0QpolJ/Z/7hZeQW5tHItj5vdRpF0zrSQ1GIqqLU8HFzc2Pt2rXc7YkdFRVVbGRrIyMjLC0tmTdvXvlXKcQ9ziVH8emvy9EUFdDZqT3BPq+jlk4FQlQppYZPo0aNdA+XTp8+nZCQEKytZcRfYVgnrp/li2OrKNIW8WyzLvzDeyjGxnpdPRZCVCJ6dTj4+OOPKSwsJCUlRTcem6IoaDQaLl68SN++fcu1SCEADsWfYNnJdWgVLS8+/QxveA6Sh0aFqKL0Cp+DBw8yffp0MjMzS6yztbWV8BHlbl/MIVad2YyCwoBWL/Cqe4AEjxBVmF7XKz777DO6du3KDz/8gJWVFevWrWPhwoXUq1dPnvMR5e7Hy2H8+8wmFBSGtunPYI9+EjxCVHF6nfnEx8ezaNEinJ2dadWqFTk5Obz44ouYmZnx1Vdf0adPn/KuU9RAiqLww8WdhF7cDcAor8H0frrsWXOFEJWfXmc+5ubmupu6TZs2JTr6zvAlrVu3LjZ/jhBPiqIorDu3hdCLuzEyMiLY53UJHiGqEb3Cx9vbm2XLlpGVlYWHhwdhYWEUFBRw8uRJ6QEnnjitVsuK09+y60o4JsYmTOocRM9m8jCzENWJXuHz3nvvceHCBbZs2ULfvn3Jzc2lffv2TJ8+neHDh5d3jaIGKdQWseTEGsJjj2BmYsbUbm/SycnL0GUJIZ4wve75FBUVsWfPHvLy8jA3N6hjdNkAACAASURBVOeHH37g5MmT1K5dm7Zt25Z3jaKG0BQV8MWxVZxOjMDcVM0032BaOzQ3dFlCiHKg15nP8OHDuXDhgm54HQsLC3r06CHBI56YvMJ8Fhz+itOJEViZWfBBz7ckeISoxvQ686lVqxa3bt0q71pEDZWjyeXjw0uJvnGVWmob3u85kSa1ZZw2IaozvcKnY8eOjBkzhk6dOvHUU0+hVquLrZ86dWq5FCeqv6z828z9ZTGxf12jrkUdPug5kYa29Q1dlhCinOkVPjExMbRr1468vDyuXLlSbJ087Cce1V+5mcz+eRHXs5JxtKrHB8+8jYNVXUOXJYSoAHqFz/r165/IzlJTU5kxYwYnTpxArVYTGBjI22+/zc2bN5kzZw7Hjh1DURS6d+9OSEgItWqVnJelrG2HDh3KhQsXMDEx0b3n9OnTxV4Lw0vLvsmsnxeRcjuNxrYNeL/nROwsahu6LCFEBdF7OODExEQWLFhAcHAwqampbN26lXPnzj3UzoKDg7G3t+fo0aN89913bN++nZ9++omQkBBMTEwIDw9n3759ZGVlMWvWrPu28aBtFUXh0qVLbNu2jbNnz+r+SPBULkm3UvgwfCEpt9NoVseJGc9OluARoobRK3wiIiLo27cv0dHRHDp0iPz8fC5dusSwYcM4cOCAXjuKiIggISGB999/H7VajZOTE+vXr8fHxwdjY2PGjx+PpaUltra2DBo0SDeR3b20Wu0Dt42Li0Or1dKsWbOH+BGIivRHxnU+OvAZN3P/okU9Fz7qOQlbtTyoLERNo1f4LFiwgLFjx7Jq1SrdhHIhISEEBwezaNEivXYUGRlJ8+bNWbJkCb6+vvj7+xMWFoajoyPLli2jSZMmum33799Py5YtSxZrbPzAbaOiorCysuKNN96gU6dOvPbaa5w9e1av+kT5i7kZz4yDn5OZl4WHoxshPSZgqbIwdFlCCAPQK3yioqJ4/vnnSywPCAggPj5erx1lZmZy5swZTE1NCQ8PZ8mSJaxevZqffvqp2HarVq0iPDycd999t8w2/75tYWEhbdq0YebMmRw6dIg+ffowevRoUlJS9KpRlJ+o1CvM+vkLsjU5eDdswzTfYMxN1WW/UQhRLekVPnXq1LlvyFy4cIG6dfXrnaRSqbC2tmbChAmoVCrc3NwYOHAgYWFhABQUFPDhhx+yZs0a1q5di4uLS6ltlbbtyy+/zFdffYWzszMqlYphw4ZRv359jh07pleNonycS77InENLyCvMp+tT3kzu+g9UJmZlv1EIUW3p1dtt6NChfPTRR7zzzjsAXLp0iYMHD7J06VJGjRql146cnZ3Jzc1Fo9GgUqkAdLOi3r59m7Fjx5KdnU1oaCiOjo6ltvOgbbdv346NjQ1+fn66ZQUFBSWeSxIV53jCbyw6vvrOtNfOXflH+yEy7bUQQr/wGTVqFFZWVnz22Wfk5uYyceJE6tWrx9ixY3njjTf02lHXrl2xs7Nj/vz5TJs2jbi4OEJDQ5kxYwaTJ09Gq9WyceNGLCwefA/gQdtmZGSwaNEiXFxcaNSoEWvWrCE7OxtfX1+9ahRP1qH4Eyw9uRZFUXix+bO80W6gPBcmhAD0DB+AwYMHM3jwYHJyctBqtQ89lYJarWbDhg3Mnj0bX19fVCoVQUFBNGnShF9++QWVSkWXLl1029vY2HDo0CGSkpLo06cPK1euxNra+oHbvvHGG2RlZfH666+TlZWFu7s7//73v2XaBwPYF/ML/z6zGYBXWr1IoHtfCR4hhI6RoiiKPhteu3aNH374gd9//x1jY2NatmzJoEGDqF+/+gyFcv36dfz8/AgPD6dxYxlb7FH9eHkfGyK2ATCsbX/6uT1n4IqEEBWtrO9TvS6+Hz16lD59+nD48GEcHByws7Nj//79vPjii5w/f/6JFy2qJkVR+O7CT7rgCWo/WIJHCHFfel12mz9/PiNGjNB1OLhrwYIFzJkzh++++65cihNVx91pr3ddCb8z7XWH1+nRrJOhyxJCVFJ6nfnEx8czYMCAEssDAwOJjo5+4kWJqkWr1fL16Y3Fpr2W4BFCPIhe4dOhQwf27NlTYvmvv/6Kl5dMcVyTFWqLWHziGw7E/oqZiRnTuo2Vaa+FEGXS67Kbu7s7y5cv59SpU7Rv3x5TU1OioqIIDw+nT58+LFiwQLetzO1TcyiKwrITa/n12mnMTdW85xtMK5l9VAihB73C58yZM7Rt25bCwkJOnDihW+7l5UVycjLJycmAzO1T0/wUvZ8j105hbqrmw55v41q3qaFLEkJUERU6n4+oPs7/eYmN5+/0ahvfcYQEjxDioej9kGlCQgJXr15Fo9EUW25kZESvXr2eeGGi8krNvsmiY6tQFIUBrZ7Hp3E7Q5ckhKhi9Aqff//73yxcuBAjIyNMTYu/xcjIiIiIiHIpTlQ+mkINC498zS1NNp4NWhPYOsDQJQkhqiC9wmfVqlVMmjSJoKAgGRSyBlMUhRWnvyUuIwFHa3smdPo/+f9BCPFI9PrmUBSFXr16yRdNDfef33/m0B8nUJuomNJ1DNYqK0OXJISoovRKk6FDh7JixYoS93tEzRGV+jvrzoUCMNbndZ6q3cjAFQkhqjK9Lrs9//zzDBs2jJ07d2Jvb1+iS3V4eHi5FCcqh5s5f/H50ZUUKVr6ufWiy1PtDV2SEKKK0yt8pkyZQqNGjXjxxRfLnG9HVC8FRQUs/HUFmfm38HBswWseLxm6JCFENaBX+MTGxrJjxw6aNWtW3vWISmbVb98Rkx6PvaUdb3UOwsTYxNAlCSGqAb3u+Xh4eBAXF1fetYhKZv/Vw7ox297t9ia2apmUTwjxZOh15hMQEMA///lP+vbti5OTU4lnfYYOHVouxQnDuXIjllW/3ZkqY4z3UJrVcTJwRUKI6kSv8FmxYgWWlpYcOHCgxDojIyMJn2rmr9xMFv66giJtES88/Qzdm3Y0dElCiGpGr/C5X+iI6qmwqJDPjq7kr7xMWto/zfB2rxi6JCFENaT3U6NarZZffvmFNWvWkJWVxfnz57l9+3Z51iYMYO25UKJvXMXOojaTugRhKh0MhBDlQK8zn9TUVEaNGkViYiJ5eXn4+fmxYsUKIiIi+Oabb3B1dS3vOkUF+DnuGHtjfsHU2JR3u46htrmtoUsSQlRTep35zJkzB1dXV44fP45arQbgk08+oV27dsydO1fvnaWmphIcHEz79u3p0qULX3zxBQA3b95k8uTJdO7cmU6dOjF16lQyMzPv24ZGo+GDDz7Ax8eHTp068fXXXxdbv379erp3746npyfvvvsuOTk5etdXk11N/4OVp78FIKj9YJkiQQhRrvQKn+PHjzNu3DhUKpVumYWFBW+//fZDjWgdHByMvb09R48e5bvvvmP79u389NNPhISEYGJiQnh4OPv27SMrK4tZs2bdt43FixcTFxdHWFgYoaGhbNu2je3btwNw+PBhli9fzqpVqzhy5AgajYbZs2frXV9NlZV3i09//ZoCbSH+Lr4869zV0CUJIao5ve/55Ofnl1j2119/YWZmptf7IyIiSEhI4P3330etVuPk5MT69evx8fHB2NiY8ePHY2lpia2tLYMGDeLMmTP3bWfbtm28+eab1KpVi8aNGzNq1Cg2b94MwPbt23nllVd4+umnsbKy4p133mHnzp1kZ2fre5g1TpG2iM+P/ZubOX/RvK4z/+c5yNAlCSFqAL3C57nnnmPevHmkpqbqxnW7fPkys2bNws/PT68dRUZG0rx5c5YsWYKvry/+/v6EhYXh6OjIsmXLaNKkiW7b/fv307JlyxJtZGVlkZaWVuweU7Nmzbhy5QoAMTExuLi46NY99dRTKIpCfHy8XjXWRBsjtnEx9Qq1zW2Z3HU0Zib6/WNCCCEeh17h895772Fvb0/37t3JycnhxRdfpH///jRu3Jj33ntPrx1lZmZy5swZTE1NCQ8PZ8mSJaxevZqffvqp2HarVq0iPDycd999t0Qbd+/fmJub65ZZWFiQl5enW3/v2HNGRkaoVCpyc3P1qrGmOfLHSXZeCcfEyJjJXf6BnUVtQ5ckhKgh9OrtZmVlxWeffcakSZO4evUqhYWFuLi4PNRYbyqVCmtrayZMmACAm5sbAwcOJCwsjICAAAoKCpg9ezYHDx5k7dq1xc5g7robLPdeAszNzcXS0lK3/m4QwZ15iDQaDVZWMu/M38X/dZ3lpzYAMMIzEDf7kj9vIYQoL3qFj5+fH1u2bMHJyQknp/8Ns5KSksLLL7/MsWPHymzD2dmZ3NxcNBqNruNCUVERALdv32bs2LFkZ2cTGhqKo6PjfduoVasW9vb2xMbG6raJi4vTXYZzdXUtNgbdtWvXUBSl2CU9Abfzs/n01+Voigro2bQzz7l2N3RJQogaptTwCQ8P1930T0xM5Msvvyx2uQvufLnrq2vXrtjZ2TF//nymTZtGXFwcoaGhzJgxg8mTJ6PVatm4cWOZUzb069ePpUuX0qJFC3Jycli1ahWvv/46AC+99BIhISH07t0bJycnFi5cSK9evXRnRuLOw8KLjq8iNfsmLnWaEOT9Won5mYQQoryVGj5ubm6sXbsWRVEAiIqKKtazzcjICEtLS+bNm6fXjtRqNRs2bGD27Nn4+vqiUqkICgqiSZMm/PLLL6hUKrp06aLb3sbGhkOHDpGUlESfPn1YuXIl3t7evPXWW8ybN4++ffui1Wp59dVXee211wDo0aMHY8eOJTg4mMzMTLp06VJql+2aanPkj0T8eQkbtTXvdP0HKulgIIQwACPlbro8wPTp0wkJCcHaunoPqX/9+nX8/PwIDw+ncePGhi7niTue8BufHV2JsZEx7/eYiLtjC0OXJISopsr6PtXrns/HH3/8xAsTFSshM4mlJ9cBMKztAAkeIYRB6f2Qqai6sjU5fHrka/IL8+n2VAf6NH/W0CUJIWo4CZ9qTqtoWXxiDcm3U2lSuzFjOgyTDgZCCIOT8KnmQi/u5rekC1irrJjSdQxqU1XZbxJCiHJW6j2fhxkVoKzu0cIwTidGEHpxF0ZGRrzVeSQO1vUMXZIQQgAPCB9PT0+9L89cunTpiRUknoykrD9ZfGINAK95vETb+q0MW5AQQtyj1PBZt25dRdYhnqDcgjw++fVrcgvy6NTYi5fcnjN0SUIIUUyp4ePj46NXAykpKU+sGPH4FEVh6cm1JGb9SWPbBoz1GS4dDIQQlY5ez/nExcWxYMECrl69qhuP7e6gnenp6URFRZVrkUJ/2y/t5eT1c1iaWTCl25tYmJmX/SYhhKhgevV2mzFjBsnJybz22mukpqYyZMgQunfvTkZGhswUWomcS77I5gs/AjCh0//RwMbBwBUJIcT96XXmExERwbp162jTpg27d++mbdu2jBo1iqZNm7Jz505eeeWV8q5TlCHldhqLjq9GQSHQvS/tG3oYuiQhhCiVXmc+iqJQr96dbrrOzs663m1+fn5yya0SyCvM59MjX5OtycG7YRsGtHrB0CUJIcQD6RU+bm5u7Nu3D7gzZ86pU6eAO50N9BiXVJSztWdD+SMzkQY2DozvOAJjI3l2WAhRuel12W38+PEEBwdjZmZGv379+Oqrrxg+fDixsbF07y4TkRnSlRuxhMcewcTYhHe7jsFSJQ/8CiEqP73Cx9fXl71791JUVISjoyObN29m+/bt+Pr66iZyExVPq9Wy6sxmAPq16IVTrYYGrkgIIfSjV/gANGz4vy+25s2bM3Xq1HIpSOgv7Oph4jISqGtZh/6tnjd0OUIIoTe9wic+Pp758+dz8eJFCgsLS9znOXbsWLkUJ0qXlXeLzRd2ADDCcxDmpmoDVySEEPrTK3w+/PBD0tLSGD16dLWfzbSq2Hh+O9kFubSt3wqfRu0MXY4QQjwUvcLnwoULrF+/Hnd39/KuR+jhyo1YDsYdxcTYhP/zCpThc4QQVY5efXIbNGhAfn5+edci9KDValn12/86GTS0cTRwRUII8fD0OvN56623mDFjBsHBwTRt2hQzM7Ni611dXculOFHS/tjDxP0lnQyEEFWb3uEDMGnSJN0yIyMjFEXByMhI7/l8UlNTmTFjBidOnECtVhMYGMjbb7+tW5+fn8/w4cP5xz/+gb+//33b6NOnD0lJSbrXRUVF5Ofns2nTJry8vFi0aBErV64sFpDLly+nY8eOetVYmWXl3WLTeelkIISo+vQKn/Dw8Ceys+DgYFq3bs3Ro0dJTU1l+PDhuLi4EBAQwOXLl/nwww+JiIh4YBu7du0q9nrSpEmYmpri5eUFQFRUFCEhIbz22mtPpObK5H+dDFpKJwMhRJWmV/g0atTosXcUERFBQkICmzZtwszMDCcnJ9avX49arSYuLo4RI0bw5ptvkpaWpnebP/74IxEREezcuVO37OLFi4wdO/ax661sincyeFU6GQghqrRSw6dz587s2rULOzs7OnXq9MAvO32e84mMjKR58+YsWbKErVu3olarGTJkCCNHjiQ7O5v9+/djbW2t9wyqeXl5LFiwgJkzZ2JpaQncuayXlpbGihUrOHfuHLVr12bUqFFVftTtezsZBLTwl04GQogqr9TwmTp1qu6ZnmnTpj32jjIzMzlz5gw+Pj6Eh4cTGxtLUFAQ9vb2BAQEPHR7W7ZswcHBAT8/P92y9PR0fHx8GDZsGIsWLeLs2bOMHTuWunXr0rNnz8c+BkO5t5OBjFgthKgOSg2fl156CWPjOz2x+/fv/9g7UqlUWFtbM2HCBODOSNkDBw4kLCzskcInNDS0xH0dNzc31q9fr3vt4+PDSy+9xL59+6ps+EgnAyFEdVTqcz6tW7fm5s2bxZadOnUKjUbzSDtydnYmNze32PvvTsn9sJKTk4mOjub554t3NT59+jRr1qwptqygoAC1uup+YX8rnQyEENVQqeFzv3l6xowZQ0pKyiPtqGvXrtjZ2TF//nw0Gg3R0dGEhobSp0+fh27r3LlzuLi4YGtrW2y5ubk5Cxcu5Oeff0ar1XL06FF27tz5RM7cDOHKjVgOSCcDIUQ19FCzjj3OxHFqtZoNGzaQkJCAr68vQUFBBAUF0bt37zLf6+npyY8//qh7nZiYiIODQ4nt3N3dWbBgAZ988gleXl7Mnj2bjz/+mDZt2jxy3YYinQyEENWZ3lMqPAlOTk6sWLHigdscOHCgxLKzZ88We303uO7nhRde4IUXqv5NeelkIISozmS+5Uro3k4Gb7QbKJ0MhBDVzgPPfLZu3ap7hgbudBDYsWMHderUKbbd0KFDy6e6GureTgYdG3sauhwhhHjiSg2fhg0bsmnTpmLL6tWrx9atW4stMzIykvB5gqSTgRCiJig1fO5370WUL+lkIISoKeSeTyUinQyEEDWFhE8lIZ0MhBA1iYRPJSGdDIQQNYmETyUgnQyEEDWNhI+BSScDIURNJOFjYNLJQAhRE0n4GFBW3i02XbgzZp10MhBC1CQSPgb07fntZGtyaOMonQyEEDWLhI+B3NvJYKRXoHQyEELUKBI+BlCik4FtfQNXJIQQFUvCxwCkk4EQoqaT8Klg0slACCEkfCqcdDIQQggJnwolnQyEEOIOCZ8KIp0MhBDifyR8Koh0MhBCiP+p0PBJTU0lODiY9u3b06VLF7744oti6/Pz8wkMDGT//v0PbKdnz560a9cOT09PPD096d27t27d8ePHCQgIoF27dgwePJhr166Vy7E8DOlkIIQQxZU6k2l5CA4OpnXr1hw9epTU1FSGDx+Oi4sLAQEBXL58mQ8//JCIiIgHtpGenk5KSgpnzpzB0tKyxLpx48Yxd+5cnn32WdauXcuoUaPYu3cvxsaGO8mTTgZCCFFchX0jR0REkJCQwPvvv49arcbJyYn169fTsWNH4uLiGDFiBC+++CINGzZ8YDsXL16kSZMmJYIHICwsjKeffprevXtjZmZGUFAQGo2GY8eOlddhlUk6GQghREkVFj6RkZE0b96cJUuW4Ovri7+/P2FhYTg4OODg4MD+/fsZMWJEmV/OUVFRKIrCwIED6dSpE6NGjeLq1asAxMTE4OLiUmz7Zs2aceXKlXI7rgeRTgZCCHF/FRY+mZmZnDlzBlNTU8LDw1myZAmrV6/mp59+wsrKCmtra73aMTY2xsPDgyVLlnDw4EFatmzJ6NGjyc3NJScnBwsLi2Lbm5ubk5ubWx6HVKb9sUekk4EQQtxHhd3zUalUWFtbM2HCBADc3NwYOHAgYWFhBAQE6N3O6NGji71+5513+Pbbb7l48SIWFhbk5eUVW5+Xl3ffS3Tl7U4ngx2AdDIQQoi/q7AzH2dnZ3Jzc9FoNLplRUVFD93OmjVrOH36dLE2ioqKUKlUuLq6EhcXV2z72NhYXF1dH73wRySdDIQQonQVFj5du3bFzs6O+fPno9FoiI6OJjQ0lD59+jxUO4mJicydO5fU1FTy8vKYN28eTZo0oXXr1vTq1YvLly+ze/duCgoKWLVqFcbGxvj4+JTTUd2fdDIQQogHq7DwUavVbNiwgYSEBHx9fQkKCiIoKKjYMzql8fT05Mcf7zwn8+6779K2bVv69+9P586dSUhIYPny5ZiYmFC3bl2WL1/OihUr8PHxYe/evSxfvhyVSlXeh6cjnQyEEKJsRoqiKIYuorK4fv06fn5+hIeH07hx40dqY1/MIf59ZhN1Lerw+Ysfyb0eIUSNVNb3qQyv8wRl5d/+XycDT+lkIIQQpZHweYKkk4EQQuhHwucJibkZz4HYX6WTgRBC6EHC5wmJ+DMKkE4GQgihjwodWLQ6e+HpZ3Cq1RDvhm0MXYoQQlR6Ej5PiKXKAp/G7QxdhhBCVAly2U0IIUSFk/ARQghR4SR8hBBCVDgJHyGEEBVOwkcIIUSFk/ARQghR4aSr9T3uzi/0559/GrgSIYSo2u5+j5Y2b5uEzz3S0tIAGDp0qIErEUKI6iEtLY0mTZqUWC5TKtwjLy+PyMhI7O3tMTExMXQ5QghRZRUVFZGWloa7uzvm5uYl1kv4CCGEqHDS4UAIIUSFk/ARQghR4SR8hBBCVDgJHyGEEBVOwkcIIUSFk/ARQghR4SR8hBBCVDgJnyfg8uXLvPrqq7Rr146AgADOnz9v6JIeS2hoKK1bt8bT01P3Z9u2bYYu66GdP3+ezp07615rNBo++OADfHx86NSpE19//bUBq3s49zsWd3f3Yr+jkSNHGrBC/fz6668MGDAALy8vevXqxebNmwHIyspi4sSJtG/fHl9fX7Zs2WLgSstW2rEkJCTg5uZW7HcTEhJi4GrLdvDgQQICAvD09MTf3193POX2uVHEY8nPz1eeeeYZ5ZtvvlE0Go2yc+dOxdvbW7l165ahS3tkM2fOVD755BNDl/HItFqt8v333yvt27dX2rdvr1v+6aefKkOHDlUyMjKUhIQEpXfv3sq2bdsMWGnZSjuW8+fPK126dDFgZQ8vKSlJ8fT0VPbt26cUFRUpERERSocOHZRDhw4pkyZNUt566y0lJydHiYqKUjp37qycOHHC0CWX6kHHsmfPHmXQoEGGLvGhpKSkKO7u7srPP/+sKIqiREZGKh4eHkpkZGS5fW7kzOcxnTx5koKCAkaMGIGZmRl9+vTB1dWV3bt3G7q0R3bx4kVatmxp6DIe2ZdffsmmTZsYO3ZsseXbtm3jzTffpFatWjRu3JhRo0bp/nVXWZV2LBcvXsTNzc1AVT2axMRE+vbtS69evTA2NqZNmzb4+Phw7Ngx9u7dy1tvvYWFhQUtW7Zk0KBBfPfdd4YuuVSlHctvv/1WJX83Dg4OHDt2jB49eqDVasnIyMDExAQrK6ty+9xI+DymmJgYXFxcii1zdnbmypUrBqro8RQVFREdHc2OHTvo1q0bvXr1YsWKFShVaBSmwYMHs3XrVtzd3XXLsrKySEtLw9XVVbesWbNmlf73dL9jAYiKiiI9PZ2AgAC6dOnCxIkTSUlJMVCV+vH29mbWrFm61xkZGZw+fRoHBweMjIxo2rSpbl1l/wyVdiytWrUiKiqK33//nd69e9OtWzdCQkLIysoyYLX6sba2Jjc3Fw8PD0aOHMnQoUOxs7Mrt8+NhM9jysnJKTFonoWFBbm5uQaq6PGkp6fj7u7Oyy+/zIEDB3T/8v72228NXZreHB0dSyzLyckBKPa7srCwIC8vr8LqehT3Oxa4U7uXlxdr167lP//5D+bm5owbN66Cq3t0t27dYuzYsbRt25bWrVujUqkwMjLSra8Kv5u77j0WPz8/bG1t6dq1K6GhoWzbto3k5GQ++OADQ5epF7VazdmzZwkNDWXLli2sW7cOKJ/PjUyp8JgsLS3Jz88vtiw3NxdLS0sDVfR47O3t2bBhg+51y5YtGTZsGPv27avSU01YWFgAFPtdVeXf0/Tp04u9fu+99+jcuTPJyck0aNDAQFXpJy4ujuDgYFxdXfn000+5evUqGo0GRVF0AVRVfjd/PxZjY2M+//xz3XobGxsmTZrEkCFDKCwsxNS0cn/lGhsbo1Kp8PDwIDAwkMjISKB8Pjdy5vOYXFxciIuLK7YsNja22GlqVfL777/z5ZdfFltWUFCAWq02UEVPRq1atbC3tyc2Nla3LC4ursr+nhYtWsTVq1d1rwsKCgAq/e/p1KlTBAYG4u/vz5dffolaraZJkyYoisL169d121WFz9D9jiU3N5cFCxZw48YN3XYFBQWYmppW6mlaTp48yYABA4ot02g02NraltvnRsLnMXXs2BFFUVizZg0FBQXs2rWL6OhoevXqZejSHomtrS3ffPMN33//PVqtlsjISNavX1/if8yqqF+/fixdupT09HSuX7/OqlWr6Nevn6HLeiTR0dHMmzePrKwssrKymDNnDj179sTOzs7QpZXq2rVrjBkzhokTJ/LOO+/oznKsrKzw9/dn4cKFZGdnc/nyZX744YdK/bsp7VgsLCw4evQon332Lh6gSAAAB31JREFUGXl5eaSmprJw4UL69+9f7LJiZdOyZUtSUlL45ptvKCoq4rfffmPLli0MHDiw/D43j91fTijR0dHK4MGDlXbt2il9+/ZVjh49auiSHsvRo0eV/v37K+3atVOeeeYZZcOGDYYu6ZEcP368WPfkvLw8ZcaMGUrnzp2Vjh07Kp999pmi1WoNWKH+/n4s6enpyuTJkxUfHx+lffv2yuTJk5WMjAwDVli2uXPnKs2bN1fatWtX7M+CBQuUjIwMZdKkSYqPj4/SrVs3Ze3atYYu94EedCwJCQnKmDFjFG9vb8XHx0eZOXOmkpeXZ+iSyxQZGakMHjxY8fLyUvr06aPs3btXUZTy+9zIZHJCCCEqnFx2E0IIUeEkfIQQQlQ4CR8hhBAVTsJHCCFEhZPwEUIIUeEkfIQQQlQ4CR9R7bz33nu0aNGi1D+LFy/mxIkTtGjRguzsbEOXe19XrlyhRYsWxZ76f5j1j9s+3BlkNjAwsMQIHvpq0aIFBw8e1GvbWbNmVck5o8Sjq9wDDQnxCEJCQnjnnXeAO0OBDB8+nB9++EE35pmlpSVmZmYcOXKkSowfZijr16/Hzc2NZs2aPdL7jxw5Qq1atfTadty4cQwcOJCePXtSp06dR9qfqFrkzEdUOzY2Ntjb22Nvb0/t2rUBsLOz0y2zsrJCpVJhb29fqYc8MaT8/Hy+/vprhg0b9sht2Nvbo1Kp9Nq2bt26dOrUqdigtqJ6k/ARNdLfL7u1aNGCXbt20a9fPzw8PHjttde4fv06M2fOxMvLi+7du7N9+3bd+2/fvq2bWrhjx45lzqcTHx/Pm2++ibe3N+7u7vTt25fw8HDd+vT0dMaPH4+npyfPPfccp06dKvb+staXVU9Z7/+7Xbv+v737C2mqjeMA/i2YsRq6i3XVxWDDZ6Rm61aJ/YHNIPDCiyhMlPRi1EUtQWFF8yqwmlAQXeywG0nRihHs3Y0wgoRqzIu0PyqrSeoc2Zyw2Goqz3tRHTbP8tSb7zHk97maz7/f7zw3D+ecx/P8g8rKSjDGAAALCwswmUyIRCJwOp2or6+Hy+VCKpWC2+2G2WyGw+HA+Pi4OEbxY7e2tjbcuXMH58+fx9GjR9HU1IQHDx6UxHQ6nRgeHhY/kkp2N1p8CPnO5/PB4/FgdHQUS0tLaGlpgUajwcOHD+F0OuH1esXF6tq1a0gkEhAEAYODg9izZw+6urqwvr4uGZdzDpfLhQMHDmBkZASPHz8GYwwejweFQgEAcPHiRXz69AlDQ0Po6+uDIAglY8jVy+Uj13+zJ0+e4Pjx45Ly27dv49atWwgEAnjx4gWam5tx7NgxPHr0CIcPH8bVq1d/Oqbf70djYyOCwSAaGhrQ19eH5eVlsb6hoQGrq6t4/fr1lrmRXeKPvw5HyF9sZmaGM8b4/Px8Sfnz5885Y4x//vyZc845Y4z7/X6x3uv1covFIn5AcX5+njPG+PT0NP/w4QNnjPFUKiW2//r1KzebzTwSiUhyyOVyXBAEnslkxLKpqSnOGOPJZJLH43HOGOOzs7NifSgUEvOWq5fLR65/OTabjQ8PD4t//7j+cDgslnV1dfHTp09L5jSbzYpz+mM+zp49y8+dOye2zWaznDHGnz59WhLX6XTywcHBsjmR3YU2HBDynV6vF3+r1WocOnRIfCf045ycQqGAZDIJADhx4kRJ/3w+j0QiAZvNVlKuVqtx5swZhEIhvHr1ColEAm/evAHwbUfZ7OwsKioqUF1dLfapr68Xf8vVx+PxLfP58uXLlv3LSafTZV/8b56j4jbFc1RO8cYFjUYDAJI7Ra1Wi3Q6vWVuZHegxYeQ7zafMrl3b/mn0hsbG1CpVAgGg5INC+V2d+VyOZw6dQr79u2Dw+GA3W7H/v370dbWVtKOF53kqVKpJOP8rF4un2fPnv3S+JuvnZf54P2vzlE5P7umYhsbG3/1oWtk+9A7H0J+k8FgwNraGvL5PPR6PfR6PQ4ePIgbN25gbm5O0j4ajWJubg5DQ0NwuVywWq3iSZecc5hMJhQKBbx9+1bsU/zeQ65eLh+5/uXodDqsrKz89tz8qUwmA51Op3hcojxafAj5TQaDAXa7HT09PYjFYnj37h16e3vx8uVLGAwGSXutVou1tTWEw2EsLi5ibGwM169fB/DtEZXBYIDFYsGVK1cwOTmJWCyGmzdvlsSTq98qH7n+5dTW1mJ6enqbZuzXZLNZJJNJHDlyRNG4ZGfQ4kPIf9Df34+6ujrxnyOz2SwCgQAqKyslbc1mM9xuN3w+H06ePIm7d++it7cXVVVV4h3IwMAAjEYj2tvb0d3dLXkkJ1cvl49c/82sVqvsduztNjExAZ1Oh5qaGkXjkp1BJ5kSQiTy+TzsdjsEQUBtba0iMS9fvgyj0YgLFy4oEo/sLLrzIYRIqNVqdHZ24v79+4rE+/jxI6LRKFpbWxWJR3YeLT6EkLI6OjowMzOD9+/f/++x7t27h0uXLomfQyK7Hz12I4QQoji68yGEEKI4WnwIIYQojhYfQgghiqPFhxBCiOJo8SGEEKI4WnwIIYQo7l9VQFGpwTfolAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(sweep, label='final temp', color='C2')\n",
    "decorate(xlabel='Time added (min)',\n",
    "         ylabel='Final temperature (C)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 分析\n",
    "我们再仔细看一下牛顿冷却定律这个公式：\n",
    "$${{dT}\\over{dt}}=-r(T-T_a)$$\n",
    "实际上这个微分方程是有解析解的，一般解为：\n",
    "$$T(t)=C_1e^{-rt}+T_{env}$$\n",
    "当$T(0)=T_{init}$时的特解为：\n",
    "$$T_{env}+(-T_{env}+T_{init})e^{-rt}$$\n",
    "我们观察到$T(t_{end})=T_{end}$，我们可以用特解来求解r。结果是：\n",
    "$$r={{1\\over{t_end}}log{T_{init}-T_{env}\\over{T_{end}-T_{env}}}}$$\n",
    "\n",
    "\n",
    "现在我们使用分析结果计算温度和时间之间的关系函数。下面的函数与run_simulation类似。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_analysis(system):\n",
    "    \"\"\"Computes temperature using the analytic solution.\n",
    "        \n",
    "    system: System object\n",
    "    \n",
    "    returns: TimeFrame\n",
    "    \"\"\"\n",
    "    T_env, r = system.T_env, system.r\n",
    "    \n",
    "    T_init = system.init.T    \n",
    "    ts = linspace(0, system.t_end)\n",
    "    \n",
    "    T_array = T_env + (T_init - T_env) * exp(-r * ts)\n",
    "    \n",
    "    # to be consistent with run_simulation,\n",
    "    # we put the array into a TimeFrame\n",
    "    return TimeFrame(T_array, index=ts, columns=['T'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "将$t_{end}=30$、$T_{end}=70$、$T_{init}=90$、$T_{env}=22$代入上面公式，可以得到r的估计值为0.0116。，这就是r_tea2的理论值。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>init</th>\n",
       "      <td>T    90\n",
       "dtype: int64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>r</th>\n",
       "      <td>0.0116102</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>volume</th>\n",
       "      <td>300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>temp</th>\n",
       "      <td>90</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_0</th>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_end</th>\n",
       "      <td>30</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>dt</th>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>T_env</th>\n",
       "      <td>22</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "init      T    90\n",
       "dtype: int64\n",
       "r                    0.0116102\n",
       "volume                     300\n",
       "temp                        90\n",
       "t_0                          0\n",
       "t_end                       30\n",
       "dt                           1\n",
       "T_env                       22\n",
       "dtype: object"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r_tea2 = 0.011610223142273859\n",
    "tea2 = make_system(T_init=90, r=r_tea2, volume=300, t_end=30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "70.0"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results = run_analysis(tea2)\n",
    "T_final_analysis = get_last_value(results.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "然后我们比较仿真的结果。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "69.99999985860761"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tea = make_system(T_init=90, r=r_tea, volume=300, t_end=30)\n",
    "results = run_simulation(tea, update_func)\n",
    "T_final_simulation = get_last_value(results.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "两者除了一些舍入误差，基本是一致的。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.4139239112864743e-07"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T_final_analysis - T_final_simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 练习\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}