laplace-fvm-inclined-poiseuille-fracture.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-04-20T04:03:24.748382Z",
     "iopub.status.busy": "2023-04-20T04:03:24.747467Z",
     "iopub.status.idle": "2023-04-20T04:03:25.158517Z",
     "shell.execute_reply": "2023-04-20T04:03:25.157738Z",
     "shell.execute_reply.started": "2023-04-20T04:03:07.085301Z"
    },
    "papermill": {
     "duration": 0.43088,
     "end_time": "2023-04-20T04:03:25.158703",
     "exception": false,
     "start_time": "2023-04-20T04:03:24.727823",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "import scipy\n",
    "from scipy.sparse import csc_matrix, linalg\n",
    "\n",
    "from math import factorial, floor\n",
    "import numpy as np\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "font = {'size'   : 14}\n",
    "plt.rc('font', **font)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-04-20T04:03:25.192011Z",
     "iopub.status.busy": "2023-04-20T04:03:25.190923Z",
     "iopub.status.idle": "2023-04-20T04:03:25.194217Z",
     "shell.execute_reply": "2023-04-20T04:03:25.193670Z",
     "shell.execute_reply.started": "2023-04-20T04:03:07.108130Z"
    },
    "papermill": {
     "duration": 0.023589,
     "end_time": "2023-04-20T04:03:25.194350",
     "exception": false,
     "start_time": "2023-04-20T04:03:25.170761",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def Vfunction(i, N):\n",
    "    sum1 = 0.\n",
    "    kmin = int(floor((i + 1.) / 2.))\n",
    "    kmax = min(i, N)\n",
    "    for k in range(kmin, kmax+1):\n",
    "        sum1 += (1.*(k**N) * factorial(2 * k) / (factorial(N - k) * factorial(k) * factorial(k - 1) * factorial(i - k) * factorial(2 * k - i)))\n",
    "    return ((-1.)**(N + i)) * sum1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-04-20T04:03:25.225272Z",
     "iopub.status.busy": "2023-04-20T04:03:25.224174Z",
     "iopub.status.idle": "2023-04-20T04:03:25.227084Z",
     "shell.execute_reply": "2023-04-20T04:03:25.227533Z",
     "shell.execute_reply.started": "2023-04-20T04:03:07.127163Z"
    },
    "papermill": {
     "duration": 0.02118,
     "end_time": "2023-04-20T04:03:25.227712",
     "exception": false,
     "start_time": "2023-04-20T04:03:25.206532",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def solve(matrixR,vectorF):\n",
    "    return linalg.spsolve(csc_matrix(matrixR, dtype=float), csc_matrix(vectorF, dtype=float)) # only for nCells <= 100 x 100\n",
    "\n",
    "    # Iterative solvers###################################\n",
    "    # Minimal residual iterative method (krylov) # can reach to nCells = 400 x 400 or more\n",
    "    \n",
    "    '''A = csc_matrix(matrixR, dtype=float)\n",
    "    A = A + A.T\n",
    "    vectorU, exitCode = scipy.sparse.linalg.minres(A, vectorF,tol=1e-100, maxiter=20000, show=False)\n",
    "    print(np.allclose(A.dot(vectorF), vectorU))\n",
    "    return vectorU'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-04-20T04:03:25.255345Z",
     "iopub.status.busy": "2023-04-20T04:03:25.254229Z",
     "iopub.status.idle": "2023-04-20T04:03:25.296773Z",
     "shell.execute_reply": "2023-04-20T04:03:25.297262Z",
     "shell.execute_reply.started": "2023-04-20T04:03:07.152662Z"
    },
    "papermill": {
     "duration": 0.058224,
     "end_time": "2023-04-20T04:03:25.297433",
     "exception": false,
     "start_time": "2023-04-20T04:03:25.239209",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def rigidity(nx, ny, s, x,y, cxx,cyy,cf, phi):\n",
    "    nCell = nx*ny\n",
    "    M = np.zeros((nCell,nCell))\n",
    "    \n",
    "    for i in range(1,nx-1):\n",
    "        for j in range(1,ny-1):\n",
    "            iCell = i*ny + j\n",
    "            iW = (i-1)*ny + j\n",
    "            iE = (i+1)*ny + j\n",
    "            iS = iCell - 1\n",
    "            iN = iCell + 1\n",
    "\n",
    "            l1 = (y[i][j+1]-y[i][j-1])/2.\n",
    "            l2 = (x[i+1][j]-x[i-1][j])/2.\n",
    "\n",
    "            c1xi  = (cxx[i][j] + cxx[i-1][j]) /2. /l2 /(x[i][j] - x[i-1][j])\n",
    "            c2eta =-(cyy[i][j+1] + cyy[i][j]) /2. /l1 /(y[i][j+1] - y[i][j])\n",
    "            c3xi  =-(cxx[i+1][j] + cxx[i][j]) /2. /l2 /(x[i+1][j] - x[i][j])\n",
    "            c4eta = (cyy[i][j] + cyy[i][j-1]) /2. /l1 /(y[i][j] - y[i][j-1])\n",
    "\n",
    "            M[iCell][iCell] = s + c1xi - c2eta - c3xi + c4eta\n",
    "\n",
    "            M[iCell][iW] = -c1xi         \n",
    "            M[iCell][iN] = c2eta\n",
    "            M[iCell][iE] = c3xi\n",
    "            M[iCell][iS] = -c4eta\n",
    "\n",
    "    \n",
    "    # Pressure boundary condition on top and bottom sides  \n",
    "    for i in range(0,nx):\n",
    "        for j in range(0,nCell):\n",
    "            M[(i+1)*ny-1][j] = 0. \n",
    "            M[i*ny][j] = 0.\n",
    "\n",
    "        M[(i+1)*ny-1][(i+1)*ny-1] = 1. \n",
    "        M[i*ny][i*ny] = 1. \n",
    "\n",
    "    # Flow rate boundary condition on left and right sides    \n",
    "    for i in range(0,nx):\n",
    "        for j in range(0,nCell):\n",
    "            M[i][j] = 0.\n",
    "            M[nCell-1-i][j] = 0.\n",
    "            \n",
    "        M[i][i] = 1.\n",
    "        M[i][i+ny] = -1.\n",
    "\n",
    "        M[nCell-1-i][nCell-1-i] = 1.\n",
    "        M[nCell-1-i][nCell-1-i-ny] = -1.\n",
    "    return M\n",
    "\n",
    "def addFracture(M, fracture):\n",
    "    # Fracture\n",
    "    for idx in range(1,len(fracture)-1):\n",
    "        # Location of a cell along the fracture\n",
    "        i,j = fracture[idx]\n",
    "        iCell = i*ny + j\n",
    "        \n",
    "        # Initialize the line of the rigidity matrix at iCell\n",
    "        for k in range(0,nCell):\n",
    "            M[iCell][k] = 0. \n",
    "            \n",
    "        iW = (i-1)*ny + j\n",
    "        iE = (i+1)*ny + j\n",
    "        iS = iCell - 1\n",
    "        iN = iCell + 1\n",
    "\n",
    "        l1 = (y[i][j+1]-y[i][j-1])/2.\n",
    "        l2 = (x[i+1][j]-x[i-1][j])/2.\n",
    "        \n",
    "        cfW  = phi*(cxx[i][j] + cxx[i-1][j]) /2. /(x[i][j] - x[i-1][j])\n",
    "\n",
    "        cfE  = -phi*(cxx[i+1][j] + cxx[i][j]) /2. /(x[i+1][j] - x[i][j])\n",
    "\n",
    "        M[iCell][iCell] = s + 2.*cf/l1**2. + cfW - cfE\n",
    "\n",
    "        M[iCell][iN] = -cf/l1**2.\n",
    "        M[iCell][iS] = -cf/l1**2.\n",
    "        M[iCell][iW] = -cfW         \n",
    "        M[iCell][iE] = cfE\n",
    "\n",
    "    # Fracture tip on the bottom side\n",
    "    # The cell corresponding to the fracture tip\n",
    "    i,j = fracture[0]\n",
    "    iCell = i*ny + j\n",
    "    \n",
    "    # Initialize the line of the rigidity matrix at iCell\n",
    "    for k in range(0,nCell):\n",
    "        M[iCell][k] = 0. \n",
    "\n",
    "    iW = (i-1)*ny + j\n",
    "    iE = (i+1)*ny + j\n",
    "    iS = iCell - 1\n",
    "    iN = iCell + 1\n",
    "\n",
    "    l1 = (y[i][j+1]-y[i][j-1])/2.\n",
    "    l2 = (x[i+1][j]-x[i-1][j])/2.\n",
    "\n",
    "    cfW  = phi*(cxx[i][j] + cxx[i-1][j]) /2. /(x[i][j] - x[i-1][j])\n",
    "    cfE  = -phi*(cxx[i+1][j] + cxx[i][j]) /2. /(x[i+1][j] - x[i][j])\n",
    "\n",
    "    M[iCell][iCell] = s + cf/l1**2. + cfW - cfE\n",
    "    M[iCell][iN] = -cf/l1**2.\n",
    "    M[iCell][iW] = -cfW         \n",
    "    M[iCell][iE] = cfE\n",
    "\n",
    "    # Fracture tip on the top side\n",
    "    # The cell corresponding to the fracture tip\n",
    "    i,j = fracture[-1]\n",
    "    iCell = i*ny + j\n",
    "    \n",
    "    # Initialize the line of the rigidity matrix at iCell\n",
    "    for k in range(0,nCell):\n",
    "        M[iCell][k] = 0. \n",
    "\n",
    "    iW = (i-1)*ny + j\n",
    "    iE = (i+1)*ny + j\n",
    "    iS = iCell - 1\n",
    "    iN = iCell + 1\n",
    "\n",
    "    l1 = (y[i][j+1]-y[i][j-1])/2.\n",
    "    l2 = (x[i+1][j]-x[i-1][j])/2.\n",
    "\n",
    "    cfW  = phi*(cxx[i][j] + cxx[i-1][j]) /2. /(x[i][j] - x[i-1][j])\n",
    "    cfE  = -phi*(cxx[i+1][j] + cxx[i][j]) /2. /(x[i+1][j] - x[i][j])\n",
    "\n",
    "    M[iCell][iCell] = s + cf/l1**2. + cfW - cfE\n",
    "    M[iCell][iS] = -cf/l1**2.\n",
    "    M[iCell][iW] = -cfW         \n",
    "    M[iCell][iE] = cfE\n",
    "    \n",
    "    return M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-04-20T04:03:25.326453Z",
     "iopub.status.busy": "2023-04-20T04:03:25.325417Z",
     "iopub.status.idle": "2023-04-20T04:03:25.349725Z",
     "shell.execute_reply": "2023-04-20T04:03:25.350263Z",
     "shell.execute_reply.started": "2023-04-20T04:03:07.196507Z"
    },
    "papermill": {
     "duration": 0.041281,
     "end_time": "2023-04-20T04:03:25.350419",
     "exception": false,
     "start_time": "2023-04-20T04:03:25.309138",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def generate_fracture(x,y,nx,ny, tip1, tip2, withoutExtension = True, verbose=0):\n",
    "    inclined_fracture = []\n",
    "    for i in range(nx-1):\n",
    "        for j in range(ny-1):\n",
    "            dx = x[i+1,j] - x[i,j]\n",
    "            dy = y[i,j+1] - y[i,j]\n",
    "            xCell = x[i,j]\n",
    "            yCell = y[i,j]\n",
    "\n",
    "            cell_in_fracture = line_passes_through_cell(tip1, tip2, xCell, yCell, dx, dy, withoutExtension, verbose)\n",
    "            if cell_in_fracture:\n",
    "                inclined_fracture.append([i,j])  \n",
    "    return inclined_fracture\n",
    "\n",
    "def line_passes_through_cell(point1, point2, xCell, yCell, dx, dy, withoutExtension, verbose):\n",
    "    # Tolerance\n",
    "    tol_x = 0.1*dx\n",
    "    tol_y = 0.1*dy\n",
    "    \n",
    "    # First, check if either point is inside the cell\n",
    "    if (xCell - dx/2 < point1[0] < xCell + dx/2 and\n",
    "        yCell - dy/2 < point1[1] < yCell + dy/2) or \\\n",
    "       (xCell - dx/2 < point2[0] < xCell + dx/2 and\n",
    "        yCell - dy/2 < point2[1] < yCell + dy/2):\n",
    "        if verbose == 1:\n",
    "            print(f'Fracture tip is inside the cell {xCell, yCell}')\n",
    "        return True\n",
    "\n",
    "    # If neither point is inside the cell, check if the line intersects the cell\n",
    "    x1, y1 = point1\n",
    "    x2, y2 = point2\n",
    "\n",
    "    # Check if line is vertical\n",
    "    if x1 == x2:\n",
    "        if x1 >= xCell - dx/2 + tol_x and x1 <= xCell + dx/2 - tol_x:\n",
    "            if withoutExtension: # If extention of the fracture line is not allowed, it is needed to verify if the cell center is between the tips\n",
    "                if min(y1,y2) + tol_y <= yCell <= max(y1,y2) - tol_y:\n",
    "                    if verbose == 1:\n",
    "                        print(f'Vertical fracture passes through the cell {xCell, yCell}')\n",
    "                    return True\n",
    "                else:\n",
    "                    if verbose == 1:\n",
    "                        print(f'Vertical fracture does not pass through the cell {xCell, yCell}, its extension passes through the cell')\n",
    "                    return False\n",
    "            else: # If extention of the fracture line is allowed, it is not needed to verify if the cell center is between the tips\n",
    "                return True\n",
    "        else:\n",
    "            if verbose == 1:\n",
    "                print(f'Vertical fracture does not pass through the cell {xCell, yCell}')\n",
    "            return False\n",
    "    \n",
    "    # Check if line is horizontal\n",
    "    if y1 == y2:\n",
    "        if y1 >= yCell - dy/2 + tol_y and y1 <= yCell + dy/2 - tol_y:\n",
    "            if withoutExtension: # If extention of the fracture line is not allowed, it is needed to verify if the cell center is between the tips\n",
    "                if min(x1,x2) + tol_x <= xCell <= max(x1,x2) - tol_x:\n",
    "                    if verbose == 1:\n",
    "                        print(f'Horizontal fracture passes through the cell {xCell, yCell}')\n",
    "                    return True\n",
    "                else:\n",
    "                    if verbose == 1:\n",
    "                        print(f'Horizontal fracture does not pass through the cell {xCell, yCell}, its extension passes through the cell')\n",
    "                    return False\n",
    "            else: # If extention of the fracture line is allowed, it is not needed to verify if the cell center is between the tips\n",
    "                return True\n",
    "        else:\n",
    "            if verbose == 1:\n",
    "                print(f'Horizontal fracture does not pass through the cell {xCell, yCell}')\n",
    "            return False\n",
    "\n",
    "    # Calculate slope and y-intercept of line\n",
    "    m = (y2 - y1) / (x2 - x1)\n",
    "    b = y1 - m*x1\n",
    "\n",
    "    # Calculate x-coordinates of intersection points with upper and lower cell boundaries\n",
    "    x_lower = (yCell - dy/2 - b) / m\n",
    "    x_upper = (yCell + dy/2 - b) / m\n",
    "\n",
    "    # Check if intersection points bound a point in the cell\n",
    "    if min(x_lower, x_upper) < xCell - dx/2 + tol_x < max(x_lower, x_upper) or \\\n",
    "       min(x_lower, x_upper) < xCell + dx/2 - tol_x < max(x_lower, x_upper):\n",
    "        if min(x1,x2) + tol_x <= xCell <= max(x1,x2) - tol_x:\n",
    "            if withoutExtension: # If extention of the fracture line is not allowed, it is needed to verify if the cell center is between the tips\n",
    "                if verbose == 1:\n",
    "                    print(f'Inclined fracture pass through the cell {xCell, yCell}')\n",
    "                    print(f'x_lower={x_lower}, x_upper={x_upper}')\n",
    "                return True\n",
    "            else: # If extention of the fracture line is allowed, it is not needed to verify if the cell center is between the tips\n",
    "                return True\n",
    "        else:\n",
    "            if verbose == 1:\n",
    "                print(f'Inclined fracture point toward the cell {xCell, yCell} but does not pass through it')\n",
    "                print(f'x_lower={x_lower}, x_upper={x_upper}')\n",
    "            return False\n",
    "    else:\n",
    "        if verbose == 1:\n",
    "            print(f'Inclined fracture does not pass through the cell {xCell, yCell}')\n",
    "            print(f'x_lower={x_lower}, x_upper={x_upper}')\n",
    "        return False\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-04-20T04:03:25.378128Z",
     "iopub.status.busy": "2023-04-20T04:03:25.377061Z",
     "iopub.status.idle": "2023-04-20T04:03:25.383086Z",
     "shell.execute_reply": "2023-04-20T04:03:25.383560Z",
     "shell.execute_reply.started": "2023-04-20T04:03:07.226414Z"
    },
    "papermill": {
     "duration": 0.021773,
     "end_time": "2023-04-20T04:03:25.383745",
     "exception": false,
     "start_time": "2023-04-20T04:03:25.361972",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Geometry\n",
    "Lx = 10.\n",
    "Ly = 10.\n",
    "nx = 101 #issue while defining nx<ny\n",
    "ny = 101\n",
    "x, y = np.mgrid[-Lx/2.:Lx/2.:nx*1j, -Ly/2.:Ly/2.:ny*1j]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "_cell_guid": "b1076dfc-b9ad-4769-8c92-a6c4dae69d19",
    "_uuid": "8f2839f25d086af736a60e9eeb907d3b93b6e0e5",
    "execution": {
     "iopub.execute_input": "2023-04-20T04:03:25.411172Z",
     "iopub.status.busy": "2023-04-20T04:03:25.410131Z",
     "iopub.status.idle": "2023-04-20T04:03:37.634773Z",
     "shell.execute_reply": "2023-04-20T04:03:37.634025Z"
    },
    "papermill": {
     "duration": 12.23973,
     "end_time": "2023-04-20T04:03:37.634926",
     "exception": false,
     "start_time": "2023-04-20T04:03:25.395196",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Apply unit pressure gradient\n",
    "pTop = Ly/2\n",
    "pBot = -Ly/2\n",
    "\n",
    "t = 1e7    # Time, must be large enough to fit with steady-state results\n",
    "\n",
    "muf = 1e-3 # Fluid dynamic viscosity\n",
    "c_fluid = 5e-10 # Fluid compressibility\n",
    "permeability = 1e-18 # Intrinsic permeability\n",
    "kappa = permeability/muf \n",
    "phi = 0.1 # Porosity\n",
    "\n",
    "# Biot's modulus\n",
    "M = 1./(phi*c_fluid) # Simple approximation for pure fluid flow problem without mechanical deformation\n",
    "\n",
    "# Diffusion coefficient\n",
    "c = kappa * M\n",
    "\n",
    "#e = 1e-5\n",
    "#cf = e**3./12./muf/c_fluid  \n",
    "ell = 1.0 # Fracture half length\n",
    "\n",
    "lambda_ratio = 0.5\n",
    "cf = lambda_ratio * 2.0*np.pi * permeability * ell/muf/c_fluid  \n",
    "\n",
    "# Diffusion matrix\n",
    "cxx = np.r_[ c*np.ones( ( int(nx/2.), ny ) ), c*np.ones( ( nx-int(nx/2.), ny ) ) ]\n",
    "cyy = cxx\n",
    "\n",
    "# Initialize the right-side vector\n",
    "nCell = nx * ny\n",
    "f = np.zeros(nCell)\n",
    "\n",
    "# Update the right-side vector with pressure boundaries condition \n",
    "for i in range(0,nx):\n",
    "    f[(i+1)*ny-1] = pTop\n",
    "    f[i*ny] = pBot\n",
    "  \n",
    "f = f.reshape(len(f),1)\n",
    "\n",
    "# Define a fracture\n",
    "tip1 = [0,-1]\n",
    "tip2 = [0,1]\n",
    "fracture = generate_fracture(x,y,nx,ny, tip1, tip2)\n",
    "            \n",
    "# Solving the linear equation system by the Laplace-FVM method\n",
    "N = 2\n",
    "sum_p = np.zeros(nCell)\n",
    "\n",
    "for j in range(1, 2*N+1):\n",
    "    s = j * np.log(2.) / t\n",
    "    M = rigidity(nx, ny, s,x,y,cxx,cyy,cf,phi)\n",
    "    M = addFracture(M, fracture)\n",
    "    \n",
    "    p = solve(M,f/s)\n",
    "    sum_p += Vfunction(j,N) * p\n",
    "\n",
    "p = sum_p * np.log(2.) / t\n",
    "p = p.reshape(nx,ny)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-04-20T04:03:37.665806Z",
     "iopub.status.busy": "2023-04-20T04:03:37.665076Z",
     "iopub.status.idle": "2023-04-20T04:03:38.047190Z",
     "shell.execute_reply": "2023-04-20T04:03:38.047713Z"
    },
    "papermill": {
     "duration": 0.398233,
     "end_time": "2023-04-20T04:03:38.047884",
     "exception": false,
     "start_time": "2023-04-20T04:03:37.649651",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUAAAAEzCAYAAABAJtrFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAew0lEQVR4nO3dXYwlR3UH8P+ZO97Zb4y9RMYGR1GQkAKIEM+DQwxZoyAckgeLEEEUjCJFXoiVD4mvyEAQEYinyMQCjDU8EIkgiwBKkJHtrAD5BZyQdURgIyxIQGBZyLDBxvsxO+OZOXno2zM9ffujuruqq6rr/5NGM3Onb3fNvX3PrVOnuq6oKoiIUrTkuwFERL4wABJRshgAiShZDIBElCwGQCJKFgMgESVr2XcDio6I6HN9N4IoQU8A51T1eV3v9yIRvdTjeD8B/lVVb+lxV6uCCoDPBfCXvhtBlKC/Bn7U536XALytx/0+CJzoczzbggqARBQXQdxBJKi2zwAc890IIjImAK7w3YgBggqARBQX9gCJJuoQgHXfjQgce4AWLSE76YhCwfOxGXuARJQs9gA9OA7gGd+NoFH5eM55nrVjD9CiGbKTzoTpdjQdPp5znmfN2AMkomSxB2iRIJxB55jf1Sgez/puwEDsARJRshgAiShpMQeRoNrepQhCNAWhDPn0xR4gESWLRRCiiToG4LzvRgSOPUCLuBoMhYbnY7PYe4BcEp+IkhVlADQZOPYxuOzjmCZFI9uFpVSOGfN5Nla78hS461coguq9dlkNJqSTwNYx+54YqQQk18esmpQc83k2RrtiT4FjbjsRecYiiEVL4DxA8if2OXk+xN4DdD4GKCLvFREVkY+7PhYRjWuMMUCXMcRp8BaRGwHcDuDbLo9D5MLxI8AzF323ImyuU2DXMcRZABSR5wD4LIA/BfABk/ssCXBoxVWL+rli5rsF5NPxI272++y2m/32drn/XV0FkT4xpCuXKfAagC+o6tccHoOIPBIAVyx3/zLkPIY4Cd4icjuAFwG4zWDbUwBOAcALDffP1ITGcOggsD6gZ5QCEWC5TxTZwgkROVO4ZU1V1/b2ax5DhrAeAEXkxQA+AuBVqrrZtv38n14DgNVlUdOUw1VqEg2m5mYGppqHDtppxsJ+3ey2v56BXqT3MNE5VV2t3me3GDKEix7gbwI4AeCsiOS3zQC8WkTeDuCIqm44OC4Rjax3D7DZaDHERQD8FwBnSrd9GsD3kUX1+og+A5B6z47Ih//rd7d8DNCy/jGkI+tNV9WnATxdvE1ELgL4uaqetX08IvJIYH04ZswYEvMkbiLyLfJLQUZpuqqeNNpQAAQ2D9C6iE+WpG35bkCgRgqAxjGkI74ciWiYiKNIlOsBEo1i6tkIBRa7l8Aq8BBhPZvjc5GmrpS+034OiiBjSv0lQ0RDRF4EiTMFdjQ7fzAfvddUjumjBxZqNhLS+Z8HwK5fgQioKcjCsemT6+MkMOnqpxKQuhzTVmraFgRdBEnbj62tVWBCCoJMgYkoSZGnwGE1nUWQTFjPyvRxjl9/DIBElDSmwEQTdAQA151sxh6gRTYvhQvrPxtXav/70BS26ZwLfUjGd/rOAEhEyWIA9ICpCY2B55kZjgFa0mVB1NBTk7KwHmkzrts81vzAIcrnme+Us8z3JXrsARJRshgALeI8wHphPVN2+Oi9uCyYpIiLIRBRstgDJKKkRRxFwmp6l8UQQhDWo7co4tTEiK2FBVwJrWBCC0J/CRNRyDgGSDRRKwCsfPz2hHEM0KJUq8AhPAuhvIuHkNYWU9chVd8UKsYMgESUtIijSJxL4qfwzhoyH4Wq0Fe9TlU+Btj1KxBhxe4uq8G4CoJhPSL1fLczliDYpxJbPLfGCoKxVoyZAhNRshgALZoBOOq7EY746PaH9ezaZ9Jrsp0phFCkCU1AKW1XU3+JEJFL7AESUbIYAC3ikviLXKcXoT5OY6S3U0hnfRdPGACJKGkcAxwZlyqnMfA8a8ceoEVdLoULcZKqj0fT5buvzyXxQ0lvu5xnPtJR3xcFMAASUdKYAlvS5UORpsTVszDGs1t38o9RYHDV4+q7X9+9MR8i7wHGeS0wEZEFEcduIvIu8h5gWE0XxLEkvs8xD1vPWFjPfD1bhZIuaS0vnzPHFaGJKFnsARJNFOcBmok4ioTV9BSWxA913l4oZ0I5VbWV5vZNa/Pzccpp7BDsARJRsjgG6AE/rSs9PtLRgwAuj3zM2LAHuJ+I3AngDQBejCxM/RuAO1X1bOudl+B/SfxcKE/qkHfXof+Dj0vhmp7X8vDIkHQZMEtrXc9K8L2aiw0OzpNBcaQDFxOhTwK4B8ArAbwG2VP8FRG5ysGxiMgndx+KdBIjxBHrsVtVX1f8XURuA/ALAL8F4P7GO0+xCBLaZW59epS2/4e819Olx9a3CBHa5XJlsV8+5ygFHhRHOhgj0TuGLLQ9NcKxiGhM440BOokjYzT9bgDfAvDICMciorGNUwV2EkecBkARuQvATQBuUtXKJEZETgE4BQDXPx/jXgoXWvnex2VuXR6Doe3rkvp2SXmncLlclRjmHvbvAZ4QkTOF39dUda3yEAZxpC9nAVBEPgrgzQBuVtUf1G03/6fXAGD1JaKu2kNEDvQPgOdUdbV194ZxpC8nAVBE7kbW6JOq+piLYxA5x0vh2jkcAxwjjriYB/gJALcBuBXAUyJyzfxPF1T1QuOdp1QFdpVed3nGLG+rA/8nMUleuqShXVLqvp8y13Y+xpCmuubgXB8URzpwMQ/wDmQVm68C+Enh610OjkVEPuU9wK5f7UaJIy7mAUrvO6f6ucAD2tq1V7ZtcKy2bbZm3d43l7d3WreZtfTmZh2vGjHqbdaJ6eqMmNrawaA40kFMYYKIQsNrgYkoaaFNJ+sgrAC4BGjERRCT9HKMfZukqNvL9TvcnjWf0dsdz/iZQaVgtl29zWxrnuM1pLsmKfbe/gw26jkMY7Rv23xfSsceIBEliwGQaJo2VpawsmHeu0wSA6A9OzPg0pH0Pqq4KR1tvW9Lugo0p6xbDX/bLp0eXVPfsnIqPCuWMEu7Xs63nR1o3c9m1bFqUmqTlHE37UYWBI35TkcH6R/oh84P9SmoAEhEcVFxO/btWpRN31xZwYENromfkg2sYGXkz0FYXzmMQxuXRj1mdBgA7VERbK6Y5RGm27k0NCU00ZSi1imnrvv/Vr2/qtvLt5mky7OGmbnLCynwYpqa31ZOa2fYxkZNjll5TIOHrdyevfvupd3rK4fbd1TDpPodjvVe91LpPjE+E8bYalABkIjioiI9x7CrRm7HF1QA3MESLqH/O25shvQgu9y3aluT3lz5fl32U9W7WiyCbBd+ru45muwHWCyU1G+7p+0lGFcPbqh+PUDArBAXqqACIBHFRSGjDAW5wgBIRL0ppNc4dSiCCoA7WMJmQzoTmikUQUzS2nJKXLxP3WOwP73dLn1fTHfrCiSbpd+rj9XvGrTaIoglqaTQTedb6OJtORF5xxSYaKJ8zD2MDQOgRawCu7lvVRptcplbOeXd2v197775kEX+tzytPFCoseYp6vJCKrxd2A8Wbivet2gxdXVXBa6be0jTEFQAJKL4xNwDDHLlgbxX0fZ97G1yl3Co8bvpNhvzfdd9L/6c94zz73nPpNhDyX8ub7tZ2GZ997YD898PLfyfe23dv5/8vudxbHfb/Oef4Np9vxe3WS/tZ+/73mNR1571iv9hyGNR9djWfe/6fNb9zfa56PLc7iqvAnf9CoWohvNRvL+2ekg/e+ZXfDejlY8nsEulrekduUsVuK76Www0+YvoaVwJALgSTwPYnwLn42hN1eC2y+Saq8Dm1da+FeMhXFebbfgN+e6jJp/TW/aS1YP6T2de0Pl4L5X/7XU825gCE9EgMafAQQVAFkHc7KPPZW/F28o9ws19aeT+IkFT0WCxGLJ3zHIRpMuiCtXHqtvGPO2zNY8vjKte3WAVmIiSpfAzJGQLAyBRjUs4hMMDFglIg/BKEFsUEvSlcCG807kqhhT/t7o5gnl1tGkeYNXzl/9tZZ4Mbleu8LI/Pc5T1ao01OTyuLr9d1Ws8toQQ0GkC6bARJQ0BsCRXcJhHAaXKk/JBg7s9iDHsokVHOClcI3YA7SoSxV4StViWyeQSYre5RK48n6rqsB18wqL2+Q/5aGkepHT/H7574uXxG20DI80p5f9qr+bli+Fm1pFmMthEVHSWASxRCGt7/I+hfRED10SP9dloYSq9QCrFkgo29wtaOTb1hdBck0LHtQVPyo/HziwooOPK1FcYgpMRMliACSipHEM0JIdLO2u/pEyWyeUScpuUgQpXxJXdQyTNLu8HmDRbOGYJpe72Uknx52bF+4QTx8a+UToIJfDIiIaQ7yhm8gxzjdtxzFAi3awFNwS5DE+ua4qxCYrx1T9Xk55yylx1f3Nqrf1l8t10TQ3z9V809Cq00PE+BrJBRUAiSgusU+EjnIMcN3gAnWTbWzzsZCDycX6tudWmrzj2+4VmPwPthcuMHk+fTznIZ3/eRGk61cowmkJ8iqw2RMX0klQ5PoFUfVu2ycI9q0Q1/2tukLcvPx+10rvkCDYt9LLINiOKTARJYlFEIuyS+FYBBnD0M8VXtxf/alUvn/eGzNpQ3l+YF+hLUIwlSIIAyARJS3mIggDIFENH2sQxib2K0GCanm2HuD4hYsU2T5pu809ND+279VTQl6dKASxp8DOpsGIyB0i8kMRuSwij4rIq1wdi4j82cas85eJMWKIkwAoIm8CcDeAjwB4BYBvAHhQRK63sX8umEBjsL0a9BTlE6G7frVxHUNyrlLgdwD4B1X91Pz3vxCRWwD8GYA76+7UZTWYtu2mUmWLje10yHd6xSDYzOEYYK8Y0pX1HqCIHABwA4DTpT+dBvBK28cjIr9sp8BjxhAXofsEgBmAJ0u3Pwngd8obi8gpAKcA4Nj1z+Gg80hsr4GX97ht99hinmKRggFFkBMicqbw+5qqruV/Q4cYMoTLKrCWfpeK2zD/p9cA4JrV6xb+TkRh6xkAz6nqass2RjFkCBdFkHMAtgFcU7r9l7AY0YmIykaLIdZ7gKq6KSKPAngtgM8X/vRaAF9sum+IS+JPtZhishBBrstcvKGXzS1uO80UeCr/l4vlsIbEkK5cpcB3AfiMiHwTwNcBvB3AtQDudXQ8IvLAYRV4lBjipOWq+jkRuRrA+wE8H8BZAK9X1R+5OB6RC1wS34yL3uxYMcRZEURV7wFwT6f7BLgaTGjG/CS0uspu07BAVW+g7gVi8sIZ8xPyqrhaEn8qXF4K1yeGdBXUtcBEFBcuie+Bj5WeTfhZGXj8HoqPAXwfV2T4eD5NhHb+c0l8S2wviW9D1yqw+xfN4v7LQdCkzfkiT03bltPt4gebd0lr6z5YvemFULWfuiDoMiC7ej6HtjmUIBj7ajBBBUAiigsDoEUhfi6wDz4KHVW35Sd23ef6Vh2jrreX/a3f5wrX8V0goQwDIBElKfYiCAMgUY1NHMABLonfiEviWxT6kvi2V1Cp1z7w7qrQASz+n7Pd4kd1SlyU9waKL4rt3dv2f6+63959uhVIbGxbxXYRJObeUh2mwESUJBZBPNjECg5gw3czaOJ4KVw7jgFapBDjCa9+Jsa6YWvFGZPq8WJ6W1WRLW+zte++TeNiVZXePI2sSo/37tdcBba9ykzzfvaOxUvh2nEMkIiSxBTYoqwIEu47ru/PqAXaCjH7B+yrepZ1hZHi/1bfS8yHHfaOk/cG8/vnL4Zi8WBr97aVfdsU1fUOm15cXV54oaRpMfeWpojPBhH1xh4gESUtlN51H0EFwKwIEuYKHBm7betT/CiXH0z20ZS67y1wUCxa7N93Ob2tKoLk2+bPX/GSxroiSJeFE/Zv037a2uiVxNyzGQsnQhNRspgCE03UBg5ghZfCtWIAtCT0KnDORzW4y2V4i2nxgdq/VX06XHne32w3da1vwyGsA9hLfYtDGRul6m9VmttW/fVVDd6wMOwRc4rYhj1AIkqWgkWQ0XGVDr/O4xiO4fyox/SRjq7jMA7xUrgWLIJY06UK7LdabOfYXarAfaq/e9uar/iSObDvtnI1eBnbOI9j8y2zluWr+OQvhqqJ0OXqr8nCqMXeRZ6OjvnB6j4+cyUmTIGJKGkMgJbEUgTJ2VrEoPkY3QsuiwWTxSJI01qBTR+GVL5P3UcYbO8rcNQXP/a26V8E6TsG5Tp1izkwmOJqMESULE6EJqKkxdzTDSoA7mAJmzshXwpn12zJJIWufjya0m+TgonZ5XHd0+Vc1zX/6tIok/00GfLi3N6J94U9FhZBiChZCon6jSLKAHjpwmEcPsr5WeTWxuUVrBzkRy80UmBriwHQip3tJVy6YFYFNt3OpdnyCFXg5R5V4IZ27aXd9YunmiyaunDM+bZdVnExmQe4e3tDL6PvC3B7q/3037jc/6MXtiMODKZUxehxDNWS7wYQEfkSVOjWHcHGejpFkNzyFebL3JdV9UIXiiBNPcKGHma5J2lWtKlX7sU19dyaehUmPashva+tZ6ffc7Ml6wHG+3gFFQCJKDIad6rPAEhEvalK1D3msALgzhJ2LvovbvTWsyiyaVJo7LDvpZnB5wNfUV0MAfZS5rq1V7oWf4akrCYvrp3tDqexy95KxD2h/qTb4x+YeFtORP4pog78DIBEdS4sA0f9fxZ00FQYAK3ZAXBB7OzLy3/W86BGd6vZaFkXbtppqhzn6W1T2t2S4pqk2PvaY5Iitb2Iur7ItgacR8V/70JYL5EFvuOzYthj7Vngzy4RBc93EB4gzgC4DswXIKZUXBbg4GJv1+0xARwc95DRyT4UJFphBcAdYP7hYu1Mt/PlCkf7XXjGGtKPyme3LpU22SYLQI0pdhOTVKnpxXS5dP8uL7y+L9LLLX9/tud+p4IBkIiSpYj6TSCsALgD4OKIx3P53/dZRMRWe7rUC7r0VJcHDnZ36Sl0eVHZWpMitJ5MaO2porD3+HtgdTEEEblKRD4mIo+JyLqIPC4inxSRq20eh4gCstXjy5KhMcf2ajDXArgOwHsAvAzAWwC8GsB9lo9DRCHIxwA9BUAMjDlWk0BVPQvgDYWb/kdE3g3gyyJyXFWfadzBDtoHnWPjKs3uu9/y/Ywuw+t5rDp9XgB9XzSu0sgY0tMxeC6CDI05Y4wBHkf2MuMSzkRTE2YV2DjmOF0QVUSuBPAhAJ9S1cqHSUROicgZETmDX/zMZXOIurnguwGTdiJ/3c+/TtnYqUnMKTLqAYrIhwG8r2Wzm1X14cJ9jgC4H8ATyPLzSqq6BmANAORXVxXnTVo0Alfz+Loa0kcf2r93nR8M7TmMcX/XQTDiKSQAhvQAz6nqat0fXcacItNT/O8B/GPLNj8uNOQogAfmv/6+qk5tZI+Icm5S4FFijlEAVNVzAM6ZbCsixwA8iOwShVtU1fw9dAf95s+54KodLntVQxblCKXHO6RH5HI+WnjjXGFwNBF6rJhj9eU4b8hpZIOQtwI4Mu+WAsDPVbVujU0iipHnidBDY47t/sgNAG6c//y90t9uBvCw5eMRkU/+q8CDYo7teYAPo/Hq/BY7iKPy5vMCQlvH7jMf0AdbLy6fL9Ipp8/+5wE+jAExJ6xrgYkoLv57gIMwABLVuQDgqO9GRIAB0JIpXgpnIrTL5WwcY4wXxRiXucUwJOMTe4BElCwGQA+YmqTnPIBjIx+TS+K344KoFnWpAoeYmvh4NF1+IqHrydFdXzhdLpO0NTety5BMxD2h3iJfEDWsAEhE8Yk48IcVALfR7V2+SSiXdg3l+hlq2r/PuYExF1HG5Dv95BggESWLAZCIksUiiEUKe2lXaJd2Ta1AEgIfg+8R93ZoUVgBkIjiwiow0UT5mHsYo4h7xWEFwG34n98XS9o4lSq3a0PHp8b6iIZYe1EsghBRsiIvgjj9VLjeLht+d7HNxZrvRXXbXOy5zUbN96a/XSh9r7qtyzZdHi+Xj/+Q/6Fqmz6P7dDns882vh//vvIxwK5fgRBV9d2GXXL1quJ1Z3w3Y3wh9MNDaAMQRjoVQhvGdp882vQpbXXk6lXF7/V4zX6m3/FsC+W0J6JYRfyGEVYADKEI0kVYj96iWAo6fQWUSlWKODAYi3wMMPSXMBGFjPMAiSaK8wDbcRqMRTsA1n03IlCc92eHyzUIU8QASETJ4higBxcBHGndiqbEx8cgbABYGfmYMeIYoCU7qJ54XMV0u1AMfaR9rG7j+uzomjr5mCFQftwjTvecYApMRMliALSoy4citQnrPxvX1Of/lUWcgg3mO/hEPgYY5rXAREQjSLmfRERDcSK0RZwHmOGcv3FFnMIFwXcaPkBYAZCI4sIiCNFEcR5gu8iLIGEFwC7zAH0Y69EK7RPtyro+DhH3EKw/FzE/FlU4BkhEyYo8BY5zGozJUt42lvvuyseVCj6OabJAgO1FBEyKY7YfC5P9+Xj8fZzbTbZ6fAUirB7gDsyfXAbBxWOOOQGaQdDdMU2EEgQ5BkhEyeIYoEXbmP76a5zjF6eIezlORT4GGFYAJKK4MAASTdQ6gEO+GxE4jgFaFvG7iZGp/39Tw0sz23EMkIiSpb4b0F+c8wCJiCxwFgAl85CIqIi80dVxiIiAfjHHZQ/wnYh6dICIItM55jgZAxSRVQB/BeAGAE+a33MbwDMumjQAJ+6RCxGXTgPUN+ZYD4AicgzAfQDepqo/FRHbhyCiYPifBzMk5rjoAd4L4CFVfcDBvokoKEHMhO4dc4wCoIh8GMD7Wja7GcALAbwcwKppA0TkFIBT2W8vQHgTr0JrD1FIevcAT4jImcLva6q6lv/iMuYUiWr7JB4ROQHgRMtmPwZwD4C3IlvXJTeb//6Iqt7UfJyXK3C6tT1EZNs1j6pq5yAi8usKfK3H8a5uPN5YMceoB6iq5wCca9tORN4H4O9KN38HwLsAfMnkWETheAbAcd+NCJybMcCxYo7VMUBVfQLAE8Xb5gOSj6vqD9r3sIPwqsBjYsXZr6oXcsrnowm/RZChMYeXwhHRQN6LIL05D4Cq6mAeTKjLdPhol480LZVj8jxr538aTFmXmBNYD7DLJ6OHWp0d0q6+9/WRpqVyzCmeZzYFMQ2mt8ACIBHFJbweYBeBBcAU1sQnmhL2AIkoWewBEk1USMWGULEHaFFIq8FwTh4B7osN8faeMuwBElGy2AP0gKkJjYGXwrVjD9AiRTjzAEOZZ0V+hTIkQy4EFgCJKD5MgS0JqQhCRO2YAhNRshgAiShZrAJb1GUxBCLyjz1AIkoWe4BEE3UewDHfjQgce4AWcTUYCg3Px2bsARJRstgD9ICXwvnl4/FP5ZixYQ/QoiksiZ8KH49/KseMCXuARJSsuHuAoqq+27BLRH4G4EcOdn0CBh+yHJjY2hxbewG2ueiXVfV5Xe8kIg8ha1NX51T1lh73syqoAOiKiJxR1VXf7egitjbH1l6AbSZgyXcDiIh8YQAkomSlEgDXfDegh9jaHFt7AbY5eUmMARIRVUmlB0hEtIABkIiSlWQAlMxDIqIi8kbf7akjIleJyMdE5DERWReRx0XkkyJyte+2FYnIHSLyQxG5LCKPisirfLepjojcKSL/ISLPiMjPROR+EXmp73aZEpH3zs/bj/tuyxQkGQABvBPZ0jOhuxbAdQDeA+BlAN4C4NUA7vPZqCIReROAuwF8BMArAHwDwIMicr3XhtU7CeAeAK8E8BpklzF8RUSu8tkoEyJyI4DbAXzbd1umIrkiiIisAvhnADcAeBLAH6rqF/y2ypyIvB7AlwFcqareP0FKRP4dwLdV9fbCbd8H8AVVvdNfy8yIyFEAvwBwq6re77s9dUTkOQD+E1kA/ACAs6r6535bFb+keoAicgxZ7+ltqvpT3+3p6TiADQCXfDdERA4geyM5XfrTaWQ9rBgcQ/Y6eMp3Q1qsIXtT+ZrvhkxJUgEQwL0AHlLVB3w3pA8RuRLAhwB8SlVDuAL9BIAZsp500ZMArhm/Ob3cDeBbAB7x3ZA6InI7gBcB+BvfbZma6AOgiHx4Pijc9HVSRG4D8HIA746lzaX7HAFwP4AnkI0JhqQ8jiIVtwVHRO4CcBOAP1DVIMeEReTFyMZX/1hVN323Z2qiHwMUkRNoX43ix8gGvt+KbNHB3Gz++yOqepObFi4ybbOqXppvfxTAA8gCy++q6gXHTTQyT4EvAfgjVf184fZPAHipqv62t8a1EJGPAngzgJtV9THf7akjIn8C4NPYX7SbIXuD2QFwRFU3PDRtEqIPgKZE5DoAzy3d/B0A7wDwJVX9wfitajcft3wQWfC7RVWD+pCKeRHkv1T1VOG27wH4YqhFEBG5G1nwO6mq3/XdnibzYY8XlG7+NIDvI+sZ/rem8iJ2IJkFUVX1CWTp4y4RAYDHAw9+p5EVPm4FcGSeCgPAzwNJie4C8BkR+SaArwN4O7LpO/d6bVWNee/0NmSP51Miko9VXgilZ12kqk8DeLp4m4hcRPb8n/XTqulIJgBG6gYAN85//l7pbzcDeHjU1lRQ1c/NJ2a/H8DzAZwF8HpVdbGwrQ13zL9/tXT73wL44LhNId+SSYGJiMqirwITEfXFAEhEyWIAJKJkMQASUbIYAIkoWQyARJQsBkAiShYDIBEliwGQiJL1/0wxJqie32zHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(Lx/2,Ly/2))\n",
    "plt.pcolor(x, y, p, cmap='jet', vmin=-5, vmax=5)\n",
    "#plt.xlim(-1,2)\n",
    "#plt.ylim(-1,1)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-04-20T04:03:38.078029Z",
     "iopub.status.busy": "2023-04-20T04:03:38.077290Z",
     "iopub.status.idle": "2023-04-20T04:03:38.083070Z",
     "shell.execute_reply": "2023-04-20T04:03:38.083559Z"
    },
    "papermill": {
     "duration": 0.022705,
     "end_time": "2023-04-20T04:03:38.083750",
     "exception": false,
     "start_time": "2023-04-20T04:03:38.061045",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# For void crack, lambda_ratio=infty\n",
    "def analytic2D(ell,x,y):\n",
    "    return np.sign(y)*np.sqrt( np.sqrt( (ell**2 + x**2 - y**2)**2 + 4*(x**2)*(y**2) ) - (ell**2 + x**2 - y**2) )  / np.sqrt(2)\n",
    "\n",
    "#plt.pcolor(x, y, analytic2D(ell,x,y), cmap='jet', vmin=-5, vmax=5)\n",
    "#plt.colorbar()\n",
    "#plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-04-20T04:03:38.114538Z",
     "iopub.status.busy": "2023-04-20T04:03:38.113756Z",
     "iopub.status.idle": "2023-04-20T04:03:38.146445Z",
     "shell.execute_reply": "2023-04-20T04:03:38.145866Z"
    },
    "papermill": {
     "duration": 0.049624,
     "end_time": "2023-04-20T04:03:38.146589",
     "exception": false,
     "start_time": "2023-04-20T04:03:38.096965",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          X         Y\n",
      "0  0.000000  0.007126\n",
      "1  0.187500  0.078385\n",
      "2  0.337500  0.135392\n",
      "3  0.487500  0.192399\n",
      "4  0.646875  0.249406\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "# Read data from text file into a pandas DataFrame\n",
    "df = pd.read_csv('/kaggle/input/poiseuillecrackpressure/pouya2010pressureAroundPoiseuilleCrack_fig5_lambda05.txt', delimiter='\\t')\n",
    "# Print the first 5 rows of the DataFrame to check if it was read correctly\n",
    "print(df.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-04-20T04:03:38.182409Z",
     "iopub.status.busy": "2023-04-20T04:03:38.181603Z",
     "iopub.status.idle": "2023-04-20T04:03:38.184469Z",
     "shell.execute_reply": "2023-04-20T04:03:38.185042Z"
    },
    "papermill": {
     "duration": 0.024684,
     "end_time": "2023-04-20T04:03:38.185229",
     "exception": false,
     "start_time": "2023-04-20T04:03:38.160545",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Extract data along fracture line\n",
    "# fracture_extention = generate_fracture(x,y,nx,ny, tip1, tip2, withoutExtension = False, verbose=0):\n",
    "\n",
    "\n",
    "# Extract a long the vertical middle line\n",
    "p_TFVM = []\n",
    "y_TFVM = []\n",
    "\n",
    "for i in range(0, ny):\n",
    "    p_TFVM.append(p[int(nx/2)][i])\n",
    "    y_TFVM.append(y[int(nx/2)][i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-04-20T04:03:38.218014Z",
     "iopub.status.busy": "2023-04-20T04:03:38.217250Z",
     "iopub.status.idle": "2023-04-20T04:03:38.399653Z",
     "shell.execute_reply": "2023-04-20T04:03:38.398966Z"
    },
    "papermill": {
     "duration": 0.200449,
     "end_time": "2023-04-20T04:03:38.399779",
     "exception": false,
     "start_time": "2023-04-20T04:03:38.199330",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAc8AAAFKCAYAAAB7BdRAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd3iUVfbA8e9NDCmEstSQQBIVFFlFkIggKCC4sgIqqCggAoIIFkREVwyKqLEsRcV1V+NPECWwLIrrioAiTZQiAQIqAaQl9CJISyHl/P54ZyYzqROYzKScz/PM48x977xzEkxO7n3ve64REZRSSinlPj9fB6CUUkpVNJo8lVJKqVLS5KmUUkqVkiZPpZRSqpQ0eSqllFKlpMlTKaWUKiVNnkoppVQpeS15GmMeM8ZsMcactj3WGGN6lPCea4wxK40x6caYA8aYF40xxlsxK6WUUoW5xIuftR/4G/AbVtIeBPzXGNNGRLbk72yMqQksAb4HrgeuBD4GzgFTvBSzUkopVYDxZYUhY8wJYJyIfFDIsZHAm0BDEUm3tY0HRgKNRUsjKaWU8hGfXPM0xvgbY+4HQoHVRXRrD6yyJ06bb4BwILpsI1RKKaWK5s1pW4wx1wBrgCDgLNBbRH4uonsY1lSvsyNOx/YUcv7hwHCA6tWrt2nevLknwlZKKVVJbNiw4biI1L/Y83g1eQLbgVZAbeBuYKYxprOI/FJE//xTs6aIdqtRJB6IB4iJiZHExMSLj1gppVSlYYxJ8cR5vJo8ReQ8sNP2MtEYcz3wFDC0kO6HsUaYzhrY/nsEpZRSykd8fZ+nHxBYxLE1wE3GmCCntluBg8DeMo5LKaWUKpI37/N8wxhzkzEm2nb/5utAZyDBdvx1Y8xSp7fMBtKAj40xVxtj+gDPAVN1pa1SSilf8ubIMwyYhXXdcynWvZt/FZFFtuONgMvtnUXkFNZIMxxIBN7Dur9zqhdjVkopVQkkJCQQHR0N0MYT5/PaNU8RGVza47aVuDeXUUhKKaWqgISEBIYPH05aWprHzunt1bblxunTpzl69ChZWVm+DkUpVYYCAgJo0KABNWvW9HUoykdiY2M9mjihiibP06dPc+TIESIiIggODkbL5SpVOYkI6enpHDhwAEATaBWVmprq8XP6erWtTxw9epSIiAhCQkI0cSpViRljCAkJISIigqNHj/o6HOUjkZGR+AHPe/CcVTJ5ZmVlERwc7OswlFJeEhwcrJdoqrC3xoxhuZ8fcR48Z5WctgV0xKlUFaI/71XYF1/Q+6WXIDfXo6etkiNPpZRSlVxaGowYAX36wMmTVpufH8AGT5xek2clsnfvXowxFFfT150+5UViYiLGGPbu3ev2ewYPHkzPnj3LLiilVPm3eTPExMAHTrtdNmkCK1Z47CM0eVYQxphiH4MHD3brPE2aNOHQoUO0atWqbANWSilvE4Fp0+CGGyA5Oa/9nnushHrTTR77KE2eFcShQ4ccjw8//LBA2zvvvOPWefz9/QkLC+OSS6rs5W6lVAVnrxbk5+dHdHQ0CQkJcOwY9OoFTz4JmZlWx5AQ+L//g//8B/70J4/GoMnzIhT6D1hGwsLCHI/atWsXaKtVq5ajb0pKCrfeeishISG0aNGCJUuWOI7ln7bNyspi1KhRhIeHExgYSJMmTXjuueeKjOPjjz8mNDSURYsW0bx5c0JCQrjjjjs4deoUn332Gc2aNaNWrVoMHDiQ9PS8fcwzMzMZPXo0DRs2JCgoiHbt2vHDDz+4nHvx4sU0b96coKAgbrrpJnbs2FHg81evXk2nTp0ctx+MHDmS06dPX9g3VSlV4dirBaWkpCAipKSk8O+hQ0lv1gy+/jqvY6tWsGEDDB0KZbFgTEQq5aNNmzZSlK1btxZ5zF2zZs2SkJAQwdpbVAAJCQmRWbNmXfS5SzJv3jyx/ulc7dmzRwC58sor5X//+5/s2LFDHnzwQalTp46cOXPGpc/69etFRGTy5MnSuHFjWblypaSkpMiPP/4o06dPL/KzZ8yYIZdccol07dpVEhMTZfXq1dKoUSPp1q2b9OzZUzZv3izLli2T2rVry+TJkx3vGzVqlISFhcmCBQtk69atMmzYMKlevbocPHhQRERSU1MlMDBQHn/8cUlOTpa5c+dKRESEALJnzx4REdmyZYtUr15dJk+eLDt27JC1a9dKu3bt5O6773Z8zqBBg6RHjx4X/T1WlY8nfu6V70VFRTl+5waAvGlN1ro+nnpKJCOj0PcDieKBHOPzJFdWj7JOns7/gM6PqKioiz53SUpKnu+//76jbf/+/QLIqlWrXPrYk+cTTzwht9xyi+Tm5rr12TNmzBBAtm3b5mh7+umnxc/PT44dO+Zoc05iZ8+elYCAAJk5c6bjeHZ2tlx22WUSGxsrIiLjxo2TZs2aucTxyiuvuCTPgQMHykMPPeQSz6ZNmwSQI0eOFPhcpZxp8qwcjDECSFOQ9fmTZoMGIosWFft+TyVPnba9QEWVeyqLMlCl1bJlS8fz8PBwgCKrqwwePJikpCSuuOIKHnvsMb7++mtyS7gfKjAwkCuvvNLxumHDhoSFhVGvXj2XNvtn7tq1i6ysLDp06OA47u/vT/v27dm6dSsAycnJtGvXzuV+vPbt27t87oYNG5g1axahoaGOh/2cu3btKjZmpVTlENmkCYOATUCMU/vKoCDYsgW6d/dKHLpq5AJFRkaSkpJSaLuvBQQEOJ7bk1FRCfG6665j7969LF68mGXLljFo0CCuvfZalixZgp9f4X9b5V9sZIxx+Ux7m/0zrT/2Cr9R3d5m71Oc3Nxchg0bxlNPPVXgWERERInvV0pVcKdOsSI8nGinQcp54IWAAFrGx0PDhl4LRUeeFyguLo6QkBCXtpCQEOLiPFkAyjtq1KjBvffey7/+9S++/vprli1bxs6dOz12/qZNm1KtWjWXBUI5OTmsWbOGFi1aANCiRQvWrVvnkkTXrl3rcp7rrruOX3/9laZNmxZ4aLlFpSq5NWugVSuinX4vbAN6N2pEyxkzGDBwoFfD0eR5gQYMGEB8fDxRUVEYY4iKiiI+Pp4BAwb4OrRSmTp1KnPmzCE5OZmdO3cye/ZsatasSePGjT32GdWrV2fkyJE899xzLFy4kOTkZEaOHMmRI0d49NFHARgxYgR79+5l9OjRbN++nc8++4z333/f5Tx/+9vf+OmnnxgxYgSbNm1i586dLFiwgEceecRjsSqlypmcHHj1VeseTeeCKUOH0vzsWb4+eNAnv3d12vYiDBgwoMIly/xq1KjBpEmT+O233zDG0Lp1axYtWlRgVH2x3nzzTQCGDBnCH3/8QevWrVm8eDGNGjUCrOnu+fPnM2bMGD744APatGnDG2+8wQMPPOA4R8uWLfn+++8ZP348nTp1Iicnh8suu4zevXt7NFalVDmxbx8MHAgrV+a11aoF8fHQt6/v4gKMO9eaKqKYmBgpqgRdcnIyV111lZcjUkr5kv7cl08JCQnExsaSmppKZGQkcXFx1qDkiy+sezTtdWkBOnSAhASIirrgzzPGbBCRmJJ7Fk+nbZVSSvlEYQUPnnz4YX7r2rVgQfeXXrJq015E4vQknbZVSinlE7GxsaSlpTletwTmpKfTbNmyvE6RkdZos2NH7wdYDB15KqWU8gnn++KfANYBLZw73HsvJCWVu8QJmjyVUkr5SGRkJPWAr4BpQJCtPc0Yq6D73LkeL+juKZo8lVJK+cRH99/Pz4DzDrybjWHpm2+WXUF3D9Frnkoppbzr/HkYP56ukya5NP9fjRpUnzaNfm7uT+xLmjyVUkp5z2+/Qb9+1nZhdg0awMyZDPNSXVpP0GlbpZRSZU8EPv4YWrd2TZzdu3u1oLunaPJUSilVtk6dgv79YcgQOHfOagsIgKlTrQ2svVjQ3VM0eaoK7bPPPit0t5YLsWLFCowxHD9+/KLO07lzZx5//HGPxOTLz7AbPHgwL7/8slc+q7QWLFhAq1atStxGT3lHQkIC0dHR+Pn5ER0dTUJCgqOgO//+d17HK6+EdevgqaesAggVUMWMuooaPHgwxhjHFmCXXXYZY8eO5Zz9LzlVpKSkJPr160d4eDiBgYFERkZy++2388UXX+gv3mL8/PPPfPnll4wePRqArKws/va3v9GyZUuqV69Oo0aN6N+/f4F9bDMzM3niiSeoV68e1atX54477mD//v0ufeLi4ujQoQPVq1cv8g+g1NRUevXqRfXq1alXrx6jRo3i/PnzjuM9e/bE39/f+iWtfCp/taB9KSnsHDKE3I4dXQu6DxtmTdu2bu2zWD1Bk2cF061bNw4dOsTu3bt59dVX+ec//8nYsWN9HVa5tmDBAm644QZOnTrFjBkz2Lp1K59//jn33HMPcXFxHDx40Nchllvvvvsud999NzVr1gQgLS2NjRs3Ehsby8aNG/nyyy/Zt28f3bt3Jzs72/G+0aNH8/nnnzNnzhxWrVrF6dOn6dmzJzk5OY4+mZmZ9OnTx5GY88vJyaFHjx6cOXOGVatWMWfOHD777DOefvppl35Dhgxh2rRpZfDVq9JwrhbUGFgKTMjKws/+x2nt2vCf/8CHH0L16j6L02NEpFI+2rRpI0XZunVrkcfKs0GDBkmPHj1c2oYNGyZhYWEiIpKRkSFPPvmkNGjQQAIDA+WGG26QVatWOfouX75cADl27Jijbc+ePQLI+vXrJTc3Vy6//HKZNGmSy2fs2LFDANmwYYOIiEyZMkWuueYaCQkJkfDwcBk6dKicPHmy2Ng//fRTiYmJkdDQUKlfv77cc889sn///gKxfffdd9K2bVsJDg6WNm3aOD7TbubMmRIZGSnBwcHSo0cP+cc//iHW/8aFO3v2rNSrV0969+5dZJ/c3Fy3Yzh+/Ljcf//9EhERIUFBQdKiRQuZPn26y/k6deokjzzyiIwaNUpq164ttWvXlrFjx0pOTo6jT2Zmpjz77LMSEREhISEhEhMTI4sXL3Y5B1DgsXz5csfxxx57zNH/xIkT8uCDD0rt2rUlKChIunbtKr/88otLXJ9//rlcffXVUq1aNWncuLG8+uqrjq+9MNnZ2VKrVi358ssvi+wjIvLrr78KIFu2bBERkT/++EMCAgJk1qxZjj6pqalijHH5Gu3mzZtX6L/hwoULxRgjqampjrZPP/1UAgMD5dSpU462lJQUAeS3334rNs6K+nNfURhjBJDeIL9by4PyHh07iqSk+DpEEREBEsUDOUZHnhVccHAwWVlZADz77LPMnTuX6dOns2nTJq655hq6d+/OoUOH3DqXMYahQ4cyffp0l/bp06fTqlUrrrvuOgD8/Px4++23+fXXX5k9ezY//fQTTzzxRLHnPn/+PBMnTmTz5s0sWLCA48eP069fvwL9xo0bxxtvvMHGjRupW7cuAwYMcGyQvW7dOgYPHszw4cNJSkqiV69evPjii8V+7rfffsvx48d59tlni/263Y0hIyOD6667jgULFvDrr7/y5JNP8sgjj7B06VKXcyQkJJCbm8uaNWv44IMPiI+P5+2333YcHzJkCCtXrmT27Nn8/PPPDBo0iF69erF582YA5s+fz6FDhxyPESNG0LBhQ5o3b17o1zB48GDWrVvHl19+yU8//URISAjdu3cnPT0dgA0bNnDvvffSp08ffv75Z9544w1ef/11/vGPfxT5fdmyZQunTp0iJqb4DShOnz4NwJ9slWA2bNhAVlYWf/nLXxx9mjRpwlVXXcXq1auLPZezNWvWcNVVV9GkSRNH22233UZmZiYbnFZrRkZG0rBhQ1Y6b1ulvO6Kxo15H5gP1LG15QBv1aoFy5dbNWorE09k4PL4KPXIM/9fSt58uCn/yHPdunVSt25d6du3r5w9e1YCAgJk5syZjuPZ2dly2WWXSWxsrIiUPPIUETl06JBccsklsmbNGsc5wsPD5d133y0yrkWLFkm1atVcRlYlSU5OFkD27dvnEpvzyOSHH35w6dOvXz/p1q2by3mGDh1a7MjzjTfeEEBOnDjhaNuyZYtUr17d8bCPkNyJoTD33XefDB061PG6U6dO0qxZM5dR3SuvvCIREREiIrJz504xxkhKvr/E77zzThk5cmSB8//73/+WoKAgx7+J/TPsI0/7zMDKlSsdx//44w+pWbOmfPjhhyIi0r9/f+nSpYvLeSdMmOCIqTBffPGFGGOK/XfNzMyUG2+8UXr16uVoS0hIEH9//wKj2i5dusjw4cMLnKOokefDDz9cIObc3Fzx9/eX2bNnu7S3bt1axo8fX2ScIjryLFNJSfJHeLjL77W9IF0DA11mIMoDdORZNS1evJjQ0FCCgoJo3749N998M++++y67du0iKyuLDh06OPr6+/vTvn17tm7d6vb5w8LC6Nmzp2P0uXjxYn7//XeXTb+XLVvGrbfeSuPGjalRowZ9+vTh/PnzHD58uMjzbty4kTvvvJOoqChq1KjhGM3kX2jSsmVLx/Pw8HAAjh49Clj7MbZv396lf/7X7rjyyitJSkoiKSkJEXGM3N2JIScnh7i4OFq2bEndunUJDQ1l/vz5Bb6Odu3auYxo27dvz4EDBzh9+jQbN25ERGjRogWhoaGOx9dff82uXbtczpOYmMhDDz3ERx99RLt27Qr9epKTk/Hz83P5XtSqVYtrrrnG8W+fnJzs8v8GQMeOHR0xFSY9PZ2AgAD8ilgNmZ2dzQMPPMAff/zBjBkzCu3jTERKvTK6qP7524ODgx2jbOVFIjBtGrRtSy2ntQP/AXo1bsyQjz5y+d1RmXgteRpjxhlj1htjThtjjhljvjLGXF3Ce6KNMVLIo2LdTetBN998M0lJSWzfvp2MjAzmz59PgwYNHNOKhf2ysbfZfwna+wIFEgfAsGHDmDt3LmlpaUyfPp0+ffo4puRSUlLo0aMHV111FfPmzWPDhg2OROu8CtLZuXPnuO222wgJCeHTTz9l/fr1LF68uND3BAQEFIjbvhrWOW53XXHFFQBs27bN0VatWjWaNm1K06ZNC/1+FRfD5MmTmTJlCs888wxLly4lKSmJu+66q8ivvTC5ubkYY1i/fr0jiSclJZGcnOwyZX7w4EHuuusuxowZQ//+/Ys8X3HfF3v8xSWuotrr1avH+fPnXbaMssvOzqZfv35s2bKFpUuXUrduXcexsLAwcnJyCtzyc/ToURqW4n6+sLCwAn+QHT9+nJycnALnOXHiBPXr13f73MoDjh2DXr3gySetcnsAISHw0Uf0zc1ly759lTZxgndX23YG/gncCNwCZAPfGWPqFPcmm+5AI6fHsuK7XwBfTtyWQkhICE2bNiUqKsrll3zTpk2pVq0aP/zwg6MtJyeHNWvW0KKFtcmP/ZeL8zXQpKSkAp/RvXt3atasyfvvv89XX33FQw895DiWmJjI+fPneeutt2jfvj1XXHFFiatVt23bxvHjx3nttde4+eabad68uWMkVxotWrRg7dq1Lm35X+f3l7/8hbp16/L666+X+vMK88MPP9CrVy8GDhxIq1atuPzyy9mxY0eBfuvWrXNJamvXriU8PJyaNWvSunVrRITDhw87krj9ERERAVjXVu+66y7atWtX4j2WLVq0cFxftTt9+jQ///yz49++RYsWLv9v2L8W++xBYVq1agVQYOYiKyuL++67jy1btrB8+XLCwsJcjrdp04aAgACWLFniaNu/fz/JycnceOONxX4tztq3b09ycrLLLS5LliwhMDCQNm3aONoyMjLYtWuX45q88oIlS6BlS6vAgV3r1rBxIzz0ULku6O4xnpj7vZAHEIp1PblXMX2isVYZxpT2/FVlta2zJ598Uho1aiRff/21bN26VR5++GGpXr26HDx4UEREzp8/L02aNJHevXvL9u3b5ZtvvpGWLVu6XPO0e/HFF6VatWoSHR3tcu1q8+bNAsjkyZNl9+7dMnv2bGnSpIkAsmfPnkLjOnr0qAQGBsqYMWNk165dsmDBAmnRooXL6lF3rseuWbNGjDHy2muvyY4dOyQ+Pl7q1q1b7DVPEZH//ve/EhAQILfddpssWrRIdu7cKVu2bJEpU6aIv7+/fPLJJ27HMGbMGImIiJBVq1ZJcnKyPProo1KzZk3p1KmT4z2dOnWS0NBQGTVqlGzbtk3mzZsntWrVclnFPGDAAImMjJR58+bJrl27ZP369TJp0iT5/PPPRUTkwQcflGbNmsmuXbvk0KFDjkdmZqbjM5xX2955553SvHlz+f7772XLli3Sq1cvady4saSlpYmIyIYNG8TPz08mTJgg27dvl1mzZkloaKhMmzat2O/dddddJ2+99ZbjdVZWltx5550SHh4uGzZscInN/lkiIiNGjJDw8HBZsmSJbNy4UTp37izXXnutZGdnO/qkpKTIpk2bZNKkSQLIpk2bZNOmTXLmzBkRsa63X3311dKlSxfZuHGjLFmyRMLDw+Xxxx93iXH58uUSGhoq586dK/Zrqag/9+VKZqbI2LEFhwBjxohkZPg6OrfgoWuevkyejWyJsWMxfezJMxU4CvwI3OPO+ati8nS+VaVatWoFblUREfnxxx/l2muvlaCgIGnXrp0sWLCg0OS5d+9eAWTixIkFPuedd96R8PBwCQoKkltuuUXmzp1bbPIUsRa9XHbZZRIYGCjXX3+9LF68uNTJU0Rk+vTp0qRJEwkKCpLu3bvLu+++W2LyFLGSR9++fSUsLEwuueQSqVOnjnTr1k0++eQTx4IYd2I4ceKE9O7d23HLzTPPPCMjR44skDwfeeQReeyxx6RWrVpSu3ZtGTNmjEviOH/+vEyYMEEuvfRSCQgIkIYNG0qvXr0kMTFRRESioqKEMrpVJSAgwK1bVURE3n//fYmJiSnw/SjsMWPGDEe/9PR0efzxx6VOnToSHBwsPXv2dLnlRMT6/7m4r1HESrA9evSQ4OBgqVOnjjz++OOSke+X9PDhw+WRRx4p9usQqbg/974ya9YsiYqKEmOMREVFyZeTJ4u0aeOaNBs0EFm0yNehlkplSJ7/ATYB/sX0qQc8DbQDYoCXbaPVB4roPxxIBBIjIyOL/ObpD1HJ1q5dK/7+/gVWhKqqJSMjQ6Kjo+X777/3dSiFOnr0qNSpU0d2795dYl/9uXffrFmzJCQkxPEHzSCQM/lHm927ixw+7OtQS81TydMnW5IZY6YCHbFGnTlF9ROR48AUp6ZEY0w94FlgViH944F4gJiYmNKvLlFkZmayb98+xo8fT+/evYmsbPdmqVIJDAxk5syZnDhxwtehFGrPnj3885//5NJLL/V1KJWKvVpQTeB9wOWO7GrV4M03YdSoCluX1hO8njyNMW8B9wNdRGT3BZxiHTDEs1Epuzlz5jB06FCuvfbaAsUSVNV08803+zqEIrVt25a2bdv6OoxKJzU1lXbAbMD5z5JtQPO1ayt8XVpP8OqfDcaYd4D+wC0isq2k/kVoBbhXMkeV2uDBg8nJyWHjxo0ulV2UUlVETg6TatZkFa6J80Ogd5MmmjhtvDbyNMa8BwwE7gJOGmPs69vPishZW5/XgbYi0tX2ehCQhXVtNBfoBTwG/M1bcSulVJWxbx8MHMjTp045mk4CDwOLQkKI99AtX5WBN6dtH7X9d2m+9onAS7bnjYDL8x0fD0RhLRTaATwkIgWud5aWSOmrnSilKiZrnYgq1vz51nZhJ086mn4KDOTezExMVBTxcXGVuuhBaXkteYpIiZlKRAbnez0TmOnpWAICAkhPTyckJMTTp1ZKlUP2UoOqEGlp1qbU8fF5bX5+MGECbZ9/npRLfLKutNyrkt+VBg0acODAASIiIggODtYRqFKVlIiQnp7OgQMHSlUasMrYvBn69YPk5Ly2yEiYPRvy1UJWrqpk8rRv7Hvw4MFCa7sqpSqPgIAAGjZs6Pi5V1h3ar77LjzzTF5dWoC+feGDD6yNq1WxqmTyBCuB6g+TUqoqSEhIIDY2ltTUVFpFRPBVvXpEONe1DgmxkumQIVWjLq0HVNnkqZRSVUFCQgLDhw8nLS2NrsCn+/fTyKnYPtddZ03TXnmlz2KsiKpueQillKoCYmNjyUhLYyLwLdYtDQ5PPw2rV2vivACaPJVSqhI7n5LCd8CL5P3CP4y1zyOTJ0NgoK9Cq9B02lYppSqrJUvY7OdHfdtm7gDfAQ8AQVFRPgurMtCRp1JKVTY5OfDii3DbbY7EmQO8ANwGnAkJIS4uzpcRVng68lRKqcrk0CHo3x9WrHA0pdeqxUNBQcw9epTIyEjitFrQRdPkqZRSlcWSJfDAA3D0aF5b164EJyQwp2FD5vguskpHp22VUqqic5qmdSROPz94+WX45hvQ6koepyNPpZSqyAqZpiUszLp3s0sXn4VV2enIUymlKpCEhASio6Px8/PjgYYNSW/e3DVxdusGSUmaOMuYjjyVUqqCsFcLykhL4yVg/NGjeSMgPz946SV4/nnw9/dZjFWFJk+llKogYmNjqZWWxtdAZ6f2o/7+NFiyREebXqTTtkopVUFckZJCEq6JcwnQMidHE6eXafJUSqnyzraadjHQwN6EVfSgO1otyBd02lYppcqzgwet1bQrVzpGO4eA/sAKIESrBfmEjjyVUqq8WrIEWrWClSsdTYeuvppejRuz0hiioqKIj4/XakE+oCNPpZQqb7KzYeJEiIsDEavNtpq20fPPk6iraX1Ok6dSSpUnTtO0Dlr0oNzRaVullCovvv22wDStFj0onzR5KqWUj9irBQUYw7TatZHu3eHYMeugvTbt4sVam7Yc0mlbpZTyAXu1oFppaSwBOp86lXdQp2nLPU2eSinlA7GxsXRIS2MWefduAqwKCuKmpCQdbZZzOm2rlFLelp3NsJSUQosedMnI0MRZAWjyVEopbzp4ELp1Yzy4FD3oCrwKNNZqQRWCJk+llPKWQlbTLgFaASvRakEViSZPpZQqa9nZ8MILkG817ea77+aRyEiOabWgCkcXDCmlVFkqqujBnDlc27kzu30XmboIOvJUSqmyUlzRg86dfRaWuniaPJVSytOys2H8+ALTtLzyihY9qCQ0eSqllAfYqwVFGMO60FDXou5hYbB0qZVQtah7paDXPJVS6iLZqwU5ih5kZuYd7IHUNHMAACAASURBVNYNZs3S0WYloyNPpZS6SC8+/zzj0tIKFD2YXKuWTtNWUl5LnsaYccaY9caY08aYY8aYr4wxV7vxvmuMMSuNMenGmAPGmBeNMcYbMSulVIkOHmR6amqhRQ+ePX1ap2krKW+OPDsD/wRuBG4BsoHvjDF1inqDMaYm1j3ER4DrgVHAM8CYsg5WKaVKZFtN28mpybnoQWRkpG/iUmXOa9c8ReQ259fGmIHAKaAD8FURbxsAhACDRCQd+MUYcxUwxhgzVcR+NV4ppbwoOxteeglee82xKCgHmAC8DuSi1YIqu1KNPI0x9YwxNxhjAj3w2TVsn3+ymD7tgVW2xGn3DRAORHsgBqWUKp2DB6Fr1wKraZc9/zyzoqIQrRZUJbg18jTG1AA+Au4BBGgG7DbGvA8cFpGXLuCz3wGSgDXF9AkD9udrO+J0bE++OIcDw0GnS5RSZeDbb+GBB/Lu3QTHatpbGzZkr440qwx3R55vAhHAdYDzKHAB0Lu0H2qMmQp0BO4WkZwSuuefmjVFtCMi8SISIyIx9evXL21YSilVuNxcmDhRix4oB3eved4B9BaRJGOMc9JKBi4rzQcaY94C7ge6iEhJZR0PY40wndlXgh9BKaXK2vHj1mjzm2/y2my1abXEXtXl7sjzT8DvhbTXwLpO7hZjzDtAf+AWEdnmxlvWADcZY4Kc2m4FDgJ73f1cpZQqDXu1oBuMYX9YmGvi7NJFa9Mqt5PneqzRp5199PkIsNqdExhj3gOGAP2Ak8aYMNsj1KnP68aYpU5vmw2kAR8bY642xvQBngN0pa1SqkwkJCQw/OGHuT0lhVVA4xyn8cG4cdZ1T52mrfLcnbZ9HvjGGPNn23vG2J63BW528xyP2v67NF/7ROAl2/NGwOX2AyJyyhhzK/AekIi1MncKMNXNz1RKqVJ5ddw4PkhP5wGntpPA2Pr1+ei113wVlipn3EqeIrLaGNMeq0DBLqziGRuB9iLys5vnKLEqkIgMLqTtZ9xP0EopdeG2beOzffv4s1PTRqzbDPYeP85HPgpLlT8lJk9jTAAwC3heRAaVfUhKKeUD8+bBQw+5JM4PgSeATCBKb39TTkq85ikiWcBfKOTWEKWUqvCysuCpp6BvXzh7FrDuxxuCddN4JlotSBXk7oKh+UCfsgxEKaW87sABa9Xs22/ntV1+Ocvj4lgeFYXRakGqCO4uGEoFxhtjbsJauHPO+aCI6AIepVTFsnQp9OvnWi3orrtgxgxur12bvc8/77vYVLnnbvIcjLXgrKXt4UzQ1a9KqYoiNxfeeANeeMF6Dta2Ya+/DmPHgu54qNzg7mrbS8s6EKWUKnMnT8KDD8KCBXltDRvC3LnQqVPR71MqH2/u56mUUl5nrxYUYwz76td3TZw33QSbNmniVKXm7q4q04o7LiKjPBOOUkp5jr1aUP/0dN4FgpyrBY0da+3HGRDgs/hUxeXuNc9r8r0OAJrb3r/RoxEppZSHvDJuHO+lpzPYqe0U8Gz9+nwwaZKPolKVgbvXPLvkb7MVa/8IWOXpoJRS6qLt3Mncffu41qlpC3A3sOv4cT7wUViqcrjga54ikgHEAbGeC0cppTzgv/+FNm1cEufHQDtgJxCp1YLURbrYBUP1gdASeymllDdkZ8Mzz0Dv3nD6NAAZwMNYFYPS0WpByjPcXTA0Jn8T1g4oA4CFng5KKaVK7dAhuP9++P77vLboaJYPHcqS//s/TGoqkZGRxMXFabUgddGMO9tiGmP25GvKBY4By4DXReRMGcR2UWJiYiQxMdHXYSilvGHlSrjvPjhyJK+tRw/45BOoU8d3calyxxizQURiLvY8WiRBKVVxicDkydYm1fbbUPz84JVX4LnnrOdKlQF3b1UpwBjTFNhvWziklFLedeoUDB5sLQ6yq18f5syBrl19FpaqGtz6s8wY85oxZpDtuTHGfAfsAA4ZY24oywCVUsrOXi2olTHsqVfPNXHeeCNs3KiJU3mFu3MaA4Dttud/Ba7FWvX9CfBGGcSllFIuEhISGD58OJ1TUlgDXJqdnXdw9GhYsQIaN/ZVeKqKcTd5NgT2257fDvxHRH4C3gVal0VgSinl7OXnn+fttDQ+BoJtbWeAx+rVg7fe0jJ7yqvcTZ6/A1G253/BWmUL1jVT3b9HKVW2du9mTmoqDzs1/QpcD/zr9999FJSqytxNnp8Ds40xS4A6wGJbeyusgh1KKVU2FiyANm24zqkpAbgB61qSVgtSvuBu8hwDTAO2AreKyDlbeyPgX2URmFKqisvJgdhY6NUL/vgDgPPAo8ADwDm0WpDyHXfv88wGphTS/pbHI1JKqaNHoV8/WLYsr61JE5YNHcrCGTO0WpDyOXfL83UCMkRkne31YGAY1mWHp0XkbJlFqJSqWn78Efr2hYMH89puuw1mzaJ7vXrsnTDBd7EpZePutO3bQBiAMeZK4AOs3X3aA7opnlLq4onA229D5855idMYeOkl+PprqFfPl9Ep5cLdCkOXAz/bnt8NLBGRR20FEj4HRpZFcEqpKuL0aRg6FD77LK+tTh2YPdsadSpVzribPAXwtz3vCnxhe34YqOvpoJRSVcjWrdCnD2zfntfWti3Mmwe6klaVU+5O264HXjDGDARuAhbZ2qOxEqhSSpVKQkICQxs04Myf/+yaOB991NpWTBOnKsfcTZ6jse7p/AcQJyK7bO33AqvLIjClVOWVMGsW24YM4cNjx6hhazsH/Pjoo/DeexAY6MvwlCqRW/t5FvlmY4KAHBHJ8lxInqH7eSpVTmVk8EW9evQ+d87RtBe4AzgdFcXevXt9FJiqCjy1n2epNrszxsQYY+4zxlS3NfljXQ9VSqmSHTwInTq5JM6VWGX2fgZSU1N9FZlSpeLulmQNjTHrgJ+A2ViF4gGmUkjxBKWUKiAxEa6/Hn76ydEUD9wKHLe91lJ7qqJwd+T5Fnkra9Oc2udhFYpXSqmizZkDN93kuH8z18+PpwMCeASwX/PRUnuqInE3eXYFYkXkZL72XYD+qaiUKlxurlWftn9/yMiw2mrXxu+bb7huxgyioqIwxhAVFUV8fLyW2lMVhrv3eQZj1WTOrz6Q4e6HGWNuBsYCbYBwYIiIfFxM/2hgTyGH/ioiiwtpV0qVF2fOwMCB8OWXeW3Nm8P//gfNmjEANFmqCsvdkef3wGCn12KM8Qf+BiwtxeeFAr8ATwLppXhfd6wdXOyPZcV3V0r51J49cOONromze3dYuxaaNfNdXEp5iLsjz2eBlcaY64FArEVCfwZqAR3c/TARWQgsBDDGfFyKOH8XES3GoFRFsHIl3H03OG9S/fTT8Oab4O9f9PuUqkDcGnmKyFbgGqyCCN8CQViLhVo7FUwoS/ONMUeNMT8aY+7xwucppUohISGB6OhoHjGGrM6d8xJntWrw8ccwebImTlWplDjyNMYEAD8AD4qIt/cCOot1jfRHIBvrPuq5xphBIjIrf2djzHBgOOiSd6W8JSEhgZEPP0xcejpPOLWn16pF8MKF1vStUpWMWxWGjDFHgY4issNjH2zMWeDx4hYMFfG+f9piaVlcP60wpJR3XNukCVP276ebU9tG4LGICNbs3++rsJQqlLcrDM0EHr7YD/OQdYCuOFCqPEhO5rN8ifM/WLtHrHPezFqpSsbdBUPVgQHGmFuBDVg1nB1EZJSnAytGK+CQFz9PKVWYhQuhXz+Xv2RfBF6xPY/SSyeqEnM3eV6FNRMDcFm+Y27XtjXGhAJNbS/9gEhjTCvghIikGmNeB9qKSFdb/0FYBUg2AblAL+AxrFtklFK+IAJTpsCzz1rPsf6afhCYb+ui1YJUZedW8hSRLh76vBhgudPribbHTKz7SBsBl+d7z3ggCsgBdgAPFbZYSCnlBRkZMGIEzJyZ19akCStHjGBDfDwmNZXIyEji4uK0AIKq1Eq1JZkxJpi85LZLREpT6MCrdMGQUh52+DD07m0VOrC78UaYPx8aNiz6fUqVI15dMGSMCTTGvA2cADYDW4ATxph3bHt6KqUqs40brR1RnBPnkCGwbJkmTlUlubva9l/APcAwrJWuTW3PewP/LJvQlFLlwrx50LEj2G878fODqVPho48gMNC3sSnlI+4uGLoX6CMiS5zadtvu//wceMjjkSmlfCs3FyZOhJdfzmurVQvmzoXbbvNdXEqVA+6OPM8BBwppP0DpCrwrpcq5hIQEroqM5HN/f9fE2awZrFuniVMp3E+e7wITbAuGAMfioRdsx5RSlUBCQgKvDBvGnH37uNup/dDVV1uJ88orfRabUuWJu9O27YBOwAFjzBZb2zW291c3xvzP3lFE7vBsiEopb5n/9NN8n5FBA6e2t4F3T59m15/+5KuwlCp33E2ex7GubTorbJNqpVRFNX06c44coZrt5XlgJDAdMPv2+S4upcohd4skDCnrQJRSPpKdDc88A2+/7UicR4G7sbZTAt2lSKn83B15KqUqoz/+gPvvh2++cTT9bAw9RUi1vdZSe0oV5O6CIaVUZbNjB9xwg0vipHdvtn74ISYqCmMMUVFRxMfHa6k9pfLRkadSVdG330LfvnDqVF7bCy/ASy9xn58f9w0d6rvYlKoANHkqVZWIwLRpMGaMVQQBIDgYZsyA++7zbWxKVSCaPJWqKjIz4bHHrLJ6dhER8OWX0KaN7+JSqgJy+5qnMeavxpgFxpitxpgmtrZhxpiuZReeUsojjh6Frl1dE+cNN8D69Zo4lboA7u6qMgD4D/AbcCkQYDvkDzxbNqEppTxi82ZrR5Qff8xrGzgQVqyARo18FpZSFZm7I89ngYdF5Ckg26l9LdDK41EppTxi5ejRpLVuDanWjSdiDPz979Zm1kG6m6BSF8rd5NkMWFNI+1mgpufCUUp5yvpBg+j0zjuE2Da8Pw3cXa0aCeHhYIxvg1OqgnM3eR4Eriik/WZgl+fCUUpdNBF49VWu/+QTR9NOrALVX2RmEhsb67PQlKos3E2e8cA0Y0wH2+smxphBwN+xNspWSpUHIjB2rHXPps1a4AYg2fY6NTW1sHcqpUrB3dq2fzfG1AKWAEHAciATmCwi75VhfEopd2Vnw/Dh1j2bNkuA3lgb8tppnVqlLl6JI09jzCXGmNuBqUA9oC3WDFB9EXmh2DcrpbwjM9MqcuCUOFOvv56+wcEuiVPr1CrlGSUmTxHJBuYDNUQkTUQSReQnETlb9uEppUp07hz06gXz5+e1DR5M5OrV/OPDD4nSOrVKeZy7FYY2A02BvWUXilKq1E6ehB49YI3TYvjRo2HKFPDzY8CAAZoslSoD7i4YegmYYoy5yxjTxBhTx/lRhvEppYpy+DB06uSaOF9+GaZOBT/dMEmpsuTuyPNr23/nA+LUbmyv/T0ZlFKqBHv3QrdusMvpTrFp0+CJJ3wWklJVibvJs0uZRqGUct/WrXDrrXDwoPXa3x+mT4cHH/RtXEpVIe7eqrKyrANRSrkhMRG6d4fff7deBwbC3Llw552+jUupKsbtLcmMMQ2Bx4AWWFO1vwL/EpEjZRSbUsrJkthY2r/+OqG2cntZQUEELFwIXXRiSClvc3dXlQ5YFb76A+lABvAA8Jsxpn3ZhaeUAljx9NPc9NprjsT5O9BVhAT71K1SyquMiJTcyZg1wM/ACBHJtbX5Ae8DV4vIjWUa5QWIiYmRxMREX4eh1MWbNYvsgQMd00QHgVuBrUBUVBR79+71WWhKVTTGmA0iEnOx53F3PXsrYIo9cQLYnk8FWl9sEEqpIrz3Hjglzl1AB6zECVqnVilfcTd5nsLaBDu/S4E/PBeOUgpw7IzC4487mn4GOuJaqUTr1CrlG+4mz38DHxljBhhjLjXGRBtjHgA+BOaUXXhKVUGF7Ixy/PLL+WtwMIedummdWqV8x93k+SzwGTAda+HQLuD/gHnAc2UTmlJVUE4ODBtmVQmy69qVeklJvKl1apUqN9xaMOTobEwIcDlWZaGdIpJWqg8z5mZgLNAGCAeGiMjHJbznGuAfWLu5nAA+AF6REgLXBUOqwsnMhAED4PPP89p694Y5c6z7OZVSF80rC4aMMSHGmPeMMQeMMUexRpuHRGRLaROnTSjwC/Ak1i0vxTLG1MTakvAIcD0wCngGGHMBn61U+ZWebhU6cE6cgwfDf/6jiVOpcqikIgkTgcFAAta9nf2AfwH3XsiHichCYCGAMeZjN94yAAgBBolIOvCLMeYqYIwxZmpJo0+lKoRz5+COO2DZsrw2p51RlFLlT0k/mX2AoSIyXERGAT2Au4wx3ioE3x5YZUucdt9gTflGeykGpcrO2bPWlmLOifOll3RnFKXKuZJ+OpsAq+wvROQnIBsreXlDGNaUrbMjTsdcGGOGG2MSjTGJx44dK/PglLooZ87AX/8KK51KR7/2GkyYAMb4Li6lVIlKSp7+wPl8bdmUoiauB+SfmjVFtCMi8SISIyIx9evXL/vIlLpQp07BX/4CP/yQ1zZpEowb57uYlFJuKykJGmCWMSbTqS0I+NAY41gwJCJ3lEVwwGEKjjAb2P6rBelVxfTHHxxv04Z6u3c7mhIfeICYsWN9GJRSqjRKGnnOxCql+bvTYxawL19bWVkD3GSMCXJqu9UW094y/FylysaJE/zeurVL4nwc6DR/PgkJCb6LSylVKqW6z/OiP8yYUKCp7eVq4A3gf8AJEUk1xrwOtBWRrrb+tYDtwArgVeAK4GNgoohMKe6z9D5PVe4cP25tYp2U5GgagXXjMmiRd6W8wduF4T0lBthkewRj3QqzCXjZdrwRVhEGAETkFNZIMxxIBN4DpmAVpFeq4jh6FG65xZE4c4Fh5CVO0CLvSlUk3lz4g4isIG/BT2HHBxfS9jNwc9lFpVQZO3LESpxbrb1QcoEhwCf5ummRd6UqDr2RTKmydOgQdO7sSJz4+bFmxAg+Cwlx6aZF3pWqWDR5KlVW9u+HTp1g2zbrtb8/JCTQ4V//Ij4+Xou8K1WBeXXBkDfpgiHlU6mp0KUL2FfVXnKJVeD9nnt8G5dSVZynFgx59ZqnUlXC3r1W4rSvnA0IgLlzrR1SlFKVgiZPpTxp924rcdpXzlarBp99Br16+TYupZRHafJUylN++81aVbt/v/U6MBDmz4fbb/dtXEopj9PkqZQnbN9ujTgPHbJeBwXBl19a9WuVUpWOJk+lLtbOndbtKIcPW6+Dg+Grr6BrV5+GpZQqO5o8lboY+/ZZSdKeOKtXhwULrGSqlKq09D5PpS7UkSPQrZtjcVAa0Dc0lIQDB3wbl1KqzOnIU6kLcfKkdT1zxw7A2vS2N/DtkSN8PXw4gBY9UKoS05GnUqV19qy1gnbLFgBygPuBb22H09LSiI2N9VV0Sikv0OSpVGlkZMCdd8LatY6mIcAX+brpDilKVW6aPJVyV1YW9O0Ly5Y5ml6oU4dPC+mqO6QoVblp8lTKHTk5MGiQdQuK3euv03zaNEJ0hxSlqhxdMKRUSURg5EirsLvduHHw3HPYlwTFxsaSmppKZGQkcXFxulhIqUpOd1VRqjgiMHYsTJ2a1/bYY/Duu2CK3NddKVVOeWpXFZ22Vao4r7zimjgffBCmTdPEqVQVp8lTqaK8/TZMmJD3uk8f+Ogj8NMfG6WqOv0toFRhfvgBxozJe33bbTB7trWptVKqytPkqVR+GRkwdKh1vROgQwdra7HAQN/GpZQqNzR5KpXfxImOsnvUqGGtss13O4pSqmrT5KmUs40bYdKkvNd//zs0aeK7eJRS5ZImT6XssrKs6dqcHOt1p05gK/KulFLONHkqZTdpEiQlWc+DguDDD3VlrVKqUPqbQSmA5GTrWqfdK69As2a+i0cpVa5p8lQqJ8earj1/3nodEwOjR/s2JqVUuabJU6n33oM1awDI9ffntgMH8KtWjejoaBISEnwcnFKqPNI7vlXVtnevVeTd5g1j+PbQIQBSUlIYblswpIXelVLOdOSpqrZHH4W0NAC2BwQwMTvb5XBaWhqxsbG+iEwpVY5p8lRV17ffwqJF1nNjGJSVxflCuqWmpno1LKVU+afJU1VNOTnwzDN5r4cO5XBUVKFdIyMjvRSUUqqi0OSpqqZPP4UtW6znISHw8svExcURkq8MX0hICHFxcT4IUClVnmnyVFVPWho4X8ccOxYaNWLAgAHEx8cTFRWFMYaoqCji4+N1sZBSqgAj9p0jKpmYmBhJTEz0dRiqPIqLg/HjrecNG8LOnRAa6tuYlFJeYYzZICIxF3ser488jTGPGmP2GGMyjDEbjDE3FdM32hgjhTy6ezNmVYkcPQpvvpn3+uWXNXEqpUrNq8nTGHMf8A7wGtAaWA0sMsaUtCKjO9DI6bGsLONUldjEiXDmjPX8qqvgoYd8G49SqkLy9shzDPCxiHwoIski8gRwCBhZwvt+F5HDTo/C7ihQqnj798MHH+S9/vvf4RKtE6KUKj2vJU9jTDWgDfBtvkPfAjeW8Pb5xpijxpgfjTH3lEmAqvKbNStvu7GOHaFHD9/Go5SqsLw58qwH+ANH8rUfAcKKeM9ZYCzQF7gdWArMNcY8UFhnY8xwY0yiMSbx2LFjnolaVQ4i1u0pdg8/DMb4Lh6lVIXmizmr/Mt7TSFtVkeR48AUp6ZEY0w94FlgViH944F4sFbbeiRaVTls2gRbt1rPQ0KgTx/fxqOUqtC8OfI8DuRQcJTZgIKj0eKsA3SjRVU6zqPO3r11ha1S6qJ4LXnaFvlsAG7Nd+hWrFW37mqFtchIKfdkZ8OcOXmvBw70XSxKqUrB29O2U4FPjTE/AT8CI4Bw4H0AY8zrQFsR6Wp7PQjIAjYBuUAv4DHgb16OW1VkS5bAEdvkRlgYdO3q23iUUhWeV5OniMw1xtQFxmPdr/kLcLuIpNi6NAIuz/e28UAU1pTvDuAhESlwvVOpIjlP2fbvr7enKKUumpbnU5XbmTNWCb70dOv1pk3QqpVvY1JK+UyFLc+nlFetWJGXOP/8Z7j2Wp+Go5SqHDR5qspt1aq85927672dSimP0OSpKrcffsh73rGj7+JQSlUqmjxV5ZWeDs7XvTt08F0sSqlKRZOnqrx++gmysqznzZtD/fq+jUcpVWlo8lSVl/P1zpuK3DZWKaVKTW94U5VLTo614fWBA7BwYV67Xu9USnmQJk9VcWRlwd69VmI8cMDan9P+3P44dChv2zFnOvJUSnmQJk9VPonArl3WdUv7Y+NGyMws/bmuugqioz0eolKq6tLkqcqHY8dcE+VPP8GJExd2rrp1ISKCA8bw3e7dvJGcTPqllxIXF8eAAQM8G7dSqkrS5Km8Ly3NGkXak+S6ddZ0rBv2AYcDA2nQqhVRN94IERGuj/BwCAoiISGB4cOHk5aWZr0xJYXhw4cDaAJVSl00rW2rykZOjnVNcudO18dvv5GbnIxfbm7J56hTB9q2ZUtwMC8tXMiqzEyO2w6FhIQQHx9fZCKMjo4mJSWlQHtUVBR73UzUSqnKx1O1bTV5qguXnQ2pqfDbbwWT5O7dcP6826fKCQjAPyYGbrgB2ra1HpddBsZcUCL08/OjsP+3jTHkupO4lVKVkqeSp07bqpL98QesXQs7drgmyD17rARaSrnAVuAnp8fpRo3YubrwPdFTU1NL1Q4QGRlZaMKNjIwsdbxKKZWfJk9V0MmTVoGBFStg5UprG68LmaFo2BCaNnV5xPTrx3bgbL6uZt++Ik9zIYkwLi7O9Zon1lRvXFxcab8KpZQqQJOnsla12pPlihWwebPbyfKgMfwmwk5gJ7CvWjX6T5jA7U88ATVqFOh//LnnOOuFRGi/FhobG0tqaiqRkZG62lYp5TF6zbMq+v13+P57a1S5YgVs2VJ8svTzg9at4brrXEaRV/XsybZCRozFXYsssAqWkhf/2N+niVApdbF0wVAJNHk6OX7cSpb2adgtW4rtngNsAH685BJajR5Nl/HjoVatAv0udFGOJkKllK9o8ixBlU6ex47lJcsVK+CXX4rv7+/PJn9/vj1/nhXAj8AZ26HiRpF6O4hSqqLR1bYqz7lz1jXL776zHps3F9/f3x+uvx46d4ZOnaBDB9rUqkVhf0YVt6JVF+UopaoqTZ4VUVYWrF8PS5dayXLNmrx9KwvrDmzw8yO0Rw+ufuwxa1Po0FCXPheyolUX5Silqiqdtq0IRGDr1rxkuWIFnDlTZPdsrHsnlwMrgNVAGmWzkEcppSoSnbat7PbvtxLl0qWkf/UVwadOFd//6quhWzfo1o06PXtSWGotbgpWR5FKKeU+TZ7lxcmT1ojSPrrcvt1xKLiQ7ufq1qX6HXdYCfOWWyAszHGsTlQUZy6gus6AAQM0WSqllBs0efpKRgasXu0YXZKYCMXc3nESWAZ8Z3tkVa/O3unTC+2rC3mUUqpsafL0lpwcSEpi06RJnP3yS2IyMgodUToEBkLHjoxbupQlwCasmrB2xZWz0ylYpZQqW7pgqKyIwK5debePLF9e/ObOxkCbNo7rltx4IwQH672USinlQbpgqDw6cgSWLctLmMUs0AHYDiwFNtevzwfbtln7V+ajU7BKKVX+aPK8QAkJCcSNG8el+/bRu0YN7q5dmz8VM5UKcBjreuVS28Pe2xw/zgeFJE7QKVillCqPdNq2NLKyYN06trz1Fmf++1/a5uYSUFz/0FCrik+3btC1K9E9epBSyGhUp2CVUso7dNrWG0SsurD220dWroSzZ2lZRPfzQLWbbsq7bnn99RCQl17jXntNp2CVUqoS0OSJ6y4f7cLDmXL77bQ/d85KmkeOFPvezeTdPvIDcOb774vsq1OwSilVOVT5adt5H3zAF6NG0fH8eboBV5T0hqgo/v377/z37FmWAcdcDun0q1JKlWc6bXuh0tPhxx8dxQnuTkzk3uL616kDXbtaj27d4LLLyJk9m690+lUppaosrydPY8yjwDNAI+BXYLSIrCqm/zXAPSz4WgAACKBJREFUP4C2wAngA+AVcXfInJMDGzfmVfL54QfIzHQc9svXPR34Hquaz5sbNkCrVuDn2kunX5VSqmrz6rStMeY+YBbwKNYlwkeBIUALESmwDNUYUxPYgZXPXgauBD4GXhKRKcV9VgNjZGZICLcAgU4jxPxygPXk3UKyBshEp2CVUqoyqqjTtmOAj0XkQ9vrJ4wx3YGRwLhC+g8AQoBBIpIO/GKMuQoYY4yZWtzoMxL4a1FJs3lzx+0jnx89ypCnntIpWKWUUm7zWvI0xlQD2gCT8x36FrixiLe1B1bZEqfdN8ArQDSwx53PPuLvT8P+/R0Jk4gIx7G+WEXWdQpWKaWUu7w2bWuMCQcOAJ1E5Hun9heBASJyZSHv+RbYLyIPObVFAinAjSKyJl//4cBw28s2+U63wSNfSMVTDzju6yB8TL8H+j0A/R6Afg8ArhSRGhd7El+sts2frU0hbSX1L6wdEYkH4gGMMYmemNeu6PT7oN8D0O8B6PcA9HsA1vfAE+fJv9i0LB3HWp8Tlq+9AVBUJYLDRfSnmPcopZRSZcpryVNEzmNNnd6a79CtwOoi3rYGuMkYE5Sv/0Fgr6djVEoppdzhzZEnwFRgsDFm2P+3d++xcpRlHMe/P48RSKrVBgmhBvqPYOVUiyn/kBYOjdwEYygYtSFQETTUoEEJSmKQQKAQCKZ/YIgNUBLuBKIYLkohgaQxYFUujQoarYoaKkisR1oK8vDHOxvHZW+zl7M77/l9kk3a2fednXnOs/vs7M4+I2mppI3AQcANAJI2SHq0NP524DVgs6RpSWuAbwMdz7Qt/GAE219HjoNjAI4BOAbgGMCQYjDn7fmKJgkXkZokbAcuaJxAJGkzMBMRS0rjlwHXk5okvEoqtJf13CTBzMxsyLLtbWtmZjYqc/2xrZmZWe3VtnhKWi/pj5L2SPqFpFVdxi+T9Lik3ZL+KukSSeo0Z9JViYGkJZKixe3EudzmYZJ0tKT7i79nSFrXw5ys8qBqDDLNg4sl/VzSLkn/kPRjSdM9zMsmF/qJQW65IOmrkp4tYrBL0s8kndxlTt85UMviWfTI3QhcCRxBOlv3oaKBQqvx7wMeIf285Ujga6Tm9N+Ykw0egaoxKDmR9H1z4/bYKLdzxBaQvjf/Oqmnf0c55gEVY1CSUx7MAN8ndSpbDbwJbJG0qN2EDHNhhooxKMklF14EvgV8AlhB2o8fSvpYq8ED50BE1O4GPAlsalr2O2BDm/HnAbuA/UrLvkPqeKRx788cxWAJqbHEinFv+4jiMQus6zImuzzoIwZZ50GxjwtIvyn/9DzOhV5iMB9y4Z/AV0aRA7U78iz1yP1p01399Mg9iJRAtdJnDBruk7RT0lZJp49kAydXVnkwoJzz4L2kT9Ve7TAm91zoJQYN2eWCpClJnye9iWjXR2CgHKhd8ST1ZpzinR2GXuKd3YgaDmwzvnFf3fQTg1ngQlIv/E+RrsB2l6QzRrWREyi3POjHfMiDjcDTpCYr7eSeC73EILtcKL7DnCVdWfIG4NSIeK7N8IFyYBy9bYdlZD1ya6TnGETEy0D5GqjbJO1P+s3traPZvImUYx70LPc8kHQdsBJYGRH/7TI8y1zoNQaZ5sLzwHLg/cBpwC2SZiJie5vxfedAHY883SO3vxi08iTw4WFtVA3klgfDkkUeSPoe8AVgdUT8ocvwLHOhYgxaqXUuRMTeiPh9RGyLiItJR98XtBk+UA7UrniGe+T2G4NWlgN/H9Z21UBWeTBEtc8DpVafa0lF47c9TMkuF/qIQSu1z4Um7wL2aXPfYDkw7rOh+jyD6nPAXuAcYCnp8/1Z4JDi/g3Ao6XxC0nvMu4EpoE1pLOsvjnufZnDGJxFemItBQ4jfdexl9Qecez702cMFpCe7MtJPZAvKf598DzKg6oxyDEPri/+jqtJRxKN24LSmKxzoc8YZJULwFXAKtLJPsuK/X0LOGkUOTD2HR4gUOtJ7w5eJx2FHV26bzOwo2n8MuAJYA/pndV3qfkp6VViUDxRfg38p0iQbcAZ496HAfd/hvTdRPNt83zJg6oxyDQPWu1/AJeWxmSdC/3EILdcKPbvT8Xr4U5gC3DCqHLAvW3NzMwqqt13nmZmZuPm4mlmZlaRi6eZmVlFLp5mZmYVuXiamZlV5OJpZmZWkYunWYYkHSPpBUlTA67nKUlrhrVdZrlw8TTL0zXAFdG9OXo3lwNXS/JrhVmJnxBmmZF0FPAR4J4hrO5B0rUhTxrCusyy4eJpNqEknSnpFUn7NC2/TdL9HaauBbZExGulOZdK2i7pLEk7JM1KulnSeyStl/SX4rGuKx9lFkeuD5Ku1GFmBRdPs8l1D+k5+pnGAkkLgVOBGzvMW0XqU9psSbGuU0jXOvws8CPgSOB40kUGzi/WX/YUcEw/O2CWqzpfDNssaxGxW9JtwNnA3cXitaQm3g90mHoIrS8rNQV8MSL+BWyX9DCpKC6OdJm730jaChwL3Fua9zdgsaR3R8SbA+2UWSZ85Gk22TYBx0n6UPH/s4FbuhSx/UhXiWj256JwNrwEvFAUzvKyA/5/GrsBAftiZoCLp9lEi4hngF8C6yRNAyuAm7pMexn4QIvlbzSvvs2y5teFRcCeiJjtaaPN5gF/bGs2+TYBFwH7A1sj4vku438FfHSIjz9NKuBmVvCRp9nkuwM4EDiPzicKNfwEWDnEx18FPDzE9ZnVnoun2YSLiH+TThjay/9OHOrkVuBQSYcP+tiSFgNHATcPui6znCgixr0NZtaFpIeAFyPi3B7HXwV8MCK+NODjXgMsjIgvD7Ies9z4O0+zCSZpEfBJ0u8wP15h6pXA+ZKmBmzRtxO4doD5ZlnykafZBJO0g3S26xURcfWYN8fMCi6eZmZmFfmEITMzs4pcPM3MzCpy8TQzM6vIxdPMzKwiF08zM7OKXDzNzMwqehvXZ5pDajoZmgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=[7.,5])\n",
    "\n",
    "plt.subplot(111)\n",
    "plt.plot(y_TFVM,p_TFVM,'ko',label='This model')\n",
    "\n",
    "y_anal = np.linspace(-Ly/2,Ly/2,1000)\n",
    "#plt.plot(y_anal, analytic2D(ell, 0, y_anal) ,'r-',linewidth=3,label='Analytic')\n",
    "plt.plot(df['X'],df['Y'],'r-',linewidth=3,label='Pouya and Ghabezloo (2010)')\n",
    "\n",
    "plt.xlabel('y (m)')\n",
    "plt.ylabel('Pore pressure')\n",
    "plt.xlim(0,3)\n",
    "plt.ylim(0,3)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-04-20T04:03:38.438093Z",
     "iopub.status.busy": "2023-04-20T04:03:38.437309Z",
     "iopub.status.idle": "2023-04-20T04:03:38.440706Z",
     "shell.execute_reply": "2023-04-20T04:03:38.440166Z"
    },
    "papermill": {
     "duration": 0.025634,
     "end_time": "2023-04-20T04:03:38.440839",
     "exception": false,
     "start_time": "2023-04-20T04:03:38.415205",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Geometry\n",
    "dx = Lx/(nx-1)\n",
    "dy = Ly/(ny-1)\n",
    "xb, yb = np.mgrid[-(Lx/2.+dx/2):(Lx/2.+dx/2):(nx+1)*1j, -(Ly/2.+dy/2):(Ly/2.+dy/2):(ny+1)*1j] # this is coordinates of cell boundaries\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-04-20T04:03:38.484745Z",
     "iopub.status.busy": "2023-04-20T04:03:38.479891Z",
     "iopub.status.idle": "2023-04-20T04:03:39.158393Z",
     "shell.execute_reply": "2023-04-20T04:03:39.157706Z"
    },
    "papermill": {
     "duration": 0.702326,
     "end_time": "2023-04-20T04:03:39.158520",
     "exception": false,
     "start_time": "2023-04-20T04:03:38.456194",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUEAAAEzCAYAAACv5LH7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOy9fXAcyXUn+KuP7uovdANoNNj4IIAGGmCTmObUYKChiJOl8drS+oNn2eexfJLPlvZsiTeOC0leR6zlOcdKjtBKO9aGjtau7Butb+4knkIT5zn7rJVlyzp7aXl5GkujYd/MQKQw5HAwJEGQIAgSbHz0Z90f4Eu+yspq8AP0EGS9CASbVe+X+fK9V1lZmS9fao7jIKCAAgroQSX9rRYgoIACCuitpKATDCiggB5oCjrBgAIK6IGmoBMMKKCAHmgKOsGAAgrogaagEwwooIAeaDLfagE4dXV1OUNDQ2+1GAEFFNB9Rj/4wQ8uOY6TUd27pzrBoaEhvPjii2+1GAEFFNB9RpqmzfrdCz6HAwoooAeagk4woIACeqAp6AQDCiigB5qCTjCggAJ6oCnoBAMKKKAHmoJOMKCAAnqgKegEAwoooAeagk4woIACeqAp6AQDCiigB5ruqR0jt0qPPPIIFhYWxP/z+TxOnjzp4gmu3bvX7gUZgmt3du2tkCGTyeDYsWPYKtrWneDCwgLK5bLrmvz/4Nq9fe1ekCG4dmfX3ioZtoq29edwPp+HbdvKe5lMBrZt45vf/KbyfjabhW3bGB8fV97v7e3FkSNHUCgUlPf7+/th2zbGxsaU93/v934Ptm0jn88r7+dyOdi2jc985jPK+9S2wcFB5f2xsTEcOXIE/f39yvuFQgG2baO3t1d5//nnn4dt28hms8r7xWIRtm3jfe97n/K+bduwbRvpdNpzT9d12LaNI0eOoKOjw3M/FAqJ+8lk0nPfsixRfiKR8NyPxWI4cuQIbNtGLBbz3E8kEgKv0s/+/fvF/VAo5Lnf0dEh5NN17yOSTqcFXkWZTAZHjhxBsVhU3t/M95588knYtr2p7z3zzDPK+4ODg5v63pEjR5DL5ZT3N/O9Z555xle3AHD48OGWvjc+Pg7btvGhD31IeZ98L5Px5jto1a7bJe1eOmhpcnLSuZUECv39/Zifn0ej0QAAGIYBAGg0GtA0DbquIx6PY2VlRfAQXzgcRrVaFfwcC2w8yG1tbWg2m7h27ZoLCwCO40DTNPEvxwJAPB6HaZq4du0ams2mRz5d16FpGiKRCNbX1z3yRaNRrK2t+ZafSCRgGIZv+fSb2imXT3pxHAfNZtPTfsMwkEgkUKlUsL6+7mk/UbPZFB0FryOVSgEAlpeXQT4my5dIJLC+vo5KpeIqX9d10SZV+Zqmic5TVT7HUFlcNsuyEIlEUC6XW/pOIpHAysoK6vW6S75QKIRarebRNZUViURgWZZv+WT7WCyG1dVVj20syxI6uRd9k+T2K7+trQ26rt+xbzabTTiO48Fms1mcPXsWt0Kapv3AcZxJ1b1tPRIMKKCAArpTCjrBgAIK6IGmoBMMKKCAHmja1p1gPp8Xc2M0x0T/tywLiUQCDz/8sIuH+AYGBpBIJBCPxz1YwzAQiURg2zb27NnjwSYSCUQiEYGXsYZhYGRkREzcq+SLRqNIJBIYGxtTyrdr1y4kEgnBJ/Ps3r0btm0jEokoy4/FYkgkEhgaGlKW/9BDDyGRSMCyLGX74/E4bNvGzp07le2nP9Izx5qmKRYOwuGwB2uapli86O7u9pTf0dEheEOhkKd8WlihhQ25fI6hsngd3d3dYtHFNE1P+8PhsJCvra3NI19vb6+LX5aPFi7It/x8s1gsKm0zODiIRCIhfMfPNwuFgtI25DPkA3Ido6OjsG0b0Wi0pe+MjIwo5SsUCq5nQOahhQ8/36Rnhp5Bufy9e/f6+mYikdjyhZFt3QmePHkS1WoV4XAY4XBYrHSGw2FomoZqtYqLFy+Ka5xvfn4e1WoVtVrNgw2Hw3AcB/Pz8zhx4oQHC2xMPlerVdTrdQ82HA6jUqmgVCqh0Wgo5Ws2m6hWq1hdXVXKR22jieVsNuvimZmZwfz8PBzHEdc4X71eR7VaxdzcnAcLAIuLi2JhSMaGw2HUajWUSiXouu7BZrNZVKtVVKtVMTHOsaFQCKVSCfPz8y488YVCIVSrVczPz2N5ednT9kqlIsqniXCuI8MwMD8/j1KpJDotrl/DMASeVqd5HbFYDKVSCdVqFaFQyGMbXdeFfM1m0yMf6Y58T5ZvYWEB8/Pzwrdk/QJAtVrF9PS00jbkG/V6XWkbx3FQKpU8uiM+wtOCg1zH6dOnPW3jfI1GA9VqFbOzs0r5lpeXUa1W4TiO0veXlpZQKpVcvsn5arWa8D2V77/yyivCt2RstVr1xBHeKW3rOEFgQzEUQzQ/Py9+h0IhWJaFTCbjURp1gNQZra6uurBUbjabxdzcnMtgxGdZFkKhEBqNhgcLbKxAAhsP5Nramke+SCQC0zQRjUY9WFrxps6sXC4LmYmoczhz5oxYSaQHt1qtIhaLuR5mjq1Wq0in0+I6/XG+eDwOAOKB4Nj5+XnxQFQqFY983HnfeOMNsbpK2Fqthng87rlPbV9dXRWhL2trax796rouXmT1el28KIiPRtoAsLS05NEvf/HQCjmvwzRNRCIRZLNZnDp1CjyCYn5+HpVKRdiX2s7roM7h/PnzWFlZ8diGHvhms6m0Db04yDdl21DdyWRSaRtd12GaJhqNhtJ3SH9nz54VK/+8Dhoh1ut1pXzJZFJ0Zirfb29vx/z8vODntgEgfJP0x4m/2CuViqcOVcjUndK27wRpGZ1+AxtOSKEfFMrAHZmHXxCeY+nfVljCqLAAxINNsvjJx/lU8nGsHM5Uq9Vc5Wua5pFJhdU0DbVaTVynUAe5PsMw0Gg0PFhe9p3It5ltuJ45n2wblX5V5RDRCIm3/3Z8x08+ajvXCW+7LJ+sX+rY6Z5sG/qX9Cvbhvvm7dhmM/l42/zafjO28fN9P1lkv94qCuIEcW/GYgVxgkGcoF/5b7VvBnGC9xAFCyPBwkiwMBIsjNwp3RcjQSKaL+DzDPfDtXtFjqBdD3a77pW2bvVIcNvPCQI3huKkJP5J4ziO53O4XC6LiX3i51gA4nNK/pQiPv7JIWOBjfkO0zRd8xicjz456vW655OjXC6LiX2/8unTwK98+oSgRQ+5fFr4ILxcB+ErlYqy/UQqrHyfXrQq+VRt5/t1VeXTJ5hf+fxzeHV11WPD9fV1RCIRV5kq3wGg9B2+31jlO3zhRFU+2Z4WPuTyCetX/lvtm7Ro5Vc+6f9OfZNsq8JuJW3rz+GAAgoooDulbd8JUkwghWUAGyMF0zQRDofR09Mjrsl8FNMmYzVNg2EYIguGCkvlU5yZzEfZVah8uQ7CdXV1bSqf3E4aCfX29sI0Tdc14uNyqbAU/2UYhgfLZUilUh4sj+syTdODpZFKb2+vKJ9jKQSklX6pfOLlfLquCyyNWjiWxybSiJrXQZ9XVL5sG5oXvBn5VL5Dbec64Tqi8v1s09HR4fJNmY90nslklLYh22/mO9w3b8V3MpmMsL3K93fs2AEALturfL+9vb2lfsm3OZZ/wW0V3RedoK7rrodF13XhaJSOh3g4HzmKjKU4q+7ubl+sqhPkfGRgKl+ug/AUzOsnH+9kOA8AdHd3wzRN1zXikx1ZxnZ1dYnyZSzhgRsrfXId9MdX+lTy0WqvXEc4HG6pX16+Sr+tsIQJh8Pis5fz0efcZr5zM/KpfIfazm3P6+Odh0pvFIdHeJmPfCKdTittw33zbvhOOp32dIKchwYAvHzOR+XTCr+ffsl35JfbVtO27wSbzSZisRhisZhYtYrFYgiFQiLWi65xPgpGBeDBxmIxmKaJ9fV18dDIdRiGIWKaZGwsFnONQFTyUegGrRTKdVD5FKpRKBRcPIZhYH19Xazm0Yoe8ZFuwuGwB5tIJNBoNNBsNmGapgcbi8XEPFs6nfZgC4WCiGOjlViOpbavr68jHA57sNFoFM1mU8zNyW2PxWKifMuyPPql0A0ALjzxWZYl8LlczqNfGuE1m01Eo1GPbcg31tfXEY/HPfLx8lW+Y1kW1tfX4TiO0jamaQrbqmyTTCZdvinrl148sm8SH/kWjXrlOlS+zflovtI0TaV8VL5hGErfJyLfknVEum9vb1f6PtUfCoU82LsRJ7itO8F8Pg9d17G6uorV1VWcOHECwMZkeK1Wc705iIf4stms+BSSsTRZHYlEUCgUPFjgxgQyAA92dXUVa2trsG0bjuMo5aOHgGIV5TqobWT0EydOuHjGxsYQiUTQaDTENc4HQHw2yljgRjxevV73YFdXV6FpGmzbxuXLlz3YEydOiDc3xctxLLU9EomIrYEcu7a2Bl3XEYlEkEwmPW23LEuUX6lUPPqtVCoiTKRSqXj0W6lUBP706dMe/c7NzcG2bei6jrW1NY9taNdGJBJx+Rfx8RGYyncymQwikYjwLVm/9Xoduq5jfHxcaZvl5WXhWyrbNBoN2LYtFlbkOvjChMp3crkcIpGIC8/5HMeBrusYHBxUykflNxoNpe+T73Df5Hw0pXHlyhWl7xeLReFbMlbX9SBEhlMQLB0ESwfB0kGw9M3QfRssHVBAAQV0pxR0ggEFFNADTXc9WFrTtKcA/BsAX3Qc53+8G3XQcFkVbU6Tw6oIdPocpk8zGQvcWC2UscTHQy44llbwKCTCTz7aQqaSr1X5VOZm5dMqZqvyW7Xfsizx6aeqgz4bVVi674cFNj4d+efmVpdPE/FynUStdAdsfDqqdEc6UemOBzu3Kp/apiqfPoe3o2+STe7UN6mcVvbfCrqrI0FN094O4MMAXr4b5efzebFaZ1mWOP3Ksix0d3ejUCjgT/7kT8Q1zpfNZlEoFDA6OurBWpaFnp4eHD582FU+5+vr60OhUMDw8LAHa1kWfvu3fxu2bWNoaEgp38DAAAqFAp566imlfNS2nTt3AtiY/+Q8Y2NjOHz4MPr6+sQ1zjcyMoJCoYDe3l4PFgC++MUvolAoIJvNerBUvm3beP/73+/B9vf3o1AooFAoiDhHjqVFi8OHDyOdTnuwtOJ3+PBhJJNJT9tpRbhQKCCVSnn0m0wmxYlmHE98qVRK4OlEOF7HT/7kT4q9t7Say+tIp9NCPl3XPfLRffINWb5MJoPDhw9jbGxMaRvyvfHxcaVtPvKRj6BQKGBkZERpm76+Pti2jc9+9rNK2+zcuROFQgFDQ0NK38nlcjh8+DAGBgaU8g0NDaFQKGBwcFAp32c/+1kUCgX09fUpff/zn/88bNtGT0+P0vdHR0dRKBTwoeunzcn6LRaLKBQK6O7u9mALhcL2WRjRNC0F4CVsdIL/GsCrm40Eb2dhpFwuu/Yx0m/KJ1gsFvHd737XhUskEp58gvIbLxwO47HHHsNLL70kVsV4HTyfoGEYnrfU+Pg4Zmdn0Wg0RD5BXgflE8zlcnjllVc88lG5tHVJli8ajeLRRx/F97//fdfCAvHxfILy2xYA9u3bh+npabF1SZYvHo9D0zRkMhmxwsrrIKKcbxyraZpIwPDCCy94RnorKyst78v5BOWci7quY2pqCqVSSay0c/lpYQnYeHAWFxdd8vf29mJ5eRkAxOIQl5/yCdq2jaNHj7pSPdFiEc8nKNsmFAph3759OHbsmMgnyOWjERAt4Mi2GRsbw9zcnPBN2Tbke2NjY3jppZc8tqFtceR7snyRSASTk5N48cUXXYtexMfzCZqm6ZFvYmICMzMzqNVqCIVCnvuTk5M4ceKEr2+Rb1LyYFl+WnBR+VYisZFMYbvsHf4SgOcdx/k7TdP+9d2qhOe7o1UmWvGkJXq6xjEc6ziOCwvA5YB+uego1o6I85Fz8b2tKvmIz08+Kl/O6wfcCLXhOeGIT65Xzgm3trYmyienk3VpmqbIIMyx8j7R25HvZmxD92TbUNgR/Vbpl36rctbxlV2V/W/Wd/zko7Zznci2oVVolW1ovzbJJvORzil6QLaN7Ju34zu8va18R9V2eum38n3HcUQHKeuXwm+IT6XrLSVS9Fb+YWP09wMA4ev/PwLgP/jwfgTAiwBeHBgYcG6F+vr6HMMwHAAOAMcwDPF/TdMcwzCcZDLp4iG+aDTq4ue/ATi6rjupVMppa2vzYA3DcHRdd/0r1xGPx51UKuXouq6Uj3DxeFwpXyKRaFl+IpFoWT79pnbK5ZNeCC/XYRiGk0qlnEgkomw//ZGe5TpSqZSTSqUcTdN85UulUo5lWZ7yQ6FQy/I1TWtZPsdQWbwOy7KcVCq1qe+kUinHNE2PfJFIpKXvRCKRluWTTdva2pS2icVi97Rvktx+5be1tW2Jb5JtZWxfX98t9RPORkfzol9/teVzgpqm7QLwGQC/4jhOdTN+x3G+5DjOpOM4k6oT5wMKKKCA7ibdjYWR/QC6ALyqaVpd07Q6gHcB+M3r/7daw2+egqSqQVLVIKlqkFT1TmnLF0Y0TWsH0C9d/t8AvIaNEeK041NpkFRVfe1ekSNo14Pdrnulrfd8UlXHca4AuCIJsALgsuM4r251fUCQVDVIquouP0iqGiRVvRXa9jtGaJM7bXgHbqSqikQiGBgYENdkvkgkIlIFySmXTNPE4OCgLzYUCgm8Kl0TpWHiqYg4H+HoTNVW8sntpIdgcHDQk66J+LhcKmx/f78I05GxFBsHAJ2dnR5sJBIRfxR46yefnK6J+CKRSEv9Uvn0+SfztcISJhKJiAeK89G+481852bkU/kOtZ0nguD1Ufl+euvq6nL5psxHnXBPT4/SNmT7zXyH++at+E5PTw8ikYh4BmQd9fX1AXCn0lL5PqXc8tMv4TmWB7pvFf2TdIKO4zzu3KXdIjypqBwJb5qmMnEj8ZmmKfAcS3jK9afCUvn0J/O1tbXdlHzE5ycfvfl4OTz5pmG4k5by9vByZGx7e7soX8ZyGSi1klwH/anko4elo6NDjCo4ljqCVvql8omX8+m67srDKGMJ45f0kx4kKl+2za3Ip/Idaruc8Jb+T7b3sw0lOKB7Mh/pXJXwlmyymf3vxHdSqZSrHllHlDyDl8/5qHx6QbXyfRV2q2lbjwTz+Tzy+Tyi0Sii0ahYKIlGo8hkMsjn8/iDP/gDcY3zDQwMIJ/PY3h42IONRqPIZrM4dOgQ9uzZ48EmEglks1nk83kMDg56sNFoFE8++aRYWFDJ19fXh3w+j4997GNK+Xbt2oV8Pi8m4Xk7o9Eodu/ejUOHDiGbzYprnI/kGhoa8mATiQQ+/elPI5/Po7u724ONRqMYHh6Gbdt44oknPFjSez6fRzqdVspn2zYOHTqEzs5OD5b+PXTokKift72jo0OUn0qlPPpNpVI4dOgQbNtGKpXy6Jcw+Xwejz76qEe/tGOEl8t/d3Z2Cvna2to88vX29rp8T5avv78fhw4dEr4l67e7uxv5fB7FYlFpmw9+8IPCt1S2yWazsG0bTz31lNI2JN/OnTuVthkdHcWhQ4fQ19enlG/nzp3I5/MYGRlRyvfUU08hn8+LZ0DW76c+9SnYtu3yTc43PDyMfD6PD3zgA0rf37t3L/L5PDKZjAdLet9Kui9SaRHdK5PFW33tXpEjaNeD3a57pa33/MLIW0HBwkiwMMLLDxZGgoWRW6Ft/TkcUEABBXSnFHSCAQUU0ANNQScYUEABPdC07ecEedzQ6OgoTp8+7UrNNDw87LlGfHzuSYXN5XI4ffq0mKPw42t1jeb2bgfLQzlyuRxee+01wWMYhueais80Tezateum+ORrAFAoFG4buxkf6ZeHPdysjgjLTze7Xf1uVocKy+emZOzN6kjXdRSLxdvW78jIyG1hb+XaZvKp9DYyMrJlvu/Ht5W0rUeC+Xwek5OTIu6I4pNM00RPTw8mJydx9OhRcY3zZbNZTE5O4uGHH/ZgTdPEwMAAjh07hkKh4MECwNDQECYnJ7Fnzx4P1jRNfO5znxP7O1XyjY6OYnJyEl/84heV8lHbRkZGANyIzaK/sbExHDt2DENDQ664LuJ76KGHMDk5id7eXg8WAL797W9jcnIS/f39HqxpmnjkkUdg2zY+/OEPe7CpVAqTk5OYnJwUB21zbDgchm3bOHbsGDKZjAcbjUYxOTmJY8eOIZlMetpuWZYov7293aPfZDKJY8eOiaSqsn7b29sFfnx83KPfn/qpn4Jt25icnEQ0GvXYJpPJCPl4zCHxpdNpUb7KdzKZDI4dO4ZHHnlEaZv+/n4hm8o2n/jEJzA5OYmHHnpIaZuhoSHYto3nnntOaZuRkRFMTk6iUCgofSeXy+HYsWMYHR1VylcoFDA5OYnBwUGlfM899xwmJycxNDSk9P2vf/3rsG0bAwMDSt9/+OGHMTk5iY997GNK3y8Wi5icnERPT48HOzk5GYTIcOrv78eVK1dE/rJoNCrOezVNE6FQCCMjI5iennblLKNzb2kFa3193YUFNhS+Z88enDp1ypUYk/jIaFQOxwIbo4hLly6hWq2KpKe8Dkos0Nvbi9dff90jH63uNRoN1Go1RCIR0U5gI4g5n8/j+PHjIj8eBQKvra2JIx9JPo7VNA27d+/G6dOnUavVUK/XXVgAYkdAIpHA+fPnXVj+lq7Vamg0Gi6spm2kXc/lcpiennat6tIJZpFIBKOjoy75qe2UCBa4sdLK9atpGorFIk6fPo1yuSyuk341TRMrrOFw2HUiHbCxI6NSqYjVT8dxXLahJAqjo6N49dVXXfn1otGoa5WSkr5y+eis55mZGZEvkuuXfFPXN07Ck23T39+PS5cuCd+UbUM7NXp6ejAzM+OxDcnXbDZRrVY9vkOxgpT4VJaPdroQyfKNjo7izJkzYpVZ9v1CoYBz585hbW1NjN64jsg30+k0zp075/F9wzCE39frdRc2Go2ivb19S0Nktn0nGBy5GRy5KZcfHLkZHLkpU6tOcFt/DgcUUEAB3Slt+4URIDhtrlX5wWlzwWlzMvZu+ybZ5E59k8ppZf+toG09EqQ9hSqi+ay9e/cq72ezWSQSCdfKLyfLssTCht/9RCLhWp3kNDg4KBJXqogSUg4PDyvvU9v8smbQaXD8geNEiTF7e3uV98fHx5FIJFw7ZzhRUlXKhiNTInEjgalMmqaJpKeqDe/08NDChkxyolaZDMMQ5fvdJ7xKP8lkUiRV5btPiOjht23b1ckS0X5pP9/LZDIiqaqKKGkrLdrIRPvF/XyLfHNsbEx5n3zLz/dyuZxIeqoi2iNMWXRkokTAfr5HCX/97lNSVcqgJBMlm1X5Fu1x3koK5gRxb867BHOCwZygX/lvtW8Gc4L3GEWjUZGLjN58lHONMlbQNZkvGo2KlTCOJTy9cVTYcDgs8DJW13UxAuM55Tgf4Sn3mp98NFLj7aQHPJ/Pe3LCER+XS4XN5XKIRqNilVLmo1FCJpPxYHnWD3pbcyw5bT6fh2maHix18q30S+UTL+czDENgqdPkWMLQiEaug9KrUfmybUzTvGn5VL5Dbee5+LiOyDf9bLNjxw6Xb8p85BMDAwNK23DfvBu+MzAw4KpH1hGNIHmuS85H5dNXhp9+Cc+xfqPbO6FtPyfI8+HxeTJ6WOgh4J89tATP851xLOCei5KxdJ9+y1gArgdSJZ9hGOKt2ko+x3FQq9Vc5RAlEgnx1pb1QQ83bycnOv+CRj5+slqW5cHy0aCqDvqX5OMhJly+VvpV6ZmXz7GyzBxDnTSvIxwOY319/bZ9p16vK+dkOR+N1FW2ofJpBC7rl5LJym3jvgncyPXIiXyG28jPd/zkI9+kdspYOvuEf0FwHp7I1s/3+YtWloHkJ7xK11tJ234kWK/XRZZfOqCGnKher2N5eVlc43zhcBj1eh3NZtODpTf40tKSMJZcB2XZcBzHg+VzLRQXJddBsjebTaV8JD8AMT/DeQzDwNLSkvhsoXkg4ms2m6jX6wiHwx4sTXbX63XRPrkO+kSJxWIe7ODgIOr1unhIZCyNQJaWlmCapgdrWRbq9TqWlpZcZVPbY7GYKD8UCnn0Gw6HsbS0BAAiSzHXbygUEnia++V1ULLUer0uPo15HaZpCvni8bhHPpKffE+Wz7IsLC0tifhJWb+6rgvdq2xD8pNvyvqlTqNSqShtQ21zHEfpO6ZpCv2p5HMcR8QAquSjGEtN05S+T3GF/NnhfOSbNPcr65eeLcMwPNi7sWtkW3eC9Lm1vr6O9fV1zM7OAriRnolHmxMP8WWzWZjmRvZaGUuBmR0dHSgUCh4sABGQreu6B0tzaLZtQ9M0pXyAO3O0XAe1jWh2dtbFMzY2ho6ODhFQS2USH33W9fb2erDARsdKwd4ylkZJtm1jbW3Ng52dnRXB4tRZcmylUoFt2+jo6BAByRxbqVRgmhuZm5PJpKftlmWJ8mu1mke/tVoNHR0dsG0btVrNo99arSbwlG+S17G0tCQWbWjOk9dBHQBlxpblS6fTonyV72QyGZG5WWWbZrMJ0zQxPj6utA3JT3XL+nUcRyw8qGxDvkWdpVxHLpcTLwKVfJTBeXBwUCkf2cdxHKXvx2Ix2Lbt8k3OZxiG0L3K94vFovAtGWuaZrAwwilYGAkWRoKFkWBh5Gbovl4YCSiggAK6Ewo6wYACCuiBpqATDCiggB5o2tYhMnSy1vHjxwFsRKrPzMyg0Wigs7MTO3bswJe+9CU8/vjjrnmH3bt3Y3l5GclkEvV6Ha+99poLCwDd3d345je/iQ9+8IOYnp52YWdmZsTkN2UR4VgA+MIXvoAvf/nLuHTpkpjY5XX09/ejra0Nv/Vbv4Unn3zSIx/Na129ehUXLlzA2NiYaCfxHD58GAcOHBCT/4ZhCL7h4WExL5VIJFxYwzDw13/91/joRz+KS5cuYWFhwYUFNnYFtLW14V3vehc+//nPu7B8p8L58+dx7do1F1bXddi2jWeeeQY/8RM/gatXr7qwr732GkZHR3H48GG8+93vFiuV1K7Z2Vmxykm/uX6j0Sj+4R/+AQcPHsQPf/hDsRpJ+rUsS+Cz2Sy+853vuPT7/ve/HzMzM1hdXcVrr72Ger3usk0qlUJvby8OHz6M/fv3u7YN7tE+WzsAACAASURBVN69G4uLi+LM3OPHj3t8J5vN4hvf+AY+8IEPiCwvXL+ZTAZdXV1iblG2zSc/+Ul87WtfQ6VSweuvv+6xTTabRV9fH37nd34H73vf+zy2uXr1KlKpFFZXV3H27FmP7+TzeTz33HP4+Z//eTG3xuugFXrK4CPL97WvfQ2f/OQnsbS0hI6ODo/vf/3rX8enPvUpnDlzBhcvXnTZptFoiBRev/ALv4Cnn37a4/s0J3rhwgVcvnzZhd29e7fvTprbpW3dCZ48eRLlcllMnJJyDMPA1atXcfXqVXz84x8X14hisRhmZmZcnQfHAsDi4iIOHjyIEydOeLAAsLCwgIWFBQAbq8ByHV/+8pdRKpVc13kd58+fx/nz5/HMM88o5SMsv8Z5ZmZmcPDgQdGByXzU8ZJ88oLGJz7xCdcDKtdx6tQpwStjb0a+UqmEgwcPuuxDfI7jCPmXl5c99yuViitFlGybarWKgwcP+uqX42Us6Y7LL9u/XC4L+fiiEZU3MzODxcVFJRbY8I2DBw/i1KlTSttcvnwZly9fBqC2zV/8xV+0tA353tNPP620zczMDC5cuOCpl+j06dM4ePAgzp8/r5SPLzqo5Hv66aeFfBTsz3k+9alPtfT9119/HQDwrW99y4OVfUvW78zMjO92xdulbb86vLS0JFYvI5EIKpWKWFGiGLOZmRlXzjKKhaIVNIq3IixwI/Hpm2++6conSHy0zE8jNo4FNqLql5aWUKlUxCo0ryMUCsEwDOzYsQNvvvmmRz6+qkkxVXyVNhaLYXBwUIxkAIg8euvr6yLImdrJsZQT7s033xTxbhwLbKyghsNhxGIx1wPFc/UBG+EczWbTI19bWxv6+/vxox/9SKwQ8joikQhyuRxOnjzpGmlFIhHUajUR5Ey5Bbl+NU1DoVDA2bNnXfkESb+apoldFaFQyMUDAJ2dnajVaiK8RrYN7crI5XIu+YmPrz5TW7h84XAYIyMjeP3118UqL287hdcYhiHCQLh+e3p6cPnyZeGbsm1oR0t3d7frZUV8PJ9grVbz2CYajWJoaMile14H7SQh/5HlGx4exrlz50Qokez7+XweFy5ccOUT5Doi3+zo6MD8/LzH9ylOk+eqJCzFeQZHbjKqVCrCScnhaGm92Wzi3Llz4hrH0INF9/hvYCPYlBQtY+l+o9FwhQlwvoWFBZimiWq1qpSPEkbS54JcRzQaFYbXNM3VTqKzZ8+KoFpg401JfNTxUjvlXRtzc3OifBrtcD7qfJaWljxYHtLCXwIq+UhHsnyVSsVXv1Sen200TXNhqXweVkK/6SHndVy7dk08WH6+00q+UCjUUj5qO7c9bzv5Do16Zf1evnzZVaZsG+pY5ufnlbahYOfNfIdeYLJ81WpVhMio5Jufnxf10Kcs5zl//rwICFfpl3yTRsOyfqlj5baQdb2VFCyMBBRQQA80betOkNJN0V5Hnn6JtuQ8/PDDLh7iGxgYEKm0ZCxt17FtG3v27PFgaQsP4WWsYRgYGRmBbdtinkWugzb3U1oiuY5du3aJdEiq8ildEW2hk8unVFpDQ0PK8h966CGRDknVfkqltXPnTmX7eaoqGWuapkh1RccIcCxPVdXd3e0pv6Ojw5WqSy4/FAqJ8mlagZfPMVQWr6O7u1uk0qLPUl4HpbqybRttbW0e+SjVFfHL8vX394tUWq18k1JGyeXTNjHyHT/fLBQKStuQz5APyHWMjo6KNG+tfGdkZEQpX6FQcD0DMs/4+HhL36Rnhp5Bufy9e/f6+mYiEaTSchHtGCGiCVNVksbtfO1ekSNo14PdrnulrVu9Y2TbzwkCN7bskJL41ifHcTzb5srlsiuZaKPRcGGBG5k65C1XxMcXVmQsADFpTHNusnw08Ux7a+U6KAuNX/k8y4yqfFpxq1arnq1J5XJZZPogvFwH4SuVirL9RCqsfJ9etCr5VG0n3fuVz7OOqMrnCxerq6seG9JiBi9T5TsAlL7Dk32qfIdCS/zKJ9s3m03PtrlyuexaeLoXfZNWa/3K54eY3Ylvkm1V2K2kbf05DGysklIuMjIOre7FYjHs2rVLXJP5KEOKjKWcb5RVWoWl8iORiAer6zp27twJ4MZKr1yHZVlihbeVfPRA8HbSQ1AoFDw54YiPy6XCjo6OIhaLufIJcj7qhHfs2OHBxmIx8UcvE44lRy0UCq6ccsRHn1yt9Evlm6bp0a9hGAJrGIYHS5hYLCb2GPM6Ojs7XbLItgmFQjctn8p3qO08Fx/XEZXvZ5ve3l6Xb8p8pPNcLqe0DfkWdfR+vsNzXd6K7+RyOWF7le/TMbHc9pyPyu/p6WmpX/Jtjt3qGEHgPhgJ0luLfgM3csxRZ0PXOEbm57/pHjmRH5bedjIWgHBUP/lkh5broDcyr0+Vd06+Tv/fDMsfAC6fLGsoFPJg+UhNhaV/Ke2Tn3yb6ZfrivPRyiX99tMvoM71SNl5NvOd25WP2r6ZbWjEJOuXv5ha2UaV65F3KDRCvl3f8ZOP+46q7fTMbeb7qlyP/JnitpB1vZW07TtBHoZAS+/0mVCtVkVAM1+GJ75qtSqG7BwLwDd8hfgoBIF/knC+K1euAIArDIHXQeEDtFtCrsMwDE/b5DCHixcvusIQdF0XfLVaTTiXjNV1XZyJTHiOJbnD4bArEwivg4jazrFU78WLF9FoNJTyVavVlvrlGX5k21DZdE3WL8dQR8Cxq6urLpkdx1H6zs3Ip/Idko/bnredQlzkECriW15edvmmbBvKs7i4uKi0DenkTnyHXsIq+RYXF8VnrEq/FEjOy+c6It9cXl5W6jcUConyZf1y39squi86QXo78NOrSGHz8/NIJNwnZtHoi59oJp981Ww2MTc353rbyXVQB5BIJDynZtFWMUqeKWPpQbhy5YqvfLQdjCaGuQNQrB85q8xH8XEqbCKRwIULF8Q1Xdc9fBQnSLkJOZbrzq8OAJibm3ONlDlftVrF3NycJ1Gm6gQy1bW5uTnxW6VfHicp65fqo3g4Cq4mLHWKc3NzntPmbsZ3LMvC3Nycyzd528k3CSfbhsuvsg353eLi4m3ZxjRNzM3NuXyT8/FkvjI2kUiIThC4kc6f18kHHirbkG+urq76+j5vv8quW0nbek4wOG0uOG3O7z7hg9PmvBScNuem+yJEJkiqGiRV5eVzDJUVJFUNkqoGSVUDCiiggBQUdIIBBRTQA01bvjCiadrvAvhvAOwCUAHwAoDfdRzn1a2uC4Br3qRYLGJ6etq12vTQQw95rhEfn69SYcfHxzE9PY22traWfK2u8XmhW8Xy+Kjx8XG88sorgudWrj3yyCO3hQWAiYmJ28Zuxkf65fOSN6sjwvJ5udvV72Z1qLB8vkrG3opt9u3bd9v6LRaLt+0TWyXf3fZ9P76tpLsxEnwcwB8BmALwzwDUAfw/mqZ1bnVF+XweExMTnrgjXdexY8cOTExM4Nvf/ra4xvmy2SwmJiZQLBY9WF3X0d/fj6NHj6JQKHiwwEaqrImJCREQLNfx+7//+2LyWiXfyMgIJiYm8LnPfU4pH7Utl8sBcAee6rqOsbExHD161HUAN+fbs2cPJiYm0Nvb68ECGznr6L6M1XUdDz/8MGzbxq/8yq94sKFQCBMTE5iYmEAmk/Fgae/w0aNHkU6nPVjLsjAxMYGjR48imUx62k73JyYmkEwmPfpNJBI4evSoWNyQ9ZtMJgWe9MfreOc73wnbtjExMSFi3ngd6XRayMdxxEf3yfdk+TKZDI4ePYqHH35YaZve3l5MTExgfHxcaZuPfvSjmJiYwJ49e5S2GRgYgG3bePbZZ5W2yeVymJiYEAsnch25XA5Hjx7FyMiIUr6xsTFMTExgcHBQKd+zzz6LiYkJDAwMKH3/ueeeg23b6O/vV/p+sVjExMQEPvzhDyt9n+5ToD7HTkxMbL+FEU3TEgCuAvh5x3H+Uyve21kYuXz5sisnHIUg0AR2f38/Xn/9dddbhK+K8XyCPCbKNE0MDw9jbm7OtYRPfLThWxXHBWwkm7x27RrW19dd+QSJjzb+d3V1YW5uziOfYWzkmqM/Clchisfj6Ovrw+nTp8UkPYUT8HyCRBxLDwLlhKNQCV6HZVnij0IeeB1EtVoNjuO4sJqmoa2tDdlsFqdOnXJN6FNqL8uyMDAwgDfeeMMV9kAhMzRKp1ASrl9d15HP5zE/Py/Cobh+gRuhFqZpiq1zRB0dHSIWkOepk31nYGDAJT/VwRcEaAsel48yW9PxorJtKJ8g5c2TbdPd3Y0rV64I35RtEwqFEI1G0dnZiTfeeMNjG55PkM6elnNR9vf3+/oO34XUbDY98g0MDIgFSTmeFdhYfV5cXMTq6qornyDxkW+mUiksLCx4fJ8Wnng+QcJGIhF0dnZuu73DbdgYcS5txng7REGl9BuAyHPmOI5ICMo7e/5g0T05uLbRaIjkDDIWUOcT5Hx08Dh1ErJ8hKecanId0WhUtE3TNFc7iebn50XQtqwPcm5qp4y9ePGiuM7lo9+UjJOvvso6AOB6CRAf/UsPiiwfj+FUtZ3Ko/JVgc+E5XuHVRh6CHkdKysrno5L5Tt+8lEwL92T5aO2c9tz/VIQMAVLy/q9evWqp0zOR226dOmS0ja8DbfrO7Q6rMJSoD3PJ8h5FhYWoOu6y/a8PeSbFEvrp185kJ3/3kr6p1gY+UMAJQDfVd3UNO0jmqa9qGnai3zEEVBAAQX0T0F3tRPUNO3zAN4B4Bcdx/GmGAHgOM6XHMeZdBxnkuaXAgoooID+qeiudYKapv3PAN4P4J85jvP63agjSKoaJFUNkqoGSVXvlO7KwoimaX8I4L8F8LjjOMc34ycKkqqqr90rcgTterDbda+09Z7fMaJp2hcB/AtsjAKXNE3LXv/b2nPyrlM8HheT+5SA0XEcmKaJeDyOPXv2iGv8r1qtIh6PIxqNerCO44ilepqclesIh8OIx+NitVCug+cTVMlnWRbi8TgGBweV8jUaDcTjcbH9SnW/WCwiHA4r5Y9Go4jH467FD85XKBQQj8cRCoWU8lOcVjqdVsofj8fF9isZS1QsFmEYhgdLdisWi1hfX/fIt7q6KsrXdd0jn6ZpKBaL5G8e+XRdF3ieFIH+6KHiiWV5HTRaKRaLYmsh/6tUKqJ8le4rlQqKxaI4XlTlW/F43Ne3stmsyzdl+Wnle2RkpKVv0dZAuY56vY5isQjLspTyRSIRxONxsXAi1zEyMiJ0q5JvdHQUAIRvyXWQb3Z3dyvlI/+ixK+y32013Y3V4d+8/u/fStd/H8Cntroynk+O/0t/tO+Qh4uo+FRl3C4W8Oan88NyPr865N9EhmF4rvvJ1QrrVx+1Q8bejKw3I99m+m2lP469WQzRzdrmduW7mbbzvb2yfknnm9mG6uB0O/a/E6xf229Gv36+z/WymSxbQVveCTqOs/VStiC/YxMpjOD8+fPiGsdomuY6clJ15KbfcZ3AjTyB9GaV+S5duuTik+WjGCjik+vgxw7K7SQ6d+6cJ2edfGwitVPOCXf+/HmROp/CXOQjNy3LwtWrVz1YnvCA2s6x5Kjnzp3z5BPkR2620q/q+EzicxxHYHkYhgqjyvVISQf8jnUk/M3Ip/Idarsc28iP3KQRoco2dOQm+aZsGzqX+eLFi0rbUBu4Xm7Fd+g4WMMwlPJdvHjR5Tty2yksrdWRm5qmKXNp8uNwCS/rmkbyW0Xbeu9wPp93ZbmlCWfatWBZlpjE5lHpFKhLwcAylqLU+/r6RNS+XAftbAiHwx6sruvo6uoSE/cq+QifyWSU8u3atQuWZYlJfjmb7+7du9HX1+cqn/ORXENDQx5sIpFAT08PLMuCaZoeLOFt20Z7e7sHywOpaZFDrsO2bfT19cEwDKV8lmWhr68P3d3dnrZ3dHSI8vlCAvGYpom+vj6RqkvWL2Eozbys37a2NpHqSWUbwvf19YmMKJyvt7e3pe/09/ejr6/PFXTM206+SZ+ksn47OzuFb6lsQwtDsu6Ij3zLz3dGR0db+g7hR0ZGlPJ1d3e76pH1m81mW/o+6Y4WrWT97t271+WbMnZbLIzcLgULI+pr94ocQbse7HbdK20NTptTUHDaXHDaHC+f5xMMTpvbet8MTpsLKKCAArqPKOgEAwoooAeatn0nSBPhFE0OwDWZTmeMyFHphOVR+ZyPznHww1K0PE28y3x09qq8Y4T4KKqfYqr85KPPYt5O+iSQo/I5H5dLhZV3jMh89EnS39/vwSYS7h0jMpYywMg7RoiP7xhp1XZagJL1SztSAIgdHxwr7xiR66BzU0gW2TZ8x8hm8ql8h9rOd4xwHZFv+tlG3jEi89GnvLxjhPsM4Vv5Dt8xciu+I+8YkXW0Z88eABC+JeuIyh8YGGipX75jhPjo91bStu8EgRsOSVuh+MPJt91wPn4usIxNJBKe+3IdNC9kGIYHKz98Kvmoo0ilUkr5yOF4h8F56D5/4DkfzRvRKq9cB8lPD7xcB8279fX1KeUjogeaY3lAK99aRXyb3ecJNaPRqEd+fp9vDSM+fsDQ7t27PfLzw61oCxevg6da49siiY/P2alsw8+6VtmG+5bKNtRJc9/jfGQ73nbOR75FvifXQffJt+Q6uO+p5CP90yq1fJ/05+db5Jt02JksA/meyrfuBm3rTjCfzyOdTqNSqaBSqeDMmTMAbsTZpdNpvOc97xHXOF82m0U6nUZ7e7sHSxPbBw4cQKFQ8GCBjYc3nU6jra3Ng61UKti1a5cwnkq+eDyOdDqNhx56SCkftY0c7syZMy6esbExHDhwAJFIRFzjfKlUCul0Gr29vR4sAPz4j/840um0iMmS6+jo6IBt2+js7PRgz5w5g3Q6jXQ6LR4Yjq1Wq7BtGwcOHBDxZxxbq9WQTqdx4MABJJNJT9styxLlO47j0W+z2cSBAwdg27aIHeP6dRxH4OlsW14HdT7pdBq1Ws1jG13XhXxcfuKjssn3ZPkymQwOHDiAjo4OpW3C4TDS6TTGx8eVtunv70c6nUYqlVLahvYOP/bYY0rbxGIxcSKeyndyuRwOHDiAeDyulC+RSCCdTmNwcFAp32OPPYZ0Oo1IJKL0/Xe84x3iZaHy/fb2dqTTaezYsUPp+8ViEel0WsTKcmw6nQ5CZDgFp80Fp80Fp80Fp83dDLUKkdnWI8GAAgoooDuloBMMKKCAHmi6L4Kl+eoW4A7opclpVQQ6fQ7zVSc5GJi2UclY4qPyZSxNPvON+H58fvK1Kp/K3Kz8cDgsMqn4ld+q/ZZliU8/VR302ajC0n0/LHDjTJG7VT7Nqcp1ErXSHbDx6ajSHelEpTu+cNKqfGqbqnz6HN6Ovkk2uVPfpHJa2X8raFuPBPP5PJLJpDIrhWVZSCaTeOyxx8Q1zpfNZpFMJtHW1ubB0nzI/v37USgUlJkvotEoksmka7WL842NjYmVUJV8sVgMyWQSu3fvVspHbePHbsrl79+/H9Fo1JVdg/gSiQSSyaQ4TU7OyPHoo48imUwiEol4sDTnRvt/Zayu60gmk0gmk66VTs5n2zb279+PcDjswRJ+//79LvtR28l2yWRS7B/lOjJNE/v37xcrobJ+TdMUeNIfryOdTsO2bXHSnWybcDgs5JPbRYsmVL7KdzKZjKdtXEeRSATJZBLj4+NK2wwNDQnfUtmGVsz37t2rtA35Fq3Cy3Xkcjns378fsVhMKV88HkcymcTg4KBSvr179yKZTCIajSp9/5FHHhHhWyrfb2trQzKZFCEysn6LxSKSyaQ4kIljk8nkli+MbOtO8OTJk66cbLTFin43m02sra15cpY1m028+eabrmwbHEv/X1tbww9/+EMPtlwui/Ipk4ZcR6VSQalU2lS+9fV1pXw/+tGPBA/Vx3mOHz+OtbU113XOR9g33njDgy2Xy1hfXxc8MtZxNvIVlkolTz5C3nY/LACUSqVN5VtbW8PFixc9bV9aWnKV72ebm9EvbTvjddRqtZvCrq2t4dq1ax756HRAP/nOnj2LtbU14Vt+tnnllVeUtqHsM376bTabKJVKnlyMsm38fOe11167Kd85derUpr6j8v3NbEPPTLVaVfr+yy+/7KvfZrOJkydPbmk/sq07wYACCiigO6VtPydIbyX6DdzIJ+g4Dt58801xjWN0XRdvUhqRcb5arYbZ2Vkllu7LR26qcqrxnHK8jmq1inq97jo6ktcRCoU8beM8ADA7O+vJCUd8lUpFfErIWF3XcebMGRH+QCElnI/i0S5fvuzB8pAZCm/gWNLH7OysK6ccr2N9fb2lfunfRqPhsY3jOALL8+apMKpcj1evXoWmaUIWx3E8vtNKPsLSPVk+ajvPxcfbTjajc35l/S4sLLh8U7YN5ROUz6vmfPw4TT/fkfMdEh/l+zNNUynf3NycuK7K10h5GLlvch2Rb6pyadKZ2dw3ZV1vddD0tu8EubH5GbP0cFy5ckVc45hwOOzil8+npc8yFRaAMBDvBDlfuVwWh2v7ydds3ojzkusIhUICSxk95JjOpaUlz7m+xEfOSe2UsVevXnWlT/fTJcWDyTonooeUY+lf+qyV5aOH1E+/VB6Vrzo7mGNlmTmGd5xE1MHzDprXIcv3mONgBcCrkm2I308+We/0f9KJn234oeUq25BtVWdCy3Xcju/U63XRCaqwy8vLrnrktl+9elXEZ/r5PrVTxtbrdViWJWwj65f73lbRtv4cDk6bC06bu9unzf1cJIK/BfBHAIzrq7nBaXPBaXN3jYKkqupr94ocD1q7fsWy8B9XV/FDAD+j67joExKy3dq13W0YJFVVUJBUNUiqysunT3Tg9pOqfhjAH6+u4rsAfhbAVabjIKlqkFT1niL+QMi/NU1zOazMx2OfOJboZrC8DE48mLmVfBS4ervyyddvFkvxdSos/61yOlXbVVg/+e62bTied6hEcmybjP1XAJ5xHPxdKISfNgxcvU35bka/fr4j+5Xqt5/vtLLNzch3s77j5/utfJqX7+dbrWRQ1XentO07QT7vQm8QmteJx+MiGJnPOxAfne0qY2neZnx83BdLZ7vSGaoy39DQEAC45l04H53tmsvlWspHIxbeTnKe8fFxV842zsflUmHp3GGeT5Dz0du+p6fHg6Uzd+ns2Vby8Tk7zhePx1vql8qnM3o5n2maAsvzCRIfYeLxuEjkwOvo6upS+47j4HOahn/rOPjTUAj/9u1vx9r1h9VPPpXvUNv5fDCvj3zHT2/9/f0u35T5aKSYz+eVtiHfoq8JP9vIuShv1nfy+bzwHZXvj42NiWdE5ftUfl9fX0v703wyx96Nc4e39Zzg448/DmAjMBfYSBRJvzOZDPr6+vDNb34Tvb29Lpxt25ifn0c2m0WtVsP09LQLCwC9vb344Q9/iLe//e04ceKEC1sqldDf34+uri6srq4iFou5sADw/PPP49Of/jTK5bII7uR15HI5pFIpfOYzn8HP/MzPeOQrl8tIJBJYWlrC7OysR76xsTF873vfw0MPPeSaHyG+QqGASCSCixcvoru72yPfyZMn8cQTT2B+fl7Mq/I66OD0973vfXjqqac88hGdOXMGi4uLLqyu69i7dy+OHDmCXC4nVlkJOz09jfHxcRw5cgQDAwMi3RXdP378uHh5nTx5Evl83iV/LBbD3NwcHn/8cczMzIhVRpKBT5739/fjG9/4hkv+J598Et/97ncBANPT06jVaph4+GF86Lvfxa/Vajjc1oZnH3kE//nv/x6dnZ2uEA7btnHmzBns3LkTwIbvybbJZDJ47bXX8GM/9mN45ZVXPLbJZrPC90KhkMc2zzzzDP74j/8Y6+vrwvd4HeR7zzzzDPbt2+exzdLSEjo6OoTvyfLlcjkcO3YMjzzyCE6fPu2RjxYcqRxZvn/8x3/EwYMHcenSJXR1dXnuv/LKK/jVX/1VXLx4EXNzcx756eX4oQ99CB//+Mc98tNn77lz57CwsODCku8dOXIEt0Kt5gS3dSfY398vYo0AoLu7G8vLy644tvvhGm2zunjxouC5H67dK/oNOw6+Zhj42bU1fNYw8PuGAWhaYK971IZdXV1bujDi2rLyVv89+uijzq1QX1+fYxiGA8AB4BiGIf6vaZpjGIaTTCZdPMQXjUZd/Pw3AEfXdSeVSjltbW0erGEYjq7rrn/lOuLxuJNKpRxd15XyES4ejyvlSyQSLctPJBIty6ff1E65fNIL4eU6DMNwUqmUE4lElO2nP9KzXEcqlXJSqZSjaZqvfKlUyrEsy1N+KBRqWb6maS3L5xgqi9dhWZaTSqU29A843wIcB3D+5fWyCJ9KpRzTND3yRSKRlr4TiURE+a1s39bWprRNLBa7p32T5PYrv62tbUt8k2wrY/v6+m6pn3A2Rnov+vU798XqcEAB3Q51OA7+E4C3AfjvNQ1f0XVAscId0P1NQScY0ANJ2WYTzzebGAXwiwC+oVhFDujBoKATDOiBo2HHwV+Vy+gE8NMA/jOArY8+C2i70LbvBCk/HQBxsA2tFgIQOc34NeKjfYi6riuxdNYDj0lT8amu0SE2PLTAj880TWV51K5YLAbLslyrrBQqsbq6KspX8YXDYbS1tbmuEd/a2pooS4UFgGQy6Qpmlfkof50Ka1mWK2ZM5qNQD776KutoM9vw8m/GNkXHwV/V6zAAvMcw8ANdR8gHS/KpyuOByzLWsixhm1Y2pNyEKttQ0tZW+o1Go64FBBWf6pppmjflO3RYmao8qlelt2g0elO+39bW5muvzZ7rraRt/Q2Qz+cxNTWFWCyGWCwmwipisRgGBwcxNTUlQgCIh/iy2Sympqbwtre9zYONxWLI5/O4cOECCoWCBwtshKhMTU3Btm0PNhaL4atf/apIfKmSb3x8HFNTU/jzP/9zpXzUNsLs3r3bxTM2NoYLFy5gbGxMXON8ExMTmJqaQm9vrwcLbISGTE1N7HkXFAAAIABJREFUYXh42IONxWLYt28fbNvGJz/5SQ929+7dmJqawtTUlAgV4Vjad3vhwgX09fV5sKlUClNTU7hw4YJIfMrbblmWKL+7u9uj366uLly4cAG2baOrq8uj3+7uboF/97vfLbA/Hg7j7+p1WPE4fuvRR2G94x1IpVIe2/T19Qn5qBPidaTTaVG+yncymQwuXLiAffv2KW0zPDyMqakpjI+PK23zhS98AVNTU5iYmFDahhL2fuc731HahuxDZ27LdeRyOVy4cEHEMsp17N27F1NTUxgcHFTK953vfAdTU1MYGxtT+v6LL74I27aRz+eVvv+2t70NU1NTePrpp5W+XywWRf0ydmpqKtg7zKm/vx/lctm1j5F+h0IhWJaFYrEoYsKIEokEqtUqwuGw2LrEscDGCOqxxx7DSy+95HoDEZ9lWQiFQiKmSd5ONj4+jtnZWTQaDTHi4nVEIhGYpolcLueKJSM+Kpe2LsnyRaNRPProo/j+97/vOq2N+GKxGHRdF+2U5du3bx+mp6fF1iVZPsqInclkXLFkxEdUqVRQq9VcWO16dmLbtvHCCy940uevrKy0vE+xlwCwtraGaDTqkl/XdUxNTaFUKoltcVx+CsgFNkZli4uL+OeOgz9tNnEWwAd7enD8uk3pxD0uv2maIknB0aNHXVlOEomEa1scxXNy+UKhEPbt24djx45hZWXFYxtKK0/by2TbjI2NiTRZ5Hu8DvK9sbExvPTSSx7b0LY48j1ZvkgkgsnJSbz44ouukSTxUZA2lSPLNzExgZmZGRHnKN+fnJzEiRMnfH2LfDObzXoSpCYSCbHtT+VblEwh2DvMiDLk0m8Arpxtr7/+urjGMYZhiCMtKYsw56tWq8JAMpbu81RXMh8FifKccrwOMjAZU64jHA572ibnhDt58qQnJxzxUd47aqecE+706dMi+zE9jJxvfX0d0WgUCwsLHiy1A7ixZ9RPPjmfINW5trbWUr/0b71e99jGcRyB5emgVJhqtYr3NZv4suPgVQA/BWD12jWhFy4/1UF4P/m4DlS+Q23nufi4fmm/LE1JyPqdn593+aZsm2q1ilAoJLKjy7Yh/+d6UdlGzndIfJTvLxQKKeV78803xXVVvkbKw8h9k+uIfJNiDWX9RiIRl2/KuuYv4a2gbd8J8oeAHJqcB3BvKucYwzBc/Pw3/euHBW5s7uadIOdbW1sTb2M/+RzHEQZW1cGTtsq536htcr4+flwAANFOGUtJJagNfrqsVCoerGqz/K3KxxMDyPepPLon24bKpmsq/dLv36jX8e8dB/8FwH8NYBmAVat59HIrvnOz8sk64fqWbcyJzzmqbEPyqXI98jpu13coiQR96chYSvrAdc95VlZWRAKFVvqlUaif75N8Kl1vJW37ThAITptrVf4De9qc4+ATjoNPOw7+2jTxRLOJdU2DgeC0ueC0OTdt+4URv6ExJdakyWGZstksEomE74Zsy7JE4kq/+4lEQsxdyTQ4OCgSV6qIElIODw8r71Pb+APLiSbH+QPHiRYo5H3TROPj40gkEq6UYpxozq67u1t5n+Zm/LK5UNJTv4witHiSTCY99+VErTIZhiHKV97XdRwyTXzacfC8ZeG/i8WwzrKP0El6iURCmZWEJ31VZaFJp9NCPhVlMhmRVFVFlLSVkkDIRElb/XyLfJMSFchEvuXne7lcTiQ9VRElZR0cHFTep0TAfr5HCX/97lNS1Ww2q7xPyWZVvsX3hW8Z+W0leSv+gm1zwba5O902F9Z159nr2+D+vaY5YdNsuW2ule8E2+YejG1z23okCLjPVOXDcIozorchz39GfMQjY+lPdWatqg4Vlt6Cm8lHb+NW8snl0OiFnxtL13h7WmHpzFgVlsvAzw3mfHLbb0e+VvptZRsqm4iuRzUNzzWb+BcAPq3r+JeGAZ19+tEfjU79bLNV8t2MflV6C4fDrvL9+Pi5vrdjGz/5NsNGIpGWz00rn+bl01fIrehXNTK/U9r2nSAZkzutrusIhUKIxWIYHR0V12S+WCwmDMqxZCD63FBhLcsSeBmr6zr6+/sB3HAYuQ7CU5ydn3zUmfJ2kiOMjY3BsizXNeLjcqmwIyMjiMVi4oGT+ejl0d3d7cHyuC76ZOFY+kQdGxtDKBTyYCk4u5V+qXzTND36NU1TYOlw9qSm4RuOg18A8NuGgX/X1obY9UPE5To6Oztdsvj5zs3Ip/Idaju9aGT9hsPhlrbp6elx+abMRz4xNDSktA33zbvhO0NDQ656ZB3RFA/vzDkflU+fw376Jd/hWL8pgjuhu9YJapr2m5qmndY0bV3TtB9omvZjd6uugB5s6nQcfLvZxH/VaOBDmob/oJgnDCggP7ornaCmab8M4A8BfAbAIwD+XwB/pWnawFbWQxHp/C2SSCRcb9vR0VFxjfMNDAx43uacz7IsjI2NYc+ePR4sTQrHYjFEo1EPlkaCNPmsko/ehnTillzHrl27xNuWJsk5z+7duz1vc85H5dNbm2MTiYRnJCjXEY1GYds2duzY4cHKI0GVfDRxz0eCxMdHgvJIMxaLoaOjwzMS5DqikaBt2xg0DPx9s4kigN/o6MD/nUwKTCwWQzKZ9Oi3s7PTdRKgbBs+Emxra/PI19vb29J3+vv7PSNBriPyzWKxqLSNPBKU9UsLI/JIkPj4SFBlm9HRUYyNjbl8U+U75CNyHfJIUNbv8PCwWBhR+T7pjhYnZf3SLivyTRm7LXaMaJr2jwBedhznw+zaawCedxznd/1wwWlz6mv3ihz3WruGGw382fVECL+g6/hBW9t90a7tdO2tqPOeP21O07QwgEcB/Dvp1t8AmNrq+oDgtLkH8bS5YqOB/2tlBRqAHwfwUrMJ4zrfVpw2R3iV7wSnzQWnzW1GXdjITHRBun4BgCcwSNO0j2ia9qKmaS8uLCzcBXECut/oHY6Db5TLqAB4J4CXNgMEFFALupurw/J3tqa4BsdxvuQ4zqTjOJOZTOYuihPQ/UA/A+BbAC7oOn46kcCPFMHOAQV0K3Q3ts1dAtCAd9TXDe/o8I6IJkhpBJnJZFAul0XmkUQigfe+97346le/6tr0n8lkxHagZrOJxcVFFxbYiGr/tV/7Nbzwwgs4duyYC1su38hUwjNt8Dre85734PLlyzh+/DiuXLnikS+ZTIqUUX/zN3/jkS+bzWJ+ft6VSYOPlCld0uHDh8XnArV5YWEBnZ2dMAwDyWQSy8vLLmw0GsUv//Iv4y//8i+xurqKlZUVFxbYOMxmdHQUuq7j6NGjLizfKXHt2jWxA4SwmnZjx8izzz4rtoARdnFxEel0Gk888QT+7M/+DBcu3HCLTCYjspMAwJUrV9De3o5yuYyfK5fxvzYaeFXT8Pyv/zq6ZmYQ/f73xbY70m+9Xkd7ezuAjfkkSkhAVCwWEY/HcfLkSSwuLqLZbLpsY1kWkskknnjiCXzlK19xZYLJZDIilyH5nuw7/f39+Nmf/Vn86Z/+qTgIjOs3Ho8jFttI9zY7O+uxzdTUFF5++WU0Gg1cvnzZY5tEIoFisYjBwUE899xzHtvQ53itVkOlUvH4zujoKN71rnfhq1/9qjjpj9fR3t6OUCiEdDqNxcVFj3zvfe978bd/+7euDD9cv7/0S7+EU6dO4eWXXxa64zpKp9PQdR179uzB9773PY/v53I5nD59WvBzbCaT2fKFkS0fCTqOUwXwAwDvlm69GxurxFtGJ0+eFPNn0WhUKIfm08rlMqanp8U1zjc/P49yuYzV1VUPNhqNotlsolQq4cSJEx4scGOuZn193YONRqO4fPkySqUSKpWKUr5arYZyuSwyach1UNtoj2o+n3fxzMzMoFQqodFoiGucb21tDeVyGXNzcx4sAPzoRz8SLwEZG41GsbKyglKpJO7LdZTLGynMKICWYyORCEqlEkqlklhp5ljLslAul1EqlbC8vOxpe6VSEeWbpol8Po/fqFTwvzca+EfDwM+1teHo9faHQiGPfumlVC6XRSwaryMej6NUKon5N9k2uq4L+eT25/N5LC4uivJVvrOwsIBSqSReLrJ+aZ5zenpaaZvV1VXXy1zWb6PRQKlUwuzsrNI25FvUUct1nD59GqVSCbVaTSkf6X92dlYp3+zsLMrljQQRKt8/deqUR3ecj9pHe/dl/b7yyiuuFzvHlstlT/qtO6W7lUDh8wAOa5r2PQBHAfwPAHoB/C93ozJSGDk2cCOfIJ9YJiqVSp58ghwLQCyc8JxuvA6eT1DGAhDZeHk+Qc5H+QSvXr2qlE/OJyjXQQ5Jb3tZPgptqFarSvmoc+U53zgf5RM8d+6cUj4iGqlyrKZpYnJ7fX3dlSCBOodW9135BFdX8e4XXsD/VKng6wA+4Diorq1h5Do/zydIMtAoHQCOHz/ukZ/SqwE38gly+SmfIL/P5ZPzCcr6pVEsLcrJtuH5BFW24SmmyPc4H/kedSaybeR8gnId1DbyLVk+nk9QJR/JJNudiMr08y3yTRqIyPLTwo7Kt/z2a98J3bWkqpqm/SaAfwWgB8CrAH7LcZzvtMLcTojMxYsXxUNkmiaazSaazSY0bSODRXt7Oy5fvuzKWWaaJkzTFPkAG42GCwvcSH1OIyqO5TxkMI4FgFQq5VmZ5nXwuKvl5WWPfPTJxevinQV9UtHnHPEQn2EYYotYvV53YXV9IyX81atXPfJwXcZiMTQaDdeDzHcOAO7VVF4Hfe5funRJdCK8DtM0hQz0oFC9RJrj4Ol6HR8H8BVNw687DhrXR56ZTAaVSsV1cDvp13EcsYpI9XH9RqNR0Ym08p2Ojg5cvnzZtXppmqYrDRu1hdufPqeXlpZcRziQLPQ7Eol4XgK6riMej4t8grRay/VLfmMYhiddv67rrtVZx3E8tonFYojH45v6TjgcFnkzeR2pVArXrl0T9ci+39nZiXq9Lr40ZP1S+ZFIxPUSIz76kqF0XRxrmia6u7vv7RAZIsdx/gjAH92t8ol4CAZPukn/p1ES7+yJj/PLWMdxWmIJQ52gzFer1VwPpYwlPHUAfvLxTlZ+YVGuPxnLN4fLOiKihJd+fCQDvShknRPdrnytbKPrOrRGA19yHHwQwBc0Db+taag7DsBGCVwGWb9yElFeR6PRQCgUum3fMQyjpe9Q22WdcH3LNubEO22VbehfnjhVpd87sQ35tApLCYv92l6tVl2dsKxf+qPOtZXvy3XIdt0Kui/yCQZ0f5HlOPg/mk28F8AnAfwbTQOCVeCA7hLdF51gq8SNtHKpikAPkqree0lVe+JxfHVlBe8E8FFNwx/dYfk0t6iaEwOCpKp3wzfJJnfqm1ROK/tvBW3rLDJBUtX7K6lql+Pg6ysrmGo08OFoFF9JJm8vqer1jiORUCf+DJKqBklVOW370+YuXbokRnThcFjMp9BG+0wmg/Pnz7vmEmhjNs071Go1FxbYeJB6enpw5coVz9urXq+LtxzVJU++d3V1oVqtuibfeR2UAiqVSrkmqIkvFAqJuRfi5wsI8XgcHR0dOH/+vGs7FfGFQiGxYMHnHokvm83i0qVLaDQaYvKd10GhJ5qm4erVqy4s79j4wgKvI5lMIplMilPTZPnC4TC6u7tx8eJFVKtV9DsOvgVgCMD7DQN/fb0Oihnk+tU0DX19fVheXsa1a9fEnBLp13Ec8QAZhuFaYAIgNuOvr6+39J3u7m6XfomPtqQBEG3h8lmWJY7d5OcHU9uN6+fx0sKDbBuKiyTflG1jGAbi8Tii0agrxpL45AUL2TaxWAydnZ2Yn593LbYQH/kmLarJ8mUyGSwtLbkiGLh+s9msiD8l3XEdkW/G43GxOMf1S4uCxM+x4XAYXV1dW7owsq1HgoA7Xx+92SgTSCQSwdDQkLgm80UiEZHpgmMJT6M0FTYcDgu8jNV1HT09PQDgyqTB+QhPIzU/+WikxttJo5Ph4WFXzjbOx+VSYQcHBxGJRESWFz9ddnV1ebCRSET8UWfDsTQ6Gx4eFg8UxxqGgUgkIvRb0DT8F2zEUP1iPI5vXq+fwohk/RqGIbCUBYbLTBgaEcn6pUBqkkW2DeFb2Z/KV/kOtZ3rhOuIfNPPNt3d3S7flPnIJ3bu3Km0DffNu+E7O3fudNUj62hgYCNZFM8gxPmofNoh5qdfwnOs3+j1Tui+mBMMaPuS3Wzi69dHAj+h6zhpmgAbeQQU0N2mbT8SXFlZEam5aR6QMmCsrKzgxIkT4hrnC4fDWFlZwdramgeraRpqtRqmp6c3wjUkbCKxcXj7ysoK1tfXPVhN08S5w3TGqlxHpVLBysoKzp49q5TPMAysrKyIbU/0CUF/hmGIw9PpGudbW1vDysqK+LTm2EQigZmZGaysrKBWq3mwmqaJgNgrV654sKFQCCsrK1hZWUG9Xvdgiaanp0UcIceS3VKlEv6mVsMagMcNA/+ffiOXH5VPn09cR47jiJ1AFAbC9dtsNl14Wb8UvE7xj7JtKDZyenpaBI3zOizLEuWrfMeyLExPT2N1dVVpm1qtJo6lVNlmYWHB5Zuyfunz9PTp00rbkG+tr68rfcc0TUxPT4vzheU61tfXsbKyAtM0lfKdPn0aKysrYtpF1u+pU6cAQPiWrCPyzcXFRaXva5omfEvG8pjVraJt3Qnm83nE43ERT0TzLxRgGY/HsXv3bnGN82WzWRFwLGMdx0E4HMb4+DgKhYIHC2x8DtO8jIx1HAe9vb1i8lkln2VZiMfjIr2+XAe1jT5pKCaM/sbGxjA+Po5wOOyKxyI+2h7W29vrwQIbk9vxeByhUMiDdRwHsVgMtm2jo6PDg63VaojH44jH42J+kGMBwLZtjI+Pi3RfHAsAvxgO4/+8dg0XTBPv1HXMXNdBrVYTuonH49B13aNfXdcxPj4uFi5k/eq67sLL+qWEsTydGK+D5tzGx8dd5RNfOp0W5at8J5PJYHx8XPiWrN9QKCTKV9kmk8kI31LZJhwOw7Zt5HI5pW1If/TpKNeRy+UwPj4Oy7KU8kUiEcTjcQwODirly+VyiMfjCIfDSt/P5/Owbdvlm5yPgrXT6bTS9///9t4/OK7rOhP8br/+3Q00GyDAJgACBNGgmoSaakOayEnFDmeTbO24sltOrdZZ2fFamyqb46kpx1lnJhntrq1UVHEp2fUynngdaVKeOHSVPbXemnLsYpxks6VKllHW0chYM7QY6AcESwQhgiQECALRjQbe/tH4Ls+773YDBECDLdyvqovg6/Pde+65591+795zz+Xebqb7ktxMJuNOm5Nwp82152lzHwH8VcD/T57n98XjofJ3ctqcyWFZsg532pw7bU5+2vpJ0KH98C/X1/E1AH8N4L/KZHDDBUE77DHcIOjwY8Nn19dx1vfxHwH8l5EIltwA6HAPoO0HwWg0GpjwBW5HqkejUR0OYU4Mk+uJRANSLhJpJBloxmX5/JhyDM3YTL+OjXMxWulnlsPFh3w+r/XnNcqx/GbcXC6nw1dMrtQhnQ6fHSzbbdOP83D5fL4R+gDgiwA+6/v4k0gEv6QU1mKxlvZl+TIuUcq14pITjUatZ9tynlW2X9ZB/lb0s/kO2y5tIm0k22TrGx4wxb415Wjzzs5Oa99I39yO72zG7dw4zEr6mJTJ5XK6H1r5PoPB7+TetAXf7xRtHSJTLBZRLBZ1aqdKpaL/7unpQX9/P86fPx/aNVGpVDA7O4tCoaBXgSUXaETt//CHP8R73/teXL58OcCdmJjAwMAADh48qNM+SS4AfO1rX8OTTz4ZyH8m6xgeHkYul8Pv/M7v4AMf+EBIv6WlJWSzWczPz2N6ejqk3/Hjx/G9730P999/fyBwlHKlUgnJZBLXrl1Db29vSL8XXngBjzzyCGZnZ/VhVbKOcrkMz/PwoQ99CI8//nhIP+L111/HjRs3AtxIJIJTp07h2WefxejRo3hyZga/tLqKL8Xj+OZ734vcD3+IsbExPPvssxgcHAxkgqlUKnjxxRf1gtbLL78c6GOgMUhMTU3h9OnTmJyc1CvZ1EHuKhgYGMB3vvOdgP6f/OQn8dxzzwForGCvrq4G9M/n8xgaGsKzzz6Lrq6uQDBvpVLB66+/rhe0JiYmQn3T09ODl156Ce973/tw8eLFUN8UCgXte7FYLNQ3f/AHf4Avf/nLWFlZ0b4n66DvPf3003j44YdDfTM/P498Pq99z9RveHgY3//+9/Ge97wHU1NTIf24W4nlmPr91V/9Fc6cOYPr16/j4MGDoe//9m//Fh/96Edx7do1HSUhdRgbG0MsFsNjjz2GT3/60yH9GYR95coVzM3NBbjS93YLbT0IMvGoDXNzc5ibmwsNMIS8+W2YmZnB6dOnAwOgxBtvvKEHH1vHPPnkkyHnkKDzmQMMsVniyMnJSZw+fbpp5LzU27b17ZFHHmmpn7x5bWjFZZ68//x978MXr1zBP6vX8duJBP7XRAKVjRCRiYkJnD59OjAAEtVqtWX5y8vLOH36dFMZ5vhrhueee67l9/Pz85ifn8fp06etWUtu3LiBGzduNOXPzc3h9OnTTW1o/vCY+PKXv9xSP/remTNnrN9PT09jenq6KZ8/IHIAlJC+x6dhiTNnzmj9Dh48GPr+ox/9aEv9Gd70x3/8x9bvW/kef+R2FXu9Iiw/21kdlqt30Wg0sMIVjUb9fD4fWoGKRqN+KpUKyMu/yc/n86EVOMqxfM/zQlyIFTJ5XcqRZ1shjEajfjab1fKUlTLZbNbP5/OhFTjK8W+2U3I9z/MPHDjgR6PRwOqwlOPqaDqdDnGj0aj+cIVQcpVS/pHOTv//iUb9NcD/5wZXKaX7JpFIhNoei8V0+ZQ1V4fz+XxodZhy5ESjUT8ej4fsy9Vblm/2jdTPtF00GvWTyaQu3+Y7yWQy5HfSRvSdjo4Oa99kMhmrr0jf5Oq4rW/kvzbfyWQyLfXj39TDrIO2k/eAlKFu0jdt7bGtPkejUb+zszPgm+Z9uturw239JAhAxy7Z/vZ933oMJOWacQnbU4CtDhtkYstWuraqg+Urpbatn3lNcs3r263D1O+g7+PbS0s4ub6OXwbwdYPLT7OyzfJ2u+2mfW1yrfRrVtdu6reTvpFlt/KdZnVspW2b+RbnhZuVsxX7blWfnaLtF0Yc7i0cAfA3AEbX1/HhdBrfcCvADvc42noQ5ASut5GVg+mNPM/Tqa4eeOCBgAzlBgcHdSotk+ttbPCvVCo4efJkiMsUV+SbXM/zMDIygkqlgnQ6bdWP6YqYlsis47777tPpkGzlM10RkwCY5TOV1tGjR63l33///Todkq39TKV15MgRa/tlqirJ/U2lcAjAb1QqmPsn/wTxeDzEjYpUVb29vaHy8/l8IFWXqVssFtOptLjFUJYvOSxL1tHb26tTaXEVUtbBVFeVSgUdHR0h/ZjqivKmfgMDA3pHSivfZMoos/yhoSGdSquVb5ZKJWvf0GfoA2Ydo6OjOs1bK98ZGRmx6lcqlQL3gCnD3TzNfJP3DO9Bs/xTp0419U256LVbaPtUWnJxgxOmtiSN7XztXtFjK9fivo9hAFc2Qn/eLe26k2v3ih7v1rYWCgWXSkuCZ0XwTAge7sKN1wwVoAzluPk7upF/TXLJ7+7uDvCkHDe/exv56sw6GAPlbeQcNOvgEw0zF5t18BwMz/OwtLSkc+Lxs7a2ps9v5TUpF93Y/F6r1ULcpaUlfbasUirEJR8IHpAj5WKxmOZL7orvg+vS3Btqcn3f1+farqyshNq+vLysywcQsq/v+4F9p6Z96Rc8DdC0r/Qd8mUd/K67u1snYZB1yPbbfKdaraK7uztgO2kj+iZlTfvyNDnyzb7hfFtHR4e1b+hbzXynXq+ju7s74JtSjvx6vW7Vr6OjQ6e5svk+4wTpW6aN6JvMKm3al/5BvuTaEq3uFG0/CMp8fQyCZTBqIpHAwMCAvmbKJRIJnVNNcjnAMRbMxmU2ESaBNOWYK03mbJNy5DN8pZl+7HTZTpnXTeZsk3JSLxu3v79fl29ypQ464Nmogx8Ols30Y0CyWUcikWhpX5bP10cp53me5noinyDlyEkkEvrHSNbBbNab+c5W9LP5Dtsuy5f1se+b2a27uzvgm6YcfeLw4cPWvpG+eTd85/Dhw4F6TBsxl6YsX8qx/K6urpb2pe9IbrNs1TtB2w+CDg4ODjvBu2JOkGEw3sZ2InkWbiaTCaT5phxTm1NecoHGrw5fN95+++0AFwiGH3A7kayDqYB4PqupH7cU8exZUz+evdqsfE4UNyuff7OdZvm0C18HzfZzErparQYO6KYcIV/PZB18JVpcXNRhDaZ+2WwWKysrgeMz+WTHNtnKV0rppzlb+ZLDsqRuzGy8tLTU0ney2azOayf149EHpq1ZFrM6NyuffZ9Op7G8vBzqG6aXb1b+Xvsm9W5WfkdHByKRyI59k1MfJne35wTbPk4QcKfNtSq/3U6b2+3y3Wlz7rS5zdDWr8PutLl312lzhBl2Y8Lz3Glz7rS53YN7Hca9+crhXofd63Cz8vfaN99tr8Nt/STo4ODgsFO0/SAoD4Lh4zPnImKxmA5VkTnLKMdYLJNLfisuy5fxbLacas30Y6yULWedqZ9ZDn99e3p64HnBnHCUY/nNuN3d3TqWzORKHTo6OkJc2W7+Sksun8Ka6ccwj836Rsqacq245MgQDinH10SWv13fYZtN/dh28/Ap/p/lN+sbxuHRN0052ryrq8vaN+z7VnX09PSE8h3K9rTidnV16b639Q1jOGXf23zflktT1ku+5N6NOMG2XxhhvBYQzFlG8IB2CVOO25BMbn9/P+bm5gLzVjY527XR0VFMTEzo82mbyY2MjAReaaScfPVspd9mutiuDQ4O4ubNm5vKnTx5smV5rWy/mX78Xs5NmeVw/qtZ2byR7qTtDzzwACYmJjb1na3YdzP95LyYTe7hhx8OXSuVSvpaK985evQo3npwtk2bAAAgAElEQVTrrU3b20o/OW9ok3vooYda1mvjDA0NYX5+PjCnaZM7ceJEU10365vdRFs/CRaLRZTL5VDGGO4oKJfL+OY3v6mvSblCoYByuaxPk5Nc3/dRKBRw/vz50GlzlOvr60O5XMbIyEiI6/s+fv3Xfx2VSgXHjh2z6jc0NIRyuYzPfvazVv3YNgbsrq2tBWSOHz+uE8bKjBuUGx0dRblcRl9fX4gLAH/yJ3+CcrmM3t7eENf3fZw8eRKVSgUf/OAHQ9y1tTWUy2WUy2X9JCu5XBg5f/48crlciBuNRlEul3H+/Hl0dnaG2p5IJHT56XQ6ZN9UKoXz58/rhSfTvul0WvM5+S7reOihh1CpVFAul/WJZrKOXC6n9YtEIiH96Fv0PVO/np4enD9/HidPnrT2TW9vL8rlMsbGxqx98yu/8isol8sYHR219g1PMjx79qy1b44cOYJyuawX3cw6hoeHcf78eQwNDVn1O3bsGMrlMoaGhqz6nT17NuBbpn2feeYZVCoVFAoFq++XSiWUy2U8+uijVt+nbbu7u0PccrnsFkYk3MKIWxhxCyNuYWQrcAsjDg4ODk3gBkEHB4d9DTcIOjg47Gu09SDokqree0lVPa9xnCJ3dLikqi6pqln+vZZUta0HwZdffhm1Wk3vQeRKYDweh1IKtVoN165d09ek3OzsLGq1GlZXV0PceDwO3/cxOzuLy5cvh7hAY/K5VquhXq+HuPF4XJ+Ytra2ZtVvfX0dtVpNHxdp1sG2cWK5UCgEZCYnJzE7Owvf9/U1KVev11Gr1TAzMxPiAo0T07gwZHLj8bg+EY6pj8w6arUaarWanhiXXB4jOTs7G+BTjnkOZ2dnsbi4GGo78z3WajU9ES5t5HkeZmdndRiRaV/P8zSfq9eyDh6RypySZt9EIhGt3/r6ekg/2o6+Z+o3NzeH2dlZ7VumfYHGvvVLly5Z+4a+Ua/XrX3j+z4mJiZCtqMc+VxwMOuYmpoKtU3Kra2toVarYXp62qrf4uIiarWajqIw7Ts/P4+JiYmAb0q51dVV7Xs237948aL2LZNbq9U2PYnxTtH2cYJyA/bs7GwgqWYikUBPT0/IaBwAORgtLy8HuCy3UChgZmYm0GGUYz61tbW1EBe4HT/neR5u3boV0i+ZTCIajSKVSoW4XPHmYLa0tKR1Jjg4vP7663olkTdurVZDOp0O3MySW6vV0N3dra/zI+W475U3hOTOzs7qG4IJLyVXOu9rr72mV1fJXV1dRSaTCX3PtvMsZwC4detWyL6RSET/kNXrdf1DQTk+aQON4zNN+8ofHq6Qyzqi0SiSySQKhQJeeeUVyAiK2dlZVKtV3b9su6yDg8PVq1fxzjvvhPqGN/z6+rq1b/jDQd80+4Z1d3Z2WvuGORHX1tasvkP7vfHGG3rlX9bBJ8R6vW7Vr7OzUw9mNt8/cOCAzvhOrpSjb9J+EvKHvVqthupotl97J2j7QZDL6PwbuB13tL6+rkMZpCPL8AvyJZf/tuKSY+MC0De2jK+y6SflbPpJrhnOtLq6GihfKRXSycZVqnH2L68z1MGsz/M8HScmubLsnei3Wd9IO0s5s29s9rWVQ/AJSbZ/O77TTD+2XdpEtt3Uz7QvB3Z+Z/YN/6V9zb6RvrmdvtlMP9m2Zm3fSt808/1muph+vWuQgYo7/QDoAvBvAVwGcAvA6wC+DKB7K/ztnDvsGWeneuLsWM/z/M7OztDZpp7n+alUKiAv/wZun+1qnjtMucjGebv816wjk8mEzl6VcuTZzl71PM/PZrMty89msy3L599sp1k+7RIR5w6btszlcn4ymbS2nx/a2ayDZ88qcS6wqV8ulwudO+x5nh+LxVqWr5RqWb7ksCxZRyKRCJ0JbfMdnq9r6pdMJlv6Ds81blY++9R25rTneX46nb6nfZN6NyufZ27v1DfZtyZ3t88d3u05wT4A/QD+NYAygF8G8H4Ej57dNbiFEbcw4hZG3MLITnHXd4wopT4A4DsADvi+v9hK1p02Z792r+jh2rW/23WvtLUdM0t3AqgCWL5bFXCOh0aSW5983w9tm1taWgokE+UEsiyLE7fmlivK+WJrkskFGvMd3JfKeQwpx61J9Xo9tDVpaWlJT+w3K59biJqV73mNrUZc9DDL58IH+WYd5PPkO7P9hI1rfs8fWpt+trbLRKa28rki3ax8uW1ueXk51IcrKys6MUEr3wFg9R2ZycTmO3LhxFY++54LH2b5MunCveibXLRqVj7tv1PfZN/auLuJuxoio5Q6AOC3Afw73/frTWQ+oZR6Xin1vJmxw8HBweFuY0uDoFLqSaWUv8nntMHJAPg2gCtozBFa4fv+M77vP+T7/kPM33YnYEwgwzI26kY0GkU8HtfH/8mcZZRjTJvJZc42pqa3cVk+48xMOeZUk/nYpBx5Bw8e3FQ/s518Eurr6wvlhKOc1MvGZfyX53khrtQhl8uFuDKui/kOJZdPKn19fbp8yWUISCv7snzKSrlIJKK5fGqRXBmbyCdqWQdfr1i+2TecF9yKfjbfYdulTaSNWH6zvsnn8wHfNOVo856eHmvfsO838x0z3+FWfaenp0f3vc33Dx06BCCYT9Dm+wcOHGhpX/q25Mo3uN3CVp8EzwI4scnnexRWSmUB/NnGf3/B9/0V3CXQkeXNwgNo4vG4TowpzzY1BxmTyzgr25nArQZBKccOlmevSjnyGczbTD85yEgZAOjt7Q2d60s505FN7sGDB3X5JlcmMmVGELMOfvh60kw/z/Os+sXj8Zb2leXb7NuKS048HtevvVKOr3Ob+c5W9LP5DttunutLeTl42OzGODx5JrSUo090d3db+0b65t3wne7u7tAgKGX4ACDLl3Isn5mAmtmXvmP+uO02tjQI+r5/3ff9y5t8lgFAKdUB4LsAPAAf8H3/7hwRtYH19XWk02mk02m9apVOpxGLxXSsF69JOQajAghx0+k0otEoVlZW9E1j1uF5no5pMrnpdDrwBGLTL7KR4kke2CTlWH5kI6VTqVQKyHieh5WVFb2axxU9ytE28Xg8xM1ms1hbW8P6+jqi0WiIm06n9Txbd3d3iFsqlXQcG1diJZdtX1lZQTweD3FTqRTW19f13JzZ9nQ6rctPJBIh+zLFE4AAn3KJRELzh4eHQ/blE976+jpSqVSob+gbKysryGQyIf1k+TbfSSQSWFlZge/71r6JRqO6b21909nZGfBN07784TF9k3L0LT71mnXYfFvKcb4yGo1a9WP5nudZfZ+gb5k2ou0PHDhg9X3WH4vFQty7ESe4q3OCGwPgXwDIA3gMQEYpVdj47PoQXiwWEYlEsLy8jOXlZVy+fBlAYzJ8dXU18MtBGcoVCgX9KmRyOVmdTCZRKpVCXOD2BDKAEHd5eRm3bt1CpVKB7/tW/XgTMKLerINtY6dfvnw5IHP8+HEkk0msra3pa1IOgH5tNLnA7bx99Xo9xF1eXoZSjcSoN2/eDHEvX76sf7kZGCu5bHsymdRbAyX31q1biEQiSCaT6OzsDLU9kUjo8qvVasi+1WpVh4lUq9WQfavVquZPTU2F7DszM6NPkrt161aob7hrI5lMBvyLcvIJzOY7PT09SCaT2rdM+9brdUQiEYyNjVn7ZnFxUfuWrW/W1tZQqVT0wopZh1yYsPnO8PAwkslkgC/lfN9HJBLB0NCQVT+Wv7a2ZvV9+o70TSnHKQ1mpzbtWy6XtW+Z3EgksushMrsdLH0aIvDR+JzejO+CpV2wtOeCpV2w9I85WHpXQ2R8338WQPggVwcHB4d7FG2dRcbBwcFhp2j7BArA7cBeW7Q5J4dtEeicj+N2HJML3F4tNLmUkyEXkssVPIZENNOPW8hs+rUqn2VuVj5XMVuV36r9iURCz/vZ6uDCjY3L75txgcbChswis9vlcyLerJNoZTugkU3HZjvaxGY7Gezcqny2zVY+MwO1o2+yT3bqmyynVf/vBtr6SbBYLOrVukQigYGBAQANJ+zt7UWpVMIf/dEf6WtSrlAooFQqYXR0NMRNJBI4fPgwzp07FyhfyvX396NUKuHYsWMhbiKRwGc+8xlUKhUcPXrUqt/g4CBKpRIef/xxq35sG0+bGxgYCMgcP34c586dQ39/f+B4QsqNjIygVCqhr68vxAWAL33pSyiVSigUCiEuy69UKnj00UdD3IGBAZRKJZRKJR3nKLlctDh37hy6u7tDXK74nTt3Dp2dnaG2c0W4VCohl8uF7NvZ2Ylz586hUqkE+JTL5XKaX6lUQvb9uZ/7Ob33lqu5so7u7m6tXyQSCenH7+kbpn49PT04d+4cjh8/bu0b+t7Y2Ji1bz7xiU+gVCphZGTE2jf9/f2oVCr4/Oc/b+2bI0eOoFQq4ejRo1bfGR4exrlz5zA4OGjV7+jRoyiVShgaGrLq9/nPfx6lUgn9/f1W3//CF76ASqWCw4cPW31/dHQUpVIJjz32mNX3eQpkb29viFsqldpv7/CdYDt7h5eWlgL7GPk38wmWy2U899xzAV42mw3lEzR/8eLxOH7iJ34CL7zwgl4Vk3XIfIKe54V+pcbGxjA9PY21tTWdT1DWwXyCw8PDuHjxYkg/lsutS6Z+qVQKDz74IP7+7/8+cFob5WQ+QfPXFmicd3vp0iW9dcnUL5PJQKnGIedcYZV1EMz5JrlKKZ2A4e/+7u9CT3rvvPNOy+/NfIJmzsVIJIKf+qmfwsTEhF5pl/p73u18golEAjdu3Ajo39fXh8XFxjZ25hOU+jOfYKVSwYULFwKpnrLZbCifoNk3sVgMDz/8ML7//e/rfIJSPz4BcZXV7Jvjx49jZmZG+6bZN/S948eP44UXXgj1DbfF0fdM/ZLJJB566CE8//zzgZMEKSfzCUaj0ZB+4+PjmJycxOrqKmKxWOj7hx56CJcvX27qW/RNJg829fc3tv3ZfCubbSRTaLe9w3cVMt+dPJ+U+yrpRNKRKSfPVJVcAAEHtHFZvi23HgDtXHJvq00/yjXTj+Wbef2A26E2vK6U0nJmvZKrlMKtW7d0+XQ605bRaFRnEJZcc5/odvTbSt/wO7Nv/I2wI/5tsy//tuWsk8dl2vp/q77TTD+2XdrE7BvuUbb1DfdrUzdTjjbnkaxm35i+uR3fke1t5Tu2tvNHv5Xv+76vB0jTvgy/oZzN1ruJtn8SdOcOu3OHzfIlh2VJ3dy5w+7cYYm2fxIE3MKIWxhpXr5bGHELI5uh7RdGpNEkuJ3r1KlT1u8LhQKy2axO22MikUjoyfNm38ttQCaGhoZ04kobmJDy2LFj1u/ZNnnDSnDhQt5wEtzixC1iJsbGxpDNZpvuxeScHffPmuDcjEwrRXDHQKVS0Tec+X02m9ULGybMRK0mPM/T5Tf7nnybfTo7O3VSVT7JSMikr3KQJbq7u3X5NvT09OikqjbE442krWNjY9bvmbS1mW/RN48fP279nr7VzPeGh4f1jh4bmJR1aGjI+j0TATfzPSb8bfY9k6ryACUTTDZr861stg2Tqt4J3Ovw7Trc67B7HW5W/l775rvtdbitnwQdHBwcdoq2HwS5yZ0b3oHbqaqSySQGBwf1NVMumUzqVEFmyqVoNKpfB2zcWCym+bZ0TXyNlKmIpBx5fCVopZ/ZTj7lDA0NhdI1UU7qZeMODAzoMB2Ty9g4AOjq6gpxk8mk/vB1t5l+ZromyiWTyZb2Zfme51nt24pLTjKZ1K+kUo5PkZv5zlb0s/kO2y4TQcj6WH4zux08eDDgm6YcXxMPHz5s7Rv2/Wa+I33zTnzn8OHDSCaT+h4wbdTf3w8gmErL5vtMudXMvuRLbrNX+J2g7QdBmVTUnPCNRqPWxI2Ui0ajmi+55DPXn43L8vkx5To6OrakH+Wa6cdXAVkOX0Hy+XwoaalsjyzH5B44cECXb3KlDkytZNbBj00/3iz5fF6/WkkuB4JW9mX5lJVykUgkkIfR5JLTLOknbySWb/bNnehn8x223Ux4y/+z75v1DV9X+Z0pR5vbEt6yTzbr/534Ti6XC9Rj2ohTIbJ8Kcfy+QPVyvdt3N1GWw+CxWIRxWIRqVQKqVRKLyakUin09PSgWCzid3/3d/U1KTc4OIhisYhjx46FuKlUCoVCAWfPnsXJkydDXE7qFotFDA0NhbipVAqf/OQn9WltNv36+/tRLBbxq7/6q1b97rvvPhSLRT1JLtuZSqVw4sQJnD17FoVCQV+TctTr6NGjIW42m8WTTz6JYrGI3t7eEDeVSuHYsWOoVCp45JFHQlzavVgs6kUCs45KpYKzZ8+iq6srxOW/Z8+e1fXLtufzeV1+LpcL2TeXy+Hs2bOoVCrI5XIh+5JTLBbx4IMPhuzLHSOyXPl3V1eX1q+joyOkX19fX8D3TP0GBgZw9uxZ7VumfXt7e1EsFlEul61987GPfUz7lq1vCoUCKpUKHn/8cWvfUL8jR45Y+2Z0dBRnz55Ff3+/Vb8jR46gWCxiZGTEqt/jjz+OYrGo7wHTvk888QQqlUrAN6XcsWPHUCwW8eEPf9jq+6dOnUKxWERPT0+IS7vvJt4VCyOEbSn93XDtXtHDtWt/t+teaauLE7SAq0s0klzh83132pw7bc6dNrebvulOm3NwcHB4F8ENgg4ODvsabT0IcsLU87zALgEGnGazWTzwwAMBGcoNDg7qHSMmlyEWlUoFJ0+eDHG5k4N8k+t5HkZGRlCpVPTBOGYdnGRm9L1Zx3333acng23lMyqfYSRm+dwxcvToUWv5999/v476t7WfO0aOHDlibT8/5u4OrhhyRwdPDZNcuSODJ9JJmXw+r2W580eWH4vFdPmxWCxUvuSwLFlHb2+v3jHCVUhZB3d0VCoVdHR0hPTjYhXlTf0GBgb0jpFWvsmdEWb5Q0NDyGaz2nea+SYP5DLroM/QB8w6RkdH9W6mVr4zMjJi1Y8HH/EeMGXGxsZa+ibvGd6DZvmnTp1q6pvZ7O7vGGnrQfDll19GvV7XMVvDw8MAoI1fr9d1RL2MbRseHsbs7Czq9TrW1tZCXB6Ss7CwgMuXL4e4APScCU80M+uIRCKYmJiA7/tW/YDG3Azns8w62DaCh+PwMzk5iYWFBR3ywTLlITr1eh0zMzMhLoDAgT8mlwc4TUxMIJPJhLjDw8Oo1+uo1+t6jkZyE4kEJiYmsLCwoGPiJDeRSKBer2NhYQGLi4uhtlerVV1+LBYL2TcWi2FhYQETExM6Vk3aNxaLaT63Dco6uru7MTExgXq9rnePyDqi0ajWb319PaTfjRs3dPk235mbm8PCwoL2LdO+kUjjgKtLly5Z+4b2Yd2mfZVSmJiYwOrqqrVvTN8y65iamsLCwkJAZynn+z7q9Tqmp6et+q2urqJer0MpZfX9lZUVTExMBHxTyq2traFer+t5b9O+Fy9e1L5lcuv1eij91k7R9gsjMt/Z1NSU/pv5BLPZbGgyf2pqKpRPUHKBRpBzLpfTRy9KLievmU/Q5AK30x0ppXRqISnHmw0ILzZMTU1hbW1N34xLS0uYmpoK6MFQEd/3A9vOKCfzCZpcAProRZnzTcpxcnppaSnEnZqa0rpzW53kMp4rl8vpwUJyq9UqMpkMcrkcVlZWQt+b+QRN+0YiER2LVqvVtK0pxydtAJiZmQnZl/kFo9Go3joo6+DAncvl9Pem/jKfoKlfLBZDLpeD53mBfIK0kcwnaOubarWqj+VkOi8px7pjsZi1b9gGLqqYdbBtAEK+vbKyEsgnaNMvFoshGo1idXXV6vvUb319PeBblKNv3rp1y+r7vu8jGo3qfIKSy5Xi3UTbD4LyF0/mjuOTEJNnSkemHLn85ZNya2tr+tfSxuUqolydlHJ0fpmPTdbBFTHKmXXwSZNysp0EnzZkTjjKyXxyJlcphcXFRV0+VxNNW8ZiMWvOOjloNauD+pk56yjHJ61WfcPyzb5ZX1/XXFm+jcOVVFkHfzSa9T/5W9HP5jtsu7SJ2Td8Crf1DZ/SZT5BKbe2toZYLIa3337b2jf0mc36ppXvSJuadbz99tu6Hq7USxm+fbXy/fX1df1wYNqXAzD71mbr3cS7Ik7QJVBwCRRk+S6BgkugYKJVnGBbzwk6ODg47BRt/zoMuKSqrcp3SVVdUlWT65KqBtHWT4IuqapLqtrsexnCY8IlVXVJVSXcnCDuzXkXNyfo5gSblb/XvunmBO8xpFIpnYuMv3zMucaMFbxmyqVSKZ3zTXLJ5y+OjRuPxzXf5EYiEf0EJnPKSTnymXutmX58UpPt5A1eLBZDOeEoJ/WycYeHh5FKpXQ+QlOOTwk9PT0hrsz6wV9ryaXTFovFUD5Bhl8wI0irtktZKed5nuZy0JRccvhEY9bB9Gos3+ybaDS6Zf1svsO2y1x80kb0zWZ9c+jQoYBvmnL0icHBQWvfSN+8G74zODgYqMe0EZ8gZa5LKcfy+ZbRzL7kS26zp9udoO3nBD3vdj48OU/Gm4U3gXzt8TxP/6qYOdooJ+eiTC6/598mF0DghrTpxwh5vvI008/3fayurgbKIbLZrP7VNu3Bm1u2U4K7GWSYg03XRCIR4sqnQVsd/Jf6ySMcpX6t7Guzsyxfck2dJYeDtKwjHo9jZWVl274jg8RtvsO2m3bn/1k+n8BN+zLY32yb9E3gdq5HCfqM7KNmvtNMP/om22lyuZNFvkFIGZnItpnvyx9aUwfqT77N1ruJtn8SlDtGuN1I7hhZXFzU16RcPB4PROVLLn/B5+fndWeZdci4KpMr51rkjhEpR93X19et+lF/AHp+Rsp4nof5+flAVL6UY6xbPB4PcTnZLXeMmHXwFSWdToe4Q0NDgR0jJpdPIPPz84EdI5Tjjoj5+flA2Wx7Op0O7Bgx7RuPxzE/Pw8AOkuxtK/cMcK5X1kHk6XKHSOyDgapz8/Ph3bMSP3pe6Z+iUQC8/PzgR0j0kaRSETb3tY31J++adqXg0a1WrX2Ddvm+77Vd6LRqLafTT/G4zG7ulkHd/Qopay+z/g/ee9IOfom535N+/Le8jwvxL0bcYJtPQjydWtlZQUrKyuYnp4GcDs9E3ct8JqUKxQKet+oyV1ZWYHv+8jn8yiVSiEuAB3VHolEQlzOoVUqFSilrPoBwczRZh1sGzE9PR2QOX78OPL5PHzfD9RJOb7W9fX1hbhAY2DlrgSTy6ekSqWCW7duhbjT09M6OzAHS8mtVquoVCrI5/M6HZPkckdEPp9HZ2dnqO2JREKXv7q6GrLv6uoq8vk8KpUKVldXQ/ZdXV3VfOablHXMz8/rRRvOeco6OAAwM7apX3d3ty7f5js9PT06c7Otb9bX1xGNRjE2NmbtG+rPuk37+r6vFx5sfUPf4mBp1jE8PKx/CGz6ccfP0NCQVT/2j+/7Vt9Pp9OoVCoB35Rynudp29t8v1wua98yudFo1C2MSLiFEbcw4hZG3MLIVvCuXhhxcHBw2AncIOjg4LCv4QZBBweHfY22DpHhyVovvvgigEak+uTkJNbW1tDV1YVDhw7hmWeewenTpwPzDidOnMDi4iI6OztRr9fx0ksvBbgA0Nvbi/Pnz+NjH/sYLl26FOBOTk7qye9bt24hlUoFuADwxS9+EV/96ldx/fp1PbEr6xgYGEBHRwd+7dd+DZ/85CdD+nFea2FhAW+++SaOHz+u20mZc+fO4Rd+4Rf05L/neVru2LFjel4qm80GuJ7n4bvf/S4+9alP4fr165ibmwtwgcaugI6ODvzMz/wMvvCFLwS4cqfC1atX8fbbbwe4kUgElUoFTz/9NH72Z39WZ2Mh96WXXsLo6CjOnTuHn//5n9crlWzX9PS0XuXk39K+qVQKf/M3f4MzZ87ghz/8oV6NpH0TiYTmFwoF/PVf/3XAvo8++igmJyexvLyMl156CfV6PdA3uVwOfX19OHfuHH7yJ38ysG3wxIkTuHHjhj4z98UXXwz5TqFQwHe+8x18+MMfxuTkZKhvenp6cPDgQT23aPbN5z73OXz9619HtVrFq6++GuqbQqGA/v5+/MZv/AY+9KEPhfpmYWEBuVwOy8vLeOONN0K+UywW8Y1vfAMf/OAH9dyarIMr9EwZZur39a9/HZ/73OcwPz+PfD4f8v0//dM/xRNPPIHXX38d165dC/TN2toaRkdHEY1G8Yu/+It46qmnQr7POdE333wTN2/eDHBPnDjRdCfNdtHWg+DLL7+s88cBt/eJep6HhYUFLCws4NOf/rS+RqTTaUxOTgYGD8kFGjnnzpw5g8uXL4e4ADA3N4e5uTkAjVVgs46vfvWrmJiYCFyXdVy9ehVXr17F008/bdWPXHlNykxOTuLMmTN6ADPlOPBSP3NB4zd/8zcDN6hZxyuvvKJlTe5W9JuYmMCZM2cC/UM53/e1/ouLi6Hvq9Wq1s20G9DY833mzJmm9pV8k0vbSf3N/l9aWtL6yUUjljc5OalzEtp8Z25uDmfOnMErr7xi7ZubN2/i5s2bAOx9861vfatl39D3nnrqKWvfTE5O4s033wzVS0xNTeHMmTO4evWqVT+56GDT76mnntL6MdhfyjzxxBMtff/VV18FAPz5n/95iGv6lmnfycnJXc8p2Parw/Pz83r1MplMolqt6hUlxphNTk4GcpYxFooraIy3IhdohBiMjo7iRz/6USAxJuW4zM8nNskFGlH18/PzqFarehVa1sG08IcOHcKPfvSjkH5yVZMxVXKVNp1OY2hoSD/JAI1VU8oxyJntlFyllG4b490kF2isoMbjcaTT6cANRTlidXUV6+vrIf06OjowMDCAf/zHfwwkmKUcMwa//PLLgSetZDKJ1dVVHeTM5LfSvkoplEolvPHGG1haWgpk565Wq1BK6V0VsVgsIAMAXV1dWF1d1eE1Zt9wV8bw8HBAf8rJ1We2ReoXj8cxMjKCV199Va/yyrbLg+sZBiLte/jwYdy8eVP7ptk33NHS29sb+LGiHFGgT/QAABTCSURBVFdR19fXsbq6GuqbVCqFo0ePBmwv6+BOEvqPqd+xY8dw5coVHUpk+n6xWMSbb76JW7duad+UNqJv5vN5zM7OhnyfcZqrq6s61pJcxnnu5uqwTty42x8ACsB3AfgAHtkK58EHH/TvBP39/b7nef5GHb7nefr/Sinf8zy/s7MzIEO5VCoVkJd/A/AjkYify+X8jo6OENfzPD8SiQT+NevIZDJ+LpfzI5GIVT/yMpmMVb9sNtuy/Gw227J8/s12muXTLuSbdXie5+dyOT+ZTFrbzw/tbNaRy+X8XC7nK6Wa6pfL5fxEIhEqPxaLtSxfKdWyfMlhWbKORCLh53K5TX0nl8v50Wg0pF8ymWzpO8lksmX57NOOjg5r36TT6XvaN6l3s/I7Ojp2xTfZtya3v7//jsYJvzEgPd9s3LmbCyOfARA+iNbBwcHhHsJdGQSVUg8B+FUA//3dKJ9wp8250+bcaXPutLmdYtcXRpRSHQC+DuCM7/vXbPnadgtcGLGBcz4/+MEPrN/Pzs62TNJYrVZDk//m95zvsU3UTk9PY3p6Wq9cmuA8CyeJTci22crnxLTcbSHBA3pmZmas31+6dKll+995552W7ZdcM2+c7/uaa9vr6fs+lpaWmpZfrVZb6ra2tqa5cmVRft+Kv7i4iImJiaYyPNyqmX43btxoWT4XLuRcsgQPt5JRBxK2w6EkNvNNOYdn8x0exmQeoETQZ+XimsTk5GTgQDMTXE1u5pu0CxcmTVy8eLFp+5eWltritLk/BPBd3/fP34WyreCNQMPJrU++74e2zS0tLQWSicqbRm5NAhDackU5XyysmFwAetLY39iWZurHiWfurTXrYBaaZuXLLDO28j2xkmpuTVpaWtKZPsg36yC/Wq1a20/YuOb3G3PEVv1sbaftm5Uvf1ht5cuFi+Xl5VAfcjFDlmnzHQBW35E3vs135Gl0tvLZ9zxNzixf/qjci77J1dpm5dP+O/VN9q2Nu5vY0uuwUupJpZS/yee0UuqjAB4A8K+2qoBS6hNKqeeVUs8z5OROwOP7IpGI7hyu7qXTadx33336minHDCkmlznfmFXaxmX5yWQyxI1EIjhy5AiA2yu9Zh2JREKv8LbSjzeEbCdvglKpFMoJRzmpl407OjqKdDodyCco5TgIHzp0KMRNp9P6wx8TyaWjlkqlQE45yvGVq5V9WX40Gg3Z1/M8zfU8L8QlJ51O6z3Gso6urq6ALmbfxGKxLetn8x22XebikzZi+c36pq+vL+CbphxtPjw8bO0b+hYH+ma+I3Nd3onvDA8P6763+f7IyAiAYD5BKcfyDx8+3NK+9G3J3e0YQWDrT4JnAXxtE5kfAXgMwEkAS8Zr8H9QSj3n+/5PmyTf958B8AzQCJHZoj4a/NXi38DtHHMcbHhNckx5+Te/oxM14/LXzuQC0I7aTD/Toc06+Iss67PlnTOv8/+bceUNIPUzdY3FYiGufFKzcfkv0z41028z+0pbSTlu7uffzewL2HM9MjvPZr6zXf3Y9s36hk9Mpn3lD1OrvrHlepQDCp+Qt+s7zfSTvmNrO++5zXzflutR3lOyL0xb7ya2NAj6vn8dwPXN5JRS/yOA/8W4fBHArwP41h1rtwXIw7cZj8fXhFqtpgOaZawX5Wq1mn5kl1yg8crAaHcbl/na5CuJlHvrrbcA3I6jM/VjvjjuljDr8Dwv1DYpAwDXrl3T5QANB6Lc6uqqdi6TG4lEcP36ddRqNc2XXOodj8cDmUBkHQTbLrms99q1a1hbW7PqV6vVWtpXZvgx+4Zl85ppX8nhQCC5nC+lLr7vW31nK/rZfIf6yb6Xbec8aTQatfbN4uJiwDfNvmGexRs3blj7hjbZie/wR9im340bN/RrrM2+DCSX5Usb0Td5Jrhp31gspss37St9b7ewq3OCvu9fAXBFXtu4IV73fd++ArBD1Go1/esgT6+iwWZnZ3USUYJPX/JEM/Pkq/X1dczMzAR+7cw6OABks9nQqVncKsbkmSaXN8Jbb73VVD9OUHNyWzqA53mYmZnRzmrKMQjWxs1ms3jzzTf1tUgkEpJjkDJzE0qutF2zOoDGBL98UpZytVoNMzMzoUSZthPIbNfkgo/NvtSFK72Sy/pqtZp+2pBcDoozMzOh0+a24juJRAIzMzMB35Rtp2+SZ/aN1N/WN/S7GzdubKtvotEoZmZmAr4p5WQyX5ObzWb1IAjcTucv65QPHra+oW8uLy839X3Zflu/7ibaOoGCO23OnTbX7HsZwmPCnTbnTpuTaPttcy6pqkuqapbvkqq6pKomWm2ba+snQQcHB4edwg2CDg4O+xptnUoLQGDepFwu49KlS4HVpvvvvz90jXJyvsrGHRsbw6VLl9DR0dFSrtU1OS90p1wZHzU2NoaLFy9qmTu59p73vGdbXAAYHx/fNnczOdpXzktu1Ubkynm57dp3szpsXDlfZXLvpG8efvjhbdu3XC5v2yd2S7+77fvN5HYTbf0kWCwWMT4+Hoo7ikQiOHToEMbHx/GXf/mX+pqUKxQKGB8fR7lcDnEjkQgGBgZw4cIFlEqlEBdopMoaHx/XAcFmHb/1W7+lJ69t+o2MjGB8fBy/93u/Z9WPbRseHgYQDDyNRCI4fvw4Lly4EDiAW8qdPHkS4+Pj6OvrC3GBRs46fm9yI5EIHnjgAVQqFXzkIx8JcWOxGMbHxzE+Po6enp4Ql3uHL1y4gO7u7hA3kUhgfHwcFy5cQGdnZ6jt/H58fBydnZ0h+2azWVy4cEEvbpj27ezs1HzaT9bx/ve/H5VKBePj4zrmTdbR3d2t9ZM8yvF7+p6pX09PDy5cuIAHHnjA2jd9fX0YHx/H2NiYtW8+9alPYXx8HCdPnrT2zeDgICqVCr7yla9Y+2Z4eBjj4+N64cSsY3h4GBcuXMDIyIhVv+PHj2N8fBxDQ0NW/b7yla9gfHwcg4ODVt//xje+gUqlgoGBAavvl8tljI+P4+Mf/7jV9/k9A/Uld3x83C2MSAwMDODmzZuBnHAMQeAE9sDAAF599dVQTjjCF/kEZUxUNBrFsWPHQvs4KccN37Y4LqCRbPLtt9/GyspKIJ8g5bjx/+DBg5iZmQnp53mNXHP8MFyFyGQy6O/vx9TUlJ6kZziBzCdISC5vBOaEY6iErCORSOiP3MlDOWJ1dRW+7we4Sil0dHSgUCjglVdeCUzoMzdgIpHA4OAgXnvttUDYA0Nm+JTOUBJp30gkgmKxqPd/8zrtC9wOtYhGo3rrHJHP53UsoMxTZ/rO4OBgQH/WIRcEuAVP6sfM1jxe1Owb5hNk3jyzb3p7e/HWW29p3zT7JhaLIZVKoaurC6+99lqob2Q+QZ49beaiHBgYaOo7chfS+vp6SL/BwUG9IGnGswKN1ecbN25geXk5kE+QcvTNXC6Hubm5kO9z4UnmEyQ3mUyiq6trVxdG2v51mEGl/BtoDGwM5GRCUDnYyxuL35nBtWtra3qDt8kFGrFmrIM3hJTjweMcJEz9yGeGYbOOVCql26aUCrSTmJ2d1UHbpj3o3Gynyb127Zq+LvXj30zGKVdfTRsACPwIUI7/8kYx9ZMxnLa2szyWbwt8JlfuHbZxeBPKOt55553QwGXznWb6MZiX35n6se2y76V9GQTMYGnTvgsLC6EypRzbdP36dWvfyDZs13e4OmzjMtCe9Zhtn5ubQyQSCfS9bA99k7G0zexrBrLLv3cTbf067ODg4LBTuEHQwcFhX6OtB0GXVNUlVXVJVV1S1Z2i7RdGZGLGVnsn2/navaKHa9f+bte90la3Y8RAJpMJHJrCDfDRaBSZTAYnT57U1+SnVqshk8kglUqFuL7v66V6Ts6adcTjcWQyGb1aaNYh8wna9EskEshkMhgaGrLqt7a2hkwmo7df2b4vl8uIx+NW/VOpFDKZTGDxQ8qVSiVkMhnEYjGr/ozT6u7utuqfyWT09iuTS5TLZXieF+Ky38rlMlZWVkL6LS8v6/IjkUhIP6UUyuUyAOjtYbKOSCSi+TIpAj+8qWRiWVkHn1bK5bLeWig/1WpVl2+zfbVaRblc1seL2nwrk8k09a1CoRDwTVN/rnyPjIy09C1uDTTrqNfrKJfLSCQSVv2SySQymYxeODHrGBkZ0ba16Tc6OgoA2rfMOuibvb29Vv3oX0z8avrdbqPtV4dlPjn5Lz/cdyjDRWxytjK2ywXC+emacaVcszrMvwnP80LXm+nVitusPrbD5G5F163ot5l9W9lPcrfKIbbaN9vVbyttl3t7TfvS5pv1DeuQ2E7/74TbrO1bsW8z35d22UyX3UDbD4LValWHOTAmS+Zsu3r1qr4mOUopHSO2vr4e4AKNMIQrV65YucDtPIH8ZTXlrl+/HpAz9WMMFOXMOnieq+SaEfNXrlwJ5ayjnEwTZXIjkQiuXr2qU+czzEXK1Wo1JBIJLCwshLgy4QHbLrl01CtXroTyCVKuWq22tK9MIGD2je/7mivDMGwcW65HJh2gLnx6oxz5W9HP5jtsuxnbyPr4hBWLxax9c/PmzYBvmn3Dc5mvXbtm7Ru2QdrlTnxndXUV9XodnudZ9bt27VrAd8y2MyxN5hOUNqJv2nJpVqtVHUtKvmlrPsnvFtr6dbhYLAay3HLCmbsWEomEnsSWUekM1GUwsMlllHp/f7+O2jfr4M6GeDwe4kYiERw8eFBP3Nv0I7+np8eq33333YdEIqEn+c1svidOnEB/f3+gfClHvY4ePRriZrNZHD58GIlEAtFoNMQlv1Kp4MCBAyGuDKTmIodZR6VSQX9/PzzPs+qXSCTQ39+P3t7eUNvz+bwuXy4kUCYajaK/v1+n6jLtSw7TzJv27ejo0KmebH1Dfn9/v86IIuX6+vpa+s7AwAD6+/sDQcey7fRNvpKa9u3q6tK+ZesbLgyZtqMcfauZ74yOjrb0HfJHRkas+vX29gbqMe1bKBRa+j5tx0Ur076nTp0K+KbJdQsjAm5hpL2v3St6uHa1V1t3e2Gk7V+HAXfanK18zsswi7M7bc6dNke40+aCaOvXYQcHB4edwg2CDg4O+xptPwhyIpzR5AACk+k8Y8SMSidXRuVLOZ7j0IzLaHlOvJtyPHvV3DFCOUb1M6aqmX58LZbt5CuBGZUv5aReNq65Y8SU4yvJwMBAiJvNBneMmFxmgDF3jFBO7hhp1XYuQJn25Y4UAHrHh+SaO0bMOnhuCnUx+0buGNlMP5vvsO1yx4i0EX2zWd+YO0ZMOb7KmztGpM+Q38p35I6RO/Edc8eIaaOTJ08CgPYt00Ysf3BwsKV95Y4RyvHv3URbzwkylx0NUywW8fLLLwdkPM8LHehiynV2dlq5/E4eBmSTs13L5/P66EGet2GTy+Vym+rX7Brr2Q5XDi6t5Pr7+/Wpdza5VrYHgIMHD7asY7Pvs9ls07I5ELXS33bt2LFjW9b/0KFD29avo6MjkJB3q/r19vbq4yjpeza5dDq9qe/kcrmm+h04cGBTXZrV28pucnBsVU6hUAjNMW/Ft3jf7xbaehD8/ve/v9cqODg4tDna/nXYwcHBYSdwg6CDg8O+hhsEHRwc9jXcIOjg4LCv4QZBBweHfQ03CDo4OOxruEHQwcFhX8MNgg4ODvsabhB0cHDY13CDoIODw76GGwQdHBz2Ne6pzNJKqTkA03utRxMcBHB9r5XYBbh23Ht4t7TlXm7HkO/71swL99QgeC9DKfV8s/Tc7QTXjnsP75a2tGs73Ouwg4PDvoYbBB0cHPY13CC4dTyz1wrsElw77j28W9rSlu1wc4IODg77Gu5J0MHBYV/DDYIODg77Gm4Q3CZUA99VSvlKqUf2Wp87gVKqSyn1b5VSl5VSt5RSryulvqyU6t5r3bYCpdS/UEpNKaVWlFL/SSn1vr3W6U6glPo3Sqm/V0otKqXmlFLfVkrdv9d67RRKqcc37oc/2Gtd7gRuENw+PgNgba+V2Cb6APQD+NcAygB+GcD7AXx9L5XaCpRSvwTg9wH8DoD3APhbAH+mlBrcU8XuDKcB/O8AfgrAfwagDuD/Ukp17aVSO4FS6r0APg7gB3uty53CLYxsA0qphwD8RwAPAngTwH/j+/4391arnUEp9QEA3wFwwPf9xb3WpxmUUv8vgB/4vv9xce0lAN/0ff/f7J1m24dSKgtgAcAHfd//9l7rc6dQSuUAvIDGIPhZAP/g+/6/3Futtg73JHiHUEp1oPHEdMb3/Wt7rc8uohNAFcDyXivSDEqpOBo/PH9hfPUXaDxVtSs60LgX5/dakW3iGTR+hP7vvVZkO3CD4J3jDwF81/f983utyG5BKXUAwG8D+He+79f3Wp8WOAjAQ+PpW+JNAIWweNvg9wFMAHhurxW5UyilPg6gCOB/3mtdtgs3CAJQSj25MaHb6nNaKfVRAA8A+Fd7rbMNW22HwckA+DaAK2jMEbYDzDkcZbnWFlBKfQHATwP4r33fb6s5ZqXUfWjMzX7E9/3aXuuzXbg5QQBKqYNoPGW0wo/QmMz+7wCsi+vexv+f833/p++OhlvDVtvh+/7yhnwWwHk0BpF/5vv+0l1WcUfYeB1eBvCo7/v/h7j+JQD3+77/M3um3DaglPrfAPy3AP6p7/uX91qfO4VS6jEA/x7BBUIPjR+kdQAZ3/ere6DaHcENgncApVQ/gLxx+SKA/wHAt3zff/XHr9X2sDG3+WdoDID/he/7b++xSlvCxsLI/+f7/ifEtUkA/2c7LYwopX4fjQHwtO/7L+61PtvBxjTKgHH53wN4CY0nxEt+Gwww0b1WoJ3g+/4VNF4bNZRSAPB6Gw6Af4HGYsgHAWQ2XosB4OY9/mrzBQDnlFLfA3ABwD9HI+TnD/dUqzvAxpPrR9Gw/bxSivOZS/f607iE7/tvAXhLXlNKvYOGD/3D3mh153CD4P7EgwDeu/H3pPHdPwXw7I9VmzuA7/v/YSOo+38CcBjAPwD4gO/792oyXhv+xca/f2Vc/y0AT/x4VXFwr8MODg77Gm512MHBYV/DDYIODg77Gm4QdHBw2Ndwg6CDg8O+hhsEHRwc9jXcIOjg4LCv4QZBBweHfQ03CDo4OOxruEHQwcFhX+P/B5FMA46wIsXhAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(Lx/2,Ly/2))\n",
    "\n",
    "plt.plot(xb,yb, 'k')\n",
    "plt.plot(xb.T,yb.T, 'k')\n",
    "plt.plot(x,y, 'k--')\n",
    "plt.plot(x.T,y.T, 'k--')\n",
    "\n",
    "# Fracture\n",
    "tip1 = [-1,-1]\n",
    "tip2 = [1.1,1]\n",
    "plt.plot([tip1[0], tip2[0]], [tip1[1], tip2[1]], 'r')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-04-20T04:03:39.204806Z",
     "iopub.status.busy": "2023-04-20T04:03:39.199691Z",
     "iopub.status.idle": "2023-04-20T04:03:39.287106Z",
     "shell.execute_reply": "2023-04-20T04:03:39.286440Z"
    },
    "papermill": {
     "duration": 0.112155,
     "end_time": "2023-04-20T04:03:39.287229",
     "exception": false,
     "start_time": "2023-04-20T04:03:39.175074",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[40, 40], [41, 41], [42, 42], [43, 42], [43, 43], [44, 43], [44, 44], [45, 44], [45, 45], [46, 45], [46, 46], [47, 46], [47, 47], [48, 47], [48, 48], [49, 48], [49, 49], [50, 49], [50, 50], [51, 50], [51, 51], [52, 51], [52, 52], [53, 52], [53, 53], [54, 53], [54, 54], [55, 54], [55, 55], [56, 55], [56, 56], [57, 56], [57, 57], [58, 57], [58, 58], [59, 58], [60, 59], [61, 60]]\n"
     ]
    }
   ],
   "source": [
    "fracture_1 = generate_fracture(x,y,nx,ny, tip1, tip2)\n",
    "print(fracture_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.017291,
     "end_time": "2023-04-20T04:03:39.320831",
     "exception": false,
     "start_time": "2023-04-20T04:03:39.303540",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.01589,
     "end_time": "2023-04-20T04:03:39.353171",
     "exception": false,
     "start_time": "2023-04-20T04:03:39.337281",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  },
  "papermill": {
   "duration": 19.596667,
   "end_time": "2023-04-20T04:03:39.477146",
   "environment_variables": {},
   "exception": null,
   "input_path": "__notebook__.ipynb",
   "output_path": "__notebook__.ipynb",
   "parameters": {},
   "start_time": "2023-04-20T04:03:19.880479",
   "version": "2.1.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}