OptimalSwitching_SHE_Multilevel.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/bingsen-wang/PowerElectronics/blob/main/OptimalSwitching_SHE_Multilevel.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jAMKoELNpsK-"
      },
      "source": [
        "<center><font size='5'><b>Optimal Switching: Selective Harmonic Elimination for Multilevel</font></b>\n",
        "\n",
        "Dr. Bingsen Wang\n",
        "\n",
        "12/3/2023\n",
        "</center>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "#Background"
      ],
      "metadata": {
        "id": "APo2BMHCTsaX"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Without loss of generality, a square waveform with firing angle $\\theta$ only has the Fourier coefficients for the sine terms.\n",
        "$$b_n = \\frac{4}{n\\pi}\\cos(n\\theta)\\tag{1}$$\n",
        "\n",
        "If the square waveform has an additional firing angle to create a \"staricase\", then the Fourier coefficient for the fundamental component is given by\n",
        "$$b_1 = \\frac{4}{\\pi}(\\cos\\theta_1 + \\cos\\theta_2)\\tag{2}$$\n",
        "Due to the extra freedom with the firing angle $\\theta_2$, an additional constraint can be satisfied, for instance, to eliminate the $n$th-order harmonic, i.e.\n",
        "$$b_n =0 = \\frac{4}{n\\pi}[\\cos(n\\theta_1) + \\cos(n\\theta_2)]\\tag{3}$$"
      ],
      "metadata": {
        "id": "nmdKU9MtTv1y"
      }
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Hv0YxDrWU4h8"
      },
      "source": [
        "#Python Code\n",
        "The Python code illustrates the waveform spectrum by varying firing angle $\\theta_1$ with keeping the constraint in (3) with n=3, 5, 7."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 49,
      "metadata": {
        "id": "XeGsHYTt8IXU",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 928
        },
        "outputId": "f3fabd60-b6c1-4940-991b-c4baf52a5a39"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 900x900 with 3 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA5cAAAOPCAYAAABfEDUTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACXfElEQVR4nOzdd3iTVePG8bu7tEALdACyEUT2FNkgAiK4UERFEfV1IOB6xYXIEBfgiz+GOBBQHCiIgIoM2SB7tgwB2UNoC6VQWrqe3x9pQkqTNm3aJinfz3X1IuQZ55zk6bhzznOOl2EYhgAAAAAAcIK3qysAAAAAAPB8hEsAAAAAgNMIlwAAAAAApxEuAQAAAABOI1wCAAAAAJxGuAQAAAAAOI1wCQAAAABwGuESAAAAAOA0wiUAAAAAwGmESwAAAACA0wiXAAAAAACnXffhctYsKSnJ1bUAAAAAAM/mUeEyI0NKTS2YcxmG9MIL0i+/SIGB+T9PUpLUuLEUFVUw9QIAAAAAT+TW4fKzz6SOHaV69aSICMnPTxoypGDOPWiQKRDOnCl5eeX/PCVKSF9+KfXsKW3aVDB1AwAUvpEjR6pmzZoqW7as0tPTXV0dAAA8nq+rK5CTsDCpeXNp+XIpJsb0XJMmzp/3009Nw2GjoiR/f+fP16KF9OGHUpcu0vbtUo0azp8TAFC4hg8froyMDC1fvlw+Pj6urg4AAB7PrXsuH3hAGjdOGjny6nNt2zp3zr//ll5+WRo/XqpY0bFjzp837b9nj/19Hn5YuvNOqW9fKS3NuToCAIrGzp071a5dO1dXAwCAYsGtw6WZebhphQpSzZrOnevZZ6XataVHH83bMa+8In30Uc77jRkj7dghvfuuU1UEABSR9evXq62zn1oCAABJHhIuV60y/evs7/+FC03nGjpU8naw5RkZ0pIlpse5DcmtXFl67jlTb+uZM87VFQBQuA4cOKDY2Fi1adPG1VUBAKBYcPtwmZx8tefS2ZFLH3wglS0r3Xef48fs3ClduGB63LFj7vsPHGiaQXbMmHxVEShS/fv3l5eXV45fQUFBql27tp599lntyWlsODzCqlWrFBYWpl69eskwjDwff+TIkWzXyMqVKwu+onk0atQo9enTR0OGDNGYMWP03//+V7Vq1dLq1ast+5w7d06DBg3SsGHD9NJLL+mnn35SgwYNFBIS4nT5ixcv1lNPPaWbbrpJISEhCgwMVOXKldWtWzeNHz9e8fHxTpfhqJUrV2Z7j44cOVJk5efFSy+9lOvPIHtf1m1ytM3OXv/OcGXZ7syTrlcAuXP7cLlhg3TliumxMz2X0dHS2rWm+zgDAhw/ztxrWqaM1LBh7vvfeKPUvr1pptsi/FsCyJf33ntPUVFRGj16tOW5xYsXKyoqSlFRUdqyZYumTZumG264QV988YWaNGmiadOmFVj51n9UFMR+yN2cOXMUFxenX375RefOncu2PbfX+oYbbrBcH4Vt8+bNuvfee9WmTRs1a9ZM7dq10zfffKPUa9akWr58uVq3bq0nnnhC48aNU9myZTVixAidP3/e0sYLFy6oU6dOatKkid59912NGzdOX3zxhdP3Wx4+fFitW7fWHXfcoWPHjum1117T77//rrVr12ry5MmqXbu23n77bVWvXl1ffvmlU2U5qkWLFoqKitLixYuLpDxnvPHGG4qKisryc2XatGmWa8zWl/XPKzNH25zb9V+YXFl2blz5M9aTrlcAuXPr2WKlq+GudGmpUaP8n+eHH0z/3nFH3o4zf+jdvr3jQ2l79jTV+9dfpccey1t5QFG64YYbdMMNN2jLli2W52rXrq1q1apZ/t+sWTP16dNHDz30kH766Sc9/fTTatSokZo1a+aCGsNZzz33nLZv36527dqpXLlyeT7ez89P9evXL4SaZfXLL79o6NCh+uGHH9Qo84f/zJkz1a9fP33xxRdaunSpSpQoIUm6cuWKunbtqo8//lgtWrTQf/7zH0lSbGys5XxvvPGGSpQooaeeekqS5OvrK39/f6fCZVRUlDp27KiEhAR9//33evjhh7Ptc/fdd2vIkCHq1q2bnnnmGR04cEBjCnloS3BwsOrXr6+SJUsWajkFoXz58ipfvnyW96p69eo5XmPWP6/MHG2zs9e/M1xZtjvzpOsVQO7cvufSHC5btXI83Nkyf77pX0eGtpoZhrRmTd6PMwfYn392/BjAnXl5een999+XJGVkZGjixIkurhHyq169elq7dq0++OADV1fFrsTERA0ePFi//vqrJVhK0mOPPaYGDRpo3bp1+t///md5vnv37vLx8dGKFSvUpUuXbOc7f/68pk+froceesjy3KlTp3Tw4MF8T+YTGxurHj166Ny5c/rwww9tBkuzKlWqaNGiRQoKCtLYsWM1Y8aMfJUJk/DwcLVs2VIBeRmGlMmV178nfO8BgLPcOlympJiGxUqm+y3N9zI2a2YafhoeLnXtKm3cmPN54uJMy4hUrWoa3mrPmTNSrVpStWqmr4gIyfxh6iefXH2+fv2clxupW1cKCpIWLza1ASgOatasqeDgYEnS9u3bXVwbFGfr1q3TyZMn1adPH12+fDnLtptvvtmyj7W0tDStXbtWHTp0yHa+9evX68qVK7rtttssz61atUo1atRQRUfXpLrGyJEjdfz4cVWqVEkvvvhirvtXrVpVgwYNkiS9+uqrunTpUr7KhdSjRw9t2LBBFSpUcHVVAADXcOtwuXmzKVBKUmSk1LKlKawtXiwdPCjt3Wtat7JTp5zXoNywwdQL2aBBzuVFRkoHDkhHjpi+Mv8OUI0aV587csR0/6ZvDgOKvb1NQ3iTk3OuF+Bp/Pz8JJn+kDc7cuSIPvzwQ3Xp0kUVKlSQn5+fSpcurWbNmmnUqFFKSEjIdh7z/T2dOnWyPGdrkhhH97NmGIZmz56t7t27Kzw8XP7+/oqIiFC3bt30zTffKD09Pcv+9iao2b17t/r27auKFSvK19c3y0QTtvbfuXOnevXqpfDwcAUGBqpu3br6+OOP8zxxh71Jlq5tZ/v27eXl5ZVtkiVbxx85ciTXiXjy81pbK6j2S6aeS0naunVrtvZdybwJv3Tp0lme37Jli5KTk9W6dets50tOTpaXl5fq1q2bpb3mIbEff/xxnup34cIFffHFF5Kk3r17yzenXwhWHnnkEUlSXFycvv76a8vzeb0GrV26dEnDhg3TTTfdpMDAQIWHh+vOO+/UihUrHKpTXr5fnKlnQejYsaM65mUYkZXcrn9727dt26a77rpLZcqUUWhoqDp16qTly5dbjtuyZYvuvPNOlSlTRiVLltRtt92m9evXF0jZef2eKsyfxdby+jPWmrPXKwAPYLix994zDFMsNIxKlQxj5crs+7zwgmn788/bP8+4caZ9BgzIW/mdOpmOe/zxvB1nGIbxyCOmY7/+Ou/HAkVt+vTphiRDknH48GGb+/z777+WfXr16mV5vkOHDoYko0OHDsbPP/9sbN682ViwYIHRt29fw8vLy6hZs6Zx+vTpLOe6dOmSERUVZUybNs1yzqioqCxfly5dcng/s+TkZOP+++83JBmtW7c2fvzxR2PDhg3Gd999ZzRr1syQZHTu3NlITEy0HJOSkmI5l7mM0aNHGzfeeKPxxRdfGJs3bzZ++uknIyIiwpBkc/93333XaNasmfHTTz8ZmzZtMr766isjLCzMkGS88cYbeXovTpw4YWzYsMHw8fExJBmffvpptnaeOnXK8Pb2NiQZI0eOzHZ8VFSUUaVKFePJJ580oqKijJSUFJv1XrFiRZ7fE2v22j99+nQjMjIyX+03DMNISkoy7rvvPuPRRx810tLSsmyrWrWqIcmYPn16luc/+OAD45ZbbrF5vpMnTxqBgYHGuXPnDMMwjOXLlxvh4eHGuHHjjH379hkzZszIU/1+/vlnS9tnzZrl8HHp6elGiRIlDEnGXXfdZXk+r9eg2ZkzZ4x69epZvv8WLFhgbNmyxZg5c6ZRp04d4913383x+zqv3y/5racjVqxYYfO6tNahQwejQ4cOOZ7n8OHDNtuc2/Vva/uoUaOM1q1bG/PmzTM2bNhgjBkzxggICDB8fX2NRYsWGatXrzbuuusuY9GiRcbatWuNIUOGGF5eXkZgYKARFRXlVNn5+ZlSmD+LzfLzM9bM2esVgGdw63DZtevVcDl/vu19Xn3VtL1bN/vnGTDAtM+oUY6XfeWKYZQoYTpu6tS81dswDOO110zHvvZa3o8Fipoj4fLVV1+17PPHH39Ynu/QoYPRrFkzIyUlJdsxb775ZrYwas36D8qcOLrfc889Z0gy2rVrly2UpKamGo0bNzYkGc8++6zN481lhIaGZnsd3nvvvWzlm/cvX768ERMTk2XbwoULDUlGiRIljISEhBzrbUunTp0MScbzNj45mzhxoqXsBg0aZNu+d+9eQ5KxZs0am+fO6Y94R19r6/PYav+SJUucar8ta9asMSQZTZs2NVJTU7NsGzBggDFu3Di7x3777bdGv379jBEjRhhTp041Fi9ebPTo0cN44YUXjMuXL+epHsOHD7e0ffPmzXk6tm7duoYko0qVKja35+UavPPOOw1JRsuWLbN9/8XExBgVK1bM8fvame+XvH6v5Kaww6Wtutsrx7z9hhtuMOLj47NsGzlypOUa7NWrV7bX/bHHHjMkGY899phTZefnZ0pR/Cx25ppx9noF4BncNlymphpGyZKmgPboo/b3e/hh0z4PPWR/nwceMO3zf//nePlr114Ntvv2OX6c2dixpmOffDLvxwJFzV64TElJMfbu3Wu8+OKLlp6y1675xGT69OnGSlvDCgzDOHLkiCHJ8PHxMS5cuJBte0GGy7179xpeXl6GJGP16tU29/nuu+8MSYafn5/x77//ZttuLmPgwIHZtl28eDHbHzzm/V999dVs+ycnJ1teM3t/SOZk8uTJhiSjQoUKRnp6epZtHTp0MNq0aWMpf//+/Vm2jxo1yqhYsaKRkZFh89wFHS5ttT8lJcXS+5qf9l/r2LFjxs0332y0adPGOHPmjNPnc8bzzz9vafvff/+dp2NvvfVWQ5IRFBRkc7uj1+DGjRst+y5YsMDmucwhz9Yf685+v+T1eyU31tedt7e34ePjk+3L3OOVk4IMl6+//nq2bevWrbNsnzJlSrbtM2fONCQZVatWdars/PxMKeyfxc5cM85erwA8h9vec7l1q2Se7+CZZ+zvZ761IaeZ8c3nCQx0vHzzLLUREdJNNzl+nFnmDPm6cCHvxwKudOONN8rX11e+vr4KCAjQzTffrEmTJqlVq1aaP3++Pvrooyz79+/f3+YkKpJpEhNJSk9P14EDBwq13rNnz5ZhGAoMDNStt95qc586depIklJTU7XavM6QDbZmEC1ZsmSWJVqstWjRIttzAQEBCgsLkyT9+++/uVU/m/vuu09eXl46ffq0/vrrL8vzZ86c0Zo1a/Tee+9Z6jN79uwsx86ZM0e9evUqsjXrbLXfz8/PstxCftpv9vDDD6t9+/Zq2LChwsPD9e233yoiIiLf5/MUuV2DCxYssDxvfb+ctZyWWSmo75e8fq84YurUqdqxY0e2r+bNm+f7nPlha7kl62vP1vby5ctLkk6fPu1U2fn5mVLYP4uduWacvV4BeA63XefSHO5KlDBN5GPLwYOmCXYk6fbb7Z/LvN62g3MuSLq6vmU+Z6knXMJjLVy4MMsMmiVKlFCFChUUFBRkc/+MjAx9//33+uGHH7Rz507FxcVlW+ReUqHPjrlz505JpslbzOsf5uTYsWN2t4WHh+epbHtr1pnrkZycnKfzSVKFChXUunVrrVu3Tj///LPlj/i5c+cqPDxc7dq10/3336+PP/5Yc+bM0VtvvSVJOnDggHbt2qUJEybkucz8Koz2m/2QuUhxRkaGPv30U910000aMmSIRo0aJW9n1qfKQceOHXXixIlsz48ZM8YywYqZrUlScmLeP7eAnNs1aJ7oKCwszO76gOagY0tBfb/k9XvFEfbWuTTPVl1UbF3X1pM35bQ9xcmp4vPzPVXYP4uduWacvV4BeA63DZfmcNeqleTvb3sf84f1FStKNj7kszD/Tezo3zfp6ZJ5lvv8fpCWn0ALuIPatWs73OuQmpqqHj16aOnSpapYsaJee+01NWnSRGXLlrXs0yBzmmYjH7OG5sWFzE9yIiMj9eeff+a6f2RkpN1tPj4+eSo7r/s76v7779e6des0d+5cjR8/XpKp9+C+++6Tt7e3HnjgAX388cfavn27Dh06pBo1amj27NmKjIws0l6Awmq/NW9vbw0aNEi//fab3nvvPV25ckVjx44tlLKOHDmio0ePZnveHAwbNmxoee6ff/5xuEctIyNDhw8fliQ1btw4x31ze03Ndcnpj3zz7M62FNT3S1G892Y5zVZcGHL78KKwPtyQ8v66FsXPYmeuGWevVwCewy2jT0aGtHat6XFOs47/+KPp3wceMC3/YU+pUqZ/HQ2X27ZdHUqb37/PzEuo2PmADigWPv30Uy1dulS+vr5asmSJ6tWr57K6hISESDJ9qm6r18MT3X///XrllVd07Ngxbd68WdWqVdPq1astvZQtW7ZUpUqVdOLECc2ePVuvv/56lvBZHN1zzz1avHixxo8fr8GDB6tKlSoFXsYR85AYOzp37ix/f3+lpKRo48aN6tOnj0PnjY6OVlLmL4c777zTqTqal2K5dh1Qa7Z6rcyK4/fL9awofhY7c804e70C8BxF8tfH7Nmz9cknnzi8/44dV4eT2hmar127pMwRGnr00avPL10qnT2bdd/KlU3/xsQ4Vr55SG7JktK1Hy6vXStNn577OS5eNP2beXsEUCwtW7ZMkqm305XBUpIaNWokyfTpek73+G3atElTp051+p6oolClShXLvVc///yz5s6dq9DQUMtaf15eXurVq5ck032W//zzj3bs2KHevXu7qsoFYtKkSWrZsqV++eWXbNvMverp6enavHlzEdfMJCQkRM9kTgYwe/bsHNf1s/b9999LkipWrKhHrX9x5YN5zc64uDhdNP/CuUZO3wee/P1y5coVnThxwrLmKYrmZ7Ez14yz1ysAz1Ho4TI1NVVPPfWUXn31VcXFxTl0jDncBQVJt9xiex9zwKtb9+qQ2IwMqX9/6do5LGrWNP17/LhjdbYeknvtyJS335Yc+aDcfLtOjRqOlQl4ooyMDEn2h1nl1gNkff+S9TmWLVumhQsX5mm/3r17W3rrfvvtN7tlDhgwQC+88EKR37+VX/fff78kU7icM2eO7r333iyvxwMPPCDJtJj7uHHjFB4ebndSD0c4+p4UpjfeeEObNm3SqFGjsm2z7vkIDQ0tkvrYMnz4cFWuXFknTpzQxIkTc93/+PHjmjRpkry8vPTJJ584dM9aTu655x7LY3sL0K9Zs8bu8Z78/bJ+/XpVrlxZ680z+qFIfhY7c804e70C8ByFHi6//fZbXbx4UYZhWO4Zyo05XOZ0v+XPP5v+vbbXsl496dr5BTI/bFPmrS65Mk+kdu1kaEuWmIa7du6c+znM4bJ2bcfKBDyR+b6+v//+22Yv0uTJk3M83noCh3PnzlkeP/vssxozZkye9qtTp46ee+45SdJ7771n88OsadOmadu2bRo8eLBlmJa7M4fLgwcPatmyZZb/m7Vp08by+nz++ee69957nboPztH3pDBVrlxZgYGB6tevX7ZtUVFRkkwTybRq1apI6mNLWFiYfvvtN5UpU0avvfaafvrpJ7v7Hj9+XN27d1diYqLef//9AulZbtGihWVo7fvvv59tSGFsbGyOobe4fr9cr4riZ7Ez14yz1ysAz1Go91ympqZq9OjRqlq1qtLS0jRhwgS9/PLLdmdBk0wrS5o/vMrpfkvzsFnzLLEZGdKwYdKHH2bft2lTKTjYNIw2IyPn+zMlKTJS2rfv6nBaSTp/XnrhhauhNjeZf/9kC6iAOzl58qTOnz+vkydPWp7bv3+/ZTbB3O6rGTRokGVmwjvvvFNvvfWWbr31ViUkJGjWrFlZhjUePnxYYWFhql69uuUT7RtvvFF169bVnj17NHr0aD3yyCNasmSJ/vnnHw0ePNhyrKP7jR8/XnFxcfrxxx/VsmVLvfXWW2rUqJFiY2M1f/58ffHFF+rWrVu2HrHo6Ogs/zfX1d/fX7VtfEJkb39z2/bv36+UlBTLH1AnT55UdHS0KlWqlOfethtvvFENGzbUrl27FBISotuvmRrb29tb9913n6ZMmSLDMCw9mddKTU3V33//bbPe1vVy5LUu7Pa/+eab+vbbb/Xss89meT4lJUUzZ86Ul5eXJk2aZHcG46LSsGFDbd26VX379lWfPn301Vdf6aGHHlLt2rUVGBioU6dOaenSpZo2bZoCAgL0ww8/6KGHHrJ5rrxeg5I0ffp0derUSRs3blTXrl313//+VxUrVtTevXv17rvvqnPnzvruu+8kXf2+rlatmmW2zvx8v+SnnjnZv3+/jh07ZpmJVJK2bt2qtLQ0u8dY72uWmJiow4cP69SpU1nOfenSJVWvXl3+/v65Xv85Xdfm422dv379+pbyD1t9im0+30033eRU2Y58TxXVz+L8/oyVnL9eAXiIwlxEc9q0aYYk4+677zZuueUWIygoyBg6dGiOx0RFGYYpYhrG2rX293v0UdM+a9YYRlKSYTz/vGEMGGB//+7dTftHR+de76+/Nu37yium/+/daxi33GIY33yT+7GGYRgxMabja9Z0bH/AVR5//HHLotW2vhxx8eJFY9iwYcbNN99sBAQEGAEBAUatWrWMgQMHWhbutv66dvHvffv2GXfeeacREhJiBAQEGLVr1zbef/99Iy0tLV/7GYZhzJ8/3+jZs6cRERFh+Pr6GqGhoUaHDh2Mr776ykhPT8+2v73257YQur22Va1a1eb26dOnO/SaXmvkyJGGJKNfv342ty9btsyQZJQtW9ZITU21uY/14vK51Su317oo2j9t2jSjSZMmxsiRI40FCxYY8+bNM7p27WpUrFjRmDdvnsPnKSqLFi0ynnzySaN27dpGqVKlDH9/f6NixYpGt27djPHjxxsJCQk5Hp/Xa9AsISHBGDp0qHHjjTca/v7+RmhoqNG+fXvj+++/t/me//HHH9nOkZfvl/zW054XX3wxx59BOX1Z/yxZsWJFjvs5cv3n9/jcyne2bMNw7HuqqH4W5/WasVYQ1ysA9+ZlGIWzPkBqaqrq1Kmjxo0bq3Tp0jpw4IDatGmjKVOm6PDhw3Z7L48dM91D2aSJtHCh/V7GCxekl16SFi2SSpeWHntMeust+/vPmCE98YQ0dqz06qu513/yZOmTT0xLilSqJI0c6dhwWMnUu/nAA6ZyCmmmfAAo9i5fvqyVK1dqz5498vb2Vt26dXX77bdnuT8MAAC4j0ILl4sWLVL37t21Y8cOffLJJzpw4IB++eUXVatWTWPGjNHAgQMLo1i7LlwwrYfZrNnVCXsKyxNPmMLs5s2Sg8ufAQAAAIBHK7QJfdq0aaM1a9ZYpq6WTBMwrFmzxu49QYUpJET6z39MS4ns31945SQnSwsWmEIlwRIAAADA9aLQwmWpUqXUtm3bbM83bdpUkZGRhVVsjl59VQoIMA15LSw//iidOye99lrhlQEAAAAA7qbQlyJxJ5Urm+7L/Pxz6ejRgj9/errpHssWLaRrVgsAAAAAgGLtugqXkvT666a1J//734I/94wZpiVMPv889+VOAAAAAKA4ue4ikL+/NH++tHKlKQwWlBMnpCFDpPfeM810CwAAAADXk+suXEpS9eqmSXfeekv66y/nz5eaKj36qHTXXaaeUQAAAAC43lyX4VKSWrc2rZH53HNSYmL+z2MYUv/+Uq1a0rRpBVY9AAAAAPAo1/VK1A0bSrt2OXeOpCTpjjukxx4rmDoBAAAAgCe6bnsuC0pQEMESAAAAAK7rnku4lmFIsbGurgXcUViY5OXl6loAAAAgLwiXcImkJKllSykqytU1gTu6+WZp82YpONjVNQEAAICjGBYLl1izhmAJ+/bulVascHUtAAAAkBf0XMIlUlKuPr7pJtNsu8A//5iCpZT1GgEAAID7I1zCJQzj6uN+/UxrjgLjxklDhpgeW18jAAAAcH8Mi4VLWAcHJm6BmfW1QLgEAADwLIRLuAThErYQLgEAADwX4RIuQbiELYRLAAAAz0W4hEsQLmEL4RIAAMBzES7hEoRL2EK4BAAA8FyES7gE4RK2EC4BAAA8F+ESLkG4hC2ESwAAAM9FuIRLEC5hC+ESAADAcxEu4RKES+SGcAkAAOBZCJdwCcIlbKHnEgAAwHMRLuEShEvYQrgEAADwXIRLuAThErYQLgEAADwX4RIuQbiELYRLAAAAz0W4hEsQLmEL4RIAAMBzES7hEoRL2EK4BAAA8FyES7gE4RK2EC4BAAA8F+ESLkG4hC2ESwAAAM9FuIRLEC5hC+ESAADAcxEu4RKES9hCuAQAAPBchEu4BMEBthAuAQAAPBfhEi5BzyVsIVwCAAB4LsIlXIJwCVsIlwAAAJ6LcAmXIFzCFsIlAACA5yJcwiUIl7CFcAkAAOC5CJdwCcIlbCFcAgAAeC7CJVyCcAlbCJcAAACei3AJlyBcwhbCJQAAgOciXMIlCJewhXAJAADguQiXcAnCJWwhXAIAAHguwiVcgnAJWwiXAAAAnotwCZcgXMIWwiUAAIDnIlzCJQiXsIVwCQAA4LkIl3AJwiVsIVwCAAB4LsIlXIJwCVsIlwAAAJ6LcAmXIFzCFsIlAACA5yJcwiUIl7CFcAkAAOC5CJdwCcIlckO4BAAA8CyES7gE4RK20HMJAADguQiXcAnCJWwhXAIAAHguwiVcgnAJWwiXAAAAnotwCZcgXMIWwiUAAIDnIlzCJQiXsIVwCQAA4LkIl3AJwiVsIVwCAAB4LsIlXKZkyYvq2HGl0tIuuroqcBN80AAAAOC5CJdwCcOQSpW6pI4dVykt7ZKrqwM3Qc8lAACA5yJcwiUYFgtbCJcAAACei3AJlyBcwhbCJQAAgOfyLYpCNm+W9u2TSpcuitLgCa5ckcqVc3Ut4G6sw+UHH0jjx7uuLnAvPj7So49KEye6uiYAAMCeQg+XZ85Iu3ebHl9k3hZYSUoK1K5dDdS2baCrqwI3UabM1cepqaYvwGzSJOn116VKlVxdEwAAYEuhh8vLl68+Ll2aPwpgrYwqVuylTp1cXQ+4i+bNpZdflhYvdnVN4E5OnJASEkyPrX+nAAAA91Lo4dL6vqkePaTvvy/sEuEp0tLSlJCQoIyM0vL2LpIR2nBzXl7S//7n6lrA3Tz+uPTNN6bH3IsLAID7KvQJfZi4BfbExMRo4sSJiomJcXVVAAAAADiJcAkAcGvMIgwAgGcgXAIAAAAAnEa4BAC4NXouAQDwDIRLAAAAAIDTinS2WMIlrFWoUEHDhw93dTUAuDl6LgEA8Az0XAIAAAAAnEa4hMvExsbqq6++UmxsrKurAsCN0XMJAIBnIFzCZVJTU3XixAmlpqa6uioAAAAAnES4BAC4NXouAQDwDIRLAAAAAIDTCJcAALdGzyUAAJ6BcAmXCQ0N1X333afQ0FBXVwUAAACAk1jnEi5TokQJNWzY0NXVAODm6LkEAMAz0HMJl0lMTNSmTZuUmJjo6qoAAAAAcBLhEi6TkJCgP/74QwkJCa6uCgA3Rs8lAACegXAJAAAAAHAa4RIA4Nb43QEAgGcgXAIAPAbDYgEAcF+ES7iMv7+/atasKX9/f1dXBQAAAICTWIoELlOuXDk9+uijrq4GADfHhD4AAHgGei7hMhkZGbpy5YoyMjJcXRUAAAAATiJcwmXOnDmjDz/8UGfOnHF1VQC4MXouAQDwDIRLAAAAAIDTCJcAALdGzyUAAJ6BcAkAAAAAcBrhEgDg1ui5BADAM7AUCVwmIiJCr776qgIDA11dFQAAAABOIlzCZXx8fBQcHOzqagBwc/RcAgDgGRgWC5c5d+6cfvjhB507d87VVQEAAADgJMIlXObKlSvav3+/rly54uqqAHBj9FwCAOAZCJcAAAAAAKcRLgEAbo2eSwAAPAPhEgAAAADgNMIlXKZUqVLq2rWrSpUq5eqqAHBj9FwCAOAZWIoELlOyZEm1atXK1dUAAAAAUADouYTLJCUlaffu3UpKSnJ1VQC4MXouAQDwDIRLuEx8fLzmzJmj+Ph4V1cFAAAAgJMIlwAAt0bPJQAAnoFwCQAAAABwGuESAODW6LkEAMAzEC7hMr6+vipfvrx8fQt90mIAAAAAhYylSOAy4eHhevbZZ11dDQBujp5LAAA8Az2XAAAAAACnES7hMqdPn9bo0aN1+vRpV1cFgBuj5xIAAM9AuIRLpaenu7oKAAAAwHXphx9+UJcuXbRo0SIZBfAJLuESAODW6LkEAKBwNGnSRImJierevbtatWrldMjM84Q+ly9fVlxcnMP7x8RIUqKkK7pw4biOH89riSiuYmNjJUn//vuv0tLSXFwbAO4qIeHq47Nnxe8RAAAKSHBwsGbNmqU1a9Zo/Pjx6t69u5o0aaIhQ4booYceklceewe9jDxG09atW2v9+vV5KgSwpUKFCnr22Wf1+eefc98lAAAA4Ea+/PJL/ec//8nTMXkOl4cPH9b+/fsd3n/dOmn06I9lGCf19NP/0/3356l+KMYyMjKUnJyswMBAeXsX+ghtAB7qq6+k2bNNj8eMkRo2dG19AAAobk6fPq1Zs2bpzz//VOnSpfXwww/rk08+kZ+fX57Ok+dhsdWrV1f16tVz3OfkyZN699139ccffyglJUWVKvkpMjJSX3zRLa/FoRi59rq45ZZbNHHiRFWpUsXVVQPgxlasuPq4RQupY0eXVQUAgGLl7Nmzeuutt/T111+rXLlyGjt2rJ599lkFBQXl63wF3l20bNkyNWrUSOXLl9fff/+tw4cPKyUlRdHR0bp8+XJBFwcPYeu6iI+P12233cZ1AQAAALjA7Nmz9ccff2jMmDE6dOiQXn755XwHS6mAw+XatWvVs2dPDRw4UCNGjFBgYKACAwNVt25dJScna/LkyQVZHDyEvevixRdf1D///MN1ASBHzDQOAIBjzp8/r/Hjx2vPnj0O7T9w4ECdPHnS6VBpVmDh8syZM+rTp4+qVq2qYcOGZdlWpkwZSdL3339fUMXBQ+R0XTTMvHGK6wKAo1iKBAAA+5599lm98sor+uijj1xSfp7vubTFMAw98sgjOnXqlL799lv5+mY9rb+/vyRp165dSkpKUokSJQqiWLi53K6LkJAQSVwXAHJGzyUAALnLyMjQkiVLJJnWr3SFAum5nDFjhpYvX67w8HD17t072/bU1FRJpgb/+++/BVEkPEBu18WlS5ckcV0AcBw9lwAA2LZz505duHBBktTRRbPfOR0uk5OTLcMdn3zySUsv5bX7mMXExDhbJDyAI9eF9bXAdQHAHnouAQDI3apVqySZbkls6KJ1u5wOl19//bVOnjwpSXrmmWds7mNO0JLkxV8J1wVHrou9e/daHnNdAHAEPZcAANi2evVqSVL79u1dtoa806VOmjRJktS0aVPVqFHD5j7nz5+3PDZP7oPizZHrYteuXZbHXBcA7OGzJwAAcmYYhtasWSPJdUNiJSfD5datWxUdHS1Juvfee+3uZz3kMTw83Jki4QEcvS42bdpkecx1AcAR9FwCADzdsmXL9MQTT+T5uBkzZqhcuXJ65JFHdObMGdWqVUvVqlVTtWrVFBERodjYWEnSJ598Ynm+fv36SktLK+gm2OVUuJw9e7blsb0QkZycbGlo1apVLTOEovhy9LrYvHmzJK4LADmj5xIAUFxcuXJFgwYN0sCBA/N87CuvvKJz587phx9+0OXLl3XgwAEdOXJER44c0aBBgyRJNWrUsDx35MgRRUdHZ1uxoTA5VdLvv/9uOomvr+666y6b+yQlJSkjI0OS1KxZM2eKg4dw9Lq4cuWKJK4LAI6j5xIA4Mk+/vhj3X333WrevHmej33kkUc0efJkSdK///6r6tWrW7aZJ/Np165dwVQ0n/Ldc3nu3DnL0Mf+/ftnScjWXw8++KDlmPbt29s819GjR1WlSpX8VgVupCCuixkzZqhr165q1aqVatWqpZdeeklnzpwp0nYAcB/0XAIAioMrV67o66+/1ltvvZWv48eNG6f69etLksqVK2d5PiUlRRs2bJDkweHS+n65W2+91e5+5hQtSffcc0+WbRcuXNCXX36p5s2b6/jx4/mtCtyIs9fFsGHDtH//fi1evFjr16/XvHnzNGfOHDVo0CDLBEAArk/0XAIAPNV3332n3r17270dzDAM/fbbb0pKSrK5PTAwUM8//7z8/PyydMxt3rzZckzbtm0LvuJ5kO9w+ffff1se33LLLTb3OXHihKKioiRJwcHBqlatmmXbXXfdpXvuuUeHDh3KbxXghvJ6XTRp0sRyXWzevFnr16/X+++/b1mapF69evrss88UExOj3r17W4ZYA7h+0HMJACgOfvjhB/Xt29fu9i+//FJ33XWXli9fbncf8xqWgYGBlufMnTYRERG66aabCq7C+ZDvey7NPY2BgYGqW7euzX3mzZtneRwZGZll26+//mp5/MMPP1gm/YFny+t1Yb0G5ueff67nn38+2/533nmnypQpo/3792vZsmXq0qVLwVYagMeg5xIA4ImSkpJ06tQp3XzzzXb3+fzzzyVJN9xwg9199u/fn+2WMvP6lq7utZSc6LlMSEiQJNWqVUs+Pj429/nmm28kSf7+/oqIiMhvUfAgebkuQkND9dhjj1me37p1q/r162eZEMjM29tbtWrVkiTt3r27MKoNwI3RcwkA8HQbN25U5cqV7W6/cOGCtm/fripVqqhRo0Z29/v555+z3FKWnp6udevWSXL9/ZaSk0uRSKZlJGzZt2+fZamJOnXq2A0aKJ4cuS6ef/55BQcHW7alpaUpMTFRCxcuzHZcenq6JDEsFrjO0XMJAPBEx48fV5kyZexu37dvnwzD0GOPPWa5Pexac+fOVWxsbJYeym3btunSpUuSPDxchoaGSpJKly5tc7u5dyoyMtIyqxGKv7xcF6+99lqWbW+//bZatmypp556Ksvzqamp2rdvnySpadOmBVxjAO6OnksAgKeLiYnRyZMn7W7fs2ePJOm5556zuf306dMaMGCAnn/++Syddub7LUuWLKnGjRtnOWbt2rWaPn26kzXPm3yHy5o1a5pO4J39FBcvXrSMGR4/frz8/f3zWww8TF6ui2tnyurTp482bNiQLUD++OOPSkxMVN26de0uZwPg+kDPJQDAE6WlpWnDhg06ffp0tm27du3SK6+8Ikk6f/58tu0nT55U586d5e3trYEDB2bZZr7fslWrVtlGir799ttFvtxjvsNlq1atJJnWNbzW2LFjde7cOT388MN6+OGH8187eJyCvi4yMjI0duxY+fj46LPPPrMZWgEUb/RcAgA8XaVKlZSamqqnn35a8fHxkkzrXn722Wdq27atZWjrSy+9pAsXLkiSLl26pMmTJ6tBgwbau3evJk6caBklaHbgwAFJ2ZcAXLJkiZKSktS5c+fCbdg18j1bbMOGDVW/fv1saw9u3bpVH374oZo3b66pU6c6XUF4loK+LsaPH69du3Zp6tSpbjGOHIBr0XMJAPBEbdu2lbe3t37//XeVL19eFSpU0KlTp5SSkiJvb2/NmTNHDz74oJYvX67IyEhFRkbq1KlTSktLkyQNGTJEDzzwQLbzRkZGat++fVkmCzp//rxeeOEF/fzzz0XWPjOnuoE++eQTnTx5Ut99950kadOmTbr77rvVvn17/fnnnwoKCiqQSsKzFNR1sWPHDr3zzjv69NNPs92HCeD6Qc8lAMDTValSRffee68kU4/lkSNHlJKSIl9fX3311Ve677779Oijj1q2Hzt2zBIsX3zxRY0ZM8bmeZ988klJssxPsm/fPt1xxx0aOnSo6tWrV8ityi7fPZeS1LlzZ82cOVMjRozQW2+9pTJlymj48OH6z3/+w/DF61hBXBdnz55Vr169NHHiRMs3DQDQcwkA8FRTp05VUFCQFixYIB8fH7Vq1UrvvPOOWrZsKUmaMGGCvLy8NG/ePKWnp6tly5Z6+eWX1b17d7vn7Nevny5evKhPPvlEP//8sypVqqT333+/yIfDmnkZRuH/qn7iiSd04MABrV271ub2atWq6ejRoyqCqsADpKSk6Pbbb9eTTz6p/v37W55PSEjQuXPnVK1aNZfVDUDRGz1aGjbM9Pi336QePVxbHwAAYBvdi3A7zz//vB5++OEswVKSli5dqpUrV7qkTgDcA59BAgDgvgiXcCsTJkxQrVq19PTTTystLc3ylZqaqmXLlqlOnTquriKAIsY9lwAAeAan7rksCAkJCbp48aIk6fjx41lmOsL15c8//9Qrr7yi9PR0vfHGGzb3ef/994u4VgDcCT2XAAC4L5f1XD799NNq3Lixypcvb1kTsWbNmqpfv756cEPNdWnQoEFKT0+3uz0yMjLb2j4Aij96LgEA8Awu67n88ssvXVU03JR5CmUAsIeeSwAA3Bf3XAIA3Bo9lwAAeAbCJQDAY9BzCQCA+yJcAgDcGj2XAAB4BsIlAMBj0HMJAID7IlwCANwaPZcAAHgGwiUAwGPQcwkAgPsiXAIA3Bo9lwAAeAbCJQAAAADAaYRLAIBbs+65ZFgsAADui3AJAAAAAHAa4RIA4NbouQQAwDMQLgEAAAAATiNcAgDcGj2XAAB4BsIlAAAAAMBphEsAgFuj5xIAAM9AuAQAAAAAOI1wCQBwa/RcAgDgGQiXAAAAAACn+bq6AgAA5MS65/Ljj6WffnJdXeBevL2lXr2khx5ydU3gLs6flz74QDp0yNU1gbv6/nvJ39/VtSi+CJcAALfma/WbatMm0xdgNneu1KmTFBnp6prAHUyZIo0d6+pawJ1lZLi6BsUbw2IBAG6tZ0+pQgVX1wLuKj1dOnPG1bWAuzhxwtU1AK5v9FwCANxatWrS0aPS2bOurgncyeuvS999Z3rMRE8ws74W/vhDatDAdXWBewoIcHUNijfCJQDA7fn5STfc4OpawJ0EB199TLiEmfW1EBnJzw2gqDEsFgAAeByWqIEt1teC9TUCoGgQLgEAgMchXMIWwiXgWoRLAADgcQiXsIVwCbgW4RIAAHgcwiVsIVwCrkW4BAAAHodwidwQLoGiR7gEAAAeh3AJW+i5BFyLcAkAADwO4RK2EC4B1yJcAgAAj0O4hC2ES8C1CJcAAMDjEC5hC+EScC3CJQAA8DiES9hCuARci3AJAAA8DuESthAuAdciXAIAAI9DuIQthEvAtQiXAADA4xAuYQvhEnAtwiUAAPA4hEvYQrgEXItwCQAAPA7hErYQLgHXIlwCAACPQ3CALYRLwLUIlwAAwOPQcwlbCJeAaxEuAQCAxyFcwhbCJeBahEsAAOBxCJewhXAJuBbhEgAAeBzCJWwhXAKuRbgEAAAeh3AJWwiXgGsRLgEAgMchXMIWwiXgWoRLAADgcQiXsIVwCbgW4RIAAHgcwiVsIVwCrkW4BAAAHodwCVsIl4BrES4BAIDHIVzCFsIl4FqESwAA4HEIl7CFcAm4FuESAAB4HMIlbCFcAq5FuAQAAB6HcAlbCJeAaxEuAQCAxyFcwhbCJeBahEsAAOBxCJewhXAJuBbhEgAAeBzCJWwhXAKuRbgEAAAeh3AJWwiXgGsRLgEAgMchXCI3hEug6BEuAQCAxyFcwhZ6LgHXIlwCAACPQ7iELYRLwLUIlwAAwOMQLmEL4RJwLcIlAADwOIRL2EK4BFyLcAkAADwO4RK2EC4B1yJcAgAAj0O4hC2ES8C1CJcAAMDjEC5hC+EScC3CJQAA8DiES9hCuARci3AJAAA8DuESthAuAdciXAIAAI9DuIQthEvAtQiXAADA4xAuYQvhEnAtwiUAAPA4hEvYQrgEXMvX1RUAAADIK+vgMGuWFB3turrAfRw+fPUx4RIoeoRLAADgcbytxl4tXGj6AqwRLoGix7BYAADgcW67TfL3d3Ut4K7uuINwCbgCPZcAAMDj1K4tnTwp7dnj6prA3QQESM2bu7oWwPWJcAkAADxSWJjUvr2rawEAMGNYLAAAAADAaYRLAAAAAIDTCJcAAAAAAKcRLgEAAAAATiNcAgAAAACcRrgEAAAAADiNcAkAAAAAcBrhEgAAAADgNMIlAAAAAMBphEsAAAAAgNMIlwAAAAAApxEuAQAAAABOI1wCAAAAAJxGuAQAAAAAOI1wCQAAAABwWpGEy9KlS6t06dJFURQAAAAAwAW8DMMwCruQixcvKiUlReXKlSvsogAAAAAALlAk4RIAAAAAULxxzyUAAAAAwGmESwAAAACA0wiXAAAAAACnES4BAAAAAE4jXAIAAAAAnEa4BAAAAAA4za3D5YI77tC/69e7uhpuh9cFAAAAgLtx63B56fhxXf73X1dXw+3wugAAAABwN76urkBBid21S+v++1+lJSXpSny8ZBhZtvuHhMi3RAk1/u9/Ve3OOyVJp9as0aYRI5SenGw6xpqXlwJCQuQTGKiWo0apQps2RdMQAAAAAPBAxSZchjVsqHuWLpUknVi+XKsHD5Ykhdx4o27/+msFhIZmO6Ziu3a6d9kySdKhefO0YehQSVLZevV021dfyb9UqaKpPAAAAAB4OLceFptfl44ftzyu1rOnzWB5rYtHj1oeV7/nHoIlAAAAAORBsem5tHZ22zbL44gWLRw7ZsuWq8c0a1bgdQKA69WOHTv0xRdfaM2aNTp69KguX76soKAgVapUSbVq1VKTJk3Utm1btWnTRiVKlHB1dQEAQD4Vy3AZkxkufUuUULn69XPdPy05WXFRUZIkv9KlFVq7dqHWDwCuFyNHjtSoUaMUFhamQYMGqUWLFoqMjNTFixe1fft2ffbZZ1qwYIEkafz48XrppZdcW2EnrVy5Up06dZIkGdfc+w8AQHFX7MLlhX/+0ZVz5yRJYU2ayNs39ybG7tihjNRUSVJ4kyby8i6Wo4UBoEj98ssvGjFihMqVK6dNmzapatWqWba3b99ezzzzjG677TZt2LDBRbUEAAAFpdilqLNbt1oeRzo6JHbzZsvjiKZNC7xOAHA9+uKLLyRJd999d7ZgaVaiRAl98MEHRVktAABQSIpduIyxCpeO3m95xipchjdvXuB1AoDr0aFDhyRJpUuXznG/1q1b64EHHlCtWrWKoloAAKCQFLtwae659HHwfsv0K1cs91v6BAaqbN26hVo/ALhelCxZUpK0aNEipaSk2N3P399fs2fPVo8ePSzPeXl5ZfmaMWOGtm7dqt69e6tChQoKCAhQxYoV9cQTT+jIkSM51mPLli3q27evKleurICAAJUpU0YtW7bU+++/rwsXLuR47Lx583TXXXepfPny8vPzU2hoqJo1a6YXX3xRq1atsuy3cuVKeXl5We63tNWGlStX2m3biRMn9Mwzz6hq1ary9/e3bNu3b1+2/a9tb2BgoM1yzKpVq5Zl+4gRIxQTE6Nnn31WFStWVHBwsBo1aqTPP//ccsy///6rp59+WhUrVlRgYKAaNmyo6dOn5/haAQBQrMLlpZMndfn0aUlSeOPG8vbzy/WY2B07lJH5R0+5Bg3k4+9fqHUEgOtFhw4dJEl///23evTooe3btzt8bFRUlKKiolSxYkVJ0h9//KFu3bqpZcuWmj9/vpYsWaIePXpoxowZatSokdavX2/zPGPGjNEtt9yiVatWadiwYVq1apVmzJihGjVqaOjQoWrcuLH+/vvvbMdduXJFvXv31n333aczZ85o4sSJ2rBhg7799lvVqVNHEyZMUMeOHTV58mRJUosWLRQVFaVp06Zla4P5q0XmaJpr23bgwAF17NhRderU0Zw5c/TLL79YenG9vb0VFRWlxYsX232ttm3bpqjMD0ltWbJkiaKiotQ8c2RObGysevTooRYtWmjBggX68ssvFRcXp+eee07vvPOOzp49q379+qlz586aP3++xo8fryNHjujJJ5/U1KlT7ZYDAECxmtDHekjsud27Nb9Ll1yPSUlIsDyOYEgsABSY1157TbNmzdLp06f1559/qmnTpmrcuLHuvfdedevWTS1atJCPj4/NY+tnjjzxy/yQcM6cOVqzZo1at25t2adDhw4KCgrShAkT9MADD2jPnj0KCQmxbP/xxx/1+uuvq0yZMtq6dasiIyMt2+655x6FhYVp0qRJ6tWrl3bt2pWlLi+++KLmzJmjJk2aaPXq1QoMDLRs69mzpyIiIvTJJ58oNXMyuODgYNWvX1+xsbHZ2pBb28aNG6fVq1erZcuWln0uXLigvn37yt/fX7Vr17b0AttSN5cRN7UzZ0APDg6WJH3++edas2aNbr31VklS8+bNValSJXXo0EFjx47V33//rS+++ELVqlWTZArOAQEBeuqpp/Tuu+/qP//5T47lAQCuX8Wq59J6rcrbZ87UPUuX5voVZPXHButbAkDBKV++vNavX59luOuOHTs0YsQItWrVSuHh4erfv7/++uuvXM/VpUuXLMHSbNiwYfLx8dGpU6csEwhJUlpamoYMGSJJevnll7MES7M333xTkrRnzx7Nnz/f8vzevXst5xo2bFiWYGlmPndB6NKlS5ZgKUm9evXS4cOHValSpQIrx6xZs2aWYGnWtm1bBQYGKjk5WV5eXpZgada1a1dJ0rFjx3IdhgwAuH4Vq3BpXt/Sr1QphdSsmev+KQkJupA54YSXr6/CGjXKc5kJhw9r74wZOrxggdJzuKcov/sDgCerWrWqfvvtN+3atUuvv/66br75Zsu28+fP6+uvv1abNm10//33K8FqJMm12rZta/P5sLAwNcr82W1eL1OS1q9fr+PHj0uSOnbsaPPYihUrWiYbWrZsmeX52bNnW9ao7Ny5s91jly9frvvuu89unR1lq22BgYGqVq2afB1YTiuvmtn4INXb21vlypWzu718+fKWx6czbz8BAOBaxSZcJsfFKeHwYUmmeye9vLxyPSZ21y4p8w+IsjffLN+goDyVeWzxYi178kmFNWyog3PmaPULL+R7/0vHj+vwb7/lqXwA8BQNGjTQhx9+qD179ujIkSOaNGlSllA1d+5cPfzww3aPtw431zL3su3Zs8fy3M6dOy2PO3XqJF9fX5tf5kB77Ngxy/67du2SZAquOc1026lTJ7tLrORFeHi40+fIC3OIvJY5yNrabh1yr1y5UjgVAwB4vGJzz6X1+pbhTZo4dEzsjh1Xj8njkNjUxERtGjFCNe67T1cuXFDM1q2qnDlsKD/7bxo1Sud271b1nj3zVA8A8DRVq1bVwIEDNXDgQG3cuFH333+/Tp48qYULF2r37t2qV69etmNy6sELyvxg0Lrn03oW2F9//VWVK1fOsU7W9zSajy1RooRjDXKSvftOC4u3d86fK+e2HQAAe4pluHR0eGus1Sfbeb3f8uSKFUpJSFB4s2a6oWNH9fj1V5WqUiVf+2ekpipm+3ZVsHE/EQB4qkuXLik9PT3LJDvXatmypcaMGaO+fftKkqKjo22Gy7S0NLvnuHz5sqSs62lal1mhQgW7k+vYYj42KSnJ4WNcLafXBwCAolJsPp60zBTr5aVyDRvmur+RkaG4zKFP8vJSeNOmeSrPHGbDGzeWl5eXQmrUkHcOn6zntH9cdLTSk5IUec2EDgDgyQYNGqTbb7891/2aWI028bezHNSZM2fsHm+eYMY6lDay+pBx3759do89fvy4pk6dqm2Z9+xbHxsbG5vjfaCJiYmWYFvYAgICLI9tDUu1nqUWAABXKRbhMvXSJcVnrlMWUrOm/EuVyvWYCwcPKvXSJcsxAaGheSozdudOlaxcWYF27l3Jy/7mWW4jb7klT3UAAHe3e/duJSYm5rjPyZMnLY8b2vlwcM2aNTafj42Ntdxfeffdd1ueb9Wqlapkjg75LYf72T/55BM9/fTTOnfunOW53r17W4aG/vnnnzaPO3DggEqWLKlnnnkmy/PWw3fNkwJJpgmDFi5caLceuSlXrpzl3KdOncq2fdOmTfk+NwAABaVIwuXs2bP1ySefFNr5Y7Zvl5GRIcnxIbExVvdbOjokNv7AAS168EEt7NVL8X//rbTLl7XowQe16MEHdcbGL/ac9j+1dq3lcfSUKfIJDNSGoUNN///sM4fqAwDuLikpybLkh73tI0eOlCR1795dNe3M9L1s2TKtX78+2/Pvvvuu0tPTVbFiRT399NOW5319fTV27FhJpvUuN27cmO3YqKgoff7552rZsmWWHtY6deroueeekySNHj1aycnJ2Y4dMWKEvLy8NHDgwCzPW088ZB1Yn332WY0ZM8Zm2xzh7++vFi1aSJKWL1+eZVtGRoYmTJiQ73MDAFBQCv2ey9TUVD311FO6fPmyHnvsMbuz1DnDen3LsMaNHTomP5P5hNaqpTt++klx0dFa3KePbn7ySd3cv3++96/Ytq3SU1I0p1Ur3dCxo9p+/LFD9QAAT+Dn5ydJmjhxonbs2KF+/fqpXr16CgoKUnx8vLZu3arPPvtMBw4c0C233KJvvvnG7rlefvllPfTQQ3rhhRfUvn17JSUlaebMmZo6dapKlSql2bNnZ7u388EHH9SJEyc0ZMgQdenSRa+//rpuu+02paWlae3atRozZowiIiI0a9asbOWNHz9esbGx+umnn9ShQwcNGTJENWrU0IkTJ/Tll1/qt99+08cff6xWrVplOe7GG29U3bp1tWfPHo0ePVqPPPKIlixZon/++UeDBw+WJO3fv18pKSlKTU2VZOq5jY6OlqQc7w0dOnSo7rrrLn300UcKDAxUp06ddPnyZU2YMEHdu3e3LKdy+PBhhYWFqVKlSgoNDdXhw4eVmJho6UE+e/asoqOjFRERoYiICLv1qV69uoKDgy11MzOf37wdAAALo5BNmzbNkGR4e3sbQ4cOzdOx39Wta/wzd26u+y3p29f4rm5d47u6dY34gwcdOveCO++0HJN4+nSe6rV/1izju7p1jX83bnR6/zObNxvf1a1r7J81y+HyHX1dAMCV0tLSjNWrVxvDhw83unfvbtx4441GyZIlDW9vbyMoKMioUaOG0atXL+P777830tPTbZ6jatWqhiRj+vTpxsGDB41HH33UqFixouHn52eUL1/eePzxx41Dhw7lWI8dO3YY/fv3N6pWrWr4+/sbQUFBRqNGjYwRI0YYFy5cyPHYX375xejRo4cRERFh+Pr6GmXLljW6d+9uLFy40O4x+/btM+68804jJCTECAgIMGrXrm28//77RlpaWpY22frKza+//mq0bt3aCAoKMkqWLGm0bt3amJv5++Dac02fPt0wDMPo0KGDzbKGDx+eY31WrFhh87zXbgcAwKxQey5TU1M1evRoVa1aVWlpaZowYYJefvnlAu29PLd7t2m9Skle3t4KymEtNLPTf/2li5kTQPiWKKGAsmXzXKa8vFS2bl2n9zcPp2UyHwDFjY+Pj9q1a6d27doVyPlq1qypmTNn5vm4Ro0aafr06fkq895779W9996bp2Nuuukm/f7773a3mycgyo+ePXuqp50lqwyrezytrVy5Msdz5lYfe+cFAOBahRouv/32Wx06dEh33323/v33X50/f17jx4/X6NGjnTrvlvfe04nly5WSkKA0q5n6jIwMzWndWgGhoQoICdEdc+bIJ3Pmwb/efFNnNm5USkKC0q2ml09LStLsW25RQGioSkRGqtusWfLy8sqx/HN79qhU1arys1oXLb/7n928WSUiI1U6cxFwAAAAAPBEhRYuzb2WvXr1UunSpRUXF6dBgwYVSO9l86FD1Xzo0Dwd0/qDD/JdnrX0lBRdOHBAlbt1c3r/9JQUxe7apcpduhRI3QAAAADAVQpttthly5bp0KFDeueddyzPvfrqq0pPT7c5eYKniP/7b2WkpTk8JDan/eN27VJ6cnKWJUj+Xb/eNIwWAK5j0dHRio6OzjbJzNmzZ11cMwAAYE+hhcs2bdpozZo1WRayDg8P15o1a/TAAw8UVrGFLi4qSpJU1mqx7vzuH5c5A19406aSpIy0NO3//nuF1q5dEFUFAI/VoEEDNWjQwLKm49tvv60GDRro008/dXHNAACAPYU2LLZUqVJq27ZttuebZgYpTxUXHS0vHx+Vvflmp/cPyJw2v0R4uAzD0PZx43TTo4/KO3P6fgC4XjGJDAAAnqfQ17ksbuKiohTetKnDk/nktH+1nj3174YN+uv11+Xj769qPXsyaywAAAAAj0S4zIOUhAQlHD6smvffXyD7e/v5qfVHHxVkFQEAAADAJQrtnsvi4vCCBZrTqpWOLV6smO3bTT2MPXoU2P4AAAAAUBy4dbgsWbmygsqXd2kdji1erLSkJAWGhenw/Pm6+cknVSI8vMD2zw93eF0AAAAAwJpbD4u9e9EiV1dBdZ9+Wt5+fjo0d65KREaq/oABBbp/frjD6wIAAAAA1ryMIpiS74knntCBAwe0du3awi4KAAAAAOACbj0sFgAAAADgGQiXAAAAAACnES4BAAAAAE4jXAIAAAAAnEa4BAAAAAA4jXAJAAAAAHAa4RIAAAAA4DTCJQAAAADAaYRLAAAAAIDTCJcAAAAAAKcRLgEAAAAATiNcAgAAAACcRrgEAAAAADiNcAkAAAAAcBrhEgAAAADgNMIlAAAAAMBphEsAAAAAgNMIlwAAAAAApxEuAQAAAABOI1wCAAAAAJxGuAQAAAAAOI1wCQAAAABwGuESAAAAAOA0wiUAAAAAwGmESwAAAACA0wiXAAAAAACnES4BAAAAAE4jXAIAAAAAnEa4BAAAAAA4jXAJAAAAAHAa4RIAAAAA4DTCJQAAAADAaYRLAAAAAIDTCJcAAAAAAKcRLgEAAAAATiNcAgAAAACcRrgEAAAAADiNcAkAAAAAcBrhEgAAAADgNMIlAAAAAMBphEsAAAAAgNMIlwAAAAAApxEuAQAAAABOI1wCAAAAAJxGuAQAAAAAOI1wCQAAAABwGuESAAAAAOA0wiUAAAAAwGmESwAAAACA0wiXAAAAAACnES4BAAAAAE4jXAIAAAAAnEa4BAAAAAA4jXAJAAAAAHAa4RIAAAAA4DTCJQAAAADAaYRLAAAAAIDTCJcAAAAAAKcRLgEAAAAATvN1dQXyyzAMbdiwQXPmzNGWLVv0zz//KC0tTeHh4WratKk6d+6s3r17q0SJEq6uKgAAAAAUex7Zc7l+/Xo1b95crVu31k8//aTIyEg9+eSTGjhwoNq0aaPdu3fr8ccf1w033KBx48YpPT3d1VUGAAAAgGLNo3ouDcPQ8OHD9d5776l58+b6448/1LVrV3l7Z8/IBw8e1Pjx4/Xaa69p7ty5mjt3rsqXL++CWgMAAABA8ecxPZeGYejll1/Wu+++q1GjRumvv/7SHXfcYTNYStKNN96oyZMna+3atTp69Kg6duyos2fPFnGtAQAAAOD64DHh8ocfftD//d//6dNPP9XQoUPl4+Pj0HGtW7fW6tWrFR8fr8cff1yGYRRyTQEAAADg+uMR4TI2NlaDBw9Wnz59NGDAgDwfX7NmTX311VdatGiRZs6cWQg1BAAAAIDrm0eEy6lTp+ry5cuaMGFCvs/Ro0cP3XPPPfroo4/ovQQAAACAAuYx4bJPnz6KiIjI8vyuyZMVNWWKzWOipkzRrsmTszw3ePBg7dmzRxs2bCi0ugIAAADA9cjtw+W///6rf/75Rz179sy2zcvbW1GTJmULmFFTpihq0iR5XTPZT8eOHRUUFKR169YVap0BAAAA4Hrj9kuR7Nq1S5LUpEmTbNsaZN5/GTVpkuX/5mDZYNAgy3YzHx8fNW7cWDt37izkWgMAAADA9cXtw2VCQoIkqWzZsja3WwfM3Z9/rozUVJvB0qxMmTKWcwIAAAAACobbD4sNDAyUJCUmJtrdp8GAAfL281NGaqq8/fzsBktJunz5suWcAAAAAICC4fbhsm7dupKkqKgou/tETZliCZYZqal2J/kxDEO7du1SvXr1CqWuAAAAAHC9cvtwWb16dYWFhWn58uU2t1vfY/nQjh1qMGiQzUl+JNP9m3FxcbrlllsKu9oAAAAAcF1x+3Dp5eWlxx57TNOnT1dSUlKWbbYm72kwYIDdgPnZZ5+pQoUK6ty5c5HVHwAAAACuB24fLiXp+eefV0JCgkaPHp3leSMjw+bkPeaAaWRkWJ7btWuXvvrqKw0cOFB+fn5FUm8AAAAAuF54GYZhFHYhTzzxhA4cOKC1a9fm+xyjRo3SqFGj9Oeff6pjx455OjYhIUHt2rVTRkaGtmzZooCAgHzXAwAAAACQnUf0XErSm2++qU6dOqlnz55atGiRw8edOXNG3bp109GjR/X9998TLAEAAACgEHhMuPTz89O8efPUoUMHde/eXc8995xOnTpld/+UlBRNnz5ddevW1aFDh/Tnn3+qQYMGRVhjAAAAALh++Lq6AnkRHBysX3/9VZ999plef/11ffXVV+rRo4dat26tOnXqyNfXV6dOndLWrVs1d+5cnT17Vn369NGkSZMUFhbm6uoDAAAAQLHlMfdcXuvChQv65ptvNHv2bG3fvl2XLl2SJHl7e+vmm2/Wbbfdpueee86yTiYAAAAAoPB4bLi0lpGRobNnzyo9PV1lypRRUFBQoZQDAAAAALDNo4bF2uPt7a3y5cu7uhoAAAAAcN3ymAl9AAAAAADui3AJAAAAAHAa4RIAAAAA4DTCJQAAAADAaYRLAAAAAIDTCJcAAAAAAKcRLgEAAAAATiNcAgAAAACcRrgEAAAAADiNcAkAAAAAcBrhEgAAAADgNMIlAAAAAMBphEsAAAAAgNMIlwAAAAAApxEuAQAAAABOI1wCAAAAAJxGuAQAAAAAOI1wCQAAAABwGuESAAAAAOA0wiUAAAAAwGmESwAAAACA0wiXAAAAAACnES4BAAAAAE4jXAIAAAAAnEa4BAAAAAA4jXAJAAAAAHAa4RIAAAAA4DTCJQAAAADAaYRLAAAAAIDTCJcAAAAAAKcRLgEAAAAATiNcAgAAAACcRrgEAAAAADiNcAkAAAAAcBrhEgAAAADgNMIlAAAAAMBphEsAAAAAgNMIlwAAAAAApxEuAQAAAABOI1wCAAAAAJxGuAQAAAAAOI1wCQAAAABwGuESAAAAAOA0wiUAAAAAwGmESwAAAACA0wiXAAAAAACnES4BAAAAAE4jXAIAAAAAnEa4BAAAAAA4jXAJAAAAAHAa4RIAAAAA4DTCJQAAAADAaYRLAAAAAIDTCJcAAAAAAKcRLgEAAAAATiNcAgAAAACcRrgEAAAAADiNcAkAAAAAcBrhEgAAAADgNMIlAAAAAMBphEsAAAAAgNMIlwAAAAAApxEuAQAAAABOI1wCAAAAAJxGuAQAAAAAOI1wCQAAAABwGuESAAAAAOA0wiUAAAAAwGmESwAAAACA0wiXAAAAAACnES4BAAAAAE4jXAIAAAAAnEa4BAAAAAA4jXAJAAAAAHAa4RIAAAAA4DTCJQAAAADAaYRLAAAAAIDTCJcAAAAAAKcRLgEAAAAATiNcAgAAAACcRrgEAAAAADiNcAkAAAAAcBrhEgAAAADgNMIlAAAAAMBphEsAAAAAgNMIlwAAAAAApxEuAQAAAABOI1wCAAAAAJxGuAQAAAAAOI1wCQAAAABwGuESAAAAAOA0X1dXAMXPldNXdGnnJaUnpMs7wFtBNwepxI0l5OXt5eqqAQAAACgkhEsUiNS4VJ368pROTz2t5H+Ss233LeOryL6RqjiwooLrBLughgAAAAAKE+ESTjEMQ2d/PKsDgw4oIzFD4X3CVfOjmirZpKT8yvkp/XK6EqMSFb8iXqenndbJKSdV5Y0qqjasmrwDGJUNAAAAFBeES+SbYRj6Z8g/OvHxCYU/GK5aE2rJP9I/yz6+Ib4KqBCgsl3LqtrIajr20TEdffeoLqy+oAa/NZBvaS5BAAAAoDig6wj5dnT0UZ34+IRunHij6v1YL1uwvJa3v7eqDaumxisb69KuS4q+L1pGulFEtQUAAABQmAiXyJeEzQk6MuKIqg6vqkqDKuXp2JDWIar/S33Fr4jXiU9OFFINAQAAABQlwiXy5cDgAyrZuKSqvl01X8eX6VRGN7xwgw6/fVipcakFXDsAAAAARY1wiTxL2JKgixsvqtrwavL2zXoJ7Zo8WVFTptg8LmrKFO2aPNny/6pDq8rIMHR62ulCrS8AAACAwke4RJ7F/hwrv0g/letRLts2L29vRU2alC1gRk2ZoqhJk+TlffWS8w/3V9i9YYqZE1Podc6vS5ekI0ekY8ekVDpYAQAAALuYqhN5dnHLRZW+tbS8fLyybWswYIAkKWrSJMv/zcGywaBBlu1mIa1DFDs/VhmpGfL2c4/POqKjpc8+k5YulQ4ckIzMOYcCAqRGjaT775eefFIKC3NtPQEAAAB3QrhEniUdSlJ473C7260D5u7PP1dGaqrNYClJQXWCZFwxlHIqRYFVAwutzo6IjZUGD5ZmzZLKl5d69ZLefFOqVElKS5P+/lv66y/pnXek4cOlUaOkV16RfHxcWm0AAADALRAukWdGmiEv3+y9ltYaDBhgCZbefn42g6Uky3lcvSTJtm3SnXeahr5Ony498ojkf83KKnfcIb34oimEfvCB9Prr0u+/S/PnSyEhrqk3AAAA4C7cYxwiPIp/hL9STqbkuE/UlCmWYJmRmmp3kp8rp65IkvzK+RV4PR21Z4/UubNUtaq0e7fUv3/2YGktLEz6+GNp5Upp505TKE1OLqraAgAAAO6JcIk8K9m0pC5uuWh3u/U9lg/t2KEGgwbZnORHMt2/WeLGEvINcU0nekqKqZfyhhukJUtMw2Ed1b69tHixtGWLNGJEoVURAAAA8AiES+RZmc5llBidqMQ9idm22Zq8p8GAATYDZkZahmJ+jlFo59Ciqno2EyaYeitnzszf0NZbbjEFy7FjTRMBAQAAANcrwiXyLOzeMPlF+On4x8ezbTMyMmxO3mMOmEZGhuW5mB9jlHIyRRWfq1jodbYlPd0ULh97TGrSJOu2vKzX+eqrph7PzAlyAQAAgOsSE/ogz7z9vVV1WFUdHHxQkX0jVea2MpZtDQcOtHucdeBMiUnRwZcPKuy+MJVqXKpQ62vP6tXS8ePSc89l32Zer1PKWm/rnlkzPz/pP/+R/vc/afJkZo8tDInp6fojLk6bL17U35cvK9UwFOrrq8YlS6pDaKhalColL6+cJ5kCAABA4SJcIl9ueP4Gxf4cqz0P7VHjVY0VfHOww8emXUxT9H3RMjIM1Z5SuxBrmbNNm6RSpaTmzbNvy+t6nZ06mZYm2bdPqlev0Kt+3biQlqZRR47oq9OndSE9XZUDAlQvOFiB3t46kpysebGxupyRocYlS2polSp6ICLC1VUGAAC4bhEukS9e3l6qO7uudnbaqe3ttuumL29S+H321740uxR9Sfse26ekQ0lquLih/CNzmJa1kO3ZYwqC3nYGh+dlvc4GDbKeE85bFR+vR/fuVXxamgbdcIOeqVBB1UuUyLJPumHoz/PnNf74cfXes0cPxMRo6k03KcSXH20AAABFjXsukW/+Yf5qvKqxQtqGaHev3dp15y7FLYpTRlpGlv0Mw9ClXZe0f8B+bW26VRkpGWqyuolCbnXt4pBXrkhBQTnv02DAAMtyKjmt1xkcfPWccN7CuDh13blTN5YooegWLfRBjRrZgqUk+Xh5qVvZsvqjYUP9WLeu/jx/Xp137FB8aqoLag0AAHB94+N9OMWvrJ/q/1JfMT/F6OgHRxXVPUregd4Kbhgsv3J+Sr+crsSoRKWdS5N/BX9VG15NlV+tLO8A13+uUbq0aRhrTmyt12krYMbFXT0nnHPg8mX13r1bd5Qtq9n16snfXteyFS8vLz0YEaGbgoLUaccOPbZvnxbUr899mAAAAEWIcAmneXl5KaJPhMIfDNfFLRd1Yd0FXdpxSekJ6fIN8VWZzmVUqkUplelcRt5+rg+VZo0aSTNmmHobAwKyb7/2Hkvz/yVlC5jbtpn+bdy4cOtc3GUYhp7Yt08V/P31fd26DgVLa41KltSMOnV0T3S0vj1zRo/lZeFSAAAAOIVwiQLj5eWl0i1Kq3QLz+i+a9tWSk2VliyR7ror6zZ763VKshkwf/tNqlRJqly5aOpeXC07f17rEhK0pGFDBedz2t27w8L0QHi4Rh45or6RkfKm9xIAAKBIuE83ElDEGjWSWrSQJk7Mvi0v63WeOyd9+6309NMSOcY5n506pQbBwbq9TJksz+dl3VFJeqVSJf2TnKw/z58vtLoCAAAgK8IlrmtvvCEtXSr99FPW5xsOHGh38p4GAwZkWc/zpZdMa10++2whVvQ6YBiGlsfHq3d4eLZ7Jc3rjl4bMM09zF7XDJ+9tXRp3eDvr+WESwAAgCLDsFhc13r1kh54QHruOdMSInldRuTLL6WZM033bkZGFkoVrxtHkpMVn5am5qVKZduW13VHvby81KxUKW27dKnwK55HuyZPlpe3t80PL6KmTJGRkZHlwwsAAABPQc8lrntffCFVqSJ17Cj9+adjx6SnS++/Lz3zjPT881K/foVaxetCbObyIRVtza6kq0OSoyZN0qzGje0GS7OKAQGWc7qTvPbCAgAAeAr+inE3ly9L589Lycmursl1o0wZaflyqWFDqUsX6YknpOho2/ump0sLF0pt2khvvy0NHSpNmsS9lgXBJ/NFTM3IsLuPo+uOms/j44ZvjHVINgfMnHphUbjyej8vAACwj2GxrnbpkvTdd9Lvv0tbtkinT1/dVrOmacaZBx6Q7r7bdGMfCkXZsqZ7Lz//XBo50jTMtX59qXlz0yyw6emmNTE3bpROnZKaNpXWrpVat3Z1zYuPGoGBkqS9ly+ruZ0FQx1dd9R8nlolShRafZ1hPcx39+efKyM1lWDpIuaeZCnrDNDWgd8TMNwaAOAOCJeukpIiffihNG6cdPGi7X3++cf0NWuWVLGi9N570uOP001WSLy9pQEDpKeekubNMw2R3bZNWrbMtK1GDemRR0xZ/5ZbPORtuBInnV4indsqJR6VZEglKkhlmkrlb5eC3WftlFA/P9UqUULrLlywuT5lXtYdTUxP1/ZLl9QrPLxI6p4fDQYMsATL3HphUXjyej+vuyouIRkA4NkIl67w99/Sgw9Ku3Y5fsypU6bxmrNnm2aQKVu28Op3nfP3N709Dz7o6po4IfGYFDVcOvKDlHFFCq4mlaolyUu6sFvaP9mUjm+4S2owSirT0NU1liT1Dg/XpJMnNa5mTZX0vfrjKa/rjv5w5oySMzJ0f1hYEdY+b/LSC4vCVRx6kotLSC5OiktvcnFpB4CiwT2XRS06WmrbNm/B0trChaaZZ+LiCrRaKEYOfSP9Xl86vVhqOFK677R0z2HptiXSbYulnnulB85LzT+VEvZJi5tL0e9JhuHqmuvZihV1OSNDY44fz/J8XtYdTUxP1wfHjumucuVUzU2HxVr/0f/Qjh3Z7sFE0cvL/bzuKq+TXqFwFZfJu4pLO4qL4nSfeHFqC66i57IoxcVJd9whxcZmfd7bW+rZ09RV1rSpFBIinTkjbdokffON9NdfWfePipLuuUdatUry8Sm6+sP97R0nbR8iVX9cavaJ5B9qez//EKnWs1KN/lL0KGnX29LlY1KLz1w63rdKYKDeqlJF7x87pp7lyumWzHsvc/pU3PoPZ8Mw9N+DB3U6JUWLG7pHb+y18toLi6JRXHqSGW7tPopLb3JxaUdxUZyGwBentuAqwmVReuEF6eTJrM/Vry99/bUpVFqrWFFq0kR69lnp11+lp582BU6zdeuk8eOlV18t/HrDM5yYbwqW9d6SGo52LCT6BEiN3pNK3ihtfFIqVVu6+b+FX9ccDK1aVUvOn9edu3ZpUcOGdif3uZZhGBp2+LA+P31an9eurRuDggq5pvmTUy+seTuKVl7u53V3xSUkFxfFYci1VHzaURwUp7BfnNqCq7wMo/DHwj3xxBM6cOCA1q5dW9hFua9160zDYa21bi0tWiTZWDQ+myNHpPbtJevhgiVKSEePSm48aQmKyJVz0u91pbItpA4L8tf7uO2/pnsxu++QQuoUeBXz4lxqqu7ctUvbL13SiGrV9N/KleWfw9CrQ0lJevrvv7U8Pl5ja9TQq1WqFGFt4cns/SHjiX/g2AvJntSG4mpW48aW0P/Qjh2urk6+FZd2FAfm72/zB0me/H1enNoC7rksOpmfyliEhpom53EkWEpStWrSDz9kDQ1JSdL06QVVQ3iyA1Ok1ATpls/zP6y14WgpMFLa80HB1i0fyvr5aWXjxnqpUiW9ffiwKq9frzf++UcL4+J0LDlZZ1NSFH3pkr7+91/dFx2tWhs36kBSkpY0bEiwRJ7k5X5ed2ZvuDX387qerd5kT1Rc2lFcFIf7xM2KU1tQiOGyuNykWyDtuHRJ+vnnrM+9/bZp6GteymjTRurbN+sOX3+de/mZeE/cS4G1w8iQDn4uVXtECqqY/zJ8S0i1B0pHfzT1hDqosN6PQB8ffVSzpqJbtFCfiAhNPX1aPaKiVHXDBkX+9ZcabNmi/vv26Uhysj6tXVt7WrRQFydmUS4u15VUNG0pLmU0HDjQ7h8yDQYMKJBZMIuiHUURkovLe16U5RT25F204/orw/p8hRX2i/r3YXFqiyfauHGj3nvvPcXExBTI+QotXBaX2cUKpB3bt0upqVf/7+srPflk/sp45pms596zR0pIKLq2uAHacY2LB6TLx6UqfZwvo2of09IlMeuKvh123BwcrAm1aimmTRv907KlFjZooHn162tV48aKb9tW25s317MVK2ZZuiQ/ist1JRVNW4pLGUWhKNpRFCG5OL3nRVFOUfQm047rrwzr8xVW2C/Kn73FqS2eKi4uTh9++KGqV6+uN954w/mQaRSiXZ9+anxXt65xz623Gm3atLH8f9ennxZmsQXu2nrnuR0TJxqGaaEH01ejRvkvIznZMPz8sp5vzZqia4uboB1WDn9vGN/JMJLjnC8jI8Mw5oQZxs7hRd8ON1Bc2mEYRdOW4lJGUaAd7lVGUZSzc9Iku+fa9emnxs5JkwqkHNpxfZVh73yFXU5RvFa5PV9Q5Xjqz9/CFBMTY7z55ptGyZIljeDgYOP11183zp49m69z5XlCn//85z9asWKFw/unJCQoNi5O6ZLK+PrKPyRE/g7O/uhOUhISlHLhgul+NsPIWzvOnzd9mQUFSeXL57+MY8ektLSr/y9f3nTOomiLG6EdmVIvmIaxlqxeMGVcPiH5BEoBYUXbDjdRXNohFU1biksZRYF2uFcZRVlOYaMd108ZKZmj1WydL6dt+S2rMF+r4tSW4iI9PV3x8fG6cOGCvLy8NG/ePN199915Okeex5F169ZN4XmcnfS7Tz7RpbQ03VqmjOpaDQf1NHumTZORni4vH5+8teOvv6TVq6/+v3x505qW+S1j4sSs4bJdO6lmTcfr42g5HoB2SDq3Vfp3mXRz7xwn83G4jINTpeCqUoUueatHXspwc8WlHVLRtKW4lFEUaId7lVGU5RQ22nH9lVEUiks7pOLVlsKQkZGhPXv26K+//pIk1apVSw0aNMjzefIcLnv37q3evXs7vH/UlCn6KyBAZ7y91ScsTA2qVPHIWaCipkxRVFjY1WmS89KOr7/OGi7T0qQPss/I6VAZcXHShx9mfe7tt6U8LBjvVFvcCO3IdOoPaeVS6a6npVK2P2RwuIyMVOmn8VLjx6Q6LxZtO9xEcWmHVDRtKS5lFAXa4V5lFGU5hY12XH9lFIXi0g6peLWloKWlpWnWrFl69913tX//ft19991655131KxZs/ydsEAH7F6Dey4zRUVlvUdSMj2XnzJmzsx6nsBAw0hNLbq2uAnaYSXpjOmey0PfOF9GzAbTuc44fh9vgbXDDRSXdhhG8bi3qKjKKAq0w73KKMpyChvtuP7KKArFpR2GUbza4qgvvvjCuPnmm41KlSoZjRs3NubNm2d33wkTJhiSjLvvvtvYsmWL02UXWri0fuP69+9vtGnTJtvznqBAbjZOSzOM8uWzhsL+/fNeRlqaYTRtmvU8XboUbVvcAO2wYWl7w1jawfkyNjxtGHNvMIz0/H9gkWsZbqq4tMMwiqYtxaWMokA73KuMoiynsNGO66+MolBc2mEYxastjnr11VcNLy8vY8GCBUZcXJxx5513Gt7e3sbu3btt7n/+/Hnj77//LrDynZu7P6ceUes1tzZtsjxv7oL2lIWpc1o7zLw9Vz4+0n/+I40effW5GTOkRx6RunRxvIzx46Vt27Ke+7nnirYtboB22FBroLSuj+ney/Kd81dGwn7pyEyp3lDJ2/EfDbwf7qco2lJcyigKtMO9yijKcgob7bj+yigKxaUdUvFqiyPWrl2rcePGqXXr1rrrrrt0//33a+HChZKkffv2qW7dutmOCQ0NVWhoaIHVIc+zxebHE088oQMHDmjt2rWFXZT7OnVKql1bSky8+lxIiLRwodS6de7Hf/ON9MQTkvU3Qa1a0u7dkp9fwdcXnsXIkJZ1lhIPS3dskwLK5u349BRpWScp+ax0507J1/HZhwEAAOB6vXv31pw5czRkyBCNGTNGt956qzZu3KiKFSsqOjpaZcqUKfQ6sHJoUalYUfroo6zPXbggtW8vDR0qnTtn+7gjR6THHpMefzxrsPTykr76imAJEy9v6dZpUtolacUd0pU4x49NTzb1ep7bIrX6mmAJAADgYa5cuaLff/9dktS8eXNJ0ty5czVr1ixt3769SIKlRM9l0crIkO67T1qwIPu2wEDp9tulpk1NPZpnzpiGE69aZbq78lrDhkmjRhV+neFZzm2XVnSVvP2lWz6XKvbIcXkSxW6SNj4hXfxHavezdEOPoqsrAAAACsSaNWvUvn17SdLBgwdVM4/LFBaUQrvnEjZ4e0uzZkkPPGAaDmstOVn67TfTV25eeUUaObJw6gjPVraJ1H27tPFpadVdUpnGUo0npHItpVK1JC8fKfGIFLdZOvKddHalVKapdMdmKTTvaxkBAADA9bZs2SJJCgkJUY0aNVxWD8JlUStRQpo3T3rvPdNXWprjx4aESP/3f1K/fjn3RuH6FlRJ6rhQOr1Y2j9Z2vZfybj2OvOSwttKrb+XqjwgeTO8GgAAwFPt2LFDktS4cWN5uTAnEC5dwc9PGjFC6tVL+vBDac4cKTXV/v6lS5vuuXzjDdO9m0BuvLykineYvtKTpfgoKfGoJEMKLC+VaSL5lXR1LQEAAFAAdu7cKckULl2Jey7dwZkz0tKl0tat0oED0pUrUnCwVK+e1Ly51KWLVJIgAAAAAEDau3evevTooYyMDKWlpenkyZOSpPLlyysgIECS5OPjo9WrV+uGG24osnoRLgEAAADAQ7nLZD5SES1FUrp0aZUuXbooigIAAACA60ZUVJQkqWTJki6dzEcqonsuR48erZSUlKIoCgAAAACuG7t27ZIkNWjQwKWT+UhFFC5LlSpVFMUAAAAAwHXF3HPZsGFDF9ekiIbFAgAAAAAKXnR0tCTCJQAAAAAgn44cOaKEhARJhEsAAAAAQD6Zh8RKhEsAAAAAQD6ZJ/OpVq2aW6zOQbgEAAAAAA/kTpP5SIRLAAAAAPBIO3fulES4BAAAAADkU1JSkg4cOCCJcAkAAAAAyKfo6Gilp6dLIlwCAAAAAPJp+/btkqRSpUqpVq1aLq6NCeESAAAAADyM+X7LW2+9Vd7e7hHr3KMWAAAAAACHmXsuO3Xq5OKaXOVlGIbh6koAAAAAAByTkZGh0qVLKzExUVFRUapfv76rqySJnksAAAAA8Cj79u1TYmKibr75ZrcJlpLk6+oKOGLTpk168MEHlZiYqLi4OF3b2Vq2bFkFBwdrzJgxeuihhyRJixYt0jPPPKPLly8rLi4uy/5eXl4qW7asgoKCNHXqVHXt2rWwGyA9+KCUmCjFxUnXdhaXLSsFB0tjxkiZ9deiRdIzz0iXL5uOydoA0zFBQdLUqVJh119SwqYE7X5wtzISM5Qalypd0wTfsr7yCfZRjTE1FPlQpCQpblGc9j+zX+mX05UWl3ZNGzKPCfLRTVNvUtmuZQu1/h7/FsRuktY9KKUlSlfilO0N8C8r+QZLjcdI1TIbcGqRtOkZKf1y5jFZWiAFlJV8gqSWU6UKhduATQkJenD3biVmZCguNfXa2qusr6+CfXw0pkYNPRRpun4WxcXpmf37dTk9XXFpWa8fr8xjgnx8NPWmm9S1bOFePwAAAK701ltvadKkSeratavmzJmjNWvWSJIGDhzo4ppl5XHDYufPn697771XklSvXj2tWrVK5cqVy/GYr7/+Wv3795ckNW/eXH/++adCQkIKuaZ2zJ8vZdZf9epJq1ZJudRfX38tZdZfzZtLf/4puar+kmLnxyr63mhJUlC9IDVZ1UR+5fxyPObfr//Vvv77JEmlmpdSoz8byTfENZ9tePxbcGK+tPpe0+OQetLtq6SAXBpw6GtpQ3/T47LNpdv+lPxd04D5sbG6N9p0/dQLCtKqJk1Uzi/n6+frf/9V/32m66d5qVL6s1Ejhfh6xGdjAAAATtm7d6/q1q0rSQoLC1NMTIxuu+02HT58WLt371ZQUJCLa3iVx/11dujQIcvjvn375hosJVkWF5Wkfv36uS5YSpJV/dW3b+6pRpKs6q9+/VwaLCUp6VCS5XFk38hcg6UkXT5w+eox/SJdFiylYvAWXLJqQLW+uQdLSbpo1YDq/VwWLCXpUNLV66dvZGSuwVKSDly+ev30i4wkWAIAgOtGmTJl5O/vr6CgIE2YMEG//PKLVq9erUWLFrlVsJQ88J5LcxewJHXs2NGhY1atWmV53L59+4KuUt5Y1V8O1l9W9Zer6y/pwpoLlsehHUMdO2aV1THtHTumsHj8W3DWqgERHR08xqoBEa5twJoLV6+FjqGhDh2zyuqY9g4eg4LTsWNHeXl5Zflq2rRpns6RnJysihUrZjvPjBkzCqfSRWzcuHEqVaqUxo0b5+qq2LRy5cpsr/2RI0dcXS0AgAPKly+v3377TY0bN9awYcP01ltvae7cubr99ttdXbVsPG5YbEREhGJiYhQcHKz4+Hj55tKDkZSUpNDQUKWkpCg0NFRxcXGuXQcmIkKKiTHd4BcfL+XWA5OUJIWGSikppn/j4iQXr2OzLmKdUmNS5R3srbbxbeXtm3N90pPStTZ0rYwUQ76hvmoT10Ze3l5FVNvsPP4t+DlCuhJjusfygXjJO5cGpCVJc0KljBTJL1R6IE7ycl0DItatU0xqqoK9vRXftq18c3kxk9LTFbp2rVIMQ6G+vopr00beXq67fq5Hhw8fVmJiojZv3qwnn3zS8vz8+fN19913O3SOCRMm6MUXX5QkVaxYUYsXL5YkVapUSaHF4AODBg0aKDo6WvXr11dUVFS27SNGjNDIkSPVoUMHrVy5ssjrl5iYqMOHD+vUqVPq1q2bJNP7Wq1atSKvCwCg+PKonsu9e/cqJiZGktS6detcg6UkrV+/XikpKZKkNm3auDZY7t1rSjWS1Lp17qlGktavN6UaSWrTxuXBMnFvolJjUiVJIa1Dcg2WkpSwPkFGiukzjNJtSrs0WHr8W3BhrylYSlJY69yDpSTFrjcFS0kKb+PSYLk3MVExqabrp3VISK7BUpLWJyQoJfMzsDalSxMsXaB69eqqX7++qlevLkmWn72jRo1y6PgrV65ozJgx8sscAu3n56f69eurfv36xSJYStI777yj5s2ba9iwYa6uik3BwcGqX7++ateu7eqqAACKMY8Kl6tXr7Y8dnRIrPUnxO3atSvgGuWRVf0dHo9p/Qm3q+sv6cLqvA+JjV8Zf/WYdo4dU1g8/i04a9WAyI4OHrPy6uMI1zZgdT6GxK6Mj7c8bldMgoine/TRRyVJW7du1W+//Zbr/lOnTlViYqJ69uxZ2FVzmd69e2vz5s168MEHXV0VAABcxqPCZX7ut7QOlx55v6V1snF1/SXFr4m3PM5PuAxp79rJiDz+LYjJx/2WZ1ZefRzu4vstrYJifsJlexdPZgWTvn37Wnoxc+u9TElJ0UcffaQXXnhBpUuXLorqAQAAF/HIcBkUFKQWLVrkun9ycrI2bdokSSpRooSaNWtWqPXLlTnZBAVJDtRfycmmBRolqUQJydX119XJfLyDvFWqRalc909PTlfCpgTTMSW8VapZ7scUJo9/C8yT+fgESeUcaEB6shSX2QCfElJZ1zbAPJlPkLe3WpTK/VpITk/XpgTT9VPC21vNHDgGhc/X11dvvvmmJGnz5s1auHCh3X2nTZum+Ph4y/2W9hiGoeXLl2vw4MFq2rSpSpcuLT8/P0VGRqpHjx6aN29ervXauXOnevfurYiICAUEBKhKlSrq37+/Dh48qBEjRmSZzMZ8r6G956dOnapmzZopODhYISEh6tq1qzZs2JCtzBkzZmSbKMea+fwjR46UZJpgztb+n332mc16mG3YsCHHcqxdunRJw4YN00033aTAwECFh4frzjvv1IoVK3J9DSXTezF79mx1795d4eHh8vf3V0REhLp166ZvvvlG6enpNo87ffq03nzzTTVq1Ehly5ZVYGCgatSooQcffFDTpk3TBauRCwCA4sljwuWRI0d07NgxSab7Lf0cWL5g/fr1unLliiSpZcuW8vf3L9Q65ujIESmz/mrdWnKg/lq/Xsqsv1q2lFxZf0lJR5J05ZipPiGtQ+Tt5+D9llcy77dsWVre/q675Dz+Lbh0RLqc2YDw1pK3Aw2IXS9lZDagXEvJx3UNOJKUpGOZL2brkBD5OXi/5ZXM+y1bli4tfxffc4yrHn/8cVWuXFmSLMHpWqmpqfrwww/1/PPPq2zZsjme7+jRo+rcubM+++wzde7cWfPmzdOaNWv04Ycf6siRI7rvvvtyXCj6l19+UYsWLTRv3jz1799fK1as0I8//qjKlSvr1ltv1bZt2ySZ1jqOiorSkiVLJEnPP/+8oqKiNHr0aMu5Bg0apL/++kv/+9//tHDhQj3xxBNaunSpOnbsqF27dmUp995771VUVJSmTZtms17m8w8YMCBL+dZfktSnT59s9bDWqFGjHMsxO3v2rG699VaNHj1aFSpU0OzZs7Vo0SI98sgjev755/Xtt9/mePyVK1fUu3dvPfjgg0pISNDkyZO1Zs0affLJJ4qLi9Pjjz+ubt266bLV8kCSFB0drXr16mnKlCnq16+fFixYoJUrV2rIkCHasGGDnnrqKb300ks5lg0A8Hwes1ic9ZDYrVu3OjTD3fnz5y2P3WpI7NatkiMz9FnV3/XjMbMuQXJx60Wtr7Y+12PSzqdZHrvTkFiPfAush8Se2yrNr5b7MSlWDXCjJUi2Xryoautzv37Op129fhgS6178/f312muvafDgwdq0aZMWLVqkO+64I8s+M2bMUExMjP773/86fN6PPvpIr7zyiuX/t956q3r37q2mTZvq008/1R133KG77roryzGHDx/WY489ptTUVE2YMEGDBw+2bGvVqpVuuukmPfbYY5KuTmxjFhERoYiICG3ZskWSdOLECSUmJmr69OmWfTp06KDTp0/rp59+0tixYzVz5kzLttDQUIWGhio2NtZme8znj4iIsFm+WZkyZVSmTBlLPa5VokQJ1a9f3245Zk888YR2796tli1baunSpZYPYps1a6Y77rhDjRo1yvH4l156ST///LPatWunFStWyMfHR5LpA9oHH3xQLVq00LJly/TKK6/os88+sxz3zjvv6Pz58/riiy/09NNPW56/9dZb1aVLF9WrV08eNjk9ACAfPCZcWk/ms2bNGtWrVy/XY+rVq6c9e/ZIcmwyn6NHj6pdu3aWHtICZT2TzJo1kgP1V716Umb9c5xJ5tdfpfHjpdRU6dw56bbbpNdekzJ7FQqK9WQ+TdY0UXC94FyP2VRvky7vMX3CHdLOdjj4541/5F/eX2W7lVVA5QAZaYaSDycrdn6sAqsGqsITFQqk/oXxFiQkSPXrS599JtWpY1qqJCREyvx7rGBZT+Zz+xop1IEG/F5PupDZgJwm84ndJEW/K6Vdkow0Kaiy1GSMFFTJuTpbsZ7MZ02TJqoXnPv1U2/TJu3J7CFp50C4PJqcrHbbt+tYq1YO1Smv+yOr//znP3rvvff077//auTIkVnCZVpamj744AM9++yzCg8Pz/VcoaGhGj58uJ566qls20qWLKm+fftqxIgRmjlzZrZwOXbsWCUmJioiIsLSQ2jt0Ucf1ahRo3TgwIFc65Genq5XX3012/NdunTRTz/9lOWDTnezadMmyxDloUOHZhvhExYWpoEDB2ro0KE2j9+3b58+//xzSdJ7771nCZZmvr6+GjJkiPr27atp06Zp5MiRioyMlCTL79qSJUtmO++NN96oZ555RjVr1nSugQAAt+cxY8zMv9BDQkJUt27dXPePj4/X3r17JZl+IbbK4Y/HCxcu6Msvv1Tz5s11/Pjxgqnwtcx/kISESA7UX/HxpnUzJNN6Gfbq/8UX0rhx0vffm8pYu1ZaulRq0EDasaMgan61SpmT+fiE+CioblCu+6fGp+ryXlMw8PL1Ukgr2+EgYUOC/nn5H22uu1lrS63VujLrtLXpVp3/87zCe+f+R6mjCuMtOHBAOn5c6tFDqllTKlfOtO+1Xw500uXO3HPpFyKFONCAlHjT0iWS5OUrhdm5ho7NkZZ1lG4aLN2+QuqyRipZXfqzg5QcUwAVNzFP5hPi46O6QblfP/GpqdqbGSx9vbzUKodweSEtTV+eOqXmW7fquHkccw7yuj9sCwwMtASxDRs2WIaaStI333yjU6dOaciQIQ6dKzQ0VCNGjFCInfe5atWqkmT5uW5twYIFkkwjVOwtUeXobOFBQUE2f8fccMMNkqR///3XofO4gvl1kKROnTrZ3Cen12H27NkyDEOBgYG69dZbbe5Tp04dSaYhz9Yf+pqXOHn99de1ePHibL2UEydOZFgsAFwHPCJcnj17Vn///bck09CcnCYyMNu4caPll1vTpk0VbKeX5K677tI999yjQ4cOFVyFr3X2rJRZf7VsKTmyTt/GjZL5l3PTppKt+ickSO++K82dK5Uvb3quTBlp4EDpwgWpAH+Rp5xNUdLfSZJM90468h5c3HhRymxCyaYl5RNsvzvPt8zVPwhL1C6hmuNrqsmqJvItWTCd64X1Fhw8aPvY9PSrX9262f9swGHJZ6WEzAaUc7ABsRtleQPKNpV8bTTg8glpwxPSDT2lCl2vPl/3DSnxmBQ13MmKm5xNSdHfSabrp2Vpx66fjRcvmmuvpiVLKthOd/BdUVG6JypKh5KTHapLXvdHzp577jmFhYVJunrvZXp6ut5//3099dRTqlDB8ZEHZ8+e1dtvv61bbrlFZcuWlb+/v3x9feXr62vp0bx06VKWYxISEnTy5ElJyvF2ifLmn5G5KFu2rM3rs0SJEpJkuY/fHZl7D8PCwmz2IEo5vw47d+6UZJoMr0SJEpbX3vrrlltusexvPcrngw8+UGRkpI4fP6477rhDVapU0fPPP6/ffvvNrV8zAEDB8ohhsdbDkFq3bu3QMX/99ZflcU6f1P7666+Wxz/88EOu97Pki/UwKgfrL6v62x0Su26ddOKE9Mgj0qJFVwNHrVqmfzduzHtd7bC+3zKktWP3vl34y+oYO0NizdrEtVHauTR5B3rnGELzq7DegoMHpalTpYcfNk324+Nz9W1ISpI6dZK+/DJ/dc7irFUDwh1sQKxVA8LtNODviaahsNdu9ytlCqSHvpaafCz5lshbfa9hfb9lawfvnfzL6pichsT+2qCB5fEPZ84ot+/gvO6PnAUHB+uVV17RW2+9pb/++kt//vmnTp06pWPHjun11193+Dzbtm1Tly5ddO7cOXXu3FlffPGFqlevroCAAEnS/Pnz9fbbb2frEUvInE1YuhoAbXFkEjhJ2YaCehLza5Hf18E8m2tkZKT+/PPPXMszD4mVTLehREdHa/Lkyfr66691+PBhTZkyRVOmTFFoaKgGDx6st99+27UT6wEACp1HhEvroTc5DW+1tt5qHKLLJ/OxvtnP0S4s63GU9upvnq1vyRIpLk7K7D1QSorp3wJcUy5+dbzlcelWjp03Yf3VP/pC24fmuK+Xl5f8yjn2x19+FNZbcPCg9OyzpqVNrvXGG9Lzz0sVKzpeT7us77e0N7z1WrFWDbA3mc+p303/lqyRfVupWqZlTM6ulCp2d6xMO1ZbrVXZysHrcr1VaGjv4JqYcI1BgwZp7NixOn/+vIYPH67Y2Fg9/vjjqlKlisPneOqpp3Tu3Dm1adNGS5Yskfc1MwPbm+jGeu3Ma2cwtZaamupwXdxZmtUkV9cyvxb5fR3MQ5KTk5NtTjqUm7CwMA0fPlzDhw/Xtm3b9PPPP2vmzJk6fvy43n33XR04cEA//PBDns8LAPAcHjEs1txz6eXlpZYtW+a6f0ZGhjZm9tp5eXmpbdu2hVq/XJm7zby8TGMyc5ORcbXX0ctLslf/7t2lnj2lF164Giwlafdu07+33Zb/Ol/D0nPpZRoWmxsjw1DCxgTLMSFt3WOm2IJ+Cx5+OOtLb7ZypWnpk3798lNbGywzxXqZhsXmxsiQ4sw9115SuJ0GJB4x/esTmH2bf+bSEee25qGitpl7Lr1kGhabmwzD0MbMcOklqS0zxbq1UqVK6YUXXpBkGjXyzz//WNbBdMS5c+e0I/Me8XvuuSdbsMxJ6dKlLfdDHjlyxO5+7nyvpDVzT629oaQ5ja4x3ysaFxenixcv2twnp9fBPJPshQsXctxv06ZNmjp1qk6fPm13n6ZNm+q9997ToUOHLNfGrFmzCm9eAwCAWyiycJnfKcgTEhIs94HUrVvX7mQP1nbv3m0ZHlS3bt1c11crVAkJUmb9VbeuaTaZ3OzebTrOfIy9+gcFmWaK/b//y/r8L79IgYHSsGH5r7eVtIQ0Xdppus8pqG6QfENy7/BO3J2o9IR0yzF+ZXPulUz6J0nRD0Rre4ft2tpyq/b226srJwvmPp3CfAu6ds3+3JUrpttdrWbpd05qghSf2YCQupK/Aw24sNt0nPmYAHvfAznc++iV+T5fPuFwVW1JSEvTzsz75OoGBSnEzoQr1nYnJiohc6H2ukFBKuvgkEa4zosvvqhSpUpJkh555BHVqGGjN9yOjIwMy2N7vytyCo733HOPJNMoF3s9e2vXrnW4PoXBPNGQdfsMw9CsWbOy9Mqa74mMjY212cu4adMmu2WYXwdJWrFihc19cprttnfv3pZg/9tvv9ndb8CAAXrhhReyzGXQokULvfXWW9n29fX11ahRoyz/zymQAgA8X5GEy9atW6tq1ar5Cpjr1q2z/OHh6JBY6/stXT4kdt06UzeY5Ph4TOub/fJa/08/lfbtk+bMcWxKVAdcWHdBymyCvRlfr5Xwl+NDYiXpwOADqjnWNIlP0/VNZaQb2tZmW4EEzKJ+C/73P1PozMM8JjmLWWfqiZQcHxIbY9WAnNa3DK5m+jc9Kfu2K5k9JCnnHCvTjnUXLpgvnxxnfLX2F0NiPU6ZMmX0/fffa+zYsVnChCPCwsIss5D++OOP2QLipUuX9M0339g9/tVXX1VwcLBiYmI0ZcqUbNu/++47HbQ3+1YRMYfGc+eufj+dPHlSDz/8cJYg16xZMwUGBiotLS3LLSGSFBMTo++++85uGS1atNCdd94pSXr//fezhdPY2FhNnDjR7vF16tTRc889J8m0FElcXFy2faZNm6Zt27Zp8ODBWYYkx8TEaM6cOUpKyv6zZN++fZJM9+c6Mts7AMBzFXq4PHjwoDZs2KDjx49r2bJleT7e2fstHZ1+vtA4e7OfI/VPS5O6dJGaNTPd6Pd//ydl/oFREKzXt3T0fssL6x2fzKdcz3K6+fubVaK6aRIKL28v3fjJjUo5maKDLzv/B2FRvAVmMTHShx9KeVgzPnfO3m9pbzIf6eq9lEk2ehPid5j+TYl3rEw7rNe3dPh+Swcn80HROHnypKKjo3X48GFJ0uHDhxUdHa2zZ89m2a9nz5569dVXs83aun//fkVHRys+897b1NRURUdHKzo6WomJiZKkSZMmyd/fX9u2bVOnTp00d+5cbd68Wd99951atmxpCTrmY/fv3285f/Xq1TVz5kz5+fnplVde0euvv67169drw4YNeueddzR8+HA988wzNtsWHx+v6Ohoy4yz154/MTExS9slWepufbyt7dbhrlu3bgoICNCePXs0ffp0bdy4Ua+88op8fHyyrNtZunRpDR48WJLUt29fzZgxQ1u2bNGcOXPUvXt3yxBTe+VMnz5ddevW1caNG9W1a1f99ttv2rZtm7777ju1bdtWnTt3zva+WM/AO378ePXp00dHjhxRy5YtNW3aNG3dulWLFy/W888/r2eeeUbdunXL9gGCl5eXDhw4oPbt2+vbb7/Vxo0b9ddff2ny5Mm6//775e3trU8//dTuLLYAgGLCKGSPP/64UaJECSM4ONho06aNkZGRkafj27RpY8i0noKxZ88eh46pXbu25Zjjx487XFbVqlWNAn9J2rQxDNOKFobhYP2N2rWvHpOH+huGYRgnTxpG/fqG0aqVYZw4kff62rC1zVZjhVYYK7TCuLTnkkPHbKi9wXJM0vGkfJW7qf4mY4XXCiP5ZHK+jjcryrfg/fcNo2HD/NXTriVtDOM7mb7iHWzAgtpXj0nMoQGJxw1jVpBhrH0o6/Pnowzj5wjT8avuy3/dDcNos3WroRUrDK1YYey55Nj1U3vDBssxx5Mcv36q/vWXoRUrCm3/69Xjjz9u+Zlq/TV8+HCHjjf/bLX1tcLq9d+xY4fRp08fIzIy0vDx8TFKly5ttGrVypg0aZLx5ZdfZjmuatWq2crZsWOHcf/99xthYWGGv7+/UbVqVeO5554zzpw5YwwbNsyQZHTp0iXLMdOnT7dZL/P5V6xYYbfuOR0vyTh8+HCWsn7//XejRYsWRmBgoFGqVCmjefPmxty5c7O1Iz093Rg7dqxRp04dw9/f3yhXrpxx7733Gjt37rRZn2vLSUhIMIYOHWrceOONhr+/vxEaGmq0b9/e+P77743Dhw9nO/6PP/7IVof58+cbPXv2NCIiIgxfX18jNDTU6NChg/HVV18Z6enp2fY/fvy48eGHHxqdO3c2KlSoYPj5+Rn+/v5G9erVjccee8zYsmWLjSsDAFDceBlGPm+GdMDBgwdVp04dNW/eXAkJCdq7d6+WLl2q22+/3aHjt27dqpYtWyo9PV3e3t6Kj4+33NNjz9KlS9U180a44OBgnTt3zuGpz6tVq6ajR4/m+/7QbLZuNc0ek54ueXtL8fFSLvXX0qVXb+QLDpbOnTOtcZEXv/9umuinfn1p82bT/Zf5dHHrRW1tuVVKl+QttY1vK99SOd8zd27pOe3qukuS5B3srbbn2srbP++d5Ns7bteFVRdU55s6Kv+YY2vUXaso34L0dKlGDalHD9Po5AJxbqu0uKVkpEte3tID8aZlQnJyeqm0IrMBvsHS/ecknxwacPRHaX0/qe1PUqV7TGtqbn1JCgiX9k+Qag2QWuSvQVsvXlTLrVvNl4/i27ZVqVzuuVx67py67jJdP8He3jrXtq38HZzgpdr69Tp65YqMjh0LZX94rhdeeEETJ05U37599e2337q6OgAAFEuFuhTJ6NGjFRERoVq1aunQoUNq0aKFRowYoc6dO9tdRH3w4MGaN2+e4uPjswzVycjIUNmyZRUWFqayZctq+/btltDYr18/LV++XOfPn88yBXtiYqJKlSqlcuXKqVKlStq4caNDi7c7ZfBgad48U4qxXuw7I8M0K0xYmOnf7duvJpZ+/aTly6Xz568uL2JqgCkJlSsnVapkmr7UkfrffrtpwcXoaGnWLKl//zw14cDgA4qdF6u0+DSlX0q3aoO0ruw6+YX5ybesr5pvb24JjXv77dX55eeVdj5NGZevTs6RkZihNaXWyK+cnwIqBajpxqZZ3oOTk0/qzA9ndPPMq8Nir5VyKiVP9XfVW7Bhg3TsmFSzZp6qm92WwdKJeabhqGn/3979x1pd33ccf/FDQCwKVAurlQtSadcfol7UmqzWH2lnJ96aJWq1q21N6VrYtP+safojpbVJ1zmDi5qsIUukLZLorHQNoa2AZdplBaGZtkPARileA0IBuQoUrtz9cctVymVwfHPvuZTHIyH3m3O+55z3NZj49PM9n+/rfoGu/cm/j02Gn969Qc9Vv3wtGv/r5mTzsmTv9uTV1/0Cna8kD45Khr85OfltyV/28gu03JCcMilZc0f3n5FnJef9Y7Lmn7ufH3thQ+P//fr1Wbh1a3Z0dublV1/7+7M/ydif/zynn3RSxg4dml9Om9YTjTevWZNl27dne2dndr1uc5dX9u/PqMcey5tPOilvGz48v7jggr7/d5jj0p133pkPfOADmTZtWq/PH9g0p+m7hwPAn7A+W7k8sGo5Z86crF69OuvXr8+XvvSlXH311Q2tXvanY75y2dfuuit54IHk298+9IuB48cnmzcns2Yl99zTlPGOxoo/X5FdT+/K2f90dib8w8H3xFt53sq88j+vZMq/Tslb//ZY3Cyyb33jG8nXvpbMm3cMb0HSTP/518nzDydtzyZvmtjsaY6KlcsT18SJE/POd74zixcvPuR/QPzsZz/L5ZdfnnHjxuWZZ57xvT8A6CN9tqHP9u3bc9lll2XGjBk9j334wx/ORz/60cPef4sGffnL3TvP3Hnnoc8d2LFvgG+GcmAVdMyVYw55bu+m7hXLUdOOcBnoAHFg46AB/o/8UHteTPb/0S0PurqSbU8k4648bsISfvKTn6StrS0LFy7M6tWr89hjj2X27NmZPn16xo4dmx/84AfCEgD6UJ9dFnvhhRdmyZIlBz02aNCgLFiwoK8+8sRz5plJe3ty3XUHP/7CC6/dpPGqq/p/rgaM+eCYvP1f3p5RFxwckK+seSX7Nu/Lm1rflFGtx0dc/va33T9HjmzuHA3Z+t/JI3+RjP9gcvni1x7fvDTZtTG5+N+aNxs04KGHHsrDDz+cpUuX5tZbb82LL76YIUOGZNKkSfnc5z6Xz3/+8znzzDObPSYA/Enr0+9cHk927tzZs6K6cePGnHXWWU2e6CjcdluyfHlyww0HP37//d0/r7++sftoNMFbZ741T9/8dN794Lsz5JQhPY9vum9TBo8YnHfMfUcTp2vMli3dP4+wX83A0rmr+7ucI1/3933vS90b+kz5u+TPPti00Rq1s7MzHX/4jufGPXty1hE2smr0fAa21tbWtLa25pvf/GazRwGAE9bx9J/BfWLGjBlZuXJl1q1b13Pz58mTJ2fKlClpaWnJokWLmjzh/2PWrO4dZ97//uRjH0vOOSdZsyb5+teTz342mTOn2RMe0bDTh2Xi7Il5avpTGfc34zLszGHZvmR7tj60NVOXTs2o84+PVcuk+0rSJBk9uqljNOYt70/edm33hj/P/0f3hkD/++3kzRcm5/dyufUANGPt2qzcuTPrdu/O7j9sBjT5F7/IlJNPTsuIEVl07rml8wEAODp9eiuSAz71qU9l/fr1efzxx/v6o05Mv/td8tOfJhs2dG9r+qEPJS0tzZ6qIZ0vd2brw1uzd9PejJwyMmP/amwGn9RnXwnuE8uXJ889l3ziE82e5A3Y/Giy7ZfdO8m+5bJk7PnNnggAgOOMuAQAAKDs+FoaAgAAYEASlwAAAJSJSwAAAMrEJQAAAGXiEgAAgDJxCQAAQJm4BAAAoExcAgAAUCYuAQAAKBOXAAAAlIlLAAAAysQlAAAAZeISAACAMnEJAABAmbgEAACgTFwCAABQJi4BAAAoE5cAAACUiUsAAADKxCUAAABl4hIAAIAycQkAAECZuAQAAKBMXAIAAFAmLgEAACgTlwAAAJSJSwAAAMrEJQAAAGXiEgAAgDJxCQAAQJm4BAAAoExcAgAAUCYuAQAAKBOXAAAAlIlLAAAAysQlAAAAZeISAACAMnEJAABAmbgEAACgTFwCAABQJi4BAAAoE5cAAACUiUsAAADKxCUAAABl4hIAAIAycQkAAECZuAQAAKBMXAIAAFAmLgEAACgTlwAAAJSJSwAAAMrEJQAAAGXiEgAAgDJxCQAAQJm4BAAAoExcAgAAUCYuAQAAKBOXAAAAlIlLAAAAysQlAAAAZeISAACAMnEJAABAmbgEAACgTFwCAABQJi4BAAAoE5cAAACUiUsAAADKxCUAAABl4hIAAIAycQkAAECZuAQAAKBMXAIAAFAmLgEAACgTlwAAAJSJSwAAAMrEJQAAAGXiEgAAgDJxCQAAQJm4BAAAoExcAgAAUCYuAQAAKBOXAAAAlIlLAAAAysQlAAAAZeISAACAMnEJAABAmbgEAACgTFwCAABQJi4BAAAoE5cAAACUiUsAAADKxCUAAABl4hIAAIAycQkAAECZuAQAAKBMXAIAAFAmLgEAACgTlwAAAJSJSwAAAMrEJQAAAGXiEgAAgDJxCQAAQJm4BAAAoExcAgAAUCYuAQAAKBOXAAAAlIlLAAAAysQlAAAAZeISAACAMnEJAABAmbgEAACgTFwCAABQJi4BAAAoE5cAAACUiUsAAADKxCUAAABl4hIAAIAycQkAAECZuAQAAKBMXAIAAFAmLgEAACgTlwAAAJSJSwAAAMrEJQAAAGXiEgAAgDJxCQAAQJm4BAAAoExcAgAAUCYuAQAAKBOXAAAAlIlLAAAAysQlAAAAZeISAACAMnEJAABAmbgEAACgTFwCAABQJi4BAAAoE5cAAACUiUsAAADKxCUAAABl4hIAAIAycQkAAECZuAQAAKBMXAIAAFAmLgEAACgTlwAAAJSJSwAAAMrEJQAAAGXiEgAAgDJxCQAAQJm4BAAAoExcAgAAUCYuAQAAKBOXAAAAlIlLAAAAysQlAAAAZeISAACAMnEJAABAmbgEAACgTFwCAABQJi4BAAAoE5cAAACUiUsAAADKxCUAAABl4hIAAIAycQkAAECZuAQAAKBMXAIAAFAmLgEAACgTlwAAAJSJSwAAAMrEJQAAAGXiEgAAgDJxCQAAQJm4BAAAoExcAgAAUCYuAQAAKBOXAAAAlIlLAAAAysQlAAAAZeISAACAMnEJAABAmbgEAACgTFwCAABQJi4BAAAoE5cAAACUiUsAAADKxCUAAABl4hIAAIAycQkAAECZuAQAAKBMXAIAAFAmLgEAACgTlwAAAJSJSwAAAMrEJQAAAGXiEgAAgDJxCQAAQJm4BAAAoExcAgAAUCYuAQAAKBOXAAAAlIlLAAAAysQlAAAAZeISAACAMnEJAABAmbgEAACgTFwCAABQJi4BAAAoE5cAAACUiUsAAADKxCUAAABl4hIAAIAycQkAAECZuAQAAKBMXAIAAFAmLgEAACgTlwAAAJSJSwAAAMrEJQAAAGXiEgAAgDJxCQAAQJm4BAAAoExcAgAAUCYuAQAAKBOXAAAAlIlLAAAAysQlAAAAZeISAACAMnEJAABAmbgEAACgTFwCAABQJi4BAAAoE5cAAACUiUsAAADKxCUAAABl4hIAAIAycQkAAECZuAQAAKBMXAIAAFAmLgEAACgTlwAAAJSJSwAAAMrEJQAAAGXiEgAAgDJxCQAAQJm4BAAAoExcAgAAUCYuAQAAKBOXAAAAlIlLAAAAysQlAAAAZeISAACAMnEJAABAmbgEAACgTFwCAABQJi4BAAAoE5cAAACUiUsAAADKxCUAAABl4hIAAIAycQkAAECZuAQAAKBMXAIAAFAmLgEAACgTlwAAAJSJSwAAAMrEJQAAAGXiEgAAgDJxCQAAQJm4BAAAoExcAgAAUCYuAQAAKBOXAAAAlIlLAACAE9Crr76aPXv2HLP3E5cAAAAnoLvvvjunn356vvjFL2bLli3l9xOXAAAAJ6Bbbrklt956a+69995MmjSpHJmDurq6uhp5werVq/PrX/+6oQ+ZO3duNm3alK9+9asNvQ4AAIC+1dHRkR//+Md55JFH0tXVlauuuirf+973MmrUqIbep+G4vPLKK7Ns2bKGPgQAAIDjx7x583LzzTc39JqG43L//v3Zt29fQx8yY8aMPPPMM3n00Ucbeh0AAAB9p7OzMw888EC+9a1vZf369Zk+fXq+8pWv5OKLL274vYY2+oLBgwdn+PDhDb1myJAhb+h1AAAAJ6r29vbcfvvtWbx4cfbu3ZuLLrood999dyZMmHBM3n/VqlW56aabsm7durS1tWXBggVpbW19w+9nQx8AAIABZunSpZk6dWrGjx+ftWvX5tlnn82OHTtyxRVXZNeuXcfkM9rb2zNt2rQ88cQT+eEPf1gKy0RcAgAADCiPP/54pk+fnlmzZmX27NkZMWJERowYkdtuuy2/+c1vcu+99x6Tz2lra8v8+fPLUXmAuAQAABggNm/enBtuuCEtLS2H3G3j3HPPTZLcf//9zRjtiBr+ziUAAADHXldXV2666aa88MIL+f73v5+hQw/OtdNOOy1J8uSTT2b37t05+eSTmzHmYVm5BAAAGADuu+++LFu2LGeccUauu+66Q55/+eWXk3TfwWPTpk39Pd4RiUsAAIAm27NnT89lsLfcckuGDRt2yDlbtmzp9XigEJcAAABNNm/evLS3tydJPvOZz/R6zpo1a3qOBw0a1C9zNUJcAgAANNk999yTJLngggty9tln93rOk08+2XM8ZsyYfpmrEeISAACgiVatWpVf/epXSZJrr732sOetWLGi5/iMM87o67EaJi4BAACa6MEHH+w5Plxc7tmzJytXrkyStLS09OwcO5C4FQkAAEATLVq0KEkydOjQXHPNNb2es3v37vz+979PkrS2tvbbbI2wcgkAANAk27Zt67kk9pOf/GSee+65Xv9cf/31Pa+59NJLe32vDRs2ZMKECf0yd2/EJQAAQJO8/nuU73vf+w573vLly3uOP/KRjxz03EsvvZS5c+dm2rRp2bhx47Ef8ii5LBYAAKBJ1q5d23N80UUX9XrO888/n6eeeipJcv7552fixIk9z11zzTXp6OjIJZdc0qdzHg1xCQAA0CQHVhpHjBiRd73rXb2es3Dhwp7jP74H5o9+9KOe4wULFmTr1q3Hfsij1C+XxZ566qk59dRT++OjAAAAjhs7d+5MkpxzzjkZMmRIr+d897vfTZKMHj06H//4x/tttkb1y8rl7bffnn379vXHRwEAABx3Wlpaen386aef7rkFycyZM3PKKaf051gN6Ze4tGoJAABwqNGjRyc5fDMdWLUcN25cvvCFL/TXWG+I3WIBAACaZPLkyUmSwYMPTbOOjo585zvfSZLMmTMnp512Wr/O1ihxCQAA0CQHdnndtm3bIc/dcccd2bZtW2688cbceOON/T1awwZ1dXV1NXsIAACAE9V73/ve7Nix46B7VK5atSqXXHJJpk6dmuXLl2fkyJFHfJ+JEydmw4YNaVbiWbkEAABoorvuuivt7e2ZP39+kmTFihVpa2vLpZdemiVLlhxVWA4EVi4BAACabP78+Zk9e3b27t2bMWPGZObMmfn0pz/d63cxD6fZK5fiEgAA4E9As+PSZbEAAACUiUsAAADKxCUAAMBxbufOneno6EiSg3ad7U/iEgAA4Dg1Y8aMnHfeeRk/fnzPvTInT56c97znPbn66qv7dRYb+gAAAFBm5RIAAIAycQkAAECZuAQAAKBMXAIAAFAmLgEAACgTlwAAAJSJSwAAAMrEJQAAAGXiEgAAgDJxCQAAQJm4BAAAoExcAgAAUPZ/I0bskazEDVgAAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        }
      ],
      "source": [
        "import matplotlib\n",
        "from os import linesep\n",
        "from numpy import linspace,cos,pi,sin,random,arccos,heaviside,array,fft,abs, angle,sqrt,c_\n",
        "from scipy.optimize import fsolve\n",
        "import matplotlib.pyplot as plt\n",
        "from matplotlib import animation,rc\n",
        "rc('animation', html='jshtml')\n",
        "plt.rcParams.update({\"font.family\" : \"serif\",\"mathtext.fontset\" : \"cm\"})\n",
        "\n",
        "#functions\n",
        "def h(theta,t): #square waveform with one firing angle theta\n",
        "  return (2*heaviside(sin(2*pi*t),0)-1)*heaviside(abs(sin(2*pi*t))-sin(theta),0)\n",
        "\n",
        "#spectrum\n",
        "def spectrum(xt): #spectrum - magnitude and angle\n",
        "  xfft = fft.rfft(xt)\n",
        "  xfft_angle = angle(xfft)\n",
        "  xfft_mag = abs(xfft)/len(xt)*2\n",
        "  xfft_mag[0]=xfft_mag[0]*0.5 #adjust the dc component\n",
        "  return xfft_mag,xfft_angle\n",
        "\n",
        "# draw dotted line between two sets of points in two different plots\n",
        "def l2Ax(fig,ax1,ax2,ps1,ps2):\n",
        "  line_set = []\n",
        "  transFigure = fig.transFigure.inverted()\n",
        "  for k in range(len(ps1[0])):\n",
        "    p1 = transFigure.transform(ax1.transData.transform([ps1[0][k],ps1[1][k]]))\n",
        "    p2 = transFigure.transform(ax2.transData.transform([ps2[0][k],ps2[1][k]]))\n",
        "    line = matplotlib.lines.Line2D((p1[0],p2[0]),(p1[1],p2[1]),\n",
        "                      transform=fig.transFigure,ls='--',color='gray', lw=1)\n",
        "    line_set.append(line)\n",
        "  fig.lines = line_set\n",
        "  return\n",
        "\n",
        "#parameters\n",
        "Nf = 200 # number of frames\n",
        "Np = 4096 # number of points\n",
        "theta1_v3 = linspace(0,pi/2 - pi/3,Nf) #firing angle 1 to eliminate 3rd-order harmonic\n",
        "theta1_v5 = linspace(0,pi/2 - pi/5,Nf) #firing angle 1 to eliminate 5th-order harmonic\n",
        "theta1_v7 = linspace(0,pi/2 - pi/7,Nf) #firing angle 1 to eliminate 7th-order harmonic\n",
        "theta2_v3 = theta1_v3 + pi/3 #firing angle 2 to eliminate 3rd-order harmonic\n",
        "theta2_v5 = theta1_v5 + pi/5 #firing angle 2 to eliminate 5th-order harmonic\n",
        "theta2_v7 = theta1_v7 + pi/7 #firing angle 2 to eliminate 7th-order harmonic\n",
        "\n",
        "#time series\n",
        "t = linspace(0,1,Np)\n",
        "\n",
        "#spectrum and frequency series\n",
        "fser = linspace(0,1,33) #frequency series for plotting\n",
        "fser_sliced = fser[1:13:2]\n",
        "\n",
        "ht3,ht5,ht7 = ([],[],[])\n",
        "hfft3, hfft5,hfft7 = ([],[],[])\n",
        "for k in range(Nf):\n",
        "  htk = (h(theta1_v3[k],t) + h(theta2_v3[k],t))/2.2\n",
        "  ht3.append(htk)\n",
        "  hfft3.append(spectrum(htk))\n",
        "  htk = (h(theta1_v5[k],t) + h(theta2_v5[k],t))/2.2\n",
        "  ht5.append(htk)\n",
        "  hfft5.append(spectrum(htk))\n",
        "  htk = (h(theta1_v7[k],t) + h(theta2_v7[k],t))/2.2\n",
        "  ht7.append(htk)\n",
        "  hfft7.append(spectrum(htk))\n",
        "\n",
        "#set up figures\n",
        "fig = plt.figure(figsize=(9,9))\n",
        "fig.tight_layout()\n",
        "ax_frame = [[[0, .7, 1, .275], [-.02,1.1], [-1.1, 1.1]], #sine1 & triangle\n",
        "            [[0, .4, 1, .275], [-.02,1.1], [-.1, 1.4]], #sine2 & triangle\n",
        "            [[0, 0.02,  1, .35], [-.02,1.1], [-.1,1.4]], # ac output spectrum\n",
        "            ]# [pos-boundary, xlim, ylim] for subplots\n",
        "ax_lbl=[[['$h(t)$'], r'$\\omega t$'],\n",
        "        [['$|H_{fft}|$'], '$f$'],\n",
        "        [[r'$H_{1}$','$H_{3}$','$H_{5}$','$H_{7}$','$H_{9}$','$H_{11}$'], r'$\\theta_1$'],\n",
        "        ] #yx axes labels\n",
        "title_lst = ['Square Waveform',\n",
        "             'Spectrum',\n",
        "             'Magnitudes',\n",
        "             ]\n",
        "clst=[['b'], ['brown'], ['k','r','m','b','orange','c']] #colors of lines\n",
        "ax_lst=[] #axis list or the list of subplots\n",
        "tx_lst=[]\n",
        "lines = [] #array for all lines\n",
        "for k in range(len(ax_frame)):\n",
        "  xn,xm,yn,ym = ax_frame[k][1]+ax_frame[k][2]\n",
        "  ax=fig.add_axes(ax_frame[k][0],xlim=[xn,xm], ylim=[yn,ym],fc='none') #no fc\n",
        "  ax.axis('off') #turn off axis frames\n",
        "  ax_lst.append(ax)\n",
        "  ax.annotate(\"\", (xm, 0), (xn, 0),arrowprops={'arrowstyle':\"->\"}) #x-axis\n",
        "  ax.annotate(\"\", (0,ym), (0,yn),arrowprops={'arrowstyle':\"->\"}) #y-axis\n",
        "  for i in range(len(ax_lbl[k][0])):\n",
        "    y_text=ax.text(0.01+0.06*i,ym,ax_lbl[k][0][i],size=24,va='center',color=clst[k][i]) #y label\n",
        "    if k==1:\n",
        "      lines.append(ax.plot([], [], clst[k][i], ls='',marker='x',markersize=6)[0])\n",
        "      lines.append(ax.scatter(fser_sliced, fser_sliced*0,s=200,ec=clst[-1],\n",
        "                              marker='o',fc='none',lw=1.5)) #circle H1,3,5,7,9,11\n",
        "    else:\n",
        "      lines.append(ax.plot([], [], clst[k][i], lw=2)[0])\n",
        "  ax.text(xm-.02,0,ax_lbl[k][-1],size=24,va='top',ha='right') #x-axis label\n",
        "  tx_lst.append(ax.text(0.55,ax_frame[k][2][1],title_lst[k],color='k',\n",
        "                        va='center',ha='center',size=18))\n",
        "tx_firingAngle1 = ax_lst[0].text(0,0,r'$\\theta_1$',size=24,va='top',ha='center') #theta1 label\n",
        "tx_firingAngle2 = ax_lst[0].text(0.2,0,r'$\\theta_2$',size=24,va='top',ha='center') #theta2 label\n",
        "\n",
        "# animation function. This is called sequentially\n",
        "def animate(i):\n",
        "  k=int(i/Nf)\n",
        "  i=i%Nf\n",
        "  if k==0:\n",
        "    tx_lst[0].set_text('Pattern with $3^{rd}$-Order Eliminated')\n",
        "    hti = ht3[i]\n",
        "    hffti = hfft3[i]\n",
        "    lines[2].set_linewidth([1,3,1,1,1,1])\n",
        "    tx_firingAngle1.set_position([theta1_v3[i]/(2*pi),0])\n",
        "    if theta2_v3[i] == pi/2:\n",
        "      tx_firingAngle2.set_text('')\n",
        "      l2Ax(fig,ax_lst[0],ax_lst[0],[[theta2_v3[i]/(2*pi)],[0]],[[theta2_v3[i]/(2*pi)],[0]])\n",
        "    else:\n",
        "      tx_firingAngle2.set_text(r'$\\theta_2$')\n",
        "      tx_firingAngle2.set_position([theta2_v3[i]/(2*pi),0])\n",
        "      l2Ax(fig,ax_lst[0],ax_lst[0],[[theta2_v3[i]/(2*pi)],[0.5]],[[theta2_v3[i]/(2*pi)],[0]])\n",
        "    for kk in range(6):\n",
        "      lines[3+kk].set_data(theta1_v3[0:i]/(pi/2-pi/3),array(hfft3)[0:i,0,2*kk+1]) #magnitudes\n",
        "  elif k==1:\n",
        "    tx_lst[0].set_text('Pattern with $5^{th}$-Order Eliminated')\n",
        "    hti = ht5[i]\n",
        "    hffti = hfft5[i]\n",
        "    lines[2].set_linewidth([1,1,3,1,1,1])\n",
        "    tx_firingAngle1.set_position([theta1_v5[i]/(2*pi),0])\n",
        "    if theta2_v5[i] == pi/2:\n",
        "      tx_firingAngle2.set_text('')\n",
        "      l2Ax(fig,ax_lst[0],ax_lst[0],[[theta2_v5[i]/(2*pi)],[0]],[[theta2_v5[i]/(2*pi)],[0]])\n",
        "    else:\n",
        "      tx_firingAngle2.set_text(r'$\\theta_2$')\n",
        "      tx_firingAngle2.set_position([theta2_v5[i]/(2*pi),0])\n",
        "      l2Ax(fig,ax_lst[0],ax_lst[0],[[theta2_v5[i]/(2*pi)],[0.5]],[[theta2_v5[i]/(2*pi)],[0]])\n",
        "    for kk in range(6):\n",
        "      lines[3+kk].set_data(theta1_v5[0:i]/(pi/2-pi/5),array(hfft5)[0:i,0,2*kk+1]) #magnitudes\n",
        "  else:\n",
        "    tx_lst[0].set_text('Pattern with $7^{th}$-Order Eliminated')\n",
        "    hti = ht7[i]\n",
        "    hffti = hfft7[i]\n",
        "    lines[2].set_linewidth([1,1,1,3,1,1])\n",
        "    tx_firingAngle1.set_position([theta1_v7[i]/(2*pi),0])\n",
        "    if theta2_v7[i] == pi/2:\n",
        "      tx_firingAngle2.set_text('')\n",
        "      l2Ax(fig,ax_lst[0],ax_lst[0],[[theta2_v7[i]/(2*pi)],[0]],[[theta2_v7[i]/(2*pi)],[0]])\n",
        "    else:\n",
        "      tx_firingAngle2.set_text(r'$\\theta_2$')\n",
        "      tx_firingAngle2.set_position([theta2_v7[i]/(2*pi),0])\n",
        "      l2Ax(fig,ax_lst[0],ax_lst[0],[[theta2_v7[i]/(2*pi)],[0.5]],[[theta2_v7[i]/(2*pi)],[0]])\n",
        "    for kk in range(6):\n",
        "      lines[3+kk].set_data(theta1_v7[0:i]/(pi/2-pi/7),array(hfft7)[0:i,0,2*kk+1]) #magnitudes\n",
        "\n",
        "  hfft_sliced = hffti[0][1:13:2]\n",
        "  lines[0].set_data(t,hti) #ac output\n",
        "  lines[1].set_data(fser,hffti[0][:len(fser)]) #spectrum\n",
        "  lines[2].set_offsets(c_[fser_sliced,hfft_sliced]) # marked components\n",
        "  return\n",
        "\n",
        "anim = animation.FuncAnimation(fig, animate, frames=Nf*3, interval=50)\n",
        "# anim #uncomment to generate animation in the output area\n",
        "# to save the animation, uncomment the following three lines\n",
        "fn = r\"OptimalSwitching_SHE_Multilevel.mp4\"\n",
        "writervideo = animation.FFMpegWriter(fps=25)\n",
        "anim.save(fn, writer=writervideo,dpi = 120)"
      ]
    }
  ]
}