bvp_triangular.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e09cfced-6a1d-425c-bb38-bab6f6417630",
   "metadata": {},
   "outputs": [],
   "source": [
    "def rmse(n,exa,num):\n",
    "    e = 0.\n",
    "    for i in range(n):\n",
    "        e += (exa[i]-num[i])**2\n",
    "    e = math.sqrt(e/n)\n",
    "    return e\n",
    "\n",
    "def contour_plot(nx,ny,x,sol,exa,unod):\n",
    "    X = np.zeros(nx)\n",
    "    Y = np.zeros(ny)\n",
    "    Z_sol = np.zeros((ny,nx))\n",
    "    Z_exa = np.zeros((ny,nx))\n",
    "    for i in range(nx):\n",
    "        X[i] = x[i][0]\n",
    "    for i in range(ny):\n",
    "        Y[i] = x[nx*i][1]\n",
    "    c = 0\n",
    "    for i in range(ny):\n",
    "        for j in range(nx):\n",
    "            Z_sol[i][j] = sol[c]\n",
    "            Z_exa[i][j] = exa[c]\n",
    "            c += 1\n",
    "            \n",
    "    #plt.contour(X, Y, Z_sol, 20, colors='black')\n",
    "    #plt.contour(X, Y, Z_exa, 20, colors='red')\n",
    "    #plt.quiver(X,Y,unod[:,0],unod[:,1], color='b', units='xy', scale=1)\n",
    "    #plt.show()\n",
    "    #plt.legend(['numerical reults','exact function'])\n",
    "    #plt.legend(['potential','|u|'])\n",
    "    \n",
    "def plotMesh(elements, nodes,l2g=False):\n",
    "    from matplotlib import pyplot as plt\n",
    "    plt.figure(figsize=(9,9)); plt.rcParams.update({'font.size': 22})\n",
    "    for eN, e in enumerate(elements):\n",
    "        element_nodes = nodes[e[0]], nodes[e[1]], nodes[e[2]]\n",
    "        X = [n[0] for n in element_nodes]; Y = [n[1] for n in element_nodes]\n",
    "        cx = sum(X)/3.; cy = sum(Y)/3.\n",
    "        X.append(element_nodes[0][0]); Y.append(element_nodes[0][1])\n",
    "        plt.plot(X, Y); plt.text(cx, cy, \"$T_{{{0:d}}}$\".format(eN))\n",
    "        for n in e:  \n",
    "            plt.text(nodes[n,0], nodes[n,1], \"$N_{{{0:d}}}$\".format(n))\n",
    "    if l2g:\n",
    "        plt.text(nodes[n,0], nodes[n,1], \"$N_{{ {0:d} = elements[{1:d},{2:d}] = elements[0,1]}}$\".format(n,eN,2))\n",
    "    plt.axis('equal')\n",
    "    return plt\n",
    "\n",
    "def mesh(L,h,he):\n",
    "    # element number and size\n",
    "    nex = int(math.ceil(L/he))\n",
    "    ney = int(math.ceil(h/he))\n",
    "    ne = 2*nex*ney\n",
    "    dx = L/nex\n",
    "    dy = h/ney\n",
    "    \n",
    "    # node number\n",
    "    nx = nex+1\n",
    "    ny = ney+1\n",
    "    n = nx*ny\n",
    "    \n",
    "    # mesh\n",
    "    c = 0\n",
    "    x = np.zeros((n,2))\n",
    "    for j in range(ny):\n",
    "        for i in range(nx):\n",
    "            x[c][0] = dx*i\n",
    "            x[c][1] = dy*j\n",
    "            c += 1\n",
    "            \n",
    "    # local-global index for triangle mesh\n",
    "    c = 0\n",
    "    ln = np.zeros((ne,3), dtype=int)\n",
    "    for j in range(ney):\n",
    "        for i in range(nex):\n",
    "            ln[2*c][0] = nx*j+i    # lower tri.\n",
    "            ln[2*c][1] = nx*j+i+1\n",
    "            ln[2*c][2] = nx*(j+1)+i\n",
    "            ln[2*c+1][0] = nx*j+i+1     # upper tri.\n",
    "            ln[2*c+1][1] = nx*(j+1)+1+i\n",
    "            ln[2*c+1][2] = nx*(j+1)+i          \n",
    "            c += 1\n",
    "            \n",
    "    # who is the node with BC\n",
    "    idir = []\n",
    "    ineu = []\n",
    "    for i in range(n):\n",
    "        if x[i][1] == h:   # free surface\n",
    "            idir.append(i)\n",
    "        elif x[i][0] == 0: # left wall\n",
    "            ineu.append(i)\n",
    "        elif x[i][0] == L: # right wall\n",
    "            ineu.append(i)\n",
    "        elif x[i][1] == 0: # bottom\n",
    "            ineu.append(i)\n",
    "    ndir = len(idir)\n",
    "    nneu = len(ineu)\n",
    "    return ne,n,nx,ny,ln,ndir,idir,nneu,ineu,x\n",
    "\n",
    "def local_a(N,invJ,detJ):\n",
    "    a = np.zeros((3,3))\n",
    "    temp = np.dot(invJ,N)\n",
    "    for i in range(3):\n",
    "        for j in range(3):\n",
    "            a[i][j] = np.dot(temp[:,i],temp[:,j])*detJ #*0.5\n",
    "    return a\n",
    "\n",
    "def loc_divN(): # N is dN/dxi\n",
    "    N = np.zeros((2,3))\n",
    "    N1 = np.zeros((2,1))\n",
    "    N2 = np.zeros((2,1))\n",
    "    N3 = np.zeros((2,1))\n",
    "    N1[0][0] = 1.\n",
    "    N2[1][0] = 1.\n",
    "    N3[0][0] = -1.\n",
    "    N3[1][0] = -1.\n",
    "    N[0][0] = 1.\n",
    "    N[1][1] = 1.\n",
    "    N[0][2] = -1.\n",
    "    N[1][2] = -1\n",
    "    return N1,N2,N3,N\n",
    "\n",
    "def jacobian(N,x):\n",
    "    jcb_tr = np.dot(N,x)\n",
    "    invJ = np.linalg.inv(jcb_tr)\n",
    "    detJ = np.linalg.det(jcb_tr)\n",
    "    return invJ,detJ\n",
    "\n",
    "def load(n,ndir,idir,phi):\n",
    "    F = np.zeros(n)\n",
    "    for i in idir:\n",
    "        F[i] = phi[i]\n",
    "    return F\n",
    "    \n",
    "def Dirichlet(n,ndir,idir,L,h,x):\n",
    "    phi = np.zeros(n)\n",
    "    for i in idir:\n",
    "        phi[i] = np.sin(math.pi*(x[i][0]-0.5*L)/L)*np.cosh(math.pi*h/L)\n",
    "    return phi\n",
    "    \n",
    "def assemble(ne,n,ln,ndir,idir,N,x):\n",
    "    Keff = np.zeros((n,n))\n",
    "    for i in range(ne):\n",
    "        # get inverse of JT and determinant of J\n",
    "        [invJ,detJ] = jacobian(N,x[ln[i][:]][:])\n",
    "        \n",
    "        # get local matrices of an element\n",
    "        loc_a = local_a(N,invJ,detJ)\n",
    "        \n",
    "        # assemble local matrices\n",
    "        for j in range(3):\n",
    "            for k in range(3):\n",
    "                Keff[ln[i][j]][ln[i][k]] += loc_a[j][k]\n",
    "        \n",
    "    # swap elements to apply Dirichlet BC\n",
    "    for i in idir:\n",
    "        Keff[i][:] = 0.\n",
    "        Keff[i][i] = 1.               \n",
    "    return Keff\n",
    "\n",
    "#************************************************************************\n",
    "#   CE 7700 project Wen-Huai Tsao 2022 Spring at LSU\n",
    "#   The ALE-FEM for 2D sloshing problem in potential flow\n",
    "#************************************************************************\n",
    "import numpy as np\n",
    "import math\n",
    "import random as rd\n",
    "import sys\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy.sparse import csr_matrix,linalg\n",
    "from scipy import interpolate\n",
    "from scipy import integrate\n",
    "#from mesh import *\n",
    "\n",
    "def main(he):\n",
    "    # physical constant\n",
    "    tho = 998.2\n",
    "    grav = 9.81\n",
    "    \n",
    "    # tank dimensions and 1st-order natural frequency (rad/s)\n",
    "    L = 2.\n",
    "    h = 1.\n",
    "    wn = math.sqrt(np.pi*grav*np.tanh(np.pi*h/L)/L)\n",
    "\n",
    "    # harmonic excitation: x = af*sin(wf*t)\n",
    "    beta = 0.95\n",
    "    af = 0.001\n",
    "    wf = beta*wn\n",
    "    \n",
    "    # mesh size and time interval\n",
    "    #he = 0.05\n",
    "    dt = 0.001\n",
    "    T = 1.\n",
    "    nt = int(T/dt)\n",
    "    \n",
    "    # mesh and l2g index\n",
    "    [ne,n,nx,ny,ln,ndir,idir,nneu,ineu,x] = mesh(L,h,he)\n",
    "    #plotMesh(ln,x)\n",
    "    \n",
    "    # build the derivative of shape function\n",
    "    [N1,N2,N3,N] = loc_divN()\n",
    "    [NN,NS] = loc_NN()\n",
    "    \n",
    "    # time-marching start\n",
    "    for i in range(1):\n",
    "        # prepare BC for weak form\n",
    "        phi = Dirichlet(n,ndir,idir,L,h,x)\n",
    "\n",
    "        # coordinate transformation and assemble local matrices\n",
    "        Keff = assemble(ne,n,ln,ndir,idir,N,x)\n",
    "        \n",
    "        # build load vector that satisfies BC\n",
    "        F = load(n,ndir,idir,phi)\n",
    "\n",
    "        # solve for potential\n",
    "        sol = np.linalg.solve(Keff,F)\n",
    "        \n",
    "        # calculate nodal velocity by averaging cell velocity\n",
    "        #ucel = v_cell(ne,ln,N,x,sol)\n",
    "        #unod = v_c2n(n,ne,ln,ucel)\n",
    "        \n",
    "        # calculate nodal velocity by applying Galerkin method to grad.phi=u\n",
    "        unod = v_Galerkin(n,ne,ln,N,NN,NS,x,sol)\n",
    "        \n",
    "        # calculate the exact solution\n",
    "        [exa_phi,exa_u] = exa_sol(L,n,x)\n",
    "\n",
    "        # calculate the RMS error of potentials\n",
    "        e1 = rmse(n,exa_phi,sol)\n",
    "        e2 = rmse(n,exa_u[:,0],unod[:,0])\n",
    "        e3 = rmse(n,exa_u[:,1],unod[:,1])\n",
    "        \n",
    "        print(he,e1,e2,e3)\n",
    "        \n",
    "        # plot contour of potentials and nodal velocity\n",
    "        #contour_plot(nx,ny,x,sol,exa_phi)\n",
    "        #contour_plot(nx,ny,x,unod[:,0],exa_u[:,0])\n",
    "        #contour_plot(nx,ny,x,unod[:,1],exa_u[:,1])\n",
    "        #contour_plot(nx,ny,x,sol,np.sqrt(unod[:,0]**2+unod[:,1]**2))\n",
    "        #contour_plot(nx,ny,x,sol,exa_phi,unod)\n",
    "        \n",
    "        return e1,e2,e3\n",
    "        \n",
    "        # solve velocity, pressure, etc\n",
    "    \n",
    "        # ALE remesh\n",
    "        \n",
    "def vel_c2n(n,ne,ln,u):\n",
    "    unod = np.zeros((n,2))\n",
    "    c = np.zeros(n)\n",
    "    for i in range(ne):\n",
    "        c[ln[i][0]] += 1.\n",
    "        c[ln[i][1]] += 1.\n",
    "        c[ln[i][2]] += 1.\n",
    "        for j in range(2):\n",
    "            unod[ln[i][0]][j] += u[i][j]\n",
    "            unod[ln[i][1]][j] += u[i][j]\n",
    "            unod[ln[i][2]][j] += u[i][j]\n",
    "    for i in range(n):\n",
    "        unod[i][0] = unod[i][0]/c[i]\n",
    "        unod[i][1] = unod[i][1]/c[i]\n",
    "    return unod\n",
    "\n",
    "def vel(ne,ln,N,x,sol):\n",
    "    u = np.zeros((ne,2))\n",
    "    phi = np.zeros(3)\n",
    "    for i in range(ne):\n",
    "        # get inverse of JT and determinant of J\n",
    "        [invJ,detJ] = jacobian(N,x[ln[i][:]][:])\n",
    "        phi = [sol[ln[i][0]], sol[ln[i][1]], sol[ln[i][2]]]\n",
    "        u[i][:] = np.dot(invJ,np.dot(N,phi))\n",
    "    return u\n",
    "\n",
    "def exa_sol(L,n,x):\n",
    "    exa_phi = np.zeros(n)\n",
    "    exa_u = np.zeros((n,2))\n",
    "    for i in range(n):\n",
    "        exa_phi[i] = np.sin(math.pi*(x[i][0]-0.5*L)/L)*np.cosh(math.pi*(x[i][1])/L)\n",
    "        exa_u[i][0] = (math.pi/L)*np.cos(math.pi*(x[i][0]-0.5*L)/L)*np.cosh(math.pi*(x[i][1])/L)\n",
    "        exa_u[i][1] = (math.pi/L)*np.sin(math.pi*(x[i][0]-0.5*L)/L)*np.sinh(math.pi*(x[i][1])/L)\n",
    "    return exa_phi,exa_u\n",
    "\n",
    "def he_converge():\n",
    "    he = []\n",
    "    e1 = []\n",
    "    e2 = []\n",
    "    e3 = []\n",
    "    for i in range(7):\n",
    "        he.append(1./2**i)\n",
    "        [a,b,c] = main(he[i])\n",
    "        e1.append(a)\n",
    "        e2.append(b)\n",
    "        e3.append(c)\n",
    "    \n",
    "    # converge speed\n",
    "    slope1 = (np.log(e1[-1])-np.log(e1[0]))/(np.log(he[-1])-np.log(he[0]))\n",
    "    slope2 = (np.log(e2[-1])-np.log(e2[0]))/(np.log(he[-1])-np.log(he[0]))\n",
    "        \n",
    "    # plot dx vs error in loglog\n",
    "    plt.loglog(he,e1,'-o',he,e2,'-x',he,e3,'-s')\n",
    "    plt.xlabel(\"he\")\n",
    "    plt.ylabel(\"RMS error\")\n",
    "    plt.legend(['potential','u','v'])\n",
    "    plt.text(0.1,0.1,'slope phi = '+str(slope1))\n",
    "    plt.text(0.1,0.05,'slope vel = '+str(slope2))\n",
    "\n",
    "def v_Galerkin(n,ne,ln,N,NN,NS,x,phi):\n",
    "    unod = np.zeros((n,2))\n",
    "    C = np.zeros((n,n))\n",
    "    D1 = np.zeros((n,n))\n",
    "    D2 = np.zeros((n,n))\n",
    "    for i in range(ne):\n",
    "        # get inverse of JT and determinant of J\n",
    "        [invJ,detJ] = jacobian(N,x[ln[i][:]][:])\n",
    "        \n",
    "        # get local matrices of an element\n",
    "        [d1,d2] = local_d(N,NS,invJ)\n",
    "        \n",
    "        # assemble local matrices\n",
    "        for j in range(3):\n",
    "            for k in range(3):\n",
    "                C[ln[i][j]][ln[i][k]] += NN[j][k]*detJ\n",
    "                D1[ln[i][j]][ln[i][k]] += d1[j][k]*detJ\n",
    "                D2[ln[i][j]][ln[i][k]] += d2[j][k]*detJ\n",
    "    # solve u\n",
    "    \n",
    "    b = np.matmul(D1,phi)\n",
    "    unod[:,0] = np.linalg.solve(C,b)\n",
    "    # solve v\n",
    "    b = np.matmul(D2,phi)\n",
    "    unod[:,1] = np.linalg.solve(C,b)\n",
    "    \n",
    "    return unod\n",
    "\n",
    "def loc_NN():\n",
    "    # NN is local integral of Ni*Nj\n",
    "    NN = np.zeros((3,3))\n",
    "    f = lambda y, x: x**2\n",
    "    s = integrate.dblquad(f, 0, 1, lambda x: 0, lambda x: 1-x)\n",
    "    NN[0][0] = s[0]\n",
    "    f = lambda y, x: x*y\n",
    "    s = integrate.dblquad(f, 0, 1, lambda x: 0, lambda x: 1-x)\n",
    "    NN[0][1] = s[0]\n",
    "    f = lambda y, x: x*(1-x-y)\n",
    "    s = integrate.dblquad(f, 0, 1, lambda x: 0, lambda x: 1-x)\n",
    "    NN[0][2] = s[0]\n",
    "    f = lambda y, x: y**2\n",
    "    s = integrate.dblquad(f, 0, 1, lambda x: 0, lambda x: 1-x)\n",
    "    NN[1][1] = s[0]\n",
    "    f = lambda y, x: y*(1-x-y)\n",
    "    s = integrate.dblquad(f, 0, 1, lambda x: 0, lambda x: 1-x)\n",
    "    NN[1][2] = s[0]\n",
    "    f = lambda y, x: (1-x-y)**2\n",
    "    s = integrate.dblquad(f, 0, 1, lambda x: 0, lambda x: 1-x)\n",
    "    NN[2][2] = s[0]\n",
    "    NN[1][0] = NN[0][1]\n",
    "    NN[2][0] = NN[0][2]\n",
    "    NN[2][1] = NN[1][2]\n",
    "    \n",
    "    # NS is local integral of Ni\n",
    "    NS = np.zeros(3)\n",
    "    NS[0] = 1./6.\n",
    "    NS[1] = 1./6.\n",
    "    NS[2] = 1./6.\n",
    "    \n",
    "    return NN,NS\n",
    "\n",
    "def local_d(N,NS,invJ):\n",
    "    d1 = np.zeros((3,3))\n",
    "    d2 = np.zeros((3,3))\n",
    "    Nxy = np.dot(invJ,N)\n",
    "    for i in range(3):\n",
    "        for j in range(3):\n",
    "            d1[i][j] = Nxy[0][j]*NS[i]\n",
    "            d2[i][j] = Nxy[1][j]*NS[i]\n",
    "    return d1,d2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "61d07b36-b7f5-4bd8-a0ba-bb8c3fa9a903",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.1 0.0016749570146842962 0.03846029275967852 0.058227104694874084\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.0016749570146842962, 0.03846029275967852, 0.058227104694874084)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "main(0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "a38e8d29-730f-420a-aad6-3e35c21fb8bf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0 0.14698716585939656 1.8425137884833362 1.6658755296381904\n",
      "0.5 0.04149373087502142 0.5316200198869317 0.5774800564562319\n",
      "0.25 0.010543911746161112 0.16696243657591237 0.22247702813317224\n",
      "0.125 0.0026226114765448222 0.05466414750036998 0.0810377407411448\n",
      "0.0625 0.0006517812384623231 0.018509631361028483 0.028904864463951664\n",
      "0.03125 0.00016232378046080107 0.006395804701894447 0.010242422268923863\n",
      "0.015625 4.04946176570193e-05 0.0022346998916171344 0.0036230394883368772\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEKCAYAAAAMzhLIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABQV0lEQVR4nO3dd1xWdf/48dcHREFQUVwMEQRcgIji3hNHKmqmdpttG3ZbVja+7bt+d/POsmyYlmWllRVortxa7r0HSwWcKMpe1+f3x9ErkSnrYryfjwcPvc51zue8D17y5pzPeCutNUIIIUR+rCwdgBBCiIpNEoUQQogCSaIQQghRIEkUQgghCiSJQgghRIEkUQghhChQDUsHUJqUUiOAEXXq1Hm4ZcuWlg5HCCEqld27d1/SWje6dbuqivMogoKC9K5duywdhhBCVCpKqd1a66Bbt8ujJyGEEAWSRCGEEKJAkiiEEEIUqEp1ZhckMzOTmJgY0tLSLB1KpWNra4ubmxs2NjaWDkUIYQHVJlHExMRQp04dPDw8UEpZOpxKQ2tNfHw8MTExeHp6WjocIcSt/voIXDuAZ+9/tkVtgtg90POpUjlFtXn0lJaWhpOTkySJ26SUwsnJSe7EhKioXDvAL/cZyQGMP3+5z9heSqrUHcWNeRTe3t75vV++AVUR8n0TouLqu+NV4pvYw6apcD1X0MQepx2vssFzQ6mco0rdUWitl2qtp9SrV8/SoZSKjz76iJSUlGIfv2HDBrZs2WJ+/cUXX/Ddd98VeMzrr7/OBx98UOxzCiHKV3xa/G1tL44qdUdRmkL3xvL+quPEJaTi4mjHjOBWhAS6lmsMH330EZMmTaJ27drFOn7Dhg04ODjQvXt3AB599NHSDE8IYUmXI2HrZ+Vyqip1R1FaQvfG8uJvB4lNSEUDsQmpvPjbQUL3xpao3ejoaFq3bs2//vUv2rRpw5133klKSgpr164lMDAQf39/HnjgAdLT05k1axZxcXH069ePfv36AfDnn3/SrVs3OnTowLhx40hKSgLAw8OD1157jQ4dOuDv78+xY8eIjo7miy++YObMmbRv357NmzfnuFv46quv6NSpEwEBAYwdO7ZEdy5CiHJ0Zgf8NAlmdSBtz/xyOWW1vKN4Y+lhjsRdy/f9vacTyMg25diWmpnNc4sPsHDH6TyPaetSl9dG+BZ67uPHjzNv3jx69OjBAw88wIcffsiXX37J2rVradmyJZMnT+bzzz/nqaee4sMPP2T9+vU0bNiQS5cu8dZbb7FmzRrs7e159913+fDDD3n11VcBaNiwIXv27OGzzz7jgw8+YO7cuTz66KM4ODjw7LPPArB27VpzHGPGjOHhhx8G4OWXX2bevHn8+9//LjR+IYQFmLLh2DLY8gk6ZgdHHBrwu19flmecg8zkMj+93FHk4dYkUdj229GsWTN69OgBwKRJk1i7di2enp7cWMTw3nvvZdOmTbmO27ZtG0eOHKFHjx60b9+eb7/9llOnTpnfHzNmDAAdO3YkOjq60DgOHTpEr1698Pf354cffuDw4cMlvjYhRCnLSIYdX8EnHUlYfC8/ZJ7jztaBTGjkQGhaDH2a9SuXMKrlHUVhv/n3eGcdsQmpuba7Otrx0yPdSnTuW0cQOTo6Eh9feKeT1ppBgwaxcOHCPN+vVasWANbW1mRlZRXa3n333UdoaCgBAQHMnz+fDRs2FB68EKJ8JJ6HHXPI3jWP7aTyW6NmrKvfnEydja9DU14JfJwhnkOoW7MuW+O25tlx7WTrVGrhVMtEUZgZwa148beDpGZmm7fZ2VgzI7hVids+ffo0W7dupVu3bvz4448EBQXx5ZdfEh4ejre3NwsWLKBPnz4A1KlTh8TERBo2bEjXrl2ZOnWqeb/k5GRiY2MpaDn1OnXqcO1a3o/YEhMTcXZ2JjMzkx9++AFX1/LtqBdC5OHCMdj6KTGHFxNWuyahTRtyjkzq1bRlvNcIQrxDaNUg58+hDeM3lHlYkijycGN0U1mMemrVqhWzZ8/mgQceoG3btsyaNYuuXbsybtw4srKy6NSpk3l00pQpUxgyZAguLi6sX7+e+fPnM3HiRNLT0wF46623CkwUI0aM4M477yQsLIxPPvkkx3tvvvkmXbp0oVGjRnTp0oXExMQSX5sQohi0hqhNpG2ZxdqzW/i9bl22uzZCoeju0plnfUbTr1k/alrXtFiIFb4ehVLKHvgMyAA2aK1/KOyYvOpRHD16lDZt2pRNkEUUHR3NHXfcwaFDhywaR3FUhO+fEFVKdib60G8c3fYRv2WcZbmDA4lWCtfaTQlpOZYQ7xCa2jct15Dyq0dhkTsKpdTXwB3ABa21303bhwAfA9bAXK31O8AYYLHWeqlS6ieg0EQhhBAVVtpVEnbOYdnB+fxeI4vjtjWpaefIQPeBjGk1jk5NO2GlKtY4I0s9epoPfAqYpwkrpayB2cAgIAbYqZRaArgBB6/vlk0l5uHhUSnvJoQQJWdKOM22TW/ye+wm1trZkGlvRVv7FrzkO5mhLYZRr1bFXVHCIolCa71JKeVxy+bOQLjWOhJAKbUIGIWRNNyAfchwXiFEJRMbsZrQ7f8jLPUUZ2vUoF5tO+5qNoDR7afk6piuqCpSZ7YrcOam1zFAF2AW8KlSajiwNL+DlVJTgCkA7u7uZRimEEIULD0zlbU7P+b3E7+yHWOofbfazjzd7kH6tb6TWta1LBzh7alIiSJPWutk4P4i7DcHmANGZ3ZZxyWEELc6en4/v+34H8vi95KowNWkeaxRR0K6/R/ODSvH3UNeKlKiiAWa3fTa7fq2IitsmXEhhChtV9Ovsuz4L/x++DuOZVyhpkkzwFSTMa3uonOXp7CysbV0iCVWkRLFTsBHKeWJkSAmAHffTgNa66XA0qCgoIfLID4hhADApE1sO7uN0MPfszbubzIw0SY9g/+zdWNY92eo13IoVKE6LpYaHrsQ6As0VErFAK9precppZ4AVmEMj/1aa31bCxCV2h1FOZQWFEJUPnFJcYSFhxF6/Cfi0uKpm53NnSlpjHbuReshz0OTtpYOsUxYatTTxHy2LweWl6Dd0rmjuFFacNx8I1ncKC04bn6Jmr11wt0HH3xAUlISr7/+eonaFUKUnfTsdNadXsfvJ39j29ntoDVd01J5Kk3T328Stbo8BnXKd2JceatIj55KrMh3FCtegHMHC96njjMsGG38mXgWGrWGDe8aX3lp6g9D3ylW3EKIiudo/FF+D/+dZZHLuJZxDRcTPHY1gVFWjrh0eQ4C/wU17S0dZrmoUomiVPsobB2NJHH1DNRrZrwWQlQZfX/qm+eqq/Y29rjXcefo5aPUxIoBqemMTrhMl4b+WA1+C1rfAVbWFojYcqpUoiiyovzmf+NxU+/nYNc86Pt8zj6LYqhRowYm0z81LdLS0krUnhCi+PKrKZ2cmYy+epoXL19jeGKi0TE9bBq4dynnCCuOKpUoSq0z++Y+Cc/e4Nkr5+tiatKkCRcuXCA+Ph4HBwf++OMPhgwZUrJYhRC3rbDFUH+JjjYeLXV9HJy8yieoCqxKJYpSe/QUuydnUvDsbbyO3VOiRGFjY8Orr75K586dcXV1pXXr1iUKUwhxe6KvRrMiegUrolYUvOP0w2BfeoV/KrsqlShKTV5DYD17l/jRE8C0adOYNm1aidsRQhTNueRzrIxayfKo5Ry9fBSFoqOTH1EFHSRJIocqlShkZrYQAuBy2mX+jP6TFVEr2HNhDwB+Tn4822IswbFHabp7Jf4eUtWxqKpUopCZ2UJUX4kZiaw9vZYVUSvYfnY72Tobr3pePOE/haGp6bjvWwwXl4Ndfeg2FaeLq4nPSs7VjpNNHQtEX7FVqUQhhKheUrNS2RizkRWRK9gcu5lMUyauDq7c73c/Qx198Tn2J2rVe5B+DZwDYNRs8BsLNnZs+Msp/xUYRA6SKIQQlUpmdiZb4rawInoF60+vJyUrhYZ2DRnfajxDmwfjfzkWtfMriHwVrGzAdzR0ngJuQTnXXyrDvsiqRhKFEKLCyzZls+v8LlZErWDN6TVcTb9K3Zp1Geo5lKGeQwmq2wLrvT/Aj/fA1dNQxwX6vQwd7wWHxpYOv9KrUolCOrOFqDq01hy8dJAVUStYFb2Ki6kXsathR3/3/gz1GEp3l+7YnDsI2+fCwcWQnQ4evSD4LWg1HKyr1I83i6pS30npzBai8jtx5QQrooy5DrFJsdhY2dDLtRdDWwylj1sf7LCCw6GwcgjE7gIbewicBJ0eqrKrt1palUoUpSW/NWCcbJ3YMH5D+QckRBV3+tppVkStYGX0SsITwrFW1nRx7sKjAY8ywH0AdWrWgasxsPED2D0fUi6BkzcMeRfaTwTbepa+hCpNEkUe8lsDJr/tRfHCCy/QrFkzpk6dCsDrr7+Og4MDzz77bLHbFKIyO598npXRK1kZtZJD8cbS+x0ad+ClLi8xqPkgnOycQGuI3gw75sCx5aBN0GoodH4YPPuClZVFr6G6qJaJ4t0d73Ls8rFiHXv/yrzLd7du0JrnOz+f73Hjx4/nqaeeMieKn3/+mVWrVhUrBiEqqytpV1h9ajUrolaw+/xuNJo2DdrwTMdnCPYIxtnB2dgxPRF2fAU758LFY2DXALr/G4IegPrNLXsR1VCVShQVuTM7MDCQCxcuEBcXx8WLF6lfvz7NmjUr/EAhKrmkjCTWn1nP8qjlbIvbRpbOwqOuB48FPMZQz6F41PP4Z+eLJ2DnV7BvIWQkgnN7GPUZ+I0BGztLXUK1V6USRVE7swv6zR/A/1v/fN/7Zsg3xYoNYNy4cSxevJhz584xfvz4YrcjREWXlpXG5tjNrIhawaaYTaRnp+Ns78w9vvcwzHMYreq3Qt2Y02DKhhMrjcdLkRvAuuY/cx9cO1ap2tOVVZVKFBXd+PHjefjhh7l06RIbN260dDhC3LaCBnqsHreabXHbWBG1gnVn1pGcmUwD2waM8RnDMM9htGvUDit1U59Ccjzs+RZ2fW0UCKvrCv1fgQ73gkOjcrwqURhJFHlwsnXK9z9DSfj6+pKYmIirqyvOzs4laksISyhooEf/n/uTkJ5AHZs6DG4+mKGeQ+nUtBM1rG75MRO7G3bMhUO/3jT34b/QapjMfaig5F8lD2U5BPbgwUJqdQtRSXVz7sZQz6H0cO1BTeuaOd/MTIMjocbjpdjdUNMBOtxjzH1o3MYi8Yqik0QhhCiSc8nnCnz/vT7v5d6YcMZ4tLTnW0iJBycfGPo+BEwA27plFKkobZIohBD5SslMYe3ptSyJWML2s9uLdpDWELXRGN56fLmxreX1uQ8t+krndCVUpRJFRR4eK0RlkW3KZuf5nSyNWMrqU6tJzUrF1cGVRwMe5fP9n+d/YHoi7F9kJIhLx425Dz2eNOY+OLqX3wWIUlelEkVhw2O11v8MyRNFVlghelE1RCREsCRiCcsil3E+5Tx1bOowzHMYI71GEtg4EKUUPx+aT3x2aq5jnbQV/K+NMffBJRBCvjCGuNrYWuBKRGmrUomiILa2tsTHx+Pk5CTJ4jZorYmPj8fWVv7DV0WX0y6zImoFSyKWcCT+CNbKmh6uPXi207P0deuLbY2c/+4benwAv9wHY+dBRhJsfA/OHQCrGuB35/W6Dx0tczGizFSbROHm5kZMTAwXL160dCiVjq2tLW5ubpYOQ5SS9Ox0Np7ZyNKIpfwV+xdZOos2DdrwXKfnGOo5lIZ2DfM/2MkbfILh+zHGukvKCgLvgQGvydwHCwrdG8v7q44Tl5CKi6MdM4JbERJYejXBq02isLGxwdPT09JhCGERWmv2X9zPkoglrIxeSWJGIo3tGnOP7z2MaDECn/o++R9sMkHkemP00vEVoLOhvgdciYaeT8OAV8rrMkQeQvfG8uJvB0nNzAYgNiGVF38zhuGXVrKoNolCiOroTOIZ/oj8g6URSzmTeAa7GnYMcB/ACK8RdGnaBWsr6/wPTroI+743lvW+Eg21naD7E9DYF1a9CL2fg13zoEUfKR9qQe+vOm5OEjekZmbz/qrjpZYoZI1eIUpB37592bVrV7mfd8OGDdxxxx05tiVmJPLriV/xCfah/yf9+Xzf57jYu/BWj7dYf9d63u71Nt1duuedJLSG6L9g8QPwYRvil7xKvy9icHgvnSeig8F7oJEkxs2H/i8Zf/5yH/uXf023bt3w9/dnxIgRXLt2DYAffviB9u3bm7+srKzYt28fALt378bf3x9vb2+mTZtmHjSxf//+PNu64fTp0zg4OPDBBx+Yt3l4eODv70/79u0JCgrKdVn/+9//UEpx6dKl65epmTZtGt7e3rRr1449e/aY9x0yZAiOjo65vq+ffvop3t7eOdq5cY3t2rXD39+f7t27s3//fvN7K1eupFWrVnh7e/POO++Yt69du5YOHTrQvn17evbsSXh4OADTp083f69atmyJo6NjjhiuXbuGm5sbTzzxBABpmdnEJuQeXAAQl8/2YtFaV7mvjh07aiHKU58+ffTOnTvL/bzr16/Xw4cP1xnZGXrjmY36mQ3P6A7fddB+8/30iN9H6K8OfKXjEuMKbyjlstZbP9P6kyCtX6ur9dvNtF7+vE6K2qM3b96sP//8cz116lStN8/UOnJjzmMjN+qg1s30hg0btNZaz5s3T7/88su5TnHgwAHdokUL8+tOnTrprVu3apPJpIcMGaKXL1+utdY6KCiowLbGjh2r77zzTv3++++btzVv3lxfvHgxz0s7ffq0Hjx4sHZ3dzfvs2zZMj1kyBBtMpn01q1bdefOnc37r1mzRi9ZskQPHz48Rzt79uzRUVFRuc71999/68uXL2uttV6+fLm5raysLN2iRQsdERGh09PTdbt27fThw4e11lr7+PjoI0eOaK21nj17tr733ntzxT1r1ix9//3359g2bdo0PXHiRP3QlEf1x2tO6A7/+VM3f/6PPL+6v702z+9HQYBdOo+fqXJHIcRtSE5OZvjw4QQEBODn58dPP/2Ua5+FCxfi7++Pn58fzz//z0rFDg4OTJ8+HV9fXwYMGGAeWBEREcGQIUPo2LEjvXr14tix3LVSXn/9de655x66deuGj48PX331FVprTl09xZGzR3Dp5sKgToP47qXvGOMzhoXDF3Lto2u0T2//T42HW2kNZ3bC74/B/1rDyhegVl1jWe+nj8HQd7D3CKRnz57/jHrr+VTux0yevTkRd5XevY3tgwYN4tdff83z+zJhwgQAzp49y7Vr1+jatStKKSZPnkxoaCgAJ06cyLet0NBQPD098fX1zfua8jB9+nTee++9HKMdw8LCmDx5MkopunbtSkJCAmfPngVgwIAB1KlTJ1c7gYGBeHh45NrevXt36tevD0DXrl2JiYkBYMeOHXh7e9OiRQtq1qzJhAkTCAsLA0ApZb5Tunr1Ki4uLnl+vyZOnGh+vXv3bqLOxJHSyJfFe2L5cPUJApo58kQ/L+xscv4ot7OxZkZwqyJ/jwojfRRC3IaVK1fi4uLCsmXLAOM/+c3i4uJ4/vnn2b17N/Xr12fw4MGEhoYSEhJCcnIyQUFBzJw5k//85z+88cYbfPrpp0yZMoUvvvgCHx8ftm/fzuOPP866detynfvAgQNs27aN6IvR9Ojcg19r/Mrx48c5ffQ098+/n4mdJ/LShJfok9UHv4Z++V7D9GlTWb9yKaRehsxUY+SSnSMTJtzHCw9/XKzvi6+vL2FhYYSEhPDLL79w5syZXPv89NNP5h+UsbGxOUbSubm5ERsbW2BbSUlJvPvuu6xevTrHYycwfvAOHjwYpRSPPPIIU6ZMAYyE4OrqSkBAQI79Y2Njc9SDuXH+ki7WOW/ePIYOHZrvObZvN2a3z507l2HDhmFnZ0fdunXZtm1bjnZOnTpFVFQU/fv3B+DAmSuMvOcRavSfRlrkXpo71eanp3rTqqmR0Lwb16neo56UUi2Al4B6Wus7LR2PqN78/f155plneP7557njjjvo1atXjvd37txJ3759adTIGCr6r3/9i02bNhESEoKVlZW5DsmkSZMYM2YMSUlJbNmyhXHjxpnbSE9Pz3XezOxM2vRuw1N/PcW2s9vQXprUyFQmt53Mpm6bmHvXXAAWBy4mOjqanj175g7+7H7Y9TUzm/wOdydDk/bQ6QHwHwe1cv8GfTu+/vprpk2bxptvvsnIkSOpWTPnooDbt2+ndu3a+Pnln8AKa+v1119n+vTpODg45Drmr7/+wtXVlQsXLjBo0CBat25NUFAQ//3vf/nzzz9LdG1FtX79eubNm8dff/1V6L4zZ85k+fLldOnShffff5+nn36auXPnmt9ftGgRY8eO5e+Iy3y5KYIVP82npkt7pg7rRO3oNMKPZJqTBBijm0ozMdyqTBOFUupr4A7ggtba76btQ4CPAWtgrtb6nXyaQGsdCTyolFpclrEKURQtW7Zkz549LF++nJdffpkBAwbw6quvFqstpRQmkwlHR0dzB+/NTNrErnO7WBKxhPmH55NpyqTDtQ487P8w693XMzloMvXq1WO33W7zMdbW1mRlZf3TSGYa7P3eGNoauxtq2DJ9R1PWRziATRJ8OQuYBcCECRN44YUXinUtrVu3Nv9APnHihPmO64ZFixbleIzi6upqfkQDEBMTg6ura4Ftbd++ncWLF/Pcc8+RkJCAlZUVtra2PPHEE+ZjGzduzOjRo9mxYwf169cnKirKfDcRExNDhw4d2LFjB66urjnuem4+f3EcOHCAhx56iBUrVuDk5GS+xrzOcfHiRfbv30+XLl0Ao07NkCFDzPtlZpv48psFNB7yOH98vYPGdWrhzVnOHt/LnMfXkJSUREZGBg4ODjk6yMtUXh0XpfUF9AY6AIdu2mYNRAAtgJrAfqAt4A/8cctX45uOW1zU80pntigrsbGxOjU1VWut9dKlS/WoUaO01v90ZsfFxZk7TbOysvSAAQN0aGio1lprQC9cuFBrrfWbb76pn3jiCa211t26ddM///yz1lprk8mk/9j8h/5498d64C8Dtd98P93lhy669wO9tY+vj05OSdaXLl3SzZo107GxsebO7BumTp2qv/nmG63PH9V9/JvpnY83MTqnP+mk9dbPjU7rYvjmm2+Mzux8nD9/XmutdXZ2tr7nnnv0vHnzzO9lZ2drFxcXHRERkeOYWzuzly1bVmhbN7z22mvmzuykpCR97do189+7deumV6xYkeuYmzuh//jjjxyd2Z06dcqx763f1/za0VrrU6dOaS8vL/3333/n2C8zM1N7enrqyMhIc2f2oUOHdGZmpnZyctLHjx/XWms9d+5cPWbMGJ2UlqnnbY7U7ad/o63rNtb9P1ivf9p5WqdlZuVot7B/i5Ign87sMr2j0FpvUkp53LK5MxCujTsFlFKLgFFa67cx7j6KRSk1BZgC4O4uC5CJsnHw4EFmzJiBlZUVNjY2fP55zkXynJ2deeedd+jXrx9aa4YPH86oUaMAsLe3Z8eOHbz11ls0btzY3BH+ww8/8NAjD/HsK8+SkJKAbZAtTUOa0t2lO093fJq+zfry7ol3icyMZED/AVy6dIlXXnkFFxcXTpw48c/Js9LhUjhs2gjRpyExBdyHwfjnoHn3Yq/a6uHhwbVr18jIyCA0NJQ///yTtm3b8tBDD/Hoo48SFBTEwoULmT17NgBjxozh/vvvNx+/adMmmjVrRosWLXK0+9lnn3HfffeRmprK0KFDzc/2C2orL+fPn2f06NHGtyAri7vvvjvHb+h5GTZsGMuXL8fb25vatWvzzTf/lDi+MaAgKSkJNzc35s2bR3BwMLNmzeK9997j3LlztGvXjmHDhjF37lz+85//EB8fz+OPPw5AjRo12LVrFzVq1ODTTz8lODiY7OxsHnjgAXMn/FdffcXYsWOxsrLCoW49ej/4Kt3fWcfV1EzsT2zmzrvu4sfpfbCyqhjLDSldxgu+XU8Uf+jrj56UUncCQ7TWD11/fQ/QRWv9RD7HOwH/DxiE8Zjq7cLOGRQUpC0xpl2IvNwoH3rkkSO0/bKtebuTrRMvd32ZpRFL2RS7iSxTFq3qt2KE1wiGtxieYymN119/HQcHB5599tncJ4iPMCbF7fvBqPlQ3xM63gft/yXLalRgkReT+GpzFL/uiSEz20Rw26ZM6dOCDu71LRaTUmq31jrXRJQK35mttY4HHi3KvrLMuKiICiofOn3DdBraNWRSm0nc0eIOWjUo4pDG7Eyj1sOub4zlNZQ1tB5mLOnt2ResZOR7RbXn9BXmbIxk1ZFz2FhbMbaDGw/38qRFo9yd9BWFJRJFLNDsptdu17eVmC5kmXEhLOnmu4kbvhj4BV2cu+SuK32L119/3fhLwhmjWtye7yDpPNR1hX4vGQvz1ZU67BWVyaRZd+wCczZFsiP6MnVtazC1rzf3dvegUZ1alg6vUJZIFDsBH6WUJ0aCmADcXRoNyx2FqEhSs1JZe3ptgfv0cO1ReEOmbAhfY4xcOvmnMVHOZzAE3Q/eg8C6wj8YqLbSs7IJ2xfHnE2RhF9IwtXRjlfvaMv4Ts2wr1V5/t3KenjsQqAv0FApFQO8prWep5R6AliFMQLqa6314dI4n9xRCEvTWrP3wl7CIsJYFb2K5Mzk4jeWeA72LoDd38LVM2Df2FitteO9UjGugruWlsmP20/z9V9RXEhMp41zXT4a357h7Zyxsa58jwXLNGKt9USttbPW2kZr7aa1nnd9+3KtdUuttZfW+v+VZQxClIfuvbrzfz/+H3f8fgf3rryXFVErGOg+kK+Dv769hkwmiFgPP90DM31h3Vvg5AXjvoWnjxhLehchScyfP9+8cFxJvPTSSzRr1izPSW43O3DgAN26dcPX1xd/f3/S0tKKdPyvv/6KUsq8oGJ0dDR2dnbmhfEeffSf7smMjAymTJlCy5Ytad26tXlpjy+++MK8IGDPnj05cuQIAPHx8fTr1w8HB4dc34v8FiS8fPkygwYNwsfHh0GDBnHlyhXAmOHdrl0786KDN0+qe/755/Hz88PPz48vv1nAf5cfpfvb63jm8Qc59smDWP32LE33fc1wv8bYWFvx/vvvm6/Pz88Pa2trLl++DBgT8Xx9ffHz82PixInm7+MN06ZNy/G9nD9/Po0aNTK3d/OkvVKV15jZyvoFjADmeHt7397gYSGKITkjWS8JX6IfXPWgtm9lr71e89IPrHxAh4WH6eSMZPN+fRb10X7z/XJ99VnU55/Gki5p/dfHWn/c3pj38I6H1qte0vpSeLFiK62x9lu3btVxcXHa3t4+330yMzO1v7+/3rdvn9Za60uXLumsrKxCj7927Zru1auX7tKli3lBxaioKO3r65vneV599VX90ksvaa2NORY35jJcvXrVvE9YWJgODg7WWhtzKnIsaHiT/BYknDFjhn777be11lq//fbb+rnnntNaa52YmKhNJpPWWuv9+/frVq1aaa2N+RgDBw7Uh2Mu62nfbdW1nH108+k/6yd+3KNnf/uzNplM2mQy6QkTJujPPvss1zUtWbJE9+vXT2utdUxMjPbw8NApKSlaa63HjRtnzIm5bufOnXrSpEk5vpelPaeC6rAooNZ6qdZ6Sr169SwdiqiikpKS6DmwJ019mtKwRUOm/m8qMYkxuDi4MKv/LOYFzyN5RzJdOnQxLwq4wTWEg71nEzU1ioH7BmJ6z0TjL2vyS+2BcGoLEbPvYkh7ZzqOm06vT6M5Fvg6PH0UBht3EyaTCQ8PDxISEsxx+Pj4cP78eS5evMjYsWPp1KkTnTp14u+//y7V6+3atWuh6x/9+eeftGvXzjwD2snJCWtr60KPf+WVV3j++eeLXGb366+/5sUXXwTAysqKhg2N4cN169Y175OcnGxe/M/e3j7ngobXFbQgYVhYGPfeey8A9957r3m7g4ODud0b59Bas3LzTuLreDHsky38eSKBNr5+POeXxicTA3l88jiUUiil6Ny5c46Z6DfcuvBfVlYWqampZGVlkZKSYl4sMDs7mxkzZvDee+8V6XtV2qpUohCirMQlxfHF/i/o/UZvjmQewf11d55Z9AzfPfUdy8csx8XBhca1G5sXBVy3bh379u1j586dhB5Ohl/uMy8KePjX9+ljd4I33p0J3wxlyodL+OS5e9m97xAfLFjJ4x8vBZt/frhZWVkxatQofv/9d8BYyqJ58+Y0adKEJ598kunTp7Nz505+/fVXHnrooQKvY/369TnqQ9z46t69e7G/NydOnEApRXBwMB06dCjSD7M9e/Zw5swZhg8fnuu9qKgoAgMD6dOnD5s3bwYwJ8lXXnmFDh06MG7cOM6fP28+Zvbs2Xh5efHcc88xa9asAs9d0IKE58+fNye2pk2b5jjH77//TuvWrRk+fDgPv/QeIZ9t4ecoa07s2swTvZqx9KF2XI3YR+qVCznOl5mZyYIFC3JNAkxJSWHlypWMHTsWMJb8ePbZZ3F3d8fZ2Zl69eoxePBgwKiFMXLkyDyT7q+//kq7du24884781yMsTRUnm73IpBRT6I0pWSmsPb0WsLCw9h+zlj1s3Xr1sQujKXn3p4Mdh9Mr5ZFWBTw8GFCnvgGqxf6Mj7hM1h0iEl+MOa3dJIGzmTLu9MYN3snzDZ+s8xrUcDx48fzn//8h/vvv59FixaZFxdcs2aN+Zk8GIVtkpKS8r2mfv365bmuVElkZWXx119/sXPnTmrXrs2AAQPo2LEjAwYMyHN/k8nE008/zfz583O95+zszOnTp3FycmL37t2EhIRw+PBhsrKyiImJoXv37nz44Yd8+OGHPPvssyxYsACAqVOnMnXqVH788Ufeeustvv322xJf1427gRuG3jGSZOcOfPDt77z08it0njqT/02fTIQvfP/SZFY1akS3bt3Md1M3PP744/Tu3TvXApJLly6lR48eNGjQAIArV64QFhZGVFQUjo6OjBs3ju+//57+/fvzyy+/sGHDhlwxjhgxgokTJ1KrVi2+/PJL7r333jxXHi6pKpUotIx6EiWktWb3+d0siVjCquhVpGSl4ObgxtT2UxnhNQJXB1cuj75c9EUBM5Lh7AFYudV4ff4guARA8HOoTa9iajceR8dXC/3h3a1bN8LDw7l48SKhoaG8/PLLgPFDd9u2bUV+fLN+/XqmT5+ea3vt2rXZsmVLkdq4lZubG7179zY/Cho2bBh79uzJN1EkJiZy6NAh+vbtC8C5c+cYOXIkS5YsISgoiFq1jHkFHTt2xMvLixMnTtCxY0dq167NmDFjABg3bhzz5s3L1faECRN47LHHCoy3oAUJmzRpwtmzZ3F2dubs2bM0btyYK8kZLNh2im+3RBOfnEGAbxDnls7k53v9aNK4EXR9hddeNeqG33333bRs2dLc9htvvMHFixf58ssvc8Vx60KJa9aswdPT0/xLxpgxY9iyZQv169cnPDycG78Ap6Sk4O3tTXh4uHkBQoCHHnqI5557rsBrLy559CQEEJsUy+f7P2fYb8O4f9X9rIpeRbBHMPOHzGf5mOU8GvAorg6uxMXFUbt2bSZNmsSMGTNylNAE6Ny5Mxs3buTS7qVk//wACz94lj56C2RnYNKwOGsgXIvlx8Vh9OzZk7p16+Lp6ckvv/wCGInq5lKaNyilGD16NE8//TRt2rQx/4AYPHgwn3zyiXm/whLOjTuKW7+KmyQAgoODOXjwICkpKWRlZbFx40bats09ufCGevXqcenSJaKjo4mOjqZr167mJHHx4kWys436z5GRkZw8eZIWLVqglGLEiBHm36rXrl1rPsfJkyfNbS9btgwfH58C43V2djbXgNBa891335nX4xo5cqT5buTjz7+ioW8Pur+zjnd/3oC/a10WTenKq93ssNZZNG7UkOzsbOLjjZn3Bw4c4MCBA+bHRXPnzmXVqlUsXLgQq1tmyl+9epWNGzeazwvGGnXbtm0jJSUFrTVr166lTZs2DB8+nHPnzpm/X7Vr1zaXTr1RbAlgyZIltGnTpsBrL7a8ergr6xcy6knchuSMZB0WHqYfWPmAeSTSgysf1EvCl+QYtXSzlStXan9/fx0QEKCDgoLMo3X69Omjd25cpfXmmfrHyc21X2Mr7du4hn5ubCetd36j9bue2r62nZ4+fbr2bemp+3nZ6gs7w7TWWkdGRurg4GDdrl073aZNG/3GG2/kee6dO3dqQM+fP9+87eLFi/quu+7S/v7+uk2bNvqRRx7RWpfeaJgZM2ZoV1dXrZTSrq6u+rXXXtNaG6OLXnnlFfN+CxYs0G3bttW+vr56xowZhR5/s5vLyC5evFi3bdtWBwQE6MDAQL1kyRLzftHR0bpXr17a399f9+/fX586dUprbZQHvXFM37599aFDh8zHNG/eXNevX1/b29trV1dXcynSnTt3al9fX93Erblu2nWUbv7cUt397bX623UHdZcevXW9Js20nUeA9py+SD/90z79zEtvmM/RtWtXvXnzZq211qmpqbpNmza6TZs2ukuXLnrv3r3mc1tbW+sWLVrogIAAHRAQkOPf9ZtvvtHjx4/P9b149dVXdatWrbSvr6+eNGmSTktLy7XPzaOeXnjhBd22bVvdrl073bdvX3306NHc/4i3gXxGPZX5ooCWIIsCivyYtInd53cTFh7G6lOrSclKoVmdZozyGsUIrxG4OOQuSVlwg9nGvIc98+H4CjBlgXs36HAvtB0FNWvDXx+Bawcc/If9038QtQli9xilRYVFhO6N5cXfDpKamW3eZqXApMGhVg3u7uLO/T08cK5nZ8Eoy1elXRRQiNIQkxjD0oilhEWEEZsUi72NPUM8hzDKaxSBjQNzdFoWydUYoyDQ3u+NWdO1naDLo9BhMjS6ZWG/vJKBZ+/ctadFuXp/1fEcSQKMJFHXtgZ/vdCfurY2Foqs4pFEIaqslMwUVp9aTVhEGDvP7USh6OzcmantpzLAfQC1bWrfXoPZmXBipbEgX/ga0CZo0Q8G/QdaD4caBS/uVtBoJFG+srJNxCak5vleYlqWJIlbVKlEIcNjxc2Plv489SepWam413Hn34H/ZkSLETg7FGOF1cuRRnLY96OxYmsdZ2PNpQ73QH2PUr8GUXaysk0s2R/Hp+vC893HxbH6PGoqqiqVKLQMj622ziSeYWnEUpZELDE/WhrmOYxR3qNo36j97T9aykyDY38YBYGiN4OyAp9gY0E+WbG10snKNhG2L45P14cTdSmZNs51ub+HB4t2nCY102Tez87GmhnBRawJUo3Ip11UWimZKfx56k/CwsPYdX4XCkUX5y48EfgEA9wHYFejGL8ZXjhq3D3sXwipV4wF+Pq/bFSLq3ubHd3C4rKyTYTui+PTdSeJjk+hrXNdvrynI4PaNMHKShHg5sj7q44Tl5CKi6MdM4JbERLoaumwKxxJFKLCulFC9FZ1a9alb7O+rD61unQeLWUkw+HfjeW8Y3aAlY3R59DxXqkWV0llZZv4fW8sn64P59T1BDHnno4Matskx91lSKCrJIYikEQhKqz8Sohey7jG2tNrGeY5jBDvEAIaBdz+oyWAuL1Gcji4GDISwcnHWIgvYCLYNyz8eFHhZN5IEOvCOX05BV+Xunw1OYiBbRoX7zMigEIShVLKGjistW5dTvGUiHRmVx/r71pfvEdLaVfhwM/G46VzB6CGLfiONuY9uHcF+WFSKWVmm/htTwyfrg/nzOVU/F3rMXdyEAMkQZSKAhOF1jpbKXVcKeWutT5dXkEVl3RmV343Ri2FhocWuN9tJQmt4cx24+7h8O+QlQpN/GHYB+A/DuwcSxSzsJzMbBO/7jYSRMyVVNq51eP1Eb70by0JojQV5dFTfeCwUmoHYK7rqLUeWWZRiWonNimWJRFLCAv/Z0JciSXHG53Se76DS8ehpgMEjDfuHlwC5e6hEsvIMvHrnhhmX08QAW71+M8oX/q1kgRRFoqSKF4p8yhEtZSalcqaU2tyLOPdpWkX84S4Lj92uf1GTSaI2mgkh2N/QHYGuHWCkZ8aj5hqFVzSU1RsGVkmFu82EkRsQioBzRx5c5QffVs1kgRRhgpNFFrrjUqpJkCn65t2aK0vFHSMEPnRWrPv4j5Cw0NZFb2K5Mxk8zLeI71G5lhrycnWKc8ObSdbp1zbuHYW9v0AexfAlWiwdYSgB40lNZrkv5KpqBwyskz8svsMn62PIDYhlfbNHPl/o/3o01ISRHkoNFEope4C3gc2AAr4RCk1Q2u9uIxjE1XIueRz5kdLpxNPY1fDjsHNBxPiHUKHJh2wUrmHoG4Yv6HgRrOzjKU09nwLJ1aBzgaPXtDvZWgzIkeVOFE5pWdl88uuGD5bH07c1TQC3R357xh/evs0lARRjory6OkloNONuwilVCNgDSCJQhQoLSuNdafXERoeyraz29BogpoE8XC7hxncfHDhay1dX3U1x+J5UZsgfC1Y28DeHyAxDuwbQ/d/G3cPTl5lek2ifKRnZfPzzjN8tiGCs1fT6ODuyDtj29FLEoRFFCVRWN3yqCkeKXgk8qG15uClg4SGh7IyaiWJmYm42LvwSMAjjGwxkmZ1mxW9MdcO8Mt9MG4+NOsKm/9nfJkyAQXeA2HYe9ByiJE4RKWXlpnNz7vO8Pn1BNGxeX3eu7MdPb0lQVhSURLFSqXUKmDh9dfjgeVlF1LxyTwKy7mQcoE/Iv8gLDyMyKuR2FrbMrD5QEK8Q+jUtFOej5YK5dkbBr0FP4wDlDGstXZD6PQQBE4Cx9tIOqJCS8vM5qedRoI4dy2NoOb1ef/OAHp4O0mCqAAKLFykjH8hN4yO7J7XN2/WWv9eDrEVmxQuKh8Z2RmsP7OesPAw/o77G5M2Edg4kFFeowj2CMahZjFHGKUnGfMd9i4w5j+gAA1+d8KYOWBlXVgLopJIy8xm0Y7TfL4xgvPX0unkUZ+nBraku5ckCEsoVuEirbVWSi3XWvsDv5VZdKLS0FpzJP4IoeGhLI9azrWMazSu3ZgH/R5kpNdIPOp5FLdhOLMD9n4Hh36HzGRjSY2O98GRMOj0MOyaB6f+loI/VUBaZjYLd5zm8w0RXEhMp7NnA2be1Z5ukiAqpKI8etqjlOqktd5Z5tGICutS6iWWRS4jNDyU8IRwalrVZEDzAYR4hdDFuQvWxf0tP+miMSlu7/fGpDgbe/AbDYGTISsNFt8Pd313vSJcr3/6LCRZVEppmdn8uN24g7iYmE4XzwZ8PCGQbl55DHkWFUahNbOVUscAb+AUxsxshXGz0a7swyseefRUOjKzM9kUs4nQ8FA2x24mW2fTrmE7RnmPYojnEOrWrFu8hrOzIGKtMSnuxEqjzrRbZ6MQkO9oqFXH2C+/UU9Sa7rSSc3I5oftp/hyUyQXE9Pp2qIBTw5oKQmigsnv0VNR+ih6YSSJHLTWubZVFJIoSubY5WOEhYexLHIZV9Kv0NCuISO8RhDiFUILxxbFbzg+wrhz2L8QEs+CfSMImACB9+SuMy2qhBsJ4ouNkVxKSqdbCyeeHOhD1xaSICqikvRRzL7eRyGqsCtpV1gWuYywiDCOXT6GjZUN/Zr1Y5T3KLq7dKeGVTFXpM9IMfoY9n4Pp/4yKsV5D4Jh78uw1iosJSOLH7ad5stNEVxKyqC7lxOz7w6kiySISkn6KKqxTFMmf8f+TWh4KBtjNpJlyqKtU1v+r8v/MdRjKI62jsVrWGuI2wN7FsChXyH9GtT3hP6vQPu7pVJcFZaSkcX3204xZ1Mkl5Iy6OHtxGcDWtLZs4GlQxMlUJRE0QWYpJSKppL0UYiCnbxykrDwMP6I/IP4tHga2Dbg7tZ3M8p7FC3rtyx+w8nxcPBnI0FcOAw17KDtKKPvoXkPWa21CkvJyGLBViNBxCdn0NO7IU8O9KGThySIqqAoiSK4zKMohFIqBBgO1AXmaa3/tGxEFVd+5UMb2Dbg0YBHCQsP43D8YWqoGvRp1odRXqPo6dYTG6tiPgIyZUPkeiM5HF9urNbqEgjDPwT/O8G2XgmvSFQ0oXtjzXWmm9azJah5ff6OiOdycga9fBry5AAfgiRBVCmFjnoCUEr1BHy01t9cX+vJQWsdVaQTKPU1cAdwQWvtd9P2IcDHgDUwV2v9ThHaqg98oLV+sKD9qnNntv+3BXcntarfihDvEIa1GEYD2xL8Z75y6vpqrT/AtRiwqw/tJhgzppv6FX68qJRC98by4m8HSc3MzrG9VdM6/He0Px2b17dQZKI0FKsz+/qBrwFBQCvgG8AG+B7oUcRzzwc+Bb67qU1rYDYwCIgBdiqllmAkjbdvOf6Bm9aaevn6caIYfhnxC60blKCqbWaaUeNh7wKI3Ghs8+oHg9+E1sOhRq3SCVRUWO+uPJYrSQAkpWVJkqjCivLoaTQQCOwB0FrHKaXqFPUEWutNSimPWzZ3BsK11pEASqlFwCit9dsYdx85XB+m+w6wQmu9p6jnrk4SMxJZGb2ywH2KnSTOHjCSw4GfIS0B6rlD3xeMjmlH9+K1KSqV1IxsFmyL5uzVtDzfj0tILeeIRHkqSqLIuD5MVgMopUqhRiWuwJmbXsdgdJrn59/AQKCeUspba/3FrTsopaYAUwDc3avHDy+TNrH97HZCw0NZe3ot6dnppdd46hU4uNiYFHfuAFjXNGo8BN4Dnn3AShYQrg7SMrP5fts/8yBq1bAiPcuUaz8Xx9uoYS4qnaIkip+VUl8Cjkqph4EHgK/KNqyctNazgFmF7DMHmANGH0V5xGUpZ66dITQilCURSziXfI46NesQ4h1CiHcIE5dNLH7DJhNEbzbuHo4uNZbQaOIPQ98D/3FQWzooq4tb12Lq7uXE55M6EHslNVcfhZ2NNTOCZcJkVVaUUqgfKKUGAdcw+ile1VqvLuF5Y4Gb14h2u76tRKryMuMpmSmsil5FaHgoey7swUpZ0c2lG890fIZ+7v2oZW30D9xW+dAbrsbAvoVGgkg4BbXqGZ3SgfeAS/syuiJREaVnGct9z14fzvlrxmJ9syYGmmdSd/Iw9rsx6snF0Y4Zwa0ICXS1XNCizBVp1FOJT2L0UfxxY9STUqoGcAIYgJEgdgJ3a60Pl8b5qsqoJ5M2sfv8bkLDQ1l9ajWpWal41PVglPcoRrQYQRP7JsVvPCvDGM66dwFErANtMsqIdph8vYyoPEqoTjKyTPy8y0gQZ6+m0cmjPtMHtpTVXKuZYo96KoUTLwT6Ag2VUjHAa1rreUqpJ4BVGCOdvi6NJFFV7ijikuIIiwhjSfgSYpJisLexZ5jnMEK8QwhoFFCy/7gXjhpzHg4sgpR4qOMCvZ6B9v+CBp6ldxGiUsjMNrF4dwyfrgsnNiGVDu6OUlFO5FIudxTlrTLeUaRmpbLm1BrCIsLYcXYHGk2Xpl0Y5T2Kgc0HYlejiL/h57Xi6vGVxlpLiXEQuxusbKDVUOPuwau/FAKqhjKzTfy2J4ZP1oUTcyWVgGaOPD2oJb2lJnW1VuI7CqWUDeAHxN5SQ7vCqGx3FFpr9l/cb9SXjl5JcmYyrg6uPNb+MUZ5jcLFoRhrIt2oM33nN2BVw6gxHbHWeK9RGwj+L7QbD/YNS/VaROWQlW3i972xfLIunNOXU2jnVo83R/nRt1UjSRAiX/kmCqXUF8AnWuvDSql6wFYgG2iglHpWa70wv2MtRWu9FFgaFBT0sKVjKcj55PMsjVxKWHgY0deisathx6DmgwjxDqFjk47Fqy99g5M3tBoGC0KMfgcU+AyGPs+Da0dZb6mayjZpwvbFMmvtSaLjU/B1qcvcyUEMaNNYEoQoVEF3FL201o9e//v9wAmtdYhSqimwAqhwiaIiS89OZ/3p9YRGhLI1bismbaJD4w484PcAgz0GY29TgukpWRlwYoXxeCl8jZEg6roZS2v0eBIGvVF6FyIqlWyT5o8DcXy85iSRl5Jp41yXL+/pyOC2TSRBiCIrKFFk3PT3QcAvAFrrcxX1A1bRHj1prTkcf5jQ8FBWRK3gWsY1mto35SH/hxjlNQr3uiWcGHjukLHe0oGfrndMO0PP6dCwJaz6P+j9nFFn2nuAlA6tZkwmzbKDZ/l47UnCLyTRqkkdvpjUgcFtm2JlVTH//4qKq6BEkaCUugNj+GoP4EEwD22tkGMnK8qjp1vrS9eyrsUA9wGM8h5Fl6YlqC8N/8yY3vs9nN1ndEy3Hm7Me/DqD6f+zllXWupMVysmk2bFoXN8vPYEJ84n4dPYgdl3d2ConyQIUXwFJYpHMGZDNwWe0lqfu759ALCsrAOrbPKsL92oHa90faVk9aXBmDEdtcFYqfXoUshOhyZ+MORdY8a0/U2T6WL35EwKnr2N17F7JFFUYSaT5s8j5/hozUmOnUvEq5E9syYGMtzfGWtJEKKEqtTw2JsePT188uTJcjnn8cvHCQ0PNdeXbmTXiBFeIxjlNapk9aUBrkTDvh+Nr6tnwNYR2t1lzHlwDpCOaYHWmtVHzvPRmpMcOXuNFg3tmTbAhxEBLpIgxG3Lb3hsvolCKVXY2krTSim2UlfW8yiupF1hedRyQsNDc9SXDvEOoZtLt+LXlwajxvTRpbDve4jaBChjKe/ASdBqONjYltp1iMpLa826Yxf4aM1JDsZepblTbZ4c4MPIABdqWMuCjaJ4ijOP4lHgEPAzEIdRArXayjJlmetLb4jZUHr1pcGoMR27x1hO40aNacfm0O8lCJgIjs0Kb0NUC1prNpy4yEerT7A/5irNGtjx/p3tGB3oKglClJmCEoUzMA4YD2QBPwGLtdYJ5RBXucqvfKiTrRPzgucRGh7K0oilpVtfGiDpojFiae/3cPHoPzWmAycZNaZlKW9xndaaTScvMXP1CfadScDV0Y53x/ozpoMbNpIgRBkrailUN2AC8DTwvNZ6QVkHVhzF7aMorHxoDVWD3m69GeU9il5uvYpfXxogOwvCVxvJ4cRKMGWBWyej38FvjNSYFjlorfk7PJ6Za06w+9QVXOrZ8kR/H+7s6EbNGpIgROkqSSnUDsBEjLkUK4DdpR9e6SiL4bEzgmYwvMVwnOwKWKa7KC6eMPod9i+CpPNg3wi6PgbtJ0HjEpQnFVXWlohLfLT6JDuiL+Ncz5Y3Q/y4K8iNWjVkbS5RvgpawuM/wHDgKLAIeFFrnVVegVUUk30nF//gtGtw+Hfj7iFmByhraBlsPFryGQzWJbgzEVXW9kjjDmJb5GWa1K3FGyN9Gd+pGbY2kiCEZRR0R/EyEAUEXP/67/UZ2QrQWut2ZR9eJaQ1nNpiJIcjoZCZAg1bwaA3jcX46pSghoSo0nZFX2bmmhP8HR5Pozq1eG1EWyZ2dpcEISyuoEQhxQlux9VY2P+jMSnuShTUrGNMhgu8B9yCZM6DyNfuU1f4aM0JNp+8REOHmrw8vA2TujaXBCEqjHwThdb6VF7blVJWGH0Web5fGRWrfChAVvr1KnHf56wS1/cFo0pczRIs9CeqlNC9sbnKh3o0tGfm6hNsPHERJ/ua/N+w1kzq2pzaNcu8npgQt6WgCXd1gamAK7AEWA08ATwD7NdajyqvIIuq3GZmnz1gJIeDPxtrL9V1g/YTof3d0KCEs7FFlRO6N5YXfztIama2eZuVApOG+rVtmNLbi8ndmmNfSxKEsKzizMwOA65g1KEYADTG6J94Umu9r+xCLbkymZmdcvn6YnwL4NwBsK4Jre8wOqZb9JUqcSJfPd5ZR2xCaq7tdW1rsOXFAThIghAVRHGGx7bQWvtfP3gucBZw11qnlVGMlpNX+dCoTRCzG5z9jbuHY8sgO8NYY2nYB+A3Fmo3sFjIovKIyyNJACSmZUmSEJVCQZ/SzBt/0VpnK6ViqmSSgH/Kh95YdXX/T7D0SaOPIeUS2NWHoAeuL8Yng71E0Rw/l8jHa0+Q35RWF8cKuVq/ELkUlCgClFLXrv9dAXbXX98YHluCdbMrmBtLcS+621ih9eoZQIFHz+uL8Q2FGrUsHKSoLMIvJPLRmpMsO3gW+5o1CG7bmI0nL5GWaTLvY2djzYzgVhaMUoiiK2jUU/V66O7Z26gMF7vbGLk0+kuo52rpqEQlEnExiVlrT7Jkfxx2NtY83teLh3q2oL59zTxHPYUEyudLVA7ygPSGqE1wOQp6z4BdX8PlCEkUokiiLiXzydqThO6LpVYNax7p7cWU3i1oYF/TvE9IoKskBlFpValEUeya2VGbjD6Ku769Xj60t5QPFYU6HZ/CrHUn+X1vLDbWiod6tWBK7xY0dJDHlKJqqVIV7m647eGx+Y16it0DPZ8q7fBEJXfmcgqfrgtn8Z4YalgpJnVtziN9WtC4jhSVEpVbsVePrRbySgY37iyEuC7mSgqz14fzy64YrKwU93RtzuN9vWhcVxKEqNokUQhRiLNXU5m9Ppyfdp5Bobi7izuP9/WmaT1JEKJ6kEQhRD7OX0vjs/XhLNxxBo3mrqBmTO3nLfMfRLUjiUKIW1y4lsbnGyP4YftpTCbNuCA3pvbzxq1+bUuHJoRFSKIQ4rqLiel8uTGCBdtOkWXSjO3gyr/7+9CsgSQIUb1JohDVXnxSOnM2RfLt1mgyskyMDnTj3/298Wgoy8QLAZIoRDV2JTmDOZsj+XZLNGmZ2Yxq78q/+3vTopGDpUMTokKRRCGqnYSUDOZujuKbv6NIycxmRDsXpg3wwbuxJAgh8lLhE4VSqg3wJNAQWKu1/tzCIYlK6mpqJvP+iuKbv6JITM9ieDtnnhzgQ8smdSwdmhAVWpkmCqXU18AdwAWttd9N24cAHwPWwFyt9Tv5taG1Pgo8er0E63eAJApxW66lZfLNX9HM/SuSxLQshvo15cmBPrRuWnUWQBaiLJX1HcV84FOMH/AAKKWsgdnAICAG2KmUWoKRNN6+5fgHtNYXlFIjgceABWUcr6hCktKzmP93FF9tjuJqaiaD2zbhyYE++LrUs3RoQlQqZZootNablFIet2zuDIRrrSMBlFKLgFFa67cx7j7yamcJsEQptQz4Ma99lFJTgCkA7u7upXMBolJKTs/i263RzNkUSUJKJgPbNOapgS3xc5UEIURxWKKPwhU4c9PrGKBLfjsrpfoCY4BawPL89tNazwHmgLEoYCnEKSqZlIwsvtt6ijmbIrmcnEG/Vo14amBLApo5Wjo0ISq1Ct+ZrbXeAGwoyr7FXmZcVGqpGdl8v+0UX26K4FJSBr1bNuKpgT50cK9v6dCEqBIskShigWY3vXa7vq3EtNZLgaVBQUEPl0Z7omJLy8zmh+2n+XxDBJeS0unp3ZDpg3zo2LyBpUMTokqxRKLYCfgopTwxEsQE4O7SaFjuKKqHtMxsFu04zWcbIriQmE63Fk589q8OdPaUBCFEWSjr4bELgb5AQ6VUDPCa1nqeUuoJYBXGSKevtdaHS+N8ckdR9dxca9q5ni3dvRvy18lLnLuWRmfPBnw8IZBuXk6WDlOIKk0q3IkKK3RvLC/+dpDUzOwc2z2davPWaH+6ezmhlLJQdEJUPflVuLOyRDBlRSk1Qik15+rVq5YORZSC91Ydy5UkADKyTfTwbihJQohyUqUShdZ6qdZ6Sr16Ml6+MsvKNvHzrjPEJaTl+X5+24UQZaPCD4+9HdKZXbllZZsI2xfHrHUnORWfgo21IjM796NRqTAnRPmSOwphcdkmze97Yxg0cxPP/LIf+5o1+GpyEO+NbYedjXWOfe1srJkR3MpCkQpRPVWpOwpRuWSbNH8ciGPW2pNEXEymddM6fDGpI8G+Tcz9D0op86gnF0c7ZgS3IiTQ1cKRC1G9SKIQ5c5k0iw7eJaP154k/EISrZrU4fN/dSDYtylWVjk7qEMCXSUxCGFhVSpRSB9FxWYyaVYePsfHa05y/HwiPo0d+PTuQIb5OedKEEKIikPmUYgyZzJp/jxyjo/WnOTYuUS8Gtnz5MCWDPd3xloShBAVRn7zKKrUHYWoWLTWrD5yno/WnOTI2Wu0aGjPR+PbMyLARRKEEJWIJApR6rTWrDt2gY/WnORg7FU8nGrz4V0BjAxwoYZ1lRpoJ0S1UKUShfRRWJbWmg3HL/LRmhPsj7mKe4PavH9nO0YHukqCEKISkz4KUWJaazadvMTM1SfYdyYBt/p2TOvvw+gOrthIghCi0pA+ClHqtNb8FW4kiD2nE3B1tOPtMf6M7eBGzRqSIISoKiRRiNumtWZrRDwz15xgZ/QVnOvZ8v9G+zGuYzNJEEJUQZIoxG3ZFhnPzNUn2B51mSZ1a/HmKF/u6tSMWjWsCz9YCFEpValEIZ3ZZWdH1GVmrj7B1sh4Gtepxesj2jKhszu2NpIghKjqpDNbFGj3qcvMXH2Sv8Iv0dChFo/39eLuLpIghKiKpDNb3JY9p68wc/UJNp+8REOHmrw8vA3/6tIcu5qSIISobiRRiBz2n0lg5poTbDh+kQb2NXlxaGvu6dac2jXloyJEdSX/+wUAB2OuMnPNCdYdu0D92jY8P6Q1k7s1x76WfESEqO7kp0A1dyj2Kh+tOcmao+epZ2fDjOBW3NvdAwdJEEKI6+SnQTV19Ow1PlpzglWHz1PXtgbPDGrJfT08qGNrY+nQhBAVjCSKaub4uUQ+WnOCFYfOUce2Bk8N9OGBnp7UlQQhhMhHlUoUMo8ifyfPJ/LR2pMsP3gWh5o1mDbAhwd7elLPThKEEKJgMo+iCgrdG2uuM92oTi3cHO3YG5NAbRtr7u/hyUO9PHGsXdPSYQohKhiZR1FNhO6N5cXfDpKamQ3AhcR0LiSmM7BNY96/M4D69pIghBC3R1Zwq2LeXnHUnCRudvRsoiQJIUSxyB1FFXEqPplP1oVz/lp6nu/HJaSWc0RCiKpCEkUld+ZyCp+sO8mve2KpYaWwr2VNcnruOwoXRzsLRCeEqAokUVRSZy6nMHt9OIt3x2BlpZjcrTmP9fFiS0R8jj4KADsba2YEt7JgtEKIykwSRSUTm5DK7PXh/LLrDArFpK7NeayvF03q2gIQEugKYB715OJox4zgVubtQghxuyRRVBJnrxoJ4qedRoKY0Mmdx/t54Vwv9yOlkEBXSQxCiFJTKRKFUsoe2Ai8rrX+w9LxlKdzV9P4bEM4i3acQaO5K6gZU/t5S5+DEKLclGmiUEp9DdwBXNBa+920fQjwMWANzNVav1NIU88DP5dZoBXQ+WtpfL4hgh93nMZk0owLasbUfl641a9t6dCEENVMWd9RzAc+Bb67sUEpZQ3MBgYBMcBOpdQSjKTx9i3HPwAEAEcA2zKOtUK4kHg9QWw/TZZJM66jG1P7edOsgSQIIYRllGmi0FpvUkp53LK5MxCutY4EUEotAkZprd/GuPvIQSnVF7AH2gKpSqnlWmtTHvtNAaYAuLu7l+JVlI+Liel8uTGCBdtOkWXSjAl05d/9fXB3kgQhhLAsS/RRuAJnbnodA3TJb2et9UsASqn7gEt5JYnr+80B5oCx1lNpBVvWLiWlM2dTJN9tjSYjy8ToQDf+3d8bj4b2lg5NCCGAStKZDaC1nl/YPpVp9djLyRl8uSmC77acIj0rm5D2rvx7gA+ekiCEEBWMJRJFLNDsptdu17eVmNZ6KbA0KCjo4dJoryxcSc5gzuZIvt0STWpmNiMDXJg2wAevRg6WDk0IIfJkiUSxE/BRSnliJIgJwN0WiKNcJaRkMHdzFN/8HUVKZjZ3tHPhyQHeeDeuY+nQhBCiQGU9PHYh0BdoqJSKAV7TWs9TSj0BrMIY6fS11vpwKZ2vwj16upqSyby/Ivnm72gS07MY3s6ZJwf40LKJJAghROUghYvKyNXUTL7+K4qv/44iMS2LoX5NeXKgD62b1rVoXEIIkZ9qUbioItxRJKZl8s3f0czdHMm1tCyCfZvw5ICWtHWRBCGEqJyqVKKwZGd2UnoW8/+O4qvNUVxNzWRQ2yY8NdAHX5d65R2KEEKUqiqVKCwhKT2Lb7dE89XmSBJSMhnQujFPDWyJv5skCCFE1VClEkV5PnpKTs/iu62nmLMpgispmfRr1YinBrYkoJljmZ9bCCHKU5VKFOXx6CklI4vvt53iy42RxCdn0KdlI54a6EOge/2yOqUQQlhUlUoUZSk1I5sftp/ii40RXErKoJdPQ54a2JKOzSVBCCGqtiqVKMri0VNaZjY/bD/NFxsjuJiYTg9vJ74Y2JIgjwaldg4hhKjIZB7FdaF7Y3OUD31qoA/J6Vl8tiGCC4npdGvhxPRBLensKQlCCFE1VYt5FMUVujeWF387SGpmNmDUpZ6x+AAAnT0bMGtiIF1bOFkyRCGEsBhJFMD7q46bk8TNGjrU5KcpXVFKWSAqIYSoGKwsHUBpUkqNUErNuXr16m0dF5eQmuf2+KQMSRJCiGqvSiUKrfVSrfWUevVub7Kbi6PdbW0XQojqpEoliuKaEdwKOxvrHNvsbKyZEdzKQhEJIUTFIX0UQEigK0COUU8zgluZtwshRHUmieK6kEBXSQxCCJEHefQkhBCiQFUqURR31JMQQoj8ValEUdxRT0IIIfJXpRKFEEKI0ieJQgghRIGq5KKASqmLwCkLh1EPsERnSVmct6RtFvf42zmuqPsWtl9h7zcELhUxporOEp/Rivj5LG4blvh8FrZPST+fzbXWjXJt1VrLVxl8AXOqynlL2mZxj7+d44q6b2H7FeH9XZb4dy2LL0t8Rivi57O4bVji81nYPmX1+ZRHT2VnaRU6b0nbLO7xt3NcUfctbD9L/btZgiWutSJ+PovbhiU+n7d73lJRJR89CVFWlFK7dB7r9QtREZTV51PuKIS4PXMsHYAQBSiTz6fcUQghhCiQ3FEIIYQokCQKIYQQBZJEIYQQokCSKIQoJUqpFkqpeUqpxZaORQgApZS9UupbpdRXSql/FbcdSRRCAEqpr5VSF5RSh27ZPkQpdVwpFa6UeqGgNrTWkVrrB8s2UlHd3eZndQywWGv9MDCyuOeURCGEYT4w5OYNSilrYDYwFGgLTFRKtVVK+Sul/rjlq3H5hyyqqfkU8bMKuAFnru+WXdwTSoU7IQCt9SallMctmzsD4VrrSACl1CJglNb6beCOcg5RCOD2PqtADEay2EcJbgzkjkKI/Lnyz29jYPyny7derlLKSSn1BRColHqxrIMT4ib5fVZ/A8YqpT6nBEt/yB2FEKVEax0PPGrpOIS4QWudDNxf0nbkjkKI/MUCzW567XZ9mxAVTZl+ViVRCJG/nYCPUspTKVUTmAAssXBMQuSlTD+rkiiEAJRSC4GtQCulVIxS6kGtdRbwBLAKOAr8rLU+bMk4hbDEZ1UWBRRCCFEguaMQQghRIEkUQgghCiSJQgghRIEkUQghhCiQJAohhBAFkkQhhBCiQJIohCgjSimPW5eCFqIykkQhhBCiQJIohChb1terix1WSv2plLJTSnkppVYqpXYrpTYrpVpbOkghCiKJQoiy5QPM1lr7AgnAWGAO8G+tdUfgWeAzy4UnROFkmXEhylaU1nrf9b/vBjyA7sAvSqkb+9Qq/7CEKDpJFEKUrfSb/p4NNAEStNbtLROOELdPHj0JUb6uAVFKqXEAyhBg4ZiEKJAkCiHK37+AB5VS+4HDGLWNhaiwZJlxIYQQBZI7CiGEEAWSRCGEEKJAkiiEEEIUSBKFEEKIAkmiEEIIUSBJFEIIIQokiUIIIUSBJFEIIYQo0P8HfRP7vnSwIJ8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "he_converge()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c96b4e5-c1c0-42c8-ae86-b5697b50ae3e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python(Proteus Dev)",
   "language": "python",
   "name": "python-proteus-dev"
  },
  "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.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}