finback-at/FMI10CS_BouncingBall.ipynb

## FMI10CS_BouncingBall.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FMI1.0 CS のCo-Simulation を行う\n",
    "上述の FMI_BouncingBall.ipynb では、BouncingBallモデルの入力信号である反発係数 $e$ をPythonの関数で作成しました。\n",
    "この反発係数を出力する部分をMSLに含まれる Modelica.Blocks.Sources.Step でモデル化し、FMUを作成します。\n",
    "\n",
    "FMI1.0 for Co-Simulationの簡単な例として、この反発係数を出力するモデル(model1)とBouncingBallモデル(model2)の２つのサブシステムモデルがデータ交換を行うCo-SimulationモデルをPythonスクリプトで構成します。\n",
    "\n",
    "\n",
    "JModelica.org で生成できる FMU はスタンドアロン型と呼ばれ、FMUの中にサブシステムモデルとソルバーが含まれています。FMI1.0 for Co-Simulation のサブシステムには次のような状態があります。\n",
    "* Instantiated (インスタンス化終了)\n",
    "* Initialized (初期化終了)\n",
    "* StepInProgress (コミュニケーションステップの計算中)\n",
    "* StepComplete (コミュニケーションステップの計算終了)\n",
    "* StepFailed または StepInComplete (コミュニケーションステップの計算がうまくいかない状態）\n",
    "* Terminated (シミュレーション終了)\n",
    "\n",
    "Pythonスクリプトで、load_fmu でサブシステムモデルをインスタンス化し、モデルの状態を遷移させる関数(model.initialize, model.do_step など)やデータ交換関数(model.get, model.set など)を使用して Co-Simulationプロセスを実現します。\n",
    "\n",
    "#### BouncingBall モデル\n",
    "BouncingBall モデルのソースコードを作成します。入力信号は反発係数$e$、出力信号は高さ$y=h$です。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting BouncingBall.mo\n"
     ]
    }
   ],
   "source": [
    "%%writefile BouncingBall.mo\n",
    "model BouncingBall\n",
    "  import SI = Modelica.SIunits;\n",
    "  Modelica.Blocks.Interfaces.RealInput e(start=0.8);\n",
    "  Modelica.Blocks.Interfaces.RealOutput y;\n",
    "  SI.Length h;\n",
    "  SI.Velocity v;\n",
    "  parameter SI.Acceleration g = Modelica.Constants.g_n;\n",
    "  parameter SI.Height h0 = 1.0;\n",
    "initial equation\n",
    "  h = h0;\n",
    "equation\n",
    "  v = der(h);\n",
    "  der(v) = -g;\n",
    "  y = h;\n",
    "  when h < 0 then\n",
    "    reinit(v, -e * pre(v));\n",
    "  end when;\n",
    "end BouncingBall;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "FMI1.0 for Co-Simualtion の規格に従う BouncingBall モデルの FMU を生成します。 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymodelica import compile_fmu\n",
    "fmu = compile_fmu(\"BouncingBall\",\"BouncingBall.mo\", version='1.0', target='cs',compile_to = \"BouncingBallCS10.fmu\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 反発係数を出力するモデル\n",
    "\n",
    "同様に、Modelica.Blocks.Souces.Step から入力信号モデルの FMU を生成します。 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "fmu = compile_fmu(\"Modelica.Blocks.Sources.Step\", version='1.0', target='cs', compile_to = \"Step.fmu\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Co-Simulation モデル\n",
    "\n",
    "シミュレーションの開始時刻、終了時刻、communication step を決めます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_start = 0.0\n",
    "t_end = 3.0\n",
    "dt = (t_end - t_start)/500"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Co-Simualation モデルを構成し、シミュレーションを実行します。このモデルでは、model1のコミュニケーションステップが終了してからmodel2のコミュニケーションステップを実行しています。より実用的なモデルでは、サブモデルのコミュニケーションステップが同時に進行するように非同期的な処理が必要です。また、StepFailed や StepIncomplete に対する処理も省いてエラーで停止します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# load models (instantiate)\n",
    "import pyfmi\n",
    "model1 = pyfmi.load_fmu(\"Step.fmu\")\n",
    "model2 = pyfmi.load_fmu(\"BouncingBallCS10.fmu\")\n",
    "t = t_start\n",
    "\n",
    "# initialize result data array\n",
    "t_res = []\n",
    "e_res = []\n",
    "h_res = []\n",
    "\n",
    "# initialize model1 with input parameters\n",
    "model1.set(\"height\", -0.3)\n",
    "model1.set(\"offset\", 1.0)\n",
    "model1.set(\"startTime\", 1.0)\n",
    "model1.initialize()\n",
    "y1 = model1.get(\"y\")\n",
    "\n",
    "# initialize model2 with initial input data\n",
    "model2.set(\"e\",y1[0])\n",
    "model2.initialize()\n",
    "y2 = model2.get(\"y\")\n",
    "\n",
    "t_res.append(t)\n",
    "e_res.append(y1[0])\n",
    "h_res.append(y2[0])\n",
    "\n",
    "while t <= t_end:\n",
    "    status = model1.do_step(t,dt)\n",
    "    if status == pyfmi.fmi.FMI_OK:\n",
    "        status = model2.do_step(t,dt)\n",
    "        if status == pyfmi.fmi.FMI_OK:\n",
    "            ## succeed to simulate communication step!\n",
    "            \n",
    "            # get data\n",
    "            y1 = model1.get(\"y\")\n",
    "            y2 = model2.get(\"y\")\n",
    "            \n",
    "            # set data\n",
    "            model2.set(\"e\",y1[0])\n",
    "            \n",
    "            # calculate next communication step\n",
    "            t += dt\n",
    "            \n",
    "            # restore result data\n",
    "            e_res.append(y1[0])\n",
    "            h_res.append(y2[0])\n",
    "            t_res.append(t)\n",
    "        else:\n",
    "            print(\"model2 error\")\n",
    "            break\n",
    "    else:\n",
    "        print(\"model1 error\")\n",
    "        break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### シミュレーション結果"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f1a3e684490>,\n",
       " <matplotlib.lines.Line2D at 0x7f1a3e6845d0>]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xd4XOWV+PHvGTV3yZIly7IlS7bcC9gYY5oD2BADCSR0EhJISIDdBZJNyP5IsmEJm2w2gW3JkmQhIZSEDkkcMM0002xs3HGVu1zl3mXLen9/vHOtsTySptw7d+bO+TyPn7FGV/eeO7LPvHPeJsYYlFJKBUvI7wCUUkq5T5O7UkoFkCZ3pZQKIE3uSikVQJrclVIqgDS5K6VUAGlyV0qpANLkrpRSAaTJXSmlAijXrwv36tXLVFdX+3V5pZTKSJ988sl2Y0xpR8f5ltyrq6uZM2eOX5dXSqmMJCLrYjlOyzJKKRVAmtyVUiqANLkrpVQAaXJXSqkA0uSulFIB1GFyF5FHRGSbiCxu4/siIr8UkToRWSgiY90PUymlVDxiabk/Ckxp5/sXA4PCf24BfpN8WEoppZLR4Th3Y8wMEalu55DLgceN3a9vpogUiUgfY8xml2I8wdJZr7Nr8RvsKahgbdEEDuUXe3GZrNf56C76755F98atNOV0YmvXoWzqPhpEoh4/ul8Rk4f3TnGUSqm2uDGJqS+wIeLr+vBzJyV3EbkF27qnqqoqoYvtWfE+Z214CIAmE+LPzefwi6br2E5RQudTJ+rJXu7KfZarQ++SJ8dO+N6q5j789NgNvN085oTnjYGKwk6a3JVKI24k92hNuai7bhtjHgIeAhg3blxCO3NP+Mp90PRDaFhO7oKnuPrjh7m6xzK49gmompDIKZWjfg48cxccaIDTvg5jvgylQ6FxH9RNZ+D7/80j2++H8bfAZ38GOfafz/dfXMT0pVt9Dl4pFcmN0TL1QGXE1/2ATS6ct225BdBnNEz5Gdz2HhR0h8cvhzUzPL1soK370L6GOXnwzTfh0gegYgzkdYZuZXDql+C29+HM2+Hjh+CFr0OzbdmHBJqbE3qvVkp5xI3kPhX4anjUzARgj1f19qjKhsHNb0DxAHjyOtiyKGWXDoyGFfDktdCjAm5+HfqcEv243Hz47E/hop/Akr/CtO8BkBMSmo0md6XSSSxDIZ8CPgKGiEi9iNwsIreJyG3hQ6YBq4E64GHg7z2Lti1dS+Arf4ZOPeDpL8OhXSkPIWMdOQBPX28/Dd3wAnQv7/hnzroDzroT5vwe5v2JkAjHtOWuVFqJZbTM9R183wD/4FpEiepeDtc8Dn+4GKb9E1z5sN8RZYbp98KOVXDTS1AURyf35Hth0zx4+buUDHkMze1KpZdgzVCtHA/n3gWLnoVl0/yOJv2tec/Wz8+4DarPie9nQzlwxUOQk8fn1/87zabZmxiVUgkJVnIHOPe7UDYcXvs+NDX6HU36ajoCf7sTetbApB8ldo4eFXDhfVTvm8ul5j1341NKJSV4yd3p9Nu11rZKVXSfPAo7V8PFP4f8romfZ+yNbOk6hG+HnoGjh10LTymVnOAld4CBF0DthfDu/XBgh9/RpJ/De+Hdn0P/c2DQRcmdKxTi3f530Fe2w2zt51AqXQQzuYMdrndkH3z0K78jST8fPQgHt8OF97W5nEA8NhaN591jo2HGA3B4jwsBKqWSFdzkXjYUhl8Os3+vCSdS436Y9VsY+jnod5orpwyFhPubroHDu225Rynlu+Amd4Bz/hEa99oEr6x5T9gkfPa3XDtlSITFZgDN1RNh5m9sZ61SylfBTu59ToGBk2Dmr+HoIb+j8d+xJvjo11A5wQ4bdUlOyJZ2jp15J+zbDIuec+3cSqnEBDu5g22hHmiAT//sdyT+W/IX2LMezr7T1dM6ZftjNedD75Hw4a/sUpFKKd8EP7nXTISSQTDnD35H4r85j0DPahh8saunzQln92MGOyGqYSmsn+nqNZRS8Ql+cheBcV+D+o+ze1Gx7Sth3Qcw9kYIuftrd8oyzcbAyCsgvzvMfczVayil4hP85A5wyvWQU5Ddrfe5j0EoF079suunlnDLvbkZOyFq9NW2DKYLuCnlm+xI7l2KbYty4bNw5KDf0aRe0xGY/xQMngLd3d8tKSdccz++7O/YG6HpMCzUjlWl/JIdyR3sZhNH9sHyLFxQbMWrdtLS2Bs9Of3x0TJOcq841Y5UmveEJ9dTSnUse5J7/3Oge0V2DtNb9Bx0LbXLMnigpSwTMUJm9HWwZaGt9SulUi57knsoBKOuhLrp2bXezOG9sPJ1GPHF43ueuq2lQzXiyRFfBAQWPe/JNZVS7cue5A4w6hpobrLjvbPF8mm2/j3yKs8uEXLGuUeObe/Rx64Rv+g5HfOulA+yK7mXj4LSodlVmln0PBRWuTojtbVQtLIMwKirYOcq2Dzfs2srpaLLruQuAiOusBNs9m/zOxrvHdwJq9+2I4VcWP2xLSeMc4807DII5cHiFzy7tlIquuxK7gBDLwUMLHvZ70i8t+JVW4Ya8QVPL+O03E/aJLtLsZ0hvPQlLc0olWLZl9x7j7BT8Je95Hck3lv2MvToC31O9fQyoWgdqo6hl8KuNdCwzNMYlFInyr7kLmLXMl/9brDXeT9yEOrehCGXeFqSgZYO1ZPKMmCvD9nxSUmpNJJ9yR1g2Oeh+SisfMPvSLyz+h1oOgRDL/H8UjltlWXAjprpe5omd6VSLDuTe7/x0LUs2KWZ5S9DQaGdvOWxUFsdqo6hl8KmubB3k+exKKWs7EzuoZDdGHrVW3YDi6Bpboblr8KgCyE33/PLhSIXDotmyKX2cfkrnseilLKyM7kDDJpsa+4b5/gdifs2z7dryQy6KCWXywn/K2qz5V46BIqq7JupUiolsje5DzgfJCeYdXcniXq0lkxrIq0WDjv5AKidbDuxdX9VpVIie5N75yI7a7MuoMm9fDR0K03J5XLamqEaqXayXZWz/uOUxKRUtsve5A424WxeAPu2+h2Jew7vhQ2zUtZqhzYWDmutZqLdLKRuemqCUirLZXdyH3ShfVz1pr9xuGnt+3ZWau2klF3y+AbZ7WX3gu5QOUGTu1IpElNyF5EpIrJcROpE5O4o368SkbdFZJ6ILBQR7wdXu6F8NHTpZceEB8WqNyGvK1SekbJLOmUZ09ESA7WT7D62+7akICqlsluHyV1EcoAHgYuB4cD1IjK81WH/DDxrjBkDXAf82u1APSECNefCmveCs/bJqrfsUru5BSm75Ek7MbXFKRWtmeFxREqpWFru44E6Y8xqY8wR4Gng8lbHGKBH+O+FQObMVqmZCPs2wc7VfkeSvJ1r7H2ksCQDEaNl2i26Y5dc7lSoyV2pFIglufcFNkR8XR9+LtK9wA0iUg9MA+5wJbpUqJ5oH9e8628cblj9tn1MYWcqtLTcO/zwE8qB/mfbfgGllKdiSe7RVp1q/d/4euBRY0w/4BLgCRE56dwicouIzBGROQ0NDfFH64WSgXZv1TXv+R1J8tZ+AN3KoaQ2pZcNxdKh6qg+164Suafe26CUynKxJPd6oDLi636cXHa5GXgWwBjzEdAJ6NX6RMaYh4wx44wx40pLUzMGu0NO3X1thtfdjYF1H0D12Z6vAtlaqKNJTJFqzrWPQXgzVSqNxZLcZwODRKRGRPKxHaZTWx2zHpgEICLDsMk9TZrmMaiZCAcaMnvN8V1rYN9m6H9Wyi/dUpaJIbmXjYDOPbU0o5THOkzuxpgm4HbgNWApdlTMpyJyn4hcFj7su8A3RWQB8BRwk4npf3qaqHZakxnc0bf2A/uYglUgW2vZiSmWg0PhunsGv9ZKZYDcWA4yxkzDdpRGPndPxN+XAGe7G1oK9ewPRf1tcj/jVr+jScy6D6BLiV2kK8U6XDistepz7XLLu9bZ114p5brsnqEaqebc8OzOWJqfaWjdB7Ykk+J6O7QMhYw5uTt1dy3NKOUZTe6OqjPh8G7YsdLvSOK3ewPsXu9LSQY62IkpmtJh9lPGWu1UVcormtwdlRPs4/qZ/saRiHXhenu1P5WxmBYOixQK2dd7wyzvglIqy2lyd5QMtK3JTEw46z6wMz/LWq8KkRpOJajdJX9bqxxvZ9Puz5xBVUplEk3uDhG72FZGJvePbFkplOPL5XM62kM1Gmdhs/rZHkSklNLkHqlyPOyogwPb/Y4kdgd32n6CyvG+hRDXJCZHxakQysvMN1OlMoAm90hO3T2TEs7GT+xjv9N9CyEUy05MreV1hj6nwAbdmUkpL2hyj1QxxrYmM6lTtX42SAgqxvoWQtwdqo7KM2DTXN1XVSkPaHKPlNfJlgsyqTW54WM7pb+gm28hxLVwWKTK8dB02G7goZRylSb31irPgE3zoKnR70g61txsyzL9xvkaRiiRDlVo6SfIpDKYUhlCk3trlWfAsUbYNN/vSDq2fQU07vW1MxVaJjHFndx7VEBhpSZ3pTygyb01J1FunONvHLGoD5ePfOxMhTgXDmutcrxN7hm0zpxSmUCTe2vdy6FHX9g41+9IOlY/GzoVQfFAX8MIxbtwWKR+4+1SxXs3uhuUUlkuplUhs07FGDuKI93Vz7Gt9pC/79FOy71hXyMrtu6L62c7dR5KFbBpyYfsH3CxB9EplX7KuhdQ1CXf02toco+m71i7JO2hXXZjiXR0eA9sWwrDv+B3JOSIkJ8b4tEP1/Loh2vj+tkCjrC4IIcXX36JB5q6ehOgUmnmJ18YyQ0TvF3uWpN7NM6Y8U3zUr7ZdMw2zQMM9DvN70gIhYTnbj2T+l2HEvr5g28O5tqCBmrO8W+svlKpNLJvD8+vock9moox9nHj3DRO7uHRPD5OXop0SmURp1QWJfbDa0+HJVO5dFS5L+vRKxVE2qEaTedwJ+WmeX5H0rZN86CoCroU+x1J8irG2rX0d631OxKlAkOTe1v6jk3vETOb50OfU/2Owh3OJ6VM6MRWKkNocm9LxVjYtwn2bfE7kpMd2mVbuRUBSe5lwyGnIL0/KSmVYTS5t6VvuJadjq33zQvsY1Ba7rn5UD4yM2YFK5UhNLm3pXw0SE56lgqOd6aO8TcON1WMtfeVqRuUK5VmNLm3Jb8LlA1Lz5b7pnlQGJDOVEfFGDiyz26WopRKmib39vQ5FbYsTL91TzbPD0693eGUwdLxk5JSGUiTe3v6jIYDDbB/q9+RtAhaZ6qj12DI66J1d6Vcosm9PeWj7OPmhf7GESlonamOUA70HgFbF/sdiVKBoMm9Pb1H2sctaZTcg9iZ6igflZ5lMKUykCb39nTqAT1r0iu5b54fvM5UR/kouyDang1+R6JUxtPk3pHyUem1x+eWxS3loqApH20f0+n1VipDaXLvSJ/RsHM1HN7rdyRw5CDsXGUn/ARR2XCQkCZ3pVygyb0jTmty66f+xgF2/XbT3NIXEDT5XeyCbZrclUpaTMldRKaIyHIRqRORu9s45hoRWSIin4rIk+6G6SOnBJIOCWdrOIagttyhpVNVKZWUDpO7iOQADwIXA8OB60VkeKtjBgHfB842xowAvu1BrP7o3ge69IItC/yOxH56yO8GRdV+R+Kd8lGwez0c2u13JEpltFha7uOBOmPMamPMEeBp4PJWx3wTeNAYswvAGLPN3TB9JJI+napbFtux4D7vmeqp42UwHe+uVDJiyRJ9gcixafXh5yINBgaLyAciMlNEpkQ7kYjcIiJzRGROQ0NDYhH7oXyUrXcfO+pfDMbYlntQ6+2OdCqDKZXBYknu0fY9az3LJBcYBJwHXA/8TkRO2nPNGPOQMWacMWZcaWlpvLH6p88pcOwINCzzL4bd66Fxj225B1n33tC1TJO7UkmKJbnXA5URX/cDNkU55q/GmKPGmDXAcmyyDwYnoW5b6l8MTpkiqGPcI2mnqlJJiyW5zwYGiUiNiOQD1wFTWx3zF+B8ABHphS3TrHYzUF+V1EIoz9/hkFsWA2LHggdd+UjYtszfMphSGa7D5G6MaQJuB14DlgLPGmM+FZH7ROSy8GGvATtEZAnwNvA9Y8wOr4JOuZw8u2qh3y334hoo6OZfDKlSNhyaj9rJY0qphOTGcpAxZhowrdVz90T83QDfCf8JprJhsGGWf9ffujj4nakO59PJtiVQOsTfWJTKUAEeU+ey3sPtglZ+LEPQuB92rsmOejvYT0kSgq1L/I5EqYylyT1Wx1uTPpRmGpYDJjvq7QB5newyBNs0uSuVKE3usTqe3H3oVG0Iv6GUDUv9tf1SNszfPg6lMpwm91gVVdmp/7603JdBTgH0rE79tf1SNtx2qB495HckSmUkTe6xErGtST/qwNuW2Tp0KCf11/ZL2TDAhEtSSql4aXKPR9lwWwdO9TZwDcugbGhqr+m3dJg4plQG0+Qej7LhcGgn7N+aums27rOjdEqzLLn3rLGlKD/6OJQKAE3u8egdMf46VZyyRDZ1pgLk5EKpzxPHlMpgmtzj4YyYSWXd3VmsLNta7hAug2lyVyoRmtzj0bWXXbEwlS33bUsht1N2jZRxlA2DvRt14w6lEqDJPV6lQ1I7gqNhGfQalF0jZRzOJyU/l1pWKkNpco9X6RDYviJ1I2a2LYPSLKu3O5x+Bp2pqlTcNLnHq3QoNO6FfVu8v9bhvbC3PvuGQTp69IO8LtCwwu9IlMo4mtzj1WuwfUxFqWB7OKlla8s9FLIlqe06kUmpeGlyj5ezBO32FLQmnZEi2bzsba8h2nJXKgExreeuInTrDQWFqelUbViWvSNlHKWDYdGzdtnjdN+o5FgTrH4H1r4HB3fYfyvVZ0PNefZTiFIppMk9XiItnape27Y0+9aUaa1X+FPLjpVQMcbfWNqz4nV49W7YuQpy8qFzTziwHd57wP4OL/0PqJnod5Qqi2hzIhGlg1NTc29Yln0zU1s73seRpqUZY2D6j+HJq+12jFc/Bnevh7tWwA82wpW/t3vBPvZ5mPFA6tclUllLk3sieg2BAw1wcKd31zi8107gycaZqZGKB4DkpGenqjEw7Xvw/n/C2BvhlndhxBcgr7P9fl5nGHUV/N2HMOpqeOtfYfq9voassoeWZRIR2alaNcGba2xfeeK1slVuvk3w6bj07wf/A7MfhjNvh4t+Ykt20eR3gSsetvsBfPDf0KMCzrg1tbGqrKMt90Q4CdfLhLMjnNxLBnl3jUyRqj6OeKx627bCR1zRfmJ3iNi6+5BL4LUfwIaPUxKmyl6a3BNRWAW5nb1NONtX2nJENo+UcfQabHdlOnbU70isgzvhL39nx+Bf/mDHid0RyoEv/AZ69IXnb7YjgJTyiCb3RIRC0KvW207VHSttYs/N9+4amaJ0CDQ3wc41fkdiTf8X2+dyxUO25BKPzkX25/ash7d+4k18SqHJPXGlQ70dwbG9zrYMVcuImXToVN00D+Y+AWfclvjQzKoJcPo3YNZvYctid+NTKkyTe6J6DbGtryMH3D93c7MdL11S6/65M9Hx4ZA+J3dj4JW7oUsJTPxecuc6/4fQqYeOnlGe0eSeqF7hxLtztfvn3rMBmg5ry91R0M0uIuZ3p+rSv8GGmTDpHlteSUaXYjjnO1D3BqyZ4U58SkXQ5J4op1W9o879c2/XkTInKR3sb8vdGJjxCygeCGNucOecZ9xqO1en36uTm5TrNLknqniAffQiuTvDILXl3qKk1n5K8isJrnwdtiyCc7/j3nIQeZ1h4l2w8RNY+74751QqTJN7ovK72lbXjlXun3v7Srs4WddS98+dqUpq7Tr6BxpSf21jYMb9UFgJo69199ynXA9desGHv3T3vCrraXJPRkmtdy33XrWxj5/OBiUD7aMXr3dHNnwM9bPh7G/Z9WPclNfZlmdWvp7ajddV4MWU3EVkiogsF5E6Ebm7neOuEhEjIuPcCzGNldTaVrbbpYLtdS0jRJR1vI/Dg09KHZn9OyjoYVvZXjj9G3bHqY8e9Ob8Kit1mNxFJAd4ELgYGA5cLyLDoxzXHbgTmOV2kGmrpBYO73Z3AbHG/bBvkw6DbK2w0i6lm+qW+/4GWPIXm9i9Wk++SzGMvgYWvwCHdntzDZV1Ymm5jwfqjDGrjTFHgKeBy6Mc96/AL4DDLsaX3rwYMeOcSztTTxTKgZ41qU/u856AY0ds69pLY2+EpkOw6Dlvr6OyRizJvS+wIeLr+vBzx4nIGKDSGPOSi7GlPy/qwM65dBjkyUpqU1uWaW6GOX+wm2yUelwmqxgD5aNh7mM6LFK5IpbkHq1X7/i/PhEJAf8FfLfDE4ncIiJzRGROQ4MPox7cVtQfQrnuJvftKwFpGWqpWpQMtMMhm5tTc711H9hZyGO+6v21ROC0G+1wy03zvL+eCrxYkns9UBnxdT9gU8TX3YGRwDsishaYAEyN1qlqjHnIGDPOGDOutDQAw/xyct0vFexYCUVVkNfJvXMGRclAONYIe+tTc72Fz9g12IdemprrjbrarjY674+puZ4KtFiS+2xgkIjUiEg+cB0w1fmmMWaPMaaXMabaGFMNzAQuM8bM8STidON2qWD7Sq23t8XLWcGtHT0ES/4Kwy6Lf+XHRHUqhCFTbAduuixvrDJWh8ndGNME3A68BiwFnjXGfCoi94nIZV4HmPZKBtpFvtwoFRhj3yi03h5dKodDLn/FTpo6xeVJSx0ZeRUc3AGr303tdVXgxLTNnjFmGjCt1XP3tHHsecmHlUFKau0iX3s3QlFlx8e3Z+8mOHqgZVEydaJuvW2ZJBUt94XPQvcKqD7X+2tFGnShnZ28+HkYNDm111aBojNUk+VmqWBnuEVaPDD5cwWRhDuavU7uh/fCqjdhxBfdW0cmVrkFMOzzsPQlWxpSKkGa3JPlanIPLx9cosm9TakYDrnydTu2fbhPVcdRV8KRfVA33Z/rq0DQ5J6s7uWQ19WdhLNztZ2F2aNvx8dmq5Ja2L0Omo54d42lU20JqN94767RnuqJ0KkIlr3sz/VVIGhyT5aIbWm71XLvWZ36UkAmKakF0wy71npz/iMHYeUbMPRzdq9cP+TkwuApsOJVONbkTwwq42lyd4Nbq0PuXKuTlzri9eqQq96Cowdt3dtPQy+FQ7tg/Uf+xqEyliZ3NxwvFTQmfg5jbMtdk3v7nNdnp0d192Uv25JI9TnenD9WAy+AnAJYPq3jY5WKQpO7G0oGhksF6xI/x/5tdhikJvf2dSm2G1R70XJvbrZ7mtZOcn/d9ngVdIOB58Oyl3StGZUQTe5u6FljH3etSfwczkiZ4prk4wk6r0bMbFlgd3qqvdD9cydiyCWwez1s0008VPw0ubvBScg73Uju2nLvUPEAb5L7yvDQw9pJ7p87EYPCbzI6JFIlQJO7G7qW2uGQybbcJcduSqHa17PGbmji9iSfujegz6nQrczd8yaqRwWUjYC6N/2ORGUgTe5ucGZOJttyL6ryv9abCZxPSrvXu3fOgzvtPqmD0qQk46idZEfMNO73OxKVYTS5u6W4OvmWu5ZkYtPThTJYa6vftp3igy5y75xuqJ1kZ8uufd/vSFSG0eTulp41drRMIqtDGmMTlSb32BS70IHd2srp0Lkn9D3NvXO6oepMu3m21t1VnDS5u6W4xm4ksW9Tx8e2dnAnNO7R5B6rLiWQ3929lrsxNnkOvCD9ZgfnFtht/jS5qzhpcndLMqUCHSkTHxG7TINbLfdtS+HANpvc09GA8+29utnHoAJPk7tbkikVaHKPX3G1e+vLOPXsVK/dHquacFxr3vM3DpVRNLm7pUc/u1l2Ii33XWsAgZ79XQ8rsJLp42ht7QworErf1790mC1FaaeqioMmd7fk5NqhjE4rPB47V9vx7bkF7scVVMn0cURqbrZJsyZNW+1gV6fsfzasfU+XIlAx0+Tupp41iZdldNmB+Lg1HHLbErv6YrqWZBw1E2HPBrtAnVIx0OTupuIau2xvvK0rHeMev57V9jHZTtW14Tq236tAdsSJT+vuKkaa3N3Us8YOaTy0K/afObTb7navyT0+hZWJ93FEWvu+faNIdnNzr5UOhS69Wt6MlOqAJnc3JbKAmNPy1LJMfHJybYJPZsSMU29P91Y72OGf1efYeLXurmKgyd1NiSz96yQnp8ygYlecYB+HY+tiOLzb7lmaCWrOhb0bE+u0V1lHk7ubnAQdV8s93EFWlKbD8NJZz5rkyjLHx7dnQMsdWjp9dUikioEmdzfld4Fu5fG1Jnevg87F0KmHd3EFVc9q2/KOp48j0oZZdnx7YV9Xw/JMr8H238qGj/2ORGUATe5uK46zNblrXfpOnkl3yWySYoxN7pXj3Y3JSyI23g2z/Ln+kQN2qYYdq6DpiD8xqJjl+h1A4BQPgFVvxX78rrVQPsqzcAItso+j79j4fnZPPezbDJVnuB+XlyrHw4pX4cAO6FqSmmtumA0z7odVb0Jzk30uvzuMuBw+8//s5D2VdjS5u61njU0aRw9BXuf2j21uthNThn0uNbEFzfGx7mvj/1mn9ZtJLXdoeTOqnw1Dpnh7rWNH4fUfwazfQLfeMOHvoc8pdn35dR/Aohdg4XNw6QMw9qvexqLipsndbZEJp2xY+8fu22z/o2hnamIKukHXssTKMvWz7TrpvUe6H5eXKsba7RjrP/Y2uTc1wtNfsksNj78VJt1jX2/HqV+C874Pf70dpt5h/71f8CNbOlJpQWvubnPq57Esz+pMJddhkInrWZ14y73vaXa8fCbJ7wJ9Rnvbqdp8DJ7/uk3sn/tvuOQXJyZ2R2E/uOEFGHsjvPcf8P5/eheTipsmd7c5rfBdMawBomPckxdvBzbYjsHNCzOv3u6oPAM2fmLLJl5452ew7CWY8nMY97X2jw3l2DeAUVfDm/fB8le9iUnFLabkLiJTRGS5iNSJyN1Rvv8dEVkiIgtF5E0Ryd46Q7cyyO0U2wJPu9YBYltAKjFF/e3KkPEkuk3zwBzL3OTe73Q4etBOwnLb6nds5+mYG+CMW2P7mVAILvuVrcf/+RbYu9n9uFTcOkzuIpIDPAhcDAwHrheR4a0OmweMM8aMBp4HfuF2oBlDxI4eiCW5714HPSp0qd9kFFXZja331Mf+M05nar9x3sTkNeee0m43AAAS2ElEQVRNye3STOM++OsdUFILlzwQX/08rzNc9Qc7RPJv39IlEtJALC338UCdMWa1MeYI8DRweeQBxpi3jTEHw1/OBLK7KVpUFVvNfdc67UxNVjx9HI4NH9sJQV2KvYnJa4X9oHuF++Pd3/l3O3rr8l93PNIrmpKBtuN15Wuw7GV3Y1NxiyW59wU2RHxdH36uLTcDr0T7hojcIiJzRGROQ0ND7FFmmqL+sdXcd6/TenuynDHWsa5zbowdKdMvw4ZARhKBytPtfbhl5xr4+CE49ctQlUS5avwtdgXLN36kE518Fktyj/bZLOpnLhG5ARgH3B/t+8aYh4wx44wx40pLS2OPMtMUVdlp8Yf3tH1MUyPs3aSzU5PVoy9IKPaW++71donlvmO8jctrfU+z93Jghzvnm36vXUL5gh8md56cXLjop3Zxs9kPuxKaSkwsyb0eiFzsuh9w0t5mIjIZ+CFwmTGm0Z3wMlQspYLdGwCjZZlk5eTZ/WtjTe6b5tnHijhntKabivCbk3M/ydj4CSz5C5x1h+0DStagyVA7Gd79ha3jK1/EktxnA4NEpEZE8oHrgKmRB4jIGOD/sIl9m/thZpjjpYL2kvta+6gt9+TF2scBNhmG8qD3CG9j8lqfU+2jG8n9/f+CToU2ubvl/B/YT69z/uDeOVVcOkzuxpgm4HbgNWAp8Kwx5lMRuU9ELgsfdj/QDXhOROaLyNQ2TpcdiqrtY3t19106gck1RVWx9XEAbJoL5SMzf4RSpx5QMij55N6wApa+ZGvlBd3diQ1s2ajmM/DRg7YEqVIupul5xphpwLRWz90T8ffJLseV2boUQ3639luTu9ZCToFdIlglp2d/u5RDU2P7Sbu5GTYtgFFXpi42L1WMSX7bvQ//x87LOOM2d2KKdM4/whNfgAVPwWk3uX9+1S6doeqFWMa6715n9+0M6a8gaUVVgOl4rPvO1XaP24oM70x19B1r39QSnTS0bysseMZOWOray93YAAacZ1/rD/9Xx737QDOLV4r6d9By1zHurol1OGRQOlMdzpvU5vmJ/fz8P0Lz0dhnosZLxJZ7dqzU3aN8oMndK04duK0Wy27dpMM1sXRgg62353a247CDoHy0HQa6cW78P9vcDJ88arfu6zXI9dCOG/FF21n7iXasppomd6/07A9H9kXfAu7wHvu8dqa6o3uFHaPdYXKfZ1dUzLSVINuS3wVKhyXWqbrqLft6dbQwWLLyOtuJUUumwv4AT1xMQ5rcvdJeqUA3xXZXTq6dzNTeiJnmY7B5QXDq7Y6+Y8ILocVZ0/7kD9ClFwz9vDdxRTrta7b8M/+P3l9LHafJ3StF7Uxkcp7T7cnc09FY94bldiXFoCX3ijFwcLtdEyZWB3bYrfpOvR5y872LzVE6GConwIKntWM1hTS5e8VJ3NFak85/RE3u7unZQQf25gX20Zn8ExTHZ6rG0am65M92L9TR13oTUzSjr4aGZbBlUequmeU0uXulc5HtSIract9gO/a6pGiD42xQ1B/2b4Gjh6N/f8si+5p72Xnoh7LhtlM1nrXdFz1va/Wp3GJw+Bdtv8iiZ1N3zSynyd1LbY1137PBjnHX/Sbd43wKaqs8sWUh9B5udw4KkrzOdvniWFvEu9fD+o9g1FWp/ffXtcSuN7PoBdv/oTynyd1LbS39u2cDFFae/LxK3PE+jiivtzE2+ZWPSm1MqVI+Kvbkvug5+zjqau/iacvoa+yuWes+SP21s1BAxoSlqaL+UPemTS6RraTdG+yWZMo97Y1131NvF7EKcnJf9Bwc3NnxBiSf/tmuZe/HHIvBF0NeF1jyV6iZmPz5GvfDnEfs+XastAuRlw2FEVfA2K9Aftfkr5HBtOXupZ79oekQHIgY33vkoB3doC13d3Uvt6s9Rvuk5LRqy0enNqZUcd60Oqq771prX4thKRj+GE1+F6idZHdpam5O7lx10+FXY+2mIKbZdg6PvgaaDsOr/w8ePMPuB5vFNLl7KVpr0ln/REfKuCuUY/sxorXctywCxHY+BlHvcHLvqDTjbH037HPextOeoZ+36+Fs/CTxc8z6P/jjVXac/s1vwC1vwyX3w6UPwK0z4Guv2E8IT1xhZ+FmKU3uXnJa55GdfHvWn/g95Z62xrpvWWj39yzolvqYUqFbqV1dtKPkvvQlKBsBxQNSE1c0gy+yo2aW/S2xn5/9O3jln2DopfCN6VAZZbvE/mfBN9+EgRfYzbrnP5lczBlKk7uXCsP7hEeuVrjbGeOuyd11bY1OCnJnqqOjTtX9DXaUjJ+tdoDOPe16Nktfin9C06q3Ydr3bO3+6kdtmactBd3h2j/alSmn3gHrZyYRdGbS5O6lzkVQ0KMloYNtxUsOdO/jX1xBVVhl+zeOHmp57tBum/CzIbk3LGt7Y4wVrwAGhvqc3MHW/HeusvHGat9WeP7rdtG3Kx+22yt2JK8TXP2Y/ZT83E3u7TebITS5e62w38kt9x59gzfeOh04n4b2bGx5buun9jGonamO8lF21mnD8ujfX/Ga3Ws2Hd7kBk+xjytfj+14Y2x55ehB22KPZ8eozkVwzeNwYDu8enfcoWYyTe5eK6xsqbNDywQm5b7jZbCIT0rOCJJ0SGpect68opVmmo7YkSODJqfHxLnCvrb2v/KN2I5f9pL95HHBj6B0SPzX6zMaJt5lZ8fGes0A0OTutWgtd+1M9Ua0DuwtC6FrKXTr7U9MqVJcA3ldoyf3DTPhyH6ovTD1cbVl0GRbBz+8t/3jmhrh9X+2yyUksxXgOd+Bklp47YdwrCnx82QQTe5eK6q0a7c37rf/qPZt0pa7V3pU2HVWIt9MtyyG3iPSo8XqpVCOXV4hWnJf+YadAzDgM6mPqy21F9plgNe82/5xs35rx+dP+bfk1uHPzYfJP4bty2He44mfJ4Nocvfa8dZkvU3spllb7l7JybMd1U4HdvMxW4MuG+FvXKlSNgy2LTl5FErdm1A1Ib5atdeqJkB+9/bLJId2w4wHYNBn7bDGZA29FKrOhLf/DRr3JX++NKfJ3WuRpQIdBum9wn4tZZlda+0M4bJhvoaUMmXD4dDOE2dE79kI2z6FQWlUkgH7RjzgM3amaVtDImf/Dhr3wgU/dOeaInDhffb1+eSxxM5x7Cjs33biiKw0pcnda5GdfE7SKdTZqZ4prGx5nbcttY9BnZnamvMmtm1Jy3Or3rKPtZNTH09HaifD3o2wfeXJ3ztyEGb+GgZd5O46TJXj7Tj7jx60Hc2xWv2OnRX7bxXwwCD4aR/47Tnw8cM24achTe5e614e3t8zouXuJHzlvsJ+trXa3NyS3BMZYZGJnDcx574B1r5nO5TT8Q3O6QOIVnef+zgc3AHnftf96579bVsidVbIbE/jPnjua/D45fZNc9zNcMkDcN7d9v/1tLvg/ybCtjjG7KeIrgrptVCO7ejbU28nVXQts4/KG0WVtqNu/1b7n7FndXCXHWitW5ldb8VpuRsDa2bYlmo6dij3rLGftNa8C+O/2fJ88zGY+aCtj1dNcP+6tZPsejwf/A+c+qW2X5sD21uS+vn/DGffCbkFLd8/725YNs2Owf/9RXD9U1B9tvvxJkhb7qlQWNVSc9d6u7ciO7C3LU3PFquXyoa1tNx3rLKLdNWc629MbRGBms/AmvdOXCWybrpdI2j8Ld5d96zb7ciZte9FP6Zxn03sO+rgy8/DZ753YmJ3DL0EvvmW/YT+5LV2s/I0ock9FZyx7rpJh/ec13fnarvGd7Z0pjrKhtvkbgysnWGfq3Zh7XSv1Ey0a+1vjRjCOecR+wnXy6UShl8OnYpgzh9O/l5zM7x4q30dr/uTbem3p6gSvvoXu27OU1+y6/ikAU3uqVBUaTuOdq3VlrvXnP6M1W/b6fjZ2HI/st82JNbMgO4VdkXMdOV8qlgTfiPavd4ulTD2q3ZsulfyOtuSzNK/2fJLpNkPw/KX4bM/jb0jukeFfSM4uANe/Gb8i6J5QJN7KhT2s+Pbm5t0pIzXOvWwG5OveM1+nY0td4CtS2y5oyZN6+2OHhVQMghWhztV5z5u4z3tJu+vfdpNtn9m/p9antu5BqbfCwMnxT8jts9omPIz27CY28ZQyz318OxX7eQ6j2lyT4XIUoy23L1XWGnHe4dybeLIJmVD7eOnL9odv6rTtN4eqWYirPvQDk1c8AwMOD81/09Kh0DVWXbMuzHhBcrutKu2XvbLxN4Ux33dvuav/wj2bj7xe4tfhP8dDyteP3FEk0diSu4iMkVElotInYictLSaiBSIyDPh788SkWq3A81okcm9OI0/IgeFU4YoqfX2o3066lRoV39c+Iz92o29Sr024DNw9IAdIbNnPZxyXequfeqX7PLDG+fanarWzIAL7018uLKIfWNoaoTXftDy/MzfwPNfg/KR8A8zYbT3G5R3mNxFJAd4ELgYGA5cLyKtC5k3A7uMMbXAfwE/dzvQjBb5D6VXlrUk/eAsN1DkwybQ6cApRRVV+bMRdrycTxfv/NxujzfkktRde/hlkFMAC56EN38MvQbD2JuSO2fxALsK5acv2qUf5j9plxse+jm48W92eG4KxNJyHw/UGWNWG2OOAE8Dl7c65nLAKTI9D0wSSedCX4o5O8bkdU3v+mdQOP95JEurjl172Udn3fR016XYLsncdMgmwFTOS+hUCEOm2KUOtq+ASf+S3AJljrO/ZT85Tr3DjoOvmQhXPRJ9OKVHYvnX3xeIWEOV+vBzUY8xxjQBe4ASNwIMjNvehzvTZwxsoPU/0z6e/g1/4/CLM8IjFZ2SbqkJz1YdfW3qr33mHfaxcoJdXMwNuQVw6X/YUXJdy+CqR1Oa2CG2GarRmpqtx/nEcgwicgtwC0BVVZaNGgn6ZhHppKgK7t3jdxT+GXmlXZOlUw+/I4nduK+HFxM7L/XXrjwd7qqzydfNT9YDzrOt9fLR0DX1bd1YWu71QGTXdT9gU1vHiEguUAjsbH0iY8xDxphxxphxpaWliUWslGqfSGYldrCd4JPvdackkohupd68ZiOv9K2fLZbkPhsYJCI1IpIPXAdMbXXMVODG8N+vAt4yJg1G8SulVJbq8G3SGNMkIrcDrwE5wCPGmE9F5D5gjjFmKvB74AkRqcO22FM4lkkppVRrMX0GMsZMA6a1eu6eiL8fBrwfuKmUUiomWTpWTCmlgk2Tu1JKBZAmd6WUCiBN7kopFUCa3JVSKoDEr+HoItIArEvwx3sB2zs8KjPovaSfoNwHBOdegnIfkPy99DfGdDgL1LfkngwRmWOMGed3HG7Qe0k/QbkPCM69BOU+IHX3omUZpZQKIE3uSikVQJma3B/yOwAX6b2kn6DcBwTnXoJyH5Cie8nImrtSSqn2ZWrLXSmlVDvSOrkHaWPuGO7lJhFpEJH54T9puY2QiDwiIttEZHEb3xcR+WX4PheKyNhUxxiLGO7jPBHZE/H7uCfacelARCpF5G0RWSoin4rIt6Ick/a/lxjvIyN+LyLSSUQ+FpEF4Xv5cZRjvM1fxpi0/INdXngVMADIBxYAw1sd8/fAb8N/vw54xu+4k7iXm4D/9TvWGO5lIjAWWNzG9y8BXsHuzjUBmOV3zAnex3nAS37HGeO99AHGhv/eHVgR5d9X2v9eYryPjPi9hF/nbuG/5wGzgAmtjvE0f6Vzyz1IG3PHci8ZwRgzgyi7bEW4HHjcWDOBIhHpk5roYhfDfWQMY8xmY8zc8N/3AUs5eZ/jtP+9xHgfGSH8Ou8Pf5kX/tO6g9PT/JXOyT1IG3PHci8AV4Y/Mj8vIpVRvp8JYr3XTHBm+GP1KyIywu9gYhH+aD8G21KMlFG/l3buAzLk9yIiOSIyH9gGvGGMafN34kX+Sufk7trG3Gkgljj/BlQbY0YD02l5R880mfI76chc7DTvU4BfAX/xOZ4OiUg34AXg28aYva2/HeVH0vL30sF9ZMzvxRhzzBhzKnbf6fEiMrLVIZ7+TtI5ubu2MXca6PBejDE7jDGN4S8fBk5LUWxui+X3lvaMMXudj9XG7kSWJyK9fA6rTSKSh02IfzLGvBjlkIz4vXR0H5n2ewEwxuwG3gGmtPqWp/krnZN7kDbm7vBeWtU/L8PWGzPRVOCr4dEZE4A9xpjNfgcVLxEpd+qfIjIe+39lh79RRReO8/fAUmPMf7ZxWNr/XmK5j0z5vYhIqYgUhf/eGZgMLGt1mKf5K6Y9VP1gArQxd4z3cqeIXAY0Ye/lJt8CboeIPIUdsdBLROqBf8F2FmGM+S12r91LgDrgIPA1fyJtXwz3cRXwdyLSBBwCrkvThgPA2cBXgEXhGi/AD4AqyKjfSyz3kSm/lz7AYyKSg30DetYY81Iq85fOUFVKqQBK57KMUkqpBGlyV0qpANLkrpRSAaTJXSmlAkiTu1JKBZAmd6WUCiBN7kopFUCa3JVSKoD+P/VBXtsxMhq1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1a3e6924d0>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.figure(1)\n",
    "plt.plot(t_res,e_res,t_res,h_res)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}