{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Tài liệu này mang giấy phép Creative Commons Attribution (CC BY).\n", "(c) Nguyễn Ngọc Sáng, Zhukovsky 12/2018.\n", "\n", "[@SangVn](https://github.com/SangVn) [@VnCFD](https://vncfdgroup.wordpress.com/)\n", "\n", "*Thực hành CFD với Python!*\n", "\n", "# Bài 13. Định lý Godunov, sơ đồ TVD, sơ đồ MUSCL" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Xét phương trình Burgers không nhớt:\n", "$$\\frac {\\partial u}{\\partial t} + \\frac{\\partial F}{\\partial x} = 0; \\ F=\\frac{u^2}{2} \\qquad (1)$$\n", "\n", "với sơ đồ sai phân phương pháp thể tích hữu hạn:\n", "$$\\frac{u_i^{n+1} - u_i^n}{\\Delta t} + \\frac{F_{i+1/2} - F_{i-1/2}}{\\Delta x} = 0 \\qquad(2)$$\n", "\n", "## 1. Định lý Godunov\n", "\n", "Định lý Godunov phát biểu về tính đơn điệu của sơ đồ sai phân.\n", "Sơ đồ sai phân gọi là đơn điệu (monotone scheme) nếu như với hai điều kiện ban đầu $u^0 > v^0$ thì $u^n > v^n$.\n", "Nếu sơ đồ sai phân cho nghiệm $u^n$ đơn điệu tăng (giảm) và $u^{n+1}$ cũng đơn điệu tương tự thì có nghĩa là tính đơn điệu được duy trì (monotonicity preserving).\n", "\n", "Godunov đã chỉ ra rằng: `sơ đồ sai phân đơn điệu có bậc không quá bậc một`.\n", "\n", "Thực vậy, bài 2 và bài 12 cho thấy, sơ đồ bậc hai làm xuất hiện các dao động khi nghiệm gián đoạn, tức là điều kiện đơn điệu không được thỏa mãn, sơ đồ cho nghiệm không chính xác. Tuy nhiên so với sơ đồ bậc một, sơ đồ bậc hai cho nghiệm chính xác hơn nếu nghiệm trơn (mượt), tại nơi nghiệm gián đoạn sơ đồ bậc hai tuy làm xuất hiện các dao động, nhưng độ dốc gián đoạn lớn hơn. Sẽ là một sơ đồ tốt nếu kết hợp được hai sơ đồ trên. Ý tưởng là ta sẽ sử dụng một 'công tắc chuyển đổi' giữa sơ đồ bậc 2 và bậc 1 tùy vào nghiệm mượt hay gián đoạn.\n", "\n", "## 2. Sơ đồ TVD (total variation dimishing)\n", "Lax (1973) đã chỉ ra rằng, định luật bảo toàn một đại lượng vô hướng được mô tả bằng phường trình Burgers, `có biến phân toàn phần` (total variation) của nghiệm vật lý không tăng theo thời gian. Với biến phân toàn phần (TV) là:\n", "$$TV = \\int |\\partial u/\\partial x| dx$$\n", "Đối với nghiệm rời rạc:\n", "$$TV(u) = \\sum_i|u_{i+1}-u_i|$$\n", "Một sơ đồ sai phân được gọi là TVD nếu như:\n", "$$TV(u^{n+1}) \\leq TV(u^n)$$\n", "\n", "Harten (1983) chứng minh:\n", "1. sơ đồ đơn điệu thì TVD\n", "2. sơ đồ TVD thì duy trì tính đơn điều (monotonicity-preserving: nếu $u_{i+1}^n > u_i^n$ thì $u_{i+1}^{n+1} > u_i^{n+1} \\forall i$ )\n", "\n", "Chúng ta có thể xây dựng sơ đồ TVD bằng phương pháp MUSCL.\n", "\n", "## 3. Sơ đồ MUSCL (Monotonic Upwind Scheme for Conservation Laws)\n", "\n", "Quay trở lại với phương pháp Godunov ở bài 12, hàm dòng F trên bề mặt thể tích được tính như sau: \n", "$$F_{i+1/2} = F(u_{i+1/2}^*(u_{i+1/2}^L, u_{i+1/2}^R))$$\n", "với $u^*$ là nghiệm bài toán Riemann $u^*(u_L, u_R)$. Phương pháp Godunov có thể chia là ba bước:\n", "* **Bước 1:** Tái cấu trúc nghiệm (reconstruction) trong từng thể tích hữu hạn. Thực chất là xác định giá trị uL, uR ở hai bên bề mặt để tìm $u^*$ ở trên bề mặt thể tích.\n", "* **Bước 2:** Giải bài toán Riemann về phân rã gián đoạn trên bề mặt thể tích $u^*(u_L, u_R)$ tìm $u^*$, tính dòng F.\n", "* **Bước 3:** Tìm nghiệm ở bước thời gian tiếp theo $u^{n+1}$\n", "\n", "Ở bước 1, để tìm $u_{i+1/2}^*(u_{i+1/2}^L, u_{i+1/2}^R)$ Godunov đề xuất cho các hằng số: $u_{i+1/2}^L = u_i, u_{i+1/2}^R = u_{i+1}$. Cách này gọi là `tái cấu trúc hằng số từng mảnh` (piecewise constant). Ngoài ra còn có `tái cấu trúc tuyến tính từng mảnh` (piecewise linear) và `tái cấu trúc parabolic từng mảnh` (piecewise parabolic) (thứ tự tương ứng như trên hình):\n", "\n", "\n", "Ví dụ tái cấu trúc tuyến tính từng mảnh, sử dụng khai triển Taylor:\n", "$$\n", "u_{i+1/2} = u_i + \\frac{\\partial u}{\\partial x} (x_{i+1/2} - x_i)\\\\\n", "u_{i+1/2} = u_i + \\frac{u_i - u_{i-1}}{\\Delta x} \\frac{\\Delta x}{2}\\\\\n", "u_{i+1/2} = u_i + 0.5(u_i - u_{i-1}) \\qquad (3)\n", "$$\n", "Công thức (3) nên sử dụng khi nghiệm mượt mà, khi nghiệm gián đoạn ta chỉ cần lấy $u_{i+1/2} = u_i$. Ta có thể thêm vào (3) `công tắc chuyển đổi` $\\phi: \\ u_{i+1/2} = u_i + 0.5\\phi(u_i - u_{i-1})$, với $\\phi$ bằng 0 khi nghiệm gián đoán, bằng 1 khi nghiệm trơn.\n", "\n", "Trên ý tưởng đó, ra đời sơ đồ **MUSCL** dùng `tái cấu trúc tuyến tính từng mảnh` và sử dụng 'công tắc chuyển đổi' gọi là `giới hạn độ dốc` (slope limiter):\n", "\n", "$$\n", "u_{i+1/2}^L = u_i + 0.5 \\phi(r_i)(u_{i+1}-u_i);\\ u_{i+1/2}^R = u_{i+1} - 0.5 \\phi(r_{i+1})(u_{i+2}-u_{i+1})\\\\\n", "u_{i-1/2}^L = u_{i-1} + 0.5 \\phi(r_{i-1})(u_i-u_{i-1});\\ u_{i-1/2}^R = u_i - 0.5 \\phi(r_i)(u_{i+1}-u_i)\\\\\n", "r_i = \\frac{u_i - u_{i-1}}{u_{i+1} - u_{i}} \\qquad (4)\n", "$$\n", "\n", "Hàm $\\phi(r)$ chính là giới hạn độ dốc. Có nhiều hàm giới hạn khác nhau [tham khảo](https://en.wikipedia.org/wiki/Flux_limiter):\n", "\n", "**minmod** - (Roe, 1986):\n", "\n", "$$\\phi_{mm}(r) = max[0, min(1, r)]$$\n", "\n", "**superbee** - (Roe, 1986):\n", "\n", "$$\\phi_{sb}(r) = max[0, min(2r,1), min(r, 2)]$$\n", "\n", "**van leer** - (van leer 1974):\n", "\n", "$$\\phi_{vl}(r) = \\frac{r+|r|}{1+|r|}$$\n", "\n", "**Koren** - (Koren 1993)\n", "\n", "$$\\phi_{kr} = max\\left[0, min\\left(2r, min\\left(\\frac{1+2r}{3}, 2\\right)\\right)\\right]$$\n", "\n", "\n", "\n", "Hàm giới hạn bằng không khi r<0 và bằng 1 khi r=1, điều này đảm bảo rằng tại những nơi nghiệm gián đoạn, sơ đồ sẽ chuyển về bậc một, nhưng về tổng thể ta có sơ đồ bậc hai.\n", "\n", "## Bài toán\n", "Xét bài lan truyền sóng xung kích với điều kiện ban đầu như ở bài 12: $u[0." ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "#chia lưới trên đoạn [0, 1]\n", "nx = 41\n", "x, dx = np.linspace(0, 1, nx, retstep=True)\n", "\n", "#bước thời gian\n", "r = 0.4\n", "dt = r*dx\n", "\n", "#điều kiện ban đầu\n", "uL = 1.2\n", "uR = 0.4\n", "\n", "u0 = np.zeros(nx)\n", "i = int(0.5/dx)\n", "u0[:i] = uL\n", "u0[i:] = uR\n", "\n", "#Phân rã gián đoạn godunov trên từng mặt thể tích hữu hạn \n", "def decay_godunov(uL, uR):\n", " ustar = 0\n", " if uL > uR:\n", " if (uL + uR)/2. > 0: ustar = uL\n", " else: ustar = uR\n", " else:\n", " if uL > 0: ustar = uL\n", " elif uR < 0: ustar = uR\n", " else: ustar = 0\n", " return ustar\n", "\n", "\n", "#tính hàm dòng F(u*)\n", "def flux(u):\n", " return u**2/2.\n", "\n", "#Ta có thể viết code sơ đồ muscl như sau đây, nhưng cách này chậm, tốn bộ nhớ\n", "def muscl(limiter, nt):\n", " u = u0.copy()\n", " un = np.zeros(nx)\n", " for n in range(nt):\n", " un = u.copy()\n", " for i in range(2,nx-2):\n", " r_im1 = (un[i-1] - un[i-2])/(un[i] - un[i-1] + 1e-6)\n", " r_i = (un[i] - un[i-1])/(un[i+1] - un[i] + 1e-6)\n", " r_ip1 = (un[i+1] - un[i])/(un[i+2] - un[i+1] + 1e-6)\n", " \n", " u_im05L = un[i-1] + 0.5*limiter(r_im1)*(un[i] - un[i-1])\n", " u_im05R = un[i] - 0.5*limiter(r_i)*(un[i+1] - un[i])\n", " u_im05 = decay_godunov(u_im05L, u_im05R) #u_{i-1/2}\n", " \n", " u_ip05L = un[i] + 0.5*limiter(r_i)*(un[i+1] - un[i])\n", " u_ip05R = un[i+1] - 0.5*limiter(r_ip1)*(un[i+2] - un[i+1])\n", " u_ip05 = decay_godunov(u_ip05L, u_ip05R) #u_{i+1/2}\n", " \n", " u[i] = un[i] - dt/dx*(flux(u_ip05) - flux(u_im05))\n", " \n", " return u\n", "\n", "# viết lại đoạn code như sau \n", "def MUSCL_scheme(limiter, nt):\n", " u = u0.copy()\n", " r = np.zeros(nx)\n", " uL = np.zeros(nx)\n", " uR = np.zeros(nx)\n", " for n in range(nt):\n", " #tìm r, để tránh chia cho 0 ta cộng vào mẫu số 1e-6\n", " r[1:-1] = (u[1:-1] - u[:-2])/((u[2:] - u[1:-1]) + 1e-6)\n", " \n", " #tái cấu trúc nghiệm u_{i+1/2}^L\n", " uL[1:-2] = u[1:-2] + 0.5*limiter(r[1:-2])*(u[2:-1] - u[1:-2])\n", " \n", " #tái cấu trúc nghiệm u_{i+1/2}^R\n", " uR[1:-2] = u[2:-1] - 0.5*limiter(r[2:-1])*(u[3:] - u[2:-1])\n", " \n", " #giải bài toán riemann \n", " ustar = np.array([decay_godunov(ul, ur) for ul, ur in zip(uL[1:-2], uR[1:-2])])\n", "\n", " #tính hàm dòng và tìm nghiệm u ở bước thời gian tiếp theo \n", " u[2:-2] = u[2:-2] - dt/dx*(flux(ustar[1:]) - flux(ustar[:-1]))\n", " \n", " #điều kiện biên không đổi \n", " return u\n", "\n", "\n", "#slope limiter: phi(r), rs - mảng r\n", "def minmod(rs):\n", " phi = np.array([max(0, min(1,r)) for r in rs])\n", " return phi\n", "\n", "def koren(rs):\n", " phi = np.array([max(0, min(2*r, min((1+2*r)/3, 2))) for r in rs])\n", " return phi\n", "\n", "def superbee(rs):\n", " phi = np.array([max(0, min(2*r, 1), min(r,2)) for r in rs])\n", " return phi\n", "\n", "def vanleer(rs):\n", " phi = np.array([(r+abs(r))/(1+abs(r)) for r in rs])\n", " return phi\n", "\n", "#tìm nghiệm\n", "nt = 25\n", "u_mm = MUSCL_scheme(minmod, nt)\n", "u_kr = MUSCL_scheme(koren, nt)\n", "u_sp = MUSCL_scheme(superbee, nt)\n", "u_vl = MUSCL_scheme(vanleer, nt)\n", "\n", "#nghiệm chính xác: vận tốc lan truyền sóng xung kích s=(uL+uR)/2. = 0.8\n", "u_exact = np.zeros(nx)\n", "i = int((0.5 + 0.8*nt*dt)/dx)\n", "u_exact[:i] = uL\n", "u_exact[i:] = uR\n", "\n", "#tải nghiệm sơ đồ godunov ở bài 12\n", "u_godunov = np.genfromtxt('data/burgers_godunov_test1.dat')\n", "\n", "#lưu giá trị u_sp vào file để so sánh ở bài 14\n", "np.savetxt('data/burgers_muscl_test1.dat', u_sp)\n", "\n", "\n", "fig = plt.figure(figsize=[15, 5])\n", "i = int(0.6/dx)\n", "j = int(0.8/dx)\n", "\n", "plt.plot(x[i:j], u_exact[i:j], x[i:j], u_mm[i:j], x[i:j], u_kr[i:j])\n", "plt.plot(x[i:j], u_sp[i:j], x[i:j], u_vl[i:j], x[i:j], u_godunov[i:j])\n", "plt.legend(['exact','minmod', 'koren', 'superbee', 'vanleer', 'godunov'])\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Kết quả cho thấy, rõ ràng sơ đồ bậc hai MUSCL cho nghiệm chính xác hơn sơ đồ Godunov bậc một. Giới hạn độ dốc superbee cho kết quả tốt nhất, minmod kém nhất." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3YAAAEyCAYAAAC2+0LeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzs3XdcVfUfx/HX4XLZKAKiDBmKA8GNK7fmStPMkbkbWplty7L6qWU2zMqZaZmVZllWNswciCjiALeCLJGlTEU2d5zfH9cs3OMiYJ/n43Ef3HvO93y/33NRHvd9z/d8v4qqqgghhBBCCCGEqL4sKrsDQgghhBBCCCFujwQ7IYQQQgghhKjmJNgJIYQQQgghRDUnwU4IIYQQQgghqjkJdkIIIYQQQghRzUmwE0IIIYQQQohqToKdEEIIIYQQQlRzEuyEEEIIIYQQopqTYCeEEEIIIYQQ1ZxlZXfgalxdXVVfX9/K7oYQQgghhBBCVIqoqKhsVVVr30jZKhvsfH19iYyMrOxuCCGEEEIIIUSlUBTl1I2WlaGYQgghhBBCCFHNSbATQgghhBBCiGpOgp0QQgghhBBCVHNV9h67K9HpdKSmplJSUlLZXRHXYGNjg5eXF1qttrK7IoQQQgghxH9CtQp2qampODo64uvri6Iold0dcQWqqpKTk0Nqaip+fn6V3R0hhBBCCCH+E6rVUMySkhJcXFwk1FVhiqLg4uIiV1WFEEIIIYS4g6pVsAMk1FUD8jsSQgghhBDizjJLsFMUZYWiKJmKohy9yv7RiqIcvvDYpShKC3O0K4QQQgghhBDCfFfsVgL9rrH/JNBNVdXmwNvAMjO1e1cJDQ1l165dt11PUlISQUFB1y03Z86ccq/vueee225bCCGEEEIIceeZZfIUVVXDFEXxvcb+f6eV3YCXOdq924SGhuLg4HDHAtacOXOYPn36xdfmCJVCCCGEqL50GRkUhIWBqlZK+6fOJ5NTnH1H2lJVMBi0pofe6uJz1Mq5paTUYCS/RA+V89abgXrxoaj/vFb+tU9R1QvnZ/zXdlAwXuW8y29Urrj9kgPVq764+FK5ypvc9YmHqN+0+RX3VQeVMSvmY8CfV9qhKMokYBKAt7f3nezTDVu1ahULFiygrKyM9u3bs2TJElJTU7n33nuJiIjA2dmZbt268eabb9KnTx8eeOABUlJSKCkp4bnnnmPSpEkAbNy4kenTp2MwGHB1deWLL75g6dKlaDQaVq1axcKFC+nSpcvFdrdv385zzz0HmO5hCwsLw8HBgVdeeYU///wTRVF44403eOihh8r1d+XKlURGRrJo0SIABg4cyNSpU9m4cSPFxcW0bNmSwMBAVq9ejYODAwUFBaiqesV6Q0NDmTlzJq6urhw9epQ2bdqwatUquadOCCGEuAuUJSdzauw49BkZldYHa8DjFo5TUdBp7S88HEw/Lf9+7kCZ1h691p6yv/dp7dFb2oFS7aabEBUo6fAhCXY3SlGUHpiCXecr7VdVdRkXhmkGBwdf8/uKWb8d43j6ebP2r6lHDWbcH3jV/dHR0Xz//feEh4ej1WqZPHkyq1evZty4cUybNo0nn3yS9u3b07RpU/r06QPAihUrcHZ2pri4mLZt2zJ06FCMRiMTJ04kLCwMPz8/cnNzcXZ25sknn8TBwYGpU6de1vaHH37I4sWL6dSpEwUFBdjY2PDTTz9x8OBBDh06RHZ2Nm3btqVr1643dK7vvfceixYt4uDBg5ftu1a9Bw4c4NixY3h4eNCpUyfCw8Pp3PmKv04hhBBCVBNlKSmcGj8BtbQUn29Xo/X0vGNtH848xAf7PiCjOJORjR/ifr/70ZVAWbGR0mIjZUXqP8+LVUqLjJSVGCkrMr0uKzZSVnL1j40WlmBta4GVrQXWtgo1bC2wsrPAyla5uN30ULCytcDiDmW99Lxifjt0mj2JOSiKQo/GtenlXxNbylAMxSi6YtAXo+iLUXRFWOhMz9EVY2Ew/VT0RShlpSiGIhRdEcqFMoqqv277qmKBammHqrVFtbQ1PddYgoUWFAUUDaqiMYVfiws/FQ1YaFAvPFf/tc300wIUC1TFEhQF9e/tioWpLgsNYHFxu2phcXG/6Y3/+81XTH3428XnCv9ct1P+VU759+W8C89N29Vy+/9V/t8HXHjq3+b+G/8FVkF3LNgpitIc+Bzor6pqzp1q15y2bt1KVFQUbdu2BaC4uBg3NzcAHn/8cX744QeWLl1aLiwtWLCAn3/+GYCUlBTi4uLIysqia9euF9d5c3Z2vm7bnTp14sUXX2T06NE8+OCDeHl5sXPnTh5++GE0Gg116tShW7du7Nu3j+bNb++bhqvVW6NGDdq1a4eXl2kkbcuWLUlKSpJgJ4QQQlRjZalpnBo/HmNRET4rv8QmIOCOtFtqKGXRgUV8dewrvJy8mNNhBWnrVX5Zl33V4YgaSwtsHbVY21thV0OLs4cWW3st1g5abB202DhosbHXYutgZXruoEVrpbkj53OjYlIyWLd5B2kJh2moyWC+33maWmWiPZsAYWevX4GlDVg7mh52jlCrBli7/rPt4qPGtbdZ2pQPT6LauyPBTlEUb+AnYKyqqrHmqPNaV9YqiqqqjB8/nnffffeyfUVFRaSmpgJQUFCAo6MjoaGhbNmyhYiICOzs7OjevTslJSWoqnrTwxdfffVVBgwYwIYNG+jQoQNbtmxBvYHx75aWlhiNxouvb2R9uWvVa21tffG5RqNBr7/+N0JCCCGEqJp06ekkjx+PsaAQ7y9X3LFQF50TzfSd04k/F8+IRiN4UDOB8GWJaDQWtOnng10NK2zsTcHsYkiz12JpZVE9bgExGuBcMuQkQE485MSTnxZNWWYsTfSZvA6gvVC2zB0c/aHpYHDyAZual4cymxqmbVYOYGlViScmqjKzBDtFUdYA3QFXRVFSgRlc+OeqqupS4H+AC7Dkwn9GvaqqweZo+07q1asXgwcP5oUXXsDNzY3c3Fzy8/Px8fFh2rRpjB49Gh8fHyZOnMjvv/9OXl4etWrVws7OjpiYGHbv3g1Ax44defrppzl58mS5oZiOjo6cP3/l4aUJCQk0a9aMZs2aERERQUxMDF27duWzzz5j/Pjx5ObmEhYWxty5c8uFN19fX5YsWYLRaCQtLY29e/de3KfVatHpdGi12nJtXa3emJiYCnhXhRBCCFEZdKdPc2r8BAznz+O9YgW2gRX/pbneqGfF0RV8evBTatnUYlH3xVjsqUtoSBx1/GrQd2IQjs42Fd4Ps1BVKMy+GNzKPXITwVB2sWiRYkeCoS6pFo1w8X6AFi2CsfNoAs71TcFNCDMw16yYD19n/+PA4+ZoqzI1bdqU2bNn06dPH4xGI1qtlsWLF5OUlMS+ffsIDw9Ho9Gwbt06vvzyS0aNGsXSpUtp3rw5jRs3pkOHDgDUrl2bZcuW8eCDD2I0GnFzc2Pz5s3cf//9DBs2jPXr1182econn3zCtm3b0Gg0NG3alP79+2NlZUVERAQtWrRAURQ++OAD6tatS1JS0sXjOnXqhJ+fH82aNSMoKIjWrVtf3Ddp0iSaN29O69atWb169cXtQ4YMuWK9EuyEEEKIu4MuI8MU6s6exXvFF9g2u/4ySbcrKS+J13e+zuHsw/T37c/zjV5m19enyDiZSvOeXtzzoD8ayyo4mUlZYbkrb+UeJXn/lLPQmoKaiz9qw74kGOuyKk7L72n2GO1cebxnfcZ28MHRRnv1toS4DcqNDOerDMHBwWpkZGS5bdHR0QTcoSEC4vbI70oIIYSomnSZmSSPHYc+OxvvLz7HtmXLCm1PVVW+O/EdH0V+hJXGijc6vEFgUTs2rziOwWCk59gA/Nu4VWgfbkhZESTvgqwT/wS37HjITy9frmY9cGkALv4XHg1Nr2vWQ7XQsD02i0Uh8USeOkttR2ue6FqfUe29sbOqjMnoRXWnKErUjY50lH9hQgghhBD/EfqsLJLHT0CflUW9zys+1J0pPMP/wv9HxOkIOnl2YlaHWSRtK+S3DYdw8bCn36RmONWxq9A+XFNRLsRuhJg/IH4r6ItN222cwLUh1O9ePsQ51wery/urqipbojNZFBLHodQ8PGra8NbgQEYE18NGW7UmbxF3Lwl2QgghhBD/AfrsbE5NeARdRgbeyz7DrnWrCmtLVVU2nNzAO3veQW/U82aHNxlQdzBbPj9OasxZmnSsS9eHG1fOjJVnT8GJDaYwd2oXqAZw9IBWY6Bxf/BoBXbXn7EcwGhU+fPoGRaGxBFzJp96zra8+2Azhrb2wqoqDisVdzUJdkIIIYQQdzl9Tg7JjzyCLj2dep8txS644uawO1dyjrd3v82mU5toUbsFczrPQZtZgx/m7KOkSE+PsU1o2ulWliG/RaoKGUdNQS7mdzhzxLS9dgB0fgGaDDCFuZuYbVNvMPL74dMs2hZPfGYB9V3tmTe8BYNbemCpkUAnKocEOyGEEEKIu5j+7FmSH3mUspRU6i1din27dhXWVlhqGDN2zeBc6Tmea/0cE5pO4EhIOhE/H8DRxYahU1pQu94dmAXSoIeU3f+EuXPJgAL12kPvt01hzqXBTVerMxj5eX8aS0LjScoponEdRxY+3Ir7mrmjsagGyzCIu5oEOyGEEEKIu9TFUHfqFPWWfop9h/YV0k6hrpC5++ayLm4dDWs1ZOm9S/G1bcDm5dEkHsyifsva9BwfgLVtBX70LCuCxG2mMHfiTyjOBY216T65LlNNwywdbm2SllK9gR8iU/k0NIG0c8UEedZg6Zg29GlaBwsJdKKKkGAnhBBCCHEXMpw7R/Jjj1GWmIjXkiXYd+xYIe1EZUTx+s7XSS9I55GgR5jScgp56aWs/XgfBTkldBrmT4te9SpmYfErTX5iXRMa9TVdlfPvdVvrxBWXGVizN5nPwhLIOF9KK28nZj8QRPfGtavHQuniP0WCXQX49ddfOX78OK+++uoda9PBwYGCgoI71p4QQgghqi5DXh7Jjz1OWVw8XksW49C5k9nbKDOUsejAIlYeW4mHgwcr+62klVsronedJmxNLDYOWh54sRXu/k7mbfhak580GQC+nUFze2vFFZbqWbX7FMt3JJJdUEZ7P2c+GtGSexq4SKATVZYEuwowaNAgBg0aVNndEEIIIcR/kOH8eZIfn0hJbCxeCxfg0KWL2duIyY3htR2vEX8unmGNhjE1eCpWqg0hX0UTs/sM9QJq0fvRQGwdrW6/sQqY/ORqzpfo+Co8iS/CT3KuSEeXhq4807Mh7fxubJZMISqTBLublJSURL9+/ejcuTO7d++mRYsWPPLII8yYMYPMzExWr17N8ePHiYyMZNGiRUyYMIEaNWoQGRnJmTNn+OCDDxg2bBihoaHMmDGDOnXqcPDgQR588EGaNWvG/PnzKS4u5pdffqFBgwacOnWKRx99lKysLGrXrs2XX36Jt7c3J0+eZNSoUej1evr161fZb4sQQgghqgBDQQHJEydSEhOD1/z5OHbvbtb69UY9K4+tZPHBxThZO7G412K6enXlXEYRvy2LJCe9kLYDfAke4Hd7955dbfIT7w7QZzY0vu+WJj+5mrOFZXwZfpIvdyWRX6KnVxM3pvT0p5V3LbO1IURFq77B7s9X//nGxlzqNoP+7123WHx8PD/88APLli2jbdu2fPvtt+zcuZNff/2VOXPm8MADD5Qrf/r0aXbu3ElMTAyDBg1i2LBhABw6dIjo6GicnZ2pX78+jz/+OHv37mX+/PksXLiQTz75hClTpjBu3DjGjx/PihUrePbZZ/nll1947rnneOqppxg3bhyLFy827/sghBBCiGrHUFBIyuMTKTl2HK/5n+DYs4dZ608+n8z0ndM5lHWIPj59eLPDmzjZOBEflUnIN9FoNBbcP6UF3oEut9aAqkJCCBxdZ/bJT64mu6CU5TsSWRVxisIyA/0C6zKlpz9BnjXN2o4Qd0L1DXaVyM/Pj2bNmgEQGBhIr169UBSFZs2akZSUdFn5Bx54AAsLC5o2bUpGRsbF7W3btsXd3R2ABg0a0KdPHwCaNWvGtm3bAIiIiOCnn34CYOzYsbzyyisAhIeHs27duovbp02bVjEnK4QQQogqz1hYSMoTT1B85AieH3+EY69eZqtbVVXWnljLvKh5WFpY8l6X97jP7z6MBpUda2M5HJJKHb8a9J0YhKOzza00AAlbYdscSIsy6+QnV3Mmr4TPwhJYszeZMr2Rgc09mNLTn0Z17sBSDEJUkOob7G7gylpFsba2vvjcwsLi4msLCwv0ev01y6uqesv1AOVu2JWbd4UQQghhLCoi5YknKT54EM95H1LjwhfF5pBRmMGMXTMITw+no3tH3ur0FnXt65KfW8Jfy4+ScfI8LXrWo+ODDdBY3uTC3KoKiaGmQJe6F2rWg/vnQ4tRYGmGe/OuIPVsEZ+GJvBDZCoGVWVIK08md29A/doOFdKeEHdS9Q12/xH33HMP3333HWPHjmX16tV07twZgE6dOvHdd98xZswYVq9eXcm9FEIIIURlMBYXk/LkUxTt34/nh3OpYcb77jckbmD2ntnoDDpeb/86DzV+CEVROHUshy0rjmMwGOk7MQj/NrcwPPJkmCnQJUdADU8Y8BG0GlthgS4pu5AlofH8tD8NRYFhbeoxuXsD6jnbVUh7QlQGCXZV3IIFC3j00UeZO3fuxclTAObPn8+oUaOYP38+Q4cOreReCiGEEOJOMxYXk/LUZIoiI/F4/31q3HefWeo9V3KOd/a8w8akjTSv3Zw5nefgU8MHo1Fl7++JRG5IwsXDgX6TgnCqc5PBKCkcQt+FpB3g6A73fQitx4Gl9fWPvQXxmfksConn10PpaDUWjOngwxPd6uNe07ZC2hOiMin/HhpYlQQHB6uRkZHltkVHRxMQEFBJPRI3Q35XQgghRMUxlpSQOnkyhRG78Xj/PWqaaZmlqIwoXt7+MmdLzjK55WQeCXoESwtLis6XsXnFMVJjztLkHne6jmyE1kpz4xUn7zZdoTu5HRzqQOcXoc0E0N7CPXk34Hj6eRZvi2fD0dPYajWM6eDD4138cHOsmPaEqCiKokSpqhp8I2Xlip0QQgghRDViLC0ldcozFEbsxn3OHLOFuviz8UzZOgUXWxcW91pMgIvpC9rT8ef4a/lRSor09BjbhKadPG680pR9EDrHNNulfW3oOweCHwVtxVwxO5RyjoUh8WyJzsDB2pLJ3RvwWOf6ONtXzBBPIaoSCXZCCCGEENWEsayM1GefpXDnTtzfmY3TkAeuf9ANyC7O5umtT2NjacPy3stxd3BHVVUObkkh4ucEHF1sGDqlBbXr3eCskWlRsO1diN8Mdi7Q+21o+xhY2Zulv5eKTMplQUg8YbFZ1LTV8sK9jZhwjy817bQV0p4QVZEEOyGEEEKIasBYVkbas89RuD2Mum/NwslM99gX64t5NuRZcktyWdlvJe4O7pQW6Qj5OobEg1nUb1WbnuMCsLa9gY+N6QdN99DFbgTbWtBrBrSbBNbmn3VSVVUiEnJYGBJPRGIOLvZWTOvXhDEdvHG0kUAn/nsk2AkhhBBCVHFqWRlpz79AQWgodWfOoNaIEWap16gamb5jOkezj/Jxj48JdA0kKyWfjcuOUpBTQqdh/rToVe/6SyydPgyh78GJP8DGCXq+Ae2eAJsaZunnv6mqyvbYLBaGxBN16ixujta8MSCAUe29sbOSj7biv0v+9QshhBBCVGGqTkfaSy9REBJCnTffoNbIkWar+5OoT9iSvIWXg1+ml3cvYvedIeSrGGwctDzwUmvcG9S8dgUZx0yBLvpX08Li3adDhyfB5jrH3QJVVdl8PINF2+I5nJqHR00b3h4cyPDgethob2IiFyHuUhLshBBCCCGqKFWnI23qy+Rv3kKd6dNxHj3abHWvPbGWL499yUONH2Js07HERJxm69fRePg70W9SELaO15hwJDMGtr8Hx34GK0foNg06TAZbJ7P1728Go8rGo2dYGBJHzJl8vJ3teH9oM4a08sLqZhdFF+IuZpZgpyjKCmAgkKmqatAV9ivAfOA+oAiYoKrqfnO0faclJSUxcOBAjh49WtldEUIIIcRdTNXrSXvlFfL/+gu3V6fhPG6s2eoOTwtnzp45dPbszKvtXiV612m2rYqhXpNa9H+q+dWXMsiKhe3vw9F1polQukyFjk+DnbPZ+vY3vcHIb4fTWRQST0JWIfVr2/PRiBYMauGBpUYCnRCXMtcVu5XAIuDrq+zvDzS88GgPfHrh53+GXq/H0lIukAohhBDi+lS9nvRpr5L/50bcXn4ZlwkTzFZ37NlYXtr+Eg2cGvBhtw85EZ5B6OoTeAc60/+JZlheKdTlJJgC3ZEfwNIWOj8PHZ8Bexez9etvZXojPx9IZUloAqdyimhS15FFo1rRP8gdjcV17vUT4j/MLElDVdUwRVF8r1FkMPC1aloNfbeiKE6KorirqnraHO1XlsTERIYOHcqCBQtYsWIFkZGRWFpa8tFHH9GjRw9WrlzJH3/8QUlJCYWFhYSEhDB37lzWrl1LaWkpQ4YMYdasWSQlJdG/f386d+7Mrl278PT0ZP369djaVswaL0IIIYSoulSDgfTXpnP+jz+o/dKLuDz2qNnqzirK4umtT2Nvac/iXotJDD9L2Hex+DRzod+kICwvvVctNxG2z4XD34HGGjpOgU7Pgb2r2fr0txKdgR+iUlkamkDauWKaedbks7Ft6B1QBwsJdEJc1526hOQJpPzrdeqFbeWCnaIok4BJAN7e3tes8P297xOTG2PWTjZxbsK0dtNuqOyJEycYOXIkX375JVu3bgXgyJEjxMTE0KdPH2JjYwGIiIjg8OHDODs7s2nTJuLi4ti7dy+qqjJo0CDCwsLw9vYmLi6ONWvWsHz5ckaMGMG6desYM2aMWc9PCCGEEFWbajRyevrrnP/tN2o//zyuEyeare4iXRFTQqaQV5rHyn4rydijY+faOHybu9JvYhAa7b+GN55NgrC5cHANaLTQ/inTVToHN7P152/FZQa+3ZvMsrAEMs6X0trbidlDgujeqPb1Z+MUQlx0p4Ldlf5XqpdtUNVlwDKA4ODgy/ZXFVlZWQwePJh169YRGBjIrFmzeOaZZwBo0qQJPj4+F4Nd7969cXY2jTvftGkTmzZtolWrVgAUFBQQFxeHt7c3fn5+tGzZEoA2bdqQlJR0509MCCGEEJXq7DffkLd+Pa7PTMH1ySfMVq/BaODVHa8SkxvDgh4LKDvgQPiPcfi1cKXvxCA0f09Ccj7dNOTywCpQNNBuInR+ARzrmq0vfyso1bNq9yk+35FIdkEZHeo78/GIlnRs4CKBTohbcKeCXSpQ71+vvYD026nwRq+sVYSaNWtSr149wsPDCQwMxDTC9Mrs7e0vPldVlddee40nnij/hzopKQlra+uLrzUaDcXFxebvuBBCCCGqLN2ZM2TNX4B91y64Tp5s1rrnRc1jW8o2Xm33Ko7RvoT/FE+DVrXp/XggGo0FqCrs/wo2vQn6EmjzCHR5EWp4mLUfAHnFOr7alcSK8JOcK9LRtVFtnunpT1tf80/AIsR/yZ0Kdr8CUxRF+Q7TpCl51fn+OisrK3755Rf69u2Lg4MDXbt2ZfXq1fTs2ZPY2FiSk5Np3Lgx+/eXn/izb9++vPnmm4wePRoHBwfS0tLQarWVdBZCCCGEqEoy3nkH1Wik7v/+Z9YrVmti1vDN8W8YHTCaJqc6EfFLAv7Bbtz7SFNTqMs9Cb89CyfDwLcLDFoAzvXN1v7fzhaW8cXOk3y1K4n8Uj33BrgxpWdDWtYz/xIJQvwXmWu5gzVAd8BVUZRUYAagBVBVdSmwAdNSB/GYljt4xBztViZ7e3t+//13evfuzRtvvMHhw4dp1qwZlpaWrFy5stwVuL/16dOH6OhoOnbsCICDgwOrVq1Co5FFNYUQQoj/svyQEPI3b6H2Sy9i5eVltnrDUsN4b+97dPfqTs+sEez+LZGGbetw74QALBQVIpZAyNumYZcDP4HW48HCvEsJZOWX8vmORL7ZfYpinYH+QXV5uoc/gR7mX8RciP8y5VrDCCtTcHCwGhkZWW5bdHQ0AQEBldQjcTPkdyWEEELcGGNhIQkD70fj4IDfT+tQzDSaJyY3hvF/jsfH0Yfn1Lc5+GcajdvXpef4ACxyYmH9FEjdCw37wMCPoab5AiXAmbwSlm5PYM3eZHQGI/e38GBKD38a1nE0aztC3M0URYlSVTX4RsrKwmpCCCGEEJUoa9Fi9KdP4/ntt2YLdRmFGTy99WkctY48qXuTg5vTaNKxLj1G+WOxc55pghQrexiyDJqPADMO/UzJLeLT7Qn8GJmKUVUZ0sqTyT388XO1v/7BQohbJsFOCCGEEKKSlERHk/v11zgNH45d61ZmqfPvZQ0KSgv4n2Yh0SFZNO3kTvdeZShf9IIzh6HpA3DfXLMuX3Ayu5Al2+L5+UAaForCsGAvnurWgHrOdmZrQwhxdRLshBBCCCEqgWowcHrGTDROTri99KJZ6jQYDbwc9jKxubFMV+eTtLOAwE516Ob+I8rnn4CtM4z4BpoOMkt7AHEZ+SzaFs9vh9LRaiwY08GHJ7rVx72mrdnaEEJcnwQ7IYQQQohKcPb77yk5fBiPuR+gcTLPzJAf7PuAsJQwni99l+wDRpq10dLl3CMoCSegxSjo+w7YmWdZgWPpeSzeFs+fR89gq9UwsUt9Huvih5ujjVnqF0LcHAl2QgghhBB3mC4zk6yPPsb+no7UGDjQLHWujl7Nt9Hf8tj51yk5bkfz+sl0Tn0epaYnjF4HDe81SzuHUs6xMCSOLdGZOFpb8nR3fx7t7IezvZVZ6hdC3BoJdkIIIYQQd1jGu++ilpWZbc26bcnbeH/PB4zMeh5tghstnLfRqXABSrvH4d6ZYH37M1HuS8plYUg8YbFZONlpebF3I8bf40tNW1mTV4iqQIJdFefg4EBBQUFld0MIIYQQZlKwYwf5f27E9dlnsPL1ve36juUcY1rYqwxKexynFD9a2f9ExzrhKIP/AN/Ot1W3qqpEJOSwICSO3Ym5uNhb8Wr/Jozp4IODtXyMFKIqkf+RVZSqqlQj+CE1AAAgAElEQVTVNQaFEEIIcWuMxcWcmfUWVn5+uDz++G3Xd6bwDM9seYbu8cNxPxNIG/t1tO9TC6VHOFjd+myUqqoSGpvFopB4ok6dxc3RmjcHNmVUO29srTS33W8hhPlZVHYHqpvCwkIGDBhAixYtCAoK4vvvv8fX15fs7GwAIiMj6d69OwAzZ85k7Nix9OzZk4YNG7J8+fKL9cydO5e2bdvSvHlzZsyYAUBSUhIBAQFMnjyZ1q1bk5KSAsBLL71E69at6dWrF1lZWQAkJCTQr18/2rRpQ5cuXYiJiQEgKyuLoUOH0rZtW9q2bUt4ePidemuEEEIIcR3Zny5Fl5pK3ZkzsbC6vXvSCsoKmLzpSVoe7YPPmdYEu4XQ/vlHUPrOvuVQZzSqbDp2hsGLw3nky32cySvh7QeCCHulB4919pNQJ0QVVm2v2J2ZM4fS6Biz1mkd0IS606dfs8zGjRvx8PDgjz/+ACAvL49p06Zdtfzhw4fZvXs3hYWFtGrVigEDBnD06FHi4uLYu3cvqqoyaNAgwsLC8Pb25sSJE3z55ZcsWbIEMAXJ1q1bM2/ePN566y1mzZrFokWLmDRpEkuXLqVhw4bs2bOHyZMnExISwnPPPccLL7xA586dSU5Opm/fvkRHR5vvTRJCCCHELSmJjSVnxQpqDhmCfft2t1WX3qhn6u/j8d7fkQbZbWkXlELbJ98ES+tbqs9gVPnz6GkWhcQTcyYfHxc73h/ajCGtvLCylOsAQlQH1TbYVZZmzZoxdepUpk2bxsCBA+nSpcs1yw8ePBhbW1tsbW3p0aMHe/fuZefOnWzatIlWrUwLkRYUFBAXF4e3tzc+Pj506NDh4vEWFhY89NBDAIwZM4YHH3yQgoICdu3axfDhwy+WKy0tBWDLli0cP3784vbz58+Tn5+Po+Pt3zQthBBCiFujGo2cmTkLjYMDbq+8fHt1nT/Du+sfxibmPhrmtKF9T3uCR4y/pbr0BiO/Hkpn8bZ4ErIKaVDbno8fasH9zT2w1EigE6I6qbbB7npX1ipKo0aNiIqKYsOGDbz22mv06dMHS0tLjEYjACUlJeXKXzrTlaIoqKrKa6+9xhNPPFFuX1JSEvb29tdsX1EUjEYjTk5OHDx48LL9RqORiIgIbG1lUVAhhBCiqji3bh3F+/fj/s47WNaqdWuVqCoc+o6vwmZyNmMSDXNa0XGwH637+910VWV6Iz8fSGVJaAKncopoUteRxaNa0y+oLhqL25+lUwhx58lXMTcpPT0dOzs7xowZw9SpU9m/fz++vr5ERUUBsG7dunLl169fT0lJCTk5OYSGhtK2bVv69u3LihUrLs52mZaWRmZm5hXbMxqN/PjjjwB8++23dO7cmRo1auDn58cPP/wAmG5wPnToEAB9+vRh0aJFF4+/UvgTQgghxJ2jz8kh88N52AUHU/PBIbdWybkUWD2MzX9N5ciZJ2mQ04qOQxvcdKgr0Rn4JiKJHh+GMm3dEWrYaFk2tg0bnu3CgObuEuqEqMaq7RW7ynLkyBFefvllLCws0Gq1fPrppxQXF/PYY48xZ84c2rdvX658u3btGDBgAMnJybz55pt4eHjg4eFBdHQ0HTt2BExLGqxatQqN5vIbku3t7Tl27Bht2rShZs2afP/99wCsXr2ap556itmzZ6PT6Rg5ciQtWrRgwYIFPP300zRv3hy9Xk/Xrl1ZunRpxb8xQgghhLiijPffx1hURN1ZM29+zTqjEaJWwOYZHNJYsSVvGvVzg2g/1JfWvX1uuJriMgOr95xiWVgimfmltPGpxTtDgujWqLZZ1tETQlQ+papOqR8cHKxGRkaW2xYdHU1AQEAl9ejmzZw5EwcHB6ZOnVrZXbnjqtvvSgghhKgIhRERJD/yKC5PPYnbc8/d3ME5CfDrs3BqJ8k+Xfks7j48chvT+kFPOvZpfENVFJTq+SbiFJ/vSCSnsIyO9V14pqc/HRu4SKATohpQFCVKVdXgGykrV+yEEEIIISqAsbSUMzNnofX2xvWS++qvfaABdi+BkHdAY8XZvvP48g8LPHJ9afqA8w2FurxiHSvDk1gRfpK8Yh1dG9Xm2Z7+BPs638YZCSGqMgl2FWjmzJmV3QUhhBBCVJKcZcspO3WKel98joWNzY0dlLYfNrwMaZHQeADF977L0oU7cc12x2egFT36tbzm4bmFZXyxM5Gvd50iv1TPvQF1eKanPy3qOZnhjIQQVZkEOyGEEEIIMytNPEnOsmXUGDgQh06drn/AuRTY+hYcWQv2tWHYCsr8B7H0/V9xzKyDU58iBg7sedXDM/NL+HzHSVbtPkWxzkD/oLpM6dGQph41zHhWQoiqTIKdEEIIIYQZqarKmZkzUWxtqfPqtGsXLjkPOz+CiCWgKNDlJej0PDrFni8+2Ij2dE3U7umMfnDcFQ8/nVfMZ9sTWbM3GZ3ByKAWHjzdw5+GdWT9WiH+ayTYCSGEEEKYUd769RTt3UvdmTOxdHW9ciGDHvavhG3vQlE2NB8Jvd6Eml6Ulej55sNt6NOsye14hDcfunzSlZTcIj7dnsCPkakYVZUHW3syubs/vq7XXg9XCHH3kmAnhBBCCGEm+rNnyXz/A2xbtsRpxPDLC6gqxP4Fm9+E7Fjw6Qx9Z4NHKwCK88v4ceFuilJUEtrs4P0xr2Oh/LPs8MnsQhZvi+fnA2loFIXhwV482a0B9Zzt7tQpCiGqKAl2lWTChAkMHDiQYcOGVXZXhBBCCGEmmR9+iCE/n7qzZqFYWJTfefoQbHoDToaBiz+MXAON+4OioBpVjoens/OnWEpL9Bxs/gcfjv0f1hprAGIz8lkUEs/vh9PRaiwY19GHJ7o2oG7NG5yURQhx1zNLsFMUpR8wH9AAn6uq+t4l+72BrwCnC2VeVVV1gzna/q/T6/VYWko+F0IIISpbUWQkeet+wuXxx7Bp3OifHefTYevbcGgN2NaC/nMh+BHQaAHITs0nZNVxspIKSa8RR0yLUD4cMgdnG2eOpuWxeFs8fx49g52Vhold6/N45/rUdrSupLMUQlRVt50IFEXRAIuB3kAqsE9RlF9VVT3+r2JvAGtVVf1UUZSmwAbA93bbrgzTpk3Dx8eHyZMnA6YlDRRFISwsjLNnz6LT6Zg9ezaDBw8mKSmJ/v3707lzZ3bt2oWnpyfr16/H1ta2XJ1RUVG8+OKLFBQU4OrqysqVK3F3dychIYGnn36arKws7OzsWL58OU2aNGHChAk4Oztz4MABWrduzbx58yrjrRBCCCHEBWpZGadnzETr4YHrhc8IlOZD+ALYtRBUA3R61jQ5ik1NAMpK9Oz9LZFDISmUWBYR4f8L7bo25ptWXxB7pozHft7H1phMHK0teaanP4928qOWvVUlnqUQoiozx6WedkC8qqqJAIqifAcMBv4d7FTg7/l2awLpt9vojrWxZKcU3G415bjWc6DLiEbXLDNy5Eief/75i8Fu7dq1bNy4kRdeeIEaNWqQnZ1Nhw4dGDRoEABxcXGsWbOG5cuXM2LECNatW8eYMWMu1qfT6XjmmWdYv349tWvX5vvvv+f1119nxYoVTJo0iaVLl9KwYUP27NnD5MmTCQkJASA2NpYtW7ag0WjM+h4IIYQQ4ublrFhBWUIC9T5bioW1FUStNC0wXpgJQUOh1wyo5QOYZs1MPJBF6PcxlOTpOe62i7MtYpnV/VUKz7vzxNdH2BGXjZOdlpd6N2LcPb7UtNVW7gkKIao8cwQ7TyDlX69TgfaXlJkJbFIU5RnAHrj3ShUpijIJmATg7e1thq6ZX6tWrcjMzCQ9PZ2srCxq1aqFu7s7L7zwAmFhYVhYWJCWlkZGRgYAfn5+tGxpWky0TZs2JCUllavvxIkTHD16lN69ewNgMBhwd3enoKCAXbt2MXz4Pzdel5aWXnw+fPhwCXVCCCFEFVCWnEz2p0tx7NsXBw8dfNYFMo9DvQ7w8BrwCr5YNi+rmO3fxZBy7Cy59ulEtPiZkV0H42s1ibd+TGTPySRcHax4tX8TxnTwwcFabrcQQtwYc/y1UK6wTb3k9cPASlVV5ymK0hH4RlGUIFVVjeUOUtVlwDKA4ODgS+so53pX1irSsGHD+PHHHzlz5gwjR45k9erVZGVlERUVhVarxdfXl5KSEgCsrf8ZA6/RaCguLi5Xl6qqBAYGEhERUW77+fPncXJy4uDBg1fsg729TGcshBBCVDZVVTkz6y0UjQV1GsfD6qFQyw9GfA0Bg0xr0wEGnZEDm5PZtyGRMrWMvT5/4BSs8pjbTFZtP8+B5Ejq1LDmfwOb8nA7b2yt5MtbIcTNMUewSwXq/eu1F5cPtXwM6AegqmqEoig2gCuQaYb277iRI0cyceJEsrOz2b59O2vXrsXNzQ2tVsu2bds4derUDdfVuHFjsrKyiIiIoGPHjuh0OmJjYwkMDMTPz48ffviB4cOHo6oqhw8fpkWLFhV4ZkIIIYS4Ged/WkNheDh12pxHez4H+r4LbR8Hy3/uhUs9cZbQb6PJyygh0fkQMU1C6eY3mu0H6/LS1mQ8nWyZ/UAQw4O9sLaUQCeEuDXmCHb7gIaKovgBacBIYNQlZZKBXsBKRVECABsgywxtV4rAwEDy8/Px9PTE3d2d0aNHc//99xMcHEzLli1p0qTJDddlZWXFjz/+yLPPPkteXh56vZ7nn3+ewMBAVq9ezVNPPcXs2bPR6XSMHDlSgp0QQghRFZQVYtgyj4y3V2PjbKDW6HHQ/WWwc75YpOh8GeHr4ojdk0Gh7Vm2N1mLu399yuIns+S4Hh8XAx8Mbc6Q1p5oNRbXaEwIIa5PUdVrjni8sUoU5T7gE0xLGaxQVfUdRVHeAiJVVf31wkyYywEHTMM0X1FVddO16gwODlYjIyPLbYuOjiYgIOC2+ysqnvyuhBBC3JWMBtOyBSGzOb2tiHMJ9viumI9txz4Xi6hGlWM709n1cxxlpXr2u28mtX4sxdmDSDnthr+bA1N6+DOwuTuWEuiEENegKEqUqqrB1y9ppnXsLqxJt+GSbf/71/PjQCdztCWEEEIIUSkStsGmNyHjCMVKc84l5OA8bmy5UJeVnE/otzFkJuVzpmYiOxr/SKG+FVlHJhDg7syS0f70C6yLhcWVpigQQohbJ1MtCSGEEEJcS2Y0bP4fxG0CJ2/UB5Zzesa3WNaxxPWZZ4ELa9L9epJD21Io0xazw/9HEh3Pk5f6MM3qNGDOuIbcG+CGokigE0JUjGoX7FRVlT+KVZw5hvcKIYQQla4gE7bNgf1fgZUj9H4L2j1B7terKT1xAq9FC7GwtyM+KpMda2MpzCslus4udntuJS/7Xpobe/DsmIZ0a1RbPrsIISpctQp2NjY25OTk4OLiIn8gqyhVVcnJycHGxqayuyKEEELcGl0xRCyGnR+DvsQ0y2W3V8HeBV1aGlmLFuPQsyfGFvfw+6JDJB/LJdf+DKFB35JmcKUZr/PCQ63oWF8+rwgh7pxqFey8vLxITU0lK6vaTqj5n2BjY4OXl1dld0MIIYS4OfpS2P817PgI8tOh8QDoPQtcGwIX1qx7622MFpakd32C32ftoUwtY4/vbxyudZxG2gm833Mgwb7O12lICCHMr1oFO61Wi5+fX2V3QwghhBB3E30ZHFwFYfPgfCrUaw9Dl4Nv53LF8jdtJvnQaeI6v0VhaDbxzofZ5fsztaw782W3NbT1rVNJJyCEENUs2AkhhBBCmI1BBwdXmwJdXjJ4tYVBC6BBT7hkCGX+6bNs/vI4Z1o+R556jp1NviCvlpG3Oy6gTyNZY1YIUfkk2AkhhBDiv8WgM61FFzYXziWDZxsY+DH497os0BmNKrv+SuLwL7FQI4DEGn+xvclOHm3+FE+2GYWFIuvQCSGqBgl2QgghhPhvMOjh8HemQHc2CTxawX3zoGHvywIdwJEjWWz5Jhqr83qcziWSqf2e4mGt2NhlPa62rne+/0IIcQ0S7IQQQghxdzPo4cha2P4BnD0J7i3g4e+gUb8rBroTKedY99VxHFOL0WtK8Ej/Du/T+/H55iNeaNK3Ek5ACCGuT4KdEEIIIe5ORgMc+RG2vw+5CVC3GYz8Fhrfd8VAdyzpHGu/j8Y+qYgaKsTU3Ytj1k8MjC3Ade671JZQJ4SowiTYCSGEEOLuYjTA0Z9MgS4nDuoEwUOroMnAKwa6/bHZrF97ghqpJbihkFErjY1e39OgpjVP/6rHvksXXAcOroQTEUKIGyfBTgghhBB3B6MBjv1sGnKZfQLcmsKIr6HJ/WBx+SQnu49msOHHWJzPlFEbSHdNJML9J/QuBTzZ4km6LtlNofEkdWf8TxYaF0JUeRLshBBCCFG9GY1w/BfTFbqsGKgdAMNXQsDgKwa67fvSCPklAZccHa6oJNc5zh739di6ang8cDxD/Idg2LGH1M1bqP3ii1h5ed35cxJCiJskwU4IIYQQ1ZPRCNG/mgJd5nFwbQzDVkDTIZcFOlVV2Ryewq7fE3E9Z8BZMRLnHkWk+wbq1nHhlaDn6efXD62FFmNREQmz38a6oT8uj0yonHMTQoibJMFOCCGEENWL0Qgxv5sCXcZRcGkIQ7+AwCFgobmkqJENIUlEbTyFa4FKTQsdR7x2sb/uJpp4+PNes7fp4tWl3Hp0WQsXoU8/jee3q1G02jt9dkIIcUsk2AkhhBCielBVOLEBQt+FM0fAuQE8uByChl4W6AwGI7/8Gc/xrWk4F6vYW5YQ6b2Nw3VC6eDbjs+CFtO6TutyxxjLysic+yFnv/kGp+HDsWtdfr8QQlRlEuyEEEIIUbWpKsRuNAW604fAuT48sBSaDQdN+Y8yOr2Rdb/EkrgjnZqlYKUtZKffX8S67aG3fy++DVxFY+fGlzVRlppK2gsvUnLkCM7jx+H20kt36uyEEMIsJNgJIYQQompSVYjbZAp06Qegli8MXgLNH7os0JWU6Fm7LobTuzNx0IFqfY4Q/z9IdjvCA40GMzfwZ7wcrzwJSv6WLaRPfx1UFa9FC3G89947cHJCCGFeEuyEEEIIUbWoKsRvMQW6tChw8oFBi6DFSNCUv+etsLCM77+PJicqGzuDQrFdJjv9fiOnbhIjmzzE6IAPcLF1uXIzZWVkzptH7ldfYxMUhOcnH8sMmEKIakuCnRBCCCGqhr+HXIZ9CGmRUNMb7l8ALUddFujyzpWwdk00+UdysTYqnHNMYavXb5TVPcu4wHEMazQMByuHqzZVlppG2osvUnL4MLXGjsXt5alYWFlV9BkKIUSFkWAnhBBCiMplNMDx9bDjI8g4YrpCN/ATaDkaLMuHreysQn5YHU3piTy0qsIZpziivH7H1gMmBj3C/Q3ux0pz7YCWv3Ur6a9NB1XFc8F8avTpU5FnJ4QQd4QEOyGEEEJUDoMOjvwIO+ZBThy4NoIhn0HQsMvuoUtPPc/P38ZgSMxHQeWU6xGiPP/EvV4tpgU9Sy/vXmgumRnzUmpZGZkffUzuypXYBAaahl7Wq1eRZyiEEHeMWYKdoij9gPmABvhcVdX3rlBmBDATUIFDqqqOMkfbQgghhKhm9KVwcDXs/BjOJUOdZjB8JQQMumzZgqT4c/z6fTRKShFGxUis2z4OeG6iqZ8/7wfNooN7BxRFuW6TurQ0Ul98kZJDh6k1ejRu016RoZdCiLvKbQc7RVE0wGKgN5AK7FMU5VdVVY//q0xD4DWgk6qqZxVFcbvddoUQQghRzZQVQtRXsGsB5J8Gz2DoPxca9YVLwtmJo1ls/CEWy4xSDIqO4x47OeK+jXsatmNZ0CKCXINuuNn8kG2kv/YaGAx4fvIJNfr1NfeZCSFEpTPHFbt2QLyqqokAiqJ8BwwGjv+rzERgsaqqZwFUVc00Q7tCCCGEqA5KzsO+5RCxGIpywLcLDFkKft3KBTrVqHJoz2lCf0tAm6tDpylhf70Qot13cV+TPkwP/Aa/mn433Kyq05H58SfkrliBddMAvD7+GCsfn4o4QyGEqHTmCHaeQMq/XqcC7S8p0whAUZRwTMM1Z6qquvHSihRFmQRMAvD29jZD14QQQghRaYpyYfensPczKMkD/97QdSp4dyhXTFdmYE9IMlGbkrEsMlCsPc9e3y2cdD/IsMAhfBDwM3Xs69xU07rTp0l74UWKDx6k1qiHcZs2DQtra3OenRBCVCnmCHZXGtiuXqGdhkB3wAvYoShKkKqq58odpKrLgGUAwcHBl9YhhBBCiOogPwMiFsK+FaArhID7octL4NGqXLHCvFLC/0wiJjwdjU4lxy6dww3/Itc9hdGBD/Npk3eoaV3z5psPDeX0tFdR9Xo8P5pHjfvuM9eZCSFElWWOYJcK/HtKKS8g/QpldquqqgNOKopyAlPQ22eG9oUQQghRFZxLgfD5sP9rMOpMs1t2eRHcAsoVy0krYMcfiaQcyEZRVU45RXPYczM2nhY81nw89/ndd90lC65E1enImj+fnM+/wDogAK+PP8LK19dMJyeEEFWbOYLdPqChoih+QBowErh0xstfgIeBlYqiuGIamplohraFEEIIUdlyEkwzXB5aAyjQ8mHo9Dy4NLhYRFVVUqPPEv5HIjkJ59Epek7U3sMRj2009m3IrObT6Oje8YZmuLwS3enTpL00leL9+3Ea+RB1XntNhl4KIf5TbjvYqaqqVxRlCvAXpvvnVqiqekxRlLeASFVVf72wr4+iKMcBA/Cyqqo5t9u2EEIIISpRZrRpDbqj60BjBcGPwj3PgtM/A3kMOiNxkRns/jOJwsxiCi0LOVpvGzF19tKzYQ9WNl+Gfy3/2+pGQVgY6a9MQy0rw+PDD6k5cMDtnpkQQlQ7iqpWzVvZgoOD1cjIyMruhhBCCCEulX4Awj6EmN9Baw9tH4OOU8DxnwlOSgp1HA1LY//WZHQFerJtMjniuZkUtzhGNh3BuKCHcbF1ua1uqHo9WfMXkLN8OdaNG+P5ycdY+934rJlCCFHVKYoSpapq8I2UNcsC5UIIIYT4D0jeDWFzIX4L2NSEbtOg/ZNg53yxyLnMIg5vTeFoeDqqXiW5RhyHAzZT6FrEE60eYVjjJVhrbn+IpC4jg7QXX6I4KgqnESOoM/01LGxsbrteIYSoriTYCSGEEOLqVBUSQ01X6E7tBDtX6DUD2j4ONjUuFFE5nZDHwS3JnDyYjUExEOsaxWGPEByc3Zja/ll6+nS55fvnLlWwYwfpr0zDWFqKx9wPqHn//WapVwghqjMJdkIIIYS4nKpC7EbTFbq0KHB0h37vQevxYGUHgNFgJOFAFge3JJOZlE+ppoSjnmEcrRNOfdc2LOmygKDaTczXJb2erIWLyPnsM6wbNsRz/idY169vtvqFEKI6k2AnhBBCiH8YDXD8F9jxMWQcAScfGPgJtBwFlqYhlGUleqLDT3NwazIFuaXkWZ/lsN8WYl2O0rZuP37u9iOejje3oPj16DIySX/pJYoiI6k5bCh1X38dC1tbs7YhhBDVmQQ7IYQQQoC+FA5+a1qH7uxJcG0EQz4zrUWnMX1cyM8t4fC2VI7tSEVXYuS0YxKHGm8m1TGbvt7DWd5tHo7WdmbvWsHOcNJfeQVjcTEe779HzcGDzd6GEEJUdxLsxP/Zu+/wqqp8/+PvfU5670ASAqTQCS1UaUJoKlJU7OM4KvYyOj+n3nvn6p0Zp1wbYBs7zthGUGw4VOklQEILkN57PWmn7t8f8TrqoKNDJAE+r+fJk2z3Wnvtb/7g5ONeey0RETmf2W2w/2XYuQJaKiF2NMxeBYMvAYsFgJpiGwfXF5O7vwqP6SEv4hBZKRtp8vHnykHXct/kBfh4df2fFKbLRc3KldQ98yy+yUnEPf44vklJ/7qjiMh5SMFORETkfNRaB3ufhT3PQkcjDJgGi5+BxBlgGJgek8JDtWRuKKb8ZCMuq5OjvbZzqPd2TEsyy1J/zo3jLsBi6ZoFUb7KWV1N+QM/oW3fPkKXLKH3f/xKUy9FRL6Bgp2IiMj5pKm08+ncgVfA2db5ZG7K/RA/FgCXw83x3ZVkbiymqaqdNl8bmf02kh2RSaAxgf83fiVLRg7pshUuv8rjcGD7+3qqfvc7PG1t9Pnd7whbvOh7GUtE5FyiYCciInI+qM2B7Y/DoTcBE0ZcARfcBzGdq1a2NTs4vKWUw5+WYm91URdUzsGU9eQGldLLks7vp/yE9MF9v7dA13HsGI3vrKb5gw9wNzXhm5JM3Csv45uc/L2MJyJyrlGwExEROZeVH4Rtj0L2+52rWqbdCJPvhrAEAOrKW8jaWMKJPZW4XR6Kw4+R2X8jJRYPyX4X88z0JUxOivpeAp2roYHm99+ncfUa7MePY/j4EJw+i9DFSwicPAnDau3yMUVEzlUKdiIiIuca04TCbZ2BLn8z+IbC1Pthwu0QFI1pmpRm15O5oZjio/V4rG6yo3ZxqPdWap19GBF0Db+bmc7YfhFdf2suFy3bttG0eg22LVvA6cRv+HB6/ed/EHrRRVjDwrp8TBGR84GCnYiIyLnC44GTH3cGurIMCIyB9F9D2k3gF4Lb6SFnVwWZG4qpK2vF4dtGZt9NHIk6gK1lOJMCHuTpWRMYER/a5bdmz8ujcfVqmtauxV1TizUigohrryV08WL8Bg3s8vFERM43CnYiIiJnO7cTjrzT+Q5dTXbnpuIXPwqjrgVvPzpanRz5uJBDm0tob3bSHFjD/qS/czK0gI7GicwM+i33XjqcIX1Cuva2bDaaP/yIxjWr6cg6BFYrQdOnE7ZkMUHTp2N4e3fpeCIi5zMFOxERkbOVsx0OvgY7n4TGYogZCkueh2GLwepFY3UbWRtPkL2rHLfDpCz8JAeHbKDYx46zYSoXh9zJnYsGkhwT1GW3ZHo8tO3eTePqNdjWr8e02/FNSSbmwQcJvXQBXlFRXTaWiIj8g4KdiIjI2aajCfY9D7ufhtYa6DsB5v8RUuZgGgYVuU1kbiim4FAtpuHhZHgvawUAACAASURBVNQ+svpsocYVhasxnSVDJ3P70iT6RQZ22S05SkpoWrOGxnffxVVegSUkhNAliwlbsgS/4cO/t9U0RUSkk4KdiIjI2aKlGnY/BfteAHszJKd37kHXbzIej0ne/hoObiimpsiG06eDQ7Gfkt17Lw0tQ/BU3cCVo1K57fok4sK6ZqNvT1sbzX//O02r19C2dy8YBoGTJxPzwAMEp6dj8fXtknFERORfU7ATERHp6RqKOqdbHnwNXHYYtgim/Bj6jMTe7iJ7YwmZG4tpbXDQElDPgQEbKIzOo6F2LJaSe7h2/CCWTUukV4jfad+KaZq0HzxI4+rV2D5eh6e1Fe+EBKLvu5fQhQvx7tOnCwoWEZHvSsFORESkp6rOhu2PweG/gWGBUVd3bioemYStvoOsv+VwdFsZLruHqtACDgxaT114OzVlE/EtXMBNk5O4acoAooJO/8mZs6qKpnffo2nNGhyFhRgBAYTMm0fYksX4jx2rqZYiIt1MwU5ERKSnKdnbGehOfATegTDxdph0J4TEUlXYTOaaI+QdqMZjesiLzCRr4CacgdGUFU0jsCmROy5I5EcX9CcswOe0bsPjcNCycSONq9fQumMHeDz4p42lzy23EDJvLpbArntHT0RETo+CnYiISE9gmpCzHnY8DkU7wD8cZvwcxi/D4xdO4aFaDq7fT2VeE24vB4d7beNk3B4My3CKCq4izLs3D0xL5PpJ/Qjx+/e3ETDdbtozM2n+8COaPvwQT1MTXr17E7nsFsIWL8anX78uLFpERLqKgp2IiEh3crvg6JrOQFd1BELiYd4jMOYHOE0/ju+qIHPjcZprOmj3a+ZA/w1UxOVgsU+kKP9WogLC+cXsRK6ZkECg77/3se6x22ndtQvbhg20bNqMu74ew8eH4PR0QpcsIXDSRAyrtYsLFxGRrqRgJyIi0h0cbZD5l3/sQRc9GBY9AyMup7XFw+GPSjm8tRRHm5u64FL2D1yPM64Ne+MUCo7OoXdIIL++OJGrxifg5/3dQ5fbZqNly6fYNm6kdetWPG1tWAIDCZo+neD0WQROm4Y1qOv2txMRke9XlwQ7wzDmAU8AVuB50zQf+Zp2lwNvA+NM08zoirFFRETOKu0NsPd52PMMtNVC/HiY/wdImUtdRRuZq3I4sa8Sj9ukMPIwWYmbiIjvja10BvkHexEfHsBvFydz2dg4fL2+W6BzVlXTsnkTtvUbaN27F5xOrNFRhCxYQHD6LAImTMDic3rv5YmISPc47WBnGIYVWAnMBkqBfYZhrDVN89hX2gUD9wB7TndMERGRs05zOexaCftfBkcLpMyFKfdh9p1IyfEGMpdnUZLdgMfq4mj0Tk7G7WJAnzG4cq9jz+5ABkQF8sfLk1g0Og5vq+VbD2svKMC2YQO2DRvoyDoEgHe/BCJ+cD3B6en4jxyJYfn21xMRkZ6pK57YjQdyTdPMBzAM4w1gIXDsK+0eBv4A/KQLxhQRETk71JyEnU9A1ptgemD4ZXDBvbgjh3JyXxWZr+ylvrwNu28LmQmbqUw4zuDImbiP3cPGPBjYK4gnrkrmktRYrJZ/vaWAaZp0HDmCbf0GbBs34sjLA8Bv2DCi77uX4Fmz8ElO1vYEIiLnmK4IdnFAyReOS4EJX2xgGMZooK9pmh8YhvG1wc4wjGXAMoCEhIQuuDUREZFuUroftj8Kxz8EL19IuxEm3UWHTyxHt5WRtWk77c0uGgOrOJC0HjOpmYFBF5F3cD5rjzgZ2ieEp69NZu6w3lj+RaAznU7a9u3DtmEjto0bcVVVgdVKwLhxhF99NcGzZmrjcBGRc1xXBLtTfdqYn580DAvwGPDDf3Uh0zSfA54DSEtLM/9FcxERkZ7FNCFvU+cedIXbwC8Upv0Ext9KU3sgWZ+UcmzHdtxOk7KwExwcspGEwVGMsCzg433BHGi2M7JvIA9dmszMwTHf+FTN09ZGy/btnStZbvkUT3Mzhp8fQVOnEDTrPoJnzMAaFnYGixcRke7UFcGuFOj7heN4oPwLx8HAcGDLZx9QvYG1hmFcqgVURETknOB2QfZ7sP1xqDwEwbEw5zcw9gYqSj1kvlZMfmYNpuHhRNQ+suO2M2HYWKbZ7+Sd3W5qW+yM7x/In64YyZTkqK8NdK6GBlo2bca2YQOtO3di2u1YQ0MJnjmT4NnpBE6ejMXf/wwXLyIiPUFXBLt9QIphGAOAMuAq4Jr/O2maZhMQ9X/HhmFsAX6iUCciImc9Z8dnWxYsh4YCiEyBhSvxDLuC/MNNHHz8ONUFNpzeHRyO20ppv8PMHzqfhIZf89fNDTS0tTElOYq7Z45mQmLkKYdwlJbRsmkjtvUbaNu/HzwevGL7ELZ0KcGzZhGQNhbDS7sXiYic7077k8A0TZdhGHcBn9C53cGLpmkeNQzjISDDNM21pzuGiIhIj9LRBPtegN1PQ2s1xI2FOQ/j6DeX7N1VZP73PlrqHLT4NXBwwAYcSdUsHryUivJrefGDcpo7qpk5OIY7L0xmbL/wf7q8o7SUprVrsW3YgP1YNgC+KSlE3XYrQbNm4Td0qBY/ERGRLzFMs2e+ypaWlmZmZOihnoiI9CC2Stj9FOx7ERw2SJoFU35Ma9g4Dm0p4/DWEpztHqpDCjnYeyO9hgWyKOkqDp3szWu7i2l1uJk7rBd3z0xheFzoly7taW2l+e/raVqzhra9e8Ew8B89muBZswhOn4VPv37dVLSIiHQXwzD2m6aZ9m3aau6GiIjIv1KXBzuegKzXweOCYYvhgnupdSeRub6Ykxk78XhM8sOzOJLyKeNSh/Ozfj9lQ6YXD7xShN1VwCWpsdx1YTKDegd/flnTNGnPyKBx9RqaP/kEs60N734JRN93L6GXXop3bGw3Fi0iImcTBTsREZGvU36wc0GUY++B1QdGX4c56S6KK8M4+GYRZcf34rY6ORa9k/yE/Vwyag4/7LOct/c088N1Jbg9JotGxXHHhUkkRQd9fllHaRlN771L07vv4SwpwRIYSMhF8wlbsgT/0aM1zVJERL4zBTsREZEvMk3I3wI7Hu/87hsCU36MO+1WThz1cHBFEY0VRbT72shK2IItuZirU5fy4/C7eXFbOUv+dgTDgMvHxnP79GQSIgOAzu0JbOvX07h6DW179oBhEDBxAtF330VwejqWgIBuLVtERM5uCnYiIiIAHnfnk7kdj0NFFgT1hvT/pmPI9RzZYyPzNyew29zUB1aQmbyR0KEGPxp+HXG+Y3lmSz6/yNqL1WJw7YQEbp2eRGyYP6Zp0paRQeOaNdg+XoenrQ3vhASi772nc6plXFx3Vy0iIucIBTsRETm/Odsh869f2LIgGRY8SWPsQrK2VHHszSw8TigOy+bw0C2MGJnIQ8MewOLsy4rNuXx0eBt+XlZ+dEF/bpmaSEyIH87ycmpff5nGNe/iLC7GEhBA8Px5hC1ejP/YsZpqKSIiXU7BTkREzk/tjbDvedjzDLTWdG5ZMPshKryncHBDKQVZ+/FYPJyM3Eduwl5mj57Gj4csp6bBn8c/ymH9sW0E+Xpx+/QkbpoygHCrB9v6v1O0Zg1tu/eAaRIwcSLRd95B8OzZmmopIiLfKwU7ERE5vzSXw66VsP9lcLRAcjqeSfdS0DSY/e8XUlOQicOrncNxW6lLyuWKUYt5JOVejpfb+dlbuWw5UUOInxf3pafww0n98TlxhKbf/w85H32Mp7UV7/h4ou66k9CFi/CJ11RLERE5MxTsRETk/FBzAnY8CYfeBNMDw5fgHHcPx/Mj2P9iIa11R2jxq+dg/014D23l+tRrmdn3N+wrbOKWlw+zM6+OiEAfHpw3iKv7+eD+5CNq/7gGZ1ExRkAAIXPnErbks6mWFkt3VysiIucZBTsRETm3lezt3LLgxIfg5Q9pN9I2/FYOZ1rIeqwIZ3st1UFFZA3cTOKoKH42fBkjo0eyNaeWq/+8l32FDUQH+/Kf6QNYYMuh/dXfUrlrd+dUy/HjibrtdkLmzMYSGNjdlYqIyHlMwU5ERM49pgk5f+8MdMU7wT8cpv+U+n4/IHOHjeO/L8bjNikMP0x2yg6mjBnD8qG/IT44no3Z1Sx6awdZpU3Ehvjy6FCD8Sc+pfVn66htacE7Lo6oO+4gdPEifOLju7tSERERQMFORETOJW4nHHkHdjwB1ccgJB5z7iOUhyxk/+YqSt48idvi4nj0bsoGHGbB2Ln8x8DnCfYOYd3RSpZt2k52RTND/Z284J1Dv62bcBYW0uLvT8jcuYQuXkzAuDRNtRQRkR5HwU5ERM5+jlY48GrnoihNJRA9BM/CZ8lzTiNjfTH1Jcfp8G7lcPxWOgaVc82YK5nf/+cYWPngUAUrNmeSV9XMPEcpv647SMiBXeBy4ZU2lqhbbiF47lysQZpqKSIiPZeCnYiInL1a62Dvc7D3WWhvgIRJONL/RHbVEPa/VUR7w3Ga/GvITNxI1Ehvbku9ngm9J+DymKw5WMZTm3NpKK/mmvpD/G/RHnwqy7CGhhJ63XWELV2Kb+KA7q5QRETkW1GwExGRs09DUefTuQOvgqsdBl1ES+o9HDoRzaHni3Hb86gIzuPI4E8ZOS6J3w59kOTwZOwuN3/ZU8zTm3OJzD3CsqoMRhVmYrhc+KeNJfyBewmeMweLr293VygiIvKdKNiJiMjZo/JI5/tzR94BwwKpS6lLup0D+w1OrqzC9BSSH5lF3rB9zB43hfsGP06UfxTtDjcv7SjgL58cIvXodn5Xto/oxiosISGEXnM14UuX4puc3N3ViYiI/NsU7EREpGczTSja0bnCZe568AnCnHA7pVE3kLG9mfJ1NbisDrJjdlGfksvlaQv5TdId+Hv502p38eyWXLa+s54LsrfxWMVhvNwu/EePJuzK+wiZNw+Ln193VygiInLaFOxERKRn8ng6957b/jiUZUBAFO7p/0GudTH7NlfSVFFKm08zhxM+xWd4G9eNvppp8f+FxbDQ3OHklY8PUPT635iRs4NptmrMwCAirrqSsCuX4jdwYHdXJyIi0qUU7EREpGdx2eHQm7DjSajLgfD+2Gc/ytHm6ez/uBRHcxEN/pUcStpM/7QIfjriZoZFDQOgodXOmlUf43j3HSYWZzLd48IzeBh9rr+XkIvmY/H37+biREREvh8KdiIi0jN0NEHGS7D7aWiphN6p2Oa+RFbJUA6/WY7HUUxp6AlODN/F5PGpPDH0IfoE9QGguqKGzU++Quimj5jUVIHdxw/vSy5lwI+ux2/w4G4uTERE5PunYCciIt3LVtkZ5jJeBHszDJhOzeSnyTgWSf6rtXgoJTfyAOWJh7lkQjq/THmKIJ8gTNOkbMdesp5+mdiD20l1O6mOTcS45RekXr0ES6D2nRMRkfOHgp2IiHSP2lzY+SRkvQ4eF+aQhZT0uYO9+6BqVwtOaxnHeu/CMbSCq9IuJ73f/XhZvHDbbOS/+jKVr71BeGURvbx8yR89jdTbfsj0KWndXZWIiEi36JJgZxjGPOAJwAo8b5rmI185fz9wM+ACaoAfmaZZ1BVji4jIWaZ0P+x4DLI/AKsP7pHXkxt4A3u227BtaqHVp4nDCZ8SMcbCslHXMLbXWAA6jhwh/5W/0PHJOryddurD4ji+8BZm3Xk9YxOiu7koERGR7nXawc4wDCuwEpgNlAL7DMNYa5rmsS80OwikmabZZhjG7cAfgCtPd2wRETlLmCbkbujcg65wG/iF4pj4IEddC8nYWo2juY56/wqOpmxl6IR4HhpxLwNCB+BuaaXxzbeo/MvrkHOCDqs32/qOwbxkEZdfO4e48IDurkxERKRH6IonduOBXNM08wEMw3gDWAh8HuxM09z8hfa7geu6YFwREenp3C44uroz0FUdgeBYWqf8nsz6KRz6uAqPvYrykBxyR+xl2uSxrBz8WyL9I+k4fpyKR39N49r3ob2N/JA+bBh9Gb0vW8iNc0YQE6K950RERL6oK4JdHFDyheNSYMI3tL8J+LgLxhURkZ7K0QoHX4OdK6CpGKIGUT/tOTIKB5OzuhbTU0l+ZBZVo4+xcNJcfpW0Al+3Bdu6dRS+/gbtmZm4vLzZFDuSLSlTmHTJNH41NZHIIN/urkxERKRH6opgZ5ziv5mnbGgY1wFpwPSvOb8MWAaQkJDQBbcmIiJnVGsd7H2u86u9HjN+IhWpf2T3kQgq3mrBZangePRu3CNquWr8Emb0/TGu4hIaHl1B0+rVuBsbqQvvzd+GX8qelAlcOWs4L08eQGiAd3dXJiIi0qN1RbArBfp+4TgeKP9qI8Mw0oFfAtNN07Sf6kKmaT4HPAeQlpZ2ynAoIiI9UEMR7FoJB14FVzuegReTH3kHu/ZBc4aDdq9KjvbdTmSawW1jriE1fBgtW7ZQ+t/LaN2xA9Ni5XjSaF4dNpbifkO5eVoiD0/sR7CfAp2IiMi30RXBbh+QYhjGAKAMuAq45osNDMMYDTwLzDNNs7oLxhQRkZ6g8nDn+3NHVoNh4Bx2Dcd9bmDPThv2BgdNvjVkJ21n6OR4Hkq9m95tvjS+9Ta5b9+Hq7ISd2Q0WyYu5MXQVKwxMdw6LZFrJiQQ4KPdeERERL6L0/7kNE3TZRjGXcAndG538KJpmkcNw3gIyDBNcy3wRyAIeNswDIBi0zQvPd2xRUSkG5gmFG6HHY93rnTpE0T76HvItC8kc3s9nvZmqoIKyRu2l+lTx7J88MN4H8ym4Zd/InfjRnC76Rg1jrfGLOFNr370Dg/kvhlJLE3ri5+3tburExEROSsZptkzZzympaWZGRkZ3X0bIiLyfzxuOP4BbH8cyg9AYDSNQ+9lX81kTu5rALeFwvDD1KSc4NIps5kfNZW2tR/Q+MabOAoLsYSGUj9tLk8Fp7KtzY+EiADumJHEkjHx+HhZurs6ERGRHscwjP2maaZ9m7aa6yIiIt/M2QGH3oCdy6EuF8IHUDluJbsKkil7vxWPUcfJ6H2YI+q4atJixtTPoOmlNyn66L8x7Xb8Ro6k4vYH+aM9nux6B4mBgfzvxcksHBWLl1WBTkREpCso2ImIyKl1NEHGi7D7aWipwuw9iqKxq9h+NJym953YrXVkx+0karyFO4YtJnZ3Pg33PkHx0aMYAQEEX3op+0fO5NECk8KKNgb18mX51cO4aEQfrJZTLagsIiIi/y4FOxER+bLmCtj9FGS8BA4b7v7pnEh+hp0ZYM8Em081JxJ3MuyCOP4rZAFeazfR9KvbqLDZ8E1JJvKXv2RD/BhW7qmkLKOV4XEhPHPdWOYM7YVFgU5EROR7oWAnIiKdanM6V7g89CZ4XNhTriDL+gP27+3A0wq1AaUUDMtgxqSR3FhxAR3Pvodt75/B25uQOXMIuHwpa1xRPLstn6qjBYxOCON/Fg1nxqBoPls4S0RERL4nCnYiIue7kn2dK1we/xC8fGkdegt7WxdwbFcrOB2UhJ6kbtxJlgwdy627+9N838vU19biHRdH9P33433Jpbye08LzG/OpbalhwoAIHl06islJkQp0IiIiZ4iCnYjI+cg0IWd9Z6Ar2gF+YTSM+hU7qidTuLEN09NKXuRBSK3lOt8ken3ioOV/H6HBNAmaPp3wq6/CNXYCq/aU8MLzWTS2OZmaEsXdM1MYPyCiu6sTERE57yjYiYicT9xOOPJO55TL6mMQEkflmEfZmj+Q6nUO3BYbJ6P30muEg9sqAvH68wGcpWtpj4wk8pZbCF96BS1h0Ty7o4CX//gpNruLWYNjuGtmMqMTwru7OhERkfOWgp2IyPnA3gIHV8GuldBUghk1hOLRL7DlWCQtH5nYrU2c6Lub0f0N7j/UivPhTZh2O95pacTc/2OC09OptZv8aVs+q3Yfps3hZv7w3tx5YTLD40K7uzoREZHznoKdiMi5rLUW9jwL+/4M7Q14+k7m5IA/sW2/F44jFlp8GihI3Mtcfw8LtxfiWHUYp78/oYsWEX7NNfgNGkhlUwePfZLD63uLcbg8LBgZy50XJjOwV3B3VyciIiKfUbATETkXNRTCzhVw8DVwteNMvpRDfj9i714nnmYvGvwrqU86yEJbK4vfy8JdVw/9+tHr5z8jdPFirCEhlDa08fSaw7ydUYrbNFk8Oo47ZiSRGB3U3dWJiIjIVyjYiYicSyoOdb4/d3QNGBY6Bl/HHvdlHN7TjtFhUhlUgm/fo1xUUIn/y4fA48F/xgzCr7mGwAsmY1gsFNa2svLtLNYcLMMw4Iq0vtw+PYm+EQHdXZ2IiIh8DQU7EZGznWlCwdbOFS7zNoFPELaRP2Zbw3Tyt9oxXE7KQk7QL+A4Vx44gbWwDGtoKKE/vIHwq6/GJz4egNxqGys25bI2qxxvq4XrJvbj1umJ9An17+YCRURE5F9RsBMROVt53JD9fmegKz8IgTHUp/2WjcXDqPrEhWnaqQg5yJjWbGZuyIK2dvyGDiX8N3cQcvFFWPz8ADhW3syKzTl8fKQSf28rN09N5OapA4gJ9uvmAkVEROTbUrATETnbODsg66+wcznU50NEIuVpK1l/vDctH1hwWtpp9tvHlLIs0rccA29vQubNI+Laa/AbOfLzTcOzShpZvimXDdlVBPl6cceMJG6akkhEoE83FygiIiLflYKdiMjZor0RMl6A3c9AazVmn9Hkjn2VzQf9cH7gjd3aimHZzZRjuwiqrMCrd2/C77uXsMsvxysq6vPLZBTW8+SmXLaerCHU35sfpw/kh5P7Exrg3Y3FiYiIyOlQsBMR6emay2H3U5DxMjhsuAekcyh0Gbv2mpgHfWi31hPRsZOpBz/F295GwIQJhP/iZwTPnInh1fnPvGma7Mqr48lNOezOrycy0IefzhvMdRMTCPZToBMRETnbKdiJiPRUbfWw8aHOLQtMN87BV7DbuJLMPQ4srd50WMpJrNnGwGPbsfr7EnbZIsKvvhrflJTPL2GaJp+erGH5plz2FzUQE+zLry4ewjUTEgjw0UeAiIjIuUKf6iIiPY1pwpF3YN3PoK2e9hHL2Nw6h9wdDqwOA4+nhKF5m+lTlonvgAGE/+LnhC5aiDU4+AuXMFl/rIoVm3M5VNpEbKgfDy8cxhVpffHztnZjcSIiIvJ9ULATEelJGkvgwwcg5xOao6bzSdAyKjf6YHGDb8dxUrPXE2YrJGjmhUQ8/AIBkyZ9vhgKgNtjsu5IJcs35XC80kZCRAC/v2wEi0fH4+Nl6cbCRERE5PukYCci0hN43LD3z7DxIeqcvfjI90kaj8ZimBDesI/BOesJ9u0gfOkVhF91Jd5xcV/q7nJ7eP9QOSs25ZJX00pidCCPLh3JpSNj8bIq0ImIiJzrFOxERLpb1VFYew+VhdV8YPycjtohYLqIrdxOYuF6QlNiifjPewmZPx+Lr++XujpcHtYcLOWpLXkU1bUxuHcwK64ZzfzhfbBajK8ZUERERM41CnYiIt3F2QFb/0jBlg/5pP1K3C0jMTwd9CvdQN/KT4lKn0rkw8/gn5r6T107nG7e3l/KM1vyKGtsZ0RcKM9eP5bZQ3phUaATERE573RJsDMMYx7wBGAFnjdN85GvnPcFXgXGAnXAlaZpFnbF2CIiZyOzYBvH3nySLZXToOO3WF0tDCj5gNiOTOKuXkr4Ff8Pr8jIf+rX7nDz173FPLc1j6pmO2MSwvifxcOZMTD6S+/aiYiIyPnltIOdYRhWYCUwGygF9hmGsdY0zWNfaHYT0GCaZrJhGFcBvweuPN2xRUTONmZbPftfe5zdh+IxXLfj7WiiX8lq+kTV0f/+Gwme9YfP9577oha7i9d2F/H8tnxqWxxMTIzgsaWjmJQUqUAnIiIiXfLEbjyQa5pmPoBhGG8AC4EvBruFwK8/+/lvwArDMAzTNM0uGF9EpMfzuD1sfWsVR7caYM7Av6OO+PK36DcmmP4/vQe/QQNP2a+p3ckrOwt5cUcBjW1OpqZEcc+sFMb1jzjDFYiIiEhP1hXBLg4o+cJxKTDh69qYpukyDKMJiARqu2B8EZEey+Vy8/c31lK4rQPT6EtAWyW9G95kyIJU+l71B6whIafs19Dq4MUdBby8oxCb3UX6kBjumpnCqL5hZ7gCERERORt0RbA71Rygrz6J+zZtMAxjGbAMICEh4fTvTESkm3S02/n4xdVUZnrjsUYQ1NJMlOt10m5aTEz60xiWU29BUGOz8/y2fFbtLqLN4Wb+8N7cNTOZYbGhZ7gCEREROZt0RbArBfp+4TgeKP+aNqWGYXgBoUD9Vy9kmuZzwHMAaWlpmqYpImed5tp61i1/i/qy3ri9ehHcmkekz7tMe+Amgsf8+Wv7VTZ18Myneby+txin28OCkbHceWEyA3sFn8G7FxERkbNVVwS7fUCKYRgDgDLgKuCar7RZC9wA7AIuBzbp/ToROZdUHc1my3Mf09CagttrIIHt2fQOep6Zd1+Oz8RV8DULnJTUt/HMp3m8nVGKxzRZPDqOOy5MZkBU4BmuQERERM5mpx3sPntn7i7gEzq3O3jRNM2jhmE8BGSYprkWeAFYZRhGLp1P6q463XFFRLqbaZoUfLyOfW9n0mAdidtrFH7OLAYGvMOMRYOwXPQyBEadsm9hbSsrN+ey5mAZFsPg8rR4bp+eRN+IgDNbhIiIiJwTjJ764CwtLc3MyMjo7tsQEfkn7pYWsl95jeytddSGjMVj8cLbzGJk+OtM6G2HBY9Bcvop++ZU2Vi5OZe1WeV4Wy1cPT6BW6cn0ifU/wxXISIiIj2dYRj7TdNM+zZtu2SDchGR84E9L48jf36ZgpP+VEeNwwxLxBqYzfSotxnuPAoT74ALfwE+/zyN8mh5Eys35/LxkUr8va3cMjWRm6YOICbYrxsqERERkXONgp2IyDcw3W5smzdz9JW3KW1Npjp6OkS7sMYWMS9hB/2LVkP4CLh0I8SN+af+WSWNLN+Uw4bsaoJ9vbhzRjI/mjKAiECfbqhG9osVggAAHRtJREFUREREzlUKdiIip+BpbaVh9WqO/+VDSoImURd1JWZgBz5Da1k4opFeu38DZQ5I/zVMugus3l/qv6+wnuWbctl6soawAG/unz2QGyb3J9Tf+5TjiYiIiJwOBTsRkS9wVldTu2oVuR/soSTmQhr6L8NttBE8oYXLpvcnZPOD8OkuGDANLnkcIpM+72uaJrvy6nhyUw678+uJDPThZ/MHc93EfgT56p9bERER+f7oLw0REcCek0PV8y9QsKuQovjZ2AbdhsurhZipThZdNB3//U/BazeBdwAsfApGXfP5FgamabLlZA0rNuWyv6iBmGBf/uOSoVwzPgF/H2v3FiYiIiLnBQU7ETlvmaZJ2+7dVL3wAsUn2inoP5e2oXNw+DbT/0KDS+bPx7tsB7yaDjXHYfhlMO8RCIoBwOMxWZ9dxYpNuRwuayIuzJ+HFw3nirHx+Hkr0ImIiMiZo2AnIucd0+mked06Kp9/gdKmCAr6z8E+LAZHYBODZ/sxe9Y0vPI3wKvzoCwDQvvCNW/DwDkAuD0mHx+pYMWmXI5X2ugXGcDvLxvB4tHx+HhZurk6EREROR8p2InIecNts9H41ttUrnqNcutA8vtdh6tPBPawRkbPD2XGBdMwjq+F52dA1WEIS4CLH4VR14K3Hy63h7VZ5azcnEteTStJ0YE8duVIFqTG4mVVoBMREZHuo2AnIuc8Z0UF9a+uovrtNZSFp1GQfA8erxAc0Y1MviSG8WOnYBx+C55eBnW5EDUQFj0DIy4HqzcOl4c1+4p5akseRXVtDO4dzMprxjBveG+sFqO7yxMRERFRsBORc1f70aPUv/Qyteu3UBI7naIxvwBLIK7YJi5cmMCIIRMxMv8Cy6+CphLoPQKueAWGLACLlQ6nm7f3FvLMp/mUNbaTGh/Kc9ePJX1ILywKdCIiItKDKNiJyDnF9Hho3baNupdeouFANoX90ymd+BCG4YunfxNzF6cwsH8YZLwIT6yA1mroO6FzymXKbDAM2h1u/rInn+e25lNtszO2Xzi/WTyc6QOjMQwFOhEREel5FOxE5JzgcThofv996l58iabSBnKT51A9+WrAijWllYsvG0NCjDfsfQ7WPA0djZA4A6a+CP2ngGHQYnexalcRz2/Lp67VwaTESB6/chSTkiIV6ERERKRHU7ATkbOaq6GBxjffpH7VazS3WTk+eD6Nk8ZiGuA31M6ll4+jV7ALdq2Av74AjhYYdBFM/QnEjwWgqd3JyzsKeXFHAU3tTqYNjOaemcmk9Y/o5upEREREvh0FOxE5KzmKi6l/5VUaV79DszWKw0Mvoz1gBB6Lm9BRbi5dMoVwrwbY+Ts48Cq4HTBsMUx9AHoNA6C+1cEL2/N5dWcRNruL9CG9uHtmMiP7hnVzdSIiIiLfjYKdiJxV2jMzqXvxJWwbNtAQ2p+sMTfh8R6Ky8tBzDiDBYumE+Qqg+0/haw3OjuNvAqm3A+RSQBU2zp4flsBr+0uot3pZv7w3tx1YQpDY0O6sTIRERGRf5+CnYj0eKbbjW3TJupfepm2Aweo6T2cQ1PuxWJJxunTQcJkLy5aMAW/lhxYfwccXQNWH0i7CSbfDWF9AahoaufZT/N5fW8xTreHS0fGcueFyaT0Cu7mCkVEREROj4KdiPRYpstF03trqX3uWRxFxZQljuPIhQ/iY/bD7ddO4nR/Zs+fhk9tFqy9AU58BD5BMPkemHQnBMUAUFLfxtOf5vG3jFI8psmSMXHcMSOZ/lGB3VyhiIiISNdQsBORHsd0uWj64ANqn34aR1EJBUOmcnzW9fi5YzEC2hicHsb0WdPwKtsFb10G+VvAPxxm/AImLOv8GSiobWXl5lzWHCzDahhckRbPbdOT6BsR0L0FioiIiHQxBTsR6TFMt5vmjz6iduVTdBQVkzdkOjkzb8TPE4NXUBuj5sUwaeoQLPkbYdVFULIHgnrB7Ich7Ubw7ZxSebLKxsrNubyfVY631cIPJvXj1mlJ9A716+YKRURERL4fCnYi0u1MjwfbunXUrFhJe2ExeYNnkHfhj/D1ROId1srES+IZMz4J48T78PztUHkIQvvCRX+C0deDd2dgO1rexIpNuXx8pJIAHyu3TEvk5imJRAf7dnOFIiIiIt8vBTsR6Tamx4Pt7+upXbmCtrxicgbPpGj6zXiboXhHtTL10v6MGBmHcXQ1PHMt1J6EyGRY+BSkLgWrNwCZJY2s2JTDhuxqgn29uHtmMj+6YADhgT7dXKGIiIjImaFgJyJnnGma2DZsoHbFSlpzi8gZlE7J9JvxMoPx7tNK+uIUBvb3wdj/MjzxIrRUQq8RcPlLMHQhWKwA7Cus58mNOWzLqSUswJsHZg/kB5P7E+rv3b0FioiIiJxhpxXsDMOIAN4E+gOFwFLTNBu+0mYU8DQQAriB35im+ebpjCsiZyfTNGnZvJma5StoySnm5KB0yqfejJUAfPq2MWfJYJKCK2DPw7BmNXickDwbJqyA5HQwDEzTZGduLU9uzGFPQT1RQT78bP5grpvYjyBf/b8qEREROT+d7l9BPwM2mqb5iGEYP/vs+KdfadMG/MA0zRzDMGKB/YZhfGKaZuNpji0iZwnTNGn59FNqV6yk6WQxOSmzqZxyCxZ88U1qY/6iwSS0b4dt10NZBvgEw7ibYNwtEJX8+TW2HK9m+aYcDhQ30ivEl/+8ZChXj0/A38fazRWKiIiIdK/TDXYLgRmf/fwKsIWvBDvTNE9+4edywzCqgWhAwU7kHGeaJq3bt1OzfDmNJ0o4mTKX6sk3YeCF/yAHF8/vS5+Kv8F790NLVef7c/P/CKOu/nyFS4/HZH12FSs25XK4rIm4MH/+Z9FwLh8bj5+3Ap2IiIgInH6w62WaZgWAaZoVhmHEfFNjwzDGAz5A3mmOKyI9mGmatO7cSc3y5dQfLyMnZS61k27GNAyChrlZMCWAqLw34K13O6dbpsyBCbdC4kywWABwe0w+OlzBys25HK+00S8ygD9clsriMXF4Wy3dXKGIiIhIz/Ivg51hGBuA3qc49cvvMpBhGH2AVcANpml6vqbNMmAZQEJCwne5vIj0EK2791D95BPUHa8gJ3kedRNvxmPxEJbqYcGwJsKzn4F3D4BvCIy/BcbdDJFJn/d3uT28l1nOyi255Ne0khwTxGNXjmRBaixeCnQiIiIip2SYpvnvdzaME8CMz57W9QG2mKY56BTtQuicpvk70zTf/jbXTktLMzMyMv7texORM6tt3z6qnnySmuxKcpLm0xCeisfqJmqkiwVxRwk+9hy01kDUQBi/DEZe9fl0SwCHy8PqA6U8tSWP4vo2hvQJ4e6Zycwb1huLxejGykRERES6h2EY+03TTPs2bU93KuZa4Abgkc++v3eKm/EB1gCvfttQJyJnj7YDB6h68gmqjtWQmziPprThuLwc9BnRxCVBGwjIewvK3TBw7mfTLS8E4x9BrcPp5q2MEp7Zkkd5Uwep8aH8xyVppA+JwTAU6ERERES+jdMNdo8AbxmGcRNQDFwBYBhGGnCbaZo3A0uBaUCkYRg//KzfD03TzDzNsUWkG7VnZlLx5BNUHqsjN3E+tjGDcHl3kDC0mPnmX/Gr2AO+oTD+Vhh/M0Qkfql/m8PFX/cU89zWfKptdtL6hfO7y1KZlhKlQCciIiLyHZ3WVMzvk6ZiivRM7YcPU/7EY1QcayJvwDxagpNw+baRmJzP3LZn8WkvhahBMGEZpF4FvkFf6m/rcLJqdxEvbCugrtXBpMRI7p6VzKTESAU6ERERkS84k1MxReQ80X70KGVPPEb5sRbyBsyjLTUBp5+NIX23MavlabzqOmDQfBi/AhJnfGm6JUBTm5OXdhbw0o5CmtqdTB8Yzd0zk0nrH9Et9YiIiIicSxTsROQb2fPzKfrDI5Qda6eg/1zah8fi9G8kNXo10+2vY3UEwsSbOle3jBjwT/3rWx28sD2fV3cWYbO7mD20F3fPTCY1PqwbqhERERE5NynYicgpuerrOf6nRyjb20xJ3CzsQ6NxBdQwNvjPTGYdltCBMOGPkHol+AT+U/9qWwd/3prPa7uL6XC5uWh4H+6amcyQPiHdUI2IiIjIuU3BTkS+xN3RwZ7H/kjVrg6qYqbhTg7AE1DKeP//Zbz3LozB82H8uzBg2j9NtwQob2znua35vL63GKfbw8JRcdx5YRLJMcGnGE1EREREuoKCnYgAYHd2sH7FSpp2Q1PYRMw+4O1/gPSgd0gOqYNRV8O4lRDe/5T9S+rbeGpLHn/bX4JpwmVj4rl9RhL9o/75aZ6IiIiIdC0FO5HzXH1rA+8/9wrOzBDs/qOxhLQR5r2OebHbiEidAsP+BP0uAIv1lP3za1pYuTmPdzPLsBoGV47ry23Tk4gPDzjDlYiIiIicvxTsRM5TOSVH+WjVWrwLBmJaU/E1qujHG8yeF4zvmEXQ7w9fG+YATlTaWLk5lw8OlePjZeGGSf25dXoivUL8zmAVIiIiIgIKdiLnFbPDxrYtr7BjezPBNaPxMiYQYjtB/6hNTHpwMdbBT39jmAM4UtbEik25rDtaSaCPlWXTkrh56gCignzPUBUiIiIi8lUKdiLnOnsLjuMfsnb7FrILBxNhG0Gox0mv6gySE+oY9uRP8Yq5/V9e5mBxAys25bLxeDXBfl7cMzOZGy8YQHigzxkoQkRERES+iYKdyLnIboOTn9BwaA1/y7VQ0ziH0PYriHHaiC/7kOS+dvr/6T58k5P/5aX2FtSzfFMO23JqCQvw5idzBvKDyf0J8fM+A4WIiIiIyLehYCdyrrDb4MQ6OPYuubkZvN8xF3fjFfi7QoixlzOgYBUJYc3E/ucDBE6e/I2XMk2THbl1PLkph70F9UQF+fDz+YO5bmI/An31z4aIiIhIT6O/0ETOZh3NcHIdHH0Xcjeww5PANvsSAhpuwMf0IsB5koFHXyLap5GY++4jdOGlGNavf4fONE02n6hm+aZcDhY30jvEj/9aMJSrxyfg5/3N796JiIiISPdRsBM523Q0w4mP4di7kLsRl8vBB75zONbyO0KbBuBvcRBuyWbonvcI9DQRedOPiLzxRiwBX7/9gMdj8vdjVazYnMORsmbiwvz5zeLhXD42Hl8vBToRERGRnk7BTuRs0FYPOevh6BrI2whuB02BA3gn6C5qi4cS2B6Bj28zUeFHGLb5HazNtYRddhlRd9+Fd0zM117W7TH58HAFKzflcqLKRv/IAP5weSqLR8fhbbWcwQJFRERE5HQo2In0NB1NUJEF5Qc7v8oOQGNR57mQOAoG38z7ZUm4cvvg6/LHE1ZNQt98Bn38Nu6SYgKnTiXmJ3/Gb9DArx3C6fbwXmY5T23OJb+2lZSYIJ64ahQXj+iDlwKdiIiIyFlHwU6kOzlaoeLQP0Jc+UGoy/nH+bAEiB0NaTeyx96frXva8dsYjQE4+1YxfqCD+Hfeo/3dTLwGDiT2+ecJmnLB1w/n8vDOgVKe2pJLSX07Q/qE8NS1Y5g3rDcWi/H91ysiIiIi3wsFO5EzxdkBVUeh/MA/QlzNcTA9neeDYztDXOqVmLGjqbD0Z9+JUgpPVGPfYyWwNRQvqw+uEdXMT0sm6I2PsK1ahyM6ij7/8zChixd/7cIoHU43b+4r4ZlP86ho6mBkfCj/dckwZg2JwTAU6ERERETOdgp2It8HtxOqj/1jKmX5wc5jj6vzfEAUxI2BIQsgdjSeXqMorjM4kHWc0t11uMrb8bUXd17KGoAzqh7fNFg6bRye196g4ebf0+LlRdSddxL5oxuxBAae8jbaHC7+sruY57blU2Ozk9YvnEcuS2VaSpQCnYiIiMg5RMFO5HR53FBz4gvTKQ9A5RFw2zvP+4V2PombfE/n99jROP36UJhbTeahk1Rub8RTmYWX2weANl8TR3QVvgMCGD48kfFDJuGHlYY33qB2yVLczc2ELllM9D334N2r1ylvydbh5NVdRbywvYD6VgeTkyJ58qrRTEyMUKATEREROQcp2Il8Fx4P1Od9+Z24iixwtnWe9wmCPqNg/C2dT+RiR0P4ANpsTgpOVHJ4Tx41+ZkYtScwTAsmHhoCbbgSyumVFEzq8BTGRA+D3CLsJ07Q8db7VJ74X+w5OZh2O4GTJxHz4IP4DR58yttranPy4o4CXtpRQHOHixmDorl7ZjJj+0WcwV+SiIiIiJxphmma3X0Pp5SWlmZmZGR0923I+cbeAi1VX/iq7vxuq+pcmbIiC+zNnW29/KFP6udP4YgdDZEpmIZBU3U7BceryD5WSF1BO5ZmPwBchoOa4BLM3q3EJoYyNiKc5Ebw5BbQceI49hMncVVWfn471ogI/AYPwnfgIAKnTiFw8uRTPnGra7Hz/PYCVu0qosXuYs7QXtw1M5nU+LAz8msTERERka5nGMZ+0zTTvk1bPbGTc5/bBa01XwlsXwlt/3fsbP3n/oYVgmIgJBZGXA6xnz2Jix4MVi/cbg+1xS0U7q/iZPanNBU5MDq8AWj3aqE6pBCflFZSAgyGO1zEVLTg2p6D/eVcTIeDGgAvL3wTEwkYN+7zIOc3eBDWqG9+F666uYPntubzlz3FdLjcXDSiD3ddmMyQPiHfz+9SRERERHqk0wp2hmFEAG8C/YFCYKlpmg1f0zYEyAbWmP+/vbuLkao+4zj+/c3smyIgsoLI2wrCWtO0tVKLFw2oJTGaaFuaVmKN9EVrFb2obaSxF01NE9Ok8crYUNPUNKnEelGstZqIYtIGEzGgVqvIi+CKLrCFNejK7sw8vZgDO7DL7oGdl52d3yfZzDlz/vM/z+6z/5l99pzzPxFrxrJfMyKK93s7VpyNVLB92gMMc2S6bSqcM7P4NfvyZHnG4OPkC4rLZ50HmcF7u/X35fhody97Nu9g59sfcuSDHMoVZ6PsbT3AkdbdTD/rMIv6P6bjYA9tb3xEvnv/8df3TZ9OW2cn026+uVjEdXbSumABamlJ/e3vO9zH71/ayfpX3ieXL3Djl2Zz11ULuXjG5DP+kZqZmZlZ/RrrEbu1wMaIeFDS2mT9vlO0fQB4aYz7s3oRAbmjkOtLHj8rTvefK/06CgN96dsN9MGnBwcLtmOTk5TKtgwWa9M6YO4VQwq1mHQ++ZZ2crSQ6y+Q68+TGzjpsbdA7mCeXP9Rcv1d5AYK9B76hPe2d9PXHShEgQKfNHeh/C4u/KybjgN7mb63Cw0kM182N9O6YAFtS5fSuriT1ks6aevspKm9/Yx/rHt7PuWRl3bw5KtdRMDKL8/hJ8sX0tE+/KyYZmZmZtYYxlrY3QgsT5YfAzYxTGEn6XJgJvAskOoc0fHo7S2b2brh2ZOePcU1iqmuXRytzdDtSp5TRLI5StrF8e1EJDEEgyfyRUlcpW0BCkm/x/ZR0m8Mvk4ERB5FAUUeIlCyTuTJEBCF4vqITjy9cOh3miGUIZQFZSgoS6iVgcwicvoiuczZ5Ggjn2klH60UaKFAM4VogsMZKGSJKD4SWVRoghCKHohexM5R4huqwFEK+fc458gu5nXv4oLunTQlxWX+3PMY6FjIkeuXMjB/AQPzFzIwex40N5/Yyf487O8+/X1H8NybH7Fh2z6yEt/9ylzuWLaQOdPOPu2+zMzMzGziGWthNzMiPgSIiA8lzTi5gaQM8DvgFuCakTqTdDtwO8C8efPGGFr5/WfjJv53YFmtw7BSUSCb7ydbGKCl0E8mP0C20E+m0Ecmef74Y6GfbL74WLp8QpvjbYe+TpEjl8mwd/JM3p5yIf+8ZAW7pl7I7qmz6G0tOQVyD7DnEDDsWclnrK05w61XdvDjZQuYOaWtrH2bmZmZWX0btbCT9DxwwTCb7k+5jzuBZyLi/dHunxUR64B1UJwVM2X/VfPVlSt5/cUXiisjfS9Dtg3XViNsO8Vrj78kUzwOJwhlgOL1X8UfWAaUKbYVQHbwxRW6fdmQbsu5n5P6yjZlyDZlaMlmaWrOks1kacq0kdUkMiquZ8gM39cYBBCtrWRmz6WjqYkO4Kqy72Vks889i2mT0l+HZ2ZmZmaNY9TCLiK+fqptkrolzUqO1s0C9g/T7Erga5LuBM4BWiQdiYi1Zxx1jcy9eDFzL15c6zDMzMzMzMxOMNZTMZ8CbgUeTB43nNwgIm4+tixpNbCkHos6MzMzMzOz8Wqs56w9CKyQ9C6wIllH0hJJj441ODMzMzMzMxudItXsjdW3ZMmS2LJlS63DMDMzMzMzqwlJr0ZEqrsKlH+WCTMzMzMzM6sqF3ZmZmZmZmZ1zoWdmZmZmZlZnXNhZ2ZmZmZmVudc2JmZmZmZmdU5F3ZmZmZmZmZ1zoWdmZmZmZlZnRu397GTdADYU+s4xqF24GCtg7CacO4bl3PfuJz7xuS8Ny7nvnGdKvfzI+L8NB2M28LOhidpS9qbFNrE4tw3Lue+cTn3jcl5b1zOfeMqR+59KqaZmZmZmVmdc2FnZmZmZmZW51zY1Z91tQ7Aasa5b1zOfeNy7huT8964nPvGNebc+xo7MzMzMzOzOucjdmZmZmZmZnXOhZ2ZmZmZmVmdc2E3Tkm6VtI7knZIWjvM9jskvSFpm6R/Sbq0FnFa+Y2W+5J235YUkjwt8gSRYtyvlnQgGffbJP2oFnFaeaUZ85K+I+ktSW9K+ku1Y7TKSDHmHyoZ79slHa5FnFZ+KXI/T9KLkrZKel3SdbWI08ovRe7nS9qY5H2TpDmp+/Y1duOPpCywHVgBdAGvAKsi4q2SNlMi4uNk+Qbgzoi4thbxWvmkyX3SbjLwD6AFWBMRW6odq5VXynG/GlgSEWtqEqSVXcq8LwKeAK6OiEOSZkTE/poEbGWT9v2+pP3dwGUR8YPqRWmVkHLcrwO2RsQjyT/vn4mIjlrEa+WTMvd/BZ6OiMckXQ18PyJuSdO/j9iNT1cAOyJiV0T0A+uBG0sbHCvqEpMAV+gTw6i5TzwA/Bb4rJrBWUWlzb1NLGnyfhvwcEQcAnBRN2Gc7phfBTxelcis0tLkPoApyfJUYF8V47PKSZP7S4GNyfKLw2w/JRd249Ns4P2S9a7kuRNIukvSTop/4N9TpdisskbNvaTLgLkR8XQ1A7OKSzXugZXJ6RlPSppbndCsgtLkfTGwWNK/Jb0syWdnTAxpxzyS5gMXAS9UIS6rvDS5/xXwPUldwDPA3dUJzSosTe5fA1Ymy98EJkuanqZzF3bjk4Z5bsgRuYh4OCIWAvcBv6x4VFYNI+ZeUgZ4CLi3ahFZtaQZ938HOiLiC8DzwGMVj8oqLU3em4BFwHKKR20elXRuheOyykv1WZ+4CXgyIvIVjMeqJ03uVwF/iog5wHXAn5O/Aay+pcn9z4BlkrYCy4APgFyazv0LMj51AaX/iZ/DyIfg1wPfqGhEVi2j5X4y8Hlgk6T3gKXAU55AZUIYddxHRE9EHE1W/wBcXqXYrHLSvN93ARsiYiAidgPvUCz0rL6dzmf9Tfg0zIkkTe5/SPHaWiJiM9AGtFclOqukNJ/1+yLiWxFxGXB/8lxvms5d2I1PrwCLJF0kqYXiG/pTpQ2Si+mPuR54t4rxWeWMmPuI6I2I9ojoSC6ifhm4wZOnTAhpxv2sktUbgP9WMT6rjFHzDvwNuApAUjvFUzN3VTVKq4Q0uUdSJzAN2Fzl+Kxy0uR+L3ANgKTPUSzsDlQ1SquENJ/17SVHZ38B/DFt5y7sxqGIyAFrgOco/uH2RES8KenXyQyYAGuSaa+3AT8Fbq1RuFZGKXNvE1DK3N+TjPvXKF5Xu7o20Vq5pMz7c0CPpLcoXkj/84joqU3EVi6n8X6/ClgfnsZ8wkiZ+3uB25L3+8eB1f4dqH8pc78ceEfSdmAm8Ju0/ft2B2ZmZmZmZnXOR+zMzMzMzMzqnAs7MzMzMzOzOufCzszMzMzMrM65sDMzMzMzM6tzLuzMzMzMzMzqnAs7MzMzMzOzOufCzszMzMzMrM79H4EutrwAiYHdAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "#xét bài toán sóng giãn\n", "uL = -0.4 # <0 - dịch chuyển sang trái\n", "uR = 1.2 # >0 - dịch chuyển sang phải \n", "\n", "u0 = np.zeros(nx)\n", "i = int(0.5/dx)\n", "u0[:i] = uL\n", "u0[i:] = uR\n", "\n", "u_mm = MUSCL_scheme(minmod, nt)\n", "u_kr = MUSCL_scheme(koren, nt)\n", "u_sp = MUSCL_scheme(superbee, nt)\n", "u_vl = MUSCL_scheme(vanleer, nt)\n", "\n", "u_exact = u0.copy()\n", "t = nt*dt\n", "i = int((0.5 + uL*t)/dx)\n", "j = int((0.5 + uR*t)/dx)\n", "for k in range(i+1, j): u_exact[k] = (x[k]-0.5)/t\n", " \n", "fig = plt.figure(figsize=[15, 5])\n", "i = int(0.3/dx)\n", "j = int(0.9/dx)\n", "plt.plot(x[i:j], u_exact[i:j], x[i:j], u_mm[i:j], x[i:j], u_kr[i:j], x[i:j], u_sp[i:j], x[i:j], u_vl[i:j])\n", "plt.legend(['exact solution', 'minmod', 'koren', 'superbee', 'vanleer'])\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So sánh với sơ đồ Godunov bậc 1 ở bài 12, ta thấy sơ đồ bậc 2 MUSCL chính xác hơn. Đối với các giới hạn độ dốc, nếu như bài toán sóng xung kích superbee cho kết quả tốt nhất thì với bài toán sóng giãn cho kết quả kém nhất, ngược lại với minmod. Vanleer có vẻ như phù hợp nếu dùng cho cả hai trường hợp.\n", "\n", "**Bài tập**: tìm hiểu về `tái cấu trúc nghiệm parabolic từng mảnh và tái cấu trúc hàm dòng`." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "#lưu giá trị u_g vào file để so sánh ở bài 14\n", "np.savetxt('data/burgers_muscl_test2.dat', u_mm)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### [Bài 14. Sơ đồ WENO5, phương pháp Runge-Kutta](Bai_14.ipynb)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.15rc1" } }, "nbformat": 4, "nbformat_minor": 2 }