check_null_curl_spec.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np, plottools as pl, local, tools_cmb, quad_func, binning as bn, analysis as ana\n",
    "from matplotlib.pyplot import *\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_spec(qobj,rlz=None,cn=2,bnum=10,lmin=40,lmax=2048,ylim=[-2,2]):\n",
    "    mb = bn.multipole_binning(bnum,lmin=lmin,lmax=lmax)\n",
    "    Mkk, Vkk, Skk, Okk, N0 = local.load_spec(qobj,mb,rlz=rlz,cn=cn,outN0=True)\n",
    "    st0 = ana.PTEs(Okk,Skk,x1pte=True)\n",
    "    st1 = ana.PTEs(Okk,Skk,x1pte=True,diag=True)\n",
    "    # plot binned spectrum lin-scale\n",
    "    pl.plot_1dstyle(fsize=[10,4],xlim=[mb.lmin,lmax],ylim=ylim,ylab=r'$10^2L^{0.5}C_L^{\\varpi\\varpi}$')\n",
    "    #errorbar(mb.bc-5,(Okk-Mkk),yerr=Vkk,fmt='o')\n",
    "    errorbar(mb.bc-5,mb.bc**0.5*(Okk-N0/100.-Mkk+N0/99.)*1e6,yerr=mb.bc**0.5*Vkk*1e6,fmt='o')\n",
    "    #errorbar(mb.bc+5,mb.bc**0.5*N0*1e6/100.,yerr=mb.bc**0.5*Vkk*1e6,fmt='o')\n",
    "    #plot(aobj.l,aobj.ckk,color='k',ls='--')\n",
    "    plot(aobj.l,aobj.l**0.5*aobj.ckk*1e6,color='k',ls='--')\n",
    "    figtext(.7,.80,r'$\\chi$-PTE:'+str(np.round(st0.px1,2)))\n",
    "    figtext(.7,.75,r'$\\chi^2$-PTE:'+str(np.round(st0.px2,2)))\n",
    "    figtext(.7,.70,r'$\\chi^2$-PTE (diag):'+str(np.round(st1.px2,2)))\n",
    "    axhline(0,color='k')\n",
    "    show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "dirs = local.data_directory()\n",
    "aobj = local.init_analysis_params(qid='boss_s15d',fltr='none',wind='com15',ivar='noivar',ptsr='base')\n",
    "bobj = local.init_analysis_params(qid='boss_s15d',fltr='cinv',wind='base',ivar='base',ptsr='base')\n",
    "#bobj = local.init_analysis_params(qid='comb_dn',fltr='none',wind='com15',ivar='noivar',ptsr='PT')\n",
    "qobj = quad_func.reconstruction(dirs['local'],aobj.ids,stag=aobj.stag,rlmin=500,rlmax=3000,qlist=['TT'])\n",
    "robj = quad_func.reconstruction(dirs['local'],bobj.ids,stag=bobj.stag,rlmin=500,rlmax=3000,qlist=['TT'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chi: 12.4 , chi (sim mean): 0.0 , chi (sim std): 6.1 , PTE: 0.04 \n",
      "chi^2: 28.3 , chi^2 (sim): 11.1 , PTE: 0.01 \n",
      "chi: 12.4 , chi (sim mean): 0.0 , chi (sim std): 6.1 , PTE: 0.04 \n",
      "chi^2: 29.3 , chi^2 (sim): 11.1 , PTE: 0.00 \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAEOCAYAAAAEzVrIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAA9hUlEQVR4nO3deXxU1f3/8dcngSRkQdlJQhApIDsBEpeqgFYRrAQRQRYVf9K64UbVuqDWrV93bal+jfDVaisCrQtQBBUqoEUtBA3IIgRFhbCDEgIJIeH8/pjJmJCFBCaZzOT9fDzmwdxzztz53Mtk8sm595xjzjlEREREJPiFBToAEREREfEPJXYiIiIiIUKJnYiIiEiIUGInIiIiEiKU2ImIiIiECCV2IiIiIiEiqBI7M0sys0VmttbM1pjZbeW0MTObbGYbzWyVmfUJRKwiIiIita1BoAOopkLgDufcF2YWB6wwswXOubUl2gwGOnofZwAvef8VERERCWlB1WPnnNvmnPvC+3w/sA5IPKrZUOBvzuNz4GQzi6/lUEVERERqXbD12PmYWTugN/Dfo6oSgc0ltrd4y7Yd9frrgOsAYmJi+nbu3LnGYhURERHxlxUrVux2zrUory4oEzsziwXeBm53zuUczz6cc1OAKQApKSkuIyPDjxGKiIiI1Awz+76iuqC6FAtgZg3xJHXTnHPvlNMkG0gqsd3GWyYiIiIS0oIqsTMzA14B1jnnnqug2Rzgau/o2DOBfc65bRW0FREREQkZwXYp9mzgKuArM8v0lt0HtAVwzqUD84CLgY3AQeD/1X6YIiIiIrUvqBI759x/ADtGGwdMqJ2IREREROqOoLoUKyIiIiIVU2InIiIiEiKU2ImIiIiECCV2IiIiIiFCiZ2IiIhIiFBiJyIiIhIilNiJiIiIhAgldiIiIiIhQomdiIiISIhQYiciIiISIpTYiYiI1GHh4eEkJyfTvXt3RowYwcGDB9mzZw/JyckkJyfTunVrEhMTfdsFBQW+1xQ/nnjiiTL7ff/99znttNPo0KFDufVVaVdUVETv3r255JJL/H7ccnzMs7Rq/ZaSkuIyMjKq1PbIkSM45wgPD6/hqERERCA2Npbc3FwAxo4dS9++ffnd737nq3/ooYeIjY3lzjvvLPc15SkqKqJTp04sWLCANm3akJqayvTp0+natWu12j333HNkZGSQk5PD3Llz/XnYUgkzW+GcSymvTj12VbB8+XJuuukmOnXqRGRkJA0aNOCZZ54B4ODBgyxevJiioqIARxlYY8fmMnZsxV8iIiJS2oEDB0hMTOS+++4DPL9rkpOTyc/Pr/A15557Lhs3bjzh9162bBkdOnSgffv2REREMGrUKGbPnl2tdlu2bOG9997jN7/5zQnHI/6jxO4Ytm/fzi9/+Utef/11unTpwp133snDDz/MoEGDAFi4cCHnnXceHTt25KmnniInJyfAEYuISDCIiYlh1apVTJ8+nYMHD3Lttdfy2muvERUVVW77wsJC5s+fT48ePY6577y8vFKXYmfOnAnAxRdfzNatW8nOziYpKcnXvk2bNmRnZ5fZT2Xtbr/9dp566inCwpRK1CUNAh1AXZWbm0tsbCytW7dm1qxZ9OvXj7i4uDLtzj//fKZPn056ejp33303Tz75JHfffTe33347ERERAYhcRESCRbNmzYiOjmb8+PFcddVVJCcnl2lTnKSBp8du/Pjxx9xvo0aNyMzMLFM+b968E4zYY+7cubRs2ZK+ffuyePFiv+xT/ENpdjm+//57TjvtNGbMmAHAr3/963KTOvDcxzBq1CgWL15MRkYGZ5xxBq+//rruwRMRkSrp2bMn27dv990j9+KLL/p62rZu3epL0jIzM/nLX/7il06DxMRENm/e7NvesmULiYmJVW63dOlS5syZQ7t27Rg1ahQfffQRV1555QnHJSdOid1R8vLyGDp0KAcOHKBnz57Vem3fvn2ZN28en376KeHh4ezbt4+xY8eyYcOGGopWRESC2a5du1i0aBHDhw/3XdKcMGGCL5FLSEiokfdNTU0lKyuLTZs2UVBQwIwZM0hLS6tyu8cff5wtW7bw3XffMWPGDM4//3zeeOONGolVqkeJ3VEeeOABVq5cWe7ooKo66aSTAFi5ciVz586lZ8+ePP3004wevU8DDERExGf8+PGcf/75rFy50q/7Pfoeu3vuuQf4+R67Bg0a8MILL3DRRRfRpUsXRo4cSbdu3Xyvr2o7qXuCaroTM3sVuATY6ZzrXk79AGA2sMlb9I5z7pFj7bd4upMVK1aQmprK9ddfz0svveSXmLdt28aNN97I7NmzadYslTPPfIm5c/v6Zd91SXHCOm1abIAjEREJDi+//DIfffQRzzzzDIMHD2b16tWBDkmCRChNd/IaMOgYbT5xziV7H8dM6krKyMggISGh0okaqys+Pp53332XadOmsX//RjIzH/DbvkVEJDhlZWXx/PPPk56eTlJSEvHx8QwcODDQYUkICKpRsc65j82sXU3t//rrr+fqq6+mUaNGft2vmTFmzBjeeusMnPPMd5ednc2hQ4do3769X99LRETqvo4dO/L111/7thcsWBDAaCSUBFuPXVWcZWYrzWy+mVX5RoB169YB+D2pK6lRo1ZER3tuhL3tttvo2bMnU6ZMIZguh4uIiEjdFWqJ3RfAKc65XsBfgFkVNTSz68wsw8wytm7dSteuXfnrX/9aW3Hy3HPPceaZZ3L99dfz61//mm3bttXae4uIiEhoCqnEzjmX45zL9T6fBzQ0s+YVtJ3inEtxzqUUFhYSFRXF0KFDay3Wtm3b8uGHHzJ58mQWL15Mt27d+PTTT2vt/UUkdGhJv9Aya9Ysfvvb33LFFVfw4YcfBjocCTIhldiZWWszM+/z0/Ec355jvW7v3r0MHTqUpk2b1nSIpYSFhXHLLbfw5ZdfMmDAgOOeXkVERELHpZdeytSpU0lPT/ctBVae8PBwkpOT6d69OyNGjCA7O9s3vUnr1q1JTEz0bRcUFPjaFz8qGij4/vvvc9ppp9GhQ4dqt7n22mtp2bIl3buXmbhCaklQDZ4ws+nAAKC5mW0B/gA0BHDOpQOXAzeaWSGQB4xyVbiBraioiKuvvrrG4j6W0047jXfeeQeAgoIChg0bxs0338zgwYMDFpOIiATWY489xoQJEyqsL7ls2NixY5k5c6Zv+6GHHiI2Nta3msXR7StSVFTEhAkTWLBgAW3atCE1NZW0tLRSHQ+Vtbnmmmu4+eabA/o7tb4Lqh4759xo51y8c66hc66Nc+4V51y6N6nDOfeCc66bc66Xc+5M51yVrm02aNCgzgwz37FjBz/88AMXX3wx119/Pbm5urwiIhKKDhw4QGJiIvfddx8Ay5cvJzk5mby8PO6++24GDx5Mnz59qrSvc889l40bN55wTMuWLaNDhw60b9+eiIgIRo0axezZs6vcpl+/frV+9UtKC6rErqZ07NiRBg3qRudlUlISGRkZ/P73v2fq1Kn07NmTTz75JNBhiYgERCjfPxgTE8OqVauYPn06Bw8e5Nprr+W1115j6tSpLFy4kLfeeov09PRj7qewsJD58+fTo0ePStsdvRpF8WXe4lUmwDMVV1JSku81bdq0ITs7u9R+qtJGAqduZDMBFh0dHegQSomMjOTJJ59kyJAhjBs3jltvvZUVK1b41hEUEZHQ0KxZM6Kjoxk/fjxXXXWVL+m69dZbj/na4kQNPD1248ePr7R9RZdi582bdzyhSx2lxK4OO+ecc1i5ciW7du0iLCyMnJwcvvnmG3r37h3o0ERExE969uzJ9u3bS90Pd7QXX3yRqVOnAp5ELCEhoUr3zFVXYmIimzdv9m1v2bKFxMTEareRwFEXUB0XGxvLqaeeCnhuhj399NN54IEHyM/PD3BkIiJyonbt2sWiRYsYPnx4pVdlJkyYQGZmJpmZmSQkJNRYPKmpqWRlZbFp0yYKCgqYMWMGaWlp1W4jgaPELojcf//9jBo1iscee4wePXqwcOHCQIckIiInYPz48Zx//vmsXLmyxt/r6Hvs7rnnHqD0PXYNGjTghRde4KKLLqJLly6MHDmSbt26lWpXWZvRo0dz1llnsX79etq0acMrr7xS48clpZmWs4KUlBSXkZFR4+9TfAPwtGmxJ7SfhQsXcuONN7Jx40b+9Kc/cdttt/kjvBPir2MTkeoL5Z+/UD62l19+mY8++ohnnnmGwYMHs3r16kCHJEHCzFY451LKq1OPXRC64IIL+Oqrr/jDH/7ApZdeCsCePXsoKioKSDxz5hSQmVnEsmVF9O+fw5w5BQGJQ0QkWGRlZfH888+Tnp5OUlIS8fHxdWbaLQluGjwRpKKionjooYcAcM4xduxYtm/fzp///Gf69+9fa3HMmVPApEl5FHhzua1bHZMm5QGQlhZRa3GIiASTjh078vXXX/u2FyxYEMBoJJSox66W1HSv1rXXXsvevXsZMGAAI0eO5Pvvv/fr/ivy7LP5HD2OIz/fUy4iIiK1S4ldLaioV8tfyZ2ZMXLkSL7++msefvhh5s6dS+fOnWtl8eht28q/R7OichEREak5SuxqQW31akVHR/Pggw+yfv16xo0bx5lnngnAN998w6FDh/z6XsXi461a5SIiUrlZs2bx29/+liuuuKJW/kCX0KLErhbUdq9WUlIS6enpNG7cmKKiIoYOHUrnzp35+9//7vcBFnfcEUVUVOmyqChPuYiIVN+ll17K1KlTSU9P9y37VZ7w8HCSk5Pp3r07I0aMIDs72zeVSevWrUlMTPRtFxQU+NoXP5544oky+8zLy6N///7l/q546KGHeOaZZ3zbv/zlL4/r+AoKCujXrx+FhYXl1r///vucdtppdOjQodwYq9uuvlFiVwsC2asVFhbGc889R5MmTbj66qvp3bs3//rXv/DXNDdpaRH88Y+NiPCOk0hIMP74x0YaOCEicoIee+wxJkyYUGF98coTq1evJiIigpkzZ/omMb7hhhuYOHGibzsiIsLXvvhRPI9dSa+++iqXXXYZ4eHhx4zv008/Pa7jioiI4Fe/+lW5SWtRURETJkxg/vz5rF27lunTp7N27drjblcfKbGrBYHs1TIzBg4cSEZGBjNmzCAvL4+0tDRmzZrlt/dIS4sgOTmc008PZ8mSxkrqRESq4MCBAyQmJnLfffcBsHz5cpKTk8nLy+Puu+9m8ODB9OnTp0r7Ovfcc9m4ceMJxzRt2jSGDh3q2/7jH/9Ip06dOOecc1i/fn2ptrGxnrkFL730Uvr27Uu3bt2YMmVKqTaPPvoop512Gueccw6jR4/29fhdeumlTJs2rcz7L1u2jA4dOtC+fXsiIiIYNWoUs2fPPu529ZESu1pQF3q1wsLCuOKKK1i7di1vvPEGQ4YMAeDNN9/k7bff5siRI7UWi4iIQExMDKtWrWL69OkcPHiQa6+9ltdee42pU6eycOFC3nrrLdLT04+5n8LCQubPn0+PHj0qbXf0yhNH95gVFBTw7bff0q5dOwBWrFjBjBkzyMzMZN68eSxfvrzc/b766qusWLGCjIwMJk+ezJ49ewBPovr222+zcuVK5s+fT8mFALp3715qf8WrWmRnZ5OUlOQrb9OmDdnZ2WXes6rt6iPNY1dL0tIimDnTMwo2kDOoN2zYkLFjx/q2X375ZT7++GO6d+/O/fffz/Dhw2nQQB8LEZHa0KxZM6Kjoxk/fjxXXXWVL+m69dZbj/na4kQNPD1248ePr7R98aXYiuzevZuTTz7Zt/3JJ58wbNgwoqOjASpcD3by5Mm8++67AGzevJmsrCyaNWvG0qVLGTp0KFFRUURFRfk6FMBzf2BERAT79+8nLi6OefPmHfN4pWrUY1fPffTRR0ybNo2ioiJGjRpFp06d+Ne//hXosERE6o2ePXuyfft27rzzzgrbvPjii76kr3hd15L3zP3lL38hIuLErgI1atSI/KOncDiGxYsXs3DhQj777DNWrlxJ7969q7yPQ4cOEXXUfUqJiYls3rzZt71lyxYSExPLvLaq7eojJXb1XHh4OGPGjOGrr77inXfeoXXr1jRq1AiAXbt2sX379gBHKCISunbt2sWiRYsYPnw4YWEV/0qeMGGCL4lLSEiokViaNGlCUVGRLzHr168fs2bNIi8vj/3795f7R/++ffto0qQJ0dHRfP3113z++ee+urPPPpt//etf5Ofnk5uby9y5c311e/bsoXnz5jRs2LDU/lJTU8nKymLTpk0UFBQwY8aMcnsKq9qupo0dm+tbz7iuUGIngCfBGzZsGJ9++im/+tWvAHj88cc55ZRT+O1vf8uaNWsCHKGISOgZP348559/PitXrqzx9zr6HrvyRsUOHDiQ//znPwD06dOHK664gl69ejF48GBSU1PLtB80aBCFhYV06dKFe+65xzd/KniSr7S0NHr27MngwYPp0aMHJ510EgCLFi3i17/+ta9t8T12DRo04IUXXuCiiy6iS5cujBw5km7dupVqA1Tart5zzgXNA3gV2AmsrqDegMnARmAV0Kcq++3bt6+rDWPG7Hdjxuyvlffyhw0bNrgbbrjBRUVFOcANGDDAzZo1q9y2wXZsIqFi9uxDrmvXn1yHDj+5fv32udmzDwU6JL8K5e+W9PR0N3LkSPfDDz+4bt26BToc55xzK1ascFdeeaXf9rd/v+f/7sCBA65v375uxYoVzjnnhg0b5tavX++39wmUQH0+gQxXQU4TbD12rwGDKqkfDHT0Pq4DXqqFmEJWx44deemll9i8eTNPPvkkmzZt4u233/bV7927N4DRiUhNL1coNScrK4vnn3+e9PR0kpKSiI+PZ+DAgYEOiz59+nDeeef5bTL76667juTkZPr06cPw4cPp06cPBQUFXHrppXTq1Mkv7yGlmfPTRLW1xczaAXOdc93LqXsZWOycm+7dXg8McM5tq2yfKSkpruQw7JpSfB0+kKNiT0RRURE5OTk0adKEzMxMTj/9dEaMGMGNN97I//5vL8wsaI9NJBj175/D1q1lv8MTEowlSxoHICL/C/bvTQltgfp8mtkK51xKeXXB1mN3LInA5hLbW7xlZZjZdWaWYWYZu3btqpXggl14eDhNmjQBPEP0b7rpJubOncu5557Le+/1Ze3aP7F///4ARylSf9T2coXif3Xx5nt/CvXjq4tCLbGrMufcFOdcinMupUWLFoEOJ+gkJSXxpz/9ia1bt/LXv/6VyMhmrF79JGaeZdK+//77CtcBFBH/CORyhbVhzpwCMjOLWLasiP79c3SJWaQKQi2xywaSSmy38ZZJDYmJieGaa67hwgsXkJa2itjYWJxzXHLJJZxyyince++9rFu3LtBhioSkQC5XWNN0/6DI8Qm1xG4OcLV5nAnsO9b9deI/UVGenk/nHI888gi9e/fmqaeeomvXrqSkpGgdPxE/qwvLFdaUZ5/N5+h5bvPzPeUiUrGgSuzMbDrwGXCamW0xs/FmdoOZ3eBtMg/4Fs90J1OBmwIUar0WFhbGsGHDmDt3LtnZ2Tz33HMcOXLEd//d9u3bfWsjisiJSUuLIDk5nNNPD2fJksYhkdSB7h8UOV5BtSioc270MeodMKGWwpEqaN26NRMnTmTixInFcw3y1ltvccsttxATE0NaWhojR45k0KBBZZaWEZH6Kz7eyh3xGyr3D4rUlKDqsZPgVjyw4qabbmLRokWMHTuWDz/8kGHDhhEfH09urkZOiYhHKN8/KFKTgqrHTkJDWFgYAwYMYMCAAbzwwgssXryYzMxMYmM98wCNHDmS6OhoLr/8ci644AL15InUQ8WXlO+91zOAIiHBuOOOqJC51CxSU9RjJwHVsGFDLrzwQu666y7AM/DipJNOYtasWQwZMoTmzZszYsQIFixYEOBIRaS2her9gyI1SYldLZo2LVazpx+DmTF16lR27NjB/Pnzueqqq1i6dCmZmZkA7Nu3j8mTJ/P9998HNlAREZE6SImd1EmRkZEMGjSIl156iS1btnDLLbcA8PHHH3PbbbfRrl07evfuzcMPP0xmZibBtjSeiIhITVBiJ35Rk72RYWFhvvvshgwZwoYNG3j66aeJjY3l4Ycfpnfv3mzYsAHwrHixZ8+eGolDRESkrlNiJ0GnY8eO3HnnnXzyySds376dN998k06dOgFw//3306JFC8444wz+8Ic/8Omnn2ppMxERqTeU2ElQa9myJaNHj/ZNpTJx4kQeeughwsPDeeyxxzj77LM5++yzfe337dsXqFBFRERqnKY7kZDSp08f+vTpw4MPPsjevXtZuHAhRUVFABQWFnLqqafSunVrBgwYwHnnnceAAQNo0aJFgKMWERHxD/XYSchq2rQpI0eOZPRoz4IlBQUF3HvvvZxyyin87W9/Y+TIkbRs2ZJnn30WgMOHD+v+PBERCWrqsZN6Izo6mrvuuou77rqLw4cPs2LFChYtWkS/fv0AWLp0Keeddx49e/ZkwIABvsu4iYmJAY5cRESkaqrcY2dmp9dkICK1qWHDhpx55pnce++9pKamAtCuXTseffRRWrRowdSpU7niiito06YNGRkZAHzzzTdkZmb6Lu2KiIjUNdXpsZtlZm8CBcAK4APnnBb3lJDRrl077r//fu6//34OHz5MZmYmS5cupWfPngC8+OKLPP/888TGxnLmmWdy9tlnc9ZZZ3HhhRcSFqa7GkREJPCqk9g54AcgArgUeNrMbnTOfVATgYkEUsOGDUlNTfX15gH87ne/IyUlhaVLl7J06VIeffRRmjVrxo4dOwBP4peXl0dqaip9+/b1rX0rIiJSW6qT2O10zk0u3jCzeGAeoMRO6oU2bdowZswYxowZA0BOTg6bNm3yTbXyzjvv8NFHHwGepdG6dOnC5ZdfzsMPPwx4Bm9ERGitSxERqTnV7bHDzK4B4oDGaFSt1GONGzemV69evu1///vf7Ny5k4yMDJYvX87y5cs5cOAAAM454uPjad68Ob169aJXr14kJyeTkpJCq1atAnUIIiJynObMKSAzs4iCAujfP4c77ogiLS3wf7wfz6jYNcBg7+Mu/4YjEtxatmzJxRdfzMUXX1yqvKCggFtvvZWVK1eSkZHBP//5TwDuvfde/ud//ofc3FwefPBBevXqRY8ePejcuTPR0dGBOAQRETmGOXMKmDQpj4ICz/bWrY5Jk/IAAp7cVSexa2xmzZxzy4HlZvYk8DnwYc2EJhI6IiMj+cMf/uDbzsnJYdWqVb7eum+//ZaXXnqJ/Px8wHMpt127dkyePJlLLrmEH3/8kaysLLp06UJcXFxAjiGYjR3rGedVU+sZixyPutrj4y+hfHzPPpuP9+vaJz/fUx7oY6xOYtccyDCzb4F1wMmA1URQlTGzQcCfgXDg/5xzTxxVfw3wNJDtLXrBOfd/tRqkyDE0btyYc845x7fds2dPcnNzycrKYs2aNaxdu5a1a9fSunVrABYtWsTw4cMBSEpKomvXrnTt2pXbb7+dtm3bUlhYSHh4uO9+PxGp2+pyj48/hPrxbdvmqlVem6qT2A0GlgPnAh2BI8AdNRFURcwsHHgRuBDYgqfncI5zbu1RTWc6526uzdhETlR4eDidO3emc+fOviSu2LnnnsusWbN8Cd/atWtJT0/npptuAjwjch944AE6duxIx44d6dSpEx07duTyyy+nUaNGgTgcEalEXe7x8YdQP774eGPr1rJJXHx84P+4rk5i9zYQ6HnsTgc2Oue+BTCzGcBQ4OjETiSktGjRgqFDhzJ06FBf2ZEjR3w9dMnJyYwbN46srCyWL1/OP//5T44cOeJLEB955BHmzZtHx44dad++Pe3atePUU0+lf//+6uUTCYC63OPjD6F+fHfcEcWkSXmlkteoKE95oAXbPHaJwOYS21uAM8ppN9zM+gEbgInOuc1HNzCz64DrANq2bVsDoYrUrJKTIvfv35/+/fv7tg8dOsQPP/zgG4DRsmVLYmJiWLJkCdOmTcM5R7Nmzdi9ezcAt99+O2vWrOHUU0/1PTp16kSfPn1q96DkuOjeweBTl3t8/CHUj6+41/Heez2XmxMSrM7cQxiK89j9C5junDtkZtcDrwPnH93IOTcFmAKQkpISGn9CiHhFRkbSsWNH3/YNN9zADTfcAHhG6P7www++pA4gNjaWnJwcZs2axa5duwDo3bs3X3zxBQBjxoxhz549JCUl0bZtW5KSkujWrRunn66VBkWOR13u8fGHUD8+8CR3M2d6biKsS39cBds8dtlAUontNvw8SAIA59yeEpv/BzxVC3GJBI2IiAg6dOhAhw4dfGWPPfYYjz32GAC5ubl89913vhG6AHFxcWzcuJFVq1axfft2ANLS0pg9ezbgGfwRFhZWKvE766yzfL2Ihw4dIjIysrYOUaTOq8s9Pv4Q6sdXl9X4PHZmNtM5d8VxvE95lgMdzexUPAndKGDMUe8X75zb5t1MwzOCV0SqKDY2lu7du5cqe/nll33PDx06RHZ2NkeOHAE8ky/369eP7777jh9++IGlS5fy448/MmHCBPr378/hw4eJjo7m5JNPJj4+noSEBBISEhg+fDhDhgyhsLCQZcuWkZCQQHx8vBJAKaUu9YT4W13t8fGXUD++uqo25rFLPUZ9lTnnCs3sZjyXf8OBV51za8zsESDDOTcHuNXM0oBCYC9wjb/eX0Q8l3nbt2/v2zYzXnjhhVJtcnNzKfDOc1BYWMgjjzzC1q1b2bp1K9u2bePrr7+me/fuDBkyhOzsbM4++2zfa5s2bUpCQgIPPPAAI0eOZM+ePbz55pu0atWKli1b0qpVK1q1akWTJk008ENE5CjVSeyaUcE8dmb2GLDa+/jaOVfo5zh9nHPz8NzbV7LswRLP7wXuran3F5Fji439+a/zRo0aMWnSpArbNm/enPnz5/sSv+Lkr3HjxgB8/fXX3HrrrWVe9+abbzJ69GhWrFjBpEmTfAlfcfI3cOBAWrVqxeHDhzlypJCwsOO5QCEiElyq8013PvAV5c9jtwfP3HIT8Vwq3Yrnku1qQP2vInVcIFdmiImJYdCgQRXWn3XWWezcuZOdO3eyY8cO3yM11XMx4ODBg+zdu5d169axY8cODh06BMDHH39Mq1at+Mc//sGMGVcRGdmUzMxWtGjRgubNm/PMM8/Qrl07Vq9ezZdffknz5s1p1qwZzZs3p3nz5sTFxalHUESCTpUTO+fcl96ni7yPknXPl9z23gPX3ftYcIIxikg9FhYWRosWLWjRogXdunUrU3/uueeybNkywHO/X05ODjt37iQxMRGA7t270737PeTn7+S00/aye/du1q79eerLefPmcffdd5fZ75YtW0hMTGTKlCm8+eabvoSvOPmbMGECERERZGdnc+jQISWDIlInVOvahJmlANcDdznnfqqonXNuE7AJz9QjIiK1wsw46aSTOOmkk3xlvXr1omfPXwDl90hOmDCByy67jN27d7N792727NnD7t27ad68OeBZEeTIkSOsXbuWPXv2sGfPHo4cOcItt9wCeEYUp6enA9CgQQOaNGlC69atWbVqFeAZeLJy5UpOPvlkmjRpwsknn0yrVq1IS0sDYNeuXTRs2JDGjRuXmptQROR4VPemk98D+eUldWbWCOiNZxBDgR9iExGpcTExMWWmfylp/PjxjB8/3rd95MgRcnJyaNDA8/X5m9/8hjPPPNOXGP70008UFRX52mdmZvLWW2/x448/+spPOeUUX2J35ZVX8uGHH/qS0iZNmtC7d2/efvttAJ5++ml27drlSwqbNGnCKaecwllnnQXAjh07iIqKIjY2lvDwcP+fIBEJKtVN7FKAq8urcM7lmdlDwFLg4ROMS0SkTgoLC+Pkk0/2bfft25e+fftW2P6ll17ipZdewjlHbm4uP/74Y6k5Am+55RYGDRrEjz/+6Hu0atXKVz9//nw+++yzUq8577zz+OijjwDPpeisrCzAM2ilcePGDBkyxNeLeN1111FYWEjjxo056aSTaNy4Mb169eKCCy4A4L///S8xMTE0btyYxo0bExcXpwRRJIhVN7GLB76rpP7/8Mxtp8RORKQEMyMuLo64uLhS5ZdcckmlrytO4PLz8/nxxx/56aefSl2yffjhh9m2bRs5OTm+R5cuXXz1K1euZOvWreTk5LB//36cc1x77bVccMEFOOc4++yzS/UwAkycOJHnnnuOgoICfvnLXxIbG1vqMXToUIYMGUJeXh6vvPKKrzwuLo7Y2Fh+8Ytf0Lp1a4qKisjPzyc6Olr3HorUkuomdrvxrPywpYL6FUD7CupEROQ4RUVFER8fT3x8fKny0aNHV/q6//73v77nR44c4cCBA6Uml543bx45OTns27fPlxgWjzg+fPgwrVq1Ijc3l23btpGbm0tubi6dO3dmyJAh7Nq1y3evYUnPPfccEydOJCsriy5dumBmxMTE+BLAxx9/nMsvv5ysrCzuu+8+YmNjiYmJITo6mpiYGEaMGEHXrl3Zvn07S5Ys8ZVHR0cTHR1N+/btiY2NpbCwEOccDRs2PNHTKxIyqpvYLQGuBT6roN4ArRciIlIHhYWFleoxDAsLY+DAgRW2j4mJ4b333quwPjExkZ07d/oSvuJH8TrFzZo146mnnipTXzwwZf/+/axZs4b9+/dz8OBBDh48SH5+Pj179qRr1658+eWXjBo1qsz7fvjhh1x44YXMmjWLESNG0KBBA1/SFxMTw1tvvUVycjIffvghkydP9tUV19955520atWKr776iuXLl5eqj4qKIjU1lcjISPbt28fBgweJioqiUaNGREZGqudR6rzqJnbP4Fl14ivn3ORy6s8GvjnxsEREpK4LDw/3TUVTnhYtWnDXXRWvPNmnT59SU8+Ap1fROQdAv379WLNmjS/pK3706tULgK5du/LYY4/5yg8cOMDBgwd9o6IPHDjAtm3byrz++uuvp1WrVrz//vv8/ve/LxNXdnY2CQkJPP/88zz8cOk7i6Kioti5cydxcXE88cQT/P3vf6dRo0a+5K9Ro0bMnj0bM+PNN99k2bJlperi4uK48cYbAVi+fDnZ2d8THh7Fp5829Q2C6dSpE+CZozEsLIyIiAiNmJYqq1Zi55zLNLPfAlPN7HLgL3guv+YC/YEn8SR/IiIi1VYygYmJiaFr164Vtu3atWul9cOGDWPYsGEV1t90001cccUVpZK+vLw8X4/ikCFDiI+PJy8vj/z8fN+/jRo1AiA+Pp4uXbr4yg8cOMBPP/3k69X7/PPPef3118nPz/ctsdekSRNfYvfkk0+yZIln9LP3VkqSkpL44YcfALjsssv44IMPAGjYsCGRkZH07NmTpUuXAnDVVVexevVqoqKiiIyMJCoqiu7du/PMM55fw4888gg7d+4kMjLS9+jUqZPv8v1bb71FQUGBr674cn/xfJEbN26kQYMGpeojIyM1uKaOq/YaO86518xsPfAsMBNw3ioD5gN/9l94IiIiNSMmJoaYmJgK64814nncuHGMGzeuwvrJkyczebLn4lbxQJLilVEAnnnmGXJzb6GoKI+77gojLy/PN40OeKba6d+/P4cOHfI9ipNOgDZt2rBv3z4OHTpEfn4+P/30E3v37vXV//vf/2b16tW+1xYWFjJo0CBfYjdx4kS2bCl9y/zll1/OP//5TwBSU1P56aefStVfc801/PWvfwWgZcuWmBmRkZFEREQQERHB1VdfzT333ENBQQGDBw9m3bpwwsIaMmJENBEREQwbNozLL7+c3NxcHnroId/rih/9+/cnNTWV/fv3M2fOHCIiIkrtv3PnziQkJJCfn8+mTZuIiIigYcOGvkfjxo2JiIjw9frWx0vnx7V4onPuM+CXZtYB6ANEA6udcxn+DE5ERCQUhIeHl0kk27VrR7NmnkRt4MCyk2ePGDGi0n0+/vjjldYvWbKk1HZRUVGpEdCfffaZ777G4uSv5FQ+U6ZM4cCBA6Xqi3vznHOMHj2agoICCgoKOHToEAUFBb6peoqKijh8+DAFBTkcOVLA2rWFFBQUkJKSAkBubi7p6ekUFBRw+PBh33s+/fTTpKamsm3bNq688soyx5Sens7111/P2rVry02633jjDcaOHcvHH3/MgAEDSiV9DRs25I033uCiiy5i8eLF3HjjjWXq//znP9OnTx8++eQTnn/++TL19913H+3atWPZsmXMmjWLVasgLKwhTz0VQ8OGDRk3bhxNmzZlzZo1fPnll2VeP2DAAKKiotiyZQvbt28vU9+2bVvCwjxJ/pEjR3zl1UlQT2hVbOfcRmDjiexDREREal54eHipy6ht2rSptH1liaWZ8ec/V3yBrlGjRnz88ccVrkPdunVrcnM9dc45Dh8+zKFDh3wjnNu1a8eGDRvKJI7FE4m3a9eOmTNn+soPHz7M4cOHfYlj27ZteeCBB3zlxY+EhATAM+djjx49ytQXn5+cnByysrJK1RUUFPhGgWdmZvLUU0/5EmXvQjNccsklNG3alPfee6/cpQq3bdtG69atmTJlCo8++miZ+v379xMbG8ukSZN4/vmfV2sNDw8nIiKCAwcOHDPJs+Luyqoys+54eunWOOdWVOvFdVRKSorLyFBno1Ssoi+nUKHjEwmMUP9shvrxjRmznyNHDvPKK5EcPnzYN8H3vn372LVrV6mksDjxbNiwIRs2bGD9+vVlEsurrrqKBg0asGjRIjIyMkolrUVFRTzxxBMAmNkK51xKeTEds8fOzP7tnPuV9/kY4D5gHnCrmb3mnHvBf6dIRMS/5swpIDOziIIC6N8/hzvuiCItTbMyiciJMzPCwyPK3Kt59JrVR+vUqZNv9HN5zjvvPM4777zjiqkql2Kblnh+G3CBc267mcUCnwJK7ESkTpozp4BJk/LwDkhk61bHpEl5AEruRCQkVWViHDOzRmYWA4Q557YDOOdygaLKXyoiEjjPPptPiSVWAcjP95SLiISiqvTYnQyswTOdiTOzeOfcNm+PXf0bRywiQWPbtvLvIa6oXEQk2B2zx84518451945d6r3323eqiNAxTM/1hAzG2Rm681so5ndU059pJnN9Nb/18za1XaMIlI3xMeX/7dnReUiIsHuuNcocc4ddM5t8mcwx2Jm4cCLwGCgKzDazI6ednw88KNzrgPwPJ7VMESkHrrjjiiiokqXRUV5ykWk5k2bFhuyI2LrqmpPd1LqxWZfOOf6+DGeY73fWcBDzrmLvNv3AjjnHi/R5gNvm8/MrAGwHWjhKjnQuLg4V9ns4iLr1nluJ+3SJTSX0gnl4ztwoB97916FWXMaNNhNXNzfiYn5ONBhiYgctyVLlhz/dCeVqc2kzisR2FxiewtwRkVtnHOFZrYPaAbsLtnIzK4DrgOIjIysqXhFJMBiYj7mhx8WAaGZuIqIlFSlxM7MmgAD8SRNAFuBD5xzP9ZUYDXNOTcFmAKeCYoXL14c2ICkTgvlSTbnzCng3ns9U4I4ZyE5z1so//+JSP1T2eoTx7zHzszGA5/h6RkL8z7OAD711tWmbCCpxHYbb1m5bbyXYk8C9tRKdCJBpqJ53ubMKQhsYCIiclyq0mP3e6CPc+5AyUIzewD4AnilJgKrwHKgo5mdiieBGwWMOarNHGAcnmT0cuCjyu6vE6nPKpvnLdR67URE6oOqJHYOiAMOHFUe562rNd575m4GPgDCgVedc2vM7BEgwzk3B0+i+Xcz2wjsxZP8iUg5NM+biEhoqUpidyewxMxW8/NlzzZAN+COmgqsIs65eXjWqi1Z9mCJ5/nAiNqOSyQYxccbW7eWTeI0z5uISHA6ZmLnnJtrZvOB04EEb/FWYJlzTkuKiQSxO+6IYtKkvFKXYzXPm4hI8KrSqFhvAvdZDcciIrWs+D664lGxCQmhOSpWRKS+OKF57EQk+KWlRTBzpmcUrKYDEREJbse9pFh5zKyzma305z5FREREpGr8mtjhGana3c/7FBEREZEq8HdiJyIiIiIBUq177MzsNWCF95HpnDtYE0GJiIiISPVVd/BEEpAGnAwUmdkGfk70vgAK/RqdSB0wZ04BmZlFFBRA//45GjUqIiJ1VrUSO+fcrwDMrD3Q1/voAzwANC1u5s8ARQKporVUASV3IiJS5xzXdCfOuW+Bb4F/FpeZWTsgBU+iJxIStJaqiIgEE7/NY+ec+w74DnjLX/sUCTStpSoiIsGkSqNizayRmSWWU97N/yGJ1B0VrZmqtVRFRKQuOmZiZ2aXA1nAe2a2yszOKFH99xqLTKQOuOOOKKKOWjZVa6mKiEhdVZUeu/uBvs65ZOD/Aa+Y2RhvnbotJKSlpUXwxz82IsJ7O11CgvHHPzbS/XUiIlInVeUeu4bOuR0AzrkVZtYPeNfMOqARsFIPaC1VEREJFlXpsdtpZj2LN5xze4ELgS5AzwpfJSIiIiK1qiqJ3VXAzpIFzrkC59xooH+NRCUiIiIi1XbMS7HOuS2V1C31bzgiIv6nS+giUl9Uex47M2sD3Aj8EmjtLd4OLAVeds5t9l94pd63KTATaIdnvryRzrkfy2lXBHzl3fzBOZdWE/GIiIiI1DVVmseumJmdA6wDRgBrgDe9jzXFZWZ2tr+D9LoH+LdzriPwb+92efKcc8neh5I6ERERqTeq22P3J+Cvzrlby6s0sz9726SeWFjlGgoM8D5/HVgM3F0D7yMiIiISlKrVYwd0A16spP4loPvxh1OpVs65bd7n24FWFbSLMrMMM/vczC6taGdmdp23XcauXbv8HauIiIhIratuj9024GxgfQX1Z3vbHBczW8jP9+2VNKnkhnPOmVlFc+id4pzLNrP2wEdm9pVz7pujGznnpgBTAFJSUjQfn4iIiAS96iZ2zwDpZnY6sADY4S1vhWduu2uA2483GOfcBRXVmdkOM4t3zm0zs3iOmoKlxD6yvf9+a2aLgd5AmcROREREJNRUK7Fzzv2vme0BJgLjgXBvVRGwArjaOfcP/4boMwcYBzzh/Xf20Q3MrAlw0Dl3yMya4+lBfKqG4hEJGZoOREQkNFR7uhPn3Exgppk1BJp7i3c75w77NbKyngD+YWbjge+BkQBmlgLc4Jz7DZ7VMF42syN47h98wjm3tobjEhEREakTqp3YFfMmcqXupzOzJOBh59y1JxpYOe+3B/hVOeUZwG+8zz8Fevj7vUVERESCQXVHxR5LUzyXSUVERESkllWrx87Mrj5Gk7YnEIuIiIiInIDqXop9DTgIVDQ9iL97AEVERESkiqqbiG3FM/I1rrwHnlGoIiIiIhIA1U3sVgB9Kql3gB1/OCIiIiJyvI5nguLKJrzaCJx3/OGIiIiIyPGq7gTFnxyj/gCw5IQiEhEREZHjosEOIiIiIiHiuCcoFqlPtOSWiIgEA/XYiYiIiIQIJXYiIiIiIUKJnYiIiEiIUGInIiIiEiKU2ImIiIiECCV2IiIiIiFCiZ2IiIhIiFBiJyIiIhIilNiJiIiIhIigSezMbISZrTGzI2aWUkm7QWa23sw2mtk9tRmjiIiISCAFTWIHrAYuAz6uqIGZhQMvAoOBrsBoM+taO+GJiIiIBFbQrBXrnFsHYGaVNTsd2Oic+9bbdgYwFFhb4wGKiIiIBFgw9dhVRSKwucT2Fm9ZGWZ2nZllmFnGrl27aiU4ERERkZpUp3rszGwh0LqcqknOudn+fC/n3BRgCkBKSorz575FREREAqFOJXbOuQtOcBfZQFKJ7TbeMhEREZGQF2qXYpcDHc3sVDOLAEYBcwIck4iIiEitCJrEzsyGmdkW4CzgPTP7wFueYGbzAJxzhcDNwAfAOuAfzrk1gYpZREREpDbVqUuxlXHOvQu8W075VuDiEtvzgHm1GJqIiIhInRA0PXYiIiIiUjkldiIiIiIhQomdiIiISIhQYiciIiISIpTYiYiIiIQIJXYiIiIiIUKJnYiIiEiIUGInIiIiEiKU2ImIiIiECCV2IiIiIiFCiZ2IiIhIiFBiJyIiIhIilNiJiIiIhAgldiIiIiIhQomdiIiISIhQYiciIiISIpTYiYiIiIQIJXYiIiIiISJoEjszG2Fma8zsiJmlVNLuOzP7yswyzSyjNmMUERERCaQGgQ6gGlYDlwEvV6Htec653TUcj4iIiEidEjSJnXNuHYCZBToUERERkTopaC7FVoMDPjSzFWZ2XaCDEREREaktdarHzswWAq3LqZrknJtdxd2c45zLNrOWwAIz+9o593E573UdcB1A27ZtjztmERERkbqiTiV2zrkL/LCPbO+/O83sXeB0oExi55ybAkwBSElJcSf6viIiIiKBFlKXYs0sxsziip8DA/EMuhAREREJeUGT2JnZMDPbApwFvGdmH3jLE8xsnrdZK+A/ZrYSWAa855x7PzARi4iIiNSuOnUptjLOuXeBd8sp3wpc7H3+LdCrlkMTERERqROCpsdORERERCqnxE5EREQkRCixExEREQkRSuxEREREQoQSOxEREZEQocROREREJEQosRMREREJEUrsREREREKEEjsRERGREKHETkRERCREKLETERERCRFK7ERERERChBI7ERERkRChxE5EREQkRCixExEREQkRSuxEREREQoQSOxEREZEQocROREREJEQosRMREREJEUGT2JnZ02b2tZmtMrN3zezkCtoNMrP1ZrbRzO6p5TBFREREAiZoEjtgAdDdOdcT2ADce3QDMwsHXgQGA12B0WbWtVajFBEREQmQoEnsnHMfOucKvZufA23KaXY6sNE5961zrgCYAQytrRhFREREAqlBoAM4TtcCM8spTwQ2l9jeApxR3g7M7DrgOu9mrpmt92uEdVdzYHegg6hjdE7K0jkpn85LWTon5dN5KUvnpHzHc15OqaiiTiV2ZrYQaF1O1STn3Gxvm0lAITDtRN7LOTcFmHIi+whGZpbhnEsJdBx1ic5JWTon5dN5KUvnpHw6L2XpnJTP3+elTiV2zrkLKqs3s2uAS4BfOedcOU2ygaQS2228ZSIiIiIhL2jusTOzQcDvgTTn3MEKmi0HOprZqWYWAYwC5tRWjCIiIiKBFDSJHfACEAcsMLNMM0sHMLMEM5sH4B1ccTPwAbAO+Idzbk2gAq6j6t3l5yrQOSlL56R8Oi9l6ZyUT+elLJ2T8vn1vFj5VzRFREREJNgEU4+diIiIiFRCiZ2IiIhIiFBiF0LMLMnMFpnZWjNbY2a3ecsfMrNs772JmWZ2cYnX3Otdfm29mV0UuOhrjpl9Z2ZfeY89w1vW1MwWmFmW998m3nIzs8nec7LKzPoENvqaYWanlfg8ZJpZjpndXh8/K2b2qpntNLPVJcqq/fkws3He9llmNi4Qx+IvFZyTcpd1NLN2ZpZX4jOTXuI1fb0/exu9580CcDh+UcE5qfbPS6gte1nBeZlZ4px8Z2aZ3vL68lmp6Hdx7XyvOOf0CJEHEA/08T6Pw7P0WlfgIeDOctp3BVYCkcCpwDdAeKCPowbOy3dA86PKngLu8T6/B3jS+/xiYD5gwJnAfwMdfy2cn3BgO54JL+vdZwXoB/QBVh/v5wNoCnzr/beJ93mTQB+bn8/JQKCB9/mTJc5Ju5LtjtrPMu95Mu95GxzoY/PzOanWz4v38Q3QHojwtuka6GPz93k5qv5Z4MF69lmp6HdxrXyvqMcuhDjntjnnvvA+349nZHBiJS8ZCsxwzh1yzm0CNuJZlq0+GAq87n3+OnBpifK/OY/PgZPNLD4A8dWmXwHfOOe+r6RNyH5WnHMfA3uPKq7u5+MiYIFzbq9z7kc8a1sPqvHga0h558RVbVlHH+95aeyc+9x5fkv9jZ/PY9Cp4HNSkYp+XkJu2cvKzou3120kML2yfYTgZ6Wi38W18r2ixC5EmVk7oDfwX2/Rzd4u3leLu38pfwm2yhLBYOWAD81shXmWkgNo5Zzb5n2+HWjlfV5fzklJoyj9xVufPyvFqvv5qG/n51o8PQzFTjWzL81siZmd6y1LxHMeioXqOanOz0t9+5ycC+xwzmWVKKtXn5WjfhfXyveKErsQZGaxwNvA7c65HOAl4BdAMrANT9d4fXKOc64PMBiYYGb9SlZ6/0Ksl/P+mGci7zTgn96i+v5ZKaM+fz7KY2WXddwGtHXO9QZ+B7xpZo0DFV8t089L5UZT+o/GevVZKed3sU9Nfq8osQsxZtYQzwdpmnPuHQDn3A7nXJFz7ggwlZ8vodWLJdicc9nef3cC7+I5/h3Fl1i9/+70Nq8X56SEwcAXzrkdoM9KCdX9fNSL82M/L+s41vuLCe/lxj3e5yvw3EPWCc/xl7xcG3Ln5Dh+XurF5wTAzBoAlwEzi8vq02elvN/F1NL3ihK7EOK9n+EVYJ1z7rkS5SXvERsGFI9emgOMMrNIMzsV6IjnBtaQYWYxZhZX/BzPDeCr8Rx78QijccBs7/M5wNXeUUpnAvtKdJ2HolJ/Udfnz8pRqvv5+AAYaGZNvJfjBnrLQoZVsKyjmbUws3Dv8/Z4Phvfes9Ljpmd6f1uupqfz2NIOI6fl/q07OUFwNfOOd8l1vryWanodzG19b1S06ND9KjVkTjn4OnaXQVkeh8XA38HvvKWzwHiS7xmEp6/mtYTxKOQKjkn7fGMPFsJrAEmecubAf8GsoCFQFNvuQEves/JV0BKoI+hBs9NDLAHOKlEWb37rOBJbLcBh/HcwzL+eD4feO472+h9/L9AH1cNnJONeO73Kf5uSfe2He792coEvgCGlNhPCp5k5xs8y0JaoI/Nz+ek2j8v3u/kDd66SYE+rpo4L97y14AbjmpbXz4rFf0urpXvFS0pJiIiIhIidClWREREJEQosRMREREJEUrsREREREKEEjsRERGREKHETkRERCREKLETERERCRFK7EREaoCZ/dvM/hboOESkflFiJyJSM/oAKwIdhIjUL0rsRET8zMx+AZyMEjsRqWVK7ERE/K8vcATPUkIiIrVGiZ2IiP/1BTY453IDHYiI1C9K7ERE/E/314lIQCixExHxPyV2IhIQSuxERPzIzE4FmgJfBDoWEal/GgQ6ABGRENPX++9hM+teorzIObcuEAGJSP2hxE5ExL+KE7ulR5WvBnrUciwiUs+Ycy7QMYiIiIiIH+geOxEREZEQocROREREJEQosRMREREJEUrsREREREKEEjsRERGREKHETkRERCREKLETERERCRFK7ERERERCxP8HkBcWm//2uP0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_spec(qobj,bnum=10,lmin=40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chi: 2.0 , chi (sim mean): 0.0 , chi (sim std): 5.9 , PTE: 0.70 \n",
      "chi^2: 5.5 , chi^2 (sim): 11.1 , PTE: 0.91 \n",
      "chi: 2.0 , chi (sim mean): 0.0 , chi (sim std): 5.9 , PTE: 0.70 \n",
      "chi^2: 4.1 , chi^2 (sim): 11.1 , PTE: 0.97 \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAEKCAYAAACfXBjeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAA0X0lEQVR4nO3deXxU1f3/8deHxABhKRAEsgphExIwTYIoIlhFFquAUBRBUEm/9vct1rZCkeq3FW2ttJWiVFqWutSlgEVUUGgFRdSCsmjYRAigIDsJKkLASDi/P2YyTSAJmZBkkpv38/G4j8yce+bmc2/uzHxy7rnnmHMOEREREan56oQ6ABERERGpGErsRERERDxCiZ2IiIiIRyixExEREfEIJXYiIiIiHhEe6gCqg+bNm7vWrVuHOgyRCnf69GmOHj1KkyZNADhy5AiNGzcmPFxvfRGRmmrdunXZzrkLi1unT3egdevWrF27NtRhiFSqPXv20Lp1axo2bMgvf/lL7r77burXrx/qsEREJEhmtqukdboUK1JLxMXF8dFHH9GzZ08mTpxI+/btefLJJzl16lSoQxMRkQqixE6kFunSpQuvvfYaK1asID4+nrvvvpvs7OxQhyUiIhVEiZ1ILdSrVy9WrlzJunXraNWqFc45MjIyWLBgAadPnw51eCIiUk5K7ERqKTPj4osvBuDw4cO89957DB06lNTUVF5++WUleCIiNZASOxGhRYsWbN68meeee47c3FyGDBlCWloa27dvD3VoIiISBCV2IgJAeHg4t956Kx9//DHPPvssjRs3JjY2FoBdu3apBU9EpAZQYiciRYSHhzNq1ChWrFhB/fr1OXXqFFdddRVJSUk888wz5OXlhTpEEREpgRI7ESmVmfHII49Qt25d7rjjDtq1a8fjjz/O8ePHQx2aiIicQYmdiJQqLCyM4cOH89FHH7F48WLatGnDz372M5YtWxbq0ERE5AxK7ESkTMyMAQMGsGLFCj744ANuuOEGACZPnsxdd91FVlZWiCMUEREldiIStEsvvZQ6dXwfHwcPHmT27Nl07NiRgQMH8tZbb+GcC3GEIiK1kxI7ETkvU6dOZdeuXfzqV7/i/fff55prruEXv/hFqMMSEamVlNiJyHlr1aoVDz74ILt37+bJJ5/klltuAWDr1q08+OCD7N+/P8QRiojUDkrsRKTC1KtXjzFjxpCWlgbAG2+8waRJk0hISOCmm25i+fLlukwrIlKJlNiJSKX5yU9+wtatW/nJT37CsmXLuPrqq+nWrZsGOxYJQlhYGCkpKSQnJzNs2DByc3PJyckhJSWFlJQUWrVqRWxsbOB5Xl5e4DUFy+TJk8/a7r/+9S86duxIu3btil2/devWItto3Lgxjz32WBXssZwP03/PkJ6e7tauXRvqMEQ87cSJE8ybN49Dhw4xYcIEnHM89NBD3HDDDaSmpoY6PJFqq2HDhhw7dgyAkSNHkpaWxj333BNYP2nSJBo2bMj48eOLfU1x8vPz6dChA0uXLiUuLo5u3boxZ84cOnfuXGL92NhYPvjgAy666KIK2jMpLzNb55xLL26dWuxEpErUr1+f22+/nQkTJgCwe/du/vCHP5CWlkb37t155plnNOix1CrHjx8nNjaW++67D4A1a9aQkpLCyZMnS3zNlVdeWSFzOK9evZp27dqRmJhIREQEw4cP59VXXy2x/ptvvknbtm2V1NUANS6xM7P+ZrbVzLab2cRS6g01M2dmxWa0IhJaF110EXv37mXatGkcPXqUO+64g+joaN5///1QhyZSJRo0aMCGDRuYM2cOubm5jBkzhmeeeYZ69eoVW//UqVMsWbKELl26nHPbJ06cKHIZdd68eQBcd9117Nu3j7179xIfHx+oHxcXx969e0vc3ty5cwM3RUn1Fh7qAIJhZmHAdOBaYA+wxswWOuc+PqNeI+CnwAdVH6WIlFWTJk34yU9+wl133cW7777L888/T9euXQF49tlnOXDgAKNGjSI6OjrEkYpUjqioKCIjI8nIyGDUqFGkpKScVacgSQNfi11GRsY5t1u/fn0yMzPPKl+8eHHQMebl5bFw4UIeeeSRoF8rVa+mtdhdCmx3zu10zuUBc4FBxdT7DfB7oOT2bBGpNsyMXr16MWvWLCIjIwFYvnw59957L/Hx8Vx//fW89NJL5OXlhThSkYrXtWtXDhw4EOgjN3369EBL2759+wJJWmZmJn/+85+JiIg4798ZGxvL559/Hni+Z88eYmNji627ZMkSUlNTadmy5Xn/Xql8NS2xiwU+L/R8j78swMxSgXjn3OulbcjM7jSztWa29vDhwxUfqYicl6effpqtW7cyYcIEPvroI37wgx8wYsSIUIclUqEOHz7M8uXLGTp0aGA2l7FjxwYSuZiYmEr5vd26dSMrK4tPP/2UvLw85s6dy8CBA4utO2fOHF2GrUFqWmJXKjOrA/wJGHeuus65Wc65dOdc+oUXXlj5wYlI0Dp06MDvfvc7du3axeLFi/nZz34GwN69ewNDOOzevTu0QYqch4yMDK6++mrWr19fods9s4/dxIm+LukFfezCw8N54okn6NevH506deKmm24iKSkp8PqCesePH2fp0qUMGTKkQuOTylOjhjsxs8uBSc65fv7nvwRwzj3if/4dYAdQcI93K+AIMNA5V+J4JhruRKRmyczM5Mc//jGrVq0CfP2ORowYwciRI2nUqFGIoxMpm5kzZ/LWW2/x6KOPMmDAADZt2hTqkKSG8NJwJ2uA9mbWxswigOHAwoKVzrmvnHPNnXOtnXOtgfc5R1InIjVPSkoKK1euZMeOHfz2t78lOzubsWPHBsbt+uyzz8jNzQ1xlCIly8rKYurUqcyYMYP4+Hiio6Pp27dvqMMSD6hRLXYAZnYd8BgQBjzlnHvYzB4C1jrnFp5R921g/LkSO7XYidRszjl27NhBu3btAOjbty8rV67kxhtvZNiwYfTt27fEISRERGqa0lrsalxiVxmU2Il4y4oVK3j++eeZP38+X375JY0aNeKee+5h0qRJoQ5NROS8eelSrIjIOfXu3ZvZs2dz8OBBlixZwk033RQYquHYsWOMGjWKBQsW6HKtiHiOWuxQi51IbbJ27VoGDBhAdnY2DRo04Pvf/z7Dhg3juuuuC4yhJxJKr7zyCq+//jpHjx4lIyNDfe/kLGqxExHxS09PZ//+/Sxbtoxbb72V5cuXM2zYMHbu3An4brw4ePBgiKOU2mzw4MHMnj2bGTNmBKYCK05YWBgpKSkkJyczbNiwwDBAKSkptGrVitjY2MDzvLy8QP2CZfLkycVu91//+hcdO3akXbt2JdZ5/PHHSU5OJikpicceeyxQPmbMGFq0aEFycvJ5HQM5D865Wr+kpaU5EamdTp065d599113+vRp55xzo0ePdmbmevTo4f7whz+4rVu3hjhCqa3uuecet27duhLXN2jQIPB4xIgRbsqUKYHnDzzwgPvjH/9YYv2SnDp1yiUmJrodO3a4b775xnXt2tVt3ry5SJ2NGze6pKQkd/z4cfftt9+6a665xmVlZTnnnFuxYoVbt26dS0pKKtM+Svngu2G02JxGLXYiUquFhYXRs2dPzAyA8ePH88ADD3DixAkmTJhAx44d6d+/f6C+U/cVqSDHjx8nNjaW++67D4A1a9aQkpLCiRMnuPfeexkwYACpqall2taVV17J9u3bzzum1atX065dOxITE4mIiGD48OG8+uqrReps2bKF7t27ExkZSXh4OL1792bBggUA9OrVi2bNmp13HFJ+SuxERArp0qULDzzwAB9++CG7du3iz3/+M4MHDwYgPz+fiy++mBEjRvD888+j6QjlfDRo0IANGzYwZ84ccnNzGTNmDM888wyzZ89m2bJlzJ8/nxkzZpxzO6dOnWLJkiV06dKl1HpnzkZRcJm3YJYJ8M3qEh8fH3hNXFwce/fuLbKd5ORk3n33XXJycsjNzWXx4sVF5p2V0AoPdQAiItVVQkICd911V+D5119/TY8ePViyZAlz5szBzLj00kt58MEH6devXwgjlZoqKiqKyMhIMjIyGDVqVCDpuvvuu8/52oJEDXwtdhkZGaXWr1+/PpmZmWeVL168OKiYO3XqxL333kvfvn1p0KABKSkphIWFBbUNqTxqsRMRKaMmTZrw9NNPs2/fPtauXcuDDz6ImQW+1NasWcPtt9/OvHnz+OKLL0IcrdQUXbt25cCBA4wfP77EOtOnTw8kfQWtawWJWmZmJn/+85+JiIg471hiY2OLtL7t2bOH2NjYs+plZGSwbt063nnnHZo2bUqHDh3O+3dLxVBiJyISpDp16pCWlsavfvUrVq1aRZ8+fQDYuXMnixYtYvjw4Vx44YVcccUVTJo0iaNHj4Y4YqmuDh8+zPLlyxk6dCh16pT8lTx27NhAEhcTE1Np8XTr1o2srCw+/fRT8vLymDt3LgMHDjyr3qFDhwDYvXs3CxYsYMSIEZUWkwRHiZ2ISAW5+eabOXToECtXrmTixIl8++23TJ06NTCd2TPPPMO0adPYsmWLbsIQwNfydfXVV7N+/fpK/11n9rGbOHEiULSPXXh4OE888QT9+vWjU6dO3HTTTSQlJZ1Vb+jQoXTu3JkbbriB6dOn06RJEwBuueUWLr/8crZu3UpcXBxPPvlkpe+XFKUBitEAxSJSeXJzcwMDHw8cOJBFixYBEBMTQ58+fRg0aBBDhgwJZYgSIjNnzuStt97i0UcfZcCAAWzatCnUIUkNoQGKRURCpPBsFgsXLmTnzp3MmjWLnj178vrrrzNnzpzA+l//+tfMnz8/cJlLvCsrK4upU6cyY8YM4uPjiY6O1gwTUiHUYoda7EQkNE6fPs1XX31F06ZNycnJISEhITB/7cUXX0zv3r2544476N69e4gjFZHqRC12IiLVUJ06dWjatCngG/biiy++YNWqVUyePJnExETmzJnDtm3bANi2bVtgnLNPP/1UffREpFhqsUMtdiJSPeXn55Ofn09ERARLlizh1ltv5ciRIwDEx8fTq1cvfve735GQkBDiSEWkKpXWYqcBikVEqqmwsLDAGHkDBgzg8OHDfPzxx6xYsYJ33nmH5cuX06hRIwD+8pe/sGzZMi6//HJ69OhBWlpa4G5cqVleeeUVXn/9dY4ePUpGRob63klQdClWRKSGqFOnDsnJyYwdO5Z58+axZ8+ewKXcEydOsGHDBiZMmEDPnj1p3Lgx3/ve9wKXbI8fPx7K0CUIgwcPZvbs2cyYMSMw7VdxwsLCSElJITk5mWHDhrF3797AUCatWrUiNjY28DwvLy9Qv2CZPHnyWds8ceIEvXv3Jj8//6x1kyZN4tFHHw0879GjR7n2Ly8vj169enHq1Kli1//rX/+iY8eOtGvXrtgYCzz++OMkJyeTlJTEY489BsDWrVuL7GPjxo0D62oLtdiJiNRQZhZ4PG7cOMaNG8ehQ4dYtWoVK1eu5Pjx44E611xzDfv37+fyyy+nW7dudOvWjdTUVBo2bBiq8OUcfvvb3zJ27NgS1xeeImzkyJHMmzcv8HzSpEk0bNiwyGwWJU0pVthTTz3FkCFDyjRF2MqVK89ZpzgRERFcc801zJs3j5EjRxZZl5+fz9ixY1m6dClxcXF069aNgQMH0rlz5yL1Nm3axOzZs1m9ejURERH079+f66+/no4dOwb2MT8/n9jYWG688cZyxVlTqcVORMRDWrRowaBBg/j973/PE088ESgfPXo0l112GatWrWL8+PH07t2b0aNHB9Y/9dRTfPDBB5w8eTIUYddKx48fJzY2lvvuuw/wTUmXkpLCiRMnuPfeexkwYACpqall2taVV17J9u3bzzumF154gUGDBgWeP/zww3To0IGePXuydevWInUL/ikYPHgwaWlpJCUlMWvWrCJ1fvOb39CxY0d69uzJLbfcEmjxGzx4MC+88MJZv3/16tW0a9eOxMREIiIiGD58OK+++upZ9bZs2UL37t2JjIwkPDyc3r17s2DBgiJ13nzzTdq2bctFF11UvoNRQymxExGpBX784x8zb948du3axcGDB3nttdf4+c9/DvimtcrIyOCyyy6jUaNGpKam8qMf/Yj//Oc/IY7a2xo0aMCGDRuYM2cOubm5gbueZ8+ezbJly5g/fz4zZsw453ZOnTrFkiVL6NKlS6n1zpx54szLvHl5eezcuZPWrVsDsG7dOubOnUtmZiaLFy9mzZo1xW73qaeeYt26daxdu5Zp06aRk5MD+BLVl156ifXr17NkyRIK36SYnJxcZHsFs1rs3buX+Pj4QHlcXBx79+4963cmJyfz7rvvkpOTQ25uLosXLy4yxy3A3LlzueWWW0o9Jl6kS7EiIrVMixYt+P73vx943rx5c3bv3s2aNWtYs2YNa9eu5cUXX6R79+5cccUVbNmyhYyMDLp160Z6ejrdunWjQ4cOpc5tKmUTFRVFZGQkGRkZjBo1KpB03X333ed8bUGiBr4Wu4yMjFLrn+tSbHZ2dmBqMIB3332XG2+8scjMKcWZNm0aL7/8MgCff/45WVlZREVF8Z///IdBgwZRr1496tWrxw033BB4TVhYGBEREXz99dc0atSIxYsXn3N/C+vUqRP33nsvffv2pUGDBqSkpBS5fJyXl8fChQt55JFHgtquFyixExGp5cyM+Ph44uPjA9ObOecCnduPHTtGeHg4Tz75JNOmTQOgUaNG/Pvf/+byyy9nz549HDx4kKSkJN2JWw5du3blwIEDRfrDnWn69OnMnj0bgMWLFxMTE1OmPnPBqF+/ftCX4t9++22WLVvGqlWriIyM5KqrrirzNr755puzzpfY2NgiLW979uwhNja22NdnZGQEktn77ruPuLi4wLolS5aQmppKy5Ytg9ofL9C/WyIichYz44ILLgCgW7duvPPOO3z11Vds2rSJp59+mtGjR5OYmAjAP/7xD9LT02nYsCHJycnceuut/PGPf9SduGVw+PBhli9fztChQ0ttAR07diyZmZlkZmYSExNTKbE0bdqU/Pz8QGLWq1cvXnnlFU6cOMHXX38dmOe4sIKZUyIjI/nkk094//33A+uuuOIKFi1axMmTJzl27BivvfZaYF1OTg7NmzcPnGMFunXrRlZWFp9++il5eXnMnTu3xJbCgqn3du/ezYIFCxgxYkRg3Zw5c2rlZVhQi52IiJRRWFgYSUlJJCUlcfvttwfKR44cSWJiIpmZmaxfv54VK1bw4osv8tOf/hTwdaBftWoVKSkpXHLJJVxyySW0a9eO8HB9BWVkZHD11Vezfv36Sv9dhS/dAvTv3/+s4UT69u3Le++9R58+fUhNTeXmm2/mkksuoUWLFnTr1u2sbfbv358ZM2bQqVMnOnbsyGWXXRZYV3BHa9euXWnZsiVdunThO9/5DgDLly8v0h3guuuu429/+xsxMTE88cQT9OvXj/z8fMaMGUNSUtJZdQCGDh1KTk4OF1xwAdOnTw9cRj5+/DhLly5l5syZFXLcahrNPIFmnhARqWhfffVV4Et88uTJzJkzh48//jhwebdNmzbs3LkTgAULFhAREUFycjIJCQm1pu/ezJkzeeutt3j00UcZMGAAmzZtCnVIfPjhh0ydOpXnnnuuQrZ37NgxGjZsSG5uLr169WLWrFmkpqYyZMgQJk+eTIcOHSrk99Q2mnlCRESqVEFSBzBx4kQmTpzIN998w5YtW1i/fj3ffPNNYP19990XGEqjYcOGdO7cmUGDBgWGAcnJyaFZs2ZFxu2r6bKyspg6dSqrVq2iadOmREdH07dvX954442QxpWamsr3vvc98vPzyzSW3bnceeedfPzxx5w8eZLbbruN1NRU8vLyGDx4sJK6SqIWO9RiJyISSl9++SUff/wxmzZtYvPmzWzatImUlBSmTJmCc47GjRtzwQUXBC4Dd+rUid69exe5rChSm6jFTkREqq0mTZrQo0ePYqeo+vbbb3nkkUcCSd+8efP48ssv+b//+z9SUlL48ssvueqqq7j44ouLLB07dqR+/foh2BuR0FJiJyIi1VZERAR33XVX4LlzjsOHDwcuyx49epSYmBjWrFnDiy++GJgbd+bMmdx5553s3LmTxx9/vEjCFx0d7anLuiKFKbETEZEaw8xo0aJF4HlCQkJgcNsTJ06wfft2Pvnkk8AdnDt27ODJJ58sMvRKwYC4PXv2ZPv27Xz00Ue0b9+e9u3b06BBg6rdIZEKpsROREQ8oX79+nTp0qXI1FrXXnstX3/9NXv37uWTTz7hk08+YevWrYFpswpPrQYQHR1N+/btmTNnDjExMXz66ad8/fXXtGvXLjADg0h1ppsn0M0TIiK11fHjx9m2bRvbt28nKysrsLzxxhtERkYyfvx4pkyZAkBMTAxt27albdu2zJ49m/DwcPbs2UPdunVp3ry5Lu9KlSnt5gkldiixExGR4m3fvp1169aRlZXFjh072LFjB1988QUbN24E4Oabb+bFF1+kUaNGJCYmkpiYSNeuXZk0aRLgm1miSZMmZ82wIHI+lNidgxI7EREpj/fee48PP/yQHTt2sHPnTnbs2MF3vvMdVq1aBfim1frggw9ISEggMTGRtm3b0r17d8aMGQPAgQMHiIqKUuInQfFUYmdm/YHHgTDgb865yWesvwf4IXAKOAyMcc7tKm2bSuxERKSiOOcCl2X/+c9/smHDhkDSt2PHDi6//HIWLlwIQHx8PPv27SMmJob4+HgSEhLo168fd9xxBwCbN28mJiaGJk2a6FKvBHhmHDszCwOmA9cCe4A1ZrbQOfdxoWofAenOuVwz+1/gD8DNVR+tiIjURoUTsGHDhjFs2LAi6/Pz8wOPH3roIT777DN2797N7t27Wbt2La1atQLg1KlTXHLJJeTn59OwYUMSEhJISEjg1ltvZeTIkeTn5/Pee++RkJBAbGwsERERVbODUq3VqMQOuBTY7pzbCWBmc4FBQCCxc84tL1T/feDWKo1QRESkFIWn6ipomSuOc4558+YFkr5du3axe/dujhw5AsC+ffu46qqrAF8yGR0dTUJCAhMmTODGG2/k6NGjLF++PJAQem1aNileTUvsYoHPCz3fA3QvpX4GsKRSIxIREakEF1xwAUOHDi1xfVRUFEuXLg0kfgVLQeK4adMmBg8eHKjfoEEDEhISmDZtGn369OHzzz8vkvjFxcWp1c8DalpiV2ZmdiuQDvQuYf2dwJ3gG+BSRESkJomMjKRPnz4lrk9JSWHNmjVnJX5NmzYFYNWqVdx2222B+mZGq1atWLRoEWlpaaxfv54VK1YQExMTWKKjo6lbt26l75uUX01L7PYC8YWex/nLijCzPsD9QG/n3DfFbcg5NwuYBb6bJyo+VBERkdCJjIwkPT2d9PRi+9gzaNAgtm3bdlbiV9DH780332TcuHFnvW7nzp20adOG+fPns3jx4iKJX0xMDKmpqYSH17T0wjtq2pFfA7Q3szb4ErrhwIjCFczsu8BMoL9z7lDVhygiIlL91a1bNzCVWnF+/vOfM3r0aPbt21dkiY6OBmD37t288cYbHDhwoMgNISdPniQ8PJwHHniARYsWFUn64uLi+OEPfwjAF198Qf369alXr17l72wtUhOHO7kOeAzfcCdPOeceNrOHgLXOuYVmtgzoAuz3v2S3c25gadvUcCciIiLlk5+fz+HDh9m3bx8HDx5kwIABAMycOZNFixYFEsJDhw7RvHlzDh3ytbkMGTKEl19+mcaNG9OiRQtatmxJ586dmTVrFgBLlizh5MmTtGzZMrA0bNgwZPtZnXhqHLvKoMRORGqykSOPAfDCC/rSk+rr22+/5ciRI7Rs2RLwzdO7YcMGDh48GFiioqKYP38+AJdeeilr1qwpso0rr7ySd955B/C1KB47dqxI4tehQwdSUlKAouMJeo1nxrETERGRmumCCy4IJHUA119/Pddff32J9RcuXMj+/fvPSvwKbN68mQ0bNnD48GFOnz4NwODBg3n55ZcB342Rzrkiid+1117LiBG+HlwrVqwgKiqKli1bEhUVRZ06dSpjt6tcmRM7M7vUObe6MoPxOv1XLSIiUjatWrUK3MhRnDfeeAPwXQrOycnh4MGDRaZmK+gfWJAUbty4kWbNmjFixAjy8vICYwAC1KlThwsvvJBx48bxi1/8gtzcXCZOnEjz5s2JioqiefPmNG/enE6dOhETE1Np+1wRgmmxe8XM/gHkAeuAfzvnjlVOWCIiIiLnFhYWRosWLWjRokWR8ocffrjE19SpU4e33367SGvgwYMHSUxMBODIkSM899xzfPnll0VeN2XKFO655x62bdtGWlpakaSvefPm/OhHP+LKK68kOzub5cuXF1kXFRVVJeMEBpPYOWA3EAEMBv5oZv/rnPt3ZQQmIiIiUhnCw8Pp3bvYYW4BiIuL44svvuDUqVMcOXKE7OxssrOzueiiiwDfYM//8z//Q3Z2Njk5OWRnZ5OVlRUYUHrjxo3cdNNNZ2331VdfZeDAgbz//vs8+OCDRRK/5s2bM2jQIFq1asXXX3/NsWPHaNasWdDjBgaT2B1yzk0reGJm0cBiQImdiIiIeE54eHixrYGxsbH86U9/KvF13bt3Z+PGjYGEsGBJSkoCIDc3l+zsbD755BOys7M5dsx3AfS73/0urVq14p///CcZGRmAL4ls1qwZzZo146WXXqJt27alxxzE/jkAM7sdaAQ0BrzR01BERESkgkRGRpKcnFzi+quvvrrIHb/ffPMNOTk5gZtDevTowV//+leOHDlCTk5O4GdZhnspz12xm4EB/uUX5Xi9iIiIiPjVrVu3yE0ZF198MRdffHG5thVMYtfYzKKcc2uANWb2e+B94I1y/WYREZFz0GgCIsEJJrFrDqw1s53AFqAJ4M2R/0RERKqAElepaMEkdgPwzdV6JdAeOA2cPTuwiIiIiIREMIndS4DGsZNi6b9OERGR0AvmrtaCceyO4BvHboOZ9auMoEREREQkeBrHTkRERKQcquPVKo1jJyIiIuIR5UnMNgNNgespwzh2ZjavHL9DRERERIIUTGIXGMfOOfcQcBXw+zK8rlu5IvOYhQvzyMzMZ/XqfHr3PsrChXmhDklEREQ8JphLsVGUMI6dmf0W2ORfPnHOnargOGu0hQvzuP/+E+T5c7l9+xz3338CgIEDI0IYmYiIiHhJMInd1cBGih/HLge4Fvg50N7M9uG7ZLsJqD49CkNkypSTnDxZtOzkSV+5EjsRERGpKGVO7JxzH/kfLvcvhddNLfzczNoAyf5l6XnGWOPt3++CKhcRKauCbh55edC791HGjaunfxhFarFgWuwws3TgR8AvnHNfllTPOfcp8Cmw6Lyi84joaGPfvrOTuOhozcgmIuWnbh4icqZg74qdANQtLqkzs/pm1sPM9GlyhnHj6lGvXtGyevV85SJS+UaOPBYYb8pLSuvmIVIdePW9V50Fm9ilA7OKW+GcOwFMAn55njF5zsCBETz8cH0i/ClvTIzx8MP19R+1iJwXdfMQkTMFm9hFA5+Vsv5v+Ma3kzMMHBhBSkoYl14axooVjZXUSbWi/6prppK6c6ibh0jtFWxilw3El7J+HZBY/nBERKSs1M1DRM4UbGK3AhhTynoD1BQlIlIF1M1DRM4UbGL3KHC7md1dwvorgB3nF5LUNJpVQyR0vNzNQ58tIsELKrFzzmUC/wNMMbN3zGyYmSWaWQszG4ZvirHnKyFOqaZKGm5BH8Aicj5qw2eLElepDMG22OGcewbohW8MvHlAFrDf/3gd8HgFxifVnIZbEJHK4PXPltqQuEpoBDVAcQHn3Cqgh5m1A1KBSGCTc25tRQYn1Z+GWxCRyuD1zxZNNSmVpVyJXQHn3HZgewXFIjWQZtUQkcrg9c8WryeuEjpBX4o1s2QzG21maZURkNQsGm5BRCqD1z9bNAahVJZzJnZm9mahxyOAuUAyMNPM7qrE2KQG0HALIlIZvP7Z4vXEtTaorje/lOVSbLNCj38K9HHOHTCzhsBK4IlKiUxqjIEDI5g3z3dCv/BCwxBHIyJe4eXPloIE9Ze/9N1AERNjjBtXzzOJq9eVdPMLEPK/YVkuxZqZ1TezBkAd59wBAOfcMSC/UqMTERHxKC+PQeh11fmu7bK02DUBNuObVcKZWbRzbr+/xU6dAYLgtf84RUREaqPqfPPLOVvsnHOtnXOJzrk2/p/7/atOAzdWbnhnM7P+ZrbVzLab2cRi1tc1s3n+9R+YWeuqjlFERES8qzrf/BL0XbEFnHO5zrlPKzKYczGzMGA6MADoDNxiZp3PqJYBfOGcawdMxTcbhoiIiEiFqM43v5hz5W82NLMPnXOpFRjPuX7f5cAk51w///NfAjjnHilU59/+OqvMLBw4AFzoStnRRo0aubQ0jd5yPrZs8XW37NQpLMSRSHl4+e93/HgvjhwZhVlzwsOzadToORo0eCfUYVUoL//9vLxv4O398/p7L5T7t2LFinXOufTi1p3vAMVVltT5xQKfF3q+B+heUh3n3Ckz+wqIArILVzKzO4E7AerWrVtZ8YpICB0/3ouvvhpLnTq+/6Lz81vw1VdjATz1BSNS3dSG916DBu+we/dyoHol5mVK7MysKdAXX9IEsA/4t3Pui8oKrLI552YBswDS09Pd22+/HdqAariRI48BukGkpvLq369376N8+WXRxnrn6hEZOZ63334oRFFVPK/+/cDb+wbe3T+99yqXWcl9+coyQHEGsApfy1gd/9IdWOlfV5X2AvGFnsf5y4qt478U+x0gp0qiE88aOfJY4A0sNUd1vnNNxMv03gudsrTYTQBSnXPHCxea2a+AD4EnKyOwEqwB2ptZG3wJ3HBgxBl1FgK34UtGfwC8VVr/OhHxLq/PNypSXem9FzpluSvWAY2KKW/kX1dlnHOngLuAfwNbgBedc5vN7CEzG+iv9iQQZWbbgXuAs4ZEEZH/qq7T4lSE6nznmoiX6b0XOmVpsRsPrDCzTfz3smcckASMq6zASuKcWwwsPqPs14UenwSGVXVcIjVRdZ4WpyLUlmmbvNY/S2q+2vLeq47Omdg5514zsyXApUCMv3gfsNo5pynFRGqw0qbF8coHsJfnGxWpzvTeC40y3RXrT+BWVXIsIlLF1MFZRMRbzmscO5EC+m+sZlIHZxERbyn3lGLFMbOLzWx9RW5TRCqPOjiLiHhLRbfYhQHJFbxNEakk6uAsIuItuhQrUsupg7NUZzonRYITVGJnZs8A6/xLpnMutzKCEhERqQ2UuEpFC7bFLh4YCDQB8s1sG/9N9D4ETlVodCIiIiJSZkElds65awDMLBFI8y+pwK+AZgXVKjJAERERESmbcvWxc87tBHYC/ywoM7PWQDq+RE9EREREqliF3TzhnPsM+AyYX1HbFBEREZGyK9M4dmZW38xiiylPqviQRERERKQ8zpnYmdkPgCzgdTPbYGbdC61+rtIiExEREZGglOVS7P8Bac65g2aWBvzdzH7nnPsHoHmHREREpFaqjsPVlCWxu8A5dxDAObfOzHoBL5tZO3QHrIiIiEi1UZY+dofMrGvBE+fcEeBaoBPQtcRXiXjEwoV5ZGbms3p1Pr17H2XhwrxQhyQiIlKssiR2o4BDhQucc3nOuVuA3pUSlUg1sXBhHvff75tHFWDfPsf9959QciciItXSOS/FOuf2lLLuPxUbjkj1MmXKSU6eLFp28qSvfODAiNAEJUGrjv1gRGoDvfeqXtDj2JlZHPC/QA+glb/4APAfYKZz7vOKC08ktPbvL74baUnlIiIioVSmcewKmFlPYAswDNgM/MO/bC4oM7MrKjpIkVCJji7+xu+SykVEREIp2Ba7x4CnnXN3F7fSzB731+l2fmGJVA/jxtXj/vtPFLkcW6+er1xERKS6CarFDkgCppey/q9AcvnDEaleBg6M4OGH6xPh704XE2M8/HB99a8TEZFqKdgWu/3AFcDWEtZf4a8j4hkDB0Ywb57vLlh1BBYRkeos2MTuUWCGmV0KLAUO+stb4hvb7nbgZxUVnIiIiIiUXVCJnXPuL2aWA/wcyADC/KvygXXAaOfcixUboohUNrVEioh4Q9DDnTjn5gHzzOwCoLm/ONs5922FRiYiIiIiQQn25okA59y3zrn9/uVbADOLN7OnKi48ERERESmrcid2JWgG3FbB2xQRERGRMgjqUqyZjT5HlYTziEVEREREzkOwfeyeAXKBkuZTqugWQBEREREpo2ATsX347nxtVNyCbxw7EREREQmBYBO7dUBqKesdoEk0RUREREKgPAMUlzbg1Xbge+UPR0RERETKK9gBit89x/rjwIrzikhEREREyqXG3OxgZs3MbKmZZfl/Ni2mToqZrTKzzWa2wcxuDkWsIiIiIqFQYxI7YCLwpnOuPfCm//mZcvHd3JEE9AceM7MmVReiiIiISOjUpMRuEPB3/+O/A4PPrOCc2+acy/I/3gccAi6sqgBFREREQqkmJXYtnXP7/Y8PAC1Lq2xmlwIRwI7KDkxERESkOgj2rthKZWbLgFbFrLq/8BPnnDOzkgZJxsyigeeA25xzp0uocydwJ0BCgibMEBERkZqvWiV2zrk+Ja0zs4NmFu2c2+9P3A6VUK8x8Dpwv3Pu/VJ+1yxgFkB6enqJSaKIiIhITVGTLsUuBG7zP74NePXMCmYWAbwMPOucm1+FsYmIiIiEXLVqsTuHycCLZpYB7AJuAjCzdOD/Oed+6C/rBUSZ2e3+193unMus+nDFS154obRxuUVERKoHc05XIdPT093atWtDHYaIiIjIOZnZOudcenHratKlWBEREREphRI7EREREY9QYiciIiLiEUrsRERERDxCiZ2IiIiIRyixExEREfEIJXYiIiIiHqHETkRERMQjlNiJiIiIeIQSOxERERGPUGInIiIi4hFK7EREREQ8QomdiIiIiEcosRMRERHxCCV2IiIiIh6hxE5ERETEI5TYiYiIiHiEEjsRERERj1BiJyIiIuIRSuxEREREPEKJnYiIiIhHKLETERER8QgldiIiIiIeocRORERExCOU2ImIiIh4hBI7EREREY9QYiciIiLiEUrsRERERDxCiZ2IiIiIRyixExEREfEIJXYiIiIiHqHETkRERMQjlNiJiIiIeIQSOxERERGPUGInIiIi4hE1JrEzs2ZmttTMsvw/m5ZSt7GZ7TGzJ6oyRhEREZFQqjGJHTAReNM51x540/+8JL8B3qmSqERERESqiZqU2A0C/u5//HdgcHGVzCwNaAm8UTVhiYiIiFQPNSmxa+mc2+9/fABf8laEmdUBpgDjz7UxM7vTzNaa2drDhw9XbKQiIiIiIRAe6gAKM7NlQKtiVt1f+IlzzpmZK6bej4HFzrk9Zlbq73LOzQJmAaSnpxe3LREREZEapVolds65PiWtM7ODZhbtnNtvZtHAoWKqXQ5caWY/BhoCEWZ2zDlXWn88EREREU+oVondOSwEbgMm+3++emYF59zIgsdmdjuQrqROREREaoua1MduMnCtmWUBffzPMbN0M/tbSCMTERERqQbMOXUvS09Pd2vXrg11GCIiIiLnZGbrnHPpxa2rSS12IiIiIlIKJXYiIiIiHqHETkRERMQjlNiJiIiIeIQSOxERERGPUGInIiIi4hFK7EREREQ8QomdiIiIiEcosRMRERHxCCV2IiIiIh6hxE5ERETEI5TYiYiIiHiEEjsRERERj1BiJyIiIuIRSuxEREREPEKJnYiIiIhHmHMu1DGEnJkdBnaFOo4q0BzIDnUQ1ZCOy9l0TIqn43I2HZPi6bicTcekeOU5Lhc55y4sboUSu1rEzNY659JDHUd1o+NyNh2T4um4nE3HpHg6LmfTMSleRR8XXYoVERER8QgldiIiIiIeocSudpkV6gCqKR2Xs+mYFE/H5Ww6JsXTcTmbjknxKvS4qI+diIiIiEeoxU5ERETEI5TYiYiIiHiEEjsPMbN4M1tuZh+b2WYz+6m/fJKZ7TWzTP9yXaHX/NLMtpvZVjPrF7roK4+ZfWZmG/37vtZf1szMlppZlv9nU3+5mdk0/zHZYGapoY2+cphZx0LnQ6aZHTWzn9W2c8XMnjKzQ2a2qVBZ0OeGmd3mr59lZreFYl8qUgnH5Y9m9ol/3182syb+8tZmdqLQOTOj0GvS/O+97f5jZyHYnQpRwjEJ+v1iZv39ZdvNbGJV70dFK+G4zCt0TD4zs0x/eW05V0r6Lq6azxbnnBaPLEA0kOp/3AjYBnQGJgHji6nfGVgP1AXaADuAsFDvRyUcl8+A5meU/QGY6H88Efi9//F1wBLAgMuAD0IdfxUcnzDgAHBRbTtXgF5AKrCpvOcG0AzY6f/Z1P+4aaj3rRKOS18g3P/494WOS+vC9c7Yzmr/sTL/sRsQ6n2r4GMS1PvFv+wAEoEIf53Ood63ij4uZ6yfAvy6lp0rJX0XV8lni1rsPMQ5t98596H/8dfAFiC2lJcMAuY6575xzn0KbAcurfxIq4VBwN/9j/8ODC5U/qzzeR9oYmbRIYivKl0D7HDOlTb7iifPFefcO8CRM4qDPTf6AUudc0ecc18AS4H+lR58JSruuDjn3nDOnfI/fR+IK20b/mPT2Dn3vvN9Sz3Lf49ljVPCuVKSkt4vlwLbnXM7nXN5wFx/3RqrtOPib3W7CZhT2jY8eK6U9F1cJZ8tSuw8ysxaA98FPvAX3eVv4n2qoPkX34n2eaGX7aH0RLCmcsAbZrbOzO70l7V0zu33Pz4AtPQ/ri3HpLDhFP3grc3nCgR/btSmY1NgDL4WhgJtzOwjM1thZlf6y2LxHYsCXj0uwbxfatu5ciVw0DmXVaisVp0rZ3wXV8lnixI7DzKzhsBLwM+cc0eBvwJtgRRgP76m8dqkp3MuFRgAjDWzXoVX+v9DrJXj/phZBDAQ+Ke/qLafK0XU5nOjJGZ2P3AKeMFftB9IcM59F7gH+IeZNQ5VfFVM75fS3ULRfxpr1blSzHdxQGV+tiix8xgzuwDfifSCc24BgHPuoHMu3zl3GpjNfy+h7QXiC708zl/mKc65vf6fh4CX8e3/wYJLrP6fh/zVa8UxKWQA8KFz7iDoXPEL9tyoNcfGzG4HrgdG+r+Y8F9uzPE/XoevD1kHfMeg8OVazx2XcrxfatO5Eg4MAeYVlNWmc6W472Kq6LNFiZ2H+PszPAlscc79qVB54T5iNwIFdy8tBIabWV0zawO0x9eB1TPMrIGZNSp4jK8D+CZ8+15wh9FtwKv+xwuB0f67lC4DvirUdO5FRf6jrs3nSiHBnhv/BvqaWVP/pbi+/jJPMbP+wARgoHMut1D5hWYW5n+ciO/c2Ok/NkfN7DL/Z9No/nssPaEc75c1QHsza+NvLR/ur+tFfYBPnHOBS6y15Vwp6buYqvpsqey7Q7RU6Z04PfE17W4AMv3LdcBzwEZ/+UIgutBr7sf3X9NWavBdSKUck0R8d56tBzYD9/vLo4A3gSxgGdDMX27AdP8x2Qikh3ofKvHYNABygO8UKqtV5wq+pHY/8C2+/isZ5Tk38PU52+5f7gj1flXScdmOr79PwWfLDH/dof73VibwIXBDoe2k40t2dgBP4J/tqCYuJRyToN8v/s/kbf5194d6vyrjuPjLnwH+3xl1a8u5UtJ3cZV8tmhKMRERERGP0KVYEREREY9QYiciIiLiEUrsRERERDxCiZ2IiIiIRyixExEREfEIJXYiIiIiHqHETkSkEpjZm2b2bKjjEJHaRYmdiEjlSAXWhToIEaldlNiJiFQwM2sLNEGJnYhUMSV2IiIVLw04jW8qIRGRKqPETkSk4qUB25xzx0IdiIjULkrsREQqnvrXiUhIKLETEal4SuxEJCSU2ImIVCAzawM0Az4MdSwiUvuEhzoAERGPSfP//NbMkguV5zvntoQiIBGpPZTYiYhUrILE7j9nlG8CulRxLCJSy5hzLtQxiIiIiEgFUB87EREREY9QYiciIiLiEUrsRERERDxCiZ2IiIiIRyixExEREfEIJXYiIiIiHqHETkRERMQjlNiJiIiIeMT/B9HFY809DfUDAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_spec(robj,bnum=10,lmin=40,ylim=[-.5,.5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py4so",
   "language": "python",
   "name": "py4so"
  },
  "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}