chapter_5_numerical_integration/numerical_integration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Numerical integration\n",
    "\n",
    "by Xiaofeng Liu, Ph.D., P.E.\n",
    "Associate Professor\n",
    "\n",
    "Department of Civil and Environmental Engineering  \n",
    "Institute of Computational and Data Sciences  \n",
    "Penn State University \n",
    "\n",
    "223B Sackett Building, University Park, PA 16802\n",
    "\n",
    "Web: http://water.engr.psu.edu/liu/\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical integration\n",
    "\n",
    "\n",
    "### Mathematical definitions\n",
    "Integration is the computation of an integral. For example, a function $f(x)$ can be integrated over an interval $x \\in [a,b]$ as\n",
    " \\begin{equation}\n",
    "    I = \\int_a^b f(x) dx\n",
    " \\end{equation}\n",
    "where $f(x)$ is called the integrand. Graphically, the integration of a function $f(x)$ over $[a,b]$ is equivalent of the \"signed\" area bounded by the curve $f(x)$, the $x$-axis, and the two vertical bounding lines. \n",
    "\n",
    "<img src=\"integration_area.png\" width=\"500\"/>\n",
    "<h3 align=\"center\">Figure. Integration over an interval.</h3> \n",
    "\n",
    "Integration can also be over a surface, which gives a volume.\n",
    " \\begin{equation}\n",
    "    I = \\int_a^b \\int_c^d f(x,y) dx dy\n",
    " \\end{equation}\n",
    "where $f(x,y)$ depends on both $x$ and $y$. Graphically, the integration of a function $f(x,y)$ over $[a,b] \\times [c,d]$ is equivalent of the \"signed\" volume bounded by the curved surface $f(x,y)$, the $x-y$-plane, and the four vertical bounding surfaces. \n",
    "\n",
    "<img src=\"integration_volume.png\" width=\"400\"/>\n",
    "<h3 align=\"center\">Figure. Integration over a surface.</h3> \n",
    "\n",
    "Integration is the reverse of differentiation. If\n",
    " \\begin{equation}\n",
    "    F(x) = \\int  f(x) dx\n",
    " \\end{equation}\n",
    "then\n",
    " \\begin{equation}\n",
    "    f(x) = \\frac{d F(x)}{dx}\n",
    " \\end{equation}\n",
    "Here, $F(x)$ is called an antiderivative of $f(x)$. If $f(x)$ is a real-valued continuous function on a closed interval $[a,b]$, then \n",
    "  \\begin{equation}\n",
    "    I = \\int_a^b f(x) dx = \\left.F(x) \\right\\vert_{a}^{b} = F(b) - F(a)\n",
    " \\end{equation}\n",
    " \n",
    "Due to historical reasons, integration is also termed quadrature. Mathematically, there are two categories of integrals, i.e., definite and indefinite integrals in the forms of \n",
    " \\begin{equation}\n",
    " \\int_a^b  f(x) dx    \\quad \\mathrm{vs.   } \\int  f(x) dx\n",
    " \\end{equation}\n",
    "In this chapter, we are only interested in definite integrals.\n",
    "\n",
    "\n",
    "### Application examples of numerical integration in CEE\n",
    "\n",
    "In CEE, numerical integration is often used. The following is a short incomplete list of examples:\n",
    "- discharge across a surface\n",
    "- length of a curve\n",
    "- surface area\n",
    "- volume calculation\n",
    "- mass of a 2D plate\n",
    "- fluid force on a 2D surface\n",
    "- center of mass and moment of inertia\n",
    "- acceleration, velocity, and distance \n",
    "- ...\n",
    "\n",
    "Next, some of the examples above are described in more details. \n",
    "\n",
    "#### Example: Total discharge at a river cross-section\n",
    "River discharge $Q$ at one particular cross-section is the volume of water flows passing the cross-section per unit time ($m^3/s$ or $ft^3/s$). To get $Q$, we need to know the flow velocity distribution across the cross section, which is usually obtained from measurements or computer simulations. The calculation of flow discharge is through the integration of velocity over the cross-sectional area.\n",
    "\\begin{equation}\n",
    "Q = \\int_S U ds\n",
    "\\end{equation}\n",
    "where $U$ is flow velocity and $S$ is the river cross-section. Flowing water in rivers also carries dissolved and particulate matters, such as oxygen and sediment. It is often of interest to their fluxes too. To do that, we also need to know the concentration distribution $C$. With that, the flux is calcualted as another surface integration\n",
    "\\begin{equation}\n",
    "Q = \\int_S U C ds\n",
    "\\end{equation}\n",
    "\n",
    "<img src=\"river_discharge.png\" width=\"200\"/>\n",
    "<h3 align=\"center\">Figure. Integration over a river cross section to get flow discharge.</h3> \n",
    "\n",
    "#### Example: Hydrostatic/hydrodynamic force on a structure or wind pressure force on a building \n",
    "To design structures like dams, bridges, or buidlings, it is important to calculate the load of water or wind on them. Regardless it is the hydrodstatic force on a dam or dynamic force on a bridge or wind load on a building, we need to know the pressure (normal) and shear stress (tangential) distribution on the surface of a structure. \n",
    "\n",
    "<img src=\"wind_wave_CFD_results.png\" width=\"600\"/>\n",
    "<h3 align=\"center\">Figure. Storm surge and wind loading on a simulated house: (a) and (b) are computer rendering of the wave impact on the house, (c) shows the pressure distribution on the building envelope, and (d) shows the pressure force history in three different directions (note the house is a scaled down model home).</h3> \n",
    "\n",
    "Typically, the pressure contribution is dominating. Then, the pressure force can be calculated with the integration as follows:\n",
    "\\begin{equation}\n",
    "\\mathbf{F} = \\int_\\mathbf{S} p d\\mathbf{s}\n",
    "\\end{equation}\n",
    "where $p$ is pressure and $\\mathbf{S}$ is area vector (magnitude = area, direction = normal).\n",
    "\n",
    "#### Example: Length of a curved road (arc)\n",
    "When we design a curved road or a highway, which is an arc, we need to calculate the total length. If we know the mathematical expression for the curve, we can calculate the length through integration as follows:\n",
    "\\begin{equation}\n",
    "L = \\int ds = \\int \\sqrt{\\Delta x^2 + \\Delta y^2} =\\int \\sqrt{1 + \\left(\\frac{\\Delta y}{\\Delta x}\\right)^2} \\Delta x = \\int_a^b \\sqrt{1+f'(x)^2} dx\n",
    "\\end{equation}\n",
    "Here $L$ is the length of the arc and slope $f'(x) \\approx \\frac{\\Delta y}{\\Delta x}$.\n",
    "\n",
    "<img src=\"arc_length.png\" width=\"400\"/>\n",
    "<h3 align=\"center\">Figure. Length of a curved road.</h3> \n",
    "\n",
    "#### Example: Center of mass\n",
    "For structural and other analyses, the center of mass is an important property. For very simple shapes and density distributions, it is easy to calculate the center of mass either by hand or through a table. However, for more complicated geometry and density distribution, numerical integration is needed. \n",
    "\n",
    "For example, a beam is located between $x\\in[a,b]$, with uniform cross-sectional shape, but nonuniform density $\\rho(x)$. Its center of mass along the $x$-axis can be calculated as \n",
    "\\begin{equation}\n",
    "x_{\\textrm{cm}} = \\frac{\\int_a^b x \\rho(x) dx}{\\int_a^b \\rho(x)dx}\n",
    "\\end{equation}\n",
    "where the numerator is the moment around $x$ = 0 and the denominator is the total mass of the beam.\n",
    " \n",
    "<img src=\"center_of_mass_beam_nonuniform_density.png\" width=\"400\"/>\n",
    "<h3 align=\"center\">Figure. Center of amss for a beam with nonuniform density.</h3> \n",
    "\n",
    "Another example is a 2D plate with uniform density. The plate is formed by two curves $f(x)$ and $g(x)$ and the two curves intersect at $x = a$ and $x = b$. \n",
    "\n",
    "\\begin{equation}\n",
    "x_{\\textrm{cm}} = \\frac{M_y}{M} = \\frac{\\int_a^b x [f(x)-g(x)] dx}{\\int_a^b [f(x)-g(x)]dx}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "y_{\\textrm{cm}} = \\frac{M_x}{M} = \\frac{\\int_a^b [f(x)^2-g(x)^2] dx}{2\\int_a^b [f(x)-g(x)]dx}\n",
    "\\end{equation}\n",
    " \n",
    "<img src=\"center_of_mass_2D_plate.png\" width=\"400\"/>\n",
    "<h3 align=\"center\">Figure. Center of amss for a 2D plate with uniform density.</h3> \n",
    "\n",
    "#### Example: Deflection of beams\n",
    "\n",
    "Beams will deflect under loading. In any structural analysis textbook, there are tables for the deflection formulas if the loading regular, such as uniform or linear. For more complicated loading, the deflection needs to be solved with computer. In this example, a nonuniform load $w(x)$ is applied to a beam and induces nonuniform moment $M(x)$. The deflection $y$ is governed by the moment equation as\n",
    "\\begin{equation}\n",
    "\\frac{d^2y}{dx^2} = \\frac{1}{EI(x)} M(x)\n",
    "\\end{equation}\n",
    "where $E$ is the Elastic or Young’s modulus of the beam material, $I(x)$ is second Moment of Area of the beam's cross-section, and $M(x)$ is the bending moment at a section $x$ from the origin. Integrate the equation once will give the slope along the beam:\n",
    "\\begin{equation}\n",
    "\\frac{dy}{dx} = \\int \\frac{1}{EI(x)} M(x) dx + C_1\n",
    "\\end{equation}\n",
    "Integrate the equation again will give the deflection along the beam:\n",
    "\\begin{equation}\n",
    "  y(x) = \\int \\left[ \\int \\frac{1}{EI(x)} M(x) dx \\right] dx + C_1x+C_2\n",
    "\\end{equation}\n",
    "where $C_1$ and $C_2$ are integration coefficients depending on boundary conditions. \n",
    "\n",
    "<img src=\"beam_deflection.png\" width=\"400\"/>\n",
    "<h3 align=\"center\">Figure. Beam deflection.</h3> \n",
    "\n",
    "\n",
    "## Numerical integration\n",
    "Numerical integration is the use of computer to calculate the approximation of an integration. The basic idea of numerical integration is to divide the integration interval into small pieces and then do summation. For example, for the following integration \n",
    "\\begin{equation}\n",
    "    I = \\int_a^b f(x) dx\n",
    "\\end{equation}\n",
    "the integration region $x\\in[a,b]$ can be divived into many subregions, whose areas are approximated with different methods. All the different methods differ in how the subregion areas are calculated. It is also good to note that the integrand $f(x)$ can be in the form of tabulated data or a function.\n",
    "\n",
    "<img src=\"numerical_integration_idea.png\" width=\"300\"/>\n",
    "<h3 align=\"center\">Figure. Idea of numerical integration.</h3> \n",
    "\n",
    "\n",
    "In this chapter, two categories of numerical integration methods are introduced. They are very popular in practice. \n",
    "- Newton-Cotes methods\n",
    "- Gauss quadrature\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Newton-Cotes methods\n",
    "\n",
    "In Newton-Cotes methods, the integrand is replaced with a simple polynomial which is easy to be integrated. It can be written as the following general equation:\n",
    " \\begin{equation}\n",
    "    I = \\int_a^b f(x) dx \\approx \\int_a^b f_n(x) dx\n",
    " \\end{equation}\n",
    " where \n",
    " \\begin{equation}\n",
    " f_n(x) =a_0 + a_1x + \\ldots + a_{n-1}x^{n-1} + a_nx^n\n",
    " \\end{equation}\n",
    " is a $n$-th order polynomial. The polynomial can be for the whole data points or as piecewise segments for each subregion. For example, the following figure shows that the integrant is approximated as a piecewise linear function.\n",
    "\n",
    "<img src=\"newton_cotes_linear.png\" width=\"300\"/>\n",
    "<h3 align=\"center\">Figure. Newton-Cotes method where the integrant is approximated as a piecewise linear function.</h3> \n",
    "\n",
    "Three popular Newton-Cotes methods will be introduced here:\n",
    "- trapezoidal rule\n",
    "- Simpson's 1/3 rule\n",
    "- Simpson's 3/8 rule"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Trapezoidal rule\n",
    "The trapezoidal rule is defined when the approximation polynomial is first order, i.e., the integrand curve on an interval is approximated as a linear line. The approximated area under the line is a trapezoid (thus the name). On an interval $[x_i,x_{i+1}]$, the integral (area) is\n",
    "\\begin{equation}\n",
    "      I_i = (x_{i+1}-x_i)\\frac{f(x_i)+f(x_{i+1})}{2}\n",
    "\\end{equation}\n",
    "\n",
    "<img src=\"trapezoidal_rule.png\" width=\"300\"/>\n",
    "<h3 align=\"center\">Figure. Trapezoidal rule.</h3> \n",
    "\n",
    "If the integral interval $[a,b]$ is divided into $n$ segments of equal size $h$, the points on the $x$ axis can be denoted as $(x_0,x_1,\\ldots,x_n)$. Then,\n",
    "   \\begin{equation}\n",
    "      h = \\frac{b-a}{n}\n",
    "   \\end{equation}\n",
    "The total integral (area) is the summation of small areas, which can be expressed as\n",
    "  \\begin{eqnarray}\n",
    "     I &=& h\\frac{f(x_0)+f(x_1)}{2} + h\\frac{f(x_1)+f(x_2)}{2} + \\ldots + h\\frac{f(x_{n-1})+f(x_n)}{2} \\\\\n",
    "     & = & \\frac{h}{2}\\left[f(x_0) + 2 \\sum_{i=1}^{n-1} f(x_i) + f(x_n) \\right]\n",
    "\\end{eqnarray}     \n",
    "   \n",
    "<img src=\"trapezoidal_rule_all.png\" width=\"300\"/>\n",
    "<h3 align=\"center\">Figure. Trapezoidal rule with the summation of all small areas.</h3> \n",
    "\n",
    "To increase the accuracy of the integration with the Trapezoidal rule, you can increase the number of segments $n$. It is equivalent to reduce the segment size $h$. Indeed, the error of the Trapezodal rule follows\n",
    "\\begin{equation}\n",
    "      E \\propto -\\frac{(b-a)^3}{12n^2}\n",
    "\\end{equation}\n",
    "When $n$ doubles, the error is reduced to a quarter. This behavior is called \\bverb|second| order convergence.\n",
    "\n",
    "Next, let's do a manual calculatio example. Using the trapezoidal rule, calculate the following integration \n",
    "\\begin{equation}\n",
    "I = \\int_0^1 (x^2+1) dx\n",
    "\\end{equation}\n",
    "with \n",
    "- only one segment \n",
    "- two segments\n",
    "- three segments\n",
    "(Note the true solution is $I=\\frac{4}{3}$)\n",
    "\n",
    "The following is an implementation of the Trapezoidal rule. It is used the solve the same example integration as above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trapezoidal rule integration result =  1.3333335000000004\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from math import exp  \n",
    "import sys\n",
    "\n",
    "#define the trapezoial rule: f is the integrand function\n",
    "#a and b are the bound so of integration\n",
    "#n is the number of segments \n",
    "def trapezoidal_rule(f, a, b, n):\n",
    "    h = (b-a)/n\n",
    "    s = (f(a) + f(b))\n",
    "    for i in range(1,n):\n",
    "        s = s + 2.0*f(a + i*h)\n",
    "    \n",
    "    return s*h/2.0\n",
    "\n",
    "#the integrand function\n",
    "def func(x):\n",
    "    return x**2+1.0\n",
    "\n",
    "a = 0\n",
    "b = 1.0\n",
    "n = 1000\n",
    "\n",
    "I = trapezoidal_rule(func, a, b, n)\n",
    "\n",
    "print(\"Trapezoidal rule integration result = \", I)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simpson's rules\n",
    "The trapezoidal rule uses a trapezoid to approximate the area under a segment. It uses a line connecting two points to approximate the integrant function. The idea of the Simpson's rules is to use higher-order polynomials to approximate the function. Two different Simpson's rules will be introduced. \n",
    "- on a segment with three points, a parabola can be fitted and this is called the 1/3 rule\n",
    "- on a segment with four points, a third-order polynomial can be fitted and it is called the 3/8 rule\n",
    "\n",
    "First, let's look at the Simpson's 1/3 rule. Let a segment be made of $[x_i,x_{i+2}]$ and $x_{i+1}$, with $x_{i+1}$ being  the mid-point. With the three points, a second-order Lagrange polynomial can be used to approximate the integrand\n",
    "  \\begin{eqnarray*}\n",
    "  f_2(x) &=& \\frac{(x-x_{i+1})(x-x_{i+2})}{(x_i-x_{i+1})(x_i-x_{i+2})} f(x_i) +  \n",
    "             \\frac{(x-x_{i})(x-x_{i+2})}{(x_{i+1}-x_{i})(x_{i+1}-x_{i+2})} f(x_{i+1}) \\\\\n",
    "         & & + \\frac{(x-x_{i})(x-x_{i+1})}{(x_{i+2}-x_{i})(x_{i+2}-x_{i+1})} f(x_{i+2})           \n",
    "  \\end{eqnarray*}\n",
    "The integration is then\n",
    "  \\begin{eqnarray*}\n",
    "    I &=& \\int_{x_{i}}^{x_{i+2}} f(x)dx \\approx \\int_{x_{i}}^{x_{i+2}} f_2(x)dx \\\\\n",
    "    &=& \\int_{x_{i}}^{x_{i+2}} [ \\frac{(x-x_{i+1})(x-x_{i+2})}{(x_i-x_{i+1})(x_i-x_{i+2})} f(x_i) +  \n",
    "             \\frac{(x-x_{i})(x-x_{i+2})}{(x_{i+1}-x_{i})(x_{i+1}-x_{i+2})} f(x_{i+1}) \\\\\n",
    "         & & + \\frac{(x-x_{i})(x-x_{i+1})}{(x_{i+2}-x_{i})(x_{i+2}-x_{i+1})} f(x_{i+2})    ] dx  \n",
    "  \\end{eqnarray*}\n",
    "\n",
    "Let $h=(x_{i+2}-x_i)/2$, the integration gives\n",
    "  \\begin{eqnarray*}\n",
    "  I &\\approx&   (x_{i+2}-x_i) \\frac{ \\left[f(x_i) + 4f(x_{i+1}) + f(x_{i+2}) \\right]}{6}  \\\\\n",
    "    & = &  \\frac{h}{{3}}\\left[f(x_i) + 4f(x_{i+1}) + f(x_{i+2}) \\right]\n",
    "  \\end{eqnarray*}\n",
    "The name \"Simpson's 1/3 rule\" is from the \"3\" in the denominator.\n",
    "\n",
    "If the integration interval $[a,b]$ is segmented with evenly distributed $n+1$ points where $n$ is an even number, there will be $n$ segments. Each segment has a length $h=(b-a)/n$. Then the total integration is \n",
    "  \\begin{eqnarray*}\n",
    "    I &\\approx& 2h \\frac{f(x_0)+4f(x_1)+f(x_2)}{6} + 2h \\frac{f(x_2)+4f(x_3)+f(x_4)}{6} \\\\\n",
    "      & &  + \\ldots + 2h \\frac{f(x_{n-2})+4f(x_{n-1})+f(x_n)}{6}\n",
    "  \\end{eqnarray*}\n",
    "It can be simplified as\n",
    "  \\begin{equation}\n",
    "  I \\approx \\frac{h}{3} \\left[f(x_0) + 4 \\sum_{i=1,3,5}^{n-1} f(x_i) + 2 \\sum_{i=2,4,6}^{n-2} f(x_i) + f(x_n) \\right]\n",
    "  \\end{equation}\n",
    "An equivalent way to write the formula is\n",
    "  \\begin{equation}\n",
    "  I \\approx \\frac{h}{3} \\sum_{i=1}^{n/2}  \\left[f(x_{2i-2}) + 4 f(x_{2i-1}) + f(x_{2i}) \\right]\n",
    "  \\end{equation}\n",
    "The error of the Simpson's 1/3 rule follows\n",
    "   \\begin{equation}\n",
    "      E \\propto -\\frac{(b-a)^5}{180n^4}\n",
    "   \\end{equation}\n",
    "which is fourth-order accurate because of the fact of $n^4$ in the denominator.  \n",
    "  \n",
    "Next, we show the use of the Simpson's 1/3 rule with a manual calculation example. Using the Simpson's 1/3 rule, calculate the following integration \n",
    "\\begin{equation}\n",
    "I = \\int_0^1 (x^2+1) dx\n",
    "\\end{equation}\n",
    "with \n",
    "- only one segment \n",
    "- two segments\n",
    "(Note the true solution is $I=\\frac{4}{3}$)  \n",
    "\n",
    "Similar to Simpson's 1/3 rule, a third-order Lagrange polynomial can be used to approximate the integrand. In this case,  Four points are needed for a segment. Let a segment be $[x_i,x_{i+3}]$, and $x_{i+1}$ and $x_{i+2}$ are two mid-points.\n",
    "Let $h=(x_{i+2}-x_i)/3$, similar analysis will give\n",
    "  \\begin{eqnarray*}\n",
    "  I &\\approx&  \\frac{{3}h}{{8}}\\left[f(x_i) + 3f(x_{i+1}) + 3f(x_{i+2}) + f(x_{i+3}) \\right]\n",
    "  \\end{eqnarray*}\n",
    "The name \"Simpson's 3/8 rule\" is from the \"3/8\" in the formula.\n",
    "\n",
    "If the integration interval $[a,b]$ is segmented with evenly distributed $n+1$ points, there will be $n$ segments. The number $n$ needs to be $3m$ to use Simpson's 3/8 rule, where $m$ is a positive integer. With these, the following statements can be made:\n",
    "- Each interval is $[x_{k-1},x_k]$ where $k$ =1, $3m$. \n",
    "- Each interval has a length of $h=\\frac{b-a}{3m}$.\n",
    "- The integration can be approximated as\n",
    "  \\begin{equation}\n",
    "  I =\\int_a^b f(x) dx \\approx \\frac{3h}{8} \\sum_{k=1}^m \\left[ f(x_{3k-3}) + 3 f(x_{3k-2}) + 3f(x_{3k-1}) + f(x_{3k}) \\right]\n",
    "  \\end{equation}\n",
    "  \n",
    "It can be shown that the error of the Simpson's 3/8 rule is \n",
    "  \\begin{eqnarray}\n",
    "  E \\propto - \\frac{(b-a)^5}{80n^4}\n",
    "  \\end{eqnarray}\n",
    "  \n",
    "Manual calculation example: Using the Simpson's 3/8 rule, calculate the following integration \n",
    "\\begin{equation}\n",
    "I = \\int_0^1 (x^2+1) dx\n",
    "\\end{equation}\n",
    "with \n",
    "- only one segment \n",
    "- two segments\n",
    "(Note the true solution is $I=\\frac{4}{3}$)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Double integration\n",
    "\n",
    "Double integral is integration over an area\n",
    "\\begin{equation}\n",
    "  \\int_c^d \\left(\\int_a^b f(x,y) dx \\right) dy\n",
    "\\end{equation}   \n",
    "which gives the (signed) volume between the $xy$ plane and the surface described by the function $f(x,y)$. The order of double integration is changeable because it will give same volume. \n",
    "\\begin{equation}\n",
    "  \\int_c^d \\left(\\int_a^b f(x,y) dx \\right) dy = \\int_a^b \\left(\\int_c^d f(x,y) dy \\right) dx\n",
    "\\end{equation}    \n",
    "\n",
    "<img src=\"integration_volume.png\" width=\"300\"/>\n",
    "<h3 align=\"center\">Figure. Double integration over an area.</h3> \n",
    "\n",
    "To perform numerical double integration, it is not that different from one-dimensional integration. The double integration can be performed as follows: \n",
    "- integration in one-dimension with methods like trapezoidal or Simpson's rule with the values of the second dimension held constant,\n",
    "- then integrate in the second dimension.\n",
    "    \n",
    "<img src=\"double_integral.png\" width=\"300\"/>\n",
    "<h3 align=\"center\">Figure. Numerical double integration.</h3>     \n",
    "\n",
    "The following table shows the tabulated value at each grid point $(x_i, y_j)$, where $i\\in [0,n]$ and $j \\in [0,m]$. The rows and columns can be used to performed the integration in one-dimension and the other. \n",
    "\n",
    "|  | $x_0$ | $x_1$ | $x_2$ |  $\\ldots$ | $x_n$\n",
    "| :---:  |  :----:   |  :----:   |  :----:   |  :----:   |  :----:   |\n",
    "|$y_0$ | $f(x_0,y_0)$ | $f(x_1,y_0)$  |  $f(x_2,y_0)$ | $\\ldots$  |  $f(x_n,y_0)$  |\n",
    "|$y_1$ | $f(x_0,y_1)$  | $f(x_1,y_1)$  |  $f(x_2,y_1)$ | $\\ldots$  |  $f(x_n,y_1)$  |\n",
    "|$y_2$ | $f(x_0,y_2)$  | $f(x_1,y_2)$  |  $f(x_2,y_2)$ | $\\ldots$  |  $f(x_n,y_2)$  |\n",
    "|$\\vdots$ |    |    | $\\vdots$   |    |\n",
    "|$y_m$ | $f(x_0,y_m)$  | $f(x_1,y_m)$  |  $f(x_2,y_m)$ | $\\ldots$  |  $f(x_n,y_m)$  |\n",
    "\n",
    "The next we will show an example: \n",
    "\\begin{equation}\n",
    "I =  \\int_{-1}^1 \\left(\\int_{-1}^1 \\sqrt{x^2 + y^2} dx \\right) dy\n",
    "\\end{equation}\n",
    "with $n$ = 10 and $m$ = 10.\n",
    "\n",
    "<img src=\"2D_3D_plot_function.png\" width=\"300\"/>\n",
    "<h3 align=\"center\">Figure. Example double integration.</h3>     \n",
    "\n",
    "The $(x,y)$ coordiantes where the function values are evaluated are given in the following table:\n",
    "\n",
    "|$x_0$ | $x_1$ | $x_2$ | $x_3$ | $x_4$ | $x_5$ | $x_6$ | $x_7$ | $x_8$ | $x_9$ | $x_{10}$ |\n",
    "| :---:  |  :----:   |  :----:   |  :----:   |  :----:   |  :----:   |  :----:   |  :----:   |  :----:   |  :----:   |  :----:   |\n",
    "|  -1. |  -0.8 | -0.6 | -0.4 | -0.2 |   0. |  0.2 |  0.4 |  0.6 |  0.8 | 1. |\n",
    "|$y_0$ | $y_1$ | $y_2$ | $y_3$ | $y_4$ | $y_5$ | $y_6$ | $y_7$ | $y_8$ | $y_9$ | $y_{10}$ |\n",
    "|  -1. |  -0.8 | -0.6 | -0.4 | -0.2 |   0. |  0.2 |  0.4 |  0.6 |  0.8 | 1. |  \n",
    "\n",
    "The process of doing the double integration can be demonstrated as in the following table. \n",
    "<img src=\"double_integration_example.png\" width=\"600\"/>\n",
    "<h3 align=\"center\">Figure. Example double integration.</h3> \n",
    "    \n",
    "Triple and multiple integrals are integration over a volume or higher dimensional space. The integration over one-dimension and followed by integration in other dimensions is no different than the doulbe integration.\n",
    "\n",
    "The following is an implementation of double integration. Trapezoidal rule is used in both directions. The result is also compared with the result from using Python's \"dblquad(...)\" function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "xs =  [-1.  -0.8 -0.6 -0.4 -0.2  0.   0.2  0.4  0.6  0.8  1. ]\n",
      "ys =  [-1.  -0.8 -0.6 -0.4 -0.2  0.   0.2  0.4  0.6  0.8  1. ]\n",
      "result from this program =  3.0824460805504503\n",
      "result from dblquad(...) =  3.0607828658561447\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy import integrate\n",
    "import matplotlib.pyplot as plt\n",
    "from math import exp  \n",
    "import sys\n",
    "\n",
    "#define the trapezoial rule: f is the integrand values\n",
    "#h is the segment size, and n is the number of segments \n",
    "def trapezoidal_rule(f, h, n):\n",
    "    s = (f[0] + f[n])\n",
    "    for i in range(1,n):\n",
    "        s = s + 2.0*f[i]\n",
    "    \n",
    "    return s*h/2.0\n",
    "\n",
    "#define the integrand fucntion\n",
    "def integrand_func(x, y): \n",
    "    return pow(pow(x, 2) + pow(y, 2), 0.5) \n",
    "\n",
    "#perform the double integration \n",
    "#The number of grid points in x and y are [0,n] and [0,m], respectively\n",
    "#The lower and upper bound in x and y are [xl,xu] and [yl,yu]\n",
    "def doubleIntegration(n, m, xl, xu, yl, yu): \n",
    "    # z stores the fucntion values at grid points \n",
    "    # g stores the integral the integration value corresponding to each y   \n",
    "    z = np.zeros((m+1,n+1))\n",
    "    g = np.zeros(m+1)\n",
    "    \n",
    "    xs = np.linspace(xl,xu,n+1)\n",
    "    ys = np.linspace(yl,yu,m+1)\n",
    "    \n",
    "    print(\"xs = \", xs)\n",
    "    print(\"ys = \", ys)\n",
    "\n",
    "    # calculate the segment size in x and y \n",
    "    delta_x = (xu - xl)/n \n",
    "    delta_y = (yu - yl)/m \n",
    "\n",
    "    # populate the integrand value table \n",
    "    for j in range(0, m+1): \n",
    "        for i in range(0, n+1): \n",
    "            z[j][i] = integrand_func(xs[i],ys[j]) \n",
    "            \n",
    "    #print(\"z = \\n\", z.round(2))        \n",
    "\n",
    "    # integrate with y held constant at each y value \n",
    "    for j in range(0, m+1): \n",
    "        g[j] = trapezoidal_rule(z[j,:],delta_x,n)    \n",
    "\n",
    "    #print(\"g = \\n\", g.reshape(m+1,1).round(2))\n",
    "    \n",
    "    # integrate in y direction with the values calculated above\n",
    "    answer = trapezoidal_rule(g,delta_y,m)\n",
    "    \n",
    "    \n",
    "    return answer \n",
    "\n",
    "\n",
    "# Application of the code above\n",
    "\n",
    "# xl and xu are upper and lower limit of x integral \n",
    "# yl and yu are upper and lower limit of y integral \n",
    "# n and m are number of segements in x and y\n",
    "xl, xu, n = -1, 1, 10\n",
    "yl, yu, m = -1, 1, 10\n",
    "    \n",
    "print(\"result from this program = \", doubleIntegration(n, m, xl, xu, yl, yu)) \n",
    "    \n",
    "print(\"result from dblquad(...) = \", (integrate.dblquad(integrand_func,xl,xu,yl,yu))[0])\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gauss quadrature\n",
    "Newton-Cotes methods introduced above use evenly distributed data points $[x_0,x_1,\\ldots,x_i,\\ldots]$. The Newton-Cotes integral is: \n",
    "- Trapezoidal rule on an interval $[x_i,x_{i+1}]$: \n",
    "   \\begin{equation}\n",
    "      I \\approx (x_{i+1}-x_i)\\frac{f(x_i)+f(x_{i+1})}{2}\n",
    "   \\end{equation}  \n",
    "- Simpson's 1/3 rule on an interval $[x_i,x_{i+2}]$:\n",
    "  \\begin{eqnarray*}\n",
    "  I &\\approx&   (x_{i+2}-x_i) \\frac{ \\left[f(x_i) + 4f(x_{i+1}) + f(x_{i+2}) \\right]}{6}  \n",
    "  \\end{eqnarray*}\n",
    "- Simpson's 3/8 rule on an interval $[x_i,x_{i+3}]$::\n",
    "    \\begin{eqnarray*}\n",
    "  I &\\approx& (x_{i+3}-x_i) \\frac{\\left[f(x_i) + 3f(x_{i+1}) + 3f(x_{i+2}) + f(x_{i+3}) \\right]}{8}\n",
    "  \\end{eqnarray*}\n",
    "  \n",
    "All these integrals are weighted-average of function values at evenly distributed points. As you recall, 1D integration is to calculate the area:\n",
    "  \\begin{eqnarray}\n",
    "  \\textrm{Area = width} \\times \\textrm{height}\n",
    "  \\end{eqnarray}\n",
    "\n",
    "<img src=\"Newton_Cotes_average_idea.png\" width=\"500\"/>\n",
    "<h3 align=\"center\">Figure. Average area for the Newton-Cotes methods.</h3>     \n",
    "\n",
    "These evenly distributed data points may not represent the integrand well. We have seen this in the interpolation chapter on the discussion of overfitting/oscillation with evenly distributed data points. Unevenly distributed data points may be better to capture the general trend of the integrand function. One such method is called Gauss quadrature.\n",
    "\n",
    "To come up with the Gauss quadrature methods, the task is to determine the data points and weighting coefficients to achieve maximum accuracy (or even exact solution). \n",
    "  \\begin{equation}\n",
    "  \\int_a^b f(x) dx \\approx \\sum_{i=0}^n a_i f(x_i)\n",
    "  \\end{equation}\n",
    "  using $n+1$ data points.\n",
    "An $(n+1)$-point Gauss quadrature rule is constructed to give an exact integration for polynomials of degree upto $2n+1$ (=$2(n+1)-1$). In this task, there are $2n+2$ unknowns:\n",
    "- $n+1$ data points $x_i$ in $[a,b]$\n",
    "- $n+1$ weighting coefficients $a_i$\n",
    "\n",
    "The $2n+2$ unknowns can be derived based on the requirement that the rule exactly integrates a polynomial of degree upto $2n+1$. Because of this reason, the Gauss quadrature rule is more accurate if the integrand $f(x)$ behaves close to a polynomial of degree upto $2n+1$.\n",
    "\n",
    "For simplicity, we will consider the integration region to be $[a,b]=[-1,1]$. Generalization will be made later because some linear transformation can be used later for different integration region. We will derive the data points and weighting coefficients with an example $n$ = 1 (two data points):\n",
    "\\begin{equation}\n",
    "  \\int_{-1}^{1} f(x')dx' \\approx a_0 f(x'_0) + a_1 f(x'_1)\n",
    "\\end{equation}\n",
    "The polynomial function $f(x')$ has a degree of $2n+1$ = 3, i.e., a cubic polynomial:\n",
    "\\begin{equation}\n",
    "  f(x') = b_0 + b_1 x' + b_2 x'^2 + b_3 x'^3\n",
    "\\end{equation}\n",
    "So plug into the definition, we have\n",
    "  \\begin{eqnarray*}\n",
    "     \\int_{-1}^{1} f(x') dx' &=& \\int_{-1}^{1} \\left(b_0 + b_1 x' + b_2 x'^2 + b_3 x'^3 \\right) dx' \\\\\n",
    "     &=& a_0 \\left(b_0 + b_1 {x'}_0 + b_2 {x'}_0^2 + b_3 {x'}_0^3 \\right) \\\\\n",
    "     & & + a_1 \\left(b_0 + b_1 {x'}_1 + b_2 {x'}_1^2 + b_3 {x'}_1^3 \\right)\n",
    "  \\end{eqnarray*}\n",
    "Collect similar terms with $b_0$, $b_1$, $b_2$ and $b_3$, we have\n",
    "  \\begin{eqnarray*}\n",
    "  b_0 \\left(a_0 + a_1 -\\int_{-1}^1 dx'\\right) + b_1 \\left(a_0 {x'}_0 + a_1 {x'}_1 -\\int_{-1}^1 {x'} d{x'}\\right) & & \\\\\n",
    "  b_2 \\left(a_0 {x'}_0^2 + a_1 {x'}_1^2 -\\int_{-1}^1 {x'}^2 d{x'}\\right) + b_3 \\left(a_0 {x'}_0^3 + a_1 {x'}_1^3 -\\int_{-1}^1 {x'}^3 d{x'}\\right)  = 0\n",
    "  \\end{eqnarray*}\n",
    "Since the polynomial $f(x')$ is arbitrary, i.e., the coefficients $b_0$, $b_1$, $b_2$, and $b_3$ are arbitrary too. Then, for the above equation to be zero, all the terms in the parentheses have to be zero. \n",
    "   \\begin{eqnarray*}\n",
    "   a_0 + a_1 -\\int_{-1}^1 dx' = 0 \\quad\n",
    "   a_0 {x'}_0 + a_1 {x'}_1 -\\int_{-1}^1 x' dx' &=&0 \\\\\n",
    "   a_0 {x'}_0^2 + a_1 {x'}_1^2 -\\int_{-1}^1 x'^2 dx' = 0 \\quad\n",
    "   a_0 {x'}_0^3 + a_1 {x'}_1^3 -\\int_{-1}^1 x'^3 dx' &=& 0\n",
    "   \\end{eqnarray*}  \n",
    "\n",
    "We have four unknowns ($a_0$, $a_1$, $x_0$ and $x_1$) with four equations. The solution is as follows:\n",
    "  \\begin{eqnarray*}\n",
    "  a_0 = 1 \\quad a_1 = 1 \\quad {x'}_0 = -\\frac{\\sqrt{3}}{3} \\quad {x'}_1 = \\frac{\\sqrt{3}}{3}\n",
    "  \\end{eqnarray*}\n",
    "Therefore\n",
    "    \\begin{equation}\n",
    "  \\int_{-1}^{1} f(x')dx' \\approx a_0 f(x'_0) + a_1 f(x'_1) = f(-\\frac{\\sqrt{3}}{3}) + f(\\frac{\\sqrt{3}}{3})\n",
    "  \\end{equation}\n",
    "This is our 2-point Gauss quadrature rule! It can exactly integrate a polynomial of degree upto 3!\n",
    "\n",
    "Manual calculation example: integrate $f(x')=x'^2$ on $[-1,1]$ with $n=1$\n",
    "\n",
    "Manual calculation example: integrate $f(x')=x'^3+2x'^2+1$ on $[-1,1]$ with $n=1$\n",
    "\n",
    "Manual calculation example: integrate $f(x')=\\cos(x')$ on $[-1,1]$ with $n=1$\n",
    "\n",
    "We derived the Gauss quadrature on a nominal integration interval $x' \\in [-1,1]$. What about a general interval $[a,b]$?\n",
    "The solution is a simple linear transformation or change of variables from $x'$ to $x$. Let $x=\\frac{b-a}{2}x'+\\frac{b+a}{2}$, it is easy to show that \n",
    "  \\begin{equation}\n",
    "  \\textrm{when } x' = -1, \\quad x = a,\\quad \\textrm{and when } x' = 1, \\quad x =b \n",
    "  \\end{equation}\n",
    "With this linear transformation, we have\n",
    "  \\begin{equation}\n",
    "  dx = \\frac{b-a}{2} dx'\n",
    "  \\end{equation}\n",
    "  then\n",
    "  \\begin{equation}\n",
    "  \\int_a^b f(x)dx = \\int_{-1}^1 f\\left(\\frac{(b-a)x'+b+a}{2} \\right) \\frac{b-a}{2} dx'\n",
    "  \\end{equation}\n",
    "  \n",
    "Manual calculation example: integrate $f(x)=x^2$ on $[a,b]=[1,3]$ with $n=1$ (note: analytical solution = $\\frac{26}{3}\\approx 8.6667$)  \n",
    "\n",
    "The previous derivation and examples use $n=1$ and we get:\n",
    "    \\begin{equation}\n",
    "  \\int_{-1}^{1} f(x')dx' \\approx a_0 f(x'_0) + a_1 f(x'_1)\n",
    "  \\end{equation}\n",
    "  with \n",
    "    \\begin{eqnarray*}\n",
    "  a_0 = 1 \\quad a_1 = 1 \\quad {x'}_0 = -\\frac{\\sqrt{3}}{3} \\quad {x'}_1 = \\frac{\\sqrt{3}}{3}\n",
    "  \\end{eqnarray*}\n",
    "This is exact for polynomials up to the degree of $2n+1=3$. Can we go for \\bverb|higher-order| polynomials?\n",
    "    \\begin{equation}\n",
    "  \\int_a^b f(x') dx' \\approx \\sum_{i=0}^n a_i f(x'_i)\n",
    "  \\end{equation}\n",
    "It turns out $x'_i$ are roots of the Legendre polynomials:\n",
    "  \\begin{eqnarray*}\n",
    "  P_0 &=& 1 \\\\\n",
    "  P_1 &=& x \\\\\n",
    "  P_2 &=& \\frac{3x^2-1}{2} \\\\\n",
    "  P_3 &=& \\frac{5x^3-3x}{2} \\\\\n",
    "  \\ldots &=& \\ldots\n",
    "  \\end{eqnarray*}\n",
    "and a recursive relationship can be used to get higher-order Legendre polynomials:\n",
    "  \\begin{equation}\n",
    "  P_n(x) = \\frac{2n-1}{n}x P_{n-1}(x) - \\frac{n-1}{n} P_{n-2}(x), \\quad n=2,3,4,\\ldots\n",
    "  \\end{equation}\n",
    "For a $n$-point Gauss quadrature, the data points should be at the roots of $P_{n}$. For example, our 2-point Gauss quadrature rule, the data points are at the roots of $P_2$. The following figure shows the Legendre polynomials. Note the roots of each polynomial labeled as scatter points.\n",
    "\n",
    "<img src=\"legendre_polynomials.png\" width=\"500\"/>\n",
    "<h3 align=\"center\">Figure. Legendre polynomials and their roots.</h3> \n",
    "The weights $a_i$ can be derived as above: to make the numerical integration accurate for polynomials upto to the order of $2n+1$. The following table shows the roots of the Legendre polynomials and the integration weights.\n",
    "\n",
    "| Number of points | Points $x_i$ |  Weights $a_i$ |\n",
    "| :---:        |    :----:   |          :---: |\n",
    "| 1 | 0  |  2 |\n",
    "| 2 | $\\pm$ 0.57735... | 1 |\n",
    "| 3 | 0, $\\pm$ 0.774597... | 0.888889... and 0.555556... | \n",
    "| 4 | $\\pm$0.339981..., | $\\pm$ 0.861136... | 0.652145... and 0.347855... |\n",
    "| $\\vdots$ | $\\vdots$ |  $\\vdots$ |\n",
    "\n",
    "Keep in mind, these points $x_i$ are for the nominal interval $[-1,1]$. For general interval $[a,b]$, a linear transformation is needed. Weights $a_i$ do not depend on the interval."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python libraries for numerical integration\n",
    "\n",
    "SciPy provides an \"integrate\" module:\n",
    "[scipy.integrate documentation](https://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html). It provides two categories of numerical integrations:\n",
    "\n",
    "* integrating functions given function values at points: It provides trapezoidal (**trapz**), Simpson's rule (**simps**), etc. The documentation of the function **simps** is not clear on what Simpson's rule it implements. From reading the source code, it seems that the Simpson's 1/3 rule is implemented if the data points are evenly distributed. For unevenly distribbuted data points, the formula is slightly different. Some weighted-average 1/3 rule between two adjacent segments is used. \n",
    "\n",
    "* integrating functions given function definitions: It provides single, double, triple, and multiple integrations (**quad**, **dblquad**, **tplquad**, etc.).\n",
    "\n",
    "The following is an example of the first category of numerical integration functionality in Python. The integration data points are generated on the interval $[1,5]$ with $x$ spacing $dx$ =1. The values $y$ are from the function $x^2$. The analytical solution is $\\frac{124}{3}\\approx41.33333...$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "trapz result =  42.0\n",
      "trapz result with refined data points =  41.5\n",
      "trapz result with very refined data points =  41.333333440021335\n",
      "simps result =  41.333333333344\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy import integrate\n",
    "\n",
    "#set up the x and y values with only five data points\n",
    "x = [1,2,3,4,5]\n",
    "y = [1,4,9,16,25]\n",
    "\n",
    "print(\"trapz result = \", integrate.trapz(y,x))\n",
    "#print(\"trapz result = \", np.trapz(y,x))   !this also works\n",
    "\n",
    "#set up the x and y values with ten data points\n",
    "x = [1,1.5,2,2.5,3,3.5,4,4.5,5]\n",
    "y = [1,2.25,4,6.25,9,12.25,16,20.25,25]\n",
    "\n",
    "print(\"trapz result with refined data points = \", integrate.trapz(y,x))\n",
    "\n",
    "#set up the x and y values with a lot more data points\n",
    "x = np.linspace(1,5,10000)\n",
    "y = x**2\n",
    "print(\"trapz result with very refined data points = \", integrate.trapz(y,x))\n",
    "\n",
    "#similary, you can also use the Simpson's rule\n",
    "print(\"simps result = \", integrate.simps(y,x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following is an example of the second category of numerical integration functionality in Python. It uses the **quad()** function of the **scipy.integrate** module. The quad function returns the integral results and an esitmate of the errror. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The integration result from quad(...) function is 41.333333 (+-4.58892e-13)\n"
     ]
    }
   ],
   "source": [
    "from scipy.integrate import quad\n",
    "\n",
    "# function we want to integrate\n",
    "def func(x):\n",
    "    return x**2\n",
    "\n",
    "a = 1.0\n",
    "b = 5.0\n",
    "\n",
    "# call quad to integrate \n",
    "res, err = quad(func, a, b)\n",
    "\n",
    "print(\"The integration result from quad(...) function is %f (+-%g)\"%(res, err))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gauss quadrature and Legendre polynomials\n",
    "\n",
    "Python's SciPy integrate module implements the Gauss quadrature in its **quadrature** function. Indeed, it implements an adaptive version of the Gauss quadrature with a specified tolerance. If the tolerance is not specified, a default value (~$10^{-8}$) is used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The integration result from quadrature(...) function is 41.333333 (+-7.10543e-15)\n"
     ]
    }
   ],
   "source": [
    "from scipy import integrate\n",
    "\n",
    "# function we want to integrate\n",
    "def func(x):\n",
    "    return x**2\n",
    "\n",
    "a = 1.0\n",
    "b = 5.0\n",
    "\n",
    "# call quadrature to integrate \n",
    "res, err = integrate.quadrature(func, a, b)\n",
    "\n",
    "print(\"The integration result from quadrature(...) function is %f (+-%g)\"%(res, err))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We know the data points $x_i$ for the Gauss quadrature are the roots of the Legendre polynomials. In Python, you can examine the Legendre polynomials, their roots, and even the Gauss quadrature weights. The following is an example Python code to plot several Legendre polynomials and their corresponding roots. In **SciPy**, the **legendre** sub-module is in its **special** module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaoAAAEQCAYAAADh3jDlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd3yN1x/A8c/J3iGJvalZtfeuaosORVGrRo3alNo1atSmNWpTs4pWbWpvFYSixEokIUQiS9bNvef3x036iwgZ7grn/XrlleYZ53yj995vnuc553uElBJFURRFsVRW5g5AURRFUV5GJSpFURTFoqlEpSiKolg0lagURVEUi6YSlaIoimLRbMwdgKXz8vKSRYsWNXcYiqIo2cr58+cfSylzGaItlajSUbRoUby9vc0dhqIoSrYihPA3VFvq1p+iKIpi0VSiUhRFUSyaSlSKoiiKRVOJSlEURbFoKlEpiqIoFi1bJCohRHUhxG0hRNd0jusohDgvhPAWQswWQogU+zyEENuFECeEEGeEEFWMHriiKIryyiw+UQkhWgJDgPB0jisPzAY+BGoAVYC+KQ5ZBFyUUtYDRgF/CiHsjRK0oiiKYjAWn6iAc1LKDkBUOsd9BeyWUj6WUuqAlcDXoL+aAtoAKwCklIeBBODj9DqXqGVQFEVRMkpqtdy/ftSgbVr8hF8pZWAGD60ObE/x8zXgbSGEI/qrq3gp5b1U+6sBW1/W6M3HwbRbcjoTESuKoryZ7HWxNPGfRdVtdw3abna4osqoPDx7ezACEIBXGvuS9+dOqyEhRK+k51zeiYSj5akx4lUURXk9SEm1uFPMCulJ7IPb+k9eA7L4KyoDeNk/WZr7pJRLgaUATsWcZO2qlxhabagxYlMURcnewu/B7uHgu4fd+UtT4FAU8WWLwvXrBuvidbqiegTkSPGzOyCBkKR97qmOd0/a/lLu9u6s/3c9QdFBhopTURQl+9Nq4MQ8WFgT7h4jvskEVtq6UiIYCr3/qUG7ep0S1TmgdIqfywFXpZSxwAXAQQhRKNX+c+k1mtspN9bCmh8v/GjQYBVFUbKte2dgSQM4MB5KNIZ+Z9mYIyd5rwQD4NqwoUG7y7aJSgjhJYQ4njSiD2A50FwI4SmEsAK6AosBpJShwGage9K5DQE7YFd6/dha2dK5XGf23N3DlcdXjPCbKIqiZBMxYbB9AKz8EOIi4YuN8MV6wh1cWfrPUt5/4Il1Li/sy5Y1aLcWn6iEEFWFEEeASsBIIcTvSbscgTKAE4CU8gowDNgPnAV80M+dStYXqCqEOAFMA1pIKeMyEsNX73yFh4MHs7xnIaUarq4oyhtGSvDZCAuqwcX1UGcA9DsLZZoDsOTyEmLjoyl1MwaX+g1IUWvBICx+MIWU8jzQKI3tAUCuVNvWA+tf0E4YkKUbp862zvSt2JfJZydzJOAI7xZ+NyvNKIqiZD8hvrDrG/A7DgVrwMdzIW/5/3YHRAbw641f+Yr6EHUIlwb1DR6CxV9RWYpWpVpRzL0Yc87PQaPTmDscRVEU49LEwqHJ8HMdCL4MH8+D7vueSVIA8y7Mw9bKlo8fFQBra5zr1DF4KCpRZZCtlS1DqgzBL9KP331/T/8ERVGU7OrWAVhUC47NhPKtof95qNYNrJ5NGT6PfNjvv5+ub3dFd8Ybx8qVsHZzM3g4KlFlQqNCjaiWpxqLLi0iOiHa3OEoiqIYVlQwbO4G61qDlQ18uR1aLQGXXM8dKqVktvdsvBy96JT7I+Kv/YtLA8OO9kumElUmCCEYVm0YYXFhrLyy0tzhKIqiGIZOC38vgwXV4fouaDQK+pyC4i9OPAfvHcQnxId+lfqhO30ewCjPp0Alqkx72+ttmhdrzpprawh+GmzucBRFUV7NfR9Y3gR2D4MCVaDvaWg0EmxevLiERqth7vm5vJXjLT576zOijx3DJndu7EuXfuE5r0IlqiwYWGUgOqljwcUF5g5FURQla+IiYc9IWPYuRARC6xXQeRt4lkj31N98f+Ne1D2GVB2CtVby9NQpnBvUN/iw9GQqUWVBAZcCdCrbie23t3Mj7Ia5w1EURck4KeHan7CwBpxdDNW6Q/9z8M7nkIFEE5kQyeJLi6mZryb1C9Qn1scHXVQULg0aGC1klaiyqEeFHrjZuzHTe6aaBKwoSvbwxA82tIXfvgRnL+hxAD6aDY450j012fJ/lhMRH8HQqkMRQhB97DjY2BhlWHoylaiyyM3OjT4V+3D2wVkOBRwydziKoigvlpgAx+fAwlrgfwo+/AF6HoGC1TLVjH+kP2uvreXTEp9S1lNfJin62DGcqlTB2sXFCIHrqUT1CtqVbsdbOd5i5rmZxGvjzR2OoijK8/xP6QvIHpwIJZtAv7+hdl+wznxhopnnZmJvbc/gqoMB0Dx8SPyNG0Yb7ZdMJapXYGNlw4gaIwiKDmLN1TXmDkdRFOX/nobCn/1gVTNIeArtN0G7deBeIEvNnQg6wdHAo/Su0BsvRy9AfzUF4GzE51OgEtUrq5WvFu8Vfo9l/yzj4dOH5g5HUZQ3nZT6wrELqsGlX6HuYOh3Bko3zXKTGq2G6X9Pp4hbETqV7fTf9qfHjmOTNy/2JUsaIvIXUonKAIZVG4ZWp2XuhbnmDkVRlDfZo+uw+iP4sy94lYLex+D9iWDn/ErNbri+Ab9IP4ZXH46ttS0AMiGBp6dO4VLfeMPSk6lEZQAFXQvStXxXdt3Zhc8jH3OHoyjKmyYhBg5MhMV14dE1+HQ+dNsDed5+5aYfxz5m8aXF1C9QnwYF/3+LL+aiD7qnT3FpaNzbfpANlvkQQjigXwCxDPp4R0sp96dx3F7AIcUmW+AdKaXbC/ZHSyk/NlScX5X/ij9v/ckPf//AhuYbsLayNlTTiqIoL+a7X19VItwfKnaADybph54byPyL84nTxjG8+vBntkcfOwq2tjjVqm2wvl7E4hMVMAEQUspaQohSwBkhRFkpZeoHQneklH2TfxBCtOTZ9aeCpZRdjRWkk60T31T9hhHHR7Dt1jZal2ptrK4URVEg8j7sHamfvOtVCrrshGKGHX139fFV/rj5B13e7kJR96LP7Ht67DhOVati7fJqtxUzwqJv/SUtKd8DWAEgpfQFLgKdUh+bMkkl6QisM3aMKTUr1owquavw08WfiEyINGXXiqK8KbSJcOZnfQFZ333QeCx8fdLgSUpKyQ9//4CHgwe9K/R+Zp/mwQPib97Epb5xh6Uns+hEBRQHPIHrKbZdA146S00I4Y5+6frDKTY7CyHWCyGOCyG2CSEqGjpYIQQja4zkSdwTFl9abOjmFUV50wWdh+WN9VdShWrqC8g2+BZs7Aze1c47O7kUcolBVQbhYvfsZN7oY8cBTPJ8Ciw/UeVJ+h6eYlsEkDud81oDW6WUuhTbbgPjpZT10V+hnRBC5DNYpEnKepalVclWbPx3I3fC7xi6eUVR3kRxEbBrGCx7D6IewueroNNW8ChulO5iNDHMOz+P8p7lafFWi+f2Rx8/hk3+fNiVSL+ArSFYeqJ6kfTGQj53209KOVJKeSvpv3egvzL7Ms3GheglhPAWQniHhIRkOriBVQbiaOPIjHMzVB1ARVGyTkq4slV/m+/ccqjRC/r/DeVbZaiAbFYt+2cZj2IfMbLmSKzEs2lCJiQQc+o0Lg0aGH1YejJLT1SPkr6nrJjonmL7c4QQBYAcUsp/0mn7HlAkrR1SyqVSympSymq5cj2/smV6PBw86FOpDyfvn+Ro4NFMn68oikLYHf1Ku1u6g2te6HkIms8AB3ejdhsQGcAvV3/hk+KfUDHX809IYi5cQBcTY9Rq6alZeqK6DYQBKVfjKgece8k57YENKTcIIXILIbqmOi4PcN8AMabpizJfUNy9ODPOzSBBm2CsbhRFed0kxsOxmbCoNgT8DU2nQ8/D+kUNTWCm90xsrWz/q+eXWvSRowhbW5xr1jRJPGDhiSrpGdMyoDuAEKIk+kES64UQZYUQB4UQqScsfUGqRAU4AUOFEM5J7VQHqgObjBW7rZUtI2qMICAqQC1bryhKxvidgMX14NBkKNVUf5uv1tdgonmZxwKPcTjgML0q9CK30/NDAaSURB04gFOtWlg5G39YejKLTlRJJgBCCHEG2Ai0l1IGo78FWAb9xF7QH1QOCJVSPkjVRjCwDfhLCHEMmA+0lFLeNGbgdfLXoWnRpiy7vIx7kfeM2ZWiKNnZ08ewra++/FFiPHTYDG1/Abf8JgshNjGWqWenUsK9BF+WS/PxPXHXrqEJDMSt6YcmiwuywYRfKWUc0DWN7WeAAqm2XQOe+xdMauO7pC+T+rb6txwPOs6Us1NY3GSxyR4+KoqSDeh04LMO/hoH8dFQ7xv9cHM7J5OHsvTyUoKig1j14ar/6vmlFrV3H1hb49K4sUljyw5XVNlabqfcDKg8gFP3T7HPf5+5w1EUxVI8vAarm8P2AZCrLHx9ApqMN0uSuh1+m9VXVtOiRAuq5U17mqqUksj9+3CuWRObnDlNGp9KVCbwRekvKOdZjhl/zyA6Idrc4SiKYk4JT+Gv8bCkPoTcgBYLodtuyF3GLOFIKZl0ZhLOds58U+2bFx4Xf+MGGv97uJr4th+oRGUS1lbWjKs1jsexj1ngs8Dc4SiKYi6++/TLwZ+cBxW/gP7eULmTUedEpWf77e2cf3ieIVWG4OHg8cLjIvftAysrXJs0MWF0eipRmcjbXm/TrnQ7Nl7fyNXQq+YOR1EUU4oIgk2dYENb/a29rrv1V1LOnmYNKzwunNnes6mUqxItS7Z84XFSSqL27sOpRg1sPF6czIxFJSoTGlhlIB4OHkw6PQmtTmvucBRFMTZtIpxeCAtrwM2/4L1x0Ps4FK1r7sgAmHdhHpEJkXxX+7vnKlCklHDrFgl37+L24QcmjO7/VKIyIVc7V76t9i1XQ6/ym+9v5g5HURRjCvSGpY1g32goXBv6nYX6Q41SQDYrfB75sPXmVjqX60ypnKVeemzkvv0ghFlu+4FKVCbXrFgzauWrxU8XfuJx7GNzh6MoiqHFhsPOb2B5E4h5DG3XQMfNkLOouSP7j0an4fsz35PXOS99KvZJ9/iofftwqloVmyyUlDMElahMTAjB2FpjSdAmMOPcDHOHoyiKoUgJlzfrC8ieXwU1v4Z+f0O5FmYdLJGWDf9u4OaTm4yqMQon25cPh4+/c4f4mzdx/dD0o/2SqURlBkXcivDVO1+x5+4eTt0/Ze5wFEV5VaG3Ye1n8HsPcC+gr83XbBo4uJk7suc8iH7AQp+FNCrYiMaF05+4G7VPP//T9QPzPJ8ClajM5qt3vqKwa2Gmnp1KvDbe3OEoipIVifFwZLq+gGzQBWg+C3ochPyVzB3ZC037exoAo2qOytDxkfv241ilCrZ50lsG0HhUojITe2t7xtQag3+kP8v/WW7ucBRFyaw7R+HnOnBkKpT9GPqfgxo9TVZANisO3zvMoYBDfF3xa/K7pF9HMMHPj/jr18022i+ZSlRmVCd/HT4q/hHLLy/nRtgNc4ejKEpGRIfA771gzaeg0+pX2v18pX7NKAsWmRDJpDOTKJWzFJ3Ldc7YOfv/AsD1/feNGVq6VKIysxHVR+Bm78a4U+NI1CWaOxxFUV5EpwPvVbCgKlz5XV88tu9peMs8Q7Yza7b3bMLiwvi+7vfYWqVddDa1qH37cKhYAdv8pqvinhaVqMwsp0NORtUcxbXQa6y5tsbc4SiKkpbgK7DyQ9g5GPJWgD6noPFYsHU0d2QZcvr+aX6/+Ttd3u7C255vZ+ichMBA4q5exe0D8432S6YSlQX4sMiHvFf4PRb5LMIvws/c4SiKkizhKewfC0saQNht+GwxdNkBuV4+QdaSxGhimHh6IkXdimZozlSy/0b7mfn5FGSDRCWEcBBCrBZCnBFCeAsh0vxXE0J0FUL4CCGOpPiql2K/hxBiuxDiRFJbplnXOQOEEIypOQY7azvGnxqPTurMHZKiKNd3w8KacGq+vnBsf2+o1N7i5kSl58cLP3I/+j7f1/0eBxuHDJ8XuW8/DuXLY1ewoBGjyxiLT1QkrfArpawFdAB+FULkecGxg6WUjVJ8nUixbxFwUUpZDxgF/CmEsDdq5JmQyykXw6sP58KjC2y6scnc4SjKmys8ADZ2gF/bg70rdN8Hn/4ETqYvxvqqLj66yMbrG2lfpj2Vc1fO8HmaoCDiLl+2iKspsPBEJYSwAnoAKwCklL7ARaBTJtvxANqkaOcwkAB8bMh4X1WLEi2om78uc8/PJSg6yNzhKMqbRavRXz0trAm3D0GTidD7GBSuZe7IsiReG8+4k+PI75KfQVUGZercyL/0o/3czDjJNyWLTlRAccATuJ5i2zUg7SUoobcQ4mjSbb++KbZXAeKllPcy2I5ZCCEYV3scAsHEUxORUpo7JEV5MwT8rS8gu38sFKuvLyBbbzC8YEn27OBnn5/xi/RjXO1x6ZZJSi1q337sy5bFrkgRI0WXOZaeqJJv8YWn2BYBpDVF+iGwH2gEtAX6CyEGpGgnPNXxL2oHIUSvpOdh3iEhIVkMPWvyu+RnSNUhnH5wmm23tpm0b0V548Q+gR2DYMX7+v9utw7a/wo5LeMDOquuhl5l9dXVtCrZijr562TqXM3Dh8RevGj2Sb4pWXqiepHnnmZKKfdIKVdJvUfAfKDv86e+vJ2ktpZKKatJKavlMkO14Lal21I1T1Vmes8kJMa0iVJR3ghSwqVNML8aXFgLtfvrr6LKfpLtBkukptFqGHdyHJ4OngytNjTT50clT/K1gGHpySw9UT1K+p4jxTb3FNtf5h6Q/GfRo6TzUspoOyZnJayYWGciCdoEJp2ZpG4BKoohPb4Jv3wCf/TSXzn1OgIfTtEPnHgNrLiyAt8nvoytNRY3u8wXxY3ctxf7kiWxL17MCNFljaUnqttAGFA6xbZywLnUBwohhqfalAe4n/TfFwAHIUSh9NqxFEXcitCvUj8OBxxmn98+c4ejKNmfJg4OT9XX53twGT6aA1/9BfkqmDsyg7n55CZLLi+hWdFmvFv43Uyfr7l/n9jzF3BtajlXU2DhiUpKqQOWAd0BhBAlgUrAeiFEWSHEQSFEcgXIZkKIhknHOQG9gbVJ7YQCm1O00xCwA3aZ8NfJtM7lOlPeszxTz05Viywqyqu4fQh+rg1Hp+vXh+p/Dqp/ZdEFZDNLo9Pw3cnvcLV1ZWTNkVlqI3zbNpAS9xafGTi6V2Nj7gAyYAKwWAhxBn287aWUwUKIokAZwBbQArOBcUlD2l2Ag8APKdrpC6wWQpwArIEWUso4U/0SWWFjZcPkepNpu6MtE05NYH7j+Yhsfv8cQBsdjeb+/f++Eh88QBN0H01wMMLODhtPT2y8PLH29MLG0xNrTw9sPL2w8fLExssLYZt9R2IpJhb1UL8U/JUt4FECOm+DEpm/0sgOll5eytXQq8xuOBsPh8zP+ZI6HRG//4FTrVrYFSxghAizzuITVVIy6ZrG9jNAgRQ/7wR2vqSdMOBTI4RoVCVylGBI1SFMPzedrTe38nmpz80dUqbF+foSuXMX0cePowkKQhcZ+ewBtrbY5suHbd68yPh4Yv/5B+3jx+hiYp5rSzg64lSjOi516+Jcrx52xYq9FslbMTCdFrxXwsFJkBgLDUdAvW/ANuOVGbKTSyGXWHZ5GZ+W+JQPimZttF7MOW80gYHkGjTQwNG9OotPVAp0KNuBI4FHmHFuBjXy1qCwW2Fzh5QuTVAQEbt2E7lzJ/G+vmBtjVP16jhVroxt/nzY5s+Pbf782OTLj00uL4TV83ehdbGxJIaGoQ19TGJoKImhocRfv8HTkyd5ePQYADb58uFct44+cdWujXWOHM+1o7xhHlyCnUMg6DwUa6h/FuX1lrmjMpoYTQyjj48mt1NuRtbI2i0/gIjft2Ll6mr2JT3SohJVNmAlrJhcdzKttrdi9InRrG66Ghsry/tfl/jkCZF79hC5cxexFy4A4Fi5Mnm+G4tb06bYeHpmqj0rR0f9LYg0bkMkBAbx9NRJnp44SdT+v4jYshWEwKHCO+T4/HPcP/kEK4fX869n5QXio+DwD3D2Z3DyhFbL4J022X64eXpme88mICqAFR+uwNUuayMXtdHRRO7bj3uLFhb5vhFq6PPLVatWTXp7e5s7DAB239nNiOMjGFB5AL0q9DJ3OP/RPHzIo5mziNy7FxITsS/5Fm4ff4LbR81NUtBSJiYSd+UK0Sf1SSv+xg2sc+QgR7t25OzQHts8LyoNqbwWpITrO2HPCIi8D9W6wXvjwDGnuSMzumOBx+h3sB9d3+6apTlTyZ789hvB48ZT9LdNOFYwzChIIcR5KaVBqv+oRJUOS0pUAMOPDucv/79Y99G6DK8rYyxSoyFs3Xoez5+P1GrJ2aED7p99hkNp8y2BIKUk5tw5wtasIfrgIbC2xq1pUzy6fInjO++YLS7FSMLvwe5vwXcv5HkHPp4LhaqbOyqTeBL3hJZ/tsTD0YNfP/oVO2u7LLfl90V7tNFRFN+xw2DPfA2ZqCzv/pHyUmNqjeH8o/OMOj6K3z7+LVNl+w0p5vx5gid+T7yvL84NG5B37FjsChVK/0QjE0LgXKMGzjVqkBAQwJN16wnfsoXInTtxrFwZjy874/r++wgb9dLP1rQaOL1QP9wcAR9Mhpp9wPrN+P8qpWTi6YlEJkSy5P0lr5Sk4m/fJtbHh9zDh1vswCSLnkelPM/d3p1JdSdxN+Iu8y7MM3n/iaGh3B85Cv+OndBGR1FwwXwKLV5sEUkqNbtChcgzaiRvHT1CntGjSQwNJWjIN9z5tAXRx4+bOzwlq+6d0S9keGA8lGisL31UZ8Abk6QAtt/ezsF7BxlQeQClPUqnf8JLhP/+O1hb4/7pJwaKzvBUosqG6uSvQ8eyHVn/73pO3T9lkj6lVsuTjRu53aw5Ebt24dmzJyV27sS1SROL/SssmbWLCx5fdqbEnt0U+OlHpDaRgJ69COj9NfF37po7PCWjYsJg+wD9kvBxkfDFRvhiPeSwvD+SjCkoOogf/v6Bqnmq8mW5L1+pLanREPHndlwaNcLGy8tAERqeSlTZ1OAqgynmXozvTnxHRHyEUfvSBAfj174DwRO/x6FsWYpv+4PcQ7/ByilzSweYm7C2xu2DDyixYwe5v/2WmPPnufPppzz8YRra1HO7FMshJfhsgAXV4OJ6qDNQfxVVprm5IzM5rU7L6OOjAZhSbwrWr1hZI/r4CbSPH5OjVUtDhGc0KlFlUw42DvxQ/wfC4sKYcmaK0fqJ8/XF74v2JNy+Tf6ZMym8ehX2JUoYrT9TEHZ2eH7VnRJ795Cj5WeErVnD7Q+b8uTXTUit1tzhKSmF3NAXkN3WR19Zovcx+GAS2LuYOzKzWHNtDRceXWBUjVEUcHn16hERf/yOtacnLg0aGCA641GJKht72/Nt+lTqwx6/PWy/vd3g7T89cxb/jp1Ap6PI+nW4f/Kxxd/mywwbLy/yTZpEsa1bsCtRnOAJE7jbqjVPz/5t7tAUTay+qsTPdSH4Mnw8T78kfN7y5o7MbK6FXmP+xfk0KdyET0u8epGdxNBQog4fwb1FC4svS6YSVTbXvXx3quapyuQzk7kTccdg7Ubs3EVAz57Y5MlN0V834lCmjMHatjQO5cpRZO1aCsybizYqkntduhD8/SR0sbHmDu3NdOsALKoFx2dB+dbQ/7x+blQa1UveFNEJ0Qw7OgwPBw/G1x5vkD8YI3bsgMREi7/tBypRZXs2VjZMrz8dB2sHvj36LXGJr1ZnV0pJ6IoV3B82DMeKFSm6fj22+fMbKFrLJYTArWlTSuzaRc4vO/NkwwbutmxF7D//mDu0N0fkA9jcFda1Bisb+HI7tFoCLqZfvNSSJA9Fvx99nxkNZpDD4dXLhEkpidi6FYeKFbB/y/LLS6lE9RrI45yHKfWm4PvElxnnZmS5HanV8nDKVB7NnIVrs6YUWrkCa/fU602+3qwcHck7ejSFV61EFxeH3xftCZm/AKnRmDu015dOC2eXwsIacH03vDsG+pyC4g3NHZlF2HpzK3v99tK/cn+q5KlikDbjrlwh/uYtcrRsZZD2jE0lqtdE/YL16Va+G5t9N7P37t5Mn6+LiyNo8GCerFuHR9euFJg9Gyu7rE8izO6ca9em+PY/cfuoOY8XLsSvfQfi7xju1qqS5P5FWP4e7PkWClSFvqeh4XCwsTd3ZBbhRtgNpv09jdr5atO9fHeDtRv+++8Ie3vcPsoeIydVonqNDKg8gAq5KjDh9AQCIgMyfJ42PJx73boTdeAgeUaNJM/IEWlWM3/TWLu5UWDGDArMm4cmIIC7LVsRtnYdUqczd2jZX1ykvjbfssYQEQStV0DnP8Aze48oNaQYTQzfHvsWVztXptafipUwzHtSFxdH5M5duH7wAdauWStia2oW/2kkhHAQQqwWQpwRQngLIdJcbEUI8bEQYp8Q4pAQ4oIQYmCq/XuFEEdSfL1w7arsytbKlpkNZmItrBl6dCgJ2oR0z9ElJBDQrz9xV65QYO4cPLp0MUGk2Ytb0w8ptmM7TjVr8HDKFAJ69EDz6JG5w8qepISr2/S3+c4ugWrd9avtvvP5a1/lPLOmnJ2CX4Qf0+pPw8vRcJNxow4cRBcVRY7W2eO2H2SDRIV+hV8hpawFdAB+FUKkVQ57HjBSStkYaA6MF0K0SLE/WErZKMXXx0aP3Azyu+RnUt1J/Bv2L3POz3npsVJKgr8bR+z58+SfPg23pk1NFGX2Y5s7N4WWLCHvxInEXPThbuvWxJw/b+6wspewu7C+DWzuAs5e0OMgfDQbHNUaYqn9eetPtt/eztcVv6ZmvpoGbTvi963YFiiAU40aBm3XmCw6USUtK98DWAEgpfQFLgKd0jh8oZTyYtJxwcBhIGtLXWZzjQs3plPZTqz/dz0H7x184XGhy5YT8eefeA3oj1vz7HGv2pyEEORs15Ziv23C2skZ/y5dCVuzFrUCQToSE+D4bP2Q83un4cMfoOcRKFjV3JFZpDvhd5hydgrV81and4XeBm1bExTE09NncG/ZMlvd3rf0SIsDnsD1FNuuAc+VjpdSzk21yQEISetZUncAACAASURBVPGzsxBivRDiuBBimxCiosGjtSBDqg6hnGc5vjv5HUHRQc/tj9y/n5A5c3D76CO8+vY1Q4TZl33JkhTdshmXhg15OHUq978dji4mxtxhWSb/U7CkPhz8Hkp+AP3+htp936gCspkRlxjH0KNDcbB2YFr9aa9cIim1J5t+AyFw/+wzg7ZrbJaeqJJv8YWn2BYB5H7ZSUIIN6A6sCrF5tvAeCllffRXaCeEEPlecH6vpOdh3iEhIWkdYvHsrO2Y1WAWOqlj+NHhaHT/H14de+Uq94ePwLFiRfJNnfJaVZswFWtXVwrO/4lcgwcTuWuXvsyUv7+5w7IcT0NhWz9Y1Qw0MdDhN2i3FtxfvezP62z6uencCr/F1PpTye300o+5TNPFxPBk0yZcmzTRr5ydjVh6onqR9D5ZpwPfSyn/++SQUo6UUt5K+u8d6K/M0iw9LKVcKqWsJqWslitX9p1sWMitEBPqTODy48vMO69fEkTz8CGBffti7ZGTggsXYGWvhgFnlbCywuvr3hRaupTEhw+5+3kbog4fNndY5qXTwYW1sKAqXP4V6g2Bvmeh1Ifmjszi7bqziy2+W+hevjv1CtQzePvhv/+BLiICj65dDd62sVl6okoeWpXyaat7iu3PEUL0AjRSyoXptH0PKPJq4Vm+pkWb0q50O9ZcW8Pef7cR2KcvuuhoCv282KLL+mcnLvXrUXTrVmwLFSSwT19CfvrpzSxu++hfWP0RbO8PucrA1yegyQSwy15V9s3hRtgNJpyaQJXcVehfub/B25daLWG//IJjxYo4Vals8PaNzdIT1W0gDEi5Mlg54FxaBwshWgLvAoOTfi6Z9D23EKJrqsPzAPcNHK9FGlF9BJW9KvFo1Fjirl8n/5zZZl0u/nVkV7AARTdswL1lSx4v+pnAvv3QRj81d1imkRADBybA4noQ8i98ugC67obcZc0dWbYQHhfOoMODcLN3Y3aj2dhaGb5AbNShQ2gCAvDo1s3gbZuCRScqKaUOWAZ0h/8STyVgvRCirBDioBDCOmlfA2Ag0B9wEkK4AGOTmnIChgohnJOOrY7+GdYmU/4+5mJrbcv3V8tQ/bqWP5q6o61VydwhvZasHBzIN3UKecZ9R/SJE/h36oTmwQNzh2VcvvthUU04MRcqtIP+3lCl8xtdQDYzEnWJfHvsWx7FPGJuo7kGnS+VUtiq1dgWLIjr+02M0r6xZYdX0wRACCHOABuB9knDz92BMkDynx8bgUbAYyAq6Sv51l4wsA34SwhxDJgPtJRS3jTR72BW4du2EbtyHdoWTdhSKZbhx4aj1b2Bt6ZMQAiBR4cOFFqyBE1gIH5t2xF75aq5wzK8iCDY1Bk2tAEbR+i6Cz5bpJ8fpWTYTxd/4syDM4ytNZYKuSoYpY/YS5eIvXABjy87I6wNO4rQVISaA/Jy1apVk97e3uYOI8sS/P250+IzHCtWpPDyZfx+dzsTTk+ge/nuDKk6xNzhvdbifH0J/LoPiU+eUGDWTFzfe8/cIb06bSKcWwaHJoMuERoMgzqDwObNrQuZVXvv7uXbY9/SrnQ7xtYam/4JWRQ4eAhPT57krcOHsXZxNlo/qQkhzkspn5tKlBXZ4YpKySKp0/FgzFiErS35Z0xH2NrSulRr2pRqw8orK9nrl/nitUrGOZQqRdHfNmFfsiSB/QcQump19p4cHHQeljeGvSOhcC3oewYafKuSVBbcCLvBuFPjqJy7MiOqjzBaPwmBQUTt30/Odm1NmqQMTSWq19iTDRuJ8fYmz8iR2Ob5f9WpUTVGUSlXJcadHIfvE18zRvj6s/Hyosgvq3H94AMeTZ9O8ISJyMREc4eVOXERsGsYLHsPoh5Cm9XQcQt4FDN3ZNlSRHwEgw4PwtXWlTmN5mBrbbzVdZ+sXQNWVuTslFYxn+wjw4lKCFFOCDEhqbLDLSFElBDiiRDiuhDiLyHEICFE9ppF9hpLCAjg0ezZONerh3uqFTxtrW2Z02gOLrYuDDo0iIj4CDNF+WawcnSkwNw5ePbsSfimTQT0/hptVJS5w0qflHBlKyyoDueWQ41e0P9veLulKiCbRVqdluHHhvMw5iFz3p1jtMETANrISMI3b8GteTNs8+Y1Wj+mkG6iEkJ4CCG2AVeAXkAicAH4A9gJXEU/qm4c4CeEmCuEUPVRzEjqdDwY+x3Cyop8k75Ps/JELqdczGk0h+CYYEYcG6EGVxiZsLIi99BvyDdlMk/PnsW/QwfLHhEYdke/0u6W7uCaF3oeguYzwOHNWkjT0OZfnM+p+6cYU3MMFXMZt4pb+ObN6GJi8MyGE3xTe2miEkK4A8fRz2V6W0qZX0r5rpSyrZTySyllZyllayllXSmlJ/o5TGWA340fuvIi4b/9RszZs+QeMRzbfGlWiQKgUu5KjK45mpP3TzL/4nwTRvjmytG6NYWXL0PzIFi/GONNCxt4mhgPR2fCwloQ8Dc0nQ49D0MBw6ws+ybb57ePFVdW0KZUGz4v9blR+5IaDWFr1+FUsyYO5coZtS9TSO+KajAwU0rZXUr5b3qNSSlPSCmbAdeEENlnsZPXiCYoiEczZuJcpzY52rRJ9/g2pdrQumRrVlxZwY7bO0wQoeJcqxZF1q0FrRa/jp2IuXDB3CHp3T2un7R7eDKUbqZfJ6rW12DgwqhvoiuPrzD2xFgq5qrIyBojjd5f5N59JAYH49Gtq9H7MgU1PD0d2Wl4upSSgK96EOPjQ/Ht2zNceFKj1fD1ga+58OgCS99fSvW81Y0cqQL6EVkBPXqgefCAAnNmm2/4+tPHsH8sXNoIOYro14gq+b55YnkNBUUH0XFXRxxsHFjXfJ1Rn0uB/nPAr/Xn6OLiKL5zh9mW8zDb8HQhRPd09v/2auEoryJ8yxaenjpFnm+HZao6cvLgikKuhRh0eBB3Iu4YMUolmV3BAhTZuAH70qUJHDBQvwSDKel0cP4XmF8V/tkC9Yfqh5yrJGUwkQmR9DvQjwRtAgvfW2j0JAUQ8/c54q5dw6NLl2y15tTLZPa3GCOEcExrhxDiI0C9ws1E8+ABj6bPwKlGDXK0a5fp893t3Vn03iJsrWzpd6AfYXFhRohSSc0mZ06KrF6Fc726BI8fT8iChaaZa/Xwqn4Jjh0DIc/b+gKy741TBWQNSKPTMPTIUPwj/Zn77lxK5Chhkn7DVq/G2sMD9xafmqQ/U8hsosoJnBZC/HePImlU4FpgB+BmyOCUjJFS8mDceKRWS74pk7P8V1RB14LMbzyfkNgQBhwaQFxinIEjVdJi5eREoYULcf/sMx4vWKCfa2Ws6usJT+GvcbCkATz2hRaL9OWPcpcxTn9vKCklk05P4syDM0yoM8Hgy8m/SPydu0QfPkzO9u2xcnAwSZ+mkNlPtA3Ap8BYIcQKIUQn9Os6fQx8Dfxj4PiUDIj4YxtPjx8n9zffYFeo0Cu1VSFXBX6o/wP/hPzD6BOj0UmdgaJUXkbY2pLvh6n/zbUKHDQIXZyB/1C4sQcW1oSTP0LF9jDgPFTuqOZEGcHyf5bzx60/6F2hNy3eamGyfsNWr0bY2ZGzQ3uT9WkKmUpUUsr+Usp76KuZfwL8AvgD5ZIWG1RluU0sMSSEh9Om4VitKjk7djBIm+8XeZ9vqn7DX/5/Me/CPIO0qaRPCEHuod+QZ/Roog8eIqBHT8NMDI4IhF87wsYvwM4Fuu2FFgvAyePV21aes/vObn66+BPNizWnX6V+Jus3wd+f8N9/J8fnrbHx9DRZv6aQ2cEU+5MGVFwErIEhQBwwUAhhL4ToY4QYlZcIWbQIXUwM+b6fZNAHp13e7kLbUm1ZdWUVm303G6xdJX0eX3Ym/6yZxPj44N+lC4mhoVlrSJsIpxbAghpw66B+EcPex6BIbUOGq6Rw4eEFxp4cS5XcVZhUd1Kak+2NJeTHHxG2tnj1ef0+hjP7yfYe+vWhDqK/ivoJ/dIaD4HT/H/9J8UEEvz8CN+8hZxt22Bf3LB114QQjKo5inoF6jHlzBROBp00aPvKy7l/9BGFfl5Ewp27+HfoSEJgUOYaCPSGpY1g/xgoWhf6ndEvC68KyBrNvch7DDo8iPwu+fnx3R+xszbdv3XslatE7t6DZ7eu2OTKZbJ+TSWziSoaaJdUjeIhgNSbB7QB1M1uI0sIDOJm4/dICAzi0Y8/IuzsjPYXlI2VDbMazqJEjhIMPTqUG2E3jNKPkjaX+vUpvHIliU+e4N+xI/G3bv23L+Xr4Bmx4bDzG1jeBGIeQ9s10OE3yFnUtMG/YZ7EPaHvwb4ALHpvETkccpik3+TXwcOpU7HOmROP7i+dQZRtZWrCrxBinZTyhWV4hRCTpZQGvaoSQjgAi9GXZrIBRksp97/g2KFA8oOaX6WUM1PsKwqsSmpDC3STUt5Nr39Lm/Dr3707MWfO4lC+PHGXL+PVtw+5Bg40ap/BT4PpuLsjWp2Wn5v8zLXQa/hH+VPEtQhNizXF2Tb7Lh+QHcTduMG9Hj0gQUOhpUtwrFgR/27diTlzBufatSm8coW+gOw/W2DfaH2Cqvk1vDsa7F3NHf5r6anmKXvv7sU/yp88jnn449Yf+EX6seyDZVTOXdlkcfh3707M6TMgJXlGj8Ljyy9N1nd6DDnh1+IrUwghpgH5pJRdhBClgDNA2eQruhTHNQXmoV+qHsAHGCql3JW0/zSwVEq5SgjRDegjpayRXv+WlKiijx4lcPAQZGwsWFlh5eTEW0cOY+3iYvS+74TfodPuTkRrorGztiNeG4+jjSNWWLGoySKq5FG14IwpISCAe92/IjE0FK+ePXi8+GdkvAZhb0uh6aNxfrQR7hyB/FXgk3mQz7gFT99kFx5eoO+BvujQEZsYixVW6NAxpMoQur9juiua6KNHCRg0GOLiQAgKLphvUYtzmqwyhRBilBCidWYbTTrvlcdkCiGsgB7ACgAppS/6gRxpXdX1BjZIKeOklHHAevRD5hFCVAQqJ20j6fs7Qoiqrxqjqeji47k/9jt9kgJ9VQEpEbbGW8smpTzOedDoNEgk8dp4AGITY3ma+JS+B/oSo4kxSRxvKrtChSiyfh22BQoQ8uNP+iRlJfF8KwzHk18jg85D81nQ44BKUkb0VKN/vT9NfEpsov69qEM/hWPp5aUmex8kfx6QPIVBSoInTEQXH2+S/k0tveU4lgEnhRANgdlSSv+XHSyEqAyMQT/xd7oB4isOeALXU2y7BqSVpasDG1Md1zfFvjtSygQAKWWCEOJWUjvnXxbAjRs3aNSo0TPb2rZtS9++fYmJiaF58+bPndO1a1e6du3K48eP+fzz56sk9+nTh3bt2hEQEEDnzp2f2z906FA++eQTbty4Qe/evQF9sVnNgweg09Hb04s6zs5cjY6ie/ny2BZ4tlzS1KlTqVOnDqdOnWL06NHPtT9v3jwqVarEgQMHmDx58nP7lyxZQunSpdmxYwezZ88GICQ2hIDIALRSS8FeBbHztCPibAShh0KxFtbU+KnGf+VhtmzZgpeXF6tXr2b16tXPtb97926cnJxYtGgRv/32fNmgI0eOADBr1ix27tz5zD5HR0f27NkDwKRJkzh48OAz+z09Pdm6dSsAo0aN4vTp08/sL1iwIOvWrQNg8ODB+Pj4PLO/VKlSLF26FIBevXrh6/vswpKVKlVi3jz9kP1OnToRGBj4zP7atWvzww8/ANC6dWtCU43Ye++99/juu+8AaNasGbHJf3gk+fjjjxk2bBjAc687gKY5c9JCCGK1iQwI9cc2Rov2lBU6Zw9s962na1dbg7/2Uho7dixNmjTBx8eHwYMHP7ffGK+9lNauXUuhQoXYtGkTP//883P7jf3a6/FjD3ToePTnI6KvRT+z397Nnr019tKqZCujv/ba1q3Hvev/6v9gTd7/5AlTli4j14D+Rnntvern3qt46RWVlPIx0BD986E7Qgg/IcQBIcSvQojVQohfhBC/CSGOCiGCgXNAKPCplAaZKZq8LG14im0RQO4XHPui41Lve1k7CCF6CSG8hRDeGo0m00Ebg+bhw2delADEJ+i3m0B8YjxamXa1BK3UqioWJhJ/yZu81cJwzhuP5qk1saG2xEfYoHn02NyhvRH8o/z/u5JKTavTci/ynkniSLh797nPA5moIWztWpP0b2oZfkaVdJvsc/TD0fOi/5DXoh+aHgjsAbZIKQ1W0VQIURc4AdgnXw0JISYDdaSUjVMdqwE+Sh5oIYRoAuyTUloLIcYk7auT4vgTwH4p5fcvi8FSnlE9+mk+oStWQIpLe+HoiGf37uQa0N/o/W/13cr0c9PTfJM6WDswquYoWpVUK7sYjU4HF35Bt2M4QiYQctWF0GuuIAUIiefn75F70kJzR/na23JjC1POTiFRJj63z9HGkZE1RprkffBw1mzCli9/ZpspPw8ywlzV0y+hv/1XW0pZTErpLKV0k1KWTFpMcYYhk1SSR0nfU471dE+xPfWxqY8LecG+l7Vjkbx693rueZS1iwuevXqapP+mxZpi9YKXi0an4f3Cqh6x0QT/Ays/hJ2DETYCv0OehF510ycpACnQ+BxGJiSYN843wMPYh2kmKQArrGhatKlJ4hC2zz+1MeXngallZCl6KyHEDPS3yh4KIR4LIYw7Hvr/bqNfXbh0im3l0N9iTO3cS447BxQTQtgBJH0v8YJ2LJLUJOpnuSdVnxAODuSbMhkre3uT9O9s68yiJotwtnHG0UZfQN/RxhE7Kzu0Usus87NMU/X7TRIfDfvGwJKG+qXha/RGaBPIVS4aYZ1028dKh1uRGCJvSgJ7tH9tH6Zbgo3XN7L40mLq5K+Dk7XTM+8DZxv9+8PJ1vjV5xNDQ3nyyxocq1RBOOpjMPXngamlN5gC9CPnhqGvPHEPKAXMFUI8kFIatbaOlFInhFiGvrbgcSFESfTDzzsKIcoCC4APpJRa9HOt5gghkudOdUiKGymljxDCB2iPvj5he+CqlPKlAyksSdiqVeiiorB/5x3ir17FqWpVXBo0MGkMVfJU4VDbQ+z128u9yHsUditM06JNWXllJUsuL8HJxonh1YebtGzMa+v6Ltg9HCIDoUoXffmjcH+wtsEFcAw+QczNEJxL5qVA37o4HrvNw62XCejZi4KLFppkysKb5M9bfzL17FQaFWrEnEZz0Gg1z70PTJGkAB7/vBhdfDz5Jk8mePIkYs6cNcvngUlJKV/6BVwGPk+1bQDwV3rnGuILcABWo58/5Y0+MQHUAoIAhxTHDk065jzwbap2igKHgOPAEaBYRvqvWrWqNDdNSIj8t3IVGTBosIwPCJS+7zaW8QGB5g7rPzqdTk47O02WX11eTj0zVep0OnOHlH098ZdywxdSjneTcmFtKf3PpHlYWq+D8O3b5bW3y8s7rVpLTWioqSJ+7W313SrfWf2O/GrfVzIuMc6sscT7+8tr5d+R98eN1/9sgZ8HyQBvaaA8kO5gCiHETSllyVTbrNBfkZR9pSyZDVjCYIrgSZN58uuvlNi1E7uiRc0ay4tIKZnlPYs119bQplQbxtYai5V4PVYXNQmtBs78DEf0Q9tpNBJq9QXrzM2Tiz56lMBBg7HNm5fCK1dgmz+/EYJ9c2y6vonJZydTN39d5r07Dwcb867xFDR0GFEHD1Ji3z5s86Q5aNlimHowRXTqDVI/9Py58chCiNUGiElJIeHePZ5s2kSOtm0sNkmBvojtsGrD6PFODzb7bmb8qfFodUZa/O91c++s/jnUX99BsYbQ7yzUHZTpJAXg0rAhhVcsJzE0FL8Oz9YHVDJn3bV1TD47mYYFG/Jj4x/NnqRivL2J3LULjy+/tPgkZWgZSVQvuuRKa3uFV4hFScPjnxcjbG3J1bdv+gebmRCCgZUH0rdiX7bd2saYk2NI1KU9QkoBYsJgxyBY+QHERcAXG6DDr5Cj8Cs161S1KkXWrkFqE/Hr2ImYCxcMFPCbY+WVlUw/N50mhZswt9Fc7K3NO0hBFxvL/TFjsC1QAK/evcwaizlkZDBFJSFEWn8aixdsVwxE8/ARETt3krNdu2xTul8IQZ9KfbCxsuGniz+h0WqY1mAatlamKfWULUgJlzfpR/TFPoHa/aHRKLA33AAIhzJlKLpxIwFf9eBet+4UmDsH18aN0z9RYfGlxSz0WUizos2YUn+KRbx2Q378CY3/PQqvXoWV85tXBDojieoJsD0Dxwn0S9IrBvJk/XrQavHoYjkVkTOqZ4We2FnbMct7FolHEpnVcBa2WbiV9dp5fBN2DgG/41CgGnT+A/IZ50aEXcGCFNm4gYDeXxPYfwB5J04gZ5s2RunrdSClZIHPApZeXsonxT9hUt1JWFtZmzssYi5cJOyXX8jxRTuca9UydzhmkZFEdU9K2S0jjQkhLr5iPEoSXUwMTzZtwrVJE+wKFTJ3OFnS5e0u2FjZMO3vaQw+Mpg5jeaY/RaK2Wji4PhsODkPbBzhozlQtdt/8+KMxcbDgyKrVxE4ZAjB340jMSQErz591BSCVKSUzL0wl1VXVtGqZCvG1RpnEUlKFxfHgzFjsMmXl9zDvjV3OGaTkXfJB5loLzPHKi8R/scf6CIi8OjW1dyhvJKOZTvyXa3vOBZ4jIGHBr6ZVdZvH4Kfa8OxGVDuMxjgDdW/MnqSSmbl7EyhhQtx/+wzHv80n+CJE5Faddc+mU7qmHFuBquurKJd6XaMrz3eIpIUwOMFC0i4e5d8kyZh7fLm3fJLlu4VlZQyJL1jsnKs8mJSqyXslzU4VqyIU2XTLcJmLG1Lt8XWypYJpyfQfV93Fry34L9K66+1qIf6hQyvbAGPEvDln1C8kVlCEba25PthKja5vAhdthxtaCj5Z816bSsZZFSCNoExJ8aw128vncp2sqgJ67GXLhG6chU52nyOS9265g7HrNREFwsUdegQmnv38OiWoTuu2ULLki2Z12get8Nv03l3Z/wi/MwdkvHotPD3MlhQHf7drh8o0eeU2ZJUMiEEuYcOJc/oUUT9dYB7X32FNiLCrDGZU0R8BL3/6s1ev70MqTrEopKULiGB+2PGYJM7N7mHDzd3OGanEpUFClu1GtuCBXF9v4m5QzGodwu/y8oPVxKTGEPnPZ3xeeST/knZzYNLsOJ92D0M8leCPqf1k3dtzTsHJyWPL78k/+xZxF66jF+HjiSkWlPrTXA/+j5d9nTBJ8SH6fWn0718d4tJUgCPFy4i4dZt8n0/EWtX1yy3c/fuXRo1aoSDgwNFixalUaNG1KpVizJlyjB8+PD/1qSKjIyka9euVK9enSpVqjBixAgSEy1oaomhSly8rl+mLqEU4+Mjr5UuI0N/+cWk/ZqSf4S/bL61uay6tqo84HfA3OEYRlyklHtGSTkhh5QzSkh5aZOUFl5KKvrMWXm9eg15o3Yd+fTCBXOHYzL/hv4r3930rqy9vrY8e/+sucN5Tsw/V+S1cm/LoJGjDNZmkSJF5Pjx4//72dfXV7q6usq+fftKKaVs2bKl7NSpk5RSyvj4eFm7dm05atSr9Y8BSyipKyoLE7p6NVaurri3am3uUIymsFth1jZfS2mP0gw5MoT1/643d0hZJyVc2w4LasCZRVC1K/Q/BxXaggX9hZ4W55o1KPrrr1i5uHCvS1cik1ZOfp2dCjpFlz1dsLay5pdmv1AjXw1zh/QMmZDAg9GjsfH0JM/IEUbrp2TJkjRu3JjNmzdz5coV/vjjD4Yn3WK0s7Nj8ODBzJs3j+jo5woTmYVKVBYkITCIqH37ydmu7Ws/wsfDwYPlHyynUaFGTPt7GrO9Z6MzyKLQJvTEHzZ+Ab91BidP+Oov+HguOOY0d2QZZl+8GEU3/YpD+fIEDfmGx4uXJBdxfu1su7WNfgf7UdC1IOuaraNkzpLpn2RijxcvId7Xl7wTJ2Dt7m7UvjQaDXZ2dhw4cAAHBwfKly//377q1asTGxvLiRMnjBpDRmVkHpViIk/WrgErK3J26mTuUEzC0caRuY3mMu3vaay+uprgp8FMrjfZ8udaaTVwegEcmQ7CCj6YAjW/Buvs+XayyZmTwqtW8mDMWELmzSPB3598Eycg7OzMHZpBSClZcnkJC30WUitfLeY2mouL3atXAZm44yrX7kcaIEI9XUwMcdfApsVY7HwdwPf0c8eUy+/G+E/efuW+Tpw4wYEDBxg5ciR37twhT548zzyjy5s3LwB37hh6LdysyZ7vrNeQNjKS8M1bcGvWDNukF8mbwNrKmtE1R5PfJT9zzs/hfvR9ZjeaTV5nC/038D8Nu76BR9egzMfQdBrkyJ4TslOysrcn/8wZ2BUpwuOFC9EEBVHwpx+xzpF6YezsJUYTw4TTE9hzdw+flviUCbUnWGSFFKnREH/zJsLGBrvCr1br8UVWr17NkSNHiI+Px87Ojrlz59K7d2969eqFfappCsk/x8RYxrxHlagsRPjmLehiYvDM5hN8s0IIQbfy3SjkWogxJ8bQbmc7ZjWcRfW81c0d2v/FhMFf4+DiWnArqC8gW+Yjc0dlUEIIcg3oj12RwjwYMxa/L9pTaOkSo31wGtu9yHsMPjKYW09uMbDyQHq808OgI/sMcWUDoIuP596XXYi7cYMi69bhWN4w7abWtWtXJkyY8Nx2Z2dn4lOtDJ38s5OTaRaDTI/FPqMSQjgIIVYLIc4IIbyFEC+seiGE+FgIsU8IcUgIcUEIMTDV/r1CiCMpvnYa/zfIOKnRELZ2LU41a+JQrpy5wzGbJkWasPGjjbjZudFzf0/WXltr/uclUoLPBlhQTf+9zkD9MhyvWZJKyf3TTym8aiXa8HD82rQl+sRJc4eUaccCj/HFzi94FPOIxU0W07NCT4safp5MSsmD0WOIvXSJ/DOmGy1JvUzx4sV59OjRM++14OBgAEqUKGHyeNJisYkKmAAIKWUt9MvK/yqEyPOCY+cBI6WUjYHmwHghRIsU+4OllI1SnbfEWQAAIABJREFUfFlU8dzIvftIDA7O9uWSDKF4jv+1d9/hUVXpA8e/JySTEFKAEFogCb0ISFUU6UXpioC6+FtRERRddFWQdbGhq2BDpNrDsgiKDQEpCRJDaBKKSq8SWgghENIzyby/P+6AA6RnJjMTzud55pm5Zc5952Tmvrn3nntOQxYPXEz3et15e9vbvLDhBed1u3TuAEQMgh+eMHqWGBcD/V63ay/nrsq3Y0fCv1qCZ82anHjsMaORhcX1G7tYxMK8XfN4at1ThPiHsGTgEm4Pud3ZYRUoae5cLq1cSfA//0lAP+f0QNe3b18yMzPZs2fPlXlxcXFUrlyZLi7SI4ZLJirrCMJjgM8AROQgsBMoqJXBHBHZaV03AViPm/Q7KCIkf/EFpgYN8OvWzdnhuAQ/kx8zes7g6fZPs/rYah5c9SDxl+LLLwBzJqx7HeZ1gbO7YfBMeGQN1G5V9HsrEFNYGOFfLSFgwADOffABJ/8xgbzUVGeHVaBLOZeY8PME5v42l0ENB7Gw/0Lq+ddzdlgFSlm5kqRZswkcOpSgsY85LY6bbrqJe+65h3feeQcwWgPOnDmTZ555Bj8/F/mnzF43ZNnzATTGGJixps28WcDiYr5/BfCqzfRSYBGwAfgBuLm4sTj6ht+0LVtlb7PmkrzkK4dux13FnoyVLou7yG2LbpNfTvzi+A0ejBT5oI3IKwEi344VSU10/DZdnMVikfMLFsjem1rJ4TvvkqyDB50d0nUOJB+Q/t/2l7YL2sqX+74Ui4vfbJ2xa5fsa91Gjv1tlORlZztsO0ePHpXu3buLt7e3hIWFSffu3eXixYvXrZeSkiJ///vfpWPHjtKuXTuZOHGimM3mMm0bO97w6/SklG9Q0MWaqEw2894A1hXjvQHAWSDMZt40oLH19WAgFahTSBljgTggLjQ0tJR/puKJf/wJOXDb7ZKXmenQ7bizE5dOyPAfh0vriNYye+dsMeeV7QeUr5TTIl8/ZCSoD9uLHC2HpOhm0rdtkwNd7pB9bdtJysqVzg7niuVHlkun/3WSHl/1kB1nXb+HjZxTp+RAlzvkUJ++Yk5OdnY4DmPPROWSp/4KUZyrodOBqSJy/PIMEZksIoetr5cDe4ECRyMUkY9FpKOIdAx24Mi65tOnSYuOpup9I/HwcZ2+4FxNPf96/Lf/fxncaDDzf5vPQ6sfst+pQEsebP0Y5twC+3+Cnv82OpBtoE/DXsu3Y0cafPstPs2bc+rZ5zj71jTEbHZaPCnZKUz6ZRL/2vAvWlRvwdeDvqZdTdcebSAvLZ0TT4xHsrKoP28untXc5+ZwZ3JKolJKrVJKpRXwWAskWle1vYkj0GZ+QeWOBcwiMqeIEOKBsFJ/ADu5+MMPIELVeytud0n2UtmzMv+54z9M7zqdYynHGL58OEsPLr18BFw6p3fCJ71g1UQI6QDjN0P3SeDp4jccO5FXrZqELYig2oMPkrxgAfGPPIr5bKE/S4fYdHoTw34cRuTxSP7R7h98dudnBPs67p9Ke5C8PE4//zzZhw8TMmMG3o0bOzskt+GURCUi/UXEr4BHP+AIkAw0s3lbS2BbQWUqpe4BegLPWKebWJ9rKqVGX7N6LeC0/T5RyYnFQsr3P+DbuTOmeq57wdfVDGg4gO+GfEeb4DZM3TyVp35+iqTMpJIVknUJVr1gJKnUMzD8c2NI+CDXaIrr6pTJRO0p/6bu29PJ/OMPjg0ZUm79BGblZjHt12mMixxHFa8q/G/g/xjbZiyeHq59S6iIcPataaRFR1Pr3y/i1/UOZ4fkXux1DtHeD4zrSl9YXzcBzgO1rdMtgHVAJet0N4yWfkGAn/WxwLosHPgDqGKd7gRkAk2KE4ejGlNcbkRxcdkyh5Rf0eVZ8uR/e/8nHRZ2kK6Lu0rU8WL0wm6xiOz+XuTdZiKvBIqseE4k44Ljg63Aso4claMjRsreZs3l5HPPS24+F+rtZXfSbhn8/WBpFdFK3tr6lmSa3eO6riU7W05NmiR7mzWXhDffcnY45YaK3pjC+Iz4ABHAFoyGDf1slnUGTgE+1ulTGI0vbB/RNuW8DmwCYqzl3VXcOByVqE5NekH2d+goeRkZDin/RnH4wmEZ8eMIaRXRSqbETpHU7NT8Vzx/VGThvUZjiXl3iJyMK99AKzCL2SyJs2fL3pY3ycFu3SVt40a7lm/OM8tHv30kbRe0lV5f9ZKNp+xbviPlpqbK8Ycflr3Nmsu5uXNdvjWiPdkzUSmjPK0gHTt2lLi4OLuWmZeWxqE7uhI4ZAh1pr5m17JvROY8M/N+m8dnuz+jTpU6TOk8hTtCrKdWcnNg8yz45W3w8IReU6DTY27bgawry/xjN6dfeIGco0ep9uCD1HzuWTwqVy5TmQcvHGTq5qn8du437gq/iymdpxDo7dhexe3FfDaRE+PGkX34MHWmTqXqsHucHVK5UkptF5GOdilLJ6rCOSJRXVi6lISXXib8qyVUvvlmu5Z9I9uVuIspG6dw/NJx+oT2YVKdXtSJeh2SDkCLIdB/OgTUdXaYFZolK4vE997nwsKFmBo2pO706VRuXfIbpdNy0pizaw6L9y/G3+TP5FsmM7Ch+3RblX34MPFjx2K5mELIzJk35DUpnajKkSMS1Z/3P0BeaioNVyx3yf7H3FlOXg4Lds7j4z2foyy5jM3y4KGe0/CqwH3zuaL0TZs4/eK/yU1KIujRR6kxbiwexejgVET46dhPvBv3LuczzzO86XAmtJtAVR/36cU9Iy6OE+OfRHmbqD9/PpVvKv/++1yBPROVu91H5fayjx4lc9cuqg4bppOUvVksmH77isfWzWDZyTPc5hvCTF8Ytm8+m09fP7aP5jhVbr+dhst+IHDgAM5/9BFH+g8g5ccfC+0v8MjFIzy69lEmb5hMLd9afDnwS16+7WW3SlKXVq8m/uFH8AwKInzxkhs2SdmbTlTlLOW776BSJQKHDHZ2KBVL4j6IGAg/PgU1mlF3TDQz71vLnN5zyJM8xkaO5flfnichPcHZkd4wKgUGUnf6dMK+XIRnjRqcnvQCxx/4G5m//37VehnmDN6Pe5/hPw7nQPIBXur8EosGLKJVDffpW1FESF6wgFP/fBafVq0I+3IRpnohzg6LY8eO0aNHD3x8fAgPD6dHjx507tyZ5s2bM2nSJDIzM6+su23bNlq2bMno0aOdF3AB9BXlciS5uVxctgy/7t3xdGCPFzeUnAyIeRs2zQJvfxj8IbT7P/Aw/gfrVq8bt9a5lc93f85nf3xGzMkYxrQew6gWo6jiVcXJwd8YfNu3J3zp16R8/wOJH8zgz5H3ETh0KNX/OYHVab8ya+csEjMSGdZkGE+3f5rqPtWdHXKJ5Pz5Jwlv/If02Fj8+/ah7jvvuExPMw0aNCA6Oprw8PCrxqM6dOgQHTp0ID09nTlz5jB9+nQiIyNd9iyPPqIqR2kbNpB3LumGa/3jMAfXwtxbIXYGtLkPnoqDDg9dSVKXeVfy5ombn+D7od/TuU5nZu2cRf9v+/PF7i/IzM0soHDNnpSHB1XvHUajVaupNuZRLq5cwf5+fdjxzr+p7RnEwv4Lee3219wqSVmysjj34SyODh5C5s6d1HrxX4R88IHLJKnCNGnShF69erF06VIAmjdvztq1a3Fkl3FloY+oylHKd99TKSgIv+7dnR2Ke0s5Basnw74foUYzGL0SwotuVVXfvz4f9vqQ38/9ztxdc3l/+/ss2LOAR1s/yoimI/DxdP0djDuziIXIpA3Ma7yB9DHweKwfD/ySiufBJKqn7CLv3kZU8vd3dpjFkvbLLyS8/gbmkycJGDSImpMm4lWzprPDKhGz2YzJZAJg6NChRaztXDpRlZPc5GRS16+n+oMPory8nB2Oe8rLhW2fwM9vgCUXer1kjLjraSpRMW2C2zC/73x2nN3BnF1zeHvb20TsjmBMmzHc2+ReTJVKVp5WOBHh5/ifmfPbHA5dOETDwIZMvPs9+j7dl8ytv3Ju1mwSp00n6cNZBA4bRvUHR2EKD3d22PkynzpFwltvkRa1jjo9wH9wDSoF7IWVD5dPALVbQ/9pZS4mNjaWqKgoJk+ebIegHE8nqnJyaflyyM0lUJ/2K52T22HFM5DwOzTuAwPeheoNylRk+1rt+ezOz/j1zK/M3jWbN7e+yRe7v2BM6zEMbjSYyp5lu1n1RpdnySP6ZDQf/fYR+5L3ERYQxrSu07gr/C4qeVQCoErnzlTp3JnM3Xu4sPC/XPjqKy4sWoRfjx5U//v/4du5s0tcN7Hk5JAcsYCkuXNBKYKfe5bAugdRiXudHVqxRUREEB0dTXZ2NiaTiRkzZjBu3Dhnh1Us+j6qItjjPioR4djQu1He3jRY+rWdIrtBZKUYo+1u+xT8asFdb8FN94Cdd14iwubTm5mzaw6/J/2Ov8mfuxvfzX3N7iMswOkd7buV5Kxkvjv0HV8f+Joz6Weo51ePx29+nIENBxbZeaw5MZGLS5ZwYfES8i5cwLtJE6o/9Hf877yz3E8LisVCRlwcl1as5NKaNVhSUvDv24da//oXXnXd68bxaxtTFKRHjx6Eh4cTERFR5m3a8z4qfURVDrL27CX74EFqv/qKs0NxHyKw5ztY/S9IPwe3jDW6P/IJcMjmlFLcHnI7t9W9jR2JO1iyfwmL9y1m4d6FdKnbhfub30/XkK5XjgS0q4kIvyf9zpL9S1jz5xrMFjO31r6ViZ0m0rN+z2L3bu5VsybBEyYQNG4cl1asJPm//+XMlJc488qrVL75Zqp0uR2/O+7Ap1UrVCX7/y1EhOz9+0lZvoJLP/1EbkICytcX/969qXrvMKp07mz3bWpF04mqHKR89y3K25uAAQOcHYp7OH8EfnoejvwMddrCA0sgpH25bFopRYdaHehQqwPnMs7xzaFv+ObAN/zj538Q4hfCiKYjGNZkGNV89IB3AJm5maw+tprF+xezL3kfVbyqMLzpcO5vdj8NqzYsdbke3t5UvXcYgcPuIXPnTtJiYkiP3UjS7DkkzZqNR0AAVW677UriKssRjiUjg5wTJ0j7+WdSVqwk58gR8PTEr2tXAiY+j3/PnsXqVUNzHJ2oHMySnU3KipX49+1LpQDHHA1UGLnZsPFDiHkHKpmg/zvQ6VFw0lFMsG8wT9z8BGNaj2F9/Hq+OvAVH+z4gLm75tIlpAt9w/rSvX53Akw31t81Oy+bzac3E3k8kvUn1pOak0rjqo2ZcusUBjUaZNf705RS+LZvj2/79vDMM+ReuEDG5s2kbdxIeuxGUtesAaBSjRp41qiBZ1AQlYKq4xlUA88aQVQKCsIzKAiPKn7kJiZiPnMG8+nTxuPMaXJPnSYvJeXK9nw7dqT6q6/if2c/PfquC9HXqIpQ1mtUKStXcvq55wn9/DOq3H67HSOrYI5tgJXPQtJB4xrUnW9BQB1nR3WdIxeP8M3Bb4g8HsnZjLN4enhya51b6Rval16hvSrskVaGOYONpzcS+WckMadiSDen42/yp2f9ntzT+B461OpQ7o0eRIScI0dIi40l+/Bh8s4nk3v+PHlJSeSeP49kZ+f7Pg9fX7xC6uJZty5edeviVcd49m3fzu2uPRXl2LFjPPzww2zZsoXatWsTHh7OsmXLCAy8ugf6jz/+mC+//JJdu3bh4+ND8+bNefbZZxkyZEipt607pS1HZU1U8Y+OIfvYURpHRaE89P3V10k7B5EvwW+LoVo4DHgPmvRxdlRFsoiF3Um7iToeReTxSE6mncRDedCpVif6hPWhW71u1PVz753ehawLbDq9iajjUcSeiiUrL4tq3tXoFdqLvmF9uaX2LXhVcs1bLUQES3o6eefPk3v+PJbUVDxr1sSrbl08AgJcoiVhRXfDJCqllA8wH2iOcZryRRFZW8C6ozGGob9oM3uKiMRal1fHGIixurWs8SKyo6gYypKozGfOcLhXb2o88QTBE/5RqjIqLIsFdv4XIl+BnHTo8jR0ex683K9JuIiwP3k/kccjiYqP4ljKMQBq+daiXc12Vx5NqzV12cYYIkJ8ajw7zu5g17ld7EzceeVzBFcOpndob/qG9aV9rfYuP+y75hpupFZ/r2Ik085KqabAFqVUCxE5W8D6z4hIdAHL5gI7ReQVpVRPYJlSqrGI5H9+wA4u/fQTiBB4t2vf9V3uzu6BFf+EE1sh7A4Y9D4EN3N2VKWmlKJFUAtaBLVgQvsJHLl4hK1ntrIrcRc7Enew+s/VAPh6+nJz8M20q9mONsFtCA8Mp7Zv7XJPXiJCUmYS8anx7E7azc7EnexM3ElyVjIAgd6BtA1uy5BGQ+hYqyNtgtvgofTZAM15XDZRKaU8gDHAMAAROaiU2gk8CLxXwrKqAyOABtay1iulcoBBwLf2jNvWpTVr8bnpJkyhoY7ahHvJSYfoabB5DvgEwt3z4OYH7H5PlLM1qtqIRlUb8bcWfwPgTNoZdiTuYGfiTnYl7mLeb/MQjDMZnh6e1POrR33/+oQGhFLfv/6VRzXvalT2qozJw1SiU1Vmi5kMcwapOamcTDvJidQTnLh0gvjUeON16omr+jis71+fO0LuuHLk1yCwgU5Mmktx2UQFNASCgP028/YChR1KjlNKvQYI8LWIzLXObw9ki0h8PmU5JFGZT50i6/ffCX72WUcU734OrIKfJkLKCaN3875Twdd9OiAtizp+dRjoN/DKCLWpOansT95P/KX4q5LH9rPbycjNuO79lVQlfD19qexV2Xj2rIyvly8eyoNMcyYZuRnGw5xBZm4mZov5ujJMHibq+dcj1D+UW2rfciUpNqvWjGBf1+yIVNMuc+VEVcv6bHvNKQVoWcD6Z4G1GNehgoFopVQlEZllLeviNeunAPn2IqmUGguMBQgt5dHQpchIAALu7Feq91cYKSdh1QuwfwUEN4eHV0HYjd360d/kT6fanehUu9NV80WE5KzkK4nrUs4lMnMzyTAbicj2dYY5gzxLHlV9qlLHsw6+nr74evle9VzFqwohfiGEBoRS07emPkrS3JYrJ6qC5HsORERW2UwmKqVmAROAWaUo62PgYzAaU5QmyNQ1a/Fu0QJT2A3a/U5eLmydD+vfBLFAn1eh85Ml7kD2RqKUIqhyEEGVg2hbs62zw9E0l+G0RKWUWgV0LWDxJuBJ6+uqQKL1daDN66LEA5ezRKL1vbYCgT3FLKtEzGfPkrlzJ8HPPO2I4l3fyThY/gyc/QOa3AkD3oFqN2jC1jStzJyWqESkf2HLrY0pkoFm/JWcWgI/FbD+JBF522ZWLeC09fUOwEcpVV9ETtiUFVG66AuXusZoQe/f705HFO+6Mi/AuqkQ9wX414GRC6HF4ArXWELTtPLlsietRcQCfAI8AqCUagK0BRZZp1sopdYppS637e2vlOpuXeYLjAMWWss6Dyy1Kas7YAJWOiL2S2vX4N2kCd4NyzYMhdsQgd+XwuxOsD0COj8BT/0KLYfoJKVpWpm5+jWqV4H5SqktGLE+ICIJ1mWBGDcCewF5GE3WX7YeifkB64C3bMoaD0QopWKBSsBQEcmyd8DmxEQyt++gxpNPFr1yRZB02Oj66NgvULc9jPoG6urrK5rmCvLrQikrK4uLFy8yZMgQXnvtNS5dusTMmTOJiYnB09OTlJQURo4cycSJE/H0dI0U4RpRFMCaSEYXsGwLEGIzvQJYUUhZyUDpO64qptSoKBAh4K4KftrPnAUbP4AN74GnjzGQYcdHnNaBrKZp12vQoAHR0dHXjUd16NAhOnToQHp6Oi1atGDNmjVER0fj7+/PiRMnaN++PZmZmUydOtW5H8DKZU/9uavUNWsxNWqEd+PGzg7FcY5Gw/wuEP2WcQ3qqW1wy2M6SWmam2jSpAm9evVi6dKl1KxZk4kTJ+JvHZiyfv36jBw5kkWLFjk5yr+49BGVu8k9f56Mbduo8bh7DO9cYmmJsObf8MfXUK0BPPgdNO7t7Kg0TSsFs9mMyWRi5MiR1y2rXLky2QX0Pu8MOlHZUWpkFFgs+N9ZwU77WSywIwKiXoWcDOg2Cbo+65YdyGqavUz/dTr7k/cXvaIdNa/enBdueaHM5cTGxhIVFcXkyZPzXb558+Z8E5iz6ERlR6lr12AKC8O7aVNnh2I/CX8YHcie3AbhXWHg+xBcgT6fpt0gIiIiiI6OJjs7G5PJxIwZMxg37vqzP1FRUcTHx7NypUMaRZeKTlR2knvhAulbfyVozJiKMdZNdppxDWrLPKhcDe75GNqM1M3NNc3KHkc25cm2MUVBjh8/zvjx4/nxxx+pWrVq+QRWDDpR2UnaunWQl1cx+vbbtwJWTYJLp6DDaOj9yg3Tgaym3agSExMZOnQon3zyCe3atXN2OFfRrf7s5NKatXjVr493ixbODqX0LsbD4gfgq1HgUxUeWQuDZ+okpWkV3IULFxg4cCBvvfUW3bt3B4zh6V2FTlR2kJeSQvrmzQTc2c89T/vlmWHjTJhzq9H0vO/rMO4XCL3V2ZFpmuZgaWlp9O/fnz59+hAcHExcXBxxcXF89NFHzg7tCn3qzw5S1/0Mubnu2dovfqvRWCJxDzTtDwPehqp6oEdNqwgu90yRkJBwpTHFsmXLCAz8q4/umTNnsnXrVrZu3cq0adOcGG3BlEipRrG4YXTs2FHi4uIKXefEuMfJPnSIRuui3OeIKiPZaG6+YwEEhED/t6HFIGdHpWlaBaGU2i4ihQ10W2z6iKqM8lJTSdu0ieoPPugeSUoEfv/KuHE38wLc9hT0+Bd4+zk7Mk3TtHzpRFVGaevXg9nsHq39kg4Zp/n+3AD1OsGgH6B2a2dHpWmaViidqMro0pq1eNaujU+bNs4OpWDmTNjwvtGJrFdlGDQD2o8GD92WRtM016cTVRnkpaWTvmEDVe+/D+WqO/3D6+Cn5yH5KLQeCXf+B/xqOjsqTdO0YtOJqgzSoqORnBwCXLG1X2oCrHkRdn8L1RvB//0AjXo6OypN07QSc9lEpZTyAeZjDI7oCbwoImsLWHc14GMzywtoLSIBBSxPE5EyN3FLjYzEMziYyq50F7clD+I+N4aEz80yGkp0eQa8fIp+r6Zpmgty2USFMbqvEpHOSqmmwBalVAsROZvPukdFZPzlCaXUPVw9SGKCiIy2Z3BiNpO+cSMB/fu7zmm/M78ZjSVObYcG3Y0OZGtU4HGxNE27IbhkorIOJz8GGAYgIgeVUjuBBzGGnL+KbZKyGgXMc2SMGTt3YklLo0q3ro7cTPFkp8L6N2HrfPANgmGfQuvhugNZTdMqBJdMVEBDIAiwHexlL1DkzWNKqUCgLbDeZnYVpdQiIBQ4D7wiIr+VJcD0mBjw8qLKbbeVpZiyEYF9y2HVC5B6Bjo+DL1fNno71zRNqyBcNVHVsj5ftJmXArQsxnvvBb4VEYvNvCPApyJyWCk1GIhVSjUVkTP5FaCUGguMBQgNzb87obSYDfi2b08lPyfdKHvhOPw0EQ6tgVqtYeR/oX4n58SiaZpLutyF0pYtW6hduzbh4eFkZWVx8eJFhgwZwmuvvYaHhwdvvPEG69evx9vbm+TkZMLCwnj33Xdp3NhFLh2IiMs9gC6AACabeW8APxfjveswGlIUts5W4IXixNKhQwe5Vs7p07K3WXNJ+vSz65Y5XG6OyIb3RV6vJfJGHZFNs0VyzeUfh6ZpbiMsLExeeeWVK9MHDx4Uf39/GT9+vJw5c0bq1KkjCQkJIiKSl5cnI0aMkPz2fSUBxImdcoJTWgEopVYppdIKeKwFEq2r2o7cFWgzv6ByQ4CqIvJHESHEA2GljT9twwYA/Mr7+tTxzTC/q9FHX+Pe8NSvcNuTUMlVD4w1TXNFTZo0oVevXixdupTq1auzcuVKatUyTmR5eHjQtWtXDh8+7OQo/+KUPZyI9C9subUxRTLQjL+SU0vgpyKKfgD48pqyagIDRCTCZnYtoNTXqNJiYvCsWwdTeR0WZyRD5MuwcyEE1ocHlkCzQqtQ0zStUGazGZPJhMlkumqgxFOnTrFgwQImTJjgxOiu5pL/iouIRSn1CfAIsEEp1QSjgcQoAKVUC2A20E9E8mzeej8w+JrifIHnlFJLRSRdKdUJ6AQ8WqrYcnLI2LSZgMGDHd8JrQjs+hLWToGsFOjyNHR/AUxVHLtdTdOKlPDmm2Tv21/0inbk3aI5tV98sczlxMbGEhUVxeTJk6/MO3XqFEOGDGHPnj0888wzRQ5bX55cMlFZvQrMV0ptwYjzARFJsC4LxLgR2AvIA1BKtQTOy/UNJBKAH4BIpVQuYALuEZFDpQkqY8cOLBkZ+HXvVpq3F9+5A8Y9Ucc3Qr1bYPAHUOsmx25T07QK6/J4VNnZ2ZhMJmbMmMG4ceOuLA8JCWH79u2cPXuW4cOHM2LECL799lsnRvwXPR5VEa4dj+rs2+9wYeFCmm7ZjEcVBxzZ5GTAhndh44fGkVPfqdDu/3QHspqmlVp4eDijR48u9lHS9u3b6dixI+vWraNXr16l2qY9x6PSe78SSt8Qg2+njo5JUoeiYG5n2PCeccPuU3HQ4SGdpDRNc5i8vDzy8vKumteypXEn0J49e5wR0nVc+dSfyzGfPk32ocMEDrvXvgVfOgOrJ8PeHyCoCTy0HBo4+NSipmkasHDhQpKSknj++eevzDt9+jRgnA50Bfpf9RJIi7E2S7fX9SlLHmz9CGZ3ggOroOcUeGKjTlKappWrTz/9lHPnzgFGa8CXXnqJ0NBQ+vTp4+TIDPqIqgTSYmLwCgnB1KBB2Qs7tcNoLHFmFzTqBQPehaBGZS9X0zTN6nLPFAkJCVcaUyxbtozAwMAr6/Tu3ZsdO3bQr18//P39SU9Pp1GjRkRFRREQEODE6P+iE1UxWXJySN+yhap3Dy1bs/SsFPj5P7DtE6gSDMM/h5uG6Q5kNU2zuwYNGhAdHV3oOvXr1+fDDz8sn4BKSSeqYsqMi0OSTL6OAAAKbUlEQVQyMqjStZS9UYgY16BWTYa0s3DLY9BrCvgEFv1eTdO0G5hOVMWUFrMBZTJR5dZbS/7m5GPGcPCHo6B2G3jgSwjpYP8gNU3TKiCdqIopLSYG306d8PD1Lf6bcnNg04cQ8w54eMFd06HTGN03n6ZpWgnoPWYx5Jw8Sc7Ro1S7/77iv+nPjUZjiaQD0HIo3DUNAuo6LkhN07QKSieqYkiLiQEo3vWp9PMQ+RLsWgRVQ+FvX0PTOx0coaZpWsWlE1UxpMdswCs0FFN4eMErWSxGcop8yRga/o5nodtEMJXgVKGmaZp2HZ2oiiJiNEu/996Cm6Un7jNO88VvhtDbYNAMqNmifOPUNE2roHSiKoIlPR1B5d8bRU4GxLwNm2aBtz8MmQ1tR+m++TRN0+xIJ6oi5KWmomoE43vLLVcvOLjGaHJ+Md5ITn1fhypBzglS0zStAtOJqgiWtDR8Bw7Ew8fHmJFyCla/APuWQ41mMHolhN/h3CA1TdMqMJc/R6WU6qSUOqKUGl2MdUcppbYrpeKUUu8pm4tKSqnqSqkflVKxSqktSqn2xdm+ZOfg17Ub5OXClnkw5xY4FAm9X4bHY3WS0jRNczCXPqJSSt0DjAAuFmPdVsB7QCsgGVgHjAfmWFeZC+wUkVeUUj2BZUqpxiKSXVTZfk2rwSc9IeF3aNzH6EC2uh06ptU0TdOK5OpHVNtE5G9AajHWfRT4SUSSRMQCfA48DsbRFEbC+wxARNYDOcCgogpVnh6YVo2C9HMwYgGM+kYnKU3TtHLk0olKRE6WYPVOwH6b6b3ATUqpykB7IFtE4q9ZXuQwyR6eZrh1HDz5K9x0t+7lXNM0rZy59Km/EqrF1acIUwAF1Mhn2eXlNfMrSCk1FhhrncxWA97eDW/bN1r7qwEkOTuIIrhDjKDjtDcdp325S5zN7FVQRUpUBSnsECjfZSLyMfAxgFIqTkSKPPJyNneI0x1iBB2nvek47cud4rRXWU479aeUWqWUSivgsbYURSYCVW2mAwEBzlmXXTvwU6B1vqZpmubCnHZEJSL97VzkNq4+1GwJ7BGRTKXUDsBHKVVfRE7YLI+wcwyapmmanbl0Y4rCKKVqKKU2WFv0AXwKDFBKBSmlPIDRwHwAETkPLAUesb63O2ACVhZjUx/bO3YHcYc43SFG0HHam47Tvm64OJWI2Kssu1NKdcC4N6otkADsFZFh1mX1gR1Au8utA5VSo4BnAQsQAzwv1g9oTWgRQHWgEvCkiOwo1w+kaZqmlZhLJypN0zRNc9tTf5qmadqN4UZonl4kpVQnYAnwuohEFLHu5dOLAvxCwacXPYHx9jq9qJTywbjm1txa9osikm/rSKXUasDHZpYX0FpEAgpYniYiRfbSYecYRwPPcPX9bVNEJNa63FXqchDwD4w6rApEiMiHNsvtWpcljO054G/WySUi8o7NsnDgC2sZecDDInKstHGVJsbyrrsyxDkaJ30XSxinU37X18RQrH2l3feTInJDP4B7gC+B7cDoItZthXGtrAbG0eh6jGtdl5cvAV6zvu4JnAC87RTnNGCB9XVTjP4MaxWw7tx8PuMXNtMRDqrLksQ4GuhRSFmuUpeHMa6DAtQGzgNDHVWXxY0NuAujJxYf62M/MNBm+WaM5ATwMPCrE2Is17orQ5xO+y6WME6n/K6v2V6R+0ocsJ902IdylwdQz/ocXVjlW9eZAXxuM/1/wB/W19Ux/nMNtVl+BLjXDjF6YNyJ3s1m3jrguWK+/xugt810hAPqsUQxFrZzcKW6BP6ZT13OcURdliQ24HvgZZvpl4Dl1tc3A1mAyTptAjKBDuUcY7nVXRnjdMp3sTTfx3zq06G/62u2V6x9pSP2kzf8NSpxgf4Ei6EhEJTPtossWykViNFqcr3N7CpKqUXW5v0/KKVudlKM45RSvyilopVS423mu0xdisiMa2b5YNxEfpk967IkseX3Xexos+yoiORYP0MOxtFNudZfOdddqeO0csZ3sTRxAuX6u76iBPtKu+8n9TWqkrFbf4Kl2C75bLtlMd57L/CtGD3KX3YE+FREDiulBgOxSqmmInKmHGM8C6zFOFcdDEQrpSqJyCxctC6VUgEYP8InbWbbsy5LElt+38WaBSy7dnlZlKr+yqHuyhKns76LJY3TVnn9rkvD7vvJG/6Iyk5K3J9gOWz3slHA/2xniMhkETlsfb0c4z+av9s/PKDg/hRXicgXYkgEZmGMH1bisuykOGVPB6aKyPHLM8qpLu31uZ1Zf86qu2tdF6cLfheLU76zf9elUer9ZIVOVO7Sn2Ax4rxcxrXbLrRspVQIUFVE/igihHggzBkxFhCDK9blWMAsInMKW49i1GUhShJbft/FcwUsK6wcR8YIlFvdXass30e7fBeLqTT1abfftYPYfT9ZoROViPQXEb8CHv1KUWSB/Qli9JLhY+0xw3b5NjvEeQSjJdC12y6q7AcwWulcoZSqaW2Oa6sWcLo8Y1RKTSokBpeqS2WMNN0TowkzSqkm1udS1WUhShJbft/FbTbLGiilTNY4TUCjAspxZIzlWXeljtNR30V7x2nDbr9rB7H/ftKerULc+UE+LVkwzqluAKpbp1sBZzAufnpgtM65ttnlK9bX3TGaXfrYKb5pWJuiAk0wmvrWtk63sMZS6Zr3xAF1rpkXDvwBVLFOd8JoEdakPGPEuAjc3fraF9h6ue5cqS6BbtZYgwA/6+NyU2K712VBseUT110Yp3YuN0/fx/XN0x+yvn4IY7Rse/1WihtjudZdGeJ02nexJHHarF+uv+t84o3GZl9JOewn7foB3PEBdLBW/EWMlirf2Syrj3G4Ws9m3iiM+wi2YfRDqGyWVQd+BGKtO4r2dozTB+Ni7xbrF7WfzbLOwCnbPzbGfylrCijndWATRn+IW4C7yjtGYJD1C7zeWpfTsDandqW6tL6Wax7RjqrLgmIr4G/8nHWd7cDEa8oJB37G2IFEAw0cXX/OrrsyxOm072Ip/ubl/ru22Ua++0rKYT+p+/rTNE3TXFqFvkalaZqmuT+dqDRN0zSXphOVpmma5tJ0otI0TdNcmk5UmqZpmkvTiUrTNE1zaTpRaZqmaS5NJypN0zTNpelEpWmaprk0nag0TdM0l6YTlaa5CevorZeUUhalVJR13hyl1AWl1DGl1Bhnx6hpjqD7+tM0N6KUGgF8DTwmIp8qpcKA5cDtIpLm3Og0zTF0otI0N6OU+g7oA7QFPgfeFJHSDASqaW5BJypNczNKqdoYY1HlAStE5GEnh6RpDqWvUWmamxGRBOA1jAHr1js5HE1zOH1EpWluRinlgTGAXWUgFLhJRJKcGpSmOZA+otI09/M0xgiud2OM7DrTueFommPpIypNcyNKqUbANxit/DKVUuOA+cBgEVnh3Og0zTH0EZWmuQml1H+AWKA28Ih19njr8yKl1DdOCUzTHEwfUWmapmkuTR9RaZqmaS5NJypN0zTNpelEpWmaprk0nag0TdM0l6YTlaZpmubSdKLSNE3TXJpOVJqmaZpL04lK0zRNc2k6UWmapmku7f8BtpqS7eSDByUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.special import legendre\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "#plt.rc('text', usetex=True)  #allow the use of Latex for math expressions and equations\n",
    "plt.rc('font', family='serif') #specify the default font family to be \"serif\"\n",
    "\n",
    "min = -1.0\n",
    "max = 1.0\n",
    "step = 0.05\n",
    "\n",
    "#generate the Legendre polynomials of degree up to 3, \n",
    "#i.e., P0, P1, P2, and P3.\n",
    "for n in range(4):\n",
    "    Pn = legendre(n)\n",
    "    x = np.arange(min,max+step,step)\n",
    "    y = Pn(x)\n",
    "    plt.plot(x, y, label='P%d'%n)\n",
    "\n",
    "#plot the roots of each Legendra polynomial above (except P0\n",
    "#which has no root)\n",
    "plt.scatter([], [], marker = 'x', s = 50)\n",
    "\n",
    "#for P1\n",
    "x = np.polynomial.legendre.leggauss(1)[0]\n",
    "plt.scatter(x, [0], marker = '^', s = 50) \n",
    "\n",
    "#for P2\n",
    "x = np.polynomial.legendre.leggauss(2)[0]\n",
    "plt.scatter(x, [0,0], marker = 'o', s = 50) \n",
    "\n",
    "#for P3\n",
    "x = np.polynomial.legendre.leggauss(3)[0]\n",
    "plt.scatter(x, [0,0,0], marker = 'd', s = 50) \n",
    "    \n",
    "#plot the x-axis    \n",
    "plt.plot([-1,1],[0,0],'k--')\n",
    "    \n",
    "plt.xlim(-1.0,1.0)\n",
    "plt.ylim(-1.0,1.1)\n",
    "\n",
    "#set x and y axes label and font size\n",
    "plt.xlabel('x',fontsize=16)\n",
    "plt.ylabel('P(x)',fontsize=16)\n",
    "\n",
    "#show the ticks on both axes and set the font size \n",
    "plt.tick_params(axis='both', which='major', labelsize=12)\n",
    "\n",
    "plt.legend(loc='lower right',fontsize=14,frameon=False)\n",
    "plt.savefig(\"legendre_polynomials.png\",dpi=300, bbox_inches='tight', pad_inches=0)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also get the Gauss quadrature nodes $x_i$ with the functions in Numpy. The following example shows how to get the nodes $x_i$ and corresponding weights $a_i$ for the 2-point Gauss quadrature rule. Note the nodes $x_i$ are defined on the interval $[-1,1]$. A linear transformation is needed if you want to use these nodes for integration on a general interval $[a,b]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2-point Gauss quadrature nodes x_i:  [-0.57735027  0.57735027]\n",
      "2-point Gauss quadrature weights a_i:  [1. 1.]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "#get the 2-point Gauss quadrature nodes $x_0$ and $x_1$\n",
    "print(\"2-point Gauss quadrature nodes x_i: \",np.polynomial.legendre.leggauss(2)[0])\n",
    "\n",
    "#and the corresponding weights\n",
    "print(\"2-point Gauss quadrature weights a_i: \", np.polynomial.legendre.leggauss(2)[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}