CEE314Problemsets4Q3_teaching.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/AlexCHEN-Engineer/CEE314HW/blob/main/CEE314Problemsets4Q3_teaching.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "mediterranean-century"
      },
      "outputs": [],
      "source": [
        "%time\n",
        "# Do not revise here. \n",
        "try:\n",
        "    import google.colab\n",
        "except ImportError:\n",
        "    import ufl\n",
        "    import dolfin\n",
        "    import mshr\n",
        "else:\n",
        "    try:\n",
        "        import ufl\n",
        "        import dolfin\n",
        "        import mshr\n",
        "    except ImportError:\n",
        "        !wget \"https://fem-on-colab.github.io/releases/fenics-install.sh\" -O \"/tmp/fenics-install.sh\" && bash \"/tmp/fenics-install.sh\"\n",
        "        import ufl\n",
        "        import dolfin\n",
        "        import mshr"
      ],
      "id": "mediterranean-century"
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Variational formulation\n",
        "\n",
        "In this part, we will solve a coupled poromechanics problems.\n",
        "\n",
        "For a domain denoted as $\\Omega$, the boundary is $\\Gamma = \\Gamma_u\\cup \\Gamma_t = \\Gamma_p\\cup\\Gamma_q$, $\\Gamma_u\\cap \\Gamma_t=\\phi$ and $\\Gamma_p\\cup\\Gamma_q=\\phi$. \n",
        "\n",
        "The strong form can be expressed as:\n",
        "\\begin{equation}\n",
        "\\begin{aligned}\n",
        "&\\nabla(\\boldsymbol{\\sigma}^\\prime-p\\boldsymbol{1})+\\rho\\boldsymbol{g}=\\boldsymbol{0},   \\\\\n",
        "&\\nabla⋅\\dot{\\boldsymbol{u}}+\\nabla⋅{\\boldsymbol{q}}=0, \\\\\n",
        "& \\boldsymbol{u} = \\overline{\\boldsymbol{u}}\\quad on \\quad \\Gamma_u, \\\\\n",
        "& \\boldsymbol{\\sigma}^\\prime\\cdot\\boldsymbol{n} = \\overline{\\boldsymbol{t}}\\quad on\\quad \\Gamma_t, \\\\\n",
        "& p = \\overline{p}\\quad on \\quad \\Gamma_p, \\\\\n",
        "& \\boldsymbol{q}\\cdot \\boldsymbol{n} = \\overline{q}\\quad on \\quad \\Gamma_q.\n",
        "\\end{aligned}\n",
        "\\end{equation}\n",
        "\n",
        "As for the constitutive equations, we need to refer to\n",
        "\\begin{equation}\n",
        "\\begin{aligned}\n",
        "& \\boldsymbol{\\sigma}^\\prime = \\mathbb{C}^e:\\boldsymbol{\\epsilon} ,\\\\\n",
        "& \\boldsymbol{q} = -\\frac{\\boldsymbol{\\kappa}}{\\mu_f}\\cdot \\nabla p.\n",
        "\\end{aligned}\n",
        "\\end{equation}\n"
      ],
      "metadata": {
        "id": "P2bXKmcEqOJN"
      },
      "id": "P2bXKmcEqOJN"
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Derive the weak form\n",
        "\n",
        "The trial space is defined as \n",
        "\\begin{equation}\n",
        "\\begin{aligned}\n",
        "& S_u = \\{\\boldsymbol{u}:\\Omega→\\mathbb{R}^{ndim}|\\boldsymbol{u}→\\boldsymbol{H}^1,\\ \\boldsymbol{u}=\\overline{\\boldsymbol{u}}\\quad on \\quad \\Gamma_u \\} ,\\\\\n",
        "&S_p = \\{p:\\Omega→\\mathbb{R}|p→H^1,\\ p=\\overline{p}\\quad on \\quad \\Gamma_p \\}\n",
        "\\end{aligned}\n",
        "\\end{equation}\n",
        "\n",
        "For the weighting space,\n",
        "\\begin{equation}\n",
        "\\begin{aligned}\n",
        "& V_u = \\{\\boldsymbol{\\eta}:\\Omega→\\mathbb{R}^{ndim}|\\boldsymbol{\\eta}→\\boldsymbol{H}^1,\\ \\boldsymbol{\\eta}={\\boldsymbol{0}}\\quad on \\quad \\Gamma_u \\} ,\\\\\n",
        "&V_\\psi = \\{\\psi:\\Omega→\\mathbb{R}|\\psi→H^1,\\ \\psi=0\\quad on \\quad \\Gamma_\\psi \\}\n",
        "\\end{aligned}\n",
        "\\end{equation}\n"
      ],
      "metadata": {
        "id": "MLHKWGvgv_Cg"
      },
      "id": "MLHKWGvgv_Cg"
    },
    {
      "cell_type": "markdown",
      "source": [
        "The weak form becomes:\n",
        "\\begin{equation}\n",
        "\\begin{aligned}\n",
        "& \\mathcal{G}_1 = \\int_{\\Omega}(\\nabla^s\\boldsymbol{\\eta}:\\boldsymbol{\\sigma}^\\prime-p\\nabla\\cdot\\boldsymbol{\\eta}-\\boldsymbol{\\eta}\\cdot{\\rho}\\boldsymbol{g})dV - \\int_{\\Gamma_t}\\boldsymbol{\\eta}\\cdot\\overline{\\boldsymbol{t}}dA = \\boldsymbol{0} \\\\\n",
        "& \\mathcal{G}_2 = \\int_{\\Omega}(\\psi\\nabla\\cdot\\dot{\\boldsymbol{u}}-\\nabla \\psi\\cdot \\boldsymbol{q})dV + \\int_{\\Gamma_q} \\psi\\bar{q}dA\\\\\n",
        "\\end{aligned}\n",
        "\\end{equation}"
      ],
      "metadata": {
        "id": "62rA3t9lxyo6"
      },
      "id": "62rA3t9lxyo6"
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Exact solution\n",
        "\n",
        "The exact solution for this problem is \n",
        "\\begin{equation}\n",
        "u(Z,T) = \\sum_{k=0}^{\\infty}\\frac{4w}{(2k+1)\\pi}\\sin\\left(\\frac{2k+1}{2}\\pi Z \\right)\\exp\\left[-{\\left(\\frac{2k+1}{2}\\pi \\right)^2T}\\right],\n",
        "\\end{equation}\n",
        "- $u$:excess pore pressure;\n",
        "- $Z=\\frac{z}{H}$;\n",
        "- $T = \\frac{c_v t}{H^2}$;\n",
        "- Coefficient of consolidation: $c_v=\\frac{K}{\\gamma_w m_v}$;\n",
        "- Coefficient of compressibility: $m_v=\\frac{1}{E}$ for $\\nu=0$."
      ],
      "metadata": {
        "id": "WyCO_zcJhSGc"
      },
      "id": "WyCO_zcJhSGc"
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yWGBMCZ5jAX2"
      },
      "source": [
        "# Import necessary library"
      ],
      "id": "yWGBMCZ5jAX2"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "infectious-train"
      },
      "outputs": [],
      "source": [
        "from fenics import *\n",
        "from mshr import *\n",
        "from ufl import nabla_div"
      ],
      "id": "infectious-train"
    },
    {
      "cell_type": "code",
      "source": [
        "# Create mesh and define expression\n",
        "# mesh = RectangleMesh.create([Point(0, 0), Point(1,1)],[2,2],CellType.Type.quadrilateral)\n",
        "# We can still create the mesh here.\n",
        "# But, we can't use plot to see the results if we use quadrilateral mesh."
      ],
      "metadata": {
        "id": "xXlk8W3E671m"
      },
      "id": "xXlk8W3E671m",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Parameters\n",
        "\n",
        "For the conversion between elastic modulus, refer to this [wiki page](https://en.wikipedia.org/wiki/Elastic_modulus)."
      ],
      "metadata": {
        "id": "PsmZYVTj5Fd3"
      },
      "id": "PsmZYVTj5Fd3"
    },
    {
      "cell_type": "code",
      "source": [
        "T = 1000                 #final time\n",
        "num_steps =100        #number of time steps\n",
        "dt = T / (1.0*num_steps) #time step size\n",
        "\n",
        "rho_l = 1000\n",
        "\n",
        "\n",
        "# E = 1000 #Pa\n",
        "# nu = 0.0\n",
        "k = 1.e-7 \t#intrinsic permeability\n",
        "mu_l = 1e-2 \t#dynamic viscosity of IF\n",
        "\n",
        "lambda_= 0.\n",
        "mu = 500"
      ],
      "metadata": {
        "id": "4Au2rvfqZMTc"
      },
      "id": "4Au2rvfqZMTc",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Create mesh and define expression\n",
        "mesh=RectangleMesh(Point(0.0,0.0), Point(1., 10.), 1, 20,'crossed')\n",
        "# mesh = RectangleMesh.create([Point(0, 0), Point(1,10)], [1, 20], CellType.Type.quadrilateral)\n",
        "top =  CompiledSubDomain(\"near(x[1], side)\", side = 10.)\n",
        "bottom =  CompiledSubDomain(\"near(x[1], side)\", side = 0.0)\n",
        "left_right = CompiledSubDomain(\"(near(x[0], 0.0) )|| (near(x[0], 1.) )\")\n",
        "\n",
        "boundary = MeshFunction('size_t', mesh, mesh.topology().dim()-1)\n",
        "boundary.set_all(1)\n",
        "top.mark(boundary, 2)\n",
        "ds = Measure('ds', domain = mesh, subdomain_data = boundary)\n",
        "\n",
        "T  = Expression(('0.0', '-10'),degree=1)  # Load on the boundary\n",
        "\n",
        "#Define Mixed Space (R2,R) -> (u,p)\n",
        "V = VectorElement(\"CG\", mesh.ufl_cell(), 2) # define the order\n",
        "W = FiniteElement(\"CG\", mesh.ufl_cell(), 1) # define the order\n",
        "L = FunctionSpace(mesh,W)\n",
        "MS = dolfin.FunctionSpace(mesh, MixedElement([V,W]))\n",
        "\n",
        "# Define boundary condition\n",
        "bcu1 = DirichletBC(MS.sub(0).sub(0), 0.0, left_right)  # slip condition\n",
        "bcu2 = DirichletBC(MS.sub(0).sub(1), 0.0, bottom)      # slip condition\n",
        "bcp = DirichletBC(MS.sub(1), 0.0, top) # drained condition\n",
        "bc=[bcu1,bcu2,bcp]\n",
        "\n",
        "# Define strain\n",
        "def epsilon(u):\n",
        "    return 0.5*(nabla_grad(u) + nabla_grad(u).T)\n",
        "\n",
        "def sigma(u):\n",
        "   return lambda_*tr(epsilon(u))*Identity(2) + 2.0*mu*epsilon(u)\n",
        "\n",
        "# Define variational problem and initial condition\n",
        "X0 = Function(MS)\n",
        "B = TestFunction(MS)\n",
        "\n",
        "e_u0 = Expression(('0.0', '0.0'), degree=1)\n",
        "e_p0 = Expression('10.0', degree=1)\n",
        "\n",
        "u0 = interpolate(e_u0, MS.sub(0).collapse())\n",
        "p0 = interpolate(e_p0, MS.sub(1).collapse())\n",
        "\n",
        "Xn = Function(MS)\n",
        "assign(Xn, [u0, p0])\n",
        "\n",
        "(u,p)=split(X0)\n",
        "(u_n,p_n)=split(Xn)\n",
        "(v,q)=split(B)\n",
        "\n",
        "F = (1/dt)*nabla_div(u-u_n)*q*dx + (k/(mu_l))*dot(grad(p),grad(q))*dx \n",
        "F += inner(sigma(u), epsilon(v))*dx -p*nabla_div(v)*dx- dot(T,v)*ds(2)\n",
        "\n",
        "#solver tuning\n",
        "dX0 = TrialFunction(MS)\n",
        "J = derivative(F, X0, dX0)\n",
        "Problem = NonlinearVariationalProblem(F, X0, J = J, bcs = bc)\n",
        "Solver  = NonlinearVariationalSolver(Problem)\n",
        "Solver.parameters['newton_solver']['convergence_criterion'] = 'incremental'\n",
        "Solver.parameters['newton_solver']['relative_tolerance'] = 1.e-11\n",
        "Solver.parameters['newton_solver']['absolute_tolerance'] = 5.e-10\n",
        "\n",
        "#if you want to store the results\n",
        "vtkfile_u = File('Terzaghi_CG/u.pvd')\n",
        "vtkfile_p = File('Terzaghi_CG/p.pvd')\n",
        "\n",
        "t = 0\n",
        "\n",
        "p_out = []\n",
        "\n",
        "for n in range(num_steps):\n",
        "\tt += dt\n",
        "\tSolver.solve()    \n",
        "\tassign(Xn,X0)\n",
        "\n",
        "\t(_u,_p)=X0.split()\n",
        "\t\n",
        "\tvtkfile_u << (_u,t)\n",
        "\tvtkfile_p << (_p,t)\n",
        "\tp_out.append(_p(0.5, 5))"
      ],
      "metadata": {
        "id": "pw-EBzjqQiha"
      },
      "id": "pw-EBzjqQiha",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# u_result, p_result = X0.split()\n",
        "# import matplotlib.pyplot as plt\n",
        "# c=plot(p_result)\n",
        "# plt.colorbar(c)\n",
        "# plt.show()"
      ],
      "metadata": {
        "id": "G0dWBTs8TB3q"
      },
      "id": "G0dWBTs8TB3q",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Postprocess"
      ],
      "metadata": {
        "id": "bZFw87laa7Tp"
      },
      "id": "bZFw87laa7Tp"
    },
    {
      "cell_type": "code",
      "source": [
        "# Use the following method to zip the file.\n",
        "!zip -r 'Terzaghi_CG.zip' '/content/Terzaghi_CG' "
      ],
      "metadata": {
        "id": "joHFy9lSTN27"
      },
      "id": "joHFy9lSTN27",
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 3(a)"
      ],
      "metadata": {
        "id": "5_CjpEE2a4n1"
      },
      "id": "5_CjpEE2a4n1"
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "depth = np.linspace(0,10,1000)\n",
        "pressure = np.array([_p(0.5, 10-d) for d in depth])\n",
        "plt.xlabel(r\"$p/w$\")\n",
        "plt.ylabel(r\"z/h\")\n",
        "plt.gca().invert_yaxis()\n",
        "plt.plot(pressure/10, depth/10, 'k')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 300
        },
        "id": "kXhNsYMdVkIu",
        "outputId": "18aabb8a-7307-43bd-b4c4-dc3b46fae9b9"
      },
      "id": "kXhNsYMdVkIu",
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[<matplotlib.lines.Line2D at 0x7f0d8669f690>]"
            ]
          },
          "metadata": {},
          "execution_count": 49
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEKCAYAAAASByJ7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAVbUlEQVR4nO3df5BddZ3m8feT7k4ggIomzmIghnWhJJPFgulldASN4s4E1FDLsFMwo6yzrllGmdpaxF22xkIW/1CkdKZkGZVBixmrFEFYK+sE44o4Kr8bRAyhgIgjCVqbCCQICEknn/3jnmAnpn8k5tx70/1+VXXVOed+++bJ7b73ued8b5+TqkKSpFm9DiBJ6g8WgiQJsBAkSQ0LQZIEWAiSpMZgrwPsrXnz5tWiRYt6HUOSDij33HPPL6pq/kRjDrhCWLRoESMjI72OIUkHlCQ/nWyMh4wkSYCFIElqWAiSJMBCkCQ1LARJEtByISRZluShJOuSXLSH2+ck+Upz+51JFrWZR5I0vtYKIckAcCVwGrAYOCfJ4t2GvRd4qqr+FfDXwGVt5ZEkTazNv0M4CVhXVY8CJLkWOANYO2bMGcAlzfJXgf+VJNXCObm///3v881vfpOhoSEGBwf3+DXebftre5L9/d+SpP2mzUJYAKwfs74B+P3xxlTVaJItwCuAX4wdlGQFsAJg4cKF+xTm9ttv56Mf/eg+fe/+MjAwsF+KpRvltb8yzZo1yyKUDhAHxF8qV9VVwFUAw8PD+7T38KEPfYgLL7yQHTt2MDo6usvXtm3bfmNbN7bv7fc8//zz+/Rv79ixY7/+PPZWv5WUe4PSnrVZCI8DR41ZP7LZtqcxG5IMAi8FnmgrUBIGBgYYGBhgzpw5bf0zfWfHjh1s3769lZJqe/vzzz+/z/fVS3u7N3gg7fW5Nzh9tVkIdwPHJDmazgv/2cCf7jZmJfAfgNuBs4BvtzF/MNPNmjWLWbNmMTQ01OsoXVNVPdkb3J9F6N7ggVGc06kEWyuEZk7gfGA1MAB8oaoeSHIpMFJVK4HPA19Msg54kk5pSL819wZ7v3c3U/YGd77ZaruMDj30UC644AJe+cpXtvZ/aXUOoapWAat223bxmOXngX/fZgZpJnFvcP+W1HPPPceWLVvYvHkzmzdvZsuWLS9+7VzfvHkzTz/9NC+88EJr/8ckvPzlL+ess846cAtBkto2OjrKs88+u8vXM888s1+2jY6O7lWWuXPncuihh3LIIYfs8vXbbjvooIO6cmjKQpDUuu3bt+/Vi/LejN26deteZTn44IP3+CK8YMGC3+oF/OCDD2bWrAP7bEAWgiSgM//Q1jvtvT2cMmfOnD2+AB9xxBG/sW1vXsDnzp3LwMBAS4/ggc9CkKa5quLDH/4wP/7xjyd8Af/Vr361V/c7e/bsPb4Az58/n0WLFu3zO+25c+cyOOhLUy/4qEvT3NatW7nhhht46KGHXtx2yimnsHDhwn1+p33IIYfMqInrmcJCkKa5OXPm8KMf/YgbbriBT37yk4yMjLBmzRrOO+88zj//fF71qlf1OqL6xIE9AyJpSoaGhjj77LO56667+N73vsfSpUv5+Mc/zqJFizj33HO57777eh1RfcBCkGaQJJx88snceOONPPLII5x33nnceOONnHDCCbz1rW/l61//es//2lm9YyFIM9RrXvMaPv3pT7N+/Xo+8YlP8Mgjj/DOd76T4447jrVr105+B5p2LARphjv88MO58MIL+djHPsbg4CA/+clP2LhxY69jqQcsBGmG+/nPf86ZZ57Ju9/9bk444QTuvfdeli5d2utY6gELQZqhqoprrrmGxYsX841vfIPLL7+c2267jSVLlvQ6mnrEj51KM9Bjjz3GihUrWL16NSeffDKf//znOfbYY3sdSz1mIUgzzC233MLy5cupKq644gre//73H/Dn4NH+YSFIM8yGDRt45pln+NznPseKFSt6HUd9JAfaBcqGh4drZGSk1zGkA9bo6CjHH38827dvZ82aNZ6CYoZIck9VDU80xv1EaYYZHBzksssu4+GHH+bqq6/udRz1EQtBmoHe8Y538KY3vYlLLrmEX/7yl72Ooz5hIUgzUBIuv/xyNm7cyBVXXNHrOOoTFoI0Qy1ZsoTZs2fz1FNP9TqK+oSFIM1Qt956K1u3buXUU0/tdRT1CQtBmqG+9a1vMTQ0xCmnnNLrKOoTFoI0Q91888284Q1v4JBDDul1FPUJC0GagZ588knuvfde3va2t/U6ivqIhSDNQOvXr6eqOO6443odRX3EQpBmoMMOOwyAZ599tsdJ1E8sBGkG2lkI/lGaxrIQpBnoJS95CQBPP/10j5Oon1gI0gw0Z84choaG3EPQLiwEaYY67LDDLATtwkKQZqhZs2axY8eOXsdQH2m1EJIsS/JQknVJLtrD7RckWZvk/iQ3J3l1m3kk/dq2bduYPXt2r2Ooj7RWCEkGgCuB04DFwDlJFu827AfAcFUdD3wV+ERbeSTtauvWrV4cR7tocw/hJGBdVT1aVVuBa4Ezxg6oqluq6rlm9Q7gyBbzSBpj27ZtFoJ20WYhLADWj1nf0Gwbz3uBm/Z0Q5IVSUaSjGzatGk/RpRmpqpidHTUQtAu+mJSOcm7gGHg8j3dXlVXVdVwVQ3Pnz+/u+GkaWh0dBTAQtAuBlu878eBo8asH9ls20WStwF/Bby5ql5oMY+kxs5PFw0MDPQ4ifpJm3sIdwPHJDk6yWzgbGDl2AFJTgA+Byyvqo0tZpEkTaK1QqiqUeB8YDXwIHBdVT2Q5NIky5thlwOHAtcnuS/JynHuTpLUsjYPGVFVq4BVu227eMyyJ2OXpD7RF5PKkqTesxAkSYCFIElqWAiSJMBCkCQ1LARJEmAhSJIaFoI0g1VVryOoj1gI0gyUpNcR1IcsBEkSYCFIkhoWgiQJsBAkSQ0LQZIEWAiSpIaFIEkCLARJUsNCkCQBFoIkqWEhSJIAC0GS1LAQJEmAhSBJalgIkiTAQpAkNSwESRJgIUiSGhaCJAmwECRJDQtBkgS0XAhJliV5KMm6JBdNMO6Pk1SS4TbzSJLG11ohJBkArgROAxYD5yRZvIdxhwH/BbizrSySpMm1uYdwErCuqh6tqq3AtcAZexj3UeAy4PkWs0iSJtFmISwA1o9Z39Bse1GSE4GjquofJ7qjJCuSjCQZ2bRp0/5PKknq3aRyklnAp4APTja2qq6qquGqGp4/f3774SRpBmqzEB4HjhqzfmSzbafDgCXAd5L8M/B6YKUTy5LUG20Wwt3AMUmOTjIbOBtYufPGqtpSVfOqalFVLQLuAJZX1UiLmSRJ42itEKpqFDgfWA08CFxXVQ8kuTTJ8rb+XUnSvhls886rahWwardtF48zdmmbWSRJE/MvlSVJgIUgSWpYCJIkwEKQJDUsBEkSYCFIkhoWgiQJsBAkSQ0LQZIEWAiSpMaUCiHJmUkeSbIlydNJfpnk6bbDSZK6Z6rnMvoE8M6qerDNMJKk3pnqIaP/ZxlI0vQ24R5CkjObxZEkXwG+Bryw8/aqurHFbJKkLprskNE7xyw/B/zhmPUCLARJmiYmK4T/C6yuqie6EUaS1DuTFcJRwPVJhoCbgZuAu6qqWk8mSeqqCSeVq+qyqnorcDrwQ+A/Avcm+VKSc5P8TjdCSpLaN9WPnX4N+FRV/WeAJL8L3Ad8G/ijlrJJkrpoqh87PRr4b0k+AlBVDwBrqsoykKRpYqqFsBk4FfidJP8nyUvpfMpIkjRNTLUQUlWjVfV+4Abg+8Ar24slSeq2qc4hfHbnQlVdk+RHwAfaiSRJ6oUpFUJVfW639XvofOJIkjRNePprSRJgIUiSGhaCJAmwECRJDQtBkgRYCJKkhoUgSQJaLoQky5I8lGRdkovGGfMnSdYmeSDJl9rMI0ka31T/UnmvJRkArgT+LbABuDvJyqpaO2bMMcD/AN5YVU8l8XQYktQjbe4hnASsq6pHq2orcC1wxm5j3gdcWVVPAVTVxhbzSJIm0GYhLADWj1nf0Gwb61jg2CS3JrkjybI93VGSFUlGkoxs2rSppbiSNLP1elJ5EDgGWAqcA/xdkpftPqiqrqqq4aoanj9/fpcjStLM0GYhPE7nmsw7HdlsG2sDsLKqtlXVT4CH6RSEJKnL2iyEu4FjkhydZDZwNrBytzFfo7N3QJJ5dA4hPdpiJknSOForhKoaBc4HVgMPAtdV1QNJLk2yvBm2GngiyVrgFuBDVfVEW5kkSeNr7WOnAFW1Cli127aLxywXcEHzJUnqoV5PKkuS+oSFIEkCLARJUsNCkCQBFoIkqWEhSJIAC0GS1LAQJEmAhSBJalgIkiTAQpAkNSwESRJgIUiSGhaCJAmwECRJDQtBkgRYCJKkhoUgSQIsBElSw0KQJAEWgiSpYSFIkgALQZLUsBAkSYCFIElqWAiSJMBCkCQ1LARJEmAhSJIaFoIkCWi5EJIsS/JQknVJLtrD7QuT3JLkB0nuT3J6m3kkSeNrrRCSDABXAqcBi4FzkizebdiHgeuq6gTgbOBv28ojSZpYm3sIJwHrqurRqtoKXAucsduYAl7SLL8U+FmLeSRJE2izEBYA68esb2i2jXUJ8K4kG4BVwF/u6Y6SrEgykmRk06ZNbWSVpBmv15PK5wDXVNWRwOnAF5P8RqaquqqqhqtqeP78+V0PKUkzQZuF8Dhw1Jj1I5ttY70XuA6gqm4HDgLmtZhJkjSONgvhbuCYJEcnmU1n0njlbmMeA04FSHIcnULwmJAk9UBrhVBVo8D5wGrgQTqfJnogyaVJljfDPgi8L8kPgS8D76mqaiuTJGl8g23eeVWtojNZPHbbxWOW1wJvbDODJGlqej2pLEnqExaCJAmwECRJDQtBkgRYCJKkhoUgSQIsBElSw0KQJAEWgiSpYSFIkgALQZLUsBAkSYCFIElqWAiSJMBCkCQ1LARJEmAhSJIaFoIkCbAQJEkNC0GSBFgIkqSGhSBJAiwESVLDQpAkARaCJKlhIUiSAAtBktSwECRJgIUgSWpYCJIkwEKQJDVaK4QkX0iyMcmacW5Pkk8nWZfk/iQntpVFkjS5NvcQrgGWTXD7acAxzdcK4DMtZpEkTaK1Qqiq7wJPTjDkDOAfquMO4GVJjmgrjyRpYr2cQ1gArB+zvqHZ9huSrEgykmRk06ZNXQknTWdJOOuss3jta1/b6yjqI4O9DjAVVXUVcBXA8PBw9TiOdMAbGhri+uuv73UM9Zle7iE8Dhw1Zv3IZpskqQd6WQgrgXObTxu9HthSVT/vYR5JmtFaO2SU5MvAUmBekg3AR4AhgKr6LLAKOB1YBzwH/HlbWSRJk2utEKrqnEluL+ADbf37kqS9418qS5IAC0GS1LAQJEmAhSBJaqQzt3vgSLIJ+Ok+fvs84Bf7Mc7+1s/5zLZv+jkb9Hc+s+27PeV7dVXNn+ibDrhC+G0kGamq4V7nGE8/5zPbvunnbNDf+cy27/Y1n4eMJEmAhSBJasy0Qriq1wEm0c/5zLZv+jkb9Hc+s+27fco3o+YQJEnjm2l7CJKkcVgIkiRgmhZCkmVJHkqyLslFe7h9TpKvNLffmWRRH2W7IMnaJPcnuTnJq7uVbSr5xoz74ySVpGsfvZtKtiR/0jx+DyT5Ur9kS7IwyS1JftD8bE/vYrYvJNmYZM04tyfJp5vs9yc5sY+y/VmT6UdJbkvyun7JNmbcv0kymuSsbmVr/t1J8yVZmuS+5vnwT5PeaVVNqy9gAPgx8C+B2cAPgcW7jXk/8Nlm+WzgK32U7S3A3Gb5L7qVbar5mnGHAd8F7gCG+yUbcAzwA+DwZv2VfZTtKuAvmuXFwD938ef6JuBEYM04t58O3AQEeD1wZx9l+4MxP8/T+inbmJ/9t+mczv+sbmWb4mP3MmAtsLBZn/T5MB33EE4C1lXVo1W1FbgWOGO3MWcAf98sfxU4NUn6IVtV3VJVzzWrd9C5kly3TOWxA/gocBnwfJ9lex9wZVU9BVBVG/soWwEvaZZfCvysS9moqu8CT04w5AzgH6rjDuBlSY7oh2xVddvOnyddfj5M4XED+EvgBqBbv2svmkK+PwVurKrHmvGTZpyOhbAAWD9mfUOzbY9jqmoU2AK8ok+yjfVeOu/cumXSfM3hhKOq6h+7mAum9tgdCxyb5NYkdyRZ1kfZLgHe1VwsahWdF5J+sbe/l73S7efDhJIsAP4d8JleZxnHscDhSb6T5J4k5072Da1dIEe/nSTvAoaBN/c6y05JZgGfAt7T4yjjGaRz2GgpnXeS303yr6tqc09TdZwDXFNVn0zyBuCLSZZU1Y5eBzsQJHkLnUI4uddZxvgb4L9X1Y7uHGDYa4PA7wGnAgcDtye5o6oenugbppvHgaPGrB/ZbNvTmA1JBunswj/RJ9lI8jbgr4A3V9ULXci102T5DgOWAN9pngD/AliZZHlVjfQ4G3Te2d5ZVduAnyR5mE5B3N0H2d4LLAOoqtuTHETnBGRdP9SwB1P6veyVJMcDVwOnVVU3nqdTNQxc2zwX5gGnJxmtqq/1NtaLNgBPVNWzwLNJvgu8Dhi3ELo2AdLFiZZB4FHgaH49wfe7u435ALtOKl/XR9lOoDNBeUw/Pna7jf8O3ZtUnspjtwz4+2Z5Hp3DIK/ok2w3Ae9plo+jM4eQLv5sFzH+5OPb2XVS+a4u/95NlG0hneuu/0E3M00l227jrqHLk8pTeOyOA25ufj/nAmuAJRPd37TbQ6iq0STnA6vpfALgC1X1QJJLgZGqWgl8ns4u+zo6kzJn91G2y4FDgeubdx6PVdXyPsrXE1PMthr4wyRrge3Ah6oL7yinmO2DwN8l+a90JpjfU82ztm1JvkznMNq8Zg7jI8BQk/2zdOY0Tqfzwvsc8OfdyDXFbBfTmd/72+b5MFpdOsvoFLL11GT5qurBJN8A7gd2AFdX1cQfoe3S76Qkqc9Nx08ZSZL2gYUgSQIsBElSw0KQJAEWgiSpYSFIkgALQZLUsBCk/STJZ5O8sdc5pH3lH6ZJ+0mS+4Dfq6rtvc4i7Qv3EKS9kOTLzdX27kry0yRvb7YfR+ekYUc1JxEjyYnNVeXmJRlorvo1t5f5pYlYCNLeeR3waFWdBPwZnfPHQOdqXt8ANtM5FxV0rnlwB50rV70d+Fb9+uJHUt+xEKQpak5ZPR/4n82mtcDhzfIf0SmEp4G5SeYBRwC3NmNW0L8XUpEAC0HaG0uAR6pq56VDTwR+2BwGellV/aw6F7wp4D/ROavuL+nsVQzUBBcmkfqBhSBN3euAhUkOSnIInT2FvwbeAtwyZtwOYDnwv+nsMXwQ6PnpkqXJWAjS1L0OuBG4k85V2D5TVbfy6/mDnbYBN1Xnet1P07k4yde7nFXaa37sVJqiJP8ErKiqh3bbfi/w+9W5dKd0wLIQpClqrkq1sJknkKYdC0GSBDiHIElqWAiSJMBCkCQ1LARJEmAhSJIaFoIkCbAQJEmN/w8Y2YVB2cZI9AAAAABJRU5ErkJggg==\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "depth = np.linspace(0,10,1000)\n",
        "pressure = np.array([_p(0.5, 10-d) for d in depth])\n",
        "plt.xlabel(r\"$p/w$\")\n",
        "plt.ylabel(r\"z/h\")\n",
        "plt.gca().invert_yaxis()\n",
        "plt.plot(pressure/10, depth/10, 'k')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 300
        },
        "id": "9yAN7gucikr2",
        "outputId": "b01a6f6e-10b0-4c4d-9e50-cc1f5c009666"
      },
      "id": "9yAN7gucikr2",
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[<matplotlib.lines.Line2D at 0x7f0d866fdfd0>]"
            ]
          },
          "metadata": {},
          "execution_count": 51
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de3RU9b338fd3kpCE3LgE5JYAchW1aAtCa6tgbQ0qUHvUQusF0XK0evr0eOnT1i61dj2tPaI9Pa3W2spptS1IrbW0gLQ9oBVFBQRUEDQQ1HAnCIGEkNvv+SOZOUlIwiSZnT175vNaa9bsy2+S72bC77Pv25xziIhI8gr5XYCIiPhLQSAikuQUBCIiSU5BICKS5BQEIiJJLtXvAjoqPz/fDRs2zO8yREQCZf369Qedc/1amxe4IBg2bBjr1q3zuwwRkUAxs/fbmqddQyIiSU5BICKS5BQEIiJJTkEgIpLkFAQiIknO0yAwsyIz22ZmxWb2rVbmp5vZ043zXzOzYV7WIyIiJ/MsCMwsBXgEmAaMA2ab2bgWzW4EPnLOjQR+DPzIq3pERKR1Xl5HcB5Q7JzbAWBmi4CZwJYmbWYC9zUOPwP8zMzMeXBv7NWrV/O3v/2NtLQ0UlNTW321NS9W080s1oslItJlXgbBYODDJuOlwKS22jjnas3sCNAXONi0kZnNA+YBFBYWdqqYNWvW8P3vf79Tn42VlJSUmARKd4RWrGoKhUIKQJE4F4gri51zjwOPA0yYMKFTWwt33XUXd955J/X19dTW1jZ71dTUnDStO6Z39DNVVVWd+t319fUx/T46Kt7CSVt/Is15GQS7gIIm40Map7XWptTMUoE8oMyrgsyMlJQUUlJSSE9P9+rXxJ36+nrq6uq6NbRiNb0z4Ree56eObv0FaStPW3+Jx8sgWAuMMrPhNHT4s4Avt2izBLgeWANcCaz04vhAsguFQoRCIdLS0vwupds45+Jq66+jAautv+AEZiKEn2dB0LjP/zZgBZACLHDObTaz+4F1zrklwBPAU2ZWDByiISxEukxbf/7vyuyOrb/wdD+FV7K8DqHs7Gxuv/12+vfvH/Nl8PQYgXNuGbCsxbR7mgxXAVd5WYNIMomHrb/6+vpIB11TU9PmcFemtRw+ceIEx48fp6qqqtl7ZWVlZLjp9PArFltQ9fX1nDhxghMnTsTgX+9kTYPgyiuvDF4QiMjJnHPU1dXFpJP0omPt6me6e+9uKBSKdJbhNeqm7+HhrKwsevXq1eb8tqadan6sPtPatJSUlG7Z9aQgkLjU3lqlH51bLH+OH7syWnY20XZEGRkZcdtxhodDId0pp6sUBAHV1lplvHWSne1Yu/ugZ9O1yo50ROnp6WRnZ3f7mmJHPtNda5USXAoCj5WUlPDAAw9QWVkZ047Vr7XKznRE6enpvqwpdmSa1iolmSkIPLZnzx7+9re/UVFRQU1NDdXV1VRXV/t+pgNARkYGPXv2JCsri549e0ZeTcfDw5mZmZEOvUePHvTo0SOmw1prFfGPgsBjn/rUpygpKTlpunMuEgxNA6Irw7H4WQcPHmT37t2ttqmp8e4iLTPzJGA6MtzRzyi8JFEoCHzStOMLCucctbW1cRFax44di7q9V8ysW4MnFkGn8JLWKAgkauGOLy0tjaysLL/LiUo4vPze0gqHV7S/z0tB2NpqOpwoV+/GMwWBJLSm4dWzZ0+/y4lK+Iwwv7e6ampqqKys5PDhw1G19/L6gc6GyNChQ3nooYdITVVX1x7964jEGTOLnKHVUeFjTydOnIh00G0Nd3ZeW+2OHz/O0aNHOXbsGBUVFRw9epSqqqqY/JuEz5irqKhodX64409PT4+EQXp6OhUVFb7ffykIFAQiUWq6pt7Vjtardl4c0E9JSWnWubY1nJmZSV5e3inbNd1dFIt2aWlp2nXURQoCiRv19fVx0Zm2N8+L3R/RdIbp6enk5ubGtAONpl34ALMkNgVBkmh6ump37SroaLu6urqYL3fLA5ZtdXg9e/akV69e3d7R6kCoxAMFQUC98cYb3H333VRWVkbV6VZXx/5MlFAoFHWHl52dHfNdAqdqpyuGRaKjIAio6upqDh8+TEVFBSdOnKCqqiry8uqWuFlZWeTm5pKXl0dubi49e/YkIyODjIwM0tPTI8OxHNcas4j3FAQBNXnyZNasWdPm/PD+9nAwNA2JlqHRlfFDhw61OT8WYRTe6vAiZKJtq4ORkugUBAkqFApFOjO/OOciYeRFCIXHP/roo3bbdpWZeRYy0Y4rjMRLCgLxjJmRnp5Oeno6eXl5vtQQDqNYhlBr8w4fPtzm/FiFkVchE23bHj16KIwSlIJAElrTMMrNzfWlhvAZW7EKnbbGjxw5wv79+9ucH4tTX9sKiqysLObPn8+nPvWpGPyLSXdTEIh4rOkNBnNzc5sFQ8tn6bY3rb154Z9//PhxUlNTSUlJiTy/GBrCqKtbJuFAzczMJCMjo9l7VlYWaWn+PSdZukZBIEmttrY26s42Fh12eLgrtz1IS0trtTMOv+fm5p40rbV27f2MltMyMjJ0Km4CUxBIXKirq+tQhxqrDrsrDwhKSUlpt5Pt27dv1B1wRzplXekrsaYgkGbq6+s5ceJEt64dHz9+vEv3yDGzSIfZWkfaq1cvBgwY0KVOubVpuqOlJAr9JSepN954g7lz51JeXt6sQ+7quf/tdZ7Z2dn069evS2vDrU3TqZUiXaMgSFK5ubmMGTOGw4cPU1lZSUVFReQ9PNyZUAjfQjn80PrwgcTws49be492XmZmpnaLiHhAQZCkRo4cydNPP91um/CB1NZCorX39ubt27fvpGnHjx/vcN0ZGRmdDpJo5ml3jyQj/dVLm1JTU8nJySEnJ8eTn19fX39S0HQkXJoOf/TRR5SWlp40r6Pnzocfw9nZcDlVG12UJfFIQSC+CYVCkU7TC+Fz5zsbLk3fjx49yr59+05q09FbZ6ekpMR0C6ble0ZGhoJGOkxBIAmr6dlEffv2jfnPD18YFotdZxUVFZSVlfHRRx9Frg7uDDOLhMJll13GggULYrzUkogUBCKdUFNTw7FjxyLP5w0Ptxzv6LyOXNcQfs5DVlYW2dnZkVd4fNKkSR7+C0giURBIQqutrW2z8+1sB37s2LEOXfeQmpraakc9YMCANjvx9obD47qlg8SKp0FgZkXAT4AU4FfOuQdazL8duAmoBQ4Ac51z73tZk8Snuro6KisrY7ZmHR7uyCmwKSkprXbA/fr1Y/jw4VF1zq0N9+jRw8N/OZGu8ywIzCwFeAT4HFAKrDWzJc65LU2abQAmOOcqzewW4D+AL3lVk3RdfX19ZL92tGvP0bTtyKmkZtZqh9unTx8KCgqiXqNuOZyenq4DrZKUvNwiOA8ods7tADCzRcBMIBIEzrlVTdq/ClzjYT3SBf/5n//Jd7/7XSoqKjr1+T59+jBw4MDIYy4HDRrUqbXr7OxsnRkjEmNeBsFg4MMm46VAe0evbgSWtzbDzOYB8wAKCwtjVZ90wOTJk7nmmmsoLy+nvLycI0eOnPTe3h01Dx06xKFDhyLPPW767OPwe8tp1dXV1NXV4ZwjFApFrlbWmrtIbMXFwWIzuwaYAFzY2nzn3OPA4wATJkzo+tM1pMMmT57M5MmT25zvnKOysrLNkGgrQMrLy9m9e3dk2tGjR09Zi5mRk5MTVZi0Ny8zM1OBIoK3QbALKGgyPqRxWjNmdjFwN3Chc67rTzsXX5hZ5KKngQMHdvrn1NfXc/To0XaDo7WgKSsro6SkJDKtsrLylL8rJSWlQ8HR1rT09PROL69IPPAyCNYCo8xsOA0BMAv4ctMGZnYu8AugyDm338NaJCBCoRB5eXldfsZxTU1NJFA6soWyZ88etm7dGplWXV19yt/Vo0ePLodJbm6u7nMkvvHsL885V2tmtwEraDh9dIFzbrOZ3Q+sc84tAR4EsoE/NG6if+Ccm+FVTZI80tLS6NOnD3369OnSzzlx4sQpd22Fhw8cOMAHH3zAli1b2L+/4+s1P/rRj/jmN7/ZpXpFOsPTVRDn3DJgWYtp9zQZvtjL3y/SESdOnODw4cOnfB05cqTV6ac6BTY1NZXevXvTq1evk155eXlccskl3bSkIs1pW1QSRlVVVZuddDSvU93fJy0t7aQOfMiQIa127K29dHBa4pWCQOJGVVVVpzvxw4cPn/Iq4rS0tJPWyAsLC8nLy1NHLklNQSAxU11dTVlZWVS7UWLVkQ8dOjTqNXJdiCbSOgWBxMywYcPYs2dPVG1TU1MZNmwYI0aMoLCwsM195+rIRbynIJCYeeihh3jzzTc5ePAgBw8epKysLPJeVlbW7CEutbW1FBcXU1xcTGpqKvn5+fTt25f8/PyohnNzcxUKIjFiHX2Un98mTJjg1q1b53cZ0kH19fWRC79aBkVbwy3Do6nU1NQOBYfCQ5Kdma13zk1obZ62CKRbhEIhevfuTe/evRk5cmRUn6mvr6e8vLzVgGg5bevWrZHhU4VHRwIkLy9P4SEJT0EgcSsUCkWOD0QbHs65qLY8ysrK2LZtW2S4rSeDpaam0qdPn6iDY+jQoXpgjASOgkASiplFwmPEiBFRfcY5F9nyaBoau3bt4p133om8tmzZcsqfddVVV7F48eKuLoZIt1IQSNKora2lrKyM/fv3R14HDhxoNt50Wnl5eas/Jz09nf79+9OvXz/69+/f7PX5z3++m5dKpOsUBBJYzjk++uijdjvzpq9Dhw7R2skRKSkp5OfnRzrziRMnttrJh6fl5OTouIEkFAWBxA3nHBUVFVGtrYeH29q336dPn0jnPW7cOKZMmXJShx5+9e7dm1Ao1M1LKxI/FATiuaqqqsgdOU/Vubd147bs7OxIxz106NB219rz8/N1wFakAxQE4rk5c+bw9NNPt9umf//+XHDBBZxzzjmMHTv2pM49MzOzm6oVST66oEw8t337dlasWMHu3bubvXbt2sWhQ4dOat+jRw8GDRoUeQ0ePLjZeHhaTk6OD0sjEkztXVCmIBBfVVVVsWfPnlZDoul4a88yzs7ObjUgmo4PHDhQWxMi6MpiiWMZGRkMHz6c4cOHt9vu6NGj7Nmz56SACL9effVVdu3a1eodTPv06XNSYLQMjdNOO03HFSRpKQgkEHJycsjJyWH06NFttgmfTtre1sWWLVvYs2fPSbehMDP69+/fLCD69u1LSkoKX/3qVyksLPR6EUV8o11DknTq6urYu3cv69ev57XXXuPVV1/ltddeo6KiotX2Y8aMYevWrd1cpUhsadeQJJX6+nr279/Phx9+2OZr9+7d1NfXN/tcdnY2BQUFkdeAAQOoqanhxhtv9GlJRLqHgkACxTlHWVlZq517aWlp5L2mpqbZ5zIyMhgyZAgFBQVcdNFFzTr88Et3GpVkpSCQuHLkyJF21+RLS0tPuugsNTU10sl/8pOfbLWT79u3rzp5kTYoCMQXf/7zn/nTn/5EKBRi9+7dkY6+5WmioVCIgQMHUlBQwDnnnMP06dNP6uRPO+003SJCpAsUBNLtDhw4wBe+8IVW5w0aNIipU6cyZcoULrzwQoYNG6bTOkU8prOGxBcvvfQSGzduZM+ePRQXF7N9+3aKi4ub3frZzBgyZAgjRoxg5MiRjBw5MjI8YsQIXVks0gG6slgCIXwguGkwNH3fv39/s/b9+/dvFgxN33VMQKQ5BYEkhPLycnbs2NEsIMLDpaWlzZ41kJeX12ZIDBw4UMcUJOkoCCThVVVVUVJS0urWxM6dO5s9tyAjI4MRI0aQn59PXV0d8+fPZ9KkST5WL+I9BYEkrbq6OoqLi/n73//O0qVLef75509qc/bZZ/Pmm2/6UJ1I99GVxZLwDh8+zLZt29i2bRtbt26NDL/33ntUV1dH2vXp04exY8cyZswYxo4dy4gRI/j0pz/tY+Ui/vM0CMysCPgJkAL8yjn3QBvt/gV4BpjonNPqvrSqrq6OnTt3Nuvow8P79u2LtEtNTWXEiBGMGTOGSy+9NNLxjxkzhvz8fB+XQCQ+eRYEZpYCPAJ8DigF1prZEufclhbtcoD/A7zmVS0SLK2t3W/dupXi4uJma/d9+/Zl7NixXHbZZZE1/DFjxnD66afr2gORDvByi+A8oNg5twPAzBYBM4EtLdp9H/gRcJeHtUgc27BhA9dddx319fWUlZW1unY/duxYLr/88mYdft++fX2sWiRxeBkEg4EPm4yXAs1OzTCzjwMFzrmlZtZmEJjZPGAeoPvCJ6B77rmHt99+u9m0j33sY8yaNYsvfvGLjBo1Sqd7injIs7OGzOxKoMg5d1Pj+LXAJOfcbY3jIWAlMMc5t9PMXgDuPNUxAp01lHgqKipYs2YNmzZtYuPGjWzYsIGtW7dGHh6Tm5vL+PHjOeecczjnnHM499xzGTduHOnp6T5XLhIcfp01tAsoaDI+pHFaWA5wFvBC4xWgA4AlZjZDB4yTS1ZWFhdffDEXX3xxZNrx48fZvHlzJBg2btzIggULIg+PSU1NZdy4cZx77rmRgBg/fjy9e/f2azFEAsvLLYJU4F3gszQEwFrgy865zW20fwFtEUg76uvrKS4uZuPGjZHXhg0b2Lt3b6TNsGHDIsEQ3nooKCjQ7SYk6fmyReCcqzWz24AVNJw+usA5t9nM7gfWOeeWePW7JTGFQiFGjx7N6NGjufrqqyPT9+7dy6ZNmyJbDhs3buTPf/5z5JYToVCI+vp6HnzwQe68806/yheJW7qyWBJKVVUVmzZtYtWqVTzxxBMUFxdH5g0ZMoQPPvhAWweSlHRlsSSkuro63nnnHV5//XXWrl3L66+/zptvvhm5r9CAAQOYMWMGEydOZOLEiUyaNEkhINIKBYEEgnOOnTt3Rjr9tWvXsn79+sjB49zcXCZOnMidd97JxIkTOe+88xg8eLA6fpEoKAgkLu3fvz+ylh/u+A8ePAhAeno65557LnPnzo10+rrWQKTzFAQSN8rLy7nssstYvXp1ZFooFGLcuHGRXTznnXceZ511Fj169PCxUpHEoiCQuFFSUtIsBPLy8rjjjjuYMWMGZ599ttb4RTyis4Ykruzbt48XXniBlStXsmrVKt577z2g4QZzU6ZMYerUqVx00UWMHTtW+/9FOkAPppHAKi0tZdWqVaxcuZKVK1fywQcfAA1nBIVDYerUqZx++ukKBpF2KAgkITjnKCkpiWwtrFy5MnJVcWFhYSQULrroIoYMGeJztSLxRUEgCck5x9atWyOhsGrVKg4dOhSZv2jRIr70pS/5WKFI/GgvCHT0TQLLzDjjjDOYO3cu//qv/8q1117bbP5LL73kU2UiwRLVWUNm9kUaHh7TH7DGl3PO5XpYm0ibtm/fzvLly3n++edZtWoVlZWVpKen8/nPf56ioiKmTZvGmDFj/C5TJBCiPX30P4Dpzrl3vCxGpC2VlZW8+OKLLF++nOXLl0fuITRixAjmzp3LtGnTuPDCC8nKyvK5UpHgiTYI9ikEpDs553j33Xcja/0vvvgiVVVVZGZmMnXqVL7+9a8zbdo0Ro4c6XepIoHXbhA07hICWGdmTwPPASfC851zz3pYmySpV199lU9+8pOR8TFjxnDzzTdTVFTEBRdcQGZmpo/ViSSeU20RTG8yXAl8vsm4AxQEEnMlJSWR4ZycHK6//npmzZrF8OHDfaxKJHG1e/qomX0ZWOGcK+u+ktqn00eTw759+1i8eDELFy5kzZo1AEyePJnZs2dz9dVXM2DAAJ8rFAmWTl9HYGb/F7gESAP+B1gOvO58vPhAQZB8du7cyaJFi1i4cCFvvvkmoVCIqVOnMnv2bL74xS/qOcUiUejyBWVmlgNcDBQB5wHvAM/TsLWwL4a1npKCILlt2bKFhQsXsnDhQrZv305aWhrTpk1j9uzZTJ8+XWcNibQhFkHwP8DDzrmljeNnAhuBlc65S2JZ7KkoCAQazipat24dCxcu5Omnn2b37t0APPnkkyddWCYisbmyeDjwTTO7F8A5txl4u7tDQCTMzBg7dixnnnkmBQUFkenHjx/3sSqRYIo2CA4DnwVOM7O/mFkeDWcNiXQr5xyrV69m7ty5DBw4kJtuuony8nLmz5/P3r17mTdvnt8ligROtBeUmXOuFviamc0BVgM6QifdZu/evfzmN79hwYIFvPvuu2RnZzN79mxuvPFGPZRepIuiDYLHwgPOuV+b2VvArd6UJNKgpqaGZcuW8cQTT7Bs2TLq6ur4zGc+w7e//W2uuuoqHRgWiZGogsA594sW4+uBuZ5UJALs2rWLT3ziE+zbt48BAwZw1113ccMNNzB69Gi/SxNJOLoNtcSllJQUUlMb1lN69+7N2Wefzemnn+5zVSKJSUEgcWnAgAGUlJTw5JNPYmZ85StfYfTo0fz85z+nqqrK7/JEEoqCQOJWWloa1157LW+99RbPPfcc/fv352tf+xrDhg3jgQce4MiRI36XKJIQFAQS90KhEDNnzmTNmjWsWrWK8ePH8+1vf5vCwkK+853vsG9ft17cLpJwFAQSGGbGlClTWLFiBevXr+eSSy7hhz/8IUOGDGHjxo1+lycSWAoCCaSKigr27t0LQF5eHvn5+T5XJBJcCgIJlPXr1zNt2jQuuOACiouLefTRR9m9ezdDhgzxuzSRwPI0CMysyMy2mVmxmX2rjTZXm9kWM9tsZr/3sh4Jri1btnDllVcyYcIEXn/9dR588EGKi4u55ZZb6NGjh9/liQRatFcWd5iZpQCPAJ8DSoG1ZrbEObelSZtRwLeB851zH5lZf6/qkWDasWMH9913H7/73e/Iysri3nvv5fbbbyc3N9fv0kQShmdBQMNzC4qdczsAzGwRMBPY0qTNV4FHnHMfATjn9ntYjwTMW2+9xbnnnktdXR133XUX3/zmN3UsQMQDXu4aGgx82GS8tHFaU6OB0Wb2spm9amZFrf0gM5tnZuvMbN2BAwc8KlfizcCBAyO3mN69ezc+PhhPJKH5fbA4FRgFTAFmA780s14tGznnHnfOTXDOTejXr183lyh+yc/PZ+vWrdxzzz0sXryYM844g6eeekqBIBJjXgbBLqCgyfiQxmlNlQJLnHM1zrkS4F0agkEEgPT0dL73ve+xYcMGRo0axXXXXUdRURElJSV+lyaSMLwMgrXAKDMbbmY9gFnAkhZtnqNhawAzy6dhV9EOD2uSgDrzzDNZvXo1P/3pT3nllVc466yzeOihh6itrfW7NJHA8ywIGh9kcxuwgoaH3S92zm02s/vNbEZjsxVAmZltAVYBdznnyryqSYItJSWF2267jS1btjB69GjuvPNOvvWtVs9KFpEO8PKsIZxzy4BlLabd02TYAbc3vkSi8sorr7Bz504APvvZz/pbjEgC8PtgsUjUjhw5wrXXXsusWbMYO3YsxcXFTJs2ze+yRAJPQSCB8NJLLzF+/HgWLlzIfffdx0svvcSIESP8LkskISgIJK5VV1dz9913M2XKFFJTU1m9ejX33ntv5OllItJ1+t8kce173/seP/jBDygqKmLx4sXk5OT4XZJIwtEWgcS1Sy+9lFAoxFtvvcWOHTqzWMQLCgKJa+effz5vvPEGZsanP/1pli1bduoPiUiHKAgk7o0fP57XXnuNUaNGMX36dB599FG/SxJJKAoCCYRBgwbxz3/+k8suu4xbb72VO+64g7q6Or/LEkkICgIJjOzsbP70pz8xadIkHn74YX7/ez3HSCQWFAQSKKtXr2bz5s0UFhbqqmKRGFEQSGA8//zzFBUVUVBQwCuvvMKgQYP8LkkkISgIJBD++Mc/MmPGDM444wxefPFFBg9u+YwjEeksBYHEvaeeeoqrr76aiRMnsnLlSvRwIpHY0pXFEtcqKiqYM2cO9fX1rFixguzsbL9LEkk42iKQuJaVlcWXv/xlAJYuXepzNSKJSUEgce+Xv/wl559/Ptdffz0vv/yy3+WIJBwFgcS9jIwMnnvuOQoLC5k5cybvvfee3yWJJBQFgQRCfn5+5D5Dl156KQcPHvS5IpHEoSCQwBg5ciQ/+clPKC4u5he/+IXf5YgkDAWBBEZlZSXz58+nV69efOUrX/G7HJGEodNHJRCcc9x8881s2rSJpUuXMmzYML9LEkkYCgIJhEceeYSnnnqK+++/Xw+sF4kx7RqSuLd69Wr+/d//nenTp3P33Xf7XY5IwlEQSNx7+OGHqa2t5Ve/+hWhkP5kRWJN/6sk7s2bNw+An/3sZz5XIpKYFAQS94qKirjuuuv44Q9/yKZNm/wuRyThKAgkEH784x/Tt29fbrjhBmpqavwuRyShKAgkEPr06cOjjz7Khg0bmD9/vt/liCQUBYEExhVXXEGvXr145pln/C5FJKEoCCQwnnvuOQ4fPsycOXP8LkUkoXgaBGZWZGbbzKzYzL7VyvxCM1tlZhvM7E0zu9TLeiS4Kisr+cY3vsHHPvYxbrnlFr/LEUkonl1ZbGYpwCPA54BSYK2ZLXHObWnS7LvAYufcz81sHLAMGOZVTRJcP/jBD/jggw/47W9/S2qqLogXiSUvtwjOA4qdczucc9XAImBmizYOyG0czgN2e1iPBFRxcTEPPvgg11xzDZ/5zGf8Lkck4XgZBIOBD5uMlzZOa+o+4BozK6Vha+DfWvtBZjbPzNaZ2boDBw54UavEsRdeeIHq6mrmzp3rdykiCcnvg8WzgV8754YAlwJPmdlJNTnnHnfOTXDOTejXr1+3Fyn+mjFjBj179mTBggV+lyKSkLwMgl1AQZPxIY3TmroRWAzgnFsDZAD5HtYkAdS/f39uvfVWfv/737Nt2za/yxFJOF4GwVpglJkNN7MewCxgSYs2HwCfBTCzM2gIAu37kZPcddddZGZmcv/99/tdikjC8SwInHO1wG3ACuAdGs4O2mxm95vZjMZmdwBfNbNNwEJgjnPOeVWTBFe/fv247bbbWLhwIe+8847f5YgkFE+PETjnljnnRjvnRjjn/l/jtHucc0sah7c45853zo13zp3jnPubl/VIsE2bNg3nHG+88YbfpYgkFL8PFotE7bHHHiMnJ4fLL7/c71JEEoqCQAJhx44dLF68mJtvvpm8vDy/yyVLgV0AAAhYSURBVBFJKAoCCYSHHnqI1NRUvvGNb/hdikjCURBI3Dtw4AALFizg2muvZdCgQX6XI5JwFAQS90pKSqiqqtLtJUQ8oiCQuDdhwgQKCwv57W9/63cpIglJQSBxLxQKcdNNN/GPf/yD7du3+12OSMJREEggzJ07l5SUFH75y1/6XYpIwlEQSCAMHjyYyy+/nP/+7/+murra73JEEoqCQAJj7NixHDx4kOPHj/tdikhCURBIINTX17No0SIuvvhiXVAmEmMKAgmEF198kffff18PrhfxgIJAAuHXv/41ubm5fOELX/C7FJGEoyCQuHfs2DGeeeYZZs2aRWZmpt/liCQcBYHEvYqKCiorKykoKDh1YxHpMAWBxL3TTjuNSZMm8eyzz/pdikhCUhBIIFx11VVs2LBBVxaLeEBBIIFw5ZVXAvCHP/zB50pEEo+CQAJh6NChTJo0SUEg4gEFgQRGr169qKio8LsMkYSjIJBAOHLkCCtXrmT69Ol+lyKScBQEEghLly6lpqaGK664wu9SRBKOgkAC4dlnn2XgwIFMnjzZ71JEEo6CQOJeVVUVy5cv54orriAU0p+sSKzpf5XEPeccdXV1pKen+12KSEJSEEjcy8zMZOrUqfzlL3/BOed3OSIJR0EggTB9+nSKi4vZtm2b36WIJBwFgQRC+LTRpUuX+lyJSOJREEggFBQUMHLkSF5++WW/SxFJOAoCCYyMjAyOHj3qdxkiCUdBIIHw/vvv8/bbb1NUVOR3KSIJx7MgMLMFZrbfzN5uY76Z2X+ZWbGZvWlmH/eqFgm+8LEB3WJCJPa83CL4NdDe6ts0YFTjax7wcw9rkYD761//yqhRoxg9erTfpYgkHM+CwDn3T+BQO01mAk+6Bq8CvcxsoFf1SHBVVFSwcuVKLr/8cr9LEUlIfh4jGAx82GS8tHHaScxsnpmtM7N1Bw4c6JbiJH4cPnyYK664QjecE/FIqt8FRMM59zjwOMCECRN0aWmSGTx4MAsXLvS7DJGE5ecWwS6goMn4kMZpIiLSjfwMgiXAdY1nD00Gjjjn9vhYj4hIUvJs15CZLQSmAPlmVgrcC6QBOOceA5YBlwLFQCVwg1e1iIhI2zwLAufc7FPMd8CtXv1+ERGJjq4sFhFJcgoCEZEkpyAQEUlyCgIRkSRnQXv0n5kdAN7v5MfzgYMxLCfeaPmCTcsXbPG+fEOdc/1amxG4IOgKM1vnnJvgdx1e0fIFm5Yv2IK8fNo1JCKS5BQEIiJJLtmC4HG/C/CYli/YtHzBFtjlS6pjBCIicrJk2yIQEZEWFAQiIkkuIYPAzIrMbJuZFZvZt1qZn25mTzfOf83MhnV/lZ0XxfLNMbMDZrax8XWTH3V2hpktMLP9ZvZ2G/PNzP6rcdnfNLOPd3eNXRHF8k0xsyNNvrt7urvGrjCzAjNbZWZbzGyzmf2fVtoE9juMcvmC9x065xLqBaQA24HTgR7AJmBcizZfAx5rHJ4FPO133TFevjnAz/yutZPLdwHwceDtNuZfCiwHDJgMvOZ3zTFevinAX/2uswvLNxD4eONwDvBuK3+fgf0Oo1y+wH2HibhFcB5Q7Jzb4ZyrBhYBM1u0mQn8pnH4GeCzZmbdWGNXRLN8geWc+ydwqJ0mM4EnXYNXgV5mNrB7quu6KJYv0Jxze5xzbzQOHwXe4eRnkQf2O4xy+QInEYNgMPBhk/FSTv6iIm2cc7XAEaBvt1TXddEsH8C/NG52P2NmBa3MD6polz/IPmlmm8xsuZmd6XcxndW4y/Vc4LUWsxLiO2xn+SBg32EiBoHAX4BhzrmPAX/nf7d+JP69QcM9YcYDPwWe87meTjGzbOCPwDecc+V+1xNrp1i+wH2HiRgEu4Cma8BDGqe12sbMUoE8oKxbquu6Uy6fc67MOXeicfRXwCe6qbbuEM33G1jOuXLn3LHG4WVAmpnl+1xWh5hZGg2d5O+cc8+20iTQ3+Gpli+I32EiBsFaYJSZDTezHjQcDF7Sos0S4PrG4SuBla7xKE8AnHL5WuxvnUHDfsxEsQS4rvHMk8nAEefcHr+LihUzGxA+XmVm59HwfzQoKyk01v4E8I5z7uE2mgX2O4xm+YL4HXr2zGK/OOdqzew2YAUNZ9gscM5tNrP7gXXOuSU0fJFPmVkxDQfuZvlXccdEuXxfN7MZQC0NyzfHt4I7yMwW0nDWRb6ZlQL3AmkAzrnHgGU0nHVSDFQCN/hTaedEsXxXAreYWS1wHJgVoJUUgPOBa4G3zGxj47TvAIWQEN9hNMsXuO9Qt5gQEUlyibhrSEREOkBBICKS5BQEIiJJTkEgIpLkFAQiIklOQSAikuQUBCIiSU5BIBIjZvaYmZ3vdx0iHaULykRipPFK00845+r8rkWkI7RFINIBZraw8el2r5vZ+2Z2WeP0M2h4SEmBmf2zcdrHzcyZWb6ZpZjZW2bW08/6RVqjIBDpmPHADufcecBXaLhXEMA04HngMJDdOO3fgFeBXsBlwD+cc5XdW67IqSkIRKJkZhlAP+B7jZO2AL0bhy+hIQjKgZ6Ntx0eCLzc2GYe8PNuLVgkSgoCkeidBbznnKtqHP84sKlxd08v59xu51w94ICbaLjL7VEatiJSnHPv+lG0yKkoCESiNx4oNLMMM8uiYcvgx8BUYFWTdvU0PAfiTzRsIdwBPNbNtYpETUEgEr3xwLM0PKN2LfBz59zL/O/xgbAaYHnj87DLgZ7AX7u5VpGo6fRRkSiZ2YvAPOfcthbT3wAmOedq/KlMpGsUBCJRanyiWGHjcQCRhKEgEBFJcjpGICKS5BQEIiJJTkEgIpLkFAQiIklOQSAikuQUBCIiSU5BICKS5P4//NMKtjZGa44AAAAASUVORK5CYII=\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 3(b)"
      ],
      "metadata": {
        "id": "rtxvdCroiuq_"
      },
      "id": "rtxvdCroiuq_"
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "time = np.linspace(0,10,100)\n",
        "pressure = np.array(p_out)\n",
        "plt.xlabel(r\"$t$\")\n",
        "plt.ylabel(r\"$p/w$ at $z/h=0.5$\")\n",
        "plt.plot(time, pressure/10, 'k')"
      ],
      "metadata": {
        "id": "BgAwNRBDfINY",
        "outputId": "b13c2345-a635-4c07-e2cf-1451897f9cfd",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 298
        }
      },
      "id": "BgAwNRBDfINY",
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[<matplotlib.lines.Line2D at 0x7efc6f28ce50>]"
            ]
          },
          "metadata": {},
          "execution_count": 25
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEICAYAAABMGMOEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3RUdd7H8fc3oaMCCrJK1wWkSQsl9BUVAixFBWkiritnZVEQLFhX1EdEmgXUdZVVqUZkAWEFQUACEiSRJgiKIEpx5bGtgijl+/yRgSci6Ixk5iaTz+ucOTNz753czxwxn9z6M3dHRETk1yQEHUBERPIGFYaIiIRFhSEiImFRYYiISFhUGCIiEpYCQQeIltKlS3vlypWDjiEikqdkZmb+r7uXOdm8uC2MypUrk5GREXQMEZE8xcx2nmqedkmJiEhYVBgiIhIWFYaIiIRFhSEiImFRYYiISFhiVhhmNsnMPjez904x38zsCTPbZmYbzKxBtnnXmtmHoce1scosIiL/L5ZbGC8A7X9hfgpQNfQYADwNYGZnA38DmgCNgb+ZWamoJhURkZ+J2XUY7r7czCr/wiJdgJc8637r6WZW0szOA9oAi9z9SwAzW0RW8UyPRs79+/czatSoU843s+PPxx4JCQkkJCSQmJhIQkICBQsWpFChQhQqVIjChQtTtGhRihYtSvHixSlRogQlS5akVKlSlChR4vjPExHJ7XLThXvlgE+zvd8Vmnaq6T9jZgPI2jqhYsWKvynEgQMHeOihh046L6fHDilYsCBlypShbNmyVKxYkUqVKlGpUiWqV69O9erVqVy5MgUK5Kb/RCKSn8XVbyN3fxZ4FiApKek3/XYvU6YMR48eDWddxx9Hjx49/jhy5AiHDh06/jh48CDff/89Bw4c4LvvvuObb77hm2++4csvv2Tfvn18/vnn7N27l48++oglS5bw7bffHl9H4cKFqVOnDg0aNKBBgwY0a9aMWrVqkZCgcxVEJPZyU2HsBipke18+NG03Wbulsk9fFrNUp3BsdxRAYmJijv3cL7/8kq1bt7JlyxY2b97M2rVreeWVV3j22WcBKFWqFC1atODSSy8lJSWFqlWr5ti6RUR+icVyiNbQMYx57l77JPM6AoOADmQd4H7C3RuHDnpnAsfOmnoXaHjsmMapJCUlebzcS8rd2b59OytWrCAtLY233nqLbdu2AXDhhRfSpUsXunfvTuPGjbX1ISKnxcwy3T3ppPNiVRhmNp2sLYXSwH/IOvOpIIC7P2NZf65PIOuA9gHgOnfPCH32T8BdoR/1P+7+z19bXzwVxsl89NFHLFiwgPnz57N48WIOHTpE+fLl6dWrF/3796dmzZpBRxSRPChXFEasxXthZPf1118zb948UlNTef311zl8+DCNGjXihhtuoHfv3hQvXjzoiCKSR/xSYWj/RRwoWbIkffv2Ze7cuezevZtx48Zx8OBBBgwYQPny5Rk6dCjbt28POqaI5HEqjDhz7rnncsstt7B+/XrS0tJo3749Tz75JFWrVqVnz568++67QUcUkTxKhRGnzIwWLVowffp0du7cya233srrr79Ow4YNadeuHenp6UFHFJE8RoWRD5x//vmMGjWKTz75hJEjR/Luu++SnJxMx44dyczMDDqeiOQRKox8pESJEgwfPpwdO3YwcuRI0tPTSUpKonfv3uzYsSPoeCKSy6kw8qEzzjjjeHHcc889zJ49m+rVqzN06FC+/vrroOOJSC6lwsjHzjrrLB588EE+/PBD+vXrx2OPPUa1atV4/vnnw7o9iojkLyoMoVy5cjz33HNkZmZSrVo1/vznP9OkSRPyy3UsIhIeFYYcV79+fdLS0pg6dSq7du2iSZMm3Hzzzfz3v/8NOpqI5AIqDPkJM6N3795s2bKFG2+8kQkTJlCjRg3mzJkTdDQRCZgKQ06qRIkSTJgwgfT0dEqXLk3Xrl3p2bMn+/btCzqaiAREhSG/qHHjxqxZs4YHHniAWbNmUaNGDWbOnBl0LBEJgApDflWhQoW49957Wbt2LVWqVKF79+706dOHL7/8xTvMi0icUWFI2GrVqsXbb7/NiBEjSE1NpU6dOixatCjoWCISIyoMiUjBggW57777WL16NSVKlODyyy9n6NChHDx4MOhoIhJlKgz5TRo0aEBGRgZ//etfGT9+PE2aNGHz5s1BxxKRKFJhyG9WrFgxJkyYwLx589i7dy9JSUlMmjSJeB2USyS/U2HIaevYsSPr1q0jOTmZ66+/nj59+uhiP5E4pMKQHHH++efzxhtv8NBDD/Hyyy+TlJTE+vXrg44lIjlIhSE5JjExkbvvvpulS5fy3Xff0bRpU55//nntohKJEyoMyXGtWrVi7dq1NG/enD//+c9cd911fP/990HHEpHTpMKQqChbtiwLFy7kvvvu48UXX6RZs2Zs37496FgichpUGBI1iYmJjBgxgvnz5/Pxxx/TsGFD/v3vfwcdS0R+IxWGRF2HDh3IzMykcuXKdOrUiYceekgDNInkQSoMiYkLLriAlStX0rt3b+69916uvPJKnXorkseoMCRmihUrxuTJkxk/fjyvvfYaycnJbNu2LehYIhImFYbElJkxZMgQFi1axH/+8x8aNWrEG2+8EXQsEQmDCkMC8Yc//IE1a9ZQoUIFUlJSGDdunK7XEMnlVBgSmCpVqvD222/TrVs3hg0bxvXXX88PP/wQdCwROYWYFoaZtTezrWa2zcyGn2R+JTN708w2mNkyMyufbd4RM1sXesyNZW6JnjPOOIPU1FTuu+8+/vnPf9K2bVs+//zzoGOJyEnErDDMLBGYCKQANYFeZlbzhMXGAC+5+8XAA8DIbPO+d/d6oUfnmISWmEhISGDEiBHMmDGDzMxMGjVqxMaNG4OOJSIniOUWRmNgm7tvd/cfgRlAlxOWqQksCb1eepL5Eseuvvpq0tLSOHz4MM2aNWP+/PlBRxKRbGJZGOWAT7O93xWalt164IrQ627AmWZ2Tuh9ETPLMLN0M+t6shWY2YDQMhn79u3LyewSI0lJSbzzzjtUq1aNzp07M378eB0MF8klcttB71uB1ma2FmgN7AaOhOZVcvckoDfwmJldeOKH3f1Zd09y96QyZcrELLTkrHLlyrF8+XK6du3K0KFDGThwIIcOHQo6lki+F8vC2A1UyPa+fGjace6+x92vcPf6wN2haV+HnneHnrcDy4D6McgsASlevDivvPIKw4cP55lnnqFTp0588803QccSyddiWRhrgKpmVsXMCgE9gZ+c7WRmpc3sWKY7gUmh6aXMrPCxZYDmgAaQjnMJCQmMHDmSSZMmsWTJEpo1a8bHH38cdCyRfCtmheHuh4FBwELgfSDV3TeZ2QNmduyspzbAVjP7ACgL/E9oeg0gw8zWk3Uw/BF3V2HkE9dddx1vvPEGe/bsoUmTJrzzzjtBRxLJlyxeDygmJSV5RkZG0DEkB23ZsoUOHTrw2WefMWXKFK644opf/5CIRMTMMkPHi38mtx30Fjmliy66iPT0dOrWrctVV12lM6hEYkyFIXnKueeey5IlS7jiiisYOnQoN998M0eOHPn1D4rIaVNhSJ5TtGhRUlNTufXWW5kwYQLdunVj//79QccSiXsqDMmTEhISGD16NBMnTmT+/Pm0adOGzz77LOhYInFNhSF52sCBA5k9ezabN28mOTmZ999/P+hIInFLhSF53h//+EfeeustDhw4QPPmzUlLSws6kkhcUmFIXEhKSiI9PZ1zzz2XSy+9lNTU1KAjicQdFYbEjWMDMjVu3Jirr76aMWPG6LRbkRykwpC4cvbZZ7No0SK6d+/ObbfdxpAhQ3TarUgOKRB0AJGcVqRIEWbMmEGFChUYN24cu3btYsqUKRQtWjToaCJ5mrYwJC4lJCQwduxYxo0bx7/+9S8uu+wyvvjii6BjieRpKgyJa7fccgsvv/wya9asoXnz5rrbrchpUGFI3OvevTuLFi3iP//5D8nJyaxduzboSCJ5kgpD8oVWrVqxYsUKChYsSKtWrVi0aFHQkUTyHBWG5Bu1atVi1apVXHDBBXTo0IHJkycHHUkkT1FhSL5ybLzwVq1a0a9fP0aOHKlrNUTCpMKQfKdEiRK8/vrr9O7dm7vuuotBgwbpWg2RMOg6DMmXChUqxOTJkylXrhyjR49mz549TJs2TddqiPyCiLcwzKyFmQ01s8ujEUgkVhISEnj00Ud5/PHHmTNnDpdeeqmu1RD5Bb9aGGb2TrbXNwATgDOBv5nZ8ChmE4mJm2++mdTUVDIzM2nRogU7d+4MOpJIrhTOFkbBbK8HAJe5+wjgcqBPVFKJxNhVV13FG2+8wWeffUZycjLr1q0LOpJIrhNOYSSYWSkzOwcwd98H4O77gcNRTScSQ8eu1UhMTKRVq1YsXrw46EgiuUo4hVECyAQygLPN7DwAMzsDsChmE4m5Y9dqVKpUiQ4dOjB16tSgI4nkGr9aGO5e2d0vcPcqoee9oVlHgW7RjScSe+XLlyctLY3mzZvTt29fRo8erWs1RDiN6zDc/YC778jJMCK5RcmSJVmwYAE9evTg9ttv17gaIpzmdRhm9jt3/yynwojkJoULF2b69OmUK1eO8ePHs2fPHiZPnkyRIkWCjiYSiNO90vv5HEkhkkslJCQwbtw4xo4dy8yZM7n88sv56quvgo4lEojTKgx375hTQURys6FDhzJ9+nRWr15NixYt+OSTT4KOJBJzYRWGmV1kZneY2ROhxx1mViPSlZlZezPbambbTnbRn5lVMrM3zWyDmS0zs/LZ5l1rZh+GHtdGum6R09WzZ08WLlzI7t27SU5OZsOGDUFHEompcK70vgOYQdYptO+EHgZMj+RKbzNLBCYCKUBNoJeZ1TxhsTHAS+5+MfAAMDL02bOBvwFNgMZkXWVeKtx1i+SUNm3asGLFCsyMli1bsmTJkqAjicRMOFsY1wON3P0Rd58SejxC1i/u6yNYV2Ngm7tvd/cfySqhLicsUxM49n/g0mzz2wGL3P1Ld/8KWAS0j2DdIjmmdu3apKenU7FiRdq3b8+0adOCjiQSE+EUxlHg/JNMPy80L1zlgE+zvd8VmpbdeuCK0OtuwJmhK8zD+SxmNsDMMswsY9++fRFEE4lM9ms1+vTpw6hRo3SthsS9cE6rHQK8aWYf8v+/tCsCvwcG5XCeW4EJZtYfWA7sBsI++d3dnwWeBUhKStL/vRJVx67V6N+/P8OHD+fTTz/l8ccfJzExMehoIlHxq4Xh7gvMrBpZu5SO/VW/G1jj7pFcybQbqJDtffnQtOzr2kNoCyN065Er3f1rM9sNtDnhs8siWLdIVBQuXJipU6dSoUIFRo8eza5du5g2bRrFihULOppIjgvrwj13Pwqkn+a61gBVzawKWUXRE+idfQEzKw18GVrfncCk0KyFwMPZDnRfHpovErhj42pUqFCBwYMH07ZtW1577TVKly4ddDSRHBWzIVrd/TBZu7AWAu8Dqe6+ycweMLPOocXaAFvN7AOgLPA/oc9+CTxIVumsAR4ITRPJNW666SZeffVV1q1bR3JyMh999FHQkURylEVyoM7MLnH3Jceeo5jrtCUlJXlGRkbQMSQfevvtt+ncuTMJCQnMmzePxo0bBx1JJGxmlunuSSebF+kWxpgTnkXkBM2aNePtt9/mjDPOoE2bNsydOzfoSCI54rfuktI4GCK/oFq1aqxatYpatWrRrVs3nn766aAjiZy2mB3DEMlvypYty7Jly+jYsSMDBw7kjjvu4OjRSC5dEsldVBgiUVS8eHFmzZrFjTfeyKOPPkrv3r05ePBg0LFEfpPTGg9DRH5dgQIFmDhxIpUrV+aOO+5gz549zJ49m7PPPjvoaCIRiXQL47vQ87c5HUQknpkZt99+OzNmzGD16tU0a9aM7du3Bx1LJCIRFYa7t8r+LCKRufrqq1m8eDGff/45ycnJvPPOO0FHEgmbjmGIxFjLli1ZtWoVxYsXp02bNsyePTvoSCJhUWGIBKB69eqkp6dTp04drrjiCh5//PGgI4n8KhWGSEDOPfdcli5dSteuXRkyZAiDBw/myJFI7ucpElthF4aZXWJmz5vZWDO7zswamlnhaIYTiXfFihXjlVdeYejQoTzxxBNcccUV7N+/P+hYIicVyRbGJOA1su5aewFwH7ApGqFE8pPExETGjh3LhAkTmDdvHq1atWLv3r1BxxL5mUiuw9jp7seOzr0SjTAi+dlf//pXKlWqRM+ePWnSpAnz58+nTp06QccSOe5XtzDM7CUzGwKkm9nQGGQSybc6depEWloaR44coXnz5ixcuDDoSCLHhbNL6gWybjZYFrjGzHaa2Vwze9DMukc1nUg+VL9+fVavXs0FF1xAx44d+fvf/x50JBEgvML4HnjM3fu7e33gQuBuYCtZw7aKSA4rX748aWlptGvXjr/85S8MGzZMZ1BJ4MIpjH5AppnNMLP+QGl33+juU9z9tujGE8m/zjzzTObMmcOgQYMYN24cV155pc6gkkD9amG4+43u3gC4HygFvGBmq8zsYTNrZWaJ0Q4pkl8VKFCAJ598kscff5zXXnuNVq1asXv37qBjST4VyWm1tcnaNdUeuARYAXQHVkcjmIj8v5tvvpm5c+fywQcf0KRJE9auXRt0JMmHIimMycA0M0t09+/d/d/Au6ca+1VEclbHjh1ZuXIlCQkJtGzZUkO/SsxFUhhbgLeAV82sYGjaTTkfSURO5eKLL2b16tXUrFmTrl27MmbMGNw96FiST0RSGO7uzwCzgLlmVhSN7S0Sc+eddx7Lli3jqquu4rbbbuOGG27gxx9/DDqW5AORXOn9FYC7v2RmB4D5QLGopBKRX1SsWDFmzJjBRRddxIMPPsi2bdt49dVXOeecc4KOJnEs7C0Md2+b7fVMYBygf50iAUlISOCBBx5gypQppKen06RJE7Zs2RJ0LIljv/n25u4+z91L52QYEYlcnz59WLp0Kd9++y1NmzZl0aJFQUeSOKXxMETiwLHhXitWrEhKSgoTJ04MOpLEIRWGSJyoVKkSK1eupGPHjgwaNIiBAwdy6NChoGNJHIlkAKVR4UwTkeCceeaZzJo1izvuuIOnn36adu3a8cUXXwQdS+JEJFsYl51kWkpOBRGRnJGYmMgjjzzCiy++yMqVK2nSpAmbN28OOpbEgXDGw7jRzDYC1c1sQ7bHDmBjJCszs/ZmttXMtpnZ8JPMr2hmS81sbWgdHULTK5vZ92a2LvR4JpL1iuRH/fr1Y9myZXz33Xc0bdqU+fPnBx1J8rhwtjCmAX8E5oaejz0aunufcFcUuknhRLK2SmoCvcys5gmL3QOkhm6j3hN4Ktu8j9y9Xujxl3DXK5KfJScns2bNGn7/+9/zxz/+kdGjR+vKcPnNwrlb7Tfu/rG79wL+S9ZASpWA2mbWKoJ1NQa2uft2d/8RmAF0OXF1wFmh1yWAPRH8fBE5iQoVKpCWlsZVV13F7bffzrXXXsvBgweDjiV5UCQHvf8MLAcWAiNCz/dHsK5ywKfZ3u8KTcvufqCvme0C/s1P71VVJbSr6i0za3mKjAPMLMPMMvbt2xdBNJH4Vrx4cV5++WUefPBBJk+eTOvWrdmzR3+PSWQiOeg9GGgE7HT3PwD1ga9zOE8v4AV3Lw90ACabWQKwF6gY2lU1lKy75p514ofd/Vl3T3L3pDJlyuRwNJG8zcy45557+Ne//sXmzZtJSkoiPT096FiSh0RSGAfd/SCAmRV29y1A9Qg+vxuokO19+dC07K4HUgHcfRVQhKwR/n5w9y9C0zOBj4BqEaxbREK6du3KqlWrKFKkCK1bt+aFF14IOpLkEZEUxi4zKwnMBhaZ2RxgZwSfXwNUNbMqZlaIrIPaJ97Q/xOgLYCZ1SCrMPaZWZljI/uZ2QVAVWB7BOsWkWxq167NmjVraNmyJddddx2DBw/WRX7yq8K+W627dwu9vN/MlpJ1UHpBBJ8/bGaDyDr2kQhMcvdNZvYAkOHuc4FhwD/M7BayDoD3d3cPHVx/wMwOAUeBv7j7l+GuW0R+7pxzzmHBggXcdtttPPbYY2zcuJHU1FRKl9Yt4uTkLF5PsUtKSvKMjIygY4jkCS+99BIDBgzgvPPOY/bs2dStWzfoSBIQM8s81UiqupeUiNCvXz+WL1/OoUOHSE5OZsaMGUFHklxIhSEiADRu3JiMjAwaNmxIr169uP322zl8+HDQsSQXieQ6jCLRDCIiwfvd737Hm2++yY033sjo0aNJSUnRzQvluEi2MN4xs7Fm9vuopRGRwBUqVIinnnqK5557juXLl5OUlMS6deuCjiW5QCSFUQ9YBow3s/lm1snMLDqxRCRo119/PWlpaRw6dIhmzZoxderUoCNJwCIpjJLAJrJuCzILeBRdCyES1xo3bkxmZiaNGjWib9++DBkyRNdr5GORFMb/ApOBHsD5wLPAg9EIJSK5R9myZVm8eDFDhgzh8ccfp23btnz22WdBx5IARFIYScAHQB1gM/CEu0+KSioRyVUKFizI+PHjmTp16vEzqd5+++2gY0mMhV0Y7v6uu18H9AV+Dyw3s7uilkxEcp3evXuTnp5O0aJFad26NRMmTND4GvlIJKfVvmVmGUAacC1ZxzSuilYwEcmdLr74YjIyMmjfvj033XQT11xzDfv37w86lsRAOEO0JofOhupH1o0Ba7l7TXev7e4Nop5QRHKdkiVLMmfOHB588EGmTZtG06ZN+fDDD4OOJVEWzhZGPyATGAV0I2vEPRHJ5xISErjnnntYsGABe/bsISkpidmzZwcdS6IonCFabwxtSdwPlAJeMLNVZvawmbU6dttxEcmfLr/8ct59912qVatGt27ddEuROBbJQe8t7j7e3dsDlwArgO7A6miFE5G8oVKlSqxYseL4LUXatm3L3r17g44lOSySg94pZrbazLYCLwJfuftNp7oNrojkL4ULF+app55i8uTJZGRk0KBBA956662gY0kOiuQ6jKfIGk+7KVkX7Y02s15RSSUieVbfvn1ZvXo1JUqU4JJLLuGRRx7h6NGjQceSHBBJYXzu7ivd/St3Xwy0A+6OUi4RycOODQHbvXt37rzzTrp06cKXX2qQzLwuksLYYWYPhcbjBjgE6MiWiJzUmWeeyfTp03nyySdZuHAh9evXZ/VqHfLMyyIpjKNknVb7qZmtALYBy8ysalSSiUieZ2YMGjSIlStXYma0bNmSJ554QleH51GRnCXV291rARWBwWSdZmvAP8zsk+jEE5F40KhRI9auXUv79u0ZPHgwV111FV9//XXQsSRCEQ/R6u4/uHumu09y98Hu3sbdK0YjnIjEj1KlSjFnzhzGjBnD3LlzadiwIZmZmUHHkghoTG8RiRkzY9iwYSxfvvz4wExPPvmkdlHlESoMEYm55ORk1q5dS7t27bj55pu1iyqPUGGISCDOOecc5syZw9ixY5k7d67OosoDVBgiEhgzY+jQoaxYsQKAFi1aMGbMGF3ol0upMEQkcE2aNGHt2rV07tyZ2267jU6dOrFv376gY8kJVBgikiuULFmSmTNnMnHiRJYsWULdunVZunRp0LEkGxWGiOQaZsbAgQNZvXo1Z511Fm3btuXee+/V7dJziZgWhpm1N7OtZrbNzIafZH5FM1tqZmvNbIOZdcg2787Q57aaWbtY5haR2Kpbty6ZmZn079+fhx56iNatW7Nz586gY+V7MSuM0EBLE4EUoCbQy8xqnrDYPUCqu9cHepJ1h1xCy/UEagHtgac0cJNIfCtevDiTJk1i2rRpbNy4kbp16/LKK68EHStfi+UWRmNgm7tvd/cfgRlAlxOWceCs0OsSwJ7Q6y7AjNBV5jvIuo9V4xhkFpGA9erVi3Xr1lG9enV69OjBDTfcwP79+4OOlS/FsjDKAZ9me78rNC27+4G+ZrYL+DdwUwSfxcwGmFmGmWXoDAuR+HHBBRewYsUKhg8fzvPPP0/Dhg1Zu3Zt0LHyndx20LsX8IK7lwc6AJPNLJIbJD7r7knunlSmTJmohRSR2CtYsCAjR45k8eLFfPvttzRt2pRx48bpmo0YimVh7AYqZHtfPjQtu+uBVAB3XwUUAUqH+VkRyQcuueQS1q9fT4cOHRg2bBgpKSkaPzxGYlkYa4CqZlYlNAhTT2DuCct8ArQFMLMaZBXGvtByPc2ssJlVAaoC78QsuYjkKqVLl2bWrFn8/e9/Jy0tjYsvvpi5c0/8dSI5LWaF4e6HgUHAQuB9ss6G2mRmD5hZ59Biw4AbzGw9MB3o71k2kbXlsRlYAPzV3Y/EKruI5D5mxoABA8jMzKR8+fJ06dKFG2+8kQMHDgQdLW5ZvN5WOCkpyTMyMoKOISIx8MMPP3DvvfcyZswYqlWrxtSpU2nYsGHQsfIkM8t096STzcttB71FRCJWuHBhHn30URYvXsx3331H06ZNGTlyJEeOaEdETlJhiEjcuOSSS9iwYQPdunXjrrvu4g9/+AMff/xx0LHihgpDROLK2Wefzcsvv8xLL73EunXruPjii3nxxRc1ql8OUGGISNwxM6655ho2bNhAvXr16N+/P927d+eLL74IOlqepsIQkbhVuXJlli5dyqOPPsrcuXOpXbs2r7/+etCx8iwVhojEtcTERG677TbWrFlD6dKl6dChA3/5y1/47rvvgo6W56gwRCRfqFu3LmvWrOHWW2/l2WefpV69eqxcuTLoWHmKCkNE8o0iRYowevRoli5dypEjR2jVqhXDhw/nhx9+CDpanqDCEJF8p3Xr1mzYsIE//elPjBo1ikaNGrFu3bqgY+V6KgwRyZfOPPNM/vGPfzBv3jz27dtHo0aNeOihhzQc7C9QYYhIvtaxY0fee+89unfvzr333ktycjKbN28OOlaupMIQkXzvnHPOYdq0aaSmprJjxw4aNGjA6NGjdWuRE6gwRERCunfvzqZNm0hJSeH222+nZcuWfPDBB0HHyjVUGCIi2ZQtW5ZZs2YxefJktmzZQt26dRk/fry2NlBhiIj8jJnRt29fNm3axGWXXcbQoUNp3bp1vt/aUGGIiJzCeeedx5w5c3jxxRfZtGkTdevWZdy4cfl2azpkqjMAAAjxSURBVEOFISLyC8yMfv36sXnzZi677DKGDRtGy5Yt2bJlS9DRYk6FISIShmNbG1OmTGHr1q3Uq1ePUaNG5avrNlQYIiJhMjP69OnDpk2b6NChA8OHDyc5OZmNGzcGHS0mVBgiIhH63e9+x6uvvsrLL7/Mzp07adiwIffffz8//vhj0NGiSoUhIvIbmBk9evRg8+bN9OjRgxEjRtCwYUPeeeedoKNFjQpDROQ0lC5dmilTpvDaa6/x1VdfkZyczLBhwzhw4EDQ0XKcCkNEJAd06tSJTZs2MWDAAMaNG0edOnV48803g46Vo1QYIiI5pESJEjz99NMsW7aMxMRELr30Uv70pz/x1VdfBR0tR6gwRERyWOvWrVm/fj133nknL730EjVq1CA1NRV3DzraaVFhiIhEQdGiRXn44YfJyMigfPnyXH311XTu3JlPP/006Gi/mQpDRCSK6tWrR3p6OmPHjmXJkiXUrFmTJ554Ik/eXkSFISISZQUKFGDo0KG89957NG/enMGDB9OsWTM2bNgQdLSIqDBERGKkSpUqvP7660ybNu34QE133HFHnjkFN6aFYWbtzWyrmW0zs+EnmT/ezNaFHh+Y2dfZ5h3JNm9uLHOLiOQUM6NXr15s2bKF/v378+ijj1K7dm0WLlwYdLRfFbPCMLNEYCKQAtQEeplZzezLuPst7l7P3esBTwKzss3+/tg8d+8cq9wiItFw9tln89xzz/HWW29RqFAh2rdvT8+ePdm7d2/Q0U4pllsYjYFt7r7d3X8EZgBdfmH5XsD0mCQTEQlIq1atWL9+PSNGjGD27NnUqFGDp59+mqNHjwYd7WdiWRjlgOznk+0KTfsZM6sEVAGWZJtcxMwyzCzdzLqe4nMDQstk7Nu3L6dyi4hEVeHChbnvvvvYuHEjSUlJDBw4kGbNmrFu3bqgo/1Ebj3o3ROY6e7Zzzur5O5JQG/gMTO78MQPufuz7p7k7kllypSJVVYRkRxRtWpVFi1axJQpU9ixYwcNGzbklltu4dtvvw06GhDbwtgNVMj2vnxo2sn05ITdUe6+O/S8HVgG1M/5iCIiwTo25saWLVu44YYbeOyxx7jooouYOXNm4FeKx7Iw1gBVzayKmRUiqxR+draTmV0ElAJWZZtWyswKh16XBpoDm2OSWkQkAKVKleKZZ55h1apVlClThu7du9OhQwc++uijwDLFrDDc/TAwCFgIvA+kuvsmM3vAzLKf9dQTmOE/rdIaQIaZrQeWAo+4uwpDROJe06ZNycjI4LHHHmPlypXUqlWLESNGcPDgwZhnsaA3caIlKSnJMzIygo4hIpJj9uzZw7Bhw5gxYwYXXnghEyZMoH379jm6DjPLDB0v/pncetBbREROcP755zN9+nQWL15MgQIFSElJ4corr+STTz6JyfpVGCIieUzbtm1Zv349I0eOZMGCBVx00UU8/PDD/PDDD1FdrwpDRCQPKly4MMOHD+f9998nJSWFu+++mzp16rBgwYKorVOFISKSh1WsWJFXX331eFGkpKTQo0ePqJyCWyDHf6KIiMRcu3bt2LhxI+PGjWP//v2YWY6vQ4UhIhInChcuzJ133hm1n69dUiIiEhYVhoiIhEWFISIiYVFhiIhIWFQYIiISFhWGiIiERYUhIiJhUWGIiEhY4vb25ma2D9h5Gj+iNPC/ORQnr8hv3zm/fV/Qd84vTuc7V3L3k45xHbeFcbrMLONU94SPV/ntO+e37wv6zvlFtL6zdkmJiEhYVBgiIhIWFcapPRt0gADkt++c374v6DvnF1H5zjqGISIiYdEWhoiIhEWFISIiYVFhnMDM2pvZVjPbZmbDg84TbWZWwcyWmtlmM9tkZoODzhQrZpZoZmvNbF7QWWLBzEqa2Uwz22Jm75tZctCZos3Mbgn9u37PzKabWZGgM+U0M5tkZp+b2XvZpp1tZovM7MPQc6mcWJcKIxszSwQmAilATaCXmdUMNlXUHQaGuXtNoCnw13zwnY8ZDLwfdIgYehxY4O4XAXWJ8+9uZuWAm4Ekd68NJAI9g00VFS8A7U+YNhx4092rAm+G3p82FcZPNQa2uft2d/8RmAF0CThTVLn7Xnd/N/T6W7J+iZQLNlX0mVl5oCPwXNBZYsHMSgCtgOcB3P1Hd/862FQxUQAoamYFgGLAnoDz5Dh3Xw58ecLkLsCLodcvAl1zYl0qjJ8qB3ya7f0u8sEvz2PMrDJQH1gdbJKYeAy4HTgadJAYqQLsA/4Z2g33nJkVDzpUNLn7bmAM8AmwF/jG3d8INlXMlHX3vaHXnwFlc+KHqjAEADM7A3gVGOLu/w06TzSZWSfgc3fPDDpLDBUAGgBPu3t9YD85tJsitwrtt+9CVlmeDxQ3s77Bpoo9z7p2Ikeun1Bh/NRuoEK29+VD0+KamRUkqyymuvusoPPEQHOgs5l9TNZux0vMbEqwkaJuF7DL3Y9tPc4kq0Di2aXADnff5+6HgFlAs4Azxcp/zOw8gNDz5znxQ1UYP7UGqGpmVcysEFkHyOYGnCmqzMzI2q/9vruPCzpPLLj7ne5e3t0rk/XfeIm7x/Vfnu7+GfCpmVUPTWoLbA4wUix8AjQ1s2Khf+dtifMD/dnMBa4Nvb4WmJMTP7RATvyQeOHuh81sELCQrDMqJrn7poBjRVtz4Bpgo5mtC027y93/HWAmiY6bgKmhP4a2A9cFnCeq3H21mc0E3iXrbMC1xOFtQsxsOtAGKG1mu4C/AY8AqWZ2PVnDPPTIkXXp1iAiIhIO7ZISEZGwqDBERCQsKgwREQmLCkNERMKiwhARkbCoMEREJCwqDBERCYsKQyTGzKy8mV0ddA6RSKkwRGKvLfF/HyeJQ7rSWySGzKwFWff1+Rr4FrjC3bcHm0okPCoMkRgzswXAre7+3q8uLJKLaJeUSOxVB7YEHUIkUioMkRgys9Jkjfx2OOgsIpFSYYjEVmXicFxpyR9UGCKxtYWscQveM7P8MvqbxAkd9BYRkbBoC0NERMKiwhARkbCoMEREJCwqDBERCYsKQ0REwqLCEBGRsKgwREQkLP8HmLQueHUS3NwAAAAASUVORK5CYII=\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        ""
      ],
      "metadata": {
        "id": "Zp1Rs8homNZ9"
      },
      "id": "Zp1Rs8homNZ9",
      "execution_count": null,
      "outputs": []
    }
  ],
  "metadata": {
    "colab": {
      "machine_shape": "hm",
      "name": "CEE314Problemsets4Q3_teaching.ipynb",
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython"
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}