{ "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)2014 Lorena A. Barba, Olivier Mesnard. Thanks: NSF for support via CAREER award #1149784." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Flow over a cylinder with source-panels" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In previous lessons, we used potential-flow singularities of defined strength to represent the shape of simple geometries, such as a [Rankine oval](02_Lesson02_sourceSinkFreestream.ipynb) or a [circular cylinder](03_Lesson03_doublet.ipynb), immersed in a free stream. We were rather lucky that when superposing a few fundamental potential-flow solutions, the stream-line pattern that resulted had a closed dividing stream line that we could interpret as a solid body.\n", "\n", "But what if we want to represent the stream lines around an *arbitrary* geometry? Would you be able to define the combination of fundamental solutions to get the expected result? *How could you do that?* Trial and error? It would take enormous luck and a lot of work to get a geometry we want.\n", "\n", "In this lesson, the objective is to calculate the source-strength distribution that can produce potential flow around a given geometry: a circular cylinder. We know that we can get the flow around a cylinder by superposing a doublet in a free stream, but here we want to develop a more general approach that can later be extended to *different* shapes.\n", "\n", "The method we will use consists of representing the geometry of the body by a series of short linear segments, called *panels*, that correspond to [source sheets](08_Lesson08_sourceSheet.ipynb) like we studied in the previous lesson.\n", "\n", "What we are aiming for is a method that starts with a defined body geometry, then works out the strength of sources needed in each panel representing that geometry to get a dividing streamline right on the body boundary. We will have to *solve* for the source strengths by specifying that the body be a streamline, i.e., the velocity be tangent there.\n", "\n", "Let's start by loading the Python libraries that we will need in this notebook." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import math\n", "import numpy\n", "from scipy import integrate\n", "from matplotlib import pyplot\n", "# embed the figures into the notebook\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will add a uniform horizontal flow of magnitude `u_inf`, so let's make that equal to 1:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "u_inf = 1.0 # free-stream speed" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Definition of the geometry" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The geometry considered here will be a circular cylinder of unit radius. We can define this geometry very easily by a set of points going around the angular range between $0$ and $2\\pi$." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAASkAAAEPCAYAAAAebDKpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJztnXl4VOX1xz8nISQIRZa4sIdVQbQg\nCKJU2WutIlpEbVWw4orWyk+rlLrTulWlahWttVC0AioKKtaqECjIjiwqZRFQEFB2CYRA4P39ce6Y\nATLJJJmZe2dyPs/zPu+de99773mzfOddzxHnHIZhGEElzW8DDMMwSsJEyjCMQGMiZRhGoDGRMgwj\n0JhIGYYRaEykDMMINCZShmEEGhMpwzACjYmUYRiBporfBvhNdna2y8nJKbXcnj17qF69evwNSgCp\nUpdUqQdUzrosXLhwq3PuuNLKVXqRysnJYcGCBaWWy83NpVu3bvE3KAGkSl1SpR5QOesiIl9F8zzr\n7hmGEWhMpAzDCDQmUoZhBBoTKcMwAo2JlGEYgSZwIiUiL4vIdyLyWYTrIiJPi8hqEVkqIqeHXRso\nIqu8NDBxVhuGES8CJ1LAaOC8Eq7/DGjppeuB5wFEpA5wH9AZ6ATcJyK142qpYRhxJ3DrpJxzM0Qk\np4QiFwH/dOr3eI6I1BKRekA34EPn3HYAEfkQFbvX4muxESvy8+Gbb2DDBs03b4bduyEv7+jkHHz/\n/WkcfzxUqVKUsrKgbl3Izj46NW4MNWr4XUujrAROpKKgAbA+7PMG71yk80chItejrTBOOOEEcnNz\nS31pXl5eVOWSAT/rUlgorF9fjbVrq7NuXXXWrq3Oxo3V2Lo1k++/zyjj0+qU+f116hTQoEE+DRrk\nU7/+Pho0yKdx4z3k5OylShX//P3b31dkklGkpJhzroTzR5907kXgRYCOHTu6aFbHVsYVwRXl4EH4\n7DP473/hk09g2TJYsQIOHCi+fEYGNGgADRtqqlcPataE6tW1BRRK1atDWhosWrSUNm1Oo7CQH1J+\nPmzbBlu3Hp62bIF162D79ky2b89k2bJah707KwvatYOOHYvSySdDenrcf0yA/X2VRDKK1AagUdjn\nhsBG73y3I87nJswqg/37Ye5cFaWZM1WYdu06ulyzZtC2bVFq2RIaNYLjjlPxiZYqVbZTlv+Fgwe1\nG7l69eFp2TLN58zRFOJHP4Jzz4XevTWdfDJIcV+FRlxJRpGaDNwiIuPQQfJdzrlNIvIB8KewwfI+\nwDC/jKwsbN8O778PkyfDv/8N339/+PUmTeAnP4Gzz4YOHaB1a//GhdLTdVyqcWPo0ePwazt2wMKF\nsGCBpvnz4euv4d13NYG28kKCdf75UKvW0e8wYk/gREpEXkNbRNkisgGdscsAcM6NAqYA5wOrgb3A\nNd617SLyEDDfe9SDoUF0I7Z89RW88Qa88462mA4eLLrWujV07w5du2pq1Cjyc4JE7drQq5emEBs2\nwIcfavroI22FjR6tKSMDfvpTGDAALrpIu6VGfAicSDnnrijlugOGRLj2MvByPOyq7OTlwZtvwpgx\nMG1a0fkqVbRV0rcvXHihduVShYYN4ZprNB06BEuXqmBNmQIzZhS1sjIz4bzz4LLLVLCOOcZvy1OL\nwImUERwOHYLp01WY3ngD9uzR81lZKkoXX6z/nJWh25OWpgPr7drBnXfq8oiJE2H8eB2DmzRJ07HH\nwsCBcOON2qo0Kk4QF3MaPrNnDzz/vP6T9eihIrVnj44rvfii/oOOHw+XX145BKo4TjwRbr5ZRXzD\nBnj6aejcWScKnn4a2rSBbt1g3DidUDDKj4mU8QMbN8Lw4TqwfPPNsHKljin94Q96PHMmXHedthaM\nIurXh1tv1ZnBRYvghht0mcT06XDFFfozfOgh2LnTb0uTExMpgy++gKuvhpwc+NOfdMauc2eYMAHW\nrNF/sJYt/bYyOWjfHkaNUsF/7jk49VT47ju4916d6fz973XNlhE9JlKVmDVrVJzatoWxY3WWrn9/\nXd80Zw5ceqkOjBtlp2ZNuOkmWLIEpk6Fnj11ecbDD+uXwdChKmRG6ZhIVUK2bKnKTTfBSSepOKWn\n60DvqlXw+uvQpYvfFqYOIrok46OPYPZsuOAC2LsXnnpKZ0Lvuqv4Ba9GESZSlYhdu3Rm6sorOzNq\nlM7eDRyo403PP59ayweCyJln6tqyTz/VFmtBATz2GLRoAZMm1aew0G8Lg4mJVCXAOXjlFW05/fnP\nsH9/Ov3767660aOhaVO/LaxctGunLdZ583TB69atMHJkK047TddgOf/2OQcSE6kU5/PPtbtx1VXw\n7bdw1lkwatQCXn/d1vH4zRln6KLQN96A+vXzWb4cfv5zTV9FFeypcmAilaLs3g133KHf2tOn6+bd\nf/xDFx6edFKe3+YZHiLwi1/AP/4xjyee0OUd778Pp5wCzz6rXfLKjolUCjJjBpx2GjzxhM7Y3XST\nukgZNKhsXgaMxFG1qmPoUPjf/3S8as8eXXt1zjl6rjJjf7IpREGBzhZ166a+k9q313GP557TDbRG\n8DnxRB2vmjhRj2fNgh//WNevVdaBdROpFGHpUujUSWeLRHSV+Jw56rzNSD4uvlgX2V57rW6rGT5c\n11p9843fliUeE6kkxznt1p1xhgpVixa6feWhh6BqVb+tMypC7drw0kvwn/+ol9IZM3SMccoUvy1L\nLCZSSczu3Tp+cccd+m17ww26BscWY6YWvXvD4sXqv2rrVp39u/POyG6YU43AiZSInCciK7y4encX\nc/0pEVnspZUisjPs2sGwa5MTa3liWbFC99dNnKhbMN5+W/eMWTSU1OT447UF9cgjukPgz39Wj6cb\nNvhtWfwJlEiJSDrwVzS2XhvgChFpE17GOXe7c66dc64d8AwwMexyfuiac65vwgxPMJMm6fjT8uXq\nEmT+fHW2ZqQ2aWk6MTJjhnqqmDtX/w7mzy/93mQmUCKFBvVc7Zxb45zbD4xD4+xF4goqUVw95+C+\n+6BfP92s2r+//qG2auW3ZUYiOessdQnTvTts2qTLFN54w2+r4kfQRKossfOaAE2BqWGns0RkgYjM\nEZF+8TMz8Rw4oG5sH3xQv1EffVRdqVj3rnJSt64Gvrj2Wti3Tz1W/PGPqbmlRlyAaiUilwI/dc4N\n9j5fBXRyzt1aTNm7gIbh10SkvnNuo4g0Q8Wrp3Puy2LuDQ8O2mHcuHGl2paXl0cNnxQhPz+NBx44\nhblz65KVdZD77vucM88sf4wJP+sSS1KlHlD+ujgHr7/ekFGjmuOc0Lv3Zu68cwUZGf4GOo2mLt27\nd1/onCt9kYxzLjAJ6AJ8EPZ5GDAsQtlPgbNKeNZooH9p7+zQoYOLhmnTpkVVLtZs2eJcp07OgXPZ\n2c7NnVvxZ/pVl1iTKvVwruJ1eftt5445Rv9Ofv5z5/LzY2NXeYi2LsACF4UuBK27Nx9oKSJNRaQq\ncDkaZ+8wROQkoDYwO+xcbRHJ9I6zgbOBLxJidZxYt079is+bp47SZs3SgVLDOJKLLtIB9bp14b33\n1G9VKHBGshMokXLOFQK3AB8Ay4EJzrnPReRBEQmfrbsCGOepcYjWwAIRWQJMAx5xziWtSK1dqwOi\nK1fqtohPPrEBcqNkOnSA3Fw44QT4+GNdV5UKDvUC5xzWOTcFDQAafu7eIz7fX8x9nwCnxtW4BPH1\n1xqlZf16ncmZMsWCHxjR0battqh69tSWd69e8MEHUKeO35aVn0C1pAzdm9Wjh3b1OndWtx0mUEZZ\naNVKXfI0a6Yh43v2TO4WlYlUgNi0SQXqyy+16f7vf1v4bqN85ORoi6plS91S07cv5Of7bVX5MJEK\nCNu26TfeypW6ifQ//6m8gTeN2NCggYaFr19fBevyy5PT3YuJVAAoKFDXHMuX65jChx8m9xiCERya\nNNEvvNq1YfJkuP765FvwaSLlM87Br3+tYwgNGmgXLzvbb6uMVOKUU3RZwjHHqAvpu+7y26KyYSLl\nM/fdB//6l25vee89FSrDiDVdusCbb2qw18cfVz9VyYKJlI+MGaPO6dLSYPx4XQ9lGPHivPPgxRf1\n+Oabde1dMmAi5ROzZsF11+nxs8/C+ef7a49RObjmGg3wcOCARqlJBnfEJlI+sGULXHaZ/qH89rca\nzcUwEsUTT2iwjs2b4ZJL1ItCkDGRSjCHDmmgzm++0X15jz3mt0VGZSMjQ938NGmi+0JvuinYM34m\nUgnm4Yd1m0J2Nowbp38whpFojjtOXU5XqwajR8Mrr/htUWRMpBLItGlw770acmrsWGjY0G+LjMpM\nu3Y6HgowZIhuag8iJlIJYutW+OUvtbs3fLjOtBiG31xzjY5L7d6twxBBXJFuIpUgfvMbHag85xy4\n/36/rTEMRUSXJdSrpzPOjzzit0VHYyKVACZNgtdeK1rxm57ut0WGUUTdurpmD/QLdN48X805ChOp\nOLN9O9x4ox4//LC6zzCMoNG7ty6HOXgQBg3SYLNBIXAiFUVw0EEisiUsCOjgsGsDRWSVlwYm1vLi\nuf127eZ17Qq33OK3NYYRmYcfVtcuy5frWqqgECiRiiY4qMd4VxQE9CXv3jrAfUBnNH7ffSJSO0Gm\nF8uUKfDPf0JWFrz8sm5/MYygkpUFzz2nxw89FJzZvqD925Q1OGg4PwU+dM5td87tAD4EfJtDKyjQ\nwXLQX3jLln5ZYhjR06sXXHGFOsi79dZgLPIMmkhFGxz0FyKyVETeEJFGZbw3ITzzjHrYbNNG+/qG\nkSw8+aR6hH3vPV3w6TdBC8QgxZw7UsvfAV5zzhWIyI3AGKBHlPfqSw4PDkpubm6phuXl5UVVDmDH\njgzuv78zUIWBA5cwc+aOqO5LFGWpS5BJlXpA8OoyaFADnn66JTfcsI9q1eaRlXUo6ntjXpdogvMl\nKlGG4KDe9XRgl3d8BfBC2LUXgCtKe2c8goNef70GaTz//KhvSSipElQzVerhXPDqUljoXPv2+nf8\n6KNlu7fSBwcVkXphH/ui8flAY/X18YKE1gb6eOcSypIl6lCsSpVgzZAYRllITy9a2PnII7Bzp3+2\nBEqkXHTBQX8jIp97QUB/Awzy7t0OPIQK3XzgQe9cQvnd73Try5AhcPLJiX67YcSO3r3VpcuOHerN\n0y+CNiZVanBQ59wwtBtY3L0vAy/H1cASmDtXnd7XqAH33OOXFYYRG0R07VSXLjBypM72nXhi4u0I\nVEsq2XnoIc1vuUW3GhhGsnPmmdCvH+zdCyNG+GODiVSMWLSoKCLH0KF+W2MYsWPEiKKNyH64GzaR\nihGhb5kbb1SHYoaRKpxyCvTvr+6un3km8e83kYoBy5bBW2/ptoI77vDbGsOIPf/3f5q/8ALk5SX2\n3SZSMeDJJzUfPFj98hhGqtG5s/rk37lT96EmEhOpCrJ9u/oqB7jtNn9tMYx4EmpNjRypLl0ShYlU\nBRkzRkMC9ekDLVr4bY1hxI++faF5c/WOMHFi4t5rIlUBDh2C55/XY4udZ6Q66enqHw1g1KjEvddE\nqgJMnQqrVmnUlwsu8Nsaw4g/v/oVZGZq5KOvvkrMO02kKkCoFXX99bpXzzBSnVq14OKL1c/U2LGJ\neaeJVDnZsQMmT1Zvm4MHl17eMFKFQYM0HzMmMU7xTKTKyVtvaYyyHj1s2YFRuejVC+rXh9Wr4ZNP\n4v8+E6lyMmGC5gMG+GuHYSSa9HQNJAoaoj3emEiVg23b4KOP9Jd18cV+W2MYiefqqzWfODH+UY9N\npMrBxIm6mK1XL8jO9tsaw0g8bdpocJHt22H27Pi+y0SqHFhXzzDgwgs1f+ed+L4ncCIVRXDQoSLy\nhRct5mMRaRJ27WBY0NDJR94bC3bvhtxcndXr1y8ebzCM5KBSilSUwUE/BTo6504D3gAeC7uW74qC\nhvYlDkyfrn3wTp2gTp14vMEwkoOzz9Z1U//7n870xYtAiRRRBAd1zk1zzu31Ps4BGibSwI8+0rxX\nr0S+1TCCR0YG/OxnehzP1lTQ1kkXF+CzcwnlrwXeD/ucJSILgELgEedcsaENKxJ3b9KkM4DqZGd/\nSm7urlLvCyJBi/FWXlKlHpC8dWna9HigDePHb6N9+2VA6sfduxR4KezzVcAzEcpeibakMsPO1ffy\nZsA6oHlp7yxL3L2NGzUO2THHOFdQENVtgSRoMd7KS6rUw7nkrcvXX+v/RK1azh08qOdSPe7eBqBR\n2OeGwMYjC4lIL2A40Nc5VxA675zb6OVrgFygfSyN+/hjzc89F6pWjeWTDSM5adRI086d8MUX8XlH\n0EQqmuCg7dHoxH2dc9+Fna8tIpnecTZwNhDTH9t//6t5jx6xfKphJDddu2o+c2Z8nh8okXLRBQd9\nHKgBvH7EUoPWwAIvaOg0dEwqpiK1cKHmZ5wRy6caRnJz9tmaz5oVn+cHbeA8muCgxc6rOec+AU6N\nl10HDgjLdFyQ9jHtRBpGchNqScVLpALVkgoy69ZVZ/9+aNUKatb02xrDCA5t22rU7rVrYcuW2D/f\nRCpKVq78EQCnn+6zIYYRMNLTdS8fwPLlsX++iVSUrFxZA4AOHXw2xDACSEik4jHDZyIVJV9+qSJl\n41GGcTQmUgFg06YsQMekDMM4HBMpn9m7F7ZvzyQjQ92mGoZxOCZSPrNuneZNmuggoWEYh9OkCVSr\nBps2QV5ebP9JTKSiYM0azZs189cOwwgqaWnQoIEeb9uWGdtnx/RpKcratZo3beqvHYYRZEJRk7Zt\ni+3GVhOpKAh190ykDCMyJlI+snWr5scf768dhhFkQiK1fbt19xLOjh2a167trx2GEWRCM9/WkvKB\nnTs1r1XLXzsMI8hYd89HrCVlGKVz7LGa79tnSxASTqglZSJlGJHJ0k0Z7N8fW1kJnEhFEXcvU0TG\ne9fnikhO2LVh3vkVIvLTWNkUaklZd88wIhMSqYICn0RKRGaJyFUhF73xIMq4e9cCO5xzLYCngEe9\ne9ug7oZPAc4DnvOeV2H279c89EswDONoqlXT3M+W1AFgDLBRRJ4UkZNjaolSatw97/MY7/gNoKeI\niHd+nHOuwDm3FljtPc8wjATge3fPOdcN9SM+Brga+FxEckXkMhHJiJE9xcXdaxCpjOcTfRdQN8p7\nDcOIE6FhkbVra8T0uWXyce6cWwEMFZFhwAA0wOa/gK0i8g/gRS+cVHmR4l4bZZlo7tUHlDE4qHPn\nAGlMnz6djIxiH5lUJGsgyiNJlXpAatRl8eJjCUWRC0xwUOB0NL7dIS8VAq8DJ5bzeV2AD8I+DwOG\nHVHmA6CLd1wF2IoK1GFlw8uVlKIJDpqRoQEQkzkgaDjJGojySFKlHs6lRl2WLNH/k+OOy4+qPPEK\nDioi1UTk1yIyD42TdxxwG1AfuAk4C3i1zGqplBp3z/s80DvuD0z1KjwZuNyb/WsKtATmldMOwzDK\nSGiCqXbtAzF9btTdPRE5FbgB+BVQHZgE3OWcmxZW7G8ishltTZUZ51yhiITi7qUDLzsv7h6qupOB\nvwNjRWQ1sB0VMrxyE9CAoIXAEOfcwfLYcSRVqsCBA/pLsMjFhlE8+/ZpXrXqoZg+tyxjUkvQkOcj\n0bGnTRHKrQZml9cgV3rcvX3ApRHu/SPwx/K+OxK1akF+vi7qrBHbMUHDSBny8zWPtUiVpbt3KdDE\nOfdACQKFc265c657xU0LDqGV5qGV54ZhHE2oJZWZGZMOzA9E3ZJyzr0Z0zcnEaGV5qEpVsMwjmb3\nbs0zM/1rSVVaQi0pEynDiMzmzZrXqbM/ps81kYqCUEvKunuGEZmNGzWvW9dEKuHUrat5yEOnYRhH\ns8kbqa5TpyCmzzWRioLGjTUP+To3DONoQiJlLSkfCIWyWlORDT+GkeKYSPlISKRCoa0Mwzgc5+Cb\nb/S4bl3r7iWcUCirtWv1l2EYxuFs2qRLEOrUgZo1C2P6bBOpKKhZE2rWPEB+Pnz7rd/WGEbw+OIL\nzdu0ASnOH0kFMJGKkvr1dc3/6tU+G2IYASRcpGKNiVSU5OTsAWDxYp8NMYwAYiIVAFq1ygNg0SKf\nDTGMAGIiFQBatdKNSQsX+myIYQQM5+Dzz/W4devYP99EKkqaN88jLU1/GSGXFIZhwKpVsH07nHAC\nNIhDVAETqSjJyjpE69Zw8CAsXeq3NYYRHGbN0vzss2M/swcBEikRqSMiH4rIKi8/Kl6wiLQTkdki\n8rmILBWRy8KujRaRtSKy2EvtYm1jhw6aL1gQ6ycbRvIyc6bmXbvG5/mBESngbuBj51xL4GPv85Hs\nBa52zoUCgI4UkfC4wnc659p5KebzcGedpXmSB/UwjJgS3pKKB0ESqfCgn2OAfkcWcM6tdM6t8o43\nAt+hgSASQq9emk+dqt0+w6jsbNkCK1Zo9OL27ePzjiCJ1Akht8RefnxJhUWkE1AV+DLs9B+9buBT\n8QgH36wZ5OToIOGnn8b66YaRfIRaUZ07Q0asQgQfQZmCg1YUEfkIOLGYS8PL+Jx6wFhgoHMu5Kt0\nGLAZFa4XgbuAByPcX6bgoKDBG6dPz+WUU1qxbl19XnxxDXl5X5fF7MCQCoEoIXXqAclbl5deagXU\nJydnLbm5XwFxqEs0wfkSkYAVQD3vuB6wIkK5msAi4NISntUNeDea90YTHNS5ouCN48drAMSePaO6\nLZCkQiBK51KnHs4lZ10OHnSuXj39f1i0qOh8tHUhXsFB40h40M+BaFy/w/AChr4F/NM59/oR1+p5\nuaDjWZ/Fw8gePTSfORP27InHGwwjOVi0SL0fNGwI7WI+l15EkETqEaC3iKwCenufEZGOIvKSV2YA\ncA4wqJilBq+KyDJgGZANjIiHkdnZ0KkTFBTA++/H4w2GkRy8847mF14Yn/VRIRI6JlUSzrltQM9i\nzi8ABnvHrwCvRLi/R1wNDGPAAJg3DyZMgP79E/VWwwgW4SIVT4LUkkoaLvXiJ7/7LuTl+WuLYfjB\n+vU6w129OnSPcyhgE6ly0LgxdOmie/jee89vawwj8bzi9WfOPx+ysuL7LhOpcjJggOYTJvhrh2Ek\nGudgjLfsetCg+L/PRKqchLp8U6bArl3+2mIYiWTuXF1lfuKJ0KdP/N9nIlVOGjTQ5Qj79sHYsX5b\nYxiJY/Roza+8EqokYOrNRKoC3Hij5s89Z1FkjMrBvn0wbpweDxxYctlYYSJVAfr10ybv8uUwY4bf\n1hhG/Jk4UYc3OnSAtm0T804TqQqQkQHXXafHzz/vry2GEW+cgyef1OPQ330iMJGqINddB2lp+g1j\nMfmMVGbGDPXxn50NV1+duPeaSFWQRo10xe2BAzBqlN/WGEb8eOIJzW++Wf1HJQoTqRhw++2a/+Uv\n8P33/tpiGPFgxQrdBpOZCUOGJPbdJlIx4Nxz4Sc/gR07dKbPMFKNp57S/Kqr4PgS3VHGHhOpGHHP\nPZo/8YS5cDFSiw0bilaYDx2a+PebSMWIXr3UherWrTY2ZaQWDz6o66MuvTQ+wT9Lw0QqRogUtaYe\nf9wCiBqpwcqV8PLLkJ4ODz3kjw0mUjHk/PPh9NN1KcLIkX5bYxgV5557NDLSNdfASSf5Y0NgRCqa\n4KBeuYNhXjknh51vKiJzvfvHe66GE4oIPPaYHv/pT+pa1TCSlUWL1MtHZibcd59/dgRGpIguOChA\nvisKANo37PyjwFPe/TuAa+NrbvH07AkXXaTO8IaXKQaOYQQH5+Bu7z/wllvUj7lfBEmkSg0OGgkv\n+EIP4I3y3B9rHn9ct8yMHq0rdA0j2Zg4ET78EI49tkis/EJcQLbvi8hO51ytsM87nHNHdflEpBBY\nDBQCjzjn3haRbGCOc66FV6YR8L5zrtgtkEfE3eswLrStuwTy8vKoUaNG1PV5/vnmTJjQiNNO28nI\nkYvj6qi+rJS1LkElVeoBwarL3r3pDBzYia1bM7nttpX067exTPdHW5fu3bsvdM51LLVgNHGvYpWA\nj9BQU0emi4CdR5TdEeEZ9b28GbAOaI6GWl8dVqYRsCwam8oady9aduxwLjtbY5K9+mqZbo07yRjj\nrThSpR7OBasuQ4fq323Hjs4VFpb9/qSOu+ec6+Wca1tMmgR8GxY7rx7wXYRnbPTyNUAu0B7YCtQS\nkZALroZA2eQ/xtSqBY88ose/+Y1tPjaSgyVLdHtXWpqu90tP99uiYI1JRRMctLaIZHrH2cDZwBee\nKk8D+pd0f6L59a91kee2bXDrrX5bYxglc+gQ3HSTLjkYMkR9RgWBIIlUNMFBWwMLRGQJKkqPOOe+\n8K7dBQwVkdVAXeDvCbW+GETgb3/TsD+vvw5vvum3RYYRmZEjYfZsdeTo18LN4ki24KCfAKdGuH8N\n0CmeNpaHnBxdOzVkiLq4OPdc9cdjGEFiyRIYNkyPX3hBZ/WCQpBaUinLjTdCt27w3Xe65iQgE6qG\nAegWrl/9Cvbvh+uvh759S78nkZhIJYC0NHjpJe32jR+vx4YRFIYNg88/h1atitwDBwkTqQTRvHmR\nd4Rbb4XFi/21xzAA/vMfnc2rUgVefVW/SIOGiVQCufJK9YleUKBuL8yLp+En69erEzuABx6AjqUv\nq/QFE6kE85e/wI9/DKtXw+DBNj5l+EN+Plx8sY6T9uwJd93lt0WRMZFKMNWq6c7yH/1IlyU884zf\nFhmVDed0gHzhQmjaVMdJg7BoMxImUj7QqlXR4Pntt8OUKf7aY1QuRo6EV16BY46Bt9+GunX9tqhk\nTKR8YsAAuPdeXeV72WU2kG4kho8+gjvu0OPRo+G003w1JypMpHzk/vt1fUpeHlxwgTq8N4x4sWQJ\n9O+vX4y//71O3iQDJlI+IgJ//7uGw/rmGw0yunu331YZqciaNXDeebBrF1xyiQZXSBZMpHwmMxPe\negtattQuX//+GpnDMGLF5s3Qu7fm3bvreqggD5QfiYlUAKhbVwfPjztOF9ddeqluUTCMirJzp7ag\n1qxRrwZvvw1ZWX5bVTZMpAJCixbw8ccqWO++C5dfDgcO+G2Vkczs2aP78JYs0Rnl99+HmjX9tqrs\nmEgFiFNPVb/StWppF/BXv4LCQr+tMpKRXbvgpz+F//4XGjTQFvpxx/ltVfkwkQoY7durUNWsqYs9\nr77ahMooG9u3q7PFWbOgUSNsjNYjAAANiklEQVSYNg2aNPHbqvJjIhVAOnaEDz7QVemvvaazMXv3\n+m2VkQx8+626BVqwAJo1gxkzdFImmQmMSEUTHFREuocFBl0sIvtEpJ93bbSIrA271i7xtYgdZ56p\nTfQ6deCdd4rcEBtGJL75Rp0qLlsGJ5+sApWT47dVFScwIkUUwUGdc9OcFxgUjbO3F/hPWJE7XVHg\n0KRfw33mmTBzpjbZZ8/W9VRff+23VUYQWbxY/15WrNBV5NOn61hUKhAkkSprcND+aGy9lO4ItW6t\nAtW2LSxfDl266DelYYSYPBm6dtUdC2edpWNQxx/vt1WxI+mCg4Zdnwo86Zx71/s8GugCFOC1xJxz\nBRHujXtw0FiTl1eF4cPbsnRpLY45ppBhw5bTtWv5+n9+1yVWpEo9oHx1cQ4mTGjECy80wzmhd+/N\n3HHHCqpW9fd/utIHB/Wu1QO2ABlHnBMgE22J3RuNTfEKDhoP8vOdGzBAAzeCc3/4Q3yDNwadVKmH\nc2WvS0GBc4MHF/0tjBjh3KFD8bGtrMQ6OGhCo8U453pFuiYi34pIPefcppKCg3oMAN5yzv2w3NE5\nt8k7LBCRfwB3xMToAJGVBePG6ezf3XfDiBE6i/Ovf0HtiG1OI9X46iu44godBsjKgrFjdTtVqhKk\nMalSg4OGcQXwWviJsOjHgo5nfRYHG31HBO68U5co1K0L//63itbSpX5bZiSCSZN0Ld3s2dCwoc7g\npbJAQbBEKprgoIhIDtAImH7E/a+KyDJgGZANjEiAzb7Rq5e2otq3131ZnTqpM7NDh/y2zIgHBQVw\n223Qrx/s2KEeMxYvhjPO8Nuy+BMYkXLObXPO9XTOtfTy7d75Bc65wWHl1jnnGjjnDh1xfw/n3KnO\nubbOuSudc3mJrkOiycnRVcXXXqt/xLffrrvd16/32zIjlqxapbN2Tz8NGRnw1FPaogq6R81YERiR\nMspHtWrqivjtt3Vv1tSpugfw1VctyEOyU1gIjz+u654WLVJ/5LNmwW9/q93+yoKJVIpw0UW6furC\nC3Vz6ZVXqlvizZv9tswoD0uW6OLM3/1O/YtddRV8+mnl6N4diYlUCnHCCdoNeOklqFFDNyifdJJ2\nE2yTcnJQUAD33KOTIQsXQuPG6mLln/+EY4/12zp/MJFKMUR0jGrpUvj5zzUA6W236Tfw7Nl+W2dE\nwjl47z2NyThihH6pDBkCn32mTusqMyZSKUrTproxedIkddOxeLEOvl57LezYkeG3eUYYS5fCnXee\nxgUX6N67Vq10acGzz6onjMqOiVQKI6KeGb/4AoYPh6pV4eWX4Ze/PJN77lHXsoZ/bN4M112ny0gW\nLqxDrVrwxBMqWj/5id/WBQcTqUrAMcdoF2LZMg2dtW9fOiNGaGvr4YfVzayROLZt05iLLVvq+GFa\nGlxyyQZWr4ahQzU4h1GEiVQlolUr7QI+++wiunXTltTvfw/Nm+vgujnWiy+bN+tugSZN4KGHNN5i\n37467nTrrasrzbqnsmIiVQk55ZTvmTpV3RSfcYZ6c7ztNp1JuuceW7YQa77+Gm65RRff/vnP2nL9\n2c/UV9ikSToDa0TGRKqSIqJba+bO1aAPHTtqN2TECP2mv+Ya2w9YEZxTv06XXaYt1b/+VZcXXHyx\nbmeaMgXOPttvK5MDE6lKjojuB5s3TyOLXHKJhtIaPVqnw3v21NXr1hWMjh07dA9l69bQowdMmKD7\nKX/5Sx0TnDhR498Z0ZNQVy1GcBFR745du8KXX+oY1d//rttspk7VqfABA2DgQC1TmbZllEZhobaa\nXn0Vxo8vikBdv77O3g0erB4LjPJhLSnjKJo3h7/8Rd3RPvccdO4Mu3eraJ1zjgYyvf9+XXtVWfcH\nFhaqeN9wA9SrB336wJgxKlB9+mgX+quv9OdkAlUxrCVlRKRWLbjpJk3/+5/+E44dq65hHnhAU6NG\nul+wb18NpZTK0+e7d0Nurm5TefNN+C7MLeNJJxW1NJs3983ElMREyoiKk0/WNVUjRmg4+Ndf13Dw\n69dra+u553S/YJ8+KlZdu+ru/fR0vy0vP4WFMH++zoJ++CHMmXP4HsgWLXRgfMAA9TxhXeD4EBiR\nEpFLgfuB1kAn59yCCOXOA/4CpAMvOedCzvGaAuOAOsAi4Crn3P4EmF6pSE9XIerTRweEFyzQtVeT\nJ+ts4MSJmkDHsc46S1dPd+2qK6tr1vTX/kg4p3HrFizQNH++znzu2lVUJi1NPRP07q2TDe3bmzAl\ngsCIFOru9xLghUgFRCQd+CvquXMDMF9EJjvnvgAeBZ5yzo0TkVHAtcDz8Te78pKWph5BO3XSxYlf\nfaVujWfO1LR2rX7+4IOiexo31vBc4alVK6hePTE2OwebNsHq1Zq+/FLFdf58XS92JC1aqCj17g3d\nu2sX2EgsgREp59xyACn5q6kTsNo5t8YrOw64SESWo8FCf+mVG4O2ykykEkiTJnD99ZpAWyYzZ+rS\nhk8+0T2EX3+tacqUw++tVUsHmMNTvXra8qpeXbuSoVS9ugrk5s1ZrF6tXbBQys/X9V5btx6etmxR\n0fzyy8jLKWrV0vVioXTGGSqqhr8ERqSipAEQ7hx3A9AZqIuGxCoMO58i8VuTlwYNdMzmssv0c2Gh\nisRnnxWlZctUPHbu1PRZmcJnnFkuu7KzdXC7RQtNJ5+sgtSsmXXfgkhCRUpEPgJOLObScOdcSdFh\nfnhEMedcCecj2REeHJTc3NxSX5yXlxdVuWTA77rUrQvnnqsJtAv2/fcZbNmS+UPaurUq27Zlkp+f\nTn5+Ovv2pf1wnJ+fjnNCWtpBqlQR0tMd6emOtDRH1aqHOPbYA9SseYBjjz08nXDCPho02EeNGkd7\nAFy/3l/f8H7/TmJJrOsSmLh7UbIBjRQToiGwEdgK1BKRKl5rKnQ+kh0vAi8CdOzY0XXr1q3UF+fm\n5hJNuWQgVeqSKvUAq0tJJNtizvlASxFpKiJVgcuByV401GlAKAJZaXH7DMNIEgIjUiJysYhsALoA\n74nIB975+iIyBcBrJd0CfAAsByY45z73HnEXMFREVqNjVH9PdB0Mw4g9gRk4d869BbxVzPmNwPlh\nn6cAU4optwad/TMMI4UITEvKMAyjOEykDMMINCZShmEEGhMpwzACjYmUYRiBRlxl9VrmISJbgK+i\nKJqNLhpNBVKlLqlSD6icdWninDuutEKVXqSiRUQWOOc6+m1HLEiVuqRKPcDqUhLW3TMMI9CYSBmG\nEWhMpKLnRb8NiCGpUpdUqQdYXSJiY1KGYQQaa0kZhhFoTKQiICKXisjnInJIRCLOVIjIeSKyQkRW\ni8jdibQxWkSkjoh8KCKrvLx2hHIHRWSxlyYn2s5IlPYzFpFMERnvXZ8rIjmJtzI6oqjLIBHZEvZ7\nGOyHnaUhIi+LyHciUqwvVVGe9uq5VEROL/fLnHOWiklo1JqTgFygY4Qy6cCXQDOgKrAEaOO37cXY\n+Rhwt3d8N/BohHJ5fttanp8xcDMwyju+HBjvt90VqMsg4Fm/bY2iLucApwOfRbh+PvA+6jX3TGBu\ned9lLakIOOeWO+dWlFLsh8AQTsNnjQMuir91ZeYiNDgFXt7PR1vKSjQ/4/D6vQH0lFIievhEsvy9\nlIpzbgawvYQiFwH/dMoc1HNuvfK8y0SqYhQXGCKIASBOcM5tAvDy4yOUyxKRBSIyR0SCImTR/Ix/\nKOPUMeIu1PFh0Ij27+UXXhfpDRFpVMz1ZCBm/xuBcXrnB3EMDJFwSqpLGR7T2Dm3UUSaAVNFZJlz\n7svYWFhuovkZB+b3UArR2PkO8JpzrkBEbkRbiD3iblnsidnvpFKLlItfYIiEU1JdRORbEannnNvk\nNbm/i/CMjV6+RkRygfboGIqfRPMzDpXZICJVgGMpuSviF6XWxTm3Lezj39Cgt8lIzP43rLtXMYoN\nDOGzTcUxGQ1OARGCVIhIbRHJ9I6zgbOBLxJmYWSi+RmH168/MNV5o7cBo9S6HDFu0xf15Z+MTAau\n9mb5zgR2hYYcyozfswRBTcDF6LdBAfAt8IF3vj4w5YhZjJVoi2O433ZHqEtd4GNglZfX8c53BF7y\njs8ClqEzTsuAa/22u6SfMfAg0Nc7zgJeB1YD84Bmfttcgbo8DHzu/R6mASf7bXOEerwGbAIOeP8n\n1wI3Ajd61wX4q1fPZUSYIY8m2YpzwzACjXX3DMMINCZShmEEGhMpwzACjYmUYRiBxkTKMIxAYyJl\nGEagMZEyDCPQmEgZhhFoTKSMpEFEqovI/0RknohkhJ3v4zknHOKnfUZ8sBXnRlIhIu2BOcBTzrm7\nReR4YCkwzznX11/rjHhgImUkHSJyO/AE0Ae4AzgV+LFzLlUiABthmEgZSYfndfM91M9SVaC3c+5j\nf60y4oWNSRlJh9Nv1rFAJrDEBCq1MZEykg4ROREYCSwCfiwit/lskhFHTKSMpMLr6o0B9gO9UbF6\nVERO89UwI27YmJSRVIjI/6Ehuno456Z7Hi7noF2/js65fF8NNGKOtaSMpMFbfvAn4GHn3HQAp6Gh\nrgBygCf9s86IF9aSMgwj0FhLyjCMQGMiZRhGoDGRMgwj0JhIGYYRaEykDMMINCZShmEEGhMpwzAC\njYmUYRiBxkTKMIxA8/8uIxybK5bDNQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# define the cylinder of untit radius centered at (0, 0)\n", "R = 1.0\n", "x_center, y_center = 0.0, 0.0\n", "theta = numpy.linspace(0.0, 2 * math.pi, 100)\n", "x_cylinder, y_cylinder = (x_center + R * numpy.cos(theta),\n", " y_center + R * numpy.sin(theta))\n", "\n", "# plot the cylinder\n", "size = 4\n", "pyplot.figure(figsize=(size, size))\n", "pyplot.grid()\n", "pyplot.xlabel('x', fontsize=16)\n", "pyplot.ylabel('y', fontsize=16)\n", "pyplot.plot(x_cylinder, y_cylinder, color='b', linestyle='-', linewidth=2)\n", "pyplot.xlim(-1.1, 1.1)\n", "pyplot.ylim(-1.1, 1.1);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Discretization into panels" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A panel, which represents a source sheet, is defined by its two end-points (`xa`, `ya`) and (`xb`, `yb`) and its strength `sigma`. We'll also need its center point (`xc`, `yc`) and its length for further calculations. The orientation of the panel is defined by the angle between the $x$-axis and its normal in the counter-clockwise sense.\n", "\n", "What information do we need to compute on each panel? First of all, we will need the strength of the source sheet that will lead to the correct streamlines. In addition, we'll also want the tangential velocity (the normal velocity on the body is zero for an inviscid flow) and the pressure coefficient.\n", "\n", "In this lesson, you'll really appreciate having learned about classes. It will make the code so much easier to manage. We create a class named `Panel` containing all the geometry data related to one panel. With a start- and end-point, the class internally calculates the center-point, length and normal vector. It also initializes to zero the source strength, tangential velocity and pressure coefficient. (These will be updated later.)" ] }, { "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, length, and angle \n", " (with the x-axis) of the panel.\n", " Initializes the strength of the source-sheet, the tangential velocity,\n", " and the pressure coefficient 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\n", " self.xb, self.yb = xb, yb\n", " \n", " self.xc, self.yc = (xa + xb) / 2, (ya + yb) / 2 # control-point (center-point)\n", " self.length = math.sqrt((xb - xa)**2 + (yb - ya)**2) # length of the panel\n", " \n", " # orientation of the panel (angle between x-axis and panel's normal)\n", " if xb - xa <= 0.:\n", " self.beta = math.acos((yb - ya) / self.length)\n", " elif xb - xa > 0.:\n", " self.beta = math.pi + math.acos(-(yb - ya) / self.length)\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": [ "To store all the discretization, we create a NumPy array of size `N_panels` where each item in the array is an object of type `Panel`." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAF8CAYAAAAdLMb/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3Xd4FMUbwPHvJLn0kIQ0EkCIoTeF\n0EEpIkgHlSLSBBEVUbEBFpqiKCo/C70IFqQrSBFBKYIgEBRRIUgvoaZALsklubv5/bFJSLmEJORK\nwnye5x5yu7O3b47cvrfvzs4IKSWKoiiKYg1O9g5AURRFKbtUklEURVGsRiUZRVEUxWpUklEURVGs\nRiUZRVEUxWpUklEURVGsRiUZRVEUxWpUklEURVGsRiUZRVEUxWpUklEURVGsxsXeAdhbYGCgrFq1\narG2TUpKwsvLq2QDKgEqrqJRcRWNiqtoympcUVFR16SUQbdsKKW8ox+RkZGyuLZt21bsba1JxVU0\nKq6iUXEVTVmNCzggC3GMVeUyRVEUxWpUklEURVGsRiUZRVEUxWpUklEURVGsRiUZRVEUxWpUklEU\nRVGsRiUZRVEUxWpUklEURVGsRiUZRVEUxWocLskIIRYJIa4IIf7OZ70QQnwqhDguhPhLCNEo27oh\nQoj/Mh5DbBe1oiiKYonDJRlgMfBQAes7A9UzHk8BswGEEOWBiUAzoCkwUQjhb9VIFUVRlAI5XJKR\nUu4E4gpo0hP4MmP4nL2AnxAiFOgEbJFSxkkp44EtFJysFMXhmM3aQ1HKCqGNc+ZYhBBVgfVSynoW\n1q0Hpkkpd2U8/xkYC7QF3KWU72QsfwtIkVJ+aOE1nkI7CyIkJCRy2bJlxYpTr9fj7e1drG2tScVV\nNCUdl5Sg17tw9aobV6+6ce2a9m98vI6UFGcMBmdSUm4+sj9PTXVCSoEQEmfn/B8uLpJy5dIpVy4d\nX9+bj9zPAwLSKFfOWGK/G9w5/48lpazG1a5duygpZeNbtSuNQ/0LC8tkAcvzLpRyHjAPoHHjxrJt\n27bFCmT79u0Ud1trUnEVTXHiSk2F6Gj4+2/49184exYuXIDz57WHTgcVK0KlStqjShVo0gR8fMDb\nW3t4ed38OfPh4QFOTmA2C37+eSetWrXBaCTPIzUV4uI8iI2Fa9dyPs6cufnz+fPg7AzVquV8RERo\n/wYHg7D0ySnh98sWVFxFY6u4SmOSOQ9Uzva8EhCTsbxtruXbbRaVUiaZTHDihJZMsj9OnYK774Z6\n9aBOHWjf/mZCqVhRSya3w8kJdDqJp+ftvY6UEBsLx4/ffPz0k/bviRNgMGjJ5t57oXFjLRE2aADu\n7re3X0XJVBqTzDrgOSHEMrSL/NellBeFEJuBd7Nd7O8IjLdXkErpFBsLu3fDr7/Crl1w6BCEhmrJ\npF496N0b3noLatQANzd7R3trQkBgoPZo3jzv+oQE+O8/+OMP2L8fFizQztBq1dKSTuajXj1wdbV9\n/Erp53BJRgjxLdoZSaAQ4jxajzEdgJRyDrAR6AIcB5KBJzLWxQkh3gb2Z7zUFCllQR0IlDLOlKQn\n/UoMuuAwnL3y1p6l1MpcW7aEsGyZlljOn9cOxq1bw3vvaQdYe5TT9Xo9586ds3o9389PO3tp0gSe\nekpblpICf/2lJZ3ffoNPP4WTJ6FhQ+jQAQIDy9GqlVYStORW77tyZ3G4JCOlfOwW6yUwKp91i4BF\n1ohLKT2k0ci5qa8Qu2IRwsUFaTQS0HcYld/4kOt6F378ETZsgO3bIT0datUKoHdv7SDboAG42PFT\nYTQaGTNmDAsXLgRg5MiRDB8+nBkzZuBio8A8PKBZM+2RSa+HPXtgyxb49NPqvPUWtGmjJZ0HH4Sa\nNQFT/u+7sOebqtiV+p9XypxzU18hduViZKoBmQpXjSH8Mv8y/y75itM3KlK/HnRpDm9Mh9Aw+Ouv\nv7inQQPQQ/Jv9o391U9m8tWGzaSkpmYtWzh/PmkxZ5n+gsXvVjbTzA2adYPud/1Flbsa8Mef8MfP\n8Np0QEDtckepdv0IdV3A3UkPQOzKxQDcNfF/9gtcsSuVZJQyxZSk58ryxZzVVyLE5TweTskEuVym\nnctq2rEa/NC6iKyC5FVwAvBCq73aW7JZsuS0kdRcfSJTUlNZ8t06Bv2xEU+nInYFswIv4Bpa75vK\nQI/MazUpQLbrNhfSq1DefAVWfkHFV95RpbM7lEoySpkgJfz5JyyZkUD52N508/kGgBSzB+4iReum\n6+yM173NcPLwyrFtfHwc/v7l7RB1Tteu63E+txuMpjzrXJydSb6nJSG+9j9QW3q/zClJJP35O5hM\npEsXnJBU1J3hQnpVvrj0JA8suUHPod633VtOKX1UklFKtUuX4JtvYMkSCLoRxQS/wfj4/Je13sMp\nJetnoXOl2qL1eb5Rb9++naYOcB9DqF6PDA4GY0qedWadKy2/2uAQN/VZer9MSXr+ahqGNJnQiZs3\nf1bUneZp37f4eUEylV5/k06dXejXD7p2zb/jgFK2ONywMopyKwYDrFwJ3bppXW3//svEnIfeZ7r7\nffjo/8PZrzy45rzRQ7h7EtDnCYcu2Xh7ezN8+HA8c33d9/T0ZPjw4Q6RYPLj7OVNQN9hCPdcpyrO\nLjgJyYM33mVn2/Z0uvckH3+s3Zz61lta7z6lbFNJRik1/v0Xnn1Wu9lxzhzo2xdO7TvD+LSOuK95\nC4xGgoaOpt6O4wT2G4Zw98DJyxvh7kFAn6FUfiPPCEMOZ8aMGQwbNgwPD4+sx7Bhw5gxY4a9Q7ul\nym98SECfoTne98ABT1Ft8UZ0IRVJ+3sv9y5tzHcvfMmWLZIbN7Ru0T16wMaN2o2vShkkpbyjH5GR\nkbK4tm3bVuxtraksxWU2S/nTT1I+9JCUISFSTpwo5enT2rrYtUvln/cGyqgInTzUvLK8vmNzjm2N\n+kSZcjJaGvWJJR6XtSUmJsovv/xSJiYWHLs93Or9svS+p8fHyhPP9ZdREToZFaGTJ0b1k+nxsVKv\nl3LhQikbN5ayalUp331XykuXrBOXvZTVuIADshDHWHUmozgkgwEWLdLuWxkzBvr0gdOnYdIkqFT+\nOqdeGszpl4ZgSryOb4fu1NlwkHL3d8zxGs5e3riH13DoEll+vL29qVy5skOXyPJj6X138StP+KdL\nqfLBApy8vEn4cQ1HukZi/msbw4ZpN36uXKkNdVOrFjz2mNaRQyn9VJJRHMrVqzBlClStqh10Pv4Y\nDh+GYcO08bT0+3dxpFtj4tctw8nDk7vemc3ds1fhUj7Q3qErtyCEIODhwdRetx+vhs1Jv3yB/wY/\nxPlp4zCnptK4sTaszalT2kgLXbtC9+6wd6+9I1duh0oyikM4dw6eflobE+zcOfj5Z9i0SbubXAiQ\n6elc+Ogtjj3egbQLZ/CsH0mtdfsI7D8cUdRhhBW7cqsSQY1vfyH0+bfAyYkrCz4m+tHWpPz3L6AN\ndfPyy9pZTZcu0L+/NrLAtm1aV3WldFFJRrGry5fhxRe1UYD9/bXBGefPh7p1b7YxnDpGdN/7uTz7\nfZCSCs+Oo+aKnbiH17Bf4MptES4uhD7/FjWXbce18t2kHDnE0V7NufrVbGRGJnF3h2ee0QbwHDgQ\nRo7UxpTbuFElm9JEJRnFLuLj4Y03tGHypYR//tEGpAwOvtlGSsm1ZQs52qMpyYejcK1YhRpLfybs\npSkIdZNFmeDVsBm1f9hPwKNDkKkGzk1+gRMjepF+7XJWG50Ohg6FI0dg9GgYOxYiI2HdOpVsSgOV\nZBSb0uvh3Xe1stilS3DwIHzyCVSokLOdMe4aJ595lLNvPoM5JRn/Ho9Re/0BvJu0tk/gitU4e/tQ\nZdp8wj/7Fmdff25s38SRLo24/suGnO2ctdLZoUMwYQKMH6/N4/PHH3YKXCkUlWQUm0hNhdWrK1K9\nujaM/K5dsHChdlNebjd2/sS/XRtxfesPOHmXo+rHSwj/eAnOPr62D1yxGf/Oj1B7fRQ+LdphjLvK\niad6c3bCc5hTknO0c3KCXr20ZNOvH3TuDE88AdeuqQlvHJFKMorV/fQT1K8P+/eXZ9MmWLYsY2j4\nXMypBs69/RLHh3XDePUS3k1aU3tDFOV7FDj7g1KGuIZWotqSTVQc9z5C58q1pfM40rMpyX/nPV1x\ncdE6i0RHQ0gIDB/ehMmTISnJDoEr+VJJRrGa8+e1+1ueflrrijxt2mHuvddy2+Sjf3G0dwuuLvkc\nXFwIe/ltqn+9BbeKFk51lDJNODkR8uQYaq7ZjXu12qSePEZ0n9ZcmjsdaWFYAF9fmDYN5s6N4uhR\n7T6bL78Es9kOwSt5qCSjlLi0NPjgA63HWJ062kX9bt0st5VmM5cXfUJ075YYjv2DW3h1aq78lQrP\njEU4O9s2cMWheNa+h1rf7yVo0LPI9HRipr/Bf4M6kRZjecCzChUMfPstrFgBs2dD06baNT/FvlSS\nUUrUtm1actm2TbuJbvJkbaZFS9Iux3D8iW5cePdVZHoagf2fpNbafXjVj7Rt0IrDcnL3oPLE/xGx\ncB0ugSHo9+3kSNdI4jasyHebFi20aaOff167XjNunDaltGIfKskoJeLiRRgwQOtqOnWqdi9DtWr5\nt4/f/B1HujYicfdWXPwDuXv2Ku56ZxbOnl75b6TcsXzbPETtDVH4tu+KKfE6p18YyOlXnsCUeMNi\neyFg8GCtk8mpU3DPPbBzp42DVgCVZJQSsHKldvZSpYo2UnLv3tqH3BJTkp4z457i1Kh+mBLiKHdf\nR2pvPIjfgz1sG7RS6ugCgrl77hoqvz0T4e5B3PffcKR7Y/RR+c+ZHRICy5dr5dsBA7RRvG9YzkuK\nlagkoxRbQoJ2J/abb8L69drNlF4FnIg4nTzK0R5NiV21GOHqRqW3PiZi4Tp0QRXy30hRshFCEPTY\nCGqv/R2Pug1JO3+aY4+1J2bGRDAa892uVy/4+29IT4d69WDDhnybKiVMJRmlWH75RRsh2c9Puxmu\nSZP820qjkYufT8Vr2ouknjmOR6361Pp+D8FDnkM4qT9BpejcI2pRc+WvhIx8FaTk0sz38Hp/DIbT\n/+W7jZ+fNmTRF19o12sGDtS+KCnWpT7hSpGkpGhD7w8erH1gP/+cAudtTz13imOPd+Di/yYjzGaC\nh71IzdW/4VGjnu2CVsokJ1dXKr46lepfb0EXWhnnU9Ec7dGUayu+yBr/zJIHHtCu1fj6QqNG2jQD\nivWoJKMU2h9/aEOwX7ig3W3dqVP+baWUxH73NUe6NSYp6jd0IWEkvTSNSq9/gJObm+2CVso8n2b3\nU3v9AdKbtsWcnMTZ10dyclRfjPGx+W7j5QUzZ8L06dqUAjNmqHHQrEUlGeWWpNRupuzUCV5/XbuQ\nGhCQf3vj9XhOvziQM68Ow5yUiF+nXtReH4WpTiPbBa3cUVx8/UkZMZ4qH36Bk3c5rv+0liNdG3Fj\n19YCt3vkEfj9d20Uih49IDb/vKQUk0oySoGSkrRZCpcu1coKjz+et+eYXq/n2LFj6PV6Evfu0CYV\n27ASJ08v7npvHuGfL8fFv4CspCglQQgCej1O7fUH8IpsSfqVixwf2oXzU1/FnGrI8XeaXXg4/Pqr\nNlJAw4bauHpKyVFJRsnX8ePQvLl2zWXXrryDWRqNRkaPHk1wcDCRkZEE+fsz4qH2pMScxfOeptRa\nt5/APkPVpGKKTblVqkqNb7YSOmYSODsTs+h/DKpZmaDAQCIjIwkODmb06NEYs/VGc3XVSmezZ8Oj\nj2o9JdWwNCVDJRnFog0boGVL7b6ChQu1CaRyGzNmDIsWLSIlJQW9Xo/BaGTtDTNzwhpRc9k23KsW\ncDemoliRcHEhdNTr1Fyxk/+l+rDm7DUMqano9XpSUlJYtGgRY8aMybNd165w4IA2K2vXrnD9uh2C\nL2McLskIIR4SQkQLIY4LIcZZWD9DCPFnxuOYECIh2zpTtnXrbBt52WA2a0PBjBwJ33+vzUxo6URE\nr9ezYMECkpNzDsNukLAs6jBJqak2ilhR8icjavNdbAqGXBf1k5OTWbhwYZ7SGUClSloX/WrVtC9a\nJ0/aKNgyyqGSjBDCGZgJdAbqAI8JIepkbyOlHCOlvFdKeS/wGbAm2+qUzHVSSnULeRElJEDPnrBl\ni3b9pWXL/Nue+ecwTqZ0i+ucnZ2JiYmxUpSKUngxMTG4uLhYXOeMzPfv1MUFPvtM+5LVqhXs3m3N\nKMs2h0oyQFPguJTypJQyDVgG9Cyg/WPAtzaJrIw7eRKaNYOqVbVvcaGh+be9vn0TiaMewWTMO+w6\ngMlkIiwszDqBKkoRhIWFYbIwPQBAusGAcf4HmJLzn4Dmuee0mzd794avv7ZWlGWboyWZisC5bM/P\nZyzLQwhRBQgHfsm22F0IcUAIsVcI0ct6YZYtBw9C69bwwgvatzfXfCYYNBtSODfpBU482RPXhGv0\nqV4Jz1xDLHt6ejJ8+HC8vb1tELmiFMzb25vhw4fjmeuOYQ9XHb18XTB89yVHezQh6a8D+b7GQw9p\no4pPmKANoaQ6BBSNKOjOWFsTQvQBOkkpn8x4PghoKqUcbaHtWKBS9nVCiDApZYwQ4m605POAlPKE\nhW2fAp4CCAkJiVy2bFmx4tXr9Q55MC1KXFFR/rzzTm3GjDnG/fdfy7ed09kTeMx/D+eLZ5HOLqT2\nGkpKh17MnD2HjRs34uTkhNlspkuXLowaNQpnC3PBlIX3y5ZUXEWTX1wmk4mZM2fm+Tt9rkdnfL6Y\njvOF00hnZ1K7DyKtSz9wsjyPUUKCjjffrEdgYCrjxh3F3b1w2aa0vV+F1a5duygpZeNbNpRSOswD\naAFszvZ8PDA+n7Z/AC0LeK3FwKO32mdkZKQsrm3bthV7W2sqbFxLl0oZHCzljh35tzGbTPLS/I/k\nwVqeMipCJ//uWE8mHT6Yo01iYqKMjo6WiYmJJRKXram4iqa0xmXp79RkSJFn335JRkXoZFSETh7t\n11Yazp3K9zVSUqQcOFDKJk2kvHq1ZOKyl9uNCzggC3Fcd7Ry2X6guhAiXAjhCvQH8vQSE0LUBPyB\nPdmW+Qsh3DJ+DgRaAf/aJOpSaMYMeO012LoV7r/fcpu0i+c5PqQzF6aNQ6anE/j409T+/nc86zXM\n0c7b25saNWo45Lc1Rclk6e/Uyc2dym9+RLUvNuASVIGkA7s50q0xcWuXWnwNd3dtaucHHoC2bbV5\nlJSCOVSSkVIageeAzcARYIWU8h8hxBQhRPbeYo8ByzKyaabawAEhxCFgGzBNSqmSTC5ms5Zc5s3T\neszUr2+5XfzGVRzpFkninm24lA8iYt533DX5U5w8ChgNU1FKqXL3PUidDQfxfbAHZv0NTr88lFNj\nBmG8kXeYZiG0mzUfewzatIFz5yy8oJLFct8+O5JSbgQ25lo2IdfzSRa2+w3I55CpgDbdxrBh2p38\nu3ZZHn/MpE/k3JQXiVvzFQDl2namyrR56AJDbBytotiWS/lA7p61ktgVX3D+nZeI/2E5SVF7qPLh\nF/g0vS9P+zfe0EbDuP9+rSIQEWGHoEsBhzqTUazHZNKG579yRftAWEow+oN7OdK9CXFrvkK4uVN5\n0qdEzP9eJRjljiGEILDfMGqt24dn/UjSYs7y3+MduDD9DcxpaXnajxkDY8dqpbOjR20fb2mgkswd\nwGTSzmCuXoXvvss7/4s0Gon5ZArHHmtH2rmTeNS+h1rf7yVo4NNq3DHljuQeXoOaK3ZS4dlxIASX\n507nWN/7MZyMztP26adh6lRo316bAkPJSSWZMs5shqeegrNnYe1ayHVbC6lnTnCsfzsuffYOmM0E\nj3iJmqt24VG9juUXVJQ7hNDpCHtpCjW+2YprxSok/32Qoz2bcfXb+XkmRRs8GD75BDp2hH377BSw\ng1JJpgyTUhvg8tgx+OGHnGcwUkpiV3/JkR5NSPrzd3QhFan+5WYqjZ2mJhVTlGy8m7Sm9voD+Pd4\nDHNKMufeGsXJpx8hPfZqjnZ9+sCCBdCtG/z5p52CdUAqyZRRUmrzmB86BBs3QvbexcaEOE6Nfowz\nY5/EnKTHr/Mj1N4QhU+LtnaLV1EcmbOPL+EfL6HqjC9x9vHl+s/rOdK1Edd3bM7Rrnt3mDVLG8H5\nRJ7bwO9MKsmUQVLCyy/D3r3w44/g43NzXeKebRzpGknCj2tw8vKmygcLCP90KS5+5e0XsKKUEuW7\n96fW+gN4N7kP47XLnBjenXNTxmA2pGS1efRRmDhRK52p+2hUkimTFi4MZ/t2+Okn8PXVlplTUzn/\n3lj+G9SJ9MsX8GrYnNo/HCDg4cHq4r6iFIFbxSpU//onwl55B1xcuPrlTI72bkHykZtX/Z96Suts\n89BDoNc73J0iNqWSTBkzbx5s3x7ETz+Bv7+2LOXYP0Q/2oorC2eAszOhL0ygxre/4HbX3fYNVlFK\nKeHsTIWnX6PWql24hVfH8N+/RD/cissLZiAzRtB8/XVo1w5ef70euaZduqOoJFOGbN6sjRQ7bdph\nAgO1i/tXvpzF0d4tSDnyF66V76bmsu2Ejn4Tkc8cG4qiFJ5nvUbUWruPwMdGINPTuDBtLMeHdiHt\n0gWEgI8/hgoVDPTrB+mWp18q81SSKSP++gsGDYJVq6BSpRTSr17ixJM9OT/lRWSqgYBHh1D7h/14\nNWxm71AVpUxx9vTirrdncvec1bj4B5L42y8c6RZJ/I9rcHKC116LxmyG4cPvzGkCVJIpA2JitF4t\nn3yizQvj8ucejnSN5MaOH3H29Sf882VUmTYfZ2+fW7+YoijF4tehO7U3HqTc/Z0wJcRx6rn+nBk3\nAhdjEitXasM5TZpk7yhtTyWZUk6v1xLMU09B355JnH1rFJ6fT8QYdxWflu2pvT4K/4cetneYinJH\n0AVVIGLhOipNmIFwdSN21RK8Jz+LjP6d776DJUtgzZpbv05ZopJMKWYywYAB0KABvNjjIEd7NePa\nt/ORLjoqjnufaos34hpayd5hKsodRQhB8OBR1Pp+Lx616uN0NYbo/m0xr3iHNSuNjBwJf/9t7yht\nRyWZUuzllyFZb+Kdxh8Q3ac1qSeP4V6tNklvfErIk2MQTuq/V1HsxaNGXWqu/o3Ujo+CycTFT6bg\n/eEDfP7WSXr1grg4e0doG+ooVEotXQp715/l08BOXJ7xJhiNBGV8ezJXVmOOK4ojcHJzI7XvU1Rb\nsgldSBhJB/dQc2ETXmz4Jf37SYxGe0dofSrJlCKmJD2GU8f452ASa15ZzizPSAwHd+ISGELEwnVU\nnjADJ3ePW7+Qoig2Va7VA9ReH4Vfp96YkxJp+ceTDLgygAkvx2d9rk1JenuHaRXqZolSQBqNnJv6\nCrErFmEQXvwe25zXvTZAMvi278pd781DFxBk7zAVRSmAi38A4Z8vI271l5yb8iINklYTtu435n1X\nnxb+vyGNRgL6DqPyGx+WqfvY1JlMKXBu6ivErlyMTDWQkgT3e20AwKtZG+6eu0YlGEUpJYQQBDw6\nhFrr9uMSGEKg80Wa6LZw9boHMtVA7MrFnJv6ir3DLFEqyTg4U5Ke2BWLkIZkDGYP/JxjMUntvy35\n0D7MyUl2jlBRlKLSBVXAeCMBAAEEuFwlzaxDGpKJXflFmSqdqSTj4NKvxCBcXLiQVgV3pxTSpCvO\nQrttWDg7k34lxs4RKopSVOlXYnDS6TKeaROgnUuvBpS9z7VKMg5OFxxGUpobf6c2BcBV3JxnXJpM\n6ILD7BWaoijFpAsOQ2Z0LRNCm57jLtfj7E1qX+Y+1yrJODhnL28+dVtFK6+ckyMJd08C+jyBs5d3\nPlsqiuKonL28Ceg7DOGuTVcrBOhEOr+mdEV2Hl2mPtdlpwtDGbVpEyRfT8Xb+QYIgfDwBLOZgD5D\nqfzGh/YOT1GUYsr8/Mau/ALMZmRaKv0qLmPS0V9ZK7XEUxaoMxkHdv26NibZ6+20wY6Ch4+h9trf\nafD7Be6a+L8y1c1RUe40wsWFuyb+jwa/X6D615vByYkqpj+JPX+DefPsHV3JUUnGgb3yCnR9yIjX\nP98DUL5Hf9zDa5SpU2lFudM5e3nj3agl3k3vRxrTmTN8HW+8AdHR9o6sZKgk46C2bNGmT57U71eM\nsVdwq1oNj9r32DssRVGsxL/LowB4/LmaKVNg8GBtENzSTiUZB5SYqJXJ5s2DtJ2rAfDr/AiirBRp\nFUXJw69TL3ByInH3Vp58LB6djjJRNlNJxgGNG6fNDd7xASMJP34HgH/nR+wclaIo1qQLCNZKZunp\nJP7yA7Nnw8SJcPmyvSO7PQ6XZIQQDwkhooUQx4UQ4yysHyqEuCqE+DPj8WS2dUOEEP9lPIbYNvKS\nsX07rF2rzQ2u3/8rxrirqlSmKHeIzJJZ/KbV1K8PQ4dq12ZLM4dKMkIIZ2Am0BmoAzwmhKhjoely\nKeW9GY8FGduWByYCzYCmwEQhhL+NQi8RqakwYgTMng1+ftofGqhSmaLcKbKXzIzX45kwAXbuhF9+\nsXdkxedQSQYtORyXUp6UUqYBy4Cehdy2E7BFShknpYwHtgAPWSlOq/jsM6hdW5tOWRpVqUxR7jTZ\nS2bXt/6Atzd88gk8+6z2JbQ0crQkUxE4l+35+YxluT0ihPhLCLFKCFG5iNs6pCtXYNo0+DDj/kpV\nKlOUO1P2khlAz55Qo8bNY0Np42h381mqCclcz38AvpVSpgohngaWAO0Lua22EyGeAp4CCAkJYfv2\n7cUKVq/XF3vb3D76qAbt2pmIiTlBTAy4f/UprsCNOk3YsWOH3eIqSSquolFxFU1ZiUt4B+EtnLj+\n6xa2b/gBvHwYMMCdkSMjCQ+PIizMYJe4ik1K6TAPoAWwOdvz8cD4Ato7A9czfn4MmJtt3VzgsVvt\nMzIyUhbXtm3bir1tdn/+KWVwsJRxcdpzc3q6PNQkTEZF6GTSv3/aLa6SpuIqGhVX0ZSluI4N7Cij\nInTy2qolWcvee0/K7t3tG1d2wAFZiOO6o5XL9gPVhRDhQghXoD+wLnsDIURotqc9gCMZP28GOgoh\n/DMu+HfMWObQpIQxY2DSJPDP6KaQo1RWq4Fd41MUxfb8Mq7DZpbMQDtOHDoEu3fbK6ricagkI6U0\nAs+hJYcjwAop5T9CiClCiB4XFW+lAAAgAElEQVQZzZ4XQvwjhDgEPA8Mzdg2DngbLVHtB6ZkLHNo\na9dq12NGjLi5LH7jKkD1KlOUO1XuXmYAbm4webJ2H520eCHAMTlUkgGQUm6UUtaQUkZIKadmLJsg\npVyX8fN4KWVdKeU9Usp2Usqj2bZdJKWslvH4wl6/Q2Glpmp94GfMgMyxLqXRSMJmbayyzAuAiqLc\nWXQBwfg0a5PVyyzToEEQG6uNzl5aOFySuZPMnKl1WX7wwZvLEvftVKUyRVEslsycnWHqVHj9dTCb\n7RVZ0agkYyfJyfDBB/DuuzmXJ2T8Qfl3eVSVyhTlDmapZAbQq5dWOlu+3I7BFYFKMnYybx60bAn1\n699clr1U5qduwFSUO1p+JTMh4L334K23ID3djgEWkkoydmAwwPTp8OabOZerUpmiKNlZKpkBtG8P\n4eGwaJE9oioalWTsYNEiaNgQGjXKuVyVyhRFyS6/khlopfYpUyAlxU7BFZJKMjaWlgbvv6+d6maX\nfawyVSpTFAXyL5kBNGkCkZGwZImdgisklWRs7MsvoVYtaNYs5/LEfTsxxl/DLby6KpUpipIlq2SW\ncf9cdi+/rN0C4cg9zVSSsSGj8eYFu9yySmXqBkxFUbLJKpn99nOektn994OPD6xfb6fgCkElGRta\nuhTuugtat865XJXKFEXJT0ElMyG0s5mPPrJTcIWgkoyNSKkN1f3663nXqVKZoigF8csc/t9CyezR\nR+HUKThwwNZRFY5KMjaye7d20b9Dh7zrVKlMUZSC+HXsCU5O3LDQy0yngxdecNyzGZVkbGT2bHj6\nae30NjtVKlMU5VYyS2YYjVzfsi7P+iefhJ9+grNn7RDcLagkYwNXr8LGjTBkSN51qlSmKEph+OWa\nMTM7X18YOlSbqtnRqCRjA4sWQe/eN+eLyS4ho8aqSmWKohSkoJIZaCWzL74Avd4OwRVAJRkrM5lg\nzhx45pm863KMVaaG9VcUpQC6gGB8mrfNt2R2113aeIhr1tg+toKoJGNlmzdDYKB2d25uOUplNevn\nbaAoipJNfmOZZRo6FBYvtl08haGSjJXNnm35LAZUqUxRlKK5Vcmse3f46y84c8YOweVDJRkrOnMG\n9uyB/v3zrlOlMkVRiupWJTM3N+jXTxu+ylGoJGNFS5dC377g6Zl3XeLvO1SpTFGUIitMyWzJEu0G\ncEegkowVrVihfauwRA3rryhKcfh17HWzZJYQl2d948baGc3u3XYIzgIXewdQVh07Bpcu5R2nDO6M\nGTCvX7/OtWvXSEtLy1rm6+vLkSNH7BiVZSquoiktcbm6uhIYGIivr68doyp5uoAgfJq3JfG3X7i+\n9QcCHs15A54Q2j15ixdbPv7YmkoyVrJihTamkLNz3nUXt/3IictXqRRRrUyWygwGA5cvX6ZSpUp4\neHhknaklJibi4+Nj5+jyUnEVTWmIS0pJSkoK58+fx83NDXd3dztHV7L8Oj9C4m+/cH7dcmIbtCAs\nLAxvb++s9QMHQr168Omnlsv1tqTKZVZiqVRmNBoZPXo0EV16MfC8kft3H+f555/HaDTaJ0gruXr1\nKkFBQXh6eqpSoGIXQgg8PT0JDAzk6tWr9g6nxHm378YH10w0+3oTkY0aERwczOjRo7OOJWFh0LQp\n/PDDLV7IBlSSsYIjRyA2VrsxKrsxY8awaNEiDEYTyRIM6UYWLVrEmDFj7BOolRgMhhzfqhTFXnx8\nfDAYDPYOo8S99s67rNVDqgR9UhIpKSl5jiW9e6skU2atXAl9+oBTtndXr9ezcOFCkpOTc7RNTk5m\n4cKF6B1tLIjbYDQacXFRlVjF/lxcXMpcpSDzWGIw5ZwOM/expFs32LRJmyzRnlSSsYLly7Wuy9nF\nxMTgbOkCDeDs7ExMTIwNIrMdVSZTHEFZ/Dss7LGkYkWoWtX+vcxUkilh//wDiYnQvHnO5WFhYfl+\nozKZTISFhdkgOkVRSruwsDBMJpPFdbmPJd27279kppJMCdu4UfuPdcr1znp7ezO4ayfcc32x8vT0\nZPjw4eoaxh1g+/btCCHYvn171rK2bdvStm1bq+5DKVu8vb0ZPnw4nrm6jVk6lvToYf8k43CFcyHE\nQ8AngDOwQEo5Ldf6l4AnASNwFRgmpTyTsc4EHM5oelZK2cNmgWfYujX/scpeqx7CdR/B2hRnXNzc\nMZlMDBs2jBkzZtg2SMVhzJo1y94hKKVQ5jFj4YIFiFQDJmDo4wPyHEsaNoSkJO2+vRo17BAoDnYm\nI4RwBmYCnYE6wGNCiDq5mv0BNJZSNgBWAR9kW5cipbw342HzBGMwwG+/gaUvptJoJGnrD7wW5MKZ\nPb8SFRXFlStX+Oyzz9RF8jtYnTp1qFMn95+4Y5BS5riZVnEcLi4ufPbZZ1y5epXverRma1UXpnRs\nnedYIoTWAcCeZzMOlWSApsBxKeVJKWUasAzomb2BlHKblDKzi9ZeoJKNY8zXnj1Qty74+eVdlzVW\n2d01CGrUjBo1aqgSWSlz6NAhevfuTUBAAB4eHtSsWZP33nuP5557jpCQENLT03O01+v1+Pj4MH78\n+HxfM3e5LLPctW7dOl5++WUCAwMJCgpi4MCBJCQk5Nj26tWrDBgwgHLlyuHn58fgwYPztMm0Zs0a\nmjdvjqenJ35+fvTp04ezuebqrVq1KgMHDmTRokXUqlULV1dXNmzYUMR3SbElb29vGvUfgqeTID5j\nVPfcuneHdXnH0rQZR0syFYFz2Z6fz1iWn+HApmzP3YUQB4QQe4UQvawRYEG2boUHH7S8Ll4N61+q\n7du3jxYtWnDixAlmzJjBhg0beOmllzh//jzPPvssV65c4bvvvsuxzTfffENSUhIjRowo8v5eeOEF\nhBAsXbqUCRMmsHr1al544YUcbR5++GHWr1/Pu+++y/Lly3FxcWH06NF5XmvOnDk88sgj1KlTh1Wr\nVjF37lz+/vtv2rRpQ2JiYo6227Zt4+OPP2bixIn8+OOPNGigpgR3dL4P9gRnZ2789rPFsczat4eD\nB+H6dTsEh+Ndk7F09LU4lqgQYiDQGGiTbfFdUsoYIcTdwC9CiMNSyhMWtn0KeAogJCSk2BdJ9Xp9\njm3XrGnEyJEn2L491/+myYT3hlU4AceDqnLMyhdlc8dla76+vnkOXqD1fLG03N4KE9eYMWMoX748\nW7Zsybrg2qRJEwYMGABA69atmTVrFp07d87aZvbs2bRv356goCASExOz7pFKTk7O2l9mL6HM55lt\nWrRowfvvv4+zszMtWrTg8OHDfPnll3z22WcIIfjll1/YtWsXixYt4tFHtakiWrZsycMPP8z58+ez\n9qHX6xk7diwDBw7kk2wTwNetW5dGjRoxc+ZMRo0aBWjlsfj4eHbs2EFISEhW29zvTWn7fzQYDHb9\nPNji8+hZowEuR/7g98+nk966U571ERH3MnfuGZo2vTkHjc2OE1JKh3kALYDN2Z6PB8ZbaNcBOAIE\nF/Bai4FHb7XPyMhIWVzbtm3L+jkuTkofHykNhrztru/aKqMidPLvB+tKs9lc7P0VJy57+Pfffy0u\n1wYft+/Dkhs3bhT4+yQlJUknJyc5duzYfNssX75cCiHksWPHpJRS7tu3TwJyzZo1WW22bdsmgRz/\nP23atJFt2rTJ02b58uU54pozZ44E5MWLF6WUUk6ePFk6OzvLtLS0HHEsXrw4xz5++uknCcitW7fK\n9PT0HI/69evL3r17Z21bpUoV2a5duwLfCylv/X7ZS35x5ff3aCu2+Dxe/Xa+jIrQyf+e6GZx/euv\nS/nWWyUbF3BAFuK47mjlsv1AdSFEuBDCFegP5KgmCiEaAnOBHlLKK9mW+wsh3DJ+DgRaAf/aKvBt\n26BVK22I7dxUqUxz40ai3dNMccTHx2M2m6lUKf/Lf71796ZChQrMnTsX0EpUYWFhdO/evVj7LF++\nfI7nbhl/WJlDpFy8eBF/f390Ol2OdtnPQACuXNE+Ih06dECn0+V4HD58mNjY2BztQ0NDixWvYl+3\nKpm1amW/mzIdqlwmpTQKIZ4DNqN1YV4kpfxHCDEFLWuuA6YD3sDKjAN2Zlfl2sBcIYQZ7VrTNCml\nzZJMftdjpNFIwk/asP7+agbMUsnf3x8nJycuXLiQbxudTseTTz7JrFmzeO2111i2bBkvv/yy1XoO\nhoaGEh8fT3p6eo5Ec/ny5RztAgICAFi8eDF169bN8zq5R1O+k78ElWZZw//v/pmELesI7DM0x/oW\nLWDfPkhPh1zfS6zO0c5kkFJulFLWkFJGSCmnZiybkJFgkFJ2kFKGyFxdlaWUv0kp60sp78n4d6Et\n4969G+67L+/yxL3bMcXH4nZ3Ddxr1LNlSEoJ8fT0pHXr1nz99dekpKTk227kyJFcv36dPn36kJqa\nWqwL/oXVokULTCYTq1fnnB1x2bJlOZ63bNkSHx8fjh8/TuPGjfM8atasabUYFdvyz5ibKsHCjJn+\n/toQM4cO2TgoHOxMprRKSYH//oP6FqaGiVczYJYJH374IW3atKFFixa8/PLLVKpUiZMnT/Lnn3/y\n2WefAVCxYkW6d+/Od999R/fu3alcubLV4nnwwQdp3bo1I0eO5Nq1a1SvXp3ly5fz999/52hXrlw5\npk+fzqhRo7h69SqdO3fG19eXCxcusGPHDtq2bZvVeUEp3Xwf7AkTR2eVzFz8cpZcW7WCXbu0mTNt\nyeHOZEqjw4ehZk3IPS9SjlJZGZ0B807RpEkTdu/eTeXKlRk9ejRdunRh+vTpea7T9OnTB9DOaqxt\nzZo1dOnShfHjx9OvXz+MRiOff/55nnYjR45k3bp1REdHM2jQIDp37szEiRMxGo3ce++9Vo9TsY3M\nkhlGIwlb8t4YY6/rMupMpgRERUFkZN7lqlRWtjRs2JAfbnHr9Pr166lSpUqOrsyZ2rZtm9nzMUvu\nLqTZ22Tvjjt06FCGDh2ao21QUBDffvttnv3k3gdAly5d6NKlS4Gxnz59usD1iuPz7/yIdl1m0+o8\n12Vat4axY7UOMLYsqqgzmRIQFQWNGuVdrkpld469e/cyZ84cli9fzksvvYRT7hFSFcUGCuplVrWq\n9q+tv0uoT0IJOHgw75mMTE9XpbI7SIsWLXj11VcZMmQIzz77rL3DUe5QBZXMhIB77oFcl+2sTiWZ\n25SaCkePQu7RNxJ/36FKZXcQKSWJiYksXLhQDXiq2FVBvczq1IF/bXZjh0Ylmdt0+DBUqwYeHjmX\nq1KZoij2UFDJTCWZUsjSRX9VKlMUxV4KKpmpJFMKHT6s1TmzyyqVRdRUpTJFUWzOv7M2ukjuklnt\n2nDkCJjNtotFJZnbdPIkRETkXKbGKlMUxZ78Oloumfn5ga8vnDtXwMYlTCWZ23TqFNx9983nMj2d\nhC1rAVUqUxTFPlzKBzpMyUwlmdtgNmt9zjP7n4MqlSmK4hiySma5ZsysW1clmVIjLs4VX1/w8rq5\nTJXKFEVxBFklsz2/5CiZ1akD//xjuzhUkrkNFy96EB5+83mOUpka1l+xke3btyOEsOvsj4rjya9k\nFh4OZ87YLg6VZG5DTIx7jusxWWOVRdTEvXreuTsURVFsyVLJLDQULl60XQwqydyGS5fcc5zJZN2A\nqUpliqI4AEsls9BQiImxXQwqydyGmBiPrDMZVSor2yZNmoQQgsOHD9OuXTs8PT0JDQ1lwoQJmDNu\nOjAYDIwZM4Z69erh7e1NhQoV6N69O0ePHs3xWosXL0YIwd69e3n88ccpV64cYWFhPP/881nTK2dK\nTk5m7NixhIeH4+rqSnh4OFOnTs3aZ342b95My5Yt8fX1xdvbm5o1azJlypSSfVMUh6eVzNrlKJmV\nL6/NgZWaapvDv0oyt+HSJfesnmWqVHZn6NWrFx06dOD7779nwIABvP3221kH79TUVBITE3nzzTfZ\nsGEDs2fPxmAw0Lx5cy5dupTntQYNGkRERARr1qzhmWeeYebMmbz33ntZ641GI506dWLBggW88MIL\nbNq0iSeffJK3336bV199Nd8YT548SY8ePQgPD2f58uWsW7eOl156iaSkpJJ/QxSHlzWWWUbJTAio\nUAFiY11tsn81kt9tuHFDR1CQ9rMqlRXsYDXb/EHfSqPjabe1/YgRIxg3bhwAHTt25MaNG3z00Ue8\n+OKL+Pn5sWDBgqy2JpOJTp06ERISwrfffsuYMWNyvNaAAQOYPHkyAB06dOD333/n22+/zVq2cuVK\ndu3axY4dO7j//vsBeOCBBwCYPHkyY8eOJTg4OE+MBw8eJC0tjdmzZ1OuXDkA2rdvf1u/t1J6+XXs\nydmJz2kls/hYXPwDCA21XZJRZzK3Qa93wd9flcruJH379s3xvH///uj1+qxpj1esWEGzZs3w8/PD\nxcUFLy8v9Ho90dHReV6ra9euOZ7Xr1+fs2fPZj3funUrVapUoWXLlhiNxqxHx44dSU9PZ+/evRZj\nvPfee9HpdPTv359Vq1Zx5cqV2/21lVLMUsksLAzi4txss3+b7KWMSkx0wc9PlcoKo9HxNBITE/Hx\n8bF3KLclJCTE4vMLFy7www8/0K9fP4YMGcLEiRMJDAzEycmJLl265LnWAlC+fM452N3c3EhNTc16\nfu3aNc6cOYNOp7MYS2xsrMXl1apVY/Pmzbz//vsMGjSI1NRUmjRpwgcffECbNm2K9PsqZYM2Y+ZW\nbcbMvk/Y9ExGJZliSk0Fo1Hg5QVns0plalj/su7y5cvcna3f+uXLlwGoWLEis2fPplq1aixevDhr\nfXp6OnFxcblfplD8/f0JDw9nxYoVFtdXzT7URC7t2rWjXbt2pKamsnv3biZMmEDXrl05ffo0gYGB\nxYpHKb1yl8xCQwM4ckSVyxxaQgL4+BjBmG1Y/y5qrLKyLvcBf9myZXh7e1OvXj2Sk5PzTFj21Vdf\nYTKZirWvBx98kHPnzuHt7U3jxo3zPAqTLNzc3Gjfvj2vvfYaSUlJnDp1qlixKKVb7pJZcDBcv275\nDLnE922TvZRB8fHg7W0kce9uTAlxqlR2h5g/fz5ms5kmTZqwefNmFixYwKRJk/Dz8+Ohhx7i+++/\nZ8yYMXTr1o2oqCg+/fRT/Pz8irWvvn378u233/LAAw/w8ssvc88995CWlsaJEydYt24d33//PZ6e\nnnm2mzNnDjt37qRLly5UrlyZa9eu8d577xEWFka9emo8vTtV9pKZR/snSE11tsl+VZIppoQELcnE\nq1LZHWXt2rWMHj2at99+G19fX958803eeustQOt5du7cORYtWsTcuXNp0qQJP/zwA7179y7WvnQ6\nHZs3b2batGnMmzePU6dO4eXlRUREBF27dsXV1XK545577mHTpk2MHz+eK1euUL58eVq3bs0333yD\nR+4pXJU7RvaSmWebWNLSbFTIklLe0Y/IyEhZHBs3Stms8SX5Z2SIjIrQyeTow8V6HWvYtm2bXff/\n77//Wlx+48YNG0dSOIWJa+LEiRKQ6enpNohIU5rfL3vIL678/h5txd6fx+yODe4soyJ0cuvri2Tz\n5tdu67WAA7IQx1h1TaaYEhLgXpedmBLicI+opUpliqI4vMxbLLwOr7bZmYxKMsVkMECjVK3PuZ+6\nAVNRlFLA78Ee4OyM7ujP6AwJNtmnwyUZIcRDQohoIcRxIcQ4C+vdhBDLM9b/LoSomm3d+Izl0UKI\nTtaM05ycSB39ZkD1KrsTTJo0CSllnt5jilKaZPYyE2YT96RuxJSkt/o+C51khBC7hRCDhBBWu01U\nCOEMzAQ6A3WAx4QQdXI1Gw7ESymrATOA9zO2rQP0B+oCDwGzMl6vREmjkbOTX8QwbQCe5gQQgitL\n5yONxpLelaIoSomSRiPSpB2rmhjX8lfTMM5OftGqx6+inMmkA0uAGCHEx0KIWlaIpylwXEp5UkqZ\nBiwDeuZq0zMjDoBVwANCq1X1BJZJKVOllKeA4xmvV6LOTX2F2JWLCRMntAVSErdqCeemvlLSu1IU\nRSlR56a+gv6P3wFo5LIdmWogduViqx6/Cp1kpJRtgdpoB/jBwD9CiO1CiH5CiJK6q6cicC7b8/MZ\nyyy2kVIagetAQCG3vS2mJD2xKxYhDckEudyckEEakold+YVNTj0VRVGKI/P4RWpKjuXWPn4VqcAs\npYwGXhJCjAf6Ak8BS4FrQogvgHlSypO3EY+lq+eykG0Ks632AkI8hRY7ISEhhZ62VqalYhg1BWk2\nIxZ+hWfiPyR06oehdkOEkxMJu3cjXG0z6FxB9Hq9Xafi9fX1JTExMc9yk8lkcbm9qbiKprTFZTAY\n7Pp5sPfnMVPW8ctopMKnrwOQ2OJBkpo9YNXjl9C6OxdzYyEaAR8D92csMgPfAaOllHkn0Lj167UA\nJkkpO2U8Hw8gpXwvW5vNGW32CCFcgEtAEDAue9vs7QraZ+PGjeWBAwcKFZ8pSc9fTcOQqQb+SmlK\nA499N2N396DB7xdw9vIu/C9sJdu3b6dt27Z22/+RI0eoXbt2nuWOOkCmiqtoSltc+f092oq9P4+Z\nsh+/AC6Z76KCkzbqd3GOX0KIKCll41u1K3LvMiGEhxBimBBiH7Af7QD/AhAGPAO0BL4p6utm2A9U\nF0KECyFc0S7kr8vVZh0wJOPnR4FfMm4MWgf0z+h9Fg5UB/ZRgpy9vAnoOwzh7slF412YM0+e3DwI\n6POEQyQYRVEUSzKPXzhrBaw9pu4ACHdPqx6/itK7rL4Q4nMgBpgDnAE6SCnrSik/k1JeklLOB54G\nWhUnmIxrLM8Bm4EjwAop5T9CiClCiB4ZzRYCAUKI48BL3DyD+QdYAfwL/AiMklIWb2TCAlR+40MC\n+gwl1cWP07qGAHg3bknlNz4s6V0piqKUqErj3kdkTB1xwKkzwt2DgD5DrXr8Kso1mUNoCeZ/aNde\nLubT7jhQYImqIFLKjcDGXMsmZPvZAPTJZ9upwNTi7rswhIsLd038H1VCDRxe8DF3xx3E2d0Doe6f\nUG7T0KFD2b59O6dPn7bpfhcvXswTTzzBqVOnCpw+IL9tzWYzw4YNs05wSolK+vN3pCEFU0AVLqZE\n0mCP9Uv8RSmX9QGqSCknF5BgkFIekVK2u/3QHJt/sDv7nHqCENzY+ROmxBv2DklRiqVr167s2bOH\n0NDQIm+7ePFiFi1aZIWoFGtIyBjQ19CwD67uwiYl/qJ0YV5tjfJTaeXnBzEpFfFu3BqZnkbCzz/Y\nOyRFKZagoCCaN2+Om5v9e0Yq1iNNJuJ/XAPAjbqP4upqtsl+HW5YmdLC3x/0epesAecyvyEoZdeh\nQ4fo0aMH/v7+eHh40KpVK3799des9UOHDqVSpUr88ccf3HfffXh6elK9enXmzJmT57V+/vlnGjVq\nhLu7OxEREcydO7fQcSxevBghBDt37qRXr154e3sTEBDAqFGjSEnJeQ/ExYsXGTx4MIGBgbi5udGg\nQQO+/vpri6+XvUxXtWpVBg4cyLJly6hduzZeXl40btyYXbt2ZbVp27YtO3bsYPfu3QghEEJk9aK6\ndOkSQ4YMISwsDDc3N0JDQ+nWrRtXrlwp9O+plCz9gd0Yr13GtfLd3AhoiJubbc4ZVJIpJn9/SEx0\nwe+h3tlKZtftHVaZZUrSYzh1zG43vB48eJCWLVsSFxfH/PnzWb16NQEBAXTo0IGoqKisdjdu3GDA\ngAEMHDiQtWvX0qRJE5555hm2bduW1ebIkSN06dIFDw8Pli1bxrvvvsv//vc/fv755yLFNHDgQKpV\nq8aaNWsYM2YM8+fP55lnnslan5SURJs2bdi0aRPvvvsu33//PfXr12fQoEHMmzfvlq//66+/8tFH\nH/H222+zfPlyTCYT3bp1IyFBG1hx1qxZNGzYkAYNGrBnzx727NnDrFmzABg0aBB79uxh+vTpbNmy\nhU8//ZRKlSqRnJxcpN9RKTnxG1cB4N/5YVIMwmZnMnafz8Xej+LOJ3P9upTu7kYppZTR/dvLqAid\nvPbd18V6rZJm7/krSnI+GXN6ujwz6QV5sI6P/KOBvzxYx0eemfSCNJfgvC6Fiat9+/ayVq1aMjU1\nNWuZ0WiUtWrVkj179pRSSjlkyBAJyF9++SWrjcFgkAEBAXLEiBFZywYMGCADAgKkXq/PWnb27Fmp\n0+lklSpVbhnXF198IQE5cuTIHMvfeecd6eTkJKOjo6WUUn722WcSyPP38MADD8igoCBpNBpzvN6p\nU6ey2lSpUkX6+fnJuLi4rGX79++XgFywYEHWsjZt2shWrVrlidHLy0t+8sknFuO3FjWfTP7MRqM8\n1KySjIrQyaTDUfLjj6V8+OFzt/WaqPlkrMvHB9LSnEhPR5XMrChzrDiZasCcpLfJWEu5paSksGPH\nDvr06YOTkxNGoxGj0YiUkg4dOrBz586stp6enrRrd7Pfi5ubG9WrV+fs2bNZy/bs2UOXLl3w8vLK\nWla5cmVatcrZ899kMmXtK3N/2fXt2zfH8/79+2M2m9m3T7s9bOfOnVSsWDHPjYADBw7k6tWr/Pvv\nvwX+3i1atMDf3z/ref369QE4f/58gdsBNGnShOnTp/PJJ59w+PDhPLErtpW9VOZRtyEXL0JAQJpN\n9q2STDEJoU2/nJCAKplZSfax4rKz9VhxcXFxmEwm3n77bXQ6XY7H559/Tnx8PGazVnrIflDO5Obm\nhsFgyHp+8eJFQkJC8rTLveyee+7Jsa8lS5YU2D7z+YULF7LittRjrEKFClnrC1K+fPk8vweQ43fJ\nz/Lly+nRowcffPABDRo0oGLFikyZMiXrfVJsK3upTAjBxYtQvnyqTfatbu64DT4+6cTH6wiqUQHv\nJveh37eThJ/XE9DrcXuHViakX4lBuLggLXwWhLMz6VdicA6vYfU4/Pz8cHJyYtSoUQwePNhiGyen\nwn9fCw0N5fLly3mW5162fPnyHPPXhIeH52lft27dHM8BKlbUxoUtX7480dHRefZz6ZI24lNAQECh\nYy6q4OBgZs6cycyZMyZCoaUAACAASURBVImOjmbJkiVMnDiRoKCgHNeNFOuTJhMJm78DblZdLl6E\nBg1scyajksxt8PY2Eh+v/ezf+REtyWxarZJMCdEFh+U7z4U0mdAFh9kkDi8vL+677z4OHTpEo0aN\nipRQLGnRogUbN24kKSkpq2R27tw5du/eTVjYzd+pbt26BY4RtmLFCtq3b5/1fNmyZTg5OdG0qTbD\nRZs2bVi5ciW7d+/OUYpbunQpwcHBJTKel5ub2y0Hy6xZsybvvvsuc+bM4e+//77tfSpFo9+/K0ep\nDCAmxnblMpVkbkNQUCrnz0OzZlrJ7NyUF7NKZs4+vvYOr9TLHGspduXiHCUzbayloTYdK+7jjz/m\n/vvvp1OnTgwfPpzQ0FCuXbvGwYMHMZlMTJs2rdCv9eabb7Jy5Uo6duzIq6++SlpaGhMnTrRYQivI\nxo0befXVV+nYsSP79u1j8uTJDB48mBo1tLO7oUOH8sknn/Dwww8zdepUKlWqxDfffMOWLVuYO3cu\nzs63P6dfnTp1mDVrFsuXLyciIgIfHx8qVKhAhw4dePzxx6lVqxY6nY61a9cSHx9Px44db3ufStHE\nZ1wrziyVATa9JqOSzG2oUMHAyYyJDXRBqmRmDZljKsWu/ALh7Iw0maw+1pIljRo1Yv/+/UyePJnn\nn3+e69evExQURKNGjXj66aeL9Fq1a9fOShD9+vWjYsWKjB07lj179hRpSPivv/6ajz76iNmzZ+Pq\n6sqIESP48MOb74uXlxc7duzgtddeY9y4cSQmJlKzZk2++uorBg4cWKSY8zN27Fiio6N58skn0ev1\ntGnThs2bN9OoUSPmz5/PmTNncHJyombNmnzzzTf07Jl7DkLFmnKWyrTRuFJSIDkZypVLt1EQDtCN\n2J6P4nZhllLK558/Jp955ubzK1/OklEROnl8RK9iv2ZJsHeXyZLswpzJqE+UKSejpVGfWOzXyM/t\nxGVNt+rC/N9//9k4Ik1pe7/u5C7MN/Zsl1EROnm4bU1pNpullFKePCll5cq3HxeqC7P1VaiQknUm\nA9l6mf26RfUyK2HOXt64h9dQ0ykoShHkVyoLs83lTEB1Yb4tYWGGHEkms2SmjWW23n6BKYpyx7NU\nKgM4f14lmVKjQgUDZ8+CKdsQQP6dHwEgIaNfuqKUtKFDhyKlpFq1/7d33/FNVe8Dxz+n6W6BlkIL\nBWRDUVFkOpAlQwFliYKAKFVxgIogiAqy/AIqVkVRfwwZoqgoioAiIKAiG2WWDTLKpkBLZ9rz++Om\nJWnTSTNanvfrlVeTe89NntwmeZLnnntOLVeHItyYba+yBpnL9+6FiAjnxSFJ5jr4+KRTtqzRHTCD\nlMyEEO7AXqkMYM8euPlm58UhSeY6Va8OR45cu21TMlspw/8LIZwvp1IZSJIpdmrUwOa4DFiVzGQs\nMyGEC+RUKjOb4cABKZcVK3XrQnS07TIpmQkhXCmzVNaxh02p7PBhqFgR/P2dF4skmevUsCFYTScC\nSMlMCOE6NqWyBx62Wbd7N1gNd+cUkmSuU8OGsG0bZB3JXIb/F0K4Qk6lMnD+8RiQJHPdKlQAPz+w\nmrkWgKAOXaVkJoRwusxh/bOUykCSTLHVqJHxa8aalMzE9Zg9ezbz5s1zdRiFZjabUUoxYcKEAm+7\nbds2xowZkznNs8g/nZbGpd9+BLKXysAolxXB4NsFIkmmCDRqlP24DEjJTBRecU8ynp6erF+/nief\nfLLA227bto2xY8dKkimE3Epl8fFGz7LbbnNuTJJkioC9g/8gJTPhXpKTnTMTYoY777wzcwI14Ry5\nlco2boQGDcDX17kxSZIpAhnlsqwH/6VkVrJs376dbt26ERISgp+fH3Xr1mXixImZ63/44QfuvPNO\n/P39CQoKomfPnhw7dszmPqpVq0bfvn1ZsGAB9erVIyAggMaNG/PXX39ltmnVqhVr165lw4YNKKVQ\nStGqVavM9UeOHKFPnz6UL18eHx8fGjRowKJFi2weZ8yYMSil2LVrFx06dCAwMJBHHnkkx+eWUd56\n6623GD9+PJUqVcLX15eWLVuyc+dOm7Zaa6ZMmUKdOnXw8fEhPDycwYMHEx8fn+3+rMtlb775Jkop\nDh06xAMPPEBAQADVqlVjwoQJmdMyz5gxg6effhowZgLNeP4nTpwAjHl96tWrh5+fH2XLlqVJkyYs\nXrw41//bjSKvUtm6dWA1d53TSJIpAuHh4OkJx49nXycls6IRHx/P/v37bT7InGnTpk3cddddHDp0\niKioKJYuXcorr7yS+eH32Wef0aNHD26++WYWLlzI559/zq5du2jZsmW2mSP//PNPpkyZwvjx4/nm\nm29IS0ujc+fOmeWhadOmcccdd3Drrbeyfv161q9fz7Rp0wBjBs1mzZqxfft2oqKiWLx4MQ0bNqRH\njx52P2y7dOlCy5YtWbx4MUOGDMnzec6aNYvffvuNTz75hC+++IKYmBjatGljU7oaPXo0w4YN4/77\n7+fnn39m2LBhzJo1i86dO2cmi9x069aNdu3a8dNPP9G5c2dGjRrF/PnzM+MdOXIkYCTtjOcfGhrK\nnDlzGDFiBH369GHZsmV8+eWXdO/enQsXLuT5mDeCjFKZz001s5XKwEgyzZu7ILD8zAdQki/XM5+M\n9XwMDz6o9TffZG+TcvaU3lrLW2+LCNDmK5cK/ViFjcsVinI+mdTUVD1o0CDt5+enAwMDtZ+fnx40\naJBOTU293jALFNe9996rK1eurK9evZptXVxcnC5durR+8sknbZYfOXJEe3l56aioqMxlVatW1UFB\nQfrixYuZyzZv3qwBPX/+/MxlLVu21HfeeWe2xxowYIAuV66cPn/+vM3ytm3b6ttvvz3z9ltvvaUB\n/cEHH+T53LQ29jOgy5cvb/McDx48qE0mkx4zZozWWuuzZ89qLy8vHRkZabN9xhw3S5cutbm/8ePH\nZ7Z54403NKDnzp1rs21ERIR+4IEHMm9Pnz5dA/rIkSM27QYOHKibNGmS43O40eeT+W/UIL21ppc+\n8c7r2daZzVqXLq31uXNFFxfFbT4ZpVRZpdQKpdQBy99gO20aKKXWK6V2K6V2KKUetVo3Wyl1RCn1\nr+WSPZU7UJs2sGpV9uVe5SsQ2LSFlMwKaciQIcyaNYvExETi4+NJTExk1qxZ+fpWXlQSEhJYt24d\nffr0wd/OqdLr16/nypUr9OnTB7PZnHmpXLkyERER/PHHHzbt77rrLoKDr72869evD5CttGbPr7/+\nSseOHSlTpozNY3Xo0IHt27dz5coVm/bdunWzuZ2enm6zXdZfHp07d7Z5jjVr1qRJkyasX78+87mm\npqZmm1mzd+/eeHh4sHbt2jyfQ6dOnWxu33rrrfl67k2aNGHr1q289NJLrFq1ioSEhDy3uVHkdgIm\nwM6dxpn+5co5OzL3Kpe9BqzSWtcGVlluZ5UAPK61vgW4H/hAKRVktf5VrXUDy+Vfx4d8Tdu2sHKl\n/XUZY5nFyvD/BRIfH8/MmTOzfZgkJCQwc+ZMp5XOYmNjSU9Pp3LlynbXnz17FoC2bdvi5eVlc9m5\nc2e2ck7ZsmVtbvv4+ACQlJSUZyxnz55l7ty52R7n1VdfBcj2WBUrVrS5PXr0aJvt2rdvb7M+LCws\n22OGhYVx8uRJAC5evGj3fn18fAgODs5cnxt7zz8/z33AgAF8/PHH/P3337Rr146QkBB69OiRrwRV\n0sVv+hPzhbPuVyoDPF3zsHZ1AVpZrs8B1gAjrBtorfdbXY9RSp0FygMu7+t4yy1w9aoxNlCNGrbr\ngjp05fjYl4j7ayVpcZcxlSrjmiCLmZiYGEwmk911JpOJmJgY6tSp4/A4goOD8fDwyPygzSokJAQw\nuh3fYmfMjlKlShVZLCEhIdx7772MGDHC7vrwLLNRZe1h9Pzzz9O1a9fM26VLl7ZZf+bMmWz3eebM\nmcxeYhkJ4vTp09StWzezTUpKCrGxsZn7whGUUjz33HM899xzXLx4keXLlzN06FB69+7NunXrHPa4\nxUHGWGVBWYb1z7BuHbRr5+yoDO6UZMK01qcAtNanlFKhuTVWSjUFvIFDVovfVkqNxvJLSGvttD6b\nShm/Zlatyp5kMkpm8RvXcmnlz4R062v/ToSN8PBw0qxnhLOSlpaW7QPVUfz9/WnevDlffvklo0eP\nxs/Pz2b93XffTalSpTh48CD9+/cvksf08fHh8uXs3d7vv/9+1q9fzy233JItjvwIDw/Pdb8tWbKE\nhISEzJLZoUOH2Lx5M6NGjQKMUp+XlxcLFiygZcuWmdt9/fXXpKen2ywrrIxfdomJiTm2KVu2LL17\n92b9+vXMmTPnuh+zOMurVKY1/PUXjB3r7MgMTk0ySqmVQAU7q94o4P1UBOYB/bXWGUXlkcBpjMTz\nfxi/gsblsP0zwDNglALWrFlTkIfPFB8fb7Nt5cphfPVVCLVr78nW1qtWffw2ruXgvM/ZGWy/7FJU\nssblbGXKlMnWowqMxGBveW769evHvHnzbD5w/Pz86NevH1rrAt+fPfmJa+zYsXTs2JFmzZoxaNAg\nKlWqxNGjR9mxYwfvvfce48ePZ+jQocTExNCuXTtKly5NTEwM69ato3nz5pndh7XWpKam2n285OTk\nzOW1atVizZo1zJ49m+rVq1OqVClq167N8OHDad26Nffccw/PPPMMVatW5dKlS+zZs4ejR49m9kLL\nOCcmLi4OT8+83+ZmsxkAb29v2rZty4svvkhiYiITJkwgKCiIAQMGEBcXh6+vL88++ywff/wxJpOJ\ntm3bEh0dzdtvv03z5s25++67iYuLy7w/6+dkHZO11NRU0tPTM5dXrVoVMLorP/roo3h5eVG/fn1e\nfvllgoKCaNq0KeXKlePAgQPMnz+fNm3aEBcXl+P/MSkpyaXvB0e/H017/yXgwlnSy4ez6dwlyPJY\n//3nT3LybZw4sQHrH+NO+5zIT+8AZ1yAfUBFy/WKwL4c2pUGtgE9c7mvVsCS/DxuUfUu01rr48e1\nDgnROi0te9uUs6f01to+TullJr3LCia/cW3btk137txZlylTRvv6+uq6devqSZMmZa5funSpbtWq\nlS5VqpT29fXVNWvW1E8++aTevXt3ZpuqVavqPn36ZLtvQL/11luZt0+dOqXbtWunAwMDNaBbtmyZ\nue748eM6MjJSh4eHay8vL12hQgXdtm1bPW/evMw2Gb3L8rufMnqDjR49Wo8bN06Hh4drHx8f3aJF\nC719+3abtpcvX9bvvvuurl27tvby8tIVK1bUgwYN0nFxcdnuz17vsqz69Omja9asabNs1KhRumLF\nitrDw0MD+vjx43rWrFm6RYsWuly5ctrHx0dXr15dv/LKK5n/vxu1d1luvcq01nryZK2fe67o4yKf\nvctcnlwyA4F3MUpcYBz0f8dOG2+MUtjLdtZlJCgFfABMys/jFmWS0VrriAitt261337fY2311ppe\n+uj86Xrfvn02b8qiVJKSTIa4uDiH7bPricuRnBlXRlKwTnQ5KW77qyQnmSuXLumf6ofqP6t76qu7\nttlt07y51r/8UvRx5TfJuFPvsklAO6XUAaCd5TZKqcZKqRmWNo8ALYAn7HRVnq+U2gnsBMoBBR+Z\nrwjk1susVPuuvHPOTN3HB9KoUSNCQ0MZPHhwZmlB5CwwMJA6deoQGBjo6lCEcDmz2czgwYMJDQul\n966ztP0vjRGfzcr2WXL+vNF92WrACKdzmwP/WusLwH12lm8BnrJc/xL4Moft2zg0wHxq1w6iomD4\n8OzrJvy5jZ/iNMlak2zpfjtr1iwApk6d6swwhRDFWMb5Y0nJKcYCre1+lixbZpzD5+zxyqy50y+Z\nEqFdO/jnH8jaEzQ+Pp4vvvqapCzjmzn7nA8h7PH09ERrzZgxY1wdishDQc4f+/lnePBBZ0doS5JM\nEfPzg06d4PssQ5Xl55wPIYTIS34/S1JSYMUK4/PIlSTJOMCjj8K339ouc5dzPoQQxVt+P0vWrjUm\nKAvN9YxDx5Mk4wDt28P27WD94yQwMJDIyMhsY1/5+/sTGRlZ4g5oG51PhHCtkvg6DAwMJHLAAHw9\nbM/sz/pZ8vPP8NBDrojQliQZB/D1NeqgWUtmUVFRDBgwAF9vL/wV+Jo8GDBgAFFRUa4J1EG8vLxy\nPVtbCGdJTEzEy8vL1WEUuXGPdqNLIPh4KAIDA/Hz87P5LDGbjWpKjx4uDhRJMg5jr2Tm6enJ1KlT\nidm3ly9v8mJlTR+i3h6frzOyi5PQ0FBOnjxJQkJCifwmKdyf1pqEhAROnjxJqKvrRQ4Qt2IRw8t7\nsmPcMLZu3crZs2eZOnVq5mfJ8uXG8FZOGNovTyXr082NtGsH/frBiROQdfDe4Go1uOWelsRvXMvl\nVT8T0q2fa4J0kIxBF2NiYkhNTc1cnpSUhK8r+1LmQOIqmOISl5eXF2FhYdkGAS3ujLHKjBkwK3ft\njb+dTDJ7NjzxhHPjyokkGQfx9oYuXWDhQnj55ezrgzs+TPzGtcQu+77EJRkwEk3WN/eaNWu44447\nXBRRziSugpG4XMtmWP+bsw/rf+EC/PYbTJ/uguDskHKZA/XqBV/aPXXUGP4fDw/i/lqB+YrLZyoQ\nQhQTsb8Y81IFdexhd1j/BQugY0cICsq2yiUkyThQ27bGt4otW7Kv8yoXZpkxM5XLq2TGTCFE3qxL\nZRmTIWY1Z477lMpAkoxDmUwwcCB8+qn99ddmzPzefgMhhLCSV6ls9244edL4gusuJMk42IABRlfm\n2Njs66RkJoQoiLxKZXPmGB2OchgQwCUkyThYaKhRH7U3eZ+UzIQQ+aXN5lxLZampxjHgIpqctchI\nknGC55+Hzz4zpkHNSkpmQoj8iN9sKZVVrWW3VPbdd8Z5MfXquSC4XEiScYJ77gEvL1i9Ovs6KZkJ\nIfIj9hfji2jQA92zlcq0hilTYOhQV0SWO0kyTqAUPPec/Q4AUjITQuQlr1LZ2rVw9arrR1y2R5KM\nk/Tta8yYaW9EfymZCSFyk1epbMoUGDIEPNzwE90NQyqZSpc2Es1HH2VfJyUzIURuciuV7d0LmzbB\n44+7IrK8SZJxoldfNYZ6uHDBdrlNyWzlYtcEJ4RwS9ps5tKviwD7pbKoKHj2WWPCRHckScaJbroJ\nuneHDz/Mvi6zZPbLD06OSgjhzuI3/4n54jm7pbJz54zR3l94wUXB5YMkGScbORKmTYNLWapiUjIT\nQtiTW6ls2jTo2dP1s1/mRpKMk9WoYfQAmTrVdrmUzIQQWdmUyjo+bLMuLs5IMkOGuCKy/JMk4wKv\nv250AIiLs12e8SKSkpkQArKUyurdbrPu/feNqd7d7eTLrCTJuEDdusakZtOm2S6XkpkQwlpOpbJz\n54wvqmPHuiqy/JMk4yJvvGH0Crl69doyr5BQKZkJIYDcS2UTJ0Lv3kb53d1JknGRW26B5s2z/5qR\nkpkQAnIulR07Zgy4++abLgyuACTJuNCECfDOO8ZP3wxSMhNCAMQuswzrn6VUNnasMU9VhQquiqxg\nJMm4UEQE9OkDo0dfW+YVEkqpZi2lZCbEDcxmrDKrUtnevbB4MQwf7qrICs5tkoxSqqxSaoVS6oDl\nb3AO7dKUUv9aLoutlldXSm20bP+NUsrbedEX3ujR8MMPsHPntWVBcmKmEDe0nEplb74Jw4ZBUJAL\ngysgt0kywGvAKq11bWCV5bY9iVrrBpbLQ1bLJwNRlu1jgUjHhls0ypY1Es2QIdfmm5GSmRA3tsxS\nmdUMmBs3wvr1MHiwKyMrOHdKMl2AjPkj5wBd87uhMv4LbYCFhdne1QYOhFOnYMkS47aUzIS4cdkb\n1j8tzZj8cOJE8Pd3ZXQF505JJkxrfQrA8jengRJ8lVJblFIblFIZiSQEuKS1NltunwAqOTbcouPp\naZxYNXQopKQYy4Jk+H8hbkhxm/7IViqbNg1KlYJ+/VwcXCEobW9OYEc9mFIrAXt9It4A5mitg6za\nxmqtsx2XUUqFa61jlFI1gN+B+4ArwHqtdS1LmyrAMq11/RzieAZ4BiAsLKzRggULCvV84uPjCQwM\nLNS29rz2Wn0aNozlkUdOoK7EEji0t1E2i/oW/PP/OEUdV1GRuApG4iqYkhKX77wP8V67lOSOvUnu\n/iQXLngTGdmYDz74l2rVElwWV1atW7feqrVunGdDrbVbXIB9QEXL9YrAvnxsMxt4GFDAecDTsvwu\nYHl+HrdRo0a6sFavXl3obe2Jjta6XDmtT582bu/v215vremlz38/x6VxFRWJq2AkroIpCXGlp6bq\n7U3C9daaXvrq7n+01lr36qX1yJGujcseYIvOx2esO5XLFgP9Ldf7Az9lbaCUClZK+ViulwPuAfZY\nnvBqjIST4/buLiICIiNh0CDjtpTMhLixZC2VrVgBGzYUnxMv7XGnJDMJaKeUOgC0s9xGKdVYKTXD\n0qYesEUptR0jqUzSWu+xrBsBvKKUOohxjGamU6MvImPGGN2ZFy606mW2bqX0MhPiBnApY6yyjj1I\nTlY8/7wxYntxO9hvzdPVAWTQWl/AOL6SdfkW4CnL9b8Bu8dZtNaHgaaOjNEZfH1h1izo0QNa7TR6\nmcWtX83llYsJ6e6m86sKIa5b1l5lkydD/frQubOLA7tO7vRLRljcfbcx+N1LL0nJTIgbRWaprFot\njptuZ+pU+7PoFjeSZNzUhAmwaRNs1FYls8uxrg5LCOEgGaWy0h160K+fYuxYqFLFxUEVAUkybsrf\nH2bOhIHDQ/FtmHFi5s+uDksI4QDWw/ovONqD8uWNky9LAkkybqxFC+jeHZbFZgz/LyUzIUqiuE1/\nYI49T3pYLd7/4XZmzQKrgZeLNUkybm7iRJi3vwtaeRD310qu7thM2tV4V4clhCgiaVfjufCN0Rl2\n0akeTJ+uCAtzcVBFSJKMmwsMhE9ml+XfpHvQ5lT292rDjqbhHBv7MtpszvsOhBBuSZvNHBv7Mtub\nVCR26XcAJJapRaf7S9b7WpJMMXDTmmEkaaOjfFpyCjo5iQvfzeb428NcHJkQorCOvz2MC9/NhpRk\nAGJSq/Jw0ogS976WJOPm0q7Gc+HbWTTzWU669kBjFGp1UgIXvvtCSmdCFEMZ72udlEBKujH1lbdK\nwjf1Yol7X0uScXOpZ2NQnp54KEjVXphU+rWVHh6kno1xXXBCiEIx3rcarcHbI4WE9ADKeZ4BQJlM\nJep9LUnGzXmFhmcee/HxSOZ4Sg2uphsjp+qEqyQd3u/K8IQQBWS+comY999CJyejFOxIaoavupq5\nXqel4RUa7sIIi5YkGTdnCggk5JEBKF/jmEwV78OsjO/O/uTbADg8sDsnJo4gPTnZlWEKIfIhbtOf\n7O3cmEu/fE+a8mbaxTFU99qDh6W7svL1J6Tnk5gC3G/KgsKSJFMMVHnjPUJ6PoHy9cMjIJAu5b9j\nbsDnbK01Gkwmzs6MYt/D95C4f7erQxVC2JGeksLJd9/gQJ+2pMQcI716I56J20qfPlA6wIxHQCDK\n14+Qnk9Q5Y33XB1ukXKbATJFzpSnJze99QGVhk0g9WwMXqHhfJUaSOPGjZjQvx31Vz9BYvQO9na7\ni0ojJkKVm10dshDCwuP0cfY/0oKEXdvAw4OAviNoP30U02Z60brT66RdfTHzfV2SfsFkkCRTjJgC\nAjFVrwNAELBkCbRs2YyvZm2i1p9DubBwNifGDcH/1sak3voDXuXtTUIqhHAGrTXnF8wgYPwrJKQk\n412pKhX/9wWdhjZn4PPQqZPRzvp9XRJJuawYi4iAb76B3k+W4urj/0f1jxdgKhOM564tRHdsyKVV\nS1wdohA3pNQL5zj8bA+Oj3oBlZJM2S6PEfHzFl6d2ZwqVWDkSFdH6DySZIq5Vq3gvfeMb0Upd3Sn\n3tJtmOvdgTn2PIcHdufYqBdIS7ia5/0IIYrG5bW/Et2pIZdXLcFUqgwJT4+k2pTZTJtdhi1b4Isv\nSs64ZPkhSaYEePxx6N8fHnwQzKUrkTBkIpVGvoPy8ub819PZ26WpUQ8WQjhMelIix8cN4VDkQ5jP\nnyGwyb1ELNmCuVlrvvsOJk+Gn34yhoq6kUiSKSHeegvq1oW+fSEdD8IiX6buD+vwrX0zyUcOsPfh\n5pz+7B10WpqrQxWixEmI3s7ebndxbu4n4OlJ+LAJ1P7yN3wqVWXLlmBeeAGWLYPq1V0dqfNJkikh\nlIIZM+DCBfj885oA+Ne7nYhF6ynffxCYzcS89yYH+rYnJeaYi6MVomTQ6emcmRHFvu73kHRgDz7V\naxOx8C8qPDscZTKxaRNMmFCP77+H2293dbSuIUmmBPHxgUWLYP36EN5/31jm4etHlVHvU3Pmz3iW\nCyN+859Ed2rExZ8XuDZYIYq5lNMnOfhER05OGoFOTaFc76eJ+GkT/rc2BCA6Gh56CIYP38e997o4\nWBeSJFPClC0L7767nY8/ho8/vra8TMsO1Fu6jTL3dSYt7jJHhzzOkVf6kxZ32XXBClFMxf76A9Gd\nGhL39+94Bpejxmffc9P4TzD5BwBw/Djcf79xHObuuy+4OFrXkiRTAoWFJfP77/Duu/D559eWe4WU\nN94ME6bh4edP7OKvie7cmPjNf7kuWCGKkbT4OP577WmODOpF2uVYSrfoQL1l2whq+2Bmm/PnoX17\nePFFo0POjU6STAlVrRr8/ju8/TbMmnVtuVKKcr2eIuKnjfjf2pCUk/+xv09bTk4ZhU5NdVm8Qri7\nq/9sZO9DTbmwcA7K24fKoz+g5szFNic9X7kCHTtCly4wdKgLg3UjkmRKsJo1YeVKGDUK5s61Xedb\noy51vv2DsGeHg9ac+XQy+x5tSdLRA64JVgg3pc1mTk2dwL5erUg+dgi/iPpE/LiB0MefR1md8HLx\nIrRtC40bG9OmC4MkmRKuTh0j0bz2Gnz9te06D29vKg2bQO35K/EOv4mEHVvY+2ATzi+YidbaNQEL\n4UaSjx1m/2P3mGv8XQAAHI9JREFUcerDcZCWRmjkEOp+/zd+dW6xaXf2LLRuDffeC598cmOdbJkX\nSTI3gHr14Lff4JVX4Lvvsq8v1dQ4aSz4oV6kJyZw7M3nOPx8T8wXzzs/WCHcgNaaCz/MJfrBJlzd\nth6vsHBqzfmFyiMn4+HjY9P25Elo2dIokb33niSYrCTJ3CBuvRV+/RUGD4bZs7Ov9ywdRPX351Jt\nymw8AktzecVi9nRqyJU/Vzg9ViFcyXw5liMvPsZ/w58i/WocQR26UW/JVkrfc1+2tkePQosWxgH+\nceMkwdgjSeYGcvvtsGYNjBkDkyaBvYpY2S6PUW/JFgIa34P53GkOPtmJ4xOGkp6c5OxwhXC6uPVr\niO7UiEu/fI9HQCBVJ02n+scL8AwOydb2wAHjF8xLLxnlaGGf2yQZpVRZpdQKpdQBy99gO21aK6X+\ntbokKaW6WtbNVkodsVrXwPnPwv1FRMDff8NXXxlvjvT07G18KlejzvyVhL8yDjw9OTd7Knu73UXi\nvp3OD1gIJ0hPTubE5Nc48HgHUk+fIKBBM+ot3kzIw/1tDu5n2L3bGJx21Cijq7LImdskGeA1YJXW\nujawynLbhtZ6tda6gda6AdAGSAB+s2ryasZ6rfW/Tom6GAoPhz/+gB07oHdvsDdzszKZqPD8a9T9\n9g98qtUiaf9u9na9i7NffIS2l5mEKKYSD0az7+HmnJ3+PihFhcFvUmfBanyq1rTb/q+/4L774J13\n4KmnnBxsMeROSaYLMMdyfQ7QNY/2DwO/aK0THBpVCRUUZByjSU83zky+nMOJ/wG3NSZi8WbK9XoK\nnZrCibeHcfDJzqSciXFuwEIUMa015778jL1dmpEYvR3vKtWps2A14S+NRnnan89x3jzo3t04rtmn\nj3PjLa7cKcmEaa1PAVj+hubRvheQpVMubyuldiilopRSPvY2Etf4+sKCBXDLLUZt+dQp++1M/gHc\nNGEaNT5diCk4hLh1K4nu3IhLv/0IQHx8PPv37yc+Pt6J0QtRMNav09TzZzj0TDeOj3kRnZxE2e79\nqLd4M4EN77K7bXo6vPkmjB4Nq1cbX8xE/ihnng+hlFoJ2JsT+A1gjtY6yKptrNY623EZy7qKwA4g\nXGudarXsNOAN/B9wSGs9LoftnwGeAQgLC2u0YEHhBouMj48n0A0nhyhoXFrD/Pk3sXRpRSZM2EXN\nmjlPcqYuXcDviyl47t6CWWve9a7M4sMxmEwm0tLS6NSpEy+88AImk+m643IWiatgiltcaWlpfPLJ\nJyxdutR4nZpT6VrGk6GlzZgCSpH4+MuYG7fI8X6TkjyYNCmC8+d9GD9+F8HBBRsZo7jtr/xq3br1\nVq114zwbaq3d4gLsAyparlcE9uXS9iXg/3JZ3wpYkp/HbdSokS6s1atXF3pbRypsXF99pXW5csbf\n3KSnpekzs6fqR4O9tK9Cw7WLv7+/HjRoUJHG5WgSV8EUt7gGDRqk/f39bV6nvgrdt04VnXzyWK73\nGROjdZMmWvftq3ViYtHG5WrXGxewRefjM9adymWLgYzh5PoDP+XStjdZSmWWXzIooytIV2CXA2Is\n0Xr3NkYHePNN48RNs9l+O+XhgX+PJ/gpwYOkLD+EExISmDlzppTOhFuIj49n5syZJCTYHrpN0vD9\n8fOklLZbLAHg33+hWTNjuP65c43ysig4d0oyk4B2SqkDQDvLbZRSjZVSMzIaKaWqAVWAtVm2n6+U\n2gnsBMoBE5wQc4lz++2webMxF0a7dsZwGfbExMTg6eVld53Jw4OYGOkYIFzv5IkTeKTbnw3WZDLl\n+Dr97jtjJOUpU4wvXXKSZeHZ70LhAlrrC0C2U2q11luAp6xuHwUq2WnXxpHx3UjKloUlS4wpnRs3\nhoULoWlT2zbh4eGk5TCVc2rCVXy2/omuXdvuOQZCOENKzHESRj+HOSXF7vq0tDTCw8NtliUmGr/i\nV6yAX36BRo2cEWnJ5k6/ZIQbMZlgwgT46CPo3BlmzrRdHxgYSGRkJP7+/jbLfT1NdC2luDDqOY4O\n6Yf5yiUnRi2E4eLSb4nu3Ai9bR3dygfgl2W8MX9/fyIjI20OfO/dC3feCbGxsG2bJJiiIklG5Kpr\nV+PEzffegwEDIC7u2rqoqCgGDBiAn58fgYGB+Pn58dSzz/HBZ/+Hh38AsUu+JbpTI+I2/uG6JyBu\nKGlxVzj66gCOvtSXtCuXKN26I7N27Cfy6adtXqcDBgwgKioqc7t584wRlAcNMkYrL13ahU+ihHGb\ncplwXxERsGmTUUZo0ADmzIHmzcHT05OpU6cyceJEYmJiCA8Pz/xmGHRXS44OfYKE7Zs50LcdYc8M\ngztaufaJiBLNdHA30WMHknL8CMrXj8oj36HcY8+glMrxdRofbySWjRuNSf7q13fxkyiB5JeMyJdS\npWD6dIiKgkcegZEjIaPUHRgYSJ06dWxKD77ValN3wRoqDHodlOLM5+8S8L+XSDq010XPQJRUOjWV\nmA/G4D95KCnHj+B3cwMiftxA+T4DbY4JZn2d7tgBTZoYB/W3bJEE4yiSZESBPPSQ0bUzOtroDLAr\nl47iysuL8JfHUOerVXhXrobp2EGiuzTj3Ff/J5OiiSKRdPQg+3q15vTH/wM0Yc8Mo+7Cv/CrVS/H\nbcxmmDzZGH/s9dfhiy8gIMB5Md9oJMmIAgsNhUWLjNFnW7eG99+3P5pzhsDG91Dv582k3NUWnZTI\n8dGDODywO6kXzjkvaFGiaK05/91s9j7UhITtm/CqUJmEoe9Qafj/8PD2znG7f/4xvhytWmV01e/X\nz4lB36AkyYhCUcroCLBpk5Fw7rsPjhzJub2pVBmSIodT7cMvMZUO4vLvS4nu1JDLa391XtCiRDDH\nXuDIoEc5NvIZ0hOuEtypJ/WWbiUt4vYct0lKMn61dOhgTHGxfDlUq+a8mG9kkmTEdale3ZgIrWNH\no7799tv2pw7IULbTI9RbspXApi0wnz/DociHOD72ZdKTEp0Wsyi+rqxbZQzOuvxHPAJKUfXdWVT7\n4Es8y+R85v5ffxknGR84YByH6d9fTq50Jkky4rqZTPDqq8bB04wDqL/9lnN77/Aq1J63nPBX30Z5\neXFu3jT2dr2ThD0yBZCwLz05mRP/G87B/g+QeiaGgEZ3U2/JFkK69c3xhN8rV+CFF4yOKhMnGmfx\nV7A3PK9wKEkyoshUq2aUzqKi4Lnn4OGH4fhx+22VyUSFga9S97u/8KlRh6SD0ezr0ZwzM96XSdGE\njcT9u9nX427OzvoATCYqvvwWdeavxKdKdbvttTamsLj1VuMM/t27jTlghGtIkhFFrlMno9fZrbfC\nHXcYMwjmMLIH/rfeQb2fNlHusYHo1BROTnqNg090JOX0SecGLdyO1pqzcz9hb9c7Sdy7E5+balL3\nm7VUHPRGjpOKrV8Pd99tvObmzoVZsyA450qacAJJMsIh/PxgzBjjJLe1a42TOLdssf9u9/Dz56Zx\nU6nx+Q94li1P3N+/E92pIbG/fO/coIXbSD13mkORD3Fi3BB0SjIhPZ8k4ufNBDRoarf90aMwbtzN\n9OwJzz5rlG1btXJqyCIHkmSEQ9WsaQy2+b//wYcf1qZNG+Pbpj1B93Wm3rJtlG71AGmXYzkyuDdH\nRzxFWnyc/Q1EiXRp5c9Ed2zIlT+WYwoqS/VPvqHqxM8xBWSfYOvKFXjtNWOcsapVr7Jvn3Fg30M+\n2dyG/CuEwylljIE2e/Zm+vSBRx81TurcsSN7W69yYdSc/iNV3voQ5ePLxe/nEv1gE67+s9H5gQun\nSku4yrFRL3D42R6YY89T6p77qLdkK8EdumVrazbD559DnTpw5kxGr7H/5KRKNyRJRjiNyaSJjIT9\n+6FtW+Ochd69jdvWlFKU7/ccET9uwK/ebaQcP8y+Xq049dF4dE4zqYliLWHXNvZ2acr5r6ejvLyp\n9Pq71PpiKd4VbGf1SEkxRgSPiDAO7i9bZpyxXynb5B/CXUiSEU7n62uMFnDgANx2G9xzDzz1FBw7\nZtvOr/bN1F24jtCnXoH0dE59NJ79vduQfOywawIXRU6npXH608nsfbg5yUcO4FvnFuou+puwAS+h\nrGpeiYnw8cdQqxZ8841xQH/1amjY0IXBi3yRJCNcJjDQGGhz/34ICzN6okVG2o6H5uHjQ+XXJlF7\n7q94hVXi6j8biH6wMRd+mCvjnxVzySf/40Df9sRMGQVmM+WfGEzEovX4R9yW2SYuDt59F2rUMCYS\nW7jQOAerRQsXBi4KRJKMcLngYGOkgH37jA+T9u2NUtry5cY5DwCl7mpNvaVbCXqgB+lX4/lv+FMc\nefExzJcuujZ4USgXf17A3s6Nid/8J57lK1Br1hKqvDkFDx9fAC5dgvHjjY4jW7car4Wffso+Q6tw\nf5JkhNsoVw7eeMMYA+2xx2D4cONcmxkzjLGnPIPKUv2jr6j6zgw8AgK59Mv3xqRo61e7OnSRT2lx\nlznySn+ODnmctLjLlGn7IPWWbKV0i/aA8St20CAjuRw+DH/+aRx7ue22PO5YuC1JMsLt+PgY3VD/\n/RemToUff4SqVY3zbs6eVYR0f5x6izcT0KAZqWdOcuDx+zkx6TXScxs0Tbhc/Oa/iO7cmNjFXxvn\nRk34lBqfLiQtoDzz5xsT4XXoACEhsH27cUC/bl1XRy2ulyQZ4baUgjZtjPNs1q6F06eND50ePWD5\n9ppUn7eaii+OAg8Pzs54n30PNyfxYLSrwxZZ6NRUTk4Zxf4+bUk5+R/+9RsRsXgTlxpHMmKE4qab\njNlWhw41TqocOxYqV3Z11KKoSJIRxUJEBHz2Gfz3H9x/vzFsSJVqnrxzZBQe41fjXaUGidHb2dul\nGefmfYrWmrSr8SQd2U/a1XhXh19g8fHxHD9+nPj44he79X5POrKffY+25Mynk0Fryj09gj2P/cFD\nA+tw991G+7//Ng7md+sGXl6ujV0UPUkyolgpUwaeftoYvn3dOmNa6C6v3slTFzdzKqI/OjmJ42Nf\nYler2mxvUpG9Xe9kR9Nwjo19uVicY2M2mxk8eDChoaEMHDiQ0NBQBg8ejLkYxK7NZo6NfZkdTcOJ\n7tKM7Y1C2XN/AxJ2bMEcXJVvb1nF7ZPG89EnXvTvb3RZf+cdo1uyKLnsjzInRDFQqxaMG2ccq1m7\nthRz5kzn4pX7GV56IP4nj5GuFaQYx2nOL5hJysljlO8zMNv9mHbt4LJyj+M5wz/8hHnLfiPR6vjS\nzOnTSTn5H++89IILI7smp/11bv7nXPlzBaQmk64VHsroGrgioScrw6bRpW0Z9nwKFSs6O2LhSpJk\nRLHn4WFMA926NVw+04EZjdpSRe2lls+ea41Sk7ny+xKu/L4k2/YBwCHnhZujhHTNnKNmkrOc/pOY\nnMycH3+m37+/4O/h+tm28rO/PJQmIT2AfxLvpkHprQz7zYRJhny5IUmSESWKT0IM94X8RvrVeNK0\nB3HpQcSmlediWnkUitBwT8IqehFcFjI+ry9evEjZsmVdGzhw4XI8puN/gTkt2zpPk4nEBvdQoUz2\nQSKdLWN/aSA+Ds5fgPNn07gUq/H3iCfEdIYg0zn8Pa5yT8AKPHwDST0bg6l6HVeHLlxAkowoUbxC\nwzOPvZhUOkGmiwSZLlJN7+OQasSOTmtY8qsP+/40hrNp3hwCAv7hmWfuwMfHtbFXiI9Hh4aCOftU\n1Ole3tw1dwmBga5NMidPwg9Tozlxoh4rVxqjNrRrB/c9ncjN4+pSynw62zY6LQ2v0HAXRCvcgRz4\nFyWKKSCQkEcGoHz9bZZ7+PlzZ79mvDnGhw0bjNEFnnzSGMH3449rERIC994Lr79uDLp46ZLzYw8M\nDCQyMhJ/f9vY/f39iYyMdHqCSU01zlWaMcOYo+WWW4yTIjdsCKFFC2PKhoMH4dNP4eHH/Kjeu3u2\n/a58/Qnp+aTdYfrFjcFtfskopXoCY4B6QFOt9ZYc2t0PfAiYgBla60mW5dWBBUBZYBvQT2udw3yM\noiSr8sZ7AFz47guUyYROSyOk5xOZywFCQ43zbXr0gC5dttKoUSs2bjTOMJ8yxZiOoEYNYxiT+vWN\nkQduvdXYzpGioqIAmDlzZuayAQMGZC53lLQ02LvXmOwr47JjhzGlduPGxiUy0hiQ8s8/99CqVfYd\nkZ/9Lm48bpNkgF1Ad+DznBoopUzAJ0A74ASwWSm1WGu9B5gMRGmtFyilPgMigU8dH7ZwN8rTk5ve\n+oBKwyaQejYGr9DwPL9JlyplTD/Qtq1xOzUVtm0zLjt3GgMz7toFnp7XEk7G5eabISioaGL39PRk\n6tSpTJw4kUWLFtGtW7ci/QUTGwuHDhm/QA4eNK7v328klIoVryWUnj2NAUtLlcr/fRdmv4uSz22S\njNY6Goy5RHLRFDiotT5sabsA6KKUigbaAI9Z2s3B+FUkSeYGZgoILPTBZi8vaNbMuGTQGk6dMpLN\nrl2wYYNRStqzB0wm4yz1rJdKla5dDw42RjHIj8DAQKpUqZLvBGM2w8WLcP687eXkyWsJ5eBBYz6W\nWrWuXZo3hyeeMMpgwfZnxy6w69nvouRxmySTT5WA41a3TwDNgBDgktbabLVcpjESRUopCA83Lu3b\nX1uutfEL4eRJOHHi2mXDBuNvxvKrV40D5bldAgKMLtlmM/z3X22++sq4nvWSnHwtqVy4YAyJHxxs\njPtVrty1S8WKxnhgL7xgJJXy5fOf6IQoCsqZc3IopVYCFeysekNr/ZOlzRpgmL1jMpbjNh201k9Z\nbvfD+HUzDlivta5lWV4FWKa1rp9DHM8AzwCEhYU1WrBgQaGeT3x8vMt7+9gjcRWMs+JKSVEkJZlI\nTLx2yXo7MdGE1sYsomZzEv7+3phM2ubi6WlcSpVKpUwZ4xIYaHbavPY3+v+xoEpqXK1bt96qtW6c\nZ0OttVtdgDVA4xzW3QUst7o90nJRwHnA01673C6NGjXShbV69epCb+tIElfBSFwFI3EVTEmNC9ii\n8/EZW9y6MG8GaiulqiulvIFewGLLE14NPGxp1x/4yUUxCiGEsHCbJKOU6qaUOoHxK2SpUmq5ZXm4\nUmoZgDaOuQwClgPRwLda692WuxgBvKKUOohxjGZm1scQQgjhXG5z4F9rvQhYZGd5DNDR6vYyYJmd\ndocxjs8IIYRwE27zS0YIIUTJI0lGCCGEw0iSEUII4TCSZIQQQjiMJBkhhBAOI0lGCCGEw0iSEUII\n4TCSZIQQQjiMJBkhhBAOI0lGCCGEwzh1qH93pJQ6B/xXyM3LYYz+7G4kroKRuApG4iqYkhpXVa11\n+bwa3fBJ5noopbbo/Myn4GQSV8FIXAUjcRXMjR6XlMuEEEI4jCQZIYQQDiNJ5vr8n6sDyIHEVTAS\nV8FIXAVzQ8clx2SEEEI4jPySEUII4TCSZPKglOqplNqtlEpXSuXYE0Mpdb9Sap9S6qBS6jWr5dWV\nUhuVUgeUUt8opbyLKK6ySqkVlvtdoZQKttOmtVLqX6tLklKqq2XdbKXUEat1DZwVl6VdmtVjL7Za\n7sr91UAptd7y/96hlHrUal2R7q+cXi9W630sz/+gZX9Us1o30rJ8n1Kqw/XEUYi4XlFK7bHsn1VK\nqapW6+z+T50U1xNKqXNWj/+U1br+lv/7AaVUfyfHFWUV036l1CWrdQ7ZX0qpWUqps0qpXTmsV0qp\njywx71BKNbRaV/T7Smstl1wuQD2gLrAGaJxDGxNwCKgBeAPbgZst674FelmufwY8V0RxvQO8Zrn+\nGjA5j/ZlgYuAv+X2bOBhB+yvfMUFxOew3GX7C6gD1LZcDwdOAUFFvb9ye71YtXke+MxyvRfwjeX6\nzZb2PkB1y/2YnBhXa6vX0HMZceX2P3VSXE8AH9vZtixw2PI32HI92FlxZWk/GJjlhP3VAmgI7Mph\nfUfgF0ABdwIbHbmv5JdMHrTW0VrrfXk0awoc1Fof1lqnAAuALkopBbQBFlrazQG6FlFoXSz3l9/7\nfRj4RWudUESPn5OCxpXJ1ftLa71fa33Acj0GOAvkebJZIdh9veQS70LgPsv+6QIs0Fona62PAAct\n9+eUuLTWq61eQxuAykX02NcVVy46ACu01he11rHACuB+F8XVG/i6iB47R1rrPzC+UOakCzBXGzYA\nQUqpijhoX0mSKRqVgONWt09YloUAl7TW5izLi0KY1voUgOVvaB7te5H9Bf625edylFLKx8lx+Sql\ntiilNmSU8HCj/aWUaorx7fSQ1eKi2l85vV7strHsj8sY+yc/2zoyLmuRGN+IM9j7nzozrh6W/89C\npVSVAm7ryLiwlBWrA79bLXbU/spLTnE7ZF95Xu8dlARKqZVABTur3tBa/5Sfu7CzTOey/Lrjyu99\nWO6nIlAfWG61eCRwGuOD9P+AEcA4J8Z1k9Y6RilVA/hdKbUTuGKnnav21zygv9Y63bK40PvL3kPY\nWZb1eTrkNZWHfN+3Uqov0BhoabU42/9Ua33I3vYOiOtn4GutdbJS6lmMX4Ft8rmtI+PK0AtYqLVO\ns1rmqP2VF6e+tiTJAFrrttd5FyeAKla3KwMxGOMCBSmlPC3fRjOWX3dcSqkzSqmKWutTlg/Fs7nc\n1SPAIq11qtV9n7JcTVZKfQEMc2ZclnIUWuvDSqk1wB3A97h4fymlSgNLgTctpYSM+y70/rIjp9eL\nvTYnlFKeQBmMEkh+tnVkXCil2mIk7pZa6+SM5Tn8T4viQzPPuLTWF6xuTgcmW23bKsu2a4ogpnzF\nZaUX8IL1Agfur7zkFLdD9pWUy4rGZqC2MnpGeWO8oBZr42jaaozjIQD9gfz8MsqPxZb7y8/9ZqsF\nWz5oM46DdAXs9kRxRFxKqeCMcpNSqhxwD7DH1fvL8r9bhFGv/i7LuqLcX3ZfL7nE+zDwu2X/LAZ6\nKaP3WXWgNrDpOmIpUFxKqTuAz4GHtNZnrZbb/Z86Ma6KVjcfAqIt15cD7S3xBQPtsf1F79C4LLHV\nxTiQvt5qmSP3V14WA49bepndCVy2fIlyzL5yRO+GknQBumFk+GTgDLDcsjwcWGbVriOwH+ObyBtW\ny2tgfAgcBL4DfIoorhBgFXDA8resZXljYIZVu2rAScAjy/a/AzsxPiy/BAKdFRdwt+Wxt1v+RrrD\n/gL6AqnAv1aXBo7YX/ZeLxjlt4cs130tz/+gZX/UsNr2Dct2+4AHivj1nldcKy3vg4z9sziv/6mT\n4poI7LY8/mogwmrbAZb9eBB40plxWW6PASZl2c5h+wvjC+Upy2v5BMaxs2eBZy3rFfCJJeadWPWa\ndcS+kjP+hRBCOIyUy4QQQjiMJBkhhBAOI0lGCCGEw0iSEUII4TCSZIQQQjiMJBkhhBAOI0lGCCGE\nw0iSEUII4TCSZIRwE0qpAKXUXqXUJqWUl9Xy9sqYNO+F3LYXwh3JGf9CuBHL2GAbgCit9WtKqVBg\nB7BJa/2Qa6MTouAkyQjhZpRSQ4ApGAMUDsOYpuF2rfV5lwYmRCFIkhHCzVhGel6KMR+KN9BOa73K\ntVEJUThyTEYIN6ONb37zAB9guyQYUZxJkhHCzSilKgAfANuA25VSL7k4JCEKTZKMEG7EUiqbA6QA\n7TCSzWSl1G0uDUyIQpJjMkK4EaXUUOAdoI3Weq1lxsUNGKWzxlrrRJcGKEQByS8ZIdyEpfvy/4CJ\nWuu1AFrrFIzps6sB77suOiEKR37JCCGEcBj5JSOEEMJhJMkIIYRwGEkyQgghHEaSjBBCCIeRJCOE\nEMJhJMkIIYRwGEkyQgghHEaSjBBCCIeRJCOEEMJh/h9YaKes3oko+wAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "N_panels = 10 # number of panels desired\n", "\n", "# define the end-points of the panels\n", "x_ends = R * numpy.cos(numpy.linspace(0.0, 2 * math.pi, N_panels + 1))\n", "y_ends = R * numpy.sin(numpy.linspace(0.0, 2 * math.pi, N_panels + 1))\n", "\n", "# define the panels\n", "panels = numpy.empty(N_panels, dtype=object)\n", "for i in range(N_panels):\n", " panels[i] = Panel(x_ends[i], y_ends[i], x_ends[i + 1], y_ends[i + 1])\n", " \n", "# plot the panels\n", "size = 6\n", "pyplot.figure(figsize=(size, size))\n", "pyplot.grid()\n", "pyplot.xlabel('x', fontsize=16)\n", "pyplot.ylabel('y', fontsize=16)\n", "pyplot.plot(x_cylinder, y_cylinder,\n", " label='cylinder',\n", " color='b', linestyle='-', linewidth=1)\n", "pyplot.plot(x_ends, y_ends,\n", " label='panels',\n", " color='#CD2305', linestyle='-', linewidth=2)\n", "pyplot.scatter([p.xa for p in panels], [p.ya for p in panels],\n", " label='end-points',\n", " color='#CD2305', s=40)\n", "pyplot.scatter([p.xc for p in panels], [p.yc for p in panels],\n", " label='center-points',\n", " color='k', s=40, zorder=3)\n", "pyplot.legend(loc='best', prop={'size':16})\n", "pyplot.xlim(-1.1, 1.1)\n", "pyplot.ylim(-1.1, 1.1);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Flow-tangency boundary condition" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In [Lesson 1](01_Lesson01_sourceSink.ipynb), you worked out the velocity potential of a single source as part of the final *Challenge Task.* It followed simply from integrating the radial velocity, $u_r=\\frac{\\sigma}{2\\pi r}$, that\n", "\n", "$$\\phi=\\frac{\\sigma}{2\\pi}\\ln r$$\n", "\n", "(The integration also gives a function of $\\theta$ that is seen to be a constant because $u_{\\theta}=0$; we take this constant to be zero.)\n", "\n", "We will use the velocity potential in this lesson to easily express that the velocity be tangent at the panel, i.e., that $u_n=0$ with:\n", "\n", "$$u_n(x,y)=\\frac{\\partial \\phi}{\\partial n}(x,y)$$\n", "\n", "at a given point on the panel. We choose the point to enforce that velocity-tangency condition as the center of the panel (and we call it the *control point*).\n", "\n", "The velocity potential in Cartesian coordinates of a [source sheet](08_Lesson08_sourceSheet.ipynb) on a panel is\n", "\n", "$$\\phi\\left(x,y\\right) = \\frac{\\sigma}{2\\pi} \\int_\\text{panel} \\ln \\sqrt{ \\left(x-x(s)\\right)^2 + \\left(y-y(s)\\right)^2 } {\\rm d}s$$\n", "\n", "where $s$ is the running coordinate along the panel and $\\left(x(s),y(s)\\right)$ are the Cartesian coordinates of $s$.\n", "\n", "Superposition of the potential of each panel gives the total potential at any point $\\left(x,y\\right)$, so we make a sum of all the panel contributions as follows (moving the $\\frac{1}{2}$ exponent in the logarithmic term as a factor outside the integral):\n", "\n", "$$\\phi\\left(x,y\\right) = \\sum_{j=1}^{N_p} \\frac{\\sigma_j}{4\\pi} \\int \\ln \\left( \\left(x-x_j(s_j)\\right)^2 + \\left(y-y_j(s_j)\\right)^2 \\right) {\\rm d}s_j$$\n", "\n", "By finally superposing the free stream, the flow around an immersed circular cylinder will be represented by the following velocity potential:\n", "\n", "$$\\phi\\left(x,y\\right) = U_\\infty x + \\sum_{j=1}^{N_p} \\frac{\\sigma_j}{4\\pi} \\int \\ln \\left( \\left(x-x_j(s_j)\\right)^2 + \\left(y-y_j(s_j)\\right)^2 \\right) {\\rm d}s_j$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Enforcing the flow-tangency condition on each *control point* approximately makes the body geometry correspond to a dividing streamline (and the approximation improves if we represent the body with more and more panels). So, for each panel $i$, we make $u_n=0$ at $(x_{c_i},y_{c_i})$:\n", "\n", "$$u_{n_i} = \\frac{\\partial}{\\partial n_i}\\left\\lbrace \\phi\\left(x_{c_i},y_{c_i}\\right) \\right\\rbrace = 0$$\n", "\n", "which leads to\n", "\n", "$$\n", "0 = U_\\infty \\cos\\beta_i + \\sum_{j=1}^{N_p} \\frac{\\sigma_j}{2\\pi} \\int \\frac{\\left(x_{c_i}-x_j(s_j)\\right) \\frac{\\partial x_{c_i}}{\\partial n_i} + \\left(y_{c_i}-y_j(s_j)\\right) \\frac{\\partial y_{c_i}}{\\partial n_i}} {\\left(x_{c_i}-x_j(s)\\right)^2 + \\left(y_{c_i}-y_j(s)\\right)^2} {\\rm d}s_j\n", "$$\n", "\n", "where $\\beta_i$ is the angle that the panel's normal makes with the $x$-axis, so\n", "\n", "$$\\frac{\\partial x_{c_i}}{\\partial n_i} = \\cos\\beta_i \\quad \\text{and}\\quad\\frac{\\partial y_{c_i}}{\\partial n_i} = \\sin\\beta_i$$\n", "\n", "and\n", "\n", "$$x_j(s_j) = x_{a_j} - \\sin\\left(\\beta_j\\right) s_j$$\n", "$$y_j(s_j) = y_{a_j} + \\cos\\left(\\beta_j\\right) s_j$$\n", "\n", "But, there is still a problem to handle when $i=j$. From the previous notebook, we have seen that the strength of the [source sheet](08_Lesson08_sourceSheet.ipynb) should be a specific value so that the streamlines do not penetrate the panel. This helps us determine that the contribution of the $i$-th panel to itself is $\\frac{\\sigma_i}{2}$.\n", "\n", "Finally, the boundary condition at the center point of the $i$-th panel gives\n", "\n", "$$\n", "0 = U_\\infty \\cos\\beta_i + \\frac{\\sigma_i}{2} + \\sum_{j=1,j\\neq i}^{N_p} \\frac{\\sigma_j}{2\\pi} \\int \\frac{\\left(x_{c_i}-x_j(s_j)\\right) \\cos\\beta_i + \\left(y_{c_i}-y_j(s_j)\\right) \\sin\\beta_i} {\\left(x_{c_i}-x_j(s)\\right)^2 + \\left(y_{c_i}-y_j(s)\\right)^2} {\\rm d}s_j\n", "$$\n", "\n", "From the equation above, we understand that we will have to compute integrals using the SciPy function `integrate.quad()`. We define a function `integral_normal()` that will do the job." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def integral_normal(p_i, p_j):\n", " \"\"\"\n", " Evaluates the contribution of a panel at the center-point of another,\n", " in the normal direction.\n", " \n", " Parameters\n", " ----------\n", " p_i: Panel object\n", " Panel on which the contribution is calculated.\n", " p_j: Panel object\n", " Panel from which the contribution is calculated.\n", " \n", " Returns\n", " -------\n", " Integral over the panel at the center point of the other.\n", " \"\"\"\n", " def integrand(s):\n", " return (((p_i.xc - (p_j.xa - math.sin(p_j.beta) * s)) * math.cos(p_i.beta) +\n", " (p_i.yc - (p_j.ya + math.cos(p_j.beta) * s)) * math.sin(p_i.beta)) /\n", " ((p_i.xc - (p_j.xa - math.sin(p_j.beta) * s))**2 +\n", " (p_i.yc - (p_j.ya + math.cos(p_j.beta) * s))**2))\n", " return integrate.quad(integrand, 0.0, p_j.length)[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solving the system of equations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We just developed an equation to enforce a flow-tangency condition on the $i$-th panel. There are `N_panels` panels $i$ and `N_panels` unknown strengths $\\sigma_i$. Therefore, the problem represents solving a linear system of equations of the form\n", "\n", "$$[A][\\sigma] = [b]$$\n", "\n", "where\n", "\n", "$$\n", "A_{ij} = \\begin{cases}\n", "\\begin{matrix}\n", "\\frac{1}{2} & \\mbox{, if } i=j \\cr\n", "\\frac{1}{2\\pi} \\int \\frac{\\left(x_{c_i}-x_j(s_j)\\right) \\cos\\beta_i + \\left(y_{c_i}-y_j(s_j)\\right) \\sin\\beta_i} {\\left(x_{c_i}-x_j(s)\\right)^2 + \\left(y_{c_i}-y_j(s)\\right)^2} ds_j & \\mbox{, if } i\\neq j\n", "\\end{matrix}\n", "\\end{cases}\n", "$$\n", "\n", "and\n", "\n", "$$b_i = - U_\\infty \\cos\\beta_i$$\n", "\n", "for $1\\leq i,j \\leq N_p$. Let's fill a matrix $A$ and a right-hand side vector $b$ with the necessary values:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# compute the source influence matrix\n", "A = numpy.empty((N_panels, N_panels), dtype=float)\n", "numpy.fill_diagonal(A, 0.5)\n", "\n", "for i, p_i in enumerate(panels):\n", " for j, p_j in enumerate(panels):\n", " if i != j:\n", " A[i, j] = 0.5 / math.pi * integral_normal(p_i, p_j)\n", "\n", "# compute the RHS of the linear system\n", "b = - u_inf * numpy.cos([p.beta for p in panels])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Hey! We just used a new Python built-in function: [enumerate()](https://docs.python.org/2/library/functions.html#enumerate). It allows us to have access to each element `panel` in the array `panels` while keeping a count `i` (that starts from `0`) to locate the element of `A` to fill." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we can easily solve the linear system of equations using the function [`linalg.solve()`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.solve.html) from NumPy, and assign each source-panel its appropriate strength:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# solve the linear system\n", "sigma = numpy.linalg.solve(A, b)\n", "\n", "for i, panel in enumerate(panels):\n", " panel.sigma = sigma[i]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Pressure coefficient on the surface" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "At this point, we have the source strength distribution required to compute the streamlines around our geometry. A very useful measurement of the results is the pressure coefficient along the surface of the geometry.\n", "\n", "From Bernoulli's equation, the pressure coefficient on the $i$-th panel is\n", "\n", "$$C_{p_i} = 1-\\left(\\frac{u_{t_i}}{U_\\infty}\\right)^2$$\n", "\n", "where $u_{t_i}$ is the tangential component of the velocity at the center point of the $i$-th panel,\n", "\n", "$$u_{t_i} = \\frac{\\partial}{\\partial t_i}\\left\\lbrace \\phi\\left(x_{c_i},y_{c_i}\\right) \\right\\rbrace$$\n", "\n", "which we can obtain as:\n", "\n", "$$\n", "u_{t_i} = -U_\\infty \\sin\\beta_i + \\sum_{j=1}^{N_p} \\frac{\\sigma_j}{2\\pi} \\int \\frac{\\left(x_{c_i}-x_j(s_j)\\right) \\frac{\\partial x_{c_i}}{\\partial t_i} + \\left(y_{c_i}-y_j(s_j)\\right) \\frac{\\partial y_{c_i}}{\\partial t_i}} {\\left(x_{c_i}-x_j(s)\\right)^2 + \\left(y_{c_i}-y_j(s)\\right)^2} {\\rm d}s_j\n", "$$\n", "\n", "with\n", "\n", "$$\\frac{\\partial x_{c_i}}{\\partial t_i} = -\\sin\\beta_i \\quad\\text{and} \\quad \\frac{\\partial y_{c_i}}{\\partial t_i} = \\cos\\beta_i$$\n", "\n", "Note that the contribution to the tangential velocity at a source panel from its own velocity potential is zero, which makes sense because streamlines go *outwards* from a source.\n", "\n", "We define a function `integral_tangential()` that will compute the integrals above using the SciPy function `integrate.quad()` once again:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def integral_tangential(p_i, p_j):\n", " \"\"\"\n", " Evaluates the contribution of a panel at the center-point of another,\n", " in the tangential direction.\n", " \n", " Parameters\n", " ----------\n", " p_i: Panel object\n", " Panel on which the contribution is calculated.\n", " p_j: Panel object\n", " Panel from which the contribution is calculated.\n", " \n", " Returns\n", " -------\n", " Integral over the panel at the center point of the other.\n", " \"\"\"\n", " def integrand(s):\n", " return ((-(p_i.xc - (p_j.xa - math.sin(p_j.beta) * s)) * math.sin(p_i.beta) +\n", " (p_i.yc - (p_j.ya + math.cos(p_j.beta) * s)) * math.cos(p_i.beta)) /\n", " ((p_i.xc - (p_j.xa - math.sin(p_j.beta) * s))**2 +\n", " (p_i.yc - (p_j.ya + math.cos(p_j.beta) * s))**2))\n", " return integrate.quad(integrand, 0.0, p_j.length)[0]" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# compute the matrix of the linear system\n", "A = numpy.empty((N_panels, N_panels), dtype=float)\n", "numpy.fill_diagonal(A, 0.0)\n", "\n", "for i, p_i in enumerate(panels):\n", " for j, p_j in enumerate(panels):\n", " if i != j:\n", " A[i, j] = 0.5 / math.pi * integral_tangential(p_i, p_j)\n", "\n", "# compute the RHS of the linear system\n", "b = - u_inf * numpy.sin([panel.beta for panel in panels])\n", "\n", "# compute the tangential velocity at each panel center-point\n", "vt = numpy.dot(A, sigma) + b\n", "\n", "for i, panel in enumerate(panels):\n", " panel.vt = vt[i]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Once we have computed the tangential velocity on each panel, we can calculate the pressure coefficient." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# calculate the surface pressure coefficient\n", "for panel in panels:\n", " panel.cp = 1.0 - (panel.vt / u_inf)**2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Alright! It is time to plot the surface pressure coefficient.\n", "\n", "Just before this, we should remember that in the lesson on the [doublet](03_Lesson03_doublet.ipynb), we found that the exact pressure coefficient on the surface on a cylinder was\n", "\n", "$$Cp = 1 - 4\\sin^2 \\theta$$\n", "\n", "i.e.\n", "\n", "$$Cp = 1 - 4\\left(\\frac{y}{R}\\right)^2$$\n", "\n", "We can use this to compare with the results obtained with our source-panel code." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnUAAAGKCAYAAABq27cFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3Xd4VFX+x/H3AQIJhE4ABWkCoSog\nIiBId0VQRERXQUAWsaGirLuyKAgW/FnWgoLYwLqKCjYQ1gKKFEVQWVoAadI7IUAgJOf3x0lCQiYh\nCZPcKZ/X8+S5M3fu3PuZk8nMN/fec66x1iIiIiIiwa2I1wFERERE5OypqBMREREJASrqREREREKA\nijoRERGREKCiTkRERCQEqKgTERERCQEq6kRCiDFmkDHGGmMOGmPKn/ZYsdTHHvEg1yOp2y5W2NvO\nC2NMEWPM88aYHcaYFGPMp15nyi9jzDxjzLxC3F5pY8wzqduNT/19d8xm2SLGmJHGmE3GmERjzO/G\nmD6FlVUkVKmoEwlNZYF/eh0iCF0H3As8DVwK/MPbOEGlIjAYOAl8fYZlHwUeAV4CugOLgY+MMVcW\nZECRUBfQ/zWLSL79F7jbGPO8tXan12EKgzGmhLX2+FmupmHq9HlrbcrZZgozm621FQCMMV2Ba30t\nZIypDPwdeNJa+0zq7LnGmLrAk8CswggrEoq0p04kND2WOh2V00Jph0V9zJ9qjNmU4X6t1MNptxtj\nxhtjdhpjDhtj3jXGlDTG1DXGzDHGJBhj1htjBmazyYbGmLnGmKOphzjHGWMyfQ4ZYyoZYyYZY7YZ\nY44bY9YYY4aetkzaYebLjDEfGWMOAj+d4bVeYYxZZIw5Zow5ZIz51BgTm+HxTbi9RwDJqesflMP6\nrDHmcWPMKGPM1tT1/mCMaXbacpcbY2alvt6jxpgVxpgRxpiipy23KbU9/2qMWW2MOWKM+cUY087H\ntjsYY75N/R0cSW37Jmd4/dHGmAnGmC2p7brLGPONMaZBTs/LLZv7yxP9BSgOvHva/HeBpsaY2v7I\nIxKOVNSJhKYduENbQ40xNf243pHAucBAYDRwA/AKMAOYCfQGlgNTjDGNfTz/U+Ab4BrgfeDh1PUA\nYIwpAywAeuAKrB7AF8AkY8zdPtb3HrARd9j0wexCG2OuSM2XkJr5DqAJ8KMxplrqYr2Bqam326T+\nzMxunakGAFcCw4BBQBXgW2NMhQzL1AG+xR2a7AG8lfraHvexvvbACFy73AAUBb40xpTL8Fp6pK4v\nAegP3ASUBuYbY87LIetzwPXAWKAbcDvwG1Auh+dkPB+yVk7L5UFj4Diw/rT5K1Onjfy0HZGwo8Ov\nIqHr/4DbgDG4gsIf/rDWpu2Fm2OMaQ/cDNxsrX0XwBjzC3A1rtBaedrzX7PWPpl6+7+pRdyI1MPE\nB3Hns9UEmlpr16Uu901qUTPGGDPJWnsyw/o+ttbm5ry3x4ANQPe05xtjFgFrcUXU/dbaX40x2wCs\ntYtz2R5RwOXW2iOp6/wJWAfchyvMsNa+krawMcYA83F7qv5ujPnXaYd5ywDNrLUHUpffCSzBFY7v\npy7zAvC9tbZXhvXOTX19I4Dh2WRtA7xnrX0jw7wZuXiNKUAy4K8LhVcADvrYs7c/w+Mikg/aUycS\noqy1+4FngQEZDzOepa9Ou78mdTonw3YPALsBX3uNpp12/wMgGrfXDOAK3GHUjcb11i1mXI/ZObgT\n8U/fi3PGosQYUwpoAXyYsSC01m7E7RXscKZ15GBWWkGXus5NuJP+22TY/jnGmMnGmM3ACSAJV2SW\nAyqftr5FaQVdqv+lTmukrqsecD7w3mntcxRYBFyWQ9YlwCBjzL+MMS1PP/ybHWvtOGttMWvt5tws\nnwsG3wWi8dP6RcKWijqR0PYcbg/IOD+t78Bp90/kMD/Sx/N3ZXM/7RBoZVxhknTaz0epj1c87fk7\nzhyZ8riCwdeyOzm7PUOnv560edXADd0BfA70xBVynYGLOXXo9fQ22p/xToaOH2nLpRWBb5C1jXqS\ntX0yuhuYjNtruwTYbYx5zhhTMofnFIT9QPnUvZYZlc/wuIjkgw6/ioQwa22CMWY8bo/d0z4WSQQw\nxhS31p7IMD+n4uBsVMEdJsx4H2Bb6nQfbi/fvdk8P+60+7k5JHggdbmqPh6rmrrN/KqSzby013M+\n0JIMh6cBjDFX5XN7aVlH4s5NPN0JH/MA915Ifd7I1PMsr8P1Nj1B4Q5/sxIogWubjOfVpe2FXVWI\nWURCivbUiYS+ibgi4zEfj6UdUkvvOZl6/lrbAspy/Wn3/4o74X9F6v3ZQANgi7X2Fx8/h/O6wdTD\no0uBvhkPOaYWNm2B7/PzQlJdmXp4N22dtYDWuEOhAGl7wZIyLBMB9Mvn9uKATUDjbNpneW5WYq3d\nbK19Fnd4N8deswVgNq6QPL0N+gMrUg+Li0g+aE+dSIiz1h43xowDXvXx8FfAIeA1Y8wY3B6Uf+AK\nrYJwa+ohySW4oS2GAI+kdpIAd7j4BlxPzudwRUwpXKHXPmPngDx6GNeT9UtjzETceXxjca/92fy+\nGOAYrsPH07i2GwvEp74OgNW4wvlxY0wyrri7L78bs9ZaY8xdwGfGmOK4cxT34vYOtsUVw//29dzU\njiGf4wq5BNy5hBfieuNmyxgzGtdD+fwznVdnjOmO+301TZ3VwRhTCThirf0q9TXsTv3djjTGHAaW\n4X7nnYH8/n5FBBV1IuFiCvAAUC/jTGvtQWNMT1wRMg3Yijv/rivQsQBy9AIm4IqsQ7i9h49myHPI\nGNMWV0T8E3du2kFccfdJfjdqrZ2dOhTIGNzrPAHMA/5hrd2e3/UCbwNHcMPHVMIVq39N7aSCtfaE\nMeaa1Mffxp0v9iawBXgtn69lljHmMtwYhK/jeuDuxHXQ+DCHp/6A21P6IO6zfwNwn7X2xTNssghu\naJXcdGSYhOu9nOaR1OlmoFaG+aNwheW9uEPgccD11tovcrENEcmGyf14kSIiksa4QZsft9Y+5HUW\nERHQOXUiIiIiISEoijpjzHnGXVpotTFmpTEmu55xIiIiImEpKA6/GmPOAc6x1i4zxpTG9WS7xlqr\nru8iIiIiBMmeOmvtDmvtstTbh3E9yqrl/CwRERGR8BEURV1GqeNANcddSkhERERECLIhTYwx0bhh\nDYZba+N9PD4UGAoQGRl5UY0aNQo5YWBLSUmhSJGgq+MLnNrFN7WLb2qXrNQmvqldfFO7+LZ27dq9\n1tqYs1lHUJxTB+mjsH8JzMlucM2MYmNjbVzc6VcUCm/z5s2jY8eOXscIOGoX39QuvqldslKb+KZ2\n8U3t4psxZqm1tuXZrCMoSuXUCz+/AazOTUEnIiIiEm6CoqgDLgVuBjobY35L/bnS61AiIiIigSIo\nzqmz1v5I7i5RIyIiIhKWgmVPnYiIiIjkQEWdiIiISAhQUSciIiISAoLinDoREcm7+Ph4du/eTVJS\nUoFup2zZsqxevbpAtxGM1C6+hVO7FCtWjMjISGJiYoiMjCz47RX4FkREpNDFx8eza9cuqlWrRlRU\nFG5kqIJx+PBhSpcuXWDrD1ZqF9/CpV2stZw8eZKEhAS2bNlClSpVKFu2bIFuU0WdiEgI2r17N9Wq\nVaNkyZJeRxEJS8YYIiIiKF++PCVKlGDnzp0FXtTpnDoRkRCUlJREVFSU1zFEBIiKiuL48eMFvh0V\ndSIiIaogD7mKSO4V1t+iijoRERGREKCiTkRERCQEqKgTERHJxiOPPJKvQ2effvopL730Upb58+bN\nwxjDvHnz/JAus02bNmGMYerUqX5ftwQHFXUiIiJ+ll1R16JFCxYtWkSLFi08SCWhTkOaiIiIFJIy\nZcrQunVrr2NIiNKeOhERCXjr16/n5ptvpnbt2kRFRVGnTh3uuOMODhw4kGm5QYMGUb16dX799Vfa\nt29PyZIlqVevHq+88kqm5fbs2cNtt91G/fr1KVmyJOeddx433XQT27ZtyzFH06ZN6d27d5b5aYdV\n58yZw6BBg3jrrbfYvn07xhiMMdSqVSvTcqcffp0xYwaXXnop0dHRlClThlatWvH555+nP/7SSy/R\npk0bKlSoQLly5WjdujUzZ87MQwtKONCeOhERCXjbt2+nevXqPP/885QvX54NGzbwxBNPcOWVV7Jo\n0aJMy8bHx3PTTTcxfPhwRo8ezZQpU7jjjjuIjY2lU6dOAOzfv5/IyEjGjx9PTEwM27dv59lnn+XS\nSy9lzZo12V7S6Y477uDee+9l+/btnHvuuenzJ0+eTO3atbn88supW7cue/bs4eeff+aLL74AoESJ\nEtm+tgkTJnDPPfdwzTXX8NZbbxEdHc2yZcvYtGlT+jKbNm1iyJAh1KpVi5MnT/LFF1/Qs2dPZs2a\nRffu3fPbrBJqrLUh+VO/fn0rmc2dO9frCAFJ7eKb2sW3YGmXVatW+ZwP3v/4Q1JSkp0/f74F7LJl\ny9LnDxw40AL2u+++S5+XmJhoK1asaG+99dZs13fy5Em7ZcsWC9jp06enzx8zZowlQ+j4+HhbunRp\nO27cuPR5e/bsscWLF7fjx4/PlOPcc8/Nsp25c+daIP19dOjQIRsdHW179+6d69eenJxsk5KSbLdu\n3ezVV1+dPn/jxo0WsFOmTMn1urwQHx/vdQRPZPc3mQb4xZ5l7aPDryIiYaQgyrT4+MN5Wj4/Tpw4\nwRNPPEGDBg2IiooiIiKC9u3bAxAXF5dp2ZIlS6bvkQO3l6xevXps2bIl03KTJk3iwgsvJDo6mmLF\nilGjRg2f68uodOnS9O/fn9dff52UlBQApkyZgrWWW265Jc+va+HChSQkJDB06NAcl1u6dCk9e/ak\nSpUqFCtWjIiICL7++uscs0r4UVEnIiIBb+TIkTzyyCP079+fmTNn8vPPPzN9+nQAEhMTMy1bvnz5\nLM8vUaJEpuUmTJjAnXfeSdeuXZk+fTo///wzixcv9rm+0915551s2bKFWbNmYa3l1VdfpXfv3lSp\nUiXPr2vfvn0AVK9ePdtl/vzzT7p06cL+/fuZMGECCxcuZMmSJVxxxRVnzCrhRefUiYhIwPvggw8Y\nMGAADz30UPq8hISEs1pfly5dePbZZ9Pnbdy4MVfPbdKkCe3bt2fy5MlERkayfv16Jk+enK8clSpV\nAmDbtm00adLE5zKzZ8/m0KFDTJs2LVPxd/To0XxtU0KX9tSJiEjAO3r0KBEREZnmTZkyxbP13Xnn\nnXz11Vc88sgj1K9fn86dO2d6vESJEhw7duyM62nbti3R0dG8+uqrOWYFMuVdu3YtCxYsyHVeCQ/a\nUyciIgHviiuu4K233qJp06bUrVuX6dOns3DhwrNa3//93//xxBNP0KpVK7777js+/vjjXD+/T58+\nDB8+nAULFmTa25emUaNGHDhwgEmTJtGyZUsiIyNp2rRpluVKly7N+PHjufvuu+nTpw/9+vWjdOnS\n/Pbbb0RGRnL33XfTtWtXihUrxoABAxgxYgQ7duxgzJgx1KhRI/28PhFQUSciIkFgwoQJWGsZNWoU\nAFdeeSX/+c9/aNWqVb7WN3r0aA4ePMhzzz1HYmIiHTp0YM6cOdSpUydXz4+IiKBXr1689dZbDBw4\nMMvjQ4YMYf78+fzrX//i4MGD1KxZM9MQJRkNGzaMqlWr8vTTT9OvXz8iIiJo2LAhDz/8MACNGzfm\nvffeY/To0Vx99dWcf/75PPnkk8yePbtALjcmwcvY/HZFCnCxsbFWvYIymzdvHh07dvQ6RsBRu/im\ndvEtWNpl9erVNGzYsFC2dfjwYUqXLl0o2woUJ0+epG7durRv35533nnH5zLh2C65Ea7tcqa/SWPM\nUmtty7PZhvbUiYiI5FJ8fDwrVqzg/fff588//2TEiBFeRxJJp6JOREQkl5YtW0anTp2oXLkyL7zw\nAs2aNfM6kkg6FXUiIiK51LFjR0L1tCUJfhrSRERERCQEqKgTERERCQEq6kRERERCgIo6ERERkRCg\nok5EREQkBKioExEREQkBKupEREREQoCKOhEREcmTWrVqMWjQIK9jAPD8888zffr0LPMfeeQRjDGc\nPHmyQLffsWPHgLl0oIo6ERERCVrZFXXhSEWdiIhILhw/ftzrCCI5UlEnIiK5knwkgcSNa0k+klDo\n2167di29e/emcuXKREZGUqNGDfr27Zvp0FpcXBy9e/emXLlyREVF0bp1a2bPnp1pPYMGDaJWrVpZ\n1n/6IbR58+ZhjGH69OnceuutxMTEUKVKlfTHf//9d3r37k3FihWJiooiNjaW8ePHZ1rn559/TuvW\nrSlZsiTlypWjb9++bNmy5YyvddOmTRhjmDhxIvfffz+VK1emZMmS9OzZk02bNmVa9oMPPqBz587E\nxMQQHR1N8+bNeeutt7Ks0xjDQw89xIsvvkjt2rUpXbo0HTp0YOXKlVmWnT59er5yny6tDT/99FNu\nu+02KlSoQPny5XnwwQdJTk5myZIltGvXjlKlStG4cWPmzJmTZR3ff/89Xbp0oXTp0pQqVYq//OUv\nrFixIv3xWrVqsXnzZt577z2MMRhjshwW3rhxIz169CA6OpqaNWsybtw4UlJSMi2Tm/cOuPZu0KAB\nJUqUoHHjxsyYMSPP7VKQVNSJiEiO7MmTbBk7nOWtzmXNNa1Z3upctowdji3gc5Uy6tmzJ9u2bWPS\npEnMmTOHJ598khIlSqR/OW/fvp127drx+++/89JLLzFt2jTKlStHjx49+Oqrr/K93bvvvhtrLe+8\n8w5Tp04F4Oeff6ZNmzb88ccfPPfcc8ycOZP777+frVu3pj/vlVdeoX///jRq1IiPP/6YyZMns2LF\nCjp06MDhw4dzte3x48ezbt06pkyZwssvv8zSpUu5/PLLSUpKSl9mw4YNXHfddbz33nt8+umnXHXV\nVQwZMoRXXnkly/reffddZs6cyQsvvMCUKVPYsmULvXr1ylQYv/LKK/Tp0+escp9u+PDhlCpVig8/\n/JBhw4YxceJEhg8fzoABAxg8eDDTp0+nQoUKXHvttezduzf9eTNnzqRLly5ER0fz7rvv8v7773P4\n8GHat2/Pn3/+CcCMGTOoWrUqf/nLX1i0aBGLFi3i4YcfzrT93r1707lzZz799FOuueYaxowZk6nw\nze1755tvvuGmm26iXr16TJ8+nQceeIB7772XuLi4fLVLgbDWhuRP/fr1rWQ2d+5cryMEJLWLb2oX\n34KlXVatWuW3dW1+5F67rHFZu/T8iPSfZY3L2s2P3GuttTY+Pt5v2/Jlz549FrCfffZZtsuMGDHC\nFi1a1K5bty593smTJ239+vVt8+bN0+cNHDjQ1qxZM8vzO3ToYDt06JB+f+7cuRaw11xzTZZl27dv\nb6tXr26PHDniM8vhw4dtmTJlbP/+/TPN37hxo42IiLDPPfdctq8jbTnANmzY0CYnJ6fP//HHHy1g\nX3/9dZ/PS05OtklJSXbIkCH2ggsuyPQYYOvWrWtPnDiRPu+jjz6ygF2wYEGm3LfccssZc9esWdMO\nHDgwx9eR1oanr+/CCy+0gJ0/f376vN9//90CdurUqenzzj//fNu5c+dMzz106JCtWLGivffeezNl\n6devX5btjxkzxgL2zTffzDS/SZMmtlu3bun3c/veadu2bZbfyeLFiy2Q6b2TnTP9TQK/2LOsfbSn\nTkREspV8JIF9097EJh7NNN8mHmXfR1MK5VBsxYoVqVOnDg8++CCvvfYa69aty7LMDz/8QOvWralb\nt276vKJFi3LjjTfy22+/ER8fn69t9+7dO9P9o0ePsmDBAvr160fJkiV9PmfRokXEx8dz/fXXc/Lk\nyfSf6tWr06BBA3744QcAUlJSMj1++iHB6667jiJFTn1NX3rppVSvXp1Fixalz1u3bh033ngj1apV\nIyIigoiICF5//XWfe4+6detGRERE+v2mTZsCpB9aTcvdr1+/HHPnVffu3TPdr1evHqVKlaJdu3bp\n8xo0aACQvgdu3bp1/PHHH1mylCxZkjZt2uQpS48ePTLdb9KkSabDybl576QdLj79d3LJJZf4PJzv\nFRV1IiKSraTd2zHFivl8zBQtStLu7QWewRjD119/TcuWLRk5ciT169enTp06TJo0KX2Z/fv3c845\n52R5btWqVbHWcuDAgXxt+/R1HjhwgJSUFKpXr57tc3bv3g3A1VdfnV5opf3873//Y9++fQAMHjw4\n02ODBw/OtJ6M5/BlnLdt2zYAEhIS6NatG7///jtPPvkk8+fPZ8mSJQwePNhnp44KFSpkul+iRAkA\nEhMTM+Xu2rVrjrnzqnz58pnuFy9enHLlymWZ5yvL3/72tyxZvvzyyzxl8fW607YDuXvv7N27l6Sk\npGx/J4HC91+qiIgIEFH53GzPnbPJyURUPpekFFvgOerUqcPbb7+NtTb93Kc777yTWrVq0b17dypU\nqMDOnTuzPG/nzp0YY9K/2CMjIzlx4kSW5fbt20fFihWzzDfGZLpfvnx5ihQpkl5Y+ZK2nkmTJtGy\nZcssj5cuXRpw46gNGzYsfX6lSpUyLbdr164sz921axfNmjUD3J61zZs3M3/+/Ex7vfI7Llta7qlT\np9K4ceNscxeGtCzjx4+na9euWR5PKwL9ITfvnZIlSxIREZHt76RmzZp+y3M2tKdORESyVbRUNBWv\nH4yJzHyo0USWpGLfWyhaKrpQ8xhjaNasGf/+978B0ntCdujQgcWLF2fqHZqcnMyHH35I8+bN0wuS\nmjVrsmvXrkwn5P/xxx+5Ptm9ZMmStGvXjnfffZdjx475XKZt27aULl2aDRs20LJlyyw/sbGxgOu5\nmXH+6YfxPv7440yHZBcsWMDWrVtp06YN4A4FA5kOqR44cIDPPvssV68lu9zr16/PMXdhiI2NpVat\nWqxcudJnlgsuuCB92RIlSmT7u8iN3Lx3ihYtysUXX5zld/LTTz9l6ZHsJe2pExGRHJ036hkA9n00\nBVO0KDY5mYp9B6XPL2jLly/n3nvv5YYbbqBu3bokJyczdepUihUrRufOnQG47777mDp1Kt26dWPs\n2LGUKVOGiRMnsnbtWmbOnJm+rr59+/Lwww/Tr18/7r//fvbu3cv48eOz7CXLyTPPPEOHDh1o06YN\nI0aMoHr16mzYsIHffvuNCRMmUKZMGZ5++mnuuusu4uPj6d69O2XLlmXbtm18//33dOzYkZtuuumM\n2zl8+DDXXHMNt912G3v27GHkyJHUq1ePAQMGAK4IK1OmDHfddRdjx47lyJEjPPbYY1SqVIlDhw7l\nsZXJlHvPnj35zu0PxhhefvllevXqxYkTJ7j++uupVKkSu3btYuHChdSoUYP7778fgEaNGjF//ny+\n/PJLqlatSqVKlfJ0nltu3ztjx47l8ssvz/Q7GTNmDFWrVvX3y8+/s+1pEag/6v2aVbD02itsahff\n1C6+BUu7+LP3a5qTCYftsQ1x9mTC4UzzC7r3665du+yAAQNsvXr1bFRUlC1fvry97LLL7OzZszMt\nt2bNGturVy9bpkwZW6JECXvJJZfYr776Ksv6ZsyYYRs3bmwjIyPtBRdcYOfMmZNt79evv/7aZ6Zl\ny5bZnj172rJly9rIyEgbGxtrn3zyyUzLfPTRR7Zjx462dOnSNjIy0p5//vn2lltusStXrszx9ab1\nfn355ZftfffdZytVqmSjoqLslVdeaTds2JBp2W+//dY2a9bMRkZG2jp16tgXXnghvddnRoAdNWqU\nz+1MmTIl0/yZM2eeMXdeer+e3oY33XSTrVatWpblfWVcuHCh7dGjhy1XrpwtUaKErVmzpr3hhhvs\nwoUL05dZvXq1bdeunY2KirJAeq60dkhKSsq0Tl89oHP73nn//fdt/fr1bfHixW2jRo3s9OnTs7x3\nslMYvV+NW0/oiY2NtQE1dkwAmDdvXsBcny6QqF18U7v4Fiztsnr1aho2bFgo2zp8+HChnm8VLPLb\nLps2baJ27dq89tprDBkypACSeStc3y9n+ps0xiy11mY9CTMPguKcOmPMm8aY3caYFWdeWkRERCT8\nBEVRB0wFrvA6hIiIiEigCoqOEtbaH4wxtbzOISIiUtBq1apFqJ4aJVlZCz6GycuXoDmnLrWo+9Ja\n2ySHZYYCQwEqVYq56KOPphVOuCCRkJBAdHThDj8QDNQuvqldfAuWdilbtmymEfILUnJyMkWLFi2U\nbQUTtYtv4dou69ev99kree1aOHwY/v73Tmd9Tl1Q7KnLLWvtq8CrAOedF2uPHu3IlVd6HCqABMsJ\n3oVN7eKb2sW3YGmX1atXEx0dnWXw3IIQrie+n4naxbdwbBdrLZGRkTRv3jzT/OXL4bbb4LXX/LOd\nYDmnLs+MgR49IDnZ6yQiIoUvIiLirAZkFRH/OXbsWPpl2TK68EI39Vcn55At6qKi3LRVK29ziIh4\noXLlymzbto2jR4/q/CwRD1hrSUpKYv/+/WzdujXLZejSCrl8XlLXp6A4/GqM+Q/QEahkjNkKjLHW\nvnGm5334IdxwA/z0E1xySUGnFBEJHGXKlAFg+/btJCUlFei2EhMTiYyMLNBtBCO1i2/h1C7FihUj\nMjKSGjVqZHrNmzfDG2/AI49A6mWJ/bM9/62q4Fhrb8zP866/Hm68EVq3hpQUd0hWRCRclClTJr24\nK0jz5s3Lcq6QqF2yo3aBtKuYjRnj3/WG7OHXNPHxbnr11d7mEBEREXnoITf980//rzvki7pSpWDi\nRPjyS1ih61GIiIiIR/bsgccfhzvvhOrV/b/+kC/qAO64w02bNvU2h4iIiISvypXd9OWXC2b9YVHU\nAezf76YDB3qbQ0RERMLPU0+56dq1BbeNsCnqypeHJ56At9+G9eu9TiMiIiLh4sAB+Oc/oV8/qFev\n4LYTNkUdwMiRblqQDSoiIiKSUdqwJe++W7DbCauiDmDXLje9/XZvc4iIiEjoSzt/rjA6a4ZdUVe5\nshvsb/Jk2LjR6zQiIiISquLjYdgwuOoqaNy44LcXdkUdnBrsr04db3OIiIhI6Cpb1k0/+6xwtheW\nRR3Ajh1ues893uYQERGR0DNxopv+9lvhXdEqbIu6qlXdqM4TJrhrsImIiIj4w6FDcNdd0KsXXHhh\n4W03bIs6gEcfddO0a7Dl1YnW4DUiAAAgAElEQVQ9Ozk0/7+c2LPTb5lERETEe2fzHV+unJvOmOHn\nUGdQrHA3F3h27IBzznFXnZg0KXfPSUlMZE2ftiTGnerKEhnbhAafLKRIZGQBJRUREZGCdrbf8RMm\nuOny5YV32DVNWO+pA3cYdtw4eOWV3A9KfPovGyAxbgVr+rQtgIQiIiJSWM7mO/7gQXeu/rXXenNp\n0rAv6gAefthNczMo8Yk9O7P8stMkxq3QoVgREZEgdbbf8eXLu+nHH/s7We6oqEu1d6+b9u+f83LH\n1iw/q8dFREQkMJ3Nd/wzz7jpihWFf9g1jYq6VBUrwrPPwnvvwapV2S8X1eCCHNdzpsdFREQkMOX3\nO373bnjgARgwoHAGGc6OiroM7r/fTRs3Bmt9L1M8piqRsU18PhYZ24TiMVULKJ2IiIgUpPx+x1ep\n4qZvvVVQyXJHRd1pDh5002uuyX6ZBp8szPJLT+sZIyIiIsErr9/x//ynmwbCpUfDfkiT05Ut64Y2\nueMO+OUXaNky6zJFIiNpNHMZJ/bs5Nia5UQ1uEB76EREREJAXr7jN2+Gp56Cv/89/2Pe+pOKOh9u\nv90VdRdfDMnJUCSb/ZnFY6qqmBMREQlBufmOTyvknn664PPkhg6/ZuPoUTdt3tzbHCIiIhJ4Bg1y\n0127PI2RiYq6bERFwaefuhGhZ83yOo2IiIgEipUrXaeIp56CypW9TnOKiroc9OoF9etDjx5w/LjX\naURERMRr1kKT1H4UDzzgbZbTqag7g7Qx69IuzisiIiLhq1MnN00bLSOQqKg7g6JFYdEiSEyE117z\nOo2IiIh45dtv4fvv4d133WgZgUZFXS60bg09e8LQoXDggNdpREREpLCdOAFdu8K550K/fl6n8U1F\nXS59/rmbVqjgbQ4REREpfGnf/5s3e5sjJyrqcskYWLfO3X7oIW+ziIiISOGZPBmOHIEFC6BYAI/w\nq6IuD+rWhREj4PHHYf16r9OIiIhIQdu3z12U4OqroW1br9PkTEVdHj3zjJvWq+e6NYuIiEjoqlTJ\nTT/7zNscuaGiLh/SujF36OBtDhERESk499zjphs2eJsjt1TU5UPZsvCf/8D8+TB7ttdpRERExN9W\nr4YJE2D0aKhd2+s0uaOiLp/++ld3CLZ7dzh2zOs0IiIi4i/WQqNG7vbYsd5myQsVdWdhzRo3LVnS\n2xwiIiLiP5dc4qaHDnmbI69U1J2FIkXgt9/c7f/7P2+ziIiIyNmbMQOWLIFPPoEyZbxOkzcq6s7S\nhRfCbbfBgw/Cli1epxEREZH8OnQIrr0WWrZ002Cjos4PXnnFTWvW9DaHiIiI5F+5cm7688/e5sgv\nFXV+snevm/bo4W0OERERybthw9x0/Xp3FalgpKLOTypWhClTYNYs+PZbr9OIiIhIbv36K7z8Mjz6\nKJx/vtdp8k9FnR8NGuQOwXbtCkePep1GREREziQ5GVq0cLeD/druKur8LG3U6VKlvM0hIiIiZ3be\neW4aCmPOqqjzsyJF3CjUAMOHe5tFREREsvfaa7BjB8ydC5GRXqc5eyrqCkCDBm4E6hdegGXLvE4j\nIiIip9u1C4YOhV69oGNHr9P4h4q6AjJ6tJtedBEkJXmbRURERDKrWtVNP/3U2xz+FDRFnTHmCmNM\nnDFmvTHmQa/z5EZiopuWL+9tDhERETnlqqvcdMcOb3P4W1AUdcaYosDLQHegEXCjMaaRt6nOrEQJ\nWLwYjhyBp5/2Oo2IiIh8+aX7mTr11N66UBEURR3QClhvrd1grT0BfAD08jhTrlxyCdx6K/zjH/DH\nH16nERERCV8HDri9dBdeCAMHep3G/4y11usMZ2SMuQ64wlo7JPX+zcAl1tphpy03FBgKEBMTc9G0\nadMKPWt2li5104su8i5DQkIC0dHR3gUIUGoX39QuvqldslKb+KZ28c3LdgmE7+LsdOrUaam1tuXZ\nrKOYv8IUMF8X7MhSjVprXwVeBYiNjbUdA6g7S4sWULYsNG4MK1Z4k2HevHkEUpsECrWLb2oX39Qu\nWalNfFO7+OZVu/Tu7TpFbN0K1aoV+uYLRbAcft0KnJfhfnVgu0dZ8qVMGZg5E1auhLfe8jqNiIhI\n+Jg1yxV0b7wRugUdBE9RtwSoZ4ypbYwpDvwV+NzjTHl25ZVwxRXucmLbtnmdRkREJPQdPAg9ekCT\nJjB4sNdpClZQHH611p40xgwD5gBFgTettSs9jpUvX30FxkD16pCS4m6LiIhIwUgbVmz5cm9zFIag\nKOoArLWzgFle5/CHAwfcm6xhQ1izxus0IiIioalvXzfdsiU8dqIEy+HXkFKuHPz3vxAXBxMmeJ1G\nREQk9MyZAx9/DK++Cuedd+blQ4GKOo906wYDBsA998C6dV6nERERCR3x8e4c9gYN3Fix4UJFnYfS\nesHWrw/Jyd5mERERCRVly7rpqlXe5ihsKuo8duSIm+r6sCIiImfv+uvddPPm8DiPLiMVdR4rWdJd\nH/bwYRgzxus0IiIiwWv6dPjoI3jzTahRw+s0hU9FXQC45BJ44AEYNw5+/dXrNCIiIsFn+3bo0wc6\ndIBbbvE6jTdU1AWIp56CIkXc5cSOH/c6jYiISPBISTl1pYh58zyN4ikVdQEkMdFNIyO9zSEiIhJM\nqld300OHvM3hNRV1ASQi4lRPnTvu8DaLiIhIMBg3DnbsgPnz3XXWw5mKugDTsCE88wy88kp470IW\nERE5kyVLXCfDBx6Adu28TuM9FXUBaMQIqFMHOnWCffu8TiMiIhJ4jhyBVq3cKBJPPeV1msCgoi5A\nrV/vppUqgbXeZhEREQk00dFuGh/vbY5AoqIuQBkDe/e627Vre5tFREQkkFx9tZtu3AhFi3qbJZCo\nqAtgFSu68+o2b4bRo71OIyIi4r0PPoAvvoC334ZatbxOE1hU1AW4Dh1g5Eh49FH48Uev04iIiHjn\nzz/hxhvh8svh5pu9ThN4VNQFgSeegHPOgfbt4cABr9OIiIgUvpSUU5f+mjPH2yyBSkVdkNi2zU0r\nVFDHCRERCT+lS7upOkZkT0VdkDAGdu92txs29DaLiIhIYRoyBI4ehV9+OVXcSVYq6oJITAx88w3E\nxcHjj3udRkREpOB99BG88QZMmAAXXeR1msCmoi7IdOniRs5+6CFYvNjrNCIiIgVn/Xq4/nro1g2G\nDfM6TeBTUReEnnrKDXfSpo0uXiwiIqEpMRHq1XO3//tfb7MECxV1QSrt/Lpy5dRxQkREQk9UlJse\nP+5tjmCioi5IFSkCO3e62xdc4G0WERERf7rsMjfdsAGKF/c2SzBRURfEqlSB2bNhxQp3jp2IiEiw\n+/e/Yf58mD5dl8nMKxV1Qe4vf3GXEHv8cZg50+s0IiIi+bd4MYwYAXfcAb17e50m+KioCwFjx0Kr\nVtCzp9tVLSIiEmz273cdAMuVg4kTvU4TnFTUhYiffnLT88+HY8e8zSIiIpIX1rpRHQD27fM2SzBT\nURdC0noIlSypHrEiIhI8Kld20717XUdAyR81XQgpXhy2bnW3Y2O9zSIiIpIbd9/tirmFC0/trZP8\nUVEXYqpVg+++g3Xr4L77vE4jIiKSvRkz4KWX3KD6bdp4nSb4qagLQZ06wZNPwvPPw8cfe51GREQk\nqxUr4Npr3XfWAw94nSY0qKgLUf/8p7tObN++sGaN12lERERO2b8fmjZ1t7/7ztssoURFXQj75hs3\nbdgQEhK8zSIiIpIm7dy55GRvc4QaFXUhLinJTUuX9jaHiIgIwNKlbnrokHq6+puaM8QVK3bqGrHL\nl3ubRUREwttFF7lpXByUKeNtllCkoi4MVKkCCxa4vXa33up1GhERCUf33QfLlkG9elC/vtdpQpOK\nujDRti3UqAGvvw5vvul1GhERCSfvvedGZHjySe2hK0gq6sJITAzceCP87W8wb57XaUREJBwsXQr9\n+8PVV7uRGaTgqKgLM++/D3XrunGB1q3zOo2IiISyXbugZUt3+crPPvM6TehTUReG0oq5+vXhwAFv\ns4iISGg6cQKqVnW3NaxW4VBRF6ZOnnTTChVODXsiIiLiLyVKuGlCAhjjbZZwUexsV2CMqQAMBSoB\nvwAfWWs1nGCAK1oU4uPdCavFi0NKiv7oRETEP+rWddONG6FUKW+zhBN/7KmbAVwExAMDgKXGmKp+\nWK8UsNKlYfNmd7tyZW+ziIhIaLj1VvjjD5g7F2rV8jpNePFHUVfJWtvXWjvOWnsl8Cjwih/WK4Wg\nRg1YvBj27oWePb1OIyIiwWzCBDd01oQJ0LGj12nCjz+Kun3GmOi0O9baT4A6flivFJJLLoEPP4SZ\nM2HUKK/TiIhIMJo+He65B26/HYYN8zpNePJHUZcIzDLGNAcwxtQFdvlhvaSur68xZqUxJsUY09Jf\n65XMrr8exo6FJ56Ad9/1Oo2IiASThQuhTx/o0gUmTfI6TfjyR1H3I3AQmGmMOQbEASuMMV2NMeX8\nsP4VwLXAD35Yl+Rg9Gi49lq4+WZ3WTEREZEzWbsWLr0Uzj0XvvnG6zThLVe9X40xUdbaY74es9aO\ny7BcNeBioCXwANAcOKtT8K21q1PXfTarkVz65BM47zxo1w42bIDatb1OJCIigWrXLoiNdbe3bfM2\ni+SiqDPGdAbmGGNuttZ+kNOy1tptwDbgUz/lEw/8+acb3qROHTh4EMqW9TqRiIgEmoSEU4MLJ2sg\ns4BgrLU5L2DMJ0CMtfayHJZpCcQCn1trD+c5hDHfAL6GQRllrf0sdZl5wN+ttb/ksJ6huDHziImJ\nuWjatGl5jRLSEhISiI6OPvOCqZYuddOLLiqgQAEir+0SLtQuvqldslKb+Bbq7ZL2HdGiRd7GOQ31\ndsmvTp06LbXWnlXfgdwcfr0UONMleFcCXwIVgRfzGsJa2zWvz8lmPa8CrwLExsbajupPncm8efPI\nS5s0bw7lUs+KDOXBifPaLuFC7eKb2iUrtYlvodou1kKR1DPy9++H8uXz9vxQbZdAkJuOEuWBDTkt\nkHq+3VuARjoLIWXLwpYt7nYRXVBORESAqCg33bw57wWdFKzcfFXvBarkYrkfcYdg/coY09sYsxVo\ng+thO8ff25DsnXcerFnjbuvcOhGR8NaiBRw/Dr/+6gavl8CSm6LuZ6BPLpY7RO6Kvzyx1s6w1la3\n1paw1lax1v7F39uQnMXGwi+/uGvFNmrkdRoREfHCX//qirk5c6BZM6/TiC+5KereAPoaY649w3J1\ncNd/lRB00UXw3XewerUbXFJERMLHqFHuykNTpsDll3udRrJzxqLOWvsl8B/gA2PMOGNM6dOXMcYU\nB4YDGrI2hHXqBDNmuOKuXz+v04iISGF49VV3taFx42DQIK/TSE5yNfgwMAg4BjwE3GeMmY670sMu\noDowEKhN6nAiErquuQbefBMGD4aYGHj+ea8TiYhIQZk5E267DQYOhIcf9jqNnEmuijprbTIw1Bjz\nITAS6EfmvXw7gL7W2p/9H1ECzS23wIEDMGKEK+xGjfI6kYiI+NvcudCzJ7RtC1Onep1GciO3e+oA\nsNZ+C3ybek3XC4CyuL11y6y1JwsgnwSo+++HvXvhoYegYkW4/XavE4mIiL/8+CN07uw6x+la4MEj\nT0VdGmvtQeAHP2eRIPPEE7B7N9xxhxur6IYbvE4kIiJn66efoH17d+3vlSu9TiN5ka+iTiTN66+7\nCzr/9a/u6hN/0YAzIiJBa+lSaN3aXdN1Q46XHZBApOsEyFn74gto2RKuuAIWL/Y6jYiI5Mfy5e6z\nvEwZ2LHD6zSSHyrqxC+WLIFq1aBNG1ixwus0IiKSF6tWwYUXQrFicOiQ12kkv1TUid9s3eqmTZvC\nxo3eZhERkdxZuxYaN3a3k5K8zSJnR0Wd+FVKipvWqQPr13ubRUREcrZhg7sUJJz6/JbgpaJO/MqY\nUx8M9erBmjXe5hEREd82b4bzz3e3U1Lc57cENxV14ncZC7uGDdUlXkQk0GzbBrVqudsq6EKHijop\nEBkLuyZN4Pffvc0jIiLOzp1Qvbq7nZysgi6UqKiTAmMMWAsREdCsmRv/SEREvLNnD5xzjrt98iQU\nURUQUvTrlAJ34gSULevGP/rpJ6/TiIiEp337oHJldzspCYoW9TaP+J+KOikUBw+6Ecpbt3bXFBQR\nkcKzfz9UquRuHz/uxqOT0KOiTgrNjh3uxNz27WHePK/TiIiEh23boGJFdzsxEYoX9zaPFBwVdVKo\nNm50PWI7dYKvv/Y6jYhIaFu79lSniKQkKFHC2zxSsFTUSaFbtQqaN4fLL4dZs7xOIyISmn799dTA\nwsnJOuQaDlTUiSeWLXPXie3RAz77zOs0IiKhZd48aNECoqPd8FLq5Roe9GsWzyxc6A7DXnMNfPyx\n12lERELDZ5+5z9bYWDh8WOPQhRMVdeKp776D7t2hb194/32v04iIBLcpU9w/yp066TKN4UhFnXhu\n1izo3Rv69YNJk7xOIyISnJ55BgYPhptucv8wS/hRUScBYfp0GDIE7rwTRozI/fNO7NnJofn/5cSe\nnQUXTkTEA3n5fPvnP+GBB+C+++C99wohnAQk9YWRgPHaa1CvnvtwWrMGZs7MftmUxETW9GlLYtyK\n9HmRsU1o8MlCikRGFkJaEZGCkdfPt1tugalT4YknYOTIQgwqAUd76iSg/OMfMG2aOySbNraSL6d/\n4AEkxq1gTZ+2BZxQRKRg5eXzrXt3V9BNnqyCTlTUSQDq2xcWLHCjoBvjxlfK6MSenVk+8NIkxq3Q\noVgRCVp5+Xy74AKYPduNHjB0aGEllECmok4CUtu2sG6du12sGBw9euqxY2uW5/jcMz0uIhKocvv5\nVr48/O9/8M030KdPYSSTYKCiTgJW3bqwZ4+7XaoU7Nrlbkc1uCDH553pcRGRQHWmz6/I2AswBg4e\nhJ9/hi5dCimYBAUVdRLQKlWCY8fc7apVYfVqKB5TlcjYJj6Xj4xtQvGYqoWYUETEf3L6fCterwmR\nVdzn2+rVcPHFhZlMgoGKOgl4kZGnzqtr1Mhd/qbBJwuzfPCl9Q4TEQlmvj7fitRuwkWz3efbrl3Q\noIEXySTQaUgTCQpFioC17rI3nTrBO+9E0n/mMk7s2cmxNcuJanCB9tCJSEgoEhlJowyfb7/uv4BO\n17jPt+PHoXhxjwNKwNKeOgkqcXHupOCbb4axY92hirLtL1dBJyIhp3hMVd75/XI6XVOVunXdP7Yq\n6CQnKuok6Hz8MYwaBY884i6HIyISigYOhLvvdlfbSRsNQCQnOvwqQemxx6BOHfjb32DlSvj9d68T\niYj4h7VQowZs3eoGFdYYdJJbKuokaA0e7D74unVzgxQnJbkx7UREgtXx465zGMD8+dCunbd5JLjo\n8KsEta5d3Z46gIgI2L7d2zwiIvm1c+epgm7zZhV0kncq6iToNWoEhw+729WqwbffeptHRCSvliyB\nc85xt48edUchRPJKRZ2EhOhoSEmBhg3d3rsxY7xOJCKSO++8A61aQZky7nMsKsrrRBKsVNRJyDAG\nVq2Cf/0Lxo2Dpk29TiQikrPhw2HAALj2Wjh0yH2OieSXijoJOY8/DnPmwIoV7gPy6FGvE4mIZNWi\nBbzwAjz1FHzyiddpJBSor6CEpMsvd8MBVK8OpUq56yTqsjoiEggy9nD96iu44gpv80jo0J46CVnV\nqsGJE+52w4bw3nve5hERWbXqVEEXF6eCTvxLRZ2EtIgIN5Bnjx7Qvz/ccovXiUQkXL34IjRu7G4n\nJkL9+t7mkdCjok7CwpdfwsSJMHWq+y85JcXrRCISLqyFCy+Ee+91/1haCyVKeJ1KQlHAF3XGmKeN\nMWuMMcuNMTOMMeW8ziTB6Y474Jdf3PksRYvC3r1eJxKRULdvHxQpAsuXw6xZ8OabXieSUBbwRR3w\nNdDEWnsBsBYY6XEeCWIXXQQHDrjbMTGwcKG3eUQkdMXHQ6VK7vbu3dC9u7d5JPQFfFFnrf2vtfZk\n6t3FQHUv80jwK1fOHX6NiYFLL3WX5hER8aehQ2HdOnfFm7TPG5GCZqy1XmfINWPMF8CH1tp3s3l8\nKDAUICYm5qJp06YVZryAl5CQQHR0tNcxAsqff0Lx4gls3x5Ns2Ya+DMjvV98U7tkpTY5xVpYtszd\nrlMngfLl1S6n0/vFt06dOi211rY8m3UERFFnjPkGqOrjoVHW2s9SlxkFtASutbkIHRsba+Pi4vwb\nNMjNmzePjh07eh0j4MycOY+ePTsC8Ouv0KyZt3kChd4vvqldslKbOGvWuOGTAP73P9i7V+3ii94v\nvhljzrqoC4jDr9bartbaJj5+0gq6gUBPoF9uCjqRvChVyg0vANC8OYwY4W0eEQk+EyeeKugSE6FJ\nE2/zSHgKiKIuJ8aYK4B/Aldba3XBJykQJUq4wyajR8O//63Li4lI7ljrOmDddRfcfLOGKxFvBXxR\nB7wElAa+Nsb8Zox5xetAErrGjoWVK93tUqVg7lxv84hI4Nq/3w1XsmwZfPEFvP2214kk3AX8tV+t\ntXW9ziDhpVEjOHkSataEzp3hhhvggw+8TiUigeTLL+Gqq9ztXbugcmVv84hAcOypEyl0RYvC1q3w\n8svw4YfucOz+/V6nEhGvpaS4S31ddZU7by4lRQWdBA4VdSI5uPNO2LLF3a5YET75xNs8IuKdJUvc\nP3yrVsHnn7serhoGSQKJijqRMzjvPPffeOvWcN110K6dOxlaRMLHtddCq1buHLqjR08dehUJJCrq\nRHLBGFi0CKZNgwUL3Af71q1epxKRgrZhg/v7nzHDDVuSnAxRUV6nEvFNRZ1IHvTtC3v3utvnnQev\nqC+2SMj6xz/g/PPd7b174Y47vM0jciYq6kTyqGJFd/j1uuvch3yxYnDkiNepRMRf9u1ze+eeftoV\ndta6v3uRQKeiTiSfPvrIHYpNToboaNdTVkSC28SJUKmSu71hA/zf/3mbRyQvVNSJnIW2bV0nih49\nYNgw99992uFZEQkex465v9+77nKdIqyF2rW9TiWSNyrqRM6SMW4g0lWr3P2YGHjoIW8ziUjuff45\nlCzpbv/8s4YukuClok7ETxo2dP/d3347PP64K/Y2bfI6lYhkJznZ/d326uUGEk5Ohosv9jqVSP6p\nqBPxs0mT4M8/3e3atWHQIE/jiIgP33zjOjmtWeOu2/q//7mhikSCmd7CIgWgenW31+6xx+Ctt9xe\nu+XLvU4lIgcOuOKtWzcoV86dS9ezp9epRPxDRZ1IARo16tQ1Yy+8ELp00dUoRLxy551QoYL7G/zx\nR1fgRUZ6nUrEf1TUiRSw8uXdl8jrr8N337m9BN9/73UqkfAxd67bWz5pEtx/v/t7vPRSr1OJ+F8x\nrwOIhIu//Q1uugnOOQc6doR69WDlSoiI8DqZSGg6dAgqV4YTJ6BUKdi+HcqU8TqVSMHRnjqRQhQV\nBQcPwmefwbp1ULy4O+dORPxr+HB3ztyJE27PeEKCCjoJfSrqRDxw9dWQlAQtWrjesca48bFE5OzM\nn+/+nl54Ae65xx1qvewyr1OJFA4dfhXxSLFisHSpG/6kRg245BI3f/t2d4hWRHIvPh7OPdddh7l4\ncdi9G8qW9TqVSOHSnjoRj5133qneeOC+mC6+GI4f9zaXSLD4+99dAXfkiOuMdPy4CjoJTyrqRALE\npZe64m7yZPjlFzfUQtrhIxHJauFCd6j12WfdlVyshU6dvE4l4h0VdSIBZuhQSEmBIUNgwgQ3BMo7\n73idSiRwbN7sirlLL3XT/fvdcCUi4U5FnUgAMgZeew2OHoVGjWDAADdvyRKvk4l4Z/t214O8Vi13\n/9tv3T9A5ct7GkskYKioEwlgUVFuLLstW9z9Vq1ccbdzp7e5RArT7t1QsSJUqwaJifDll+5Qa+fO\nXicTCSwq6kSCQFpnih9+cPfPOQdat3ZjcImEqv373Xu/ShV3++OP3d9Bjx5eJxMJTCrqRIJI+/bu\nS23SJPjpJyhRAu66yx2CEgkVhw5B/fpu79zWrfDee+5936eP18lEApuKOpEgdPvtrpAbPBgmToSi\nRaFnT3doSiRYJSRAs2buShDr1rnrJVvrLq8nImemok4kSBkDb7wByclw990wc6Y7B69ePdizx+t0\nIrl37Bi0bQulS8Pvv8NLL7li7m9/8zqZSHBRUScS5IoUgRdfdF+CL74I69e7i5gbA6tXe51OJHvH\nj0OXLlCyJCxaBE895d7Hd93ldTKR4KSiTiSE3H23+1L88kt3v1EjV9x99523uUQySkx0pwtERrr3\n5tix7n37wANeJxMJbirqREJQjx7uS/LXX939Ll1ccTdlire5JLz9+qu7LmtUlDtd4MEH3bmho0d7\nnUwkNKioEwlhzZq54m7rVjcMyuDBrrgbOVKXH5PCYS089ph737VoAUlJ8Pbbbv748W6+iPiHijqR\nMFCtmhuN//Bh6NgRnnzSnYvXp487r0nE37ZtgyZN3Pvs4YfdECVbtrhi7uabvU4nEppU1ImEkeho\nmDsXTp5015idPt2d13T++epUIf7x4Ydu71v16u5qKA8/7Hpox8W5gYRFpOCoqBMJQ0WLwuTJbq/J\nM8/Ahg2nOlX07+/26InkVkqK2+trDPz1r27eTz+599e4cW5vnYgUPP2piYS5ESPcl++2bdChgxu9\nv0wZWLoUnn9e595J9hYvdoXcr7+6vb7XXw9Hjrj3TKtWXqcTCT8q6kQEgHPPhXnzTl1jtkgRuO8+\nN42MhO+/9zqhBIKTJ11HG2OgTRs3r04d97758EM35pyIeENFnYhk0b49NG/uDqtNmOA6U3Ts6L7I\n27VzvWklfBw5An//u/v9R0S4jjYXXug631gL5ct7nVBEQEWdiOTAGBg2zH1xHz4MAwfCggXuhHdj\nYPhw9Z4NVXv3ut+3Ma6DzbPPQqVKMGuWez/89psbJkdEAoeKOhHJlehomDrVfaHHxUHjxvDCC+7Q\nbNrAxikpXqeUs7FpE7njVg0AAA4+SURBVHTv7n6fMTFuPLlGjU51etizxz0uIoFJRZ2I5Fn9+rBi\nhfui/+wzN2/wYNer1hi4/PJTV7OQwPb7725QYGOgdm2YPdsdao+Lc7/flSvV6UEkWKioE5GzcvXV\n7svfWvj6a3cu3tdfnyoUjIF773V7eSQwzJ176hB6s2auAL/+etixw/0e5851hbuIBBcVdSLiN127\nwrJlrjA4fhxefBFKlHDTypVdEVG2LLz6qutFKYVj9Wq47bZTRXbnzq6zy7BhcOjQqZ6rVat6nVRE\nzoaKOhEpEMWLw913Q2KiKxq2b4fbb4f4eFdgREScGhZj/nyv04aOxER3fmPTpqeKuEaNXCENMHbs\nqd/JhAluTEIRCQ0q6kSkUJxzDkyadOpQ7aJFcNllbgDbyy47VYBccYXrkHHokNeJg8OqVXDrrafa\nLyrKnd+4YgVccgl88AGcOHGq3UePdntPRST0qKgTEU+0bu0GNLbWHYp9803X43LOHLjlFihX7lSh\nUqQI3HQTfP65t0OonNizk6IrfuHEnp2ebD8xEd54w/U8Tmubxo3h9dfd4/fd5y75llbALV4MN9zg\n9oqKSOgL+KLOGPOoMWa5MeY3Y8x/jTHnep1JRPyraFFXyO3efaogOXwY/vMf6NXL3U+7nTaEijFu\n3LS774aFCwt2OJWUxERW9WjBijY1KPn8v1jRpgarerQgJTHR79tKTHRX9Hj0UXeOYrFimffCDRni\n9s61bQsffQRJSafa7N//dj1YRSQ8FfM6QC48ba19GMAYcw8wGrjd20giUtCio93F4dMuEJ9m1y74\n5BNX5P34I7z0kvvJqFYtaNLE9eDM+HPuua44yqs1fdqSGLcCgLSnJ8atYE2ftjSauSxP6zpyxBWh\n33/vfn788czPueQS6NLFHWatVStv2UUkfAR8UWetjc9wtxSgy4uLhLEqVeDOO91PRmvXuh6c//mP\n6+25aVPe1lupUtYisGJFOLp9J1XiVuCrFjwWt4Ibr9zJrsSqHD7s9i4mJJB+2+by06p9e+jQwf20\naQOlSuUtu4gIgLG5/dTxkDHmcWAAcAjoZK31OeKVMWYoMBQgJibmomnTphVeyCCQkJBAdHS01zEC\njtrFt1Bsl5QUd3jz+HE3zXg7Odn3c8pv+YULZ/7LZ1Fngd97PMGBGi19PrdIEfcTFQWlS7ufUqXy\nt7cwkIXie8Uf1C6+qV1869Sp01Jrre8Pk1wKiKLOGPMN4GuEpFHW2s8yLDcSiLTWjjnTOmNjY21c\nXJwfUwa/efPm0bFjR69jBBy1i29qF+fEnp2saFMj28ebLNpC8ZjwHuBN7xXf1C6+qV18M8acdVEX\nEIdfrbVdc7no+8BM4IxFnYiIPxSPqUpkbJP0c+oyioxtEvYFnYgEjmDo/Vovw92rgTVeZRGR8NTg\nk4VExjYBTp3UGxnbhAafLPQulIjIaQJiT90ZPGmMiQVSgM2o56uIFLIikZE0mrmME3t28tOH73LJ\nDf21h05EAk7AF3XW2j5eZxARAXcoNrlJSxV0IhKQAv7wq4iIiIicmYo6Efn/9u42Rq6yCuD4/wi2\nRI3SgkhRtG2CKAnykoYQTUQqFPRDixG1JMSiGAISY1SiRfxgTBQh0RqjCaJBEBNAaog1SEhpAb9Q\nkBhKebG0gFGgUBHBGLQgHD/cZ/GyvbszW3Zndp7+f8lk7zz3Jc9z9tw7Z+69MyNJqoBFnSRJUgUs\n6iRJkipgUSdJklQBizpJkqQKWNRJkiRVwKJOkiSpAhZ1kiRJFbCokyRJqoBFnSRJUgUs6iRJkipg\nUSdJklQBizpJkqQKWNRJkiRVwKJOkiSpAhZ1kiRJFbCokyRJqoBFnSRJUgUs6iRJkipgUSdJklQB\nizpJkqQKWNRJkiRVwKJOkiSpAhZ1kiRJFbCokyRJqoBFnSRJUgUs6iRJkipgUSdJklQBizpJkqQK\nWNRJkiRVwKJOkiSpAhZ1kiRJFbCokyRJqoBFnSRJUgUs6iRJkipgUSdJklQBizpJkqQKWNRJkiRV\nwKJOkiSpAhZ1kiRJFbCokyRJqoBFnSRJUgUs6iRJkipgUSdJklSBkSnqIuKCiMiIOHDYfZEkSZpt\nRqKoi4hDgZOBvwy7L5IkSbPRSBR1wBrgq0AOuyOSJEmz0awv6iJiOfB4Zm4edl8kSZJmq32H3QGA\niLgFOLhj1kXA14FlfW7nHOCc8nRXRNw3PT2sxoHA08PuxCxkXLoZl27GZXfGpJtx6WZcuh3+WjcQ\nmbP3imZEHAlsAJ4vTe8AngCOy8wne6x7d2YumeEujhRj0s24dDMu3YzL7oxJN+PSzbh0m464zIoz\ndRPJzC3AQWPPI+LPwJLMtMKXJElqmfX31EmSJKm3WX2mbrzMXDiFxS+fqX6MMGPSzbh0My7djMvu\njEk349LNuHR7zXGZ1ffUSZIkqT9efpUkSarAyBZ1EfGJiLg/Il6OiAk/LRIRp0bE1ojYHhGrW+2L\nIuLOiNgWEddFxJzB9HxmRcT8iFhfxrU+IuZ1LHNiRNzTevwnIk4r866MiEdb844e/CimXz9xKcu9\n1Br7ulb73pwvR0fEHWV/uzciPtWaV02+THSsaM2fW/7320suLGzNu7C0b42IUwbZ75nWR1y+HBEP\nlNzYEBHvas3r3J9q0EdczoqIv7XG/7nWvFVln9sWEasG2/OZ00dM1rTi8VBEPNuaV3OuXBEROyf6\nmrVo/LDE7d6IOLY1b2q5kpkj+QDeS/OdLrfRfCK2a5l9gIeBxcAcYDNwRJn3K2Blmb4MOG/YY5qm\nuFwKrC7Tq4FLeiw/H3gGeEN5fiVw+rDHMay4AP+aoH2vzRfg3cBhZfoQYAewf035MtmxorXM54HL\nyvRK4LoyfURZfi6wqGxnn2GPaYBxObF1/DhvLC7leef+NOqPPuNyFvCjjnXnA4+Uv/PK9Lxhj2kQ\nMRm3/BeAK2rPlTK2DwLHAvdNMP+jwE1AAMcDd+5prozsmbrMfDAzt/ZY7Dhge2Y+kpkvANcCKyIi\ngKXA2rLcVcBpM9fbgVpBMx7ob1ynAzdl5vM9lht1U43LK/b2fMnMhzJzW5l+AtgJvHVgPRyMzmPF\nuGXasVoLfLjkxgrg2szclZmPAtvL9mrQMy6ZeWvr+LGJ5vtEa9dPvkzkFGB9Zj6Tmf8A1gOnzlA/\nB2mqMTkDuGYgPRuyzPw9zcmTiawAfpGNTcD+EbGAPciVkS3q+vR24K+t54+VtgOAZzPzv+Paa/C2\nzNwBUP4e1GP5ley+Y327nAJeExFzZ6KTQ9BvXPaLiLsjYtPYJWnMl1dExHE078IfbjXXkC8THSs6\nlym58BxNbvSz7qia6tjOpjnjMKZrf6pBv3H5eNk31kbEoVNcd9T0Pa5yiX4RsLHVXGuu9GOi2E05\nV2b1V5rEJD8flpm/6WcTHW05SftImCwuU9zOAuBI4OZW84XAkzQv3JcDXwO+tWc9Haxpiss7M/OJ\niFgMbIyILcA/O5bbW/PlamBVZr5cmkc2X8bp55hQ5fGkh77HFhFnAkuAE1rNu+1Pmflw1/ojpp+4\n/Ba4JjN3RcS5NGd5l/a57iiayrhWAmsz86VWW6250o9pO7bM6qIuM096jZt4DDi09XzsZ8aepjm9\nuW95xz3WPhImi0tEPBURCzJzR3kR3jnJpj4J3JCZL7a2vaNM7oqInwMXTEunB2A64lIuL5KZj0TE\nbcAxwK/Zy/MlIt4M3Ah8o1weGNv2yObLOBMdK7qWeSwi9gXeQnNJpZ91R1VfY4uIk2jeJJyQmbvG\n2ifYn2p4oe4Zl8z8e+vpT4FLWut+aNy6t017DwdvKvvBSuD8dkPFudKPiWI35Vyp/fLrH4DDovnk\n4hyaRFqXzR2It9LcTwawCujnzN8oWEczHug9rt3uaSgv7GP3kZ0GdH5aZwT1jEtEzBu7fBgRBwIf\nAB7Y2/Ol7Ds30Nzzcf24ebXkS+exYtwy7VidDmwsubEOWBnNp2MXAYcBdw2o3zOtZ1wi4hjgJ8Dy\nzNzZau/cnwbW85nVT1wWtJ4uBx4s0zcDy0p85gHLePXVklHVzz5ERBxOc9P/Ha22mnOlH+uAT5dP\nwR4PPFfeME89V4b9qZA9fQAfo6lidwFPATeX9kOA37WW+yjwEE3Ff1GrfTHNgXc7cD0wd9hjmqa4\nHABsALaVv/NL+xLgZ63lFgKPA68bt/5GYAvNi/MvgTcNe0yDigvw/jL2zeXv2eZLApwJvAjc03oc\nXVu+dB0raC4lLy/T+5X//faSC4tb615U1tsKfGTYYxlwXG4px+Cx3FhX2ifcn2p49BGXi4H7y/hv\nBd7TWvezJY+2A58Z9lgGFZPy/JvAd8etV3uuXEPzrQEv0tQtZwPnAueW+QH8uMRtC61v9JhqrviL\nEpIkSRWo/fKrJEnSXsGiTpIkqQIWdZIkSRWwqJMkSaqARZ0kSVIFLOokSZIqYFEnSZJUAYs6SZKk\nCljUSdIEIuKNEfGniLgrIl7fal8WES9HxPmTrS9Jg+QvSkjSJMpvm24C1mTm6og4CLgXuCszlw+3\nd5L0fxZ1ktRDRHwJ+B7ND2pfABwJHJWZTw+1Y5LUYlEnST1ERAA3AkuBOcDJmblhuL2SpFfznjpJ\n6iGbd79XA3OBzRZ0kmYjizpJ6iEiDgZ+APwROCoivjjkLknSbizqJGkS5dLrVcALwMk0xd0lEfG+\noXZMksbxnjpJmkREfAW4FFiambdHxByaT8POBZZk5r+H2kFJKjxTJ0kTKF9n8h3g4sy8HSAzXwDO\nABYC3x9e7yTp1TxTJ0mSVAHP1EmSJFXAok6SJKkCFnWSJEkVsKiTJEmqgEWdJElSBSzqJEmSKmBR\nJ0mSVAGLOkmSpApY1EmSJFXgf1SdISr0TFrCAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# calculate the analytical surface pressure coefficient\n", "cp_analytical = 1.0 - 4 * (y_cylinder / R)**2\n", "\n", "# plot the 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(x_cylinder, cp_analytical,\n", " label='analytical',\n", " color='b', linestyle='-', linewidth=1, zorder=1)\n", "pyplot.scatter([p.xc for p in panels], [p.cp for p in panels],\n", " label='source-panel method',\n", " color='#CD2305', s=40, zorder=2)\n", "pyplot.title('Number of panels : %d' % N_panels, fontsize=16)\n", "pyplot.legend(loc='best', prop={'size':16})\n", "pyplot.xlim(-1.0, 1.0)\n", "pyplot.ylim(-4.0, 2.0);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Challenge task" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have computed the pressure coefficient on the surface of the cylinder, it will be interesting to visualize what the streamlines look like.\n", "\n", "To do that, we use the function `streamplot()` from Matplotlib, requiring the Cartesian velocity components (`u`,`v`) on a mesh grid (`X`,`Y`). Therefore, the first step is to derive the equations for the velocity components.\n", "\n", "The potential at point $\\left(x,y\\right)$ of the $N_p$ source sheets in a uniform horizontal flow $U_\\infty$ is\n", "\n", "$$\\phi\\left(x,y\\right) = U_\\infty x + \\sum_{j=1}^{N_p} \\frac{\\sigma_j}{4\\pi} \\int \\ln \\left( \\left(x-x_j(s_j)\\right)^2 + \\left(y-y_j(s_j)\\right)^2 \\right) {\\rm d}s_j$$\n", "\n", "And the velocity field at point $\\left(x,y\\right)$ is\n", "\n", "$$u\\left(x,y\\right) = \\frac{\\partial}{\\partial x}\\left\\lbrace \\phi\\left(x,y\\right) \\right\\rbrace$$\n", "\n", "$$v\\left(x,y\\right) = \\frac{\\partial}{\\partial y}\\left\\lbrace \\phi\\left(x,y\\right) \\right\\rbrace$$\n", "\n", "Your task will be to:\n", "\n", "* derive the Cartesian velocity components\n", "* create a mesh grid\n", "* compute the velocity field on the mesh grid\n", "* plot the results\n", "* change the number of panels to improve the visualization" ] }, { "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": 13, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 13, "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 }