test_jacobian.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Verify the Jacobian Approximation of Measurement Equation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 277,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pypower.api import case14,ext2int, makeYbus, runopf\n",
    "from pypower.idx_brch import *\n",
    "from pypower.idx_bus import *\n",
    "import scipy\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 332,
   "metadata": {},
   "outputs": [],
   "source": [
    "case = case14()\n",
    "# No transformer\n",
    "case['branch'][:,TAP] = 0\n",
    "case['branch'][:,SHIFT] = 0\n",
    "#case['branch'][:,BR_B] = 0\n",
    "case['bus'][:,GS] = 0\n",
    "case['bus'][:,BS] = 0\n",
    "case_int = ext2int(case)\n",
    "ys = 1/(case['branch'][:,BR_R]+1j*case['branch'][:,BR_X])\n",
    "gs = np.real(ys)\n",
    "bs = np.imag(ys)\n",
    "gsh = case['bus'][:,GS]/case['baseMVA']\n",
    "bsh = case['bus'][:,BS]/case['baseMVA']\n",
    "bc = case_int['branch'][:,BR_B]\n",
    "\n",
    "Ybus, Yf, Yt = makeYbus(case_int['baseMVA'], case_int['bus'], case_int['branch'])\n",
    "# To numpy\n",
    "Ybus = scipy.sparse.csr_matrix.todense(Ybus).getA()\n",
    "Yf = scipy.sparse.csr_matrix.todense(Yf).getA()\n",
    "Yt = scipy.sparse.csr_matrix.todense(Yt).getA()\n",
    "# Branch Incidence Matrix\n",
    "f_bus = case_int['branch'][:, 0].astype('int')        # list of \"from\" buses\n",
    "t_bus = case_int['branch'][:, 1].astype('int')        # list of \"to\" buses\n",
    "\n",
    "no_bus = len(case['bus'])\n",
    "no_brh = len(case['branch'])\n",
    "\n",
    "Cf = np.zeros((no_brh,no_bus))         # \"from\" bus incidence matrix\n",
    "Ct = np.zeros((no_brh,no_bus))         # \"to\" bus incidence matrix\n",
    "for i in range(no_brh):\n",
    "    Cf[i,f_bus[i]] = 1\n",
    "    Ct[i,t_bus[i]] = 1\n",
    "A = Cf-Ct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 333,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PYPOWER Version 5.1.4, 27-June-2018 -- AC Optimal Power Flow\n",
      "Python Interior Point Solver - PIPS, Version 1.0, 07-Feb-2011\n",
      "Converged!\n",
      "\n",
      "Converged in 0.76 seconds\n",
      "Objective Function Value = 8082.58 $/hr\n",
      "================================================================================\n",
      "|     System Summary                                                           |\n",
      "================================================================================\n",
      "\n",
      "How many?                How much?              P (MW)            Q (MVAr)\n",
      "---------------------    -------------------  -------------  -----------------\n",
      "Buses             14     Total Gen Capacity     772.4         -52.0 to 148.0\n",
      "Generators         5     On-line Capacity       772.4         -52.0 to 148.0\n",
      "Committed Gens     5     Generation (actual)    268.3              89.8\n",
      "Loads             11     Load                   259.0              73.5\n",
      "  Fixed           11       Fixed                259.0              73.5\n",
      "  Dispatchable     0       Dispatchable           0.0 of 0.0        0.0\n",
      "Shunts             0     Shunt (inj)              0.0               0.0\n",
      "Branches          20     Losses (I^2 * Z)         9.30             40.67\n",
      "Transformers       0     Branch Charging (inj)     -               24.4\n",
      "Inter-ties         0     Total Inter-tie Flow     0.0               0.0\n",
      "Areas              1\n",
      "\n",
      "                          Minimum                      Maximum\n",
      "                 -------------------------  --------------------------------\n",
      "Voltage Magnitude   0.977 p.u. @ bus 14         1.060 p.u. @ bus 1   \n",
      "Voltage Angle     -14.85 deg   @ bus 14         0.00 deg   @ bus 1   \n",
      "P Losses (I^2*R)             -                  2.89 MW    @ line 1-2\n",
      "Q Losses (I^2*X)             -                  8.82 MVAr  @ line 1-2\n",
      "Lambda P           36.72 $/MWh @ bus 1         41.34 $/MWh @ bus 14  \n",
      "Lambda Q           -0.17 $/MWh @ bus 1          0.73 $/MWh @ bus 14  \n",
      "\n",
      "================================================================================\n",
      "|     Bus Data                                                                 |\n",
      "================================================================================\n",
      " Bus      Voltage          Generation             Load          Lambda($/MVA-hr)\n",
      "  #   Mag(pu) Ang(deg)   P (MW)   Q (MVAr)   P (MW)   Q (MVAr)     P        Q   \n",
      "----- ------- --------  --------  --------  --------  --------  -------  -------\n",
      "    1  1.060    0.000*   194.24      0.00       -         -      36.716  -0.171\n",
      "    2  1.040   -4.000     36.72     16.92     21.70     12.70    38.358     -\n",
      "    3  1.016   -9.922     28.83     24.89     94.20     19.00    40.577     -\n",
      "    4  1.015   -8.671       -         -       47.80     -3.90    40.202   0.141\n",
      "    5  1.020   -7.473       -         -        7.60      1.60    39.666   0.148\n",
      "    6  1.016  -13.242      0.00     24.00     11.20      7.50    39.768   0.047\n",
      "    7  1.019  -11.391       -         -         -         -      40.177   0.260\n",
      "    8  1.059  -10.594      8.52     24.00       -         -      40.170   0.241\n",
      "    9  0.997  -13.376       -         -       29.50     16.60    40.175   0.337\n",
      "   10  0.993  -13.669       -         -        9.00      5.80    40.352   0.451\n",
      "   11  1.001  -13.594       -         -        3.50      1.80    40.200   0.332\n",
      "   12  1.000  -14.149       -         -        6.10      1.60    40.474   0.284\n",
      "   13  0.995  -14.187       -         -       13.50      5.80    40.682   0.449\n",
      "   14  0.977  -14.849       -         -       14.90      5.00    41.340   0.728\n",
      "                        --------  --------  --------  --------\n",
      "               Total:    268.30     89.81    259.00     73.50\n",
      "\n",
      "================================================================================\n",
      "|     Branch Data                                                              |\n",
      "================================================================================\n",
      "Brnch   From   To    From Bus Injection   To Bus Injection     Loss (I^2 * Z)  \n",
      "  #     Bus    Bus    P (MW)   Q (MVAr)   P (MW)   Q (MVAr)   P (MW)   Q (MVAr)\n",
      "-----  -----  -----  --------  --------  --------  --------  --------  --------\n",
      "   0      1      2    129.41     -4.54   -126.52      7.54     2.889      8.82\n",
      "   1      1      5     64.83      4.54    -62.78     -1.41     2.047      8.45\n",
      "   2      2      3     55.54     -0.44    -54.19      1.46     1.342      5.65\n",
      "   3      2      4     48.84     -1.61    -47.56      1.91     1.282      3.89\n",
      "   4      2      5     37.16     -1.27    -36.43     -0.18     0.728      2.22\n",
      "   5      3      4    -11.18      4.43     11.28     -5.50     0.098      0.25\n",
      "   6      4      5    -49.94      4.74     50.26     -3.71     0.326      1.03\n",
      "   7      4      7     23.48     -1.18    -23.48      2.30    -0.000      1.12\n",
      "   8      4      9     14.94      3.93    -14.94     -2.64    -0.000      1.29\n",
      "   9      5      6     41.35      3.71    -41.35      0.47    -0.000      4.17\n",
      "  10      6     11      5.62      5.19     -5.57     -5.08     0.054      0.11\n",
      "  11      6     12      7.62      2.76     -7.54     -2.60     0.078      0.16\n",
      "  12      6     13     16.91      8.08    -16.69     -7.64     0.225      0.44\n",
      "  13      7      8     -8.52    -22.98      8.52     24.00    -0.000      1.02\n",
      "  14      7      9     32.00     20.68    -32.00    -19.14    -0.000      1.54\n",
      "  15      9     10      6.96      2.60     -6.94     -2.55     0.018      0.05\n",
      "  16      9     14     10.48      2.59    -10.33     -2.27     0.149      0.32\n",
      "  17     10     11     -2.06     -3.25      2.07      3.28     0.012      0.03\n",
      "  18     12     13      1.44      1.00     -1.43     -1.00     0.007      0.01\n",
      "  19     13     14      4.62      2.83     -4.57     -2.73     0.051      0.10\n",
      "                                                             --------  --------\n",
      "                                                    Total:     9.305     40.67\n",
      "\n",
      "================================================================================\n",
      "|     Generation Constraints                                                   |\n",
      "================================================================================\n",
      " Gen   Bus                Active Power Limits\n",
      "  #     #    Pmin mu    Pmin       Pg       Pmax    Pmax mu\n",
      "----  -----  -------  --------  --------  --------  -------\n",
      "   3     6    0.232      0.00      0.00    100.00      -  \n",
      "\n",
      "Gen  Bus              Reactive Power Limits\n",
      " #    #   Qmin mu    Qmin       Qg       Qmax    Qmax mu\n",
      "---  ---  -------  --------  --------  --------  -------\n",
      "  0    1   0.171      0.00      0.00     10.00      -  \n",
      "  3    6     -       -6.00     24.00     24.00    0.047\n",
      "  4    8     -       -6.00     24.00     24.00    0.241\n",
      "[ 1.29408306  0.648307    0.55536726  0.48840338  0.37158018 -0.1117923\n",
      " -0.49935361  0.23480407  0.14935962  0.41353356  0.05624287  0.07618532\n",
      "  0.16910727 -0.08523486  0.32003893  0.06959517  0.10480338 -0.02058136\n",
      "  0.01440359  0.0461929 ]\n"
     ]
    }
   ],
   "source": [
    "result = runopf(case)\n",
    "pf_result = result['branch'][:,PF]/100\n",
    "print(pf_result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 345,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Active power\n",
    "# From complex power\n",
    "vmag = result['bus'][:,VM]\n",
    "vang = result['bus'][:,VA]*np.pi/180  # In radian\n",
    "vcomplex = vmag*np.exp(1j*vang)\n",
    "Sf = np.diag(Cf@vcomplex)@np.conjugate(Yf)@np.conjugate(vcomplex)\n",
    "Yf1 = np.diag(ys)@A+1j*np.diag(bc/2)@Cf\n",
    "Yf2 = np.diag(ys)@A\n",
    "Sf2 = np.diag(Cf@vcomplex)@np.conjugate(Yf2)@np.conjugate(vcomplex)\n",
    "\n",
    "Ybus1 = A.T@np.diag(ys)@A + 1j*Cf.T@np.diag(bc/2)@Cf + 1j*Ct.T@np.diag(bc/2)@Ct\n",
    "Ybus2 = A.T@np.diag(ys)@A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 346,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.+0.051j , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ],\n",
       "       [0.+0.j    , 0.+0.0826j, 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ],\n",
       "       [0.+0.j    , 0.+0.j    , 0.+0.0283j, 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ],\n",
       "       [0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.0234j, 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ],\n",
       "       [0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.0419j,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ],\n",
       "       [0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ],\n",
       "       [0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ],\n",
       "       [0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ],\n",
       "       [0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ],\n",
       "       [0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ],\n",
       "       [0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ],\n",
       "       [0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ],\n",
       "       [0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ],\n",
       "       [0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ,\n",
       "        0.+0.j    , 0.+0.j    , 0.+0.j    , 0.+0.j    ]])"
      ]
     },
     "execution_count": 346,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ybus1 - Ybus2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 337,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.02966304, -0.02764056, -0.02367559, -0.01837831, -0.01870263,\n",
       "       -0.00661157,  0.        ,  0.        ,  0.        ,  0.        ,\n",
       "        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,\n",
       "        0.        ,  0.        ,  0.        ,  0.        ,  0.        ])"
      ]
     },
     "execution_count": 337,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.imag(Sf) - np.imag(Sf2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 340,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "19.4980702"
      ]
     },
     "execution_count": 340,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "15.26308652 + 4.23498368"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 339,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-19.44707021,  15.26308652,   0.        ,   0.        ,\n",
       "          4.23498368,   0.        ,   0.        ,   0.        ,\n",
       "          0.        ,   0.        ,   0.        ,   0.        ,\n",
       "          0.        ,   0.        ],\n",
       "       [ 15.26308652, -30.2721154 ,   4.78186315,   5.11583833,\n",
       "          5.1939274 ,   0.        ,   0.        ,   0.        ,\n",
       "          0.        ,   0.        ,   0.        ,   0.        ,\n",
       "          0.        ,   0.        ],\n",
       "       [  0.        ,   4.78186315,  -9.82238013,   5.06881698,\n",
       "          0.        ,   0.        ,   0.        ,   0.        ,\n",
       "          0.        ,   0.        ,   0.        ,   0.        ,\n",
       "          0.        ,   0.        ],\n",
       "       [  0.        ,   5.11583833,   5.06881698, -38.31973174,\n",
       "         21.57855398,   0.        ,   4.78194338,   0.        ,\n",
       "          1.79797907,   0.        ,   0.        ,   0.        ,\n",
       "          0.        ,   0.        ],\n",
       "       [  4.23498368,   5.1939274 ,   0.        ,  21.57855398,\n",
       "        -34.93350411,   3.96793905,   0.        ,   0.        ,\n",
       "          0.        ,   0.        ,   0.        ,   0.        ,\n",
       "          0.        ,   0.        ],\n",
       "       [  0.        ,   0.        ,   0.        ,   0.        ,\n",
       "          3.96793905, -17.34073281,   0.        ,   0.        ,\n",
       "          0.        ,   0.        ,   4.09407434,   3.17596397,\n",
       "          6.10275545,   0.        ],\n",
       "       [  0.        ,   0.        ,   0.        ,   4.78194338,\n",
       "          0.        ,   0.        , -19.54900595,   5.67697985,\n",
       "          9.09008272,   0.        ,   0.        ,   0.        ,\n",
       "          0.        ,   0.        ],\n",
       "       [  0.        ,   0.        ,   0.        ,   0.        ,\n",
       "          0.        ,   0.        ,   5.67697985,  -5.67697985,\n",
       "          0.        ,   0.        ,   0.        ,   0.        ,\n",
       "          0.        ,   0.        ],\n",
       "       [  0.        ,   0.        ,   0.        ,   1.79797907,\n",
       "          0.        ,   0.        ,   9.09008272,   0.        ,\n",
       "        -24.28250638,  10.36539413,   0.        ,   0.        ,\n",
       "          0.        ,   3.02905046],\n",
       "       [  0.        ,   0.        ,   0.        ,   0.        ,\n",
       "          0.        ,   0.        ,   0.        ,   0.        ,\n",
       "         10.36539413, -14.76833788,   4.40294375,   0.        ,\n",
       "          0.        ,   0.        ],\n",
       "       [  0.        ,   0.        ,   0.        ,   0.        ,\n",
       "          0.        ,   4.09407434,   0.        ,   0.        ,\n",
       "          0.        ,   4.40294375,  -8.49701809,   0.        ,\n",
       "          0.        ,   0.        ],\n",
       "       [  0.        ,   0.        ,   0.        ,   0.        ,\n",
       "          0.        ,   3.17596397,   0.        ,   0.        ,\n",
       "          0.        ,   0.        ,   0.        ,  -5.42793859,\n",
       "          2.25197463,   0.        ],\n",
       "       [  0.        ,   0.        ,   0.        ,   0.        ,\n",
       "          0.        ,   6.10275545,   0.        ,   0.        ,\n",
       "          0.        ,   0.        ,   0.        ,   2.25197463,\n",
       "        -10.66969355,   2.31496348],\n",
       "       [  0.        ,   0.        ,   0.        ,   0.        ,\n",
       "          0.        ,   0.        ,   0.        ,   0.        ,\n",
       "          3.02905046,   0.        ,   0.        ,   0.        ,\n",
       "          2.31496348,  -5.34401393]])"
      ]
     },
     "execution_count": 339,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.imag(Ybus)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 351,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f9368b59b40>]"
      ]
     },
     "execution_count": 351,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAzOklEQVR4nO3deXxjdbnH8c+TtGnatE33dTobzL4xM50VGUSQRRTEizqAsogiKle9ilfUK1cRFcQFr6KAsowKw66igmyioMy+Mxuzb+10uqZ7muV3/0hm7HTSmXZykrTJ8369+po05+ScZ9L02+Sc3+85YoxBKaVU8rMlugCllFLxoYGvlFIpQgNfKaVShAa+UkqlCA18pZRKEWmJLqA/RUVFZvTo0YkuQymlhpU1a9Y0GGOKIy0bsoE/evRoVq9enegylFJqWBGRff0t00M6SimVIjTwlVIqRWjgK6VUitDAV0qpFKGBr5RSKUIDXymlUoQGvlJKpQgNfKWUpbYs/yu7Ni1PdBkqAksCX0QuFpHtIrJTRG6LsPxLIrJFRDaKyGsiMsqK/Sqlhp6sl2+l+0+3JroMFUHUgS8iduA+4BJgMnCViEzus9o6oNoYMx14BvhBtPtVSg1N7mALI707CQYCiS5F9WHFO/y5wE5jzG5jTA/wBHB57xWMMa8bYzrD3y4HRliwX6XUEBPw+3GbdnKki5o9WxNdjurDisCvBA70+v5g+L7+3Ai8GGmBiNwkIqtFZHV9fb0FpSml4qml8TA2CV02te6dFQmuRvUV15O2IvIxoBq4J9JyY8yDxphqY0x1cXHEZm9KqSGsvanu2O2eA+sTV4iKyIpumYeAql7fjwjfdxwRuQD4BnCuMcZrwX6VUkNMe3Mo8INGyG56O8HVqL6seIe/ChgnImNExAEsBp7vvYKIzAQeAC4zxhyxYJ9KqSHI6wkF/s70cYzw7sAEgwmuSPUWdeAbY/zALcBLwFbgKWPMZhG5Q0QuC692D5ANPC0i60Xk+X42p5QaxnxtoXNvTWXnkE8bdQd3Jbgi1ZslF0AxxrwAvNDnvtt73b7Aiv0opYa2YEcjAO5pF8HBh6jdtoKykeMSXJU6SmfaKqUsY+tsoJUsxkw7G7+x0b1/TaJLUr1o4CulLJPW3USr5OLMyuaAfSRZjXridijRwFdKWSajp4kOex4ADTkTqex6J7EFqeNo4CulLJPp99DlyAcgUDadIlpoqOn3mtoqzjTwlVKWyQl46AkHvntsNQCHti1LZEmqFw18pZQlTDBInvEQyCwAoGrSXIJG6Ny3LsGVqaM08JVSlmhrbcYhAcRVBEB2bj4H7RU46zcluDJ1lAa+UsoSbU2HAbBn/7sP1hHXRMo79cTtUKGBr5SyxNHAd+QWHbvPXzqNMupprq9NVFmqFw18pZQlultCbbIy88qO3Zc9ejYAB7fqJQ+HAg18pZQlfK2hwM/O/3fgV01ZCED7ntUJqUkdTwNfKWWJQHsDAO6ifwe+O7+IQ1JKhp64HRI08JVS1uhsxGvSyXLlHnd3nWsCJR3bE1SU6k0DXyllCXtXIy2Si9iOjxVv8TRGmMN4mhsSVJk6SgNfKWUJR08TbeE+Or25RoVP3G7Ra9wmmga+UsoSmb4WOtPzTri/ctI8ANr26onbRNPAV0pZIrtXH53eCktHUEchaXUbE1CV6k0DXyllCXfQg99ZEHFZTdYEitv1xG2iaeArpaLW3dWBS7oxWYWRlxdNpSpwkM52T5wrU71p4CulouZpDLVVsLmKIy7PHDkLmxj2b1kZz7JUHxr4SqmotTWe2Eent4pJ8wFo2bUqbjWpE2ngK6Wi1tVSB0CGuzTi8uLyUTTixq4nbhNKA18pFTVvaz0Arl59dHoTm42DzvEUtm6NZ1mqDw18pVTU/G2hwHcXRg58gM7CqYwM7Ke7qyNeZak+NPCVUlEzHQ0EjJCbH/mkLYBz5EzSJMj+rToBK1E08JVSUbN1NeKRHGx2e7/rlE4Inbht1hO3CaOBr5SKWrq3mVab+6TrlI8chwcX1G6IU1WqLw18pVTUnD3NdKTlnXQdsdk4kDGOAj1xmzAa+EqpqLkCLXjTT+yj01d7/hRG+fbg6/HGoSrVlwa+UipquUEPvn766PSWVjUTh/jZv31dHKpSfWngK6WiEvD7cZt2TGbkPjq9lY6fC0DTDu2Nnwga+EqpqHia6rCJAVfktgq9VY6dSodxEqxZH/vC1Ak08JVSUTnaRyct59SBb7Pb2ec4E3eLnrhNBA18pVRU2pvDfXRyI/fR6as1fwqjfLsI+P2xLEtFoIGvlIqKt/UIAFn5Awt8e+VZZEoPB3foePx408BXSkXFF26clnuSPjq9FYdP3Nbv0N748WZJ4IvIxSKyXUR2ishtEZYvEpG1IuIXkSut2KdSamgIdjQAkFtQMqD1R5w5nS7jwH9ofQyrUpFEHfgiYgfuAy4BJgNXicjkPqvtB64HHo92f0qpocXW2UCbySTDmTWg9dPSHexPH0tO8+YYV6b6suId/lxgpzFmtzGmB3gCuLz3CsaYvcaYjUDQgv0ppYYQe3cznlP00emrJW8yI707CQYCMapKRWJF4FcCB3p9fzB8n1IqBWT0NNFhzxvUY6R8BjnSRc0eHZ4ZT0PqpK2I3CQiq0VkdX19faLLUUoNQJa/ha70vEE9puDMOQDUvaMzbuPJisA/BFT1+n5E+L5BM8Y8aIypNsZUFxf3fyEFpdTQkRPw0JNx6j46vY2cOJseY6fnwNoYVaUisSLwVwHjRGSMiDiAxcDzFmxXKTXEmWCQPNNKIHNwge/IcLI/bTTZTXriNp6iDnxjjB+4BXgJ2Ao8ZYzZLCJ3iMhlACIyR0QOAh8GHhAR/SkrlQTa21pwiB8ZQB+dvppyJ1Hl3YEJ6liOeEmzYiPGmBeAF/rcd3uv26sIHepRSiWR1sZacgDbaQS+KZ9BXvOfqT2wg/JRE6wvTp1gSJ20VUoNL21NocZpGe7Bn3PLPzM047Z2m864jRcNfKXUaev2hEbTZeYNrK1CbyMnVuM3Nrx64jZuNPCVSjKb33qB/e+sj8u+esKN07LzBx/4zqxsDthHktX4ttVlqX5o4CuVRPZsWcW4lz5Gw/PfjMv+gm1HG6cNrFNmXw05E6nsesfKktRJaOArlSR8PV78z96MQwLkdNfGZ6edjXhNOq7swbVWOCpQNp0iWmio2WdxYSoSDXylksSax7/NuMBOaqSEfH98Zqrbu5tokVzEdnpR4h5bDcChbcusLEv1QwNfqSSwZ8sqZu25n7XZ57Kv6gqKaMHb3Rnz/Tq8TbQNso9Ob1WT5hI0QudePXEbDxr4Sg1zfl8P/uc+Q4e4GPXxX2DPD015icdhkkxfM11pp3c4ByA7N5+D9gqcDZssrEr1RwNfqWFu1ePfYpx/B7vnfovC0hFkFo4CoOXwnpjv2xXw4B1kH52+jrgmUt6pJ27jQQNfqWFs79bVzN79AGuzFzH7fTcC4C4LBX5XQ+zf4buDrfid0QW+v3QaZTTQdOS0ei6qQdDAV2qY8vt68D17Mx2SxaiP//LY/cWVZ4SWNx3o76GW8HZ3ki1dmKzCqLaTPXo2AIe2aqvkWNPAV2qYWvX4t487lHNUpiuHFrKRtti+Y25pCA39PJ0+Or1VTVkIQPveNVHXpE5OA1+pYWjf1jXM3n0/a12LmHXxDScsb7SX4OyM7Vj89qY6ABy50V27wp1fxCEpJaNeT9zGmga+UsOM39eD99mb6ZBMRn78FxHHwLdllJLjrYtpHZ0toe1nuE9vlm1vda4JlHRsj3o76uQ08JUaZlYt/Q7j/e+we863KCqririON6ucwmBsJ195PaE+Oq7T6KNzwraKpzHCHMbT3BD1tlT/NPCVslAwEGDbipepO7grJtvft3UNs3f9grWuc5h1ySf6ryOnAjcddLS1xKQOAP/RPjoF0b/Dd40Knbg9uEVP3MaSBr5SFvA01bP8sTuouXMKE1/8MDm/WsCKp39k6dWcjj+U88uTtjNIzw+98284tNuy/fdlOhsJGMFdUBL1tionzQOgbc+qqLel+mfJFa9UajPBIG1tLXjqa2hvqqXbc4Sxsy7AfZodFIeTXRvfovH1+5jW9DLzpYdt6ZNZNekzZG5/jnmb72DDrhepuPbXFFeMjnpfq5d+h/n+d1g954dU93Mo56is4tBY/Na6vTBxVtT7jsTW1YhHciiw26PeVmHpCOooJO2InriNJQ38Plq7fdz5i4dYsPDdXDF/YqLLSRi/z0dzYy1tDbV0NNXi9dThb63DtNdj72rA4W0ky9dMrr+JPOMhV3zk9nr8hrVzmf6Vl067qdZQ1uPtZuPLS8je+CgTfVuoMA42FVxI4XmfY+L00BDDYOAWVjz9A6Zv/TE9Dy5k9ZxvU33pp057n/u2rWXmrl+yLvtdzL7kxlOun1c+FoCuhv2nvc9TSe9uotXmJrppV/9WkzWB4vZtFm1NRaKB38ebb77OD1q/yhsvzmDnmL9wZmlOokuKq462Fg789CLG+7ZTLIa+A+56jJ0WcdNqz6czvQCPawx7M4sguxh7dgkZeWV07HiDBYceZf1rT3DWe69OyP8jFuoO7mL3X3/O+IPPUo2Hg1LO8nFfZtIln2FuwfHPlM1uZ97ir3Fgx/vofPJTVK+6lTVb/8wZ1z9AXtHgTnL6fT14n/k0XZJB1cfvH9Af0eKK0QSNEGg5OKh9DYazp5nOtDzLttddNJWqfcvoaGvBlWPddtW/aeD3Ydb8BoBFsoH7l3yHkV+6C0da8r1L7c/G393GAv82lldcg+SNxJFbijOvDFdhOe7iSnLdhZTYbJzsqK1v4fvZe9drlLz1v3SffRnOrOy41W81Ewyy+a0/07PsQaa3/4tiDBuz5nFo3k1MPeeDjDjF4YyqcTPwf/WfLH/sW8zacz+tP5/HhkX3MOM9HxlwDaufuDN0KKf6nlMeyjkq3ZFBveRhj+HkK1fAQ2PmGMu2lzlyFrb9hgNbVjJx3oWWbVf9mwZ+L9sO1HFO12vsLruIbJuXG2oe5rHnF3HDh96f6NLiYufGZcw5/CQrCz/A/E//4rS3k+7IoOP8uxn9ytUsW/otFtz4QwurPDkTDLLisW9ha9xBMCMXk+FGMt3YM92kZeXhyC7AmVNAVm4+LncR2Tl52CKEdpuniS0vPkDZO48xNXiAZnJYVXE1oy78T84aM7hDfWnpDuZf/z12bXw/tj/czIw3PsXKTX9k8g0/Jzs3/6SP3bdtLTN3/iJ0KOd9nxzUfpvTisnsOjyoxwxGbtDD4Sj76PRWMWk+/BNadq8GDfyY0MDvZfNrjzNROrGdexM5I2fQdu9cFqy/jVXTZzPnzPJElxdTAb+fwPNfoFWymXjNj6Pe3pSzL2XN8vcwa/+jHNp9I5VjJ1lQ5amt+uPPmb/rpzSQh9N4yZauk64fNEKrZNEuLrps2XTbs/HbMzmzcwPzpJt30sazavp3mXbRDSzIdEVV2xnTF+Idv5xlS77KvJrfcvgnK9n33p8yZeH7Iq4f8Pvpfubm8KGck4/KiaQ9o4yirtgMDw34/bhNGyYzuj46vRWXj6IRN/a6jZZtUx1PAz/M6w9QtfdpGtLLKZr4HrDZSPvQ/Ux86iM8tfRWJnzlUXKd6YkuM2ZWPvMjFvi3s3b23cwqjH6YHUDV4h/jf3A+9U9/kcqvvmTJNk+mZs82Jq//HpszpjHpq//AZrfj9/XQ0dpMu6eJrrZGutua8HU04+9oIdDVAl0exOvB1tNKuq8Nh6+NHG8dW/LejfvczzJ+1rmW1pjhzGLBp3/GthUfIPuv/8mkl65m+YbFnHX9j3D2+YOyaukdzPdvZ/XsH1BdNnLQ++pxlVPUvgITDFp+8tzTVEeBGHBZF/his3HQOZ7C1q2WbVMdTwM/7J8rVnI+m9kz6UsUhX85MidfxJHJ1/ORLY/y4GMPc9ONn05wlbFRV7OPKVvvZbNzJjMvvcmy7ZZUjmH5mZ9m/q6fhk7gnr/Ysm33FfD7aXn8RnKBgmseOnaYJi3dgbuwdNBDRMfGoMbeJs67kI7JK1j16BeYX7eUvff8E99lv2TcWecAsG/7+tChHNfZzD7d0T3uEWQd8eJprrd8iGxbUx0FQFpOdH10+uosnMqUg0vo7uo44Q+gil7qnI08hY7ljxLAxqj3HH+ctOSKu2nIOoMr9n+Xl1Ym5xjh/Y9/ASc95H/4Z5a/E5z1ka+zzzaC4n/eTndXh6Xb7m3V0juY7HubrTP/h/JRE2K2Hyu5cvKY959L2Pjuh8kKdjD695ez7OGv4O3upPvp8Kicawc2KicSR0F48lWN9RdC6WgO99HJtfYPiXPkTNIkyP6tqy3drgrRwAcONrYyr/Wv7MtfiC2v8viF6U7cH1uCWzpxvvB5altif53QeFrzt2eZ0/4668fcSMWZ0yzfviPDSdt536PS1LFu6bct3z7A7rdXMGvnfax1nUP1ZZ+NyT5iafq7/4OMz69kg/s8Fux/kNa7JjPBv42ds2+n6DQO5RzlCk++ajuy16JK/63bEwr8rHxrA790wnwAmnel5oxbT2Mdyx78PMt/eXNMtq+HdIA1rz7F5dJC47sij4JIr5hG4znf5Nw3b+fhR77L9V+4E5tN4lyl9Tra2yh742scsFUy86rYhDHA1HMuZ+2Kc5m572Fq9nyCikGOcjkZb3cn8txNtEo2Y67/1bCd6OUuKKb6S8+y9sVHGLPim6ErWEUxUQugoCJ0YMrbaP3kK19rqI9OjgV9dHorHzkODy7Ktj7K6h8vAwQjAoS/RDCEf/ekz33H1gPXzP9g6rsus7S2WGpvbWbTs3cxde8S5tHNmrwLY3LuJeUDPxA0FGx/ghZbPoVn9T/8svA9n6dm20tcfeQBnnvp3Vx5yXvjWGVsrHvsf3gXdWy/8HHSMzJjuq+Kj/6Y4K8XUvf0f1Hx3y9att21S77CguBeNix6gBnFw38k1axLbqDnPVcxw54W9S97QckIfMZOsMX6K18FO0JdLa0+NyA2G1vKr6Dq8KuUt20Kx7dBjPn37fAXEPG203gxr/yVI2OmUVJp3TyBWOjqaGPD73/IhJ0PsYA21rkWknfpt5gzZV5M9pfygb9y42YWBNawd8InyLOfZBSOCOXXPkzbT+YxZfmX2TbtNSaOsPaEVTzteHs1c2t+y7qCi5g5/9KY76+s6kyWjf0UC/b8nA2vP82M8z4c9Ta3LHuReTWPsaLwMua9J3YnhOPNkeG0ZDv2tDSOSAHp7TWWbK836WykzWSS48yyfNsLPn1fVI8/tHsrhUvO4cDSL1Jy6x8tqspa3u5O1v/x/xi79X7m08xGZzX1F97OTItHhfU1PD//WujIm4+QJkFGXnDqETiSUwaX/5xJso+3f3Mr3b5AHCq0XiAQxPuHL9AlmZxxzb1x2+/sxd/kgFRQ+MY38XZHdy6kzdNE3sufp8ZWytQbfmZRhcmnJb2ErG7rL4SS1t2Ex+a2fLtWqBw7iXWjPsHs9r+z6Y2hFfh+Xw+rnvspTXdNZ97W71OfXsmWi55g+m2vWT4EOJKUDvym9m7OavgTe7Nn4igZP6DH5M74AAfHXcOVPX/gmad+G+MKY2PZsz9lqv9t9s78b3KLKuK2X0eGk5Z3f5cRppa1T9wR1ba2PvI5SoP1tF/yc+27chIdmeXk+Y5Yvt2MniY67EMz8AFmXvW/HJQy8l6/Leo3F1YIBgKs/vOD1H5vBnM23k5bWj6bznuESV97k8kLLolbHSkd+Mv/9jyjpI6MudcP6nEjPvxDjjhH8953/pflm96JTXExcvjwQaZs/hHvZExh+gduifv+p537Ida6FnHWnoeo3Xd6l7Rb9/LvmNvyAitHXMfEucP/XEos+VwVFAUbCAas/TSa5ffQlX7ythCJ5Mx00bTou1SZmqjfXETDBIOse/l37PvuTKpXfwW/OFi38D7GfX0F0879UNwHGaRs4BtjcG76He3ionzBRwf3YEcWudcsIV/a8T73WVo6vLEpMgZ2P/Ylsukk58r7EFv0fcxPR8VHf4xBqH3qy4N+bMPhA4x662vstJ/B7GvvjkF1ycWWV4lDAjTVW9tELSfQQk+GdX10YmH6eVcee3NRsze+18s1wSAb//4sO783l5lvfY4042P1nB8y6htrmXnhxxI2mixlA3/Tzn2c3fMWNVUfgPTBj1BxVp1F47yvca5ZxV+W3IUxJgZVWmvV359nYdtLbBp1LeXjZiasjrKR49g45pPM6niTTf94bsCPM8EgB3/zKVymi7QrH7Ts5GYyyygIjeNvqrHuylcmGCTPtBLIHNqBD1Cx+CcYhLqnvhi3fe5Y/yZbv38O0//+CbL9Layc8R3Kv76B6ks/FbFRXzylbODvff0RMsRH5fmn3y6h/KIvcSB/HlfU3cfLb/zTwuqs197RQfE/bqNWSpl29Z2JLoeZi/+HA1JB3t+/MeBjrKt+/3+c1bmMdeM/z+hJ1TGuMDlkl4QmX7VbOPmqo92DQ/xIlnV9dGKlrOpMNp7xaWZ2vsX6V5fGfH/7tq2l7PcfpsR3kBWTvk7h195m7hWfJy3dEfN9D0RKBn6n18eEQ7/ngHM8rlFRXP7NZqPi+kcJ2ByM+Nt/cqC+xbIarbbqsW8z2hyi/YK7SHcmvj99hjOL5kXfCR1jffK7p1z/0O6tTNn4fTY7ZjB38TfiUGFyKAxPvuppsm4svqch1HLZlj08hiXP+ug32GurouRft9PV0Raz/bQ0HMb+5FV4JQP/J15l3ke/OuQ+hVoS+CJysYhsF5GdInJbhOUZIvJkePkKERltxX5P11tvvsoE2UfgrI9HvS27uwLv++5liuxh9SO3EggOvUM772zdwIJDD7PJfR7jzv5Qoss5Zvp5V7Iu62xm7P4Vhw/s7He9gN9P69JPEEQo/NhDCf9YPJzkFZbSZRzgse4YfntTLQAZ7uER+I4MJx3n302FOcL6pbfHZB++Hi+HHvwwJcEGGi59iLKR42Kyn2hFHfgiYgfuAy4BJgNXicjkPqvdCDQbY84EfgIk9GxbcM1v6MbBqHOjD3yAojlXsnfUh7m84xn+9IcnLNmmVQKBIO3PfRG/pDHqmv9LdDknKP3IT7ARpObJL/W7zsrHv80k3xa2zbp9yP4iDVVis9FgK8LRYd3kqy5PaJhnZt7gLtWYSFPOvpTVuRcw+8BvOLBjg6XbNsEga+//JFN6NrJx9p1MnHOBpdu3khXv8OcCO40xu40xPcATwOV91rkcWBK+/QxwvogkpBnN7pojLOh8nX1lFyKZ1g0rG3X1vdQ7RjB/w9fZvHOvZduN1pt/eIBZvrXsmf4lcktOvxFXrFSMnsC60Tcyq/0fESfJ7Nr4FrN33cfa7EVUfyA2DaWSncdRgsvCyVc94T46rjxr2yrE2uirf4KXdFqe+SImGLRsuyue/D7zmp5nWcW1VF/2Gcu2GwtWBH4l0PsA4cHwfRHXMcb4AQ9wwhkfEblJRFaLyOr6+noLSjvRlleWkCNdlJwbXWOqviQjm6yrHqFIPDQsvZn2bp+l2z8dh+sOM3Xj99ntGM/Uy/t/B51oMxffzkEpI/fvX6fH233s/u6uDmx/uJlWyWHsMG6MlmhdmeXk+62bfBVsD/1uugd5MfZEKyobyeaJn2eady1r/7rk1A8YgI1/f5Y52+5hXdbZzLvxXku2GUtD6jfIGPOgMabaGFNdXGz98UFfIMiIPc9wOH0E+ROtn8acM3YOh2Z+mXMDy1j9/Qv5zQP38PqGXXj9iWnBsO2xr5BPK1kf+hliH7ptk5yZLhrP+TajggdZ+9T3jt2//tEvMya4j0Pn/pC8YRYuQ0kgu4Ii04yvx6L5Ih0NeE06ruyhO9O2P9VX3sou+1iqVt5Be2tzVNvat309o/9+C3vTRjP+M48Pi3NLVgT+IaCq1/cjwvdFXEdE0gA30GjBvgdlxcplnMU22iZdHW6lar1RH7iNmmmfZWb6fq6tvZOFz83jX9+5kN/96h7+sWkXPX7rPkqezIo3XmKR509srlpM2cT5cdlnNGa8ZzHrsxYwfef9HDm0h81vvcDcw0+wovCDljRaS2W2vBHYxNBQu9eS7dm7m2iR3GH5iSst3YHv4nsooYm3Hz/90V4tDYexP/FRfKTjuu7pYdPew4qf2CpgnIiMEREHsBh4vs86zwPXhW9fCfzNJGCmUvuyR/BjZ8z5N8ZuJzY7Ff/xfdxf34H/uheon3AVs9P38rFDdzL/mXn86873svRX9/Dmpl34ArEJ//aubvJf/28abQVMunr4zEYt+fC92Ahy6LHPUvDy56mxlTHthqF3onm4cRaFxuK3HN5ryfbSvc20D+E+Oqcycc4FrMx/P9W1S9mzZfAXWuk9Iqf+0oeH1UCCqD/nG2P8InIL8BJgBx42xmwWkTuA1caY54GHgN+KyE6gidAfhbiqa26l2vMSuwvOYbw7DocHbDbSxpzNiDFnQ/Cn+PYt58iyJ5i5+wXyDt2J95m7+JfM4HDlxYxceCVzJ44mzX76f399gSCN7T00tLRy4IV7uMTsZdd591OcNXx+MSvGTGTZqBtYsP9BAkbY8f5nGTEMDxsMNe7S0QB01O+zZHuZvmY60/Is2VaijLv6h7Tf93e6/vBfmIlvDPjTigkGWfvAp5jXs5FVs77PnCE8IicSSw7sGmNeAF7oc9/tvW53Awn9XL72laVcIq0Ezo7hu/v+2Gykj1lI1ZiFELyXnn3LqXvrSWbs+QvvPvRdvE/dzVsyg8NVFzNywZVUTwiFvzEGT5eP+tZumpoaaWs4RHdzDf7WWmirI63zCE5vPTm+RvKCzZRIC1OlnanA9tyzmbBo+PWIn7n4f9ly7zJaR17A/OrzE11OUiisDE2+8ls0+So70EJtVtWpVxzC8ovLWTn1y8x9+9us+tP9zLl8YJfGXPHkXcxv/CPLKq5lwQAfM5QM3TN5FgqGr2rVaCuidFbsL/ZxUjYbjjELGTlmIQR/gnffCmrfeoKpe15g0YHv0rP/LpbLdHpsTnL9TRTTTKV4GCcnnnDzkYYnrZDOrCJ6MsdR7yqhMaeMjIJKxp99VczOU8SSMyubyV8f2m0qhpvs3HxayUJarZl85Q62csA59PvonEr1FV9g+7YnGLvuLjyLPoI7v+ik62/6x3PM2fYD1rkWDosROZGkROCve3sTc/zr2D7h0xQmqENkRDYbGWMWMHrMAgj+hO59K6l9aynj970K2Oh2FeHPGkNtdilp7nKcBRVkF40gK78cySkjPTOfomEY6ir+Gm3FODpro96Ot7uTbOnCDIPGaadis9uxf+An5D13Kasf+wrzbnmk33X3bV/PqNc/xz77KMZ/ZumwGJETSUoEfv0bDwMw5oIhPHHHZsM5Zj5jxgz9ETVq+Gl1lJDjjX7ylaexjhKGTx+dUzlzxtmsePNDVNc/x84NN3LmjHedsI6nsQ77E4vxkU7mdU8NmxE5kQy/cVWD5OnoZlr9n9iVU42zZGhf0FipWOnOKqcwEP1kxrbG0KcER25yBD7AxGt+QLO4Cf7pv064UIyvx8uBBz9CSbCe+vc9RPmoCQmq0hpJH/irXv89ldIw6KtaKZVMgrmV5NNKd2d7VNvpbAl9SshwD6+2Cifjzi9i76zbGO9/h1XP3Xvs/tCInJuY6l3Phpl3JMXV1ZI+8J0bf4dHcqlacGWiS1EqYdLyQqNq6mv2RLUd77E+OiVR1zSUzH7/p9nsmMaEzT+muT70KWblU3czr/EPLCu/ljkf/FyCK7RGUgf+tl27metdxsGqDyDpQ6svtVLxlFkcapznORxd4PvbQoGfW1gedU1DidhsuK64F5fpYsfjX2bTP56jeuvdrMtayLxP3pvo8iyT1IG/928P45AAVVFc1UqpZJBXGjp/1Rnl5CvT0UDACO6C5HqHDzB6UjVryhczt/kvjP3bZ9hvH8W4m4dHj5yBStrA7+7xM+7Qc+xxTiZ31IxEl6NUQhVVhgI/0HIwqu3YuhrxSE5ShWBv0675HnUU4hUHzuueIjvXuhbqQ0HSDstc+eZfWcQhdsy4JdGlKJVwzkwXTeRib4tu8lV6dzNtNjfDfxR+ZK6cPNo/+Qp+sVFemXyj+pI28INrltCJkzPOs+aqVkoNd432Epydh6PahtPXTMcw76NzKqUjzkh0CTGTlId0DtbWMbfjH+wuuxibMyfR5Sg1JLRnlJLbE93kK5e/he705DrMkUqSMvC3vPooWeKl9NybEl2KUkNGjyv6yVe5QQ++JOijk6qSLvADQUPl7qc5mD6a4okLE12OUkOGya0kR7pobTm9aw8F/H5yTRvBJOijk6qSLvCP7N3CFLOD1knDs1ukUrGSVhCafNVYs/u0Ht/aXI9dDOI6eVdJNXQl3Unb8rFTCH52JROyk2+csFLRyA5f+aq1bi9MnjPox7c21pIPpOUkTx+dVJN0gQ9gKxneDY6UioW8itCFULob9p/W4zuaw310cpOnj06qSbpDOkqpyIrKRhIwQtBzepOvvK1HAMjK03f4w5UGvlIpIi3dQYMUkNZWc1qP7wk3TstJsj46qUQDX6kU0pxWQlbX6V35KtgeCnx3oR7SGa408JVKIR3OUty+I6f1WOlqos1kkuHMsrgqFS8a+EqlEJ+rnKJgAyYYHPRj07qbaLXlxqAqFS8a+EqlkrwqnOKjuWHwh3UyvE202/Osr0nFjQa+UinEURC6EEpT7eAvhJLlb6FL++gMaxr4SqWQnJLQ5Ku2ur2Dfmx2wENPhrZVGM408JVKIfnloR7v3sbBTb4ywSB5ppWgU9/hD2ca+EqlkILiCnpMGgxy8lVHu4cM8YH20RnWNPCVSiE2u516WyHpHYM7aetpCF04xZats2yHMw18pVJMS3oJWd2Du/JVR0u4j45bA38408BXKsV0ZpaT7xvchVC6WkJ/IJy52oV2ONPAVyrF+LMrKDKNBPz+AT/G6wn9gcgu0D46w5kGvlIpxuYeQZoEaTi8b8CPOdZHp6gsVmWpONDAVyrFOItCV75qHsTkK9PZiNek48p2x6osFQca+EqlmJyS0Fj8jvqBv8O3dzXhkRzEppExnOlPT6kUUxi+8pWv6cCAH+PwNtGmfXSGPQ18pVJMbl4hHcYJrYcG/JhMXzOdaXmxK0rFRVSBLyIFIvKKiOwI/xtx3rWI/FVEWkTkz9HsTykVPbHZaLAX4xjE5CtXwINX++gMe9G+w78NeM0YMw54Lfx9JPcAH49yX0opi3gcJeR4Bz75yh304M/QPjrDXbSBfzmwJHx7CfDBSCsZY14D2qLcl1LKIt2ZZeT7Bzb5ytvdSY50YbIKY1yVirVoA7/UGHP0c+FhIKqLXYrITSKyWkRW19cPbiagUmrgAjmVFNGCt7vzlOt6GkNtFbSPzvB3ysAXkVdF5O0IX5f3Xs8YYwATTTHGmAeNMdXGmOriYn1xKRUraXkjAGioOfXQzLamUOA7cvV3crhLO9UKxpgL+lsmInUiUm6MqRWRcuD0ro6slIorZ9FoAFoO76Fy7KSTrtvZHG6cpn10hr1oD+k8D1wXvn0d8Mcot6eUigN3WejKV10Np36H720NBb4rP6ojtmoIiDbw7wLeKyI7gAvC3yMi1SLy66MricibwNPA+SJyUEQuinK/SqkoFFeeAYCv+dSTr/xtofNpuYXaOG24O+UhnZMxxjQC50e4fzXwyV7fnxPNfpRS1sp05dBCNraBTL7qaCRghNx8PYY/3OlMW6VSVKO9GGfnqSdfSVcjrZKDPS2q94dqCNDAVypFtTlKyfHWnXK99O4mWm3aJTMZaOArlaK8rgoKg6ee7+L0NdOhfXSSgga+UikqmFOBmw462z0nXS/L30J3urZVSAYa+EqlqPT80IVQ6g+d/EIouUEPPqcGfjLQwFcqRWUVh8bitx7e3e86wUAAt2kjmKl9dJKBBr5SKSqvLHTlq66G/f2u09pcj10M4iqKV1kqhjTwlUpRRRWjCRoh0HKw33U8jaEWymk5OgY/GWjgK5WiHBlOmsSNva3/yVcdzaHAz9DGaUlBA1+pFNaUVkJmV/8XQvF6QuP0s/K0j04y0MBXKoW1Z5SR6+u/yW1PawMAOdpHJylo4CuVwnpc5RQH6jHBYMTlwY7QxCx3ob7DTwYa+EqlMnclWeKltTnyjFvpbKTdZJLhzIpzYSoWNPCVSmHp+SMBaKiJPPkqrbsJjy03niWpGNLAVyqFuUpCgd92ZG/E5Y6eZtrtefErSMWUBr5SKaywInQhFG9j5MlXWb4WurSPTtLQwFcqhRWUjMBn7AT7mXyVE2jBl6GBnyw08JVKYfa0NBqkgPT2EydfmWCQPNNKwFmQgMpULGjgK5XiWtJLyOo+8UIonR2tZIgPtI9O0tDAVyrFdTjLyIsw+crTEJqBa8vWtgrJQgNfqRTnyy6nKNhAMBA47v72cB8dh/bRSRoa+EqlOFteFQ4J0FR//HH8rpZQ4Ge6SxJRlooBDXylUlxGQWgsflOfyVdH++hkF2gfnWShga9UissuCV35qv3I8YEfaA/30Skqi3tNKjY08JVKcYUVYwHoaTpw3P2mo4Eek4Yr252IslQMaOArleLyCkvpNungOf4Yvr2riRbJRWwaE8lCf5JKpTix2ai3FePoqDnufkdPM23aRyepaOArpfA4SnD1mXyV2dNMZ1peYgpSMaGBr5SiK7OcfP/xk69cgRa8Du2jk0w08JVSBLIrKDLN+Hq8x+5zBz34tY9OUtHAV0phyxuBTQwNtXsB6PF2kyNdmKzCxBamLKWBr5TCWRiafNVyeC8Ankbto5OMNPCVUuSWjQago34fAK3hwE/P0U6ZyUQDXyl1bPKVPzz5qrM5NGLH6S5NWE3Kehr4Sily3AW0koW0hiZfedtCI3Zc+Rr4yUQDXykFQKOtGEdnLQD+1lAfndxCbZyWTKIKfBEpEJFXRGRH+N8TBu2KyFkiskxENovIRhH5aDT7VErFRqujhBxvePJVRwNBI+Tm60nbZBLtO/zbgNeMMeOA18Lf99UJXGuMmQJcDNwrInlR7lcpZbHurHIKA6F39tLViEdysKelJbgqZaVoA/9yYEn49hLgg31XMMa8Y4zZEb5dAxwB9G2DUkNMMKeCfFrp7mwn3dtEq027ZCabaAO/1BhTG759GDjpGR4RmQs4gF39LL9JRFaLyOr6+vooS1NKDUZafhUA9TV7yOhppkP76CSdU35eE5FXgUhXQPhG72+MMUZEzEm2Uw78FrjOGBOMtI4x5kHgQYDq6up+t6WUsl5mcehCKJ7De8j2t9CUOTqxBSnLnTLwjTEX9LdMROpEpNwYUxsO9CP9rJcL/AX4hjFm+WlXq5SKmbzSMQB0NeynPNhKnfbRSTrRHtJ5HrgufPs64I99VxARB/B74DfGmGei3J9SKkaKKkOB72/ah9u0EczUPjrJJtrAvwt4r4jsAC4If4+IVIvIr8PrfARYBFwvIuvDX2dFuV+llMWcmS6ayCWjcRtpEkRc2lYh2UQ15soY0wicH+H+1cAnw7d/B/wumv0opeKjyV5MWed2ANKyNfCTjc60VUod05ZRSoUJTb7KcJckuBplNQ18pdQxPa6KY7ez8rSPTrLRwFdKHWNyK4/dzi7QwE82GvhKqWPSCqqO3c4r0sZpyUYDXyl1THZRaPJVu8kkw5mV4GqU1TTwlVLH5JWHxuJ7bLkJrkTFgga+UuqYovJRBIzQYdfGaclIA18pdUxauoMGKaAr/YRLW6gkoM2ulVLH2T/zVpz5esI2GWngK6WOM+fyzya6BBUjekhHKaVShAa+UkqlCA18pZRKERr4SimVIjTwlVIqRWjgK6VUitDAV0qpFKGBr5RSKUKMMYmuISIRqQf2RbGJIqDBonJiQeuLjtYXHa0vOkO5vlHGmOJIC4Zs4EdLRFYbY6oTXUd/tL7oaH3R0fqiM9Tr648e0lFKqRShga+UUikimQP/wUQXcApaX3S0vuhofdEZ6vVFlLTH8JVSSh0vmd/hK6WU6kUDXymlUsSwDnwRuVhEtovIThG5LcLyDBF5Mrx8hYiMjmNtVSLyuohsEZHNIvKFCOu8W0Q8IrI+/HV7vOrrVcNeEdkU3v/qCMtFRP4v/BxuFJFZcaxtQq/nZr2ItIrIF/usE9fnUEQeFpEjIvJ2r/sKROQVEdkR/jfi9QFF5LrwOjtE5Lo41nePiGwL//x+LyJ5/Tz2pK+FGNb3LRE51Otn+L5+HnvS3/cY1vdkr9r2isj6fh4b8+cvasaYYfkF2IFdwFjAAWwAJvdZ57PA/eHbi4En41hfOTArfDsHeCdCfe8G/pzg53EvUHSS5e8DXgQEmA+sSODP+zChSSUJew6BRcAs4O1e9/0AuC18+zbg7giPKwB2h//ND9/Oj1N9FwJp4dt3R6pvIK+FGNb3LeDWAfz8T/r7Hqv6+iz/EXB7op6/aL+G8zv8ucBOY8xuY0wP8ARweZ91LgeWhG8/A5wvIhKP4owxtcaYteHbbcBWoDIe+7bY5cBvTMhyIE9EEnHB0/OBXcaYaGZfR80Y8wbQ1Ofu3q+zJcAHIzz0IuAVY0yTMaYZeAW4OB71GWNeNsb4w98uB0ZYvd+B6uf5G4iB/L5H7WT1hbPjI8BSq/cbL8M58CuBA72+P8iJgXpsnfAL3gMUxqW6XsKHkmYCKyIsXiAiG0TkRRGZEt/KADDAyyKyRkRuirB8IM9zPCym/1+0RD+HpcaY2vDtw0BphHWGyvP4CUKf2CI51Wshlm4JH3J6uJ9DYkPh+TsHqDPG7OhneSKfvwEZzoE/LIhINvAs8EVjTGufxWsJHaKYAfwM+EOcywN4lzFmFnAJ8DkRWZSAGk5KRBzAZcDTERYPhefwGBP6bD8kxzqLyDcAP/BYP6sk6rXwS+AM4CygltBhk6HoKk7+7n7I/y4N58A/BFT1+n5E+L6I64hIGuAGGuNSXWif6YTC/jFjzHN9lxtjWo0x7eHbLwDpIlIUr/rC+z0U/vcI8HtCH517G8jzHGuXAGuNMXV9FwyF5xCoO3qYK/zvkQjrJPR5FJHrgfcD14T/KJ1gAK+FmDDG1BljAsaYIPCrfvab6OcvDfgQ8GR/6yTq+RuM4Rz4q4BxIjIm/A5wMfB8n3WeB46OhrgS+Ft/L3arhY/3PQRsNcb8uJ91yo6eUxCRuYR+HvH8g+QSkZyjtwmd3Hu7z2rPA9eGR+vMBzy9Dl/ES7/vrBL9HIb1fp1dB/wxwjovAReKSH74kMWF4ftiTkQuBv4buMwY09nPOgN5LcSqvt7nhK7oZ78D+X2PpQuAbcaYg5EWJvL5G5REnzWO5ovQCJJ3CJ29/0b4vjsIvbABnIQOA+wEVgJj41jbuwh9tN8IrA9/vQ+4Gbg5vM4twGZCIw6WAwvj/PyNDe97Q7iOo89h7xoFuC/8HG8CquNco4tQgLt73Zew55DQH55awEfoOPKNhM4LvQbsAF4FCsLrVgO/7vXYT4RfizuBG+JY305Cx7+Pvg6PjlyrAF442WshTvX9Nvza2kgoxMv71hf+/oTf93jUF77/0aOvuV7rxv35i/ZLWysopVSKGM6HdJRSSg2CBr5SSqUIDXyllEoRGvhKKZUiNPCVUipFaOArpVSK0MBXSqkU8f8EPV8qhstdqgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Active power\n",
    "# From complex power\n",
    "vmag = result['bus'][:,VM]\n",
    "vang = result['bus'][:,VA]*np.pi/180  # In radian\n",
    "vcomplex = vmag*np.exp(1j*vang)\n",
    "Sf = np.diag(Cf@vcomplex)@np.conjugate(Yf)@np.conjugate(vcomplex)\n",
    "\n",
    "# From active power\n",
    "Pf = -(np.diag(Cf@vmag)@np.diag(Ct@vmag)@(np.diag(gs)@np.cos(A@vang) + np.diag(bs-bc/2)@np.sin(A@vang)) - gs*(Cf@vmag**2))\n",
    "Qf = -(np.diag(Cf@vmag)@np.diag(Ct@vmag)@(np.diag(gs)@np.sin(A@vang) - np.diag(bs-bc/2)@np.cos(A@vang)) + bs*(Cf@vmag**2))\n",
    "plt.plot(Qf)\n",
    "plt.plot(np.imag(Sf))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 352,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-6.37582737e-04, -1.26220922e-03, -6.55299381e-04, -4.89311899e-04,\n",
       "       -3.84250664e-04, -7.79263340e-06, -7.10542736e-15, -8.88178420e-16,\n",
       "       -2.22044605e-16,  0.00000000e+00,  0.00000000e+00, -5.23886490e-16,\n",
       "        0.00000000e+00, -8.88178420e-16, -2.80331314e-15,  0.00000000e+00,\n",
       "        4.44089210e-16,  0.00000000e+00,  0.00000000e+00,  4.44089210e-16])"
      ]
     },
     "execution_count": 352,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.imag(Sf) - Qf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 298,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.00000000e+00,  5.55111512e-17,  2.77555756e-16, -3.33066907e-16,\n",
       "        1.11022302e-16,  8.88178420e-16, -2.77555756e-17,  4.30211422e-16,\n",
       "       -8.67361738e-19,  6.76542156e-17,  3.33066907e-16,  0.00000000e+00,\n",
       "       -2.22044605e-16,  4.44089210e-16,  5.55111512e-17, -4.44089210e-16,\n",
       "        1.11022302e-16, -3.33066907e-16, -6.66133815e-16,  5.55111512e-17])"
      ]
     },
     "execution_count": 298,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Pf - pf_result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 258,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f9367f02fb0>"
      ]
     },
     "execution_count": 258,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAx1UlEQVR4nO3deVzVVf7H8de5cAHZlE1BQUHEJVQUEZdcyrJFK9utpsU2f9M47U3Tb/pNU04zUzPTMlNN5ZhpZeloZlo22b6YC6iIuAIKsqmALCKy3Mv5/cHN0EDBu3wv3M/z8eDBvfd7vt/z9svl4+Wc76K01gghhOj6TEYHEEII4RpS8IUQwkNIwRdCCA8hBV8IITyEFHwhhPAQ3kYHaEt4eLiOjY01OoYQQnQqmzdvLtNaR7S2zG0LfmxsLOnp6UbHEEKITkUpld/WMhnSEUIIDyEFXwghPIQUfCGE8BBuO4YvhPAcjY2NFBYWUldXZ3SUTsPPz4/o6GjMZnO715GCL4QwXGFhIUFBQcTGxqKUMjqO29NaU15eTmFhIXFxce1eT4Z0hBCGq6urIywsTIp9OymlCAsL6/BfRFLwhRBuQYp9x5zN/upyBb+q/BDr3/wt2RnfGR1FCCHcSpcr+CZvM2PyXqdsy2qjowghRKtmzZrF8uXLXd5vlyv4Qd1DOeAVg3/pVqOjCCGEW+lyBR/gcPdh9D2+C93UZHQUIUQn8tZbbzF8+HCSkpK45ZZbyMvLY8qUKQwfPpwLLriAAwcOAM2f0O+55x7Gjh1L//79+frrr7njjjsYMmQIs2bNOrG9wMBAHnzwQRITE7ngggsoLS39WZ+bN29m8uTJjBo1iosvvpiSkhKqqqoYNGgQe/bsAeDGG2/k3//+t93/PocclqmUWgBcBhzWWg9tZfl5wIfAfttLK7TWcx3Rd6v6pBBS8TGF+3YSPeBncYQQbuyp1TvYWVzt0G2e0zuYP1yeeNo2O3bs4Omnn+aHH34gPDycI0eOcNttt534WrBgAffddx8rV64EoKKigvXr17Nq1SquuOIK1q1bx/z58xk9ejQZGRmMGDGCY8eOkZKSwgsvvMDcuXN56qmnePnll0/02djYyL333suHH35IREQES5cu5fHHH2fBggW8/PLLzJo1i/vvv5+Kigruvvtuu/eDo47DXwi8DLx1mjbfaa0vc1B/pxUx5FzIgoM7v5OCL4Roly+//JLrrruO8PBwAEJDQ1m/fj0rVqwA4JZbbuHRRx890f7yyy9HKcWwYcPo1asXw4YNAyAxMZG8vDxGjBiByWRi5syZANx8881cffXVJ/W5Z88esrKymDp1KgBWq5WoqCgApk6dyrJly5gzZw7btm1zyL/RIQVfa/2tUirWEdtyhL6DRnFM+2E9sAm4x+g4QogOONMncXfh6+sLgMlkOvH4x+cWi6XVdU49lFJrTWJiIuvXr/9Z26amJnbt2oW/vz8VFRVER0fbndmVY/jjlFLblFKfKKVa/YkqpWYrpdKVUumtjXW1l5e3N3m+gwit3H7W2xBCeJYpU6awbNkyysvLAThy5Ajjx49nyZIlACxevJiJEyd2aJtNTU0njsZ59913mTBhwknLBw0aRGlp6YmC39jYyI4dOwB44YUXGDJkCO+++y633347jY2Ndv37wHWXVtgC9NNa1yilpgErgYRTG2mt5wHzAFJSUrQ9HVaHjyCl6B3qamvw8w+0Z1NCCA+QmJjI448/zuTJk/Hy8mLkyJG89NJL3H777fztb38jIiKCN998s0PbDAgIYNOmTTz99NP07NmTpUuXnrTcx8eH5cuXc99991FVVYXFYuGBBx7A29ub+fPns2nTJoKCgpg0aRJPP/00Tz31lF3/RqW1XXX1pw01D+l81NqkbStt84AUrXVZW21SUlK0PTdA2br2HUb+MIfdly5j8JiLzno7Qgjn27VrF0OGDDE6hsMFBgZSU1PjtO23tt+UUpu11imttXfJkI5SKlLZBq+UUqm2fsud2WfM8EkAVOb8fGxMCCE8kaMOy3wPOA8IV0oVAn8AzABa69eAa4F7lFIW4Dhwg3bUnxZtCI/sS7HqiU/JZmd2I4QQbXLmp/uz4aijdG48w/KXaT5s06VKAhPpc1QmboUQArrombY/aowaRSRllBbnGR1FCCEM16ULfo+EcQAUbJcrZwohRJcu+LFDx9GgvajP22R0FCGEMFyXvsWhX7cA9prjCS6XK2cKITrmySefJDAwkOrqaiZNmsSFF15o1/YyMjIoLi5m2rRpDkrYcV36Ez5ARchw4ur3YmlsMDqKEKITmjt3bqvF3mq1dmg7GRkZrFmzxlGxzkqXL/hefVPxV/Xk795idBQhhJv705/+xMCBA5kwYcKJSxO3vFlJbGwsv/3tb0lOTmbZsmWsXbuWcePGkZyczHXXXXfiMMy0tDTGjx9PUlISqampVFVV8cQTT7B06VJGjBjxszNuXaVLD+kARCVOgs1Qtvt74oeNNTqOEOJMPnkMDjr4cOrIYXDpM6dtsnnzZpYsWUJGRgYWi4Xk5GRGjRr1s3ZhYWFs2bKFsrIyrr76aj7//HMCAgJ49tlnef7553nssceYOXMmS5cuZfTo0VRXV+Pv78/cuXNJT08/6fLIrtblC37v2EFUEIyp6Owv0yCE6Pq+++47rrrqKvz9/QG44oorWm334+WON2zYwM6dOzn33HMBaGhoYNy4cezZs4eoqChGjx4NQHBwsAvSt0+XL/jKZCK/2zn0rM4yOooQoj3O8EncaAEBAUDzpY2nTp3Ke++9d9Ly7dvd92TPLj+GD3C85wj6NRVQVdHmtdqEEB5u0qRJrFy5kuPHj3P06FFWr1592vZjx45l3bp15OTkAHDs2DH27t3LoEGDKCkpIS0tDYCjR49isVgICgri6NGjTv93nI5HFPzA+PEAHJATsIQQbUhOTmbmzJkkJSVx6aWXnhiSaUtERAQLFy7kxhtvZPjw4YwbN47du3fj4+PD0qVLuffee0lKSmLq1KnU1dVx/vnns3PnTkMnbR12eWRHs/fyyC1VV5YT+EI8G2P/h3G3P+uQbQohHKerXh7Z2dzy8shGC+4RxgGvaPwPy6GZQgjP5REFH+Bw8DD6Ht+FbmoyOooQQhjCYwo+0aMJ4ShF+3YanUQI0Qp3HV52V2ezvzym4IcPbj5W9uBOmbgVwt34+flRXl4uRb+dtNaUl5fj5+fXofW6/HH4P+o3eBS12hfrgTTgHqPjCCFaiI6OprCwkNLSUqOjdBp+fn5ER0d3aB2PKfhe3t7s9x1EaGWm0VGEEKcwm83ExcUZHaPL85ghHYDq8BHENu6jrta97jMphBCu4FEF3y92DGZlJS9rvdFRhBDC5Tyq4McMnwRAZfYPBicRQgjX86iCHx7ZlxIiMJfICVhCCM/jUQUfoDhoKH1qdhgdQwghXM7jCn5jZDKRlFJWnG90FCGEcCmPK/g9BjZfObMgS07AEkJ4Fo8r+LFDx9Ggvajbv9HoKEII4VIeV/D9ugWQZ44nuHyr0VGEEMKlPK7gA1SEDCeufi+WxgajowghhMt4ZMH3ihmNv6onf7ccnimE8BweWfCjhk4EoGz3OoOTCCGE63hkwe8dO4QKglFFjrmFohBCdAYeWfCVycSBbkPoVb3d6ChCCOEyHlnwAWp7jqRfUwFVFWVGRxFCCJfw2IIfGD8WgAPb5QQsIYRn8NiC32/4JJq0oiZ3g9FRhBDCJTy24Af3COOAVzT+h+UELCGEZ/DYgg9wOHgYfY/vRDc1GR1FCCGcziEFXym1QCl1WCmV1cZypZT6p1IqRymVqZRKdkS/9tJ9UgjhKMV5u4yOIoQQTueoT/gLgUtOs/xSIMH2NRt41UH92iV88LkAlMiVM4UQHsAhBV9r/S1w5DRNZgBv6WYbgB5KqShH9G2PfoOTqdW+WA9sMjqKEEI4navG8PsABS2eF9peO4lSarZSKl0plV5aWur0UN5mH/b7DiKkUk7AEkJ0fW41aau1nqe1TtFap0RERLikz+qwEcQ15lJXW+OS/oQQwiiuKvhFQEyL59G21wznFzcGs7KSl7Xe6ChCCOFUrir4q4BbbUfrjAWqtNYlLur7tGJsV86szJaCL4To2rwdsRGl1HvAeUC4UqoQ+ANgBtBavwasAaYBOUAtcLsj+nWE8N79OEgE5pLNRkcRQginckjB11rfeIblGpjjiL6coSgwkT41O4yOIYQQTuVWk7ZGaYxKJpJSyorzjY4ihBBOIwUf6JEwHoACOQFLCNGFScEHYoeOo0F7UZe30egoQgjhNFLwAT//QPLN/QkuyzA6ihBCOI0UfJsjPYYTV78Hq8VidBQhhHAKKfg2Xn1T8Vf15O+WwzOFEF2TFHybyHMmAFC663uDkwghhHNIwbfp0/8cKghCFaUbHUUIIZxCCr6NMpk40G0IPavlyplCiK5JCn4LtT2T6WstpLqy3OgoQgjhcFLwWwiMH4tJafIzZRxfCNH1SMFvoe+w5itn1uT+YHASIYRwPCn4LXQPCSffFIP/4a1GRxFCCIeTgn+KQ8HD6Ht8F7qpyegoQgjhUFLwT6H7jCKEaorzdhkdRQghHMoh18PvSsIHT4AdUP3eXRSEDIOQfvj1jCekdwI9+w6kW0CQ0RGFEOKsSME/ReyQFDb1mEZE9U76HVqJ/+F62PPT8jJ6UGbuTU23PjR274dXaCyBkQMIixlIRFQsJi8v48ILIcRpqOabUbmflJQUnZ5u7FmvuqmJI6XFlBbsoeZgLo1l+/CqzCegtpDQhhJ66jK81E/7r0F7c8jUkyN+0fic9whDxlxsYHohhCdSSm3WWqe0tkw+4Z+GMpkI6xVNWK9o4IKfLW+or6OkIIeK4mxqD+XSdGQ/vkcPEF2TRdCam9le9zrDJl/t+uBCCNEKKfh28PH1I3rAUKIHDD3p9fJDhZS8Pp1BX95NhqWBERfcYFBCIYT4iRyl4wRhvaIJm7OWfHMcid/+ii3/XWh0JCGEkILvLN3DetHr15+S6zOI4esfJH3160ZHEkJ4OCn4ThTcI4yY+z5hj+9QktN/S9qKfxgdSQjhwaTgO1lAUA/637+GrG7JjM58go1LnzU6khDCQ0nBd4FuAUEMeuAjtvqPZ8yuP7Nh8VNGRxJCeCAp+C7i6+fP0AdWsiVwMmOzn2f9wseMjiSE8DBS8F3I7OPL8PuXk9b9Isblvcr6+Q/KRdqEEC4jBd/FvM0+jLpvCZtCLmNc4QI2vv4rKfpCCJeQgm8Ak5cXKb9+i43h1zD20Hts+tedNFmtRscSQnRxUvANYvLyIvVX89kQ+QvGlK0g/eVbsFosRscSQnRhUvANpEwmxsx+mfUxd5Fa8TFb/zkTS2OD0bGEEF2UFHyDKZOJcXc+x/q4OaRUf07mi9fQUF9ndCwhRBckBd9NjLvtz2wY+AjJx75l54szqDt+zOhIQoguRgq+Gxl70+/ZeM7jjDi+gf3PX0j5oUKjIwkhuhAp+G5mzPWPsjn1eeIasml89Tz2ZW00OpIQoouQgu+GRk27k4IrV2DCSuSyy9m69h2jIwkhugCHFHyl1CVKqT1KqRyl1M+uGaCUmqWUKlVKZdi+7nJEv11ZwshJqNlfUWjux8gf5rBh4e/kBC0hhF3sLvhKKS/gFeBS4BzgRqXUOa00Xaq1HmH7mm9vv54goncsfR/6ivTgCxmb9wqbX7iWutoao2MJITopR3zCTwVytNb7tNYNwBJghgO2KwA//0BGPbCs+bDNo19w4PnzKS3OMzqWEKITckTB7wMUtHheaHvtVNcopTKVUsuVUjGtbUgpNVspla6USi8tLXVAtK5BmUyMu+3PbB3/CtGN+eh555O99VujYwkhOhlXTdquBmK11sOBz4BFrTXSWs/TWqdorVMiIiJcFK3zGHnRzZRcu4omvIhZeTWbP5aRMSFE+zmi4BcBLT+xR9teO0FrXa61rrc9nQ+MckC/Hil+2FjM93zNfp8ERqU9zIb5D8mF14QQ7eKIgp8GJCil4pRSPsANwKqWDZRSUS2eXgHsckC/HiusVzT9H/6CTT2mMbbwDTKen0FtTZXRsYQQbs7ugq+1tgC/Bj6luZD/R2u9Qyk1Vyl1ha3ZfUqpHUqpbcB9wCx7+/V0vn7+jL5vMRsSHiap5ntKXjiPgweyjY7VJVTWNmCxyiGwoutRWmujM7QqJSVFp6enGx2jU9j21TL6f30v9cqXsulvMHj0hUZH6rQOV9dx0Yvfct7ACF68YaTRcYToMKXUZq11SmvL5EzbLiDp/Os4cuMa6pQf/T+aSdrKV4yO1Gk9tXonlbWNfLitmNxSOedBdC1S8LuIfoOTCZzzDXv9hjI643ekr3rN6Eidzhe7DvHx9hL+nFhEgvdhXv061+hIQjiUFPwupEd4JIMeXkuhisS8632j43Qqx+otPPHhDq4PzeWm3N/wVuDLfLi1gMKKWqOjCeEwUvC7GLOPL0Vh40mo3UZ9nRSr9nr+s70crSzjaf4FPkFEHs9humkDr3+zz+hoQjiMFPwuyHfQVPxVPdmbvzA6SqewvbCKN9ftZ1HkMnyOH4ZbV0LPRP7P/wPeT9/P4Wq5A5noGqTgd0HxqZdg0SaO7vjM6Chuz2Jt4rEVmVzvv5mRlWth8qMQnQJT/o/whkKu0N8w//v9RscUwiGk4HdBQd1DyfYZQvihdUZHcXsLf8jjcHE+c73mQ+9kmPhw84JBl0KfFH7bbSX/2ZBNxTG5ubzo/KTgd1GVvScSb8mlorTE6Chuq+BILc+t3cOCkIWYmxrg6nngZW5eqBRc8AQhllKutq7lzXXyKV90flLwu6iQYRdhUpp9mz42Oopb0lrz+w+zuEF9xrDjaaipcyE84eRG/SdD3GQe8F3Nf37YxdG6RmPCCuEgUvC7qAFJE6kmAGu2TNy25qPMEvL2ZvI778UQPwVGt3ETtgueILipkmsbP+LtDfmuDSmEg0nB76K8zT7kBIyib+UmuTXiKapqG3l61XZeC5iHt48fzHgFTG38KkSnwKDp/MrnY5Z9m8nxBrkyqei8pOB3YY2xk4mkjAPZmUZHcSvP/HcXM+uXMdiyBzX9OQjuffoVpvwf3fRxrm9YwZK0A64JKYQTSMHvwmJSLgOgZIuM4/9o0/4jbE/7hvvNK2DoNTDs2jOv1Osc1PDrucN7Le9/k06DRf5iEp2TFPwurHfcYApVFH4HOtftEP+1NpOlG/fj6Cu51lusPLkinZf9XkUF9oRpf2//yuc9hllZub52KSu2FDo0V1ezec2bbP/2A6NjiFZIwe/iikLHMLB2Kw31neNs0araeq5YdxXDP76cJ15/j5Kq4w7b9mtf7+PaijeI1UWYrvwX+Ie2f+XQ/qjkW7nJ+0tWfrVOrpffhpxt60ja+BADv7ibnG3fGx1HnEIKfhfnc+IyC18aHaVdcndnEq3KGGgq5omSX7PihXv5IN3+T/u5pTVs+foD7vD+L6TObj4yp4PU5N+gTN5cV7OYjzLl/IZTWS0W9Or7qVTBVKpg/FfeTlX5IaNjiRak4Hdx8amXYtEmqnesNTpKu1TkbALg2Mzl1A+awRyWkbDqSv74xn8oq6k/w9qt01rz9PL1POv9GpaQAXDhU2cXLrg3pjGzucrrez76/Euamtzz5kFGSfvPMyRYsskf/XsqL5tPeFM5ef/+hdxz2Y1Iwe/ignuEkeMzmLDOcpmFkgzqMRM8cCKBN72J9fq3ifet4rGCe1j63P18mlnQ4U0uSy/kiuIX6Kkq8b52Hvj4n3U8NeFBrN7+XFO9iLU75dPrjw4eyGb4nn+yzW80yZfewaCUKWxN/F+S6tLYuOgxo+MJGyn4HqAiagIDGrOpLDtodJQzCq3eSbFv/IlLHHidcwXdHkinfsClzNHvEbn8Cv6y6AOqatt31mtZTT2bPl7AVV7rYNJvoM8o+wIGhGEafy+XeqWx9rM1Dp9Y7ox0UxMHl9wLQMQNr6Bs5zSkXvswad0vZkz+v9n21TIjIwobKfgeIGRo82UWcjd9YnSU0zp6vJ54Sy41oUNPXhAQTtAti7FcvYCBPuU8tO9u3nruQb7adeZx9H988C2P63nU9UzCNOkRh+T0Gj+HenMPrjyygG+zyxyyzc5s69q3GFG7nsyEe+gdO+jE68pkYujsN9jvHUvsN/dTvH+3gSkFSMH3CANGTuao7oY1+3Ojo5xW7p7tBKvjmGOSW13uPfwauj2QTl3sFO61vk3we5fz3HsfUVNvabX917sPceHeuQR6NeJ3/Rs/XRjNXn7BeE1+mEle2/nqE8++s1h1ZTkxG54k16s/KTMf/9nybgFB+P1iMQpN7Ts3UVcr9wk2khR8D+Bt9iEncBQxFRvd+jILR7KbJ2x7DR7TdqPAnnSftZSGGa9zjvkQc3bP4o2//4YNuaUnNTveYGXz+88x2SsTWrswmp28x9zNMd8ILit/g425nvspf9c7jxCqK7FOfxFvs0+rbfr0T2TfxBcYYM0lc97dbv0e7Oqk4HuIhn6TiaKUwtztRkdpk7ZN2Ib0Szp9Q6XwGXkD3e7fRF3MRO5vXIBadBkvLV9LXWPzESGLPv6cexoWUhk1EfOY2Y4Pa+6Gz5T/JcW0l+8+edfx2z9LlsYG0la+zO405/81tzvtc0aXfkBar+sYmDz5tG1HXHAD66PvILVyDWkrXnR6NtE6KfgeIjplOgDFm9cYnKRtoVU7Kfbp3/6hl+Aoety5gvrLXiLJu4A7t9/Ma3//X5Zt3MeYrb8Dbx963Div7Quj2cmccitVftFMOzyfbQeOOKWPjtj+7QcU/mUUozMeJ/ajG8j6fpXT+mpsqMf3k4coVaEMvfmv7VonddbfyPQbxYjtf2Lvlm+clk20TQq+h+jTP5Ei1QsfN73MQk1dI/GWXI6eOmF7Jkrhm3Irfvdvor53Kg/Uv864NVMZacpBT2vHhdHs4WXGZ+r/cY4pn00fzXdeP2dQkLOdjL9ewrAvZ+Gj60kb+RdKvHrT/7O72LneORP16e89RVxTPiXnPk1gcEi71vHy9qbvXe9yRIUQvOpOuTmPAaTge5DC0LEMPLaFxoazO4HJmZonbGsxx4w8uw10jyZk9mqOX/I8EV61HI67koBRMx0bshXdRs6k1H8AFx6cz54i137Kr64sZ8Orv6TX25NJOLaV9f3vI/TRrYye8SuCZn9MqVcEsf+9jd2bHHtv46J9Oxi5bx5bAiYyYupNHVq3R3gkR2csIFRXUPDGL7BaWp9wF84hBd+D+Ay8gABVR/aWr4yO8jMnJmwHjT37jShFt7F34vtYLj1vedNByc7AZKLbxX8gznSILatfcUmXVouFjcuew/LiCFIPLiEj5GKO35PGuFv/iF+3AADCI2MIuGsN5aYwoj++hb1bvnZI37qpiSNL59CIN9E3vXRW20gYOYmMYY8zvG4zmxb+xiG5RPtIwfcg/VOnY9WK6iw3vMxC8VYa8CY09gwTtu3h4++0cfvWBA6/nKKARCaXLCD/YLlT+9qx7mPy/pLCmB1zOWSOIfeq1aQ+8B7hkX1/1ja8dz987vyIKlMwkatuJDvjO7v73/zR6wyr38rOcx6kZ5+4s97O6KsfYFOPaYwrXEDGF0vsziXaRwq+B+keEk6OeRAhB93vKoY9qnZQ5NMfvFs/tM+tKYX/tKforY6Q+eELTumiaN8utvztMhI/u4kAaw2bRz/H4P/9noQRE0+7Xq/oeLxu/5gaFUjPlTPJ3b7hrDNUlh0kfsuf2eM9mNHX2ncSmzKZGD773+R4xdP/uwcp2rfDru2J9pGC72GORE1gQONeqo6Unrmxi9TW2yZsQzo4YetGQhKnkhuYwvjiRZQcdty+ramuYP28e4lYNIHBNZtY3++X9Hg0g1HT7zpxCYMzieybALeu5jjdCH3/OvbvTDurLNnvPECgPobPVf/E5OV1Vttoyc8/EP+b30WjqFv8C44fO2r3NsXpScH3MCFDL8ZLaXI3uc9dsPbtbZ6w9Y4+ywlbNxE0fS5hqppdHzxr97aarFY2ffBP6p4fwbjit9jWYwo1szcy7vZn8fMP7PD2escNxnLzhzTiTfB/riF/95YOrZ+1bjWjKz8hvc8viEs8zYlxZ5Erb/I/iLPkkTXvzrM+Kavu+DF2bfyU9YseZ+tfL2XDa7+i7vgxh+XsKryNDiBcK37kZGrWdKNx7xdwySyj4wBQtrd5wranPRO2bqDnkHPJCprA6OJ3KDv8COE9IwGwWhppqK+joaGOxvo6LA11WBrqsTT+9L2psQFrYx1Nlnoaa6vpvvV1Uq057PEeTPklbzI6pePX7z9V9IChHLjpA/zfnUHAkqso+MUqYhLOPGdSd/wY3T9/lCLVi5E3/8XuHKdKOv861u/fxLgD89i4/DnGXH/midzKsoPkZXzF8dzv6VG2hfiGvQxRzUf8FKpIomt/IPfv6/G58S1iBgxzeObOSgq+hzH7+JIVMJKYig3opqZ2Dws4ky7OoAFvwvo7YMLWYN2nP0nAe1PxemU4dTRhxoKX0nQDunVgO4cJJT35WZKn3+2Q4ZMf9R04gvyZKwhaeiW+i6+k6NaP6NM/8bTrZCz+PWN1MdunLKTPWfx10R5jbnuGbX/PYOSOv7AnfRSDWvwHp5uaKM7bQ3HmlzQd2ECvyq3ENhUwAmjQXuwzD2RL1Ez84ifQb8T5REdEkfHFEmK/exjvt6eyOfVpRk2/yym5Oxvlrpd3TUlJ0enp6UbH6JI2Ln2GMbv+QuHN64geYPy4+danJ9LDdJy433WNn3fa+y9iPZgFXj5oLx/w9kV5+YC3D8r22GT2RXn7YPL2xWT2xcvsi7fZD5PZF28fP3r3TzyroZv22r9jIz2WXUM9vjTNWnPSVS5byt+1maglU8nsPoWUh5Y7LQ9AVfkhjr08AS9tpeyilzmavxWfoo3E1GQSQQUA1fizv9tQjvcaTfCgifRPmtjmfjp4IJvKt25hsGUXG8OuJOmuf504bNVdNdTXkbdjIw3Hqxl67uVntQ2l1GatdUqry6Tge56CnO3EvDOBjUN+x5iZvzU0S12Dhfo/9SUv8iKS7lloaBZPk5v5A+ErruWYCsB0+5rmyd0WmqxW9jwzkajGfJp+tYnQnn2cniln2zqiV8zATzXf76CECIqCk7BGj6Fn4mT6DU7p0F88jQ31bF7wIGMPLibXKw6fG992myEe3dTEocJcirK+pTF/Ez2OZBLXkI2vaiTXqz/xv996Vts9XcGXIR0PFN0/kWLVE5/8bwBjC37O3iyGqmOdfsK2M4ofPp5s6xJ6fTiT6jcv4/Bd/z3p2Pq0FS8ypnEHm5L+SKoLij3AgKRz2VW3mGNlB4gefj5RMQOIsmN7Zh9fxv7yX2z7chL9vn0Y89sXkp76NCnT73ZY5vaqralif+Y6qnPW43doC9HHdhBJBZFAnTaz3yeBrZHXYu43hj5DJzglg0MKvlLqEuAfgBcwX2v9zCnLfYG3gFFAOTBTa53niL5FxymTicKQMQwp/5zGhnrMPr6GZTmSvRGAiIGOO/JDtF/CyEnssS6mz+qbqJg/jbK7/0t4736UHTzAkB1/Z4fPcEbP+LVLMw0Zc7HDt5k05QYOJqRwcNHNpKQ9wsZ935F057+cNmzWZLVSkLOdw7u+p6kgjfDKTPpZ8khUzUchFaoo8oNT2Nd7FKGDzqXfOakM8fVzSpaW7C74Sikv4BVgKlAIpCmlVmmtd7ZodidQobUeoJS6AXgWcP6FTkSbvBMuIOjIanZv/YbBYy4yLIe1qPkM24j+IwzL4OkGpUxht/Ut+q65mdL502H2J+Qvvo9huoGga19yi4l9R4iMGUDYo9+w/s2HGFfyDrnPZeJzw6J2HanUHhWlJWR/vxxz9hria7fRj2P0A47qbuT5DSE96na69U+l77DJREdEEe2QXjvGEZ/wU4EcrfU+AKXUEmAG0LLgzwCetD1eDryslFLaXScQPEB86jSsGx6kImstGFjwe1TupMgnjjiz8z/diLYNHnMRO6wL6f/pbVS9dh6jKGN97C8ZN3CE0dEcyuzjy7j/eYVtX05sHuJ55yLSR/+RlMvO7p4JxXl7OPDDfwjK+5TB9VmkKs1BwtkdOgViUuk1ZAIxCUkMc+CRVvZwRMHvAxS0eF4InPr3+Yk2WmuLUqoKCANOulWQUmo2MBugb9+fXxtEOE73sF7sNScQUmLcZRbqGizENeaQ12uqYRnETxLHTyPLOp8Bn99JvlcMyTf+wehITnNiiOetW0hJ/03zEM9dr55xiEc3NbFvxyYOp71Pz6LPibfuozew39SPtJg7CB99NfHDxhPppn8VudWkrdZ6HjAPmo/SMThOl1ceOYHRBQupqiije0i4y/vfn72TIeoY3n1kwtZdDJ04g4LIWAKCQvH18zc6jlNFxgwg7Ddfnxji2ff3TLxvWETfU/6qsVos7En7jOqMlfQ9/BXx+hBxWrHH5xw2xD1I9NhriRswlLO/lJzrOKLgFwExLZ5H215rrU2hUsob6E7z5K0wUPehF+NduIDcTZ+QfPEtLu+/NHsjQ4DwQTJh604cNabdGfw0xDOZvt8+hM/ii0kf/UeGTrmJ3T+spiFrFQmV33MO1TRob3b5j6JwwD30P/dahkTGnLkDN+OIgp8GJCil4mgu7DcAp94VYRVwG7AeuBb4UsbvjZeQfD7HPvGjce/nYEDBbyraSiNe9IqXT/jCWElTrudgQjKHbEM8dWm/Y4RqpBp/9gaPRw2+jEETriSpnXf3cld2F3zbmPyvgU9pPixzgdZ6h1JqLpCutV4FvAG8rZTKAY7Q/J+CMJjZx5cdASPpc+TsL5lrj+6VOyg0y4StcA8nhnjeewrT0RL8h13OoLHTSHHB4ZKu4pAxfK31GmDNKa890eJxHXCdI/oSjlXfdzLRu9dTtG/HGa+p4tB+Gy3ENWST1/NCl/UpxJmYfXwZd9ufjY7hNO45lSxcJir5UgAK09ecoaVj7c/eSQ91DC+ZsBXCZaTge7iYAcM5SATm/G9c2m+Z7QzbcDnDVgiXkYLv4ZTJxIGQVAYc24ylscFl/VoLmydsIxOSXdanEJ5OCr7AK+ECgqklZ5v9N7lurx8nbJVM2ArhMlLwBfGp02nSiorMT13SX6PFSmxDNpU9znFJf0KIZlLwBT3CI8k1D6BHiWs+4e/PaZ6wNcmErRAuJQVfAFDW61wSGnZztOqI0/sq3WubsE2QCVshXEkKvgAgOPFivFUTORudf3impXArjdqLKJmwFcKlpOALABJGTaFW+9Kw9wun99W9YgeF5lhMPh25rbcQwl5S8AUAPr5+ZPuPoHe5cy+zYLFY6ScTtkIYQgq+OOF438nE6GKK9+92Wh95+3YTompkwlYIA0jBFydEjWy+zEJB+sdO6+PwnuYJ29ABqU7rQwjROin44oS+A0dwiDDM+V87rQ9L4RYatRd9BqU4rQ8hROuk4IsTlMlEfo8xxNdsxmqxOKWPIJmwFcIwUvDFSUwJU+jOMXIzHX+vW6u1idj6bCq6y4StEEaQgi9OEp96GU1aUb75A4dvO3//HkLUUVTvEQ7fthDizKTgi5OERESxLWAcQ4qWU1tT5dBtH9q1HoCwBJmwFcIIUvDFz3Q7/2F6UEPmqpccul1LUfMZtn0Gj3bodoUQ7SMFX/zM4NEXstNnGLF736Sxod5h2w06kkWhuR9eMmErhCGk4ItWNY67n0jKyFgz3yHbs1qb6CcTtkIYSgq+aNXwydewzxRLROarNFmtdm+vIK95whaZsBXCMFLwRauUycSR5DnENhWQ+eVSu7d3aLecYSuE0aTgizaNuHgWxaoX3Tb+A93UZNe2Ggu2YNEmouUMWyEMIwVftMnb7EPB4DsZZNnNro323f4wsCKLAnM/vP0CHJROCNFRUvDFaSVdPocjBNP4zXNnvY0maxN96/ZyJFgmbIUwkhR8cVp+/oHsjb2FpLo0crdvPKttFORnE6qOoqJGODacEKJDpOCLMxoy4yFqdDcqPvvrWa1/aHfzTVVCZMJWCENJwRdn1D0knKyoqxlZ9QVF+3Z1eP3GQtuE7RA5w1YII0nBF+0Sf8WjWPGicE3HP+UHlGdR6N0Xs0zYCmEoKfiiXSJ6x5IReglJpaspO1jQ7vV0UxMx9XsplwlbIQwnBV+0W9S0R/HBQvbqv7d7ncL8HMKoRkclOTGZEKI9pOCLdotJSCIjaCKJRcs4WnWkXesc3CMTtkK4Cyn4okOCLvgNwRwja9U/2tW+4cAWrFoRM0QKvhBGk4IvOiRh5CS2+45kQO4i6utqz9g+sHw7B7z74dMt0AXphBCnIwVfdJia8CARVLDt49dP2+6nCdshLkomhDgdKfiiwxLPvZxs7wSisuZhtVjabFdSkEso1Wg5w1YItyAFX3SYMpk4OmoOMbqYbZ+93Wa7kl22Cdt4OeFKCHdgV8FXSoUqpT5TSmXbvoe00c6qlMqwfa2yp0/hHpKm3kKB6k1Q+ittXjq5vmAzVq2IlglbIdyCvZ/wHwO+0FonAF/YnrfmuNZ6hO3rCjv7FG7Ay9ubkqGzSbBmk7VudattAsqzKPDqi59/kIvTCSFaY2/BnwEssj1eBFxp5/ZEJ5I0/X8oJQS+f+Fny3RTE9F1eymTCVsh3Ia9Bb+X1rrE9vgg0KuNdn5KqXSl1Aal1JVtbUwpNdvWLr20tNTOaMLZfP38yR1wG8Pqt5K99duTlh0s3EcYVTRFyhm2QriLMxZ8pdTnSqmsVr5mtGyntdaAbmMz/bTWKcBNwItKqfjWGmmt52mtU7TWKRERER39twgDJF5+P9UEcPSLv530eontksg94mX8Xgh34X2mBlrrC9tappQ6pJSK0lqXKKWigMNtbKPI9n2fUuprYCSQe3aRhTsJ6h7K+j7XMaZwEQXZ24hJaP5EX5/ffIZt38QxBicUQvzI3iGdVcBttse3AR+e2kApFaKU8rU9DgfOBXba2a9wIwmXP0ID3hSv+elTfrfy7RR4xciErRBuxN6C/wwwVSmVDVxoe45SKkUpNd/WZgiQrpTaBnwFPKO1loLfhYRHxrAt4nJGHvmEw0V5P03YBsmErRDu5IxDOqejtS4HLmjl9XTgLtvjH4Bh9vQj3F/M9N9iWriS3NV/RV/+CL2oJFcmbIVwK3KmrXCI3nGDyeg+heEl75Of9gkA3WXCVgi3IgVfOEzoRY8SoOron/FXrFrR7xyZsBXCnUjBFw7Tf+gYtnVLJZxKCryi6RYYbHQkIUQLUvCFQ5knPQRAqUzYCuF27Jq0FeJUQ8ZczIY9cwgdfrHRUYQQp5CCLxxKmUyMnfVno2MIIVohQzpCCOEhpOALIYSHkIIvhBAeQgq+EEJ4CCn4QgjhIaTgCyGEh5CCL4QQHkIKvhBCeAjVfGdC96OUKgXy7dhEOFDmoDjOIPnsI/nsI/ns4875+mmtW71HrNsWfHsppdJt99F1S5LPPpLPPpLPPu6ery0ypCOEEB5CCr4QQniIrlzw5xkd4Awkn30kn30kn33cPV+ruuwYvhBCiJN15U/4QgghWpCCL4QQHqJTF3yl1CVKqT1KqRyl1GOtLPdVSi21Ld+olIp1YbYYpdRXSqmdSqkdSqn7W2lznlKqSimVYft6wlX5WmTIU0ptt/Wf3spypZT6p20fZiqlkl2YbVCLfZOhlKpWSj1wShuX7kOl1AKl1GGlVFaL10KVUp8ppbJt30PaWPc2W5tspdRtLsz3N6XUbtvP7wOlVI821j3te8GJ+Z5UShW1+BlOa2Pd0/6+OzHf0hbZ8pRSGW2s6/T9Zzetdaf8AryAXKA/4ANsA845pc2vgNdsj28AlrowXxSQbHscBOxtJd95wEcG78c8IPw0y6cBnwAKGAtsNPDnfZDmk0oM24fAJCAZyGrx2l+Bx2yPHwOebWW9UGCf7XuI7XGIi/JdBHjbHj/bWr72vBecmO9J4JF2/PxP+/vurHynLH8OeMKo/WfvV2f+hJ8K5Git92mtG4AlwIxT2swAFtkeLwcuUEopV4TTWpdorbfYHh8FdgF9XNG3g80A3tLNNgA9lFJRBuS4AMjVWttz9rXdtNbfAkdOebnl+2wRcGUrq14MfKa1PqK1rgA+Ay5xRT6t9VqttcX2dAMQ7eh+26uN/dce7fl9t9vp8tlqx/XAe47u11U6c8HvAxS0eF7IzwvqiTa2N3wVEOaSdC3YhpJGAhtbWTxOKbVNKfWJUirRtckA0MBapdRmpdTsVpa3Zz+7wg20/Ytm9D7spbUusT0+CPRqpY277Mc7aP6LrTVnei84069tQ04L2hgSc4f9NxE4pLXObmO5kfuvXTpzwe8UlFKBwPvAA1rr6lMWb6F5iCIJeAlY6eJ4ABO01snApcAcpdQkAzKcllLKB7gCWNbKYnfYhyfo5r/t3fJYZ6XU44AFWNxGE6PeC68C8cAIoITmYRN3dCOn/3Tv9r9LnbngFwExLZ5H215rtY1SyhvoDpS7JF1zn2aai/1irfWKU5drrau11jW2x2sAs1Iq3FX5bP0W2b4fBj6g+U/nltqzn53tUmCL1vrQqQvcYR8Ch34c5rJ9P9xKG0P3o1JqFnAZ8Avbf0o/0473glNorQ9pra1a6ybg3230a/T+8wauBpa21cao/dcRnbngpwEJSqk42yfAG4BVp7RZBfx4NMS1wJdtvdkdzTbe9wawS2v9fBttIn+cU1BKpdL883Dlf0gBSqmgHx/TPLmXdUqzVcCttqN1xgJVLYYvXKXNT1ZG70Oblu+z24APW2nzKXCRUirENmRxke01p1NKXQI8Clyhta5to0173gvOytdyTuiqNvptz++7M10I7NZaF7a20Mj91yFGzxrb80XzESR7aZ69f9z22lya39gAfjQPA+QAm4D+Lsw2geY/7TOBDNvXNOCXwC9tbX4N7KD5iIMNwHgX77/+tr632XL8uA9bZlTAK7Z9vB1IcXHGAJoLePcWrxm2D2n+j6cEaKR5HPlOmueFvgCygc+BUFvbFGB+i3XvsL0Xc4DbXZgvh+bx7x/fhz8eudYbWHO694KL8r1te29l0lzEo07NZ3v+s993V+Szvb7wx/dci7Yu33/2fsmlFYQQwkN05iEdIYQQHSAFXwghPIQUfCGE8BBS8IUQwkNIwRdCCA8hBV8IITyEFHwhhPAQ/w8c72jtH2sAmgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot(np.real(Sf), label = 'complex')\n",
    "plt.plot(Pf, label = 'direct')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 259,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-8.88178420e-16,  0.00000000e+00, -7.21644966e-15, -9.54791801e-15,\n",
       "        1.11022302e-15,  8.88178420e-16,  3.01980663e-14, -6.17353793e-03,\n",
       "       -4.97473850e-03, -2.99547367e-02,  3.06421555e-14,  1.12354570e-13,\n",
       "       -2.44249065e-14,  2.84876289e-14,  4.44089210e-15, -3.90798505e-14,\n",
       "        5.10702591e-15,  4.08562073e-14,  8.88178420e-16,  7.99360578e-15])"
      ]
     },
     "execution_count": 259,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Pf - np.real(Sf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 262,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-15.26308652,  -4.23498368,  -4.78186315,  -5.11583833,\n",
       "        -5.1939274 ,  -5.06881698, -21.57855398,  -4.78194338,\n",
       "        -1.79797907,  -3.96793905,  -4.09407434,  -3.17596397,\n",
       "        -6.10275545,  -5.67697985,  -9.09008272, -10.36539413,\n",
       "        -3.02905046,  -4.40294375,  -2.25197463,  -2.31496348])"
      ]
     },
     "execution_count": 262,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 266,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 984.73691348,   15.26308652,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ],\n",
       "       [ 995.76501632,    0.        ,    0.        ,    0.        ,\n",
       "           4.23498368,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ],\n",
       "       [   0.        ,  995.21813685,    4.78186315,    0.        ,\n",
       "           0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ],\n",
       "       [   0.        ,  994.88416167,    0.        ,    5.11583833,\n",
       "           0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ],\n",
       "       [   0.        ,  994.8060726 ,    0.        ,    0.        ,\n",
       "           5.1939274 ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ],\n",
       "       [   0.        ,    0.        ,  994.93118302,    5.06881698,\n",
       "           0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ],\n",
       "       [   0.        ,    0.        ,    0.        ,  978.42144602,\n",
       "          21.57855398,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ],\n",
       "       [   0.        ,    0.        ,    0.        , 1040.49629332,\n",
       "           0.        ,    0.        ,    4.88951266,    0.        ,\n",
       "           0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ],\n",
       "       [   0.        ,    0.        ,    0.        , 1063.0920996 ,\n",
       "           0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           1.85549956,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ],\n",
       "       [   0.        ,    0.        ,    0.        ,    0.        ,\n",
       "        1146.67803439,    4.25744534,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ],\n",
       "       [   0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,  995.90592566,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ,    4.09407434,    0.        ,\n",
       "           0.        ,    0.        ],\n",
       "       [   0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,  996.82403603,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ,    0.        ,    3.17596397,\n",
       "           0.        ,    0.        ],\n",
       "       [   0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,  993.89724455,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           6.10275545,    0.        ],\n",
       "       [   0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ,  994.32302015,    5.67697985,\n",
       "           0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ],\n",
       "       [   0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ,  990.90991728,    0.        ,\n",
       "           9.09008272,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ],\n",
       "       [   0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ,    0.        ,    0.        ,\n",
       "         989.63460587,   10.36539413,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ],\n",
       "       [   0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ,    0.        ,    0.        ,\n",
       "         996.97094954,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    3.02905046],\n",
       "       [   0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,  995.59705625,    4.40294375,    0.        ,\n",
       "           0.        ,    0.        ],\n",
       "       [   0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ,    0.        ,  997.74802537,\n",
       "           2.25197463,    0.        ],\n",
       "       [   0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ,    0.        ,    0.        ,\n",
       "           0.        ,    0.        ,    0.        ,    0.        ,\n",
       "         997.68503652,    2.31496348]])"
      ]
     },
     "execution_count": 266,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.imag(Yf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.12 ('event_trigger_mtd')",
   "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.9.12"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "7dac1cb1ce86535b23df7051fbd88cb254e53683203cb38e019b53c60e510804"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}