optimized_cicle.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8899ca7e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "%pylab is deprecated, use %matplotlib inline and import the required libraries.\n",
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "from scipy import optimize\n",
    "rc('font',**{'family':'sans-serif','sans-serif':['Latin Modern']})\n",
    "rc('text', usetex=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ffeb5148",
   "metadata": {},
   "outputs": [],
   "source": [
    "datafile=\"g-10.00_h20.00/lapse_5.53e-01.dat\" # this one has omega = 7.044948446677619103e-01\n",
    "\n",
    "### OJO #################################################\n",
    "#omegaB=0.82666    #DONT FORGET CHANGING THE VALUE OF OMEGA, g AND h HERE\n",
    "\n",
    "omegaB=7.044948446677619103e-01\n",
    "g=-10.0\n",
    "h=20.0\n",
    "#########################################################\n",
    "\n",
    "import numpy.polynomial.chebyshev as pch\n",
    "\n",
    "N=18; M=int(N/2); N0=N\n",
    "r_m1=0.1; r_m2=2.5; r_m2m=5.0; r_m3=25.0\n",
    "\n",
    "G=np.loadtxt(datafile); M0=int(N0/2)\n",
    "u10=np.zeros(N0+1); v10=np.zeros(N0+1); w10=np.zeros(N0+1); z10=np.zeros(N0+1)\n",
    "i1,i2=0,M0+1\n",
    "for ii in np.arange(M0+1):\n",
    "    u10[2*ii]=G[ii+i1]\n",
    "i1,i2=i2,i2+N0+1; u20=G[i1:i2]; i1,i2=i2,i2+N0+1; u20m=G[i1:i2]; i1,i2=i2,i2+N0+1; u30=G[i1:i2]; i1,i2=i2,i2+N0+1; u40=G[i1:i2]\n",
    "i1,i2=i2,i2+M0+1\n",
    "for ii in np.arange(M0+1):\n",
    "    v10[2*ii]=G[ii+i1]\n",
    "i1,i2=i2,i2+N0+1; v20=G[i1:i2]; i1,i2=i2,i2+N0+1; v20m=G[i1:i2]; i1,i2=i2,i2+N0+1; v30=G[i1:i2]; i1,i2=i2,i2+N0+1; v40=G[i1:i2]\n",
    "i1,i2=i2,i2+M0+1\n",
    "for ii in np.arange(M0+1):\n",
    "    w10[2*ii]=G[ii+i1]\n",
    "i1,i2=i2,i2+N0+1; w20=G[i1:i2]; i1,i2=i2,i2+N0+1; w20m=G[i1:i2]; i1,i2=i2,i2+N0+1; w30=G[i1:i2]; i1,i2=i2,i2+N0+1; w40=G[i1:i2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b7a55d95",
   "metadata": {},
   "outputs": [],
   "source": [
    "if omegaB!= G[-1]:\n",
    "    print(\"STOP!!!! ! ! !\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ad0c6d54",
   "metadata": {},
   "outputs": [],
   "source": [
    "def phi_interpol(ri):\n",
    "    if ri<r_m1:\n",
    "        x=ri/r_m1\n",
    "        val=pch.Chebyshev(u10)(x)\n",
    "    elif ri>=r_m1 and ri<r_m2:\n",
    "        x=(2*ri-r_m2-r_m1)/(r_m2-r_m1)\n",
    "        val=pch.Chebyshev(u20)(x)\n",
    "    elif ri>=r_m2 and ri<r_m2m:\n",
    "        x=(2*ri-r_m2m-r_m2)/(r_m2m-r_m2)\n",
    "        val=pch.Chebyshev(u20m)(x)\n",
    "    elif ri>=r_m2m and ri<r_m3:\n",
    "        x=(2*ri-r_m3-r_m2m)/(r_m3-r_m2m)\n",
    "        val=pch.Chebyshev(u30)(x)\n",
    "    else:\n",
    "        x=1.0-2.0*r_m3/ri\n",
    "        val=pch.Chebyshev(u40)(x)\n",
    "    return val\n",
    "def N_interpol(ri):\n",
    "    if ri<r_m1:\n",
    "        x=ri/r_m1\n",
    "        val=pch.Chebyshev(v10)(x)\n",
    "    elif ri>=r_m1 and ri<r_m2:\n",
    "        x=(2*ri-r_m2-r_m1)/(r_m2-r_m1)\n",
    "        val=pch.Chebyshev(v20)(x)\n",
    "    elif ri>=r_m2 and ri<r_m2m:\n",
    "        x=(2*ri-r_m2m-r_m2)/(r_m2m-r_m2)\n",
    "        val=pch.Chebyshev(v20m)(x)\n",
    "    elif ri>=r_m2m and ri<r_m3:\n",
    "        x=(2*ri-r_m3-r_m2m)/(r_m3-r_m2m)\n",
    "        val=pch.Chebyshev(v30)(x)\n",
    "    else:\n",
    "        x=1.0-2.0*r_m3/ri\n",
    "        val=pch.Chebyshev(v40)(x)\n",
    "    return val\n",
    "def Psi_interpol(ri):\n",
    "    if ri<r_m1:\n",
    "        x=ri/r_m1\n",
    "        val=pch.Chebyshev(w10)(x)\n",
    "    elif ri>=r_m1 and ri<r_m2:\n",
    "        x=(2*ri-r_m2-r_m1)/(r_m2-r_m1)\n",
    "        val=pch.Chebyshev(w20)(x)\n",
    "    elif ri>=r_m2 and ri<r_m2m:\n",
    "        x=(2*ri-r_m2m-r_m2)/(r_m2m-r_m2)\n",
    "        val=pch.Chebyshev(w20m)(x)\n",
    "    elif ri>=r_m2m and ri<r_m3:\n",
    "        x=(2*ri-r_m3-r_m2m)/(r_m3-r_m2m)\n",
    "        val=pch.Chebyshev(w30)(x)\n",
    "    else:\n",
    "        x=1.0-2.0*r_m3/ri\n",
    "        val=pch.Chebyshev(w40)(x)\n",
    "    return val\n",
    "\n",
    "def dN_interpol(ri):\n",
    "    if ri<r_m1:\n",
    "        x=ri/r_m1\n",
    "        val=pch.Chebyshev.deriv(pch.Chebyshev(v10),m=1)(x)/r_m1\n",
    "    elif ri>=r_m1 and ri<r_m2:\n",
    "        x=(2*ri-r_m2-r_m1)/(r_m2-r_m1)\n",
    "        val=pch.Chebyshev.deriv(pch.Chebyshev(v20),m=1)(x)/(0.5*(r_m2-r_m1))\n",
    "    elif ri>=r_m2 and ri<r_m2m:\n",
    "        x=(2*ri-r_m2m-r_m2)/(r_m2m-r_m2)\n",
    "        val=pch.Chebyshev.deriv(pch.Chebyshev(v20m),m=1)(x)/(0.5*(r_m2m-r_m2))\n",
    "    elif ri>=r_m2m and ri<r_m3:\n",
    "        x=(2*ri-r_m3-r_m2m)/(r_m3-r_m2m)\n",
    "        val=pch.Chebyshev.deriv(pch.Chebyshev(v30),m=1)(x)/(0.5*(r_m3-r_m2m))\n",
    "    else:\n",
    "        x=1.0-2.0*r_m3/ri\n",
    "        val=pch.Chebyshev.deriv(pch.Chebyshev(v40),m=1)(x)*(2.0*r_m3)/(2*r_m3/(1-x))**2\n",
    "    return val\n",
    "def dPsi_interpol(ri):\n",
    "    if ri<r_m1:\n",
    "        x=ri/r_m1\n",
    "        val=pch.Chebyshev.deriv(pch.Chebyshev(w10),m=1)(x)/r_m1\n",
    "    elif ri>=r_m1 and ri<r_m2:\n",
    "        x=(2*ri-r_m2-r_m1)/(r_m2-r_m1)\n",
    "        val=pch.Chebyshev.deriv(pch.Chebyshev(w20),m=1)(x)/(0.5*(r_m2-r_m1))\n",
    "    elif ri>=r_m2 and ri<r_m2m:\n",
    "        x=(2*ri-r_m2m-r_m2)/(r_m2m-r_m2)\n",
    "        val=pch.Chebyshev.deriv(pch.Chebyshev(w20m),m=1)(x)/(0.5*(r_m2m-r_m2))\n",
    "    elif ri>=r_m2m and ri<r_m3:\n",
    "        x=(2*ri-r_m3-r_m2m)/(r_m3-r_m2m)\n",
    "        val=pch.Chebyshev.deriv(pch.Chebyshev(w30),m=1)(x)/(0.5*(r_m3-r_m2m))\n",
    "    else:\n",
    "        x=1.0-2.0*r_m3/ri\n",
    "        val=pch.Chebyshev.deriv(pch.Chebyshev(w40),m=1)(x)*(2.0*r_m3)/(2*r_m3/(1-x))**2\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5abda7e7",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d85d33e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "xf=480/2\n",
    "num=4*2400+1\n",
    "\n",
    "xx=np.linspace(0,xf,num)\n",
    "xx_double=np.linspace(0,xf,2*(num)-1)\n",
    "\n",
    "phi_array=zeros(2*(num)-1)\n",
    "for oo in arange(2*(num)-1):\n",
    "    phi_array[oo]=phi_interpol(xx_double[oo])\n",
    "metric_term_array=zeros(2*(num)-1)\n",
    "for oo in arange(2*(num)-1):\n",
    "    metric_term_array[oo]=dN_interpol(xx_double[oo])/N_interpol(xx_double[oo])+2.0*dPsi_interpol(xx_double[oo])/Psi_interpol(xx_double[oo]) \n",
    "metric_term2_array=zeros(2*(num)-1)\n",
    "for oo in arange(2*(num)-1):\n",
    "    metric_term2_array[oo]=Psi_interpol(xx_double[oo])**4  \n",
    "metric_term3_array=zeros(2*(num)-1)\n",
    "for oo in arange(2*(num)-1):\n",
    "    metric_term3_array[oo]=N_interpol(xx_double[oo])**2  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "140cb899",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 15.0)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGbCAYAAABklPKCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAvLUlEQVR4nO3deXzcd33n8bek0YxujUaWb/kYOXbixJjoSAzhCCCnBUJLWJtQWloKa/nRY2G3tFa9tMsjD9qqdg8WemzttOzSgxJLQHcpbcFDCCQhcWSNc+DEjq2fD/mSLY1mRveMZn77x0hjSZbs0TW/OV7Px2MeM/PTz85Hv0i/eft75pimaQoAACBJcq0uAAAAZBfCBwAASCrCBwAASCrCBwAASCrCBwAASCrCBwAASCrCBwAASCqb1QVMF41GdeXKFZWWlionJ8fqcgAAQAJM01R/f79Wr16t3Nzbt22kXPi4cuWKqqurrS4DAADMQ1dXl9auXXvbc1IufJSWlkqKFV9WVmZxNQAAIBHBYFDV1dXxz/HbSbnwMdHVUlZWRvgAACDNJDJkggGnAAAgqQgfAAAgqQgfAAAgqQgfAAAgqQgfAAAgqQgfAAAgqQgfAAAgqQgfAAAgqeYUPrxer+rq6u54nmEYOnjwoNra2nTw4EH5/f751gcAADJMwiuctrW1ye12y+v13vHc3bt3q6OjQ1IsiOzZs0etra3zrxIAAGSMhMPHrl27EjrPMIwp791utzwez9yqAgAAGWvRx3x4PB65XK4px1wuV0ItJgAAIPMt+sZys43v8Pl8Mx4fHR3V6Oho/H0wGFzskgAAQApJ2myX2UJJS0uLysvL44/q6upklQQAACyw6OHD6XTe0srh8/nkdDpnPH///v0KBALxR1dX12KXBAAAUsiih4/GxsYZj9fX18943OFwqKysbMoDAABkrnmFj+ldKF6vNz7Lxe12T/maYRiqr6+fteUDAABkl4TDh8fjUXNzs6TYOI22trb416a/b21tVXNzs9ra2nTo0CHW+AAAAHE5pmmaVhcxWTAYVHl5uQKBAF0wAACkibl8frO3CwAASCrCBwAASCrCBwAASCrCBwAASCrCBwAASCrCBwAASCrCBwAASCrCBwAASCrCBwAASCrCBwAASCrCBwAASCrCBwAASCqb1QUAAID0FYmaGgqN6Ub/SMJ/hvABAEAGM01ToUhUQ6MRDYbGNByKaDAU0dDoWOw5NKahUETDoYiGw7HnoVBEw+GxSa9n+nrsdSgSlSRFR4cSronwAQBAihiLRDUUjsSDwtTAcPP90HhoGBwdfx4PE0PjwWDi3IlgMRY1l7z2nJzEzyV8AAAwT6ZpanQsqv6RMQ2Ojmlg4jEypsHQ2JTjt5wzGns/OUCExqJLWq/dlqtie56K7DYVO2LPRfY8FdnzVGi3qTA/V0V2mwry8+LHJ14X5uepcPy5yG6LvbbnqWj8+OjQgJxfSqwOwgcAIOtEo6YGxsNBcDis/pEx9Y+Ep4SCgZExDYxGNDAa1uBoRP1Tjt98RJagVcGWm6Mie56KHbb4c2H+pPd2m4ocEwHBFgsUk7828Wcmn5ufJ1ve0s0zCc2h6YPwAQBIOyPhSCw4jISnBIjY+7CCw7EwERy5+Tz5nIHRMZmLmBlycqRiu00ljliLQklBvkoceSpx2FTiGH9dYFOxw6ZSR+y5xDFxfuzPFE6ECLtNdltmT0YlfAAAks40TQ2GIvIPhRQYDiswFJZ/OKzAcFj+odhzYDj2tenhIjgytmjdE3ZbrsoKbCoryFdJwc1AUOKwxcNCybSgUFow9XWxw6ai/Dzl5s5h0EOWI3wAAOYtHIlOCgyh+OubAWIiUITi4SIw/rWFDoLMyZFKHLHgUFpgU1lhfjxITLwvjb/PV1mhLfZcEHsuLYiNbUDyET4AAJJiMy38w2H1DYbkGwypbygk32B4/DkUOz7+euL9YCiyoP+mPS9XzqJ8lRfmx5/LC+2TXo+HBkf+zTAx/lxit9HakKYIHwCQgUzT1MDomHoHQuodHI2FiPHwMDVchNQ3FJZvMNZqMV9lBTY5i+zTQsTN185Cu8onhQxnYezcgvxc5cxljiYyAuEDANJEJGrKPxRSz0BIvQOjujEwGn/dMzCq3oGQesaP9QyManSe4yKcRflyFdlVUWxXRZFdruJ8VRTb48dufi1frmK7SgvylUcLBOaA8AEAFjJNU4HhsLqDo7reP6Lu4Khu9E+EiZtBomcgJN/gqOY6TKLInqfKErtcxQ65iqaFiHi4GA8Y4y0XSzkdE5AIHwCwJEzTVN9QOB4orgdHdL3/5nP3xPv+0TnP3KgoyldliUPLSuyqLHGoqsShymK7lpU6tKzEocoSe+xYiV1Fdm7zSD38VALAHI2ORdQdGNWVwLCuBUbiz/FAMd56MbHnRSKcRflaUVqg5WUOVZU64uEhFiZiQWNZiUOuYrvyaZlAmiN8AMAkE8HiamBYVwMj44/h+PO1wIh6BkIJ/30VRflaUVagqlKHVpQVaPmk5+Xjz1WlDqZ8IqsQPgBklcBwWJf7hnWpb0iX+obHH0PxoNEzMJrQ3+Ow5WpVeYFWlRdqVXmBVo4/pocKh41QAUxH+ACQUQLD4VuCxeTX/SNjd/w77PFgcTNcrHIWalVZgVY5Y8cqivKZIgrME+EDQFoZi0R1NTCi872DutA7pAvjz11zCBeVxXatrSjUmopCra0o0hpnoVY7C+OBw1VsJ1gAS4jwASDljI5F1OUb1kXfoM73DOmibygeNi71DSkcuf1804lwsbaiaPz55us1FYXMAAEsxm8gAEtEo6auBIbVeWNQxo0Bdd4YkHEjFjCuBIZvu+OoPS9X6yqLtN5VpPWVxVpfWaR1LsIFkC74DQWwpIZCYzJuDKrzxsCkoDGocz0DGgnPPhW1yJ6n9ZXF2lBZpHWVRdowHjLWVxZrZVkBK2oCaYzwAWBRBEfCOtPdr9PXBvRmd7/OXh+QcWNAVwIjs/6Z/Lwcra8sVk1VsdxVJXIvK9bGZcVaX1msZSWMuwAyFeEDwJyMhCM6e31Ap6/1683ufp3u7teb1/pvGzJcxXa5lxWrpqpENcuL5V5WoprlJaquKGQpbyALET4AzMg0TXX5hnXySkBvXA3GQkb3gC70Ds66v8jKsgJtXlmqzctLtHlFqdzjLRquYntyiweQ0ggfABSORHX2+oBOXgnq5JWAXr8S1OtXg7NOW60oytfmFaXasrJ0ynN5YX6SKweQjggfQJYZCUfiIePk5aBOXg3ozWsDM+5DYs/L1eaVJdq6qkxbVpZpy4pSbV5ZoqoSB+MxAMwb4QPIYJGoqbPXB/RKl18vX/LrlS6/Tl/r19gM/SalDpvuWV2me1eX6d7V5bp3dZlqqkpktzEmA8DiInwAGcI0TV32D+uVroBeueTXy11+/fRyQEOhyC3nLitxaNuamyHj3tXlqnYV0poBICkIH0CaCkeiOnklqOPnfeq40KfjF/p0o//WTdGK7XnatrZc26udeutap7ZXO7WqvICgAcAyhA8gTQSGw/Je7FPH+T4dv+DTy13+WxbpsuXm6J5VZXrLRNiodqqmqoQFuQCkFMIHkKJ8gyEdM3r1k85etZ/36XR3/y1LjpcX5qt+fYXqNlSofr1Lb1lbroJ8tnAHkNoIH0CKCAyH9dI5n37S2aMXOnt16lr/LeesryxS3foKNWxwqX59hWqqSpRLqwaANEP4ACwyHIro2LlevdAZa904eSVwy+Jdm1eU6G3uSj3orlT9hgotLy2wplgAWESEDyBJTNPUmesD+tHpG/rxmRs6ds6n0NjUMRvuZcXaUVOpt7krtcNdqapSh0XVAsDSIXwASygwHNbzZ3vigePqtP1P1jgL9Y5Ny/S2mljYWFlOywaAzEf4ABZZ540BeV7vlueNbnkv+hWZ1JfisOVqh7tS795cpXdtrlJNVTFTXgFkHcIHsECRqKkTF/t09I1uHX29W8aNwSlf37S8RO+6q0rv3lKlBze6mI0CIOsRPoB5GAlH9OM3b+jo6916+tR19Q6G4l/Lz8vRDneldm5doffevVxrK4osrBQAUg/hA0jQSDiiZ8/06LuvXpHnjesaGL2542tZgU3vuXu5Gu9ZoXdvqVJZAbu7AsBsCB/AbYyORfTsmz367mtXdfT17imBY3V5gR65d6Ue2bpCDRtdys9jAzYASAThA5gmEjX1otGrb3kv6/snr6l/UuBYWVagD2xbpQ++ZZXur3aywBcAzAPhAxh39nq/vum9rH85cXnKlNjlpQ59YNsqPfqWVapdV0HgAIAFInwgq/kGQ/rOK1f0Te8lvXopED9eVmDTo9tX6+e3r1bDBheBAwAWEeEDWScaNfV8Z4++fuyijr7erbHxdTjycnP0ni1V+kjtWr337uVMiQWAJUL4QNboGRhV6/FL+kb7RV3oHYof37amXB+pXaMPbV+tZSUsZw4AS43wgYxmmqZe6OzVP710Ud8/eU3hSKyVo9Rh02O1a/QLD6zTPavKLK4SALIL4QMZaTgU0bdPXNZXnz+ns9cH4se3Vzv1iw+s06PbV6nIzo8/AFiBuy8yyrXAiP7hxfP6+rGL6hsKS5KK7Xn68P1r9PEH1+ne1eUWVwgAIHwgI7x6ya+vPndO//rq1fgA0rUVhfrk2zfoow3VrDgKACmE8IG0ZZqmXjB69dc/7NRzZ3vixx/Y6NKnHtqonVtXKI8psgCQcggfSDumaerpU9f1lz88qxMX/ZIkW26OPrR9tT710EZtW0vXCgCkMsIH0kYkaurfXruqv/rhWZ261i9Jstty9bGGajW9y83usQCQJggfSHmmaep7J7v150dP683u2MyVYnuefult6/Xpd2zU8tICiysEAMzFnMKHYRhqa2uT2+2WYRhqamqS0+mc9dxDhw6ppqZGnZ2d2r9//6znAjMxTVM/PtOjP/v+6fjS52UFNn36HW598u0bVF7EIFIASEdzCh+7d+9WR0eHpFi42LNnj1pbW2c8d+fOnero6JDT6ZTX61Vzc7MOHTq08IqRFdrP+/Qn3zutl875JElF9jx9+h0b9Z/f6VZ5IaEDANJZwuHDMIwp791utzwez4znThyfaOmora3V4cOHCR+4o/M9g2r59zf0vZPdkmJjOj6xY71+7eEalj4HgAyRcPjweDxyuVxTjrlcLnm9XtXW1k457vf7Z/w7Zjp3dHRUo6Oj8ffBYDDRkpBBAkNh/cXTZ/S1F84rHDGVmyM93lCtz7zvLq0qL7S6PADAIko4fMwWKHw+3y3Hamtrp7SUeL3eWc9taWnRE088kWgZyDDhSFRfP3ZR/9PzZnxF0ndvrtLnP3iPNq8otbg6AMBSyF3oXzBTKHG73Tpw4IAOHz4sv98fDyLTW04kaf/+/QoEAvFHV1fXQktCmjh+3qdHv/KcvvD/TqpvKKy7lpfo//xqg772qQcIHgCQwRJu+XA6nbe0XPh8vllnsOzbt0+GYcgwDDU2NkqKhZLpHA6HHA768rOJbzCkln97Q60dlyRJFUX5+twjW/SxhmrZ8hachwEAKS7h8NHY2DjjgNH6+voZzzcMIx42JsZ6MNU2u0Wjpo4c79If/8cp+ce7WD7WUK3mn71bFcV2i6sDACRLwuFjequFYRiqr6+PBwqv1yun0xk/r66uTufOnZPT6dShQ4d04MCBxasaaafzxoCa217V8Qt9kqS7V5bqDx+7T3Xrb+2KAwBkthzTNM1ET55YOKyhoUHt7e1TFg7bvXu3GhoatG/fPknS4cOH5XK55PP55Ha7410vdxIMBlVeXq5AIKCysrK5f0dIKZGoqa8+d05/+v3TGh2Lqtiep/+2c7M++fYNdLEAQAaZy+f3nMJHMhA+MkfnjQH9Tusr8o5v/vbOu5ap5SPb2IMFADLQXD6/2dsFiy4aNfXV58/pT74Xa+0ocdj0ex+8R483VCsnhy3uASDbET6wqK4HR/S51lf07JkeSdK7Nlep5SPbtMbJQmEAgBjCBxbN06e69dutr8o3GFJBfq7+x6P36hceoLUDADAV4QMLNhKO6I///ZT+z0/OS5K2rirTV37hfm1aXmJtYQCAlET4wIJc6B3Ur/2jV69fje3J86mHNqr5/VvksOVZXBkAIFURPjBvT5/q1n/9xssKjoypstiuP/3odr1ny3KrywIApDjCB+YsGjX15R+c0Zd/cEaSVLvOqb/+xTqtLC+wuDIAQDogfGBOAkNhffapE3rm9A1J0i+/bb1+74NbZbexYBgAIDGEDyTsfM+gPvW1dhk3BlWQn6s/emybPlK71uqyAABphvCBhBwzerX3HzvkHwprdXmBnvyVet27utzqsgAAaYjwgTv6Zscl/e63XlU4Ymr72nI9+Sv1Wl7K+A4AwPwQPjAr0zT150ff1F88fVaS9IFtK/Vnu9+qQjvTaAEA80f4wIzGIlF9/ts/1VPHuyRJv/GeGn1u5xbl5rJaKQBgYQgfuMVIOKLPfuOEvneyW7k50h89tk0fe2Cd1WUBADIE4QNT9I+E1fT3HXrB6JXdlquvfOx+/ex9K60uCwCQQQgfiOsbDOkTXz2mn14OqsRh0+FfrtPba5ZZXRYAIMMQPiBJ8g2G9It/e0xvXA2qstiur33qAd23hqm0AIDFR/iAfIMhffzJF3XqWr+WlTj0jaYHtWl5qdVlAQAyFOEjy/UOjOoX//aYTl3rV1WpQ/+8Z4c2LS+xuiwAQAYjfGSxWIvHMZ3uJngAAJKH8JGl+kfC+uT/fkmnu/u1vNShf27aoZoqggcAYOmxFWkWGglHtOfvj+vVSwG5iu36+p4HCR4AgKQhfGSZsUhUv/n1E3rR8KnEYdPXfvUBBpcCAJKK8JFFTNPUvm++Ks8b3bLbcvW3v1KvbWuZTgsASC7CRxb50tE39S3vZeXl5uivP16rHe5Kq0sCAGQhwkeWaD3epa+M7077R4/dp8atKyyuCACQrQgfWeAnZ3u0/1uvSZJ+/eEaPd7AJnEAAOsQPjLc2ev92vuPHRqLmnr0Lav0249ssbokAECWI3xksMBQWJ/+2nH1j4ypbn2F/nT3duXm5lhdFgAgyxE+MlQkauqzT53Qhd4hrXEW6vAn6lSQn2d1WQAAED4y1ZeOvqlnTt+Qw5arQ5+oU2WJw+qSAACQRPjISP/x06v6yx/GZrb88X/apvvWsJYHACB1ED4yzNnr/frckVckSZ96aKMeu3+txRUBADAV4SODjIQj+o1/OqHBUEQ73C799w/cbXVJAADcgvCRQb74r6/rdHe/lpU49Be/UCtbHv97AQCph0+nDPFvr13VPx27KEn60uPbVVXKAFMAQGoifGSALt+Qmr/5qiTp1x6u0TvvqrK4IgAAZkf4SHPhSFSf+cYJ9Y+MqXadU7+1c7PVJQEAcFuEjzT3v57p1ImLfpUW2PTlj92vfMZ5AABSHJ9UaezklYC+8oMzkqQ/+PB9qnYVWVwRAAB3RvhIU6NjEX3uyCsai5p6/30r9XPbV1tdEgAACSF8pKkve87o1LV+VRbb9Qcfvk85OWwYBwBID4SPNHTiYp/+5kedkqQ/fGwb+7YAANIK4SPNhCNR/e43X1PUlD781tX62ftWWl0SAABzQvhIM08+a+h0d79cxXZ94UP3Wl0OAABzRvhIIxd6B/VlT2x2y+c/cI8qiu0WVwQAwNwRPtKEaZr6/Ld/qtGxqB7aVKmP1K6xuiQAAOaF8JEm/uXly3rubI/stlz9wYe3MbsFAJC2CB9poH8krD/87ilJ0mfeu0kblxVbXBEAAPNH+EgDf/n0WfUMjGrjsmI1vavG6nIAAFgQwkeKM24M6KvPn5Mk/f6j98hu438ZACC98UmW4v7gu28oHDH18JYqvffuFVaXAwDAghE+UtgPT1/X06euy5abo99/dKvV5QAAsCgIHykqHInqi//6uiTpk2/foJqqEosrAgBgcRA+UtQ32rtk3BhUZbFdn2m8y+pyAABYNISPFDQUGtNXfhBbyfS/vHeTygryLa4IAIDFQ/hIQV997pxu9I+q2lWojz+43upyAABYVISPFNM3GNKhHxmSpN9+ZAtTawEAGYdPthTzVz88q/7RMW1dVaYPvWW11eUAALDoCB8p5LJ/WH//wgVJ0r6f3aLcXPZvAQBkHsJHCvmrH55VKBLVDrdL795cZXU5AAAsCcJHirjsH1br8S5J0m/t3MKutQCAjEX4SBF/80ynwhFTO9wuPbDRZXU5AAAsGcJHCrgWGNFT7bFWj8++b7PF1QAAsLRscznZMAy1tbXJ7XbLMAw1NTXJ6XTOeq7H45HL5ZJhGNq1a5fcbvdi1Jxx/uZHnQpFonpgg0s73LR6AAAy25zCx+7du9XR0SEpFi727Nmj1tbWGc9ta2vTvn374u/37t2rQ4cOLaDUzHQ9OKJ/fumiJOkz77uLsR4AgIyXcLeLYRhT3rvdbnk8nlnPf+qpp+ZfVRY5/GNDo2NR1a2v0EObKq0uBwCAJZdw+JjoQpnM5XLJ6/XOeL7L5VJdXV28+2Xnzp0LqzQDBYbD8VaP33zvJlo9AABZIeHw4ff7Zzzu8/lmPD7RHVNTU6PW1lbt2rVrxvNGR0cVDAanPLLFP790UYOhiLasKNXDrOsBAMgScxrzMZPZQonH49GBAwdkGIb27t0rSTOO+WhpadETTzyx0DLSTmgsqv/9/DlJ0p53uWn1AABkjYRbPpxO5y2tHD6fb8bZLoZhqL29XY2NjWpqalJnZ6eOHDlyy7gRSdq/f78CgUD80dXVNffvIg1955Ur6g6OakWZQz+3nT1cAADZI+Hw0djYOOPx+vr6W455vV41NDTE37vdbu3fv3/GVhKHw6GysrIpj0xnmqaefDYWxD759o3sXAsAyCoJf+pNX6PDMAzV19fHWz68Xm+8ZaO2tlbt7e1Tzu/t7VVtbe0Cy80Mz57p0alr/Sq25+njD66zuhwAAJJqTmM+Wltb1dzcrIaGBrW3t09Z46OlpUUNDQ3at2+f3G63du7cqYMHD8bDycS4Dyje6vF4wzqVF+ZbXA0AAMmVY5qmaXURkwWDQZWXlysQCGRkF8yb3f165Es/Vm6O9KPfeY+qXUVWlwQAwILN5fObwQZJ9g8vXJAk7dy6guABAMhKhI8k6h8J61veS5KkX37bBmuLAQDAIoSPJPr2icsaDEVUU1Wst9ewlDoAIDsRPpLENE39/XiXyyd2rGdRMQBA1iJ8JMkLRq/OXh9QkT1PH6lba3U5AABYhvCRJBMDTR+7f43KCpheCwDIXoSPJLjeP6Kjr3dLkj7xtvUWVwMAgLUIH0nwLe9ljUVN3b/OqbtXZt7aJQAAzAXhY4mZpqkj7bHN8h6vr7a4GgAArEf4WGLHL/TJ6BlUkT1Pj7J7LQAAhI+l9tR4q8cHt61SiWNOW+kAAJCRCB9LqH8krO++elWS9HgDXS4AAEiEjyX1nVeuajgckbuqWHXrK6wuBwCAlED4WEJHjse6XD7WUM2KpgAAjCN8LJFzPYN6ucuv3Bzpw/evsbocAABSBuFjifzLicuSpHfcVaXlpQUWVwMAQOogfCwB0zT1Ly/Hwsdj9zO9FgCAyQgfS+BEl18XeodUmJ+nR7autLocAABSCuFjCUx0ufzMvStUzNoeAABMQfhYZOFIVP86vrbHY7VrLa4GAIDUQ/hYZM+euSHfYEjLShx6qKbS6nIAAEg5hI9F9u0TVyRJP7d9tWx5XF4AAKbj03ERjYQj+sEb3ZKkD21fZXE1AACkJsLHInrm9HUNhSJa4yzUW6udVpcDAEBKInwsou++dk2S9IFtK1lOHQCAWRA+FsnkLpcPbKPLBQCA2RA+FgldLgAAJIbwsUjocgEAIDGEj0VAlwsAAIkjfCyCZ07foMsFAIAEET4WwX/8NLac+vvvo8sFAIA7IXwsUDgS1dOnrkuSfuY+drAFAOBOCB8L1H7Op+DImCqL7apdV2F1OQAApDzCxwJ9//XYQNP33r1cebl0uQAAcCeEjwUwTVNHx8PHzq0rLK4GAID0QPhYgDeu9uuyf1gF+bl6511VVpcDAEBaIHwsgGd8bY93bKpSoT3P4moAAEgPhI8FmOhyeYQuFwAAEkb4mKergWG9djmgnBzpPXcvt7ocAADSBuFjnjzjrR616ypUVeqwuBoAANIH4WOevs8sFwAA5oXwMQ+Do2M6ZvgkSY33ED4AAJgLwsc8vNDZq1AkqmpXoWqqiq0uBwCAtEL4mIdn3ozt5fLw5uVsJAcAwBwRPubINE09c/qGJOnhLSwsBgDAXBE+5sjoGdSlvmHZ83L1tppKq8sBACDtED7maKLV44GNLhXZbRZXAwBA+iF8zNEzp8fHe9DlAgDAvBA+5mA4FNGxc7EptoQPAADmh/AxBy8YPQqNRbXGWaiaqhKrywEAIC0RPuZg8iwXptgCADA/hI8ETZ1iy0ZyAADMF+EjQed7h3TRN6T8vBym2AIAsACEjwQ9d7ZHklS3vkIlDqbYAgAwX4SPBD1/JhY+3rFpmcWVAACQ3ggfCYhETf2kMxY+HiJ8AACwIISPBJy8ElBwZEylBTZtW1NudTkAAKQ1wkcCJsZ77HBXypbHJQMAYCH4JE3A82cZ7wEAwGIhfNzBSDii9vN9khjvAQDAYiB83EHHhT6FxqJaUeZQTVWx1eUAAJD2CB93MDHe46FNy1hSHQCARUD4uAPGewAAsLgIH7cRGArrtcsBSYz3AABgsRA+buMFo0emKW1aXqIVZQVWlwMAQEaY0yYlhmGora1NbrdbhmGoqalJTqdzxnPb2trU2NgoSbOek+p+0tkrSXqIjeQAAFg0cwofu3fvVkdHh6RYENmzZ49aW1tnPXe6AwcOaN++ffMo0xovGrHwscNN+AAAYLEk3O1iGMaU9263Wx6PZ8Zz/X6/WltbZZpm/JFuwcM3GNKb3QOSpAc2uiyuBgCAzJFw+PB4PHK5pn4Iu1wueb3eGc/ftWtX/HVbW9uU95ONjo4qGAxOeaSCl87FWj02ryhRZYnD4moAAMgcCYcPv98/43Gfz3fLscljPPx+v3w+n9xu94x/vqWlReXl5fFHdXV1oiUtqReN2PdFlwsAAItrwbNdZgslE5qbm/XRj3501q/v379fgUAg/ujq6lpoSYtiYrzHgxsJHwAALKaEB5w6nc5bWjl8Pt9tZ7L4/X55PJ7bnuNwOORwpFa3Rt9gSKeu9UuSHnQz3gMAgMWUcMvHxLTZ6err62f9M8ePH0/LabYvnY+FrE3LS7SM8R4AACyqhMPH9DEbhmGovr4+Hi68Xu8tM2K8Xu8tg1TTwbH4eI/0qx0AgFQ3pzEfra2tam5uVltbmw4dOjRljY+Wlha1tbXd8mdmG2iayhjvAQDA0skxTdO0uojJgsGgysvLFQgEVFZWlvT/fmAorLd+8fsyTemlz79Py0tZVh0AgDuZy+c3e7tM89J5n0xTclcVEzwAAFgChI9pjtHlAgDAkiJ8THPsHINNAQBYSoSPSYIjYZ28EpDEyqYAACwVwsckHRf6FDWlDZVFWlHGeA8AAJYC4WOS4+OLizVsoMsFAIClQviY5Pj5PklS/YYKiysBACBzET7GhcaiernLL0mqW0/LBwAAS4XwMe7klYBGx6KqKMpXTVWx1eUAAJCxCB/jOi7Eulzq1ruUk5NjcTUAAGQuwse49vHBpoz3AABgaRE+JJmmGW/5aCB8AACwpAgfki70DqlnICS7LVf3rSm3uhwAADIa4UM3u1zesqZcDluexdUAAJDZCB+6Odi0nsXFAABYcoQPTRpsup7xHgAALLWsDx99gyF13hiUJNURPgAAWHJZHz4mulw2LS9RRbHd4moAAMh8WR8+2i/Q5QIAQDJlffjoOM9gUwAAkimrw8dIOKJXLwUk0fIBAECyZHX4+OnlgEKRqJaV2LW+ssjqcgAAyApZHT6OT6zvwWZyAAAkTXaHj/h4D7pcAABIlqwNH7HN5GIzXVjfAwCA5Mna8NF5Y1B9Q2EV5Ofq3tVsJgcAQLJkbfg4Pr6k+va1TtltWXsZAABIuqz91I0PNmW8BwAASZW14YOdbAEAsEZWho8b/aM61zOonBypdh0tHwAAJFNWho+JVo/Ny0tVXphvcTUAAGSXLA0f45vJMd4DAICky8rw0c7iYgAAWCbrwsdwKKKTVyY2k2OwKQAAyZZ14eOVS36FI6ZWlDm0tqLQ6nIAAMg6WRc+OthMDgAAS2Vd+JhY2ZTxHgAAWCOrwkc0ak5p+QAAAMmXVeHjzPUBBUfGVGTP0z2rSq0uBwCArJRV4eP4+Poe969zypaXVd86AAApI6s+gY+Pr+9RR5cLAACWyarw0T4+2LSBwaYAAFgma8LH1cCwLvUNKzdHup/N5AAAsEzWhI+JLpetq8tU4rBZXA0AANkra8IHU2wBAEgNWRM+bo73IHwAAGClrAgf/SNhvXE1KImVTQEAsFpWhI8TF/2KmlK1q1ArygqsLgcAgKyWFeFjYj+XBsZ7AABguawIH+3jM13qGe8BAIDlMj58hCNRneiKhQ8WFwMAwHoZHz5OXglqJByVsyhfNVUlVpcDAEDWy/jwMTHeo25dhXJzcyyuBgAAZEH4YLwHAACpJKPDh2maOn6BzeQAAEglGR0+zvcOqWcgJLstV9vWlltdDgAAUIaHj4kl1bevLZfDlmdxNQAAQMr08HFufLApi4sBAJAyMjp8HBsPHw+6CR8AAKSKjA0fl/3DuugbUl5ujurXM9gUAIBUkbHh45jRK0m6b025SgvyLa4GAABMyNjw8eJ4+NixkS4XAABSScaGj4nxHjvclRZXAgAAJsvI8HHFP6wLvUPKzZHqWVwMAICUYpvLyYZhqK2tTW63W4ZhqKmpSU6nc9bzPR6PDMOQ2+2WJDU2Ni6o2EQdOxfrctnGeA8AAFLOnMLH7t271dHRISkWRPbs2aPW1tYZz/V4PGptbdWhQ4dkGIZ27typzs7OhVecgBc76XIBACBVJRw+DMOY8t7tdsvj8cx6/t69e+NBxe126+jRo/Msce5eHG/5IHwAAJB6Eh7z4fF45HJNnTnicrnk9XpvOdcwDPl8PjmdTnm9Xvn9/njXy1K7GmC8BwAAqSzh8OH3+2c87vP5bjnm9Xrlcrni40MOHz6stra2Gf/86OiogsHglMdCHDNi9bC+BwAAqWlOYz5mMlMo8fl8MgxDjY2NcjqdampqUkVFhUzTvOXclpYWPfHEEwstI+75sz2S6HIBACBVJdzy4XQ6b2nlmOhamc7tdsvpdMa/NvE8UxfN/v37FQgE4o+urq7Eq5/GNE09eyYWPt5517J5/z0AAGDpJBw+ZpsmW19ff8uxuYzvcDgcKisrm/KYr84bA7oWHJHDlquGDaxsCgBAKko4fEwPFIZhqL6+fkqrxsSMGLfbrfr6+niXzMRaH7W1tYtT9Sx+/Gas1eOBjS4V5Oct6X8LAADMz5zGfLS2tqq5uVkNDQ1qb2+fssZHS0uLGhoatG/fvinn1tXVqaOjIylTbZ89c0MSXS4AAKSyHHOmUaAWCgaDKi8vVyAQmFMXzOhYRG994qiGwxH9+2ffqXtWzb/7BgAAzM1cPr8zZm+Xjgt9Gg5HtKzEobtXllpdDgAAmEXGhI/nxme5vOuuZcrJybG4GgAAMJuMCR/xKbabGe8BAEAqy4jwcb1/RD+9EpAkPbSJ8AEAQCrLiPDxw1PXZZrS9rXlWl5aYHU5AADgNjIifBx9/bokqfGeFRZXAgAA7iTtw8dwKKLnzsbW93gf4QMAgJSX9uHj+bM9GglHtcZZqHtWMcUWAIBUl/bhw/NGtySp8Z7lTLEFACANpHX4iEZNed4YH++xlS4XAADSQVqHj+MX+tQzMKpSh00Pbqy0uhwAAJCAtA4f33nliiTpZ+5bKbstrb8VAACyRtp+Yo9Fovq3165Kkj60fbXF1QAAgESlbfh4wehV72BIFUX5ensNXS4AAKSLtA0f/+/lWJfL+7etUn5e2n4bAABknbT81B4YHdN3x7tcfp4uFwAA0kpaho/vvHJFQ6GI3MuK9cBGl9XlAACAOUjL8PGN9i5J0uMN1SwsBgBAmkm78PH6laBe6fIrPy9H/6lurdXlAACAOUq78PHV589Jkh7ZulLLShwWVwMAAOYqrcLHtcCI/u/LlyVJn37nRourAQAA85FW4ePvnjMUjphq2FCh2nUVVpcDAADmIW3CxxX/sL72wgVJ0q8/vMniagAAwHylTfj486NvKjQW1YMbXXp4S5XV5QAAgHlKi/Dxk7M9auu4JElqfv/dTK8FACCNpXz48A+F1PytVyVJv7RjHWM9AABIcykdPkbCETX9Q4e6fMNaW1Go333/PVaXBAAAFshmdQGz6bzery9+/zV5L/pV6rDp736lQSWOlC0XAAAkKGU/zX/+r36iXEeRSgtsevKX67VlZanVJQEAgEWQsuEjN0d6ZOsK/d4Ht2pdZZHV5QAAgEWSsuHjx/veo3Url1ldBgAAWGQpO+DUWWS3ugQAALAEUjZ8AACAzET4AAAASUX4AAAASUX4AAAASUX4AAAASUX4AAAASUX4AAAASUX4AAAASUX4AAAASUX4AAAASUX4AAAASUX4AAAASUX4AAAASWWzuoDpTNOUJAWDQYsrAQAAiZr43J74HL+dlAsfvb29kqTq6mqLKwEAAHPV39+v8vLy256TcuHD5XJJki5evHjH4jNdMBhUdXW1urq6VFZWZnU5luE6xHAdYrgOMVyHGK5DTCpcB9M01d/fr9WrV9/x3JQLH7m5sWEo5eXlWf2DNFlZWRnXQlyHCVyHGK5DDNchhusQY/V1SLTRgAGnAAAgqQgfAAAgqVIufDgcDn3hC1+Qw+GwuhTLcS1iuA4xXIcYrkMM1yGG6xCTbtchx0xkTgwAAMAiSbmWDwAAkNkIHwAAIKkIHwAAIKlSKnwYhqGDBw+qra1NBw8elN/vt7okS3i9Xh08eFAHDx7U7t27s/Y6TNbc3JzV18Hj8ejw4cPyeDzyeDxWl2MJwzB0+PDh+P3BMAyrS0oar9erurq6W45n2z1ztuuQbffM2a7DZCl/zzRTSG1tbfx1Z2enuWvXLgursc6BAwemvJ58XbJRR0eHKcns6+uzuhRLHD161GxqajJNM/Z74Xa7La7IGpN/L0zTjF+TTNfa2hr/HZgum+6Zt7sO2XTPvN11mJAO98yUCR+dnZ23/MA4nU6LqrFOR0fHlO+7s7PTlGR2dnZaWJW1WltbTbfbndK/SEtp+veerT8L0+8P2RI+Jkz/sMnWe+b065Ct98zbhY90uGemTLeLx+OJ7+syweVyyev1WlSRNWpra/Xkk0/G3080m02/Ntmira1Nu3btsroMyxiGIZ/PJ6fTKa/XK7/fL7fbbXVZlnC5XKqrq5NhGPJ4PNq5c6fVJVmKe2YM98yp0uWemTLhY7a+KZ/Pl9xCUsDkH5ynnnpKjY2Ncjqd1hVkEb/fn5Xf92Rer1cul0ttbW1yu93xMQ/ZqLW1VZJUU1Oj1tbWtLjBLiXumTdxz4xJp3tmym0sN11KD5hZYn6/X21tbero6LC6FEscOXJETU1NVpdhKZ/PJ8Mw4jfTpqYmVVRUyMzCtQE9Ho8OHDggwzC0d+9eSdKhQ4csrir1cM/knpkOUqblw+l03pLYJ5qbs1Vzc7OOHj2aldfA4/Hoox/9qNVlWM7tdsvpdMZ/Biaes61p3TAMtbe3q7GxUU1NTers7NSRI0eyasbLdNwzb8U9M33umSnT8tHY2Djjv2Lq6+stqMZ6Bw8eVHNzs9xud/xfMtn2C3XkyJH4a8Mw1NLSoscff1y1tbUWVpVc2Tq+Yzqv16uGhob4e7fbrf3792f1v/K5Z07FPTPN7plWj3idbPq0scbGRgursU5ra6t59OhR0zRNs6+vzzx06JDFFVlPWTB6fTaNjY3xUevZOtW2s7PT3Ldv35Rj099nOs0wdTIb75kzXYdsvGfOdB2mfz2V75kp0/IhxQaUNTc3q6GhQe3t7fEBZtnEMAzt3r17yrGJvv5s5Pf7dfjwYUnSgQMHtHfv3tRM8Uto4veirq5OHR0dOnr0qNUlJZ3b7dbOnTt18ODB+L9mJ8Z9ZDqPxxP/f97S0qKGhob4AMtsumfOdh2y7Z55u58HKX3umexqCwAAkiplBpwCAIDsQPgAAABJRfgAAABJRfgAAABJRfgAAABJRfgAAABJRfgAAABJRfgAAABJRfgAAABJRfgAAABJRfgAAABJ9f8Bm0k0iJQYiYIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(xx_double,sqrt(metric_term3_array))\n",
    "xlim(0,15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4baa7ff5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.04499027, 2.04492344, 2.04472301, ..., 1.00336196, 1.00336179,\n",
       "       1.00336161])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "metric_term2_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "0c6fe81d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(yy,t):\n",
    "    Rp,Ip,dRp,dIp,Rm,Im,dRm,dIm=yy\n",
    "    \n",
    "    RKelement=round((t-xx_double[0])/(xx_double[-1]-xx_double[0]) * (len(xx_double)-1))\n",
    "    phi=phi_array[RKelement]\n",
    "    U=1.0+2.*g*phi**2+3.*h*phi**4\n",
    "    W=g*phi**2+2*h*phi**4\n",
    "    metric_term=metric_term_array[RKelement]\n",
    "    metric_term2=metric_term2_array[RKelement]\n",
    "    metric_term3=metric_term3_array[RKelement]\n",
    "    \n",
    "    dRpdt=dRp\n",
    "    dIpdt=dIp\n",
    "    ddRpdt = -2.0*dRp/(t+1e-10) - (metric_term)*Rp - metric_term2*((omegap**2/metric_term3-U)*Rp-W*Rm)       \n",
    "    ddIpdt = -2.0*dIp/(t+1e-10) - (metric_term)*Ip - metric_term2*((omegap**2/metric_term3-U)*Ip+W*Im)       \n",
    "    dRmdt=dRm\n",
    "    dImdt=dIm\n",
    "    ddRmdt = -2.0*dRm/(t+1e-10) - (metric_term)*Rm - metric_term2*((omegam**2/metric_term3-U)*Rm-W*Rp)       \n",
    "    ddImdt = -2.0*dIm/(t+1e-10) - (metric_term)*Im - metric_term2*((omegam**2/metric_term3-U)*Im+W*Ip)       \n",
    "    \n",
    "    return array([dRpdt,dIpdt,ddRpdt,ddIpdt,dRmdt,dImdt,ddRmdt,ddImdt])\n",
    "\n",
    "def RK(f,t,U0):\n",
    "    N=len(f(U0,t[0]))\n",
    "    U=zeros((N,len(t)))\n",
    "    k1,k2,k3,k4=zeros((4,N,len(t)))\n",
    "    U[:,0]=U0\n",
    "    for ii in range(len(t)-1):\n",
    "        Uii=U[:,ii]\n",
    "        h=t[ii+1]-t[ii]\n",
    "        k1[:,ii]=h*f(Uii,t[ii])\n",
    "        k2[:,ii]=h*f(Uii+0.5*k1[:,ii],t[ii]+0.5*h)\n",
    "        k3[:,ii]=h*f(Uii+0.5*k2[:,ii],t[ii]+0.5*h)\n",
    "        k4[:,ii]=h*f(Uii+k3[:,ii],t[ii]+h)\n",
    "        U[:,ii+1]=Uii+k1[:,ii]/6.+k2[:,ii]/3.+k3[:,ii]/3.+k4[:,ii]/6.\n",
    "    return U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "dc07dcce",
   "metadata": {},
   "outputs": [],
   "source": [
    "def integrate_THETA(Rm00,Im00):\n",
    "    Rm0=Rm00\n",
    "    Im0=Im00\n",
    "    U0=Rp0,Ip0,dRp0,dIp0,Rm0,Im0,dRm0,dIm0\n",
    "\n",
    "    SOL=RK(f,xx,U0)\n",
    "\n",
    "    Rp_sol=SOL[0]\n",
    "    Ip_sol=SOL[1]\n",
    "    Rm_sol=SOL[4]\n",
    "    Im_sol=SOL[5]\n",
    "    dRp_sol=SOL[2]\n",
    "    dIp_sol=SOL[3]\n",
    "    dRm_sol=SOL[6]\n",
    "    dIm_sol=SOL[7]\n",
    "\n",
    "    kp=sqrt(omegap**2-1.)\n",
    "    km=sqrt(omegam**2-1.)\n",
    "\n",
    "    chip=Rp_sol+1j*Ip_sol\n",
    "    dchip=dRp_sol+1j*dIp_sol\n",
    "    chim=Rm_sol+1j*Im_sol\n",
    "    dchim=dRm_sol+1j*dIm_sol\n",
    "\n",
    "    Ap1=0.5*xx[-1]*exp(-1j*kp*xx[-1])*(kp*chip[-1]-1j*dchip[-1])\n",
    "    Am1=0.5*xx[-1]*exp(-1j*km*xx[-1])*(km*chim[-1]-1j*dchim[-1])\n",
    "    Bp1=0.5*xx[-1]*exp(1j*kp*xx[-1])*(kp*chip[-1]+1j*dchip[-1])\n",
    "    Bm1=0.5*xx[-1]*exp(1j*km*xx[-1])*(km*chim[-1]+1j*dchim[-1])\n",
    "    \n",
    "    second_point = -100   #this will take also this point and take the average -10 and the trivial -1 also work\n",
    "    Ap2=0.5*xx[second_point]*exp(-1j*kp*xx[second_point])*(kp*chip[second_point]-1j*dchip[second_point])\n",
    "    Am2=0.5*xx[second_point]*exp(-1j*km*xx[second_point])*(km*chim[second_point]-1j*dchim[second_point])\n",
    "    Bp2=0.5*xx[second_point]*exp(1j*kp*xx[second_point])*(kp*chip[second_point]+1j*dchip[second_point])\n",
    "    Bm2=0.5*xx[second_point]*exp(1j*km*xx[second_point])*(km*chim[second_point]+1j*dchim[second_point])\n",
    "    \n",
    "    #Ap2=0.5*xx[-1]*exp(-1j*kp*xx[-1])*(kp*chip[-1]-1j*dchip[-1])\n",
    "    #Am2=0.5*xx[-1]*exp(-1j*km*xx[-1])*(km*chim[-1]-1j*dchim[-1])\n",
    "    #Bp2=0.5*xx[-1]*exp(1j*kp*xx[-1])*(kp*chip[-1]+1j*dchip[-1])\n",
    "    #Bm2=0.5*xx[-1]*exp(1j*km*xx[-1])*(km*chim[-1]+1j*dchim[-1])\n",
    "    \n",
    "    Ap=(Ap1+Ap2)/2.     #here we take the average\n",
    "    Am=(Am1+Am2)/2.\n",
    "    Bp=(Bp1+Bp2)/2.\n",
    "    Bm=(Bm1+Bm2)/2.\n",
    "\n",
    "    #SANITY CHECK  (check that every solution satisfy eq. 44 of BS superradiance manuscript)\n",
    "\n",
    "    check1=abs(Ap)**2/kp+abs(Bm)**2/km\n",
    "    check2=abs(Am)**2/km+abs(Bp)**2/kp\n",
    "\n",
    "    if (check1-check2)/(check1+check2)>1e-3:\n",
    "        print(\"there was a problem \"+str((check1-check2)/(check1+check2)))\n",
    "    return array([Ap,Am,Bp,Bm])\n",
    "\n",
    "\n",
    "def one_mode(xxxx):\n",
    "    Ap,Am,Bp,Bm = integrate_THETA(xxxx[0],xxxx[1])\n",
    "    if omega>0:\n",
    "        return array([Am.real,Am.imag])\n",
    "    if omega<0:\n",
    "        return array([Bm.real,Bm.imag])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "f6af9b2c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0[ 0.0478368  -0.06608159]\n",
      "1[ 0.05057904 -0.07013735]\n",
      "2[ 0.05357091 -0.0751408 ]\n",
      "3[ 0.05713006 -0.07974438]\n",
      "4[ 0.06024529 -0.08509872]\n",
      "5[ 0.06404532 -0.09098932]\n",
      "6[ 0.06809022 -0.09716148]\n",
      "7[ 0.07250435 -0.10420999]\n",
      "8[ 0.07634266 -0.11185212]\n",
      "9[ 0.08188575 -0.12017299]\n",
      "10[ 0.0870658  -0.12922834]\n",
      "11[ 0.09232201 -0.13904743]\n",
      "12[ 0.09867978 -0.15037507]\n",
      "13[ 0.10550573 -0.16315489]\n",
      "14[ 0.11210739 -0.17743673]\n",
      "15[ 0.11986433 -0.19413188]\n",
      "16[ 0.1259468  -0.21429612]\n",
      "17[ 0.13322899 -0.23608931]\n",
      "18[ 0.13774995 -0.26308278]\n",
      "19[ 0.14290943 -0.29252454]\n",
      "20[ 0.14290943 -0.29252454]\n",
      "21[ 0.141711   -0.29556817]\n",
      "22[ 0.14211921 -0.29526783]\n",
      "23[ 0.14295111 -0.29865276]\n",
      "24[ 0.1414312  -0.29919996]\n",
      "25[ 0.14340059 -0.30085914]\n",
      "26[ 0.14185439 -0.30353221]\n",
      "27[ 0.14264777 -0.30344511]\n",
      "28[ 0.14304536 -0.30695275]\n",
      "29[ 0.14173277 -0.30730936]\n",
      "30[ 0.14358805 -0.30929076]\n",
      "31[ 0.1419291  -0.31178787]\n",
      "32[ 0.1428085 -0.3118653]\n",
      "33[ 0.14303969 -0.31538416]\n",
      "34[ 0.14176234 -0.31578238]\n",
      "35[ 0.14351404 -0.31780406]\n",
      "36[ 0.14188238 -0.32036123]\n",
      "37[ 0.14259717 -0.32049177]\n",
      "38[ 0.14291791 -0.32396535]\n",
      "39[ 0.14150148 -0.32459939]\n",
      "40[-0.04497692  0.04890472]\n",
      "41[-0.0475505   0.05131215]\n",
      "42[-0.05014778  0.05380966]\n",
      "43[-0.05287951  0.05627307]\n",
      "44[-0.05578486  0.05875187]\n",
      "45[-0.0588225   0.06149252]\n",
      "46[-0.06239651  0.06403185]\n",
      "47[-0.0657971   0.06718593]\n",
      "48[-0.06957245  0.06998108]\n",
      "49[-0.07325948  0.07251253]\n",
      "50[-0.07750287  0.07553632]\n",
      "51[-0.0819438  0.078462 ]\n",
      "52[-0.08674251  0.08129579]\n",
      "53[-0.09128557  0.08451766]\n",
      "54[-0.09709706  0.08676191]\n",
      "55[-0.10212836  0.08888404]\n",
      "56[-0.10765497  0.09165653]\n",
      "57[-0.11481286  0.09297893]\n",
      "58[-0.12146904  0.09381168]\n",
      "59[-0.1289325  0.0924598]\n",
      "60[-0.1289325  0.0924598]\n",
      "61[-0.12972722  0.09212694]\n",
      "62[-0.13056511  0.09207311]\n",
      "63[-0.13128335  0.09220707]\n",
      "64[-0.13184737  0.09238986]\n",
      "65[-0.13231414  0.09252056]\n",
      "66[-0.13277108  0.09255153]\n",
      "67[-0.13329288  0.09245512]\n",
      "68[-0.1339097   0.09217547]\n",
      "69[-0.13454509  0.0916042 ]\n",
      "70[-0.1349013  0.0906739]\n",
      "71[-0.13446737  0.08972974]\n",
      "72[-0.13327236  0.08998419]\n",
      "73[-0.1336266   0.09224633]\n",
      "74[-0.1375141   0.09197547]\n",
      "75[-0.13588874  0.08785337]\n",
      "76[-0.13573204  0.0933685 ]\n",
      "77[-0.13662288  0.08623851]\n",
      "78[-0.14259078  0.09032635]\n",
      "79[-0.14413062  0.09360818]\n"
     ]
    }
   ],
   "source": [
    "### NOT CHANGE\n",
    "Rp0=1.0    #\n",
    "Ip0=0.0    #\n",
    "dRp0=0.0   #\n",
    "dIp0=0.0   #\n",
    "dRm0=0.0   #\n",
    "dIm0=0.0   #\n",
    "###########\n",
    "\n",
    "OMEGAm=concatenate((geomspace(-5,-omegaB-1.1,20), geomspace(-omegaB-1.1,-omegaB-1.01,20)))\n",
    "OMEGAp=concatenate((geomspace(5,omegaB+1.1,20), geomspace(omegaB+1.1,omegaB+1.01,20)))\n",
    "OMEGA=concatenate((OMEGAm,OMEGAp))\n",
    "AE=zeros(len(OMEGA))\n",
    "AE_Cardoso=zeros(len(OMEGA))\n",
    "\n",
    "semilla = array([0.04623308, -0.01241643])\n",
    "\n",
    "\n",
    "for cc in arange(len(OMEGA)):\n",
    "    omega=OMEGA[cc]\n",
    "    omegap=omegaB+omega\n",
    "    omegam=omegaB-omega\n",
    "\n",
    "    \n",
    "    root_ = optimize.root(one_mode, x0=semilla)\n",
    "    \n",
    "    Ap,Am,Bp,Bm = integrate_THETA(root_.x[0],root_.x[1])\n",
    "    \n",
    "    kp=sqrt(omegap**2-1.) # require to be recalculated here, sorry\n",
    "    km=sqrt(omegam**2-1.)\n",
    "    \n",
    "    temp1=(omegap**2/kp**2)*abs(Ap)**2+(omegam**2/km**2)*abs(Bm)**2\n",
    "    temp2=(omegap**2/kp**2)*abs(Bp)**2+(omegam**2/km**2)*abs(Am)**2\n",
    "\n",
    "    AE[cc]=(temp1/temp2)**(sign(omega))\n",
    "\n",
    "    temp1=(omegap/kp)*abs(Ap)**2-(omegam/km)*abs(Bm)**2\n",
    "    temp2=(omegap/kp)*abs(Bp)**2-(omegam/km)*abs(Am)**2\n",
    "\n",
    "    AE_Cardoso[cc]=abs(temp1/temp2)**(sign(omega))\n",
    "\n",
    "    print(cc , end='')\n",
    "    semilla = root_.x\n",
    "    print(semilla)\n",
    "    #if sign(OMEGA[cc])!=sign(OMEGA[cc+1]):\n",
    "    #    semilla = array([0.04623308, -0.01241643])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "0413fcbd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x1276450f0>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiYAAAGfCAYAAABr4xlmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAABBV0lEQVR4nO3de3xUd53/8feZXElCMgxQEioFBi81sLXmgtuqWyvBokhXd6G6Wi8/dwGru176WEjRn4tZd6Xg/tTuaruUXVdrV1eS30MFW9Mlu1ZdW0vKWCvN+tNlKEUI0JLLkISEkDm/P4YZcpnLmZkzlzN5PR+P6MyZb5Jv8zhM3vlePl/DNE1TAAAAecCV6w4AAACEEUwAAEDeIJgAAIC8QTABAAB5g2ACAADyBsEEAADkDYIJAADIGwQTAACQNwgmAAAgbxQn09jn82nz5s06cuRIwnZdXV2SpO7ubu3bt09ut1uS5Pf71dHRIa/XK7/fry1btkReAwAAs5vlYBIOEz6fL2Hbrq4ubd++XZK0Z88erVmzJhJmNm3aFHns9/u1efNmtbe3p9J3AABQYIxkz8oxDEPxPsXn82nNmjXq7++XFAofK1as0LFjxyRNDSaSNG/evEjbRILBoE6fPq25c+fKMIxkug0AAHLENE1duHBBixcvlssVfxVJUlM5VjQ0NGjfvn2R5wMDA5Ikj8ej/fv3y+PxTGnv8Xjk8/nU0NAw42uNjY1pbGws8vzUqVOqr6+3u8sAACALTp48qZe97GVx29geTCRp48aNkcff+c531NLSIrfbHQkp0/X19UW9vmvXLrW1tc24fvLkSVVXV9vSV8Buw8PDWrx4sSTp9OnTqqyszHGPUIi4z+AkgUBAS5Ys0dy5cxO2zUgwCRsYGFBHR0fCxbKxAsuOHTt09913R56H/8Oqq6sJJshbRUVFkcfV1dX8wkBGcJ/Biawsw8hoMGltbdWhQ4ciu27cbveM0ZG+vr6Yu3LKyspUVlaWyS4CAIA8krE6Jnv27FFra6u8Xq8GBgY0MDCglpaWqG2bmpoy1Q0AAOAgKY2YDAwMTBnl8Pl8crvd8nq9kkJbixsaGiKhZP/+/VHrlfj9fjU1NVHHBAWluLhYH/jAByKPgUzgPkOhsrxduKurS4cOHdKePXu0fft2NTc3Rxa5btq0Sc3Nzdq+fXtke/Bkbrd7yvbhvXv3qrm5Wd3d3dqxY4flYBIIBFRTU6PBwUHWmAAA4BDJ/P5Ouo5JLhFMAABwnmR+fzP+B9jMNE2NjIxIkioqKigGiIzgPkOh4hA/wGYjIyOqqqpSVVVV5BcHYDfuMxQqggkAAMgbBBMAAJA3CCYA4HA/ePa0njx2XhNBx+xlAGJi8SsAONCh585EHm9rf1au0nLV1ZRr54Z6rVtVl8OeAelhxAQAHKbzaK8+/m/PzLh+ZnBUdz3sU+fR3ux3CrAJwQQAHGQiaKrtYI+iTdqEr7Ud7GFaB47FVA5gs6KiokhV5MknwAJ2OHy8T72DozJcLlW86vWSJMN19W9MU1Lv4KgOH+/TTSvm56iXQOoIJoDNysvL1d7enutuoECduzAqSTKKS7XwHTsStgOchqkcAHCQa+aW29oOyDcEEwBwkNXLPaqrKVesAvSGpLqacq1e7slmtwDbEEwAmw0PD8swDBmGoeHh4Vx3BwWmyGVo54Z6BS+N6sTut+vE7rcreOnK9M6VNjs31KvIxdk5cCaCCQA4zLpVdbrv3TfOuF5bU64H7mygjgkcjcWvAOBAa1fWRh5/YdMNuu4aj1Yv9zBSAscjmACAw739hsWqrKzMdTcAWzCVAwAA8gbBBAAA5A2CCQAAyBusMQFsVlRUpLe97W2Rx0AmcJ+hUBmmaTrmpKdAIKCamhoNDg6quro6190BAAAWJPP7m6kcAACQNwgmAAAgbxBMAJsNDw+rsrJSlZWVlKRHxnCfoVCx+BXIgJGRkVx3AbMA9xkKESMmAAAgbxBMAABA3iCYAACAvEEwAQAAeYNgAgAA8ga7cgCbuVwu3XLLLZHHQCZwn6FQUZIeAABkFCXpAQCAIxFMAABA3iCYADYbHh7WwoULtXDhQkqFI2O4z1CoWPwKZMBLL72U6y5gFuA+QyFixAQAAOQNggkAAMgbBBMAAJA3CCYAACBvEEwAAEDeYFcOYDOXy6WmpqbIYyATuM9QqChJDwAAMoqS9AAAwJEIJgAAIG8QTACbjYyMaNmyZVq2bJlGRkZy3R0UKO4zFCoWvwI2M01TJ06ciDwGMoH7DIWKERMAAJA3CCYAACBvEEwAuwUnrj5+/mdTnwMA4iKYAHbqOSB9pfnq83/dKH15Veg6ACAhgglgl54D0v73Sxd6p14P9IauE04AICF25QB2CE5Ina2STBmGVL8wlPkNQ5JMSYbUeY90/XrJVZTDjqJQGIah+vr6yGOgUBBMADuceEIKnJYkVZQYeu4jVdMamFLgVKjd8jdmv38oOBUVFXruuedy3Q3AdklN5fh8PjU2NqbV1ufzyefzSZL8fn/kMeBoQ2ftbQcAs5TlYNLR0SFJloJEvLZ79+5VY2OjDMPQ1q1b5fV6rXYByF9Vi+xtBwCzlOWpnI0bN1r+ovHaNjY2qr+/X5Lkdrstf00grw2fjzwcGTfVvG9YktS9uVIVJVfm/6uvlZbenIveoQCNjIyouTm0A6y7u1sVFRU57hFgj5ysMSGQoKAEJ6R/3xF5appSz4vByOOI2z7PwlfYxjRN9fT0RB4DhSLrwWRgYCAy1dPd3R13OmdsbExjY2OR54FAICt9BJIyaeFrXBXzM98XAHC4rAeTLVu2REZMvF6v1q5dq2PHjkVtu2vXLrW1tWWxd0AKrC5o/X+PsiMHABLIeoE1v98feez1euX3+6dcm2zHjh0aHByMfJw8eTJb3QSss7qg9dn9lKcHgASyGkx8Pp/WrFkz47rH44navqysTNXV1VM+gLyz9GZr0zQjL4WmfQAAMaUUTAYGBqY89/l8MUc9Jrf1er3avXt35HlXV5c2btzIYlg4m6tIuuFd1tpSxwQA4rK8xqSrq0uHDh2SFFr70dzcHNkWHH6+ffv2uG3dbreampq0Z88eud1uHTt2TO3t7Xb/NwHZ96q3ST+/X1KoDP3SGiPyeArqmMAmhmFo6dKlkcdAoTBMB+0zCwQCqqmp0eDgINM6yC+PfVp68ivx21RfK33iV2wZBjDrJPP7m9OFgXRdvhQZLYnrLX9DKAGABAgmQLq690lmMHG7C72Z7wsAOBzBBEhX//NTnl4cN9W8b0jN+4Z0cdyM2Q5Ix8WLF9Xc3Kzm5mZdvHgx190BbJOTkvRAQZm3bMrToCk9fToYeRyrHZCOYDCop59+OvIYKBSMmADpat4sGYn+KblC7QAAcRFMgHQVl0qvXJegUVD6TWdWugMATkYwAdIVnJD8P07crvMeStIDQAIEEyBdj++RxocTtwucoiQ9ACRAMAHSEZyQfvoF6+0pSQ8AcbErB0jHsR9J5szpmQUVMUqEU5IeNlqwYEGuuwDYjmACpOPQzhmXKksNvbht7sy2FQtCJxEDNqisrNSLL76Y624AtiOYAKkKTkjneqy3X/9/KEkPW00ETR0+3qdzF0Z1zdxyrV7uUZGLA/3gbAQTIFXP/5cki4WtFq2SVr4jk73BLNN5tFdtB3vUOzgauVZXU66dG+q1blVdDnsGpIfFr0Cqjv806uWL46be9PVhvenrw1dL0v9pVxY7hkLXebRXW//lSf3igU/ozLfuUXB8TJJ0ZnBUdz3sU+dRzmWCczFiAqTMjHo1aEo/PjEReSzPCql0Thb7hUI2ETTVdrBHpmlq7OTR0EUzdC+akgxJbQd7tLa+lmkdOBIjJkCq5rittWv6XxntBmaXw8f7pkzfTGdK6h0c1eHjfdnrFGAjggmQqsqF9rYDLDh3IXYoSaUdkG8IJkCqhl+ytx1gwTVzy21tB+QbggmQqhGLQ+VW2wEWrF7uUV1NuWKtHjEU2p2zerknm90CbEMwAVJ1/CfW2rn4Zwb7FLkM7dxQH/W1cFjZuaGeha9wLN4xgVRcviSd6o75ckVJ6EOStPQN2ekTZo11q+p037tvlKukXEZJWeR6bU25HrizgTomcDS2CwOp6N4X86XKUkPDn6oOPSmpkJa/MUudwmzyjtUrdGl0hMqvKDgEEyAV549Za3fdTZShR8YUuQzdtGJ+rrsB2IqpHCAVFyxW1iwuS9wGABBBMAFSMbc25kujl02t/9aI1n9rRKNlHEuPzBgdHdX69eu1fv16jY5SswSFg6kcIBXzXx7zpYmg9OhvL4ceu5dnq0eYZSYmJvToo49GHgOFghETIBXNmyXDwj+fRsrRA0AyCCZAKopLpZv+3Fo7AIBlTOUAqXrL50L//+RXJDN49brBLhwASBXBBEjHWz4nvfkzobom/c9L85ZJ9e+RPk85cABIBcEESFdxqXTTR68+Hx7OXV8AwOFYYwIAAPIGIyaAzSorK2WaZq67gQLHfYZCxYgJAADIGwQTAACQNwgmgM1GR0e1adMmbdq0iVLhyBjuMxQqw3TQJGUgEFBNTY0GBwdVXV2d6+4AUQ0PD6uqqkqSNDQ0pMrKyhz3CIWI+wxOkszvb0ZMAABA3iCYAACAvEEwAQAAeYNgAgAA8gbBBAAA5A2CCQAAyBuUpAdsVlFRoaGhochjIBO4z1CoCCaAzQzDoKYEMo77DIWKYAIkIzghnXhCGjorVS2Slt4suYpy3SsAKBgEE8CqngNSZ6sUOH31WvViad1uqf72yKWxsTFt3bpVkrR3716VlZVlu6eYBbjPUKgoSQ9Y0XNA2v9+SdP/uRih/7vjoUg4oVQ4soH7DE5CSXrATsGJ0EjJjFCiq9c67wm1AwCkhWACJHLiianTNzOYUuBUqB0AIC0EEyCRobP2tgMAxEQwARKpWmRvOwBATAQTIJGlN4d234QXus5gSNXXhtoBANJCMAEScRWFtgRLmhlOrjxfdy/1TADABkkFE5/Pp8bGxrTa+v1+7dmzRx0dHdqzZ48GBgaS6QKQG/W3h7YEV9dNvV69eMpWYSlUHvzcuXM6d+4cpcKRMdxnKFSW65h0dHTI6/WqsbFRiT4lXtvGxkYdOXJEUiiktLa2qr293VJnqWOCnKPyKwAkLZnf30kXWDMMI2EwidXW7/dr06ZNkWAiSfPmzVN/f7+lr0cwAYCrJoKmDh/v07kLo7pmbrlWL/eoyBVrLRSQO8n8/s5qSfquri55PJ4p1zwej3w+nxoaGrLZFSBjxsbGdPfdd0uSvvjFL1IqHBlxwPe8tnzk4xq+NCHPm/9MRnGJ6mrKtXNDvdatqkv8BYA8ldXFr7HWk/T19UW9PjY2pkAgMOUDyHeXL1/W/fffr/vvv1+XL1/OdXdQgDqP9uovHj6is08d0NAvHpF5perwmcFR3fWwT51He3PcQyB1ebErJ1Zg2bVrl2pqaiIfS5YsyW7HACDPTARNtR3siXdAgtoO9mgi6Jhj0IApshpM3G73jNGRvr4+ud3uqO137NihwcHByMfJkyez0EsAyF+Hj/epd3A05uumpN7BUR0+Hn0kGsh3WQ0mLS0tUa83NTVFvV5WVqbq6uopHwAwm527EDuUpNIOyDcpBZPpUy8+n09+vz9hW6/XO+U1v9+vpqammCMmAICprplbbms7IN9YDiZdXV1qbW2VFFr70dHREXlt+vN4bdvb29Xa2qqOjg7t3bvXcg0TAIC0erlHdTXl8Q5IUF1NaOsw4ERJ1zHJJeqYwAmGh4dVVVUlSRoaGlJlZWWOe4RC03m0V1u/9oRe+NJGSdKST3bIVXo1rDxwZwNbhpFX8raOCTAbzJkzR8ePH488Buy2blWdHvjg7+vTxQ/rxQtjMkpKJUm11DFBAWDEBAAcisqvcApGTABgFihyGbppxfxcdwOwVV4UWAMKyaVLl7Rt2zZt27ZNly5dynV3UKC4z1ComMoB4knhNGEWvyIbuM/gJEzlAHboOSB1tkqB01evVS+W1u2W6m/PXb8AoIAxlQNE03NA2v/+qaFEkgK9oes9B3LTLwAocAQTYLrgRGikJN4xaZ33hNoBAGxFMAGmO/HEzJGSKUwpcCrUDsgDP3j2tJ48dp4ThVEQWGMCTDd01t52QAYceu5M5PG29mflKi1XHQXWUAAYMQGmq1pkbzvAZp1He/Xxf3tmxvUzg6O662GfOo/2Zr9TgE0YMQGmW3pzaPdNoFfR15kYodeX3hz10+fMmaOjR49GHgN2mgiaajvYI5WUqu5DX5WkSEl6U6FD/NoO9mhtfS1VYOFIjJgA07mKQluCJWnGGa5Xnq+7N2Y9E5fLpZUrV2rlypVyufgnBnsdPt6n3sFRGYZLpQuXqnThUhnG1fvMlNQ7OKrDx/ty10kgDbxrAtHU3y7d8ZBUPW2uvnpx6Dp1TJAj5y6M2toOyDdM5QCx1N8uXb8+6cqvly5d0uc//3lJ0qc+9SmVlpZmo7eYJa6ZWy5JMifGNfjkfklSzU13yCgqidoOcBqCCRCPq0ha/sakPmV8fFxtbW2SpG3bthFMYKvVyz2qqynX6RdHNfizb0uSqlf/cSSYGJJqa0InDQNOxFQOADhIkcvQzg31UV8Lr4jauaGeha9wLIIJADjMulV1uu/dN864XltTrgfubKCOCRyNqRwAcKC1K2sjj7+w6QZdd41Hq5d7GCmB4xFMAMDh3n7DYlVWVua6G4AtmMoBAAB5g2ACAADyBlM5gM3Ky8t1+PDhyGMgE7jPUKgIJoDNioqK1NzcnOtuoMBxn6FQEUyAsOBE0lVeAQD2IpgAktRzQOpslQKnr16rXhw6zC/Jc3EuXbqk++67T5L08Y9/nMqvyAjuMxQqwzTNaOe656VAIKCamhoNDg6quro6191Boeg5IO1/v0Lnsk52pR5Ekof2DQ8Pq6qqSpI0NDTENk5kBPcZnCSZ39/sysHsFpwIjZTMCCW6eq3znlA7AEDGEUwwu514Yur0zQymFDgVagcAyDiCCWa3obP2tgMApIVggtmtapG97QAAaSGYYHZbenNo941iHXxmSNXXhtoBADKOYILZzVUU2hIsaWY4ufJ83b3UMwGALKGOCVB/e2hLcNQ6JvcmXcekvLxcP/rRjyKPgUzgPkOhoo4JEEblVzjMRNDU4eN9OndhVNfMLdfq5R4VuWJNSwK5k8zvb0ZMgDBXkbT8jbnuBWBJ59FetR3sUe/gaORaXU25dm6o17pVdTnsGZAe1pgANhsfH9dXv/pVffWrX9X4+Hiuu4MC1Hm0Vx/+xmH95kcduuD7gcyJy5KkM4OjuuthnzqP9ua4h0DqmMoBbEapcGTSRNDUG3b/p069OKCTX9ooSVryyQ65SkPrTAxJtTXl+q/WNzOtg7xBSXoAKFCHj/dNmb6ZzpTUOziqw8f7stcpwEYEEwBwkHMXYoeSVNoB+YZgAgAOcs1ca1uDrbYD8g27cjD7sC0YDrZ6uUd1NeU6/WL0EZHwGpPVyz3Z7RhgE4IJZpeeAzEKqe1OupAakAtFLkM7N9Rr69dmnngdXuq6c0M9C1/hWEzlYPboOSDtf//UUCJJgd7Q9Z4DuekXkKR1q+p037tvnHG9tqZcD9zZQB0TOBojJpgdghOhkRJF2x1vSjKkznuk69enPa1TVlamH/zgB5HHQCa8vWGpvn/goH5z9oK8Nzapbl4VlV9REAgmmB1OPDFzpGQKUwqcCrVLs/prcXGx1q9fn9bXABIpLi7W7RvenutuALYjmGB2GDprbzsgT3BeDgoNwQSzQ9Uie9vFMT4+rn/913+VJL33ve9VSUlJ2l8TmG58fFz37P6qvvfMaV1e/noZRaG3c87LgdNRkh6zQ3BC+vKq0ELXqOtMjNDunE/8Ku01JpSkRzZ87/AxvfN1L5c0syS9JBbBIq9Qkh6YzlUU2hIs6epbt6Y+X3cv9UzgCBNBU5//4X9HfS0cu9sO9mgi6Ji/O4EIgglmj/rbpTsekqqn/RVZvTh0nTomcIjDx/t0ZnAs5uuclwMnY40JZpf620Nbgqn8CgfjvBwUMoIJZh9XUdpbgoFc4rwcFDKmcgDAYVYv96i2JnbxPkOh3TmclwMnIpgAgMMUuQx96q2vjvoa5+XA6ZKayvH5fNq8ebOOHDkSt53f71dHR4e8Xq/8fr+2bNkit9sd+RqS1NDQIL/fr4GBATU0NKTWeyCRHJwkXFZWpv3790ceA5nw9oal+tTf7VX7kd9prPhqrZxa6pjA4SzXMQkHjcbGRiX6lMbGxkh48fv9am1tVXt7uyRp69atevDBByVJLS0tam9vj4SWRKhjgqRwkjBmASq/wgmS+f1tecRk48aNltr5/f4pz71er7q6uiLPGxsb1d/fL0mWAwmQtPBJwtOLqYVPEmZ7MApEkcvQTSvm57obgG1s35XT1dUlj2fqgiuPxyOfzxeZsrEaSMbGxjQ2dnWvfiAQsK2fKGBZPEk4msuXL+u73/2uJOmd73yniovZ/Ab7he+zYNDU4htvUd/Fy4yYoCDY/o45MDAQ9XpfX1/k9Y6ODklSd3e3tm7dKq/XG/Vzdu3apba2Nru7iEKXxZOEoxkbG9Mdd9whKVSSnmCCTJh8n00uSc9ZOXC6rL1jhgPL5IWwXq9Xa9eu1bFjx6J+zo4dO3T33XdHngcCAS1ZsiTTXYXTcZIwZoFDz52Jev3M4KjuetjHWTlwLNu3C7vd7sjoSFhfX18kjExegxLetTN9XUpYWVmZqqurp3wACWXxJGEgFzgrB4XM9mDS0tIS9XpTU5N8Pp/WrFkz47Xpa1KAtCy9ObT7ZsZhfWGGVH1tqB3gQJyVg0KWUjCZvo7E5/NFRj2mrxfx+/1qamqS2+2W1+vV7t27I691dXVp48aN7M6BvThJGAWOs3JQyCyvMenq6tKhQ4ckhRalNjc3R7YQh59v375dktTe3q7W1lY1Nzeru7s7UsPE7XarqalJe/bskdvt1rFjxyKvAbYKnyQctY7JvWwVhqNxVg4KmeUCa/mAAmuIKVaF1xxUfh0eHlZVVZWk0K6cysrKjH4/zD4TQVM3fe4RdX92g6Spu3Kk0LhgbU25/qv1zWwdRl7ISIE1IG8lqvCa5ZOES0tL9S//8i+Rx4DdilyG/ur2G/TBw5+QJBlFV9/KOSsHTseICZwtVoXX8NszFV5RwDqP9qrtYI96B6+uJaGOCfJRMr+/CSZwruCE9OVVcYqpGaGRk0/8ioWuKFiclQMnYCoHs0OOK7zGcvnyZT322GOSpNtuu43Kr8iImfcZ5+WgMPCOCefK0wqvY2Njevvb3x761pSkR4ZMvs/+49kTujBRxIgJCgLvmHAumyu8MiQOp/rgv3RzVg4KBsEEzhWu8BroVfSThK+sMbFQ4ZVFhHAazspBoSKYwLnCFV73v1+hXTiTw0nsCq/TR0b6hy/po9/yzYg2ZwZH9eGHffpkyyu0bEEloyjIG4nOyjEUOitnbX0t9ysch2ACZ0uywmu0kRGXEX28JXztS12/jVxjFAX5IJmzcm5awaJYOAvBBM4SrZJr/e3S9esTVnjtPNqrux6eOTKSzAGsDJMjH3BWDgoZwQTOkUaF14mgqbaDPVFHRpLBMDnyAWfloJARTOAMsSq8BnpD12NUeA2vJ/nZ/7w4ZfomHYmGyUtLS/WVr3wl8hiw2+rlHtV5qjS89sOSppakl66elbN6uScHvQPSQzBB/gtOhEZKYq4EMaTOe0LTOZOmb6KtJ7FTrGHykpISffSjH83I9wSk0Fk5be+4QXcNXZYUddk3Z+XAsVy57gCQUDIVXq8IryfJVCiRGCZHbq1bVactf7BcxrTsYRjSlj9YzhooOBbBBPkvyQqvqa4nsfrHpaHQ7pxYw+QTExN6/PHH9fjjj2tiYiLJXgDWPPLL3+m+h7+vkRPPygxevc+CpvTgT46r82hvDnsHpI6pHOQ/i5VbnwvM0f88c0ovXRhLaqQknEe+8iev1bzKMp27MKrnXxrRl7t+Iyn5YfLR0VHdeuutkkIl6SsrKy33BbBiImiq7XvP6Oy3PyVJWvLJDhmlU3ehsUAbTkUwQf6ZviV4yeviVng1ZeisPNpwMKignkn629XGqE3yqtqqGWtUYrUFsok6JihkBBPkl1hbgldtlJ74B02v8GrKkGma2jn+PgWTnJn881tfrte/fEHMaq7rVtVpbX0t5+cg71DHBIWMYIL8EW9L8BP/IN38F9LRjimh5aw82jn+Pj0WXG3524S3Un5y7SsThowilzHlL86JoKknj50nqCCnqGOCQkYwQX6wsiX46P+VPvZL6eRT0tBZPReYc2X6xvpISTpbKTnoD/li9XKPamvKdDLG69QxgZOxKwf5weKW4IkXfq4ng/X6/sRN+uHQy5OevqmtKU+pnHys7cfhEvXsgEA2FbkMfeqtr475uinp9tfUMZoHR2LEBPnB4pbgtm/9px4aupDUl/7M+ldrwdyylKde4m0/pkQ9cmXtytq4rz/4k+N67XXzGM2D4xBMkBvTd95ULrT0ab8Zsb71Njyc/cHXL08rMBw+3hd3+/H0HRAlJSXas2ePpFAVWCATXEXFum7dZgVGL8soKorahsAMJyKYIPui7byZWyfN8UgX+xVtnUlQ0hlzvg4Hr7f0Lewsy53sDojS0lJt27Ytre8JJPLMqSEZr/lD1cR4nS3DcCqCCbIr1s6bC2ci10wZMqZtCZZpqi2JLcF21hthBwTyEVuGUagIJsgeCztvxkprNHCpSIt0PvLKGXn0WQtbgv/81hV6xaK5tm/jXb3co7qacp0ZHI3a8+k7ICYmJuTz+SRJDQ0NKooxzA6kY35FicZ6Q9WJSxetkOGKfp8RmOE0BBNkj4WdN2WXBvSJS5+SKZeu0YDOya3DwestjZS8/uULMzJkXeQytHNDve562DetvFv0KaPR0VGtXh0KUZSkR6bcUFehMw/dLSl6SXop/plOQL4imCBzpi9wvWBtS+1CBXQgeLPlb5ONmg3rVtXpgTsbKFGPvGFlRJAtw3AiggkyI9oC1wproxnn5Lb8bexc5JoIJeqRTyaCic/PPvDLXm1f92ruUTgKwQT2i7XAdSS0biRc+2O6oCmdUfydN+45JRq4OB55nu0Ri+kl6oFcefr5voRt2JUDJyKYwF5xF7heuWqGtv9O/iMu/Mdfop03X31Pg1wugxELzHovDsU+XXgyduXAaQgmsFeCBa7Glf/pM+dqvq5WcD2j+WqLs/MmvI7k91fMJ4gAkhZWlVlqx64cOA3BBOlJcYHrX4+/T2flsbTzJpvrSACnaFqWeLE3u3LgRAQTpC6NBa5n5dHPg/VRX8v1OpJ0lZSUaOfOnZHHQCaUl5Xqze/5iI6c6I9Zkp5dOXAigglSE2eBa3gdiRHl/dDKAlenryMpLS3VZz/72Vx3AwWuqLhEI6v+SO4lsdeQsCsHTkQwQWzTp2mW3iy5iiwtcDUUCiHJLHBlHQlgXaLDJSV25cCZCCaILto0TfViad1uac68uAtcXVLSC1wLaR1JMBjUf//3f0uSXv3qV8vlsna+D5CMM4MjuvTiCUlSyYIlMozo9xm7cuA0BBPMdGWaxgwdnxdhBnpl7H+/9Pt3Wfoy8Ra4On0dSTwXL17UqlWrJFGSHplTU2Kq92sflRQuSR999w27cuA0BJPZLNpUjSR1ts4IJZJkhK8+uz9qgbTp4i1wdfo6EiDXrOzKmVdRwq4cOA7BZLaKNVXT8EEpcDpm8DBkSiMv6bw5V/N0QdGyRLwFrqwjAbIncdF6IP8w+T0bhadqpq0TMQO9Mh//vKUv8d3Lb5B0dUFrWLwFroW0jgTINSsl6QdGxnX4eOJ2QD4hmBSq4IR0/KfSrzpC/x+cuHo9zlSNVV1mo+4a/4TOaOow8RnN113jn9BjwdXyVJZOea22plwP3NlQEOtIgFyjJD0KFVM5hajngMzOVhmTRkTM6sUyJu2oiT1VEzJ9q68mXT+j+fpN6Sr1jwZ1aKxJq12/nrLA1ZRLdTXl+vG2W3XkRD/rSIAMoCQ9ChXBxGli1RYJ6zkgM+qOmtOh6ZvXfdjSMFmiOiQfePPL9eWu38iUa8oC18nTNaXFLuonABny2uvmJWzjMqTGpYnbAfmEYJIvEgUOKf5ISP3tUnBCFw9uU5lpzhjtcEkKmqZGj/ybKix054uXN+pPiv9Ti3V1fvqM5uuvx9+nZ+f+ge5/88v1qtoqtR3smVLkqZC2/aaqpKREf/mXfxl5DGTCr04PqXr1H0lSzJL0QVM6cqKfPxDgKASTfJAocITbxBkJMe54SBPlbs25eEax5mlchlRxud/SjpqvTrxD90+8Q82Tpmm6r9QheeDK4tV1q+q0tr5Wh4/3MV0zSWlpqb7whS/kuhsocP1jQc279UMJ27HGBE5DMMmwicuX9eunHtPF/lOaM+9aXf+621RUPOnHbiFw6Pr1iUdCDm7TyYbteqWFPn338hv0oeIfxp2q+XjL9fq37hf088Gr0zR1UUZDilwGf40BOWB17cjzLw1nuCeAvQgmKUgYNq74xWPf0OIn27RS5yPXzh6ar9M37dRrb/uAtamXg9tUWlqdcCRkzsUzGh08a6n/Py95nZ6+9Cr9VclDcadq/vzNL2c0JAXBYFAvvPCCJOm6666jJD0yommpW/MmBvTi0JiKqhfGLEn/7cMv6M/f/Ar+7cIxCCayHjQkC2FjUrvXPPGx0JNJ7wcLzfNa+MTH9AtJN7xiuaXA8cIz/67rLPx3nA/O1WnTo1r1xZ2mueENb9WX/uOYDo01xZ2qkcRoSAouXryo5cuXS6IkPTLn0tionvm7OyXFL0l/JjDGQX5wlFkfTKwGjXDbRGHjtbd9QBOXL2vxk22SZm65dRmhgFD3ZJv+x/i0XmWhj6f6L1oKJsXua/X3JX+mz4/viTlN8/clf6q/XfMqvbKuRm0HexJO1SCxS5eD+uaTz+tE34iWeir0RzcszHWXgClYZwInmdXBxGrQkGQ5bEysea9+/dRjoaATZySkVufVc+qkpWBysrpBS3/3vYQjIcXL36A3Lb5JH/nWpZjTNO/Y9CEWrtpo16M92vfT41Mq4H7ue7/IXYeAKBZYrHkC5INZG0ySCRpFxcWWw8ZzV6aErBgpqbE09XLtjW/R3//2hcQjISsWhoLFez6sTQderyVDv4xM05yseo0+s+n3WLhqo12P9mjvT47PuD69TD+Qc9yTcJBZG0ySCRorX7/ectgIr1OxYumyV+jvn7cw9fKKa3ThHR+yNBIiadJoSCOjIRly6XJQ+346M5REa8cKE+TaS8PWytcD+WDWBpNkgoYky2EjvHj27KH5WmiejzkScs6Yr/rfX6c3uRstT71YHQmRGA3JtG8++bylkZFvP3VCH3nLqsx3CIiDqRw4SVLBxOfzafPmzTpy5Ejcdn6/Xx0dHfJ6vfL7/dqyZYvcbnfC17IpmaAhyXLYCO/oOX3TTi184mMxR0J6b9qp2uLipAIHIyH54/h5a7Uhnu+jhgTyAFM5cBDLwSQcJnw+X8K2mzZtioQXv9+vzZs3q729PeFr2ZRM0JCUVNiQQotmfyFp8ZNtWjRpx885Y756p+34SSZwMBKSH+JFQcNVpKrXrpckFRXN2kFJZFhxcbHeescH9NPfviRj+vEV0zCVAyex/K65ceNGS+38fv+U516vV11dXQlfy7Zkg4aUXNgIt59Y897IgtjwNE9tlBopBA5nee2Sefrmz1+I+ppRXKL5b7lLktTsXZTNbmEWKSsr06f/9u/03n96KmFbpnLgJLb/OdfV1SWPxzPlmsfjkc/n09NPPx3ztYaGBru7klCyQSP8OVbDhhQKQCtfvz5T/wnIkUXV1sqBW20HpMTqFA1TOXAQ24PJwMBA1Ot9fX1xX4tmbGxMY2NXhyADgUC63Zsh2aAhETaguHM5pmkqeDF0r5r8RkCGmKapYydPa2JkUK451TKM2DclUzmzUzJVzVNpnylZ+46xQkm813bt2qW2trbMdGgSggaS9dJQ7Dd6c3xMv/uH90qSOt94WG94BZVgYb+RkRF9YM1rJMUvSS8xlVNIbDurLc32mWR7MHG73TNGQPr6+uR2u+O+Fs2OHTt09913R54HAgEtWbLE7i4DSbN6susPnu1V2x+b7JxCTn37qef1+pcvyHU3EIPdYSOZquaptM802489bWlpiXq9qakp7mvRlJWVqbq6esoHkA9WL/fIU1mSsF3f8LgOH48+VQlkyyO/OqtLl4O57sasM3H5sp772SN6+gcP6rmfPaKJy5dntPnFY9/QS3/zSq089B41Pb1NKw+9Ry/9zSv1i8e+MaPda574mBaa56dcX2ie12ue+FikfaKq5tKVquZX+pJs+2xIKZhMn3rx+XyRHTder3fKa36/X01NTXK73XFfA5ykyGXonTdaq4XDAWrINVOhooCwz0TQ1JPHzuv7z5zSk8fOa2JaxUUrgSMTYePXTz2mRYpeBiPcvlbn9eunHpOkpNtng+WpnK6uLh06dEhSaO1Hc3NzZAtx+Pn27dslSe3t7WptbVVzc7O6u7un1CmJ9xrgJG++fpH++WfPJ2y3oJL5feTeib6RXHfBMSaCZtwDTjuP9upzB341syjm7aGimFamRm5Y817bD4ZN5qy2cLtk22eD5WDS0tKilpYW7d69e8Zr08OF1+uNtJte/yTea4CjWF02wvIS5IGlnopcd8EREoWOzqO9+t63/lHtJQ9pcenVadrTYx799bfer+C7NqvRQuDoqZin38tA2Ei2qnmy7bOBspRAiuLtzEmlHZCqRNnXZUjvu2lZNrqS9+KNhiQMHe/+sH5y4Gu6v+TLM75urfp0f8mXtfd7p/U2C4Hj+K//01J/kw0byVY1T7Z9NhBMgBTFmqIxXEWqXLUm8pipHGRCcXGxPvCB0E6JFbes0D8/Gfuv6s1vXK7SYtv3OuSlRMEj1mjI2vpaPf69+KHjnu+69Ing1yXFHg15d/BRi6Ok1oZSkw0byVY1T6UKeqYRTIBUxXhfMYpLtGD9JxO2A9JRVlamr3/965HnxSWl2vfT41NOvXYZoVCy42312e9gDnQe7VXbwR71Dl5dcF5XU66dG0L//fFGQ35z65/oY+P/JCl26Nge/CctMGIX+nQZ0jwNWepr1fW36uyp72ckbKRyfEqyVdAzyTBN0zGlKQOBgGpqajQ4OMjWYeTc9585pY//2zMJ29337hv1hxZ38ADpuHQ5qG8++bxO9I1oqadC77tpWcGNlMQaEek82qu7HvbJUFCrXb+OjIh0B69XUC69c45Pfxf8O0nRf7H/o3GHPqL9tvRxQFWqNofiBo6F//s3evY//jWySDZan35589/PqDcyPTycUezwkE+VX5P5/c2ICZCiWFM0pmnKHA+tKzFKypjKQUaYpqmRkdBOm4qKChmGodJil/70jd4En+lcsUZEPrO+Xp97pEdvcR3WzpKHtNiYNCJievS58ffpL4PflBR7NORPTKtTMImdffUHVd3zlYSjG5k8GFZKvqp5vlRBZ8QESNHP/uelqCe7Bi+N6uSXQjvOlnyyQ9/+yC1U3YTthoeHVVVVJUkaGhpSZWVljntkn2ijIod6zsQcEZmQS7e5DuuBK+tDpocBQ1Kco4SSMlY6TyVj/TFHQ8YqajVnW49+cehhy6Mb+XJGTSYxYgJkAbtyAPtFGxWprS7X6OWJuCMinymJPSJi9c/vi0XVKrsciB861u+W2fFBBWVOqVAalGQYhuZs+ILkKkpqdCNfRiryBcEESJHV83KstgNmu/A6kek54kxgdMqIyGShHTP3xR0RsTpaUvqGj8r48b3xQ0f97TJcD0mdrVLg9NXvUX2tjHX3SvW3R64ROFJDMAFS1Lh0XmR+OhaXEWoHYKrp0zWNS+ep7WCPTEmuadM1TwdfqZ0lD0lKb0Rk+pqPydfHKmo155Zt0qL6xKGj/nYZ16+XTjwhDZ2VqhbJWHqz5CpK4SeB6QgmQIqOnOiPG0qk0BvekRP9umnF/Ox0CnCAaNM1nspS9Q1f0m1RpmvOm3M137gQ8+tZHRExjNDoR7wpGMuhw1UkLX+jtW+MpBBMgBRZPZzvUM8ZgglmrekjI/3DY/rot34xY7omHEqiTdd4FDuUTGYq+sYaU4aM6sUy3vJ56d93JJyCIXTkFsEESJHVtSPff+a0Pr2+fsohYMBsEG1kxGUo6ekayyMiuhJCJsWe0HNJV8KHUb+BKZg8RzABUrR6uUeeyhL1DY9PuW64XKp41esjj88PX9Lh432MmsBWRUVFkYNQi4ry7xdrrIWsQVMpTdfEEx4R0Vs+L2PGiMjiSCiRxGiIAxBMgBQVuQy988Zr9c8/e37KdaO4VAvfsWPKNavTPoBV5eXlM052z6XJUzYLqsr02QPPzQglkmyfrpk+IqJpIyJiRMRxCCZAGlrqa2cEk2jYMoxCFm3KJmzylM2LqtbOkm+Erqc6XVMxXxq5WrSMEZHCQzAB0sCWYcx2saZspOhTNqkzpOrF0seekU4+xYhIASOYAGmItmV4ekl6lZazZRi2s1KSPtaBd+lId8omNVf6vO5eqbiUEZECRzAB0mB17QhrTJBtsQ6827mhXutW1dn2NaNxKRhzh40l06ZrNH26BgWNYAKkgbL0yEcxS7sPjuquh3164M6GpMNJvCmb6Va7fp3i9A3TNSCYAGlhjQnyzUTQjJR2ny68o6XtYI/W1tfGndaxOmUTNnmR68uN36XQc6ZrEEIwAdJAWXrkm8PH++JOtZiSegdH49bWsTplE2bLIlema3AFwQRIA2tMkG+s3ms/PNorKVQoUFJkdOT5l4b1pa7fWv5+sRa5hg/Wi74N+MqUzR/eL428xHQNpiCYAGlgjQnyjdV77aEnT+ihJ0/IXVEiSRoYGU/wGTPFW+RqXCk9P9OkKZsVb0r6e6LwEUyANERbY2K4XJrjbYo8Zo0JMqGoqEhve9vbIo/DVi/3qK6mXGcGRy0tVE0lkES+V4JFrlEHS5iyQQIEEyAN0daYGMWlumbTZyPPWWOCTCgvL9cjjzwy43qRy9DODfW662HflUPt7JXSItc3bpOuuZ4pG1hCMAHSwBoT5KN1q+r0wJ0NSS1gtSLlRa7eW9hlA8sIJkAaWGOCfLVuVZ3W1tfq8PE+/fBorx568kRaXy+yyHXa/Mz0Q/WmurLIdenNaX1vzC6uXHcAcLLwfP5kwUujeuGLf6wXvvjHCl4K/bXaP3wpF91DARseHlZlZaUqKys1PDwctU2Ry9BNK+brrSlWeg2bssh12mtxQ4kUWk/C1A2SQDAB0lDkMvSZ9a+ecd0cH5M5PhZ5/rlHejSRqOAJkKSRkRGNjIwkbBcO0KmekhNe5JpUefnqxdIdD7HIFUljKgdI07zKsoRtEhW0AjIp2QWxhkKjJLubhrSs/IIWjfZLRy18Ixa5wgYEEyBNLICFE8RaEDuvokSmpm4bflfVM9pZ8pDmHD2T3DdhkStsQDAB0sQCWDjF5AWx5y6M6pq55TMqv17f/7he+eMvyLiczNQji1xhH4IJkKbw/H28bZl1NVd/AQC5FF4QO91NK+ZLwQnpy3+j5KqfsMgV9mLxK5CmIpeh218Tf9fD7a+pi3uSK5AXTjwhBU4n9zkscoXNGDEB0jQRNHXgl71XLxiGypasijyWpAO/7NX2da8mnMA2LpdLt9xyS+SxLYbOWmvHIldkEMEESNP0Y+ZdJWWqfc+9U9qwKwd2mzNnjh5//HF7v2jVImvtWOSKDGIqB0gTu3JQMJbeHJqaiVnxxJCqr2WRKzKKYAKkiV05cLzghHT8p9Jz35UaPqjoheZZ5IrsYCoHSNP0XTnBS6M69Y8fkiRd++GvyVVazq4c2G54eFjLli2TJD3//POqrKxM7Qv1HJA6W6cuep3jkWRKF/uvXqteHAolLHJFhhFMgDSFd+Xs/cnxyLXgxcCUNuzKQSa89NJL6X2BngPS/vdrxvbgcCB506ek+StY5IqsYioHSNOMXTlRHPhlL2flIL8EJ0IjJVFrlly55vuGtPKdoYWuhBJkCcEESNP0XTnRhHflAHkjYc0SUwqcCrUDsohgAqSJXTlwJKs1S6y2A2xCMAHSxK4cOJLVmiVW2wE2IZgAaWpcOk+J1rW6jFA7IOfCW4Mv9EoV8Qr+UbMEucGuHCBNR070a8q6VsNQae0rIo8lKWiG2lH5FXZxuVxqamqKPLYk2tbgqKhZgtwhmABpmr52xFVSproPfClhOyAdc+bMUXd3t/VPiLU1OBpqliCHCCZAmlhjgrwXd2uwJBmhaZ11u6S5ddQsQU6xxgRIU7jya5zTRaj8ityysjV45KVQKKFmCXKMYAKkqchlaOeGekmhEBIcH9XvHviQfvfAh2SOh6Zvdm6op/IrbDUyMqJly5Zp2bJlGhkZid+YrcFwEIIJYIN1q+r0wJ0Nqq0pl0xpInBOE4FzWlRdpgfubNC6VXW57iIKjGmaOnHihE6cOCHTTLBuhK3BcBCCCWCTdavq9ONtt+qeda+KXHvsE7cQSpB7S28OLWiNN+HI1mDkCYIJYJPOo7265Qs/0r2d/y9y7bYv/1idR+OfowNknKtIWrf7ypPp4YStwcgvBBPABp1He3XXw74ZZ+acHRzTXQ/7CCfIjXAxtV91SHPmSRu/LlVPG8GrXizd8RBbg5E3ktou7Pf71dHRIa/XK7/fry1btsjtdsdsu3fvXq1YsULHjh3Tjh07Im19Pp8kqaGhQX6/XwMDA2poaEjrPwTIlYmgqbaDPTHPaDUltR3s0dr6WhbAInuiFVOrXizdtiu0NXjobGhNCVuDkWeSCiabNm3SkSNHJIWCx+bNm9Xe3h617dq1a3XkyBG53W75fD61trZq7969kqS9e/fqwQcflCS1tLTE/BqAEyRzujCVX5EVsYqpBXql9g+GRkh+b2MuegYkZHkqx+/3T3nu9XrV1dUVtW34eniEpKGhIRJEJKmxsVH9/f3q7+/XoUOHYo66AE4wo6KrIZXMv04l86+bMp1/qOdMdjuGgmYYhurr61VfXy/DmHSjxS2mduVa5z2hdkAeshxMurq65PFMLRDl8Xgi0zKTDQwMRP0ak9u63W4CCQrC9IqurpJyLf6z+7X4z+6Xq+Tqa99/5rQmghbKgQMWVFRU6LnnntNzzz2nioqKqy9YKaYWOBVqB+Qhy1M5scJGX1/fjGvhtSNh4UASbjswMKCOjg5JUnd3t7Zu3Sqv1zvj64yNjWlsbCzyPBAIWO0ukDWrl3vkqSxR3/B43Hbnhy8xnYPMo5gaHC7tXTnRAovX69Xu3bv14IMPamBgIBJSwiMuW7Zs0caNG7Vx40a9613v0tq1a6N+7V27dqmmpibysWTJknS7C9iuyGXonTdea6ktB/kh4yimBoezHEzcbveM0ZG+vr6Y0zHbt29XS0uL/H6/WlpaJCkyKjJ5NCW8w2f6GhZJ2rFjhwYHByMfJ0+etNpdIKta6msjj4Pjozr9Tx/R6X/6iILjU4MIB/nBLiMjI1q5cqVWrlw5tSQ9xdTgcJaDSThcTNfU1BT1ut/vl9frjUzrNDQ0RHborFmzZkb76etXJKmsrEzV1dVTPoB8NOUgP1MaP/+Cxs+/EFlryEF+sJtpmurp6VFPT8/VkvSXL0lP/aN0zUpFX/xKMTXkP8vBZPoaEL/fr6ampim1SSaPejQ2Nkamefbu3avdu3dHvk74sRRaVLtx40YWwsLRph/kN1n4OQf5IaM6d0h/s1B67FPS/xyK3oZianCApOqYtLe3q7W1Vc3Nzeru7p5Sf2TXrl1qbm7W9u3bJUm7d+9WV1eX+vr6tGnTpsiIi9vtVlNTk/bs2SO3261jx45RxwQFIXyQ31/93yN6YdL12ppy7dxQz5k5yJyvvVU6/8vYr7/yNummv6CYGhzBMBMeS5k/AoGAampqNDg4yLQO8lbgwpBqqudKkv7j2RO6ZeUSRkpgu+HhYVVVVUmShnbMVWVpnHvMKJI+fUYqLs1S74Cpkvn9zVk5gM0mh5DXeecTSpB75oTUvS/XvQAsIZgANptcRO0p/3mKqiEz/uOvk2vf/3xGugHYjWAC2KjzaK9avvhjFVVfo6Lqa/TBr3frDbv/k9OFYa9//4yMpx7Q0hpDS2sMGVYG5eYty3SvAFuwxgSwSefRXt31sG/GJs3w74wH7mxgASzSd/mS9LeLJDOY3Of97xdZY4KcYY0JkGUTQVNtB3viHZumtoM9TOsgfd37kg8l819JKIFjEEwAGxw+3qfewdjl5k1JvYOjOnx85tlSQFJSWSvy1j22dwPIFIIJYIPJZ+AEx8fU+41Pqvcbn1RwfCxmOyAlV9aKXBw31bxvSM37hnRxPN5InCF5/yArXQPsQDABbDDlDBzT1KUzv9WlM7+Vpi3h4qwcpK3xQ5KkoCk9fTqop08HFXeG8B3/SFE1OArBBLDBlLNyouCsHNjm1NPW2xbPkW58d+b6AmQAwQSwAWflIGsuJLH1fM1nMtcPIEMIJoBNwmflLKopm3K9tqacrcKwz9BZ622bN2euH0CGJHWIH4D41q2q081L36Saz4aef/1/NXNWDux1ccBau8VNbBGGIzFiAtiMs3KQWRbvpxW3ZrYbQIYQTIAMWLBggRYsWJDrbqAQLX9j5OGCCkMLKmIElUntACdhKgewWfmcCh34+a917sKonj0zqtXLKxg1gX2WvUGa41Gl+vTitrnR28zxhNoBDkQwAWzUebRXbQd7plSBrasp184N9Sx+hT1cRdKG+6T974vdZsN91C6BYzGVA9gkfIjf9NL0ZwZHddfDPk4Yhn3qb5fu+KZUvXjq9eprQ9frb89NvwAbMGIC2GDyIX7B8TGda98pSbpmU5tcJWUyFDrEb219LdM6sMXF5Wv11gP/Rxot1w//4S81Z8F10tKbGSmB4xFMABtMOcTPNDV28mjksTT1EL+bVszPTSdRUILBoH7845+EHte/Q6qszG2HAJswlQPYwOrhfBziBwDxEUwAG1g9nI9D/AAgPoIJYAMO8QMAexBMABtwiB8A2INgAtiEQ/wAIH3sygFsFD7Eb9HuCpmmySF+yKiKiopcdwGwHcEEsFn13CpdHBnOdTdQ4CorKzU8zH2GwsNUDgAAyBsEEwAAkDcIJoDNRkdHtX79eq1fv16joxRUQ2Zwn6FQscYEsNnExIQeffTRyGMgE7jPUKgYMQEAAHmDYAIAAPIGwQQAAOQNggkAAMgbBBMAAJA3HLUrxzRNSVIgEMhxT4DYJlfjDAQC7JhARnCfwUnCv7fDv8fjMUwrrfLE7373Oy1ZsiTX3QAAACk4efKkXvayl8Vt46hgEgwGdfr0ac2dO1eGwaFogUBAS5Ys0cmTJ1VdXZ3r7hQsfs7Zwc85e/hZZwc/56tM09SFCxe0ePFiuVzxV5E4airH5XIlTFqzUXV19ay/6bOBn3N28HPOHn7W2cHPOaSmpsZSOxa/AgCAvEEwAQAAeYNg4mBlZWXauXOnysrKct2VgsbPOTv4OWcPP+vs4OecGkctfgUAAIWNERMAAJA3CCYAACBvEEwAAEDecFQdE1jT2tqqHTt2yO1257orBcfn86mrq0uS1N3drX379vFztoHf71dHR4e8Xq/8fr+2bNnCzzUDuH+zj/fj5LH4tcD4fD41Njaqv7+ffwgZsGfPHm3fvj3y+Dvf+Y6OHDmS4145X2NjY+Tn6Pf71draqvb29hz3qvBw/2YX78epYSqnwPj9fnm93lx3oyD5fD7t2rUr8nzjxo3y+Xzy+/057JXzTf/5eb3eyF/1sA/3b/bxfpwagkkB6ejo0MaNG3PdjYLV0NCgffv2RZ4PDAxIkjweT456VBi6urpm/Aw9Ho98Pl+OelSYuH+zi/fj1BFMCsTAwABDhVkw+Y3mO9/5jlpaWvi5pyn8C3K6vr6+7HZkFuD+zQ7ej9NDMCkQ+/fvV0tLS667MWsMDAyoo6ODdRAZFCuwIH3cv5nF+3F62JWTxx588EEdO3Ys5utr165VS0uLurq6dMcdd2SxZ4XF6s95stbWVh06dIi/imzgdrtnjI709fXxs80g7t/M4f04fezKKQBdXV1TFrBt3bpV27dv17ve9S41NDTksGeFac+ePdq4caO8Xm/kr3re4FPn9/u1adOmKbtD5s2bp+PHj/NzzQDu38zi/Th9BJMCZBiGjh07xmrwDOjo6JDb7VZLS4sGBga0f/9+bdmyJdfdcrzp24W3bt2qQ4cO5bhXhYf7N/t4P04ewaSADAwM6MEHH1Rra6u2bNmirVu3ktBt5Pf7tWLFiinX3G63+vv7c9SjwuH3+7V37141Nzeru7ubglQZwP2bXbwfp45gAgAA8ga7cgAAQN4gmAAAgLxBMAEAAHmDYAIAAPIGwQQAAOQNggkAAMgbBBMAAJA3CCYAACBvEEwAAEDeIJgAAIC8QTABAAB54/8DCR+AplyGPgMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "scatter(OMEGA,AE)\n",
    "scatter(OMEGA,AE_Cardoso)\n",
    "axvline(x=-omegaB-1, color=\"k\", linestyle = \"--\")\n",
    "axvline(x=omegaB+1, color=\"k\", linestyle = \"--\")\n",
    "#ylim(0,2.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "1c81965f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "511f3ae3",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd00fec0",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6de7b270",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03cb99a5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd0561b3",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44e95c12",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47b6d205",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fac586ac",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01918047",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3d960aa",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c080b68",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ca9188d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7dc02d0b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa472711",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "402f6158",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAF8CAYAAADRg4HlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABHGklEQVR4nO3de3hb9YEn/K/kixwnto7lEAjgFMvEhVJaIjsEttPS1nZTprPdncTGnWm5zDCxmWk72+2AhcnMdum+qZGgs/vsvtO3UkiXlrYzjgUzZdphQKKFtlNw7ShhaLkk8QlNQgxJLB85N1913j8UnfhY5+hmH12/n+fJg/U750g/IUv6+nc1ybIsg4iIiKiAmHNdASIiIqJ0McAQERFRwWGAISIiooLDAENEREQFhwGGiIiICg4DDBERERUcBhgiIiIqOAwwREREVHDKc12BfBWJRHDixAnU1NTAZDLlujpEREQFQ5ZlnDlzBldeeSXMZmPaShhgdJw4cQINDQ25rgYREVHBOnbsGK6++mpD7psBRkdNTQ2A6P/82traHNeGiIiocExNTaGhoUH5LjUCA4yOWLdRbW0tAwwREVEGjByCwUG8REREVHAYYIiIiKjgMMAQERFRweEYGCIiStnCwgLm5uZyXQ3KsYqKCpSVleW0DgwwRESUlCzLePfddyFJUq6rQnlCEARcccUVOVsrjQGGiIiSioWXdevWobq6mgt8ljBZlnH+/HmcPHkSALB+/fqc1IMBhoiIElpYWFDCS319fa6rQ3lg1apVAICTJ09i3bp1OelO4iBeIiJKKDbmpbq6Osc1oXwS+33I1ZgoBhgiysj06f2YOvI0pk/vz3VVKEvYbXQJxwLl/veBXUhElLZT+x+B9IZHuS1c34vLNj2YwxoRZdfAwABcLlfWH1cURfh8PtjtdoiiiJ6eHgiCoHt+IBCAKIqw2+0AgPb29izV1HgMMESUlunT+1XhBQCkNzyoadiKqrWbclQrotLQ1dWFffv2AYiGmR07dmBoaEjz3EAggKGhIXg8HoiiiI6ODoyNjWWzuoZigCGitMyeOaJbzgBDpcDn86G7uzvrjyuKouq23W5HIBDQPb+3t1cJO3a7HX6/39D6ZRvHwBBRWiprGtMqJyo2fr8fDocj648bCARgs9lUZTabDcFgMO5cURQRCoUgCAKCwSAkSVK6kYoFAwwRpaVq7SYI1/eqyoTr72PrC6VseBh48snof43mdrs1yzs6OjS/+JORJAlNTU2aj7O4NcTpdK74QF+9+wuFQnFlwWAQNptNGS/j9Xrh8/lWtD65xi4kIkrbZZseRE3DVsyeOYLKmkaGF0qZ0wkszhR9fYBRY2FFUdQdtNrR0aH5xZ+M1+tFT09PXLnH41G6a2LnaQ3ylSQJAwMDCR+jvr4efX19KddJK9iEQiHl+QuCgJ6eHtTV1UGW5ZTvN98xwBBRRqrWbmJwobQMD6vDCxC9vW0bsGXLyj9eMBhEZ2cngOiXvNfrVYKBw+FAa2ur6ty2tja4XC7YbDaEQiGMjY3FhZCxsbG4WT+xABErDwQCusFJEISMZy8JghAXumLdREvZ7XYIgqAci/03GAzmpPvLCAwwRESUFQcP6pcbEWAWf9mPjo6qvriXts7Y7XbY7XZV60pdXZ0qbASDQXR0dMQ9ztLA4vf7Nc8DltcC097eDo/HE1e+OIgtfj7FjgGGiIiyork5vfLl8vv9aG9vhyRJcLlcSquL1+uNa4UYHR1VQkistWb37t2qcwYHBzVbT/x+v2pcjM/ng9/vh8/nU1qAYpbTArM0lIiiiNbWVlXriiAIShhrbW2FJEkQBEFZC6ZYWl8ABhgiSsP48DBCBw/C1tyM9Ub8yUxFbcuW6JiXxd1ITqcxrS8A0N3djZaWFrS3t2NoaAg7duxAY2Mj+vv747p4/H4/6uvrEQgEIEkSJiYmNMe6aBkdHQUA5dre3l4EAgHNlpHlGhoagtPpxObNmzEyMqJaA2ZgYACbN29WWm9i57a0tGDfvn1FN43aJBfTiJ4VNDU1BavVinA4jNra2lxXhyjnXnI6MbLom2dzXx9uy8FKpJR909PTOHLkCBobG1FVVbXs+xsejnYbNTcbF17SFfuSj4nNXoqFgdhsnqUtGJIkoa2tTXVtqUj0e5GN71BOoyaipMaHh1XhBQBG3G6MZ2MeLBWdLVuAO+/Mn/CiZWRkRBVW9NZ+MaqlhZJjFxIRJRXSGX0ZOniQXUlU8JxOJwAo66SMjIygo6NDNSZGa+0XIDoOpaurKzsVJRV2IelgFxLRJePDw/jBLbfElX/+lVcYYErASnchUXFgFxIR5b2DTz8dV3az08nwQkQ5wwBDRAlpjX8BgLMnTuSgNkREUQwwRJSQ3viX1598koN4iShnGGCIKCFbglXG9MINEZHRGGCIKKH1W7bg+i98QfNYonBDRGQkBhgiSmrNlVfGlXEQLxHlEgMMESWkN4h34x/+YQ5qQ5QfYjtQU+4URICJbaylt7unHqfTCa/XC6/XqyxQRETpEZ99Nq1yolKQbEdpo4iiCLfbDZ/PB7fbnXKQcjqdqnN9Ph8kSdK8PhgMIhgMKo8X+znf5P1KvMFgEKOjo5AkSbU1eiKxvSleeOEFCIKAYDCIlpYWcM0+ovSdOX48rXIiMk5XV5ey75IoitixY4dqQ0ctwWAQbrcb/f39qvtZyuVyoa+vDx6PB16vFwCUjTDzUd4HGIfDAYfDkVYLitPpRHd3t7LFuMPhKLpdOImIKDd8Ph+6u7uz/riiKKpu2+12BAKBlK6z2+3KbUmSMDQ0hM7OTqXM7XYrG1e2tLRgcnISAJTv0XxUEF1I6fJ6vejs7IQoisqLu3TrdCJKTc3VV6dVTlTs9DZ2NFogEIDNZlOV2Wy2hF08Pp9PFVRiFpdpnSMIQl6HF6AIA0wsoQaDQUiSBLvdjt7e3qQpdWZmBlNTU6p/RATUNjSkVU6UzPjwMH6bpYUQ3RoD0AGgo6Mjo7Edehs7ut1u1ffM0jEnK0Hv/vSGV0iSpBlCFpfFhmcsbaHx+Xzw+XxwOp1xLT/5Iu+7kNIV+x8tCIKSkF0uFxobG5UmMS0DAwN4+OGHs1JHokJirqxMq5wokZecTtWsts19fbjN5TLksURR1G197+joSHlc5WJerxc9PT1x5R6PRxmbEjvPpfG8JElKOgC4vr5e6c5JhV6w2bt3r2ZdF3M6nXH17OnpUUKO3W5HR0cHxsbGUq5PthRdgIlpbW1VfhYEAZIkIRAI6P4y9/f346tf/apye2pqCg38C5MIB771Lc3y0JtvZrkmVOi0puSPuN1o3rbNkDWFgsGg0jUSm80aCwYOh0P1PREMBtHW1gaXywWbzYZQKISxsbG4L/exsbG4Vo1YgIiVJ/quEQRBM9ikQhCEuNAVCoU0W1kCgQDuuOOOhPcX+15cer0oikoDgN1uhyiKceNo8kHRBRi9/8GCICRsBrNYLLBYLEZVi6ggjQ8PY/yVVzSP/fZ738NHd+3Kco2okOltPRE6eNCQALP4y350dFQ1bmVp64zdbofdble1WNTV1anCRjAY1FzOY2lg8fv9ust+LKcFpr29HR6PJ658cRBbbO/evcrPoihiYGAA3d3dyv+H0dHRuPASC3JLeyyWjr3JB0UZYGKJcfEvqyRJui8yEWlLtNfR2ePHMT48zNV4KWV6W08YtSWF3+9He3s7JEmCy+VSWl28Xm/cINzR0VElhMRaa3bv3q06Z3BwULP1xO/3q8bF+Hw++P1+3cGxmbbALP0DXRRFtLa2KiEkGAxCEATY7fa4FqDe3l709vaq7iMYDMYFE7vdrqpfIBBAZ2dnXg7oLZhBvHp9lbFFfRZzuVwYHBxUbvt8PrS3t+dk1DhRIYvMziY8zs0cKR3rt2zB5iUtC0ZuSdHd3Y2WlhYMDAxgaGgIoiiisbERQPzMVL/fj/r6egQCAQQCAUxMTKQ8e3V0dBRjY2MIBALw+XzKxBEjulyGhobgdDrh8/ng8XhUa7QMDAzELTkiSZLyHelyueIGLi+toyAIaG1thdvthtfrxcjISN6uA2OS83x1N1EU4fP5MDg4iGAwiL6+PmzevFlJtbGBUksHGHm9XqVfcmJiIu3EOzU1BavVinA4jNra2hV5LkSF5rdPPoln77pL9/jnX3mFLTAlYHp6GkeOHEFjYyOqqqqWfX/jw8MIHTwIW3Nz3vz+tLS0qAbhxr70Y105Pp8Pdrs97g/h2MKpi68tFYl+L7LxHZr3XUh2ux19fX26I7J7eno0R1knG3lNRMklaoG58iMfyZsvHyos67dsyfvfnZGREfT29iq3/X6/5viTQCDA4Qk5kvcBhohy552XX9Y99olvfjOLNSEyjtPpBACl+2VkZAQdHR2qMTFaa78A0V4CrWX5yXh534WUK+xCIgK+u2kTTh04EFduqavDlzNYQ4MK00p3IVFxyHUXUsEM4iWi/DEzOYmXLv7VSkSUCwwwRKSrYs0a3WOHn/pfOBp4FNOn92exRkREUQwwRKSreu1azfIPdNfhY1+/EjMnv4Vjz2/Dqf2PZLlmRFTqGGCISNcJjc326poqsfGzVlWZ9IaHLTFElFUMMESk6bU9e3B+fDyufPX6Cs3zZ88cMbpKREQKBhgi0vTb731Ps/zc+JxmeWVNo5HVISJSYYAhorRMjs3i0DNhVZlw/X2oWrspRzUiyr7YSu+UO1zIjog0XfV7v4fjP/+55rHXBydRd/123HDX7aisaWR4oZIzMDCQ8aaMyxHbXie2aXFPT4/uRouBQABANGyNjIyodqKO7YnkcDggiiIkSSq4/QIZYIhIk7lCe6wLAGx9/HHceO+9WawNEQFAV1eXsu+SKIrYsWOH7maLXV1deOGFF9De3o5QKISuri5l30CPxwOv1wsgurFlvm7YmAgDDBFpOnP8uGb52g99iOGFSprP50N3d3fWH1cURdVtu92utLJoGRoaUrWqLG6paWlpweTkZFx5IeEYGCLSFHrrLc3ySm6tQSXO7/fnpLslEAjAZrOpymw2m9IdtFRsLycgGmYWb04JRINLoYYXgAGGiHRcOHUqYfn06f2YOvI013+htGXzd8ftdmuWd3R06H7xJ6K3saPb7Va1hjidzhUf6Kt3f6EE+5IFg0E4nU50dHSgp6dHdV8+nw8+nw9OpzOudacQsAuJiDSFj2iv61LX3IxT+x+B9IZHKROu78Vlmx7MVtWogGXzd0cURVUrxGIdHR0Jv/j1eL1eVRCI8Xg8ytiU2Hlag3wlScLAwEDCx6ivr0dfX1/KdUoUlBwOB+x2O5xOJ3w+Hzo7OwFANfjXbrejo6NDGR9TKBhgiCjOL3buRGR2VvPYzX91B6Q3/kZVJr3hQU3DVs5GooSmT+9XhRfA2N+dYDCofGFLkgSv16sEA4fDgdbWVtW5bW1tcLlcsNlsCIVCGBsbiwshY2Njcd0usQARKw8EArrBSRCEjGcvCYIQF7pCoVDSbiBBENDV1YWOjg5MTk5CEASIoqh0g8VmNImiCLvdnlHdcoEBhojinHj5Zc3ymg0bULuhGhfeiT82e+YIAwwlpLdas1G/O4u/7EdHR1XjVpa2ztjtdtjtdlXrSl1dnSpsBINBdHR0xD3O0sDi9/s1zwOW1wLT3t4Oj8cTV744iC2uU1dXlzJQNxZMYl1FbW1tyrGYpeNr8h0DDBHFqd2wQbN8wyc+obviLlfipWSy/bvj9/vR3t4OSZLgcrmUVhev1xs3CHd0dFQJIbHWmt27d6vOGRwc1Gw98fv9qnExPp8Pfr9f1WUTs5wWmKWtI6IoorW1VWmBCQaDEAQBdrsdNptNFapixxwOh/L/IyYQCKCzs7PgBvQywBBRnAs6YwMuhEKoWrsJwvW9S8YxcCVeSi7bvzvd3d1oaWlR1jnZsWMHGhsb0d/fH9fF4/f7UV9fj0AgAEmSMDExoTnWRcvo6CgAKNf29vYiEAhotows19DQEJxOJzZv3oyRkRHV+i0DAwPYvHkz+vr64HA40N3draz14vf7lTE6giCgtbUVbrcbgiBgbGysINeBMcmyLOe6EvloamoKVqsV4XAYtZw2SiVmz3XXYVJjGnXd+9+Pe998E9On9+PciRcBAKuv/DjDS5Gbnp7GkSNH0NjYiKqqquXf3+n9mD1zJK9WcW5paVENwo3NXop15cRWv13aciNJEtra2lTXlopEvxfZ+A5lCwwRxSnT+ZIqq6qKm0USWZjJmy8hKgxVazfl/e/MyMiIat0Uv9+vOf7EqJYWSo4BhojiWKxWzXLbxqqsziIhygan0wkg2soCRMNLR0eHakyM1tovQHQcSldXV3YqSioMMEQUZ2bJ7IQY2/vrAbwXV84ZSFTIlg6q1Rp4q7cuSzrrtdDK4kq8RKQyPjyM06+9pnnsmo4/0iznDCQiyjYGGCJSeXnXLs1ya1MTrr7tC4B5lfqAeRVbX4go6xhgiEjl7Dsaq9QhOi7m9IFHgcgF9YHIhWg5EVEWMcAQkUqiGUgXTmtvfqdXTsUlEonkugqUR3L9+8BBvESkUm6x6JavWuvA9MlX4o6tWuvQuIKKRWVlJcxmM06cOIHLLrsMlZWVMJlMua4W5Ygsy5idncWpU6dgNptRWVmZk3owwBCRyvmTJzXLazdswNqbHsDk69+KO7b2pgeMrhblkNlsRmNjI8bHx3HixIlcV4fyRHV1NTZs2ACzOTedOQwwRKQYHx7GxG9/q3nMVF6uO9bl9IFHGWKKXGVlJTZs2ID5+XksLCzkujqUY2VlZSgvL89pSxwDDBEpQgcP6h6bnpjQHesydeQfGWBKgMlkQkVFBSoqKnJdFSIO4iWiSyKzs7rHmv7gD1BRfZXmsYUL45g+vd+oahERxWGAISLFoR/9SLO8UhBw47334uyJn+peG9vckYgoGxhgiEihtwaMYLcjfHgQ8qz2FgMAMH/+XaOqRUQUhwGGiBSJ1oA5984LCa8tr77CiCoREWligCEixZljxzTLyy0WzJ7VPhaz+sqPG1AjIiJtnIVERACiU6jP6gWYVatgMutPl6y55g+5HxIRZRVbYIgIQOIp1GWVlbBYr9M+Vt2AK/7D3xpVLSIiTQwwRAQAsDU36x6rqq/HuXd/qXmsvPpyo6pERKSrIAKMJEnwer3o6OjI6PpMryMqJYefeUb32Afv/Dgi0+9pHps/P25UlYiIdOX9GJhgMIjR0VFIkoRQKJT29T6fD4FAwICaERWXI//6r5rlqy6/HKvqTuGC9gxrWKz6LTdEREbJ+wDjcDjgcDjg8/nSvjbT0ENUisqrqzXL6zZuxLl3f6V7Xc01f2BUlYiIdBVEF1Km9u7dizvuuCPX1SAqCLb3v1+3XF64oHtdZU2jUVUiItKV9y0wmQoEAmhvb0/5/JmZGczMzCi3p6amjKgWUd469uKLmuXnT56ExXodZifjd6muqG3m9GkiyomibYGRJAl2uz3l8wcGBmC1WpV/DQ0NBtaOKL+MDw8jPDameWwmHMb5UyOaxyrW8H1CRLlRlAHG6/Wis7MzrWv6+/sRDoeVf8d0FvQiKkbis8/qHruiZR0Wzh3VPDYjvWFUlYiIEiq6LqRgMIjW1ta0r7NYLLBYLAbUiCj/nT1+XPdYY8f7MDc5qnnMZK40qkpERAkVXYAJhUIIBoPK1Omxi83ibrcbdrs97ZYZolJw/uRJzXKhqQk1V12FkM4m1NXrbjawVkRE+gomwOhNhxZFET6fD319fQCA9vZ21eDdYDAIr9erHCeieFM6XaYVNTWQI/O611mv/ZxRVSIiSijvx8CIogi32w2Px4NgMAin06laEyYQCMDj8Whe6/P5MDAwAABwOp1c0I5Ix4XTpzXLK2trce7Ei5rHzJZ6zkAiopwxybIs57oS+WhqagpWqxXhcBi1tbW5rg6RYcaHh/GDW27RPHbD3Xfj2tuHgYXpuGNmy1o0bdeenUREpS0b36F53wJDRMZKtAu1ueyEZngBAHPFGqOqRESUFAMMUYmLzM7qHhMa9ce/VNY2GVEdIqKUMMAQlbhDP/qR7rHLPnyN7rEq2w0G1IaIKDUMMEQl7uw72ttMr7r8cpRbzuhet/rKjxtUIyKi5BhgiEpcWVWVZrnQ1ISZ0G+0LzJXcQYSEeUUAwxRiTt54IBm+SrhLBDRHsBrKuOq1USUWwwwRCXstT17sHD+vOaxdR/SbpkBgIrVVxlVJSKilDDAEJWwsR//WPfYFa11useE5ruMqA4RUcoYYIhKmN4KvGuuugry/O90riqD9dpu4ypFRJQCBhiiEjZ79qxm+WUfqtNdwM7EBeyIKA8wwBCVML0p1ELjnO41ZRb9riUiomxhgCEqUePDw7hw6pTmMVuzfitL9bqbjaoSEVHKGGCISpT47LO6xyy1M7rHrNd+zojqEBGlhQGGqESdOX5c91hZhfbYmAvTlXh1jAvYEVHuMcAQlajQW29pln/8G+th1vlkOD1VjwQNN0REWcMAQ1SitMa/1DVVwvo+/VV2//nf/tDIKhERpYwBhqhEXQiF4srWfXiV7vmRCPCU72NYd/xJjA8PG1k1IqKkGGCIStD48DCmNVpg6j+gv33AyXdX4y9xC87vuQs/uOUWvOR0GllFIqKEGGCIStDLu3Zpltc363cfTb85rro94nazJYaIcoYBhqgETR48GFf2e//tCpjL9D8Sjr4YPzMppHE/RETZwABDVIIiCwuq23VNlah/v373UUiqxeTYbFy5rbl5xetGRJQKBhiiEjQ9MaG6/aE/rU94/qN7H8JP0acqu9npxPotW1a8bkREqSjPdQWIKLvGh4cxMzmp3K5rqoRwjf7Yl5m5cjz1UjeAbvwG2/D/fu0gNt/ezPBCRDnFAENUYpYO4N3052sTnv9Pv9im/HwMWxDEFnyW2YWIcoxdSEQl5tSrryo/1zVVomZ9ZcLzh17k3kdElH8YYIhKzPSiBexuca5LeO6L+2/Dvy/Z+6ihwZBqERGlhQGGqISMDw9j7mx0OvSG21ajcrV+L/KZ89XoeeyJuPJjx4yqHRFR6hhgiErIq7t3Kz8nmnkky0DLjt9mo0pERBlhgCEqIcd/8QsAwHWdVpSV67/9A/vadI/dfvuKV4uIKG0MMEQl5Nx4dDuAjf/JmvA8zzNf1Cx3OgHOniaifMAAQ1QixoeHMXfmDG7aYYPZrP/WD761KW7gLgA89BDwyCNG1pCIKHUMMEQlIrb+S8PHanTPmZ834XNffzquvKIC0Nn/kYgoJxhgiErEO7/8JdoeuxJms0n3nL/5zoBm+ac/bVStiIgywwBDVALGh4dRbTuHNQkWrZuYqru4ZUC8nTuNqhkRUWYYYIhKgPjss7j2s/oDd2UAt/55UPNYVRUH7hJR/mGAISoB//7447jCsUr3+C/f2K577LrrjKgREdHyMMAQFbnx4WFs/P0LCWce9Qw8pnvsS18yolZERMvDAENU5MRnn0XDx9boHh+f246FBe1jZjNw770GVYyIaBn0N0LJI5IkYe/evRgaGoLf70/pGrfbDQAYGxsDAHg8HsPqR5TPTAt7dVtfZBn4m+/ot75s3GhUrYiIlifvA0wwGMTo6CgkSUJo0S66iTidTrhcLuV2b28vOjo6Ug4/RMWk/v3noNfYujD/IVxcnFfTdv2hMUREOZX3XUgOhwM9PT2w2+0pnS9JEoLBICRJUsp6e3sRCAQgiqJBtSTKT2989z/pt75EZFx/949w+LD2tSYTF68jovyV9wEmE6Ojo6qwEgs/i0MNUSkwl/277rFTb16OnTuj3UhaGhoMqhQR0QrI+y6kdAmCgMnJSVVZIBAAgIStODMzM5iZmVFuT01NGVNBoiw5feBR6E08kiMyaht78J2/1r/+wx82pl5ERCsh5RaY7u5uvP322wZWxTgDAwPweDwQBCHhOVarVfnXwD8/qcBNvv4t3WPj+y/gxnvvxbvv6l/P1XeJKJ+lHGCGhoYyGkOydetWdHd34/HHH8eBAwfSvn65nE4nuru70dPTk/C8/v5+hMNh5d+xY8eyVEOilff2jz+V8Pj5iRsTBpTqaq6+S0T5LeUuJEEQ8MlPflJVduDAAdx0000Jr5uYmMBzzz2neWz37t0wmUy44447UFtbm2pVUubz+dDU1JQ0vACAxWKBxWJZ8ToQZdv06f2Ymzqkezz8uxnc9Oe70PlJ3VPwla+sfL2IiFZSyi0wNpstrmzHjh1Jr2tvb9c9tmPHDnR1daGzsxN/8Rd/kWpVUhIb9xILL5IkcRYSlYSJ3/yd7jFZlvHL/yeMo9iC8+f174Ozj4go36UcYEwmU1xZXV1d0uvq6+uVn1944QU89thjqq4kq9WK559/Hv/wD/+Q8H701oARRVFZtC4mGAwiGAzC4XBAFEWIogiv16sZwoiKzfkTL+geG993Hi1f+Qp279a/nsO/iKgQpBxgxsbG8NBDD6nCx9LZPloWB5+2tja0tbWhq6sLjz/+uGqmzx133KF5fSygeDweBINBOJ1O+Hw+5XggEFCtsitJEtra2uB0OtHU1KT8czqdCQfxEhWDZGNfDj8Txkd37cLPf65/zqcS3wURUV5IawzM4cOH8clPfhImkwnt7e0pdclMTEyobm/atAnbt2/Hn/3Zn6nKm5qaNK+32+3o6+tDX1+f5vGenh7VGBetadREpSDZ2JeJty4A5Y0AgEP6pyGFnmEiopxLawzM3r17EQqFMDo6itbWVjQ2NqK+vh7f/OY3da/zeDzYunUrHnvsMbz66qsAgLVr1y6/5kSkcvynd+seiyzI+OXX30Pz9u245x79+6ip4ewjIioMKQcYh8Oh/NzY2IgHHngAo6OjGBkZwaFDh3SDjN1uR2NjI7797W9j06ZNKCsrg8fjQX9/P372s58p52mNsSGi1IQPD0KeP6N7/NAzEgDgo7t24ckn9e/ny19e4YoRERnEJMt6C4mrPfXUU2hqatKdNi1JEh555BFce+21qu6hp556Ctsv7ggXDofh9/sxOjqKQCCAYDAIk8kEh8MBSZJwKFG7dpZNTU3BarUiHA4bMsWbaCUd+qEdgPZbObIg45/v+h1uuPtunPjoE1jSe6uS2qcBEVFi2fgOTTnApCIcDqO1tTWtIBIIBOD3+/H444/HjZfJJQYYKhSnDzyacNXdt/5xEm/6wrhflrFpE6C3nuTGjcDBg8bUkYhKSza+Q1d0L6S2tra0r2lvb0d7ezvHxRBlKFF4WZiN4E1fGFseeggA8Npr+vfjdK50zYiIjLOiAaa3txetra0ZXdvZ2bmSVSEqCYf+/jrdY7Is48d/chQwmfDRXbuwZw+wsKB9bkUFcO+9BlWSiMgAKQ/iTcWOHTuwadOmhOf89Kc/1SxvbGxcyaoQFb3TBx4F5Bnd40dfig7q3dLfDwDo7dW/r1tvXdGqEREZbkUDTCoWLzpHRJlL1HUkR2Qc2B1KqfUFAO66y4AKEhEZaEW7kIDowJ1QKARRFJX9hyYmJpQl/YPB4Eo/JFHJOfTDxC2WBx4/DeBS68sDD+ifu3o1u4+IqPBkFGCefvppiKKIsbExiKKoCixLLZ3kxPVeiJbn0OCNCY9Ph+dx9KVzSuvL8DCQaHHqs2dXuIJERFmQUYB5/vnn4fV6YTKZlIASW7DObrfDZrOhqakJg4ODcLlcSpnVas14kC8RAcd/djewoJ845IiM5/7iOIBLrS+JJgdefvmKVo+IKGsyCjDf/va30dPTgwcffBD33Xcftm3bpnneyMhI3NTq9vb2TB6SqORNn96PC+MJdmHEpa4jmM3K2Jdz5/TPZ9cRERWqjAfxOhwOPP/885BlGf39/aqdpWO0uoseeeSRTB+SqKQde74r4fGz785Fu44AbPV6AQBf/GLi+9y1a0WqRkSUdcsexLt9+3Zs374djz76KEwmE+6///6VqBcRLXJ4aBMA/WlE87MRvPBX7wAAqtauxY333os9e4AZ/VnWuFt/70ciory3YtOoH3jgAWzfvh333Xef7lovRJS+4z+7G/KclPCc1564tA3Hl06dAoCEex4BwBNPLLNiREQ5tKLrwMR2nZZlGQ8++OBK3jVRSQofHkw67kXVdfT44wCAK69MfL8XTyMiKlgrupnjUoXcrcTNHCnX3v3Vf8WZt/8p4TkzZ+fxr73RWUdVa9fiS6dOYc+exK0vq1YB58+vYEWJiJYouN2otYTDYQwMDODmm2/Gtm3bsHXrVjz33HNGPuSKYIChXHrnZ3fjfJKWl/nZCH7yJ0eV2/dffCsnW2rplVeALVuWXUUiIl3Z+A41fCsBq9WKRx55BFarFddeey0CgYDRD0lU0N7+8dak4SWyIKvCy+dfeQUAYLcnvu+rrmJ4IaLikLW9kNra2vD8889n6+GICtKhv2/G3NTBhOfIERn/fNfvlNsfuPNOrN+yBXfeCRw5kvj+n3pqJWpJRJR7Wd3M0W63xy1sR0TRReoO/bARkOcSnidHZDxz56XwUv/BD+L3v/c9DA8D3/9+4se4/Xa2vhBR8TB8DMxS4XAYVqs1mw+ZEY6BoWw5+tw2zEzsT3re0vBirqzEVy8u9JJs3IvVCmhsVUZEZIhsfIeu+G7UyRRCeCHKlkN/35y01QWIDy8AUg4vAMMLERWfrHYhEVFU+PBgSl1GALAwH4kLL6nOOAKis46IiIpN1ltgiErd4b0fgjx/JqVz52fm8ZM/PX6poKwM98/PA0gtvHDcCxEVKwYYoiyYPr0fJ4PfwMzp0ZSvmXjrAn759feU2+mMeQGA9euBf/mXtKtKRFQQGGCIDPb2j7cmnRq9WGQhgl8+/C4mx2aVssraWvxlOIydO4FvfCP5fVitwIkTmdSWiKgwMMAQGeSdn/0pzo//LK1r4rqMANQ2NqJHFGGxALOzOhcuUlnJQbtEVPw4iJdoBU2f3o+3f3I7Dv2wMe3wMvHWhbjwsvXxx/Grj4kwmVIPLxd7mYiIihpbYIhWwLu/uv/ixosLaV+7MLeAf/sf76m6jMwVFbD+f7P4UIJNGZeqrwdOn0774YmIChIDDFEGpk/vx8Rv/g7n33sZWMhsa+dIRMahH0l40xdWlc9bG7HzjIhIGuHl058Gnn02o2oQERUkBhiiFIUPD2Li9W9j4exRAJGM7ycSieC9/Rfw6789pS4HMITHMRq+N6374+7SRFSKGGCIdChTn0O/BSIXVuQ+w7+bxosPvasqiwAYwd3w4Ym07y+7G4EQEeUPBhiii04feBThsUFE5s4CkTksp5VlsUgkgtBbM3j97ydV41yWE1xuuAH4zW9WpHpERAWJAYZKjtKyMvkmEJkFIKe0pH+6IpEIjv38LA7sDqnLkXlwMZmAl19mlxEREQMMFa24oGIyrWjLipZIJIK5sxG8/cKZuMG5ywkuAPDQQ8CuXcuvIxFRMWCAoYKlDKq98F50MIipLHpAXrjYsrIkqBg4XmRhLoLj/6bd2nIeNryM+/A8Mksfd98NPPHE8utIRFRMGGAor4QPD2Ly4Pcwf+4dyJF5/RMXzsPQRJJEZD6C2XMRTB6awaFnwnFjW2axGq+iM+PWFpMJ6O9niwsRkR4GGFpRsfVRpicOQI7Mw2QuV4KI1s9yZB6RuVnMXZhDhUWGuTyFXQpzIBKJIDIr48LEAg7/JIyjL51THV8AMI9VOIRP4Lv4ScaP09AAfO1rwL3pzaQmIio5BRFgJEnC3r17MTQ0BL/fn9I1brcbgiAo1/f19RlYw/wXPjyI8JGnMX/2GOTIHEzmCpjKLJAXZpTbADR/nj1zDnMXpoBIbJVZGbJsUu2ILMtARZUJ5VXq3SnkFH42mwHLagDIn/AyP7cALMiYOw8cfSl+PMsCgAgqcAZXwI+vYRSZJ47Vq4HOTnYTERGlI+8DTDAYxOjoKCRJQigUSn4BouEFAHp6egAAgUAAvb298Hg8aT/+sRfuwhrLHCLz52Eur4apvBry/HlE5s9jfnoa0xNnMBO+AFN5BGXR730szAELM0CZBUoZzCaYTWZE5AgWZmTVcVmOAHIEMJlhMiXenmrpueby6NAPeT6C+dkIIhEzZJgBRGBCBOYKM6oFE8oqMg8HZeVAWQ0AlGV8H/lqYT4CeSECwIT5aVmzSwgA5gFEUIlJXIMX0beswFJWBlx5JVtaiIiWwyTLhbEUls/nw8DAAPbt25f03Lq6Ohw5ckRpgQEAk8mEdJ7q1NQUrFYrgrs3YE0197wsdAvzEcjzkWjLjwxMhyJxXUERAPOogAkyZJgxizV4Gx/BT7ETx5D5vGWzGRAE4L77OKaFiDIzPAzs3g3s2wdMTEQ3d62sBKqqgOlp9W1Au2yxZMdTtXo1sGYNsO53e3B56Meomj+Ncvks5rGA/zn7GsLhMGprazN/gATyvgUmXaIoQpIkVXiJCQQCaG9vz36lyHDzswtYmI0GVJM5GhoW5oELJxdwJDClBJUFAAuoAGCCDBNMqMA8qvEa/nPGA24XM5mA6mp2CREZbc8ewO0Gji/awH1hIdrCWV4OzM+rb2eqvPzS/c3Pq2+vmx7GR2d3YV3kVZRFpgGYon8AmUwwyTJkmADIQOy2yYQIKgEAZswqZY23VcL+kTKsEi62nJdHW9flhQjm56Kt6mZE0HtdBGUfNMNUbsbCQrS134wIon9+mRFB6n9sz88A89NAeRVQbknt+KUyM8osJpgQwfyMDExPo6IKKLOYsDAjY35axrnpCP7nY5n/f09FUQYYLYIgQJIk3etmZmYwMzOj3J6amlrpqlGa5mcWsDCn32oWmQUkMdrlc3psFvOwQL7YzWVatCu0jDKcxwbswxcynsoMRD8IFzfiVVYCra3AY49xYTmiRO65B/infwLOnYs/ZjJF/+CQZSASUd+W5ejPSy36qNbUgGF8ErvQgF+jAqltAxJBORZQjnVNMj742XJYN5TDfPEbMtoqq+6GN8WVVF66r4UIIJtgLjdp/By5+HMZFuZkVNeXo9wS6543X/yHRbehcyxZuZGWrqUV34Sz6rxx623FFF2A0WOz2RKOoRkYGMDDDz+cxRoVv/npBcydl1XLsyz9efFteQGIzMuYPBrB689MY2LsUggxYx4mLEBGGSKLfm2nIWAfdmQUTJL9VWaxAJs2MaAQLTY8DPT2Am++CcwtWcDaZIr+i4UPkynaCpKJT2EnWvB9VEFSvf8BaP4c+1y49jZgY/tqVNWXwVRWAaBCuc/Y54zWZ4+pLNrSUF5VMl+LBa9kXqlkA4D7+/vx1a9+Vbk9NTWFhoYGo6uVdfMzC5g7JyMyL2NhPoKyChPMZWZE5mXAJCf8WQYudrvE/zUCRP8iWZg3QTq6gDeeuYBTYyYsoBxlmIcZ0enTsb9yFpfNYRWO4eaEY01MJqCm5lKzcEzZxQ8hy5IPydhfciZT9GdZjgaWDRuAvj4OniXSs2cP8PDDwIkTl95ri99Ly9WJe3ADfowKnIXp4qeKDLPq5wrMKG0KdU2V2PhZK+o2WlBuufh+1qlHWZUZ5eUcs1gqii7A2O12zXJJknSPAYDFYoHFEt8RePrgBcxZy1R9e+VVJpRZol/gsTJzVQVMluj/Tu2+RdPFN2m0z3Dx8XT6MJeeW1ZxsWtjPoKF+Uvl5kWzkMxlwAVJxlsBYOzngBlziJgqAZhhlmcRQRkAE8yYBeR5LO3HnTbV4a3K7fiVdRfmL64tp9cnvPh4+Tp1WTLz89EgEuuzjgWOL32JgYNopcWCyjvvJA8msRaVVLRiDz4ONwQcvdg6Yr7Y3RK5+IeL2nWdVmy4rQYV1dHlFUy4VB9zJQMJ6SvKACMIAkRRjAssmQzg/fH/ugJlZissOIsKnIMsAwsRYN68GufKGzBR5cCx1bfjVFV6fQyZjBBfWrb4NhD9C+nDHwZ27gQ2a1Tn1r9O88kTUVEYHga+8AVAFFemFQVQB5VYa6oZ85qLLWy4bTWu/YwVq+rLAJMc7V4CLv6BxYBCmSmYAKPXBSSKInw+n2qhuv7+fgQCAWUdGJ/Pp/ycrv9x+peGTQEjIjLCSgeWBgxjG3qxDm/CjLmEbcU37bDhCkc1zJXR5SvMZSaY2YqSksi8jDMnF+JmIeFi6z0QQbnFjLIK88Uxg+mtI7aY5nplSY7HyspXmVFeaVKta1a5uhxVtlqsWnsFTOXVmDu/AODoSvxv0ZX368DEAsrg4CCCwSD6+vqwefNmdHZ2AgC8Xi9cLhfGxsZU17ndbqUFZmRkBC6XK63Hja0DY+QcdiKilRILLYcPL/++WrEHHXgYNXgXJszp/qWralkxA+ZyE8xa04YKVGRexszZCGKdW0tXHzeZAFNFJSxrVsNcUXnpvAQrmy/+ObYauqm8CjUNt2PtTQ9k9fkZKRvfoXkfYHKFAYaICsFHPgL86lfLv59O3IMb8U+wIKy75rYSWNaWRcNKnnf/yLIJpvJq5bbJXK7ah21pGQCYy6tRvqYB1sZtsF7bnZN6F4NsfIcWTBcSERFF7dwJ/J//A5w5s7z7uRufQTMCKMesbpfQdZ1WXNNWg4rVeTJexVx1aQ601uHyKlTV34T6D34RVWs3ZbFilG0MMEREBeKee4Dvfnd59/Ep7MSt+DaqEUo4jmV962qUV2e5S8hUGe2XUZZpkAFTOSpq3oe65rvYIkIqDDBERHlueBi45Zbl3Ucn7sFmfDd5S8sas8GhxQSYKpQV70zl1VhzVRuu+A8GrztPRYcBhogoT+3ZA3z5y8CF1FbD13Q3PoMP4F90g8vNX70Ml29aZUxoMVVEZ8iYK2Gx3YB1jofYrUMrhgGGiCgPrV0b3XU4U924Bw6dFpe6pkp86E/qUfu+SpjN8atqZ64MKFuFmoatbFEhwzHAEBHlkeV0F5WVAX95xT248h3tgTIbbluNGz5vQ+Vq/UGwaTFVoGz1Vaj/wH0cn0JZxwBDRJQnbr8d+Nd/Tf+6G28E/u4/7sTwN74BvBN/fMNtq3HjPfUor1xuN5EJZWvex8BCeYEBhogoD2TSZXT33cATTwDfrKzE8Gtzmud82nM1LGuW81FvQoX1/bhiyzc4foXyCgMMEVGOlZerd1lP5qGHgF27oj8/ZtIew3LTDhsaPlaT8RgXy9pWDrqlvMYAQ0SUQzr5Q9OVV0Z3jwaAZ++5B7/VWBSmrqkSH/naFSjLYNG5sjXXsHuICgYDDBFRjqQaXszm6HYBWy7uMv9YWZnmLo1b7l+HKzZVx5UnU73+E7jqE99J+zqiXGKAISLKgVTDy+LuovHhYfxAZ4rS7buvRmV1eh/pNdds53RnKlgMMEREWZZqeFm81e5P7rwTb3z/+3Hn1DVV4vf++/q0xrrUfeAvimrnYypNDDBERFmUSXj5l7vu0gwv6XYZscWFigkDDBFRlpSluH7c4vAyPjyM1598Mu6crd+6GlXW1D7CTeU1uPaOf0/twYkKBAMMEVEWrFqlOe42zuLwAkBzzMtnvtOAcktqaWjdzY9wVhEVJQYYIiKD3XorMD2d/Lyl4UVrjZfPPLEB5RWpTZHe+MdHUjqPqBAZuWc6EVHJGx4GXnkl8TkmU2rh5Q+++76UwktFbTPDCxU9tsAQERkk1Y0Zl3Yt6bW8lJUnHwHc8KmnuXoulQS2wBARGcDpTC28LG2d+VuLJe6cz/zf1LqNNv7xEYYXKhkMMEREK2x4GHC7k5/ndF5aXRcAfnDrrYjMzqrO+cz/3ZDCLtImdhlRyWEXEhHRChgeBnbvBg4dAn7+8+TnL15hF4hOlx5f0hzzme+kEF5MFdj4RwczqDFRYWOAISJaJqcztRaXmMZGdXgB4qdLb/3W1Si3JAkvZgs2fu7N1B+YqIiwC4mIaBlS7S6Kqa4GRFFdtnTcyy0PrEthkToTwwuVNAYYIqJlOJhm7825c+rbv9i5UzXupa6pEpfflGx7ABM2/rGY5Byi4sYAQ0S0DM3NqZ+7dK0XABj+xjdUtz/68Pqk98PwQsQAQ0S0LFu2AFddlfw8rcXs/rcgqG7/x++9D6Ykuz1ythFRFAMMEdEyORyJj995p3q6NAC8tmcPZsNh5fZn9myAuYzhhShVDDBERMu0KIfEueUW4Hvfiy/3f/GLys/XdVpRXpX445jhhUiNAYaIaBnGh4excepJNGA47tiaNcDLL8df89qePYjMzCi3m/+zkPAxGj719HKrSVR0GGCIiDL0ktOJH9xyC647cBf+ErfgdjhVx//qr7Sve/7P/1z5+bPfTzzuxVK/idsDEGlggCEiysD48DBGliwA80m4VS0xt98ef91re/ZAnpsDkDy8AGXYsJWtL0RaGGCIiDIQ0lkA5jJEy5fucxTz4v33A0h1xtHh5VWSqIhxKwEiogzYdBaA+eLXmnHD7drhZXx4GDOShNt3X80ZR0TLxBYYIqIMrN+yBZv7+lRlNzud+NP/vkUzvADAz+6/H5/2XI3K6sR/O3LQLlFybIEhIsrQbS4XmrdtQ+jgQdiam7FeL7lctOne36G8MvHHLgftEqWGAYaIaBnWb9mSNLgAwFvfb0R5ZVnik0wVHLRLlCJ2IRERGWj69H4c+mEjzCl82m78ozR3hiQqYWyBISIyyNHn/hAzEwdSOpeDdonSUzABxu12Q7i48ZkkSehbMnhOi9frhSRJEAQBY2Nj6O/vV+6DiMhIh/7hOiAyk/Q8WQaaP8/wQpSugggw7ouLRfX09AAAAoEAent74fF4El7T09OjCj07duzA0NCQ4fUlotJ26IfXAlhIep4sy2j+/NuG14eoGBXEGJiBgQElvABAe3s7vF5vwmv8fr+qtUUQBEiSZFANiYiiDv3DdUgpvERk/OLrFuMrRFSk8j7AiKKodAMtFQgEdK8TBAEdHR1KaBFFEXa73aBaEhEBR5/bllq3UUTGM3f+Djen0BVORNryvgtJFEXN8mQtKrt370ZLSwvq6urQ19eHpqamhF1OMzMzmFm0O+zU1FTGdSai0jN9ej9mJvYnP0+ax3NfPA5TRQVuvPfeLNSMqDjlfQuMHpvNhlAopHtcEAQ4nU50dnbC7XZjaGgoYeAZGBiA1WpV/jU0NBhQayIqVsee3570nPDRGTz3xeMAgKtuvdXoKhEVtYINMInCCwA4nU7Y7XYMDQ1hbGwMoVAILS0tuuf39/cjHA4r/44dO7bSVSaiInX6wKMAZN3jsixjv/cUXuwfV8puuOuuLNSMqHjlfYDRG7ciSZLusdi4mfb2duU+9u3bB0EQ4PP5NK+xWCyora1V/SMiSsXk699OePzA7tM4+tK5SwVlZew+IlqmgggwgiBojoWJBZSlRFHUHPTb29u70tUjohIXPjwIIKJ7fFqaV4cXAPbf/32Da0VU/PI+wADR7p3FM458Pp9qWrUoispaMUA02ASDwbgxL/v27UNnZ6fh9SWi0nHy13+d8Piv//ZkXNmtO3caVR2ikmGSZVm/4zaPuN1upctoZGQELpdLOeb1euFyuTA2NqaUSZKEgYEB1NfXKzOWFi9sl8zU1BSsVivC4TC7k4hIU/jwIE7++kHd4+/uP4fhx06pyipqavBfOMuRilw2vkMLJsBkGwMMESVz5JlPYP7s25rHIgsR/PNdR+PKr/7Yx/C5l14yuGZEuZWN79CC6EIiIspHeuEFAI794qxmOWcfEa0MBhgiogxEp05ri0QiOLBbY6kHzj4iWjEMMEREGZh88zu6x8K/m9Usb9y61ajqEJUcBhgiojRNn94PRKZ1j7/tP6NZft3nPmdUlYhKDgMMEVGaJn7zd7rHIvORuHVfYmzNzUZViajkMMAQEaXp/Lu/0D323qsXNMsramqwfssWo6pEVHIYYIiI0vD2jzuAiPYYFwA48evzmuWOL3/ZqCoRlSQGGCKiFIUPD2Ju6nDCc86Nz2mWf3TXLiOqRFSyGGCIiFI0deTphMdnzsxjciy+dUbYuNGoKhGVLAYYIqIURea0F6eLef3vJzXLGz72MSOqQ1TSGGCIiFJUZfug7rH5uQXd2Ucf2rHDqCoRlSwGGCKiFFmv1V/H5b0D2uvClFVVcfYRkQEYYIiIUjQjHdQ9dnK/9uyjy1tbjaoOUUljgCEiStF06FXdY2eOa88++sRjjxlVHaKSxgBDRJSisso63WOr11fElZkqKth9RGQQBhgiohRVWpt0j2mt/2K95hoDa0NU2hhgiIhSVFnTqFl+9OdnNNd/2eJ0Gl0lopLFAENEtExvBzR2nzaZcOO992a/MkQlggGGiChFs2eOaJZrjX+pXrfO6OoQlTQGGCKiFOl1IWmNf1l/881GV4eopDHAEBGlqGrtJgjX96rKDj4jaY5/uWXnzmxVi6gklee6AkREheSyTQ+isqYRh/7xf0P8yaua2wdYbDZOnyYyGAMMEVEaTu1/BNIbHlivAjb1XIY1V1Ti9UH1Jo5X3nJLjmpHVDrYhURElKLp0/shveFRlW38rBV1TZWqspqGhmxWi6gkMcAQEaUo1VlI6zdvzkZ1iEoaAwwRUYpSmYVUabVy/ReiLGCAISJKUSqzkFq+8pUs14qoNDHAEBGl4bJND6LhU09jdvaz+Pl/O4E3BiXV8VqOfyHKCs5CIiJKU9XaTUDZbzTXfzFXVmpcQUQrjS0wREQZsDU3p1VORCuLAYaIKAPrt2zB5r4+VdnNTicXsCPKEpMsy3KuK5GPpqamYLVaEQ6HUVtbm+vqEFGeGh8eRujgQdiamxleiC7Kxncox8AQES3D+i1bGFyIcoBdSERERFRwGGCIiIio4DDAEBERUcFhgCEiIqKCwwBDREREBadgZiG53W4IggAAkCQJfUvWX9DjdDrR1NQEALDZbOjs7DSqikRERJQlBRFg3G43AKCnpwcAEAgE0NvbC4/Ho3uNJEloa2vDCy+8AEEQEAwG0dLSAi57Q0REVPgKYiG7uro6HDlyRGmBAQCTyZQwjPT29qKpqUnVUhMIBNDe3p7SY3IhOyIiosxk4zs07wOMKIpoamqKCysmkwl+v183kJhMJoyNjSn3kWpwiWGAISIiykw2vkPzfhCvKIqa5YIgQJKkhNcEg0FIkgS73Y7e3l4EAgHdx5mZmcHU1JTqHxEREeWnghgDo8VmsyEUCmkeiwUYQRDgcDgAAC6XC42NjZicnNS8ZmBgAA8//LAxlSWiojJ9ej9mzxxBZU0jqtZuynV1iEpS3rfA6NELL4u1trYqP8dabPRaYfr7+xEOh5V/x44dW7G6ElHxOLX/ERx7fhvee/mvcOz5bTi1/5FcV4moJOV9C4zdbtcsj3UNpXONIAi6XVIWiwUWiyWzShJRSZg+vR/SG+rZj9IbHtQ0bGVLDFGW5X0LjN1u1w0eegNz7XY77HZ73DWSJKlaZYiI0jF75kha5URknLwPMEC0e2dx14/P51PWhAGiY15ia8XEuFwuDA4Oqq5pb29XxsQQEaWrsqYxrXIiMk7eT6OOcbvdStfQyMgIXC6Xcszr9cLlcinTpheXx2YqTUxMqK5JhtOoiUjLqf2PqLqRhOvvw2WbnDmsEVH+4TowOcQAQ0R6OAuJKLFsfIfm/SBeIqJ8U7V2E4MLUY4VxBgYIiIiosUYYIiIiKjgMMAQERFRwWGAISIiooLDAENEREQFhwGGiIiICg4DDBERERUcBhgiIiIqOFzITkdsgeKpqakc14SIiKiwxL47jVzsnwFGx8TEBACgoaEhxzUhIiIqTBMTE7BarYbcNwOMDpvNBgA4evSoYf/z88HU1BQaGhpw7Nixot7zic+zuPB5Fp9Sea6l8jzD4TA2bNigfJcagQFGh9kcHR5ktVqL+pcspra2ls+ziPB5FpdSeZ5A6TzXUnmese9SQ+7bsHsmIiIiMggDDBERERUcBhgdFosFX/va12CxWHJdFUPxeRYXPs/iUirPEyid58rnuXJMspFznIiIiIgMwBYYIiIiKjgMMERERFRwGGCIiIio4JT8OjCSJGHv3r0YGhqC3++PO+52uyEIgnJuX19f0vvM5Jps6urqQnd3N+x2u1LPGLvdrnlNIBCAx+NBR0cH7HY7/H4/Nm/ejM7OzizUODOZ1jnfXz8tbrcbADA2NgYA8Hg8Cc8vhNezGN97WorxtVuqVN6LxfrZmrffk3IJ27dvn+zxeGSXyyU7HI644y6XS3a5XMptv98v9/T0JLzPTK7JNrvdLgOI+9fZ2al7zdDQkCwIggxAttvtssfjyWKNM5NJnQvh9Vuqr69Pdbunp0dub29PeE2+v57F+t5bqhhfOy2l8l4sxs/WfP6eLOkAEzM0NKT5wgiCIE9OTqrKkmW+TK7JtsW/ODHJ3jRDQ0NxzyvfZVLnQnj9FpucnJTb29tVdd63b58MQB4bG9O9Lt9fz2J97y1WrK+dllJ4L8pycX+25uP3JMfA6BBFEZIkxTUDAtEmv5W6JheWNk0GAgG0trbmqDb5o1Bev6VGR0chiqJyO9ZULUlSjmq0PMX83luq2F67lVKor2epfbbm+r1a8mNg9Cz+UFlMEATdD5dMrsmFxX2xoihCFEW0t7cnvW7v3r2w2WwIhUIYGxuDy+UysporIp06F8rrt5ggCJicnFSVxT4E9PrcY/L19Szm995ixfjaJVLs70WgtD5bgdy/Vxlg0hT7JTP6mmxxuVxJBw0CgMPhAHDpDer1etHV1YWhoSFD67ccK1XnfH79tAwMDMDj8Wj+hRNTiK9nsb33tBTra1eK78Vi/mxNJlvv1aIJMD6fD4ODg0nP6+/vV35hMpHJm8fIN9xynncwGEz5cZb+RXjHHXegt7dXtylwpWXyPFeqztn8wFzu77HT6UR3dzd6enoSXp/r1zMT+fbeW2nF/NoV4ntxOQrps9UIWXuvpj1qpghpDU4aGxvTHFQEQPb7/Zr3k8k1udTT06M56EzL0NBQXBkAed++fStdrRWTbp0L7fVbamhoKOUZDPn8epbCe2+pYnnt9JTae7EYP1vz8XuSAUZOPLp66WyAZJkvk2tyJdUpe5OTk3EzI2Jl+Tp6PtM6F9Lrt5jf71d9EE5OTurOZCmE17PY33uLFdtrt1SpvRdluTg/W/Pxe5KzkKDfdNXf368aFe3z+VTNu6IoKotQpXpNPhFFUXOw4NLnJQgC+vr6VOd6vV50dnbmbRNnKnUu9NcvJhgMIhgMwuFwKAMHvV4vbDYbgMJ8PYv9vRdTjK/dUqX0Xowpxs/WfPyeLOndqEVRVMYcBINB9PX1xa2A6Ha7lV+ukZER1ehwr9cLl8ulrKCZyjX5pKmpCUNDQ3FjKbSelyRJ8Hq9yu2JiYm8fV4xyepc6K8fEH2OjY2NmqP3Y2/tQn09i/m9BxT3a7dUKbwXFyumz9Z8/p4s6QBDREREhYldSERERFRwGGCIiIio4DDAEBERUcFhgCEiIqKCwwBDREREBYcBhoiIiAoOAwwREREVHAYYIioKoijmugpElEUMMERU8Nxut2ppciIqflyJl4gKXktLC1544YW82z+GiIzDFhgiKmiiKMJmszG8EJUYBhgiKmgejwe9vb25rgYRZRm7kIiooDU1NcXtdEtExY8tMERUsILBINrb23NdDSLKAQYYIso6p9OZ9By3242Ojo6E53g8HnR1dSU8p6mpCcFgULOcU6+JChcDDBFlVTAYRH19fdLzJiYmEAgEIEmS7jmBQCBhC0wgEIAoirDb7apyn8+nDP4losLEAENEWTU4OIi+vr6k57lcLgDQnV3k8/nQ2dmZ8D78fj/sdnvcfeiVE1HhYIAhoqyRJAlNTU1x5T6fT7OlJVFAGRwcTDr7SK+FJlnLDRHlPwYYIsqavXv34o477ogr37Fjh2ZryNKunxhJkjS7hpaeEwwG48bRxK5NNr6GiPIbAwwRZc2+ffvigookSZpBJBAI6IaMvXv3oru7O+FjjY6OAkBcS0tsywG2wBAVNgYYIsqaUCgUV6bXnePxeHRDhsfjQU9PT8LH8vv9EAQhLjCNjIwo419irTREVHgYYIgoa5YGBlEU4XQ64fV6Vec5nU5s3rxZ9z5S2TogEAhotux4vV44HA4A0ZacRN1QRJS/ynNdASIqHR0dHWhra0NPTw8kSYLX68W+ffuwY8cOtLS0oL29HcFgEKFQCPv27dO8D6/Xm3TtF72WFbfbDZvNpoSWsbExzkQiKlDcSoCIsqqjowOBQAAOhwO7d++Gw+FAMBhEV1cXRFFEZ2cndu/erRssUtl52ufzoaurC319faqZT52dnZAkCQMDA9i8eTM6OzvZAkNUoBhgiKhgBINBDAwMYGhoKOF5vb29CAQC3COJqIhxDAwRFYxU1n4BuM4LUSlggCGiguHz+ZIGE67zQlQaGGCIqCCk2qqit/4LERUXBhgiKggejyel7qPW1lZlDRgiKl4cxEtEBSEYDCrrtxARMcAQERFRwWEXEhERERUcBhgiIiIqOAwwREREVHAYYIiIiKjgMMAQERFRwWGAISIiooLDAENEREQFhwGGiIiICg4DDBERERWc/x+XHZRZJGGeugAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 600x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "scatter(OMEGA1,AE_Cardoso1, color=\"blue\", linestyle=\"-\", label=r\"$\\omega_B/\\mu=0.6$\", s=10)\n",
    "scatter(OMEGA2,AE_Cardoso2, color=\"darkred\", linestyle=\"-\", label=r\"$\\omega_B/\\mu=0.475$\", s=10)\n",
    "scatter(OMEGA3,AE_Cardoso3, color=\"goldenrod\", linestyle=\"-\", label=r\"$\\omega_B/\\mu=0.35$\", s=10)\n",
    "#plot(omegaNu10,MNu10/Mp**2, color=\"red\", linestyle=\"-\",label=r\"$h=10$\")\n",
    "\n",
    "#scatter(omega_0[argmax(M_0)],M_0[argmax(M_0)], color=\"k\", marker=\"s\", s=10)\n",
    "\n",
    "rcParams[\"figure.figsize\"] = (6,4)\n",
    "plt.xlabel(r\"$\\omega/\\mu$\", fontsize=16)\n",
    "plt.ylabel(r\"$A^\\mathcal{F}_E$\", fontsize=16)\n",
    "\n",
    "#plt.legend(loc=\"center left\", fontsize=14)\n",
    "\n",
    "xlim(-10,10)\n",
    "#ylim(0.9,1.2)\n",
    "plt.xticks(fontsize=12)\n",
    "plt.yticks(fontsize=12)\n",
    "\n",
    "plt.legend(loc=\"best\", fontsize=10)\n",
    "\n",
    "#plt.savefig('AFlinear_gm10_h10.pdf',bbox_inches='tight')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}