Galerkine_1D.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/usr/bin/env python3\n",
    "# -*- coding: utf-8 -*-\n",
    "\n",
    "\"\"\"\n",
    "Méthode Galerkine point par point\n",
    "\"\"\"\n",
    "\n",
    "DataFolder=\"C:/Users/Julie_000/Desktop/Stage/Stage_l3\"\n",
    "#DataFolder=\"E:/Clément/Julie/Stage_L3\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Importations des librairies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mat4py import loadmat #pour charger des .mat\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt # pour tracer les figures\n",
    "import os # pour modifier le dossier de travail\n",
    "import numpy as np \n",
    "from pathlib import Path # Pour rendre les Path compatibles entre Mac et Windows\n",
    "from matplotlib import rc #Les trois prochaines lignes pour que Ã§a ressemble Ã  latex\n",
    "rc('font', size=16)\n",
    "rc('text', usetex=True)\n",
    "import scipy as sc\n",
    "import matplotlib.cm as cm #colormaps\n",
    "from matplotlib.colors import Normalize #Pour l'utilisation des couleurs dans quiver\n",
    "from matplotlib.ticker import MaxNLocator #pour les courbes de niveau\n",
    "from scipy.ndimage import gaussian_filter #Filtrage gaussien\n",
    "\n",
    "from scipy.linalg import pinv as scpinv\n",
    "from numpy.linalg import pinv as nppinv\n",
    "\n",
    "import fonctions_projet_divergence as fpd #pour les données\n",
    "\n",
    "colormap2=cm.PiYG#colormap\n",
    "\n",
    "from mpl_toolkits.mplot3d import Axes3D#plot3D\n",
    "import matplotlib.patches as mpatches #pour le nageur\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_champ_vitesse(u,v,x,y,title):\n",
    "    \"Fonction de tracé du champ de vitesse\"\n",
    "    #norme\n",
    "    f=2\n",
    "    u1=np.copy(u)[::f,::f]\n",
    "    v1=np.copy(v)[::f,::f]\n",
    "    x1=np.copy(x)[::f,::f]\n",
    "    y1=np.copy(y)[::f,::f]\n",
    "    \n",
    "    velocity=np.sqrt(u1**2+v1**2)\n",
    "    velmax=np.mean(velocity)+3*np.std(velocity)#norme maximale représentée sur la colormap borne sup arbitraire\n",
    "    colors = np.copy(velocity)\n",
    "    colors[velocity>velmax]=velmax\n",
    "    norm = Normalize()\n",
    "    norm.autoscale(colors)\n",
    "\n",
    "    #plot   \n",
    "    fig,ax=plt.subplots(num=title, figsize=figsize)\n",
    "    im=ax.quiver(x1,y1,u1/velocity,v1/velocity,colors,headwidth=8)\n",
    "    cbar=plt.colorbar(im, ax=ax)  #pour la barre de couleur\n",
    "    cbar.set_label(r'vitesse en $mm\\cdot s^{-1}$',fontsize=fontsize)\n",
    "    circle = mpatches.Circle((0,0), 6, color='gainsboro')#patch\n",
    "    ax.add_patch(circle)\n",
    "    \n",
    "    fig.suptitle(title,fontsize=18)\n",
    "    plt.xlabel(r\"$x (mm)$\",fontsize=fontsize)\n",
    "    plt.ylabel(r\"$y (mm)$\",fontsize=fontsize)\n",
    "    plt.show()\n",
    "    \n",
    "def plot_divergence(div,dx,dy,title):\n",
    "    \"Tracé des divergences\"\n",
    "    #norme\n",
    "    vmax,vmin=np.max(div),np.min(div)\n",
    "    norm=Normalize(vmin=vmin, vmax=vmax)\n",
    "    levels = MaxNLocator(nbins=25).tick_values(vmin,vmax)\n",
    "    \n",
    "    fig, l0 = plt.subplots(num=title,figsize=figsize)\n",
    "    im0 = l0.contourf(dx,dy,div,levels=levels, cmap=colormap2,norm=norm)\n",
    "    l0.set_xlabel(r\"$x (mm)$\",fontsize=fontsize)\n",
    "    l0.set_ylabel(r\"$y (mm)$\",fontsize=fontsize)\n",
    "    circle = mpatches.Circle((0,0), 6, color='gainsboro')#patch\n",
    "    l0.add_patch(circle)\n",
    "    cbar=plt.colorbar(im0,ax=l0)\n",
    "    cbar.set_label(r'$s^{-1}$',fontsize=fontsize)\n",
    "   \n",
    "    fig.suptitle(title,fontsize=20)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chargement des données\n",
    "## Choix des variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "prof=15 # Choisi sur quelle profondeur d'eau travailler (5, 10 ou 15 mm)\n",
    "m=5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Constitution des tableaux\n",
    "### Données brutes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(49, 60, 60) (49, 60, 60) (49, 60, 60) (49, 60, 60) (49,)\n",
      "(49, 59)\n",
      "(49, 59, 59)\n"
     ]
    }
   ],
   "source": [
    "os.chdir(Path(DataFolder)) # se place dans le dossier où est rangée la manip\n",
    "manips=loadmat('PIVnagseul.mat') # Charge les données dans la variable manip\n",
    "\n",
    "piv=fpd.PIV(prof,manips)#données pour la profondeur\n",
    "\n",
    "P=fpd.nb_plan(prof)#nombre de plan dans la profondeur choisie\n",
    "\n",
    "u,v,x,y,z=fpd.donnees2(P,piv)\n",
    "print(np.shape(u),np.shape(v),np.shape(x),np.shape(y),np.shape(z))\n",
    "\n",
    "#tableau des z:\n",
    "z=np.reshape(z,(P,1))\n",
    "z_2D=np.hstack((z,z))\n",
    "for i in range(57):\n",
    "    z_2D=np.hstack((z_2D,z))\n",
    "print(np.shape(z_2D))\n",
    "\n",
    "z_3D=np.hstack((z_2D,z_2D))\n",
    "for i in range(57):\n",
    "    z_3D=np.hstack((z_3D,z_2D))\n",
    "z_3D=np.reshape(z_3D,(P,59,59))\n",
    "print(np.shape(z_3D))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Divergence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(49, 59, 59)\n"
     ]
    }
   ],
   "source": [
    "div_2D=np.zeros((P,59,59))\n",
    "dx=np.zeros((P,59,59))\n",
    "dy=np.zeros((P,59,59))\n",
    "\n",
    "for plan in range(P):\n",
    "    div_2D[plan],dx[plan],dy[plan]=fpd.divergence2D_gauss(u[plan],v[plan],x[plan],y[plan],1)\n",
    "    \n",
    "print(np.shape(div_2D))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Manque de données"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "176400 (49, 60, 60)\n",
      "36706 36706 10638\n",
      "0.2 % des valeurs sont égales à 0\n"
     ]
    }
   ],
   "source": [
    "print(np.size(u),np.shape(u))\n",
    "nu=np.sum(u==0)\n",
    "nv=np.sum(v==0)\n",
    "ndiv=np.sum(div_2D==0)\n",
    "print(nu,nv,ndiv)\n",
    "print('{} % des valeurs sont égales à 0'.format(np.around(nu/np.size(u),decimals=1)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "77970 170569\n"
     ]
    }
   ],
   "source": [
    "v_nan=np.where(v==0,float('NaN'),v)\n",
    "u_nan=np.where(u==0,float('NaN'),u)\n",
    "\n",
    "div_2D_nan=np.zeros((P,59,59))\n",
    "\n",
    "for plan in range(P):\n",
    "    div_2D_nan[plan]=fpd.divergence2D_gauss(u_nan[plan],v_nan[plan],x[plan],y[plan],1)[0]\n",
    "\n",
    "print(np.sum(np.isnan(div_2D_nan)),np.size(div_2D_nan))\n",
    "div_2D_nan0=np.where(np.isnan(div_2D_nan),4,div_2D_nan)\n",
    "#print(div_2D_nan0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhAAAAGnCAYAAADmJ27ZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdb2wcZ54n9u9PFEmJMu1mi+y2lxBokbI0A0rJSKYsYwJH6xXpkxAdsLNLWXHyIgFmhpokrxLsSefNiwMOyU2oS4Ak92KXnPhe3JuFJc34cufD+obyOopxznhMy7Nrc2dWHtIWtFov/5hqi7YoSqKevKinWsVidXdVd/15quv7AQipu6urq7vJrm8/f36PKKVAREREFMSWpA+AiIiI0ocBgoiIiAJjgCAiIqLAGCCIiIgoMAYIIiIiCowBgoiIiALbmvQBEBGRGURkAkAewAdKqfNJH48XERkH0A9gWSl1JunjyTK2QERARA6JiPL4uSUiEyKSS/oYm4WIjItIKMVMRGRMv098fyKk37NblS5TovJKqVOmhgcAUEqdU0qdghV0KEEMENE6B6BL/zwL4McAXgbwmYgccm8sIhfDOhnSRvW+tml7T9J2vJQNOiRu+syrsv2oDvQTIjIlIqNRHh/Vh10Y0ZpTSpX0/68CuCoikwDe1j9dru2n4jy4jKn3tU3be5K246UmJSL9sL5ElQCMwefvpg4LV5VSc/pyDsCHItJvcstIFrEFImY6UPwQQE5Ezrpum9RNcxSyel/btL0naTteal5KqTml1Bml1DkAywHumrfDg95PCVYQGQ/7GKkxDBAJUEpdBXAZwKtJHwsRkSl0q8WER3fHVX27724Qih4DRHKmYLVClAfsOQeT6b6/TX3ZXgP9dH/hh/r6D51/ZCJyVkQ+dOx/1nXbrL7fhP6ZdW1Tc9960OiU1zaObe1tbumfKde+Kj6OHyIy7Nj/rIgMu27fNFBPH7+9/QSATYMnPQb8OZ/zhuMVkX7H6zDr1W8bxuup92G/b7f0uIf+Ks/Tvb8xj9fBz+MGeb8bej8b2Yf+G3Eep+fvhGP7Wq/PpoG64v13WPFvzcfzvCUiU47rxqrdp1nplofzAOZcN9mvs/v6MnGMm3BdPqvf35z+Oat/Lrr/RsPYR6YopfgT8g+AQwAUgNEq24zqbYYd140DuFVtHwA+BPCh4/JZvd1ZAMMAJvTlnHOf+n6zAM7q68fs/ev7zept+gEcCrBvpe83qvd5C8Cs65iH9XYX9Xaj+v9n/TxOjdfafQyj+rkoAP1er63rMS/qx7SPfcPjetzPfrxZx+Pdsp+3+xhcx9rw66nfH6W3PaS3m3C8Z+7jHXVs73yeE1Vew0rvo9/tfP1OVnqNQ/6dGHO9t/2ubX2/Pq77jbmPBxX+1mocaz8e/V1M6J9xOD4XYv7suhjRfmcbeU76NfnQdd1F1+Ux/e8t59+Eva2+bsxx3SGP37uG95Gln8QPoBl/4C9A2CdVzwChL99y/pHASuHK8UtuXz7r2veHeHRytj9Mpzy2cX5Ijjo/JAPu2/0clOs+s+7HD/I4NV5rr8crn2RrvLbu18Q+afkJEJues/N48ejEdMjv8/Tzejr263kirfA8x13bDLuOze/76Of4/D7PigEiot8J+29ywnW979fHtU2lALHpb63GsboDzSFU+dxwbXsRVkum359+P/v0e+xBftBAgNC/D7c8XivnZ+MwrL97+3fH/cVrChW+2IS5j6z9sAsjOXZzbMUmOQCTsE4YNrtZ84L+d0j/Oy6OehN63yOufQUtuBJk39OO/88C5ZHTdp9mPyoPgAryOL4oqxn0Mqw/7k10U3gO1rcJp5LH5pVses76MW32+2rPVQ/l9XQ8xtu6CbVis75ust/0PJVSl2E919MBHtfvdmG8n1H8TlyF1Y9u77ue18cv339ryjFYUL+X/UqpSz7ve0opNRLgp9pnjckuAjhV7fiVUpf17UP6svs19PoMOgTH32wY+8gaTuNMzmFg4weIh9cBnBWRUf3LfBrAZfVoaqj9oT2AGqOcPR7nMoAxsaq6LcP60HP+IQTZd7UTb62g5PtxAirB+oP3Yl9f9wdqhedcbX+hvJ5KqZKIjMD6IBsHABG5rJTyOrHaz9Pr8ebw6L2p+bgBtgvj/Yzqd8L9nAO9Pn7Vc6K2x2f4DQ9x0V8ALga4y2VlzboI6/HHYbUQ+T1Jj8B1Qtfhtt99PazP09cj2kcmMEAkQP9RjsIaLFSRUuqqiMzB+iW9BOsDzTlF76r+N1fnt4tlPPr2fBXAsRD3bbPv2w/vE2xYj+NW6fGcxxRnJbvQnqf+MH3W8Xs0LiJn1eY58s7X/qrrtn5sbEkISxjPM67fiSRen030ILw53UoS5H5BTuwAcC7o66m3fzbg44RCDySdChAeAKvV0X1CHwZQcrX25LD58zTMfWQCuzBipn/pLsL6hvxjH3eZADBsj/R1fkPRv8xz8JgOWqHp2WkYwBmllADoUko96/xm2eC+y/SHYgkezboikgvrcVz364f1h+35gez4oHYf00A9j+dHFM9TWfPsz8N6noc9NpmGx2uvf5fs38NQhfE8I/qdOITNvxONvD6hnFS9woNYM4pqtn7oLowgP6npwnC8Lpcd1x3Sf9vVeHUpbGpRgFUReE4pNaf36+zuDGMfmcAWiGj1Oz7w8rBO2uf0/4/5bDK+BKup+lX9f7czAKb0txF7KuIZWB/A1fpi52B9c30dQElEAGDa9S2o3n27nfLYz2n970gYj6OnwI3r+/4EtQPaeVjdQxOwThSHYA2ijFIYz3NMbzuh72eHJfd4Dru744cALur3136e4wAuBfxmF0QYvzdh/k7YfdcbficCvD5f6v2N4dHYmoanWeoTzmkAcyLyKqxvvYdhjYPIxLda/fn4NoAf21+O9OuSh1W51w4MOVhfeCq+945uIPcXh2FsHrvwLB4FgtN2t0sY+8iUuEdtZuEHj0Z8u3/sqUG+RtA7rv8QVWZ16Meb8noMeIwg19dPVDjGD53HF3Tf8Bid7rGfWegPaz+PU+O1HtP7G3a8TptGnHu9tvZ1+v7jjn3VnIVR6zk7fgeGq7wOgV9PPDoZ2lNVZ1Fltom+zvnabJpeGOA51ft+ez3PqtM4G/ydKM8WqfY7EeD1yblur/a7sulvrcpxnnXdd8OsqyR+wnx8/bqNw/pbtz9bNkxRxaMZFu6ZZV4/7hkQm6ZxwjXLRl+/4X1y/G7Z02f7w9xHln5EvxCUITplT8Ea2Wyn/hys0cdTACYVl8mllNID784qq3uOAhCRiyolrR9pOtZmxTEQ2dQPbBpPUVJWk+1lVJ69QEREBIABIqsuwBr3MKVLtfbrgVsT8O7rIyIi2oABIoOUNXhzN/RASjwak5AH8KyKbnAdERE1CY6BICIiAIBYC4BdBfCB2lxXxAh6jEs/rIGLidSoIAsDBBEREQXGLgwiIiIKjAGCiIiIAmOAICIiosAYIIiIiCgwBggiIiIKjAGCiIiIAmOAICIiosAYIIiIiCgwBggiIiIKjAGCiIiIAmOAICIiosAYIIiIiCgwBggiIiIKjAGCiIiIAmOAICIiosAYIIiIiCgwBggiIiIKjAGCiIiIAmOAICIiosAYIIiIiCgwBggiIiIKjAGCiIiIAtua9AEkYev2Lart8Uw+dWpCjz3WityOHLZtfRytajtw/yEe3r0feD8P7z2I4Og22tK2+e9uy7ZWqHbBA3UXt+8u4fY3K7h792Hkx0LmWF24v6SU6oli351929R6A79Pqwv3/71S6niIh9Q0MnkWbXt8K545XUj6MIhC8cLRp3By6HsY3DmC4vogtvz9HazMzAfez+qNxQiObqPtuzafIzoHi1jra8Xyg2v4+W8n8db7V3Bt5k7kx0Lm+Kt/cfN6VPtev/uwoc/7v/oXN7tDPJymwi4MIiIiCowBgoiIiAJjgCAiIqLAGCCIiIgoMAYIIoplACURNRcGCKKMY3ggonowQBBlGMMDEdWLAYIooxgeiKgRmSwkRZR1DA9E4RCRHICX9cUBAFBKnauy/UUAEwDmACw7b1NKlURkFMAIgIv66lMAxpVScyEfesMYIIgyhuGBKFTjAM4ppUoAICIfishZpdT5CtsfAjDlcf0pAJcA5AEMAxgDcBXAD00MDwC7MIgyheGBKHRDsE74tjkAh6tsP6GUEucPrAByyd5AKTWgb3tWKXU1ouNuGAMEEcWGAYaajT7JX3JcVamFwebcFjVaK4zGLgyijODJm6gu3SIy7bg8qZSa9NpQRM4CuFzpdgBwdkeIyDCAyx77GYM1PiKv71Nxf0ligCCiWK3eWPRclZMoCi1bBflCSyO7WFJKDVXbwDWQcjbAvs8opU65rpsGULKDhohMiMioq5XDCOzCIMoAtj4QRUcpVVJKTequiMN6pkVVerbFBx77uuoaNDkF4NXwjjY8DBBEGWDaN35noFmZmce2ZUFHSzcO9h7DgWf6Gv3GSBQLEcnpbgunKQCjPu5+BtYsC/c+h11XlWCNqzAOAwRRRpgcIm7/4gZyt3ZiX+cJnDzwfbz43b3YO9iR4NER+TIEYFx3YdhylTZ2GYarDoSI9AOY8tifkdM4OQaCKEPqCRFxdX/c/sUNPP78LuzLnwAOAL35twFcwbWZO7E8PlFQSqnLIlKuAaGNACjPqtCh4JBzDIMjIDjvB6XUnIicce3vNKxaE8ZhgCCiqpyhI+ww4R5QefsXN9A5WMS+vhNAL4AjAEMEGe6SoxtjJ4Ap17RMu7KkexDkpkqU2mWP/XEWBhGlW5RhwrYyM49OWCGiZ8+AvpYhgsykBzxWrOOgw8R513Ul6LLXQfdnEgYIoiZzbymeE22UYcIOEcUnB/HSnjEAQLH7V3j3yhehPg4R1S8VAUL3F72KR/Nrp53lPXVzzxwML7pBFKZ8oQWD3y7g5ND3cLjwCjrvFfHgkyWsza/EfizusRVhBAp3iOjtmgLwBkMEkSGMDxA6PFxUSo3oy2OwwsQpfXkcwAf2ABURGTe16AZRWPKFFrz43b14buDFxMODl2qtE0EGcq7MzKN9qRPF/YNWb/AQALyBmV8vYHlhPZyDJaK6GB8gAPwE1tKntgvYWPpzzLV06hSAc9g8YIWoKbjDQ+7WTqx9Zk54cGt0+qj9vIr7B9FR6AaGgELXO3jnvWsMEUQJSkOAGAVwTk+FyemuC3vZVK/iGsvYuDIaUdOww8PJA99H32NHkLu1E7d/cSPpw4qcHSJyu3ficOGV8vUMEUTJMTpAOALCIeiKXbpE6A/1KNY8Nk+DKYGoCe0d7MCBZ/pw8sD3sa/zBLYtSybCg60cImCFiEJHP4DX8PGn1zlDgygBRgcIAP363znHwiKvw+rWOAXvil/LerucsxiHHjsxBgCtnSyTS2aqVMK5u6cdx48cxcHeY5kMD7a1+RWsza8g9/wubNMFp4DXADBEEMXN9AAxB1iLi7ius+uM260QTu7L0PuYBDAJAB3FNhXuYRI17oWjT1W87eDAd8rhof36fdyemY/xyMzDqpVEyTM9QHh1R9jjH3KwWhvcrRA5oFyogygVXjj6FA4OfAe9Xbs9by909JfDw0rGw4PNWbWSBaeoktatgu6e9qQPoykZHSB0XfCSiPQ7ljfNwVorvQTgqoi4g0IeG2dpEBnLWcthcOcIOlq6PbfbJnmGBw92rYh8314WnCKKmdEBQvsxrFkVdnGo0/o62wVX3YcRbJz2SWQkr0JQLbcfeG6rVu4xPFSwuWql9VHBEEEULeMDhFLqvIicdSwu8qVzoRKl1Bl9+yisQZezLCJFpqtUCOobQ2s5mK5S1UoWnCKKjvEBAigvRlL37USmsftkramIQMs3D+Dd9kB+VapaOd/zVaD9LC2uMXQQ+ZCKAEFE5IdX1cpfzr4TbCfPsEAVkR8MEETUVNxVK+1WnmBYoIqoFgYIImo6zqqV2/Ingu+AtSWIamKAIKKmZFet7BwsBr7vvr4TQC+AIwBDBJE3BgiihPTmnwZg1XgAYOxqmlGxl/ludLXOWuzpr+5lxSvZvqsHnWCBKqJaGCCIIrJ3sKPibVld18LrJL56YzHyEBHE6o1FrN5YROH4fkdtCRaoInJjgCCKgF2aupKX9owhv3VvJta18PvN3zQLb32yIUSwtgTRRgwQRCF74ehT5dLUlRTXB7Hl5p2mrS5pYmjYvqsn8HEtvPUJckP9G6aFAm9gBgwRRAwQRCFxr2tRXB+suO2Wv2++8NBIaDCtG8OpND2HHB5NC8UQUOh6h7UiKPMYIIhC4LWuxYNPlipu30wDJk1sbQibFSL6kYOztsRrDBGUaQwQRCHo7mnHwYHvoNDRj22Sx5a/v9NUIcFLFoKDU2l6Dtvne5B7fpdVW+IAwIJT5tu6tQXF7ieSPoymxABhsHyhBd097fxwIjKEHZoef34X9ukQwYJTlFUMEIaym8Stkfz8cCKzZK31wcl+7p2DrBVB2cYAYSDnUs/2SH7OQScyh7NWRL5vL2tFUCYxQBjGDg8nD3wffY8dQee94oY56PxwIjKHu1YEMAkA/DulTGCAMMjewQ4ceKavHB5yt3Zi7bOlzXPQWciGyBjOWhEM+5QlDBCGcq+P0Lm/WJ6DzkI2lLR6ijI1M7tWhDvszy99FWg/S4tr/Lum1GCASAEWsiETpS1ERH2s9t+pM+zfvPVZoH3cXP6c00IpNRggUsJdyMbGEEFJsqtHpilIRMkd9gd3Vi4m5mVxdZbTQik1GCBSxFnIhtXwyCRpa42IkjPsd+4oBrpvvnMv0AvgCMAQkQ4ikgPwsr44AABKqXM17nMWwByAvN5+MsjtpmCASBn7Qzr3/C6g6whOshoeGaKR1ghT18Golx32OweDBYh2gLUl0mccwDmlVAkARORDETmrlDrvtbGIjAP4QCl1yb4sIqPOy9VuNwkDRAo5Q4SzpC7AEJGUYvcT6O3ajZ7tA2hdawVwP+lDSkzQIBFHeEiidcSuFREUa0ukzhCAYQD2CX4OwOEq24+5WiimAJxz3L/W7cZggEgpltSNxt7Bjrrq5jtX4HzwyVLTr4PhR61ujWZrdQhLpdoSUZtf+oqfHXVQSj3ruuoQrFaJTUTkkMfVy7ACSM3bTcMAkWLukrrsO23M3sEOHD9yFL1duwPf17kCJ8PDI+7WCIYGf9whwqotES1rxgg/Oxqhxy5crjJmIQ8rEDiVAtxuFAaIlCuHCLAufyPs8HCw9xh6tg/UvoNDR0s3w0MNDA7BOQtUdRS6I3+8RzNGmuuzo3VrCwpduUZ20S0i047Lkx6DHp0DKWer7MvrQJYd+6h6uz3OwhQMEE2AdfkbY4eHl/aMIb91rx7D4F/L7QdY+4zhgcLnrC0RtW3teX52eFtSSg1V20Cf2CcBQEQuishFpdQpj01L0DMrHPIBbjcKA0QTYV3+4F44+hQODnynHB7ar9/HvYDVAx8ADA8UGTtERK29u4OfHQHpVoMx14yLKQATFe6yjM2tDDnACiEiUvX2xo84XAwQTca775TrZ3jJF6ymzd6u3eho6UbrWivuLX3FMEDGKU3PRf4YuaF+tHe2oaOrG71du3Gz63PkC/zcqGEIwLiITDpO8BXznlLqqoi4g0AewGU/t5tmS9IHQOFbeOsTPPhkCcX1QRwuvGLNEPh2AflCS9KHRkTUNJRSl+GoAaGNACi3SIhIv4iMOm6/4Lo8go0tFrVuNwZbIJoU188gIorFJT37AgB2AphydWmMwgoBlwBAKXVGRM7qkNAPYNZZJKrW7SZhgGhi7vUz7NLXrFpJRBQOpdQcHC0OHrefd99eqUql39tNwQDR5JzrZ7BqJRERhYUBIgNYtZKIiMKWugAhIhNKqTOu61KxclkQd9UytmFnaPtzV61kwSnKCq4SShSNVM3C0KuU9XtcN6eUuqSDw4BrBCtpqzcWsfDWJ2i/fh/FdWua5/EjR/HC0aeSPjQiIkqZ1ASICouMAFYRD+cI1SkAZypsS7CmeW75+zvlEHFy6HsMEUREFEiaujCGYIWDcpBI28plJtlUZ38IYMEpIiLyKxUBQndJXIA+zTmkauWyILbJo/LntZZFrpezzr5dKwJ4AzNgiKDmwPEPRNExPkDoWuMlXSfcfbPvlctEZAzAGAC0P9FiZFXGYvcTeG7gRfRsH0DrWivUSvSDG1lwiogAaz2Xtu4ObMvnMbhzBDdvfYbBb5fw7gLXwyBvxgcIAC9XmVXhe+UyvY9JACgM7FCD3y6Ed4QhKHTl8NzAizhceCX2paErFZxiiCDKlpWZeXSi2FSLarW2tKM3/3TSh9GUjA4QeoxDtUVE6lq5LLcjh5ND32v8AEOWRHiwOQtOoesITuqCUwwRRNnSjCGComF0gIDVmjDs6Lo4DKBf1324VO/KZR1b81ZzvWFyt3Zi7bP4w4PN7i92hwiWvqY04vgH/7bv6tlweWVmHu1LnSju37iqL0MEORkdIPRKZ+UwoMcx9LvqhF8QkVHHVM6aK5dtXd+K3K3wijSFZe2zUuJLSTtDBEtfEzU/d3iw2Z9Fxf2D1hJRnKlFLkYHCCcdHs7AaoEYB/BjpVSpnpXL1LqCWrkXw1H7d2/pTuLhwZbV0tct3zzAg6QPgsggzhCxYbo3Z2oRUhQgnIMgPW4LtnJZ6xY8fLIjjMMKzdYnO9DW3YGVmfmkDwXA5tLX6AVwBACuJHlYobJnvQzuHEHnvWKi3UdESajU+uBk/004p3sXut7Bx59ex9Li2qbtu3vaAx/HXwW+B5kgNQEiTPdlFfMtM0kfxib5vr3oRNG8EIGN62fcXP48waMKT1KzXohM4Cc82Oy/jfJ0bwDAO8Az4RzLJXwazo4oVpkMELfvLuHnvzVvva2X9owZGSJWbyyicHw/8n178dKeMSyuziZ9WKHoe+wIwwORT+UQsWG6dzj+SfVx72SobAaIb1bw1vtmNsUf7D2GfX0n8HjnLtz+xY2kD6ds4a1PUDi+H8UnB5Hv3Jv04YRi27Kw26JJcQZGdUFaH5zW5lewNr/yaJA1ZVomA8Tduw8NHQx4xRpn0Avsy5/A48+bGSKC93Caac2ggatEcak3PDjd/sUNdA4WQzgaSrNMBgiT3Vz+HAd7gbtqGdtg3lTThbc+CW1fYXyQEZF/lf7m2oudaOvugHS2BdqfabPZKF4MECkR1YJaSVq9scgQQZSw9mIntu7vxlr7fQAq2J3zrZEcE6UDAwQliiGCwtZsQTssXn9ndniYb5nB4kpzDI6m+DBAEBE1Oa/w0DlYxMMnOzDfMoOZL6fwy9l3EjgySjMGiBRpxm4MgK0QFJ5m/PuIQudgEWt9rVh+MIOf/3YSH83+CjO/Xkj6sChlGCCIiFJq+64etBc7A92nrbsDa32t+JuVP8dHN9/GW+83d4n6tq1t6O3anfRhNCUGiJSxv6k32zcttkJQo5rtb6KW3FA/2nfnsL4j2Mf4Wvv9zIQHihYDREo1Y5BgiCDyJzfUj637u1Fqe1Sx9s76kq/7Lq7M4s2PX8PHn3KFXWoMA0TKNVuQYIigejTL778fdniYb5nBzMJU4Pv/cvYdvPPeNa6mSQ1jgGgSzRQkGCKIvG0ID19O4c3pNwLvY+bXXIqbwsEA0WQqnXjTFiwaDRHO58sw0tzS9rtdr8Lx/eVpl/bMiXevfJH0YVGGMUAYxOqPvI7e/NtAL7BNr4cRRrnYzsEi7i3dQWl6rvEDjUk9IcLrZMIWDUpCbqgfbd0doe3POe2Sgx/JBAwQhrE+FK7g5jOfAweAvq4j2JbPh7Lv1ie7kQOaNkRU+ybKEEFxsrsarPLQ4Vh+cI3hgYzCAGEguyUCeA0nD4S3357tAyjuH0ShuyPURbGixpM/pYkdHlba5nF95f3Q9stpl2QaBghDXZu5g6XFawBeQ2/+6ZD2+jZe2jOG4pODKBzfn7kQwSBCUXPPkLh567PQ9s3wQKZhgDDY8sI63nnvGrp7roeyv2L3EwAm8dKeMeT79qJwfD9WZuZTMwgtjONkiGguJv3ues2QmF/6KrT9MzyYSURyAMb0xcMAppRSk/VuLyKjAEYAXNRXnQIwrpQyru+ZAcJwywvrWF4I54PjGu6UP9DsENGJIgCzPoijxhBBYfMKD5whkRmvKqXO2RdEZFZEUCVE1No+D2AYVsi4CuCHJoYHgAEic+xBmgBwsPcY9vWdwOOduwBkK0QQETVKtyb0u66eAHAOwKYA4Xd7pdRAuEcaDQaIDHLP9Ninp4sCDBFERAHkAQyLSL+jlaCEzSGh3u2NxgCRUc6ZHvZ00VwThwh2WRBRnbpFZNpxedLubtAhoMu1/QiAy1478ru9iIwBWIYVOKp1hySKASLDnDM9Th4AoENEe7EzVbUiKmFoIKKtW9pR6GjoC/6SUmrIz4a6i2IYwLEGtp8GULJbKERkQkRGlVKXgh129BggMs6e6QG8hucG5nC48Apy2Ikc+lMbIhgcKqvUusTXjCgUPwFwTCl1td7tPe47BeBVAMYFiC1JHwAlb3lhHR9/eh03b32G61+/j7t5hbbujlSeVNJ4zCZoxm4rojiJyFkAE37DQ6XtRWTYtWkJwKFwjjJcDBBEGVErJKQtRKTteKl56doNV5VSl/Vldwjwtb2I9AOY0l0bthwAI5uDGSCImtjqjcXyj9/ticg/ffLPA5gWkX4dAkYct/frwFBzez3u4YxSquR4iNMAxmN4KoFxDARRE2okCLDQFpE/uqVgSl+ccNzkHK9gV5a85HP7y7p7AwB2okZlyyQxQFBV9omE30zTIwslv/n7SCbQLQVSY5vzAM4H2H7O3t507MIgX7bv6jH6hAJwACXAEysRxYcBgqhJhB0eGEaIqBrjuzD8rHSm+4vmYHjVrrQz/YRierN7VKJ8X+x9Z/F1JaLqjA8QqLFymYiMA/jArtIlIuOmVu0iClOcgc6kcGZ6kCXKCqO7MGqsXGYbc4WFKQBnoj42oiQlcRINMh00ymMgIjOY3gJRdeUyEfGqzrUMq7Y4hYgf3GYw4X1IqlvDhOdORI8Y3QKhlJpTSnU5wgOwceWyPKzA4FQCBba0uIaby58DAO6qZUhnG8t/jf8AACAASURBVNqLnckeFG1g2gnUhBYJE5jStUMUN6MDhJtj5TK7CyPnsdmyY1vnfcdEZFpEph+sPoz2QFPIXlTrzY9fwwcLf4ZS15do350rh4i0fEg24wnN9BN1HMdm6vNPy98FURRM78Jwc69cVoKeeeHgvgygPDNjEgA6im0qsiNMMTtELNwqAUPA4cIr6NxRBACsza8kfHT+mTTgrx6mniyribJbw9TXI82/Y1mydUs7erYPJH0YTSk1AaLCymXL2NwKkQPKFb8ooOWFdcxgAcAbwBAwuHMExf2D5dtN/TB3CztE1Pu8ax1DWl5Pv0x53aPG8ECUkgDhtXKZUuqyUuqqiLiDQh6PxkhQHZYX1vHuwhewQwR2AsX9g2jr7gBg7oe6W1gns0bXlciatLcAVdOsz4uoHsYHCMfKZZf1qmXAxoGUF1x1H0awcZESqtO7V6wQcXPgM7y0ZwzFJwfRiWKqTorOYw364Z+m52maRkOEia89wwPRRkYHCD8rlymlzojIWd1K0Q9glkWkwjO/9BUwACyuziLfuRftsD5ITfyAr8VvmAjy3ArH9zd0TGm2MjNf9bWqd1yEib9bDA9EmxkdIPysXKa3S8XKZWSOME5SheP78fDJDtxvvx/CEaXP4527ANR+Lf22RpgYHIioMqMDBJGJtu/qQedgEWt9rVh+MAM8SPqIktHR1Y3c87vQXuxEaXqu6rZpDgdsfSDyxgBBFIAzPPzNyp/jo5tvJ31Iient2o3DhVeQw07k0F8zRKQRwwNRZQwQRD5t39WDx5/fhbt5VQ4Pb71/JenDSsyBZz4H8KheSA5oqhDB8EBUHQMEkQ/u8PDmx6/h40+v49rMnaQPLTFLi4+Kjtn1QpolRDA8ENXGAEFUQ26oH+27cyh1fYnrK+8zPGjuomN2vZBCdwcW3vok6cMjoogxQBBV4QwPHyz8GX45+w7eee8alhfWkz40IziLjjnrhRSO748lRHi1FKR5wCZRmjBAEFWQG+rH1v3dKLXNMzzUYBUdA4BJvLRnDPm+vSgc31+zVkRQfroW3NswUBBFgwGCyIMdHlZ0eHhz+g3M/HqB4aGKd698YRUeA8ohohPWYmxhnMTrHZfAVgqiaDBAELnYBaLmW2YwszDF8BCANS7EmplysPcY9vWd8F1wqpqwBzVykGR2bMFWdLR0J30YTWlL0gdAZJLcUD+ksw0rbfNYXJ3FL2ffYXgI6NrMHXw0+yss3JnD8oNrWN+xFZ2Dxbr3x5M9kZnYAkFEpDWy+BpR1jBAEFFm+VkMzMZAQbQRAwQRZVLQMRkMFEQbMUAQUeaEMQvDuY/2YmfD+yNKGwYIIsoMTt8kCg9nYRBRJjA8EIWrrhYIEXla/3dZKXU7tKMhInKxT/z1jDlgaCCKjq8AoQPDjwAcAqAAfOW4LaevuwjgAgMFEUWhVhjYvquHgYEoRjUDhIj8I1gBYUIp9VmV7Q4C+GMR+aVS6mchHiMRUU0MD5QEEekHMApgp1LqXI1tRwGMwPrCDQCnAIwrpeYc25wFMAcgDwBKqckojjsMVQOEDg+TSqmvqm0HAEqpjwB8JCIHReQHSqn/K6yDJCIiMo2IDAPIARjweZc8gGEAYwCuAvihKzyMA/hAKXXJviwio/Zl01QNEEqpfx50h3aQqPuIiIiIUkApdRkAROQwrCDh5z7VwsaYqxVjCsA5AEYGiIZmYYjIEyLyeFgHQ0RElEUicsjj6mVYLRZGCjwLQw+oXFZK3VZKfSUifygit5RSfxH60RER9g52oNj9RKD7vHvli4iOhogaISJjsIKBe4xDXl/vVIrx0AILFCBE5E9hPclDIqJgNau8DmAIAAMEpd7a/Arad/tqiYzF3sEOHD9yFL1duwEAolrQhifQqnZA0AJAYR1rWJMS1nEXEPueXIKcCAC2yFZsk3wju+gWkWnH5ckGBjZOAyjZ4x5EZMIxxsHrg2dZb5dTShkXJoK2QEwppX4KWN0XAF4G8McAfh72gRFlnR0eRp7+Ebas5XDvDqAewhESHPT1W9uA9h0CPCsAfoYZMEQQNWhJKTUUxo6UUlddV00BeBXWl/ESdKuEQ0PJJ2p1l7LWMzN+on+IKEQvHH0K/8kzx3Ak91/j/petG29UFe6kgAdrwPoasFdexn9xcDcud/1LvPPe3zBEEBlARIbtgZdaCVZ9JcBqbXC3QuQAwMTWByDgIEql1E9F5I9EhCvHEEXkd393F0a//T/guR0/gLrfWvsOLgoKSgFFHMZ/vvt/w+8f/V3kCy0RHKl/95buJPr4REnT9SKmdPFFWw5WzQe7dcIdFPIALsNQgQKEHgPxEoDrIvKpiPyJiPwBZ2IQNS5faMHw7+3B2H/8f2JX+2F491UEoICWhx34z37nf8R/eewPsHewI5TjJCJ/RKRfF4+CHvdwxtWacBrAuOPyBXt7bQTARPRHWp9Gx0AMwwoUp/UPEdUhX2jBof278YMD/zvaEG4Dn6AFv1v4b7DlP9oK4P/GtRm2BhCFQU+9PA2rEqVdCOp1x1gHu/KkXcfhsq40CQA7YZ1TywMylVJnROSsDhH9AGZNLSIFND4G4qf6h6ipOEdtd/e0Y3khupNuvtCCF7+7F7/f909CDw82QQv+08IPse15wRv41wwRRCHQQeEqrGJPXrefB3DecXnOebnKfWIhIo83sn5V0EJSl0Xkj+p9MCLayA4PI0//EDvwVKSPJWjBUOd/hRNHfo/dGUQEWDMp6xa0BWIcQL+I/DGAWViDO6YATHMVTqJg9g524MAzfTi5fwxPfO1VhC58LWobXnzqv4U6so5i969YcIooA0Tk59g8f0sA7AZQ97pVQQPEh0qpH+kD2g1rDMSPADwB4B/UexBEWVMODwe+j6I6hDVVeXZm6O48hpGBMQBW12uUIaKjpRvI0AzSjpbu8v/zhRZOnyVTjCul3nZfKSJ/2MhOAwUIpdRPROT3AMzppb0TrwORpqVPKR3Uyj0gH3z6ZH0PBqx9HWN4AKCUgroX0/PLiNL0HHIAOvcXcbjwilWbF2+wkBcZwSs86OsbGsMYeDEtpdRfKKU+B5JfTEuPeJ1TSl3SwWHANQWGyLfVzxdjf8wtqh0qzvQAQCng4Vrd46epgtL0HB58soTcrZ04XHgFJ4e+hxe/uzfxGhxEUQkcIETkaTs06JkYI7pVIgljrikuUwDOJHQslGJ2eFi9EW+IaFnf0XC5h3o8jLEFouWbBwCSCWhxK03PYe2z0qMQceD7DBFkHBH5QRj7Se1iWmlc+pTMlOSJbcv69nj7L2zryZzQVj9fxPanexJ57LiUpuewfb4Hued3AV1HcPIAALyGjz+9zumzZIpQvrakeTGt1C19SuapFR5680/jY1yP7PFFBW4ETKW1+ZXy/+3XvJmDhN2SlXt+F7blTwA6RAAMEWSEUL62pHkxrUBLn+o12McAoLWTzYlUOzz0bB8A4Dn2KDRKkmh+MEOzt0bYIaJzsIh9fVaI6M2/DeAKQ0SM5KGgdY2DhqMQdBbGT3UhqQml1ErNO0Qr0NKnepDlJAB0FNuy+6lNAKqHh3tLd9D6ZDcQQ858uGUt+gfxsuVhMo/rkpkQAStE9OwZ0LcwRFD6pXkxrdQtfUpmqBQevvnNYuzjIdZbvk5kEOWWtgfxP2gFzT64cvXGIhbe+gTt1+8jv3UvXtozhuNHjuKFo9FWHiWqIpRPnaAdsFNKqZeUUnlYAycvwwoUsXdjpHHpU0qe18nqm98s4pvfJHMSe7hlNf5BlAJI2/2YHxRVX+NmDxEAyiGiuD6Il/aM4eTQ9xgiKBFKqVDO2WlfTOuCiIw6pnIavfQpJatSeEiUAG0dwL0YW7MFwJZt9yLbf77QgkJXDoWOfmyTPNTKvXJT/je/WcSOb3l3WTR7dwZghYjcUD+K+wettRjtglO/ZsEpSp9UL6allDoDa22OUV2R0uilTyk5fsLD1zPzsQWKazN38PGn1/HRzbextOUvY3lMm2xfxV8vT+Gj2V9h5tcLsT42UL3FJwstEXbBqeL6IAtOUWJE5A9F5HX9/z8QkX8vIj8Oso/UL6YV59KnlE5+w4NtbX4FW/d3u+8SOmsQ3RXcfOZzjO75p2i/X4j8MSEP8enaW3jzwzeMXUgrCy0RRAYoKaVO63WtLgI4BeAjEfmBUsrXAltBWyA+dIyBeBnWGhQ/0g9OZJyg4cFmV0+Mmt0S8bPf/lOsy92IH03hL795Hf/mw9cTCQ/u1znrYyKIEvaE/ncUwFdKqZ/pNa4+87uDQAHCXkxLRJ5WSn2mlPqJUuplpRRX4iTj+DkJeYUHm3NlxShdm7mDt9/7a/y7v/0xFKIKLgqfr72HP3vvXyXa8hAkRDSD7bvYkkLG+kwv8z0OwFna2vew7oYW0yIyVbWpmjb3yeyL20tVA0WUlhfW8W/+33fx7/7un2Ed4daGUFjH52v/AX/67v+UWO0B5/vhN0SkuRVi+66ecniw/88wQSZRSn2kexS2KKV+phfH/C28izR6qhgg9M7qWiRLRP6gnvsRRalWeEja8sI6/vWV/wf/cua/x118GcpMbYV1/Hr1Tfzpu/9zYuHh3tLmx23WEFErKDBIkKnsxTGVUj/ze5+KAULv7JaI/JGIPO1nZyJyTM/SYC0GSkzQ6Zru8BD3ipxOywvreP+vfoM/+dV/hy/k/wPkIaTOJHFv6xIu3TiLf/Uf/tSIqofu17lZQwRRWukxEL5VnYWhlPoI1qjMfyQiQ7AGTc7CqgJZgtXUYReV6gLwulLqf63nwInCEHTQpAktD27LC+t4d+ELAP8H/uGzf4h9j50A7jyG9fuwWiW8eih1xhBdV2IRf4l/O/MneOe9a4nWF2hda90wquOL20t46vHKY0sq1YlIy8yM1RuLNVsYVm8sor3YGdMREUXH1zROpdQ/BwAROQigH8AArDIoX8IKFed0iwWR0aqNcfi7rxc3ndx6u3YDuBLxUXmzBjv+FH87MIuDvceAbS1oWd+Olocd2LK+DaK2QImCkvtYb7mD9ZZVKFnDwuocfjn7TuLhwY+vZ+bx2GAx6cOITZKtW0RhC7qY1kcAPoroWIga0mjXhYnevfIF5pe+ws1nPg90P5PDg7sVwh0imrUVguGBmk2gACEijydVMIqomka7Lv7u643bqpV7VqecAa7N3MHS4rVA9zEtPLjfi6yGCIqfPFSx1XXJmqCVKP9CV62yK1Be5pROMlG94eGL20t4OtIjq49pgSAM9YaINGLrAzWjoIWkhmCNgZgEsAfAJRH5UtfQ5tRNSkSYo/SrDfCjxrhbeYDa3UheXVBpmZWxemOx/EPUjOopJPWVUuptpdQ/1oFiANa4iOd0kHg89KMkqiDsrgubdLaFcHQUlN9CXmkJEUTNLFCAEJGDIvIDZ10IpVQJwC+VUv8YwGkAY6EeIRE1tXrrQxBRsoKOgRiGNX3zRyLyBKyxELP6up8ppUoiEqgQBVG9omp9cLp5i7/OpkrLgEqiZhW0C2MOwD/TXRdDsALEVwD+F136ehoBFuIgqlejC2WRWVilkih9gg6i/CmAw3o1zq+UUj/VK3J+pQtJnQpSR5soTI3UfPidx3rwO49t/jZrFZKiRu34lvfr65SVrozcUD/ad+ew0jaP61+/n5qiX0RuQbswoJR6u8ptbO+lyIXV+lDphLb01l9j+64e5J7fhcOFV3DzyGcArhixngRtlqaujNxQP7bu70apbR4fLPwZwwOlWuBZGEQmCvvb6Y2fXMHtX9xA7tZO/P6+P8bxI0exd7Aj1MfIgvvt9wNt38ytEHZ4WNHh4c3pNxgeKNUCt0AQJSnowEkgeMlqu6DRjZ9cwS4cRe75XXhpzxh6u3bjze439BoVZBLTWyHs8DDfMoOZhSm8Of0GZn69wPDQJETkLKwxgnkAUEpNVtn2IoAJvf2y8zY9EWEUwAiAi/rqUwDGlVJzERx6QxggKDWSGDBnh4jfeX4fOgrd1tBhMEQE0dYdrOWmmSpUbt/Vg87BIh4+2YH5lhn8/LeT+Gj2V/z9aSIiMg7gA6XUJfuyiIzalz0cglXJ2e0UgEuwQsgwrJIIVwH80MTwALALg1IurNaHpx7v3nTSsjm7Mw4XXsHJoe/hhaNPNXDU2eR+jRtheoXK7bt6sH1XD9qLnXj4ZAdW2uYx8+UUPpr9FWZ+vZD04VG4xlxhYQrAmSrbTyilxPkDa0Xr8j6UUgP6tmeVUlejOvBGsQWCUiHKk4PXSc29xPTyu38NAOWBlRgCCl0cAOelu6cdvfmn0bN9AK1rrViZ+dsNtz/1eHfFYFcpYPhd9jvprgwuoJUtInLI4+plWC0IlWxomRCRs0qp86EeWEwYICi1arU++OE+YXmdpHZ8qwfbn+4pr2lQDhEaQ8Qjewc7cPzIUby0ZwzF9UE8+GQJy+/+9ab3yn7dnUEiaOuESV0ZDA5NrVvXOLJNOsY45OEaxwCgVG1nzu4IERmGVU9pAxEZ0/utOaYiSQwQZLyoWh/8hocNx8IQUZFXePji0vtV7xMkNJjaCsHwYDa1rqBW7jWyiyVdPNFLzuO6ZQAQkZxe6qGaM0qpU67rpgGU7KAhIhM1xlQkhgGCjFYpPPhpfag2/qFWeKj2zdYdIgod/QBey3SIeOHoUzg48J2q4cF+jcOsEGpSKwRlUgm6lcDBfdmTnm3xgft6jzEPUwBehavrwwQcREmZ4jWQL0h4sK3eWCwPrOx77AhOHvg+XvzuXuQLLaEebxq8cPQpnBz6XsXw4A57floRKuFqnWSYZWxuhcgB5YUmqzkDa5bFBrpbw6kEa+aGcRggyFh+Wx/88jNY0is8VGoOd4aIfZ0nyiEiSwWn7PBwuPBKxZYHr9e0kRDhxnUyKCm6tcAdFPLwGNfgYRiu8RMi0g9gSkScoSQHq2aEcdiFQalX6Vupc7S/n772IOHBtnpjEe3FTmxDDvvyJ4ADAPAalhav1Xw8P8LqEomiZWTw24VyeMjd2om1z5ZQmvb+nNvxrR7Ploi0L3jG8Q8E4IJrjMIIrEJRAMqh4JBzDIMjIGwIH0qpORE542q9OA1gPJpDbwwDBDW9SuGh1rdgvwPx1uZXAADbkENf1xGc1CFi4Vb1Fsz5pa9q7rs74vNTsfuJqrcXurzGiFmeG3jRER5KFcODLawQ4XcwJZD8tE5qfkqpMyJyVo9p6Acw6xrwaFeWdI9h2FSJUrusK1sCwE4AU5yFQRRAWN0XfsODu/WhnpPO2mcl5LAT6DqC5wask+nNW97ry91c/rzqydmtVhipxs/j9Oaf3njZYxVSa7DoI32PHUHu1k7c/sWN8sBS4NFr5/UeRtUSwcGUlKRqdRz0bedd15UADFTYfs69vakYICjV6jnxhBke3E3YauUeOncUcbjwCu6sL2FwJ7C4Orthm4U7cxWXCa8UOHp9jev2x0846NlufbZ1tGwMYNvk0YFsW5ZN4cFp+9M9FUMEsDEMRjFDw8ZWCKJoMECQccIePOkUVnhwB4e1+RW0Fztxb+kO2gB0PlnEtnbrZJvv3Iu76lFLpX1ydrJDhvtEblu4U/8Yqkr7dB6HMyjYIaF1rRXQQzBavnlQvt2eU397Zr5ieLBF1Rrh1Y3BVgiieDFAUGaEER78DJqzQ0Q7AOlsAwC07rAe+377fWzburE54a5aRt9jG7/p31nfWMPCK3TUwyso2FrXWstBwQoJ93Fv6Q4A4AEejfUAUDM4uFULEmGqFCLYCkEUPk7jJKMEaX0Is7m71snFXhwpKPvbun1ibl1r3bSN+0QOWCd6d/dBI9z78woP1TjDQyPcr3PUUzyJKDrGtkDoaS5j+uJheIxEDbIGO2VbPcWibGFP1Wtda8X99vsbrtsm+Q3dHDb7pO9ukfDLTwipFR6iFtcUT7ZCEIXL5BaIV5VS5/XPKQDn9AIjAMprsM8ppS7p4DCgp9FQSkXVvF1v10W9rQ4Ayk3/Ts5xBH5bImx2C0KQVolK2zofx3kcG7svvJ9DGLxe73paIioFjDDGyhBRbUYGCN364B75NQHgnONy0DXYKaXCPCH4CQ+NBIcgKoWIakEC2Bgmqv3U2neSLQ9hhYigWJ2SKDxGBghYXRLDuoKXrQQdKupcg50MFlfrQy1RBIdqKwFWOonXChFB1dqfs3UkLn67E+oJEWyFIIqekWMgdCGNLtfVI3hUXzzwGuy6+2MMAFo7s7fgUVo1ciII2nURR6tDyzcPsL5j45+d15gIoPK4iKC8wkOl4OIVdsIaQOnFXSvCazwEEO6YiDDHQqzeWGQ5a8OpBw8j647LOiMDhJvu0hgGcExfFXgNdj1OYhIAOoptKqJDJUOEVaY6CkFDBIC6gkSlVgd3eAir9aHWidTv1M9KIcKt1nvMuhBE0TK1C8PtJwCOOdZJr3sNdjJP2N0XXieWWieSKL5FOr/1uL/Ze520q41J8DM2wr29H5XCQ9BvbH5ev0rbRNmVQUTRia0FQs+QOF1js2Wl1IaBkHqq5oQjPACNrcFOKVFP94Wf8GDyVL5KLRE2ZzCwWyX8hoVagyarjdWoJkj42r6rx7MlIu6uDHZjEDUutgChZ0y4VyOrSoeOq0qpy/rysFLqslLqqojUuwY7Nal6v6HG9eGvVu6VK1MC3l0ZQO0QYQvSIuEVHpIYOBkWv+91XN0YDBGURcZ2YYjIMKxQMC0i/XpGxohjkwuuug8b1mCndGik+8LPSSTO1od6BhtWOomHOcWykX3Vek71nDT9dmVUOvGH1ZURdtdZ0PLeRGlnZIDQgyanYAWCWwBm9U95Wqfu6ugXkVHdzeFeg51SzG/3xWODxfJPPaL+1uhnLEESLQHux3R2X8QxYt3v6x5G60GcUzoZIihLjJyFoccxiI/tUrFmOiXDxLEP7m6Mavx2ZdTah5spXRde4yEqLQGeJuzOoKwwsgWCsiHKE4Wfb65RfMjXWzMhiq4Mv/etd/BkXKIcw5D2sEKUJAYIMk4UTc5JFI2yubsEgp6w6wkRle5Tq/UhSPdFGK+h1z78lrkOIu7KlOzKoCxggKBEJN36kDSvEFHt5N661lr+aWQbP8dhi7ICpVPSzf1shSCqDwMENb24xz40cuL1Mz7BGRKChAZTxj74EUUrRNzYCkHNjgGCjNJoU3NSYx9q8eoaaHTsQZCWBq/w4H78JNcLiOM94QJbROFigKDYxdlknNTMi6hbIbIgzlYIdmMQBWfkNE6iepja+mC7t3QHbd0dG66rNK2zUpXKevhpfSBqVurBemzjebKGLRDUtJKu+5DmD61qxx5V335cMzKIKBwMEBSrqJqK03JSCTIWIoyuDL/7SHL8Q5yqjYNgNwZRMAwQZIyoB7nV6r5oL3aWf6IUV4iodF+Tuy/YCkGUHgwQ1JQa7b4IK0SE0Y0R9aDKrLQ+EFG4GCAoNiZ3X3gFhihbIoJO6wwaIhppfYhi7EajLTtshSAyDwMEpZrXScTrZJPk7IsgJ+RaIcJPkEj7NNAg71U9IYLjIIjCwQBBRkiyyE+1b8Zxt0L44RUQ7HBRLTx4hZO4uy+iHl9CRPFhgKDU8vvts9o3Wj8ntDBOepVaIeqtUOkMDKa3ONTz+vkdTAmwK4MoKQwQVNHewQ4cP3IUB3uPYV/nCbRfv4+Vmfm66gDE1TScdO2HasLqykhSGDUg2ApB1BxYiZI82eHhpT1jyG/di/br97Hw1ieRPFY93RdRDZystm1UhaG8KlSGzW8giav4lZ/Xc/uunk2BZfvTPZ5hdMe3egL9Hn3zm8WKv0Orny8aHUTJPCJyFsAcgDwAKKUmG9k+6P6SwhYI2uSFo0+Vw0NxfTDS8BC1St0XSX0LNqkVIuj4B64uSbSZiIwDmFNKXdIn+gERGa13+6D7SxIDBG3wwtGncHLoe+Xw8OCTpYbDQ9jdF5W+Oab5W2OUgxnDCCJRhId6QxzHQpBhxpRSlxyXpwCcaWD7oPtLDAMEldnhYXDnSDk8lKbnkj4sYyRRXMrUsRBhqfWaBp1+G1aI4HRO8kNEDnlcvQxguJ7tg+4vaQwQBMAa83Bw4DsY3DmC/Na92PL3d2LpDw86/iFI64Np3Re1xDml0uuxKr3fjbQ+RDXLJS2tTWvzK2j55gG2SR6Fjn4UunLo7mlP+rAoPHlYJ3inUgPbB91fojiIkiLFb3Le1uZXfJ84Ky357fe+WeN3QGW1gZRhscNX52AR+/pOAAeA3vzbAK7g2gxLiMfh4b0HjXbBdYvItOPypGNQY85j+2UAEJGcUsp98q+6fR37SxQDBJFBKs3IaCRERCHOyp5eszGAyjMygOCzMryENRujHCJghYiePQP6FoaIlFhSSg1VuK0EPVPCwX05yPZB95codmFQapjQfRFm90fQLqKgrQmVtjdx8SxTu5XCsnpjEQtvfYL26/eR37oXL+0Zw/EjR/HC0aeSPjRqzDI2txrkAKBCa0Gt7YPuL1EMEJQKWRphX+0EH1WXRFz1H+pVKRBWayEw8XfGDhHF9UG8tGcMJ4e+xxCRYkqpq9g8RiEP4HI92wfdX9IYIKiiRqfuxTH+IS2D6cLkJ0SEFTSyXvshit/hhbc+wYNPllBcH8TgzpFyiMgXWkJ/LIrFBVedhhEAE/YFEel33V51ex+3G4MBghLjt4866DfJJFfeDCrIGhlOlQKCWrmX2oGTYdeFAMxshXDqaOkGABS6vMbOURoopc4A6BeRUV1BctZVx2EUjjoOtbb3sT9jcBAlpVbQ1odm62evJyjEOf4h7Ne70mBKoqQppc7XuO28x3V17c8kbIEgo5n+DTIM9bZCRP34laSphQeo/jvkpxUsqq4401dRJaqFAYIikVT9h7Sd3GxxDGIMI5CY+vrWao3KQhAlihsDBBmr2od+kt0XcXaFmDTlMo7wEHZpa6dGQgQLohFtxgBBRE3BT6g0pSXC2eJU6OgHAJa4a2vq4wAAE7xJREFUptThIEoKnZ9va41WCfRiavO6X5XKW1eqThlEpZYME+s/tBc7qx5XMwymtI8/9/wubMtbJa6B1wBcZ3VKSo3UBAgRmdDTW5zXnQUwB13q01GfnFLO1O4LMlu18ta2estch1Xaurw/HSIef34X9uW5TgalTyoChIiMA+j3uO4De36siIyLyKip82WJ/IiiFSKNgydrtUI0C/diW1wng9LE+DEQFdZHB4AxV1iYgqNYByWDsy+iE/aAyihP0FG3+lR7v9M0FgKwQsTKzDzXyaDUSUMLxBCscFAOEhVCxTKA4bgOiurXyPiHtHdfrN5YrBl2giz1XYtJsziCMrEVIuxujPJ+byxi9cYiCsf3o/iktU4GYPXIvnvli9AfL0serj3gLJqIGN0CoeuBX/C4KQ+9RrpD1ZXKRGRMRKZFZPrB6sOwDrEp5AstOPBMH3q7dqOjpRvbliXRE49J3w6j0MgAwDQHAtPs+FaPcb9rKzPzaPnmAfJb9+Jg7zEUunLYO9jYAFqiqBgbIEQkB6BUYQlTr8Lxy477baKUmlRKDSmlhrZuN/Zpxy5faMGL392L5wZexOHCK8jd2om1z0ooTc8F3lcYKb+ewZNp6b4IEhyqffP2GyKqbRfkm72pr2+j3RhE1BiTuzBerjKrogQ988LBfZlq2DvYgQPP9OHkge+j77EjyN3aidu/uJH6KXI207ovnPx0ZZCZourGIEqb2AKE7o44XWOzZaXUGT3Godr658vY3AqRA4AKLRbk4gwP+zpPYNuyxBIeqo1/MK05OUxer2utEFFtLEStWRlhtT4kzcRxEERkiS1A6BkTfqdY5gEMi4h9+TCs5U3PAriklLoqIu6gkEf10EHa3sEOHD9yFAd7j2Ff5wm0X7+P2zPzDYWHRrsvaoWHtHdfVNJISwTHQ1QvKuWnJgQR1c/ILgyl1GU4woCIjAHody1xesFV92EEwESMh5lKdnh4ac8Y8lv3ov36faw0GB5MZHL3RRBhzsig8LAbg8jgQZQ2HR7OwGqBGLcHSeqqlP0iMqpbJmZZRKq2YvcT5f+3rrUCaGxWQBiy2vpgi/P1Z3cAEYXFyBYIJz2Q0nMwpatFghLit5k4ivUvmkW1roykWiHqDTZhH2sja2OwG4MoOsa3QBBlRbUTdhgtB2x9IKIwMUBQoqLovojj23pWxiU0SzdRFNiyQVnHAEGUEo20ILD1gYjCxgCRIflCCwpduXLJ6pZvHjQ8FTDItzC7dLDd6lBv60MzqzXuIM4gkKbWB1alJIofA0RG2CWrTx74fsMlq8PQSNGopLsv0qYZWh9MfV/ZjUFZxgCRAXsHO8rhwVmyOqnwQNWF2QrRDOGBiMxk/DROakxSJasbxWbn6lhgisifh3cfcAp5RNgC0cTsqpN2eGi/fj/U8JBE822a+uUb4ec9qtW6EKT1YfXGovGh0lTsxqCsYgtEk8pKyeqsc4cEu1Wi3q6LRlcJNbVVhAWliMLHFogmlC+0oNj9RHm2RetaK+4t3UlNeKi3+8LUk1e96nm/1uZXGh73kJbfEyJKFgMEpUZWui+c4jiZpyUwmBwQ2bpBWcQAQXXhByaZJosBkyhJDBBkFM6+2CzKFoK0tD4QkXkYIIhSgCf6xkUdTtkqR1nDAEGpUKt52uT+8bCEHSLC3h8XMSPKFk7jbDJ2yernBl7E4cIr6LxXxINPlkKtOhnVNy12X8SHLRpE1Ci2QDQRr/Uuwg4PVJ+wBviFceJneCCiMLAFokk4w4NzvQueLJpPI8Wesv77wIJSlDQROQtgDkAeAJRSk1W2vQhgQm+/7LxNKVUSkVEAIwAu6qtPARhXSsXyrZEtEE3AuVjWvs4TqQwP1bovOP5hs6DvLUtVx4PhhKoRkXEAc0qpSzo4DOgQUMkhAFMAZgHccvwM69vz+v9TAMYBTMQVHgAGiNSLer0LN35AmsMOBbXe6ziCQ5whLouBkZrGmFLqkuPyFIAzVbafUEqJ8wfAOec+lFID+rZnlVJXozpwLwwQKZYvtODAM33o7dqNvseOYNuypKpkNYXHK0yw1SF9Vm8sQq3cQ+taK3q2D+C5gRdR7H4i6cPKum4RmXb8jNWzExE55HH1Mh61Jnhxhg2IyFml1Pl6Hj8KHANBiWP3RbgYGqozfRzEwlufIDfUj+L+QXQUuoEhAHgDM79ewPLCetKHlzrrq/fx9cx8I7tYUkoNhXAoebjGMQAoVbuDsztCRIYBXHZvowPNMnyMqQgbAwQRNY3tu3qaIkCVpueQA9C5v4jDhVeAIaDQ9Q7eee8aQ0R65TyuWwYAEckppaqGCQBnlFKnXNdNAyjZQUNEJkRk1NVNEhkGCCLKHNNbIYBHISK3e6cVIjSGiNQqQbcSOLgve9IDLT9wX+8x5mEKwKtwdX1EhQGCEpWl7otm+XYcprbujvL/7y3dSfBIzFSansP2+R7knt+Fw4VXUOjoB/AaQ4Qh9In9dI3NlpVSZ2C1NrhbIXKANSWzxj7OwJpl4X78YaWUs1ujBGvmRiwYIMg307+xUbo4w4N9uVaIaOvugHS2bbqvU+dg0VcYyQ31b7pubX6l5v2AaMPr2vzKpsGwAJB7fhe25U8ABwDgNXz86XVcm2HoSpLuKvD1bV8pdVVE3EEhD49xDR6GAZxzXiEi/QCmRKTLEUBysGpGxIIBgoiMUS1EtHV3YK2vFYAC8q2e27SuWddvfbJywKhm6/5utHzzoOo26zui/dhs351De7FzUwVZtXKv/Lx780/jY1yP9DgoEhdcYxRGYBWKAlAOBYecYxhExG612BA+lFJzInLG1XpxGh4tFVFhgKDEZKn7opnFsojW7hzu5hWWH1zD4urshtt6tg88utASwoM9HsI+GtDR1Y0cdiKHfpahbzJKqTMiclZ3ffQDmHUNeLQrS7pbNTZVotQu68qWALATwBRnYRBRU6vWBeFmh4e/WflzfHTzbY8tvK5Lr96u3dZCeDuKyAEMEU2mWh0Hfdt513UlAAMVtp9zbx8nBggyTlgLT1E6Obsx2nfnUOr6EtdX3sdHN9/GR7O/wvzSVwkfYbSK3b8ChoDBnSMo7h9EDv7HZhDFiQGCEtHI0t1p7r7I8kyMIK0Obd0dePhkB0pt87j+tRUe3nr/SiYGDV7DHQBvWAWkdgLF/YOBXjuiuDBAkC+cgUH1ssPAw4D3W9Hh4c2POeOAyEQMEER1ai92smm5Bjs8rLRZpYTvrC/5vu/i17MMD0QGMzpA6Okrr8JayhQApp2Vt4Ksq07maGT2BW1ktww10iUUFTs8zLfMYGZhKvD9fzn7DsMDkcGMDRA6PFxUSo3oy2OwwsQpfXkcwAf2FBgRGY+zBjglI83jH8Lm7FZydzElHSg2hIcvp/DL2Xc23L5wq3rhvUJXjuGByHDGBggAP4GjwAaAC9hYsWtMKeWszDUFq1IXAwQ1vVpjUlY/X0wsRNgFn5YfWOHhzWlrJckgZsCVJ4lMZ3KAGAVwTlfmyumuixJQ37rqugVjDABaO8OoNpO85YX18je5O+tL6NxRjORxwhxAmfQ34yxJIkTYLQ931Xy55WF+6SuGAUrM/YcP8MVt/2NvyL8tSR+AF0dAOOS47qKjpGc966pPKqWGlFJDW7cb+bTr8u6VL/Dm9Bv4+W8nMd8yg637u1E4vj/pw6IqGh3nESTQ1Rv+sjrVlIj8M/VMaq9yM6eUmtOtD6/D6tYAaqyrHsPxGeXdK1/grfevlEPEWl8rQwSVcQouEUXB1C6MOWDTWudzsLo1gAbWVW9W1mCzKwCAg73HsK/vBArH92NlZp7fJomIKHSxBYiA66Z7dUfY4x9yaGxd9aZlh4ibz3wOHAD29Z3A4527ANTfJB3n+Ic0TuFMSy2IJAdVElFzii1ABFw3fU5ESiLSrxcLAayAUNIBoZF11ZuaFSKuA3gNOAD0dR1B7vnGQgQ1B4YIIgqTqWMgAODH2Dir4rS+znZBt2rYNqyrnmVLi2tYuFXCwh0re63v2Mr6CUREFCpjA4Re1jSn104/C+BL5zKouqujX0RG9e3uddWJqAFssSKiakwdRAmg+rrpfm4nc7DpfKOkVuUMuxuDLVtE2WVsCwQlj9P/iIioEgYIMkIaZ2AQEWUZAwRFLqwmczaXh4MtS0QUBgYIIiIiCowBgoiIiAIzehYGJSfOZm6Of4gfi0pRVtxff4C/+5rddlFgCwRFKgsnqSyPzUhDGW8iigYDBFFC2PJCRGnGAEFEFbEapXl6tg8kfQhEABggiFInrG4hTudMl5WZeWxbFnS0dONg7zEceKYP+UJL0odFGcZBlJQoNuMT+WO3BuWe3wV0HcHJAwDwGt557xqWF9YTPTbKJgYIikwWBlASxckZIrblTwA6RHz86XVcm7mT6LFR9jBAEGUYp3Omjx0iOgeL2NdnhYje/NsArjBEUKwYIIiIUqYcIqBDRC+AIwBDBMWJAYKIKIVWbyxi9cYiCsf3Y1/fCfTssWdnMERQPDgLgygEcReTYrcDESWNLRBEREQxEZF+AKMAdiqlztXYdhTACICL+qpTAMaVUnOObc4CmAOQBwCl1GQUx+2FLRBEREQxEJFhAIcADADI+bhLHsAwgCkA4wAmXOFhHMCcUuqSDg4DOnTEggGCKKXYjUGULkqpy0qpSwBKAe4zoJQSpdSzSqmrrpvH9P5sUwDOhHGsfjBAUGpw4SYiIouIHPK4ehlWi0U8x6CUiuuxjCEiiwCuJ30cEeoGsJT0QcQsi88ZyObzzuJzBpr7efcppSJpUhORt2C9dvXaBuCu4/Jko+MMdNdDTilVtbVARMb0f5fhGuOgu0MmlFIDju37AcwqpaSR4/Mrk4Moo/pFNYWITCulhpI+jjhl8TkD2XzeWXzOQHafd6OUUseTPoYGTAMo2eMeRGRCREZ1t4XXGIplvV1OKeW7m6Re7MIgIiIykFLqqnPQJKwxDq/q/5egWyUc3JcjlckWCCIiokbpGQ+na2y2XKurosr+h5VSlx1XlWDN4gCs1gZ3K0QOAOJofQAYIJpVbPOADZLF5wxk83ln8TkD2X3extJdCZdqblgHPZ5hSkS6HIEgB6vmA5RSV0XEHRTyAC4jJuzCaEJxFhIxRRafM5DN553F5wxk93lniYj023UcdNfFGVdrwmlY9SBsF1x1H0YATER/pJZMzsIgIiKKm556eRpWJUrAar143a7voKtKjiilRvTlfse2O2HNsJh07dOuRNkPa8BlbEGTAYKIiIgC4xgIohRKsv59HEQkB8CeA38YwJTzOfpZIyCNTF/7gMiJAaKJZPHDJ4snGl2E5gO7hK2IjDvmhjeLV50LDYnIrIg4f1/tNQLGAFwF8MM0v6cOVZ9XRt57SgkGiOaSxQ+fLJ5oxlyr+E0BOIeIRoPHTYfCftfVE7CeYzkcOivwNZMaz6up33tKF87CaDImL7wSthonmrIar0mqmFD/PgZ5AMN6AJmthM3vdaZk5L2nFGELREY06YdP+UTjaFVo9hNNHrpcrUMsRWPiot/LLtfVI3DNb9frBGxaIyDtqjyvpn/vKV0YIJpMlj58MnqiSbz+fdx0S9MwgGOOq6utEZBmRq99QOTEANFcMv3hk5ETTeL17xPwEwDHnN1PHl1R9hoBaX1fAdR8Xll878lgHAPRRExfeCUGnieaKq9JGiVe/z5OetbQhPvEqpcydnKuEZBaNZ5Xpt57Mh9bIAwWdKEW0xde8aPexWmqnWiqvCapY0L9+7jo34Wr9vtnv5e11ghIqzSsfUDkxABhsCALtTTLh089i9Nk7UQDXf/e0Q0Ta/37OOhv4nkA9vsI6PEtSqk5Eam1RkDq+HxeTf/eU3qwlHUTEZExVxGli3AUVhKRCX25XAcCjroQaaRPNP0ALuBRl8wZe658rdckrZKsfx81PZbllsdNl5RSp/Q2NdcISCPT1z4gcmKAaCJZ+/DJ8omGiChpDBBEREQUGGdhEBERUWAMEERERBQYAwQREREFxgBBREREgTFAEBlGzy4hIjIaAwSRQfSqqUMJPO7ZuB+TiNKNAYLILO7S23GZZIggoiAYIIgMoU/giVQFtcsnO8pGExFVxQBBZI7DrpVD43YJwJmaWxERgYtpEcVCD4x8GdZiXlfhWBlVKXVJj33YFB5EZAzAswDOAbCXeh5RSp3Ri4iVLzdyH30cc/o4iIhqYgsEUTxeVkpNKqXOA7gIK0T06x/of2edd9Crik7q28aUUpf0wmf9InLWcTlvdz3Ucx8ionqwBYIoHvZS44cAzOmuivOO2/uxuQXCvjyklBpxXD8E4JTj8iEAyw3cx6kkIv0Jd6UQUQqwBYIoBo4T8jB0mHDZCaDkvo9uJSifzHVXyLJj0GNOb1uq9z5ERPVggCCK1wiAKWDTjIcv4RgX4eAOHO7LL0PP3HAUoKrnPk5eLRNERBswQBBFTETGRGRcn6yH8Kh1YNSx2RyAvMfdy4GjyuXXHYMj672PLceWCSLyg2MgiKJnf6MfAnAMwKiIlABMOra5DGDc4745ANOOy3lsbE2Y0vudc5z467kPEVEgopRK+hiICICIXFRKnaq9ZWSPfwjW4MvJmhsTUeaxC4PIHK+LyHDtzSJzGsCFBB+fiFKEAYLIELo+QyKFnHTrwwfs0iAivxggiMwyWWFwY9SGdIAhIvKFYyCIiIgoMLZAEBERUWAMEERERBQYAwQREREFxgBBREREgTFAENH/v1EwCkbBKCAZAACoJc7uPdmlIQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plan=2\n",
    "fontsize=15\n",
    "figsize=(8,6)\n",
    "plot_divergence(div_2D_nan0[plan],dx[plan],dy[plan],r'Divergence bidimensionnelle pour $z = {} mm$'.format(np.around(-z[plan],decimals=1)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(49, 59, 59)\n",
      "[       nan        nan        nan        nan        nan 0.01602563]\n",
      "5\n"
     ]
    }
   ],
   "source": [
    "print(np.shape(div_2D))\n",
    "c=np.shape(div_2D)[1]\n",
    "print(div_2D_nan[:6,0,0])\n",
    "print(np.sum(np.isnan(div_2D_nan[:6,0,0])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "count=0\n",
    "for i in range(c):\n",
    "    for j in range(c):\n",
    "        if np.sum(np.isnan(div_2D_nan[:,j,i]))>20:\n",
    "            count+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5357655846021259\n"
     ]
    }
   ],
   "source": [
    "print(count/(c*c))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "count=0\n",
    "for i in range(c):\n",
    "    for j in range(c):\n",
    "        if np.sum(np.isnan(v_nan[:,j,i]))>20:\n",
    "            count+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.1373168629704108\n"
     ]
    }
   ],
   "source": [
    "print(count/(c*c))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "count=0\n",
    "for i in range(c):\n",
    "    for j in range(c):\n",
    "        if np.sum(div_2D[:,j,i]==0)>10:\n",
    "            count+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.06434932490663603\n"
     ]
    }
   ],
   "source": [
    "print(count/(c*c))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Moyennage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Mean_V = np.zeros([len(x_bins), len(y_bins)])\n",
    "\n",
    "for i, x_bin in enumerate(x_bins[:-1]):\n",
    "    bin_x = (x > x_bins[i]) & (x <= x_bins[i+1])\n",
    "    for j, y_bin in enumerate(y_bins[:-1]):\n",
    "        bin_xy = (y[bin_x] > y_bins[j]) & (y[bin_x] <= y_bins[j+1])\n",
    "        if (sum(x > 0 for x in bin_xy) > 0) :\n",
    "            Mean_V[i,j]=np.nanmean(V[bin_x][bin_xy])\n",
    "        else:\n",
    "            Mean_V[i,j]=np.nan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Galerkine 1D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def methode_Galerkine(u,v,x,y,z,P,m,h):\n",
    "    \"Applique la méthode Galerkine1D\"\n",
    "    #divergence bidimensionelle\n",
    "    div_2D=np.zeros((P,59,59))\n",
    "    dx=np.zeros((P,59,59))\n",
    "    dy=np.zeros((P,59,59))\n",
    "    for plan in range(P):\n",
    "        div_2D[plan],dx[plan],dy[plan]=fpd.divergence2D_gauss(u[plan],v[plan],x[plan],y[plan],2)\n",
    "    #Matrice pinv\n",
    "    pinv=nppinv(np.transpose(np.array([n*np.pi/h*np.cos(n*np.pi/h*z) for n in range(1,m+1)])))\n",
    "    #Calcul des coefficients\n",
    "    div_2D=np.reshape(div_2D,(P,59*59))\n",
    "    a=np.matmul(pinv,div_2D)\n",
    "    #calcul de vz\n",
    "    sin=np.transpose(np.array([np.sin(m*np.pi*z/h) for m in range(1,m+1)]))\n",
    "    vz=np.matmul(sin,a)\n",
    "    vz=np.reshape(vz,(P,59,59))\n",
    "    return(vz)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def deriv(y,x,axis):\n",
    "    \"\"\"Dérive un tableau 2D y par rapport à x suivant un axe donné\"\"\"\n",
    "    if axis==0:\n",
    "        return((y[1:,:]-y[:-1,:])/(x[1:,:]-x[:-1,:]))#Donc là le tableau est de taille (n-1,n)\n",
    "    if axis==1:\n",
    "        return((y[:,1:]-y[:,:-1])/(x[:,1:]-x[:,:-1]))#Donc là le tableau est de taille (n,n-1)\n",
    "    return(\"deriv ne fonctionne que pour des tableaux 2D\")\n",
    "\n",
    "def divergence2D_gauss(u,v,x,y,sigma):\n",
    "    \"Applique le programme de divergence2D avec filtre gaussien sur la vitesse\"\n",
    "\n",
    "    #Filtrage\n",
    "    u_filtre=gaussian_filter(u,sigma)\n",
    "    v_filtre=gaussian_filter(v,sigma)\n",
    "    #dérivation\n",
    "    du_filtre,dv_filtre=deriv(u_filtre,x,1),deriv(v_filtre,y,0)\n",
    "    dx,dy=abcisse(x,1),abcisse(y,0)\n",
    "    #shape\n",
    "    du_filtre,dv_filtre=good_shape(du_filtre),good_shape(dv_filtre)\n",
    "    dx,dy=good_shape(dx),good_shape(dy)\n",
    "    #somme\n",
    "    div = somme(du_filtre,dv_filtre)\n",
    "    #masque\n",
    "    div_masque=masque(div)\n",
    "    return(div_masque,dx,dy)\n",
    "\n",
    "def methode_Galerkine(u,v,x,y,z,P,m,h):\n",
    "    \"Applique la méthode Galerkine\"\n",
    "    #divergence bidimensionelle\n",
    "    div_2D=np.zeros((P,59,59))\n",
    "    dx=np.zeros((P,59,59))\n",
    "    dy=np.zeros((P,59,59))\n",
    "    for plan in range(P):\n",
    "        div_2D[plan],dx[plan],dy[plan]=fpd.divergence2D_gauss(u[plan],v[plan],x[plan],y[plan],2)\n",
    "    #Matrice pinv\n",
    "    pinv=nppinv(np.transpose(np.array([n*np.pi/h*np.cos(n*np.pi/h*z) for n in range(1,m+1)])))\n",
    "    #Calcul des coefficients\n",
    "    div_2D=np.reshape(div_2D,(P,59*59))\n",
    "    a=np.matmul(pinv,div_2D)\n",
    "    #calcul de vz\n",
    "    sin=np.transpose(np.array([np.sin(m*np.pi*z/h) for m in range(1,m+1)]))\n",
    "    vz=np.matmul(sin,a)\n",
    "    vz=np.reshape(vz,(P,59,59))\n",
    "    return(vz)"
   ]
  }
 ],
 "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}