{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Content provided under a Creative Commons Attribution license CC-BY 4.0; code under BSD 3-Clause license. (c)2015 L.A. Barba, Olivier Mesnard, Pi-Yueh Chuang, Natalia Clementi."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Vortex-source panel method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In [Lesson 9](09_Lesson09_flowOverCylinder.ipynb) of _AeroPython_, you learned to use a _source panel method_ to represent a circular cylinder, and in [Lesson 10](10_Lesson10_sourcePanelMethod.ipynb) we used it for a symmetric airfoil at zero angle of attack. But what if we want the airfoil to generate some lift? If we place the airfoil at a non-zero angle of attack, we _should_ get lift, but will a source-panel representation be able to give you lift? Remember the [_Kutta-Joukowski theorem_](http://en.wikipedia.org/wiki/Kutta%E2%80%93Joukowski_theorem)?\n",
    "\n",
    "\n",
    "Historically, the first panel method ever developed was a source-sheet method. At the time, Douglas Aircraft Company was concerned with calculating the flow around bodies of revolution, and it was only later that the method was extended to lifting surfaces. (See the reference below for a nice historical account.)\n",
    "\n",
    "A *source-panel method* leads to a solution with no circulation, therefore no lift. The objective of this lesson is to start with the source panel method we implemented in the previous lesson and add some *circulation* so that we may have a lift force. We introduce an important concept: the **Kutta-condition** that allows us to determine what the right amount of circulation should be."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Reference\n",
    "\n",
    "* Smith, A.M.O., The Panel Method: Its Original Development. In _Applied Computational Aerodynamics_, Vol. 125, edited by P.A. Henne, published by AIAA (1990). [Read it on Google Books.](http://books.google.com/books?id=5Ov2tHj0wxoC&lpg=PA3&ots=SnUiqcdEnb&dq=The%20Panel%20Method%3A%20Its%20Original%20Development&pg=PA3#v=onepage&q&f=false)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A lifting-body panel method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we were to simply increase the angle of attack in the freestream and calculate the flow with a source sheet only, the rear stagnation point will not be located at the trailing edge. Instead, the flow will bend around the trailing edge and the stagnation point will be somewhere on the top surface of the airfoil. This is not a physically possible solution.\n",
    "\n",
    "For example, using the source-sheet panel method of [Lesson 10](10_Lesson10_sourcePanelMethod.ipynb) with an angle of attack $\\alpha=4^\\circ$ (using 40 panels), and plotting the streamlines in an area close to the trailing edge, we get the following plot:\n",
    "\n",
    "<center><img src=\"./resources/StreamlinesTrailingEdge.png\" width=\"600\"></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the streamlines behave strangely at the trailing edge. We know experimentally that the flow leaves the trailing edge of an airfoil smoothly, so this must be wrong. What's wrong?  We can't exclusively use sources to calculate the potential flow of an airfoil at non-zero angle of attack--we need circulation. But how do we obtain circulation?\n",
    "\n",
    "### The Kutta condition\n",
    "\n",
    "The *Kutta-condition* states that the pressure below and above the airfoil trailing edge must be equal so that the flow does not bend around it and instead leaves the airfoil tangentially. The rear stagnation point must be exactly at the trailing edge.\n",
    "\n",
    "It's natural to be a little perplexed by this. How can we justify this seemingly arbitrary condition? Remember that potential-flow theory completely ignores fluid viscosity, so if we are leaving out this physical effect, we shouldn't be surprised that the theory needs some adjustment for those situations when viscosity does play a role. A real viscous fluid is not able to turn around a sharp corner like an airfoil trailing edge without separating there. The Kutta condition allows us to correct potential-flow theory so that it gives a solution closer to reality.\n",
    "\n",
    "Remember [Lesson 6](06_Lesson06_vortexLift.ipynb), where we studied lift on a cylinder by combining a doublet and a freestream, plus a vortex. That's when we learned that **lift always requires circulation**. If you experimented with the circulation of the point vortex (which you *did*, right?), you found that the stagnation points moved along the cylinder.\n",
    "\n",
    "Like for the circular cylinder, the amount of circulation we add to an airfoil will move the stagnation points along the surface. And if we add just the right amount, the rear stagnation point can be made to coincide with the trailing edge. This amount of circulation makes the flow a physically relevant solution. And this amount gives the correct lift!\n",
    "\n",
    "To implement the Kutta-condition in our panel method we need to add one more equation to the system, giving the circulation that moves the stagnation point to the trailing edge. By placing a vortex-sheet with the same constant strength at every panel, we can add the circulation to the flow with just one more unknown.\n",
    "\n",
    "How do we enforce this in our code? We can re-use most of the code from [Lesson 10](10_Lesson10_sourcePanelMethod.ipynb), and enforce the Kutta-condition while adding circulation to the flow. Previously, we discretized the geometry into `N` panels, with a constant source strength on each one (varying from panel to panel), and applied a Neumann boundary condition of flow tangency at the `N` panel centers. This led to a linear system of `N` equations and `N` unknowns that we solved with the SciPy function `linalg.solve`. In the lifting-body case, we will instead have `N+1` equations and `N+1` unknowns. Read on to find out how!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Discretization into panels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's get the preliminaries out of the way. We need to import our favorite libraries, and the function `integrate` from SciPy, as in Lesson 10."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# import libraries and modules needed\n",
    "import os\n",
    "import numpy\n",
    "from scipy import integrate, linalg\n",
    "from matplotlib import pyplot\n",
    "\n",
    "# integrate plots into the notebook\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by importing the NACA0012 geometry from a data file, and we plot the airfoil:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# load geometry from data file\n",
    "naca_filepath = os.path.join('resources', 'naca0012.dat')\n",
    "with open(naca_filepath, 'r') as infile:\n",
    "    x, y = numpy.loadtxt(infile, dtype=float, unpack=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAngAAACTCAYAAAAdi2RSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl8nNV97/HPT7tla7EWa7XlHezE\ngMFg08S9Nlto8A2khYLbOEANJq8CDTcNMQm5SUNpWJIWmq0pLjTgCzUJ3FscbKDEsRNK7IKxSbjG\n2PGizdosybY2S9Zy7h8zmjuSZqSRLY1Gj77v1+u8nmXO88yZwzD+6jybOecQEREREe+IG+sGiIiI\niMjIUsATERER8RgFPBERERGPUcATERER8RgFPBERERGPUcATERER8ZiYC3hmdq2ZHTCzQ2b2QIjX\n/9DM9phZl5nd2O+1W83s9/5ya/RaLSIiIhI7LJbug2dm8cBB4GqgEngXWO2c+zCozkwgHfgysNk5\n95J/fRawG1gCOOA94BLn3IkofgQRERGRMRdrI3iXAYecc0ecc2eATcD1wRWcc6XOud8BPf22/RTw\npnOu0R/q3gSujUajRURERGJJrAW8IqAiaLnSv260txURERHxjISxbkA/FmJdpMeQI97WzNYB6wAm\nTZp0yfTp0yN8i9jR09NDXFys5XNvU59Hn/o8+tTn0ac+j77x3OcHDx6sd87lDlUv1gJeJRCctoqB\nqmFsu6LftjtCVXTOPQU8BbBkyRK3e/fu4bZzzO3YsYMVK1aMdTMmFPV59KnPo099Hn3q8+gbz31u\nZmWR1Iu1+PouMM/MZplZEnALsDnCbd8ArjGzqWY2FbjGv05ERERkQompgOec6wLuwRfM9gM/dc7t\nM7OHzOwzAGZ2qZlVAjcB/2xm+/zbNgJ/iy8kvgs85F8nIiIiMqHE2iFanHNbga391n0jaP5dfIdf\nQ237DPDMqDZQREREJMbF1AieiIiIiJw7BTwRERERj1HAExEREfEYBTwRERERj1HAExEREfEYBTwR\nERERj1HAExEREfEYBTwRERERj1HAExEREfEYBTwRERERj1HAExEREfEYBTwRERERj0kY6waIeFlP\nTw/t7e20trbS1tbG6dOn6ejoGLScOXOGzs5Ourq6hiwjzcxITEwMlISEhD7Lg60PVzcxMZGUlBQm\nTZoUmCYmJmJmI95+ERHxUcAT6aerq4sTJ05w4sQJGhsbOXHiBE1NTezevZs9e/bQ1NREc3Nzn2lT\nUxNtbW0DyunTp8f648SkuLg4Jk2aNKD0BsDe0tTUxPPPPx+yTmpqKpMnTyYtLY0pU6YMmE6ZMoX4\n+Pix/qgiImMi5gKemV0L/CMQD/yLc+7Rfq8nA88BlwANwM3OuVIzmwnsBw74q+5yzn0hWu2W2NTd\n3U1DQwPHjx+nrq6uz7R3vjfE9U6bm5tHtA3BYSQlJYXk5ORBS1JSEklJSSQkJAxa4uPjR3wUrLu7\nm66ursAIYmdnZ5/Sf10kdTo7O2lvb+f06dOB0tXVRWtrK62trSPa/v4mTZo0IPQNFgh759PT08nM\nzCQjIyMwTUlJ0aijiIwbMRXwzCwe+CFwNVAJvGtmm51zHwZVWwuccM7NNbNbgMeAm/2vHXbOXRTV\nRsuYcM5RX19PRUUFVVVVIUt1dTW1tbU454a1bzNj6tSpfUpGRgYtLS2cd955pKenk5aWRnp6emC+\nNxikpqYGwlxqaiopKSnExelU1/66uroGhL5QZc+ePcyaNSvka21tbbS2ttLS0kJzc3PIaW/d48eP\nn3Obk5KSyMjI6BP6+k8Hey09PZ2EhJj6yRURD4u1X5vLgEPOuSMAZrYJuB4IDnjXA3/jn38J+IHp\nz2rPcc5RXV3NoUOHKCsro6ysjPLy8j7zkR7+zM7OJjc3l2nTppGbmztgPjs7m6ysrECYS09PDxnK\nduzYwYoVK0b4k05MCQkJgVGzweTk5Jx1nzvnOH36dNjwFy4Y9h52P3XqFKdOneLkyZOcPHmSM2fO\nBEZ+z1ZaWhpZWVnDLikpKWf9niIyMcVawCsCKoKWK4Gl4eo457rM7BSQ7X9tlpntBZqArzvn3hrl\n9so5cM5RU1PD/v37OXToUJ9y+PBh2traBt0+MzOT6dOnU1RURGFhYcgybdo0EhMTo/SJJJaYWWBE\nNS8v75z3197eHgh8g03DvdZ7zmZzczNlZWXDeu/U1NQhQ2BOTk6g5ObmMnXqVJ2DKDKB2XAPX40m\nM7sJ+JRz7g7/8hrgMufcvUF19vnrVPqXD+Mb+WsBpjjnGszsEuDfgY8555pCvM86YB1AXl7eJZs2\nbRrlTzbyWlpahhz9iBW9h1OPHj0aGIErLS2lrKyMlpaWsNulp6dTVFREQUEB06ZNIy8vr0+ZPHly\nFD/F+Opzr/BSn/f09NDW1jbgAp3gUcNw07O5YtrMSEtLG3AIObj0P7SckpJCa2urZ/p8vPDS93y8\nGM99vnLlyvecc0uGqhdrI3iVwPSg5WKgKkydSjNLADKARudLqh0Azrn3/MFvPrC7/5s4554CngJY\nsmSJG4+H3WL1cGF3dzcHDhxg7969vP/++4FpQ0NDyPqZmZksWLCA+fPnM3fu3ECZM2cOU6dOjXLr\nBxerfe5l6nPfH0itra00NjaGLQ0NDTQ0NFBfX8/x48epr6+nsbExcIV3RUXF0G8EJCcnk5aWRnFx\ncWA0cNq0aYGSl5fXZz7af2R5lb7n0TcR+jzWAt67wDwzmwUcA24B/qxfnc3ArcBO4Ebgl845Z2a5\n+IJet5nNBuYBR6LX9Inp2LFj7Ny5k127drFr1y727NkT8ty4rKwsFi1axMKFC1m4cCELFixg4cKF\n5Ofn68pEkUGYWeB8xRkzZkS8XVdXF42NjdTX1/cJfoPNt7W10dHRQX19fUTvkZqaGjb89Z/Pzs7W\nIWORKIqpgOc/p+4e4A18t0l5xjm3z8weAnY75zYDTwMbzewQ0IgvBAL8IfCQmXUB3cAXnHON0f8U\n3uWcY9++fWzfvp233nqLnTt3UllZOaBeSUkJixcv5qKLLmLx4sUsXryY4uJiBTmRKEpISAiErEi1\ntbXx85//nHnz5gUuKOm9nVBdXR21tbV95tva2igtLaW0tHTIfcfFxfUZEczPz6egoCDkdOrUqfq9\nEDlHMRXwAJxzW4Gt/dZ9I2i+HbgpxHYvAy+PegMnEOccv//979m2bRvbt29nx44dA64gzMjIYOnS\npSxbtoxly5Zx2WWXkZ2dHWaPIhLLei9Iufjii4es65yjpaUlZPALNd/Q0BCYH0pSUlKfwBcuDObl\n5ZGUlDQSH13Ec2Iu4MnYOnPmDL/+9a959dVX2bJlC4cOHerzemFhIStXrmTFihX8wR/8Aeeff77u\n8yYyAfVexJGWlsacOXOGrN/Z2Ul9fT11dXXU1NRQW1tLdXU1NTU1A6ZNTU2Ul5dTXl4+5H6zs7Mj\nCoMZGRkaFZQJRQFPOH36NFu2bOHFF1/k9ddf73Nla1ZWFldddRVXXHEFK1euZN68efqRFJFhS0xM\npKCggIKCAi688MJB67a1tfUJfKFCYO+NzHsvMNm3b9+g+0xJSQkEwMLCwsDtlYqKigKlsLCQtLS0\nkfzYImNGAW+C6uzs5I033mDTpk288sorfULdokWLuO6661i1ahVLly7V3fdFJKpSU1OZPXs2s2fP\nHrRed3c39fX1YQNgcDhsaWmJ6HzBtLS0PoGvfwAsKioiPz9fv4sS8/QNnWDKysrYsGEDTz/9NDU1\nNYH1l156KatXr+aP//iPKSkpGcMWiohEJj4+PnBfzKFGBVtaWgJhr6qqimPHjgWmvaWqqorm5mY+\n+ugjPvroo7D7iouLIy8vL2wA7C06LCxjSQFvAnDO8ctf/pInnniCrVu3Bp7Nev7557NmzRpuvvnm\niM6hEREZr6ZMmRK4z2Y4zjlOnDgRNvz1ztfV1VFdXU11dTXvvfde2P2lpqYGQl9xcTHTp08PTHuL\nLkqT0aKA52HOObZu3crDDz/Mrl27AN/VaTfeeCN33XUXy5cv11+XIiJ+ZhZ49NuiRYvC1uvs7KSm\npiZsAOydb2lpCTx+MZyUlBSys7OZP39+IPT1D4GZmZn6rZZhU8DzqLfffpsvfvGLgb8us7Ozue++\n+7jrrrvIzc0d49aJiIxfiYmJgfA1mKampkDoq6yspKKiIjDtLadOnQrUCWfy5MkDQl//5fT09JH+\nmDLOKeB5TFVVFV/5yld4/vnnAcjPz+f+++9n3bp14/a5eyIi41F6ejrp6eksWLAgbJ3m5mZefvll\nCgsL+wS/4CDY0tLCgQMHOHDgwKDvVVxczIwZMygpKaGkpKTPfGFhoZ4kMsEo4HnIyy+/zNq1azl1\n6hTJycmsX7+e9evXk5qaOtZNExGRENLS0pg5c2bY56I65zh16lTI0b/gdU1NTXz44Yd8+OGHIfeT\nkJBAcXFxIPD1D4AzZswgJSVlFD+pRJsCngd0dHTwpS99iR/96EcAXHfddXz/+99n1qxZY9wyERE5\nF2ZGZmYmmZmZYc8LdM7R2NhIRUUF5eXllJWVBaa9pba2dsjbxOTl5YUNgCUlJWRmZo7Sp5TRoIA3\nzrW1tXHDDTfw5ptvkpSUxHe/+13uuecenZArIjJBmBnZ2dlkZ2dz0UUXhazT3t7eJ/T1D4CVlZXU\n1tZSW1vLO++8E3If6enpYcPfzJkzycvL0789MUQBbxxra2tj1apVbN++nWnTprFlyxaWLFky1s0S\nEZEYk5KSwvz585k/f37I17u7u6mqqgoZ/npLU1MTH3zwAR988EHIfUyaNImZM2cya9askEUjgNGl\ngDeO3XfffWzfvp2CggK2bds26Im8IiIi4cTHxw96ZbBzjoaGhrDhr7S0lMbGRvbv38/+/ftD7iMz\nMzNs+Js5cyaTJk0azY844SjgjVNvvfUWGzZsIDk5mddff13hTkRERo2ZkZOTQ05ODhdffHHIOqdO\nnaK0tJSjR48GypEjRzh69CilpaWcPHmSvXv3snfv3pDb5+fnhwx/c+bMobi4WFcBD5MC3jjU09PD\nhg0bAHjssce44IILxrhFIiIy0WVkZHDhhReGfGycc466uro+4S+4lJeXB54dvHPnzgHbJyYmMmvW\nLGbPns2cOXP6lNmzZ2v0L4SIA56ZvQ38GPipc65jtBpkZtcC/wjEA//inHu03+vJwHPAJUADcLNz\nrtT/2leBtUA38FfOuTdGq51jadu2bVRUVFBcXMzdd9891s0REREZlJkFnhu8bNmyAa93dXVx7Nix\nkOHv8OHD1NTUcPDgQQ4ePBhy/4WFhX0CX3AAzM7OnpAXfwxnBK8TeBZ40syeBZ5yzoV/GvNZMLN4\n4IfA1UAl8K6ZbXbOBd/YZy1wwjk318xuAR4DbjazhcAtwMeAQuAXZjbfOdc9km2MBVu2bAHg9ttv\nJyFBg7AiIjK+JSQkBK7IDXVPwNbWVo4cOcKRI0c4fPhwn1JaWkpVVRVVVVW89dZbA7ZNT08fMOKX\nlpYWhU81tiJOB865FWZ2HnAX8Hngi2b2FvBPwP92znWOQHsuAw45544AmNkm4HogOOBdD/yNf/4l\n4Afmi+bXA5v8o4tHzeyQf38Dx3rHud7nGoa7HF5ERMRLJk+ezKJFi0LeC7Crq4uKiooBwa83DDY1\nNYU89++zn/2sp2/uPKzhH+fcAeBL/kOhfwqsA14A6s3sX/GN6h05h/YUARVBy5XA0nB1nHNdZnYK\nyPav39Vv26JQb2Jm6/xtJy8vjx07dpxDk6Ov90aV5eXl467t41lLS4v6O8rU59GnPo8+9fnISEhI\n4LzzzmP+/PmB5wBXVVWxf/9+du/eTXl5eZ/6jz76aNgniHjBWR3f84+SbTSzfcA/AH8IfAX4spn9\nH+Be51zNWew61EFyF2GdSLb1rXTuKeApgCVLlrjx9h/4ggsuYN++fWRmZnr6yxlrduzYof6OMvV5\n9KnPo099Pnzd3d1UVlYOGLXrLU1NTUPu4+tf/7qnT3Ma9iczs0nAauAL+C50+Aj4IvAz4L/jO3z6\nPHDlWbSnEgi+CU8xUBWmTqWZJQAZQGOE23rCxz72MQD+4z/+g9tuu21sGyMiIjIKWltbKS0tDRng\nSktL6ewMf2bYlClTBlxt21umT5/Of/7nf3o63MHwrqJdhO/8uz8HJgOvAOudc9uDqm0wsxp8Ye9s\nvAvMM7NZwDF8F038Wb86m4Fb8Z1bdyPwS+ecM7PNwAtm9g/4LrKYB4R+3so4t2bNGr7xjW/w0ksv\n8Z3vfIeiopBHokVERGJWR0cH5eXlgfvkBV85W1paSl1d3aDb5+fnhw1xOTk5E/LK2WDDia+/xTci\n9iS+c+2qw9Q7xFle2OA/p+4e4A18t0l5xjm3z8weAnY75zYDT+M7PHwI38jdLf5t95nZT/FdkNEF\n3O3FK2gBZsyYwfLly/nVr37FnXfeyZYtWyb8F1lERGJL/1uf9A9xVVVVOBfyTCrAd++7kpKSkLc/\nmT17NpMnT47ipxl/hhPwbgL+fajQ5JzbD6w82wY557YCW/ut+0bQfLu/LaG2/Tvg7872vceTu+++\nmw8++IDXXnuNRx55hK997Wtj3SQREZlAenp6qKmp6RPcgucrKiro6uoKu31cXBwzZswIPKos+LFl\ns2bNorCwkLi4uCh+Im8Zzm1SXh7Nhsjw5Obm8vTTT/PZz36WBx98kLi4OB544IGxbpaIiHhET08P\ntbW1lJWVUV5ePiDIlZaW0tEx+HMPCgoKBgS33vnp06eTmJgYpU8z8Xj7DEOPu+GGG3j66ae54447\n+OpXv8qxY8f47ne/S3Jy8lg3TUREYlxbWxsVFRWUl5cHQlzwfEVFxaAXMgDk5OSEHYErKSnx9H3m\nYp0C3jj3F3/xFyQmJrJ27Vp+8IMf8Pbbb/Piiy8yb968sW6aiIiMkd5nv/YPbcHL9fX1Q+4nNzeX\nGTNmMGPGDEpKSvqEuJkzZ06IJ0KMVwp4HrBmzRoWLFjAzTffzN69e7ngggtYv34969ev1wOYRUQ8\nqL29ncrKypAjb71lqMOniYmJTJ8+nZKSkj4hrnd++vTppKamRukTyUhTwPOIJUuWsGfPHu699142\nbtzIt771LX7yk5/w6KOPctNNNxEfHz/WTRQRkQh0dnZSVVVFZWVln9J7OLW8vJza2toh95OVlRUy\nuPUu5+Xl6SIGD1PA85CMjAyee+457rzzTu69915++9vfsnr1ar71rW/xta99jdWrV3v+xo4iIrGs\no6ODY8eO9Qluu3bt4nvf+15guaamZtDbhwDEx8dTXFw86OibDp9ObPrX3oOWL1/O7t27eeaZZ/j2\nt7/NRx99xOc//3m++c1v8pd/+Zfcdttt5OTkjHUzRUQ85fTp0wNG3fqXoW7eC2BmFBYWUlxcHCjT\np0+nqKgoEOQKCgp0ZEYGpYDnUQkJCaxbt47bb7+d559/nkceeYSDBw9y//338+CDD3LjjTdy1113\nsXz5ct0kWURkED09PRw/fpyqqqoBJXg0rqGhYch9xcfHDwhv7e3trFy5MrCcn5+v24fIOVPA87jE\nxERuu+021qxZw2uvvcaPf/xjtm7dygsvvMALL7xASUkJN998M7fccgsXXXSRwp6ITBjOORobG0MG\nt+BSXV1Nd/fQD0ZKTEykqKioT3jrPwqXl5c3YORtx44drFixYpQ+pUxUCngTRHx8PKtWrWLVqlWU\nlZWxYcMGnn32WcrKynj88cd5/PHHmT9/Pn/yJ3/CqlWrWLp0qYb/RWRccs7R1NQ0ZHCrqqrizJkz\nEe0zOzubwsLCkKU3wE2bNk0XLUjMUMCbgEpKSnj44Yd56KGH+M1vfsOmTZv42c9+xsGDB3nkkUd4\n5JFHyM7O5tOf/jTXXXcdV199NVlZWWPdbBGZ4Lq6ujh+/Dg1NTXU1tb2mdbU1PQJbm1tbRHtMyMj\nI2xw6y35+fm6Ya+MOwp4E1hcXByf/OQn+eQnP8mTTz7Jjh07+PnPf86rr77KkSNH2LhxIxs3bsTM\nuPDCC1m5ciVXXHEFy5cvJyMjY6ybLyIe0NPTQ0NDQ5+gFm6+vr5+yKtLe6WmplJUVDRocCsoKNAD\n68WzFPAE8F2UcdVVV3HVVVfx5JNPcuDAAV599VW2bNnCb37zG95//33ef/99nnjiCeLi4rj44ou5\n/PLLufzyy1m2bBkzZ87U+XsiAvgOkZ48eTKi0FZXVxfR+W3gu7o0NzeX/Px88vLyyM/PD8zn5eX1\nCXRpaWn6TZIJTQFPBjAzzj//fM4//3y+/OUv09bWxs6dO9m+fTvbt2/nnXfeYffu3ezevZvvf//7\nAOTl5bFs2TKWLVvG4sWLueiii8jLyxvjTyIiI6W9vZ3jx49TV1fH8ePHQ87X1tYGwluk57YBTJ06\nNWRo6z+fm5ure3mKRChm/k8xsyzgRWAmUAr8qXPuRIh6twJf9y8+7Jx71r9+B1AAnPa/do1zbugb\nDsmQUlNTufLKK7nyyisBaGlpYdeuXX1KbW0tr7zyCq+88kpgu4KCgkDYW7x4MR//+MeZM2eOLv8X\niQHhAtt7773Hxo0bB6xvaWkZ1v7T0tIiCm3Tpk0jOTl5lD6lyMQVMwEPeADY5px71Mwe8C+vD67g\nD4HfBJYADnjPzDYHBcE/d87tjmajJ6IpU6YEDueC73DM4cOH2bVrF++8807gcG51dTXV1dVs3bo1\nsG1iYiLz589n4cKFgbJgwQLmzp2r5+aKnKW2tjYaGhqor6+noaEhUMItHz9+fNiBLTExkdzcXHJz\nc5k2bVrI+d7QlpeXp2eYioyxWAp41wMr/PPPAjvoF/CATwFvOucaAczsTeBa4N+i00QJxcyYO3cu\nc+fO5XOf+xzgO3H6yJEj7N27l7179/L++++zf/9+SktL2bdvH/v27Ruwn+Li4sB+esucOXOYNWuW\nLuqQCaGzs5OTJ09y4sQJTpw40SecDRbY2tvbh/1e4QJbS0sLS5cuHbA+IyND57SJjCOxFPDynHPV\nAM65ajObFqJOEVARtFzpX9frX82sG3gZ3+HbyC63khEXFxcXCGk33XRTYH1LSwsHDhzgww8/DJT9\n+/dz9OjRwN3gd+zYMWB/6enplJSU9Cm9z1wsKioiPz+fpKSkKH5CkdDa29sDAS04rEWybrijar2S\nk5PJzs4mJyeH7OzsQAm3nJOTEzaw6aa7It5g0cxAZvYLID/ESw8CzzrnMoPqnnDOTe23/f1AsnPu\nYf/y/wTanHN/b2ZFzrljZpaGL+D9L+fcc2HasQ5YB5CXl3fJpk2bRuLjRVVLSwtTpkwZ62aMmO7u\nbmprazl27Fig9D4GqKamho6OjiH3kZmZSVZWVp9/xLKyspg6dSoZGRlMnTqVzMxM0tLSzupmpF7r\n8/Egmn3e09PD6dOnaWtro62tjdbW1pDT/utaWlpoaWmhubmZlpaWYV1c0F9cXBxTpkwJlPT0dDIy\nMkJOe0tGRgYpKSkjNrqm73n0qc+jbzz3+cqVK99zzi0Zql5UR/Ccc1eFe83Mas2swD96VwCEukCi\nkv9/GBegGN+hXJxzx/zTZjN7AbgMCBnwnHNPAU8BLFmyxI3Hv1Yn0l/ZzjkaGhooKyujrKyM8vLy\nPvPV1dXU1NRw8uRJTp48yZEjRwbdX1xcHDk5OYFDUFlZWYEg2DvtP5+RkcGePXsmTJ/HinDf856e\nnkDACg5b4ZZbW1tpbm6mqakpbGlubh6RNicmJga+N71/VAy2HLzubP/4GEkT6bclVqjPo28i9Hks\nHaLdDNwKPOqfvhKizhvAt82sd2TvGuCrZpYAZDrn6s0sEVgF/CIKbZYoMDNycnLIycnhkksuCVmn\nu7s75MPAq6ur+1wpWFdXx8mTJ6mrq6Ouri7kuYCDSUlJIT09nbS0tMAISu/85MmTSU1NHVCC16ek\npJCcnDxkGet/5IfinKOrq4uOjg7OnDlDR0fHsMpg2wQHtGPHjpGYmDggsJ3NOWeRmDx5cp/RseBR\nslDre0twUEtNTdW5aiIy5mIp4D0K/NTM1gLlwE0AZrYE+IJz7g7nXKOZ/S3wrn+bh/zrJgNv+MNd\nPL5wtyH6H0HGSnx8fOD2CxdffPGgdc+cOUN9fX0g8DU2NgbOhQo33zvK097eTnt7O3V1o3sHHjMj\nPj6euLg44uPj+8z3n0ayDnwhuKen56yn/deNtd7w3Bugh5oPDuWhSlpamp6/LCKeETMBzznXAFwZ\nYv1u4I6g5WeAZ/rVaQVCD+2I9JOUlBS42/1wbN++naVLlwYO5wVPm5qaAudn9T9fK3jd6dOnIxrd\n6h0hi2VxcXERjUYOtwSHsoMHD/KJT3xiQFhLSUmJ+VFOEZGxFDMBTyTWmVngUGt+fqhrhUaGc67P\niFmoEbXBRttCrTOzsCOAkUzDrRttqampXHrppaP+PiIiXqOAJxJjeg/PxsfH66kfIiJyVnSMQ0RE\nRMRjFPBEREREPEYBT0RERMRjFPBEREREPEYBT0RERMRjFPBEREREPEYBT0RERMRjFPBEREREPEYB\nT0RERMRjFPBEREREPEYBT0RERMRjFPBEREREPCZmAp6ZZZnZm2b2e/90aph6r5vZSTN7td/6WWb2\nX/7tXzSzpOi0XERERCS2xEzAAx4Atjnn5gHb/MuhfAdYE2L9Y8AT/u1PAGtHpZUiIiIiMS6WAt71\nwLP++WeBG0JVcs5tA5qD15mZAVcALw21vYiIiIjXxVLAy3POVQP4p9OGsW02cNI51+VfrgSKRrh9\nIiIiIuNCQjTfzMx+AeSHeOnBc911iHVukHasA9b5F1vM7MA5vv9YyAHqx7oRE4z6PPrU59GnPo8+\n9Xn0jec+L4mkUlQDnnPuqnCvmVmtmRU456rNrACoG8au64FMM0vwj+IVA1WDtOMp4Klh7D/mmNlu\n59ySsW7HRKI+jz71efSpz6NPfR59E6HPY+kQ7WbgVv/8rcArkW7onHPAduDGs9leRERExEtiKeA9\nClxtZr8HrvYvY2ZLzOxfeiuZ2VvAz4ArzazSzD7lf2k98CUzO4TvnLyno9p6ERERkRgR1UO0g3HO\nNQBXhli/G7gjaHl5mO2PAJcXR69RAAAFRElEQVSNWgNjz7g+xDxOqc+jT30eferz6FOfR5/n+9x8\nRzdFRERExCti6RCtiIiIiIwABbwYZ2bXmtkBMztkZgOe7mFmyf5Hsx3yP6ptZvRb6S0R9PmXzOxD\nM/udmW0zs4guWZfwhurzoHo3mpkzM09f/RYNkfS5mf2p/7u+z8xeiHYbvSaC35YZZrbdzPb6f18+\nPRbt9Aoze8bM6szs/4Z53czse/7/Hr8zs4uj3cbRpIAXw8wsHvgh8EfAQmC1mS3sV20tcMI5Nxd4\nAt8j2+QsRdjne4ElzrkL8D095fHottJbIuxzzCwN+Cvgv6LbQu+JpM/NbB7wVeATzrmPAfdFvaEe\nEuH3/OvAT51zi4FbgB9Ft5We8xPg2kFe/yNgnr+sA/4pCm2KGgW82HYZcMg5d8Q5dwbYhO+RbsGC\nH/H2Er6ri0Pd+FkiM2SfO+e2O+fa/Iu78N13Uc5eJN9zgL/FF6bbo9k4j4qkz+8EfuicOwHgnBvO\nvUlloEj63AHp/vkMBrmfqwzNOfdroHGQKtcDzzmfXfjup1sQndaNPgW82FYEVAQth3oEW6CO/ybP\np/DdJkbOTiR9Hmwt8Nqotsj7huxzM1sMTHfOvRrNhnlYJN/z+cB8M3vbzHaZ2WAjITK0SPr8b4DP\nmVklsBW4NzpNm7CG+3s/rsTMbVIkpEgewTasx7TJkCLuTzP7HLAE+G+j2iLvG7TPzSwO3+kHt0Wr\nQRNAJN/zBHyHrlbgG6V+y8w+7pw7Ocpt86pI+nw18BPn3N+b2eXARn+f94x+8yYkT//7qRG82FYJ\nTA9aDvUItkAdM0vAN6w/2JC0DC6SPsfMrsL3DOXPOOc6otQ2rxqqz9OAjwM7zKwUWAZs1oUW5yTS\n35ZXnHOdzrmjwAF8gU/OTiR9vhb4KYBzbieQgu+ZqTI6Ivq9H68U8GLbu8A8M5tlZkn4Trrd3K9O\n8CPebgR+6XRzw3MxZJ/7Dxf+M75wp/OSzt2gfe6cO+Wcy3HOzXTOzcR33uNn/DdBl7MTyW/LvwMr\nAcwsB98h2yNRbaW3RNLn5fhv+G9mC/AFvONRbeXEshn4vP9q2mXAKedc9Vg3aqToEG0Mc851mdk9\nwBtAPPCMc26fmT0E7HbObcb3SLaN/ke0NeL70ZCzFGGffweYAvzMfz1LuXPuM2PW6HEuwj6XERRh\nn78BXGNmHwLdwP3+Jw7JWYiwz/8a2GBm/wPfocLb9Af72TOzf8N3ikGO/7zGbwKJAM65H+M7z/HT\nwCGgDbh9bFo6OvQkCxERERGP0SFaEREREY9RwBMRERHxGAU8EREREY9RwBMRERHxGAU8EREREY9R\nwBMRERHxGAU8EREREY9RwBMRERHxGAU8EZFzZGaTzewjM3vHzBKD1l9jZj1mdvdYtk9EJh49yUJE\nZAT4n1G8C3jCOfeAmU0Dfge8o0fZiUi0KeCJiIwQ/zNE/x64BvgysAi40DlXP6YNE5EJRwFPRGSE\nmJkBW4ArgCTgaufctrFtlYhMRDoHT0RkhDjfX8wbgWTgtwp3IjJWFPBEREaImeUDTwJ7gAvN7Itj\n3CQRmaAU8ERERoD/8OyzwBnganxB7zEzu2BMGyYiE5LOwRMRGQFm9tfA48AVzrlfmVkSvqtqk4El\nzrnTY9pAEZlQNIInInKO/LdI+TbwiHPuVwDOuTPAamAm8A9j1zoRmYg0giciIiLiMRrBExEREfEY\nBTwRERERj1HAExEREfEYBTwRERERj1HAExEREfEYBTwRERERj1HAExEREfEYBTwRERERj1HAExER\nEfGY/wcGO1ReCcYaXgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1133052b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot geometry\n",
    "width = 10\n",
    "pyplot.figure(figsize=(width, width))\n",
    "pyplot.grid()\n",
    "pyplot.xlabel('x', fontsize=16)\n",
    "pyplot.ylabel('y', fontsize=16)\n",
    "pyplot.plot(x, y, color='k', linestyle='-', linewidth=2)\n",
    "pyplot.axis('scaled', adjustable='box')\n",
    "pyplot.xlim(-0.1, 1.1)\n",
    "pyplot.ylim(-0.1, 0.1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The contour defining the airfoil will be partitioned into `N` panels, using the same method as in [Lesson 10](10_Lesson10_sourcePanelMethod.ipynb).\n",
    "\n",
    "We define a class `Panel`  that will store all information about one panel: start and end points, center point, length, orientation, source strength, tangential velocity and pressure coefficient. We don't save the vortex-sheet strength because all panels will have the same value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Panel:\n",
    "    \"\"\"\n",
    "    Contains information related to a panel.\n",
    "    \"\"\"\n",
    "    def __init__(self, xa, ya, xb, yb):\n",
    "        \"\"\"\n",
    "        Initializes the panel.\n",
    "        \n",
    "        Sets the end-points and calculates the center-point, length,\n",
    "        and angle (with the x-axis) of the panel.\n",
    "        Defines if the panel is located on the upper or lower surface of the geometry.\n",
    "        Initializes the source-strength, tangential velocity, and pressure coefficient\n",
    "        of the panel to zero.\n",
    "        \n",
    "        Parameters\n",
    "        ---------_\n",
    "        xa: float\n",
    "            x-coordinate of the first end-point.\n",
    "        ya: float\n",
    "            y-coordinate of the first end-point.\n",
    "        xb: float\n",
    "            x-coordinate of the second end-point.\n",
    "        yb: float\n",
    "            y-coordinate of the second end-point.\n",
    "        \"\"\"\n",
    "        self.xa, self.ya = xa, ya  # panel starting-point\n",
    "        self.xb, self.yb = xb, yb  # panel ending-point\n",
    "        \n",
    "        self.xc, self.yc = (xa + xb) / 2, (ya + yb) / 2  # panel center\n",
    "        self.length = numpy.sqrt((xb - xa)**2 + (yb - ya)**2)  # panel length\n",
    "        \n",
    "        # orientation of panel (angle between x-axis and panel's normal)\n",
    "        if xb - xa <= 0.0:\n",
    "            self.beta = numpy.arccos((yb - ya) / self.length)\n",
    "        elif xb - xa > 0.0:\n",
    "            self.beta = numpy.pi + numpy.arccos(-(yb - ya) / self.length)\n",
    "        \n",
    "        # panel location\n",
    "        if self.beta <= numpy.pi:\n",
    "            self.loc = 'upper'  # upper surface\n",
    "        else:\n",
    "            self.loc = 'lower'  # lower surface\n",
    "        \n",
    "        self.sigma = 0.0  # source strength\n",
    "        self.vt = 0.0  # tangential velocity\n",
    "        self.cp = 0.0  # pressure coefficient"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Like before, we call the function `define_panels` to discretize the airfoil geometry in `N` panels. The function will return a NumPy array of `N` objects of the type `Panel`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def define_panels(x, y, N=40):\n",
    "    \"\"\"\n",
    "    Discretizes the geometry into panels using 'cosine' method.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    x: 1D array of floats\n",
    "        x-coordinate of the points defining the geometry.\n",
    "    y: 1D array of floats\n",
    "        y-coordinate of the points defining the geometry.\n",
    "    N: integer, optional\n",
    "        Number of panels;\n",
    "        default: 40.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    panels: 1D Numpy array of Panel objects.\n",
    "        The list of panels.\n",
    "    \"\"\"\n",
    "    \n",
    "    R = (x.max() - x.min()) / 2.0  # circle radius\n",
    "    x_center = (x.max() + x.min()) / 2.0  # x-coordinate of circle center\n",
    "    \n",
    "    theta = numpy.linspace(0.0, 2.0 * numpy.pi, N + 1)  # array of angles\n",
    "    x_circle = x_center + R * numpy.cos(theta)  # x-coordinates of circle\n",
    "    \n",
    "    x_ends = numpy.copy(x_circle)  # x-coordinate of panels end-points\n",
    "    y_ends = numpy.empty_like(x_ends)  # y-coordinate of panels end-points\n",
    "    \n",
    "    # extend coordinates to consider closed surface\n",
    "    x, y = numpy.append(x, x[0]), numpy.append(y, y[0])\n",
    "    \n",
    "    # compute y-coordinate of end-points by projection\n",
    "    I = 0\n",
    "    for i in range(N):\n",
    "        while I < len(x) - 1:\n",
    "            if (x[I] <= x_ends[i] <= x[I + 1]) or (x[I + 1] <= x_ends[i] <= x[I]):\n",
    "                break\n",
    "            else:\n",
    "                I += 1\n",
    "        a = (y[I + 1] - y[I]) / (x[I + 1] - x[I])\n",
    "        b = y[I + 1] - a * x[I + 1]\n",
    "        y_ends[i] = a * x_ends[i] + b\n",
    "    y_ends[N] = y_ends[0]\n",
    "    \n",
    "    # create panels\n",
    "    panels = numpy.empty(N, dtype=object)\n",
    "    for i in range(N):\n",
    "        panels[i] = Panel(x_ends[i], y_ends[i], x_ends[i + 1], y_ends[i + 1])\n",
    "    \n",
    "    return panels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use our new function to define the geometry for the airfoil panels, and then plot the panel nodes on the geometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# discretize geoemetry into panels\n",
    "panels = define_panels(x, y, N=40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAngAAACTCAYAAAAdi2RSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl841W9//HXJ1vTpvu+z3RWFhGE\nAVFQB1lEvcK9/uQK1wUVQdwVQdyvIiLX3euCouIFBUcEr46Il4cgIyoMMIAwDDDQ6XTf23RJ0zTb\n5/dHMqXTpp3OlqaZz/PxyCPJN99vcnImk75zzvecI6qKMcYYY4zJHo6lLoAxxhhjjDm0LOAZY4wx\nxmQZC3jGGGOMMVnGAp4xxhhjTJaxgGeMMcYYk2Us4BljjDHGZJmMC3gicq6I7BSRZhH5VIrHXy0i\nj4tIVETeMuuxi0XkheTl4vSV2hhjjDEmc0gmzYMnIk7geeBsoBN4FLhIVZ+Zsc9KoBC4Etisqnck\nt5cC24ANgAKPASepqj+Nb8EYY4wxZsllWgveKUCzqraoahjYBJw/cwdVbVXVp4D4rGNfB/xZVYeT\noe7PwLnpKLQxxhhjTCbJtIBXB3TMuN+Z3Ha4jzXGGGOMyRqupS7ALJJi22L7kBd9rIhcBlwGkJub\ne1JDQ8MiXyJzxONxHI5My+fZzeo8/azO08/qPP2sztNvOdf5888/P6iqFfvaL9MCXicwM23VA937\ncezGWcduSbWjqt4I3AiwYcMG3bZt2/6Wc8lt2bKFjRs3LnUxjihW5+lndZ5+VufpZ3Wefsu5zkWk\nbTH7ZVp8fRRYKyJNIuIBLgQ2L/LYe4BzRKREREqAc5LbjDHGGGOOKBkV8FQ1CnyIRDB7FrhdVXeI\nyDUich6AiJwsIp3ABcCPRWRH8thh4MskQuKjwDXJbcYYY4wxR5RM66JFVe8G7p617Qszbj9Kovs1\n1bE3ATcd1gIaY4wxxmS4jGrBM8YYY4wxB88CnjHGGGNMlrGAZ4wxxhiTZSzgGWOMMcZkGQt4xhhj\njDFZxgKeMcYYY0yWsYBnjDHGGJNlLOAZY4wxxmQZC3jGGGOMMVnGAp4xxhhjTJaxgGfMMjG0+Vds\nf/UaHl+bw/ZXr2Fo86+WukjGGGMylAU8Y/bDUoWsoc2/ov0z7yfS3Q6qRLrbaf/M+5ck5FnQNMaY\nzOda6gIYs1y4tt5H+y+/h4aCANMhC6DsvItSHhOPxwmFQkxMTBAMBpmcnGRqamrBSzgcJhKJEI1G\npy+n3/YNcpOvu4eGguz47Af53788nNggctDvUURwu93TF5fLtdf9qhf+SdM9m3BEw9N1sPvqS3ny\niX8SfcWZe+078+L1esnNzZ2+drvdyCEo79DmX9H9jc8T6enAXdNA7ZVfnvffwhhjjiQW8IyZJRqN\n4vf78fv9DA8P4/f7GRsbo/bWH06Huz00FOSpqy7jB1/5OrHQJLHQJPHQJPGpEISnIBImR8AjQo6D\n5O2Z1zJ92yPgdSSu82dszxHwOkkZ4HzBcd5x1w/nbI+rooDuKeeei764ba/te/bVWfdnXRc6wDGr\nHI5ImLybv8XTN36TKYUphXDyMhXX6ft7XSPgcoPHA+4cJCcHyfEiHi9Oby7O3FxcuXn4J4LcfnMT\nbp8Pb56P3Ly86ZBYt/tpVt7zKxyRF8Nm26ffh394mIp/ezv5+fk4nc5F/qsbY0x2ybiAJyLnAt8F\nnMBPVfX6WY/nALcAJwFDwFtVtVVEVgLPAjuTu25V1cvTVW6zNPbVghOLxRgaGmJgYID+/v69rgcG\nBujv62NiaICof4j46DASGCNnapJip1DkhGIHFDuFYies8UrKkFUYDfHWnu3TgSasMCUw5YGw20HM\n4STucqEuD3GXG3W5iLvc4PIw5fEw5c5BPB7E450OOo4cL7hzwONBPDmcff8m8ibH57z2VEEJD1/8\n6cQdnRHFkmlN9qS25OOy134ztqGgSjwWS7Qa7mlBjEQSrYnJy39s+WXKfwe3Q3ik6QQkGsYRjSLR\nCI5oBEcsiiMWxRkJ44tGKIhGcMaiuDROjiOCJx4hJzxBTgQ8E8mQK5IIvHsCceuD5AgILwbEKYUy\n59ywyVSInms+xvev/AjBOEw53cQ8Oag3D7x5kJePM78AZ0ERnqJicoqKyS8sIj8/n/z8fAoKCsjP\nz6ewsJDi4mKKioqmr71e76JaHa1V0RiTCTIq4ImIE/gBcDbQCTwqIptV9ZkZu10C+FV1jYhcCPwX\n8NbkY7tU9YS0FtosCVWl7dYbGb7uKgiHgEQLTvOV7+GH11/H4xNRwkP96NgIxQ4ockCx88Wwdowj\ncV3khIjCSAxGY8qIB0acDkLuHMLePKby8unNL6KvsIR1ux+lMBaeU5ZYaSWV37iN/Px88vLyyMvL\nw+fzkZdsbXI4Dv5U16HNp9P+mffv1YIo3jzWX/NdXpnG8LD91Q8kzgOcxVPbyI0PPLro54lGo4RC\nISYnJ+e9BCcn+fvjj9PU1MTk5CShiQBTgQDhiQCRiQAX3J86bOaJMOHKwR0JU6QRfJEoebEJ8ifB\nNwo+h5DngHyBXEciLE7E91yUoTi0x2FCE/f3PBZyOFFvHpLrw5FfiKugEE9xCd7iUnJLy8kvLWNN\nXwvr7r9zTqviRCBA7b+/G5cro75yjTFZLNO+bU4BmlW1BUBENgHnAzMD3vnAF5O37wC+L4fiZB6T\nUVSVnp4empubaWtro62tjfa2NgZ3NxPu2I1zsJcrCmMUOPf+p3fHY7yu7zlWh2EkrozkCiGPl1iu\nDy0oJlBcSriskkB1DaO19RTXr6C0uoa60lJeUlJCSUkJhYWFKUPZ36//PDLjHDxIhKzVn/8GZa95\nzWGtjz0tQEvdMlR75ZdTBs3aK7+8X8/jcrmmW80WUl5ezsaNG1M+Nm/YrGvkJw80o6pMTk4yPj5O\nIBCYc90XCDA+NsrkiJ+Qf5jIqJ/I2AjR8TFigXHiE+MQDOAMBSmaClGtMXzRAL6JAL7JPnxDgk/A\n50hcXJIYtZaqVbH1sx/gix+8nCl3DnFfAVJQhLOolJyycrzlVRRWVFJaVkZpaWnKi9fr3a/6ncla\nFI05MmVawKsDOmbc7wRePt8+qhoVkVGgLPlYk4g8AYwBn1PVvx3m8ppFmO8PjKrS29vLs88+S3Nz\nM83Nzex64QWGXniOWFcrFbEp6t1CvVtY5RZe7Ya4QkdE6XRB/jwtY05xsPa3f6e2tpbKykrcbvch\neR/RU8+k8ZhjluyPZdl5Fy35H+ZMCZqw77ApItMtqlVVVQf9eqFQiNHRUUZGRva6bk1ejw4P8a+3\nfzPlsT6H4HMItYQpnByiKDxEoX83Re2JVmSXwGgMRuMwFFNa4jAWg7G4MhqDSZebeN6eYFiCp7Qc\nb3klBRVV08GwvLx8+lJRUUFJSQkjf7x9rzpazMAgY0x2EJ15Ps4SE5ELgNep6nuT998BnKKqH56x\nz47kPp3J+7tItPwFgHxVHRKRk4DfAceq6liK17kMuAygqqrqpE2bNh3md3boBQKBfbZ+ZALX1vvI\nveU7SHhqeltYHNzpKOWfgyOJEOdiOsjVuBJ/5DojSh9uxvOLCJdUoBW1uBuaKGlYQVVVVeLypffh\nGO6f85rx0koCX0vdfXcwlkudZ5N91blr6314f/tzZHgALa0g9OZ3Ez31zDSWcG/5n3z7vJ/Jsetv\nIRgMMj4+ztjY2F7XwZERwiNDxEf96PgoMjGOc3IC99QkOZEQ+cQpcgqFDihyCkUOKHQmzlEcjcFY\nPHGKQeIaRuPKeBzeVuykMMXvoHFvPv/3pvdNn1+451xDr9fLxMSEfc7TzL5b0m851/kZZ5zxmKpu\n2Nd+mRbwXgF8UVVfl7z/aQBV/eqMfe5J7vOQiLiAXqBCZ70REdkCXKmq2xZ6zQ0bNui2bQvukpG2\nbNkyb9fVUorFYuzcuZMnnniCfz7xBG/4w48oSnHeWlSVRyaVAYeHaHk1noZVFK07hprjT2L1Mcey\nevVqSkpKFnytPXPDzW7BabzuhsPSOpGpdZ7NlludH47PpKoyMTHB8PDwnIt/oJ9AXw+T/b1MDQ0Q\n8Q8RG/UjE+O4poJcUuxIOTBEVemIgD+u+GMwHAN/TBkXF5MeL+6ySpyl5eRUVJNfU09FdTWVlZVU\nVVVRWVk5fdvn8x1QHWVCC3AmWW6f82ywnOtcRBYV8DKti/ZRYK2INAFdwIXAf8zaZzNwMfAQ8Bbg\nL6qqIlIBDKtqTERWAWuBlvQV/cjU1dXFQw89xNatW3n8wX8w9cwTHCURXuoVXu9NtDikGnnqFAf/\n77FWqqurD3g+tEzqLjQGDs9nUkSmz1dsbGxc9HHRaJQdr1lLrK9rzmPB3HzuPeHVhAf7iPsHYWwU\n9+Q41bEIxbEJSod2UzrSSnFbYnBSIA7+GPTHlOeSYdAfg4DTTTy/GEdJKTkV1eRW1VJY10BVMhDO\nDIVlZWXWZWxMGmVUwEueU/ch4B4S06TcpKo7ROQaYJuqbgZ+BvxCRJqBYRIhEODVwDUiEgViwOWq\nOpz+d5EdUv3KLn3ThezYsYP777+fvz3wALu2/oOqkV6O8woneh38mxuez1daPYX0rlzP5Mmn89q/\n3IbLPzjn+T21DdTU1Bx0OTPhvDRjZsqUz6TL5aL+6utStigec90P+X6KMgaDQf7whz+wdu3axDRC\nAwPs6O9jpLOdie4uQv3dRIcHYdSPc2KMBo1QMjlISXiQ0sEXKHke8h2JbmJ/DJrjyiPRRCAcUeE/\nipwUOvbuNdJQkBe+dAX/9FVSU1NDdXU1JSUlh2QibGOOZBkV8ABU9W7g7lnbvjDjdgi4IMVxdwJ3\nHvYCHgFmdzNFuttp/sR72PSByxgKBHmpV3ivV1AnPFOcQ6B2Jf0nvZKSc8/nvNNOo6ysbMZznXRI\nRl0aY/bf/rYo7hmQcuKJJ+7zuVWVQCBAf38/fX199Pf309zfT39PD2Nd7Uz0dhHu7yXmH8IxPoJn\nKkiBxIG5wc09MkjLO1/HozEYjCp+nMTyi5DSCtyV1eTVNlDS2ER1Tc10CKypqaGqqgqPx3NQdWRM\ntsq4gGeWVjgcpuXaq3DOWrHBrTHemhPlHpePyaNOoPfsN3LSuf/C2UcfveA8b9aNaszSOlwtiiJC\nQUEBBQUFrF69ep/7RyIRdmxcRzxFl3HAk8vWitXEhvpxjo+QHwlRNumnvNdP+cALlO9MTEUzHIOh\nqNISg8Hk7UmvDy0qwVVejbe6jvyGFVTVN0yHwD3XRUVF062Cdh6gORJYwDNMTk7yxz/+kV9v2kTP\nX/7Ed0vCKc+b8zgcfKF5dL+7TjKly8oYs3TcbjcN83QZv+S6G/jZjO+IYDBIb28vPT09iamUenvp\n7exkvKOVYHcH0cE+GBnEMzFObWSC8pEg5eNdlLc/TunjEIonAuDOqPJgDAZjyqi4iBYUc5TPw5mh\nAdwaBxI9FK2fuoye7h5WvO1SCgoK0l43xhwOFvCOUJFIhHvuuYdNmzbx2F3/yxmOEBcXOAjmQ1Cc\n+IjPOcZd02DnxRhjDthiW/Tz8vJYtWoVq1atWvD5YrEYg4ODe4XBJ7u78He0E+hsI9zfgw4P4Bzz\nUxwLUz4xxBkIbsfe32MSnmL0a1fzn5+7ihFnDvHiMpwV1eTWNlDYtIbahkbq6uqora2lrq6O6upq\nW5XEZDz7hGa52V0Rnos/yi929/Hrn/2EE4KDvLFAuLhY+Kevmp7zL+Lsyz5EwVMP2XlzxpjD4lC2\n6Dudzul5MY8//vgF9w0EAvT29jJ67rFMr888g1eEUXFRFpuiYrSHyokeKjv+SdljMB6H/ihsjSr9\nUWUwLoTyCqG0Ak91HfmNKylfsYraujrqZlxmdgvvkap7mMKDH3BmzGwW8LJYqsESE1+5gnVh5Wde\nYXt+EeGzzqPpY5/i7PVHvXjgypWAnTdnjMke+fn5rFmzhu21DfMucffT5hfw+/10dXXR3d1NV1cX\nj3V2Mry7mWBnG5G+bhgeIHcyQNPUGJUDY1T4W6h8/m+4BQZi0BZVtkUTU8qMiJtYUWmiNbB+BS91\nxdjwxH04oxHgxWliXG//MCzTOdlM5rKAl8W6v/H5vVrhADwOYbUvh9h3b+cDr3v9vF2udt6cMSYb\nLbTEnYhMr/973HHHzfsckUiE3t5eurq66Orq4smuLvra2xhr3UWou534YB/O0WEqY2EqAv1Uhvqp\n6NlOlZM537kaChL5ydd47+a7yW1YSfGqddStWEFDQ8P0pbi42E6PMfvNAl6W+sc//oG3uz3FhASQ\nG49w4rlvSHuZjDFmqR2Kkf1ut3s6fC1kbGxsOgR2dnZSde1lKfcrJM4r256kquspyrfBSAz6osqO\naOLa7/AQKy7DWVWHr7GJ0lVrqW9s3CsEFhYWpnxuGzF85LKAl2W6u7v55Cc/ydN33saNtU5cKX71\nuWsW/lIyxphslq4eisLCQgoLCzn66KMB2H7TtSm7hyPFZdR/9RY6Ojp4tK2N4ZbnmexoJdrbhWNk\nkNJYmOqxXqqCvVR1PE7hVhiMQldUeTwKvVFl1OUlVlKOp6aBgqY1VK9Zx0tGuqm96+bptcBt5ZAj\niwW8LHLnnXdy+SXv4W2uAF+rdtF61MmsbXsaDU1O72ODJYwxZmnM1z0cveBSzjnnnJTHqCqjo6N0\ndHTQ2dlJa0cHXa278be8wGTHbuL9PbhGh6mLh6ga7aJqopuq1ofJ2QIuSd0lvOPqy/nLw9upb1rF\nihUrWLFiBY2NjXi93sP59k2aWcBbxmY2vY/n+PhT1yg3FDoIltWy+qbfsfplJ1nzvDHGZIj5uoe3\nLzCKVkQoLi6muLh43vMCVZXh4WE6Ojpob2/n6bY2uluaecsfbki5vy8c5NxNX2c0lmgF3BaBnqgy\nkVsAFdXk1K+kaPV6alevmQ6AK1asoLi4eM5z2d+YzGUBb5lybb2P9l9+b/qXYEEowPtLnfQcfxpv\nvuO+6dUlbLCEMcZkjpTfyVu2HNRzighlZWWUlZVxwgknTG/f/tgfU3YJh4vK+NMb3stQ83ME21qI\n9XXhGR2memqcmr4ANUO7qNlxH1MKPVG4L6p0RxS/y0u8tBJ3bQP5TWs5nkmO/ftmHJEwYF3AmcYC\n3jLl/e3P54yQFYEVgx0LLh1mjDHmyDBfl/C6L36bV8wKYLFYjO7ubtra2mhvb+evra30Ne8ksLuZ\nSE8HzuF+yqIhavwdVI93UrPrIQoc4EjRBfz0VZey+a77qF13FE2rVtHU1ERTU1PKFsCZrDXw0LKA\nt0zJ8EDK7ZGejjSXxBhjTCbanxHDTqdzwZHBqsrQ0BDt7e20tbXxZFsbr/r+J1Pumx+d4rw/34zc\nCz0RuD+q3BpR/O5ctLyanIYmitYeTcO69dPhr2THI/R88aN7zdtqrYEHxwLeMjXpKyRvYnTOdhsh\na4wxZo9DdZqOiFBeXk55eTknnngiANt/+/2UXcCxkgoef9tVdD3/HCMvPEuksw3ncB8V0UnqBlup\nGW2j7rkthBS6Isr2KLwiV8h3zm0NbLn2KgLHv5L6+nqcTmfKslnLX2oW8JaheDzODaPKR5261zQo\nNkLWGGNMuszXBbz6C9/klFkBS1Xp7+9n9+7d7N69my0tLfQ+/yyBXc8T7WnnLIZSvoZjqI+fnbya\n3riDcEkF7voVFK49hvqjjmH1mjWs6NxJ/IfXTs8Wkarl70hdHm7RAU9E/gH8CLhdVacOV4FE5Fzg\nu4AT+KmqXj/r8RzgFuAkYAh4q6q2Jh/7NHAJEAM+oqr3HK5yLpWhzb+i5dor+bhzFBCksBgdH7Vf\nLcYYY9Jqf7qARWR63eBTTz11zuPbX72aSPfcU4zGnW46cnIoCI6xeryP+uZ+6lsfRe6FrggE3JDr\nmNvy1/ylK2iuWkXVC0/iv/6Te3X9tl1xMQXA404nZRe+lxVf+t4hqI3Msz8teBHgZuA7InIzcKOq\nPncoCyMiTuAHwNlAJ/CoiGxW1Wdm7HYJ4FfVNSJyIfBfwFtF5BjgQuBYoBa4V0TWqWrsUJZxKe1Z\nW9YZCiZGVACEw6z45v9YsDPGGJN2h6oLuPbKa1O2Bh533Q38/LyLmJiYoKWlhZaWFh7btYuO555h\n9Pln+GDntpTP5xoZZOu/ncZLvYLPMXfCfwGIxRi69ccAWRnyFj3cUlU3AkeTCHnvBHaIyBYReauI\nuA9ReU4BmlW1RVXDwCbg/Fn7nJ8sA8AdwJmSmMnxfGCTqk6p6m6gOfl8WSPV2rIaCtL9jc8vUYmM\nMcaYg1d23kU0XncD7tpGEMFd20jjdTdMh0efz8dxxx3H+eefzxVXXMG3b/wpN215EE9tY8rnC7i9\nPFy2krxFrOE7tOmnh/S9ZIr9OgdPVXcCVyS7Qv8duAy4DRgUkZ+TaNVrOYjy1AEz22g7gZfPt4+q\nRkVkFChLbt8669i6VC8iIpcly05VVRVbDnIOonQp6O5IubZsuLtj2byH5SwQCFg9p5nVefpZnaef\n1XlSYQ1c82LYGoJ9zhGYs+6leGatu64AJ2/k1Fe8Hv32xxHVBZ9DY7GsrP8DGmSRPAfvFyKyA/gW\n8Grgk8CVIvK/wIdVtfcAnjpVfpn9LzPfPos5NrFR9UbgRoANGzboxo0b96OIS2d7bUPKEUue2gaW\ny3tYzrZs2WL1nGZW5+lndZ5+VucJ+zMaNhaL0dnZSf9Tj8z54y+A/P3/mLr3T0ievHhK0zzE6czK\n+t/vgCciucBFwOUkBjo8B3wU+A3wJuCLwK3AmQdQnk5g5jwf9UD3PPt0iogLKAKGF3nssjbfiCUb\nOWuMMWY523OO+d7z4F1OZ2cnbfXr2dXcTMdzzzD2wrNEOltx+weoccT514LUAc7nELaXNXBidICC\n6MLjQssufO9heU9LbX9G0R4HvA94G+ADfg9crar3z9jtJyLSSyLsHYhHgbUi0gR0kRg08R+z9tkM\nXAw8BLwF+IuqqohsBm4TkW+RGGSxFnjkAMuRkfb8kum4/tNE+7qIIxS+9VIbYGGMMWbZmZqaor29\nnd27d1PwxY+RM+cc80kC3/gMgQic6Eq0KHVGoTOidKnSlVPIuHOKwnh0znN76hq56YHmOcFx+rlJ\ntNzZKNqEJ0m0iH2HxLl2PfPs10wifO235Dl1HwLuITFNyk2qukNErgG2qepm4GckuoebSbTcXZg8\ndoeI3A48A0SBD2bTCNo99oxY2rhxIxuf+Tvn3vI9Hr/lv22aFGOMMWm3ULdqNBqlq6treu671pYW\nBp5/honWZuI9neQFRqh1CbVuONGbuiXO6xDuW72BorVHU7f+aFavWcPrV69m1apV+Hy+lAFuZs/W\nfFO5bC+sycpu2Zn2J+BdAPxuX6FJVZ8FzjjQAqnq3cDds7Z9YcbtULIsqY79CvCVA33t5eTqV55A\nUceDeDQO2LIuxhhj0mvwd7fR/tn3w9SLkwzvuvI9/PiLn+PZkQk8o0NUO+LUuYVal3CUC8bj0BVV\nuiPQ7RDafSV0VtezfrSF/PDknNfw1Dbysy0PzluGxczFl3IqlywcVDHbogOeqt55OAti9k/FvXcw\ne2qfPVOmWMAzxhhzsOLxOH19fbS1tdHe3k7r7t307nyGsZbniXR38O5IP4WzVg9zxWO8dqQdV1Dp\nFqUvt5D+ylqeXbGakrVHs2LtOpqamnjNypU0NDTgdidmWdtXS9xCDtVcfNnGlipbpmR4IOX2SM/c\nmcCNMcYcmRbqQg0Gg3R0dNDe3k5bWxsdra0MN+8k2LaLeH83ntFhqhxxql1Q4xJOd8GkQndE6YlC\ngU9INYGF2+Hg/L8+xYoVK/B6vYsq5/6simEWxwLeMqWlFchw/5ztYYWdf9/C+tM3pr9QxhhjMoKq\nsvuXP8b/1auQcGIUaaS7neZPvJsffeoT7ByZIH8qQI1LqEkGuJe4YDgGvckA16NKuzuftvJqPHWN\nFK9eT0OyBe6ElSvRD5xPtKdzzmt7ahtYv379fpfZWuIOLQt4y1Toze/G98vv7dWcHVLh/kCMU99x\nDv945Ru48Ce3MXnv7+0XkTHGZJlQKJSYQiTZfdre2spg804m2lqI9HXh8g/ynoI4PufeLWxujfPa\n0ADhqNIvTlrzS2mvqiW3cRUla9bT2LSKxsZGTmtspKGhgby8vHnLMHTVV2zqrgxmAW+Zip56Jo3H\nHLNXeKv8wKfZefdf+fmvf8GXHvgjd7+kjJUeB45oBLCBGMYYs9RSdZlSWLPXPpFIhO7ubjo7OxOX\njg4GWpoZb20m3N2BDvWRFxyn2iVUu6DKJRztgrEY9EaVvij0oeQ5nCnL4HE4uPzJVqqqqnA4Fr1i\n6RzWrZrZLOAtY6mas2+58BL+dumlfOxDH+JbY8/giO79680GYhhjzNJIDCS4HA29OOq05ZPv5a9V\nR3Gbr5Rwdwcy1E9ucIwqJ1Qlw9srXBBSEsEtqvRFlX6E5/JKaa6qJbehiZLV62lclWh925BsfWt9\n48vmXf2opqZmzvYDYd2qmcsCXhZ61atexcOPPcaTR6VuWreBGMYYc+hNTk6+2Oq2p+Vt9y4C7S1M\n9XTy7kA7hY69V9B0RiO8sfMp2iKJ8NYbhT5VWrxFtJRVkVPbQMHKNdQ0raKuro6TGxtZsWIFNTU1\nOJ2pW+j2sNWPjmwW8LKUy+XCU9uY8tdbOK5848xTOOXjn+NVb3wTso91+owxJpvsz5qnkJguZGBg\ngO7u7ulLb0c7o627CHa2EenrxjEyiG8qSIVTqHRBhUtY40y0vA1EoT+mFOSmHnUqIsS/eSsvq6+n\nvr6e6urq6elDDoZ1oR7ZLOBlsVS/3sLi4FZ/lMqnH0c//P+48SNewq84i9M+cjUvO+Xl02Fvf78A\njTFmOUi15mnbp9/H9iefYnf9Orq7uhhq281ERyuRvi4YHsATGKXcCZXORHBb5YLjHDCUDG4DUeiP\nKkPqZMBXws7yanLrGilauZqalU3U19dzcn09zisuJN7XNadMWlrJBReknL//oFkX6pHLAl4WS/Xr\nbcWVX+bK41/JT37yEy69+eesH+3hjff/gZEH7+LrrkKcZ/wLp61dRc7N35616LMNzjDGLA+qytjY\n2F4tbt1dXQy0t/Gmv9xKQTQ6Us2qAAAUM0lEQVS09wFTIXJ+/g3KYrDemZhuqj8GA1GlPwoDqnS4\nfLQXl+GqrCW3vpGSFauoraujtraW45Itb5WVlQsOWhi6+rqUXaahN7/7cFWFOYJZwMtyqX69lQHX\nXnst11xzDQ8++CCbNm3i+3f8mlOmhnnj/92G4z5QscEZxphD72B6B6LRKAMDA/T29tLX15e47ulh\npH03wa4OQn1dxIcGcI75KdQI5U6h3AnlLmGVE+JAnpByzVOPQ7j3pDdSsGIVVY0rqK2t5WW1tdTW\n1lJdXb3oCXsXstC6qMYcahbwjmAOh4PTTz+d008/neh3vsOWLVv4w+bNvPOPP0q5f7i7nbt+/zte\ntfEMioqK0lxaY8xyl6p7tP0z72d8bIzACafR29s7Hd4GOjsIdLYx1ddNbLAPGfWTG5qg1Ml0aKt1\nwjFOGI3BUAwGY8pgFIZiSq/DTW9BGc7yKnJr6ilsXEll40pOu/XruEeH5pTNU9vIf//mfw97HRyp\n66Ka9LOAZ4DEoIyzzjqLs846i+3/vDvl4IyoKsUfv4DfTEFfSRXuY09k5dn/wslnns3KlSttsIYx\nBkh0kY6MjEy3su0JbSf+9Bp8M7onIdE70P75D/GPoFLmgnKncLIT3AKDydA2FE1eK7S689ldXIaz\nrBJvTT2F9Y1U1tRSVVXFMcku09raWgoKClJ+Jw2tqbGRpeaIYAHPzJFqcIZ6cnjixLO5Y1cXoacf\n59jJHo4fvpvqR/7EP6+B2yWHqRXrKD31Naw942xOOPEkqqqqABuwYcxSOlT//0KhEAMDA/T39zMw\nMJC43dfLaFcnwZ4uQgM9RIcG0NFhHIExiohR6hSKnVDqFI5zQp6DlN2jPofQXFzNrtJycqpq8dWv\noLSugeqaGmqqqzmhqorq6moqKipwuQ7uz5aNLDVHiowJeCJSCvwaWAm0Av+uqv4U+10MfC5591pV\nvTm5fQtQA0wmHztHVecu1mr2ab4vwJPOu4jLgUAgwNatW9m6dStff+gheh59kBVTY7x08inqWreT\n+5vvc+eU0uL2UV9ZySvHunDFooAN2DAmnebrEgXwnfNvKQJbH9sfeZjNN/6AUF8PseFBdMyPc2KM\n/GiYEqdQkgxsVc7EgISJOPhj4I8pw8lrfwx63Dn0FJTgKCnDW1lDfm0D5z/8e/ImRueU01PXyM8e\naE5bvdjIUnMkyJiAB3wKuE9VrxeRTyXvXz1zh2QI/E9gA6DAYyKyeUYQfJuqbktnobPVQl+A+fn5\n0925kOiO2bVrF1u3buX+Rx7hhcceJb7zKdZNBXi5M4jLMXfAxjOfvIwHH3iYpg2ncvSxx7JmzRpy\nc3MP+/syJhsFg0GGhoYYHBxkaGho+vLSH3+BvBRdok9/7GJ+P/Z2ip1CaTKwlTlhrRNOVhiJwXAy\nsO253YuTrpwCKCzBXVZJTmU1+TV1lFfXUFFRQUVFBcckW9qqqqpSrmE6tPlM6x41Jk0yKeCdD2xM\n3r4Z2MKsgAe8Dvizqg4DiMifgXOBX6WniCYVEWHNmjWsWbOGt7/97UBiYtCWlhZGzz0m5TF5kRAv\nveMHFP/2B7wQhfsjyojXR7SsGk9DE8Xrj6X2hA2sOupompqa5h3UYd2/JtMczGcyEokwMjKC3+/H\n7/e/GNYGBhjr7Wair4fQYD8R/yCx0REIjOKcnMAXj1LoFIocUOiEQodwlBNy5+kSzRfFr04Gcgqg\nsBhXWSU5FdX4auoYnQzx8pe/nIqKCtZVVk6Ht6KiooM+z9a6R41Jn0wKeFWq2gOgqj0iUplinzpg\n5jpbnclte/xcRGLAnSS6b/deE8akjcPhYM2aNWyfZzWNSHE5D5zzLp5/ejv+556GkS5qJgPUB5qp\n6N5F/WP3Uf0b6I7BIxGlDzfBwlK0spacxlUUrz+WdRODVP/+50h4KvGc1v1rlliqtUZ3X30Z//j7\n32lrOAr/0BCTA32EBvsJ+weJjfrR8VEcwQDuUBBvLEyRQyh0kghrDqHOCesdEIzDWBxGY5q8htG4\nMhaDPoeLnpx8HAVFOItL8ZSWk1tRxVu23UVecHxOOd01DdzUvCtlYNuyZQsbN248bHVk3aPGpIek\nMwOJyL1AdYqHPgvcrKrFM/b1q2rJrOOvAnJU9drk/c8DQVX9pojUqWqXiBSQCHi/VNVb5inHZcBl\nAFVVVSdt2rTpULy9tAoEAuTn5y91MfbJtfU+cm/5znQIg8SAjcl3fozoqWdOb4vFYvT19dHV1TV9\n6e3qYqq7HY9/gCqiNLih3iXUu4V6N+QIOFP8gRoTF7c2vIyc8ipyK6rJr6ikqLiYkpISiouLKSgo\nWHAy0vkEAgGKn34Y729/jgwPoKUVhN787r3ehzm09vU5d22975D9e8TjcSYnJwkGgwSDQSYmJva+\nDgSIBMaIjo8RC4yjwQAEA+hkEEcoiGNqkotzpihIsTxoVJVgPDHIYCIZzsbiymgcxmYEtTEVptw5\nRHJyieX6IL8IZ1EJ7uJSCouLKSwspKioiMLCwulLUVERXq83ZVhb7P+/mZbLd0s2sTpPv+Vc52ec\nccZjqrphX/ulNeAtRER2AhuTrXc1wBZVXT9rn4uS+7wvef/Hyf1+NWu/dwEbVPVD+3rdDRs26LZt\ny++0vcP9K/tQOthuVFVlaGiItrY22traaG9vp621lbff9cMUqzom9n92SilOjuBzS+I8opEYjMSV\nkRhMeRJ/QLWgCCkswVlShqesgtzKGnzVtRRVVlFaVkZJSQklJSUUFRXx9A3fpOBXP5hz/lDjdTek\npUUiU7qj01mO+T7n8Xic3jtupu9LH0OnJqe3qycH/wWXM3TMBiYmJhLhbGyU8KifqZHEJTo2QjQw\nRnxiHIITSGgSZ3gSdyRMnkPId4AveclzCPmy5zZEFALxREibiCsTuud24v4FhY6UQUuBX57xDnwV\nVRSXlk5/rvZcipM/QA70x8dC9vffazl9t2QLq/P0W851LiLLLuB9HRiaMciiVFU/OWufUuAx4MTk\npseBk4AxoFhVB0XETeKcvHtVNfWMvTNYwFu+tr96Ter5+koq+Oel/0l3dzc9PT0M9/US7OshPNRP\nzD+MazIwHf6KHbx42ynJ+4nlwEfie4JhIhSe7hN8jrl/vANOD3eteyVuXz7uPB8eXz6e/AJy8gvI\nKSgkt6iIvPwC8vLy8Hq95OTk7PMy+4/87NGQkN5wuXA5cim86nqcr3k9U1NTi7qEp6aIBCeIBCeI\nTgaJTQaJTgaJhyaJTASITgaJBoOMDQ3gdTrQqSk0HELDUxAO44xFubDIQb5z7r9HRJXeaDKkCTgl\nEcqCM4JYIA5B3fv+RByibg/qzYNcHw5fAa6CQtyFReQUl+ItKaMg2Yo2+7InqA1d+CqiPZ1zyuSu\nbeS4NI4SPRj23ZJ+Vufpt5zrfLEBL5POwbseuF1ELgHagQsARGQDcLmqvldVh0Xky8CjyWOuSW7z\nAfckw50TuBf4SfrfgkmnVPP1iTeP1Z//BqcsEHrC4TCDg4PTU0QMDw/j9/tp8/t5Inl7fGiA8OAA\nsZEhGB/FNTnBORJK+Xy+6BTrH/8LOZJY7ihHwCOJLuQ9lzgwpYk1LkcUwvHE7ZAq4eT2qVnXEYQw\nQkSEC3xKwayGHQ0F2f7xd3HrRz+MQxw4HIJDJHnbkbid3Aag8TioovE4qopq4n58znYFjaNx3Ws/\nVLnQp3O6IDU0Sc+XPsIDE4pnz3t3kKyHF+vDJ7xYNw5hKq5z3nfi9ovbp+Iz60Rf3AdShm0Alwh/\nPP4sXAWFuAqK8BYU4svPp6CgYDqQ1aYIaQUFBTidKfpX94P3qq/YKFFjTEbImICnqkPAnJNCktOe\nvHfG/ZuAm2btM0GiJc8cQQ50RJ7H45me7X5/bDulHhmeO7VirLSCnI9+jWAwyOCs87eCwSDBiQlC\nwQmmAuNEg0F0KpRojZqagkiiRYpoGEckgkQjOOPRZDDS6UCULw5I0SFdQJyqiZFEEJuxXWdd0MTR\nqdrrHTOP0Rf3UebeznekLkeuCM/mlhB3u8HtAXcO4vEgOV4kx4sjx4vDm4szeXHl5pEzT2tmXl4e\nJT4fPp+P559/ntecdho+n4+8vDx8ye1er5cdG9elbMH11DZy4+a7U7zTw89GiRpjMkXGBDxjDkQ6\nR+SF3vxufL/8XooWw28u2GK4vzTZqhaPx4nFYsRiMXade1zKrj9ndR2X/uWvxGKxvfbfc3vPtYjg\ncDhwOp0HdL3ndvM5LyHS0zGnHJ66Rn5xGLog8/LyOPnkk1M+Nl8L7lK3ltkoUWNMJrCAZ8wiRU89\nk8ZjjjnsrTMigtPpxOl04na7Aaibp+uv4eqvUrZixSF9/YXUXnVtxoQqay0zxpj5WcAzZj8sVetM\npoSZTCnHzPJYoDPGmLks4BmzTGRKmMmUchhjjJnfoZ1wyRhjjDHGLDkLeMYYY4wxWcYCnjHGGGNM\nlrGAZ4wxxhiTZSzgGWOMMcZkGQt4xhhjjDFZxgKeMcYYY0yWsYBnjDHGGJNlLOAZY4wxxmQZC3jG\nGGOMMVkmYwKeiJSKyJ9F5IXkdck8+/2fiIyIyF2ztjeJyMPJ438tIp70lNwYY4wxJrNkTMADPgXc\np6prgfuS91P5OvCOFNv/C/h28ng/cMlhKaUxxhhjTIbLpIB3PnBz8vbNwL+m2klV7wPGZ24TEQFe\nC9yxr+ONMcYYY7JdJgW8KlXtAUheV+7HsWXAiKpGk/c7gbpDXD5jjDHGmGXBlc4XE5F7geoUD332\nYJ86xTZdoByXAZcl7wZEZOdBvv5SKAcGl7oQRxir8/SzOk8/q/P0szpPv+Vc5ysWs1NaA56qnjXf\nYyLSJyI1qtojIjVA/3489SBQLCKuZCtePdC9QDluBG7cj+fPOCKyTVU3LHU5jiRW5+lndZ5+Vufp\nZ3WefkdCnWdSF+1m4OLk7YuB3y/2QFVV4H7gLQdyvDHGGGNMNsmkgHc9cLaIvACcnbyPiGwQkZ/u\n2UlE/gb8BjhTRDpF5HXJh64GrhCRZhLn5P0sraU3xhhjjMkQae2iXYiqDgFnpti+DXjvjPuvmuf4\nFuCUw1bAzLOsu5iXKavz9LM6Tz+r8/SzOk+/rK9zSfRuGmOMMcaYbJFJXbTGGGOMMeYQsICX4UTk\nXBHZKSLNIjJndQ8RyUkuzdacXKptZfpLmV0WUedXiMgzIvKUiNwnIosasm7mt686n7HfW0RERSSr\nR7+lw2LqXET+PflZ3yEit6W7jNlmEd8tjSJyv4g8kfx+ecNSlDNbiMhNItIvIk/P87iIyH8n/z2e\nEpET013Gw8kCXgYTESfwA+D1wDHARSJyzKzdLgH8qroG+DaJJdvMAVpknT8BbFDVl5JYPeVr6S1l\ndllknSMiBcBHgIfTW8Lss5g6F5G1wKeB01T1WOBjaS9oFlnk5/xzwO2q+jLgQuCH6S1l1vkf4NwF\nHn89sDZ5uQy4IQ1lShsLeJntFKBZVVtUNQxsIrGk20wzl3i7g8To4lQTP5vF2Wedq+r9qhpM3t1K\nYt5Fc+AW8zkH+DKJMB1KZ+Gy1GLq/FLgB6rqB1DV/Zmb1My1mDpXoDB5u4gF5nM1+6aqDwDDC+xy\nPnCLJmwlMZ9uTXpKd/hZwMtsdUDHjPuplmCb3ic5yfMoiWlizIFZTJ3PdAnwp8Naouy3zzoXkZcB\nDap6VzoLlsUW8zlfB6wTkX+IyFYRWaglxOzbYur8i8DbRaQTuBv4cHqKdsTa3+/7ZSVjpkkxKS1m\nCbb9WqbN7NOi61NE3g5sAF5zWEuU/RascxFxkDj94F3pKtARYDGfcxeJrquNJFqp/yYiL1HVkcNc\ntmy1mDq/CPgfVf2miLwC+EWyzuOHv3hHpKz++2kteJmtE2iYcT/VEmzT+4iIi0Sz/kJN0mZhi6lz\nROQsEmson6eqU2kqW7baV50XAC8BtohIK3AqsNkGWhyUxX63/F5VI6q6G9hJIvCZA7OYOr8EuB1A\nVR8CvCTWTDWHx6K+75crC3iZ7VFgrYg0iYiHxEm3m2ftM3OJt7cAf1Gb3PBg7LPOk92FPyYR7uy8\npIO3YJ2r6qiqlqvqSlVdSeK8x/OSk6CbA7OY75bfAWcAiEg5iS7blrSWMrssps7bSU74LyJHkwh4\nA2kt5ZFlM/DO5GjaU4FRVe1Z6kIdKtZFm8FUNSoiHwLuAZzATaq6Q0SuAbap6mYSS7L9IrlE2zCJ\nLw1zgBZZ518H8oHfJMeztKvqeUtW6GVukXVuDqFF1vk9wDki8gwQA65KrjhkDsAi6/wTwE9E5OMk\nugrfZT/YD5yI/IrEKQblyfMa/xNwA6jqj0ic5/gGoBkIAu9empIeHraShTHGGGNMlrEuWmOMMcaY\nLGMBzxhjjDEmy1jAM8YYY4zJMhbwjDHGGGOyjAU8Y4wxxpgsYwHPGGOMMSbLWMAzxhhjjMkyFvCM\nMcYYY7KMBTxjjDlIIuITkedE5BERcc/Yfo6IxEXkg0tZPmPMkcdWsjDGmEMguUbxVuDbqvopEakE\nngIesaXsjDHpZgHPGGMOkeQaot8EzgGuBI4DjlfVwSUtmDHmiGMBzxhjDhEREeCPwGsBD3C2qt63\ntKUyxhyJ7Bw8Y4w5RDTxi/kXQA7wpIU7Y8xSsYBnjDGHiIhUA98BHgeOF5GPLnGRjDFHKAt4xhhz\nCCS7Z28GwsDZJILef4nIS5e0YMaYI5Kdg2eMMYeAiHwC+BrwWlX9q4h4SIyqzQE2qOrkkhbQGHNE\nsRY8Y4w5SMkpUq4DvqqqfwVQ1TBwEbAS+NbSlc4YcySyFjxjjDHGmCxjLXjGGGOMMVnGAp4xxhhj\nTJaxgGeMMcYYk2Us4BljjDHGZBkLeMYYY4wxWcYCnjHGGGNMlrGAZ4wxxhiTZSzgGWOMMcZkGQt4\nxhhjjDFZ5v8DN/43tQSZlzsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1162fea90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot discretized geometry\n",
    "width = 10\n",
    "pyplot.figure(figsize=(width, width))\n",
    "pyplot.grid()\n",
    "pyplot.xlabel('x', fontsize=16)\n",
    "pyplot.ylabel('y', fontsize=16)\n",
    "pyplot.plot(x, y, color='k', linestyle='-', linewidth=2)\n",
    "pyplot.plot(numpy.append([panel.xa for panel in panels], panels[0].xa),\n",
    "            numpy.append([panel.ya for panel in panels], panels[0].ya),\n",
    "            linestyle='-', linewidth=1, marker='o', markersize=6, color='#CD2305')\n",
    "pyplot.axis('scaled', adjustable='box')\n",
    "pyplot.xlim(-0.1, 1.1)\n",
    "pyplot.ylim(-0.1, 0.1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Free stream conditions\n",
    "\n",
    "The airfoil is immersed in a free-stream $(U_\\infty$,$\\alpha)$ where $U_\\infty$ and $\\alpha$ are the velocity magnitude and angle of attack, respectively. Like before, we create a class for the free stream, even though we will only have one object that uses this class. It makes it easier to pass the free stream to other functions later on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Freestream:\n",
    "    \"\"\"\n",
    "    Freestream conditions.\n",
    "    \"\"\"\n",
    "    def __init__(self, u_inf=1.0, alpha=0.0):\n",
    "        \"\"\"\n",
    "        Sets the freestream speed and angle (in degrees).\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        u_inf: float, optional\n",
    "            Freestream speed;\n",
    "            default: 1.0.\n",
    "        alpha: float, optional\n",
    "            Angle of attack in degrees;\n",
    "            default 0.0.\n",
    "        \"\"\"\n",
    "        self.u_inf = u_inf\n",
    "        self.alpha = numpy.radians(alpha)  # degrees to radians"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# define freestream conditions\n",
    "freestream = Freestream(u_inf=1.0, alpha=4.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Flow tangency boundary condition\n",
    "\n",
    "A constant vortex strength $\\gamma$ will be added to each panel (all panels have the same, constant vortex-sheet strength). Thus, using the principle of superposition, the velocity potential becomes:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\phi\\left(x_{c_i},y_{c_i}\\right) &= V_\\infty x_{c_i} \\cos \\alpha + V_\\infty y_{c_i} \\sin \\alpha \\\\\n",
    "&+ \\sum_{j=1}^N \\frac{\\sigma_j}{2\\pi} \\int_j \\ln \\left(\\sqrt{(x_{c_i}-x_j(s_j))^2+(y_{c_i}-y_j(s_j))^2} \\right) {\\rm d}s_j \\\\\n",
    "&- \\sum_{j=1}^N \\frac{\\gamma}{2\\pi} \\int_j \\tan^{-1} \\left(\\frac{y_{c_i}-y_j(s_j)}{x_{c_i}-x_j(s_j)}\\right) {\\rm d}s_j\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "The flow tangency boundary condition is applied at every panel center:\n",
    "\n",
    "$$0 = \\vec{V}\\cdot\\vec{n}_i = \\frac{\\partial}{\\partial n_i} \\left\\{ \\phi\\left(x_{c_i},y_{c_i}\\right) \\right\\}$$\n",
    "\n",
    "i.e.\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "0 &= V_\\infty \\cos \\left(\\alpha-\\beta_i\\right) + \\frac{\\sigma_i}{2} \\\\\n",
    "&+ \\sum_{j=1,j\\neq i}^N \\frac{\\sigma_j}{2\\pi} \\int_j \\frac{\\partial}{\\partial n_i} \\ln \\left(\\sqrt{(x_{c_i}-x_j(s_j))^2+(y_{c_i}-y_j(s_j))^2} \\right) {\\rm d}s_j \\\\\n",
    "&- \\sum_{j=1,j\\neq i}^N \\frac{\\gamma}{2\\pi} \\int_j \\frac{\\partial}{\\partial n_i} \\tan^{-1} \\left(\\frac{y_{c_i}-y_j(s_j)}{x_{c_i}-x_j(s_j)}\\right) {\\rm d}s_j\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "We already worked the first integral in the previous lesson:\n",
    "\n",
    "$$\\frac{\\partial}{\\partial n_i} \\ln \\left(\\sqrt{(x_{c_i}-x_j(s_j))^2+(y_{c_i}-y_j(s_j))^2} \\right) = \\frac{\\left(x_{c_i}-x_j\\right)\\frac{\\partial x_{c_i}}{\\partial n_i} + \\left(y_{c_i}-y_j\\right)\\frac{\\partial y_{c_i}}{\\partial n_i}}{\\left(x_{c_i}-x_j\\right)^2 + \\left(y_{c_i}-y_j\\right)^2}$$\n",
    "\n",
    "where $\\frac{\\partial x_{c_i}}{\\partial n_i} = \\cos \\beta_i$ and $\\frac{\\partial y_{c_i}}{\\partial n_i} = \\sin \\beta_i$, and\n",
    "\n",
    "$$x_j(s_j) = x_{b_j} - s_j \\sin \\beta_j$$\n",
    "\n",
    "$$y_j(s_j) = y_{b_j} + s_j \\cos \\beta_j$$\n",
    "\n",
    "We now need to derive the last integral of the boundary equation:\n",
    "\n",
    "$$\\frac{\\partial}{\\partial n_i} \\tan^{-1} \\left(\\frac{y_{c_i}-y_j(s_j)}{x_{c_i}-x_j(s_j)}\\right)= \\frac{\\left(x_{c_i}-x_j\\right)\\frac{\\partial y_{c_i}}{\\partial n_i} - \\left(y_{c_i}-y_j\\right)\\frac{\\partial x_{c_i}}{\\partial n_i}}{\\left(x_{c_i}-x_j\\right)^2 + \\left(y_{c_i}-y_j\\right)^2}$$\n",
    "\n",
    "where $\\frac{\\partial x_{c_i}}{\\partial n_i} = \\cos \\beta_i$ and $\\frac{\\partial y_{c_i}}{\\partial n_i} = \\sin \\beta_i$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Enforcing the Kutta condition\n",
    "\n",
    "To enforce the *Kutta-condition*, we state that the pressure coefficient on the first panel must be equal to that on the last panel:\n",
    "\n",
    "$$C_{p_1} = C_{p_{N}}$$\n",
    "\n",
    "Using the definition of the pressure coefficient $C_p = 1-\\left(\\frac{V}{U_\\infty}\\right)^2$, the Kutta-condition implies that the magnitude of the velocity at the first panel center must equal the magnitude of the last panel center:\n",
    "\n",
    "$$V_1^2 = V_N^2$$\n",
    "\n",
    "Since the flow tangency condition requires that $V_{n_1} = V_{n_N} = 0$, we end up with the following *Kutta-condition*:\n",
    "\n",
    "$$V_{t_1} = - V_{t_N}$$\n",
    "\n",
    "(the minus sign comes from the reference axis we chose for the normal and tangential vectors).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's derive the tangential velocity for every panel, since it will be useful to compute the pressure coefficient.\n",
    "\n",
    "$$V_{t_i} = \\frac{\\partial}{\\partial t_i} \\left(\\phi\\left(x_{c_i},y_{c_i}\\right)\\right)$$\n",
    "\n",
    "i.e.,\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "V_{t_i} &= V_\\infty \\sin \\left(\\alpha-\\beta_i\\right) \\\\\n",
    "&+ \\sum_{j=1,j\\neq i}^N \\frac{\\sigma_j}{2\\pi} \\int_j \\frac{\\partial}{\\partial t_i} \\ln \\left(\\sqrt{(x_{c_i}-x_j(s_j))^2+(y_{c_i}-y_j(s_j))^2} \\right) {\\rm d}s_j \\\\\n",
    "&- \\sum_{j=1,j\\neq i}^N \\frac{\\gamma}{2\\pi} \\int_j \\frac{\\partial}{\\partial t_i} \\tan^{-1} \\left(\\frac{y_{c_i}-y_j(s_j)}{x_{c_i}-x_j(s_j)}\\right) {\\rm d}s_j \\\\\n",
    "&-\\frac{\\gamma}{2}\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "\n",
    "which gives\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "V_{t_i} &= V_\\infty \\sin \\left(\\alpha-\\beta_i\\right) \\\\\n",
    "&+ \\sum_{j=1,j\\neq i}^N \\frac{\\sigma_j}{2\\pi} \\int_j \\frac{\\left(x_{c_i}-x_j\\right)\\frac{\\partial x_{c_i}}{\\partial t_i} + \\left(y_{c_i}-y_j\\right)\\frac{\\partial y_{c_i}}{\\partial t_i}}{\\left(x_{c_i}-x_j\\right)^2 + \\left(y_{c_i}-y_j\\right)^2} {\\rm d}s_j \\\\\n",
    "&- \\sum_{j=1,j\\neq i}^N \\frac{\\gamma}{2\\pi} \\int_j \\frac{\\left(x_{c_i}-x_j\\right)\\frac{\\partial y_{c_i}}{\\partial t_i} - \\left(y_{c_i}-y_j\\right)\\frac{\\partial x_{c_i}}{\\partial t_i}}{\\left(x_{c_i}-x_j\\right)^2 + \\left(y_{c_i}-y_j\\right)^2} {\\rm d}s_j \\\\\n",
    "&-\\frac{\\gamma}{2}\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "where $\\frac{\\partial x_{c_i}}{\\partial t_i} = -\\sin \\beta_i$ and $\\frac{\\partial y_{c_i}}{\\partial t_i} = \\cos \\beta_i$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Building the linear system\n",
    "\n",
    "Here, we build and solve the linear system of equations of the form\n",
    "\n",
    "$$[A][\\sigma,\\gamma] = [b]$$\n",
    "\n",
    "where the $N+1 \\times N+1$ matrix $[A]$ contains three blocks: an $N \\times N$ source matrix (the same one of Lesson 10), an $N \\times 1$ vortex array to store the weight of the variable $\\gamma$ at each panel, and a $1 \\times N+1$ Kutta array that repesents our Kutta-condition.\n",
    "\n",
    "We are going to re-use the function `integral` from [Lesson 10](10_Lesson10_sourcePanelMethod.ipynb) to compute the different integrals with the SciPy function `integrate.quad`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def integral(x, y, panel, dxdk, dydk):\n",
    "    \"\"\"\n",
    "    Evaluates the contribution from a panel at a given point.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    x: float\n",
    "        x-coordinate of the target point.\n",
    "    y: float\n",
    "        y-coordinate of the target point.\n",
    "    panel: Panel object\n",
    "        Panel whose contribution is evaluated.\n",
    "    dxdk: float\n",
    "        Value of the derivative of x in a certain direction.\n",
    "    dydk: float\n",
    "        Value of the derivative of y in a certain direction.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    Contribution from the panel at a given point (x, y).\n",
    "    \"\"\"\n",
    "    def integrand(s):\n",
    "        return (((x - (panel.xa - numpy.sin(panel.beta) * s)) * dxdk +\n",
    "                 (y - (panel.ya + numpy.cos(panel.beta) * s)) * dydk) /\n",
    "                ((x - (panel.xa - numpy.sin(panel.beta) * s))**2 +\n",
    "                 (y - (panel.ya + numpy.cos(panel.beta) * s))**2) )\n",
    "    return integrate.quad(integrand, 0.0, panel.length)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first define the function `source_contribution_normal` to build the source matrix, i.e. the weight from each source on the normal velocity of each panel:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def source_contribution_normal(panels):\n",
    "    \"\"\"\n",
    "    Builds the source contribution matrix for the normal velocity.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    panels: 1D array of Panel objects\n",
    "        List of panels.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    A: 2D Numpy array of floats\n",
    "        Source contribution matrix.\n",
    "    \"\"\"\n",
    "    A = numpy.empty((panels.size, panels.size), dtype=float)\n",
    "    # source contribution on a panel from itself\n",
    "    numpy.fill_diagonal(A, 0.5)\n",
    "    # source contribution on a panel from others\n",
    "    for i, panel_i in enumerate(panels):\n",
    "        for j, panel_j in enumerate(panels):\n",
    "            if i != j:\n",
    "                A[i, j] = 0.5 / numpy.pi * integral(panel_i.xc, panel_i.yc, \n",
    "                                                    panel_j,\n",
    "                                                    numpy.cos(panel_i.beta),\n",
    "                                                    numpy.sin(panel_i.beta))\n",
    "    return A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we define `vortex_contribution_normal` to build the vortex matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def vortex_contribution_normal(panels):\n",
    "    \"\"\"\n",
    "    Builds the vortex contribution matrix for the normal velocity.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    panels: 1D array of Panel objects\n",
    "        List of panels.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    A: 2D Numpy array of floats\n",
    "        Vortex contribution matrix.\n",
    "    \"\"\"\n",
    "    A = numpy.empty((panels.size, panels.size), dtype=float)\n",
    "    # vortex contribution on a panel from itself\n",
    "    numpy.fill_diagonal(A, 0.0)\n",
    "    # vortex contribution on a panel from others\n",
    "    for i, panel_i in enumerate(panels):\n",
    "        for j, panel_j in enumerate(panels):\n",
    "            if i != j:\n",
    "                A[i, j] = -0.5 / numpy.pi * integral(panel_i.xc, panel_i.yc, \n",
    "                                                     panel_j,\n",
    "                                                     numpy.sin(panel_i.beta),\n",
    "                                                     -numpy.cos(panel_i.beta))\n",
    "    return A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's call both of these functions to create the source matrix `A_source` and the vortex matrix `B_vortex`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A_source = source_contribution_normal(panels)\n",
    "B_vortex = vortex_contribution_normal(panels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We still have to include the Kutta condition in our system.\n",
    "\n",
    "Following this lesson, there is an [exercise](11_Lesson11_Exercise.ipynb) for you to complete, which will help you to show that:\n",
    "\n",
    "* the matrix resulting from the contribution of the vortices on the normal velocity, $B^n$, is **exactly** the matrix resulting from the contribution of the sources on the tangential velocity, $A^t$, i.e.\n",
    "\n",
    "$$B^n_{ij}=A^t_{ij} \\quad \\forall (i, j) \\in \\{1, \\cdots, N\\}^2$$\n",
    "\n",
    "* the matrix resulting from the contribution of the vortices on the tangential velocity, $B^t$, is the **opposite** of the matrix resulting from the contribution of the sources on the normal velocity, $A^n$, i.e.\n",
    "\n",
    "$$B^t_{ij}=-A^n_{ij} \\quad \\forall (i, j) \\in \\{1, \\cdots, N\\}^2$$\n",
    "\n",
    "where the superscripts $n$ and $t$ denote the normal and tangential velocity components, respectively.\n",
    "\n",
    "Thus, the Kutta condition can be cast in the following form\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "(A^t_{11}+A^t_{N1}) \\\\\n",
    "\\vdots \\\\ \n",
    "(A^t_{1N}+A^t_{NN}) \\\\\n",
    "\\left( \\sum_{j=1}^N (B^t_{1j}+B^t_{Nj}) \\right)\n",
    "\\end{bmatrix}^T\n",
    "\\begin{bmatrix}\n",
    "\\sigma_1 \\\\\n",
    "\\vdots \\\\\n",
    "\\sigma_N \\\\\n",
    "\\gamma\n",
    "\\end{bmatrix} = -(b^t_1+b^t_N)\n",
    "$$\n",
    "\n",
    "Let's define a function `kutta_condition` that does the job:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def kutta_condition(A_source, B_vortex):\n",
    "    \"\"\"\n",
    "    Builds the Kutta condition array.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    A_source: 2D Numpy array of floats\n",
    "        Source contribution matrix for the normal velocity.\n",
    "    B_vortex: 2D Numpy array of floats\n",
    "        Vortex contribution matrix for the normal velocity.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    b: 1D Numpy array of floats\n",
    "        The left-hand side of the Kutta-condition equation.\n",
    "    \"\"\"\n",
    "    b = numpy.empty(A_source.shape[0] + 1, dtype=float)\n",
    "    # matrix of source contribution on tangential velocity\n",
    "    # is the same than\n",
    "    # matrix of vortex contribution on normal velocity\n",
    "    b[:-1] = B_vortex[0, :] + B_vortex[-1, :]\n",
    "    # matrix of vortex contribution on tangential velocity\n",
    "    # is the opposite of\n",
    "    # matrix of source contribution on normal velocity\n",
    "    b[-1] = - numpy.sum(A_source[0, :] + A_source[-1, :])\n",
    "    return b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have all the ingredients we need.  The function `build_singularity_matrix` will gather together the source matrix, the vortex array, and the Kutta array to build the matrix needed to solve the linear system. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def build_singularity_matrix(A_source, B_vortex):\n",
    "    \"\"\"\n",
    "    Builds the left-hand side matrix of the system\n",
    "    arising from source and vortex contributions.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    A_source: 2D Numpy array of floats\n",
    "        Source contribution matrix for the normal velocity.\n",
    "    B_vortex: 2D Numpy array of floats\n",
    "        Vortex contribution matrix for the normal velocity.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    A:  2D Numpy array of floats\n",
    "        Matrix of the linear system.\n",
    "    \"\"\"\n",
    "    A = numpy.empty((A_source.shape[0] + 1, A_source.shape[1] + 1), dtype=float)\n",
    "    # source contribution matrix\n",
    "    A[:-1, :-1] = A_source\n",
    "    # vortex contribution array\n",
    "    A[:-1, -1] = numpy.sum(B_vortex, axis=1)\n",
    "    # Kutta condition array\n",
    "    A[-1, :] = kutta_condition(A_source, B_vortex)\n",
    "    return A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On the right-hand side, we store everything that does not depend on the unknown strengths."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def build_freestream_rhs(panels, freestream):\n",
    "    \"\"\"\n",
    "    Builds the right-hand side of the system \n",
    "    arising from the freestream contribution.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    panels: 1D array of Panel objects\n",
    "        List of panels.\n",
    "    freestream: Freestream object\n",
    "        Freestream conditions.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    b: 1D Numpy array of floats\n",
    "        Freestream contribution on each panel and on the Kutta condition.\n",
    "    \"\"\"\n",
    "    b = numpy.empty(panels.size + 1, dtype=float)\n",
    "    # freestream contribution on each panel\n",
    "    for i, panel in enumerate(panels):\n",
    "        b[i] = -freestream.u_inf * numpy.cos(freestream.alpha - panel.beta)\n",
    "    # freestream contribution on the Kutta condition\n",
    "    b[-1] = -freestream.u_inf * (numpy.sin(freestream.alpha - panels[0].beta) +\n",
    "                                 numpy.sin(freestream.alpha - panels[-1].beta) )\n",
    "    return b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A = build_singularity_matrix(A_source, B_vortex)\n",
    "b = build_freestream_rhs(panels, freestream)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As in the previous lesson, we solve the linear system using the Scipy function `linalg.solve`. Then we store the results in the attribute `sigma` of each `Panel` object. We also create a variable `gamma` to store the value of the constant vortex strength."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# solve for singularity strengths\n",
    "strengths = numpy.linalg.solve(A, b)\n",
    "\n",
    "# store source strength on each panel\n",
    "for i , panel in enumerate(panels):\n",
    "    panel.sigma = strengths[i]\n",
    "    \n",
    "# store circulation density\n",
    "gamma = strengths[-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Surface pressure coefficient\n",
    "\n",
    "The pressure coefficient at the $i$-th panel center is:\n",
    "\n",
    "$$C_{p_i} = 1 - \\left(\\frac{V_{t_i}}{U_\\infty}\\right)^2$$\n",
    "\n",
    "So, we have to compute the tangential velocity at each panel center using the function `compute_tangential_velocity`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def compute_tangential_velocity(panels, freestream, gamma, A_source, B_vortex):\n",
    "    \"\"\"\n",
    "    Computes the tangential surface velocity.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    panels: 1D array of Panel objects\n",
    "        List of panels.\n",
    "    freestream: Freestream object\n",
    "        Freestream conditions.\n",
    "    gamma: float\n",
    "        Circulation density.\n",
    "    A_source: 2D Numpy array of floats\n",
    "        Source contribution matrix for the normal velocity.\n",
    "    B_vortex: 2D Numpy array of floats\n",
    "        Vortex contribution matrix for the normal velocity.\n",
    "    \"\"\"\n",
    "    A = numpy.empty((panels.size, panels.size + 1), dtype=float)\n",
    "    # matrix of source contribution on tangential velocity\n",
    "    # is the same than\n",
    "    # matrix of vortex contribution on normal velocity\n",
    "    A[:, :-1] = B_vortex\n",
    "    # matrix of vortex contribution on tangential velocity\n",
    "    # is the opposite of\n",
    "    # matrix of source contribution on normal velocity\n",
    "    A[:, -1] = -numpy.sum(A_source, axis=1)\n",
    "    # freestream contribution\n",
    "    b = freestream.u_inf * numpy.sin([freestream.alpha - panel.beta \n",
    "                                      for panel in panels])\n",
    "    \n",
    "    strengths = numpy.append([panel.sigma for panel in panels], gamma)\n",
    "    \n",
    "    tangential_velocities = numpy.dot(A, strengths) + b\n",
    "    \n",
    "    for i, panel in enumerate(panels):\n",
    "        panel.vt = tangential_velocities[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# tangential velocity at each panel center.\n",
    "compute_tangential_velocity(panels, freestream, gamma, A_source, B_vortex)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we define a function `compute_pressure_coefficient` to compute the surface pressure coefficient:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def compute_pressure_coefficient(panels, freestream):\n",
    "    \"\"\"\n",
    "    Computes the surface pressure coefficients.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    panels: 1D array of Panel objects\n",
    "        List of panels.\n",
    "    freestream: Freestream object\n",
    "        Freestream conditions.\n",
    "    \"\"\"\n",
    "    for panel in panels:\n",
    "        panel.cp = 1.0 - (panel.vt / freestream.u_inf)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# surface pressure coefficient\n",
    "compute_pressure_coefficient(panels, freestream)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Time to plot the result!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAncAAAGKCAYAAABuLmc4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3Xd4VMXi//H3BBIIoYUWeoIooAgo\nRRCkioDY9XotwYIoyrV87aCgoAgWrvXqFQuCJYoNywV+CgoBFRAQEdEIinSQGmmhBDK/P2Y3dVNJ\ntuXzep7z7O7Z2bNzDiF8mDkzY6y1iIiIiEh4iAh0BURERESk9CjciYiIiIQRhTsRERGRMKJwJyIi\nIhJGFO5EREREwojCnYiIiEgYUbgTKceMMdcbY6wx5m9jTGyu9yp63hsTgHqN8Xx3RX9/d3EYYyKM\nMc8ZY7YaYzKMMZ8Guk4lZYxJNsYk++m7vvD8+T7m471YY8zrxpidxpgDxpivjDFt/FEvkXChcCci\nADWA4YGuRAj6B/B/wASgG3B/YKsT/IwxVwHt8nnPAJ8DA4DbgcuASGCuMaax3yopEuIU7kQEYBZw\nuzGmfqAr4i/GmEqlcJiTPY/PWWsXWmtXl8Ixw5YxpibwLHB3PkUuBM4CrrHWvmet/cKzLwIFZ5Ei\nU7gTEQBv99jIggp5u0t97J9ijFmX7XWCp9vtFmPM48aYv4wx+4wx7xhjqhhjTjTGfGmM2W+M+cMY\nc10+X3myMWauMSbN0/X5qDEmx+8tY0wdY8zLxpjNxpjDxpjfjDFDc5Xxdj/3MMZ8aIz5G/i+kHMd\nYIxZaIw5aIzZY4z51BjTMtv764AxnpfHPMe/voDjWWPMOGPMSGPMJs9x5xtjTstVrp8xZqbnfNOM\nMSuNMfcYYyrkKrfOcz2vNMakeLowlxpjzvLx3T2NMV97/gwOeK79qYWcf1VjzH+MMRs813Wbp4u0\nVUGfK8RTwC/W2vfyef9CYIu1dq53h7V2D/A/4KLj+F6RckXhTkQAtgIvAkONMfGleNwHgIbAdcDD\nwBXAROATYAZwCbACmGyMae3j858CXwEXA+8CD3mOA4AxpjrwHXAeLmidhwsCLxtjbvdxvCRgLa47\ndUR+lTbGDPDUb7+nzsOAU4FvjTGNPMUuAaZ4np/p2Wbkd0yPa4GBwG3A9UAc8LUxpla2MicAXwM3\neM7nTc+5jfNxvO7APbjrcgVQAZjuaSHznst5nuPtBwYBVwPVgG+MMU0KqOuzwD+BR4BzgFuA5UD2\nY0/xFfZ98YTOa4F/FVCsNbDSx/5fgKbGmKpF+S6Rcs9aq02btnK64QKGBU4EagF/A2943qvoeW9M\ntvJj3K+NPMeZAqzL9jrB89k5ucpN8+wflG1fLHAUGJ37e4ARuT7/GrAPqOl5/RBwCDjJR7mdQMVc\n5/lsEa/LUuB37+c9+5oB6cAz2fY95ut65HNM66lTTK7rlA6MzeczxvPnMBJIBSKyvbfOsy82276O\nnu+5Otu+P4Cvcx23uqcuz2XblwwkZ3u9Mvu55lO/ScDRIpx7JC6gPZbrejyWq9xqYKqPz9/oKd8k\n0H9ntGkLhU0tdyICgLV2N/A0cG327sfj9P9yvf7N8/hltu9NBbYDvlqRPsj1eipQFdeKBu7G+++B\ntZ7RvRU9I2y/BGoDp+T6/CeFVdgYEwO0B9631h7NVs+1uFbCnoUdowAzrbUHsh1zHbAI1+rn/f4G\nxphXjDHrgSO48PcYrsWsXq7jLfRcP6+fPY9NPcc6CWgOJOW6PmnAQqBHAXVdAlxvjHnQGNMxd7ew\np/5DrLVFGdE8HIjGd+tjdgYX4nztF5EiUrgTkeyeBXYDj5bS8VJzvT5SwP7KPj6/LZ/X3q7ReriA\nkp5r+9Dzfu1cn99aeJWJxYUJX2X/wrVwllTu8/HuawRuahXcaNHzcYGuD9CJrFCU+xrtzv7CWns4\nVzlvGJxE3mt0PnmvT3a3A6/guoeXANuNMc8aY6oU8Jk8jDFNcS2PDwGVjDE1s3Ube197g+NufF9f\n7zQ9uX9uRMSHoJ5DSkT8y1q73xjzOK4Fb4KPIocAjDFR1toj2fYXFBKORxzwZ67XAJs9j7twrX7/\nl8/nV+V6XZT7w1I95XyNHK7v+c6Sistnn/d8muO6Vq+x1r7jLWCMuaCE3+et6wO4exdzO+JjH+B+\nFjyfe8BzH+Y/gCc8nynOtDkn4MLmOz7eu9eznY67n+8XoJ+PcqcAGzx1EpFCqOVORHL7Ly5s5Jlg\nFljvecwcaelphelaRnX5Z67XV+IGBnhvuv8CaIX7h3+pj21fcb/Q0236A3B59q5IT8DpCswryYl4\nDPR0+3qPmQB0wXWRAnhbxdKzlYkEEkv4fatw9+a1zuf6rCjKQay16621T+O6fQscZevDcqC3jw1c\n4OuNuy8QXKtlI2NMZte3Z9DMBZ73RKQI1HInIjlYaw8bYx4FXvXx9v8D9gCvGWNGA5Vw84+VVYvK\nTZ6uyiVAf9yN9WOstX973n8WN0r0G2PMs7gwE4MLfN2ttSWdPuMh3MjX6caY/+Lu83sEd+5Pl/Rk\ngIPALGPMBNy1ewTY6zkPgBRcgB5njDmGC3l3lfTLrLXWGHMr8JkxJgp3D+NOXGthV1wofsbXZ40x\nC3GB6mfcn29P3OTDb2YrMwm4rqD77jx/Vsk+jg+w3lqb/b3PcUH3HWPMfbhW1Adw3eRPFemkRUQt\ndyLi02TcaNEcPP9Qnw9k4ILC48B/gLm5y5aSi3DTcHyOm8bjMWBstvrswYWUmbiuwi+BNzyfK3Gd\nrJs89zzcIIYPcNO3pABnWWu3lPS4wFu40PgiLiTtAM72DGbB09V9Me7evreAl4D5uO7QErHWzsTd\nlxgDvI67Rk/hupgXFvDR+biW0yRPnf8B3GWtfT5bmQqerVRYazNwP1+zcS3InwDHgN7W2o2l9T0i\n4c5YW6QpikRE5Dh45oMbZ60dFei6iEh4U8udiIiISBgJ+nBnjJngWU5ohTHmk+wzr+cqN8AYs8q4\npYzynXleREREJJwFfbesMaYfbpb7o8aYJwGstcNzlamAm9n8HGAT7ubrq6y1v/q7viIiIiKBFPQt\nd9baWdlmiV8ENPZR7AzgD2vtn54bkqeiRaZFRESkHAr6cJfLDeRdzgjc7O7ZR1JtImsGexEREZFy\nIyjmuTPGfIXv2eBHWms/85QZiVtcPMnXIXzs89nfbIwZCgwFiI6O7tCkia/lLINfRkYGERGhls1D\nm665/+ma+5+uuf/pmvtfqF7z1atX77TW1i2sXFCEO2tt34LeN8Zch5v76Gzr+ybBTeRcdLwx4HMu\nKmvtq3gmZ+3YsaNdunRpieocaMnJyfTq1SvQ1ShXdM39T9fc/3TN/U/X3P9C9ZobY9YXXioEumWN\nMQNwk5NeaK1Ny6fYEuAkY0wzzyzsV6KlakRERKQcCvpwh5vJvRow2xiz3BgzEcAY09AYMxPAM+Di\nNtzM6ynAB9baXwJVYREREZFACYpu2YJYa0/MZ/8WYGC21zNxSxCJiIiIlFuh0HInIiIiIkWkcCci\nIiISRhTuRERERMJI0N9zJyIiUlr27t3L9u3bSU9PD3RVMtWoUYOUlJRAV6NcCbZrHhkZSb169ahe\nvXqpHE/hTkREyoW9e/eybds2GjVqRHR0NMb4mv/e//bt20e1atUCXY1yJZiuubWWgwcPsnnzZoBS\nCXjqlhURkXJh+/btNGrUiCpVqgRNsBMxxlClShUaNWrE9u3bS+WYCnciIlIupKenEx0dHehqiPgU\nHR1darcLKNyJiEi5oRY7CVal+bOpcCciIiISRhTuRERERMKIwp2IiIiUGxkZGQwfPpwGDRoQERHB\nxRdfHOgqlTpNhSIiIiLlxkcffcTLL7/M008/zZlnnknt2rUDXaVSp5Y7ERGRkkhKgoQEiIhwj0lJ\nga5R0LPWcuTIkYB89+HDhwEyJy++8847OfPMM2nRokVA6lOWFO5ERESKKykJhg6F9evBWvc4dKhf\nA971119PQkJCnv29evWiV69ema+Tk5MxxvDxxx9z/fXXExsbS/Xq1UlMTGTXrl05PmuMYeTIkYwb\nN47GjRsTHR1Njx49WL58eZ7vmTZtGl26dKFKlSrUrFmTyy+/nA0bNuQok5CQwKBBg3jjjTdo1aoV\nUVFRzJgxI99zev755zn55JOJjo4mNjaWjh078sknn+Q43vXXX5/nc8YYxowZk/l6zJgxGGNYuXIl\n/fv3p2rVqvzzn/8kISEhs1yFChUwxjBlyhQARo8eTfv27alRowZ16tShT58+LFq0KM937dixg3/9\n6180adKESpUq0aRJE6655prM8Ajw008/ceGFFxIbG0t0dDTdunXjm2++yfe8S5u6ZUVEpPwqzalR\n0tJg0CC3FUM1cAGxjN1555307duX9957j99//50HH3yQLVu2MHfu3Bzl3nrrLZo2bcqLL77I4cOH\nefjhhzn77LP5/fffqVWrFgATJ05k2LBhDB48mIcffph9+/YxZswYevbsyYoVK3Ks/jB37lyWL1/O\n6NGjqVevns9ACpCUlMQ999zDww8/TPfu3Tl48CArVqxg9+7dJT7niy66iCFDhjB8+HAiIiKoUaMG\nL7zwAlOmTGHhwoUANG/eHIDNmzdz11130bhxYw4cOMA777xDjx49WLp0KW3btgUgNTWVrl27snv3\nbkaNGkXbtm3Zvn07n332GUeOHKFSpUosW7aM7t27c/rpp/Paa69RpUoVJk6cSN++fVmwYAEdOnQo\n8fkUlcKdiIhIOdC6dWsmT54MwIABA6hVqxaDBg3i66+/5uyzz84sd/DgQWbNmkVMTAwAnTt35qST\nTuLZZ59l7Nix7N+/n+HDhzN48GDeeOONzM917tyZFi1aMGnSJO68887M/ampqfzwww/Ur1+/wPot\nXLiQtm3b8vDDD2fuGzhw4HGd8x133MH//d//5djXqFEjALp06ZJj/+uvv575/NixYwwYMIDWrVsz\nadIknn/+eQCeffZZ/vzzT5YuXcrpp5+eWf6qq67KfH7ffffRtGlT5syZQ1RUFAD9+/fn1FNPZezY\nsXz66afHdU5FoW5ZEREpv6wt2RYf7/t48fHFPta+vXv9cqr//Oc/c7y+/PLLiYiIyGzB8ho4cGBm\nsAPXFdqlS5fMcgsXLmTv3r0kJiZy9OjRzK1x48a0atWK+fPn5zhely5dCg12AJ06dWL58uXcfvvt\nfPXVV6SlpZX0VDNdcsklRS771Vdf0bt3b2rXrk3FihWJjIxk9erVrFq1KrPMrFmz6NSpU45gl93B\ngweZN29e5rX1XhtrLX379s1zbcqKwp2IiEhxjRsHVark3FelitsfpOLi4nK8joqKIjY2NnPB+vzK\nefd5y3nXP+3bty+RkZE5tp9//jnPfXwNGjQoUv2uvfZaXn75Zb7//nv69+9PrVq1uPTSS1m3bl1R\nTzGPon73smXLGDhwIFWrVmXSpEksWrSIJUuW0K5dOw4dOpRZbteuXTRu3Djf4+zevZtjx44xduzY\nPNfmxRdfJDU1lYyMjBKfT1GpW1ZERKS4EhPd48iRsGEDNG3qgp13vx9UrlzZ58jTXbt2+ZzeY9u2\nbTleHzlyhNTU1MxuyvzKefd5y3mPPWXKFFq3bp2nbPb77aDoy2oZY7j55pu5+eabSU1NZdasWdxz\nzz1cccUVfP/994Dvcy7onryifvfHH39MxYoVmTZtGpGRkZn7U1NTqVmzZubrOnXq5AnD2dWsWZOI\niAhuvfVWrr32Wp9lIiLKvl1N4U5ERKQkEhP9GuZyi4+PZ9u2bezcuZM6deoAsGbNGlatWkXXrl3z\nlP/ggw+44YYbMl9/+OGHZGRkcOaZZ+YoN3PmTA4cOJDZNbtu3ToWLVrEiBEjAOjatSvVqlXjjz/+\n4LrrriuTc4uNjc0Mda+88krm/vj4eFauXJmj7PTp04/7+9LS0jJHz3rNmTOHDRs20KxZs8x9/fr1\n47HHHuOnn36iXbt2eY4TExND9+7d+emnn2jfvr1fgpwvCnciIiIh6PLLL+ehhx4iMTGRu+++m507\nd/L4449nBr3cfvnlFwYPHsyVV17J6tWrGTlyJD179swxmAIgOjqafv36cd9993H48GFGjx5N9erV\nueuuuwCoXr06EyZM4NZbb2XHjh2ce+651KhRg82bNzNv3jx69erF1VdfXezzGTp0KNWqVePMM8+k\nXr16rF69mrfffpt+/fpllrnyyiu54YYbuOuuuzj//PP56aefMqcyOR4DBgzgueee4/rrr2fw4MGs\nXr2asWPH5mnVvOuuu3j33Xfp27cvo0aNok2bNuzcuZPPPvuMiRMnUq1aNZ555hl69OhB//79GTJk\nCA0aNGDnzp0sW7aMY8eO8cQTTxx3fQujcCciIhKCTjzxRD766CNGjRrFxRdfTIsWLXjmmWcYP368\nz/LPP/88n3/+OVdccQXHjh3jggsu4IUXXshT7tprryUmJobbbruNnTt30qlTJ6ZOnZo5DQrAzTff\nTJMmTZgwYQLvvvsu6enpNGrUiB49enDaaaeV6Hy6devG5MmTefvtt9mzZw8NGzZk0KBBPPLII5ll\nrrvuOjZu3MikSZN45ZVX6N69O5988gknnnhiib7Tq3///rzwwgs888wzfPzxx5x66qm89dZbPPbY\nYznK1axZk++++45Ro0bxxBNPsGvXLuLi4ujTp0/myNj27duzZMkSHnnkEe644w727NlD3bp1ad++\nPbfccstx1bOojPXD3DrBqmPHjnbp0qWBrkaJJCcn55ikUsqerrn/6Zr7Xzhf85SUFE4++eRAVyOP\nffv25blPrTQlJyfTu3dvZs+eTd++fQss653EOHeoCTdlfc1LqrCfUWPMD9bajoUdR6NlRURERMKI\nwp2IiIhIGNE9dyIiImGsV69eFPUWrPJ8q1Y4UcudiIiISBhRuBMREREJIwp3IiIiImFE4U5EREQk\njCjciYiIiIQRhTsRERGRMKJwJyIiIhJGFO5ERERC0JgxYzDGBLoaIeG3336jT58+VK9eHWMM06dP\nD3SVypQmMRYREZGwdvfdd/Pnn3/ywQcfULNmTRo2bBjoKpUptdyJiIiUQFISJCRARIR7TEoKdI2C\n0+HDhwPyvdZajhw5AkBKSgo9evRgwIABdOnShdjY2IDUyV8U7kRERIopKQmGDoX168Fa9zh0aOAD\n3t69e7ntttto2LAhlSpVomXLljz77LOZy4odO3aMmjVr8thjj2V+5ueff8YYw1lnnZXjWI0bN+b+\n++/PfJ2Wlsbw4cNp1qwZUVFRNGvWjHHjxpGRkZFZJjk5GWMM06ZN46abbqJu3brExcXlW9+//vqL\n6667LrO+DRo04Pzzz2f79u05jpecnJzjc1OmTMEYw7p16zL3JSQkMGjQIN544w1atWpFVFQUM2bM\nyCz39ttvY4zJ7Mr+448/uOaaa2jWrBnR0dGccMIJDBs2jNTU1Dz1nDdvHueccw41atQgJiaGdu3a\nMWnSpBxlXnvtNdq1a0flypWpU6cOQ4YMYffu3fmee1lSt6yIiEgxjRwJaWk596Wluf2JiYGpU0ZG\nBueddx7Lli3j0UcfpU2bNsyYMYO7776bHTt2MH78eCpUqECPHj2YM2cOo0aNAmDOnDlER0ezePFi\nDhw4QExMDKtWrWLz5s307t0bgKNHj9K/f39+/fVXHnroIdq0acOiRYsYO3Ysu3fv5umnn85Rl9tv\nv51zzz2Xt99+m0OHDuVb52uuuYb169czYcIEmjRpwrZt2/j6669Jy31xi2ju3LksX76c0aNHU69e\nPerUqcPChQu58MIL6dSpEw899FBm2S1bttC4cWOee+45YmNj+fPPPxk/fjwDBw5k4cKFmeU+++wz\nLrvsMrp168Yrr7xCnTp1+OWXX1i/fn1mmREjRvD0009zxx13MGHCBDZv3syoUaNYuXIlCxYsoEKF\nCiU6nxKz1pbbrUOHDjZUzZ07N9BVKHd0zf1P19z/wvma//rrrz73u7a3wG4lMXr0aEu2D//vf/+z\ngJ08eXKOckOGDLFRUVF2x44d1lprn3nmGVu5cmV76NAha621F110kb3llltsTEyM/eKLL6y11r78\n8su2YsWKdt++fdZaa9966y0L2Hnz5uU49mOPPWYjIyPttm3brLXu5wewF198cZHOISYmxj7//PP5\nvu89Xu6fy8mTJ1vArl27NnNffHy8jY6Otlu3bs1znEaNGtnrrrsu8/XevXvzlElPT7fffPONBeyy\nZcustdZmZGTY+Ph426FDB3vs2DGfdVy7dq2NiIiwjzzySI793377rQXsJ598ku/55Zbfz6gXsNQW\nId+oW1ZERMq1ksSx+Hjfx4qPL/6x9u7dVyrnMX/+fCIiIrjqqqty7B80aBBHjhzJbI3q3bs3hw4d\nYsGCBWRkZDBv3jz69+9Pt27dmDNnDuBa8zp16kTVqlUB+OKLL4iPj6dr164cPXo0c+vXrx/p6eks\nWrQox3decsklRapzp06dmDBhAs8//zw///wzLr+UXJcuXahfv36Ryh45coTx48fTqlUroqOjiYyM\npHv37gCsWrUq83H9+vXceOONRET4jkyzZ88mIyODxMTEHNemc+fOVK9enfnz5x/XOZWEwp2IiEgx\njRsHVark3FelitsfKLt376ZWrVpUqlQpx35v2PHe/9WuXTtq167N3Llz+fHHH9m7dy89e/akd+/e\nzJ07F2stycnJmV2yANu3b2f9+vVERkbm2M444wwAdu3aleM7GzRoUKQ6v//++1x44YU89dRTtG3b\nlkaNGvHoo4/muI+vOIr6vQAPPPAAY8aMYdCgQcyYMYPFixczbdo0gMyuZO95NW7cON/jeO8PPPHE\nE/Ncn7179+a5Nv6ge+7CSVKSu+FjwwZo2tT9lgnUzR8iImHM+6s1mH7l1qpVi927d3PkyBGioqIy\n9//1118A1K5dGwBjDD179mTOnDlUq1aN0047jdjYWPr06cOoUaP47rvv2LFjR45wV7t2bZo1a8YH\nH3zg87sTEhJyvC7q/Hv16tXjpZde4qWXXmLVqlW8+eabjB49mrp16zJs2DAqV64MkDnq1Su/wFSc\nef+mTp3Ktddem3nvIcD+/ftzlKlTpw4Amzdvzvc43us6a9Ysn6Nwve/7k8JduPAO3fLehOodugUK\neCIiZSAxMbh+vfbs2ZMJEybw4YcfkpitYklJSURFRdGlS5fMfb179+buu++mQoUK9OnTB4AOHToQ\nExPDmDFjiIqKolu3bpnlBwwYwMcff0zVqlVp1apVmdS/ZcuWjB8/nokTJ7Jy5UoA4j393ytXrqRf\nv36ZZWfOnHnc35eWlkZkZGSOfZMnT87xukWLFiQkJPD6668zdOhQn+HxnHPOISIigg0bNnDOOecc\nd71Kg8JduAjGoVsiIuI35557LmeddRa33HILO3bsoHXr1sycOZPXX3+dBx54ILMVCqBPnz6kp6cz\nf/58hg8fDpA5knb69On06NGD6OjozPKJiYlMnjyZs88+m3vuuYd27dpx5MgR1qxZw+eff86nn35K\nldz91IXYs2cPffv2JTExkVatWhEZGclnn31GampqZpBr0KABPXv25PHHH6dOnTrUq1ePd955hzVr\n1hz39RowYABvvvkmbdq04cQTT2TatGksWLAgRxljDM899xyXXnopffr04ZZbbqFu3bqkpKSwfft2\nHnnkEZo3b87w4cO57bbbWLVqFT179qRy5cps3LiR2bNnc+ONN+ZoBfUHhbtwsWFD8faLiEhYiYiI\nYMaMGTz44IM8+eST7Nq1i4SEBJ555hnuvPPOHGVPOeUU4uLi2LVrV+YgAnChb/r06XnCSGRkJF9+\n+SVPPPEEr776KmvXriUmJobmzZtz3nnn5egGLqrKlSvTvn17XnvtNdavX09ERAQtW7YkKSmJiy66\nKLPcO++8w7Bhw7jjjjuoXLkyN9xwA6NGjeKmm24q9ndm95///AdrLSNHjgRg4MCBvPfee5n3EXpd\ndNFFzJ49m7FjxzJkyBAAmjdvnuOajh8/npNPPjmzi9kYQ5MmTTj77LM56aSTjqueJWGOd2RKKOvY\nsaNdunRpoKtRIsnJyfTq1StrR0KC64rNLT4esk3yKCWX55pLmdM1979wvuYpKSmcfPLJga5GHvv2\n7aNatWqBrka5EqzXvLCfUWPMD9bajoUdJ+hHyxpjLjfG/GKMyTDG5HtCxph1xpifjTHLjTGhmdiO\nx7hxkPt/TpUrB3boloiIiPhd0Ic7YCVwKVCUiWJ6W2tPK0qqDTuJiXDhhTn3nXACXH11YOojIiIi\nARH04c5am2KtXRXoeoSERo3c46hRULMm/PorzJoV2DqJiIiIX4XTgAoLzDLGWOAVa+2rvgoZY4YC\nQwHi4uLyLEYcKvbv35+n7ievXEkckHLsGFFXXEHzV15h/623svTVVyGfmbWl6Hxdcylbuub+F87X\nvEaNGuzbVzqrQZSmY8eOBWW9wlmwXvNDhw6Vyt+/oAh3xpivAF/rhYy01n5WxMN0s9ZuMcbUA2Yb\nY36z1ubpyvWEvlfBDagI1RuHfd70PH48ACf36AG9esHMmVRds4ZeW7bAoEF+r2O4CecbzYOVrrn/\nhfM1T0lJoWrVqsWa6NYfgvXm/nAWjNfcWkvlypU5/fTTj/tYQdGcY63ta6091cdW1GCHtXaL53E7\n8AlwRsGfCEM7drjHunXdYIpHH3WvR40Cz1IqIiLlVWRkJAcPHgx0NUR8OnjwYJ5JlUsqKMLd8TLG\nxBhjqnmfA/1wAzHKl+zhDuCaa+DUU90UKS+/HLh6iYgEgXr16rF582bS0tKOe4F6kdJirSUtLY3N\nmzdTr169UjlmUHTLFsQYcwnwH6AuMMMYs9xa298Y0xB43Vo7EIgDPvE0tVcE3rXWfhGwSgeCtVnh\nzjsLeYUK8MQTcP758NhjMHiwG2ghIlIOVa9eHYAtW7aQnp4e4NpkOXToUOYaquIfwXbNIyMjiYuL\ny/wZPV5BH+6stZ/gullz798CDPQ8/xNo5+eqBZd9++DIEahSxW1eAwdCz54wbx48+SQ8/njg6igi\nEmDVq1cvtX9AS0tycnKp3Gfg1kqtAAAgAElEQVQlRRfu1zwsumWFvF2yXsa4UAfw3HOwebN/6yUi\nIiJ+pXAXLvILdwCdO8M//uEGVYwe7d96iYiIiF8p3IWLgsIduGXIKlSAyZPd5MYiIiISlhTuwkVh\n4a5FCxg6FDIy4IEH/FcvERER8SuFu3BRWLgDePhhiImBzz+Hb7/1T71ERETErxTuwsXOne6xoHBX\nvz7cc497fv/9bvoUERERCSsKd+GiKC13APfe68osXAifFXkBEBEREQkRCnfhoqjhrlo11z0LcPnl\nEBEBCQmQlFSm1RMRERH/CPpJjKWIihruwAU8Y+DoUfd6/Xo32AIgMbFs6iciIiJ+oZa7cFGccDd6\ndN777dLSYOTI0q+XiIiI+JXCXbgoTrjbsKF4+0VERCRkKNyFg4MH4cABiIpyXa6Fadq0ePtFREQk\nZCjchQNvq12dOu5eusKMGwdVquTdf8cdpVsvERER8TuFu3BQnC5ZcIMmXn0V4uNdGPQGvfffh/T0\nsqmjiIiI+IXCXTgobrgDF/DWrXPLkW3aBE2awOLF8OijZVJFERER8Q+Fu3BQknCXXWwsvP22a8Ub\nPx6++ab06iYiIiJ+pXAXDo433AH07AkjRriWvEGD4O+/S6duIiIi4lcKd+GgNMIdwCOPQMeObkqU\nYcO09qyIiEgIUrgLBzt3usfjDXeRkfDuu26AxdSp8M47x183ERER8SuFu3BQWi13ACedBC+84J7f\neiv8+efxH1NERET8RuEuHJRmuAO44Qa47DLYt8/df+ddg1ZERESCnsJdOCjtcGeMmwevUSNYuBAe\ne6x0jisiIiJlTuEuHJR2uAOoVStrepSxY2HBgtI7toiIiJQZhbtQl57upi2pUMHNV1eaeveG++5z\n06MkJsLevaV7fBERESl1CnehzjtStnZtiCiDP86xY6F9e7eaxa23lv7xRUREpFQp3IU6b5dsnTpl\nc/yoKDc9SnS0mxrl3XfL5ntERESkVCjchbqyuN8ut5Yt4bnn3PNhw1wrnoiIiAQlhbtQ549wB3DT\nTXDxxe6+u2uu0fQoIiIiQUrhLtT5K9wZA6+9Bg0awLffwhNPlO33iYiISIko3IU6f4U7cPf1vfmm\nez5mDHz/fdl/p4iIiBSLwl2o82e4AzjnHLjnHjh2DC64AJo2daN0ExIgKck/dRAREZF8KdyFOu9U\nKP4KdwDjxrlQt2MHbNwI1sL69TB0qAKeiIhIgCnchTp/t9wBVKrkJk/OLS0NRo70Xz1EREQkD4W7\nUBeIcAfw11++92/Y4N96iIiISA4Kd6EuUOGuaVPf+5s08W89REREJAeFu1CWkQG7drnntWv797vH\njYMqVfLuj4lxa92KiIhIQCjchbLdu13Ai42FyEj/fndiIrz6KsTHuznw6teH6tUhJQW6dXMDLERE\nRMTvFO5CWaC6ZL0SE91SZBkZsHUrrFgBp5wCv/4KXbrAsmWBqZeIiEg5pnAXygId7nKLj4fvvoPe\nvd2Aix49YMaMQNdKRESkXFG4C2XecFenTmDrkV3NmvDFF2792QMH4MILYeLEQNdKRESk3FC4C2XB\n1nLnFRXllil76CHXZTtsGAwf7p6LiIhImVK4C2XBGu7ADbJ49FGYNAkqVoSnnoKrr4ZDhwJdMxER\nkbCmcBfKgjnced1wA8ycCdWqwfvvu7VpvdO3iIiISKlTuAtloRDuwAW6b7+FRo3cY9eu8Oefga6V\niIhIWFK4C2U7d7rHYA93AG3bwvffQ7t2sHq1myrl++8DXSsREZGwE/ThzhhzuTHmF2NMhjGmYwHl\nBhhjVhlj/jDGjPBnHQMmVFruvBo1gvnzoX9/V/feveHTTwNdKxERkbAS9OEOWAlcCszPr4AxpgLw\nEnAucApwlTHmFP9UL4BCLdyBW8Xif/+DIUPg4EG49FJ4/vlA10pERCRsBH24s9amWGtXFVLsDOAP\na+2f1tojwFTgorKvXQBZG5rhDtxSaa+9Bo895s7jzjvdduxYoGsmIiIS8ioGugKlpBGwMdvrTUBn\nXwWNMUOBoQBxcXEkJyeXeeXKwqHt2yE9naPR0Xy7aFGgq1My3bpR78EHafXUU0Q8/zw7li0j5cEH\nyahcOdA182n//v0h+/MSqnTN/U/X3P90zf0v3K95UIQ7Y8xXQH0fb4201n5WlEP42Gd9FbTWvgq8\nCtCxY0fbq1evolYzqHyflARAxfr1CdVzAKBXLzea9pJLqPvNN9S97TY4fBg2b4amTWHcOLeGbRBI\nTk4O7WsdgnTN/U/X3P90zf0v3K95UIQ7a23f4zzEJqBJtteNgS3HecygFvn33+5JqHXJ+tKrl1uT\ntkcPWLMma//69TB0qHseJAFPREQk2AX9PXdFtAQ4yRjTzBgTBVwJfB7gOpWpsAp3AKecAtHRefen\npcHIkf6vj4iISIgK+nBnjLnEGLMJOBOYYYz50rO/oTFmJoC19ihwG/AlkAJ8YK39JVB19ofMcFen\nTmArUpo2b/a9f8MGN/BCREREChX04c5a+4m1trG1tpK1Ns5a29+zf4u1dmC2cjOttS2stc2tteMC\nV2P/iNqzxz0Jl5Y7cPfY+WItDBwIqwobNC0iIiJBH+7Et7DrlgU3eKJKlZz7IiNdd+0XX0CbNnD/\n/bBvX2DqJyIiEgIU7kJUWIa7xER49VWIjwdj3OPkybBuHdx4Ixw9ChMmQIsW8Pbb6qoVERHxQeEu\nFCUlUXfePPd8+HDwTIsSFhITXZjLyHCPiYlQr56b9Pj776FzZ/jrL7j2WjjrLFi2LNA1FhERCSoK\nd6EmKQmGDqXCkSPu9Y4dbrqQcAp4+enUCRYsgClTIC7OPe/YEW6+GXbuDHTtREREgoLCXagZOdJN\nD5JdeZouJCICrrvODa645x6oUMF15bZoAS+95LpuRUREyjGFu1CzYUPx9oerGjXg3/+GFSvcChep\nqXDbbdChA8yfH+jaiYiIBIzCXajJb7qQ/PaHu5NPhi+/hE8+gYQEF/Z69oSrroJNmwJdOxEREb9T\nuAs1vqYLqVLF7S+vjIGLL4Zff4VHHoHKlWHqVGjZEh5/3K1VKyIiUk4o3IUaz3Qhh+LisqYLefVV\nrb0Kbj68hx+G336Dyy5z9yI++CC0bg3Tpwe6diIiIn6hcBeKEhNZNHVqzulCJEt8PHz0EXz1lVuz\nds0auOACOO88WL060LUTEREpUwp3Er7OPhuWL4dnn4Xq1WHmTDj1VBgxAiZNcvfoRUS4x/IwlYyI\niJQLCncS3iIj4c474fff4YYbID0dnnwSbroJ1q93q1ysX19+5goUEZGwp3An5UO9eq617vvvISoq\n79Jl5WmuQBERCWsKd1K+nHGGa73zZf162LLFv/UREREpZQp3Uv4UNCdgfLwboLJkif/qIyIiUooU\n7qT88TVXYKVKrlUvIwPefdc979YNPvhAS5qJiEhIUbiT8sczVyDx8VlzBXrvx/vzT7jvPqhZExYs\ngCuugBNOgKeeouLevYGuuYiISKEU7qR8Skx0cwTmniswPh6eego2boSXXnKrXGzcCMOHc+YVV8Cw\nYZCSEsiai4iIFEjhTsSXqlXhX/9yS5rNnAn9+lHh0CGYONFNjHzuuW5N29yjbkVERAJM4U6kIBER\nmUFu8eTJcPPNbpmzL76AAQNc0Js4EQ4cCHRNRUREAIU7kSJLS0hwQW7jRnj8cWjUyK1jO2wYNGni\nVr7YuDHQ1RQRkXJO4U6kuGrXdkFu7VqYOhW6dIHUVLfyRbNmbhDGwoXqshURkYBQuBMpqcjIrCC3\naBFcdZUbffvBB9C1K3Tu7KZVOXIk0DUVEZFyROFOpDR4g9zatfDAA1CrlpsIOTHRteaNH++6dBMS\n3H18CQlay1ZERMqEwp1IaWrc2AW5jRvdXHqnnOKWNBs50t2bt369665dvx6GDlXAExGRUqdwJ1IW\nqlSBm26ClSth1iyoXDlvmbQ0uP9+/9dNRETCmsKdSFkyBs45Bw4f9v3+li3u/rz//Ae2bfNv3URE\nJCwp3In4Q9Omvvcb4wZk3HEHNGwI/frBlCmwZ49fqyciIuFD4U7EH8aNc1212VWpAq+/7gZiXHAB\nVKgAs2fD4MEQFweXXQYffwwHDwamziIiEpIU7kT8ITHRDbCIj3etdfHx7vUNN7gpVD7/HP76y+3r\n1ctNnzJtGvzjHy7oXX+9u3fv6NFAn4mIiAQ5hTsRf0lMhHXrICPDPSYm5ny/Vi03CGPuXDfa9umn\noUMH2LcP3nwT+vd3q2LcfjssWKBJkkVExCeFO5Fg1KgR3H03LF0Kq1bBmDHQsiVs3w4vvgjdusEJ\nJ7g59X7+OdC1FRGRIKJwJxLsWrSA0aMhJQV++AHuvdeFv3Xr4IknoG1baNPGza+3dm2gaysiIgGm\ncCcSKoyB9u1hwgTYsAGSk+Hmm1137sqVbqLkE06AM8/U1CoiIuWYwp1IKIqIgJ493ZJmW7fC9Olw\n9dUQE+PWufU1tUpSkpY/ExEpByoGugIicpyiouC889x24AD8739uepUvvnBTq8yeDTfe6MoeO+Ye\nvcufQd6BHSIiEtLUcicSTmJi4Morc06t0ru3C3XeYOeVlgYjRgSmniIiUmYU7kTClXdqlTlz3P16\nvmzaBGed5QZm/PKLplcREQkDCnci5UF+y58BfPedm1Ll1FOheXN3v96sWfmvhysiIkFN4U6kPCho\n+bOPPnIrYNSt66ZS+c9/3ITJdeq4FTKmTHHz64mISEjQgAqR8sA7aGLkSDeNStOmLvB59192mbsn\nb8kSNyBj+nRYscKtbfvxx65bt3Nntwbu+ee7efXy6+oVEZGAOu5wZ4ypBQwF6gBLgQ+ttccK/pSI\n+F1iYsEjYytUgC5d3DZunBtRO3262+bMcVOsLFrkAmLTpi7kXXCBWwu3cmW/nYaIiBSsNLplPwE6\nAHuBa4EfjDH1S+G4IhJI8fFw663w//4f7NoFn3wCN9wAcXGu9e+//4Vzz3Xdt5dcAm+84UboiohI\nQJVGt2wda21P7wtjzGXARODiUji2iASDqlXh4ovdlpHh1rydPt114S5fDp9+6jaAM87IatVr107d\ntyIiflYaLXe7jDFVvS+stR8DJ5TCcUUkGEVEuAD36KPw44+wcSO8/LKbRLlyZVi8GB5+GE4/3XXf\nDhsGM2fCwYNZx9BqGSIiZaY0Wu4OATONMf9nrf3RGHMiUGqLWnru6XsfSADWAf+01qb6KHcM+Nnz\ncoO19sLSqoOIFKBxY7jlFrcdOABff511r96mTW6JtIkT3ejcvn3dqNx3380Ke1otQ0SkVJVGy923\nwN/ADGPMQWAVsNIY09cYU7MUjj8C+NpaexLwtee1Lwettad5NgU7kUCIiYELL3QrY2za5LpvR4+G\nDh3cihiffw6TJuVsxQP33siRgamziEiYKVK4M8ZE5/eetfZRa+2F1tqGwInAZcAB4D5gdSnU8SLg\nTc/zN9G9fCKhISLChboxY1zI27TJhb78rF8PTz7ppmPJvVSaiIgUWaHhzhjTB9hrjLmysLLW2s3W\n2k+ttaOstf2ttfVKoY5x1tqtnuNvBfI7ZmVjzFJjzCJjjAKgSLBp1MgthxYfn3+ZESPc/Xy1a8PF\nF9Po44+1LJqISDEZW8gvTWPMx0Bda22PAsp0BFoCn1tr9xW7EsZ8BfiaPmUk8Ka1tma2sqnW2lgf\nx2hord1ijDkBmAOcba1d46PcUNy8fMTFxXWYOnVqcasbFPbv30/VqlULLyilRte8dNT76ita/vvf\nVMi2vNmxSpXYeu65RKSnE/vjj0Rv2ZLjM0diY0k9/XT+bt+e1NNP51DDhv6udrmhn3P/0zX3v1C9\n5r179/7BWtuxsHJFCXd/AcOttW8WUCYaWAuMt9a+UNzKFvL9q4Be1tqtxpgGQLK1tmUhn5kCTLfW\nflRQuY4dO9qlS5eWXmX9KDk5mV69egW6GuWKrnkpSkrKf7UMcF20c+bw13vvUX/lSti6NefnExKg\nT5+srUEDv1Y/nOnn3P90zf0vVK+5MaZI4a4o99zFAn8WVMBaexB3P9z5RatesXwOXOd5fh3wWe4C\nxphYY0wlz/M6QDfg1zKoi4iUhsREWLfOzZm3bl3eUbLx8TB4ML89+CBs3gwpKfDSS3DppRAb6z7z\nxhswaBA0bAinnAK33eYmWt69OwAnJCISPIoyFcpOIK4I5b4FCr0vrwSeAD4wxgwBNgCXQ2ZX8C3W\n2huBk4FXjDEZuMD6hLVW4U4kHBgDrVq57V//coHwp5/clCtz5sD8+S78eQOgMW6OvbPPdq163bu7\nUbwiIuVEUcLdYtwI2AK7OIE9FC0EFou1dhdwto/9S4EbPc8XAG1K+7tFJAhFRLjwdvrpcO+9kJ7u\nJk6eM8cFvoULYdkyt02YAJGR0LmzC3pnn+2eV6oU6LMQESkzRemWnQRcboy5tJByJ+DWlxUR8Z/I\nSOjWDR56CJKTITUVZs/OGnl77Bh8+61bUaNnT9et269fzmlXtGKGiISRQlvurLXTjTHvAVONMU8A\nE3KPiDXGRAF3At+VTTVFRIrIuxJG377u9Z49MG9eVsveypUu/M2e7d6PjoYjR7Lm1tOKGSIS4oq6\n/Nj1wEFgFHCXMWYasBK3zFhj3ECHZnimGBERCRo1arhVMy70LFyzbZtr4fPes7cmz4xJbsWMW2+F\n+vVdN24ITpkgIuVXkcKdtfYYMNQY8z7wAJBIzi7drcDl1trFpV9FEZFSFBcHV1zhNnBdsb6mhNqz\nx7X+VagA7dpB166u+7dbN2jSxL91FhEphqK23AFgrf0a+NqzZmxboAau9W6ZtfZoGdRPRKRsNW3q\numJzq1bNjdD98cesARovvujea9w4K+h17erCX8Vi/ToVESkzJfptZK39G5hfynUREfG/cePcPXZp\naVn7qlSBl19299wdOOAGXnz3HSxY4LZNm+D9990GbqqVzp2zWve6dIGaNX1/n4hIGdN/NUWkfPMO\nmshvxYyYGOjVy23g5tlLSXEh77vv3PbHH+7+vTlzXBljoHXrnK17J5zg9ouIlDGFOxGRxMSij4yN\niHDBrXVruOkmt2/79qywt2ABLF3qRuWuXAmvvOLKxMVlBb1u3aB9e4iKKpvzEZFyTeFOROR41asH\nF1/sNoBDh+CHH7LC3nffuVG606a5DaByZejUKSvsde0KtWu79wpbe1dEpAAKdyIipa1y5awuWXCj\ncX//PWdXbkoKfPON27xatnTTryxc6ObeA827JyLFpnAnIlLWjIEWLdx2/fVu3+7dLsR5W/cWL4ZV\nq9yWW1oa3HWXm5olrtRXeRSRMKNwJyISCLVqwXnnuQ1cS93y5W7UrS87drhWvaZNXXfuGWe4xw4d\noHp1/9VbRIKewp2ISDCIinKBLT7e97x7lSq5dXQ3bHDbxx+7/ca4+fi8Ye+MM6BtW1deRMolhTsR\nkWCS37x7r74KV17pum0XL3Zz7y1eDD/95O7fS0mBN9905aOi3MTK2Vv4WrVyI31FJOwp3ImIBJPC\n5t075RS3ee/dO3wYVqzIGfh++809X7IE/vtfV65aNdeFm72Fr0kTzb0nEoYU7kREgk1x5t2rVMmF\ntU6dsvbt3eumYvGGvSVLXFBMTnabV716OcNe9mOISMhSuBMRCTfVq0Pv3m7z2rYtZ9hbvNhNvjx9\nuts8OjdsCN27Z4W99u3h0081755ICFG4ExEpD+Li4Pzz3QZu7r21a3OGvWXLiN6yJee6ueC6bq11\nzzXvnkjQU7gTESmPjHHr3Z5wghuoAXD0KEvefJNOkBX6fvwxK9h5paXB4MGui7dDB7dphK5I0FC4\nExERp2JFDjRvDr16wZAhbl9ERN5wB5CeDq+/7jbPZzn11Kyw5w18lSv7rfoi4ijciYhI/po29T3v\nXlwcjBjhBm788IMbobt8udsmTXJlKlaE1q1zBr527RT4RMqYwp2IiOQvv3n3nn465z13+/e7YOcN\ne97A99NPbnvjDVeuQgUX+Dp2zNnCFx3t3/MSCWMKdyIikr/C5t3zqloVzjrLbV4HDuQNfCkpbl6+\nFSvyBr7cLXy+Al9SkkbuihRC4U5ERApWnHn3souJgW7d3OZ14IBrycse+H79NSvwTZ7sylWo4CZr\nzh74fvsNbrstqxVRI3dFfFK4ExER/4mJga5d3eaVlpYz8C1d6gLfzz+7bcqU/I+XluZa8hTuRDIp\n3ImISGBVqQJnnuk2r9yB74cfXNDzZf16uPbarBa+005z3cQi5ZTCnYiIBB9fga9pU9i40Xf5t992\nG7g5/Fq1cqtreAPf6ae79XVFygGFOxERCQ2PP+575O7w4VC3rmvdW7bMtfClpLgtKcmVMwZatMgb\n+GrUCMy5iJQhhTsREQkNRR25e/iwC3i5u3RXrXLbe+9llT3pJBf0vKGvfXuoWdN/5yRSBhTuREQk\ndBRl5G6lSm4evY4ds/YdOQIrV+YMfCtWwO+/u23q1KyyzZvnDXy1arn3NBWLhACFOxERCX9RUS6k\ntW8PN93k9qWnwy+/5Ax8P/0Ea9a47YMPsj7frJnr+v3xR/c50FQsErQU7kREpHyKjHQja087LWst\n3fR0d6+eN+wtW+YmYl671m25paXB7be7Vrz27d1ULyIBpnAnIiLiFRnplkNr2xYGD3b7jh51Eyi3\naeP7M6mp0KMHRES4lTbOOAM6dXKPp57qjiniRwp3IiIiBalY0YW0+HjXFZtb1apuYMaKFVkTL0+a\n5N6rXNmNyvWGvU6d4MQTXRAUKSMKdyIiIkUxbpzvqVgmTnT33B086LpwFy+GJUvc4++/w8KFbvOq\nWdMN9vCEvaijR/1/LhLWFO5ERESKorCpWKKj8068nJrqllPLHvi2boWvvnIb0BWgYcOc3bkdO+ac\nkkWjdKUYFO5ERESKqihTsWQXGwvnnOM2r82bc4S9o4sWUXHLFvj0U7d5tWjhwl5EhBu5e/iw269R\nulIIhTsRERF/atQILrnEbcC3c+bQq1GjrJa9JUvclCurV7vNl7Q0GDFC4U580h2dIuJ3SUmQkOAa\nJBISslaIEimXIiKgZUsYNAheeMHdn7d3r5uK5eWX8//cpk3uc0OGwOTJ7v4+a/1XbwlaCnci5USw\nBKqkJNejtH69+3fI28MUyIAXLNdGJJN30uVbbnGjdH0xxrXsvfEG3HCD68atXx8uuwyefdbd66fB\nGuWSumVFygFvoPIO8ivqLTvWui0jw23HjmU9L+nr++7LOdgQ3Ovhw+Hss92/aVFRbgUpfzRClPTa\niPhNfqN0X34ZTj4Zvv02a9u+HaZNcxu4SZW7dIGzznJbly5u6hYJawp3In5S2oPdDh2Cv/8ueNuz\nxz1++qkrn11aGlx7Ldx6a/5hLCPDNQ5ERECFCu4x+5Z7X2GvIyLcQEFfNm92CwUcPuyWAT1yBI4d\n65kZ9ryBL/tjaex74AHfYfO++6B7d/dvaEyMm67MmJL/eRWHBkZKDoWN0u3UCe66y/1v6I8/coa9\n1avh66/dBu4v5WmnuaDXvTt06+Za+0A/eGFE4U7ED3y1Dt10kwteffrkH8oK2qx1A/Fq1sx/S0hw\nj++/77te1sK6dfkHN2NKP9AkJPieBzY+3tUlu6+/nke3br1yBD7v85Lu+/vvnPt27fJdz61b3b97\nBw64P7f09Kyg533MvvnaV5yy3kUM1JIoPhVllK4xbjLlk07KWl1j2zZYsCAr7C1blrW02vPPuzIn\nnggNGsD337u/FKAfvBCncCdSRtLS3P3Nq1a51rHcrUMHD8Idd7jfw4WFM19b5cpFr8uIEb4DVdOm\nOafS8of8epjGjctbtkIFd57FOdfimj27aGHz6FFX5wMHsrbcr3Pv27mzaOUOHHD/LsfEwP79rgU1\nu7Q0d+vVjz9C9epQrZp7zG+rWtVdu9KiBp0QFheXY2QuBw64EbnffOPC3sKFrrXvjz/yflYjckOW\nwp3IcbAWtmxxAe6333I+btsGzZu7wWx79uT/+d9+K/t6FidQlbXCepj8rajXpmLFrPBU2qx1DSZp\naVC7tu8y+/e73rO9e10Y3bvX97ZvnysbHV1wAMwvJP7+e1WaNMl6/dFHakkMKzEx0Lu328D9r2XF\nCujQwXf5TZvcwI7+/WHAADdBc1SU/+orJaJwJ+VScVsiDh6ENWti2L49Z4BbvdoFgZYtoVUr99i/\nv3tMSMhqPcmvK7Jp07I4u7yCLVAVdx7YshQM18YYd/9fpUru+/NrSbz33qIdLyPDNdAUFAD37nVd\n0mvX5nxvy5ZWPP541uv09LzHT0uDm292PXyxsVm3B3ifZ39dqdLxXZvc1IpYyipWdOEtv3VzjXFN\nxj/+CE884f5H0KePC3r9+0OzZv6vsxQqZMKdMWYA8DxQAXjdWvtErvcrAW8BHYBdwBXW2nX+rqcE\nv4LuaerdO28L3G+/uVa4uLhT6NAhK8DdcYd7XpRuzWBoOQumQBVsgunalMbPSkSE+ze4WjU3X25x\nJCcvpVevXjmO5WvU8oEDriUxNdW1Xqemuu3vv7Oep6a6ewkLCn8FvRcTk/OeT92PWIby+8F76SX3\nQ/Tll/DFF/DLL/DZZ24DN/3KgAFu69nTfUYCLiTCnTGmAvAScA6wCVhijPncWvtrtmJDgFRr7YnG\nmCuBJ4Er/F9bCXYjR/oeHXnNNVC3blYLXKtW0K9fVivct98uyfGPXnEEQ+uQhIZg+1kpqCXxvvsK\n/qy17u9WfsHv779dy2F+76WnZ91jGhvreg99jfq++243I0jdulCnjuuSlmIq7AfvnHPg3/+GjRth\n1iwX9GbPzlpF44UXXDNtjx5ZXbinnOLSuZpb/S4kwh1wBvCHtfZPAGPMVOAiIHu4uwgY43n+EfCi\nMcZYq+m6xUlNhf/9z/c/VF7btpXd9wdT65AEt2D6WTmelkTvIJGYGGjcuPjf7R3d7A1+Xbv6Lrd9\nu1ukYccOt0VGupDnDRmIEukAABZzSURBVHuFPcbGuhbKogjrnFKUH7wmTdzFHjLE3a+3eLELel98\n4SZNnj3bbffe6/7Qmzd3gzY0CtevTChkH2PMP4AB1tobPa+vATpba2/LVmalp8wmz+s1njI7cx1r\nKDAUIC4ursPUqVP9dBala//+/VTVRJSF2rUriu++q8M339QhJaU67dun8tNPNdi7N+8NwXFxh5g6\ndVG+x9I19z9dc//zdc2/+qoer79+Atu3V6JevcPceOOf9O273e91u/LKLmzblnfodPa/u9bCwYMV\n2LMnkr//jsx83Ls3kr//jsqx37ulpVWgWrWj1KiRTo0a6dSsmZ75vEaNI5n7fvutGu+8E8+RI1lD\nkStVOsa99646rusRLj/nkXv2ELt0KbUWL6bWkiVEpabmW/ZQXByLAvjvb6he8969e/9gre1YWLlQ\nCXeXA/1zhbszrLW3Zyvzi6dM9nB3hrU2n1msoGPHjnbp0qVlW/kykpycXOIuwnD355/wySdugvaU\nFBg4EC691PUUxMTkvW8HXEvEq68W/B9JXXP/0zX3v2C+5iX9u1uY9HTYvdu1+u3cmf/jvHm+B5hU\nqQLDhkHDhjm3Bg3c75zCBPM1L7GMDNePfvrp+ZeZONH9cq5b13/18gjVa26MKVK4C5Vu2U1Ak2yv\nGwNb8imzyRhTEagB7PZP9SSQrIVff81acWfLFrj4YnjoITeoK/eo/WC7p0lEiqas/u5GRrrp4OLi\nCi6XX9dtWhrUq+fq9P337neQd4uKyhv6soe/hg3h8OHC+4RDrjs4IsKthJHfKFxwkzf+61/uF/Xl\nl7ugV6eOf+sZpkIl3C0BTjLGNAM2A1cCV+cq8zlwHbAQ+AcwR/fbhS9rYcmSrEB3+LD7vfDCC+6+\nnMImcA2me5pEpOgC+Xe3oMEl99+fd7+17n5Bb9DbutU9rlnj5hDOCoFnUbVq3tDn3X7+GZ56yk3J\nBCF225qvmzajo90Itk2b3OCMr75ymzfo/fOfbtLl/CZ9lEKFRLiz1h41xtwGfImbCuUNa+0vxphH\ngaXW2s+BScDbxpg/cC12VwauxlJSBf3v9OhRN6H6tGmu27VqVRfo3nvPTdPkr3U/RaR8Ku7gEmOy\npndp3Tr/486dO5+2bXvlaPHzTo4+dy7MmJE1HsErLc2tfFOhghuz0Lw51Kp1/OdY6gprbk1NdYtf\nf/CBC3jeARm33AJ9+7oWvWPHYPz4EGq2DLyQCHcA1tqZwMxc+x7O9vwQcLm/6yWlJ785rH78MWuk\na9OmLtDNmuWmPhAR8Zey6hY2xjVS1a4NbdrkfT+/7uA9e+Djj11L4Jo1rtyJJ2aFPe/mXTq2sBHB\nZdb1W1Bza2ysWwd38GA3q/ann8KHH7qg9+WXbssupJotAydkwp2Ev/zmn3vxRTcx+kMPufnmREQC\nJRDdwgV1B3/4oXturctGa9a4ZWLXrHEDQN54wz3fs8ctJpE98Hmfe48T8Amia9fOmmZl504X9G6/\n3ffkhtdf7yZSbtMma2vWzCXYkLtBsfQp3EnQ2LDB9/4jR+DOO/1bFxGRYFGU7mBj3FiEOnWgc+e8\nx9i/380k4G3lW7nSZaM//oDNm13P57FjOT+TluYyUkByUZ06cOONWQkzt6NHXSL1pltwQ5Pr14d1\n67JOppy29CncSVA4cMD9vdy/P+97/lp/VUQkGJVGd3DVqtC2rdtyS0/Pfw3g9evd2IdOndx22ml+\nXgEkv2bL+vVdl87PP2dtW7e65JpbQFNqYBRxTm6RsrNggfuF0a5d3l8a/l5/VUQkGCUmugapjAz3\nWJo5JTIy//9EN2gAvXq56aZuv931nJ52Gtx0E7zyCixblnewR25JSe6WmogI95iUVIzKjRuXd73a\nKlXcUmjXXecev/zSjUDZuTP/kXX5dQ2FKYU7CZjDh2HECLjsMnjySTcS9rXX3P0fxrjH452cVERE\n/n979x+kV1Xfcfz9zQKJSRAYUEoJSXBErTAqEiNC0ASLUhzDVJlOmAWaRg1aEAn+xDgGcLZWGccf\nTNCEEgEnEMDpYFpAIC2ZkkioICqEiI2SX42RlCASItGQ0z/OrmyS3eyzP55797n3/ZrZyT73uYHv\nnll2P5xzz/n2rbcMdfXV+RG473wndxfbti3/nD7xRFi1Kuerww7LS8EXXww33ACrV7+8Ktq1UW79\n+vxcYNcqacMBr709/yJo5BfD4Yf3nlJrtgTksqxK8eijcMEF+aHen/0sHwAKnj8nSWVodOl31KiX\nl2i7bN+ef6Y//HA+yaCjA7ZsyQHwpz/teaNcv1ZJ+/OLYTDNkCvEcKdC7doFX/4yXHMNfO1rcN55\nnk8nScPBQP/neuxYOO20/NHl2WfhkUfgjDN6/jtNWyXt+gI+9jF4/nk45BCYP792swYuy6owa9bk\n7hErVuTnNM4/32AnSVV02GH5DOIJE3p+v60NvvQl+OUvm/Avb2+Hm27Kn594Yu2CHRjuVIDdu+Hr\nX4d3vhNmzYIf/hDGjSu7KklSs/X2LN/cubB1K7zrXbnD0Fe/mjeKDJlTT81/PvRQ3zs+Kshwp6Z6\n6imYNi23DFu1KneUcbZOkuqht/0QV1yRe4Fv2pQf0fnVr/JzfO94B3zzm3nza5cB7bZ91avgDW/I\nDXl/8pPmfHHDmOFOTZFS3lE1eTJMnw7Ll+eT0CVJ9bK/Y1za2vIEwIIFOdDNm5c3Z5xwQj6CZdas\nfOzKgHbbdj0E+MADQ/9FDXOGOw25zZvhfe/LW+eXL4dPfjL/ByxJUm8OPBDOPDMfp7J5M8yZA7fe\nmiffuuvabdsnw500eCnBzTfn51ff/va8DHv88WVXJUlqNaNGwdln7xvsujS023bKlPznypV52rBG\nDHcasO7PQRxzTA50HR1w1115av3AA8uuUJLUyno7e3jECFi0KLdO69XEiXD00fnk5TVrmlHesGW4\n04Dsfer4pk35OYlPfxpOOqns6iRJVdDbbtvPfz6vFL3+9fsJeRG1XZo13GlA5s7d99TxXbvyDihJ\nkoZCb7ttr7oKli2DG2/cN+Ttsbv2ngUs5tx8wGqN2KFCA9Lb8w41680sSWqy/XXOOO20HPIeeACu\nvBI++1n4/e9fPtpu/bOvZDbXwd2XU6ejjJ2504AcdVTP12vWm1mSNAx0hby2tn3PLN7BGOZuu6xW\nsw+GO/VbSnDooftumKhhb2ZJ0jDy9NM9X9/A+Fo9d2e4U79ddx2MGQPXX7/vcxA1bOEnSRomels9\nGs+GWj135zN36pcNG/JmivvvzyeIn39+2RVJkpR1dOSTHLpv+HvFyJfo2Pl5eODn5RVWMGfu1LCU\n4MIL4ROfyMFOkqThZO/dtWPGwJQpQfvIf4XVq+GZZ8ousRCGOzXsxhthy5a8G0mSpOGoey/bjRvh\n8SdGsPL1s/KbK1eWWltRDHdqyObN8JnPwHe/a+cJSVJrOOww+Na34CP/ewU7Oag2z90Z7tSnlOCj\nH80fb3lL2dVIktS4D34QXtG2kyP5LSOu/mcmHrCJ//5GL9tqK8Jwpz7dfDM89RR84QtlVyJJUv/c\nfNEK1jx9OM9xKIkRrH9pHFf+4H0s/sfqzuIZ7rRfW7bAZZfl5diDDiq7GkmS+mfuwon8gT0b1O5g\nDHMXTiynoAIY7tSrlOCii2DWLJg0qexqJEnqvw0v/WW/rleB4U69uv12WLMG5s0ruxJJkgZmfNvm\nfl2vAsOderR1az7PbtEiGDWq7GokSRqYjtnrGM0Le1wbzQt0zF5XTkEFMNypR5dcks8KOvnksiuR\nJGng2q+dwsKPPcqEtk0Eu5nQtol5Z99J+7VTyi6taQx32scdd8Ajj8BVV5VdiSRJg9d+7RTW7RrH\n7jSCdbvGMfnSV5ddUlPZW1Z72LYtb6JYsgRGj+77fkmSNLw4c6c9XHppPvDxtNPKrkSSJA2EM3f6\nszvvzJ1ZHnus7EokSdJAGe4EwO9+BxdeCN/7HowZU3Y1kiRpoFyWFQCf+hS8//0wbVrZlUiSpMFw\n5k7cey8sW+ZyrCRJVeDMXU0tXgwTJ8KIEXDWWTBjBhx8cNlVSZKkwTLc1dDixTB7Nqxfn/vHvvQS\nXHNNvi5Jklqb4a6G5s6FHTv2vLZjR74uSZJam+GuhjZs6N91SZLUOgx3NTR+fP+uS5Kk1mG4q6GO\nDhg1as9ro0fn65IkqbUZ7mqovR0+8AEYOxYiYMIEWLgwX5ckSa2tZcJdRJwZEU9GxNqI+FwP78+M\niK0R8dPOjw+XUWer2LYNbroJdu+GdesMdpIkVUVLHGIcEW3AfOAMYBPw44hYmlJ6Yq9bb00pXVx4\ngS3mxRdh5Uq45ZayK5EkSUOtVWbuJgNrU0q/Tin9EVgCnF1yTS1r5Uo44QQ49NCyK5EkSUMtUkpl\n19CniDgHODOl9OHO1+cDb+8+SxcRM4EvA1uBXwJzUkobe/hnzQZmAxx55JEnLVmypPlfQBNs376d\nsWPHDujvLljwGkaO3M3MmeuGtqiKG8yYa2Ac8+I55sVzzIvXqmM+bdq0R1JKk/q6ryWWZYHo4dre\nqfTfgFtSSjsj4qPAjcDp+/yllBYCCwEmTZqUpk6dOsSlFmP58uUMtPY5c2D+fDjllIlDWlPVDWbM\nNTCOefEc8+I55sWr+pi3yrLsJuCYbq/HAZu735BSeialtLPz5XXASQXV1lJ++9u8gWLy5LIrkSRJ\nzdAq4e7HwHERcWxEHATMAJZ2vyEijur2cjqwpsD6WsayZTBtGhzQKnO2kiSpX1riV3xKaVdEXAzc\nA7QBi1JKqyPiKuDhlNJS4JKImA7sArYBM0sreBi79154z3vKrkKSJDVLS4Q7gJTSXcBde137YrfP\nLwcuL7quVpJSDnfz5pVdiSRJapZWWZbVEHj8cRgzBl7zmrIrkSRJzWK4q5H77oMzzii7CkmS1EyG\nuxrxeTtJkqrPcFcTXS3Hpk0ruxJJktRMhruaWLEC3vQmW45JklR1hruacElWkqR6MNzVhOFOkqR6\nMNzVwJYtsH49vO1tZVciSZKazXBXA8uWwemn23JMkqQ6MNzVgEuykiTVh+Gu4lLy8GJJkurEcFdx\nthyTJKleDHcV55KsJEn1YrirOMOdJEn1YrirsD/8AX70I1uOSZJUJ4a7CluxAt78ZjjkkLIrkSRJ\nRTHcVZhLspIk1Y/hrsIMd5Ik1Y/hrqJ+8xvYuBEmTSq7EkmSVCTDXUUtW5Y3UthyTJKkejHcVdR9\n97kkK0lSHRnuKigln7eTJKmuDHcV9NhjcPDBcOyxZVciSZKKZrirIGftJEmqL8NdBRnuJEmqL/dS\nVsjixXD55fkIlF/8ArZvh/b2squSJElFMtxVxOLFMHs27NiRX2/cmF+DAU+SpDpxWbYi5s59Odh1\n2bEjX5ckSfVhuKuIDRv6d12SJFWT4a4iXv3qnq+PH19sHZIkqVyGuwp48EF44QUYOXLP66NHQ0dH\nOTVJkqRyGO5a0OLFMGPGyYwYAUcdBe99L9x+O1x/PUyYABH5z4UL3UwhSVLduFu2xby8K3YUAFu2\n5Bm7Z57JQc4wJ0lSvTlz12J62hW7c6e7YiVJUma4azHuipUkSftjuGsxve1+dVesJEkCw13L6ejI\nu2C7c1esJEnqYrhrMe3teRfskUe+6K5YSZK0D3fLtqD2djj66FVMnTq17FIkSdIw48ydJElShRju\nJEmSKsRwJ0mSVCGGO0mSpAox3EmSJFWI4U6SJKlCWiLcRcSiiHg6Ih7v5f2IiG9FxNqI+HlEvLXo\nGiVJkoaDlgh3wA3Amft5/2+A4zo/ZgPfLqAmSZKkYaclwl1K6b+Abfu55WzgppStAg6NiKOKqU6S\nJGn4aIlw14CjgY3dXm/qvCZJklQrVWk/Fj1cSz3eGDGbvHQLsD0inmxaVc11BPB/ZRdRM4558Rzz\n4jnmxXPMi9eqYz6hkZuqEu42Acd0ez0O2NzTjSmlhcDCIopqpoh4OKU0qew66sQxL55jXjzHvHiO\nefGqPuZVWZZdClzQuWv2ZOC5lNJvyi5KkiSpaC0xcxcRtwBTgSMiYhMwDzgQIKX0HeAu4CxgLbAD\n+IdyKpUkSSpXS4S7lNK5fbyfgIsKKme4aPml5RbkmBfPMS+eY148x7x4lR7zyLlIkiRJVVCVZ+4k\nSZKE4W5Yi4gzI+LJzrZqn+vh/ZERcWvn+w9FxMTiq6yWBsb8soh4orPN3X9EREPb0tW7vsa8233n\nRESKiMrucCtKI2MeEX/X+b2+OiJuLrrGqmngZ8v4iLg/Ih7t/PlyVhl1VkmdW5ca7oapiGgD5pNb\nq70RODci3rjXbR8Cnk0pvRb4OvCVYquslgbH/FFgUkrpTcD3ga8WW2W1NDjmRMTBwCXAQ8VWWD2N\njHlEHAdcDpyaUjoeuLTwQiukwe/zLwC3pZROBGYA1xZbZSXdQE1blxruhq/JwNqU0q9TSn8ElpDb\nrHV3NnBj5+ffB94dET0d6KzG9DnmKaX7U0o7Ol+uIp+pqIFr5Psc4EvkIP1ikcVVVCNj/hFgfkrp\nWYCU0tMF11g1jYx5Al7Z+fkh9HJWqxpX59alhrvhq5GWan++J6W0C3gOOLyQ6qqpv23sPgTc3dSK\nqq/PMY+IE4FjUkr/XmRhFdbI9/nrgNdFxMqIWBUR+5v9UN8aGfMrgPM6j/u6C/h4MaXVWmVbl7bE\nUSg11UhLtYbbrqkh/Wljdx4wCXhXUyuqvv2OeUSMID9yMLOogmqgke/zA8hLVVPJs9MPRMQJKaXf\nNbm2qmpkzM8FbkgpfS0i3gF8r3PMdze/vNqq7O9QZ+6Gr0Zaqv35nog4gDyVv78paO1fQ23sIuKv\ngbnA9JTSzoJqq6q+xvxg4ARgeUSsA04GlrqpYlAa/dnyg5TSn1JKTwFPksOeBqaRMf8QcBtASulB\nYBS5/6map+HWpa3GcDd8/Rg4LiKOjYiDyA/YLt3rnqXA33d+fg7wn8mDCwejzzHvXCJcQA52Poc0\nePsd85TScymlI1JKE1NKE8nPOU5PKT1cTrmV0MjPljuAaQARcQR5mfbXhVZZLY2M+Qbg3QAR8Vfk\ncLe10Crrp7KtS12WHaZSSrsi4mLgHqANWJRSWh0RVwEPp5SWAteTp+7XkmfsZpRXcetrcMyvBsYC\nt3fuXdmQUppeWtEtrsEx1xBqcMzvAd4TEU8ALwGfTik9U17Vra3BMf8kcF1EzCEvDc70f9YHp86t\nS+1QIUmSVCEuy0qSJFWI4U6SJKlCDHeSJEkVYriTJEmqEMOdJElShRjuJEmSKsRwJ0mSVCGGO0mS\npAox3EnSIEXEayPiTxFx5V7Xvx0Rz9sLV1KRDHeSNEgppbXAvwBzOnuxEhFfBGYBf2svXElFsv2Y\nJA2BiPgL4FfAtcAvgIXAuSml20otTFLtHFB2AZJUBSmlLRHxDXID+AOASwx2ksrgsqwkDZ3/AUYC\nD6aU5pddjKR6MtxJ0hCIiNOBBcCDwKkR8eaSS5JUU4Y7SRqkiHgrcAd5U8VUYAPwT2XWJKm+DHeS\nNAgR8VrgbuBe4OMppT8CVwJnRcQ7Sy1OUi25W1aSBqhzh+yPyDN1700p7ey83gY8DjybUjqlxBIl\n1ZDhTpIkqUJclpUkSaoQw50kSVKFGO4kSZIqxHAnSZJUIYY7SZKkCjHcSZIkVYjhTpIkqUIMd5Ik\nSRViuJMkSaqQ/weqgfYYJ57PpQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1160b3ba8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot surface pressure coefficient\n",
    "pyplot.figure(figsize=(10, 6))\n",
    "pyplot.grid()\n",
    "pyplot.xlabel('$x$', fontsize=16)\n",
    "pyplot.ylabel('$C_p$', fontsize=16)\n",
    "pyplot.plot([panel.xc for panel in panels if panel.loc == 'upper'],\n",
    "            [panel.cp for panel in panels if panel.loc == 'upper'],\n",
    "            label='upper surface',\n",
    "            color='r', linestyle='-', linewidth=2, marker='o', markersize=6)\n",
    "pyplot.plot([panel.xc for panel in panels if panel.loc == 'lower'],\n",
    "            [panel.cp for panel in panels if panel.loc == 'lower'],\n",
    "            label= 'lower surface',\n",
    "            color='b', linestyle='-', linewidth=1, marker='o', markersize=6)\n",
    "pyplot.legend(loc='best', prop={'size':16})\n",
    "pyplot.xlim(-0.1, 1.1)\n",
    "pyplot.ylim(1.0, -2.0)\n",
    "pyplot.title('Number of panels: {}'.format(panels.size), fontsize=16);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Accuracy check\n",
    "\n",
    "For a closed body, the sum of all the source strengths must be zero. If not, it means the body would be adding or absorbing mass from the flow! Therefore, we should have\n",
    "\n",
    "$$\\sum_{i=1}^{N} \\sigma_i l_i = 0$$\n",
    "\n",
    "where $l_i$ is the length of the $i^{\\text{th}}$ panel.\n",
    "\n",
    "With this, we can get a measure of the accuracy of the source panel method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sum of singularity strengths: 0.004606\n"
     ]
    }
   ],
   "source": [
    "# calculate the accuracy\n",
    "accuracy = sum([panel.sigma * panel.length for panel in panels])\n",
    "print('sum of singularity strengths: {:0.6f}'.format(accuracy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lift coefficient\n",
    "\n",
    "The lift  is given by the Kutta-Joukowski theorem, $L = \\rho \\Gamma U_\\infty$, \n",
    "where $\\rho$ is the fluid density. The total circulation $\\Gamma$ is given by:\n",
    "\n",
    "$$\\Gamma = \\sum_{i=1}^N \\gamma l_i$$\n",
    "\n",
    "Finally, the lift coefficient is given by:\n",
    "\n",
    "$$C_l = \\frac{\\sum_{i=1}^N \\gamma l_i}{\\frac{1}{2}U_\\infty c}$$\n",
    "\n",
    "with $c$ the chord-length of the airfoil."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lift coefficient: CL = 0.506\n"
     ]
    }
   ],
   "source": [
    "# compute the chord and lift coefficient\n",
    "c = abs(max(panel.xa for panel in panels) -\n",
    "        min(panel.xa for panel in panels))\n",
    "cl = (gamma * sum(panel.length for panel in panels) /\n",
    "      (0.5 * freestream.u_inf * c))\n",
    "print('lift coefficient: CL = {:0.3f}'.format(cl))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Challenge task\n",
    "\n",
    "Based on what has been done in the previous notebook, compute and plot the streamlines and the pressure coefficient on a Cartesian grid."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Please ignore the cell below. It just loads our style for the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: rgb(245,245,245);\n",
       "}\n",
       "\n",
       "div.cell { /* set cell width */\n",
       "    width: 750px;\n",
       "}\n",
       "\n",
       "div #notebook { /* centre the content */\n",
       "    background: #fff; /* white background for content */\n",
       "    width: 1000px;\n",
       "    margin: auto;\n",
       "    padding-left: 0em;\n",
       "}\n",
       "\n",
       "#notebook li { /* More space between bullet points */\n",
       "margin-top:0.8em;\n",
       "}\n",
       "\n",
       "/* draw border around running cells */\n",
       "div.cell.border-box-sizing.code_cell.running { \n",
       "    border: 1px solid #111;\n",
       "}\n",
       "\n",
       "/* Put a solid color box around each cell and its output, visually linking them*/\n",
       "div.cell.code_cell {\n",
       "    background-color: rgb(256,256,256); \n",
       "    border-radius: 0px; \n",
       "    padding: 0.5em;\n",
       "    margin-left:1em;\n",
       "    margin-top: 1em;\n",
       "}\n",
       "\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Alegreya Sans' sans-serif;\n",
       "    line-height: 140%;\n",
       "    font-size: 125%;\n",
       "    font-weight: 400;\n",
       "    width:600px;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Alegreya Sans', sans-serif;\n",
       "    font-style:regular;\n",
       "    font-weight: 200;    \n",
       "    font-size: 50pt;\n",
       "    line-height: 100%;\n",
       "    color:#CD2305;\n",
       "    margin-bottom: 0.5em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Fenix', serif;\n",
       "    font-size: 22pt;\n",
       "    line-height: 100%;\n",
       "    margin-bottom: 0.1em;\n",
       "    margin-top: 0.3em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Fenix', serif;\n",
       "    margin-top:12px;\n",
       "\tfont-size: 16pt;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: regular;\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {    /*Use this for captions*/\n",
       "    font-family: 'Fenix', serif;\n",
       "    font-size: 2pt;\n",
       "    text-align: center;\n",
       "    margin-top: 0em;\n",
       "    margin-bottom: 2em;\n",
       "    font-style: regular;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {  /*Use this for small titles*/\n",
       "    font-family: 'Alegreya Sans', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 16pt;\n",
       "    color: #CD2305;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .5em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 { /*use this for copyright note*/\n",
       "    font-family: 'Source Code Pro', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 9pt;\n",
       "    line-height: 100%;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "\t\n",
       "    \n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"], \n",
       "                           equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling(filepath):\n",
    "    styles = open(filepath, 'r').read()\n",
    "    return HTML(styles)\n",
    "css_styling('../styles/custom.css')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}