03_NonUniformFlow03/Q2DNonUniformFlow-text-3_code_CalSectprof.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3d68231e",
   "metadata": {},
   "source": [
    "# 分割断面法のプログラミング"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51c9f370",
   "metadata": {},
   "source": [
    " - 前回定義したクラスsectionを利用する。\n",
    " - 分割断面法は、各断面を複数の領域に分割するため、クラスsectionをクラスsubsectionと再定義して、複数のsubsectionの集合が新たなクラスsectionとなるように定義する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "62bc4703",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c0d1c066",
   "metadata": {},
   "outputs": [],
   "source": [
    "class subsection(object):\n",
    "    def __init__(self, l, z, n):\n",
    "        self.l = l\n",
    "        self.z = z\n",
    "        self.n = n\n",
    "\n",
    "    def H2ABSKn(self, H): #変更:引数は水位のみ\n",
    "        A, B, S, SN = float(0), float(0), float(0), float(0)\n",
    "        l = self.l #変更:メンバ変数を参照\n",
    "        z = self.z #変更:メンバ変数を参照\n",
    "        n = self.n #変更:メンバ変数を参照\n",
    "        \n",
    "        for i in range(1, len(l)):\n",
    "            dx = l[i] - l[i-1]    \n",
    "            dy = z[i] - z[i-1]    \n",
    "            hb, hf = H - z[i-1], H - z[i]\n",
    "            \n",
    "            if hb <= float(0) :\n",
    "                if hf > float(0) :\n",
    "                    dx_dh = dx / (hf - hb)\n",
    "                    B += hf * dx_dh\n",
    "                    A += 0.5 * hf * hf * dx_dh\n",
    "                    Sp = hf * np.sqrt( dx_dh * dx_dh + 1.0)\n",
    "                    S +=  Sp\n",
    "                    SN += Sp * n[i-1]**1.5\n",
    "            elif hf <= float(0) :\n",
    "                if hb > float(0) :\n",
    "                    dx_dh = dx / (hf - hb)\n",
    "                    B -= hb * dx_dh\n",
    "                    A -= 0.5 * hb * hb * dx_dh\n",
    "                    Sp = hb * np.sqrt(dx_dh * dx_dh + 1.0)\n",
    "                    S += Sp\n",
    "                    SN += Sp * n[i-1]**1.5\n",
    "            else :\n",
    "                B += dx\n",
    "                A += 0.5 * dx * (hf + hb)\n",
    "                Sp = np.sqrt(dx**2 + dy**2)\n",
    "                S += Sp\n",
    "                SN += Sp * n[i-1]**1.5\n",
    "                \n",
    "        if S <= float(0):\n",
    "            nd = float(0)\n",
    "            K = float(0)\n",
    "        else:\n",
    "            nd = (SN/S)**(2.0/3.0)\n",
    "            K = A**(5.0/3.0)/nd/S**(2.0/3.0)\n",
    "            \n",
    "        return A, B, S, K, nd"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ffb6f4cd",
   "metadata": {},
   "source": [
    "下図のように分割断面を設定する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "20b605f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "l=np.array([0, 5,93,100,200,206,294,300])\n",
    "z=np.array([6, 3.5, 3.5,  0,  0,  3,  3,  6])\n",
    "n=np.array([0.041, 0.041,0.030,0.030,0.030,0.040,0.040])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7e356726",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhYAAAGdCAYAAABO2DpVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAj+klEQVR4nO3dfXBU5dnH8d/mbRPzVkAhSYmUp059C1AFq8GXR1GZZorV2jrqWEvfnNICI8XWGv0DaKeGdqaOPkOlSh2rMypOR7HOqFgcJegwtEBxoFYtVlqigIyo2Wwgm5Dczx/NrhvIQjbsOec+9/l+ZnaGbDbk8vbMnh/XfZ2zMWOMEQAAQAEUBV0AAABwB8ECAAAUDMECAAAUDMECAAAUDMECAAAUDMECAAAUDMECAAAUDMECAAAUTInfv3BgYEB79uxRdXW1YrGY378eAACMgjFGXV1damhoUFFR7r6E78Fiz549amxs9PvXAgCAAujo6NDEiRNzft/3YFFdXS3pv4XV1NT4/esBAMAoJBIJNTY2Zs7jufgeLNLbHzU1NQQLAABC5nhjDAxvAgCAgiFYAACAgiFYAACAgiFYAACAgiFYAACAgiFYAACAgiFYAACAgiFYAACAgiFYAACAgsk7WLz//vv65je/qXHjxumkk07SF7/4RW3dutWL2gAAQMjkdUvvjz/+WBdeeKEuu+wyvfDCCxo/frz+9a9/6TOf+YxH5QEAgDDJK1j86le/UmNjox5++OHMc5/73OcKXRMAAAipvLZCnn32Wc2YMUPXXXedxo8fr3POOUerVq065s+kUiklEokhDy8sWbJEP/rRj/TBBx948vcDAGC7W2+9VT/96U/10UcfBVZDzBhjRvri8vJySdLixYt13XXX6a9//asWLVqkBx54QN/61reG/ZmlS5dq2bJlRz3f2dlZuE837e6WqqokSds3btTU5ubC/L1wU9bxomRSqqwMth5gOBynyNNAV5eKBs+rH/zrX5rwP/9T0L8/kUiotrb2uOfvvDoWAwMDOvfcc3X33XfrnHPO0Q9+8APdcsstWrlyZc6faW1tVWdnZ+bR0dGRz6/Mm1cdEQAAbNbd3Z35c3V1dWB15BUs6uvrddZZZw157swzz9Tu3btz/kw8HldNTc2Qh5eSyaSnfz8AADbKPv9VVFQEVkdeweLCCy/U22+/PeS5f/7zn5o0aVJBizoRdCwAAFGUff6LxWKB1ZFXsPjxj3+sTZs26e6779Y777yjxx9/XA8++KDmz5/vVX156+rqCroEAAB8Z8v5L69gcd5552nNmjV64okn1NTUpF/84he69957ddNNN3lVX95sWVgAAPxkS8c+r/tYSNKcOXM0Z84cL2opCIIFACCKbJkxdO6zQggWAIAosuX851ywsKUVBACAn2w5/zkXLGxJbAAA+MmW8x/BAgAAB9hy/iNYAADgAFvOfwQLAAAcYMv5z7lgYcvwCgAAfrLl/OdcsLDlOl4AAPxky/nPuWCRSCSUxyfBAwDgBDoWHjEa+tGxAABEATMWHrIltQEA4BeChYdsWVwAAPxiy7nPyWBBxwIAECV9fX061NMTdBmSCBYAAISeLd0KydFgYdMCAwDgNZv+Qe1ksLBpgQEA8JpN/6AmWAAAEHI2nfcIFgAAhJxN5z0ng4VNLSEAALxm03nPyWBhU3IDAMBrNp33CBYAAIScTec9ggUAACFn03nPyWBh014TAABes+m852SwsCm5AQDgNZvOewQLAABCzqbzHsECAICQs+m852SwsGmvCQAAr9l03nMyWBw8eFCHDx8OugwAAHxBx8IHNqU3AAC8RLDwULysTJJdiwwAgJdsOuc5Fyxqamok0bEAAESDMcaqc55zwaKqqkqSXekNAACv9PT0WDVX6FywqK6ulkSwAABEg23nO+eCBVshAIAoSQeL6sGOfdCcCxZ0LAAAUZL+h3T6/Bc0ggUAACGW6VgQLLxBsAAARAnBwmPphWXGAgAQBelgkZ4xDJpzwSK9sHQsAABRwIyFx7iPBQAgStLnuyquCvEGHQsAQJSwFeIxZiwAAFHCVojHuCoEABAldCw8RrAAAEQJMxYeI1gAAKIk1PexWLp0qWKx2JBHXV2dV7WNSvZnhRhjAq4GAABvpWcsbNkKKcn3B84++2y99NJLma+Li4sLWtCJSreC+vr6lEqlVF5eHnBFAAB4x7aORd7BoqSkxLouRbbsPaZEIkGwAAA4zbZgkfeMxc6dO9XQ0KDJkyfrhhtu0LvvvnvM16dSKSUSiSEPLxUXF3OTLABAZIQ6WJx//vl69NFH9eKLL2rVqlXat2+fZs6cqQMHDuT8mba2NtXW1mYejY2NJ1z08WTPWQAA4KqBgQElk0lJIQ0WLS0t+vrXv64pU6boiiuu0HPPPSdJeuSRR3L+TGtrqzo7OzOPjo6OE6t4BLgyBAAQBelQIYV4eDNbZWWlpkyZop07d+Z8TTweVzweP5Ffkzdu6w0AiIL0ea6kpMT3c20uJ3Qfi1QqpTfffFP19fWFqqcg2AoBAERB9l03Y7FYwNX8V17B4ic/+Yna29u1a9cu/eUvf9E3vvENJRIJzZ0716v6RoWOBQAgCmy7h4WU51bIe++9pxtvvFEffvihTjnlFF1wwQXatGmTJk2a5FV9o8KMBQAgCmy7IkTKM1isXr3aqzoKio4FACAKbPsAMsnBzwqRmLEAAEQDwcIndCwAAFFg44yFk8GCGQsAQBTYOGPhZLCgYwEAiAK2QnzCjAUAIAoIFj6hYwEAiAJmLHzCjAUAIAqYsfAJHQsAQBSwFeKT9AInk0kNDAwEXA0AAN5gK8Qn6QU2xqi7uzvgagAA8AYdC5+Ul5eruLhYEtshAAB3MWPhk1gsxpwFAMB5dCx8xL0sAAAu6+3tVSqVkkSw8AWXnAIAXJb9D2e2QnzAVggAwGXp81tFRYVKSkoCruZTBAsAAELIxvkKKQLBghkLAICLbLyHheRwsGDGAgDgMjoWPmMrBADgMhvvYSFFIFiwFQIAcBEdC5/RsQAAuIwZC58xYwEAcBkdC5/RsQAAuIwZC58xYwEAcBkdC5/RsQAAuIwZC58xYwEAcBkdC5/RsQAAuIwZC5+lg0VPT4/6+voCrgYAgMKiY+Gz7ATHACcAwDXMWPistLRU5eXlktgOAQC4h62QADBnAQBwkTGGrZAgcC8LAICLDh06pIGBAUkEC1/RsQAAuCh9XovFYqqsrAy4mqGcDhbcywIA4KLs+YpYLBZwNUM5HSzoWAAAXGTrfIUUkWDBjAUAwCW2XmoqRSRY0LEAALiEjkVAmLEAALjI1ntYSI4HC7ZCAAAuomMRELZCAAAuYsYiIGyFAABcRMciIHQsAAAuYsYiIMxYAABcRMciIHQsAAAuYsYiIMxYAABc5GzHoq2tTbFYTIsWLSpQOYWV3bEwxgRcDQAAheHkjMXmzZv14IMPaurUqYWsp6DSwaK/v189PT0BVwMAQGE417FIJpO66aabtGrVKo0ZM6bQNRVMZWVl5lPf2A4BALjCuRmL+fPn6ytf+YquuOKK4742lUopkUgMefilqKhIVVVVkggWAAB32NyxKMn3B1avXq2tW7dqy5YtI3p9W1ubli1blndhhVJTU6Ouri6CBQDACf39/eru7pbkwIxFR0eHbr31Vj322GMqLy8f0c+0traqs7Mz8+jo6BhVoaPFvSwAAC7JPp+FvmOxdetW7d+/X9OnT88819/frw0bNmjFihVKpVIqLi4e8jPxeFzxeLww1Y4C97IAALgkHSzKysoCPb/mklewuPzyy7Vjx44hz33nO9/RGWecoZ/97GdHhQobcC8LAIBLbJ6vkPIMFtXV1WpqahryXGVlpcaNG3fU87agYwEAcInN97CQHL/zpsSMBQDALU51LIazfv36ApThHToWAACX2HwPCykCHQtmLAAALrG9Y+F8sGArBADgEmYsAsZWCADAJWyFBIxgAQBwCVshAWPGAgDgErZCAsaMBQDAJXQsAsZWCADAJcxYBIxgAQBwCR2LgKX3oJLJpAYGBgKuBgCAE8OMRcCyE10ymQywEgAAThwdi4DF43GVlpZKYjsEABB+zFgELBaLcckpAMAZdCwswAAnAMAFqVRKvb29kpixCBT3sgAAuCD7H8gEiwDRsQAAuCD9D+TKykoVFxcHXM3wIhEsmLEAALjA9vkKKSLBgo4FAMAFtt/DQopYsGDGAgAQZnQsLEHHAgDgAtvvYSFFJFgwYwEAcEEYOhYlQRfgh/T/gPb2di1atCjYYhC4sr4+/TroIgBgFMIwYxGJYFFfXy9Jeuutt/TWW28FXA2CdpJEsAAQSnQsLPG1r31N9957r/bv3x90KbCASSal//s/SVJvb6/KKisDrggARiYMMxaRCBbxeFy33npr0GXAEoc7OzPBoqurS+PGjAm4IgAYmTB0LCIxvAlkKyn5NE9zCTKAMAnDjAXBApFGsAAQJmHYCiFYINIIFgDChK0QwHIECwBhQrAALMdN0wCECTMWgOXoWAAIE2YsAMsRLACEhTGGrRDAdgQLAGHR3d0tY4wkggVgLYIFgLBIdyuKiopUUVERcDW5ESwQaQxvAgiL7PmKWCwWcDW5ESwQaclkMugSAGBEwjBfIREsEHF0LACERRguNZUIFog4ZiwAhAUdCyAECBYAwiIM97CQCBaIOIIFgLCgYwGEAMECQFgwYwGEQCKRyNxwBgBsRscCCIEBY3Tw4MGgywCA42LGAggJtkMAhAEdCyAkuJcFgDBgxgIICYIFgDBwsmOxcuVKTZ06VTU1NaqpqVFzc7NeeOEFr2oDfEGwABAGTs5YTJw4UcuXL9eWLVu0ZcsWzZo1S1dffbXeeOMNr+oDPMeMBYAwCEvHoiSfF1911VVDvv7lL3+plStXatOmTTr77LMLWhjgFzoWAMIgLDMWeQWLbP39/frjH/+o7u5uNTc353xdKpVSKpXKfM2bOGzDMQkbDQwMZFrKt99+u3pLSwOtB8H76KOPJDnWsZCkHTt2qLm5WT09PaqqqtKaNWt01lln5Xx9W1ubli1bdkJFAl5iKwQ22rhxoy4a/PNv779f3G0FklRSUqJx48YFXcYx5R0sTj/9dL3++uv65JNP9NRTT2nu3Llqb2/PGS5aW1u1ePHizNeJREKNjY2jrxgoMDoWsNH+/fszf/7pT36ivrKyAKuBLS644AL3OhZlZWU67bTTJEkzZszQ5s2bdd999+mBBx4Y9vXxeFzxePzEqgQ8RLCAjbI7aUuXLpUqK4MrBsjDCd/HwhgzZIYCCBuCBWzEcYmwyqtjceedd6qlpUWNjY3q6urS6tWrtX79eq1du9ar+gDPMWMBG3FcIqzyChYffPCBbr75Zu3du1e1tbWaOnWq1q5dqyuvvNKr+gDP8S9D2CiZTAZdAjAqeQWLhx56yKs6gMAQLGAjjkuEFZ8VgsjjDRw2YisEYUWwQOTxBg4bEXgRVgQLRB5v4LARMxYIK4IFIq+7u1v9/f1BlwEMQScNYUWwAMSbOOxDJw1hRbBApJUNfrATwQK24ZhEWBEsEGnpe+7zr0PYxBhDsEBoESwQadXV1ZIIFrBLKpVS3+HDQZcBjArBApFWVVUliWABu3A8IswIFoi09FYIbWfYhGCBMCNYINLYCoGNCLoIM4IFIo3hTdiI4xFhRrBApDFjARtxPCLMCBaItPRWCK1n2IRggTAjWCDS2AqBjQi6CDOCBSKNrRDYiOMRYUawQKTRsYCNOB4RZgQLRBozFrARxyPCjGCBSKNjARtxPCLMCBaINGYsYCOOR4QZwQKRxp03YSOOR4QZwQKRxmeFwEYcjwgzggUiLd2x6O3tVSqVCrga4L/oWCDMCBaItPSMhcSbOezBsYgwI1gg0oqLi1VZWSmJN3PYg2MRYUawQOQxZwGbDAwMKJlMBl0GMGoEC0Qe97KATbq7u2WMCboMYNQIFog8LjmFTdLHYUlxccCVAKNDsEDksRUCm6SDRTrwAmFDsEDksRUCm6QDLsECYUWwQOQRLGCT9HGYPi6BsCFYIPKYsYBN2ApB2BEsEHnMWMAm6WCRffM2IEwIFog8tkJgk3TAZSsEYUWwQOQRLGATtkIQdgQLRB4zFrAJwQJhR7BA5DFjAZsQLBB2BAtEHlshsAn3sUDYESwQeQQL2IT7WCDsCBaIPGYsYBMuN0XYESwQedkzFnyqJIJGxwJhR7BA5KXfwI0x6u7uDrgaRB0zFgg7ggUir6KiQsWDH1HNdgiCxlUhCDuCBSIvFosxZwFrECwQdgQLQNzLAnbo6+tTT0+PJIIFwiuvYNHW1qbzzjtP1dXVGj9+vK655hq9/fbbXtUG+IZLTmGD7GBLsEBY5RUs2tvbNX/+fG3atEnr1q3T4cOHNXv2bAbeEHoEC9ggffxVVFSotLQ04GqA0SnJ58Vr164d8vXDDz+s8ePHa+vWrbrkkksKWhjgJ2YsYAPmK+CCvILFkTo7OyVJY8eOzfmaVCqlVCqV+Zo3btiIGQvYgI9MhwtGPbxpjNHixYt10UUXqampKefr2traVFtbm3k0NjaO9lcCnmErBDbg5lhwwaiDxYIFC7R9+3Y98cQTx3xda2urOjs7M4+Ojo7R/krAM2yFwAYEC7hgVFshCxcu1LPPPqsNGzZo4sSJx3xtPB5XPB4fVXGAX+hYwAbMWMAFeQULY4wWLlyoNWvWaP369Zo8ebJXdQG+YsYCNmDGAi7IK1jMnz9fjz/+uP70pz+purpa+/btkyTV1taqoqLCkwIBP9CxgA3YCoEL8pqxWLlypTo7O3XppZeqvr4+83jyySe9qg/wBTMWsAHBAi7IeysEcBEdC9iAGQu4gM8KAcSMBezAjAVcQLAARMcCdmArBC4gWABixgJ2IFjABQQLQJ++kR86dEiHDx8OuBpEFTMWcAHBAtDQN3LmLBAUZizgAoIFIKmsrEzl5eWS2A5BcNgKgQsIFsAg5iwQJGMMwQJOIFgAg7jkFEE6dOiQ+vv7JTFjgXAjWACDuOQUQUoH2lgspsrKyoCrAUaPYAEMIlggSNlXhBQV8daM8OLoBQYxY4EgcakpXEGwAAYxY4EgMbgJVxAsgEFshSBI3MMCriBYAIMIFggSHQu4gmABDGLGAkFixgKuIFgAg5ixQJDoWMAVBAtgEFshCBIzFnAFwQIYRLBAkOhYwBUEC2AQMxYIEjMWcAXBAhjEjAWCxFYIXEGwAAaxFYIgsRUCVxAsgEHZwcIYE3A1iBqCBVxBsAAGpfe2Dx8+rJ6enoCrQdQwYwFXECyAQVVVVZk/M2cBvzFjAVcQLIBBRUVFXBmCwLAVAlcQLIAsBAsEob+/X8lkUhLBAuFHsACycGUIgpAOFRIzFgg/ggWQhXtZIAjp4620tFTxeDzgaoATQ7AAstCxQBCy5ytisVjA1QAnhmABZGHGAkFgcBMuIVgAWehYIAjcwwIuIVgAWZixQBC4hwVcQrAAstCxQBDYCoFLCBZAFmYsEASCBVxCsACy0LFAEJixgEsIFkAWZiwQBGYs4BKCBZCFjgWCwFYIXEKwALIwY4EgsBUClxAsgCxshSAIdCzgEoIFkIWtEASBGQu4hGABZMnuWAwMDARcDaKCjgVcQrAAsmTvcWd/lDXgJWYs4BKCBZClvLxcJSUlkpizgH/oWMAlBAsgSywWY84CvmPGAi4hWABHIFjAT6lUSqlUShLBAm7IO1hs2LBBV111lRoaGhSLxfTMM894UBYQHO5lAT9lb7lVVVUFWAlQGHkHi+7ubk2bNk0rVqzwoh4gcNzLAn5KH2cnnXRSZr4HCLO8j+KWlha1tLR4UQtgBbZC4CcGN+Eaz+Nx9v6hxJs17EewgJ8IFnCN58ObbW1tqq2tzTwaGxu9/pXACWHGAn7iHhZwjefBorW1VZ2dnZlHR0eH178SOCHMWMBPXGoK13i+FRKPxxWPx73+NUDBsBUCP7EVAtdwHwvgCAQL+IlgAdfk3bFIJpN65513Ml/v2rVLr7/+usaOHatTTz21oMUBQWDGAn5ixgKuyTtYbNmyRZdddlnm68WLF0uS5s6dqz/84Q8FKwwICjMW8BMzFnBN3sHi0ksvlTHGi1oAK7AVAj+xFQLXMGMBHIFgAT8RLOAaggVwBGYs4CdmLOAaggVwBGYs4CdmLOAaggVwhPQbfE9Pj3p7ewOuBq5jKwSuIVgAR8huSdO1gNcIFnANwQI4QklJiSoqKiQxZwHvMWMB1xAsgGEwZwE/GGOYsYBzCBbAMLjkFH44ePCgBgYGJBEs4A6CBTAMLjmFH9LHV1FRkU466aSAqwEKg2ABDIOOBfyQPV8Ri8UCrgYoDIIFMAxmLOAH5ivgIoIFMAw6FvADl5rCRQQLYBjMWMAPXGoKFxEsgGHQsYAf6FjARQQLYBjMWMAPzFjARQQLYBh0LOAHOhZwEcECGAYzFvADMxZwEcECGAZbIfADHQu4iGABDIOtEPiBGQu4iGABDINgAT/QsYCLCBbAMJixgB+YsYCLCBbAMLJnLIwxAVcDV7EVAhcRLIBhpN/o+/v7dejQoYCrgavYCoGLCBbAMCorKzOfNsl2CLxCsICLCBbAMGKxGHMW8BwzFnARwQLIgXtZwEv9/f06ePCgJDoWcAvBAsiBS07hpezASscCLiFYADkQLOCl9HEVj8cVj8cDrgYoHIIFkAMzFvAS8xVwFcECyIEZC3iJe1jAVQQLIAe2QuAlLjWFqwgWQA4EC3iJYAFXESyAHJixgJeYsYCrCBZADsxYwEvMWMBVBAsgB7ZC4CW2QuAqggWQA1sh8BLBAq4iWAA50LGAl5ixgKsIFkAOzFjAS8xYwFUECyAHOhbwElshcBXBAsiBGQt4ia0QuIpgAeSQ/pdkd3e3+vv7A64GrqFjAVcRLIAcst/wk8lkgJXARcxYwFUECyCHeDyusrIySWyHoPDoWMBVBAvgGJizgBeMMcxYwFkEC+AYuDIEXkilUurr65NExwLuGVWwuP/++zV58mSVl5dr+vTpevXVVwtdF2AF7mUBL2QfT1VVVQFWAhRe3sHiySef1KJFi3TXXXdp27Ztuvjii9XS0qLdu3d7UR8QKDoW8EL6eKqqqlJxcXHA1QCFlXewuOeee/S9731P3//+93XmmWfq3nvvVWNjo1auXOlFfUCgmLGAF5ivgMtK8nlxb2+vtm7dqjvuuGPI87Nnz9bGjRuH/ZlUKqVUKpX5mjdohEm6Y/Hoo49q+/btAVcDV+zZs0cS8xVwU17B4sMPP1R/f78mTJgw5PkJEyZo3759w/5MW1ubli1bNvoKR6KyUjLG298Bd+RxvNTX10uS2tvb1d7e7mVViKD08TUs3teQL0uOmbyCRVosFhvytTHmqOfSWltbtXjx4szXiURCjY2No/m1gO9uv/12jRs3TgcPHgy6FDimuLhYN9xwQ9BlAAWXV7A4+eSTVVxcfFR3Yv/+/Ud1MdLi8bji8fjoKwQCVFdXp7vuuivoMgAgNPIa3iwrK9P06dO1bt26Ic+vW7dOM2fOLGhhAAAgfPLeClm8eLFuvvlmzZgxQ83NzXrwwQe1e/duzZs3z4v6AABAiOQdLK6//nodOHBAP//5z7V37141NTXp+eef16RJk7yoDwAAhEjMGH9HSBOJhGpra9XZ2cmlVgAAhMRIz998VggAACgYggUAACgYggUAACgYggUAACgYggUAACgYggUAACgYggUAACgYggUAACgYggUAACiYUX1s+olI3+gzkUj4/asBAMAopc/bx7tht+/BoqurS5LU2Njo968GAAAnqKurS7W1tTm/7/tnhQwMDGjPnj2qrq5WLBYr2N+bSCTU2Niojo4OPoPkOFirkWOt8sN6jRxrNXKs1ch5uVbGGHV1damhoUFFRbknKXzvWBQVFWnixIme/f01NTUceCPEWo0ca5Uf1mvkWKuRY61Gzqu1OlanIo3hTQAAUDAECwAAUDDOBIt4PK4lS5YoHo8HXYr1WKuRY63yw3qNHGs1cqzVyNmwVr4PbwIAAHc507EAAADBI1gAAICCIVgAAICCIVgAAICCcSZY3H///Zo8ebLKy8s1ffp0vfrqq0GXFKilS5cqFosNedTV1WW+b4zR0qVL1dDQoIqKCl166aV64403AqzYXxs2bNBVV12lhoYGxWIxPfPMM0O+P5L1SaVSWrhwoU4++WRVVlbqq1/9qt577z0f/yv8cby1+va3v33UsXbBBRcMeU0U1qqtrU3nnXeeqqurNX78eF1zzTV6++23h7yG4+pTI1kvjq3/WrlypaZOnZq56VVzc7NeeOGFzPdtO66cCBZPPvmkFi1apLvuukvbtm3TxRdfrJaWFu3evTvo0gJ19tlna+/evZnHjh07Mt/79a9/rXvuuUcrVqzQ5s2bVVdXpyuvvDLzWS6u6+7u1rRp07RixYphvz+S9Vm0aJHWrFmj1atX67XXXlMymdScOXPU39/v13+GL463VpL05S9/ecix9vzzzw/5fhTWqr29XfPnz9emTZu0bt06HT58WLNnz1Z3d3fmNRxXnxrJekkcW5I0ceJELV++XFu2bNGWLVs0a9YsXX311ZnwYN1xZRzwpS99ycybN2/Ic2eccYa54447AqooeEuWLDHTpk0b9nsDAwOmrq7OLF++PPNcT0+Pqa2tNb/73e98qtAeksyaNWsyX49kfT755BNTWlpqVq9enXnN+++/b4qKiszatWt9q91vR66VMcbMnTvXXH311Tl/JqprtX//fiPJtLe3G2M4ro7nyPUyhmPrWMaMGWN+//vfW3lchb5j0dvbq61bt2r27NlDnp89e7Y2btwYUFV22LlzpxoaGjR58mTdcMMNevfddyVJu3bt0r59+4asWTwe1//+7/9Gfs2kka3P1q1b1dfXN+Q1DQ0NampqiuQarl+/XuPHj9cXvvAF3XLLLdq/f3/me1Fdq87OTknS2LFjJXFcHc+R65XGsTVUf3+/Vq9ere7ubjU3N1t5XIU+WHz44Yfq7+/XhAkThjw/YcIE7du3L6Cqgnf++efr0Ucf1YsvvqhVq1Zp3759mjlzpg4cOJBZF9ZseCNZn3379qmsrExjxozJ+ZqoaGlp0WOPPaaXX35Zv/nNb7R582bNmjVLqVRKUjTXyhijxYsX66KLLlJTU5MkjqtjGW69JI6tbDt27FBVVZXi8bjmzZunNWvW6KyzzrLyuPL90029cuRHsBtjCvqx7GHT0tKS+fOUKVPU3Nysz3/+83rkkUcyw0+s2bGNZn2iuIbXX3995s9NTU2aMWOGJk2apOeee07XXnttzp9zea0WLFig7du367XXXjvqexxXR8u1Xhxbnzr99NP1+uuv65NPPtFTTz2luXPnqr29PfN9m46r0HcsTj75ZBUXFx+Vuvbv339UgouyyspKTZkyRTt37sxcHcKaDW8k61NXV6fe3l59/PHHOV8TVfX19Zo0aZJ27twpKXprtXDhQj377LN65ZVXNHHixMzzHFfDy7Vew4nysVVWVqbTTjtNM2bMUFtbm6ZNm6b77rvPyuMq9MGirKxM06dP17p164Y8v27dOs2cOTOgquyTSqX05ptvqr6+XpMnT1ZdXd2QNevt7VV7eztrJo1ofaZPn67S0tIhr9m7d6/+/ve/R34NDxw4oI6ODtXX10uKzloZY7RgwQI9/fTTevnllzV58uQh3+e4Gup46zWcqB5bwzHGKJVK2XlcFXwcNACrV682paWl5qGHHjL/+Mc/zKJFi0xlZaX597//HXRpgbntttvM+vXrzbvvvms2bdpk5syZY6qrqzNrsnz5clNbW2uefvpps2PHDnPjjTea+vp6k0gkAq7cH11dXWbbtm1m27ZtRpK55557zLZt28x//vMfY8zI1mfevHlm4sSJ5qWXXjJ/+9vfzKxZs8y0adPM4cOHg/rP8sSx1qqrq8vcdtttZuPGjWbXrl3mlVdeMc3Nzeazn/1s5Nbqhz/8oamtrTXr1683e/fuzTwOHjyYeQ3H1aeOt14cW59qbW01GzZsMLt27TLbt283d955pykqKjJ//vOfjTH2HVdOBAtjjPntb39rJk2aZMrKysy555475JKlKLr++utNfX29KS0tNQ0NDebaa681b7zxRub7AwMDZsmSJaaurs7E43FzySWXmB07dgRYsb9eeeUVI+mox9y5c40xI1ufQ4cOmQULFpixY8eaiooKM2fOHLN79+4A/mu8day1OnjwoJk9e7Y55ZRTTGlpqTn11FPN3Llzj1qHKKzVcGskyTz88MOZ13Bcfep468Wx9anvfve7mfPbKaecYi6//PJMqDDGvuOKj00HAAAFE/oZCwAAYA+CBQAAKBiCBQAAKBiCBQAAKBiCBQAAKBiCBQAAKBiCBQAAKBiCBQAAKBiCBQAAKBiCBQAAKBiCBQAAKBiCBQAAKJj/B9uWO6mThcBIAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "plt.plot(l,z, c='k')\n",
    "plt.vlines([0,93,206,300], ymin=0, ymax=6, color='r')\n",
    "plt.show()                # 描画"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1790d5af",
   "metadata": {},
   "source": [
    "左岸側の高水敷のsubsectionは次のとおりに定義できる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e798aa57",
   "metadata": {},
   "outputs": [],
   "source": [
    "sub = subsection( l[0:3], z[0:3], n[0:2] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "dc318ecb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "河積：134.25\n",
      "水面幅：91.0\n",
      "潤辺：91.35410196624969\n",
      "通水能：4232.411951420601\n",
      "合成粗度係数：0.04100000000000001\n"
     ]
    }
   ],
   "source": [
    "A, B, S, K, nd = sub.H2ABSKn(5.0)\n",
    "print('河積：{}\\n水面幅：{}\\n潤辺：{}\\n通水能：{}\\n合成粗度係数：{}'.format(A, B, S, K, nd))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4dce28d3",
   "metadata": {},
   "source": [
    "同様にすべてのsubsectionを含めた配列を作成する。\n",
    "\n",
    "その際に境界の座標点は重複させる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f863eb28",
   "metadata": {},
   "outputs": [],
   "source": [
    "section = [ subsection( l[0:3], z[0:3], n[0:2] )\n",
    "            , subsection( l[2:6], z[2:6], n[2:5] )\n",
    "            , subsection( l[5:], z[5:], n[5:] )\n",
    "          ]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45c2d615",
   "metadata": {},
   "source": [
    "次にsubsectionの配列となるクラスsectionを定義する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "d7eab3e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "class section(object):\n",
    "    def __init__(self, section):\n",
    "        self.subsections = section"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08e1f768",
   "metadata": {},
   "source": [
    "となるが、より抽象度を高めて以下のように定義する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "12d410e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "X = [l[0:3], l[2:6], l[5:]]\n",
    "Y = [z[0:3], z[2:6], z[5:]]\n",
    "N = [n[0:2], n[2:5], n[5:]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "41030594",
   "metadata": {},
   "outputs": [],
   "source": [
    "class section(object):\n",
    "    def __init__(self, L, Z, N):\n",
    "        self._subsections = [subsection(l,z,n) for l,z,n in zip(L,Z,N)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb2b322c",
   "metadata": {},
   "source": [
    "※ クラスのメンバ変数名の最初の_はクラスの外からアクセスして欲しくないことを示す暗黙のルール。pythonではprivate変数は設定できない。\n",
    "\n",
    "前回はLと定義した下流端からの追加距離をdistanceとして追加"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "d0530754",
   "metadata": {},
   "outputs": [],
   "source": [
    "class section(object):\n",
    "    def __init__(self, L, Z, N, distance=np.nan):\n",
    "        self._subsections = [subsection(l,z,n) for l,z,n in zip(L,Z,N)]\n",
    "        self.distance = distance"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "634466bc",
   "metadata": {},
   "source": [
    "次に断面の諸元を計算するメソッドを定義する。\n",
    "\n",
    "A,B,S,Kはsubsectionごと値の総和をとる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "2cc3e156",
   "metadata": {},
   "outputs": [],
   "source": [
    "class section(object):\n",
    "    def __init__(self, L, Z, N, distance=np.nan):\n",
    "        self._subsections = [subsection(l,z,n) for l,z,n in zip(L,Z,N)]\n",
    "        self.distance = distance\n",
    "        \n",
    "    def H2ABSKnSub(self, H):\n",
    "        num = len(self._subsections)\n",
    "        A = np.zeros(num)\n",
    "        B = np.zeros(num)\n",
    "        S = np.zeros(num)\n",
    "        K = np.zeros(num)\n",
    "        n = np.zeros(num)\n",
    "            \n",
    "        for nump in range(num) : \n",
    "            A[nump], B[nump], S[nump], K[nump], n[nump] = self._subsections[nump].H2ABSKn(H)\n",
    "            \n",
    "        return A, B, S, K, n\n",
    "    \n",
    "    def H2ABSK(self, H):\n",
    "        A, B, S, K, n = self.H2ABSKnSub(H)\n",
    "        return np.sum(A), np.sum(B), np.sum(S), np.sum(K)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a41f6fc",
   "metadata": {},
   "source": [
    "水位が5のときのA,B,S,Kを計算すると以下のとおりとなる。\n",
    "\n",
    "断面全体の粗度係数は総和では計算できないため、計算に含めていない。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "63643f8d",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "河積：858.0\n",
      "水面幅：296.0\n",
      "潤辺：298.3606797749979\n",
      "通水能：62445.88519359645\n"
     ]
    }
   ],
   "source": [
    "X = [l[0:3], l[2:6], l[5:]]\n",
    "Y = [z[0:3], z[2:6], z[5:]]\n",
    "N = [n[0:2], n[2:5], n[5:]]\n",
    "\n",
    "s = section(X,Y,N,0)\n",
    "A,B,S,K = s.H2ABSK(5.0)\n",
    "print('河積：{}\\n水面幅：{}\\n潤辺：{}\\n通水能：{}'.format(A, B, S, K))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b80ab81",
   "metadata": {},
   "source": [
    "課題と同様に、水位を0.01～6.0まで変化させて場合の以下の諸元を計算する。\n",
    " - 通水能（y軸：水位、x軸：通水能）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "587b45d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "Harr = np.arange(0.01,6.001,0.01)\n",
    "Karr = np.zeros_like(Harr)\n",
    "\n",
    "for i, H in enumerate(Harr) : \n",
    "    _,_,_,Karr[i] = s.H2ABSK(H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "1b0d39cf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxYAAAEiCAYAAABkykQ1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABGIklEQVR4nO3deVhUZf8/8PewDavjgqCIGu4aaiZJoOZCmiAuqYk7aj2PpeZWfUtzz8LHFtNMSzPNEjUXUFM0XHBJ0MRwq8zKNUFT2UQYGLh/f/ibk5OgDMzMmTPzfl0X1yXnHJjPOUc4vOdz7vuohBACREREREREleAgdwFERERERKR8DBZERERERFRpDBZERERERFRpDBZERERERFRpDBZERERERFRpDBZERERERFRpDBZERERERFRpDBZERERERFRpDBZERERERFRpDBZECvHee+8hPj7+geVJSUlQqVRISkqyeE1ERPZs9erVUKlUOH78uMHymzdvIigoCJ6enkhMTJSpOiLLY7AgUoiygsWTTz6J5ORkPPnkk5YvioiIDFy9ehUdO3bEn3/+iT179qBbt25yl0RkMU5yF0BUXvn5+XB1dYVKpZK7FKtSpUoVPP3003KXQURk986fP49nn30WRUVFOHDgAFq2bCl3SUQWxY4FyebXX3/F4MGD4evrC7VajXr16mHEiBHQarVSe/n777/H6NGjUbNmTbi7u0Or1aKkpAQLFixAs2bNoFar4ePjgxEjRuDq1asG3/+nn35CZGQkfHx8oFar4efnh549expst3HjRgQHB0Oj0cDd3R0NGjTA6NGjy1X/smXL0Lp1a3h6esLLywvNmjXDtGnTDLbJyMjAmDFj4O/vDxcXFwQEBGDOnDnQ6XQG22m1WsydOxfNmzeHq6sratSogS5duuDIkSMAAJVKhby8PHz11VdQqVRQqVTo3LkzgLJvhdq2bRtCQkLg7u4OLy8vdOvWDcnJyQbbzJ49GyqVCmfPnsXgwYOh0Wjg6+uL0aNHIzs7u1zHgYiIgLS0NHTo0AFOTk44fPgwQwXZJXYsSBYnT55Ehw4d4O3tjblz56Jx48ZIT0/Htm3bUFhYKG03evRo9OzZE19//TXy8vLg7OyMV155BcuXL8f48eMRGRmJixcvYsaMGUhKSsKJEyfg7e2NvLw8dOvWDQEBAfj000/h6+uLjIwM7N+/H7m5uQCA5ORkREVFISoqCrNnz4arqysuXbqEffv2PbL+9evXY+zYsXj11VfxwQcfwMHBAb///jt+/vlnaZuMjAy0a9cODg4OmDlzJho2bIjk5GTMmzcPFy9exKpVqwAAOp0O4eHhOHToECZNmoSuXbtCp9MhJSUFly9fRmhoKJKTk9G1a1d06dIFM2bMAHCvU1GW2NhYDB06FN27d8e6deug1WqxYMECdO7cGXv37kWHDh0Mtu/fvz+ioqLw4osv4vTp05g6dSoA4MsvvyznGSUisl+HDx/G7NmzUbduXXz//feoXbu23CURyUMQyaBr166iatWq4saNG6WuX7VqlQAgRowYYbD8l19+EQDE2LFjDZYfPXpUABDTpk0TQghx/PhxAUDEx8eXWcMHH3wgAIisrCyj6x8/fryoWrXqQ7cZM2aM8PT0FJcuXSr1dc+ePSuEEGLNmjUCgFixYsVDv5+Hh4eIjo5+YPn+/fsFALF//34hhBDFxcXCz89PtGzZUhQXF0vb5ebmCh8fHxEaGiotmzVrlgAgFixYYPA9x44dK1xdXUVJSclDayIismf6axUAodFoyrymEdkL3gpFFnf37l0cOHAAAwcORM2aNR+6bf/+/Q0+379/PwBg5MiRBsvbtWuH5s2bY+/evQCARo0aoVq1anjzzTfx2WefGXQS9J566ikAwMCBA/Htt9/ir7/+emCb4uJi6HQ66aOkpER6vaysLAwePBhbt27FzZs3H/ja7777Dl26dIGfn5/B9wgPDwcAHDhwAACQkJAAV1fXct+C9Sjnzp3DtWvXMHz4cDg4/PMj7unpif79+yMlJQV37941+JrevXsbfN6qVSsUFBTgxo0bJqmJiMiW9e7dG9nZ2Zg0aRKKi4vlLodINgwWZHGZmZkoLi6Gv7//I7f9dzv51q1bpS4HAD8/P2m9RqPBgQMH8MQTT2DatGl4/PHH4efnh1mzZqGoqAgA8MwzzyA+Ph46nQ4jRoyAv78/AgMDsW7dOul7NmzYEM7OztLH3LlzAQDDhw/Hl19+iUuXLqF///7w8fFBcHCwwbSC169fx/bt2w2+3tnZGY8//jgASGHk77//hp+fn0EIqIxHHaOSkhJkZmYaLK9Ro4bB52q1GsC9AfNERPRwM2bMwMyZMxEbG4thw4YxXJDd4hgLsrjq1avD0dHxgcHWpfn3DFD6P4DT09MfCCbXrl2Dt7e39HnLli2xfv16CCFw6tQprF69GnPnzoWbmxveeustAECfPn3Qp08faLVapKSkICYmBkOGDMFjjz2GkJAQbN++HVqtVvqefn5+0r9HjRqFUaNGIS8vDwcPHsSsWbMQGRmJ3377DfXr14e3tzdatWqFd999t9R903+vmjVr4vDhwygpKTFJuLj/GP3btWvX4ODggGrVqlX6dYiI6B9z5syBSqXCnDlzUFJSgrVr18LJiX9mkX1hx4Iszs3NDZ06dcLGjRtLvYXoYbp27QoA+OabbwyW//jjj/jll18QFhb2wNeoVCq0bt0aCxcuRNWqVXHixIkHtlGr1ejUqRP+97//Abg3oxRwL5wEBQVJH/cHCz0PDw+Eh4fj7bffRmFhIc6ePQsAiIyMxJkzZ9CwYUOD7/Hv7xUeHo6CggKsXr36ofuuVqvL1UFo2rQp6tSpg9jYWAghpOV5eXnYvHmzNFMUERGZ1uzZszFnzhx8++23GDJkyAMzABLZOkZpksVHH32EDh06IDg4GG+99RYaNWqE69evY9u2bfj888/L/LqmTZviv//9Lz755BM4ODggPDxcmhWqbt26mDx5MoB74xuWLl2Kvn37okGDBhBCYMuWLcjKypIeVjRz5kxcvXoVYWFh8Pf3R1ZWFhYtWgRnZ2d06tTpofX/5z//gZubG9q3b4/atWsjIyMDMTEx0Gg00tiNuXPnIjExEaGhoZgwYQKaNm2KgoICXLx4ETt37sRnn30Gf39/DB48GKtWrcLLL7+Mc+fOoUuXLigpKcHRo0fRvHlzDBo0CMC9kJOUlITt27ejdu3a8PLyQtOmTR+ozcHBAQsWLMDQoUMRGRmJMWPGQKvV4v3330dWVhbmz59foXNGRESPNnPmTDg4OGDGjBkQQmDdunXsXJDd4P90kkXr1q1x7NgxzJo1C1OnTkVubi5q1aqFrl27wsXF5aFfu2zZMjRs2BArV67Ep59+Co1Ggx49eiAmJka6Dahx48aoWrUqFixYgGvXrsHFxQVNmzbF6tWrER0dDQAIDg7G8ePH8eabb+Lvv/9G1apVERQUhH379knjIMrSsWNHrF69Gt9++y0yMzPh7e2NDh06YM2aNdKA9Nq1a+P48eN455138P777+Pq1avw8vJCQEAAevToId2O5OTkhJ07dyImJgbr1q3Dxx9/DC8vL7Ru3Ro9evSQXnPRokUYN24cBg0ahLt376JTp04PPLtCb8iQIfDw8EBMTAyioqLg6OiIp59+Gvv370doaGi5zhEREVXM9OnT4eDggLfffhslJSVYv349nJ2d5S6LyOxU4v57JYiIiIiIiCqAYyyIiIiIiKjSGCyIiIiIiKjSGCyIiIiIiKjSGCyIiIiIiKjSGCyIiIiIiKjSGCyIiIiIiKjSLP4ci5KSEly7dg1eXl5QqVSWfnkiIjKCEAK5ubnw8/ODg0Pl34viNYCISFmMuQ5YPFhcu3YNdevWtfTLEhFRJVy5cgX+/v6V/j68BhARKVN5rgMWDxZeXl4A7hVXpUoVS788EREZIScnB3Xr1pV+d1cWrwFERMpizHXA4sFC3/quUqUKLypERAphqtuWeA0gIlKm8lwHOHibiIiIiIgqjcGCiIiIiIgqjcGCiIiIiIgqzehg8ddff2HYsGGoUaMG3N3d8cQTTyA1NdUctRERERERkUIYNXg7MzMT7du3R5cuXZCQkAAfHx/88ccfqFq1qpnKIyIiIiIiJTAqWPzvf/9D3bp1sWrVKmnZY489ZuqaiIiIiIhIYYy6FWrbtm0ICgrCCy+8AB8fH7Rp0wYrVqwwV22lWrt2LRYsWIBr165Z9HWJiJTm4MGD+Oijj3D48GG5SyEiIgvT6XQ4fPgwpk+fjj/++MMir2lUx+LPP//EsmXLMGXKFEybNg3Hjh3DhAkToFarMWLEiFK/RqvVQqvVSp/n5ORUquB58+bh119/RXBwMPz8/Cr1vYiIbNmOHTuwYMECTJkyBR06dJClBlNfA4iIqGzXrl3D7t27kZCQgMTERGRlZQEAvL29MWnSJLO/vlHBoqSkBEFBQXjvvfcAAG3atMHZs2exbNmyMoNFTEwM5syZU/lK/z/9eI7MzEyTfU8iIluk/z0p5zg4U18DiIjoH0VFRThy5Ah27dqFhIQEnDx50mB99erV0b17dwQGBlqkHqOCRe3atdGiRQuDZc2bN8fmzZvL/JqpU6diypQp0uf6x4JXVLVq1QAwWBARPYr+96T+96YcTH0NICKyd1evXkVCQgJ27dqFPXv2GHSCVSoVgoKCEB4ejvDwcDz11FNwdHS0WG1GBYv27dvj3LlzBst+++031K9fv8yvUavVUKvVFauuFPp33vStHSIiKp3+96ScHQtTXwOIiOxNYWEhDh8+LHUlzpw5Y7De29sbzz33HMLDw9G9e3fUrFlTpkqNDBaTJ09GaGgo3nvvPQwcOBDHjh3D8uXLsXz5cnPV9wB2LIiIyscaOhZERGS8q1evYufOndi5cyf27t2LO3fuSOtUKhWCg4MRHh6OHj16oG3bthbtSjyMUcHiqaeeQlxcHKZOnYq5c+ciICAAH3/8MYYOHWqu+h7AjgURUflYQ8eCiIgerbi4GCkpKdixYwd27NiBU6dOGaz38fFBjx490KNHD3Tv3h01atSQqdKHMypYAEBkZCQiIyPNUUu5sGNBRFQ+7FgQEVmvW7duYffu3dixYwd27dqF27dvS+v0XYmePXsiPDwcbdq0gYODUU+JkIXRwUJu+gskOxZERGUTQki/JxksiIjkJ4TAqVOnpK5ESkoKSkpKpPVVq1ZFjx49EBERgR49esg6VqKiFBcsON0sEdGj5ebmShcs3gpFRCSPvLw87N27Fzt27MDOnTtx9epVg/WBgYHo2bMnevbsiZCQEDg5Ke5PcwOKq563QhERPZr+d6RarYabm5vM1RAR2Y8//vhD6kokJSWhsLBQWufm5oawsDD07NkTERERqFevnoyVmp7iggUHbxMRPRoHbhMRWYZ+Olh9mPj3oxkCAgKkING5c2ebfrNHccGCHQsiokfjwG0iIvO5ffs2EhISsG3bNuzatcvgIXVOTk7o0KGDdItTs2bNoFKpZKzWchQXLPTvvuXn50Or1fLBS0REpWDHgojItM6fP4/t27dj27ZtOHz4MIqLi6V1Pj4+iIiIQEREBLp37w6NRiNjpfJRXLDQaDRQqVTSjCe+vr5yl0REZHXYsSAiqpzi4mIkJydj27Zt2L59O3799VeD9YGBgejduzd69eqFdu3aKWI6WHNTXLBwcHCARqNBVlYWgwURURk41SwRkfFyc3Oxe/dubN++HTt27MCtW7ekdU5OTujUqZMUJgICAmSs1DopLlgA91r7WVlZHGdBRFQG/e9H3gpFRPRwly9flm5x+vcsTtWqVUNERAR69eqFHj162O0tTuWlyGBRrVo1XLx4kTNDERGVgbdCERGVrqSkBKmpqVKYOHnypMH6Ro0aoXfv3ujduzfat2+v+GdLWJIijxQfkkdE9HAcvE1E9A+tVot9+/YhPj4e27dvR3p6urTOwcEBoaGh6NWrF3r37o2mTZvazSxOpqbIYMEpZ4mIHo4dCyKyd9nZ2UhISEB8fDx27tyJ3NxcaZ2npyeee+459O7dGxEREfD29paxUtuhyGDBh+QRET0cOxZEZI8yMjKwdetWxMfHY+/evSgqKpLW1a5dG3379kXv3r3RpUsXPrLADBQZLNixICJ6OHYsiMhenD9/HnFxcYiPj0dKSgqEENK6pk2b4vnnn0ffvn3x1FNPcUpYM1NksGDHgojo4TjdLBHZKiEEUlNTER8fj7i4OPz8888G69u1ayeFiWbNmslUpX1SZLBgx4KI6OE43SwR2ZKioiIcPHgQ8fHxiI+Px9WrV6V1Tk5O6NKlC/r27Ys+ffqgTp06MlZq3xQdLNixICJ6kFarRX5+PgB2LIhIufLz87Fr1y7ExcXhu+++M3hD2cPDA+Hh4ejbty969uzJN1GshCKDBaebJSIqm/5NF5VKhSpVqshbDBGREe7cuYOdO3di06ZN2LlzJ/Ly8qR1NWvWRO/evdG3b188++yzcHV1lbFSKo0igwVvhSIiKpv+d6NGo+FARSKyetnZ2fjuu++wadMm7Nq1CwUFBdK6evXqoX///nj++ecRGhoKR0dHGSulR1FksODgbSKisnGqWSKydrdv38bWrVuxefNmJCYmorCwUFrXqFEj9O/fHwMGDEDbtm35sDoFUWSwuH+MRUlJCd+RIyK6D6eaJSJrdOPGDcTFxWHz5s3Yv38/dDqdtK558+YYMGAA+vfvj1atWjFMKJQig4X+XTghBHJzc6HRaOQtiIjIinCqWSKyFteuXcOWLVuwadMmHDp0CCUlJdK61q1bo3///ujfvz9atGghY5VkKooMFq6urnB1dUVBQQEyMzMZLIiI7sOpZolITpcvX8bmzZuxadMmHDlyxGBdUFCQ1Jlo1KiRTBWSuRgVLGbPno05c+YYLPP19UVGRoZJiyqPatWqIT09neMsiIj+hbdCEZGlXbt2DRs3bsSGDRuQnJxssC4kJAQDBgxAv3798Nhjj8lTIFmE0R2Lxx9/HHv27JE+l2t0ftWqVZGens6ZoYiI/oWDt4nIEq5fv45Nmzbh22+/xaFDhyCEAHBvquuOHTtKYYIPrLMfRgcLJycn1KpVyxy1GIVTzhIRlY4dCyIyl5s3b2LLli3YsGEDkpKSDMZMhIaGIioqCgMGDICfn5+MVZJcjA4W58+fh5+fH9RqNYKDg/Hee++hQYMG5qjtoTjlLBFR6dixICJTyszMRFxcHDZs2IC9e/eiuLhYWteuXTsMHDgQL7zwAurVqydjlWQNjAoWwcHBWLNmDZo0aYLr169j3rx5CA0NxdmzZ1GjRo1Sv0ar1UKr1Uqf5+TkVK7i/48dCyKi0llTx8Jc1wAiMq/s7Gxs3boVGzZsQGJiIoqKiqR1bdq0QVRUFAYOHIiAgAAZqyRrY1SwCA8Pl/7dsmVLhISEoGHDhvjqq68wZcqUUr8mJibmgQHfpsCOBRFR6axpullzXQOIyPTy8vKwfft2rF+/HgkJCQYPrWvZsqUUJho3bixjlWTNKjXdrIeHB1q2bInz58+Xuc3UqVMNQkdOTg7q1q1bmZcFwI4FEVFZrGm6WXNdA4jINIqKirBnzx7ExsYiLi4OeXl50rpmzZohKioKUVFRaN68uYxVklJUKlhotVr88ssv6NixY5nbqNVqqNXqyrxMqe5/+jYREf3Dmm6FMtc1gIgqTgiB5ORkxMbG4ttvv8Xff/8trQsICMDgwYMxaNAgBAYG8gnYZBSjgsXrr7+OXr16oV69erhx4wbmzZuHnJwcREdHm6u+MunfiWPHgojoHyUlJcjOzgZgHR0LIrIeZ8+eRWxsLGJjY3Hx4kVpuY+PD6KiojBkyBAEBwczTFCFGRUsrl69isGDB+PmzZuoWbMmnn76aaSkpKB+/frmqq9MvBWKiOhBOTk50lzyDBZEdPnyZaxfvx5r167FqVOnpOWenp7o168fhgwZgrCwMDg5VeomFiIARgaL9evXm6sOo3HwNhHRg/S/E11dXeHq6ipvMUQki1u3bmHTpk1Yu3YtDh06JC13dnZGREQEhgwZgsjISLi7u8tYJdkixcZTdiyIiB5kTeMriMhyCgoKsH37dqxZswa7du2CTqcDcO8p2J06dcKQIUPQv39/VK9eXeZKyZYpNliwY0FE9CBrmmqWiMxLCIEjR45gzZo12LBhgzS+Crj3rIkhQ4Zg0KBB8Pf3l7FKsieKDRb6i2ZBQQEKCgrY8icignVNNUtE5vHnn3/i66+/xtdff40//vhDWl6vXj0MGzYMw4YN4/SwJAvFBgsvLy84ODigpKQEWVlZqFWrltwlERHJjrdCEdmm7OxsbNy4EWvWrDEYN+Hp6YkXXngBI0aMwDPPPAMHBwcZqyR7p9hg4eDgAI1Gg8zMTGRmZjJYEBHhn1uh2LEgUj6dTofvv/8ea9asQXx8PLRaLYB7fwM9++yzGDFiBPr27QsPDw+ZKyW6R7HBArj3jlxmZibHWRAR/X/sWBApX1paGtasWYPY2Fhcv35dWv74448jOjoaQ4YMQZ06dWSskKh0ig4WfEgeEZEhdiyIlOn27duIjY3Fl19+iZ9++klaXrNmTQwZMgQjRoxAmzZt+PA6smqKDhaccpaIyBA7FkTKUVJSgr179+LLL79EXFycdKuTi4sLevfujejoaDz33HNwdnaWuVKi8lF0sOCUs0REhjjdLJH1u3jxIlavXo1Vq1bh8uXL0vLWrVvjxRdfxJAhQ1CjRg0ZKySqGEUHC3YsiIgMcbpZIutUUFCAuLg4fPnll9i7dy+EEADu/awOHToUo0eP5q1OpHiKDhbsWBARGeKtUETW5cSJE/jyyy+xdu1ag79XwsLC8OKLL6Jv375wc3OTr0AiE1J0sGDHgojIEAdvE8kvOzsba9euxYoVK5CWliYtr1evHkaNGoWRI0fisccek60+InOxiWDBjgUR0T3sWBDJQwiBH3/8EZ9//jnWr1+Pu3fvArg3ELtfv34YPXo0wsLC+AA7smmKDhacbpaI6B8FBQXSrDLsWBBZRk5ODmJjY/H5558bdCdatGiBMWPGYNiwYahevbp8BRJZkKKDBW+FIiL6h/53oYODA7y8vGSuhsi2HT9+HMuXL0dsbCzy8vIAAGq1GgMHDsSYMWMQGhrKgdhkdxQdLDh4m4joH/ePr+DtFkSml5ubi3Xr1uHzzz/HiRMnpOXNmjXDmDFjMGLECHYnyK4pOliwY0FE9A9ONUtkHqdOncLSpUuxdu1a3LlzB8C9sRMDBgzAmDFj0LFjR3YniKDwYKG/eObk5KCkpITv0BGRXePAbSLTKSoqQlxcHJYsWYJDhw5Jy5s0aSJ1J7y9vWWskMj6KDpY6C+eQghkZ2fzYkpEdo1TzRJVXnp6OlasWIHPPvsM6enpAABHR0f069cPr7zyCjp37szuBFEZFB0sXFxc4O7ujrt37yIrK4vBgojsGjsWRBUjhMCRI0ewZMkSbNq0CTqdDgDg6+uLMWPG4L///S/q1Kkjc5VE1k/RwQK4987c3bt3kZmZiYCAALnLISKSDTsWRMa5e/cuYmNjsWTJEpw8eVJa3r59e4wfPx79+vWDi4uLjBUSKYvig0W1atVw7do1DuAmIrvHjgVR+Vy+fBlLlizBihUrpEDu5uaGoUOHYty4cXjiiSdkrY9IqRQfLDjlLBHRPfrfgwwWRKVLSUnBwoULsXnzZhQXFwMAGjRogLFjx2LUqFGcKpaokhQfLDjlLBHRPZxuluhBOp0OW7ZswcKFC5GSkiItDwsLw6RJkxAREcFZJYlMpFI/STExMVCpVJg0aZKJyjEeOxZERPfwViiif2RlZeGDDz5Aw4YNERUVhZSUFLi4uGDUqFE4efIk9uzZg8jISIYKIhOqcMfixx9/xPLly9GqVStT1mM0diyIiO7h4G0i4Pfff8eiRYuwatUq5OXlAQBq1qyJsWPH4pVXXoGvr6/MFRLZrgoFizt37mDo0KFYsWIF5s2bZ+qajKIPFuxYEJG9Y8eC7NnRo0exYMECxMXFQQgBAAgMDMTkyZMxZMgQuLq6ylwhke2rULAYN24cevbsiWefffaRwUKr1UKr1Uqf5+TkVOQly6R/Z44dCyKyd9bYsTD3NYDsmxACu3btwv/+9z8cOHBAWh4eHo4pU6YgLCyMD7MjsiCjg8X69euRmpqK48ePl2v7mJgYzJkzx+jCyou3QhERAcXFxcjOzgZgXR0Lc18DyD7pdDps2LABCxYswKlTpwAATk5OGDZsGN544w20aNFC5gqJ7JNRI5auXLmCiRMnYu3ateVuKU6dOhXZ2dnSx5UrVypUaFk4eJuIyLATYE0dC3NfA8i+5OXl4ZNPPkGjRo0wbNgwnDp1Cp6enpgyZQouXLiAVatWMVQQyciojkVqaipu3LiBtm3bSsuKi4tx8OBBLFmyBFqtFo6OjgZfo1aroVarTVNtKdixICL653egu7u7VT0p2NzXALIPmZmZ+OSTT7B48WLcunULwL0B2RMnTsTYsWOtqktHZM+MChZhYWE4ffq0wbJRo0ahWbNmePPNNx8IFZbAjgUREQduk226efMmFi5ciE8++QS5ubkA7j3Q7vXXX8fIkSPh5uYmc4VEdD+jgoWXlxcCAwMNlnl4eKBGjRoPLLeU+zsWQggO0iIiu2SNA7eJKur69ev48MMPsXTpUmnK2FatWmHatGkYMGCALG9kEtGjKf7J2/qLaGFhIQoKCvjuBRHZJXYsyBakp6fj/fffx2effYb8/HwAwJNPPokZM2agd+/efJgdkZWrdLBISkoyQRkV5+XlBUdHRxQXFyMzM5PBgojsEjsWpGRXrlzBggULsGLFCml64nbt2mHmzJmIiIjg3QhECqH46K9SqfgsCyKye+xYkBJlZGRgwoQJaNSokTQJTPv27bF7926kpKSgZ8+eDBVECqL4W6GAe+/Q3bp1iwO4ichu6X//MViQEty+fRsLFizA4sWLpVueOnXqhFmzZqFz584ME0QKZRPBglPOEpG90//+461QZM1yc3Px8ccf44MPPpCevRIcHIx3330XYWFhMldHRJVlE8GCU84Skb3jrVBkzfLz87F06VLMnz8fN2/eBHBvlqd58+YhMjKSHQoiG2ETwYIdCyKydxy8TdaouLgYX3/9NaZPn46//voLANC4cWPMnTsXAwcO5CxPRDbGJoIFOxZEZO/YsSBrs2fPHrz++us4efIkAKBu3bqYPXs2RowYAScnm/jzg4j+xSZ+stmxICJ7x44FWYuzZ8/ijTfeQEJCAgBAo9Fg+vTpGD9+PFxdXWWujojMyaaCBTsWRGSv2LEguWVkZGDWrFn44osvUFJSAicnJ4wbNw4zZsxAjRo15C6PiCzAJoIFn2NBRPZMCMHpZkk2Op0OS5YswcyZM5GbmwsA6N+/P2JiYtC4cWOZqyMiS7KJYKG/kO7evRtNmjSRuRr7Vr9+fcTFxcHT01PuUojsRn5+PgoLCwHwViiyrB9++AFjx47FqVOnAABPPfUUPvroI3To0EHmyohIDjYRLFq2bAng3sX1/PnzMldj386fP49Dhw4hPDxc7lKI7Ia+W+vo6MhQTxbx999/480338SqVasAANWrV0dMTAxeeuklzvREZMdsIli0aNECf/75J65duyZ3KXZt0qRJOH78OG9JI7Kw+wdu83kAZE5CCKxcuRL/93//J/2uf/HFFzF//nx4e3vLXB0Ryc0mggUABAQEICAgQO4y7Fq9evVw/PhxDqInsjAO3CZL+P333/Hf//4X+/fvBwA88cQTWLp0KUJCQmSujIisBfuVZDKc9pdIHpxqlsxJp9Ph/fffR8uWLbF//364ubnhww8/xI8//shQQUQGbKZjQfLjgwqJ5MGOBZlLWloaXnrpJaSmpgIAwsLCsHz5cjRo0EDmyojIGrFjQSbDjgWRPDjVLJlaQUEBpk2bhqCgIKSmpqJq1ar48ssvkZiYyFBBRGVix4JMhg8qJJKHPszzVigyhaNHjyI6Ohrnzp0DAAwYMACffPIJatWqJXNlRGTt2LEgk+GDConkwVuhyBQKCwsxffp0hIaG4ty5c6hduza2bNmCjRs3MlQQUbmwY0Emw1uhiOTBwdtUWWfOnMHw4cORlpYGABg8eDCWLFmC6tWry1sYESkKOxZkMhy8TSQPdiyoooqLi/H++++jbdu2SEtLQ40aNfDtt98iNjaWoYKIjMaOBZkMOxZE8mDHgiriwoULGDFiBA4fPgwAiIyMxIoVK3jbExFVGDsWZDL6P2qys7NRUlIibzFEdoQdCzLW+vXr8cQTT+Dw4cPw8vLCypUrsW3bNoYKIqoUBgsyGX2wEEIgJydH3mKI7Ag7FlReeXl5GD16NAYPHoycnBx06NABp06dwujRo6FSqeQuj4gUzqhgsWzZMrRq1QpVqlRBlSpVEBISgoSEBHPVRgrj6uoKV1dXABxnQWRJ7FhQeaSlpaFt27ZYtWoVHBwcMHPmTOzfvx+PPfaY3KURkY0wKlj4+/tj/vz5OH78OI4fP46uXbuiT58+OHv2rLnqI4XhOAsiy9LpdMjNzQXAYEGlE0Lg008/RXBwMM6dO4c6depg3759mDNnDpycONSSiEzHqGDRq1cvREREoEmTJmjSpAneffddeHp6IiUlxVz1kcIwWBBZVnZ2tvRvjUYjYyVkjfLy8jBs2DCMHz8ehYWF6NOnD06ePIlOnTrJXRoR2aAKv1VRXFyMjRs3Ii8vDyEhIWVup9VqodVqpc95771t45Sz9iU5ORnz5s0z+BknyyooKAAAeHp6wtnZWeZqHsRrgHzOnz+Pfv364cyZM3B0dMSCBQswefJkjqUgIrMxOlicPn0aISEhKCgogKenJ+Li4tCiRYsyt4+JicGcOXMqVSQpBzsW9mXhwoXYuXOn3GUQgEaNGsldQql4DZDH3r17MWDAAGRlZcHX1xfffvstnnnmGbnLIiIbZ3SwaNq0KdLS0pCVlYXNmzcjOjoaBw4cKDNcTJ06FVOmTJE+z8nJQd26dSteMVk1dizsy82bNwEAr7766kM7l2ReKpUKHTt2lLuMUvEaYHnLly/HuHHjoNPp8PTTT2Pz5s3w8/OTuywisgNGBwsXFxfpnbGgoCD8+OOPWLRoET7//PNSt1er1VCr1ZWrkhSDHQv7og+Q4eHhCA8Pl7cYskq8BlhOcXExXnvtNSxatAgAMGTIEKxcuVKarY+IyNwqPR2EEIL3V5OEHQv7wmlOiaxDTk4OBg0aJE0BP2/ePEybNo3jKYjIoowKFtOmTUN4eDjq1q2L3NxcrF+/HklJSdi1a5e56iOFYcfCvvDBbETyu3DhAnr16oWzZ8/Czc0Na9aswYABA+Qui4jskFHB4vr16xg+fDjS09Oh0WjQqlUr7Nq1C926dTNXfaQwDBb2o6SkRJrqlB0LInn88MMP6Nu3L27evAk/Pz9s3boVQUFBcpdFRHbKqGCxcuVKc9VBNoK3QtmP7OxsCCEAsGNBJIfY2FiMGjUKhYWFePLJJ7Ft2zbUqVNH7rKIyI4Z9YA8okdhx8J+6MOjm5sbB+cSWZAQAjExMRg6dCgKCwvRr18/HDx4kKGCiGTHYEEmxY6F/eDAbSLL0+l0eOWVVzBt2jQAwOuvv46NGzfCw8ND5sqIiEwwKxTR/dixsB8cuE1kWXfu3MGgQYOwY8cOqFQqLF68GOPHj5e7LCIiCYMFmZT+j0ytVouCggLOn27D2LEgspyMjAxERkYiNTUVrq6uiI2NxfPPPy93WUREBngrFJmUl5cXHBzu/bdi18K2sWNBZBm//vorQkJCkJqaCm9vb+zfv5+hgoisEoMFmZSDg4P0hyaDhW1jx4LI/H744QeEhobi4sWLaNiwIZKTk/H000/LXRYRUakYLMjkOIDbPujPL4MFkXns2LED3bp1Q2ZmJoKDg5GcnIxGjRrJXRYRUZkYLMjkOIDbPujPL2+FIjK9b775Bn369EF+fj4iIiKwb98+1KxZU+6yiIgeisGCTI4dC/vAW6GIzGPx4sUYPnw4iouLMWzYMMTHx8Pd3V3usoiIHonBgkyOHQv7wMHbRKYlhMDMmTMxceJEAMCECRPw1VdfwdnZWebKiIjKh9PNksmxY2Ef2LEgMp3i4mKMHz8en332GQDgnXfewdtvvw2VSiVzZURE5cdgQSbHjoV9YMeCyDQKCwsxYsQIbNiwASqVCp9++ileeeUVucsiIjIagwWZHIOFfWDHgqjyCgoK8MILL+C7776Ds7MzvvnmGwwcOFDusoiIKoTBgkyOt0LZB3YsiCrn7t27eP755/H999/D1dUVcXFx6NGjh9xlERFVGIMFmRw7FrYvPz8fWq0WADsWRBVx584d9OrVC0lJSfDw8MD27dvRpUsXucsiIqoUBgsyOXYsbJ8+NDo4OMDLy0vmaoiUJTs7GxEREThy5Ai8vLyQkJCA9u3by10WEVGlMViQybFjYfvuvw2Ks9YQld/t27fx3HPP4fjx46hatSp2796Ndu3ayV0WEZFJMFiQybFjYfs4cJvIeDdv3sSzzz6LkydPwtvbG4mJiXjiiSfkLouIyGT4gDwyOf0fm9nZ2SguLpa5GjIHDtwmMs7t27fRrVs3nDx5Er6+vkhKSmKoICKbw2BBJnf/H5s5OTnyFUJmw44FUfllZWWhe/fuSEtLk0LF448/LndZREQmx2BBJufi4gJ3d3cAHGdhq9ixICqf7OxsPPfcc0hNTUXNmjWxd+9eNGvWTO6yiIjMgsGCzIIDuG0bOxZEj5abm4vw8HAcO3YMNWrUwJ49e9ipICKbxmBBZsEB3LZNf14ZLIhKd+fOHURERCA5ORnVqlXDnj170KpVK7nLIiIyK6OCRUxMDJ566il4eXnBx8cHffv2xblz58xVGykYOxa2TX9eeSsU0YPy8/PRq1cvHD58GBqNhrM/EZHdMCpYHDhwAOPGjUNKSgoSExOh0+nQvXt35OXlmas+Uih2LGwbb4UiKl1RURGioqKQlJQELy8v7N69G23btpW7LCIiizDqORa7du0y+HzVqlXw8fFBamoqnnnmGZMWRsrGjoVt4+BtogeVlJTgxRdfxPbt2+Hq6orvvvsOwcHBcpdFRGQxlRpjkZ2dDQCoXr26SYoh26EPFuxY2CZ2LIgMCSEwadIkfP3113B0dMTGjRv5hhsR2Z0KP3lbCIEpU6agQ4cOCAwMLHM7rVYLrVYrfc7nGtgH/TvZ7FjYJnYsqLzs5Rowd+5cfPLJJwCAr776CpGRkTJXRERkeRXuWIwfPx6nTp3CunXrHrpdTEwMNBqN9FG3bt2KviQpCG+Fsm3sWFB52cM1YPHixZg9ezYA4JNPPsHQoUPlLYiISCYVChavvvoqtm3bhv3798Pf3/+h206dOhXZ2dnSx5UrVypUKCkLB2/bruLiYuldZ3Ys6FFs/Rqwbt06TJw4EcC9rsX48eNlroiISD5G3QolhMCrr76KuLg4JCUlISAg4JFfo1aroVarK1wgKRM7FrZLP7YKYMeCHs2WrwEHDhzAyJEjAQATJkzA9OnT5S2IiEhmRgWLcePGITY2Flu3boWXlxcyMjIAABqNBm5ubmYpkJSJHQvbpQ+LHh4ecHZ2lrkaInn8/PPP6Nu3LwoLC9G/f38sXLgQKpVK7rKIiGRl1K1Qy5YtQ3Z2Njp37ozatWtLHxs2bDBXfaRQ7FjYLg7cJnuXnp6OiIgIZGVlITQ0FF9//TUcHCo1ySIRkU0w+lYoovK4v2MhhOA7eTaEA7fJnt25cweRkZG4dOkSGjdujK1bt7JjT0T0//EtFjIL/R+dhYWFyM/Pl7kaMiV2LMhe6XQ6DBw4ECdOnEDNmjWRkJAAb29vucsiIrIaDBZkFp6ennB0dATA26FsDTsWZK9ef/11JCQkwM3NDdu3b0fDhg3lLomIyKowWJBZqFQqDuC2UexYkD1auXIlFi1aBAD4+uuvERwcLHNFRETWh8GCzIYDuG0TOxZkbw4fPoxXXnkFADBnzhz0799f5oqIiKwTgwWZDTsWtkl/PhksyB5cunQJ/fr1Q1FREV544QXMmDFD7pKIiKwWgwWZDTsWtkl/PnkrFNm6O3fuoE+fPvj777/Rpk0brF69mjPcERE9BIMFmQ07FraJt0KRPRBCYPTo0Th58iR8fX2xdetWuLu7y10WEZFVY7Ags2HHwjZx8DbZg0WLFmHjxo1wdnbGli1bULduXblLIiKyegwWZDYMFraJHQuydYcPH8Ybb7wBAPjoo48QGhoqc0VERMrAYEFmw1uhbBM7FmTLrl+/joEDB0Kn02Hw4MEYN26c3CURESkGgwWZDTsWtkcIwY4F2SydTodBgwYhPT0dLVq0wPLlyzlYm4jICAwWZDbsWNie/Px8FBUVAWCwINszc+ZMJCUlwdPTE5s3b4anp6fcJRERKQqDBZkNOxa2R38uHR0d4eHhIXM1RKazb98+zJ8/HwDwxRdfoFmzZjJXRESkPAwWZDbsWNie+2+D4i0iZCtu3bqF4cOHQwiBl156CVFRUXKXRESkSAwWZDbsWNgeDtwmW6MPE9euXUPTpk3x8ccfy10SEZFiMViQ2eiDRW5uLnQ6nczVkClw4DbZmuXLlyM+Ph7Ozs5Yt24db/EjIqoEBgsyG41GI/07OztbxkrIVNixIFvy66+/YvLkyQCA+fPno02bNjJXRESkbAwWZDbOzs7SrCq8Hco2sGNBtqK4uBgjR45Efn4+unfvjkmTJsldEhGR4jFYkFlxALdtYceCbMWHH36Io0ePQqPRYOXKlXBw4OWQiKiy+JuUzIoDuG0LOxZkC3755RfMnDkTALBw4UL4+/vLXBERkW1gsCCzYsfCtjBYkNIVFxdj1KhR0Gq1CA8Px8iRI+UuiYjIZjBYkFmxY2FbeCsUKd1HH30k3QK1fPlyPo+FiMiEGCzIrPTBgh0L28COBSnZhQsXMGvWLAD3AgZvgSIiMi2jg8XBgwfRq1cv+Pn5QaVSIT4+3gxlka3Qv7PNjoVtYMeClEoIgQkTJiA/Px+dO3fGqFGj5C6JiMjmGB0s8vLy0Lp1ayxZssQc9ZCN4a1QtoUdC1KqrVu34rvvvoOzszOWLl3KW6CIiMzAydgvCA8PR3h4uDlqIRvEwdu2hR0LUqK8vDxMmDABAPD666+jefPmMldERGSbOMaCzIodC9uh0+mQm5sLgB0LUpZ33nkHV65cQf369TF9+nS5yyEisllGdyyMpdVqodVqpc9zcnLM/ZJkRdixsB3Z2dnSv9mxoPKS+xpw4cIFLFy4EACwePFiuLu7W/T1iYjsidk7FjExMdBoNNJH3bp1zf2SZEXYsbAd+nPo6ekJJyezvydBNkLua8DUqVNRWFiIsLAw9OrVy6KvTURkb8weLKZOnYrs7Gzp48qVK+Z+SbIinG7WdnDgNlWEnNeAlJQUbNiwASqVCh9++CEHbBMRmZnZ33ZUq9VQq9XmfhmyUvdPNyuE4IVdwThwmypCrmuAEAKvvfYaAGDkyJFo3bq1xWsgIrI3RgeLO3fu4Pfff5c+v3DhAtLS0lC9enXUq1fPpMWR8unf3dbpdMjLy4Onp6fMFVFFsWNBShIfH48jR47A3d0d77zzjtzlEBHZBaODxfHjx9GlSxfp8ylTpgAAoqOjsXr1apMVRrbB3d0dTk5O0Ol0yMrKYrBQMHYsSClKSkowc+ZMAMDkyZNRp04dmSsiIrIPRgeLzp07QwhhjlrIBqlUKlSrVg1///03MjMz4e/vL3dJVEHsWJBSbN68GWfOnEGVKlWk26GIiMj8+BwLMjtOOWsb2LEgJSgpKcGcOXMAAJMmTWIQJiKyIAYLMjtOOWsb2LEgJdi0aRPOnj0LjUaDyZMny10OEZFdYbAgs+OUs7aBwYKsnRACMTExAO51K9hdIyKyLAYLMrv7p5wl5eKtUGTtkpKSkJaWBjc3N7z66qtyl0NEZHcYLMjseCuUbWDHgqzdhx9+CAAYNWoUatSoIXM1RET2h8GCzI6Dt20DOxZkzX755Rfs2LEDKpUKkyZNkrscIiK7xGBBZseOhW1gx4Ks2eLFiwEAvXv3RuPGjWWuhojIPjFYkNmxY6F8Qgh2LMhq5eXlYe3atQCACRMmyFwNEZH9YrAgs2PHQvny8vKg0+kAsGNB1mfjxo3Izc1Fw4YN0blzZ7nLISKyWwwWZHbsWCif/tw5OzvD3d1d3mKI/uWLL74AALz44otwcOBljYhILvwNTGbHjoXy6c9d1apVoVKpZK6G6B/nzp3DDz/8AEdHR0RHR8tdDhGRXWOwILNjsFA+Dtwma7VhwwYAQPfu3eHn5ydzNURE9o3BgsxOfytUXl4eioqK5C2GKoQDt8labdy4EQAwcOBAmSshIiIGCzI7jUYj/ZvjLJSJHQuyRr/88gvOnDkDZ2dn9OnTR+5yiIjsHoMFmZ2TkxO8vLwAMFgoFTsWZI02bdoEAOjWrRtDLxGRFWCwIIvgOAtlY8eCrFFCQgIA4Pnnn5e5EiIiAhgsyEI45ayy6c8bgwVZi6ysLBw9ehTAvYHbREQkPwYLsgh2LJTt/ulmiazBvn37UFJSgqZNm6JevXpyl0NERGCwIAthsFA23gpF1mbv3r0A2K0gIrImDBZkEbwVStk4eJusTUpKCgCgY8eOMldCRER6DBZkEexYKBs7FmRNCgoKcOrUKQBAu3btZK6GiIj0GCzIItixUDZ2LMiapKWlQafTwcfHh+MriIisCIMFWQQ7FsrGjgVZk9OnTwMAnnzySahUKpmrISIiPQYLsgh2LJSrqKgIeXl5ANixIOtw/vx5AECTJk1kroSIiO5XoWCxdOlSBAQEwNXVFW3btsWhQ4dMXRfZGHYslOv+MMhgQdZAHywaNWokcyVERHQ/o4PFhg0bMGnSJLz99tv46aef0LFjR4SHh+Py5cvmqI9sBIOFcunPWZUqVeDo6ChzNUSQrjcBAQEyV0JERPdzMvYLPvroI7z44ot46aWXAAAff/wxdu/ejWXLliEmJsbkBZJt0L/TnZmZifz8fHmLIaNcv34dALsVZD1u374NAPD29pa5EiIiup9RwaKwsBCpqal46623DJZ3794dR44cKfVrtFottFqt9HlOTk4FyiSl03csbt++DXd3d5mroYrgwG2qCHNcA/TBonr16pX+XkREZDpG3Qp18+ZNFBcXw9fX12C5r68vMjIySv2amJgYaDQa6aNu3boVr5YUy9fXF08++aTcZVAlREREyF0CKZCprwFCCPTo0QNhYWHsWBARWRmjb4UC8MD0fkKIMqf8mzp1KqZMmSJ9npOTw3BhhxwdHfHjjz9KswuRsjg4OMDDw0PuMkiBTH0NUKlU2LBhgylKIyIiEzMqWHh7e8PR0fGB7sSNGzce6GLoqdVqqNXqildINsPBwQFeXl5yl0FEFsRrABGR/TDqVigXFxe0bdsWiYmJBssTExMRGhpq0sKIiIiIiEg5jL4VasqUKRg+fDiCgoIQEhKC5cuX4/Lly3j55ZfNUR8RERERESmA0cEiKioKt27dwty5c5Geno7AwEDs3LkT9evXN0d9RERERESkABUavD127FiMHTvW1LUQEREREZFCGf3kbSIiIiIion9jsCAiIiIiokqr0K1QlSGEAMAncBMRKYH+d7X+d3dl8RpARKQsxlwHLB4scnNzAYAPySMiUpDc3FxoNBqTfB+A1wAiIqUpz3VAJUz1NlQ5lZSU4Nq1a/Dy8irzad0Po39q65UrV1ClShUzVCgvW94/7pty2fL+2fK+AZXfPyEEcnNz4efnBweHyt89W9lrAGD756wyeGzKxmNTOh6XsvHY3GPMdcDiHQsHBwf4+/tX+vtUqVLFpk+yLe8f9025bHn/bHnfgMrtnyk6FXqmugYAtn/OKoPHpmw8NqXjcSkbj035rwMcvE1ERERERJXGYEFERERERJWmuGChVqsxa9YsqNVquUsxC1veP+6bctny/tnyvgG2uX+2uE+mwmNTNh6b0vG4lI3HxngWH7xNRERERES2R3EdCyIiIiIisj4MFkREREREVGkMFkREREREVGkMFkREREREVGmKCxZLly5FQEAAXF1d0bZtWxw6dEjukow2e/ZsqFQqg49atWpJ64UQmD17Nvz8/ODm5obOnTvj7NmzMlZctoMHD6JXr17w8/ODSqVCfHy8wfry7ItWq8Wrr74Kb29veHh4oHfv3rh69aoF96Jsj9q/kSNHPnAun376aYNtrHH/YmJi8NRTT8HLyws+Pj7o27cvzp07Z7CNks9defZPqedu2bJlaNWqlfTAppCQECQkJEjrlXzeyssWrgN6lvxZzMzMxPDhw6HRaKDRaDB8+HBkZWUZbHP58mX06tULHh4e8Pb2xoQJE1BYWGiWfTdGTEwMVCoVJk2aJC2z5+Py119/YdiwYahRowbc3d3xxBNPIDU1VVpvr8dGp9Nh+vTpCAgIgJubGxo0aIC5c+eipKRE2sZej43FCAVZv369cHZ2FitWrBA///yzmDhxovDw8BCXLl2SuzSjzJo1Szz++OMiPT1d+rhx44a0fv78+cLLy0ts3rxZnD59WkRFRYnatWuLnJwcGasu3c6dO8Xbb78tNm/eLACIuLg4g/Xl2ZeXX35Z1KlTRyQmJooTJ06ILl26iNatWwudTmfhvXnQo/YvOjpa9OjRw+Bc3rp1y2Aba9y/5557TqxatUqcOXNGpKWliZ49e4p69eqJO3fuSNso+dyVZ/+Ueu62bdsmduzYIc6dOyfOnTsnpk2bJpydncWZM2eEEMo+b+VhK9cBPUv+LPbo0UMEBgaKI0eOiCNHjojAwEARGRkprdfpdCIwMFB06dJFnDhxQiQmJgo/Pz8xfvx4yxyMMhw7dkw89thjolWrVmLixInScns9Lrdv3xb169cXI0eOFEePHhUXLlwQe/bsEb///ru0jb0em3nz5okaNWqI7777Tly4cEFs3LhReHp6io8//ljaxl6PjaUoKli0a9dOvPzyywbLmjVrJt566y2ZKqqYWbNmidatW5e6rqSkRNSqVUvMnz9fWlZQUCA0Go347LPPLFRhxfz7D+/y7EtWVpZwdnYW69evl7b566+/hIODg9i1a5fFai+PsoJFnz59yvwapezfjRs3BABx4MABIYTtnbt/758QtnPuhBCiWrVq4osvvrC581YaW7kOlMVcP4s///yzACBSUlKkbZKTkwUA8euvvwoh7r2R4uDgIP766y9pm3Xr1gm1Wi2ys7PNt9MPkZubKxo3biwSExNFp06dpGBhz8flzTffFB06dChzvT0fm549e4rRo0cbLOvXr58YNmyYEMK+j42lKOZWqMLCQqSmpqJ79+4Gy7t3744jR47IVFXFnT9/Hn5+fggICMCgQYPw559/AgAuXLiAjIwMg/1Uq9Xo1KmT4vazPPuSmpqKoqIig238/PwQGBiomP1NSkqCj48PmjRpgv/85z+4ceOGtE4p+5ednQ0AqF69OgDbO3f/3j89pZ+74uJirF+/Hnl5eQgJCbG58/ZvtnYdKI25fhaTk5Oh0WgQHBwsbfP0009Do9EYbBMYGAg/Pz9pm+eeew5ardbgNhtLGjduHHr27Ilnn33WYLk9H5dt27YhKCgIL7zwAnx8fNCmTRusWLFCWm/Px6ZDhw7Yu3cvfvvtNwDAyZMncfjwYURERACw72NjKU5yF1BeN2/eRHFxMXx9fQ2W+/r6IiMjQ6aqKiY4OBhr1qxBkyZNcP36dcybNw+hoaE4e/astC+l7eelS5fkKLfCyrMvGRkZcHFxQbVq1R7YRgnnNTw8HC+88ALq16+PCxcuYMaMGejatStSU1OhVqsVsX9CCEyZMgUdOnRAYGAgANs6d6XtH6Dsc3f69GmEhISgoKAAnp6eiIuLQ4sWLaQLmi2ct9LY0nWgNOb8WczIyICPj88Dr+nj42Owzb9fp1q1anBxcZHl+K5fvx6pqak4fvz4A+vs+bj8+eefWLZsGaZMmYJp06bh2LFjmDBhAtRqNUaMGGHXx+bNN99EdnY2mjVrBkdHRxQXF+Pdd9/F4MGDpXoB+zw2lqKYYKGnUqkMPhdCPLDM2oWHh0v/btmyJUJCQtCwYUN89dVX0uBRW9hPvYrsi1L2NyoqSvp3YGAggoKCUL9+fezYsQP9+vUr8+usaf/Gjx+PU6dO4fDhww+ss4VzV9b+KfncNW3aFGlpacjKysLmzZsRHR2NAwcOSOtt4bw9jC39fryfuX8WS9u+IttYwpUrVzBx4kR8//33cHV1LXM7ezsuAFBSUoKgoCC89957AIA2bdrg7NmzWLZsGUaMGCFtZ4/HZsOGDfjmm28QGxuLxx9/HGlpaZg0aRL8/PwQHR0tbWePx8ZSFHMrlLe3NxwdHR9IeTdu3HggESqNh4cHWrZsifPnz0uzQ9nCfpZnX2rVqoXCwkJkZmaWuY2S1K5dG/Xr18f58+cBWP/+vfrqq9i2bRv2798Pf39/abmtnLuy9q80Sjp3Li4uaNSoEYKCghATE4PWrVtj0aJFNnPeymLL1wFz/yzWqlUL169ff+B1//77b4Nt/v06mZmZKCoqsvjxTU1NxY0bN9C2bVs4OTnByckJBw4cwOLFi+Hk5CTVY2/HBbj3u6pFixYGy5o3b47Lly8DsN//MwDwxhtv4K233sKgQYPQsmVLDB8+HJMnT0ZMTIxUL2Cfx8ZSFBMsXFxc0LZtWyQmJhosT0xMRGhoqExVmYZWq8Uvv/yC2rVrIyAgALVq1TLYz8LCQhw4cEBx+1mefWnbti2cnZ0NtklPT8eZM2cUt78AcOvWLVy5cgW1a9cGYL37J4TA+PHjsWXLFuzbtw8BAQEG65V+7h61f6VRyrkrjRACWq1W8eftUWzxOmCpn8WQkBBkZ2fj2LFj0jZHjx5Fdna2wTZnzpxBenq6tM33338PtVqNtm3bmn7nHyIsLAynT59GWlqa9BEUFIShQ4ciLS0NDRo0sMvjAgDt27d/YEri3377DfXr1wdgv/9nAODu3btwcDD809bR0VGabtaej43FmHt0uCnppxlcuXKl+Pnnn8WkSZOEh4eHuHjxotylGeW1114TSUlJ4s8//xQpKSkiMjJSeHl5Sfsxf/58odFoxJYtW8Tp06fF4MGDrXa62dzcXPHTTz+Jn376SQAQH330kfjpp5+kqR/Lsy8vv/yy8Pf3F3v27BEnTpwQXbt2tZqpLx+2f7m5ueK1114TR44cERcuXBD79+8XISEhok6dOla/f6+88orQaDQiKSnJYLrVu3fvStso+dw9av+UfO6mTp0qDh48KC5cuCBOnTolpk2bJhwcHMT3338vhFD2eSsPW7kO6FnyZ7FHjx6iVatWIjk5WSQnJ4uWLVuWOj1mWFiYOHHihNizZ4/w9/e3mukx758VSgj7PS7Hjh0TTk5O4t133xXnz58Xa9euFe7u7uKbb76RtrHXYxMdHS3q1KkjTTe7ZcsW4e3tLf7v//5P2sZej42lKCpYCCHEp59+KurXry9cXFzEk08+aTB9pFLo50x2dnYWfn5+ol+/fuLs2bPS+pKSEjFr1ixRq1YtoVarxTPPPCNOnz4tY8Vl279/vwDwwEd0dLQQonz7kp+fL8aPHy+qV68u3NzcRGRkpLh8+bIMe/Ogh+3f3bt3Rffu3UXNmjWFs7OzqFevnoiOjn6gdmvcv9L2CYBYtWqVtI2Sz92j9k/J52706NHS78CaNWuKsLAwKVQIoezzVl62cB3Qs+TP4q1bt8TQoUOFl5eX8PLyEkOHDhWZmZkG21y6dEn07NlTuLm5ierVq4vx48eLgoICc+2+Uf4dLOz5uGzfvl0EBgYKtVotmjVrJpYvX26w3l6PTU5Ojpg4caKoV6+ecHV1FQ0aNBBvv/220Gq10jb2emwsRSWEEJbojBARERERke1SzBgLIiIiIiKyXgwWRERERERUaQwWRERERERUaQwWRERERERUaQwWRERERERUaQwWRERERERUaQwWRERERERUaQwWRERERERUaQwWRERERERUaQwWRERERERUaQwWRERERERUaQwWRERERERUaf8PMwtV9ZrFESQAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x300 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "fig, axes = plt.subplots(1, 2, figsize=(8,3), tight_layout=True, sharey=\"all\")\n",
    "axes[0].plot(l, z, c='k')\n",
    "axes[1].plot(Karr, Harr , c='k')\n",
    "axes[0].set_title('cross-section')\n",
    "axes[1].set_title('K')\n",
    "# plt.savefig('tmp.png', bbox_inches=\"tight\", facecolor=\"white\")\n",
    "plt.show()                # 描画"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ffedbe58",
   "metadata": {},
   "source": [
    "さらに、流量$1000 \\mathrm{m^3/s}$、河床勾配$1/300$における$\\dfrac{\\tau}{\\rho g A} = \\dfrac{Q^2}{(\\sum K)^2}$と$H$の関係は以下となる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "752b8f8a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAHECAYAAAAXlgHcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA260lEQVR4nO3de1xUdeL/8fcAgoiAikoqaJqmJuIFvJalZZbalmFeyoxqsyyt1K02a2t37buxPqxd27y3bVpmaipeSjNM0kotZSUtXdO8a2qgMYA63Ob3hz9mY72BDPM5M/N6Ph7n4XDmMLzng8Lbcz7nHJvT6XQKAADAggJMBwAAALgYigoAALAsigoAALAsigoAALAsigoAALAsigoAALAsigoAALCsINMBKqOkpERHjx5VeHi4bDab6TgAAKAcnE6ncnNz1bBhQwUEXHqfiVcXlaNHjyo2NtZ0DAAAcAUOHTqkmJiYS27j1UUlPDxc0rk3GhERYTiNb8vPlxo2PPf46FEpLMxsHgCA97Lb7YqNjXX9Hr8Ury4qpYd7IiIiKCpVLDDwv48jIigqAIDKK8+0DSbTAgAAy6KoAAAAy6KoAAAAy6KoAAAAy6KoAAAAy6KoAAAAy6KoAAAAy6KoAAAAy6KoAAAAy6KoAAAAyzJeVI4cOaL7779fUVFRqlGjhtq3b6+MjAzTsQAAgAUYvdfPqVOndP3116tXr15atWqV6tevrx9//FG1atUyGQsAAFiE0aIyceJExcbG6p133nGtu/rqq80F+pVffvlFdrtdjRs3Nh0FAAC/ZfTQz/Lly5WYmKhBgwapfv366tChg956662Lbu9wOGS328ssVWHGjBmKiorS+PHjq+T1AQBA+RgtKnv37tX06dPVokULrV69WiNHjtRTTz2ld99994Lbp6SkKDIy0rXExsZWSa7rrrtOJSUlSktLU0lJSZV8DQAAcHk2p9PpNPXFg4ODlZiYqA0bNrjWPfXUU9q8ebM2btx43vYOh0MOh8P1sd1uV2xsrHJychQREeG2XAUFBYqKilJeXp4yMjLUsWNHt722t8rPl2rWPPc4L08KCzObBwDgvex2uyIjI8v1+9voHpUGDRrouuuuK7OudevWOnjw4AW3DwkJUURERJmlKgQHB+vmm2+WJH366adV8jUAAMDlGS0q119/vXbt2lVm3Q8//KAmTZoYSvRfffr0kURRAQDAJKNFZezYsdq0aZNeffVV7dmzR/PmzdOsWbM0atQok7Ek/beofPnll8rLyzOcBgAA/2S0qHTq1Empqan64IMPFBcXp1deeUWTJ0/WsGHDTMaSJDVv3lxNmzZVYWGh1q1bZzoOAAB+yfiVae+44w5t375dZ8+e1c6dOzVixAjTkSRJNpuNwz8AABhmvKhY2W233SZJWr16teEkAAD4J4rKJfTq1UuBgYHatWuXDhw4YDoOAAB+h6JyCbVq1VKXLl0kcfgHAAATKCqXUXr4h6ICAIDnUVQuo3RC7Zo1a1RcXGw4DQAA/oWichmJiYmqVauWfvnlF23evNl0HAAA/ApF5TKCgoLUu3dvSZz9AwCAp1FUyqF0nsonn3xiOAkAAP6FolIOffv2lSR9/fXXysrKMpwGAAD/QVEph0aNGqldu3ZyOp0c/gEAwIMoKuXUr18/SdLHH39sOAkAAP6DolJO/fv3l3RungqnKQMA4BkUlXLq0qWLateurVOnTunrr782HQcAAL9AUSmnoKAg19k/K1euNJwGAAD/QFGpgNJ5KhQVAAA8g6JSAbfddptsNpu2bt2qo0ePmo4DAIDPo6hUQP369dWpUydJXPwNAABPoKhUEId/AADwHIpKBZUWlU8//VQFBQWG0wAA4NsoKhWUkJCg+vXrKzc3V1999ZXpOAAA+DSKSgUFBAS47v3D4R8AAKoWReUKlB7++eijjwwnAQDAt1FUrkCfPn0UFBSk//znP9q9e7fpOAAA+CyKyhWoVauWevbsKUlasWKF2TAAAPgwisoVuvPOOyVJy5YtM5wEAADfRVG5Qr/5zW8kSV9++aWys7MNpwEAwDdRVK7Q1Vdfrfj4eJWUlHD2DwAAVYSiUgmlh3+WL19uOAkAAL6JolIJd911l6Rz9/1xOByG0wAA4HsoKpXQsWNHNWzYUHl5eUpPTzcdBwAAn0NRqYSAgADXpFoO/wAA4H4UlUr69TwVp9NpOA0AAL6FolJJN998s8LCwnTkyBFt3brVdBwAAHwKRaWSqlevrttuu00SF38DAMDdKCpuwGnKAABUDYqKG/Tv318BAQHKzMzUgQMHTMcBAMBnUFTcoG7duurRo4ckacmSJYbTAADgOygqbpKUlCSJogIAgDtRVNzk7rvvliR99dVXOnbsmOE0AAD4BoqKm8TGxqpz585yOp1aunSp6TgAAPgEioobDRw4UBKHfwAAcBeKihuVzlNJT0/XyZMnDacBAMD7UVTcqHnz5oqPj1dRUZFWrFhhOg4AAF6PouJmpXtVFi9ebDgJAADej6LiZqXzVD799FPl5uYaTgMAgHczWlT+9Kc/yWazlVmuuuoqk5EqrU2bNmrRooUcDodWrlxpOg4AAF7N+B6VNm3a6KeffnIt27dvNx2pUmw2G2f/AADgJsaLSlBQkK666irXUq9ePdORKq10nsrHH3+sM2fOGE4DAID3Ml5Udu/erYYNG6pp06YaOnSo9u7de9FtHQ6H7HZ7mcWKEhMTFRsbq/z8fH366aem4wAA4LWMFpUuXbro3Xff1erVq/XWW2/p2LFj6t69u7Kzsy+4fUpKiiIjI11LbGyshxOXz68P/3z44YeG0wAA4L1sTqfTaTpEqfz8fF1zzTV67rnnNG7cuPOedzgccjgcro/tdrtiY2OVk5OjiIgIT0a9rI0bN6p79+6qWbOmTpw4odDQUNORKiU/X6pZ89zjvDwpLMxsHgCA97Lb7YqMjCzX72/jh35+LSwsTG3bttXu3bsv+HxISIgiIiLKLFbVtWtXNW7cWHl5eVq1apXpOAAAeCVLFRWHw6GdO3eqQYMGpqNUms1m0+DBgyVJCxYsMJwGAADvZLSoPPPMM1q3bp327dunr7/+Wvfcc4/sdruSk5NNxnKbIUOGSJI++ugj5efnG04DAID3MVpUDh8+rHvvvVctW7ZUUlKSgoODtWnTJjVp0sRkLLdJSEhQs2bNdPr0aX300Uem4wAA4HWMFpX58+fr6NGjKigo0JEjR7R48WJdd911JiO5lc1mc+1V4fAPAAAVZ6k5Kr6otKisXLnSstd9AQDAqigqVSw+Pl7XXnutHA6Hli9fbjoOAABehaJSxTj8AwDAlaOoeEBpUVm9erVOnTplOA0AAN6DouIBbdq0UZs2bVRYWKhly5aZjgMAgNegqHhI6V6V+fPnG04CAID3oKh4yNChQyVJa9as0fHjxw2nAQDAO1BUPKRFixbq1KmTiouLmVQLAEA5UVQ86P7775ckvf/++4aTAADgHSgqHjRkyBAFBgbqm2++uegdogEAwH9RVDwoOjpat956qyT2qgAAUB4UFQ8bNmyYJGnu3LlyOp2G0wAAYG0UFQ8bMGCAatSooR9//FHffPON6TgAAFgaRcXDatasqQEDBkg6t1cFAABcHEXFgNKzfxYsWKDCwkLDaQAAsC6KigG33nqr6tWrp59//llr1qwxHQcAAMuiqBgQFBTkulIth38AALg4ioohpWf/LF26VHl5eYbTAABgTRQVQzp37qzmzZvr9OnTSk1NNR0HAABLoqgYYrPZNHz4cEnSnDlzDKcBAMCaKCoGPfDAA5KktWvX6sCBA4bTAABgPRQVg66++mr16tVLTqdT7777ruk4AABYDkXFsAcffFCSNHv2bC6pDwDA/6CoGDZw4EDVrFlTe/fu1Zdffmk6DgAAlkJRMSwsLEyDBw+WJL3zzjuG0wAAYC0UFQsoPfyzcOFCrqkCAMCvUFQs4IYbbtA111yj/Px8LVmyxHQcAAAsg6JiATabrcykWgAAcA5FxSIeeOAB2Ww2paena9++fabjAABgCRQVi2jcuLFuueUWSeKaKgAA/H8UFQspPfwzZ84clZSUmA0DAIAFUFQs5O6771Z4eLj27dundevWmY4DAIBxFBULqVGjhu69915J0ltvvWU4DQAA5lFULGbEiBGSpMWLFys7O9twGgAAzKKoWExCQoLat2+vgoICvffee6bjAABgFEXFYmw2m2uvyltvvcWNCgEAfo2iYkHDhg1TaGioduzYoY0bN5qOAwCAMRQVC4qMjHTdqJBJtQAAf0ZRsajSwz8LFixQTk6O4TQAAJhBUbGo7t27q3Xr1jpz5ozmzZtnOg4AAEZQVCzKZrPp0UcflcThHwCA/6KoWNjw4cMVHBysrVu3KiMjw3QcAAA8jqJiYVFRURo4cKAk9qoAAPwTRcXiSifVzps3T3l5eYbTAADgWRQVi+vZs6eaN2+u3NxczZ8/33QcAAA8yjJFJSUlRTabTWPGjDEdxVJsNpsee+wxSdLUqVO5Ui0AwK9Yoqhs3rxZs2bNUnx8vOkolvTwww+revXqyszM5Eq1AAC/Yryo5OXladiwYXrrrbdUu3Zt03EsqU6dOrrvvvskSVOmTDGcBgAAzzFeVEaNGqX+/furd+/el93W4XDIbreXWfzFqFGjJEmLFi3SsWPHDKcBAMAzjBaV+fPnKyMjQykpKeXaPiUlRZGRka4lNja2ihNaR8eOHdWtWzcVFhZyqjIAwG8YKyqHDh3S008/rffff1/Vq1cv1+eMHz9eOTk5ruXQoUNVnNJaRo8eLUmaOXOmCgsLDacBAKDq2ZyGTiNZunSp7r77bgUGBrrWFRcXy2azKSAgQA6Ho8xzF2K32xUZGamcnBxFRERUdWTjHA6HGjdurBMnTujDDz/UPffc47GvnZ8v1ax57nFenhQW5rEvDQDwMRX5/W1sj8ott9yi7du3KzMz07UkJiZq2LBhyszMvGxJ8UchISGu+/9MnTrVcBoAAKqesaISHh6uuLi4MktYWJiioqIUFxdnKpblPfbYYwoMDNTnn3+u7777znQcAACqlPGzflAxMTExGjBggCRp2rRpZsMAAFDFjM1RcQd/m6NSKj09XTfffLPCwsJ05MgRRUZGVvnXZI4KAMBdvGKOCq5cz549dd111yk/P1/vvvuu6TgAAFQZiooXstlsrgvAcf8fAIAvo6h4qeHDhys8PFy7du3SZ599ZjoOAABVgqLipcLDw5WcnCyJ+/8AAHwXRcWLlR7+WbFihQ4cOGA4DQAA7kdR8WKtWrXSLbfcopKSEs2YMcN0HAAA3I6i4uVK7//zz3/+U2fPnjWcBgAA96KoeLk77rhDsbGxysrK0ocffmg6DgAAbkVR8XJBQUF6/PHHJTGpFgDgeygqPuCRRx5RcHCwvvnmG23evNl0HAAA3Iai4gPq1aunIUOGSOKuygAA30JR8RGlpyrPnz9fWVlZhtMAAOAeFBUf0blzZyUmJsrhcOjtt982HQcAALegqPiIX9//Z/r06SouLjacCACAyqOo+JAhQ4YoKipKBw4c0Mcff2w6DgAAlUZR8SGhoaH67W9/K4lJtQAA30BR8TGPP/64bDabPv30U/3www+m4wAAUCkUFR9z9dVX64477pAkTZs2zXAaAAAqh6Lig0rv//POO+8oLy/PcBoAAK4cRcUH9e7dWy1atJDdbtfcuXNNxwEA4IpRVHxQQECA61TlqVOnyul0Gk4EAMCVoaj4qOTkZIWFhem7777T+vXrTccBAOCKUFR8VK1atTR8+HBJ3FUZAOC9KCo+rPTwT2pqqg4fPmw4DQAAFUdR8WFxcXHq2bOniouLNXPmTNNxAACoMIqKjys9VXnWrFlyOByG0wAAUDEUFR931113KSYmRidOnNCiRYtMxwEAoEIoKj4uKChII0eOlMSkWgCA96Go+IERI0YoODhYmzZt0pYtW0zHAQCg3CgqfqB+/foaPHiwJO6qDADwLhQVP1E6qfaDDz5QVlaW4TQAAJQPRcVPdO7cWYmJiXI4HHr77bdNxwEAoFwoKn7CZrO59qpMmzZNRUVFhhMBAHB5FBU/MmTIEEVFRengwYP66KOPTMcBAOCyKCp+pHr16hoxYoQkTlUGAHgHioqfGTlypAICAvTZZ59p586dpuMAAHBJFBU/06RJE915552SOFUZAGB9FBU/VDqpds6cObLb7YbTAABwcRQVP3TzzTerVatWysvL07vvvms6DgAAF0VR8UO/PlV5ypQpcjqdhhMBAHBhFBU/9cADDyg8PFy7du3SZ599ZjoOAAAXRFHxU+Hh4UpOTpbEqcoAAOuiqPixUaNGSZJWrFih/fv3mw0DAMAFUFT8WKtWrdS7d2+VlJRoxowZpuMAAHAeioqfK51U+89//lNnzpwxnAYAgLIoKn7ujjvuUOPGjZWdna0FCxaYjgMAQBlGi8r06dMVHx+viIgIRUREqFu3blq1apXJSH4nMDBQTzzxhCTpzTff5FRlAICl2JwV+M1U3quYRkRElGu7FStWKDAwUM2bN5d07kqpkyZN0tatW9WmTZty5YmMjFROTk65vybOl5WVpZiYGDkcDm3cuFFdu3Y9b5v8fKlmzXOP8/KksDAPhwQA+IyK/P6uUFEJCAiQzWa76PNOp1M2m03FxcXlT/s/6tSpo0mTJum3v/3tec85HA45HA7Xx3a7XbGxsRQVN3jooYc0e/ZsDRs2THPnzj3veYoKAMBdKlJUgirywunp6a7HTqdT/fr10z//+U81atToypL+SnFxsT788EPl5+erW7duF9wmJSVFf/7znyv9tXC+0aNHa/bs2Vq4cKFef/11RUdHm44EAEDF9qj8r/DwcH377bdq1qzZFQfYvn27unXrprNnz6pmzZqaN2+e+vXrd8Ft2aNStbp166ZNmzbplVde0R/+8Icyz7FHBQDgLhXZo2L8rJ+WLVsqMzNTmzZt0uOPP67k5GTt2LHjgtuGhIS4Jt6WLnCfJ598UpI0Y8YMFRYWGk4DAIAFikpwcLCaN2+uxMREpaSkqF27dnrjjTdMx/JL99xzj6Kjo3XkyBEtXbrUdBwAACpfVC41ufZKOJ3OMod34DnBwcF69NFHJXH/HwCANVRoMm1SUlKZj8+ePauRI0cq7H8mLCxZsqRcr/fCCy+ob9++io2NVW5urubPn6/PP/9cn3zySUViwY0ee+wxvfrqq1q/fr22bdum+Ph405EAAH6sQkUlMjKyzMf3339/pb748ePHNXz4cP3000+KjIxUfHy8PvnkE916662Vel1cuUaNGikpKUkffvihpkyZolmzZpmOBADwY5U668c0LvhWNdavX6+bbrpJoaGhOnTokKKiojjrBwDgNl511g+sp0ePHurQoYPOnDnDHhUAgFEUFZzHZrNp7Nixks5Nqi0oKDCcCADgrygquKDBgwfrqquu0tGjR7Vo0SLTcQAAfoqiggsKCQlx3VX573//O3dVBgAYQVHBRY0cOVIhISHasmWLNm3aZDoOAMAPUVRwUfXq1XOdgj516lTDaQAA/oiigksaM2aMJGnZsmVmgwAA/BJFBZcUFxen3r17y+ksMR0FAOCHKCq4rNK9KgAAeBpFBZfVt29ftWrV2nQMAIAfoqjgsgICAlwXgJPO3YwSAABPoKigXAYNGuR6PG/ePINJAAD+hKKCcgkODnY9njx5soqLiw2mAQD4C4oKKmzv3h+VmppqOgYAwA9QVHBF/vrXv3JZfQBAlaOooMKqVw9VRkaG1q5dazoKAMDHUVRQYcnJyZLO7VUBAKAqUVRQYU899ZQCAwO1Zs0affPNN6bjAAB8GEUFFdakSRPXzQonTJhgOA0AwJdRVHBFXnzxRQUEBOjjjz/Wli1bTMcBAPgoigquSIsWLTRs2DBJ7FUBAFQdigquWOlelRUrVujf//636TgAAB9EUcEVa9mypYYOHSpJeuWVVwynAQD4IooKKuUPf/iDbDabli5dqm+//dZ0HACAj6GooFJat26tIUOGSGKuCgDA/SgqqLSXXnpJNptNS5YsYa8KAMCtKCqotOuuu06DBw+WdO5QEAAA7kJRgVtMmDBBgYGB+uijj7RhwwbTcQAAPoKiAre49tpr9dBDD0mSXnjhBe6sDABwC4oK3Obll19WSEiI1q1bp7S0NNNxAAA+gKICt4mNjdUTTzwhib0qAAD3oKjArcaPH6+aNWsqIyNDixcvNh0HAODlKCpwq3r16mncuHGSzp22XFRUZDgRAMCbUVTgdr/73e9Up04d/ec//9GcOXNMxwEAeDGKCtwuIiJCL774oqRz11XJy8sznAgA4K0oKqgSo0aNUrNmzXTs2DG99tprpuMAALwURQVVIiQkRBMnTpQkTZo0SUePHjWcCADgjSgqqDIDBw5U9+7ddfr0aS6tDwC4IhQVVBmbzabXX39dkjR79mxlZmaaDQQA8DoUFVSprl27asiQIXI6nXrmmWe4CBwAoEIoKqhyKSkpCg4O1meffaaPP/7YdBwAgBehqKDKNW3aVE8//bQkacyYMTp79qzhRAAAb0FRgUf84Q9/UIMGDfTjjz/qb3/7m+k4AAAvQVGBR0RERGjSpEmSpL/85S86dOiQ4UQAAG9AUYHH3Hfffbrhhht0+vRpPfPMM6bjAAC8gNGikpKSok6dOik8PFz169fXgAEDtGvXLpORUIVsNpumTJmigIAALVy4UGvXrjUdCQBgcUaLyrp16zRq1Cht2rRJaWlpKioqUp8+fZSfn28yFqpQu3bt9Pjjj0uSnnrqKRUWFhpOBACwMpvTQhe2+Pnnn1W/fn2tW7dON95443nPOxwOORwO18d2u12xsbHKyclRRESEJ6P6nfx8qWbNc4/z8qSwsCt/rZMnT6ply5bKysrS3/72N40dO9Y9IQEAXsFutysyMrJcv78tNUclJydHklSnTp0LPp+SkqLIyEjXEhsb68l4cJM6deooJSVFkvTyyy8zsRYAcFGW2aPidDp111136dSpU/riiy8uuA17VMxx5x4VSSopKVGPHj20YcMG3XXXXVq6dGmlMwIAvINX7lEZPXq0tm3bpg8++OCi24SEhCgiIqLMAu8UEBCgmTNnKigoSMuWLVNqaqrpSAAAC7JEUXnyySe1fPlypaenKyYmxnQceEhcXJx+//vfSzpXVEsP/QEAUMpoUXE6nRo9erSWLFmitWvXqmnTpibjwIAXX3xRzZs319GjR/Xiiy+ajgMAsBijRWXUqFGaO3eu5s2bp/DwcB07dkzHjh3TmTNnTMaCB4WGhmrGjBmSpGnTpmnTpk2GEwEArMRoUZk+fbpycnLUs2dPNWjQwLUsWLDAZCx42C233KLk5GQ5nU6NGDFCBQUFpiMBACzC+KGfCy0PPvigyVgw4LXXXlPdunX13Xff6f/+7/9MxwEAWIQlJtMCdevW1bRp0yRJr776qjIyMgwnAgBYAUUFljFo0CANHjxYxcXFevDBB8tcMwcA4J8oKrCUKVOmqF69evruu+80YcIE03EAAIZRVGAp9erV0/Tp0yVJf/3rX7V582bDiQAAJlFUYDkDBw7U0KFDVVJSouTkZE5XBwA/RlGBJU2ZMkXR0dHauXOnnn/+edNxAACGUFRgSVFRUXrnnXckSf/4xz+0atUqw4kAACZQVGBZffv21ZNPPilJevDBB3XixAnDiQAAnkZRgaVNnDhRcXFxOnHihB566CE5nU7TkQAAHkRRgaWFhoZq3rx5CgkJ0cqVK10XhQMA+AeKCiyvbdu2mjhxoiTpmWee0ffff284EQDAUygq8ApPPfWUbr/9dp09e1aDBg1Sfn6+6UgAAA+gqMAr2Gw2zZ49Ww0aNNDOnTv1+OOPM18FAPwARQVeIzo6WvPnz1dAQIDee+89vf3226YjAQCqGEUFXuXGG2/UX/7yF0nS6NGjlZmZaTYQAKBKUVTgdZ577jn1799fDodDgwYNUk5OjulIAIAqQlGB1wkICNCcOXPUuHFj7dmzRw8//DDzVQDAR1FU4JWioqK0cOFCVatWTUuWLHEdDgIA+BaKCrxWly5dXBeAe+mll7R8+XLDiQAA7kZRgVd75JFH9MQTT0iS7r//fu3cudNwIgCAO1FU4PUmT56sG2+8Ubm5ubrrrrv0yy+/mI4EAHATigq8XrVq1bRo0SI1btxYu3fv1r333qvi4mLTsQAAbkBRgU+oV6+eli5dqtDQUH3yySf63e9+ZzoSAMANKCrwGR06dNDs2bMlSW+88YbefPNNs4EAAJVGUYFPGTx4sFJSUiRJY8aM0YoVKwwnAgBUBkUFPuf3v/+9HnnkEZWUlGjo0KHKyMgwHQkAcIUoKvA5NptN06ZNU58+fXT69GndcccdOnjwoOlYAIArQFGBT6pWrZoWLlyouLg4HTt2TLfffruys7NNxwIAVBBFBT4rMjJSH3/8sRo1aqSdO3eqf//+ysvLMx0LAFABFBX4tMaNG+vTTz9VnTp19PXXXyspKUkOh8N0LABAOVFU4POuu+46rVy5UmFhYUpLS9MDDzzABeEAwEtQVOAXunTpotTUVNfclVGjRsnpdJqOBQC4DIoK/Matt96q999/XzabTTNnztSzzz5LWQEAi6OowK8MGjRIM2fOlCS9/vrrGj9+PGUFACyMogK/M2LECE2ZMkWSNHHiRL300kuUFQCwKIoK/NKoUaM0efJkSdJf/vIXTZgwwWwgAMAFUVTgt55++mm9/vrrkqQ//elPeuWVV9izAgAWQ1GBXxs3bpwmTpwoSXr55ZeZswIAFkNRgd977rnnXHtWJk6cqNGjR6ukpMRwKgCARFEBJJ3bszJz5kzXDQ0ffvhhFRUVmY4FAH6PogL8f48++qjee+89BQYGas6cORo6dKgKCgpMxwIAv0ZRAX5l2LBhWrRokYKDg7V48WL1799fdrvddCwA8FsUFeB/DBgwQCtWrFBYWJjWrFmjHj166OjRo6ZjAYBfoqgAF9CnTx+tX79e0dHR2rZtm7p27arvv//edCwA8DsUFeAiOnbsqI0bN6ply5Y6dOiQbrjhBq1bt850LADwK0aLyvr16/Wb3/xGDRs2lM1m09KlS03GAc7TtGlTffXVV+revbt++eUX9enTR/PmzTMdCwD8htGikp+fr3bt2rnuuwJYUVRUlNasWaOkpCQVFBRo2LBheuGFF7jWCgB4gM1pkctw2mw2paamasCAARfdxuFwyOFwuD622+2KjY1VTk6OIiIiPJDSf+XnSzVrnnuclyeFhZnNY0JxcbFefPFF15Vs77zzTs2dO1fh4eGGkwGAd7Hb7YqMjCzX72+vmqOSkpKiyMhI1xIbG2s6EvxIYGCg/vrXv+q9995TSEiIli9fru7du2vfvn2mowGAz/KqojJ+/Hjl5OS4lkOHDpmOBD90//33a926dWrQoIG+++47derUSWlpaaZjAYBP8qqiEhISooiIiDILYEKXLl20efNmJSYmKjs7W7fddpsmTJjAvBUAcDOvKiqAlTRq1Ejr16/XiBEj5HQ69cc//lH9+vVTVlaW6WgA4DMoKkAlhIaGatasWZo9e7ZCQ0O1evVqdejQQZs2bTIdDQB8gtGikpeXp8zMTGVmZkqS9u3bp8zMTB08eNBkLKDCkpOTtWnTJrVo0UKHDx/WjTfeqEmTJnEoCAAqyWhR2bJlizp06KAOHTpIksaNG6cOHTro5ZdfNhkLuCLx8fHasmWL7rnnHhUWFuq5557TrbfeqsOHD5uOBgBey2hR6dmzp5xO53nL7NmzTcYCrlhERIQWLlyoWbNmqUaNGlq7dq3i4+O1aNEi09EAwCsxRwVwM5vNphEjRmjr1q1KSEjQqVOnNGjQID388MPKzc01HQ8AvApFBagi1157rTZs2KDx48fLZrPpnXfeUbt27ZSenm46GgB4DYoKUIWCg4P16quvKj09XY0bN9a+fft088036/HHH5fdbjcdDwAsj6ICeMBNN92k7du3a+TIkZKkGTNmKC4uTqtWrTKcDACsjaICeEhERISmT5+u9PR0NWvWTIcOHVK/fv2UnJysn3/+2XQ8ALAkigrgYT179tS2bds0duxY2Ww2vfvuu2rZsqVmzpyp4uJi0/EAwFIoKoABYWFh+tvf/qYNGzaoffv2OnXqlEaOHKlu3bopIyPDdDwAsAyKCmBQ165dtXnzZv3jH/9QRESENm/erE6dOumJJ57QyZMnTccDAOMoKoBhQUFBevLJJ/Wf//xHw4YNk9Pp1PTp09W8eXNNnjxZBQUFpiMCgDEUFcAiGjRooLlz52rt2rVq06aNTp06pbFjx6pNmzZasmSJnE6n6YgA4HEUFcBievXqpczMTM2aNUvR0dHas2ePBg4cqJtuukmbN282HQ8APIqiAlhQUFCQRowYod27d+ull15SaGiovvjiC3Xu3FmDBw/Wzp07TUcEAI+gqAAWFh4ergkTJuiHH37Qgw8+KJvNpg8//FBxcXEaPny49uzZYzoiAFQpigrgBWJiYvTOO+/o22+/1d13362SkhLNnTtXrVq10iOPPKIDBw6YjggAVYKiAniRtm3basmSJcrIyFD//v1VXFyst99+W82bN9dDDz3EISEAPoeiAnihjh076qOPPtKGDRvUu3dvFRUVafbs2WrTpo2SkpKYdAvAZ1BUAC/WrVs3paWl6euvv9aAAQPkdDqVmpqqzp07q3fv3kpLS+O0ZgBejaIC+IDOnTsrNTVVO3bsUHJysoKCgvTZZ5+pT58+iouL08yZM5Wfn286JgBUGEUF8CGtW7fW7NmztWfPHj311FOqWbOmduzYoZEjRyomJkbPPvus9u/fbzomAJQbRQXwQU2aNNEbb7yhw4cP6+9//7uaNWumX375Ra+99pquueYaJSUlKS0tTSUlJaajAsAlUVQAHxYZGakxY8bohx9+0IoVK9S7d2+VlJQoNTVVffr0UbNmzTRhwgQdOnTIdFQAuCCb04tn2tntdkVGRionJ0cRERGm4/i0/HypZs1zj/PypLAws3lw5Xbs2KGpU6fq/fffV05OjiTJZrPp9ttv1yOPPKI77rhDwcHBhlMC8GUV+f1NUUG5UFR8z5kzZ7R48WK9/fbb+vzzz13r69Wrp2HDhmnYsGFKSEiQzWYzFxKAT6KowO0oKr5t9+7d+te//qXZs2fr2LFjrvXXXnut7rvvPt13331q0aKFwYQAfAlFBW5HUfEPhYWF+uSTT/T+++9r+fLlOnPmjOu5Tp066d5771VSUpKaNGliMCUAb0dRgdtRVPxPbm6uli1bpvfff19paWkqLi52PZeQkKCkpCQlJSWpVatWBlMC8EYUFbgdRcW/nThxQgsXLtSiRYv0xRdflDmtuXXr1kpKStKAAQPUsWNHBQRwMiGAS6OowO0oKih14sQJLV++XEuWLNGaNWtUWFjoei46Olq33367+vXrp1tvvVW1a9c2mBSAVVFU4HYUFVxITk6OVq5cqcWLF2v16tXKy8tzPRcYGKhu3bqpb9++6tevn9q1a8cZRAAkUVRQBSgquJyCggJ9+eWXWrVqlVauXKkdO3aUeT46Olo9e/ZUr1691KtXL7Vo0YLiAvgpigrcjqKCijpw4ICrtHz22Wc6ffp0mecbNmzoKi29evVS06ZNKS6An6CowO0oKqgMh8Ohr7/+Wunp6UpPT9fGjRtVUFBQZpvGjRvrxhtvVPfu3dW9e3fFxcUpMDDQUGIAVYmiArejqMCdzpw5ow0bNriKyzfffKOioqIy29SsWVNdu3ZV9+7d1a1bN3Xt2lW1atUyExiAW1FU4HYUFVSlvLw8bdiwwbVs2rRJubm5523XvHlzJSYmKiEhQYmJierQoYMiIyMNJAZQGRQVuB1FBZ5UXFys77//Xhs3bnSVlz179lxw2xYtWrjKS7t27dS2bVtFR0d7ODGAiqCowO0oKjAtOztb//73v7VlyxZlZGRoy5YtOnDgwAW3rV+/vtq2bVtmadOmjWrUqOHh1AAuhKICt6OowIqysrKUkZHhWrZv3649e/boQj/WbDabmjdvrjZt2qhly5ZllqioKAPpAf9FUYHbUVTgLfLz87Vjxw5t3769zHLixImLfk5UVFSZ4nLttdeqZcuWuuaaaxQSEuLB9IB/oKjA7Sgq8HbHjx/X9u3btWPHDu3atUs//PCDdu3apUOHDl30c2w2mxo1aqSmTZvq6quvPu/PmJgYBQUFefBdAL6BogK3o6jAV+Xn52v37t3atWtXmQKza9euC5559GuBgYGKjY29YIFp1KiRGjVqpDD+sQDnoajA7Sgq8DdOp1MnTpzQ/v37tX//fu3bt6/Mn/v37z/vonUXUqtWLVdpKV1Ki0zDhg0VHR2t+vXrq1q1ah54V4A1UFTgdhQVoKySkhIdO3bsgiXmyJEjOnz4cJmbNF5OnTp1FB0dfcHlqquucj2Oiori7CV4PYoK3I6iAlSc3W7XkSNHXMXlfx//9NNPOnHihIqLiyv0utWrV1dUVJTq1KmjqKioiy6/fr527drMp4FlVOT3N39rAaCKREREKCIiQq1bt77oNiUlJTp58qSOHTum48ePX3YpKirS2bNnXaWnIsLDwxUREaHIyEhFRkZe9vGF1lWvXr2ywwJUCEUFAAwKCAhQ3bp1VbduXcXFxV1yW6fTqdzcXGVnZ+vkyZPKzs4+b7nQ+pycHElSbm6ucnNzK1xwfi04ONhVwMLCwlSjRg3Xn79+fKF15XmeG1HifxkvKtOmTdOkSZP0008/qU2bNpo8ebJ69OhhOhYAWI7NZnOVhKZNm5b784qKinTy5EnZ7Xbl5OQoJyenwo9Lz4AqKChQVlaWsrKyquQ9hoSEXLDIhIaGKiQkRCEhIapevXqlH1erVk3BwcGqVq3aeUvpepvNViXvERVjtKgsWLBAY8aM0bRp03T99ddr5syZ6tu3r3bs2KHGjRubjAYAPiMoKEj169dX/fr1r/g1SkpKlJub6yowdrtdp0+fVn5+vk6fPn3Rx+V9vpTD4ZDD4dCpU6fc8dYrJTAw8IIF5mLF5lLrAwMDq2QJCgqqstcuXapXr67Q0FBj3wejk2m7dOmijh07avr06a51rVu31oABA5SSknLZz2cyrecwmRZAVXE6nTp79uwlS82ZM2dcJebs2bOVflxYWKiCggIVFhaqsLBQRUVFpofBsoYOHaoPPvjAra/pFZNpCwoKlJGRoeeff77M+j59+mjDhg0X/JzSv2il7HZ7lWYEAFQ9m82m0NBQo/9rdzqdKioqKlNefr1cbP2lnvv1+uLiYrcvRUVFVfK6pUsp0/OGjBWVrKwsFRcXn3c79ujoaB07duyCn5OSkqI///nPnogHAPAjNpvNdZgG55SUlFhiT1OA6QD/O1nJ6XRedALT+PHjXRO7cnJyLnmPDgAAcOUCAgIUHBys4OBgozmM7VGpW7euAgMDz9t7cuLEifP2spQqnbUNAAD8g7E9KsHBwUpISFBaWlqZ9WlpaerevbuhVAAAwEqMnp48btw4DR8+XImJierWrZtmzZqlgwcPauTIkSZjAQAAizBaVIYMGaLs7GxNmDBBP/30k+Li4rRy5Uo1adLEZCwAAGAR3JQQ5cJ1VAAA7lKR39/Gz/oBAAC4GIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLIoKAACwLKNXpq2s0mvV2e12w0l8X37+fx/b7VJxsbksAADvVvp7uzzXnPXqopKbmytJio2NNZzEvzRsaDoBAMAX5ObmKjIy8pLbePUl9EtKSnT06FGFh4fLZrOZjmOU3W5XbGysDh06xO0EPIyxN4NxN4NxN8PXxt3pdCo3N1cNGzZUQMClZ6F49R6VgIAAxcTEmI5hKRERET7xl9gbMfZmMO5mMO5m+NK4X25PSikm0wIAAMuiqAAAAMuiqPiIkJAQ/fGPf1RISIjpKH6HsTeDcTeDcTfDn8fdqyfTAgAA38YeFQAAYFkUFQAAYFkUFQAAYFkUFQAAYFkUFQAAYFkUFT+Um5urTp06qX379mrbtq3eeust05H8xt13363atWvrnnvuMR3FpzHOnsfPFbOCgoLUvn17tW/fXo888ojpOG7F6cl+qLi4WA6HQzVq1NDp06cVFxenzZs3KyoqynQ0n5eenq68vDzNmTNHixYtMh3HZzHOnsfPFbPq1q2rrKws0zGqBHtU/FBgYKBq1KghSTp79qyKi4vLdattVF6vXr0UHh5uOobPY5w9j58rqCoUFQtav369fvOb36hhw4ay2WxaunTpedtMmzZNTZs2VfXq1ZWQkKAvvviiQl/jl19+Ubt27RQTE6PnnntOdevWdVN67+WJccfl8X0wwx3jzs+VK+OOsbfb7UpISNANN9ygdevWeSi5Z1BULCg/P1/t2rXTlClTLvj8ggULNGbMGL344ovaunWrevToob59++rgwYOubRISEhQXF3fecvToUUlSrVq19O2332rfvn2aN2+ejh8/7pH3ZmWeGHdcnju+D6g4d4w7P1eujDvGfv/+/crIyNCMGTP0wAMPyG63eyp+1XPC0iQ5U1NTy6zr3Lmzc+TIkWXWtWrVyvn8889f0dcYOXKkc+HChVca0SdV5binp6c7Bw4cWNmIfqEy3wfG+cq54+8/P1eujDvG/vbbb3du3ry5qiJ6HHtUvExBQYEyMjLUp0+fMuv79OmjDRs2lOs1jh8/7mrbdrtd69evV8uWLd2e1Ze4Y9xReXwfzCjPuPNzpWqUZ+xPnTolh8MhSTp8+LB27NihZs2aeTxrVQkyHQAVk5WVpeLiYkVHR5dZHx0drWPHjpXrNQ4fPqzf/va3cjqdcjqdGj16tOLj46sirs9wx7hL0m233aZ///vfys/PV0xMjFJTU9WpUyd3x/VZ5f0+MM7uVZ5x5+dK1SjP2O/cuVOPPfaYAgICZLPZ9MYbb6hOnTom4lYJioqXstlsZT52Op3nrbuYhIQEZWZmVkEq31eZcZek1atXuzuSX7rc94FxrhqXGnd+rlStS4199+7dtX37dhOxPIJDP16mbt26CgwMPO9/8SdOnDivccN9GHdr4PtgBuNuDmNPUfE6wcHBSkhIUFpaWpn1aWlp6t69u6FUvo9xtwa+D2Yw7uYw9hz6saS8vDzt2bPH9fG+ffuUmZmpOnXqqHHjxho3bpyGDx+uxMREdevWTbNmzdLBgwc1cuRIg6m9H+NuDXwfzGDczWHsL8PY+Ua4qPT0dKek85bk5GTXNlOnTnU2adLEGRwc7OzYsaNz3bp15gL7CMbdGvg+mMG4m8PYXxr3+gEAAJbFHBUAAGBZFBUAAGBZFBUAAGBZFBUAAGBZFBUAAGBZFBUAAGBZFBUAAGBZFBUAAGBZFBUAAGBZFBUAAGBZFBUAlnLnnXfKZrNdcFm+fLnpeAA8jHv9ALCU7OxsFRYWKi8vTy1atNDKlSvVoUMHSVLdunUVFMRN3wF/QlEBYEkbN27U9ddfr5ycHIWHh5uOA8AQDv0AsKRt27bp6quvpqQAfo6iAsCStm3bpvj4eNMxABhGUQFgSfv371fLli1NxwBgGEUFgCWVlJTowIEDOnz4sJhKB/gvJtMCsKRVq1bp0Ucf1alTp2S32xUQwP+rAH9EUQEAAJbFf1EAAIBlceUkAB5hs9k8/jXZYQx4P4oKAI+gNAC4Ehz6AeC1tm7dqmrVqqlHjx6mowCoIkymBeC1evTooRtuuEFTp05VTk6OkcNLAKoWe1QAeNyXX36pzp07q3r16qpbt67+/ve/u577/vvvdeONNyo0NFTt27fXV199JZvNpm+//bbMa8ybN0+1a9fWqFGjlJubq71793r6bQDwAIoKAI9auXKl7r77bj3xxBPatm2bHnvsMY0bN04//vijvv/+e3Xt2lU9evTQ1q1b9fLLL+uee+5RtWrV1Lp1a9dr5Ofn64UXXtDEiRMVExOjyMhIZWZmmntTAKoMh34AeMzZs2fVokULTZw4Uffdd58kqbi4WLVr19bkyZP13nvvqX79+lqwYIHrc5KSkrR3794yReSFF17QyZMnNWPGDElS9+7ddcstt+iVV17x6PsBUPU46weAx6xdu1ZnzpzRkCFDXOsCAwMVFBSkn376SZ9//rm+++67Mp8TEhKidu3auT7eu3evZs2aVWa7uLg49qgAPoqiAsBj0tPT1a5dOwUGBrrW7du3T6dOnVLNmjUVHBysNm3alPmcnTt3Kjk52fXx2LFjlZ2drZiYGNe6kpISNWrUqOrfAACPo6gA8JitW7eqoKCgzLo333xTHTt21DXXXKOioiKdPXtW1atXlyStW7dO3377rWuPSlpamr766itt3bpVQUH//fG1efNmPfzww8rOzlZUVJTn3hCAKkdRAeAxmZmZKikp0Zw5c9StWzctXLhQ06dP11dffaUGDRqoWrVqevbZZzV27Fjt2LFDY8aMkSS1b99eRUVFevrpp/Xss8+qffv2ZV43IiLC9fq33HKLh98VgKrEWT8APOLgwYPKzs7W3Llz9dprr6lt27ZatmyZPvnkE3Xs2FENGjTQv/71Ly1btkzx8fH617/+pYceekjNmzdXnTp19Oabbyo7O1ujR48+77VjY2NVo0YN5qkAPoizfgB4xPLly/XQQw8pOzu7XNuXlJSoV69euv766/Xqq69WcToAVsWhHwAesXXrVrVt2/aiz69fv14///yzOnTooKysLE2aNEn79+9XamqqB1MCsBqKCgCP2Lp1q+Lj4y/6/PHjx/X888/ryJEjio6OVu/evfXNN9+oTp06HkwJwGo49AMAACyLybQAAMCyKCoAAMCyKCoAAMCyKCoAAMCyKCoAAMCyKCoAAMCyKCoAAMCyKCoAAMCyKCoAAMCyKCoAAMCyKCoAAMCy/h8GrdJcMcT4JwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "Q = 1000\n",
    "plt.plot(Q**2/Karr**2, Harr, c='k')\n",
    "plt.xscale('log')\n",
    "plt.axvline(x=1/300, c='b')\n",
    "plt.xlabel(r'$\\dfrac{\\tau}{\\rho g A}$')\n",
    "plt.ylabel(r'H')\n",
    "plt.show()                # 描画"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.12"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}