plot_distorted_test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: No release version set. Setting default to DR15\n",
      "筛选fracDeV，得到late-type星系个数：  2238\n",
      "筛选fracDeV+初步b/a，得到late-type+不是很倾斜的星系个数：  2103\n",
      "不考虑nair catalog时： \n",
      "已从 manga 4675 个星系中找出  1400  个星系 \n",
      "距论文中的 1222 个星系还差   178 个\n",
      "\n",
      "替换nair_ba之前，b/a平均值 =  0.6989220308384543\n",
      "替换nair_ba之后，b/a平均值 =  0.6989220308384543\n",
      "\n",
      "考虑nair catalog时： \n",
      "已从 manga 4675 个星系中找出  1400  个星系 \n",
      "距论文中的 1222 个星系还差   178 个\n",
      "\n",
      "nair与manga重叠匹配到  0  个星系\n"
     ]
    }
   ],
   "source": [
    "# 结合筛选出的星系，做出spaxel分布图\n",
    "#######################################\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.io import fits\n",
    "import copy\n",
    "from os.path import join\n",
    "import math\n",
    "import numpy as np\n",
    "from matplotlib.colors import ListedColormap\n",
    "from marvin.tools import Maps\n",
    "import marvin.utils.plot.map as mapplot\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')    # 不显示warning信息\n",
    "\n",
    "new_ra_dec = np.loadtxt('/root/cosmos/data/nair_ra_dec_ba.txt')\n",
    "\n",
    "# 读取nsa与manga的匹配数组\n",
    "nsa_vs_manga = np.loadtxt(\"/root/.jupyter/cosmos/code_git/Metallicity-deviation/tmp_storage/nsa101_vs_manga.dat\")\n",
    "\n",
    "#read in the vac file\n",
    "vac0 = fits.open('/root/cosmos/galaxies-mzr/data/manga_firefly-v2_4_3-STELLARPOP.fits') \n",
    "basic_manga = copy.copy(vac0[1].data)\n",
    "\n",
    "# load the data of NASA-Sloan catalog v1 0 1\n",
    "nsa101 = fits.open('/root/cosmos/data/nsa_v1_0_1.fits') \n",
    "basic_nsa = nsa101[1].data\n",
    "nsa_ba = copy.copy(basic_nsa['PETRO_BA90'])\n",
    "nsa_mass1 = copy.copy(basic_nsa['ELPETRO_MASS'])    #Stellar mass from K-correction fit (use with caution) for elliptical Petrosian fluxes\n",
    "nsa_mass2 = copy.copy(basic_nsa['SERSIC_MASS'])    #Stellar mass from K-correction fit (use with caution) for Sersic fluxes\n",
    "#ELPETRO_BA\n",
    "# http://skyserver.sdss.org/dr15/en/tools/search/IQS.aspx\n",
    "# 衡量星系的年龄大小的量 fracDeV<0.7\n",
    "fracDeV = np.loadtxt('lalala_manga_fradeV.txt')\n",
    "\n",
    "# https://www.sdss.org/dr12/algorithms/classify/#photo_adaptive\n",
    "# http://skyserver.sdss.org/dr15/en/help/docs/QS_UserGuide.aspx#SAS\n",
    "# 筛选倾斜程度，椭圆的长轴短轴为a b\n",
    "# 计算出每个星系的a b\n",
    "# b/a>=0.3\n",
    "#q_u = np.loadtxt('manga_Q_U_4675.txt')\n",
    "#b_over_a = 2/((q_u[:,0]**2+q_u[:,1]**2)**0.5+1)-1\n",
    "\n",
    "\n",
    "#########################################\n",
    "# 去重\n",
    "tmp_num_ij = nsa_vs_manga    # 临时存储一下，防止改动\n",
    "#tmp_num_ij = tmp_num_ij[np.where((old_same_gal[:,0]==1) & (~np.isnan(nsa_ba)) & (tmp_num_ij[:,0]!=-999))]    # 筛选出未重复的星系编号 & 剔除无效数据nan\n",
    "tmp_num_ij = tmp_num_ij[np.where(~np.isnan(nsa_ba))]    # 剔除nsa_b/a无效数据nan\n",
    "tmp_num_ij = tmp_num_ij[np.where(tmp_num_ij[:,0]!=-999)]    # 将匹配到的有效位置提出，-999为初始值\n",
    "a,last_num = np.unique(tmp_num_ij[:,0],return_index=True)    # 第一列nsa数据去重\n",
    "tmp_num_ij = tmp_num_ij[last_num]\n",
    "#tmp_num_ij = tmp_num_ij[np.where(~np.isnan(nsa_ba))]    # \n",
    "#a,num_last_same_gal = np.unique(last_same_gal[np.where(old_same_gal[:,0]==1),0],return_index=True)    # 第一列sdss数据去重\n",
    "a,last_num = np.unique(tmp_num_ij[:,1],return_index=True)    # 第二列manga数据去重\n",
    "tmp_num_ij = tmp_num_ij[last_num]\n",
    "tmp_num_ij = tmp_num_ij.astype(np.int64)    #最后得到了最终版的nsa_vs_manga数组位置\n",
    "######################################\n",
    "\n",
    "b_over_a = nsa_ba[tmp_num_ij[:,0]]    # 根据数组位置提出nsa_b/a\n",
    "#nsa_mass = nsa_mass1[tmp_num_ij[:,0]]\n",
    "\n",
    "# 星系质量筛选\n",
    "# stellar mass > 10^9 太阳质量\n",
    "#glo_mass_id = basic_manga['PHOTOMETRIC_MASS'] > 10e8\n",
    "manga_mass = basic_manga['PHOTOMETRIC_MASS']\n",
    "manga_mass = manga_mass[tmp_num_ij[:,1]]    # 匹配数组维度\n",
    "#glo_mass_id = glo_mass_id[tmp_num_ij[:,1]]    # 匹配数组维度\n",
    "glo_mass_id = np.where(manga_mass>1e9)\n",
    "\n",
    "# 星系star-forming区域筛选\n",
    "# sfr spaxels >= 20\n",
    "#plfu_sfr = np.loadtxt(\"/root/.jupyter/cosmos/code_git/Metallicity-deviation/plfu_sfr_spaxel.dat\")\n",
    "plfu_sfr = np.loadtxt(\"/root/.jupyter/cosmos/code_git/Metallicity-deviation/tmp_storage/plate_ifu_spaxel_v2.0.dat\")\n",
    "plfu_sfr = plfu_sfr[tmp_num_ij[:,1]]    # 匹配数组维度\n",
    "sfr_id = np.where(plfu_sfr[:,2]>=20)\n",
    "\n",
    "fra_b_over_a = np.zeros([4675,3])\n",
    "fra_b_over_a[:,0] = fracDeV\n",
    "fra_b_over_a = fra_b_over_a[tmp_num_ij[:,1]]\n",
    "fra_b_over_a[:,1] = b_over_a\n",
    "manga_id = tmp_num_ij[:,1]\n",
    "\n",
    "# 开始筛选，在第三列进行判断计数\n",
    "fra_b_over_a[np.where(fra_b_over_a[:,0]<0.7),2] += 1    # fracDeV筛选\n",
    "sel_galaxy = np.where(fra_b_over_a[:,2]==1)\n",
    "print(\"筛选fracDeV，得到late-type星系个数： \",np.shape(sel_galaxy)[1])\n",
    "\n",
    "fra_b_over_a[np.where(fra_b_over_a[:,1]>=0.3),2] += 1    # b/a筛选\n",
    "sel_galaxy = np.where(fra_b_over_a[:,2]==2)\n",
    "print(\"筛选fracDeV+初步b/a，得到late-type+不是很倾斜的星系个数： \",np.shape(sel_galaxy)[1])\n",
    "\n",
    "fra_b_over_a[np.where(glo_mass_id),2] += 1    # 星系总质量筛选\n",
    "fra_b_over_a[sfr_id,2] += 1    # 星系形成区域筛选\n",
    "\n",
    "sel_galaxy = np.where(fra_b_over_a[:,2]==4)    # 四次筛选结果\n",
    "\n",
    "print(\"不考虑nair catalog时：\",\"\\n已从 manga 4675 个星系中找出 \",np.shape(sel_galaxy)[1],\" 个星系\",\"\\n距论文中的 1222 个星系还差  \",np.shape(sel_galaxy)[1]-1222,\"个\")\n",
    "print(\"\\n替换nair_ba之前，b/a平均值 = \",np.mean(fra_b_over_a[:,1]))\n",
    "\n",
    "\n",
    "#########################################################\n",
    "# 第五次筛选：manga与nair匹配的星系中，nair_b/a < manga_b/a\n",
    "\n",
    "# 读取manga vs nair 匹配数组\n",
    "last_same_gal = np.loadtxt(\"/root/.jupyter/cosmos/code_git/Metallicity-deviation/tmp_storage/manga_vs_nair.dat\")\n",
    "# 去重\n",
    "tmp_num_ij = last_same_gal    # 临时存储一下，防止改动\n",
    "tmp_num_ij = tmp_num_ij[np.where(tmp_num_ij[:,0]!=-999)]    # 筛选出未重复的星系编号\n",
    "#tmp_num_ij = tmp_num_ij[np.where(old_same_gal[:,0]==1)]    # 筛选出未重复的星系编号\n",
    "#a,num_last_same_gal = np.unique(last_same_gal[np.where(old_same_gal[:,0]==1),0],return_index=True)    # 第一列sdss数据去重\n",
    "a,last_num = np.unique(tmp_num_ij[:,1],return_index=True)    # 第二列nair数据去重，last_num为最后编号\n",
    "tmp_num_ij = tmp_num_ij[last_num]\n",
    "tmp_num_ij = tmp_num_ij.astype(np.int64)\n",
    "\n",
    "# manga数组编号为tmp_num_ij[:,0]\n",
    "# nair catalog 的b/a为nair_ba\n",
    "nair_ba = new_ra_dec[tmp_num_ij[:,1],2]    #nair catalog 中的b/a\n",
    "# 剔除999999坏数据\n",
    "tmp_num_ij = tmp_num_ij[np.where(nair_ba!=999999)]\n",
    "manga_vs_nair_id = tmp_num_ij    # 最终版的manga_vs_nair\n",
    "nair_ba = nair_ba[np.where(nair_ba!=999999)]    # nair catalog 的 b/a\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "##################################################################################\n",
    "# nair_b/a中较小的值\n",
    "manga_vs_nair_id = manga_vs_nair_id[np.where(nair_ba<=0.1)]    # 数组位置进行同步筛选\n",
    "nair_ba = nair_ba[np.where(nair_ba<=0.1)]\n",
    "\n",
    "manga_ba = fra_b_over_a[:,1]    # 原始b/a，数据来源于nsa_v101\n",
    "#manga_id前面已赋值4561个星系，即所在的数组位置信息\n",
    "k=0\n",
    "# for 循环将manga_ba的值替换成nair中b/a较小的值\n",
    "for i in range(0,manga_ba.size):    # 循环nsa\n",
    "    for j in range(0,nair_ba.size):    # 循环nair\n",
    "        if (manga_id[i]==manga_vs_nair_id[j,0])&(manga_ba[i] > nair_ba[j]):\n",
    "        #if (manga_id[i]==manga_vs_nair_id[j,0]):\n",
    "            #manga_ba[i] = nair_ba[j]\n",
    "            k+=1\n",
    "            fra_b_over_a[i,2]+=66\n",
    "        #if manga_id[i]==manga_vs_nair_id[j,1]:\n",
    "            #fra_b_over_a[i,2]+=1\n",
    "\n",
    "fra_b_over_a[:,1] = manga_ba    # 替换完成\n",
    "fra_b_over_a[np.where(fra_b_over_a[:,1]>=0.3),2] += 1    # 加入nair的b/a之后的筛选\n",
    "sel_galaxy = np.where(fra_b_over_a[:,2]==5)    # 五次筛选结果\n",
    "\n",
    "# 筛选出的星系的manga数组位置\n",
    "num_sel_gal = manga_id[sel_galaxy]    # 筛选出的星系的manga数组位置\n",
    "\n",
    "print(\"替换nair_ba之后，b/a平均值 = \",np.mean(fra_b_over_a[:,1]))\n",
    "print(\"\\n考虑nair catalog时：\",\"\\n已从 manga 4675 个星系中找出 \",np.shape(sel_galaxy)[1],\" 个星系\",\"\\n距论文中的 1222 个星系还差  \",np.shape(sel_galaxy)[1]-1222,\"个\")\n",
    "print(\"\\nnair与manga重叠匹配到 \",k,\" 个星系\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# num_sel_gal：筛选得到的接近1222个星系的数组位置\n",
    "# spaxel_data：获取的所有的spaxel的数据\n",
    "# 0: spaxel所在的星系的数组位置\n",
    "# 1: spaxel在2800数组中的位置\n",
    "# 2: 星系中的y位置\n",
    "# 3: 星系中的x位置\n",
    "# 4: spaxel's surface mass density\n",
    "# 5: nii 6584\n",
    "# 6: oiii 5007\n",
    "# 7: oii 3729\n",
    "# 8: oii 3726\n",
    "# 9: log(Oiii/Oii)\n",
    "# 10: log(nii/Oii)\n",
    "# 11: E(B-V)\n",
    "\n",
    "# 载入数据\n",
    "spaxel_data = np.loadtxt(\"/root/.jupyter/cosmos/code_git/Metallicity-deviation/tmp_storage/galid_ij_den_NO_logOiiioii_logniioii_star_forming_v2.0.dat\")\n",
    "# spaxel_data = galid_den_NO\n",
    "# spaxel_ha_hb = np.loadtxt(\"/root/.jupyter/cosmos/code_git/Metallicity-deviation/tmp_storage/manga_sf_spaxel_flux.dat\")\n",
    "# spaxel_ha_hb = galid_den_NO[:,[5,6]]\n",
    "\n",
    "lala = np.where(spaxel_data[:,0]!=-9999)\n",
    "spaxel_data = spaxel_data[lala]    # 剔除无效数据【spaxel数据为空】\n",
    "\n",
    "sel_spaxel = np.zeros([spaxel_data[:,0].size,1])    # 计数筛选数组\n",
    "for i in range(0,num_sel_gal.size):\n",
    "    id_sel = (num_sel_gal[i]==spaxel_data[:,0])\n",
    "    #id_sel = (lalala[:,0]==num_sel_gal[i])\n",
    "    sel_spaxel[id_sel] = +1\n",
    "\n",
    "# 按照筛选出的1222个星系挑选数据\n",
    "spaxel_data = spaxel_data[np.where(sel_spaxel>=1)[0]]\n",
    "# spaxel_ha_hb = spaxel_ha_hb[np.where(sel_spaxel>=1)[0]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAETCAYAAAD+spv+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsnXm4VlX1xz9LkCugV0EQQSVAMqcU4yqmPwmHVDQ0TZscIicsslGTHEJTS9PMLC0ckJw1h4KcyBG0RC+JpqKGQIoggyigIAis3x/77LPXe+95J7jD+3L393ne593vPufsc8477b3W+q7vElUlIiIiIiKiKbBRa19ARERERMSGgzipREREREQ0GeKkEhERERHRZIiTSkREREREkyFOKhERERERTYY4qURERERENBnipBIRERER0WSIk0pERERERJMhTioRERERbRgispOI7G1edxSRX4rIX0XkjHLHi5NKRERERNvGtcAw8/oK4AfAJsBlInJWOYPFSSUiIiKibWNX4F8AIrIxcDzwQ1U9FDgHOKmcweKkEhEREdG20RlYmrT3Tl7fl7z+N/CpcgaLk0pERERE28ZM3GQCcBTwgqq+l7zuBiwrZ7D2TXhhFY1u3bppnz59WvsyIiIiKhxTp05dpKrd12eM/iK6vMR958EjiauptfBb4I8iciywB/Bts20I8FI5g7WZSaVPnz7U19e39mVERERUOETkf+s7xnJgRIn7XuCsgVaDqt4oIv8F9gRGqepjZvNi4Kpyxmszk0pERERES0Gorj9XVZ0ETMrov6DcsarpviMiIiKqAhsBHVv7IgpARHYuZ39VfbXUfeOkEhEREdHEEGDj1r6IwngZKKXsryT7tSt14DipRERERDQxqsD9tX9zDVzh9x0RERFRfah0S0VVn2quseOkEhEREdHEqAJLpRFEZChQB2wHXKyqb4nIYGCGqs4tdZxqu++IiIiIikelWyoWItIDGA8MBGYDfYE/AW/hclY+Br5T6nhxUomIiIhoYlQ6+6sBfg9sCuyIm1RWmW2PAqPLGSxOKhERERFNjGqyVIBDgW+p6gwRacjymgNsU85gcVKJiIiIaAZU2Z/rmjz93YAV5QwUBSUjIiIimhjeUinlUQGYDJzRwErxOSwnAY+XM1iVTaYRERERlY+WZn+JyE64wlrdgMdU9Y9lHH428DQuIfJ+3IRyqojsiqu1sneBYxshWioRERERTQwfqC/lkQ8iMlZEFojIyw36DxWR10VkhoiMAlDV6ap6OvBVHC24ZKjqyzjmVz0wHOcKOxp4Gxikqm+UM16cVCIiIiKaGE3k/hqHC6KHcZ2L6hpgKLAz8A2v4yUiR+AsjscoE6r6pqqeoKq9VLWDqm6tqsep6n/LHStOKhERERFNDO/+KuWRD4ly8OIG3XvhkhFnquoq4E7gyGT/8aq6D3BcWdfqXGeFth9ezngxplIFEJmf81q1RytdSURERCkoRim+G/hLeFlOPZVtcG4pjznAIBEZgnNZ1QAPljEewL9F5HzgN6qaikyKyGbA1cCJREHJiIiIiNZDsUD9N5MHwA6wqMyhG0JV9UngyTLGsTgPuAg4UkS+paozReSLwA242xhWzmBxUmlFiFyY0Xt6i19HRERE06IZkx/n4LS5PLYFStblyoKq/kZEHgT+DLwkIv8AjgDuAM5Q1ffLGa/iYipZzIYG24eLyEIRmZY8TmmN62xpqPZIHxEREZUNoTjrazJwvtt98zKGfh74tIj0FZEOwNdxul3rBVWdDozCubmOBF4AvlvuhAIVNqkUYjY0wF2qOiB53NCiFxkRERFRBAJs3L7w4/D28DvnK1qSOYbIHcC/gM+IyBwROVlVVwPfAx4BpgN3q+or63WtIjUiciUwEZgAfAUX53k5cYOVhUpzf6XMBgAR8cyGkktZVjqyXV4REREbEkSgfan/rquzu1X1G3n6H6T8YHwhTMNNIser6p0AIvIo8DvgYRG5TlWrVqU4k9mQsd9XEp3/N4AfqerbGftUBVRH5040J+S6t/TmFr6giIiI9YYIbFwyX6rVMQPYX1Xf9R2qugw4SUTuA8ZQxdL3mcyGBq8nAHeo6koROR0XXDogczCR04DTAHr37t2U19mikLNyX+vlrXMdERERpaEUS+Xvn8ADnwDlxVSaHKqal92lqn9P5FpKRqVNKkWZDar6nnl5PXBZvsFU9TrgOoC6urqGk1OrQLVIaYJbcnNSuDkG5iMiqg0isHFN4X2OqoGjgBsXZMdUKgXlBusrbVJJmQ3AOzhmwzftDiLSU1XnJS+PwAWrIpoB1i2XNRlK99DWhS1xRRERVYIqqycsIl8DTgV2ADZpuF1Vtyp1rIq6bVVdLSKe2dAOGKuqr4jIL4B6VR0PfD/RuFmNkzAY3moX3EQoZr1Ed1dERJWhhEllwnKY4CqVtKr7S0S+CYzFaY0dkLQ3wi3aPwDKiuyKycrfoFFXV6f19fVNOmaxlfz6jx9cYeuTn9KQcVbqteZjqvnjsyyVYsc0xXVFRDQnRGSqqpal9NsQdTWi9duWeL6ZrPf51gci8gJwD3Ap8AlQp6r/TmRa/gHco6pXlDpeReWpRFQWVEcX/KPXheERERFhIDhfSymP1sengWdUdQ1O9r4WUgbYZbi8mJJRUe6viMbikaUd0zorfhmcNCaXnnsjJzbPtUREVBSqK6ayBCdECS6WvRNBR0yALcsZrHpuuwLREn/eTSHLsr7Xue7HF9ExO2F0zMOJ2DAhhL/pykc9sBsulj0e+LmIrAZWAT8HppQzWJxUykQ5VkGzx1y65+nPiGusXxwld3JYl4kuvdahZR8aEVF9KCVQ/4F70MqBeuBXwKeS9s+T9rU459zzJLl+pSJOKhWGdbJMTjATxi2tLQPTeALKmfxi2k1EW0AJk8qwbu5x/aLWzVNR1WeBZ5P2BzgJ/BqgRlWXljtenFQiIiIimgOVEYQvCyIiOB2wRaq6cl3GiJNKBSPLfZYGx5vpPE0+9kwlW30nImIDRnUF6hGRw3DFugbirny1iEwFLlHVB8oZq4puuzJQTmyk2QP5vrL0bNBJrim3rH/8ozH+1OB1cl++cvXho9Gzk/NbV9cUTZ+1n5tYGk1gl8fclIgNEFU0qYjICFwM5THgB8ACYCtceeLxIvJdVR1T6nhVctsR6Z/xftl/wvmoyP64zOTDl9bxWjwteJeM6+uW5/rykAqaGnJf7ms9umXOGxGRgxLYXxPmw4QFQOsH6s8BsuTt/yQifwLOxSkVl4Q4qVQwUpeXXd1PvrDRBJEzoezXI2ffkmAngkW5k5Nqj9zxrTT/qIzxF5m+7TPG7zY6JktGbPgoJVC/jXtc/3arC0puCdyXZ9u9wPHlDBYnlQpGfuvD93u3VOl17XMmqBeLu55aysKIiNigUEXuL+AJ4As4SZaG+AIwqZzBque22zIKWBJue37rpGBcZ/cMSyOP+ypFHspyplX1bDLW3g2vufl4xdHdFVER8DIt1YGrgRtEZEvgr4SYylG4zLJTbFl3VS1YiTdOKmWiIftKy5rDWwfrIv1S9jmiRRMREVBdlsojyfOI5NGQsvlw8izJtoLTZfXcdhtENnMr9In/2BeVYZ1wejpuOtnYOEwfs+vI5Pkds71/7ti6W8MJxbji7ilwGRERGzKqS6Zl/6YcLE4qLYyyZF7WsQhWeo4T8jHF/DWUHovJe64oEBkR0RhVZKmo6lNNOV6V3HbloJLcXZlxjKIwOScvJhaIrZRwQUZ9nVeMJXyg6fc6XlPsMVtl5qTEWikRbQqlaH/Nhgn/A1q/SNdWQGdVnZW8FlwVyJ2Bx1R1QjnjxXoqEREREc2B9oUfw/rDdW6R1tqU4nHAj8zrC3HJkIcC94vI8HIGi5ZKHqyrGnGxfdf7eh4xYx/S2uKRlQcnCxPgraaIiBZFdbG/PgdcByAiGwHfAc5R1V+L++P5IW7iKQlxUmlhlCXzUqREb4CJjdig+y0N81mAS835d0+232iOedL8Cfs4zqLQtdPCF3LO/Cp7IPKLBhfuzhFdXhFtFlUUU8G5395L2gOBrsBtyevHgZ+UM1j13HYrY33iAxtyLfYsunKMpUS0eVQX+2sOLn4yGafo95qqvpNs2xz4uJzB4qSSB/bPsJil0KJ/nHPsi8RCOTOXeqyXu2cn20OOy+yqg0ek7R++4uR8Nj/+3bRv1cfhl7BiehfXOCuMPf3EPcKL1BKyLLKG4pMti+juiqgIVJelMhb4tYgchJtUfma27Q1ML2ew6rntNoh1laRvDnn8dYW3ZLJybhpaOU2jqBwRUQGooklFVX8lIu8AewJn4CYZj67ADeWMJ6oZFNINEHV1dVpfX9/al1EWipf2TawCK62yk9k8OfnTftbGTMx2nwczJHwHduo3LW1PH5tYJduaY241ba9S/IrpuyU7ez9OKhHVAhGZqqp16zNG3bai9SOL7wcg57De56skVMlc2rxY14REKC8psSnOX+3Ivdfc5EvpTgNts2xXWtb702np+2l7eW2X9bjCiIgmQhH214RXYYJzLLW29H2TIk4qlQxfO6VP6Drw5lCEbTbHAPDmFHOMYbzfffCpAHx17Pi077mzP5u236Q/ACctDdZuap0AHY9xf9QdNglVRZdsvXU4gU+aHGLOf0uee8mCjQWNK+O4iIhKRynS97u5x/XPtXqeSpMiTioViGIFuaoRzVmyOCKi4lBd7K8mRYypVAByXGoZ4pByVu7+Kbsr0d3a6PKP0m1rd+0cdkzkUw7uF1QWJp51RNrueL6zRF6r3THtOzQVJIXZS/sAsOKe4E46+KRg9UycOcw1phnG1b/MhV6RMZHcmGeiHJc89zF9t5Qec+lvAjszbEnKiIgy0SQxlU+J1p9T4vlOjzGViPVATnD6zOKBaRkMDMrog9w/4IhWRbZIZ25MSHV0m4qftWlUEfurqZH3tkWk3zqO+baqfrKOx7ZNDLe5Jabfu7+6mb4HTHu2e1r7Uae0a9jCv6TtCYOPBWAiwTrZaVLIiH99/g4A3FN7TNpnYyrbn+RW/8NOCuWp+zArba/p5yKRj3FYuKb/mhtIinR13CkE0VdYxruVxvcT5xXr9qcbrZOIikKFy7SIlKcvrqo3l7pvobl0Bq4gS7nYE/j3OhwXERERsWGg8i2VcQ1e+/96yegDaJJJBeAS4M0Sx2oHXF/qifNBRA4FfpeMd4OqXtpgew3uBgfi9Gq+pqqzSx+/5SREyjlXDGRnw7v6tp00I+17O2GtFcOVhESBn8hWmftklQ8IfYVjOuV+ZvEzbkMQYJPWvoiC2My0dwTuBm4E7iOUE/4KcBLw1XIGLjap/F1VnytlIBFpR5mZl3nGuAb4Ik6Q5HkRGd+gJvLJwPuq2l9Evg5cBnxtfc7barhiPqo9cv9sXjSTz/f8ftmHb7vQ/dHOeSn8yU6YGVxZem1S12RuWHBswQdpe0aPTwPQ9/B5ad+jD4SCKQuSKpP3EMa07q8VOLfbOf3CNV86/Oy0vX+PJwGYbYI/AwY9mLbvvea4cDNfTp6rnvH2JzMptfKlRLQeKtz9paopu0dEfgNco6pXml0WA5eIyMfAlcAXSh270KSyH7m50sUuco2I7Ae8VuoxGdgLmKGqMwFE5E7gSMBOKkcCFyTte4A/iIhoC9HY1jXQugfPJq2+BceLKA7rDS7d05tnrExBzGxVADmRvIoBGzoiwaBMVL77y2Iv4Fd5tr0MXFTOYHlvW1WfKWegdT2mAbYB3jav59CI+xT2UdXVIrIE2JIcgXYHETkNOA2gd+/eSe/6l9AtFwP0ENOexQvsDeRzuYS+zVc68cdBHYKxOPGlEHS/kh8D0Gu3uWnf4PmhMuiF/RKrYW24lmcWH5C2j+96o2uY4P+FhD+LR3DXPY5vp32zzaS4ig4ADOGJtK9djzU0RP2qgWn71Q47p+2tbn4rbS8Y25tGODw050iixJ2nRHIWGrq8GrOvzHfBj2snjW7G3TWUDJyeusSyFgfxj7eNo3omlbeBbwOPZGw7mQYytsVQabed5TBoaIGUso/rVL2OpPhMXV1dSZZMa8m2y2UtdqoNCqlV4evITM5v+ZXHdynz/GXsFzXO2gAq3P3VAOcAd4rIy8B4QkzlCFy8pazwQiFK8dh82zKgqnpyOSfOgznAdub1tsDcPPvMEZH2ON2cxaWeoGV/0G4lPE3ynTdZKV9qtl06Gk0MDJ8UOXGnI8jCV6e7RMTFCwOl2NKLfzH/fAC+uVv4KHcxxLzrcTIud5nvzM8IvAgff9necDW2INCDt0uMytGEIl19TcxldfKrOqRDWAC9aYLs3tIBWOBFK3/Y8C4T3JtM8H81fUVcUTtpeN9Syf4TRme6zDLzTGwi6uxiC4yWt4BbEtHqKhOl1Kj/J0xwXvFW1f5S1XtFZBAwCvgGsDXwLvA88C1VnVrOeIVu+2BKpxQ3VTzjeeDTItIXeAf4OvDNBvuMB76Fy90+Bni8ueMpheIehaTds/bLTYjL/iNKz9ct/pCbA1bBoCkEQRuNX1RduroRi7CVgBJkWobt7x7XP9D62l+q+m/KZHnlQ6GYyrb5tjUXkhjJ93C+vXbAWFV9RVy92npVHY+jvd0iIjNwFsrXm/YaivxI7B/9omKuDz+B5PtDSbaPyo6p8FjybEUajXjkwQudpfJF/pH2/bBf8KMtS1iDrxLiGDZm4lldlhF2uanIdS4XA7DG2PF9fcalgbVO9mNy2n4wSYocRpCJ6UVgmh3Eo2n7koPPBeAUQyA8eeJtaTv16g4BPck1ZXaYyGdO6mWuxxm3Oa4pb4DlYdKFzyjPpO8je5PN5hNs6eZihIs/xT/gtoTqCtQDICJdgF1xnqCHVPV9EdkEWKWqawsfHVDI/XVSGdejqnpTGfsXGuhB4MEGfT837Y+BY5viXB5NWddjnQtrleibj4BL+VHS+m3m9h54AkAbVfSLaH1U0aSSpHL8ChgJdMR5nvYE3gfuBeqBkldEeQUlRaTkmQk3qVR0WKqQoOT6JLk1XH3m2zdrlZrpx89iHFkFkmPD57VR5+VAbhxl823DvXTr8B4AZ3B15jX9h90AOIIgEvkSQRrfWyXW0vmALdL2NRf9pNHl22rW47c7GIB9+Gfa5xllAPdzVNruhLuXDqxK+26aPzxtj+rhLLDNWJZ5LzdwChCsM4AFgw2jbHKW+7H0zy/NnzGMNEY1/t7kfpcKnyuiMtEkgpK7iNbfXuL5BrSuoKSIXIZjyf4YeAKYCdSp6r8TBu3pqvq5UscrNJf2LbCt6jB16ielxz+6F9y8Xoh5KZWFVJyzrCrc2bCFwiIiqoj9dSIwSlVvSqwWizeBsnQgC8VU/rcOF1f96JYx4XQbnQZ0n6HhgqLh6jO/Sm0jpLkRZqIZGsYbeLNL+5l61r5h+8LAqF47zcncb3t0kDA5zCSdHJHEMn7A79K+QwwV/YZZLmX/mr6BuLevsSreT6ySL3N/2rfnr19O28+ePwDItV4sPDvsbUPoO9DEURaxZdo+Dresm094/5f12DRtfxvnXR3JNWnfAhpLr1im2oLpxlLxn+uUn4e+4RnsdBszG276vdBln2xro+OSDo07X4yWSZtFFbm/gC3IL8fVgTKnx5JuW0Q2Bj5HoPu+Dfw7qhFHREREZKC6inS9jFMqeTRj21DKFAguOKmIyEbAubjsgS0IiYcKLBGRq4BLVLVxGnUFo1Q3WO4xbqX6dEYIqlDd9YgID/89aY3YSpRZaWFUl6VyMXCviHQE/oL7fx8gIkcBI4DsRLk8KMT+2giXavYl4CFcfshs3Nv1KYIG114ickQ5lLPWwUKcOyrrT9+4qRZlF1Py+D8ZVvoPMivPxPY9lDyfafpMPZGpu7j+4ZeH6xt3VuPrH3R04Bl/aALVr/MZAGa8vGvaN2rXC9L2QX3/DsD2BPdZH0MZ7sECAPa8Mbi8PvrBRmm7HauB4OYC6LPSCE7WOAJBt7c/TPtWBo8X/TsFi9u7wt5k+3BfBHma8cn3emcjA2fPeyauIM3xGBry30Mzrd3yF+PyskF337b1amzM3QfqbcLlfmFRsuJLvmW+Nwda+nmMpbUpVNGkoqp/E5FvAr/GqRKDEwd+BzhBVbPkW/Ki0G2PwCVAfjnJD2mIMSJyJHBXsu8fyzlxJWBdabzrIrnRsCRwREQx6yFaF1WMKppUAFT1buBuEdkBVxZwMfD6uiSWF7rtbwO/zzOh+Av5m4j8Idm3wieV7jgrJSsh8fSik0O+GhuNkWH15JlQfrjQZeRdJStD573mj+Mr7k9l3ChzUMaaYS4h8c/KqHj67c93DcWyL156Qdq+pvY0IJeG+wID0vagJNPylZOD9bDLpGBdfGbfNwDo9FGgAbe/I1xXzTcSC2WW6TM1Jg5/7LHwIrmFNQPDV/KIeRPT9oieVyXXFKyXf7JP2p6Gk2GxRIQPBgUCwZtLEl52yBOFnUx7XPJsZGI6HmQqVl7UxTXOM98To9yPzzmdnEfm5czWnRDihNTyqOwkiwAR+TmudtVcVX0DeMNs6wmcqqq/yDtAA2xUYNtOwMMljPEQuT/PqoR0D4/M7XJhpksspaTmG/e+Jri4iA0Keb9Lefojqg+6EazapLRHBWA0TmcxC70oI/ERClsqSrYicENURymi9hvDFj2w74+nCWdOJFZi3VJ+9ysj295bHeNMn5F2SS0UG2e5yuzrV7efD11nHHx52n4viUPcPjNI0+tDYZ3wykhnYVhK8Y9rf522fSzlM2FhwrFfDZIq79ztxu+x9L1wAT8JzdobnIXy7u5BD2/rXkHGqL3P/TCH58QsQj0w1LGj2eWZYAn9Y9//S9t1OE27k14OGWVzd+2Ztm9K5Pmt1WXlZzYf8q4b5+CgjffYRBNUSSyUHOvk0S5p+5zLHRX5l4PNgs0YWmlY6kbz/XjSxsqKab81pp9H66J6oQKr2xVas1u0ejhayK/fuC1QVgJWoUllOnAI2TQzi6HkFtGqeojMz9V1iohoQZQq2JifdZidHxUnqZaDirCmfalBlVXFd2liiMi3cMK84CaUP4rI0ga7bQJ8FphIGSh01+OA34rIU6r696wdROQI4LuQVIuqZKyem/q4G/24FmX8OG+ZHwowWUFHW68jtTDMD/peMxklMZGcErl9zPYeSb9hfDHd7uueNv/Su2nX7webAE2yurbJjzuPDJTyVyc6ZYUBB7+Q9u3BtLR93JX3AnDbj7+S9t12d2gfd6Pb/sDJwaQ4/NNmeZ68L1vXGpFV6+47wT2tNhZJ+/+Y7QtCU5L2RweE1d0Xv/l02GHvpB2UXfj+mt+n7XntXFDGJk/aRMkXOriYy1RCwTCLrY52emELZprKCwPC4u2XY52Fsu2k8F4fZ5hml73rPjf7WS052RrxGX/0OQsXJ8kf3V8bDta0q+igynKCD0GAJTQuIbIKF964tpyBC00qY3B04r+JyAPABEj5pn1w3OXDcHGXMeWctLXRlsvCRlQuSp1Qmkp8tCnZZZGplgtFctS9mxsi8mUcMX4rXL35gtaFqv4Fl5OCiNwE/EJVZxU6puRrKcQYS4pgnQf8AFdIxu8swFLg6uRiVjfFxTQnRHZXmBhWh3ZSecSsGENKBXq0Pzb/j7iRiKCNuQxPnp8MXX+8OcQ/vnOWkx7Z9vKw+p0zNhSxIon5bH/2K5nnfnNmIvQ4LayIdZvQXlrnpENmt+uT9k1mv7Q9cmFSTviwMOZrz38qbdckZnnfUUGuPidL2BsoNthoLAl8ZWPDCLNxFA42bW94m/FXfyG02/83aRhNUDVa1T7g+WqnwBnZ49dB0GvWT138ZQwj0j4r6TItYb1ZQU2b8zMgsfDmG2mYR+YHcczNujimXU2H4MpYMNHIxFycPFuNMbsU+0ppLi9ouiqScVLJRlMISu5e114fqd+0+I5AT1mSeb6kUOKXgAWquqvpPxT4HU4+5QZVvdRs6wJcUU7RRBGpBTqr6ryMbT2BZar6YeMjs1HQ6ZdMFheIyC+BgTiZFsHJtNSrass7A1sQmSrCDfdpRvHJiLaJcgq6hWMa/6nHSaP1oAir1l+nZRzwByCtVZoIPl4DfBFXZeh5ERmvqj6ufV6yvRzcgFsinpqx7QKcQVFy3aqSIknJ5PGv5FHd8LGSm7NXeBdzJgDny2aZ23PYPUnZ3xzfuB12XPLcJ3R9R8yLxFCYc5axTr5ojj/Z/bnsfHbgQVgZem+hnHZ0YHd1Whrcosu/1xWAu/4YygV/k8CemtXdrd7XPB/M9F5rwmKldlayZrDWh4mDpIr5nTP6gBXJZXU0jDGjTZlbRNWPu2Poam+r6vgiWcZ6EZv/klhNKw8IP+R//DSwx744z8VkDuwZeCcLzIflM/WtpWKz96ckF2AFO3v1CO/VjTO/C8BF/ULM65d7n5u2VxyeMMmGm3vKXMi2nMxPnBiaD03h/lLVSSLSp0H3XsAMVZ0JICJ3AkeKyHRcKbqHkiqO5WAw+b94D1JmDmIhmZbBONHIks2e5JipqvpRORcRERERsaGh0KRy53Ufctd16d9ktzKG3QaMPpGzVgYBZwAHAZuLSH9VLSKRnoPNISlo1BgfA13ybMtEIUvlCVyGxHMF9kmRmGVP4CqGlTtTtihE5nOjBsLayd0Di+eiZqhZHhHRXMgX3C8U9M+KyeQrLhatmXWDIqwuMKkcc9rmHHOaM9V3kv8tKmPorLxAVdWrIU81vuL4Ly7InxXcP4z8sviZKDSpCHCYiOxYYB+LUjN9WgmFBCUDzveaKibgvv2kECh/87LGlRnvuDn4h75xmfHvTG/wDLn04gzo5uE7I484XkRHs4hYtMooMia5fbZu/MzaUFvtmj+6WN2PTNndowzn17t3xtxqtEmsjZlc9zNXBRruvr8OyYOpOoy5JAxluKMPPVqionWlPWtOlbi6dvqB2W5Z897FZvkptqh0UnFy73mBMr3SuNdWJ8d3YkXat58pON/75cT/FrQ3c6pQDkySL5fTkSz8rd+hAFxK0NTpUzs7ba86ey4Ab040ZTytD8An237Z9I3I/13RhaWxxZpqUoiTS3lw7q9mEf+aA6ZAkUtOnLueY/4e+JOIrMI57eecdjidAAAgAElEQVQBPXF5LCOB75QzWLG7/nmR7REREa2AUmnFTcUUiygPLlCfUbjN4KkJy5g04UPIjS4Ww/PAp0WkL05F+OvAN9f1OgFU9XoR6QH8jNycw4+B81T1+nLGK1Sj/lOZG4pjbiUW75LN6pQ96kOg9GyzcVEZIpE5go8ZLB1riXhFlb3z1C1PVqcd/2CkQaYb9+U77skm1FnK6pZJ7tL0KXukfWcMCjIuXvr+kYeOTPvGDw08Xp8oeOrnTXanKTL5/hVuVd5ldFjdY/IceS15/q/ps1poPhBtb99aGobK3CjtqiF8/N1ob9qgvb3uFFYCz9+2IW4+u10Qz+yUWIOfXRwka6Z1NfTkhc5se7j7kLTPVrw8KBGeuItAiviFWZPtwOsAvJF8JgALLguU442GOxNx7dZXZNxIQyTWtqWvW2t4UTZrMbi6Sqcvt0U0BaV4p7qOenN9/+I7AnvJy/koxXcAQ3Axl/nAaFW9UUQOwwk6tQPGquol63Ot5nybA/sAXXGJkf9S1SWFj2qMQpbKPsDfVXXZul1iREREJaFwnCV7W6QhryvW3/2lqt/I0/8gjpXVpEgmkIeK7lgEhe76j8BNIvIUjgg6XlXX13fXahj4GaifZNwBtha9pQn7QP0Qc/BZWdYJ8Gxy3N5mRTgoNPGFm+y5HgtjdezjLJQVi0yN9ydDc8DZLuhwmPn+HGSk2Paf4hjeBw0KKjperh5CzOCIoXenfeNv/Wo4gU8ktFEzYwl0WZhYKDa19RDT9haKZbDfYNrekjBWyGyTvNjH0pP9Os2KT1qOi7dU7LW8Zdre2LIWi6U6e5if4t4LQvwljfWY+JAvHQAwqfte7pqNqWX/NCYlSaWTTHLpH42l8FzyxZhnyhRwdriBBWMTqyWfkKmVzvdfwTwFw4L8fn4CUNv6g295lEIpfnrCBzwzYQmU5/5qFojIbrgqv3W4OM3nVfXfInIJ8LSqljzZFAqud8f91GYB5wNvi8gUEfmZiFSd1P3UqZ80mbxFRES1w8nsz08f6z9eGKtQCYm2hDW0K/j4/LAtOfO6fpDrVG5xiMhQYCqwNS7RcmOzeSWOrlwy8loqSVzkoeRxuoh8HldC+ETgYhGZgSs3fL+qPptvnEqEag+n/+UxzrS9YKRZGQ6YFG5v2mV7h32fTJ7NKvGpy4Op8oUrkkxAa6kERRY2282thC1LaIuzg1z7v2YOAWBZv5CIaQtyHdTHVZz6trmBkav+kLYHdXBscGvd5MikJAvZRU+HLLxu1xhK0inJcyCUYfQowRvn1uL4bON955uYS58TzHYbB0gslHeeCl2bGfmXWs/Ksp5nK15pDIAU1q72yyD7TbXlhP+SPBujs/+MOaE9M2l/2hxjrLotd3Ss0A8MpX8FndL2ZxNanGWU2UTWB49xAaYVu9iUgPAdvOjyM9P2Fatce8nIrcOuOYXgClOCi7PGmigBs9votLxEW0NLa3+tJ34FjFPVUxNpLvulmUaZX4iSnX6q6jPqRyU046Nwk8yZIrIA+Juqtlw68HqgB28BvYvuFxFR7UgthkXF6ccFx2lgzZTLJGuOmE0lQxFWrr9MS0thR8CvWhoyt5biAvclY50iSar6Gm52+1UiOPZlnGpx5aLHxnBCD7a6PPFjm9/Ed2++Mm1f273xl3zaxL0b9QFpwSn9rcktOct8Jv5jCmEOzjk6MIJuSEwB66d/aMrRafvWQccAcPxZ96R9218e8pBe6rE7AMemy2z4R4eg83JRwj6y0iw2ZvHC0275vsfEYDL4wl4Au8xKzmXjGHebtrc6rPFu4xiJdua/QzVjVhui2TC76k/m+G2sJfWxaXurwFgnS00eS20/f4I81+KtpXkZfRBEMe2yyApe+riQsaTUyM948ck9CGUGLLz8S51VxDTYotZZqK8OCtbLm4PCZ3H+RMMK6598x8LXIkcKKE2/6ZbnD9v3L2q6QHqkK+eiFEtlyoSFPDdhEbR+TGUB0C/Ptl3IjV4WxXonLKrqPFX9o6oOXd+xIiIiIjYE+Eml0KNu2NZ897pdoZVjKsCdwC9E5P9Mn4rIDrjki9uyD8tGIe2vSet2fZysqv8tvltEREQlIXVReUumAB1H5MJGllBanjvmwQAUlGmpMJwP7Iyzw31S3N9wgfuJwC/LGaxQ8uNa4AVyxTIKjgXsB+y5DiqZzQ6RXgoj+I26qPJPBudRh04C9Yv1srRrBsENsddZRofEy8DZuutG8WSjfZOEtvnBD6MbGemexOiVc81nMDI0nx7kotL1plrhXYa/+895BwDwVs9Q46PX0hA1b+/Zw8eHMc88/qK0PQp3j93GmOC8MXRXJiK7NVZl2CYveu3Sg0yfdTl5V5uVccm3JkuOm22C932yguImID/FBN0HeSqxvT7rCvNrMEs6sPDntQmVa0zbv8X2f8KSEhJSgK0Bs6g2eDVmJyeuMdH9GYTkOE8F7/12+Py+v92v07ZNtPTHvWm+lwu6mxihnwwsEcKuNY+jMawrLM+koskys9Hk47dvIJNKUyQ/bl+3hV5aP7j4jsBXZcJ6n68pICIH4qg83XAO38dU9R/ljlMspvIdVS1VULI961FsWUS6AnfhvMOzga+q6vsZ+60hKEy9paqVHcuJiNhAIEX+I7OC8YEo0DjQX4iFVo2TkUWVsb8AUNXHgMeK7lgEhSaVG8klixbDmuSYdSURjsLNjJeKyKjk9dkZ+61Q1QEZ/QXRfWA7jq3fjJ94Q2E/s3Fy+MJ//NEFAGwiQRpFbzDUTSNSnUprHBuooxP7hV/ewYOTiOl54ZjtDg6ewTmXuRXnnJvDoNtODIKl1+JqdPx5aSjidkjtI2l7ZM/fAPCrNT9L+2bUBnWdux52Vs1vVwXz6YPHe6btOw9w6oVvjwj6dCeEekBsPTExKyyJxf5O9kyerXXSt/G+0w0NeadTzHZbW8UG5T3suInQ5XxrnVgasQ+0G+tm8TOh3dWXu59pjrHSMN7SsZTpl03bW3u2yqWlJyfWWHvj+N26azDLtt7UDfxA98BEsFaLr8h53Mv3pn2dtgtCojYB1su/LLgvWCebvxO+r0tqki+5pbLfaq41iwnWQOhUJzVgfO2XEYgvxijLRxRoA3Dsr8LaXy9MmMu0CXOhFQL1ItKp+F4BqppPGr8RCuWpZFUBK3RSJbtyWKk4kpDH/mdcFkjWpBIREdHMiMmL64dSVIp3G9ab3Yb15qnrZ7ZGoP5DGtOHC6Fks6tZtJnXET18jWRVnSciW+XZbxMRqcd5zC9V1b/mG1BETgNOA6CmN9cO/jFe4+Kbk8am+80yy+tN5AIAdjKKNHKZee+PDe2145Kl9JCw+eALgpz6Hya5Zfn3XgraJYcQLI2bhjurYZsxgec7b0TwnW99kfuuTTp/r7RvPMPS9haJ9v177ULQYseD/pe2t3vU6Zx88GywTmxM4shN3Vt364chI3Hrl833O7FTP7oskAQ7n742bHchHeYHRjM9DGWYhB2dIxZvrIf5ht7cI4lF9LFrNpu8mNxiznrZfHsXJ1Tfjsaq6mqslqU+1pInQtgxMQQ2tjawjen44+31WUvq8eTZipJbaf7kwg9fHbwLr/Tcnoa4behX0nZP8wb8xwRw9ucJAHodHbbfe58JlHgLKlRsyE0aHZlYEPavzHrOxyXP1tKZbK2SDEHLyYXpySINJWNG5+xTiox/taHC3V8nUd6kUjIKsb/KklNW1duL7SMij+IYBQ1xbkZfPvRW1bki0g94XET+o6qZRWRU9TrgOkhUiiMiIloEfpIoFBvZECcSj0qPqajquOYau5ClcmuBbQ2hQNFJRVUPyrdNROaLSM/ESulJnniOF7VU1Zki8iSwB6VUJvsYx4Z5xK2ubpegVvhDDaf6VyLPMn1mmIM6fifwBZaa5NKa4Y41dVSPYCz1HxR0WLzvWxcFxteFxqO3Rw/nZz9uxI1p3xn8Plzz+W55fMSqv6Vdz3QIVHKfPLnI0Kt+8WhIrrxpe1dbZ9GbQYblEaMIuehDd9wPjrwunNMk/71yvFtJ7/KQeXsNu+m2JLnxOFOX3ic8AqmMSx8TL1k8LrRr7UrfxypsnMOIc/47sSQ+Z5hkn5iCYnMSptduxpJ4x7CffKhjiDn+HWMp1XoLw67e7bV4Vpy9ZhvT8e+LIQfm/IL8vZhMgF1mhPf13f7uwvuYgmu+iBrAHg+Fm3lgqIvLtDP0tM8eHYqDvZrQtt4bFGJ19asCg9DL90y8z3BcvhO+71+93MXV7pZwLTklH/6VPF9RjmZYcbGNag/OW1T6pJIFEdkZGIizt8eq6rsi0h+YX45afaFJpaVFI8fjKo1dmjz/reEOItIFWK6qK0WkGy68+uuG+0VERES0JkqRaXllwmxenTAbWjmjXkQ2BcYCxwCf4OaFh3E5K7/ELanOzDtAw/EK5KmUJY6lqmWl8mecb0ucCEhv3E0cq6qLRaQOOF1VTxGRfYAxwFqcGsBVqnpj3kHt+L5Il2di2VXklxvv/5WjA7F/GsHR/ubMIKOh81ysQcZkSLMAJ+/mcmFGcWnaZ9u/w9XO/TG/Sfvuem142j5sR8cEevDU4Gcfe33wSnr2kP3yLiCEorx4oZXD/+yasJSe2s6tXrczGvM9VwY//XM1bnm9/7n/Svus5Mz7jyZFvG41RbxsbsidyXNYJDPdpFH1Mav+jokBtdhYPV1NVG1xht0625zLr46MXifbmPagJD6ywsRpOtqfcjLWCmP9rDCCkV29gWqPscIW/lpsTMZGBb2FZNOCrQyM/7U9HbosWb6+a6hz/F5CQWxn3uy/mjrNWyVGvs1jsbhnqZP/GVYb3mwrhDlllYvhLfmt8VTb+IwPpeQzVG5pnKeyvtphLYmmyFPZtq6H/qD+a8V3BH4qv2/VPBURuQ5XMu8EXNTzY6Aukb4fDpypqrsWGCIHhSyV2TQTOyALqvoeuRq6vr+eRC9XVf9JbspZREREBSOfrL7vz6pG6fpHV32Rrypyfx0N/EBVnxCRhhf9P6CsKsCFJpXTaCZ2QKugM86v/Q/Qy0Gs7p9dPSYr2inGoT9nbMh8PvCkkD7/Wj/3Xh+4b+iz0vQrEt6TlTv/NjelbW8BWYn023YMVsmda1xQ4q3rw5LXChL6c1kJdcsY6pW09104lSys7t74S9/5jsDu2v8pZ6GsMOyujoZS0eVTzkJZbeII7X9nBvPJ/zZPJSjzc+H3QvtHScykq3W9W6vFWwgmjrGZsY29VfGOid8MMpbC7OQzNjYVS41VMshm/ftzWm1WH5aylootKOaNhnz2umeP2TweG7/xIRPDOBNzfXsuCUkzS3u7/Idl7UJJhB4mBOnLSNuYjP2D61Drvo9vG6qatVa36+DaT549JFz+/D5pO1WIsG/mO6b9kBeszFdEbMOva6RINcm0dCT322yxGbnaEkVRKE/lhnzbIiIiIiLyo5Q8lQrC87g6WQ9nbDsG+Gc5gxWiFL9N6ZaKqmpZJlJr4npOBJM5HhERUT2oFm2xKnJ/nQc8mqR8/AX3v3+YiPwIN6mUJmKWoNBU+hgbkPtrs22XsNflD/A17nIdJolxYL8wEQ9M3Eu3LA2lISeeFDRdvjjFRFITt84/v7FP2rVlbbAivRviLENQqzGusPlJxHPilBCR/e6ga9N27ctu31d3Dz6fvXcLvqTXX3LjH2X0TvqOMQVDktDaon0DpdjSj33yZP9xocKhdS9pwh/4900h0r7vD40rLXFVtbfqazZQ7ZVFjBtqhbF/R99k9vV1Smy9lqDpyVLHjqbWBNo3NlSSjRP301DjsppvEi29p2YnSwPOqmFvxs+p15Ic94lxlS427qse/loeD33saNrebWaJDFbyxbvtrBtuQfa+tYvd96LTpy0nOjT3SRaWlrQxyxRc8ZVA7+GYtM/W9PErbJuoO79HcF+92cMRACyleepZ+6Zt/p48n5UvGG8TId2kUMmTw7pAEVYVkWmpFKjq04mY5KXAH3DiwBfi0mgPUtXnyxmvkPtr+HpcZ0RERESbRSkxlTcnvMabE16H1i/Shao+A+wnIh2BLsAH5eh9WeSlFG9o2KRuF922/m7efGkXAE7Y7fp02y2XBcmyO8521MxeeSQyJhklysFJib2RhwdW84kPmETCBHaV+MTSwCO9ptYpyFjfqzWZvSCdPf9dVw9P21d+3wlObkbISzr1/4xOSkLffWWwqeb4mklk9Fw7U03xowFGkqU+CdoHQyc3uc8Hki0d/07TPix5Psr0/ca07U/JG1tGpp8xpu0JjYbmOz0spNnJK81YQUizen8nsVq2MZbKdGOVbJNYArXmmCxYynFHSyn2sWcr7WItjSyCvrWU/HvYL6MPcgP83mrpmtEHqVTMR73CZ/lmTSCbeBn95YYgYqX130vMJUs6yWovNwI8T08JFUd5MnkeZQLyVvLFBPDXhVrc3O6vpqAU96jbVr9Z/73iOwJXyc8qQvreIskJ/BQwXVVXFtvfomgkKRn8FJzSk6eLvI1zj92YJU8fERER0dZRLTEVcbN0jaqOSl4fgEs+7wS8KyIHq+orhcawKDipJIPfiRN8n4/LXRGc5/wQ4GwR+YaqProO99KiWK3tWbRqSzbq4Zaat4wN1sm2Z4eUuW9c5pbM55wd5E6sdMZ/2C1teyrwyMOCpWK/SC+wBwCvXv25tO+27wfKsKd0XnHl+eFCzer9zpNdVuaopWF5//D3h6Tt2YkF9CN+m/a98nSwSnzy2xETJ4ZBrfTI8OTZSNAvGBRWjn03cfGZ93cPK9IPdg8r2r43JvEbmzlk4hgpSdHUdffSLUCutLy3ZmxMwqZbeQvEfGN3sms7fw4rrW+o0Nt4S8FYKvb4pclbNN9YFz2s1ZFYAh2tdWCp6N6as78oG7/x+CijD4LVYmnG1lKxVo9fN9prsQv+JG7TeXGgh+/W+420/VZ3F2t5LyO+BsHytZbM2zlKmQ5PzwyqS9saeaI5ryRWUQO5fE1rya5f4mM1xF+qTKblOOBi8/o3uDTcC4FLgF8BJdetKsT+2gGXKfAicHjDYI2I7AVcBfxVROpU9bVSTxoREdH2UC2sraZAleWp9CKpNCQi2wG7AyNU9TkRuRK4qdDBDVHIUjkPeA3YX1UbVXRMTrg/jsN8Do7nXLHYSNayWYcP+WwPFxSYekxYpo4klBa+9EdO8NEmP95vggJ7ELImvSVyysiQ0TckdSjDhQmz5eLv23q8AafgqFCv/LixdQFwzFInVNk+5Dty6H/D+H1GOMaO9YfvMS4ID94+3Em6HNHPWCpmxfvaJY4FvuOtQS6/70WGPZb48bt0DVluXRabjDe/Ercrdpu86CXY85UTtuwoT1TKkLsHgmCjLZFr4wj2OA+bzuXjR/ZazTG1SaJmrbn9nIJe3mqwvxirjuctkCLlknP+Z2xMxbPCrCWTVbgMQiqavZbOedoZ6D0v+RL0zN4+N9lgi0zZGGMqSjojCKXOmRZiNqlmarEiXhswHPursPZXBWEZwUdyAPC+qfj7MVBWQa9Ck8oBwKisCcUjEXa8CoygVUREREQRlCKNX82oMvfXU8AoEVmLUy+0Yr47gJFbKAGFJpXuYIIJ+TGL3DVlReKTFTXMeak/p+zmrIMetWHJ/lt+lLa/1sHlsTxnLBVrnViT1lsals9vS8Q+wf4A7G+sFxuf8UKSae4MMNVUUzriwMTCsCtiIz3v4zt2FWn3HeHpU2bFftvxIaazX8JeyyngZPNMfKzFMo+sO9zHKX5q+uwq2TtEjzZ9NuZiV+X+OHuvljF1R/JsrQcjI5/eg13d2/iKl3yxjDMr6OjjO/b+bU6JP87KsFjBSP8eWUvI/rr8ojWf4IXvt/9D1hKz1+2/YlbtJF9JOw/7vvXypwwnszETLy75IUEGxua0LJiYfDC2sJcRGg3srtMz2V0bch0ViyqaVH4E3IKLn08jt77VicCkrIPyodCkshDoQ45uaib6AouK7BMRERHRZlBKTGXuhBeYO2EatHKeiqq+Q1rHtREOIb8jNhOFpO//jOP17J3PBSYiNbiSPf9R1W+Vc+KWhny2Tvnb82zU2eXzrH0mLKn/fnRgsfjVhS1mtaVZ6v+TkD3/6NtfAmC1+Uq0N+GTZ653megDVgZLZ2pNMAsGT0jclja3wlgF9x0/FIA6QhZ76g+HkF/y1dD1Qvew1O+Iu9cdzwwxE2sdfHSVy2PoPMWUCLaCkH7VbvJBUosBSENR9rdjGV3J8e9fGNhjXV40MRlbetiv8C3HxMYnfMK2ZZJZJpvP3v+56bNxFn8uY+nljLV38mwtEcuu8tn1lhFmLa0Pk2f792CdBv5ztSyuLEvIHm9d8ptk9Ofb7o0OYyR8tGPIWVlV4w70ig/usoKpcxdOst0KlS4zVsubg12uV64wZMNywQ7V6N5qijyVLnX99MD6i0ra9145vuLyVNYHGxXYdgnOn/aYiHyu4UYRGYjLVdmBNM0uIiIiIsLLtJTy2NBQSKblDRH5Mm5t+ryIvEuIsfTB1Zr/ADhaVV9v5uuMiIiIqBpUGaW4SVFUpkVEugKnEjLqhZBRf0NSXKvisXNdR72tvi83JUXUr34mRJev3Pe7aftaXHvGlJB595YpzDGGEWHMpF6FdzMBHL78obS9dydXj+SHJjnRBtX3X+rCVe2NG2XpVmHlUjsx8ToaN8w1g05O275OyzBTeGTre4LPaNYxjhpqE9u6jDbuJx9A/zB0YbONPOXXMJJT6RUAr0JxlemzrjwP6yY6xLSt5It3D5lEzBz6h79G636yZWK8R9jUYMkhIHj3Vb7fuXe/ZdGEIbvahL2+rOWZdb/5z9AKXmS5v6z32hIV7LV4t1c+cczkHlaa65vVKYiIe4LHFPZK+25YFVgNNR3c9yoNyINLj/PIQxWuRldXFprC/bV5XX/du/7KkvadKEduUO6vojItqroYpxd7WbF9IyIiIiKqjlLcpKiaKjLri/fpwl18jatPdRbKFteELLdehLavnCjdQ/D684YSbOuCe0qxTY58oNPQtO1rhXcyJfK+eFAg02ki835nzy+nfTZ4OmKoC35uPStYHyOnBEmYGYO2ddufMRFtIx3fd6vkviyd1CiUp0UkbRKiCZR7OY/em9uMydBc+qyzqmqvMTyOrMqGNjidL3nRWx3WUrDX7ddxJhE0R8bFn8MG4i08Zdeu6K0l4S2wWabPWg3eQrL3Ysfyv6QM6i4QxC2tpWLbWdaHDdrb99XDGBLv7hh29vJBM0wi7ZMJvR2CIORdS0MN9RWLQgItFyRJjcHoBkAXumcR07eBWCdNjTipZEBEHge+W6r8iohsBDyKS+//b7H9IyIiIjZkVHJMRUQWUka9LFUtlgmVopClMgQMj7A4ZB2OaTF8wBaM5wjqr3dL3kE8l247wjjiR05wlsDFwwI3+Ptrfp+2Z7frk7Zvw8mg2MSx19khHMfVSV+wPpY+EmImC9q5z8muaEa/HLyMmqxafU1ygLf7hnP5Gvdb7hsc/l1eNjGTZHU7Y/C2add7EgpyDUj89Is6Bef7NlPCWL2fTZbldvVuVue+iBgDzXYrc+JX3fnkVCxl139l7Uo9KznSYknGdttnYx5+/Pcy+iw2z9POUtywlGb/sdiYj73mpRl9W2Xs2z7PdnMt/udd3zWYavY7+CZOMmU8w9I+Wx7hoZmJZf0XY3JYDdpbCmt0ReukONayUaXLtFxDMxVhLOb++quIlKOl3zaKs0REtBHky36PE0txVLL7S1UvaK6xC00qf17HMSsyu36X+dOpv/JzHPFjF3SwjKicUqvDXC1UK/s9rd0eaftsI3N2HLcnz7elfWMTdhnAhCRAYYX5au8I8Yfbj3dJl9aSOabfvWm7Jlmp1y4Ix+zyYCiytexkZxTaAkpddgrb39nOLdX7nxGsk54fh9SkmmkubrTNuWb5bn33nt1lF1zWkvBVku1vx0qy+IREa+nYIly2338TLSMqix1lYQtaeXaYlUHJinnYcTbJaNvt9vz+Lconbe+l7+3xCzL2te9lltVi+t7vG5JG5xuzZVoiZGqLt9nkRB9L+c+qsH3Vx+bEwxMLJSd50RTMutlNGG1FTqU5EGMqGVDVb+fbFhERERGRH0rFx1TuLr5XgKp+tfheydhtpZywbFmnHF7Pr2/+PpArw9LOLG8fWXgkAJ1qgsP8lNob0rbNMxl1a5KgYXIjXro7xFS8vMujqe46zDamwBa4opnfXxxKEF/W9Ydp+7trrgWgdrFhV5nV81s93erVSt9/ZnkoxlSzpPExOTIkf2nwDCx9K1hVmy1x513e2ZQYXmAkXbyIoJVQt6WHvfT93qbP5qbYJU2WupCNqfgx7Orfyqx49pYd3zLJ/HH2nLZ0sHfyWuvJWiJZ0vRZcvMr82z3/y/GuvpoO2M1rnTv6we14Q205RfmG80Vn18ymcFpXyeTKzV1SkLx+5I5v30vJpdfzreh1bIhu7+aIk+lY93O2rf+toL7LJvwFB9OmMQH198/Q1U/XXDnJoaIPFHO/qq6f/G9HNoMpTgiIiKipVCK+6vTsAPoNOwAPrj+/nxVeJoN5UwS5aLNTCrt+66i681vpdnndnVvP/wru7uM+uVjAjXptREhG3nHCUGc8dnjHftq2vFBL/0gQmVlv6K0LBCr9eNjMjd3DZblIkNZerSds3DWdA8f07EPTUjbvY91y+/evwvL97cGBt9778tc/4wLA/tr+84hviLewWlWsbVzjVWUZNJ37mqsExvH8KwvSyC3jKiDM7bb0sM2Z8UT161Wqo15+JiGXf1b9tjqjO3WQvMWRpaEPATWWOOqubmwMSd7r/68+Rhjyfum5n9mWk2I1a2pcRusdWLL/f5pVVBy8PGRFdt3CYONMefyu9rM98nBuhCx4o+lWRwbsmXSHFAkJ5baltBmJpWIiIiIloKzVKrn71VENqdbMwsAACAASURBVAOOxAkEb9Jwu6r+tNFBeVA9dx0RERFRRagW9peIbI8rn9cJZ3MvxPkB2gPv4+z4OKk0xKZ8yH5MTgP0XmIFcl1Ovl7KESMCOWL8mOCeOnfE+WnbC0me9+Jv0r5ndg+ZgG8kVOEzVl6d9j1WE2q3bLbc8WA/0ykE120VSX/89gSacA58gN0EuXs/biLZSaB5u+XB5fV21+Ae22JzR6uubRdcXh9tZYLyX0jcXva3catpe8kXu66xAWEfSLfCiFbGxLqvRibPlkiQVTnRBsoPNG3vHrPnsu4tf9vWJWbH3zyjz5IHajO22yqRyfUt7ZXt8nivnfuOWemUNwyV3AudepFQgGlLwwlW/NG4uh5Inq17KxT0BE5v8JyL6MpqflQZpfi3OAGkY3G/sMOAF4GvAb9KnktGSZOKiJxYYPNa3Ez2gqrOKbBfsXMcC1yA+1vaS1Xr8+x3KK6UVDucSvKlWftFREREtBYUYc3aqplU9sIV3/bLsA6quga4XUS64f5v98l3cEOUaqmMI2TLG22HnL61InIX8O18lSKL4GVc6tyYfDuISDucvMAXgTm4Oi/jVfXVYoN/wsbMpRe34OZHmyx23kPB0rh9qNP4HmJEJC8ccXbavmRqqOb2p4HDAbhy9yCdP3L5tWl7Rie3Kr2oJpQjvHThBeGiLnFPfa4KPFYb1PfS+t0mBm36SUODXPl2SWnBvmOCNsq7I0KkeOvOLvpcYyos9v7IWDKJVbHaiEx2nmaC8t6AszTeQabtrQZLhrT7ekvCyuHbQHYWfTefTIsPpFtLxlpIWcmFtu1/LpZnY3/zPtCedQyEAL6hTC/tGqwSb4lYuRTbfjDhPFtLxCat+n0XnGhMLSOXwn7GumiQtAhAN7N9UWGZlYjmh64VVn5c0TItFpsAS1V1rYgsJvdX+DKwezmDFar8aLEv8D/gDzgd2B2T52txP/PDgZ8BR+GsjbKhqtNLKPa1FzBDVWcmE9eduOBSRERERMVAVVizul1JjwrAG4CnuL4AnC4im4jIxsDJ5Cr2FUWplsqZwJ2qek6DC5ksIsuA01T1KBGpxZXzOSdrkCbANuRW/p5D7to5ByJyGnAaQG3vzdmfJxiQLJs/y0vpfmcODdaHXzH+9safpX2jTr4gbT8zMMRMtsKtGK1Y3/hOQcSvb1Io81s3muRVm7x4lXP027r3w1aGTMrOYxrHNAZvGYQwfUzi/RFBzqPb0rAU/2iAWzO0/2ywPmpuCof71Xn78Y37gFAk2tZ9NyKNSwck0vcPGsPUWip+0W1puFlxDttvLQn77dyywXPDc2UJSmZ9u/NRfv35TZzHpqP5BND3arqlfTbh0FOBreROlozK5KX7pX0r7jFxknHJs6VZ25jI5Kzkw2wrxErTR7QSlEqZMErBnbgI4S3A+cAjOAnUtbhf0fByBit1UvkizirJwuMElahJFGAJiMijuDLEDXGuqv6thOvI+rnklQRQ1euA6wB61vVqG9IBERERrQ5VYfUn1TGpqOqVpv2siOwKHAp0BB5X1ZfLGa/USWUVLtXtsYxtA5Pt4NxpH2XsA4CqHpRvW4mYQ2562raUaJpts3Ael4y5iK+NGAfAXYbQYMvxnogLQNx3cii29RXuSdt7jjbvbyId8uehgR1mrRYv6fL+ycGSsH72FxL6kI+dAHQeHayKZy8d0GhMyyjxMZUuE43cvflEFx3gVsc5Rbys4KOPtdgSwdYS8DEPGxM5ODTTREnrgbXiFt71b2Mf1jqw/UsytlurJuubai0g/y2w12KP90yurDgMsDKxdJZ1CkGTleZi/Oc224iP2s/Sy8xbEcclT5r1k7eQwlcpWCcAi/ybZRMTI6oXwto11UmuVdW3gevX9fhS7/ovwIUisgb3s1iAc154xtbYZL8BQLG4yPrgeeDTItIXeAf4OiRFTSIiIiIqBQpUj/urSVGSoKSIdMTNXN/I2Hw7cKqqfiwihwPLVHVS2RcichTwe6A78AEwTVUPEZFeOOrwYcl+hwFX4SINY1X1klLG36Zua/1O/fEsSIT5rp4RvHRW2K/z4sRSeDYc+9pRQablHo5J27541+R2wU9++JRgzP15kLNg2hvBSssq22ZSklxh9aDN+mDSAY7pNWh5iKPUGKL1ysQqeL+TYXwtNFaJf2d+QCYm9XXj2zyYbWYZGXy/+g8amaw2lkR7bwkEcloud++U5NkytmxMw1oVWZaElXfxsQ4b87ExFR+fsdL3xmZemcRHVrcLn/W8mnABHZOSz2+aPJIpGeG60UtDbGOz2mBBLngpCSBZIvw40zaxEl+WN6Iy0RSCkrJrnXJfZlZEY3xm/c9XSSjJUlHVFcDxIvILnNNna1yNvymq+obZ74E8Q5RyjvuB+zP652IcNKr6IKFSR0RERERlYnXxXTZElOX0SyaQN4ruGBERURRtSU6+zcEVVKlYiEhvYJ6qftLUY5c8qYhIJ+AkXH5KV5wwxpPAOFVdXuDQisBaNmIFnVIasA3vd34qBMdfOdm5P1YdFRLbbLGd86aGRMkZA536r02kXGlUeL/1WkIlvsNcyPdC863BzmfT+5zgx3ntgOBqG7DGRchrQu4lz/4+SHfsPcFt73SYCdRbi9sfZ+RGbL37fZY6t1p7k4T41u5G5XjT5LqMR6x91jfGjL/S0I/bJT+q9sa1rMb9lVOo2rctzcMmVb6Xcbx1lSXXsNK41Kxb0KtSv0egBOfU0Unke6wysE1O9JUXO2wSLnrBYJOomFJ+rTRKvqB7nDw2eLTwpCIi/YBzgc1V9Zhi++NoOJ8HnhORx4HvquprRY4pCaXKtGyNm0B2wCVBvosT8z4GOENEhqhqRppvREREQ+Qr0+v718ViiVZPhUGB9bQBRGQsrtTaAlXd1fQ3kqpS1ZnAySJyT/ZojbAC0prpQ8gtWbdeKNVS+TXQBdhPVZ/xnSKyD3AvcBllJsi0NLZZNo9LHr+Ihw8YAsBLg0P0+cHBgVM76kRXzfHZm41F8CvDqTWr57kD3UrWikC+1Cn97Jm8o6vM9+N9Q4rP893D9l44eZV3TrY8XjN+O1dSsXZkqOGSI1KXDFU7zSQfmkC6rwzZe1awhGpvMeKR33BB6//sHhL29pgXIsrv93RU6A61YXXeeVaw6lYmC/XlnQJlusu8YDW929NZClu/F8gDYi2RrORGg49qDYGCtY22L+oZ6L/eWrR1ciy81fGCUYFckf6mQv9Uw4lOg+8Q6Nd26ZQjl3J6RhXFxn/s+f78Yz34DQxKLqV93TAOp2Jys+9YH6mqBngB+J2I/CN5fYaIzMuzr6rq2Xm2NUKpk8pQ4Gw7oSRn+qeInAdEUceIiCaCnWCixVGlaAL3l6pOEpE+DbpTqSoAEfFSVeVOKqcClyfHKk6pL980qARnelGUOqlsSv4kwznkViavSMzdbGsuPOBbacKgxTBCNUVPv7Ur3kU/C7f3KCF/01sots9aEt8kqVFtKgR2MwEK77/f7cXAfXh49xCU8cmTL+wY9ORtHOC1vi7+0qlvsA7mG+2THp5za5YCs4aHgvJb4KTvBywO1old/Xeb92Gjvs6dA3+4Jk04DOdXYyn1WOwsFBvH+aBdeF+3WxwsqFVJfMRSflfUBEuiJjnvq7WG32zgExFtLfe5hHv1VolNep01v0/aXvtRcq7hRrQhS7gRU9f95tLqu+cemz1JrK/FEiefCkOxSeXB69zDoVuBPRsiU6pKRLbEJRHsISI/U9VfFbw8Fz8ZBiAia4Evq+pzhY4pFaVOKq8DJwAPZ2w7nlAMNiIiIgPFrI/o/trAUGxSOfg09wA4XBaVMXKmVJWqvke+AjrF0RfI5/oqG6VOKlcAN4tID1yy4zxcrsrXgYNwE05Fo4aV9GFWyv5ZYFb0h054MuyYuNytxXE/R6XtU68MOvLv/tjFDI6dEiyd1aZI1ZRaJz659YFT0z4rOFiTWJurjdxIf2akbb/q3t702cJO3hKx1knd4iAjM7trslL/Qhjf3tfyZPXeZUmwNLqtCZbIyoQ81e3l0PfuroFRVWNk3APCWD7WstnS0NeuJvx+vEijxeyavo36AKbUDmp0/TPob+7FnWuBsVTs+3LPfEeISS0SgGnm9/lB8pwj3Jj7G20cM2keRKtjA0Dzsb/WWaoqH1T1fyLSXkS+Bvwfjt27GJgM3KeqZd1JqcmPtyaU4l+AKZnowpanq+rt5Zw0IiIiYoNHsb/i5ydA/QTI1ZkohiaXqhKRrYCJwG7AbNz/+udx9VhfFJGDVUvXgShJpsWcfCPgM4SZ7HVVbUzLqUD0q+uiv6wfwmf5D5ArW26LXL0wwpkaezwU4gz3DQ3iktZSGJgkhbyXxyXqYy4+B6LhebdKLI0VBPbU6gzByB3vDOyvSV8PRbp8zKX/vFBw0zOuIMSFLMtph5VBms1L0thjrPXhGVXWOuiQaZ0EbLMwxIze6e5iRlusfD+MWRNyetaYNY2XR2mXo7MSUI+z+qwlkiPouMoJOi57P4y/9kWj6XJx8pwTJ8nNI1EdHam5EU0j09KvTrm4RJmW47LPJyJ34Oi+3XB/9KNV9cZ1larKe60it+L8GUer6vOmf08cu/cpVS3ZG1VuRv1aGlR8EJGDgCtVdbdyxoqIiIjYYNE07K8srcXmkKo6DPienVCS8zwvIj/DaTKWjKbQZt4c2KUJxmlWKMJKOqQy8raAUt8lxlJhDwBeHbpz2mcFF2cbvXW/qrdZ2CNn3Ji2J/V3VoWVrrcrcR9T6UaIM2y1JjCi3mznVu8Pf31I2mdl8v1Yk3oG68VaEt4q2mNGWAc83z/kyXTr+V6ja+qyPOSUvNdpy0ZjWqvFt7daGVb/73a3Vo+7vzdqgnVn3wubse7bNuZkLbhnkkJmObkl80POydpxiVVi1ecG5bb1chBpbJ0Ueh0RsU4oZVKZNgFeLNv91RyoAfPDzMUyoEOebZkotZxwRERERESp8JNKoceuw+C46yC3Xmlr4FngbBGxGuAkr88mR7O9OKqzikxExAaAGL/ZwFHBgpIN8BPgCeBtEZmIi99sBRyCozAPKWewNjOp1LKUg3iMbV5zLp+f7WgIE6YAsg+ub3ONUVE0TNK5x4SEOk/53SLlo8K7/YMl6ym7vQzjb49ZxhXV17minmD/tG9Nu+BeOv21cQA8vGO4AOtq84H0fMHtlYnV+lr/IFJp3UfLE/rvcuNmWtOpcWGhLdaE+3u7XQiO+/MvN0mK1j1mXVkeNqnUVlH07kgbfJ+1Kmxf8tukimLgTOSIc6ZD9TF9VzTI/bg8/mlHtBAqXKXYQlWnicingTOBPXEssHk4JsuVqlpOHk3+SSVRvSwFWTXnIyIiItou1mJTtrLx6gR4rSJiKiQTx6imGCsvpThJ3S+FbyzumrSia2fuXtdeH6zfIrVUntjx8+m2fZb/K237GuWW2mtXz3veGJILR5zsxCfHTPlh2vfwoCFp+9AxTwJw54gvp32WUnzEwokA/Ll7qHHf2ygweEvHHmNlZnxwe6WJo3nKNMC0hHTQw6ggWpn+HZLKz6tyCsM3hrWEVplzeasjVxrFBt+DVefxzyTgDtDJ/Orun+/eo5zkRItBSaKiSS7NoQd3S65hUXQpRawfmoRSvE2dMrJESvG5bafy47cLbIuIiIiIKIQqcX81NfJOKqr655a8kObGh2zKP9mH13d0goQ2tjClU+CeDiApjGUEO/dfGCyZv5w8LG37+Mszgwamfa+bgu7LR7hYhbUOrHTIw92HANDfUJatVeDb9nhrNW2asACtpWIpz97CsWNaenCWJZMFe4ylYvuYiLWObMzEJn1mXd+yVUGocu1ZCfHEktNH5V6Xao+SpOOjdRLR6qiimEpTo80E6iM2TEQhxoiKRJxUNnwspxNTGcgP+S0AfzUikYsMo2pKkjG3BUFaZH73bPbVCWucuORl7UKpgY4m/uGtil5GALSDsYC8NWQtkVcJSZf78M/k+MAeW24srA9TGZXwMVqhTG9hWOuhD7PTtk9EtNaVteC2TGT6p5gsQjuWZ41Z6f/6VcFq85Ipax8w9PcnycYtjS0NGdWwHO/oHCskUnIjKhZrySmznYk3J8DM1g3Ui0gNjvX1d1V9sSnGjMmPEREREc2BYsmPnxoG+7du8qOqrsTVts8umboOaDOWymYs40AeTVf3W5k4wjGEss5e5n4Y49O+GzglbXcxOSkPtnNliL/M/WmftRo8+8mu/i3jya/wBzEl7bNWibd0rIyLtVS8pWHZYTYW5BlsVmPBsrN8LMVaRxZeWt4yvqzV4mX6bb7JkicNw/zW5PmhzOFhUbFyvIURLZOIikV1ub+mAAOBp5pisDYzqURERES0GBT4pLUvomT8FLhdRFbhhCrn0yCdRFWXZx2YhTYzqaygI/9hN7748tMAjNz1mnSbjWn41f/tHJf2Weukp7EkvFXyoTneFo7y8QfLvvoMQXreM828BD7kxmy8VWLZV4uMzL7PD+mVp0aPvy8bB7HX4q0OK/JoYzKevWVzY/41JWT//+vJ0E7xQOMuhue+1Mvdc0NxR4jijhEbCBTyCF1UIryr5Grgd3n2KTkPsc1MKhEREREthlLcX29PgDkVkVF/EqUlupeEOKlEtDgiDThig4dSXKal2zD3mHF9q6oUq+q4phyvrMqP1Yxd6jbRO+o/lSbvZbmBAAbwAhDqv0NuINsm//nkR+uyWmkkT7zbyLrXbKKhp+xOICRU2tot2/FWo/uwY/l6LvkqM05mv0bbO5pv+hsJldi6ty57ybicNk2+Gx+aWu4Wu/vJ4fTs7WQnLOZujy6uiMpCk8i0dKlT9i9RpuX+ypBpEZGdcQH77XAVJd8Vkf7AfFXNV2+lEaKlEhEREdHUqCL2l4hsCowFjsHRC9oDDwPvAr8E3sLlspSEiplURORY4AKcZOBeqpo5zYvIbFw1sjXA6lJn+I7LV7Lb1DcYM3AEALMMDXYk16btGUmtdFuL3iY07sfktO2l3esIl2olUz7DG0C2XAkEq8NaCr5WO4QA+po8MTJfr91Kt9hr9VbLXwmClrZy5FvJcZMYHAZ915zg4sRC6WP6ZmddyZ8KWhzRGoloc6iiSQW4EtgHOBB4hty0zQdxE0r1TSrAy8DRwJgS9t2/XI3/iIiIiBZDdVGKjwZ+oKpPiEjDFez/gE9lHJMXFTOpqOp0AJE8/vv1xLxOPbh44PHp61O4IW3bmIqPtVjroMbEKaYQ6sEP4jkg11KwiYDekrCUXR/ngCDEaIt8WUvCWy0+9tJw33rqkr4gKWOtIk9Jnj5xj7Rv1pBwfWmi4hwCFpr25Av9RedFtEIiIvKgeijFHcH8yeRiM8q8k2qUaVFgoohMFZHTWvtiIiIiIhrBa38VesydAK+cBq1PKX4eODHPtmMgESEsES1qqYjIo2RXijxXVf9W4jD7qupcEdkK+IeIvKaqk/Kc7zTgNIDa3puzgk6pJTLZxBHs6n9FIpJoWVJ2u4+TQG78w2OIUUz0Vost12vjL3ZcD5uo6GVWcmVSgqXkExVt4avX5wdxyLXjEiHHT4fxl2xj3v5FWeyt3ITEaIlERKwDSnF/bTbMPRa0LqUYOA94NPl//gvu6g8TkR/hJpXBhQ5uiBadVFT1oOJ7FR1jbvK8QETuB/YCMicVVb0OuA6gZ12vtsGdjoiIaH1UUUa9qj4tIgcClwJ/wFXzvRB4FjhIVZ8vZ7yKiamUAhHpDGykqsuS9sHAL0o5diU1zGD7NP9kSyPSuMCUw/Xxj2WEAlJWhNFaDV6a3loKNxlNkjqmArly81Zm3sus5CvBu1sSc7Hj22v1sZwFY3uHG7UxEW+hXEURROskIqLJUT3sL1T1GWA/EekIdAE+KEfvy6JiYioicpSIzAE+DzwgIo8k/b1E5MFktx7A0yLyIvAc8ICqPtw6VxwRERGRB55SXMqjsvAxznFXTA8gLyrGUlHV+8FoyIf+ucBhSXsmsHsLX1pEREREeSilSFcFQUQOw8VWBuLmhdUiMhW4RFWzZGLzomImleZGB1axHW+nNF4bEO9m2HTeFTXPuKT6GvLDpoYe/ARDgNyAu3d5QaAn2+D6KiPjkqVibOuxPJO4vWw9lqkv7Rtu6pbkuUGJeb3ZPUt305e4xWK1xIiIFkAVJT+KyAjgWuAx4AfAAmArXP7KeBH5rqqWkj8ItKFJpS3BTiYRERGthCqZVIBzgOtU9TsN+v8krj7FuZSWlA60oUlFEVZRkwa37erfijB6cUhb48RaErPpm7Z97Xmb/GjhKzLaY2ygP9RLCdUapzEgbS+YmYx7gUkIvcWYJfslQfvJxvro1qAmyUJyX0fLJCKi+VEKpXjlBFhVEdL3WwL35dl2L3B8nm2ZqJhAfURERMQGA08pLvRoPww6tW6N+gRPAF/Is+0L5EnZyIc2Y6ksYzMmsx97JUXOrLS9pfH2Tei/1jqZb2i8uVUUnX1rLRVbL76egY2uw1otc2YmyZN/NJZIN7PzFcnzojzS8ZOd5ZGjbLPowmiNRES0Nio8ppLI3HtcDdwgIlsCfyXEVI4Chv5/e/ceLGdd33H8/eEmisiAKRAgEjEWLIwdB4QpTgcwYKkQDCgVHGaINoOM2tG2ysW0AmorKDOoBSuIFgqKQDq0wYarkUumgATKRQyRBFIICeGixkQgkfDtH7/f5jxnz57dZ8959nbO5zWzs88+z2+f/f02l+/+7sDcdu49aYKKmVnXvM44BuV2xS8YvtujgE/mR+TXNTfj7YRH2oHfczD3bRldVZxwWFQ7X6ypFPdtr40eA/gZaY/24uKTxQUla301D/1uqJ/klaU7D33YJfn5pkIGXmw0Oqt5zcM1E7M+1N8z6o/o1I0nTVAxM+uqPl4YKiLu7NS9J01Q2ZrNw5agL47+eq1Q01ibayXF68V5KMWthWsLSj61dvqWc6+v3WHoQ2vzSO4rZKQ4Uqu2kOOUoT4bppwzYsSWmVk3SNoGCkNUs3aWbPHoLzOzSUzSTpK+I2kNaR2A9Q0epU2amsomtuMZprFNbugszowvznKvXS/WahZsmrXleMZ2K7Ycr3hk/3RQ7BMpLmiwND8XN+I85hzizOF5k+qmxBdGm5mZddgVpKHD3wOWQ2Hi3hhMmqBiZtY9/T/8q2Am8MmIuKaKmzmoDACv12U2aAZqk/qngTEtc9/IpAkqqaN+A7vl1ReXF3ZtXMGMLce//F3qiN/4ylBf1RveOFQbfOCiwoKOtaau0ZZJqQ0PPqtVIPhu3WsHDbPB1uezH4c7AzhP0v9GxNPjvdmkCSr9qL4GYmYTRZmayq350du1vyJioaQjgeWSVsLIfc4j4uARbxzFpAkqG3gz/8OhW3Z8HLac/JNDOyvyQp5IesnQqVcoDBNeWbjp3bUO9sIe7w0nL7bWLJ2bu8wGTZmgckR+XN3Ttb8kXQh8Drgfd9SbmfWjYIA66ucC8yLia1XcbNIElT/EtqzeNJWld7wnndhQuFjcw31Oft6/cO5ChhmanNho6G/5WoVrIGYT1UD1qbwMhd0Fx2nSBJVO8ugsMxtuoEZ/fQs4TdJtETHuxWUmTVDZvG5b1v1kd7gnn7iwfsJhdnf9SKzEgcLMyhuomsoU4BBgmaQ7GNlRHxH1U7ZHN2mCiplZ9wxUTeUjpAi4LXBUg+sBOKiMsGI1fPg8hkZqNa6RQG3jq/LDfV2LMbPhBqemEhFvb52qvMkTVMzMumaglmmp1KQJKgceuAdLlhRrFN74ysw6ZXCavyR9qlWaiPhO2ftNmqBiZtZdg9H8BVzc5FptNFjpoDJp9lN54IHVSOd5aRQz64JaTaXMo7ciYqv6B7ALcDLwMBR2JizBNRUzs8oNTvNXIxHxW+BaSTsBlwKHl32vg4qZWeUGZ/RXC08BB7XzhkkTVEZ21JuZdcrgj/6SNBX4e1JgKW3SBBUzs+4ZnOYvSS8w1CFfsx2wI2nP+hPauV/fBBVJ3wBmkZZdXgF8PLfr1ac7mrRWzdbA5RFxflczambWUnebvyTtQBqhtQm4IyJ+2MbbL2FkUHkVWAXcHBEvtZOXvgkqwG3A2RHxmqQLgLOpWxpA0takL+AoUoHvl7QgIn7Z9dyamY1q/DUVST8AjgWej4gDCucb/bA+AZgfETdKuhYoHVQi4txxZbRO3wwpjohbI6IW2u8F9mqQ7GBgeUQ8GRGbgB8DH+pWHs3MyqnVVMo8RnUFcHTxROGH9V+ShvqeLOlPSP9fPpOTba6oEGPSTzWVok8A1zY4vydDXxyk2sohXcmRmVlp4++oj4i7JE2vO73lhzWApNoP61WkwPIQJSoLkha1l5WYWTZxV4OKpNuB3RtcmhcR/5XTzCOF70bVNzU4N+r6/5JOA07LLzdIWtZejisxBfIexhPDRCrPRCoLuDxV2Xv8t1hzC5w7pUmCKcAf5ePt27jxaD+svw1cLOkY4MYS9ynTTzIVOJQm/8c20tWgEhFHNrsu6VRSG+LMUTaLWQVMK7zeC1jd5PMuAy4bQ1YrI2lJRLQ1zrufTaTyTKSygMvTTyLi6NapxqThD+uI+D3w8bI3iYgTR/0A6W2k/uxjSUH9onYy2DfNX7nz6UzgsIh4eZRk9wPvlPR24FngJOBjXcqimVmvtfXDuh2SZpAGSJ0CPJ+PL42Ittrx+qajnrSo2Y7AbZIekvRdAEl7SFoIkDvyPwPcAiwFrouIx3qVYTOzLtvyw1rSdqQf1gvGc0NJ+0v6Een/1COAzwLviIhvthtQoI9qKhExY5Tzq4EPFl4vBBZ2K18V6GnzWwdMpPJMpLKAyzOhSLqGtObWFEmrgHMi4vuSaj+stwZ+MNYf1pIOBOaROvp/BcwFro6IcY0eUwX73JuZ2QCRdBPwAeAR4J8j4vrK7u2gYmY2uUh6PR/+mjT+uamI2LXsvfum+WuikHQicC7wLuDgiFgyLYlVawAAByZJREFUSrqVwHrSRKXX+nWUSxvl6fvlcyTtQpr/NB1YCfxVRPymQbrNwKP55dMRcVy38lhGq+9a0huAfwcOJA0d/WhErOx2PssoUZY5wDdIA3MALo6Iy7uayYmpYxtLuaZSMUnvIkX+S4HPtwgqB0VEX88rKFOePMv3VxSWzwFO7rflcyR9Hfh1RJwv6Sxg54g4s0G6DRHx5u7nsLUy33XeHvbdEXG6pJOA4yPioz3JcBMlyzKH9O/kMz3JpLWtn0Z/TQgRsTQiejHJsiNKlmdQls/5EHBlPr4SmN3DvIxVme+6WM75wExJjeY39Nqg/L2xNjio9E4At0p6IM/8H2SNZvnu2aO8NLNbRKwByM+jtRNvL2mJpHsl9VvgKfNdb0mTh+GvA97aldy1p+zfmw9LekTSfEnTGly3PuI+lTEos9xMCe+LiNWSdiXNzXk8Iu6qLpflVVCetpbP6aRmZWnjNm/Lfzb7AIskPRoRK6rJ4biV+a775s+jhTL5vBG4JiI2SjqdVAN7f8dzZmPmoDIGrZabKXmP1fn5eUk3kJoCehJUKihPx2b5tqtZWSStlTQ1ItbkXe2eH+UetT+bJyXdAbyHtMdPPyjzXdfSrJK0DbATaZRPv2lZlrq9PL4HXNCFfNk4uPmrByTtIGnH2jFpvPgvepurcal8lm+HLABOzcenAiNqYZJ2zqOnkDQFeB/QTwMOynzXxXJ+BFg0ylp6vdayLDn41xxHmvVt/Swi/KjwARxP+gW2EVgL3JLP7wEszMf7AA/nx2OkZqae532s5cmvP0gaybOiX8tD6lf4KfBEft4lnz+INJwV0qqsj+Y/m0eBv+51vhuUY8R3DXwZOC4fbw9cDywHfg7s0+s8j6MsX8v/Rh4Gfgbs1+s8+9H84SHFZmZWGTd/mZlZZRxUzMysMg4qZmZWGQcVMzOrjIOKmZlVxkHF2iJpjqTIW4/2Mh9X5HxEnqBYOz+9cH7EzGtJi4vp87mQ9NVG1yXtVbhfSBoxuVLS2ZKeaJavQtp9JV0p6VlJm/LzVZL2HaWMq/Lx3HzP6YXrtxc+6+rW35pZ5zmo2CB7Dvgz4FOjXP+nCj7j+fwZn26SZjbwn63ylQPSg8CfAl8EjiTtA74/8GCjgNXC3+TPea7N95l1jJdpsUG2MSLuHeXarcAHJM2KiBvH+gGRVs+9V9L2ja7nGd/vBf6uWb4kvZW0Cu/DwPsj4tV86S5J1wGLgB9L2jeGL03SLG9L8703tlMms05yTcU6QtIpkh6W9KqkF3MTz9S6NG+S9K+SXpK0XtINkg7NzTlzxpmF+aRawVc7vOz7bOAF4J4W6eaSZvR/thBQAMivP5evz+1EJs26xUHFKpeX8r+KtE7TCcBZwF8Ad0oqbn51GfAJ4MKcbhnww4qyEcA/AO8GOrlB1WxgQUS02pJ1JvBcRNzf6GJE/Jy0DI5X4LWB5uYvq1Teze8rwB0RcVLh/OPA3aQg8u3cMf0x4KyI+HpOdpukN5H6CsYtIm6StBj4sqT5kfYWqYyktwCHkwJiK9NIWxg3s5Lhq/aaDRzXVKxq+5I2vxpW44iIxcD/AYflU4eQ9tO4vu798yvOzxeBdwJzKr4vwDHAJuD2EmnLNMH14+6MZm1xULGq7ZKf1zS49lzheq1/pX5Pk7VVZiYi7gZuBr5UW9K+QrOBmyOiTEf5M8D0Fmn2ZvhOiGYDx0HFqlbbDKrR7ou7A7WRTbWgU7+l724dyNM80gZQp1d1w7z/x9EMH0rczE+B3SW9d5T7HUwq+6JqcmjWGw4qVrVlpNrGScWTkg4l/RK/M5+6j9SZfmLd++tfj1tEPAj8B6kpbIeKbjsTeCPw3yXTXw78BvhW/fDk/PqbpIB8eUX5M+sJd9TbWB0tqX7S3bqIuE3Sl4BL8yzvq4E9SRMRnwD+DSAilkn6EfAVSVsBD5BGPs3K92o1mqpd/0jaXXNXhgLbeMwG7oyI35ZJHBEvSjoZuAG4R9JFwFOkJrG/BfYDji87R8WsXzmo2Fj9S4NzjwEHRMRlkl4GvkDasncDsBA4IyI2FNKfBqwHzgC2IzX9fBr4CbCuysxGxOOSrqKCDvs872UWbc7Yj4hbJB1IqjGdD0whNQcuAk6JiH7atthsTLzzo/UVSV8ALgCmR8TTTdJdQRrOOwOIiNjcwTxtQxq1djtwFClI3gNMi4hVPczXVqQm7OXA4og4pVOfZVaWayrWM5KOBQ4AHiI1d/058HngumYBpWBv4A+k5qzDO5THvagbkZWXYGk2/Lfj+cpuJfXtACzu4OeYleaaivWMpMNItZL9SB3ozwLXAufUL2XS4L3TSc1HAOsjYlmH8rgtaQHImmURsb7X+cqf9cfAW/LLlyLiqU59lllZDipmZlYZDyk2M7PKOKiYmVllHFTMzKwyDipmZlYZBxUzM6uMg4qZmVXm/wHw3j+RzXa1NAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "############################################\n",
    "####### 消光修正 extinction correction #####\n",
    "############################################\n",
    "\n",
    "# 元素的波长设定\n",
    "# ha_6564    hb_4862    nii_6585    oii_3727    oii_3729    oiii_5008    \n",
    "ha_lam = 0.6564\n",
    "hb_lam = 0.4862\n",
    "nii_lam = 0.6585\n",
    "oii7_lam = 0.3727\n",
    "oii9_lam = 0.3729\n",
    "oiii_lam = 0.5008\n",
    "\n",
    "# 获取元素的flux\n",
    "# %xdel ha\n",
    "# %xdel hb\n",
    "# %xdel nii\n",
    "# %xdel oii7\n",
    "# %xdel oii9\n",
    "# %xdel oiii\n",
    "ha = spaxel_data[:,5]\n",
    "hb = spaxel_data[:,6]\n",
    "nii = spaxel_data[:,7]\n",
    "oii7 = spaxel_data[:,10]\n",
    "oii9 = spaxel_data[:,9]\n",
    "# oii9 = spaxel_data[:,7]*0\n",
    "oiii = spaxel_data[:,8]\n",
    "\n",
    "E_BV = spaxel_data[:,11]    # 获取E(B-V)\n",
    "Es_BV = 0.44*E_BV    # 计算Es(B-V)\n",
    "\n",
    "# 判断是否 ha/hb >= 2.86，若是，则修正；若不是，则不管\n",
    "# %xdel ha_over_hb\n",
    "ha_over_hb = ha/hb    # 计算ha/hb\n",
    "tmp_id = np.where(ha_over_hb>=2.86)    # 找出 ha/hb>=2.86 的数组位置\n",
    "Rv = 4.05    # 论文中的Rv取4.05\n",
    "def flux_coeff(tmp_lamda):\n",
    "    \n",
    "    k_lamda = 2.659*(-2.156+1.509/tmp_lamda-0.198/tmp_lamda**2+0.011/tmp_lamda**3)+Rv\n",
    "    # 修正系数\n",
    "    x2 = 0.4*Es_BV*k_lamda\n",
    "    flux_coff = np.power(10,x2)\n",
    "    \n",
    "    return flux_coff\n",
    "\n",
    "# 进行flux的修正！！！\n",
    "nii[tmp_id] = nii[tmp_id]*flux_coeff(nii_lam)[tmp_id]\n",
    "oii7[tmp_id] = oii7[tmp_id]*flux_coeff(oii7_lam)[tmp_id]\n",
    "oii9[tmp_id] = oii9[tmp_id]*flux_coeff(oii9_lam)[tmp_id]\n",
    "oiii[tmp_id] = oiii[tmp_id]*flux_coeff(oiii_lam)[tmp_id]\n",
    "\n",
    "\n",
    "# 计算元素的flux之比，并log\n",
    "# %xdel lg_oiii_oii\n",
    "# %xdel lg_nii_oii\n",
    "lg_oiii_oii = np.log10(oiii/(oii7+oii9))\n",
    "lg_nii_oii = np.log10(nii/(oii7+oii9))\n",
    "############################\n",
    "####### 消光修正结束 #######\n",
    "############################\n",
    "\n",
    "# 画图\n",
    "# %xdel y\n",
    "# %xdel x\n",
    "y = lg_oiii_oii\n",
    "x = lg_nii_oii\n",
    "\n",
    "fig = plt.figure()    # 设定图形大小\n",
    "ax1 = fig.add_subplot(111)    # 设定子图，调整位置\n",
    "\n",
    "# hist2d画图\n",
    "xbins,ybins,hist,cax = ax1.hist2d(x,y,bins=100,range=[[-1.8,0.8],[-2,0.8]],cmap='jet',norm=mpl.colors.LogNorm()) \n",
    "# xbins,ybins,hist,cax = plt.hist2d(x,y,bins=100,range=[[-1.8,0.8],[-2,0.8]],cmap='jet',norm=mpl.colors.LogNorm()) \n",
    "# ax1.set_aspect(1)\n",
    "plt.xlabel('Log [NⅡ]/[OⅡ]',fontsize=16)    # x轴标签\n",
    "plt.ylabel('Log [OⅢ]/[OⅡ]',fontsize=16)    # y轴标签\n",
    "# colorbar调整\n",
    "cbar = fig.colorbar(cax)\n",
    "cbar.ax.tick_params(which='both',direction='in')\n",
    "cbar.set_label(\"Number of all selected spaxels\",fontsize=15,labelpad=8)\n",
    "plt.show() \n",
    "# b = np.argwhere(np.isnan(a[:,8]))\n",
    "# c = np.unique(a[b,0]).astype(int)\n",
    "# print(lalala.shape)\n",
    "# print(\"max = \",np.max(lalala[:,2]))\n",
    "# print(\"min = \",np.min(lalala[:,2]))\n",
    "# print(lalala[20:40,2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAETCAYAAAD+spv+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsnXe8HFX1wL8nIYQEAiYEQigxSCihCSSAINJFQENRECkqCgJSpag0BdQIWEBpIu0XUSlKM5EqvZfQ0UCMEJJACEJCMxAgOb8/7szcs+/N7s6+t+/t7nvn+/nsZ++7M3Pnzr7dOXPuaaKqOI7jOE496NPoCTiO4zg9BxcqjuM4Tt1woeI4juPUDRcqjuM4Tt1woeI4juPUDRcqjuM4Tt1woeI4juPUDRcqjuM4Tt1woeI4jtOLEZHRIvIZ8/cAEfm5iNwgIkfUOp4LFcdxnN7NBcA48/evgKOAJYAzReT7tQzmQsVxHKd3sy7wEICI9AP2A76nqjsCJwLfrmUwFyqO4zi9myWBd5L2Z5K/r0v+fgL4ZC2DuVBxHMfp3bxIECYAuwNPquqbyd9DgXdrGWyxOk6sqRk6dKiOHDmy0dNwHKfJefzxx99Q1eU6M8YoEZ1fcN/ZcGuy1NQozgZ+JyJ7AhsC3zLbtgaeqWWwXiNURo4cyeTJkxs9DcdxmhwRebmzY8wHDq6yzwvJazas1tnzdQZVvVRE/g1sDByvqneYzXOB39QyXq8RKo7jON2FUP3muk7yegLe7voZVUZV7wXuzek/tdaxXKg4juPUmT7AgEZPogIisnYt+6vqv4ru60LFcRynzgjQr9GTqMxzQJGyv5Ls17fowC5UHMdx6kyR5a8Gs01XDdzk1+04jtN6NLumoqr3dNXYLlQcx3HqTAtoKu0QkZ2AscAqwM9UdYaIbAlMU9VXi47TatftOI7T9DS7pmIRkWHARGAMMB1YFbgQmEGIWfkA+G7R8VyoOI7j1Jlm9/5qw7nAUsBaBKHyodl2O3BKLYO5UHEcx6kzraSpADsC31TVaSLS1strFrBSLYO5UHEcx+kCWuzmurBM/1Dg/VoG8oSSjuM4dSbVVIq8moD7gCPaaClpDMu3gTtrGazFhKnjOE7z02LeXz8E7icERF5PECjfEZF1CbVWPlPh2Ha4puI4jlNnUkN9kVejUdXnCJ5fk4H9CUthXwZmApuq6tRaxmshYeo4jtMatJihHlX9D/D1eozlmorjOE6dSZe/irwajYiMrrL9i7WM50KlBRCZU/JyHKe56W5DvYiMFpELReQaESkcqJjwhIgcJyLSZsxBIvJ/hMDIwrhQcRzHqTP10FRE5DIReV1EnmvTv6OIvCAi00TkeABVnaKqhwBfJaRaqYWTgZ8A94rIp5JzfJ5guN8BGFfLYM2gffV4RE4r+Vv1lNz+wCHdMCPHcbqSOtlUJgDnAZdn4wa33/OBzxMCEx8TkYmq+i8R2QU4PjmmMKr6axG5CfgD8IyI/APYBbgSOEJV59UyXtNpKnlSuM32/UXkvyLyVPI6sBHz7G5Uh2Uvx3GaG6Hz3l9JNca5bbo3ISR4fFFVPwSuAnZN9p+oqpsD+9Y6X1WdQhBIfZPxngQOrVWgQJNpKpWkcJtdr1bVw7t9go7jOAUQoF+Fu+sfFoVXwtAahl6J4OqbMgvYVES2JrgB9wduqmE8RKQ/cDpwJHAdQUM5G3hORA5Q1X/UMl5TCRWMFAYQkVQKFy5l2QrkL3s5jtNTEIHFKtxdD0heAMt8wBu1DJ3Tp6p6N3B3DeNYniIItv1U9SoAEbkd+C1wi4hcpKotm6U4Vwrn7PeVJM//VOBoVZ2Zs0/TkNpQoIBA2bt0eUuv6IoZOY7TlYhAv8IFeGtiFqHeScrKQOFaJ2WYBmyjqq+lHar6LvBtEbkO+D0tnPo+Vwq3+XsScKWqLhCRQwjGpW1zBxM5CDgIYMSIEfWcZzJ+dY3DCpQOn+eYNmOe1ekhS8cv40jgOE7HqKapAPz9I7jxIwCWqWHox4DVRWRV4BXga8A+HZtlQFXLenep6t+TdC2FaTahUlUKq+qb5s+LgTPLDaaqFwEXAYwdO7atcGoI1W7YMqSbJuI4TpchAv36V95n9/6wO3Dp67ydP4ZcCWwNDBWRWcApqnqpiBwO3Eowql+mqv+s59zbUquxvtmESlUpLCLDVXV28ucuwJTunWL9sBpCOWEjMgeObozHV7X5WQGobX1UHKc3U4eMkqq6d5n+m6jRGF8NEdkL+A6wBrBEzjmXLzpWUwkVVf04TwqLyE+Ayao6ETgy8cf+mOBut3/DJmyo15KRvTln0fNnz+lSV2Jf7nKcOlNAqEyaD5NCpZJalr/qjojsA1xGiIvZNmn3ITy0v4WJkylCUwkVyJfCqvpj0z4BOKG759XVxPQrF5re+gRC1sNmUkSryjtXpWPcluP0aKrcXcctHV4Xv5i//NWNfB/4KXAGwQZ9gao+ISKDgH8A82sZrOmESm8iRtY3Zz6vypH/vuTlOGURwlpLa7A68ICqLhSRhcDSEDzARORMQszKr4oO5kKlxch7uq/khVZEQyi3X9W5bJ40HioedyOd8lNxnBahtap0vU0ImoRgyx5NjHkRYNlaBmudy25C6rZcs5mxlzxUuuTVDGlZOn6dVZbv9j7F43CcnokQb9NlmPQ2THoHaLBNhVCca32CLXsi8GMR+Rj4EPgx8Egtg7lQaXEybaFe4xVIctkRQZd5iu1Y+5wcp+UooKmMWza8Ln6z4TaV04FPJu0fJ+0LCAt4j5HE+hXFhUoToA/av2q4YW+WaBBHG03i7PosbXWc9gKoJPZmhW6ciuM0ihZa/lLVh4GHk/ZbwK5JPrD+qvpOreO1yGU7juO0GK1jqM9ICnUNBd5Q1QUdGcOFShdivbo6tGRU0I23yBgdcQOuBzK1y4Z2nOalSJzKGzApeFA22qaCiOxMKNY1hjDzj0XkcWC8qt5Yy1guVFqNzXKEw8NxCU0mmO3zujgb8pfS91PQE5Pz26Wuh8381ki2txVgZ3lsitMDKWJTWSG8Lp7dWJuKiBxMsKHcARwFvA4sT0ilP1FEDlXV3xcdz4VKJ7BG8lK7SIVjuijgL9OKBudrRHlajzzRwXOlCS5NFYhs/MFl0s10U04zaeNNpu7C7DSCAt5fTcSJQF56+wtF5ELgJEKm4kK4UOlCOusOnCdw8oRXSfCk9a66MjWa2yj9HKwgmFcaiKk6rHR8m4fspBxNyGpHa+aMP/gUD5p0ej4tZKgnxKFcV2bbtcB+tQzWOpfdg2n7FJ/edHO1j5KbfiosakvnkmkVj1fXkjxrsuN0gNYSKncBWxFSsrRlK+DeWgZrnctuQooueZUc0xXuvVboXJm3vFamSNiYHE2jzPJVRo7LcjxPm/HvScbaqm0ami5MjunLXU4zUCBNy6RXw4vGG+rPAS4RkWWBG4g2ld2BnYADRWTtdOec8u4luFBpINXsIHU/T1eewzUax4kUMdSPCK+LX2p48OOtyfvByUspLZh4S/IuybaK4tKFSjOwVrkNyfLWPKsJGPtIqlXMq0VoHJLZenKj8Uea9qnJ+2tG6C1VqsnoRm0FilmKu6GGaTlOT6K1DPXb1HMwFypdSGfjVOoyh2Ny+tJ5bdb5OXmCSMfJoYVsKqp6Tz3Ha5HL7plkGkO5paNq9o15pyVZistsP7rK8Wck71ajyLP5P2fafzLtkcn7w6aPYbkxKV4rxelVtJBQEZHlgSVV9aXkbyFUgVwbuENVJ9UyXp/6T9FxHMdhsYKvxjMBONr8fRohGHJH4HoR2b+WwVyolEH2KX1167mHFDd8d2WKlVZEppa+HKchpN5fRV6NZyPgTgAR6QN8FzhRVdcCxgPfq2Ww5pCTPZSidpSSuvRlYlYybjTLSF+skqY+jXi33mXnmV2PT95taNM0057evm+7B0vTAN3OF9sLtmSpy5e8nF5Lkdxf/wkvGu9SvAzwZtIeAwwB/pz8fSdwbC2DuVBpctrdsG9srht1nruy21KcXk8B769xa4fXxc823KV4FsF+ch/wReB5VX0l2bYM8EEtg7lQKUPDKhLmRswbppv23skNe13TZ1OnPJdsN4b4X295WNY+dvL5AKx8aFRF3pr/iaz93qxE1THG+TtmfCH+8cn04cZa96ukhOliUicBx2koLWSoBy4DfiEi2xOEyglm22eAKbUM1jqX7bQkqSaTtxTYVstphtLJjlMXWkioqOrpIvIKsDFwBEHIpAwBLqllvBa57N6DvbHK5omm8ZC5+RqtYdn7g4b65hYrZX19Xjsuay8anTSMInHsFefHP3YLb7NmjMy6VhvxQtZ+7/ZEUzF2mKWGvhW3n5/M9X5zAVfWlofMcXokLSRUAFT1cuDynP6af9AtdNk9j3KeWz3RDpF/raXfVxlC9eU/8j+fvnP+l7UXDluyAzN0nDrTHJ5d3Y4LlVbApJvf4ayJWftDFgfghQc/jH2LFs/a580Nbl17X3B91nfPoZtm7ZmsAsCFHJz13X/T5+N5t28/lVUGzszaU1L3MmvTyc2oXMbOMt7s+6v8XRynJWkxTaWe9NLLbm5+zElJa3xD5+E4RckrPtdVBelagtbK/VVXXKg0kLx08T/SqHVktpT94hP9bTftErdPD29L7f9G1vXehFiOce/tg4ayxaGxTMJWNz2StbfbOcScXM43s75ddv5bHH7+qgAcPfDsrO8RoqYzZfsNQ2OyuajjTDuVjXubm8kWZvsnTDsrLmbT+Jcuj1Uy5K84bLb5a1TZ/RynWygSp/JMeNH4OJW64kKlm7ECZKn/HVZhT6eViP/X8u7Vvf7pvTdRJPX9RuF18f0Nj1OpK2UvW0Q+1cExZ6rqRx08tndiEj/+1KSEeUJDXZyNhsSaOH2mRIP0otuDQfq956N2YrWS+y8I9pHvrxENFoN3jt5b/yKMv9BYFJ97YuOsvd1GpdHzAHtwTdZeuEY47o4lYuzKUrvF8ecfsBQAw4a9lPXNnrpqHCwOBSsk72d37KY707WThpJb+ro3C8wCRboaiYh8o5b9E++wQlSSpdMIBVlqZWPgiQ4c5ziO0zNofkP9hDZ/p/d6yemDHHfjclS77PHAfwqO1Re4uOiJyyEiOwK/Tca7RFXPaLO9P+ECxxDy1eylqtM7e16nOZGbwvvwnaO28yqrltm7lJ8ZA8+PZFDuPnl2rdhXOTiz1mSenvyzFyHAEo2eREXsD2It4C/ApcB1xHLCXwG+DXy1loGrCZW/q+qjRQYSkb7UGHlZZozzgc8T8tE8JiIT29REPgCYp6qjRORrwJnAXp05b6N4b8k0ENGsw98Tb1wbyV2hsVncvPaw+FEMPOt9AOYzIOu7/4noEvzyoWGsT06NN8f914jr/Lck1vFR+8zK+na74qp4LsK5HiSWiOzLwqydLpudPuL7Wd/4+Sdl7S2H3QfAx2YdYOwa0ao/6bk944WlSS0fNksmJ9OCXGiEUoOn4jSOJl/+UtVsHV1Efg2cr6pnmV3mAuNF5APgLGCromNXEiqfA/5ZwyQXisjngOeLHpPDJsA0VX0RQESuAnYFrFDZlVjo9hrgPBERVe3IUl2XYsv16oONm0dP4yR+nLXH85NOjZWfEDO/PLPsA1xZS+lmp9fS/Mtflk2A08tsew74aS2Dlb1sVX2gloE6ekwbVgJmmr9ngfFhbbOPqn4sIm8DywJvtNkPETkIOAhgxIgRnZxafbDGy+zmtXe+q+zK+i4Asy6I2/9B1ETuTkpLr8irWd+XN7gua09iHAAHrBFTs5w///Cs/Z2BFwEw6Ip3s75fErWOJwkuwy+wZtZ3NNG9+BTCcs56PJv1HTnwnHbXcTnRJvgW87L2alfEZ5b/PLFOaHzGHDg9NmdLWNI9URdSlLZLXu29r4yGmLo9W6Fhyy2PzDvDIbF6Z87SVq82VDutJFRmAt8Cbs3ZdgDhPlyYZrvsvAWDthpIkX1Cp+pFwEUAY8eObagmU02AOB0j+1wTAbDJg/eW3/eYLjx/Dfu1euJML21QgCZf/mrDicBVIvIcMJFoU9mFYG+pybxQyaX4snLbclBVPaCWE5dhFiS5QwIrg3kML91nlogsRggcalvKqilIl7zK3Xju1N0B2HaIWRubHpuzHkqetP8c+9ZdFJ/u39wvJJKcfUWMInzzqZhc8siVfgvAxGG7Zn37DowZKY9LcqOkGg1EN2OAWwmuwn35OOu7hAOz9u3//RIAey03IesbxHtZO00jsz23x/mxbO65MnYz7fdM+55w8/q5rcJZZSlqtMZA0SnHJIGae5feBNMSB7lxJg9Z7aPazdMTaTqGFlr+UtVrRWRTQtm+vQkO/q8BjwHfVNXHaxmv0mXvQHGX4nppAY8Bq4vIqsArwNeAtsV8JwLfBB4C9gDubEZ7itN7KZI80+nhtFiaFlV9ghq9vMpRyaaycj1OUAuJjeRwwtpeX+AyVf2niPwEmKyqEwlub38UkWkEDeVr3T3PmjEeXfpUXL0TeS009jb7XmluSOcnT8drxa43r4qayNZX3ALAnvw165uwUdRUv3nvXwD43LDbsr6vcXXWfpwxACxPfOJf2/hEpLaQVOMA2JL7svbPlgtVRkcZr3M7Vuo1dtXc/bO+Y4ZEe+CPjJH9uxsFr7RXWTHrW8WY1x69YMvQsNrJ+fFznXLoyKy9VqLulWiIibbXZ7sYPBpIMhoPTj7reWVsLqmt5yGz2S5lXlnNXfhCXyrqTbSQppIiIoMJ6WFXAW5W1XkisgTwoaouKjpOpeWvb9cwH1XV/6th/0oD3QTc1Kbvx6b9AbBn2+McJ2VU5rQ4tOJ+YLSKwV1/w5fE3KNbdvmpasbTx9SZIrm/7oVJ4Rmtobm/klCO04HDgAGElaeNgXnAtYTsfoW/EFJu5UhECksmglBparPU2LFjdfLkydV37AJkiPnDaB12zT6LbZgR7RebjAi2Fvv0blOqzJmzPACDl4upUQ7u8/usfS+fA0q1iwHMz9qvJ8kb9yG/dvLrhPFfMq5Pr5uEj+NPSTwNDzQHxelz16ohwGYDnsr6UjtN2/YnCNdgY26mGq+zrbkbgIFm/rZ9BscD0J8FWd9/fr5OnMxJidZSkpq/embdjNzSze2DI0vtZyb31z3xN+lCpbkRkcdVdWxnxhi7jujkgiXJZQM6fb7OICJnErxkjwHuAl4ExqrqE4kH7SGqulHR8SrJ0mJhy05dyATPUxV3c7oQuaBzx9tCYb0BF0RVaOrH7BK+ARyvqv+XaC2W/wA15YGsZFN5uQOTc/Kw1Qy/ZJ6Uz0t+hLYY1idjDZVhSULJR5+Ij7arbRS9v84YdgIAN7Fz1jf+KzFO6flrPwnAh8ZieF+ivQCc/NKvAfjDqtE+N4gYs5LyRbMaOeoB47KexNjOWGX5rGu60WqGJ4571uPrC8YVfnGjVWzKo+32PYcjs/ZehEj/NDYG4H0GmuNDSv9rb9o367v2xPi5fOVXyTU8bC5sf3MTTCP67TLYjmbfK3Ncwo3WkyX6tLEt1v6y1WlNfdNt5rm1JK1lU/kE5dNxLU6N4rHQZYtIP2AjorvvTOAJz0bsOI6TQ2t5fz1HyFRye862nagxQXBFoSIifQillr5HkGap65ICb4vIb4DxqjWEOTtOL0U2p03si2sIPZbW0lR+BlwrIgOAvxLu7xuIyO7AwYQgyMJUMtT3AW4AvgTcTIgPmU74uD5JkGw7Ejy1dqnF5awRNNJQX45sTdoE5M29Ii7p/D6pHW+DBP/48+9k7a+e+AcgLh1BaaDiYrSX9cvzetbehpCwcujNMcpwxk5xKWtOYqjfeMJzcQBjTpyxbtj+rkl4Omp+ey26vy1BZH5oDy+3QdYemmTZsdc6xzgFpAb8N41Hl3VaSJ0N1p0RU8+tO+LJrP3cz2OdmIwNTDu1Zf3d9I1sf4gNTi1JKZN4h5e4FttAy/Q4FypNT10M9Z8WnXxLwfOt2FhDPYCIfBX4BWDzWb0CHKuqf6llrEqy9GBCAORuSXxIW34vIrsCVyf7/q6WE/cGMuP7n0xnnoLp9Fp6SsoTN9q3obU0FRLB8RcRWYPgiz8XeKEjgeWVLvtbwLllBEo6kb+JyHnJvi5UaiYJrjNPtw9eEVMbn3BBSN547KE/y/rWPfGxrJ0+te9C/BfNNFlu0kSP1iV5/dlTs/Zlw0OygjV3eiHre8sUjk/dfF/ZPxrPV7rmzay9yopB65G3o/bDnbG5IA3qfDH2GUWJzyw0rm7J+vPwHWJWniXvicrvkzuMBko1GZuS/6lE7bBp+M8nlmte6ntBE3pvgoldsQI+7d7P9H0pan0rj5gOwKtzhmd9i3ZfMu57VbLvdHMztVrL+KTfGu+dHk1zB1lEROTHhNpVr6rqVGCq2TYc+I6qFk4HXkmojAZOKDDGzZCs0ziO02FCjMuFZbd399N/ifZhveLm9QztqivRPvBhcxfpspwC3EL7PIsAKybb6yJUlPyMwG3xUkRlOG9uiAo8XEyOzM3ij3C0hv/hlGNi3wLjHtznK8FN9dfHxGpVb/8yupTs1zesq/3E1Be5fMJB8VyJKeTqX0U34ruGb521j3o8pL6/d8wmWd8Xd70jHn9s8m6DN38dm5IoLQu+Hvv6m337pxnxY9HGEk0lz/y35GRjmjPe8Wn6/X2fvzbr+9daUWt5MtFUbHCmdV++aWD4XP98aHQ5PsTcwGfPCP+jVCMBmDU11r1PXZp/veD4OKm8ImLnmfZ75oa7VRp8aW/OPac2iwuXUlTg4759Cu7dcHO0UD5/48pg6lUUoJJQmQJ8gepWgJ0oLaLlOE4XUjVhZRo/k5MxIByfCrMLy2yvvexxuWN6q7BRERYuVtSo8mGXziUPEfkmITEvBIHyOxF5p81uSwDrAbdRA5WuegJwtojco6p/z9tBRHYBDiWE9zttOCyprnxYWVNXSMcuJiDvm/Mvz9q/HPYDAH51Vqy1/n1+kbWveycELZ6y9I+yvjf2XyprDz0/eHVtTkytbwtmpV+VFcdErfe1v8U0RCsckbhtnWimvLxpJxpIf5v1zRYhSEu2GU2GD0zbPqgnX+cFpmhp/2gS4WufuSE0oiLH+Jkx0POuVUJKGJv6ZTvzPJQGR15tSkO8Oz96rW09Iuz77KL1sr7ha0QV69czgoZiPcoOHBGrZ39vakiPM3yjeMzsC0xSivRGH/+V8AkTKDkteT+b5qRcok2nLAv7NrVRZT6QGkgFeJv2JUQ+JJg3aso1UUmo/J7gTvw3EbkRmER0jBxJWLzYmbAW9/uc450EsaFDY8yPc7Pe+RTnNIZqBcVKctTVaczi4/Qs7zFFSlzemw1V/SshJgUR+T/gJ6r6UuWjilEpTcuixGX4ZOAogoBJn7mF8Gw5PplMwxcFW5qropfRewdG76RjjwtlgPWNaLbaZoebs3b/94MmMm/JeDdY+hKjSidZUFZbGGNHzlhgfC8SrWPUF03qlWNjM0t+aewoJQkV09I9NnLYKsq7J+9WU7Gajkk+mSal7D/e9NkH4vSbavNXG5vMNrODW9Wyw6N32vqXRk83VrkfgPk7xISV4wZOytqph9wGfaJHmk1ps/uI68PxJjXMTxdFW9YWa/wDgKfmb5j1pTYxgEVvJZ5iW8QpLbvFK1n7zaFJSYMi2ZKzNP1lbug5hvTcRJf2oeZ5ciinkVyYs/2QnO29F0X4uImFShuOIqsBUUri/fWuqr6Xtz2Piot+qvoxcKqI/BwYQ0jTIoQ0LZNVtfsXAx2nlyKbV9+n3TF10iQ6Sk+Jw6kVRUpy7jU5lxCWv76Ts+1UQmr+wnWrClmSEuHxEO5l3zHGlFmPTmIi7Dr9mlfEp+trLwieSjLSGGXM6uYPDw3j2iJdww6ON5HLjwieYP0XRNn/6sAYZ7HqVrNDw373bdKB1DxhtQOb5j715LrH9Nnc1lcm79bh/OembewjpE5nG5o+q/Vsm7yvbvqeNe1EKxs4LqbDL/EuSz4Wm30gL42+TYhpU/andikb87Nan2lZ+zeLjgbgxwOj5+U1A/fI2o/ul0iE6fEnt+AD88GnWuFDZYTAZjl17QcPQ3MLaVfQTiBqOmWi+2XzZHuu9hLGb3vevOWrWgz+PU3gdPfyl4jsBnyRsBZwvqrWYlzfkvJq6U3UGINY1udNRLYUkaXKba9wTK4a5TiO05tYSN9Cr3KIyGUi8rqIPNemf0cReUFEponI8QCqeoOqfgfYH4w3SjGWAfOEVcoHwOBaBqukqdwFbAbm0a4CSR7+uwgVw2rKatlKtF1OUC19gnzePOn+gF92x5ScXkze8lJHDO4lx61VcbeS83alhlHtt9bM1MmmMoEQ+ZS5hCb32fOBzwOzgMdEZKKqpmEdJyfba+HfBC0nT7vZmfJp8XOpJFQE2FlECnzFgApaT69n79IfniYV4dI18uemxWSH039m1o+SjCmzdo7G++8aTXRk4oy3LNE4/f3Z52bt354blr+O+u9FWd+qb8+O4/82ef+smdxM006N9meaPhtXm3of22+RjclN+23hse1M+4HYnD4lvI+0OTBXNO1/J+922cUutSV+K6MeN04H9vhkXoPWjfViRi6Izi7rPxWWHRffNC4VbkBclkxrtzzCplmfrf3ypz5hLfNuts76UjdmgFVGhA/2+v67kUsaSLluHJO/mw92mtn3ofw4k8rkrW7kxLZATKRZbimu4Hl72pJWLYTlr/K3179e9A5/vSgLC8mte62q94rIyDbdmwDTVPVFABG5CthVRKYAZxBqy9f6UH8ucKGIfEgQZLOB4YQ4lsOA79YyWDWbyo+rbHcofaKaUiD9WvZ06S7FTh3pbOXKDp2zjPNAZzSZRjsX1INgqF+87PZdDxrKrgcFWbKRTHmjhqFXovTRbxawKXAEodzfMiIySlULP3Go6sUiMoyQlsvGHH4AnKyqF9cwvy4pJ5yXP6bHUE0FHy2PZPtkaR5NreoS42X6JLhfHHP3gddn7T+eF5wxVt43VntcdmH0w5104Z5hTndEl+NR18bl1+v4MgDzlotutIMfeD9rv3RuMNqverrRXmJmE+Yk2sGwk2LfO8Yo/0oSyDjaPFzPNWnuh6QQUMlPAAAgAElEQVSBjDYH0iOmbarvjvxy0njTbLerB2m//cYaTSfLo2nT7NufVZLmfp0VoyY/bd2Vs/aoyUHD2XHFu3Pnl+bvH7RqDKi0yS3XICTl3HFaPH6PUTE99bDEU2DtYTH5xHNTo4baZ0w42aIzjUky9/mVaLT/jLlhjzLb09Oa71VJFcsvJe/7kk+asj3POQDg+eS8dg1js1PQJMa2IxH5+bS9L7bOQ5hCV7kU56XFUlU9Bzino4Oq6s9E5Fxgc0JipjeBh1T17cpHtqeSUNkc+Luqtq8v61RknollGMz8Dq9xO04rkSdMYl9OHMvRbZaFzyo+fvMvrVVe/gK4b9LbPDDpHYgLyUWYBSYVecjNVZcH+USA3Fx1xypUKtL1FuEZ8x7gemCiqrasFtJdRbqsQIH2QmXluXFxfJYMoi0rGxl+ZPLgcfyc07O+RdfGJ9lZh4ZHWesaO8dEF64592UALhhyQNZ32PcuzdrvJHkQlr7OTMDaP9Khppg+k28ys1kYZfyjGbHdLyli9YDRfizrGQ1m6XSs6PHMXKPVDEi8bwdYrWcH007dm8uFBmzVZj8oDfdKC3ZZl2XjUquJS/OcIfH3b28a7xIcJacmiS8BBpj/y7OsD8ALrJH1TTeLAbfNSPyrPyhjR7F+mGkYmv1fmaDKLHnl3kbTWMFsT9MCfcn0jTTtNCnm4abP/nTS4x8qo5GkNkSb+r8GoWKXv1Ktvzsj7utRpGutsUvqpZNHF9p3C3m87PkSm8rfVXXd5O/FCKnptyMU0XoM2EdV/5l3fFFEZH1Cld+xBEG1mao+ISLjgftVtbCwqWRcX44QE/0S8CNgpog8IiIniEixT6u38CeJrzbUbynAcXo2Iqe1+73k9VU6tuj+3UEdXIqvJMQGrikis0TkgCQg/XDCQvUU4C91ECg7EfJjrEDwNOtnNi8g2GsKUylNy0cEVehm4BAR2YxQQvgbwM9EZBqh3PD1qvpwuXF6G4P/9D7sF7S/7Mtt1p5nyZ8rHr+LKc27fZIQcfByP8z6PnHo9Ky90nPB0PDO6GgQvLVvTKj44ZC7AThsjaid3DV1s6y9zegkltUmebTt9InVPkOZgMWPEsfFj026lQHGvpIq5TYzy+p2LOPJ9VFiM5luCnqNNIsC/dK2Geyjm3K2lwuOTEv/Wq3LlhNO7TN2IcKcSxKtboVtzRLzAtMcGy5grX+/nPU9/+lPZu1BiXoxzGTRHGoMSCuOCB+WTQ3zH9bJ2r9ZI0aQfm9G+OC33jnad+6+1xhNUluI1V6s1pNqKDb6wbZTDcXaTLYw/+Tjk9vG82W0hXT8K/M3Zw4qZ5cJzuTCmoMnm416BD+q6t5l+m+iNCS5s5wOTFDV7ySakP3HPkWNGUQLF7xU1TSi/vjEzXh3gpA5TkReB/6mqi2TvrSnJbBznN5EswscRVhQJU3Lw5Pe4JFJb0JtNpWuYC1i/uy29pB3KK2oVJUOVVFW1ecJ0u30JOHYbuSWXOol7Jdvl0oF1cFZQAhcJIe02WdYyRry+VfFxW85JxnXrH2/sWz0WFrx4PBY/3D2GF6ahmS9uUnKF/NMs83T7TPtPPPl2F7/frNhQvJuvaBMeFS/PZN3a5vIKchls308a9bmlzb96TO9DUN5xSgF6ydzeMfYbIaZ8/ZLtQo7V+u8lGpARrvAVo9Iz2UfmO1PfbGc7UYT6J8uCJs5r/XfqLU8s1y4waxn1CebnDIt3dyXGKiz+RqxZMElJj/OV0cEbfeORdtnfcO3NCn3T05sNdPNXK132Fu0x9pc0v/RJ0zfe+ZWkV63TS5qSTWdx09p3wdRa3yovoknm+nBsIimsvG4YWw8bhi3XDy7Zg+rOvM6JSXxSlgHmFFmWy6dDlhU1dmq+jtV3amzYzmO4/QEUqHSGZtKN3IV8BMRsQumKiJrAD8EKq/Zt6GspiIi93Zsfhygqv+uvlvrIG1KkFVyfax1WS3s3zKrho7TLfSEAMgWSn3/I2Btgqfva0nf3wj6622UpoGtSqXlry2AJyldJKiEAJ8D2vvJNiHdqSpfJHa94cKcsq2x715ivfhZXwsuw0PnR4PuO/2jUX72vWGZY8SU6Cfb9+C4fCKp8dnWSPmjaSdLRusfZfpsubV0xcW6rtpl4sSV+M/GZXhfa4hPlOZxxuB9lnHp3d98+z5K7MA24fG+ZnlrTrKsZW81c8xS16bp8phZfptmnA42Sve1Hv7WvSS1Q4/O6YPoPm0XKuxnkToA2PuI+VzWXzUsRb7y2ejJYNO4rDk/bO9vnPavGhVTuthls8XTNTyzznDNnJgROQukXCF+gCtvZFzZ97FrYckxP4wf5jbD7gZK3dPfMjkFZ60QjreVMWeLcaU/P1l3vCV2layBpkG/tnbMWjnbOSS6FLdZ1de5zW1XqZamBWDypNd4fNIcaLBNRVUXAF8Ske0IrspDCSvRd6jqP2odr5pN5buqWjSh5GJ0otiyiAwBriZ4zE8Hvqqq83L2W0j065mhqt1uy7FJI9diel2eqtIfyD0F0rw4Tm+hXOBws2syRWwqG45biQ3HrcQdF89otE0FAFW9g9JItA5RSahcSmmoWDUWJsf8t4NzOZ4gGc9I0jkfT1jPa8v7qrpBTn+XYZe7ni+JEsvjkHapXC7Vqfm75gSJ/cSkWzuOXwFw9cCYyfr/7oy53Q7aNjgAPLllfLze8DbjM5sYUp+ZGgPu1r8zzuXjJLnjYuVSn6TmuT1Nn0nZkgbh7WtNfDlJHK8zT98HmuDFB4wmsVPyAL+RSdPyjtFEhiTHzTXHrGRONTdnwdVeVuqy/IpxWV7a7JDFNNp/3ZOmnWpg9pnSuienxm2r6Yww7eSXtOJcc4HmhpnWuVn12ZgyZ+CoGDw5Lib9yZJa2uDJPYZdk7X/MvWbAKy2UQxfmDPfXFjq52M0tcHLRW363WSxYdUSS39sD1ojBOhOeSIWvxluKtHOThxDlto6prR6bw+Tc+bCnPQvr5l26ktik1xW0WSajeD9VT73V6MRkYHV94qoarnU+O2oFKeSVwWs0kmV/MphRdkVshSvfwDuJl+otAxrJ3el46rs5zhOz6LI8leDeY/27sOVKGwgKpumpbsRkbdU9RPm73mq2q44jIh8TFjl/xg4Q1VvqDDmQcBBACNGjBjz8ssvl9u1UwyaX5pk9N2B4alsbfOo+6+SkoaljDXqweNXxDz0ukKI0JdL4v9I94xR+5vuHiwQjzy9VdaHeWL/37iw6L7k64tiZywSGQ0Y1s3WphFNfxPWodCmyU80iekmHf5I4578fmJTeMYmjjSHW2Pd6uuF97nm6f9Na9NIKHE5br+5ZEwbh7lSzna7srJycq1Wa1vaal2pLcYOaueXmh+spmfTyKRjWa3OpKRJNaBXTOr794mJQN81pso3kkk8Zb5TM42xKF12GW5SQtl90+y5c4xaNspERw7g/Xbnf9Icn85lG+7K+u5im6zdP/lC2ZIMj84w6Yx/k3zI1qbVJs1+cLUvbzNp9jQtI8cuqydN/mLFfZ6eNJNnJs3ivounTVPV1SvuXGdEZH9qECqq+oei+1by/tqn6CDJSa+oto+I3E6pR3zKSTl95Rihqq+KyKeAO0XkWVXNLSKjqhcBF0HI/VXDOTqFJOr/aBN9/kuOBOAHYi81+OmP0e1xHCfS7DaTahSxqaw7biTrjhvJfRdP63abiqpO6KqxK+lnf6qwrS1KSYL3Mjtp+buniMwRkeGqOjsJqMy156RJLVX1RRG5m1DVvKbKZPXmvSVNZakb49PflH3C090vrjiy4vGPb24e/42XjOwY5OBqV8S18R9zYtY+Lcmm8M9Pr5b13ffpmObjkHMmADDvSJP6fsWY+v7Cv+0f9vvqhKzv3t9F77PU42jwKfEY7ozNVxJN5CN7MeZ3lHpnbWq8vyaZ/6r98q2eeFI9Yp7+R5rtaUISa7KwPkypDrqR6bNfoNTPqtwKfL/kvEub888x9pdhyUP/AOPx1c9OJu236o9NfZJ+RWxCzFfbt1fa29hcbMoco8G98/WgaSzWN37vtjZaw1tJ1OJMY9QZRExUmqaHsan7rXdZmkrmGqJH2WGmmOALSdJMazOwWkt6M711UUwZVEJqXrEeYUeb/8zZ0WYSqW+gZFfT3TXq64GIrA2MIejll6nqayIyCphTS7b6SkKlu5NGTiRUGjsjef9b2x1EZDAwX1UXiMhQwmLML7p1lo7jOFUokqalWRCRpYDLgD0Iz4iLERzCXyPEqMygBtNwpdT3I3I3lEFVawrlzznfssBfCD4zM4A9VXWuiIwFDlHVA0Vkc0IkxSKCl/5vVPXSsoMauiv1fYkLZOqx8hnTd7ZZJzapv/PSfqeJAfXUaEdZZYdoNEmfhJYl2nSevSZqGlftEeIcvvaSMTvZIlnJKu6fx3wl69p3n2uz9tzE/jLElhO2/+XUwGE0kX//2gyfLP+9b4JPrjP2lX1NaeGbE0dGq0nYBZBUgbAuKDbNS6pLWTvLSNNObSn2GJvoMtXljE5W0k7HsulbBxmtI9Vg+pWLOEhS52PLJdsJpHYtm/LG2rdsSpTEPqNlzjVzSBjYVh58yXwaHyY3O6tpbGiCkVKbiU0jM90cnx73oblpTmZMu/GtHejaK2JFsD7bJXE0u5uLtWn00/iVea1rU1l57DA9avJe1XcEfiDndvp8nUFELiKkiv06wSr4ATA2SX2/P3Bcmnq/CJU0lel0kXdAHqr6JqUVzNP+ySRheKr6ILBeZ87TrHghL8cpTmfKFXcX1Za/pkx6kSmTXoTGJ5T8MnCUqt4lIm0n/TIxNV8hKgmVg6hNqDgA8+Lz9egHw0K59efn7NjUeYl319/Mxzx4GJpoAJnWYopgPb8gFoH6dP9nAHjT1p01/vzjFoTYhtdWjd/ZFZaINsHXhof+fWdG7cRG3w9JHQ3sAqP9yqVGi/+LXasb54TUO2yAeTrf2kSZn2vCrI5IntPmmuj9UTnfThvbYqNyv5q8W+0mLxXELNO2YSYb5WzPw7rovPtBzg7G5NrPeoqlST2tR1mev49dHygXhpaEhIjVasxtacR6yT/GzG/UYvHK5g0Petkn5kZd7NUhcbJrvxNimd5aOlYGG2kyhabxMa8bVcvGtKQ2F5vmf+V9TET/FYk1zNpUBp8Sv/dpWqQJRmDMs//Z5revKFI1Tcvq41Zn9XGr89jFzzU6+HEApYW8LYMo1a+rUilO5ZJaBnIcx3ECLRCnYnmMUCfrlpxtewAP5vSXpZJL8UyKayqqqjWpSE55Xsk0j04VdHOcliKLSzm6eZe0aqGFvL9OBm5PQj7+Srjv7ywiRxOEypa1DFZJlN6BL3/VTGnaiKRt/VxL1oCTtsmCPOsss5SVIAvMv2GJuAzwU/0lAGfOPz7uu3jc9+b+ISBtm/l3x+ONn+0dw4OH99arxO0PrhKD1Pb89SQAnn8mPi/YaoVDn0/ytMTMMnxk662kDgomzctKJjjzCBv8l/gaDjHRjVOsS2/yTbVuxNYMlbo198vpg2igX7rM9nSpzB5vx08XiuwijM0MMzpZano9L4klMChdqjIJL0vq0Cxs8w6lS40TTTtNUmTr1tvPMnVVtuObX/rgJZOrMY4AKz3VfvVj6JLvxT+MU8GgFcOy6+v943fdGuXTZZ/HTcnQWU+Y/9z05P040MRDviTQ0VSEjDXqm3/Jy6JIiaNEM6Oq9yfJJM8g1HsV4DRCeOr2qvpYLeNVWv7avxPz7BXY8NDqoZ8Jp5g69ukNfoDLbqd30lMdVIrYVJoJVX0A+JyIDAAGA2/Vku/L0jKLfj0ZfV/gd0GwHGaCGw/VP3M+x5T88I7Vi7P2t5LSjCe/Hf141So6yVP/I0Ni/s311noma6fBbZubJdM9J0zK2jOuCI+naz1t0tvYB8b06dUYn/vZ7enDr7WIHxabU4xTwOjkqXuucVm2WsWQxKHxCWPIn262H5jM4X0TULi+eVK/JNEarNKYZ5S3emaeJmON/yNNO1Ww7A/qHRNIOevt9uP3M+7daXqYpW3Epg1zsDffVEWyPkM2+WX6/7CainUaSHkupw/iB28na9ynBy4T0v4sMGUYVjPxx4snycon2WKwNmFk+nXcz560XKBj0Oatl1czp7xPaTGbSoaqvg+8LyKDkyJdU5LU+IWpetVJwOGBBE/7NMHQTMLy2KV56emdwIpmrePVkqCD4vTUJznH6elUs6m8NOlfvDTpX9Bgl2IJUrq/qh6f/L0tIfh8IPCaiOygqoUNvBUTSiaDX0VIrDCH8HAoBL/lYYRn0b1V9fZyYzQL3RX8aMkTKiVPWS+Y4MckO33J9iQITK8xS2bm6VXGJf+7p+Ij8Xkj4hNfqoncfcyOWd9qZ8Xvxl5c3W7Otm/944Nr6ZNnmNT6M43akShIU34bu0bbNPnpk7JxiS4pGHalaafBgTY5kHW0TFUBI2Q/MkaNLOjQaAfTjdaS2kqsdjLStNOgyemmb0BO22pPNtByUJv92tIvZ7t9ohtm07ek+9qUMFZrSfe1mohtp/ta7cSebPk2+0HpbS19/rHHmGIT/1s9JCqd2T8msbyDmIEpTc3/x3tN0nKrfKQ+Rm2CG5sl5qQewY/Lj11F95x8dKF9L5BjGx38OA34WZoPTESeJOiWpwHjgf/VUreqkvfXGgTz4NPAF9saa0RkE+A3wA0iMlZVn88ZxnEcp9fRYjaVFYEXAURkFeDTwMGq+qiInEVJJFp1Ki1/nUwIT9pGVdtVdExOuA3Bh/lEgp+zY8hb8qr6NDa4fcCX7GG0SVusaC3QB+EBkwfmE8RiS4c9FzLYyIR4/H9GrZO1Nz00LOrfSIxYTAPXANb/X9BUbjdPoRv+L2oq/zszPLEu/RtToeB840WUBvzZp2hbmnhT07Yp41NMGn3S2mIm5U0/kxLmo8QWMd1oN8vaIlxJUOG/jUeZNfVsmsxxkPHYysugZzUN60OfeoeV8z5LsTYZO9acxGvMZtNf2niSDcgpvTzEaBclmkxqCzHXUpLIMv2MViyzPa0TZksvGz7sH052n/E0tanx/zwn8WCx3mlX5ttBmkU7qTfB+6s1cn8Rvurpt2lbYJ6p+PsBUFNBr0pCZVvg+DyBkpIkdvwNwRXNob0Rsat/NCKncb87jzlOU9FiWYrvAY4XkUWExJE2me8alC5gV6WSUFmO0iXmcrxEvm+J0xFK0lEkfvqTok1lxXHxUXu2XA7AFnJjPMCk3k+1AjXHyxanZu1rDg02l8sfPyjrO2vMoVn7wHPPC30LYyDNLmv9JWtPvCYkR1ly+aidnHVYPH6Xw0JwxajjjCXjfnNxNtFmus5vE15aRS/1CrN+KGadv1/iybS6zR5nftPvJ1rT1nbl2niavZ881a9u6p29b0xwc+xTf4LVNEYnz3m2MkZemn2rqeTFzFhNp9yPs1+y4XVzriFGGxyQ/jus9mEH69/mvW07/dzN5ztvrXi1qa3uEaNqPkJMZLrozGQyZ5d+l3sbLSRUjgb+SLCfP0VpfatvAPfWMlglofJfgi3z/gr7QPgKvlFlnx5NNRfHzmgvInPQidX3cxyneWglm4qqvkJ0lWnLFyit7FOVahH1x4jIX8otgYlIf4KUK5f6rlfSmSWvhzR6am0mwU1GdnnNROq/lHNUdK1ZausYCLLXzsGTaz2bevH8jbPmcan9xHhRfX/lmD3ylmFhLhv1fSLru4f4KP+PPbYASj1/bLr085OglLP/ekI8QVR0SqtoJd/EN8bHhfihp5iI7jTRotVe7KP+Mm3eocRoMsDaZ1KMd9iAVGuyT/zmSX1kOlerk1vtJdEaRpvt75sH9QGJ1rC8MZrMNZpGajOx5YzfNz/lxXJ+qSvZa81LLmm1D5uoMv3cbEJLW+Y4iQl6eK34AbxrDCSppnL7nKgWLhptJpDa/TYbhmYhUD3TdlKOVo1TaYuq5uVlrUilqx4PPAHcISJHqeoTdqOIjAF+S1hz27vWEzuO4/RUiqRpeW3SZOZMehwan/q+rlRK0zJVRHYjRBM8JiKvEW0sIwm15t8CvqyqL3TxPHskNrBR55bfD+C1pERsaei04zjNSJHlr6HjNmXouE2ZcfEdjU59X1cq6meqeruIrAl8hxhRL4T0uecAlyTFtZwK1LIctpnsFv/4c1jymr3PJ7Ku1DgPwI1h3Jt33saMcDc7JvXCD0us2rMXRd/RD/Y3gZTzof9Aha1i38JHjR9o4ub7ny2jG/JbO8e5rJ2sL33+0mh2093j4ZKWNTd50R7bNBaQ2/icmCdk3pHBEGxdmod+5vF4YOqfYFOL2CSKW+Vst4GYqf+KzVBgjPKZ+7O9D9iC2nnPkrZyY7o8Zpa3BthzJf39TMLIEkP+x23egaUX5GyHuNRVbvkrnZf9ddtlw9RZwSx5vfbpONh9OUlpLyA6YNw9NVmiXdOs7w22VRytgT7PXaF30BOWvzpC1atW1bnAmcnLaQXOSYTEkb+uvB/wS47k+108HcfpbbSYS3Fd6Z2itM501DCft+RVmjgvPPENPzwGNP5J98ja+yXbd+KqaMg/J2odFwxJXIGj7Z/Drjgva69JWLWUFWO2vx/pOVn7p6ePB+ArZ/w565tpIuLWeSRJIrhDHH/jIVFrmbxVMORb7WH6Z+Mj88bbxg2DXwrhg5+9JGonM8ZHVWDFd4Iq8OzSJmXMA8YSf2d4W2DS8Pe3iwqp+3COazAA+yfv9hibWiQ1apsiniU+ManR32oX1hEh/djs+e2+6XmtRmKxWklqgLfBhXmainE0+Ni4Us9fMqz1P9431pW3Aa6TE1XmwfmxDMJ7N5hMpfummoj5gNpkAOypQY1FcaGSg4jcCRxaNP2KiPQBbieE9/+72v5OffkGIXvx5VX2c5zuIM/NvrcJmmZ2KRaR/1JDvSxVXb76XoFKmsrWYCrvVEc6cIxTgVigKP5A95P12u33dRPIIqeWfk90LvBfY0f5L7Bc2EfSB1HjZjF+Tky9P+HBvQDYf/OYZPKMB2NBsJc2DUaNczgy65v89BZxLrPCefrsFh/PX7B5Zox39MeJfeOu8fF4m+ZixMTgtb7hF6J28sxn18jacz4bvvODiG7Iaw6MFzb4f7Eee3ZOYzNZLLUMWk3CaiqpUlXuF5PaJ8p59KeaiA1ItFrJ3Jy+PO3DzmF0/r7/WyWkz3mqf1Sr5ptQzdQF/F+snfW9anK2PP7zz4aGjT6zgYybpXaSUzKX4VZIR9+dLKJPs6dpOZ8uKsJYbfnrBhGpJZe+JwxpMkTmoK9X389xOoILk/I08/KXqp7aVWNXEip/6OCYvTq6vssxCSePmBvKCZ+7jzG17xirUKY2GVn+1JIh0moHLz4Yn05XTYwCq2TZBGH/fYKGsvWDt2R9q4vJXZKmhDHFln5zQ9Rk7t8nrM3bJ+LVn4hphC4dt2/WPmBqsNvogqhVXbZcLK25zX43A3CaCaLbcmYM6py4ykgA3iImt1zTqGATRwXDz/BRMeJxtnk6X2/pZwFY9fl4/dPWWjlrr/KpkGpm/sD4xD94dtR+FiQ2j8WNpiJ5GZOWzumD+Es0mswCo328OjC6uqV2rWVNaee72Tprv0Gwf9iARXutty76AgBv3r9SPIGtCXBSjs3E8lBYyuruPHethNtUclDVb3XnRHoysnnp3zHKuHlIBdDKrms6TqdRmt6m8pfqe0VU9atF93XvrxagXCnVc7OH+jaJ+64YlhyX2mQOIY9PHZM8tZv18lnHjIo7rAB6FvDITlnXn/UrWfsSDgTg7vOie9nft4xeROk6/THzz8r6VtsoFgl7lmgfGrNGCIrZlHuyvm8Yt4M5SbzDRMZlfYetEj3ZLuS7AFxPjPN5kCjN+yeuVsNNbpbFiEEj6Q3g/LViPMamJr3NIwND8sS0VC7AisPjWNOTkl/LDozaQ98h0UCS2noGmYT6M3Nyy79FjAN60+SEedeYKtPPIj0nwH2LPpe1F+8T5jj7JhOcYksOpKc4rPzSVdBE2vdV+tuxNH2aluW6auCmvmrHcZxWpMjy13uT7uZ/k+6FBqRpUdVtqu/VMSqWE+5JNKKc8DTzJDoqpySBRJNBZgepRslSmolczrQSE8X9u7lhBfO7Q0zhthpKuKZLYiVr6+PN/ifFsdJxDiNqJVmczM/MoIcZreroGG3d54fB7WrR28aQYE671M+Cqe69N+KTPNPjM9EmW4bs3JPnxNiLRa/EsQ7aKNQ8fsoEmliby5NJUMeW3Jf1zTHR4KOYBsCCMh49qSYxkPlZn23PS9QDqx1ZrWQaQUMcZrTOu+bH3/0qA+P3Z+b88L2y3m2P3/TZOBmrlaSc1P7/7jaRfOpRTrj/2HV1+ORrCu37soxuaDnheuOaiuM4Tp1ptSzFIjII2JWQIHiJtttV9QdFx2qdq+4l5D09lvRt1vWVJCP5tpjc4xIN6dAqiTEdp7fQKt5fIrIaIdPfQELE038JWfIWI+RKeBsoLFR8+avJqCpUcmrYh/3aB0qmx3f8/FGolBu/3ZyOK92kJ7ZdsisjINO8DSY2soQvhbc+B8ToxEV3mKWyLYJRfPiIuEw0+4loqB6+UYi0nH2FMV6bBYfha4TtHy6K6crfnBWXv5ZdOXzWn+gTU+ZYo32a5nzO/HjMyIExuvO5e5M6NqZETJ8xOdcShy9t5wUi7l2arDG6kveuevD1ph7LX/3Grq9DJv+90L6vyycbuvwlIhOBPoQUrP8j/DKeBvYCTgf2UNXHio5XSFMRkW9U2LyIIMmeVNVZFfardo49gVMJccKbqGquBBCRHQl1XPoSsiSf0dFzNhr74++qH3x6Dr+h9Hw8ELF5UISFi1pDUwE2AQ4kZqNbXFUXAleIyFDC/Xbzcge3pZCmIiKLiNHy1tHQ9i0Crga+Va5SZJVzjE7G+D1wXJ5QEZG+wFTg84RwrceAvVX1X9XGb8h1YQ4AABdbSURBVEZNpahQscZ3qF57pdbxO8umpoT1o0NK06brXJCppftrkl3FOipwZTrX/CW30Rrcd6dsbjI6mqDLzD32aHOdE8z2/ZN3+zxoNYEcxhz6QNbODOHljpmevBs/AuPxG7efnDMngDSn5y2mr4xTRd7/1Y3u9aMemkrfjTbQgfffXmjf95ZcrtGaylvALqp6r4i8ARykqtcl27YFJqnqkhUHMfQpuN9ngZeB8whVKNZK3i8g5H/9InACsDtB26gZVZ1SoNjXJsA0VX0xEVxXEYxLjuM4TYOqsPDjvoVeTcBU4JNJ+0ngEBFZQkT6AQdQUni7OkU1lWuBF1T1xJxtPwdGq+ruIvJTYF9V/VS7QYpOSORuymsqewA7quqByd9fBzZV1cPLjHUQcBDAiBEjxrz88ssdnZbTCaJrciBL01/TGJW1riJLP5XcaDO7z0NVxrF2oM+Y/ofNmGmSxRJNrLr9q+2cnMZQD02lzwYbar8776m+I/Dhsss0WlM5BlhJVY8Vkc8AtwIDCCtHiwH7q+qfio5X1Pvr8wStJI87gfSmfi8VvARE5HZCGeK2nKSqfyswD8npKysVVfUi4CIIy18Fxm9Kqi1t1LL0UV9DfmNugHEOxb3Tmom2aXvK7ldBUBb9H7uAagyqwscfNYUWUhVVPcu0HxaRdQlVmAYAd6rqc2UPzqGoUPkQGAPckbNtTLIdwnJauTJIqOr25bYVZBaU5LZYmRpVM6f76Yhm0n6MfNtB5JBC5+nz2nH5Gx7KsVNUsWVVFbBX5sy1xCW8VCjW43NymgVh0cLWjNhQ1ZmQFGjqAEWv+q/AaSKyELiGUNNueYIL2qnAZcl+G1BSnaPuPAasLiKrAq8AXwP2qXyI01to6+nWUQeHIudwnIoo0Bz2km6nqFA5hlB86xfJy3IFcGzSfg54qCMTEZHdgXMJic5uFJGnVPULIrIiwXV4Z1X9WEQOJ6z59QUuU9V/Vhi2R1BtCaOWJY7OLoc0fjnFPN0PNk/284rd7BcOK+zEUhch1I6HTjOfYWXNpMNlqhv+P3JYJPBBa2oqnaXQVavq+8B+IvITgnlyBWA28IiqTjX73djRiajq9cD1Of2vAjubv28CburoeRzHcbqFj6vv0hOpSZQmAmRq1R2dpqIW4/pebWqzXc03u2ROXUl2vYPzHRpa8Um+GRwknBoIBVWaFhEZAcxW1Y/qPXZhoSIiA4FvE+JThgBvAncDE1R1foVDHadulDdm5xjy551WwMCfjlvsJl2utk1nxnR6IE0uVICXgM2AR0XkTuBQVX2+yjGFKJqmZQWCAFmDEAT5GvApYA/gCBHZWlXnlB/BcepPVzy9N0oj6KwbsGsyTYYCddcByiMinwJOApZR1T0KHPI+IYEkwNaUL3RdM0U1lV8Ag4HPqWqWu0JENgeuBc6kNOmE06K04nKXpdzNtCtust1543Yh0WIoMZNWBxGRywipVF9X1XVNf7v8h6r6InCAiBQr4hIi538rIv9I/j5CRGaX2VdV9YdF511UqOwE/NAKlORMD4rIyUDLJnVsBvwpszhd5dLbTK7CtWQP8O9Kk1Kf5a8JhNRYWV3tJP/h+Zj8hyIysUj+wzZ8B/glIc2VAttRXgwqUHehshTlgwxnJdudJqUn33haXQNpe86OCree/D9uSeogVJIEjyPbdGf5DwFEJM1/WJNQSewn45IxFgG7qeqjnZtxoKhQeQH4OqU5VFP2I1bDcJympCNP965BOh2mmlC5+SK45aL0r6E1jLwSlNQmnwVsKiLLAuOBDUXkBFU9vYYxVyWEiNSFokLlV8DlIjKMEOw4mxCr8jVge4LAcTqI36yK01WfVbP8D4rMo1nm6lSgmlD5/EHhBbCrvFFhz7bk5j9U1TfpYDI8VX1ZRBYTkb2ALQjevXOB+4DrVLUmnato8OOfEpfinwCXmE1zgENU05pzjuM4DtBVLsV1z38oIssDtwHrEyr/zCG4Gx8GPC0iO6jqf4uOVzhORVUvEpFLgDWJkuwFVV1UfPqO0xg68nTvGoHTYRYBH1TZ54lJ4QXL1DByV+Q/PAtYllBGJCsbLCIbE7x7z6KG1ahaI+oXAVNsn4hsD5ylquvXMpbjdDVuE3EaRhFD/frjwuuui9/O2ywiVxJiSIaKyCzgFFW9tAvyH+4MHN62Dr2qPiYiJxByMhamHhnPlgHWqcM4juM4PYP6eH/tXaa/3vkP+wPvltn2LrB4LYP1zjSajuM4XUkRofL0JHim5uWvruBh4IcicqeqZvWwRGRJQnzKw2WPzMGFiuM0CF+e6+FUEyrrjAuv+/OXv7qRY4G7gJkichvBUL888AWCt9nWtQzmQsXpsfhN2mkYzZ9QMkNVnxKR1YHjgI0JXmCzgQsJ9vJaXJ7LC5UkQVkR8mrOO05D8Kd/pylYREjZ2CIkguP4eoxVSVOZRpC31ZCC+zmOY3CB14NRYGGVfaZMguebwqZSVyoJlW912ywcx3F6GtWWv1YfF16TG25TqStlhYqq/qHcNsdpVvzp32kKWsimUm/cUO84jlNvXKg4juM4daNImpYmQET6E7y+/q6qT9djTBcqjuM4XUE1TeWlSeHVQEO9qi4QkZOA++s1Zp96DeQ4juMkpMtflV6rjIMtLwJotKH+EWBMvQZzTcVxHKfeKPBRoydRmB8AV4jIh4ScYnNoEyaiqvOLDuZCxXEcp94UiVNpHh5J3s8Bfltmn75FB3Oh4jiOU29ay/vr29QxgN2FiuM4Tr1RqqdpeXUSzG58RL2qTqjneG6odxzHqTfp8lel17BxsEFTGOoBEJG1ReTrInKiiKyQ9I0SkUG1jOOaiuM4Tr1poeUvEVkKuAzYg+BesBhwC/Aa8HNgBiGWpRBNo6mIyJ4i8k8RWSQiYyvsN11EnhWRp0RkcnfO0XEcpxBFXIrTV+M5C9gc2A4YREgSnHITsGMtgzWTpvIc8GXg9wX23abWHP+O4zjdRmu5FH8ZOEpV7xKRtl5eLwOfrGWwphEqqjoFQESq7eo4jtP8tI5L8QDgzTLbBlHjlTTN8lcNKHCbiDwuIgc1ejKO4zjtSHN/VXq9Ogn+eRA0vp7KY8A3ymzbA3iwlsG6VVMRkdvJrxR5kqr+reAwn1XVV0VkeeAfIvK8qt5b5nwHAQcBjBgxokNzdhzHqZkiy1+DxoXX6w2vp3IycHtyf/4rYfY7i8jRBKGyZS2DdatQUdXt6zDGq8n76yJyPbAJkCtUVPUi4CKAsWPHenVKx3G6hxaKqFfV+0VkO+AM4DyCof404GFge1V9rJbxmsamUgQRWRLoo6rvJu0dgJ80eFqO4zjtaQ7PrkKo6gPA50RkADAYeKuWfF+WprGpiMjuIjIL2Ay4UURuTfpXFJGbkt2GAfeLyNPAo8CNqnpLY2bsOI5ThtZyKbZ8QFi4q5YPoCxNo6mo6vXA9Tn9rwI7J+0XgU9389Qcx3Fqo0WKdKWIyM4E28oYglz4WEQeB8ar6o21jNU0morjOE6PoYU0FRE5GJgEvAccBeyZvL8HTEy2F6ZpNBXHcZweRRMIjIKcCFykqt9t03+hiFwInESxoHTAhYrjOE79KeJSvGASfNj4LMXAssB1ZbZdC+xXy2C+/OU4jlNvimQpXmwcDGyKLMV3AVuV2bYVZUI2yuGaiuM4Tr1p8izFIrK2+fMc4BIRWRa4AXgdWB7YHdgJOLCWsV2oOI7j1JtFdMIpt1t4jtJqjwIcnLyU0kzFt+DlhB3HcRpMc0fUb9NVA7tQcRzH6QqaODGUqt7TVWO7UHEcx3EAEJHFgMXb9teSssW9vxzHcXoxIrKMiFwgIrMJeQDezXkVxjUVx3Gc3s0EguvwxcA04MPODOZCxXEcp+40v/uXYTvgYFW9sh6DuVBxHMepOy1VpH4G0KE093m4TcVxHKfutFBGSfgBcLKI1KU8rmsqjuM4daeIpnJb8mps7i9VvUlEtgemich04K2cfTYpOp4LFcdxnLpTRKhsk7z+1NDcXyLyK+B7wGO4od5xHKcZUVrIUH8gcJKqnl6PwVyoOI7j1J0mzyhZynzg8XoN5oZ6x3GcupMufxV5NZzfAgeJiFTdswCuqTiO49SdltJUhgKbAi+IyN20N9Srqv6w6GAuVBzHcepOS8Wp7EGQgP2Az+dsV8CFiuM4TuNoHU1FVVet53guVBzHcepOS6VpqSsuVBzHcepO6yx/icih1fZR1QuKjudCxXEcp0tojeUv4LwK29JSY4WFirsUO47j1J3WcSlW1T5tX8AQYG/gaWDtWsZzTcVxHKfutM7yVx6q+hZwtYgsA/we2LrosS5UHMdx6k7reH9V4SVgbC0HuFBxHMepO63v/SUiw4FjCYKlMC5UHMdx6k7rLH+JyH+JBvmUxYFBhJr1X65lvKYRKiLyS2AcIe3yf4BvJet6bffbkZCrpi9wiaqe0a0TdRzHqUr3Ln+JyJIED60PgbtV9c81HH4+7YXKB8As4BZVfbOWuTSNUAH+AZygqh+LyJnACbRJDSAifQkfwOcJF/yYiExU1X91+2wdx3HK0nlNRUQuA74EvK6q65r+vAfrLwPXqOokEbkaKCxUVPXUTk20DU3jUqyqt6lqKtofBlbO2W0TYJqqvqiqHwJXAbt21xwdx3GKUZdywhOAHW2HebDeieDqu7eIrE24X85MdltYp4voEM2kqVi+DVyd078S8YODoK1s2i0zchzHKUznDfWqeq+IjGzTnT1YA4hI+mA9iyBYnqKAsiAid9Y2Fd2u6M7dKlRE5HZghZxNJ6nq35J9TiKI7zz1LS/ff9u1QHu+g4CDkj/fE5EXaptxXRgKvNGA83YVPel6etK1gF9Pvfhk54eYfSucOrTCDkOB5ZL2EjUMXO7B+hzgPBH5IjCpwDhF7CTDgc2pcI/No1uFiqpuX2m7iHyTsIa4narmXcgsYBXz98rAqxXOdxFwUQemWjdEZLKq1uTn3cz0pOvpSdcCfj3NhKruWH2vDpH7YK2q/wO+VXQQVd2z7AlERhDs2V8iCPWza5lg0yx/JcanHwJbqer8Mrs9BqwuIqsCrwBfA/bppik6juM0mpoerGtBREYRHKT2A15P2r9X1ZrW8ZrGUE9IajYI+IeIPCUiFwKIyIoichNAYsg/HLgVmAL8RVX/2agJO47jdDPZg7WILE54sJ7YmQFFZB0RuYJwT90GOApYTVV/U6tAgSbSVFR1VJn+V4Gdzd83ATd117zqQEOX37qAnnQ9PelawK+nRyEiVxJybg0VkVnAKap6qYikD9Z9gcs6+mAtImOAkwiG/qnAgcCfVLVT3mOSb7pwHMdxeioicjOwA/AM8HNV/Wvdxnah4jiO07sQkUVJcy7B/7kiqrp80bGbZvmrpyAiewKnAqOBTVR1cpn9pgPvEgKVPm5WL5carqfp0+eIyBBC/NNIYDrwVVWdl7PfQuDZ5M8ZqrpLd82xCNU+axHpD1wOjCG4ju6lqtO7e55FKHAt+wO/JDjmAJynqpd06yR7Jqd11cCuqdQZERlNkPy/B46rIlTGqmpTxxUUuZ4kyncqJn0OsHezpc8RkV8Ac1X1DBE5Hhisqj/M2e89VV2q+2dYnSKfdVIedn1VPUREvgbsrqp7NWTCFSh4LfsTfieHN2SSTs00k/dXj0BVp6hqI4Isu4SC19Mq6XN2Bf6QtP8A7NbAuXSUIp+1vc5rgO1EJC++odG0yvfGqQEXKo1DgdtE5PEk8r+VyYvyXalBc6nEMFWdDZC8l1snXkJEJovIwyLSbIKnyGed7ZO44b8NLNsts6uNot+br4jIMyJyjYiskrPdaSLcptIBiqSbKcBnVfVVEVmeEJvzvKreW79ZFqcO11NT+pyupNK11DDMiOR/8yngThF5VlX/U58Zdpoin3XT/D+qUGSek4ArVXWBiBxC0MC27fKZOR3GhUoHqJZupuAYrybvr4vI9YSlgIYIlTpcT5dF+dZKpWsRkTkiMlxVZydV7V4vM0b6v3lRRO4GNiTU+GkGinzW6T6zRGQxYBmCl0+zUfVa2tTyuBg4sxvm5XQCX/5qACKypIgMStsEf/HnGjurTlH3KN8uYiLwzaT9TaCdFiYigxPvKURkKPBZoJkcDop81vY69wDuLJNLr9FUvZZE+KfsQoj6dpoZVfVXHV/A7oQnsAXAHODWpH9F4Kak/Sng6eT1T8IyU8Pn3tHrSf7emeDJ859mvR6CXeEO4N/J+5CkfyzBnRVCVtZnk//Ns8ABjZ53znW0+6yBnwC7JO0lgL8C04BHgU81es6duJbTk9/I08BdwFqNnrO/Kr/cpdhxHMepG7785TiO49QNFyqO4zhO3XCh4jiO49QNFyqO4zhO3XCh4jiO49QNFypOTYjI/iKiSenRRs5jQjIPTQIU0/6Rpr9d5LWI3G/3T/pURH6Wt11EVjbjqYi0C64UkRNE5N+V5mX2XVNE/iAir4jIh8n7H0VkzTLXOCtpH5iMOdJsv92c60/VPzXH6XpcqDitzGvAZsChZbaPr8M5Xk/OcViFfXYDbqg2r0QgPQF8GjgR2J5QB3wd4Ik8gVWFI5LzvFbjcY7TZXiaFqeVWaCqD5fZdhuwg4iMU9VJHT2Bhuy5D4vIEnnbk4jvjYFjKs1LRJYlZOF9GthWVT9INt0rIn8B7gSuEpE1tTQ1SaW5TUnGXlDLNTlOV+KaitMliMh+IvK0iHwgIm8kSzzD2+wzUER+JyJvisi7InK9iGz+/+3dTYhVdRjH8e8jFVTgIqIEE2fRiwtpExUIoSCRUMIUBBYuJCSCFhWkSJlBtqgIslolQoEVVBNB9GJp0uDAUFG0CbwU2CsaJSUjLYr4tXj+l+5cbvfemfmf7gz8Phs55/zvOc+sHv+vTxnO2bbAECbIXsHjDR/7Pg78AkwPaLed3NF/X0dCAaBc31+eb28iSLP/i5OKVVeO8j9IntN0G7ALuAmYjIjO4lf7gbuAp0u7FvBKpTAE7AauBposUDUOvC1pUEnWjcApSZ/1eijpU/IYHJ/Aa0uah7+sqlLNby/wsaQtHfePA8fIJPJcmZi+E9gl6anS7HBEXEDOFSyYpPcjYgp4LCImlLVFqomI5cAGMiEOsoosYdzPt8w+tddsyXFPxWq7iix+NavHIWkK+A5YX25dT9bTeKPr9xOV43kIuALYVvm9ADcDfwJHhmg7zBDcYqzOaDYnTipW20Xl35M9np3qeN6eX+muafJzzWAkHQMOAXvaR9pXNA4ckjTMRPkPwNiANquZXQnRbMlxUrHa2sWgelVfXAG0Vza1k053Sd9LG4jpYbIA1D21Xljqf2xi9lLifj4CVkTEtf/xvuvIv/1onQjNRsNJxWprkb2NLZ03I2Id+T/xyXLrE3Iy/fau33dfL5ikL4A3yaGwCyu9diNwPvDukO0PAL8Bz3YvTy7X+8iEfKBSfGYj4Yl6m69NEdG96e6MpMMRsQd4oezyfhlYSW5E/Bp4EUBSKyJeBfZGxDLgc3Ll0+byrkGrqebqEbK65iX8m9gWYhyYlPT7MI0l/RoRdwBvAdMR8QxwghwSewBYA9w67B4Vs8XKScXm6/ke974C1kraHxF/ADvIkr1ngfeAnZLOdrS/G5gBdgLnkUM/9wLvAGdqBivpeEQcpMKEfdn3spk57tiX9EFEXEP2mJ4ALiaHA48CWyUtprLFZvPiyo+2qETEDuBJYEzS933avUQu570ckKS/G4zpHHLV2hHgRjJJTgOrJP04wriWkUPY3wBTkrY29S2zYbmnYiMTEbcAa4EvyeGuG4AHgdf7JZQOq4G/yOGsDQ3FeBldK7LKESz9lv82HlfxITm3AzDV4HfMhuaeio1MRKwneyVryAn0n4DXgEe7jzLp8dsxcvgIYEZSq6EYzyUPgGxrSZoZdVzlW1cCy8vlaUknmvqW2bCcVMzMrBovKTYzs2qcVMzMrBonFTMzq8ZJxczMqnFSMTOzapxUzMysmn8A9xtFuiTi85oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}