analytical_dmb.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "169171.7559705463\n",
      "zth: 169171.75597054456\n",
      "zcrit: 95615.49032033558\n",
      "kcrit: 0.4288036022565983\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "\n",
    "N0 = (2**(7/2))/(3*np.sqrt(math.pi))\n",
    "c = 3e5 # speed of light in km/s^2\n",
    "Mpc2cm = 3.086e24 # 1 Mpc = 3.086e24 cm\n",
    "omega_b = 0.044 # baryon density\n",
    "rho_crit = 4.755e-6 # critical density in (GeV/c^2)*cm^-3\n",
    "YHe = 0.24 # helium mass fraction\n",
    "Yp = 1.0-YHe # proton mass fraction\n",
    "k = 8.617e-14 # Boltzmann constant in GeV/K\n",
    "T0 = 2.73 # temperature of CMB today in K\n",
    "kT0 = k*T0\n",
    "\n",
    "B = c*Mpc2cm*N0*omega_b*rho_crit*Yp*(kT0**0.5)\n",
    "\n",
    "omega_rad = 1.0e-4 # radiation density\n",
    "H0 = 67.27 # Hubble constant in km/s/Mpc\n",
    "mp = 1.0 # proton mass in GeV/c^2\n",
    "\n",
    "mx = 1.0 # dark matter mass in GeV/c^2\n",
    "sigma0 = 1.8e-25 # cross section in cm^2\n",
    "\n",
    "zth = np.power((H0*np.sqrt(omega_rad)*np.sqrt(mx+mp)*np.sqrt(mx*mp)*((mx+mp)/mx))/(B*sigma0), 2./3.)\n",
    "print(zth)\n",
    "\n",
    "from scipy.optimize import fsolve\n",
    "\n",
    "def func(x):\n",
    "    return [(((mx/(mx+mp))*(B*sigma0*(1/(mx+mp))*(((x[0])*((1/mp)+(1/mx)))**(0.5)))) - (H0*math.sqrt(omega_rad)*(x[0]**2))),\n",
    "            ((B*sigma0*(1/(mx+mp))*((((x[1])/mp)+(((x[1])**2)/(mx*x[0])))**(0.5))) - (H0*math.sqrt(omega_rad)*(x[1]**2)))]\n",
    "\n",
    "root = fsolve(func, [1, 1])\n",
    "zth = 1/root[0]\n",
    "zcrit = 1/root[1]\n",
    "print(\"zth:\", zth)\n",
    "print(\"zcrit:\", zcrit)\n",
    "\n",
    "kcrit = (2*H0*(math.sqrt(omega_rad))*zcrit)/c # 1/Mpc\n",
    "print(\"kcrit:\", kcrit)\n",
    "\n",
    "k = 0.0857 # 1/Mpc\n",
    "bool(kcrit > k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mx: 1.000000e-06\n",
      "zth: 19724774.973122437\n",
      "zcrit: 19719.840696770367 \n",
      "\n",
      "mx: 1.000000e-03\n",
      "zth: 629432.5588197407\n",
      "zcrit: 19743.921030704474 \n",
      "\n",
      "mx: 1.000000e+00\n",
      "zth: 29848.361615531358\n",
      "zcrit: 19971.735603187735 \n",
      "\n",
      "mx: 1.000000e+03\n",
      "zth: 19973.340487138656\n",
      "zcrit: 19960.040439232096 \n",
      "\n",
      "[4.33148322e-26 2.40647515e-25 2.08122157e-24 1.34626058e-21]\n",
      "Percent errors:\n",
      "1.000000e-06 : 88.41043646918044\n",
      "1.000000e-03 : 89.81913442993874\n",
      "1.000000e+00 : 90.58245393827596\n",
      "1.000000e+03 : 95.81061787983339\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAEfCAYAAABxkKByAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXzU1b3/8dcn+75A2BNJsCAEAiirBVtblWJFpGqV1mpBK7cL2sUutvVXvbW9tb3aW1FvvSpqWxFrrVK12tqq1GqtigrI4sIqYV8CBMg2M5/fH+ebZBKSkElmTT7Px2MeyXznO9/5zBDyzjnfc85XVBVjjDEmWpJiXYAxxpjexYLHGGNMVFnwGGOMiSoLHmOMMVFlwWOMMSaqLHiMMcZElQWPMb2ciNwtIv8vgsd/RUROjdTx23i92SLySLRez4TOgseElYhsEZEaEakWkYMi8i8R+bKIJAXt86CIqIjMbvXcX3nb50Wp1nki8nI0XiscROQmEXmo1bblIvKl7hxXVb+sqjd3r7q2icj5QLWqvh3m4/YRkSdE5KiIbBWRzzc+pqpPAmNEZGw4X9OEjwWPiYTzVTUXGArcAnwPWNxqn/eBLzbeEZEU4LPAxnAV4R3TdEBEkiP8El8GfheB494F1AMDgMuAX4vI6KDHlwILIvC6JhxU1W52C9sN2AKc3WrbZCAAjPHuPwjcCuwCCr1ts4BngZeBeR0cX4FrgU3APuC/gSTvsXnAK8D/AAeAnwD5wG+BvcBW4AbcH1yjgFrADxwBDnrHaHP/oOO/7NVeBWwGzu2g1u8B24Fq4D3gLG97MvADXMhWA28CJd5jtwPbgMPe9jO87TNxv2gbvHpXAT/16q/1tt3p7TsS+Jv3GbwHXBJU04PAr4FngKPA2d62n3iPnwlUAtcBe4CdwPyg5/cFnvLqe8P7jF9u5/2nATVAcdC2m4BHvc+4GlgLTAzxZyzb+yxGBG37HXBL0P1pwOZY/3+wW9s3a/GYiFPV13G/zM4I2lwLPAnM9e5fgftl1BmfASYCpwEXAFcGPTYFF0r9cb+Y78CFyTDg497rzFfV9bi/xl9V1RxVLfCe3+b+rY7/HlAE/AJYLCLSukAROQVYCExS1/r7FC6UAb4FfA74NJDn1X/Me+wNYDzQB3gY+IOIZKjqX4D/An7v1TtOVX8I/BNY6G1bKCLZuNB52PsMPgf8b6vWwOe9zyYXF6StDfQ+gyHAVcBdIlLoPXYXLrAG4lqsX2zj+Y2GAwFVrWy1fTbwCFCA+xm4M+hze9rrom3r9rS32wjAr6rvBx1zFRD8HtcDpSKS10F9JkYseEy07MD9Mg32W+AKEcnH/ZJf1slj/VxVD6jqh8CvcL9cm15HVe9QVR/ur+JLge+rarWqbgFuAy5v66Bet9OJ9t+qqveqqh/4DTAI193Tmh9IB8pFJFVVt6hqYzfil4AbVPU9dVap6n4AVX1IVferqk9Vb/OOcUonPxdwLcctqvqAd4y3gD8CFwft8ydVfUVVA6pa28YxGoAfq2qDqj6Da02d4n0+FwE3quoxVV3nfQbtKcC1alp7WVWf8T7D3wHjGh9Q1VmqWtDObZa3Ww5wqNUxD+GCtFHj6xZg4o4Fj4mWIbiunyaq+jLQD9ed9bSq1gQ/LiJrReSIdwtuLW0L+n4rMLidx4pw3T1bW+0/pJ0aO7P/rqD6G1spOa0PpKobgG/gupb2iMgjItJYZwntnMsSketEZL2IHBKRg7iWR1E79bZlKDAluKWAOwcyMGifbW0/tcl+L7gbHcO9x35ASqvnd3SsKlqGQaNdQd8fAzJCPB93BNdSDJZHy5BrfN2DIRzXRIkFj4k4EZmE++XdVrfOQ7jzCcd1s6nqaK8LKUdV/xn0UEnQ9yfhWlNNTwv6fh/ur/ehrfbf3sa+ndk/JKr6sKpO946nwM+9h7YBJ7fe3wvX7wGX4M59FeD+km/symtrKfnW27YB/2jVUshR1a908JzO2gv4gOKgbSXt7AvwASAi0l7QH0dEng36Y6P17Vlvt/eBFBEZHvTUcbjzRY1G4Vp+hzv72iZ6LHhMxIhInojMwvXnP6Sq77Sx2yLgHOClEA79HREpFJES4OvA79vayevKeRT4qYjkishQ3PmVxiHJu4FiEUnr5P6dJiKniMgnRSQddz6rBtf9BnAfcLOIDBdnrIj0xf2V7sP9gk8RkR/R8i/73bjzFkmttg0Luv80MEJELheRVO82SURGhfoeWvM+n8eBm0QkS0RG4s6Btbd/A/B3XDdqZ1/j3KA/NlrfzvX2OerV8WMRyRaRabhzfcGj5z6OG6xi4pAFj4mEp0SkGvfX9w+BX9LyBH0T71zN86oayl/hf8KN+FoJ/Jnjh2oHuwZ3MnwTrsX1MHC/99gLuL+Sd4nIvk7sH4p03FDyfbiupf64kWzgPo9Hgedwo8MWA5nAX3G/LN/HdfHV0rIr6w/e1/0i8pb3/e3AxSJSJSKLVLUamIEbtLHDe+2fe/WEw0Jc998u3C/6pUBdB/v/H+2cU+umr+I+sz1eDV9R1eAWz+e81zZxSEL7/25MbImIAsO9cygmxkTk58BAVW13dJs3SfcaDfMk0g5e73zgclW9JBqvZ0JnwWMSigVPbHnda2nAO8Ak3HygL6lqZ0ckGkPCzuwWkTnAebgujLtU9TlvDsP/4obRLlfVJbGs0ZgeKBfXtTUY1811G67r05hOi0mLR0Tux8032KOqY4K2z8T1WScD96nqLZ04ViFwq6peJSKX42agPyUiv1fVSyP0FowxxnRRrFo8D+JmKzcNofUmp92FG+FUCbwhIk/iQuhnrZ5/paru8b6/wXseuGGejSOn/BhjjIk7MQkeVX1JREpbbZ4MbFDVTQDiljW/QFV/hmsdteAtU3IL8Kw3OxtcYBXjRjvZiD1jjIlD8XSOZwgth45W4tbFas81uAUO80XkI6p6N25s/50ich5uIcPjiMgCvFVrszIzJ5w8rDQMpRtjTO8gksTqNWv3qWq/rh4jnoLnuIUW6WCGtaouwk0+DN52lHbmiwTtcw9wD0BF+Uj940M21N8YYzorLTOX0pGnbj3xnu2Lp+6oSlouv1FMy6VQjDHG9ADxFDxvAMNFpMxbwmQubsl0Y4wxPUhMgkdElgKv4pZarxSRq7zVcBfilg1ZDzzaagkMY4wxPUCsRrV9rp3tz+BmQhtjjOmh4qmrzRhjTC9gwWOMMSaqLHiMMcZElQWPMcaYqLLgMcYYE1UWPMYYY6LKgscYY0xUWfAYY4yJKgseY4wxUWXBY4wxJqoseIwxxkSVBY8xxphOW7OnjuScvgO7c4x4uhCcMcaYOPavygb+85XDJGcXDunOcSx4jDHGdOhog/Lo+jp+v76egNL29aJDYMFjjDGmTQ1+5c8bG3hobR2H6pTpJWm8saMBVLU7x7XgMcYY04Kq8nKlj8Wr6th+JMDUIWn8YHoeYwek8ebeJKYsOrijO8e34DHGGNNkzV4f96ysY/1+P8P7pPDA7ALOHJqOiOtfmzA4E/+R/bu68xoWPMYYY9h22M/iVXW8st3HgOwkfn5WPhePyiI5qZsndNpgwWOMMb1YVW2A362p45mNDWSkCN8+PZcrx2eTlRq52TYWPMYY0wvVNCiPvVfPH96toyEAX6jI5prJORRlJUf8tS14jDGmF/EHlL9sauC3a+o4UKuce3IG352WR1lB9OLAgscYY3oBVeXfO3zct6qODw8HmDAojXun5zFhUFrUa7HgMcaYHu7d/X7uXVnL6r1+ygqS+b/zCpkxLKNppFq0WfAYY0wPtaM6wAPv1LL8Qx99M5O4+cx85o7OIjU5NoHTyILHGGN6mEN1AZasreepDfWkJgnXTs5hwWk55KTFx7rQFjzGGNND1PmUJ96v55H1ddT44NLyLL45NZf+2ZEfqRYKCx5jjElw/oDy9y0N/GZNHXuPKWeVpXP9R/MY3jc11qW1yYLHGGMSlKqyYpef+1bVsulggLH9U7ljZh5Ti9NjXVqHLHiMMSYBbahyI9Xe2u2nJC+ZO2cWct7w2I1UC4UFjzHGJJDdRwM8+E4dz29pID9d+NHH8rhsTDbpKfEfOI0seIwxJgFU1yuPrKtj2fv1iMCXJ+Tw5Yk55KfHx0i1UFjwGGNMHKv3K09+UM/D6+o5Uq9cOCqT66bmMTg3vkaqhcKCxxhj4lBAleUf+nhgdS27jiofOymd66flUd4vPkeqhcKCxxhj4szK3T7uXVnL+1UByotS+O9z8jjjpIxYlxU2CRs8IjIHOA/oD9ylqs+1tS2WNRpjTCg2H/SzeHUdr+3wMSgniV+eU8CckZkkJcBItVDEJHhE5H5gFrBHVccEbZ8J3A4kA/ep6i3tHUNVlwHLRKQQuBV4rq1tEXwbxhgTFvuOBfjNmjqe29xAdqrw/Wl5fHFcNhkJNFItFLFq8TwI3An8tnGDiCQDdwHnAJXAGyLyJC6Eftbq+Veq6h7v+xu85wVra5sxxsSVow3Ko+vr+ON79fgV5o/LZuGkXAozE2+kWihiEjyq+pKIlLbaPBnYoKqbAETkEeACVf0ZrnXUgrhZUrcAz6rqW+1tM8aYeOMLKH/e2MBDa+o4WKfMHpHJd07PpSQ/Yc9+hCSe3uUQYFvQ/UpgSgf7XwOcDeSLyEdU9e52trUgIguABQCDBw0IV+3GGHNCqsrLlT7uX11HZXWAqUPS+MH0PMYOiP7F2GIpnoKnrc5MbW9nVV0ELDrRtjaedw9wD0BF+ch2j2+MMeG0Zq+Pe1fWsW6/n+F9UnhgdgFnDk1PiCVuwi2egqcSKAm6XwzsiFEtxhgTFtsOu5Fqr1T66J+VxM/PyueiUVmkJPW+wGkUT8HzBjBcRMqA7cBc4POxLckYY7qmqjbA79bU8czGBjJShG+fnsuV47PJSk3wgQO7VjMkVwZ25xCxGk69FDgTKBKRSuBGVV0sIguBv+JGst2vqmtjUZ8xxnRVjU957N16/vBuHQ0BuKwii2sn51KUlbhL3DT54DlYfgsDc2RIdw4Tq1Ftn2tn+zPAM1Euxxhjus0fUP66uYHfvFPHgVrl3JMz+M5H8xhWGE8dS120Zz28/RBsfSUsh+sBn4gxxsSOqvLvHT7uW1XHh4cDTBiYyr1n5DNhUIKPVFOFnStd4Gx/EzIKYPzlsOYPKEe6NTDLgscYY7rovf1+7llZy+q9fsoKkrn7vEI+NSwxLsbWLlXY9m8XOLvXQnY/OOfHMPFKSM+FCVew+xdTuzXwy4LHGGNCtKM6wAPv1LL8Qx99M5O4+cx85o7OIjU5gQMn4IfNL8HKJbB/A+QXw6dvhVO/AKmZzfuVTGZ7te7qzktZ8BhjTCcdrgvw0Np6ntpQT2qScO3kHBaclkNOWgKPVAv44IO/waqlcPBD6Dsc5vwaKj4LyZG5BIMFjzHGnECdT3nig3p+v66OYz64tDyLb0zJZUBOAo9U89XBe8/AqkfgyG4YWAGffRBGzYakyL4vCx5jjGmHP6A8v7WBB9+pY+8x5ayydK7/aB7D+ybwxdjqj8G6P8E7f4CaA1A8GWbfCcPPgSidm7LgMcaYNqzY6eO+VbVsPBhgbP9U7piZx9Ti9FiX1XW1h2DN47D2Cag7DMM+AR/7NgydFrXAaWTBY4wxQTZU+blvZS1v7vZTkpfMHTMLmTU8gUeqHdsPqx+FdU+CrwZGzoIzvgVDJsSsJAseY4wB9hwN8MA7dTy/pYH8dOFHH8vjsjHZpCfqxdiqd7rzN+896wYQjLkIpn8LBpTHujILHmNM73akXlm6ro5l79eDwH9MyOErE3PIT0/QkWoHt8LbD8OGv4MkwfjPw/RvQJ9hsa6siQWPMaZXqvcrT22o5+G19VTXKxeOyuRbU3MZkpugvxb3fQBv/w42/xNS0mHKf8DpCyG/W8uqRUSCfsLGGNM1AVWWf+jjgdW17DqqfOykdK6flkd5vwQdqbbrHbfKwLbX3MoCZ3wLpn4VsotiXVm7LHiMMb3Gqt0+7llVy/sHApQXpfDf5+RxxkkZsS4rdKqwfYULnJ2rIKsvfPL/weSrISM/1tWdkAWPMabH23LIz32r6nhth49BOUn88pwC5ozMJCnRRqppALa84pa12fsu5A6CmbfAaVdAWnasq+s0Cx5jTI+1rybAb9+p46+bG8hOFa6flsu8cTlkJNpItYAPNr4IKx+Gqs1QWArn3w7jPufO5yQYCx5jTI9ztEH5w/o6HnuvHr/C/HHZLJyUS2Fmgo1U89fD+39166gd3gH9RsKF98Hoz0By4v76TtzKjTGmFV9A+fPGBh5aU8fBOmX2iEy+c3ouJfkJ9quuoQbWPw3v/B6O7oPBp7mVokecC0kJFp5tSLB/DWOMOZ6q8nKlj/tX11FZHWDqkDR+MD2PsQMS7GJsddVuSZs1f3RL3JROhwu/DcPOjPqyNpFkwWOMSWhr9vq4d2Ud6/b7Gd4nhfvPL+ATpemJtcRNTRW88xisWwb1R2H4DDjj23DSlFhXFhEWPMaYhLTtsJ/Fq+t4pdJH/6wkfn5WPheNyiIlKYEC58geWP17ePdp8NVD+QVwxnUwaGysK4soCx5jTEKpqg3wuzV1PLOxgfQU4bqpuVx1ajZZqQl07uNQpRuh9sFz7v7YS2H6N6FoeGzrihILHmNMQqjxKX98t55H362jIQCXVWRx7eRcirIS6GJsBzbB20tg04uQlAoT5sO0a6HgpFhXFlUWPMaYuOYPKH/d3MBv19Sxv0Y59+QMvvPRPIYVJtCvrz3r3SoDW19xEz1PX+huuQNiXVlMJNC/nDGmN1FVXtvh475VdWw9HGDCwFTumZXPhEEJMlJNFXaudIGz/U3IKIAzvw+TF0BWn1hXF1MWPMaYuPPefj/3rKxl9V4/ZQXJ3H1eIZ8aliAXY1OFD191y9rsXgs5/eGcm2HifLeIp7HgMcbEj51HAty/upblH/rom5nEzWfmM3d0FqnJCRA4AT9sfskFzv4NkF8C590G478AqQm4EGkEWfAYY2LucF2AJWvreXJDPalJwrWTc1hwWg45aQkwUi3ggw/+5kapHdoGfYfDnF9DxWchOUEvtRBhFjzGmJip8ylPfFDP79fVccwHl5Zn8Y0puQzISYCRar46ePfPbh7Okd0wsAJm/ARGnQ9JCVB/DFnwGGOiat0+H2/v9tHgh79taWDPMeWssnSu/2gew/smQAuh/iisexLe+QPUHICSKTD7Thh+To9a1iaSLHiMMVFxpF55YWs9v36rDp+6bcMKk7l9ZgGnFyfA0v61h2DN47D2cbem2smfdKsMDJ1mgRMiCx5jTFj5A8rOIwE2Hgyw6aCfTQcDbD7oZ/cxbbFfEnDhyKz4D51j+1132rqnwFcDI2e5y0sPmRDryhKWBY8xpsuO1iubD/ldyFT52XTQz5ZDAWr97vFkgWGFKUwpzmBUUQqpSfDfr1bjC0BqMvEdOod3uuvgvP8XN4BgzMVuWZsB5bGuLOFZ8BhjTiigys4j6rVgXCtm00E/u442t2IK0oVR/VI5/aRURhWlUl6Uwkf6pB53tc9TB6Xz78o6phanx+dk0KqtboTahr+7QQLjPw/Tvg59hsW6sh4jYYNHROYA5wH9gbtU9TlvezbwEnCjqj4dwxKNSUjHGlwrZlNVgI1eyGw55KfG5x5PEigrSGbi4HRGFaU23QbmJHVqgueEQWnxGTj73nerDGz+p7uc9JT/gI9eA3mDY11ZjxOT4BGR+4FZwB5VHRO0fSZwO5AM3Keqt7R3DFVdBiwTkULgVsBb5pXvAY9GqnZjeoqAKruPKhsP+tl8MND0dceRQNM+eenCqKJULvW6ykYVpTK8TwqZibQS9InsWu0CZ9vrbmWBM66DqV+B7KJYV9ZjxarF8yBwJ/Dbxg0ikgzcBZwDVAJviMiTuBD6WavnX6mqe7zvb/Ceh4icDawDbJqwMUFqGlsxwSf8D/k51uAeF6C0IJlTB6bx+X6pjOybyqh+KQzOSU6MZWpCpQqVb7jA2bUasorgrB/BpC9BRn6sq+vxYhI8qvqSiJS22jwZ2KCqmwBE5BHgAlX9Ga511IK4/w23AM+q6lve5k8A2UA5UCMiz6hqoPVzjemp1GvFBJ+H2eS1YhrPxuSkCSP7pnLxqOaushF9UxLrejZdpQHY8ooLnH3vQe4gmHkLnPZFSMuKdXW9Rjyd4xkCbAu6Xwl0dN3Xa4CzgXwR+Yiq3q2qPwQQkXnAvrZCR0QWAAsABg/qnUuSm56h1qdsOeTCpbGbbPNBP0eCWjEn5SczdkAal45JZWTfFEb1S6U4t4e2YjoS8MHGF9yggaotUFgG5y+CcXPd+RwTVfEUPG39T9A2trkHVBcBi9p57MEOnncPcA9ARfnIdo9vTLxQVfYe06YT/Y2tmO3Vza2YrFRhVFEKc0Y2t2JO6ZtCdiKsdRZJ/np47y9uWHT1Tug/Ci5aDOVzIDmefv31Lp365EWkMxePCKjqwW7UUgmUBN0vBnZ043jGJJw6n7LlcOOcGBcymw8FqK5v/hupJC+Z0f3TuGiUa8GUF6VSnJdMUm9rxXSkoQbWPwXvPApH98HgCXDeL2HETEjq5WEcBzob+Tu8W0c/2clAd67f+gYwXETKgO3AXODz3TieMXFLVdlXc/y5mMrqAAEvY7JShFOKUjh/hGvFlPdLZUSfFHLT7Rdnu+qqYe0TsOaPbomb0ulw4bdh2Jm2rE0c6WzwrFfVUzvaQUTe7uyLishS4EygSEQqcXNuFovIQuCvuBC7X1XXdvaYxsSrer+y9VDznJjN3tfDQa2Y4txkRvVLY87IlKauspPyrRXTaTVVbtHOtcug4RgM/5QbFn1SR6eJTax0NnhOD9M+AKjq59rZ/gzwTGePY0w8UVX216rrJjvU3F1WWR3A72VMRgqc0jeV84Z752L6pXBK31TyrBXTNUf2wKpH3OUJ/PUweg5M/xYMGhvrykwHOhU8qlobjn2M6Snq/cqHhwPHdZUdqmtuxQzOSaa8Xyrnn9I8u39ofjLJSdaK6bZDlW6E2gfevPGxl7p11IqGx7Yu0yknDB4ROQe4BLcszUoRWeCNDDOmVzhQE2gRLpsO+vnwcHMrJj0ZRvRNZebJ6Yzql+ImXxalkp9hrZiwO7DJzcHZtBySUmHCfJh2LRR05/SyibbOtHi+CswHbvBGt42PbEnGxEaDX9lWHWBjqxFlVbXNrZiBOUmUF6Vy7vDGVkwKZQUp1oqJtD3r4K2H4MN/QVq2W0Nt6tcg1+biJaLOBM9eb5j0t0XkFmBShGsyJuKqagNN65MFt2J83pTjtGQY0SeVs4elMaootakVU5hprZioUYUdb8PKJbD9TcgshDN/AJOvhqzOzPAw8aozwfPnxm9U9XoRuSaC9RgTVr6Asu1wcFeZC5oDQa2Y/tmuFTPjZNeCKS9KpawwhRRrxcSGKnz4Kry9BPashZz+cM7NMHG+W8TTJLwTBo+q/qnV/TsiV44xXXeorjlgNla5RTC3HgrQ4LViUpNgeJ8UPlGWEbScfwp9MpNjW7hxAn7Y/A8XOAc2Qn4JnHcbjP8CpNq6vz1JSGtGiMhE4IfAUO+5Aqiq2thFEzX1fuWf2xp4dYePVIFD9bD5oJ99Nc2tmH5ZSYwqSuWsMm9eTL9UhhWkkJpsrZi442+AD/4Gqx52o9WKhsOcu6HiYkhOjXV1JgJCXaxoCfAd4B3AVn02EeMLKLuOBNh+xK1JVlntvu44EmD3UW2xiN/QvGQ+NtS1YkZ614wpyrJWTNzz1bn5N6t/D0d2w8CxMOOnMOp8d+VP02OFGjx7VfXJiFRieh1/QNlzTJtCpTFkdlQH2HW0ebgyQG6aUFqQwuQhaew75ufVynoUSBa4ZHQWX5tkff8JQRW2/NOto7bnXaivhpKpcMFd8JGzbVmbXiLU4LlRRO4DngfqGjeq6uNhrcr0GAFV9jWGixcsjS2XnUeaz7+AW5ustCCZ8QPTKCtMobTADVUuzU+mT2bzZZXf3FnPZU/so8EPqckwtdiWtY9LddVu3k3TbTPs3wA+b665JMGnb4PJX4ptnSbqQg2e+cBIIJXmrjYFLHh6scalYhpDpbH1ssP7Wu9v3jc9GYbmpzCqXxqfHp5CWUFyU8D0y0rq1HViJgxKY8lnivh3ZR1Ti9OZMCgtgu/OnJC/Hg5+6MJl/yao8kLm6N7mfTLyoX85DBjjhkajgEDdoVhVbWIo1OAZp6oVEanExDVV5WCdtgiWxm6x7UcC1Pia901NchcgO7lPKmcPC2q5FKQwMCcpLAtfThiUZoETbRqA6t0tWzFVm+HgNlDvr4ukVOh3Cpz8SS9oRruveYNdN9q21+E3s11YJadB6RmxfU8mJkINnn+LSLmqrotINSbmDtcp24/422y9NF7ZEty5leK8ZMoKUjmjNIWy/BRKC5IpK0hhcK6tR5bwag81d481hcwWt/Jzo4KToP9oGHNRc8j0/UjHI9FKJsMXn3TneUrPcPdNrxNq8EwHvigim3HneGw4dQI62nB8t1jjeZfgRS4FGJLrwmVqcXKLlktxbrINTe4JfHVwcGvL8zAHNsOxfc37ZBa6YCmdHtSKGdX1yZwlky1werlQg2dmRKowYVfj06ZusOCQ2XGk5ax9cOuPlRWkcOqg9KZgKStIpiQvhfQUC5ceQQPu0s/7W3WTHap0jwEkp7tusuHneAFT7lo0uQNttJkJq1CD58fA1xsvcS0ihcBtwJXhLsycWL1f2dEqWBqDJngyJbgJlaUFKZw9LJ0yr0ustCCFofnJZKba+mM9Ss3B48/DHNgCvhpvB4HCUjdvZuzc5oDpMwySQ/2VYEzoQv0pG9sYOgCqWiUiHV6Z1HRPg1/ZdbRlqDR+v6fVRMrCDKGsIIWPl6Z5w5BduJQWJJOTZuHS4/hq3XmX4PMwBzZDzYHmfbL6utZL2cebu8n6jYT0nJiVbUyowZMkIoWqWgXgXSbB/kTqJn9A2X1UWwWLn+3Vyu5WEynz0oXS/BROL05r6hIr9Vov+XYVy54p4IfDO1qGS9UmOLQdGv/0SMlwgXLKuS27yXL6WzeZiTuhhsZtwL9E5DHcT/wlwE/DXjWoly4AABl+SURBVFUPFFBl77HjhyNvrw6w82jzcvwAWamu5XLaoOYuscaQKczo3FwXk6COHWg5mqxqs7v5Gudri+sSG3yqWzyzqZuszJaZMQkjpOBR1d+KyArgk7hBTxfa0Opmqsr+Gm2zW2xHdYD6oHDJSHETKcf0T2NW02gx13rp7ERKk8AaarxusuBWzGaoqWreJ7uf6xo7+SwvYMpdqyYtK2ZlGxMOIXeTeUHTa8MmeCJlZatg2dHmRMoUPtInjRkne11i+S5kBoRpIqWJcwEfHN5+/HmYwzto6iZLzYR+5TByVvOEywGjIbsopqUbEymdCh4ReUtVT+vuPonkcF2gzcUrtx8JcDRoImVKEhTnJlNWmMrHS1u2XAbn2ETKXkMVju0//jxM1VY3Sx/c2mR9hkHxRNc91tiKKSyDJDs/Z3qPzrZ4RonI6g4eFyA/DPVE1b6aAM9uqic9WY7rFjtc33xGP0lgsDeR8vSS5mApK0ihOC/ZrlTZ29Qf84Yob245ZLk2aN2xnAGu1TL8U82tmH6nuNaNMb1cZ4NnZCf28Z94l/hyoFb55eu1TfcHeRMpJwxO95aAcSFjEyl7qYDPTbBsPau/ekfzPmnZLlTK57TsJsvqE7u6jYlznQoeVd0a6UJiJQn44rhsvjctjwwLl95J1a2k3LoFU7UVAl6/qiS7dchOmto8kmxAOeSfZN1kxoSo18/BSUuBWSMyLXR6i/qjrVow3uKXdYeb98kd5Fotp3zafR0wGopGQIpd98eYcOjVwTMgO5klnymy5fV7ooDPu0ZMq1ZM9a7mfdJzXdfYmIuCusnK3aKYxpiI6dXB0z87yUInUe1eCztWwqBxkNPv+FbMwQ9d+AAkpUDRcBg6vVU3WYnN6jcmBrocPCLyeWA2blCBAE+p6tJwFWZMC/XH3En9wzvdFSzXP9m8qnKwvGLXeim/oDlg+g6HFPsDw5h40Z0Wz8dVdW7jHRG5C7DgMV2jATi6z02srN7pvjZ+X72z5Yz+FsStTzbt6+4aMRkJN6rfmF6nO8GTLiLnAduAYsAmKJiONdR4oeIFS2MLpvFrIGhmriRDfrFbg6xkslvGv7DMfT22Hx65rPnyydO/aRcWMyaBhBQ8InKjqv6nd/erwIVABS58Foa5NpNoVN2S/Ie3N4dLcKvl2P6W+6fnujAZNB5Gf8YLl1K3Lb+440so2+WTjUlYobZ4bhSRLKAP8BbwSOMlEkwv4atzI8OqdzQHS2OrpXpn0CrKAAJ5Q1yrZfCpLlT6lDWHS2Zh10/u2+WTjUlYoQaPArXAX4HTcJdImKuqq8JemYkNVag92LK10tiCqd7pJloGS81yYTJgjFvksjFU+nitFpv7YoxpJdTgeVdVb/S+f0xEHgTuxl0mIapEZA5wHtAfuEtVnxORJOBmIA9Yoaq/iXZdCcFfD9W7g1otO1t+33SJZE/uIBcow2e07A4rLHUrKNuQZGNMCEINnn0iMkFV3wRQ1fdFpF+oLyoi9wOzgD2qOiZo+0zgdiAZuE9Vb2nvGKq6DFgmIoXArcBzwAXAEOAAUBlqXT2GqpuJ36LVEvT9kT0QfNHslAwXIkWnuEUtm7rDSqFgKKRmxOZ9GGN6pFCD51rgERF5E3gHGAts7sLrPgjcCfy2cYOIJAN3AefgQuMNEXkSF0I/a/X8K1V1j/f9Dd7zAE4BXlXV//Oukvp8F2pLDAEfHNntncTf3nJ0WPUOtzRMsJz+rpUy7BPNodIYMDkDrNVijImaUK9AukpExgNnA2OAF+nC3B1VfUlESlttngxsUNVNACLyCHCBqv4M1zpqQdwlOm8BnlXVt7zNlYB38ZPEWy37OHXVbbRYdsDhXS50NOgtJqdBwUnuei8nf6Jld1jhULeKsjHGxIGuXIG0DvizdwunIbhh2Y0qgSkd7H8NLgDzReQjqno38Dhwh4icAbzU1pNEZAGwAOCkIQPDUXfXBXzuZH3reS2NYRO8cCVAVl8XJkM/2rI7rLDMnYexVZKNMQkgntZqa6uvR9vY5h5QXQQsarXtGHBVRy+iqvcA9wBMHFfe7vHDpv5o2+dZDu+EIzshENRqSUqBghIoHAal01t2hxUMhYy8iJdrjDGRFk/BUwmUBN0vBna0s2/8CF7qpa3Z+LUHW+6fWeiCpGRSUFdYqQuYvCGQlByDN2GMMdETT8HzBjBcRMqA7cBc4POxLcnTtNTLjuaT+Y2z8dtb6qWwFEqmtDrXUgqZBbF5D8YYEydiEjwishQ4EygSkUrgRlVdLCILcZNTk4H7VXVtRAs5dsAtr99/lPs+eC5LcOul5kDL5wUv9VI+p+X5lvySjpd6McaYXk5UI3+aI15NHJysKxbkunMrwa2W4KVemk7gl3r3u7nUizHGJDgReVNVJ3b1+fHU1RYjCsUTYMzFzV1iBSW21IsxxkRILw8egZRMOOdmW3DSGGOipHdP/Mgb5JbXt9Axxpio6d3BkzPAQscYY6KsdwePMcaYqLPgMcYYE1UWPMYYY6LKgscYY0xUWfAYY4yJKgseY4wxUWXBY4wxJqoseIwxxkSVBY8xxpiosuAxxhgTVRY8xhhjosqCxxhjTFRZ8BhjjIkqCx5jjDFRZcFjjDEmqix4jDHGRJUFjzHGmKiy4DHGGBNVFjzGGGOiyoLHGGNMVFnwGGOMiSoLHmOMMVFlwWOMMSaqLHiMMcZElQWPMcaYqLLgMcYYE1UWPMYYY6LKgscYY0xUWfAYY4yJqoQNHhGZIyL3isifRGSGt+0kEXlSRO4XketjXaMxxpjjxSR4vGDYIyJrWm2fKSLviciGEwWHqi5T1auBecCl3uYRwJ9V9UqgPBK1G2OM6Z5YtXgeBGYGbxCRZOAu4FxcaHxORMpFpEJEnm516x/01Bu85wG8DcwVkReAFyP+LowxxoQsJRYvqqoviUhpq82TgQ2quglARB4BLlDVnwGzWh9DRAS4BXhWVd/yNs8HbvSO/xjwQITegjHGmC6Kp3M8Q4BtQfcrvW3tuQY4G7hYRL7sbfsLcK2I3A1saetJIrJARFaIyIq9e/d2v2pjjDEhiUmLpx3SxjZtb2dVXQQsarVtDXBxRy+iqvcA9wBMnDix3eMbY4yJjHhq8VQCJUH3i4EdMarFGGNMhMRT8LwBDBeRMhFJA+YCT8a4JmOMMWEWq+HUS4FXgVNEpFJErlJVH7AQ+CuwHnhUVdfGoj5jjDGRE6tRbZ9rZ/szwDNRLscYY0wUxdPgAmOM6ZKGhgYqKyupra2NdSk9SkZGBsXFxaSmpob1uBY8xpiEV1lZSW5uLqWlpbgpfqa7VJX9+/dTWVlJWVlZWI8dT4MLjDGmS2pra+nbt6+FThiJCH379o1IK9KCxxjTI1johF+kPlMLHmOMCZMnnngCEeHdd9/t8jHmzZvHY4891uE+//Vf/9Xi/kc/+tEuvdZNN93Erbfe2qXndocFjzHGhMnSpUuZPn06jzzySERfp3Xw/Otf/4ro64WbBY8xpld6c2sVd724gTe3VoXleEeOHOGVV15h8eLFTcGzfPlyzjzzTC6++GJGjhzJZZddhqpbqevHP/4xkyZNYsyYMSxYsKBpe6Pnn3+ez3zmM033//a3v3HhhRdy/fXXU1NTw/jx47nssssAyMnJadrvF7/4BRUVFYwbN47rr3dXl7n33nuZNGkS48aN46KLLuLYsWNhec9dZaPajDE9yn8+tZZ1Ow53uE91bQPv7qomoJAkMHJgLrkZ7Q8ZLh+cx43nj+7wmMuWLWPmzJmMGDGCPn368NZbbtH8t99+m7Vr1zJ48GCmTZvGK6+8wvTp01m4cCE/+tGPALj88st5+umnOf/885uO98lPfpKvfe1r7N27l379+vHAAw8wf/58zj//fO68805Wrlx5XA3PPvssy5Yt47XXXiMrK4sDBw4AcOGFF3L11VcDcMMNN7B48WKuueaaDt9PJFmLxxjT6xyu9RHwGhgBdfe7a+nSpcydOxeAuXPnsnTpUgAmT55McXExSUlJjB8/ni1btgDw4osvMmXKFCoqKnjhhRdYu7blQi0iwuWXX85DDz3EwYMHefXVVzn33HM7rOHvf/878+fPJysrC4A+ffoAsGbNGs444wwqKipYsmTJca8VbdbiMcb0KCdqmYDrZrvsvn/T4AuQmpLE7XNPZcLQwi6/5v79+3nhhRdYs2YNIoLf70dE+PSnP016enrTfsnJyfh8Pmpra/nqV7/KihUrKCkp4aabbmpz2HJjCycjI4PPfvazpKR0/CtbVdsciTZv3jyWLVvGuHHjePDBB1m+fHmX32s4WIvHGNPrTBhayJIvTeVbM05hyZemdit0AB577DGuuOIKtm7dypYtW9i2bRtlZWW8/PLLbe7fGDJFRUUcOXKk3VFsgwcPZvDgwfzkJz9h3rx5TdtTU1NpaGg4bv8ZM2Zw//33N53Daexqq66uZtCgQTQ0NLBkyZLuvNWwsOAxxvRKE4YW8rVPfKTboQOumy14IADARRddxMMPP9zm/gUFBVx99dVUVFQwZ84cJk2a1O6xL7vsMkpKSigvL2/atmDBAsaOHds0uKDRzJkzmT17NhMnTmT8+PFNQ6VvvvlmpkyZwjnnnMPIkSO7+jbDRlqPpOhNJk6cqCtWrIh1GcaYblq/fj2jRo2KdRkRsXDhQk499VSuuuqqmLx+W5+tiLypqhO7ekw7x2OMMXFqwoQJZGdnc9ttt8W6lLCy4DHGmDj15ptvxrqEiLBzPMYYY6LKgscYY0xUWfAYY4yJKgseY4wxUWXBY4wxYSAiXHfddU33b731Vm666aao1rBixQquvfbaLj33zDPPJFrTSyx4jDEmDNLT03n88cfZt29fTF7f5/MxceJEFi1aFJPXD4UFjzGmd9r2OvzzNvc1DFJSUliwYAH/8z//c9xjrS/u1ngZg+XLl/Pxj3+cSy65hBEjRnD99dezZMkSJk+eTEVFBRs3bgRg7969XHTRRUyaNIlJkybxyiuvAO5CbgsWLGDGjBlcccUVLF++nFmzZgHuMg3z58+noqKCsWPH8sc//hGAr3zlK0ycOJHRo0dz4403huW9h8rm8RhjepZnr4dd73S8T91h2L0GNACSBAPGQHpe+/sPrIBzbznhS3/ta19j7NixfPe73+10uatWrWL9+vX06dOHYcOG8aUvfYnXX3+d22+/nTvuuINf/epXfP3rX+eb3/wm06dP58MPP+RTn/oU69evB9xcn5dffpnMzMwWi3/efPPN5Ofn88477rOoqnLXHfrpT39Knz598Pv9nHXWWaxevZqxY8d2ut5wsOAxxvQ+tYdc6ID7Wnuo4+DppLy8PK644goWLVpEZmZmp54zadIkBg0aBMDJJ5/MjBkzAKioqODFF18E3OUO1q1b1/Scw4cPU11dDcDs2bPbfK2///3vLa6EWljo1qR79NFHueeee/D5fOzcuZN169ZZ8BhjTLd0omXCttfhN7PBXw/JaXDRfVAyOSwv/41vfIPTTjuN+fPnN21LSUkhEHBBp6rU19c3PRZ82YSkpKSm+0lJSfh87jpBgUCAV199tc2Ayc7ObrOOti6RsHnzZm699VbeeOMNCgsLmTdvXpuXY4g0O8djjOl9SibDF5+ET/7QfQ1T6IC7+Noll1zC4sWLm7aVlpY2LX/zpz/9qc1LGnRkxowZ3HnnnU3327r66ImeU1VVxeHDh8nOziY/P5/du3fz7LPPhlRHuFjwGGN6p5LJcMZ1YQ2dRtddd12L0W1XX301//jHP5g8eTKvvfZau62U9ixatIgVK1YwduxYysvLufvuu0/4nBtuuIGqqirGjBnDuHHjePHFFxk3bhynnnoqo0eP5sorr2TatGkhv7dwsMsi2GURjEl4PfmyCLEWicsiWIvHGGNMVFnwGGOMiSoLHmOMMVFlwWOM6RF68/nqSInUZ2rBY4xJeBkZGezfv9/CJ4xUlf3795ORkRH2Y9sEUmNMwisuLqayspK9e/fGupQeJSMjg+Li4rAfN2GDR0RGAV8HioDnVfXXIpIN/C9QDyxX1SWxrNEYEx2pqamUlZXFugzTSTHpahOR+0Vkj4isabV9poi8JyIbROT6jo6hqutV9cvAJUDjePILgcdU9WpgdkSKN8YY0y2xOsfzIDAzeIOIJAN3AecC5cDnRKRcRCpE5OlWt/7ec2YDLwPPe4cpBrZ53/uj8D6MMcaEKCZdbar6koiUtto8GdigqpsAROQR4AJV/Rkwq53jPAk8KSJ/Bh4GKnHhsxIbOGGMMXEpns7xDKG5tQIuRKa0t7OInInrWksHnvE2Pw7cKSLnAU+187wFwALvbl3r7r4O5AOHOvlYR/cbv2+9De9+EdDZSxiGUlN7dZyovnwgNYSaOqornDWF67M6UU2drSWcNZ2ohkSsqXEbEa6pvTpOVJ/9nIf2M3VKJ+tpm6rG5AaUAmuC7n8WuC/o/uXAHRGuYUUI+97T2cc6ut/4fettQdsjUlN7dZyoPq+2TtfUUV3hrClcn9WJagqhloj++yV6TUE/S/Zz3st+ztu6xVOLpxIoCbpfDOyIUS1tabMF1c5jHd1/qoNtkaypo9fsqL6ngNPCVFc4awpVV2vqbC3hrOlENSRiTV2tJ9SaOnpN+znveFukfqaOE7PVqb1zPE+r6hjvfgrwPnAWsB14A/i8qq6NYA0rtBsrrEaC1dR58ViX1dQ5VlPnxWNd3a0pVsOplwKvAqeISKWIXKWqPmAh8FdgPfBoJEPHc0+Ej98VVlPnxWNdVlPnWE2dF491daumXn09HmOMMdFnQ46NMcZElQWPMcaYqLLgaYeIJInIT0XkDhH5YqzrATd3SUT+KSJ3e/OY4oKIZIvImyLS5kTfaBORUd5n9JiIfCXW9TQSkTkicq+I/ElEZsS6HgARGSYii0XksRjXkS0iv/E+n8tiWUujePlsgsXpz1DI/996ZPCEYy044ALcpNYG3FDveKhJgSNARhzVBPA94NHu1hOumrTtdfzioa5l6tYRnAdcGic1bVLVq7pbSxjqi8o6i6HUFMnPphs1hfVnKEw1hf7/rTuTgOL1BnwMNyY/eIJqMrARGAakAatwa8JVAE+3uvUHrgf+w3vuY3FSU5L3vAHAkjip6WxgLu4/wqx4qMl7zmzgX7gh+XHxMxX0vNuA0+Kspm7/jHezvu8D4719Hg53LV2pKZKfTRhqCsvPULhqCvX/WzxNIA0bDcNacCJSibu8AoRhwdFw1BSkCrdUUMxrEpFPANm4Xx41IvKMqgZiWZN3nNbr+HVLmD4rAW4BnlXVt+KhpkgKpT6itM5iiDWti1QdXa1JRNYTxp+hcNQErAv1/1uPDJ52hLQWHG7dtztE5AzgpXioSUQuBD4FFAB3xkNNqvpDr7Z5wL7uhE64apK21/GLhFB/pq7BtRDzReQjqnp3rGsSkb7AT4FTReT7XkBFUnv1LeIE6yxGu6YYfDYnrIno/AyFVFNX/r/1puCRNra1O4lJVY8Bke7fDbWmx3GBGEkh1dS0g+qD4S+lSaif03JgeaSKCRJqXYtwv2AjKdSa9gNfjlw5x2mzPlU9CsyPYh3B2qsp2p9NsPZqisbPUHvaq2k5If5/65GDC9oRj2vBWU2dE481QXzWFY81BYvH+qymzglbTb0peN4AhotImYik4U6IP2k1WU3dEI91xWNNweKxPqsp2jVFcrRGrG7AUmAnzUOhr/K2fxq3EOlG4IdWk9WUyHXFY03xXp/VFB812Vptxhhjoqo3dbUZY4yJAxY8xhhjosqCxxhjTFRZ8BhjjIkqCx5jjDFRZcFjjDEmqix4jDHGRJUFjzHGmKiy4DEmxkSkVERqRGRl0LYBIvKwiGwSd3XXV0XkMyc4znIR+VSrbd8Qkf8VkUwRWSki9SJSFKn3YkxnWPAYEx82qup4aLpuzzLgJVUdpqoTcOtiFZ/gGEu9/YLNBZaqao13/FgvNGmMBY8x3SUifxCRO0XkZRHZKiLTReS3IvK+iCzuwiE/CdRr0LVWVHWrqt4R9JpfEJHXvVbM/4lIMvAYMEtE0r19SoHBwMvdeX/GhJsFjzHdVwFsUtXpwG+AxcD3gDHAhY1BEILRQLtXlxSRUcClwDSvFeMHLlN3/ZjXgZnernOB36styGjiTG+6EJwxYSciGbgrwv7K21QDLFbVnd7jx2i+hHpXX+MuYDquFTQJOAuYALzheuXIBPZ4uzd2t/3J+3pld17bmEiwFo8x3TMaeEubL/k9DngNQEQaL5Q1VkSeEJFUEXmoEyf31wKnNd5R1a/hwqaft0mA36jqeO92iqre5D22DDhLRE4DMlW13ZaTMbFiwWNM91QAq4LujwVWe9+PA1ar6ipgJfAA8CtV3XeCY74AZIjIV4K2ZQV9/zxwsYj0BxCRPiIyFEBVj+AuQ3w/rvVjTNyx4DGmeypwodLY7ZapqlXeY8EhVAWsV9UVJzqgd05mDvBxEdksIq/jzh19z3t8HXAD8JyIrAb+BgwKOsRSXOg90s33ZkxE2IXgjIkwESkHrgP2q+p323i8FHhaVcdEoZYtwMROtLqMiRhr8RgTQSKSD/wQ+CpQJiJT2tjND+QHTyCNQB2Z3vFTgcCJ9jcmkqzFY4wxJqqsxWOMMSaqLHiMMcZElQWPMcaYqLLgMcYYE1UWPMYYY6LKgscYY0xUWfAYY4yJKgseY4wxUWXBY4wxJqr+P+4g55o9kajtAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import fsolve\n",
    "\n",
    "k = 0.0857 #1/Mpc\n",
    "\n",
    "def k_crit(mx,sigma0,omega_b=0.05,omega_rad=9.23640e-5,\n",
    "                 H0=67.27,kT0=2.353e-13,rho_crit=4.755e-6,\n",
    "                 cc=3e5,Mpc2cm=3.086e24,reduction_factor=1.,mp=1.,YHe=0.25,n=0.):\n",
    "    N0 = (pow(2.,(7./2.)))/(3.*np.sqrt(np.pi))\n",
    "    B = cc*Mpc2cm*N0*omega_b*rho_crit*(1.-YHe)*np.sqrt(kT0)\n",
    "    zth = np.power((H0*(np.sqrt(omega_rad))*(np.sqrt(mx+mp))*(np.sqrt(mx*mp))*((mx+mp)/mx))/(B*sigma0), 2./3.)\n",
    "    #print(zth)\n",
    "    H = H0*np.sqrt(omega_rad)\n",
    "    def f(z):\n",
    "         return(B*(z**3.)*sigma0*(1./(mx+mp))*(((z/mp)+((z**2.)/(mx*zth)))**(0.5))) - (H*(z**2.))\n",
    "    zcrit = fsolve(f, zth)\n",
    "    #print(zcrit)\n",
    "    kcrit = (2*H0*(np.sqrt(omega_rad))*zcrit)/cc\n",
    "    return zth, zcrit, kcrit\n",
    "\n",
    "masses = np.array([1e-6, 1e-3, 1., 1e3])\n",
    "sigmas = np.zeros(len(masses))\n",
    "sigma_array = np.logspace(-18,-30,1000)\n",
    "\n",
    "for i in range(len(masses)):\n",
    "    mx=masses[i]\n",
    "    print(\"mx:\", format(mx, 'e'))\n",
    "    \n",
    "    zdec, zdecprime, kdec = k_crit(mx=mx, sigma0=sigma_array)\n",
    "    sigmas[i] = sigma_array[np.argmin(np.abs(kdec-k))]\n",
    "    print(\"zth:\", zdec[np.argmin(np.abs(kdec-k))])\n",
    "    print(\"zcrit:\", zdecprime[np.argmin(np.abs(kdec-k))], \"\\n\")\n",
    "    \n",
    "print(sigmas)\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "proton_chains = np.array([5.02e-27, 2.45e-26, 1.96e-25, 5.64e-23])\n",
    "\n",
    "f,ax = plt.subplots()\n",
    "\n",
    "ax.plot(masses, sigmas, marker='.', label='Analytical')\n",
    "ax.plot(masses, proton_chains, marker='.', label='Numerical')\n",
    "ax.set_title('DM-proton scattering (n=0)')\n",
    "ax.set_xlabel('$m_\\chi$ $[\\mathrm{GeV}]$')\n",
    "ax.set_ylabel('$\\sigma_0$ $[\\mathrm{cm}^2]$')\n",
    "ax.set_xscale('log')\n",
    "ax.set_yscale('log')\n",
    "ax.set_xlim([1.0e-6, 1.0e3])\n",
    "ax.set_ylim([1.0e-30, 1.0e-20])\n",
    "ax.fill_between(masses, sigmas, 1.0e-20, alpha=0.2)\n",
    "ax.fill_between(masses, proton_chains, 1.0e-20, alpha=0.2)\n",
    "plt.legend(loc=\"lower right\")\n",
    "\n",
    "print(\"Percent errors:\")\n",
    "for i in range(len(sigmas)):\n",
    "    error = abs((proton_chains[i]-sigmas[i])/sigmas[i])*100\n",
    "    print(format(masses[i], 'e'), \":\", error)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mx: 1e-06\n",
      "zth: 456671.82075219066\n",
      "zcrit: 19972.333828000406 \n",
      "\n",
      "mx: 1e-03\n",
      "zth: 25637.637148333786\n",
      "zcrit: 19980.04483011891 \n",
      "\n",
      "mx: 1e+00\n",
      "zth: 19948.280674775753\n",
      "zcrit: 19941.489008883043 \n",
      "\n",
      "mx: 1e+03\n",
      "zth: 19846.40375042482\n",
      "zcrit: 19846.396989420653 \n",
      "\n",
      "[5.91250841e-27 7.12611543e-26 5.59432571e-23 5.63314267e-20]\n",
      "Percent errors:\n",
      "1e-06 : 83.3%\n",
      "1e-03 : 91.9%\n",
      "1e+00 : 90.1%\n",
      "1e+03 : 84.5%\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAEfCAYAAABxkKByAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd3zV9fX48dfJIDshYQkEBARkBVCmAm4RFdSKRdSqqJXW1tFdv9VWWmvV/tQqoiLKcCAuFEex1ao4EJG9QZEhAWSPMHKTe+/5/fH5BC4x415yV5LzfDzySO77fu7nc+7lcs99b1FVjDHGmGhJiHUAxhhj6hdLPMYYY6LKEo8xxpiossRjjDEmqizxGGOMiSpLPMYYY6LKEo+JCBGZIiJ/j3UcdZWIrBCRsyJ07iYiskZEUiNx/ngQifeniLwhIkPCec66yhJPHSYiG0TksIgUicheEflCRH4uIgkBx0wRERWRS8o99lG3fFTUAz8awygR+TxW1w8HEZklIj8tV6Yi0r4m51XVrqo6q0bBVe5OYLKqFofzpCLSU0QWiMgh93fPIB+3WkRurKD8DhGZH84Yg4jlXDeeQyLysYicGHD3A8B90YyntrLEU/cNU9Us4ESc/xh/BCaWO+Zr4PqyGyKSBPwY+DZaQR4vEUmMdQzR5P7bRPL8KTjvhRfDfN4GwFvueXOB54C33PLqPAdcV0H5te595a8lgV+uwkVEGgNvAH8G8oD5wCtl96vqV0C2iPQO97XrGks89YSq7lPVt4ErgetFpFvA3e8AA0Qk1709BFgKfF/VOUWkk4h8ICK73aaZEVUcO1REFgfUvLoH3NfKbabYISK7RGSciHQGxgOnicgBEdnrHjtFRJ4SkZkichA4W0Q6uzWLvW4T1CUB554iIk+IyL/dmt9cETmpkhhTReRFN4a9IjJPRJq59+WJyGQR2SIie0RkhlueKyLvurHvcf/Od++7DxgEjHOfwzgR+dS93BK37MogXp8NIvJHEVkKHBSRJLfsPPf+MSLyqog87z7HFYEffiJyqogscu97TURekcqbmfoBe1W1MODxs0TkXhGZ7Z7jffdDOBRnAUnAo6rqUdWxgADnBPHYF4CBgbUL9/3RHZgWEON9IjIbOAS0C+X9GaTLgRWq+ppbGxwD9BCRTgHHzAIuruF16jxLPPWM+62sEOcDsUwx8DYw0r19HfB8VecRkQzgA+AloClwFfCkiHSt4NhTgUnAz4BGwNPA2yKS4tZY3gU2Am2AlsDLqroK+DkwR1UzVbVhwCmvxmnSyALm4iTO9904bgOmisjJAcdfBfwV55v2WipvDrkeyAFauXH+HDjs3vcCkA50da/zL7c8AZiMU6Ns7R4/DkBV7wI+A251n8OtqnqG+7gebtkrVb0+5Z7DxUBDVfVWEPslwMtAQ5x/y3FwpKbxJjAF51v6NOBHlTx/gAJgTQXlVwM3uM+9AfC7sjvcZFnZz53uYV2BpXrsGl1L3fIquUnwY5waTpnrgJmqujOg7FpgNM77YgfBvz9bV/Mcrg54DksC4jqI0yoQeM5VQI/qnlN9Z4mnftqC8yEU6HngOhHJAc4EZlRzjqHABlWdrKpeVV0ITAeuqODYm4GnVXWuqvpU9TnAA/QH+gItgN+r6kFVLVbV6vp13lLV2arqB3oCmcADqlqiqh/hJLKrAo5/Q1W/cj+wp7qPqUgpzgd/ezfOBaq6X0SaAxcCP1fVPapaqqqfAKjqLlWdrqqHVLUIJ6mdWU38obw+Zcaq6iZVPVzxKfhcVWeqqg8nSZZ9+PXHqWmMdeN+A/iqilgaAkUVlE9W1a/d679KwGuoqg2r+HnAPSwT2FfunPtwkkQwnsNNPOI0o13DD5vZpqjqCvffeQhBvj9V9btqnsNLITyHIpzX0FQhou3FJm61BHYHFqjq5yLSBLgbeFdVD4vIkftFZAXOt3pwPoRPBPqVNYG5knA+9Mo7Ead577aAsgY4CccHbKzkW3xlNgX83QLY5CahMhtxnmOZwCbDQzgfIBV5Aae287KINMTpj7jLLdutqnvKP0BE0nFqP0NwalQAWSKS6CaBYFT1+pTZRNXKP8dUcfqDWgCby9U0qjrXHipOBsG+hpU5AGSXK8um4iRXkTdwaiz9cWqe6cC/yx0T+LxCeX8GK5jnkAXsxVTJajz1jIj0wflQrqhW8SLwWypoZnNHUWW6P5/h/Cf/pNw3w0xVvaWC824C7it3bLqqTnPvay0Vd5pXtnR6YPkWoJUc25ncGthcyWMr5dYI/qqqXYDTcWp117kx5rnJqLzfAicD/VQ1GyhrSivL2sEs/17V63MkvFCfj2sr0FICv0U4ibQyS4GOoVzA7auq7OdP7mErgO7l4ujulldLVQ8Br+P8e1yL0xxbUv6wgL+Dfn+6TW1VPYdrAp5Dj4DHZQAnlXsOnQlojjMVs8RTT4hItogMxekHeFFVl1Vw2FjgfODTCu4r712go4hcKyLJ7k8ft9O3vGeAn4tIP3FkiMjFIpKF0+yzFXjALU8VkQHu47YB+VL1yKe5wEHgD24MZwHD3OcZEhE5W0QK3H6n/ThNbz5V3Qq8h/ONO9e9TlmCycLp19krInnAPeVOuw1oV01ZVa9PTc3BqVXeKs6ghEtxmjcr8xXQUERaVnHMMQK+kFT08w/3sFluHLe7fXu3uuUfwZGh8xuqudRzOINjhlPBaLZygn5/uk1tVT2Hqe6hbwLdRGS4OHOc/oLTb7U64HRn4rxXTBUs8dR974hIEc43wLuAR3A6iX9AVXer6oflmmUq5PZnDMYZkLAFpynmQSClgmPn4/RjjMNpylkLjHLv8+EkivbAdzgDH650H/oRzrfJ70VkJxVwv/VegtP8txN4Eriu3IdBsE7A+Va9H6eT+BOODiu+FicRrQa2A79yyx8F0txrfwn8p9w5HwOuEGfE21i3bAzwnNtxPaKq16em3NfncuAmnCagn+B8KHuqOH6Ke1zYuOe9DKfGshe4EbgsoNbSCphdzWk+xelT2ayq86q5XtDvzxCeww6cpHcfzr9TP44OyClrTTiozgAeUwUJ4jPGGFOHiMhcYLyqTq7k/iY4o/FOqWIwQ7hjeh+4Q53RjLWSiEwHJqrqzFjHEu8s8RhTx4nImThDpHfijAYbD7RzmxCNibpaO6pNRC7DmdfQFHhCVd93O/ueBEqAWQFts8bUZyfjDIHOxJl3coUlHRNLManxiMgknBFD21W1W0D5EJw28UTg2YA5AFWdKxd4SFVvEpFrcWZdvyMir6jqldU93hhjTHTFqsYzBacj9ciwXXck0RM4o6oKgXki8jZOErq/3ONvVNXt7t93u48DyAfKRmsFO4fCGGNMFMUk8ajqpyLSplxxX2Ctqq4DEJGXgUtV9X6c2tEx3PkADwDvubOSwUlY+cBibMSeMcbEpXjq42nJsTOPC3GGK1bmNuA8IEdE2qvqeJzZzeNE5GKc9bt+QERG46znRHpaWq+T2rUJQ+jGGBM9inKgRNlx0EepH7JTEjghM5EGiVL9g2ui9DCy7zs27PGx85Ae98XiKfFU9CQq7YBSZ3XbseXKDlLJHJWAYyYAEwAKunTS6S8+HXqkxhgTI+v3+nhqYTGLtvsY0CiJMWfmcFr+cU9PCt7hvTDz97BrF70nHKjRqeIp8RRy7FIe+TgTv4wxpt4rKlGeX+7hnW9KyGwg/O3MHK4uSCcpIcK1HFVY8x7MHQ8lByEhCT3+JZyA+Eo884AOItIWZ52tkThLsRtjTL3l8yv/WVfK5GUeikqUq7ul85v+WeSlRWEPxD0b4fNHYOsSaH0aDH0UPPvZ9lD/GlUKYpJ4RGQazsZQjUWkELhHVSe66zf9F2ck2yRVDWoBQWOMqYuW7/DyxMJi1u7x06dFA8acmUPXJsmRv7DXA4tegCUvQ4NMuORx6PkTSHDGbG0u0io3iaxOrEa1XVVJ+UzAlpswxtRrOw/5eXaJhw83ltIsI4GxQ3IZ1iGVYxf3jpDC+fD5v2D/Zug+Egb/HTKbhPUS8dTUZowx9VqJT5m+poRpKz34FG7rk8ktvTNJT47C7JBDu+HLJ2Ht/yDvJLjubWgX6p6GwbHEY4wxMaaqfLnFy/hFHrYc8DO4XSp3D8qmdU4UPqLVD6v/DV9NcJrYzrwTBv4aklMjdklLPMYYE0Ob9vt4alEx87b6OCk3iecvbcgZJ0buQ/8Yu9fBZ4/AtuXQZqAzeKBxh4hf1hKPMcbEwMFS5cXlHt78uoS0ZOHPg7K5rnsGyZGeBArgLYYFz8Gy1yAlGy4bDz1GQjT6kLDEY4wxUeVX5YP1pUxc6mFvsTKiSzq/Pz2LxulRGB4N8N1cmP0oFG2FU34C598L6XnRubbLEo8xxkTJ6l0+nlhwmNW7/ZzSLJkxl+bQo1lVO7uH0aFd8MU4WPex05w2fCa0GVD94yLAEo8xxkTY7sN+Ji718P76UpqkJ/DI+Q25rFMaCdFo2vL7YNU7MO8Z8Hnh7LthwO2QFIVldiphiccYYyKk1KfM+KaEF1d4KPHBz3plclufTDIbRGnx/J3fOCsPbF8F7c6Cix+BRidF59pVsMRjjDERMG+rl6cWFrOpyM/ZJ6bw5zNyaJcbpY/c0kOwYAosmw7puXD5s1BwRdQGD1THEo8xxoTR5iI/4xcV8+UWL21yEpk0LI9z2kZpeDTAxi9g9mNwYBv0GgXnjYG03OhdPwiWeIwxJgwOlyovrfQwfU0JyYnCnQOyuKFHJilJUaplHNjuDB7Y8Ck06QQjnofW/aNz7RBZ4jHGmBpQVT7a6OXZJcXsPKxc3imNP56eTbPMKA2P9vtg5QyYN9FZheDce+C0WyEpSqPljoMlHmOMOU7f7PbxxMJiVuz00a1JMk8PzaFX8yh+4O9Y46w8sHMNtD8PLnoI8tpG7/rHyRKPMcaEaG+xnynLPMz8tpTctAQePDeHH3dJj87waICSQzB/Iqx4EzIawxWToeuP4mbwQHUs8RhjTJB8fuWdtaU8t9zD4VLlxp4Z3N4vi5yUKA2PVoUNn8MXY+HgTuhzE5zzZ0hrGJ3rh4klHmOMCcLC7708taiYDfv8DGjVgDFn5NChURQ2ZStzYBt8/hh89wU06wpXvQz5vaN3/TCyxGOMMVX4/oCfCYuL+azQS35WIk9fnMvgdlHalA3A74Xl02H+FBCctdX63wKJUUx6YRZU4hGRYFaQ86vq3hrGY4wxcaHYq7yyysOrq0tIFOG3/bO4+dRMUqM1PBqcFQc+ewR2fQMdL3AGDzRsHb3rR0iwNZ4t7k9Vr3giUPtfEWNMvaaqfLrJy4TFxWw/pAzrmMb/DcimRVaUhkcDlByAr56FlW9B1gkw4gXoPKzWDB6oTrCJZ5WqnlLVASKyKAzxGGNMzKzf6wyPXrLdR+fGSTx+YQ79WkZxMU1VWP+JMxH08G7o9zM4+y5IzY5eDFEQbOI5LUzHGGNM3NnvUZ5bXsy7a0vJbiD8/ewcruqaTmJCFGsY+7c6++RsmgvNe8A1r0HLU6N3/SiqNvGIyPnACBF5QlUXi8hoVZ1Q/jhVLY5IhMYYEyE+vzJzXSlTlnk4UKL8pCCd3/TPpmFqlIZHgzN4YOmrsPA5SEiCIQ9An5shse6O/Qrmmf0CuAG42x1k0DOyIRljTOQt2+7liYXFfLvXT/+WDbjnzBw6N47ySLHvlzvbFuxeB52GwoX/hJyW0Y0hBoJJPDvc0Wq/E5EHgD4RjskYYyJm+0E/zywpZtZ3XppnJvDEhblc1D6Kw6MBPEXw1QRng7bsljByGnS6KHrXj7FgEs+/y/5Q1TtF5LYIxmOMMRFR4lNeW13Cyys9KHB730xu6ZVJWnIUm9VU4duPYM44KN7vLOZ51v9BSmb0YogD1SYeVX2r3O3HIxeOMcaEl6ryxWYvTy8qZutB5cKTUvnToGxaZUe5D2X/ZvjsX7B5PrToBcMehebdoxtDnAjplReR3sBdwInuYwVQVa2fr54xJq5t3OfjqYXFLNjmo0NeElMH5zCgVRSHRwP4SmHJy7DoRWe1gYsegt43QkIU5wXFmVBT/lTg98AywB/+cIwxpuYOligvrPAw4+sS0pOFe87I5icFGSQnRnkC5tal8NnDsHcjdLnMGbGW3Ty6McShUBPPDlV9OyKRGGNMDflVeX99KZOWethbrIzsls7v+mfRKD3KtYvifTD3aVgzE3JawdWvQcfB0Y0hjoWaeO4RkWeBDwFPWaGqvhHWqIwxJkQrdzrDo7/e7afXCcn89awcujWN8i6cqvDN+/DlU87ItQF3wJl/hAYZ0Y0jzoWaeG4AOgHJHG1qU8ASjzEmJnYd9jNxiYcPNpTSLCOBRwc35NKT06I7PBpg7yb4/F+wZSG07OMMHjihW3RjqCVCTTw9VLUgIpEYY0wISnzKm1+XMHWFB68fftE7k1/2ziSjQRSHRwP4SmDxS7BoKiSnwdB/wamjICHKcdQioSaeL0Wki6qujEg0xhgThLlbShm/yENhkZ/z2qZw96Ac2jSMwRIzWxY52xbs2wTdroAL/gFZzaIfRy0T6r/UQOB6EVmP08djw6mNMVFTWORj/CIPc7d4adswkcmX5HF2m9ToB3J4r9OP881/IbcN/OQNaH9u9OOopUJNPEMiEsVxEJHLgIuBpsATqvp+RWWxjNEYEx6HSpWXVnqYvqaElEThroHZXN8jgwbRHh6tCmveg7njofQwDPodnPE7p4nNBC2kxKOqG8NxURGZBAwFtqtqt4DyIcBjOJvKPauqD1QRywxghojkAg8B71dUFo54jTGx4Vflww2lTFzqYddh5YrOafzh9GyaZsRg8uWejc6CnluXQOvTYOij0LRT9OOoA0JdueA54I6yLa7dD/iHVfXGEK87BRgHPB9w7kTgCeB8oBCYJyJv4ySh+8s9/kZV3e7+fbf7uEAVlRljapE1u308uaCYlbt8dG+azLPDcjjlhCgPjwbwemDRC87qAw0y4ZLHoedPbPBADYTa1Na9LOkAqOoeEalyZ9KKqOqnItKmXHFfYK2qrgMQkZeBS1X1fpza0THEGSv5APCeqi6srMwYU7vsKfYzaamH/64rJS8tgf93XkOGd04jIRbbPhfOd4ZI798M3UfC4L9DZpPox1HHhJp4EkQkV1X3ALj784RrKElLYFPA7UKgXxXH3wacB+SISHtVHV9J2TFEZDQwGqBFcxt9Yky88PqVt78p4YXlHop98NNTMritbxbZKTGoWRzaDV8+CWv/B3knwXVvQ7szox9HHRVq0ngY+EJEXseZODoCuC9MsVT0dUYrO1hVxwJjqyur4HETgAkABV06VXp+Y0z0LPjey5MLi/luv58zWqfwlzOyaZ8X5U3ZANQPq//t7JXj9cCZd8LAX0NyDEbO1WGhDi54XkTmA+fgJIrLwzinpxBoFXA7H9gSpnMbY+LQ1gN+nl5UzOzNXlplJ/Ls0DzObZsS/VUHwNkF9LNHYNtyaDPQGTzQuEP046gHQm4mcxNNJCaQzgM6iEhbYDMwErg6AtcxxsTQyp1e5n/vZdsBP7O+85KYIPz+tCxuOiWT1KQYJBxvMSx8Hpa+CinZcNl46DESYpH86omgEo+ILFTVU2t6TMCx04CzgMYiUgjco6oTReRW4L84I9kmqeqKYM5njKkdVuzw8ruPD+F1V3oc1LoB/zw3l+ZZMdqb5ru5MPsxKNrijFQ7/2+Q0Sg2sdQjwdZ4OovI0iruFyAn2Iuq6lWVlM8EZgZ7HmNM7bFyp5f75hw+knQSgP4tU2KTdA7tgi/GwbqPnea04f92mtdMVASbeIKZJeWrSSDGmLrp+wN+Ji4tZtZ3XnJThOQE8CskJ0L//CjvBur3wap3YN4z4PPC2XfDgNshKcpx1HNBJZ5wrVhgjKk/DpYqL6/08MaaEhIS4Pa+mfzs1ExW7/LyZaGH/vkp9GoexQmhO79xVh7YvgranQUXPwKNTore9c0RMVjO1RhTl/n8yn/WlfLccg97ipXLO6Xx+9OyjzSp9WreILoJp/QQLJgCy6ZDei5c/gwU/NgGD8SQJR5jTNgs+N7L04uKWb/PT58WDfjzoGy6N4vBMjdlNn7hDB44sA1OvR7OGwPpebGLxwA1TDwi8gSQoaqjRGSwrQZtTP20cZ+PCYs9fLXVmY/z1EW5DDkpNTbzcQAObHcGD2z4FJp0ghHPQ+v+sYnF/EBNazwlwDb373Ow1aCNqVf2efw8v9zDv9eWkpYs/GlgNtd3zyAlFvNxwBk8sHIGzJvorEJw7j1w2q2QFMNal/mBmiaeQzjroiUDrcMQjzGmFijxKW99U8JLKzwc9sLV3dL5Vb8sGqXHaD4OwI41zsoDO9dA+/Pgoocgr23s4jGVqmniuQdnwc0ngJdqHo4xJp6pKp8XenlmiYetB/ycfWIKfxqYTYdGMVhXrUzJIZg/EVa8CRmN4YrJ0PVHNnggjoW6H889qvrXstuq6gWeDHtUxpi4s2a3j6cXFbNsh48OeUk8f2lDzjgxhotnqsKGz+GLsXBwJ/S5Cc75M6Q1jF1MJiih1njuEZF0IA9YCLxctkWCMaZu2nHIz+SlHj7YUEqjtAT+cU4OI7qkk5QQwxrFgW3w+WPw3RfQrCtc9TLk945dPCYkoSYeBYpx1lM7FWeLhJGquiTskRljYuqwV3l1lYfXVpegwC29MvlF70yyYrE/Thm/F5ZPh/lTnIW6zr8X+t8CiTFs6jMhCzXxrFbVe9y/XxeRKcB4nBFtxpg6wK/KB+tLmbzMw67DytAOqfxxQDatsmM87W/7KmfwwK5voOMFzuCBhjamqTYK9Z20U0R6qeoCAFX9WkRsH1hj6ogl27yMX1zM2j1+ejRLZsLQnOiuMlCRkgPw1bOw8i3IOgFGvACdh9nggVos1MRzO/CyiCwAlgHdgfVhj8oYE1WFRT6eXexh9mYvzTMTeOyChgzrmEZCLD/cVWH9J85E0MO7od/P4Oy7IDU7djGZsAh1B9IlItITOA/oBnwMTItEYMaYyCsqUaau8PDWNyU0SIzxhmyB9m+F2Y/CprnQvAdc8xq0DGq7L1MLHM8OpB7g3+6PMaYW8vqVd9eW8sJyD0UlypVd0/lN/yyaZsRwAig4gweWvgoLn4OEJLjgfug7GhJtWcm6xP41jalHVJW5W7xMWOxhU5GfAfkNuGtQDl2axMGosG0r4LOHYfc66DQULnwQcvJjHZWJAEs8xtQT3+7x8fTiYhZt89GuYRIThzXknDYpsVvIs4ynCL6a4GzQlt0SRr4EnS6ObUwmomq6OnVzYLfb/GaMiUO7D/uZsszDf9aVkpMijDkzm2u6ZZCcGOOEowrffgRzxkHxPmcxz7P+D1IyYxuXibia1nheAE4Skemq+rtwBGSMCQ+PV5m+poSXV3ko9cONPTO4vW8WOakxnAAKTpPa+k9gy1LYuRpanArDHnUGEZh6oUaJR1XPE6ee3iVM8RhjakhV+Xijl4lLi9l+SLngpFTuHJBN24Zx0LK+ZRH8+3egPud2/1/C4HshIcaDGkxUhbpIaAowHGgT+FhV/Vt4wzLGHI8VO72MX+Rh9S4fXZsk89iQbE7LT4l1WEfn5Hz60NGkI4mQ0ciSTj0U6legt4B9wALA+nWMiRPfH/Dz7JJiPtnkpWl6Av/vvIYM7xzjCaBldn0LXzwOWxdDbjvwlTgbtiU2gDaDYh2diYFQE0++qg6JSCTGmJAdLFWmrfTw5poSEhLgjr6Z/KxXJunJMe7HASjeC/Mmwep3IbUhXPwI9BoFmxfAhs+cpNOqb6yjNDEQauL5QkQKVHVZRKIxxgTF51feW1fKc8s87PUol3dK4/enZdM8Kw6arfxeWPk2LJjsbNLWdzScdSek5Tr3t+prCaeeCzXxDARGich6nKY2AVRVu4c9MmNMheZv9TJhcTHr9/np06IBfx6UTfdmMV7Is0zhfGd49J4N0O4sGPIANO0c46BMvAk18VwYkSiMMdXauM/HhMUevtrqpVV2IuMvyuWCk1JjPwEUYP9mmPMkbJwNuW2cSaAnX2QrSJsKhbpI6MZIBWKMqdjeYj/PL/cw89tS0pOFuwZmc133DFJivZAnQOkhWDTVWV8tMRnO/YszRDo5hltim7gX8sB+EekBlA1F+cx2HzUmMkp8ylvflPDSCg+HvXBNQTp39M2iUXoc9OOoH775n7PUzaGd0H0knDcGspvHOjJTC4Q6j+cO4GbgDbfoRRGZoKqPhz0yY+opVeWzQi/PLi5m60Hl7BNTuGtQNu3z4mAhT4Dtq53h0dtXOKsOXPUytOoT66hMLRJqjecmoJ+qHgQQkQeBOYAlHmPCYM0uH+MXFbN8p4+OeUk8eF42Z5wYJ81Wh3bBV8/A1/+BjCZw6ZPQ4ypIiIOh26ZWCTXxCOALuO1zy4wxNbDjkJ+JSzx8uLGURmkJ/OOcHEZ0SScpIQ7+e/lKYPkbsPB58JXCgDtg0O9sJ1Bz3EJNPJOBuSLyJk7CuQyYFPaojKknDpcqr6z28PrqEhT4Re9MbumVSVZKHNQiVOG7Oc5otf2F0PECZ2O2RifFOjJTy4U6qu0REZkFDMBJPNer6uJIBFYdEbkMuBhoCjyhqu+75RnAp8A9qvpuLGIzpjp+Vd5fX8qUZR52HVaGdUzjD6dn0So7DhbyBNi7EeY8AZu+gsYd4Jrp0OG8WEdl6oig3uUi8rmqDhSRIkAJaF4TEVXVkOrcIjIJGApsV9VuAeVDgMeAROBZVX2gsnOo6gxghojkAg8B77t3/RF4NZR4jImmJdu8jF9czNo9fno2S2bC0Bx6NY+TCaCeIqdJbcWbkJwGF/zD3Xo6TgY2mDohqMSjqgPd31lhuu4UYBzwfFmBiCQCTwDnA4XAPBF5GycJ3V/u8Teq6nb377vdxyEi5wErgTjpjTXmqMIiH88u9jB7s5cWWYk8dkFDLumYFh8TQP0+WPMezJvorLF26nVwzp8hs0msIzN1UKjDqR9U1T9WV1YdVf1URNqUK+4LrFXVde55XwYuVdX7cWpH5WMR4AHgPVVd6BafDWTg7A90WERmqqo/lNiMCbeiEuXF5R7e+qaElCTh96dlcdMpmaTGw+r+Se8AACAASURBVARQgO+XOsOjd34DrU9zlrlp0TPWUZk6LNQG5fNxmrICXVhB2fFoCWwKuF0I9Kvi+NuA84AcEWmvquNV9S4AERkF7Kwo6YjIaGA0QIvmzcIQtjEV8/qVd9aW8OLyEg6UKld2SefX/bNomhEHE0ABDmyHueOd7aezW8DwidBtuC1zYyIu2D6eW4BfAO1EZGnAXVnA7DDFUtG7XSs7WFXHAmMruW9KFY+bAEwAKOjSqdLzG3O8VJUvt3iZsNhDYZGfAfkNuPuMHDo3jpN+Eq8Hlr4Ci19ybp/xBxj4K2iQEdu4TL0RbI3nJeA9nL6WOwPKi1R1d5hiKQRaBdzOB7aE6dzGRMXaPT4mLCpm0XYf7RomMWlYQ85ukxIf/Thlu4B+OR4OfA9dLoXz74XcE2Mdmalngh1csA9n59GrIhjLPKCDiLQFNgMjgasjeD1jwmbXYT9Tlnn477pSclKEv56Zw9Xd0klOjIOEA8fuAtq0Cwx/BtqeEeuoTD0V6uCC54A7VHWvezsXeFhVbwzxPNOAs4DGIlKIM+dmoojcCvwXZyTbJFVdEcp5jYk2j1d5fU0Jr6zyUOqHm07J4LY+WeSkxsEEUIDifTB/Eqx6B1Jz4OKH4dRRkBgn84VMvRTqu697WdIBUNU9InJKqBdV1QprTqo6E5gZ6vmMiTa/KrM2epm4tJjth5QhJ6Vy54Bs2jSMkw/0I7uAToGSg9DnZmcX0PS8WEdmTMiJJ0FEclV1D4CI5B3HOYyp1Zbv8PL0omJW7/bTrUkyY4dk0z8/JdZhHbV5AXwxDvash7ZnOsOjm3WJdVTGHBFq0ngYmCMir+GMOBsB3Bf2qIyJQ1sP+Jm4pJhPNnlplpHAQ+c35PJOaSTEw8ABgP1b4MunYMNn0PBEuHIqdLrYhkebuBPqWm3Pi8h84Byc4c+Xq+rKiERmTJw4WKK8tNLDm1+XkJgAd/TN5Ge9MklPjpN+nNJDsOglWPYqJCQ5Kw6cdqvtAmriVqiDCwQ4FchT1b+JSGsR6auqX0UmPGNix+dXZq4r5fllHvZ6lOGd0/j9admckBknE0BVYe0Hzi6gB3dC9yvdXUBbxDoyY6oUalPbk4Afp8bzN6AImA7Y9oOmTpm31cuExcVs2Oenb4sG/PmMbAqaxslCngA73F1At62AFqfAyGnQqm+sozImKKEmnn6qeqqILIIjo9ri6H+jMTWzcZ+PpxcXM2+rj9Y5iYy/OJcL2qXGxwRQcHcBfdbdBbQxXPoE9LjadgE1tUqoiafUXUVaAUSkCU4NyJhabW+xn+eXe5j5bSnpycLdA7O5tnsGKfGykKevFJZPP7oL6Om3Okvd2C6gphYKNfGMBd4EmorIfcAVONsSGFMrlfiUGV+X8NJKD8Ve+ElBBnf0yyQvLU76ccDdBfQJ2FcIHS5w9shp3D7WURlz3EId1TZVRBYA57pFl6nqqvCHZUxkqSqfFXp5dnExWw8q57RJ4U8Ds2mfFycLeYK7C+iTsGkuNGoP17wOHc6PdVTG1Fiwq1P3ATap6vequlpE+gHDgZNFZEwYFwo1JuLW7PIxflExy3f6OLlREv88P5tBreNo6HHJAVjwPKyYDsnpMPg+ZxfQJOtONXVDsDWep3H2vkFEzsBZpfo2oCfOFgNXRCQ6Y8Jo+0E/k5Z6+HBjKY3SErj/nBxGdEknMSFO+nH8PmfQwLxn4fBeOPVadxfQprGOzJiwCjbxJAbUaq4EJqjqdGC6iCyOTGjGhMfhUuWV1R5eX12CAr/oncktvTLJSomjkWDfL3N3Af0aWvWHCx9whkkbUwcFnXhEJElVvTj9O6OP4xzGRJXPr3ywoZTJSz3sLlaGdUzjj6dnkZ8dR2/ZA9th7tPw7YeQ1dx2ATX1QrD/A6cBn4jITuAw8BmAiLTH2afHmLiyeJuX8YuK+Xavn1OaJfPMsBx6NY+jPpLAXUBVbRdQU68EuxHcfSLyIdAceF9Vy7aMTsDp6zEmLhTu9zFhiYc5m720yEpk7JBchnWIowmgqrD+U5g7Hoq2QudLYPDfbRdQU68E3eagql9WUPZ1eMMx5vjs9yhTV3h465sSUpOEP5yexY09M0mNlwmgALvXOf04WxZB087wo6eh3ZmxjsqYqIujxm5jQuf1K29/U8KLK0o4WKpc2SWd35yWRZP0OJoAWrwP5k+GVW87u4Be9BD0usF2ATX1lr3zTa20YoeXd74tYel2HzsOKQNaNeDuQTl0bhxHE0B/sAvoT+Gs/7NdQE29Z4nH1BqHSpVvdvv4ZFMp764tRXE2hfrTgCxuPjUzfvpxADYvhDnjnOa1tmfAkAdtF1BjXKHux/Nj4D+qWiQid+PszfN3VV0YkehMveX1K+v3+lm928eaXT7W7PaxcZ8fLXdcgkCpn/hJOvu3wpdPuruAtoYrX4ROQ214tDEBQq3x/FlVXxORgcAFwEPAU0C/sEdm6g1VZfMBP2t2HU003+7xUeKue56XlkDPZsn8qFMDejRzmtJ+PnM3pT5IToT++SkxjN5VehgWT4Wlr0JCIpxzN5x2m+0CakwFQk08Pvf3xcBTqvqWiIwJb0imrtt92M8aN8Gs3u3j691+ikqcukxaklDQNJlRPVPp0cxJNC2zEn9Qo5n6o8Z8Weihf35KbOfnqMLa/7m7gO6AghFw/l9tF1BjqhBq4tksImXrtj0oIik4c3mMqVBZv8ya3W6S2eVj2yEnySQKnNwoiaEd0+jZLJkezRrQPi+JpCDWTuvVvEHsJ4TuWA1fjINty6F5T7hyKrS2yr8x1Qk18YwAhgAPqepeETkB+H34wzK1UXX9Mq1zEumbn0qPZsn0bJZM1ybJpCXXwu8th3Y7C3mueQ8yGsEl46DnNbYLqDFBCjXxXEwFgwvCH5aJd6rKlgN6pLlszS4fa/f6KHEbY8v3y3Rvlhxfm6sdj2N2AS2B034JZ/7BmZtjjAmaDS4wQdlT7Gf1rsr7Zbo1Teb67m6/zAnJ5FfQL1OrffeluwvoJugw2N0FtEOsozKmVrLBBeYHDpcqX+8JSDKV9Mv0aJpMjxMa0CHIfplaae937i6gX0LeSXD1a9BxcKyjMqZWs8EF9ZzXr2zYd7Q2s2a3j437/fjdjplW2XWkXyZUJQecJrXlZbuA/h36/sx2ATUmDGo6uKA5Nrig1qiuXyY3Veh5QgMudftletSFfplQqd8ZNFC2C+gpP4Fz/2K7gBoTRiElHlU9JCLfAheIyAXAZ6r6fmRCMzW1p9h/TJJZU0m/TPdmDejZLJn87DrWLxOq75e7u4Cugfy+cOGD0PLUWEdlTJ0T6pI5dwA3A2+4RS+KyARVfTzskZmQHC5Vvtnjzpdxm8y2HTzaL9OxURJDO6Q5NZm63i8TqgPbnQmga//n7AJ6+bNQcIUtc2NMhITa1HYT0E9VDwKIyIPAHMASTxRV1y+Tn51InxYp9HRHmHVtkkx6feiXCZXX4yxxs3iqswLBoN/BwF9DSmasIzOmTgs18QhHR7bh/m1fCyNIVdl6QI+ZlPnNnsr7Zbo3TaZRPO1FE49UnUU8v3zK3QV0mLsLaJtYR2ZMvRBq4pkMzBWRN93blwETwxtS/RbYL/P1bh9rdvnZ7/bLpCZBtyYNuM6dL2P9Msdh9zpnmZstC20XUGNiJOjEI86n22vALGAgTk3nBlVdFJnQqo3nMpz5RE2BJ1T1fRFJAO4FsoH5qvpcLGILVlX9MgkCHfOSuKisX6ZZAzo2sn6Z41a8H+ZPcnYBTcm2XUCNiaGg/9epqorIDFXtBdRo/x0RmQQMBbarareA8iHAY0Ai8KyqPlBFPDOAGSKSi7OCwvvApUBLYDdQWJMYw62sXyZwlJn1y0SB3wur3nWSTskB6H0TnP0n2wXUmBgK9evelyLSR1Xn1fC6U4BxwPNlBSKSCDwBnI+TNOaJyNs4Sej+co+/UVW3u3/f7T4O4GRgjqo+LSKvAx/WMM7jUlG/zNo9Pjxuv0zDVKFnswZccnIyPU9oYP0ykbJlkTM8evc6aDPIGR7drGusozKm3gs18ZwN/ExENgIHcZrbVFW7h3ISVf1URNqUK+4LrFXVdQAi8jJwqarej1M7Oobb9PcA8F7ADqiFQIn7t6/8YyJlT7GfrwOayyrql/lJQSo9TnD6ZVpZv0xk7d/qDBzY8KmzC+iIF5wBBPaaGxMXQk08F0YkCkdLYFPA7UKqXnz0Npyle3JEpL2qjseZX/S4iAwCPq3oQSIyGhgN0KJ5s5CDPOxV1u72HVOb+d76ZeJD6WFY/BIsfcXZBfTsu+H0WyE5LdaRGWMCBJV4RKQ90ExVZ5crHwRsCVMsFX06awVlzh2qY4Gx5coO4cw1qpSqTgAmABR06VTp+QF8fmV9Vf0yWYn0KuuXaZZMt6bWLxMTqvDthzD3aXcX0B/DeX+FnJaxjswYU4FgazyPAn+qoPywe9+wMMRSCLQKuJ1P+JJatVSV7w/qkeay1bsq75fp0awB3Zsl09j6ZWJvxxp3F9Bl0LwHXPkitO4f66iMMVUINvG0UdWl5QtVdX4FfTXHax7QQUTaApuBkcDVYTp3hb4/6OfhuYfYVQxf7/axz+NUZVISoaCp9cvEte/mwPzJsPNryGgMlzzu7gJqXwaMiXfBJp7UKu4LuQFdRKYBZwGNRaQQuEdVJ4rIrcB/cUayTVLVFaGeOxT7S5T/rPfSOjuRISel0eMEt18mL4nkREsycaHkAOwrdH72bnI2Ytv5jfMbICEJLn8GTjontnEaY4IWbOKZJyI3q+ozgYUichOwINSLqupVlZTPBGaGer6aSBS4sms6v+yTFc3LmkC+Emck2r5NbpLZ5CaZQji8O+BAcUapJTbAHVDp9O9sWWSJx5haJNjE8yvgTRG5hqOJpjfQAPhRJAKLluRE6J+fEusw6j71w8GdAUll09GaTNFW5/4y6Y2gcUfodBE0an/0J7ctJKfCpq/guUuchJXYwJmjY4ypNUS1yoFdxx4scjZQttLAClX9KCJRRUnr9p30zVeeo1dz21UybIr3H1tzKWsi21/orAZdJjnt2KTSqIP7ux2k5VZ/nU1fOQt9thkErfpG7vkYY35ARBaoau/jfnwoiaeu6d2ji85/b2qsw6h9vB7Yv/loc1hg7aV479HjJBFyW0Ojjm5SOeloosluYRM6jamlapp4bIVEUzG/z9kgrXy/y/5CKNrGMVOsMps5yaRlL+d3Y7f20vBESLLapDHmWJZ46jNVKN5Xcb/L/s1OH0qZBplOMmkzKKB5zK3BpNjADGNM8Czx1Aelh48mlMCay75C8BQdPS4h2dkMrWln6Dw0oN+lPWQ2taYxY0xYhJR4RCQFGA60CXysqv4tvGGZkPm9UPS9m1jK9bsc3H7ssVktnOaw1qcfTSyN20NOa9ufxhgTcaF+yrwF7MMZUu2p5lgTbqrOvJa9FfS77N/s9MuUSc1xaiwnnXM0sTRqD3ntoEFG7J6DMabeCzXx5KvqkIhEYo4qOfjDiZRlt0sPHT0uMQXy2jprlHUbfuzQ5PQ8axozxsSlUBPPFyJSoKrLIhJNfeIrhaItAc1iAQnm0K6AAwVyWjk1lnZnuv0ubqd+Tr6tTWaMqXVCTTwDgVEish6nqe24NoKrNwJn6x9Zayxwtn5A01h6IyepdLzwaGJp3OHobH1jjKkj4mkjuNrLU1SuUz8gwXiLjx6XlOoklPzex853yWvnNI0ZY0w9EFLiUdWNItIDKFsc6zNVXRL+sOKQ1wP7txxNLHsD+l0qnK3fAToMDpit3wGymkOCbRRnjKnfQh1OfQdwM84W0wAvisgEVX087JFFw6HdsG0FNOvq3Pb7nKHHgUvwH2ka+55jZ+s3dWfrX3LsfJfcNjZb3xhjqhDqIqFLgdNU9aB7OwOYU1v7eHq3SNL5P8uGZgXg2VfBbP2MY5NK4Gz91OzYBW6MMTEU7bXaBAjoEcfnltVS6nTw7/8O8vtBp4uP9rs0au+sQWZDko0xJqxCTTyTgbki8iZOwrkMmBT2qKJGICkNrnrFltY3xpgoCXVwwSMiMgsYgJN4rlfVxZEILCqym8P10y3pGGNMFAWVeETkc1UdKCJFOD3sEnCfqmrt7PDIbGZJxxhjoiyoxKOqA93ftv69McaYGglpUomIPBhMmTHGGFOZUGcznl9Bma1mYIwxJmjB9vHcAvwCaOfO5SmTBcyORGDGGGPqpmBHtb0EvAfcD9wZUF6kqrvDHpUxxpg6K9jBBftwNoC7KrLhGGOMqetCHVzwnIg0DLidKyK1eAKpMcaYaAt1cEF3VT2yFLOq7gFOCW9Ixhhj6rJQE0+CiOSW3RCRPEJfdscYY0w9FmrSeBhn++vX3ds/Bu4Lb0jGGGPqslDXanteROYD5+Asm3O5qq6MSGTGGGPqpOPZDnMr8BWwBGgsImeENyRjjDF1Wag7kP4UuAPIBxYD/YE5ODUgY4wxplqh1njuAPoAG1X1bJwRbTvCHpUxxpg6K9TEU6yqxQAikqKqq4GTwx9W9UTkMhF5RkTeEpHBbllrEXlbRCaJyJ3VncMYY0z0hZp4Ct0JpDOAD0TkLWBLqBd1E8N2EVlernyIiKwRkbXVJQ5VnaGqNwOjgCvd4o7Av1X1RqBLqHEZY4yJvKD7eEREgNvdCaRjRORjIAf4z3FcdwowDng+4PyJwBM4K2AXAvNE5G0gEWeNuEA3qup29++73ccBLALuEpErgReOIy5jjDERFnTiUVUVkRlAL/f2J8d7UVX9VETalCvuC6xV1XUAIvIycKmq3g8MLX8ONxE+ALynqgvd4huAe9zzvw5MPt4YjTHGREaoTW1fikifiEQCLYFNAbcL3bLK3AacB1whIj93y/4D3C4i44ENFT1IREaLyHwRmb9jh42LMMaYaAt15YKzgZ+LyAbgIM4kUlXV7mGIRSoo08oOVtWxwNhyZcuBK6q6iKpOACYA9O7du9LzG2OMiYxQE08kdxstBFoF3M7nOAYuGGOMiW+hNrV9BwwCrlfVjTg1kmZhimUe0EFE2opIA2Ak8HaYzm2MMSZOhJp4ngRO4+iGcEUcHVEWNBGZhrPiwckiUigiN6mqF7gV+C+wCnhVVVeEem5jjDHxLdSmtn6qeqqILAJnPx63dhISVa1wJ1NVnQnMDPV8xhhjao9QE0+pO99GAUSkCeAPe1TGGBOC0tJSCgsLKS4ujnUodUpqair5+fkkJyeH9byhJp6xwJtAUxG5D2cE2Z/DGpExxoSosLCQrKws2rRpgzPFz9SUqrJr1y4KCwtp27ZtWM8d6n48U0VkAXAuzvDny1R1VVgjMsaYEBUXF1vSCTMRoVGjRkRivmOo2yI8qKp/BFZXUGaMMTFjSSf8IvWahjqq7fwKyiI5t8cYY2qNN998ExFh9erV1R9ciVGjRvH6669Xecw//vGPY26ffvrpx3WtMWPG8NBDDx3XY2siqMQjIreIyDKc4c9LA37WA0sjG6IxxtQO06ZNY+DAgbz88ssRvU75xPPFF19E9HrhFmyN5yVgGM6EzmEBP71U9ScRis0YYyJmwcY9PPHxWhZs3BOW8x04cIDZs2czceLEI4ln1qxZnHXWWVxxxRV06tSJa665BlVnpa6//e1v9OnTh27dujF69Ogj5WU+/PBDfvSjHx25/cEHH3D55Zdz5513cvjwYXr27Mk111wDQGZm5pHj/vnPf1JQUECPHj24805nd5lnnnmGPn360KNHD4YPH86hQ4fC8pyPV1B9PKq6D9jH0YmjxhgTl/76zgpWbtlf5TFFxaWs/r4Iv0KCQKcTsshKrXzIcJcW2dwzrGuV55wxYwZDhgyhY8eO5OXlsXChs2j+okWLWLFiBS1atGDAgAHMnj2bgQMHcuutt/KXv/wFgGuvvZZ3332XYcOGHTnfOeecwy9/+Ut27NhBkyZNmDx5MjfccAPDhg1j3LhxLF68+AcxvPfee8yYMYO5c+eSnp7O7t27Abj88su5+eabAbj77ruZOHEit912W5XPJ5JC6uMRkRQRuVpE/iQifyn7iVRwxhgTCfuLvfjdCoZfnds1NW3aNEaOHAnAyJEjmTZtGgB9+/YlPz+fhIQEevbsyYYNGwD4+OOP6devHwUFBXz00UesWHHsQi0iwrXXXsuLL77I3r17mTNnDhdeWHWX+v/+9z9uuOEG0tPTAcjLywNg+fLlDBo0iIKCAqZOnfqDa0VbqPN43sKp+SwAPOEPxxhjaqa6mgk4zWzXPPslpV4/yUkJPDbyFHqdmHvc19y1axcfffQRy5cvR0Tw+XyICBdddBEpKSlHjktMTMTr9VJcXMwvfvEL5s+fT6tWrRgzZkyFk1/Lajipqan8+Mc/Jimp6o9sVa1wJNqoUaOYMWMGPXr0YMqUKcyaNeu4n2s4hDqqLV9Vr1TVf6rqw2U/EYnMGGMipNeJuUz9aX9+M/hkpv60f42SDsDrr7/Oddddx8aNG9mwYQObNm2ibdu2fP755xUeX5ZkGjduzIEDByodxdaiRQtatGjB3//+d0aNGnWkPDk5mdLS0h8cP3jwYCZNmnSkD6esqa2oqIjmzZtTWlrK1KlTa/JUwyLUxPOFiBREJBJjjImiXifm8suz29c46YDTzBY4EABg+PDhvPTSSxUe37BhQ26++WYKCgq47LLL6NOn8v01r7nmGlq1akWXLl2OlI0ePZru3bsfGVxQZsiQIVxyySX07t2bnj17Hhkqfe+999KvXz/OP/98OnXqdLxPM2yk/EiKKg8WWQl0ANbhNLWFcyO4qOvdu7fOnz8/1mEYY2po1apVdO7cOdZhRMStt97KKaecwk033RST61f02orIAlXtfbznDLWPZwhusjneCxpjjAlOr169yMjI4OGH61aPRlCJR0SKqDjZlCWh7HAGZYwxBhYsWBDrECIi2Hk8WZEOxBhjTP0Q6uACY4wxpkYs8RhjjIkqSzzGGGOiyhKPMcaEgYjw29/+9sjthx56iDFjxkQ1hvnz53P77bcf12PPOussojW9xBKPMcaEQUpKCm+88QY7d+6MyfW9Xi+9e/dm7NixMbl+KCzxGGPqp01fwWcPO7/DICkpidGjR/Ovf/3rB/eV39ytbBuDWbNmceaZZzJixAg6duzInXfeydSpU+nbty8FBQV8++23AOzYsYPhw4fTp08f+vTpw+zZswFnI7fRo0czePBgrrvuOmbNmsXQoUMBZ5uGG264gYKCArp378706dMBuOWWW+jduzddu3blnnvuCctzD1WoE0iNMSa+vXcnfL+s6mM8+2HbclA/SAI06wYpVUxHPKEALnyg2kv/8pe/pHv37vzhD38IOtwlS5awatUq8vLyaNeuHT/96U/56quveOyxx3j88cd59NFHueOOO/j1r3/NwIED+e6777jgggtYtWoV4Mz1+fzzz0lLSztm8c97772XnJwcli1zXos9e5x9h+677z7y8vLw+Xyce+65LF26lO7do7v4jCUeY0z9U7zPSTrg/C7eV3XiCVJ2djbXXXcdY8eOJS0tLajH9OnTh+bNmwNw0kknMXjwYAAKCgr4+OOPAWe7g5UrVx55zP79+ykqKgLgkksuqfBa//vf/47ZCTU311mT7tVXX2XChAl4vV62bt3KypUrLfEYY0yNBFEzYdNX8Nwl4CuBxAYw/Flo1Tcsl//Vr37Fqaeeyg033HCkLCkpCb/fSXSqSklJyZH7ArdNSEhIOHI7ISEBr9fZJ8jv9zNnzpwKE0xGRkaFcVS0RcL69et56KGHmDdvHrm5uYwaNarC7Rgizfp4jDH1T6u+cP3bcM5dzu8wJR1wNl8bMWIEEydOPFLWpk2bI8vfvPXWWxVuaVCVwYMHM27cuCO3K9p9tLrH7Nmzh/3795ORkUFOTg7btm3jvffeCymOcLHEY4ypn1r1hUG/DWvSKfPb3/72mNFtN998M5988gl9+/Zl7ty5ldZSKjN27Fjmz59P9+7d6dKlC+PHj6/2MXfffTd79uyhW7du9OjRg48//pgePXpwyimn0LVrV2688UYGDBgQ8nMLh5C2RahrbFsEY+qGurwtQqxFYlsEq/EYY4yJKks8xhhjosoSjzHGmKiyxGOMqRPqc391pETqNbXEY4yp9VJTU9m1a5clnzBSVXbt2kVqamrYz20TSI0xtV5+fj6FhYXs2LEj1qHUKampqeTn54f9vLU28YhIZ+AOoDHwoao+JSIZwJNACTBLVafGMkZjTHQkJyfTtm3bWIdhghSTpjYRmSQi20VkebnyISKyRkTWisidVZ1DVVep6s+BEUDZePLLgddV9WbgkogEb4wxpkZi1cczBRgSWCAiicATwIVAF+AqEekiIgUi8m65n6buYy4BPgc+dE+TD2xy//ZF4XkYY4wJUUya2lT1UxFpU664L7BWVdcBiMjLwKWqej8wtJLzvA28LSL/Bl4CCnGSz2Js4IQxxsSleOrjacnR2go4SaRfZQeLyFk4TWspwEy3+A1gnIhcDLxTyeNGA6Pdm57yzX1VyAH2BXlfVbfL/i5fhnu7MRDsFoahxFRZHNXFlwMkhxBTVXGFM6ZwvVbVxRRsLOGMqboYamNMZWVEOKbK4qguPnufh/aeOjnIeCqmqjH5AdoAywNu/xh4NuD2tcDjEY5hfgjHTgj2vqpul/1dviygPCIxVRZHdfG5sQUdU1VxhTOmcL1W1cUUQiwR/fer7TEFvJfsfV7P3ucV/cRTjacQaBVwOx/YEqNYKlJhDaqS+6q6/U4VZZGMqaprVhXfO8CpYYornDGF6nhjCjaWcMZUXQy1MabjjSfUmKq6pr3Pqy6L1HvqB2K2OrXbx/OuqnZzbycBXwPnApuBecDVqroigjHM1xqssBoJFlPw4jEuiyk4FlPw4jGumsYUq+HU04A5wMkiUigiN6mqF7gV+C+wCng1kknHNSHC5z8eFlPwe4RSEQAACqlJREFU4jEuiyk4FlPw4jGuGsVUr/fjMcYYE3025NgYY/5/e2ceu0dRxvHPl0ppqVpCOSKaWCpEgZbTAxUMpyRKpKjhEMQKCVJQQjACDSYaEqIGBeSoGCiHCE2wwZaAmGq1FETC0btcBqQGA7FgJfbQYvv4x8xbttud9333fffdd9s+n+TN77fzzvGdZ2Z3dubdncepFR94HMdxnFrxgSeBpJ0kXS3pRklfG7YeCO8uSXpU0i3xPaZGIGmMpGckFb7oWzeSDog2miVp6rD1tJA0WdKtkuZI+uyw9QBImiBphqRZQ9YxRtJd0T5nDVNLi6bYJktD+1Dp8227HHiq2AsOOIXwUuvbhEe9m6DJgDXAqAZpArgcuK9fPVVpsuJ9/Jqga7aFfQSnAKc3RNPLZnZev1oq0FfLPotlNA3SNn1oqrQPVaSp/PnWz0tATf0AnyE8k599QXUE8BIwARgJLCHsCTcJeDD32Qu4AvhGTDurIZp2iun2Bu5piKYTgDMIJ8LJTdAU03wBeJzwSH4j+lQm3U+Awxumqe8+3qe+acChMc69VWvpRdMgbVOBpkr6UFWayp5vTXqBtDKsgr3gJL1KcK8AFWw4WoWmDKsJWwUNXZOkY4ExhIvHekm/MbNNw9QU88nv49cXFdlKwA+Bh81sYRM0DZIy+qhpn8WSmp4dlI5eNUl6jgr7UBWagGfLnm/b5cCToNRecIR9326UdDSwoAmaJH0ROAnYDbipCZrM7MqobQrwRj+DTlWaVLyP3yAo26e+RZghjpW0n5ndMmxNksYBVwOHSZoWB6hBktJ3Ax32Waxb0xBs01ET9fShUpp6Od92pIFHBWHJl5jMbB0w6PXdspruJwyIg6SUps0RzO6sXspmytppPjB/UGIylNV1A+ECO0jKanoTuGBwcraiUJ+ZrQW+XqOOLClNddsmS0pTHX0oRUrTfEqeb9vlwwUJmrgXnGvqjiZqgmbqaqKmLE3U55q6ozJNO9LA8xSwv6R9JY0k/CD+gGtyTX3QRF1N1JSlifpcU92aBvm0xrA+wEzgNd55FPq8GP45wkakLwFXuibXtC3raqKmputzTc3Q5Hu1OY7jOLWyIy21OY7jOA3ABx7HcRynVnzgcRzHcWrFBx7HcRynVnzgcRzHcWrFBx7HcRynVnzgcRzHcWrFBx7HcRynVnzgqQFJGyUtlrRC0hJJl0rq2vaSxucdMiXimaS7M8fvkrRK0oMd0u0m6cLUca9k6r1c0q8k7Voy/VY6JD3eIc2akmX0ZLNekTRa0iOSRsTjtvUpmXepusc0lbT1toSkkZIWSKplk+R4/q6XtDgTtrekeyW9rOC998+STu2Qz3xJJ+XCLpE0PfarxZI2SNpjUHWpCh946mG9mR1qZgcBJxK2nfheNwmjD5du22ktMFHS6Hh8IvD3LtLtBlzY5rgjCuR1tuo9keDbqOudfmO9d8/rMLNPldHVBb3arFfOBe43s40wkPqUJdnWiTbd5jGzDcA8BujBs4CXzOxQ2Ny3ZwMLzGyCmR1B2PfsAx3ymBnjZTkDmGlm62P+w95ItCu2u07VdMzsH8D5wDfjiT073vGskHQ+bL5Dek7SdGAhmR1hFfzAL5L0sUQRDwOfj/+fSeismykqj+BY6kPxjumagmMknS3pyRj2c0kj2uks4FFgvzYaiuo9o0DHmvj3HElL4wzy7q2LK9Zclc0kjZH0UCx/uaTTi8IKyjoLmJPJu1WfVt1vjeXMzQyGpeuVilNgty3auqhNFWboy+PnkjJ6Y7znJd0W098j6QRJf5L0F0kfr8rOiT6aapPZsS06ojBbv0nSY5JWSjpK0i8kvShpRjd55DgO2GAZXzpmttLMbsyUWdR+s4CTJe3Ssi2wD/BYDxqGS50bz+2oH2BNQdhqggvr3ePxaGA5MA4YD2wCjozfjY/ffRhYRHQRXFQOcDChg44ieHM8BngwEydVXtbFbf74AIJzrp3j8XTgnLzOVL0Jfp/mAFNTGjLlblXvgjoeBLwA7JHNL1dmoeYKbfYl4NZMnLFFYbmyRgKvJ2w0Hvgf77h/vg84u0Bvsl6d6l5kt0TbZ9vgCGAZwcvsu4EVwGEl9LbiTSLc6D4D3E7w7XIKMLsKO7epc2GbENw4r+ry/H0euDT+f1W04ftie64GdumQPm/ji4Hr2sRv18YPETzKAlwBXJNL+0qrfZv88RnP8Gg5VbpY0hLgCcKMYf8YvtLMnsjE35Nw8T7bzBaTwMyWEjr6mRR7A0yV147jCRegpxTWqY8n+F0v0plldIz/NPA3wgymk4Z2+bU4DphlZm8AmNk/S2regh5ttgw4QdKPJB1tZm8lwrLsAfyrTb3+mmnbZ6KmXuqVitON3WDLNjgK+LWZrTWzNQRHhEeX0NuKt8yCd9oVwDwLV8lluTT92DlV58I2sbDUuUHSexKaAZA0irAceX0MWg/MMLPXLCzZrSMsI/eMpJvjjOypGNSujbPLbWeQm51vK+xIHkgbg6QJwEbgQIIb20+a2TpJ8wl33RB+e8jyFsHt7KcJJ287HgB+TLhzH5cp95g25bWVDNxlZtNy9RhfoDNLa905m6aThnb5ZfV02la9UHMbStnMzF6UdATh97ofSJprZlcVhWXKWE97e/838/9Gwp1/L/VKtdfFdOFNli3boMjrZItu9ObjbcocbyJeg/q1M2HmUWiXNm2yC/CfNvWDMEtcaO+4dD8E+FnMt+UI7WBJ3wdOA+4ALmkN7glWEGZiAJjZRQoPBDzdkpyqC2GJ8FpJhwOjzWxhB/2NxGc8NSNpT+AW4CbgvcDqeKJ9BDiyTdINwGTgHElf6VDM7cBVZrYsFz42Ud6/geydX/54HvBlSXvFOuwu6YMdNKRIaSgiryOr5zRJ41p6EnHKaC5lM0n7AOvM7JeEAevworBsRma2GhgR76J7pZt6peIU2S1l4xYLgMmSdpU0BjiV8Htd1fRr58I6p9ok2mCVmb0dj+dJen+BrknAkszxwcDS+P8hwFIzW0JYor0DuL7DoAPwB2CUpKmZsOwTn8k2jrPO+YT+uk3OdsBnPHXRWnLambDefTdwbTy+QNJSwrpx2yUmM1sr6WTgd5LWmtmcRLxXgZ8WfPXbovLM7E2FH3uXAw+b2XcKjr8LzFV4yult4CLg9bKGSGlI1GMrXTF8haSrgUckbST87jUll/bZhOaVibJK2YxwQbpG0qaY99REWJ65hOWr36fq3Y5u6pWKY2ZP5O1mZlOyNgZuzpW3UNKdwJMx6DYzWxRnu1XSl53b2GVsQXqAY4nLqjH+fkDR0uMkYt3jDcPoeAMBWw5Cq4HnzOzprbPYEjMzSZOB6yRdBqwizDIvj993auOZhCXP/BNu2wzuCM5xakTSYYQfqr86bC07MpLuB6aZ2QuSJgLnmtmlPeZ1IPBt4E0zu6zg+/GEh1Um9iG5Wy2vAB/tYtY1VHypzXFqxMwWAX9U+tFuZ8BIGkl4mu4FADNb3segMxa4kvAu1L6SPlEQbSMwVpkXSKtG8QVSwirKpk7xh43PeBzHcZxa8RmP4ziOUys+8DiO4zi14gOP4ziOUys+8DiO4zi14gOP4ziOUys+8DiO4zi14gOP4ziOUys+8DiO4zi14gOP4ziOUyv/B2JAEUMESxA0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import fsolve\n",
    "\n",
    "k = 0.0857 #0.0857\n",
    "\n",
    "def k_crit(mx,sigma0,omega_b=0.05,omega_rad=9.23640e-5,\n",
    "                 H0=67.27,kT0=2.353e-13,rho_crit=4.755e-6,\n",
    "                 cc=3e5,Mpc2cm=3.086e24,reduction_factor=1.,mp=0.511e-3,YHe=0.25,n=0.):\n",
    "    N0 = (pow(2.,(7./2.)))/(3.*np.sqrt(np.pi))\n",
    "    B = cc*Mpc2cm*N0*omega_b*rho_crit*((1.-YHe)/1836.152)*np.sqrt(kT0)\n",
    "    zth = np.power((H0*(np.sqrt(omega_rad))*(np.sqrt(mx+mp))*(np.sqrt(mx*mp))*((mx+mp)/mx))/(B*sigma0), 2./3.)\n",
    "    #print(zth)\n",
    "    H = H0*np.sqrt(omega_rad)\n",
    "    def f(z):\n",
    "         return(B*(z**3.)*sigma0*(1./(mx+mp))*(((z/mp)+((z**2.)/(mx*zth)))**(0.5))) - (H*(z**2.))\n",
    "    zcrit = fsolve(f, zth)\n",
    "    #print(zcrit)\n",
    "    kcrit = (2*H0*(np.sqrt(omega_rad))*zcrit)/cc\n",
    "    return zth, zcrit, kcrit\n",
    "\n",
    "masses = np.array([1e-6, 1e-3, 1., 1e3])\n",
    "sigmas = np.zeros(len(masses))\n",
    "sigma_array = np.logspace(-18,-30,1000)\n",
    "\n",
    "for i in range(len(masses)):\n",
    "    mx=masses[i]\n",
    "    print(\"mx:\", format(mx, '.0e'))\n",
    "    \n",
    "    zdec, zdecprime, kdec = k_crit(mx=mx, sigma0=sigma_array)\n",
    "    sigmas[i] = sigma_array[np.argmin(np.abs(kdec-k))]\n",
    "    print(\"zth:\", zdec[np.argmin(np.abs(kdec-k))])\n",
    "    print(\"zcrit:\", zdecprime[np.argmin(np.abs(kdec-k))], \"\\n\")\n",
    "    \n",
    "print(sigmas)\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "electron_chains = np.array([9.90e-28, 5.78e-27, 5.53e-24, 8.71e-21])\n",
    "\n",
    "f,ax = plt.subplots()\n",
    "\n",
    "ax.plot(masses, sigmas, marker='.', label='Analytical')\n",
    "ax.plot(masses, electron_chains, marker='.', label='Numerical')\n",
    "ax.set_title('DM-electron scattering (n=0, Vrel=0)')\n",
    "ax.set_xlabel('Dark Matter Particle Mass (in electron masses), $m_\\chi$ $[\\mathrm{GeV}]$')\n",
    "ax.set_ylabel('Interaction Cross Section, $\\sigma_e$ $[\\mathrm{cm}^2]$')\n",
    "ax.set_xscale('log')\n",
    "ax.set_yscale('log')\n",
    "ax.set_xlim([1.0e-6, 1.0e3])\n",
    "ax.set_ylim([1.0e-30, 1.0e-20])\n",
    "ax.fill_between(masses, sigmas, 1.0e-20, alpha=0.2)\n",
    "ax.fill_between(masses, electron_chains, 1.0e-20, alpha=0.2)\n",
    "plt.legend(loc='lower right')\n",
    "plt.savefig('dme_analytical.pdf')\n",
    "plt.savefig('dme_analytical.png', dpi=300)\n",
    "\n",
    "print(\"Percent errors:\")\n",
    "for i in range(len(sigmas)):\n",
    "    error = abs((electron_chains[i]-sigmas[i])/sigmas[i])\n",
    "    print(format(masses[i], '.0e'), \":\", format(error, '.1%'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mx: 1e-06\n",
      "zth: 2267950.5272753322\n",
      "zcrit: 99187.78207406001 \n",
      "\n",
      "mx: 1e-03\n",
      "zth: 123469.78142315909\n",
      "zcrit: 96223.05494561423 \n",
      "\n",
      "mx: 1e+00\n",
      "zth: 98612.69568184736\n",
      "zcrit: 98579.12163640605 \n",
      "\n",
      "mx: 1e+03\n",
      "zth: 98109.07542996648\n",
      "zcrit: 98109.04200749476 \n",
      "\n",
      "n = 0.0\n",
      "[5.34229330e-28 6.74262224e-27 5.08987019e-24 5.12518693e-21] \n",
      "\n",
      "mx: 1e-06\n",
      "zth: 459876.015831241\n",
      "zcrit: 96411.42398807018 \n",
      "\n",
      "mx: 1e-03\n",
      "zth: 111001.37269913143\n",
      "zcrit: 96663.29929870977 \n",
      "\n",
      "mx: 1e+00\n",
      "zth: 97341.62881250696\n",
      "zcrit: 97321.745389593 \n",
      "\n",
      "mx: 1e+03\n",
      "zth: 97062.85883189054\n",
      "zcrit: 97062.8389922553 \n",
      "\n",
      "n = 2.0\n",
      "[8.99402217e-27 1.70699493e-23 1.92879151e-20 1.94217468e-17] \n",
      "\n",
      "mx: 1e-06\n",
      "zth: 278193.6235244337\n",
      "zcrit: 98209.98040013728 \n",
      "\n",
      "mx: 1e-03\n",
      "zth: 107824.5749749279\n",
      "zcrit: 97996.78993551355 \n",
      "\n",
      "mx: 1e+00\n",
      "zth: 96813.78429869135\n",
      "zcrit: 96799.6592799096 \n",
      "\n",
      "mx: 1e+03\n",
      "zth: 98233.78236843963\n",
      "zcrit: 98233.76802631743 \n",
      "\n",
      "n = 4.0\n",
      "[6.02254120e-26 3.05833803e-20 5.47947234e-17 5.20854855e-14] \n",
      "\n",
      "mx: 1e-06\n",
      "zth: 213473.0131681883\n",
      "zcrit: 97779.01264134399 \n",
      "\n",
      "mx: 1e-03\n",
      "zth: 104920.3330685932\n",
      "zcrit: 97535.14653993762 \n",
      "\n",
      "mx: 1e+00\n",
      "zth: 97925.39349195482\n",
      "zcrit: 97914.28140438127 \n",
      "\n",
      "mx: 1e+03\n",
      "zth: 97791.67238872261\n",
      "zcrit: 97791.66128394287 \n",
      "\n",
      "n = 6.0\n",
      "[3.02329468e-25 4.60960449e-17 1.16698982e-13 1.17508713e-10] \n",
      "\n",
      "[[5.34229330e-28 6.74262224e-27 5.08987019e-24 5.12518693e-21]\n",
      " [8.99402217e-27 1.70699493e-23 1.92879151e-20 1.94217468e-17]\n",
      " [6.02254120e-26 3.05833803e-20 5.47947234e-17 5.20854855e-14]\n",
      " [3.02329468e-25 4.60960449e-17 1.16698982e-13 1.17508713e-10]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import fsolve\n",
    "from scipy.special import gamma\n",
    "\n",
    "k = 0.42\n",
    "\n",
    "def k_crit(mx,sigma0,n,omega_b=0.05,omega_rad=9.23640e-5,\n",
    "                 H0=67.27,kT0=2.353e-13,rho_crit=4.755e-6,\n",
    "                 cc=3e5,Mpc2cm=3.086e24,reduction_factor=1.,mp=0.511e-3,YHe=0.25):\n",
    "    N0 = pow(2,(5.+n)/2.)*gamma(3.+n/2.)/(3.*np.sqrt(np.pi))\n",
    "    B = cc*Mpc2cm*N0*omega_b*rho_crit*((1.-YHe)/1836.152)*pow(kT0, (n+1.)/2.)\n",
    "    zth = np.power((H0*(np.sqrt(omega_rad))*(pow((mx*mp)/(mx+mp), (n+1.)/2))*((pow(mx+mp, 2.))/mx))/(B*sigma0), 2./(n+3.))\n",
    "    #print(zth)\n",
    "    H = H0*np.sqrt(omega_rad)\n",
    "    def f(z):\n",
    "         return(B*(z**3.)*sigma0*(1./(mx+mp))*(pow((z/mp)+((z**2.)/(mx*zth)), (n+1.)/2)) - (H*(z**2.)))\n",
    "    zcrit = fsolve(f, zth)\n",
    "    #print(zcrit)\n",
    "    kcrit = (2*H0*(np.sqrt(omega_rad))*zcrit)/cc\n",
    "    return zth, zcrit, kcrit\n",
    "\n",
    "masses = np.array([1e-6, 1e-3, 1., 1e3])\n",
    "sigma_array = np.logspace(-5,-30,1000)\n",
    "ns = np.array([0., 2., 4., 6.])\n",
    "sigmaout = np.array([np.zeros(len(masses)), np.zeros(len(masses)), np.zeros(len(masses)), np.zeros(len(masses))])\n",
    "\n",
    "for j in range(len(ns)):\n",
    "    sigmas = np.zeros(len(masses))\n",
    "    for i in range(len(masses)):\n",
    "        mx=masses[i]\n",
    "        print(\"mx:\", format(mx, '.0e'))\n",
    "    \n",
    "        zdec, zdecprime, kdec = k_crit(mx=mx, sigma0=sigma_array, n=ns[j])\n",
    "        sigmas[i] = sigma_array[np.argmin(np.abs(kdec-k))]\n",
    "        print(\"zth:\", zdec[np.argmin(np.abs(kdec-k))])\n",
    "        print(\"zcrit:\", zdecprime[np.argmin(np.abs(kdec-k))], \"\\n\")\n",
    "    \n",
    "    print(\"n =\", ns[j])\n",
    "    print(sigmas, \"\\n\")\n",
    "    sigmaout[j] = sigmas\n",
    "\n",
    "print(sigmaout)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7ffb8c0b5970>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAEPCAYAAAByRqLpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd3Sc5Zn4/e89XTOj3qxmW82yZVuuuMpgTDPYGGIgGAgQIDiwgSS/vNnfJpu8C5vdDWx29z1v8iNvNmwKJwkLu4cAIYmzMUlMsSmxwb3bsnrXSKPp7bnfP57xuMQ2llVmJN2fczjHetA8zyXL0jV3uy4hpURRFEVRxooh2QEoiqIok4tKPIqiKMqYUolHURRFGVMq8SiKoihjSiUeRVEUZUypxKMoiqKMKVOyA0imrAynnFKQm+wwFEVRxg2D0cThYw29Usr8K73HpE48hXnZ/H/fejLZYSiKoowbtsxcVq5/uGk491BTbYqiKMqYUolHURRFGVMq8SiKoihjSiUeRVEUZUxNmMQjhKgVQvy3EOIHQog7kx2PoiiKcmEpkXiEED8RQnQLIQ6cd32tEOKoEOKEEOJrn3Cbm4H/I6V8HHhg1IJVFEVRhiVVtlO/ADwH/Oz0BSGEEfg+cAPQCuwUQrwBGIFnznv9w8DPgaeEEBsAdThHURQlRaVE4pFSviOEmH7e5SXACSllA4AQ4mXgNinlM8D6i9zqC/GE9epoxaooiqIMT0oknosoAVrO+rgVWHqxT44nrr8FHMC/XOLzNgObAQryskYgTEVRFGUoUjnxiAtcu2i7VCllI/GEcilSyueB5wFqKspU+1VFUZQxlhKbCy6iFSg76+NSoD1JsSiKoigjJJUTz06gWghRLoSwAJuAN5Ick6IoijJMKZF4hBAvAe8DNUKIViHEI1LKKPAE8HvgMPDfUsqDyYxTURRFGb6UWOORUt5zketbgC1jHI6iKIoyilJixKMoiqJMHirxKIqiKJfN2NhGock0ZTj3SImpNkVRFCXFaRppe4+T/tv3yDWaSoZzK5V4FEVRlL8kJcZ+D9aGNiyn2rE2dmAIhIALH7IcCpV4FEVRFAAMvoCeZE61YznVjmnAC0A0OwP/oplEcjPJ2vIeEjmsw/cq8SiKokxSIhzB0tyFpaEN66l2zF0uALQ0K8FZ03HXlhOsLSdamANCH+cYFs2j774vDeswv0o8iqIok4WmYW7vTYxoLC1diJiGNBkJVZXiXVlHsLac8LQiMF5471l0xnS6otHO4YShEo+iKMpEJSXGPjfWU+1YG9qxNHVgCIYBCJcVMnjDEoK15YSqpyKt5jELSyUeRVGUCcTgDWA51Ya1QV+rMQ76AIjmZuJbUkuwtpzgzOloGY6kxagSj6IoyjgmwhEsTZ1n1mm6+wGI2W16kqmdrq/T5Gcn1mmSTSUeRVGU8SSmYW7v0afOTrVhae1BaPo6TXDGVLxXL9DXaaYWgiE1awSoxKMoipLKpMTU606MaCxNHRhCEaSA8NQpDN60jGDtdELVZUjL2K3TDIdKPIqiKCnG4PHpSeb0Oo3HD0AkPxvf8rnxdZppaE57kiO9MirxKIqiJJkIhbE0dmI9pVcJMPcMABBzpsXXacoJzpqur9NMACrxKIqijLWYhrmtO77zrA1zWw9Ck2hmE6EZU/GuXkSgtpxIWSEYUmNDwEhSiUdRFGW0SYmppz8xdWZp6sQQjiCFIDy9iMGbVxCoLSdUVQrmif9reeJ/hYqiKElgGPRhPb0h4FQ7Rm8AgEhhDr6VdXqimTkNzZGW5EjHnko8iqIoI0AEw1iaOvRqzg3tmPvcAMTS7QRnVxCIr9XEcjOTHGnyqcSjKIpyJaIxLK3diWrO5rYehJRoFjOhmql4rr+KYG05kZKCCblOMxwq8SiKolwOKTF1uRLbnC3NnRgiUaRBECovxr1+pV73rLIUTMZkR5vSxm3iEUJUAN8AMqWUd8av3Q6sAwqA70sptyYxREVRxjmD26vvPGtow9LYjtEXBCAyJTdRISBYMxVptyU50vElKYlHCPETYD3QLaWcc9b1tcB3ASPwIynlsxe7h5SyAXhECPHKWddeB14XQmQD/wqoxKMoymUTgRDWxo5EkU2TaxCAaIaDQF1V/DxNObGcjOQGOs4la8TzAvAc8LPTF4QQRuD7wA1AK7BTCPEGehJ65rzXPyyl7L7E/b8Zv5eiKMrFRaNYWrrjdc/aMXf06us0VgvBmdPi5WjKiRTnpUyBzYkgKYlHSvmOEGL6eZeXACfiIxmEEC8Dt0kpn0EfHX0iIYQAngV+J6X8eOQiVhRlQpASU2dfvD9NG5bmLkQ0hjQYCFWW4N6wSl+nKS9W6zSjKJXWeEqAlrM+bgWWXuyThRC5wD8BC4QQX48nqCeB64FMIUSVlPLfL/C6zcBmgIK8rBEMX1GUVGTs9ySmziyNHRj9+jpNuDgPz+pFZ9Zp0qxJjnTySKXEc6FxrLzYJ0sp+4DHzrv2PeB7l3qIlPJ54HmAmoqyi95fUZTxSfiDWBs7EtWcTf0eAKJZTgILZiTqnsWy0pMb6CSWSomnFSg76+NSoD1JsSiKMl5EolhauvSps1PtmDv7EBI0m4XgrOkM3qKXo4lOyVXrNCkilRLPTqBaCFEOtAGbgHuTG5KiKClH0zB39p2pe9YSX6cxGghVluJeOodAbTnh8mIwpmYjtMkuWdupXwJWA3lCiFbgKSnlj4UQTwC/R9/J9hMp5cFkxKcoSgqREmO/JzGisTZ2YAiEAAiXFuBZs1ivezZjKtJmSW6symVJ1q62ey5yfQuwZYzDURQlxRh8ASyNHYlkYxrwAhDNzsC/aKZe92zWdLRMZ5IjVa5EKk21KYoyWUWiWJo69W3O8XUaAC3NSnDWdNyzKwjWlhMtyFbrNBOASjyKoowpc0sXlsYONIcNoy+I5fQ6TUxDmoyEqkrxrtDbO4enF4FBrdNMNCrxKIoyJoQ/iOPPB3G+uxekTJyfCE0tZPCGpQRrpxOqnoq0mpMapzL6VOJRFGV0SImp1431WDO2482YW7oRUiLRD+1JAe5bVuK+49pkR6qMMZV4FEUZOTENS3OnnmyOtWDq14tshqYW4r61nmhuJjm/+B+IxZBGI4F51UkOWBmqBl8TpmzTlOHcQyUeRVGGRQRCWE+0YjvWjPVkK4ZgGGky6psC1q8kMK+KWM6ZrpuRojxsR5sI1kwjXFWaxMiVoTria+CFjl9iyjSVDOc+KvEoijJkxj63nmiOt2Bp7kRokli6Hf/iWfjn62VpLnamJlxVqhLOOBHWIpwKtHI80MgJfxNtoc4Rua9KPIqifDJN08vSHGvBeqwZc58bgHBJPoM3r8A/v5pweYlq8TzOxaRGa7CDE4EmjvsbaQq2EZUxjMLITHsFsxyVvN3/IciL19G8HCrxKIpyQSIYxnqyVV+vOdGKIRBCGg0EZ07Hs3YZ/nnVxFSF93FNSkl3uI/jgUaO+xtpCLQQ1PSqEOW2UtblXUudcyazHJXYjHr17puKruPWwXuHVUdTJR5FURKM/YOJjQGWpk6EphFzpOFfMIPA/BkEZleo9gHj3EBkMDGiOeFvYjCmV4UotOSxKusq5qXPZI5zBhmmC1eFmJVeTbQ/Oqw5N5V4FGUy0zTMbT36es2xFsw9/QCEi3IZvGkZgfnVhCpL1CHOccwfC3Iy0MRxfxMn/I30RFwAZBid1KXPZK6zhjpnDYXWvDGLSSUeRZlkRCiMtaFNX6853oLRH0QaDARnTMVz3WIC86uJFuQkO0zlCkW0KI3B1sSIpjXUiURiNViY7ahmXf611DlrmGorxiCS84ZCJR5FmQQMbm98VNOMtbEDEdOI2W0E6yrxz59BYE4l0m5LdpjKFdCkRluoi+N+fZ2mMdiqbwjAwAx7OXcV3sw850yq7NMxG1LjV35qRKEoysiSEnNbT7xqQAvmLn16JVKYw+D1S/QptKoy1a9mHJJS0hNxccKvr9OcDDQT0PR23tNsJdycew116TOpdVSRZkzNNxMq8SjKBCHCESyn2hPna4zeANIgCFWV4blmoT6FNiU32WEqV2Aw6k2MaE4GmhiI6u288805rMhayDznTOY4a8gyj4923irxKMo4Zhj0nTuFFo2h2SwE6qoIzKsmMLcKzZmW7DCVIQrEgjQEWvR1mkATXeFeANKNDuY6a5ibXsM850wKLXmIcdgmQiUeRRlPpMTU0Rcf1TRj6dD71kTys/CsXkRgfjXBGVPBZExyoMpQRLUoTcE2fedZoJGWYAcaEoswU+us5qbcVdQ5ZzI9rSRpGwJGkko8ipLqIlG9QVp8vcbo8SOFIFRZQv+dawjMqyZSnKcapI0jmpS0xzcEnAg0cSrQQkRGMWCgyj6NjQVrmZdewwx7OWbDxGsToRKPoqQgg9eP9VgLtmPNWE61Y4hE0awWAnMqCMyPT6FlOJIdpnKZpJT0RQbiiaaRE/5m/FoAgDJrETfm1lPnnMlsZzV248SfGh23iUcIUQF8A8iUUt4Zv7Ya+AfgIPCylPKtpAWoKEMhJaYuV+Igp6W9B4Bobibeq+cTmFdNsGYamMftj+yk44n6zqoQ0Eh/VG8RkWvOYmnmvMThzRxz5ifcaeJJyr9iIcRPgPVAt5RyzlnX1wLfBYzAj6SUz17sHlLKBuARIcQrZ18GvIANaB2N2BVlxERjWBs79I0Bx5sxuX1IAeHyYvo/tZrAghlESvLVFNo4EdRCnDq9IcDfREdYf/PgMNqZ65xBnXMmdek1FFkKxuWGgJGUrLdPLwDPAT87fUEIYQS+D9yAnjR2CiHeQE9Cz5z3+oellN0XuO+7Usq3hRCFwP8D3DcKsSvKFTP4AliPx6fQGtowhKNoFjPB2RUMzK8mUFeFlnnhGllKaonKGM3B9kSiaQ62o6FhFmZmOSq5Lncldc4aytPKME6ADQEjKSmJR0r5jhBi+nmXlwAn4iMZhBAvA7dJKZ9BHx1dzn21+B/7AVXJUEk+KTH19CfWa8xt3QgJ0ax0fCvq9N41s6arKbRxQJOSznB3ItE0BFoIywgCQWXaVG4vuIF5zpnUOCqwTMANASMplf61lwAtZ33cCiy92CcLIXKBfwIWCCG+LqV8RgixEbgJyEIfUV3odZuBzQAFqqS7MhpiMSxNnYmDnKZ+/bBfaHoR7tuuxj9vBpGphWoKbRxwxTcE6Nucm/DF/ACUWAu5LmcFdek1zHHOwGG0JznS8SWVEs+Ffgov2mxIStkHPHbetVeBVy/1ECnl88DzADUVZcNqZqQopwl/ENuJFr3w5sk2DKEwmtlEsLYc962rCMyrJpY9Pk6VT2beqD9eybmR44EmXJEBALJNmSzOmBNfp5lJrlm9aR2OVEo8rUDZWR+XAsNqNqQoo0bKM+2fjzVjaelGSEk0w4Fv6Wz9IOescqRVTbmksrAWpiHQmjhP0x7qAsBusDHHOYO6ghupc9ZQYp0y6TcEjKRUSjw7gWohRDnQBmwC7k1uSIpylpiGpbkzcZDT5NK3x4bLCnGvX0lg/gzC04pU++cUFpMxWoIdiUTTFGgjhoZJGJlpr2R19lLqnDVU2qdiFKr6w2hJ1nbql4DVQJ4QohV4Skr5YyHEE8Dv0Xey/URKeTAZ8SnKaSIQwnqyVR/ZnGjFEAwjTUaCM6fjXreSwLwqYjmT7xzGeCGlpCvcG586a6TB30JIhhEIytNK2ZB/PXPTa5jlqMRqsCQ73EkjWbva7rnI9S3AljEOR1HOkZhCO96CpbkToUli6Xb8i2fhn1dNcHYF0qZ+SaWq/og7MaI54W/CE/MBUGTJZ3XOUuqcemvndJOq/JAsqTTVpijJoWmYW7oTycbcqy8oh0vyGbx5Bf751YTLS9QUWoryxQKc9DdxPKC3du6N6O27M03pzE+vpS59JnXOGvItqqtqqlCJR5mURDCM9WSrvl5zohVDIIQ0GgjWTMNz41L886uJqe32KSmsRTgVaOVEQN/m3B7qRAJpBiuzndVsyL+Ouc6ZTLUVqQ0BKUolHmXSMPZ74hsDmrE0dertnx1peuvn+dV6++c0de441cSkRmuwI1H3rCnYprd2FkZq7OWsyl5PnbOGKvt0TGpDwLigEo8ycWka5raeROFNc48+BRMuymXwxqUE5lUTqioFgypnkkqklHRH+hIVAk4GmglqIQCm20pZl3ctc5011DqqsBnVG4XxSCUeZUIR4QjWk2164c0TLRh9Qb3984ypeK5bTGBeNdFCNdefagYig2dVcm5iMOYFoNCSx6qsqxIVAjJN6hDuRKASjzLuGdzec9s/xzRidivBuir882YQmFuJtNuSHaYS1xho44jvJBaDCXfUy3F/Iz0RFwAZRidz02v0CgHOGgqteUmOVhkNKvEo45JtzzHs+05gdHsTtdAiBTl4rrsK//xqQlVlqv1zCjndcfP9gd3s9OxDxqthmYWJuc4a1uXr02fTbMUTorWzcmkq8SjjSyRK1i+3kXasWf/VJQSDaxbjue4qokW5yY5OOYs/FuS4/xSHfSc55j+VOE9zmkBwR8FaPj3lliRFqCSLSjzKuGHqcpH1y22YeweQ6FVlpYBYdrpKOilASklHuJsjvgaO+E7SFGxDQ+Iw2lmYXsvCjNk4jQ7+pfE/iMoYJmFkXvrMZIetJIFKPErqkxL7zkNk/GEnMbsN1z03kPXKNojFkEaj3hJaSYpgLMTxQCNHfA0c9TfgjurTnhVpZWwsuImFGbOptk8/p+7Z31d+iQPeY8xxzqDGUZGs0JUrZOk9QEm6mDKce6jEo6Q04Q+S9ca72I41459bSd8jG9AyHITKS7AdbSJYM41wVWmyw5w0pJR0h/s47D/JEd9JTgVa0dBIM9iYnz6LhRlzWJBeS4754vXrahwVKuGMR1LibNnGlN3fY4pTlAznVirxKCnLcqqdrNffxhAI4brnRjzXX5VonhauKlUJZ4yEtTAn/E0c9jdw1HeS/qhelXuqrZjb8q9nYcZsahwV6vDmBGTy95DWsxd7zx7sPXsxhfpH5r4jchdFGUkxjfS3Psbx3l6ihbl0/V/3Epk6rJG9MgRSSnoj/RzxneSw/ySnAi1EZQybwUqdcyaLMmazIH02eZbsZIeqjDBD2Eta7/5EorF4WwGI2nLxll5DOL2UvP3/gcQ7rCaaKvEoKcXYP0jWq29haevBc/V8+u+5EWlVlaBHW0SLcDLQzBHfSY74T9EXL7RZYp3CLXmrWZQ+h5mOSswG9StjIhGxCDbXYew9e0nr2YOt/zgCjZjJjq9oGb2zH8ZbUk8ou+bMbEPFDXT9y5phNem8rH9FQojLOeqtSSkHhhOMMrnZ9p8kc8sOMBjoeXwj/qtqkx3ShNYXGdATja+Bk4EmIjKKRZipS6/hjvjGgAKL2i04oUgNy2Aj9u492Hv2kNZ3EEMshBRG/Pnz6V7wJN6SegL5C5DGC7/hC05ZTJtHdg4njMt9+9Ie/+9SpV6NwNThBKNMTiIUJuN/PsC+9zjBylJ6P3+7qgw9CqJalIZAC0f8+nbn09UCiiz53Ji7ioXps5ntrMZiUO26JxKTv/tMoundhymkjw+CWVW4Zt6Lt7geX9EyNMvYlSO63MRzWEq54FKfIITYPQLxKJOMqb2X7Fe3Yez3MLBhFe5bV4FRnVwfKf0RdyLRnPA3EZYRzMLEbOcMNuRfx4KM2RRbC5IdpjKCDGEP9p59pJ1ep/Hps2KRtAI8ZWvwltTjLa4n6kjeuunlJp7lI/Q5iqKTEscHB0j/0y5iGQ66/vdnCKnzOMMWkzEaA62JHWid4V4A8s05rMlZzsKM2cx11qg2zxOIiIWxuQ7FRzV7sQ6cQCCJmR34ipbTO3ezvk6TVZ1Yp0m2y0o8UsrgSHyOogAYvH6yfvUO1pNt+BfW0PfZ9WjOtGSHNW65ox6O+ho47DvJ8UAjIS2MURipdVSxNu8aFqbPpsRaqJqiTRRSw+pu0KfOuvfq6zRaGClM+AsW0L3wf+EtqcefPw9SdNr0ExOPEOIG4NPA96WUe4QQm6WUz49+aMpEZD3eQuYb7yLCEfoeuBnvNQtT5l3YeBGTGs3BNr00jf8k7aFuAHLNWVydtYRF8VFNmlFV5J4oTL5OfYtz9x7svfswhvWzVMHsGly19+vrNFOWolmcSY708lzOiOevgIeAb8Z3t80f3ZAujxCiAvgGkCmlvDN+bSrwHNALHJNSPpvEEJWzRWOk/2kXzg8OEC7Jp/exjURK8pMd1bjhifo4Gl+rOeZvJKAFMWBglqOS63JWsDB9NlNtxWpUM0EYQm7svfsS02dmv76JLGIvxD3tJnwl9XiLVxC1FyY50itzOYmnJ75N+qtCiGeBq4b7UCHET4D1QLeUcs5Z19cC30XfIfejSyUOKWUD8IgQ4pWzLs8Afiul/KEQ4mfDjVMZGcY+N1m/3Ials4/BNYsZ+PR1SEtqTgGkCk1qtAQ7EhsDWkP6L54sUwYrshawIH0289Jn4jDakxypMhJELERa3yF9Q0D3Xqzuk/F1Gife4pV4S76At7iecGbFhJghuJzE89vTf5BSfk0I8eQIPPcF9JFJIjkIIYzA94EbgFZgpxDiDfQk9Mx5r39YStl9gfvuBr4hhLgb+PkIxKkMh5Sk7T1Oxv+8jzSb6H7yLgILapIdVcryxfwc853isL+BY74GfFoAgWCGvZx7p9zKwvQ5TE8rUf1qJgIZwzpwMj59theb6xAGLYJmMOMvWMRA1e36eZq8OpiAh3Y/8SuSUv7qvI//z3AfKqV8Rwgx/bzLS4AT8ZEMQoiXgduklM+gj44ux0PAU/H7vwL8dLixKldGBMNk/nYHaQcbCNZMo3fzbcSyM5IdVko53RztsO8kR/0naQ62I9G7cC7OnMui9DnMS59FusmR7FCV4ZISs68de89e7N36eRpjRG/vHciZRd/sh+LrNEuQ5ok/ih1SKhVCLEZfV5kWf60ApJSybgRiKQFazvq4FVh6iVhygX8CFgghvh5PUP8DPC2EuBdovMjrNgObAQrUIcVRYW7pIuu1tzAO+ujfuJrBW1aAQb1Lh3Obox31n8Ib8yEQVKZN5dOF61iYMZvKtKlqVDMBGEMDeoHN0+s0AX2SJuwoxl2+Pn6eZgWxtMnX3nuoY7gXgb8G9gPaCMdyoYnLixaik1L2AY+dd+0AcOelHhLfkfc8QE1F2bAK3Snn0TScO/bhfPtjYtkZdH79QcKVk7uC9OU0R5ufXkumaexOjSujQ0SDpPUdTBTYtLobAIhZMvR1muJ6vCWrCGdMmxDrNMMx1MTTI6V8Y1Qi0Uc4ZWd9XIpepkcZBwyDPrJefxtrYwe+JbX0PXAL0j45t/OeaY52kqO+BtwxfUpFb462loUZtVTbyzGqUc34psWwDRxPbAhIcx1GyKi+TlN4Ff3Vd+ItWUUgdw4YVMuIsw018TwlhPgR8EcgdPqilPLVEYhlJ1AthCgH2oBNwL0jcF9llFmPNpH1xrugafQ+fCu+lXWT6h2dlJKucG9iB9rp5mh2g4356bUsyJj9ic3RlHFASsze1sSGgLTefRijfiSCYO5seud+Lr5OcxXSpA5EX8pQE89DwEzAzJmpNgkMKfEIIV4CVgN5QohW9A0BPxZCPAH8Hn0n20+klAeHGJ8yliJRMt78M45dhwlNLaT3sY1Ep0yOasaheHO0I+c1R5tmK1HN0SYQY7A/MXWW1rMHc0AvQRR2luKuvC2+TrOSmO1yCvgrpw018cyTUs4d7kOllPdc5PoWYMtw76+MPlN3P1mvbsPc3c/gjUvpv+NaME+8bZ+nSSnpibj00jSqOdqEJSL+c9dpBhsBiFqz4udp9AKbkQxVV3A4hvqb4gMhRK2U8tCoRKOkPimxf3SEjK0foqVZ6fpf9xCcW5nsqEbFuc3RGuiL6OXkVXO0CUSLYus/lkg0NtcRhIyhGa34CpfgqrkHb0k9wdzZoNbkRsxQf2LqgQeFEKfQ13hGcju1kuJEIETWG+9iO9pEYE4FvY9sQMscH7WhLteZ5mgnORloPq852lrVHG28kxKLpzleYHMP9r4DGKIBJIJA3lx66h7TC2wWLEKaJufmmLEw1MSzdlSiUFKepbGDrNffxuAL4Pr09XhuXAqG8b+B4ExzNL0Tp2qONvEYA736wc34qMYU1L/HoYxp9FfdgbekHl/RCmI2da5vrAw18XwL+NLpFtdCiGzg34CHRzowJUVoGs63d+PcvpdofjbdX95EeHpRsqMalv6IO1HZ+ezmaHNUc7QJwRDxkdZ7QB/V9OzB6tHPpUet2XhLViU2BETSyz7hTspoGWriqTuddACklP1CiEt2JlXGL+OAh6zX3sLS0o13ZR2u+9YibeOvgVg03hztiGqONjFpEWyuo/ERzR5s/ccQUkMz2vAVLcM160F9nSZnplqnSRFDTTwGIUS2lLIfIN4mQa2sTkC2gw1k/nYHAD2bb8e/bM4nvCK1qOZoE5iUWAab4olmN2m9BzHEgkhhIJBXR888vZKzv3Ah0mhNdrTKBQw1afwb8F68AKdEbxD3TyMelZI0Ihwh4/cfYN99jFBFMb2bP0W0IPW3CKvmaBObyd+TmDqz9+zDFOoHIJhZQX/N3Xo5mqJlaFZ1SHc8GFLikVL+TAixC1iDvqNto9paPXGYOvvI+uU2TC437nUrGLjtGjCl5gHIxkAbh30nEAh6Iy7VHG2CMYS9pPXuS2wKsHjbAIik5eEpXa1vCCheScRZnORIlSsx5GmyeKJRyWYikRL7nw+S8YedxNLtdH/1PoKzypMd1QVpUvIn13tsdb2bqCDrNNpZkbWAhelzqEuficOoypWMNyIWweY6nFinsfafQKARM9nxFS2jd/YjeEvqCWXXTKpyTBPVZSUeIcTHUsqFw/0cJfUYfAEy33gH2/FW/POq6Xv4VrT01OsHoknJAe9R3nTtoDPck7guEGzIv447C29OYnTKUNhch0nr2UfElocp7Nan0PoOYoiFkMKIP38+3Que1Buh5S9AGtWmj4nmckc8s4QQ+y7x/wWgJlfHGcvJNrJ+9TaGYBjXfTfhWbM45d5NSik56DvOm67ttIe6KbYWcnfhOl7r3j5uHF4AACAASURBVEpUxjAJI3OdqqvpeGDy95DZuIXs478EqSX6oASzqnHNvFcvsFm0DM2iWkRMdJebeGZexufEhhOIMoZiMdK3fYTzvf2Ei3Lp+uvPECkrTHZU55BScsh3gjdd22kLdTHFks+Xpj5IfdZVGIWB+emzOOA9xhznDGocFckOV7kAQ8RPWu9+7N274+s0rYC+K0kAEkFP3eN0LflaUuNUxt5lJR4pZdNoB6KMDaNrkKxX38LS3oPnmgX0b7oRaU2dU/lSSo74T7K1bzutoU4KLXk8WfYAV2dfhfGsSs81jgqVcFKNFsM2cCzecXPPWXXPbPF1moeJ2nIp3f43oEWQBjOeaTckO2olCdQZnEkkbd9xMra8jzQZ6P7CnQQWXc5AdmxIKTnqP8XWvndpCXWQb87hC2X3c032EtVaIFVJidnXEe9Ps1vvTxPxnVf3bBX+wkXnnKeJZEzF0fE+vqLl+AsXJfELUJJFJZ5JQITCZGx5D/v+kwSry+jdfDux3NRYkpNSctzfyFbXuzQF28kzZ/N46X1cm7NMJZwUZAh7ztQ9696N2d8FQNhRjLviVr0kTXE9MdvFz375CxephDPJqcQzwZnbesh69S2MAx4Gbrsa9/p6MCa/bIiUkhOBJrb2bacx2EquOZvPl97Dmuzlqs1AKtEipLmO6Os03XuwDhxHIImZnXiLV8Rrn60inFGechtTlNHx51MujM7cKcO5xxX/hAsh7gU2oG8qEMCvpZQvDScYZQRJieO9/aRv+4hYlpOuv7mf0IypyY4KgJP+Jra6ttMQaCHHnMWjJXdzfc4KzKoCdPJJicXTEt8QsJu03gPxcjSntzl/SZ8+K5gP6vs1KQQjMQ51DLK3ZYAPGvpodwcxOrJLhnPP4by1vEZKuen0B0KI7wMq8aQAg8dP1q/extrQjm/xLFwP3oLmSP6hyoZAC1v73uVkoJlsUwafK/k01+esVC0HkizR3jm+KcAU7AMglDGd/hmfxlN6dXybc0aSI1XGQkyTnOzxsqdlgL2tAxzuGCQSk5iNgikZNn0b/DAHt8NJPFYhxDqgBSgFkv+bTcF6rJnMN95FRKL0fXYd3lXzkz4FcirQyta+dzkRaCLTlM7DxXdyQ269qgadJCIW0ts7xxON1d0AnG7vXJ+YPouklyY5UmUsSClpHwiyp3WAvS0D7G9z4w1FAagtyuChleWsrMrjqunZHO7wcN+PPtBfNAxiKK8XQjwlpfz7+J/twEb0pNMCvCal9A8nmLFWU1Emf/CPX0p2GCMjGiPjjztxfHiQcGkBPY9tJFqcl9SQmgJtbHVt55j/FJkmJ58quImbclephDPWpIbVfQp7j75OY+s7iEGLoBnM+AsXJxJNIHcOGNSGjsmg3x9mX6ubPS397G1x0+MNAVCSlcaq6jxWVuWxojKXXOdfVvf+qKmfpXOq26Ke3it+ZzLUEc9T8YSTA3wMvHy6RcJYE0LcDqwDCoDvSym3CiEqgG8AmVLKO5MRVzIYewfI/uU2zF0uBq+/iv67rgNz8hboW4IdbO17lyP+BtKNDh4o+hRrc6/GpkrUjxlToDdxcDOtZy+mkN5GK5hdg6v2QTwlq/BNWYo0p155JGXkBcIxDna42RufPjvVq48RMtPMrKjMpb46j/qqPKbm2D+xqO6iadnEvH2dw4lnqL+dJBAEfg8sRG+RsElKuXcoNxFC/ARYD3RLKeecdX0t8F3ACPxISvnsRQOR8nXg9XgX1H8FtkopG4BH4m0bJj4pSdt9jIzff4C0mun+4t0E5lcnLZzWYCdb+97lsP8kTqOdzxTdzs25V6tWBGNARPzY+w4kdp9ZvHrXzUhaHp6yaxPbnKOOYW1GUsaJmCY53uVJTJ8d6fQQ1SQWk4GrpmVz1+Iy6qvymF2ciTEJLeyHmniOSCmfiv/5FSHEC8C/o7dJGIoXgOeAn52+IIQwAt8HbgBagZ1CiDfQk9Az573+YSlld/zP34y/blIRwRCZv9lB2qFTBGZNp+9ztxHLTk6Nq7ZgJ1td2znkO4HDaOfeKRtYl7daJZzRJGNY+0/g6P74AlUCltJX+yCe0qtVNedJQkpJ60CAPc36iGZ/mxt/OIYAZpdk8LlVFdRX5bF4ejY2c/KnU4eaeHqFEIuklB8BSCmPCSHyh/pQKeU7Qojp511eApyIj1oQQrwM3CalfAZ9dHQOoY8HnwV+J6X8eKgxjGfm5i6yXnsLo8dP/51rGFy7HJLwrqU91M2bfds54DuGw5DGpinrWZ93LXbVlmBUmHydeqLp3hOvEuBFIgjmzqan7vPxbc6LkCaV8CcDly/M3tYB9rQMsK91gF5vGICynDRum1/Cquo8llfkku1IvTXVoSaeLwIvCyE+AvYDdcCpEYqlBH2TwmmtwNJLfP6TwPVAphCiSkr570KIXPSOqAuEEF+PJ61zCCE2A5sBCvKyRij0MaJpOLfvxfn2bqK5mXR+/UHClcPaTn9FOkM9vOnazj7vUdIMNj5deAu35q/BYVTrBSPJEPbGqwTsJq1nLxZfBxCvElC+Hm9JPd6SemK2nCRHqowFfzjKgbZB9sanz5pc+jpNlt3Myip9jWZlZR5Tc1P/53CoHUj3CiHmo//CnwNsY+TO7lzoLftFt9xJKb8HfO+8a33AY5d6iJTyeeB50He1DT3M5DC4vWS9/jbWpk68y+bguv9mZNrYLtZ3hXoTCcdmsHBX4c3cmncdTlPq/0MfF7QINtfRxKYAW/9xvRma2YGvaAW9dY/jKVlFOLNCTZ9NAtGYxtEuT3xDgJujXR5imsRqMrCkPId7l05lZVUetUUZGJIw4zEcV9KBNAT8Nv7fSGoFys76uBRoH+FnjEu2w41k/mY7aJLez23At6JuTJ/fHe7jzb4d7PUexmqwsLHgRjbkX0+6yTGmcUw4UmL2tsY3BOzG3rv/rCoB8+LN0K5WVQImCSklzS5/YvrsQNsggYi+TjO3NJPPX11BfXUeC6emxjrNcKRSUaydQLUQohxoAzYB9yY3pCSLRMnY+iGOj44QmjaF3sc2Ei0cu2mVnrCLP7h2sNtzCIvBzO35N3BbwfVkmJxjFsNEYwy5SYsX2LT37MEc6AUglDGN/hl34S25Gm/xclUlYJLo9YbY2zKQ2H3W748AMD3Xzh2LSqivymN5RR6Z9on1xiMpiUcI8RKwGsgTQrQCT0kpfyyEeAJ9q7YR+ImU8mAy4ksFpi4XWa9uw9wzgHvtcgY2rgbT2LzL6Q338wfXDj72HMQizNyWfx23FdxApkl1hhwqEQtj6zuE/fTus9NVAiyZ8TWa01UCyj7hTspE4AtF2d925jxNS38AgByHhVXV+dRX5bGiKpfS7Ik9fZ2UxCOlvOci17cAW8Y4nNQiJfZdh8l4889odhtdX7mH4JzKMXl0X2SAP7p28NHgAYzCxPr8NXwq/wayzOrd92WTGpbBxsR5mrS+gxi0cLxKwCI6q+9UVQImkUhM42infp5mT8sAx7s8aBLSzEaWlOfwwPLp1FfnUVOYPu7WaYYjlabaJj3hD5L1xrvYjjUTmFtJ7yMb0DJGfx3FFRngj6732TW4H6MwcHPeaj5VcCM55tTo2ZPqjIHeeN2z3djPrhKQNQNX7f1441UCNLNaE5voNClp6vMnps8OtrkJRjUMAuaVZvGFa6tYWZXHgqlZWMdoBiMVqcSTIiyn2sl6/W0M/iCuTTfguX7JqJ/N6Y+4+ZPrfXYO7kMIwdq8q9lYcCM55nG2zXyMiWiAtN4DiXUaq6cZiFcJKL02MYWmqgRMDt2eoJ5oWtzsaxtgIL5OU5nv4O6rylhZlcfSilwy0ybWOs1wqMSTbDGN9Lc/xrFjL9GCHLq+ci+RaaP7C2sgMsif+j/gz+49IAQ35NazseAm8iwX7xo5qckY1oET+qimew9prsMIGUUzWvFNWYpr1gN4S1YRzJmptjlPAt5glP1tA+yOb3NuH9DXafKcFlbPyKe+Op+VVbkUZaqD1BejEk8SGfs9ZL32FpbWbrz183DdexPSNnqnjN1RD9tc7/Ph4F4kcF3uCu4oWEu+RR1APJ/J15lo72zv2Ycx4gEgkDubnrmP6lUCCherKgGTQCSmcbhjMNGf5kS3F02C3WJkWUUuj9SXU1+Vx4xC5ycW2FR0KvEkie3ASTJ/+x4I6HnsU/iXzB61Zw1GvWzr/4AP3LuRUnJtznLuLFxLgSV31J453hjCXtJ69yWmz85UCSjCXX5LvMjmSmJp6u9sotOk5FSvLz59NsChjkFCUQ2jQbCgLIsn11RTX53HvNIsLKbkt5Efj1TiGWMiHCHjd+9j33ucYGUJvZ//FLFRKt3jifp4q/8D3nfvJiY1rs1Zyh0FNzPFmtw+PSlBi2JzHUl03rT1H9OrBJjs+IpX0Fv3GN6SVYQyK9X02STQNRhkTzzR7G9z4w7o6zTVBU7uWTKVVdV5LCnPId2m1mlGgko8Y8jU0Uv2L7dh7B9kYH097tuuBuPIv2PyRv281f8B77k/JipjrM5eyp2FaymyFoz4s8aNRJUAffdZWu8BjFE/UhgI5J2uErAKf8ECVSVgEhgMRNjf5k5Mn3W4gwAUpFu5blaBXvesKo/CDDWVOhpU4hkLUuL44ADpf9pFLN1B119/htDM6SP+GF/Mz9v9f2bHwEdEZJSrs6/irsKbKbYWjvizxoMzVQL2xKsE9AAQSp/GQPUd+vRZ0XI0q9o2PtGFojEOd3gSieZktxcJOKxGllfkJcrRVOardZqxoBLPKDN4A2S+8Q62E634F8yg76H1aM6RPZXsjwUSCScsI9RnLeauwpsptU2u7bwiFsbmOhTffbYbm/skoFcJ8JWsxJOoEjA1yZEqoy2mSRp6vIlSNIc6BonEJCaDYMHULL58/Qzqq/OoK83EPAqzDsqlqcQziiwnWsn61TuIcIS++2/Gu3rhiK4X+GNB3hn4M9sHdhHWIqzIWsinC2+hzFY0Ys9IVTbXYdJ69hNJy8MUGtCnz/oOYYiF0Awm/AWL6azaqFcJyJurqgRMcFJKOuPrNHtbBtjX5sYTjAJQMyVdrxBQpa/TOKzq116yqe/AaIjFSP/jLpwfHCBcnEfvYxuJlI7c+kogFuTdgZ28O7CLoBZieeYCPl24jmlpxSP2jFRlCA+SdfI35Bx7GaSW6KURzKrGNfM+vUpA0TJVJWAScAci7IuXotnbOkDXYAiAKZk21s6eQn11HssrcylIV+s0qUYlnhFm7HOT9epbWDp68Vy7iP67r0daRmaxOhgL8e7ALt4d2ElAC7I0cx53F65jelrpiNw/JUmJ1d2AvWsXjq5d2FxHEWhI9AZOEkFP3eN0LflasiNVRlkwEuNQ+2Bi+qyh1wdAus3E8opc6qv1ZmjleQ61TpPiVOIZKVKStvc4Gf/zPtJkovuJuwgsrBmRWwe1EDsGPuKd/j/j14Isyajj7inrKE+bmBWNRcSPvWcPjniyMQVdAPjz5tG94ItEnCUUv/d/gxZBGsx4pt2Q5IiV0RDTJCd7vIltzkc64+s0RsHiadlsXFjCyqo85pZkYlLrNOOKSjwjQATDZG7ZQdqBBoI1U+l99HZiOcOv6BzSwomE49MCLEqfw6Yp66m0T7DF8fhWZ0fXLhydu0jrO4iQUWLmdDyl1+ApuxZv6Wqi9vzES0JZVTg63sdXtBx/4aIkBq+MFCkl7QPBxIhmf5sbb0hfp6ktyuChleWsrMrjqunZ2C3qV9d4pr57w2Ru7Sbr1bcwDnoZ+NQ1uNetBMPw3n2FtTDvuXfzVv+H+GJ+FqTXsmnKeqrt00cm6BQgYiHSevfj6NyJo+sjzP5OAALZNfTMfRRP2Rr8hQsveqbGX7hIJZxx6kjHIPvb3MwtyaQw05boTbO3xU2PV1+nKc6ysW5uESur81hRmUuec2zbvCujS0gpkx1D0tRUlMkf/OOXruzFmobjvf2kv/URsewMej9/O6Gq4U19hbUI78cTjjfmY55zFpumrKfGUT6s+6YKk68rPn22k7SefXqfGlMa3uJ6PGXX4im7loizJNlhKldAalGktw9ikUt+Xjim4fZHzlqj0xkEWE1GrGYDNpNBTZ2lCJvNRmlpKWbzuW8AhRAfSSkXX+l91YjnChgGfWS9/jbWxg58V9XS9+AtSPuV75yJaBE+cO9hW/8HeGI+6pwz2TRlHTMdY9MAbtRoEdL6DuPo2oW9aydWTwugt3l2zfoMnrI1+KYsUYU2JwDp7SMvJ4us7OxzFvalhEAkRiAcIxCJISIxzt7fmW41UZhpI81sVBsCUoyUkr6+PlpbWykvH9k3vyrxDJH1aDOZv34HEdXofWg9vvp5V3w2J6JF+XBwD9tcHzAY8zLHMYNNU9ZT66wa4ajHjjHoiieaXdi792CM+tEMZnxTltE3+2E8ZWsIZ06MEZxylliErOxsQCSSTCAcIxjRkEgEgjSLkWy7hYFABKRECEFBhk2t16QoIQS5ubn09PSM+L3Vd/xyRaNkvLkTx85DhMsK6XlsI9GiK6tUHNWi/HlwH3/qfx931EOto4pNU9YzxzljhIMeAzKGrf8Y9s5dOLo+wuY+AUDEPoWBqtv1UU3xSnWuZoKKxDSOdXmYZozSPhAkGNWQUiIAm8VIXroFp9WE3WLCGG9smBOK4gtFcVhN6jBnihutUei4/a4LIW4H1gEFwPellFvj1x3AO8BTUsrfjMSzTD39ZP1yG+bufgZvWEL/nWvAPPS/uqiMsXNwH39yvc9AdJBZ9kq+MvVh5jhnjKtpBkN4EEfXx/rZmu6PMYYHkcKIv2ARnTP+hsGy6whl16iqzhNQTJOc6Payr22Afa1ujnQMEoxq/GhDEXlArkNPNA6rEeNFNtmMZsJ57bXX2LhxI4cPH2bmzJlXfJ/PfvazrF+/njvvvPOin/Ptb3+bv/3bv018vGLFCt57770hP+vpp5/G6XTy1a9+9YpivVJSSr70pS+xZcsW7HY7L7zwAgsXLhyTZycl8QghfgKsB7qllHPOur4W+C5gBH4kpXz2YveQUr4OvC6EyAb+Fdga/19/A/z3iAQqJfaPj5K+9UOkzULXlzcRrBv6NFhMxtg1eIA/ut6jP+pmhr2cL059kHnOmeMj4UiJxX0KR9fOcw5xRm05DE67AU/ptXhKr0azqpbZE01M03vT7I8nmkMdg/jDMUBvGbBpyVSWVeRSJFxUF6YnOVp46aWXqK+v5+WXX+bpp58e1Wedn3iuJOkk0+9+9zuOHz/O8ePH+fDDD3n88cf58MMPx+TZyRrxvAA8B/zs9AUhhBH4PnAD0ArsFEK8gZ6Enjnv9Q9LKbvjf/5m/HUIIa4HDgHDXq0WgRCZv95O2pFGArXl9H1uA7Gsof1gxWSMjwYP8Mf+93FFBqhKm8YXyj7DgvTalE845x7i/AhTsA84c4jTU7aGQH4dCLX7aCLRpKS5z8++Njf7Wgc42D6YOEtTnufg9gUlrKjMZVnFuVucDx/uH/KzPmrq54OGPpZV5LJo2vDbrnu9Xnbs2MG2bdvYsGFDIvG89dZbPP300+Tl5XHgwAEWLVrEL37xC4QQfOtb3+LXv/41gUCAFStW8MMf/vCcn80//vGPPPfcc7z22msAvPnmm/zgBz9gxowZBAIB5s+fz+zZs3nxxRdxOp14vV4AvvOd7/Dzn/8cg8HAzTffzLPPPst//Md/8PzzzxMOh6mqquLnP/85dvvlFQx++umnaW5upqGhgebmZr785S/zxS9+cVh/X7/61a944IEHEEKwbNkyBgYG6OjooKho9Gs9JiXxSCnfEUJMP+/yEuCElLIBQAjxMnCblPIZ9NHROYT+r+NZ4HdSyo/jl68FHEAtEBBCbJFSakONz9LUQdZrb2PwBui/6zoGb1qm7/e8TDGpsdtzkD+43qMv0k9FWhmPl97LwvTZqZtwruAQpzL+SSlpHQiwv1VPNAfaBxNN0Eqz07hl7hRWVOaxrCKXKZmX937u7399kEPtg5f8HE8wwpFOD5rUf7RmTkm/ZJO12uIMnrr10l16X3/9ddauXcuMGTPIycnh448/Tkwd7d69m4MHD1JcXMzKlSvZsWMH9fX1PPHEE/zd3/0dAPfffz+/+c1vuPXWWxP3XLNmDV/4whfo6ekhPz+fn/70pzz00EPceuutPPfcc+zZs+cv4vjd737H66+/zocffojdbsfl0itvbNy4kUcffRSAb37zm/z4xz/mySefvOTXdLYjR46wbds2PB4PNTU1PP7443+xzfnuu+/m6NGjf/Har3zlKzzwwAPnXGtra6Os7MwRkNLSUtra2iZu4rmIEqDlrI9bgaWX+PwngeuBTCFElZTy36WU3wAQQnwW6L1Q0hFCbAY2AxSc3/lT03C+swfnu3uI5mXR/aW7CZdffuFNTWrs9hziD6736I24KLeVsrnkbhZnzE3JhDPcQ5zK+HO6ivO+Vjf7Wt0caHfj8oUBmJJh47qZBSyvzGV5ZS6l2SPbvuNsg8EoWvwQjyb1j4fb3fOll17iy1/+MgCbNm3ipZdeSiSeJUuWUFqq1zScP38+jY2N1NfXs23bNr7zne/g9/txuVzMnj37nMQjhOD+++/nF7/4BQ899BDvv/8+P/vZz/7y4Wf5wx/+wEMPPZQYzeTk5ABw4MABvvnNbzIwMIDX6+Wmm24a0te3bt06rFYrVquVgoICurq6El/Taf/1X/912fe70BnOsfo9lUqJ50Jf8UVPt0opvwd87yL/74VLvO554HnQD5Cevm4c8JD12ttYWrrwrpiL6761yLTLOy2tSY29nsO86dpBT8TFNFsJj5RsZknGvJRLOBc/xLmS7rIvqkOcE1C3Jxgf0bjZ33amOkCe00J9lV7BeXlFLtNy7SPy7/WTRiagT7Pd96MPiEQ1zCYD3920YFjTbX19ffzpT3/iwIEDCCGIxWIIIfjOd74DgNV65mfZaDQSjUYJBoP81V/9Fbt27aKsrIynn36aYDD4F/c+PcKx2WzcddddmEyX/rUp41vFz/fZz36W119/nXnz5vHCCy/w1ltvDelrvNDXcL6hjHhKS0tpaTnzXr+1tZXi4rGpcJ9KiacVOPvofynQPhYPth06ReZvtoOE3kdvw7d87mW9TpOSfd4j/MG1na5wH2W2Iv538aMsyZyHIVXWPrQoaX2HEmdrrJ5mQB3inMhcvjD7WvWeNPtb3XQO6r9Ms+xmVsSTzPLK3KR221w0LZsXP7dsxNZ4XnnlFR544AF++MMfJq5dc801bN++/aKvOZ1k8vLy8Hq9vPLKKxfcxVZcXExxcTH/+I//yJtvvpm4bjabiUQifzHddeONN/Ktb32Le++9NzHVlpOTg8fjoaioiEgkwosvvkhJyV++wXvuuecAeOKJJ4b2FxA3lBHPhg0beO6559i0aRMffvghmZmZYzLNBqmVeHYC1UKIcqAN2ATcO5oPNHj8ZL/8JrZjzYTKi+n9/O1EC3I+8XWalBzwHuVN1w46wz2UWqfw1WmfY1nm/JRIOJc+xPmQOsQ5wbgDEfbHNwPsb3PT2h8A9HYBS8tz2Xx1Bcsrc6kpTMcwhLXK0bZoWvaIbCoAfZrta187tzXGHXfcwX/+539y9913X/A1WVlZPProo8ydO5fp06dz1VVXXfT+9913Hz09PdTW1iaubd68mbq6OhYuXMiLL76YuL527Vr27NnD4sWLsVgs3HLLLXz729/mH/7hH1i6dCnTpk1j7ty5eDyev3jOkSNHWLly5VC//Ctyyy23sGXLFqqqqrDb7fz0pz8dk+dCkmq1CSFeAlYDeUAX+pmbHwshbgH+X/SdbD+RUv7TaMYxx5Ym/3v6dHzL59L30HowXbpLpSYlB33HeLNvBx3hbkqshdxduI7lWQsxJjPhXOIQ5+DUNeoQ5wTjDUbZ334m0TT1+QGwW4wsKc+Jj2ryqC3OSBzaHG2HDx9m1qxZY/KsZHjiiSdYsGABjzzyyKg+Z/369bz66qtYLJZRfc5QXOh7Oy5rtUkp77nI9S3AljENRggixXmXTDpSSg76jvOmawftoS6KLAV8eepDrMxalLSE80mHOD2lawjmzFSHOCcAfzjKgbZB9rfpiaahx4cEbGYDi6Zl8+nFZSyvzGVuSSZmVVxzxC1atAiHw8G//du/jfqzfvObETnznvJSaapt7AmQJiPBmmkX/N9SSg77TvKmazutoU6mWPL5YtmDrMpejFFcenQ04tQhzkkjGIlxqGMwvhlggBPdXjQJZqNg0bRsNswrYXllLvPKMrF+wihdGb6PPvoo2SFMOJM68USz0un6688Qrjp3S6KUkqP+Brb2bacl1EGhJY8nyu7nmuwlY5pw1CHOySEUjXGk08P++K6zY10eoprEZBDML8viC9dWsbwil4XTsrGZVaJRxr9JnXi0TOc5SUdKyTH/Kba6ttMcbCffnMNfld7H6pxlmMYi4ahDnJPC6cKap7c3n27pbBAwtySTz63SNwMsnpatimgqE5L6V42ecI4Hmniz710ag23kmbN5vPReVmcvw2wY3b+iM4c4d+Ho2pU4xBlUhzgnjERhzfgW59OFNQX6ifwHl09neWUuV5XnkDHMQ5SKMh5M6sQzGPOyvf8j9nkPcyrYSq45i8+XbGJNzopRTThnDnHuIq13H4ZYSB3inEBOF9Y8vevs7MKaMwrPFNZcVpFDlj11di8pyliZ1InHHfXwq943STc6eLTkbq7PWYF5NEYWlzrEOfM+dYhznDtTWFOv4Hx2Yc2KPAefWqBvBji/sKYy8lRbhMv34osv8s///M8AOJ1OfvCDHzBv3v/f3v0HR13feRx/viEhISEgBEKp/Eqq/CpI+CnSqrTWFn8MInhXWu6uFiltqR3t3cxhp50p50yttt5MD7TXirH1PLVWpljqwZwOlmrvWiVWof6gXoutjVqRGAgQICG874/vrEt3iAAAEkFJREFUN2EJ2U02+93vbpLXY2Ynu99897uvXb7LO99f78/MWF67XxceCPr0XDXyI1wx8tJIl5v6Is7rOTzuMl3E2Uu5O3UNx8LOAGc21hw3YjBXzRjT3u9s9FD9MREnDYvQfZWVlfzyl79k+PDhbN++nTVr1vT5YRHygmEUWiEzy3r+l1E7b6W44f8oCRtuaiTOvsPdefvQ8bA7QHCKc0NTUGjGDCvmsqkV7W1ostlYs8/5y3Pwp2dg4sUwbn7Gi9OwCOlZuHBh+/0FCxZQV1eX0fLS0a8Lz4jCYfzLB25icmlVj54/oLmRkv0vBN2ddRFnn7K/8Xh7r7PExpqjhhRx8fmjIm+s2adsvwX++rvU85xohHdeAj8VXA4wejoUDU0+//tmwBVJx4UENCwCpNckNFFNTQ1XXHFFt7Nkql8XnuEFw9IrOrqIs8+qP3Ii2KLp0FhzeElhe5HJdWPNPuX4oaDoQPDz+KHUhacbNCxCek1C2/ziF7+gpqYmZUPVqPXrwtMdyS/ivEAXcfZiB5ua+d2bh9pvbY01hxYXcGFVOWuq8rOxZq/QxZYJEOxmu38JtDbDwEGw/N6MdrdpWIRAuls8e/bsYfXq1Wzfvp3y8vK08mRChaejxIs436ll8AFdxNkXHD7ewkttWzQdGmteWDmCf7hoQuyNNfu1cfPhM1sjO8ajYREC6WzxvPHGGyxbtowHHniASZMm9ej1ekqFhzMv4izZ/zyDjr4N6CLO3ixVY825E0aosWY+GDc/kpMKQMMi9MStt95KfX09a9euBaCgoIDa2tpYXjsnwyLki9kThvr//OM0ihr3nXER5+FxH9VFnL3M8ZZWXnmrMdyiOd1Yc1DBAGaPP4eLqkaqsWYWaViEaGhYhH5gQMtRihv2cqjyShomf0oXcfYie+oOsvO1dzl1KjjVWY01JVs0LEL0+nXhAcAGcLz8gxwZe0mukwhBA82GpmYajrbwXlMzDUeb2382NDXT0NTC/sbjNB4/fWD1vFGlrL64ioUfKGfuxOGUDNJqLdHRsAjR6+ffUMMHFHJ0zEW5DtLnHW9p5b2weJz+2XK6sITFJbGgtDGD8tJBjCor4txzBlNcMICX32rEgYEG184ey5c+cl78b0pEeqRfF56W0tG8fuUDNI2ek+sovZK7c/RE61lbJmcWmBYamprbm2QmKhhojBpSREVZEZNGl1FRVkRFWTEVQ4vOuF9eOoiChBMAnv9zAyvv/Q0tJ09RWDCABVXxnQYqIpnr14Xn5OAKFZ1OtJ5yGo+dvavrvaaWs3Z9tbSefXLK4MKBjCoromJoEdXlpe33K8qKg4IS3j9ncGGPrpGZM2E4D65ewG/21bOgqpw5E4ZH8bZFJCb9uvD0N6mOn7yXcAzlYFMzpzo52XHY4EJGlRUxemgR08YMpaKsiFHhLXFLZUhRQdav7p8zYbgKjkgv1WsLj5ktBa4CKoC73f0JM7sYWEnwvqa5+8JUy+grjjW3nnHspCfHTyrKihk3ooQ5E84sIqPCrZRRZUU6O0zynoZFSN+uXbtYsGABjzzySMr3G6WcFB4zuw+4Gtjv7tMTpi8G/g0YCNzr7kl7b7j7Y8BjZjYcuBN4wt2fAZ4Ji9KubL6HbHN3jpw42X6c5OwD880cDKcfa4nu+IlIb6ZhEdLT2trKunXr0u4bl6lcbfH8CLgLaO+2Z2YDgbuBy4E6YJeZbSUoQt/q8PxV7r4/vP/18HmJPg2sjj525jI9flIyKDx+UlbExJGlp4tIh+Mo55QUqpml5LUX979I7Tu1zB09l+qK6oyXp2ER0rdx40aWL1/Orl3x/p2ek8Lj7k+b2cQOk+cDf3D3fQBm9mPgGnf/FsHW0RksWDtuB7a7+28Tpo8HDrl7Y5bid6ql9VSXhaQnx08qhoZFJbw/pKjX7h2VfuKO5+5g73t7U85zpPkIv2/4PY5jGJOHT2bIoCFJ558yYgrr5q9LuUwNi5Bek9A333yTLVu28NRTT/WPwpPEucBfEh7XARemmP/LwMeAYWZ2nrt/P5x+A/DDZE8yszXAGoAx547rMlTH4ydnnSocFpfDnRw/GWBQXhpsiYzX8RORdodbDuMEf4E5zuGWwykLT3doWIT0moTefPPN3HHHHQwcGP//PflUeDrbL5S0kZy7bwA2dDL9G6lexN3vAe4BqKia5lt3v8WQooK0jp8Uth0/GVrM5PeVdXqqcEVZESN0/ET6oa62TCDYzfa5Jz5Hy6kWCgcUcvvFt2e0u03DIgTS2eKpra1lxYoVABw4cIBt27ZRUFDA0qVL08rVE/lUeOqAxE2QscBb2XzBQ8da2PTMvvbHOn4iEo/qimo2fXxTZMd4NCxCIJ0tntdff739fttZfHEUHcivwrMLON/MKoE3gRUEJwlk1QCDVR+q5ObLJ+n4iUiMqiuqIzmpADQsQm+Tk2ERzOxhYBEwEngH+Ia715jZlcB3Cc5ku8/dv5nNHEVjzvfK1Rt4cPUCXYwokgENixANDYuQRe7+qSTTtwHb4soxemixio6IpKRhEaLXr/ctVZQVqeiISEoaFiF6OuVKRERipcIjIiKxUuEREZFYqfCIiEisVHhEpM/YsmULZsbeval7xXXl+uuvZ/PmzSnnue222854vHBhz0ZhWb9+PXfeeWePnpupnTt3tjc6vfTSS2N7XRUeEekzEodFyLaOhae3DYtw8OBB1q5dy9atW3n55Zd59NFHY3ttFR4RyYmmF17gwA/uoemFFyJZXtuwCDU1NWcUnp07d7Jo0SKuu+46pkyZwsqVK2m7cP7WW29l3rx5TJ8+nTVr1tDxgvodO3Zw7bXXtj9+8sknWbZsGbfcckv7sAgrV64EYMiQ001Ov/3tbzNjxgxmzpzZ3lFh06ZNzJs3j5kzZ7J8+XKampq6/d7Wr1/PqlWrWLRoEVVVVWzYcFabyrQ99NBDLFu2jPHjxwNQUVGR8TK7q19fxyMi0fvrbbdx4tXUu7pajxzhxN694A5mFE2ZwsAhybtTF02dwvsSBl3rjIZFSK9J6GuvvUZLSwuLFi3i8OHD3HTTTWfNky0qPCISu1ONjUHRAXDnVGNjysLTHRoWIb0moSdPnuT5559nx44dHDt2jIsuuogFCxYwadKktHL1hAqPiESqqy0TCHazvfHZVXhLC1ZYyPvv/A4ls2b1+DU1LEIgnS2esWPHMnLkSEpLSyktLeWSSy5h9+7dKjwi0jeVzJrF+B/eR9NzuyiZPy+jogMaFqFNOls811xzDTfeeCMnT56kubmZZ599lq985Ss9et106eQCEcmJklmzGPn5NRkXHQh2syWeBACnh0VIJnFYhKVLl3Y5LMK4ceM6HRah7eSCNosXL2bJkiXMnTuX6urq9lOl24ZFuPzyy5kyZUqnr7N3717Ky8u7fL9RmDp1KosXL+aCCy5g/vz5rF69munTp8fy2jkZFiFfzJ0712tra3MdQ6TX07AI0dCwCCIiomERskCFR0QkBQ2LED0d4xERkVip8IhIJPrz8eK+Klv/pio8IpKx4uJi6uvrVXz6EHenvr6e4uLiyJetYzwikrGxY8dSV1fHu+++m+soEqHi4uKzuiNEodcWHjObCtwEjAR2uPu/m1kp8D2gGdjp7g/mMqNIf1FYWEhlZWWuY0gvkZNdbWZ2n5ntN7OXOkxfbGa/N7M/mNktqZbh7q+6+xeAvwXazidfBmx2988BS7ISXkREMpKrYzw/AhYnTjCzgcDdwBXANOBTZjbNzGaY2eMdbhXhc5YAvwJ2hIsZC/wlvN8aw/sQEZE05WRXm7s/bWYTO0yeD/zB3fcBmNmPgWvc/VvA1UmWsxXYamb/BTwE1BEUnxfRiRMiInkpn47xnMvprRUIisiFyWY2s0UEu9aKgG3h5J8Cd5nZVcDPkzxvDbAmfHii4+6+FIYBh7r5u1SP2+53nEb4eCRwIAuZkuXoKt8woDCNTKlyRZkpqs+qq0zdzRJlpq4y9MZMbdPIcqZkObrKp/U8vXVqcjfzdM7dc3IDJgIvJTz+G+DehMd/D2zMcobaNOa9p7u/S/W47X7HaQnTs5IpWY6u8oXZup0pVa4oM0X1WXWVKY0sWf336+2ZEtYlref9bD3v7JZPWzx1wLiEx2OBt3KUpTOdbkEl+V2qxz9PMS2bmVK9Zqp8PwdmR5Qrykzp6mmm7maJMlNXGXpjpp7mSTdTqtfUep56WrbWqbPkrDt1eIzncXefHj4uAF4DLgPeBHYBn3b3l7OYodYz6LCaDcrUffmYS5m6R5m6Lx9zZZopV6dTPwz8GphsZnVmdoO7nwRuBP4beBX4STaLTuieLC+/J5Sp+/IxlzJ1jzJ1Xz7myihTvx6PR0RE4qdTjkVEJFYqPCIiEisVniTMbICZfdPMNprZZ3KdB4Jrl8zsGTP7fngdU14ws1Ize97MOr3QN25mNjX8jDab2RdznaeNmS01s01m9jMz+3iu8wCYWZWZ1ZjZ5hznKDWz+8PPZ2Uus7TJl88mUZ6uQ2l/3/pk4YmiFxxwDcFFrS0Ep3rnQyYHjgDFeZQJYB3wk0zzRJXJO+/jlw+5HvOgj+D1wCfzJNM+d78h0ywR5Iulz2I6mbL52WSQKdJ1KKJM6X/fMrkIKF9vwCUE5+QnXqA6EPgjUAUMAnYT9ISbATze4VYB3AJ8Pnzu5jzJNCB83mjgwTzJ9DFgBcEX4ep8yBQ+ZwnwvwSn5OfFOpXwvH8FZudZpozX8QzzfRWoDud5KOosPcmUzc8mgkyRrENRZUr3+5ZPF5BGxiPoBWdmdQTDK0AEDUejyJSggaBVUM4zmdlHgFKC/zyOmdk2dz+Vy0zhcjr28ctIRJ+VAbcD2939t/mQKZvSyUdMfRbTzPRKtnL0NJOZvUqE61AUmYBX0v2+9cnCk0RaveAI+r5tNLOLgafzIZOZLQM+AZwD3JUPmdz9a2G264EDmRSdqDJZ5338siHdderLBFuIw8zsPHf/fq4zmVk58E1glpl9NSxQ2ZQs3wa66LMYd6YcfDZdZiKedSitTD35vvWnwmOdTEt6EZO7NwHZ3r+bbqafEhTEbEorU/sM7j+KPkq7dD+nncDObIVJkG6uDQT/wWZTupnqgS9kL85ZOs3n7keBz8aYI1GyTHF/NomSZYpjHUomWaadpPl965MnFySRj73glKl78jET5GeufMyUKB/zKVP3RJapPxWeXcD5ZlZpZoMIDohvVSZlykA+5srHTInyMZ8yxZ0pm2dr5OoGPAy8zelToW8Ip19J0Ij0j8DXlEmZenOufMyU7/mUKT8yqVebiIjEqj/tahMRkTygwiMiIrFS4RERkVip8IiISKxUeEREJFYqPCIiEisVHhERiZUKj4iIxEqFRyTHzGyimR0zsxcTpo02s4fMbJ8Fo7v+2syu7WI5O83sEx2m3Wxm3zOzwWb2opk1m9nIbL0Xke5Q4RHJD39092poH7fnMeBpd69y9zkEfbHGdrGMh8P5Eq0AHnb3Y+Hyc91oUkSFRyRTZvaomd1lZr8ysz+b2YfN7D/M7DUzq+nBIj8KNHvCWCvu/md335jwmn9nZs+FWzE/MLOBwGbgajMrCueZCLwf+FUm708kaio8IpmbAexz9w8D9wM1wDpgOrCsrRCk4YNA0tElzWwq8EngQ+FWTCuw0oPxY54DFoezrgAecTVklDzTnwaCE4mcmRUTjAj73XDSMaDG3d8Of9/E6SHUe/oadwMfJtgKmgdcBswBdgV75RgM7A9nb9vd9rPw56pMXlskG7TFI5KZDwK/9dNDfs8EngUws7aBsi4wsy1mVmhm/9mNg/svA7PbHrj7lwiKzahwkgH3u3t1eJvs7uvD3z0GXGZms4HB7p50y0kkV1R4RDIzA9id8PgCYE94fyawx913Ay8CPwS+6+4HuljmU0CxmX0xYVpJwv0dwHVmVgFgZiPMbAKAux8hGIb4PoKtH5G8o8IjkpkZBEWlbbfbYHdvCH+XWIQagFfdvbarBYbHZJYCl5rZ62b2HMGxo3Xh718Bvg48YWZ7gCeBMQmLeJig6P04w/cmkhUaCE4ky8xsGvBPQL27/3Mnv58IPO7u02PI8idgbje2ukSyRls8IllkZsOArwFrgUozu7CT2VqBYYkXkGYhx+Bw+YXAqa7mF8kmbfGIiEistMUjIiKxUuEREZFYqfCIiEisVHhERCRWKjwiIhIrFR4REYmVCo+IiMRKhUdERGKlwiMiIrH6fxb3Aj6IoiHXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "f,ax = plt.subplots()\n",
    "\n",
    "ax.plot(masses, sigmaout[0], marker='.', label='Analytical, n = 0')\n",
    "ax.plot(masses, sigmaout[1], marker='.', label='Analytical, n = 2')\n",
    "ax.plot(masses, sigmaout[2], marker='.', label='Analytical, n = 4')\n",
    "ax.plot(masses, sigmaout[3], marker='.', label='Analytical, n = 6')\n",
    "#ax.set_title('DM-electron scattering (n=0, Vrel=0)')\n",
    "ax.set_xlabel('$m_\\chi$ $[\\mathrm{GeV}]$')\n",
    "ax.set_ylabel('$\\sigma_0$ $[\\mathrm{cm}^2]$')\n",
    "ax.set_xscale('log')\n",
    "ax.set_yscale('log')\n",
    "ax.set_xlim([1.0e-6, 1.0e3])\n",
    "ax.set_ylim([1.0e-30, 1.0e-8])\n",
    "ax.fill_between(masses, sigmaout[0], 1.0e-8, alpha=0.2)\n",
    "ax.fill_between(masses, sigmaout[1], 1.0e-8, alpha=0.2)\n",
    "ax.fill_between(masses, sigmaout[2], 1.0e-8, alpha=0.2)\n",
    "ax.fill_between(masses, sigmaout[3], 1.0e-8, alpha=0.2)\n",
    "plt.legend(loc='lower right')\n",
    "#plt.savefig('dme_analytical.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Analytical:\n",
      "TeV n2/n0 = 3789.4709190746717\n",
      "TeV n4/n2 = 2681.812609453029\n",
      "TeV n6/n4 = 2256.074066496864\n",
      "\n",
      "Numerical:\n",
      "TeV n2/n0 = 11037.234042553191\n",
      "TeV n4/n2 = 2771.084337349397\n",
      "TeV n6/n4 = 9391.304347826088\n",
      "\n",
      "Percent errors:\n",
      "191.3%\n",
      "3.3%\n",
      "316.3%\n"
     ]
    }
   ],
   "source": [
    "print(\"Analytical:\")\n",
    "print(\"TeV n2/n0 =\", sigmaout[1][3]/sigmaout[0][3])\n",
    "print(\"TeV n4/n2 =\", sigmaout[2][3]/sigmaout[1][3])\n",
    "print(\"TeV n6/n4 =\", sigmaout[3][3]/sigmaout[2][3])\n",
    "\n",
    "n0 = 3.76E-21\n",
    "n2 = 4.15E-17\n",
    "n4 = 1.15E-13\n",
    "n6 = 1.08E-09\n",
    "\n",
    "print(\"\\nNumerical:\")\n",
    "print(\"TeV n2/n0 =\", n2/n0)\n",
    "print(\"TeV n4/n2 =\", n4/n2)\n",
    "print(\"TeV n6/n4 =\", n6/n4)\n",
    "\n",
    "print(\"\\nPercent errors:\")\n",
    "print(format(abs(((n2/n0)-(sigmaout[1][3]/sigmaout[0][3]))/(sigmaout[1][3]/sigmaout[0][3])), \".1%\"))\n",
    "print(format(abs(((n4/n2)-(sigmaout[2][3]/sigmaout[1][3]))/(sigmaout[2][3]/sigmaout[1][3])), \".1%\"))\n",
    "print(format(abs(((n6/n4)-(sigmaout[3][3]/sigmaout[2][3]))/(sigmaout[3][3]/sigmaout[2][3])), \".1%\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Analytical:\n",
      "TeV n2/n0 = 3789.4709190746717\n",
      "TeV n4/n2 = 2681.812609453029\n",
      "TeV n6/n4 = 2256.074066496864\n",
      "\n",
      "Numerical:\n",
      "TeV n2/n0 = 4764.638346727899\n",
      "TeV n4/n2 = 2771.084337349397\n",
      "TeV n6/n4 = 5982.608695652175\n",
      "\n",
      "Percent errors:\n",
      "25.7%\n",
      "3.3%\n",
      "165.2%\n"
     ]
    }
   ],
   "source": [
    "print(\"Analytical:\")\n",
    "print(\"TeV n2/n0 =\", sigmaout[1][3]/sigmaout[0][3])\n",
    "print(\"TeV n4/n2 =\", sigmaout[2][3]/sigmaout[1][3])\n",
    "print(\"TeV n6/n4 =\", sigmaout[3][3]/sigmaout[2][3])\n",
    "\n",
    "n0 = 8.71E-21\n",
    "n2 = 4.15E-17\n",
    "n4 = 1.15E-13\n",
    "n6 = 6.88E-10\n",
    "\n",
    "print(\"\\nNumerical:\")\n",
    "print(\"TeV n2/n0 =\", n2/n0)\n",
    "print(\"TeV n4/n2 =\", n4/n2)\n",
    "print(\"TeV n6/n4 =\", n6/n4)\n",
    "\n",
    "print(\"\\nPercent errors:\")\n",
    "print(format(abs(((n2/n0)-(sigmaout[1][3]/sigmaout[0][3]))/(sigmaout[1][3]/sigmaout[0][3])), \".1%\"))\n",
    "print(format(abs(((n4/n2)-(sigmaout[2][3]/sigmaout[1][3]))/(sigmaout[2][3]/sigmaout[1][3])), \".1%\"))\n",
    "print(format(abs(((n6/n4)-(sigmaout[3][3]/sigmaout[2][3]))/(sigmaout[3][3]/sigmaout[2][3])), \".1%\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Analytical:\n",
      "TeV n2/n0 = 3789.4709190746717\n",
      "TeV n4/n0 = 10162650.893930014\n",
      "TeV n6/n0 = 22927693128.656673\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAEUCAYAAAAbV1CxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de3RU5fX/8feGYIIUUW6KBAQRVC7hFpGiVaqiWCkgWuXSUhBloUXF+m3Fn36XeCvIwm8VsSIWBRWxylJACrUgIFUREq0gAioCQtCK3AQqCIH9+2OSGMKEzGRmzkzC57XWrOQ885zn7JEjm+ecM882d0dERCQoVZIdgIiIHF+UeEREJFBKPCIiEiglHhERCZQSj4iIBEqJR0REAqXEIyIigVLiERGRQKUlO4B4MbPGwARgG/CZu49JckgiIhJGSsx4zOxZM9tqZqtKtHc3s0/NbJ2ZjSxjmBbA3939BqBlwoIVEZGYWCosmWNmFwF7gefdvXVBW1XgM6AbkAfkAP2AqsDoEkPcABwCZgAOvODuzwUTvYiIRCMlEg+AmTUB5hRLPD8FRrn7FQXbdwO4e8mkU7j//wDL3X2Jmc1w92tL6TcUGApQo0aNjuecc068P4qISKX2wQcfbHP3euXdP5Xv8TQENhfbzgPOP0b/fwCjzKw/sLG0Tu4+CZgEkJ2d7bm5ubFHKiJyHDGzL2PZP5UTj4VpK3V65u6rgLCzHBERSR0p8XBBKfKARsW2M4GvkhSLiIjESSonnhyguZk1NbMTgL7A7CTHJCIiMUqJxGNm04GlwNlmlmdmQ9w9HxgOvAmsAV5x90+SGaeIiMQuJe7xuHu/UtrnAnMDDkdERBIoJWY8IiJy/FDiERGRQCnxiIhIoJR4REQkUEo8IiISKCUeEREJlBKPiIgESolHREQCpcQjkmr2/AeeuxL2fJPsSEQSQolHJNW8PRY2vQ9vP5LsSEQSIiWWzBER4KH6kP/Dj9u5k0OvtHS4d2vy4hKJM814RFLF7Suh9a8grXpoO606tPkV3P5xcuMSiTMlHpFUUfM0SK8Jh36AtIzQz/SToOapyY5MJK50qU0klfx3K3QcDNmDIfc52KsHDKTyUeIRSSV9p/34e4//S14cIgmkS20iIhIoJR4REQmUEo9Iitm6ez/XPb2UrXv2JzsUkYRQ4hFJMePf+pycjTsYv+DzZIcikhB6uEAkRZx97zx+yD9ctP3isk28uGwT6WlV+PShK5MYmUh8acYjkiL+9cef07Pd6WRUC/1vmVGtCr3anc6/7vp5kiMTia8Km3jM7Ewzm2xmM4q11TCzqWb2jJkNSGZ8ItGqf1IGNdPT+CH/MOlpVfgh/zA109OoXzMj2aGJxFVSEo+ZPWtmW81sVYn27mb2qZmtM7ORxxrD3de7+5ASzX2AGe5+E9AzzmGLJNy2vT8w4PwzeP2WCxhw/hl8u/eHsncSqWCSdY9nCjABeL6wwcyqAk8C3YA8IMfMZgNVgdEl9r/B3cOtmpgJFC5sdSjOMYsk3NO/yS76/aHerZMYiUjiJCXxuPsSM2tSorkTsM7d1wOY2ctAL3cfDfSIcOg8QsnnI0qZzZnZUGAoQOPGjaOOXUREYpNK93gaApuLbecVtIVlZnXMbCLQ3szuLmh+DbjGzJ4C3gi3n7tPcvdsd8+uV69enEIXEZFIpdLj1BamzUvr7O7bgWEl2v4LDI5zXCIiEkepNOPJAxoV284EvkpSLCIikiCplHhygOZm1tTMTgD6ArOTHJOIiMRZsh6nng4sBc42szwzG+Lu+cBw4E1gDfCKu3+SjPhERCRxkvVUW79S2ucCcwMOR0REApRKl9pEROQ4oMQjIiKBUuIREZFAKfGIiEiglHhERCRQSjwiIhIoJR4REQmUEo+IiARKiUdERAKlxCMiIoFS4hERkUAp8YiISKCUeEREJFBKPCIiEiglHhERCZQSj4iIBEqJR0REAqXEIyIigVLiERGRQCnxiIhIoJR4REQkUBU28ZjZmWY22cxmFGvrbWbPmNksM7s8mfGJiEh4SUk8ZvasmW01s1Ul2rub2admts7MRh5rDHdf7+5DSrTNdPebgEHA9XEPXEREYpaWpONOASYAzxc2mFlV4EmgG5AH5JjZbKAqMLrE/je4+9ZjjH9vwVgiIpJikpJ43H2JmTUp0dwJWOfu6wHM7GWgl7uPBnpEMq6ZGTAGmOfuH5bSZygwFKBx48blil9ERMovle7xNAQ2F9vOK2gLy8zqmNlEoL2Z3V3QfCtwGXCtmQ0Lt5+7T3L3bHfPrlevXpxCFxGRSCXrUls4FqbNS+vs7tuBYSXaxgPj4xyXiIjEUSrNePKARsW2M4GvkhSLiIgkSColnhyguZk1NbMTgL7A7CTHJCIicZasx6mnA0uBs80sz8yGuHs+MBx4E1gDvOLunyQjPhERSZxkPdXWr5T2ucDcgMMREZEApdKlNhEROQ4o8YiISKCUeEREJFBKPCIiEiglHhERCZQSj0iqeOcx2LDkyLYNS0LtIpWIEo9IqmjYAV4d9GPy2bAktN2wQzKjEok7JR6RFDFxU0NWXfB4KNksfBheHcSqCx5n4qZS18oVqZCUeERSRFZmLQYuzGBzs36wZCybm/Vj4MIMsjJrJTs0kbhS4hFJEV2a1eX5S/bzk4+f5/3MIfzk4+d5/pL9dGlWN9mhicSVEo9IqtiwhNbv3s6CVmPou+5SFrQaQ+t3bz/6gQORCk6JRyRVbPmQVRc8zui19bntkrMYvbZ+6J7PlrDFdEUqLCUekRTxXoNfM3BhBhP6t+f3l5/NhP7tGbgwg/ca/DrZoYnElRKPSIpYmfcdE/q3L7qn06VZXSb0b8/KvO+SHJlIfJl7qdWlK73s7GzPzc1NdhgiIhWKmX3g7tnl3V8zHhERCZQSj4iIBEqJR0REAqXEIyIigVLiERGRQKVF0snMakfQ7bC774oxHhERqeQiSjzAVwUvO0afqkDjmCOKkJmdCdwD1HL3a4u11wCWAPe5+5yg4hE5Hh08eJC8vDz279+f7FAkATIyMsjMzKRatWpxHTfSxLPG3dsfq4OZ/TvSg5rZs0APYKu7ty7W3h14nFAS+6u7jyltDHdfDwwxsxkl3roLeCXSWESk/PLy8qhZsyZNmjTB7Fj/LpWKxt3Zvn07eXl5NG3aNK5jR3qP56dx6lNoCtC9eIOZVQWeBK4EWgL9zKylmbUxszklXvXDDWpmlwGrgW+iiEVEymn//v3UqVNHSacSMjPq1KmTkNlsRDMedy/zyJH0KdZ3iZk1KdHcCVhXMJPBzF4Gern7aEKzo0j8HKhBKHHtM7O57n64eAczGwoMBWjcOLArgyKVlpJO5ZWoP9syZzxm1s3MnjGzdgXbQxMSCTQENhfbzitoKy2uOmY2EWhvZncDuPs97j4CeAl4pmTSKegzyd2z3T27Xr168f0EIiJSpkgutd0C/AH4tZldArRLUCzhUmupC8m5+3Z3H+buzQpmRcXfm6IHC0SOH6+//jpmxtq1a2MaZ9CgQcyYUfK28ZH+9Kc/HbHdpUuXch1r1KhRjBs3rlz7xsLdue222zjrrLPIysriww+DL7sRSeL51t13ufv/AJcD5yUoljygUbHtTEJP0olIJfLBlzt5ctE6PvhyZ9zGnD59OhdeeCEvv/xy3MYsTcnE89577yX8mPE0b948Pv/8cz7//HMmTZrEzTffHHgMkdzj+XvhL+4+0sxuTVAsOUBzM2sKbAH6Av0TdCwRibP73/iE1V/tPmafPfsPsvY/ezjsUMXgnNNqUjOj9Ed1W55+Evf9stUxx9y7dy/vvvsuixYtomfPnowaNQqAxYsXM2rUKOrWrcuqVavo2LEjL774ImbGAw88wBtvvMG+ffvo0qULTz/99BH3M9566y0mTJjA66+/DsD8+fN56qmnaNGiBfv27aNdu3a0atWKadOm8ZOf/IS9e/cCMHbsWF544QWqVKnClVdeyZgxY3jmmWeYNGkSBw4c4KyzzuKFF17gxBNPjOQ/KaNGjWLTpk2sX7+eTZs2MWLECG677baI9i3NrFmzGDhwIGZG586d2bVrF19//TUNGjSIadxolDnjcfdZJbafiPWgZjYdWAqcbWZ5ZjbE3fOB4cCbwBrgFXf/JNZjiUjq2L0/n8MFF9APe2g7VjNnzqR79+60aNGC2rVrH3Hp6N///jePPfYYq1evZv369bz77rsADB8+nJycHFatWsW+ffuYM+fIK/OXXHIJa9as4dtvvwXgueeeY/DgwYwZM4bq1avz0UcfMW3atCP2mTdvHjNnzmTZsmWsWLGCP/7xjwD06dOHnJwcVqxYwbnnnsvkyZOj+nxr167lzTffZPny5dx///0cPHjwqD7XX3897dq1O+r1/PPPH9V3y5YtNGr048WlzMxMtmzZElVMsYr0ezwAmFk2oS9tnlGwrwHu7lnRjOPu/UppnwvMjWYsEUkNZc1MIHSZbcBf3+dg/mGqpVXh8b7t6XjGKTEdd/r06YwYMQKAvn37Mn36dDp06ABAp06dyMzMBKBdu3Zs3LiRCy+8kEWLFjF27Fi+//57duzYQatWrfjlL39ZNKaZ8Zvf/IYXX3yRwYMHs3Tp0rB/iRe3YMECBg8eXDSbqV07tODLqlWruPfee9m1axd79+7liiuuiOrzXXXVVaSnp5Oenk79+vX55ptvij5Tob/97W8RjxeuBlvQTyZGlXiAaYQeNPgYOOqJMRGRY+l4xilMu7Ez76/fTucz68ScdLZv387ChQtZtWoVZsahQ4cwM8aOHQtAenp6Ud+qVauSn5/P/v37ueWWW8jNzaVRo0aMGjUq7HdVBg8ezC9/+UsyMjL41a9+RVrasf+6dPewf4EPGjSImTNn0rZtW6ZMmcLixYuj+ozhPkNJ119/PZ9++ulR7b///e8ZOHDgEW2ZmZls3vzjA8R5eXmcfvrpUcUUq2gTz7fuPjshkYjIcaHjGafEnHAKzZgxg4EDB/L0008XtV188cW88847pe5TmGTq1q3L3r17mTFjBtdee+1R/U4//XROP/10HnroIebPn1/UXq1aNQ4ePHjUMjKXX345DzzwAP379+fEE09kx44d1K5dmz179tCgQQMOHjzItGnTaNjw6G+JTJgwAQhdAiyPaGY8PXv2ZMKECfTt25dly5ZRq1atQO/vQPSrU99nZn81s35m1qfwlZDIRETKMH36dK6++uoj2q655hpeeumlUvc5+eSTuemmm2jTpg29e/fmvPNKf1B3wIABNGrUiJYtWxa1DR06lKysLAYMGHBE3+7du9OzZ0+ys7Np165d0aPSDz74IOeffz7dunXjnHPOCXuctWvXUqdOnTI/bzz84he/4Mwzz+Sss87ipptu4i9/+Usgxy3Owl3vK7Wz2YvAOcAn/Hipzd39hgTElnDZ2dmem5ub7DBEKqw1a9Zw7rnnJjuMhBk+fDjt27dnyJAhCT1Ojx49eO211zjhhBMSepzyCPdnbGYfuHt2eceM9lJbW3dvU96DiYhUFB07dqRGjRo8+uijCT9WyafqKrtoE8/7ZtbS3VcnJBoRkRTxwQcfJDuESivaxHMh8Fsz2wD8QDkfpxYRkeNXtImne9ldREREShftU20PAN+5+5fu/iWwG7gv/mGJiEhlFW3iyXL3XYUb7r4TOGZlUhERkeKiTTxVzKzom19mVpvoL9eJyHHuz/M/i+t4KosQuWnTppGVlUVWVhZdunRhxYoVgccQbeJ5FHjPzB40sweA94Cx8Q9LRCqzx9/6PK7jqSxC5Jo2bcrbb7/NypUr+d///V+GDk1Ubc/SRTVbcffnzSwXuITQE2199Gi1iEBkZRGKu/7ppWX2UVmE+JdFKD5D69y5M3l5eTGNVx7Rznhw99XuPsHdn1DSEZFI5e38nmUbdrBsww6Aot/zdn4f07gqixBdWYTiJk+ezJVXXhlVPPEQ0YzHzD509w6x9hGRyiuSsgiFmoz8OxvHXBWX46osQnSLhBZatGgRkydPPuaCqokS6aW2c81s5THeN6BWHOIREYmYyiKERFMWAWDlypXceOONzJs3L7DFSYuLNPGEX1L1SIdiCUREjh+3X9o8LuOoLEJINDOeTZs20adPH1544QVatGhRruPFKqJ7PIVfGC3jFfwdKhGpkO7oFp+/8FQWIXoPPPAA27dv55ZbbqFdu3ZkZ5d7kelyi6osQmWjsggisVFZhPhQWQQREVFZhARS4hERCUNlERKn3InHzPoDPQk9VGDAG+4+PV6BRXD8M4F7gFrufm1BWxXgQeAkINfdpwYVj4iIRCbqL5AWc7G793X3Ae7en1CtnoiY2bNmttXMVpVo725mn5rZOjMbeawx3H29u5e88NoLaAgcBPSwg4hICorlUlu6mV0FbAYygepR7DsFmAAUfSPLzKoCTwLdCCWNHDObDVQFRpfY/wZ33xpm3LOBpe7+tJnNAN6KIiYREQlAVDMeMytee+cW4BTgFwU/I34A3d2XADtKNHcC1hXMZA4ALwO93P1jd+9R4hUu6UAoYe0s+D3s94rMbKiZ5ZpZbuFyGCIiEpxoL7XdZ2aPmNkzwG+Bv7v7GHef5u6xLbgUukS2udh2XkFbWGZWx8wmAu3N7O6C5teAK8zsCWBJuP3cfZK7Z7t7dr169WIMWUQiNfHtL3jvi21HtL33xTYmvv1FzGOrLEL0cnJyqFq1apmfNxGiTTwO7AfeBBoRKpHQNk6xHL3WROh44QNx3+7uw9y9mbuPLmj73t2HuPut7v5knOISkTjIyqzF8Jf+XZR83vtiG8Nf+jdZmbGvtqWyCNE5dOgQd911V9TrxsVLtPd41rp74eW2GWY2BZhIqExCrPIIJbNCmcBXcRhXpOJZNBp+fnfZ/VJIJGUR6tdMZ+Dk5Zx6Ujrf7P6Bs+r/hMcXfM7jC8LX51FZhPiXRQB44oknuOaaa8jJyYl5rPKIdsazzcw6Fm64+2dAvK5X5QDNzaypmZ0A9AVmx2lskYrl7THJjiAhalWvxqknpbNl135OPSmdWtWrlb1TGVQWIbqyCFu2bOH1119n2LBhUcURT9HOeG4DXjazD4CPgSxgQ7QHNbPpQFegrpnlAfe5+2QzG07oMl5V4Fl3/yTasUUkOSIpi1B4ee22S87ixWWbuP2y5nRpVjem46osQnSLhI4YMYJHHnmEqlWrRhVHPEVbgXSFmbUDLgNaA4uAqL806u79SmmfC8yNdjyRSmHR6CNnOqMK7n1cPLLCXXYLpzDpTOjfni7N6tK5WZ0jtstDZRFCoimLkJubS9++fQHYtm0bc+fOJS0tjd69e0cVVyyi/h6Pu/8A/L3gJSLx8vO7f0wwo2rBqO+SG0+crcz77ogk06VZXSb0b8/KvO/KnXhUFiEkmhnPhg0/XqQaNGgQPXr0CDTpQGwrF4hIHP15/mc0Gfl3mowM/Zuu8Pc/z/8syZHFx7CLmx2VYLo0q8uwi5uVe0yVRaiYVBZBZREkBT1272BGPPRcssMok8oixMfxVhZBMx6RFPRY/tGXfiRYHTt2ZOXKlfz6179O+LHmzJmTkkknUVQWQSQFxas0tJSfyiIkjmY8IikoXqWhRVKREo+IiARKiUdERAKlxCMiIoFS4hGRQG3dvZ/rnl7K1j1HrxZQXiqLEJ3FixcXLXR68cUXB358JR4RCdT4tz4nZ+MOxpeyInV5qCxC5Hbt2sUtt9zC7Nmz+eSTT3j11VcDj0GPU4tIXJRVFmH5xh0U/776i8s28eKyTZhBpya1w+6jsgjxL4vw0ksv0adPHxo3bgxA/fr1YxqvPDTjEZFAtMs8mTo1TqBKwd/vVQzq1DiBdpknxzSuyiJEVxbhs88+Y+fOnXTt2pWOHTuWuep2ImjGIyJxEUlZhHte/5iXlm8iPa0KBw4d5srWp/HQ1W1iOq7KIkS3SGh+fj4ffPABb731Fvv27eOnP/0pnTt3pkWL4L47psQjIoHZtvcHBpx/Bv07Neal5Zv4NsYHDFQWISSasgiZmZnUrVuXGjVqUKNGDS666CJWrFihxCMildPTv/lxXcmHereOeTyVRQiJZsbTq1cvhg8fTn5+PgcOHGDZsmXccccd5Tpueekej4hUWCqLEL1zzz2X7t27k5WVRadOnbjxxhtp3Tr2fwREQ2URVBZBpNxUFiE+jreyCLrUJiISRseOHalRowaPPvpowo9V8qm6yk6JR0QkDJVFSBzd4xERkUBV2MRjZmea2WQzm1GsrbGZzTazZ81sZDLjExGR8JKSeAoSw1YzW1WivbuZfWpm68pKHO6+3t1L3vFrAfzd3W8AWobZTUREkixZM54pQPfiDWZWFXgSuJJQ0uhnZi3NrI2ZzSnxKm1xoX8Dfc1sIbAogfGLiEg5JSXxuPsSYEeJ5k7AuoKZzAHgZaCXu3/s7j1KvLaWMvRg4D53vwS4KlwHMxtqZrlmllu4DpOIVFxmxp133lm0PW7cuKKFQoOSm5tb7sU7u3btSjK+1rFjxw66detG8+bN6datGzt37gzs2Kl0j6chsLnYdl5BW1hmVsfMJgLtzezuguZ/ALcVtG8Mt5+7T3L3bHfPrlevXnwiF5HIbV4O/3o09DMO0tPTee2119i2bVtcxotWfn4+2dnZjB8/PinHL68xY8Zw6aWX8vnnn3PppZcyZsyYwI6dSo9TH73IEZT67VZ33w4MK9G2Cjh67QsRSbx5I+E/Hx+7zw+74ZtV4IfBqsCprSH9pNL7n9YGrjz2X4hpaWkMHTqUP//5zzz88MNHvDdo0CB69OhRtCROYQmDxYsXc99993Hqqafy0Ucf0adPH9q0acPjjz/Ovn37mDlzJs2aNePbb79l2LBhbNq0CYDHHnuMCy64gFGjRvHVV1+xceNG6taty9ChQxk3bhxz5sxh79693HrrreTm5mJm3HfffVxzzTXcfPPN5OTksG/fPq699lruv//+sv+bFujatSvnn38+ixYtYteuXUyePJmf/exnEe8fzqxZs4rWjfvtb39L165deeSRR2IaM1KplHjygEbFtjOBr5IUi4gkwv7vQkkHQj/3f3fsxBOh3/3ud2RlZRWVIojEihUrWLNmDbVr1+bMM8/kxhtvZPny5Tz++OM88cQTPPbYY9x+++3ccccdXHjhhWzatIkrrriCNWvWAKHv+bzzzjtUr179iIU/H3zwQWrVqsXHH4eScOElrIcffpjatWtz6NAhLr30UlauXElWVlbE8ebn57N8+XLmzp3L/fffz4IFC454f8+ePaUmo5deeumIZX8AvvnmGxo0aABAgwYN2Lq1tDsY8ZdKiScHaG5mTYEtQF+gf3JDEpGIlTEzAUKX16b2hEMHoOoJcM1foVGnmA990kknMXDgQMaPH0/16tUj2ue8884r+ou3WbNmXH755QC0adOGRYtCzyYtWLCA1atXF+2ze/du9uzZA0DPnj3DHmvBggVHVEI95ZRTAHjllVeYNGkS+fn5fP3116xevTqqxNOnTx8gtKLCxo0bj3q/Zs2afPTRRxGPl0xJSTxmNh3oCtQ1szxCDwRMNrPhwJtAVeBZd/8kGfGJSII06gS/nQ0b/wVNfhaXpFNoxIgRdOjQgcGDBxe1paWlcfhwaIbl7hw4cKDoveLlBqpUqVK0XaVKlaLSA4cPH2bp0qVhE0yNGjXCxhGuPMKGDRsYN24cOTk5nHLKKQwaNChsKYZjKYyvtNII0c54Tj31VL7++msaNGjA119/HWgl0mQ91dbP3Ru4ezV3z3T3yQXtc929hbs3c/eHyxpHRCqgRp3gZ3fGNelAqPDaddddd0SFzyZNmhQtfTNr1qyw1TuP5fLLLy8qWQBENKMouc/OnTvZvXs3NWrUoFatWnzzzTfMmzcv7L4DBw5k+fLyPXRROOMJ9yqZdCA0Y5s6dSoAU6dOpVevXuU6bnmk0lNtIiIxufPOO494uu2mm27i7bffplOnTixbtqzUWUppxo8fT25uLllZWbRs2ZKJEyeWuc+9997Lzp07ad26NW3btmXRokW0bduW9u3b06pVK2644QYuuOCCsPuuXLmy6PJfoo0cOZL58+fTvHlz5s+fz8iRwS32orIIKosgUm6VvSxCkHbv3s2QIUN49dVXkx3KERJRFkEzHhGRFHDSSSelXNJJFCUeEREJlBKPiIgESolHREQCpcQjIiKBUuIREZFAKfGISPAWjY7bUCqLUD5/+MMfOOecc8jKyuLqq69m165dgR1biUdEgvd2/JbgV1mE8unWrRurVq1i5cqVtGjRgtGj4/ePgbKk0iKhIlKRRVIWobjnwtZqPJLKIgCJKYtQuCgqQOfOnZkxY0ZM40VDiUdEgrHrS/iuWK3HL98J/azVCE4+I6ahVRYh+kVCi3v22We5/vrrI44lVko8IhIfkZRFKDSqFoz6Lm6HVlmE8pdFePjhh0lLS2PAgAFR71teSjwiUimoLEL0M56pU6cyZ84c3nrrraNiTiQ9XCAiwbs4/ishqyxCdGUR/vGPf/DII48we/ZsTjzxxHIds7yUeEQkeD+/OyHDqixC5IYPH86ePXvo1q0b7dq1Y9iwYYEcF1QWQWURRGKgsgjxo7IIIiISKJVFEBERSRAlHhERCZQSj4iIBKrCfo/HzHoDVwH1gSfd/Z9mVgP4C3AAWOzu05IZo4iIHC0pMx4ze9bMtprZqhLt3c3sUzNbZ2bHfNDf3We6+03AIKBwrYc+wIyC9p6JiF1ERGKTrEttU4DuxRvMrCrwJHAl0BLoZ2YtzayNmc0p8apfbNd7C/YDyAQKF4M6lNBPIJIga3MWsHTq/2NtzoKyO1ck7zwGG5Yc2bZhSag9BiqLEJtx48ZhZoGu7p2UxOPuS4AdJZo7Aevcfb27HwBeBnq5+8fu3qPEa6uFPALMc/cPC8bII5R8QPevpAJam7OAM+b047z1T3HGnH6VK/k07ACvDvox+WxYEtpu2CGmYVUWofw2b97M/Pnzady4caDHTaV7PA35cbYCoSRy/jH63wpcBtQys7PcfSLwGjDBzK4C3gi3k5kNBYYCgf/HFinLztULOYt80uwweD47Vy+E8y5LdliRiaQsQs0G8ALJVHUAAAmoSURBVMLVoZ97voZ658DiR0KvcFQWAUhMWQSAO+64g7Fjx9KrV6+Yx4pGKiWecCvUlbqsgruPB8aXaPsvMDj8HkV9JgGTILRyQfRhiiTOKS0v4eD6Z8DzOUgap7S8JNkhxVfGyaGk893mUDmEjJPjMqzKIkS/SOjs2bNp2LAhbdu2jTiGeEmlxJMHNCq2nQl8laRYRJLinPMuYy3T2bl6Iae0vIRzKspsByIri1B4ee2iP0LuZOh6FzS9KOZDqyxCdGURvv/+ex5++GH++c9/Rnz8eEqlxJMDNDezpsAWoC/QP7khiQTvnPMuqziX16JRmHR+NSWUbJr+7MjtGKksQuQzni+++IINGzYUzXby8vLo0KEDy5cv57TTTosqrvJI1uPU04GlwNlmlmdmQ9w9HxgOvAmsAV5x90+SEZ+IJMCWD49MMk0vCm1v+fBYe0VMZREiL4vQpk0btm7dysaNG9m4cSOZmZl8+OGHgSQdSN5Tbf3cvYG7V3P3THefXNA+191buHszd3+4rHFEpAK5cMTRM5umF4Xa40RlESoGlUVQWQSRclNZhPhRWQQREQmUyiKIiIgkiBKPiIgESolHREQCpcQjIiKBUuIREZFAKfGISLD2/AeeuxL2fBOX4VQWofyeeOIJzj77bFq1ahXVOnexSqUlc0TkePD2WNj0Prz9CPT4v5iHKyyLcPfdd1O3bt04BBidwrII2dnl/lpLUixatIhZs2axcuVK0tPT2bp1a2DHVuIRkfgoqyzCpneh+BfWcyeHXmbQOPw3+VUWISQRZRGeeuopRo4cWbQGXP369cvYI36UeEQkGKefBzs3wL7t4IfBqsCJdeCUpjEPrbII0ZdF+Oyzz/jXv/7FPffcQ0ZGBuPGjeO8886LOJ5YKPGISHxEUhbhjTvgwymQlgGHDsC5PeNyuU1lEaIriwChRLZz507ef/99cnJyuO6661i/fv1RK2snghKPiATnv1uh42DIHgy5z8He+DxgACqLEO2MJzMzkz59+mBmdOrUiSpVqrBt2zbq1asXVVzloafaRCQ4faeFZjintQn97DstbkOrLELkZREAevfuzcKFC4HQZbcDBw4E9nCGEo+IVBoqixC5G264gfXr19O6dWv69u3L1KlTA7nMBiqLoLIIIjFQWYT4UVkEEREJlMoiiIiIJIgSj4jE5Hi+XF/ZJerPVolHRMotIyOD7du3K/lUQu7O9u3bycjIiPvY+h6PiJRbZmYmeXl5fPvtt8kORRIgIyODzMzMuI9bYROPmfUGrgLqA0+6+z/DtSUzRpHKrlq1ajRtGvuSN3J8ScqlNjN71sy2mtmqEu3dzexTM1tnZiOPNYa7z3T3m4BBwPWltYmISGpJ1oxnCjABeL6wwcyqAk8C3YA8IMfMZgNVgdEl9r/B3QvX8L63YL/iwrWJiEgKSEricfclZtakRHMnYJ27rwcws5eBXu4+GuhRcgwLfcV2DDDP3T8srU1ERFJLKt3jaQhsLradB5x/jP63ApcBtczsLHefWErbEcxsKDC0YPOHkpf7KqhawHeV5Lixjlme/aPZJ9K+kfQ7Vp+6wLZS3qtoknF+VpZzM9r94nV+lvX+2RHGE567J+UFNAFWFdv+FfDXYtu/AZ5IcAy5yfr8cf4ckyrLcWMdszz7R7NPpH0j6XesPpXl3IzHn2mqHDMZ52a0+8Xr/Izg/ZjOz1T6Hk8e0KjYdibwVZJiqWjeqETHjXXM8uwfzT6R9o2kX7L+3IKWjM9ZWc7NaPeL1/mZ0D+zpC0SWnCPZ467ty7YTgM+Ay4FtgA5QH93/ySBMeR6DAvdiSSKzk1JZbGen8l6nHo6sBQ428zyzGyIu+cDw4E3gTXAK4lMOgUmJXh8kfLSuSmpLKbz87guiyAiIsFLpXs8IiJyHFDiERGRQCnxiIhIoJR4ijGz3mb2jJnNMrPLkx2PSCEzO9fMJprZDDO7OdnxiBRnZjXM7AMzO2qVmXAqfeKJZkFS1yKjEqAoz8017j4MuA7QY9aSUOVYyPku4JVIx6/0iYfQgqTdizcUW5D0SqAl0M/MWhbrokVGJQhTiOLcNLOewDvAW8GGKcehKUR4bprZZcBq4JtIB6/0icfdlwA7SjQXLUjq7geAl4FeFvIIWmRUAhDNuVnQf7a7dwEGBBupHG+iPDd/DnQG+gM3mVmZeSWVFgkNUmkLkpa5yKhIgoU9N82sK9AHSAfmJiEukbDnprsPBzCzQcA2dz9c1kDHa+KxMG3u7uOB8UEHI1JMaefmYmBxsKGIHCHsuVn0i/uUSAeq9JfaSqEFSSVV6dyUVBW3c/N4TTw5QHMza2pmJwB9gdlJjkkEdG5K6orbuVnpE08KLUgqcgSdm5KqEn1uapFQEREJVKWf8YiISGpR4hERkUAp8YiISKCUeEREJFBKPCIiEiglHhERCZQSj4iIBEqJR0REAqXEI5JkZtbEzPaZ2UfF2k41s5fMbH1BZcelZnZ1GeMsNrMrSrSNMLO/mFl1M/vIzA6YWd1EfRaRSCjxiKSGL9y9HYCZGTATWOLuZ7p7R0LrYmWWMcb0gn7F9QWmu/u+gvG14KgknRKPSIzM7FUzm2Bm75jZl2Z2oZk9b2afmdnkcgx5CXCgeD0od//S3Z8odsxfm9nyglnM0wXVIWcAPcwsvaBPE+B0QlVLRVKGEo9I7NoA6939QmAqMJlQDfrWQJ/CRBCFVkCpFXDN7FzgeuCCglnMIWCAu28HlvNjyeK+wN9cCzJKijleC8GJxIWZZQAnA48VNO0DJrv71wXvfw8ciPEYTwIXEpoFnQdcCnQEckJX5agObC3oXni5bVbBzxtiObZIImjGIxKbVsCHxcr9tgWWAZhZYaGsLDN73cyqmdmLEdzc/wToULjh7r8jlGzqFTQZMNXd2xW8znb3UQXvzQQuNbMOQHV3L3XmJJIsSjwisWkDrCi2nQWsLPi9LbDS3VcAHwHPAY+5+7YyxlwIZJjZzcXaTiz2+1vAtWZWH8DMapvZGQDuvpdQiexnCc1+RFKOEo9IbNoQSiqFl92qu/vOgveKJ6GdwBp3zy1rwIJ7Mr2Bi81sg5ktJ3Tv6K6C91cD9wL/NLOVwHygQbEhphNKei/H+NlEEkKF4EQSzMxaAncC2939j2HebwLMcffWAcSyEciOYNYlkjCa8YgkkJnVAu4BbgGamtn5YbodAmoV/wJpAuKoXjB+NeBwWf1FEkkzHhERCZRmPCIiEiglHhERCZQSj4iIBEqJR0REAqXEIyIigVLiERGRQCnxiIhIoJR4REQkUP8fELOH6x8cqZoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "an2n0 = sigmaout[1][3]/sigmaout[0][3]\n",
    "an4n0 = sigmaout[2][3]/sigmaout[0][3]\n",
    "an6n0 = sigmaout[3][3]/sigmaout[0][3]\n",
    "\n",
    "print(\"Analytical:\")\n",
    "print(\"TeV n2/n0 =\", an2n0)\n",
    "print(\"TeV n4/n0 =\", an4n0)\n",
    "print(\"TeV n6/n0 =\", an6n0)\n",
    "\n",
    "n0 = 8.71E-21\n",
    "n2 = 4.15E-17\n",
    "n4 = 1.15E-13\n",
    "n6 = 6.88E-10\n",
    "\n",
    "nn2 = n0*an2n0\n",
    "nn4 = n0*an4n0\n",
    "nn6 = n0*an6n0\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "f,ax = plt.subplots()\n",
    "\n",
    "ax.plot(1.0e3, n0, marker='.', label='Analytical, n = 0', color='tab:blue')\n",
    "ax.plot(1.0e3, nn2, marker='+', label='Analytical, n = 2', color='tab:blue')\n",
    "ax.plot(1.0e3, nn4, marker='x', label='Analytical, n = 4', color='tab:blue')\n",
    "ax.plot(1.0e3, nn6, marker='*', label='Analytical, n = 6', color='tab:blue')\n",
    "\n",
    "ax.plot(1.0e3, n0, marker='.', label='Numerical, n = 0', color='tab:orange')\n",
    "ax.plot(1.0e3, n2, marker='+', label='Numerical, n = 2', color='tab:orange')\n",
    "ax.plot(1.0e3, n4, marker='x', label='Numerical, n = 4', color='tab:orange')\n",
    "ax.plot(1.0e3, n6, marker='*', label='Numerical, n = 6', color='tab:orange')\n",
    "\n",
    "#ax.set_title('DM-electron scattering (n=0, Vrel=0)')\n",
    "ax.set_xlabel('$m_\\chi$ $[\\mathrm{GeV}]$')\n",
    "ax.set_ylabel('$\\sigma_0$ $[\\mathrm{cm}^2]$')\n",
    "ax.set_xscale('log')\n",
    "ax.set_yscale('log')\n",
    "ax.set_xlim([1.0e2, 1.0e4])\n",
    "ax.set_ylim([1.0e-22, 1.0e-8])\n",
    "plt.legend(loc='lower right')\n",
    "plt.savefig('dme_analytical_TeV_comparison.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}